毛聰聰，阮文俊，孫雪明，賈昊楠

(1.南京理工大學 能源與動力工程學院，江蘇 南京 210094；2.江西星火軍工工業(yè)有限公司,江西 南昌 331709；3.中國人民解放軍32181部隊,河北 石家莊 050000)

火箭彈發(fā)射過程中，發(fā)射筒會經(jīng)受高溫、高壓、高速燃氣射流的沖擊。在對火箭發(fā)射筒進行力學分析時，也需將發(fā)射筒內(nèi)表面壓力數(shù)據(jù)作為邊界條件，為不同狀態(tài)下的發(fā)射筒有限元分析提供載荷數(shù)據(jù)[1-4]。尤其對帶導轉(zhuǎn)槽的火箭武器系統(tǒng)進行發(fā)射過程分析時，需要考慮由于導轉(zhuǎn)槽存在而導致的不對稱力帶來的影響，因此有必要對帶導轉(zhuǎn)槽的火箭發(fā)射筒三維內(nèi)流場及發(fā)射裝置內(nèi)的不對稱力進行研究。

國內(nèi)外許多科研人員對燃氣射流的數(shù)值模擬做了大量研究。張磊等[5]通過大渦模擬研究了不同噴管口徑對燃氣射流流場的影響。鄧春麗等[6]建立了導彈發(fā)射過程二維流場計算模型，研究了不同時刻發(fā)射筒的內(nèi)流場和外流場特征。張俊等[7]通過研究某固體火箭發(fā)動機燃氣流場的分布，得到了流場軸向與徑向的燃氣射流參數(shù)分布規(guī)律。王華等[8]建立了多噴管流場計算模型，研究了燃氣射流流場的速度與溫度分布。Sinha等[9]對燃氣射流沖擊發(fā)射筒過程進行數(shù)值模擬，將不同網(wǎng)格數(shù)量和湍流模型得到的計算結(jié)果進行了對比分析。

為了使數(shù)值模擬更加接近實際的發(fā)射過程，筆者采用動網(wǎng)格技術(shù)對火箭彈發(fā)射時的運動過程進行模擬[10-12]。通過建立火箭彈發(fā)射過程發(fā)射筒內(nèi)流場計算模型，使用SSTk-ω湍流模型對非穩(wěn)態(tài)流場進行數(shù)值模擬，并對帶導轉(zhuǎn)槽的火箭發(fā)射筒內(nèi)流場和非對稱力進行研究。

1 理論基礎(chǔ)

1.1 控制方程

質(zhì)量守恒方程：

(1)

動量方程：

(2)

能量方程：

(3)

1.2 湍流模型

筆者采用的湍流模型為SSTk-ω湍流模型，與標準k-ω模型相比，該模型在湍流粘度定義中考慮了湍流剪切應(yīng)力的輸運過程，是對標準k-ω模型的增強。表達式為

(4)

(5)

式中：GK表示湍流的動能；Gω為ω的方程；Γk和Γω分別代表k與ω的發(fā)散項；Dω代表正交發(fā)散項。

2 計算模型及網(wǎng)格劃分

2.1 發(fā)射筒實體模型及計算模型

某單管火箭武器發(fā)射筒的實體模型如圖1所示。將發(fā)射筒尾部的凸臺結(jié)構(gòu)簡化，并將導轉(zhuǎn)槽簡化為連續(xù)的螺旋槽后，根據(jù)實體模型建立的三維仿真計算物理模型如圖2所示。

因為不需要考慮火箭彈受到燃氣射流的沖擊狀況，所以將火箭彈簡化為噴管。在發(fā)射筒力學分析時由于需要考慮導轉(zhuǎn)槽處的應(yīng)力集中，因此可以通過三維帶導轉(zhuǎn)槽的發(fā)射筒內(nèi)流場計算，得到螺旋導轉(zhuǎn)槽處不同時刻表面的壓力分布。

2.2 初始及邊界條件

