壩體與庫水的流固耦合分析

龔亞琦，蘇海東，崔建華

壩體與庫水的流固耦合分析

龔亞琦，蘇海東，崔建華

（長江科學(xué)院 材料與結(jié)構(gòu)研究所，武漢 430010）

基于流固耦合分析方法，分別應(yīng)用勢流模式和位移模式2種流體單元模擬庫水，并與傳統(tǒng)的附加質(zhì)量方法模擬庫水的方法做了比較，結(jié)果顯示該方法可行；相比附加質(zhì)量方法，流固耦合方法與理論解更為接近，其中位移模式流體單元對單元形態(tài)要求較高。對庫水截取區(qū)域的敏感性分析表明，流體區(qū)域選取3～5倍壩高可以得到較為滿意的結(jié)果。以重力壩的不同壩頂寬度為例，討論了需要考慮庫水可壓縮性的條件，這些成果為后續(xù)工作奠定了基礎(chǔ)。

流固耦合；動水壓力；自振頻率；可壓縮性

1 概 述

壩體-庫水之間的相互作用一直是水工結(jié)構(gòu)動力分析中的一個重要課題。以混凝土重力壩為例，當(dāng)大壩受到地震激勵時，庫水的作用對壩體的動力響應(yīng)有著顯著的影響。在地震作用下，壩體的振動帶動了庫水的運(yùn)動，庫水狀態(tài)變化產(chǎn)生的附加動水壓力，反過來又作用于壩體，影響壩體的動力特性，這種情況下，結(jié)構(gòu)受到的動壓力是不能事先確定的，它與結(jié)構(gòu)變形密切相關(guān)，因此，這種作用具有耦合特性［1］。

壩-水耦聯(lián)振動問題，最早由 Westergaard［2］在30年代提出。當(dāng)時為了滿足建壩過程中計算壩體上游面承受的附加動水壓力的需要，提出了剛性重力壩在水平地震荷載作用下的動水壓力分布公式。Westergaard公式采用剛性壩面假設(shè)，忽略了壩體的彈性變形和庫水的狀態(tài)對動力響應(yīng)的影響，該成果方法簡便、計算量小，被水工結(jié)構(gòu)抗震規(guī)范所廣泛采用并沿用至今。近年來，隨著計算機(jī)技術(shù)的發(fā)展和各種數(shù)值計算方法的提高，壩體-庫水的耦合問題可以采用流固耦合方法來求解，通過建立流體單元模擬庫水運(yùn)動狀態(tài)并考慮實際的邊界條件，對壩-水系統(tǒng)作更為真實的仿真模擬。

2 壩體-庫水?dāng)?shù)值模型

2.1 附加質(zhì)量公式

1933年，Westergaard將壩面動水壓力分布近似地用拋物線來表示，根據(jù)力矩相等條件，給出了水平地震動激勵下，具有鉛直上游面的剛性重力壩沿壩高的壩面動水壓力分布，即

式中：Pw（h）為水深h處的動水壓力；ah為水平設(shè)計地震加速度代表值；ρ為水體質(zhì)量密度；H0為庫水深。

動水壓力的特點(diǎn)與慣性力相似，所以可以用附著在壩面的一定質(zhì)量的水體來代替水的動力學(xué)效應(yīng)。為了適用任意形狀的壩體和河谷形狀，1982年，Clough［3］推廣了韋斯特伽特公式，建立了適用于有限元分析的廣義質(zhì)量表達(dá)式

式中：Mai為附加質(zhì)量；λi是壩面點(diǎn)的法向矢量；Ai為該點(diǎn)在壩面上的隸屬面積。

有限元中，流固耦合的求解有2種求解模式：①勢流模式，流體采用歐拉坐標(biāo)系下的標(biāo)量函數(shù)作為場變量，而固體采用拉格朗日坐標(biāo)系下的位移作為場變量；②位移模式，應(yīng)用位移流體單元，流體和固體均采用位移作為場變量，其好處是結(jié)構(gòu)和流體具有相同的自由度，自動滿足流固耦合邊界上的變形協(xié)調(diào)條件。

2.2 勢流模式

對于勢流單元，基于流體無旋、無粘和均質(zhì)假定，可以得到以壓力P為目標(biāo)的波動方程為流體中的聲速，kw為流體體積模量，ρ為流體的密度。

不考慮水體可壓縮性時，（3）式簡寫為▽2P＝0，其中，邊界條件如下。

流固交界面：

式中：n為壩體-庫水交界面的法向；¨Un為法向加速度。

