使用拍照手機(jī)的數(shù)字圖像相關(guān)測量系統(tǒng)在實(shí)驗(yàn)力學(xué)教學(xué)中的應(yīng)用

1、基于Matlab求解時(shí)諧復(fù)系數(shù)彈性波方程的超弱變分方法摘要 眾所周知，彈性波方程被用來描述工程應(yīng)用中的彈性波傳輸問題。高效數(shù)值求解彈性波方程是一個(gè) 重要的研究課題。提出了基于平面波離散的間斷有限元超弱變分法。首先，選取滿足剖分單元上齊次偏微分方 程和其對(duì)偶方程的解析解分別作為檢驗(yàn)空間和測試空間，得到等價(jià)的間斷有限元變分形式；然后分別定義二維 空間下的矢量平面波檢驗(yàn)離散空間和測試離散空間，進(jìn)而得到離散變分形式?；贛atlab軟件的數(shù)值實(shí)驗(yàn)驗(yàn) 證了新算法具有較高的精度。關(guān)鍵詞 彈性波方程，超弱變分方法，平面波AN ULTRA WEAK VARIATIONAL FORMULATION FOR SO

2、LVINGTIME-HARMONIC ELASTIC WAVE EQUATIONS WITH COMPLEXWAVE COEFFICIENTS BASED ON MATLAB1）Abstract As is known, some engineering problems, such as the elastic wave transmission, can be described by the elastic wave equation. The efficient numerical algorithms for solving elastic wave equations are im

3、portant. An ultra weak variational formulation based on the plane wave discretizations is proposed. First, by choosing the analytic solutions satisfying the homogeneous elastic wave equations and its adjoint equations to form the trial space and the test space, respectively, we obtain the equivalent

4、 Trefftz variational formulation; then we define the two-dimensional vector plane wave discretized trial space and test space to discretize the proposed continuous variational formulation. Numerical results based on Matlab have validated the proposed method.Key words elastic wave equations, ultra we

5、ak variational formulation, plane wave“偏微分方程數(shù)值解”課程是數(shù)學(xué)與應(yīng)用數(shù)學(xué)、 信息與計(jì)算科學(xué)等專業(yè)的核心課程。該課程針對(duì)刻 畫物理模型的數(shù)學(xué)方程，利用有限元、有限差分法、 有限體積法和邊界元法等方法數(shù)值求解方程，并理 論分析相應(yīng)數(shù)值解的誤差收斂階。該課程涉及知識(shí) 面廣泛，不但要講授變分方法、索伯列夫空間、離散 空間（如多項(xiàng)式空間）、線性方程組的迭代求解及數(shù) 值逼近，還要培養(yǎng)學(xué)生的編程驗(yàn)證能力【IT。因此,學(xué)生系統(tǒng)地掌握該課程的知識(shí)結(jié)構(gòu)是至關(guān)重要的。本文基于間斷有限元方法數(shù)值求解描述一大類 物理問題模型(比如聲波散射、彈性波散射問題)的 彈性波方程。在數(shù)

6、值求解彈性波方程的發(fā)展過程中, 有限元I4-5、邊界元法6-7和有限差分方法I起到 了很重要的作用。彈性波方程解析解的一個(gè)重要參 數(shù)是波數(shù)。它描述了解析解的振蕩行為。波數(shù)越大, 振蕩越頻繁。這個(gè)特點(diǎn)在有限元分析中可以通過經(jīng) 驗(yàn)法則來刻畫一一單個(gè)波長內(nèi)的網(wǎng)格點(diǎn)數(shù)。因而由 此經(jīng)驗(yàn)法則很難獲得大波數(shù)情形下令人滿意的數(shù)值 結(jié)果。為了克服此缺點(diǎn)，很多改進(jìn)方法10-12被陸續(xù) 提出。最近用來求解中高頻彈性波方程的平面波方 法13-15在工業(yè)界和數(shù)學(xué)界比較流行。平面波方法的 核心是選取滿足無約束方程的精確解作為逼近基函 數(shù)，如平面波函數(shù)、貝塞爾函數(shù)、格林函數(shù)等，并且 未知量是定義在剖分單元上的。大量的數(shù)值實(shí)

