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

邢浩潔 李鴻晶

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

固體力學

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

邢浩潔 李鴻晶2)

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

將邢浩潔和李鴻晶提出的多次透射公式(multi-transmitting formula,MTF)的譜元格式應用于均勻介質中線彈性SH波動問題的譜元模擬.假定緊鄰人工邊界的一層譜單元為具有直線邊界的四邊形單元，以保證每個人工邊界節點都唯一對應一條指向內域的離散網格線.人工邊界節點在某時刻的位移由該離散網格線上的節點在前若干時刻的位移確定，按照MTF譜元格式進行計算.通過平面波以一定角度傳播的外源問題算例和點源脈沖自由擴散的內源問題算例，驗證了方法的可行性以及對實際復雜波動問題的適用性.通過不同類型初值問題算例，在時域內分析了插值多項式階次、人工波速和透射階次三個參數對反射誤差的影響.結果表明：插值多項式階次較高的格式會表現出更好的精度，但總體上對反射誤差的影響較小；人工波速對反射誤差具有顯著影響，當人工波速小于介質物理波速時反射誤差較大，而當人工波速等于或稍大于介質物理波速時反射誤差處于較低水平；透射階次對反射誤差具有決定性影響，表現在不失穩的情形下提高透射階次能夠迅速降低反射誤差，但內源問題從三階MTF開始出現飄移失穩，外源問題從二階MTF開始出現輕微的飄移失穩.

多次透射公式,譜元法,SH波動,MTF譜元格式,飄移失穩

引言

局部場地的地震波動數值模擬是地震學和地震工程等領域的重要課題，這類復雜介質中的大規模開放系統問題，對計算模型的空間離散以及人工邊界條件都有很高要求.譜元法[2](spectral element method,SEM)兼具了有限元法的靈活性和譜方法的高精度，具有高效處理復雜幾何或物理模型的能力，成為近年來備受關注的空間離散方法.從20世紀90年代開始，譜元法被用于求解地震波傳播問題，在地球介質中的聲波和彈性波模擬方面取得了進展[310].近十年來對譜元法的研究逐漸增多，在譜元法基本理論[11]、大規模地震動模擬[1214]、三角網格劃分[1518]、并行計算[19]、非線性分析[20]以及與有限元法的對比[21]等方面都有成果發表.可以說，在處理大規模復雜波動問題方面，譜元法比傳統有限差分法[22]或有限元法[23]具有顯著優勢.但由于譜元法在有限單元上采用了基于不等距節點的高階插值格式，傳統的以有限元法等采用低階插值格式構造的人工邊界條件計算公式在波動譜元模型中無法被直接使用，因此需要專門研究適合譜元模型的人工邊界條件的具體實現方式.

Clayton-Engquist(CE)邊界[24]是早期在譜元法地震波動模擬中使用最多的人工邊界條件，它與譜元法結合的方式有兩種，一種需將邊界表達式寫成沿邊界的等效積分并納入到譜元法的等效積分方程中[4]，另一種則將邊界表達式寫成關于邊界上黏性阻尼力以及法向、切向應力的形式[79].作為一種時、空局部化的近似邊界，它僅能吸收傳播方向與邊界法線夾角較小的外行波動.全局化的DtN邊界在譜元法中有一定應用[2527]，該邊界由波傳播的遠場解得到，是一種精確的人工邊界.但全局邊界的計算量很大，許多實際波動問題的遠場解比較復雜，而在數值計算中需要對級數表達式進行截斷[28]，這些都限制了該邊界的推廣.完美匹配層(perfectly matched layer,PML)邊界是目前波動譜元模擬中比較熱門的人工邊界[2933]，其優點在于外行波在內域和吸收層的交界面上不發生反射，而且能夠迅速地被吸收層內的衰減函數和阻尼作用消減掉.但無論是經典PML，還是復頻移PML，形式都比較復雜，而且數值離散后的PML仍然存在精度和穩定性問題.最近，多次透射公式(multi-transmitting formula,MTF)[3435]被嘗試引入波動譜元模擬中，文獻[36]從時間回退和空間回退兩個角度，給出了在譜元法中使用多次透射公式的方法并探討了其穩定性問題.

MTF通過透射階次和人工波速兩個參數控制透射精度，譜元法則通過單元階次和單元尺度來模擬復雜波場，兼具高精度和靈活性是二者的相似之處.因此，在波動譜元模擬中使用MTF作為人工邊界條件，可能是提高波動模擬精度和計算效率的重要途徑之一.MTF易于與內域離散格式結合，但關于這方面的研究不多.目前在波動有限差分或有限元模擬中廣泛使用的是一種基于三點拋物線內插和齊次內插遞推過程的MTF數值格式[37].在文獻[1]中，我們詳細探討了在譜元離散網格中施加MTF時所面臨的新問題，提出了一組基于拉格朗日多項式插值和簡單內插方案的MTF譜元格式，并在一維波動問題中研究了這組格式的精度和穩定性特點.本文將根據文獻[1]提出的MTF譜元格式，進一步探討在二維波動譜元模擬中實現MTF的具體方法，分析方法的可行性以及復雜波動問題的物理特征與MTF計算參數之間的聯系，并初步探討與之相關的穩定性問題.

1 二維波動譜元模型

以SH型波動問題說明譜元法的基本公式以及MTF的實現過程.SH波動方程為如下的二維標量波方程

其中,u(x,y,t)為出平面波動位移，f(x,y,t)為外荷載，c為介質物理波速.模型所在區域為(x,y)∈Ω，計算時間t∈(0,T].

