張 木，譚俊杰，易文俊，蔡曉偉

（南京理工大學(xué)1.能源與動(dòng)力工程學(xué)院；2.瞬態(tài)物理重點(diǎn)實(shí)驗(yàn)室，南京210094）

空化是流動(dòng)中的常見(jiàn)現(xiàn)象，往往能導(dǎo)致空蝕、噪聲和機(jī)械效率下降等負(fù)面現(xiàn)象的產(chǎn)生，因此，在早期工程設(shè)計(jì)中，人們對(duì)空化現(xiàn)象的分析與研究主要是為了避免空化現(xiàn)象的產(chǎn)生[1].但隨著俄羅斯超空泡魚雷“暴風(fēng)”的問(wèn)世，人們發(fā)現(xiàn)利用超空泡技術(shù)可使水下航行體被包裹在空泡內(nèi)，僅頭部與水體接觸，從而減少航行體在水中的表面摩擦阻力與壓差阻力，提高水下航行體的運(yùn)動(dòng)速度和航行距離[2，3].因此，國(guó)內(nèi)外許多學(xué)者在空泡流動(dòng)的數(shù)值模擬方面開(kāi)展了大量的研究工作，如KINNAS和FINE[4]發(fā)展了一種基于速度勢(shì)的邊界元方法，用于二維水翼的局部和超空化流動(dòng)分析；文獻(xiàn)[5]采用回射流模型對(duì)三維水翼的局部空泡進(jìn)行了非線性分析；文獻(xiàn)[6]基于均相N－S方程，采用高雷諾數(shù)k－ξ湍流模型，對(duì)空化數(shù)為0.15條件下繞三維尖頭回轉(zhuǎn)體的空泡流動(dòng)進(jìn)行了數(shù)值模擬，其結(jié)果與實(shí)驗(yàn)數(shù)據(jù)吻合較好；文獻(xiàn)[7]基于simple算法，采用帶低雷諾數(shù)修正k－ξ湍流模型，模擬了在空化數(shù)為0.01條件下的繞三維軸對(duì)稱航行體的空泡流動(dòng)，得到航行體運(yùn)動(dòng)阻力系數(shù)與空化數(shù)之間的變化關(guān)系；文獻(xiàn)[8]利用CFD商業(yè)軟件對(duì)空化器線形與超空泡形態(tài)的關(guān)系及其阻力特性進(jìn)行了數(shù)值仿真，分析了影響超空泡航行體減阻效果的因素.以上這些數(shù)值仿真研究主要針對(duì)空化數(shù)大于0.01的條件下，利用雷諾平均湍流模型模擬水下航行體的空泡流動(dòng).

針對(duì)此種情況，本文基于N－S方程及PISO算法，引入考慮氣液間相變的空化模型，采用大渦模擬湍流模型和VOF方法，在空化數(shù)為0.78×10－3～4.88×10－3的條件下對(duì)繞三維水下自然超空泡射彈的空泡流動(dòng)進(jìn)行了數(shù)值模擬.

1 控制方程

1.1 大渦模擬

為求解非定常超空泡流動(dòng)，近似認(rèn)為氣液兩相流動(dòng)是等溫均相流動(dòng).另外根據(jù)文獻(xiàn)[9]的研究，在等溫均相流動(dòng)的前提下，空化流動(dòng)中氣相與液相的可壓縮性影響很小，而氣液混合區(qū)域，其可壓縮性影響需要考慮.因此，對(duì)于氣液混合區(qū)域，建立速度、壓力和密度的耦合關(guān)系式，考慮可壓縮性的影響.對(duì)于氣液兩相采用VOF方法捕捉氣液兩相間的界面，以及采用大渦模擬湍流模型對(duì)流場(chǎng)的湍流運(yùn)動(dòng)進(jìn)行模擬.關(guān)于氣液兩相之間相變問(wèn)題，本文采用KUNZ[6]空化模型進(jìn)行模擬.基于以上分析，控制方程可寫為

