陸辰昱，王 鑫，周志壇，梁曉揚(yáng)，樂貴高

(1. 南京理工大學(xué)機(jī)械工程學(xué)院，南京 210094；2. 中國人民解放軍96901部隊，24分隊，北京 100094)

0 引 言

潛射導(dǎo)彈的發(fā)射方式目前主要有兩種：濕式發(fā)射和干式發(fā)射。采用干式發(fā)射時，運(yùn)載器包裹著導(dǎo)彈在水下航行。期間導(dǎo)彈不與水直接接觸，并且大部分水動力載荷由運(yùn)載器承擔(dān)。當(dāng)運(yùn)載器穿過水面時，導(dǎo)彈助推器點(diǎn)火推動其與運(yùn)載器分離，轉(zhuǎn)入空中飛行彈道[1]。彈器分離階段雖然時間短，卻提供了導(dǎo)彈空中飛行的初始條件，是發(fā)射成功的關(guān)鍵[2-3]。然而分離過程在自由液面進(jìn)行，受海浪和多體相互作用等因素影響，分離運(yùn)動呈現(xiàn)很強(qiáng)的非線性。因此，有必要對運(yùn)載器式潛射導(dǎo)彈水面分離運(yùn)動特性加以研究。

目前國外關(guān)于多體出水分離運(yùn)動的研究較少，關(guān)注更多的是單體跨介質(zhì)運(yùn)動[4-5]和空中多體分離問題[6-8]。而國內(nèi)則有較多學(xué)者開展了飛行器出水分離的相關(guān)研究。彭正梁等[9]建立了高速航行體與運(yùn)載器水面分離的線性多體動力學(xué)方程，并基于勢流理論計算了水動力，進(jìn)而對該過程進(jìn)行了二維彈道仿真；顧媛媛等[10]基于切片法建立了航行器水面分離的動力學(xué)模型，并通過Matlab進(jìn)行數(shù)值模擬，研究了分離過程中航行器的運(yùn)動參數(shù)變化；李晶等[11]對筒式導(dǎo)彈出水分離過程進(jìn)行了研究。通過受力分析建立了多體動力學(xué)方程，進(jìn)而計算了分離過程中姿態(tài)和速度等參數(shù)的變化；劉曜等[12]通過構(gòu)建動力學(xué)方程組對潛射導(dǎo)彈的水下彈道和水面分離彈道進(jìn)行了仿真研究，并通過與水洞試驗(yàn)對比驗(yàn)證了計算模型的精確性。以上研究為運(yùn)載器式潛射導(dǎo)彈水面分離過程的預(yù)測提供了一定的參考依據(jù)。然而，目前絕大多數(shù)研究缺乏對液體相的模擬，且在處理水動力問題時均采取如勢流理論等簡化的液體流動模型。實(shí)際情況下，運(yùn)載器出水分離伴隨著海浪的擾動以及波浪的破碎，其受力情況較為復(fù)雜。因此，需要在計算運(yùn)載器與導(dǎo)彈分離運(yùn)動的同時，精確地模擬自由液面的狀態(tài)。

本文建立了自由液面模型和六自由度運(yùn)動方程，并將以上多介質(zhì)耦合模型的計算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)對比。在驗(yàn)證了計算模型的精度和有效性基礎(chǔ)上，對不同海況以及出水角度下運(yùn)載器式潛射導(dǎo)彈水面分離過程中的質(zhì)心位移，姿態(tài)角及軸向速度等進(jìn)行了三維仿真計算。本文的數(shù)值方法和仿真結(jié)果可以為水面分離系統(tǒng)的設(shè)計提供指導(dǎo)。

1 分離過程與坐標(biāo)系建立

