丁少超，黃 青

(浙江省水利水電技術(shù)咨詢中心，浙江 杭州 310020)

千百年來堰壩作為一種壅水建筑物，主要用于灌溉，是古代先民治水歷史中的重要智慧結(jié)晶[1-2]。它通過攔截水流，抬高水位來滿足引水灌溉的需求，同時在一定程度上也可對整個河道的坡降進行調(diào)整，以起到鞏固泥沙的作用；另一方面它還能夠改善和治理區(qū)域內(nèi)的生態(tài)環(huán)境和生態(tài)質(zhì)量，從而促進人與自然的和諧發(fā)展[3]。在全國中小河流中有著不計其數(shù)的堰壩，其中多為實體固定堰。近年來較多堰壩并未達到使用年限就發(fā)生局部沖毀甚至整體破壞的現(xiàn)象[4-5]。因此，在設(shè)計階段完善堰壩流場分析與穩(wěn)定性分析顯得尤為重要。

目前堰壩穩(wěn)定計算主要采用結(jié)構(gòu)靜力學的方法，在堰體穩(wěn)定分析時常選取堰壩剖面進行單寬計算，而不規(guī)則堰壩各剖面均不同，采用該方法無法準確反映堰壩及基礎(chǔ)受力情況[6]。此外，在受力分析時常將堰體、基礎(chǔ)和水體作為單獨對象，沒有考慮結(jié)構(gòu)之間的相互影響以及應力、變形的分布情況。

近些年，隨著計算機輔助工程(CAE)軟件的完善，有限元法在水利工程領(lǐng)域已經(jīng)得到廣泛應用，但考慮自由表面外部流場與應力場耦合的堰壩穩(wěn)定分析還較為欠缺[7-16]。本文基于ANSYS Workbench仿真平臺對不規(guī)則實體堰進行流固耦合仿真計算。首先，考慮結(jié)構(gòu)之間的相互作用關(guān)系，建立堰壩和基礎(chǔ)的材料模型；然后，將模型導入ANSYS CFX中引入氣液兩相流來模擬堰壩的流場以及水壓力分布；最后，將流體計算結(jié)果耦合到結(jié)構(gòu)力學分析中完成堰壩的應力及變形分析。

1 工程概況

設(shè)計中的堰壩位于浙江省東陽市南江水庫下游湖溪鎮(zhèn)境內(nèi)，堰壩寬77.3 m，長15 m，消力池長度為11.5 m，堰壩底部最低點高程為133.52 m，堰頂高程為139.5 m，設(shè)計20年一遇洪水位高程為141.1 m。堰壩上下游齒墻深入弱風化粉砂巖0.5 m，堰壩上游與消力池下游設(shè)有拋石防沖槽。堰體采用金包銀結(jié)構(gòu)，堰體內(nèi)為C20埋石混凝土，外部采用50 cm厚C20混凝土。消力池為C25F50鋼筋混凝土結(jié)構(gòu)。基礎(chǔ)材料從上至下分別為含細粒土礫、強分化粉砂巖和弱風化粉砂巖。圖1給出堰壩的整體模型和材料分區(qū)情況。

圖1 堰壩模型和材料分區(qū)

由于堰壩左右兩岸基巖面高程相差較大，弱風化粉砂巖高程相差1.6 m，因此，堰壩底部高程需根據(jù)基巖面進行調(diào)整，整體為左高右低的傾斜面。堰壩底部高程的改變將會造成堰體內(nèi)部結(jié)構(gòu)的變化，而采用傳統(tǒng)的斷面單寬進行計算顯然并不準確。本文通過建立各材料的三維模型并設(shè)置對應的材料參數(shù)以及材料之間的接觸關(guān)系，來模擬實際工程中堰壩的情況。

2 基本假定及控制方程

2.1 基本假定

數(shù)學模型建立時采用如下基本假定：①水體為不可壓縮黏性流體；②堰體為具有小變形理想彈性體；③水體與堰體接觸界面上法向位移是相等的，即水體與堰體在法向不發(fā)生直接脫離，只沿切向發(fā)生相對滑動。

