基于Simulink的延伸噴管燃?xì)庹归_過程聯(lián)合仿真*

陳 鵬，宋學(xué)宇，曹濤鋒，尤軍峰

(1.中國航天科技集團(tuán)公司四院四十一所，西安 710025；2.固體火箭發(fā)動機燃燒、熱結(jié)構(gòu)與內(nèi)流場國防科技重點實驗室，西安 710025；3.西北工業(yè)大學(xué) 航天學(xué)院，西安 710072)

0 引言

延伸噴管技術(shù)是提高固體火箭發(fā)動機性能的關(guān)鍵技術(shù)之一。對于采用作動筒推動方式展開的延伸噴管，作動筒展開力是決定延伸噴管動力學(xué)特性的關(guān)鍵因素。作動筒的展開方式包括燃?xì)馐?、氣瓶式、產(chǎn)氣式等。

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

現(xiàn)有的仿真計算局限于單系統(tǒng)的建模仿真，如燃?xì)獍l(fā)生器的內(nèi)彈道仿真計算，氣體充壓過程仿真計算[5]，延伸噴管展開動力學(xué)、展開碰撞、燃?xì)馕擦鞣抡嬗嬎鉡6-11]等。但由于燃?xì)獍l(fā)生器內(nèi)彈道參數(shù)與下游過程相關(guān)，且受延伸噴管展開系統(tǒng)機構(gòu)的影響，下游氣體做功過程無法通過簡單試驗?zāi)M，必須聯(lián)合整個點火展開系統(tǒng)進(jìn)行試驗，周期較長且成本較高。在噴管外載荷不確定，且展開機構(gòu)受加工工藝影響較大的條件下，僅通過單個系統(tǒng)的仿真計算，無法給出藥柱與展開動力學(xué)參數(shù)的關(guān)系。而如何確定燃?xì)獍l(fā)生器藥形與藥量，是在工程實際中必須予以解決的問題。

針對作動筒式雙級延伸噴管采用燃?xì)獍l(fā)生器展開的工況，本文通過在Matlab/Simulink中構(gòu)建燃?xì)獍l(fā)生器藥柱燃燒與作動筒內(nèi)燃?xì)獾呐蛎涀龉δＰ停Y(jié)合Ansys的動力學(xué)仿真模塊，對從燃?xì)獍l(fā)生器點火開始的延伸噴管展開全過程進(jìn)行聯(lián)合仿真，以期在燃?xì)獍l(fā)生器藥型與藥量發(fā)生變化時，可以對延伸噴管展開過程進(jìn)行預(yù)測。通過與已有試驗測試數(shù)據(jù)的對比，驗證仿真結(jié)果的準(zhǔn)確性，建立完整可靠的延伸噴管點火展開模型。

1 燃?xì)庹归_模型

在延伸噴管的展開過程中，作動筒的充氣過程從燃?xì)獍l(fā)生器點火開始。不計點火時間，主裝藥開始燃燒并產(chǎn)生高焓值的氣體，經(jīng)過濾網(wǎng)后進(jìn)入管路，最終進(jìn)入作動筒并膨脹做功，推動作動筒展開，帶動延伸噴管展開到位。通過建立藥柱燃燒及燃?xì)馀蛎涀龉Φ哪Ｐ?，可以為延伸噴管動力學(xué)模型提供輸入條件，再經(jīng)動力學(xué)仿真，最終得到延伸噴管展開的動力學(xué)參數(shù)。

1.1 燃燒過程

由維耶里模型可知，燃速rc與燃燒室內(nèi)壓pc的關(guān)系為

rc=a0pcn

(1)

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

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

故燃面Ab為

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

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

(2)

(3)

由質(zhì)量守恒可知，藥柱燃燒減少的質(zhì)量，一部分成為殘渣被濾網(wǎng)阻攔，一部分作為燃?xì)膺M(jìn)入作動筒內(nèi)膨脹做功，則有

(4)

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

pc=k1pzdtek2pzdt

(5)

注意到此關(guān)系內(nèi)不含時間項t，且對于同一個燃燒室內(nèi)壓pc可能對應(yīng)不同的作動筒內(nèi)壓pzdt。

1.2 膨脹做功過程

理想氣體方程:

(6)

結(jié)合質(zhì)量守恒方程，有

(7)

開口系統(tǒng)能量守恒方程:

(8)

由于每個燃?xì)獍l(fā)生器為2個作動筒供氣，有

(9)

1.3 動力學(xué)展開過程

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

