透射邊界條件在波動譜元模擬中的實現：一維波動1)

邢浩潔 李鴻晶

(南京工業大學土木工程學院，南京 211816)

透射邊界條件在波動譜元模擬中的實現：一維波動1)

邢浩潔 李鴻晶2)

(南京工業大學土木工程學院，南京 211816)

多次透射公式(multi-transmitting formula，MTF)是一種具有普適性的局部人工邊界條件，但其在近場波動數值模擬中一般與有限元法結合.由于波動譜元模擬的數值格式與有限元格式有極大的不同，傳統的MTF在譜元離散格式中無法直接實現.為了使物理概念清楚、精度可控的多次透射人工邊界條件能夠適應波動譜元模擬的需求，首先指出多次透射邊界與譜元離散格式結合的基本問題，并分析了空間內插和時間內插兩種方案的可行性.然后從空間內插角度出發，提出基于拉格朗日多項式插值模式的MTF譜元格式，并采用一種簡單內插方法實現高階MTF.最后通過一維波動數值試驗檢驗這些MTF譜元格式的精度，并討論其數值穩定性.結果表明：對于一、二階MTF，幾種格式的精度相當；對于三、四階MTF，基于譜單元位移模式插值的格式精度最高.相反，隨著插值多項式階次的升高，不同MTF格式的穩定臨界值逐步降低，但是所有格式均在人工波速大大超過物理波速時才可能發生失穩.

波動數值模擬，譜元法，多次透射邊界，MTF，數值穩定性，數值精度

引言

近場波動數值模擬需要從實際的無限介質中截取有限模型進行分析.為保證波動模擬的有效性，通常假定所有波源或散射體、不規則和不均勻區域都被包括在有限模型內，而僅用模型邊界上的人工邊界條件來模擬外部無限介質對計算區域的影響.因此，它需要處理兩個重要問題：一是計算區域內的波場分布和介質幾何、物理參數的空間離散；二是模型邊界上的人工邊界條件[1-2].針對具體問題采用適當的空間離散方法和人工邊界條件.在保證計算精度和數值穩定性的前提下，最大程度地提高計算效率，是近場波動數值模擬技術追求的目標.

在空間離散方法中，有限差分法[3]和有限元法[4-5]是近場波動數值模擬中使用最廣泛的兩種.但由于數值穩定性方面的限制，它們均采用了低階格式的數值方案，數值精度和計算效率都不是很高.譜方法 (spectral method)[6]具有無窮階收斂性的優點，結合有限單元概念發展起來的譜元法 (spectral element method,SEM)[7]，既保持了譜方法的指數收斂率，又體現了有限元法的靈活性.它不僅能夠以較少的節點模擬復雜波場和不規則幾何形狀，還可以有效地處理不均勻介質，并能夠有效地實現并行計算，為求解大型波動問題提供了強有力的支持[8-9].近年來，應用譜元法實現波動數值模擬開始受到國內外學者的重視，在地震波傳播問題研究方面取得了一些進展[10-23].在上述波動譜元模擬研究中，人工邊界條件一般采用黏性邊界[24]、CE邊界[11,13,25]或完美匹配層邊界(perfectly matched layer，PML)[26-30]，而多次透射邊界作為一種具有普適性的人工邊界卻并未引起重視.由多次透射公式(multi-transmitting formula,MTF)定義的透射邊界條件[31-32]是一種直接以離散形式給出的人工邊界，它通過模擬各種單向波動的共同運動學特征建立邊界條件，具有公式簡單、施加方便、復雜波場模擬精度高等優點，且其高階形式易于與內域離散格式相結合.但目前MTF一般只能與內域的有限元格式相結合完成近場波動的數值模擬，而譜元法的單元模型數值格式與有限元法完全不同，在譜元格式中不能直接套用同有限元法相結合的MTF格式，因而需要按照譜元格式的要求專門研究與其相適應的MTF的具體格式.

實際上，透射邊界條件所具有的精度可控、易于實施以及善于處理復雜波場的優點，與譜元法的理念是高度切合的，因此，在波動譜元模擬中引入透射邊界條件，可能是提高近場波動數值模擬精度和計算效率的一條值得重視的途徑.本工作旨在解決透射邊界條件與譜元離散格式結合的關鍵問題，即提出一種適用的MTF譜元格式，確保其能夠在離散意義上穩定地實現多次透射以達到提高精度的效果. MTF與譜元離散格式結合的方式，是將MTF公式中涉及到的計算點位移用鄰近的譜元離散點位移進行插值表示的過程.最直接的考慮是將目前廣泛使用的MTF有限差分或有限元格式推廣到譜元法中，例如文獻[33]所做的工作.這一推廣包括兩個部分，一是將實現一階MTF的三點拋物線內插公式由等距節點推廣到不等距節點；二是將原來實現高階MTF的齊次內插遞推方案，替換為一種簡單內插方案.研究表明，這種基于譜元不等距節點進行三點拋物線內插的MTF譜元格式，在通常情況下不失為一種可供選擇的方案，但是對于較為復雜的波動情形，其模擬效果往往并不理想.本文將從透射邊界條件的數學本質和譜元離散格式的特定形式出發，系統探討MTF與譜元離散格式結合的相關問題，包括空間內插方案與時間內插方案的含義、特點及適用性，MTF插值多項式階次的選取，高階MTF的實現方案等.通過數值試驗和理論分析總結高精度方案應當具有的特點，在幾種不同的數值方案中，推薦一種能夠較好地適應復雜波動情形的MTF譜元格式，并分析該格式在一維波動問題譜元求解時的精度和穩定性.

1 基本問題

本文以一維波動問題為例，探討在譜元離散格式中實現透射邊界條件的方法.之所以選擇從一維波動問題開始研究，主要基于如下兩點考慮.

(1)一維條件下得到的多次透射邊界數值格式可以不加修改地直接用于二維、三維模型，因為后兩者實現MTF的方式也是“一維化”的，即只涉及到每個邊界節點上一條指向內域的離散網格線上的節點.

