水下爆炸沖擊波數值仿真精度研究

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

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

0 引言

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

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

1 水下爆炸數值仿真

1.1 數值模型

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

圖1 數值計算模型Fig.1 Numerical calculation model

1.2 數值模型

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

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

表1 TNT 狀態方程參數Table 1 State equation parameters for TNT

當水介質處于壓縮狀態時,其狀態方程為

當水介質處于膨脹狀態時,其狀態方程為

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

表2 水狀態方程參數Table 2 State equation parameters for water

空氣使用POLYNOMIAL 狀態方程進行描述,其形式為

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

表3 空氣狀態方程參數Table 3 State equation parameters for air

1.3 沖擊波峰值壓力計算結果

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

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

實時監測柴油發電機的運行狀態、報警信息及運行參數，并對柴油發電機的輸出電源的質量進行監測。在動力機房室，通過發電機廠家提供的專用智能通訊接口，及相應的通訊協議，軟件工程師將根據此編集驅動，導入系統平臺進行實時監測；實時地監視柴油機的起停狀態，各相電壓，電流，頻率等參數進行監測。

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

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

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

2 網格尺寸及人工粘性對Pm 的影響

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

式中:ρ為材料密度;Q1為二次項人工粘性系數;Q2為一次項人工粘性系數;l為特征長度;C為當地聲速;ε為體積變化率。

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

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

網格尺寸及Q2對Pm結果如圖2～圖6 所示。由圖2 和圖3 可知,當λ較小時,Pm較低,與經驗公式值偏差較大,隨著λ的增加,Pm顯著增大且近場峰值壓力與經驗值吻合較好。但隨著λ增加,Pm趨于穩定,不會隨著λ的變化出現顯著變化,這表明隨著網格密度的增加,Pm對網格的敏感性越來越低。與此同時,網格尺寸對沖擊波超壓爬升至峰值所需時間影響很大,λ越小,沖擊波超壓爬升至峰值時間越長且峰值壓力越低,如圖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]中有所不同。在實際計算中發現,當網格密度較大時,過小的Q2反而導致Pm數值解與經驗值的偏差增大,如表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

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

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

3 計算誤差預估模型

為了便于應用,實現在預定精度下快速確定網格尺寸和一次項系數,需構建出關于λ和Q2的誤差EP預估模型?？紤]到實際工程需要,誤差應控制在20%以內。同時,由前面分析可知,較大的網格密度以及過小的一次項系數均可能導致計算誤差增大。因此,將誤差預估模型的變量區間范圍限定在0.03≤Q2≤0.10 和3≤λ≤8,則圖7 中的數據在該區間范圍顯示為圖8 中曲線。

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

將圖8 中的EP和λ取對數得到lgEP和lgλ的關系曲線如圖9 所示,k為斜率,b為截距?？梢钥吹?各曲線呈近似平行的線性關系,因此可構造線性表達式為

圖9 lgEP 和lgλ 關系曲線Fig.9 The relationship curves between lgEP and lgλ

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

表6 lgEP 關于lgλ 的擬合參數Table 6 Fitting parameters for λ in EP

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

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

表7 b 關于Q2 的擬合參數Table 7 Fitting parameters for b in Q2

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

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

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

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

4 結論

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

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

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

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