李治濤， 趙世平， 盧丙舉， 于 勇

(1. 中國船舶集團有限公司第七一三研究所，鄭州 450015；2. 河南省水下智能裝備重點實驗室，鄭州 450015； 3. 北京理工大學 宇航學院，北京 100081)

射彈高速入水會產(chǎn)生空化現(xiàn)象，在真實的海況條件下，還會受到波浪等因素的影響。射彈在波浪中入水，除受浮力、自身重力、砰擊載荷外，還要考慮波浪載荷對射彈的影響。波浪載荷對射彈入水過程結(jié)構(gòu)的安全性以及彈體入水彈道的穩(wěn)定性均有非常重要的影響。

儲慧林[1]基于流體體積函數(shù)(volume of fluid, VOF)法建立了二維數(shù)值水槽，研究了波浪相位、波高對魚雷入水過程砰擊載荷、彈道和速度的影響。王平等[2]基于VOF法，建立了楔形體自由入水模型，研究了入水時楔形體所受的載荷、自由液面的變化以及物面壓強分布等。楊衡等[3]采用非線性雙漸進法，研究了三維圓柱體在波浪、洋流下入水運動速度及軌跡。王文華等[4]采用二維計算流體動力學(computational fluid dynamics,CFD)方法，數(shù)值模擬了圓柱在規(guī)則波中的入水過程，研究了波浪參數(shù)(波高、周期和入水相位)對圓柱入水動力學特性的影響。金禹彤等[5]用CFD方法，數(shù)值模擬了楔形體波浪入水過程，研究了規(guī)則波波形位置對楔形體受力特性和運動姿態(tài)變化的影響。鄒麗等[6]通過具有造波功能的水槽，開展了兩種曲面楔形體分別在靜水和波浪中入水砰擊試驗，記錄了入水過程楔形體的加速度和壓力變化，用高速攝像系統(tǒng)捕捉砰擊過程的射流飛濺現(xiàn)象。侯昭[7]采用大渦模擬方法，數(shù)值模擬了圓柱低速波浪入水，研究了入水空泡演變機理以及漩渦特征等。郭張霞等[8]用數(shù)值方法對回轉(zhuǎn)體波浪入水進行了研究，波浪對回轉(zhuǎn)體的空泡特征、頭部壓力以及誘導氣相速度產(chǎn)生影響。朱仁慶等[9]用數(shù)值方法對不同剛度三維楔形體波浪入水進行了研究，彈性效應在一定程度上會減緩砰擊應力。劉邦威[10]基于勢流理論和邊界元法開發(fā)了物體三自由入水過程的全非線性數(shù)值模擬程序，對三維圓錐波浪入水進行了研究，分析了波高對物體入水過程的影響。史崇鑌等[11]用CFD法對跨介質(zhì)航行器波浪入水過程進行了數(shù)值模擬，分析了有無波浪和波高對航行器入水運動特征的影響。趙慶新等[12]用CFD法對水下運載器波浪入水過程進行了數(shù)值模擬，研究了波浪參數(shù)和風速等對水下運載器受力的影響。陳晨等[13]用能夠考慮氣、汽、液三相可壓縮性數(shù)值方法對射彈跨聲速入水初期多相流場進行了數(shù)值模擬，并分析了入水高度的影響。

對波浪入水問題，目前研究較多的是楔形體、圓柱、航行器、圓錐等，高速旋轉(zhuǎn)射彈波浪入水研究較少。本文用定義輸入邊界造波法、重疊網(wǎng)格技術(shù)和6自由度算法對某一高速旋轉(zhuǎn)射彈波浪入水過程進行數(shù)值模擬，研究了波峰、波節(jié)1、波谷、波節(jié)2等波面不同位置入水彈體空化形態(tài)、彈道特性和流體動力載荷特性的變化規(guī)律，研究結(jié)果對復雜海況下射彈入水問題具有一定的參考價值。

1 數(shù)值計算方法

本文采用均質(zhì)多相流理論中的流體體積函數(shù)模型對空氣、水蒸氣、水組成的多相流動系統(tǒng)進行描述，流動被看作單一介質(zhì)的流動系統(tǒng)，不同流動組分共享一套動量方程，通過計算得到單元內(nèi)各相流體介質(zhì)的體積分數(shù)，從而確定不同相在流動系統(tǒng)中的分布。

1.1 控制方程組

描述混合介質(zhì)流動的N-S方程組為

(1)

式中：v為速度分量；ρm=αlρl+αvρv+αgρg和μm=αlμl+αvμv+αgμg分別為混合介質(zhì)密度和動力黏度；ρl，ρv和ρg分別為水、水蒸氣和空氣的密度；μl，μv和μg分別為水、水蒸氣和空氣的動力黏性系數(shù)；μt為湍流黏性系數(shù)；αl，αv和αg分別為水、水蒸氣和空氣相的體積分數(shù)。對于任意控制體均有如下關(guān)系成立

αl+αv+αg=1

(2)

體積分數(shù)的變化滿足如下輸運方程

(3)

數(shù)值計算中采用剪切應力輸運k-ω湍流模型(shear-stress-transportk-ω, SSTk-ω)。該模型考慮了湍流剪切應力效應，對很多湍流流動問題的求解具有很高的精度和可信度。其數(shù)學方程為

(4)

式中：k為湍動能；ω為湍動能耗散率；υt=μt/ρ為湍流動力黏性系數(shù)；σk，γ，σω和σω2為湍流模型常數(shù)。

