透射邊界高頻失穩機理及穩定實現
——P-SV波動

章旭斌，廖振鵬，謝志南

中國地震局工程力學研究所，中國地震局地震工程與工程振動重點實驗室，哈爾濱 150080

0 引言

近場波動數值模擬是地震工程、地球物理、電磁學等學科共同關注的重要研究方向，其涉及內域數值算法和人工邊界條件，其中人工邊界用來截斷外部無限域，并模擬外行波無反射地穿過邊界.迄今常用的內域算法有差分法、有限元法(Finite Element Method, FEM)和譜元法(Liu et al., 2019)等，而人工邊界主要有黏性及黏彈性邊界(Liu and Li, 2005; Zhao et al., 2011)，Higdon 邊界及其輔助變量實現(Higdon, 1987; Hagstrom et al., 2008)，透射邊界(廖振鵬等,2002)，物理阻尼層(Semblat et al., 2011)，加權混合人工邊界(Liu and Sen, 2012)和完美匹配層(Bérenger, 2007)等.由內域和邊界算法的兩者組合將形成諸多近場波動數值模擬方案，其穩定性是獲得可靠模擬結果的前提.

透射邊界亦稱多次透射公式(Multi-Transmitting Formula, MTF)，是一類典型的具有高階精度、低計算量、易于實現等特點的人工邊界條件，其已應用于諸多波動模擬領域，如電磁學(Taflove and Hagness, 2005)、地震工程(陳厚群,2006)、流體力學(Xu et al., 2016)等.然而在波動模擬中MTF存在數值失穩現象.下面結合穩定性分析方法加以概述已有失穩機理的研究結果.

從大量失穩的數值實驗中可以發現，失穩的肇始具有局域性，即由邊界開始向內域擴散，可見失穩僅是邊界與相鄰的運動方程不當耦合所致，而與遠處的運動方程和邊界無關.基于正則模態分析建立的GKS(Gustafsson, Kreiss, Sundstr?m)定理及其群速度解釋(Trefethen, 1983)是分析局域失穩的重要理論.周正華和廖振鵬(2001)基于GKS定理分析了MTF引發的飄移失穩，指出數值解中的零頻零波數諧波從邊界進入內域引發失穩.景立平等(2002)基于群速度解釋分析了SH波動有限元模擬中MTF引發的高頻失穩，指出失穩是由MTF與內域格式組成的半無限模型支持群速度向內而相速度向外的高頻平面諧波所致.Duru等(2015)基于正則模態分析了不同截斷邊界的PML穩定性，指出采用自由邊界截斷的PML存在隨時間增長的模態.然而，正則模態分析亦有其局限性，除局域引發的失穩外，遠處邊界和運動方程亦可能導致全局性失穩.基于人工邊界反射系數可揭示MTF引發的反射放大失穩，其失穩機理為在有限區域內MTF對外行波的反復反射放大(Liao and Liu, 1992; 章旭斌等, 2019).

基于對失穩現象和機理的認識，MTF引發的數值失穩可歸納為局域(Locality)耦合失穩和全局(Globality)反射放大失穩，前者依據失穩形態可區分為高頻失穩和飄移失穩.對于上述失穩問題，研究人員已提出多種穩定措施.針對反射放大失穩，選擇合理的數值模擬參數可避免這類失穩；針對飄移失穩，通過小量修正MTF或內域格式，可有效抑制飄移失穩(周正華和廖振鵬, 2001; 李小軍和楊宇, 2012).針對高頻失穩，曾采用阻尼濾波來消除(廖振鵬等, 2002)，但會影響數值模擬的精度；對于波導模型，可以設置合理的長方形單元長寬比來避免高頻失穩(謝志南和廖振鵬, 2012)，但該措施并不適用全空間模型和實際場地采用的半空間模型.

P-SV彈性矢量波動方程是地震波動模擬常用的計算模型，在地形效應、盆地效應以及全球地震動模擬(Li et al., 2014)都存在應用.MTF在矢量波動模擬中相比標量波動更易引發失穩，而針對標量波動模擬中得到的MTF穩定條件(章旭斌等, 2015)并不能用于更為復雜的矢量波動模擬.另外，GKS定理的群速度解釋雖然是針對一維和標量波動，但其失穩模態能自然地推廣到P-SV波動中的P波和SV波.本文將在以往研究基礎上，將GKS定理及其群速度解釋應用到P-SV波動模擬中MTF的穩定性分析，揭示MTF引發的失穩機理，并采用修改的數值積分方法(Yue and Guddati, 2005)修改有限元剛度來調整內域格式的頻散，消除引發MTF失穩的波動成分，從而穩定實現MTF.同時利用長持時數值實驗加以驗證.

