馬慶鵬 魏英杰 王 聰 曹 偉

(哈爾濱工業(yè)大學(xué) 航天學(xué)院，哈爾濱 150001)

航行體以一定的初速度穿越氣水界面進(jìn)入水中這一過程稱之為入水過程.關(guān)于入水問題的研究，最早可以追溯到19世紀(jì)末，A.M.Worthington和 R.S.Cole［1］利用高速閃光相機(jī)、玻璃水槽等設(shè)備開展了球體垂直入水實(shí)驗(yàn)，闡述并分析了球體入水各個(gè)階段的流動(dòng)現(xiàn)象.二戰(zhàn)后，美國海軍基于對(duì)空投魚雷研制的需求，召集諸多學(xué)者開展了大量的入水實(shí)驗(yàn)研究.美國海軍軍械實(shí)驗(yàn)室(NOL)的 D.Gilbarg 和 R.A.Anderson［2］通過實(shí)驗(yàn)研究了大氣壓力對(duì)球體入水噴濺和空泡的影響規(guī)律，分析了弗勞德數(shù)在入水空泡和噴濺模擬中的影響.A.May 和 J.C.Woodhull［3－5］在此基礎(chǔ)上開展了大量的球體及錐體、圓盤等軸對(duì)稱航行體的入水實(shí)驗(yàn)，得到了入水空泡的發(fā)展規(guī)律以及入水過程的流體動(dòng)力、入水彈道以及阻力系數(shù)等實(shí)驗(yàn)數(shù)據(jù).公開文獻(xiàn)顯示，這一階段的入水實(shí)驗(yàn)大多以低速(小于空氣中聲速)為主.在高速入水方面，1946年普林斯頓大學(xué)的 E.N.Harvey和 J.H.Mc-Millen［6］在美國國防研究委員會(huì)的報(bào)告中研究了鋼制球體在約243～914m/s速度范圍內(nèi)的入水問題，并通過實(shí)驗(yàn)觀察到了入水時(shí)產(chǎn)生的彈性激波.在此之后，公開文獻(xiàn)沒有400m/s以上入水速度的實(shí)驗(yàn)開展.在理論研究方面，M.Lee等人［7］從能量守恒的角度，建立了結(jié)構(gòu)入水過程能量轉(zhuǎn)換方程，得到了空泡半徑及空泡壁面擴(kuò)張速度的表達(dá)式，在此基礎(chǔ)上求解了球體以500m/s初始速度入水的空泡形態(tài)及擴(kuò)張速度.

國內(nèi)，何春濤、王聰?shù)龋?－9］開展了錐頭圓柱體在2～56m/s范圍內(nèi)垂直、傾斜入水的模型實(shí)驗(yàn)，研究了入水空泡生成、發(fā)展和閉合流動(dòng)過程以及空泡穩(wěn)定性，并分析了初始速度對(duì)空泡閉合方式和形態(tài)的影響.同時(shí)對(duì)低速入水問題開展了相關(guān)的數(shù)值模擬研究，分析了空泡內(nèi)的壓強(qiáng)分布及對(duì)空泡形態(tài)、閉合方式的影響.

綜上所述，國內(nèi)外對(duì)于入水問題的實(shí)驗(yàn)及數(shù)值模擬研究大多為低速狀態(tài)(低于空氣中聲速)下，對(duì)于高速入水的研究很少.本文采用有限體積法和流體體積法(VOF，Volume of Fluid)多相流模型，考慮入水過程的空化及湍流現(xiàn)象，對(duì)半錐角為63.5°的錐頭圓柱體高速狀態(tài)下垂直自由入水問題開展了數(shù)值計(jì)算研究，并將計(jì)算結(jié)果與文獻(xiàn)中基于能量守恒的求解結(jié)果進(jìn)行對(duì)比，驗(yàn)證了本文數(shù)值計(jì)算方法的合理性.進(jìn)一步開展不同初始速度下的數(shù)值計(jì)算，分析了入水空泡的生成、發(fā)展過程及流場(chǎng)壓力分布規(guī)律.

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

