超空泡射彈高速傾斜入水的空化流動數(shù)值模擬

秦 楊,易文俊,管 軍

(南京理工大學(xué) 瞬態(tài)物理國家重點(diǎn)實(shí)驗(yàn)室， 南京 210094)

航行體高速入水時(shí)，在流體力的作用下,在航行體表面一部分形成低壓區(qū),當(dāng)航行體最小壓力點(diǎn)處的壓力降低到飽和蒸汽壓時(shí),液體介質(zhì)會發(fā)生汽化而產(chǎn)生空泡[1]?？栈鲃拥牧黧w界面上會出現(xiàn)較大的密度比和流場參數(shù)梯度，這給數(shù)值求解這類問題帶來了很大的困難[2]。入水過程伴隨著大量復(fù)雜的流動現(xiàn)象如湍動、相變、可壓縮等，具有非定常、強(qiáng)瞬時(shí)及高載荷等特性，會對航行體的運(yùn)動、結(jié)構(gòu)性能產(chǎn)生嚴(yán)重影響[3]。

入水角度是影響航行體入水空化流場特性的一個(gè)重要因素，過小的入水角度會造成航行體忽撲或者冒出水面，不能按照近似于空中軌跡延長線在水下運(yùn)動。過大的入水角度有可能使航行體不能最快到達(dá)目的點(diǎn)。因此選擇合適的入水角度是保證航行體優(yōu)良性能的重要前提。尤其近年來隨著宇宙飛船和火箭發(fā)動機(jī)等在水面回收等技術(shù)的應(yīng)用以及水下武器的發(fā)展，掌握入水角度對航行體入水過程流場特性的影響規(guī)律顯得更加迫切和重要。

20世紀(jì)60年代Logvinovich基于大量實(shí)驗(yàn)和理論研究，提出了經(jīng)典的基于勢流理論的空泡截面獨(dú)立膨脹原理，描述了運(yùn)動體入水空泡生成和演化過程，為分析空泡壁面運(yùn)動貢獻(xiàn)了重要的理論成果。Savchenko等[4]對不同形狀航行體空泡形態(tài)、空化器阻力特性等開展了研究。Sun等[5]采用具有完全非線性邊界條件的不可壓縮速度勢理論，對多種不同斜角錐體的傾斜入水進(jìn)行了仿真。Yves-Marie S[6]基于Wagner theory發(fā)展了數(shù)值求解三維物體入水的理論方法，通過與實(shí)驗(yàn)結(jié)果對比表明該方法準(zhǔn)確可靠。熊天紅等[7]對小攻角對水下高速射彈空泡形態(tài)的影響的研究表明攻角會影響空泡的對稱性，從而使得航行體失去平衡。宋武超等[8]研究了不同頭型回轉(zhuǎn)體入水過程中空泡形態(tài)的規(guī)律、運(yùn)動特性及流體動力特性。胡青青等[9-10]對不同頭型鈍體的超空泡流動特性在不同速度條件下進(jìn)行了實(shí)驗(yàn)以及數(shù)值模擬。

根據(jù)已有文獻(xiàn)，對入水問題空泡變化及流場特性的數(shù)值仿真研究，主要以垂直入水為主，然而實(shí)際中入水問題幾乎都是傾斜入水。本研究采用CFD軟件Ansys fluent18.2模擬二維情況下初速度為500 m/s的射彈以30°、45°和60°三個(gè)角度的傾斜入水過程，得到了不同角度下彈體入水空泡形態(tài)發(fā)展規(guī)律、彈道特性及流體動力特性變化規(guī)律，研究結(jié)果可為工程實(shí)踐提供理論參考。

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

一般默認(rèn)對于初始速度不是太快的入水情況，流體可壓縮性[11]可參與不考慮。本研究將流體介質(zhì)視作不可壓縮，同時(shí)不計(jì)流體黏性所導(dǎo)致的熱傳遞。VOF多相流模型中分別用αl，αg，αv表示液體、氣體和水蒸氣的體積分?jǐn)?shù)，它們的關(guān)系為:

αl+αg+αv=1

(1)

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