(2)一維波動模型摒除了二維、三維模型的諸多復雜性，有助于獲得對MTF離散格式的精度和穩定性較為清晰的認識，而一維波動問題的求解方法和數值規律，可為二維、三維波動的模擬奠定基礎.

1.1 一維波動譜元模型

對于一維無限或半無限長直桿中無阻尼波的傳播問題，在適當的位置截取人工邊界得到有限計算模型.模型中質點的運動由如下波動方程控制

式中，u=u(x,t)為待求的波動位移函數，c為介質波速，f(x,t)為輸入場.

采用具有高精度集中質量矩陣的勒讓德譜元法(Legendre spectral element method,LSEM)，對波動方程進行空間離散.主要步驟包括從等效積分“弱”形式出發、劃分互不重疊的單元、將關于每個物理單元的計算變換到標準的參考單元上進行，最后利用伽遼金原理得到空間離散的運動方程

式中，u為譜元離散模型的節點位移向量，Me,Ke,fe分別為單元的質量矩陣、剛度矩陣和等效節點向量.

譜元法與有限元法的區別在于：單元位移模式和計算單元特性矩陣的數值積分方法不同.波動模擬選擇的有限單元通常為線性單元，只用到兩個端節點，數值積分采用高斯積分.勒讓德譜單元則為高階單元，除了兩個端節點之外還有內節點，單元節點位置根據GLL(Gauss Lobatto Legendre)積分點的相對位置來確定，數值積分采用GLL積分.具體過程見文獻[1,34]，這里僅討論與文中主題高度相關的譜單元位移模式，以及譜單元和有限單元的空間離散尺度.

譜單元位移模式定義在標準的參考單元 Λ ∈[-1,1]上，單元內任意位置的位移由下式表示

其中，ng為單元節點數，ng=NE+1，NE為單元階次；ξi為GLL節點坐標，Ni(ξ)為節點i的形函數.表1列出了部分單元階次NE下的GLL節點坐標.形函數為定義在GLL節點上的拉格朗日插值函數

表1 LSEM在參考單元上的節點坐標Table 1 Node coordinates of LSEM in reference element

由于GLL節點為勒讓德多項式的極值點，是根據對勒讓德多項式的數值求解確定的，因此該形函數也被稱為勒讓德基函數.

各個物理單元的節點位移與相應參考單元的節點位移ue(ξi)是一一對應的，單元節點坐標通過轉換關系來確定，為1號單元節點的坐標，?xE為譜單元尺寸.在等效積分方程的實際計算中，通過雅可比行列式|J|=dx/dξ=?xE/2實現積分區間的變換.

高精度譜單元相對于有限單元的優勢在于數值頻散低，能夠以較少節點模擬波場.對于感興趣的最短波長λmin=c/fmax的空間離散，有限元通常采用10個單元左右[1]，而四階譜單元則只需采用一個單元[10]，如圖1所示.

圖1 有限單元與譜單元對波場的空間離散Fig.1 Spatial discretization of wave fiel using finit element or spectral element

1.2 多次透射公式

有限模型人工邊界節點的位移由多次透射公式計算，這是一種直接以離散形式給出的一維化的時、空外推公式.MTF涉及的外推節點如圖2所示的空心圓點1a,2a,3a,···，它們以時、空步距(?t,ca?t)從邊界節點逐步向內域延伸.其中ca為人工波速，這些外推節點被稱之為邊界計算點.MTF定義式為

圖2 有限元、譜元離散網格中的MTFFig.2 MTF in the discrete grid of finit element and spectral element

圖2(a)給出了目前廣泛使用的MTF有限元格式[1,31-32]所涉及的離散網格點，1a,2a,3a,···點的位移分別由3,5,7,···個有限元節點的位移插值得到.這是一種基于三點拋物線內插公式計算1a點位移，再由一種齊次內插遞推過程計算2a,3a,···各點位移[1]的數值格式.這種格式是由時、空移動算子表示的MTF二項式形式[37]得到的，其優點在于各階MTF的反射系數為R=-fN，f為一階MTF的反射因子[38]，使得在滿足穩定條件|f|≤1時，隨著透射階次N的增加，反射系數能夠穩步地減小.不過，二項式遞推只能在等距節點下實現，因此該MTF數值格式要求涉及到的有限元節點為等距離分布.

圖2(b)表示譜元離散網格中的MTF，與有限元相比，譜元節點之間的間距要大得多且為不等距分布，此時點1a,2a,3a,···的位移可能需要采用與有限元不同的插值方式來計算.可供選用的插值方式有多種，對應不同的MTF譜元格式，關鍵在于確定能夠較好地滿足精度和穩定性要求的格式.關于不同MTF數值格式對精度和穩定性的影響，相關研究幾乎未曾見到，目前所有討論都是基于上述格式進行的[35-45]，它僅是一種有限元等距節點下的MTF數值格式.研究發現，不同MTF數值格式會受到插值公式本身和內域網格頻散效應兩方面的影響，在精度和穩定性方面表現出一定差異.與內域離散格式結合的MTF數值格式和MTF定義式的區別在于，前者并不總能穩步地實現多次透射，達到提高精度的效果，而這將是判斷MTF數值格式好壞的重要標準.

2 MTF譜元格式

Liao等[32]在提出多次透射公式時就指出：多項式、三次樣條、三角函數等各種插值方法，甚至時間步內插法，都可以嘗試作為MTF的數值格式，只是每種格式都存在需要進一步研究的數值誤差問題.文獻[1]對于MTF的有限元格式，除了上述常用的三點拋物線內插和齊次內插的MTF有限元格式之外，還提出了簡單內插、時間內插兩種思路.文獻[33]對于MTF的譜元格式，提出了空間域回退、時間域回退兩種方式.我們研究發現，在譜元離散網格下選擇MTF插值方式是有規律可循的，如：從MTF的本質以及波動數值計算的具體實踐考慮，時間內插方案并不適用；從單元位移模式考慮，多項式插值是最簡單可能也是最好的選擇.

2.1 時空內插方案的選取

理論上，在離散格式中實現MTF可以采用空間內插與時間內插兩類方案.針對譜元離散格式，我們考察了這兩類方案的具體含義，并探討了它們的可行性.