式中，φl(shuí)為液相的體積分?jǐn)?shù)，ρl為液相密度，δij為Kronecker函數(shù)，μeff為有效動(dòng)力粘性系數(shù)為相變質(zhì)量傳遞源項(xiàng)，ρm為氣液兩相的混合密度為濾波速度，為濾波壓力，ksgs為亞網(wǎng)格湍動(dòng)能，g為重力.變形張量1，2，3.為了求解μeff＝μm＋μsgs，需要求解μm和ksgs.由于ρm＝φl(shuí)ρl＋（1－φl(shuí)）ρv，ρv為氣相密度，因此μm＝φl(shuí)μl＋（1－φl(shuí)）μv.關(guān)于ksgs求解，采用k－Δ模型進(jìn)行封閉，其方程可以寫為

1.2 空化模型

式中，為濾波壓力，pv為空化壓力，t∞為特征時(shí)間，u∞為來(lái)流速度，Cp和Cd為經(jīng)驗(yàn)系數(shù)，本文中Cp和Cd分別為1 000和800，t∞＝L∞/u∞，L∞為特征長(zhǎng)度.

2 離散方法

本文采用有限體積法對(duì)控制方程進(jìn)行離散，其中時(shí)間項(xiàng)采用二階精度的三點(diǎn)向后隱式差分格式，對(duì)流項(xiàng)采用二階精度的線性迎風(fēng)格式，擴(kuò)散項(xiàng)采用二階中心差分格式.另外基于PISO算法，建立壓力、密度與速度的耦合算法，考慮氣液混合相可壓縮性的影響.在預(yù)測(cè)步驟中，將上一時(shí)刻的壓力分布作為當(dāng)前時(shí)刻動(dòng)量方程中的壓力分布p＊，然后求解如下方程：

式（8）中方程右端項(xiàng)為混合密度在p0點(diǎn)處的一階泰勒展開(kāi)式，而混合密度可通過(guò)式（1）、式（2）以及混合密度表達(dá)式導(dǎo)出.通過(guò)式（8）得出壓力修正量p′，進(jìn)而計(jì)算出速度修正量u′.之后利用關(guān)系式u＊＊＝u＊＋u′和p＊＊＝p＊＋p′，求解出速度與壓力的修正值u＊＊和p＊＊.在第二校正步中，將速度與壓力的修正值代入第二次壓力修正方程得：

通過(guò)式（9）得出第二次壓力修正量p″，進(jìn)而計(jì)算出第二次速度修正量u″，然后根據(jù)關(guān)系式u＊＊＊＝u＊＊＋u″和p＊＊＊＝p＊＊＋p″，求解出改進(jìn)的速度值和壓力值.

3 計(jì)算結(jié)果與分析

采用圖1所示的圓盤空化器多錐體射彈模型，進(jìn)行了自然超空泡流動(dòng)數(shù)值計(jì)算.

圖1 圓盤空化器多錐體超空泡射彈模型

在圖1中，空化器直徑為3.5mm，射彈長(zhǎng)徑比為1∶10.計(jì)算域采用圖2（a）所示的分塊結(jié)構(gòu)化網(wǎng)格對(duì)超空泡射彈進(jìn)行三維網(wǎng)格建模，網(wǎng)格總數(shù)為1 602 000.為了保證1＜y＋＜3，第一層網(wǎng)格點(diǎn)距離壁面為6.5μm.在非定常計(jì)算中，為獲得較好計(jì)算結(jié)果，通常以較小時(shí)間步長(zhǎng)計(jì)算，文中計(jì)算時(shí)間步長(zhǎng)為4.3ns，而且當(dāng)流場(chǎng)變量的殘差降到5×10－4以下時(shí)，認(rèn)為計(jì)算過(guò)程收斂.由于網(wǎng)格數(shù)較大及時(shí)間步長(zhǎng)較小，采用MPI并行算法提高計(jì)算效率.

超空泡流場(chǎng)計(jì)算邊界條件如圖2（b）所示.計(jì)算域頂部、底部與左端均為速度入口，右端為壓力出口，射彈表面采用固壁邊界條件.關(guān)于物性參數(shù)的選取，水的密度為1 000kg/m3，動(dòng)力粘性系數(shù)為1×10－6m2/s；水蒸氣的密度為23.08g/m3，動(dòng)力粘性系數(shù)為4.27×10－4m2/s.定義無(wú)量綱量空化數(shù)σ＝（p∞－pv）/（0.5ρu2∞），其中，p∞為來(lái)流壓力，pv為空化壓力，ρ為水的密度，u∞為來(lái)流速度.在自然超空泡條件下，本文通過(guò)變化來(lái)流速度來(lái)調(diào)整空化數(shù).

