劉慧開(kāi)，張勸華，趙小龍，王星宇，楊 立

(1.海軍后勤部，北京 100036; 2.海軍裝備部，北京 100036;3.海軍工程大學(xué)動(dòng)力工程學(xué)院，武漢 430033)

水下航行器在通氣管航態(tài)下，需要將發(fā)動(dòng)機(jī)產(chǎn)生的高溫廢氣排入海水中，廢氣進(jìn)入海水將破裂成大小不一的氣泡，氣泡與海水混合形成具有復(fù)雜結(jié)構(gòu)的氣液兩相流。為了掌握水下航行器通氣管航行時(shí)水下排氣的熱特征，有必要對(duì)氣泡在海水中的運(yùn)動(dòng)特性以及氣泡與海水的氣液兩相流進(jìn)行深入研究。當(dāng)連續(xù)地將發(fā)動(dòng)機(jī)廢氣排放到液體中去，由于氣體自身的壓力與周?chē)h(huán)境壓力的不平衡以及氣液交界面上的表面張力等作用，氣體將破裂成大小不一的氣泡，在浮力力作用下，氣泡將向上浮升。上升的氣泡誘導(dǎo)周?chē)囊后w向上運(yùn)動(dòng)，形成氣泡羽流流動(dòng)［1-3］。對(duì)于氣泡羽流的研究，國(guó)內(nèi)外都開(kāi)展了一些研究工作，但還很不成熟。1996年哈爾濱工程大學(xué)張文平［4］等利用攝像機(jī)拍攝和重放技術(shù)，觀察研究了柴油機(jī)水下排氣從管口到水面氣液兩相流動(dòng)狀態(tài)以及管口氣泡形成、上浮和破碎的過(guò)程。發(fā)現(xiàn)整個(gè)兩相流動(dòng)狀態(tài)可明顯的分為氣泡形成、自由氣泡、氣泡云以及水面噴涌4 個(gè)區(qū)域。氣泡羽流的運(yùn)動(dòng)場(chǎng)通?？梢詣澐譃? 個(gè)區(qū)域:形成區(qū)、形成后區(qū)和表面流區(qū)。在形成區(qū)，氣流破碎成氣泡，并與環(huán)境液體混和，羽流寬度和軸線(xiàn)流速快速增長(zhǎng); 而形成后區(qū)羽流寬度和軸線(xiàn)流速的增加要緩慢得多; 在液體表面附近，羽流轉(zhuǎn)向水平方向流動(dòng)，形成具有垂向速度梯度的表面流區(qū)［5］。早期，國(guó)內(nèi)外對(duì)氣泡羽流的研究主要集中在實(shí)驗(yàn)研究和積分理論分析階段［6，7］。近年來(lái)，對(duì)氣泡羽流的數(shù)值研究逐漸增多，采用的數(shù)值方法主要有有限差分法、有限分析法和混合有限分析法。1994年林永明［8］等采用有限差分法對(duì)氣泡羽流紊流場(chǎng)進(jìn)行了數(shù)值分析和模擬，與Kobus 的實(shí)驗(yàn)結(jié)果符合較好。1999年Robert［9］基于完整的兩相流模型對(duì)氣泡羽流進(jìn)行了二維和三維的數(shù)值模擬。2000年武漢水利電力大學(xué)的王雙峰、槐文信［10］等從兩相流理論中已被廣泛接受的雙流體概念出發(fā)，建立了描述兩相湍流流動(dòng)兩方程模型，用于均勻環(huán)境中圓形氣泡羽流的數(shù)值計(jì)算，所得結(jié)果初步證實(shí)了其正確有效性。2001年馬霞［11］等假設(shè)氣泡在水中的運(yùn)動(dòng)僅僅是一種紊動(dòng)擴(kuò)散現(xiàn)象，提出了一個(gè)預(yù)測(cè)圓形氣泡特性的k-ε-E湍流模型，對(duì)靜止環(huán)境中的圓形氣泡羽流場(chǎng)進(jìn)行數(shù)值模擬，模擬結(jié)果與實(shí)測(cè)資料吻合較好。2006年重慶大學(xué)劉洪濤［12］建立了氣液兩相流動(dòng)傳熱的三維穩(wěn)態(tài)數(shù)學(xué)物理模型，給出了描述不同模型下的換熱及流動(dòng)控制方程。2012年中國(guó)船舶研究中心的喻太君等［13］對(duì)水下沖壓發(fā)動(dòng)機(jī)噴管內(nèi)的氣液兩相流進(jìn)行了研究。

