固定界面Level Set方法在爆轟波陣面?zhèn)鞑ブ械膽?yīng)用

張莉，吳開騰

(1.內(nèi)江師范學(xué)院 數(shù)學(xué)與信息科學(xué)學(xué)院，四川 內(nèi)江 641110；2.內(nèi)江師范學(xué)院 四川省高等學(xué)校數(shù)值仿真重點(diǎn)實(shí)驗(yàn)室，四川 內(nèi)江 641110；3.內(nèi)江師范學(xué)院 四川省數(shù)據(jù)恢復(fù)重點(diǎn)實(shí)驗(yàn)室，四川 內(nèi)江 641110)

0 引言

爆轟是一類特殊的爆炸現(xiàn)象，主要包括起爆、傳播和近區(qū)效應(yīng)等基本問題。爆轟波傳播是炸藥爆轟過程的中間環(huán)節(jié)，是爆轟近區(qū)效應(yīng)研究的基礎(chǔ)。爆轟波傳播問題一直是爆轟物理研究領(lǐng)域的重要研究內(nèi)容，主要涉及到爆轟波陣面及反應(yīng)區(qū)的運(yùn)動規(guī)律、波后爆轟產(chǎn)物的流場分布、炸藥與惰性介質(zhì)或真空界面處爆轟波與介質(zhì)的相互作用。由于爆轟問題的復(fù)雜性，很難用理論分析方法得到解析解，實(shí)驗(yàn)研究方法是目前最有效的研究手段，但往往存在實(shí)驗(yàn)設(shè)備要求過高、實(shí)驗(yàn)數(shù)據(jù)采集有限、實(shí)驗(yàn)準(zhǔn)備周期較長和存在安全隱患等問題，因此采用數(shù)值模擬方法研究爆轟問題具有重要的意義。隨著計(jì)算機(jī)技術(shù)的迅速發(fā)展，國內(nèi)外爆轟物理領(lǐng)域工作者進(jìn)行了大量爆炸力學(xué)的數(shù)值計(jì)算工作，取得了重要的研究成果。Mader等[1]用二維拉格朗日反應(yīng)流程序模擬了炸藥的拐角繞射現(xiàn)象，認(rèn)為未分解的炸藥總量主要依賴于爆轟波拐角前的曲率。Hallquist等[2]用二維有限元分析程序DYNA2D模擬了爆轟波拐角繞射現(xiàn)象，發(fā)現(xiàn)反應(yīng)延遲現(xiàn)象與拐角處的材料有關(guān)。Bdzil等[3-5]引入沖擊波動力學(xué)思想來解決爆轟波傳播計(jì)算，提出了計(jì)算多維爆轟波傳播的方法，即爆轟波沖擊動力學(xué)(DSD)方法。DSD理論指出爆速曲率關(guān)系可以將爆轟波陣面?zhèn)鞑栴}轉(zhuǎn)化為曲率相關(guān)的界面運(yùn)動問題，而Level Set方法是描述運(yùn)動界面使用最廣泛的方法之一，適用于處理曲率相關(guān)和拓?fù)浣Y(jié)構(gòu)變化的界面運(yùn)動問題。1996年，Aslama等[6]成功地將Level Set方法應(yīng)用到爆轟波沖擊動力學(xué)，發(fā)展了一種適用于任意復(fù)雜二維邊界條件的爆轟波傳播數(shù)值計(jì)算方法。符松等[7]將Level Set方法應(yīng)用到近水面水下爆炸現(xiàn)象的數(shù)值模擬，取得了較好的效果。陳森華等[8]研究了貼體坐標(biāo)系下非理想爆轟波陣面?zhèn)鞑ビ?jì)算的Level Set方法，給出了穩(wěn)定的本質(zhì)無振蕩(ENO)差分格式及邊界處理方法，并對圓弧炸藥進(jìn)行了數(shù)值計(jì)算。文尚剛[9]研究了三維爆轟波傳播的Level Set方法，鐘敏等[10]采用貼體坐標(biāo)系中非正交網(wǎng)格的Level Set方法，研究和發(fā)展了二維和三維爆轟波傳播計(jì)算問題。柏勁松等[11]以DSD理論和Level Set方法為基礎(chǔ)，利用程序燃燒法實(shí)現(xiàn)了波陣面和波后流場的耦合。Wiggins等[12]利用惠更斯原理計(jì)算了二維核爆波在Ia型超新星中的傳播，并利用流體力學(xué)計(jì)算了核爆波在Ia型超新星中的傳播。Clavin[13]綜述了激波前沿三相點(diǎn)的形成及氣體爆轟胞狀結(jié)構(gòu)的分析研究進(jìn)展。Ciccarelli等[14]數(shù)值模擬了爆轟波在帶孔板圓管中傳播的三維結(jié)構(gòu)。Avdeev等[15]對氣泡爆轟波的結(jié)構(gòu)進(jìn)行了數(shù)值模擬。馬天寶等[16]提出了一種弧長算法，計(jì)算了氣相爆轟波在二維管道中的傳播問題，偽弧長算法可以使網(wǎng)格按照一定的形式自適應(yīng)移動，從而達(dá)到在強(qiáng)間斷區(qū)域自動加密的效果，提高了網(wǎng)格的分辨率。