1.1 流體動(dòng)力學(xué)控制方程

本文數(shù)值計(jì)算假設(shè)流體介質(zhì)為不可壓縮，且忽略流體黏性產(chǎn)生的熱傳導(dǎo)效應(yīng)，即不求解能量方程，同時(shí)不考慮航行體在空氣中以超聲速運(yùn)動(dòng)產(chǎn)生的脫體激波對(duì)航行體姿態(tài)及水域的影響.VOF多相流模型將多相流體看作單一的流體介質(zhì)混合物，針對(duì)本文所求解的水、空氣及水蒸氣三相流動(dòng)問題，分別用 αl，αg，αv表示水、空氣和水蒸氣的體積分?jǐn)?shù)，在整個(gè)流場(chǎng)的任一區(qū)域，三者滿足關(guān)系式:αl+αg+αv=1.

混合流體介質(zhì)的連續(xù)性方程:

混合相的動(dòng)量守恒方程為

其中，ui為速度分量;ρm=(1 － αg－ αv)ρl+αgρg+αvρv代表混合介質(zhì)密度;μm=(1 － αg－αv)μl+ αgμg+ αvμv代表混合介質(zhì)的動(dòng)力黏性系數(shù);μt= ρmCμk2/ε 為湍流粘性系數(shù).

1.2 空化及湍流模型

空化是入水問題最重要的流動(dòng)現(xiàn)象之一，本文采用Schnerr and Sauer空化模型計(jì)算液態(tài)水的空化問題，該模型的守恒方程基于水蒸氣相建立，水蒸氣相的體積控制方程為

其中，RB=1 ×10－6m為氣核半徑;αnuc=5 ×10－4為不可凝結(jié)氣體體積分?jǐn)?shù);Fvap=50和Fcond=0.001為經(jīng)驗(yàn)常數(shù).

湍流是流體的一種隨機(jī)的非線性運(yùn)動(dòng)，也是空化流場(chǎng)的主要特征之一，對(duì)于高速入水所涉及的這種強(qiáng)非線性的流動(dòng)問題，本文采用基于k-ω的SST(Shear Stress Transport)湍流模型進(jìn)行模擬計(jì)算.SST模型綜合了常見的k-ε模型邊界層外部獨(dú)立性和k-ω模型近壁穩(wěn)定性的優(yōu)點(diǎn).

2 數(shù)值計(jì)算

2.1 計(jì)算模型及邊界條件

由于錐頭圓柱體為軸對(duì)稱體，所以本文數(shù)值計(jì)算采用二維軸對(duì)稱模型，航行體半錐角63.5°，柱段直徑D=10mm，柱段長度L=50mm.計(jì)算域示意圖及航行體附近網(wǎng)格劃分如圖1所示.

圖1 計(jì)算域及模型表面網(wǎng)格劃分示意圖

計(jì)算域以x軸為對(duì)稱軸，正向?yàn)榱鲌?chǎng)深度增加方向，y軸為流場(chǎng)寬度方向，坐標(biāo)原點(diǎn)O取航行體頂部所在初始位置.流場(chǎng)空氣域高度為32D，流場(chǎng)直徑30D，初始時(shí)刻錐角頂端距水面距離為2D.航行體徑向6D范圍內(nèi)采用三角形網(wǎng)格加密，外場(chǎng)采用四邊形網(wǎng)格.其中加密區(qū)域近壁面第1層網(wǎng)格高度為D/1000，最外層網(wǎng)格高度為D/40，網(wǎng)格總數(shù)量為589752.

2.2 數(shù)值計(jì)算方法

