水下爆炸問題的數(shù)值模擬及離心機(jī)試驗(yàn)驗(yàn)證

劉世聰，王秋生，婁浩然

(北京工業(yè)大學(xué)建筑工程學(xué)院，北京 100124)

0 引 言

水下爆炸研究在軍事和民用建設(shè)方面都具有重要作用，其載荷主要包括兩部分：沖擊波以及氣泡脈動(dòng)，兩者大約各占爆炸總能量的50%。水下爆炸研究初期，試驗(yàn)多以大當(dāng)量炸藥爆炸[1-2]和機(jī)理試驗(yàn)[3-4]為主，但在大多數(shù)機(jī)理試驗(yàn)過程中要忽略重力的影響，Takayama等[5]提出在沖擊波階段可以不考慮重力的影響，但氣泡脈動(dòng)階段必須考慮。因此，考慮重力效應(yīng)的縮比模型試驗(yàn)成為水下爆炸試驗(yàn)研究的主要方式，離心機(jī)水下爆炸便是縮比模型試驗(yàn)的一種，其利用離心機(jī)在水下形成超重力場(chǎng)，達(dá)到極少量炸藥模擬原型大當(dāng)量炸藥的爆破效果。中國(guó)水利水電科學(xué)研究院已在離心機(jī)水下爆炸領(lǐng)域做了一些探索[6-7]。

隨著計(jì)算機(jī)技術(shù)的不斷進(jìn)步，數(shù)值計(jì)算已成為研究炸藥水中爆炸的重要手段之一。水的狀態(tài)方程、人工粘性系數(shù)、網(wǎng)格尺寸都是影響數(shù)值計(jì)算精度的重要因素，針對(duì)這幾個(gè)影響因素的分析，國(guó)內(nèi)外已經(jīng)進(jìn)行了一系列的研究。李曉杰等[8]對(duì)不同水狀態(tài)方程進(jìn)行理論分析并提出了不同水狀態(tài)方程的適用范圍，然后通過數(shù)值模擬驗(yàn)證了理論分析的結(jié)果。Huang等[9]建立一維楔體中心對(duì)稱模型，對(duì)不同藥量的球形炸藥水下爆炸進(jìn)行了計(jì)算，研究了人工粘性系數(shù)對(duì)于水下沖擊波的影響。Shi[10]等對(duì)比了空中爆炸模型在不同網(wǎng)格尺寸下的沖擊波峰值、沖擊波比沖量等數(shù)值計(jì)算結(jié)果及誤差，提出了減小網(wǎng)格尺寸帶來的誤差的方法。但目前許多數(shù)值模擬研究缺少試驗(yàn)數(shù)據(jù)的驗(yàn)證，僅僅依靠數(shù)值計(jì)算本身很難具有普遍的適用性。

本文利用RDX（黑索金）炸藥離心機(jī)水下爆炸相對(duì)應(yīng)的數(shù)值模擬研究，首先驗(yàn)證了相似準(zhǔn)則在數(shù)值模擬中的合理性，同時(shí)建立具有不同水狀態(tài)方程、人工粘性系數(shù)、網(wǎng)格尺寸的模型，分析其對(duì)水下爆炸沖擊波及氣泡脈動(dòng)數(shù)值計(jì)算結(jié)果的影響并提出了一些參數(shù)取值的建議，最后將數(shù)值計(jì)算結(jié)果與離心機(jī)水下爆炸試驗(yàn)數(shù)據(jù)對(duì)比，驗(yàn)證了水下爆炸沖擊波及氣泡脈動(dòng)規(guī)律。

1 離心機(jī)水下爆炸試驗(yàn)

離心機(jī)水下爆炸試驗(yàn)是在中國(guó)水利水電科學(xué)研究院IWHR試驗(yàn)室的LXJ-4-450離心機(jī)（見圖1）上進(jìn)行的。隨著離心機(jī)的轉(zhuǎn)動(dòng)，在水下形成一個(gè)比重力加速度g大很多倍的離心加速度N，忽略重力及假設(shè)離心加速度保持不變，小模型便實(shí)現(xiàn)了與原型一樣的重力場(chǎng)環(huán)境，能以極少量炸藥模擬原型大當(dāng)量炸藥水下爆炸效果。

圖1 LXJ-4-450離心機(jī)

