基于殘值函數(shù)法的雙燃速藥柱內(nèi)彈道分析*

褚佑彪，趙天泉，董新剛，任 萍

(中國航天科技集團有限公司四院四十一所，固體火箭發(fā)動機燃燒、熱結(jié)構(gòu)與內(nèi)流場國防科技重點實驗室,西安 710025)

0 引言

燃面計算用于確定裝藥在燃燒過程中燃燒表面積隨燃燒時間的變化規(guī)律，直接影響發(fā)動機內(nèi)彈道性能預(yù)示精度，是發(fā)動機內(nèi)彈道設(shè)計的基礎(chǔ)，在固體火箭發(fā)動機的設(shè)計中一直占有重要地位[1]。為滿足先進固體動力技術(shù)的發(fā)展要求，燃燒室的藥型設(shè)計越來越復(fù)雜。在單室雙推發(fā)動機的設(shè)計過程中，為保證助推段與續(xù)航段的推力比，常采用高低燃速搭配的思路進行裝藥設(shè)計。當(dāng)高低燃速推進劑在交界面處燃燒時，燃面會出現(xiàn)分離、交匯等復(fù)雜的拓?fù)浣Y(jié)構(gòu)變化，是燃面推移計算的一個難點。

目前，可實現(xiàn)燃面不等速推移的有實體造型法[2]、網(wǎng)格推移法[3-4]、最小距離函數(shù)法[5-6]、Level Set方法[7-8]和殘值函數(shù)法[9]等。其中，實體造型法通過特征造型或驅(qū)動尺寸實現(xiàn)燃面推移過程的模擬。但是，對于結(jié)構(gòu)復(fù)雜的藥型，實體造型法的推移過程十分繁瑣，并可能出現(xiàn)奇點，導(dǎo)致計算無法繼續(xù)；對于多燃速藥柱，燃面的分離需要在推移過程中構(gòu)建新的幾何特征，較難實現(xiàn)。網(wǎng)格推移法在初始藥型上生成非結(jié)構(gòu)網(wǎng)格，利用當(dāng)?shù)厝妓龠M行推移，再通過網(wǎng)格重構(gòu)獲得新的燃面，該方法通用性較好，但其處理燃面交匯、分離等拓?fù)浣Y(jié)構(gòu)變化的穩(wěn)定性較差。最小距離函數(shù)法通過計算藥柱內(nèi)部各點到初始燃面的距離，即最小距離函數(shù)(MDF)，進行燃面推移分析，可以實現(xiàn)不等速推移，但是其無法自然捕捉到由于燃速間斷的存在而導(dǎo)致的燃面拓?fù)浣Y(jié)構(gòu)的變化。Level Set方法采用初值形式的偏微分方程將一個純幾何問題轉(zhuǎn)變?yōu)橛闷⒎址匠堂枋龅臄?shù)學(xué)問題，對于推進劑內(nèi)部存在燃速間斷的工況，由于燃?xì)馀c推進劑之間也存在間斷，在數(shù)值求解微分方程組時，會遇到間斷相交的情況，進而導(dǎo)致燃面計算精度及穩(wěn)定性較差。殘值函數(shù)法[9]在笛卡爾網(wǎng)格上，采用殘值函數(shù)記錄燃面位置，通過惠更斯原理模擬燃面的推移過程。此算法可以準(zhǔn)確捕捉由于燃速間斷的存在導(dǎo)致燃面的交匯、分離、消失等復(fù)雜拓?fù)浣Y(jié)構(gòu)變化，且不需要進行網(wǎng)格重構(gòu)，穩(wěn)定性好。

本文應(yīng)用殘值函數(shù)法模擬了多燃速藥柱的燃面推移過程，并與試驗結(jié)果進行對比分析。在此基礎(chǔ)上，開展參數(shù)研究，分析了高低燃速推擠劑交界面的位置對藥柱燃面演化過程和內(nèi)彈道的影響。

1 內(nèi)彈道計算方法

為準(zhǔn)確分析雙燃速藥柱在界面處的耦合燃速過程，本文采用殘值函數(shù)法[9]進行燃面計算分析。該方法基于燃面推移的一般性原理——惠更斯球面波傳播原理，不僅可以處理燃面的等速推移，還可以處理燃面的不等速推移，包括燃面上存在燃速間斷的工況。由于殘值函數(shù)法采用笛卡爾網(wǎng)格離散藥柱，對燃面的推移區(qū)域可以有效地提前預(yù)估，因此可以將三維體循環(huán)縮減為具有一定厚度的曲面循環(huán)，大大降低了計算量。此外，在燃面計算過程中引入殘值函數(shù)，不僅有效地控制燃面推移的累計誤差，還可準(zhǔn)確捕捉燃面推移過程中形狀變化，自然處理界面拓?fù)浣Y(jié)構(gòu)變化，避免燃面的網(wǎng)格重構(gòu)，提高算法的穩(wěn)定性。

