毛德龍，王康碩，沙 江，何 凜，任海剛，王怡豪，王 奔

(中國船舶工業(yè)系統(tǒng)工程研究院，北京 100094)

0 引 言

水下航行體因其極強的隱蔽性與突襲性而在軍用領(lǐng)域占有重要地位，其中以潛艇最具代表性[1]。近年來，受益于光電技術(shù)的飛速發(fā)展，紅外探測設(shè)備的分辨率、精度和抗干擾能力大大提高，紅外探潛作為一種重要的非聲探潛措施，以其探測范圍大、晝夜工作、被動探測等優(yōu)點而被廣泛關(guān)注[2-4]。水下航行體在航行過程中，其動力系統(tǒng)產(chǎn)生的廢熱會被冷卻水吸收，并隨冷卻水排放到海水中。這些冷卻水的溫度明顯高于周圍水體的溫度，在密度差的驅(qū)使下向上浮升，有可能在海面形成溫度異常的紅外特征。

目前，相關(guān)學(xué)者針對密度分層對熱射流流動的影響進行了大量研究，并取得了一系列成果。江傳富等[5]提出了熱尾流不同發(fā)展階段的物理模型，分析了其在密度均勻和分層流體中的浮升規(guī)律；張昊春等[6]建立了溫度分層環(huán)境和均勻介質(zhì)環(huán)境中的射流模型，計算得到了熱尾流的溫度分布特性及浮升擴散過程中的衰減規(guī)律；楊立等[7]通過對線性密度分層流體中熱尾流的浮升過程研究，發(fā)現(xiàn)了流體的分層會抑制熱尾流的浮升。

海水的溫度分布受諸多因素的影響，如深度、地理位置、季節(jié)、太陽輻射和海洋大氣熱交換等。海水溫度隨深度的增加而顯著遞減，根據(jù)NOAA（National Oceanic and Atomspheric Administration）數(shù)據(jù)庫WOD(World Ocean Database)可得海水溫度沿著深度方向的變化規(guī)律[8]，不同緯度地區(qū)海水總體溫度變化值約在30℃左右；同一地區(qū)，在50 m深度范圍內(nèi)，海水溫度變化范圍約在1℃～2℃，靠近海平面處（深度20 m以內(nèi)）溫度沿深度變化不明顯，原因為海上風(fēng)浪大，海平面處海水運動劇烈，不同深度海水之間對流混合，溫度趨于一致；深度大于20 m后水溫迅速下降。海水鹽度是指海水中全部溶解固體與海水質(zhì)量之比，也會隨深度的變化而發(fā)生改變。根據(jù)海水溫度、鹽度在深度方向上的變化情況，總結(jié)了40° N附近50 m深度內(nèi)，密度變化規(guī)律，如圖1所示。在深度方向上，海水的密度為線性增加，50 m深度下密度的增長幅度為0.35‰。需要指出的是，大量文獻中對表層海水密度分層的處理方式也是線性的，與圖1中的規(guī)律有一致性[9,10]。

本文在海水密度沿深度方向變化規(guī)律的基礎(chǔ)上，使用數(shù)值仿真手段模擬熱射流在密度分層海水中的運動，探究密度分層對熱射流流動的具體影響關(guān)系，矯正未分層情況下數(shù)值仿真的預(yù)測結(jié)果，增加數(shù)值仿真預(yù)測結(jié)果的可靠性。

1 數(shù)學(xué)模型

數(shù)值仿真采用流體連續(xù)介質(zhì)模型假設(shè)，其流動受納維斯托克斯（N-S）方程控制[11]。N-S方程的張量形式如下：

式中：fi為流體受到的體積力，此處為熱射流受到周圍冷流體浮力與重力的合力。

以溫度形式表達的能量方程為：

其中 α為熱擴散率。根據(jù)Boussinesq假設(shè)，忽略密度變化在非穩(wěn)態(tài)項、對流項和擴散項中的影響，即在上述3項中密度取常數(shù)，只考慮密度變化在浮力項中的影響。

密度與溫度的關(guān)系式為：

式中：TRef為參考溫度，取環(huán)境的水溫，其他所有物質(zhì)特性均取參考狀態(tài)下的介質(zhì)屬性。