物體在水下高速運動時會引起空化現(xiàn)象，需要引入空化模型，用于描述水相與水蒸氣相之間的質(zhì)量輸運。本文選用Schnerr-Sauer空化模型，其表達式為

(5)

式中：RB=1×10-6為空化氣核半徑；αnuc=5×10-4為不可凝結(jié)氣體的體積分數(shù)；Fvap=50和Fcond=0.001為經(jīng)驗常數(shù)；飽和蒸汽壓pv取3 540 Pa。式(5)等式右側(cè)第一項是水蒸氣的生成項，第二項是水蒸氣的凝結(jié)項。

1.2 重疊網(wǎng)格技術(shù)

彈體入水是一個極強的非定常、非線性問題。彈體入水運動過程中彈體與水域接觸面的邊界位置，需要根據(jù)之前彈體邊界所受的載荷計算獲得。彈體每個時刻在水域中的位置均不同，用重疊網(wǎng)格技術(shù)可以模擬彈體實時運動引起的流域形狀隨時間變化。

重疊網(wǎng)格法將彈體入水過程復雜的流動區(qū)域分成流場區(qū)域和彈體周圍區(qū)域等兩個簡單的子區(qū)域。各子區(qū)域中的網(wǎng)格單獨生成，彼此存在著重疊或嵌套關(guān)系，流場信息通過插值在重疊區(qū)邊界進行匹配和耦合，各子區(qū)域無須進行繁瑣的拓撲分區(qū)，降低了網(wǎng)格生成難度[14-15]。本文在后續(xù)計算中流場區(qū)域即背景網(wǎng)格，彈體周圍區(qū)域即部件網(wǎng)格，二者均采用結(jié)構(gòu)網(wǎng)格。

1.3 6自由度算法

計算過程中部件網(wǎng)格區(qū)域隨彈體同步運動，部件和背景網(wǎng)格之間的挖補插值操作會在每一個時間步進行，彈體運動采用6自由度算法計算。

采用地面坐標系和彈體坐標系來描述彈體的運動規(guī)律。地面坐標系OXYZ與CFD計算采用的坐標系一致如圖1所示。其中：OX軸在水平面內(nèi)；OY軸在縱平面內(nèi)，方向垂直向上，重力為-Y方向；OZ軸滿足右手坐標系。彈體坐標系原點位于彈體質(zhì)心，坐標系OX1軸與彈體軸線相一致，指向彈體尾端；OY1軸垂直于彈體橫向平面，方向垂直向上；OZ1滿足右手坐標系。

圖1 地面坐標系和彈體坐標系Fig.1 Ground coordinate system and body coordinate system

射彈質(zhì)心處的速度為V，根據(jù)動量定理在地面坐標系中有

(6)

式中，F(xiàn)為作用于射彈上的合外力。依據(jù)慣性系(地面坐標系)與動坐標系(彈體坐標系)中速度導數(shù)關(guān)系有

(7)

將式(7)代入式(6)，得

(8)

式中，F(xiàn)x，F(xiàn)y和Fz為作用于彈體的外力矢量在彈體坐標系中的分量。

慣性系下彈體旋轉(zhuǎn)運動的動量矩方程為

(9)

式中：M為合外力對彈體質(zhì)心的力矩；J為射彈的轉(zhuǎn)動慣量。

根據(jù)慣性系(地面坐標系)中與動坐標系(彈體坐標系)中動量矩導數(shù)的關(guān)系有

(10)

將式(10)代入式(9)，得

(11)

式中，Mx，My和Mz為外力對射彈質(zhì)心力矩M在彈體坐標系中的分量。

2 數(shù)值方法驗證

用Schnerr-Sauer空化模型、SSTk-ω湍流模型和重疊網(wǎng)格技術(shù)，建立了高速射彈跨介質(zhì)入水空化多相流場數(shù)值仿真模型。為了驗證該空化多相流場數(shù)值仿真模型的可靠性和準確性，首先對其進行了驗證。用文獻[16]中的平頭圓柱高速入水問題作為驗證算例。

2.1 平頭圓柱模型

采用郭子濤研究中所用的平頭圓柱，其參數(shù)如表1所示。

表1 平頭圓柱參數(shù)Tab.1 Parameters of flat-headed cylinder

2.2 計算網(wǎng)格及邊界條件

部件網(wǎng)格為包含平頭圓柱的圓柱形區(qū)域，平頭圓柱位于該計算域的中心，計算域直徑為15倍平頭圓柱長度即0.381 m，高度為11倍平頭圓柱長度即0.279 4 m，網(wǎng)格數(shù)目為310 374，如圖2(a)所示；背景網(wǎng)格區(qū)域為長方體，長為0.76 m，寬為0.76 m，高為1.20 m，網(wǎng)格數(shù)目為1 110 262，其中水域高為1.00 m，空氣域高為0.20 m；背景網(wǎng)格與部件網(wǎng)格初始化后在Z=0截面分布，如圖2(b)所示，Z軸滿足右手坐標系。

邊界條件設(shè)置如圖3所示。計算域X軸最大值與最小值的邊界面、Z軸最大值與最小值的邊界面均設(shè)為對稱邊界條件；Y軸最大值邊界面設(shè)置為壓力出口，壓強大小與標準大氣壓一致；Y軸最小值邊界面設(shè)置為壓力出口，壓強為對應水深的靜水壓力；平頭圓柱邊界設(shè)置為Wall；部件網(wǎng)格外邊界設(shè)置為Overset邊界；計算時間步長為1×10-6s。

