傅慧萍

(上海交通大學海洋工程國家重點實驗室，上海200240)

微氣泡減阻的隱含機理仍然是研究熱點［1］。Legner認為近壁流體可以被認為是空隙率可變的各向同性混合物，宏觀量和湍流輸運的變化導致摩擦阻力的下降［2］。Madavan等的結論類似，二者的分析都驗證了空隙率的存在導致了粘性底層的增厚從而降低了摩擦阻力，并得出“當氣泡停留在湍流邊界層(TBL)的過渡層時最有效”的結論［3］。Kunz等采用歐拉雙流體計算了BDR，得到了與減阻率密切相關的近壁氣泡群的演化［4］。多相流在高精度計算方面的進展使氣液相互作用的細觀研究成為可能。Ferrante等采用歐拉－拉格朗日模型揭示了氣泡與湍流在近壁區(qū)的相互作用細節(jié)［5］。Lu等采用直接數(shù)值模擬揭示了:當大的可變形氣泡出現(xiàn)在過渡區(qū)，會與流向渦產(chǎn)生強烈的相互作用［6］。Mohanarangam等基于ANSYS CFX中的多尺寸組模型對二維平板BDR進行了數(shù)值模擬，并全面借鑒了試驗模型的試驗參數(shù)及豐富的試驗數(shù)據(jù)進行了數(shù)學模型的校驗［7］。丁力等采用混合物模型對二維簡化船型進行了模擬，得到了噴氣參數(shù)對船舶阻力的影響規(guī)律，在噴氣量的相似換算方面取得了一些成果［8］。陳顯文等同樣采用混合物模型對某回轉體及SUBOFF進行了模擬。由于采用的也是混合物模型，無法考慮微氣泡運動過程中的變形、破裂和聚合，也忽略了氣泡尺度的影響，而且數(shù)值計算結果沒有試驗數(shù)據(jù)進行支撐與校驗［9］。王炳亮等分別采用Euler雙流體模型(均一氣泡直徑)和混合物模型對三維平板及某散貨船模型進行了模擬［10］。與前兩者的區(qū)別是參照了Madavan等［11］的試驗參數(shù)進行數(shù)值建模。但并沒有采用后者的試驗結果進行對比，而是采用了一系列的平板阻力經(jīng)驗公式進行阻力評估獲得CFD結果的參考數(shù)值進行誤差分析。我國目前所擁有的微氣泡減阻技術離實用還有一段距離。根據(jù)現(xiàn)代測試手段獲取的試驗數(shù)據(jù)，建立更加符合物理本質(zhì)的微氣泡減阻模型;發(fā)展高效、高精度的高雷諾數(shù)網(wǎng)格生成技術與數(shù)值分析方法，用以分析微氣泡減阻的流動特性和機理，是微氣泡減阻數(shù)值模擬的發(fā)展趨勢。

1 研究對象

本文的研究對象為平板，為與試驗對比，依文獻［11］進行建模。圖1給出了平板底部通氣的二維數(shù)值模擬示意圖。計算域總長0.61 m，高0.057 m。沿長度方向，將平板分為3段:wall-1、wall-2、wall-3，長度分別為:0.178 m、0.178 m、0.254 m，其中wall-2為噴氣入口。

圖1 計算域Fig.1 Computational field

圖2為三維建模中的底板部分，采用對稱邊界條件。陰影部分為所考察的平板，寬為0.051 m，計算域寬度為0.250 m。

圖2 三維建模Fig.2 3D model

網(wǎng)格劃分的邊界層最小網(wǎng)格尺度為Δyp按下式計算:

式中:y+為當?shù)乩字Z數(shù)，L為特征長度，此處取計算域總長。由于本文的研究對象為5～500 μm的微氣泡，故邊界層網(wǎng)格要求較高，當?shù)乩字Z數(shù)要求為個位數(shù)的量級。計算中采用第1層網(wǎng)格高度是0.001 5 mm，圖3所示為邊界層網(wǎng)格劃分。

圖3 邊界層網(wǎng)格Fig.3 Mesh in the boundary

邊界條件的設置如圖1所示，左端和wall-2分別為水流和噴氣入口，均設置為速度入口;右端為流動出口，設置為outflow類型;wall-1、wall-3設置為無滑移壁面;計算域最大高度處設置為水流速度入口條件。當計算無通氣平板的摩擦阻力時，wall-2設為壁面條件即可。

2 多相流模型

2.1 混合物模型

混合物模型用來模擬2種流體混合之后的兩相流，它將各相設置為相互貫穿的連續(xù)體，從而只對混合物求解共同的控制方程。連續(xù)方程:

式中:Um為混合相速度，ρm為混合相密度。動量方程:

