基于流固耦合的氣液兩相流豎直彎管振動特性數(shù)值仿真

張宇祥，夏利娟，李琳依

(1.上海交通大學(xué) 海洋工程國家重點實驗室，上海 200240；2.高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海 200240)

0 引 言

輸流管道系統(tǒng)廣泛應(yīng)用于船舶海洋工程領(lǐng)域，其振動特性對于保證船舶正常工作而言是不可忽視的關(guān)鍵部分，由液體壓力脈動等交變載荷引發(fā)的管道振動，極有可能導(dǎo)致管道疲勞破壞，影響船舶設(shè)備的安全運行。對于石油鉆井平臺、LNG-FPSO（浮式天然氣生產(chǎn)儲卸設(shè)備）等結(jié)構(gòu)物，其管道系統(tǒng)中存在氣液兩相流，相較單相流其運動過程更為復(fù)雜，產(chǎn)生的流體激振力會導(dǎo)致管道顯著的振動。因此，對典型管道結(jié)構(gòu)中兩相流誘導(dǎo)振動特性進(jìn)行深入研究對于工程應(yīng)用十分必要。

目前國內(nèi)外針對氣液兩相流誘導(dǎo)振動主要研究管道外橫向兩相流和軸向兩相流，對管道內(nèi)流的研究相對較少。張濤等計算了不同含氣率下氣液兩相流對管道內(nèi)壁產(chǎn)生的流場作用力，表明當(dāng)含氣率大于50%時，長圓管內(nèi)兩相流流型轉(zhuǎn)變?yōu)閷訝盍?，流場作用力隨之減小。Riverin等測量了多種兩相流流型下U 型管道結(jié)構(gòu)的流場激勵力，表明兩相流激勵力由窄帶和周期性分量組成，且激勵力的主頻隨著兩相流折算流速的增加而成比例增加。張紅艷等基于流固耦合方法研究了管道振動峰值與不同氣液組分參數(shù)的段塞流的關(guān)系。馬曉旭等利用ADINA 對水平管內(nèi)氣液兩相流誘導(dǎo)振動進(jìn)行數(shù)值仿真，研究了兩相流流型、體積含氣率和分相折算速度對激振力的影響，總結(jié)了水平管道多相流流激振動的機理。蔡標(biāo)華等基于流固耦合算法，應(yīng)用有限體積和有限元方法研究了注水管道的流激振動特性，對振動產(chǎn)生的噪聲進(jìn)行了數(shù)值模擬。

氣液兩相流誘導(dǎo)管道振動由于流動的復(fù)雜性以及流固耦合的特殊性，對計算資源有較高的要求，數(shù)值仿真的相關(guān)研究相對較少。因此，本文使用Ansys Fluent 與Ansys Mechanical 進(jìn)行豎直彎管內(nèi)的氣液兩相流誘導(dǎo)振動分析，主要研究豎直彎管的轉(zhuǎn)彎角度對氣液兩相流流型變化以及管道動態(tài)響應(yīng)的影響，為涉及氣液兩相流的船舶管道系統(tǒng)設(shè)計提供參考。

1 兩相流管道數(shù)學(xué)模型

1.1 流固耦合控制方程

流固耦合問題的控制方程包括流體域、固體域、流固耦合交界面3 個部分。

對于流體域，其遵循的守恒定律包括質(zhì)量守恒定律、動量守恒定律以及能量守恒定律。在研究氣液兩相流誘導(dǎo)管道振動時，不考慮氣液兩相的相變過程，即在整個流體域內(nèi)不發(fā)生能量傳遞，因此僅給出粘性不可壓縮流體質(zhì)量和動量的控制方程，即Navier-Stokes 方程。

質(zhì)量控制方程：

動量控制方程：

式中：為流體速度矢量；ρ為流體密度；為流體壓力；為體積力矢量；μ為流體動力粘度。

氣液兩相流在運動過程中產(chǎn)生的激勵力作用在管道內(nèi)壁，并誘導(dǎo)管道振動。固體域即管道在該激勵力作用下的結(jié)構(gòu)動力學(xué)基本方程如下：

式中：為管道結(jié)構(gòu)質(zhì)量矩陣；為阻尼矩陣；為剛度矩陣；{} 為結(jié)構(gòu)運動的位移矢量；{˙}為速度矢量；{¨} 為加速度矢量；{}為作用在管道結(jié)構(gòu)上的外界動態(tài)載荷，即流體激勵力。

在流固耦合交界面處遵循最基本的物理守恒原則，在不考慮熱效應(yīng)的情況下，應(yīng)滿足流體域與固體域應(yīng)力 τ 和位移的守恒，即滿足以下方程：

