許 賀 ，鄒德高 ，孔憲京 ，劉京茂

（1.大連理工大學(xué) 海岸與近海工程國家重點實驗室，遼寧 大連 116024；2.大連理工大學(xué) 水利工程學(xué)院，遼寧 大連 116024）

1 研究背景

比例邊界有限元方法（SBFEM）是一種近年新興起來的半解析的數(shù)值技術(shù)［1］，其集成了有限元法和邊界元法各自的優(yōu)勢。因此，SBFEM既適用于有限域分析，例如斷裂問題［2-3］、飽和土液化分析［4］、跨尺度分析［5-7］、磁電或壓電材料分析［8］等；又適用于無限域分析，例如結(jié)構(gòu)與地基相互作用［9］和壩-水相互作用［10-14］等。采用比例邊界有限元法模擬庫水，可建立一種高效的壩前動水壓力計算模型［10］，即直接離散壩水交界面來模擬半無限域庫水，不但將該問題的求解維數(shù)降低一維，而且自動滿足流體域無窮遠(yuǎn)處的輻射條件。該動水壓力計算模型可考慮庫水可壓縮性、地震激勵方向和庫底淤沙的波能吸收效應(yīng)等因素。由于動水壓力對于大壩抗震安全評估的重要性，地震荷載條件下作用于壩面上的動水壓力問題在試驗分析與數(shù)值計算方面一直備受關(guān)注［15-21］。在壩-庫水動力相互作用分析中，庫水可壓縮性是一個重要的影響因素與研究熱點［19-21］。Saini等［19］采用有限元法對重力壩進(jìn)行了頻域動力流固耦合分析，其結(jié)果表明，在壩-庫耦合系統(tǒng)的基頻附近，忽略庫水可壓縮性可使壩頂動位移產(chǎn)生超過30%的誤差。Porter等［20］研究了動水壓力對拱壩動力響應(yīng)的影響，結(jié)果顯示，在壩-庫耦合系統(tǒng)的基頻附近，不考慮庫水可壓縮性時，壩頂加速度響應(yīng)的誤差高達(dá)40%以上。以往的研究表明，數(shù)值分析時忽略庫水可壓縮性可能對壩體動力響應(yīng)造成明顯的誤差，不利于合理評價大壩的抗震性能，因此，進(jìn)行壩-庫水動力相互作用分析時，需要考慮庫水可壓縮性。

當(dāng)采用SBFEM離散壩前半無限域可壓縮庫水時，需要首先進(jìn)行一系列離散的動水壓力頻域求解，再經(jīng)過傅里葉逆變換將頻域解轉(zhuǎn)化為動水壓力時域脈沖響應(yīng)函數(shù)，方可計算時域的動水壓力。但是，在時域計算前的預(yù)備工作中，動水壓力頻域求解會消耗許多時間。本文建立基于SBFEM的面板壩（CFRD）與可壓縮庫水動力耦合彈塑性分析方法，在此基礎(chǔ)上，為了減少頻域求解的計算量，對動水壓力計算過程進(jìn)行簡化處理，并且驗證該簡化方案的精度和效率。

2 計算方法及其簡化方案

彈塑性面板壩采用有限元方法（FEM）離散，壩前可壓縮庫水的動水壓力采用基于SBFEM的計算模型。現(xiàn)將壩前動水壓力計算方法［10］、本文建立的面板壩與庫水動力耦合彈塑性分析方法以及動水壓力計算的簡化方案介紹如下。

2.1 基于SBFEM的壩前動水壓力計算方法 假定庫水為可壓縮、小擾動理想流體，地震荷載下，壩前動水壓力滿足方程：

忽略庫水自由表面微幅重力波：

壩水交界面上的邊界條件：

庫底和岸坡迎水面上的邊界條件：

由于采用SBFEM離散了壩前半無限域庫水，庫尾無窮遠(yuǎn)處的輻射條件自動滿足。

圖1 壩前庫水單元比例邊界坐標(biāo)

