基于多物理場耦合的雙脈沖發動機點火過程數值模擬

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

1.南京理工大學 機械工程學院， 南京 210094

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

基于多物理場耦合的雙脈沖發動機點火過程數值模擬

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

1.南京理工大學 機械工程學院， 南京 210094

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

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

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

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

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

圖1 雙脈沖發動機結構簡圖
Fig.1 Schematic diagram of dual pulse motor

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

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

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

本文發展一套多物理場耦合求解器，對Ⅱ脈沖啟動過程中點火藥氣體沖擊過程、推進劑裝藥受熱點火過程、燃氣的非定常流動及金屬膜片機械響應過程進行數值模擬，獲得金屬膜片破裂時間和壓強，分析點火具質量流率和金屬膜片厚度對發動機點火瞬態過程的影響。

1 控制方程和計算方法

1.1 流體區域

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

?ΩHdV+?ΩSdV

(1)

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

(2)

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

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

1.2 固體結構區域

基于有限元方法的基本原理和基本步驟，結構動力學運動方程可以寫為

(3)

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

結構動力學運動方程的求解采用經典的Newmark 方法[22]。給定時間間隔[t,t+Δt]，已知t時刻的q(t)及其導數，t+Δt時刻的位移、速度和加速度由下式計算：

(4)

式中：α和β為兩常數，下標t和t+Δt分別代表時刻t、t+Δt的值。

1.3 固體推進劑點火模型

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

固體推進劑內部區域非定常熱傳導方程為

(5)

(6)

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

(7)

式中：Tw為固體推進劑表面溫度；σ為Stefan-Boltzman常數；Csf為考慮熱傳導的經驗系數，取為0.25。上述耦合傳熱計算方法驗證見文獻[25]。

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

1.4 流固耦合過程

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

結合固體推進劑的點火過程，給定所有計算區域初始值和邊界條件，計算步驟如下：

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

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

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

4) 由流固耦合界面上的外力載荷，結合t時刻結構的狀態進行t+Δt時刻結構動力學數值模擬。

5) 將流固耦合界面上的信息傳遞給流場，調整流體區域網格，計算網格移動速度。

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

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

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

2 計算物理模型及邊界

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

點火具入口給定質量流率和點火燃氣總溫為2 590 K，點火燃氣沿軸線流動，且質量流率隨時間變化，最大質量流率為1.0 kg/s，點火持續時間為15 ms，噴管出口為壓力出口邊界條件，其余表面均為無滑移絕熱壁面，當金屬膜片破裂后將金屬膜片壁面邊界條件改為內部邊界。流體計算域初始值為壓力101 325 Pa、溫度300 K、速度為0，固體推進劑初始溫度為300 K。仿真中所用燃氣和推進劑的物性參數見表1。

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

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

表1 仿真參數Table 1 Parameters of simulation

圖4 膜片結構簡圖
Fig.4 Schematic diagram of burst diaphragm

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

3 算例驗證

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

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

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

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

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

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

4 計算結果與分析

4.1 Ⅱ脈沖點火過程內流場特性

Ⅱ脈沖點火過程等壓線分布如圖8所示(圖中y方向尺寸放大2.5倍，下同)。由圖可見，Ⅱ脈沖點火具啟動后，點火燃氣膨脹并向燃燒室下游移動，如圖8(a)所示。與推進劑裝藥內表面碰撞后以正激波形式繼續向前推移，且在裝藥通道中沿軸線傳播，如圖8(b)所示。點火激波沿裝藥通道向下游推進，波后出現高壓區域，形成較大的增壓速率，使Ⅱ脈沖裝藥表面各點的壓強依次驟升。裝藥通道中正激波繼續向前推移，遇到金屬膜片后形成反射激波，并向Ⅱ脈沖燃燒室上游流動，與點火藥氣體相互作用形成更加復雜的波系，如圖8(c)所示，反射激波造成推進劑裝藥表面二次升壓。

在點火藥氣體和加質燃氣共同作用下，Ⅱ脈沖燃燒室內的壓強逐漸升高，金屬膜片破裂，氣體膨脹并進入Ⅰ脈沖燃燒室，并以正激波繼續向前推移，如圖8(d)所示。由于Ⅱ脈沖燃燒室尾部壓強(2.98 MPa)高于Ⅰ脈沖燃燒室壓強(0.1 MPa)，在Ⅰ脈沖燃燒室內形成管內約束高度欠膨脹射流現象，如圖9所示。不斷向外擴張的膨脹波在Ⅰ脈沖燃燒室壁面反射形成入射激波，入射激波遇到馬赫盤后再次發生反射，產生反射激波。從圖9中可以清晰看出，在馬赫盤邊緣位置處入射激波和反射激波交匯形成三叉激波結構。

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

圖9 Ⅰ脈沖燃燒室內馬赫數云圖
Fig.9 Mach number contours in first pulse combustion chamber

隨著時間的推移，壓力波以正激波形式向Ⅰ脈沖燃燒室下游移動，部分燃氣在遇到噴管收斂段后形成反射激波，并最終從噴管中排出，如圖8(e)所示。隨后燃氣不斷填充Ⅰ脈沖燃燒室，燃燒室內壓強逐漸升高，Ⅰ脈沖和Ⅱ脈沖燃燒室壓強趨于一致，燃燒室內高度欠膨脹射流退化為弱欠膨脹射流，最終欠膨脹射流現象消失，燃燒室壓強逐漸趨于平緩，上述壓強驟然上升過程逐漸消失，發動機進入穩定工作狀態，如圖8(f)所示。

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

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

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

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

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

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