空間內插方案如圖2(b)所示，其基本特點是：MTF計算點1a,2a,3a,···與內域網格點的時間坐標相同，空間坐標不同，此時插值只需在空間方向進行.觀察MTF的定義式(5)不難發現，采用這類插值方案是十分自然的，如圖2(a)作為目前唯一廣泛使用的MTF有限元格式，就是一種空間內插方案.隨后是譜元離散網格下具體插值格式的選擇，這將在下一節進行專門探討.

時間內插方案如圖3所示，其基本特點為：MTF計算點1a,2a,3a,···與內域網格點的空間坐標相同，時間坐標不同，插值只需在時間方向進行.這種方案與文獻[1]介紹的有限元網格下的時間內插思路類似，它無法從MTF定義式(5)得到，而是基于另一種形式的MTF表達式

式中各個 MTF計算點所在的時刻為p+1-sj/ (ca?t)，sj為第j個譜元節點與邊界節點的距離，顯然它們內域節點對應的時刻不一致.

圖3 譜元離散網格的MTF時間內插方案Fig.3 Temporal interpolation scheme of MTF in the discrete grid of spectral element

目前關于MTF時間內插方案的研究很少，僅在文獻[33]中見到，其模擬結果顯示時間域插值的精度和穩定性不如空間域插值.我們對MTF時間內插方案進行了研究，從波動數值模擬的具體實踐并結合對多次透射公式本質的認識，認為譜元離散網格的MTF時間內插方案存在以下問題：

(1)式(6)的有效性缺乏理論支持.與圖2相比，圖3中的計算點1a,2a,3a,···間隔要大得多且為不等距分布.間隔大導致的問題是計算點所在范圍內的波形起伏變化加大(如圖1)，此時通過少量計算點推算邊界節點位移的準確性會下降.不等距分布導致的問題是導出的MTF定義式(5)時所使用的時空差分原理不再適用于止，因此，僅僅仿照式(5)的形式寫出的式(6)，在增加透射階次以提高精度方面缺乏嚴格的理論支持.

(2)邊界計算點1a,2a,3a,···跨越的時間步較多.如圖3所示，1a點跨越3個時間步，2a點跨越7個時間步，3a點跨越11個時間步，這相當于三階MTF計算公式采用的時間步長約為3?t、7?t和11?t，計算精度較低.另外，時間步數大大超過了內域時間積分通常采用的一個時間步(Newmark法)或兩個時間步(中心差分法)，不僅需要在初始時刻對邊界進行專門處理，還可能導致較為嚴重的穩定性問題.

(3)具體的插值公式難以確定.圖3中計算點1a,2a,···附近的內域節點比較密集，可以用來插值的內域節點數目以及節點位置存在多種組合，插值公式難以確定.

上述問題表明，時間內插方案不宜作為研究MTF譜元格式的出發點.接下來，將從空間內插方案出發，導出具體的MTF譜元格式.

2.2 插值格式的確定

為了從空間內插角度得到MTF的譜元格式，首先考慮套用傳統MTF有限元格式的可行性.根據1.2節內容不難得知，直接套用的做法行不通，理由是：(1)等距節點的三點拋物線內插可以推廣到譜元不等距節點，但二者的插值系數并不相同；(2)高階MTF計算點的齊次內插遞推過程無法在不等距節點下實現.

于是，我們決定采用不等距節點的三點拋物線內插公式在譜元網格中實現一階MTF，并借鑒文獻[1]介紹的簡單內插思路，即對點2a,3a,···采用與點1a完全相同的插值格式，在譜元網格中實現高階MTF.由此得到一種基于三點拋物線內插的MTF譜元格式，如圖4所示.

圖4 基于三點拋物線內插的MTF譜元格式Fig.4 MTF scheme based on parabolic interpolation applied in spectral element method

該格式人工邊界節點0在p+1時刻的位移由下式計算

其中，上標T表示向量的轉置.式(8)中的插值系數tj,0,tj,1,tj,2與MTF計算點ja(j=1,2,···,N)以及式(9)中的譜元網格點是一一對應的，其數值由拉格朗日插值公式確定.

其中，s0(s0=0)，s1和s2為離散網格點與邊界節點的距離(圖4)，可根據譜單元尺寸和表1給出的GLL節點坐標進行計算.

從圖4可以看出，式(7)～式(10)表示的MTF譜元格式需要滿足各個計算點1a,2a,3a,···位移均由三點拋物線“內插”公式計算的條件，即

上式對能夠實現的MTF階次作了一定限制.根據內域時間積分的穩定條件以及譜元節點位置關系不難得出，當人工波速等于物理波速時(實際計算時通常如此)，能夠實現的MTF階次一般能夠達到3以上，基本滿足使用要求.但是，若人工波速取值較大且時間步長接近穩定臨界值時，則需要驗算MTF階次是否滿足上式要求.

數值試驗表明，一般情形下該MTF譜元格式能夠有效地模擬外行波在人工邊界上的透射過程，但對于比較復雜的情形，如人工波速與介質物理波速差異較大時，該格式往往難以隨著透射階次的提高而較好地收斂到理想結果.分析其原因，可能是三點拋物線內插格式的插值多項式階次僅為2，與常用的譜單元階次(如4～8)差距較大，導致精度不足.考慮到復雜二維、三維波動問題的模擬要求，有必要開發精度更高的MTF譜元格式.

為改善邊界精度，提高插值多項式階次，即利用更多的譜元節點進行插值，得到新的MTF譜元格式.若采用的插值多項式階次為M，則稱之為基于M次多項式插值的MTF譜元格式，如圖5所示.

圖5 基于M次多項式插值的MTF譜元格式Fig.5 MTF scheme based onMth-order polynomial interpolation applied in spectral element method

此時人工邊界節點0在p+1時刻的位移計算公式為

式(13)中的插值系數tj,0,tj,1,···,tj,M可根據拉格朗日插值多項式計算

離散網格點與邊界節點的距離sk(k=0,1,···,M)，根據譜單元尺寸和表1給出的GLL節點坐標進行計算.

