劉海軍，王 聰，鄒振祝，王本利，王柏秋

(哈爾濱工業(yè)大學(xué)航天學(xué)院，150001哈爾濱)

水下航行體的發(fā)射方式可分為有動(dòng)力發(fā)射和無動(dòng)力發(fā)射.根據(jù)布置方式可分為垂直發(fā)射、水平發(fā)射和傾斜發(fā)射;根據(jù)發(fā)射筒內(nèi)航行體周圍的環(huán)境分為干式發(fā)射和濕式發(fā)射;按其發(fā)動(dòng)機(jī)的點(diǎn)火時(shí)機(jī)分為筒內(nèi)點(diǎn)火、水下點(diǎn)火和水面點(diǎn)火等［1］.無論采用哪種水下發(fā)射方式，航行體都要經(jīng)歷出筒、水中航行、出水和空中飛行4個(gè)階段，盡管前3個(gè)階段的運(yùn)動(dòng)距離和時(shí)間比較短，但會(huì)形成復(fù)雜多相流場(chǎng)，對(duì)航行體流體動(dòng)力、彈道及載荷等產(chǎn)生直接影響，從而影響航行體運(yùn)動(dòng)的穩(wěn)定性［2］.

潛射航行體在水下運(yùn)動(dòng)的時(shí)間很短，卻存在復(fù)雜的流場(chǎng)情況，殷崇一等［3］運(yùn)用有限體積法和SIMPLE算法計(jì)算了航行體水下發(fā)射內(nèi)流場(chǎng)的數(shù)值研究.陳雄州等［4］根據(jù)固體發(fā)動(dòng)機(jī)推力最大原則，并結(jié)合水下航行體外彈道運(yùn)動(dòng)方程，求出了不同使用環(huán)境下的發(fā)動(dòng)機(jī)最佳出口壓強(qiáng).盧立秀等［5］運(yùn)用MATLAB對(duì)飛機(jī)內(nèi)彈道彈射機(jī)構(gòu)進(jìn)行了研究.王漢平等［6］采用多相流模型和動(dòng)網(wǎng)格技術(shù)研究了均壓氣體對(duì)航行體射后效的影響.仲維國(guó)等［7］基于魚雷的航行力學(xué)導(dǎo)出水下航行體的彈道方程.Acosta等［8］采用數(shù)值模擬的方法對(duì)重力場(chǎng)中的軸對(duì)稱細(xì)長(zhǎng)體沿軸線的定?？张蓍L(zhǎng)度、空泡截面面積和阻力系數(shù)進(jìn)行了研究，得到了重力場(chǎng)中弗洛德數(shù)對(duì)空泡的長(zhǎng)度和面積影響較大，對(duì)軸對(duì)稱細(xì)長(zhǎng)體的阻力系數(shù)影響較小.Basharova等［9］采用基于細(xì)長(zhǎng)體空泡理論的近似方程組計(jì)算了有重力影響無限流場(chǎng)中沿軸線方向上空泡的形態(tài)和基本尺寸.J.J.Yagla等［10］對(duì)水下垂直發(fā)射系統(tǒng)的動(dòng)力環(huán)境進(jìn)行了研究，提出了筒內(nèi)發(fā)動(dòng)機(jī)點(diǎn)火的方式可以為航行體的出水提供1個(gè)從出筒到出水整個(gè)過程的氣幕，并進(jìn)行了相關(guān)的實(shí)驗(yàn).Chris J.Weiland等［11］研究了不同水深和橫向流對(duì)氣幕的影響.劉筠喬等［12］研究了航行體垂直發(fā)射出筒過程的通氣空泡流.張紅軍等［13］對(duì)不帶空化模型的潛射航行體垂直出筒過程進(jìn)行了研究.劉樂華等［14］采用TVD有限體積法研究了深水發(fā)射的內(nèi)流場(chǎng)，得到燃?xì)饬鲊姵鲋罂焖倥蛎?，進(jìn)而形成高溫燃?xì)馍淞?，?duì)反射板有極大的破壞作用.朱明駿等［15］對(duì)潛射航行體的充氣均壓系統(tǒng)進(jìn)行了研究，并利用PID控制方法對(duì)充壓系統(tǒng)進(jìn)行了設(shè)計(jì).

