王桂軍，閆雪飛，金良安，苑志江

(海軍大連艦艇學(xué)院 航海系，遼寧 大連 116018)

水面艦船在航行過(guò)程中會(huì)在其尾部一定區(qū)域內(nèi)留下一條長(zhǎng)達(dá)數(shù)千米的具有特殊性質(zhì)的尾流，其中的氣泡尾流由于可以被尾流自導(dǎo)魚(yú)雷探測(cè)到，具有重要的軍事和民用價(jià)值，是國(guó)內(nèi)外相關(guān)領(lǐng)域的一個(gè)研究熱點(diǎn)［1］。氣泡尾流可以按離艦船尾部的距離長(zhǎng)短劃分為近程尾流和遠(yuǎn)程尾流，其中關(guān)于遠(yuǎn)程靜尾流中氣泡動(dòng)力學(xué)特性的研究屢見(jiàn)不鮮［1］，而專門(mén)針對(duì)近程湍流尾流的研究報(bào)道卻很少見(jiàn)諸文獻(xiàn)。

目前對(duì)艦船尾流的氣泡特性研究還主要考慮氣體的傳質(zhì)因素，而很少有系統(tǒng)考慮聚并影響的。Stewart采用大渦模擬技術(shù)(LES)和拉格朗日粒子跟蹤方法(LPD)對(duì)尾流環(huán)境下的氣泡動(dòng)力學(xué)進(jìn)行了研究，重點(diǎn)分析了氣體傳質(zhì)與擴(kuò)散對(duì)氣泡群尺寸演變的影響［2］。陳圣濤采用氣液兩相流模型分別對(duì)近程和遠(yuǎn)程氣泡尾流含氣量的時(shí)間分布特征做了比較全面的數(shù)值模擬和實(shí)驗(yàn)研究，指出進(jìn)程氣泡尾流含氣量以乘冪形式衰減，而遠(yuǎn)程尾流基本以指數(shù)形式衰減，然而陳的方法過(guò)于復(fù)雜，體現(xiàn)在占有大量運(yùn)算內(nèi)存和時(shí)間，一個(gè)算例需要在集成計(jì)算機(jī)上運(yùn)算3 400個(gè)CPU小時(shí)［3］。朱東華［4］曾以MINER［5］的初始值生成方法為依據(jù)，對(duì)尾流的氣泡數(shù)密度分布進(jìn)行了三維數(shù)值模擬，基于波爾茲曼型氣泡輸運(yùn)方程和雷諾平均方程計(jì)算出了尾流中不同尺寸的氣泡數(shù)密度分布特征，結(jié)果表明氣泡尾流主要位于10～200 um之間。劉慧開(kāi)建立了海洋表層氣泡運(yùn)動(dòng)和半徑變化的數(shù)學(xué)模型，發(fā)現(xiàn)小于臨界半徑(74.1 um)的氣泡會(huì)迅速溶解，而大于臨界半徑的氣泡會(huì)上升到水面破裂，由此得出海水中的氣泡主要以50～74 um之間居多［6］。Carrica采用氣液兩相流技術(shù)對(duì)艦船尾流中不同位置處的氣含量及氣泡數(shù)密度分布進(jìn)行了數(shù)值模擬［7］，但沒(méi)有對(duì)聚并這一重要因素進(jìn)行必要的深入分析，而且Carrica的算法同樣要占大量存儲(chǔ)和運(yùn)算時(shí)間(200CPU h)，普通電腦上根本無(wú)法實(shí)現(xiàn)，難以達(dá)到實(shí)時(shí)應(yīng)用的要求。