采用具有集中質量矩陣的高精度勒讓德譜元法(Legendre spectralelementmethod,LSEM)對波動方程進行空間離散，通過以下3個步驟：(1)導出原方程的等效積分“弱”形式；(2)將整個計算區域的積分劃分為各個單元積分的疊加；(3)將各個單元的積分和微分轉換到標準的參考單元上.得到波動方程的等價形式如下

其中，ue(ξ,η),ve(ξ,η)分別為參考單元上的波動位移函數和測試函數，頂標“¨”表示對時間的二階偏導數，上標T表示向量或矩陣的轉置.ne為模型的單元數；Λ為參考單元，ξ,η為參考單元坐標；|J|和J-1為雅可比行列式以及雅可比矩陣的逆陣.

參考單元上的波動位移函數，即LSEM的單元位移模式由下式定義

式中，單元的總體節點編號i以及ξ,η方向的局部節點編號j,k如圖1所示.單元節點數ng=nξ×nη，這里nξ,nη分別為ξ,η方向的節點數.

圖1 參考單元的節點編號Fig.1 Nodenumbering of a referenceelement

勒讓德基函數為定義在GLL(Gauss Lobatto Legendre)節點上的拉格朗日插值函數，具有表達式

式中，GLL節點ξi,ξj是數值分析和工程計算中的一類重要節點，很多文獻都給出了其數值，這里不再贅述.η方向的勒讓德基函數定義與上式相同.

方程(2)再經過如下步驟：利用伽遼金原理，選取與波動位移函數相同的測試函數；將波動位移函數、測試函數、雅可比行列式以及雅可比逆矩陣代入方程(2)，根據二維GLL數值積分計算各個積分項，得到各個單元的運動方程；將單元方程組裝為總體方程.得到空間離散后的譜元節點運動方程

式中，M為總體質量矩陣,K為總體剛度矩陣，F為外載荷向量，u為被求的節點位移.

這個運動方程是關于模型所有內域節點和自由邊界節點的(譜元法基本公式已將自由邊界條件包含在內)，再引入物理邊界條件和人工邊界條件，就能通過時域逐步積分對上述方程進行求解.其中，采用適當的人工邊界條件，降低或消除人工邊界反射波的干擾，是保證波動模擬精度的關鍵.

2 多次透射公式

有限模型人工邊界節點的位移由MTF進行計算.MTF是一種直接以離散形式給出的一維化的時空外推公式，它由下式定義

由于MTF計算點與內域離散點位置不一定重合，上式無法直接用于波動數值計算.只需通過適當的插值方法，由臨近內域離散點位移插值得到各個MTF計算點位移，就能得到適用于波動數值模擬的MTF數值格式.

2.1 人工邊界附近的譜單元類型

在波動有限元模擬中使用MTF時，通常要求每個邊界節點的內法線上有一系列等距分布的有限元節點，然后利用這些節點進行三點拋物線內插和齊次內插遞推方案，得到實用的MTF有限元格式[3435].其中包含兩個要點，一是邊界法線，二是等距分布的有限元節點.然而研究發現：“法線”的條件不必嚴格滿足，但它必須是一條將人工邊界節點和一組內域有限元節點串連起來的直線，這條直線通常是從邊界節點指向內域的離散網格線；等距分布的有限元節點是保證齊次內插過程得以實現的條件，因此在譜元不等距節點分布的條件下無法使用齊次內插方案.

為了在譜元離散網格中實現一維化的MTF，首先需要確保每個人工邊界節點都對應一條上述“直線”.如圖2所示的三角形譜單元和曲邊譜單元顯然不滿足這一要求，因為前者有部分邊界節點上沒有指向內域的離散網格線，而后者有大量的離散網格線為曲線.

如圖3所示的具有直線邊界的四邊形譜單元能夠滿足上述要求.同時由于譜單元尺寸較大，對于常用階次的MTF，其計算點一般不會超過一個譜單元范圍.于是對邊界附近的譜單元類型作如下假定：緊鄰人工邊界的一層譜單元為具有直線邊界的四邊形單元.

圖2 三角形譜單元和曲邊譜單元Fig.2 Triangleand curvilinear spectralelements

圖3 直邊譜單元中的MTF計算點Fig.3 Computation nodesofMTF in a spectralelementw ith straightedges

基于該假定，每個人工邊界節點的MTF計算點都位于一條指向內域的離散網格線上，該點位移可由文獻[1]給出的MTF譜元格式進行計算.

2.2 邊界位移計算公式

以圖3中A點為例來說明在二維譜元離散網格中施加MTF的具體過程，其他人工邊界節點與之類似.A點指向內域的離散網格線為一條直線，線上有不均勻分布的譜元節點A,A1,A2,···,A′，以及均勻分布的MTF計算點1a,2a,···,Na.各個MTF計算點分別位于不同的時刻，如圖4所示.

根據文獻[1]提出的MTF譜元格式，對各個MTF計算點的位移均采用同樣的拉格朗日多項式插值公式進行計算，得到A點在p+1時刻的位移計算公式為

圖4 計算A點位移的MTF譜元格式Fig.4 SpectralelementschemeofMTF for computing the displacement of node A

式中，上標T表示向量轉置，M為插值多項式階次.插值系數tj,0,tj,1,tj,2,···,tj,M與MTF計算點ja(j=1,2,···,N)和譜元網格點0,1,2,···,M一一對應.s0,s1,···,sM分別為節點A與A,A與A1,···,A與AM之間的距離.

插值多項式階次M的取值范圍為2,3,···,NEn；NEn為MTF計算方向的譜單元階次，因此式(8)～式(11)表示的是一系列MTF譜元格式.不同格式的區別在于采用的插值多項式階次M不同，即涉及的譜元離散點數目不同.其中，M=2時為三點拋物線內插格式，M=NEn時為與譜單元位移模式一致的格式.式(8)～式(11)表示的MTF譜元格式能夠實現的MTF階次范圍，由保證各個插值均為“內插”的條件確定，為

