基于水平集方法和最小距離函數(shù)法的復(fù)雜裝藥燃面退移問題研究

王革， 韓萬之， 李冬冬， 郜冶

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院, 黑龍江 哈爾濱 150001)

基于水平集方法和最小距離函數(shù)法的復(fù)雜裝藥燃面退移問題研究

王革， 韓萬之， 李冬冬， 郜冶

(哈爾濱工程大學(xué) 航天與建筑工程學(xué)院, 黑龍江 哈爾濱 150001)

利用水平集(Level-set)方法和最小距離函數(shù)(MDF)方法，對固體火箭發(fā)動機裝藥燃面退移問題進行研究?；贚evel-set方法和MDF方法的相互耦合，建立了一種可以準確預(yù)測固體推進劑裝藥燃燒表面退移的方法；對不同幾何藥型、不同燃速組合推進劑裝藥及帶侵蝕效應(yīng)的裝藥進行計算分析，其結(jié)果表明，該組合方法及所采用的網(wǎng)格界面分割和幾何重構(gòu)計算燃面面積的思想，在預(yù)測固體推進劑裝藥燃面瞬態(tài)退移問題上具有良好的適應(yīng)性和較高的可信度。

兵器科學(xué)與技術(shù)； 裝藥燃面退移； 燃面面積； 水平集方法； 最小距離函數(shù)方法

0 引言

準確的固體火箭發(fā)動機內(nèi)彈道計算對于預(yù)測固體火箭發(fā)動機工作過程有著重要意義，而裝藥燃面退移過程直接決定著其內(nèi)彈道變化特性。

目前國內(nèi)外燃面退移實現(xiàn)方法主要有以下5種：

1)通用坐標法。該方法于1968年，由Peterso等[1]提出，主要思想是：將初始空腔分解為基本幾何元素，通過幾何元素的擴大，實現(xiàn)燃燒表面退移。1975年，Nickson等[2]將通用坐標法應(yīng)用到內(nèi)彈道計算中，設(shè)計開發(fā)了固體火箭發(fā)動機性能程序(SPP)，經(jīng)過后來一系列改進，成為了固體火箭發(fā)動機性能評估的標準。侯曉等[3]和周華盛[4]通過對通用坐標法計算燃面結(jié)果跳動的分析，分別采用直接求解燃面和通過將肉厚等份數(shù)設(shè)置為積分區(qū)間長度的倍數(shù)減小體積分誤差的方法，消除了計算結(jié)果的波動。鮑福廷等[5]在通用坐標法的基礎(chǔ)上，結(jié)合計算機圖形處理直接獲得幾何數(shù)據(jù)，研制了三維裝藥設(shè)計程序，并建立了一套圖形庫。這種處理方式不依賴于復(fù)雜的求交運算，可以獲得更為完備的燃燒室內(nèi)通道參數(shù)，為一維內(nèi)彈道計算及內(nèi)流場計算提供詳盡的數(shù)據(jù)。

2)實體造型法。隨著幾何造型技術(shù)的發(fā)展，出現(xiàn)了很多依托計算機輔助設(shè)計CAD軟件平臺，進行二次開發(fā)的方法，可以得到任意時刻燃面和與藥柱體積相關(guān)的信息。1993年，田維平等[6]借助計算機輔助設(shè)計與輔助制造I-DEAS軟件和實體造型方法設(shè)計開發(fā)了藥柱CAD燃燒模擬分析系統(tǒng)——CAPBAS. 1993年，方蜀州等[7]利用實體造型法思想，借助AutoCAD11.0的繪圖功能，在其高級造型擴展軟件包AME的基礎(chǔ)上開發(fā)了用于燃面計算的Regress3D程序。2003年，顏仙榮等[8]以計算機輔助設(shè)計與輔助制造UG軟件為開發(fā)平臺借助UG的實體造型功能，完成了固體火箭發(fā)動機裝藥的設(shè)計、燃面退移模擬和燃面計算，與實驗結(jié)果的對比表明該方法具有較高的精度。2014年，劉偉[9]借助建模軟件Pro/E，使用Qt平臺進行軟件界面的開發(fā)，實現(xiàn)了發(fā)動機藥型設(shè)計與內(nèi)彈道計算，并且可與通用計算流體力學(xué)CFD軟件進行數(shù)據(jù)交流，實現(xiàn)藥型設(shè)計、內(nèi)彈道計算與流場計算一體化。