(2)

其中：ρm為混合物密度；ρl為水的密度;ρv為水蒸汽的密度;Vm為混合物速度矢量。

混合相的動量守恒方程為

(3)

其中：ui和uj為速度分量;μm為混合介質(zhì)的動力黏度;μt=ρmCuk2/ε為湍流動力黏度。其中Cu=0.09是常數(shù)，k為湍動能，ε是湍流耗散率。

ρm=αlρl+αgρg+αvρv

(4)

μm=αlμl+αgμg+αvμv

(5)

本研究采用Schnerr and Sauer空化模型對流動中的空化問題進(jìn)行求解，在這個(gè)模型里水蒸汽相體積分?jǐn)?shù)的輸運(yùn)方程為

(6)

其中：Fvap=50和Fcond=0.001為經(jīng)驗(yàn)常數(shù)，αnuc=5×10-4為不可凝結(jié)氣體積分?jǐn)?shù);RB=1×10-6m為瑞利方程中的氣核半徑。對于流動中的湍流現(xiàn)象，本研究采用k-ωSST湍流模型對流體控制方程進(jìn)行封閉求解。該模型能夠更恰當(dāng)?shù)拿枋鐾牧骷羟袘?yīng)力的傳輸，在預(yù)測近壁區(qū)繞流和旋流有優(yōu)勢。

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

2.1 計(jì)算模型

本文的數(shù)值模擬在二維情況下進(jìn)行。針對截錐形頭部的彈體入水問題，開展了入水角為30°、45°、60°3種工況下的高速傾斜入水?dāng)?shù)值模擬研究。射彈模型采用圓盤空化器，頭部為截錐體，后體部分為圓柱體。其中，彈體材料為普通鋼，密度為ρ=5 g/cm3。超空泡射彈的頭部是唯一穩(wěn)定的沾濕區(qū)域，還是90%以上航行阻力的來源。根據(jù)文獻(xiàn)[12]的實(shí)驗(yàn)數(shù)據(jù)結(jié)果，彈體空化器直徑與后體直徑的比值應(yīng)大于0.26。無附體超空泡射彈的尺寸參數(shù)如圖1所示。

圖1 射彈尺寸參數(shù)

2.2 計(jì)算域及邊界條件

計(jì)算域?yàn)閳A柱體，圖2為其對稱面示意圖。根據(jù)文獻(xiàn)[13]的數(shù)值模擬結(jié)論，流域的徑向尺寸應(yīng)大于46倍彈體最大直徑，選取計(jì)算域直徑為2 000 mm，高度為2 500 mm。氣水交界面取在坐標(biāo)原點(diǎn)下方25 mm處，空氣域高為500 mm，水域高為2 000 mm。初始狀態(tài)，彈體的軸線與x軸夾角為α，坐標(biāo)原點(diǎn)取在頭部中點(diǎn)處。外流域邊界條件均采用壓力出口邊界(pressure outlet)，采用用戶自定義函數(shù)(udf)對邊界上的壓力進(jìn)行定義；計(jì)算環(huán)境壓力P0=101 325 Pa。

圖2 計(jì)算域?qū)ΨQ面示意圖

2.3 動網(wǎng)格技術(shù)與網(wǎng)格劃分

在彈體入水的數(shù)值模擬中引入動網(wǎng)格技術(shù)，從而實(shí)現(xiàn)彈體高速傾斜入水運(yùn)動。由于采用三角形網(wǎng)格，所以選用彈性光順模型以及局部重劃模型更新網(wǎng)格。通過設(shè)置模型表面網(wǎng)格的幾何尺寸，并設(shè)置尺寸變化范圍，網(wǎng)格被壓縮或被拉升超出設(shè)定網(wǎng)格尺寸時(shí)，就會被合并或分裂出新網(wǎng)格層。首先設(shè)移動邊界附近的網(wǎng)格理想高度為hideal，網(wǎng)格分裂因子為αs，網(wǎng)格層潰滅因子為αc。當(dāng)網(wǎng)格被拉伸，且被拉伸的網(wǎng)格高度hi滿足式(8)[14]時(shí)，網(wǎng)格將根據(jù)指定的網(wǎng)格層高度分割網(wǎng)格；在網(wǎng)格被壓縮，且網(wǎng)格高度滿足式(9)時(shí)，被壓縮層網(wǎng)格將與其相鄰網(wǎng)格合并為新一層網(wǎng)格。

