透射邊界條件在波動(dòng)譜元模擬中的實(shí)現(xiàn)：一維波動(dòng)1)

邢浩潔 李鴻晶

(南京工業(yè)大學(xué)土木工程學(xué)院，南京 211816)

透射邊界條件在波動(dòng)譜元模擬中的實(shí)現(xiàn)：一維波動(dòng)1)

邢浩潔 李鴻晶2)

(南京工業(yè)大學(xué)土木工程學(xué)院，南京 211816)

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

波動(dòng)數(shù)值模擬，譜元法，多次透射邊界，MTF，數(shù)值穩(wěn)定性，數(shù)值精度

引言

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

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

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

1 基本問題

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

(1)一維條件下得到的多次透射邊界數(shù)值格式可以不加修改地直接用于二維、三維模型，因?yàn)楹髢烧邔?shí)現(xiàn)MTF的方式也是“一維化”的，即只涉及到每個(gè)邊界節(jié)點(diǎn)上一條指向內(nèi)域的離散網(wǎng)格線上的節(jié)點(diǎn).

(2)一維波動(dòng)模型摒除了二維、三維模型的諸多復(fù)雜性，有助于獲得對(duì)MTF離散格式的精度和穩(wěn)定性較為清晰的認(rèn)識(shí)，而一維波動(dòng)問題的求解方法和數(shù)值規(guī)律，可為二維、三維波動(dòng)的模擬奠定基礎(chǔ).

1.1 一維波動(dòng)譜元模型

對(duì)于一維無限或半無限長直桿中無阻尼波的傳播問題，在適當(dāng)?shù)奈恢媒厝∪斯み吔绲玫接邢抻?jì)算模型.模型中質(zhì)點(diǎn)的運(yùn)動(dòng)由如下波動(dòng)方程控制

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

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

式中，u為譜元離散模型的節(jié)點(diǎn)位移向量，Me,Ke,fe分別為單元的質(zhì)量矩陣、剛度矩陣和等效節(jié)點(diǎn)向量.

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

譜單元位移模式定義在標(biāo)準(zhǔn)的參考單元 Λ ∈[-1,1]上，單元內(nèi)任意位置的位移由下式表示

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

表1 LSEM在參考單元上的節(jié)點(diǎn)坐標(biāo)Table 1 Node coordinates of LSEM in reference element

由于GLL節(jié)點(diǎn)為勒讓德多項(xiàng)式的極值點(diǎn)，是根據(jù)對(duì)勒讓德多項(xiàng)式的數(shù)值求解確定的，因此該形函數(shù)也被稱為勒讓德基函數(shù).

各個(gè)物理單元的節(jié)點(diǎn)位移與相應(yīng)參考單元的節(jié)點(diǎn)位移ue(ξi)是一一對(duì)應(yīng)的，單元節(jié)點(diǎn)坐標(biāo)通過轉(zhuǎn)換關(guān)系來確定，為1號(hào)單元節(jié)點(diǎn)的坐標(biāo)，?xE為譜單元尺寸.在等效積分方程的實(shí)際計(jì)算中，通過雅可比行列式|J|=dx/dξ=?xE/2實(shí)現(xiàn)積分區(qū)間的變換.

高精度譜單元相對(duì)于有限單元的優(yōu)勢(shì)在于數(shù)值頻散低，能夠以較少節(jié)點(diǎn)模擬波場.對(duì)于感興趣的最短波長λmin=c/fmax的空間離散，有限元通常采用10個(gè)單元左右[1]，而四階譜單元?jiǎng)t只需采用一個(gè)單元[10]，如圖1所示.

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

1.2 多次透射公式

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

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

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

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

2 MTF譜元格式

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

2.1 時(shí)空內(nèi)插方案的選取

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

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

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

式中各個(gè) MTF計(jì)算點(diǎn)所在的時(shí)刻為p+1-sj/ (ca?t)，sj為第j個(gè)譜元節(jié)點(diǎn)與邊界節(jié)點(diǎn)的距離，顯然它們內(nèi)域節(jié)點(diǎn)對(duì)應(yīng)的時(shí)刻不一致.

圖3 譜元離散網(wǎng)格的MTF時(shí)間內(nèi)插方案Fig.3 Temporal interpolation scheme of MTF in the discrete grid of spectral element

目前關(guān)于MTF時(shí)間內(nèi)插方案的研究很少，僅在文獻(xiàn)[33]中見到，其模擬結(jié)果顯示時(shí)間域插值的精度和穩(wěn)定性不如空間域插值.我們對(duì)MTF時(shí)間內(nèi)插方案進(jìn)行了研究，從波動(dòng)數(shù)值模擬的具體實(shí)踐并結(jié)合對(duì)多次透射公式本質(zhì)的認(rèn)識(shí)，認(rèn)為譜元離散網(wǎng)格的MTF時(shí)間內(nèi)插方案存在以下問題：

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

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

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