綜上所述，國(guó)內(nèi)外只是對(duì)潛射航行體的垂直發(fā)射出筒過程的內(nèi)流場(chǎng)進(jìn)行了初步研究，沒有對(duì)航行體出筒過程中的流體動(dòng)力進(jìn)行研究，由于航行體在出筒過程中受到的流體動(dòng)力對(duì)航行體的出筒及出筒后飛行的彈道都有影響，故對(duì)航行體在出筒過程中的流體動(dòng)力特性進(jìn)行研究是必要的.本文采用多相Mixture模型并結(jié)合動(dòng)網(wǎng)格技術(shù)對(duì)運(yùn)動(dòng)流場(chǎng)中混合介質(zhì)RANS方程進(jìn)行求解，實(shí)現(xiàn)了重力場(chǎng)中運(yùn)動(dòng)固體邊界與氣(汽)水流場(chǎng)的耦合求解，針對(duì)采用熱發(fā)射方式的航行體垂直發(fā)射出筒過程的多相流場(chǎng)進(jìn)行了數(shù)值模擬，研究了重力影響下的水下航行體垂直出筒過程和流體動(dòng)力特性.

1 數(shù)學(xué)模型

1.1 控制方程

描述水下航行體垂直出筒過程的控制方程主要包括混合物的連續(xù)方程、動(dòng)量方程、能量方程、湍流模型和空化模型.

混合物連續(xù)性方程為

式中:ui為i方向的速度分量;ρm為混合物的密度，ρm=αlρl+αvρv+αngρng，ρi和αi(i=l，v，ng)分別為水、水蒸氣及不可結(jié)氣體的密度和體積分?jǐn)?shù);fi=αiρi/ρm(i=l，v，ng)分別為液體、水蒸汽及不可凝結(jié)氣體的質(zhì)量分?jǐn)?shù);不可凝結(jié)氣體為空氣;m-是水的汽化源相，m+是水蒸汽液化過程的源相.

混合物動(dòng)量方程為

式中:p是混合介質(zhì)的壓力;μm是運(yùn)動(dòng)粘性系數(shù)，其中n是構(gòu)成混合介質(zhì)的相數(shù)，μk是每個(gè)相的湍流粘性系數(shù);ρm為混合物的密度;g是重力加速度.

1.2 湍流模型

標(biāo)準(zhǔn)的k-ε湍流模型為

湍動(dòng)能k和湍動(dòng)耗散率ε的輸運(yùn)方程如下:

式中:k是混合介質(zhì)的湍動(dòng)能;ε是混合介質(zhì)湍流耗散率;Gb是浮力產(chǎn)生湍動(dòng)能項(xiàng);Gk是速度梯度產(chǎn)生的湍動(dòng)能項(xiàng);μt為湍流粘性系數(shù);ρm為混合物的密度;σk和σε分別為k和ε的普朗特?cái)?shù);可以調(diào)整的經(jīng)驗(yàn)常數(shù)取值分別為:Cε1=1.44，Cε2=1.92，σk=1.0，σε=1.3，Cμ=0.09.

1.3 空化模型

基于均質(zhì)平衡流理論的Singhal空化模型，表達(dá)式如下:

式中:Vch是液相和汽相之間轉(zhuǎn)化的特征速度;τ是液體表面張力;pv是飽和蒸汽壓;fl=1-fv-fng，其中不可凝結(jié)氣體分?jǐn)?shù)fng=1.5×10-5;Cc=0.01和Ce=0.02是兩個(gè)經(jīng)驗(yàn)常數(shù).

1.4 能量方程

能量方程為

