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

(1.中國航天科技集團(tuán)公司四院四十一所,固體火箭發(fā)動機(jī)燃燒、熱結(jié)構(gòu)與內(nèi)流場國防科技重點(diǎn)實(shí)驗(yàn)室,西安　710025；2.中國航天科技集團(tuán)公司四院，西安　710025)

?

固體火箭發(fā)動機(jī)裝藥燃面計(jì)算方法

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

(1.中國航天科技集團(tuán)公司四院四十一所,固體火箭發(fā)動機(jī)燃燒、熱結(jié)構(gòu)與內(nèi)流場國防科技重點(diǎn)實(shí)驗(yàn)室,西安710025；2.中國航天科技集團(tuán)公司四院，西安710025)

通過引入殘值函數(shù)記錄燃面位置，在笛卡爾離散網(wǎng)格上運(yùn)用惠更斯原理進(jìn)行燃面顯式推移，實(shí)現(xiàn)了固體火箭發(fā)動機(jī)三維復(fù)雜藥柱的燃面推移模擬和燃面計(jì)算功能；采用此算法分別對具有復(fù)雜幾何構(gòu)型的藥柱和由不同燃速特性的推進(jìn)劑構(gòu)成的藥柱進(jìn)行燃面推移模擬。結(jié)果表明，此算法精度高、穩(wěn)定性好，可準(zhǔn)確捕捉燃面交匯、分離、消失等復(fù)雜拓?fù)浣Y(jié)構(gòu)變化，適用于復(fù)雜裝藥的燃面計(jì)算。

燃面推移；殘值函數(shù)；笛卡爾離散網(wǎng)格；惠更斯原理

0　引言

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

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

針對現(xiàn)有燃面模擬方法的不足，本文采用笛卡爾網(wǎng)格離散藥柱，通過引入殘值函數(shù)記錄燃面位置，運(yùn)用惠更斯原理進(jìn)行燃面顯式推移，實(shí)現(xiàn)了三維復(fù)雜藥柱的燃面推移模擬和燃燒面積計(jì)算功能。其基本思想是采用笛卡爾網(wǎng)格和殘值函數(shù)描述當(dāng)前燃面，通過燃面上球面波影響區(qū)域的疊加,模擬燃面的推移過程。此算法的優(yōu)點(diǎn)是在無需進(jìn)行網(wǎng)格重構(gòu)的前提下，自然處理界面拓?fù)浣Y(jié)構(gòu)變化，準(zhǔn)確模擬燃面的推移過程，且計(jì)算量較小。

1　實(shí)現(xiàn)過程

藥柱燃面推移的實(shí)現(xiàn)過程簡單描述如圖1所示。

圖1　燃面推移流程圖

1.1初始化

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

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

1.2燃面推移

燃面推移模擬的理論基礎(chǔ)是惠更斯原理，如圖2所示，燃面可以理解為“波陣面”，面上的任一點(diǎn)都可以看作是新的“次波源”，由這些“次波源”發(fā)出的“波”所形成的包絡(luò)面，就是藥柱退移到的新燃面，“波”的傳播速度對應(yīng)著推進(jìn)劑的燃速。由此可知，運(yùn)用惠更斯原理可實(shí)現(xiàn)燃面的不等速推移，燃面推移模擬的一般性原理。

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

圖3　燃面推移示意圖

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

第一步：預(yù)估燃面上的離散點(diǎn)Pi的影響區(qū)域，并向四周推移，在推移的區(qū)域內(nèi)，推進(jìn)劑轉(zhuǎn)化為燃?xì)猓埔莆灰痞由式(1)給出。

(1)

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

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

(2)

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

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

(3)

獲得。

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

(4)

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

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

(5)

式中dio為Pi到O點(diǎn)的直線距離；doj為O到Pj點(diǎn)的直線距離。

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

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

1.3燃面計(jì)算

積分獲得當(dāng)前時(shí)間步每種推進(jìn)劑的燃?xì)怏w積增量。當(dāng)前時(shí)間步內(nèi)每種推進(jìn)劑的平均燃面值為其燃?xì)怏w積增量與推進(jìn)距離之商。當(dāng)單一推進(jìn)劑的推進(jìn)距離具有空間不均勻性時(shí)，可采用平均距離代替，這是因?yàn)椋糜趦?nèi)彈道預(yù)示的物理量是燃?xì)獾捏w積增量[1]，即燃面與推進(jìn)距離的乘積，而不是單一的燃面或推進(jìn)距離。