圖2 平頭圓柱計算域網(wǎng)格分布Fig.2 Grid distribution in computational domain of flat-headed cylinder

圖3 邊界條件設(shè)置Fig.3 Boundary condition setting

2.3 入水速度與位移驗證

圖4分別給出了平頭圓柱體跨介質(zhì)入水過程中，彈體入水位移曲線和速度衰減曲線。從圖4中可以看出數(shù)值計算結(jié)果與試驗具有較好的一致性。驗證了用Schnerr-Sauer空化模型、SSTk-ω湍流模型和重疊網(wǎng)格技術(shù)所建入水空化多相流場數(shù)值仿真模型的可靠性和準確性。

圖4 平頭圓柱入水的位移與速度變化曲線Fig.4 Displacement and velocity curves of flat-headed cylinder entry water

3 數(shù)值造波

定義輸入邊界造波是根據(jù)行波的解析解，給定造波邊界處的值，從而實現(xiàn)造波的一種方法[17-22]。它適用于規(guī)則波、非規(guī)則波等多種波浪模型[23]，并且還能夠適應波浪和洋流耦合，在波浪和海洋工程領(lǐng)域應用廣泛。

3.1 入射邊界條件

采用定義輸入邊界造波法考慮到第4章射彈入水的三維數(shù)值模擬，造波采用三維計算，忽略Z向(側(cè)向)的波浪運動，即只考慮二維波浪情形。

對于二階Stokes波浪，入射邊界處的速度和波高滿足式(12)。其X和Y方向速度與波高分別為

(12)

式中,T，k，H和d分別為波浪周期、波數(shù)、波高和靜水深。

定義輸入邊界造波法的邊界條件設(shè)置，如圖5所示。用Fluent中的宏DEFINE_PROFILE(name,t,i)定義入射邊界處流體的速度；用宏DEFINE_PROFILE(name,t,i)定義入射邊界處流體的體積分數(shù)，達到控制入射邊界處流體高度，液相狀態(tài)下水的體積分數(shù)為1，空氣的體積分數(shù)為0。

圖5 定義輸入邊界造波法示意圖Fig.5 Schematic diagram of wave-generation method of defining inlet boundary conditions

3.2 海綿層阻尼消波

波浪傳播至水槽末端，如不進行消波處理，會造成波浪的反射。反射波浪的存在會與入射波浪形成波浪疊加，影響數(shù)值水槽的波形和流場。海綿層阻尼消波是在特定的計算區(qū)域內(nèi)添加阻尼項以達到消弱或消除該區(qū)域波動。本文將計算域邊界前1～2倍波長的范圍設(shè)置為消波區(qū)域，在動量方程中添加阻尼項，通過Fluent中的宏DEFINE_SOURCE(mom_source,c,t, dS, equ)實現(xiàn)消波。動量方程中添加人工黏性為

(13)

式中：α(x)為消波系數(shù)，采用線性分布的形式；x1和x2分別為海綿層消波的起始坐標和末端坐標；β為經(jīng)驗常數(shù)，與波浪特征等因素相關(guān)，一般取10。

3.3 造波計算及結(jié)果

數(shù)值模擬5級浪二階Stokes波，該波的波長、波高和周期分別為40 m，3 m和5.06 s。造波區(qū)域長160 m、寬5 m、高20 m，其中水深15 m、空氣域高5 m，坐標原點位于靜液面對稱中心處。計算采用的網(wǎng)格分別在靜液面上下半波高附近進行了加密，網(wǎng)格數(shù)目為2 381 727，如圖6所示；造波采用的邊界條件與圖5相同，入口邊界條件的速度剖面、壓強分布、體積分數(shù)等參數(shù)均通過自定義函數(shù)(user-defined function, UDF)載入。計算采用隱式VOF模型、SSTk-ω湍流模型、含分離算子的隱式壓力(pressure-implicit with splitting of operators, PISO)算法，時間步長為0.01 s；當有較為穩(wěn)定的波形出現(xiàn)時，逐步降低時間步長一直到1×10-6s，以適應第4章射彈高速入水的時間步長。采用水相體積分數(shù)為0.5的位置作為波面位置，二階Stokes波面歷時演化過程如圖7所示。

圖6 數(shù)值水槽網(wǎng)格Fig.6 Numerical flume grid

圖7 二階Stokes波數(shù)值生成歷時演化過程Fig.7 Evolution of second-order Stokes wave generation

仿真所得T=50 s波形與理論波形對比，如圖8所示。圖中X坐標從20～60 m數(shù)值計算時設(shè)置為消波區(qū)即圖中灰色區(qū)域，受消波與反射波浪的影響X坐標從20～80 m內(nèi)仿真所得波形與理論波形差別較大；圖中X坐標從-80～0 m內(nèi)數(shù)值計算所得波形與理論波形吻合較好。計算所得波形能較好地反應波浪波面變化，可用于第4章高速射彈波浪入水仿真分析。

圖8 數(shù)值波形和理論波面對比Fig.8 Comparisons of numerical and theoretical wavefront

4 高速旋轉(zhuǎn)射彈波浪入水仿真分析

選取高速射彈在波浪波面四個不同位置入水進行數(shù)值仿真研究。波面入水位置如圖9所示，分別為波峰、波節(jié)1、波谷、波節(jié)2。

計算所用射彈模型如圖10所示，彈體長為Lt，頭部直徑為Dn，彈體最大直徑為Dt，彈體細長比為8.746，彈體質(zhì)心到頭部距離與彈體長度的比值為48.37%。彈體初速度為408 m/s，名義入水角為25°(彈體軸線與水平面的夾銳角為名義入水角)，初始旋轉(zhuǎn)速度為7 500 r/min。

