郝亮，劉濤濤,2，黃彪,2，王國玉,2

（1.北京理工大學 機械與車輛學院，北京 100081；2.北京理工大學 重慶創(chuàng)新中心，重慶 401120）

航行體入水流動特性及彈道特性等問題的研究一直是國內外在發(fā)展相關裝備研制過程中的焦點。其入水過程涉及到介質突變、湍動、相變等大量復雜流動，具有強瞬時性和非定常運動特性，是物體從空中彈道過渡到水下彈道的一個重要環(huán)節(jié)。入水過程包含的非定常流體力學問題，涉及到激波、自由液面和超空泡相互作用，還有跨音速跨介質運動造成的結構動態(tài)響應力學問題，對入水過程彈道穩(wěn)定性和超空泡射彈結構具有很大的影響。入水過程是整個航行過程中流動現(xiàn)象和載荷最復雜、航行體姿態(tài)和彈道最不確定的非受控階段，尤其入水初期，介質突變、空泡演變過程和彈道特性的研究尤為重要。

Lundstrom對穿甲彈以跨音速800~1070 m/s入水過程開展了試驗研究，分析了入水空泡發(fā)展、潰滅的規(guī)律，并基于能量守恒推導出空泡預測半徑的公式。之后，Hrubes對超空泡射彈開展了以亞音速、跨音速以及超音速入水試驗，獲得了射彈在水中航行時空泡的形態(tài)，以及超音速航行時的水中激波現(xiàn)象。Owis等開發(fā)了一種基于Navier-Stokes方程求解可壓縮流動的數值計算方法，通過對繞半球體和圓錐體的空化流動進行數值模擬計算，驗證了數值方法的可行性。Gekle等對圓盤勻速垂直入水過程進行了試驗和數值模擬研究，針對伴隨空泡深閉合產生的跨音速射流形成過程建立了數學模型，數值模擬計算結果與試驗結果取得了較好的一致性。Jiang等采用有限體積法對水平放置的圓柱體入水過程進行了數值模擬，研究了湍流減阻添加劑在自然超空泡流過程中的作用，結果表明，當加入減阻添加劑時，空泡的長度和直徑大于在水中的尺寸。葉取源基于非線性自由面理論，采用了歐拉–拉格朗日混合邊界元方法對錐頭和圓盤平頭物體垂直入水過程進行了數值計算，分析了Froude數等因素對入水空泡發(fā)展的影響。李達欽等對繞跨音速射彈可壓縮超空泡流動進行了數值模擬，研究了馬赫數對超空泡流動特性的影響，推導了激波前后流動參數的表達式，同時分析了沖擊波上壓力與密度之間的關系。

早期航行體入水彈道研究主要針對二戰(zhàn)中魚雷失靈等問題，具有很強的軍事應用背景。1945年，Mason和Slichter開始了對空中發(fā)射導彈后入水的現(xiàn)象的研究。May等對不同材料的球體進行了垂直入水的試驗，分析了不同密度入水時的附加質量力。Waugh等在《水彈道學模擬》中總結了水彈道學入水試驗方法和研究現(xiàn)狀。當射彈以跨音速傾斜入水時，彈體更容易由于不穩(wěn)定性而失穩(wěn)?？缫羲購椀纼A斜入水的研究對于實現(xiàn)水下精確打擊具有重要研究意義。Truscott對特定射彈斜射入水過程進行了試驗研究，記錄了射彈以0.3~1.0馬赫數的入水速度，在入水角度為5°~15°系列工況下的入水空泡形態(tài)演變過程及彈道特征。當射彈跨音速入水時，流體的可壓縮性對入水沖擊載荷造成了很大影響。Eroshin等考慮了沖擊過程中流體的可壓縮性，研究表明，考慮可壓縮性的入水沖擊載荷偏低。

由于跨音速入水試驗難度大，得到的觀測數據有限，不少學者通過理論和數值計算方法研究了跨音速入水過程。隨著計算流體力學的發(fā)展，基于N-S方程求解的計算模型結合動態(tài)網格技術的數值模擬方法，可以考慮湍流特性對空穴形成和發(fā)展的影響，并可以實現(xiàn)非定常流動、多自由度水彈道的計算，數值計算逐漸成為研究入水問題的重要方法。何春濤以RAMICS為研究對象，對超空泡形態(tài)和阻力系數進行了數值模擬研究，分析了彈身母線形狀、彈身長細比以及空化器直徑對空泡形態(tài)和阻力的影響，并且對帶和不帶尾翼的超空泡射彈進行了三維數值模擬，分析了尾翼穿透現(xiàn)象及其對超空泡穩(wěn)定和阻力的影響。秦陽等對截錐頭彈體以500 m/s跨音速斜射入水過程進行了數值模擬研究，結果表明，入水角度越大，壓力峰值越高。從國內外的研究現(xiàn)狀來看，對于入水過程的研究主要集中在低速垂直入水，對于跨音速尤其是跨音速斜射入水的研究較少。