如圖2所示，試驗(yàn)中，模型箱的內(nèi)部尺寸為1 280 mm×720 mm×900 mm，模型箱內(nèi)水體深度為600 mm，為提高試驗(yàn)精度，炸藥采用特制的質(zhì)量W為1 g 的RDX（黑索金）球形炸藥，半徑Re為5.25 mm，密度ρ為1.65 g/cm3，并將其至于水下300 mm處。在水下與炸藥相同深度的平面內(nèi)布置水壓力傳感器以獲取沖擊波數(shù)據(jù)，炸藥與傳感器之間的距離即爆距為L(zhǎng)。側(cè)向布置高速攝像機(jī)以捕捉水下爆炸氣泡。試驗(yàn)分別在超重力環(huán)境20，30，40，50g下進(jìn)行。左側(cè)鋼板的尺寸為700 mm×600 mm×50 mm，距離炸藥300 mm，由厚為50 mm的混凝土支座固定，鋼板用于研究結(jié)構(gòu)的動(dòng)態(tài)響應(yīng)，不在本文詳細(xì)描述，本文只是將其看成與模型箱外殼一樣的邊界。各試驗(yàn)工況如表1所示。

圖2 試驗(yàn)?zāi)Ｐ拖洌▎挝唬簃m）

表1 試驗(yàn)工況

2 相似準(zhǔn)則在水下爆炸數(shù)值計(jì)算中的驗(yàn)證

相似準(zhǔn)則[11]是水下爆炸模型試驗(yàn)的理論基礎(chǔ)。文獻(xiàn)[6-7]對(duì)相似準(zhǔn)則內(nèi)容，及相似準(zhǔn)則在離心機(jī)水下爆炸試驗(yàn)中的合理性進(jìn)行了詳細(xì)描述，本文不再贅述。表2列出了離心相似關(guān)系。

利用表2中的相似關(guān)系，將下文中的二維數(shù)值模型Num-2D-40g、三維數(shù)值模型Num-3D-40g換算至原型進(jìn)行計(jì)算，并將原型計(jì)算結(jié)果與兩模型利用相似關(guān)系換算至原型條件下的計(jì)算結(jié)果進(jìn)行對(duì)比，如圖3（a），圖3（b）所示，發(fā)現(xiàn)每幅圖中的兩條曲線都基本重合，因此相似準(zhǔn)則在數(shù)值計(jì)算中也是符合的。這為以離心機(jī)水下爆炸試驗(yàn)為基礎(chǔ)而進(jìn)行的數(shù)值計(jì)算提供了理論支持。

3 數(shù)值模擬

3.1 狀態(tài)方程

狀態(tài)方程是表征物體壓強(qiáng)、密度、比內(nèi)能3個(gè)參量的函數(shù)關(guān)系式。

空氣采用理想氣體狀態(tài)方程[12]：

式中：Pa——空氣壓力，MPa；

γ——絕熱指數(shù)，取1.4；

ρ0——空氣的初始密度，取1.225 kg/m3；

ρ——空氣密度，kg/m3；

E——空氣的比內(nèi)能，GJ/ m3。

表2 離心相似關(guān)系

RDX炸藥采用標(biāo)準(zhǔn)的JWL狀態(tài)方程[12]：

式中，A，B，R1，R2，ω為狀態(tài)方程參數(shù)，A、B單位為GPa，R1，R2，ω為無量綱參數(shù)；Pe為爆轟產(chǎn)物的壓力，MPa；V為相對(duì)體積，m3；E0為爆轟產(chǎn)物的初始比內(nèi)能，GJ/m3；試驗(yàn)所用的炸藥是專業(yè)研制的，因此需要利用一種半經(jīng)驗(yàn)的方法[6]確定JWL狀態(tài)方程的各參數(shù)。最終各參數(shù)取值A(chǔ)=565.295 GPa；B=1.217 11 GPa；R1=4.213 9；R2=1.137 75；ω=0.35；E0=9.536 GJ/m3。

水的狀態(tài)方程有兩種：Shock狀態(tài)方程，多項(xiàng)式狀態(tài)方程。Shock狀態(tài)方程是以沖擊Hugoniot點(diǎn)（PH，μ）表示的Mie-Grüneisen高壓狀態(tài)方程[8]：