自由液面邊界，忽略水表面重力波影響：

采用Galerkin法分別建立庫水和壩體的動平衡方程，并組合成統(tǒng)一的矩陣形式：

式中：U和P分別為位移和動水壓力；¨U和¨P分別為位移和動水壓力對時間的二階導(dǎo)數(shù)；U和P分別為位移和動水壓力對時間的一階導(dǎo)數(shù)；M，C，K分別為結(jié)構(gòu)的質(zhì)量矩陣、阻尼矩陣和剛度矩陣；M′，C′，K′分別為流體的質(zhì)量矩陣、阻尼矩陣和剛度矩陣；Q為壩-庫交界面上的耦合矩陣；F為結(jié)構(gòu)荷載向量。從公式（1）來看，由于流體的引入，合成矩陣為非對稱的。

2.3 位移模式

對于位移流體單元，流體單元與固體一樣采用位移模式，可以直接用結(jié)構(gòu)有限元方法來建立流體單元物理方程為［4，5］

式中：σx，σy，σz分別為流體單元的正應(yīng)力；εx，εy，εz分別對應(yīng)為流體單元的正應(yīng)變；τxy，τxz，τyz分別為流體單元的剪應(yīng)力；γxy，γxz，γyz分別對應(yīng)為流體單元的剪應(yīng)變；S表示流體的等效剪切模量；kw為流體的體積模量。

與固體單元不同的是，流體與流體之間，以及流體與固體單元之間并不傳遞剪力，僅只傳遞法向壓力，所以物理方程中剪切模量應(yīng)該為0，但這樣會導(dǎo)致流體矩陣的奇異［6］，為了簡便起見，可以將流體設(shè)置為相對體積模量很小的一個數(shù)值［5］，如kw·10－9。

為了檢驗上述方法的可靠性，下面將用具體算題加以討論。

3 數(shù)值算例

3.1 懸梁一側(cè)受液體作用時的自由振動

為與理論解對比，本文選用文獻(xiàn)［1］中懸臂梁一側(cè)受液體作用時的自由振動的耦聯(lián)系統(tǒng)的經(jīng)典算例（如圖1所示）。

圖1 懸臂梁與液體耦聯(lián)系統(tǒng)Fig.1 Fluid-solid coupling system of cantilever beam and water

模型中梁寬1 m，流場與梁等高（h＝10 m），梁的密度為mb，梁的彈性模量為20 GPa，泊松比為0.2；流體密度 ρ取 1 000 kg／m3，流場在順流方向截取到距梁30 m的區(qū)域范圍。理論解［1］表明，濕模態(tài)頻率與 ρh／mb的值有關(guān)，因此計算過程中，ρh／mb分別取為0.5，1，2，3。不考慮水體的可壓縮性，分別采用勢流模式、位移模式2種流體單元和附加質(zhì)量方法，進(jìn)行系統(tǒng)模態(tài)分析，計算有水頻率如表1所示。

從表1可以看出，附加質(zhì)量方法基于剛性壩面假設(shè)，計算的附加質(zhì)量偏大［7］，增大了系統(tǒng)的整體質(zhì)量，因此與流固耦合方法相比，其計算出的基頻偏小。

采用流固耦合方法計算的頻率與理論解比較接近，相對誤差在2%以內(nèi)，但是隨著ρh／mb的增大，流體的質(zhì)量相對固體質(zhì)量的比重逐漸增大，有水頻率與理論級數(shù)解的相對誤差有增大的趨勢。

需要說明的是，采用位移模式時，計算過程中出現(xiàn)了很多小頻率，需要預(yù)先判斷結(jié)構(gòu)的自振頻率區(qū)間，進(jìn)行過濾處理。這主要是由于流體剪切模量很小，形成了病態(tài)方程，Wilson［8］提出在計算中引入一個罰參數(shù)項，并對流體進(jìn)行無旋運(yùn)動限制，通過改變罰參數(shù)來解決這一問題。

表1 耦聯(lián)系統(tǒng)中梁的基頻Table 1 Natural frequency of cantilever beam in the fluid-solid coupling system

3.2 流體作用域的影響

采用流固耦合方法，流體作用域的范圍直接影響計算的精度，下面，分別采用2種流體單元模式，選用實際大壩模型，對水體區(qū)域的截取范圍作敏感性分析。