本文基于數值計算方法，對跨音速多自由度入水展開分析，主要研究了傾斜入水過程的可壓縮流體動力特性與彈道特性，討論了入水攻角對跨音速入水過程空泡形態(tài)、流場結構以及彈道特性的影響。

1 數值計算方法

1.1 基本控制方程

本文采用VOF均相流模型對航行體跨音速入水過程多相流動進行求解。在整個計算域內，VOF多相流模型通過對互不相溶的流體求解同一個動量方程組，并追蹤每種流體的體積分數來模擬多相流。VOF模型對引入模型的每一相，引入一個稱為單元相體積分數的變量，在每個控制容積中，所有相的體積分數之和為1，混合介質的連續(xù)方程和動量方程見式（1）、（2）。

式中：和分別為混合介質的黏度和密度；、（=1,2,3；=1,2,3）表示3個正交坐標方向。

1.2 湍流模型

本文采用SST-湍流模型來封閉RANS方程。Menter提出的SST-模型，在近壁面利用-模型，在遠壁面利用-模型求解，并且考慮到了湍流剪切應力的傳輸，得到了壁面分離流動問題的精確預測。湍流動能和湍流頻率輸運方程見式（3）、（4）。

式中：Г、Г為湍流擴散系數；G、G為湍流生成項；Y、Y為湍流動能耗散項；S、S為自定義項；G為交叉擴散項；為基于平均體積分數的密度；（= 1,2,3）表示3個正交坐標方向。

1.3 空化模型

Zwart空化模型結合了泡間兩相流動理論，認為空穴的內外為連續(xù)體，將空穴看作一種密度發(fā)生劇烈變化的可壓縮黏性流體，很大程度上考慮了空化生長以及潰滅時空泡體積變化的影響，其適用于模擬空化的非定常流動過程。在該模型中，單位體積內的相間傳輸速率為：

式中：為汽核體積分數；為汽泡半徑；為汽化壓強；和分別是蒸發(fā)和凝結經驗系數。在計算中，=5×10，=1×10m，=50，=0.01。

1.4 液相可壓縮性修正

水的彈性模量為2.2 GPa，一般情況下可以忽略其可壓縮性，但當彈體航行速度達到跨音速或者超音速階段時，需要考慮水的可壓縮性。介質的可壓縮性對射彈空泡形態(tài)和流體動力均有顯著影響，在強可壓縮環(huán)境下，基于常規(guī)計算無法捕捉密度和溫度隨壓力場動態(tài)變化的過程。為了解析上述復雜耦合問題，在N-S方程的計算框架內，必須加入能量方程，其本質為熱力學第一定律，見式（7）。

式中：為內能，=（為比熱容，為溫度）；為流體的傳熱系統(tǒng)；S為熱源項。

本文數值計算選擇Tait方程作為可壓縮流體的狀態(tài)方程，忽略溫度影響的可壓縮液相水的Tait方程見式（8）。

其中，和是可以通過假設體積模量是壓力的線性函數確定的系數，其值基于參考狀態(tài)壓力、密度和體積模量求得。

Tait方程的簡化形式可以寫成：

式中：為參考壓力，取101 325 Pa；為參考壓力下的密度，取998.2 kg/m；為參考壓力下的體積彈性模量，取2.2 GPa；為密度指數，取7.15；為液體壓強；分別為壓力下的密度和體積彈性模量。

1.5 邊界條件設置與網格劃分

本文計算采用的航行體模型如圖1所示，為航行體圓柱段長度，為航行體頭部長度，為航行體直徑。航行體總長度為240 mm，其中=20 mm，=220 mm，L=20 mm。

圖1 航行體模型Fig.1 Model of vehicle

對于航行體斜射入水過程，本文只考慮縱向平面的運動，即、方向的位移和方向的偏轉3個自由度。計算域三維示意圖及邊界條件設置如圖2所示。航行體頭部距水面的高度為1.25，計算域高度為100，長度為100，寬度為25。航行體表面和Wall設置為無滑移固壁，上下邊界設置為壓力出口，并且以彈體頭部剛觸及水面為0時刻。文中采用混合密度網格進行網格劃分，航行體附近網格采用非結構網格局部加密，全局采用六面體結構網格。全局網格和局部加密網格如圖3所示。為了減少網格，節(jié)省計算時間，航行體采用半模，航行體切面與流域前面設置為Symmetry，以減少網格數量。