提高后的插值多項式階次M的取值可以為3,4,···,NE，NE為譜單元階次(這里不考慮超過譜單元階次的情形)；當M取值為2時，式(12)～式(15)與式(7)～式(10)相同.因此，本文MTF譜元格式可概括為基于M次多項式插值的格式，M取值從2到NE.

對于式(12)～式(15)表示的MTF譜元格式，保證各個計算點1a,2a,3a,···位移的插值公式均為“內插”的條件是

上式與式(11)類似，對通常使用的MTF階次不構成限制，只是在人工波速大大超過介質物理波速或透射階次很高時，才需要進行驗證.

數值試驗表明，提高插值多項式階次M能夠有效地提高本文MTF譜元格式的模擬精度，尤其是能夠改善低階插值格式在高階MTF方面的不足，增強模擬復雜波場的能力.那么，在具體的波動譜元模擬中，插值多項式階次M取多少最為合適？這個問題難以單純地從對數值試驗結果的總結來回答，還必須結合數值分析的基本理論進行分析.分析認為：當階次很高時等距節點下的高階多項式插值會出現龍格現象，導致精度降低；譜單元采用不等距的GLL節點，能夠有效地避免龍格現象的出現，具有很高精度.不過，插值精度是以單元為單位來確定的，由NE階譜單元組成的離散網格中節點的最高插值精度為2NE-1階.因此，插值多項式階次取M=NE最為合適.

對比式(15)與式(4)不難發現，當插值多項式階次取為M=NE時，式(12)～式(15)的MTF譜元格式與內域的單元位移模式是一致的，這似乎暗示了某種一般性原則.需要指出的是，一般情況下，插值多項式階次低于譜單元階次的邊界格式同樣具有較好模擬效果，不同格式之間的差異只有在模擬復雜波動問題時才表現得較為明顯.

3 精度

通過一維半無限均勻彈性直桿中波傳播問題的數值試驗來檢驗本文MTF譜元格式的精度.計算模型如圖6所示，左端為輸入端，人工邊界位于距桿端L=200 m處，桿中波速為c=200 m/s.內域采用勒讓德譜單元離散，單元階次NE=5，時間積分采用中心差分法，人工邊界采用式(12)～式(15)表示的MTF譜元格式.

圖6 半無限均勻彈性直桿的譜元離散模型Fig.6 Discretized spectral element model for a semi-infinit straight uniform elastic rod

輸入波采用三次樣條脈沖波S(t)，幅值1 m，脈沖非零段時間寬度Ts，表達式為計算時取T=0.2 s，根據傅里葉分析可知其上限頻率約為14 Hz，此時桿中最短波長λmin≈14.3 m.模型的單元數目取為14，則單元尺寸?xE與λmin一致，符合譜元離散的精度要求[10,13].時間步長取為?t=0.002 s，滿足時域積分穩定條件[34].分別使用插值多項式階次為M=2～5對應的四種MTF譜元格式進行模擬，并將它們記為interp 2，interp3，interp 4和interp spec.

人工邊界條件是通過一種假定的外行波動模式來推算邊界節點位移的，該波動模式與實際外行波動的接近程度決定了邊界的精度.多次透射邊界假定的波動模式由人工波速(ca)和透射階次(N)兩個基本參數確定，其優點在于高精度和靈活性，具體為：人工波速取值適當時，較低的透射階次就能達到很好的模擬效果；人工波速與介質物理波速差異較大時，增加透射階次同樣能夠逐步提高模擬精度.為不失一般性，這里考慮人工波速等于、大于和小于介質物理波速三種情形.

對于人工波速等于介質物理波速的情形，即ca=c，四種MTF譜元格式計算得到的邊界節點位移時程如圖7所示.

此時對于這四種MTF譜元格式中的任何一種，采用一階MTF得到的邊界位移結果已經十分接近精確解.給出高階MTF計算結果的目的在于觀察當透射階次N增加時，每種MTF譜元格式能否進一步提高或者至少保持相當于一階MTF的精度.圖7反映了MTF數值格式與MTF定義式的不同，前者并不總是能夠穩步地達到增加透射階次提高邊界精度的效果.從圖7(a)可以看出，在提高MTF階次后插值階次與單元階次差異較大的interp 2格式，不僅沒有提高精度，反而誤差越來越大.圖7(b)和圖7(c)表明，插值階次接近單元階次的interp 3和interp 4格式，其高階MTF在提高精度方面仍未達到效果，但是在控制誤差方面有所改善.圖7(d)結果表明，與單元位移模式一致的interp spec格式，始終能夠保持數值解十分接近于精確解，它在控制高階MTF誤差方面效果最好.

多次透射公式的優點在于：先以一個統一的人工波速[1,31-32]來描述外行波沿邊界“法線”的視傳播，再通過多次透射的方法逐步消除由于人工波速與實際視傳播速度不同而造成的誤差.就一維波動模擬而言，外行波的視傳播速度就是介質物理波速，實際計算時人工波速應當取為介質物理波速.這里純粹從研究角度出發，采用不同人工波速進行模擬，檢驗本文MTF譜元格式的精度.

圖7 人工邊界節點位移時程(ca=c)Fig.7 Displacement history of the artificia boundary node(ca=c)

對于人工波速大于介質物理波速的情形，取ca=2c，四種MTF譜元格式計算的邊界節點位移時程如圖8所示.

圖8結果顯示四種MTF譜元格式在一階和二階MTF時的模擬結果非常接近，在三階和四階MTF時的模擬結果則表現出明顯差異.一階MTF時幾種格式的模擬結果都存在較大誤差，完全不能滿足要求，二階MTF能夠迅速減小誤差，使計算結果接近于精確解.三階和四階MTF時，圖8(a)中interp 2格式呈現誤差增大的趨勢，圖8(b)和圖8(c)中interp 3和interp 4格式體現了提高插值階次能夠減小高階MTF誤差的作用，圖8(d)中interp spec格式在高階MTF時的誤差最小，能夠保證高階MTF計算結果的收斂性，相比于前三種格式的優勢比較明顯.

圖8 人工邊界節點位移時程(ca=2c)Fig.8 Displacement history of the artificia boundary node(ca=2c)

圖9 人工邊界節點位移時程(ca=0.5c)Fig.9 Displacement history of the artificia boundary node(ca=0.5c)

