李映坤， 韓珺禮,2， 陳雄,*， 周長省， 鞏倫昆

1.南京理工大學(xué) 機械工程學(xué)院， 南京 210094

2.北京機電研究所， 北京 100083

基于多物理場耦合的雙脈沖發(fā)動機點火過程數(shù)值模擬

李映坤1， 韓珺禮1,2， 陳雄1,*， 周長省1， 鞏倫昆1

1.南京理工大學(xué) 機械工程學(xué)院， 南京 210094

2.北京機電研究所， 北京 100083

為研究雙脈沖固體火箭發(fā)動機Ⅱ脈沖點火瞬態(tài)過程，發(fā)展一套多物理場耦合求解器。流體控制方程基于有限體積法求解，時間推進(jìn)采用雙時間步LU-SGS(Lower Upper Symmetric Guass Seidel)方法；固體推進(jìn)劑表面溫度基于耦合傳熱方法計算；結(jié)構(gòu)動力學(xué)運動方程基于有限元方法離散，采用經(jīng)典的Newmark格式進(jìn)行時間推進(jìn)，流固耦合采用松耦合算法，并通過算例驗證求解器的可靠性。計算結(jié)果表明：該求解器能夠數(shù)值模擬Ⅱ脈沖啟動過程中的點火藥氣體沖擊、燃?xì)夥嵌ǔＡ鲃蛹敖饘倌て瑱C械響應(yīng)過程，獲得金屬膜片的破裂時間和壓強；且隨著點火質(zhì)量流率增加，推進(jìn)劑裝藥首次點燃時間和金屬膜片破裂時間變短，膜片破裂壓強降低；金屬膜片破裂時間和壓強不僅與作用在其表面的壓強載荷大小相關(guān)，而且與壓強載荷加載的過程相關(guān)；金屬膜片厚度越薄，膜片破裂時間越短，膜片軸向位移越大，膜片破裂壓強越低。

多物理場耦合； 流固耦合; 耦合傳熱; 點火; 雙脈沖發(fā)動機; 固體火箭發(fā)動機; 數(shù)值模擬

雙脈沖固體火箭發(fā)動機(以下簡稱雙脈沖發(fā)動機)由兩級燃燒室組成，兩級燃燒室之間由脈沖隔離裝置連接，且共用一個噴管，兩級脈沖點火時間間隔可以根據(jù)戰(zhàn)術(shù)指標(biāo)進(jìn)行調(diào)整。與傳統(tǒng)的固體火箭發(fā)動機相比，雙脈沖發(fā)動機具有可多次點火、提供不連續(xù)推力、推進(jìn)系統(tǒng)能量可控等優(yōu)點[1]，是現(xiàn)有飛行系統(tǒng)中的一種先進(jìn)動力裝置。

脈沖隔離裝置是雙脈沖發(fā)動機的關(guān)鍵部件之一，各國學(xué)者對不同類型的脈沖隔離裝置進(jìn)行了全面系統(tǒng)的研究。Nishii[2]、Carrier[3]、Dahl[4]、Wang[5]、Schilling[6]和Stadler[7]等分別針對隔塞式、陶瓷艙蓋式、金屬膜片式、軟質(zhì)隔層式脈沖隔離裝置展開深入研究，并對其工作特性進(jìn)行實驗驗證。其中，金屬膜片式脈沖隔離裝置兼具有結(jié)構(gòu)設(shè)計簡單、研制周期短和可靠性高等優(yōu)點，在國內(nèi)外被廣泛應(yīng)用，其結(jié)構(gòu)簡圖如圖1所示，該發(fā)動機由金屬膜片式脈沖隔離裝置、兩級脈沖點火具、兩級燃燒室、兩級固體推進(jìn)劑裝藥和噴管組成。

圖1 雙脈沖發(fā)動機結(jié)構(gòu)簡圖
Fig.1 Schematic diagram of dual pulse motor

金屬膜片式脈沖隔離裝置主要由支撐件、金屬膜片和隔熱層組成。為控制金屬膜片破裂時間和壓強，膜片一側(cè)一般設(shè)置有缺陷槽。

近年來，國內(nèi)外學(xué)者相繼對金屬膜片式雙脈沖發(fā)動機工作過程進(jìn)行了仿真和實驗研究。Javed等[8]對雙脈沖發(fā)動機工作過程中三維內(nèi)流場進(jìn)行了數(shù)值模擬，分析速度、壓力和溫度沿著軸線的分布；孫娜等[9]指出由于脈沖隔離裝置的存在，使得燃?xì)饬髟冖衩}沖燃燒室內(nèi)出現(xiàn)后臺階流動，氣流發(fā)生分離再附著過程，氣流再附著點附近為絕熱層燒蝕較為劇烈的部位；王春光等[10]設(shè)計一種金屬膜片式隔艙結(jié)構(gòu)，基于ABAQUS商業(yè)軟件，利用脆性斷裂模型模擬了膜片的破壞過程，并將有限元數(shù)值計算結(jié)果、理論公式計算結(jié)果和單項實驗驗證結(jié)果進(jìn)行對比；石瑞等[11]提出了含缺陷槽鋁膜隔板及其組件的設(shè)計方法，對鋁膜隔板破裂特性進(jìn)行了數(shù)值模擬，并對不同厚度和刻痕深度的鋁膜隔板進(jìn)行耐壓和破裂實驗；劉偉凱等[12-13]通過三維虛擬裂紋閉合法數(shù)值計算預(yù)制缺陷處的應(yīng)力強度因子，通過多孔圓板強度理論建立支撐件強度校核方法，基于ABAQUS商業(yè)軟件，選取Ductile damage 模型和Brittle cracking 模型，對膜片靜態(tài)和動態(tài)打開過程進(jìn)行仿真，最后通過膜片冷流靜態(tài)和熱流動態(tài)實驗進(jìn)一步驗證計算方法的準(zhǔn)確性；王偉等[14]通過圓板大撓度理論和斷裂力學(xué)理論推導(dǎo)出金屬膜片預(yù)制缺陷處應(yīng)力強度因子計算公式，提出金屬膜片在內(nèi)壓作用下的設(shè)計方法，并通過實驗驗證其正確性；李映坤等[15]采用經(jīng)驗公式計算再附著點處的對流換熱系數(shù)，分析脈沖隔離裝置通道孔徑、寬度和角度對再附著點位置和對流換熱系數(shù)的影響；陳雄等[16]基于耦合傳熱算法，分析隔離裝置通道孔徑對雙脈沖發(fā)動機Ⅰ脈沖燃燒室熱防護(hù)層受熱的影響。

