章旭斌，廖振鵬，謝志南

中國地震局工程力學(xué)研究所，中國地震局地震工程與工程振動(dòng)重點(diǎn)實(shí)驗(yàn)室，哈爾濱 150080

0 引言

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

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

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

基于對失穩(wěn)現(xiàn)象和機(jī)理的認(rè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或內(nèi)域格式，可有效抑制飄移失穩(wěn)(周正華和廖振鵬, 2001; 李小軍和楊宇, 2012).針對高頻失穩(wěn)，曾采用阻尼濾波來消除(廖振鵬等, 2002)，但會(huì)影響數(shù)值模擬的精度；對于波導(dǎo)模型，可以設(shè)置合理的長方形單元長寬比來避免高頻失穩(wěn)(謝志南和廖振鵬, 2012)，但該措施并不適用全空間模型和實(shí)際場地采用的半空間模型.

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

1 數(shù)值格式

均勻、各向同性P-SV彈性波動(dòng)方程

(1)

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

(2)

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

dΩ=Jdξdη,

(3)

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

(4)

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

(5)

(6)

(7)

(8)

其中算子矩陣H由矩陣L作鏈?zhǔn)椒▌t得到，即

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

(9)

由此(5)式可以寫成

(10)

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

(11)

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

mi=ρΔxΔy/4

(12)

以及單元?jiǎng)偠汝?/p>

(13)

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

(14)

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

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

(15)

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

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

(16)

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

(17)

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

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

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

(18)

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

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

從而可得群速度

(19)

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

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

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

(20)

其時(shí)空離散形式為

(21)

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

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

(22)

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

(23)

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

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

ωΔt=kxΔx,

(24)

ωΔt=kyΔy.

(25)

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

圖3為連續(xù)模型和傳統(tǒng)有限元的頻散曲線.其中連續(xù)模型的頻散改寫為

(26)

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

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

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

(27)

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

(28)

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

(29)

以及(28)式的充要條件

(30)

上式表明積分點(diǎn)選取與介質(zhì)波速比和單元長寬比有關(guān).具體推導(dǎo)過程見附錄.

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

(31)

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

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

(32)

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

(33)

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

3.1 基巖半空間

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

圖6 傳統(tǒng)有限元結(jié)合MTF計(jì)算的接收點(diǎn)位移時(shí)程

圖7 修正有限元結(jié)合MTF計(jì)算的接收點(diǎn)位移時(shí)程

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

3.2 盆地-基巖半空間

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

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

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

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

圖10 修正有限元結(jié)合MTF和黏彈性邊界以及遠(yuǎn)置邊界時(shí)計(jì)算的接收點(diǎn)位移時(shí)程

4 結(jié)論

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

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

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

對(23)式求偏導(dǎo)可得

(A1)

(A2)

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

Q′1≥0,

(A3)

(A4)

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

下面先證明(A3)和(A4)式的必要條件為(29)式.取值kyΔy=π，此時(shí)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)式可以推導(dǎo)得到

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

(A10)

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

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

(A11)

(A12)

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

(A13)

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

(A14)

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

(A15)

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