hi>(1+αs)hideal

(7)

hi<αchideal

(8)

為了在實(shí)際計(jì)算中消去網(wǎng)格的對流效應(yīng)，引入動網(wǎng)格后控制體V對變量φ的守恒方程如下：

(9)

式中：ρm是流體混合物密度；u是流體速度矢量；ug為動網(wǎng)格的運(yùn)動速度；Г為擴(kuò)散系數(shù)；Sφ為標(biāo)量φ的源項(xiàng)；?V表示控制體積V的邊界。

由于非定常問題的動網(wǎng)格計(jì)算常常需要計(jì)算巨量的數(shù)據(jù)，對計(jì)算機(jī)資源要求較高。因此，為了提高計(jì)算速度、改善計(jì)算精度，將計(jì)算域劃分為彈體運(yùn)動路徑上的加密區(qū)域和受彈體運(yùn)動影響較小的稀疏區(qū)域。為保證空氣域、汽水交界面以及空泡區(qū)域的計(jì)算結(jié)果精度，加密區(qū)域采用較密集的網(wǎng)格，且對彈體周圍網(wǎng)格進(jìn)一步加密。兩個(gè)區(qū)域的網(wǎng)格通過一組網(wǎng)格界面mesh interface滑移，在fluent中，通過mesh interface工具將兩個(gè)區(qū)域的網(wǎng)格界面對接起來。網(wǎng)格劃分結(jié)果如圖3所示。

圖3 流場網(wǎng)格

2.4 數(shù)值方法

在入水過程中，射彈的運(yùn)動軌跡是由彈體自身慣性以及彈體上的作用力(重力、水動力、氣動力等)和力矩共同決定的，因此彈體的運(yùn)動軌跡是與流場的計(jì)算相互耦合的。本文采用Ansys Fluent18.2提供的6DOF求解器，計(jì)算傾斜入水過程中彈體表面上的作用力和力矩，然后根據(jù)力的平衡(重力、水動力、氣動力等)，計(jì)算出平移加速度，再積分得到平移速度；基于力矩，計(jì)算出角加速度，再積分得到角速度，然后計(jì)算得出新的重心位置和歐拉角，最終解算出彈體傾斜入水的運(yùn)動軌跡。

采用基于VOF多相流模型的有限體積法對流體控制方程進(jìn)行時(shí)間和空間上的離散，在瞬態(tài)計(jì)算過程中速度與壓力的耦合計(jì)算采用PISO(Pressure Implicit with Splitting of Operators)算法；時(shí)間離散為一階精度，對流項(xiàng)采用Quick離散格式；壓力插值采用PRESTO！離散格式；綜合考慮計(jì)算時(shí)間與收斂性，采用一階迎風(fēng)格式離散動量方程。耗散項(xiàng)和湍流采用了二階迎風(fēng)格式；各相體積率離散采用CICSAM格式?；贑語言程序，采用UDF自編程定義入水彈體質(zhì)量、慣性矩及計(jì)算域邊界壓力，最終實(shí)現(xiàn)彈體的入水運(yùn)動。

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

湍流一直是流體力學(xué)實(shí)驗(yàn)的難點(diǎn)，想要通過實(shí)驗(yàn)獲得流場結(jié)構(gòu)非常困難。隨著現(xiàn)在計(jì)算流體力學(xué)的蓬勃發(fā)展以及計(jì)算機(jī)性能的顯著提升，使用數(shù)值方法計(jì)算入水空泡流場變成了一種可行途徑。

