水下航行體垂直出筒流體動力特性研究

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

(哈爾濱工業大學航天學院，150001哈爾濱)

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

潛射航行體在水下運動的時間很短，卻存在復雜的流場情況，殷崇一等［3］運用有限體積法和SIMPLE算法計算了航行體水下發射內流場的數值研究.陳雄州等［4］根據固體發動機推力最大原則，并結合水下航行體外彈道運動方程，求出了不同使用環境下的發動機最佳出口壓強.盧立秀等［5］運用MATLAB對飛機內彈道彈射機構進行了研究.王漢平等［6］采用多相流模型和動網格技術研究了均壓氣體對航行體射后效的影響.仲維國等［7］基于魚雷的航行力學導出水下航行體的彈道方程.Acosta等［8］采用數值模擬的方法對重力場中的軸對稱細長體沿軸線的定常空泡長度、空泡截面面積和阻力系數進行了研究，得到了重力場中弗洛德數對空泡的長度和面積影響較大，對軸對稱細長體的阻力系數影響較小.Basharova等［9］采用基于細長體空泡理論的近似方程組計算了有重力影響無限流場中沿軸線方向上空泡的形態和基本尺寸.J.J.Yagla等［10］對水下垂直發射系統的動力環境進行了研究，提出了筒內發動機點火的方式可以為航行體的出水提供1個從出筒到出水整個過程的氣幕，并進行了相關的實驗.Chris J.Weiland等［11］研究了不同水深和橫向流對氣幕的影響.劉筠喬等［12］研究了航行體垂直發射出筒過程的通氣空泡流.張紅軍等［13］對不帶空化模型的潛射航行體垂直出筒過程進行了研究.劉樂華等［14］采用TVD有限體積法研究了深水發射的內流場，得到燃氣流噴出之后快速膨脹，進而形成高溫燃氣射流，對反射板有極大的破壞作用.朱明駿等［15］對潛射航行體的充氣均壓系統進行了研究，并利用PID控制方法對充壓系統進行了設計.

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

1 數學模型

1.1 控制方程

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

混合物連續性方程為

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

混合物動量方程為

式中:p是混合介質的壓力;μm是運動粘性系數，其中n是構成混合介質的相數，μk是每個相的湍流粘性系數;ρm為混合物的密度;g是重力加速度.

1.2 湍流模型

標準的k-ε湍流模型為

湍動能k和湍動耗散率ε的輸運方程如下:

式中:k是混合介質的湍動能;ε是混合介質湍流耗散率;Gb是浮力產生湍動能項;Gk是速度梯度產生的湍動能項;μt為湍流粘性系數;ρm為混合物的密度;σk和σε分別為k和ε的普朗特數;可以調整的經驗常數取值分別為:Cε1=1.44，Cε2=1.92，σk=1.0，σε=1.3，Cμ=0.09.

1.3 空化模型

基于均質平衡流理論的Singhal空化模型，表達式如下:

式中:Vch是液相和汽相之間轉化的特征速度;τ是液體表面張力;pv是飽和蒸汽壓;fl=1-fv-fng，其中不可凝結氣體分數fng=1.5×10-5;Cc=0.01和Ce=0.02是兩個經驗常數.

1.4 能量方程

能量方程為

式中:keff為傳導率，keff=(k+kt);kt為湍流的熱傳導率為物質j的擴散通量.上式中等號右邊的最后3項分別代表由熱傳導、物質擴散和粘性耗散引起的能量交換.Sh包含化學反應熱在內的一切熱源.

1.5 航行體的運動學方程

航行體的運動學方程由牛頓第二定律得到，對航行體運動學方程使用下式進行離散:

式中:Vi是航行體在第i個時間步的速度;Vi-1是航行體在第i-1個時間步的速度;Δt是時間步長;Fi-1是航行體在第i-1個時間步合外力，m是航行體的質量.合力由下式算得:

式中:FG為航行體受到的重力;FB為航行體的底部受到的推力;FG和FB由初始條件給出;FP是航行體前后壓強差產生阻力;FV是航行在粘性流體運動受到的粘性力阻力，FP和FV均是通過對流場的計算得到.

2 數值模擬

2.1 數值求解及動網格

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

采用動網格技術進行計算時，動網格區域中的網格具有一定對流運動，在實際計算過程中要去掉動網格對流運動的影響，動網格守恒方程統一形式如下:

式中:ρm為混合物的密度;u為流體的速度矢量;ug為網格的運動速度;Γ為擴散系數;Sφ為標量φ的源相;?V為控制體積V的邊界.

2.2 數值求解方法及正確性

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

圖1 實驗結果與仿真結果的比較

在計算過程中采用網格重構法，結合多相流Mixture模型和湍流模型，對流場混合介質RANS方程進行求解，進而求解固體運動邊界與汽水流場的流-固耦合問題，采用此方法對水下航行體垂直發射出筒過程流體動力特性展開數值模擬研究.

