(北京宇航系統(tǒng)工程研究所，北京 100076)

0 引言

潛射航行體的發(fā)射是一個涉及多因素干擾、多相介質(zhì)參與的復(fù)雜物理過程，目前為應(yīng)對各種復(fù)雜環(huán)境需求，航行體朝著變深度發(fā)射、簡化結(jié)構(gòu)、提高發(fā)射可靠性等趨勢發(fā)展[1]。其中發(fā)射可靠性一直是其重點(diǎn)研究方向之一，由于航行體需經(jīng)歷穿越水層、空中點(diǎn)火等對其結(jié)構(gòu)和強(qiáng)度有極高要求的過程，且以海洋為發(fā)射環(huán)境使得諸多不可預(yù)知因素成為抑制航行體彈道穩(wěn)定性的技術(shù)瓶頸，又考慮水下發(fā)射點(diǎn)火升空時間極短，需要完成多個復(fù)雜的控制動作及過程，因此如何能在可預(yù)知和不可預(yù)知的多因素作用下開展航行體發(fā)射可靠性、提高其彈道穩(wěn)定性的控制技術(shù)研究顯得極為必要。

水下平臺垂直冷發(fā)射[2]的航行體水下運(yùn)動過程可分為三個階段：出筒段、水中航行段和出水段[3]。發(fā)射時平臺須保持一定的航行速度來維持其操縱性，平臺航速的存在破壞了流場的對稱性，使得航行體所受流體動力呈現(xiàn)顯著的三維特征，進(jìn)而出現(xiàn)較大的俯仰姿態(tài)[4-5]。變深度發(fā)射有利于航行體發(fā)射的靈活性，將是水下發(fā)射技術(shù)的發(fā)展趨勢之一，因而探討有平臺航速條件下航行體的流體動力特征，研究平臺航速和發(fā)射水深對航行體流體動力和運(yùn)動速度的影響具有重要的工程價值。

隨著多相流模型、動網(wǎng)格技術(shù)、湍流模型等CFD技術(shù)的不斷發(fā)展，目前數(shù)值模擬已成為研究水下發(fā)射技術(shù)的常用手段。多相流模型根據(jù)界面處理的不同大致可分為兩類，一類是自由界面追蹤模型，如VOF[6]、MAC[7]模型，可以處理自由界面變化較劇烈但界面基本清晰的流動問題；另一類是混合模型，如Mixture模型等，可模擬各相之間強(qiáng)烈摻混、界面基本不清晰的流動問題。Dyment[8]，劉志勇等[9]分別采用VOF方法和改進(jìn)的MAC方法模擬了尾空泡發(fā)展過程；李杰等[10]采用VOF方法對細(xì)長回轉(zhuǎn)體出水過程進(jìn)行了數(shù)值模擬；張紅軍等[11]采用Mixture模型對航行體出筒過程進(jìn)行了三維非定常數(shù)值模擬研究；楊曉光[12]等采用VOF模型和動網(wǎng)格技術(shù)對航行體水下運(yùn)動及出水過程進(jìn)行了三維流場仿真。Ma 等采用MIX模型對出口航行體排氣氣膜孔參數(shù)進(jìn)行了三維仿真[13]，后續(xù)又進(jìn)行了考慮平臺速度和偏航角不確定性的水下航行器彈道和姿態(tài)的數(shù)值研究[14]。

本文研究將采用Mixture模型、標(biāo)準(zhǔn)k-ε湍流模型以及動網(wǎng)格技術(shù)求解水下航行體非定常發(fā)射過程，通過對流場計(jì)算獲得航行體所受水動力，而后將水動力作為力學(xué)參數(shù)輸入以求解航行體的運(yùn)動過程，從而實(shí)現(xiàn)航行體水下發(fā)射多相流場和彈道耦合求解。通過分析航行體表面壓力以及流體動力的三維特性，以探討發(fā)射平臺航速和發(fā)射水深對航行體所受流體動力和運(yùn)動速度的影響。