式中，μ為壓縮比，其值為ρ/ρ0-1；ρ0為水的初始密度，為1 g/cm3；E為水單位質(zhì)量的內(nèi)能，GJ/m3；PH為沖擊Hugoniot壓力，單位MPa，為含μ的未知量函數(shù)，由沖擊Hugoniot曲線確定；Г為Grüneisen系數(shù)，無量綱參數(shù)，同樣是含μ的未知量函數(shù)，其值由常溫常壓下的Grüneisen系數(shù)Г0和Г的高壓物理性質(zhì)確定。Г表達(dá)式如下：

式中a為常數(shù)，是Г0的一階體積修正常數(shù)。

圖3 原型與模型計(jì)算結(jié)果對(duì)比

在Shock狀態(tài)方程中，沖擊波速度U（km/s）與波后質(zhì)點(diǎn)速度μp（km/s）看作線性關(guān)系，即

式中，k1為沖擊波速度與波后質(zhì)點(diǎn)速度成線性關(guān)系時(shí)的斜率，是經(jīng)驗(yàn)系數(shù)；C0為常溫常壓無擾動(dòng)狀態(tài)聲速，km/s。Shock狀態(tài)方程參數(shù)取值為：Г0=0.28，a=0，k1=1.75，C0=1.483 km/s。

水的多項(xiàng)式狀態(tài)方程[12]根據(jù)壓縮狀態(tài)的不同具有不同的形式。

當(dāng)水壓縮時(shí)（μ＞0），其狀態(tài)方程：

當(dāng)水膨脹時(shí)（μ＜0），其狀態(tài)方程：

式中，E為水單位質(zhì)量的內(nèi)能，GJ/ m3；A1=2.2 GPa，A2=9.54 GPa，A3=1.457 GPa，T1=2.2 GPa，T2=0為壓強(qiáng)量綱常數(shù)；B0=B1=0.28為無量綱常數(shù)。

3.2 數(shù)值模型

基于離心機(jī)水下爆炸試驗(yàn)，建立如圖4、圖5所示的模型進(jìn)行計(jì)算。圖4為二維模型Num-2D，用來研究歷經(jīng)時(shí)間短，精度要求高的沖擊波，其采用軸對(duì)稱方式建立。炸藥質(zhì)量為1 g，半徑為5.25 mm，炸藥、測(cè)點(diǎn)置于水下300 mm，在空氣頂端施加流出邊界，以剛性邊界模擬試驗(yàn)的模型箱，水體、炸藥、空氣均采用Euler單元體建立。圖5為三維模型Num-3D，用來研究歷經(jīng)時(shí)間長(zhǎng)的氣泡脈動(dòng)，在建模時(shí)應(yīng)用網(wǎng)格映射技術(shù)，將二維模型中沖擊波剛要碰觸邊界時(shí)的計(jì)算結(jié)果映射到三維模型中繼續(xù)計(jì)算。三維模型采用面對(duì)稱方式建立，重力環(huán)境、單元體、邊界的選取與二維模型一致。表3列出了數(shù)值模型工況。

圖4 二維模型Num-2D示意圖（單位：mm）

圖5 三維模型Num-3D示意圖（單位：mm）

3.3 影響數(shù)值計(jì)算精度的因素

3.3.1 水的狀態(tài)方程

將二維模型Num-2D-40g、三維模型Num-3D-50g中的水體分別設(shè)置Shock與多項(xiàng)式兩種狀態(tài)方程來研究水的狀態(tài)方程對(duì)沖擊波及氣泡脈動(dòng)的影響。如圖6所示，利用Shock狀態(tài)方程計(jì)算出的沖擊波峰值較多項(xiàng)式狀態(tài)方程要大，且更貼近試驗(yàn)值。由于試驗(yàn)數(shù)據(jù)的限制，圖中的試驗(yàn)曲線是利用Cole[1]提出的沖擊波峰值Pm計(jì)算公式擬合得來，公式如下：

表3 數(shù)值模型工況

圖6 兩種狀態(tài)方程下的沖擊波峰值壓力

式中：W——炸藥的質(zhì)量，kg；

L——爆距，mm；

kp，αp——經(jīng)驗(yàn)系數(shù)，與炸藥種類有關(guān)，無量綱。

