低速船舶微氣泡減阻數(shù)值研究

趙曉杰，宗 智，王加夏，洪智超，4，胡俊明

（1.江蘇科技大學船舶與海洋工程學院，江蘇鎮(zhèn)江 212000；2.大連理工大學船舶工程學院，遼寧大連 116024；3.工業(yè)裝備結構分析國家重點實驗室，遼寧大連 116024；4.江蘇省船舶與海洋工程裝備技術創(chuàng)新中心，江蘇南通 226100）

0 引 言

航行中的船舶會受到阻力。為了實現(xiàn)減阻，研究者們提出了許多減阻技術[1]，如為了減少興波阻力，采用球鼻艏或優(yōu)化船體型線；為了降低摩擦阻力，在船舶表面設置超疏水涂層，或在船底噴入高分子聚合物、微氣泡等。其中，微氣泡減阻技術受到人們的廣泛關注，因為它的減阻效率較高并且無污染。對于低速船舶，摩擦阻力占主要部分，微氣泡減阻技術更具應用前景。

幾十年來，國內(nèi)外學者對微氣泡減阻進行了許多試驗研究，研究對象主要有平板、回轉體、船體等。McCormick 等[2]首先通過試驗發(fā)現(xiàn)了微氣泡減阻現(xiàn)象，然后研究學者又相繼對平板和回轉體的微氣泡減阻進行了試驗研究[3]。對于船舶微氣泡減阻，減阻情況更加復雜，因為減阻效果會受船型、噴氣形式、復雜工況等多個因素影響。Sayyaadi等[4]對一雙體船進行了試驗研究，發(fā)現(xiàn)噴氣量對減阻效果有顯著影響，并存在一個臨界噴氣量；王家楣等[5]采用多孔硅板生成微氣泡，對不同吃水的船模微氣泡減阻進行了試驗研究，發(fā)現(xiàn)船舶在淺吃水下減阻更好。

前人也對微氣泡減阻做了許多數(shù)值研究。Kunz 等[6]利用兩相歐拉模型研究了平板的微氣泡減阻，通過和試驗結果對比驗證了模型的準確性；傅慧萍[7]利用不同多相流模型對平板微氣泡減阻進行了數(shù)值模擬并對影響因素進行了分析；對于船舶微氣泡減阻，鄭曉偉等[8]利用相間滑移算法研究了二維船的微氣泡減阻，所得結果與試驗定性符合；吳浩與歐勇鵬[9]將船底簡化為平板并利用Mixture模型研究噴氣減阻，較好地模擬出了氣液混合流。可以看出，目前船舶微氣泡減阻模擬都對數(shù)值模型進行了一定簡化，幾何模型簡單，并且沒有充分考慮氣液間及氣泡間相互作用。

因此，為了進一步研究船舶微氣泡減阻規(guī)律，本文基于OpenFOAM 中兩相歐拉模型對低速散貨船的微氣泡減阻展開數(shù)值研究，考慮五種相間作用力以及氣泡的聚合和破碎，并考慮氣泡對液體湍流的影響。區(qū)別于二維平板，船舶具有不同阻力成分，本文研究氣泡對船舶摩擦阻力、粘壓阻力和總阻力的減阻率影響，同時研究噴氣量、氣泡直徑、航速及吃水等因素對減阻率的影響，分析船舶附近氣體體積分數(shù)、湍流粘度和氣泡直徑的分布等，為微氣泡減阻技術在船舶上的應用提供理論指導。

1 數(shù)值模型

1.1 兩相歐拉模型

在OpenFOAM 的兩相歐拉模型中，對氣液兩相（本文指水和空氣）分別進行建模，建立的質量和動量守恒方程[10]分別為

式中：i為相指數(shù)，這里i=l或g分別表示液體和氣體；αi表示i相體積分數(shù)；Ui表示i相流體速度；ρi表示i相流體密度；p為壓力；g為重力加速度；∑Fi為i相的相間作用力，用于考慮相間動量傳遞；νi表示i相運動粘度系數(shù)；νt,i表示i相湍流運動粘度系數(shù)。

