黃海龍，王柏秋

(1.哈爾濱工業(yè)大學(xué)土木工程學(xué)院，哈爾濱150090;2.北京宇航系統(tǒng)工程研究所，北京100076;3.哈爾濱工業(yè)大學(xué)建筑學(xué)院，哈爾濱150001)

導(dǎo)彈等物體從空中高速進(jìn)入水中的現(xiàn)象稱(chēng)為入水問(wèn)題，關(guān)于入水問(wèn)題的研究，最早可以追溯到19 世紀(jì)末，A·M·Worthington[1－2]利用高速照相機(jī)對(duì)液滴和球體入水后所產(chǎn)生的流動(dòng)現(xiàn)象進(jìn)行了試驗(yàn)觀察，由此開(kāi)始，各種針對(duì)入水空泡問(wèn)題的探索不斷的展開(kāi).從文獻(xiàn)[3－4]中可以看出，在第二次世界大戰(zhàn)前后，多個(gè)國(guó)家相繼利用試驗(yàn)方法和數(shù)值方法開(kāi)展了針對(duì)各種導(dǎo)彈、空投魚(yú)雷等武器的入水問(wèn)題的研究.其中，關(guān)于球體入水流場(chǎng)的試驗(yàn)研究[5－9]和基于勢(shì)流理論研究[10－11]是較多的，其次是關(guān)于圓柱體等平頭物體[12－14]的研究以及具有圓錐頭的柱體模型[15]和射彈模型[16]入水流場(chǎng)的分析.

以上試驗(yàn)研究和理論研究，分別得到了多種航行體的入水阻力系數(shù)與各種無(wú)量綱參數(shù)之間的關(guān)系，以及模型彈道等信息.但是，這些試驗(yàn)研究并未討論模型后體參數(shù)的變化情況，例如錐柱體入水問(wèn)題中，柱體段長(zhǎng)度對(duì)模型入水空泡的影響.并且，由于研究方法和測(cè)試手段的局限性，這些研究成果只能反映入水空泡流場(chǎng)的部分參數(shù)信息.因此，隨著計(jì)算機(jī)技術(shù)的發(fā)展，利用數(shù)值計(jì)算，直接求解NS方程的方法逐漸得到了蓬勃發(fā)展[17－18]，這使得所論目標(biāo)更接近真實(shí)，能夠獲得的流場(chǎng)信息更豐富.

綜上所述，帶有后體的航行體入水空泡問(wèn)題研究成果較少，而具有后體的航行體入水產(chǎn)生的空泡流場(chǎng)與不帶后體又具有一定差異.因此，針對(duì)后體對(duì)入水空泡流場(chǎng)的影響問(wèn)題，本文利用動(dòng)網(wǎng)格技術(shù)，采用VOF多相流模型和有限體積法求解RANS方程，對(duì)具有不同柱段長(zhǎng)度的錐頭圓柱體勻速垂直入水過(guò)程進(jìn)行數(shù)值計(jì)算，分析模型勻速垂直入水過(guò)程中，空泡的生成和發(fā)展過(guò)程，以及有、無(wú)后體和后體長(zhǎng)度對(duì)入水空泡的影響.

1 控制方程

1.1 試驗(yàn)方法連續(xù)性方程和動(dòng)量方程

混合物連續(xù)性方程:

空氣相連續(xù)性方程:

水蒸氣相連續(xù)性方程:

混合相動(dòng)量方程:

其中:ui表示混合物速度分量，m/s;下標(biāo)m、l、g和v分別表示混合相、水、空氣和水蒸氣;ρm為混合物密度(kg/m3);μm是混合物動(dòng)力黏性系數(shù)，m2/s.

其中:ρm和μm分別為:

其式中:αg、αv和 －αg、－αv分別為空氣、水蒸氣和水的體積分?jǐn)?shù).

1.2 空化模型