圖9 高速旋轉(zhuǎn)射彈波面不同入水位置Fig.9 Different water-entry positions of wave surface of high-speed spinning projectile

圖10 彈體幾何示意圖Fig.10 Geometric sketch of projectile

4.1 射彈波浪入水彈體網(wǎng)格與計算方法

高速旋轉(zhuǎn)射彈波浪入水，采用1.2節(jié)所述重疊網(wǎng)格技術(shù)，網(wǎng)格分為背景網(wǎng)格和部件網(wǎng)格。背景網(wǎng)格見圖6；部件網(wǎng)格區(qū)域是包圍彈體周圍一定區(qū)域，使彈體周圍有較好質(zhì)量且足夠精細的網(wǎng)格用于計算。采用三套不同大小流場區(qū)域的部件網(wǎng)格，它們的幾何尺寸如圖11所示，圖中D為彈體最大直徑。雖然三套部件網(wǎng)格流場區(qū)域大小不同，但它們的網(wǎng)格拓撲是一樣的。這里僅給出第3套網(wǎng)格分布如圖12所示，第1套和第2套網(wǎng)格與其類似不再附圖。圖12分別給出了過彈體重心軸截面的網(wǎng)格分布和其他兩個剖面的典型網(wǎng)格，這樣劃分的結(jié)構(gòu)化網(wǎng)格質(zhì)量很高。

圖11 部件網(wǎng)格流場流域示意圖Fig.11 Schematic diagram of flow field area of component grid

圖12 部件網(wǎng)格局部截面分布示意圖Fig.12 Schematic diagram of component grid local section distribution

3套不同大小流場區(qū)域的部件網(wǎng)格信息，如表2所示。背景網(wǎng)格尺寸大小的選取應與部件網(wǎng)格最大網(wǎng)格尺寸保持基本一致。這樣可以保證在背景網(wǎng)格和部件網(wǎng)格交界面上插值的精度，在提高部件(彈體)周圍流場分辨率的同時降低總的網(wǎng)格數(shù)目。

表2 部件網(wǎng)格信息Tab.2 Component grid information

圖13給出了3套不同部件網(wǎng)格，射彈以初速度為408 m/s、名義入水角7.2°和初始旋轉(zhuǎn)速度為7 500 r/min平靜水面工況入水，彈體質(zhì)心運動軌跡在XOY面的投影。從圖13中可以看出，受部件網(wǎng)格區(qū)域大小的影響第1套網(wǎng)格彈體在XOY平面比第2套和第3套網(wǎng)格更早緩慢抬頭，第2套和第3套網(wǎng)格彈體運動軌跡幾乎沒有差別。增加部件網(wǎng)格區(qū)域的尺寸勢必增大網(wǎng)格的數(shù)量，加大計算負擔，所以部件網(wǎng)格尺寸的選取需要綜合平衡計算精度和計算資源的取舍。事實上，當部件網(wǎng)格尺寸達到一定值后，再增加部件網(wǎng)格尺寸對計算結(jié)果基本不會有影響。經(jīng)過上述分析，又考慮到第3套網(wǎng)格僅比第2套增加了5萬多網(wǎng)格量，計算量增加不是太大，后續(xù)仿真計算中均采用第3套部件網(wǎng)格。

圖13 3套不同部件網(wǎng)格，彈體質(zhì)心運動軌跡在XOY面的投影Fig.13 Three sets of different component grid, projection of projectile centroid trajectory on XOY plane

將背景網(wǎng)格和部件網(wǎng)格同時導入Fluent軟件中，設(shè)置彈體初始位置和姿態(tài)，彈體不動，用第3章所述定義輸入邊界造波法進行數(shù)值造波，在造出穩(wěn)定波形后逐步減小時間步長，使其與射彈波浪入水時間步一致為1×10-6s；之后設(shè)置彈體6自由度運動，打開空化模型、三相流體等，進行射彈波浪入水仿真，直到計算結(jié)束。高速旋轉(zhuǎn)射彈入水位置處的波形(波峰、波節(jié)1、波谷和波節(jié)2等)可用式(12)中的波高預示獲得。

4.2 空化效應對高速射彈入水特征的影響

圖14給出了是否考慮空化效應，高速射彈波峰入水Z=0平面不同相的體積分數(shù)云圖。圖14(a)未考慮空化效應，即計算中未開啟1.1節(jié)Schnerr-Sauer空化模型，該圖中深灰色為水相、白色為空氣相；圖14(b)考慮空化效應，即計算中開啟1.1節(jié)Schnerr-Sauer空化模型，該圖中淺灰色為水相、深灰色為空氣相、白色為水蒸氣相。通過對比圖14中是否考慮空化效應可發(fā)現(xiàn)：未考慮空化效應射彈波峰高速入水，剛?cè)胨畷r彈尾攜帶部分空氣圖中白色，隨著入水深度增加彈尾空氣逐漸從彈尾斷裂溢出，入水9 ms以后彈尾幾乎沒有空氣相，入水后受水動力的作用彈體從6 ms后逐漸緩慢抬頭(見圖14(a))；考慮空化效應，射彈剛?cè)胨畷r彈體頭部產(chǎn)生了空化效應，即水相(圖中淺灰色)汽化為水蒸氣(圖中白色)，同時彈尾攜帶部分空氣(圖中深灰色)，由于空化效應在12 ms內(nèi)彈體仍保持在超空泡中直線航行。

