復雜固-液介質的速度應力方程三角形單元間斷伽遼金地震波模擬算法

曹文忠,張偉

1 中國地震局地球物理研究所,北京 100081 2 南方科技大學地球與空間科學系,廣東深圳 518055

0 引言

地震波數值模擬方法是地球物理學研究的重要領域之一,尤其是計算機技術的飛速發展,使其發揮了越來越重要的作用.近些年來,隨著海底地震儀和光纖地震儀技術的發展,海洋地震學受到了越來越多的關注.構造活動頻繁的區域是地震學研究的重點地區,構造活動會導致海底地形較為復雜.因此,準確高效的大規模復雜固-液介質地震波模擬方法在海洋地震學研究中將會具有非常大的潛在應用價值.

目前已經有多種數值模擬方法應用于地震波傳播模擬問題中,不同的數值模擬方法各自具有優勢,根據需要模擬的實際問題選擇合適的方法.對于復雜固-液介質地震波數值模擬問題,需要顯式施加邊界條件才能確保模擬結果的準確,因此所使用的方法可以靈活的刻畫復雜的固-液界面的形狀,同時還應該易于施加固-液邊界條件.

有限差分法(Graves,1996; Yang et al.,2006; 祝賀君等,2009; Lisitsa and Vishnevskiy,2010; Sun et al.,2019; Zhang and Zhang,2022) 計算效率高,易于并行求解,是目前地震波數值模擬中使用最廣的方法之一.但是其采用結構化網格,存在復雜形狀固液界面時,貼合界面生成網格的難度較大.同時,高階格式下需要較寬的差分模板,在施加邊界條件時,采用降階處理會降低整體格式精度 (Zhang et al.,2014; Sun et al.,2016).Simultaneous-Approximation-Term Summation-By-Parts (SAT-SBP)格式通過單邊差分施加邊界條件,避免了整體格式精度降低的問題(Sun et al.,2020; 楊在林等,2020),但是仍然存在復雜結構網格生成耗時的問題.有限元法(Bao et al.,1998)采用非結構網格,可以靈活處理復雜幾何固液邊界.高階格式時的有限元法求解需要對大型稀疏矩陣求逆,當網格數量較大時,內存需求高且并行求解存在難度.使用四邊形或六面體單元的譜元法(Komatitsch and Vilotte,1998; Komatitsch et al.,2000,2002)解決了有限元法中高階格式時并行求解的問題,但是只能使用四邊形或六面體網格,對于復雜水體結構,仍然存在網格生成困難的問題.有限體積法(Tadi,2004)可以使用非結構網格單元,相鄰網格單元直接使用通量進行連接,這種方式易于施加復雜邊界條件.但是,有限體積法的模擬精度較低(Dumbser et al.,2007a),目前在地震波數值模擬中使用的較少.間斷伽遼金法(K?ser and Dumbser,2006)可以看作是有限元法和有限體積法的結合,單元之間采用通量的方式進行波場的傳遞,雖然增加了部分計算量,但是使得每個單元的求解可以獨立,易于施加邊界條件和并行求解.單元內部使用有限元的基函數,可以在三角形或四面體單元時采用高階格式.間斷伽遼金法非常適用于含有復雜水體的地震波傳播模擬問題.

K?ser和Dumbser(2006)最早將間斷伽遼金法應用于二維地震波數值模擬中,并擴展到三維彈性波模擬問題(Dumbser and K?ser,2006),黏彈性模擬問題(K?ser et al.,2007),各向異性模擬問題(de la Puente et al.,2007)和變時間步長的模擬中(Dumbser et al.,2007b).隨后,間斷伽遼金法又被用于地震波震源動態破裂模擬問題(de la Puente et al.,2009; Pelties et al.,2012),孔彈性模擬問題(de la Puente et al.,2008; Boxberg et al.,2017; Zhan et al.,2018c; Shukla et al.,2019; Zhang et al.,2019)以及固液介質模擬問題(K?ser and Dumbser,2008; Wilcox et al.,2010; Ye et al.,2016; Zhan et al.,2018a; Shukla et al.,2020).

間斷伽遼金法中通量是非常重要的,目前間斷伽遼金法中應用于地震波模擬問題中的數值通量主要有中心通量(Delcourte et al.,2009),Lax-Friedrich通量(廉西猛和張睿璇,2013; 張金波等,2018; 薛昭等,2014; 賀茜君等,2014,2021),Upwind通量(K?ser and Dumbser,2006)和基于Rankine-Hugoniot跳躍條件的通量(RH-condition通量)(Wilcox et al.,2010)四種.Ye等(2016)在固體和液體內部采用中心通量,固-液邊界基于中心通量增加額外的懲罰項實現固-液邊界條件,但是中心通量會導致模擬的精度較低.K?ser和Dumbser(2008)采用Upwind通量實現了復雜固-液介質地震波傳播模擬.Wilcox等(2010)基于二維四邊形和三維六面體網格在一階速度-應變方程中采用RH-condition通量實現了固-液介質間斷伽遼金地震波傳播模擬.Zhan等(2018b)將該通量擴展到各向異性固-液介質的模擬中.

