三角網(wǎng)格間斷有限元法彈性波模擬精度分析

韓德超 劉衛(wèi)華 司文朋*

(①中國石化石油物探技術(shù)研究院，江蘇南京 211103；②中國石化地球物理重點(diǎn)實(shí)驗(yàn)室，江蘇南京 211103)

0 引言

地震波數(shù)值模擬對認(rèn)識地震波傳播規(guī)律和分析干擾波形成機(jī)理具有重要意義，因此常被用于優(yōu)化觀測系統(tǒng)和指導(dǎo)處理解釋[1-5]。此外，地震波數(shù)值模擬也是逆時(shí)偏移和波形反演的重要內(nèi)核[6]。目前已有許多基于波動(dòng)方程的數(shù)值模擬方法，如：有限差分法(FDM)[7-10]、偽譜法(PSM)[11-12]、有限體積法(FVM)[13-14]、無網(wǎng)格法[15]、有限元法(FEM)[16-18]，譜元法(SEM)[19-23]和間斷Galerkin有限元方法(DGFEM)[24-27]。FDM因容易編程實(shí)現(xiàn)、易于并行等特點(diǎn)，目前應(yīng)用最廣泛。但是FDM對起伏地表的適應(yīng)性有限，且準(zhǔn)確模擬地表自由邊界較為困難。PSM同樣也存在邊界條件難以準(zhǔn)確模擬的問題，且三維計(jì)算需要全局通信，難于實(shí)現(xiàn)并行。FVM和FEM都可以使用非結(jié)構(gòu)化的三角形(三維是四面體)網(wǎng)格以適應(yīng)復(fù)雜的地表結(jié)構(gòu)，且天然滿足自由邊界條件，但是FVM的高階計(jì)算難于實(shí)現(xiàn)，較少應(yīng)用于地震波模擬。FEM則需要存儲大型稀疏矩陣并求逆，對計(jì)算資源要求極高。SEM和DGFEM是對FEM的改進(jìn)。SEM使用GLL(Gauss Lobatto Legendre)插值點(diǎn)得到對角形式的質(zhì)量矩陣，從而避免了矩陣求逆過程。但是經(jīng)典的SEM使用的是四邊形(三維是六面體)網(wǎng)格，對復(fù)雜地表的刻畫缺乏靈活性。DGFEM在單元內(nèi)部使用有限元方法處理，而在單元邊界上使用FVM的數(shù)值流通量的方式。DGFEM集合了有限元方法的網(wǎng)格靈活性和高精度，同時(shí)又可以逐元求解，降低了對計(jì)算資源的消耗且便于并行計(jì)算。

數(shù)值穩(wěn)定性、數(shù)值頻散以及耗散問題是研究波動(dòng)方程數(shù)值解的重要內(nèi)容。有限元方法因使用的網(wǎng)格形態(tài)比較靈活，所以一般需構(gòu)建周期性網(wǎng)格進(jìn)行分析。四邊形網(wǎng)格的模擬精度分析因?yàn)橹芷诰W(wǎng)格構(gòu)建容易、基函數(shù)可以解耦等原因，已有大量研究[16，28-32]。而對于三角形網(wǎng)格的數(shù)值精度研究較少，劉少林等[33]和曹丹平等[34]分析了經(jīng)典FEM在不同形態(tài)的周期性三角形網(wǎng)格中的聲波方程頻散和穩(wěn)定性條件，也有少量的關(guān)于三角形網(wǎng)格的SEM頻散分析研究[35-37]。對于四邊形網(wǎng)格的DGFEM，De Basabe等[38]分析了不同基函數(shù)對數(shù)值精度的影響；賀茜君等[39]分析了Runge-Kutta(RK)時(shí)間格式的頻散和穩(wěn)定性；Meng等[40]分析了不同基函數(shù)對頻散和穩(wěn)定性的影響。對于三角形網(wǎng)格的DGFEM，Hu等[41]分析了半離散聲波方程的頻散和耗散，但未考慮時(shí)間格式；K?ser等[42]從實(shí)驗(yàn)角度分析了彈性波任意高階導(dǎo)數(shù)法(Arbitrary high-order derivatives，ADER)時(shí)間格式的數(shù)值收斂率，但未分析頻散和耗散；Delcourte等[43]分析了蛙跳(Frog-leap，F(xiàn)L)時(shí)間格式的基于中心數(shù)值流的收斂率；He等[44]分析了聲波方程RK和ADER時(shí)間格式的數(shù)值精度。

