基于Simulink的延伸噴管燃氣展開過程聯合仿真*

陳 鵬，宋學宇，曹濤鋒，尤軍峰

(1.中國航天科技集團公司四院四十一所，西安 710025；2.固體火箭發動機燃燒、熱結構與內流場國防科技重點實驗室，西安 710025；3.西北工業大學 航天學院，西安 710072)

0 引言

延伸噴管技術是提高固體火箭發動機性能的關鍵技術之一。對于采用作動筒推動方式展開的延伸噴管，作動筒展開力是決定延伸噴管動力學特性的關鍵因素。作動筒的展開方式包括燃氣式、氣瓶式、產氣式等。

采用燃氣發生器提供展開力的展開方式屬于產氣式，通過燃氣發生器燃燒室內固體藥柱的燃燒，產生大量高溫高壓的氣體，經濾網與管路進入作動筒內，推動延伸噴管展開。整個點火展開過程包括點火-燃燒-傳遞-做功-展開五個步驟，且氣體在作動筒內的展開做功過程(下游)對燃燒室內藥柱燃燒(上游)存在影響，為了對延伸噴管展開進行準確預示，必須要對整個點火展開過程進行完整建模分析。

現有的仿真計算局限于單系統的建模仿真，如燃氣發生器的內彈道仿真計算，氣體充壓過程仿真計算[5]，延伸噴管展開動力學、展開碰撞、燃氣尾流仿真計算[6-11]等。但由于燃氣發生器內彈道參數與下游過程相關，且受延伸噴管展開系統機構的影響，下游氣體做功過程無法通過簡單試驗模擬，必須聯合整個點火展開系統進行試驗，周期較長且成本較高。在噴管外載荷不確定，且展開機構受加工工藝影響較大的條件下，僅通過單個系統的仿真計算，無法給出藥柱與展開動力學參數的關系。而如何確定燃氣發生器藥形與藥量，是在工程實際中必須予以解決的問題。

針對作動筒式雙級延伸噴管采用燃氣發生器展開的工況，本文通過在Matlab/Simulink中構建燃氣發生器藥柱燃燒與作動筒內燃氣的膨脹做功模型，結合Ansys的動力學仿真模塊，對從燃氣發生器點火開始的延伸噴管展開全過程進行聯合仿真，以期在燃氣發生器藥型與藥量發生變化時，可以對延伸噴管展開過程進行預測。通過與已有試驗測試數據的對比，驗證仿真結果的準確性，建立完整可靠的延伸噴管點火展開模型。

1 燃氣展開模型

在延伸噴管的展開過程中，作動筒的充氣過程從燃氣發生器點火開始。不計點火時間，主裝藥開始燃燒并產生高焓值的氣體，經過濾網后進入管路，最終進入作動筒并膨脹做功，推動作動筒展開，帶動延伸噴管展開到位。通過建立藥柱燃燒及燃氣膨脹做功的模型，可以為延伸噴管動力學模型提供輸入條件，再經動力學仿真，最終得到延伸噴管展開的動力學參數。

1.1 燃燒過程

由維耶里模型可知，燃速rc與燃燒室內壓pc的關系為

rc=a0pcn

(1)

由于藥柱采用內孔+端面燃燒的方式，且由多個藥片組合成，單個藥片形狀如圖1所示。

圖1 藥片尺寸圖Fig.1 Size of the charge pellet

故燃面Ab為

Ab=2π[D2-(d+2rct)2]/4+

π(d+2rct)(h-2rct)

(2)

(3)

由質量守恒可知，藥柱燃燒減少的質量，一部分成為殘渣被濾網阻攔，一部分作為燃氣進入作動筒內膨脹做功，則有

(4)

受過濾網與氣體管路影響，燃氣發生器燃燒室內壓與作動筒內壓不一致，但兩者間存在一定的關系。一般來說，可建立起作動筒內壓pzdt(因變量)相對于燃氣發生器內壓pc、時間t(自變量)的單值函數關系。但由于試驗條件限制，燃氣發生器內壓無法測量，只能由試驗數據近似建立燃氣發生器內壓pc(因變量)相對于作動筒內壓pzdt(自變量)的函數關系如下：

pc=k1pzdtek2pzdt

(5)