流場初始化壓力大小為101.325 kPa,溫度為300 K。因涉密原因，文中噴管的入口壓力和火箭彈運動速度都經(jīng)過無量綱化處理。將噴管進口設(shè)置為壓力入口，入口燃氣溫度為2 500 K?；鸺龔棸l(fā)射時的入口壓力隨時間變化曲線由實驗測定，選取實驗曲線中壓力到達0.036 6的時刻作為計算的起始時刻，火箭彈運動出筒時刻作為計算的終止時刻，共歷時130 ms。從開始計算到火箭彈運動出筒的噴管入口壓力時間曲線如圖3所示。將發(fā)射筒壁面設(shè)定為絕熱無滑移壁面，發(fā)射筒尾端設(shè)置為壓力出口，壓力大小為101.325 kPa，溫度為300 K。根據(jù)噴管入口壓力隨時間變化曲線得到火箭彈運動速度隨時間變化的曲線，如圖4所示。

2.3 計算條件設(shè)定

流場計算選擇基于壓力的瞬態(tài)求解器，將燃氣當作理想氣體可壓縮流處理。不考慮火箭彈發(fā)射時的空氣阻力；不考慮燃氣化學反應(yīng)和固體顆粒的影響；不考慮燃氣與發(fā)射筒內(nèi)壁及周圍空氣的傳熱；不考慮火箭彈的旋轉(zhuǎn)，僅考慮火箭彈沿水平方向的運動。流體粘性系數(shù)采用Sutherland公式計算。

2.4 網(wǎng)格劃分

在劃分網(wǎng)格時，由于較難生成單塊的高質(zhì)量結(jié)構(gòu)化網(wǎng)格，因此采用分塊網(wǎng)格對計算域進行劃分，4個網(wǎng)格劃分區(qū)域分別為噴管、尾部計算域、導轉(zhuǎn)槽和筒體。筆者采用Icem軟件劃分六面體結(jié)構(gòu)網(wǎng)格，對邊界區(qū)域添加了邊界層網(wǎng)格后,初始時刻的網(wǎng)格總數(shù)為60萬個。初始時刻網(wǎng)格劃分如圖5所示。

2.5 動網(wǎng)格設(shè)置

為了處理火箭彈發(fā)射時的網(wǎng)格變化，在流場計算時通過動網(wǎng)格對火箭彈運動進行模擬。噴管網(wǎng)格水平移動，尾部計算域網(wǎng)格被拉長，筒體與導轉(zhuǎn)槽網(wǎng)格固定不動。動網(wǎng)格更新方法包括彈性光順法、動態(tài)層法和局部重構(gòu)法，根據(jù)本文的網(wǎng)格劃分方式，采用動態(tài)層法模擬火箭彈運動。該方法的基本思想是：根據(jù)與運動邊界相鄰網(wǎng)格的高度，判斷網(wǎng)格需要分割還是合并，當網(wǎng)格高度滿足式(6)時，網(wǎng)格將會被分割；而當網(wǎng)格高度滿足式(7)時，兩層網(wǎng)格將會合并。

h≥(1+α)h0，

(6)

h≤βh0，

(7)

式中：α為網(wǎng)格層的分裂因子；β為網(wǎng)格層的合并因子；h0為理想單元網(wǎng)格高度。

筆者取h0=3 mm,α=0.4，β=0.2。通過Profile文件將火箭彈的運動速度賦給噴管及尾部計算域，噴管不斷向前運動，尾部計算域網(wǎng)格通過動態(tài)層法進行更新生成。初始時刻的運動域網(wǎng)格如圖6所示，火箭彈發(fā)射時運動域的網(wǎng)格變化如圖7所示。

3 數(shù)值模擬結(jié)果及分析

3.1 不同時刻流場縱截面壓力云圖

計算時刻分別為50、90和130 ms時的流場縱截面壓力云圖，如圖8所示。

