基于三維動網(wǎng)格技術的潛艇熱尾流浮升規(guī)律及水面溫度特征研究

來慶志，王成安，譚建宇，周志權

(1.哈爾濱工業(yè)大學（威海）汽車工程學院，山東 威海 264200；2.哈爾濱工業(yè)大學（威海）信息與電氣工程學院，山東 威海 264200)

0 引 言

潛艇，特別是續(xù)航能力強的核潛艇，在海洋中航行時，由于需要冷卻動力裝置，因此會向海洋中排放大量的冷卻水而形成熱尾流[1]。在浮升力作用下熱尾流漂浮至海面，引起海面溫度分布產(chǎn)生變化。通過分析海面溫度特征，理論上可以獲得潛艇航向、位置等重要信息，從而實現(xiàn)對潛艇的跟蹤與打擊。因此，開展?jié)撏嵛擦鞲∩?guī)律的研究具有重要意義。

早在1972年Merrit等[2]便開始了對潛艇熱尾流溫度分布特性的研究，他們利用熱敏電阻對模擬的潛艇無動量尾流進行探測。1988年，俄羅斯科研人員，研究指出對于潛艇而言，熱尾流紅外探測是一種有效的非聲探測方法[3]。之后國內(nèi)外學者又進行了大量的實驗研究，并取得了一系列的重要研究成果[4–7]。但熱尾流浮升擴散是一個極為復雜的過程，實驗研究成本高、周期長，且易受測量方法及儀器精度的限制。隨著計算機技術的發(fā)展，數(shù)值方法以其低成本、高效率的優(yōu)勢，而越來越成為一種重要的科研手段。為方便研究，早期研究人員大都忽略海面空氣流動的影響，將海面簡化為壁面或壓力邊界。但真實環(huán)境中海面是大氣和海洋的交界面，海面上空氣流動、交界面內(nèi)氣液兩相能量、動量交換都會對研究結(jié)果產(chǎn)生重要影響。吳猛猛等[8]和張昊春等[1]采用VOF法在二維空間內(nèi)對潛艇熱尾流浮升擴散過程進行仿真研究，指出在水下熱尾流浮升主要受浮力和海水垂直溫度梯度的影響。到達水面后熱尾流將與氣液交界面發(fā)生摻混，逐漸擴散和衰減。研究均采用相對運動法模擬潛艇航行，即將潛艇模型視為靜止，根據(jù)相對運動原理，讓海水以與航速相同的速度向模型流動。這與潛艇航行的實際情況仍有一定差距。戴天奇等[9]基于動網(wǎng)格技術，在二維平面內(nèi)對多種工況下潛艇冷卻水排放過程進行仿真。研究表明動網(wǎng)格技術對于潛艇航行過程中熱尾流浮升過程有較高的計算精度。Shi Zhi-guang等[10]采用動網(wǎng)格技術，對不考慮海面空氣流動條件下潛艇熱尾流海面溫度分布特征及其紅外識別過程進行仿真分析，并獲得了熱尾流水面溫度分布及其紅外成像特征。

綜上所述，國內(nèi)外學者雖然對熱尾流浮升規(guī)律已經(jīng)進行了廣泛的研究，但在三維空間內(nèi)考慮海面空氣流動影響的熱尾流水面溫度分布特征方面的研究尚不多見。本文基于動網(wǎng)格技術，分別采用壁面法與VOF法對三維海域內(nèi)潛艇熱尾流浮升過程進行仿真研究。獲得了潛艇熱尾流的浮升、擴散規(guī)律及其水面溫度分布與紅外識別特征，分析了水面空氣流動及大氣-海水交界面對熱尾流發(fā)展過程的影響，并通過對不同冷卻水排放方案下熱尾流水面溫度分布特征的變化規(guī)律的分析，提出了抑制熱尾流浮升的方法。本文研究結(jié)果可以為潛艇紅外探測及紅外隱身方面研究提供理論基礎。

1 數(shù)理模型

1.1 計算模型

圖1為本文研究模型示意圖，模型長1 200 m，水面寬200 m，水深60 m（在VOF法中增加10 m空氣層，模型總高度為70 m）。潛艇在水下40 m深度處沿海水流動方向航行。