圖14 是否考慮空化效應，高速射彈波峰入水Z=0平面不同相體積分數(shù)云圖Fig.14 Whether cavitation effect is considered, volume fraction contour of different phases in the Z=0 plane of wave crest entering water of high-speed projectile

圖15給出了是否考慮空化效應，彈體波峰入水其質(zhì)心運動軌跡。從圖15中可以看出，未考慮空化效應由于彈體高速旋轉(zhuǎn)其受到較大的側(cè)向水動力載荷作用，使其產(chǎn)生了側(cè)向位移，12 ms內(nèi)彈體質(zhì)心產(chǎn)生了約40 mm的側(cè)向位移；考慮空化效應，12 ms內(nèi)彈體質(zhì)心幾乎未產(chǎn)生側(cè)向位移。

圖16給出了是否考慮空化效應，彈體波峰入水其質(zhì)心運動軌跡在XOY平面的投影。從圖16中可以看出，未考慮空化效應，彈體入水一段時間后由于水動力載荷的作用其緩慢抬頭；考慮空化效應，彈體在12 ms內(nèi)沿彈體初始對稱軸直線航行。

圖15 是否考慮空化效應，彈體波峰入水其質(zhì)心運動軌跡Fig.15 Whether cavitation effect is considered, centroid trajectory of projectile of wave crest entering water

圖16 是否考慮空化效應，彈體波峰入水其質(zhì)心軌跡在XOY平面投影Fig.16 Whether cavitation effect is considered, projection of centroid trajectory of projectile of wave crest entering water on XOY plane

仿真計算時保存彈體三個方向的受Fx，F(xiàn)y和Fz，以及力矩中心在原點三個坐標方向的力矩Mx，My和Mz，輸出的力和力矩是關(guān)于全局坐標系即地面坐標系下的值，需要變換到彈體速度坐標系計算其阻力、升力和側(cè)向力，變換到彈體坐標系計算滾轉(zhuǎn)力矩、偏航力矩和俯仰力矩。

過彈體質(zhì)心全局坐標系與彈體坐標系的變換矩陣[24]為

(14)

式中，γ，ψ和?為彈體的三個歐拉角，即繞X軸的滾轉(zhuǎn)運動、繞Y軸的偏航運動和繞Z軸的俯仰運動。

通過全局坐標系下彈體質(zhì)心受力FG和變換矩陣Ω，可求得彈體坐標系下的受力FB

FB=ΩFG

(15)

將彈體坐標系下的受力FB向彈體速度坐標系變換可得彈體運動過程中受到的阻力、升力和側(cè)向力。將彈體坐標系下的載荷用0.5ρV2S進行無量綱處理可得運動過程中彈體的阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)，其中：ρ為水的密度；V為彈體初速；S為彈體最大橫截面積。

全局坐標系下力矩中心在原點時彈體受到力矩Mo=[MxMyMz]T，全局坐標系下過彈體質(zhì)心的彈體力矩為MG，過質(zhì)心彈體坐標系下的力矩為MB=[MrollMyawMpitch]T，其中，Mroll，Myaw和Mpitch分別為彈體的滾轉(zhuǎn)力矩、偏航力矩和俯仰力矩。通過過原點全局坐標系下的力矩可得過質(zhì)心全局坐標系下的力矩，然后通過變換矩陣可得彈體坐標系下的力矩，即

MG=Mo-RoG×FG
MB=ΩMG

(16)

式中,RoG為全局坐標系中坐標原點到彈體質(zhì)心矢量。將力矩用0.5ρV2SLt無量綱處理可得力矩系數(shù)。

圖17～圖22分別給出了是否考慮空化效應，波峰入水彈體阻力系數(shù)、升力系數(shù)、側(cè)向力系數(shù)、滾轉(zhuǎn)力矩系數(shù)、偏航力矩系數(shù)和俯仰力矩系數(shù)隨時間變化曲線。從圖17～圖22中總體可以看出，未考慮空化效應彈體入水后直接與水接觸，其受到的水動力載荷幅值整體大于考慮空化效應的仿真結(jié)果，且在數(shù)值仿真計算12 ms內(nèi)未考慮空化效應載荷幅值變化較為劇烈。這是由于彈體以相同速度在水中航行和在水蒸氣與空氣相混合氣體中航行所受的流體動力載荷幅值有非常大的差別，其在水蒸氣與空氣相混合氣體中所受載荷較小，特別是阻力，極大的提高了射彈在水下航行速度，受到國內(nèi)外水動力領(lǐng)域的廣泛關(guān)注和重視。綜合本節(jié)是否考慮空化效應對射彈高速波峰入水不同相體積分數(shù)云圖、運動軌跡以及載荷特征分析可知，空化效應是射彈高速入水過程一個非常重要的物理現(xiàn)象，對射彈高速入水過程的運動軌跡以及載荷有非常重要的影響，故在后續(xù)射彈高速入水數(shù)值仿真過程均考慮了空化效應。

圖17 是否考慮空化效應，彈體阻力系數(shù)變化曲線Fig.17 Whether cavitation effect is considered, curves of drag coefficient of projectile

圖18 是否考慮空化效應，彈體升力系數(shù)變化曲線Fig.18 Whether cavitation effect is considered, curves of lift coefficient of projectile

圖19 是否考慮空化效應，彈體側(cè)向力系數(shù)變化曲線Fig.19 Whether cavitation effect is considered, curves of lateral force coefficient of projectile

