張榮培 王迪 蔚喜軍 溫學(xué)兵?

1) (沈陽師范大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院,沈陽 110034)

2) (北京應(yīng)用物理與計算數(shù)學(xué)研究所,北京 100088)

波的傳播往往在復(fù)雜的地質(zhì)結(jié)構(gòu)中進行,如何有效地求解非均勻介質(zhì)中的波動方程一直是研究的熱點.本文將局部間斷 Galekin(local discontinuous Galerkin,LDG)方法引入到數(shù)值求解波動方程中.首先引入輔助變量,將二階波動方程寫成一階偏微分方程組,然后對相應(yīng)的線性化波動方程和伴隨方程構(gòu)造間斷Galerkin格式;為了保證離散格式滿足能量守恒,在單元邊界上選取廣義交替數(shù)值通量,理論證明該方法滿足能量守恒性.在時間離散上,采用指數(shù)積分因子方法,為了提高計算效率,應(yīng)用Krylov子空間方法近似指數(shù)矩陣與向量的乘積.數(shù)值實驗中給出了帶有精確解的算例,驗證了LDG方法的數(shù)值精度和能量守恒性;此外,也考慮了非均勻介質(zhì)和復(fù)雜計算區(qū)域的計算,結(jié)果表明LDG方法適合模擬具有復(fù)雜結(jié)構(gòu)和多尺度結(jié)構(gòu)介質(zhì)中的傳播.

1 引 言

波的傳播是能量傳輸?shù)囊环N基本形式,出現(xiàn)在許多科學(xué)、工程和工業(yè)領(lǐng)域.波動傳播問題的研究對地球科學(xué)、石油工程、電信和國防工業(yè)具有重要意義[1?4].波動方程是描述波傳播的數(shù)學(xué)模型,例如聲波方程、地震波方程等.本文考慮如下二維非均勻介質(zhì)下的二階波動方程.由于所求解的區(qū)域比較復(fù)雜,以及傳播介質(zhì)的復(fù)雜和不均勻性,大多數(shù)波動方程不能精確求解,對于波動方程數(shù)值方法的研究就顯得非常重要.

以往二維波動方程數(shù)值解的計算方法主要是有限差分方法 (finite difference method,FDM)[5,6]、有限體積方法 (finite volume method,FVM)[7]和有限元方法 (finite element method,FEM)[8]等.局部間斷 Galerkin (local discontinuous Galerkin,LDG)方法是Cockburn和舒其望[9]在1998年提出的.與傳統(tǒng)的間斷Galerkin (discontinuous Galerkin,DG)方法[10]比較,LDG 方法通過引入輔助變量將含有高階導(dǎo)數(shù)的微分方程寫成只含有一階導(dǎo)數(shù)的偏微分方程組,然后用DG方法進行空間離散.Chung和Engquist[11]在交錯三角網(wǎng)格上提出了能量守恒的DG格式,但是交錯網(wǎng)格在高維情形格式構(gòu)造復(fù)雜.Chou等[12]發(fā)展了一類能量守恒的LDG方法求解二維非均勻介質(zhì)中的線性二階波動方程,但是他們的工作只限于矩形網(wǎng)格.本文將在三角網(wǎng)格上構(gòu)造能量守恒的LDG格式.首先將二階標(biāo)量波動方程寫成一階速度-應(yīng)力形式,通過選取廣義交替數(shù)值通量,對相應(yīng)的線性化波動方程和伴隨方程構(gòu)造DG格式.對于高維問題,空間離散后得到的是一組規(guī)模非常大的常微分方程組.顯式時間離散格式不需要求解代數(shù)方程組,但是時間步長有嚴格的限制;隱式時間離散允許比較大的時間步長,但是需要求解大規(guī)模代數(shù)方程組.本文采用指數(shù)積分因子方法求解常微分方程組,該方法是在2013年由Nie等[13]提出的求解剛性常微分方程組的有效數(shù)值方法.隱式積分因子方法的發(fā)展分為兩個方向:其一是緊致隱積分因子格式[14,15],該格式將解寫成矩陣形式并在每個方向計算指數(shù)矩陣;其二是Krylov積分因子格式[16],該格式針對指數(shù)矩陣與向量的乘積運算,應(yīng)用Krylov子空間方法近似.