式中:μm為混合相黏度，p為壓力，j為多相流的相數(shù)(j=1，2)，αj、ρj、Udr，j分別為第j相的體積分數(shù)、密度和漂移速度。

氣相體積分數(shù)方程:

式中:Sg為氣相的源項。相對速度(或滑移速度)，是指次要相(第2相)相對于主要相(第1相)的速度:

設任意相(j相)的質(zhì)量分數(shù)定義為

則漂移速度與相對速度通過下式關聯(lián):

FLUENT里的混合物模型采用代數(shù)滑移公式，基本假設是為了能對相對速度進行代數(shù)表達，要求在一個較短的空間尺度上達到局部相間平衡。本文采用Manninen相對速度形式。

2.2 相群平衡模型

相群平衡模型是在歐拉雙流體框架下開啟，它將第二相(空氣)進行尺寸分組，對每一組的體積分數(shù)進行求解。假設將氣泡分成N個尺寸組，各組氣泡之間存在聚并和破碎作用。那么第i組氣泡的連續(xù)性方程為

式中:Si為第i組氣泡由于聚并和破碎產(chǎn)生的源項，本文采用Luo模型;fi為第i組的體積分率:

假設ni和vi分別為第i組氣泡的數(shù)密度和氣泡體積，則有如下關系:

假設氣泡是球形的，各組氣泡體積比為

若擬在5～500 μm內(nèi)對氣泡直徑進行分組(N =5)，則q=4.985?？梢缘玫礁鞒叽绶纸M如下:D0= 500 μm、D1=158 μm、D2=50 μm、D3=16 μm、D4= 5 μm。索特平均直徑Dm是聯(lián)系流場與微觀氣泡群的一個變量，即歐拉雙流體模型和相群平衡模型通過它聯(lián)系起來。該變量與各氣泡分組的直徑及其數(shù)密度有關，定義如下:

至于湍流的模擬，本文3種多相流求解方法均是對混合相進行湍流方程的求解，所采用的湍流模型為標準k-ε兩方程模型。

3 計算結果與分析

3.1 重力的影響

為考察重力的影響，計算時先不考慮重力，待計算收斂后，再考慮重力至再度收斂。圖4給出了相同水流速度和噴氣速度下(Uw=9.6 m/s，Ua= 0.165 2 m/s)，氣泡直徑D分別采用5、50和500 μm得到的平板阻力系數(shù)收斂歷程。圖中曲線基本呈2段式:前一段不考慮重力，后一段考慮重力后，阻力明顯呈上升趨勢。這是由于考慮重力后，氣泡在水中受到向上的浮力作用。對于底部通氣方式，通氣方向與浮力方向一致，浮力作用使氣泡離開壁面，氣體的潤滑作用減小，阻力增加。從前后2段曲線的提升程度看，小氣泡的阻力增加幅度小于大氣泡，大氣泡阻力系數(shù)大。

圖4 重力的影響Fig.4 Effect of gravity on Cd

3.2 通氣方式比較

由上節(jié)的研究可知重力對大氣泡的影響顯著，為此選擇 500 μm泡徑，在 Uw=9.6 m/s，Ua= 0.165 2 m/s時，將二維平板模型旋轉180°，使底部通氣變?yōu)轫敳客?，探討通氣方式對BDR的影響。計算時仍然先不考慮重力，待第1階段收斂后再考慮。圖5給出了2種通氣方式的比較。由圖5可見:重力對2種通氣方式的影響是完全相反的。對于底部通氣，如前節(jié)所述，重力使阻力增大。對于頂部通氣，通氣方向與浮力方向相反，浮力是使氣泡趨向壁面并停留在邊界層內(nèi)的，氣體的潤滑作用得以保持所以使阻力減小。對于船舶微氣泡減阻，應該是適用于頂部通氣方式，重力的影響是有利于減阻的。

圖5 通氣方式的影響Fig.5 Effect of injection method on Cd

3.3 減阻機理研究

3.3.1 氣層厚度與減阻率

設Cd為通氣情形下平板“wall-3”段的無量綱摩擦阻力系數(shù);Cd0為不通氣情形下的“wall-3”段無量綱摩阻系數(shù);則減阻率DR由下式定義:

圖6給出了Uw=9.6 m/s時，均一氣泡直徑D= 50 μm時的采用混合物模型得到的減阻率與氣層厚度的關系。減阻率先是隨噴氣速度增大而增大，當噴氣速度達到一定數(shù)值后，減阻率不增反降，即減阻率曲線存在一個峰值。氣層厚度卻不同，隨噴氣速度增加而單調(diào)增大。因此減阻率存在一個峰值，對應著一個最佳氣層厚度或最佳噴氣速度。這是因為噴氣速度過小時，微氣泡還不足以形成足夠的氣層，甚至類似于表面粗糙度，產(chǎn)生負增長;如果噴氣速度過大，減阻效果也會下降，這是因為垂直向上的氣流對水平方向的來流造成的擾動影響了流動的光滑性，導致了減阻效果的下降。