以上對于金屬膜片式雙脈沖發(fā)動機的研究，并沒有考慮Ⅱ脈沖點火瞬態(tài)過程，也未見有關(guān)Ⅱ脈沖點火過程仿真研究或成果的公開報道。實際上，Ⅱ脈沖發(fā)動機點火過程非常復(fù)雜，為典型的多物理場耦合問題，按照時間順序可以分為：點火具噴射點火藥氣體、推進(jìn)劑裝藥表面受熱著火(流熱耦合)、火焰?zhèn)鞑?、Ⅱ脈沖燃燒室增壓、脈沖隔離裝置金屬膜片受壓破壞(流固耦合)、燃?xì)馔ㄟ^脈沖隔離裝置通道、Ⅱ脈沖燃燒室泄壓和Ⅰ脈沖燃燒室增壓。Ⅱ脈沖點火過程關(guān)系整個發(fā)動機工作的可靠性，但是由于實驗費用、實驗測試設(shè)備及測試方法的限制，使用數(shù)值模擬方法進(jìn)行Ⅱ脈沖點火過程的研究顯得更為突出。

本文發(fā)展一套多物理場耦合求解器，對Ⅱ脈沖啟動過程中點火藥氣體沖擊過程、推進(jìn)劑裝藥受熱點火過程、燃?xì)獾姆嵌ǔＡ鲃蛹敖饘倌て瑱C械響應(yīng)過程進(jìn)行數(shù)值模擬，獲得金屬膜片破裂時間和壓強，分析點火具質(zhì)量流率和金屬膜片厚度對發(fā)動機點火瞬態(tài)過程的影響。

1 控制方程和計算方法

1.1 流體區(qū)域

雙脈沖發(fā)動機點火過程中流體計算區(qū)域不斷變化，控制方程中需考慮網(wǎng)格運動。采用ALE(Arbitrary Lagrangian Eulerian)方法描述的可壓縮非定常雷諾時均Navier-Stokes方程為

?ΩHdV+?ΩSdV

(1)

式中：V為體積；s為面積；t為時間；?Ω為某一固定區(qū)域Ω的邊界；U為守恒變量；Fc為無黏通量；Fv為黏性通量；H為軸對稱幾何源項；n為控制體表面法向量，以上各式的具體形式見文獻(xiàn)[17]；S為固體推進(jìn)劑表面燃燒加質(zhì)源項[18]，包含質(zhì)量、動量和能量源項，只發(fā)生在靠近推進(jìn)劑表面的第1層流體網(wǎng)格，其具體表達(dá)式為

(2)

相對于k-ε和k-ω湍流模型，k-ωSST(Shear Stress Transport)湍流模型對近壁處和后臺階流場預(yù)測更加準(zhǔn)確。因此，本文采用Menter提出的k-ωSST湍流模型，具體方程見文獻(xiàn)[19]。

采用基于格心的多塊結(jié)構(gòu)網(wǎng)格有限體積法求解流體控制方程，對流通量離散采用三階MUSCL(Monotone Upstream centered Schemes for Conservation Laws)重構(gòu)和AUSMPW+(Advection Upstream Splitting Method by Pressure-based Weight functions)格式[20]，黏性項離散采用二階中心差分格式，時間推進(jìn)采用雙時間步LU-SGS(Lower Upper Symmetric Guass Seidel)時間離散方法[21]。

1.2 固體結(jié)構(gòu)區(qū)域

基于有限元方法的基本原理和基本步驟，結(jié)構(gòu)動力學(xué)運動方程可以寫為

(3)

式中：q(t)是單元節(jié)點的位移矢量，q(t)對時間t的一階和二階導(dǎo)數(shù)分別表示單元節(jié)點的速度和加速度；p(t)為作用在結(jié)構(gòu)上的外載荷；M為質(zhì)量矩陣；C為阻尼矩陣，一般采用瑞利阻尼方法進(jìn)行近似；K為剛度矩陣。

結(jié)構(gòu)動力學(xué)運動方程的求解采用經(jīng)典的Newmark 方法[22]。給定時間間隔[t,t+Δt]，已知t時刻的q(t)及其導(dǎo)數(shù)，t+Δt時刻的位移、速度和加速度由下式計算：

(4)

式中：α和β為兩常數(shù)，下標(biāo)t和t+Δt分別代表時刻t、t+Δt的值。

1.3 固體推進(jìn)劑點火模型

固體推進(jìn)劑的點火過程采用廣泛使用及認(rèn)可的Zeldovich-Novozhilov (ZN)點火模型[23]，并通過耦合傳熱的方法計算推進(jìn)劑表面溫度，相比于采用經(jīng)驗公式的方法，提高了推進(jìn)劑裝藥表面溫度計算的準(zhǔn)確性。

固體推進(jìn)劑內(nèi)部區(qū)域非定常熱傳導(dǎo)方程為

(5)

(6)

式中：下標(biāo)f和s分別代表流體和固體區(qū)域的物理量；qr為輻射熱流密度，文中使用的輻射傳熱模型[24]為

(7)

式中：Tw為固體推進(jìn)劑表面溫度；σ為Stefan-Boltzman常數(shù)；Csf為考慮熱傳導(dǎo)的經(jīng)驗系數(shù)，取為0.25。上述耦合傳熱計算方法驗證見文獻(xiàn)[25]。