如圖1所示，運(yùn)載器在水下發(fā)射出筒后依靠慣性和浮力上升，隨后沖出水面。此時運(yùn)載器蓋中傳感器獲得信號，繼續(xù)上升一段距離后運(yùn)載器頂蓋拋落，隨后導(dǎo)彈助推器點(diǎn)火產(chǎn)生推力，推動導(dǎo)彈向上運(yùn)動。同時，運(yùn)載器受反推力作用落入水中，完成運(yùn)載器式潛射導(dǎo)彈的水面分離[13]。在此過程中如圖2建立固定坐標(biāo)系o0x0y0z0，運(yùn)載器移動坐標(biāo)系o1x1y1z1，導(dǎo)彈移動坐標(biāo)系o2x2y2z2。初始時刻(t=0 s)移動坐標(biāo)系與固定坐標(biāo)系保持平行。其中o0位于初始時刻運(yùn)載器頂端中心，o1位于運(yùn)載器質(zhì)心，o2位于導(dǎo)彈質(zhì)心。o1z1與o2z2分別沿著運(yùn)載器與導(dǎo)彈軸線方向。記運(yùn)載器動坐標(biāo)系及導(dǎo)彈動坐標(biāo)系與固定坐標(biāo)系之間的歐拉角為滾轉(zhuǎn)角φ，俯仰角θ，偏航角ψ。

圖1 運(yùn)載器式潛射導(dǎo)彈水面分離過程示意圖Fig.1 Schematic diagram of separation process of sub-marine launched missile near the free surface

圖2 坐標(biāo)系示意圖Fig.2 Schematic diagram of coordinate systems

2 計算模型

如圖3所示，運(yùn)載器式潛射導(dǎo)彈采用單筒單彈結(jié)構(gòu)。運(yùn)載器直徑為1 m，長度為10 m，材料密度為500 kg/m3；導(dǎo)彈直徑為0.6 m，長度為6 m，材料密度為1000 kg/m3；運(yùn)載器初始出水速度為10 m/s；導(dǎo)彈在運(yùn)載器出水后0.5 s開始點(diǎn)火產(chǎn)生推力，推力大小為100 kN，同時，運(yùn)載器受到同樣大小的反向推力。為了適應(yīng)運(yùn)載器與導(dǎo)彈的姿態(tài)分布并充分模擬波浪的傳遞過程，需要較大的計算域。

圖3 水面分離計算模型Fig.3 Computational model of separation near the free surface

本文為了研究海面波浪以及出水角度對運(yùn)載器式潛射導(dǎo)彈水面分離運(yùn)動特性的影響，共設(shè)立三組計算工況如表1所示。其中工況2與工況3均采用5級海況下的浪高。

表1 計算工況Table 1 Computational conditions

3 數(shù)值方法

3.1 自由液面模型

基于分?jǐn)?shù)容積障礙法(FAVOUR)[14]，建立笛卡爾坐標(biāo)系下的三維不可壓縮Navier-Stokes方程，形式如下：

(1)

(2)

(3)

式中：ρ為流體密度；u，v，w分別為流體速度在x，y，z方向上的分量；Ax，Ay，Az分別為流體在x，y，z方向上的面積分?jǐn)?shù)；VF為流體的體積分?jǐn)?shù)；p為流體壓強(qiáng)；Gx，Gy，Gz為重力加速度；fx，fy，fz為黏滯力加速度。

湍流模型采用RNGk-ε湍流模型[15-16]，與標(biāo)準(zhǔn)k-ε模型相比，它更適于描述有強(qiáng)剪切區(qū)域的流動，因此可以更精確地計算水面波浪破碎過程。其中湍動能k方程為：

(4)

式中：ui，uj為xi，xj方向上的速度分量；μ為流體黏度；μt為動力黏度；Gk為平均速度梯度引起的湍動能k的產(chǎn)生項(xiàng)；常數(shù)系數(shù)σk=1.39。

湍動耗散率ε方程為：

(5)

(6)

