高速航行體齊射出水過程的空化與運動特性研究

周東輝，賈會霞，施紅輝，王焯鍇

(1.浙江工業職業技術學院，紹興 312000；2.西北工業大學 寧波研究院，寧波 315100；3.浙江理工大學 機械與自動控制學院，杭州 310018)

0 引言

從水下發射超空泡射彈攻擊水面目標或者低空飛行目標，可提高武器發射的隱蔽性，增強其作戰效果。例如挪威 DSG 公司研制的多環境槍彈[1]，可實現水下5 m 深度發射，打擊1 000 m 飛行高度的直升機。德國研制的反氣墊船超空泡射彈，可預先在水下將射彈垂直布置，當目標出現后在不到1 s 時間內完成攻擊任務[2]。在實際應用中，為了提高射彈的命中概率和毀傷效果，往往需要采用齊射方式，即以一定發射空間間距、在極短時間間隔內連續發射多發射彈。齊射出水過程中，不僅會出現空化、湍動、穿越自由面等復雜的三相耦合流動現象，還存在多彈體超空泡流場的相互擾動，導致射彈的超空泡流動和運動特性非常復雜。因此，開展相關的研究具有重要的學術研究和工程應用價值。

針對帶空泡航行體的出水問題，其中出水空泡形成、發展、潰滅以及影響因素分析等方面一直是研究重點，這主要是由于空泡的形態、潰滅關系到航行體的受力及沖擊載荷，進而影響航行體出水過程的彈道穩定性和結構安全。在國外，Waugh 等[3]通過實驗研究了發射角度和空泡對導彈出水姿態的擾動影響，給出了出水空泡的形態。Nguyen 等[4]利用數值模擬方法研究了射彈勻速出水過程的超空泡流動，給出了射彈出水過程的阻力變化和超空泡形態。Logvinovich等[5]提出了空泡截面獨立膨脹原理，給出了計算超空泡特征直徑和特征長度的半經驗公式。Savchenko等[6]修正了Logvinovich 提出的計算超空泡輪廓的半經驗公式。在國內，褚學森等[7]利用均質多相流方法對圓柱體垂直出水進行了數值模擬，研究發現空泡在圓柱體的肩部和尾部形成，在穿越自由面的過程中潰滅。顏開等[8]對空化器朝水面高速運動產生的非定常空泡進行了理論研究，計算了空化器出水后空泡從脫落、收縮到潰滅的時間。魏海鵬等[9]分析了空化數及非凝結性氣體含量對潛射導彈空化流動的影響。施紅輝等[10-12]對射彈高速出水的超空泡流動進行了實驗和數值模擬研究，獲得了出水超空泡的演化過程，分析了超空泡發展對射彈出水姿態的影響以及射彈頭部形狀對射彈運動、超空泡形態等的影響。王一偉等[13-14]采用實驗和數值模擬相結合的方法研究了射彈垂直出水過程的非定常空化流動，給出了出水空泡的演化和射彈表面壓力分布特征，揭示了空泡發生潰滅的機制。付國強等[15]對細長體垂直出水過程的空化流動進行了實驗研究，分析了水深和發射速度對空泡發展的影響。陳瑛等[16]采用大渦模擬的方法對潛射航行體出水過程的空化流動進行了數值模擬，給出了出水空泡的演化過程、湍流結構和壓力特征，討論了潛射深度、發射速度、迎角的影響。Chen 等[17]利用自行設計的水下發射平臺開展了射彈垂直出水實驗，研究了射彈頭部形狀和發射速度對出水空泡演化和射彈運動的影響，分析了空化數、再進入射流與空泡穩定性之間的關系。