3)動網(wǎng)格方法。1995年，Hejl等[10]，用自適應(yīng)網(wǎng)格方法解決軸對稱和二維裝藥燃面計算問題，實現(xiàn)了不等速燃燒燃面計算。1998年，Breton等[11]用非結(jié)構(gòu)網(wǎng)格方法進行二維裝藥內(nèi)彈道計算，通過求解漢密爾頓一雅克比方程，解決不等速燃燒燃面計算問題。2014年，Gueyffier等[12]使用水平集(Level-set)方法和傅里葉變換等多種復(fù)雜運算及高階算法計算界面節(jié)點移動，并配合使用體網(wǎng)格適應(yīng)技術(shù)避免網(wǎng)格的重新劃分。用來計算復(fù)雜藥型的內(nèi)流場，得到了比較好的結(jié)果。2005年，沈偉等[13]根據(jù)惠更斯波傳播原理，構(gòu)建了一種在CFD軟件非結(jié)構(gòu)網(wǎng)格系統(tǒng)上直接計算燃面推進的數(shù)值方法并實現(xiàn)了燃面推進。2006年，Jiao[14]提出了一種新的面網(wǎng)格移動技術(shù)——face offsetting。該技術(shù)提供了一個動網(wǎng)格界面移動的統(tǒng)一框架，即使燃面存在奇點和大曲率半徑也可以得到比較精準的結(jié)果。2008年，賀征等[15]使用動態(tài)層方法對星形裝藥燃面退移及內(nèi)流場進行了一體化計算，計算結(jié)果和準穩(wěn)態(tài)計算結(jié)果以及實驗數(shù)據(jù)的對比表明，動網(wǎng)格方法可以得到優(yōu)于準穩(wěn)態(tài)的結(jié)果。2013年，LI等[16]使用了face offsetting方法進行藥柱燃面退移計算，配合網(wǎng)格光順和網(wǎng)格重構(gòu)技術(shù)處理流體域網(wǎng)格變形和扭曲，模擬了三維復(fù)雜藥柱燃面退移內(nèi)流場，結(jié)果和實驗吻合較好。

4)界面追蹤法。界面追蹤法將發(fā)動機燃面視為流場中的自由邊界，通過一定手段追蹤或者捕獲燃面位置，然后求得有關(guān)燃面面積和裝藥體積。其中用的最多的是Level-set方法。1996年，Sethian等[17]利用LSM，建立運動界面的控制方程，通過求解實現(xiàn)了對自由界面的追蹤。2001年，HEGAB等[18]發(fā)展了一種異質(zhì)推進劑的非穩(wěn)態(tài)燃燒模型，同時模擬氣相的燃燒和固相的傳熱，并且在界面處施加階躍。在模擬界面退移時使用了Level-set方程推導(dǎo)出的漢密爾頓-雅克比方程，并且通過Level的不同標識出高氯酸銨(AP)和黏結(jié)劑。2003年，秦飛[19]將Level-set方法用于燃面計算。將裝藥燃面看作不同物質(zhì)的分界面，建立了燃面控制方程，通過求解控制方程計算出燃面變化規(guī)律，并對含缺陷裝藥和變?nèi)妓傺b藥進行了計算。2005年，Yild等[20]用LSM分別對二維和三維裝藥下不等燃速燃燒情況進行計算。2011年，Cavallini等[21]用準一維非穩(wěn)態(tài)內(nèi)彈道耦合Level-set界面退移計算了發(fā)動機從開始點火到推進劑燃燒殆盡整個過程，該方法特別適合計算侵蝕燃燒下發(fā)動機的內(nèi)彈道，且可以給出比零維內(nèi)彈道更準確的壓力變化曲線。2014年，Wang等[22]使用Level-set方法和零維內(nèi)彈道耦合計算發(fā)動機燃面退移，并將其應(yīng)用于固體火箭發(fā)動機的裝藥設(shè)計。2014年，Hwang等[23]發(fā)展了3種使用固定網(wǎng)格的界面追蹤法：界面追蹤、射線法和最小距離函數(shù)(MDF)法。界面追蹤是基于拉格朗日視角的方法，而射線法和MDF法是基于歐拉視角的方法。使用這3種方法進行模擬，發(fā)現(xiàn)MDF法比其他兩種方法要更準確。2015年，何濤等[24]以FLUENT為平臺，利用UDF進行二次開發(fā)，完成了星型裝藥燃面點火和燃面退移的建模， 研究了星型裝藥固體火箭發(fā)動機點火和燃面退移過程。

5)MDF方法。2005年，Willcox等[25]提出了MDF燃面計算方法，處理燃面計算普遍存在的通用性和穩(wěn)定性問題。2007年，馬長禮[26]對MDF迭代過程進行改進，提高了MDF方法燃面計算的精度，通過大量算例驗證了MDF方法的通用性，對不等燃速燃面計算也進行了嘗試。2009年，熊文波等[27]在單元法中，使用了網(wǎng)格到燃面的最小距離判斷某一時刻單元內(nèi)裝藥是否燃盡，通過相鄰肉厚的體積微分獲得燃面面積。這種方法避免了大量的數(shù)學(xué)推導(dǎo)，也避免了內(nèi)型面標準幾何體的復(fù)雜定義。