式中:keff為傳導(dǎo)率，keff=(k+kt);kt為湍流的熱傳導(dǎo)率為物質(zhì)j的擴(kuò)散通量.上式中等號(hào)右邊的最后3項(xiàng)分別代表由熱傳導(dǎo)、物質(zhì)擴(kuò)散和粘性耗散引起的能量交換.Sh包含化學(xué)反應(yīng)熱在內(nèi)的一切熱源.

1.5 航行體的運(yùn)動(dòng)學(xué)方程

航行體的運(yùn)動(dòng)學(xué)方程由牛頓第二定律得到，對(duì)航行體運(yùn)動(dòng)學(xué)方程使用下式進(jìn)行離散:

式中:Vi是航行體在第i個(gè)時(shí)間步的速度;Vi-1是航行體在第i-1個(gè)時(shí)間步的速度;Δt是時(shí)間步長(zhǎng);Fi-1是航行體在第i-1個(gè)時(shí)間步合外力，m是航行體的質(zhì)量.合力由下式算得:

式中:FG為航行體受到的重力;FB為航行體的底部受到的推力;FG和FB由初始條件給出;FP是航行體前后壓強(qiáng)差產(chǎn)生阻力;FV是航行在粘性流體運(yùn)動(dòng)受到的粘性力阻力，F(xiàn)P和FV均是通過對(duì)流場(chǎng)的計(jì)算得到.

2 數(shù)值模擬

2.1 數(shù)值求解及動(dòng)網(wǎng)格

使用SIMPLEC方法進(jìn)行計(jì)算，采用PRESTO對(duì)控制方程的壓力項(xiàng)進(jìn)行離散，采用二階迎風(fēng)格式對(duì)控制方程中的能量方程和湍流耗散項(xiàng)及對(duì)流項(xiàng)進(jìn)行離散.利用動(dòng)網(wǎng)格技術(shù)解決運(yùn)動(dòng)邊界和計(jì)算域變化的問題.編制控制航行體運(yùn)動(dòng)的程序，通過自定義函數(shù)(UDF)嵌入FLUENT，求解出航行的速度和位移.

采用動(dòng)網(wǎng)格技術(shù)進(jìn)行計(jì)算時(shí)，動(dòng)網(wǎng)格區(qū)域中的網(wǎng)格具有一定對(duì)流運(yùn)動(dòng)，在實(shí)際計(jì)算過程中要去掉動(dòng)網(wǎng)格對(duì)流運(yùn)動(dòng)的影響，動(dòng)網(wǎng)格守恒方程統(tǒng)一形式如下:

式中:ρm為混合物的密度;u為流體的速度矢量;ug為網(wǎng)格的運(yùn)動(dòng)速度;Γ為擴(kuò)散系數(shù);Sφ為標(biāo)量φ的源相;?V為控制體積V的邊界.

2.2 數(shù)值求解方法及正確性

本文采用間接證明的方法，對(duì)比分析數(shù)值模擬與實(shí)驗(yàn)結(jié)果，得到數(shù)值求解方法的正確性.驗(yàn)證算例采用模型為:半球頭型航行體直徑為50 mm、總長(zhǎng)度為200 mm.計(jì)算條件是:未受擾動(dòng)水流場(chǎng)壓強(qiáng)p∞=378 955 Pa，采用非定常過程計(jì)算.以航行體直徑D計(jì)算雷諾數(shù)為Re=106，空化數(shù)σ=0.31，圖1給出了算例中的航行體在尾部離開筒口時(shí)刻的物面壓力系數(shù)分布曲線與文獻(xiàn)［16］中實(shí)驗(yàn)數(shù)據(jù)對(duì)比曲線圖.H為距離航行體的頭部的軸向距離，從圖中看出仿真結(jié)果和實(shí)驗(yàn)數(shù)據(jù)吻合較好.

圖1 實(shí)驗(yàn)結(jié)果與仿真結(jié)果的比較

