固體火箭發動機裝藥燃面計算方法

褚佑彪，張　崗，苑　博，武　淵

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

?

固體火箭發動機裝藥燃面計算方法

褚佑彪1，張崗2，苑博1，武淵1

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

通過引入殘值函數記錄燃面位置，在笛卡爾離散網格上運用惠更斯原理進行燃面顯式推移，實現了固體火箭發動機三維復雜藥柱的燃面推移模擬和燃面計算功能；采用此算法分別對具有復雜幾何構型的藥柱和由不同燃速特性的推進劑構成的藥柱進行燃面推移模擬。結果表明，此算法精度高、穩定性好，可準確捕捉燃面交匯、分離、消失等復雜拓撲結構變化，適用于復雜裝藥的燃面計算。

燃面推移；殘值函數；笛卡爾離散網格；惠更斯原理

0　引言

燃面計算用于確定裝藥在燃燒過程中燃燒表面積隨燃燒時間的變化規律，直接影響發動機內彈道性能預示精度，是發動機內彈道設計的基礎，在固體火箭發動機的設計中一直占有重要地位[1]。傳統的燃面計算方法包括解析法、作圖法及通用坐標法，均因各自的缺點較明顯已鮮有應用。隨著計算機相關技術，特別是大型三維工程設計軟件的發展，實體造型法應運而生[2-3]。目前，實體造型法主要包括特征造型法和驅動尺寸法。盡管實體造型法在工程中已得到廣泛應用，并已成為藥型設計的主要手段，但這種基于商用軟件的燃面計算程序仍有很大的局限性，比如對于某些特殊幾何形狀的藥型，在CAD軟件的推移過程中會出現奇異點，使得燃面推移不能繼續下去，只能通過構造近似的幾何形狀進行燃面計算。

為滿足先進固體動力技術的發展要求，藥型設計也越來越復雜，對燃面模擬方法不僅提出了適應復雜幾何構型的要求，還提出可處理燃面不等速推移的要求，包括處理不同燃速的推進劑在其交界面處同時燃燒的情況，比如單室雙推發動機。鑒于實體造型法的局限性，一些新型的燃面計算方法相繼出現，其中可實現復雜燃面不等速推移的有網格法[4-5]、Level Set方法[6-7]和最小距離函數法[8-9]等。早期的網格法僅適用于軸對稱和二維藥柱[3]。2005年，沈偉和邢耀國根據惠更斯原理，構建一種在通用CFD軟件非結構網格系統上直接計算燃面推移的數值方法，可以進行三維藥柱的燃面分析，但是算法的穩定性有待進一步研究[5]。Level Set方法采用初值形式的偏微分方程將一個純幾何問題轉變為用偏微分方程描述的數學問題，但由于在燃面推移過程中需解微分方程組,計算量非常大[6-7，10]。最小距離函數法通過計算藥柱內部各點到初始燃面的距離，即最小距離函數(MDF)，進行燃面推移分析，可以實現不等速推移[9]，但是無法處理不同燃速的推進劑在交界面處同時燃燒的情況。

針對現有燃面模擬方法的不足，本文采用笛卡爾網格離散藥柱，通過引入殘值函數記錄燃面位置，運用惠更斯原理進行燃面顯式推移，實現了三維復雜藥柱的燃面推移模擬和燃燒面積計算功能。其基本思想是采用笛卡爾網格和殘值函數描述當前燃面，通過燃面上球面波影響區域的疊加,模擬燃面的推移過程。此算法的優點是在無需進行網格重構的前提下，自然處理界面拓撲結構變化，準確模擬燃面的推移過程，且計算量較小。

1　實現過程

藥柱燃面推移的實現過程簡單描述如圖1所示。

圖1　燃面推移流程圖

1.1初始化

為了便于工程應用，同時充分利用CAD軟件的強大的實體造型功能，本文首先在CAD軟件中建立藥柱的幾何模型，然后將藥柱的面信息(燃面或限燃面)保存為STL格式的文件，文件中記錄曲面小三角平面的頂點坐標及其外法向向量。

