基于FSI泄洪流動與壩體穩定性耦合分析

武錫穎，徐繼偉，洪 波，鄒 渝

(1. 河南水利與環境職業學院，河南 鄭州 450008; 2. 中國農業科學院農田灌溉研究所，河南 新鄉 450000; 3. 四川大學水力學與山區河流開發保護國家重點實驗室，四川 成都 610065; 4. 四川省生態環境科學研究院，四川 成都 610041)

0 引 言

針對夏季暴雨增多入庫導致下游壩體水位增高的情況可能出現壩體失穩潰壩等情況，合理地分析泄洪過程中水流特性流動，有助于重力壩體結構穩定性的研究。大型重力壩工程的泄水建筑物通常由溢流表孔與底孔（或深孔、中孔）組成，溢流排水孔具有泄流能力強、閘門運行調度靈活、安全可靠性高等特點，是確保大壩安全最重要的泄水建筑物，很多學者進行了大量的研究。盧吉等人通過有限單元法進行考慮壩體自重的數值模擬[1]。王京等人對重力壩的排沙底孔運行方式進行研究，通過模型試驗對運行方式進行調整[2]。王亮等人對于淤地壩進行滲流模擬分析來加固壩體[3]。蔣京名等人采用理論分析及 Midas GTS NX計算機數值模擬計算方法，研究了該礦尾礦庫下壓覆礦體在地下開采過程中產生的巖層移動和變形對尾礦庫及壩體的影響[4]。同時，需要對水流壓力下的混凝土壩體力學性能進行測試[5]。但目前對于泄洪過程重力壩的結構力學特征以及水流特征作用耦合的應用案例，沒有詳細的工程背景研究，因此采用FSI(fluid-structure interaction)方法基于某水庫重力壩體在泄洪消能過程實時流動水流特性與應力耦合分析壩體的瞬態穩定性，整個過程包括壩體泄洪過程中的水流效能效果分析以及泄流過程的壩體力學特性研究。為施工中混凝土材料性能的選擇提供進一步指導[6-8]。

本研究基于FSI方法對泄水建筑物中的高速泄水水流進行模擬分析，并計算排水口射流時流體與重力壩相互作用力。通過對不同工況下（正常蓄水位、設計洪水位）泄水建筑物的泄量、水流銜接狀態、流速及水力學特征的參數進行分析，進一步了解不同工況下的消能效果，以及對蓄水池、閘體及泄水槽邊墻等動力學特性校核對工程的設計與優化提出進一步的指導和改進，保證水工構筑物的工程穩定。

1 FSI耦合理論

相比于實驗研究，數值模擬具有成本低、操作簡便、流場信息豐富等優點，使得數值模擬方法逐漸成為近年來評價水工構筑物受力下變形的主流手段。然而，數值計算中如何深層次的剖析相互作用機理以及模擬工況的參數設置決定了模擬的準確性，也成為了泄洪過程中壩體穩定性研究的關鍵技術要點。

流固耦合相互作用理論目前主要是依靠有限元與有限體積聯合求解，既能實現流場中流體的流動，同時還能考慮流體用結構作用時結構體的大變形。通過偏微分方程建立整合控制體上的強形式，甚至進一步部分積分形成弱耦合形式的系統[9-10]。具體可以通過流體求解器、結構求解器、耦合求解器三方面理解耦合過程。

1.1 流體中的作用力

不可壓縮流動制方程中連續性方程表明對于流場中的物理量φ（速度、壓力等）以及給定的無窮小體積單元δV和通量F而言，物理量φ隨時間的變化以及通量變化之和等于源項S。

通過應用質量守恒S=0，設ρf為流體密度，令φ=ρf，F=φU，其中U是通過給定的流體單元 δV的速度，由于水流過程涉及氣液兩相因此考慮VOF方法，因此質量的連續性方程變為

其中α為相的體積分數，由于不可壓縮流動的特征是密度為固定值，這意味著速度散度為0，如下所示：

將動量守恒應用于方程(1)，定義φ=ρfUi和通量F=φU，其中Ui表示速度的分量，利用質量守恒和牛頓第二定律可使等式成立，因此源項就變成單元上的力fi。如下所示：