圖2 計算模型Fig.2 Calculation model

圖3 網格劃分Fig.3 Schematic diagram of meshing

采用上文建立的數值計算方法進行模擬，選取入水角度為45°工況。計算采用基于壓力的求解器，采用壓力–速度耦合，利用SIMPLEC求解策略，利用二階迎風格式求解密度與動量項，利用一階迎風格式求解湍動能與湍動能耗散率，利用一階隱式方程進行瞬態(tài)求解。除此之外，數值計算方法基于有限體積法，網格質量和數量對計算方法影響很大，但網格數量過大會影響計算效率，因此在保證計算精度的情況下要盡量減少網格數量。對航行體周圍網格以及加密區(qū)域邊界節(jié)點進行加密，建立3種不同網格密度的網格，網格數量分別為160萬、400萬、600萬。3種網格下的阻力系數變化如圖4所示。由圖4可知，網格為160萬時，阻力系數峰值和穩(wěn)定后的阻力系數偏低，隨著網格加密，阻力系數變化趨于一致。為了節(jié)省計算時間，選用400萬網格密度進行后續(xù)計算分析。

圖4 不同網格數量下阻力系數變化曲線Fig.4 Variation curves of resistance coefficients under different grid numbers

1.6 數值計算方法驗證

為了驗證數值計算方法的準確性和可行性，對比了利用上述數值計算方法的計算結果與魏卓慧和Hurbes的試驗結果。入水角度=45°斜射入水過程中，航行體速度和彈道軌跡隨時間變化的試驗結果和數值結果對比如圖5所示，其中航行體運動軌跡以航行體質心位置表示。由圖5a可知，在入水過程中，隨著入水深度的增加，速度不斷降低，入水速度數值計算結果和試驗結果吻合較好。由圖5b可知，數值計算結果彈道軌跡與試驗結果基本吻合。綜上所述，數值計算結果與試驗結果取得了良好的一致性。

圖5 數值計算結果與試驗結果對比Fig.5 Comparison between numerical results and experimental results: a) speed; b) projectile displacement

為進一步驗證數值計算方法對跨音速水中激波現(xiàn)象的模擬，圖6給出了=1.03（1540 m/s）下數值計算結果與Hurbes的試驗圖像結果。從數值計算結果的壓力云圖可以看出，射彈頭部出現(xiàn)的高壓現(xiàn)象與試驗結果一致，形成了弓狀激波，并且激波面與試驗中捕捉到的激波面基本吻合，充分驗證了數值計算方法的準確性與可行性。

圖6 Ma=1.03時射彈超空化試驗圖像與數值計算結果對比Fig.6 Comparison between projectile supercavity test image(a) and numerical results (b) when Ma=1.03

2 結果及分析

2.1 液相可壓縮性對入水流動特性的影響

本文模擬計算了不同時刻下航行體參考線位置的壓力分布，參考線如圖1所示。當入水角度為45°，入水速度為1 500 m/s時，其壓力分布如圖7所示。由圖7可知，在0.025 ms時，航行體周圍壓力分布具有很強的不對稱性，迎流面?zhèn)群叫畜w附近壓力最大值明顯低于背流面。隨著入水深度的增加，這種不對稱性逐漸降低，航行體兩側壓力分布趨于一致。對于不可壓縮工況，航行體兩側壓力最大值隨入水時間的增加而增大，并且在初始時期迎流面?zhèn)葔毫υ鲩L明顯，背流面?zhèn)葔毫τ绊懛秶静蛔?，而迎流面壓力分布范圍增加。在可壓縮流場中，背流面?zhèn)葔毫ψ畲笾迪葴p小、后增大，而迎流面?zhèn)葔毫ψ畲笾笛杆僭龃蠛缶徛黾樱⑶覊毫Ψ植佳貜较蚍较蜓杆僭黾雍蟪霈F(xiàn)陡降現(xiàn)象。航行體兩側最大壓力差值較為明顯，并且隨著入水深度的增加，由于航行體速度衰減，激波角增大，激波強度減弱，壓力發(fā)生陡降的轉折點逐漸遠離航行體，陡降趨勢減弱。

圖7 航行體附近壓力分布Fig.7 Pressure distribution near the vehicle: a) incompressible condition; b) compressible condition

