高毅超，徐艷杰，金 峰，王 翔

1 清華大學(xué)水沙科學(xué)與水利水電工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100084

2 華僑大學(xué)土木工程學(xué)院，廈門(mén) 361021

3 長(zhǎng)江勘測(cè)規(guī)劃設(shè)計(jì)研究有限責(zé)任公司，武漢 430010

1 引 言

隨著我國(guó)水電開(kāi)發(fā)建設(shè)的快速發(fā)展，大壩工程的規(guī)模越來(lái)越大.西部地區(qū)許多高壩壩高達(dá)300米級(jí)，水庫(kù)向上游延伸幾十甚至幾百公里.如此巨型的水利水電工程，在遭遇強(qiáng)烈地震時(shí)，大壩結(jié)構(gòu)、河谷地基、水庫(kù)庫(kù)水以及庫(kù)底的淤砂層之間都會(huì)發(fā)生顯著的動(dòng)力相互作用，由此給大壩的地震響應(yīng)帶來(lái)巨大影響，多年來(lái)許多學(xué)者在該領(lǐng)域開(kāi)展了廣泛研究［1－6］.其中，如何準(zhǔn)確合理地確定壩面動(dòng)水壓力一直是這個(gè)領(lǐng)域的研究熱點(diǎn)和難點(diǎn).早在20世紀(jì)30年代，Westergaard［7］基于剛性直立壩面和庫(kù)水不可壓縮假設(shè)，提出了附加質(zhì)量法，該方法簡(jiǎn)單易用，至今仍被工程界廣泛接受和應(yīng)用.然而，壩面動(dòng)水壓力問(wèn)題十分復(fù)雜，大量的研究表明，庫(kù)水的可壓縮性、大壩結(jié)構(gòu)特性、淤砂層特性、地基特性以及壩體、地基、庫(kù)水和淤砂層之間的動(dòng)力相互作用等因素都會(huì)對(duì)壩面動(dòng)水壓力產(chǎn)生影響［8］.對(duì)于大壩－庫(kù)水這個(gè)耦合系統(tǒng)，其動(dòng)力分析必須考慮動(dòng)水壓力波在半無(wú)限庫(kù)水中傳播引起的輻射阻尼影響.有限單元法［9－11］、有限差分法［12－14］是模擬波動(dòng)傳播的常用數(shù)值方法，但是需要在計(jì)算域的截?cái)噙吔缟鲜┘尤斯み吔纾?5－19］.Zienkiewicz等［20］以 及 Kü?ükarslan 等［21］采 用Sommerfeld輻射邊界提出了結(jié)構(gòu)－流體動(dòng)力相互作用的數(shù)值分析方法.由于Sommerfeld輻射邊界只能吸收垂直入射波，在理論上只有零階精度，需要將其置于遠(yuǎn)離近場(chǎng)結(jié)構(gòu)的地方以獲得較好的精度.邊界元法［22－23］能夠自動(dòng)滿足無(wú)窮遠(yuǎn)處的輻射邊界條件，并且使問(wèn)題維數(shù)降低一維，是模擬動(dòng)水壓力的有效方法.然而在實(shí)際應(yīng)用中，求解邊界元法所需的基本解可能十分復(fù)雜或者不存在，這大大限制了邊界元法在實(shí)際工程中的應(yīng)用.

比例邊界有限單元法（Scaled boundary finite element method，SBFEM）［24］是一種半解析的數(shù)值方法，它兼具有限單元法和邊界元法的優(yōu)點(diǎn).該方法只需離散邊界而在徑向保持解析，使問(wèn)題維數(shù)降低一維，大大減少計(jì)算量；對(duì)于無(wú)限域的模擬，與邊界元法相比，不需要基本解即可自動(dòng)滿足無(wú)窮遠(yuǎn)處的輻射邊界條件.林皋等［25－26］將綜合考慮庫(kù)水可壓縮性和庫(kù)底吸收條件，推導(dǎo)了動(dòng)水壓力的比例邊界有限單元方程，并在頻域和時(shí)域內(nèi)進(jìn)行了求解，將其應(yīng)用于大壩－庫(kù)水動(dòng)力相互作用分析.Li［27］針對(duì)半無(wú)限等截面層狀庫(kù)水，對(duì)比例邊界有限單元方程采用對(duì)角化技術(shù)，簡(jiǎn)化卷積計(jì)算，并用于時(shí)域大壩－庫(kù)水動(dòng)力相互作用分析.

