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

(1.哈爾濱工業(yè)大學航天學院，黑龍江哈爾濱150001;2.河北大學建筑工程學院，河北保定071002)

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

圓柱體在水下運動的時間很短，卻存在復雜的流場情況.Jon J.Yagla等［3］對水下垂直發(fā)射系統(tǒng)的動力環(huán)境進行了研究，提出了筒內發(fā)動機點火的方式發(fā)射導彈，可以為導彈出水提供一個從出筒到出水整個過程的氣幕，并進行了相關的實驗研究.Chris J.Weiland等研究了不同發(fā)射深度和橫向流對垂直發(fā)射導彈的彈道的影響［4］.殷崇一等運用有限體積法和SIMPLE算法計算了導彈水下發(fā)射內流場的數值研究［5］.劉樂華等采用TVD有限體積法，對導彈水下大深度發(fā)射內流場進行了數值模擬，研究流場波系結構、壓力、溫度等參數隨時間變化的規(guī)律［6］.朱明駿等對潛射導彈的充氣均壓系統(tǒng)進行了研究，利用 PID控制方法對充壓系統(tǒng)進行了設計［7］.王漢平等采用多相流模型和動網格技術研究了均壓氣體對彈射后效的影響［8］.劉筠喬等研究了導彈垂直發(fā)射出筒過程的通氣空泡流［9］.張紅軍等對不帶空化模型的潛射導彈垂直出筒過程進行了研究［10］.

國內外主要對圓柱體外流場進行了研究，而對出筒過程中內彈道和流體動力特性的研究相對較少.本文采用動網格技術，并結合求解考慮質量輸運空化模型的混合介質的RANS方程，實現(xiàn)了重力影響下的圓柱體邊界運動與氣(汽)水流場的耦合求解，針對采用筒內點火的方式對圓柱體出筒過程的多相流場進行了數值模擬，研究了重力影響下的3種頭型的圓柱體出筒過程和流體動力特性.

1 圓柱體出筒過程數學模型

1.1 流體控制方程

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

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

式中:ui為i方向的速度分量;ρm為混合物的密度，ρm=αlρl+αvρv+αngρng，ρi和αi(i=l，v，ng)，分別為水、水蒸氣及不可凝結氣體的密度和體積分數.

混合物動量方程:

1.2 湍流模型

標準的k－ε湍流模型:

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

式中:k和ε分別為湍動能和耗散率，Gb由浮力引起的湍動能k的產生項，Gk為由平均速度梯度引起的湍動能k的產生項，μt為湍流粘性系數，σk和σε分別為k和ε的普朗特數.5個經驗常數取值分別為:Cε1= 1.44，Cε2=1.92，σk=1.0，σε=1.3，Cμ=0.09.

1.3 空化模型

式中:Vch為特征速度，表示當地液相和汽相之間的相對速度，本文中取，其中k為當地湍動能;τ為表面張力系數;pv為指定溫度下水的飽和蒸汽壓力，經驗常數: Cc=0.01，Ce=0.02，f1=1－fv－fng，其中不可凝結氣體分數fng=1.5×10－5，fi=αiρi/ρm(i=l，v，ng)為各相的質量分數，不可凝結氣體定義為理想氣體;m·－和m·+分別表征汽化和液化過程的源相.

1.4 能量方程

能量方程為

式中:keff為傳導率，keff=(k+kt)，kt為湍流的熱傳導率;J為物質j的擴散通量;Sh為一切熱源.方程(9)右邊的后3項分別代表由熱傳導、物質擴散和粘性耗散引起的能量交換.

其中，h為焓，理想氣體的焓為

不可壓縮流體的焓為

式中:Yj為物質j的質量分數，

1.5 圓柱體的運動方程

假定水下圓柱體為剛體，由牛頓第二定律計算圓柱的運動，圓柱體邊界運動方程采用如下離散形式:

式中:vi和vi－1分別代表第i和i－1個時間步圓柱體的運動速度，Δt為時間步長，F(xiàn)i－1為第i－1個時間步圓柱體所受的合力，m為圓柱體的質量.合力由下公式算得

式中:FG為圓柱體受到的重力，F(xiàn)B為圓柱體的底部受到的推力，F(xiàn)G和FB由初始條件給出，F(xiàn)P為圓柱體受到的壓差阻力，F(xiàn)V為圓柱受到的粘性阻力，F(xiàn)P和FV從非穩(wěn)態(tài)流場計算得出.

2 數值仿真方法

采用有限體積法對控制方程進行離散，利用Simple方法求解混合介質RANS方程，采用二階迎風格式對動量方程、對流項和湍流耗散項進行離散，采用PRESTO!對壓力項進行離散.利用FLUENT中的自定義函數(UDF)和動網格技術，求解圓柱體出筒過程中的運動速度及位移.