圖13為Ⅰ脈沖和Ⅱ脈沖燃燒室內監測點壓力隨時間的變化曲線。由圖可見，在金屬膜片破裂前，Ⅱ脈沖燃燒室內壓強迅速上升，在t=4.9 ms時，金屬膜片破裂，由于Ⅱ脈沖燃燒室內高壓燃氣的突然釋放，壓強急劇下降，且Ⅱ脈沖燃燒室尾部壓強相比于頭部壓強，降壓幅度相對較大。

圖12 推進劑裝藥內通道壓強建立過程
Fig.12 Pressurization process of internal channel in propellant charge

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

另外，在Ⅱ脈沖燃燒室高壓燃氣釋放初期，由于燃燒室內激波的傳播，激波、膨脹波及渦之間相互干擾使得各監測點的壓強均出現一定程度的振蕩，且Ⅰ脈沖燃燒室頭部和尾部分別出現極大值。這是因為釋放的高壓燃氣壓縮Ⅰ脈沖燃燒室內低壓氣體，形成激波，且沿軸向朝噴管處傳播，激波依次掃過Ⅰ脈沖燃燒室頭部和尾部，使得各監測點處壓強先后劇增。激波傳播至噴管收斂段后發生碰撞并反射，反射激波沿負方向朝Ⅰ脈沖燃燒室頭部傳播，重新依次掃過各監測點，導致各監測點壓強先后再次劇增，出現明顯的極大值，隨后激波衰減為壓縮波，在燃燒室內來回運動造成壓強持續振蕩。隨著裝藥表面燃燒加質，高溫高壓燃氣持續填充Ⅰ脈沖燃燒室，燃燒室內壓強逐漸上升到發動機穩定工作壓強。由于金屬膜片破裂后形成級間通道，燃氣流流經級間通道后壓力有所損失，因此Ⅱ脈沖燃燒室內壓強整體高于Ⅰ脈沖燃燒室。

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

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

圖14 不同時刻膜片von Mises等效應力分布
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 之前，Ⅱ脈沖點火初期，壓力峰尚未移動至金屬膜片處，或者作用在金屬膜片上的壓強較小，金屬膜片軸向位移較小；而當時間大于1.3 ms 后，作用在金屬膜片表面的壓力迅速增大，金屬膜片開始變形，軸向位移不斷增大，3.0 ms時Mises等效應力最大值為44 MPa，4.9 ms時金屬膜片Mises等效應力最大值為350 MPa，達到金屬膜片材料的強度極限，金屬膜片失去承壓能力，膜片破裂。另外，膜片軸向位移隨時間變化與作用在膜片上的壓強變化趨勢一致，且在壓力載荷作用下膜片出現了小幅振動現象。

4.4 點火具質量流率對點火過程影響

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

表2 不同點火質量流率下Ⅱ脈沖點火特征參數

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

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

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

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

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

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

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

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

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

5 結 論

1) 通過經典算例驗證可知本文發展的多物理場耦合求解器可信度較高，建立的流固耦合方法和流熱耦合方法是正確的；該求解器能夠對雙脈沖發動機Ⅱ脈沖點火過程進行數值模擬。

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

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

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

[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] 孫娜, 婁永春, 孫長宏, 等. 某雙脈沖發動機燃燒室兩相流場數值分析[J]. 固體火箭技術, 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] 王春光, 任全彬, 田維平, 等. 脈沖發動機中金屬膜片式隔艙動態破壞過程研究[J]. 固體火箭技術, 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] 石瑞, 王長輝, 萇艷楠. 雙脈沖固體發動機鋁膜隔板設計和試驗研究[J]. 固體火箭技術, 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] 劉偉凱, 惠博. 雙脈沖發動機中金屬膜片式隔艙設計方法[J]. 固體火箭技術, 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] 劉偉凱, 何國強, 王春光. 雙脈沖發動機中金屬膜片動態與靜態打開對比分析[J]. 推進技術, 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] 王偉, 李江, 王春光, 等. 隔艙式雙脈沖發動機金屬膜片設計與實驗研究[J]. 推進技術, 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] 李映坤, 韓珺禮， 陳雄， 等. 級間通道構型對雙脈沖發動機燃燒室局部受熱的影響[J]. 推進技術， 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] 陳雄, 李映坤, 劉銳, 等. 基于耦合傳熱的雙脈沖發動機熱防護層受熱分析[J]. 推進技術, 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] 李映坤, 韓珺禮, 陳雄, 等. 基于嵌套網格的脈沖發動機噴管內流場數值模擬[J]. 固體火箭技術, 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] 鞏倫昆, 陳雄, 周長省, 等. 來流條件對固體燃料沖壓發動機燃速及自持燃燒性能影響的仿真研究[J]. 航空學報, 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] 閻超. 計算流體力學方法及其應用[M]. 北京： 北京航空航天大學出版社, 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]. 航空學報, 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">網絡出版時間

時間：2016-08-011034

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

江蘇省普通高校學術學位研究生科研創新計劃 (KYZZ15_0113)

＊

.E-mailchenxiong@njust.edu.cn

李映坤， 韓珺禮， 陳雄， 等． 基于多物理場耦合的雙脈沖發動機點火過程數值模擬J． 航空學報,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

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

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