高速入水物體在入水過(guò)程中產(chǎn)生空泡，僅有物體頭部與水接觸，在物體頭部邊緣的分離處會(huì)出現(xiàn)較大的壓強(qiáng)梯度變化，當(dāng)壓強(qiáng)降低到飽和蒸汽壓強(qiáng)之下時(shí)，可能會(huì)出現(xiàn)自然空化現(xiàn)象.本文將這一復(fù)雜過(guò)程通過(guò)式(7)、(8)中的源項(xiàng)來(lái)表示(即Zwart－Gerber－Belamri空化模型).源項(xiàng)表達(dá)式如下:

其中:Re、Rc表示水蒸氣產(chǎn)生和凝結(jié);RB=10－6m為氣泡半徑;αnuc=5×10－4，為不可凝結(jié)氣體體積分?jǐn)?shù);Fvap、Fcond為經(jīng)驗(yàn)常數(shù)，分別取 50 和 0.001.

1.2 湍流模型

本文計(jì)算中采用Realizable k－ε模型模擬湍流，該湍流模型中，湍動(dòng)能k和湍能耗散率ε的輸運(yùn)方程分別為:

其中:σk和σε分別為k和ε的Prandtl數(shù)，取值分別為 σk=1.0，σε=1.2;為由平均速度梯度所引起的湍動(dòng)能生成項(xiàng);經(jīng)驗(yàn)常數(shù)C1和C2取值分別為1.44 和1.9.

2 數(shù)值求解

2.1 計(jì)算域和邊界條件

本文選取130°圓錐頭型為研究對(duì)象，模型后體圓柱直徑 D=10 mm.根據(jù)A·M·Worthington[1]的研究，軸對(duì)稱(chēng)物體垂直入水所產(chǎn)生的空泡流場(chǎng)同樣具有軸對(duì)稱(chēng)性質(zhì)，因此，本文數(shù)值計(jì)算均采用二維軸對(duì)稱(chēng)計(jì)算模型，計(jì)算域設(shè)置和模型表面網(wǎng)格分布示意圖如圖1所示.圖1中，空氣域高度為60 D，水域深度為400 D，計(jì)算域徑向?yàn)?0 D，使用四邊形網(wǎng)格劃分計(jì)算域.不失一般性，在保證計(jì)算精度和正確性的基礎(chǔ)上，選取網(wǎng)格密度和計(jì)算量適中的網(wǎng)格分布作為本文數(shù)值模型.

圖1 圓錐及計(jì)算域分布示意圖

2.2 數(shù)值方法

本文采用有限體積法對(duì)控制方程進(jìn)行空間離散，應(yīng)用PISO方法對(duì)動(dòng)量方程方程和連續(xù)性方程聯(lián)立求解.計(jì)算過(guò)程中，對(duì)流項(xiàng)的離散采用二階迎風(fēng)格式;水蒸氣相的離散采用QUICK格式;湍流輸運(yùn)方程的離散采用二階迎風(fēng)格式;對(duì)時(shí)間的離散采用一階顯示格式.計(jì)算時(shí)間步長(zhǎng)的選取根據(jù)網(wǎng)格尺度、運(yùn)動(dòng)速度等因素綜合確定.通過(guò)動(dòng)網(wǎng)格技術(shù)實(shí)現(xiàn)模型的運(yùn)動(dòng)過(guò)程與計(jì)算域內(nèi)的網(wǎng)格更新.

由于本文計(jì)算采用結(jié)構(gòu)化網(wǎng)格，網(wǎng)格主要是被拉伸或被壓縮變形，因此，采用動(dòng)態(tài)層分裂法進(jìn)行網(wǎng)格更新.動(dòng)態(tài)層分裂法的基本思想是當(dāng)網(wǎng)格被拉伸到大于設(shè)定的理想最大尺度((1+αs)hideal)時(shí)，該網(wǎng)格分裂成兩層新的網(wǎng)格，當(dāng)網(wǎng)格被壓縮到小于設(shè)定的理想最小尺度(αchideal)時(shí)，網(wǎng)格與其相鄰的網(wǎng)格合并成新的網(wǎng)格，即:

引入動(dòng)態(tài)網(wǎng)格后，相當(dāng)于網(wǎng)格按一定速度作對(duì)流運(yùn)動(dòng)，在實(shí)際計(jì)算中應(yīng)減掉網(wǎng)格的對(duì)流效應(yīng).因而引入動(dòng)網(wǎng)格后控制體V對(duì)變量φ的統(tǒng)一形式的守恒方程可寫(xiě)為:

其中:ρm為混合物密度為速度矢量為動(dòng)網(wǎng)格的運(yùn)動(dòng)速度;Γ為擴(kuò)散系數(shù);Sφ為標(biāo)量φ的源項(xiàng);?V表示控制體積V的邊界.

3 數(shù)值方法對(duì)比與驗(yàn)證

為驗(yàn)證方法的正確性和可信性，首先對(duì)直徑為25.4 mm球體勻速垂直入水過(guò)程進(jìn)行數(shù)值計(jì)算，并將計(jì)算結(jié)果的空泡輪廓與Albert May[6]等人提出的勻速垂直入水理想空泡模型進(jìn)行對(duì)比如下.

在二維軸對(duì)稱(chēng)坐標(biāo)系下，垂直入水早期空泡形態(tài)的理想空泡模型表達(dá)式如式(14)所示，該入水空泡模型主要用于描述入水過(guò)程中，空泡閉合之前發(fā)展較為充分時(shí)期的外部輪廓.

其中:CD為入水物體的阻力系數(shù);d為入水物體的特征尺寸.

模型從空氣域勻速垂直入水，入水速度30 m/s，空氣域環(huán)境壓強(qiáng)為101 325 Pa，空氣、水、和水蒸氣的密度分別為 1.225、998.2、0.554 kg/m3.選取入水空泡充分發(fā)展時(shí)的空泡形態(tài)，與相同時(shí)刻下，由式(14)的理想空泡模型擬合結(jié)果進(jìn)行對(duì)比，結(jié)果如圖2所示，其中，阻力系數(shù)取CD=0.20.

圖2 無(wú)后體圓錐入水空泡輪廓計(jì)算結(jié)果與經(jīng)驗(yàn)公式對(duì)比

由圖2可以看出，數(shù)值計(jì)算結(jié)果與理想空泡模型擬合結(jié)果具有較好的一致性.在空泡輪廓兩端具有一定的誤差，產(chǎn)生該誤差的主要原因是由于理想空泡模型本身忽略了自由表面處和分離點(diǎn)處的空泡輪廓的不規(guī)則性.通過(guò)圖2的對(duì)比驗(yàn)證可見(jiàn)，本文所采用的數(shù)值計(jì)算方法是正確可信的，在此基礎(chǔ)之上，下面針對(duì)帶有不同尺度圓柱后體的圓錐模型，展開(kāi)勻速垂直入水空泡的數(shù)值計(jì)算.

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

計(jì)算模型外形如圖3所示，頂部錐角 a=130°，直徑 D=10mm，后體長(zhǎng)度 L分別為:0(圓錐)，0.5D，1.0D，2.0D 和 4.0D.為考察后體對(duì)入水空泡的影響，將各模型入水計(jì)算結(jié)果的某一時(shí)刻流場(chǎng)分布顯示于圖4中.

圖3 不同長(zhǎng)度后體模型

圖4 入水空泡內(nèi)部組分分布

有后體后，可以發(fā)現(xiàn)入水空泡的最大直徑隨著后體尺度有一定變化.圖5表示垂直入水空泡最大直徑隨后體長(zhǎng)度的變化過(guò)程，入水速度為200 m/s.從圖5可以看出，模型后體長(zhǎng)度較小時(shí)(例如圖5中的無(wú)后體模型和后體長(zhǎng)度為L(zhǎng)=0.5D模型)，入水空泡最大直徑略高于其余帶后體模型，隨著模型后體長(zhǎng)度的增加，入水空泡最大直徑逐漸降低，當(dāng)后體長(zhǎng)度增加到一定程度后，入水空泡最大直徑趨于穩(wěn)定，不再隨后體長(zhǎng)度變化而變化.

