基于殘值函數法的雙燃速藥柱內彈道分析*

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

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

0 引言

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

目前，可實現燃面不等速推移的有實體造型法[2]、網格推移法[3-4]、最小距離函數法[5-6]、Level Set方法[7-8]和殘值函數法[9]等。其中，實體造型法通過特征造型或驅動尺寸實現燃面推移過程的模擬。但是，對于結構復雜的藥型，實體造型法的推移過程十分繁瑣，并可能出現奇點，導致計算無法繼續；對于多燃速藥柱，燃面的分離需要在推移過程中構建新的幾何特征，較難實現。網格推移法在初始藥型上生成非結構網格，利用當地燃速進行推移，再通過網格重構獲得新的燃面，該方法通用性較好，但其處理燃面交匯、分離等拓撲結構變化的穩定性較差。最小距離函數法通過計算藥柱內部各點到初始燃面的距離，即最小距離函數(MDF)，進行燃面推移分析，可以實現不等速推移，但是其無法自然捕捉到由于燃速間斷的存在而導致的燃面拓撲結構的變化。Level Set方法采用初值形式的偏微分方程將一個純幾何問題轉變為用偏微分方程描述的數學問題，對于推進劑內部存在燃速間斷的工況，由于燃氣與推進劑之間也存在間斷，在數值求解微分方程組時，會遇到間斷相交的情況，進而導致燃面計算精度及穩定性較差。殘值函數法[9]在笛卡爾網格上，采用殘值函數記錄燃面位置，通過惠更斯原理模擬燃面的推移過程。此算法可以準確捕捉由于燃速間斷的存在導致燃面的交匯、分離、消失等復雜拓撲結構變化，且不需要進行網格重構，穩定性好。

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

1 內彈道計算方法

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

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

第一步：如圖1所示，Pi為當前燃面上的網格點。預估燃面上的離散點Pi的影響區域，并向四周推移，在推移的區域內，推進劑轉化為燃氣，推移距離Δi由式(1)給出。

Δi=dt×ri+δi

(1)

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

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

(2)

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

(1)如圖1(a)所示，若Pi與Pj之間沒有燃速間斷，當Pi與Pj之間的距離dij趨于0時，相應的燃速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點的直線距離。

第三步：確定新的燃面，即更新后的燃氣區域與藥柱區域的交界面。

結合內彈道基本方程，可以計算燃燒室內壓強，獲得發動機內彈道性能。

2 某雙燃速藥柱內彈道分析

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

表1 推進劑性能參數

圖2 藥柱結構圖Fig.2 Structure diagram of the grain

圖3給出采用殘值函數法獲得的發動機內彈道數據，與試驗結果進行對比可知，殘值函數法準確地模擬了雙燃速藥柱燃燒過程燃面的變化過程。

圖3 數值仿真與試驗結果對比Fig.3 Comparison of results of numerical simulation and experiment

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

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

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

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

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

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

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

圖6 算例A在6.3 s時的燃面拓撲結構 圖7 算例B在10 s時的燃面拓撲結構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所示，當算例B燃面推移至高低燃速交界面時，已趨于平面，此時燃面面積最小且燃速較低，使得燃燒室內的壓強降低至0.58 MPa，在此壓強下，推進劑不能維持正常燃燒，影響發動機工作可靠性。由以上分析可知，雙燃速藥柱交界面的位置直接影響著發動機內彈道性能的穩定性，因此需要在藥型參數設計及藥柱成型過程中嚴格控制。

4 結論

(1)與試驗結果對比可知，殘值函數法準確捕捉了高低燃速推進劑交界面處燃面復雜拓撲結構的變化。

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

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