針對多個航行體出水的研究，Mnasri 等[18]建立雙圓柱體低速并聯出水的二維數值模擬，分析了流場的相互擾動。盧佳興等[19]開展了發射速度為15 m/s的兩發回轉體齊射出水實驗，獲得了齊射過程中回轉體的空泡演化特性和運動特性。施紅輝等[20]進行了兩發超空泡射彈連續出水過程的數值模擬，分析了出水過程超空泡流場的干擾規律以及超空泡演化過程對射彈運動的影響。畢鳳陽等[21]建立了多細長體水下齊射多相流動與多體運動耦合數值模擬的計算模型，研究了細長體筒中段及離筒初期的流動、運動及干擾特性。Xu 等[22]進行了兩發射彈齊射出筒的六自由度數值模擬，研究了尾流的演化及其與射彈的相互作用。Gao 等[23]基于VOF 方法研究了射彈水下齊射過程中的流動干擾特性，獲得了射彈表面的壓力分布、水動力特性、彈道特性以及俯仰角的變化，由于射彈的發射速度較低，不考慮自然空化作用。

綜上所述，目前的研究集中在航行體單獨出水的空泡演化、流體動力及運動特性，對多發航行體齊射出水的研究較少，仍處于探索階段，而且已有齊射出水的研究涉及自然空化作用的較少。因此，本文基于VOF 多相流模型，通過重疊網格技術開展了雙發射彈齊射出水問題的數值模擬研究，獲得了齊射出水過程的超空泡演化特性、射彈的運動軌跡、偏轉角變化、速度衰減曲線和阻力特性，分析了超空泡流場的干擾機理以及發射時差的影響。研究結果可為水下多彈齊射安全性分析和方案設計提供相關的參考。

1 數學模型及數值方法

1.1 控制方程

采用均質多相流模型來處理射彈出水過程誘導的多相流動，即將水、空氣、水蒸氣三相當作單一介質的混合相，混合相的密度由各相的密度和體積率共同決定，各相之間不存在速度滑移，并且具有相同的壓力場和速度場，流場計算時只對混合相求解N-S 方程和湍流方程。

混合相的連續性方程和動量方程分別為：

式中：ρm為 混合相密度；xi、xj為笛卡爾坐標分量；ui、uj分別為笛卡爾坐標系中的速度分量；p為壓力；μm為混合相動力黏度；gi為重力加速度的分量。

混合相密度的表達式為：

式中：αv、αg分 別為水蒸氣、空氣的體積分數；ρv、ρl、ρg分別為水蒸氣、水和空氣的密度。

混合相的動力黏度的表達式為：

式中：μv、μl、μg分別為水蒸氣、水和空氣的動力黏度。

采用SSTk-ω湍流模型[24]對雷諾平均方程提供湍流封閉，該模型考慮了湍流剪切應力傳輸，綜合了k-ε模型在外部區域模擬和標準的k-ω湍流模型在近壁計算的優點，能夠精確模擬逆壓梯度下出現的流動分離現象。

渦黏性系數 μt的表達式為：

湍動能k和比耗散率 ω的輸運方程如下：

式中：σk、σω是 湍動能k、比耗散率 ω的湍流普朗特數；Gk、Gω分別為湍動能、比耗散率產生項；Sk與Sω為自定義項；Yk、Yω分別為湍動能、比耗散率的發散相；Dω為正交發散項。

空氣和水蒸氣相的體積分數的輸運方程為：

空化問題的求解采用基于Rayleigh-Plesset 氣泡方程建立的Schnerr-Sauer 空化模型：

式中：RB為 氣核的半徑；pv為水的飽和蒸氣壓力；n0為單位液體體積空泡個數，取值為1×1011。

1.2 6DOF 剛體運動模型

射彈在流場的運動通過6DOF 剛體運動模型[25]來求解。射彈運動可分為質心的平動和繞質心的轉動，在慣性參考坐標系下建立射彈的平移運動控制方程：

式中：VG為射彈質心的平移速度矢量；m為射彈的質量；FG是射彈受到的外力矢量。

射彈的轉動在剛體坐標下的控制方程為：

式中：ωp為射彈角速度矢量；I為射彈的慣性張量；M為射彈外力矩矢量。

1.3 計算模型、邊界條件及數值方法

本文數值計算采用的射彈模型為參考文獻[10]的射彈模型，射彈為截錐型回轉體，由頭部錐臺段和圓柱段組成，彈體總長L=48 mm，頭部空化器直徑D0=3 mm，錐臺段L0=18 mm，如圖1 所示，質量m為2.98 g。