對于人工波速小于介質物理波速的情形，取ca=0.5c，四種MTF譜元格式計算得到的邊界節點位移時程如圖9所示.

圖9結果顯示四種MTF譜元格式在一階和二階MTF時的模擬結果基本看不出差別，在三階和四階MTF時的模擬結果則出現一定差異.對于四種MTF譜元格式，一階MTF數值解都大幅度低于精確解，誤差很大，無法滿足要求；二階MTF能夠正常地減小誤差，使結果靠近精確解；但到三階和四階MTF，interp 2格式的精度不再提高，interp3和interp 4的結果慢慢地向精確解靠近，interp spec能夠穩步地逼近精確解.

四種MTF譜元格式在一維波動模擬中的精度可總結為：(1)在一階和二階MTF時，四種格式的精度相當，對于人工波速與實際“法向”傳播速度(一維情形為介質物理波速)差別不大的波動問題，通常二階MTF已能夠滿足精度要求，四種格式都可以使用.(2)在三階和四階MTF時，能夠很好地實現增加MTF階次以提高精度效果的是interp spec格式，即與譜單元位移模式一致的格式.因此，對于人工波速與實際“法向”傳播速度差別較大的情形，如一維波動問題中人為選取大于或小于介質物理波速的人工波速，或二維、三維波動問題中受透射角度或介質交界面影響導致視傳播速度難以準確定義時，推薦采用與單元位移模式一致的MTF譜元格式.

關于外行波“法向”視傳播速度與介質物理波速的關系，一維波動情形下二者相同，而在二維或三維波動情形下二者明顯不同，如：存在透射角度會導致視傳播速度大于介質物理波速，且透射角度越大二者差異越大；P-SV波動問題中P波與SV波波速的差異巨大，加上透射角度等因素的影響，導致視傳播速度變得比較復雜，難以用一個值來描述.因此，盡管一維波動算例的模擬結果對二維或三維波動問題有一定參考價值，但它們并不是完全等同的.

4 穩定性

局部人工邊界條件用于時域逐步積分計算時存在穩定性問題[35]，式(12)～式(15)的MTF譜元格式也不例外.討論人工邊界的穩定性，必須與具體的內域算法相結合[36-37]，波動問題的類型、維數、MTF數值格式、透射階次等多種因素都可能對穩定性造成影響.為使問題不至過于復雜，這里僅就簡單的一維波動情形，初步討論上述MTF譜元格式的穩定性.

文獻[38]在頻域內推導了一維波動有限元(或有限差分)離散模型中傳統MTF譜元格式的反射系數精確解，并據此闡述了失穩機理為：當反射系數大于1時才可能發生失穩；失穩來自邊界對波動有限元(或有限差分)模擬無意義的高頻段的反射放大；高頻誤差經由離散網格，在人工邊界和物理界面之間不斷反射，每至人工邊界就被放大一次，當波動循環次數足夠多的，累積效應導致出現振蕩失穩.后來，文獻[39]提出一種直接在時域內分析MTF穩定性的方法：基于傳遞矩陣譜半徑的穩定性判別法，該方法得到的結果與頻域方法幾乎完全一致.對于本文MTF譜元格式而言，使用以上兩種分析方法都比較困難，前者是因為難以得到不等距節點和高階單元位移模式的頻域解，后者因為當傳遞矩陣的最大特征值的模是1或非常接近于1時[40]，會引起較大的誤差，導致判斷結果不夠準確.但是，從離散網格相當于低通濾波器的觀點來看，譜元網格和有限元網格對波動的作用相似，區別僅在于截止頻率不同.由此不難推斷，本文MTF譜元格式的失穩機制應當與有限元或有限差分類似，為一種高頻振蕩失穩，數值試驗結果驗證了該推斷.

為消除MTF的高頻振蕩失穩，人們提出了多種方法，主要包括：對邊界計算區內的節點位移進行三點平滑[38]；在整個計算區內施加與應變速度成正比的黏性阻尼[39]；在邊界區設置黏性阻尼[41-42]；調整內域離散格式的網格參數[43]；采用低階MTF與高階MTF組合的形式[44]；在不降低精度的前提下修改內域有限元格式的剛度矩陣[45].但值得注意的是，上述措施都被用于二維或三維模型.那么，一維波動問題是否需要采用穩定實現MTF的措施？

一維波動數值模擬中極少出現失穩現象，因此一般認為不需要采取穩定措施，但這樣的直觀判斷并不能令人完全滿意.實際上，有研究已經十分接近從理論上回答這個問題，只是之前側重于從頻域反射系數的角度解釋MTF的失穩機理，而忽視了失穩現象與時域計算參數之間的聯系.文獻[46]在文獻[38]工作的基礎上，解析地證明了如下命題：對于一維波動有限元離散模型，若透射邊界的人工波速取值大于1.5倍的空間步距與時間步距的比值，則在某一高頻段內其穩態波動解在人工邊界上反射系數的模大于1.這一命題的含義為：就一維波動有限元離散模型而言，若ca為人工波速，?x,?t分別為空間步距和時間步距，則只有在滿足條件

時，MTF才可能出現高頻失穩.反過來說，當人工波速小于或等于1.5倍空間步距與時間步距的比值時，MTF不會發生失穩.對應于MTF的頻域穩定條件(即反射系數的模大于1)，將式(18)稱為MTF的時域穩定條件.

為論述方便，定義兩個參數：波速比，指人工波速與實際法向透射速度(一維為物理波速)的比值，用符號α表示，α=ca/c；無量綱時間步距，指時間步長與物理波速的乘積再除以空間網格尺寸，用符號?τ表示，?τ=c?t/?x.這樣，上述MTF時域穩定條件可表示為一種更為簡潔的形式

上式的含義如圖10所示(?τ的取值范圍由內域計算的穩定條件確定，為?τ≤1)，其直觀呈現出MTF時域穩定條件的價值在于：在一維波動有限元離散模型中，可以通過控制人工波速來保證透射邊界的穩定性.因此，可以從理論上解釋一維波動數值模擬很少出現失穩現象的原因，即人工波速取值不夠大.例如，當?τ=1時，取α＞1.5才可能失穩；當?τ=0.5時，取α＞3才可能失穩；當?τ=0.2時，至少取α＞7.5才可能失穩.