數值方法正確性驗證，將半球頭型彈體的仿真結果與已有的實驗數據進行對比，驗證了數值求解方法的正確性.驗證算例采用模型為:半球頭型彈體直徑為50 mm、總長度為200 mm.計算條件是:未受擾動水流場壓強p∞=378 955 Pa，采用非定常過程計算.以彈體直徑D計算雷諾數為Re=106，空化數σ=0.31，圖1是圓柱體表面壓力分布與實驗數據對比曲線.H為距離彈頭頂點處的軸向距離，從圖中看出仿真結果和實驗數據吻合較好.因此，采用動態(tài)分層方法的動網格技術，結合求解混合介質RANS方程，可以實現(xiàn)彈體邊界運動與汽水流場的耦合求解，并以此為基礎對潛射圓柱體垂直出筒過程空化多相流場開展數值研究.

圖1 數值模擬結果驗證Fig.1 Validity of numerical simulation results

3 計算模型及邊界條件

圖2是圓柱體出筒過程示意圖，整個計算區(qū)域分為2塊:中間為動網格區(qū)域，周圍為靜網格區(qū)域.

圖2 計算區(qū)域與模型Fig.2 Computional zone and model

圓柱體的長細比L/D=10，L為圓柱體的長度，D為圓柱體的直徑，水深為50 m.計算域由筒內的氣體和筒外的水域組成.圖2給出了流場計算區(qū)域、網格和邊界條件的示意圖，在圓柱體頭部附近的網格進行了加密處理.流場的長度是50D，寬度是12D.圓柱體噴管入口是本模型的壓力入口.

圖3 計算區(qū)域與網格Fig.3 Computational domain and grid

圖4是總壓曲線，特征時間T=L/Vmax給出，式中Vmax為出筒的最大速度，Pmax為最大壓力值.

圖4 燃氣總壓隨時間變化的曲線Fig.4 Toll pressure of gas versus time

4 計算結果及分析

對3種頭型的圓柱體出筒過程中的流場、表面壓力、水下彈道和水下流體動力進行計算及分析.

4.1 圓柱體出筒過程流場分析

水下圓柱體運動產生空化的判據用無量綱的量空化數σ來表示，其表達式為，式中v為圓柱體的速度，p∞為未擾動流場的壓力，pv為飽和蒸汽的壓力，ρ為水的密度.圖5是圓柱體空化云圖，模型Ⅰ是小鈍角頭型，模型III是大鈍角頭型.模型III肩部空泡最大，模型Ⅰ肩空泡最小.

圖5 不同結構圓柱體的空化云圖Fig.5 Contours of cavitation to different structare cylinder

圖6是氣相體積分數和汽相體積分數的變化曲線，3種頭型的氣相體積分數曲線，t/T≤1，都平行于軸線且?guī)缀踔睾?，t/T≥1，模型Ⅲ的下降速度最快，模型I則最慢.

圖6 氣相體積分數的變化Fig.6 Gas volume fraction versus time

圖7是汽相體積分數曲線，在t/T≤1時，都平行于軸線且接近0，t/T≥1，模型III的增加速度最快，模型Ⅰ則最慢.在圓柱體頭部由運動產生的動壓較大，不滿足空化的條件，故沒有出現(xiàn)空化現(xiàn)象，在大鈍角頭型的圓柱體的肩部因液體的繞流較容易產生低壓區(qū)，滿足了變化空化的條件，故模型III的肩部空化就容易產生.

圖7 汽相體積分數的變化Fig.7 Volume fraction of vapor versus time

4.2 圓柱體表面壓力分析

圖8是圓柱體出筒時的表面壓力比值的分布曲線，左側是圓柱體頭部，右側是圓柱體尾部，h為距離圓柱體頭部的距離.模型III的肩部出現(xiàn)的繞流低壓區(qū)最大，模型Ⅰ則最小;在圓柱體的中段，模型I的壓力最大，模型III的壓力最小;在圓柱體的尾部，模型I的壓力最大，模型III則最小，3種頭型都有小幅波動，在出筒時刻，圓柱體的尾部還在發(fā)射筒口的高壓氣體之中，氣團的膨脹和收縮是波動的原因.

圖8 圓柱體表面環(huán)向壓力比值沿軸向分布Fig.8 Cylinder surface pressure ratio along axis

4.3 圓柱體出筒彈道分析

圖9 圓柱體速度變化曲線Fig.9 Cylinder velocity versus time

圖9是圓柱體速度變化曲線，從圖中可以看出，模型III的速度最大，模型Ⅰ最小，這是由于航行體出筒過程中模型III的肩部空化最大，圓柱體的沾濕面積最小，摩擦阻力最小，所以速度最大.

圖10是位移變化曲線，從圖中看出，模型III的位移最大，模型Ⅰ則最小，這是由于相同時間內速度越大圓柱體上升的位移越大，模型III的速度最大，故模型III的位移最大.

圖10 圓柱體位移變化曲線Fig.10 Cylinder displacement versus time

圖11 圓柱體加速度變化曲線Fig.11 Cylinder acceleration versus time

圖11是加速度變化曲線，模型III的加速度最大，模型I則最小，3種頭型的加速度在初始刻的波動是由于發(fā)動機的推力逐漸達到峰值和氣團的膨脹和收縮的原因，隨著時間的推移3種頭型的加速度抖動下降，模型III的波動最小，模型I則最大.故模型III的水下彈道變化較小且比較穩(wěn)定，模型I的彈道變化較大穩(wěn)定性較差.通過增加頭型的角度可以得到較理想的彈道，為工程應用提供了理論參考.