1.2 相間作用力

為了充分考慮氣液間動量傳遞作用，本文考慮五種相間作用力，公式如下：

式中：Fd為相間阻力，F(xiàn)l為升力，F(xiàn)vm為虛質量力，F(xiàn)td為湍流離散力，F(xiàn)wl為壁面潤滑力。

相間阻力是氣泡在液體中運動時受到的阻力，與相對速度方向相反，計算公式為

式中：αg為氣體體積分數(shù)，dg為氣泡直徑，ρl為液體密度，Ul和Ug分別表示液體速度和氣體速度，Cd為相間阻力系數(shù)。本文采用Ishii和Zuber[11]提出的阻力模型，該模型考慮三種氣泡形態(tài)，對于球形氣泡，阻力系數(shù)計算公式為

對于橢球型氣泡，阻力系數(shù)計算公式為

式中，Re是氣泡雷諾數(shù)，Eo表示無量綱E?tv?s數(shù)，f(αg)表達式為

對于帽子形氣泡，阻力系數(shù)計算公式為

不同氣泡形狀區(qū)域劃分基于以下公式：

升力是由于液體速度梯度引起的剪切力，它是除阻力外最重要的作用力，可以影響氣體垂向分布，計算公式為

式中：Cl是升力系數(shù)。采用Tomiyama升力模型[12]，可以考慮不同直徑的氣泡，計算公式為

式中，Cvm為虛擬質量力系數(shù)，這里假設氣泡是球形，根據(jù)勢流理論取Cvm=0.5。

湍流耗散力為流體湍流脈動作用對氣泡產(chǎn)生的作用力，采用Bertodano模型[13]，計算公式為

壁面潤滑力是為了防止氣泡接觸壁面而產(chǎn)生的使氣泡遠離壁面的作用力，計算公式[14]為

式中，壁面潤滑力系數(shù)Cwl的計算公式為

式中：yw為氣泡離壁面的垂直距離；Cw1和Cw2是校準系數(shù)，根據(jù)文獻[14]推薦，取值分別為Cw1=-0.01 和Cw2=0.05，表示壁面潤滑力在距離壁面5倍直徑距離時變?yōu)?。

1.3 氣泡群體平衡模型

為了考慮氣泡直徑變化，采用基于界面面積密度方程的氣泡群體平衡模型，控制方程[15]為

式中：a表示單位體積氣泡界面面積密度；ψ表示氣泡形狀因子，忽略氣泡變形，ψ取值為1/36π；Sj表示氣泡聚并和破碎的源項，可以分成三部分：

其中，SRC和SWE分別表示隨機碰撞和尾跡夾帶作用，它們會引起氣泡聚并，表達式分別為

式中：ng為單位體積氣泡數(shù)密度，ut與ur分別為湍流速度和氣液相對速度。根據(jù)Ishii等[15]實驗研究，對于隨機碰撞，常數(shù)取值CRC=0.04和C=3。對于尾跡夾帶，常數(shù)取值CWE=0.04和αmax=0.75。STI為湍流作用，它導致氣泡發(fā)生破碎，表達式為

式中，CTI為湍流作用參數(shù)，Wecr為臨界韋伯數(shù)，取值分別為CTI=0.06和Wecr=6.0。

1.4 湍流模型

對于湍流模型，本文采用Lahey[16]提出的改進k-ε模型。為了考慮氣泡影響，在控制方程中加入源項，并增加額外湍流粘度來計算氣泡導致的湍流增加。湍流模型的控制方程為

式中，湍流粘度系數(shù)νt的計算公式為

源項Φkl與Φεl的計算公式為

式中，常數(shù)Cp取值為0.25，其它參數(shù)為Cμ=0.09，Cε,1=1.44，Cε,2=1.92，σk=1.0，σε=1.3。

2 計算設置

2.1 幾何模型