如圖1所示，根據(jù)比例邊界有限元的離散方式，將相似中心O置于水庫下游無限遠(yuǎn)處，如此可將利用二維離散網(wǎng)格來模擬三維半無限域庫水，直接利用壩水交界面上壩體的有限元網(wǎng)格生成一些由壩面（ξ=0）延伸至庫水上游無窮遠(yuǎn)處的（ξ=+∞）柱狀水體單元（ξ為解析求解的徑向坐標(biāo)）。整體坐標(biāo)（X，Y，Z）和局部坐標(biāo)（ξ，η，ζ）之間的比例邊界坐標(biāo)變換為：

對以上動水壓力的控制方程和邊界條件通過比例邊界坐標(biāo)變換和加權(quán)余量法求解，可以得到頻域動水壓力控制方程（關(guān)于ξ的二階線性常微分方程）和邊界條件，分別如下式所示：

聯(lián)立式（6）和式（7），將其轉(zhuǎn)化為特征值問題［9］，解得頻域壩前（ξ=0）動水壓力為：

其中：

2.2 時域面板壩與庫水動力耦合計算方法 壩前時域動水壓力可由頻域解通過傅里葉逆變換計算而得：

其中：

根據(jù)積分變換理論中的卷積定理，由式（10）可推得式（12），即時域動水壓力是通過時域單位脈沖響應(yīng)函數(shù)（見式（13））與壩水交界面?zhèn)鬟f而來的時域加速度卷積而得的。

程相關(guān)，因此它在tm時刻是可通過卷積計算出的已知量。

考慮地震慣性力輸入，在tm時刻壩-庫水相互作用系統(tǒng)的運動方程為：

將式（15）代入式（16），進(jìn)一步整理得：

由式（17）可知，在該強耦合計算方法中，只需將附加質(zhì)量陣［Mp］引入到壩體系統(tǒng)，即可實現(xiàn)壩與庫水系統(tǒng)的分區(qū)求解，而且不需要兩個系統(tǒng)之間的迭代計算過程，因此該方法具有高效的優(yōu)勢。該時域耦合方法可考慮壩體非線性，彈塑性面板壩與庫水系統(tǒng)動力流固耦合計算流程如圖2所示。

圖2 壩與庫水動力耦合計算流程

2.3 動水壓力計算方法的簡化方案 本文選取地震條件下壩體動力計算的常用時間步0.01 s，則動水壓力頻域求解的頻率ω范圍為0～100π。一方面，地震動的卓越頻率范圍為0～20π（0～10 Hz）；另一方面，壩體動力響應(yīng)的高階（高頻）成分可能很有限。因此，實際計算時，庫水動水壓力的高頻地震響應(yīng)的貢獻(xiàn)可能并不大。若可以縮小頻率求解范圍、截斷高頻部分，則可相應(yīng)地降低計算量。設(shè)動水壓力頻域求解的截斷頻率為ωT（ωT≤100π），則相應(yīng)的頻域求解范圍調(diào)整為0～ωT。

僅需改動上述求解動水壓力單位脈沖響應(yīng)函數(shù){hp(t)}具體流程中的第（2）步，即可改進(jìn)動水壓力計算方法，提高計算效率，進(jìn)而實現(xiàn)該動力耦合計算方法的化簡：當(dāng)0≤ωi≤ωT時，通過式（11）計算每個頻率點上的動水壓力頻域響應(yīng)函數(shù)當(dāng)ωT＜ωi≤100π時，?。憔仃嚕?/p>

在保證較高計算精度的前提下，確定ωT的取值以盡可能多的節(jié)省計算量是該簡化方案的關(guān)鍵點。本文基于面板壩的非線性動力流固耦合分析，對截斷頻率為ωT的取值進(jìn)行了敏感性分析，以下算例的截斷頻率ωT分別取為10π、20π、30π、40π、50π、100π（未截斷）。

3 面板壩算例驗證與分析

采用大連理工大學(xué)工程抗震研究所自行研發(fā)的巖土工程非線性有限元（與比例邊界有限元）程序GEODYNA［22］對典型的面板壩與可壓縮庫水進(jìn)行動力流固耦合驗證計算，通過比較截斷頻率ωT對壩前動水壓力分布規(guī)律和面板應(yīng)力極值（最大值）及其分布的影響，建議了不同壩高時，截斷頻率ωT的取值，且分析了相應(yīng)的取值依據(jù)。

