周 超,吉耿杰,魏繼鋒

(1. 92228部隊, 北京 100072;2.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國家重點實驗室, 北京 100081)

0 引言

水中爆炸載荷及其特性研究是水下爆炸技術(shù)領(lǐng)域研究的核心,對水下爆破技術(shù)應(yīng)用具有重要意義。早在1948年,Cole就根據(jù)大量的水下爆炸試驗,分析了水下爆炸過程中的基本物理現(xiàn)象、水下爆炸載荷傳播特性,并提出半理論半經(jīng)驗計算公式[1]。在此基礎(chǔ)上,Zamyshlyayev研究了沖擊波在自由液面及水底的強非線性效應(yīng)[2]。而在氣泡載荷方面,在考慮了流體的可壓縮性后,Prosperetti[3]和Lezzi[4]給出了氣泡動力學(xué)一階理論模型和二階理論模型?；跉馀輨恿W(xué)方程以及試驗數(shù)據(jù),學(xué)者們給出計算氣泡脈動參量的經(jīng)驗公式,例如,Swift等[5]給出了計算TNT和特屈兒爆炸氣泡最大半徑及脈動周期的經(jīng)驗公式;Slifko[6]給出了適用于計算遠場二次壓力波正壓持續(xù)時間和沖量的經(jīng)驗公式;Van Aanhold[7]給出了計算氣泡在第一個脈動周期內(nèi)鉛直位移、氣泡能量損失的經(jīng)驗公式。王振宇等[8]根據(jù)一維不可壓縮流體力學(xué)理論,建立了水下爆炸氣泡脈動規(guī)律和水中壓力分布規(guī)律的基本方程,并引入虛擬力和氣泡能對方程進行了改進。

由于水下爆炸試驗耗資巨大且測試困難,數(shù)值仿真技術(shù)已成為研究其的重要手段。數(shù)值仿真是研究水下爆炸載荷的重要方法之一。常見的數(shù)值方法包括有限差分法[9-11]、有限體積法[12-13]、有限元法[14-15]和SPH法[16]等。Barras等[17]運用ALE法對空間收斂性、水域尺寸和邊界效應(yīng)對爆炸氣泡的影響進行了研究;Hao Huang、劉科種、丁寧、方斌等[18-21]對網(wǎng)格密度、狀態(tài)方程、人工粘性系數(shù)對水下爆炸沖擊波載荷特性的影響進行了研究,但僅僅是定性描述;辛春亮等[22]采用了ALE和Euler算法,發(fā)現(xiàn)網(wǎng)格尺寸及網(wǎng)格質(zhì)量對模擬近場水下爆炸的計算結(jié)果影響很大。Sang-GabLee[23]仿真研究了艦船在水下爆炸載荷作用下的動態(tài)響應(yīng),探討了流體網(wǎng)格尺寸、流體邊界形狀、自由表面尺度等一些重要參數(shù)的選取對計算結(jié)果的影響。

本文中擬運用AUTODYN軟件研究水下爆炸載荷特性,分析數(shù)值仿真參數(shù)設(shè)置對水下爆炸氣泡脈動與沖擊波載荷計算的影響,提出數(shù)值模擬時水域尺寸及網(wǎng)格尺寸的取值方法,并與已有實驗數(shù)據(jù)對比以驗證其準確性,以期為水下爆炸研究提供技術(shù)基礎(chǔ)和分析依據(jù)。

1 水下爆炸載荷理論計算方法

水下爆炸載荷研究主要是指沖擊波和氣泡脈動載荷的研究。

1.1 沖擊波載荷

Cole基于大量水下爆炸實驗結(jié)果擬合出了水下爆炸載荷參量計算公式,Zamyshlyayev在此基礎(chǔ)上做了進一步細化。TNT炸藥水下爆炸沖擊波峰值壓力隨爆距和藥量的變化公式為

式中:r為爆距,m;W為藥量,kg;r0為半徑,m。

Zamyshlyaev提出的衰減常數(shù)為

水下爆炸沖擊波是一種寬頻波,沒有固定的波長。但工程應(yīng)用中,指數(shù)波的衰減常數(shù)是一個重要的時間特征參數(shù),是優(yōu)勢波分量的周期,通常以其作為沖擊波的周期,對應(yīng)的波長就定義為沖擊波的波長。單元網(wǎng)格的粗細與波長密切相關(guān)[24],因此,單元尺寸須慎重選擇。采用一個沖擊波波長內(nèi)的網(wǎng)格數(shù)量λ=cθ/d,作為網(wǎng)格密度的特征參量,研究其對沖擊波壓力峰值的影響。

1.2 氣泡脈動載荷

Cole提出的TNT炸藥氣泡最大脈動半徑與脈動周期計算表達式為

