孫　凱，黨建軍，郝維敏，蔣　彬

（西北工業(yè)大學(xué) 航海學(xué)院，陜西 西安，710072）

回轉(zhuǎn)體超音速入水沖擊數(shù)值仿真

孫凱，黨建軍，郝維敏，蔣彬

（西北工業(yè)大學(xué) 航海學(xué)院，陜西 西安，710072）

入水沖擊問題是現(xiàn)代水中兵器的研究熱點之一。采用商用CFD軟件，使用流體體積函數(shù)（VOF）模型并結(jié)合動網(wǎng)格技術(shù)，仿真了平頭回轉(zhuǎn)體模型在亞音速、跨音速以及超音速狀態(tài)下的入水過程，得到了模型在不同速度下的入水阻力特性、頭部直徑、液體可壓縮性及空氣激波對入水過程的影響。仿真結(jié)果表明，入水阻力系數(shù)峰值隨著速度的增加先降低后平緩，峰值寬度逐漸降低； 增大頭部直徑，阻力系數(shù)峰值增加，峰值寬度增加； 水的壓縮性會使阻力系數(shù)峰值出現(xiàn)的時間延后，數(shù)值降低； 空氣激波會使接觸水面時刻阻力系數(shù)降低。

平頭回轉(zhuǎn)體； 入水過程； 流體體積函數(shù)（VOF）模型； 阻力系數(shù)

0　引言

高速射彈運用超空泡技術(shù)使得自身航行阻力減小一個數(shù)量級，可用于完成攔截魚雷、擊毀水雷、破除水下障礙和獵殺蛙人等攻擊和防御任務(wù)。高速射彈的入水過程是從空中彈道轉(zhuǎn)入水下彈道的過度環(huán)節(jié)，受到巨大的沖擊載荷，可能會造成彈道失穩(wěn)甚至結(jié)構(gòu)失效。因此，入水過程對射彈的整個運動過程非常重要。

國外關(guān)于入水理論的相關(guān)研究始于20世紀20～30年代。G. V. Logvinovich提出的基于能量守恒的空泡截面擴張過程是最早描述空泡發(fā)展過程的理論，并成為后續(xù)理論研究的重要參考［1］。國內(nèi)入水過程的理論和數(shù)值研究始于上20世紀90年代。王永虎對雷彈頭型與入水沖擊阻力特征進行了分析，得到頭型長細比對入水沖擊的影響［2］；陳宇翔等人使用流體體積函數(shù)（volume of fluid，VOF）方法研究水平圓柱砰擊入水問題，得到了湍流粘性對圓柱入水的影響［3］； 何春濤等人基于RNS方程和VOF模型對回轉(zhuǎn)體垂直勻速入水進行了數(shù)值仿真，得到了空泡形態(tài)隨時間的變化規(guī)律［4］； 邱海強等人使用混合模型對不同頭型回轉(zhuǎn)體入水過程進行了數(shù)值仿真，得到不同頭型對入水沖擊載荷和空泡形態(tài)的影響［5］。

目前，對于射彈入水問題的研究主要采用試驗和數(shù)值仿真2種方法。相較試驗研究成本高、周期長、結(jié)果不精確等缺點，數(shù)值仿真方法具有優(yōu)勢。這里采用計算流體力學(xué)（computational fluid dynamics，CFD）方法，仿真了平頭回轉(zhuǎn)體模型在亞音速、跨音速以及超音速狀態(tài)下的入水過程，分析了水可壓縮性和空氣可壓縮性對入水沖擊的影響。

1　模型建立與驗證

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

1）Tait狀態(tài)方程

當水受到較大的壓力時，需要考慮水的壓縮性影響，文章應(yīng)用的水可壓縮模型建立在Tait方程基礎(chǔ)上。Tait方程在仿真水下爆炸等方面應(yīng)用廣泛。不含溫度修正的Tait方程為

其中： p0是水的參考絕對壓力；0ρ是在參考壓力下水的參考密度； K0是在參考壓力下的參考體積彈性模量； n是密度指數(shù)； p是水受到的絕對壓力； ρ是在壓力為p時的密度。

為了使可壓縮流動中的壓力計算結(jié)果一定，需要在壓力修正方程中加入與聲速有關(guān)的項。利用可壓縮流體的聲速理論可得簡單的聲速方程