點火準(zhǔn)則采用推進(jìn)劑表面溫度進(jìn)行判斷，當(dāng)某單元溫度達(dá)到推進(jìn)劑點火的臨界溫度時，推進(jìn)劑被點燃，該單元開始加質(zhì)，加質(zhì)源項見式(2)。推進(jìn)劑表面加質(zhì)單元的絕熱燃燒溫度達(dá)3 400 K，由于本文忽略推進(jìn)劑表面復(fù)雜的化學(xué)反應(yīng)和傳質(zhì)過程，此時基于耦合傳熱方法計算的推進(jìn)劑表面溫度不合理。實際上，推進(jìn)劑點燃后計算的推進(jìn)劑表面溫度已沒有意義，對整個Ⅱ脈沖點火過程流場沒有影響，因此，本文基于耦合傳熱方法計算的推進(jìn)劑表面溫度僅適用于未點燃的推進(jìn)劑。由于發(fā)動機點火過程持續(xù)時間較短，忽略固體推進(jìn)劑燃面的移動。

1.4 流固耦合過程

Ⅱ脈沖點火過程中金屬膜片變形運動為典型的雙向流固耦合問題，本文采用分區(qū)迭代耦合算法。該算法是指在每一個時間步對單物理場一次求解，并通過界面物理量傳遞方法交換信息，在時間步推進(jìn)過程中實現(xiàn)耦合問題的求解。依據(jù)在每個時間步推進(jìn)過程中是否進(jìn)行子迭代可以分為松耦合和緊耦合2種方法。然而緊耦合算法在反復(fù)迭代過程中由于計算量較大，導(dǎo)致計算效率降低，在某些工況下多次迭代后，仍然可能達(dá)不到收斂的效果，在一定程度上限制了其工程應(yīng)用，因此本文采用松耦合算法[26]，其示意圖如圖2所示。

結(jié)合固體推進(jìn)劑的點火過程，給定所有計算區(qū)域初始值和邊界條件，計算步驟如下：

1) 運行流體求解器得到t時刻流體區(qū)域的流場分布。

2) 基于耦合傳熱方法計算推進(jìn)劑裝藥表面溫度，如果某點溫度達(dá)到推進(jìn)劑點火的臨界值，則單元進(jìn)行加質(zhì)，否則不加質(zhì)。

3) 將流固耦合界面上的氣動壓強，轉(zhuǎn)換為固體邊界的外力載荷，以供固體位移場進(jìn)行計算。

4) 由流固耦合界面上的外力載荷，結(jié)合t時刻結(jié)構(gòu)的狀態(tài)進(jìn)行t+Δt時刻結(jié)構(gòu)動力學(xué)數(shù)值模擬。

5) 將流固耦合界面上的信息傳遞給流場，調(diào)整流體區(qū)域網(wǎng)格，計算網(wǎng)格移動速度。

6) 重復(fù)步驟1)，運行流體求解器，將時間推進(jìn)到t+Δt時刻。

流固耦合界面上物理量的插值采用梯度法，另外由于點火過程持續(xù)時間較短，金屬膜片最大運動位移不超過1 mm，因此流體區(qū)域網(wǎng)格只需根據(jù)運動后的邊界進(jìn)行調(diào)整。

圖2 松耦合算法求解時序圖
Fig.2 Generic cycle of loosely coupled algorithm

2 計算物理模型及邊界

根據(jù)Schilling等[6]的雙脈沖發(fā)動機實驗，建立Ⅱ脈沖點火過程計算區(qū)域及邊界如圖3所示，其中Ⅱ脈沖燃燒室長度為600 mm，Ⅱ脈沖裝藥采用單孔管狀藥，內(nèi)徑為60 mm，點火具出口直徑為20 mm，Ⅰ脈沖發(fā)動機長度為1 050 mm，噴管喉徑為30 mm。對上述計算區(qū)域進(jìn)行多塊結(jié)構(gòu)化網(wǎng)格劃分，流體計算區(qū)域網(wǎng)格單元總數(shù)為252 894，物面附近網(wǎng)格最小尺寸為0.001 mm，以保證壁面處的y+≤1，固體推進(jìn)劑計算區(qū)域網(wǎng)格總數(shù)為38 360，全局時間步長為1.0×10-6s。

點火具入口給定質(zhì)量流率和點火燃?xì)饪倻貫? 590 K，點火燃?xì)庋剌S線流動，且質(zhì)量流率隨時間變化，最大質(zhì)量流率為1.0 kg/s，點火持續(xù)時間為15 ms，噴管出口為壓力出口邊界條件，其余表面均為無滑移絕熱壁面，當(dāng)金屬膜片破裂后將金屬膜片壁面邊界條件改為內(nèi)部邊界。流體計算域初始值為壓力101 325 Pa、溫度300 K、速度為0，固體推進(jìn)劑初始溫度為300 K。仿真中所用燃?xì)夂屯七M(jìn)劑的物性參數(shù)見表1。

脈沖隔離裝置包括1個支撐件和1個高強度易變形的金屬膜片。當(dāng)Ⅰ脈沖工作時，金屬膜片蓋在支撐件后面，防止燃?xì)馔ㄟ^，并保證膜片不發(fā)生破壞；Ⅱ脈沖工作時，膜片破裂，燃?xì)馔ㄟ^支撐件從噴管中排出。其中金屬膜片的結(jié)構(gòu)形式主要為薄板結(jié)構(gòu)，且在一側(cè)設(shè)有“米”字型的預(yù)制缺陷槽，以方便控制膜片的破壞壓強和破壞時間，其結(jié)構(gòu)如圖4所示。金屬膜片外徑為70 mm，厚度h=3 mm，刻槽深度為1 mm，材料選擇硬鋁合金2A12，其密度為2 700 kg/m3，彈性模量為68 GPa，泊松比為0.33，強度極限為350 MPa，當(dāng)金屬膜片受到的最大應(yīng)力超過該值，金屬膜片失去承載能力，膜片破裂。

圖3 Ⅱ脈沖點火過程計算模型及邊界
Fig.3 Computational model and boundary conditions for ignition processing of second pulse