殘值函數(shù)法中燃面的推移流程可簡要描述為：

第一步：如圖1所示，Pi為當(dāng)前燃面上的網(wǎng)格點。預(yù)估燃面上的離散點Pi的影響區(qū)域，并向四周推移，在推移的區(qū)域內(nèi)，推進劑轉(zhuǎn)化為燃?xì)猓埔凭嚯xΔi由式(1)給出。

Δi=dt×ri+δi

(1)

式中 dt為時間步長；ri為網(wǎng)格點Pi處推進劑的燃速；δi為Pi處的殘值函數(shù)。

第二步：計算并更新殘值函數(shù)。Pj為藥柱內(nèi)的網(wǎng)格點且處于Pi的影響范圍內(nèi)。Pj處的殘值函數(shù)δj由式(2)獲得。

(2)

其中，m為當(dāng)前燃面上能夠影響到Pj的網(wǎng)格點數(shù)目，εij的值根據(jù)Pi與Pj之間是否存在燃速間斷進行分類計算：

(1)如圖1(a)所示，若Pi與Pj之間沒有燃速間斷，當(dāng)Pi與Pj之間的距離dij趨于0時，相應(yīng)的燃速ri和rj趨于相等，則推進劑沿直線燃燒，所以Pj處的εij可以由

εij=Δi-dij

(3)

獲得。

(2)如圖1(b)所示，若Pi與Pj之間存在燃速間斷，類比折射定律可知，推進劑燃燒路徑滿足：

(a) No discontinuity in burning rate (b) Discontinuity in burning rate

risinθj=rjsinθi

(4)

式中θi和θj分別為兩種推進劑燃燒路徑與推進劑交界面法向量的夾角。

O點為推進劑燃燒路徑與交界面的交點。因此，Pj處的εij可由式(5)獲得：

εij=(Δi-diO)rj/ri-dOj

(5)

式中diO為Pi到O點的直線距離；dOj為O點到Pj點的直線距離。

第三步：確定新的燃面，即更新后的燃?xì)鈪^(qū)域與藥柱區(qū)域的交界面。

結(jié)合內(nèi)彈道基本方程，可以計算燃燒室內(nèi)壓強，獲得發(fā)動機內(nèi)彈道性能。

2 某雙燃速藥柱內(nèi)彈道分析

單室雙推發(fā)動機的設(shè)計過程中，為保證助推級與續(xù)航級的推力比，常采用高低燃速搭配的思路進行助推級和續(xù)航級裝藥設(shè)計。圖2為某發(fā)動機藥柱結(jié)構(gòu)，綠色表示助推級藥柱，藍(lán)色表示續(xù)航級藥柱。助推級采用低燃速大燃面藥型，藥柱結(jié)構(gòu)為盲孔和翼柱相結(jié)合結(jié)構(gòu)，續(xù)航級采用小燃面高燃速端燃藥柱結(jié)構(gòu)。推進劑性能參數(shù)見表1。

表1 推進劑性能參數(shù)

圖2 藥柱結(jié)構(gòu)圖Fig.2 Structure diagram of the grain

圖3給出采用殘值函數(shù)法獲得的發(fā)動機內(nèi)彈道數(shù)據(jù)，與試驗結(jié)果進行對比可知，殘值函數(shù)法準(zhǔn)確地模擬了雙燃速藥柱燃燒過程燃面的變化過程。

圖3 數(shù)值仿真與試驗結(jié)果對比Fig.3 Comparison of results of numerical simulation and experiment