注意到此關系內不含時間項t，且對于同一個燃燒室內壓pc可能對應不同的作動筒內壓pzdt。

1.2 膨脹做功過程

理想氣體方程:

(6)

結合質量守恒方程，有

(7)

開口系統能量守恒方程:

(8)

由于每個燃氣發生器為2個作動筒供氣，有

(9)

1.3 動力學展開過程

延伸噴管展開的動力直接來源于作動筒內的高壓氣體。因此，對于延伸噴管，作動筒內壓產生的推力是運動輸入，各運動部件的位移、速度、加速度等動力學參數是運動輸出。

展開阻力來源有三部分：(1)作動筒展開內部阻力，由預緊力產生的摩擦力與展開過程中筒間正壓力產生的摩擦力兩部分組成；(2)展開即將到位時產生的阻力，由延伸錐側向密封圈摩擦力、延伸錐端面密封圈阻力、折轉片到位后的拉力三部分組成；(3)熱試中尾部燃氣對延伸錐產生的作用力，根據燃氣不同的膨脹狀態，可能產生引射力或者阻力，其值可由流場計算得到。

由展開過程可知，作動筒容積Vzdt，與作動筒上端沿作動筒展開方向的位移szdt關系為

Vzdt=V0+Azdtszdt

(10)

由運動關系可知：

(11)

在噴管展開方向(軸向)上，作動筒內壓(考慮作動筒內摩擦fzdt)的分量克服軸向摩擦力ff、尾部燃氣阻力Fe使噴管展開。由于展開機構限制，延伸噴管I(下標1)、II(下標2)級之間的位移關系固定為1∶2，故在軸向有

(pzdtAzdt-fzdt)szdt=(ff1/2+ff2+Fe1/2+Fe2)sysz+Ek

(12)

式中sysz為II級延伸錐位移；Ek為動能項，展開到位后動量為0。

不考慮碰撞損失，延伸錐展開到位后，作動筒內燃氣做功完全用來克服摩擦力及尾部燃氣阻力。

文獻中對延伸噴管的動力學展開進行了理論分析。但針對雙級延伸噴管，在不大量簡化模型的情況下，直接通過微分方程進行理論計算是非常復雜的。為了更加詳細的反應模型參數，采用多體動力學仿真軟件Ansys對延伸噴管展開模型進行建模分析。

2 Ansys多體動力學模型

在Ansys中對1/4延伸噴管進行建模，給定作動筒上的摩擦阻力、到位等效阻力以及燃氣對延伸錐的作用力，以作動筒內壓作為輸入量，輸出作動筒上端在作動筒展開方向上的的速度和加速度。三維模型如圖2所示。

圖2 延伸噴管1/4模型Fig.2 Quarter model of the extendible nozzle

在模型中，作動筒阻力通過設置作動筒上的預緊力和摩擦系數給出；延伸錐展開到位阻力的產生包含了與密封圈摩擦和碰撞過程，模型中采用近似曲線進行模擬，圖3中縱軸為延伸錐展開到位阻力(等效在II級延伸錐上)，橫軸為II級延伸錐質心軸向坐標。亦可采用顯式動力學計算出阻力隨展開位移變化曲線，作為參數輸入模型中；尾部燃氣產生的阻力與發動機工作狀態以及延伸噴管展開位置有關，在發動機工作狀態穩定的情況下，采用流場計算得到尾部燃氣阻力隨延伸錐展開位移的變化曲線，并作為參數輸入模型中。圖4中縱軸為兩級延伸錐所受的燃氣阻力，橫軸為II級延伸錐質心軸向坐標。

圖3 延伸錐展開到位阻力曲線Fig.3 Mechanical resistance of the extendible nozzle

圖4 噴管燃氣阻力曲線Fig.4 Gas resistance of the extendible nozzle

為模擬棘輪式折轉片的鎖緊過程，防止仿真中展開到位后延伸錐回彈，在折轉片、插銷上定義接觸，并減小計算步長。仿真中發現，棘輪式轉動鎖緊機構可以相當可靠地保證鎖緊，即使未展開至預定的最終位置，也可防止延伸錐在外力作用下回彈。

