劉世聰，王秋生，婁浩然

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

0 引 言

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

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

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

1 離心機水下爆炸試驗

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

圖1 LXJ-4-450離心機

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

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

表1 試驗工況

2 相似準則在水下爆炸數(shù)值計算中的驗證

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

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

3 數(shù)值模擬

3.1 狀態(tài)方程

狀態(tài)方程是表征物體壓強、密度、比內(nèi)能3個參量的函數(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炸藥采用標準的JWL狀態(tài)方程[12]：

式中，A，B，R1，R2，ω為狀態(tài)方程參數(shù)，A、B單位為GPa，R1，R2，ω為無量綱參數(shù)；Pe為爆轟產(chǎn)物的壓力，MPa；V為相對體積，m3；E0為爆轟產(chǎn)物的初始比內(nèi)能，GJ/m3；試驗所用的炸藥是專業(yè)研制的，因此需要利用一種半經(jīng)驗的方法[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)方程，多項式狀態(tài)方程。Shock狀態(tài)方程是以沖擊Hugoniot點（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ì)確定。Г表達式如下：

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

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

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

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

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

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

當(dāng)水膨脹時（μ＜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為壓強量綱常數(shù)；B0=B1=0.28為無量綱常數(shù)。

3.2 數(shù)值模型

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

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

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

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

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

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

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

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

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

L——爆距，mm；

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

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

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

3.3.2 人工粘性系數(shù)

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

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

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

l——計算單元長度；

ρ'——材料密度；

c——聲速；

υ——體積變化率。

在二維模型Num-2D-40g中改變C1，C2取值來研究人工粘性系數(shù)對沖擊波峰值的影響。取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以此研究一次項人工粘性系數(shù)對計算結(jié)果的影響；對于二次項人工粘性系數(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ù)對氣泡脈動的影響。

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

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

3.3.3 網(wǎng)格尺寸

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

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

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

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

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

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

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

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

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

4 離心機水下爆炸試驗驗證

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

圖10 不同λ下的第二周期氣泡面積與脈動時間曲線

圖11 試驗沖擊波壓力曲線

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

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

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

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

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

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

5 結(jié)束語

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

1）水的多項式狀態(tài)方程計算出的沖擊波峰值較Shock狀態(tài)方程??；水狀態(tài)方程的不同對氣泡最大半徑、氣泡脈動周期影響不大。

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

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

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