水下爆炸沖擊波數(shù)值仿真精度研究

敖啟源 ,盧 熹 *,姜智雅 ,康珀閣

(1.沈陽理工大學(xué) 裝備工程學(xué)院,遼寧 沈陽,110159;2.山西江陽化工有限公司,山西 太原,030041)

0 引言

水下武器作為艦船生命力的主要威脅之一,其爆炸沖擊波及氣泡載荷會對艦船造成嚴(yán)重的局部和總體破壞[1]。各國相繼開展了大量水下爆炸試驗,但實彈實測試驗安全風(fēng)險較高、成本高昂、重復(fù)性低且觀測范圍有限[2]。隨著計算機性能和仿真技術(shù)的發(fā)展,數(shù)值仿真以其較高計算精度、低成本和高可重復(fù)性等優(yōu)點成為繼實驗、理論研究后第3 種水下爆炸研究方法[3]。因此開展水下爆炸數(shù)值仿真精度研究具有十分重要的意義。

數(shù)值仿真研究水下爆炸沖擊波載荷問題時,為了處理沖擊波的強間斷面,抑制波陣面前后的數(shù)值振蕩,引入了人工粘性。由于人工粘性的引入會在幾個網(wǎng)格寬度上光滑沖擊波波陣面,因此水下爆炸沖擊波載荷計算結(jié)果直接受人工粘性的影響。同時，人工粘性的引入要求網(wǎng)格尺寸不能過大,否則計算過程中網(wǎng)格會忽略部分波陣面信息,致使峰值過低影響計算精度。Huang 等[4]通過典型三硝基甲苯(Trinitrotolue,TNT)炸藥水下爆炸數(shù)值分析,探討了一次、二次粘性系數(shù)對數(shù)值計算結(jié)果的影響,并給出了一定比例爆距范圍內(nèi)的建議值。Wang 等[5]研究了不同裝藥質(zhì)量下網(wǎng)格尺寸對沖擊波峰值的影響,并引入與裝藥半徑和單元邊長相關(guān)的無因次量表征網(wǎng)格尺寸。胡亮亮等[6]以常規(guī)TNT 水下爆炸為例,對水的狀態(tài)方程、人工粘性系數(shù)和網(wǎng)格尺寸對于數(shù)值仿真結(jié)果的影響進行了研究。張社榮等[7]基于有限元軟件AUTODYN建立了不同炸藥量的水下爆炸數(shù)值模型,對比分析了網(wǎng)格尺寸對不同爆距處沖擊波峰值壓力及比沖量的影響。此外還有其他學(xué)者[8-13]討論了網(wǎng)格尺寸及粘性系數(shù)對水下爆炸計算結(jié)果的影響,但現(xiàn)有研究無法在預(yù)定精度下快速確定網(wǎng)格尺寸和人工粘性系數(shù)。因此文中以TNT 水下自由場爆炸數(shù)值計算為例,探究網(wǎng)格尺寸和一次項人工粘性系數(shù)對水下爆炸沖擊波峰值壓力的影響,綜合分析不同工況下網(wǎng)格密度因子和一次項系數(shù)與峰值壓力平均誤差間的關(guān)系,構(gòu)建出普適性較高的水下爆炸數(shù)值誤差預(yù)估模型,為預(yù)定精度的仿真模型設(shè)計提供依據(jù)。

1 水下爆炸數(shù)值仿真

1.1 數(shù)值模型

基于文獻[14]中開展的爆炸水井試驗,建立二維軸對稱計算模型。試驗水域尺寸4.5 m×9 m,空氣域尺寸4.5 m×0.1 m,炸藥為直徑2 cm 的等高藥柱,炸藥質(zhì)量78 g,放置在水深4 m 位置處,采用中心起爆方式。網(wǎng)格的排列走向和過渡方式會影響計算結(jié)果,因此采用均勻網(wǎng)格劃分方式。在距裝藥中心0.4～2.8 m/kg1/3比例爆距內(nèi)選定16 個觀測點,采用關(guān)鍵字*INITIAL_HYDROSTATIC_ALE初始化靜水壓力,設(shè)置水域壓力梯度,以模擬真實條件下的水下壓強環(huán)境。炸藥、空氣和水介質(zhì)均選用ALE（Arbitrary Lagrangian-Eulerian）算法,數(shù)值計算模型如圖1 所示。