上述MTF時域穩定條件是針對一維波動有限元模型的，其關鍵在于數字1.5的確定，我們稱之為MTF的穩定臨界值.文獻[46]通過求解頻域反射系數的模大于1的不等式來確定有限元的MTF穩定臨界值.對于一維波動譜元模擬而言，大量數值試驗表明也存在類似的MTF穩定臨界值.但是，若要從頻域角度推導其MTF反射系數，并進一步求解不等式，是極其困難甚至難以實現的，本工作暫不討論.我們采用另外一種方法來確定MTF譜元格式的穩定臨界值，即根據高頻失穩現象比較顯著，能夠在波動數值模擬結果中輕易地被觀察的特點，采用數值算例進行大量試算，確定MTF譜元格式的穩定臨界值.試算方法的準確性已通過有限元模型得到驗證.

對于一維波動譜元模型，波速比α的定義與前文相同，而無量綱時間步距則定義為?τ=c?t/s1，s1為譜單元端部兩個節點之間的距離.若用 γc表示MTF譜元格式的穩定臨界值，則有相應的MTF時域穩定條件為

采用上一節的一維波動算例進行試算，計算時間設定為300 s，以觀察到明顯的高頻振蕩現象作為失穩標準.對于上一節的四種MTF譜元格式，考慮一階MTF，在計算過程中嘗試α和?τ的不同組合，發現總是當α?τ超過一定值時，對應的MTF才出現失穩.得到四種MTF譜元格式下，一階MTF的穩定臨界值如表2所示.

表2 不同MTF譜元格式的穩定臨界值(一維波動)Table 2 Stability thresholds of several MTF spectral element schemes(1-D wave motion)

表2顯示從interp 2格式到interp spec格式，即MTF譜元格式的插值多項式階次從M=2～5逐步升高，其一階MTF穩定臨界值逐步降低.但總體而言，MTF譜元格式的穩定臨界值均高于其有限元格式.將表2數值代入式(20)，結果如圖11所示.(根據文獻[34]給出的內域時間積分穩定條件確定?τ取值范圍為?τ≤0.86).

圖11 MTF譜元格式的時域穩定條件Fig.11 Time-domain stability condition of the spectral element scheme of MTF

圖11結果表明，不同MTF譜元格式的穩定性，隨著插值多項式階次的升高而逐步降低，變化趨勢與精度結果正好相反.對此可作如下解釋：當插值多項式階次升高(即越來越接近譜單元階次)時，相應的MTF譜元格式與內域的單元位移模式匹配得更好，表現為精度提高；但是高精度插值格式對人工波速的變化更為敏感，在人工波速增大時更有可能發生失穩.

5 結論

本文提出應用譜元法和透射邊界條件實現高效近場波動數值模擬的思路，探討了MTF與譜元離散模型結合的關鍵問題，得到的主要結論如下：

(1)MTF與譜元離散格式的結合，采用空間內插方案比較合適，不宜采用時間內插方案.

(2)傳統的MTF有限元格式適用于等距節點情形，不能直接套用到譜元法中.其中，傳統格式實現一階MTF的三點拋物線內插法可以借鑒，但需改變插值系數.高階MTF的齊次內插方法不能使用，本文提出的簡單內插方法可以作為一種替代方案.

(3)根據插值多項式階次的不同，可得到不同的MTF譜元格式，插值階次越接近譜單元階次，其精度越高.理論分析和數值試驗都表明，基于譜單元位移模式插值的MTF譜元格式具有最高的精度，原因是它不僅插值階次與單元階次一致，而且插值基函數也與單元形函數一致.

(4)不同MTF譜元格式的穩定臨界值隨著插值多項式階次的提高而降低，表明插值多項式階次較高的MTF譜元格式對人工波速的敏感性更強，在較大人工波速時相對較容易失穩.不過，所有格式均在人工波速大大超過介質物理波速時才可能發生失穩.

1 廖振鵬.工程波動理論導論.第二版.北京:科學出版社,2002 (Liao Zhenpeng.Introduction to Wave Motion Theories in Engineering(second edition).Beijing:Science Press,2002(in Chinese))

2 廖振鵬.近場波動的數值模擬.力學進展,1997,27(2):193-216 (Liao Zhenpeng.Numerical simulation of near-fiel wave motion.Advances in Mechanics,1997,27(2):193-216(in Chinese))

3 Chen YS,Yang GW,Ma X,et al.A time-space domain stereo finit di ff erence method for 3D scalar wave propagation.Computers&Geosciences,2016,96:218-235

4 Liu SL,Li XF,Wang WS,et al.A mixed-grid finit element method with PML absorbing boundary conditions for seismic wave modelling.Journal of Geophysics&Engineering,2014,11(5):1-13

5 曹丹平,周建科,印興耀.三角網格有限元法波動模擬的數值頻散及穩定性研究.地球物理學報,2015,58(5):1717-1730(Cao Danping,Zhou Jianke,Yin Xingyao.The study for numerical dispersion and stability of wave motion with triangle-based finit element algorithm.Chinese Journal of Geophysics,2015,58(5):1717-1730(in Chinese))

6 CanutoC,HussainiMY,QuarteroniA,etal.SpectralMethods:Fundamentals in Single Domains.Berlin:Springer,2006

7 Patera AT.A spectral element method for flui dynamics:laminar fl w in a channel expansion.Journal of Computational Physics, 1984,54(3):468-488

8 Seriani G,Su C.Wave propagation modeling in highly heterogeneous media by a ploy-grid Chebyshev spectral element method.Journal of Computational Physics,2012,20(2):345-351

9 林燈,崔濤,冷偉等.一種求解地震波方程的高效并行譜元格式.計算機研究與發展,2016,53(5):1147-1155(Lin Deng,Cui Tao, Leng Wei,et al.An efficient parallel spectral element scheme for solving seismic wave equation.Journal of Computer Research and Development,2016,53(5):1147-1155(in Chinese))