可以看出，插值多項式階次M越高，能夠實現的MTF階次N越高.對于最低的插值多項式階次M=2，通常N取為2和3，所以上式對常用的MTF階次一般不構成限制.只有在少數極端情形下，如人工波速大大超過物理波速或透射階次很高時，才需要驗證.

上述MTF譜元格式是在較為稀疏的不等距譜元節點上進行插值，用于計算各個MTF計算點1a,2a,···,Na位移的譜元節點數均為M+1個；而傳統MTF有限元格式是在相對密集的等距有限元節點上進行插值，用于計算各個MTF計算點1a,2a,···,Na位移的有限元節點數分別為3,5,7,···,2N+1個.二者的主要區別在于實現高階MTF的方式不同，以及插值多項式階次不同.MTF有限元格式實現高階MTF的方式為齊次內插遞推過程，這樣的好處是反射系數表達式為R=-fN，f為一階MTF的反射因子.此時各階MTF的穩定界限相同，且在不失穩的情況下，隨著透射階次N的增加，反射系數能夠穩步地減小.MTF譜元格式使用簡單內插方案實現高階MTF，其模擬效果比齊次內插方案有所降低，但是譜單元的高精度特性能夠較好地彌補這方面不足，而高階譜單元下插值多項式階次M的取值則是MTF譜元格式面臨的新問題，對此文獻[1]已有初步討論，本文還將進一步研究.

需要指出的是，上述MTF譜元格式和傳統MTF有限元格式都來源于MTF定義式(7)，因此它們的模擬效果首先是由人工波速ca和透射階次N這兩個基本參數決定的.

3 方法可行性的驗證

分別以無限均勻介質中的外源和內源SH波動問題為算例，驗證本文方法的可行性.計算模型如圖5所示，模型所在區域為寬600m、高400m的矩形，介質物理波速c=400m/s，模型四邊均為MTF人工邊界.外源問題的輸入波為從模型左下方斜入射的平面波，斜入射方向與豎軸夾角為α.內源問題的輸入波為作用在模型內部某點的一個脈沖波源，波源位置距模型左側、下側邊界各100m.

輸入波位移時程采用Ricker型地震子波，具有表達式

式中，f為中心頻率，t0為脈沖波形的起始時間.由此不難得出脈沖寬度為2/f，峰值頻率約為3f.計算時取f=10Hz，對應峰值頻率約為30Hz.對于內源問題，上式為點波源的輸入波位移時程；對于外源問題，上式僅為模型左下角點的輸入波位移時程，模型左側、下側邊界上其他節點的輸入波位移時程，需在上式的基礎上考慮相應的滯后時間.

圖5 外源和內源SH波動問題計算模型Fig.5 Computationmodelof SH-typewavemotion problem w ith outer or innerwave source

采用四階勒讓德譜單元對模型進行空間離散，單元節點數為25，模型的單元數為1247(43×29=1247).兩個方向的譜單元尺寸分別為13.95m和13.79m，接近于λmin=400/30≈13.33m，符合譜元離散的精度要求.時間步長取為0.003s，此時，與時域積分穩定性有關的庫朗數(Courantnumber)為q=c?t/?xmin=400×0.003/2.38≈0.50，小于臨界值(略小于1的某個常數)，符合穩定性要求.

3.1 外源問題

首先進行外源問題的模擬.平面波入射角度設置為α=30°，此時外行波在上邊界和右邊界的透射角度(相對于邊界法線)分別為θ=30°和θ=60°.邊界條件采用式(8)～式(11)表示的MTF譜元格式，插值多項式階次取為M=4(與譜單元階次一致).采用不同邊界參數進行數值模擬，選取其中比較有代表性的兩組結果，它們在0.3s,0.6s,1.1s,1.4s的波場快照如圖6所示.

圖6(a)為第一組模擬結果，使用的邊界參數為：二階MTF，人工波速ca取為介質物理波速c，這是MTF最常用的參數設置.前兩幅子圖顯示，在模型的左側和底部順利地實現了平面波的輸入，此時沒有外行波，邊界上不會發生反射.后兩幅子圖顯示，當入射波到達右邊界和上邊界時，MTF邊界能夠將絕大部分外行波動能量透射出去，僅有很小一部分能量被反射回計算區域.從圖中還可以看出，右邊界的反射明顯大于上邊界的反射，其原因在于右邊界為大角度透射(θ=60°)，而上邊界的透射角度相對要小得多(θ=30°).這組模擬表明了本文提出的在二維譜元離散網格中使用MTF作為人工邊界條件的方法是可行的.對于透射角度不大的外行波，常用的MTF參數設置就能取得較好的模擬效果；對于透射角度較大的外行波，模擬效果有所降低.

圖6外源問題波場快照Fig.6 Wave fiel snapshotsof outer-source problem

圖6 (b)為第二組模擬結果，使用的邊界參數為：一階MTF，人工波速ca取為各邊的法向透射速度c/cosθ，即左、右邊界ca=2c，上、下邊界前兩幅子圖顯示平面波能夠在模型的左側和底部順利地輸入計算區域，與上一組相同.后兩幅子圖則顯示當入射波到達模型右邊界和上邊界時，能夠完全透射出去，不發生任何反射.這種近乎完美的模擬結果表明,通過調整人工波速取值，可以很好地克服僅由透射角度帶來的影響.這組模擬不僅驗證了本文方法的可行性，還進一步顯示了它在處理簡單波動問題時的精確性.