本文以 “海狼”級核潛艇為研究對象，表1為其基本參數(shù)，研究重點在于潛艇熱尾流在水下浮升過程及其在水面溫度分布特征。為方便計算，文中忽略了潛艇瞭望塔及方向舵等結(jié)構，將其簡化為長110 m，最大直徑10 m，排水口直徑2 m的長柱體結(jié)構，如圖1所示。

圖1 計算模型示意圖Fig.1 Schematic diagram of simulation model

表1 “海狼”級潛艇參數(shù)表Tab.1 Parameters of Seawolf-class submarine

1.2 網(wǎng)格技術

圖2為計算模型網(wǎng)格局部示意圖，為提高計算精度，對潛艇周圍及氣液交界面處網(wǎng)格進行了局部加密處理。

網(wǎng)格劃分質(zhì)量及網(wǎng)格數(shù)量對求解結(jié)果的精度有重要的影響。為保證本文模擬結(jié)果的精準性，文中網(wǎng)格質(zhì)量均大于0.4，角度質(zhì)量均大于36°，最大縱橫比均小于30。為驗證網(wǎng)格無關性，本文分別劃分了278萬、516萬、900萬、1 022萬和1 796萬5種數(shù)量的網(wǎng)格，分別在距離潛艇尾部20 m，40 m，60 m，80 m，100 m處測量熱尾流中心溫度，如圖3所示?？梢钥闯?，當網(wǎng)格數(shù)量大于516萬時，尾流中心溫度不再隨網(wǎng)格數(shù)量增加而改變。因此，為節(jié)約計算時間，本文選擇數(shù)量為516萬的網(wǎng)格進行計算。

1.3 流動傳熱控制方程

質(zhì)量方程

圖2 局部網(wǎng)格示意圖Fig.2 Schematic diagram of local mesh

圖3 網(wǎng)格無關性驗證Fig.3 Validation of mesh independence

動量方程

能量方程

其中：

第p相的體積分數(shù)滿足方程：

式中：ρ為密度，kg/m3；ui（i=1，2，3）分別為x，y，z方向速度分量，m/s；p為流體靜壓，Pa；gi為重力加速度分量，m/s2；cp為定壓比熱容，J/（kg·K）；T為溫度，K；k為導熱系數(shù)，W/m·K；μ為粘性系數(shù)，Pa·s；αp為第p相體積分數(shù)。

VOF法是目前求解分層流動、有自由表面流動等存在2種或多種互不相溶流體的復雜多相流動問題中最常用的方法。它在計算域內(nèi)對互不相溶流體求解同一個動量方程，并追蹤每種流體的體積分數(shù)來模擬多相流動。計算單元內(nèi)各相體積分數(shù)之和為1，在單元內(nèi)，若αp=0，則單元中沒有第p相流體；αp=1，單元中充滿第p相流體；0＜αp＜1，單元中有第p相流體與其他相流體的界面。根據(jù)局部αp值，計算域內(nèi)每一控制容積被賦予適當?shù)奈镄院妥兞恐?。湍流流動是非常復雜的流動，在計算湍流流動時需要附加湍流方程，本文研究均采用標準k-ε湍流模型。

1.4 物性參數(shù)及邊界設置

參照海水物性參數(shù)隨溫度變化數(shù)據(jù)，擬合得計算所用相應物性參數(shù)關系式如下：

潛艇排放冷卻水與螺旋槳的攪動作用效果可視為“熱射流”[1]。針對目標潛艇，其熱功率為220 MW，若熱量全部通過冷卻水散出，假定海水溫度為284 K，冷卻水溫度為318 K，海水熱容取4 000 J/（kg·K），則排水速度為0.52 m/s；假定潛艇航速為1 m/s，海水入口速度為0.1 m/s；為減弱壁面影響，將海水出口及兩側(cè)面定義為自由流；壁面法中將海面與海底均定義為壁面條件；VOF法中海面上空增加10 m空氣層，海面為空氣與海水交界面，空氣層上表面給定4 m/s水平速度。