綜上可知，前人對(duì)于靜止環(huán)境中氣泡羽流的速度場(chǎng)研究較多，而對(duì)于水下高溫排氣的氣液兩相流速度場(chǎng)和溫度場(chǎng)的數(shù)值模擬研究幾乎為空白。因此，有必要開(kāi)展水下航行器水下高溫排氣的氣液兩相流數(shù)值模擬研究。本文通過(guò)建立水下排氣的連續(xù)性方程、動(dòng)量守恒方程和能量守恒方程，采用氣液兩相流混合模型和k -ε 紊流模型，利用有限體積數(shù)值計(jì)算方法，對(duì)水下排氣的速度場(chǎng)和溫度場(chǎng)進(jìn)行了研究，得到了水下排氣在海洋表面形成的溫度分布，研究成果對(duì)水下航行器排氣的紅外抑制有重要的理論和應(yīng)用價(jià)值。

1 水下排氣的氣液兩相流模型

目前多相流研究通常有四種模型:Lagrangian 離散模型、歐拉模型(Eulerian Model)、VOF 模型和混合模型(Mixture Model)。模擬水下排氣與海水熱交換的速度場(chǎng)與溫度場(chǎng)時(shí)，可以采用多相流中的混合模型。

連續(xù)性方程

能量方程

其中keff是有效熱傳導(dǎo)率，keff=k +kt，kt代表湍流熱傳導(dǎo)率。SE代表所有體積熱源。

第二相體積分?jǐn)?shù)方程

2 網(wǎng)格劃分與數(shù)值計(jì)算方法

由于氣液兩相流模擬計(jì)算的復(fù)雜性，有必要對(duì)水下航行器排氣口模型進(jìn)行一定的簡(jiǎn)化。由于航行器外形主要對(duì)其后的流場(chǎng)產(chǎn)生擾動(dòng)，而排氣口大多布置在指揮臺(tái)上，所以外形對(duì)排氣與海水熱交換所形成的溫度場(chǎng)影響不是很大，本文忽略航行器外形對(duì)水下排氣的影響。假設(shè)排氣口半徑為15 cm，航行深度為排氣口距離水面4 m，計(jì)算區(qū)域海水溫度為20℃。主流場(chǎng)采用笛卡爾坐標(biāo)系:以排氣口中心為坐標(biāo)原點(diǎn)(0，0，0)，排氣管長(zhǎng)0.5 m。其中y 軸正方向?yàn)槲擦髋艢怆x開(kāi)艇體的方向，y 軸負(fù)方向?yàn)楹叫蟹较颍瑉 軸正方向?yàn)楹Ｋ疃确较?，主流?chǎng)尺寸3 m×10 m×5 m，其數(shù)值計(jì)算實(shí)體模型如圖1 所示。排氣口采用四面體網(wǎng)格，在y =0 面采用三角形網(wǎng)格結(jié)構(gòu)，流場(chǎng)區(qū)域采用混合網(wǎng)格結(jié)構(gòu)。通過(guò)網(wǎng)格無(wú)關(guān)性檢驗(yàn)，計(jì)算域共生成網(wǎng)格1 675 608 個(gè)網(wǎng)格單元，其網(wǎng)格劃分情況如圖2 所示。