國內(nèi)外研究工作表明，界面追蹤方法已經(jīng)比較好地運用到固體火箭發(fā)動機裝藥的燃面退移計算中，其中Level-set方法得到了比較好的結(jié)果。但是Level-set使用和計算也存在著一些問題，通常使用的都是漢密爾頓- 雅克比形式的Level-set方程，而這種方法在退移界面的曲率過大或過小時，求解存在問題，計算存在凸尖點的界面退移時需要特殊處理。本文主要結(jié)合MDF燃面計算方法，對Level-set方法進行改進，使其在燃面退移的計算上更加方便和準確，發(fā)展精確的燃面面積計算方法，為燃面退移和流場的耦合計算提供良好的基礎(chǔ)。

1 數(shù)值方法與計算方案

1.1 計算方法

采用Level-set方法和MDF方法組合技術(shù)計算燃面退移。Level-set方法求解需要初場和燃面(界面)退移速度場，求解后可以得到當前燃面位置和流體與固體(簡稱流固)區(qū)域劃分；MDF方法初始化需要燃面位置和流固區(qū)域劃分信息，變換后可以為Level-set方程的求解提供所需的初場和燃面退移速度場，這兩種方法可以互相配合，其整個過程如圖1所示。

圖1 燃面退移計算方案Fig.1 Calculation scheme of burning surface regression

主要流程如下：

步驟1 根據(jù)計算域中燃面位置和區(qū)域劃分(流體區(qū)域和固體區(qū)域)初始化帶符號最小距離函數(shù)Φe；

步驟2 變換Φe得到Level-set方程求解所需的水平集函數(shù)初場φt和界面退移速度場Vφ；

步驟3 求解Level-set方程，得到t+Δt時刻的水平集函數(shù)場φt+Δt；

步驟4 使用t+Δt時刻水平集函數(shù)場φt+Δt中的零等值面對計算域中的所有網(wǎng)格進行網(wǎng)格分割和幾何重構(gòu)，得到燃面位置信息，計算燃面面積和固體體積分數(shù)；

步驟5 根據(jù)t+Δt時刻的水平集函數(shù)場φt+Δt確定計算域中流體區(qū)域和固體區(qū)域的劃分；

步驟6 重復(fù)步驟1～步驟5，直到計算結(jié)束。

1.2Level-set方法

使用Level-set方法求解燃面退移問題，主要就是利用一個水平集函數(shù)φ來表征出界面在不同時刻的位置，通常這樣的函數(shù)φ需要滿足兩個條件：1)能夠表征出流體區(qū)域和固體區(qū)域；2)能夠表征出不同時刻界面的位置。

在此使用φ>0部分表征固體區(qū)域，φ<0部分表征外部流體區(qū)域，φ=0表征界面，如圖2所示。

圖2 Level-set水平集函數(shù)示意圖Fig.2 Schematic diagram of Level-set function

傳統(tǒng)level-set方程形式為

(1)

(1)式為非守恒形式，將其轉(zhuǎn)換為守恒形式

(2)

(2)式為由非穩(wěn)態(tài)項、對流項和源項組成的標量方程。

1.3MDF方法

MDF方法基本思想是：計算藥柱內(nèi)部各點到初始燃面的位置，即最小距離函數(shù)Φe，按照平行層退移規(guī)律，選取Φe的值等于已燃厚度e的點，即可組成已燃厚度e時的燃面，如圖3所示，Φe=0是初始燃面位置，Φe=Δe是燃層厚度為Δe時的燃面位置。

圖3 MDF方法示意圖Fig.3 Schematic diagram of MDF

將MDF在整個計算域上按照如下方式擴展得到帶符號的最小距離函：

(3)

式中：d為計算域中某點與初始燃面的距離。

可以說最小距離法可以獨立地完成按照平行層退移規(guī)律的燃面計算，而且經(jīng)過其他學(xué)者證明，精度和穩(wěn)定性都很高。當裝藥燃速大小在空間存在差異時，MDF方法難以用來計算燃面退移，但是可以在燃面附近為燃面提供一個合理的界面退移方向，即MDF的法向。

1.4 界面分割和幾何重構(gòu)

t時刻的水平集函數(shù)場φt經(jīng)過Level-set方程的求解，得到t+Δt時刻的水平集函數(shù)場φt+Δt. 根據(jù)φt+Δt中的零等值線，即界面，對計算域進行分割，在此使用基于等值面的方法進行界面分割，這種方法可以使界面在網(wǎng)格交界面處連續(xù)。為了方便計算界面，在網(wǎng)格內(nèi)簡化為一線段(二維)或平面(三維)，以二維為例，如圖4所示。

圖4 基于等值線方法的界面分割示意圖Fig.4 Interface segmentation based on contour