2 計算結(jié)果

2.1 壁面法熱尾流浮升擴散規(guī)律

圖4為縱截面上熱尾流溫度分布云圖。以潛艇尾部排水口為原點，沿潛艇航行反方向，測量距潛艇尾部不同距離上熱尾流中心區(qū)溫度及其所在高度，可得尾流中心溫度與浮升高度隨其在水平方向擴散距離變化曲線如圖5和圖6所示。為方便分析，將圖6中浮升高度曲線與水平方向（圖中X軸方向）的夾角定義為熱尾流的“浮升角”。顯然，浮升角越大，熱尾流浮升速率越高，當浮升角為90°時熱尾流垂直向上浮升，此時浮升速率最高；當浮升角為0°時浮升速率最低為0。

觀察圖5和圖6可得，在距離小于15 m區(qū)域內(nèi)，由于熱尾流剛從潛艇尾部排出，與周圍環(huán)境存在較大的溫差，會迅速與環(huán)境發(fā)生劇烈的能量、動量交換，尾流溫度迅速降低；這階段浮升角處于最大，浮升高度迅速增加；隨著距離增加，尾流與海水間溫差逐漸降低，熱量交換逐漸減弱、浮升角逐漸減小，熱尾流溫度衰減、浮升速率逐漸降低；到達水面后，浮升角度減小為0°，熱尾流不再浮升，在水面逐漸擴散、衰減，最終形成水滴形溫度輪廓，如圖7所示。據(jù)此，可以將熱尾流水下發(fā)展過程大致分為：近尾流快速擴散、浮升區(qū)、過渡區(qū)、遠尾流水面擴散區(qū)3部分。

圖4 縱截面溫度分布云圖Fig.4 Temperature distribution in the longitudinal section

圖5 熱尾流中心溫度隨其水平擴散距離變化曲線Fig.5 Variation of central temperature with horizontal diffusion distance of thermal wake

圖6 熱尾流浮升高度隨其水平擴散距離變化曲線Fig.6 Variation of rising height varied with horizontal diffusion distance of thermal wake

圖8分別為 3～5 μm 及 8～12 μm 波段下水面紅外識別圖像，可以看出，當溫差大于0.06 K時，熱尾流具有較為明顯紅外識別輪廓；在8～12 μm波段下，熱尾流紅外成像更加清晰，識別效果更佳。

圖7 海面溫度分布云圖Fig.7 Temperature distribution on sea surface

圖8 水面紅外圖像Fig.8 Infrared image of sea surface

2.2 VOF法熱尾流浮升擴散規(guī)律

海面是大氣與海洋相接形成的一個動態(tài)變化的氣液交界面，它對大氣與海洋之間熱量、動量交換等具有重要的影響。為研究獲得更接近于真實海況條件下的熱尾流浮升規(guī)律及水面溫度分布特性，本文采用VOF法對海面空氣流動、氣液交界面對熱尾流的影響進行了研究。

圖9為采用VOF方法獲得的縱截面熱尾流溫度分布云圖。圖10為2種方法下尾流浮升高度隨水平擴散距離變化曲線，圖中給出了VOF法中熱尾流浮升過程中3個區(qū)域劃分?？梢钥闯?，在近尾流區(qū)2種方法熱尾流浮升速率基本相等，但是由于VOF方法中不再將海面簡化為壁面邊界，而是隨海水一起流動的氣液交界面，受海面空氣流動及交界面作用，海水流動加強，換熱量增加，導致過渡區(qū)熱尾流浮升角減小，浮升速率降低，從而使得浮升過渡段增長，浮升過程能量衰減增加；當浮升至水面后，受空氣流動及交界面的影響，水面溫度場出現(xiàn)明顯抖動失穩(wěn)現(xiàn)象，導致水面溫度輪廓不再連續(xù)，如圖11所示。

圖9 縱截面溫度分布云圖Fig.9 Temperature distribution in the longitudinal section

圖10 熱尾流浮升高度隨擴散距離變化曲線Fig.10 Variation of rising height with horizontal diffusion distance of thermal wake