3 計算模型及邊界條件

航行體的長細比L/D=10，L為航行體的長度，D為航行體的直徑，頭為半球型，發射深度為50 m.計算域由筒內的氣體和筒外的水域組成.圖2給出了流場計算區域、網格和邊界條件的示意圖，在航行體附近的網格進行了加密處理.圖3給出了靜網格和動網格示意圖，整個計算域分為兩塊:中間為動態區域，網格隨時間變化，周圍為靜態區域，網格固定不動.流場的長度是航行體長度的5倍，寬度是航行體直徑的12倍.

航行體噴管入口是本模型的壓力入口，噴管入口的壓力采用圖4給出隨時間變化燃氣總壓曲線.

圖2 計算區域與網格

圖3 計算區域與模型

圖4 燃氣總壓隨時間變化曲線

4 計算結果及分析

4.1 流場分析

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

圖5 不同時刻流場體積分數ai云圖

4.2 流場壓力分布

圖6給出了不同時刻流場的壓力云圖，所選的時間點與圖5一一對應，從壓力云圖中可以看出，水域流場的壓力分布主要受重力的影響.初始時刻，筒口附近水域壓力比較大，發動機點火筒內的壓力開始增加.隨著出筒高度的增加周圍環境壓力開始減小.從IV時刻后，在航行體的肩部出現了繞流引起的低壓區，隨出筒高度的增加，航行體的速度增大，在肩部繞流引起的低壓區繼續增大，這與圖5的壓力云圖中的肩部空泡的產生是對應的.航行體尾部離開發射筒后，筒口附近的海水壓力沿徑向往筒口逐漸增加，但筒口氣團的壓力仍高于環境壓力，此時筒口周圍的水暫時沒有向筒內回灌.

圖6 不同時刻流場壓力云圖

4.3 航行體表面壓力分析

圖7給出了不同時刻航行體表面的環向壓力沿軸向的分布，橫坐標代表的是航行體的特征長度，縱坐標代表的是表面壓力與最大壓力的比值，h是距離航行體頭部的高度，L是航行體的長度，T=L/V0，V0是出筒的速度，T是特征時間.左側是航行體的頭部，右側是航行體的尾部.從圖中看出，隨著出筒時間的增加，航行體表面的壓力系數是逐漸降低的，在航行體的頭部出現了由繞流引起的低壓區，在航行體的尾部由于處在發射筒內高壓燃氣之中，出現了高壓區，隨出筒時間的增加航行體尾部壓力逐漸減小并趨緩，這與圖6中壓力場的分布一致.

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

4.4 航行體出筒彈道分析

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

圖9是加速度變化曲線.加速度曲線在初始刻有上下波動，這是由于發動機的推力在逐漸增加至峰值時，加速度也迅速增大，然而在筒口，由于氣團的膨脹和收縮導致加速度又有所波動，最后加速度抖動下降.

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

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

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

圖10 總系數變化曲線

5 結論

采用數值模擬的方法研究了有重力影響的水下航行體垂直發射出筒過程流體動力特性，得到如下結論:

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

2)氣團的形狀從半球型變化到橢球型，最后變成細長型;航行體頭部在海水靜壓作用下出現了高壓和在肩部出現了繞流引起的低壓區;

3)航行體頭部穿透氣團之前，氣團的膨脹和收縮引起航行體的流體動力系數較大的波動，縮短航行體穿透氣團的時間可以減弱航行體載荷的波動，進而優化航行體的水彈道.

［1］顏開，王寶壽.出水空泡流動的一些研究進展［C］//第二十一屆全國水動力學研討會.北京:海洋出版社，2008:9-16.

［2］黃壽康.流體動力彈道載荷環境［M］.北京:宇航出版社，1991:165-182.

［3］殷崇一，張宇文，劉樂華，等.導彈水下發射內流場的數值研究［J］.彈箭與制導學報，2003，23(3):56-58.

［4］陳雄州，舒旭光，廉斌，等.水下火箭發動機噴管出口壓強研究［J］.艦船科學技術，2004，26(6):38-40.

［5］盧立秀，湯軍社，劉永超，等.導彈彈射機構的建模及優化設計［J］.航空計算技術，2007，37(6):58-61.

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

［7］仲國維，張嘉鐘.潛射航行體的水下彈道模擬［J］.彈道學報，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］劉筠喬，魯傳敬，李杰，等.導彈垂直發射出筒過程通氣空泡流研究［J］.水動力學研究與進展(A輯)，2007，22(5):549-554.

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

［14］劉樂華，張宇文，袁緒龍.水下大深度垂直發射內流場的數值研究［J］.水動力學研究與進展(A輯)，2005，20(1):90-94.

［15］朱明駿，楊軍.潛射導彈充氣均壓系統建模與控制系統設計［J］.戰術導彈控制技術，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.