基于上文提出的數(shù)值模擬方法，采用不可壓縮液體作為介質(zhì)，將數(shù)值模擬結(jié)果與文獻(xiàn)[15]的試驗(yàn)結(jié)果進(jìn)行對比，從而檢驗(yàn)數(shù)值模型的可靠程度。根據(jù)該文獻(xiàn)的實(shí)驗(yàn)結(jié)果，獲得射彈以初速為440 m/s、入水角度為10.7°的入水瞬間試驗(yàn)照片，根據(jù)文獻(xiàn)中的射彈外形和運(yùn)動參數(shù)，進(jìn)行數(shù)值模擬，對比空泡輪廓、位移曲線和速度變化規(guī)律的數(shù)值計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果。

從圖4可以看出，入水0.1～1.0 ms，除了噴濺形態(tài)不太一致，入水的空泡云圖與試驗(yàn)照片具有較好的一致性。圖5給出了高速射彈傾斜入水位移曲線和速度變化規(guī)律的數(shù)值計(jì)算結(jié)果，lx、vx分別表示高速射彈水平方向位移與速度，ly、vy分別表示高速射彈豎直方向位移與速度。從圖5可以看出，仿真結(jié)果和文獻(xiàn)[15]的實(shí)驗(yàn)數(shù)據(jù)的位移和速度變化規(guī)律吻合度也很好，這驗(yàn)證研究中所采用的數(shù)值模擬方法對傾斜入水過程的計(jì)算結(jié)果是正確有效的。觀察圖4與圖5不難發(fā)現(xiàn)，數(shù)值計(jì)算結(jié)果能夠較好地反映入水過程中空泡形態(tài)及彈體彈道特性的變化過程。

圖4 空泡形態(tài)

圖5 位移和速度曲線

4 不同角度入水的計(jì)算結(jié)果分析

彈體浸水時(shí)的低壓效應(yīng)和浸水阻力不僅與頭型和長徑比有關(guān)，還與入水角大小有關(guān)系，本文使用數(shù)值計(jì)算的方法探討不同角度傾斜入水過程彈體的空泡形態(tài)、彈道特性及流體動力特性變化規(guī)律。

4.1 高速傾斜入水空泡形態(tài)及流場分析

對截錐頭彈體500 m/s 初始速度，入水角度分別為30°、45°、60°的傾斜入水問題開展數(shù)值模擬計(jì)算，分析不同入水角度對其空泡形態(tài)及流體動力特性等影響規(guī)律。

圖6給出了彈體45°入水過程0～1 ms空泡變化過程。彈體傾斜入水時(shí)，因?yàn)轭^部下側(cè)首先跟水接觸，在彈體前方引起濺水，來不及逃出的空氣在隆起的水堆和彈體間形成一個(gè)低壓空泡區(qū)。同時(shí)，水動壓力垂直于彈體的沾濕面，由于沾水面不是圓心在彈體重心的球面的一部分，因此作用于彈體沾濕面的合力將產(chǎn)生一個(gè)繞彈體重心的力矩，使得文中具有圓盤空化器的彈體頭部產(chǎn)生向下的偏離。這兩個(gè)因素引起的不平衡力矩，造成彈體在這期間發(fā)生一個(gè)角速度的階躍，即忽撲現(xiàn)象。在頭部完全沾水后，由于彈體繼續(xù)前進(jìn)時(shí)攻角為零，水動壓力產(chǎn)生的俯仰力矩變?yōu)榱?，因此彈體逐漸恢復(fù)之前的入水角度。

入水瞬間空腔和大氣相通，所以稱之為開空泡。從圖6中可以看出，彈體入水形成開空泡，由于水界面空氣卷入空泡中，空泡得以持續(xù)發(fā)展。并且由于入水速度快，彈體表面空泡不易閉合。

圖6 45°入水過程0～1 ms空泡形態(tài)變化

圖7是彈體不同角度入水0.5 ms時(shí)刻頭部附近的壓力等值線分布圖，從圖7中可以看出不同入水角度的壓力分布趨勢是基本一致的，最大壓力均集中在彈體頭部。

圖7 不同角度入水0.5 ms時(shí)刻彈體頭部的壓力等值線圖