基于VOF多相流模型采用有限體積法離散、求解流體控制方程，對(duì)運(yùn)動(dòng)航行體垂直自由入水空泡的發(fā)展過程進(jìn)行數(shù)值模擬.壓力場(chǎng)與速度場(chǎng)的耦合求解選用PISO算法;壓力場(chǎng)的空間離散采用PRESTO!格式;各相體積率離散采用CICSAM格式;綜合考慮收斂性與計(jì)算時(shí)間，對(duì)動(dòng)量方程的離散采用一階迎風(fēng)格式.計(jì)算過程中，引入動(dòng)網(wǎng)格技術(shù)，利用C++語言編制UDF定義運(yùn)動(dòng)區(qū)域的運(yùn)動(dòng)參數(shù)，網(wǎng)格更新方法采用動(dòng)態(tài)層法.為保證求解精度，更新網(wǎng)格高度均為D/40.

3 數(shù)值計(jì)算結(jié)果驗(yàn)證

3.1 基于能量守恒定律的求解方法

圖2為入水空泡發(fā)展的示意圖.文獻(xiàn)［7］從能量守恒的角度給出了一種計(jì)算入水空泡形態(tài)發(fā)展規(guī)律的方法.

圖2 入水空泡示意圖

考慮不可壓縮流體，忽略入水過程中的熱效應(yīng)，根據(jù)牛頓第二定律，可以得到

式中，A0為航行體截面積;Vp為入水速度;ρl為水的密度;)為空化數(shù)，其中，p0為參考?jí)簭?qiáng)，pv為空化壓強(qiáng);Cdx為阻力系數(shù)，其數(shù)值采用Sedov［10］得到的公式Cdx=Cd0+σ來確定，Cd0為空化數(shù)等于0時(shí)的阻力系數(shù).在較低速入水情況下，重力對(duì)入水空泡的影響是不能忽略的，但是對(duì)于本文計(jì)算的500m/s左右的速度區(qū)間，可以忽略運(yùn)動(dòng)航行體重力的影響，因此，式(4)可以簡化為

在入水過程中，運(yùn)動(dòng)航行體動(dòng)能的損失率可以表達(dá)為

式中，β=ρlA0Cdx/(2m);xb為航行體頭部所在位置.

空泡壁面的徑向速度可以表示為軸線處點(diǎn)源源強(qiáng)的表達(dá)式:

式中ξ為入水深度.

取軸向微小位移d x，排開水的動(dòng)能為

將式(7)代入式(8)可以得到

式中N=ln(Ω/a)，其取值一般在15～30之間.文獻(xiàn)［11］對(duì)于亞音速入水問題研究采用N=15求解，考慮到本文所研究的問題初始入水速度較高，流場(chǎng)橫向擾動(dòng)較大，因此采用N=30求解.

空泡內(nèi)所積蓄的壓力勢(shì)能可以表示為

令pg=p0(x)－pc(x)，由于空泡內(nèi)部為水蒸氣和低壓氣體的混合物，因此，這里假設(shè) pc(x)=3540Pa.

綜上所述，根據(jù)能量守恒原理，可以得到

綜合式(12)、式(13)，并考慮到當(dāng)航行體在tb時(shí)刻到達(dá)深度xb時(shí)，空泡半徑為D/2，可以得到空泡半徑的表達(dá)式:

在該表達(dá)式右側(cè)，只有一個(gè)未知量Vp，且可以通過求解微分方程(5)得到，因此，對(duì)于給定的時(shí)間t，便可以得到該時(shí)刻各深度的空泡半徑，當(dāng)d t足夠小時(shí)，根據(jù)各深度空泡半徑便可以描繪出空泡輪廓.

3.2 數(shù)值計(jì)算結(jié)果驗(yàn)證分析

針對(duì)本文計(jì)算模型，由式(5)、式(14)分別得到航行體以Vp=500m/s初始速度垂直入水的速度衰減曲線及空泡擴(kuò)張過程.其中，對(duì)于半錐角63.5°的 錐 體，Cd0=0.637［12］;水 的 密 度 為998.2 kg/m3;航行體材料為合金鋁，密度為2700 kg/m3.采用相同的邊界條件開展數(shù)值計(jì)算.

圖3為2種方法得到的入水過程速度衰減曲線，可以看出，數(shù)值計(jì)算結(jié)果與理論解具有較高的一致性.