重新初始化是Level Set方法的一個(gè)重要特征，也是Level Set方法數(shù)值實(shí)現(xiàn)的一個(gè)關(guān)鍵問題[17]。重新初始化問題是將隨著演化后不再保持符號距離的Level Set函數(shù)重新變回符號距離函數(shù)，使其Level Set函數(shù)的梯度既不能過大、也不能過小，從而避免數(shù)值振蕩，保證數(shù)值計(jì)算精度?；谄⒎址匠痰碾[式重新初始化方法數(shù)值實(shí)現(xiàn)簡單，界面正則性好，是經(jīng)常采用的一種重新初始化方法。但是，基于偏微分方程的隱式方法很難保證Level Set函數(shù)在重新初始化過程中界面位置不動。這意味著將Level Set方法用于爆轟波陣面計(jì)算時(shí)，容易導(dǎo)致爆轟波陣面位置產(chǎn)生偏離。這種偏離在一次計(jì)算中可能不會造成太大的影響，但是經(jīng)過較長時(shí)間的演化后，誤差的累積會使爆轟波陣面位置產(chǎn)生嚴(yán)重的偏離和失真，計(jì)算也就失去了意義。

在DSD理論的基礎(chǔ)上，為了避免由于誤差累積導(dǎo)致的爆轟波陣面偏離現(xiàn)象，本文提出一種固定界面Level Set方法，其基本思想是在重新初始化相鄰兩次迭代過程中，以固定界面為基本原則，構(gòu)建保持界面不動應(yīng)滿足的等式條件，進(jìn)而推導(dǎo)光滑參數(shù)滿足的計(jì)算公式，主要改進(jìn)在于將傳統(tǒng)Level Set方法中固定的光滑參數(shù)修正為變化的光滑參數(shù)。

1 基于DSD理論的Level Set方程

根據(jù)爆轟沖擊動力學(xué)理論，爆轟波陣面?zhèn)鞑栴}可以轉(zhuǎn)化為計(jì)算界面運(yùn)動速度與界面當(dāng)?shù)厍氏嚓P(guān)的界面運(yùn)動問題，計(jì)算所需要的參數(shù)由實(shí)驗(yàn)給出，且計(jì)算波陣面位置時(shí)不涉及具體狀態(tài)方程形式，具有較強(qiáng)的通用性。在計(jì)算區(qū)域Ω中，假設(shè)爆轟波陣面為一強(qiáng)間斷面，構(gòu)造一個(gè)Level Set函數(shù)φ(x,t)，使得在任意時(shí)刻，爆轟波陣面Γ(t)恰是φ(x,t)的零等值面，即