上述問題表明，時(shí)間內(nèi)插方案不宜作為研究MTF譜元格式的出發(fā)點(diǎn).接下來，將從空間內(nèi)插方案出發(fā)，導(dǎo)出具體的MTF譜元格式.

2.2 插值格式的確定

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

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

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

該格式人工邊界節(jié)點(diǎn)0在p+1時(shí)刻的位移由下式計(jì)算

其中，上標(biāo)T表示向量的轉(zhuǎn)置.式(8)中的插值系數(shù)tj,0,tj,1,tj,2與MTF計(jì)算點(diǎn)ja(j=1,2,···,N)以及式(9)中的譜元網(wǎng)格點(diǎn)是一一對(duì)應(yīng)的，其數(shù)值由拉格朗日插值公式確定.

其中，s0(s0=0)，s1和s2為離散網(wǎng)格點(diǎn)與邊界節(jié)點(diǎn)的距離(圖4)，可根據(jù)譜單元尺寸和表1給出的GLL節(jié)點(diǎn)坐標(biāo)進(jìn)行計(jì)算.

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

上式對(duì)能夠?qū)崿F(xiàn)的MTF階次作了一定限制.根據(jù)內(nèi)域時(shí)間積分的穩(wěn)定條件以及譜元節(jié)點(diǎn)位置關(guān)系不難得出，當(dāng)人工波速等于物理波速時(shí)(實(shí)際計(jì)算時(shí)通常如此)，能夠?qū)崿F(xiàn)的MTF階次一般能夠達(dá)到3以上，基本滿足使用要求.但是，若人工波速取值較大且時(shí)間步長接近穩(wěn)定臨界值時(shí)，則需要驗(yàn)算MTF階次是否滿足上式要求.

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

為改善邊界精度，提高插值多項(xiàng)式階次，即利用更多的譜元節(jié)點(diǎn)進(jìn)行插值，得到新的MTF譜元格式.若采用的插值多項(xiàng)式階次為M，則稱之為基于M次多項(xiàng)式插值的MTF譜元格式，如圖5所示.

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

此時(shí)人工邊界節(jié)點(diǎn)0在p+1時(shí)刻的位移計(jì)算公式為

式(13)中的插值系數(shù)tj,0,tj,1,···,tj,M可根據(jù)拉格朗日插值多項(xiàng)式計(jì)算

離散網(wǎng)格點(diǎn)與邊界節(jié)點(diǎn)的距離sk(k=0,1,···,M)，根據(jù)譜單元尺寸和表1給出的GLL節(jié)點(diǎn)坐標(biāo)進(jìn)行計(jì)算.

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

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

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

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

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

3 精度

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

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

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

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

對(duì)于人工波速等于介質(zhì)物理波速的情形，即ca=c，四種MTF譜元格式計(jì)算得到的邊界節(jié)點(diǎn)位移時(shí)程如圖7所示.

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

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

圖7 人工邊界節(jié)點(diǎn)位移時(shí)程(ca=c)Fig.7 Displacement history of the artificia boundary node(ca=c)

對(duì)于人工波速大于介質(zhì)物理波速的情形，取ca=2c，四種MTF譜元格式計(jì)算的邊界節(jié)點(diǎn)位移時(shí)程如圖8所示.

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

圖8 人工邊界節(jié)點(diǎn)位移時(shí)程(ca=2c)Fig.8 Displacement history of the artificia boundary node(ca=2c)

圖9 人工邊界節(jié)點(diǎn)位移時(shí)程(ca=0.5c)Fig.9 Displacement history of the artificia boundary node(ca=0.5c)

對(duì)于人工波速小于介質(zhì)物理波速的情形，取ca=0.5c，四種MTF譜元格式計(jì)算得到的邊界節(jié)點(diǎn)位移時(shí)程如圖9所示.

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

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

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

4 穩(wěn)定性

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

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

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

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

時(shí)，MTF才可能出現(xiàn)高頻失穩(wěn).反過來說，當(dāng)人工波速小于或等于1.5倍空間步距與時(shí)間步距的比值時(shí)，MTF不會(huì)發(fā)生失穩(wěn).對(duì)應(yīng)于MTF的頻域穩(wěn)定條件(即反射系數(shù)的模大于1)，將式(18)稱為MTF的時(shí)域穩(wěn)定條件.

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

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

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

對(duì)于一維波動(dòng)譜元模型，波速比α的定義與前文相同，而無量綱時(shí)間步距則定義為?τ=c?t/s1，s1為譜單元端部兩個(gè)節(jié)點(diǎn)之間的距離.若用 γc表示MTF譜元格式的穩(wěn)定臨界值，則有相應(yīng)的MTF時(shí)域穩(wěn)定條件為