2 局部間斷有限元方法

在二維區(qū)域Ω求解如下二階波動方程:

其中,u 是位移;ρ 是物體的密度;A 是波速,本文取常數(shù);f (x,y,t) 是外部源力項.引入兩個輔助變量:p=A?u,w=ut,(1)式可寫成如下的等價一階偏微分方程組:

邊界條件考慮齊次 Dirichlet邊界:u=0,(x,y)∈?Ω×[0,T].所對應(yīng)的初始條件為 :u(x,y,0)=u0(x,y),p (x,y,0)=p0(x,y),? (x,y)∈ Ω.當(dāng)f=0時,系統(tǒng)(2)和(3)式具有能量守恒性:

下面針對系統(tǒng)(2)和(3)式構(gòu)造LDG格式.首先將二維計算域Ω離散為有限個三角形單元用表示的所有邊界的集合,其中是所有內(nèi)部邊緣的集合.定義LDG近似空間為

矢量函數(shù)q的平均值和跳躍項可以類似定義為

要注意的是,標(biāo)量函數(shù)v的跳躍項 [v]是與法線平行的向量,向量函數(shù)q的跳躍項 [q]是標(biāo)量.由于討論的是齊次 Dirichlet邊界條件,即在?Ω上時,u=0.因此對于這種邊界條件,將[v]和{q} 設(shè)為[v]=vn,{q}=q,其中n是指向外側(cè)的單位法向量.

2.1 LDG方法

在(2)式和(3)式兩側(cè)同時乘以試探函數(shù)vh,qh,并在每個單元上進行積分,通過分部積分可以得到問題(2)式和(3)式的LDG格式,即找到wh∈Vh,ph∈Σh,使得任意試探函數(shù)vh,q,對所有的K∈ Γh滿足

2.2 能量守恒性

在構(gòu)建二維波動方程的數(shù)值方法時,能量守恒通常被考慮在內(nèi).數(shù)值方法是否能夠保持能量守恒可以判斷該方法是否能更好地模擬長時間的波傳播.本節(jié)證明LDG方法可以保持能量守恒.考慮到數(shù)值通量(10)式和(11)式的定義,在所有單元上對(8)式和(9)式求和可得

對(12)式和(13)式左右兩端分別求和,可得

為了得到(12)式和(13)式的能量守恒性,首先證明如下引理.

引理對于所有的試探函數(shù)

是恒成立的.

證明對(15)式左端項進行分部積分可得

利用平均值和跳躍項的定義,把各個單元的總和改寫為邊界形式,通過簡單的計算可得到

將數(shù)值通量(10)式和(11)式代入(15)式右端項得到

即(15)式成立.

利用引理可以得到(8)式和(9)式的能量守恒性.

定理 (能量守恒性)在齊次Dirichlet邊界條件和數(shù)值通量(10)式和(11)式的定義下,系統(tǒng)(8)式和(9)式是能量守恒的,即當(dāng) f=0 時,

證明選取試探函數(shù)vh=wh代入(8)式中可得到

在 (9)式中,選取試探函數(shù)qh=ph,能夠得到

對(19)式和(20)式左右兩邊分別求和,并考慮到(15)式,可得到

3 時間離散

二維波動方程通過LDG方法進行空間離散后,與一維二階波動方相類似,也得到一組非線性常微分方程組,為了節(jié)約這種復(fù)雜代數(shù)方程組求解的計算成本,本文利用指數(shù)積分因子方法來進行時間離散.

首先將LDG格式(8)式和(9)式對于所有單元聯(lián)立,可得到全局非線性常微分方程組的矩陣方程形式,即

其中,W=(w1,w2,···,wJ)T和P=(p1,p2,···pJ)是所有單元上的自由度,wj和pj表示在第j個單元上的自由度,M1是由 3×3 矩陣組成的對角分塊矩陣,M2為6×6矩陣組成的對角分塊矩陣,其逆矩陣容易求解.將(22)式和(23)式進一步寫成如下形式:

在方程(24)式兩端同乘積分因子e-Bt,關(guān)于時間 tn到tn+1進行積分,可得

采用梯形積分公式計算(25)式中的積分,得到二階指數(shù)積分因子格式:

