透射邊界高頻耦合失穩(wěn)機(jī)理及穩(wěn)定實(shí)現(xiàn)
——SH波動

章旭斌， 廖振鵬， 謝志南

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

?

透射邊界高頻耦合失穩(wěn)機(jī)理及穩(wěn)定實(shí)現(xiàn)
——SH波動

章旭斌， 廖振鵬， 謝志南

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

就大型近場波動的高效數(shù)值模擬而言，穩(wěn)定實(shí)現(xiàn)高階人工邊界是一個尚未圓滿解決的問題.本文針對使用多次透射公式的SH波動集中質(zhì)量有限元模擬，依據(jù)GKS定理的群速度解釋，進(jìn)一步闡明了人工邊界與內(nèi)域離散格式耦合所導(dǎo)致高頻失穩(wěn)的機(jī)理，即兩者支持群速度指向內(nèi)域的外行高頻平面諧波，波動能量自發(fā)地從人工邊界進(jìn)入內(nèi)域，從而導(dǎo)致失穩(wěn),而這類諧波是由集中質(zhì)量有限元離散引入的.本文提出了消除此種耦合失穩(wěn)的一種方法：通過修改有限元剛度陣來改變內(nèi)域離散格式，并保證修改格式的精度不低于原有格式的精度.理論分析和數(shù)值實(shí)驗(yàn)表明此法能穩(wěn)定實(shí)現(xiàn)透射邊界.本文研究結(jié)果具有推廣應(yīng)用前景.

耦合失穩(wěn)； GKS定理； 群速度解釋； 異常頻散； 數(shù)值積分

1 引言

近場波動的數(shù)值模擬需引入人工邊界將無限域截斷成有限的計算區(qū)，并設(shè)置人工邊界條件(Artificial Boundary Condition, ABC)模擬外行波穿過邊界進(jìn)入外部無限域.關(guān)于ABC的研究已有半個世紀(jì)的歷史，研究思路主要建立在單向波概念之上，即ABC應(yīng)保證計算區(qū)內(nèi)的外行波在人工邊界上無反射，由此建立了眾多ABC.建立ABC的方法主要有兩類：其一由滿足場方程的外行波假定導(dǎo)出，如粘性ABC(Lysmer and Kulemeyer, 1969)，黏彈性ABC(Deeks and Randolph, 1994; 劉晶波等, 2005；杜修力等, 2006)，高階彈簧-阻尼-質(zhì)量ABC(Du and Zhao, 2010)，Clayton-Engquist ABC(Clayton and Engquist, 1977)，Higdon ABC(Higdon, 1986; 1990)，輔助變量實(shí)現(xiàn)的Higdon ABC (Givoli and Neta, 2003; Hagstrom and Warburton, 2004; Baffet et al., 2012).其二在人工邊界外設(shè)置吸收層，通過在吸收層內(nèi)設(shè)置場方程，衰減掉進(jìn)入層內(nèi)的外行波，以消除人工邊界的反射.代表性工作為完美匹配層(Bérenger, 1994; Chew and Weedon, 1994; Xie et al., 2014).

以上ABC的建立均遵循傳統(tǒng)的研究思路：ABC具有與場方程相匹配的連續(xù)形式，然后離散化實(shí)現(xiàn)數(shù)值模擬.Liao等(1984a，1984b)提出了可描述體波、面波和導(dǎo)波的單向波形式表達(dá)式，以離散形式直接模擬外行波穿過邊界，由此建立的ABC稱為多次透射公式(Multi-Transmitting Formula, MTF).MTF與前述ABC的區(qū)別如下：第一，MTF與內(nèi)域場方程無關(guān)，可直接用于各種場方程，而PML則需針對特定場方程建立相應(yīng)的ABC；第二，MTF由單向波形式表達(dá)式導(dǎo)出，可適用于人工邊界上的各節(jié)點(diǎn)，如角點(diǎn)和分層界面點(diǎn)，而應(yīng)用輔助變量實(shí)現(xiàn)的Higdon ABC則較難處理；第三，MTF是一維表達(dá)式，其數(shù)值實(shí)現(xiàn)方式靈活、易于實(shí)現(xiàn)且計算量低；第四，MTF的精度階含義明確且可控，易與內(nèi)域離散格式的精度階一致.但是，在數(shù)值模擬中MTF存在數(shù)值穩(wěn)定問題.文獻(xiàn)(Liao et al., 1984a; 廖振鵬等, 2002)中的算例清晰地顯示了MTF的潛在應(yīng)用前景，但更仔細(xì)的數(shù)值實(shí)驗(yàn)結(jié)果表明延長計算時間仍將發(fā)生失穩(wěn).已有的MTF失穩(wěn)機(jī)理研究結(jié)果可歸納為兩類：計算區(qū)內(nèi)外行波的多次反射放大(Liao and Liu, 1992)和局部區(qū)域的耦合失穩(wěn)(李小軍和廖振鵬, 1996; 周正華和廖振鵬, 2001; 景立平等, 2002).前者是人工邊界對高頻行波的反射放大，及在有限計算區(qū)內(nèi)反復(fù)放大所致.后者則為ABC與相鄰內(nèi)域離散格式的耦合所致，可由GKS定理(Gustafsson et al.,1972)的群速度解釋(Trefethen, 1983)加以闡明，即模型(人工邊界與內(nèi)域離散格式)支持群速度指向內(nèi)域的外行平面諧波，從而波動能量將自發(fā)地從邊界進(jìn)入內(nèi)域，導(dǎo)致失穩(wěn).由于前者可通過保證反射系數(shù)小于1避免(Xie and Liao, 2008)，下面概述消除局部耦合失穩(wěn)的研究結(jié)果.