綜上可知，當(dāng)前關(guān)于氣泡尾流的特征研究還存在很多缺陷，對(duì)氣泡聚并因素的認(rèn)識(shí)不足。1946年，美國(guó)國(guó)防研究委員會(huì)第六局用聲吶測(cè)量了以航速15 kn航行的驅(qū)逐艦產(chǎn)生的尾流氣泡分布規(guī)律，發(fā)現(xiàn)尾流中氣泡直徑為0.08～1.07 mm的氣泡數(shù)密度高達(dá)5.98×106/m3，如此高的數(shù)密度會(huì)不可避免導(dǎo)致聚并的發(fā)生［8］。此外，Smirnov在對(duì)對(duì)非動(dòng)力船只的近程尾流的研究中，也著重分析了氣泡聚并在近程尾流中發(fā)生的可能性，提出了氣泡聚并的重要作用，但并未對(duì)此進(jìn)行深入研究［9］。這里在對(duì)BPBE方程修正的基礎(chǔ)上，首次將其運(yùn)用于對(duì)氣泡尾流數(shù)密度和氣含量的研究上，重點(diǎn)探究了聚并因素的作用，且模型具有簡(jiǎn)單，可操作性強(qiáng)的特點(diǎn)，還可節(jié)省大量的計(jì)算時(shí)間(完整模型對(duì)3 000 m長(zhǎng)度尾流的推演時(shí)間不超過(guò)4 s)。

1 模型的構(gòu)建過(guò)程

1.1 模型構(gòu)建的思想與基礎(chǔ)

1.1.1 基本思想

傳統(tǒng)對(duì)氣泡尾流的數(shù)值模擬將尾流速度與氣泡運(yùn)動(dòng)耦合在一起進(jìn)行求解，大大的增加了模型求解的復(fù)雜性和難度，無(wú)法滿足軍事應(yīng)用的實(shí)時(shí)性要求。鑒于此，采用BPBE方程進(jìn)行相關(guān)求解，可將尾流速度的求解和氣泡特性的計(jì)算分開(kāi)而來(lái)，只要有一組準(zhǔn)確的速度場(chǎng)，即可模擬其中的氣泡動(dòng)力學(xué)特性，實(shí)踐表明，這種思想是可行的，并且可大大的節(jié)省時(shí)空復(fù)雜度。

1.1.2 基礎(chǔ)模型

氣泡群模型是尾流流場(chǎng)的氣泡群聚并模型構(gòu)建的重要基礎(chǔ)，一個(gè)典型的波爾茲曼型氣泡群平衡方程［10］:

式中:n(z，d，t)為與空間(z)、氣泡尺寸(d)、時(shí)間(t)有關(guān)的單位體積氣泡數(shù)密度函數(shù);u是對(duì)應(yīng)的(z)方向的氣泡速度;d(z，d)為與空間(z)、尺寸(d)有關(guān)的表示氣泡尺寸因氣體傳質(zhì)而引起變化的函數(shù);與氣泡聚并和破裂有關(guān)的源項(xiàng)S(z，d，t)可以如下表示:

其中，

分別表示體積小于Vi的氣泡與體積為(Vi-Vj)的氣泡聚并為體積為Vi的氣泡速率，氣泡Vi與其它氣泡聚并而導(dǎo)致氣泡Vi消失的速率，體積為2Vi的氣泡分裂產(chǎn)生兩個(gè)氣泡Vi的速率，氣泡Vi本身分裂消失的速率。

與普通的波爾茲曼型輸運(yùn)方程［4］相比，氣泡群平衡方程的優(yōu)勢(shì)在于不僅包含傳質(zhì)與擴(kuò)散項(xiàng)，還包括聚并與分裂項(xiàng)，因此對(duì)于氣泡群尺寸的分布與變化描述將更加精確，考慮到實(shí)際情況，尾流氣泡尺寸太小幾乎不會(huì)分裂，因此文中將忽略分裂的影響。下面將基于氣泡尾流模型對(duì)方程(1)進(jìn)行適當(dāng)?shù)幕?jiǎn)和離散處理，使其可以用于對(duì)艦船尾流的數(shù)值模擬。

1.2 模型建立的相關(guān)求解

1.2.1 模型的化簡(jiǎn)和離散