界面分割之后，需要對重構(gòu)的網(wǎng)格進行幾何運算，求解網(wǎng)格中的界面面積和固體體積分數(shù)，以二維四邊形網(wǎng)格為例，如圖5所示，流程如下：

1)使用泰勒展開求得網(wǎng)格頂點處的水平集函數(shù)φ，如圖5(a)中的網(wǎng)格頂點A、B、C和D的水平集函數(shù)φA、φB、φC和φD；

2)使用零等值線分割網(wǎng)格，如圖5(b)中零等值線交網(wǎng)格界面于點E和點F，其中EF即為網(wǎng)格內(nèi)的界面面積；

3)進行幾何重構(gòu)，排除φ<0的流體區(qū)域，剩余固體區(qū)域，如圖5(c)中剩余多邊形ABEFD；

4)對多邊形ABEFD的頂點進行順時針或者逆時針排序，如圖5(d)中得到頂點序列P1、P2、P3、P4和P5；

圖5 三角分割求解界面面積和體積分數(shù)Fig.5 Triangular segmentation for solving interface area and volume fraction

這種網(wǎng)格界面分割和幾何重構(gòu)的方法，適合各種二維凸多邊形網(wǎng)格或者三維凸多面體網(wǎng)格，適用性好，方便編程處理；另外，對網(wǎng)格進行界面分割時，可以將重構(gòu)后得到的退移界面上的節(jié)點存儲起來，為帶符號的距離函數(shù)重新初始化提供所需的界面信息。

1.5Level-set方程求解

Level-set方程的求解需要為其指定初場φt和界面退移速度場Vφ.

1.5.1 Level-set初場設(shè)置

對帶符號的距離函數(shù)進行如下變換可以得到一個水平集函數(shù)φ，如圖6所示。

(4)

圖6 水平集函數(shù)φ示意圖Fig.6 Level-set function

式中：L表示燃面附近的一條窄帶的寬度，通常選取網(wǎng)格平均尺寸的4～5倍距離。

選取這樣的水平集函數(shù)主要有以下優(yōu)點：

1)能夠簡單的標識出流體區(qū)域、藥柱區(qū)域及燃面，φ<0的區(qū)域為流體區(qū)域，φ>0的區(qū)域為藥柱區(qū)域，φ=0的等值面為燃面；

2)在界面附近有1階線性分布的水平集函數(shù)，方便確定界面的準確位置；

3)在遠離界面的位置處φ值方便確定，在藥柱區(qū)域內(nèi)φ=1.0，在流體區(qū)域內(nèi)φ=-1.0，便于給定邊界值，或可以直接使用源項固定部分區(qū)域的值作為內(nèi)邊界。

1.5.2 Level-set速度場Vφ的確定

圖7 界面退移示意圖Fig.7 Interface regression

2 界面退移計算驗證

2.1 裝藥燃面退移分析

推進劑藥柱根據(jù)燃燒表面的位置可以分為：端面燃燒藥柱、側(cè)面燃燒藥柱和側(cè)端面同時燃燒藥柱。側(cè)面燃燒藥柱可分為：外燃藥柱、內(nèi)燃藥柱和內(nèi)外燃藥柱。不管藥柱如何復(fù)雜，根據(jù)平行層燃燒規(guī)律，再考慮其他因素的影響，燃面退移主要包含以下6種方式，如圖8所示。

1)燃面平面退移，如端面燃燒藥柱，在退移過程中，燃面始終保持為一平面，如圖8(a)所示；

2)燃面凸面退移，如星形藥柱的和車輪形藥柱等，都存在裝藥部分凸向燃燒的燃面，在退移過程中，界面當?shù)厍什粩嘣龃螅詈蟪霈F(xiàn)凸尖點，甚至其后退移過程一直存在凸尖點，如圖8(b)和圖8(c)所示；

3)燃面凹面退移，如星形藥柱和車輪形藥柱等存在內(nèi)燃面的藥柱都存在凹向燃燒室的燃面，且初始當?shù)厍识急容^大(極限情況為凹尖點)，在退移過程中，當?shù)厍什粩鄿p小，形成曲率比較小的曲面，如圖8(d)和圖8(e)所示；

4)燃面部分燒完，燃面到達包覆層等位置，燃面與包覆層相交，如圖8(f)所示；

5)裝藥可能是由不同燃速推進劑組成的，燃面以不同燃燒速度退移，如圖8(g)所示；

6)裝藥部分燃燒表面出現(xiàn)侵蝕燃燒，下游燃速比上游高，如圖8(h)所示。

圖8 裝藥燃面退移示意圖Fig.8 Schematic diagram of grain burning regression

為了驗證Level-set方法和MDF方法組合技術(shù)是否能夠用來計算上述燃面退移，采用以下4種計算模型進行驗證：