圖1 射彈的幾何模型Fig.1 Geometric model of projectile

計算域和邊界條件如圖2 所示，計算域長35D、寬20D、高130D，空氣域高度為50D，水域高度為80D，首發射彈和次發射彈的質心距離自由液面的高度h均為71.2D，兩射彈中軸線之間的距離為4D。計算域底部為壓力入口，壓強為105 704 Pa，計算域頂部為壓力出口，壓強為標準大氣壓，4 個側面均為壓力出口，壓強依據水深環境指定，其壓力梯度分布由用戶自定義場函數定義。首發射彈和次發射彈的表面均設定為無滑移壁面條件。定義處于兩射彈的中間區域為射彈的內側，其余為射彈的外側。

圖2 計算域及邊界條件設置Fig.2 Computational domain and boundary condition setting

兩發射彈齊射出水過程中，首發射彈和次發射彈的發射時間間隔為發射時差 Δt。當 Δt=0 時，被稱為同步發射；當 Δt≠0 時，被稱為異步發射，異步發射過程中首發射彈先發射，隨后次發射彈發射。在數值模擬中通過UDF 控制次發射彈在水下運動的起始時間從而實現兩發射彈以不同發射時差出水，出水過程中，首發射彈和次發射彈的初速度均為150 m/s。

計算域網格劃分使用重疊網格技術，該技術對計算流場中多物體運動比較有效。計算區域分為背景網格區域和子網格區域，整個流場為背景網格區域，由于首發射彈和次發射彈分別獨立運動，因此需要劃分兩個子網格區域，兩個子網格區域大小相同，均為包裹射彈的圓柱，其長度為10D，直徑為3D。圖3 為計算域的網格劃分示意圖，均采用結構化網格進行劃分，背景區域中對航行體運動區域和自由液面附近的網格進行局部加密，以便精確捕捉超空泡界面和自由液面的變形，子區域中對射彈壁面附近的網格進行局部加密。計算中，設置子網格1 計算域與首發射彈同步運動以及子網格2 計算域與次發射彈同步運動，然后結合 6DOF 算法可以實現兩個射彈各自獨立運動的六自由度運動求解。計算模型中壓力與速度耦合的求解采用Coupled 算法，壓力場和空間離散采用PRESTO!格式，數值計算中時間步長選擇為1×10-6s。

圖3 網格劃分示意圖Fig.3 Schematic of computational domain grids

1.4 無量綱數

對文中的物理量進行無量綱化：

式中：t為時間；v0為 射彈的初速度；Δt為兩發射彈的發射時差；vz為射彈的豎直速度；z為射彈質心的豎直位移；x為射彈質心的水平位移。

定義無量綱超空泡半徑、無量綱位置：

式中：r為在zs位置處的超空泡輪廓半徑，設定空泡外側輪廓到中軸線距離為負，空泡內側輪廓到中軸線距離為正。

式中：p0為標準大氣壓；Fz為射彈在z軸方向受到的阻力；A0為空化器的沾濕面積，這里取射彈頭部的橫截面積。

2 數值方法驗證

先進行網格無關性驗證，針對發射時差 Δ=0 的兩發射彈齊射出水采用了3 種不同數量的網格進行數值模擬，網格數分別為157 萬（粗糙尺度網格）、274 萬（中等尺度網格）、362 萬（精細尺度網格），不同網格尺度下首發射彈的無量綱豎直速度的變化如圖4 所示。從圖中可以看出，隨著網格密度的增加，中等尺度網格和精細網格計算所得的結果基本一致，綜合考慮計算精度和計算效率，故選取了中等尺度網格開展數值計算。

圖4 不同網格數量下無量綱豎直速度衰減曲線Fig.4 Dimensionless vertical velocity attenuations of different grid numbers