其中 ρf=α1ρ1+(1-α2)ρ2，fi在連續介質力學的應用中，將定義作用于每體積的力和應力作用在邊界上的張量σij，其中ij表示xyz三個方向，上述公式表現流體中VOF 的守恒過程以及流體單元上的切應力的計算，為與固體域的耦合提供良好的基礎。

1.2 結構求解器的控制方程

對于固體部分的求解主要通過流體施加在固體邊界處的作用力來完成，目前主要考慮固體變形是線彈性的。其控制方程可表示為：

其中q為擴展位移場(m)，是材料參數，nu位移方向，u位移，ρs固體密度，μ剪切模量，vs固體體積，ss表示固體的表面積，fb作用在表面上的力。

整個流體與固體耦合計算過程通過上述公式在FSI求解器內計算步驟來實現。FSI求解器采用分區方法將流體和結構通過時間推進來實現耦合求解。 在時間步的處理上有顯式或隱式兩種求解方法。顯示方法具有一定的時間順初始時刻給流體與固體賦值然后先求解固體進一步將流體的作用傳遞給固體求解之后在這一部求解流體和固體的基礎上得到下一步的固體流體數據，嚴格依靠上一步的求解次序所以被描述為顯示方法，如圖1所示。

圖1 流固耦合顯示求解方法

顯式耦合方法適用于弱相互作用。OpenFOAM包含一個較弱的FSI解算器。然而，弱FSI解算器因為它相對簡單在小變形中具有重要作用，對于大變形而言他就是無法處理，需要一個強耦合的方法來實現耦合計算。

隱式方法是專門處理大變形的方法，隱式耦合方法適用于強相互作用。重點是大結構變形的有限狀態識別。他的求解次序基于上一步的計算結果同時求解當前時間步內流體與固體相互作用的過程，這一步不僅需要依靠上一步的計算結果還需要當前時間步內相互作用計算結果。既需要上一步的計算結果還依靠該時間步的計算結果的方法被稱之為隱式求解，如圖2所示。

圖2 流固耦合顯示求解方法

強耦合過程屬于求解一個大型的稀疏矩陣，這樣的隱式求解方法對于時間步的影響較小相比與弱耦合方法對時間步要求較高的限制下，隱式方法的強耦合性質非常明顯。

具體的耦合是通過在上述耦合接口中交換傳遞變量，如流體域的壓力（pr）和粘滯力（tr），結構中的位移增量（ur）和速度（vr）等變量來實現固體域變形與流體域網格運動的耦合，以及對流場以及結構變化的影響，如圖3所示。

圖3 流體與固體計算數據傳遞域求解流程

2 水庫泄洪建筑物數值模擬分析

根據某農村飲水安全水庫工程的設計方案可確定水庫正常蓄水位高程為1 410.00 m，設計洪水位高程為1 410.96 m，校核洪水位高程為1 413.15 m，總庫容為42×104m3，庫容38×104m3，最大壩高 30 m。重力壩體采用C25鋼筋混凝土，蓄水池采用C20混凝土擋墻，閘門采用鋼結構。壩體結構示意圖如圖4所示。

圖4 壩體結構示意圖

流體參數：水的密度為998 kg/m3，黏度為1×10-6Pa·s. 空氣相的粘度和密度分別為 1×10-5Pa·s和1 kg/m3，各參數均按照該重力壩泄洪的實際參數設置。

2.1 流動過程與流量校核

工況1：當壩前水位為正常蓄水位1 410.00 m時，排沙底孔開1孔并且下游無水時，設計流量為131 m3/s。進行泄流過程的數值模擬。該過程中流體域流體流動過程如圖5所示。

圖5 工況1泄流過程

圖5展示了通過排水口1泄洪的過程，通過圖5可以清楚顯示整個水流泄洪的過程。 通過圖6排水口流量曲線可以看到數值模擬的結果與泄流工況1的設計最大排水流量131 m3/s基本一致，模擬結果與設計結果誤差僅為6.1%。這表明了數值模擬的結果具有較高的準確性。

工況2：當上游為設計洪水位高程為1 410.96 m，下游設計洪水位高程為1 389.50 m時，排沙底孔全開，過閘流量為406 m3/s。水庫中水流通過閘門流動到下游的過程如圖7所示。

圖7 工況2泄流過程