3.1 大壩有限元模型與地震動輸入 采用4種不同壩高H（50 m、100 m、200 m、300 m）的規(guī)則面板堆石壩作為計算模型，其有限元網(wǎng)格如圖3所示。面板壩的上、下游面坡度均為1∶1.4和1∶1.6，梯形河谷兩岸坡度均為1∶1，其他參數(shù)詳見表1。庫水網(wǎng)格由壩體迎水面的有限元網(wǎng)格直接生成，如圖4所示。假定面板堆石壩坐落在剛性地基上（地震采用一致輸入法），采用一組《水工建筑物抗震設(shè)計規(guī)范》中的規(guī)范譜生成的人工地震波和一組Koyna實測地震波進(jìn)行計算，地震波時程如圖5和6所示，兩組地震波順河向和豎向的加速度峰值分別為1.5 m/s2和1.0 m/s2。

3.2 本構(gòu)模型及參數(shù) 土體的靜、動力本構(gòu)關(guān)系具有非線性特點［23-27］，本文筑壩堆石料采用廣義塑性本構(gòu)模型［23］，該模型可合理地反應(yīng)堆石料的材料非線性和高圍壓條件下的顆粒破碎等特性，17個參數(shù)取值見表2。面板與墊層之間的接觸面采用理想彈塑性模型，參數(shù)選取如表3所示?；炷撩姘宀捎镁€彈性模型，彈性模量E=30 GPa，泊松比ν=0.167。豎縫和周邊縫單元［28］的法向壓縮與拉伸剛度分別為25 GPa/m和5 MPa/m，切向剛度為1 MPa/m。水中波速c=1440 m/s，水體密度ρ=1.0 g/cm3，庫底和岸坡入射波的反射系數(shù)α=0.6。

圖3 面板壩有限元網(wǎng)格

圖4 庫水網(wǎng)格示意

3.3 動水壓力分布規(guī)律 以迎水面中線的動水壓力絕對值包絡(luò)為例說明分布規(guī)律。在Koyna地震作用下，100 m和300 m壩高面板壩的動水壓力分布見圖7，圖8為50 m和200 m高的面板壩在人工地震波作用下的動水壓力分布。可見，隨著截斷頻率減小，壩前動水壓力逐漸減小（精度降低），但動水壓力分布保持相似的形狀；壩高越低，動水壓力隨著截斷頻率減小而降低的幅度就越顯著；當(dāng)30π≤ωT≤50π時，動水壓力的誤差很?。划?dāng)ωT≤20π時，動水壓力的誤差則相對明顯。

圖5 人工地震波時程

圖6 Koyna地震波時程

表1 大壩有限元網(wǎng)格參數(shù)

表2 堆石料廣義塑性模型參數(shù)

表3 理想彈塑性接觸面模型參數(shù)

圖7 壩面中線動水壓力包絡(luò)（Koyna地震波）

圖8 壩面中線動水壓力包絡(luò)（人工地震波）

3.4 面板動應(yīng)力分析 混凝土面板作為核心防滲體，對面板壩的安全運行起著至關(guān)重要的作用［29］。通過對比面板動應(yīng)力極值及其分布差異，分析截斷頻率（10π，20π，30π，40π，50π，100π）對面板動應(yīng)力計算精度的影響。以下計算面板動應(yīng)力的誤差時，以截斷頻率ωT=100π（未截斷）的結(jié)果為基準(zhǔn)。

表4 面板動應(yīng)力極值（壩高50m，人工地震波）

表5 面板動應(yīng)力極值（壩高200m，人工地震波）

表6 面板動應(yīng)力極值（壩高100m，Koyna地震波）

表7 面板動應(yīng)力極值（壩高300m，Koyna地震波）