7、驗(yàn)表 明，相比其他離散化方法，在滿足相同精度的條件 下，平面波方法使用了相對(duì)較少的自由度；即，在選 取相同自由度的條件下，平面波離散化方法具有較 高的精度。為了讓學(xué)生直觀感受平面波方法數(shù)值求解偏微 分方程的過程，本文從方程組本身出發(fā)，分別構(gòu)造滿 足彈性波方程和其對(duì)偶方程的測試函數(shù)空間和檢驗(yàn) 函數(shù)空間，并基于間斷有限元思想，結(jié)合分部積分公 式，推導(dǎo)出復(fù)波數(shù)情形的超弱變分形式；分別構(gòu)造二 維空間下的平面波檢驗(yàn)離散函數(shù)空間和測試離散函 數(shù)空間，進(jìn)一步得到等價(jià)的變分形式。數(shù)值實(shí)驗(yàn)驗(yàn)證 了新方法的數(shù)值解具有較高的精度。1控制方程與數(shù)值方法1.1彈性波方程考慮邊值問題16 TOC o 1-5 h z IA

8、u + (A + ) V (V u) + pu = 0(1)及最低階吸收邊界條件Tn (u) iau = g(2)其中，拉梅常數(shù)入和產(chǎn)借助于泊松比和楊氏模量 E表達(dá)E 、Ev必=2(1 + v)= (1 + v)(1 2v) 材料密度p 0是常數(shù)，3 0是角頻率，n是邊界 的單位外法向，i是虛部單位，截?cái)嗨阕覶n (u)為uTn (u) = 2# + XnV-u + 訴n x 如乂 u)on定義縱波和橫波的波速分別為Cp/a + 2#Cp/a + 2#a = 3p (CP | nnT + |Cs| ssT)(3)其中s是邊界的單位切向量，|Cp|與|CS|分別表 示Cp與Cs的模，T代表的轉(zhuǎn)置

9、。假定求解區(qū) 域是有界多面體區(qū)域。Th表示區(qū)域的三角剖 分N=1，并滿足“形正則”和“擬一致”假設(shè)17o 定義rtj = n 阻(i = j)，Yk =頁n dd 為簡單 起見，我們只考慮二維區(qū)域模型，但方法可以直接推 廣到三維空間。1.2彈性波方程的變分形式對(duì)于每一個(gè)剖分單元d，令H 1(dk)為標(biāo)準(zhǔn)的 向量索伯列夫空間，并定義分片間斷Trefftz變分空 間18V(dk) = uk e H 1(dk)； uk 滿足(1) 和區(qū)域間斷Trefftz變分空間V(Th) = I!N=1 V (d) 類似地，定義分片間斷Trefftz測試空間W (風(fēng))為 滿足方程(1)的對(duì)偶方程#Aw + (A

10、+ #) V (V w) + 32pw = 0和區(qū)域間斷Trefftz測試空間NW (Th) = H W (dk)k=1對(duì)于每個(gè)單元d，定義跡算子L+Vk = Tnk (vk) iavk L-Vk = Tnk (vk) iavk M+Vk = Tnk (vk) iavk M-vk = Tnk (vk) iavk利用等距引理19-20，原問題(1)(2)等價(jià)于變分形式: 尋找u & V(Th)滿足A (u, v) = l (v), Vv e W(Th)(4)其中半雙線性型N fA (u, v) = /。 I a-1 (L+u -(昭F)ds,l=3rVv G W(Th)右端項(xiàng)為N ，- I (v

11、) = a-19 - (M-勿)ds k=J Y1.3平面波離散化本節(jié)構(gòu)造間斷Trefftz變分空間V(Th)的平面 波離散空間Vp(Th)。首先考慮沿單位方向d傳播的 時(shí)諧彈性平面波，其可分解為縱波與橫波之和v = odexp (iKpd - x) + /3e exp (iKsd - x)其中 Kp = 3/Cp, Ks = 3/Cs， Q, 是系數(shù)，d - e = 0o 縱波 Vp = odexp(i%d - x)和橫波 Vs = 3e exp (iKsd .x)。簡單計(jì)算可知，縱波和橫波均滿足 原始方程(1)，且Nx Vp = 0,Vs =0。理論分析 和實(shí)際應(yīng)用表明，借助于選取單位圓周

