章旭斌，謝志南

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

對于大多數(shù)只關(guān)心廣義結(jié)構(gòu)及其鄰近介質(zhì)中波動的問題，其本質(zhì)是無限域波動問題。對無限域波動數(shù)值模擬的關(guān)鍵在于將其轉(zhuǎn)換成有限域的數(shù)值計(jì)算，因此需引入人工邊界截取有限計(jì)算區(qū)，并設(shè)置人工邊界條件模擬外部無限域?qū)ν庑胁ǖ哪芰枯椛?。無限域波動數(shù)值模擬涉及內(nèi)域數(shù)值算法和人工邊界條件，迄今常用的內(nèi)域算法有差分法、有限元法和譜元法[1]等，而人工邊界有粘彈性邊界[2-3]、連分式邊界[4-6]、透射邊界[7-8]和完美匹配層[9-10]等。由內(nèi)域和邊界算法兩者組合，將形成多種多樣的波動數(shù)值模擬方案。

譜元法(Spectral element method, SEM)是一種節(jié)點(diǎn)非均勻分布的高階有限元，其必然存在虛假模態(tài)。虛假模態(tài)也即高階離散格式中存在的虛假解，描述了離散網(wǎng)格中允許的虛假運(yùn)動形式，可由頻散分析得到。MULDER[11]分析了譜元法中的虛假模態(tài)對數(shù)值模擬精度的影響并不明顯。DE BASABE 等[12]分析了譜元頻散和穩(wěn)定條件，闡述了4 階以上的譜元一波長只要4 個(gè)～5 個(gè)節(jié)點(diǎn)，就基本消除了數(shù)值頻散和各向異性。KOMATITSCH和TROMP[13]將譜元法推廣到三維波動數(shù)值模擬，并編寫了SPECFEM3D 程序，推動了譜元法在地震波動數(shù)值模擬中的應(yīng)用[14-16]。由廖振鵬等[7-8]提出的透射邊界亦稱多次透射公式(Multi-transmitting formula, MTF)，是一類典型的具有高階精度、低計(jì)算量、且易于實(shí)現(xiàn)的人工邊界條件，其已在地震工程、電磁學(xué)、流體力學(xué)等數(shù)值模擬領(lǐng)域得到應(yīng)用[17-19]。已有研究將MTF 應(yīng)用于波動譜元模擬中，戴志軍等[20]在譜元法中應(yīng)用了時(shí)間插值和空間插值形式的MTF，顯示后者具有很好的穩(wěn)定性。于彥彥等[21]對比了粘性邊界和二階MTF，MTF 明顯提高了數(shù)值模擬精度。XING 等[22]基于數(shù)值實(shí)驗(yàn)討論了MTF 的精度和穩(wěn)定性，針對不同入射角設(shè)置不同人工波速的MTF，其吸收精度可接近完美匹配層的精度。上述研究顯示了MTF 具有較好的模擬精度和穩(wěn)定性，然而，在波動譜元模擬中MTF 仍可能引發(fā)數(shù)值失穩(wěn)現(xiàn)象，其失穩(wěn)機(jī)理和穩(wěn)定條件尚不明確。下面梳理MTF 穩(wěn)定性研究結(jié)果及其穩(wěn)定性分析方法。

依據(jù)已有的研究結(jié)果，由MTF 引發(fā)的失穩(wěn)形式表現(xiàn)為零頻飄移失穩(wěn)、反射放大失穩(wěn)和高頻振蕩失穩(wěn)，后兩者為高頻失穩(wěn)。LIAO 和LIU[23]針對一維波動集中質(zhì)量有限元模擬，基于邊界反射系數(shù)揭示了MTF 引發(fā)的反射放大失穩(wěn)機(jī)理，即在有限區(qū)域內(nèi)人工邊界對高頻行波的反復(fù)反射放大所致。景立平等[24]針對二維SH 波動集中質(zhì)量有限元模擬，基于GKS 定理的群速度解釋[25]揭示了高頻振蕩失穩(wěn)的機(jī)理，即邊界和內(nèi)域格式支持群速度指向內(nèi)域的高頻諧波。依據(jù)對高頻失穩(wěn)的認(rèn)識，已提出諸多穩(wěn)定措施[26-28]。零頻飄移失穩(wěn)是多種數(shù)值格式(如粘性邊界和Higdon 邊界等)常見的失穩(wěn)現(xiàn)象，通常采用小量修正人工邊界格式[29]或降階方法[30]抑制飄移失穩(wěn)。

