透射邊界與譜元法的結合及對波動模擬精度的改進

于彥彥， 丁海平，2， 劉啟方

(1. 中國地震局 工程力學研究所，哈爾濱 150080；2. 江蘇省結構工程重點實驗室(蘇州科技學院)，江蘇 蘇州 215011)

透射邊界與譜元法的結合及對波動模擬精度的改進

于彥彥1， 丁海平1，2， 劉啟方1

(1. 中國地震局 工程力學研究所，哈爾濱 150080；2. 江蘇省結構工程重點實驗室(蘇州科技學院)，江蘇 蘇州 215011)

實現與高精度的譜元法具有一致精度的人工邊界條件，對基于譜元法的波動數值模擬具有重要意義。給出了具有多階精度的多次透射公式(MTF)與節點不等間距分布的譜元法的結合方法，并在通用的采用顯式Newmark預估校正式和并行計算技術的譜元程序SPECFEM2D和SPECFEM3D中實現。利用改進后的譜元程序模擬了二維和三維模型下波源問題的響應，模擬結果與程序自帶的黏性邊界的結果及遠置解的對比表明，結合MTF的譜元程序顯著提高了對邊界處大角度入射的體波及面波的吸收效果，消除了黏性邊界存在的漂移現象。高階MTF是與高精度的譜元法相匹配的邊界條件。

譜元法；多次透射公式；一致精度；波動數值模擬; 體波；面波

在區域或全球地震學中，對地震波傳播過程的精確模擬一直是研究的重點。為實現這一目的，一系列的數值模擬方法被相繼發展和應用，如有限差分法[1-4]、有限元法[5]、偽譜法[6-7]、譜元法[8-10]等。近年來，隨著KOMATITSCH等將譜元法推廣到地震波場的模擬并推出適用于二維和三維問題的譜元程序SPECFEM2D及SPECFEM3D，譜元法(SEM)作為一種高精度的數值方法得到越來越廣泛的應用[11-15]。譜元法的基本思想是將計算區域分解成有限個單元，在每個單元上配置Gauss-Lobatto-Legendre數值積分點(GLL點)為節點，然后在每個單元內使用譜方法或偽譜法。在保證積分點與插值點在同一位置的基礎上，以截斷的Lagrange正交多項式作為插值基函數來獲得單元的近似解。它結合了譜方法精度高、收斂速度快及有限元法幾何適應性好的優點。與傳統的有限差分、有限元等方法相比，譜元法能更好地模擬地表地形及面波的傳播，同時可以考慮各向異性及流-固耦合的情況[16]。在相同的精度下，譜元法所需劃分的單元更稀疏，避免了高階有限元可能存在的Runge現象。此外，由于單元質量陣為對角陣，因此大大減少了計算存儲，并易于在多節點上實現并行計算。

通過在單元上配置GLL點，譜元法可以實現空間上的高階精度。在時間積分方案上也有高階的Runge-Kutta算法[17]與之對應。在波動數值模擬中，為減少計算量，常常需要截取有限的計算區域，并在邊界處設置某種邊界條件以等效外部無限域的作用，其中理想的邊界條件應使得入射到邊界上的地震波能完全透射出去。對于SPECFEM2D及SPECFEM3D譜元程序，其采用的是一階精度的黏性邊界，這與譜元法的高階精度不匹配。有研究者將完美匹配層(PML)結合到譜元程序中并取得較好的吸收效果[17-18]，但它需要單獨設置邊界層和拆分場變量，施加過程較為繁瑣。因此，尋求一種與高精度的譜元法相匹配且易于實現的邊界條件具有理論和實際意義。