高階局部透射邊界可以在保證較高計(jì)算效率的同時(shí)達(dá)到任意階的計(jì)算精度，成為近十多年來(lái)無(wú)限域模擬中的一個(gè)研究熱點(diǎn)［28］.然而，目前大部分高階透射邊界基本屬于高頻單向漸近邊界，只能模擬行波的傳播，無(wú)法模擬快衰波的能量傳遞［29］.對(duì)于半無(wú)限庫(kù)水，存在一個(gè)截?cái)囝l率，當(dāng)激振頻率低于截?cái)囝l率時(shí)，水庫(kù)中存在快衰波，這使得部分波動(dòng)能量在未來(lái)得及向遠(yuǎn)處傳播之前已經(jīng)衰減.Prempramote和Song等［30］基于動(dòng)力剛度的連分式雙漸近解提出了一類高階雙漸近透射邊界，該透射邊界隨著階數(shù)的增加能夠在全頻范圍內(nèi)快速收斂到準(zhǔn)確解，具有很高的計(jì)算精度和計(jì)算效率.王翔等［31］將該透射邊界推廣到邊界離散的單元節(jié)點(diǎn)上，提出了動(dòng)水壓力波高階雙漸近透射邊界，并基于通用大型商業(yè)有限元軟件ABAQUS的重啟動(dòng)分析功能實(shí)現(xiàn)了該透射邊界與有限元單法的順序耦合分析模型，進(jìn)行了重力壩－庫(kù)水動(dòng)力相互作用分析，結(jié)果表明該透射邊界具有很高的精度，但由于模型基于重啟動(dòng)，無(wú)法體現(xiàn)高階雙漸近透射邊界在計(jì)算效率上的優(yōu)勢(shì).

本文將動(dòng)水壓力波高階雙漸近透射邊界直接嵌入到近場(chǎng)有限元方程中，得到一個(gè)整體控制方程，建立了混凝土壩－庫(kù)水動(dòng)力相互作用的直接耦合分析模型.基于有限元開(kāi)源軟件框架體系OpenSees編程實(shí)現(xiàn)了直接耦合分析模型，并將其應(yīng)用于重力壩和拱壩與庫(kù)水的動(dòng)力相互作用分析.

2 混凝土壩－庫(kù)水系統(tǒng)運(yùn)動(dòng)方程

典型的混凝土壩－庫(kù)水系統(tǒng)如圖1所示，其中向上游無(wú)限延伸的庫(kù)水通過(guò)一個(gè)豎直截?cái)噙吔绫环指畛山鼒?chǎng)和遠(yuǎn)場(chǎng)兩個(gè)部分.具有不規(guī)則形狀的近場(chǎng)庫(kù)水與混凝土大壩構(gòu)成廣義結(jié)構(gòu)，通過(guò)有限單元離散；遠(yuǎn)場(chǎng)規(guī)則庫(kù)水可以簡(jiǎn)化成等高或者等截面的半無(wú)限層狀介質(zhì)，如圖2所示，通過(guò)比例邊界有限單元模擬，其相似中心取在下游無(wú)窮遠(yuǎn)處.

圖1 混凝土壩－庫(kù)水系統(tǒng)Fig.1 Concrete dam－reservoir system

圖2 比例邊界有限元法應(yīng)用于半無(wú)限水庫(kù)Fig.2 SBFEM for semi－infinite reservoir

鑒于本文主要研究透射邊界的影響，不考慮材料阻尼，則混凝土壩和近場(chǎng)庫(kù)水系統(tǒng)的有限元耦合方程可以寫(xiě)成如下分塊形式

式中，u和p分別表示位移和動(dòng)水壓力，和分別表示位移和動(dòng)水壓力對(duì)時(shí)間的兩次導(dǎo)數(shù)；［M］為質(zhì)量矩陣，［K］為剛度矩陣，［Q］為混凝土壩與庫(kù)水之間的耦合矩陣，｛f｝為外荷載向量.下標(biāo)s表示壩體結(jié)構(gòu)節(jié)點(diǎn)自由度，下標(biāo)f表示不包含截?cái)噙吔绲膸?kù)水節(jié)點(diǎn)自由度，下標(biāo)b表示截?cái)噙吔绻?jié)點(diǎn)自由度.遠(yuǎn)場(chǎng)半無(wú)限庫(kù)水對(duì)近場(chǎng)的影響通過(guò)作用力－｛r｝表現(xiàn).