就MTF與相鄰內(nèi)域離散格式的耦合所造成的局部失穩(wěn)而言，可通過消除失穩(wěn)耦合達(dá)到(廖振鵬和謝志南, 2011)，即改變內(nèi)域離散格式，或MTF，或兩者.局部耦合失穩(wěn)的形態(tài)可區(qū)分為高頻失穩(wěn)(景立平等, 2002)和零頻飄移失穩(wěn)(李小軍和廖振鵬, 1996; 周正華和廖振鵬, 2001)，飄移失穩(wěn)可采用加小量修改MTF(周正華和廖振鵬, 2001)或內(nèi)域離散格式(李小軍和楊宇, 2012)來消除.高頻失穩(wěn)常出現(xiàn)在對近場波動數(shù)值模擬無意義的高頻段，曾采用阻尼濾波來消除(廖振鵬等, 2002)，但阻尼大小的選取無確定準(zhǔn)則，過小仍可能失穩(wěn)，過大則影響計算精度.謝志南和廖振鵬(2012)針對波導(dǎo)的有限元數(shù)值模擬，通過調(diào)整長方形單元的長寬比來改變內(nèi)域離散格式的頻散關(guān)系，使之不再支持群速度指向內(nèi)域的外行平面諧波，穩(wěn)定實(shí)現(xiàn)了MTF.但對于半空間或全空間模型，此法失效.

由上述討論可知，高頻耦合失穩(wěn)目前仍無有效的消除方法，這是本文的論題.本文的研究思路是通過改變內(nèi)域離散格式來消除失穩(wěn)耦合，達(dá)到穩(wěn)定實(shí)現(xiàn)MTF的目的.關(guān)于改變內(nèi)域離散格式的工作應(yīng)提及Yue 和Guddati(2005)針對標(biāo)量波動的數(shù)值模擬提出了一種修改的數(shù)值積分方法，即對質(zhì)量和剛度矩陣積分表達(dá)式進(jìn)行數(shù)值積分時，通過對頻散誤差的控制來選取積分點(diǎn)，使正方形單元頻散誤差由原來的二階精度提高到四階精度.Idesman和Pham(2014)則將這一方法推廣到彈性波情形.我們將進(jìn)一步研究這一改變內(nèi)域離散格式的方法，達(dá)到消除高頻耦合失穩(wěn)，同時保證內(nèi)域離散格式的精度不低于原有格式的精度.

本文第二節(jié)給出SH波動數(shù)值模擬的有限元離散格式.第三節(jié)分析離散格式的頻散關(guān)系，進(jìn)而闡明MTF高頻耦合失穩(wěn)的機(jī)理，并分析剛度陣對有限元頻散的影響.第四節(jié)基于修改的數(shù)值積分方法來修改剛度，提出一個穩(wěn)定實(shí)現(xiàn)MTF的集中質(zhì)量有限元格式.第五節(jié)利用數(shù)值實(shí)驗(yàn)驗(yàn)證本文措施的有效性.最后，對本文工作做了小結(jié).

2 有限元離散格式

各向同性、均勻、線彈性SH波動方程為

(1)

采用線性有限元(廖振鵬, 2002)離散上述連續(xù)模型可得

Mü+KU=F,

(2)

其中U為節(jié)點(diǎn)位移矢量，ü為其二階時間導(dǎo)數(shù).系統(tǒng)矩陣M，K分別由單元矩陣Me，Ke組裝而成，系統(tǒng)荷載向量F由單元荷載矩陣fe組裝而成.