圖8和圖9給出了0.5 ms時(shí)不同入水角度的空泡形態(tài)對比，可以看出，由于傾斜入水過程中彈體對水面的不對稱撞擊，其產(chǎn)生的噴濺很大部分出現(xiàn)在彈體的水平速度分量正方向。并且空泡右邊液面抬升隨著入水角度的增加而減小，而空泡左邊的液面抬升則增大。這是因?yàn)槿胨嵌容^大的時(shí)候，水平正方向的速度分量較小，在水平前方傳遞給液體的動能較小，因此右邊液面抬升較小；而豎直方向速度分量較大，傳遞給液體的動能較大，導(dǎo)致了左邊液面噴濺較明顯，液面抬升較大。

圖8 0.5 ms時(shí)不同入水角度的空泡云圖

圖9 0.5 ms時(shí)不同入水角度的空泡形態(tài)

4.2 傾斜入水彈道特性及流體動力特性分析

為進(jìn)一步探究初速度500 m/s彈體以30°、45°和60°3個(gè)角度傾斜入水的彈道特性及流體動力特性，對不同入射角度的高速彈體入水進(jìn)行了數(shù)值仿真，仿真結(jié)果如圖10、圖11和圖12所示。

圖10 速度衰減曲線

圖11 位移變化曲線

圖12 彈體壓力變化曲線

圖10為不同入水角度射彈總的速度以及重心水平方向、豎直方向速度曲線，觀察圖10(a)可知，入水角度為60°時(shí)，2 ms時(shí)間內(nèi)的速度衰減大于其他兩種情況，這是由于入水角度較大時(shí)入水瞬間彈體頭部的沾濕面積較大，受到了較大的阻力，動能損失較大。圖10(b)中高速射彈由于初始運(yùn)動方向的不同，30°入水角彈體的水平速度衰減速度較60°入水角彈體的衰減速度要快，豎直方向情況則正好相反。總體看來入水角度對總速度的影響不大，這主要是因?yàn)閺楏w在x、y方向受到阻力的合力大致相等。速度衰減曲線的斜率是逐漸減小的，表明彈體受到的阻力在入水瞬間達(dá)到最大值，然后逐漸減小。

圖11為不同角度入水射彈的水平方向位移lx和豎直方向位移ly曲線。從圖11中可以看出，入水角度為30°時(shí)，2 ms時(shí)間內(nèi)彈體在水平方向達(dá)到最大位移僅僅約為65D，入水角度為60°時(shí)，豎直方向有最大位移約60D。這說明彈體在入水瞬間受到了巨大的阻力。同一時(shí)刻，入水角度越小，則水平方向位移越大，入水深度越小。

圖12為不同角度入水射彈壓力變化曲線，由圖12可見，在入水之前，彈體受到的壓力近似為零，保持恒定；彈體撞擊自由液面后在極短時(shí)間內(nèi)其壓力出現(xiàn)峰值，此時(shí)彈體表面流場壓力最高可達(dá)大氣壓的千倍量級。入水角度越大，彈體表面壓力峰值越大，且壓力衰減速度越快。這是由于入水角度較大時(shí)入水瞬間彈體頭部的沾濕面積較大，持續(xù)受到了較大的阻力。在彈體觸水后，隨著入水深度的增大，壓力峰值逐漸下降，壓力數(shù)值緩慢減小，逐漸趨于穩(wěn)定，不同入水角度的彈體壓力值差距越來越小，但在2 ms內(nèi)依然保持較高的水平。

5 結(jié)論

1) 本研究采用的數(shù)值計(jì)算方法能夠有效模擬射彈入水過程中空泡形態(tài)的變化過程。彈體以不同角度入水產(chǎn)生的空泡形態(tài)差異較大，隨著入水角度的增加，空泡右邊液面抬升減小，而空泡左邊的液面抬升則增大。

2) 不同入水角度射彈的速度衰減曲線不太一致，隨著入水角度的增大，總速度的衰減率呈現(xiàn)微小增加的趨勢?？傮w來說，入水角度對射彈的總速度變化影響不大。

3) 入水初期，彈體受到較高沖擊載荷的作用，其壓力峰值可達(dá)數(shù)千倍大氣壓。入水角度越大，壓力峰值越大，入水瞬間壓力衰減速度越快。