因涉密原因，流場計算得到的壓力、速度及扭轉(zhuǎn)力都經(jīng)過了無量綱化處理。由于噴管內(nèi)壓力變化較大，因此為了更清楚顯示發(fā)射筒內(nèi)的壓力分布，將無量綱壓力在0.366以上的部分用相同顏色顯示。燃氣在噴管內(nèi)由亞音速加速到超音速后,在噴管出口形成了超音速燃氣射流。按膨脹程度劃分，燃氣射流為欠膨脹射流。燃氣射流從噴管出口流出時的壓力高于周圍環(huán)境壓力，所以一出噴管首先會膨脹加速。燃氣經(jīng)過噴管出口兩側(cè)的膨脹扇區(qū)后,繼續(xù)膨脹加速同時壓力不斷降低。當噴管出口的膨脹波系到達射流邊界后，向內(nèi)反射為一系列的壓縮波，壓縮波足夠強時會匯合形成斜激波。燃氣經(jīng)過激波后，壓力增加、溫度升高同時速度降低。當燃氣壓力增加到高于環(huán)境壓力后，燃氣膨脹加速同時壓力降低。隨著這樣膨脹壓縮過程的反復循環(huán)，發(fā)射筒內(nèi)形成了一系列膨脹波與壓縮波結(jié)構(gòu)。

3.2 中心軸線上的速度與壓力分布

圖9和圖10分別為不同時刻流場中心軸線上的速度與壓力分布。從圖9和圖10中可以看出，燃氣射流從噴管流出后，噴管出口附近的膨脹波和壓縮波強度是最大的，因此燃氣在噴管附近的速度與壓力波動也最為明顯。在50 ms時，燃氣在中心軸線上經(jīng)歷了1～2個膨脹-壓縮波；在90 ms時，燃氣在中心軸線上經(jīng)歷了5～6個膨脹-壓縮波；在130 ms時，燃氣在中心軸線上經(jīng)歷了12～13個膨脹-壓縮波。因發(fā)射筒長度有限，所以當燃氣到達發(fā)射筒出口附近時，燃氣的速度與壓力仍有較大波動。由于粘性效應(yīng)和周圍空氣的作用產(chǎn)生了湍流剪切層，湍流剪切層內(nèi)的粘性阻尼使得膨脹波與壓縮波結(jié)構(gòu)不斷瓦解，因此中心軸線上的燃氣速度與壓力呈現(xiàn)振蕩衰減的趨勢。

3.3 不同時刻發(fā)射筒內(nèi)表面壓力分布

分別取計算時刻為50、90和130 ms時的發(fā)射筒內(nèi)表面壓力云圖，如圖11所示。燃氣從噴管流出后，位于噴管出口后的發(fā)射筒區(qū)域直接受到燃氣射流的強烈沖擊，因此其表面壓力明顯高于噴管前部未直接受到燃氣射流沖擊部分的壓力。

當燃氣經(jīng)歷了在壁面附近的激波后壓力增加，其對發(fā)射筒壁面的壓力也隨之增加。50 ms時，燃氣對發(fā)射筒壁面的最大壓力為0.078；90 ms時，燃氣對發(fā)射筒壁面的最大壓力為0.065；130 ms時，燃氣對發(fā)射筒壁面的最大壓力為0.049。當燃氣壓力增加到高于周圍空氣的壓力后，燃氣膨脹加速并且壓力不斷降低，其對壁面的壓力也逐漸降低。燃氣在發(fā)射筒內(nèi)經(jīng)歷了一系列膨脹壓縮波，所以燃氣對發(fā)射筒內(nèi)表面的壓力在整體上呈現(xiàn)振蕩衰減的趨勢。

3.4 不同時刻導轉(zhuǎn)槽內(nèi)的壓力及扭轉(zhuǎn)力分布

不同時刻導轉(zhuǎn)槽表面壓力隨距離變化曲線如圖12所示。