2.2 基本方程

流固耦合中，控制方程包括流體控制方程和結(jié)構(gòu)運動方程。

2.2.1流體控制方程

求解自由邊界問題常采用拉格朗日、歐拉和任意拉格朗日-歐拉(ALE)3種方法[17]。流體的控制方程建立在拉格朗日法下可較好地模擬流體自由面和邊界情況，但流體變化較大時會導致網(wǎng)格畸變嚴重，求解不收斂。采用歐拉法計算時收斂情況較好，但由于網(wǎng)格不允許變形，難以模擬流體的變化情況。ALE坐標系為拉格朗日坐標系和歐拉坐標系的組合，獨立于材料和有限元網(wǎng)格區(qū)域的坐標系統(tǒng)，有限元網(wǎng)格和材料允許在該坐標系下任意移動[18]。因此，本文采用ALE法模擬堰壩水體自由面與流固耦合面的流速、流線及水壓力的分布情況。在ALE坐標系下求解流體的控制方程由連續(xù)方程和動量方程組成[19-21]：

(1)

(2)

2.2.2結(jié)構(gòu)運動基本方程

對于彈性小變形結(jié)構(gòu)，在xi坐標中的運動方程為

(3)

(4)

式中：σij為應力；εij為應變率；Fsi為固體體積力在xi方向上的分量；ρs為固體密度；si、ai分別為xi方向上的位移分量和加速度分量。

2.2.3流固耦合中的有限元方程

假設(shè)耦合體系的求解向量為X=(Xf,Xs)，其中Xf和Xs分別為定義在水體和結(jié)構(gòu)節(jié)點上的求解向量。將流體方程和結(jié)構(gòu)方程統(tǒng)一到剛度矩陣中，則流固耦合體系的有限元方程可寫為

(5)

式中：Aff、Bf分別為流場的系統(tǒng)矩陣、外部作用力矩陣；Ass、Bs分別為固體區(qū)域的系統(tǒng)矩陣、外部作用力矩陣；Asf、Afs為流固耦合矩陣；ΔXf,k、ΔXs,k分別為Xf、Xs在第k個迭代步的變化量。

2.3 邊界方程

在求解堰壩表面流時，流體域是由空氣和水體組成的兩相流，初始條件下考慮流體域自由面氣體的體積分數(shù)為1，即自由面為氣相。入口邊界條件為流體沿豎直方向的流速方程且入口水體的體積分數(shù)為1，入口高程為堰上水位高程，入口為液相；出口邊界條件為水體自由面相對壓強為零，壓強沿豎直方向由靜水壓強公式控制，出口為氣液兩相邊界。通過兩相流基本方程，求解水體自由面，即計算出水體的體積分數(shù)為1的臨界面。同時在計算流固耦合面上流體相對于固體的運動關(guān)系時，須滿足運動學和動力學兩種邊界條件[22-23]。

氣液兩相流動基本方程為

(6)

式中：vk、Vk、Ak、nk分別為第k相(k=1表示液相，k=2表示氣相)的速度、體積、表面面積、表面的單位外法向矢量；ρ2為氣體密度。

運動學邊界條件為

(7)

式中：us、vs、ws分別為水體自由面沿x方向(壩寬方向)、y方向(壩軸線方向)、z方向(豎直方向)的速度分量。

動力學邊界條件為

p0=p氣z=h(x,y,t)

(8)

式中：p0為流固耦合面上的壓強；p氣為氣體壓強。

流體與固體接觸面會發(fā)生相互作用，流體會使堰壩發(fā)生應變，而堰壩變形又會改變流體形態(tài)，因此在流固接觸面上需滿足作用力平衡(動力學條件)和位移一致性(運動學條件)：

(9)

(10)

3 計算模型

3.1 計算工況