采用文獻中[10]中的超空泡射彈出水的實驗數據進行數值方法的有效性驗證，實驗中的射彈模型與本文數值計算采用的射彈模型相同（見圖1），實驗中射彈的出水初速度為98 m/s。圖5 表示的是超空泡射彈出水過程的實驗和數值模擬結果對比，圖中相鄰圖片的無量綱時間間隔為24.5，通過對比可以發現，整體上數值模擬結果和實驗結果吻合度較好，但數值模擬中水在射彈尾部兩側的附著現象與實驗結果不太一致，其原因為射彈尾部兩側的水破碎成液滴，繼而受到重力作用下降，而數值模擬中對此破碎過程難以準確模擬。另外在數值模擬中尾空泡的長度要長于實驗觀測，這是由于超空泡尾部閉合區域是充滿蒸氣、液滴和旋渦的多相流湍流區，數值模擬很難精確模擬超空泡尾部滯止及潰滅。為了進一步驗證數值模擬方法的準確性，將射彈的無量綱豎直位移的數值結果與實驗數據進行對比，如圖6 所示。從圖中可以看出射彈無量綱豎直位移的數值結果與實驗數據具有較好的一致性，其最大誤差約為3.9%，誤差結果在可接受的范圍。因此，經過上述的比較，說明本文采用的數值模擬方法是有效的。

圖5 超空泡射彈出水過程的實驗[10]和數值模擬結果對比Fig.5 Comparison of the water-exit process between experimental and numerical simulation results

圖6 無量綱豎直位移變化的實驗和數值模擬結果對比Fig.6 Comparison of dimensionless vertical displacement between experimental and numerical simulation results

3 結果與討論

3.1 空化流場特性分析

圖7 給出了4 種發射無量綱時差條件下雙發射彈齊射出水過程的水相圖。為便于討論，定義處于兩射彈中間區域的超空泡壁面為超空泡內側壁面，其余為超空泡外側壁面。齊射過程中，兩發射彈都經歷了超空泡的生成、發展、潰滅，在出水階段，伴隨著自由面的抬升和隆起，劇烈噴濺出水花等現象。對于同步發射，兩發射彈的超空泡內側壁面同時擴張時，由于超空泡內側壁面附近的流體形成了對流，抑制了超空泡向兩射彈中間區域的擴張，而外側壁面自由擴張，導致單個超空泡失去了對稱性；隨著時間的發展，兩個超空泡在尾部區域相互靠攏，但兩個超空泡在空間上呈現了良好的鏡面對稱特征。不同于同步發射的工況，Δ=25 的異步發射過程超空泡形態有較大差別，從圖7（b）中可知，首發射彈先進入水中運動，其誘導的超空泡迅速擴張，空泡輪廓正常發展；隨后次發射彈緊隨而來，在次發射彈流場的干擾下，首發射彈超空泡尾部區域的內側壁面發生收縮并且空泡尾部潰滅的速度加快，而次發射彈超空泡內側壁面發生膨脹。結合圖8，根據Logvinovich 的空泡截面獨立擴張原理[5]，空泡的每個橫截面的膨脹擴張只取決于空泡外界和空泡內部的壓差，以及頭部空化器經過該截面時的速度、空化器的尺寸及阻力。對于本工況，影響空泡截面擴張的決定性因素是空泡外界和空泡內部的壓差。射彈在水下高速運動，射彈頭部區域生成高壓區，誘導生成的超空泡內部為低壓區，其值等于水的飽和蒸氣壓（約為3 540 Pa）。對于首發射彈超空泡，在次發射彈頭部高壓區持續壓縮下，處于高壓區附近的首發射彈空泡內側壁面率先收縮；對于次發射彈超空泡，由于內側區域的空泡內外壓力差相對外側較小，導致次發射彈空泡內側壁面擴張的幅度更大。對于 Δ=50 和Δ=75，它們的空泡演化過程和 Δ=25相似，這里不再贅述。

圖7 不同發射無量綱時差下射彈齊射出水的水相圖Fig.7 Supercavity evolutions of two projectiles exiting water in underwater salvo for different launch dimensionless time intervals

圖8 Δ =25 時，超空泡演化示意圖Fig.8 Diagram of supercavity evolutions for Δ =25

