激波繞射動(dòng)態(tài)裂紋的數(shù)值模擬*

羅 雷，孫振生，張世英

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

0 引言

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

激波繞射裂紋計(jì)算的一個(gè)主要難點(diǎn)是對(duì)邊界的處理。文獻(xiàn)［5］在假設(shè)裂紋壁面靜止的情況下對(duì)強(qiáng)激波繞射裂紋所形成的流場(chǎng)進(jìn)行了詳細(xì)分析，得出了許多有價(jià)值的結(jié)論。但是隨著推進(jìn)劑的燃燒，計(jì)算區(qū)域是不斷變化的。假設(shè)裂紋壁面靜止對(duì)流場(chǎng)結(jié)構(gòu)的影響是無法估計(jì)的，而引入移動(dòng)邊界進(jìn)行計(jì)算可有效地解決這一問題。

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

1 控制方程與數(shù)值計(jì)算方法

1.1 流場(chǎng)控制方程

二維非定?？蓧嚎sNavier-Stokes方程可寫成如下守恒形式:

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

1.2 界面的追蹤

Level Set的控制方程如下:

其中，Level Set函數(shù)φ如下定義:

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

1.3 Ghost Fluid 方法

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

而界面上滿足條件I=Iw。文中將變量從φ＞0的區(qū)域擴(kuò)展到φ＜0的區(qū)域，因此在上式中取減號(hào)。由于對(duì)上述方程處理方式的不同，將會(huì)得到不同的結(jié)果，在文中，按下列方式定義ghost網(wǎng)格上的變量:

式中:下標(biāo)G表示ghost網(wǎng)格上的變量，下標(biāo)E表示歐拉實(shí)網(wǎng)格上的變量，下標(biāo)w表示界面上的值，τ是與n垂直的單位向量。

1.4 邊界條件和初始條件

未點(diǎn)燃推進(jìn)劑表面采用無滑移邊界條件:u=0，v=0，?p/?n=0，T=Tw;已點(diǎn)燃推進(jìn)劑表面則采用如下邊界條件:Vτ=0，Vn=Vb，?T/?n=0。其中，Vτ是壁面處的切線速度，Vn是壁面處法向速度，Vb為推進(jìn)劑燃燒對(duì)裂紋的噴射速度，按式Vb=rbρs/ρ計(jì)算，其中ρs為推進(jìn)劑密度。因?yàn)樵谡麄€(gè)過程中，只存在燃?xì)饬魍鸭y壁面的傳熱，所以認(rèn)為在燃?xì)饬魍七M(jìn)劑之間的耦合邊界上溫度和熱流密度均連續(xù)，即:TwⅠ=TwⅡ，qwⅠ=qwⅡ。其中，Ⅰ代表固體推進(jìn)劑區(qū)域，Ⅱ代表燃?xì)饬鲃?dòng)區(qū)域。在具體計(jì)算中，首先假定耦合邊界上的溫度分布，對(duì)推進(jìn)劑內(nèi)部進(jìn)行導(dǎo)熱計(jì)算，得到推進(jìn)劑內(nèi)部的溫度分布，然后通過計(jì)算得到邊界上的熱流密度，再將熱流密度應(yīng)用于裂紋內(nèi)燃?xì)饬鲃?dòng)的計(jì)算之中，得到耦合邊界上新的溫度分布，以此分布作為推進(jìn)劑邊界的溫度分布，重復(fù)上述計(jì)算過程直到邊界上溫度分布收斂到一個(gè)確定值。

2 算例與結(jié)果分析

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

圖1 發(fā)動(dòng)機(jī)裝藥裂紋模型

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

圖2 激波繞射動(dòng)態(tài)裂紋等壓線圖

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

圖4 不同邊界條件下超壓計(jì)算結(jié)果對(duì)比圖

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

圖5 t=0.5時(shí)流場(chǎng)中的流線圖

3 結(jié)論

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

［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］崔小強(qiáng)，白曉征，周偉，等.點(diǎn)火沖擊波與藥柱裂紋流固耦合作用機(jī)理數(shù)值模擬研究［J］.固體火箭技術(shù)，2010，33(1):9－12.

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

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

［6］劉儒勛，王志峰.數(shù)值模擬方法和運(yùn)動(dòng)界面追蹤［M］.合肥:中國(guó)科學(xué)技術(shù)大學(xué)出版社，2001.

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

［8］張文普.固體推進(jìn)劑裂紋內(nèi)燃燒流動(dòng)規(guī)律研究［D］.西安:西北工業(yè)大學(xué)，2000.

［9］葛愛學(xué).固體火箭發(fā)動(dòng)機(jī)點(diǎn)火過程與裝藥裂紋相互作用機(jī)理研究［D］.長(zhǎng)沙:國(guó)防科學(xué)技術(shù)大學(xué)，2004.