1)三維星形裝藥，星形裝藥燃面在退移過程中，退移形式包含了圖8(a)～圖8(f)的6種形式；

2) 復(fù)雜三維裝藥燃面退移，退移形式包含了圖8(a)～圖8(f)的6種形式；

3) 不同燃速推進劑組成的單管內(nèi)外燃裝藥燃面退移，可以用來檢驗上述中的退移形式圖8(g)；

4) 發(fā)生侵蝕燃燒的單管內(nèi)外燃裝藥燃面退移，可以用來檢驗上述中的退移形式圖8(h)。

2.2 星形裝藥燃面退移

圖9是計算采用的星形裝藥模型，裝藥幾何參數(shù)如表1所示，裝藥兩端面包覆不燃燒，星形界面按照平行層規(guī)律退移。

圖9 星形裝藥模型Fig.9 Star grain modle

參數(shù)數(shù)值外徑D/m0.269藥柱長度L/m0.937藥柱肉厚e1/m0.0565星角數(shù)n6星邊夾角θ/(°)64.18角度系數(shù)ε0.88過渡圓弧半徑r/m0.008星角圓弧半徑r1/m0.008

提取各時刻垂直于裝藥軸線的橫截面，可以看到LSM和MDF方法組合技術(shù)計算得到的燃面退移過程，如圖10所示。從圖10中可以看到：

1)星尖從凸曲面退移為凸尖點，且在初期星尖的角度不變，后期逐漸變大；

2)槽部的過渡圓弧逐漸變大；

3)星邊初期一直保持為平面；

4)燃燒后期部分裝藥燒完，燃面與包覆層相交。

圖10 星形裝藥燃面變化Fig.10 Change of burning surface in star grain

圖11 星形裝藥燃面面積變化Fig.11 Change of burning surface in star grain

圖12 燃面面積誤差Fig.12 Deviation of burning area

對Level-set方法計算得到的Level-set函數(shù)場進行處理，采用1.4節(jié)描述的方法計算得到燃面面積S，如圖11所示，將其與建模軟件Pro/E輸出的面積值相比，其面積誤差εS曲線如圖12所示。從圖12中可以看到：Level-set方法計算燃面退移得到的燃面面積大小和建模軟件輸出的面積值是基本一致的，且誤差最大不超過5%.

2.3 復(fù)雜三維裝藥燃面退移

圖13為計算所采用的復(fù)雜三維藥柱結(jié)構(gòu)示意圖，此藥柱部分被包覆，部分裸露于燃氣中，屬于內(nèi)外燃藥柱，在此計算此藥柱的燃面退移變化規(guī)律，檢驗在裝藥幾何比較復(fù)雜的情況下Level-set和MDF組合方法計算燃面退移的準確性。

圖13 復(fù)雜裝藥幾何模型Fig.13 Geometrical model of complex grain

提取燃面退移過程中的燃面幾何，如圖14所示。從圖14中可以看出藥柱燃面的退移過程：1)初始時在藥柱的內(nèi)側(cè)、兩端及部分外側(cè)都有燃面，此時燃面面積較大；2)隨著藥柱尾部燃燒殆盡，藥柱變成近似短管形三維裝藥。

圖14 裝藥燃面變化示意圖Fig.14 Change of burning surface

計算Level-set和MDF組合方法得到的燃面面積，如圖15所示。從圖15中可以看出Level-set方法計算藥柱退移的準確性：使用該方法計算得到的燃面面積與文獻[26]給出的燃面變化曲線基本一致，說明Level-set和MDF組合方法可以準確地計算三維復(fù)雜裝藥的燃面退移。

圖15 燃面面積變化曲線Fig.15 Curves of burning surface area

2.4 不同燃速推進劑組成裝藥的燃面退移計算

在此使用如圖16所示的內(nèi)外燃管形裝藥來驗證Level-set方法和MDF方法組合技術(shù)是否能用來計算變?nèi)妓偃济嫱艘?，圖16中藥柱區(qū)域分為四部分，各個區(qū)域燃速各不相同，燃速之比分別為1.05∶1.10∶1.00∶1.15.

圖16 不同燃速裝藥示意圖Fig.16 Schematic diagram of grains with different burning rates

提取燃面退移過程中的燃面幾何，如圖17所示，圖中Z、R分別為藥柱長度和半徑。從圖17中可以看出燃面的變化過程：1)各區(qū)域處的燃面分別以各自對應(yīng)的燃速退移，維持局部燃面的半徑一致；2)在各區(qū)域交界處，由于兩側(cè)燃速不同，出現(xiàn)近似階梯狀曲面；3)燃速大的區(qū)域推進劑燃燒完后，總的燃面驟然減少，藥柱可能存在脫落未燃盡的碎塊。