式中：d和d分別為流體和管道結(jié)構(gòu)的位移；τ和τ分別為流體和管道結(jié)構(gòu)的剪切應(yīng)力。式(4)為位移協(xié)調(diào)方程，在流固耦合計算中用于確定耦合交界面處流體節(jié)點的位移；式(5)為應(yīng)力守恒方程，可由力平衡方程通過積分將流體域的均布壓力轉(zhuǎn)化為集中力施加到管道結(jié)構(gòu)節(jié)點上。

目前流固耦合問題的主流求解方法為分離解法，即無需求解耦合控制方程，而是按照規(guī)定好的順序在求解器中分別求解流體域和固體域的控制方程，并在流固耦合交界面處實現(xiàn)計算數(shù)據(jù)的互相傳遞，待當(dāng)前時間步計算達(dá)到收斂，再進(jìn)行下一時間步的計算，依次而行，直至得到最終結(jié)果。對于氣液兩相流誘導(dǎo)管道振動這一典型的流固耦合振動問題，本文采取分離解法中的單向流固耦合分析方法。

1.2 兩相流相關(guān)參數(shù)

氣液兩相流的流動形態(tài)有別于單相流動，其兩相交界面的形狀不斷變化，因此引入描述氣液兩相流流動特征的相關(guān)參數(shù)。

表1 氣液兩相流的相關(guān)參數(shù)Tab.1 Related parameters of gas-liquid two-phase flow

上述參數(shù)間存在以下關(guān)系：

2 兩相流管道數(shù)值模型

2.1 兩相流管道系統(tǒng)參數(shù)

選取豎直平面內(nèi)的彎管作為研究對象，將彎管轉(zhuǎn)彎角度記作 α，分別取 α的值為90°，105°，120°，135°，4 種不同布置形態(tài)管道的物理模型及相關(guān)尺寸參數(shù)如圖1 所示，尺寸參數(shù)參考文獻(xiàn)[7]中的水平彎管。本文所研究豎直彎管的外徑為30 mm，壁厚為2.3 mm，其在軸方向上的跨度為305 mm，在軸方向上的跨度為400 mm，彎頭半徑均為50 mm。管道兩端和為簡支約束，端為空氣和水的混合流入口，端為混合流的出口，選取管壁外側(cè)彎頭中間一點作為仿真計算的監(jiān)測點。

圖1 不同轉(zhuǎn)彎角度的豎直彎管示意圖Fig.1 The diagrammatic drawing of vertical elbows with different turning angel α

豎直彎管及管道內(nèi)流的物性參數(shù)如表2 所示。

表2 豎直彎管及管道內(nèi)流的物性參數(shù)Tab.2 Physical parameters of vertical elbows and internal flow of pipes

2.2 數(shù)值仿真方法和驗證

進(jìn)行氣液兩相流仿真與管道瞬態(tài)響應(yīng)的耦合計算，分別使用掃掠網(wǎng)格和四面體網(wǎng)格進(jìn)行結(jié)構(gòu)和流體的網(wǎng)格剖分。在進(jìn)行數(shù)值仿真之前，為獲取本文建立的數(shù)值模型求解精度受網(wǎng)格密度的影響，以轉(zhuǎn)彎角度α為90°的豎直彎管為例進(jìn)行網(wǎng)格無關(guān)性驗證。對該豎直彎管的網(wǎng)格剖分采用3 種依次加密的方案，在相同的數(shù)值仿真工況下進(jìn)行計算。

數(shù)值仿真工況設(shè)置如下：

1）入口采用速度入口邊界，氣相速度為2 m/s，液相速度為0.8 m/s，氣相的入口截面為與管道橫截面圓同心且直徑為內(nèi)徑一半（即12.7 mm）的圓截面，其體積分?jǐn)?shù)的值為1.0，液相的入口截面為氣相入口圓周和管道內(nèi)徑圓周之間的圓環(huán)截面，其體積分?jǐn)?shù)的值為1.0，圖2 為氣液兩相入口截面示意圖。可計算得到仿真入口處氣相折算速度u為0.5 m/s，液相折算速度u為0.6 m/s，體積含氣率 β為0.455。

圖2 氣液兩相入口示意圖Fig.2 The diagrammatic drawing gas-liquid two-phase flow inlet

出口采用壓力出口邊界；施加沿軸負(fù)方向的重力加速度，其值取為9.81 m/s；管道兩端為簡支約束，限制其在空間3 個自由度的平動運動；管道和流體相互接觸的壁面設(shè)置為流固耦合交界面。