利用Fluent 軟件進(jìn)行計(jì)算，通過(guò)3D 穩(wěn)態(tài)分離求解器求解，重力方向?yàn)閦 軸負(fù)方向，采用k-ε 紊流模型，定義海水為第一相，排氣為第二相。主相、氣相進(jìn)口為速度邊界條件，出口為Outflow，其他面為壁面條件。壓力速度耦合格式為壓力隱式分裂算子(PISO)，壓力采用PRESTO 分離，動(dòng)量與能量方程采用QUICK 離散格式，相體積分?jǐn)?shù)采用一階離散方式。通過(guò)殘差值檢測(cè)器和兩相進(jìn)出口質(zhì)量流量檢測(cè)器來(lái)判斷計(jì)算是否收斂，當(dāng)殘差值達(dá)到指定差值、兩相進(jìn)出口質(zhì)量流量趨于相等時(shí)，認(rèn)為計(jì)算已經(jīng)收斂。

圖1 數(shù)值計(jì)算實(shí)體模型

圖2 網(wǎng)格劃分

3 結(jié)果分析與討論

3.1 排氣速度場(chǎng)

圖3 為中心軸線(xiàn)垂直縱截面上氣相的速度分布圖，由圖可以看出水下排出氣體在排氣口大約7.5 m 處浮出水面。排氣速度在排氣口附近最大值為45.2 m/s，在離開(kāi)排氣口之后氣體速度迅速減小，至排氣口3 m 遠(yuǎn)時(shí)排氣速度的最大值已經(jīng)降至12.4 m/s，至浮出水面時(shí)氣相速度基本上已經(jīng)降至海水速度4.12 m/s。

圖3 中心軸線(xiàn)垂直縱截面上氣相的速度分布

圖4為不同垂直截面上氣相的速度分布，從圖4 可以看出，離排氣口由遠(yuǎn)及近，不同垂直截面上的氣相速度逐漸減小，氣相不斷向四周環(huán)境水體中擴(kuò)散，氣相的分布區(qū)域逐漸增大，離排氣口不遠(yuǎn)處的氣體擴(kuò)散表現(xiàn)為圓形擴(kuò)散，后來(lái)逐漸擴(kuò)散成馬蹄形;每個(gè)截面都存在一個(gè)最大速度區(qū)域，且最大速度區(qū)域都位于擴(kuò)散區(qū)域的中心區(qū)域。在y=0.6m 處，中心最大速度可達(dá)44 m/s，氣體擴(kuò)散至y =9m 處，速度降至4 m/s 左右，基本上已經(jīng)達(dá)到海水的速度。

圖4 不同垂直截面上氣相的速度分布

圖5為不同深度水平截面上氣相的速度分布，從圖5 可知，不同深度水平截面的氣相速度分布都類(lèi)似于“流星狀”，“流星狀的尾巴”在離開(kāi)排氣口的方向不斷變寬，在“流星”的根部中心區(qū)存在一個(gè)最大速度，且越靠近水面，這個(gè)速度越小，當(dāng)氣相浮升至垂直高度4 m(浮出水面)時(shí)，氣相速度基本已經(jīng)降至海水的速度。在z =0 處，可以看到氣體在排氣口附近的最大速度為45 m/s，之后不斷的衰減，形成明顯的速度梯度。

圖5 不同深度水平截面上氣相的速度分布

綜合圖3 ～圖5 可以發(fā)現(xiàn)排氣在離開(kāi)排氣口之后速度迅速減小，排氣在進(jìn)入環(huán)境水體之后不斷向水中擴(kuò)散上浮，在浮出水面時(shí)，氣泡速度已經(jīng)降至環(huán)境水體速度。

3.2 排氣溫度場(chǎng)

圖6 為中心軸線(xiàn)垂直縱截面上氣相的溫度分布，由圖6可以看出水下高溫排氣在排氣口大約7.5 m 處浮出水面，溫度降至與海水的溫差只有0.05 K。在排氣出口煙氣的溫度為400 K，在離開(kāi)排氣口之后迅速減小，至排氣口2.5 m 遠(yuǎn)時(shí)排氣溫度已經(jīng)降至294 K。由于溫度差異太大，為了能更好的表現(xiàn)出氣體在浮出水面的溫度場(chǎng)，圖6 中的溫度云圖最大值為294 K，從中也可以看出氣體在400 K 降至294 K 時(shí)所經(jīng)歷的時(shí)間是非常的短的。