計算Level-set和MDF組合方法求得的燃面面積，與燃速一致的裝藥燃面面積及文獻[28]中提出的方法計算結(jié)果進行對比，如圖18所示。圖19為本文方法與文獻中解析解的誤差曲線。從圖19中可以看到燃面面積的變化過程：1)在隨著各區(qū)域裝藥燃燒完，燃面面積出現(xiàn)了多次明顯降低，使得發(fā)動機工作末段與預(yù)期工作狀態(tài)相差很大。如果不同批次推進劑的燃速存在明顯差異將對發(fā)動機工作過程產(chǎn)生較大影響；2)在裝藥燃燒的前期，Level-set方法與解析解的方法相比誤差在0.5%以下，在后期由于各區(qū)域推進劑先后燃盡，誤差有所上升，但最大誤差不超過5%.

圖17 不同燃速裝藥燃面變化Fig.17 Change of burning surfaces of grains with different burning rates

圖18 變?nèi)妓傺b藥燃面面積對比曲線(燃層厚度以燃速最低的區(qū)域計算)Fig.18 Comparison of burning surface areas with different rates (calculation of burning layer is based on minimum burning rate)

圖19 燃面面積誤差Fig.19 Deviation of burning area

2.5 侵蝕燃燒燃面退移計算

侵蝕燃燒效應(yīng)指的是平行于裝藥燃燒表面的氣體流動狀態(tài)對裝藥燃速的影響。侵蝕燃燒效應(yīng)使裝藥燃速增大，燃氣生成率也隨之增大，從而使燃燒室內(nèi)燃氣壓強升高。在發(fā)動機工作初始階段，裝藥通道截面積最小，氣流速度大，侵蝕燃燒效應(yīng)最嚴重。

侵蝕比ε通常被總結(jié)為關(guān)于喉通比J、密流G或燃通比等參數(shù)的經(jīng)驗公式，一般可以表示為如下形式：

(5)

式中：B表示喉通比J、密流G或者燃通比等參數(shù)；Bth表示發(fā)生侵蝕燃燒時參數(shù)B的臨界值；K表示經(jīng)驗系數(shù)。

計算采用如圖20所示的內(nèi)外燃單管形裝藥，裝藥兩端包覆，燃面為內(nèi)外圓柱面，假定侵蝕燃燒只發(fā)生在燃通比比較大的管內(nèi)燃氣通道下游，且只發(fā)生在前期。侵蝕燃燒的計算需要和流場計算耦合，在此為了方便設(shè)置，在此并不進行與流場耦合的侵蝕燃燒燃面退移計算，而是根據(jù)侵蝕燃燒的特點，只進行類似于侵蝕燃燒形式的燃面退移計算。

圖20 侵蝕燃燒模型示意圖Fig.20 Schematic diagram of erosion combustion model

在此假定無侵蝕燃速r0為0.1 m/s，侵蝕比為采用如下形式：

(6)

提取各個時刻的燃面幾何圖，如圖21所示。從圖21中可以看到與實際侵蝕燃燒相類似的燃面變化趨勢：1)管內(nèi)燃氣通道從某處開始一直到末端，因為侵蝕燃燒的影響，退移速度比上游大，且越靠近末端退移速度越大；2)裝藥后段比前段燃面退移速度大，形成一個近似錐面的曲面；3)內(nèi)外燃面相遇，裝藥后段先燒完。

提取各個時刻的燃面面積，繪制曲線，與無侵蝕燃燒工況相比較，如圖22所示。從圖22中可以看到：1)在裝藥燃燒的前期，侵蝕燃燒使得裝藥燃面面積比無侵蝕工況要大，主要原因是侵蝕燃燒使得裝藥后段從圓柱面變?yōu)榱隋F面；2)在裝藥燃燒的后期，侵蝕燃燒使得裝藥燃面面積比無侵蝕工況小，主要原因是侵蝕燃燒使得裝藥后段提前燒完。

圖21 侵蝕燃燒燃面變化圖Fig.21 Change of burning surface with erosion combustion

圖22 侵蝕燃燒燃面面積對比曲線Fig.22 Comparison of burning surface areas in erosion combustion and normal combustion

綜合侵蝕燃燒對燃速和裝藥燃面面積的影響，可以發(fā)現(xiàn)侵蝕燃燒會使發(fā)動機工作的初期壓力明顯高于無侵蝕時的壓力預(yù)測值，研究侵蝕燃燒對發(fā)動機內(nèi)彈道的影響有重大意義。

3 結(jié)論

通過上述對Level-set方法、MDF方法及其耦合算法和界面面積計算方法的介紹，對不同形式和不同情況下的裝藥進行計算分析可以得到以下結(jié)論：