密度不分層情況下，Boussinesq假設(shè)參考密度ρRef是一個定值，其在空間中的分布是常數(shù)，而密度分層情況下 ρRef沿著深度方向是變化的。本文研究目的在于探究密度變化對熱射流的影響，因此實際海水中的鹽度與溫度分層將會對密度分層產(chǎn)生影響。為了控制變量，海水溫度變化也僅僅體現(xiàn)在密度項中，參考溫度TRef依然取常數(shù)。需要補充說明的是，式（4）描述了溫度與密度的關(guān)系，在此基礎(chǔ)上可以得出單位體積海水所受的重力為：

所受浮力為：

式（6）與式（5）相減，得到浮力與重力的合力大小為：

由于數(shù)值仿真軟件設(shè)置的實際需要，參考密度ρRef仍取常數(shù)，而其影響完全放在膨脹系數(shù) β 中考慮。為了使浮力與重力的合力保持不變，定常膨脹系數(shù)β變?yōu)?β′，兩者關(guān)系如下：

其中： ρRef為常數(shù)； ρ (h)為密度分層在深度h方向上的變化情況。

2 CFD模擬

2.1 計算區(qū)域與邊界條件

2.1.1 計算區(qū)域

數(shù)值仿真的計算區(qū)域如圖2所示，整個計算區(qū)域x方向（航行方向）長為120 m，z方向?qū)挒?0 m，y方向（垂直方向）高為70 m。噴嘴（即熱射流出口）離頂端自由液面高度為50 m，距離入口平面長度為5 m。整個計算區(qū)域劃分為1，2，3三個區(qū)域。根據(jù)計算精度的需要，分別對這3個計算區(qū)域劃分了疏密不同的網(wǎng)格，其中，噴嘴位于計算區(qū)域1中，其網(wǎng)格質(zhì)量將會影響整個熱射流流場的分布，因此對該區(qū)域的網(wǎng)格進行了加密處理，尤其在噴嘴周圍使用了高密度的O型網(wǎng)格。

圖2 計算區(qū)域Fig.2 Computational region

2.1.2 邊界條件

計算區(qū)域左側(cè)面設(shè)置為入口邊界條件，入口海水溫度為15℃；噴嘴也設(shè)置為入口邊界條件，熱射流的速度大小為0.33m/s。右側(cè)面設(shè)置為出口邊界條件，其余外側(cè)面均為自由滑移壁面。計算區(qū)域內(nèi)部的接觸面設(shè)置為interfaces。

2.2 工況計算

圖1定量描繪了40° N海域，50 m水深海水密度的典型分布，不難發(fā)現(xiàn)其變化范圍極小。因此，定義海水的密度分層值φ為：50 m深度海水密度 ρ50與表面海水密度 ρ0的差值 (ρ50-ρ0)，與表面海水密度 ρ0的比值（表面海水密度 ρ0為定常數(shù)），即

圖3給出了相同時刻，密度不分層φ=0%與密度分層值φ=1%的溫度云圖對比，可以看出兩者幾乎沒有差別。這說明，φ 值較小的密度分層對熱射流流場的影響并不明顯。

圖3 t=16min，xy截面溫度云圖Fig.3 xy section temperature nephogram in 16 min

為了在數(shù)值計算結(jié)果中體現(xiàn)密度分層對熱射流的影響，需使密度分層的影響遠大于數(shù)值計算誤差，因此數(shù)值仿真計算選取水下航行體航速V為0 kn，0.1 kn，0.2 kn，密度分層值φ為0%，50%，100%，熱射流與海水的溫差 ΔT為30℃，40℃，50℃。顯然數(shù)值仿真中φ值50%，100%的計算結(jié)果并不能直接應(yīng)用于實際環(huán)境，為了解決φ值偏小無法體現(xiàn)密度分層影響以及φ值偏大與實際不符的問題，將對仿真結(jié)果進行曲線擬合，得出任意φ 值情況下，密度分層對熱射流流場的影響規(guī)律。

3 密度分層對熱射流影響的仿真分析

3.1 密度分層影響的定性說明

采用熱射流高度H作為衡量熱射流影響的指標，假設(shè)海水溫度為15℃，若受熱射流影響使溫度上升0.01℃，圖4給出了熱射流某時刻溫度15.01℃的等值面圖，任意時刻熱射流高度H定義為此時噴嘴中心到T=15.01℃溫度等值面頂端的垂直距離。

圖4 15.01℃等值面Fig.4 15.01 °C isosurface