由表4對(duì)比發(fā)現(xiàn)，Rm，T幾乎不受水狀態(tài)方程的影響。因?yàn)槎囗?xiàng)式計(jì)算的水體壓力略小于Shock狀態(tài)方程，因此氣泡最大半徑Rm會(huì)比Shock計(jì)算的略大一些，但狀態(tài)方程的不同并不影響氣泡脈動(dòng)周期的大小。

表4 不同水狀態(tài)方程的Rm，T值

3.3.2 人工粘性系數(shù)

人工粘性系數(shù)是為了解決沖擊波波陣面前后壓力、密度等突變問題而引入的。有限元軟件給出的人工粘性項(xiàng)系數(shù)如下式[13]所示：

式中：C1——一次項(xiàng)人工粘性系數(shù)；

C2——二次項(xiàng)人工粘性系數(shù)；

l——計(jì)算單元長(zhǎng)度；

ρ'——材料密度；

c——聲速；

υ——體積變化率。

在二維模型Num-2D-40g中改變C1，C2取值來研究人工粘性系數(shù)對(duì)沖擊波峰值的影響。取C2=0.1，將C1分別設(shè)置為0.002 5、0.005、0.01、0.02、0.05、0.10、0.15、0.20、0.25、0.30以此研究一次項(xiàng)人工粘性系數(shù)對(duì)計(jì)算結(jié)果的影響；對(duì)于二次項(xiàng)人工粘性系數(shù)的研究，取C1=0.02，將C2分別設(shè)置為0.025、0.05、0.1、0.25、0.5、0.75、1.0、1.25、1.5。在三維模型Num-3D-50g中分別設(shè)置C1=1.0，C2=0.02；C1=0.1，C2=0.2；C1=1.0，C2=0.2；C1=0.1，C2=0.02 4種情況研究人工粘性系數(shù)對(duì)氣泡脈動(dòng)的影響。

從表5可以看出，不同的C2對(duì)于沖擊波峰值影響不大。如圖7所示，當(dāng)C1＞0.02，隨著一次項(xiàng)粘性系數(shù)的增大，沖擊波峰值逐漸減小，且各點(diǎn)與試驗(yàn)值之間的誤差都大于10%。當(dāng)C1＜0.02，沖擊波峰值與一次項(xiàng)粘性系數(shù)之間并不成一定的規(guī)律性且計(jì)算值也非常接近，各測(cè)點(diǎn)誤差均小于8%，因此在數(shù)值計(jì)算時(shí)建議C1＜0.02，以提高計(jì)算的精確度。

表6給出了三維模型Num-3D-50g取不同的人工粘性系數(shù)計(jì)算得來的Rm，T值。如表所示，只有當(dāng)一次項(xiàng)人工粘性系數(shù)與二次項(xiàng)人工粘性系數(shù)都取較小的值時(shí)氣泡最大半徑Rm會(huì)略有減小，但是氣泡脈動(dòng)周期T沒有變化。因此，建議在研究氣泡脈動(dòng)Rm及T時(shí)，可忽略人工粘性系數(shù)的影響。

3.3.3 網(wǎng)格尺寸

計(jì)算單元的網(wǎng)格尺寸d與炸藥半徑Re是兩個(gè)重要的長(zhǎng)度參數(shù)。因此，定義無參數(shù)量λ來描述單元網(wǎng)格的疏密程度，λ=Re/d。二維模型Num-2D-40g、三維模型Num-3D-40g的基礎(chǔ)之上，改變網(wǎng)格尺寸以研究其對(duì)沖擊波及氣泡脈動(dòng)計(jì)算結(jié)果的影響。將二維模型的網(wǎng)格尺寸分別取2，1.5，1.25，1，0.8，0.5，0.25，0.2，0.15 mm，對(duì)應(yīng)的λ為2.625、3.5、4.2、5.25、6.562 5、10.5、21、26.25、35；三維網(wǎng)格尺寸分別取20，15，12，10，7.5，5，2.5 mm，對(duì)應(yīng)的λ為0.262 5、0.35、0.437 5、0.525、0.7、1.05、2.1。為使本文對(duì)網(wǎng)格尺寸的研究工作更具適用性，對(duì)數(shù)值計(jì)算結(jié)果中的氣泡半徑Rn與脈動(dòng)時(shí)間Tn采用無量綱化處理，使其與試驗(yàn)值進(jìn)行對(duì)比，即R′=Rn/Rm，T′ =Tn/T，其中R′和T′ 為無量綱氣泡半徑與無量綱脈動(dòng)時(shí)間。