圖9 表示的是 Δ=75、次發射彈典型時刻的無量綱壓力云圖。結合圖9（d）可以看出，首發射彈出水后剝離在自由面下的空泡發生收縮潰滅，在潰滅的過程中空泡上端和下端產生了局部高壓，此時次發射彈處于該空泡的附近，引起次發射彈頭部附近的壓力場與空泡潰滅產生的局部高壓發生耦合作用。在首發射彈脫體空泡潰滅造成的局部高壓作用下，次發射彈超空泡在相同水深的空泡壁面率先發生收縮，最終出現了“頸縮”現象”，如圖7（d）中=175 所示。

圖9 Δ =75，次發射彈典型時刻的無量綱壓力云圖Fig.9 Dimensionless pressure distributions at typical time of the second projectile for Δ =75

由于射彈在水下高速航行時的空化數很小，超空泡的全局尺寸很大，但只有射彈附近的空泡形態對流體動力特性有影響。因此，本文選擇超空泡前沿輪廓（2 倍彈長）作為研究對象。圖10 給出了不同發射時差下射彈在特征位置處的前沿輪廓對比，其中圖10（a）為首發射彈在無量綱豎直位移=61 時的超空泡前沿輪廓對比，圖10（b）為次發射彈在=23 時的超空泡前沿輪廓對比。在圖10（a）中，同步發射條件下的首發射彈超空泡外側前沿輪廓曲率大于內側的，而異步發射條件下的首發射彈的超空泡前沿輪廓對稱性較好，這是由于同步發射時，超空泡的內側受到相鄰射彈的排擠作用，空泡內側壁面的擴張受到抑制，而空泡外側壁面自由擴張，異步發射時，首發射彈超空泡的前沿部分不受次發射彈流場的干擾，空泡前沿的外側和內側壁面均自由發展，超空泡前沿輪廓左右基本對稱。從圖10（b）中可知，對于異步發射，超空泡外側前沿輪廓的曲率小于內側前沿輪廓的曲率，其原因上文已給出分析，并且次發射彈超空泡內側前沿輪廓曲率隨著發射時差的增大而減小。

圖10 特征位置處的射彈超空泡前沿輪廓對比Fig.10 Comparison of the front part of supercavity profiles at the feature position

圖11 Δ =0，出水階段特征位置的射彈表面無量綱壓力分布及超空泡輪廓圖Fig.11 Dimensionless pressure distributions on the projectile surface and supercavity profiles at the feature position of water exit stage for Δ =0

3.2 運動特性分析

圖12 為出水過程射彈的運動軌跡。從圖12（a）中可知，異步發射的3 種發射時差下，首發射彈的運動軌跡幾乎為一條豎直線，這是由于異步發射時，次發射彈主要干擾的是首發射彈超空泡的尾部部分，對射彈附近的空泡形態幾乎無影響，而決定射彈流體動力的是射彈附近的空泡形態，這也表明異步發射條件下次發射彈流場對首發射彈的運動基本無干擾。對于同步發射，首發射彈的質心先向內側偏移，隨后向外側偏移，當完全出水后，首發射彈質心的軌跡基本為直線變化。由圖12（b）可知，對于同步發射，次發射彈和首發射彈的運動具有對稱性。對于異步發射，由于受到首發射彈超空泡流場的干擾，次發射彈在內外兩側壓差作用下質心先向內側偏移，在穿越自由面階段有所波動。異步發射條件下，次發射彈質心的最大無量綱水平位移隨發射時差的增大而減小。

圖12 出水過程射彈的運動軌跡Fig.12 Trajectories during the projectiles exiting water

圖13 為不同發射無量綱時差下射彈偏轉角隨無量綱豎直位移的變化，定義偏轉角θ為當前射彈軸線與豎直方向z軸的夾角，取射彈發生逆時針偏轉時為正值。從圖中可知，同步發射條件下，由于兩發射彈頭部附近區域均產生了高壓區，壓力場的作用下產生了方向相反的力矩，首發射彈逆時針偏轉，次發射彈順時針偏轉，即兩發射彈的運動向著兩者頭部遠離、尾部靠近的方向偏轉。此種工況下，首發射彈和次發射彈的偏轉角大小隨著無量綱豎直位移的增大而逐漸增大，最大值達到了3.1°；當射彈出水后，由于超空泡的不對稱潰滅（圖11 可知），使射彈產生了相反方向的俯仰力矩，其偏轉角開始減小。對于異步發射，出水過程中首發射彈的偏轉角基本為0°，而次發射彈由于受到首發射彈超空泡尾部低壓流場的影響發生逆時針偏轉，隨著次發射彈向上運動距離逐漸增大，在出水后，其偏轉角發生波動。次發射彈的最大偏轉角隨發射時差的增大而減小，當發射時差 Δ=75 時，次發射彈的偏轉角幾乎維持在0°附近。

