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

秦 楊,易文俊,管 軍

(南京理工大學 瞬態物理國家重點實驗室， 南京 210094)

航行體高速入水時，在流體力的作用下,在航行體表面一部分形成低壓區,當航行體最小壓力點處的壓力降低到飽和蒸汽壓時,液體介質會發生汽化而產生空泡[1]。空化流動的流體界面上會出現較大的密度比和流場參數梯度，這給數值求解這類問題帶來了很大的困難[2]。入水過程伴隨著大量復雜的流動現象如湍動、相變、可壓縮等，具有非定常、強瞬時及高載荷等特性，會對航行體的運動、結構性能產生嚴重影響[3]。

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

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

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

1 數學模型

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

αl+αg+αv=1

(1)

混合物的連續性方程為

(2)

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

混合相的動量守恒方程為

(3)

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

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

(4)

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

(5)

本研究采用Schnerr and Sauer空化模型對流動中的空化問題進行求解，在這個模型里水蒸汽相體積分數的輸運方程為

(6)

其中：Fvap=50和Fcond=0.001為經驗常數，αnuc=5×10-4為不可凝結氣體積分數;RB=1×10-6m為瑞利方程中的氣核半徑。對于流動中的湍流現象，本研究采用k-ωSST湍流模型對流體控制方程進行封閉求解。該模型能夠更恰當的描述湍流剪切應力的傳輸，在預測近壁區繞流和旋流有優勢。

2 數值計算

2.1 計算模型

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

圖1 射彈尺寸參數

2.2 計算域及邊界條件

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

圖2 計算域對稱面示意圖

2.3 動網格技術與網格劃分

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

hi>(1+αs)hideal

(7)

hi<αchideal

(8)

為了在實際計算中消去網格的對流效應，引入動網格后控制體V對變量φ的守恒方程如下：

(9)

式中：ρm是流體混合物密度；u是流體速度矢量；ug為動網格的運動速度；Г為擴散系數；Sφ為標量φ的源項；?V表示控制體積V的邊界。

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

圖3 流場網格

2.4 數值方法

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

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

3 數值計算結果驗證

湍流一直是流體力學實驗的難點，想要通過實驗獲得流場結構非常困難。隨著現在計算流體力學的蓬勃發展以及計算機性能的顯著提升，使用數值方法計算入水空泡流場變成了一種可行途徑。

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

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

圖4 空泡形態

圖5 位移和速度曲線

4 不同角度入水的計算結果分析

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

4.1 高速傾斜入水空泡形態及流場分析

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

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

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

圖6 45°入水過程0～1 ms空泡形態變化

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

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

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

圖8 0.5 ms時不同入水角度的空泡云圖

圖9 0.5 ms時不同入水角度的空泡形態

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

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

圖10 速度衰減曲線

圖11 位移變化曲線

圖12 彈體壓力變化曲線

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

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

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

5 結論

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

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

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