以上研究工作基于反射系數(shù)分析和GKS 定理的群速度解釋，揭示了波動有限元模擬中MTF 引發(fā)的反射放大失穩(wěn)和高頻振蕩失穩(wěn)機(jī)理。而關(guān)于穩(wěn)定性分析方法，應(yīng)提及BEAM 等[31]提出的P 穩(wěn)定(Practical stability)，即在保證GKS 穩(wěn)定條件下，仍需保證數(shù)值解不允許隨時(shí)間增長。TREFETHEN[25]論證了P 失穩(wěn)(P-instability)往往表現(xiàn)為反射放大失穩(wěn)，并進(jìn)一步闡明了GKS 定理群速度解釋所針對的失穩(wěn)模態(tài)的反射系數(shù)為無窮大，即對失穩(wěn)模態(tài)而言，其僅有反射波而無入射波。因此，MTF引發(fā)的2 種高頻失穩(wěn)可統(tǒng)一采用反射系數(shù)來分析。

基于反射系數(shù)分析得到的線性有限元中MTF穩(wěn)定條件[23,32]并不適用于譜元法。譜元法是高階有限元，其單元內(nèi)多個(gè)非等間距節(jié)點(diǎn)的組合呈現(xiàn)周期延拓，故而存在多條包括真實(shí)模態(tài)和虛假模態(tài)在內(nèi)的頻散曲線，而線性有限元是單一節(jié)點(diǎn)的周期延拓，僅存在一條逼近波動方程的頻散曲線。因此，在譜元法中分析人工邊界反射系數(shù)時(shí)需考慮這一特點(diǎn)。本文將依據(jù)譜元節(jié)點(diǎn)周期延拓性質(zhì)，通過構(gòu)建譜元和MTF 數(shù)值格式的向量形式，進(jìn)而基于向量形式推導(dǎo)MTF 反射系數(shù)，從而分析人工邊界穩(wěn)定性，揭示邊界引發(fā)失穩(wěn)的機(jī)理，并給出穩(wěn)定條件，為實(shí)際波動模擬參數(shù)選取提供參考。同時(shí)其研究思路對于其他類型波動和人工邊界穩(wěn)定性分析具有參考價(jià)值。

1 數(shù)值計(jì)算格式

進(jìn)而結(jié)合時(shí)間中心差分法，可得時(shí)空解耦的數(shù)值格式：

2 透射邊界反射系數(shù)

在分析離散格式中人工邊界的反射系數(shù)時(shí)，需要知道內(nèi)域格式的頻散關(guān)系，也就是內(nèi)域格式在無限域上支持的諧波頻率和波數(shù)之間的關(guān)系，可通過頻散分析得到。而頻散分析是將內(nèi)域節(jié)點(diǎn)看作在無限域上的周期延拓。由于譜元是一種高階有限元，其以單元內(nèi)多個(gè)節(jié)點(diǎn)的組合呈現(xiàn)周期延拓，如四階譜單元具有5 個(gè)節(jié)點(diǎn)，以1 個(gè)邊界節(jié)點(diǎn)和3 個(gè)內(nèi)節(jié)點(diǎn)進(jìn)行周期延拓(如圖1 所示)，其有別于線性有限元以單一節(jié)點(diǎn)的周期延拓。因此，分析譜元頻散需考慮節(jié)點(diǎn)周期延拓這一性質(zhì)。

圖1 四階譜單元節(jié)點(diǎn)周期延拓示意圖Fig. 1 Periodic extension of four nodes in fourth-order SEM

同樣的，在分析邊界反射系數(shù)時(shí)也應(yīng)考慮節(jié)點(diǎn)周期延拓性質(zhì)，故而構(gòu)建數(shù)值格式的向量形式，即將內(nèi)域譜元周期延拓的節(jié)點(diǎn)組裝成向量形式的內(nèi)域格式，將邊界節(jié)點(diǎn)與鄰近譜單元內(nèi)的節(jié)點(diǎn)組裝成向量形式的邊界格式，從而構(gòu)建了與原格式等價(jià)的向量形式格式，同時(shí)這一向量形式規(guī)避譜元非等間距節(jié)點(diǎn)分布給分析造成的困難。實(shí)際波動數(shù)值模擬通常采用四階譜元和二階MTF，圖2 為邊界節(jié)點(diǎn)附近譜單元節(jié)點(diǎn)分布。依據(jù)MTF格式(8)和內(nèi)域計(jì)算格式(3)，組裝后的邊界格式為：

圖2 人工邊界節(jié)點(diǎn)和四階譜單元節(jié)點(diǎn)組合示意圖Fig. 2 Combination of nodes in artificial boundary and fourth-order SEM