圖20 是否考慮空化效應，彈體滾轉(zhuǎn)力矩系數(shù)變化曲線Fig.20 Whether cavitation effect is considered, curves of rolling moment coefficient of projectile

圖21 是否考慮空化效應，彈體偏航力矩系數(shù)變化曲線Fig.21 Whether cavitation effect is considered, curves of yawing moment coefficient of projectile

圖22 是否考慮空化效應，彈體俯仰力矩系數(shù)變化曲線Fig.22 Whether cavitation effect is considered, curves of pitch moment coefficient of projectile

4.3 不同波面位置入水空化特性分析

圖23分別給出了波峰、波節(jié)1、波谷和波節(jié)2四種工況下彈體入水過程不同相體積分數(shù)云圖，圖中淺灰色為水相、深灰色為空氣相、白色為水蒸氣相。從圖23中可知：由于波面的存在，彈體不同位置入水的有效入水角(彈體入水位置處波面的切線與彈體軸向的夾銳角為有效入水角；彈體波峰、波節(jié)1、波谷與波節(jié)2入水的有效入水角分別為25.003 2°，39.662 0°，25.009 8°與10.309 4°)發(fā)生了變化，與名義入水角相比，波節(jié)2有效入水角減小了14.690 6°，而波節(jié)1有效入水角增加了14.662 0°，波節(jié)2工況彈體失穩(wěn)最快。

圖23 高速射彈波面不同位置入水Z=0平面不同相體積分數(shù)云圖(圖中淺灰色為水相,深灰色為空氣相,白色為水蒸氣相)Fig.23 Contour of different phases volume fraction of high-speed projectile entry water at different positions of wave surface in Z=0 plane (in the figure, light gray is water phase, dark gray is air phase and white is water-vapor phase)

圖24分別給出了波峰、波節(jié)1、波谷和波節(jié)2四種工況下彈體的空泡形態(tài)演化過程，其中空泡界面取水蒸氣相體積分數(shù)為0.5的等值面。從圖23和圖24中可以看出，彈體入水經(jīng)歷了撞擊、空泡形成、超空泡航行和沾濕階段；對比不同波面位置入水彈體產(chǎn)生的空泡形態(tài)變化過程可以發(fā)現(xiàn)，波節(jié)2工況彈體姿態(tài)最先發(fā)生翻轉(zhuǎn)，波節(jié)1工況彈體下表面沾濕和姿態(tài)變化最慢，相同時間內(nèi)彈體姿態(tài)發(fā)生很大變化，由大到小依次是波節(jié)2、波峰、波谷和波節(jié)1。

圖24 高速射彈波面不同位置入水空泡形態(tài)圖Fig.24 Cavitation shape map of high-speed projectile entry water at different positions of wave surface

4.4 波面不同位置入水彈道特性分析

圖25和圖26分別給出了波峰、波節(jié)1、波谷和波節(jié)2入水彈體質(zhì)心在三維空間中的運動軌跡和彈體質(zhì)心運動軌跡在XOY平面上的投影。從圖25和圖26中可以看出，入水初期彈體質(zhì)心運動軌跡保持在彈體初始對稱軸方向，隨后發(fā)生偏轉(zhuǎn)；波節(jié)1工況彈體質(zhì)心軌跡側(cè)向偏移和縱對稱面內(nèi)的偏轉(zhuǎn)均很小，而其他三種工況在入水達到一定深度后彈體質(zhì)心位移變化較大；四種入水工況彈體側(cè)向偏轉(zhuǎn)程度由大到小依次為波節(jié)2、波峰、波谷和波節(jié)1，圖中仿真結(jié)束時波節(jié)2工況側(cè)偏了0.321 6 m、波峰工況側(cè)偏了0.061 6 m、波谷工況側(cè)偏了0.006 17 m、波節(jié)1工況側(cè)偏了0.000 638 m。

圖27和圖28分別給出了彈體運動過程中X向、Y向、Z向的速度衰減曲線和總速度大小衰減曲線，其中速度分量的正負與坐標軸方向一致。從圖27和圖28中可以看出，入水前期彈體基本處于超空泡航行，各方向速度變化基本一致，且側(cè)向Z向速度分量近似為0；入水中期不同工況彈體沾濕出現(xiàn)的時間和變化趨勢不同，速度變化量由大到小依次是波節(jié)2、波峰、波谷、波節(jié)1，在此階段逐漸產(chǎn)生了側(cè)向Z向速度分量；波節(jié)1工況彈體速度變化相對于其他三種工況得到了很大的延緩，同一時刻該工況彈體軌跡和速度變化量均小于其他工況。

圖25 高速射彈波面不同位置入水，彈體質(zhì)心運動軌跡圖Fig.25 Centroid trajectory of projectile entry water at different positions of wave surface

圖26 高速射彈波面不同位置入水，彈體質(zhì)心運動軌跡在XOY平面投影Fig.26 Projection of centroid trajectory of projectile on XOY plane entry water at different positions of wave surface

圖27 高速射彈波面不同位置入水，彈體質(zhì)心在X，Y，Z方向速度分量變化曲線Fig.27 Velocity component curves of projectile centroid in X,Y,Z direction entry water at different positions of wave surface

圖28 高速射彈波面不同位置入水，彈體質(zhì)心速度變化曲線Fig.28 Curves of centroid velocity of projectile entry water at different positions of wave surface

