金良安，閆雪飛，遲 衛(wèi)，陳明榮

（海軍大連艦艇學院 航海系，遼寧 大連 116018）

0 引 言

艦船尾流具有存活時間長，目標特征明顯和易被探測等特點，因此一直備受關注，是軍事領域的一個重要研究熱點[1]。實驗探測、理論計算和數(shù)值模擬是當前對尾流研究的三種主要方法，其中數(shù)值模擬由于具有無需設備、經(jīng)濟適用和節(jié)省人力物力等諸多優(yōu)點而倍受青睞。

隨著計算機性能的飛速提升，對尾流的數(shù)值模擬方法也在不斷地進步和發(fā)展。艦船尾流中可供探測的特征參數(shù)包括速度、湍動能、耗散率及氣泡尾流等，其中主要以氣泡尾流的研究為主。田恒斗[2]通過建立艦船遠程尾流場中單個氣泡上浮運動與傳質速率的數(shù)學耦合模型，分析研究了氣泡尾流的存活時間、氣泡數(shù)密度的衰減規(guī)律以及不同尺寸氣泡相對數(shù)密度的分布特征。Stewart[3]采用大渦模擬技術（LES）和拉格朗日粒子跟蹤方法對尾流環(huán)境下的氣泡動力學進行了研究，重點分析了氣體傳質與擴散對氣泡群尺寸演變的影響。朱東華[4]在艦船尾流相關特征參數(shù)求解的基礎上對尾流的氣泡數(shù)密度分布進行了三維數(shù)值模擬，該方法基于波爾茲曼型氣泡輸運方程和雷諾平均方程，計算出了尾流中不同尺寸的氣泡數(shù)密度分布特征，發(fā)現(xiàn)尾流中的氣泡尺寸主要以10～200 μm為主。其它針對艦船尾流氣泡分布特征的數(shù)值分析研究更是屢見不鮮[5-9]。

當前針對艦船尾流的數(shù)值計算方法主要以直接數(shù)值模擬（DNS）和求解復雜的氣液兩相流模型為主，要消耗巨大的存儲空間和計算時間，不是特定性能的巨型機根本無法實現(xiàn)，更不能達到實時性監(jiān)測的軍事斗爭需求，為此，本文首次利用LES和BPBE方程對流場速度和氣泡運動分開求解，大大地節(jié)省了時空復雜度，實現(xiàn)了艦船氣泡尾流的大范圍數(shù)值模擬（模擬長度為3 000 m）。

1 艦船尾流LES數(shù)值模型

1.1 模型建立的需求

針對當前氣泡尾流主要數(shù)值研究方法的復雜性以及要占用大量的運行時間和存儲空間的缺陷，將流體參量與氣泡運動分開再來進行求解，不僅可以滿足誤差要求，還可節(jié)省大量的時空復雜度，具有事半功倍的效果。且流體計算和氣泡運動方程都分別具有較好的移植性，便于用于其它兩相流領域的應用和求解。為此，先建立尾流LES控制方程，以計算BPBE方程所需的一組速度參數(shù)，是氣泡尾流求解的基礎。

1.2 模型建立的過程

首先建立尾流模型的三維直角坐標系，沿尾流速度方向為x軸，縱向為y軸，垂直水面向下為z軸，對應各坐標軸的速度分量分別為u、v、w。使用格子濾波函數(shù)對N-S方程進行濾波，濾波后得到完整的LES控制方程如下式（其中i，j=x，y，z，展開方式滿足求和約定）：

式中：ν是運動粘度系數(shù)，等式右邊出現(xiàn)不封閉項τij，即為亞格子尺度應力，需要構建模型予以封閉，常用的Smagorinsky渦粘模型如下所示：

式中：cs為Smagorinsky常數(shù)，其值在0.1～0.27之間，Δ是格子濾波的長度量綱。在給出初始邊界條件的情況下，采用數(shù)值離散方法即可對上述LES控制方程實現(xiàn)求解。

2 BPBE數(shù)值模型

2.1 模型構建的基礎

BPBE方程通常用于鼓泡塔等湍流流場的氣泡群運動特征的求解，因此可以移植過來應用于近尾流湍流流場的氣泡群求解，而且包含的影響因素較為全面，下面將給出尾流中BPBE方程的詳細構建過程。

一個典型的BPBE方程如下式所示[10]：

式中：n（z,d,t）為與空間 （z）、氣泡尺寸 （d）、時間 （t）有關的單位體積氣泡數(shù)密度函數(shù)，u是對應的 （z）方向的氣泡速度，d（z,d）為與空間（z）、尺寸（d）有關的表示氣泡尺寸因氣體傳質而引起變化的函數(shù)，此外，與氣泡聚并有關的源項S（z,d,t）可以如下表示：

其中，

分別表示體積小于Vi的氣泡與體積為（Vi- Vj）的氣泡聚并為體積為Vi的氣泡速率，氣泡Vi與其它氣泡聚并而導致氣泡Vi消失的速率。這里假設聚并過程沒有體積虧損[10]。