本文研究的是95 000 DWT 散貨船，縮尺比為1:50，實船垂線間長為231 m,型寬為38 m，型深為20.7 m，設計吃水為12.5 m，排水體積為91 893.0 m3，濕表面積為12 543.0 m2，設計航速為14.5 kn，船體型線見圖1。在船模17站處設置長方形噴氣板，長400 mm，寬90 mm，如圖2所示。

圖1 散貨船橫剖線圖Fig.1 Sectional curves of the bulk carrier

圖2 船舶噴氣板示意圖Fig.2 Sketch of the jet plate on the ship model

2.2 計算域和網(wǎng)格劃分

計算區(qū)域如圖3 所示，不考慮自由面影響，采用疊模模型，將自由面當作對稱滑移壁面。為了減少計算量，采用一半船體進行計算。速度入口距離船體1倍船長，出口邊界距離船體3倍船長，底面和側面距離船體1倍船長。

圖3 計算區(qū)域Fig.3 Computational domain

網(wǎng)格劃分利用切割體網(wǎng)格，網(wǎng)格分布如圖4所示。最大網(wǎng)格尺寸為231 mm，在船體附近進行網(wǎng)格加密，最小網(wǎng)格尺寸約為7.2 mm，船體表面設置三層邊界層，第一層網(wǎng)格厚度滿足y+>30。

圖4 網(wǎng)格分布Fig.4 Mesh distribution

2.3 邊界條件和計算方法

對于邊界條件，船體對液體速度采用無滑移壁面條件，對于氣泡速度采用滑移邊界。入口為速度入口，出口為壓力出口。對于噴氣面，本文假設采用微米多孔板產(chǎn)生微氣泡，為了匹配實際物理條件，假定孔隙率為0.3，與Kunz等[6]數(shù)值研究中條件一致。對于其他邊界，設置為對稱邊界條件。

控制方程采用有限體積法進行離散。對于空間離散，梯度項和拉普拉斯項采用高斯線性格式（Gauss linear），對流項采用限制高斯線性格式（Gauss limited linear），在相輸運控制方程中，采用Gauss vanLeer離散格式?？臻g離散為二階精度。對于時間離散，采用一階精度的隱式歐拉格式。

控制方程求解采用隱式速度壓力耦合算法（PISO），在OpenFOAM[17]中利用PIMPLE 循環(huán)進行設置，外迭代次數(shù)設為1，壓力矯正次數(shù)為2。時間步采用自適應時間步，最大庫朗數(shù)設為0.9，從而保證計算收斂。

2.4 計算工況設置

為了研究航速、噴氣量、氣泡直徑和吃水等因素對船舶微氣泡減阻的影響，設置了不同計算工況?？紤]三種航速，分別為0.946 m/s、1.055 m/s 和1.164 m/s，其中1.055 m/s 為設計航速。考慮三種吃水，分別為0.138 m、0.25 m 和0.29 m，其中0.25 m 為模型船的設計吃水?？紤]三種氣泡直徑，分別為200 μm、500 μm和800 μm。計算工況設置總結見表1。

表1 計算工況總結Tab.1 Summary of calculation conditions

采用無量綱噴氣量Cv作為參數(shù)，表達式為Cv=Qg/(Qg+Ql)，式中Qg是氣體噴氣流量，Ql是邊界層內(nèi)液體流量。Ql計算公式為Ql=U( )δ-δ?b，式中U為船舶航速，b為船寬，δ為邊界層厚度，δ?為邊界層排擠厚度，表達式為δ?=0.125δ，根據(jù)邊界層理論[18]，δ計算公式為δ=0.37x/Re0.2。無量綱噴氣量Cv范圍為0～0.05，間隔為0.01。

本文研究三種阻力成分的減阻率，計算公式為ΔRi=(Ri-R′i)/Ri，式中Ri表示未噴氣時阻力，R′i表示噴氣后阻力，下標i取f、p和t時分別表示摩擦阻力、粘壓阻力和總阻力。

3 計算結果與分析

3.1 計算模型驗證