不考慮和考慮可壓縮性下，入水過程中空氣相與蒸汽相云圖見圖8，其中黑色實線為空泡輪廓。由于計算深度有限，僅能觀察到航行體帶空泡航行階段，入水空泡始終與大氣相連?？梢钥闯?，入水空泡由空氣相和蒸汽相構成。在入水初期，背流面?zhèn)瓤张輲缀跞烧羝鄻嫵桑諝庀嘀饕獜挠髅鎮(zhèn)裙嗳?。隨著入水深度的增加，蒸汽相不斷增加，并且分布區(qū)域緊鄰兩側空泡壁面，航行體兩側存在的空氣不斷減少。對比不考慮和考慮可壓縮性工況可以看出，兩者水汽分布不同主要體現(xiàn)在迎流面。可壓縮工況下，空泡內部的空氣相明顯多于不考慮可壓縮工況（如6.5 μs時刻），迎流面蒸汽相僅沿空泡壁面存在。不考慮可壓縮性時，蒸汽相存在范圍寬度更大。

圖8 不同入水時刻空氣相和水蒸氣相分布云圖Fig.8 Cloud diagram of air phase and water vapor phase distribution at different water entry times:a) incompressible condition; b) compressible condition

2.2 入水攻角對入水流動特性的影響

本小節(jié)主要討論入水攻角對跨音速入水過程流體動力特性和彈道特性的影響。入水速度為1 500 m/s，入水角度為45°，入水攻角分別為：–3°、–1°、0°、1°、3°。相同入水深度下，不同攻角的空泡輪廓以及對應時刻的水相分布云圖見圖9?？梢钥闯?，不同攻角下，空泡輪廓存在很強的相似性；在相同入水深度下，空泡直徑趨近一致，背流面空泡長度始終大于迎流面，迎流面產生的飛濺水冠高度始終高于背流面，水冠飛濺高度也不斷增長。不同攻角下，迎流面空泡直徑趨于一致，水冠飛濺高度和空泡半徑近乎一致，而背流面空泡輪廓存在較大差別。對于背流面空泡，在航行體頭部附近，空泡直徑幾乎無差別，但結合水相云圖可以看出，隨著正負攻角的增加，航行體錐段背流面或迎流面會逐漸出現(xiàn)沾濕，空泡完整性遭到破壞。在自由液面附近，空泡輪廓差異較大。隨著入水正攻角的增加，背流面空泡半徑增加；隨著負攻角增加，背流面空泡半徑減小。

圖9 不同攻角下的空泡輪廓Fig.9 Cavity profile at different angles of attack.

為了進一步研究攻角對航行體沾濕面積的影響，圖10給出了航行體無量綱沾濕面積隨入水時間的變化曲線?？梢钥闯觯礉衩娣e變化基本可以分為3個階段：1）在入水撞擊瞬間，航行體頭部以及航行體沾濕面積迅速增加，隨后航行體周圍開始形成入水空泡，沾濕面積迅速減??；2）由于在入水初期航行體始終與自由液面水相接觸，沾濕面積不斷增加；3）隨著入水深度繼續(xù)增加，航行體完全被空泡包裹，沾濕面積再次減小，并且趨于穩(wěn)定。在階段1）中，隨著入水攻角的增大，沾濕面積峰值增大，并且隨負攻角的增大而減小，隨正攻角的增加先減小、后增大。這是由于正負攻角的存在改變了航行體撞擊水面過程中背流面?zhèn)鹊恼礉衩娣e。在階段2）中，隨著負攻角的增大，航行體沾濕面積大幅增加，并且沾濕面積增長幅度增大；隨著正攻角的增加，航行體沾濕面積先減小、后增大。在階段3）中，–1°~1°攻角下，航行體完全被空泡包裹，航行體沾濕面積趨于一致，但是由于–3°和3°攻角航行體錐段存在沾濕，在穩(wěn)定階段沾濕面積高于其他工況。

圖10 不同攻角下航行體沾濕面積隨時間的變化曲線Fig.10 Variation curves of wetted area of vehicle with time under different angles of attack.

不同攻角下彈體俯仰角及其角速度的變化曲線如圖11所示。由圖11可知，在–1°~1°攻角內，彈體俯仰角及其角速度變化呈現(xiàn)相同的變化規(guī)律。航行體撞擊自由液面后，航行體俯仰角沿逆時針旋轉，并且隨著入水深度的增加而增大，角速度在入水撞擊后迅速增加。在帶空泡航行階段，角速度趨于定值，俯仰角基本呈線性增加。俯仰角偏轉值隨負攻角的增大而增加，隨正攻角增加，航行體俯仰角沿逆時針偏轉值減小。由于較大攻角工況下航行體與空泡壁面有接觸，造成這2個工況具有不同的變化規(guī)律。對于–3°攻角，在入水過程中，彈體錐段始終與背流面空泡接觸，航行體受力不穩(wěn)定性加劇，角速度隨入水深度的增加而增大，俯仰角變化不再呈線性增長。同理，3°攻角工況下，彈體錐段始終與迎流面空泡接觸，導致彈體偏轉一定角度后，開始反向偏轉，角速度在入水撞擊增至峰值后開始逐漸降低，最終反向增加。