圖29～圖31分別給出了四種工況彈體入水運動過程中的滾轉(zhuǎn)角、偏航角與俯仰角變化規(guī)律。從圖29～圖31中可以看出：彈體滾轉(zhuǎn)角持續(xù)近似線性增長，不同工況彈體滾轉(zhuǎn)角的變化在入水前期幾乎一樣，入水后期產(chǎn)生少許差別；在計算時間內(nèi)彈體偏航角逐漸增大，偏航角增大量由大到小依次為波節(jié)2、波峰、波谷和波節(jié)1，圖中仿真結(jié)束時刻波節(jié)2偏航角增大了36.062 2°、波峰增大了16.144 8°、波谷增大了5.148 7°、波節(jié)1增大了2.979 2°；彈體俯仰角在彈體接觸水面前近似不變，接觸水面后俯仰角逐漸減小，在相同時間內(nèi)俯仰角減小量由大到小依次為波節(jié)2、波峰、波谷、波節(jié)1，且波節(jié)2工況彈體俯仰角變化最快；仿真結(jié)束時波節(jié)2工況彈體俯仰角減小了75.877 8°，彈體由入水姿態(tài)變?yōu)槌鏊藨B(tài)如圖23(d)和圖24(d)時刻23 ms所示，波峰工況彈體俯仰角減小了17.433 9°，波谷工況彈體俯仰角減小了6.185 4°，波節(jié)1工況彈體俯仰角減小了3.005 9°。

圖32～圖34分別給出了四種工況彈體入水運動過程的滾轉(zhuǎn)角速度、偏航角速度和俯仰角速度變化曲線。從圖32～圖34中可以看出，彈體滾轉(zhuǎn)角速度逐漸減小，減小程度由大到小依次為波節(jié)2、波峰、波谷、波節(jié)1，仿真結(jié)束時刻波節(jié)2工況滾轉(zhuǎn)角速度減小了142.586 6 rad/s，波峰工況減小了27.884 1 rad/s，波谷工況減小了5.496 3 rad/s，波節(jié)1工況減小了2.390 2 rad/s；彈體偏航角速度逐漸增大，仿真結(jié)束時刻波節(jié)2工況出現(xiàn)了極大值；同時可以看出彈體撞擊水面產(chǎn)生的俯仰角速度擾動比較明顯，而滾轉(zhuǎn)角速度和偏航角速度幾乎沒有受到影響。

圖29 高速射彈波面不同位置入水，彈體滾轉(zhuǎn)角變化曲線Fig.29 Curves of roll angle of projectile entry water at different positions of wave surface

圖30 高速射彈波面不同位置入水，彈體偏航角變化曲線Fig.30 Curves of yaw angle of projectile entry water at different positions of wave surface

圖31 高速射彈波面不同位置入水，彈體俯仰角變化曲線Fig.31 Curve of pitch angle of projectile entry water at different positions of wave surface

圖32 高速射彈波面不同位置入水，彈體滾轉(zhuǎn)角速度變化曲線Fig.32 Curves of rolling angular velocity of projectile entry water at different positions of wave surface

圖33 高速射彈波面不同位置入水，彈體偏航角速度變化曲線Fig.33 Curves of yaw angular velocity of projectile entry water at different positions of wave surface

圖34 高速射彈波面不同位置入水，彈體俯仰角速度變化曲線Fig.34 Curves of pitch angular velocity of projectile entry water at different positions of wave surface

4.5 波面不同位置入水載荷特性分析

圖35～圖37分別給出了四種工況彈體入水運動過程的阻力系數(shù)、升力系數(shù)和側(cè)向力系數(shù)的變化曲線。撞擊水面時彈體阻力系數(shù)產(chǎn)生瞬間的峰值，隨后趨于平緩；入水后期由于彈體表面的沾濕，彈體阻力系數(shù)迅速增大，其中波節(jié)2工況阻力系數(shù)最先增大，這與圖23和圖24波節(jié)2工況彈體最先失穩(wěn)結(jié)論一致；撞擊水面時彈體升力系數(shù)產(chǎn)生瞬態(tài)峰值后趨于平緩，入水后期由于彈體下表面的沾濕，彈體升力迅速增大，由于波節(jié)2工況彈體最先沾濕，故其升力系數(shù)最先增大；撞擊水面時彈體側(cè)向力系數(shù)不受影響為0，入水后期彈體側(cè)向力系數(shù)逐漸增大，由大到小依次是波節(jié)2、波峰、波谷和波節(jié)1，這與圖25所示彈體側(cè)向位移結(jié)論是一致的。

圖35 高速射彈波面不同位置入水，彈體阻力系數(shù)變化曲線Fig.35 Curves of drag coefficient of projectile entry water at different positions of wave surface

圖36 高速射彈波面不同位置入水，彈體升力系數(shù)變化曲線Fig.36 Curves of lift coefficient of projectile entry water at different positions of wave surface

圖37 高速射彈波面不同位置入水，彈體側(cè)向力系數(shù)變化曲線Fig.37 Curves of lateral force coefficient of projectile entry water at different positions of wave surface

圖38～圖40分別給出了四種工況彈體入水運動過程的滾轉(zhuǎn)力矩系數(shù)、偏航力矩系數(shù)和俯仰力矩系數(shù)的變化曲線。在彈體帶空泡運動初期，彈體的力矩系數(shù)趨于平緩且量很小；彈體表面沾濕后，彈體力矩系數(shù)發(fā)生很大的變化，由于彈體旋轉(zhuǎn)運動的存在，彈體的偏航力矩系數(shù)和俯仰力矩系數(shù)呈現(xiàn)正負交替變化。

