秦 楊,易文俊,管 軍

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

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

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

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

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

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

一般默認對于初始速度不是太快的入水情況，流體可壓縮性[11]可參與不考慮。本研究將流體介質(zhì)視作不可壓縮，同時不計流體黏性所導(dǎo)致的熱傳遞。VOF多相流模型中分別用αl，αg，αv表示液體、氣體和水蒸氣的體積分數(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空化模型對流動中的空化問題進行求解，在這個模型里水蒸汽相體積分數(shù)的輸運方程為

(6)

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

2 數(shù)值計算

2.1 計算模型

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

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

2.2 計算域及邊界條件

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

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

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

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

hi>(1+αs)hideal

(7)

hi<αchideal

(8)

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

(9)

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

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

圖3 流場網(wǎng)格

2.4 數(shù)值方法

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

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

3 數(shù)值計算結(jié)果驗證

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

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

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

圖4 空泡形態(tài)

圖5 位移和速度曲線

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖10 速度衰減曲線

圖11 位移變化曲線

圖12 彈體壓力變化曲線

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

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

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

5 結(jié)論

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

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

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