圖1 數(shù)值計算模型Fig.1 Numerical calculation model

1.2 數(shù)值模型

裝藥選用典型的單質(zhì)TNT 炸藥,采用JWL(Jones Wilkins Lee)狀態(tài)方程描述其爆炸過程,具體形式為

式中:V為相對體積;A、B、R1、R2為常數(shù),取值如表1[15]所示;P為壓力;ω為藥量;E為單位體積內(nèi)能。

表1 TNT 狀態(tài)方程參數(shù)Table 1 State equation parameters for TNT

當(dāng)水介質(zhì)處于壓縮狀態(tài)時,其狀態(tài)方程為

當(dāng)水介質(zhì)處于膨脹狀態(tài)時,其狀態(tài)方程為

式中:μ為水的壓縮比;C為水中聲速;S1、S2和S3為常數(shù);γ0為GRUNEISEN 系數(shù);α為體積修正系數(shù);V0為初始相對體積。以上參數(shù)取值如表2所示。

表2 水狀態(tài)方程參數(shù)Table 2 State equation parameters for water

空氣使用POLYNOMIAL 狀態(tài)方程進行描述,其形式為

式中:C0、C1、C2、C3、C4、C5和C6為常數(shù)。取值如表3 所示。

表3 空氣狀態(tài)方程參數(shù)Table 3 State equation parameters for air

1.3 沖擊波峰值壓力計算結(jié)果

根據(jù)文獻[10],長徑比為1∶1 的柱形裝藥可以近似為球形裝藥,對于球形裝藥水下自由場爆炸沖擊波的傳播,Cole[16]通過大量試驗標(biāo)定了水下爆炸沖擊波相似律的公式系數(shù),獲得了不考慮水深影響的TNT 炸藥水下爆炸沖擊波峰值壓力計算公式,Zamyshlyaev[17]在其基礎(chǔ)上將經(jīng)驗公式修正為

式中:Pm為沖擊波峰值壓力;R為爆距;W為裝藥質(zhì)量;Re為裝藥半徑。

實時監(jiān)測柴油發(fā)電機的運行狀態(tài)、報警信息及運行參數(shù)，并對柴油發(fā)電機的輸出電源的質(zhì)量進行監(jiān)測。在動力機房室，通過發(fā)電機廠家提供的專用智能通訊接口，及相應(yīng)的通訊協(xié)議，軟件工程師將根據(jù)此編集驅(qū)動，導(dǎo)入系統(tǒng)平臺進行實時監(jiān)測；實時地監(jiān)視柴油機的起停狀態(tài)，各相電壓，電流，頻率等參數(shù)進行監(jiān)測。

根據(jù)式(5)可計算出水中一定范圍內(nèi)的沖擊波峰值壓力,即表4 中的經(jīng)驗值。為獲得水下爆炸沖擊波傳播演化規(guī)律,文獻[14]在爆炸水井中開展了水下爆炸沖擊波試驗,獲得了78 g TNT 裝藥不同爆距Z處的沖擊波峰值壓力,即表4 中的試驗值。基于爆炸水井試驗,利用LS-DYNA 有限元軟件對78 g TNT 炸藥水下爆炸過程進行仿真,網(wǎng)格尺寸0.25 cm,一次項系數(shù)取值0.06,獲得不同爆距處的沖擊波峰壓力,即數(shù)值解。將試驗值、數(shù)值解及經(jīng)驗值進行對比。

表4 不同爆距處峰值壓力對比Table 4 Comparison of peak pressure at different scaled blast distances MPa