對于高維情況,矩陣指數(shù)的計算將遇到巨大的困難,因為指數(shù)矩陣是大而稠密的.對于這種困難,可以通過使用Krylov子空間方法近似指數(shù)矩陣和向量的乘積來解決[16].

4 數(shù)值算例

首先給出帶有精確解的波動方程測試方法的精度和能量守恒性,然后應(yīng)用LDG方法求解波的傳播問題.第二個算例是求解非均勻介質(zhì)介質(zhì)中波的傳播問題,第三個算例是L形區(qū)域中波的傳播問題.在下面的計算中,考慮線性元,邊界條件為齊次 Dirichlet 邊界,A=I 是單位矩陣,f(x,y,t)=0.

算例1考慮具有常系數(shù)的二維波動方程utt=?2u,(x,y)∈[0,2]×[0,2],初始條件為u(x,y,0)=sin(πx)sin(πy),ut(x,y,0)=0,這個問題的精確解為時間步長取Δt=0.001,表1列出了網(wǎng)格數(shù)分別為16×16,32×32,64×648和1 28×12 時數(shù)值解 wh和p 的誤差和收斂階數(shù),通過表格可以發(fā)現(xiàn)w的精度接近2,p的精度接近1,通常觀察到這種數(shù)值精度的振蕩行為就可以說明所設(shè)計的LDG方法是能量守恒的.

算例2考慮波動方程(1)在單位正方形Ω=[0,1]2上的傳播,將其剖分為64×64個均勻大小的正方形,然后使用一個對角線將每個正方形細分為兩個三角形,如圖1(a)所示.本算例考慮非均勻介質(zhì),密度函數(shù)ρ定義為

表1 數(shù)值解 wh和p 的誤差和收斂階數(shù)Table 1.Error and convergence order of numerical solution w h and p.

圖1 (a) 算例 2 的網(wǎng)格剖分和數(shù)值解 w h 在不同時刻 (b) t=0.2,(c) t=0.3,(d) t=0.4 時的波傳播Fig.1.(a) The triangulation mesh of Example 2;contour plot of solution w h at different time:(b) t=0.2 ;(c) t=0.3 ;(d) t=0.4.

圖2 算例 2 的能量隨時間的演化Fig.2.Energy evolution with time for Example 2.

初始條件w0(x,y)=2exp[-500(x-0.5)2+(y-0.5)2],p0(x,y)=0.圖2 給出了數(shù)值解 wh在時間 t=0.2,0.3,0.4 時的等高線,可以看出,波大約在t=0.3時觸及界面 x=0.65,在那之后,波以更快的速度通過界面.這說明傳播系數(shù)不同會導(dǎo)致波的傳播速度不同.同時也給出離散能量的時間演變圖(圖2).從圖2中可以發(fā)現(xiàn),LDG方法可以保證能量守恒.

算例3考慮圓形波在L形區(qū)域Ω=[0,1]2[0.7,1]2上的傳播,初始條件同算例2.采用非一致三角剖分將其剖分為13578個三角單元,如圖3(a)所示.數(shù)值解 wh在 t=0.1,0.3和0.45 的波形傳播圖在圖3(b)—(d)給出.從圖3可以看出,波在 t=0.3 時刻到達了拐角點.在此之后,波反射得到另外一個圓形波.通過與文獻[11]比較,發(fā)現(xiàn)本文的數(shù)值結(jié)果與其數(shù)值結(jié)果是吻合的.

圖3 (a) 算例 3 的網(wǎng)格剖分和數(shù)值解 w h 在不同時刻 (b) t=0.1,(c) t=0.3,(d) t=0.45 時的波傳播Fig.3.(a) The triangulation mesh of Example 3;contour plot of solution w h at different time:(b) t=0.1 ;(c) t=0.3 ;(d) t=0.45.

5 結(jié) 論

研究了二階波動方程在二維計算區(qū)域的傳播問題,通過選取廣義交替數(shù)值通量,建立了LDG方法并分析了格式的能量守恒性.時間離散上利用指數(shù)積分因子方法.數(shù)值實驗驗證了能量守恒性質(zhì)和最優(yōu)誤差收斂階.下一步的工作考慮將守恒的LDG格式應(yīng)用到求解非線性波動方程,例如Sine-Gordon等方程上.