(3)

(4)

(5)

其中N是線性形函數(shù)矩陣，Ωe代表單元所在空間區(qū)域.

用中心差分對時間離散，可得全離散格式

M(Up+1-2Up+Up-1)+Δt2KUp=Δt2Fp,

(6)

其中Δt為時間步距，Up和Fp為pΔt時刻的位移和荷載向量.

全離散格式(6)可以寫成與節(jié)點(diǎn)(m,n)相關(guān)的局部節(jié)點(diǎn)系格式，即四個相鄰單元組成的9節(jié)點(diǎn)系(如圖1)，即

圖1 長方形單元節(jié)點(diǎn)編號和二維局部節(jié)點(diǎn)系Fig.1 Node number and 2D local grid system for rectangle element

(7)

采用長方形單元(如圖1)，其邊長為Δx和Δy，則由(3)和(4)可得單元質(zhì)量和剛度矩陣

(8a)

(8b)

其中，q=Δy/Δx.一致質(zhì)量矩陣元素

對一致質(zhì)量矩陣作行和集中得到集中質(zhì)量矩陣，其元素

剛度矩陣元素

ABC采用MTF(Liaoetal., 1984a)，

(9)

3 高頻耦合失穩(wěn)機(jī)理

GKS定理的群速度解釋表明，若內(nèi)域離散格式和人工邊界支持群速度指向內(nèi)域的外行平面諧波，則將導(dǎo)致失穩(wěn).內(nèi)域離散格式支持的平面諧波由頻散分析可知，而對于人工邊界支持的平面諧波，可將人工邊界進(jìn)行無限延拓當(dāng)成內(nèi)域離散格式進(jìn)行頻散分析，本文將無限延拓的MTF稱之為MTF格式.因此，若內(nèi)域與人工邊界格式頻散曲線的交點(diǎn)對應(yīng)群速度向內(nèi)的外行平面諧波，則導(dǎo)致耦合失穩(wěn).

首先對連續(xù)模型、一致質(zhì)量有限元、集中質(zhì)量有限元和MTF格式作頻散分析，利用GKS定理的群速度解釋來揭示MTF高頻耦合失穩(wěn)的機(jī)理.

考慮平面諧波

(10a)

其時空離散形式為

(10b)

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

作頻散分析時，可忽略荷載項(xiàng)，將平面諧波(10a)代入(1)式，可得連續(xù)模型的頻散關(guān)系

).

(11)

將平面諧波(10b)代入局部節(jié)點(diǎn)系格式(7)，可得一致質(zhì)量有限元的頻散關(guān)系

(12)

(13)

將平面諧波(10b)代入(9)，當(dāng)人工邊界垂直x方向時，可得MTF格式的頻散關(guān)系

ωΔt=kxΔx,

(14a)

同理，當(dāng)人工邊界垂直y方向時，可得MTF格式的頻散關(guān)系

ωΔt=kyΔy,

(14b)

由于頻散關(guān)系的對稱性，可只考慮第一象限，即ωΔt∈[0,π],kxΔx∈[0,π],kyΔy∈[0,π].

圖2給出q=1, Δτ=0.4時，連續(xù)模型和離散格式在第一象限的頻散，其中連續(xù)模型的頻散關(guān)系式(11)改寫為

(ωΔt)2=Δτ2((kxΔx)2+(kyΔy)2/q2).

(15)

圖3給出q=1, Δτ=0.4時，在給定不同取值kyΔy下，離散格式在第一象限的頻散曲線.

主要說明兩點(diǎn)：第一，從圖3(a、b)可以看出一致質(zhì)量有限元與MTF格式頻散曲線的交點(diǎn)斜率均為正，而集中質(zhì)量有限元與MTF格式頻散曲線的交點(diǎn)斜率在高頻段為負(fù)(高頻指ωΔt接近內(nèi)域離散格式的截止頻率).由此可知，一致質(zhì)量有限元結(jié)合MTF不會導(dǎo)致高頻耦合失穩(wěn)，而集中質(zhì)量有限元結(jié)合MTF將導(dǎo)致高頻耦合失穩(wěn).

第二，由連續(xù)模型的頻散式(15)，即圖2a可知，其x和y方向頻散曲線的割線和切線斜率均為正，表明第一象限的波動相速度和群速度均為正.一致質(zhì)量有限元的頻散圖2b或圖3a也表明其具有相同性質(zhì)，而從集中質(zhì)量有限元頻散圖2c或圖3b可知，