圖11 不同攻角下彈體俯仰角及其角速度變化曲線Fig.11 Variation curves of pitching angles (a) and angular velocity (b) under different angles of attack

不同攻角下彈道軌跡和彈道偏移量隨時間的變化曲線如圖12所示?？梢钥闯?，不同攻角下，彈道軌跡接近直線。在入水初期，彈道均向上發(fā)生偏移，隨著入水深度的增加，彈道偏移量出現(xiàn)不同的變化規(guī)律。當入水攻角為–1°~1°時，彈道軌跡偏移量隨著入水深度逐漸增大，基本呈線性增長；在攻角為3°時，彈道軌跡偏移量隨入水深度的增加而急劇增加，彈道偏移量遠高于其他工況；當入水攻角為–3°時，彈道軌跡在入水初期由向右側偏移逐漸變換，甚至發(fā)生反向偏轉，彈道逐漸向左側偏移，且偏轉量增長速度逐漸加快。整體來看，隨著正攻角的增大，彈道軌跡正向偏轉量急劇增加；隨著負攻角增加，彈道偏移量先略有增加，隨后彈道偏移發(fā)生負向偏轉。

圖12 不同攻角彈道軌跡和彈道偏移量曲線Fig.12 Trajectory and trajectory offset curves at different angles of attack.

不同攻角下航行體水平方向和垂直方向受力系數的變化如圖13所示。由圖13可知，不同入水攻角下，垂直方向受力系數的峰值均高于水平方向。隨著正攻角的增加，垂直和水平方向的受力峰值減小；隨著負攻角增加，垂直和水平方向的受力峰值先增大、后減小。這是因為隨著正攻角的增加，在入水撞擊過程中，航行體頭部與水面接觸面積減小，受力系數峰值降低；隨負攻角的增加，雖然在入水沖擊過程航行體頭部與水面接觸面積增加，但速度與航行體軸線夾角不斷變大，法向速度減小，可能造成航行體受力降低。在帶空泡航行階段，隨著入水深度的增加，航行體逐漸被入水空泡包裹，航行體受力系數迅速降低，并趨于平穩(wěn)。在小幅度攻角變化范圍（–1°~1°）內，穩(wěn)定階段受力系數相差很小，但隨著攻角進一步增加，水平方向的受力系數隨負攻角的增大而增加，隨正攻角增加而減小，垂直方向受力系數變化與之相反。因為在–1°和1°攻角下，航行體在帶空泡穩(wěn)定航行過程中始終處于空泡包裹，但當入水攻角進一步增加，在航行體錐段將會出現(xiàn)沾濕，造成了垂直和水平方向受力變化。以–3°攻角為例，航行體錐段在帶空泡航行階段與迎流面始終接觸，航行體受到空泡壁面的反作用力，既增加了航行體水平方向受力，也減小了垂直方向受力。

圖13 不同入水攻角下水平方向和垂直方向受力系數變化Fig.13 Variation of force coefficients in horizontal direction(a) and vertical direction (b) under different angles of attack

3 結論

本文應用能量方程和水的Tait方程進行了液體可壓縮修正，基于SST湍流模型、Zwart空化模型和動網格技術，建立了航行體跨音速斜射入水數值模擬計算方法。數值計算結果與試驗結果對比取得了較好的一致性，驗證了數值計算方法的可行性。基于該數值計算方法，研究了航行體跨音速斜射入水過程的流動特性及其彈道特征，得出如下主要結論：

1）航行體跨音速入水過程中，流體的可壓縮性對航行體受力特性和流場結構有較大影響。在考慮流體可壓縮性時，航行體兩側壓力分布呈現(xiàn)出明顯的不對稱性，隨著深度的增加，不對稱性減弱，并且壓力峰值先增加、后減小。除此之外，空泡內部空氣明顯多于不考慮可壓縮工況。

2）航行體以不同攻角入水時，流場結構的發(fā)展具有相似性，但隨著攻角減小，入水空泡直徑略有減小。攻角為3°時，航行體沾濕面積最大，隨后逐漸與–3°攻角情況相同。

3）航行體以–3°/3°攻角入水時，航行體姿態(tài)發(fā)生明顯偏轉，并隨著入水深度的增加，彈道偏移顯著變化。