本文針對彈性波的數(shù)值模擬精度問題，通過構(gòu)造任意等腰三角形單元組成的周期性網(wǎng)格，可以方便地分析各種形態(tài)的三角形單元對數(shù)值模擬精度的影響。從理論和數(shù)值計(jì)算上分析了基于局部Lax-Friedrichs數(shù)值流和3階總變差不增的Runge-Kutta(簡稱TVD RK)時(shí)間格式的三角形網(wǎng)格DGFEM彈性波模擬的數(shù)值頻散、耗散和穩(wěn)定性，可為基于三角形網(wǎng)格的DGFEM彈性波模擬提供參數(shù)選取指導(dǎo)。

1 數(shù)值模擬方法

二維各向同性介質(zhì)中的彈性波速度—應(yīng)力方程為[42]

(1)

式中：λ、μ為拉梅常數(shù)；ρ是密度；σxx和σzz為正應(yīng)力，σxz為剪切應(yīng)力；vx和vz分別為質(zhì)點(diǎn)速度的x和z分量；S1～S5表示外力震源。將式(1)改寫成一階線性雙曲偏微分方程矩陣形式，有

(2)

式中：u=(σxx，σzz，σxz，vx，vz)T；P、B是式(1)中彈性參數(shù)的矩陣形式。

1.1 空間離散

DGFEM的空間離散需要將計(jì)算空間劃分成不重疊的子區(qū)域，即單元。本文討論二維三角形網(wǎng)格，直接給出式(2)的DGFEM空間離散格式[36]

(3)

(4)

1.2 時(shí)間離散

間斷有限元求解有多種時(shí)間離散方法，如ADER、LF和RK法。ADER法將時(shí)間導(dǎo)數(shù)轉(zhuǎn)化到空間導(dǎo)數(shù)的計(jì)算上，可以達(dá)到時(shí)間和空間任意精度，但是實(shí)現(xiàn)復(fù)雜[42]；LF時(shí)間格式通常和中心數(shù)值流一起使用[43，46]，且可以使用更大的時(shí)間步長從而提高計(jì)算效率，但是常用的LF時(shí)間格式精度只有二階。本文使用RK法。

式(3)可簡寫為關(guān)于時(shí)間的常微分方程

(5)

式中忽略了震源項(xiàng)，L是空間離散的線性系統(tǒng)。常用的3階TVD RK時(shí)間離散格式[25，39]為

(6)

2 數(shù)值穩(wěn)定性及精度分析

2.1 基于三角形網(wǎng)格的理論分析方法

穩(wěn)定性和數(shù)值精度的理論分析通常使用Von Neumann方法：假設(shè)研究區(qū)域?yàn)闊o界各向同性的均勻介質(zhì)且無外力震源，分析平面波在研究區(qū)域的傳播問題。

首先假設(shè)計(jì)算空間剖分為以頂角為θ、腰長為h的等腰三角形為基本單元的周期性網(wǎng)格(圖1)，其中單元2-i是由單元1-i旋轉(zhuǎn)所得，除此之外，單元完全相同。因此將單元1-1和單元2-1稱為母單元，其他單元稱為子單元。為了簡化，將相鄰的兩類三角形組合在一起分析[41]，即圖1中單元1-1和單元2-1作為整體進(jìn)行分析。假設(shè)母單元中的邊界順序如圖1b所示，子單元邊界順序同母單元。

圖1 周期性三角形網(wǎng)格剖分示意圖

考慮平面簡諧波

Aexp(-iωt)exp(ikx)=A(t)exp(ikx)

(7)

則各子單元內(nèi)的平面波與母單元平面波的相位關(guān)系如表1所示。假設(shè)波的傳播方向與x軸夾角為γ，則kxh=2πδcosγ，kzh=2πδsinγ，其中k為波數(shù)矢量，δ=h/W是采樣率，W是波長。

表1 各子單元與母單元的平面波相位關(guān)系

忽略震源項(xiàng)，單元1-2和單元2-1的間斷Galerkin有限元法的半離散方程(式(3))可以寫為

(8)

式中

(9)

