激波繞射動態裂紋的數值模擬*

羅 雷，孫振生，張世英

(第二炮兵工程大學，西安 710025)

0 引言

固體火箭發動機在生產、儲存以及使用過程中，裝藥表面或內部不可避免的形成裂紋。在點火沖擊波的作用下，這些裂紋有可能變形、擴展，進而會導致藥柱結構的破壞，甚至造成發動機爆燃等惡性事故。國內外學者在點火沖擊波與裂紋的相互作用以及裂紋的擴展等方面做了大量的工作，認為點火階段壓力梯度的迅速升高是導致裂紋擴展的重要因素［1-4］。點火過程中，激波在裂紋內的繞射和反射是壓力梯度變化的重要原因，但對這一非定常過程的研究還非常少。因此，研究激波在裂紋中的傳播過程，揭示點火沖擊波與裝藥裂紋相互作用的規律和機理，對發動機判廢準則的制定具有重要意義。

激波繞射裂紋計算的一個主要難點是對邊界的處理。文獻［5］在假設裂紋壁面靜止的情況下對強激波繞射裂紋所形成的流場進行了詳細分析，得出了許多有價值的結論。但是隨著推進劑的燃燒，計算區域是不斷變化的。假設裂紋壁面靜止對流場結構的影響是無法估計的，而引入移動邊界進行計算可有效地解決這一問題。

LevelSet作為一種有效的移動界面追蹤方法，已廣泛應用于多個領域［6］。文中利用高精度加權本質無振蕩(weighted essentially non-oscillatory，WENO)格式［7］同Level Set方法相結合，考慮了由于燃燒造成裂紋邊界的移動，對激波繞射裂紋形成的非定常流場進行了計算，并對這種激波與移動邊界相互作用的非定常流場進行了詳細分析，探討了點火期間裂紋內壓力突升的主要機理。并通過對靜態和動態邊界條件下所形成的流場結構進行比較，探討了引入動態邊界條件計算的必要性。

1 控制方程與數值計算方法

1.1 流場控制方程

二維非定常可壓縮Navier-Stokes方程可寫成如下守恒形式:

其中:U=(ρ，ρu，ρv，e)T;E、F 為無粘通量，Ev、Fv為粘性通量。為了使方程封閉，補充完全氣體方程p=(γ－1)ρe。流場控制方程的求解采用高精度WENO格式，具體的離散方法可參考文獻［7］。

1.2 界面的追蹤

Level Set的控制方程如下:

其中，Level Set函數φ如下定義:

V是界面的運動速度，定義Level Set函數的單位外法向N=?φ/|?φ|，那么界面運動的速度V可以由式V=rbN給出。其中，rb為由于裂紋表面燃燒造成的邊界點推移速度，用rb=αpn計算。其中，p為該處的靜壓，α與n均為與推進劑有關的常數。

1.3 Ghost Fluid 方法

Level Set最大優點的是避免了顯式的追蹤運動界面，大大提高了對復雜界面形狀以及界面復雜運動的處理能力，但是在求解物理方程時，需要對界面附近的點做特殊的處理。為此，引入ghost網格方法，在φ(x，y，t)＞0的地方，認為是流體的真實網格，在φ(x，y，t)＜0的地方，認為是流體的 ghost網格。在ghost網格區域，采用Ghost Fluid插值。例如，標量I通過以下方程常數擴展:

而界面上滿足條件I=Iw。文中將變量從φ＞0的區域擴展到φ＜0的區域，因此在上式中取減號。由于對上述方程處理方式的不同，將會得到不同的結果，在文中，按下列方式定義ghost網格上的變量:

式中:下標G表示ghost網格上的變量，下標E表示歐拉實網格上的變量，下標w表示界面上的值，τ是與n垂直的單位向量。

1.4 邊界條件和初始條件

未點燃推進劑表面采用無滑移邊界條件:u=0，v=0，?p/?n=0，T=Tw;已點燃推進劑表面則采用如下邊界條件:Vτ=0，Vn=Vb，?T/?n=0。其中，Vτ是壁面處的切線速度，Vn是壁面處法向速度，Vb為推進劑燃燒對裂紋的噴射速度，按式Vb=rbρs/ρ計算，其中ρs為推進劑密度。因為在整個過程中，只存在燃氣流同裂紋壁面的傳熱，所以認為在燃氣流同推進劑之間的耦合邊界上溫度和熱流密度均連續，即:TwⅠ=TwⅡ，qwⅠ=qwⅡ。其中，Ⅰ代表固體推進劑區域，Ⅱ代表燃氣流動區域。在具體計算中，首先假定耦合邊界上的溫度分布，對推進劑內部進行導熱計算，得到推進劑內部的溫度分布，然后通過計算得到邊界上的熱流密度，再將熱流密度應用于裂紋內燃氣流動的計算之中，得到耦合邊界上新的溫度分布，以此分布作為推進劑邊界的溫度分布，重復上述計算過程直到邊界上溫度分布收斂到一個確定值。

2 算例與結果分析

文中采用圖1所示的物理模型，并根據文獻［8］，給定裂紋處激波前后壓力比值p2/p1=10.87，波前參數取為大氣參數:p1=1.01 ×105Pa，ρ1=1.225 kg/m3，T1=288.15 K，v1=0。然后根據激波關系式計算出裂紋入口處的流場條件。裂紋的寬度取為Wc=1 mm，深度h=4 mm。超壓 Cp定義為 Cp=(pmaxp2)/(0.5ρ2v22)，其中，pmax為裂紋壁面處最大靜壓，p2為入口靜壓，v2為入口氣流速度，ρ2為入口氣流密度。