在動力學模型建立完成后，通過輸入實測的作動筒內壓曲線，可較為準確地仿真出延伸錐展開曲線，如圖5所示；若給出期望的展開位移曲線，也可計算出需求的作動筒內壓曲線。即作動筒內壓曲線與延伸錐展開曲線之間的動力學模型是完整可靠的。但如果燃氣發生器藥量發生改變，很難通過對內壓曲線的線性變化來準確預示新的展開曲線，必須結合燃氣發生器內部的燃燒與燃氣膨脹做功的過程進行求解，進行多系統的聯合仿真。

圖5 展開速度仿真與實測對比Fig.5 Comparison between simulation and test results of velocity

3 Simulink聯合仿真平臺

由于在Ansys中無法對燃氣發生器內部的燃燒過程進行模擬，而Matlab可以采用數值求解的方式，通過建立燃燒與膨脹做功的熱力學模型，對燃燒過程進行計算。利用Ansys中的CoLink模塊建立與Simulink之間的聯系，將多體動力學的求解作為模塊嵌入Simulink中，建立燃燒-展開的聯合仿真模型。

3.1 微分方程組的建立

根據上文所述質量守恒方程、能量守恒方程，得可求解的偏微分方程組：

(13)

式(13)中，燃面Ab、作動筒位移szdt可直接求出；作動筒內溫度Tzdt、作動筒內壓強pzdt須通過方程組解出；如果采用簡化的延伸噴管模型，作動筒速度vzdt、作動筒加速度azdt也可以通過微分方程解出；如果采用聯合仿真，則可以考慮更詳細的動力學模型，通過多體動力學求解器得出vzdt、azdt。

3.2 聯合仿真平臺的建立

將求解偏微分方程組的.m文件寫成函數形式。由于Ansys求解模塊無法嵌入偏微分方程求解中，故采用時間離散的方法逐步求解，框圖如圖6所示。聯合仿真模型如圖7所示。模型中，B為燃燒的肉厚，AAB為實際燃面面積，rm為單藥片計算燃面面積，m為已燃燒藥柱質量，其余變量與前文意義相同；main為微分方程求解主程序，rm/m為藥柱燃燒參數求解程序，Ansys Client Block為多體動力學求解器，Memory1～9為對應參數初始值；由于現有裝藥形式為多片裝藥，Constant2表示裝藥片數；Switch為判斷程序，當裝藥燒完后,燃面面積設置為Constant1；scope為監視器。

圖6 聯合仿真求解框圖Fig.6 Co-simulation flow

圖7 聯合仿真模型Fig.7 Co-simulation model

由于是離散求解，無需考慮連續性問題，且為了避免求解中發生代數環(algebraic loop)錯誤，在每個循環中都加入Memory塊。

4 仿真結果與分析

仿真中微分方程求解為龍格-庫塔法，error tolerance設置為1e-7；整體仿真設置為discrete，步長為0.001，進行聯合求解。

4.1 試驗與仿真結果對比

裝藥片數為7、4時，仿真獲得作動筒內壓力曲線與實測壓力曲線對比如圖8所示。作動筒展開速度曲線與實測速度曲線(高速錄像獲得)對比如圖9所示。

(a)7 charge pellets

(a)7 charge pellets

如表1所示，由實測與仿真結果對比可以看出，不同藥量下作動筒內壓曲線與延伸錐展開曲線吻合程度較好，聯合仿真模型可以較好的模擬出延伸噴管實際的點火-展開過程。

表1 藥片數分別為7和4時仿真值與實測值對比

由數據對比可見，加入已燃燒藥量的判定可以較為準確地預測出作動筒內壓強的峰值、展開到位時壓強值以及燃燒完畢時作動筒殘壓值，最大誤差約10%。

4.2 仿真結果預測與分析

在不同藥量下的延伸錐展開預測中，關鍵點在于作動筒內壓曲線的獲得。不同藥量下作動筒內壓數據如圖10所示(8 pcs為預測數據)?？梢钥闯觯S著藥量的上升，壓強峰值、到位壓強與殘壓均有提升，僅到位壓強可以近似為比例關系；同時，作動筒到位時間、達到最大工作壓強的時間、藥柱燃燒時間均有縮短，但到位后燃氣發生器繼續工作的時間(藥柱燃燒時間-作動筒到位時間)并不會隨之縮短。所以即使在相同的外部邊界條件下，不同藥量的燃氣發生器曲線變化規律也并非簡單的線性關系。