圖2 計(jì)算網(wǎng)格與邊界條件

3.1 超空泡形成過(guò)程

在空化數(shù)σ＝0.002 17的條件下，計(jì)算了繞超空泡射彈的空泡流動(dòng)，并給出了不同時(shí)刻下超空泡形態(tài)的計(jì)算結(jié)果，如圖3所示.

由圖3可以看出，在圖3（a）～圖3（d）階段，射彈頭部空泡與尾部空泡隨時(shí)間的增加而不斷變大、變長(zhǎng)；在圖3（d）～圖3（g）階段，射彈頭部形成的空泡與尾部空泡相連接，并繼續(xù)向遠(yuǎn)處發(fā)展，形成完全包裹彈身的超空泡；在圖3（h）階段，在射彈表面出現(xiàn)穩(wěn)定的空化區(qū)域，尾部被一個(gè)拖拽的巨大氣泡包圍著，形成一種穩(wěn)態(tài)的開(kāi)式尾流超空泡流型.

3.2 空化數(shù)對(duì)超空泡形態(tài)的影響

定義超空泡最大長(zhǎng)度與空化器直徑的比值c＝Lc/Dn，超空泡最大直徑與空化器直徑的比值c＝Dc/Dn.

為研究空化數(shù)對(duì)超空泡形態(tài)特性的影響，在空化數(shù)分別為4.88×10－3，3.12×10－3，2.17×10－3，1.59×10－3，1.22×10－3，0.96×10－3，0.78×10－3的條件下，對(duì)無(wú)量綱化空泡特征長(zhǎng)度c和無(wú)量綱化空泡特征直徑c進(jìn)行了計(jì)算，計(jì)算結(jié)果如圖4所示.

圖4 超空泡特征長(zhǎng)度與直徑隨空化數(shù)變化關(guān)系

由圖4可以看出，超空泡特征長(zhǎng)度與特征直徑隨空化數(shù)的增大而減小.為了驗(yàn)證計(jì)算結(jié)果的正確性，將其分別與LONGVINOVICH經(jīng)驗(yàn)公式計(jì)算結(jié)果及SAVCHENKO經(jīng)驗(yàn)公式計(jì)算結(jié)果進(jìn)行對(duì)比，結(jié)果如圖4所示.由對(duì)比可知，關(guān)于特征直徑，本文數(shù)值計(jì)算結(jié)果與2種經(jīng)驗(yàn)公式吻合較好；關(guān)于特征長(zhǎng)度，本文計(jì)算結(jié)果與2種經(jīng)驗(yàn)公式變化趨勢(shì)大體一致，但本文計(jì)算結(jié)果與LONGVINOVICH[10]經(jīng)驗(yàn)公式吻合較好，而與SAVCHENKO[11]經(jīng)驗(yàn)公式的吻合較差，這是由于SAVCHENKO經(jīng)驗(yàn)公式主要是適用于空化數(shù)在0.012～0.057的范圍內(nèi)，而LONGVINOVICH經(jīng)驗(yàn)公式是適用于空化數(shù)在0～0.25的范圍內(nèi).

3.3 攻角對(duì)超空泡形態(tài)的影響

在空化數(shù)為0.002 17的條件下，計(jì)算了攻角α分別為1°，3°，5°時(shí)水下超空泡射彈的空泡流動(dòng)，計(jì)算結(jié)果如圖5所示.

圖5 不同攻角下超空泡射彈的空泡形態(tài)

從圖5中可以看出，在有攻角的條件下，射彈表面上、下空泡厚度出現(xiàn)不對(duì)稱性，且隨著攻角的增加，其不對(duì)稱性增強(qiáng).為了分析攻角對(duì)空泡形態(tài)的影響規(guī)律，給出了不同攻角下的射彈上、下表面空泡厚度分布曲線，如圖6所示，圖中hu，hd分別為上、下表面空泡厚度.