密度不分層情況下，隨著航速增加，相同時刻熱射流的高度H逐漸降低；隨著熱射流溫度的提高，相同時刻熱射流的高度H逐漸升高。圖5與圖6對比了相同密度分層值=50%情況下，水下航行體航速、熱射流溫度對熱射流流場的影響，結(jié)果表明上述規(guī)律依然成立。由圖5可知，在熱射流出口溫度為45℃情況下，航速0.2 kn相對于航速0.1 kn，其熱射流15.01℃等值面頂端較低，即熱射流高度H較低；由圖6可知，在航速0.1 kn情況下，熱射流出口溫度55℃相較于出口溫度45℃，其熱射流高度H較高。

圖5 t=4 min，熱射流15.01℃等值面Fig.5 Hot jet 15.01 °C isosurface at 4 min

圖6 t=4 min，熱射流15.01℃等值面Fig.6 Hot jet 15.01 °C isosurface at 4 min

由φ的定義可知，φ值越大，則某一固定深度處的海水密度ρRef越大。式說明此時熱射流所受浮力與重力的合力Fi也越大。因此φ值越大，熱射流漂浮的速度越快，相同時刻熱射流的高度H也越高。

3.2 密度分層影響的修正

為了體現(xiàn)密度分層影響對熱射流流場的影響，引入密度分層影響的修正系數(shù) ω ，定義式如下：

其中：H0為密度不分層（即φ =0），某一確定的時間t、航速V、溫差 ΔT工況下的熱射流高度；Hφ為密度分層值φ時，相同時間t、航速V、溫差 ΔT工況的熱射流高度。若是能夠得出修正系數(shù) ω ，則任意密度分層值φ時的熱射流高度為：

圖7（a）列出了 ΔT=30?C ，V=0.1 kn，不同φ值的熱射流高度曲線；圖7（b）列出了 ΔT=30?C，V=0.2 kn，不同φ 值的熱射流高度曲線。通過圖7發(fā)現(xiàn)，在溫差ΔT、航速V相同時，隨著密度分層值φ的增加，熱射流的漂浮速度是逐漸增加的，即相同時間t熱射流高度Hφ增大。密度分層值φ表征了相同深度海水的密度大小，φ 值越大，表征海水密度越大；而熱射流所受浮力大小Fi與海水密度正相關(guān)，即相同時刻熱射流的高度Hφ更大。

圖7 ΔT=30?C ，不同φ 值的熱射流高度曲線Fig.7 At ΔT=30?C, Heat jet height curves with differentφ values

圖8列出了 ΔT=30?C，φ =50%，不同航速V的熱射流高度曲線。由圖8可知，在密度分層φ =50%情況下，隨著水下航行體航速V由0增加到0.2 kn，熱射流相對于周圍海水的速度變大，兩者的對流換熱更加劇烈，導(dǎo)致熱射流受到的浮力變小，熱射流高度Hφ逐漸降低

圖8 ΔT=30?C ，φ =50%，不同航速V的熱射流高度曲線Fig.8 At ΔT=30?C andφ =50%, Hot jet height curves with different speedV

圖9列出了V=0.1 kn，φ =50%，不同溫差 ΔT熱射流高度曲線。由圖9可知，在密度分層φ =50%情況下，隨著溫差 ΔT由30℃上升到50℃，熱射流在與周圍海水混合過程中溫度較高，其上浮過程中受到的浮力增大，熱射流高度Hφ逐漸升高。

圖9 V=0.1 kn，φ =50%，不同溫差 ΔT熱射流高度曲線Fig.9 AtV=0.1 kn andφ =50%, Hot jet height curve with different temperature difference ΔT

圖10～圖14分別列出了改變密度分層值φ 、航速V或溫差 ΔT其中一個變量時，密度分層修正系數(shù)ω與時間t的關(guān)系曲線，不難發(fā)現(xiàn)圖中各密度分層修正系數(shù)ω均是隨著時間t的增加而下降的。圖10～圖12顯示在溫差 ΔT與航速V固定時，隨著密度分層值φ的增大，密度分層修正系數(shù)ω曲線向上移動，這是由于分層值φ 的增加，使得熱射流所受浮力增大，相同時間內(nèi)熱射流高度Hφ越大，相應(yīng)的密度分層修正系數(shù)ω值越大。從圖13可以看出，隨著航速V的增大，密度分層修正系數(shù)ω曲線向上移動。由于在密度不分層情況下，航速V的增大使得熱射流高度H0減小，根據(jù)定義式(10)有 ω =Hφ/H0，因此在航速V較大情況下，H0較小使得密度分層的影響更加突出，密度分層修正系數(shù)ω曲線上移。從圖14可以看出，隨著溫差 ΔT的提高，在密度不分層情況下，熱射流高度H0是逐漸增大的，使得密度分層的影響減小，密度分層修正系數(shù)ω曲線下移。

