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

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

0 引 言

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

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

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

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

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

圖1 運(yùn)載器式潛射導(dǎo)彈水面分離過(guò)程示意圖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 計(jì)算模型

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

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

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

表1 計(jì)算工況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ū)域的流動(dòng)，因此可以更精確地計(jì)算水面波浪破碎過(guò)程。其中湍動(dòng)能k方程為：

(4)

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

湍動(dòng)耗散率ε方程為：

(5)

(6)

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

(7)

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

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

(8)

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

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

VB=V+ω×rB

(9)

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

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

(10)

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

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

(11)

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

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

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

3.4 模型驗(yàn)證

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

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

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

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

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

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

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

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

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

4.1 分離狀態(tài)

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

4.2 質(zhì)心位置

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

圖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)彈在分離過(guò)程中的姿態(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)生推力開(kāi)始才逐漸增大。工況3下，由于傾斜出水，受波浪干擾力大，運(yùn)載器與導(dǎo)彈的滾轉(zhuǎn)角φ相對(duì)較大。而運(yùn)載器與導(dǎo)彈的俯仰角θ從最開(kāi)始的-10°偏轉(zhuǎn)回0°左右。各工況下姿態(tài)角的最大絕對(duì)值如表3所示。其中工況2下運(yùn)載器與導(dǎo)彈俯仰角最大絕對(duì)值達(dá)到22°左右。雖然仍在發(fā)射允許的安全范圍內(nèi)，但若浪高進(jìn)一步增加，導(dǎo)彈有傾覆的危險(xiǎn)。

圖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)角最大絕對(duì)值Table 3 Maximum absolute value of attitude angles under each working condition

4.4 軸向速度

圖14與圖15為各工況下運(yùn)載器與導(dǎo)彈在分離過(guò)程中的軸向速度。如圖14所示，各工況下，0.5 s即導(dǎo)彈點(diǎn)火前，運(yùn)載器軸向速度均保持在10 m/s左右。隨后運(yùn)載器受反推力作用軸向速度降低到-7 m/s左右。整個(gè)分離過(guò)程中，運(yùn)載器軸向速度由正到負(fù)，最終落入水中。其中工況2由于運(yùn)載器俯仰角相對(duì)較大，受波浪的推動(dòng)，反向速度增加相對(duì)較慢。由圖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)動(dòng)方程，并通過(guò)與試驗(yàn)數(shù)據(jù)對(duì)比驗(yàn)證了模型的精度。在此基礎(chǔ)上對(duì)不同海況以及出水角度下運(yùn)載器式潛射導(dǎo)彈水面分離的運(yùn)動(dòng)特性問(wèn)題進(jìn)行了三維數(shù)值模擬，主要結(jié)論如下：

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

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

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