下面我們將基于艦船尾流實際情形對方程（3）進行適當?shù)幕喓碗x散處理，以方便程序更好更快地計算和求解。

2.2 模型建立的相關求解

僅考慮氣泡數(shù)密度在軸向x和豎直方向y的變化，空間項可以如下表示：

注意到以下兩個關系式：

則有：

對于方程（3）中的擴散與傳質，有：

由（9）式可知，時間項可以被消去，因此（11）式可以化簡為如下形式：

其中：ψ（d）表征擴散與傳質，是關于直徑d的一元函數(shù)。聯(lián)立各式且考慮到氣泡尺寸小到可以忽略分裂項，最終得到如下模型：

碰撞概率θij由下式給出[12]：

聚并概率由下式給出[12]：

其中： 初始膜厚度 h0＝5×10-4， 破裂液膜厚度 hf=1×10-8，d 為氣泡直徑，ρg為氣體密度，ρl為液態(tài)密度，u是運動粘度，ξ為湍動能耗散率，Patm為標準大氣壓，CA為海水中的氣體質量濃度，CI為氣液界面處的氣體質量濃度，具體公式見參考文獻[1]。n（x,d）表示x處尺寸為d的單位體積氣泡數(shù)密度。u（x,d）表示x方向尺寸為d的氣泡速度。氣泡群平衡方程最終化為只與坐標x和氣泡尺寸d有關的二維偏微分方程，在求解前需要對軸坐標x和氣泡尺寸d進行離散化處理。

3 計算與討論

首先采用前文提出的LES模擬方法對流場速度進行求解，在已知初始條件的情況下即可實現(xiàn)對流場速度的數(shù)值模擬。流場初始流動參量的方法可采用MINER[13]提出的方法，湍動能耗散率的解算方法將結合MINER與LES計算結果進行求解。得到速度場后（展向、垂向、流向尺寸為：5 m×1 m×3 000 m），將其代入到BPBE方程求解BND的分布。本文將以文獻[4]給出的美國某型驅逐艦為實船模型，各參數(shù)具體如下：艦長為155 m，艦寬為20 m，吃水6 m，排水量8 000 t，航速15 kns，螺旋槳直徑5 m。

圖（1）是x-y平面的位于流向位置（x=20 m）處的速度u等值線圖，從圖中可以看出速度在對稱面處的值最大，沿展向和垂向逐漸減小，且呈對稱分布。同時由圖形的形狀可以看出其分布與Kelvin尾流橫截面結構相似。圖（2）反映了某一位置處（y=0 m、z=1 m）速度u沿流向的變化，可以看出尾流平均速度非常小，在500 m以內以光滑的指數(shù)曲線迅速衰減，之后變化緩慢，與文獻[14]測量結果的變化規(guī)律基本一致，這表明了本文LES計算方法的準確性，雖然在大小上有一定的偏差，造成誤差的主要原因是算例中選取的目標船與文獻[12]的實測目標船2者船型結構參數(shù)的不同引起的。

圖1 速度u的等值面圖Fig．1 The isoline-section of velocity u

圖2 速度u沿流向x的變化Fig．2 The change of velocity u with x

3.2 尾流場BND計算與分析

為了進行數(shù)值計算，氣泡尺寸需要按照大小進行分組，根據(jù)理論分析和實測經(jīng)驗，分組方法將如表1所示，分別給出了各組氣泡的尺寸以及氣泡數(shù)密度及單位體積分數(shù)，這樣共得到15個相互耦合的差分方程組。模型將在普通臺式機上運行，運行時間不超過5 s，相對以往文獻的數(shù)值模擬方法大大提高了計算效率，可以達到實時監(jiān)測的要求。

dbi /μm ni (1/m3) εdi (%) dbi /μm ni (1/m3) εdi (%) dbi /μm ni (1/m3) εdi (%)10 21 35 53 77 2.61e+10 6.31e+9 2.97e+9 1.69e+9 1.04e+9 0.010 9 0.024 5 0.052 4 0.106 0.204 109.6 151 203.5 268 346.5 6.01e+8 2.95e+8 1.51e+8 8.08e+7 4.5e+7 0.332 0.426 0.531 0.651 0.784 440.5 551.5 680.5 829.5 1 000 2.6e+7 1.56e+7 9.56e+6 6.06e+6 5.13e+6 0.931 1.09 1.26 1.45 2.14

本文重點研究氣泡群BND沿流向的變化及相對分布的特征，取LES計算結果具有代表性的某一展向位置處（y=0 m）、某一垂向深度處（z=1 m）的一組流向速度參量作為氣泡群平衡方程的流場速度。

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

圖4 不同位置處氣泡相對數(shù)密度Fig.4 Relative BND at different location