1)Level-set方法和MDF方法組合技術(shù)可以比較準確地計算各種形式的燃面退移，包括平面燃面退移、凸曲面(凸尖點)燃面退移、凹曲面(凹尖點)燃面退移、裝藥部分包覆燃面退移、燃面以不同燃速退移和侵蝕燃燒燃面退移。

2)采用的面積計算方式可以比較準確地計算燃面面積，為瞬態(tài)內(nèi)流場的計算提供了一個良好的基礎(chǔ)。

目前本文燃面的耦合計算可以實現(xiàn)，并對不同結(jié)構(gòu)的推進劑裝藥燃面退移過程進行有效的預(yù)測，但是由于所采用的MDF方法計算量較大，還需要對MDF計算進行一些技術(shù)上的改進，以期更快速、高效地進行燃面的耦合計算。

References)

[1] Peterson E G, Nielsen G C. Generalized coordinate grain design and internal ballistics evaluation program[C]∥Solid Propulsion Conference. Atlantic City, US: AIAA,2013:68-490.

[2] Coats D E, Levine J N, Nickerson G R. A computer program for the prediction of solid propellant rocket motor performance [R]. California: Ultra Systems Inc,1975.

[3] 侯曉,蹇澤群.三維藥柱燃面的通用積分計算法[J]. 固體火箭技術(shù),1993,9(3):1-6. HOU Xiao, JIAN Ze-qun. General integration calculation for burning surface of three dimensional grain[J]. Journal of Solid Rocket Technology, 1993, 9(3):1-6. (in Chinese)

[4] 周華盛.藥柱通用坐標計算法計算結(jié)果跳動的原因分析及解決途徑[J]. 固體火箭技術(shù),1994,9(3):8-16. ZHOU Hua-sheng. Analysis and solution approach about pulsation cause of calculation results of the general coordinate calculation method of the grain[J]. Journal of Solid Rocket Technology, 1994,9(3):8-16. (in Chinese)

[5] 鮑福廷,李逢春.固體發(fā)動機裝藥CAD[J].固體火箭技術(shù),1994,9(3):1-7. BAO Fu-ting, LI Feng-chun. Computer aided design of propellant grains for solid rocket motors[J]. Journal of Solid Rocket Technology, 1994, 9(3):1-7. (in Chinese)

[6] 田維平,王錕.固體發(fā)動機藥柱CAD及燃燒模擬分析[J]. 固體火箭技術(shù),1993,12(3):14-22. TIAN Wei-ping, WANG Kun. Computer aided design and burning simulation analysis of solid rocket motor grains[J]. Journal of Solid Rocket Technology, 1993,12(3):14-22. (in Chinese)

[7] 方蜀州,胡克嫻.固體火箭發(fā)動機三維藥柱燃面退移仿真技術(shù)及燃面通用計算方[J]. 固體火箭技術(shù), 1993,12(4):10-20. FANG Shu-zhou, HU Ke-xian. The general technique of emulation and calculating of burning surface of 3D grain of solid rocket motors[J]. Journal of Solid Rocket Technology, 1993,12(4):10-20. (in Chinese)

[8] 顏仙榮,于勝春，王慶官.實體造型技術(shù)與固體發(fā)動機裝藥燃面計算[J].固體火箭技術(shù),2003,26(2):20-22. YAN Xian-rong, YU Sheng-chun, WANG Qing-guan. The solid modeling technology and the grain-burning-area calculation[J]. Journal of Solid Rocket Technology, 2003, 26(2): 20-22. (in Chinese)

[9] 劉偉.基于Pro/E和Qt平臺的固體火箭發(fā)動機內(nèi)彈道性能計算[D].哈爾濱：哈爾濱工程大學(xué), 2014. LIU Wei. Performance computing of solid propellant rocket motor[D]. Harbin: Harbin Engineering University,2014. (in Chinese)

[10] Voller V R, Brent A D. The modeling of heat, mass and solute transport in solidification systems[J]. International Journal of Heat and Mass Transfer, 1989, 32(9) : 1719-1731.

[11] Breton P L, Ribéreau D, Godfroy F. SRM performance analysis by coupling bidimensional surface burnback and pressure field computations[C]∥34th AIAA JPC Conference. Cleveland, US:AIAA, 1998.

[12] Gueyffier D, Roux F X. High-order computation of burning propellant surface and simulation of fluid flow in solid rocket chamber[C]∥50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference.Cleveland, US: AIAA, 2014.

[13] 沈偉,邢耀國.基于非結(jié)構(gòu)網(wǎng)格的燃面推進算法[J]. 固體火箭技術(shù), 2005, 28 (3):176-179. SHEN Wei, XING Yao-guo. A simulation method of burning surface regression based on unstructured mesh[J]. Journal of Solid Rocket Technology, 2005, 28 (3):176-179. (in Chinese)