1 數值格式

均勻、各向同性P-SV彈性波動方程

(1)

其中位移u=[uxuy]T，荷載f=[fxfy]T，ρ為密度，μ為彈性模量，λ為拉梅常數，偏導?x=?/?x，?y=?/?y.(1)式對應的弱形式為

(2)

其中w為試函數，Ω為單元e上的積分變量.(2)式引入等參變換

dΩ=Jdξdη,

(3)

其中J為雅克比行列式.位移場和試函數采用線性形函數插值

(4)

其中N為形函數矩陣，ue和we分別為單元節點位移向量和試函數向量.由此可得單元節點格式

(5)

(6)

(7)

(8)

其中算子矩陣H由矩陣L作鏈式法則得到，即

由于一致質量的空間耦聯特征將極大地增加數值模擬的計算量，通常將一致質量作行和集中，從而得到空間解耦的集中質量(廖振鵬, 2002)

(9)

由此(5)式可以寫成

(10)

剛度陣(7)式通常采用Gauss數值積分計算，而本文將采用修改的數值積分方法(Yue and Guddati, 2005)

(11)

邊界區域通常采用規則的長方形單元，采用數值積分方法(11)式分別計算(6)和(7)式可以得到長方形單元節點質量

mi=ρΔxΔy/4

(12)

以及單元剛度陣

(13)

其中α1=1+α2，α2=1-α2，e1=(ε2β2+1)/2β，e2=(ε2+β2)/2β，g=(3-ε2)/2，h=(ε2-1)/2.進而對(10)式采用時間中心差分離散，可得長方形單元的局部節點系格式

(14)

圖1 長方形單元和局部節點系

Lu=O(Δt3+Δx3+Δy3),

(15)

可知無論α取何值，該格式都具有二階精度.

人工邊界采用MTF，在垂直x方向的人工邊界上的節點格式為

(16)

在垂直y方向的人工邊界上的節點格式為

(17)

2 失穩機理及穩定實現

根據對數值失穩現象的觀察，失穩擾動從一個或若干個節點開始出現，然后隨時間增長并向周圍擴散.也就是說失穩的肇始具有局域性，即僅與鄰近的節點運動方程之間的耦合相關，而與遠處的節點運動方程和邊界條件無關.因此，對于內域數值格式的穩定性，可以在無限模型中采用傅里葉模態來分析局域失穩，要求內域格式不支持隨時間步數無限增長的傅里葉模態，即von Neumann穩定條件(Strikwerda, 2004).

邊界引發的失穩是從邊界節點開始，然后逐步向內域擴散.可見，失穩是邊界節點與內域鄰近節點的相互作用所致，其同樣具有局域性.因此，邊界穩定性分析是針對由邊界和內域組成的半無限模型，并采用更為一般的正則模態

(18)

當z=eiωΔt，κ=e-ikxΔx時，正則模態退化為傅里葉模態.數值穩定也就是要求半無限模型不支持隨時間步數無限增長的正則模態(|z|>1,|κ|>1)，即Godunov-Ryabenkii穩定條件(Gustafsson et al., 1972).其中|κ|>1的含義為：在內域穩定條件下，數值格式已不再支持|z|>1，|κ|=1的模態，并考慮到正則模態在無窮遠處應為有限值，因此對右邊界而言，|κ|>1.圖2給出了正則模態示意圖.

圖2 正則模態示意圖

從而可得群速度

(19)

如果Cx<0，意味著z受到擾動向單位圓外移動變為|z|>1時，κ也將向單位圓外移動變為|κ|>1.因此，GKS定理補充條件的物理含義為：半無限模型不支持群速度指向內域的諧波模態.由于群速度是描述波動能量的行進速度，如果半無限模型支持群速度指向內域的諧波模態，波動能量將不斷地從邊界進入內域，從而導致失穩.

GKS定理的群速度解釋物理概念清晰，易于操作.邊界和內域頻散曲線的交點就是兩者共同支持的諧波模態，其交點切線斜率就是諧波的群速度.對于右邊界而言，如果群速度為負，則表示向內行進，邊界將會引發失穩.值得說明的是，GKS定理的群速度解釋雖是針對一維和SH波動的，但其失穩模態能自然地推廣到P-SV波動中的P波和SV波.考慮到P波和SV波的頻散是完全獨立的，因此分別針對P波和SV波，驗證邊界和內域格式是否支持群速度指向內域的失穩諧波模態.

頻散分析所采用的平面諧波形式為

(20)

其時空離散形式為

(21)

其中ω為圓頻率，k為波數，θ為諧波入射方向與x軸正向的夾角，沿x和y軸正向的視波數分別為kx=kcosθ和ky=ksinθ.