在計(jì)算過程中采用網(wǎng)格重構(gòu)法，結(jié)合多相流Mixture模型和湍流模型，對(duì)流場(chǎng)混合介質(zhì)RANS方程進(jìn)行求解，進(jìn)而求解固體運(yùn)動(dòng)邊界與汽水流場(chǎng)的流-固耦合問題，采用此方法對(duì)水下航行體垂直發(fā)射出筒過程流體動(dòng)力特性展開數(shù)值模擬研究.

3 計(jì)算模型及邊界條件

航行體的長(zhǎng)細(xì)比L/D=10，L為航行體的長(zhǎng)度，D為航行體的直徑，頭為半球型，發(fā)射深度為50 m.計(jì)算域由筒內(nèi)的氣體和筒外的水域組成.圖2給出了流場(chǎng)計(jì)算區(qū)域、網(wǎng)格和邊界條件的示意圖，在航行體附近的網(wǎng)格進(jìn)行了加密處理.圖3給出了靜網(wǎng)格和動(dòng)網(wǎng)格示意圖，整個(gè)計(jì)算域分為兩塊:中間為動(dòng)態(tài)區(qū)域，網(wǎng)格隨時(shí)間變化，周圍為靜態(tài)區(qū)域，網(wǎng)格固定不動(dòng).流場(chǎng)的長(zhǎng)度是航行體長(zhǎng)度的5倍，寬度是航行體直徑的12倍.

航行體噴管入口是本模型的壓力入口，噴管入口的壓力采用圖4給出隨時(shí)間變化燃?xì)饪倝呵€.

圖2 計(jì)算區(qū)域與網(wǎng)格

圖3 計(jì)算區(qū)域與模型

圖4 燃?xì)饪倝弘S時(shí)間變化曲線

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

4.1 流場(chǎng)分析

圖5給出了水下航行體不同時(shí)刻的流體相體積分?jǐn)?shù)分布云圖.I至III是航行體頭部沒有穿透氣團(tuán)之前的氣團(tuán)的生成過程，此時(shí)氣團(tuán)是半球型.III至VI是航行體頭部穿透氣體時(shí)氣團(tuán)的演化過程，此時(shí)氣團(tuán)演變成橢球型，因航行體的速度逐漸增大和環(huán)境的靜壓力逐漸減小，空化數(shù)逐漸減小，故在航行體的肩部產(chǎn)生空化，并逐漸發(fā)展成包裹航行體肩部的空泡.航行體離開發(fā)射筒后，筒內(nèi)的壓力減小，但是仍高于周圍環(huán)境的壓力，筒口的氣團(tuán)在海水的作用下開始徑向回縮.

圖5 不同時(shí)刻流場(chǎng)體積分?jǐn)?shù)ai云圖

4.2 流場(chǎng)壓力分布

圖6給出了不同時(shí)刻流場(chǎng)的壓力云圖，所選的時(shí)間點(diǎn)與圖5一一對(duì)應(yīng)，從壓力云圖中可以看出，水域流場(chǎng)的壓力分布主要受重力的影響.初始時(shí)刻，筒口附近水域壓力比較大，發(fā)動(dòng)機(jī)點(diǎn)火筒內(nèi)的壓力開始增加.隨著出筒高度的增加周圍環(huán)境壓力開始減小.從IV時(shí)刻后，在航行體的肩部出現(xiàn)了繞流引起的低壓區(qū)，隨出筒高度的增加，航行體的速度增大，在肩部繞流引起的低壓區(qū)繼續(xù)增大，這與圖5的壓力云圖中的肩部空泡的產(chǎn)生是對(duì)應(yīng)的.航行體尾部離開發(fā)射筒后，筒口附近的海水壓力沿徑向往筒口逐漸增加，但筒口氣團(tuán)的壓力仍高于環(huán)境壓力，此時(shí)筒口周圍的水暫時(shí)沒有向筒內(nèi)回灌.