(10)

其中的子矩陣Γ11、Γ12、Γ21、Γ22表達(dá)式見附錄A。

2.2 穩(wěn)定性分析

模擬時(shí)在網(wǎng)格尺寸h和介質(zhì)速度cP已定的情況下，時(shí)間步長需滿足Δt≤αmaxh/cP，其中α為Courant數(shù)，αmax=ΔtmaxcP/h是最大Courant數(shù)。DGFEM的3階TVD RK方法可以寫成[33]

(11)

(12)

式中ωh是與三角形腰長h有關(guān)的數(shù)值圓頻率?？梢钥闯觯?exp(-iωhΔt)是矩陣G的特征值。為了使模擬穩(wěn)定，矩陣G的譜半徑必須小于等于1，即|β|≤1[34]。需要注意的是，求解式(13)的特征值，系統(tǒng)會得到20N個(gè)特征值，要求每個(gè)特征值都滿足|β|≤1，可用迭代方法求得最大Courant數(shù)αmax[40]。αmax越小說明模擬過程的穩(wěn)定性越差，反之穩(wěn)定性越強(qiáng)。

αmax隨θ的變化曲線如圖2a所示，可見，當(dāng)單元過于狹長或扁平時(shí)模擬的穩(wěn)定性較差。圖2b為一、二、三階單元αmax隨r/h(r為內(nèi)切圓半徑)變化曲線，可知，各階單元的αmax與r/h滿足線性關(guān)系，其斜率分別為0.6137、0.3535、0.2683，約等于2/(2p+1)，其中p為單元階次。由Δt=αh/cP可知，Δt與三角形內(nèi)切圓半徑r也滿足線性關(guān)系，即Δt∝2/(2p+1)r/cP。

圖2 不同階單元αmax隨θ(a)和r/h(b)的變化曲線

2.3 頻散和耗散分析

(13)

將ωh代入式(12)左側(cè)，可得

exp(ωiΔt)[cos(ωrΔt)-i sin(ωrΔt)]=βr+iβi

(14)

進(jìn)一步可得

(15)

則數(shù)值頻散系數(shù)d和耗散系數(shù)d′可分別表示為

(16)

(17)

當(dāng)d大于1時(shí)，相速度大于介質(zhì)速度，相反則說明相速度小于介質(zhì)速度。d越接近1，模擬精度越高。由于P波和S波具有相同的頻率、不同的空間采樣率，因此分別對P波和S波進(jìn)行分析。

當(dāng)Courant數(shù)α=0.03、平面波傳播方向角γ=0°時(shí)，計(jì)算θ分別為30°、60°、90°、120°時(shí)的二階和三階單元的P波和S波頻散系數(shù)d隨δ的變化曲線(圖3)。由圖3可以看出，二階單元在δ小于0.4時(shí)，P、S波頻散系數(shù)的相對誤差小于0.9%，三階單元的相對誤差則小于0.05%，體現(xiàn)了DGFEM的弱頻散特性。二階單元一個(gè)波長內(nèi)有2～3個(gè)網(wǎng)格即可滿足頻散要求。此外，由二階單元的頻散結(jié)果可以看出，在γ=0°方向上傳播的平面波，θ為60°時(shí)頻散最弱，其次是30°、90°、120°；S波的頻散強(qiáng)弱相對關(guān)系與P波一致。不同形態(tài)(頂角θ不同)三階單元的縱波頻散相對關(guān)系與二階單元一致。

圖3 不同形態(tài)二階(a)和三階(b)單元的的P波和S波頻散曲線

圖4是一階單元和二階單元頻散系數(shù)隨傳播方向的變化曲線，由圖中可以看出，當(dāng)單元為等邊三角形時(shí)，頻散系數(shù)隨傳播方向變化的周期為60°，且沿單元邊界方向傳播時(shí)頻散最小，而垂直網(wǎng)格線傳播時(shí)頻散最大。而經(jīng)典有限元方法在等邊三角形網(wǎng)格中基本沒有數(shù)值各向異性。二者的區(qū)別可能是由DGFEM的數(shù)值流形式引起的。其他形態(tài)三角形單元的頻散系數(shù)隨傳播方向的變化以180°為周期。等邊三角形單元的頻散曲線更接近圓形，數(shù)值各向異性最弱，而其他形態(tài)的三角形單元的頻散多呈條帶狀，數(shù)值各向異性更強(qiáng)。在圖4b中，當(dāng)θ=30°時(shí)頻散曲線在γ=15°和195°處頻散最大，即垂直于等腰三角形底邊的方向，最小頻散的方向?yàn)檠氐走叿较?；?dāng)θ=90°時(shí)，最大頻散的方向?yàn)棣?135°和225°，當(dāng)θ=120°時(shí)，最大頻散的方向?yàn)棣?150°和330°，即沿底邊的傳播方向頻散最大，垂直于底邊方向頻散最小。