對比表4 數(shù)據(jù),可以發(fā)現(xiàn)數(shù)值解與試驗值最小誤差僅2.6%,平均誤差6.9%;與經(jīng)驗值最小誤差4%,平均誤差9.2%。因此網(wǎng)格尺寸和一次項粘性系數(shù)取值合理時,數(shù)值模型可獲得較高的計算精度。經(jīng)驗值與試驗值峰值壓力平均誤差7.2%,表明經(jīng)驗公式可以較為準(zhǔn)確地預(yù)估水下爆炸沖擊波峰值壓力。因此,在探究網(wǎng)格尺寸及人工粘性對數(shù)值計算精度的影響時,僅比較數(shù)值解與經(jīng)驗值。

2 網(wǎng)格尺寸及人工粘性對Pm 的影響

在有限元計算中為了反映沖擊波波陣面的強間斷,引入了人工粘性來光滑沖擊波,這使得數(shù)值計算峰值壓力低于真實值。一次項人工粘性系數(shù)對沖擊波峰值壓力影響較大,二次項人工粘性系數(shù)主要用于抑制沖擊波衰減過程中的虛假振蕩,但對水中爆炸數(shù)值仿真計算中沖擊波峰值影響較小。參閱文獻[8],二次項人工粘性系數(shù)取定值1.0。LS-DYNA 中人工粘性形式為

式中:ρ為材料密度;Q1為二次項人工粘性系數(shù);Q2為一次項人工粘性系數(shù);l為特征長度;C為當(dāng)?shù)芈曀?ε為體積變化率。

在水下爆炸數(shù)值計算中,網(wǎng)格密度和排列方式對計算結(jié)果影響很大,過大的網(wǎng)格尺寸在計算過程中會忽略沖擊波波陣面信息,沖擊波爬升至峰值所需時間變長;而過小的網(wǎng)格尺寸對計算資源帶來的壓力也不可忽視。同時,水域網(wǎng)格長寬比盡量接近于1,尤其是炸藥附近的網(wǎng)格,否則爆炸沖擊波易出現(xiàn)失真[3]。在探究網(wǎng)格尺寸對數(shù)值計算結(jié)果的影響時,引入與裝藥半徑R0和網(wǎng)格尺寸L0相關(guān)的無因次量,因此文中為研究網(wǎng)格尺寸對數(shù)值計算結(jié)果的影響[5,7-8],其具體形式為

為探究水下爆炸數(shù)值仿真中網(wǎng)格尺寸和Q2對沖擊波峰值壓力的影響,對78 g TNT 水下自由場爆炸過程進行數(shù)值仿真計算。當(dāng)Q2＞0.1 時,計算峰值壓力誤差較大[8],因此調(diào)整分別取Q2=0.02、0.03、0.04、0.05、0.06、0.07、0.08 和0.10。對于二維模型而言,網(wǎng)格密度因子λ=8 時便滿足大部分二維數(shù)值模型,因此取λ分別取λ=1、2、3、4、5、6、7 和8,共計算64 個工況。為了從全局意義上分析網(wǎng)格尺寸和Q2對Pm的影響,引入平均誤差,定義為比例爆距0.375～2.8 m/kg1/3范圍內(nèi)選定觀測點峰值壓力相對于經(jīng)驗公式誤差的平均值。

網(wǎng)格尺寸及Q2對Pm結(jié)果如圖2～圖6 所示。由圖2 和圖3 可知,當(dāng)λ較小時,Pm較低,與經(jīng)驗公式值偏差較大,隨著λ的增加,Pm顯著增大且近場峰值壓力與經(jīng)驗值吻合較好。但隨著λ增加,Pm趨于穩(wěn)定,不會隨著λ的變化出現(xiàn)顯著變化,這表明隨著網(wǎng)格密度的增加,Pm對網(wǎng)格的敏感性越來越低。與此同時,網(wǎng)格尺寸對沖擊波超壓爬升至峰值所需時間影響很大,λ越小,沖擊波超壓爬升至峰值時間越長且峰值壓力越低,如圖4 所示。

圖2 Q2=0.08 時,不同λ 時Pm 隨比例爆距變化曲線Fig.2 Variation of Pm with scaled blast distance for different λ at Q2=0.08