文中研究的工況涉及到跨超聲速入水，因此流體的材料屬性均使用可壓縮流體。

2）6自由度動網(wǎng)格模型

6自由度求解器在計算平移運動時加速度計算較精確，因而文中采用了該求解器模型。在慣性參考坐標系下建立剛體平移運動控制方程為

式中： v.G是質(zhì)心平移加速度矢量； m是剛體質(zhì)量； fG是作用在物體上的外力。文中回轉(zhuǎn)體運動通過UDF定義6自由度運動網(wǎng)格實現(xiàn)。

1.2計算模型簡化及邊界條件

建立Φ200 mm×1 000 mm的圓柱形計算域，為減小網(wǎng)格量，加快收斂速度，計算時選取流場的1/2。模型初始位置為水面上一定距離，在空氣中飛行一段時間整個流場穩(wěn)定后再與水面接觸。

經(jīng)典的高速射彈模型大多是長細比很大的回轉(zhuǎn)體，文中仿真的垂直入水工況主要研究的是入水沖擊的影響，從節(jié)省計算時間成本考慮，計算模型僅保留頭部空化器尺寸，后體采用等直徑圓柱代替，簡化后的計算模型為Φ3 mm×15 mm的平頭圓柱，并將其視為剛體?；剞D(zhuǎn)體頭部阻力系數(shù)Cd計算公式為

式中：FD為射彈入水過程中收到的阻力，F(xiàn)D= FD（t）； v為運動體入水過程中的瞬時速度，v= v（t）；ρ為水密度；A為運動體特征面積。計算流場為整個流場的1/2，該特征面積取3.532 5×10-6m2。

圖1表明了計算域各部分所采用的邊界條件：計算域下界采用壓力出口，上界采用壓力入口；遠場以及回轉(zhuǎn)體均為無滑移壁面； 計算域上部為空氣域，下部為水域； 回轉(zhuǎn)體運動使用動網(wǎng)格方法實現(xiàn)，靠近回轉(zhuǎn)體及其運動路徑附近的區(qū)域為運動區(qū)域，其余部分為靜止區(qū)域，兩部分采用interface邊界條件連接。

圖1　計算流場及邊界條件Fig. 1　Calculation of flow field and boundary conditions

Fluent仿真設(shè)置： 1）多相流模型選擇VOF模型，主相為水，副相為空氣和水蒸氣； 2）湍流模型選擇Realizable k-ε模型； 3）壓力速度場耦合方式選用PRESTO??； 4）水使用可壓縮流體模型； 5）開啟空化模型。

1.3模型驗證

為驗證上述計算模型的合理性，選取相關(guān)論文中的典型試驗工況進行數(shù)值仿真分析，評估其合理性。M.Schaffar等對圓柱體在1 036 m/s速度下入水進行了試驗研究，得到了該試驗工況下的流場分布［6］。

圖2給出了回轉(zhuǎn)體在入水25 cm深度后的流場分布，可以看出回轉(zhuǎn)體在水中形成了完全包裹自身的超空泡。圖3給出了回轉(zhuǎn)體入水后距離頭部不同位置處的空泡直徑，其中Serebryakov、Levinson是2種數(shù)值計算結(jié)果，Measurement是試驗得到的結(jié)果。從圖3中可以看出，文中使用的計算模型得到的結(jié)果與Serebryakov的數(shù)值仿真結(jié)果基本一致，與試驗的最大誤差不超過12.5%。

圖2　試驗與仿真結(jié)果流場對比圖Fig. 2　Comparison between experimental and simulated flow field

圖3　空泡外形對比圖Fig. 3　Comparison of cavitation outlines

2　仿真結(jié)果與分析

2.1回轉(zhuǎn)體高速入水流場分析

圖4和圖5顯示的是回轉(zhuǎn)體以1 000 m/s工況下不同位置處的水相云圖和蒸汽相云圖。

從圖4可看出，回轉(zhuǎn)體在入水過程中各個階段的空泡發(fā)展過程?；剞D(zhuǎn)體從距離水面20 mm處開始下落，接觸水面后，水面受到?jīng)_擊向四周飛濺； 隨著回轉(zhuǎn)體繼續(xù)下落，形成開口空泡并不斷發(fā)展； 在入水深度達到2倍回轉(zhuǎn)體長度后，空泡上端依然沒有閉合趨勢。