表1 仿真參數(shù)Table 1 Parameters of simulation

圖4 膜片結(jié)構(gòu)簡圖
Fig.4 Schematic diagram of burst diaphragm

Ⅱ脈沖點火過程中，金屬膜片承壓變形較小(最大位移不超過1 mm)，因此本文將二維軸對稱計算得到的壓力，繞軸旋轉(zhuǎn)后施加在膜片表面，作為壓力載荷。同時考慮到膜片結(jié)構(gòu)的對稱性，文中選取四分之一模型進(jìn)行計算，膜片外表面為固支邊界。另外，由于燃燒室內(nèi)溫度較高，膜片破裂后形成的碎片被消融，文中忽略碎片對流場的影響。

3 算例驗證

法國普羅旺斯大學(xué)的Giordano等[27]使用激波管研究沖擊載荷作用下平板的變形運動過程，實驗中拍攝到了平板運動過程陰影圖，得到了平板頂端位移隨時間變化的曲線。該實驗為典型的流固耦合問題，可以作為流固耦合計算平臺的驗證算例。

計算區(qū)域及模型參數(shù)如圖5所示，平板厚度為1 mm，高度為50 mm。初始流場壓力為101 325 Pa，溫度為293 K，馬赫數(shù)為1.21的激波從左向右沖擊豎直平板。平板的彈性模量E=220 GPa，密度ρ=2 700 kg/m3，泊松比ν=0.33。流體區(qū)域網(wǎng)格總數(shù)為106 547，固體區(qū)域為891。

圖5 激波沖擊豎直平板實驗裝置
Fig.5 Experiment setup of vertical plate exposed to shock

圖6給出了實驗紋影圖與本文計算結(jié)果對比，由圖可知，本文所發(fā)展的多物理場耦合求解器可以準(zhǔn)確計算出平板附近的激波結(jié)構(gòu)。圖7為平板頂點水平位移隨時間的變化曲線，從圖中可以看出：計算結(jié)果與實驗值吻合得較好，表明本文建立的流固耦合方法具有較高的可信度。

圖6 實驗激波結(jié)構(gòu)(上)與仿真結(jié)果(下)對比
Fig.6 Comparison of experimental shockwave structures (upper) with numerical result (lower)

圖7 平板頂點水平位移隨時間變化曲線
Fig.7 Time-history curves of panel tip horizontal displacement

4 計算結(jié)果與分析

4.1 Ⅱ脈沖點火過程內(nèi)流場特性

Ⅱ脈沖點火過程等壓線分布如圖8所示(圖中y方向尺寸放大2.5倍，下同)。由圖可見，Ⅱ脈沖點火具啟動后，點火燃?xì)馀蛎洸⑾蛉紵蚁掠我苿樱鐖D8(a)所示。與推進(jìn)劑裝藥內(nèi)表面碰撞后以正激波形式繼續(xù)向前推移，且在裝藥通道中沿軸線傳播，如圖8(b)所示。點火激波沿裝藥通道向下游推進(jìn)，波后出現(xiàn)高壓區(qū)域，形成較大的增壓速率，使Ⅱ脈沖裝藥表面各點的壓強依次驟升。裝藥通道中正激波繼續(xù)向前推移，遇到金屬膜片后形成反射激波，并向Ⅱ脈沖燃燒室上游流動，與點火藥氣體相互作用形成更加復(fù)雜的波系，如圖8(c)所示，反射激波造成推進(jìn)劑裝藥表面二次升壓。

在點火藥氣體和加質(zhì)燃?xì)夤餐饔孟拢蛎}沖燃燒室內(nèi)的壓強逐漸升高，金屬膜片破裂，氣體膨脹并進(jìn)入Ⅰ脈沖燃燒室，并以正激波繼續(xù)向前推移，如圖8(d)所示。由于Ⅱ脈沖燃燒室尾部壓強(2.98 MPa)高于Ⅰ脈沖燃燒室壓強(0.1 MPa)，在Ⅰ脈沖燃燒室內(nèi)形成管內(nèi)約束高度欠膨脹射流現(xiàn)象，如圖9所示。不斷向外擴張的膨脹波在Ⅰ脈沖燃燒室壁面反射形成入射激波，入射激波遇到馬赫盤后再次發(fā)生反射，產(chǎn)生反射激波。從圖9中可以清晰看出，在馬赫盤邊緣位置處入射激波和反射激波交匯形成三叉激波結(jié)構(gòu)。

圖8 Ⅱ脈沖點火過程等壓線分布
Fig.8 Pressure contour lines during ignition processing of second pulse

圖9 Ⅰ脈沖燃燒室內(nèi)馬赫數(shù)云圖
Fig.9 Mach number contours in first pulse combustion chamber

隨著時間的推移，壓力波以正激波形式向Ⅰ脈沖燃燒室下游移動，部分燃?xì)庠谟龅絿姽苁諗慷魏笮纬煞瓷浼げ?，并最終從噴管中排出，如圖8(e)所示。隨后燃?xì)獠粩嗵畛洧衩}沖燃燒室，燃燒室內(nèi)壓強逐漸升高，Ⅰ脈沖和Ⅱ脈沖燃燒室壓強趨于一致，燃燒室內(nèi)高度欠膨脹射流退化為弱欠膨脹射流，最終欠膨脹射流現(xiàn)象消失，燃燒室壓強逐漸趨于平緩，上述壓強驟然上升過程逐漸消失，發(fā)動機進(jìn)入穩(wěn)定工作狀態(tài)，如圖8(f)所示。