Γ(t)={x∈Ω|φ(x,t)=0}.

(1)

取φ(x,0)為x點(diǎn)到爆轟波陣面Γ(0)的符號距離函數(shù)，當(dāng)φ=0時(shí)(1)式表示爆轟波陣面坐標(biāo)滿足的方程，當(dāng)φ>0時(shí)(1)式表示炸藥未反應(yīng)區(qū)，當(dāng)φ<0時(shí)(1)式表示爆轟產(chǎn)物區(qū)和化學(xué)反應(yīng)。爆轟波陣面?zhèn)鞑r(shí)，Level Set函數(shù)φ(x,t)應(yīng)當(dāng)滿足

(2)

假設(shè)爆速曲率關(guān)系Dn(κ)近似呈線性關(guān)系，即

Dn(κ0)=DCJ(1-ακ0)，

(3)

式中：DCJ為CJ爆速；α為與炸藥有關(guān)的參數(shù)；κ0為φ=0時(shí)波陣面的中值曲率，即

(4)

式中：φx、φy、φz分別為函數(shù)φ(x,t)對x、y、z的1階偏導(dǎo)數(shù)；φxx、φyy、φzz分別為函數(shù)φ(x,t)對x、y、z的2階偏導(dǎo)數(shù)；φxy、φxz、φyz分別為函數(shù)φ(x,t)對x、y，x、z，以及y、z的混合2階偏導(dǎo)數(shù)。

同樣地，流場中其他部分(φ>0和φ<0)的Dn(κ)可類似處理。若所在點(diǎn)的Level Set曲面中值的曲率為κ，則將κ替換為κ0即可得到所在點(diǎn)的Dn(κ)。

事實(shí)上，(2)式是一個(gè)Hamilton-Jacobi類型方程，只要給定適當(dāng)?shù)某跏紬l件和邊界條件，通過求解(2)式就可以捕捉界面的運(yùn)動和確定界面的位置，從而掌握爆轟波陣面的傳播規(guī)律。

2 固定界面Level Set方法

Level Set函數(shù)φ(x,t)的初值滿足符號距離函數(shù)，可以克服由Level Set函數(shù)梯度過大或過小所造成的數(shù)值振蕩，然而在進(jìn)行了幾個(gè)時(shí)間步求解后，由于數(shù)值方法的內(nèi)在效應(yīng)，Level Set函數(shù)φ(x,t)將不再滿足符號距離函數(shù)，因此需要對Level Set函數(shù)φ(x,t)重新初始化，即求解初值問題(5)式的穩(wěn)態(tài)解來實(shí)現(xiàn)：

(5)

式中：φt為φ(x,t)對t的1階偏導(dǎo)數(shù)；φ0為初始Level Set函數(shù)值；sign(φ0)為符號函數(shù)，即

(6)

ε為光滑參數(shù)。

在Level Set方法重新初始化過程中，最理想的情況是既使Level Set函數(shù)重新變成符號距離函數(shù)的同時(shí)，又保證界面在這過程中不能偏離原來位置。文獻(xiàn)[18]指出，為了減少在界面附近產(chǎn)生的數(shù)值振蕩，光滑參數(shù)一般取ε=Δx(Δx為空間步長)，但這種取法并不能避免界面位置的偏離。

以一維重新初始化問題為例，即(5)式退化為

(7)

(8)

對(8)式進(jìn)行一系列恒等變形后發(fā)現(xiàn)，光滑參數(shù)ε與固定界面位置之間存在數(shù)量關(guān)系[19]，即

(9)

顯然，有

(10)

因此，有

(11)

3 控制方程

爆轟產(chǎn)物流場內(nèi)動力學(xué)方程組為

(12)

p=f(ρ,e).

(13)