“人設(shè)”先行，架空了歷史，扭曲了情節(jié)。熱播的《延禧攻略》中，主角魏瓔珞的“人設(shè)”可謂鮮明至極：她聰明果敢、為人剛強(qiáng)，見招拆招、快意恩仇，以披靡之勢走向人生巔峰。“底層逆襲”“人人都愛我”等現(xiàn)實(shí)中無法實(shí)現(xiàn)的白日夢在其中得到代償，使得這樣的“人設(shè)”很是討喜。網(wǎng)劇情節(jié)有所夸張無可厚非，但如果為了凸顯“人設(shè)”而不顧歷史常識，違背起碼的生活邏輯，生造出在現(xiàn)實(shí)中不可能發(fā)生的情節(jié)，“人設(shè)”就失去了堅(jiān)實(shí)的支撐與豐滿的內(nèi)里，雖滿足了受眾一時(shí)的感官刺激，卻成了虛假而空洞的存在。

圖4 不同形態(tài)單元縱波頻散系數(shù)隨傳播方向角γ的變化曲線

θ分別為30°、60°、90°、120°時(shí)，計(jì)算二階和三階單元的網(wǎng)格P波和S波的耗散系數(shù)隨采樣率的變化曲線(圖5)，其中α=0.03。為了明顯區(qū)分曲線，圖5是模擬了200個(gè)時(shí)間步長后的耗散曲線(即exp(200ωiΔt))。由圖可見，隨著δ的增加，數(shù)值耗散越強(qiáng)；對于相同的δ，S波的耗散略小于P波，通常情況下，同一介質(zhì)內(nèi)的S波采樣率大于P波的采樣率；另外，隨著單元階數(shù)的增大，耗散的影響減弱。通常地震波需要模擬幾千甚至幾萬個(gè)時(shí)間步長，數(shù)值耗散將對模擬結(jié)果產(chǎn)生較大影響。

圖5 不同形態(tài)二階(a)和三階(b)單元的的P波和S波耗散曲線

圖6為θ不同時(shí)，一階單元網(wǎng)格的exp(ωiΔt)隨傳播方向的變化曲線，其中α=0.03、δP=0.11、δS=0.20。由圖可見，耗散曲線的對稱性與頻散曲線(圖4)相同。單元頂角θ=60°時(shí)，耗散系數(shù)整體上更接近1。

圖6 不同形態(tài)一階單元的耗散系數(shù)隨傳播方向角γ的變化曲線

表2梳理了上述周期性網(wǎng)格的理論分析結(jié)果。表中數(shù)值各向異性是基于已比較過的單元形態(tài)。

表2 不同形態(tài)三角形單元的最大αmax及數(shù)值各向異性特征

3 數(shù)值實(shí)驗(yàn)

首先通過模擬二維均勻各向同性介質(zhì)平面波單頻成分分析不同形態(tài)單元的誤差。不考慮震源項(xiàng)，均勻各向同性介質(zhì)中彈性平面波應(yīng)力和速度分量的解析解表達(dá)式分別為

(18)

(19)

上式每個(gè)分量都包含一個(gè)沿著(nx，nz)T方向傳播的平面P波和S波。取t=0時(shí)刻的解析值作為數(shù)值模擬的初始條件。

分別用θ為30°、60°、90°的周期性網(wǎng)格進(jìn)行均勻模型模擬，并對每個(gè)θ取其數(shù)值各向異性的兩個(gè)對稱軸的方向作為平面波傳播方向，即θ=30°時(shí)，取γ=15°或105°；θ=60°時(shí)，取γ=0°或30°；θ=90°時(shí)，取γ=45°或135°。等腰三角形單元的腰長h分別設(shè)置為10.0、12.5、20.0、25.0m。均勻模型尺寸為5000m×5000m，cP=4000m/s、cS=2000m/s、ρ=2000kg/m3；模擬頻率f0=40Hz，Δt=2ms，模擬時(shí)長tmax=250ms。

