費 陽，胡 凡，張為華，孫明波

(國防科技大學(xué)航天與材料工程學(xué)院，長沙 410073)

0 引言

裝藥燃面計算是固體火箭發(fā)動機(jī)內(nèi)彈道性能預(yù)示的重要內(nèi)容，其計算精度直接影響發(fā)動機(jī)性能預(yù)示的效果。

目前，已有多種裝藥燃面計算方法[1]，應(yīng)用較廣的有作圖法、通用坐標(biāo)法、實體造型法、解析法、網(wǎng)格推移法。作圖法需在圖紙上作出初始燃面形狀，繼續(xù)作出燃去肉厚后的燃面圖紙，使用測量器具測出燃面面積，該方法不適用于構(gòu)型復(fù)雜的藥型，且計算誤差較大;通用坐標(biāo)法通用性較好，但計算某些特殊裝藥構(gòu)型時，有燃面跳動，且不能計算含缺陷裝藥燃面;實體造型法目前被工業(yè)部門廣泛采用，充分利用CAD軟件，與藥型設(shè)計結(jié)合緊密且精度高，但對于結(jié)構(gòu)復(fù)雜的藥型，推移造型過程十分繁瑣，并可能出現(xiàn)奇點，計算無法繼續(xù);解析法建立在精確的燃面數(shù)學(xué)公式上，但只針對形面簡單的燃面計算;網(wǎng)格推移法在初始藥型上生成三角網(wǎng)格，再利用初始燃面按燃去肉厚生成新燃面，該方法通用性較好，但它處理燃面交匯、分離效果不佳。另外，以上方法無法計算含裂紋裝藥燃面。事實上，固體火箭在長時間儲存過程中，受到各種因素影響，裝藥可能產(chǎn)生裂紋，若裝藥存在裂紋，就將整個發(fā)動機(jī)報廢，勢必造成大量損失。實踐證明，有些含裂紋裝藥仍能滿足發(fā)動機(jī)任務(wù)要求。因此，對含裝藥裂紋固體發(fā)動機(jī)進(jìn)行性能預(yù)示是十分必要的，燃面計算是其中的關(guān)鍵與難點。

針對現(xiàn)有燃面計算方法的不足，本文應(yīng)用Orsher Stanely與 Sethian J A提出的 Level Set界面追蹤方法[1－3]，并結(jié)合孫明波提出的改進(jìn) Sub cell fix重新初始化新方法[4]進(jìn)行燃面計算。其優(yōu)點是可準(zhǔn)確捕捉界面運(yùn)動過程中形狀變化，自然處理界面拓?fù)浣Y(jié)構(gòu)變化(交錯、分離、消失)。但其缺點在于初始燃面定義需通過修改程序定義較為繁瑣，結(jié)合實體造型軟件進(jìn)行二次開發(fā)，可解決這一困難。另外，Level Set方法基于計算網(wǎng)格，計算量較大，計算速度與網(wǎng)格數(shù)量有關(guān)。對于主頻2.0 GHz的計算機(jī)，本文算例均耗時5 min左右?？赏ㄟ^并行計算技術(shù)，縮短計算時間。

1 基于Level Set方法的燃面計算模型

Level Set方法采取符號距離函數(shù)作為高維等值面，其核心思想是將N維曲面利用N+1維方程零等值面來表示。例如，二維平面圓x2+y2=1，可用三維曲面的z=0等值面來代替，z值便是二維平面上任一點到圓x2+y2=1的符號距離。

1.1 控制方程

根據(jù)Level Set方法數(shù)學(xué)思想，作如下定義:

設(shè)藥柱、包覆層內(nèi)部為區(qū)域Ω1，外部區(qū)域為Ω2，區(qū)域之間界面即為燃面，記為Γ(t)。記為t時刻為空間點到燃面的符號距離開始計算時，只需將設(shè)置區(qū)域Ω1內(nèi)φ值設(shè)為某一負(fù)值、區(qū)域Ω2內(nèi)φ值設(shè)為某一正值即可，經(jīng)過一次重新初始化，即可將φ值轉(zhuǎn)化為初始距離場，操作簡單。

Level Set法捕捉界面的思想是以適當(dāng)速度推動符號距離變化，得到每一點的符號距離，即可確定界面位置。對燃面法向推移控制方程式(2)，空間離散采用五階WENO格式，時間離散采用三階TVD-Runge-kutta格式。