圖2 q=1, Δτ=0.4時，連續(xù)模型和離散格式在第一象限的頻散(a) 連續(xù)模型; (b) 一致質(zhì)量有限元; (c) 集中質(zhì)量有限元; (d) 人工邊界垂直x方向時的MTF格式.Fig.2 The dispersion of continuous model and discrete schemes at the first quadrant when q=1, Δτ=0.4.(a) Continuous model; (b) Consistent mass FEM; (c) Lumped mass FEM; (d) MTF scheme when ABC is perpendicular to x-axis.

圖3 q=1, Δτ=0.4時，在給定kyΔy不同取值下，離散格式在第一象限的頻散曲線(a) 一致質(zhì)量有限元(實(shí)線)與人工邊界垂直x方向時的MTF格式(虛線)的頻散曲線; (b) 集中質(zhì)量有限元(實(shí)線)與人工邊界垂直x方向時的MTF格式(虛線)的頻散曲線.Fig.3 The dispersion curves of discrete schemes at the first quadrant under the given different value of kyΔy when q=1, Δτ=0.4(a) The dispersion curves of consistent mass FEM (solid lines) and MTF scheme (dash line) when ABC is perpendicular to x-axis; (b) The dispersion curves of lumped mass FEM (solid lines) and MTF scheme (dash line) when ABC is perpendicular to x-axis.

其x和y方向頻散曲線在高頻段的割線斜率為正，而切線斜率為負(fù)，表明這類高頻波動相速度為正，而群速度為負(fù).因此，集中質(zhì)量有限元在高頻段出現(xiàn)了有別于連續(xù)模型和一致質(zhì)量有限元的異常頻散.這異常頻散對應(yīng)的平面諧波就是集中質(zhì)量有限元允許的群速度向內(nèi)的外行高頻平面諧波.可見，消除集中質(zhì)量有限元的異常頻散可穩(wěn)定實(shí)現(xiàn)MTF.

雖然一致質(zhì)量有限元與MTF結(jié)合不會導(dǎo)致高頻耦合失穩(wěn)，但一致質(zhì)量有限元空間耦聯(lián)，計算量大.因此，對于大型波動數(shù)值模擬有必要采用集中質(zhì)量有限元.

本文不探討零頻飄移失穩(wěn)，但同樣可以采用GKS定理的群速度解釋，即內(nèi)域和邊界支持群速度向內(nèi)的外行零頻平面諧波，可參看文獻(xiàn)(周正華和廖振鵬, 2001).本文的數(shù)值實(shí)驗(yàn)采用該文獻(xiàn)提出的小量處理MTF來消除零頻飄移失穩(wěn).

以上簡化分析表明，有限元質(zhì)量集中造成其異常頻散，進(jìn)而導(dǎo)致MTF高頻耦合失穩(wěn).由于有限元的頻散由其質(zhì)量和剛度陣決定，而集中質(zhì)量可以極大地提高計算效率，因此可修改剛度陣以消除有限元與MTF共同支持的群速度向內(nèi)的外行平面諧波.下面在一般情形下分析剛度對有限元頻散的影響，進(jìn)一步論證MTF的失穩(wěn)機(jī)理.

將平面諧波(10b)代入局部節(jié)點(diǎn)系格式(7)，把剛度當(dāng)成未知量，僅利用關(guān)系

K0+Kx+Kxy+Ky=0,

(16)

(17)

式(17)中的各項(xiàng)對kxΔx求偏導(dǎo)，可得

(18)

考慮到MTF格式的頻散關(guān)系(14a)和(14b)，其為三維斜面，分布在頻率[0,π]和波數(shù)[0,π]內(nèi)，必然與有限元的x和y方向上任意一條頻散曲線相交.因此，為保證集中質(zhì)量有限元和MTF格式的頻散曲線交點(diǎn)斜率為正，就要要求有限元的頻散曲線單調(diào)遞增.

由式(18)可得有限元x方向頻散曲線單調(diào)遞增的條件：

Kx≤Kxy,

(19a)

同理可得有限元y方向頻散曲線單調(diào)遞增的條件：

Ky≤Kxy,

(19b)

式(19)是保證集中質(zhì)量有限元不存在異常頻散，與MTF結(jié)合不會導(dǎo)致高頻耦合失穩(wěn)的條件.

由此可知，對于正方形單元，Kx=Ky=-1/6，Kxy=-1/3，不滿足式(19).因此，垂直x或y方向的人工邊界MTF均會導(dǎo)致高頻耦合失穩(wěn).

4 修正的集中質(zhì)量有限元