首先建立尾流場(chǎng)二維坐標(biāo)系，與航向平行的方向?yàn)閤，垂直水面向下的方向?yàn)閥。僅考慮氣泡數(shù)密度在軸向x和豎直方向y的變化，空間項(xiàng)可以表示:

注意到以下兩個(gè)關(guān)系式:

則有:

對(duì)于方程(1)中的擴(kuò)散與傳質(zhì)，有:

由式(7)可知，時(shí)間項(xiàng)可以被消去，因此式(9)可以化簡(jiǎn)為:

其中，ψ(d)表征擴(kuò)散與傳質(zhì)，是關(guān)于直徑d的一元函數(shù)。

聯(lián)立式(1)、(4)、(7)、(8)、(10)，且考慮到氣泡尺寸小到可以忽略分裂項(xiàng)，最終得到如下模型:

氣泡群平衡方程最終化為只與坐標(biāo)x和氣泡尺寸d有關(guān)的二維偏微分方程，在求解前需要對(duì)軸坐標(biāo)x和氣泡尺寸d進(jìn)行離散化處理，采用差分法離散后的數(shù)值模型:

其中，n(i，j)表示x=iΔx處尺寸為d=jΔd的單位體積氣泡數(shù)密度;u(i，j)表示x方向尺寸為d=jΔd的氣泡速度。下面分別給出聚并模型Ω與傳質(zhì)模型ψ的求解。

1.2.2 聚并模型及其修正

體積分別為Vi和Vj的氣泡集合聚并速率Ωc等于各自密度乘以碰撞概率乘以聚并概率］11］:

其中，碰撞概率θij等于湍流誘導(dǎo)的碰撞概率θijT加上浮力誘導(dǎo)的碰撞概率θijB。上述模型只適用于氣泡數(shù)密度小的實(shí)驗(yàn)室水箱或鼓泡塔，而對(duì)氣泡數(shù)密度巨大的尾流則無(wú)能為力，因?yàn)閚inj后其結(jié)果之大完全可以超乎想象，導(dǎo)致碰撞的發(fā)生會(huì)遠(yuǎn)遠(yuǎn)大于實(shí)際存在的氣泡數(shù)，這是與實(shí)際明顯相悖的，為此，將其修改為:

湍流誘導(dǎo)的碰撞概率:

浮力誘導(dǎo)的碰撞概率:

聚并概率:

其中接觸時(shí)間［11］:

聚并時(shí)間［11］:

上式思想來(lái)源于液膜排水模型，其中初始膜厚度h0=5×10-4，破裂液膜厚度hf=1×10-8。兩氣泡的特征尺寸dij:

1.2.3 傳質(zhì)模型

應(yīng)用田恒斗［1］建立的艦船尾流氣泡上浮傳質(zhì)模型:

其中，CA為海水中的氣體質(zhì)量濃度，CI為氣液界面處的氣體質(zhì)量濃度。

當(dāng)深度h固定時(shí)，對(duì)式(23)化簡(jiǎn)并對(duì)d求偏導(dǎo)即可得到ψ(d)的一元函數(shù)表達(dá)式。

2 計(jì)算與討論

求解尾流場(chǎng)的流動(dòng)參量的方法可采用MINER［5］提出的方法，在給定船型參數(shù)和附件參數(shù)(如船長(zhǎng)、船寬、吃水、船速、螺旋槳盤(pán)面尺寸、槳軸位置等)的情況下，可近似得到尾流空間各點(diǎn)的初始速度、湍動(dòng)能、湍動(dòng)能耗散率。根據(jù)理論分析和測(cè)量經(jīng)驗(yàn)，湍流尾流的速度發(fā)展變化可以用以下乘冪曲線近似表達(dá)［13］:

圖1 模擬的尾流速度沿x方向的變化Fig.1 The change of simulated wake speed along x

其中，u0為初始尾流速度，k為待定系數(shù)(文中設(shè)為5.0)。使用上述方程擬合的尾流速度如圖1所示。由圖可知擬合曲線光滑、連續(xù)、穩(wěn)定，與實(shí)際情形較為接近。