展開阻力來源有三部分：(1)作動筒展開內(nèi)部阻力，由預(yù)緊力產(chǎn)生的摩擦力與展開過程中筒間正壓力產(chǎn)生的摩擦力兩部分組成；(2)展開即將到位時產(chǎn)生的阻力，由延伸錐側(cè)向密封圈摩擦力、延伸錐端面密封圈阻力、折轉(zhuǎn)片到位后的拉力三部分組成；(3)熱試中尾部燃?xì)鈱ρ由戾F產(chǎn)生的作用力，根據(jù)燃?xì)獠煌呐蛎洜顟B(tài)，可能產(chǎn)生引射力或者阻力，其值可由流場計算得到。

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

Vzdt=V0+Azdtszdt

(10)

由運動關(guān)系可知：

(11)

在噴管展開方向(軸向)上，作動筒內(nèi)壓(考慮作動筒內(nèi)摩擦fzdt)的分量克服軸向摩擦力ff、尾部燃?xì)庾枇e使噴管展開。由于展開機構(gòu)限制，延伸噴管I(下標(biāo)1)、II(下標(biāo)2)級之間的位移關(guān)系固定為1∶2，故在軸向有

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

(12)

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

不考慮碰撞損失，延伸錐展開到位后，作動筒內(nèi)燃?xì)庾龉ν耆脕砜朔Σ亮拔膊咳細(xì)庾枇Α?/p>

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

2 Ansys多體動力學(xué)模型

在Ansys中對1/4延伸噴管進(jìn)行建模，給定作動筒上的摩擦阻力、到位等效阻力以及燃?xì)鈱ρ由戾F的作用力，以作動筒內(nèi)壓作為輸入量，輸出作動筒上端在作動筒展開方向上的的速度和加速度。三維模型如圖2所示。

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

在模型中，作動筒阻力通過設(shè)置作動筒上的預(yù)緊力和摩擦系數(shù)給出；延伸錐展開到位阻力的產(chǎn)生包含了與密封圈摩擦和碰撞過程，模型中采用近似曲線進(jìn)行模擬，圖3中縱軸為延伸錐展開到位阻力(等效在II級延伸錐上)，橫軸為II級延伸錐質(zhì)心軸向坐標(biāo)。亦可采用顯式動力學(xué)計算出阻力隨展開位移變化曲線，作為參數(shù)輸入模型中；尾部燃?xì)猱a(chǎn)生的阻力與發(fā)動機工作狀態(tài)以及延伸噴管展開位置有關(guān)，在發(fā)動機工作狀態(tài)穩(wěn)定的情況下，采用流場計算得到尾部燃?xì)庾枇﹄S延伸錐展開位移的變化曲線，并作為參數(shù)輸入模型中。圖4中縱軸為兩級延伸錐所受的燃?xì)庾枇?，橫軸為II級延伸錐質(zhì)心軸向坐標(biāo)。

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

圖4 噴管燃?xì)庾枇η€Fig.4 Gas resistance of the extendible nozzle

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

在動力學(xué)模型建立完成后，通過輸入實測的作動筒內(nèi)壓曲線，可較為準(zhǔn)確地仿真出延伸錐展開曲線，如圖5所示；若給出期望的展開位移曲線，也可計算出需求的作動筒內(nèi)壓曲線。即作動筒內(nèi)壓曲線與延伸錐展開曲線之間的動力學(xué)模型是完整可靠的。但如果燃?xì)獍l(fā)生器藥量發(fā)生改變，很難通過對內(nèi)壓曲線的線性變化來準(zhǔn)確預(yù)示新的展開曲線，必須結(jié)合燃?xì)獍l(fā)生器內(nèi)部的燃燒與燃?xì)馀蛎涀龉Φ倪^程進(jìn)行求解，進(jìn)行多系統(tǒng)的聯(lián)合仿真。

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

3 Simulink聯(lián)合仿真平臺

由于在Ansys中無法對燃?xì)獍l(fā)生器內(nèi)部的燃燒過程進(jìn)行模擬，而Matlab可以采用數(shù)值求解的方式，通過建立燃燒與膨脹做功的熱力學(xué)模型，對燃燒過程進(jìn)行計算。利用Ansys中的CoLink模塊建立與Simulink之間的聯(lián)系，將多體動力學(xué)的求解作為模塊嵌入Simulink中，建立燃燒-展開的聯(lián)合仿真模型。

3.1 微分方程組的建立

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

(13)

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