圖10給出了Ⅱ脈沖點火過程中燃燒室內(nèi)溫度演化歷程。從圖中可以看出，點火具啟動初期點火藥氣體還未接觸到裝藥內(nèi)表面，對推進(jìn)劑的加熱效果不明顯，如圖10(a)所示。隨著時間推移，點火藥氣體質(zhì)量流率逐漸增大，高溫點火藥氣體進(jìn)入Ⅱ脈沖裝藥通道，流經(jīng)推進(jìn)劑裝藥表面并將能量傳遞給推進(jìn)劑，推進(jìn)劑裝藥表面的溫度逐漸升高，如圖10(b)所示。圖11所示的是不同時刻未點燃的推進(jìn)劑裝藥表面溫度分布，由圖可見，在t=2.99 ms時刻，位于x=0.163 m處的裝藥表面溫度首先達(dá)到臨界點火溫度750 K。該位置位于點火藥氣體形成的回流區(qū)再附著點上游，根據(jù)回流區(qū)流動傳熱特性，此處熱流密度最大，從而導(dǎo)致此處裝藥表面首先被點燃。推進(jìn)劑表面點燃后，燃燒加質(zhì)產(chǎn)生高溫高壓燃?xì)?，燃燒火焰沿推進(jìn)劑壁面分別向頭部和尾部傳播，如圖10(c)所示。在t=4.0 ms和t=5.0 ms，約61.4%和89.5%推進(jìn)劑裝藥表面被點燃，最終在t=5.3 ms推進(jìn)劑表面全部點燃，如圖10(d)所示。

圖10 不同時刻燃燒室溫度分布
Fig.10 Temperature contours of combustion chamber at different times

圖11 不同時刻推進(jìn)劑表面溫度分布
Fig.11 Temperature distributions on propellant surface at different times

4.2 Ⅱ脈沖點火過程壓強沖擊特性

圖12所示的是Ⅱ脈沖推進(jìn)劑裝藥內(nèi)通道壓強建立過程，結(jié)合點火增壓過程等壓線圖可知，點火具啟動后，點火激波正向傳播，壓強逐漸升高，在約1.3 ms遇到金屬膜片形成反射激波，并迅速向Ⅱ脈沖燃燒室上游移動。壓力波沿裝藥內(nèi)通道的反射傳播過程造成推進(jìn)劑裝藥表面二次增壓。在約3.0 ms推進(jìn)劑裝藥表面開始加質(zhì)，但是由于裝藥通道容積較大，裝藥通道表面壓強并沒有迅速增加，而是與點火激波相互作用造成壓力波動。

隨著加質(zhì)燃面的增大，通道內(nèi)壓強逐漸升高，在約4.9 ms金屬膜片破裂，通道內(nèi)壓強迅速下降，且靠近Ⅱ脈沖尾部的壓強首先下降。壓力波繼續(xù)沿Ⅰ脈沖燃燒室下游移動，在約6.4 ms從噴管喉道中噴出，隨后燃燒室內(nèi)壓強逐漸升高，兩級燃燒室壓強趨于一致，發(fā)動機進(jìn)入穩(wěn)定工作狀態(tài)。

圖13為Ⅰ脈沖和Ⅱ脈沖燃燒室內(nèi)監(jiān)測點壓力隨時間的變化曲線。由圖可見，在金屬膜片破裂前，Ⅱ脈沖燃燒室內(nèi)壓強迅速上升，在t=4.9 ms時，金屬膜片破裂，由于Ⅱ脈沖燃燒室內(nèi)高壓燃?xì)獾耐蝗会尫?，壓強急劇下降，且Ⅱ脈沖燃燒室尾部壓強相比于頭部壓強，降壓幅度相對較大。

圖12 推進(jìn)劑裝藥內(nèi)通道壓強建立過程
Fig.12 Pressurization process of internal channel in propellant charge

圖13 燃燒室內(nèi)壓力隨時間變化
Fig.13 Time-histories of instantaneous pressures of combustion chamber

另外，在Ⅱ脈沖燃燒室高壓燃?xì)忉尫懦跗?，由于燃燒室?nèi)激波的傳播，激波、膨脹波及渦之間相互干擾使得各監(jiān)測點的壓強均出現(xiàn)一定程度的振蕩，且Ⅰ脈沖燃燒室頭部和尾部分別出現(xiàn)極大值。這是因為釋放的高壓燃?xì)鈮嚎sⅠ脈沖燃燒室內(nèi)低壓氣體，形成激波，且沿軸向朝噴管處傳播，激波依次掃過Ⅰ脈沖燃燒室頭部和尾部，使得各監(jiān)測點處壓強先后劇增。激波傳播至噴管收斂段后發(fā)生碰撞并反射，反射激波沿負(fù)方向朝Ⅰ脈沖燃燒室頭部傳播，重新依次掃過各監(jiān)測點，導(dǎo)致各監(jiān)測點壓強先后再次劇增，出現(xiàn)明顯的極大值，隨后激波衰減為壓縮波，在燃燒室內(nèi)來回運動造成壓強持續(xù)振蕩。隨著裝藥表面燃燒加質(zhì)，高溫高壓燃?xì)獬掷m(xù)填充Ⅰ脈沖燃燒室，燃燒室內(nèi)壓強逐漸上升到發(fā)動機穩(wěn)定工作壓強。由于金屬膜片破裂后形成級間通道，燃?xì)饬髁鹘?jīng)級間通道后壓力有所損失，因此Ⅱ脈沖燃燒室內(nèi)壓強整體高于Ⅰ脈沖燃燒室。

4.3 Ⅱ脈沖點火過程中金屬膜片力學(xué)特性

由Ⅱ脈沖點火過程流場分析可知，點火過程的瞬態(tài)沖擊作用對金屬膜片表面產(chǎn)生瞬變的沖擊載荷，金屬膜片開始變形。圖14中所示是 3.0 ms 和4.9 ms時刻金屬膜片Mises等效應(yīng)力分布，由圖可見，由于金屬膜片一側(cè)存在“米”字型凹槽，因此中心刻槽處金屬膜片的應(yīng)力最大。

圖14 不同時刻膜片von Mises等效應(yīng)力分布
Fig.14 von Mises stress distributions of diaphragm at different times

圖15 金屬膜片最大軸向位移和壓強隨時間變化
Fig.15 Variation of the maximum horizontal displacement and pressure of diaphragm with time