對自由液面采用由Hirt等[17]提出的流體體積函數(shù)(VOF)進(jìn)行求解。通過計算網(wǎng)格單元中流體和網(wǎng)格體積比函數(shù)F(x,y,z,t)來確定自由液面邊界。函數(shù)F滿足如下控制方程：

(7)

3.2 六自由度運(yùn)動方程

假設(shè)運(yùn)載器和導(dǎo)彈均為剛體，則其一般運(yùn)動可分為平移與轉(zhuǎn)動。因此，分離運(yùn)動方程形式如下：

(8)

式中：F為剛體所受的總外力，m為剛體質(zhì)量，V為質(zhì)心速度，T為作用在質(zhì)心的總力矩，ω為剛體角速度矢量，J為剛體動坐標(biāo)系下的慣性張量。

若定義B為運(yùn)動剛體表面上的任意一點(diǎn)，則該點(diǎn)速度可由下式計算：

VB=V+ω×rB

(9)

式中：rB為質(zhì)心到B點(diǎn)的矢量。

基于FAVOR法和流體連續(xù)性方程，剛體運(yùn)動的控制方程如下[18]：

(10)

式中：Vc為網(wǎng)格單元的體積；SB，n，VB分別為網(wǎng)格單元內(nèi)運(yùn)動物體的表面積、單位法向量和速度。

對于剛體之間的碰撞運(yùn)動采用Stronge[19]提出的多體動力學(xué)方程進(jìn)行求解。假設(shè)碰撞過程無摩擦，兩碰撞體的質(zhì)量為M和M′，其質(zhì)心分別位于點(diǎn)G和G′，碰撞接觸點(diǎn)分別為C和C′。在接觸點(diǎn)建立碰撞參考系n1n2n3。其中,n1和n2位于公共切面，n3沿著公共切面的法線方向。根據(jù)愛因斯坦求和公式，接觸點(diǎn)C和C′之間的相對速度vi可由下式計算：

(11)

3.3 數(shù)值模型關(guān)系與計算精度

通過上述模型可在每一個時間步內(nèi)完成流體與剛體的計算。各模型之間的數(shù)值傳遞關(guān)系如圖4所示。數(shù)值模擬過程采用有限差分法對控制方程進(jìn)行離散處理，并在時間和空間上保持一階計算精度。

圖4 數(shù)值模型關(guān)系Fig.4 Relationship between numerical models

3.4 模型驗(yàn)證

根據(jù)上述自由液面模型與六自由度運(yùn)動方程，如表2所示，對不同初始傾斜角度的圓柱入水過程進(jìn)行計算，并與實(shí)驗(yàn)結(jié)果[20]對比以驗(yàn)證模型的精確性。圓柱的長徑比為4，密度為900 kg/m3。圓柱初始入水速度垂直向下。受實(shí)驗(yàn)誤差影響，其數(shù)值略有波動。

圖5顯示了初始傾角為55.6°，5個不同時刻下，計算結(jié)果的速度云圖與實(shí)驗(yàn)高速攝影圖的對比情況。從圖5可以看出,兩者各時刻圓柱的傾角基本一致，且圓柱兩端由于高速入水產(chǎn)生的空泡外形也較為相似。圖6顯示了初始傾角為55.6°，圓柱入水過程中質(zhì)心位置和質(zhì)心垂直速度的計算值與實(shí)驗(yàn)測量值的變化曲線。從圖6可以看出，計算結(jié)果的質(zhì)心位置與實(shí)驗(yàn)結(jié)果較為吻合。同時，兩者的質(zhì)心垂直速度也基本一致，僅在末端出現(xiàn)輕微誤差。這主要是受計算網(wǎng)格限制，無法完全模擬空泡拉斷過程中的細(xì)小氣泡而造成。

表2 圓柱入水初始條件Table 2 Water entry conditions of the cylinder

圖5 速度場計算結(jié)果與實(shí)驗(yàn)高速攝影圖像(θ0=55.6°)Fig.5 Calculated results of velocity magnitude and high-speed photographs(θ0=55.6°)