從圖6（a）可以看出，在攻角為0°和1°的情況下，射彈上表面的空泡厚度從射彈前緣起先逐漸增加，隨后變小，再逐漸變大；在攻角為3°和5°的情況下，射彈上表面的空泡厚度從射彈前緣起逐漸變大；另外隨著攻角的增加，射彈上表面的空泡厚度不斷變大.從圖6（b）可以看出，在攻角為0°和1°的情況下，射彈下表面的空泡厚度從射彈前緣起先逐漸增加，隨后逐漸變小，再逐漸變大；在攻角為3°和5°的情況下，射彈下表面的空泡厚度由前緣起略有增大后逐漸變小，甚至在射彈的下表面，彈體有部分處于粘濕狀態(tài)，之后再略有增加.綜合上述分析可知，攻角對(duì)于超空泡射彈的空泡形態(tài)影響很大，特別是當(dāng)攻角較大時(shí)，射彈有部分處于粘濕狀態(tài)，這將影響超空泡對(duì)于射彈的減阻效果及射彈的水下彈道.

圖6 不同攻角下射彈上、下表面空泡厚度分布曲線

4 結(jié)論

本文基于PISO算法及N－S方程，利用LES湍流模型和VOF方法，實(shí)現(xiàn)三維非定常水下超空泡射彈形態(tài)特性的數(shù)值模擬，得出如下結(jié)論：

①通過(guò)分析超空泡形成過(guò)程以及對(duì)比計(jì)算結(jié)果與經(jīng)驗(yàn)公式，驗(yàn)證了本文計(jì)算結(jié)果的正確性；

②研究了不同空化數(shù)對(duì)超空泡形態(tài)的影響，計(jì)算結(jié)果表明，當(dāng)空化器直徑一定時(shí)，超空泡特征長(zhǎng)度和特征直徑隨空化數(shù)的增加而減小；

③射彈航行攻角對(duì)水下超空泡形態(tài)特性影響較大，其顯著改變超空泡形態(tài)的軸對(duì)稱性.

[1]FAROUK M，NAYFEH H.A compressible multi－phase flow solver for the computation of the supercavitation over highspeed torpedo[J].AIAA，2002，14（17）：1－14.

[2]LOGVINOVICH G V.Some problems of supercavitating flows[C].Proceedings of NATO－AGARD.Ukraine：NAS－IHM，1997：1－5.

[3]SAVCHENKO Y N，SEMENENKO V N.The gas absorption into supercavity from liquid－gas bubble mixture[C].Proceedings of 3rd International Symposium on Cavitation.Grenoble，F(xiàn)rance：Universite Joseph Fourier，1998：49－53.

[4]KINNAS S A，F(xiàn)INE N E.A numerical nonlinear analysis of the flow around two－and three－dimensional partially cavitating hydrofoils[J].J Fluid Mech，1993，254：151－181.

[5]DANG J，KUIPER G.Re－entrant jet modeling of partial cavity flow on three－dimensional hydrofoils[J].Journal of Fluids Engineering，1999，121：781－787.

[6]KUNZ R F，BOGER D A，STINEBRING D R，et al.A preconditioned Navier－Stokes method for two－phase flows with application to cavitation prediction[J].Computers ＆ Fluids，2000，29：849－875.

[7]陳瑛，魯傳敬，吳磊.三維小空化數(shù)空泡流數(shù)值方法[J].計(jì)算物理，2008，25（2）：163－171.CHEN Ying，LU Chuan－jing，WU Lei.Numerical method for three－dimensional cavitation flows at small cavitation numbers[J].Chinese Journal of Computational Physics，2008，25（2）：163－171.（in Chinese）

[8]賈力平，張嘉鐘，于開(kāi)平.空化器線性與超空泡減阻效果關(guān)系研究[J].船舶工程，2006，28（2）：20－23.JIA Li－ping，ZHANG Jia－zhong，YU Kai－ping.Research on the relation of cavitator profile and supercavity drag reduction[J].Ship Engineering，2006，28（2）：20－23.（in Chinese）

[9]QIAO Qin.A virtual single－phase natural cavitation model and its application to CAV2003 hydrofoil[C].5th International Symposium on Cavitation.Japan：Nova Science，2003：1－8.

[10]LOGVINOVICH G V.Hydrodynamics of flow with free boundaries[M].Kiev：Naukova Dumka，1969.

[11]SAVCHENKO Y N，VLASENKO Y D，SEMENENKO V N.Experimental study of high－speed cavitation flow[J].Hydromechanics，1998，（72）：103－111.