基于Matlab 求解時(shí)諧復(fù)系數(shù)彈性波方程的超弱變分方法1)

袁 龍 趙茂先

(山東科技大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，山東青島266590)

“偏微分方程數(shù)值解” 課程是數(shù)學(xué)與應(yīng)用數(shù)學(xué)、信息與計(jì)算科學(xué)等專(zhuān)業(yè)的核心課程。該課程針對(duì)刻畫(huà)物理模型的數(shù)學(xué)方程，利用有限元、有限差分法、有限體積法和邊界元法等方法數(shù)值求解方程，并理論分析相應(yīng)數(shù)值解的誤差收斂階。該課程涉及知識(shí)面廣泛，不但要講授變分方法、索伯列夫空間、離散空間(如多項(xiàng)式空間)、線性方程組的迭代求解及數(shù)值逼近，還要培養(yǎng)學(xué)生的編程驗(yàn)證能力[1-3]。因此，學(xué)生系統(tǒng)地掌握該課程的知識(shí)結(jié)構(gòu)是至關(guān)重要的。

本文基于間斷有限元方法數(shù)值求解描述一大類(lèi)物理問(wèn)題模型(比如聲波散射、彈性波散射問(wèn)題) 的彈性波方程。在數(shù)值求解彈性波方程的發(fā)展過(guò)程中，有限元[4-5]、邊界元法[6-7]和有限差分方法[8-9]起到了很重要的作用。彈性波方程解析解的一個(gè)重要參數(shù)是波數(shù)。它描述了解析解的振蕩行為。波數(shù)越大，振蕩越頻繁。這個(gè)特點(diǎn)在有限元分析中可以通過(guò)經(jīng)驗(yàn)法則來(lái)刻畫(huà)——單個(gè)波長(zhǎng)內(nèi)的網(wǎng)格點(diǎn)數(shù)。因而由此經(jīng)驗(yàn)法則很難獲得大波數(shù)情形下令人滿(mǎn)意的數(shù)值結(jié)果。為了克服此缺點(diǎn)，很多改進(jìn)方法[10-12]被陸續(xù)提出。最近用來(lái)求解中高頻彈性波方程的平面波方法[13-15]在工業(yè)界和數(shù)學(xué)界比較流行。平面波方法的核心是選取滿(mǎn)足無(wú)約束方程的精確解作為逼近基函數(shù)，如平面波函數(shù)、貝塞爾函數(shù)、格林函數(shù)等，并且未知量是定義在剖分單元上的。大量的數(shù)值實(shí)驗(yàn)表明，相比其他離散化方法，在滿(mǎn)足相同精度的條件下，平面波方法使用了相對(duì)較少的自由度；即，在選取相同自由度的條件下，平面波離散化方法具有較高的精度。

為了讓學(xué)生直觀感受平面波方法數(shù)值求解偏微分方程的過(guò)程，本文從方程組本身出發(fā)，分別構(gòu)造滿(mǎn)足彈性波方程和其對(duì)偶方程的測(cè)試函數(shù)空間和檢驗(yàn)函數(shù)空間，并基于間斷有限元思想，結(jié)合分部積分公式，推導(dǎo)出復(fù)波數(shù)情形的超弱變分形式；分別構(gòu)造二維空間下的平面波檢驗(yàn)離散函數(shù)空間和測(cè)試離散函數(shù)空間，進(jìn)一步得到等價(jià)的變分形式。數(shù)值實(shí)驗(yàn)驗(yàn)證了新方法的數(shù)值解具有較高的精度。

1 控制方程與數(shù)值方法

1.1 彈性波方程

考慮邊值問(wèn)題[16]

及最低階吸收邊界條件

其中，拉梅常數(shù)λ和μ借助于泊松比ν和楊氏模量E表達(dá)

材料密度ρ>0 是常數(shù)，ω >0 是角頻率，n是邊界的單位外法向，i 是虛部單位，截?cái)嗨阕覶n(u) 為

定義縱波和橫波的波速分別為

定義在邊界上的正定矩陣σ為

其中s是邊界的單位切向量，|CP|與|CS|分別表示CP與CS的模，·T代表·的轉(zhuǎn)置。假定求解區(qū)域?是有界多面體區(qū)域。Th表示區(qū)域?的三角剖分并滿(mǎn)足“形正則”和“擬一致”假設(shè)[17]。定義為簡(jiǎn)單起見(jiàn)，我們只考慮二維區(qū)域模型，但方法可以直接推廣到三維空間。

1.2 彈性波方程的變分形式

對(duì)于每一個(gè)剖分單元?k，令H1(?k) 為標(biāo)準(zhǔn)的向量索伯列夫空間，并定義分片間斷Trefftz 變分空間[18]

和區(qū)域間斷Trefftz 測(cè)試空間

對(duì)于每個(gè)單元?k，定義跡算子

利用等距引理[19-20]，原問(wèn)題(1)(2)等價(jià)于變分形式：尋找u ∈V(Th) 滿(mǎn)足

其中半雙線性型

右端項(xiàng)為

1.3 平面波離散化