為了驗證本文模型，采用兩種算例：一是計算不噴氣下船舶摩擦阻力，并和ITTC 推薦的摩擦阻力計算公式進行對比；二是計算噴氣情況下平板微氣泡減阻，并和前人試驗數(shù)據(jù)進行對比。其中ITTC推薦的摩擦阻力系數(shù)公式為

表2 表示船舶在三種速度下摩擦阻力和通過ITTC 推薦公式計算得到的摩擦阻力結果對比，可以看出計算誤差在2%之內(nèi)，說明本文的數(shù)值模型可以準確計算船舶阻力。

表2 船舶計算摩擦阻力和ITTC公式計算摩擦阻力對比Tab.2 Comparison of calculated ship resistance and the values from ITTC formula

平板微氣泡減阻試驗驗證采用Merkle 等[3]的實驗，通過多孔介質板產(chǎn)生微氣泡實現(xiàn)減阻，采用實驗工況，考慮流速為9.6 m/s和14.2 m/s，本文氣泡群體平衡模型需要先給定平均氣泡直徑，根據(jù)實驗觀測[3]以及Mohanarangam 等[19]在相同工況計算所得氣泡直徑分布，兩種流速下平均氣泡直徑分別取600 μm和400 μm。圖5和圖6分別表示流速為9.6 m/s和14.2 m/s時的數(shù)值結果，橫坐標為噴氣量，縱坐標為平板摩擦阻力減阻率，可以看出數(shù)值結果和試驗值吻合良好，說明本文微氣泡減阻數(shù)值模型是可行的和準確的。

圖5 速度為9.6 m/s時數(shù)值結果Fig.5 Numerical results at a speed of 9.6 m/s

圖6 速度為14.2 m/s時數(shù)值結果Fig.6 Numerical results at a speed of 14.2 m/s

3.2 不同氣泡直徑下船舶微氣泡減阻率

氣泡直徑是影響微氣泡減阻率的重要因素。圖7（a）～（c）表示在不同氣泡直徑下船舶摩擦阻力、粘壓阻力和總阻力的微氣泡減阻率，對應計算工況1～3。從圖中可以看出噴氣情況下，摩擦阻力、粘壓阻力和總阻力都實現(xiàn)了減阻，并且隨著噴氣量的增大，減阻率也都隨之增大，說明噴氣量是影響船舶微氣泡減阻率的決定因素。另外，摩擦阻力和總阻力的減阻率增大速度一開始很快，但隨后會變得緩慢，直到達到最大值后幾乎不再改變，說明對于摩擦阻力和總阻力，存在一個最佳臨界噴氣量，這和Sayyaadi等[4]試驗研究所得的結論相一致。對于粘壓阻力，減阻率增加存在波動，這是因為微氣泡對流場存在擾動、局部壓力出現(xiàn)不穩(wěn)定波動造成的，但總體上減阻率單調(diào)增加。

圖7 不同氣泡直徑下不同阻力成分微氣泡減阻率Fig.7 Microbubble drag reduction rates of different drag components for different bubble sizes

對于不同氣泡直徑，從圖7（a）～（c）可以看出，較小的氣泡可以獲得較大的摩擦阻力、粘壓阻力和總阻力減阻率，說明小氣泡可以更有效地減阻，這和Song等[20]三維回轉體的微氣泡減阻試驗結論相一致。由于船舶是低速散貨船，摩擦阻力占比較大（約80%），而粘壓阻力占比較?。s20%），盡管粘壓阻力減阻率存在不穩(wěn)定波動，但是摩擦阻力減阻率和總阻力減阻率的變化趨勢還是彼此接近。

3.3 不同氣泡直徑下船舶附近氣體分布

圖8（a）～（c）給出了無量綱噴氣量為0.05時不同氣泡直徑下船舶附近氣體體積分數(shù)分布?？梢钥闯鲈谧钚馀荩╠=200 μm）時，船底氣體分布均勻，計算得到表面氣體體積分數(shù)均值為0.327，此時減阻率也較大；在最大氣泡（d=800 μm）時，靠近噴氣板處氣體覆蓋較多，隨著遠離噴氣板氣體逐漸減少，導致船底氣體分布不均勻，計算得到表面氣體體積分數(shù)均值為0.308，此時減阻率也較小，說明近壁面氣體體積分數(shù)分布是影響船舶微氣泡減阻的關鍵因素。