人工邊界對入射波的吸收精度主要體現(xiàn)在反射系數(shù)的模，下文所指反射系數(shù)均表示反射系數(shù)的模。值得說明的是，式(10)中的入射波和反射波由其群速度方向確定，群速度向外的為入射波，群速度向內(nèi)的為反射波[25]。因此，若諧波的相速度向外而群速度向內(nèi)(由頻散關(guān)系確定)，依據(jù)式(12)計(jì)算的反射系數(shù)取其倒數(shù)。

依據(jù)譜單元節(jié)點(diǎn)周期延拓性質(zhì)(圖1)和內(nèi)域節(jié)點(diǎn)格式(3)，組裝后的內(nèi)域格式為：

圖3(a)為 Δτ=0.1時(shí)四階譜元的頻散，其存在4 條頻散曲線，僅其中1 條為具有物理意義的真實(shí)模態(tài)頻散(實(shí)線)，即對應(yīng)式(14)的最小特征值[12]。其余3 條均為虛假模態(tài)的頻散，這是高階格式所固有的特性，頻散曲線的數(shù)量由周期延拓的節(jié)點(diǎn)個(gè)數(shù)所決定。圖3(b)為MTF 的反射系數(shù)，包括邊界對真實(shí)模態(tài)和虛假模態(tài)的反射系數(shù)。從圖中可以看出，虛假模態(tài)的反射系數(shù)較大。雖然虛假模態(tài)對譜元模擬精度影響很小[11]，然而其對人工邊界穩(wěn)定性的影響卻不可忽略。圖3(c)為不同邊界參數(shù)取值下，二階MTF 對真實(shí)模態(tài)的反射系數(shù)，可以看出當(dāng)S=Δτ時(shí) (即ca=c，人工波速取為物理波速)，MTF 的吸收精度最好。

圖3 四階譜元的頻散曲線和二階MTF 的反射系數(shù)Fig. 3 The dispersion curves of fourth-order SEM and the reflection coefficients of second-order MTF

下面采用數(shù)值實(shí)驗(yàn)驗(yàn)證本文得到的反射系數(shù)。計(jì)算模型長1000 m，密度為2.2 g/cm3，波速為500 m/s。距離左端250 m 處施加Ricker 子波荷載，主頻5 Hz。左端設(shè)置MTF，右端為自由邊界，觀測點(diǎn)設(shè)置在距離左邊界0 m、50 m 和100 m 的位置。采用四階譜元離散，Δx=10 m ，Δt=0.002 s，此時(shí) Δτ=0.1。MTF 采用二階格式，邊界參數(shù)S=0.05。

圖4(a)為觀測點(diǎn)位移時(shí)程，100 m 處觀測點(diǎn)入射波和反射波完全分離且間隔時(shí)間較短，可由其入射波和反射波的傅里葉譜比得到MTF 反射系數(shù)。從圖4(b)可以看出，數(shù)值實(shí)驗(yàn)和理論計(jì)算的真實(shí)模態(tài)的反射系數(shù)在0 Hz～10 Hz 范圍內(nèi)基本完全一致。

圖4 三個(gè)觀測點(diǎn)的位移時(shí)程以及基于理論計(jì)算和數(shù)值實(shí)驗(yàn)的二階MTF 反射系數(shù)Fig. 4 The displacement time history of three receivers and the reflection coefficient of second-order MTF computed by theoretical calculation and numerical experiment

3 透射邊界穩(wěn)定條件

波動有限元模擬中MTF 引發(fā)反射放大失穩(wěn)和高頻振蕩失穩(wěn)。前者為在有限的計(jì)算區(qū)內(nèi)邊界對外行波的反復(fù)發(fā)射放大所致，后者為邊界和內(nèi)域共同支持了群速度指向內(nèi)域的平面諧波。實(shí)際上這兩類失穩(wěn)可統(tǒng)一采用人工邊界反射系數(shù)來分析，兩者僅在反射系數(shù)上表現(xiàn)差異，前者為反射系數(shù)大于1，后者為反射系數(shù)無窮大，即僅有反射波，而無入射波[25]。因此，MTF 引發(fā)的兩種高頻失穩(wěn)可統(tǒng)一采用反射系數(shù)來分析。

另外，譜元法中存在虛假模態(tài)，雖然對數(shù)值模擬的精度影響很小，但對數(shù)值穩(wěn)定性不可忽略。因此，在分析邊界穩(wěn)定性時(shí)，必須同時(shí)避免虛假模態(tài)和真實(shí)模態(tài)的反射放大。由反射系數(shù)式(12)，可以得到MTF 在譜元法中的穩(wěn)定條件為：