平面波傳播問題是一種理想化的波動模型，波動能量沿單一方向傳播的特性，使得每個人工邊界上外行波的視傳播速度都可以確定，這為設置人工邊界帶來了方便.實際波動問題則要復雜得多，其復雜性表現為由透射角度或波動能量空間擴散的衰減效應導致的在人工邊界不同位置處，外行波的視傳播速度各不相同.對于這種復雜的實際波動問題，本文MTF譜元格式同樣能夠很好地模擬.

3.2 內源問題

接下來進行內源問題的模擬.為顯示本文方法適應復雜波場的能力，將波源設置在模型左下角以增加邊界附近波場的復雜性.此時對于左側和下側邊界，在靠近波源的位置會受到較為明顯的波場幾何衰減效應的影響，遠離波源的位置則存在外行波透射角度較大的問題；上側和右側邊界受到的影響相對較小.采用式(8)～式(11)給出的MTF譜元格式，插值多項式階次取為M=4.由兩組不同邊界參數得到的0.3s,0.6s,0.9s,1.5s的波場快照如圖7所示.

圖7內源問題波場快照Fig.7 Wave fiel snapshotsof inner-source problem

圖7 (a)使用的是MTF常用的邊界參數：二階MTF，ca=c，它對弱衰減、小角度透射的外行波模擬效果較好.圖中結果顯示在人工邊界各個位置處，大部分波動能量都能夠順利透出，表明了本文MTF譜元格式模擬復雜波場問題的有效性.僅有少量反射出現在左側、下側邊界距波源較近的波場衰減區和遠端的大角度透射區.如：0.6s子圖顯示了由于波場衰減效應導致的未被完全透出的波；0.9s子圖中A點和B點反射波強度對比，表明本例的邊界參數設置對大角度外行波動的模擬效果不如小角度外行波動；1.5s子圖中下邊界的反射波強度大于上邊界，因為下邊界的透射角度更大.

為了保證水利工程的順利招標，我國出臺了各種法律法規，以此來確保招標工作開展的科學性和合理性。但是，在實際的招標過程中，仍然存在著不按照國家相應的法律法規進行招投標的現象。尤其是有些小型企業在招標的過程中進行相應的暗箱操作，這樣就在一定程度上打破了整個行業的競爭平衡。而有些企業則是依靠人脈、金錢等進行非法的招標，這樣不僅給公平、公正以及公開的招標環境造成了影響，而且還會影響招標工作的順利展開。

圖7(b)使用的邊界參數為：二階MTF，ca=c/cos45°，此時它對透射角度為45°左右的外行波模擬效果最好.由于本例中外行波在大部分人工邊界區域都是以一定角度透射的，所以不難預測采用45°角為透射中心，能夠顯著改善模型4個角點附近，以及左側、下側邊界上遠離波源部分的透射效果.圖中結果顯示這組邊界參數的模擬效果總體上要優于上組邊界參數.如：0.9s子圖中B點的法向反射波依稀可見，而A點附近接近45°角的反射波則幾乎被完全吸收；1.5s子圖中C點附近區域與上一組模擬結果相比，小角度反射波幾乎看不出差別，但是大角度反射波被吸收得更加充分.

內源問題算例不僅再次證明了本文方法的可行性，而且表明了它對復雜波動問題的適應性，這種適應性主要通過選取適當的邊界參數來實現.就工程應用而言，常用的MTF參數配置(N=2，ca=c)對很多實際波動問題已具有足夠的模擬精度.然而從研究角度看，深入研究MTF計算參數與波動物理特征之間的聯系，根據外行波的視傳播規律有針對性地設定MTF計算參數，對更好地求解復雜波動問題具有重要價值.

4 邊界反射誤差的時域分析

式(8)～式(11)描述的MTF譜元格式對外行波的透射效果受到插值多項式階次M、人工波速ca和透射階次N三個參數的影響，通過一組精心設計的SH波動數值試驗，在時域內定量地研究不同MTF參數取值下，本文MTF譜元格式在人工邊界上的反射誤差.

如圖8所示，采用兩個計算模型：半空間模型Ω2，左側使用式(8)～式(11)的MTF譜元格式，其余三側為自由邊界；全空間模型Ω3，四邊均為自由邊界.反射誤差的分析區域為Ω1，這樣半空間模型的計算結果為受到左側人工邊界影響的數值解，全空間模型計算結果為不受任何邊界影響的大區域數值解(即精確解).介質物理波速為1m/s，計算時間為2s，這組參數能夠保證計算時間內自由邊界的反射波不會返回分析區域.

圖8 人工邊界反射誤差分析模型Fig.8 Numericalmodel for analyzing the reflectio errorof artificia boundary

為便于分析，本文模擬初始波場的自由擴散問題.由于沒有外部波動能量輸入以及介質阻尼等因素的干擾，模型中波動能量的變化只受邊界影響.初始波場所在區域為圖中圓形區域：圓心為(0.5,0)，半徑為0.45m.分別采用三種具有不同傳播特征的初始波場：圓形擴散(circularspread)、法向透射(normal transmission)和45°角透射(45°transmission)，其中圓形擴散初始波場的表達式為如下高斯型脈沖

其中，r2=(x-0.5)2+y2.法向透射的初始波場由在上式右端乘以cos(8πx)得到，45°角透射的初始波場則通過在上式右端乘以cos(8πx-8πy)得到.根據公式λmin=1/kmax來確定3種初始波場中最短波長成分，其中kmax為波數的上限值.由于初始波場是二維的，為簡化分析，取波場變化最劇烈的方向進行一維傅里葉分析，圓形擴散情形計算exp（-30r2）的傅里葉譜，后兩種情形計算exp（-30r2）·cos(8πr)的傅里葉譜.得到圓形擴散情形kmax的值約為4m-1，后兩種情形約為8m-1.3種情形的初始速度都為零.