庫(kù)水按理想聲學(xué)流體介質(zhì)處理，不考慮庫(kù)底的吸收作用，豎直截?cái)噙吔绲陌霟o(wú)限等高或等截面層狀介質(zhì)動(dòng)水壓力波的比例邊界有限單元方程如下式［31］

［E0］｛p｝，ξξ－ ［E2］｛p｝－ ［M0］｛p｝，tt＝0，（2）式中 ｛p｝為動(dòng)水壓力，｛p｝，ξξ表示動(dòng)水壓力對(duì)比例邊界坐標(biāo)ξ的二次偏導(dǎo)，｛p｝，tt表示動(dòng)水壓力對(duì)時(shí)間的二次導(dǎo)數(shù)，系數(shù)矩陣 ［E0］、［E2］和［M0］的定義參見(jiàn)文獻(xiàn)［31］.值得指出的是，對(duì)于半無(wú)限等高或等截面層狀介質(zhì)，［M0］和［E0］存在如下關(guān)系

式中，c為動(dòng)水壓力波波速.在豎直截?cái)噙吔缟蠎?yīng)用虛功原理，可得

式中 ｛R（ω）｝和 ｛P（ω）｝分別為截?cái)噙吔缟系牡刃Ч?jié)點(diǎn)荷載和節(jié)點(diǎn)動(dòng)水壓力在頻域內(nèi)的幅值.根據(jù)動(dòng)力剛度 ［S∞（ω）］的定義，有

綜合考慮式（4）和式（5），式（2）可以改寫(xiě)成頻域動(dòng)力剛度表示的比例邊界有限單元方程

3 直接耦合分析模型的建立

3.1 動(dòng)力剛度的連分式雙漸近解

式（6）表示的比例邊界有限單元方程可以通過(guò)模態(tài)變換的方式解耦.考慮如下廣義特征值分解

式中［Λ2］為對(duì)角矩陣，其對(duì)角元素為λ2j，h為層狀介質(zhì)高度.對(duì)式（2）左乘h［Φ］T，右乘h［Φ］，并利用式（3）和式（7）可以得到如下解耦的比例邊界有限單元方程

式中a0＝ωh／c為無(wú)量綱頻率，模態(tài)動(dòng)力剛度［~S（a0）］為對(duì)角矩陣，其定義如下

根據(jù)文獻(xiàn)［31］的結(jié)果，單一模態(tài)的動(dòng)力剛度的連分式雙漸近解可以表示成

式中MH和ML分別為連分式高頻和低頻漸近階數(shù).其中式（10）和式（11）為高頻漸近解，它通過(guò)使式（8）在高頻極限a0→ ∞ 成立得到；式（13）和式（14）為低頻漸近解，它通過(guò)使高頻漸近解的殘余項(xiàng)（式（12）表示）在低頻極限a0→0成立得到.該連分式雙漸近解能在全頻范圍內(nèi)逼近動(dòng)力剛度的準(zhǔn)確解，因而能夠同時(shí)模擬行波和快衰波的傳播.

為了驗(yàn)證連分式雙漸近解的精度及其快速收斂到動(dòng)力剛度準(zhǔn)確解的特性，在頻域內(nèi)將模態(tài)動(dòng)力剛度的連分式雙漸近解與準(zhǔn)確解進(jìn)行對(duì)比，其結(jié)果見(jiàn)圖3.從中可以看出，除了在截?cái)囝l率附近的精度稍差外，MH＝ML＝2的雙漸近解基本能在全頻范圍內(nèi)逼近動(dòng)力剛度的準(zhǔn)確解；繼續(xù)增加連分式階數(shù)至MH＝ML＝4，其漸近解在全頻范圍內(nèi)快速收斂到準(zhǔn)確解.綜上，連分式雙漸近解能夠快速收斂到動(dòng)力剛度的準(zhǔn)確解，僅需要較少的階數(shù)就能達(dá)到良好的結(jié)果.根據(jù)圖3給出的結(jié)果，在時(shí)域分析中，可以選取階數(shù)為MH＝ML＝4的透射邊界進(jìn)行計(jì)算.

3.2 高階雙漸近透射邊界