圖11 海面溫度分布云圖Fig.11 Temperature distribution on sea surface

為方便分析，將圖11中由熱尾流引起的水面高溫區(qū)域稱為“溫度異常區(qū)”。圖12為溫度異常區(qū)與環(huán)境海水溫度之差大于某一值所占區(qū)域面積（AT）的折算尺寸（將溫度異常區(qū)轉(zhuǎn)換為與其面積相等的正方形的邊長長度（Lc），即Lc2=AT）變化曲線?？梢钥闯觯捎谑芎Ｃ婵諝饬鲃蛹昂?氣交界面的作用，熱尾流浮升速率降低，浮升過程溫度衰減增加，從而導致VOF法中浮升至水面的溫差特征減弱，在壁面法中溫差大于0.06 K的區(qū)域折算尺寸可達201 m，而在VOF法中減小了約55%，僅為91 m?？梢姡Ｃ婵諝饬鲃蛹昂?氣交界面對熱尾流的浮升過程及表面溫度分布特征具有重要的影響，VOF法對于有自由表面的多相流動具有較高的計算精度，仿真研究中可采用VOF法研究海面空氣流動及交界面對熱尾流的影響。

圖12 溫差異常區(qū)域的折算尺寸Fig.12 Conversion size of temperature abnormal area

圖13 為 VOF 法在 3～5 μm 及 8～12 μm 波段下水面紅外識別圖，可以看出，熱尾流依然具有較好的紅外識別性。

圖13 水面紅外圖像Fig.13 Infrared image of sea surface

2.3 不同排水方案熱尾流浮升擴散規(guī)律

潛艇冷卻水出口溫度、出口速度對熱尾流浮升、擴散過程及表面溫度分布特征都具有重要影響。潛艇正常航行時，冷卻水出口速度與溫度應滿足關系式：

式中：T1和T0分別為冷卻水出口溫度與環(huán)境溫度。

可以看出，相同熱工況下，冷卻水排出速度與溫度成反比，增加冷卻水出口溫度，則出口速度會減小。當冷卻水出口溫度增加時，熱尾流與環(huán)境溫差增加，熱尾流溫度信號增強；當冷卻水排出速度減小時，則單位時間內(nèi)排出的熱量減小，熱尾流信號減弱。因而，理論上存在一組冷卻水排出溫度與速度的最佳值，可以使?jié)撏嵛擦饕鹚鏈囟刃盘栕钊?、表面熱紅外特征最不明顯。如表2所示，本文選取以下5組冷卻水排放方案，研究冷卻水排出速度與溫度對熱尾流浮升及水面溫度分布特征的影響規(guī)律。

表2 不同排水方案出口參數(shù)Tab.2 Parameters of different discharge programs

圖14為不同排水方案下熱尾流浮升高度隨距離的變化曲線，可以看出，隨著排水溫度升高，熱尾流浮升速率逐漸增加。這是由于，冷卻水出口溫度增加后，其與周圍環(huán)境的溫差增大，密度差增大，產(chǎn)生的浮升力增大，從而加快了熱尾流浮升速率。

圖14 不同排水方案熱尾流浮升高度隨擴散距離變化曲線Fig.14 Variation of the rising height with diffusion distance of thermal wake for different discharge programs

圖15 溫度異常區(qū)域的折算尺寸Fig.15 Conversion size of temperature abnormal area

圖15為不同排水方案下溫度異常區(qū)的折算尺寸變化曲線?？梢钥闯?，隨著排水溫度升高，熱尾流引起的水面溫度異常區(qū)面積逐漸增加，排水溫度由313 K增加至333 K，溫差大于0.06 K區(qū)域折算尺寸由60.5 m增大至109.5，增加約81%。這是由于排水溫度升高后，熱尾流浮升速率增加，浮升過程能量衰減減小，從而導致水面的熱尾流溫度信號增強。

通過以上分析可得，排水溫度變化對熱尾流浮升過程及表面溫度分布特性的影響要大于速度變化的影響，潛艇航行時可通過適當增加冷卻水排出速度，降低排出溫度，以減弱潛艇尾跡信號強度，增加自身隱身性能。