圖3 入水速度變化數(shù)值計(jì)算與理論結(jié)果對(duì)比

取入水深度分別為2D，4D，6D和10D這4種情況下的空泡輪廓，如圖4所示.

圖4 入水空泡輪廓數(shù)值計(jì)算與理論結(jié)果對(duì)比

由圖4可以看出，2種方法得到的空泡輪廓基本一致，具有較好的一致性.但是由于理論計(jì)算沒有考慮入水過程的液面波動(dòng)及入水噴濺等因素，因此在自由液面以上，二者有較為明顯的差別.此外，從圖4中可以看出數(shù)值計(jì)算得到的空泡半徑略小于理論解，這是由于理論方法沒有考慮水橫向擴(kuò)張時(shí)的黏性阻力.進(jìn)一步分析在自由液面處2種算法得到的空泡半徑誤差，如表1所示.

表1 自由液面處空泡半徑誤差

由表1可知，在航行體入水初期，空泡半徑誤差較大，隨著入水深度的增加，誤差逐漸減小，這是因?yàn)樵谌胨跗冢鲌?chǎng)的運(yùn)動(dòng)較為復(fù)雜，尤其是自由液面處會(huì)產(chǎn)生水面抬升及噴濺等現(xiàn)象，這一時(shí)期流體與結(jié)構(gòu)之間互相作用力也尤為復(fù)雜.另一方面，在實(shí)際入水過程中，入水噴濺在表面張力、大氣壓力等因素的作用下會(huì)呈現(xiàn)向軸線收縮的趨勢(shì)，進(jìn)而帶動(dòng)自由液面處空泡口的收縮，而理論方法并未考慮這一點(diǎn)，因此在空泡口處的誤差要明顯大于遠(yuǎn)離自由液面處.

通過以上分析和對(duì)比驗(yàn)證，可以看出本文采用的數(shù)值計(jì)算方法能夠較好地模擬入水運(yùn)動(dòng)參數(shù)及空泡的發(fā)展過程，其計(jì)算結(jié)果是可信的.在此基礎(chǔ)上，開展不同初始速度條件下，航行體自由垂直入水空泡形態(tài)發(fā)展規(guī)律的數(shù)值計(jì)算研究.

4 高速入水空泡發(fā)展規(guī)律研究

4.1 高速入水空泡發(fā)展分析

對(duì)錐頭圓柱體以600m/s初始速度自由垂直入水問題開展數(shù)值模擬計(jì)算，分析其運(yùn)動(dòng)參數(shù)、空泡形態(tài)、入水空化及流體動(dòng)力等發(fā)展規(guī)律.

圖5給出了航行體入水后的速度及入水深度的變化規(guī)律.可以看出，其速度由600m/s衰減到約150m/s僅運(yùn)動(dòng)2ms，在這段時(shí)間內(nèi)，航行體共在水中前行約55D的距離.這表明，航行體在入水初期受到阻力非常大.

圖5 600m/s初始速度入水運(yùn)動(dòng)參數(shù)變化曲線

提取不同時(shí)刻對(duì)稱軸及航行體表面的壓力，取全局坐標(biāo)系，得到不同時(shí)刻的壓力曲線，如圖6所示.

由圖6可以看出，壓力峰值出現(xiàn)在航行體頭部，不同時(shí)刻壓力峰值的推移表征了航行體入水深度的變化.可以看出，在航行體觸水后，頭部流場(chǎng)壓力可達(dá)到大氣壓的千倍量級(jí)，且在頭部頂點(diǎn)處最大.隨著入水深度的增大，航行體速度迅速降低，壓力峰值也逐漸下降，但在2ms內(nèi)依然保持較高的水平.這表明在高速入水的前期，航行體將持續(xù)受到較高沖擊載荷的作用，這對(duì)入水結(jié)構(gòu)安全性的設(shè)計(jì)提出了較高的要求.

圖6 對(duì)稱軸及航行體表面壓力分布曲線