工程中藥柱表面一般為曲面，本文為了降低笛卡爾網格捕捉曲面的數值誤差，引入殘值函數δi，用于描述真實燃面的幾何特征。初始時刻，燃面附近區域采用網格點到初始燃面的最小距離函數作為殘值函數，根據網格離散點相對于三角平面的空間位置[9]，進行點到面、點到線、點到點3種情況分類計算。

1.2燃面推移

燃面推移模擬的理論基礎是惠更斯原理，如圖2所示，燃面可以理解為“波陣面”，面上的任一點都可以看作是新的“次波源”，由這些“次波源”發出的“波”所形成的包絡面，就是藥柱退移到的新燃面，“波”的傳播速度對應著推進劑的燃速。由此可知，運用惠更斯原理可實現燃面的不等速推移，燃面推移模擬的一般性原理。

圖2　惠更斯傳播原理示意圖

圖3　燃面推移示意圖

因此，本文采用惠更斯原理實現燃面的推移，并通過殘值函數記錄當前燃面的位置，如圖3所示。燃面的推移流程可簡單描述如下：

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

(1)

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

第二步：計算并更新殘值函數。如圖4所示，Pi為當前燃面上的網格點，Pj為藥柱內的網格點且處于Pi的影響范圍內。Pj處的殘值函數δi由式(2)獲得：

(2)

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

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

(3)

獲得。

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

(4)

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

所以，Pj處的εij可由式(5)獲得：

(5)

式中dio為Pi到O點的直線距離；doj為O到Pj點的直線距離。

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

(a)無燃速間斷　　　　(b)有燃速間斷

1.3燃面計算

積分獲得當前時間步每種推進劑的燃氣體積增量。當前時間步內每種推進劑的平均燃面值為其燃氣體積增量與推進距離之商。當單一推進劑的推進距離具有空間不均勻性時，可采用平均距離代替，這是因為，用于內彈道預示的物理量是燃氣的體積增量[1]，即燃面與推進距離的乘積，而不是單一的燃面或推進距離。

1.4內彈道計算，確定燃速

由上述燃面計算過程及內彈道學基本方程可知，在給定燃速ri的情況下，可獲得對應的燃面，進而獲得對應的燃燒室壓強p。而在已知p的前提下，ri可由p直接確定，兩者之間一般滿足：

(6)

式中ai和ni分別為Pi點處推進劑的燃速系數ai和燃速的壓強指數ni。

根據推進劑燃速ri與燃燒室壓強p的這兩種關系，在數值求解過程中，可在推進劑的燃速ri與燃燒室壓強p之間建立循環，直至壓強的殘差小于設定的預示誤差pres，即認為燃燒室內壓強的求解過程完成。當前時間步內的燃速由式(6)確定，并可作為下一時間步內燃速的預估值。

2　算例

2.1具有復雜幾何構型的藥柱

為驗證上述算法的正確性，本節選取圖5所示具有復雜幾何構型的φ400 mm藥柱作為驗證算例。圖6給出計算獲得的燃面的變化過程，并與基于Pro/E平臺的驅動尺寸方法[2]的結果進行對比。采用主頻為2.9 GHz 的計算機進行單線程計算，計算結果如表1所示。基于網格1和網格2的計算機耗時分別為27 s和96 s，誤差分別不超過1.2%和0.5%。由此可知，此算法精度較高，計算量較小，滿足工程實際應用。

圖5　三維藥柱結構圖

表1　不同計算網格的計算結果對比

圖6　燃燒面積隨時間的變化規律

2.2具有不同燃速特性的推進劑構成的藥柱

在單室雙推發動機的設計過程中，為保證助推段與續航段的推力比，常采用高低燃速搭配的思路進行

助推段和續航段裝藥設計。當燃面與高低燃速推進劑交界面相交時，推進劑的燃燒速度在燃面上出現強間斷，呈現明顯的差異性，是燃面退移計算的一個難點。

如圖7所示，本節以內孔燃燒藥柱作為數值驗證算例，圖中給出不同時刻藥柱過軸的剖面圖，灰色區域為高燃速推進劑，黑色區域為低燃速推進劑，白色區域為燃氣填充區域。為充分驗證當前算法準確處理燃面上存在燃速間斷工況的能力，模擬了3個藥柱的燃面推移過程，其交界面分別與初始燃面成45°、90°和135°夾角。