不同時間點的水面流速及高度變化的情況可以看到剛開始三個排水孔同時排水水流從大壩上游水庫中的水通過排水口進入蓄水池之后在消力墩進行消能之后進入下游。其流量如圖8所示。

圖8 排水孔的總流量變化

三個排水孔全開之后每一秒的流量約為425 m3/s，與設計流量的誤差為4.6%，在該工況下的數值結果誤差是非常低的。這表明該方法模擬的結果是非常準確的能夠應用于工程計算。這樣的準確性能夠為進一步分析流動過程中壩體受力的情況提供進一步的保證。

2.2 消力墩消能效果

對于排水孔泄洪之后需要進行消能設計，按照工程的實際模型建立消能仿真模型，對工況1與工況2進行消能效果分析，可以看到工況1下進過進入蓄水池前后速度的變化。

從圖9可以看到在進入排水孔1之前到大壩排水孔射流到蓄水池這一段過程處于一個速度上升的過程，在40～60 m這一段處由于消力墩的存在到速度下降的變化，進一步通過流動過程中東壓力變化來說明其他消能效果，可以看到在經過消能墩之后流速出現明顯的下降，具體消能效果可以參考圖10的動壓變化云圖，這表明消力墩有一定的消力功能。

圖9 工況1水流沿x方向流速變化

圖10 工況1動壓力云圖變化

對工況2的模擬結果也進行相應的分析。速度沿x方向的變化如圖11所示。

圖11 工況2 沿x方向的速度變化

可以看到由于工況2下游有水，在工況2中能夠明顯看到消力墩有降低水流速度的作用消能效果相較于工況1非常明顯，結合水流動壓力變化云圖12進一步了解。

圖12 工況2動壓力變化云圖

通過不同時刻動壓力的變化云圖可以看到水流經過消力墩之后動壓力明顯降低，該過程實現了消能的效果，也表明該消力墩有消能功能能夠滿足工程中消能的設計要求。

在對工況1和工況2的動壓力的情況進行對比如圖13所示。

圖13 工況1和2 流動過程動壓力變化情況

可以展示沿水流在經過消力墩（x位置40 m處）后工況1和工況2都有動壓力衰減的過程，但工況1和2的流量與流速不同導致衰減程度不一樣，這完全符合工程消力的設計要求。

2.3 壩體穩定性分析

對比兩種工況下的應力應變以及位移變化深度了解改變工況對于壩體穩定性的影響。

在排水孔開1孔的情況下可以看到應力應變位移最大處均發生在閘門處，這與相關閘門應力的研究結果相似[11-13]，但同時壩體頂部也出現了較大的位移；而在工況2的情況下閘門全部開啟，受力都在壩體頂部上面。可以看到工況1和2的壩體頂部位置處的受力基本呈現一致的趨勢。整體受力的情況可以明顯的了解到當水位超過正常工況時要及時排水這將極大地減小閘門的受力保證其穩定性。詳細的可以從兩種工況的不同時刻平均應力與最大應力的變化進一步了解。

圖15(a)可以看到不同時間最大應力的變化是工況1遠大于工況2，這表明工況1局部出現非常大的應力值，通過圖14應力應變位移云圖可以了解到其工況1的閘門處應力較大，可知工況1的最大應力遠超過工況2。最大應力的選取考慮了壩體以及閘門為了確定壩體的應力對比，通過統計去除閘門之后壩體平均應力如圖15(b)所示工況1要仍大于工況2的應力，這表明閘門關閉后在水位升高時會導致壩體的穩定性變差，可以通過開啟開門泄洪來實現減小壩體應力的效果，有助于保證壩體的穩定性。

圖14 工況1和2穩定情況下的應力應變位移云圖

圖15 壩體應力云圖

3 結束語

對壩體泄流的穩定性進行數值分析，具體研究流體域水流的流動以及結構部分的受力穩定性。該數值模擬方面能夠定量分析實際中排孔流量大小，這也表明了FSI中速度的求解結果較為準確，從而保證流體對結構體作用力的準確提供良好的基礎，這使得流固耦合結果變得可信。

FSI流固耦合方法對于工程涉及結構以及帶有氣液界面的兩相流體運動具有重要作用，能夠快速的工程力學特性評價與流動規律研究，有助于指導工程的設計與選擇。