圖4　1 000 m/s速度入水階段水相體積分數(shù)云圖

圖5　1 000 m/s速度入水階段蒸汽相體積分數(shù)云圖Fig. 5　Contour of vapor volume fraction in water entry phase at a speed of 1 000 m/s

從圖5可以看出，在空泡內(nèi)氣相和水相的臨近區(qū)域產(chǎn)生了大量的水蒸氣。這是因為回轉(zhuǎn)體入水速度超過聲速，空氣不能及時補充進空泡內(nèi)，導(dǎo)致空泡內(nèi)的壓力低于空化壓力，產(chǎn)生了空化現(xiàn)象。正是由于空化現(xiàn)象，入水空泡內(nèi)氣相和水相的交界面變得“凹凸不平”。

2.2速度對沖擊阻力系數(shù)的影響

圖6顯示的是回轉(zhuǎn)體以1 000 m/s速度入水過程阻力系數(shù)在峰值附近的變化過程。在空中飛行階段，由于空氣的阻力很小，阻力系數(shù)很小； 回轉(zhuǎn)體接近水面，阻力系數(shù)迅速增大，在2×10-5s時回轉(zhuǎn)體與水面接觸，此時阻力系數(shù)為2.5； 回轉(zhuǎn)體繼續(xù)運動，阻力系數(shù)增大，在2.02×10-5s時阻力系數(shù)達到最大值3.23； 阻力系數(shù)在達到最大值后逐漸降低，整個沖擊峰值持續(xù)約0.12×10-5s。

圖6　1 000 m/s速度入水阻力系數(shù)峰值放大圖Fig. 6　Enlarged view of peak value of drag coefficient in water entry phase at a speed of 1 000 m/s

圖7 顯示的是回轉(zhuǎn)體以不同速度入水時沖擊阻力系數(shù)峰值及其寬度對比。從圖7可以看出，阻力系數(shù)峰值隨著速度增加呈現(xiàn)先減小后平穩(wěn)的趨勢。在100～300 m/s速度范圍內(nèi)阻力系數(shù)峰值隨著入水速度的增加而迅速降低； 而在500 m/s及以上速度入水，沖擊阻力系數(shù)峰值降低得不明顯。入水峰值寬度時間處于10-5s量級。峰值寬度隨著速度的增大逐漸減小。在1 000 m/s及以上速度，峰值寬度改變很小。

圖7　不同速度工況下阻力系數(shù)峰值和寬度對比圖Fig. 7　Comparison of peak values and peak widths of drag coefficient at different speeds

2.3頭部直徑對沖擊阻力系數(shù)的影響

計算選取速度為1 000 m/s工況，原有計算模型頭型及長度不變，頭部直徑分別擴大為原來的2倍和5倍，研究頭部直徑對阻力系數(shù)的影響。

Q 1:燃燒室輸出到水浴的熱量;Q 2:水浴傳遞給管程內(nèi)LNG的熱量;C:水的熱容;T 1:水浴溫度;T 2:管程LNG的溫度;R:水浴傳熱熱阻;Gp 1:主被控對象(NG出口溫度);Gp 2:副被控對象(水浴溫度);Gc1:主回路控制器;Gc2:副回路控制器;τ:副被控對象時滯;τ1:主控對象存在時滯;J(s):燃料氣的干擾;E(s):入口LNG的干擾;GE(s):入口LNG的傳遞函數(shù);Ge(s):前饋補償通道函數(shù);Y 1(s):主回路的輸出;Y 2(s):副回路的輸出。

圖8顯示的是不同頭部直徑以1 000 m/s速度入水的阻力系數(shù)變化。從圖中可以看出3種直徑的阻力系數(shù)依次增大，其值分別為3.38，4.42和5.39，峰值寬度分別為1.6×10-6s，2.9×10-6s和4.5×10-6s。這是因為隨著頭部直徑變大，側(cè)向排開迎流面水的困難增加，阻力系數(shù)峰值增加，入水以后阻力系數(shù)下降緩慢。

圖8　不同頭部直徑阻力系數(shù)對比Fig. 8　Comparison of different drag coefficients for different head diameters

2.4水可壓縮性對沖擊阻力系數(shù)的影響