x和y方向的譜單元階次均取為5，即采用36節點的單元.對于圓形擴散初始波場，單元尺寸取為1/4，對于法向透射和45°角透射初始波場，單元尺寸取為1/7.區域Ω1內的誤差波場定義為由半空間模型Ω2計算的數值解減去由全空間模型Ω3計算的精確解.衡量各個時刻MTF邊界反射誤差的指標為：區域Ω1內誤差波場的弗羅貝尼烏斯范數與初始波場的弗羅貝尼烏斯范數的比值.其中，各個時刻波動位移矩陣U的弗羅貝尼烏斯范數的計算公式為

4.1 插值多項式階次的影響

3種初始波場情形下，MTF邊界上的反射誤差隨插值多項式階次M的變化規律如圖9所示(MTF計算參數：N=2，ca=c，M=2,3,4,5).

圖9 插值多項式階次M對反射誤差的影響Fig.9 Reflectio errorofMTFa ff ected by interpolation polynomial order M

圖9插值多項式階次M對反射誤差的影響(續)Fig.9 Reflectio errorofMTFa ff ected by interpolation polynom ial order M(continued)

圖9 結果顯示，人工邊界反射誤差隨著插值多項式階次M增加而逐步降低，其中，當M=2時反射誤差最大，3種情形的峰值誤差均接近4%;當M=5時反射誤差最小，3種情形下的峰值誤差分別為3.5%,0.9%和2.7%.顯然，當插值多項式階次與譜單元階次一致時，模擬效果最好.另外圖中結果還有兩點值得重視：一是所有反射誤差均在5%以內，這說明如果從工程精度考慮，由插值多項式階次不同引起的誤差變化幾乎可以忽略不計；二是不同插值多項式階次對透射方向比較單一的波動問題，如法向透射或45°角透射，具有層次較為分明的影響，而對透射方向比較豐富的波動問題，如圓形擴散，其影響程度較小.

總之，插值多項式階次對透射效果的影響是在較高精度水平上的微小變化，它并不是決定式(8)～式(11)的MTF譜元格式模擬效果的關鍵因素.從實際使用角度出發，考慮到邊界附近的譜單元階次會受到特定波動問題的模擬精度和計算效率的影響而有所不同，若存在多種不同插值多項式階次的邊界格式可供選擇，會大大增加問題的復雜性，于是筆者建議統一采用插值多項式階次與譜單元階次一致的邊界格式.

4.2 人工波速的影響

3種初始波場情形下，MTF邊界上的反射誤差隨人工波速ca的變化規律如圖10所示(MTF計算參數：N=2，M=5，α=ca/c，取值為0.5,0.75,1.0,1.25,1.5).

圖10 人工波速ca對反射誤差的影響Fig.10 Refection errorofMTFa ff ected by artificia wave speed ca

對比圖9和圖10結果可知，不同人工波速ca取值導致的邊界反射誤差變化幅度，大大超過不同插值多項式階次M導致的變化.這點不難理解，因為人工波速是MTF定義式中的一個基本參數.較大反射誤差出現在人工波速取值小于物理波速的情形，如ca=0.5c時，峰值誤差分別達到10%,6.2%和13%，ca=0.75c時，峰值誤差分別為5.3%,1.6%和6%.而當人工波速取值等于或大于物理波速時，反射誤差均處于較低水平，其峰值分別不超過3.5%,1.7%和2.7%.反射誤差越小，表明人工波速取值越接近外行波沿邊界法向的視傳播速度.

波動問題的復雜性在于受到透射角度、波幅衰減或多個子波疊加等因素的影響，確定外行波的視傳播速度常常比較困難，而且很多時候難以用一個值來精確地描述.但這部分研究能夠幫助我們認識不同因素對視傳播速度的具體影響，進而總結一些有助于優化人工波速取值的結論.圖10(c)結果表明：外行波透射角度對視傳播速度的影響是較為確定的，以θ角透射的外行波的視傳播速度為c/cosθ.圖中模擬的是45°角透射波，其視傳播速度為α=1.5曲線與此最接近，因此反射誤差最小，其他α取值偏離越遠，反射誤差越大.圖10(b)模擬的是法向透射波，但反射誤差最小值沒有出現在α=1時，這表明了初始波峰附近波動能量幾何擴散導致的波幅衰減效應帶來的影響.從圖中結果不難看出，采用稍大于介質物理波速的人工波速能夠有效地弱化這一影響，而本例中最優人工波速取值在c和1.25c之間.圖10(a)結果顯示了單一人工波速的MTF譜元格式對實際復雜波動問題的模擬效果，它與不同時刻波動能量在邊界上的透射區域和透射方向有關.在1.0 s之前，法向和小角度透射波先到達邊界，此時人工波速等于c的模擬效果最好；在1.0s之后，大角度透射波成為主要成分，此時人工波速等于1.25c或1.5c的模擬效果要好于前者；而對于人工波速等于0.5c或0.75c的格式，其模擬效果自始至終都不理想.

以上分析表明，實際波動問題與理想的平面波法向透射問題相比，存在透射角度或波幅衰減，二者都會導致視傳播速度增大，因此采用大于介質物理波速的人工波速能夠達到更好的模擬效果.這一現象可理解為：從人工邊界節點的視角來看，存在透射角度或波幅衰減加速了波動能量向人工邊界的“積聚”過程.由此可以推斷對于本例未涉及的多個子波疊加效應，其作用與前兩者相似.

至此，可將均勻介質中標量波問題的人工波速取值規律總結為：不應取小于介質物理波速的值；對于大多數實際波動問題，取等于或稍大于介質物理波速的值就能滿足要求；只有在少數特殊情形下，如透射角度較大、邊界距離點波源很近或較多子波場在邊界附近疊加時，才可能需要采用大大超過介質物理波速的值，但這往往會帶來諸如穩定性等方面的問題.

