傅慧萍

(上海交通大學(xué)海洋工程國(guó)家重點(diǎn)實(shí)驗(yàn)室，上海200240)

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

1 研究對(duì)象

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

圖1 計(jì)算域Fig.1 Computational field

圖2為三維建模中的底板部分，采用對(duì)稱邊界條件。陰影部分為所考察的平板，寬為0.051 m，計(jì)算域?qū)挾葹?.250 m。

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

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

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

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

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

2 多相流模型

2.1 混合物模型

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

式中:Um為混合相速度，ρm為混合相密度。動(dòng)量方程:

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

氣相體積分?jǐn)?shù)方程:

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

設(shè)任意相(j相)的質(zhì)量分?jǐn)?shù)定義為

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

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

2.2 相群平衡模型

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

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

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

假設(shè)氣泡是球形的，各組氣泡體積比為

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

至于湍流的模擬，本文3種多相流求解方法均是對(duì)混合相進(jìn)行湍流方程的求解，所采用的湍流模型為標(biāo)準(zhǔn)k-ε兩方程模型。

3 計(jì)算結(jié)果與分析

3.1 重力的影響

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

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

3.2 通氣方式比較

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

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

3.3 減阻機(jī)理研究

3.3.1 氣層厚度與減阻率

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

圖6給出了Uw=9.6 m/s時(shí)，均一氣泡直徑D= 50 μm時(shí)的采用混合物模型得到的減阻率與氣層厚度的關(guān)系。減阻率先是隨噴氣速度增大而增大，當(dāng)噴氣速度達(dá)到一定數(shù)值后，減阻率不增反降，即減阻率曲線存在一個(gè)峰值。氣層厚度卻不同，隨噴氣速度增加而單調(diào)增大。因此減阻率存在一個(gè)峰值，對(duì)應(yīng)著一個(gè)最佳氣層厚度或最佳噴氣速度。這是因?yàn)閲姎馑俣冗^(guò)小時(shí)，微氣泡還不足以形成足夠的氣層，甚至類似于表面粗糙度，產(chǎn)生負(fù)增長(zhǎng);如果噴氣速度過(guò)大，減阻效果也會(huì)下降，這是因?yàn)榇怪毕蛏系臍饬鲗?duì)水平方向的來(lái)流造成的擾動(dòng)影響了流動(dòng)的光滑性，導(dǎo)致了減阻效果的下降。

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

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

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

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

圖7 wall-3處各組氣泡對(duì)應(yīng)的體積分率Fig.7 fion wall-3 to bubble size Di

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

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

3.4 三維效應(yīng)

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

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

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

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

3.5 側(cè)壁通氣

側(cè)壁通氣BDR模擬必須是三維全流場(chǎng)模擬，由此也可進(jìn)一步完善重力影響的研究。本節(jié)對(duì)平板側(cè)壁通氣BDR進(jìn)行基于PBM的三維數(shù)值模擬。圖11給出了相同水流速度和噴氣速度下(Uw=9.6 m/s，Ua=0.055 1 m/s)，采用5～500 μm的尺寸分組(分10組)進(jìn)行計(jì)算，得到的側(cè)壁通氣與底部通氣的Cd收斂歷程對(duì)比。由圖11可見，側(cè)壁通氣經(jīng)過(guò)了一個(gè)偽收斂平臺(tái)再下降至最后的收斂值?？赡芤?yàn)樗沙谝蜃犹停諗亢苈?/p>

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

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

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

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

3.6 計(jì)算方法驗(yàn)證

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

表2 計(jì)算值與試驗(yàn)值的比較Table 2 Comparison between CFD and EFD

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

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

4 結(jié)論

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

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

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

4)三維效應(yīng)使減阻率更接近試驗(yàn)值:故三維建模與計(jì)算是必要的。

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

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

［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ù)換算關(guān)系研究［D］.武漢:武漢理工大學(xué)，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］陳顯文.回轉(zhuǎn)體微氣泡減阻和噪聲的數(shù)值研究［D］.武漢:華中科技大學(xué)，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ù)值模擬及機(jī)理研究［D］.哈爾濱:哈爾濱工程大學(xué)，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.