圖3 Q2=0.08 時,不同比例爆距處Pm 隨λ 變化曲線Fig.3 Variation of Pm with λ for different scaled blast distances at Q2=0.08

圖4 λ 對Pm 上升速度影響Fig.4 Effect of λ on the rate of increase of Pm

由圖5 和圖6 所見,隨著Q2逐漸變小,Pm逐漸增大,且近場沖擊波峰值壓力受Q2影響較大。由于Q2取值范圍較小,不同比例爆距處Pm隨Q2的變化曲線與文獻[8]中有所不同。在實際計算中發(fā)現(xiàn),當(dāng)網(wǎng)格密度較大時,過小的Q2反而導(dǎo)致Pm數(shù)值解與經(jīng)驗值的偏差增大,如表5 所示。通過計算,Q2=0.02 時不同爆距處的平均誤差為11.2%,而Q2=0.06 時平均誤差僅4.2%,因此有必要探究λ和Q2對峰值壓力平均誤差的影響。

表5 λ=6 時,不同比例爆距處沖擊Pm 對比Table 5 Different scaled blast distances peak pressure of shock waves at λ=6 MPa

圖5 λ=3 時,不同Q2 時Pm 隨比例爆距變化曲線Fig.5 Variation of Pm with scaled blast distance for different Q2 at λ=3

圖6 λ=3 時,不同比例爆距處Pm 隨Q2 變化曲線Fig.6 Variation of Pm with Q2 for different blast scaled distances at λ=3

對比數(shù)值解和經(jīng)驗值在不同爆距處的沖擊波峰值壓力,獲得不同Q2下峰值壓力平均誤差EP隨λ的變化關(guān)系,如圖7 所示。整體來看,隨著λ的增大和Q2的減小,EP逐漸降低。但以Q2=0.10為例,隨著λ的增大,EP先減小后增大,這是因為網(wǎng)格密度較大時,過小的Q2會加大偽振蕩,致使Pm數(shù)值解與經(jīng)驗的相對誤差增大,因此在數(shù)值計算中不能一味地減小網(wǎng)格尺寸和一次項系數(shù)。同時,當(dāng)λ=1 時,EP很大,即使調(diào)整Q2也未能使EP滿足工程精度。因此在水下爆炸數(shù)值計算中應(yīng)首先確定λ,同時調(diào)整Q2方能得到精度較高的計算結(jié)果。

圖7 峰值壓力平均誤差Fig.7 Average errors of peak pressure

3 計算誤差預(yù)估模型

為了便于應(yīng)用,實現(xiàn)在預(yù)定精度下快速確定網(wǎng)格尺寸和一次項系數(shù),需構(gòu)建出關(guān)于λ和Q2的誤差EP預(yù)估模型。考慮到實際工程需要,誤差應(yīng)控制在20%以內(nèi)。同時,由前面分析可知,較大的網(wǎng)格密度以及過小的一次項系數(shù)均可能導(dǎo)致計算誤差增大。因此,將誤差預(yù)估模型的變量區(qū)間范圍限定在0.03≤Q2≤0.10 和3≤λ≤8,則圖7 中的數(shù)據(jù)在該區(qū)間范圍顯示為圖8 中曲線。

圖8 限定區(qū)間后峰值壓力平均誤差曲線Fig.8 Average errors of peak pressure curves after limited interval

將圖8 中的EP和λ取對數(shù)得到lgEP和lgλ的關(guān)系曲線如圖9 所示,k為斜率,b為截距。可以看到,各曲線呈近似平行的線性關(guān)系,因此可構(gòu)造線性表達(dá)式為

圖9 lgEP 和lgλ 關(guān)系曲線Fig.9 The relationship curves between lgEP and lgλ

