章旭斌，廖振鵬，謝志南

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

0 引言

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

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

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

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

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

1 數(shù)值格式

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

(1)

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

(2)

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

dΩ=Jdξdη,

(3)

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

(4)

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

(5)

(6)

(7)

(8)

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

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

(9)

由此(5)式可以寫成

(10)

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

(11)

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

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)式采用時間中心差分離散，可得長方形單元的局部節(jié)點系格式

(14)

圖1 長方形單元和局部節(jié)點系

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

(15)

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

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

(16)

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

(17)

2 失穩(wěn)機理及穩(wěn)定實現(xiàn)

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

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

(18)

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

圖2 正則模態(tài)示意圖

從而可得群速度

(19)

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

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

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

(20)

其時空離散形式為

(21)

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

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

(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為連續(xù)模型和傳統(tǒng)有限元的頻散曲線.其中連續(xù)模型的頻散改寫為

(26)

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

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

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

(27)

同理，在垂直y方向人工邊界上的MTF穩(wěn)定實現(xiàn)的條件為

(28)

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

(29)

以及(28)式的充要條件

(30)

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

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

(31)

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

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

(32)

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

(33)

3 數(shù)值實驗

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

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

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

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

3.2 盆地-基巖半空間

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

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

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

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

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

4 結論

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

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

附錄 證明(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)式，經(jīng)整理可知，要證明(A11)式僅需證明不等式

(A13)

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

(A14)

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

(A15)

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