選取一重力壩［9］，如圖2所示，壩高105 m，壩頂寬10 m，壩底寬80 m，上游壩面豎直，沿壩軸線方向取寬度100 m。壩體混凝土：彈模為1.96×1010Pa，密度為 2 400 kg／m3，泊松比為0.2；水體：密度為1 000 kg／m3，水中縱波傳播速度 c＝1 440 m／s，水深為100 m。計算過程中，不考慮庫水表面重力波和壩基的影響，分別選取不同的流體域長度：1～6倍壩高，對其數(shù)值結(jié)果的穩(wěn)定性進(jìn)行比較分析。由于勢流模式和位移模式計算的壩體振型不盡相同，選取不同工況下相同振型對應(yīng)的三階壩體頻率列于表2。

圖2 重力壩-庫水耦合有限元模型Fig.2 Finite element model of gravity dam-reservoir

表2 不同流體域長度壩體自振頻率Table 2 Natural frequencies of dam by different reservoir models with different water areas rad／s

水體區(qū)域敏感性分析表明：采用勢流模式流體單元時，上游取3倍壩高可以得到穩(wěn)定的數(shù)值結(jié)果；采用位移模式流體單元時，計算結(jié)果隨流體區(qū)域的選擇有波動，5倍壩高基本可以得到滿意的結(jié)果；2種流體單元模式計算的前二階頻率相對誤差在4%以內(nèi)。

3.3 水體可壓縮性的影響

在壩水相互作用過程中，是否可以忽略庫水可壓縮性影響一直是一個有爭議的問題。與不考慮庫水可壓縮性相比，考慮庫水可壓縮性時的分析要復(fù)雜的多。考慮庫水可壓縮性的影響之后，將使得原來所需求解的拉普拉斯方程轉(zhuǎn)為相應(yīng)的波動方程，在求解方面增加了很大的難度。本文將從結(jié)構(gòu)的角度，就庫水的壓縮性在壩-水相互作用中的具體影響作進(jìn)一步探討。

模型仍然選用上述重力壩模型，壩頂寬度B作為變量，分別采用不同壩頂寬度（10，20，30，40，50，60 m），比較附加質(zhì)量法、不可壓縮庫水單元以及可壓縮庫水單元對結(jié)構(gòu)頻率振型的影響。

從圖3可以看出，隨著壩頂寬度的增加，壩體質(zhì)量逐漸增大，流固耦合方法與附加質(zhì)量方法計算的自振頻率的差距越來越大，表1的結(jié)果也說明了這一點(diǎn)，此時用流固耦合方法計算系統(tǒng)的動力特性更準(zhǔn)確。

同時，隨著壩體質(zhì)量的增加，結(jié)構(gòu)的自振頻率逐漸降低，考慮和不考慮水體可壓縮性，兩者計算的結(jié)果越來越接近。

圖4是不同壩體彈性模量條件下，附加質(zhì)量法、不可壓縮庫水單元以及可壓縮庫水單元對結(jié)構(gòu)頻率振型的影響。模型中壩體彈性模量（E）分別選取為10，20，30 GPa時，結(jié)構(gòu)的自振頻率逐漸增大，相比不可壓縮庫水單元，考慮庫水可壓縮性時求得的自振頻率增加的幅度較慢。

圖3 不同壩頂寬度（B）對自振頻率的影響Fig.3 Effect of different top dam widths（B）on natural frequency

圖4 不同壩體動彈性模量（E）對自振頻率的影響Fig.4 Effect of different dynamic elastic moduli of dam（E）on natural frequency

對比圖3和圖4可看出，當(dāng)結(jié)構(gòu)自振頻率較高時，庫水壓縮性對結(jié)構(gòu)動力特性影響較為顯著。庫水可壓縮性的影響還與地形、淤泥、基巖和共振頻率等綜合條件有關(guān)［7，9，10］

4 結(jié) 論

分別采用勢流模式和位移模式2種流體單元模擬庫水，并與壩體結(jié)構(gòu)單元耦合模擬壩體-庫水的相互作用。

（1）經(jīng)典算例結(jié)果顯示，流固耦合方法與理論解更為接近，采用流體單元模擬庫水的方法是可行的；采用位移模式時，需要加入無旋約束，過濾小頻率。

（2）庫水截取區(qū)域的敏感性分析表明：流體區(qū)域選取3～5倍壩高可以得到較為滿意的結(jié)果；2種流體單元模式計算的前二階頻率相對誤差在4%以內(nèi)。

（3）當(dāng)壩體結(jié)構(gòu)自振頻率比較高時，庫水壓縮性對結(jié)構(gòu)動力特性的影響較顯著。