4.3 透射階次的影響

3種初始波場情形下，MTF邊界上的反射誤差隨透射階次N的變化規律如圖11所示(MTF計算參數：M=5，ca=c，N=1,2,3).

圖11 透射階次N對反射誤差的影響Fig.11 Reflectio errora ff ected by transm itting order N

將圖11和圖10、圖9對比之后可以看出，透射階次N對邊界反射誤差的影響是決定性的，它大大超過了人工波速ca或插值多項式階次M的影響.其中，一階MTF的反射誤差較大，如圓形擴散和45°角透射情形，峰值誤差分別為8.7%和10.2%，不滿足工程精度要求，說明一階MTF難以用來模擬復雜波動問題.二階MTF的反射誤差顯著降低，三種情形下的峰值誤差分別為3.5%,0.9%和2.7%，都小于5%，滿足工程精度要求，說明通常采用二階MTF進行波動模擬的做法是可取的.三階MTF的反射誤差沒有像預想的那樣進一步降低，反而出現了不合理的增長，如圓形擴散和法向透射情形，三階MTF的反射誤差大大超過一階MTF，45°角透射情形，三階MTF的反射誤差與二階MTF相當，但前者隨時間增長的趨勢一直在延續.我們還采用四階MTF進行了模擬，其反射誤差已經大到無法在圖中顯示的地步，更高階次的MTF同樣如此.這種現象表明N≥3時MTF邊界發生了某種失穩.

觀察三階MTF模擬的各個時刻波場變化，我們發現這是一種飄移失穩[38]，如圖12所示.飄移現象從MTF邊界中部開始，逐漸向兩端擴展，形成中間大、兩端小的拱形位移，拱的高度隨時間不斷上升，而實際邊界位移在主波形透過之后應當趨近于零.顯然，當邊界出現這種嚴重的失穩問題時，波動模擬過程失敗.從實用角度看，不使用N≥3的MTF譜元格式就可以回避這個問題，因為這里二階MTF沒有出現失穩，而且通常具有足夠的模擬精度.關于本文MTF譜元格式的飄移失穩機理及其抑制措施，我們進行了初步研究，已得到一些基本認識.

圖12 區域Ω1在2s時的波場Fig.12 Wave fiel of domainΩ1at2s

飄移失穩現象不是本文MTF譜元格式所特有的，傳統MTF有限元格式可能出現飄移失穩已是眾所周知[3841].而其他類型的人工邊界，如文獻[42]提出的幾種Higdon邊界形式，同樣可能出現飄移失穩.關于飄移失穩機理目前還沒有統一的解釋，如：李小軍等[3839]認為飄移失穩是由于外源問題采用連續介質中的波動理論解作為輸入波場，而它與有限元離散解之間存在差異所造成的.周正華等[40-41]認為是由于MTF不滿足GKS準則導致數值解中的零頻和零波數成分可以通過邊界進入計算區，從而引起飄移失穩.Higdon[42]則將其解釋為邊界格式允許波動模擬中出現關于零頻成分的廣義特征值.我們認為這可能是由于人工邊界條件作為計算模型的支承條件本身就存在一定缺陷所致.

數值試驗結果表明，本文MTF譜元格式可能出現飄移失穩的透射階次，對于內源問題為N≥3，對于外源問題為N≥2.其中，前者與文獻[42]的結果一致，后者則與文獻[38]結論相同.由此不難看出，飄移失穩可能是邊界高階形式固有問題，而當邊界上存在外部波動輸入時，這一問題更為嚴重.已有研究提出了幾種抑制MTF飄移失穩的方法，如降低MTF階次[38]、在邊界區添加彈簧和阻尼元件[39]，以及在MTF公式中增加一個微小正數[40].我們認為這些方法都可以應用到本文的MTF譜元格式當中，這部分內容值得深入探討，但已超出本文的范疇.

最后，對本文MTF譜元格式的高頻振蕩失穩特性做一簡要說明.研究過程中我們進行了大量SH波動模擬，邊界均未出現高頻振蕩失穩.結合文獻[1]關于一維波動問題的研究，我們將其解釋為：本文MTF譜元格式的穩定條件比內域離散格式寬松，所以一般情況下不會發生高頻振蕩失穩，只有在人工波速大大超過介質物理波速或透射階次很高時，才有可能超出穩定界限，在邊界上出現高頻振蕩失穩.

5 結論

(1)一維化的MTF公式需要通過一條直線上的離散網格點來實現，本文通過假定緊鄰人工邊界的一層譜單元為直線邊界的四邊形單元，保證每個人工邊界節點都唯一地對應一條指向內域的離散網格線，用于實現文獻[1]所提出的MTF譜元格式.外源和內源SH波動問題算例證實了方法的可行性以及對復雜波動問題的適應性.

(2)插值多項式階次M對邊界反射誤差的影響較小，M取值較高的MTF譜元格式精度略好于M取值較低的格式.從使用方便的角度，我們建議M取值統一采用MTF計算方向的譜單元階次，這同時也是精度最好的選擇.

(3)人工波速ca取值對邊界反射誤差具有顯著影響.就本文研究的均勻介質中SH波動問題而言，外行波存在透射角度、波幅衰減或多個子波場疊加，都會導致波動能量更快地向邊界“積聚”，因此采用等于或稍大于介質物理波速的人工波速能得到較好的模擬效果.采用小于或等于介質物理波速的人工波速，推算邊界節點位移的準確性較差.其他類型波動問題，如涉及到成層場地、P-SV波動，或Rayleigh面波等，情況較為復雜，可在此基礎上進行專門討論.