分2種工況進行計算：①設(shè)計洪水位工況，堰上水位高于堰頂1.6 m，堰上流速1.38 m/s，荷載組合為自重、水重、靜水壓力、揚壓力；②常水位工況，堰上水位高于堰頂0.3 m，堰上流速0.34 m/s，荷載組合為自重、水重、靜水壓力、揚壓力。

3.2 網(wǎng)格模型

流體模型和堰壩模型建立后，根據(jù)求解順序先對流體模型進行網(wǎng)格劃分。在求解堰壩流固耦合問題時忽略熱傳導作用，因此，在對流體進行網(wǎng)格劃分時可先抑制堰壩模型，定義流體與堰壩接觸面，根據(jù)流場求解器選擇對應的流體網(wǎng)格劃分方式。此外，由于需要分析接觸面上的流場和水壓力分布，因此考慮對接觸面網(wǎng)格進行加密處理。流體全局網(wǎng)格單元尺寸控制為0.5 m，流體與固體接觸面網(wǎng)格單元尺寸控制為0.3 m。同時對堰頂水深處局部網(wǎng)格進行加密處理，設(shè)計洪水位工況下，局部接觸面尺寸設(shè)置為0.15 m，常水位工況下，局部接觸面尺寸設(shè)置為0.1 m。圖2(a)給出設(shè)計洪水位工況的流體計算網(wǎng)格，其中單元節(jié)點數(shù)為548 899，網(wǎng)格單元數(shù)為2 947 440。圖2(b)給出常水位工況的流體計算網(wǎng)格，其中單元節(jié)點數(shù)為431 331，網(wǎng)格單元數(shù)為2 215 298。

圖2 流體計算網(wǎng)格

在流場計算完成后進行堰壩結(jié)構(gòu)穩(wěn)定分析，需對堰壩及基礎(chǔ)模型進行網(wǎng)格劃分。在進行流固耦合分析時流體與固體接觸面會發(fā)生應力傳遞，因此需要加密流體與固體接觸面的網(wǎng)格，具體通過設(shè)置固體接觸面網(wǎng)格尺寸和膨脹層來完成。在固體表面設(shè)置3層膨脹層，增長率為1.2。堰壩外層混凝土結(jié)構(gòu)采用四面體網(wǎng)格劃分，堰體、消力池和基巖采用六面體網(wǎng)格劃分，全局網(wǎng)格尺寸控制為0.5 m。圖3給出堰壩及基礎(chǔ)的整體網(wǎng)格模型，其中單元節(jié)點數(shù)為1 552 462，網(wǎng)格單元數(shù)為501 523。

圖3 堰壩及基礎(chǔ)網(wǎng)格模型

3.3 計算參數(shù)

堰壩及基礎(chǔ)材料力學計算參數(shù)見表1。

表1 材料力學計算參數(shù)

4 計算結(jié)果分析

4.1 流場分析

圖4 堰壩縱剖面流速分布(單位：m/s)

圖5 自由表面流速分布(單位：m/s)

圖4給出堰壩縱剖面流速分布。由圖4可見，堰壩上的水體由上游至下游勢能轉(zhuǎn)化為動能，水體的流速由初始值逐漸增大，并在消力池內(nèi)發(fā)生水躍，隨后流速下降。對比圖4(a)與圖4(b)可見，相較于低流速水流，高流速水流發(fā)生摻氣的現(xiàn)象明顯，能量消散較快，相應流速變化梯度也較大。由于水流阻力，液體流速從固體壁面上零值增加到主流流速，形成一定的流速梯度。因此，水體沿堰壩下泄時最大流速發(fā)生在水體內(nèi)部。需要注意的是剖視圖中的流速涵蓋了水氣混合的流速分布。圖5為相對壓強為零的自由表面流速分布。該圖可以直觀地表現(xiàn)出流體沿堰壩下泄時的特性。從圖5(a)可以看出，堰壩下游水深小于共軛水深，發(fā)生遠離水躍，且在堰壩曲線凹側(cè)共軛水深大于其他位置處。