Rmax=3.383W1/3/((h+9.8)1/3)

Tmax=2.064W1/3/((h+9.8)5/6)

式中:Rmax為最大氣泡半徑,m;Tmax為脈動周期,s;h為爆炸水深,m。

2 數(shù)值仿真模型

2.1 模型構(gòu)建方法

采用一維楔形網(wǎng)格構(gòu)建水下爆炸計算模型,單位為mm-mg-ms,炸藥為球形裝藥,如圖2所示。爆心設(shè)于原點即炸藥中心,炸藥和水均采用Euler網(wǎng)格剖分。

2.2 狀態(tài)方程

水采用多項式狀態(tài)方程,通過改變比內(nèi)能來改變靜水壓力,從而改變水域深度。當(dāng)水壓縮時,水的狀態(tài)方程為

p=A1μ+A2μ2+A3μ3+(B0+B1μ)ρ0e,μ>0

當(dāng)水拉伸時,水的狀態(tài)方程為

p=T1μ+T2μ2+B0ρ0,μ<0

式中:μ為壓縮比,μ=ρ/ρ0-1;A1、A2、A3、B0、B1、T1、T2為材料常數(shù);e為水的比內(nèi)能,e=(p0+ρ0gh)/B0ρ0。水的狀態(tài)方程參數(shù)如表1所示。

表1 水的狀態(tài)方程參數(shù)

TNT炸藥采用JWL狀態(tài)方程

式中:V為相對比容;E為單位體積內(nèi)能;A、ω、R1、B、R2為材料參數(shù)。具體參數(shù)如表2所示。

表2 TNT炸藥的狀態(tài)方程參數(shù)

3 數(shù)值仿真參數(shù)設(shè)置

3.1 邊界條件的選擇

采用一維楔形網(wǎng)格建模進行水下爆炸計算時,不同邊界條件對沖擊波特性基本無影響,對氣泡脈動特性影響較大,有些學(xué)者選擇flow-out邊界條件模擬無限水域,也有選用transmit邊界條件[25-27]。

這里分別選用transmit、flow-out和None邊界條件,靜水壓力為0.1 MPa,分析1 g TNT球形炸藥、水域大小為1 m水下爆炸計算時,邊界條件對氣泡脈動的影響,計算結(jié)果如表3所示。

表3 邊界約束的選擇

由表3可以看出,采用flow-out與None邊界條件時得到的氣泡脈動特性完全一致,說明flow-out邊界條件在一維楔形網(wǎng)格計算中不起作用;采用transmit邊界條件計算得到的氣泡脈動特性與理論值也有較大偏差,因此該邊界條件也不夠準確。在實際計算中,應(yīng)設(shè)置足夠大的水域來模擬無限水域環(huán)境,從而減少水域邊界反射對氣泡脈動特性的影響。

3.2 水域尺寸的影響

水域尺寸對沖擊波載荷影響小,對氣泡脈動載荷影響大,在此研究水域尺寸對氣泡脈動半徑的影響十分必要。在1 g藥量TNT下,水域尺寸分別為1、2、5、10 m,將水域尺寸與裝藥半徑之比作為橫坐標(biāo),得到的氣泡脈動特性曲線與理論對比如圖3所示?？梢钥闯?氣泡半徑與氣泡脈動周期隨著水域尺寸的增大呈冪指數(shù)增長趨勢,在600倍裝藥半徑之內(nèi)增長比較迅速,之后趨于平穩(wěn),綜合準確度與計算量,在本裝藥條件下水域尺寸選擇1 200倍裝藥半徑較為適合。

3.3 網(wǎng)格密度對沖擊波特性的影響

改變網(wǎng)格密度,研究1 g球形TNT裝藥在測點0.388 m處的沖擊波特性,波形曲線如圖4所示?？梢钥闯?隨著網(wǎng)格密度的增大,沖擊波壓力上升為迅速,沖擊波壓力峰值更高,且波形下降更加平緩。

1) 不同比例距離沖擊波壓力峰值的影響。

對于1 g球形TNT裝藥,在不同網(wǎng)格密度時,各觀測點得到的沖擊波壓力峰值結(jié)果如圖5所示。隨著特征參量的增大,各測點處的沖擊波壓力峰值均成指數(shù)增大趨勢,在100 內(nèi)迅速增加,在100 后逐漸平穩(wěn)并保持不變。根據(jù)理論經(jīng)驗公式,250、500、750、1 000 mm處的沖擊波壓力峰值分別為18.61、8.50、5.38、3.88 MPa。仿真值與理論值相比誤差小于5%,說明仿真模型計算可靠。