計(jì)算實(shí)例取文獻(xiàn)［3］給出的美國(guó)某實(shí)船模型，其中各參數(shù)具體:船長(zhǎng)為86.9 m，船寬為14.6 m，吃水3.6 m，船速為12節(jié)，螺旋槳直徑2.44 m。為了進(jìn)行數(shù)值計(jì)算，氣泡尺寸需要按照大小進(jìn)行分組，采用文獻(xiàn)［7］的分組方法，如表1所示。

這樣共得到15個(gè)相互耦合的差分方程組。由于要求解氣泡群尺寸沿x方向的分布，又由于氣泡的速度可以用流場(chǎng)速度U近似，因此有u(x，d)=U。模型將在普通臺(tái)式機(jī)上運(yùn)行，完整模型的運(yùn)行時(shí)間不超過(guò)4 s，相對(duì)以往文獻(xiàn)的數(shù)值模擬方法大大提高了計(jì)算效率，可以達(dá)到實(shí)時(shí)監(jiān)測(cè)的要求。

表1 各組分氣泡數(shù)密度的初始化Tab.1 Initialization of the bubble-group BND

2.1 沿尾流流向的氣含量分布

首先對(duì)沿尾流流向的氣含量分布特征進(jìn)行求解分析。用于計(jì)算氣含量的公式:

這里，由于差分步長(zhǎng)Δx設(shè)為1 m，因此單位體積VL=1 m3，尾流中總的氣含量表達(dá)式:

圖2為只考慮傳質(zhì)、聚并及二者綜合情形下的尾流中氣含量沿尾流流向的分布?？梢钥闯鋈N情形下的氣含量總體分布均以乘冪形式進(jìn)行衰減。相對(duì)于傳質(zhì)因素，聚并作用在近尾流區(qū)(小于1 800 m)更加明顯，這是由于近尾流區(qū)的船殼、螺旋槳攪動(dòng)、湍流等因素作用非常強(qiáng)烈，大大增加了氣泡碰撞的幾率，促進(jìn)了聚并的發(fā)生，在聚并的影響下，氣含量衰減曲線更顯陡峭。當(dāng)進(jìn)入遠(yuǎn)程趨于靜止的尾流區(qū)(大于1 800 m)后，氣泡間的影響減弱，聚并作用不再明顯，這時(shí)傳質(zhì)因素開(kāi)始發(fā)揮主要影響，使得氣泡含量繼續(xù)以直線方式緩慢減少。

將完整模型的計(jì)算結(jié)果與文獻(xiàn)［3］的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行對(duì)比，在比較前，需要將長(zhǎng)度單位換算為時(shí)間(min)，近似為:

其中，L為尾流長(zhǎng)度，U為船速。

此外，由于初始值不同，文中的結(jié)果大小會(huì)與文獻(xiàn)［3］有所不同(相差約0.1)，因此，為了在同一幅圖中比較，還需各自進(jìn)行歸一化，歸一化后將文中計(jì)算結(jié)果與文獻(xiàn)［3］進(jìn)行對(duì)比，如圖3所示。由于文獻(xiàn)［3］的研究分為近程尾流(小于1 800 m)和遠(yuǎn)場(chǎng)尾流(大于1 800 m)，因此測(cè)量曲線沒(méi)有連續(xù)，但仍可將其結(jié)果與文中氣含量計(jì)算值的相應(yīng)區(qū)間進(jìn)行比較。由圖可知文中計(jì)算結(jié)果與文獻(xiàn)［3］中氣含量實(shí)驗(yàn)結(jié)果曲線的走向吻合較好。

圖4為文中計(jì)算的BND沿尾流流向分布與文獻(xiàn)［4］的計(jì)算結(jié)果的對(duì)比。可以看出，雖然增加了考慮聚并因素，仍與文獻(xiàn)［4］吻合較好，進(jìn)一步表明了文中計(jì)算模型的有效性，同時(shí)說(shuō)明聚并對(duì)BND的變化影響不大。