4.4 圓柱體出筒過程受力分析

圖12～14分別是壓差系數、粘性系數和總動力系數的變化曲線.圓柱體運動方向為坐標軸x的正向，力的系數用公式Cx=Fx/(0.5ρV2S)給出，式中ρ為水的密度;V為圓柱體的速度;S為圓柱體的橫截面面積;Fx=Fp+Fd，F(xiàn)p為推力，F(xiàn)d為阻力，F(xiàn)d= fg+fw，fg為空氣的粘性力，fw為水的粘性力.出筒過程中的壓差力、粘性力和總動力都為正值.圖12是壓差系數變化曲線，在初始刻有波動，曲線下降且趨緩，模型III的壓差力最大而且波動較小，模型Ⅰ的壓差力則最小而且波動較大.

圖12 圓柱體壓差系數變化曲線Fig.12 Cylinder pressure coefficient versus time

圓柱體的粘性力系數的變化規(guī)律，見圖13，從圖中看出，3種頭型圓柱體的粘性力都有相同的波動然后下降且趨緩，模型III的粘性摩擦力最小而且波動較小，模型I粘性摩擦力較大而且波動較大，肩部空泡的出現(xiàn)增大了圓柱體的壓差阻力系數，減小了粘性摩擦阻力系數.

圖13 圓柱體摩擦系數變化曲線Fig.13 Cylinder friction coefficient versus time

圖14 圓柱體總流體動力系數變化曲線Fig.14 Cylinder toll hydrodynamic coefficient versus time

圓柱體總流體動力系數變化規(guī)律，見圖14，從圖中可以看出，模型III的總流體動力系數最大而且波動較小，模型I的總流體動力則最小且波動較大.由于模型III的總流體動力系數的波動較小，對圓柱體水下彈道穩(wěn)定的影響較小.壓差力、粘性摩擦力和總流體動力系數出現(xiàn)的波動都是由發(fā)動機的推力逐漸增加到峰值和氣團出筒時的膨脹和收縮導致.

5 結論

本文基于混合多相流模型，采用k-ε湍流模型，Singhal空化模型和動網格技術，結合求解考慮質量輸運空化模型的混合介質RANS方程，實現(xiàn)了圓柱體邊界運動與氣水流場的問題的求解.并對重力影響下的不同頭型的水下圓柱體出筒過程進行了數值仿真研究，得到的主要結論如下:

1)在出筒的過程中模型III比模型II和模型Ⅰ產生空化較早并且形成肩部空泡最大，大鈍角易發(fā)生空化，則小鈍角不易發(fā)生空化.

2)在出筒過程中模型III的位移、速度和加速度比模型II和模型Ⅰ都大而且波動較小;模型III出筒過程中水彈道比較穩(wěn)定.

3)模型III流體動力較大而且波動較小，模型Ⅰ流體動力較小而且波動較大.模型III出筒過程中流體動力對其水彈道的穩(wěn)定影響較小.

［1］顏開，王寶壽.出水空泡流動的一些研究進展［C］//第二十一屆全國水動力學研討會暨第八屆全國水動力學學術會議贊兩岸船舶與海洋工程水動力學研討會文集.濟南，中國，2008.

YAN Kai，WANG Baoshou.Some of research progress of water-exit cavity［C］//Proceedings of the 21st National Conference on Hydrodynamics and 8th National Congress on Hydrodynamics and the Cross-Strait Conference on Ship and Ocean Engineering Hydrodynamics.Jinan，China，2008.

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

［3］YAGLA J J.Launch dynamics environment of a water piercing missile launcher［C］//24th International Symposium on Ballistics.New Orleans，2008.

［4］PAVLOS W，YAGLA J J.Concept analysis and laboratory observations on a water piercing missile launcher［J］.O-cean Engineering，2010，3(9):959-965.

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

YIN Chongyi，ZHANG Yuwen，LIU Lehua，et al.A numerical study of the interior flow for underwater launched missile［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2003，23(3):56-58.

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

LIU Lehua，ZHANG Yuwen，YUAN Xulong.Numerical simulation of the interior flow for launch under deep water［J］.Journal of Hydrodynamics:Ser A，2005，20(1):90-94.

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

ZHU Mingjun，YANG Jun.Modeling and design of pressure servo system used in underwater launched missile［J］.Control Techonology of Tactical Missile，2010，27(1):26-30.

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

WANG Hanping，WU Yousheng，CHENG Dong，et al.Analysis of ejection after-effect for underwater-launched emulating missile［J］.Journal of Ship Mechanics，2010，14 (10):1112-1128.

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

LIU Yunqiao，LU Chuanjing，LI Jie，et al.An investigation of ventilated cavitating flow in vertical launching of a missile［J］.Chinese Journal of Hydrodynamics:Ser A，2007，22 (5):549-554.

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

ZHANG Hongjun，LU Hongzhi.Numerical investigation to vertical launching of submarine launching missile［J］.Chinese Journal of Hydrodynamics:Ser A，2010，25(3): 405-415.