1.4內(nèi)彈道計(jì)算，確定燃速

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

(6)

式中ai和ni分別為Pi點(diǎn)處推進(jìn)劑的燃速系數(shù)ai和燃速的壓強(qiáng)指數(shù)ni。

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

2　算例

2.1具有復(fù)雜幾何構(gòu)型的藥柱

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

圖5　三維藥柱結(jié)構(gòu)圖

表1　不同計(jì)算網(wǎng)格的計(jì)算結(jié)果對比

圖6　燃燒面積隨時(shí)間的變化規(guī)律

2.2具有不同燃速特性的推進(jìn)劑構(gòu)成的藥柱

在單室雙推發(fā)動機(jī)的設(shè)計(jì)過程中，為保證助推段與續(xù)航段的推力比，常采用高低燃速搭配的思路進(jìn)行

助推段和續(xù)航段裝藥設(shè)計(jì)。當(dāng)燃面與高低燃速推進(jìn)劑交界面相交時(shí)，推進(jìn)劑的燃燒速度在燃面上出現(xiàn)強(qiáng)間斷，呈現(xiàn)明顯的差異性，是燃面退移計(jì)算的一個(gè)難點(diǎn)。

如圖7所示，本節(jié)以內(nèi)孔燃燒藥柱作為數(shù)值驗(yàn)證算例，圖中給出不同時(shí)刻藥柱過軸的剖面圖，灰色區(qū)域?yàn)楦呷妓偻七M(jìn)劑，黑色區(qū)域?yàn)榈腿妓偻七M(jìn)劑，白色區(qū)域?yàn)槿細(xì)馓畛鋮^(qū)域。為充分驗(yàn)證當(dāng)前算法準(zhǔn)確處理燃面上存在燃速間斷工況的能力，模擬了3個(gè)藥柱的燃面推移過程，其交界面分別與初始燃面成45°、90°和135°夾角。

由圖7可知，在交界面處，藥柱幾何特征的演化過程均被準(zhǔn)確捕捉。驗(yàn)證結(jié)果表明，當(dāng)前算法可準(zhǔn)確捕捉由于燃速間斷的存在導(dǎo)致燃面的交匯、分離、消失等復(fù)雜拓?fù)浣Y(jié)構(gòu)變化，適用于復(fù)雜裝藥的燃面計(jì)算。

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

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

3　結(jié)論

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

(2)采用笛卡爾網(wǎng)格離散藥柱，對當(dāng)前燃面的影響區(qū)域可有效地提前預(yù)估，將三維體循環(huán)縮減為具有一定厚度的曲面循環(huán)，大大降低了計(jì)算量。

(3)殘值函數(shù)的引入，不僅有效地控制燃面推移的累計(jì)誤差，準(zhǔn)確計(jì)算燃面，還可自然處理界面拓?fù)浣Y(jié)構(gòu)變化，避免燃面的網(wǎng)格重構(gòu)，提高算法的穩(wěn)定性。

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

[2]董新剛,陳林泉,侯曉.基于Pro/E平臺下的固發(fā)裝藥CAD軟件[C]//2002年中國宇航學(xué)會固體推進(jìn)專業(yè)委員會年會論文集(上),昆明:2002:109-114.[3]于勝春,趙汝巖,周紅梅，等.基于Pro/E特征造型技術(shù)的固體發(fā)動機(jī)裝藥燃面計(jì)算[J].固體火箭技術(shù),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]沈偉,邢耀國.基于非結(jié)構(gòu)網(wǎng)格的燃面推進(jìn)算法[J].固體火箭技術(shù),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]秦飛.固體火箭發(fā)動機(jī)復(fù)雜裝藥燃面算法研究[D].西安：西北工業(yè)大學(xué),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]馬長禮.固體火箭發(fā)動機(jī)MDF燃面計(jì)算方法研究[D].長沙：國防科學(xué)技術(shù)大學(xué),2007.

[10]費(fèi)陽,胡凡,張為華,等.基于平行層推移的含表觀裂紋缺陷固體發(fā)動機(jī)裝藥燃面計(jì)算[J].固體火箭技術(shù),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—)，男，博士，研究方向?yàn)楣腆w火箭發(fā)動機(jī)設(shè)計(jì)。E-mail：cyber@mail.ustc.edu.cn

V435

A

1006-2793(2016)04-0488-04

10.7673/j.issn.1006-2793.2016.04.007