圖10 ΔT=30?C ，V=0.1 kn，不同φ 值的密度分層修正系數(shù)ω曲線Fig.10 At ΔT=30?C andV=0.1 kn, Density stratificationcorrection coefficientω curves with differentφ values

圖11 ΔT=30?C ，V=0.2 kn，不同φ 值的密度分層修正系數(shù)ω曲線Fig.11 At ΔT=30?C andV=0.2 kn, Density stratificationcorrection coefficientω curves with differentφ values

圖12 ΔT=40?C ，V=0.1 kn，不同φ 值的密度分層修正系數(shù)ω曲線Fig.12 At ΔT=40?C andV=0.1 kn, Density stratification correction coefficientω curves with differentφ values

圖13 ΔT=30?C ，φ =50%，不同航速V下密度分層修正系數(shù)ω曲線Fig.13 At ΔT=30?C andφ =50%, Density stratification correction coefficientω curves with different speedsV

圖14 V=0.1 kn，φ =50%，不同溫差 ΔT下密度分層修正系數(shù)ω曲線Fig.14 AtV=0.1kn andφ ==50%, Density stratification correction coefficientω curves with different temperature differenceΔT

由圖10～圖14可知，密度分層值φ 、航速V、溫差ΔT以及時間t均會影響密度修正系數(shù) ω。假設(shè) ω關(guān)系式如下：

當φ =0時，熱射流高度不需要修正，即 ω =1，因此關(guān)聯(lián)式中第一項為常數(shù)1；其次，第二項設(shè)置成指數(shù)函數(shù)的原因在于，隨著密度分層值φ 、航速V（kn）、溫差 ΔT（℃）以及時間t（min）的變化，密度修正系數(shù) ω 均是呈現(xiàn)單調(diào)變化規(guī)律的，這與指數(shù)型函數(shù)相符；顯然當航速V=0 kn，密度分層值φ ≠0時，熱射流高度依然需要修正，即 ω ≠1。若航速影響項為Vb，則有 ω =1，因此指數(shù)函數(shù)項出現(xiàn)了 (1 +V)b。 ΔT=0表征熱射流與周圍環(huán)境不存在溫差，不會上?。籺=0表征熱射流還沒有排入海水中，因此溫差與時間兩項直接寫成指數(shù)的形式即可。通過Matlab對數(shù)值仿真中8個工況，共66個數(shù)據(jù)點進行擬合，得出密度修正系數(shù) ω 擬合式如下：

該式的驗證范圍是：0 kn≤V≤ 0.2 kn，0≤ φ ≤100%，30℃≤ ΔT≤50℃，1 min≤ t≤ 12 min。

該擬合關(guān)聯(lián)式的計算結(jié)果與數(shù)值仿真結(jié)果對比，誤差在6%以內(nèi)，且絕大部分數(shù)據(jù)點的誤差都在1%以下。因此在確定密度分層值φ 、航速V、溫差 ΔT以及時間t的情況下，可以方便準確地通過式(13)求得密度修正系數(shù) ω 。再利用式(11)即可求得對應(yīng)的密度分層值φ 、航速V、溫差 ΔT以及時間t下，熱射流的高度Hφ。

4 結(jié) 語

水下航行體航行過程中，其動力系統(tǒng)的廢熱以冷卻水的形式向外界直接排出，在密度差的作用下向海面浮升。實際海洋環(huán)境下，海水密度隨水深是變化的，而熱射流浮升力的大小受海水密度的影響?；诖?，本文研究了密度分層對熱射流流場的影響，結(jié)論如下：

1）隨著海水深度的增加，海水的鹽度上升，溫度下降，這兩者均會造成海水密度隨深度的增加而增大；

2）所得到的擬合關(guān)聯(lián)式計算結(jié)果與數(shù)值仿真結(jié)果對比，誤差在6%以內(nèi)，且絕大部分數(shù)據(jù)點的誤差都在1%以下。