(4)透射階次N對本文MTF譜元格式的模擬效果具有決定性影響.不發生失穩時，提高透射階次對降低反射誤差的效果最為明顯.但是，高階MTF容易出現飄移失穩，其中，內源問題和外源問題開始出現飄移失穩的MTF階次分別為3和2.

1邢浩潔,李鴻晶.透射邊界條件在波動譜元模擬中的實現：一維波動.力學學報,2017,49(2)：367-379(Xing Haojie,LiHongjing.Implementation ofmulti-transm itting boundary condition for thewave motion simulation by spectralelementmethod：onedimension case.Chinese JournalofTheoreticaland Applied Mechanics,2017,49(2)：367-379(in Chinese))

2 Patera AT.A spectral elementmethod for flui dynam ics：Lam inar fl w in a channelexpansion.JournalofComputational Physics,1984,54(3)：468-488

3 Priolo E,SerianiG.A numerical investigation of Chebyshev spectral elementmethod for acoustic wave propagation//Proceedings of the13th IMACSConferenceon Comparative Applied Mathematics,Dublin,Ireland,1991,2：551-556

4 SerianiG,Priolo E.Spectralelementmethod foracousticwavesimulation in heterogeneousmedia.Finite Elements in Analysis and Design,1994,16(3)：337-348

5 Seriani G.A parallel spectral elementmethod for acoustic wave modeling.JournalofComputationalAcoustics,1997,5(1)：53-69

6 SerianiG,Su Chang.Wave propagationmodeling in highly heterogeneousmedia by a ploy-grid Chebyshev spectralelementmethod.JournalofComputationalphysics,2012,20(2)：345-351

7 Komatitsch D,Vilotte JP.The spectralelementmethod：an e ffi cient tool to simulate theseism ic responseof 2D and 3D geologicalstructures.Bulletin ofthe SeismologicalSociety ofAmerica,1998,88(2)：368-392

8 Komatitsch D,Vilotte JP,VaiR,etal.The spectralelementmethod for elastic wave equations—application to 2-D and 3-D seism ic problems.International Journal for Numerical Methods in Engineering,1999,45(9)：1139-1164

9 Komatitsch D,Barnes C,Tromp J.Wave propagation near a flui solid interface：a spectral-element approach.Geophysics,2000,65(2)：623-631

10 Komatitsch D,Liu Qinya,Tromp J,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

11王秀明,SerianiG,林偉軍.利用譜元法計算彈性波場的若干理論問題.中國科學：G輯,2007,37(1)：1-19(Wang Xium ing,SerianiG,Lin Weijun.Several theoretic aspects for the calculation of elasticwave fiel using spectralelementmethod.Science China(G series),2007,37(1)：1-19(in Chinese))

12嚴珍珍,張懷,楊長春等.汶川大地震地震波傳播的譜元法數值模擬研究.中國科學：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 elementmethod.Science China(D series),2009,39(4)：393-402(in Chinese))

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

14李孝波,薄景山,齊文浩等.地震動模擬中的譜元法.地球物理學進展,2014,29(5)：2029-2039(LiXiaobo,Bo Jingshan,QiWenhao,etal.Spectralelementmethod in seism ic groundmotion simulation.Progress in Geophysics,2014,29(5)：2029-2039(in Chinese))

15劉有山,滕吉文,徐濤等.三角網格譜元法地震波場數值模擬.地球物理學進展,2014,29(4)：1715-1726(Liu Youshan,Teng Jiwen,Xu Tao,etal.Numericalmodeling of seism ic wave fiel w ith the SEM based on triangles.Progress in Geophysics,2014,29(4)：1715-1726(in Chinese))

16李琳,劉韜,胡天躍.三角譜元法及其在地震正演模擬中的應用.地球物理學報,2014,57(4)：1224-1234(Li Lin,Liu Tao,Hu Tianyue.Spectral elementmethod w ith triangularmesh and its application in seismicmodeling.Chinese Journal ofGeophysics,2014,57(4)：1224-1234(in Chinese))

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

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

19 He Chunhui,Wang Jinting,Zhang Chuhan.Nonlinear spectralelementmethod for 3D seismic-wave propagation.Bulletin of the Seismological Society ofAmerica,2016,106(3)：1074-1087

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

21 Liu Youshan,Teng Jiwen,Lan Haiqiang,etal.A comparativestudy of finit elementand spectralelementmethods in seism ic wavefiel modeling.Geophysics,2014,79(2)：T91-T104

22 Chen Yushu,Yang Guangwen,Ma Xiao,et al.A time-space domain stereo finit di ff erencemethod for3D scalarwavepropagation.Computers&Geosciences,2016,96：218-235

23 Liu Shaolin,LiXiaofan,WangWenshuai,etal.Am ixed-grid finit elementmethod w ith PML absorbing boundary conditions for seism ic wavemodelling.Journal ofGeophysics&Engineering,2014,11(5)：1-13

24 Clayton R,EngquistB.Absorbing boundary conditions for acoustic and elastic waveequations.Bulletin ofthe Seismological Society of America,1977,67(6)：1529-1540

25 Keller JB,Dan G.Exactnon-reflectin boundary conditions.JournalofComputational Physics,1989,82(1)：172-192

26 Chaljub E,Komatitsch D,Vilotte JP,etal.Spectral-elementanalysis in seismology.Advances in Geophysics,2007,48：365-419

27 He Ying,M in M isun,Nicholls DP.A spectralelementmethod w ith transparentboundary condition for periodic layered media scattering.JournalofScientifi Computing,2016,68(2)：772-802