圖8 不同氣泡直徑下船舶附近氣體體積分數(shù)分布Fig.8 Distribution of air volume fraction around the ship for different bubble sizes

不同氣泡直徑導致船體附近氣體體積分數(shù)分布不同，原因是不同氣泡直徑的氣液相互作用力不同。由公式（4）可以看出，在相同氣體體積分數(shù)時，氣泡直徑越小，氣液間相間阻力越大，使得小氣泡從噴氣板噴出后難以遠離壁面，更容易沿著液體流動方向延伸至船尾，從而造成氣體分布均勻，實現(xiàn)較高減阻率；而對于大氣泡，氣液相間阻力較小，使得大氣泡難以沿著液體流動方向向船尾延伸，從而導致氣體在噴氣口附近聚集，船尾附近氣體體積分數(shù)較低，導致減阻率也較低。

3.4 不同氣泡直徑下湍流運動粘度分布

為了進一步分析不同氣泡直徑下船舶微氣泡減阻率不同的原因，圖9（a）～（d）給出了不噴氣時及噴氣時（Cv=0.05）不同氣泡直徑的船舶附近湍流運動粘度分布?？梢钥闯?，噴氣后船底液體湍流運動粘度都會變大，說明氣泡對液體湍流具有增強作用，這和Mohanarangam 等研究得到的結論相一致。根據(jù)公式（23），可以得到在相同氣體體積分數(shù)下，氣泡直徑越大湍流運動粘度也越大的結論。圖9也表明，800 μm、500 μm 和200 μm 直徑氣泡產(chǎn)生的湍流運動粘度依次減小，因此這也是小氣泡可以更有效減阻的原因之一。

圖9 不同氣泡直徑下船舶附近湍流運動粘度分布Fig.9 Distribution of turbulent kinematic viscosity around the ship for different bubble sizes

3.5 船舶附近氣泡直徑變化

圖10（a）～（c）為無量綱噴氣量為0.05 時不同初始氣泡直徑下船舶附近氣泡直徑分布。從圖中可以看出，在船底沿著流動方向氣泡直徑都有所增大，這和Merkle 等[3]平板微氣泡減阻實驗的氣泡變化趨勢相同，原因是氣泡沿流向會發(fā)生隨機碰撞，再加上相鄰氣泡尾跡夾帶的作用，氣泡會發(fā)生聚并，導致氣泡直徑變大。在船舶尾部附近，由于存在較大的幾何形狀突變，船尾湍流作用劇烈，就會容易導致氣泡破碎，氣泡直徑又會變小。

圖10 Cv=0.05時不同初始氣泡直徑下的船舶附近氣泡直徑分布Fig.10 Distribution of bubble size around the ship for different initial bubble sizes when Cv=0.05

另外，從圖中可以看出不同初始氣泡直徑下的氣泡直徑變化情況不同。對于較小初始直徑氣泡如200 μm直徑氣泡，在噴氣板不遠處就開始發(fā)生聚并，最大氣泡直徑可達約260 μm，氣泡直徑增加約30%，大氣泡在船底分布也多；而對于較大初始直徑氣泡如800 μm 直徑氣泡，在距離噴氣板較遠處才開始發(fā)生聚并，最大氣泡直徑約為910 μm，氣泡直徑增加約13%，增大的氣泡僅分布在船尾附近?？梢娸^小氣泡的聚并作用更強，原因在于不同氣泡直徑導致的氣泡聚并和破碎作用不同，根據(jù)公式（18）～（19），氣泡隨機碰撞和尾跡夾帶作用與氣泡數(shù)密度有關，相同噴氣量下小氣泡的氣泡數(shù)密度更大，所以氣泡聚并作用更強。又根據(jù)3.4 節(jié)可得，小氣泡的湍流增強作用較小，所以湍流作用導致的氣泡破碎作用也較弱，較小氣泡的聚并作用更加劇烈。