4.2 應力分析

在進行堰壩穩(wěn)定分析時，需導入流體分析的計算結(jié)果進行耦合。圖6為耦合后堰壩及基礎(chǔ)表面的水壓強分布，可以看出模擬的水壓強分布與水力學公式計算值基本一致。

堰壩及基礎(chǔ)最小主應力分布如圖7所示，其中拉應力為正值，壓應力為負值。由圖7可以看出，壓應力自上而下遞增，設(shè)計洪水位工況下堰壩在含細粒土礫層、強分化粉砂巖層和弱風化粉砂巖層的壓應力最大值分別為55 kPa、155 kPa和180 kPa，常水位工況下各層壓應力最大值分別為45 kPa、140 kPa和171 kPa，均能滿足各層地基承載力要求。

圖6 水壓強分布(單位：kPa)

圖7 最小主應力分布(單位：MPa)

圖8給出兩種工況下堰壩剖面最小主應力分布，設(shè)計洪水位和常水位工況下堰壩基底面壓應力最大值與最小值之比分別為1.20和1.28，均能滿足規(guī)范要求。

圖8 堰壩剖面最小主應力分布(單位：MPa)

圖9 最大主應力分布(單位：MPa)

堰壩最大主應力分布如圖9所示，可以看出堰體內(nèi)部最大主應力值為負，說明堰體內(nèi)部未出現(xiàn)拉應力。在設(shè)計洪水位工況下，堰壩表面出現(xiàn)少量拉應力，主要分布在堰壩最后一個臺階與消力池內(nèi)。此外，由于堰壩兩岸受到位移約束的影響，使得拉應力極值出現(xiàn)在堰壩與兩岸擋墻交接面上。堰壩表面和堰壩上下游位置處擋墻均為C20混凝土材料，抗拉強度為1.3 MPa，模擬的拉應力極值小于混凝土的抗拉強度。

4.3 變形分析

位移約束考慮基礎(chǔ)在弱風化基巖面下1.5 m處不發(fā)生相對位移，將堰壩與擋墻結(jié)構(gòu)視為整體進行分析，堰壩模型上下游邊界視為無窮遠，相對位移為零。堰壩整體變形分布如圖10所示，可以看出整體變形自上而下遞減，變形極值出現(xiàn)在上游拋石防沖槽右岸，由材料模型可看出該處為含細粒土礫，右岸至左岸基巖面逐漸上升，整體變形也逐漸變小，符合材料模型的變形規(guī)律。

圖10 整體變形(單位：mm)

圖11給出堰壩和基礎(chǔ)沉降變形，負號代表變形方向。設(shè)計洪水位工況下最大豎向變形為0.10 mm，常水位工況下最大變形為0.048 mm，遠小于規(guī)范允許的地基最大沉降量。最大沉降均發(fā)生在含細粒土礫層，與實際情況相符。

圖11 沉降變形(單位：mm)

5 結(jié) 論

a. 基于ANSYS CFX的有限體積法求解流場，通過引入氣液兩相流求解堰壩的表面流場和水壓力。計算結(jié)果表明該方法可有效解決自由表面外部流動問題。

b. 設(shè)計洪水位工況下堰壩下游水深小于共軛水深，發(fā)生遠離水躍，且在堰壩曲線凹側(cè)共軛水深大于其他位置處，因此堰壩下游需設(shè)消能工。

c. 設(shè)計洪水位與常水位工況下堰壩在含細粒土礫層、強分化粉砂巖層和弱風化粉砂巖層的壓應力均能滿足各層地基承載力要求，且堰體內(nèi)部未出現(xiàn)拉應力。

d. 堰壩及基礎(chǔ)變形計算結(jié)果合理，且能很好地反映堰體與基礎(chǔ)非連續(xù)性的變形情況。兩種工況下最大豎向變形小于規(guī)范允許的地基最大沉降量，且沉降極值均發(fā)生在含細粒土礫層，與實際情況相符。