波動譜元模擬中透射邊界穩(wěn)定性分析

章旭斌，謝志南

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

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

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

依據(jù)已有的研究結果，由MTF 引發(fā)的失穩(wěn)形式表現(xiàn)為零頻飄移失穩(wěn)、反射放大失穩(wěn)和高頻振蕩失穩(wěn)，后兩者為高頻失穩(wěn)。LIAO 和LIU[23]針對一維波動集中質量有限元模擬，基于邊界反射系數(shù)揭示了MTF 引發(fā)的反射放大失穩(wěn)機理，即在有限區(qū)域內人工邊界對高頻行波的反復反射放大所致。景立平等[24]針對二維SH 波動集中質量有限元模擬，基于GKS 定理的群速度解釋[25]揭示了高頻振蕩失穩(wěn)的機理，即邊界和內域格式支持群速度指向內域的高頻諧波。依據(jù)對高頻失穩(wě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)機理。而關于穩(wěn)定性分析方法，應提及BEAM 等[31]提出的P 穩(wěn)定(Practical stability)，即在保證GKS 穩(wěn)定條件下，仍需保證數(shù)值解不允許隨時間增長。TREFETHEN[25]論證了P 失穩(wěn)(P-instability)往往表現(xiàn)為反射放大失穩(wě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]并不適用于譜元法。譜元法是高階有限元，其單元內多個非等間距節(jié)點的組合呈現(xiàn)周期延拓，故而存在多條包括真實模態(tài)和虛假模態(tài)在內的頻散曲線，而線性有限元是單一節(jié)點的周期延拓，僅存在一條逼近波動方程的頻散曲線。因此，在譜元法中分析人工邊界反射系數(shù)時需考慮這一特點。本文將依據(jù)譜元節(jié)點周期延拓性質，通過構建譜元和MTF 數(shù)值格式的向量形式，進而基于向量形式推導MTF 反射系數(shù)，從而分析人工邊界穩(wěn)定性，揭示邊界引發(fā)失穩(wěn)的機理，并給出穩(wěn)定條件，為實際波動模擬參數(shù)選取提供參考。同時其研究思路對于其他類型波動和人工邊界穩(wěn)定性分析具有參考價值。

1 數(shù)值計算格式

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

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

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

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

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

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

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

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

圖3(a)為 Δτ=0.1時四階譜元的頻散，其存在4 條頻散曲線，僅其中1 條為具有物理意義的真實模態(tài)頻散(實線)，即對應式(14)的最小特征值[12]。其余3 條均為虛假模態(tài)的頻散，這是高階格式所固有的特性，頻散曲線的數(shù)量由周期延拓的節(jié)點個數(shù)所決定。圖3(b)為MTF 的反射系數(shù)，包括邊界對真實模態(tài)和虛假模態(tài)的反射系數(shù)。從圖中可以看出，虛假模態(tài)的反射系數(shù)較大。雖然虛假模態(tài)對譜元模擬精度影響很小[11]，然而其對人工邊界穩(wěn)定性的影響卻不可忽略。圖3(c)為不同邊界參數(shù)取值下，二階MTF 對真實模態(tài)的反射系數(shù)，可以看出當S=Δτ時 (即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ù)。計算模型長1000 m，密度為2.2 g/cm3，波速為500 m/s。距離左端250 m 處施加Ricker 子波荷載，主頻5 Hz。左端設置MTF，右端為自由邊界，觀測點設置在距離左邊界0 m、50 m 和100 m 的位置。采用四階譜元離散，Δx=10 m ，Δt=0.002 s，此時 Δτ=0.1。MTF 采用二階格式，邊界參數(shù)S=0.05。

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

圖4 三個觀測點的位移時程以及基于理論計算和數(shù)值實驗的二階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)。前者為在有限的計算區(qū)內邊界對外行波的反復發(fā)射放大所致，后者為邊界和內域共同支持了群速度指向內域的平面諧波。實際上這兩類失穩(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)定性時，必須同時避免虛假模態(tài)和真實模態(tài)的反射放大。由反射系數(shù)式(12)，可以得到MTF 在譜元法中的穩(wěn)定條件為：

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

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

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

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

圖7 失穩(wěn)的觀測點位移時程和模擬參數(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)定的觀測點位移時程和模擬參數(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

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

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

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

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

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

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

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

圖11 右邊界上三個觀測點的位移時程Fig. 11 The displacement time history of three receivers located at the right boundary

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

4 結論

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

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

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