圖6 不同時(shí)刻流場(chǎng)壓力云圖

4.3 航行體表面壓力分析

圖7給出了不同時(shí)刻航行體表面的環(huán)向壓力沿軸向的分布，橫坐標(biāo)代表的是航行體的特征長(zhǎng)度，縱坐標(biāo)代表的是表面壓力與最大壓力的比值，h是距離航行體頭部的高度，L是航行體的長(zhǎng)度，T=L/V0，V0是出筒的速度，T是特征時(shí)間.左側(cè)是航行體的頭部，右側(cè)是航行體的尾部.從圖中看出，隨著出筒時(shí)間的增加，航行體表面的壓力系數(shù)是逐漸降低的，在航行體的頭部出現(xiàn)了由繞流引起的低壓區(qū)，在航行體的尾部由于處在發(fā)射筒內(nèi)高壓燃?xì)庵校霈F(xiàn)了高壓區(qū)，隨出筒時(shí)間的增加航行體尾部壓力逐漸減小并趨緩，這與圖6中壓力場(chǎng)的分布一致.

圖7 不同時(shí)刻航行體表面壓力比值沿軸向分布曲線

4.4 航行體出筒彈道分析

圖8 航行體速度和位移變化曲線

圖9是加速度變化曲線.加速度曲線在初始刻有上下波動(dòng)，這是由于發(fā)動(dòng)機(jī)的推力在逐漸增加至峰值時(shí)，加速度也迅速增大，然而在筒口，由于氣團(tuán)的膨脹和收縮導(dǎo)致加速度又有所波動(dòng)，最后加速度抖動(dòng)下降.

4.5 航行體出筒過程受力分析

圖10給出了壓差力系數(shù)、粘性力系數(shù)和總流體動(dòng)力系數(shù)變化曲線.由圖9可得出，在t＜0.375T時(shí)壓差力系數(shù)、粘性力系數(shù)和總流體動(dòng)力系數(shù)的都在對(duì)應(yīng)的時(shí)間上有波動(dòng).壓差系數(shù)和總流體動(dòng)力系數(shù)波動(dòng)的趨勢(shì)一樣，產(chǎn)生的原因是:1)發(fā)動(dòng)機(jī)的推力隨著時(shí)間的增加逐漸增大，并且到達(dá)峰值;2)氣團(tuán)在沒有被航行體穿透時(shí)，氣團(tuán)在水靜壓的作用下收縮和膨脹，產(chǎn)生瞬間壓力的變化.t＞0.375T時(shí)，壓差系數(shù)曲線逐漸下降且趨緩.粘性力系數(shù)在t＜0.375T時(shí)的波動(dòng)，是由于海水靜壓作用下氣團(tuán)的收縮和膨脹影響高壓燃?xì)獾乃俣?，此時(shí)氣體和航行體相對(duì)速度發(fā)生改變，進(jìn)而影響粘性力系數(shù)，粘性力系數(shù)為正值.航行體受到的粘性力由高速氣體對(duì)航行體的粘性力和水對(duì)航行體的粘性力組成，兩者的作用力方向相反，水的粘性力隨著沾濕面積的增加而增大，因此總粘性力系數(shù)逐漸減小.這說明了壓差系數(shù)是影響總流體動(dòng)力系數(shù)的主要因素，粘性系數(shù)對(duì)總系數(shù)的影響是次要因素.

圖9 航行體加速度變化曲線

圖10 總系數(shù)變化曲線

5 結(jié)論

采用數(shù)值模擬的方法研究了有重力影響的水下航行體垂直發(fā)射出筒過程流體動(dòng)力特性，得到如下結(jié)論:

1)航行體頭部距離筒口0.01D處，氣團(tuán)的膨脹和收縮變化劇烈，加速度出現(xiàn)劇烈波動(dòng)，距離筒口高度的增加，筒口氣團(tuán)的膨脹和收縮變化趨緩;