圖7為θ=90°、γ=45°時(shí)，tmax時(shí)刻的模擬波場與解析解的對比。由于沒有使用吸收和周期性邊界，所以取模型中沒有受到邊界反射干擾的區(qū)域。由圖7可以看出，使用一階單元且h=25.0m時(shí)的模擬波場與解析解在能量上有很大差異，模擬誤差主要來源于是數(shù)值耗散(縱波空間采樣率為0.25，橫波采樣率為0.5)。使用二階單元，h=10.0m時(shí)的模擬波場與解析解接近(縱、橫波的空間采樣率分別為0.1和0.2)。圖8為圖7深度為2000m處的波場，可見：二階單元h=10.0m和h=12.5m兩種網(wǎng)格的計(jì)算結(jié)果與解析解對應(yīng)較好；當(dāng)使用一階單元時(shí)，即使空間采樣率取0.1，也會產(chǎn)生很強(qiáng)的數(shù)值耗散。

圖7 θ=90°、γ=45°時(shí)tmax時(shí)刻的不同網(wǎng)格模擬的40Hz單頻波場與解析解的對比

圖8 θ=90°、γ=45°時(shí)tmax時(shí)刻、深度2000m處不同網(wǎng)格模擬的波場與解析解的對比

(20)

式中Ω為選取的計(jì)算區(qū)域。具體地，當(dāng)θ=30°時(shí)，Ω為4600m≤x≤4900m，1200m≤z≤1400m；當(dāng)θ=60°時(shí)，Ω為3500m≤x≤4500m，2000m≤z≤3000m；當(dāng)θ=90°時(shí)，Ω為2000m≤x≤3000m，2000m≤z≤3000m。

圖9、表3給出了不同θ的一階和二階單元網(wǎng)格在不同傳播方向上的模擬誤差。可以看出，由于誤差主要受數(shù)值耗散的影響，對于相同θ，兩個(gè)對稱軸方向上的誤差相對關(guān)系與耗散隨傳播方向的變化曲線(圖5和圖6)一致；對于相同階次的單元，其網(wǎng)格尺寸與誤差在雙對數(shù)坐標(biāo)系下呈線性關(guān)系：一階單元時(shí)，斜率約為2；二階單元時(shí)，斜率約為4；單元從一階提高到二階，網(wǎng)格越小，誤差減小越明顯。

表3 不同形態(tài)的一階和二階單元在不同傳播方向上的誤差統(tǒng)計(jì)

同時(shí)，計(jì)算了在一般性非規(guī)則網(wǎng)格中的3階TVD RK時(shí)間格式的DGFEM的誤差(圖9、表4)。圖10是誤差統(tǒng)計(jì)區(qū)域內(nèi)的非結(jié)構(gòu)性網(wǎng)格，呈無規(guī)律排列。圖9中非規(guī)則網(wǎng)格誤差曲線同樣表現(xiàn)為線性，而且在45°和135°兩個(gè)方向上沒有出現(xiàn)明顯的數(shù)值各向異性。

表4 非結(jié)構(gòu)網(wǎng)格中一階和二階單元在45°和135° 方向上的誤差

圖9 不同θ的時(shí)一階和二階單元網(wǎng)格在不同傳播方向上的模擬誤差隨網(wǎng)格尺寸的變化曲線藍(lán)色、紅色和紫色線分別代表頂角為60°、30°、90°的周期性網(wǎng)格結(jié)果，綠色為非規(guī)則網(wǎng)格結(jié)果

圖10 模擬使用的非結(jié)構(gòu)性網(wǎng)格