入水空泡的形成與尾端流場(chǎng)的壓力密切相關(guān)，根據(jù)伯努利方程，當(dāng)航行體高速穿過水域，在排開水的同時(shí)亦在接觸面附近形成低壓區(qū)，由此產(chǎn)生空化.提取上述時(shí)刻尾端流場(chǎng)對(duì)稱軸上的壓力，如圖7所示，其中在0.1ms時(shí)刻，航行體尾端尚未到達(dá)自由液面處.

圖7 航行體尾端流場(chǎng)壓力分布曲線

圖8給出了空泡發(fā)展過程中水蒸氣相和空氣相的體積分?jǐn)?shù)云圖.

圖8 入水過程水蒸氣及空氣相分布體積分?jǐn)?shù)云圖

結(jié)合圖7、圖8，由壓力曲線可以看出，在航行體完全入水前，尾端空氣域壓力有一個(gè)明顯的先上升后下降的過程，這是由于航行體在空氣中超音速運(yùn)動(dòng)時(shí)在尾部流場(chǎng)產(chǎn)生的壓力波動(dòng).由圖8可以看出，在入水初期，空泡內(nèi)即出現(xiàn)空化現(xiàn)象，航行體排開前端水域后，在空泡壁及靠近自由液面位置產(chǎn)生了大量水蒸氣.在空泡擴(kuò)張階段，隨著入水空泡的生長，空泡內(nèi)水蒸氣不斷增多并與空泡內(nèi)空氣混合，阻止了空氣從空泡口向空泡前端的運(yùn)動(dòng)，此時(shí)空泡前端主要由水蒸氣和少量的空氣混合組成，空泡內(nèi)部壓力迅速降為飽和蒸氣壓，并產(chǎn)生了零散分布的壓力峰值.

4.2 初始速度影響分析

考慮航行體以V1=400m/s，V2=500m/s以及V3=600m/s3種初始速度垂直自由入水，開展入水空泡發(fā)展規(guī)律的數(shù)值模擬研究.

圖9給出了不同初始速度入水過程中，空泡最大無量綱直徑(Dc/D)隨時(shí)間的變化規(guī)律.

圖9 不同入水速度下最大空泡直徑變化曲線

由圖9可以看出，速度較高時(shí)，最大空泡直徑更大，隨著航行體速度的降低，空泡直徑的擴(kuò)張速度逐漸放緩.這一規(guī)律完全符合能量守恒原理，在速度較高時(shí)，周圍流場(chǎng)可以獲得較大的動(dòng)能向四周運(yùn)動(dòng)，因此高速時(shí)的空泡直徑更大;流體周向運(yùn)動(dòng)的過程中，動(dòng)能逐漸轉(zhuǎn)化為附近流體介質(zhì)的壓力勢(shì)能，其擴(kuò)張速度也就隨之下降，當(dāng)擴(kuò)張速度降為零時(shí)便開始反向運(yùn)動(dòng)，進(jìn)而導(dǎo)致空泡的閉合或潰滅.

進(jìn)一步分析3種入水時(shí)刻最大空泡直徑的相對(duì)增長率，如表2所示.

表2 不同時(shí)刻最大空泡直徑相對(duì)增長率

由表2可以看出，各時(shí)刻最大空泡直徑的增長率都小于速度增長率，這表明流體黏性對(duì)空泡擴(kuò)張有著重要的影響.

圖10給出了入水時(shí)刻為0.5ms時(shí)3種速度下的流場(chǎng)壓力曲線.可見，3種狀態(tài)下壓力曲線規(guī)律一致，速度越高，壓力峰值越大.同時(shí)在航行體肩部有明顯的壓力震蕩，可達(dá)到參考?jí)毫?0倍量級(jí)，如圖10b所示.這表明錐頭圓柱體入水過程中，肩部受到流場(chǎng)干擾較大，對(duì)于帶有一定攻角的入水問題，這一作用極有可能導(dǎo)致彈道的失穩(wěn).

圖10 不同速度入水狀態(tài)下對(duì)稱軸壓力分布曲線

