基于FSI泄洪流動與壩體穩(wěn)定性耦合分析

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

(1. 河南水利與環(huán)境職業(yè)學(xué)院，河南 鄭州 450008; 2. 中國農(nóng)業(yè)科學(xué)院農(nóng)田灌溉研究所，河南 新鄉(xiāng) 450000; 3. 四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護國家重點實驗室，四川 成都 610065; 4. 四川省生態(tài)環(huán)境科學(xué)研究院，四川 成都 610041)

0 引 言

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

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

1 FSI耦合理論

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

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

1.1 流體中的作用力

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

通過應(yīng)用質(zhì)量守恒S=0，設(shè)ρf為流體密度，令φ=ρf，F(xiàn)=φU，其中U是通過給定的流體單元 δV的速度，由于水流過程涉及氣液兩相因此考慮VOF方法，因此質(zhì)量的連續(xù)性方程變?yōu)?/p>

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

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

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

1.2 結(jié)構(gòu)求解器的控制方程

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

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

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

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

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

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

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

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

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

圖3 流體與固體計算數(shù)據(jù)傳遞域求解流程

2 水庫泄洪建筑物數(shù)值模擬分析

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

圖4 壩體結(jié)構(gòu)示意圖

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

2.1 流動過程與流量校核

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

圖5 工況1泄流過程

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

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

圖7 工況2泄流過程

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

圖8 排水孔的總流量變化

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

2.2 消力墩消能效果

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

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

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

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

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

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

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

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

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

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

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

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

2.3 壩體穩(wěn)定性分析

對比兩種工況下的應(yīng)力應(yīng)變以及位移變化深度了解改變工況對于壩體穩(wěn)定性的影響。

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

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

圖14 工況1和2穩(wěn)定情況下的應(yīng)力應(yīng)變位移云圖

圖15 壩體應(yīng)力云圖

3 結(jié)束語

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

FSI流固耦合方法對于工程涉及結(jié)構(gòu)以及帶有氣液界面的兩相流體運動具有重要作用，能夠快速的工程力學(xué)特性評價與流動規(guī)律研究，有助于指導(dǎo)工程的設(shè)計與選擇。