應(yīng)用上述均勻各向同性介質(zhì)模型進(jìn)行點(diǎn)源激發(fā)模擬，以更加直觀地說明網(wǎng)格對DGFEM的波場特征的影響。震源為主頻40Hz的Ricker子波，時(shí)間步長Δt=8ms。圖11為h取不同值時(shí)一階直角三角形周期性網(wǎng)格與一般非結(jié)構(gòu)化網(wǎng)格模擬得到的0.3s時(shí)刻的vz波場快照對比。圖12為二階直角三角形周期性網(wǎng)格與一般非結(jié)構(gòu)化網(wǎng)格模擬得到的0.3s時(shí)刻的vz波場快照對比。由圖11可以看出，當(dāng)網(wǎng)格階數(shù)為一時(shí)，隨著網(wǎng)格尺寸的增加，兩種網(wǎng)格的波場能量總體都有所下降。由h=20.0m的直角三角形周期性網(wǎng)格模擬的波場快照(圖11a右和圖12a右)可以看出，在135°和315°方向，縱波和橫波的波形變寬，旁瓣更加明顯，能量比45°方向弱很多，這是數(shù)值耗散所致。而在圖11b右和圖12b右所示的h=20.0m的一般非結(jié)構(gòu)化網(wǎng)格模擬的性網(wǎng)格的波場快照上，觀察不到明顯的方向差異，說明非結(jié)構(gòu)數(shù)值各向異性不明顯。圖12b右中有部分雜亂的波場，這是由于網(wǎng)格尺寸增大導(dǎo)致的點(diǎn)加載的震源噪聲。通?？梢詫⒄鹪匆愿咚狗植嫉姆绞郊虞d在震源點(diǎn)附近的一定范圍內(nèi)以減弱這一噪聲[27]。

圖11 h取不同值時(shí)一階直角三角形周期性網(wǎng)格(a)與一般非結(jié)構(gòu)化網(wǎng)格(b)模擬的0.3s時(shí)刻的vz波場快照對比

圖12 h取不同值時(shí)二階直角三角形周期性網(wǎng)格(a)與一般非結(jié)構(gòu)化網(wǎng)格(b)模擬的0.3s時(shí)刻的vz波場快照對比

4 結(jié)論

本文通過構(gòu)造周期性三角形網(wǎng)格模型，分析了基于Runge-Kutta的間斷Galerkin有限元方法的穩(wěn)定性條件、數(shù)值頻散以及數(shù)值耗散，得到以下結(jié)論和認(rèn)識。

(1)DGFEM穩(wěn)定模擬的最大時(shí)間步長Δtmax與單元內(nèi)切圓半徑存在線性關(guān)系。分析結(jié)果表明，等邊三角形的穩(wěn)定性條件最寬松。

(2)從理論和實(shí)驗(yàn)角度說明了基于局部Lax-Friedrichs數(shù)值流的DGFEM的弱頻散和強(qiáng)耗散特性，且在不同傳播方向上，二者都表現(xiàn)出方向各向異性。數(shù)值各向異性呈雙軸對稱性。等邊三角形網(wǎng)格的數(shù)值各向異性最弱，而直角三角形和鈍角三角形網(wǎng)格表現(xiàn)出的數(shù)值各向異性較強(qiáng)。一般的非結(jié)構(gòu)性網(wǎng)格由于網(wǎng)格排列沒有一定規(guī)律，因此沒有表現(xiàn)出明顯的數(shù)值各向異性，說明在網(wǎng)格生成時(shí)，應(yīng)盡量避免鈍角三角形和呈大面積排列的直角三角形，以減小數(shù)值各向異性的影響。

(3)誤差分析顯示，模擬誤差與網(wǎng)格尺寸在雙對數(shù)坐標(biāo)系下呈線性關(guān)系。數(shù)值實(shí)驗(yàn)表明，模擬誤差主要來源是基于局部Lax-Friedrichs數(shù)值流的數(shù)值耗散。實(shí)驗(yàn)中，使用一階單元模擬的波場均有較強(qiáng)耗散。而選取二階單元，一個(gè)橫波波長內(nèi)包含四個(gè)單元?jiǎng)t耗散明顯減弱。

關(guān)于其他數(shù)值流形式的頻散和耗散特性，以及其對模擬精度的影響仍需進(jìn)一步研究。

附錄A 矩陣Γ11、Γ12、Γ21、Γ22的具體表達(dá)式

式(11)中子矩陣Γ11，Γ12，Γ21，Γ22分別為

(A-1)

(A-2)

(A-3)

(A-4)

其中

(A-5)

(A-6)

(A-7)

(A-8)

(A-9)

(A-10)

(A-11)

(A-12)

(A-13)

(A-14)