采用上一節(jié)的一維波動(dòng)算例進(jìn)行試算，計(jì)算時(shí)間設(shè)定為300 s，以觀察到明顯的高頻振蕩現(xiàn)象作為失穩(wěn)標(biāo)準(zhǔn).對(duì)于上一節(jié)的四種MTF譜元格式，考慮一階MTF，在計(jì)算過程中嘗試α和?τ的不同組合，發(fā)現(xiàn)總是當(dāng)α?τ超過一定值時(shí)，對(duì)應(yīng)的MTF才出現(xiàn)失穩(wěn).得到四種MTF譜元格式下，一階MTF的穩(wěn)定臨界值如表2所示.

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

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

圖11 MTF譜元格式的時(shí)域穩(wěn)定條件Fig.11 Time-domain stability condition of the spectral element scheme of MTF

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

5 結(jié)論

本文提出應(yīng)用譜元法和透射邊界條件實(shí)現(xiàn)高效近場波動(dòng)數(shù)值模擬的思路，探討了MTF與譜元離散模型結(jié)合的關(guān)鍵問題，得到的主要結(jié)論如下：

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

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

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

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

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

2 廖振鵬.近場波動(dòng)的數(shù)值模擬.力學(xué)進(jìn)展,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 曹丹平,周建科,印興耀.三角網(wǎng)格有限元法波動(dòng)模擬的數(shù)值頻散及穩(wěn)定性研究.地球物理學(xué)報(bào),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 林燈,崔濤,冷偉等.一種求解地震波方程的高效并行譜元格式.計(jì)算機(jī)研究與發(fā)展,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,林偉軍.利用譜元法計(jì)算彈性波場的若干理論問題.中國科學(xué):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 嚴(yán)珍珍,張懷,楊長春等.汶川大地震地震波傳播的譜元法數(shù)值模擬研究.中國科學(xué):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 李洪建,韓立國,鞏向博.復(fù)雜構(gòu)造網(wǎng)格化及高精度地震波場譜元法數(shù)值模擬.石油物探,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 李琳,劉韜,胡天躍.三角譜元法及其在地震正演模擬中的應(yīng)用.地球物理學(xué)報(bào),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 李孝波.基于譜元法的玉田震害異常研究.[博士論文].哈爾濱:中國地震局工程力學(xué)研究所,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 李孝波,薄景山,齊文浩等.地震動(dòng)模擬中的譜元法.地球物理學(xué)進(jìn)展,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 廉西猛,單聯(lián)瑜,隋志強(qiáng).地震正演數(shù)值模擬完全匹配層吸收邊界條件研究綜述.地球物理學(xué)進(jìn)展,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 廖振鵬,黃孔亮,楊柏坡等.暫態(tài)波透射邊界.中國科學(xué):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 戴志軍,李小軍,侯春林.譜元法與透射邊界的配合使用及其穩(wěn)定性研究.工程力學(xué),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 關(guān)慧敏,廖振鵬.局部透射邊界和疊加邊界的精度分析與比較.力學(xué)學(xué)報(bào),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 景立平,吳照營,鄒經(jīng)相.近場波動(dòng)數(shù)值模擬穩(wěn)定性問題分析.地震工程與工程振動(dòng),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 景立平.多次透射公式實(shí)用形式穩(wěn)定性分析.地震工程與工程振動(dòng),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 關(guān)慧敏,廖振鵬.局部人工邊界穩(wěn)定性的一種分析方法.力學(xué)學(xué)報(bào),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 關(guān)慧敏,廖振鵬.時(shí)域局部人工邊界的穩(wěn)定性分析方法概述.世界地震工程,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 廖振鵬,周正華,張艷紅.波動(dòng)數(shù)值模擬中透射邊界的穩(wěn)定實(shí)現(xiàn).地球物理學(xué)報(bào),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 楊宇,李小軍,賀秋梅等.散射問題中消除多次透射邊界高頻振蕩失穩(wěn)措施比較分析.地震工程學(xué)報(bào),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 謝志南,廖振鵬.透射邊界高頻失穩(wěn)機(jī)理及其消除方法——SH波動(dòng).力學(xué)學(xué)報(bào),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 章旭斌,廖振鵬,謝志南.透射邊界高頻耦合失穩(wěn)機(jī)理及穩(wěn)定實(shí)現(xiàn)——SH波動(dòng).地球物理學(xué)報(bào),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 謝志南,廖振鵬.人工邊界高頻振蕩失穩(wěn)機(jī)理的一點(diǎn)注記.地震學(xué)報(bào),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網(wǎng)絡(luò)版發(fā)表.

1)國家自然科學(xué)基金資助項(xiàng)目（51278245）.

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

邢浩潔,李鴻晶.透射邊界條件在波動(dòng)譜元模擬中的實(shí)現(xiàn)：一維波動(dòng).力學(xué)學(xué)報(bào),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