表4和表5分別為壩高50 m和200 m的面板壩在人工地震波作用下，面板動應(yīng)力極值（最大值）結(jié)果；表6和表7分別為壩高100 m和300 m的面板壩在Koyna地震波作用下，面板動應(yīng)力極值結(jié)果。表8—表11分別為相應(yīng)的面板動應(yīng)力極值的最大誤差，可見，在不同壩高條件下，隨著截斷頻率減小，面板動應(yīng)力的誤差逐漸增大；當(dāng)截斷頻率較小時（ωT≤30π），隨著壩高增大，面板動應(yīng)力極值的誤差會減小，這主要是由于壩越高其基頻越低，截斷頻率可適當(dāng)減小。

圖9為壩高50 m的面板壩在人工地震波作用下，面板壩軸向壓應(yīng)力的極值分布情況；圖10為壩高300 m面板壩在Koyna地震波作用下，面板順坡向拉應(yīng)力極值分布?？梢姡S著截斷頻率的減小，較低壩的面板高應(yīng)力區(qū)范圍變化更明顯，這主要是由于低壩的基頻較高，對截斷頻率的降低比較敏感。

表8 面板動應(yīng)力極值的最大誤差（壩高50m，人工地震波）

表9 面板動應(yīng)力極值最大誤差（壩高200m，人工地震波）

表10 面板動應(yīng)力極值的最大誤差（壩高100m，Koyna地震波）

表11 面板動應(yīng)力極值的最大誤差（壩高300m，Koyna地震波）

表12 截斷頻率及相應(yīng)計算效率

圖9 壩軸向壓應(yīng)力極值分布（壩高50m，人工地震波）

圖10 順坡向拉應(yīng)力極值分布（壩高300m，koyna地震波）

由以上分析可知，隨著面板壩高度H增加，可相應(yīng)地降低截斷頻率ωT，仍然保持具有較高計算精度，更多地節(jié)省頻域求解的計算量。根據(jù)計算結(jié)果的分析給出以下建議：當(dāng)100 m＞H≥50 m時，取ωT=40π；當(dāng)200 m＞H≥100 m時，取ωT=30π；當(dāng)H≥200 m時，取ωT=20π。根據(jù)建議，表12列出了相應(yīng)計算效率的提高程度（以ωT=100π為基準(zhǔn)），面板壩越高，則計算效率提高的越多。

4 結(jié)論

本文在面板壩與可壓縮庫水動力流固耦合時域分析方面開展了研究：（1）采用SBFEM模擬壩前庫水，并且將其與FEM離散的面板壩耦合，進(jìn)而建立了面板壩與可壓縮庫水動力耦合彈塑性分析方法，并根據(jù)壩與庫水動力耦合響應(yīng)的特點，對動水壓力計算過程進(jìn)行了簡化處理，僅需確定截斷頻率ωT，即可大幅度地降低動水壓力頻域求解的計算量，還可保證較高的計算精度。（2）地震的卓越頻率在較低的范圍內(nèi)（0～10 Hz），而且地震條件下壩體動力響應(yīng)的高階（高頻）成分對動水壓力的影響也有限。因此，庫水的高頻地震響應(yīng)貢獻(xiàn)并不大，可截斷高頻部分減少計算量。（3）隨著截斷頻率ωT減小，壩前動水壓力逐漸減?。ㄕ`差增大），但動水壓力分布保持相似的形狀；壩高越低，動水壓力隨著截斷頻率減小而降低的幅度就越顯著。在不同壩高條件下，隨著截斷頻率減小，面板動應(yīng)力的誤差逐漸增大。當(dāng)截斷頻率較小時（ωT≤30π）可發(fā)現(xiàn)明顯的規(guī)律，即隨著壩高增大，面板動應(yīng)力極值的誤差會減小；這主要是由于壩越高其基頻越低，截斷頻率可適當(dāng)減小。（4）根據(jù)計算結(jié)果分析給出以下建議：當(dāng)100 m＞H≥50 m時，取ωT=40π；當(dāng)200 m＞H≥100 m時，取ωT=30π；當(dāng)H≥200 m時，取ωT=20π。面板壩越高，則計算效率提高的越多。

本文建立的面板壩與可壓縮庫水動力耦合分析方法以及動水壓力計算的簡化方案，也適用于混凝土壩的動力計算。