圖8給出了二維模型下λ與沖擊波峰值Pm的關(guān)系曲線。圖中與各測(cè)點(diǎn)Pm-λ曲線接近的4條虛線是該測(cè)點(diǎn)的試驗(yàn)值。可以看出，隨著λ的增大，即網(wǎng)格尺寸的減小，Pm逐漸增大且趨于穩(wěn)定，在λ=35時(shí)測(cè)點(diǎn)1（L=175 mm）已出現(xiàn)Pm減小的趨勢(shì)。測(cè)點(diǎn)1在λ=10.5時(shí)與試驗(yàn)值最接近，測(cè)點(diǎn)2（L=250 mm）在λ=21時(shí)，而測(cè)點(diǎn)3（L=300 mm）、測(cè)點(diǎn)4（L=350 mm）就目前計(jì)算情況來看，則在λ=35時(shí)與試驗(yàn)值最接近。因此，爆距L的不同對(duì)網(wǎng)格尺寸取值要求也不一樣。爆距近的點(diǎn)，在網(wǎng)格相對(duì)稀疏時(shí)便可使計(jì)算誤差達(dá)到允許范圍，爆距越遠(yuǎn)，要求網(wǎng)格尺寸則越小。因此，建議二維網(wǎng)格尺寸d=（1/30~1/10）Re。

圖9給出了三維模型下R′、T′ 與λ的關(guān)系曲線。可以看出，隨著λ的增大，即網(wǎng)格尺寸d逐漸減小，氣泡最大半徑Rm、脈動(dòng)周期T與相應(yīng)試驗(yàn)值逐漸接近，但并非網(wǎng)格尺寸越小越接近。當(dāng)網(wǎng)格尺寸與炸藥半徑大小相近，即λ=1.05時(shí)，數(shù)值計(jì)算的氣泡最大半徑和脈動(dòng)周期與試驗(yàn)值最接近，之后隨著網(wǎng)格尺寸的減小，如圖中λ=2.1，即網(wǎng)格尺寸約取炸藥半徑的1/2時(shí)，數(shù)值計(jì)算的氣泡半徑與脈動(dòng)周期雖都有所增加，但增幅不大，且數(shù)值計(jì)算時(shí)間會(huì)大大增加。

圖7 不同 C1 計(jì)算結(jié)果對(duì)比

表6 50 g下不同C1，C2的氣泡最大半徑與脈動(dòng)周期值

圖10給出了不同λ下的第二脈動(dòng)周期氣泡面積S與時(shí)間的關(guān)系曲線。由圖可以看出，網(wǎng)格尺寸d越小，第二周期的氣泡最大面積Sm也越大，當(dāng)λ＞0.525時(shí)，即網(wǎng)格尺寸大于炸藥半徑的2倍時(shí)，整個(gè)第二脈動(dòng)周期氣泡面積變化不大，且第二、三次脈動(dòng)周期之間的臨界點(diǎn)不明顯；λ＜0.525時(shí)，可以觀察到氣泡大小變化，且可找出不同周期之間的臨界點(diǎn)；在與炸藥半徑相近時(shí)（ λ= 1），計(jì)算結(jié)果與試驗(yàn)結(jié)果誤差最小。

圖8 不同網(wǎng)格尺寸下各測(cè)點(diǎn)的Pm- λ曲線

圖9 不同網(wǎng)格尺寸下的R'/ T'-λ曲線

綜合上述分析，若要利用三維模型研究氣泡脈動(dòng)，建議網(wǎng)格尺寸近似取與炸藥半徑Re相近的尺寸。

4 離心機(jī)水下爆炸試驗(yàn)驗(yàn)證

依據(jù)前文的總結(jié)，利用Num-2D，Num-3D進(jìn)行各個(gè)重力加速度下的數(shù)值模擬計(jì)算。通過以上分析，水狀態(tài)方程選取Shock狀態(tài)方程，令C1=0.02，C2=0.1，二維網(wǎng)格尺寸選取0.5 mm，三維網(wǎng)格尺寸選取5 mm。水壓力傳感器測(cè)得的各試驗(yàn)工況下的沖擊波壓力曲線[6]如圖11所示。數(shù)值模擬各試驗(yàn)工況得來的沖擊波壓力曲線如圖12所示，因?yàn)槿斯ふ承韵禂?shù)的設(shè)置，導(dǎo)致沖擊波后出現(xiàn)偽振蕩[13]。