圖1 發動機裝藥裂紋模型

圖2為不同時刻裂紋內流場的等壓線圖譜變化。圖2(a)為正激波繞射裂紋左側拐角時的等壓線分布圖，激波S1繞過拐角時產生膨脹波，導致激波在裂紋進口處發生彎曲。主激波繼續向右運動，在裂紋右側壁面發生反射，在圖2(b)中可以清晰地看到激波S3和壁面反射形成的反射激波S2。其中S3向裂紋頂端運動，S2向裂紋左側壁面運動。隨著時間的推移，激波S2分裂為兩部分，即圖2(c)中的S2和 S4。S4也向裂紋頂端運動，并且由于S4的運動速度大于激波S3的運動速度，二者逐漸發生融合形成一道更強的激波S5，同時，S4在左側壁面反射形成S6，如圖2(d)所示。S5繼續向裂紋頂端運動并在裂紋壁面反射后形成激波S7，S7向裂紋入口處運動并與S2相互作用形成復雜的流場結構，造成裂紋內的壓力振蕩。從圖中還可以看出，在激波到達裂紋頂端后，裂紋頭部的燃速明顯加快，說明裂紋頭部壓力較大，這與文獻［9］的實驗事實相符，也說明在裂紋不擴展的情況下，頭部的鈍化是必然的。

圖2 激波繞射動態裂紋等壓線圖

圖3 靜止邊界條件下激波繞射裂紋等壓線圖

圖4 不同邊界條件下超壓計算結果對比圖

圖3是靜止邊界條件下計算得到的裂紋內流場等壓線圖譜變化情況，裂紋深度仍為h=4 mm，圖4是兩種邊界條件下計算所得的裂紋頂端最大壓力處超壓對比圖。從圖中可以看出，兩種邊界條件下壓力開始突升的時間基本相同，但在靜止邊界條件下計算的超壓值明顯的大于移動邊界條件下所得的超壓值。主要原因是在移動邊界條件下，隨著時間的推移，裂紋空腔不斷增大，使得激波強度有所減弱。由于超壓隨著時間的變化呈現強烈的非線性，使得這種差別是無法事先估計的，而準確估計裂紋內的超壓對發動機判廢準則的制定具有十分重要的意義，因此必須引入動邊界進行計算。圖5是無量綱時間t=0.5時兩種邊界條件計算所得流線圖，二者形成的流場結構明顯不同，在靜止邊界條件下，裂紋內出現一系列渦結構，而在移動邊界條件下，只得到兩個比較明顯的大的渦結構。因此，為了得到更為真實的流場結構圖譜，也必須引入動邊界進行計算。同時，這種引入移動邊界進行計算的界面追蹤方法也可以有效耦合裝藥燃燒和發動機內流場的計算，為全發動機三維內流場的數值模擬提供了一種很好的研究思路。

圖5 t=0.5時流場中的流線圖

3 結論

通過對計算結果的分析可知，用Level Set方法配合高精度WENO格式求解激波繞射動態裂紋可得到比較理想的結果。通過對激波繞射裂紋的研究，得到了激波在裂紋內部繞射和反射的流場變化，這在幫助我們認識點火沖擊波與裂紋相互作用規律的同時，也可以為發動機的實際判廢準則提供參考。通過對靜態和動態邊界條件下所形成的流場結構進行比較可知，不同的邊界條件所得到的計算結果相差甚遠，邊界條件對流場結構的影響是無法預計的，要得到更加接近實際的裂紋內流場，必須考慮裂紋由于燃燒而產生的區域擴展。同時由于Level Set方法對復雜界面的追蹤能力，也為進一步考慮裂紋的擴展提供了一種可能的發展方向。

［1］Kuo K K，Morcei J.Crack propagation and branching in burning solid propellants［C］∥Twenty-first Symposiums on Combustion/The Combustion Institute，1986:1933－1941.

［2］Lu YC，Kuo K K.Modeling and numerical simulation of combustion process inside a solid propellant crack［J］.Propellants Explosives Pyrotechnics，1994(19):217－226.

［3］崔小強，白曉征，周偉，等.點火沖擊波與藥柱裂紋流固耦合作用機理數值模擬研究［J］.固體火箭技術，2010，33(1):9－12.

［4］趙汝巖，于勝春，李昊，等.點火升壓階段藥柱裂紋變形研究［J］.固體火箭技術，2009，32(1):43－47.

［5］徐春光，劉君.激波在狹縫中繞射過程研究［J］.固體火箭技術，2003，26(3):39－41.

［6］劉儒勛，王志峰.數值模擬方法和運動界面追蹤［M］.合肥:中國科學技術大學出版社，2001.

［7］Jiang G S，Shu C W.Efficient implementation of weighted ENO schemes［J］.Journal of Computational Physics，1996，126:202－228.

［8］張文普.固體推進劑裂紋內燃燒流動規律研究［D］.西安:西北工業大學，2000.

［9］葛愛學.固體火箭發動機點火過程與裝藥裂紋相互作用機理研究［D］.長沙:國防科學技術大學，2004.