從理論研究和生活經(jīng)驗都可以知道，在沖擊速度達到一定大小后，水不能按照低速下的不可壓縮流體對待，需要考慮水的壓縮性。

圖9顯示的是不同速度工況下的水可壓縮性對阻力系數(shù)的影響。從圖中可以看出，在速度為200 m/s時，考慮水可壓縮性情況下阻力系數(shù)峰值為5.8，不考慮水可壓縮性情況下為6.2，兩者相對誤差為6.9%。同時，由于接觸水面時對水有壓縮作用，考慮水壓縮性情況下入水沖擊峰值出現(xiàn)的時間發(fā)生后移； 在速度為500 m/s時，考慮水壓縮性情況下阻力系數(shù)峰值出現(xiàn)的時間同樣發(fā)生延后，其值為3.2，而在水不可壓縮情況下沖擊阻力系數(shù)為4.4，兩者的相對誤差達到了37.5%。可見在入水速度為500 m/s時，水的可壓縮性對入水沖擊阻力影響非常巨大，不能忽略不計。

圖9　不同速度工況下水可壓縮性對比Fig. 9　Comparison of water compressibility at different speeds

2.5空氣壓縮性對沖擊阻力系數(shù)的影響

當回轉(zhuǎn)體在空氣中以超聲速運動時，會壓縮空氣產(chǎn)生激波。回轉(zhuǎn)體接近水面時，壓縮空氣與水相互作用是影響沖擊力的一個重要因素。計算采用如下策略： 射彈在空氣中運動時空氣和水同時使用可壓縮模型； 當射彈與水面接觸阻力系數(shù)達到峰值，關(guān)閉空氣的可壓縮性，保留水的可壓縮性直到阻力系數(shù)趨于某一穩(wěn)定值，計算結(jié)束。

圖10顯示的是回轉(zhuǎn)體以1 000 m/s速度在空氣中飛行的壓力云圖。從圖中可以看出，回轉(zhuǎn)體在空氣中飛行，產(chǎn)生了脫體激波。當回轉(zhuǎn)體接近水面時，激波與水面相互作用使臨近的水相區(qū)域產(chǎn)生了壓力梯度。

圖11和圖12為速度1 000 m/s和2 000 m/s時，考慮與不考慮空氣壓縮性情況下入水阻力系數(shù)峰值的對比。從圖中可知，回轉(zhuǎn)體接近水面過程中，考慮空氣可壓縮情況下的阻力系數(shù)開始增大的時間比空氣不可壓縮要晚，阻力系數(shù)峰值出現(xiàn)的時間延后； 速度為1 000 m/s時，在回轉(zhuǎn)體與水接觸瞬間，考慮空氣可壓縮情況的阻力系數(shù)為2.55，不考慮空氣可壓縮情況的阻力系數(shù)為3.25，兩者相差27.5%。這是因為可壓縮空氣在回轉(zhuǎn)體接近水面過程中起到了“空氣墊”作用，使得阻力系數(shù)降低。速度為2 000 m/s時，在回轉(zhuǎn)體與水接觸瞬間，考慮和不考慮空氣可壓縮情況的阻力系數(shù)分別為2.2和2.7，兩者相差22.7%。

圖10　1 000 m/s速度空中飛行壓力云圖Fig. 10　Pressure contour in air for flight at a speed of 1000 m/s

圖11　1 000 m/s速度工況空氣可壓縮性對比Fig. 11　Comparison of air compressibility at a speed of 1 000 m/s

圖12　2 000 m/s速度工況空氣可壓縮性對比Fig. 12　Comparison of air compressibility at a speed of 2 000 m/s

3　結(jié)論

采用CFD方法模擬了平頭回轉(zhuǎn)體模型在亞音速、跨音速以及超音速狀態(tài)下的入水過程，分析了水可壓縮性和空氣可壓縮性對入水沖擊的影響，得到了如下結(jié)論。

1）平頭回轉(zhuǎn)體以1 000 m/s速度入水，在入水深度達到2倍體長時，入水空泡依然沒有閉合的趨勢，且在氣液兩相的交界區(qū)域附近產(chǎn)生了大量水蒸氣；

2）平頭回轉(zhuǎn)體入水速度在100～2 000 m/s范圍內(nèi)，入水沖擊阻力系數(shù)峰值先減小后平穩(wěn)。峰值寬度處于10-5s量級，其大小隨著速度的增大呈現(xiàn)逐漸減小趨勢； 在1 000 m/s速度工況下，增加頭部直徑，阻力系數(shù)峰值增加，峰值寬度增加；