4 爆轟波陣面與爆轟產(chǎn)物流場的耦合計(jì)算

從(2)式中容易看出，方程不包含其他流體動力學(xué)參數(shù)，僅僅包含1個(gè)未知量φ，因此一種簡便可行的方法是，將(2)式獨(dú)立在只包含整個(gè)炸藥的區(qū)域內(nèi)求解。波后流場的流體動力學(xué)方程(12)式和(13)式中，如果反應(yīng)進(jìn)程變量λ是關(guān)于時(shí)間和坐標(biāo)的已知函數(shù)，則波后流場的流體動力學(xué)方程組也可以獨(dú)立求解，可以根據(jù)包括惰性介質(zhì)的整個(gè)流場特征來選取坐標(biāo)系及其差分格式。這兩部分獨(dú)立求解結(jié)果的耦合問題，可以通過程序燃燒法[11]來實(shí)現(xiàn)。程序燃燒法是根據(jù)已知的波陣面位置構(gòu)造1個(gè)人工反應(yīng)速率代替真實(shí)反應(yīng)速率。為了避免對時(shí)間的積分，以化學(xué)反應(yīng)分?jǐn)?shù)作為反應(yīng)進(jìn)程變量，其數(shù)值由以下反應(yīng)速率函數(shù)給出：

(14)

式中：d為所在點(diǎn)到波陣面的符號距離，d<0表示點(diǎn)在反應(yīng)區(qū)內(nèi)，d=0表示點(diǎn)在爆轟波陣面上，d>0表示點(diǎn)在未反應(yīng)區(qū)炸藥內(nèi)；dc>0為人工反應(yīng)速率確定的反應(yīng)區(qū)寬度(一般為經(jīng)驗(yàn)值)；對于給定的x、x1、x2，函數(shù)Λ滿足性質(zhì)

(15)

因此，利用(14)式便可計(jì)算反應(yīng)進(jìn)程量λ.

5 數(shù)值算例

空間采用5階加權(quán)無本質(zhì)振蕩格式離散，時(shí)間采用3階全變分下降龍格- 庫塔格式離散，求解(2)式和(12)式。為了敘述方便，本節(jié)將傳統(tǒng)的Level Set方法簡稱為LS方法，將固定界面Level Set方法簡稱為FLS方法。

5.1 剪切速度場

剪切速度場是測試流體介質(zhì)發(fā)生拓?fù)渥冃螘r(shí)運(yùn)動界面捕捉的經(jīng)典算例?？紤]剪切速度分布為

(16)

圖1 t=2.0 s時(shí)刻剪切流場界面位置Fig.1 Interface location of shear flow field for t=2.0 s

計(jì)算區(qū)域?yàn)閇0,1]×[0,1]，初始界面為圓心在(0.5,0.3)處、半徑為0.2的圓周，計(jì)算網(wǎng)格為100×100. 若流場在t=2.0 s時(shí)刻后反轉(zhuǎn)剪切速度，則理論上在t=4.0 s時(shí)刻流體界面應(yīng)返回到初始時(shí)刻界面狀態(tài)。從圖1中可以看出，F(xiàn)LS方法能很好地捕捉到狹窄、細(xì)長的流帶，而LS方法無法捕捉到狹長、細(xì)長的流帶，在流體頭部和尾部均出現(xiàn)嚴(yán)重的體積損失。從圖2中可以看出，采用LS方法得到的返回界面已經(jīng)嚴(yán)重失真，出現(xiàn)了較嚴(yán)重的體積損失，而采用FLS方法得到的返回界面非常接近初始圓周。這是因?yàn)榧羟兴俣葓鲈谡麄€(gè)流場分布固定，LS方法和FLS方法都沒有反作用于速度場，在界面變形過程中，F(xiàn)LS方法在重新初始化過程中能夠很好地固定界面位置，使得在其剪切流場中邊界能夠很好地返回初始狀態(tài)。