圖2 氣含量沿x方向的變化Fig.2 The change of air content with x

圖3 氣含量隨時(shí)間的變化Fig.3 The change of air content with time

圖4 氣泡數(shù)密度沿x方向的變化Fig.4 The change of bubble BND with x

2.2 氣泡尾流尺寸的相對(duì)分布

下面對(duì)氣泡尾流不同尺寸的相對(duì)氣含量及BND分布特征進(jìn)行計(jì)算分析。圖5是在分別考慮傳質(zhì)因素、聚并因素及二者綜合情形下，位于尾流充分發(fā)展的遠(yuǎn)場(chǎng)尾流區(qū)(x=3 000 m)的不同尺寸氣泡BND的相對(duì)值分布。由圖可知只在聚并作用下的氣泡數(shù)密度值與氣泡尺寸成反比，這主要是由于大氣泡在湍流誘導(dǎo)下的碰撞概率更高，因此聚并消失的速率更大，密度下降也最快。而在傳質(zhì)作用下，小氣泡由于迅速溶解而消失，大氣泡由于迅速上浮至水面破裂而消失，致使尺寸約為70～80 um的氣泡數(shù)密度最大。根據(jù)以上結(jié)論，綜合情形下的相對(duì)BND分布與只有傳質(zhì)影響下的相對(duì)BND分布區(qū)別不大。

圖6是位于尾流充分發(fā)展的遠(yuǎn)場(chǎng)尾流區(qū)(x=3 000 m)的不同尺寸氣泡氣含量的相對(duì)值分布。由圖可知尺寸200 um的氣泡氣含量最大，這與文獻(xiàn)［7］針對(duì)氣泡氣含量相對(duì)分布的研究結(jié)論是基本一致的。與圖5相比，可知對(duì)于只有聚并因素的情形，雖然小氣泡的密度較大，但是由于體積較小，因此氣含量就會(huì)相對(duì)變小(由表1也可明顯看出)。同時(shí)可以看出雖然只在傳質(zhì)作用與綜合情形下的氣泡相對(duì)BND分布比較接近，但是二者氣含量的相對(duì)分布卻有較大差異，這說(shuō)明聚并在氣泡尾流的氣含量分布與變化中發(fā)揮著主要作用，不能忽視。

圖5 不同尺寸氣泡相對(duì)數(shù)密度Fig.5 Relative BND of bubble with different dimensions

圖6 不同尺寸氣泡相對(duì)氣含量Fig.6 Relative air content of bubble with different dimensions

3 結(jié)語(yǔ)

氣泡尾流以其強(qiáng)大的制導(dǎo)作用而一直被各國(guó)軍界作為研究的重點(diǎn)來(lái)對(duì)待，鑒于以往對(duì)氣泡尾流特征參數(shù)的分布研究方面存在的不足，基于BPBE方程提出了可以同時(shí)將傳質(zhì)與聚并作用包含在內(nèi)的用于氣泡尾流特征分布研究的新思想，并使用數(shù)值方法對(duì)模型進(jìn)行了仿真求解，分別對(duì)傳質(zhì)因素、聚并因素以及二者綜合因素情形下氣泡氣含量以及BND的特征進(jìn)行了分析研究，結(jié)果表明聚并在近尾流區(qū)作用強(qiáng)烈，而在遠(yuǎn)程尾流區(qū)傳質(zhì)起主要作用。聚并對(duì)BND的分布變化影響不大，但對(duì)氣含量的影響則作用明顯，不可忽視。模型大大提高了計(jì)算效率，達(dá)到了軍事應(yīng)用的實(shí)時(shí)性需求，且結(jié)果與相關(guān)文獻(xiàn)及實(shí)驗(yàn)數(shù)據(jù)吻合較好，表明了該模型的有效性。

