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

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

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

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

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

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

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

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

0 引言

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

目前國內外燃面退移實現方法主要有以下5種：

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

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

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

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

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

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

1 數值方法與計算方案

1.1 計算方法

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

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

主要流程如下：

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

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

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

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

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

步驟6 重復步驟1～步驟5，直到計算結束。

1.2Level-set方法

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

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

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

傳統level-set方程形式為

(1)

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

(2)

(2)式為由非穩態項、對流項和源項組成的標量方程。

1.3MDF方法

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

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

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

(3)

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

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

1.4 界面分割和幾何重構

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

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

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

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

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

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

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

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

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

1.5Level-set方程求解

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

1.5.1 Level-set初場設置

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

(4)

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

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

選取這樣的水平集函數主要有以下優點：

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

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

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

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

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

2 界面退移計算驗證

2.1 裝藥燃面退移分析

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

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

2)燃面凸面退移，如星形藥柱的和車輪形藥柱等，都存在裝藥部分凸向燃燒的燃面，在退移過程中，界面當地曲率不斷增大，最后出現凸尖點，甚至其后退移過程一直存在凸尖點，如圖8(b)和圖8(c)所示；

3)燃面凹面退移，如星形藥柱和車輪形藥柱等存在內燃面的藥柱都存在凹向燃燒室的燃面，且初始當地曲率都比較大(極限情況為凹尖點)，在退移過程中，當地曲率不斷減小，形成曲率比較小的曲面，如圖8(d)和圖8(e)所示；

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

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

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

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

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

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

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

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

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

2.2 星形裝藥燃面退移

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

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

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

提取各時刻垂直于裝藥軸線的橫截面，可以看到LSM和MDF方法組合技術計算得到的燃面退移過程，如圖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函數場進行處理，采用1.4節描述的方法計算得到燃面面積S，如圖11所示，將其與建模軟件Pro/E輸出的面積值相比，其面積誤差εS曲線如圖12所示。從圖12中可以看到：Level-set方法計算燃面退移得到的燃面面積大小和建模軟件輸出的面積值是基本一致的，且誤差最大不超過5%.

2.3 復雜三維裝藥燃面退移

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

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

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

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

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

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

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

在此使用如圖16所示的內外燃管形裝藥來驗證Level-set方法和MDF方法組合技術是否能用來計算變燃速燃面退移，圖16中藥柱區域分為四部分，各個區域燃速各不相同，燃速之比分別為1.05∶1.10∶1.00∶1.15.

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

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

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

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

圖18 變燃速裝藥燃面面積對比曲線(燃層厚度以燃速最低的區域計算)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 侵蝕燃燒燃面退移計算

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

侵蝕比ε通常被總結為關于喉通比J、密流G或燃通比等參數的經驗公式，一般可以表示為如下形式：

(5)

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

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

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

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

(6)

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

提取各個時刻的燃面面積，繪制曲線，與無侵蝕燃燒工況相比較，如圖22所示。從圖22中可以看到：1)在裝藥燃燒的前期，侵蝕燃燒使得裝藥燃面面積比無侵蝕工況要大，主要原因是侵蝕燃燒使得裝藥后段從圓柱面變為了錐面；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

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

3 結論

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

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

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

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

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]. 固體火箭技術,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] 周華盛.藥柱通用坐標計算法計算結果跳動的原因分析及解決途徑[J]. 固體火箭技術,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] 鮑福廷,李逢春.固體發動機裝藥CAD[J].固體火箭技術,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] 田維平,王錕.固體發動機藥柱CAD及燃燒模擬分析[J]. 固體火箭技術,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] 方蜀州,胡克嫻.固體火箭發動機三維藥柱燃面退移仿真技術及燃面通用計算方[J]. 固體火箭技術, 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] 顏仙榮,于勝春，王慶官.實體造型技術與固體發動機裝藥燃面計算[J].固體火箭技術,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平臺的固體火箭發動機內彈道性能計算[D].哈爾濱：哈爾濱工程大學, 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] 沈偉,邢耀國.基于非結構網格的燃面推進算法[J]. 固體火箭技術, 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] 賀征,郜冶.動態網格在固體火箭發動機內流場計算中的應用研究[J].彈箭與制導技術, 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] 秦飛.固體火箭發動機復雜裝藥燃面算法研究[D]. 西安：西北工業大學,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] 何濤,孫振華.星型裝藥固體火箭發動機工作過程仿真研究[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] 馬長禮.固體火箭發動機MDF燃面計算方法研究[D].長沙：國防科學技術大學, 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]. 航空學報, 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] 林小樹, 王寶山. 雙燃速固體火箭發動機內彈道計算方法[J]. 固體火箭技術, 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—), 男, 教授, 博士生導師。 E-mail: wangge@hrbeu.edu.cn

V435+.12

A

1000-1093(2017)02-0280-12

10.3969/j.issn.1000-1093.2017.02.011