將(20)式代入(1)式，并忽略荷載項，可得連續模型的頻散關系

(22)

將(21)式代入(14)式，可得有限元的頻散關系

(23)

其中s1=sin(kxΔx/2)，c1=cos(kxΔx/2)，s2=sin(kyΔy/2)，c2=cos(kyΔy/2).當Δx→0,Δy→0，Δt→0時，(23)式可化簡成(22)式，可知(23)式的兩個式子分別是P波和SV波的頻散.

將(21)式分別代入(16)和(17)式，可得MTF的頻散關系

ωΔt=kxΔx,

(24)

ωΔt=kyΔy.

(25)

由于頻散關系的對稱性及混淆效應(廖振鵬, 2002)，可僅考慮第一象限的頻散，即ωΔt∈[0,π],kxΔx∈[0,π],kyΔy∈[0,π].

圖3為連續模型和傳統有限元的頻散曲線.其中連續模型的頻散改寫為

(26)

首先，從圖3c可以看出傳統有限元P波頻散存在單調遞減的曲線，其與MTF頻散曲線在高頻段(ωΔt>1)存在負切線斜率的交點.由于頻散曲線上一點的切線斜率表示該點對應的平面諧波群速度，對于右邊界而言，群速度為負表明諧波向內域行進.由此可知，MTF和傳統有限元共同支持了群速度指向內域的平面諧波.依據GKS定理的群速度解釋，MTF會引發數值失穩.另外，從圖3a可以看出連續模型P波頻散曲線始終單調遞增，切線斜率始終為正，因此并不存在群速度向內的諧波.這意味著傳統有限元離散引入了在連續模型中并不存在的群速度向內的諧波，這正是引發MTF高頻失穩的諧波成分.

圖3 連續模型和傳統有限元的頻散曲線，及在垂直x方向人工邊界上的MTF頻散曲線，其中參數

基于對失穩機理的認識，若有限元頻散曲線單調遞增，則可消除P-SV波動模擬中MTF引發的高頻失穩.由于有限元頻散由其質量陣和剛度陣決定，(12)和(13)式表明數值積分方法并不會改變長方形單元的集中質量大小，但會改變剛度陣，可見數值積分方法可以改變有限元頻散.因此，在垂直x方向人工邊界上的MTF穩定實現條件，即為有限元P波和SV頻散曲線在x方向上單調遞增

(27)

同理，在垂直y方向人工邊界上的MTF穩定實現的條件為

(28)

這里直接給出(27)式的充要條件

(29)

以及(28)式的充要條件

(30)

上式表明積分點選取與介質波速比和單元長寬比有關.具體推導過程見附錄.

對于正方形單元(β=1)而言，MTF穩定實現條件由(29)和(30)式可簡化為

(31)

圖4 修正有限元和在垂直x方向人工邊界上的MTF頻散曲線，其中參數

內域格式的穩定條件由Von-Neumann穩定性分析方法可以得到.在(29)式條件下，容易得到修正有限元的穩定條件為

(32)

同理，在(30)式條件下，修正有限元穩定條件為

(33)

3 數值實驗

3.1 基巖半空間

考慮基巖均勻半空間模型，介質密度ρ=2.7 g·cm-3，cp=5.5 km·s-1，cs=3.2 km·s-1.計算模型如圖5所示，Xb=2Yb=10 km.在自由地表(2 km,5 km)處作用豎直方向的瑞克子波荷載

圖5 基巖半空間模型

f(t)=a0[1-2(πf0)2(t-t0)2]e-(πf0)2(t-t0)2,

(34)

從圖6可以看出傳統有限元結合MTF引發數值失穩現象，MTF階數越高，失穩增長越劇烈，失穩出現時間越早.圖7顯示修正有限元結合MTF可以穩定實現，各階MTF均未引發高頻失穩，其中三階MTF存在飄移失穩，已通過小量修正MTF加以消除(周正華和廖振鵬, 2001)，小量取為0.005.圖8a為計算模型參數條件下，傳統有限元和MTF的頻散曲線，圖中紅線為一條水平頻散曲線，因此傳統有限元和MTF的負切線斜率交點將落在紅色方框內，從而可以得到引發MTF失穩的諧波頻率范圍，即為圖中紅色方框內的頻帶ωΔt∈[0.91,1.07]，換算為頻率f∈[144.8 Hz,170.3 Hz].圖8b為截取的接收點失穩增長數值解，其表現為節點位移高頻振蕩，通過快速傅里葉變換可以得到其失穩主頻為147.3 Hz，處在理論失穩頻率范圍內.

圖6 傳統有限元結合MTF計算的接收點位移時程

圖7 修正有限元結合MTF計算的接收點位移時程