3）考慮水壓縮性的影響，模型入水沖擊阻力系數(shù)峰值減小，峰值的時間延后。在速度為200 m/s時，考慮水可壓縮性和不考慮水可壓縮情況下，阻力系數(shù)相對誤差為6.9%； 在速度為500 m/s時，兩者的相對誤差達到了37.5%；

4）考慮空氣可壓縮性，模型接近水面時，激波與水面相互作用使臨近的水相區(qū)域產(chǎn)生壓力梯度。模型在空氣中運動的阻力系數(shù)增大的時間延后； 在接觸水面時刻的阻力系數(shù)減??； 阻力系數(shù)出現(xiàn)峰值的時刻延后。

［1］Logvinovich G V. Hydrodynamics of Free-boundary Flows［J］. Jerusalem，Israel Program for Scientific Translation，1972： 103-113.

［2］王永虎. 空投雷彈入水沖擊頭型特性參數(shù)分析［J］. 航空計算技術(shù)，2010，40（11）： 14-17.

Wang Yong-hu. Nose Performance Description Coefficience of Airborne Torpedo and Deep-mine during Water-entry Impact［J］. Aeronautical Computing Techbique，2010，40（11）： 14-17.

［3］陳宇翔，郜冶，劉乾坤. 應(yīng)用VOF方法的水平圓柱入水數(shù)值模擬［J］. 哈爾濱工程大學(xué)學(xué)報，2011，32（11）： 1439-1442.

Chen Yu-xiang，Gao Ye，Liu Qian-kun. Numerical Simulation of Water-entry in a Horizontal Circular Cylinder Using the Volume of Fluid（VOF）Method［J］. Journal of Harbin Engineering University，2011，32（11）： 1439-1442.

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

He Chun-tao，Wang Cong，Min Jing-xin，et al. Numerical Simulation of Early Air-Cabity of Cylinder Cone With Vertial Water-entry［J］. Engineering Mechanics，2012，29（4）： 237-243.

［5］邱海強，袁緒龍，王亞東，等. 回轉(zhuǎn)體高速垂直入水沖擊載荷和空泡形態(tài)仿真［J］. 魚雷技術(shù)，2013，21（3）： 161-164.

Qiu Hai-qiang，Yuan Xu-long，Wang Ya-dong，et al. Simulation on Impact Load and Cavity Shape in High Speed Vertical Water Entry for an Axisymmetric Body［J］. Torpedo Technology，2013，21（3）： 161-164.

［6］Schaffar M，Rey C，Boeglen G. Behavior of Supercavitating Projectiles Fired Horizontally in a Water Tank： Theory and Experiments［C］//The 35th AIAA Fluid Dynamics Conference and Exhibit，Toronto，Ontario Canada，2005.

（責任編輯： 陳曦）

Numerical Simulation on Vertical Water Entry Impact of Axisymmetric Body at Supersonic Speed

SUN Kai，DANG Jian-jun，HAO Wei-min，JIANG Bin
（School of Marine Science and Technology，Northwestern Polytechnical University，Xi′an 710072，China）

Water entry impact is one of the hot topics in the modern underwater weapons research. With the volume of fluid （VOF）model and dynamic mesh in commercial CFD software，the water entry processes of an axisymmetric body with flat top are simulated at subsonic，transonic and supersonic speeds，respectively. The drag coefficient，head diameter，water compressibility and the effect of air shock wave on the water entry processes of the simulated model are acquired at three different speeds. Simulation results show that： 1）with the increase of speed，the peak value of drag coefficient first decreases then remains stable in the end，and the width of peak value decreases gradually； 2）the peak value and its width of drag coefficient increases with the increase in head diameter； 3）the water compressibility delays the appearance of the peak value of drag coefficient and reduces the peak value； and 4）air shock wave reduces the value of drag coefficient at water surface.

axisymmetric body； water entry process； volume of fluid（VOF）model； drag coefficient

TJ630.1

A

1673-1948（2015）01-0002-05

2014-11-17；

2014-12-08.

孫凱（1989-），男，在讀碩士，主要研究方向為水下航行器設(shè)計.