圖15中給出了金屬膜片最大軸向位移和壓強隨時間的變化曲線，從圖中可以看出，在約1.3 ms 之前，Ⅱ脈沖點火初期，壓力峰尚未移動至金屬膜片處，或者作用在金屬膜片上的壓強較小，金屬膜片軸向位移較??；而當(dāng)時間大于1.3 ms 后，作用在金屬膜片表面的壓力迅速增大，金屬膜片開始變形，軸向位移不斷增大，3.0 ms時Mises等效應(yīng)力最大值為44 MPa，4.9 ms時金屬膜片Mises等效應(yīng)力最大值為350 MPa，達(dá)到金屬膜片材料的強度極限，金屬膜片失去承壓能力，膜片破裂。另外，膜片軸向位移隨時間變化與作用在膜片上的壓強變化趨勢一致，且在壓力載荷作用下膜片出現(xiàn)了小幅振動現(xiàn)象。

4.4 點火具質(zhì)量流率對點火過程影響

圖16和圖17所示為點火質(zhì)量流率對Ⅱ脈沖燃燒室尾部壓強和金屬膜片軸向位移的影響。由圖可見，不同點火質(zhì)量流率下，燃燒室內(nèi)壓強和金屬膜片軸向位移變化趨勢相同，質(zhì)量流率越大，燃燒室壓力越大，膜片破裂時間越短。這是因為在金屬膜片表面壓強載荷最終上升階段，不同質(zhì)量流率下的平均升壓速率分別為1.80、1.90、2.54、2.61 GPa/s，增壓速率越大，金屬膜片的破裂時間越短。因此可以得到以下結(jié)論，金屬膜片破裂時間不僅與作用在其表面的壓強載荷大小相關(guān)，而且與壓強載荷加載的過程相關(guān)，采用靜態(tài)分析方法并不能準(zhǔn)確得到金屬膜片破裂壓強和時間，這與文獻(xiàn)[13]中的實驗結(jié)果相吻合。

表2 不同點火質(zhì)量流率下Ⅱ脈沖點火特征參數(shù)

圖16 點火質(zhì)量流率對Ⅱ脈沖燃燒室尾部壓強的影響
Fig.16 Effects of ignition mass flow rates on aft-end pressure of second pulse combustion chamber

圖17 點火質(zhì)量流率對金屬膜片軸向位移的影響
Fig.17 Effects of ignition mass flow rates on the maximum horizontal displacement of diaphragm

4.5 金屬膜片厚度對點火過程影響

不同金屬膜片厚度(h)下Ⅱ脈沖點火過程特征參數(shù)計算結(jié)果如表3所示。由計算結(jié)果可見，金屬膜片厚度對推進(jìn)劑首次點燃和全部點燃時間影響可以忽略；金屬膜片越薄，膜片有效厚度越小，膜片破裂時間越短，膜片破裂壓強越小。相比于點火質(zhì)量流率，金屬膜片厚度對膜片破裂時間和壓強影響更為明顯。

圖18給出了金屬膜片厚度對Ⅱ脈沖燃燒室尾部壓強和金屬膜片軸向位移的影響。由圖可見，不同膜片厚度下，燃燒室內(nèi)壓強和金屬膜片軸向位移變化趨勢相同。但是，隨著膜片厚度的減小，膜片有效厚度降低，膜片承壓能力減弱，膜片破裂時間越短，最大軸向位移越大。膜片厚度為2.0 mm時，金屬膜片的最大軸向位移達(dá)到0.6 mm。

圖19給出了金屬膜片厚度為3.5 mm時Ⅱ脈沖發(fā)動機頭部壓力隨時間變化曲線與文獻(xiàn)[6]中實驗值的對比，由圖可見，計算得到的膜片破裂時的壓強比實驗值略大，破裂時間較早，泄壓過程持續(xù)時間較短，泄壓幅度較大。這是因為文獻(xiàn)[6]中裝藥結(jié)構(gòu)和性能、金屬膜片厚度、脈沖隔離裝置通道孔徑等參數(shù)并未明確給出，本文計算參數(shù)與文獻(xiàn)中的實驗有一定的差異。但是，計算結(jié)果與實驗值總體趨勢是一致的，表明本文多物理場耦合求解器的能夠應(yīng)用于金屬膜片式雙脈沖發(fā)動機Ⅱ脈沖點火過程的研究。

表3 不同金屬膜片厚度下Ⅱ脈沖點火特征參數(shù)

圖18 金屬膜片厚度對Ⅱ脈沖燃燒室尾部壓強和金屬膜片軸向位移的影響
Fig.18 Effects of diaphragm thickness on aft-end pressure of second pulse combustion chamber and maximum horizontal displacement of diaphragm

圖19 Ⅱ脈沖發(fā)動機頭部壓力計算結(jié)果與實驗值對比
Fig.19 Comparison of head-end pressure of second pulse motor present results with experimental values

5 結(jié) 論

1) 通過經(jīng)典算例驗證可知本文發(fā)展的多物理場耦合求解器可信度較高，建立的流固耦合方法和流熱耦合方法是正確的；該求解器能夠?qū)﹄p脈沖發(fā)動機Ⅱ脈沖點火過程進(jìn)行數(shù)值模擬。

2) 隨著點火質(zhì)量流率增加，推進(jìn)劑首次點燃時間和推進(jìn)劑全部點燃時間變短，燃?xì)庠谘b藥通道內(nèi)流速增加，膜片破裂時間變短，膜片破裂壓強降低；金屬膜片破裂時間和壓強不僅與作用在其表面的壓強載荷大小相關(guān)，而且與壓強載荷加載的過程相關(guān)，采用靜態(tài)分析方法并不能準(zhǔn)確得到金屬膜片破裂壓強和時間。

3) 隨著金屬膜片厚度降低，膜片破裂時間變短，膜片軸向位移增大，膜片破裂壓強降低，且相對于點火質(zhì)量流率，金屬膜片厚度對膜片破裂壓強和時間的影響更為明顯。