5 結(jié)論

本文對(duì)帶有63.5°半錐角的錐頭圓柱體高速垂直自由入水問題開展了數(shù)值模擬研究.通過對(duì)初始速度分別為400，500，600m/s這3種工況下的數(shù)值計(jì)算結(jié)果分析得到以下結(jié)論:

1)得到了初始入水速度為500m/s的空泡形態(tài)發(fā)展規(guī)律，并與文獻(xiàn)理論計(jì)算結(jié)果進(jìn)行了對(duì)比，二者具有較好的一致性，驗(yàn)證了本文采用的數(shù)值計(jì)算方法的正確性;

2)初始速度600m/s條件下，在入水初期航行體頭部受到千倍大氣壓力量級(jí)的沖擊載荷作用，速度迅速下降，在2ms內(nèi)衰減至初始速度的25%，同時(shí)，強(qiáng)烈的沖擊載荷對(duì)航行體的結(jié)構(gòu)安全性設(shè)計(jì)提出了較高的要求;

3)航行體入水后排開周圍流體，形成入水空泡，在空泡擴(kuò)張階段，空泡直徑及長度隨入水深度的增大而增大;初始入水速度越高，相同時(shí)刻下入水空泡的最大直徑越大;

4)在空泡分離點(diǎn)附近，航行體肩部受到較強(qiáng)的隨機(jī)作用力，對(duì)于有攻角的結(jié)構(gòu)入水，極有可能導(dǎo)致其彈道失穩(wěn).

References)

［1］ Worthington AM，Cole R S.Impact with a liquid surface studied by the aid of instantaneous photography［J］.Philosophical Transactions of the Royal Society，1900，194(A):175 －200

［2］ Gilbarg D，Anderson R A.Influence of atmospheric pressure on the phenomena accompanying the entry of spheres into water［J］.Journal of Applied Physics，1948，9(2):127 － 139

［3］ May A，Woodhull JC.Drag coefficients of steel spheres entering water vertically［J］.Journal of Applied Physics，1948，19(12):1109－1121

［4］ May A，Woodhull JC.The virtual mass of a sphere entering water vertically［J］.Journal of Applied Physics，1950，21(12):1285－1289

［5］ May A.Effect of surface condition of a sphere on its water-entry cavity［J］.Journal of Applied Physics，1951，22(10):1219 －1222

［6］ National Defense Research Committee.Mathematical studies relating to military physical research［R］.AD221604，1946

［7］ Lee M，Longoria R G，Wilson D E.Cavity dynamics in high speed water entry［J］.Physics of Fluids，1997，9(3):540 －550

［8］何春濤，王聰，魏英杰，等.圓柱體垂直入水空泡形態(tài)試驗(yàn)研究［J］.北京航空航天大學(xué)學(xué)報(bào)，2012，38(11):1542 －1546 He Chuntao，Wang Cong，Wei Yingjie，et al.Vertical water entry cavity of cylinder body［J］.Journal of Beijing University of Aeronautics and Astronautics，2012，38(11):1542 － 1546(in Chinese)

［9］王聰，何春濤，權(quán)曉波，等.空氣壓強(qiáng)對(duì)垂直入水空泡影響的數(shù)值研究［J］.哈爾濱工業(yè)大學(xué)學(xué)報(bào)，2012，44(5):14－19 Wang Cong，He Chuntao，Quan Xiaobo，et al.Numerical simulation of the influence of atmospheric pressure on water-cavity formed by cylinder with vertical water-entry［J］.Journal of Harbin Institute of Technolody，2012，44(5):14 －19(in Chinese)

［10］ Sedov L I.Two-dimensional problems in hydrodynamics and aerodynamics［M］.New York:John Wiley ＆ Sons Inc，1965:1－427

［11］ Lundstrom E A.Fluid dynamic analysis of hydraulic ram［R］.NWC TP 5227，1971

［12］ May A.Water entry and the cavity-running behavior of missiles［R］.AD A020429，1975