圖10 不同λ下的第二周期氣泡面積與脈動(dòng)時(shí)間曲線

圖11 試驗(yàn)沖擊波壓力曲線

圖12 數(shù)值沖擊波壓力曲線

從試驗(yàn)及數(shù)值結(jié)果可以看出，隨著爆距L的增大，沖擊波峰值Pm是逐漸減小的。Exp-04到Exp-07是在不同的超重力加速度下進(jìn)行的，但四者的沖擊波波型基本一致，且數(shù)值計(jì)算中四者的數(shù)值模擬結(jié)果相差也很小，這說明在研究沖擊波峰值時(shí)可以忽略重力影響，同時(shí)也證明了文獻(xiàn)[5]的結(jié)論。

高速攝像機(jī)拍攝以及利用三維模型數(shù)值計(jì)算的40g下氣泡第一周期過程圖如圖13 所示，表7為試驗(yàn)所得數(shù)據(jù)[6]及數(shù)值計(jì)算的氣泡最大半徑Rm以及脈動(dòng)周期T的對(duì)比。由于高速攝像機(jī)拍攝范圍的限制，在氣泡膨脹過程中，有些氣泡已經(jīng)超出拍攝范圍，但本文的氣泡外形仍以完整圓處理。

圖13 40 g下第一周期氣泡運(yùn)行形態(tài)

表7 各重力條件下Rm及T試驗(yàn)與數(shù)值結(jié)果對(duì)比

如圖13 所示，40g下試驗(yàn)（Exp-4）氣泡在9.5 ms膨脹至最大，數(shù)值在9.58 ms膨脹至最大，膨脹過程中呈球形，收縮時(shí)逐漸變?yōu)榧忓N形，試驗(yàn)在18.0 ms收縮至最小，數(shù)值在17.98 ms收縮至最小。數(shù)值計(jì)算得到的Rm和T與各自的試驗(yàn)值誤差都較小，均在3% 以內(nèi)?？梢缘贸鼋Y(jié)論：氣泡脈動(dòng)必須考慮重力作用，且隨著重力加速度的增大，氣泡最大半徑Rm以及脈動(dòng)周期T都是在減小的，也驗(yàn)證了文獻(xiàn)[5]的結(jié)論。

5 結(jié)束語

本文基于離心機(jī)水下爆炸試驗(yàn)，建立有限元數(shù)值模型，首先證明了相似準(zhǔn)則在數(shù)值計(jì)算上也是成立的，隨后對(duì)影響數(shù)值計(jì)算精度的因素-水的狀態(tài)方程、人工粘性系數(shù)、網(wǎng)格密度進(jìn)行了研究，并提出取值建議，最后根據(jù)數(shù)值計(jì)算結(jié)果及離心機(jī)試驗(yàn)結(jié)果得出了水下爆炸沖擊波及氣泡脈動(dòng)的規(guī)律。主要結(jié)論有：

1）水的多項(xiàng)式狀態(tài)方程計(jì)算出的沖擊波峰值較Shock狀態(tài)方程?。凰疇顟B(tài)方程的不同對(duì)氣泡最大半徑、氣泡脈動(dòng)周期影響不大。

2）對(duì)沖擊波影響較大的是一次項(xiàng)人工粘性系數(shù)C1，建議取C1取值小于0.02；人工粘性系數(shù)C1，C2對(duì)氣泡脈動(dòng)幾乎無影響。

3）建議二維模型取炸藥半徑的1/30~1/10；在三維模型研究中，在氣泡脈動(dòng)第二周期內(nèi)，當(dāng)網(wǎng)格尺寸大于炸藥半徑的2倍時(shí)，已無法清楚看清氣泡外形動(dòng)態(tài)變化，建議三維模型取與炸藥半徑相近的尺寸。

4）對(duì)于沖擊波研究，可以忽略重力的影響；對(duì)于氣泡脈動(dòng)研究，必須要考慮重力的影響。