吳霆鋒, 杜永成, 楊立, 王保霖, 張澤華
(海軍工程大學(xué) 動力工程學(xué)院,湖北 武漢 430033)
潛艇航行過程中排出大量溫度較高的冷卻水,冷卻水浮升擴(kuò)散,可能使海面上尾流區(qū)域出現(xiàn)明顯的冷熱特征,將這種尾流稱之為熱尾流。理論上,基于紅外技術(shù)對潛艇尾流水面冷熱特征進(jìn)行分析,可以獲得潛艇的航向、航速、艇型等至關(guān)重要的作戰(zhàn)信息。同時(shí),海洋環(huán)境中存在鹽度環(huán)境,即海水鹽度可能隨深度不同而變化的現(xiàn)象,鹽度環(huán)境對于冷卻水的浮升擴(kuò)散可能存在一定影響,進(jìn)而影響海面尾流區(qū)域的冷熱特征。因此,開展鹽度環(huán)境下的潛艇熱尾流研究有重要意義。
楊立等[1]通過實(shí)驗(yàn)探究了穩(wěn)定溫度密度分層環(huán)境對于水下航行器熱尾流浮升的抑制作用;劉興泉[2]整理得到了東海P-N斷面夏季溫鹽度的分布特征;Bonnier等[3]對線性穩(wěn)定鹽層中球體尾跡的擴(kuò)散過程進(jìn)行了實(shí)驗(yàn)研究;Fritts等[4]基于DNS方法模擬了潛艇尾跡演化過程;Meunier等[5]基于流體的自保性,建立了線性分層流體中拖曳尾跡的一般模型;戴天奇等[6]基于Fluent動網(wǎng)格法對潛艇水下航行的二維計(jì)算模型進(jìn)行了仿真,結(jié)果表明,動網(wǎng)格方法具有較高精度,能反映尾流浮升過程中的瞬時(shí)變化;王平等[7-8]基于重疊網(wǎng)格技術(shù)進(jìn)行了均勻水體中的潛艇熱尾流仿真模擬,并與傳統(tǒng)來流法進(jìn)行了對比,發(fā)現(xiàn)了重疊網(wǎng)格技術(shù)具有更高精度。
綜上所述,國內(nèi)外對于潛艇尾流在海洋環(huán)境下的發(fā)展過程已經(jīng)有了一些理論和實(shí)驗(yàn)研究,但大多采用傳統(tǒng)來流法,未采用更為先進(jìn)與精確的重疊網(wǎng)格方法;且多數(shù)未考慮海洋鹽度對尾流的影響研究,缺少系統(tǒng)性研究海洋鹽度對于潛艇熱尾流的影響。因此,本文將基于重疊網(wǎng)格技術(shù),開展不同鹽度環(huán)境下的潛艇尾流仿真,采用有限體積法進(jìn)行模擬計(jì)算,探究海洋鹽度對于潛艇尾流浮升擴(kuò)散與水面冷熱特征的影響。
對海洋鹽度環(huán)境中的水下航行潛艇尾流浮升擴(kuò)散過程進(jìn)行數(shù)值模擬,需滿足如下控制方程。
連續(xù)性方程:
(1)
動量方程:
(2)
能量方程:
(3)
組分傳輸方程:
(4)
式中:T為溫度,K;ρ為流體密度,kg/m3;ui、uj(i、j=1,2,3)為x、z、y方向的速度分量,m/s;gi為重力加速度分量,m/s2;μ為粘性系數(shù),Pa·s;λ為導(dǎo)熱系數(shù),W/(m·K);cp為定壓比熱容,J/(kg·K);cs為組分的體積濃度;ρcs為組分的質(zhì)量濃度;Ds為組分的擴(kuò)散系數(shù),m2/s。
采用物性多項(xiàng)式構(gòu)建海水密度與溫度的關(guān)系[9]:
ρ(T)=A+BT+CT2+DT3+ET4+FT5
(5)
式中各項(xiàng)系數(shù)如表1所示。
表1 海水物性多項(xiàng)式參數(shù)[9]Table 1 Polynomial parameters of seawater physical properties
以縮比潛艇模型為研究對象,按照1∶100的比例建立計(jì)算模型。模型中計(jì)算域的長×寬×高為10.32 m×1.0 m×1.126 2 m,分別沿x、z、y方向,艇長為86 cm,艇體直徑為10 cm,兩側(cè)冷卻水排放口直徑為10 mm,以初始時(shí)潛艇中心為模型原點(diǎn)。潛艇航速為0.1 m/s,冷卻水排放速度為0.05 m/s,溫度為350 K。冷卻水排放口設(shè)置為速度入口,艇體周圍與背景水域交界面設(shè)置為重疊界面,頂部水面出口為壓力出口,四周壁面為對稱壁面,底部設(shè)置為靜止壁面。
圖1為基于重疊動網(wǎng)格的潛艇尾流三維計(jì)算模型示意圖。
圖1 基于重疊動網(wǎng)格的潛艇尾流三維計(jì)算模型示意Fig.1 Schematic diagram of the submarine 3-D calculation model based on overset grid
在Fluent Msehing平臺進(jìn)行模型網(wǎng)格劃分。潛艇及近艇區(qū)域的重疊部分采用Poly-Hexcore(多面體-六面體混合網(wǎng)格)進(jìn)行網(wǎng)格劃分,背景海水區(qū)域采用Cartesian(笛卡爾網(wǎng)格)方法進(jìn)行結(jié)構(gòu)化網(wǎng)格劃分。為進(jìn)行網(wǎng)格無關(guān)性驗(yàn)證,將模型網(wǎng)格數(shù)目分別劃分為101萬、202萬、301萬、407萬,令潛艇進(jìn)行t=9 s的航行仿真,選取距離潛艇后方0.2 m處截面,截面垂直于潛艇的運(yùn)動方向。測得截面上熱尾流中心位置處的溫度,結(jié)果如圖2所示,可以看出,當(dāng)網(wǎng)格數(shù)目達(dá)到202萬時(shí),截面上溫度不再隨網(wǎng)格數(shù)的增加而變化,故選定總網(wǎng)格數(shù)目為202萬。
圖2 網(wǎng)格無關(guān)性驗(yàn)證Fig.2 Mesh amount independence verification
湍流模型選用Realizablek-ε模型。壓力-速度耦合模式選用Coupled算法,壓力項(xiàng)選用Body Force Weighted,動量方程、能量方程、湍動能和耗散率均選用二階迎風(fēng)格式。采用瞬態(tài)計(jì)算,時(shí)間步長設(shè)置為0.01 s,計(jì)算時(shí)長為60 s。
海洋鹽度指的是海水中全部溶解固體與海水重量之比,通常以每千克海水中所含的克數(shù)表示。在海水溶解物中,NaCl占據(jù)主要部分。因此,可將鹽度近似表示為海水中鹽類物質(zhì)的質(zhì)量分?jǐn)?shù)。據(jù)資料顯示,世界大洋的平均鹽度為35‰。
一方面,海洋鹽度受緯度影響,因海域所處的緯度位置不同而具有不同的鹽度;另一方面,鹽度還主要受到蒸發(fā)量和降水量之差的影響。水分蒸發(fā)使得海水濃縮,從而鹽度上升,降水使海水稀釋,鹽度下降。在本文中,計(jì)算域尺寸為10.32 m×1.0 m×1.126 2 m,緯度變化較小,因此主要考慮蒸發(fā)量和降水量之差的影響。
當(dāng)海水蒸發(fā)量大于降水量時(shí),海水上層鹽度上升,下層海水鹽度維持穩(wěn)定,故呈現(xiàn)出上層鹽度高、下層鹽度低的正梯度鹽度分層分布,此類正梯度鹽度分層屬于不穩(wěn)定分層,海水會受重力作用產(chǎn)生上下對流,分層存在時(shí)間短暫,故不予考慮。當(dāng)海水蒸發(fā)量小于降水量時(shí),海水上層鹽度下降,下層海水鹽度變化緩慢,基本維持原有鹽度,因此呈現(xiàn)出上層鹽度低、下層鹽度高的負(fù)梯度鹽度分層分布,目前海洋環(huán)境中以該類分布為主,這也是本文主要探討的鹽度分層環(huán)境。
研究鹽度分層環(huán)境時(shí),簡化為海水鹽度沿深度方向線性變化,當(dāng)鹽度分層的上下邊界密度均給定時(shí),即可計(jì)算得出表征分層流場特征的浮頻率數(shù)N,其計(jì)算式為:
(6)
式中:ρ0為上邊界流體密度,kg/m3;Δρ為上下邊界流體密度差,kg/m3;H為流場深度,m。浮頻率數(shù)作為分層流場的特征數(shù),當(dāng)其較大時(shí),代表流場密度分層跨度較大,密度分層梯度大;反之,當(dāng)浮頻率數(shù)較小時(shí),代表流場的密度分層梯度較小。
因NaCl溶質(zhì)為固體,其物性與液態(tài)水存在較大差距,模擬其混合溶解過程存在較大難度,且誤差較大。因此,考慮將NaCl溶液作為替代性介質(zhì),引入體積數(shù)20%NaCl溶液,其主要物性參數(shù)如下[10]。
20%NaCl溶液的密度多項(xiàng)式為:
ρNaCl(20%)=-0.001 3T2+0.28T+1 176.37
20%NaCl溶液的粘度多項(xiàng)式為:
μN(yùn)aCl(20%)=-4.80×10-9T3+4.88×10-6T2-
1.66×10-3T+0.19
20%NaCl溶液的導(dǎo)熱系數(shù)多項(xiàng)式為:
kNaCl(20%)=-5.75×10-6T2+
4.81×10-3T-0.33
20%NaCl溶液的比熱容多項(xiàng)式為:
cpNaCl(20%)=0.94T+3148.68
在Fluent中引入體積分?jǐn)?shù)20%NaCl溶液后,通過初始化Ptach功能設(shè)置其分布,模擬形成海洋鹽度環(huán)境,具體設(shè)置情況見表2。
表2 鹽度分層參數(shù)設(shè)置Table 2 Parameter of temperature stratification
本節(jié)將對均勻鹽度與鹽度分層2類海水環(huán)境數(shù)值仿真的結(jié)果進(jìn)行分析。一方面,海面溫度特征是潛艇熱尾流探測最直接的依據(jù),明顯的海面溫度特征容易被紅外探測捕捉到,進(jìn)而暴露潛艇蹤跡,因此,潛艇航行時(shí)的海面溫度場是結(jié)果分析的重點(diǎn)。另一方面,對潛艇冷卻水浮升擴(kuò)散過程的研究,能更深刻地認(rèn)識潛艇熱尾流的形成機(jī)理。
因此,可分別通過對海水溫度場與冷卻水浮升擴(kuò)散過程進(jìn)行分析,研究鹽度對于潛艇尾流的影響。
2.1.1 海面溫度特征
選取上界面溫度場進(jìn)行分析,圖3為t=60 s時(shí)不同鹽度的均勻海水環(huán)境下的海洋表面溫度云圖。觀察圖可以看到,不同鹽度下均勻海水環(huán)境的溫度分布基本一致,合理推斷均勻海水環(huán)境下,環(huán)境鹽度的改變對潛艇尾流形成的海洋表面溫度場影響非常微弱。
圖3 均勻鹽度下的海洋表面溫度云圖Fig.3 Temperature cloud of sea surface in uniform sea water
2.1.2 冷卻水浮升擴(kuò)散
在研究潛艇冷卻水的浮升擴(kuò)散過程中,冷卻水中心的浮升高度與溫度變化是重要指標(biāo)。冷卻水中心浮升高度即冷卻水中心的y坐標(biāo),可以直觀反映冷卻水的浮升軌跡;冷卻水中心的溫度變化則能清楚反映冷卻水與海水的熱交換程度。
圖4為潛艇冷卻水中心溫度隨航行距離變化圖,圖5為冷卻水中心浮升高度(即冷卻水中心的y坐標(biāo))隨航行距離變化的規(guī)律。觀察圖4,可以發(fā)現(xiàn),3種不同均勻鹽度環(huán)境下冷卻水溫度變化與浮升過程基本一致,均勻海水下整體鹽度的改變對于潛艇冷卻水的浮升擴(kuò)散影響微弱,與上述溫度場分布結(jié)果吻合,與上述結(jié)論相互印證。需要說明的是,在圖5中出現(xiàn)冷卻水中心短暫下降的現(xiàn)象,該時(shí)段為冷卻水浮升至潛艇艉翼附近,向上浮升受到艉翼阻擋后反彈,故而短暫下沉。
圖4 均勻鹽度環(huán)境下冷卻水中心溫度Fig.4 Central temperature of cooling water in uniform seawater
圖5 均勻鹽度環(huán)境下冷卻水中心浮升高度Fig.5 Buoyancy height of cooling water center in uniform seawater
2.2.1 海面溫度特征
圖6所示為t=60 s時(shí)的負(fù)鹽度梯度分層環(huán)境下的海面溫度云圖。觀察圖6可以看到,在鹽度分層環(huán)境下,當(dāng)鹽度梯度增大時(shí),浮頻率隨之增大,海面尾跡區(qū)域明顯縮小,溫度降低。因此可以得出結(jié)論,鹽度分層環(huán)境下海面溫度特征隨鹽度梯度的增大而減弱。
圖6 鹽度分層環(huán)境下的海面溫度云圖Fig.6 Temperature cloud of sea surface in salinity stratification environment
為了更清晰地研究海面溫度變化,選取海面上直線AB與直線CD,進(jìn)行兩直線上的溫度分布分析,直線AB、CD分布如圖6所示。圖7、圖8分別為直線AB與CD上的溫度分布。觀察發(fā)現(xiàn),當(dāng)鹽度分層梯度增大,即浮頻率增大時(shí),海面溫度分布趨勢保持不變,尾流區(qū)域整體溫度降低。因此可以得出結(jié)論:鹽度分層梯度的增大對潛艇熱尾流形成的海面溫度場具有一定的抑制作用。
圖7 直線AB溫度分布Fig.7 Temperature distribution on line AB
圖8 直線CD溫度分布Fig.8 Temperature distribution on line CD
圖9所示為海面最高溫度隨浮頻率數(shù)的變化規(guī)律,可以發(fā)現(xiàn),浮頻率數(shù)N=0的3種均勻鹽度環(huán)境的海面最大溫度基本相同,最大溫度均處于291.024 K左右,而當(dāng)鹽度存在分層,且分層梯度產(chǎn)生變化時(shí),隨著浮頻率數(shù)N逐漸增大,海面最大溫度先略微上升,達(dá)到291.026 K后降低,降低速度隨N的增大而減小,整體變化趨勢接近正弦函數(shù)。
圖9 海面最高溫度隨浮頻率數(shù)N變化規(guī)律Fig.9 Variation of maximum sea surface temperature with buoyancy frequency N
2.2.2 冷卻水浮升擴(kuò)散
圖10所示為t=60 s時(shí)不同鹽度分布下的中性面溫度云圖。分析圖10可以發(fā)現(xiàn),在鹽度分層環(huán)境下,潛艇冷卻水浮升受到明顯抑制。當(dāng)浮頻率數(shù)增大時(shí),冷卻水被釋放后向后方擴(kuò)散的趨勢加劇,向上浮升的趨勢受阻。
圖10 鹽度分層下的中性面溫度云圖Fig.10 Temperature cloud of neutral surface in salinity stratification environment
圖11所示為潛艇冷卻水中心浮升高度,可以發(fā)現(xiàn),隨著浮頻率數(shù)的增大,潛艇冷卻水浮升受到明顯抑制,最大浮升高度有所降低,甚至出現(xiàn)到達(dá)浮升最高點(diǎn)后下沉的現(xiàn)象。圖12為鹽度分層環(huán)境下潛艇冷卻水中心溫度隨航行距離的變化規(guī)律??梢钥闯?當(dāng)浮頻率較小時(shí),冷卻水中心溫度較低,而冷卻水排放之初溫度均為350 K,說明浮頻率數(shù)較低時(shí),冷卻水與周圍海水的熱交換更加劇烈,而當(dāng)浮頻率數(shù)增大時(shí),熱交換減弱,從而導(dǎo)致尾流浮升受到抑制,進(jìn)而使得海洋表面的尾流區(qū)域溫度整體降低。
圖11 鹽度分層下冷卻水中心浮升高度Fig.11 Central buoyancy height of cooling water under salinity stratification
圖12 鹽度分層下冷卻水中心溫度Fig.12 Central temperature of cooling water under salinity stratification
2.2.3 旋渦結(jié)構(gòu)
潛艇熱尾流問題本質(zhì)上接近于橫流環(huán)境下的熱射流排放問題,研究潛艇尾流所引起的渦結(jié)構(gòu),有助于更深入認(rèn)識潛艇熱尾流的發(fā)展演變規(guī)律。圖13所示為N=0.076 87時(shí)的鹽度分層下的潛艇后方y(tǒng)z平面速度矢量圖,位置為x=-5.4 m。觀察發(fā)現(xiàn),艇體上下均存在一對反向旋渦。潛艇兩側(cè)排水口向外排放冷卻水時(shí),冷卻水受排放慣性作用向下運(yùn)動,便產(chǎn)生了艇體下方的旋渦對;繼而高溫冷卻水受到浮升力作用向上浮升,在艇體上方形成了另一對旋渦。
圖13 N=0.076 87時(shí)潛艇后方y(tǒng)z平面速度矢量圖Fig.13 Velocity vector diagram of yz plane behind submarine when N=0.076 87
圖14所示為上下旋渦對的旋渦半徑隨浮頻率數(shù)的變化規(guī)律,可以發(fā)現(xiàn),旋渦半徑隨浮頻率增大而減小,減小速度隨浮頻率增大而逐漸降低,且艇體下方旋渦半徑的減小速度比上方旋渦更快。因此,鹽度分層梯度的增大對于潛艇尾流所引起的旋渦具有抑制作用,促使其影響面積收縮,抑制作用隨著鹽度分層梯度的增大逐漸減弱。同時(shí),鹽度分層對于艇體下方的旋渦抑制作用更為明顯。
圖14 上下旋渦對半徑隨浮頻率數(shù)N的變化規(guī)律Fig.14 Variation of vortex radius of upper and lower vortex pairs with buoyancy frequency N
1)均勻海水環(huán)境下,環(huán)境鹽度的改變對冷卻水的浮升擴(kuò)散影響非常微弱。
2)鹽度分層環(huán)境對潛艇冷卻水浮升具有明顯抑制作用,海面溫度特征隨鹽度梯度的增大而減弱,尾流區(qū)域最高溫度與浮頻率基本呈正弦函數(shù)關(guān)系。
3)鹽度分層對于潛艇航行產(chǎn)生的反向旋渦具有抑制作用,抑制作用隨著鹽度分層梯度的增大逐漸減弱。