多次透射公式(MTF)由廖振鵬等在20世紀80年代提出[19-21]，隨后得到了極大的發展和應用[22-31]。它利用單側外行波解的一般表達式直接模擬外行波在人工邊界點上穿越邊界的過程，用與人工邊界點相鄰的若干內節點的位移或其組合來表示邊界點的位移，無需專門設置邊界層，物理意義清晰且易于實現。由于邊界點求值的計算量不多于內點求值的計算量，因此基本不增加計算時間。隨著引入內點數量的增加，可以建立高精度的高階透射公式。實踐表明二階透射公式已具有較高的工程可用精度。此外，由于透射邊界是從一般無限域模型導出，模擬的是各種單向波動的共同運動學特征，因此具有較強的普適性。

本文首先以三維問題為例對譜元法的相關理論作了總結，然后給出了多次透射公式與空間上節點不等間距分布的譜元法的結合方法，并在開源的采用Newmark時間積分格式和并行計算技術的二維和三維譜元程序SPECFEM2D及SPECFEM3D中實現應用。最后利用改進后的譜元程序模擬了二維和三維內源問題的響應，通過不同模型下采用透射邊界和遠置邊界及程序自帶的黏性邊界的模擬結果對比，研究了透射邊界與譜元法的結合對波動模擬精度的改進。

1 譜元法理論

1.1 彈性波動方程及其等效積分弱形式

介質中彈性波傳播的波動方程可表示為：

(1)

式(1)也稱為波動方程的強形式，對其直接進行求解時(如有限差分方法等)需施加自由地表應力為零的邊界條件，在地表有起伏的情況下這一過程的實現較為困難。譜元法是基于波動方程的弱形式，將式(1)左右兩邊同乘以任意的試函數w并在計算域內積分得：

(2)

對式(2)右邊第一項做分部積分，

∫Ωw▽σdΩ=∫ΩwσdΩ-∫Ω▽wσdΩ=

∫ΓwσΓdΓ+∫ΩintwσintdΩint-∫Ω▽wσdΩ

(3)

式中Ωint和Γ分別表示計算區域的內域和邊界面。對于內域，各單元間相互作用的力大小相等，方向相反，因此∫Ωintw·σintdΩint=0，則式(2)為

∫Ω▽w·σdΩ+∫ΩwfdΩ

(4)

式(4)即為譜元法的波動方程。對于自由地表，滿足應力為零的邊界條件，即σΓ=0，因此包括自由地表在內的內域節點的運動方程為

(5)

對于其余的側邊界和底邊界節點，σΓ的值不為零，根據施加的邊界條件確定。

從以上推導過程可見，譜元法自動滿足自由地表應力為零的邊界條件，無需專門施加，因此更易于模擬地表地形。

1.2 單元上的插值函數和數值積分

將計算域Ω分解為ne個單元，Ωe表示第e個單元，在每個單元上波動方程的弱形式為

(6)

譜單元內任意點的函數值用單元節點上的值插值得到，插值函數常取高階的Lagrange多項式。單元節點，即插值點，取為Gauss-Lobatto-Legendre積分點(GLL點)，即以下方程的根

(7)

式(4)中的積分采用Gauss-Lobatto-Legendre數值積分方法，積分點為GLL點。

1.3 單元質量陣

(8)

1.4 時間離散

由式(4)得全局運動方程

(9)

(10)

(11)

定義預估式

(12)

(13)

取γ=1/2,β=0，則式(10)、(11)可寫為

(14)

(15)

只有當γ=1/2時這一時間離散方案才具有二階精度，當γ為其它值時僅具有一階精度，這和譜元法空間上的高精度不匹配，因此已有研究者將高階的Runge-Kutta等時步遞推格式引入到譜元法。本文是基于上述Newmark-β方案完成MTF與譜元法的結合。

2 多次透射公式及與譜元法的結合

(16)

(17)

○代表離散網格節點，?代表透射公式用到的計算點圖1 多次透射公式原理示意圖Fig.1 Schematic of multi-transmitting formula

對于譜元法，各單元內配置的節點(GLL點)是不均勻分布的，如圖2所示，其位置由Legendre多項式一階導數的零點確定，因此文獻[33]中的插值公式不能直接使用，需稍做修改以便應用于譜元程序。筆者建議的方法有兩種：