圖（3）是BND總量沿x方向的分布，可以看出BND在離船尾500 m距離以內下降速度很快，之后衰減速度趨于緩慢，總體衰減規(guī)律與冪乘曲線近似。圖中黑方塊為文獻[4]中的實測數(shù)據(jù)，可以看出兩者形狀基本相同。圖（4）反映了不同位置處（x=500 m、x=1 500 m、x=3 000 m）的不同尺寸氣泡相對BND的分布情況。由圖可知尾流中的氣泡成分主要以尺寸小于200 μm為主。在離船尾越近的地方，小氣泡的相對BND越大，隨著距離的增加，峰值開始向200 μm靠攏，但這種變化不大。在x=3 000 m的位置處，BND最大的氣泡尺寸約為70 μm。這些結論都與相關文獻和實驗測量的結論[4-5]比較接近，進一步表明了本文使用BPBE方程對艦船氣泡尾流進行求解的可行性。

4 結 語

艦船尾流一直是國內外研究的熱點，尤其是氣泡尾流的研究，在海洋工程、艦船安全、雷達監(jiān)測以及武器制導等諸多領域，都具有重要的理論應用價值。針對當前對艦船尾流進行數(shù)值模擬的現(xiàn)有方法的復雜性及計算效率低下的缺陷，本文基于LES模擬和修正的BPBE方程，實現(xiàn)了尾流速度和氣泡運動的分開求解，大大提高了計算效率，且具有普遍適用性，在普通臺式機上即可實現(xiàn)模型的求解。計算結果表明：BND沿尾流流向的分布基本符合指數(shù)分布，在500 m以內的區(qū)域衰減速度很快，500 m以外的區(qū)域衰減速度緩慢；在離船尾越近的地方，小氣泡的相對BND越大；在3 000 m距離處BND最大的氣泡尺寸約為70 μm。所得結論與相關文獻比較接近，可望為艦船氣泡尾流特征研究提供新的有效手段。

[1]田恒斗,金良安,遲 衛(wèi).船舶遠程尾流場中氣泡上浮運動模型[J].船舶力學,2010,14(11):1195-1201.Tian Hengdou,Jin Liangan,Chi Wei.Rising process model of bubble in the far wake of ship[J].Journal of Ship Mechanics,2010,14(11):1195-1201.

[2]田恒斗,金良安,王 涌,等.考慮單氣泡運動特性的艦船尾流氣泡分布研究[J].兵工學報,2011,32(9):1126-1130.Tian Hengdou,Jin Liangan,Wang Yong,et al.Study on wake bubble distribution based on single bubble motion characteristics[J].Acta Armamentarii,2011,32(9):1126-1130.

[3]Stewart M B,Miner E W.Bubble dynamics in a turbulent ship wake[R].NRL Memorandmm Report,6055,1987.

[4]朱東華,張曉暉,顧建農,等.艦船尾流及其氣泡數(shù)密度分布的數(shù)值計算[J].兵工學報,2011,32(3):315-320.Zhu Donghua,Zhang Xiaohui,Gu Jiannong,et al.Numerical calculation of ship wake and its bubble number density distribution[J].Acta Armamentarii,2011,32(3):315-320.

[5]顧建農,田雪冰,張志宏,等.艦船主尺度參數(shù)對艦船氣泡尾流幾何特性的影響[J].艦船科學技術,2012,34(1):6-10.Gu Jingnong,Tian Bingxue,Zhang Zhihong,et al.The effect of the ship main scale parameter upon ship bubble wake geometry characters[J].Ship Science and Technology,2012,34(1):6-10.

[6]陳圣濤,王慧麗,王運鷹,等.艦船氣泡尾流特性的數(shù)值模擬和實驗研究[J].船舶力學,2012,16(4):342-349.Chen Shengtao,Wang Huili,Wang Yunying,et al.Numerical simulation and experimental study of a ship’s bubble wake[J].Journal of Ship Mechanics,2012,16(4):342-349.

[7]楊 力,劉慧開,等.艦船尾流氣泡特征的運動規(guī)律研究[R].武漢:海軍工程大學,2005.Yang Li,Liu Huikai,et al.Movement regularities study of Bubbles in Ship Wakes[R].Wuhan:Naval University of Engineering,2005.

[8]Smirnov A,Celik I,Shi S.LES of bubble dynamics in wake flows[J].Computers&Fluids,2005(34):351-373.

[9]Carrica P M,Drew D,Bonetto F,et al.A polydisperse model for bubbly two-phase flow around a surface ship[J].International Journal of Multiphase Flow,1999,25:257-305.

[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]徐麥容,劉成云.水中浮升氣泡的半徑和速度變化[J].大學物理,2008,27(11):14-17.Xu Mairong,Liu Chengyun.The change of radius and velocity of the rising bubble in water[J].College Physics,2008,27(11):14-17.

[12]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.

[13]Miner E W,Ramberg S E.Method for approximating the initial data plane for surface wake simulation[R].US:RL-MR-6376,1988.

[14]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.

[15]石晟瑋,蔣興舟,石 敏,等.艦船尾流氣泡運動特性研究[J].武漢理工大學學報（交通科學與工程版）,2007,31(5):764-767.Shi Shengwei,Jiang Xingzhou,Shi Min,et al.Movement characteristic of bubbles in ship wakes[J].Journal of Wuhan U-niversity of Technology(Transportation Science&Engineering),2007,31(5):764-767.