2)氣團(tuán)的形狀從半球型變化到橢球型，最后變成細(xì)長(zhǎng)型;航行體頭部在海水靜壓作用下出現(xiàn)了高壓和在肩部出現(xiàn)了繞流引起的低壓區(qū);

3)航行體頭部穿透氣團(tuán)之前，氣團(tuán)的膨脹和收縮引起航行體的流體動(dòng)力系數(shù)較大的波動(dòng)，縮短航行體穿透氣團(tuán)的時(shí)間可以減弱航行體載荷的波動(dòng)，進(jìn)而優(yōu)化航行體的水彈道.

［1］顏開，王寶壽.出水空泡流動(dòng)的一些研究進(jìn)展［C］//第二十一屆全國(guó)水動(dòng)力學(xué)研討會(huì).北京:海洋出版社，2008:9-16.

［2］黃壽康.流體動(dòng)力彈道載荷環(huán)境［M］.北京:宇航出版社，1991:165-182.

［3］殷崇一，張宇文，劉樂華，等.導(dǎo)彈水下發(fā)射內(nèi)流場(chǎng)的數(shù)值研究［J］.彈箭與制導(dǎo)學(xué)報(bào)，2003，23(3):56-58.

［4］陳雄州，舒旭光，廉斌，等.水下火箭發(fā)動(dòng)機(jī)噴管出口壓強(qiáng)研究［J］.艦船科學(xué)技術(shù)，2004，26(6):38-40.

［5］盧立秀，湯軍社，劉永超，等.導(dǎo)彈彈射機(jī)構(gòu)的建模及優(yōu)化設(shè)計(jì)［J］.航空計(jì)算技術(shù)，2007，37(6):58-61.

［6］王漢平，吳友生，程棟，等.潛射模擬彈彈射后效分析［J］.船舶力學(xué)，2010，14(10):1122-1128.

［7］仲國(guó)維，張嘉鐘.潛射航行體的水下彈道模擬［J］.彈道學(xué)報(bào)，2005，17(1):8-12.

［8］ACOSTA A J.The effect of a longitudinal gravity field on the supercavitating flow over a wedge［J］.Trans ASME，1961，28(2):188-192.

［9］BASHAROVA V N，BUIVOL V N，SEREBRYAKOV V V.Slender axisymmetric cavities in the flow around bodies in a longitudinal gravity force field［J］.International Applied Mechanics，1983，19(4):359-366.

［10］YAGLA J J.Launch dynamics environment of a water piercing missile launcher［C］//Proceedings of the 24th international symposium on ballistics.New Orleans，LA:NDIA，2008:1-17.

［11］WEILAND C J，VLACHOS P P.Concept analysis and laboratory observations on a waterpiercing missile launcher［J］.Ocean Engineering，2010，3(9):959-965.

［12］劉筠喬，魯傳敬，李杰，等.導(dǎo)彈垂直發(fā)射出筒過程通氣空泡流研究［J］.水動(dòng)力學(xué)研究與進(jìn)展(A輯)，2007，22(5):549-554.

［13］張紅軍，陸宏志.潛射導(dǎo)彈出筒過程三維非定常數(shù)值模擬研究［J］.水動(dòng)力學(xué)研究與進(jìn)展(A輯)，2010，25(3):405-415.

［14］劉樂華，張宇文，袁緒龍.水下大深度垂直發(fā)射內(nèi)流場(chǎng)的數(shù)值研究［J］.水動(dòng)力學(xué)研究與進(jìn)展(A輯)，2005，20(1):90-94.

［15］朱明駿，楊軍.潛射導(dǎo)彈充氣均壓系統(tǒng)建模與控制系統(tǒng)設(shè)計(jì)［J］.戰(zhàn)術(shù)導(dǎo)彈控制技術(shù)，2010，27(1):26-30.

［16］ROUSE H，MCNOWN J S.Cavitation and pressure distribution，head forms at zero angle of yaw［R］.Iowa City:State University of Iowa，1948:23-26.