袁 龍 趙茂先

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

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

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

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

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

1.1 彈性波方程

考慮邊值問題[16]

及最低階吸收邊界條件

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

材料密度ρ>0 是常數(shù)，ω >0 是角頻率，n是邊界的單位外法向，i 是虛部單位，截斷算子Tn(u) 為

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

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

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

1.2 彈性波方程的變分形式

對于每一個剖分單元?k，令H1(?k) 為標準的向量索伯列夫空間，并定義分片間斷Trefftz 變分空間[18]

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

對于每個單元?k，定義跡算子

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

其中半雙線性型

右端項為

1.3 平面波離散化

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

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

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

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

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

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

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

1.4 算法實現(xiàn)

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

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

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

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

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

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

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

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

取E= 2.0×1011,ν= 0.3,ρ= 7800。將解析解u代入式(2)得到邊界條件函數(shù)g，并將以上參數(shù)代入到式 (3) 得到邊界正定矩陣σ。所有運算均在軟件Matlab 上實現(xiàn)。

表 1、圖 1 和圖 2 分別給出了數(shù)值解關(guān)于單元基函數(shù)個數(shù)p的誤差。數(shù)值結(jié)果表明數(shù)值解誤差關(guān)于單元基函數(shù)個數(shù)p是指數(shù)收斂的。

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

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

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

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

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

3 結(jié)語

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

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

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

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

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

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

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