其中，vn為燃面法向推移速度，可為空間、時間、壓強(qiáng)及流速等流場參數(shù)的函數(shù)。求解式(2)，可得各點φ值。事實上，只需準(zhǔn)確求出界面附近各點φ值便可。利用Narrow Band[1]方法，將求解區(qū)域限制在界面附近，可顯著減少計算量，時間復(fù)雜度為O(N)，N為3個坐標(biāo)方向上平均網(wǎng)格點數(shù)。

燃面面積如式(3)計算:

其中，δ為delta函數(shù)，標(biāo)記燃面點。

式中 ε =1.5Min(dx，dy，dz)。

1.2 重新初始化

因此，計算過程中每一時間步都對φ進(jìn)行重新初始化。采用國防科技大學(xué)孫明波提出的改進(jìn)Sub cell fix重新初始化格式，與傳統(tǒng)格式[5－7]相比，效果更佳、收斂速度更快[8]。

2 星形裝藥算例仿真

算例采用圖1所示星形裝藥，兩端包覆，內(nèi)表面燃燒。主要參數(shù)如表1所示。不考慮侵蝕燃燒效應(yīng)，利用Level set方法對其推移過程進(jìn)行仿真。

圖1 某星形裝藥構(gòu)型圖Fig.1 Structure of tested star grain

表1 算例裝藥結(jié)構(gòu)參數(shù)Table 1 Structure parameters of tested grain

2.1 算例1——裝藥無裂紋情況

與實體造型法進(jìn)行對比，驗證Level Set法的準(zhǔn)確性。計算網(wǎng)格數(shù)量為71×71×101，程序中將包覆層轉(zhuǎn)化為積分限，如果燃面超出包覆層，則不對其積分。圖2為裝藥燃面推移過程，星根部分逐漸消失，星尖不斷擴(kuò)展鈍化消失，最終趨于圓面，符合平行層推移規(guī)律?？梢?，Level Set法處理燃面交匯消失的能力較強(qiáng)。

圖3為Level set法與實體造型法燃面計算結(jié)果對比曲線。二者吻合較好，最大相對誤差為4.5%，平均相對誤差僅為1.3%。燃去肉厚達(dá)50 mm時，部分燃面抵達(dá)側(cè)面包覆層，燃面面積迅速下降。

圖2 無缺陷裝藥燃面推移過程Fig.2 Burning surface regression process of a perfect tested grain

圖3 燃面曲線對比Fig.3 Com parison of burning surface area lines

2.2 算例2——裝藥含單條裂紋情況

計算含單個裂紋時裝藥燃面，考察Level Set方法對復(fù)雜構(gòu)型裝藥燃面推移的處理能力。

裂紋位于裝藥正中，寬度3 mm，扇角120°，半徑95 mm。計算網(wǎng)格大小為71×71×101。圖4為含單裂紋裝藥燃面擴(kuò)展過程。

在圖4(c)中，可清楚看到裂紋頂部燃面消失，這是燃面“跑出”計算域的結(jié)果。實際上，計算域涵蓋整個裝藥，圖4(d)中方框內(nèi)部即為計算域。

圖5為裝藥含1條裂紋與無裂紋情況下燃面肉厚對比曲線。裂紋產(chǎn)生附加燃面，初始時燃面偏大，隨著燃燒過程的進(jìn)行，裂紋頂部燃面最先到達(dá)包覆層，使燃面增加速度在5 mm肉厚時放緩;裂紋擴(kuò)展產(chǎn)生附加燃面，同時裂紋到達(dá)包覆層，使燃面減少;二者作用不斷抵消，接近10 mm肉厚時，后者作用超過前者，使燃面增速低于無裂紋情況。燃去肉厚30 mm時，大部分裂紋燃面到達(dá)包覆，使總?cè)济娴陀跓o裂紋時的情況。

圖4 裝藥含1條裂紋時燃面推移過程Fig.4 Burining surface regression process of a one-cracked tested grain

圖5 裝藥含1條裂紋時燃面肉厚曲線Fig.5 Burning surface area of a one-cracked grain

2.3 算例3——裝藥含2條相鄰裂紋情況

利用Level Set方法計算裝藥含多裂紋時燃面肉厚曲線，考察Level Set方法對多裂紋復(fù)雜構(gòu)型裝藥燃面推移、交匯過程的處理能力。