由試驗結(jié)果可知，在助推級與續(xù)航級轉(zhuǎn)級過程中，發(fā)動機壓強由7.8 MPa急劇下降至1.5 MPa，然后逐漸增加至局部極大值6.4 MPa，最終逐漸減少至2.6 MPa，構(gòu)成一個非常明顯的壓強起伏。圖4給出典型時刻藥柱結(jié)構(gòu)圖，圖中給出藥柱過軸線的剖面圖，紅色區(qū)域代表燃?xì)馓畛鋮^(qū)域，綠色代表助推級低燃速推進劑，藍(lán)色代表續(xù)航級高燃速推進劑。如圖4(a)所示，燃面與高低燃速推進劑交界面剛剛接觸，此時，低燃速藥柱側(cè)面尚有3.5 mm厚推進劑。隨著燃面的推移，在t=7.97 s時刻，燃面推移至高燃速藥柱內(nèi)，燃面迅速“膨大”，但是由于低燃速藥柱燃面減少量顯著，導(dǎo)致燃燒室壓強迅速降低。在t=7.97～10.28 s范圍內(nèi)，燃面在高燃速藥柱中不斷“膨大”，且低燃速藥柱的兩側(cè)同時燃燒，燃燒室壓強逐漸回升。隨著燃面的不斷推移，藥柱燃面逐漸趨于平面，面積減少，低燃速藥柱燃盡，導(dǎo)致燃燒室壓強趨于平穩(wěn)。由以上分析可知，此雙燃速藥柱發(fā)動機的內(nèi)彈道呈現(xiàn)顯著起伏現(xiàn)象，主要是由于燃面由低燃速藥柱向高燃速藥柱過渡時，明顯“膨大”所致。

(a) t=6.82 s (b) t=7.97 s

3 雙燃速藥柱交界面位置影響分析

雙燃速藥柱成型過程中，首先進行續(xù)航級藥柱澆注，待續(xù)航級藥柱預(yù)固化后，再澆注助推級藥柱。兩級藥柱的交界面一般通過續(xù)航級裝藥量進行控制，控制精度不高，本節(jié)通過參數(shù)研究，進行交界面位置對發(fā)動機內(nèi)彈道性能的影響評估，進而指導(dǎo)藥型優(yōu)化設(shè)計及藥柱成型要求。

圖5基于圖2所示藥柱將高低燃速交界面向兩側(cè)分別移動20 mm進行對比分析，推進劑性能參數(shù)保持不變，算例A高低燃速交界面位置為x=320 mm，算例B高低燃速交界面位置為x=280 mm。結(jié)果表明，算例A和算例B內(nèi)彈道特性與原藥柱存在顯著差異。算例A在轉(zhuǎn)級過程中，燃燒室壓強陡升至14.7 MPa，超出助推級初始高壓段的最大壓強。算例B在轉(zhuǎn)級過程中，燃燒室壓強大幅降低至0.58 MPa。

圖5 不同交界面位置下的壓強-時間曲線Fig.5 Pressure-time curves at different interface positions

圖6給出算例A在t=6.3 s時的燃面拓?fù)浣Y(jié)構(gòu)，此時低燃速藥柱推移肉厚為70 mm，與圖4(b)基本一致。對比可知，與原藥柱相比，算例A在從低燃速藥柱向高燃速藥柱過渡時，高燃速燃面 “膨大”現(xiàn)象更加明顯。與低燃速燃面疊加后燃面增加顯著，使得發(fā)動機壓強出現(xiàn)陡增，在6.3 s時發(fā)動機壓強達(dá)到14.7 MPa。若高低燃速的交界面繼續(xù)向右移動，發(fā)動機轉(zhuǎn)級過程產(chǎn)生的壓強峰會進一步增大，最終導(dǎo)致發(fā)動機殼體結(jié)構(gòu)失效。

圖6 算例A在6.3 s時的燃面拓?fù)浣Y(jié)構(gòu) 圖7 算例B在10 s時的燃面拓?fù)浣Y(jié)構(gòu)Fig.6 Topological structure of burning surface at t=6.3 s for Case A Fig.7 Topological structure of burning surface at t=10 s for Case B

如圖7所示，當(dāng)算例B燃面推移至高低燃速交界面時，已趨于平面，此時燃面面積最小且燃速較低，使得燃燒室內(nèi)的壓強降低至0.58 MPa，在此壓強下，推進劑不能維持正常燃燒，影響發(fā)動機工作可靠性。由以上分析可知，雙燃速藥柱交界面的位置直接影響著發(fā)動機內(nèi)彈道性能的穩(wěn)定性，因此需要在藥型參數(shù)設(shè)計及藥柱成型過程中嚴(yán)格控制。

4 結(jié)論

(1)與試驗結(jié)果對比可知，殘值函數(shù)法準(zhǔn)確捕捉了高低燃速推進劑交界面處燃面復(fù)雜拓?fù)浣Y(jié)構(gòu)的變化。

(2)雙燃速藥柱發(fā)動機內(nèi)彈道呈現(xiàn)顯著起伏現(xiàn)象，主要是由于燃面由低燃速藥柱向高燃速藥柱過渡時，明顯“膨大”所致。

(3)雙燃速藥柱交界面的位置直接影響著發(fā)動機內(nèi)彈道性能的穩(wěn)定性，需要在藥型參數(shù)設(shè)計及藥柱成型過程中嚴(yán)格控制。