圖5 空泡最大直徑隨模型后體的變化

圖4中的入水空泡尾部的閉合情況清晰表明入水空泡閉合時(shí)間隨后體的出現(xiàn)而有所變化.圖6表示模型勻速垂直入水過(guò)程中，入水空泡尾部無(wú)量綱閉合時(shí)間隨后體長(zhǎng)度的變化情況，模型入水速度為200 m/s.

圖6 入水空泡閉合時(shí)間隨后體的變化

從圖6可以看出，模型后體長(zhǎng)度較小時(shí)，空泡閉合時(shí)間較長(zhǎng)，尤其當(dāng)后體長(zhǎng)度L=0時(shí)，閉合時(shí)間最長(zhǎng);當(dāng)模型后體長(zhǎng)度較大時(shí)，空泡閉合時(shí)間趨于一致，例如當(dāng)L＞2.0時(shí)，空泡閉合時(shí)間不再隨后體長(zhǎng)度變化而變化.由圖5、6可以看出，帶后體模型產(chǎn)生的入水過(guò)程，當(dāng)模型后體長(zhǎng)度大于某一數(shù)值時(shí)，流場(chǎng)是相近的.為考察穩(wěn)定的入水過(guò)程，下面以L=4.0D的模型為例，進(jìn)行分析.

圖7為長(zhǎng)后體模型(L=4.0D)勻速垂直入水過(guò)程中，水蒸氣和空氣在空泡內(nèi)部的空間分布的變化情況，其中，模型入水速度為200 m/s.從圖7可以看出，在入水初期，模型肩部的流動(dòng)分離點(diǎn)附近和空泡入口處的噴濺內(nèi)部邊緣有少量的水蒸氣.隨著入水過(guò)程的進(jìn)行，空泡入口的噴濺液體迅速向內(nèi)彎曲，空泡開(kāi)始趨于閉合.此時(shí)，在空泡入口處的彎曲液壁內(nèi)側(cè)，由于液面迅速?gòu)澢浇黧w流速增高，靜壓力急劇下降，導(dǎo)致液體迅速空化，從而使得該處水蒸氣含量隨著空泡閉合而迅速增加，隨著時(shí)間的推移，空化逐漸沿著空泡壁面向空泡內(nèi)部傳播.在模型附近，既空泡頭部附近，空泡輪廓線穩(wěn)定，不隨時(shí)間變化，因此，該處流場(chǎng)內(nèi)部壓力和流速均穩(wěn)定，發(fā)生空化的概率相對(duì)較低，該處只存在少量入水沖擊產(chǎn)生的水蒸氣或者不存在水蒸氣.

另外，由圖7的入水過(guò)程可知，入水空泡內(nèi)部流體壓力低于外界大氣壓和周?chē)黧w內(nèi)部靜壓，當(dāng)空泡壁運(yùn)動(dòng)能量不足時(shí)，將導(dǎo)致空泡壁向軸線坍縮，從而擠壓內(nèi)部空氣，最終可能導(dǎo)致一股沿空泡軸線向上的液體噴濺的產(chǎn)生.

圖7 不同時(shí)刻入水空泡組分分布

5 結(jié)論

本文對(duì)具有不同長(zhǎng)度圓柱型后體的圓錐勻速垂直入水過(guò)程進(jìn)行了數(shù)值模擬研究，分析了有、無(wú)后體和不同長(zhǎng)度后體模型入水后，空泡的發(fā)展過(guò)程及空泡內(nèi)部流場(chǎng)的分布，得到了以下結(jié)論:

1)模型在勻速垂直入水后形成的入水空泡內(nèi)，無(wú)后體模型流場(chǎng)內(nèi)的空化作用弱于有后體模型流場(chǎng)，入水空泡內(nèi)的水蒸氣含量明顯偏低.

2)帶圓柱后體模型，后體的存在有利于模型肩部流動(dòng)分離點(diǎn)附近自然空化的發(fā)生，并有少量的水蒸氣分布在分離點(diǎn)附近的空泡壁面.