由于有限差分等其他數值模擬方法大部分都采用一階速度-應力方程作為控制方程,開發基于速度-應力方程的間斷伽遼金法更易于同這些數值方法聯用實現混合數值模擬方法.此外,相比于三角形和四面體網格,四邊形和六面體網格單元在復雜結構中應用時網格生成過程復雜.因此,本文針對二維三角形網格單元,發展了基于一階速度-壓力聲波方程和一階速度-應力彈性波方程、采用RH-condition通量的二維復雜固-液介質間斷伽遼金地震波傳播模擬算法,并通過一系列的數值模擬驗證了模擬結果的準確性并分析了間斷伽遼金法空間格式階數與計算效率的關系.本文中使用的時間積分為低存儲的4階Runge-kutta格式(Carpenter 和Kennedy,1994).

1 原理

1.1 聲波及彈性波方程

進行固-液介質地震波傳播模擬時,液體介質用聲波方程來表示,固體介質用彈性波方程來表示.一階波動方程可以表示為

(1)

(1)對于聲波方程:

(2)

其中,P表示壓力,v表示質點振動的速度,ρ表示介質的密度,c表示波在介質中傳播的速度.

(2)對于彈性波方程:

(3)

其中,τ表示應力,v表示質點振動的速度,ρ表示介質的密度,λ和μ是應力-應變關系中介質的拉梅常數.

1.2 節點間斷伽遼金法

間斷伽遼金法根據基函數可以分為模態和節點兩種表示方式,他們之間可以互相轉化.間斷伽遼金法采用的是積分形式的方程來進行求解,其中又可以分為弱形式和強形式兩種表示形式(Hesthaven and Warburton,2007),本文采用強形式方程的節點間斷伽遼金法來進行求解,強形式方程的表達式為

=0.

(4)

其中,D表示一個網格單元,?D表示網格單元的邊界,上標k表示空間中第k個單元,表示拉格朗日插值多項式,其表達式為根據其性質可以建立Gauss-Legendre-Lobatto(GLL)點與基函數之間的關系(Hesthaven and Warburton,2007).u表示待求的波動方程中的變量,表示通量的部分,表示單元邊界的外法線方向.

在進行求解時,如圖1所示,先利用坐標變換將任意形狀的三角形變換為等腰直角的標準參考三角形,其中Ψ表示任意三角形D與標準參考三角形I之間的映射關系.通過該變換每個單元的求解都可以統一起來.

圖1 任意三角形和標準三角形變換示意圖

一個網格單元內的GLL 點分布和階數有關.圖2中N表示階數,從圖中可以看到,隨著階數的增加,一個三角形單元內的GLL點的個數也隨之增加.其對應的基函數為

圖2 不同階數下GLL點在標準三角形中的分布

(5)

有限元法的基函數可以看作在整個計算空間內是連續的,在相鄰單元的邊界隱含了變量的連續性條件,間斷伽遼金法的基函數只作用于單元內部,相鄰單元之間采用通量來進行波場的傳遞.這樣的方式使得每個單元都可以獨立進行求解,因此易于并行處理.此外,對于固-液邊界等非連續性邊界,施加邊界條件也更加靈活.

1.3 RH-condition通量

基于Rankine-Hugoniot跳躍條件推導的數值通量(RH-condition通量),考慮了相鄰單元之間的波阻抗差異,更加符合物理規律.曹文忠(2022)通過一系列的數值模擬發現,在聲波和彈性波介質中,存在強波阻抗差異時,只有RH-condition通量可以始終保持穩定,Upwind通量,Lax-Friedrich通量都會出現穩定性問題.

求解方程是在笛卡爾坐標系下來表示各個分量的,但是在進行通量計算時,使用單元邊界法向分量的波場值.首先對變量進行坐標旋轉得到對應的法向分量,計算出局部的通量后,再進行一次逆變換得到笛卡爾坐標下對應的通量值.定義單元邊界的法線方向n=(nx,nz).

(1)聲波方程坐標旋轉:

(6)

其中,旋轉矩陣T及其逆矩陣為

(7)

(2)彈性波方程坐標旋轉:

(8)

其中,旋轉矩陣T及其逆矩陣為

(9)

對于式(4)中的通量部分,其離散后的形式為

(10)

其中,Ne表示網格單元中邊的個數.

根據聲波方程中的連續性條件:

(11)

可以得到對應的通量(具體推導過程見附錄):

(12)

其中,[[·]]符號表示單元邊界兩側的波場差,例如[[a]]=a+-a-.類似地,彈性波方程的RH-condition通量為