以上研究對金屬膜片式雙脈沖發(fā)動機設(shè)計具有一定的參考意義，同時需基于擴展有限元方法展開金屬膜片破裂過程的數(shù)值模擬，并考慮多種裝藥類型，如星孔裝藥、輪轂式裝藥，將本文的流場求解器擴展至三維。

[1] NAUMANN K W, STADLER L. Double pulse solid rocket motor technology applications and technical solutions: AIAA-2010-6754[R]. Reston: AIAA, 2010.

[2] NISHII S, FUKUDA K, KUBOTA N. Combustion tests of two stage pulse rocket motors: AIAA-1989-2426[R]. Reston: AIAA, 1989.

[3] CARRIER J L C, CONSTANTINOU T, HARRIS P G, et al. Dual-interrupted-thrust pulse motor[J]. Journal of Propulsion and Power, 1987, 3(4): 308-312.

[4] DAHL H, JONES B. Demonstration of solid propellant pulse motor technologies: AIAA-1996-3157[R]. Reston: AIAA,1996.

[5] WANG C H, LIU Y, LIU Y B. Design and experimental studies on ceramic port cover for dual pulse motor[J]. Acta Astronautica, 2011, 68(11): 1881-1890.

[6] SCHILLING S, TROUILLOT P, WEIGAND A. On the development and testing of a 120 mm caliber double pulse motor (DPM): AIAA-2004-3387[R]. Reston: AIAA, 2004.

[7] STADLER L J, HOFFMANN S, HUBER J, et al. The flight demonstration of the double pulse motor demonstrator MSA: AIAA-2010-6756[R]. Reston: AIAA, 2010.

[8] JAVED A, MANNA P， DEBASIS C. Numerical simulation of a dual pulse solid rocket motor flow field[J]. Defence Science Journal, 2012, 62(6): 369-374.

[9] 孫娜, 婁永春, 孫長宏, 等. 某雙脈沖發(fā)動機燃燒室兩相流場數(shù)值分析[J]. 固體火箭技術(shù), 2012, 35(3): 335-338.

SUN N, LOU Y C, SUN C H, et al. Numerical analysis of two-phase flow in combustion chamber of dual-pulse motor[J]. Journal of Solid Rocket Technology, 2012, 35(3): 335-338 (in Chinese).

[10] 王春光, 任全彬, 田維平, 等. 脈沖發(fā)動機中金屬膜片式隔艙動態(tài)破壞過程研究[J]. 固體火箭技術(shù), 2013, 36(1): 22-26.

WANG C G, REN Q B, TIAN W P, et al. Research on the process of dynamic failure of metal diaphragm pulse separation device in pulse motor[J]. Journal of Solid Rocket Technology, 2013, 36(1): 22-26 (in Chinese).

[11] 石瑞, 王長輝, 萇艷楠. 雙脈沖固體發(fā)動機鋁膜隔板設(shè)計和試驗研究[J]. 固體火箭技術(shù), 2013, 36(2): 190-194.

SHI R, WANG C H, CHANG Y N. Design and experimental study on aluminum clapboard of dual-pulse rocket motor[J]. Journal of Solid Rocket Technology, 2013, 36(2): 190-194 (in Chinese).

[12] 劉偉凱, 惠博. 雙脈沖發(fā)動機中金屬膜片式隔艙設(shè)計方法[J]. 固體火箭技術(shù), 2013, 36(4): 486-491.

LIU W K, HUI B. Research on designing method of metal diaphragm PSD in double pulse solid rocket motor[J]. Journal of Solid Rocket Technology, 2013, 36(4): 486-491 (in Chinese).

[13] 劉偉凱, 何國強, 王春光. 雙脈沖發(fā)動機中金屬膜片動態(tài)與靜態(tài)打開對比分析[J]. 推進(jìn)技術(shù), 2014, 35(9): 1259-1264.

LIU W K, HE G Q, WANG C G. Research on dynamic and static opening of metal diaphragm in double pulse solid rocket motor[J]. Journal of Propulsion Technology, 2014, 35(9): 1259-1264 (in Chinese).

[14] 王偉, 李江, 王春光, 等. 隔艙式雙脈沖發(fā)動機金屬膜片設(shè)計與實驗研究[J]. 推進(jìn)技術(shù), 2013, 34(8): 1115-1120.

WANG W, LI J, WANG C G, et al. Study on metal diaphragm of pulse separation device in dual pulse solid rocket motor[J]. Journal of Propulsion Technology, 2013, 34(8): 1115-1120 (in Chinese).

[15] 李映坤, 韓珺禮， 陳雄， 等. 級間通道構(gòu)型對雙脈沖發(fā)動機燃燒室局部受熱的影響[J]. 推進(jìn)技術(shù)， 2014， 35(11)： 1503-1510.

LI Y K， HAN J L， CHEN X, et al. Effects of two-stage pulse channel configurations on local heat transfer characteristics in combustion chamber of dual pulse motor[J]. Journal of Propulsion Technology, 2014, 35(11): 1503-1510 (in Chinese).

[16] 陳雄, 李映坤, 劉銳, 等. 基于耦合傳熱的雙脈沖發(fā)動機熱防護(hù)層受熱分析[J]. 推進(jìn)技術(shù), 2016, 37(1): 83-89.

CHEN X, LI Y K, LIU R, et al. Heating study of thermal protection layer in dual pulse motor based on the conjugate heat transfer method[J]. Journal of Propulsion Technology, 2016, 37(1): 83-89 (in Chinese).

[17] 李映坤, 韓珺禮, 陳雄, 等. 基于嵌套網(wǎng)格的脈沖發(fā)動機噴管內(nèi)流場數(shù)值模擬[J]. 固體火箭技術(shù), 2014, 37(2): 178-183.

LI Y K, HAN J L, CHEN X, et al. Numerical simulation research of nozzle inner flow field for pulse motor based on dynamic chimera grid[J]. Journal of Solid Rocket Technology, 2014, 37(2): 178-183 (in Chinese).

[18] 鞏倫昆, 陳雄, 周長省, 等. 來流條件對固體燃料沖壓發(fā)動機燃速及自持燃燒性能影響的仿真研究[J]. 航空學(xué)報, 2016, 37(5): 1428-1439.