方法1 采用式(17)，考慮到GLL點不等間距分布且在各單元中的相對位置固定的特點，若仍僅用邊界節點及相鄰的2個或4個內節點的位移插值得到計算點位移，則不能體現整個單元的運動特征，因此插值點數應取為譜單元垂直邊界方向的一條邊所包含的GLL點的個數(NGLL)。以圖2所示一維模型為例，每個單元包含有5個GLL點，則應使用這5個點進行插值；若包含6個GLL點，則使用6個點進行插值，以此類推。筆者也嘗試過只用邊界節點及相鄰的i(i

(18)

式中L為Lagrange插值多項式的次數，L=NGLL-1。式(18)中節點i(i= 0,1,…,L) 的位移對計算點j位移的貢獻，即其對應的插值系數li為

(19)

式中xd,j=x0-jcaΔt為計算點j的坐標；ca為人工波速，可根據情況調整，一般取為介質的剪切波速；xi(i=0,1,…,L)為有限元節點的坐標，Δt為時間步距。

(20)

(21)

式中ai為待定系數，ui為節點i(i=0,1,…,NGLL-1)的位移。以NGLL=5為例，將ui在x1處作泰勒展開

(22)

(23)

由于方法1相比方法2公式簡單，更易于編程，因此本文采用方法1實現多次透射公式與譜元程序的結合，并采用廖振鵬[33]建議的在透射公式中引入遠小于1的正數的方法來壓制透射邊界可能存在的漂移失穩現象。

上述結合方案在通用的二維和三維譜元程序SPECFEM2D和SPECFEM3D中完成了實現。多維情況下只需對邊界節點各個方向的運動分別應用多次透射公式即可。

?表示透射公式所需的計算點表示GLL點(NGLL=5)圖2 一維譜單元中不均勻分布的節點:Fig.2 The non equidistant nodes of 1D spectral element for NGLL=5

3 結合透射邊界的譜元程序的精度

本節將利用結合透射邊界(以二階MTF為例)的譜元程序模擬內源激勵下二維和三維模型的響應，并通過和相同模型下遠置邊界及譜元程序自帶的黏性邊界的模擬結果對比，研究改進后的譜元程序在二維及三維波動模擬問題中的精度。本文算例采用的GLL積分配置點的個數均等于5。

3.1 二維問題——均勻半空間模型

假定無限半空間均勻介質的力學參數為：cp=1 500 m/s，cs=866 m/s，ρ=1 635 kg/m3，不考慮介質阻尼。計算區域如圖3所示，其尺寸為2 000 m×2 000 m。觀測點A、B、C的坐標分別為(100,2 000)、(100,1 000)、(100,0)。采用正方形網格劃分模型，網格大小為50 m，模型共包含1 600個譜單元，26 244個節點。計算時間步長為0.001 s。

圖3 二維均勻半空間模型及觀測點的分布(m)Fig.3 2D homogeneous halfspace model and distribution of observation points(m)

波源為垂直向下的具有Riker子波形式的集中力荷載，作用于模型的(1 000,2 000) m處。荷載持時為1.0 s，其時程及頻譜如圖4所示，可見其包含的最高頻率為5 Hz左右。按照譜元法中最短波長應至少包含5個節點(GLL點)的精度條件，以上劃分的網格滿足精度要求。取模型尺寸為40 km×22 km的遠置解作為理論解，可以保證在研究時段內無邊界反射波到達。

二維譜元程序分別施加二階透射邊界、黏性邊界及遠置邊界時，模擬得到的觀測點沿Z方向的位移時程如圖5所示，其中遠置邊界的結果視為精確解。對于二維均勻彈性半空間模型在地表施加力荷載的情況(即經典Lamb問題)，空間中產生的除有P波、S波等體波外，還有首波及沿地表傳播的無頻散的Rayleigh面波。A點處面波和體波的疊加使得其響應幅值相比B、C兩點大得多，B、C兩點面波和體波分離，時程后半段為面波，可以看出，盡管沿深度方向Rayleigh波的幅值呈指數衰減，B點處面波的幅值也明顯比體波大。C點由于位于模型底部，該處面波的幅值小于體波。

圖4 施加荷載的時程及其頻譜Fig.4 The load time history and its spectrum

對A、B、C點有影響的邊界反射波主要來自左側邊界，觀測點位置的不同反映了波源入射到左邊界的角度的變化。其中經由左邊界反射到C點的波的入射角約為72°，為大角度入射的情形；而對應A點的波相對左邊界的入射角為0°，理論上講兩種邊界都應使其完全透射出去。從圖5可見對于A點，透射邊界的結果與遠置解基本完全吻合，而黏性邊界的結果存在明顯的邊界反射波，后面的波場快照顯示其為反射的Rayleigh波；對于C點，透射邊界明顯提高了透射大角度入射波的能力，對比黏性邊界，其結果與遠置解擬合更好。B點的結果也顯示了透射邊界優于黏性邊界的模擬精度。

圖5 二維均勻半空間模型，不同邊界條件下觀測點Z分量的位移時程對比Fig.5 The Z-component displacement time histories of observation points for the 2D homogeneous halfspace model under different boundary conditions

波場快照能更直觀的看出人工邊界對外行波的吸收效果。圖6給出了分別施加黏性邊界和二階透射邊界時Z方向位移的波場快照，可以看出，由于Rayleigh波的波速小于S波(c≈0.919 4cS)，因此隨著波向外傳播，面波和體波逐漸分離。T=1.8 s左右面波通過側邊界，此前兩種邊界條件下體波在人工邊界上引起的反射波均不明顯，但之后黏性邊界的結果可見明顯的面波反射波(T=2.2 s)并向計算區域內傳播，而透射邊界中未見此現象。T=3.6 s時所有波場均已通過計算區域，透射邊界的結果顯示此時整個波場的值基本全為零，但黏性邊界的結果中依然可見顯著的邊界反射波(主要是反射的面波及大角度入射時反射的體波，如圖中箭頭所示)。此外，T=2.2 s和T=3.6 s時黏性邊界對應的體波反射波的幅值也比透射邊界的大。以上結果表明黏性邊界對于面波的吸收效果較差，對體波的吸收效果也不如二階透射邊界，而且，原則上黏性邊界不適于處理面波[21]。

圖6 二維均勻半空間模型，分別施加黏性邊界和二階透射邊界時的Z向位移波場快照Fig.6 The displacement wavefield snapshots of 2D homogeneous halfspace model implemented with the viscous boundary and the second-order transmitting boundary, respectively

3.2 二維問題——成層半空間模型

模型如圖7所示，尺寸與3.1節相同，仍為2 000 m×2 000 m，其中z=1 000 m以上為介質1，以下為介質2。介質參數為：介質1：cP=1 500 m/s，cS=866 m/s，ρ=1 635 kg/m3；介質2：cP=3 000 m/s，cS=1 732 m/s，ρ=2 208 kg/m3。網格大小為50 m，波源仍為如圖4所示的Riker子波形式的集中力荷載，作用于模型的(1 000,2 000)m處。

A、B、C三個觀測點的坐標仍分別為(100,2 000) m、(100,1 000) m、(100,0) m，在除自由地表外的三個邊界上分別施加黏性邊界和二階透射邊界，將其模擬得到的三個觀測點Z分量的位移時程與遠置邊界(模型尺寸為40 km×22 km，其結果視為精確解)的結果對比示于圖8。

圖7 二維成層半空間模型及觀測點的分布(m)Fig.7 2D layered halfspace model and distribution of observation points(m)

可以看出，由于介質分界面的存在，觀測點時程明顯比均勻模型下的時程復雜，其后半段存在持續的振蕩。二維成層半空間模型下不同邊界條件的結果對比和均勻半空間的結果類似，黏性邊界在各分量均存在低估或高估觀測點響應的現象，且邊界反射波也比二階透射邊界的明顯，而透射邊界明顯改善了對面波和大角度入射體波的吸收效果，結果精度明顯提高，與遠置解符合很好。

3.3 三維問題——均勻半空間模型

計算模型如圖9所示，其尺寸為4 000 m×4 000 m×4 000 m，介質參數示于圖中，不考慮介質阻尼。震源深度1 km，其坐標為(2.0,2.0,-1.0) km，為一雙力偶點源，走向90°，傾向90°，滑動角90°。震源時間函數如圖10所示，為寬頻帶的高斯函數，上升時間為0.8 s。同時由圖10可見待模擬的最高頻率約為3 Hz，根據最短波長至少包含5個節點的精度要求，網格大小取為100 m。

圖8 二維成層半空間模型，不同邊界條件下觀測點Z分量的位移時程對比Fig.8 The Z-component displacement time histories of observation points for the 2D layered halfspace model under different boundary conditions

2個觀測點均位于地表，其y坐標為2 100 m，x坐標分別為2 300 m和3 400 m。對觀測點A，人工邊界反射P波到時為1.4 s，反射S波到時為2.6 s；對觀測點B，反射P波到時為1.0 s，反射S波到時為1.9 s。在除自由地表外的五個邊界面均施加黏性邊界或透射邊界(二階，γ2=0.001)，并取尺寸為16 km×16 km×10 km的模型的解作為精確解，這樣可以保證在5.2 s之前觀測點處無反射P波到達，9.8 s之前無反射S波到達。遠置邊界模型共包含256萬個譜單元，在30節點240線程的并行計算系統上完成求解。模擬中取計算步長為0.001 s，共模擬5 000步。三種邊界條件下2個觀測點三分量的位移時程對比如圖11，圖中人工邊界反射P波和S波到時分別標記為ARP和ARS。

坐標系的原點及三個坐標軸的正方向如圖所示震源坐標為(2, 2, -1)km，如圖中星號所示圖9 三維均勻半空間模型Fig.9 3D homogeneous halfspace model

圖10 三維模型中所用的震源時間函數(高斯函數)及其頻譜Fig.10 Source time function (Gauss function) used in the 3D problem and its spectrum

由圖11可知，對于簡單的三維均勻半空間模型，在反射波到達前，透射邊界、黏性邊界的計算結果均與遠置解相同。但在人工邊界反射波到達后，尤其是反射S波到達后，透射邊界相比黏性邊界對于入射到人工邊界上的外行波的吸收效果明顯要好。黏性邊界的計算結果主要存在兩個問題：一是某些分量的邊界反射波較為明顯，反射波的幅值甚至接近理論的反應峰值，尤其以Z分量最為顯著。觀測點A、B的徑向-垂直向運動軌跡圖顯示這部分反射波為Rayleigh面波。二是存在明顯的漂移現象，一些分量的永久位移值明顯偏離理論值，這是由于黏性邊界不能考慮介質的彈性恢復能力所致。而透射邊界的結果精度有明顯改善，邊界反射波的幅值很小，與遠置解的符合程度也更高。

(a) 觀測點A，X分量 (b) 觀測點A，Y分量 (c) 觀測點A，Z分量

(d) 觀測點B，X分量 (e) 觀測點B，Y分量 (f) 觀測點B，Z分量圖11 三維均勻半空間模型，分別施加二階透射邊界、黏性邊界時觀測點A,B三分量的位移時程與遠置解對比Fig.11 Comparisons of three-component displacement seismograms of points A and B for 3D homogeneous halfspace model under different boundary conditions

3.4 三維問題——成層半空間模型

考慮均勻半空間上有一覆蓋層的情形，模型尺寸、震源及觀測點的位置均與3.3節相同，只是包含兩層介質。介質分界面位于地表下400 m處，上層為介質1，下層為介質2，介質參數見圖12。網格大小為100 m，在除自由地表外的其余邊界面分別施加黏性邊界和二階透射邊界(γ2= 0.001)，并選擇模型尺寸為16 km×16 km×10 km的結果作為遠置解，對比不同邊界條件下A、B兩個觀測點三分量的位移時程如圖13所示。

星號代表震源所在的位置圖12 考慮覆蓋層的三維成層介質模型，其中覆蓋層厚度為400 mFig.12 3D layered halfspace model with overburden thickness equal to 400 m

可以看出，對于有覆蓋層的三維模型，透射邊界的吸收效果依然明顯優于黏性邊界。和三維半空間模型的結果類似，黏性邊界的垂直分量等存在顯著的邊界反射波，容易造成對模擬結果的誤判；而透射邊界盡管在某些分量上也有邊界反射波存在，但幅值較小，影響不大。

4 結 論

實現與譜元法的精度相匹配的人工邊界條件，對基于該方法的波動數值模擬精度的提高具有重要意義。本文利用Lagrange插值和待定系數法給出了具有多階精度的多次透射公式與網格點不等間距分布的譜元法的結合方法，并將其實現在基于顯式Newmark時間積分方案及并行計算技術的二維和三維譜元程序中。以二階透射邊界為例，利用改進的譜元程序模擬了二、三維均勻半空間及成層半空間模型在內源激勵下的響應。通過和程序自帶的黏性邊界的結果及擴展解的對比，研究了透射邊界與譜元程序結合后的精度，主要結論如下：

(1)MTF中利用GLL點的位移插值計算點的位移時，插值點數應等于單元垂直邊界方向上的邊上GLL點的個數。

(2)結合方法在二維和三維問題中均明顯提高了波動模擬的精度，改善了原有的黏性邊界對體波，尤其是對面波的吸收效果，提高了透射大角度入射波的能力，消除了黏性邊界存在的漂移現象。

(3)高階MTF是與高精度的譜元法相匹配的邊界條件。

(a) 觀測點A，X分量 (b) 觀測點A，Y分量 (c) 觀測點A，Z分量

(d) 觀測點B，X分量 (e) 觀測點B，Y分量 (f) 觀測點B，Z分量圖13 三維有覆蓋層模型，分別施加透射邊界、黏性邊界時觀測點A,B三分量的位移時程與遠置解對比Fig.13 Comparisons of three-component displacement seismograms of points A and B for 3D layered halfspace model under different boundary conditions

此外，本文模擬的主要是內源問題，對于地震波入射等外源問題，可采用波場分離法完成透射邊界與譜元法的結合。同時，本文是以二階透射邊界為例來研究結合后的精度，而和譜元法具有同階精度的更高階透射邊界的施加過程與之類似。但數值實驗表明，透射邊界的階數越高越可能出現失穩，且三維問題相比二維問題更容易失穩，因此如何在多維問題中穩定高效地實現高階透射邊界將是下一步研究的重點。

[1] VIRIEUX J. SH-wave propagation in heterogeneous media: velocity-stress finite-difference method[J]. Geophysics, 1984, 49(11): 1933-1942.

[2] VIRIEUX J. P-SV wave propagation in heterogeneous media: velocity-stress finite-difference method[J]. Geophysics, 1986, 51(4): 889-901.

[3] GRAVES R W. Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences[J]. Bull Seism Soc Am, 1996, 86(4): 1091-1106.

[4] LEE S J, CHEN H W, HUANG B S. Simulations of strong ground motion and 3D amplification effect in the Taipei basin by using a composite grid finite-difference method[J]. Bull Seism Soc Am, 2008, 98(3): 1229-1242.

[5] BAO H, BIELAK J, GHATTAS O, et al. Large-scale simulation of elastic wave propagation in heterogeneous media on parallel computers[J]. Comp Meth Appl Mech Eng, 1998, 152(1): 85-102.

[6] TESSMER E, KOSLOFF D. 3D elastic modeling with surface topo-graphy by a Chebychev spectral method[J]. Geophysics, 1994, 59(3): 464-473.

[7] FURUMURA T, KENNETT B L N, FURUMURA M. Seismic wave-field calculation for laterally heterogeneous whole earth models using the pseudo spectral method[J]. Geophys J Int, 1998, 135(3): 845-860.

[8] PATERA A T. A spectral element method for fluid dynamics: laminar flow in a channel expansion[J]. J Comput Phys, 1984, 54(3): 468-488.

[9] MADAY Y, PATERA A T. Spectral element methods for the incompressible Navier-Stokes equations, in State-of-the-art surveys on computational mechanics[M]. New York:ASME, 1989:71-143.

[10] KOMATITSCH D, VILOTTE J P. The spectral element method: An efficient tool to simulate the seismic response of 2D and 3D geological structure[J]. Bull Seism Soc Am, 1998, 88(2): 368-392.

[11] LEE S J, CHEN H W, LIU Q, et al. Three-dimensional simulations of seismic- wave propagation in the Taipei Basin with realistic topography based upon the spectral-element method[J]. Bull Seism Soc Am, 1998, 88(2): 368-392.

[12] KOMATITSCH D, LIU Q, TROMP J, et al. Simulations of ground motion in the Los Angeles basin based upon the spectral-element method[J]. Bull Seism Soc Am, 2004, 94(1): 187-206.

[13] 嚴珍珍，張懷，楊長春，等．汶川大地震地震波傳播的譜元法數值模擬研究[J]．中國科學D輯:地球科學，2009, 39(4): 393-402. YAN Zhenzhen, ZHANG Huai, YANG Changchun, et al. The numerical simulation of seismic wave propagation in Wenchuan great earthquake by the spectral element method[J].Scientia Sinica Terrae, 2009, 39 (4): 393-402.

[14] 胡元鑫，劉新榮，羅建華，等.汶川震區地震動三維地形效應的譜元法模擬[J]. 蘭州大學學報(自然科學版), 2011, 47(4): 24-32. HU Yuanxin, LIU Xinrong, LUO Jianhua, et al.Simulation of three-dimensional topographic effects on seismic ground motion in Wenchuan earthquake region based upon the spectral-element method[J]. Journal of Lanzhou University (Natural Sciences), 2011, 47(4): 24-32.

[15] 劉啟方，于彥彥，章旭斌. 施甸盆地三維地震動研究[J]. 地震工程與工程振動，2013, 33(4): 54-60. LIU Qifang, YU Yanyan, ZHANG Xubin. Three-dimensional ground motion simulation for Shidian Basin[J]. Journal of Earthquake Engineering and Engineering Vibration, 2013, 33(4): 54-60.

[16] KOMATITSCH D, TROMP J. Introduction to the spectral element method for three-dimensional seismic wave propagation[J]. Geophys J Int, 1999, 139: 806-822.

[17] XIE Z N, KOMATITSCH D, MARTIN R, et al. Improved forward wave propagation and adjoint-based sensitivity kernel calculations using a numerically stable finite-element PML[J]. Geophys J Int, 2014, 198(3): 1714-1747.

[18] KOMATITSCH D, TROMP J. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[J]. Geophys J Int, 2003, 154:146-153.

[19] LIAO Z P, WONG H L. A transmitting boundary for the numerical simulation of elastic wave propagation[J]. Soil Dyn Earthq Eng, 1984, 3:174-183.

[20] LIAO Z P, WONG H L, YANG B P, et al. A transmitting boundary for transient wave analyses[J]. Scientia Sinica (Series A), 1984, 27(10): 1063-1076.

[21] 廖振鵬, 黃孔亮, 楊柏坡,等. 暫態波透射邊界[J]. 中國科學A輯, 1984, 26(6):50-56. LIAO Zhenpeng, HUANG Kongliang, YANG Baipo, et al. A transmitting boundary condition for transient wave[J].Scientia Sinica Mathematica, 1984, 26(6): 50-56.

[22] WOLF J P. A comparison of time-domain transmitting boundaries[J].Earthquake Engng Struct Dyn,1986,14(4): 655-673.

[23] KAUSEL E. Local transmitting boundaries[J]. J Eng Mech, 1988, 114(6): 1011-1027.

[24] CHEW W C, WAGNER R L. A modified form of Liao’s absorbing boundary condition[C]//IEEE Antennas and Propagation Society International Symposium. Chicago, IL, 1992:536-539.

[25] CHENG N, CHENG C H. Relationship between Liao and Clayton-Engquist absorbing boundary conditions: acoustic case[J]. Bull Seism Soc Am, 1995, 85(3): 954-956.

[26] WAGNER R L, CHEW W C. An analysis of Liao's absorbing boundary condition[J]. J Electromagn Waves Appl, 1995, 9: 993-1009.

[27] YEO T S, KOOI P S, LEONG M S, et al. A performance assessment of Liao's absorbing boundary conditions for FDTD method[J]. Microwave and Optical Technology Letters, 1997, 16(3): 186-197.

[28] ZHANG L, YU T. A method of improving the stability of Liao’s higher-order absorbing boundary condition[J]. Progress in Electromagnetics Research M, 2012, 27: 167-178.

[29] 陳少林，廖振鵬.多次透射公式在衰減波場中的實現[J]. 地震學報, 2003, 25(3):272-279. CHEN Shaolin, LIAO Zhenpeng. Multi-transmitting formula for attenuation waves[J].Acta Seismologica Sinica, 2003, 25(3):272-279.

[30] 丁海平，宋貞霞. 震源參數對地震動相干函數的影響[J]. 振動與沖擊，2010, 29(7): 16-19. DING Haiping, SONG Zhenxia. Effect of source parameter on seismic coherency function[J]. Journal of Vibration and Shock, 2010, 29(7): 16-19.

[31] 劉啟方，金星，丁海平. 復雜場地條件下震源參數對斷層附近長周期地震動的影響[J]. 地球物理學報，2008, 51(1): 186-196. LIU Qifang, JIN Xing, DING Haiping.Effects of the source parameters on long period near-fault ground motion in the case of complex site condition[J]. Chinese J Geophys, 2008, 51(1): 186-196.

[32] HUGHES T J R. The finite-element method, linear static and dynamic finite element analysis[M]. Upper Saddle Riverm NJ:Prentice-Hall International, 1987.

[33] 廖振鵬.工程波動理論導論[M]. 北京：科學出版社, 2002.

Integration of transmitting boundary and spectral-element method and improvement on the accuracy of wave motion simulation

YU Yanyan1, DING Haiping1,2, LIU Qifang1

(1. Institute of Engineering Mechanics, China Earthquake Administration, Harbin 150080, China；2. Key Laboratory of Structure Engineering of Jiangsu Province, Suzhou University of Science and Technology, Suzhou 215011, China)

The implementation of accuracy-matched artificial boundary condition to match the high-order spectral element method (SEM) has great significance for wave motion simulations based on SEM. A method was introduced, integrating the multi-transmitting formula (MTF) with multi-order accuracy and the spectral-element method characterized by nonuniformly distributed nodes. The technique was implemented in the widely used software packages called SPECFEM2D and SPECFEM3D that adopt the explicit Newmark time scheme and parallel computing technique. The responses of 2D and 3D wave source problems were simulated by using the modified SEM packages. By comparing the results with that based on viscous boundary and the extended solution, it is shown that SEM combined with MTF significantly improves the absorbing efficiency of large-angle incidence body waves and surface waves at the boundary, and removes the drifting phenomenon of the viscous boundary. High-order MTF is a proper artificial boundary condition which can match well with the high accuracy high-order SEM.

spectral-element method; multi-transmitting formula; matched accuracy; wave motion simulation; body wave; surface wave

國家自然科學基金(51278323;51378479)

2015-07-02 修改稿收到日期：2016-01-05

于彥彥 男，博士生，1985年9月生

丁海平 男，博士，教授，1966年2月生

P315.9

A

10.13465/j.cnki.jvs.2017.02.003