裂紋尺寸與上節(jié)相同，位于裝藥中部，兩裂紋起始間距34 mm。計算網(wǎng)格大小為71×71×101。圖6為含兩相鄰裂紋裝藥等速平行推移過程。燃去肉厚接近17 mm時，兩裂紋互相交匯，形成1條更大的裂紋繼續(xù)擴(kuò)展下去。由此可見，Level Set方法可自然處理多裂紋擴(kuò)展、交匯等復(fù)雜結(jié)構(gòu)變化，效果較好。

圖7為裝藥含1條、2條裂紋與無裂紋情況下的燃面肉厚對比曲線。2條裂紋帶來更多的附加燃面，導(dǎo)致起始階段燃面面積比無裂紋情況時高很多。接近17 mm肉厚時，兩裂紋交匯，使附加燃面迅速下降。燃去肉厚25mm時，大部分裂紋燃面到達(dá)包覆，總?cè)济娴陀跓o裂紋情況。

圖6 裝藥含2條相鄰裂紋時燃面推移過程Fig.6 Burning surface regression process of a two-near-cracked tested grain

圖7 裝藥含2條相鄰裂紋時燃面肉厚曲線Fig.7 Burning surface area of a two-nearcracked grain

3 結(jié)論

跟蹤Level Set方法最新發(fā)展動態(tài)，采用全新的改進(jìn)Sub cell fix重新初始化方法，將Level Set界面追蹤方法與固體火箭發(fā)動機(jī)裝藥燃面計算相結(jié)合，利用Level Set方法，對固體火箭發(fā)動機(jī)裝藥燃面推移過程進(jìn)行了數(shù)值模擬，較準(zhǔn)確地捕捉燃面變化情況，并計算了裝藥燃面肉厚曲線，精度較高。采用Level Set方法，對裝藥含單條、多條裂紋情況進(jìn)行了燃面計算，準(zhǔn)確捕捉裂紋擴(kuò)展、交匯的過程?？梢姡琇evel Set方法對燃面交匯、消失等復(fù)雜拓?fù)浣Y(jié)構(gòu)變化處理能力較強(qiáng)。

實際上，固體發(fā)動機(jī)藥柱常見缺陷(裂紋、脫粘、氣泡)都可看作是裂紋。Level Set燃面計算方法是含缺陷裝藥燃面計算的有力工具。Level Set燃面計算程序易與發(fā)動機(jī)內(nèi)流場計算相結(jié)合，可實現(xiàn)三維內(nèi)彈道的模擬。

綜上所述，Level Set燃面計算方法有以下3個特點:

(1)可準(zhǔn)確計算復(fù)雜構(gòu)型裝藥燃面，精度高，通用性較好;

(2)自然處理推移過程中燃面交匯、消失等情況，無需特殊處理，應(yīng)用簡便;

(3)可對含裂紋等缺陷裝藥進(jìn)行燃面計算，彌補(bǔ)現(xiàn)有燃面計算方法的不足。

[1] 馬長禮.固體火箭發(fā)動機(jī)燃面計算方法研究[D].國防科技大學(xué)，2007.

[2] Osher S，F(xiàn)edkiw R.Level setmethods and dynamic implicit surfaces[D].Spring-Verlag，New York，2003.

[3] 秦飛.固體火箭發(fā)動機(jī)復(fù)雜裝藥燃面算法研究[D].西北工業(yè)大學(xué)，2003.

[4] 孫明波.超聲速來流穩(wěn)焰凹腔的流動及火焰穩(wěn)定機(jī)制研究[D].國防科技大學(xué)，2008.

[5] Sussman M，Smereka P，Osher S.A level set approach for computing solutions to incompressible two-phase flow[J].J.Comp.Phys.，1994，114:146-159.

[6] Russo G，Smereka P.A remark on computing distance functions[J].Journal of Computational Physics，2000，163:51-67.

[7] Hartmann D，Meinke M，SchroderW.Erratum to“differential equation based constrained reinitialization for level set methods”[J].Journal of Computational Physics，2008，227(22):9682-9696.

[8] Sun Ming-bo，Wang Zheng-guo，Bai Xue-song.Assessment and modification of sub-cell-fixmethod for re-initialization of level set distance function[J].International Journal for Numerical Methods in Fluids，2010，62:211-236.