由上節(jié)可知，修改剛度可消除集中質(zhì)量有限元的異常頻散，能穩(wěn)定實(shí)現(xiàn)MTF.因此，對單元剛度的積分計算式(4)，不采用精確積分或是常用的Gauss數(shù)值積分，而采用修改的數(shù)值積分方法(YueandGuddati, 2005)，即

(20)

其中積分點(diǎn)αk∈[0,1]，det(J)為Jacobi行列式，對于長方形單元，det(J)=Δx·Δy/4.

利用數(shù)值積分(20)式，可得修正剛度

將修正剛度代入式(19a)，可得穩(wěn)定實(shí)現(xiàn)垂直x方向的人工邊界MTF的條件

(21a)

同理，將修正剛度代入式(19b)，可得穩(wěn)定實(shí)現(xiàn)垂直y方向的人工邊界MTF的條件

(21b)

令qr為長方形單元長邊與短邊的比值，則穩(wěn)定實(shí)現(xiàn)垂直x和y方向的人工邊界MTF的條件為

(21c)

式(21)表明選取合適的數(shù)值積分點(diǎn)αk可以穩(wěn)定實(shí)現(xiàn)MTF.對于積分點(diǎn)αk的取值，將依據(jù)格式精度和穩(wěn)定性來確定.

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

(22)

可知無論αk取何值，該格式具有二階精度，但也無法找到合適的αk來提高格式的精度.對于波動數(shù)值模擬，格式的頻散誤差與精度階一樣能衡量格式精度.本文同文獻(xiàn)(YueandGuddati, 2005;IdesmanandPham, 2014)一樣采用頻散誤差進(jìn)一步衡量格式精度，以其截斷系數(shù)確定較優(yōu)的αk取值以保證修正有限元的精度不低于傳統(tǒng)有限元的精度.

對修正有限元的頻散式(17)引入余量

(23)

對R中的各項(xiàng)作Taylor展開(ωΔt,kxΔx,kyΔy? 1)，并代入連續(xù)模型的頻散式(15)，則余量R可衡量有限元格式的頻散誤差.由此可得

R=a4(kΔx)4+O((kΔx)6),

(24)

(25a)

對于傳統(tǒng)有限元，有

+2q2sin2θcos2θ}.

(25b)

當(dāng)q≥1時，并由下文穩(wěn)定性條件(30)式可知Δτ≤1，故式(25a)和(25b)的共同項(xiàng)1-Δτ2+(q2-1)sin4θ≥0.若對于兩式中的sin2θcos2θ項(xiàng)系數(shù)有

(26a)

(26b)

(26c)

當(dāng)q<1時，先利用(sin2θ+cos2θ)2=1改寫式(25a)，即

(27a)

對于傳統(tǒng)有限元，

+2sin2θcos2θ}.

(27b)

由下文穩(wěn)定性條件(30)式可知Δτ≤q，故式(27a)和(27b)中的q2-Δτ2+(1-q2)cos4θ≥0.同理可得

(28a)

即

(28b)

(28c)

可以發(fā)現(xiàn)(26c)和(28c)是一樣的，故無論q取何值，滿足式(28c)的αk能保證修正有限元的精度不低于傳統(tǒng)有限元的精度.

下面給出幾個滿足式(28c)和MTF穩(wěn)定實(shí)現(xiàn)條件式(21c)的αk取值，可取

(29a)

此時(25a)和(27a)中的sin2θcos2θ項(xiàng)系數(shù)為0.

或取

(29b)

這是滿足式(21c)的最小取值.在有限元穩(wěn)定條件方面，由下文的穩(wěn)定性條件(30)容易驗(yàn)證(29b)取值優(yōu)于(29a)取值.

或取

αk=1,

(29c)

此時數(shù)值積分(20)式即為兩點(diǎn)Gauss-Lobatto積分，修正有限元格式等同中心差分格式.

內(nèi)域離散格式的穩(wěn)定性可由Von-Neumann穩(wěn)定性分析得到，即要求波數(shù)kx和ky為實(shí)數(shù)的平面諧波(10b)式不能隨時間無限增長.由格式的頻散關(guān)系(17)式可知，kx和ky為實(shí)數(shù)時，ω必為實(shí)數(shù)，故sin2(ωΔt)≤1，進(jìn)而由(17)式可得修正有限元的穩(wěn)定性條件

(30)

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

Liao等(1984a)中的算例由于計算時間不長，未出現(xiàn)MTF與傳統(tǒng)有限元耦合所致的高頻失穩(wěn).為了揭示MTF的高頻數(shù)值失穩(wěn)現(xiàn)象，以及驗(yàn)證MTF結(jié)合修正有限元可以穩(wěn)定實(shí)現(xiàn)，本文數(shù)值實(shí)驗(yàn)采用Liao等(1984a)的均勻彈性半空間模型.