取各網(wǎng)格密度下沖擊波壓力峰值與最小網(wǎng)格尺寸下沖擊波壓力峰值的相對誤差Er,得到結(jié)果如圖6所示?？梢钥闯鲭S著網(wǎng)格密度的增大,不同爆距處的峰值壓力誤差均呈減小趨勢,且爆距越小,誤差越小,網(wǎng)格密度對沖擊波峰值壓力的影響越小,計算結(jié)果越準確。

對上圖平滑后取5%誤差與10%誤差點擬合得到圖7,橫坐標(biāo)為比例距離Z=r/W1/3:

由此得到了5%誤差和10%誤差時最小網(wǎng)格密度隨比例距離變化的取值公式:

5%誤差時的擬合方程:λ=84.35+4.168 97Z

10%誤差時的擬合方程:λ=44.045+2.401 96Z

2) 某比例距離下不同藥量時的沖擊波壓力峰值。

觀測點取相同比例距離,TNT藥量取1、10、100、1 000 g,研究網(wǎng)格密度對沖擊波壓力峰值的影響,結(jié)果如圖8所示?？煽闯霾煌幜肯嗤壤嚯x得到的沖擊波壓力峰值基本一致,且隨著網(wǎng)格密度的增加而增大。

對比壓力峰值的相對誤差,結(jié)果如圖9所示。

取5%誤差與10%誤差點如圖10所示。橫坐標(biāo)為炸藥藥量,縱坐標(biāo)λ。由此圖可以看出,炸藥藥量在相同比例距離時對于網(wǎng)格密度變化造成的沖擊波壓力峰值計算誤差影響很小,因此,綜合之前得到的5%誤差和10%誤差時相同藥量不同比例距離網(wǎng)格尺寸最小取值計算公式,可以得到AUTODYN計算TNT裝藥水下爆炸時考慮計算精度的網(wǎng)格尺寸最小取值計算公式為

圖1 氣泡脈動示意圖

圖2 一維深水爆炸仿真模型

圖3 氣泡脈動特性變化規(guī)律

圖4 不同網(wǎng)格密度下的沖擊波波形

圖5 不同比例距離下的沖擊波壓力峰值

圖6 各比例距離與最小網(wǎng)格尺寸下沖擊波壓力峰值的相對誤差曲線

圖7 5%及10%誤差時網(wǎng)格尺寸與比例距離的關(guān)系

圖8 不同藥量下的沖擊波壓力峰值

圖9 各比例距離與最小網(wǎng)格尺寸下沖擊波壓力峰值的相對誤差曲線

圖10 5%及10%誤差時網(wǎng)格尺寸與炸藥藥量的關(guān)系

5%誤差時:λ=141.7+38.462Z

10%誤差時:λ=83.6+17.486Z

4 試驗驗證

對得到的計算公式與試驗結(jié)果進行對比,鐘帥[28]在加壓水罐中測過1.07 g當(dāng)量炸藥在0.388 m遠處的沖擊波壓力峰值,水深為0.605 m時,測試沖擊波壓力為11.796 MPa。根據(jù)本文中得到的計算公式,設(shè)置的網(wǎng)格尺寸如表4所示。

表4 5%及10%誤差下的網(wǎng)格尺寸

根據(jù)計算得到的網(wǎng)格尺寸建立仿真模型,仿真得到的結(jié)果如表5所示。仿真結(jié)果與實驗測得數(shù)據(jù)相符,驗證了網(wǎng)格尺寸最小取值公式的準確性。

表5 仿真結(jié)果與實驗數(shù)據(jù)的對比

5 結(jié)論

本文中運用AUTODYN數(shù)值模擬軟件,對水下爆炸仿真計算方法進行研究,并與試驗結(jié)果進行對比,主要結(jié)論如下:

1) 采用一維楔形網(wǎng)格進行水下爆炸計算時,邊界條件的選擇對沖擊波特性基本無影響,但對氣泡脈動特性的影響較大;flow-out和transmit邊界條件都無法同時兼顧沖擊波和氣泡2個載荷的計算精度,應(yīng)設(shè)置足夠大的水域才能模擬無限水域環(huán)境。

2) 氣泡脈動周期與氣泡半徑隨水域尺寸的增大呈冪指數(shù)增長趨勢,綜合準確度與計算量,在本裝藥條件下水域尺寸宜選擇1 200倍裝藥半徑。

3) 選取一個沖擊波波長內(nèi)的網(wǎng)格數(shù)量作為網(wǎng)格密度的特征參量,研究得到了不同比例距離和不同藥量時網(wǎng)格密度對計算精度的影響規(guī)律,給出了TNT裝藥水下爆炸在5%誤差和10%誤差下網(wǎng)格尺寸最小取值的計算公式,并經(jīng)實驗驗證了其合理性。