由圖7可知，在交界面處，藥柱幾何特征的演化過程均被準確捕捉。驗證結果表明，當前算法可準確捕捉由于燃速間斷的存在導致燃面的交匯、分離、消失等復雜拓撲結構變化，適用于復雜裝藥的燃面計算。

(a)推進劑交界面與初始燃面成45°夾角　　　　　　　　(b)推進劑交界面與初始燃面成90°夾角

(c)推進劑交界面與初始燃面成135°夾角

3　結論

(1)此算法基于燃面推移的一般性原理——惠更斯原理，不僅可處理燃面的等速推移，還可處理燃面的不等速推移，包括燃面上存在燃速間斷的工況。

(2)采用笛卡爾網格離散藥柱，對當前燃面的影響區域可有效地提前預估，將三維體循環縮減為具有一定厚度的曲面循環，大大降低了計算量。

(3)殘值函數的引入，不僅有效地控制燃面推移的累計誤差，準確計算燃面，還可自然處理界面拓撲結構變化，避免燃面的網格重構，提高算法的穩定性。

[1]陳汝訓,劉銘初,李志明.固體火箭發動機設計與研究(上)[M].北京:中國宇航出版社,1991.

[2]董新剛,陳林泉,侯曉.基于Pro/E平臺下的固發裝藥CAD軟件[C]//2002年中國宇航學會固體推進專業委員會年會論文集(上),昆明:2002:109-114.[3]于勝春,趙汝巖,周紅梅，等.基于Pro/E特征造型技術的固體發動機裝藥燃面計算[J].固體火箭技術,2005,28(2):108-111.

[4]Hejl R J,Heister S D.Solid rocket motor grain burnback analysis using adaptive grids[J].Journal of Propulsion and Power,1995,11(5):1006-1011.

[5]沈偉,邢耀國.基于非結構網格的燃面推進算法[J].固體火箭技術,2005,28(3):176-179.

[6]Sethian J A.Theory,algorithms,and applications of level set methods for propagating interfaces[J].Acta Numerica,1995,5(5):309-395.

[7]秦飛.固體火箭發動機復雜裝藥燃面算法研究[D].西安：西北工業大學,2003.

[8]Willcox M A,Brewster M Q,Tang K C,et al.Solid propellant grain design and burnback simulation using a minimum distance function[J].Journal of Propulsion Power,2007,23(2):465-475.

[9]馬長禮.固體火箭發動機MDF燃面計算方法研究[D].長沙：國防科學技術大學,2007.

[10]費陽,胡凡,張為華,等.基于平行層推移的含表觀裂紋缺陷固體發動機裝藥燃面計算[J].固體火箭技術,2011,34(4):466-469.

(編輯：呂耀輝)

Burning surface computing method for solid rocket motor grain

CHU You-biao1，ZHANG Gang2，YUAN Bo1，WU Yuan1

(1.The National Key Laboratory of Combustion,Thermostructure and Flow of SRM,The 41st Institute of the Fourth Academy of CASC,Xi'an710025,China；2.The Fourth Academy of CASC,Xi'an710025,China)

A new computational method, based on a Cartesian grid and Huygens' principle of wave propagation,is presented for simulating the evolution of the burning surface regression of a complex,three-dimensional solid rocket motor propellant grain.The residual displacement function (RDF) is introduced to describe the location of complex burning surface on the Cartesian grid.Burning surface calculation of grain with complicated geometry or multi-burning rate was performed.The results indicate that the method has high precision and good stability,and it is capable of dealing with complicated structure changes such as burning surface intersection,separation and vanishing.It is demonstrated that the method provides a new technical measure for burning surface calculation of complicated propellant grain.

burning surface regression;residual displacement function;Cartesian grid;Huygens' principle of wave propagation

2015-08-31；

2015-11-05。

褚佑彪(1987—)，男，博士，研究方向為固體火箭發動機設計。E-mail：cyber@mail.ustc.edu.cn

V435

A

1006-2793(2016)04-0488-04

10.7673/j.issn.1006-2793.2016.04.007