圖38 高速射彈波面不同位置入水，彈體滾轉(zhuǎn)力矩系數(shù)變化曲線Fig.38 Curves of rolling moment coefficient of projectile entry water at different positions of wave surface

圖39 高速射彈波面不同位置入水，彈體偏航力矩系數(shù)變化曲線Fig.39 Curves of yawing moment coefficient of projectile entry water at different positions of wave surface

圖40 高速射彈波面不同位置入水，彈體俯仰力矩系數(shù)變化曲線Fig.40 Curves of pitch moment coefficient of projectile entry water at different positions of wave surface

4.6 高速射彈入水不同階段的載荷與彈道特征

不同學者對于運動體入水各階段的劃分不盡相同，May[25]將整個過程分為六個階段：沖擊階段、流動形成階段、空泡敞開階段、空泡閉合階段、空泡潰滅階段以及完全沾濕階段。本文高速射彈波面入水僅研究前五個階段，由圖23可知高速射彈波面不同位置入水四個階段(這里將空泡閉合與潰滅合為一個階段)時間劃分如表3所示。

文獻[26]指出：當t?V/g時，重力在許多高速入水問題中通常被忽略。高速射彈波面不同位置入水，僅分析前23 ms，由0.023 s?(409/9.80)s=41.735 s，故在高速射彈波面入水整個過程中重力因素可以被忽略。事實上在本文數(shù)值仿真過程中均考慮了重力效應。

通過圖25～圖40可知彈體沖擊和流動形成這兩個階段，彈體所受重要載荷以及由此引起彈道變化如圖41所示。這兩個階段彈體阻力急劇增大產(chǎn)生局部極大值，然后減小趨于穩(wěn)定值如圖41(a)所示，且波面不同位置入水四種工況穩(wěn)定值幾乎完全相同，約為44.73 kN；這兩個階段產(chǎn)生了沖擊升力載荷然后趨于穩(wěn)定值，如圖41(b)所示；這兩個階段使彈體產(chǎn)生了初始俯仰角速度擾動如圖41(e)所示。

表3 高速射彈波面不同位置入水四個階段時間劃分Tab.3 Time division of four stages of high-speed projectile entry water at different positions of wave surface ms

圖41 高速射彈波面不同位置入水，沖擊和流動形成階段彈體載荷與彈道特征Fig.41 Load and ballistic characteristics of projectile entry water at different positions of wave surface in the formation stage of impact and flow

高速射彈波面不同位置入水空泡敞開階段、空泡閉合階段以及空泡潰滅階段，彈體所受載荷與彈體特征詳見4.4節(jié)和4.3節(jié)。

高速射彈波面不同位置入水彈體所受的波浪載荷已包含于其水動力載荷中，由式(12)二階Stokes波的最大橫向速度和垂向速度分別為3.08 m/s和3.06 m/s，遠小于射彈入水速度409 m/s；同時僅分析高速射彈入水前23 ms，遠小于二階Stokes波的周期5.06 s。故高速射彈入水過程中波浪載荷對彈體運動軌跡影響不大，可以忽略。

5 結(jié) 論

本文針對某一型艦載射彈在波峰、波節(jié)1、波谷及波節(jié)2下高速旋轉(zhuǎn)入水過程進行了數(shù)值模擬，分析了不同波面位置下彈體入水過程中的空泡形態(tài)、彈道特征和流體動力載荷特征的變化規(guī)律，得到如下結(jié)論：

(1) 波面不同位置入水，入水初期產(chǎn)生覆蓋彈身的超空泡，空泡不對稱，隨著入水深度的增加，空泡的不對稱性逐漸減??；入水中后期彈體尾部下垂，刺破空泡壁面，空泡形態(tài)產(chǎn)生變化，波節(jié)2工況彈體姿態(tài)最先發(fā)生翻轉(zhuǎn)，波節(jié)1工況彈體下表面沾濕和姿態(tài)變化最慢，相同時間內(nèi)彈體姿態(tài)發(fā)生很大變化依次是波節(jié)2、波峰、波谷和波節(jié)1。

(2) 波面不同位置入水，入水后期彈體質(zhì)心均會產(chǎn)生不同程度的側(cè)向Z向偏移；仿真結(jié)束時刻彈體側(cè)向位移、彈體質(zhì)心速度變化量、姿態(tài)角變化量和角速度變化量由大到小工況依次為波節(jié)2、波峰、波谷和波節(jié)1。

(3) 撞擊水面時，彈體阻力和升力系數(shù)產(chǎn)生瞬間峰值隨后趨于平緩；入水后期由于彈體表面的大量沾濕，阻力和升力系數(shù)迅速增大。入水初期彈體在空泡中航行，彈體力矩系數(shù)趨于平緩且量較小；由于彈體旋轉(zhuǎn)運動的存在彈體偏航力矩和俯仰力矩系數(shù)呈現(xiàn)正負交替變化，仿真結(jié)束時彈體偏航力矩系數(shù)、俯仰力矩系數(shù)的量級從大到小依次是波節(jié)2、波峰、波谷和波節(jié)1。

(4) 彈體不同波面位置入水，實際上改變了彈體入水的有效入水角，雖然四種工況彈體名義入水角均為25°，但是波峰、波節(jié)1、波谷與波節(jié)2入水彈體的有效入水角分別為25.003 2°，39.662 0°，25.009 8°與10.309 4°，因此在相同仿真時間內(nèi)彈體最先發(fā)生失穩(wěn)的工況依次為波節(jié)2、波峰、波谷和波節(jié)1。