吝 琳，尚永騰，李沖沖，白 濤，郭夢飛

(1.航空工業(yè)第一飛機設(shè)計研究院，西安 710089；2.中國空空導(dǎo)彈研究院，河南 洛陽 471009；3.西安航空學(xué)院 飛行器學(xué)院，西安 710077；4.西安理工大學(xué) 機械與精密儀器工程學(xué)院，西安 710048)

0 引 言

固體發(fā)動機在深空探測、重載運輸和導(dǎo)彈武器等方面得到了廣泛的應(yīng)用[1-4]，其內(nèi)彈道性能參數(shù)直接影響發(fā)動機以及飛行器的性能，且裝藥推移過程中燃面的精確計算是固體發(fā)動機燃燒室內(nèi)彈道性能計算的基礎(chǔ)，是發(fā)動機設(shè)計中的關(guān)鍵參數(shù)。燃面計算的目的是獲得發(fā)動機燃燒室內(nèi)固體推進劑燃燒的面積。發(fā)動機工作過程中，推進劑會發(fā)生燃燒推移，導(dǎo)致燃燒面積變化，需要精確計算出發(fā)動機的燃面面積隨工作時間的變化關(guān)系。

隨著復(fù)雜藥型以及高填充裝藥發(fā)動機的使用[5-6]，燃面面積、燃面推移經(jīng)歷了二維至三維的變化，計算方法也由之前的“解析法” “作圖法”等方法發(fā)展到目前的“Level-Set” “最小距離函數(shù)法(MDF)”。解析法是指使用數(shù)學(xué)模型對藥柱進行建模，具有高的精度；作圖法即依據(jù)藥柱的幾何形狀進行裝藥推移模擬。解析法和作圖法不適用于三維復(fù)雜的藥柱。為了解決三維復(fù)雜藥柱的燃面計算問題，Osher等[7]提出使用Level-Set界面追蹤方法解決活動邊界問題，具有較好的精度。秦飛等[8]使用Level-Set方法計算了發(fā)動機燃燒室的燃面面積。然而，Level-Set方法在遇到曲率過小或過大時求解存在問題[9]。Willcox等[10-11]提出最小距離函數(shù)(MDF)燃面計算方法，并使用該方法開展了發(fā)動機裝藥燃面推移和燃面面積計算。郭夢飛等[12]使用MDF方法計算星孔裝藥的燃面面積，發(fā)現(xiàn)MDF方法穩(wěn)定性好，精度高。馬長禮等[13]優(yōu)化了MDF方法的計算過程，提高了計算效率，并使用不同藥型驗證了MDF方法的通用性。彭博等[14]使用復(fù)雜藥柱，驗證了MDF方法的通用性和高精度，認為其適用于復(fù)雜裝藥的燃面計算。王革等[9]將MDF方法的計算思路引入Level-Set方法中，提高了Level-Set方法的精度。Hwang等[15]對比了Level-Set方法和MDF方法的區(qū)別，結(jié)果表明，MDF方法的精度高于Level-Set方法。

固體發(fā)動機面臨多種工作環(huán)境[16-17]，包括高的橫向過載。例如在第四代空空導(dǎo)彈中，其主動段的橫向過載達到了60～100。高橫向過載下發(fā)動機燃燒室內(nèi)推進劑出現(xiàn)不等速燃燒現(xiàn)象，燃面相應(yīng)地發(fā)生變化，影響發(fā)動機的壓強和推力；同時，推進劑的不等速燃燒可能導(dǎo)致部分裝藥提前燃燒完畢，暴露出發(fā)動機的熱防護層[18-19]，過載導(dǎo)致的粒子射流會惡化熱防護層的工作環(huán)境，嚴重時甚至引起發(fā)動機燒穿爆炸[20-22]，目前關(guān)于過載情況下發(fā)動機內(nèi)燃面面積計算的研究較少，限制了過載工況下發(fā)動機內(nèi)彈道性能的精確計算和發(fā)動機的優(yōu)化設(shè)計。

為了解決高橫向過載時發(fā)動機燃面精確計算的難題，本文使用MDF方法，針對高橫向過載的發(fā)動機工況，計算獲得推進劑裝藥燃燒時的燃面面積、燃面推移以及發(fā)動機內(nèi)彈道性能，為高過載發(fā)動機的設(shè)計優(yōu)化提供理論基礎(chǔ)。

1 計算模型

1.1 MDF方法

MDF是一種可以快速計算固體發(fā)動機復(fù)雜三維裝藥燃面推移的方法。該方法需要一個包含藥柱的三維最小距離網(wǎng)格和初始燃燒面，以此計算出每個網(wǎng)格點至初始面的距離，同時判斷出該點的距離符號，按照正號代表在藥柱上、負號代表在燃燒室空腔的原則，燃面自然就處于正負最小距離的網(wǎng)格點之間，即最小距離為零的點。依據(jù)裝藥的燃燒速度進行推移，得到不同時刻每個網(wǎng)格點的推移量，依據(jù)推移量重新計算最小距離函數(shù)，并抽取出發(fā)動機下一個時刻的燃面點，將燃面點組成燃面并進行面積計算，以此獲得實時的燃面面積。