圖6 質(zhì)心位置與質(zhì)心垂直速度對比(θ0=55.6°)Fig.6 Comparison of mass center position and comparison of mass center vertical velocity(θ0=55.6°)

在此基礎(chǔ)上，對其余初始傾角的工況也進(jìn)行了對比計算。如圖7所示，不同初始傾角下圓柱入水，在100 ms時兩端的空泡形狀仍較為相似。同時，從圖8可以看出，計算得到的圓柱入水角位移曲線與實(shí)驗(yàn)值保持較好的吻合。綜上，仿真模型具有較高的精度。

圖7 速度場計算結(jié)果與實(shí)驗(yàn)高速攝影圖像(t=100 ms)Fig.7 Calculated results of velocity magnitude and high-speed photographs (t=100 ms)

圖8 不同初始傾角下圓柱入水角位移對比Fig.8 Comparison of angular displacement for the water entry with different initial inclined angles

4 計算結(jié)果及分析

4.1 分離狀態(tài)

圖9為各工況下y+方向的分離狀態(tài)圖像序列。由圖9可知，從t=0.5 s導(dǎo)彈產(chǎn)生推力開始，工況1～3均在1.03 s完成分離過程。其中，工況1由于海面無波浪且海水橫向流速為0 m/s，運(yùn)載器與導(dǎo)彈在分離后仍能保持相對垂直的姿態(tài)。而工況2在5級海況下進(jìn)行分離，雖然分離時間較短，運(yùn)載器和導(dǎo)彈仍由于波浪沖擊產(chǎn)生較大的傾角。由圖9(c)可知，工況3雖然與工況2的分離環(huán)境相同，但由于運(yùn)載器初始時刻傾斜10°出水，在分離結(jié)束后運(yùn)載器與導(dǎo)彈最終偏轉(zhuǎn)到垂直姿態(tài)。而垂直姿態(tài)更有利于導(dǎo)彈的彈道控制，實(shí)施全方位攻擊。因此，若能控制運(yùn)載器以一定角度出水，則在較高海況等級下導(dǎo)彈仍能在分離后獲得較好的初始狀態(tài)。

4.2 質(zhì)心位置

圖10和圖11顯示了各工況下運(yùn)載器與導(dǎo)彈在分離過程中的質(zhì)心位置。其中，工況2下由于受波浪推動，運(yùn)載器與導(dǎo)彈質(zhì)心均在x方向有較大位移。而由于出水瞬間外部載荷急劇變化，工況1,2,3中運(yùn)載器質(zhì)心均在y方向上有輕微晃動。其中工況3由于運(yùn)載器傾斜出水，受波浪擾動力大，運(yùn)載器與導(dǎo)彈質(zhì)心在y方向位移相對較大。此外，各工況下運(yùn)載器受反推力作用，從0.8 s開始其質(zhì)心在z方向高度不再上升，轉(zhuǎn)而落入水中。且同一時刻下，工況3中運(yùn)載器與導(dǎo)彈質(zhì)心在z方向位置較工況1與工況2略高?？梢酝茢?，此時波浪對運(yùn)載器有一定的向上推動作用。

圖9 工況1(a)，工況2(b)，工況3(c)分離狀態(tài)圖像序列Fig.9 The image sequence of separating state in working condition 1(a), 2(b) and 3(c)

圖10 運(yùn)載器質(zhì)心位置Fig.10 Mass center position of launch canister

圖11 導(dǎo)彈質(zhì)心位置Fig.11 Mass center position of missile

4.3 姿態(tài)角

