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

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

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

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

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

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

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

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

0　引言

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

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

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

1　模型建立與驗(yàn)證

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

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

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

其中： p0是水的參考絕對(duì)壓力；0ρ是在參考?jí)毫ο滤膮⒖济芏龋?K0是在參考?jí)毫ο碌膮⒖俭w積彈性模量； n是密度指數(shù)； p是水受到的絕對(duì)壓力； ρ是在壓力為p時(shí)的密度。

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

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

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

6自由度求解器在計(jì)算平移運(yùn)動(dòng)時(shí)加速度計(jì)算較精確，因而文中采用了該求解器模型。在慣性參考坐標(biāo)系下建立剛體平移運(yùn)動(dòng)控制方程為

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

1.2計(jì)算模型簡(jiǎn)化及邊界條件

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

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

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

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

圖1　計(jì)算流場(chǎng)及邊界條件Fig. 1　Calculation of flow field and boundary conditions

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

1.3模型驗(yàn)證

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

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

圖2　試驗(yàn)與仿真結(jié)果流場(chǎng)對(duì)比圖Fig. 2　Comparison between experimental and simulated flow field

圖3　空泡外形對(duì)比圖Fig. 3　Comparison of cavitation outlines

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

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

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

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

圖4　1 000 m/s速度入水階段水相體積分?jǐn)?shù)云圖

圖5　1 000 m/s速度入水階段蒸汽相體積分?jǐn)?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)生了大量的水蒸氣。這是因?yàn)榛剞D(zhuǎn)體入水速度超過(guò)聲速，空氣不能及時(shí)補(bǔ)充進(jìn)空泡內(nèi)，導(dǎo)致空泡內(nèi)的壓力低于空化壓力，產(chǎn)生了空化現(xiàn)象。正是由于空化現(xiàn)象，入水空泡內(nèi)氣相和水相的交界面變得“凹凸不平”。

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

圖6顯示的是回轉(zhuǎn)體以1 000 m/s速度入水過(guò)程阻力系數(shù)在峰值附近的變化過(guò)程。在空中飛行階段，由于空氣的阻力很小，阻力系數(shù)很小； 回轉(zhuǎn)體接近水面，阻力系數(shù)迅速增大，在2×10-5s時(shí)回轉(zhuǎn)體與水面接觸，此時(shí)阻力系數(shù)為2.5； 回轉(zhuǎn)體繼續(xù)運(yùn)動(dòng)，阻力系數(shù)增大，在2.02×10-5s時(shí)阻力系數(shù)達(dá)到最大值3.23； 阻力系數(shù)在達(dá)到最大值后逐漸降低，整個(gè)沖擊峰值持續(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í)沖擊阻力系數(shù)峰值及其寬度對(duì)比。從圖7可以看出，阻力系數(shù)峰值隨著速度增加呈現(xiàn)先減小后平穩(wěn)的趨勢(shì)。在100～300 m/s速度范圍內(nèi)阻力系數(shù)峰值隨著入水速度的增加而迅速降低； 而在500 m/s及以上速度入水，沖擊阻力系數(shù)峰值降低得不明顯。入水峰值寬度時(shí)間處于10-5s量級(jí)。峰值寬度隨著速度的增大逐漸減小。在1 000 m/s及以上速度，峰值寬度改變很小。

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

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

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

Q 1:燃燒室輸出到水浴的熱量;Q 2:水浴傳遞給管程內(nèi)LNG的熱量;C:水的熱容;T 1:水浴溫度;T 2:管程LNG的溫度;R:水浴傳熱熱阻;Gp 1:主被控對(duì)象(NG出口溫度);Gp 2:副被控對(duì)象(水浴溫度);Gc1:主回路控制器;Gc2:副回路控制器;τ:副被控對(duì)象時(shí)滯;τ1:主控對(duì)象存在時(shí)滯;J(s):燃料氣的干擾;E(s):入口LNG的干擾;GE(s):入口LNG的傳遞函數(shù);Ge(s):前饋補(bǔ)償通道函數(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。這是因?yàn)殡S著頭部直徑變大，側(cè)向排開(kāi)迎流面水的困難增加，阻力系數(shù)峰值增加，入水以后阻力系數(shù)下降緩慢。

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

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

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

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

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

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

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

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

圖11和圖12為速度1 000 m/s和2 000 m/s時(shí)，考慮與不考慮空氣壓縮性情況下入水阻力系數(shù)峰值的對(duì)比。從圖中可知，回轉(zhuǎn)體接近水面過(guò)程中，考慮空氣可壓縮情況下的阻力系數(shù)開(kāi)始增大的時(shí)間比空氣不可壓縮要晚，阻力系數(shù)峰值出現(xiàn)的時(shí)間延后； 速度為1 000 m/s時(shí)，在回轉(zhuǎn)體與水接觸瞬間，考慮空氣可壓縮情況的阻力系數(shù)為2.55，不考慮空氣可壓縮情況的阻力系數(shù)為3.25，兩者相差27.5%。這是因?yàn)榭蓧嚎s空氣在回轉(zhuǎn)體接近水面過(guò)程中起到了“空氣墊”作用，使得阻力系數(shù)降低。速度為2 000 m/s時(shí)，在回轉(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速度工況空氣可壓縮性對(duì)比Fig. 11　Comparison of air compressibility at a speed of 1 000 m/s

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

3　結(jié)論

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

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

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

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

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

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

［2］王永虎. 空投雷彈入水沖擊頭型特性參數(shù)分析［J］. 航空計(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方法的水平圓柱入水?dāng)?shù)值模擬［J］. 哈爾濱工程大學(xué)學(xué)報(bào)，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］邱海強(qiáng)，袁緒龍，王亞?wèn)|，等. 回轉(zhuǎn)體高速垂直入水沖擊載荷和空泡形態(tài)仿真［J］. 魚(yú)雷技術(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.

（責(zé)任編輯： 陳曦）

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-），男，在讀碩士，主要研究方向?yàn)樗潞叫衅髟O(shè)計(jì).