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

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

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

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

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

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

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

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

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

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

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

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

1 譜元法理論

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

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

(1)

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

(2)

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

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

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

(3)

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

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

(4)

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

(5)

對于其余的側(cè)邊界和底邊界節(jié)點，σΓ的值不為零，根據(jù)施加的邊界條件確定。

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

1.2 單元上的插值函數(shù)和數(shù)值積分

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

(6)

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

(7)

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

1.3 單元質(zhì)量陣

(8)

1.4 時間離散

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

(9)

(10)

(11)

定義預估式

(12)

(13)

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

(14)

(15)

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

2 多次透射公式及與譜元法的結(jié)合

(16)

(17)

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

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

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

(18)

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

(19)

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

(20)

(21)

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

(22)

(23)

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

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

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

3 結(jié)合透射邊界的譜元程序的精度

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

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

假定無限半空間均勻介質(zhì)的力學參數(shù)為：cp=1 500 m/s，cs=866 m/s，ρ=1 635 kg/m3，不考慮介質(zhì)阻尼。計算區(qū)域如圖3所示，其尺寸為2 000 m×2 000 m。觀測點A、B、C的坐標分別為(100,2 000)、(100,1 000)、(100,0)。采用正方形網(wǎng)格劃分模型，網(wǎng)格大小為50 m，模型共包含1 600個譜單元，26 244個節(jié)點。計算時間步長為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個節(jié)點(GLL點)的精度條件，以上劃分的網(wǎng)格滿足精度要求。取模型尺寸為40 km×22 km的遠置解作為理論解，可以保證在研究時段內(nèi)無邊界反射波到達。

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

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

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

圖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左右面波通過側(cè)邊界，此前兩種邊界條件下體波在人工邊界上引起的反射波均不明顯，但之后黏性邊界的結(jié)果可見明顯的面波反射波(T=2.2 s)并向計算區(qū)域內(nèi)傳播，而透射邊界中未見此現(xiàn)象。T=3.6 s時所有波場均已通過計算區(qū)域，透射邊界的結(jié)果顯示此時整個波場的值基本全為零，但黏性邊界的結(jié)果中依然可見顯著的邊界反射波(主要是反射的面波及大角度入射時反射的體波，如圖中箭頭所示)。此外，T=2.2 s和T=3.6 s時黏性邊界對應的體波反射波的幅值也比透射邊界的大。以上結(jié)果表明黏性邊界對于面波的吸收效果較差，對體波的吸收效果也不如二階透射邊界，而且，原則上黏性邊界不適于處理面波[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節(jié)相同，仍為2 000 m×2 000 m，其中z=1 000 m以上為介質(zhì)1，以下為介質(zhì)2。介質(zhì)參數(shù)為：介質(zhì)1：cP=1 500 m/s，cS=866 m/s，ρ=1 635 kg/m3；介質(zhì)2：cP=3 000 m/s，cS=1 732 m/s，ρ=2 208 kg/m3。網(wǎng)格大小為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，其結(jié)果視為精確解)的結(jié)果對比示于圖8。

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

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

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

計算模型如圖9所示，其尺寸為4 000 m×4 000 m×4 000 m，介質(zhì)參數(shù)示于圖中，不考慮介質(zhì)阻尼。震源深度1 km，其坐標為(2.0,2.0,-1.0) km，為一雙力偶點源，走向90°，傾向90°，滑動角90°。震源時間函數(shù)如圖10所示，為寬頻帶的高斯函數(shù)，上升時間為0.8 s。同時由圖10可見待模擬的最高頻率約為3 Hz，根據(jù)最短波長至少包含5個節(jié)點的精度要求，網(wǎng)格大小取為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節(jié)點240線程的并行計算系統(tǒng)上完成求解。模擬中取計算步長為0.001 s，共模擬5 000步。三種邊界條件下2個觀測點三分量的位移時程對比如圖11，圖中人工邊界反射P波和S波到時分別標記為ARP和ARS。

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

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

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

(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節(jié)相同，只是包含兩層介質(zhì)。介質(zhì)分界面位于地表下400 m處，上層為介質(zhì)1，下層為介質(zhì)2，介質(zhì)參數(shù)見圖12。網(wǎng)格大小為100 m，在除自由地表外的其余邊界面分別施加黏性邊界和二階透射邊界(γ2= 0.001)，并選擇模型尺寸為16 km×16 km×10 km的結(jié)果作為遠置解，對比不同邊界條件下A、B兩個觀測點三分量的位移時程如圖13所示。

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

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

4 結(jié) 論

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

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

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

(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

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

[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] 嚴珍珍，張懷，楊長春，等．汶川大地震地震波傳播的譜元法數(shù)值模擬研究[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] 胡元鑫，劉新榮，羅建華，等.汶川震區(qū)地震動三維地形效應的譜元法模擬[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] 廖振鵬, 黃孔亮, 楊柏坡,等. 暫態(tài)波透射邊界[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] 陳少林，廖振鵬.多次透射公式在衰減波場中的實現(xiàn)[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] 丁海平，宋貞霞. 震源參數(shù)對地震動相干函數(shù)的影響[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] 劉啟方，金星，丁海平. 復雜場地條件下震源參數(shù)對斷層附近長周期地震動的影響[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