1 坐標(biāo)系及無量綱化定義

圖1所示為發(fā)射系和航行體體系坐標(biāo)系示意圖，本文將在發(fā)射系中討論航行體的運(yùn)動速度特征，在航行體體系中討論流體動力特征。坐標(biāo)系的定義如下:

發(fā)射系：以航行體初始時刻質(zhì)心位置為原點(diǎn)O，以發(fā)射平臺航速方向?yàn)閄軸(橫向)，以垂直出筒方向?yàn)閅軸(縱向)，Z軸方向由右手法則確定。航行體體系：以航行體實(shí)時的質(zhì)心位置為原點(diǎn)O1，法向方向?yàn)閄1軸，體軸方向?yàn)閅1軸，Z1軸方向由右手法則確定。

本文中壓力p和時間t分別以發(fā)射水深靜壓P和出筒時刻T進(jìn)行無量綱化；各速度分量均以水深H、平臺航速U下計(jì)算獲得的最大值為參考進(jìn)行無量綱化；對力和力矩的無量綱化定義為

Cd=F/(0.5ρV2S)

(1)

Cm=M/(0.5ρV2SL)

(2)

式中ρ——水的密度;

V——選取水深H、平臺航速U下的航行體最大縱向速度;

S——航行體特征橫截面積;

L——航行體軸向長度。

2 數(shù)值計(jì)算模型

本文采用動網(wǎng)格技術(shù)將多相流流場和航行體彈道進(jìn)行耦合計(jì)算。首先利用多相流模型計(jì)算當(dāng)前時刻瞬態(tài)流場，通過對航行體表面壓力積分，獲得當(dāng)前時刻航行體所受流體力和力矩，在此基礎(chǔ)上計(jì)算航行體的加速度、速度和位移等；然后利用動網(wǎng)格技術(shù)獲得下一時刻的航行體邊界和網(wǎng)格，進(jìn)而求解下一時刻的瞬態(tài)流場，接著計(jì)算下一時刻航行體的加速度、速度和位移等。上述過程不斷循環(huán)，直至航行體出水為止。

2.1 流體控制方程

本文采用商用數(shù)值計(jì)算軟件Ansys Fluent,采用Mixture模型模擬氣液兩相流動，忽略相間的滑移速度，計(jì)算中考慮水體為不可壓縮相，氣體為可壓縮相。采用有限體積法離散控制方程，結(jié)合Simple算法來求解下列方程。

連續(xù)性方程

(3)

動量方程

(4)

能量方程

(5)

次相體積分?jǐn)?shù)方程

(6)

狀態(tài)方程

pg=ρgRT

(7)

ρw=const

(8)

αk——第k相的體積分?jǐn)?shù)；

keff——有效傳熱系數(shù)；

SE——包含了所有的體積熱源。

上述方程采用標(biāo)準(zhǔn)k-ε湍流模型封閉。

2.2 航行體動力學(xué)方程

本文研究中將航行體當(dāng)作剛體，尾出筒后開始計(jì)算其俯仰運(yùn)動(忽略出筒過程的俯仰運(yùn)動)。通過下式計(jì)算航行體質(zhì)心的合力和繞質(zhì)心的合力矩

(9)

(10)

2.3 計(jì)算域、邊界條件、網(wǎng)格劃分及更新方法

本文建立如圖2所示的三維對稱模型，計(jì)算域包括發(fā)射筒、水域和大氣域。發(fā)射筒底燃?xì)馊肟跒閴毫θ肟谶吔纾叫畜w尾出筒后改為壁面邊界；水域和空氣域外圍均為壓力出口邊界。為計(jì)算尾出筒后航行體的俯仰姿態(tài)，在航行體運(yùn)動區(qū)域內(nèi)劃分非結(jié)構(gòu)網(wǎng)格，采用網(wǎng)格當(dāng)?shù)刂貥?gòu)更新動網(wǎng)格，為保證航行體壁面附近網(wǎng)格質(zhì)量，對航行體表面附近進(jìn)行了網(wǎng)格加密，該區(qū)域內(nèi)的網(wǎng)格質(zhì)量在航行體運(yùn)動過程不發(fā)生變化。數(shù)值計(jì)算網(wǎng)格采用Gambit生成，初始計(jì)算網(wǎng)格數(shù)量為600萬。