模型的幾何關(guān)系如圖4所示，剪切波速c=1 m·s-1.波源為在自由表面處沿z方向作用的暫態(tài)荷載：

F(t,x)=T(t)S(x)，

(31)

其中輸入時間函數(shù)為三角形脈沖：

(32)

荷載空間分布函數(shù)S(x)由下式給出：

(33)

其余三邊設(shè)置人工邊界MTF，觀測點(diǎn)設(shè)置在自由地表(1,0)處.

本文的人工邊界穩(wěn)定性分析是針對透射公式(9)，首先用數(shù)值實(shí)驗(yàn)驗(yàn)證其穩(wěn)定性.

圖5中的(a)和(b)分別為采用傳統(tǒng)和修正有限

圖4 彈性半空間SH波動模型Fig.4 Elastic half-space model of SH wave motion

圖5 (a)和(b)分別為采用傳統(tǒng)有限元和修正有限元結(jié)合MTF(式(9)，N=1,2,3)計算的觀測點(diǎn)(1,0)位移時程，離散方案為Δx=0.05 m，Δy=0.025 m，Δt=0.02 sFig.5 (a) and (b) are the displacement time history of the observation point (1,0) computed by classic FEM and modified FEM with MTF(formula (9), N=1,2,3) separately. The discretization mesh is Δx=0.05 m, Δy=0.025 m, Δt=0.02 s

圖6 (a)和(b)分別為采用傳統(tǒng)有限元和修正有限元結(jié)合MTF(空間內(nèi)插公式，人工波速ca=c，N=1,2,3)計算的觀測點(diǎn)(1,0)位移時程，離散方案為Δx=Δy=0.05 m，Δt=0.025 sFig.6 (a) and (b) are the displacement time history of the observation point (1,0) computed by classic FEM and modified FEM with MTF(spatial interpolation formula with artificial wave velocity ca=c, N=1,2,3) separately. The discretization mesh is Δx=Δy=0.05 m, Δt=0.025 s

對于MTF的計算點(diǎn)和離散網(wǎng)格節(jié)點(diǎn)不重合時，MTF的計算點(diǎn)采用空間內(nèi)插(Liao et al., 1984a)，因?yàn)閮?nèi)域格式已消除異常頻散，故不會導(dǎo)致高頻耦合失穩(wěn).下面也通過數(shù)值實(shí)驗(yàn)驗(yàn)證.

6 結(jié)論

對于大型近場波動數(shù)值模擬，建立高精度、低計算量、穩(wěn)定的人工邊界是一個尚未圓滿解決的問題.MTF具有高精度、低計算量，但其存在穩(wěn)定問題.本文針對SH波動集中質(zhì)量有限元模擬，闡明了MTF高頻耦合失穩(wěn)機(jī)理，穩(wěn)定實(shí)現(xiàn)了MTF.本文的研究結(jié)果對于地球物理中的波動正反演問題具有應(yīng)用前景.依據(jù)作者已獲得的研究結(jié)果，對于其他內(nèi)域離散格式(如差分、譜元等)，只要其不出現(xiàn)異常頻散，均不會導(dǎo)致MTF的高頻耦合失穩(wěn).本文的研究工作是針對SH波動進(jìn)行的，但其研究思路適用于更廣類型波動數(shù)值模擬中人工邊界的穩(wěn)定問題.關(guān)于在P-SV和三維彈性波數(shù)值模擬中MTF的穩(wěn)定實(shí)現(xiàn)，將在后續(xù)報告中提出.

Baffet D, Bielak J, Givoli D, et al. 2012. Long-time stable high-order absorbing boundary conditions for elastodynamics.ComputerMethodsinAppliedMechanicsandEngineering, 241-244: 20-37.

Bérenger J P. 1994. A perfectly matched layer for the absorption of electromagnetic waves.JournalofComputationalPhysics, 114(2): 185-200.

Chew W C, Weedon W H. 1994. A 3-D perfectly matched medium from modified Maxwell's equations with stretched coordinates.MicrowaveOpt.Technol.Lett., 7(13): 599-604.

Clayton R, Engquist B. 1977. Absorbing boundary conditions for acoustic and elastic wave equations.Bull.Seism.Soc.Am., 67(6): 1529-1540.

Deeks A J, Randolph M F. 1994. Axisymmetric time-domain transmitting boundaries.JournalofEngineeringMechanics, 120(1): 25-42.