3)無(wú)后體模型入水空泡的閉合明顯晚于帶圓柱后體模型入水空泡的閉合，當(dāng)圓柱后體長(zhǎng)度小于1倍圓柱體直徑時(shí)，空泡最大直徑隨著后體長(zhǎng)度的增加而減小，且空泡最大直徑略大于具有更長(zhǎng)帶后體的模型;當(dāng)后體長(zhǎng)度大于2倍以上圓柱體直徑時(shí)，入水空泡最大直徑與后體長(zhǎng)度無(wú)關(guān).

4)無(wú)后體時(shí)，入水空泡的閉合時(shí)間較長(zhǎng);有后體時(shí)，后體長(zhǎng)度大于兩倍圓柱直徑時(shí)，入水空泡的閉合時(shí)間與后體長(zhǎng)度無(wú)關(guān).

通過(guò)本文研究可以看出，模型后體的存在對(duì)入水空泡形態(tài)及相關(guān)參數(shù)有重要影響，為了使入水問(wèn)題研究結(jié)果在魚(yú)雷和導(dǎo)彈等武器中得到更好的應(yīng)用，應(yīng)加強(qiáng)對(duì)模型后體影響的研究.

[1]WORTHINGTON A M.A study of splashes[M].New York:Longmans Green and Company，1908.

[2]WORTHINGTON A M，COLE R S.Impact with a liquid surface[J].Phil.Trans.Roy.Soc.，1897，189A:137 －148.

[3]WHALLEY I A.Project upkeep－a review of the WWⅡdambubuster weapon[C]//Indianapolis:38th AIAA Joint propulsion Conf，2002.

[4]ALBERT M.Water entry and the cavity－running behavior of missiles[R].AD A020429，1975.

[5]MALLOCK A.Sounds produced by drops falling on water[J].Pro.R.Soc.Lond.，1918，A95:138 －143.

[6]BELL G E.On the impact of a solid sphere with a fluid surface

[J].Phil.Mag.j.Sci，1924，48:753 － 765.

[7]ASHLEY S.Warp drive underwater[J].Sci.Amer，2001，284(5):70－79.

[8]MAY A，WOODHULL J C.Drag coefficients of steel spheres entering water vertically[J].Journal of Applied Physics，1948，19(12):1109－1121.

[9]TRUSCOTT T T，TECHET A H.Water entry of spinning spheres[J].J.Fluid Mech，2009，625:135 －165.

[10]Yan H，Liu Y，KOMINIARCZUK J，et al.Cavity dynamics in water entry at low Froude numbers[J].J.Fluid Mech，2009，641:441－461

[11]ARISTOFF J M，BUSH J W M.Water entry of small hydrophobic spheres[J].J.Fluid Mech，2009，619:45 －78.

[12]SHI H，ITOH M，TAKAMI T.Optical observation of the supercavitation induced by high－speed water entry[J].Journal of Fluids Engineering，2000，122:806－810.

[13]何春濤，王 聰，何乾坤，等.圓柱體低速入水空泡試驗(yàn)研究[J].物理學(xué)報(bào)，2012，61(13):1－8.

[14]陳學(xué)農(nóng)，何友聲.平頭物體三維帶空泡入水的數(shù)值模擬[J].力學(xué)學(xué)報(bào)，1990，22(2):129－138.

[15]王 聰，何春濤，權(quán)曉波，等.空氣壓強(qiáng)對(duì)垂直入水空泡影響的數(shù)值研究[J].哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2012，44(5):14－19.

[16]MAY A.Effect of surface condition of a sphere on its waterentry cavity[J].Journal of Applied Physics，1951，22:1219－1222.

[17]HE C，WANG C，WEI Y，et al.Numerical simulation of pressure distribution in vertical water－ entry cavity[J].Journal of Ship Mechanics，2011，15(9):960 －968.

[18]何春濤，王 聰，閔景新，等.回轉(zhuǎn)體勻速垂直入水早期空泡數(shù)值模擬研究[J].工程力學(xué)，2012，29(4):237－243.