由于四階譜元并未有解析形式的頻散關(guān)系，故而通過數(shù)值方法給出譜元法中MTF 的穩(wěn)定條件。圖5(a)為一階MTF 穩(wěn)定條件S≤0.195，與Δτ無關(guān)，其中 Δτ取值限定在四階譜元的穩(wěn)定條件上限0.147。圖5(b)為二階MTF 穩(wěn)定條件，其表現(xiàn)為無量綱邊界參數(shù)S=caΔt/Δx和譜元參數(shù)Δτ=cΔt/Δx之 間的關(guān)系，當(dāng) Δτ取 值變大時(shí)，S取值上限變小。一維線性有限元中MTF 的穩(wěn)定條件[23]為S≤1.5，MTF 在高階譜元中的穩(wěn)定性較線性有限元更為復(fù)雜。從圖5 中S=Δτ曲線可以發(fā)現(xiàn)，對于一階MTF 而言，這一取值方式在整個(gè)區(qū)間內(nèi)都是穩(wěn)定的，而對于二階MTF 而言，這一取值方式存在失穩(wěn)的情況，即當(dāng) Δτ取值接近四階譜元穩(wěn)定條件上限0.147 時(shí)。

圖5 四階譜元法中MTF 的穩(wěn)定條件Fig. 5 The stability condition of MTF in SEM

下面通過第2 節(jié)的數(shù)值算例加以驗(yàn)證，取S=0.4，其他參數(shù)不變。圖6(a)中3 個(gè)觀測點(diǎn)的位移時(shí)程在較短時(shí)間內(nèi)出現(xiàn)數(shù)值失穩(wěn)，其中100 m處觀測點(diǎn)的失穩(wěn)時(shí)程與自由面反射波疊加，因此顯得不對稱。圖6(b)為模擬參數(shù)取值下二階MTF 反射系數(shù)，有兩條反射系數(shù)曲線均存在大于1 的值。

圖6 三個(gè)觀測點(diǎn)的位移時(shí)程和模擬參數(shù)條件下二階MTF 的反射系數(shù)Fig. 6 The displacement time history of three receivers and the reflection coefficient of second-order MTF under the given simulation parameters

實(shí)際波動數(shù)值模擬中MTF 人工波速通常取為物理波速，此時(shí)S=Δτ。圖5(b)顯示二階格式S=Δτ=0.147時(shí)，處于失穩(wěn)區(qū)域。同樣采用第2 節(jié)算例加以驗(yàn)證。圖7(a)中的觀測點(diǎn)位移模擬結(jié)果顯示失穩(wěn)。圖7(b)為截取的20 s 數(shù)值結(jié)果顯示了失穩(wěn)形式，其失穩(wěn)峰值在逐步放大。圖7(c)為模擬參數(shù)取值下二階MTF 的反射系數(shù)，其中一條虛假模態(tài)的反射系數(shù)大于1。因此，數(shù)值失穩(wěn)機(jī)理為MTF 對虛假模態(tài)的反射放大及在有限區(qū)域內(nèi)的反復(fù)反射放大所致。該數(shù)值實(shí)驗(yàn)中失穩(wěn)出現(xiàn)的時(shí)間很晚，其原因在于上述取值情況下虛假模態(tài)對應(yīng)的反射系數(shù)沒有遠(yuǎn)大于1，并且數(shù)值模擬中的虛假模態(tài)成分很小，以致其失穩(wěn)過程緩慢。保持其他參數(shù)取值不變，取S=0.137，處于穩(wěn)定范圍。圖8(a)數(shù)值模擬結(jié)果未出現(xiàn)失穩(wěn)，圖8(b)為計(jì)算了500 萬步的最后20 s 位移時(shí)程，未出現(xiàn)微小擾動，沒有失穩(wěn)跡象。圖8(c)為模擬參數(shù)取值下二階MTF 的反射系數(shù)，其值均小于等于1。

圖7 失穩(wěn)的觀測點(diǎn)位移時(shí)程和模擬參數(shù)條件下二階MTF 的反射系數(shù)Fig. 7 The unstable displacement time history of receiver and the reflection coefficient of second-order MTF under the given simulation parameters

圖8 穩(wěn)定的觀測點(diǎn)位移時(shí)程和模擬參數(shù)條件下二階MTF 的反射系數(shù)Fig. 8 The stable displacement time history of receiver and the reflection coefficient of second-order MTF under the given simulation parameters