MDF方法計算步驟分為三步：(1)裝藥初始化。將發(fā)動機裝藥使用CAD軟件進行STL離散，使用C++程序讀取發(fā)動機的網(wǎng)格數(shù)據(jù)和初始燃燒曲面，并進行最小距離計算。(2)燃面推移。依據(jù)每個點的燃燒速率進行發(fā)動機裝藥的推移模擬，得到不同時刻的燃面外形。(3)燃面面積計算。根據(jù)每個時刻得出的燃面推移結(jié)果，得出每個時刻的燃面網(wǎng)格點，進而求出燃面面積。國內(nèi)外學(xué)者的研究均證實，MDF方法在處理燃面推移和一維內(nèi)彈道性能時具有精度高、通用性強和穩(wěn)定性好的特點[12-17]。

1.2 燃燒速率計算模型

燃燒速率的計算是推進劑推移過程的關(guān)鍵參數(shù)。首先，計算出推進劑的初始燃燒速率r0:

r0=apn

(1)

其次，引入Lenoir-Robillard模型[23-24]，計算侵蝕對燃燒速率的影響:

(2)

最后，建立過載工況下推進劑穩(wěn)態(tài)燃速增大模型[25-26]，獲得過載工況下推進劑的燃燒速率，采用steffensen方法進行迭代求解:

(3)

式中：r0為基本燃速；a為燃速系數(shù)；p為燃燒室壓強；n為壓強指數(shù)；ρ為燃氣密度；u為燃氣速度；L為藥柱軸向距離；α主要取決于氣流中心溫度和燃燒產(chǎn)物特性；rg為過載下推進劑燃速；β為熱流系數(shù)；Ga為加速度質(zhì)量通量;ρp為推進劑密度；cp為燃氣定壓比熱容；δ0為燃燒區(qū)厚度；λ為氣相熱導(dǎo)率。

研究表明，Lenoir-Robillard燃速模型和過載燃速增大模型均具有較高的精度[23-26]。

過載方位角φ的計算如下：

(1)由燃面三角形網(wǎng)格的三個頂點確定出三角形的平面方程ax+by+cz+d=0，其法向矢量為(a，b，c)；

(2)通過已知的過載條件, 得到過載的方向矢量(l，m，n)；

(3)計算出過載方位角φ：

(4)

1.3 內(nèi)彈道性能計算模型

在計算過程中，如果滿足燃氣流速遠小于當?shù)芈曀?，燃氣密度遠小于推進劑密度，燃氣通道面積變化遠小于燃氣通道面積，則燃氣流動參數(shù)隨時間的變化很小，這時認為燃氣在裝藥通道中的運動是準定常的流動，控制方程中對時間的偏導(dǎo)數(shù)各項可以忽略。此外，忽略發(fā)動機工作期間裝藥通道內(nèi)燃氣向燃燒室壁的散熱損失，通道內(nèi)燃氣的一維準定常流動可以認為是絕能流動，此時得到絕能條件下的一維準定常控制方程組為[27-28]

(5)

由式(1)～(5)獲得一維內(nèi)彈道準定常絕能流動條件下的控制方程組，具體的推導(dǎo)過程見文獻[29]：

(6)

式中：ρ為燃氣密度；Ap為通氣截面積；Ab為燃燒表面積；Π為藥柱燃燒周長；r為推進劑燃燒速率；p為靜壓；T為靜溫；R為氣體常數(shù)；Tf為總溫。

計算時，將導(dǎo)彈在彈身長度方向平均分成11個截面，令第一個截面的流速為0，使用平衡壓強作為壓強迭代的初始值p1，迭代計算每個截面的燃燒面積和通道面積，通過截面裝藥基礎(chǔ)燃速的線性插值，得到各個點的基礎(chǔ)燃燒速度。

2 結(jié)果與分析

2.1 燃燒速率

將推進劑在無過載時(即過載值=0)的燃燒速率定義為1，使用式(3)計算獲得不同過載值時推進劑燃燒速率的比值。過載方位角和過載值對發(fā)動機裝藥燃燒速率的影響如圖1所示。可以看出，當過載方位角較小時，過載值對燒蝕速率影響較大。當方位角為0°，過載值為100g時，發(fā)動機裝藥的燃速增大系數(shù)為1.34，燃燒速率顯著增大，影響發(fā)動機燃燒室內(nèi)的燃面面積，導(dǎo)致發(fā)動機的內(nèi)彈道性能發(fā)生較大變化。當過載方位角較大時，過載值對燒蝕速率影響較小，過載繼續(xù)增大，燃速的增加不顯著，此時過載對燃速的影響可以忽略不計。

圖1 過載方位角和過載值對發(fā)動機裝藥燃燒速率的影響

2.2 過載工況下燃面面積計算