圖10 不同藥量下作動筒內壓Fig.10 Pressure in the actuator under different charge pellets

在保持燃氣發生器藥量的情況下，若外部條件發生變化，如燃氣尾流阻力變化等，延伸錐的展開曲線會受到較大影響，展開的最大速度也會發生變化；同時，作動筒內壓曲線也會相應變化。由仿真分析可以得出，不同展開阻力下延伸噴管的展開情況如圖11、圖12所示。

圖11 不同阻力下作動筒內壓Fig.11 Pressure in the actuator under different resistance

圖12 不同阻力下延伸錐展開速度Fig.12 Extending velocity of the exit cone under different resistance

可以看出，在燃氣尾流的影響下，展開速度明顯下降；同時作動筒內壓升高，到位時間增加，展開到位時作動筒壓強也有所升高。因此，在外部邊界條件改變的時候，作動筒內壓曲線也會相應產生變化。又由于燃氣發生器內壓與作動筒內壓相關，故燃氣發生器內的燃燒情況也會隨外部邊界條件的變化而改變。

5 問題與討論

模型中較為模糊的點在于燃燒室壓強pc與作動筒壓強pzdt的關系。如前文所述，作動筒內壓pzdt(因變量)應為燃氣發生器內壓pc、時間t(自變量)的單值函數，即：

pzdt=f(pc,t)

(14)

但由于試驗條件原因，無法準確測得燃氣發生器內壓pc的值，只能獲得作動筒內壓pzdt。為反應出pzdt對pc的影響，在仿真中采用近似的關系式：

pc=k1pzdte-k2pzdt

(15)

二者隨時間變化的壓強曲線對比如圖13所示。

圖13 燃氣發生器內外壓強對比Fig.13 Comparision between the pressure in the gas generator (inner pressure) and the actuator (outer pressure)

由于最終的仿真曲線重合度較高，可以近似認為燃氣發生器內壓曲線與實際狀況一致。可以看出，在燃氣發生器點燃的初期，作動筒內壓較小，燃氣發生器內壓快速上升；隨著作動筒內壓逐漸上升，作動筒展開加速度變快，同時燃氣發生器內壓下降；在展開即將到位時，作動筒內壓下降，燃氣發生器內壓再次上升，在藥柱燃燒完畢后停止。

最終單位時間藥柱燃燒質量曲線如圖14所示?？梢钥闯?，不同藥量下質量流率變化趨勢基本一致。在展開初期，質量流率快速上升，在作動筒展開時先降后升，在展開到位后緩慢下降。另外，隨著藥量的下降，燃面面積與燃速均有下降，質量流率降低，燃燒時間上升。

圖14 不同藥量下燃氣發生器質量流率曲線Fig.14 Mass flow rate under differentcharge pellets

為進一步完善仿真模型，后續需要對燃氣發生器內壓曲線進行準確的測量，以期獲得燃氣發生器內壓、作動筒內壓(燃氣發生器外壓)與時間三者的函數關系，或者通過流場計算得出 ，才能對延伸噴管點火-展開全過程進行更加準確的模擬。

6 結論

本文通過分析延伸噴管燃燒-傳遞-做功-展開的全過程，建立了燃氣發生器裝藥點火燃燒與燃氣做功的數學模型，以及延伸噴管展開的動力學模型，并采用聯合仿真平臺將進行耦合仿真。結論如下：

(1)聯合仿真計算結果與試驗數據一致性較好，表明聯合仿真模型合理，仿真計算方法正確，仿真計算結果可靠；

(2)隨著燃氣發生器藥量的提高，作動筒內壓峰值提高，到位時間提前且殘壓上升；

(3)若延伸噴管外部阻力上升，作動筒內壓會相應提高，但展開最大速度下降，到位時間增加；

(4)仿真模型可以較好地預估出在不同燃氣發生器藥量、不同發動機工作條件下的延伸噴管展開位移曲線，為在不同工況下的藥量調節提供的可靠的參考，具有較強的工程應用價值；

(5)仿真計算中發現，燃氣發生器內壓、作動筒內壓、燃氣發生器工作時間三者間存在較強的關系，研究其變化規律有助于進一步完善仿真模型，以獲得更準確，適用范圍更廣的結果。