利用式(8)對圖9 中的數(shù)據(jù)進行擬合得到擬合系數(shù)如表6 所示。其中,復(fù)合相關(guān)系數(shù)R1均值0.979,決定系數(shù)R2均值0.959,可見數(shù)據(jù)線性擬合精度較高。表6 中各曲線斜率k值比較接近,為了得到歸一化的誤差預(yù)估模型,可以取k值平均值,并且將截距b看作關(guān)于Q2的函數(shù)。如圖10 所示為截距b與lgQ2的關(guān)系,可以看到二者近似呈線性關(guān)系,因此可構(gòu)建線性關(guān)系式

表6 lgEP 關(guān)于lgλ 的擬合參數(shù)Table 6 Fitting parameters for λ in EP

圖10 截距-對數(shù)粘性系數(shù)線性擬合Fig.10 Linear fitting of intercept-logarithm viscosity coefficient

利用式(9)對圖10 中數(shù)據(jù)進行線性擬合,結(jié)果如表7 所示,可以看出數(shù)據(jù)具有較高擬合精度。將式(9)代入式(8),整理得

表7 b 關(guān)于Q2 的擬合參數(shù)Table 7 Fitting parameters for b in Q2

從式(10)可以看到,誤差EP是關(guān)于的函數(shù)。利用圖8 中的數(shù)據(jù)繪制EP與的關(guān)系曲線如圖11 所示。可以看到,EP與近似呈線性關(guān)系。因此,為了進一步提高誤差預(yù)估模型的擬合精度,對圖11 中數(shù)據(jù)重新進行線性擬合,其中,擬合模型的復(fù)合相關(guān)系數(shù)R=0.998,決定系數(shù)R2=0.996,說明擬合精度較高,則式(11)為最終得到的誤差預(yù)估模型

圖11 峰值壓力平均誤差擬合結(jié)果Fig.11 Fitting results of peak pressure average error

為驗證誤差預(yù)估模型具有較高的普適性,需要對不同藥量和裝藥形狀的炸藥進行仿真計算,同時結(jié)合研究背景,對0.2、5、500、1 500 和5 000 kg TNT 柱形裝藥(長徑比為1)和球形裝藥水下爆炸進行計算。在距裝藥中心0.327～2.8 m/kg1/3比例爆距間設(shè)置16 個觀測點,取Q2=0.06,λ=6。不同工況下的峰值壓力平均誤差及預(yù)估誤差如表8 所示。可以看到,預(yù)估誤差與實際計算誤差相近,因此誤差預(yù)估模型對于不同裝藥量和裝藥形狀也具有很高的適用性,有助于建立水下爆炸模型對整體計算精度和網(wǎng)格數(shù)量進行綜合分析,為預(yù)定精度的仿真模型設(shè)計提供依據(jù)。

表8 不同工況下峰值壓力誤差Table 8 Errors of peak pressure for different operating conditions

4 結(jié)論

文中基于LS-DYNA 有限元軟件,分析了網(wǎng)格尺寸和一次項系數(shù)對中近場沖擊波峰值壓力和整體計算誤差的影響,主要得到如下結(jié)論:

1) 在水下爆炸數(shù)值仿真中,近場沖擊波峰值壓力受一次項系數(shù)和網(wǎng)格尺寸影響較大,隨λ 的增大計算峰值壓力對網(wǎng)格的敏感性降低,且網(wǎng)格密度較大時,過小的一次項系數(shù)會加大偽振蕩,致使計算峰值壓力與經(jīng)驗公式值的相對誤差增大;

2) 通過研究網(wǎng)格尺寸及一次項人工粘性對沖擊波峰值壓力的影響,獲得了20%范圍內(nèi)的平均誤差EP與網(wǎng)格密度因子λ和一次項系數(shù)Q2之間的關(guān)系,并進一步擬合,獲得預(yù)定精度下快速確定網(wǎng)格尺寸和一次項人工粘性的預(yù)估模型為

3) 通過0.2～5 000 kg 范圍內(nèi)的TNT 柱形裝藥(長徑比為1)和球形裝藥的水下爆炸計算,驗證了誤差預(yù)估模型可適用于不同裝藥量的柱形(長徑比為1)和球形裝藥二維中近場范圍內(nèi)的水下爆炸計算問題。