基于平行層推移的含表觀裂紋缺陷固體發動機裝藥燃面計算①

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

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

0 引言

裝藥燃面計算是固體火箭發動機內彈道性能預示的重要內容，其計算精度直接影響發動機性能預示的效果。

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

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

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

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

1.1 控制方程

根據Level Set方法數學思想，作如下定義:

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

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

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

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

其中，δ為delta函數，標記燃面點。

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

1.2 重新初始化

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

2 星形裝藥算例仿真

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

圖1 某星形裝藥構型圖Fig.1 Structure of tested star grain

表1 算例裝藥結構參數Table 1 Structure parameters of tested grain

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

與實體造型法進行對比，驗證Level Set法的準確性。計算網格數量為71×71×101，程序中將包覆層轉化為積分限，如果燃面超出包覆層，則不對其積分。圖2為裝藥燃面推移過程，星根部分逐漸消失，星尖不斷擴展鈍化消失，最終趨于圓面，符合平行層推移規律。可見，Level Set法處理燃面交匯消失的能力較強。

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

圖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方法對復雜構型裝藥燃面推移的處理能力。

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

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

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

圖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方法對多裂紋復雜構型裝藥燃面推移、交匯過程的處理能力。

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

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

圖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 結論

跟蹤Level Set方法最新發展動態，采用全新的改進Sub cell fix重新初始化方法，將Level Set界面追蹤方法與固體火箭發動機裝藥燃面計算相結合，利用Level Set方法，對固體火箭發動機裝藥燃面推移過程進行了數值模擬，較準確地捕捉燃面變化情況，并計算了裝藥燃面肉厚曲線，精度較高。采用Level Set方法，對裝藥含單條、多條裂紋情況進行了燃面計算，準確捕捉裂紋擴展、交匯的過程。可見，Level Set方法對燃面交匯、消失等復雜拓撲結構變化處理能力較強。

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

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

(1)可準確計算復雜構型裝藥燃面，精度高，通用性較好;

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

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

[1] 馬長禮.固體火箭發動機燃面計算方法研究[D].國防科技大學，2007.

[2] Osher S，Fedkiw R.Level setmethods and dynamic implicit surfaces[D].Spring-Verlag，New York，2003.

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

[4] 孫明波.超聲速來流穩焰凹腔的流動及火焰穩定機制研究[D].國防科技大學，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.