3 航行體表面壓力和流體動力三維特征分析

為獲得航行體表面壓力的三維特征，計(jì)算時在航行體迎背水面布置了多個壓力監(jiān)測點(diǎn)，測點(diǎn)1為航行體駐點(diǎn)，測點(diǎn)2～4分布于頭錐段，測點(diǎn)5～8分布在柱段，測點(diǎn)9位于尾段。各測點(diǎn)距離航行體實(shí)際頂點(diǎn)的長度如表1所示(以航行體軸向長度L為參考)。

圖3對比了測點(diǎn)1和迎水面測點(diǎn)3、5、7、9處的計(jì)算與試驗(yàn)壓力值。從曲線分布可知，計(jì)算值與試驗(yàn)值吻合較好，驗(yàn)證了計(jì)算模型的合理性。航行體出筒及水中航行過程中，駐點(diǎn)壓力隨著航行體速度增加而升高、下降而遞減；頭錐段壓力因?yàn)槔@流作用和環(huán)境壓力的下降而呈現(xiàn)出下降趨勢；柱段脫離發(fā)射筒后，壓力隨著環(huán)境水壓的降低而降低，柱段靠近尾段部分(如測點(diǎn)7)在尾出筒時受燃?xì)饣厣渥饔眯纬闪嗣黠@的壓力脈沖，且尾出筒后受尾空泡生成演化的影響，呈現(xiàn)出明顯的壓力振蕩特性。

表1壓力監(jiān)控點(diǎn)布置

測點(diǎn)序號123456789測點(diǎn)位置00.034 L0.099 L0.198 L0.297 L0.481 L0.741 L0.916 L0.962 L

如圖4所示，發(fā)射平臺航速導(dǎo)致航行體迎背水面流場不對稱，航行體以一定的俯仰姿態(tài)出水。圖5給出了測點(diǎn)2～9迎背水面壓差時變曲線，從中可以看出，迎水面壓力高于背水面壓力，壓差主要集中在航行體頭錐段。同時頭錐段壓差在出水前隨著航行體速度增大而升高、下降而遞降；柱段壓差主要出現(xiàn)在水中航行段，兩級上明顯小于頭錐段最大壓差，且受燃?xì)夂笮в绊懀拷捕尾糠?如測點(diǎn)8)壓差在出筒時刻附近形成了脈沖變化；尾段壓差則主要受尾空泡演化影響，除在尾空泡拉斷后一段時間形成了一定量值的壓差外，大部分時間幾乎為零。

迎背水面壓差的存在導(dǎo)致航行體承受較大的法向力和俯仰力矩。從圖7給出了航行體所受流體法向力和俯仰力矩的分布可知，整個發(fā)射過程航行體所受法向力、俯仰力矩主要由頭錐段和柱段組成，尾段貢獻(xiàn)幾乎為零：頭錐段法向力量值在出筒過程中不斷增大，并成為法向合力的主要組成部分，在尾出筒后逐漸減?。浑S著柱段進(jìn)入水中，柱段法向力量值逐漸增大，水中航行段繼續(xù)保持增大趨勢，柱段出水過程逐漸減小至零；出水前，頭錐段力矩在出筒前后呈現(xiàn)先增大后減小的變化趨勢，并且決定了此階段俯仰合力矩的變化趨勢；柱段力矩在出水前量值較小，隨著柱段出水逐漸減小至一個較大的負(fù)值后逐漸增大至零，本文中稱之為俯仰力矩的“反號”變化，并決定了此階段俯仰合力矩的變化趨勢，這是由于法向合力作用點(diǎn)移動到質(zhì)心之后所致。

4 發(fā)射平臺航速影響分析