圖2 t=4.0 s時(shí)刻剪切流場界面位置Fig.2 Interface location of shear flow field for t=4.0 s

若記剪切運(yùn)動前界面所圍的面積為S0、剪切運(yùn)動后界面所圍的面積為S，則定義運(yùn)動前、后界面面積的比值為err，即err=S/S0，表1給出LS方法和FLS方法計(jì)算的剪切前、后界面所圍面積的對比結(jié)果。從表1中可以看出，隨著網(wǎng)格數(shù)的增加，傳統(tǒng)LS方法和FLS方法的質(zhì)量守恒性都會提高，但在同等網(wǎng)格條件時(shí)，F(xiàn)LS方法質(zhì)量守恒性的優(yōu)勢更明顯。

表1 剪切流場界面面積比值對比結(jié)果

5.2 中心一點(diǎn)起爆

假定爆轟波在炸藥中是理想傳播的，即以不變的理想爆速DCJ在炸藥中進(jìn)行穩(wěn)定傳播。正方形炸藥的邊長為16 mm，起爆點(diǎn)在炸藥的中心，坐標(biāo)為(8 mm,8 mm,8 mm)，設(shè)炸藥的理想爆速為8.0 km/s.

圖3～圖5分別給出了LS方法和FLS方法在不同時(shí)刻爆轟波的傳播情況。從圖3～圖5中可以看出，爆轟波以球形向四周發(fā)散，與文獻(xiàn)[9]計(jì)算結(jié)果是一致的，表明LS方法和FLS方法計(jì)算的正確性。在此基礎(chǔ)上，圖6分別給出了兩種方法在不同時(shí)刻的爆轟波陣面，其中截面為y=8 mm的Oxz平面。從圖6中也可以看出爆轟波陣面的發(fā)展演化過程。根據(jù)中心一點(diǎn)起爆的爆速曲率關(guān)系，如果忽略爆速曲率的高階項(xiàng)，則爆轟波在沒有達(dá)到炸藥邊界的情況下是能夠得到解析式[20](理想爆轟的CJ值)的，這可以作為檢驗(yàn)LS方法和FLS方法計(jì)算精度的參考依據(jù)。為了簡化，考慮一維情形，設(shè)初始時(shí)刻的爆轟波陣面半徑為r0，根據(jù)爆轟曲率關(guān)系，有

(17)

式中：r為爆轟波陣面半徑；r(0)為初始時(shí)刻爆轟波陣面半徑；β為與反應(yīng)區(qū)長度和外殼材料性質(zhì)相關(guān)的系數(shù)。

圖3 t=0.35 μs時(shí)爆轟波的計(jì)算結(jié)果Fig.3 Calculated results of detonation wave for t=0.35 μs

圖5 t=1.05 μs時(shí)爆轟波的計(jì)算結(jié)果Fig.5 Calculated results of detonation wave for t=1.05 μs

圖6 在不同時(shí)刻爆轟波陣面?zhèn)鞑サ挠?jì)算結(jié)果Fig.6 Calculated results of detonation wavefront propagation at different times

求解(17)式可得到爆轟波陣面半徑r(t)的解析式：

(18)

圖7 中心一點(diǎn)起爆爆轟波陣面半徑的對比結(jié)果Fig.7 Comparison results of radii of detonation wavefront

圖7給出了兩種方法計(jì)算的爆轟波陣面半徑與CJ理論值的對比情況。從圖7中可以看出，LS方法與CJ理論值的差距較大，F(xiàn)LS方法的計(jì)算結(jié)果比較接近于CJ理論值。隨著時(shí)間的演化，LS方法的計(jì)算結(jié)果與CJ理論值的偏離程度越來越大，而FLS方法的計(jì)算結(jié)果始終保持接近于CJ理論值。