在發(fā)射筒右視圖中，位于導轉(zhuǎn)槽左側(cè)的壁面為本文的導轉(zhuǎn)槽左側(cè)面，位于導轉(zhuǎn)槽右側(cè)的壁面為本文的導轉(zhuǎn)槽右側(cè)面。50 ms時，燃氣對導轉(zhuǎn)槽頂面、右側(cè)面和左側(cè)面的最大壓力分別為0.072，0.070和0.065；90 ms時，燃氣對導轉(zhuǎn)槽頂面、右側(cè)面和左側(cè)面的最大壓力分別為0.065，0.060和0.055；130 ms時，燃氣對導轉(zhuǎn)槽頂面、右側(cè)面和左側(cè)面的最大壓力分別為0.045，0.044和0.042。因50 ms時噴管入口壓力是3個時刻中最大的，故燃氣對導轉(zhuǎn)槽3個表面的最大壓力也是3個時刻中最大的。在不同時刻，導轉(zhuǎn)槽頂面受到的最大壓力值始終是3個表面中最大的，而導轉(zhuǎn)槽右側(cè)面受到的最大壓力略大于導轉(zhuǎn)槽左側(cè)面。

從不同時刻導轉(zhuǎn)槽表面的壓力分布曲線中可以看出，導轉(zhuǎn)槽表面的壓力沿著中心軸線方向振蕩變化，但并未呈現(xiàn)振蕩衰減的趨勢。這是由于導轉(zhuǎn)槽表面會多次反射壓縮波，當這些反射的壓縮波疊加后，燃氣再次經(jīng)過這些壓縮波會增加對導轉(zhuǎn)槽表面的壓力，因此導轉(zhuǎn)槽表面的壓力并未沿著軸線方向振蕩衰減。

導轉(zhuǎn)槽內(nèi)扭轉(zhuǎn)力隨時間變化的曲線如圖13所示。導轉(zhuǎn)槽兩個側(cè)面受到燃氣射流直接作用面積隨著火箭彈向前運動而逐漸增加，故導轉(zhuǎn)槽兩個側(cè)面受到的燃氣作用力也在不斷增加。導轉(zhuǎn)槽兩個側(cè)面上受到的燃氣射流作用力差值使得發(fā)射裝置內(nèi)產(chǎn)生了不對稱力即扭轉(zhuǎn)力，而右側(cè)面受到的燃氣射流作用力始終大于左側(cè)面，因此導轉(zhuǎn)槽內(nèi)扭轉(zhuǎn)力的方向始終垂直于右側(cè)面。扭轉(zhuǎn)力隨著時間呈現(xiàn)振蕩上升的趨勢，在火箭彈即將出筒時刻達到最大，最大的扭轉(zhuǎn)力值為769。

4 結(jié)論

筆者通過動網(wǎng)格技術(shù)模擬了火箭彈從發(fā)射到出筒過程燃氣射流對發(fā)射筒的沖擊。通過流場計算得到不同時刻發(fā)射筒內(nèi)表面的壓力分布，將其加載到包含導轉(zhuǎn)槽結(jié)構(gòu)的三維發(fā)射筒有限元分析模型中進行計算?？梢缘玫桨l(fā)射筒在發(fā)射壓力作用下的應(yīng)力應(yīng)變情況，從而確定發(fā)射筒在發(fā)射壓力作用下的相對薄弱部位以及安全儲備。從本文的數(shù)值模擬中可以得到以下結(jié)論：

1)超音速燃氣射流與環(huán)境大氣之間的壓力差，導致發(fā)射筒內(nèi)產(chǎn)生了一系列膨脹壓縮波。由于流體粘性效應(yīng)和周圍大氣的作用，膨脹波與壓縮波結(jié)構(gòu)在發(fā)射筒內(nèi)逐漸崩解。

2)燃氣在發(fā)射筒內(nèi)經(jīng)歷了一系列膨脹波和壓縮波，因此不同時刻中心軸線上的燃氣壓力、燃氣速度和發(fā)射筒內(nèi)壁面上的壓力，在整體上呈現(xiàn)振蕩衰減的趨勢。

3)導轉(zhuǎn)槽內(nèi)壓縮波的反射疊加，使得不同時刻導轉(zhuǎn)槽表面上的壓力并未沿著中心軸線振蕩衰減。導轉(zhuǎn)槽的存在導致火箭發(fā)射筒內(nèi)產(chǎn)生了扭轉(zhuǎn)力，并且隨著火箭彈向前運動，扭轉(zhuǎn)力振蕩增加。