參照文獻(xiàn)［31］的構(gòu)建高階透射邊界的思路，首先將模態(tài)動(dòng)力剛度的連分式雙漸近解式（10）～（14）改寫(xiě)成矩陣形式，然后對(duì)其左乘 ［Φ］T／h、右乘［Φ］后依次代入式（5），并引入輔助變量，最后再進(jìn)行傅里葉逆變換即可得時(shí)域動(dòng)水壓力波高階雙漸近透射邊界的表達(dá)式為

式中系數(shù)矩陣［C∞］和［A］定義如下

3.3 直接耦合分析模型

考慮MH階和ML階的雙漸近透射邊界，令高階殘余項(xiàng)取零.式（15）～式（19）表示的高階雙漸近透射邊界直接建立在豎直截?cái)噙吔珉x散節(jié)點(diǎn)之上，可以直接通過(guò)右端作用力｛r｝同近場(chǎng)有限元方程耦合.將式（15）～（19）代入近場(chǎng)有限元方程（1），即可得到考慮遠(yuǎn)場(chǎng)庫(kù)水輻射阻尼的混凝土壩－庫(kù)水系統(tǒng)動(dòng)力相互作用的直接耦合分析模型的控制方程

式中廣義系數(shù)矩陣 ［Mh］、［Ch］和［Kh］定義如下（分塊矩陣用黑體字母表示）

圖3 模態(tài)動(dòng)力剛度的雙漸近解（a）實(shí)部；（b）虛部.Fig.3 Continued fraction solution for modal dynamic stiffness（a）Real part；（b）Imaginary part.

其中分塊矩陣定義為

廣義未知向量與廣義外荷載向量分別為

4 算 例

OpenSees是由美國(guó)加州伯克利大學(xué)主研發(fā)的有限元開(kāi)源軟件框架體系，它基于面向?qū)ο蟮木幊趟枷?，具有豐富的單元庫(kù)和材料庫(kù)，并且可以方便地添加自定義單元和材料.借助于OpenSees開(kāi)源代碼的優(yōu)勢(shì)，作者在OpenSees中嵌入了流固耦合分析程序和高階雙漸近透射邊界，實(shí)現(xiàn)了大壩－庫(kù)水動(dòng)力相互作用直接耦合分析程序.采用以下算例進(jìn)行模型驗(yàn)證，計(jì)算中，時(shí)步積分法采用Newmark－β法，其中γ＝0.5，β＝0.25，即平均常加速度法.

4.1 重力壩算例

圖4 重力壩－庫(kù)水系統(tǒng)（a）幾何模型；（b）有限元網(wǎng)格.Fig.4 A gravity dam－reservoir system（a）Geometory；（b）Finite element mesh.

選取文獻(xiàn)［31］中的重力壩算例，重力壩和半無(wú)限庫(kù)水耦合分析系統(tǒng)如圖4所示.重力壩同近場(chǎng)不規(guī)則庫(kù)水統(tǒng)一采用8節(jié)點(diǎn)二次單元離散，其中壩體單元52個(gè)，庫(kù)水流體單元156個(gè)，豎直截?cái)噙吔绮捎?3個(gè)二次線單元離散.材料參數(shù)取值如下：壩體混凝土彈性模量E＝35GPa，泊松比ν＝0.2，質(zhì)量密度ρ＝2400kg／m3；動(dòng)水壓力波波速c＝1438.7m·s－1，庫(kù)水密度ρf＝1000kg·m－3.

根據(jù)前面分析的結(jié)果，選取透射邊界階數(shù)MH＝ML＝4.為了驗(yàn)證透射邊界的精度，采用有限元擴(kuò)展網(wǎng)格解作為參考解.在擴(kuò)展網(wǎng)格模型中，遠(yuǎn)場(chǎng)庫(kù)水有限單元離散的長(zhǎng)度為7200m，動(dòng)水壓力波在該范圍內(nèi)往返傳播的時(shí)間約為10s（對(duì)應(yīng)的無(wú)量綱時(shí)間tc／h約為90），即擴(kuò)展網(wǎng)格模型的前10s計(jì)算結(jié)果有效.計(jì)算時(shí)步取0.01s，對(duì)應(yīng)的無(wú)量綱時(shí)間約為0.09.

