基于SBFEM的面板壩與可壓縮庫水動(dòng)力耦合彈塑性分析方法

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

（1.大連理工大學(xué) 海岸與近海工程國家重點(diǎn)實(shí)驗(yàn)室，遼寧 大連 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］等。采用比例邊界有限元法模擬庫水，可建立一種高效的壩前動(dòng)水壓力計(jì)算模型［10］，即直接離散壩水交界面來模擬半無限域庫水，不但將該問題的求解維數(shù)降低一維，而且自動(dòng)滿足流體域無窮遠(yuǎn)處的輻射條件。該動(dòng)水壓力計(jì)算模型可考慮庫水可壓縮性、地震激勵(lì)方向和庫底淤沙的波能吸收效應(yīng)等因素。由于動(dòng)水壓力對(duì)于大壩抗震安全評(píng)估的重要性，地震荷載條件下作用于壩面上的動(dòng)水壓力問題在試驗(yàn)分析與數(shù)值計(jì)算方面一直備受關(guān)注［15-21］。在壩-庫水動(dòng)力相互作用分析中，庫水可壓縮性是一個(gè)重要的影響因素與研究熱點(diǎn)［19-21］。Saini等［19］采用有限元法對(duì)重力壩進(jìn)行了頻域動(dòng)力流固耦合分析，其結(jié)果表明，在壩-庫耦合系統(tǒng)的基頻附近，忽略庫水可壓縮性可使壩頂動(dòng)位移產(chǎn)生超過30%的誤差。Porter等［20］研究了動(dòng)水壓力對(duì)拱壩動(dòng)力響應(yīng)的影響，結(jié)果顯示，在壩-庫耦合系統(tǒng)的基頻附近，不考慮庫水可壓縮性時(shí)，壩頂加速度響應(yīng)的誤差高達(dá)40%以上。以往的研究表明，數(shù)值分析時(shí)忽略庫水可壓縮性可能對(duì)壩體動(dòng)力響應(yīng)造成明顯的誤差，不利于合理評(píng)價(jià)大壩的抗震性能，因此，進(jìn)行壩-庫水動(dòng)力相互作用分析時(shí)，需要考慮庫水可壓縮性。

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

2 計(jì)算方法及其簡化方案

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

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

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

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

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

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

圖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)變換為：

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

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

其中：

2.2 時(shí)域面板壩與庫水動(dòng)力耦合計(jì)算方法 壩前時(shí)域動(dòng)水壓力可由頻域解通過傅里葉逆變換計(jì)算而得：

其中：

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

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

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

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

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

圖2 壩與庫水動(dòng)力耦合計(jì)算流程

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

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

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

3 面板壩算例驗(yàn)證與分析

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

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

3.2 本構(gòu)模型及參數(shù) 土體的靜、動(dòng)力本構(gòu)關(guān)系具有非線性特點(diǎn)［23-27］，本文筑壩堆石料采用廣義塑性本構(gòu)模型［23］，該模型可合理地反應(yīng)堆石料的材料非線性和高圍壓條件下的顆粒破碎等特性，17個(gè)參數(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 動(dòng)水壓力分布規(guī)律 以迎水面中線的動(dòng)水壓力絕對(duì)值包絡(luò)為例說明分布規(guī)律。在Koyna地震作用下，100 m和300 m壩高面板壩的動(dòng)水壓力分布見圖7，圖8為50 m和200 m高的面板壩在人工地震波作用下的動(dòng)水壓力分布?？梢姡S著截?cái)囝l率減小，壩前動(dòng)水壓力逐漸減小（精度降低），但動(dòng)水壓力分布保持相似的形狀；壩高越低，動(dòng)水壓力隨著截?cái)囝l率減小而降低的幅度就越顯著；當(dāng)30π≤ωT≤50π時(shí)，動(dòng)水壓力的誤差很??；當(dāng)ωT≤20π時(shí)，動(dòng)水壓力的誤差則相對(duì)明顯。

圖5 人工地震波時(shí)程

圖6 Koyna地震波時(shí)程

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

表12 截?cái)囝l率及相應(yīng)計(jì)算效率

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

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

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

4 結(jié)論

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

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