圖6 中心軸線(xiàn)垂直縱截面上氣相的溫度分布

從圖7 可以看出，離排氣口由遠(yuǎn)及近，不同垂直截面上的氣相溫度逐漸減小，氣相不斷向四周擴(kuò)散，分布區(qū)域逐漸增大，離排氣口不遠(yuǎn)處的擴(kuò)散表現(xiàn)為圓形擴(kuò)散，后來(lái)逐漸擴(kuò)散成馬蹄形，且每個(gè)截面都存在一個(gè)最大溫度區(qū)域。對(duì)比圖4，我們發(fā)現(xiàn)雖然溫度場(chǎng)與速度場(chǎng)的變化有類(lèi)似之處，但溫度場(chǎng)的擴(kuò)散區(qū)域要明顯小于速度場(chǎng)擴(kuò)散的區(qū)域，這說(shuō)明溫度的降低速率要比速度降低的速率要快得多。水下航行器排氣在擴(kuò)散至y=3 m 處時(shí)，其最大溫度已經(jīng)降至僅有293.5 K。

圖7 不同垂直截面上氣相的溫度分布

圖8為不同深度水平截面上氣相的溫度分布，從圖8 可知，不同深度水平截面的氣相溫度分布都類(lèi)似于“流星狀”，“流星狀的尾巴”在離開(kāi)排氣口的方向不斷變寬，在“流星”的根部中心區(qū)存在一個(gè)最大溫度，且越靠近水面，這個(gè)溫度越小，當(dāng)氣相浮升至離排氣口4 m(浮出水面)時(shí)，氣相溫度與海水的溫度差降至0.05 K。對(duì)比圖5，我們同樣可以發(fā)現(xiàn)氣泡的溫度梯度要比速度梯度大很多，氣泡從400 K 降至394 K 所需要的時(shí)間是非常短的。

綜合圖6 ～圖8 可以發(fā)現(xiàn):水下航行器排氣在離開(kāi)排氣口之后溫度迅速降低，氣相溫度在離開(kāi)排氣口水平距離3 m時(shí)，其最大溫度已經(jīng)降至僅有293.5 K;排氣在進(jìn)入環(huán)境水體后不斷向水中擴(kuò)散上浮，在浮出水面時(shí)，氣泡溫度已經(jīng)降至與環(huán)境溫度水體相差0.05 K;氣相擴(kuò)散的范圍在剛離開(kāi)排氣口時(shí)基本為圓形，但是在浮出水面之后基本充滿(mǎn)了整個(gè)計(jì)算區(qū)域。

3.3 不同排氣溫度下的水下航行器排氣特征

在不改變其他初始條件下，改變水下航行器排氣在排氣口的初始溫度，可以發(fā)現(xiàn)降低水下航行器排氣在排氣口的初始溫度，可有效降低水下航行器排氣浮升至水面時(shí)與周?chē)Ｋ臏囟炔?。初始溫度?23 K 的排氣在離開(kāi)排氣口浮升至水面時(shí)，氣相與周?chē)Ｋ淖畲鬁夭顬?.11 K。當(dāng)初始溫度從523 K 降至400 K 時(shí)，浮升至海面時(shí)與海水的溫度差從0.11 K 降至0.05 K，進(jìn)一步降低初始溫度至333 K 時(shí)，浮升至海面時(shí)與海水的溫度差將降至0.013 K?？梢?jiàn)，降低水下航行器水下排氣口的初始溫度，可有效降低排氣浮升至水面時(shí)與周?chē)Ｋ臏囟炔睢?/p>

圖8 不同深度水平截面上氣相的溫度分布

4 結(jié)論

運(yùn)用有限體積數(shù)值計(jì)算方法對(duì)水下航行器水下排氣形成的氣泡羽流的浮升過(guò)程進(jìn)行了數(shù)值模擬，得到了水下排氣形成的羽流速度場(chǎng)和溫度場(chǎng)。