考慮重力壩基底受到順河向的三角脈沖加速度荷載作用，圖5是按無(wú)量綱時(shí)間定義的三角脈沖時(shí)程曲線.重力壩壩底動(dòng)水壓力時(shí)程如圖6所示，可以看出，透射邊界階數(shù)為MH＝ML＝4的直接耦合分析模型的計(jì)算結(jié)果與擴(kuò)展網(wǎng)格模型提供的參考解幾乎完全吻合，表明直接耦合分析模型具有很高的計(jì)算精度.從計(jì)算耗時(shí)上看，直接耦合分析模型耗時(shí)6s，而擴(kuò)展網(wǎng)格模型耗時(shí)約258s，表明該模型具有很高的計(jì)算效率.

圖5 三角脈沖荷載時(shí)程曲線Fig.5 Time history of triangular impulse

4.2 拱壩算例

考慮如圖7所示的拱壩－庫(kù)水系統(tǒng)，混凝土拱壩壩高H＝265m.計(jì)算中，分別將透射邊界置于離壩體1倍壩高和0.5倍壩高處，即近場(chǎng)庫(kù)水離散范圍L分別取H和0.5H.壩體和近場(chǎng)庫(kù)水均采用20節(jié)點(diǎn)六面體單元模擬，其中壩體單元192個(gè)，近場(chǎng)庫(kù)水單元分別為640個(gè)和320個(gè)；壩體和庫(kù)水的交界面采用64個(gè)8節(jié)點(diǎn)界面單元模擬，豎直截?cái)噙吔绮捎?4個(gè)8節(jié)點(diǎn)二次單元模擬.計(jì)算中，壩體混凝土和庫(kù)水的材料參數(shù)同重力壩算例.

計(jì)算中，選取透射邊界階數(shù)為MH＝ML＝4，計(jì)算步長(zhǎng)為0.02s.依然采用有限元擴(kuò)展網(wǎng)格模型作為參考解，遠(yuǎn)場(chǎng)庫(kù)水的離散范圍為7200m（約27倍壩高），即擴(kuò)展網(wǎng)格的前10s計(jì)算結(jié)果有效.

考慮拱壩基底受到順河向的El Centro地震波荷載作用，其中El Centro地震波的加速度時(shí)程曲線如圖8所示.圖9給出了直接耦合分析模型與擴(kuò)展網(wǎng)格模型前10s的計(jì)算結(jié)果，可以看出，即使將透射邊界置于0.5倍壩高處，直接耦合模型給出的計(jì)算結(jié)果仍同擴(kuò)展網(wǎng)格模型的解基本吻合.因此，在直接耦合分析模型中，可以減少近場(chǎng)庫(kù)水的離散范圍，節(jié)約計(jì)算時(shí)間.從計(jì)算耗時(shí)看，近場(chǎng)庫(kù)水離散范圍分別為1倍壩高和0.5壩高時(shí)，計(jì)算耗時(shí)為65s和51s，而擴(kuò)展網(wǎng)格模型計(jì)算耗時(shí)1847s，拱壩算例同樣體現(xiàn)了直接耦合模型具有很高的計(jì)算效率.

圖7 拱壩庫(kù)水系統(tǒng)Fig.7 An arch dam－reservoir system

圖8 El－Centro地震波時(shí)程曲線Fig.8 Time history of El Centro earthquake

5 結(jié) 論

高階雙漸近透射邊界能夠在低頻和高頻范圍內(nèi)迅速地逼近準(zhǔn)確解，具有良好的收斂性能和計(jì)算效率.本文將高階雙漸近透射邊界同有限單元法直接耦合起來(lái)，建立了時(shí)域大壩－庫(kù)水動(dòng)力相互作用的直接耦合分析模型，由于該模型的整體控制方程保留了近場(chǎng)有限元方程系數(shù)矩陣對(duì)稱稀疏的特點(diǎn)，可以方便地利用通用有限元的求解器.基于開(kāi)源軟件框架體系OpenSees，作者編寫(xiě)了時(shí)域大壩－庫(kù)水相互作用的直接耦合分析程序，并且應(yīng)用于大壩－庫(kù)水動(dòng)力相互作用分析.重力壩和拱壩的算例表明，該模型具有很高的精度和較高的計(jì)算效率，并且可以減少近場(chǎng)庫(kù)水的離散范圍，適用于模擬長(zhǎng)時(shí)間地震荷載作用下的動(dòng)水壓力波響應(yīng).

圖6 三角脈沖荷載作用下重力壩壩底動(dòng)水壓力時(shí)程Fig.6 Hydrodynamic pressure at heel of the gravity dam under triangular impulse