Du X L, Zhao M, Wang J T. 2006. A stress artificial boundary in FEA for near-field wave problem.ChineseJournalofTheoreticalandAppliedMechanics(in Chinese), 38(1): 49-56.

Du X L, Zhao M. 2010. A local time-domain transmitting boundary for simulating cylindrical elastic wave propagation in infinite media.SoilDynamicsandEarthquakeEngineering, 30(10): 937-946.

Givoli D, Neta B. 2003. High-order non-reflecting boundary scheme for time dependent waves.J.Comput.Phys., 186(1): 24-46.

Gustafsson B, Kreiss H O, Sundstr?m A. 1972. Stability theory of difference approximations for mixed initial boundary value problems II.MathematicsofComputation, 26(119): 649-686.

Hagstrom T, Warburton T. 2004. A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to first order systems.WaveMotion, 39(4): 327-338.

Higdon R L. 1986. Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation.MathematicsofComputation, 47(176): 437-459.Higdon R L. 1990. Radiation boundary conditions for elastic wave propagation.SIAMJournalonNumericalAnalysis, 27(4): 831-869.

Idesman A, Pham D. 2014. Finite element modeling of linear elastodynamics problems with explicit time-integration methods and linear elements with the reduced dispersion error.ComputerMethodsinAppliedMechanicsandEngineering, 271: 86-108.

Jing L P, Liao Z P, Zou J X. 2002. A high-frequency instability mechanism in numerical realization of multi-transmitting formula.EarthquakeEngineeringandEngineeringVibration(in Chinese), 22(1): 7-13.

Liao Z P, Wong H L, Yang B P, et al. 1984a. A transmitting boundary for transient wave analyses.ScientiaSinica(SeriesA), 27(10): 1063-1076.

Liao Z P, Wong H L. 1984b. A transmitting boundary for the numerical simulation of elastic wave propagation.InternationalJournalofSoilDynamicsandEarthquakeEngineering, 3(4): 174-183.Liao Z P, Liu J B. 1992. Numerical instabilities of a local transmitting boundary.EarthquakeEngineeringandStructuralDynamics, 21(1): 65-77.

Liao Z P, Xie Z N. 2011. Stability of numerical simulation of wave motion.JournalofHarbinEngineeringUniversity(in Chinese), 32(9): 1254-1261.Liao Z P. 2002. Introduction to Wave Motion Theories in Engineering (Second edition) (in Chinese). Beijing: Science Press.Liao Z P, Zhou Z H, Zhang Y H. 2002. Stable implementation of transmitting boundary in numerical simulation of wave motion.ChineseJournalofGeophysics(in Chinese), 45(4): 554-568, doi: 10.3321/j.issn:0001-5733.2002.04.011.

Li X J, Liao Z P. 1996. The drift instability of local transmitting boundary in time domain.ActaMechanicaSinica(in Chinese), 28(5): 627-632.

Li X J, Yang Y. 2012. Measures for stability control of transmitting boundary.ChineseJournalofGeotechnicalEngineering(in Chinese), 34(4): 641-645.

Liu J B, Wang Z Y, Du X L, et al. 2005. Three-dimensional visco-elastic artificial boundaries in time domain for wave motion problems.EngineeringMechanics(in Chinese), 22(6): 46-51.

Lysmer J, Kulemeyer R L. 1969. Finite dynamic model for infinite media.JournalofEngineeringMechanicsDivision, 95(4): 859-878.

Trefethen L N. 1983. Group velocity interpretation of the stability theory of Gustafsson Kreiss and Sundstr?m.JournalofComputationalPhysics, 49(2): 199-217.

Xie Z N, Komatitsch D, Martin R, et al. 2014. Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML.Geophys.J.Int., 198(3): 1714-1747.Xie Z N, Liao Z P. 2008. A note for the mechanism of high-frequency oscillation instability resulted from absorbing boundary conditions.ActaSeismologicaSinica, 21(3): 306-310.Xie Z N, Liao Z P. 2012. Mechanism of high frequency instability caused by transmitting boundary and method of its elimination—SH wave.ChineseJournalofTheoreticalandAppliedMechanics(in Chinese), 44(4): 745-752.Yue B, Guddati M N. 2005. Dispersion-reducing finite elements for transient acoustics.J.Acoust.Soc.Am., 118(4): 2132-2141.

Zhou Z H, Liao Z P. 2001. A measure for eliminating drift instability of the multi-transmitting formula.ActaMechanicaSinica(in Chinese), 33(4): 550-554.

附中文參考文獻(xiàn)