圖12和圖13為各工況下運(yùn)載器與導(dǎo)彈在分離過程中的姿態(tài)角。從圖12～13可以看出，工況1下，運(yùn)載器與導(dǎo)彈各姿態(tài)角變化數(shù)值較小。而工況2受海面波浪影響，各姿態(tài)角顯著增大。其中運(yùn)載器與導(dǎo)彈的滾轉(zhuǎn)角φ和俯仰角θ從0.2 s受到波浪沖擊力后逐漸增大，而偏航角ψ從0.5 s導(dǎo)彈點(diǎn)火產(chǎn)生推力開始才逐漸增大。工況3下，由于傾斜出水，受波浪干擾力大，運(yùn)載器與導(dǎo)彈的滾轉(zhuǎn)角φ相對較大。而運(yùn)載器與導(dǎo)彈的俯仰角θ從最開始的-10°偏轉(zhuǎn)回0°左右。各工況下姿態(tài)角的最大絕對值如表3所示。其中工況2下運(yùn)載器與導(dǎo)彈俯仰角最大絕對值達(dá)到22°左右。雖然仍在發(fā)射允許的安全范圍內(nèi)，但若浪高進(jìn)一步增加，導(dǎo)彈有傾覆的危險。

圖12 運(yùn)載器姿態(tài)角Fig.12 Attitude angles of launch canister

圖13 導(dǎo)彈姿態(tài)角Fig.13 Attitude angles of missile

表3 各工況下姿態(tài)角最大絕對值Table 3 Maximum absolute value of attitude angles under each working condition

4.4 軸向速度

圖14與圖15為各工況下運(yùn)載器與導(dǎo)彈在分離過程中的軸向速度。如圖14所示，各工況下，0.5 s即導(dǎo)彈點(diǎn)火前，運(yùn)載器軸向速度均保持在10 m/s左右。隨后運(yùn)載器受反推力作用軸向速度降低到-7 m/s左右。整個分離過程中，運(yùn)載器軸向速度由正到負(fù)，最終落入水中。其中工況2由于運(yùn)載器俯仰角相對較大，受波浪的推動，反向速度增加相對較慢。由圖15可知，導(dǎo)彈受運(yùn)載器保護(hù)，各工況下其軸向速度基本一致。因此，導(dǎo)彈軸向速度不受波浪與出水角度的影響，且導(dǎo)彈分離完成后的軸向速度均達(dá)到35 m/s左右，滿足出筒速度的要求。

圖14 運(yùn)載器軸向速度Fig.14 Axial speed of launch canister

圖15 導(dǎo)彈軸向速度Fig.15 Axial speed of missile

5 結(jié) 論

本文建立了自由液面模型與六自由度運(yùn)動方程，并通過與試驗(yàn)數(shù)據(jù)對比驗(yàn)證了模型的精度。在此基礎(chǔ)上對不同海況以及出水角度下運(yùn)載器式潛射導(dǎo)彈水面分離的運(yùn)動特性問題進(jìn)行了三維數(shù)值模擬，主要結(jié)論如下：

1)各工況下運(yùn)載器式潛射導(dǎo)彈均在出水后1.03 s左右完成分離。分離過程中，運(yùn)載器在y方向上存在輕微晃動。各工況下運(yùn)載器與導(dǎo)彈的姿態(tài)角與軸向速度變化較為穩(wěn)定，且均在安全范圍內(nèi)，證明了該種分離方式的優(yōu)越性。

2)波浪除了推動運(yùn)載器與導(dǎo)彈在橫向流速的方向產(chǎn)生一定位移外，其影響主要體現(xiàn)在對姿態(tài)角散布范圍的增大。其中，增加最為明顯的為俯仰角。在5級海況下，運(yùn)載器與導(dǎo)彈俯仰角的最大絕對值達(dá)到22°。

3)控制運(yùn)載器以一定角度出水可以保證導(dǎo)彈分離后的相對垂直姿態(tài)，但由于傾斜出水受波浪擾動力大，運(yùn)載器與導(dǎo)彈的滾轉(zhuǎn)角增大，相應(yīng)的其在y方向上的位移也有所增大。