10 Priolo E,Seriani G.A numerical investigation of Chebyshev spectral element method for acoustic wave propagation//Proceedings of the 13th IMACS Conference on Comparative Applied Mathematics, Dublin,Ireland,1991,2:551-556

11 Seriani G,Priolo E.Spectral element method for acoustic wave simulation in heterogeneous media.Finite Elements in Analysis and Design,1994,16(3):337-348

12 Seriani G.A parallel spectral element method for acoustic wave modeling.Journal of Computational Acoustics,1997,5(1):53-69

13 Komatitsch D,Vilotte JP.The spectral element method:an efficient tool to simulate the seismic response of 2D and 3D geological structures.Bulletin of the Seismological Society of America,1998,88(2): 368-392

14 Komatitsch D,Liu QY,TrompJ,et al.Simulations of ground motion in the Los Angeles basin based upon the spectral-element method.Bulletin of the Seismological Society of America,2004,94(1):187-206

15 Komatitsch D,Tsuboi S,Tromp J.The spectral-element method in seismology.Seismic Earth:Array Analysis of Broadband Seismograms.American Geophysical Union,2005:205-227

16 王秀明,Seriani G,林偉軍.利用譜元法計算彈性波場的若干理論問題.中國科學:G輯,2007,37(1):1-19(Wang Xiuming,Seriani G,Lin Weijun.Several theoretic aspects for the calculation of elastic wave fiel using spectral element method.Science China(G Series),2007,37(1):1-19(in Chinese))

17 嚴珍珍,張懷,楊長春等.汶川大地震地震波傳播的譜元法數值模擬研究.中國科學:D輯,2009,39(4):393-402(Yan Zhenzhen, Zhang Huai,Yang Changchun,et al.Numerical investigation of seismic wave propagation of Wenchuan Earthquake using spectral element method.Science China(D Series),2009,39(4):393-402(in Chinese))

18 李洪建,韓立國,鞏向博.復雜構造網格化及高精度地震波場譜元法數值模擬.石油物探,2014,53(4):375-383(Li Hongjian,Han Liguo,Gong Xiangbo.High precision spectral element method based on grid discretization of complicated structure for seismic wavefiel numerical simulation.Geophysical Prospecting for Petroleum,2014,53(4):375-383(in Chinese))

19 李琳,劉韜,胡天躍.三角譜元法及其在地震正演模擬中的應用.地球物理學報,2014,57(4):1224-1234(Li Lin,Liu Tao,Hu Tianyue.Spectral element method with triangular mesh and its application in seismic modeling.Chinese Journal of Geophysics, 2014,57(4):1224-1234(in Chinese))

20 李孝波.基于譜元法的玉田震害異常研究.[博士論文].哈爾濱:中國地震局工程力學研究所,2014(Li Xiaobo.The investigation of seismic damage anomaly in Yutian based on spectral element method.[PhD Thesis].Harbin:Institute of Engineering Mechanics,Chinese Earthquake Administration,2014(in Chinese))

21 李孝波,薄景山,齊文浩等.地震動模擬中的譜元法.地球物理學進展,2014,29(5):2029-2039(Li Xiaobo,Bo Jingshan,Qi Wenhao, et al.Spectral element method in seismic ground motion simulation.Progress in Geophysics,2014,29(5):2029-2039(in Chinese))

22 Liu YS,Teng JW,Lan HQ,et al.A comparative study of finit element and spectral element methods in seismic wavefiel modeling.Geophysics,2014,79(2):T91-T104

23 He CH,Wang JT,Zhang CH.Nonlinear spectral-element method for 3D seismic-wave propagation.Bulletin of the Seismological Society of America,2016,106(3):1074-1087

24 Lysmer J,Kuhlemeyer R L.Finite dynamic model for infinit media.Journal of the Engineering Mechanics Division,1969,95(4): 859-878

25 Clayton R,Engquist B.Absorbing boundary conditions for acoustic and elastic wave equations.Bulletin of the Seismological Society of America,1977,67(6):1529-1540

26 Berenger JP.A perfectly matched layer for the absorption of electromagnetic waves.Journal of Computational Physics,1994,114(2): 185-200

27 KomatitschD,TrompJ.Aperfectlymatchedlayerabsorbingboundary condition for the second-order seismic wave equation.Geophysical Journal International,2003,154(1):146-153

28 廉西猛,單聯瑜,隋志強.地震正演數值模擬完全匹配層吸收邊界條件研究綜述.地球物理學進展,2015,30(4):1725-1733(Lian Ximeng,Shan Lianyu,Sui Zhiqiang.An overview of research on perfectly matched layers absorbing boundary condition of seismic forward numerical simulation.Progress in Geophysics,2015,30(4): 1725-1733(in Chinese))

29 Xie ZN,Matzen R,Cristini P,et al.A perfectly matched layer for fluid-soli problems:Application to ocean-acoustics simulations with solid ocean bottoms.Journal of the Acoustical Society of America,2016,140(1):165-175

30 He Y,Min MS,Nicholls D P.A spectral element method with transparent boundary condition for periodic layered media scattering.Journal of Scientifi Computing,2016,68(2):772-802

31 廖振鵬,黃孔亮,楊柏坡等.暫態波透射邊界.中國科學:A輯, 1984,6:556-564(Liao Zhenpeng,Huang Kongliang,Yang Baipo, et al.A transmitting boundary for transient wave.Scientia Sinica (Series A),1984,6:556-564(in Chinese))

32 Liao ZP,Wong HL.A transmitting boundary for the numerical simulation of elastic wave propagation.International Journal of Soil Dynamics and Earthquake Engineering,1984,3(4):174-183

33 戴志軍,李小軍,侯春林.譜元法與透射邊界的配合使用及其穩定性研究.工程力學,2015,32(11):40-50(Dai Zhijun,Li Xiaojun,Hou Chunlin.A combination usage of transmitting formula and spectral element method and the study of its stability.Engineering Mechanics,2015,32(11):40-50(in Chinese))

34 Cohen GC.Higher-order Numerical Methods for Transient Wave Equations.Berlin:Springer,2002