2）采用標(biāo)準(zhǔn)-ε模型作為流體計算的湍流模型，該模型求解穩(wěn)定性強，應(yīng)用廣泛；采用幾何重構(gòu)的VOF 模型來進(jìn)行氣液兩相流相界面的追蹤，VOF 模型實現(xiàn)氣液運動界面追蹤的方法有多種，而幾何重構(gòu)（Geo-Reconstruct）是大多數(shù)瞬態(tài)VOF 仿真推薦的方法，是目前最精確的界面跟蹤方法；采用基于壓力耦合的求解器，壓力與速度耦合方式采用PISO 算法（Pressure-Implicit with Splitting of Operators），該算法是CFD 計算主流分離解法中比較適合瞬態(tài)不可壓縮流求解的方法。

在3 種網(wǎng)格劃分方案下，得到管道結(jié)構(gòu)在彎頭監(jiān)測點的位移響應(yīng)、速度響應(yīng)以及加速度響應(yīng)，并將三者的均方根值作為網(wǎng)格無關(guān)性驗證的衡量參數(shù)。位移響應(yīng)、速度響應(yīng)和加速度響應(yīng)的均方根值分別依照式(9)～式(11)進(jìn)行計算，90°豎直彎管的網(wǎng)格無關(guān)性驗證結(jié)果如表3 所示。

表3 轉(zhuǎn)彎角度 α為90°豎直彎管的網(wǎng)格無關(guān)性驗證Tab.3 Verification of grid independence of vertical elbows withα=90°

根據(jù)表3，采用3 種不同的網(wǎng)格剖分方案，管道結(jié)構(gòu)在兩相流激勵力作用下的位移響應(yīng)和速度響應(yīng)的均方根值差別不大，相對誤差均落在13%以內(nèi)。網(wǎng)格方案2 和網(wǎng)格方案3 在位移響應(yīng)這一計算結(jié)果上相差僅有1%～2%左右，即當(dāng)前網(wǎng)格密度下數(shù)值計算已達(dá)收斂。綜合考慮數(shù)值計算精度和計算資源，最終采用網(wǎng)格剖分方案2，在該方案下流體域和固體域的有限元模型如圖3 所示。對于其他轉(zhuǎn)彎角度的豎直彎管均采用了相同方法進(jìn)行網(wǎng)格無關(guān)性驗證來選擇相對合適的網(wǎng)格密度。

圖3 90°豎直彎管有限元模型Fig.3 Finite element model of vertical elbows with α=90°

3 氣液兩相流管道仿真結(jié)果

3.1 不同轉(zhuǎn)彎角度下管道內(nèi)氣液兩相流的流型分析

對于氣液兩相流誘導(dǎo)振動而言，流型是重要的影響因素，通常情況下氣相和液相的速度以及管道實際的布置形式?jīng)Q定了管道內(nèi)氣液兩相的分布以及特定流型的產(chǎn)生和發(fā)展。為了更直觀地觀測管道內(nèi)氣液兩相流的流型特征，提取不同時刻下管道中剖面上液相的相體積分?jǐn)?shù)云圖，如圖4 所示。

可以看出，在氣液兩相折算速度確定的情況下，相界面的波動與管道布置形式即豎直彎管的轉(zhuǎn)彎角度密切相關(guān)。仿真工況的入口邊界為環(huán)形流邊界，由于氣相進(jìn)入管道的實際速度相對較大，是液相進(jìn)入管道實際速度的2.5 倍，在初始時刻產(chǎn)生較大的兩相間沖擊力，導(dǎo)致相界面產(chǎn)生波動變化。以90°彎管在0.3 s時刻的相體積分?jǐn)?shù)云圖（見圖4（a））為例，管道豎直管段中產(chǎn)生了多段頭部較鈍、尾部扁平形如子彈的氣泡，由于氣相速度較快，在通過彎頭區(qū)域時快速擠壓液相，形成一段較為完整的氣泡，氣泡在管壁之間則形成很薄的一層液膜。這種氣液兩相的流動特征在其他轉(zhuǎn)彎角度的豎直彎管中同樣清晰可見，不同的是隨著轉(zhuǎn)彎角度 α的增大，入口邊界處的環(huán)形流快速向不穩(wěn)定分層流轉(zhuǎn)變，產(chǎn)生這種差異的主要原因是氣液兩相密度差導(dǎo)致二者所受重力不同，由于轉(zhuǎn)彎角度增大，相同時間內(nèi)氣液兩相通過管道的方向跨距減小，在重力的作用下氣相上浮，液相則向下聚集在管壁下方。