GONG L K, CHEN X, ZHOU C S, et al. Numerical investigation on the effect of inlet flow condition on regression rate and self-sustained combustion of solid fuel ramjet[J]. Acta Aeronautica et Astronautica Sinica, 2016, 37(5):1428-1439 (in Chinese).

[19] MENTER F R. Two equation eddy viscosity turbulence models for engineering application[J]. AIAA Journal, 1994, 32(8): 1598-1605．

[20] KIM K H, KIM C, RHO O H. Methods for the accurate computations of hypersonic flows: I. AUSMPW+ scheme[J]. Journal of Computational Physics, 2001, 174(1): 38-80.

[21] 閻超. 計算流體力學(xué)方法及其應(yīng)用[M]. 北京： 北京航空航天大學(xué)出版社, 2006: 160-164.

YAN C. Computational fluid dynamics method and its application[M]. Beijing: Beihang University Press, 2006: 160-164 (in Chinese).

[22] SMITH I M, GRIFFITHS D V. Programming the finite element method[M]. 4th ed. New York: John Wiley & Sons, 2005: 465-480.

[23] LI Q, LIU P, HE G Q. Fluid-solid coupled simulation of the ignition transient of solid rocket motor[J]. Acta Astronautica, 2015, 110: 180-190.

[24] JOHNSTON W A. Solid rocket motor internal flow during ignition[J]. Journal of Propulsion and Power, 1995, 11(3): 489-496.

[25] LIU R, CHEN X, ZHOU C S, et al. A couple approach for a conjugate heat transfer investigation of the shape-change effects in a composite nozzle[J]. Numerical Heat Transfer, Part A: Applications, 2015, 68(11): 1280-1305.

[26] 李偉, 馬寶峰. 一種改進(jìn)型松耦合方法在機翼搖滾計算中的應(yīng)用[J]. 航空學(xué)報, 2015, 36(6): 1805-1813.

LI W, MA B F. A modified loosely-coupled algorithm for calculation of wing rock[J]. Acta Aeronautica et Astronautica Sinica, 2015, 36(6): 1805-1813 (in Chinese).

[27] GIORDANO J, JOURDAN G, BURTSCHELL Y, et al. Shock wave impacts on deforming panel, an application of fluid-structure interaction[J]. Shock Waves, 2005, 14(1-2): 103-110.

Numericalsimulationoftheignitiontransientofdualpulsemotorbasedonmulti-physicscoupling

LIYingkun1,HANJunli1，2,CHENXiong1,*,ZHOUChangsheng1,GONGLunkun1

1.SchoolofMechanicalEngineering,NanjingUniversityofScienceandTechnology,Nanjing210094,China2.BeijingInstituteofElectromechanicalTechnology,Beijing100083,China

Inordertostudythesecondpulseignitiontransientofadualpulsesolidrocketmotor,amulti-physicsolverisdeveloped.ThegoverningequationsforunsteadycompressiblefluidflowaresolvedwithdualtimeLU-SGS(loweruppersymmetricGuassseidel)iterativealgorithmbyfinitevolumemethod.Theconjugateheattransferstrategyisemployedtocalculatethepropellantsurfacetemperature.Afiniteelementmethodisusedtodiscretizethestructuraldynamicequationinspace,whereasthetemporaltimeintegrationisachievedwiththeclassicNewmarkalgorithm.Alooselycoupledalgorithmisusedforfluidstructureinteractionproblems,andthereliabilityofthenumericalapproachisvalidatedbyacomparisonwithexperimentalcases.Resultsshowthatthemulti-physicssolvercansimulatetheimpactofignitiongas,strongunsteadyflow,andmechanicalresponseofmetaldiaphragm.Thebursttimeandburstpressureofmetaldiaphragmcanbealsoacquired.Meanwhile,withtheincreaseoftheignitionmassflowrate,thefirstignitiontimeofpropellantandthebursttimeofthediaphragmbecomeshorterandtheburstpressureofthediagramdecreases.Thebursttimeandburstpressureofmetaldiaphragmarenotonlyrelatedtothepressureloadonthesurfaceofdiaphragm,butalsotothehistoryofthepressureloadonit.Withthedecreaseofthicknessofmetaldiaphragm,thebursttimeofthediaphragmgoesshorter,theburstpressureofdiaphragmdecreases,andthemaximumhorizontaldisplacementofthediaphragmincreases.

multi-physicscoupling;fluidstructureinteraction;conjugateheattransfer;ignition;dualpulsemotor;solidrocketmotor;numericalsimulation

2016-05-09；Revised2016-06-20；Accepted2016-07-17；Publishedonline2016-08-011034

URL：www.cnki.net/kcms/detail/11.1929.V.20160801.1034.002.html

TheResearchInnovationProgramforCollegeAcademicGraduatesofJiangsuProvince(KYZZ15_0113)

2016-05-09；退修日期2016-06-20；錄用日期2016-07-17； < class="emphasis_bold">網(wǎng)絡(luò)出版時間

時間：2016-08-011034

www.cnki.net/kcms/detail/11.1929.V.20160801.1034.002.html

江蘇省普通高校學(xué)術(shù)學(xué)位研究生科研創(chuàng)新計劃 (KYZZ15_0113)

＊

.E-mailchenxiong@njust.edu.cn

李映坤， 韓珺禮， 陳雄， 等． 基于多物理場耦合的雙脈沖發(fā)動機點火過程數(shù)值模擬J． 航空學(xué)報,2017,38(4):120409.LIYK,HANJL,CHENX,etal.Numericalsimulationoftheignitiontransientofdualpulsemotorbasedonmulti-physicscouplingJ.ActaAeronauticaetAstronauticaSinica,2017,38(4):120409.

http://hkxb.buaa.edu.cnhkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0212

V235

A

1000-6893(2017)04-120409-12

(責(zé)任編輯： 彭健， 張晗)

＊Correspondingauthor.E-mailchenxiong@njust.edu.cn