杜修力, 趙密, 王進(jìn)廷. 2006. 近場波動模擬的人工應(yīng)力邊界條件. 力學(xué)學(xué)報, 38(1): 49-56.

景立平, 廖振鵬, 鄒經(jīng)湘. 2002. 多次透射公式的一種高頻失穩(wěn)機(jī)制. 地震工程與工程振動, 22(1): 7-13.

廖振鵬. 2002. 工程波動理論導(dǎo)論(第二版). 北京: 科學(xué)出版社.

廖振鵬, 周正華, 張艷紅. 2002. 波動數(shù)值模擬中透射邊界的穩(wěn)定實(shí)現(xiàn). 地球物理學(xué)報, 45(4): 554-568, doi: 10.3321/j.issn:0001-5733.2002.04.011.

廖振鵬, 謝志南. 2011. 波動數(shù)值模擬的穩(wěn)定性. 哈爾濱工程大學(xué)學(xué)報, 32(9): 1254-1261.

李小軍, 廖振鵬. 1996. 時域局部透射邊界的計算飄移失穩(wěn). 力學(xué)學(xué)報, 28(5): 627-632.

李小軍, 楊宇. 2012. 透射邊界穩(wěn)定性控制措施探討. 巖土工程學(xué)報, 34(4): 641-645.

劉晶波, 王振宇, 杜修力等. 2005. 波動問題中的三維時域粘彈性人工邊界. 工程力學(xué), 22(6): 46-51.

謝志南, 廖振鵬. 2012. 透射邊界高頻失穩(wěn)機(jī)理及其消除方法-SH波動. 力學(xué)學(xué)報, 44(4): 745-752.

周正華, 廖振鵬. 2001. 消除多次透射公式飄移失穩(wěn)的措施. 力學(xué)學(xué)報, 33(4): 550-554.

(本文編輯 胡素芳)

Mechanism of high frequency coupling instability and stable implementation for transmitting boundary—SH wave motion

ZHANG Xu-Bin, LIAO Zhen-Peng, XIE Zhi-Nan

InstituteofEngineeringMechanics,ChinaEarthquakeAdministration,KeyLaboratoryofEarthquakeEngineeringandEngineeringVibrationofChinaEarthquakeAdministration，Harbin150080,China

Stable implementation of high order artificial boundary conditions has not been solved satisfactorily for wave simulation in infinite domain. The instability problem is studied for the multi-transmitting formula, one of the promising boundary schemes. Focusing on the lumped-mass finite element simulation of SH wave motion, mechanism of high frequency instability caused by coupling the formula and the finite element scheme is further clarified based on the group velocity interpretation of GKS theorem. The mechanism is that both the formula and the finite element scheme support out-going high frequency waves whose group velocity points to interior domain, thus wave energy is allowed to enter spontaneously from the boundary into interior domain. Such type of waves do not exist in the corresponding continuous model, it is solely induced by the finite element. A method for destroying the coupling is presented by changing the interior scheme via modifying the stiffness matrix of FEM using proper integration rule, meanwhile the accuracy of the modified scheme is no less than that of the original one. Theoretical analysis and numerical experiments have shown this method can ensure MTF to be implemented stably. The research approach developed in this study can be applied to cope with the boundary instable problems for other types of wave simulation, and the research results have application prospect for forward and inverse problems in geophysics.Keywords Coupling instability; GKS theorem; Group velocity interpretation; Abnormal dispersion; Modified integration rules

10.6038/cjg20151017.

Zhang X B, Liao Z P, Xie Z N. 2015. Mechanism of high frequency coupling instability and stable implementation for transmitting boundary—SH wave motion.ChineseJ.Geophys. (in Chinese),58(10):3639-3648,doi:10.6038/cjg20151017.

中央級公益性科研院所基本科研業(yè)務(wù)費(fèi)專項(xiàng)(2014B10),國家重大科技專項(xiàng)“cap1400關(guān)鍵審評技術(shù)研究”(2013zx06002001-09),國家自然科學(xué)基金(51378479,51108431)項(xiàng)目聯(lián)合資助.

章旭斌，男，1986年生，博士研究生. 研究方向?yàn)椴▌訑?shù)值模擬.E-mail:zhangxvbin.hebi@163.com

10.6038/cjg20151017

P315,P631

2015-01-12，2015-03-27收修定稿

章旭斌， 廖振鵬， 謝志南. 2015. 透射邊界高頻耦合失穩(wěn)機(jī)理及穩(wěn)定實(shí)現(xiàn)——SH波動.地球物理學(xué)報，58(10):3639-3648,