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

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

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

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

1 引 言

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

以往二維波動(dòng)方程數(shù)值解的計(jì)算方法主要是有限差分方法 (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方法進(jìn)行空間離散.Chung和Engquist[11]在交錯(cuò)三角網(wǎng)格上提出了能量守恒的DG格式,但是交錯(cuò)網(wǎng)格在高維情形格式構(gòu)造復(fù)雜.Chou等[12]發(fā)展了一類能量守恒的LDG方法求解二維非均勻介質(zhì)中的線性二階波動(dòng)方程,但是他們的工作只限于矩形網(wǎng)格.本文將在三角網(wǎng)格上構(gòu)造能量守恒的LDG格式.首先將二階標(biāo)量波動(dòng)方程寫成一階速度-應(yīng)力形式,通過選取廣義交替數(shù)值通量,對(duì)相應(yīng)的線性化波動(dòng)方程和伴隨方程構(gòu)造DG格式.對(duì)于高維問題,空間離散后得到的是一組規(guī)模非常大的常微分方程組.顯式時(shí)間離散格式不需要求解代數(shù)方程組,但是時(shí)間步長(zhǎng)有嚴(yán)格的限制;隱式時(shí)間離散允許比較大的時(shí)間步長(zhǎng),但是需要求解大規(guī)模代數(shù)方程組.本文采用指數(shù)積分因子方法求解常微分方程組,該方法是在2013年由Nie等[13]提出的求解剛性常微分方程組的有效數(shù)值方法.隱式積分因子方法的發(fā)展分為兩個(gè)方向:其一是緊致隱積分因子格式[14,15],該格式將解寫成矩陣形式并在每個(gè)方向計(jì)算指數(shù)矩陣;其二是Krylov積分因子格式[16],該格式針對(duì)指數(shù)矩陣與向量的乘積運(yùn)算,應(yīng)用Krylov子空間方法近似.

2 局部間斷有限元方法

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

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

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

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

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

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

2.1 LDG方法

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

2.2 能量守恒性

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

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

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

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

是恒成立的.

證明對(duì)(15)式左端項(xiàng)進(jìn)行分部積分可得

利用平均值和跳躍項(xiàng)的定義,把各個(gè)單元的總和改寫為邊界形式,通過簡(jiǎn)單的計(jì)算可得到

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

即(15)式成立.

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

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

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

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

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

3 時(shí)間離散

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

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

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

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

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

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

4 數(shù)值算例

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

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

算例2考慮波動(dòng)方程(1)在單位正方形Ω=[0,1]2上的傳播,將其剖分為64×64個(gè)均勻大小的正方形,然后使用一個(gè)對(duì)角線將每個(gè)正方形細(xì)分為兩個(gè)三角形,如圖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 在不同時(shí)刻 (b) t=0.2,(c) t=0.3,(d) t=0.4 時(shí)的波傳播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 的能量隨時(shí)間的演化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在時(shí)間 t=0.2,0.3,0.4 時(shí)的等高線,可以看出,波大約在t=0.3時(shí)觸及界面 x=0.65,在那之后,波以更快的速度通過界面.這說明傳播系數(shù)不同會(huì)導(dǎo)致波的傳播速度不同.同時(shí)也給出離散能量的時(shí)間演變圖(圖2).從圖2中可以發(fā)現(xiàn),LDG方法可以保證能量守恒.

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

圖3 (a) 算例 3 的網(wǎng)格剖分和數(shù)值解 w h 在不同時(shí)刻 (b) t=0.1,(c) t=0.3,(d) t=0.45 時(shí)的波傳播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é) 論

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