本節(jié)構(gòu)造間斷 Trefftz 變分空間V(Th) 的平面波離散空間Vp(Th)。首先考慮沿單位方向d傳播的時(shí)諧彈性平面波，其可分解為縱波與橫波之和

其中κp=ω/Cp,κs=ω/Cs，α,β是系數(shù)，d ·e= 0。縱波vp=αdexp(iκpd·x) 和橫波vs=βeexp(iκsd·x)。簡(jiǎn)單計(jì)算可知，縱波和橫波均滿(mǎn)足原始方程 (1)，且?×vp= 0,?·vs= 0。理論分析和實(shí)際應(yīng)用表明，借助于選取單位圓周上分布均勻的一簇互不相同的彈性波傳播方向

我們得到定義在每個(gè)單元?k上的2p個(gè)復(fù)平面波基函數(shù)

記Qq(q= 2p) 表示單元?k上的上述 2p個(gè)復(fù)平面波基函數(shù)張成的空間，則平面波離散空間Vp(Th) 為

定義解析解uk=u|?k的平面波逼近uh,k

其中αk,l,βk,l(k=1,2,··· ,N;l=1,2,··· ,p) 為待求解的系數(shù)。對(duì)應(yīng)于連續(xù)變分形式(4) 的離散形式為：尋找uh ∈Vp(Th) 滿(mǎn)足

我們利用直接法，即高斯消去法[21]，求解上述離散線性方程組。

1.4 算法實(shí)現(xiàn)

步驟 1. 初始化輸入?yún)?shù)：計(jì)算區(qū)域?,λ,μ,ω,ρ,σ；

步驟2. 剖分區(qū)域?形成三角網(wǎng)格Th；

步驟3. 基于變分形式 (5) 的間斷有限元離散，形成剛度矩陣A及右端項(xiàng)b；

步驟4. 利用線性解法器求解線性方程組Ax=b；

步驟5. 利用解向量構(gòu)造有限元解uh并計(jì)算誤差；

步驟6. 分析誤差關(guān)于p及h的變化。

2 數(shù)值算例與結(jié)果分析

取ξ1= 1,ξ2=?1，計(jì)算區(qū)域?=[0,1]2。取Rayleigh 波作為邊值問(wèn)題(1) 和(2) 的解析解[22]

取E= 2.0×1011,ν= 0.3,ρ= 7800。將解析解u代入式(2)得到邊界條件函數(shù)g，并將以上參數(shù)代入到式 (3) 得到邊界正定矩陣σ。所有運(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ù)收斂的。

表1 數(shù)值解關(guān)于單元基函數(shù)個(gè)數(shù)p 的誤差

圖1 數(shù)值解關(guān)于單元基函數(shù)個(gè)數(shù)p 的誤差

圖2 對(duì)數(shù)下數(shù)值解關(guān)于單元基函數(shù)個(gè)數(shù)p 的誤差

表2、圖 3 和圖 4 分別給出了數(shù)值解關(guān)于網(wǎng)格步長(zhǎng)h的誤差。數(shù)值結(jié)果表明數(shù)值解誤差關(guān)于網(wǎng)格步長(zhǎng)h是代數(shù)收斂的。

表 2 數(shù)值解關(guān)于網(wǎng)格步長(zhǎng)h 的誤差

3 結(jié)語(yǔ)

基于間斷有限元方法構(gòu)造數(shù)值求解彈性波方程的平面波超弱變分方法，實(shí)現(xiàn)了模型參數(shù)輸入、網(wǎng)格剖分、平面波離散系統(tǒng)的生成及求解、誤差計(jì)算和輸出功能，具有良好的算法遷移、推廣能力。該算法在“偏微分方程數(shù)值解”的課程教學(xué)中達(dá)到了以下目的：

圖3 數(shù)值解關(guān)于網(wǎng)格步長(zhǎng)h 的誤差

圖4 對(duì)數(shù)下數(shù)值解關(guān)于網(wǎng)格步長(zhǎng)h 的誤差

(1)掌握間斷有限元方法的思想，能夠區(qū)分其與其他離散化方法的不同，使學(xué)生對(duì)算法加深理解。

(2)借助于代碼實(shí)現(xiàn)，掌握間斷有限元離散化方法的實(shí)現(xiàn)細(xì)節(jié)，特別是單剛合成總剛、線性方程組的求解等。

(3)借助于數(shù)值解誤差的圖形分析，使學(xué)生掌握數(shù)值解誤差的收斂階驗(yàn)證方法。

將該數(shù)值離散化方法應(yīng)用于課程教學(xué)和實(shí)驗(yàn)教學(xué)中，對(duì)學(xué)生深入理解基礎(chǔ)理論知識(shí)和掌握相應(yīng)的實(shí)驗(yàn)方法有積極的作用，可以激發(fā)學(xué)生的科研興趣和創(chuàng)新意識(shí)，提高學(xué)生的創(chuàng)新能力。目前，在此教學(xué)實(shí)驗(yàn)的基礎(chǔ)上，我校學(xué)生正著力于實(shí)現(xiàn)復(fù)雜的三維區(qū)域網(wǎng)格剖分，取得了較好的實(shí)驗(yàn)成果。