符合說(shuō)明:ξ為湍動(dòng)能耗散率;g為引力常數(shù)，10 m/s2;σ為表面張力，0.074 N/m;ρL為海水密度，1 200 kg/m3;Patm為標(biāo)準(zhǔn)大氣壓，105Pa;R為氣泡半徑，m;ρg為氣體密度，1 kg/m3;vb為氣泡上浮速度，m/s;DAB為氣體擴(kuò)散系數(shù)，1.8×10-9m2/s;v為流體粘度，0.001 kg/ms。

［1］ 田恒斗，金良安，遲衛(wèi).船舶遠(yuǎn)程尾流場(chǎng)中氣泡上浮運(yùn)動(dòng)模型［J］.船舶力學(xué)，2010，14(11):1195-1201.(TIAN Heng-dou，JIN Liang-an，CHI Wei.Rising process model of bubble in the far wake of ship［J］.Journal of Ship Mechanics，2010，14(11):1195-1201.(in Chinese))

［2］ Stewart M B，Miner E W.Bubble Dynamics in a Turbulent Ship Wake［R］.NRL Memorandmm Report，6055，1987.

［3］ 陳圣濤，王慧麗，王運(yùn)鷹，等.艦船氣泡尾流特性的數(shù)值模擬和實(shí)驗(yàn)研究［J］.船舶力學(xué)，2012，16(4):342-349.(CHEN Sheng-tao，WANG Hui-li，WANG Yun-ying，et al.Numerical simulation and experimental study of a ship's bubble wake［J］.Journal of Ship Mechanics，2012，16(4):342-349.(in Chinese))

［4］ 朱東華，張曉暉，顧建農(nóng)，等.艦船尾流及其氣泡數(shù)密度分布的數(shù)值計(jì)算［J］.兵工學(xué)報(bào)，2011，32(3):315-320.(ZHU Dong-hua，ZHANG Xiao-hui，GU Jian-nong，et al.Numerical calculation of ship wake and its bubble number density distribution［J］.Acta Armamentarii，2011，32(3):315-320.(in Chinese))

［5］ Miner E W，Ramberg S E.Method for Approximating the Initial Data Plane for Surface Wake Simulation［R］.US:RL-MR-6376，1988.

［6］ 劉慧開(kāi)，張勸華，楊 立.海洋表層氣泡運(yùn)動(dòng)規(guī)律的研究［J］.海洋科學(xué)，2009，33(1):34-38.(LIU Hui-kai，ZHANG Quan-hua，YANG Li.The study on the motion of bubbles near the ocean surface［J］.Marine Sciences，2009，33(1):34-38.(in Chinese))

［7］ Carrica P M，Drew D，Bonetto F，et al.A poly disperse model for bubbly two-phase flow around a surface ship［J］.International Journal of Multiphase Flow，1999，25(2):257-305.

［8］ 高 江，張靜遠(yuǎn)，楊 力.艦船氣泡尾流特性研究現(xiàn)狀［J］.艦船科學(xué)技術(shù)，2008，30(4):27-32.(GAO Jiang，ZHANG Jing-yuan，YANG Li.The present situation of research on shipwake characteristic［J］.Ship Science And Technology，2008，30(4):27-32.(in Chinese))

［9］ Smimov A，Celik I，Shi S.LES of bubble dynamics in wake flows［J］.Computers＆ Fluids，2005(34):351-373.

［10］Ryszard P W，Moniuk W l，Bielski P，et al.Modelling of the coalescence/redispersion processes in bubble columns［J］.Chemical Engineering Science，2001，56:6157-6164.

［11］Prince M J，Blanch H W.Bubble coalescence and break-up in air-sparged bubble columns［J］.A.I.Ch.E.Journal，1990，36(10):1485-1499.

［12］徐麥容，劉成云.水中浮升氣泡的半徑和速度變化［J］.大學(xué)物理，2008，27(11):14-17.(XU Mai-rong，LIU Cheng-yun.The change of radius and velocity of the rising bubble in water［J］.College Physics，2008，27(11):14-17.(in Chinese))

［13］James K E，John R D，Brian A M.The radar image of the turbulent wake generated by a moving ship［J］.London Research and Development，1993，755:343-346.