對于實(shí)際波動模擬而言，本文給出的穩(wěn)定條件過于嚴(yán)格。因?yàn)閷?shí)際波動模擬通常在較短的時(shí)間內(nèi)，而對于不滿足穩(wěn)定條件的S取值，若其反射系數(shù)并未遠(yuǎn)大于1，失穩(wěn)誤差在較短的時(shí)間內(nèi)并不顯著。因此，對于實(shí)際波動模擬，S取值可適當(dāng)放寬。

本文給出了一維波動譜元模擬中MTF 的穩(wěn)定條件，可作為多維波動穩(wěn)定條件的參考。下面以二維SH 波動模擬為例，闡述MTF 的穩(wěn)定性。對于二維平面波入射問題(如圖9)，為了吸收大角度入射波，通常需要增大MTF 的人工波速取值。若平面波入射角為θ，SH 波速為cs，則右側(cè)MTF 的人工波速可取為：

圖9 平面波入射示意圖Fig. 9 Diagram of plane wave incident

由式(16)可知，當(dāng)入射角很大時(shí)，ca會取得很大，以致S=caΔt/Δx取值過大可能引發(fā)數(shù)值失穩(wěn)。

圖10 為盆地-基巖半空間模型，計(jì)算模型長度4 km、深度2 km，盆地長度2 km、深度0.3 km，基巖介質(zhì)密度 ρ=2.5 g/cm3，cs=2.5 km/s，盆地介質(zhì)密度 ρ=1.5 g/cm3，cs=0.5 km/s。采用四階譜元離散，單元平均尺寸 Δx=Δy=50 m，時(shí)間步長Δt=0.0014 s ， 此時(shí) Δτ=0.07。除自由地表外，其余三邊均采用二階MTF，左邊和底邊MTF 的人工波速取為基巖SH 波速，右邊MTF 的人工波速按式(16)取值。3 個(gè)觀測點(diǎn)(圖10 中菱形)設(shè)置在右邊邊界上。入射平面波時(shí)間函數(shù)采用Ricker 子波，主頻10 Hz。

圖10 盆地-基巖半空間模型Fig. 10 Basin-rock half space model

圖11(a)為 θ=80°時(shí)計(jì)算的觀測點(diǎn)位移時(shí)程，MTF 引發(fā)數(shù)值失穩(wěn)，其中MTF 人工波速依據(jù)式(16)計(jì)算為ca≈14.4 km/s ，此時(shí)S≈0.403。參考一維穩(wěn)定條件S≤0.243 ， 可得ca＜8.6 km/s。圖11(b)為 θ=80°時(shí)，取ca=8.0 km/s計(jì)算的觀測點(diǎn)位移時(shí)程未出現(xiàn)數(shù)值失穩(wěn)。圖12 為各個(gè)時(shí)刻的波場圖，從1.2 s 波場圖可以看出右側(cè)邊界MTF 基本沒有反射誤差。對于大角度平面波入射，MTF 人工波速取值可參考一維穩(wěn)定條件，可適當(dāng)調(diào)大以提高對大角度入射波的吸收效果，但取值不宜過大，以免引發(fā)數(shù)值失穩(wěn)。

圖11 右邊界上三個(gè)觀測點(diǎn)的位移時(shí)程Fig. 11 The displacement time history of three receivers located at the right boundary

圖12 譜元和MTF 模擬的各個(gè)時(shí)刻波場Fig. 12 Propagation of wave simulated by SEM with MTF

4 結(jié)論

本文分析了譜元法中MTF 的反射系數(shù)，揭示了MTF 引發(fā)失穩(wěn)的機(jī)理，同時(shí)給出了MTF 穩(wěn)定條件，并通過數(shù)值實(shí)驗(yàn)加以驗(yàn)證。主要結(jié)論有：

(1) 譜元法中虛假模態(tài)雖然對數(shù)值模擬的精度影響很小，但會引發(fā)人工邊界數(shù)值失穩(wěn)。譜元法中MTF 引發(fā)高頻失穩(wěn)的機(jī)理為在有限計(jì)算區(qū)內(nèi)邊界對虛假模態(tài)的反復(fù)反射放大所致；

(2) 一維波動譜元模擬中MTF 穩(wěn)定條件表現(xiàn)為譜元參數(shù)和邊界參數(shù)之間的關(guān)系，可作為多維波動穩(wěn)定條件的參考；(3)通過構(gòu)建高階譜元格式和MTF 格式的向量形式，進(jìn)而基于邊界反射系數(shù)是分析MTF 穩(wěn)定性的可行途徑。關(guān)于多維彈性波動譜元模擬中MTF 的穩(wěn)定性，以及其他類型內(nèi)域格式和人工邊界的穩(wěn)定性仍需深入研究。