圖8 (a)計算模型參數條件下，傳統有限元和MTF的頻散曲線，紅線為一條水平頻散曲線，紅色方框為引發MTF失穩的諧波頻率范圍；(b)截取的接收點失穩增長數值解；(c)失穩增長數值解的頻譜

3.2 盆地-基巖半空間

本文得到的MTF穩定性條件是針對規則的長方形單元，然而實際場地存在地形起伏、介質分界面等復雜情形，劃分的有限元網格必然存在非規則單元.由于局域高頻失穩是邊界與鄰近內域節點的相互作用所致，因此只需保證邊界區域為規則的長方形單元且滿足穩定條件，遠處的離散形式并不會影響MTF的穩定性.

在上述基巖半空間上設置深度0.5 km，長度6 km的盆地(如圖9)，介質密度ρ=2.6 g·cm-3，cp=4.5 km·s-1，cs=2.6 km·s-1，其他參數不變，接收點設置在盆地地表距離右邊界1 km處.該模型可用于研究盆地對面波的放大效應(Bowden and Tsai, 2017).本文用于驗證內域存在介質非均勻和網格單元非規則時MTF的穩定性.

圖9 盆地-基巖半空間模型

圖10為修正有限元分別結合MTF和黏彈性邊界計算的接收點位移時程，MTF未出現失穩.與遠置邊界數值解相比，二階和三階MTF的吸收效果較好，一階MTF和黏彈性邊界存在較大的反射波.由于MTF人工波速取為S波波速，提高了對面波的吸收效果，因此顯示一階MTF的精度略高于黏彈性邊界.

圖10 修正有限元結合MTF和黏彈性邊界以及遠置邊界時計算的接收點位移時程

4 結論

P-SV彈性波動有限元模擬中MTF引發的高頻失穩機理可用GKS定理群速度解釋加以闡明，即有限元離散引入了在連續方程中并不存在的群速度指向內域的高頻P波或SV波，并得到了MTF的支持，從而導致數值失穩.可見數值模擬中由于引入人工邊界所導致的失穩并非僅是邊界的問題.本文進而采用修改的數值積分方法修改有限元剛度來調整內域格式的頻散，消除了引發邊界失穩的波動成分，穩定實現了MTF，從而構建了穩定的近場P-SV波動模擬方案.對于實際復雜問題，通過保證邊界區域為規則的長方形單元且滿足穩定條件，數值實驗顯示MTF同樣能夠穩定實現.本文研究思路可進一步推廣到三維彈性波動數值模擬，對于其他類型數值格式中人工邊界的穩定問題亦具有參考價值.

致謝感謝審稿專家對本文提出的寶貴意見.

附錄 證明(27)式的充要條件為(29)式

對(23)式求偏導可得

(A1)

(A2)

其中C0=Δτ2/(4sin(ωΔt))，Q′i=?Qi/?(kxΔx)，i=1,2,3.由(A1)和(A2)式可得(27)式的等價形式

Q′1≥0,

(A3)

(A4)

因此，要證明(27)式的充要條件為(29)式，即證明(A3)和(A4)式的充要條件為(29)式.

下面先證明(A3)和(A4)式的必要條件為(29)式.取值kyΔy=π，此時s2=1，c2=0，Q3=0，分別代入(A3)和(A4)式得到

1-(1-α2)(1+1/β2)≥0,

(A5)

(ε2+1)2[1-(1-α2)(1+1/β2)]2

≥(ε2-1)2[1-(1-α2)(1-1/β2)]2.

(A6)

由于ε>1且α∈[0, 1]，由(A5)和(A6)整理可得

(A7)

(Q′1)2≥(Q′3)2,

(A8)

β≥1,

(A9)

(A7)和(A9)式的組合即(29)式，因此(A3)和(A4)式的必要條件為(29)式.

下面證明(A3)和(A4)式的充分條件亦為(29)式.由(29)式可以推導得到

(1-α2)(1+1/β2)≤(β2+1)/(β2+ε2)<1,

(A10)

由于s1,s2,c1,c2∈[0,1]，因此可得(A3)式成立.

為了證明由(29)式可以得到(A4)式，可以通過證明由(A4)式分解得到的如下兩個不等式成立

(A11)

(A12)

將Qi,Q′i代入(A11)式，經整理可知，要證明(A11)式僅需證明不等式

(A13)

對(A13)式進一步整理可得

(A14)

由(29)式容易得到(A14)式成立.將Qi,Q′i代入(A12)式, 并利用(A13)式，經整理可得

(A15)

由(29)式容易得到(A15)式成立，因此(A3)和(A4)式的充分條件亦為(29)式.綜上可知，(27)式的充要條件為(29)式.同理可以證明(28)式的充要條件為(30)式.