圖9 地震波作用下拱壩壩底動(dòng)水壓力時(shí)程Fig.9 Hydrodynamic pressure at heel of the arch dam under El Centro earthquake

（

）

［1］ Zhang C H，Pan J W，Wang J T.Influence of seismic input mechanisms and radiation damping on arch dam response.SoilDynamicsandEarthquakeEngineering，2009，29（9）：1282－1293.

［2］ Humar J L，Jablonski A M.Boundary element reservoir model for seismic analysis of gravity dams.Earthquake Engineering＆StructuralDynamics，1988，16（8）：1129－1156.

［3］ 徐艷杰，張楚漢，金峰.非線性拱壩－地基動(dòng)力相互作用的FE－BE－IBE模型.清華大學(xué)學(xué)報(bào)（自然科學(xué)版），1998，38（11）：100－104.Xu Y J，Zhang C H，Jin F.Model of FE－BE－IBE coupling for soil－structure interaction on nonlinear response of arch dams.JournalofTsinghuaUniversity（Sci＆Tech）.（in Chinese），1998，38（11）：100－104

［4］ Liu X，Xu Y，Wang G，et al.Seismic response of arch dams considering infinite radiation damping and joint opening effects.EarthquakeEngineeringandEngineeringVibration，2002，1（1）：65－73.

［5］ 徐艷杰，牟海磊，張楚漢，et al.汶川地震中寶珠寺重力壩地震響應(yīng)的三維有限元模擬.地球物理學(xué)報(bào)，2012，55（1）：293－303.Xu Y J，Mu H L，Zhang C H，et al.3Dfinite element modeling of seismic response of Baozhusi gravity dam inMs8.0Wenchuan Earthquake.ChineseJ.Geophys.（in Chinese），2012，55（1）：293－303

［6］ 程惠紅，張懷，朱伯靖，et al.新豐江水庫(kù)地震孔隙彈性耦合有限元模擬.中國(guó)科學(xué)：地球科學(xué)，2012，42（6）：905－916.Cheng H H，Zhang H，Zhu B J，et al，F(xiàn)inite element investigation of the poroelastic effect on the Xinfengjiang Reservoir－Triggered earthquake.SciChinaEarthSci，（in Chinese），2012，42（6）：905－916.

［7］ Westergaard H M. Water pressure on dams during earthquakes.Transactions of the American Society of Civil Engineering，1933，98：418－433.

［8］ 杜修力，王進(jìn)廷.動(dòng)水壓力及其對(duì)壩體地震反應(yīng)影響的研究進(jìn)展.水利學(xué)報(bào)，2001，（7）：13－21.Du X L，Wang J T.Review of studies on the hydrodynamic pressure and its effects on the seismic response of dams.JournalofHydraulicEngineering（in Chinese），2001，（7）：13－21.

［9］ 徐世浙.地球物理中的有限單元法.北京：科學(xué)出版社，1994：260－277.Xu S Z，The finite element method in Geophysics （in Chinese），Beijing：Science Press，1994：260－277.

［10］ 馮德山，陳承申，王洪華.基于混合邊界條件的有限單元法GPR正演模擬.地球物理學(xué)報(bào)，2012，55（11）：3774－3785.Feng D S，Chen C S，Wang H H，F(xiàn)inite element method GPR forward simulation based on mixed boundary condition.ChineseJ.Geophys.（in Chinese），2012，55 （11）：3774－3785.

［11］ 李偉華，劉清華，趙成剛.飽和多孔介質(zhì)三維時(shí)域黏彈性人工邊界與動(dòng)力反應(yīng)分析的顯式有限元法.地球物理學(xué)報(bào)，2010，53（10）：2460－2469.Li W H，Liu Q H，Zhao C G，Three－dimensional viscousspring boundaries in time domain and dynamic analysis using explicit finite element method of saturated porous medium.ChineseJ.Geophys.（in Chinese），2010，53（10）：2460－2469.

［12］ Lan H，Zhang Z.Three－Dimensional Wave－Field Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface.BulletinoftheSeismologicalSociety ofAmerica，2011，101（3）：1354－1370.

［13］ Lan H Q，Zhang Z J.Comparative study of the free－surface boundary condition in two－dimensional finite－difference elastic wave field simulation.JournalofGeophysicsandEngineering，2011，8（2）：275.