[14] Jiao X M. Face offsetting: a unified approach for explicit moving interfaces[J]. Journal of Computational Physics, 2007, 220(2):612-625.

[15] 賀征,郜冶.動態(tài)網(wǎng)格在固體火箭發(fā)動機內(nèi)流場計算中的應(yīng)用研究[J].彈箭與制導(dǎo)技術(shù), 2008,28(1):164-166. HE Zheng, GAO Ye. The study of dynamic grids application to calculate the inner flow field of solid rocket motor[J]. Journal of Projectiles, Rockets, Missiles and Guidance, 2008, 28(1):164-166. (in Chinese)

[16] Li Q, He G Q. Coupled simulation of fluid flow and propellant burning surface regression in a solid rocket motor[J]. Computers & Fluids, 2014, 93(8):146-152.

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

[18] Hegab A, Jackson T L, Buckmaster J, et al. Nonsteady burning of periodic sandwich propellants with complete coupling between the solid and gas phases[J]. Combustion and Flame, 2001, 125(1): 1055-1070.

[19] 秦飛.固體火箭發(fā)動機復(fù)雜裝藥燃面算法研究[D]. 西安：西北工業(yè)大學(xué),2003. QIN Fei. Method research for burning surface calculation of solid rocket motor with complicated grain[D]. Xi’an: Northwestern Polytechnical University, 2003. (in Chinese)

[20] Yildirim C, Aksel M H. Numerical simulation of the grain burnback in solid propellant rocket motor[C]∥41st AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. Arizona, US: AIAA,2005:3995-4160.

[21] Cavallini E, Favini B, Di Giacinto M, et al. Internal ballistics simulation of a NAWC tactical SRM[J]. Journal of Applied Mechanics, 2011, 78(5): 748-760.

[22] Wang D H,Yang F, Fan H, et al. An integrated framework for solid rocket motor grain design optimization[J]. Proceedings of the Institution of Mechanical Engineers, Part G: Journal of Aerospace Engineering, 2014, 228(7): 1156-1170.

[23] Yao H H, Chiang C H. Tracking methods to study the surface regression of the solid-propellant grain[J]. International Journal of Engineering and Technology Innovation, 2014, 4(4):213-222.

[24] 何濤,孫振華.星型裝藥固體火箭發(fā)動機工作過程仿真研究[J].航空兵器,2015(1):61-64. HE Tao, SUN Zhen-hua. Numerical simulation of the working process in star grain solid rocket motor[J]. Aero Weaponry, 2015(1):61-64. (in Chinese)

[25] Michael A Willcox. Solid rocket motor internal ballistics simulation using three-dimensional grain burnback[J]. Journal of Propulsion and Power, 2007, 23(3):575-584.

[26] 馬長禮.固體火箭發(fā)動機MDF燃面計算方法研究[D].長沙：國防科學(xué)技術(shù)大學(xué), 2007. MA Chang-li. Research of MDF burning surface CaleUlation method for solid rocket motor[D].Changsha : National University of Defense Technology, 2007. (in Chinese)

[27] 熊文波,劉宇. 基于單元法的三維裝藥通用燃面計算[J]. 航空學(xué)報, 2009, 30(7):1176-1180. XIONG Wen-bo, LIU Yu. Generalized burning surface calculation of three dimensional propellant based on element method[J]. Acta Aeronautica et Astronautica Sinica, 2009, 30(7):1176-1180. (in Chinese)

[28] 林小樹, 王寶山. 雙燃速固體火箭發(fā)動機內(nèi)彈道計算方法[J]. 固體火箭技術(shù), 1991 (4): 12-18. LIN Xiao-shu, WANG Bao-shan. Calculation method of solid rocket motor internal ballistic with dual burning rate[J]. Journal of Solid Rocket Technology, 1991(4):12-18. (in Chinese)

Research on Grain Burning Surface Regression Based on Level-set Method and Minimum Distance Function

WANG Ge, HAN Wan-zhi, LI Dong-dong, GAO Ye

(College of Architecture and Aerospace Engineering, Harbin Engineering University, Harbin 150001, Heilongjiang， China)

The combustion surface regression is studied using Level-set method and minimum distance function (MDF) in SRM. A way that accurately predicts burning surface regression of solid propellant is established based on Level-set method and MDF. The grains with different geometries, the propellant grains with different burning rates and the grains with erosion combustion are calculated and analyzed. The results show that the combined method which calculates the burning area based on interface segmentation and geometric reconstruction has a good adaptability and high reliability on combustion surface regression.

ordnance science and technology; combustion surface regression; combustion surface area; level-set method; MDF method

2016-07-06

王革(1966—), 男, 教授, 博士生導(dǎo)師。 E-mail: wangge@hrbeu.edu.cn

V435+.12

A

1000-1093(2017)02-0280-12

10.3969/j.issn.1000-1093.2017.02.011