圖6 減阻率與氣層厚度的關系Fig.6 Relation between DR and gas layer thickness

3.3.2 體積分率以及氣泡數(shù)密度

考慮到采用均一氣泡直徑的混合物模型得到的減阻率與試驗值對比偏高，而PBM模型可以對氣泡進行尺寸分組，更好地模擬了利用金屬燒結板形成微氣泡的氣泡生成方式，本節(jié)將采用PBM模型對平板底部通氣BDR進行數(shù)值模擬。利用PBM模型進行計算時，邊界條件中通氣口wall-2的各分組的體積分率均設為0.2，即假設通氣口處，各種直徑的氣泡對于氣相分數(shù)的貢獻一致。氣泡尺寸分組是PBM模型特有的，可以通過對體積分率以及氣泡數(shù)密度的分析來觀察通氣口下游平板wall-3處各組氣泡的分布及演化。本節(jié)給出了二維平板底部通氣當Uw=9.6 m/s，Ua=0.055 1 m/s時，采用5～500 μm的尺寸分組(分10組)的PBM模擬結果。各組氣泡的直徑和平板wall-3處的體積分率如表1所示，對應的曲線見圖7。由表1及圖7可以看出:在平板wall-3處，直徑最大和最小的氣泡都有明顯減少，而中間尺寸的氣泡體積分率有所增大，各組分的體積分率不再相等。

表1 各組氣泡的直徑和wall-3處體積分率Table 1 fion wall-3 to bubble size Di

圖7 wall-3處各組氣泡對應的體積分率Fig.7 fion wall-3 to bubble size Di

去流向各個截面上不同尺度的氣泡數(shù)密度及體積分率變化如圖8所示。由圖8可見，大部分尺度的氣泡數(shù)密度和體積分率沿流向都呈減小趨勢，只有大氣泡組呈增大趨勢，兩者的變化趨勢是一致的。即小氣泡在流動方向上發(fā)生聚合變成大氣泡，大氣泡發(fā)生破裂，但由于聚合的影響大于破裂，總體上大尺度分組的數(shù)密度和體積分率增大。

圖8 氣泡數(shù)密度與體積分率沿流向的變化Fig.8 niand fialong x axis

3.4 三維效應

3.3 節(jié)采用PBM模型對二維平板底部通氣BDR進行模擬，得到的減阻率與試驗值對比仍然偏高。為此，本節(jié)開展基于PBM的三維數(shù)值模擬。圖9給出了相同水流速度和噴氣速度下(Uw=9.6 m/ s，Ua=0.055 1 m/s)，采用5～500 μm的尺寸分組(分10組)進行計算，得到的阻力系數(shù)收斂歷程。由圖9可見，三維較之于二維，阻力系數(shù)提升一個量級。多相流中α表示第二相即氣相的體積分數(shù):α= 1，表示純氣流;α=0，表示純水流;α為0～1的任意值時，表示氣水混合流。本節(jié)定義的氣層厚度δα表示氣相分數(shù)為α的等值面上的最大厚度。

圖9 三維效應對阻力系數(shù)的影響Fig.9 Effect of 3-D on Cd

圖10給出了當α分別為0.1、0.25、0.5時的等值面厚度分布。由圖10可見，由于氣泡的上浮運動，α等值面沿流向厚度逐漸增加，最大厚度出現(xiàn)在出口靠近中心線處;隨著α的增大，等值面最大厚度減小，外緣沿流向向中心線收縮。

圖10 三維BDR模擬Fig.10 3-D simulation of BDR

3.5 側壁通氣

側壁通氣BDR模擬必須是三維全流場模擬，由此也可進一步完善重力影響的研究。本節(jié)對平板側壁通氣BDR進行基于PBM的三維數(shù)值模擬。圖11給出了相同水流速度和噴氣速度下(Uw=9.6 m/s，Ua=0.055 1 m/s)，采用5～500 μm的尺寸分組(分10組)進行計算，得到的側壁通氣與底部通氣的Cd收斂歷程對比。由圖11可見，側壁通氣經(jīng)過了一個偽收斂平臺再下降至最后的收斂值。可能因為松弛因子太低，收斂很慢。

圖11 側壁通氣的阻力系數(shù)Fig.11 Cdwith plate vertical

圖12給出了當α分別等于0.1、0.25、0.5時的等值面厚度分布。由圖12可見，側壁通氣的氣層厚度要大于底部通氣的。在α=0.1的等值面上最大厚度為9.83 mm，為底部通氣(3.86 mm)的2.55倍。但此最大厚度已偏離出所考察的平板(wall-3)。