采用工程中常用的圓管+星孔裝藥，前段裝藥為圓柱狀，后段裝藥為星孔狀，見圖2。推進劑裝藥參數(shù)：推進劑密度1 730 kg/m3，推進劑比熱1 553.4 J/(kg·K)，燃面絕熱燃燒溫度3 400 K，燃速壓強指數(shù)0.35，燃速參考壓強6 MPa，推進劑參考燃速9 mm/s，燃氣平均氣體常數(shù)290.5 J/(kg·K)。

圖2 裝藥結(jié)構(gòu)示意圖

根據(jù)藥柱尺寸，將其劃分為50×50×201個基礎(chǔ)網(wǎng)格點，同時在圓管與星孔的結(jié)合部進行局部加密，選取藥柱推移步長為2 mm，推移步數(shù)為n。在燃面推移計算時，考慮侵蝕燃燒和橫向過載(80g)對發(fā)動機燃燒速率的影響，使用MDF方法計算裝藥的推移過程，進而獲得不同時刻的燃面面積，圖3～6給出了星孔藥柱的燃面推移情況。

圖3 橫向過載下的裝藥燃面推移(推移步數(shù)1～12)

圖3和圖4中推移步數(shù)0～12時使用藥柱末端為主視角，隨著步數(shù)的推移，藥柱內(nèi)徑逐漸增大，空腔也隨之增大。在推移步數(shù)為12時，橫向過載下裝藥已經(jīng)燃燒至星孔裝藥的邊界，裝藥末端的一部分已經(jīng)燃盡(見圖3)，而在無橫向過載時裝藥并沒有燃盡(見圖4)。由于過載的存在，使得藥柱與過載方向一致的部分燃燒速度增大，藥柱較早地出現(xiàn)了偏燒現(xiàn)象。

圖4 無橫向過載時的裝藥燃面推移(推移步數(shù)1～12)

接下來將視角切為藥柱整體。無橫向過載時，隨著燃燒繼續(xù)，裝藥也燃燒至星孔裝藥的邊界，星孔段大面積減少，t=19 s時，星孔段燃盡，一級發(fā)動機工作完成，如圖5所示。再次將視角切回藥柱尾部，此時圓柱段內(nèi)腔持續(xù)增大，燃燒面積也增大，最后剩下薄壁圓筒迅速燃燒完畢，發(fā)動機進入二級低壓強、低推力的巡航階段。

圖5 無橫向過載時的裝藥燃面推移(推移步數(shù)13～24)

在橫向過載下，發(fā)動機出現(xiàn)了明顯的偏燒，裝藥下部快速地燃燒完畢，而裝藥的上部燃燒速率相對較慢，如圖6所示。由于橫向過載的影響，使得藥柱上部的燃燒時間明顯長于藥柱下部，導(dǎo)致星孔段藥柱和圓柱段藥柱燃盡的時間變長，且燃燒后期出現(xiàn)了大量余藥。此時，余藥燃面面積迅速下降，發(fā)動機工作壓強降低，藥柱的燃燒速度低，發(fā)動機的拖尾段時間增長。

圖6 橫向過載下的裝藥燃面推移(推移步數(shù)12～24)

依據(jù)裝藥燃面的推移過程可以獲得裝藥燃面隨時間的變化關(guān)系，如圖7所示。與無橫向過載相比，有橫向過載時藥柱的燃面面積出現(xiàn)了明顯變化?？梢钥闯觯?/p>

圖7 發(fā)動機燃面面積及壓強隨時間的變化關(guān)系

(1)一維內(nèi)彈道壓強的總體變化趨勢與燃面的變化趨勢基本一致；

(2)藥柱在過載一側(cè)的燃速增加較為明顯，但并沒有引起燃面面積的大幅變化以及最大壓強的過高增加；

(3)有過載時，藥柱出現(xiàn)了偏燒現(xiàn)象，發(fā)動機具有明顯的拖尾段，導(dǎo)致發(fā)動機的總工作時間變長，此時，藥柱的燃燒面積小，壓強低。

由MDF計算結(jié)果可以看出，橫向過載不會引起較為明顯的燃燒面積增加和壓強的增大。其主要的影響是改變?nèi)紵覂?nèi)藥柱的燃燒速率，造成藥柱的偏燒，從而導(dǎo)致部分發(fā)動機殼體過早地暴露在高溫高壓的兩相燃氣中，對發(fā)動機的熱防護系統(tǒng)設(shè)計提出了很高的要求。

3 結(jié) 論

針對固體火箭發(fā)動機過載這一特殊工況，首次使用MDF方法計算了發(fā)動機裝藥的燃面面積以及發(fā)動機的藥柱推移過程，并將發(fā)動機一維內(nèi)彈道與燃面計算耦合，獲得了過載工況下發(fā)動機裝藥燃面面積、燃燒室內(nèi)壓強等參數(shù)隨時間的變化曲線。研究結(jié)果表明，橫向過載會引起發(fā)動機藥柱出現(xiàn)偏燒現(xiàn)象，改變了發(fā)動機裝藥的燃面面積和燃燒室內(nèi)的壓強，并延長了發(fā)動機的工作時間。同時，偏燒會導(dǎo)致部分藥柱提前燃燒完畢，需要對發(fā)動機熱防護系統(tǒng)進行針對性的加強。