波動譜元模擬中透射邊界穩定性分析

章旭斌，謝志南

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

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

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

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

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

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

1 數值計算格式

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

2 透射邊界反射系數

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

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

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

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

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

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

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

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

下面采用數值實驗驗證本文得到的反射系數。計算模型長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 采用二階格式，邊界參數S=0.05。

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

圖4 三個觀測點的位移時程以及基于理論計算和數值實驗的二階MTF 反射系數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 透射邊界穩定條件

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

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

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

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

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

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

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

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

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

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

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

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

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

圖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 中菱形)設置在右邊邊界上。入射平面波時間函數采用Ricker 子波，主頻10 Hz。

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

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

圖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 的反射系數，揭示了MTF 引發失穩的機理，同時給出了MTF 穩定條件，并通過數值實驗加以驗證。主要結論有：

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

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