35 關慧敏,廖振鵬.局部透射邊界和疊加邊界的精度分析與比較.力學學報,1994,26(3):303-311(Guan Huimin,Liao Zhenpeng.A comparison between the discrete local transmitting boundary and the superposition boundary in wave propagation.Acta Mechanica Sinica,1994,26(3):303-311(in Chinese))

36 景立平,吳照營,鄒經相.近場波動數值模擬穩定性問題分析.地震工程與工程振動,2002,22(2):17-21(Jing Liping,Wu Zhaoying,Zou Jingxiang.Stability analysis of numerical simulation ofnear-fiel wave motion.Earthquake Engineering and Engineering Vibration,2002,22(2):17-21(in Chinese))

37 景立平.多次透射公式實用形式穩定性分析.地震工程與工程振動,2004,24(4):20-24(Jing Liping.Stability analysis of practical formula for multi-transmitting boundary.Earthquake Engineering and Engineering Vibration,2004,24(4):20-24(in Chinese))

38 Liao ZP,Liu JB.Numerical instabilities of a local transmitting boundary.Earthquake Engineering and Structural Dynamics,1992, 21(1):65-77

39 關慧敏,廖振鵬.局部人工邊界穩定性的一種分析方法.力學學報,1996,28(3):376-380(Guan Huimin,Liao Zhenpeng.A method for the stability analysis of local artificia boundaries.Acta Mechanica Sinica,1996,28(3):376-380(in Chinese))

40 關慧敏,廖振鵬.時域局部人工邊界的穩定性分析方法概述.世界地震工程,1997,13(2):1-7(Guan Huimin,Liao Zhenpeng.A summary on the methods for the stability analysis of artificia local boundaries in time domain.World Information on Earthquake Engineering,1997,13(2):1-7(in Chinese))

41 廖振鵬,周正華,張艷紅.波動數值模擬中透射邊界的穩定實現.地球物理學報,2002,45(4):533-545(Liao Zhenpeng,Zhou Zhenghua,Zhang Yanhong.Stable implementation of transmitting boundary in numerical simulation of wave motion.Chinese Journal of Geophysics,2002,45(4):533-545(in Chinese))

42 楊宇,李小軍,賀秋梅等.散射問題中消除多次透射邊界高頻振蕩失穩措施比較分析.地震工程學報,2014,36(3):476-481(Yang Yu,LiXiaojun,HeQiumei,etal.Comparisonofmeasuresforeliminatinghigh-frequencyinstability of amulti-transmittingboundaryin scattering problems.China Earthquake Engineering Journal,2014, 36(3):476-481(in Chinese))

43 謝志南,廖振鵬.透射邊界高頻失穩機理及其消除方法——SH波動.力學學報,2012,44(4):745-752(Xie Zhinan,Liao Zhenpeng.Mechanism of high frequency instability caused by transmitting boundary and method of its elimination——SH wave.Chinese Journal of Theoretical and Applied Mechanics,2012,44(4): 745-752(in Chinese))

44 Zhang L,Yu TB.A method of improving the stability of Liao’s higher-order absorbing boundary condition.Progress in Electromagnetics Research M,2012,27:167-178

45 章旭斌,廖振鵬,謝志南.透射邊界高頻耦合失穩機理及穩定實現——SH波動.地球物理學報,2015,58(10):3639-3648(Zhang Xubin,Liao Zhenpeng,Xie Zhinan.Mechanism of high frequency coupling instability and stable implementation for transmitting boundary——SH wave motion.Chinese Journal of Geophysics,2015, 58(10):3639-3648(in Chinese))

46 謝志南,廖振鵬.人工邊界高頻振蕩失穩機理的一點注記.地震學報,2008,30(3):302-306(Xie Zhinan,Liao Zhenpeng.A note for the mechanism of high frequency instability induced by absorbing boundary conditions.Acta Seismologica Sinica,2008,30(3):302-306(in Chinese))

IMPLEMENTATION OF MULTI-TRANSMITTING BOUNDARY CONDITION FOR WAVE MOTION SIMULATION BY SPECTRAL ELEMENT METHOD: ONE DIMENSION CASE1)

Xing Haojie Li Hongjing2)
(College of Civil Engineering,Nanjing Tech University,Nanjing211816,China)

Multi-transmitting formula(MTF)is considered to be a universal local artificia boundary condition,which is generally employed in finit element simulation of near-fiel wave motion.Due to the great di ff erence between spectral element method(SEM)and finit element method(FEM),the traditional numerical scheme of MTF cannot be simply adopted in SEM without any change.In order to make use of the advantages of MTF,i.e.,clear physical mechanism and controllable accuracy,basic problems involved in the combination of MTF and SEM are discussed in this paper, then the feasibility of spatial or temporal interpolation schemes are investigated,respectively.From the view of spatial interpolation scheme,a set of numerical formulas of MTF based on Lagrange polynomial are proposed,where the higherorder MTF is implemented via a simple iteration process.The accuracy and stability of the above MTF schemes are examined by a standard 1-D spectral element model of wave motion.The numerical results show that all schemes havecomparable accuracy for 1st-and 2nd-order MTF,and the MTF scheme based on spectral element displacement mode is superior to others for 3rd-or 4th-order MTF.On the contrary,the stability threshold descends with the growth of interpolation polynomials’order of di ff erent MTF schemes,but instabilities only occur under the unusual condition that artificia wave speed is far beyond the physical wave speed.

numerical simulation of wave motion,spectral element method,multi-transmitting boundary,MTF,numerical stability,numerical accuracy

P315

A

10.6052/0459-1879-16-282

2016–10–11收稿，2016–12–23錄用,2016–12–27網絡版發表.

1)國家自然科學基金資助項目（51278245）.

2)李鴻晶，教授，主要研究方向：地震工程學.E-mail:hjing@njtech.edu.cn

邢浩潔,李鴻晶.透射邊界條件在波動譜元模擬中的實現：一維波動.力學學報,2017,49(2):367-379

Xing Haojie,Li Hongjing.Implementation of multi-transmitting boundary condition for wave motion simulation by spectral element method:one dimension case.Chinese Journal of Theoretical and Applied Mechanics,2017,49(2):367-379