12、上分布均勻 的一簇互不相同的彈性波傳播方向定義解析解Uk = unk的平面波逼近Uh,k TOC o 1-5 h z ppuh,k =ak,lvk,l +Bk,lvk,li=1i=1其中 ak,i, Pk,i (k = 1, 2,N; I = 1,2,，p)為待 求解的系數(shù)。對(duì)應(yīng)于連續(xù)變分形式(4)的離散形式 為：尋找Uh G Vp(Th)滿足A(uh, vh) = l (vh), Vvh G Wp(Th)(5)我們利用直接法，即高斯消去法少,求解上述 離散線性方程組。1.4算法實(shí)現(xiàn)步驟1.初始化輸入?yún)?shù)：計(jì)算區(qū)域Q,入心,P, b；步驟2.剖分區(qū)域Q形成三角網(wǎng)格Th；步驟3.基于變分形式(5

13、)的間斷有限元離散, 形成剛度矩陣A及右端項(xiàng)b；步驟4.利用線性解法器求解線性方程組Ax =b；步驟5.利用解向量構(gòu)造有限元解Uh并計(jì)算誤 差；步驟6.分析誤差關(guān)于p及h的變化。2數(shù)值算例與結(jié)果分析cos(2%(/ _ 1)/p) sin(2%(/ 1)/p)取& = 1, & = 1，計(jì)算區(qū)域Q =0,12o取 Rayleigh波作為邊值問題和的解析解22我們得到定義在每個(gè)單元Qk上的2p個(gè)復(fù)平面波基 函數(shù)uxuyJ vp t = di exp (iKpdi x)vkj = ei exp (iKsdi x)(l = 1, 2, ,，)記Qq(q = 2p)表示單元Qk上的上述2p個(gè)復(fù)平面波基

14、函數(shù)張成的空間，則平面波離散空間Vp(Th)為Vp(Th) = v G H1 (Q) : vK G Qq, VK G Thas (e-asy 22Kr 2 e-ap八 eiKRX (KR + aSiKR(e-as y _ 2ap aseapy ) eiKx kR + a%)其中0.87+1.123 KR = cRkR Kp, askR _ KS類似地，為了滿足方程的對(duì)偶方程，我們只需要 將vpi和vkj中的Kp和Ks分別替換為Kp和瓦 即可,進(jìn)而得到單元Qk上的上述2p個(gè)復(fù)平面波基 函數(shù)張成的空間Sq(q = 2p)。平面波離散測試空間 Wp(Th)為Wp(Th) = w G H1 (Q) :

15、 wk G Sq, VK G Th 取 E = 2.0 x 1011, = 0.3, p = 7800。將解析解 u 代入式(2)得到邊界條件函數(shù)g，并將以上參數(shù)代入 到式(3)得到邊界正定矩陣b。所有運(yùn)算均在軟件 Matlab上實(shí)現(xiàn)。表1、圖1和圖2分別給出了數(shù)值解關(guān)于單元 基函數(shù)個(gè)數(shù)p的誤差。數(shù)值結(jié)果表明數(shù)值解誤差關(guān) 于單元基函數(shù)個(gè)數(shù)p是指數(shù)收斂的。p171921232527誤差/%0.7980.6930.4310.2210.1200.059表1數(shù)值解關(guān)于單元基函數(shù)個(gè)數(shù)P的誤差UWVF00.100.150.20h0.25圖1數(shù)值解關(guān)于單元基函數(shù)個(gè)數(shù)p的誤差7.5一一澄一 一釁7.06.56

16、.05.55.04.51.81.92.02.1lg(p/lg(p)IUWVF-* 一斜率=10.69圖2對(duì)數(shù)下數(shù)值解關(guān)于單元基函數(shù)個(gè)數(shù)p的誤差表2、圖3和圖4分別給出了數(shù)值解關(guān)于網(wǎng)格 步長h的誤差。數(shù)值結(jié)果表明數(shù)值解誤差關(guān)于網(wǎng)格 步長h是代數(shù)收斂的。h1/41/61/81/101/121/14誤差/%0.8800.4700.1500.0490.0190.008表2數(shù)值解關(guān)于網(wǎng)格步長h的誤差3結(jié)語基于間斷有限元方法構(gòu)造數(shù)值求解彈性波方程 的平面波超弱變分方法，實(shí)現(xiàn)了模型參數(shù)輸入、網(wǎng) 格剖分、平面波離散系統(tǒng)的生成及求解、誤差計(jì)算圖3數(shù)值解關(guān)于網(wǎng)格步長h的誤差圖4對(duì)數(shù)下數(shù)值解關(guān)于網(wǎng)格步長h的誤差和輸出功能，具有良好的算法遷移、推廣能力。該算 法在“偏微分方程數(shù)值解”的課程教學(xué)中達(dá)到了以下 目的：掌