［1］ 黃玉盈.結(jié)構(gòu)振動分析基礎(chǔ)［M］.武漢：華中工學(xué)院出版社，1988.（HUANG Yu-ying.Structural Vibration Analysis［M］.Wuhan：Huazhong University of Science＆Technology Press，1988.（in Chinese））

［2］ WESTERGAARD H M.Water Pressures on Dams During Earthquakes［J］.Transactions of the American Society of Civil Engineers，1933，98：418－472.

［3］ CLOUGH H M.Reservoir Interaction Effects on the Dynamic Response of Arch Dam［J］.Proceedings of China-US Bilateral Workshop on Earthquake Engineering，1982：58－84.

［4］ 李遇春，樓夢麟，尚 偉，等.大型渡槽抗震分析中流體的位移有限元模式［J］.水利學(xué)報，2003，2：93－97.（LI Yu-chun，LOU Meng-lin，SHANGWei，et al.Displacement-Based Fluid Finite Element for Seismic-Resistant Analysis of Large-Scale Aqueduct［J］.Journal of Hydraulic Engineering，2003，2：93－97.（in Chinese））

［5］ 李遇春，樓夢麟，周 成.大型渡槽的豎向地震效應(yīng)分析［J］.地震工程與工程振動，2008，28（2）：102－107.（LI Yu-chun，LOU Meng-lin，ZHOU Cheng.Analysis of Vertical Seismic Effects for Large-Scale Aqueduct Bridges［J］.Journal of Earthquake Engineering and Engineer Vibration，2008，28（2）：102－107.（in Chinese））

［6］ AKKAS N，AKAY H U，YILMAZ C.Applicability of General Purpose Finite Element Programs in Solid-Fluid Interaction Problems［J］.Computers＆Structures，1979，10：773－783.

［7］ 陳厚群，侯順載，楊大偉.地震條件下拱壩庫水相互作用的實驗研究［J］.水利學(xué)報，1989，7（7）：29－39.（CHEN Hou-qun，HOU Shun-zai，YANG Da-wei.Study on Arch Dam-Reservoir Water Interaction Under Earthquake Condition［J］.Journal of Hydraulic Engineering，1989，7（7）：29－39.（in Chinese））

［8］ WILSON E L，KHALVATI M.Finite Elements for Dynamic Analysis of Fluid-Solid Systems［J］.International Journal of Numerical Methods in Engineering，1983，19：1657－1668.

［9］ 高瑞強(qiáng)，宮必寧，郭建平.壩-水相互作用中水的壓縮性影響研究［J］.水利與建筑工程學(xué)報，2008，6（1）：106－108.（GAO Rui-qiang，GONG Bi-ning，GUO Jianping.Study on Effects of Compressibility of Water in Dam-Water Interaction［J］.Journal of Water Resources and Architectural Engineering，2008，6（1）：106－108.（in Chinese））

［10］劉浩吾.混凝土壩動水壓力與庫水可壓縮性效應(yīng)［J］.水利水電科技進(jìn)展，2002，22（2）：10－13.（LIU Haowu.Hydrodynamic Pressure on Concrete Dams and Reservoir Water Compressibility Effect［J］.Advances in Science and Technology of Water Resources，2002，22（2）：10－13.（in Chinese） ）

Fluid-Solid Coupling Analysis of Dam-Reservoir Interaction

GONG Ya-qi，SU Hai-dong，CUI Jian-hua
（Yangtze River Scientific Research Institute，Wuhan 430010，China）

Dam-reservoir interaction was analyzed in this paper using fluid-solid coupling method.The reservoir water was modeled respectively by two types of finite element，namely pressure-based element and displacementbased element.The fluid-solid coupling method was compared with traditional Westergaard additional mass method and the comparison showed that the former was closer to the theoretical solution than the latter.The displacementbased element demanded more regular element shape than pressure-based element.Moreover，sensitivity analysis of selected areas of the reservoir showed that satisfactory results can be obtained if the length of reservoir area to be modeled was designed only 3－5 times the height of the dam.Finally，the compressibility of water was also explored taking the gravity dam structured model with different dam crest widths as an example.The present results will be helpful for the following research.

fluid-solid coupling；hydrodynamic pressure；natural frequency；compressibility

TV312

A

1001－5485（2011）06－0063－04

2010-05-24

長江科學(xué)院博士啟動項目（CKSQ2010082）

龔亞琦（1980-），男，湖北武漢人，工程師，主要從事水工結(jié)構(gòu)數(shù)值分析與計算工作，（電話）027-82829754（電子信箱）gong＿yaqi＠163.com。

（編輯：劉運(yùn)飛）