(13)

式(12)和式(13)得到的是沿著單元邊界外法線方向的通量,波動方程是在笛卡爾坐標系進行求解的,將其分別代入式(10),即可得到聲波和彈性波在笛卡爾坐標系下的通量.

1.4 固-液邊界條件

在固-液邊界,如圖3所示,固體介質中界面法向方向的正應力和液體介質中的壓力的相反數存在連續關系,其法向的速度分量也存在連續關系.固體介質中的剪切應力在固-液邊界處為0,其余分量沒有約束.

圖3 固-液邊界條件示意圖

根據圖3中固-液邊界不同分量之間的關系,可以得到固-液邊界條件為

(14)

(1)聲波方程固-液邊界通量處理:

(15)

(2)彈性波方程固-液邊界通量處理:

(16)

實際計算時,在聲波介質中的固-液界面,將式(15)代入式(12),在彈性波介質中的固-液界面,將式(16)代入式(13)即可實現固-液邊界條件.

1.5 吸收邊界條件

本文中使用吸收邊界條件來對避免波場在邊界位置的反射.吸收邊界條件滿足在邊界位置流入的波場為0.

(1)聲波方程吸收邊界通量處理:

(17)

(2)彈性波方程吸收邊界通量處理:

(18)

實際計算時,在計算區域的邊界,將式(17)代入式(12),式(18)代入式(13)既可以實現吸收邊界條件.

2 固-液介質模擬結果

本節中,首先采用水平層狀固-液模型驗證了模擬結果的準確性,并分析了不同網格單元大小與空間格式下模擬效率和模擬精度的關系.然后通過sin型起伏固-液模型進一步驗證在復雜情況下的模擬結果的準確性.最后,對不規則形狀固-液界面的模型進行模擬來說明本文提出的方法可以用于復雜固-液邊界的地震波傳播過程模擬問題.

2.1 水平層狀模型

首先,使用水平層狀固-液模型來驗證模擬結果的準確性.圖4展示了模型的大小,介質參數及震源和檢波器的分布.震源位于固-液界面上方100 m處,檢波器共有兩組,分別位于固-液界面上下各1 m的位置.生成網格的平均邊長約為70 m,最小的邊長51 m,最大的邊長99 m.震源時間函數是中心頻率8 Hz的Ricker子波,震源類型為壓力源.聲波和彈性波介質參數參考了Komatitsch等(2000)文章中使用的介質參數.計算時采用5階格式進行計算,時間步長為0.00085 s,并且使用譜元法作為參考解進行對比,譜元法采用邊長40 m的正方形網格,4階格式,時間步長為0.00085 s.圖5展示了模擬的波形結果,可以看到本文提出的間斷伽遼金法與參考解譜元法差異非常小,這說明施加的固-液邊界條件是準確的.圖6是模擬得到的波場快照,從圖中也可以看到,在固-液界面位置,Vx分量存在不連續的現象,而Vz分量始終保持連續,這也與固-液邊界的物理特征相符.

圖4 水平層狀固-液模型的介質參數,震源和檢波器分布

圖5 水平層狀固-液模型波形對比

在間斷伽遼金法中,單元的邊界存在兩組波場值,尤其是在低階格式時,由于網格單元數量增加,會有大量的節點存在兩組波場值,降低了計算效率.在本文中,單元網格大小是指單元的邊長,由于采用的是三角形網格,在劃分時難以生成大小一致的網格,因此本文中的單元大小指的是平均的網格邊長.表1展示了在每波長網格點基本一致時的對比結果,其中在分析絕對誤差最大值時使用的參考解為120 m單元大小10階格式得到的波形.從中可以看到,1階格式的計算時間為2階格式的2.58倍,且絕對誤差更大.4階格式相比2階格式計算精度更高,且計算時間更少.8階格式計算精度最高,但計算量也有了明顯的增加.在低階格式時,由于存在大量重復的波場值導致計算量增大,導致計算時間增加.隨著階數的增加和網格單元邊長的變大,重復的節點數大量的減少,使得計算量也顯著減小.在計算空間偏導時,需要進行矩陣乘法的運算,隨著階數的繼續增加,矩陣的規模會迅速增加,導致計算量增大,計算時間也隨之增加.更高階的格式理論上使用更少的每波長網格點即可達到低階格式相同的誤差,因此,這里的對比不夠嚴格.不過在特別高階格式時,為了提高計算效率網格單元會更大,但是對于一些復雜結構問題,網格劃分時需要準確刻畫復雜界面形狀,網格單元邊長不能過大,此時選取高階格式會存在一定的問題.