另外，側壁通氣的減阻效果低于底部通氣，這是因為:由于氣泡的上浮運動導致大部分區(qū)域氣相分數(shù)降低，甚至局部區(qū)域出現(xiàn)全沾濕狀態(tài)，即未被氣體潤滑。

圖12 側壁通氣BDR模擬Fig.12 Simulation of BDR with plate vertical

3.6 計算方法驗證

基于PBM方法，采用50～500 μm的5分組氣泡群，對Uw=9.6 m/s時的平板底部噴氣進行數(shù)值模擬。表2和圖13給出了數(shù)值結果與試驗值［3］的比較。

表2 計算值與試驗值的比較Table 2 Comparison between CFD and EFD

圖13 數(shù)值模擬校驗Fig.13 Validation of numerical method

數(shù)值模擬所得變化趨勢與實驗所得變化趨勢基本相同，而計算所得與試驗值相比偏高，這個誤差并不影響本文對微氣泡減阻的機理及尺度效應的定性研究。這可能是由于二維模擬所帶來的系統(tǒng)誤差所致。如果將二維模型變?yōu)槿S模型，并進一步增大氣泡直徑尺寸分組，可以預見，阻力系數(shù)的預報精度將提高，但隨之而來的是計算量增大和計算穩(wěn)定性問題。

4 結論

1)重力對大氣泡影響較大:考慮重力的影響后，氣泡就受到了水中的浮力，浮力使氣泡上升。越大的氣泡所受浮力也越大。

2)重力對通氣方式的影響:對于頂部通氣方式，重力使氣泡停留在壁面邊界層，從而使減阻率相較于底部通氣方式得到改善。

3)減阻率與氣層厚度之間存在一定的關系:噴氣速度增加，氣層厚度也增加，但減阻率在達到一個峰值后開始下降。

4)三維效應使減阻率更接近試驗值:故三維建模與計算是必要的。

5)側壁通氣的減阻效果不佳:重力效應使流場上下不對稱，氣泡沿流向向上漂移，考察對象局部區(qū)域全沾濕導致。

進一步的研究將拓展到三維船體底部通氣建模與計算。

［1］CECCIO S L.Friction drag reduction of external flows with bubble and gas injection［J］.Annual Review of Fluid Mechanics，2010，42(1):183-203.

［2］LEGNER H H.A simple model for gas bubble drag reduction［J］.Physics of Fluids，1984，27(12):2788-2790.

［3］MADAVAN N K，MERKLE C L，DEUTSCH S.Numerical investigations into the mechanisms of micro-bubble drag reduction［J］.Journal of Fluids Engineering，1985，107(3): 370-377.

［4］KUNZ R F，GIBELING H J，MAXEY M R，et al.Validation of two-fluid Eulerian CFD modeling for microbubble drag reduction across a wide range of Reynolds numbers［J］.Journal of Fluids Engineering，2007，129(1):66-79.

［5］FERRANTE A，ELGHOBASHI S.Reynolds number effect on drag reduction in a microbubble-laden spatially developing turbulent boundary layer［J］.Journal of Fluid Mechanics，2005，543:93-106.

［6］LU Jiacai，F(xiàn)ERNáNDEZ A，TRYGGVASON G.The effect of bubbles on the wall drag of a turbulent channel flow［J］.Physics of Fluids，2005，17(9):095102.

［7］MOHANARANGAM K，CHEUNG S C P，TU J Y，et al.Numerical simulation of micro-bubble drag reduction using population balance model［J］.Ocean Engineering，2009，36(11):863-872.

［8］丁力.微氣泡減阻噴氣參數(shù)換算關系研究［D］.武漢:武漢理工大學，2012:1-47.

DING Li.The conversion relation research about jet parameters of micro-bubbles drag reduction［D］.Wuhan:Wuhan University of Technology，2012:1-47.

［9］陳顯文.回轉體微氣泡減阻和噪聲的數(shù)值研究［D］.武漢:華中科技大學，2012:1-77.

CHEN Xianwen.Numerical simulation of axisymmetric body's drag reduction and noise by micro bubbles［D］.Wuhan: Huazhong University of Science and Technology，2012:1-77.

［10］王炳亮.船舶微氣泡減阻數(shù)值模擬及機理研究［D］.哈爾濱:哈爾濱工程大學，2012:1-70.

WANG Bingliang.Numerical simulation and mechanism research on drag reduction of ship by microbubbles［D］.Harbin:Harbin Engineering University，2012:1-70.

［11］MADAVAN N K，DEUTSCH S，MERKLE C L.Reduction of turbulent skin friction by microbubbles［J］.Physics of Fluids，1984，27(2):356-363.