［14］ 蘭海強(qiáng)，劉佳，白志明，VTI介質(zhì)起伏地表地震波場(chǎng)模擬，地球物理學(xué)報(bào)，2011，54（8）：2072－2084，Lan H.Q.，Liu J.，Bai Z.M.，Wave－field simulation in VTI media with irregular free surface.ChineseJ.Geophys.（in Chinese），2011，54（8）：2072－2084.

［15］ Givoli D.Non－reflecting boundary conditions.Journalof ComputationalPhysics，1991，94（1）：1－29.

［16］ 邵秀民，藍(lán)志凌.各向異性彈性介質(zhì)中波動(dòng)方程的吸收邊界條件.地球物理學(xué)報(bào)，1995，38（S1）：56－73.Shao X M，Lan Z L.Absorbing boundary conditions for anisotropic elastic wave equations.ChineseJ.Geophys.（in Chinese），1995，38（S1）：56－73.

［17］ 廖振鵬，周正華，張艷紅.波動(dòng)數(shù)值模擬中透射邊界的穩(wěn)定實(shí)現(xiàn).地球物理學(xué)報(bào)，2002，45（4）：533－545.Liao Z P，Zhou Z H，Zhang Y H.Stable implementation of transmitting boundary in numerical simulation of wave motion.ChineseJ.Geophys.（in Chinese），2002，45（4）：533－545.

［18］ Chen J Y，Bording R P，Liu E R et al.The application of the nearly optimal sponge boundary conditions for seismic wave propagation in poroelastic media.JournalofSeismic Exploration，2010，19（1）：1－19.

［19］ Yang D H，Wang S Q，Zhang Z J et al.n－Times Absorbing Boundary Conditions for Compact Finite－Difference Modeling of Acoustic and Elastic Wave Propagation in the 2DTI Medium.BulletinoftheSeismologicalSocietyofAmerica，2003，93（6）：2389－2401.

［20］ Zienkiewicz O C， Bettess P. Fluid－structure dynamic interaction and wave forces：An introduction to numerical treatment.InternationalJournalforNumericalMethodsin Engineering，1978，13（1）：1－16.

［21］ Kü?ükarslan S，Coskun S B，TaskIn B.Transient analysis of dam－reservoir interaction including the reservoir bottom effects.JournalofFluidsandStructures，2005，20（8）：1073－1084.

［22］ Antes H，Von Estorff O.Analysis of absorption effects on the dynamic response of dam reservoir systems by boundary element methods.EarthquakeEngineering＆Structural Dynamics，1987，15（8）：1023－1036.

［23］ Soares Jr D，von Estorff O，Mansur W J.Efficient non－linear solid－fluid interaction analysis by an iterative BEM／FEM coupling.InternationalJournalforNumericalMethodsin Engineering，2005，64（11）：1416－1431.

［24］ Song C M，Wolf J P.The scaled boundary finite－element method－alias consistent infinitesimal finite－element cell method－for elastodynamics.ComputerMethodsinApplied MechanicsandEngineering，1997，147（3－4）：329－355.

［25］ Lin G，Du J G，Hu Z Q.Dynamic dam－reservoir interaction analysis including effect of reservoir boundary absorption.ScienceinChinaSeriesE：TechnologicalSciences，2007，50：1－10.

［26］ Lin G， Wang Y，Hu Z Q.An efficient approach for frequency－domain and time－domain hydrodynamic analysis of dam－reservoir systems.EarthquakeEngineering＆Structural Dynamics，2012，41（13）：1725－1749.

［27］ Li S M.Diagonalization procedure for scaled boundary finite element method in modeling semi－infinite reservoir with uniform cross－section.InternationalJournalforNumerical MethodsinEngineering，2009，80（5）：596－608.

［28］ Givoli D. High－order local non－reflecting boundary conditions：a review.WaveMotion，2004，39（4）：319－326.

［29］ Geers T. Singly and doubly asymptotic computational boundaries.In：TL G，editor.Proceedings of the IUTAM Symposium on Computational Methods for Unbounded Domains.Dordrecht：Kluwer Academic Publishers；1998.p 135－141.

［30］ Prempramote S，Song C M，Tin－Loi F，et al.High－order doubly asymptotic open boundaries for scalar wave equation.InternationalJournalforNumericalMethodsinEngineering，2009，79（3）：340－374.

［31］ Wang X，Jin F，Prempramote S，et al.Time－domain analysis of gravity dam－reservoir interaction using high－order doubly asymptotic open boundary.Computers＆Structures，2011，89（7－8）：668－680.