基于高階雙漸近透射邊界的大壩－庫水動力相互作用直接耦合分析模型

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

1 清華大學水沙科學與水利水電工程國家重點實驗室，北京 100084

2 華僑大學土木工程學院，廈門 361021

3 長江勘測規(guī)劃設計研究有限責任公司，武漢 430010

1 引 言

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

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

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

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

2 混凝土壩－庫水系統(tǒng)運動方程

典型的混凝土壩－庫水系統(tǒng)如圖1所示，其中向上游無限延伸的庫水通過一個豎直截斷邊界被分割成近場和遠場兩個部分.具有不規(guī)則形狀的近場庫水與混凝土大壩構成廣義結構，通過有限單元離散；遠場規(guī)則庫水可以簡化成等高或者等截面的半無限層狀介質，如圖2所示，通過比例邊界有限單元模擬，其相似中心取在下游無窮遠處.

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

圖2 比例邊界有限元法應用于半無限水庫Fig.2 SBFEM for semi－infinite reservoir

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

式中，u和p分別表示位移和動水壓力，和分別表示位移和動水壓力對時間的兩次導數(shù)；［M］為質量矩陣，［K］為剛度矩陣，［Q］為混凝土壩與庫水之間的耦合矩陣，｛f｝為外荷載向量.下標s表示壩體結構節(jié)點自由度，下標f表示不包含截斷邊界的庫水節(jié)點自由度，下標b表示截斷邊界節(jié)點自由度.遠場半無限庫水對近場的影響通過作用力－｛r｝表現(xiàn).

庫水按理想聲學流體介質處理，不考慮庫底的吸收作用，豎直截斷邊界的半無限等高或等截面層狀介質動水壓力波的比例邊界有限單元方程如下式［31］

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

式中，c為動水壓力波波速.在豎直截斷邊界上應用虛功原理，可得

式中 ｛R（ω）｝和 ｛P（ω）｝分別為截斷邊界上的等效節(jié)點荷載和節(jié)點動水壓力在頻域內的幅值.根據(jù)動力剛度 ［S∞（ω）］的定義，有

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

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

3.1 動力剛度的連分式雙漸近解

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

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

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

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

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

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

3.2 高階雙漸近透射邊界

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

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

3.3 直接耦合分析模型

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

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

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

其中分塊矩陣定義為

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

4 算 例

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

4.1 重力壩算例

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

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

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

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

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

4.2 拱壩算例

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

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

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

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

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

5 結 論

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

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

圖9 地震波作用下拱壩壩底動水壓力時程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］ 徐艷杰，張楚漢，金峰.非線性拱壩－地基動力相互作用的FE－BE－IBE模型.清華大學學報（自然科學版），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.汶川地震中寶珠寺重力壩地震響應的三維有限元模擬.地球物理學報，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.新豐江水庫地震孔隙彈性耦合有限元模擬.中國科學：地球科學，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］ 杜修力，王進廷.動水壓力及其對壩體地震反應影響的研究進展.水利學報，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］ 徐世浙.地球物理中的有限單元法.北京：科學出版社，1994：260－277.Xu S Z，The finite element method in Geophysics （in Chinese），Beijing：Science Press，1994：260－277.

［10］ 馮德山，陳承申，王洪華.基于混合邊界條件的有限單元法GPR正演模擬.地球物理學報，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］ 李偉華，劉清華，趙成剛.飽和多孔介質三維時域黏彈性人工邊界與動力反應分析的顯式有限元法.地球物理學報，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］ 蘭海強，劉佳，白志明，VTI介質起伏地表地震波場模擬，地球物理學報，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］ 邵秀民，藍志凌.各向異性彈性介質中波動方程的吸收邊界條件.地球物理學報，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］ 廖振鵬，周正華，張艷紅.波動數(shù)值模擬中透射邊界的穩(wěn)定實現(xiàn).地球物理學報，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.