通過上面的對比表明,選取合適的階數才能更好的提高間斷伽遼金法的計算效率.表2分析了在4階格式下,不同網格單元大小時,絕對誤差最大值與計算時間的關系,參考解與表1中的保持一致.可以看到,隨著網格單元的變大,絕對誤差也隨之增加,計算時間與單元數成正比.通過分析可以看到,在4階格式下,每波長5～6個格點時,與參考解的誤差可以控制在1%左右.

表2 4階格式不同單元大小下,節點數、計算時間域絕對誤差最大值對比

2.2 sin型起伏模型

下面進一步對存在固-液界面起伏模型的模擬結果準確性進行驗證,采用的是sin型起伏固-液模型.圖7展示了sin型起伏界面模型的大小,介質參數及震源和檢波器的分布.生成網格單元的最小邊長為25 m,最大的邊長為61 m,平均網格單元的邊長約為40 m.選取空間4階格式進行計算,時間步長為0.00085 s,震源參數與水平層狀模型中的一致.參考解同樣為譜元法模擬得到的波形記錄,其網格單元邊長為50 m,空間4階格式進行計算.圖8展示了模擬得到的波形記錄,與水平層狀固-液模型類似,本文提出的方法得到的結果和參考解的差異非常小,說明了該方法可以用于復雜固-液介質地震波傳播過程模擬問題.

圖7 sin型起伏固-液模型的介質參數,震源和檢波器分布

圖8 sin型起伏固-液模型波形對比

2.3 復雜結構模型

通過前面的對比,說明本文提出的方法可以用于計算復雜固-液模型地震波傳播模擬問題.圖9展示了一個含有不規則復雜固-液邊界的模型.有限差分類方法難以生成對應的結構化網格,傳統的譜元法由于需要使用四邊形網格單元,針對此類問題其網格生成不夠靈活,基于三角形網格的間斷伽遼金法可以靈活的處理此類復雜結構.圖10展示了間斷伽遼金法模擬得到的波場快照,從圖中可以看到,在水平固-液界面位置處,Vx分量存在不連續現象,Vz分量保持連續.在豎直固-液界面位置處,Vx分量保持連續,Vz分量存在不連續現象,這與物理規律相符,也說明了本文提出的間斷伽遼金法可以在復雜結構時準確的施加對應的固-液邊界條件,能夠模擬復雜固-液介質下的地震波傳播.

圖9 復雜固-液模型的介質參數,震源和檢波器分布

圖10 復雜固-液模型波場快照(第一列對應Vx分量,第二列對應Vz分量)

3 結論

本文基于Rankine-Hugoniot跳躍條件得到的RH-condition通量發展了復雜二維固-液介質間斷伽遼金地震波傳播模擬算法,該數值通量考慮了單元兩側的波阻抗的差異,更符合物理規律.使用該算法模擬了水平層狀固-液模型和sin型起伏固-液模型波形記錄,并與譜元法模擬結果進行對比驗證了該方法的準確性.通過復雜模型波場快照說明該方法在復雜模型時依然可以準確的施加固-液邊界條件,能夠對復雜固-液模型地震波傳播過程進行準確模擬.在水平層狀固-液模型中通過一系列的數值實驗表明,低階格式計算效率和精度都較低.在高階格式時,相同精度條件下,網格單元的邊長更大,但是過大的網格邊長在網格劃分時不能準確刻畫出復雜界面的形狀,這會導致模擬結果出現一定的誤差.因此在實際模擬時,應該根據實際問題來選擇階數和網格單元邊長,但建議避免采用低階格式來進行計算.

附錄A RH-condition通量推導過程

Rankine-Hugoniot跳躍條件是計算通量的一種方法(Toro,1997).目前也越來越多的應用于地震波數值模擬問題中(Wilcox et al.,2010; Zhan et al.,2017,2018b; Duru et al.,2022; 曹文忠,2022).曹文忠(2022)給出了該通量在一階速度-壓力聲波方程和一階速度-應力彈性波方程中的詳細介紹及推導過程,本文以聲波方程為例進行簡單的介紹.如附圖1所示,橫軸中0的表示相鄰單元的邊界,其兩側分別對應不同的單元,其介質參數和所需求解的波場值分別用上標“-”和“+”來表示,上標“*”表示加上通量后的波場值.

附圖1 二維聲波RH-condition示意圖

本文在計算通量時,將x方向旋轉到法向方向,將公式(2)中矩陣A進行特征值和特征向量分解可以得到:

(A1)

根據Rankine-Hugoniot跳躍條件可以得到:

(A2)

(A3)

(A4)

根據聲波方程中的連續性條件:

(A5)

根據上式關系,通過求解2元1次方程組可以得到:

(A6)

其中,[[·]]符號表示單元邊界兩側的波場差,例如[[a]]=a+-a-.將式(A6)代入式(10)和(12)即可得到聲波方程的通量.同樣地,可以很容易得到彈性波方程中:

(A7)

將公式(A7)代入公式(10)和(13)即可得到彈性波方程的通量.