3 結(jié) 語

本文采用動網(wǎng)格技術對潛艇熱尾流浮升過程及其水面溫度分布與紅外識別特征進行了仿真分析，研究表明：

1）潛艇熱尾流浮升、擴散過程可分為：近尾流快速擴散、浮升區(qū)、過渡區(qū)、遠尾流表面擴散區(qū)三部分；當溫差大于0.06 K，熱尾流具有較為明顯紅外識別輪廓。

2）海面空氣流動及海-氣交界面對熱尾流浮升過程及其表面溫度分布具有重要影響：受空氣流動及交界面影響，熱尾流浮升速率降低，水下漂浮距離增加約68%，表面溫差大于0.06 K區(qū)域尺寸減小約55%。仿真研究中可采用VOF法研究海面空氣流動及交界面的對熱尾流的影響。

3）排水溫度對熱尾流浮升過程及其表面溫度分布具有重要影響，排水溫度增加20 K，表面溫差大于0.06 K區(qū)域尺寸增大約81%。潛艇航行時可通過適當增加冷卻水排出流量，降低排出溫度，以減弱潛艇尾跡信號強度，增加自身隱身性能。

[1]張昊春, 吉宇, 馬銳, 等.水下航行體熱尾流浮升特性研究[J].艦船科學技術, 2015, 07: 24–28.ZHANG Hao-chun, JI Yu, MA Rui, et al.Buoyant characteristics of thermal wakes discharged by underwater vehicles[J].Ship Science and Technology, 2015, 07: 24–28.

[2]MERRITT G E.Wake growth and collapse in stratified flow[J].AIAA Journal, 1973, 12(7): 940–949.

[3]張健, 楊立, 袁江濤, 等.水下航行器熱尾流試驗研究[J].實驗流體力學, 2008, 22(3): 7–14.ZHANG Jian, YANG Li, YUAN Jiang-tao, et al.Experimental research of underwater vehicle thermal wakes[J].Journal of Experiments in Fluid Mechanics, 2008, 22(3): 7–14.

[4]WREN G G, MAY D.Detection of submerged vessels using remote sensing techniques[J].Australian Defense Force Journal, 1997, 127(11): 9–15.

[5]MARION B, OLIVIER E.Experimental investigation of the collapse of a turbulent wake in a stratified fluid physics of fluids[J].Physics of fluids, 2002, 14(2): 791–801.

[6]THOMAS S L, DAVID C F, JOSEPH A W.Application of direct and large-eddy simulation methods to late wakes of submerged bodies[C]//IEEE Proceedings of the Users Group Conference, 2004: 112–118.

[7]張健, 陳翾, 楊立, 等.水下航行器排放冷卻水溫度分布特性研究[J].船舶力學, 2009, 04: 533–539.ZHANG Jian, CHEN Xuan, YANG Li, et al.Study of temperature characteristic of cooling water discharged by underwater vehicle[J].Journal of Ship Mechanics, 2009, 04:533–539.

[8]吳猛猛, 陳伯義, 張修峰, 等.溫度分層海洋中水下航行體引起水面特征參數(shù)變化規(guī)律的研究[J].紅外技術, 2010, 04:242–246.WU Meng-meng, CHEN Bo-yi, ZHANG Xiu-feng, et al.The study on the varied laws of surface feature parameters caused by a going body underwater in the temperature stratification ocean[J].Ship Science and Technology, 2015, 07: 24–28.

[9]戴天奇, 姚世衛(wèi), 魏志國.基于動網(wǎng)格技術的潛艇熱尾流浮升規(guī)律研究[J].艦船科學技術, 2015, 05: 86–89.DAI Tian-qi, YAO Shi-wei, WEI Zhi-guo.Numerical simulation of thermal wake buoyant law based on dynamic meshing technique[J].Ship Science and Technology, 2015, 05:86–89.

[10]SHI Zhi-guang, LI Ji-cheng.Detectability of the infrared surface features of the wake behind a moving underwater body[C]//AOPC2015: Optical and Optoelectronic Sensing and Imaging Technology, 2015.