在溫度均勻的海水中水下排氣羽流浮升時(shí)，氣相的速度不斷減小，最后與海水的速度大小基本相等; 在氣相速度減小的過(guò)程中，每一個(gè)縱截面都存在一個(gè)速度最大的中心區(qū)域，中心區(qū)域的形狀由圓形慢慢發(fā)展成馬蹄形。氣相速度場(chǎng)的水平截面區(qū)域類(lèi)似于“流星狀”，且可看出明顯的速度梯度存在。

水下排氣羽流浮升時(shí)，氣相的溫度在海水中不斷降低，由出口溫度400 K 降至與海水的溫差為0.05 K;在溫度降低的過(guò)程中，每一個(gè)縱截面都存在一個(gè)溫度最大的中心區(qū)域，中心區(qū)域的形狀由圓形慢慢發(fā)展成為馬蹄形;氣相溫度場(chǎng)的水平截面區(qū)域也類(lèi)似于“流星狀”，且可看出明顯的溫度梯度存在。與氣相的速度場(chǎng)，溫度場(chǎng)發(fā)展的區(qū)域要明顯比速度場(chǎng)發(fā)展的區(qū)域“要窄”。

降低水下航行器排氣的初始溫度，可有效降低排氣浮升至水面時(shí)與周?chē)Ｋ臏囟炔睢．?dāng)初始溫度從523 K 降至400 K 時(shí)，浮升至海面時(shí)與海水的溫度差從0. 11 K 降至0.05 K，進(jìn)一步降低初始溫度至333 K 時(shí)，浮升至海面時(shí)與海水的溫度差將降至0.013 K。

［1］吳鳳林，Tsang G.關(guān)于氣泡羽流的研究［J］.水動(dòng)力學(xué)研究與進(jìn)展，1989，4(1):104-112.

［2］吳鳳林，Tsang G.三維氣泡羽流建立區(qū)流動(dòng)形態(tài)的觀察［J］.力學(xué)學(xué)報(bào)，1991，23(1):1-8.

［3］Friedl M J，F(xiàn)annelop T K.Bubble plumes and theirs interaction with the water surface［J］.Applied Ocean Research，2000(22):119-128.

［4］張文平，張?zhí)煸?，劉志剛，?柴油機(jī)水下排氣口氣液兩相流動(dòng)狀態(tài)研究［J］.內(nèi)燃機(jī)工程，1996，17(3):14-18.

［5］王雙峰，槐文信，李煒.靜止均勻環(huán)境中氣泡射流特性的研究——數(shù)學(xué)模型及實(shí)驗(yàn)資料的分析［J］.水科學(xué)進(jìn)展，2000，11(1):25-31.

［6］Kobus H E.Proc.of 11th Conf.on coastal Eng［C］.London，ASCE.New York，p 106，1968.

［7］Fannelop T K，Sjoen K.AIAA 18th Aerospace Sci.Meeting［C］.Pasadena，CA.，1980.

［8］林永明，李煒.氣泡羽流紊流的應(yīng)力代數(shù)模型及數(shù)值分析［J］.水道港口，1994，21(2):32-36.

［9］Robert F.Mudde，Oliver simonin.Two-and three-simulations of a bubble plume using a two-fluid model［J］.Engineering Science，1999，54:5061-5069.

［10］王雙峰，槐文信，李煒.靜止均勻環(huán)境中氣泡射流特性的研究——形成區(qū)及形成后區(qū)特性的數(shù)值研究［J］.水科學(xué)進(jìn)展，2000，11(1):32-37.

［11］馬霞，李建中，魏文禮，等.氣泡羽流的數(shù)值模擬［J］.西安理工大學(xué)學(xué)報(bào)，2001，17(1):86-89.

［12］劉洪濤. 有內(nèi)熱源熔池內(nèi)的氣液兩相流流動(dòng)與傳熱［D］.重慶:重慶大學(xué)，2006.

［13］喻太君，王寶壽.水下沖壓發(fā)動(dòng)機(jī)噴管內(nèi)兩相流動(dòng)［J］.四川兵工學(xué)報(bào)，2012，33(3):4-7.