圖13 不同發射無量綱時差下射彈偏轉角隨無量綱豎直位移的變化Fig.13 Variations of deflection angles of projectiles with dimensionless vertical displacement for different launch dimensionless time intervals

圖14 表示出水過程中射彈無量綱豎直速度衰減曲線。圖14（b）中定義為當地無量綱時間，時間零點為次發射彈發射時刻。從圖中可以看出，運動過程中4 種發射時差的首發射彈和次發射彈的無量綱豎直速度之間的差異較小。從放大圖中可知，同步發射條件下首發射彈和次發射彈的無量綱豎直速度衰減略快，這是因為同步發射相對異步發射，首發射彈發生了較大偏轉，阻力增加，產生了更大的能量損失，這表明同步發射會對超空泡射彈的減阻產生不利影響。

圖14 出水過程射彈無量綱豎直速度衰減曲線Fig.14 Attenuation curves of dimensionless vertical velocity during the projectiles exiting water in vertical direction

圖15 表示出水過程中射彈的阻力系數變化。射彈起始運動時，射彈大部分還處于沾濕狀態，受到水的阻力較大，隨著超空泡的生成，射彈被超空泡包裹，只有頭部沾濕，射彈受到水的黏性阻力幾乎為0，主要為壓差阻力，此過程射彈的阻力系數比較穩定。隨著射彈穿越自由面進入空氣中，射彈頭部周圍的流體介質由水變成空氣，密度值顯著降低，導致射彈受到的阻力變成了量級很小的空氣阻力值，這時阻力系數不妨標記為0。首發射彈和次發射彈出水過程的阻力系數曲線的變化趨勢基本一致；相比異步發射，同步發射條件下射彈阻力系數的平均值略大。需要注意的是，Δ=75 時，次發射彈的阻力系數從=70（對應=145）時刻起出現增大，其原因為次發射彈受到附近潰滅空泡產生的局部高壓影響，引起次發射彈受到的壓差阻力增大（見圖9）。

圖15 出水過程中射彈的阻力系數變化Fig.15 Variations of drag coefficients during the projectiles exiting water

4 結論

本文基于求解N-S 方程的VOF 方法，采用重疊網格技術對高速射彈齊射出水問題進行了數值模擬研究，獲得的主要結論如下：

1）同步發射出水時，兩發射彈的超空泡流場互相影響，超空泡內側擴張受到抑制，引起射彈超空泡外側前沿輪廓曲率大于內側，在出水階段出水超空泡發生了非對稱性潰滅。異步發射出水時，首發射彈超空泡前沿輪廓基本對稱，而次發射彈超空泡前沿輪廓內側壁面發生膨脹，失去了對稱性，隨著發射無量綱時差的增大，次發射彈超空泡內側前沿輪廓曲率變小。

2）出水過程中，同步發射條件下兩發射彈的彈道穩定性較差，兩射彈先發生向外側的偏轉，隨后向內側偏轉，其偏轉角的最大值達到了3.1°。異步發射條件下首發射彈有很強的彈道穩定性，運動軌跡沿豎直方向向上運動，偏轉角基本為0°；次發射彈在壓差作用下向內側偏轉，運動軌跡也向內側偏移，運動過程中次發射彈的最大無量綱水平位移和最大偏轉角隨發射時差的增大而減小，當發射無量綱時差為70 時，次發射彈的偏轉角幾乎為0°。

3）相比異步發射出水，同步發射出水射彈的無量綱豎直速度衰減較快，降低了射彈的減阻性能。