如表2所示，將一個(gè)計(jì)算時(shí)間(t=1.4 μs)分成4個(gè)相同的計(jì)算周期(Δt=0.35 μs)，觀察在每一段計(jì)算周期內(nèi)LS方法和FLS方法計(jì)算的爆轟波陣面半徑與CJ理論值的絕對誤差之和變化情況。從表2中可以看出，隨著時(shí)間的演化，LS方法的絕對誤差之和越來越大，而FLS方法的絕對誤差之和變化不大，與圖7的結(jié)論是一致的。這是因?yàn)長S方法不重視重新初始化過程中的界面偏離，誤差的累積導(dǎo)致爆轟波陣面產(chǎn)生了偏離，而FLS方法在重新初始化過程中，通過固定界面有效地抑制了爆轟波陣面偏離。

表2 爆轟波陣面半徑絕對誤差之和對比結(jié)果

5.3 圓弧炸藥平面起爆

設(shè)圓弧片材料為JB9014炸藥，其內(nèi)徑、外徑分別為70 mm和100 mm，圓心角為直角，主裝藥內(nèi)、外界面均為5 mm厚的鋼外殼封閉。爆轟曲率關(guān)系見表3，空間步長為0.5 mm，時(shí)間步長為0.01 μs，dc=4 mm，流體動力學(xué)計(jì)算取穿透邊界條件，其他部分取剛性壁邊界條件。

表3 炸藥的爆速曲率關(guān)系

圖8 爆轟波陣面的對比結(jié)果Fig.8 Comparison results of detonation wavefront

圖8給出了LS方法、FLS方法計(jì)算的爆轟波陣面與實(shí)驗(yàn)結(jié)果[11]的對比情況。從圖8中可以看出，由于約束條件處理的不同，LS方法和FLS方法的計(jì)算結(jié)果比實(shí)驗(yàn)結(jié)果[11]總體上要小些，LS方法計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的平均相對誤差大約為9.625 7%，F(xiàn)LS方法計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果的平均相對誤差大約為2.792 3%. 因此，F(xiàn)LS方法比LS方法的計(jì)算結(jié)果與實(shí)驗(yàn)結(jié)果更加吻合。

6 結(jié)論

本文針對誤差累積導(dǎo)致的爆轟波陣面偏離問題，根據(jù)DSD理論，在傳統(tǒng)Level Set方法的基礎(chǔ)上，提出了一種固定界面Level Set方法，其基本思想是在重新初始化相鄰兩次迭代過程中，以固定界面為基本原則，構(gòu)建保持界面不動應(yīng)滿足的等式條件，進(jìn)而推導(dǎo)出光滑參數(shù)滿足的計(jì)算公式，主要改進(jìn)在于將傳統(tǒng)Level Set方法中固定的光滑參數(shù)修正為變化的光滑參數(shù)。

為了驗(yàn)證固定界面Level Set方法在描述爆轟波陣面?zhèn)鞑ビ?jì)算的有效性和適用性，主要討論了剪切速度場、中心一點(diǎn)起爆和圓弧炸藥平面起爆的驗(yàn)證模型。剪切速度場數(shù)值模擬結(jié)果表明，固定界面Level Set方法在重新初始化過程中能夠很好地固定流場界面位置，能夠較為準(zhǔn)確地捕捉流體介質(zhì)發(fā)生拓?fù)渥冃螘r(shí)的運(yùn)動界面，使其流場邊界能夠很好地返回初始狀態(tài)。中心一點(diǎn)起爆和圓弧炸藥平面起爆數(shù)值模擬結(jié)果表明，與傳統(tǒng)Level Set方法相比，固定界面Level Set方法能夠有效地避免由于誤差累積導(dǎo)致的爆轟波陣面偏離現(xiàn)象，使其計(jì)算結(jié)果更接近于CJ理論值和實(shí)驗(yàn)結(jié)果，能夠更準(zhǔn)確地捕捉到合理的爆轟波陣面。