本文在相同發(fā)射水深H條件下，研究了發(fā)射平臺航速(0.4U、0.7U和1.0U)對航行體流體動力和運(yùn)動速度的影響。圖7所示為，不同航速下航行體所受流體動力的時變曲線。由圖可知，流體軸向力與尾段壓力的變化規(guī)律基本一致。如圖3中所示測點(diǎn)9壓力在尾出筒后呈現(xiàn)周期性振蕩變化；隨著發(fā)射平臺航速增加，流體法向力和俯仰力矩量值增大，但軸向力幾乎不變。

圖8比較了不同平臺航速下航行體運(yùn)動速度的時變曲線。航行體尾出筒后，橫向速度逐漸減小，且出現(xiàn)了反號變化(即沿相反的方向運(yùn)動)，尾出水后不再變化；縱向速度逐漸減?。桓┭鼋撬俣戎饾u增大，受俯仰力矩“反號”變化的影響，頭出水后逐漸減小，尾出水后不再變化；隨著發(fā)射平臺航速增大，橫向速度下降加快，俯仰角速度增大，但縱向速度幾乎不受影響。圖9所示為尾出水后橫向速度和俯仰角速度與初始運(yùn)動時刻的絕對變化量，該變化量與平臺航速基本呈線性關(guān)系。

5 發(fā)射水深影響分析

下文在相同發(fā)射平臺航速U條件下，探討了不同發(fā)射水深(0.8H、1.0H和1.2H)對航行體流體動力和運(yùn)動速度的影響。圖10所示為不同發(fā)射水深下航行體所受流體動力的時變曲線。由圖可知不同發(fā)射水深條件下，航行體所受外力變化趨勢基本一致；隨發(fā)射水深增加，航行體水中段運(yùn)動時間增長，軸向力經(jīng)歷周期性振蕩次數(shù)增多，導(dǎo)致出水引起的法向力和俯仰力矩變化發(fā)生時間延遲。

由圖11所示不同發(fā)射水深下航行體的運(yùn)動速度時變曲線可知，不同發(fā)射水深條件下，航行體運(yùn)動速度的變化歷程幾乎一致；隨著發(fā)射水深增加，出水后橫向速度和俯仰角速度量值增加，而縱向速度變化較小。

6 結(jié)論

本文采用數(shù)值模擬的方法對航行體水下垂直發(fā)射多相流場和彈道進(jìn)行耦合求解，獲得了航行體表面壓力和流體動力三維特征，研究了發(fā)射平臺航速和發(fā)射水深對航行體所受流體動力和運(yùn)動速度的影響，主要結(jié)論如下：

(1)基于Mixture模型和動網(wǎng)格技術(shù)實(shí)現(xiàn)了三維非定常多相流場和航行體彈道的耦合計(jì)算，測點(diǎn)壓力計(jì)算結(jié)果與試驗(yàn)結(jié)果吻合良好，驗(yàn)證了本文所建計(jì)算模型的正確性。

(2)發(fā)射平臺航速的存在導(dǎo)致航行體迎水面壓力高于背水面，且主要體現(xiàn)在頭錐段和柱段，進(jìn)而使得航行體承受較大的流體法向力和俯仰力矩，法向力變化趨勢在尾出筒前后分別由頭錐段和柱段決定，俯仰力矩變化趨勢在頭出水前后分別由頭錐段和柱段決定。

(3)相同發(fā)射水深條件下，隨著發(fā)射平臺航速增加，航行體所受流體法向力、俯仰力矩增大，但軸向力幾乎不變；航行體橫向速度下降加快，俯仰角速度增大，且橫向速度和俯仰角速度的變化量與平臺航速成線性正比關(guān)系，但縱向速度變化較小。

(4)相同發(fā)射平臺航速、不同發(fā)射水深條件下，航行體所受流體動力和運(yùn)動速度的變化趨勢基本保持一致；隨著發(fā)射水深增加，航行體出水時間延遲，導(dǎo)致出水后橫向速度、俯仰角速度量值增大，但縱向速度變化較小。