3.6 航速與吃水對船舶微氣泡減阻的影響

圖11 為不同船舶航速下不同阻力成分的微氣泡減阻率，對應工況1、4 和5?？梢钥闯觯煌剿傧履Σ磷枇涂傋枇p阻率都隨噴氣量的增加而增加，到最后增加速度會變慢，減阻率都會接近一個恒定值；粘壓阻力減阻率也隨噴氣量的增加而增加，但也存在波動，這是因為氣泡導致流場局部壓力不穩(wěn)定造成的。此外，航速較高時船舶減阻率較高，說明此時船舶附近氣體體積分數(shù)較大。這是因為不同航速下氣液相互作用力不同，根據(jù)公式（4）可知流速較高時，氣泡和水的相對速度較大，此時氣液相間阻力較大，氣泡噴出后遠離壁面趨勢較小，氣泡在較強的氣液相互作用下更容易延伸至船尾，造成船底氣體體積分數(shù)較大且氣體分布均勻，因此船底獲得更大氣體覆蓋面積，從而獲得較大減阻率。

圖11 不同航速下不同阻力成分微氣泡減阻率Fig.11 Microbubble drag reduction rates of different drag components for different ship speeds

圖12 為不同船舶吃水下不同阻力成分的微氣泡減阻率，對應工況1、6 和7?？梢钥闯觯瑢τ谀Σ磷枇?、粘壓阻力和總阻力，小噴氣量時減阻率差別較小，隨著噴氣量增加，減阻率差別逐漸變大。另外，小吃水時減阻率較大，這和王家楣等[5]試驗的結論相一致，這是因為船舶減阻主要原因是船底摩擦阻力的減少，小吃水時船底摩擦阻力在整個船體總阻力中占比較大，更容易獲得較大的減阻率。

圖12 不同船舶吃水下不同阻力成分微氣泡減阻率Fig.12 Microbubble drag reduction rates of different drag components for different ship drafts

4 結 論

本文基于OpenFOAM 中兩相歐拉模型對三維低速散貨船的微氣泡減阻進行了數(shù)值研究，考慮五種相間作用力、氣泡聚合和破碎以及氣泡對液體湍流的影響，研究了氣泡對不同阻力成分減阻率的影響，研究了噴氣量、氣泡直徑、航速及吃水等因素的影響，分析了船舶附近氣體體積分數(shù)、湍流粘度和氣泡直徑的分布，得出的主要結論如下：

（1）氣泡可以同時減少低速船舶摩擦阻力、粘壓阻力和總阻力，并且隨著噴氣量增大，減阻率也都隨之增大，說明噴氣量是影響船舶微氣泡減阻率的決定因素。

（2）船舶附近氣體體積分數(shù)分布是影響船舶微氣泡減阻的關鍵因素。氣泡較小時船舶附近的平均氣體體積分數(shù)較大，且氣體分布更加均勻，減阻率較高。另外，較小氣泡產(chǎn)生的湍流運動粘度也較小，這也是小氣泡可以更有效減阻的原因之一，因此在工程應用中應盡量保證氣泡直徑足夠小。

（3）氣泡會沿著液體流動方向發(fā)生聚并，導致氣泡直徑增大。較小氣泡的聚并作用更劇烈，原因是相同噴氣量下它的氣泡數(shù)密度更大，隨機碰撞和尾跡夾帶作用更強，另外較小氣泡的湍流增強作用較小，導致氣泡破碎作用也較弱。

（4）較高航速和小吃水更有利于船舶微氣泡減阻。較高航速時，氣液相間阻力較大，船底氣體體積分數(shù)較大且氣體分布均勻，從而獲得較大減阻率。小吃水時，船底摩擦阻力在整個船體總阻力中占比較大，更容易獲得較大的減阻。