28趙密.近場波動有限元模擬的應力型時域人工邊界條件及其應用.[博士論文].北京：北京工業大學,2009(Zhao M i.Stress-type time-domain artificia boundary condition for finite-elemen simulation ofnear-fiel wavemotion and itsengineering application.[PhD Thesis].Beijing：Beijing University of Technology,2009(in Chinese))

29 Berenger JP.A perfectlymatched layer for theabsorption ofelectromagnetic waves.Journal ofComputational Physics,1994,114(2)：185-200

30 Komatitsch D,Tromp J.A perfectlymatched layerabsorbingboundary condition for thesecond-order seism icwaveequation.Geophysical Journal International,2003,154(1)：146-153

31 Martin R,Komatitsch D,Gedney SD.A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seism ic wave equation.Computer Modeling in Engineering&Sciences,2008,37(3)：274-304

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

33 Xie Zhinan,Matzen R,Cristini P,etal.A perfectlymatched layer for fluid-soli problems：Application to ocean-acousticssimulations w ith solid oceanbottoms.Journalofthe AcousticalSociety ofAmerica,2016,140(1)：165-175

34廖振鵬,黃孔亮,楊柏坡等.暫態波透射邊界.中國科學：A輯,1984,(6)：556-564(Liao Zhenpeng,Huang Kongliang,Yang Baipo,et al.A transm itting boundary for transientwave.Scientia Sinica(SeriesA),1984,(6)：556-564(in Chinese))

35 Liao Zhenpeng,Wong Hongliang.A transm itting boundary for the numerical simulation of elastic wave propagation.International JournalofSoilDynamicsand Earthquake Engineering,1984,3(4)：174-183

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

37陳少林,趙宇昕.一種三維飽和土--基礎--結構動力相互作用分析方法.力學學報,2016,48(6)：1362-1371(Chen Shaolin,Zhao Yuxin.A method for three-dimensional saturated soil-foundationstructure dynam ic interaction analysis.Chinese JournalofTheoreticaland Applied Mechanics,2016,48(6)：1362-1371(in Chinese))

38李小軍,廖振鵬.時域局部透射邊界的計算飄移失穩.力學學報,1996,28(5)：627-632(LiXiaojun,Liao Zhenpeng.Thedriftinstability of local transmitting boundary in time domain.Acta Mechanica Scinica,1996,28(5)：627-632(in Chinese))

39李小軍,楊宇.透射邊界穩定性控制措施探討.巖土工程學報,2012,34(4)：641-645(LiXiaojun,Yang Yu.Measures for stability controlof transmitting boundary.Chinese Journal ofGeotechnical Engineering,2012,34(4)：641-645(in Chinese))

40周正華,廖振鵬.消除多次透射公式飄移失穩的措施.力學學報,2001,33(4)：550-554(Zhou Zhenghua,Liao Zhenpeng.A measure for eliminating drift instability of the multi-transmitting formula.Acta Mechanica Scinica,2001,33(4)：550-554(in Chinese))

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

42 Higdon RL.Absorbing boundary conditions for di ff erence approximations to themulti-dimensionalwave equation.Mathematics of Computation,1986,47(176)：437-459

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

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

IMPLEMENTATIONOFMULTI-TRANSM ITTING BOUNDARY CONDITION FOR WAVEMOTION SIMULATION BY SPECTRAL ELEMENTMETHOD:TWO DIMENSION CASE1)

Xing Haojie LiHongjing2)
(College ofCivilEngineering，Nanjing Tech University，Nanjing 211816，China)

Thespectralelementschemeofmulti-transm itting formula(MTF)which proposed by Xing and Liisdeveloped and expanded to numerical simulation of SH-typewavemotion in homogeneous linearmedium.The spectral elements adjacent to the artificia boundariesare supposed to be rectilinear quadrilateralelements,in order thateach node on the artificia boundary locates on a unique grid line pointing to the interior domain,hence the displacementof a particular artificia boundary node at a certain time can be inferred from displacements of nodes on the grid line atearlier times.Two numericalexamples,the planewave propagationw ith a certain angle(outer-source problem)and the free spread ofa pulsewaveoriginated from a point(inner-source problem),areemployed for verificatio of thiswavemotion simulationprocedurewhich combines the spectral elementmethod w ith MTF boundary.Themain parameters a ff ecting reflectio error of the MTF scheme,such as order of interpolation polynom ial,artificia wave speed and transm itting order,are investigated in time domain via a seriesof initial-value problems.The results show that the order of interpolation polynomial influence the reflectio error very little,while higher interpolation ordermay lead to better accuracy.For the choice of artificia wave speed,it is preferable to choose values equal or slightly greater than the physicalwave speed,otherw ise bigger reflectio errors come about.Transm itting order exerts themost significan impacton reflectio error,which would be reduced greatly with the increase of transm itting order of MTF,but the undesired drift instabilitymay arisewhen transmitting order reaches three or two for inner-source and outer-source problems,respectively.Although thework combining MTF boundary w ith spectralelementmethod is conducted on simulation of SH wave propagation in infinit homogeneousmedium in thispaper,ithas laid the foundation of research onmore complicated situations,such as issuesof layered soilsites,propagation of P-SV waveorRayleighwave,etc.

multi-transm itting formula,spectral elementmethod,SH-type wavemotion,spectral element scheme of MTF,drift instability

P315

A

10.6052/0459-1879-16-393

2016-12-25收稿，2017-03-01錄用,2017-03-06網絡版發表.

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

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

邢浩潔,李鴻晶.透射邊界條件在波動譜元模擬中的實現：二維波動.力學學報,2017,49(4)：894-906

Xing Haojie,LiHongjing.Implementation ofmulti-transm itting boundary condition forwavemotion simulation by spectral element method：two dimension case.Chinese JournalofTheoreticaland Applied Mechanics,2017,49(4)：894-906