3.2 聯(lián)合仿真平臺的建立

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

圖6 聯(lián)合仿真求解框圖Fig.6 Co-simulation flow

圖7 聯(lián)合仿真模型Fig.7 Co-simulation model

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

4 仿真結(jié)果與分析

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

4.1 試驗與仿真結(jié)果對比

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

(a)7 charge pellets

(a)7 charge pellets

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

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

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

4.2 仿真結(jié)果預(yù)測與分析

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

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

在保持燃?xì)獍l(fā)生器藥量的情況下，若外部條件發(fā)生變化，如燃?xì)馕擦髯枇ψ兓?，延伸錐的展開曲線會受到較大影響，展開的最大速度也會發(fā)生變化；同時，作動筒內(nèi)壓曲線也會相應(yīng)變化。由仿真分析可以得出，不同展開阻力下延伸噴管的展開情況如圖11、圖12所示。

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

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

可以看出，在燃?xì)馕擦鞯挠绊懴?，展開速度明顯下降；同時作動筒內(nèi)壓升高，到位時間增加，展開到位時作動筒壓強也有所升高。因此，在外部邊界條件改變的時候，作動筒內(nèi)壓曲線也會相應(yīng)產(chǎn)生變化。又由于燃?xì)獍l(fā)生器內(nèi)壓與作動筒內(nèi)壓相關(guān)，故燃?xì)獍l(fā)生器內(nèi)的燃燒情況也會隨外部邊界條件的變化而改變。

5 問題與討論

模型中較為模糊的點在于燃燒室壓強pc與作動筒壓強pzdt的關(guān)系。如前文所述，作動筒內(nèi)壓pzdt(因變量)應(yīng)為燃?xì)獍l(fā)生器內(nèi)壓pc、時間t(自變量)的單值函數(shù)，即：

pzdt=f(pc,t)

(14)

但由于試驗條件原因，無法準(zhǔn)確測得燃?xì)獍l(fā)生器內(nèi)壓pc的值，只能獲得作動筒內(nèi)壓pzdt。為反應(yīng)出pzdt對pc的影響，在仿真中采用近似的關(guān)系式：

pc=k1pzdte-k2pzdt

(15)

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

圖13 燃?xì)獍l(fā)生器內(nèi)外壓強對比Fig.13 Comparision between the pressure in the gas generator (inner pressure) and the actuator (outer pressure)

由于最終的仿真曲線重合度較高，可以近似認(rèn)為燃?xì)獍l(fā)生器內(nèi)壓曲線與實際狀況一致?？梢钥闯?，在燃?xì)獍l(fā)生器點燃的初期，作動筒內(nèi)壓較小，燃?xì)獍l(fā)生器內(nèi)壓快速上升；隨著作動筒內(nèi)壓逐漸上升，作動筒展開加速度變快，同時燃?xì)獍l(fā)生器內(nèi)壓下降；在展開即將到位時，作動筒內(nèi)壓下降，燃?xì)獍l(fā)生器內(nèi)壓再次上升，在藥柱燃燒完畢后停止。

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

圖14 不同藥量下燃?xì)獍l(fā)生器質(zhì)量流率曲線Fig.14 Mass flow rate under differentcharge pellets

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

6 結(jié)論

本文通過分析延伸噴管燃燒-傳遞-做功-展開的全過程，建立了燃?xì)獍l(fā)生器裝藥點火燃燒與燃?xì)庾龉Φ臄?shù)學(xué)模型，以及延伸噴管展開的動力學(xué)模型，并采用聯(lián)合仿真平臺將進(jìn)行耦合仿真。結(jié)論如下：

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

(2)隨著燃?xì)獍l(fā)生器藥量的提高，作動筒內(nèi)壓峰值提高，到位時間提前且殘壓上升；

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

(4)仿真模型可以較好地預(yù)估出在不同燃?xì)獍l(fā)生器藥量、不同發(fā)動機工作條件下的延伸噴管展開位移曲線，為在不同工況下的藥量調(diào)節(jié)提供的可靠的參考，具有較強的工程應(yīng)用價值；

(5)仿真計算中發(fā)現(xiàn)，燃?xì)獍l(fā)生器內(nèi)壓、作動筒內(nèi)壓、燃?xì)獍l(fā)生器工作時間三者間存在較強的關(guān)系，研究其變化規(guī)律有助于進(jìn)一步完善仿真模型，以獲得更準(zhǔn)確，適用范圍更廣的結(jié)果。