圖4 不同轉(zhuǎn)彎角度豎直彎管在不同時刻下管道中剖面液相體積分?jǐn)?shù)云圖Fig.4 Contour plots of water volume fraction in the middle section of vertical elbows with different turning angles at different times

氣液兩相在完全流經(jīng)豎直彎管后呈現(xiàn)的流型特征也較為統(tǒng)一，在水平管段隨著時間推移在重力的作用下逐漸形成波動相界面的分層流；在豎直管段，隨著轉(zhuǎn)彎角度 α的增大，管道在傾斜方向的實際長度增大，氣液兩相的發(fā)展更為充分，液相對氣相的擠壓作用更強烈，因此轉(zhuǎn)彎角度較大的管道更易產(chǎn)生不連續(xù)的氣泡，且液相相對而言占據(jù)著較大的管道空間，當(dāng)轉(zhuǎn)彎角度增大到135°時，氣液兩相在整個管段均相對連續(xù)。

3.2 不同轉(zhuǎn)彎角度下管道的動態(tài)響應(yīng)分析

對于不同轉(zhuǎn)彎角度下的豎直彎管，選取管壁彎頭區(qū)域中間一點作為監(jiān)測點，提取其振動響應(yīng)時域數(shù)據(jù)并得到時域分布曲線。不同轉(zhuǎn)角豎直彎管的振動響應(yīng)時域分布曲線相類似，以90°豎直彎管為例，其振動位移響應(yīng)的時域曲線如圖5 所示。

圖5 90°豎直彎管彎頭監(jiān)測點處的振動位移時域分布曲線Fig.5 Curve of the displacement and time at the monitoring point of vertical elbows with α=90°

可知，氣液兩相流在豎直彎管中誘導(dǎo)產(chǎn)生的振動具有高度不穩(wěn)定特性，在彎頭區(qū)域會產(chǎn)生較為劇烈的振動，且主要表現(xiàn)為垂向（方向）振動。

計算不同轉(zhuǎn)彎角度豎直彎管監(jiān)測點處的振動位移、振動速度以及振動加速度的均方根值，結(jié)果如表4所示。圖6 給出了管道振動響應(yīng)均方根值隨管道轉(zhuǎn)彎角度變化的趨勢。

根據(jù)表4 和圖6 的仿真計算結(jié)果可以得出，隨著轉(zhuǎn)彎角度的增大，豎直彎管的振動響應(yīng)即位移、速度和加速度響應(yīng)幅值的均方根值均呈現(xiàn)出先增大后減小的趨勢，這也與氣液兩相流在不同轉(zhuǎn)彎角度豎直彎管中的流型分布特征相吻合。

表4 不同轉(zhuǎn)彎角度豎直彎管的振動響應(yīng)均方根值Tab.4 Root mean square value of vibration response of vertical elbows with different turning angles

圖6 管道振動響應(yīng)隨轉(zhuǎn)彎角度變化曲線Fig.6 Curve of the flow-induced vibration response and time at the monitoring point

4 結(jié) 語

本文基于流固耦合方法針對氣液兩相流在豎直彎管內(nèi)的誘導(dǎo)振動進(jìn)行數(shù)值仿真和研究，探討管道布置形式即管道轉(zhuǎn)彎角度對于氣液兩相流流型以及管道振動響應(yīng)的影響，可以得到以下結(jié)論：

1）通過CFD 仿真以及結(jié)構(gòu)瞬態(tài)響應(yīng)仿真的耦合計算對氣液兩相管道內(nèi)流誘導(dǎo)振動進(jìn)行模擬是可行的，流固耦合計算可以實現(xiàn)對兩相流的流型進(jìn)行判斷以及對結(jié)構(gòu)的振動響應(yīng)進(jìn)行前期預(yù)判。

2）豎直彎管轉(zhuǎn)彎角度增大幅度較小時，如從90°增大到105°和120°，氣液兩相流相界面的波動會更為劇烈，在傾斜管段更易產(chǎn)生不連續(xù)的氣泡，在彎頭處氣液兩相之間的擠壓沖擊作用更加顯著；當(dāng)轉(zhuǎn)彎角度繼續(xù)增大，氣液兩相在整個管段均相對連續(xù)，彎頭處的相間擠壓沖擊作用相對平緩。

3）氣液兩相流在豎直彎管中誘導(dǎo)產(chǎn)生的振動具有高度隨機性和不穩(wěn)定性，隨著轉(zhuǎn)彎角度的增大，管道振動響應(yīng)均方根值均呈現(xiàn)出先增大后減小的趨勢，這與氣液兩相流的流型轉(zhuǎn)化密切相關(guān)。