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

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

0 引言

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

數(shù)值仿真研究水下爆炸沖擊波載荷問題時,為了處理沖擊波的強間斷面,抑制波陣面前后的數(shù)值振蕩,引入了人工粘性。由于人工粘性的引入會在幾個網(wǎng)格寬度上光滑沖擊波波陣面,因此水下爆炸沖擊波載荷計算結果直接受人工粘性的影響。同時，人工粘性的引入要求網(wǎng)格尺寸不能過大,否則計算過程中網(wǎng)格會忽略部分波陣面信息,致使峰值過低影響計算精度。Huang 等[4]通過典型三硝基甲苯(Trinitrotolue,TNT)炸藥水下爆炸數(shù)值分析,探討了一次、二次粘性系數(shù)對數(shù)值計算結果的影響,并給出了一定比例爆距范圍內(nèi)的建議值。Wang 等[5]研究了不同裝藥質(zhì)量下網(wǎng)格尺寸對沖擊波峰值的影響,并引入與裝藥半徑和單元邊長相關的無因次量表征網(wǎng)格尺寸。胡亮亮等[6]以常規(guī)TNT 水下爆炸為例,對水的狀態(tài)方程、人工粘性系數(shù)和網(wǎng)格尺寸對于數(shù)值仿真結果的影響進行了研究。張社榮等[7]基于有限元軟件AUTODYN建立了不同炸藥量的水下爆炸數(shù)值模型,對比分析了網(wǎng)格尺寸對不同爆距處沖擊波峰值壓力及比沖量的影響。此外還有其他學者[8-13]討論了網(wǎng)格尺寸及粘性系數(shù)對水下爆炸計算結果的影響,但現(xiàn)有研究無法在預定精度下快速確定網(wǎng)格尺寸和人工粘性系數(shù)。因此文中以TNT 水下自由場爆炸數(shù)值計算為例,探究網(wǎng)格尺寸和一次項人工粘性系數(shù)對水下爆炸沖擊波峰值壓力的影響,綜合分析不同工況下網(wǎng)格密度因子和一次項系數(shù)與峰值壓力平均誤差間的關系,構建出普適性較高的水下爆炸數(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)格的排列走向和過渡方式會影響計算結果,因此采用均勻網(wǎng)格劃分方式。在距裝藥中心0.4～2.8 m/kg1/3比例爆距內(nèi)選定16 個觀測點,采用關鍵字*INITIAL_HYDROSTATIC_ALE初始化靜水壓力,設置水域壓力梯度,以模擬真實條件下的水下壓強環(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

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

當水介質(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 沖擊波峰值壓力計算結果

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

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

實時監(jiān)測柴油發(fā)電機的運行狀態(tài)、報警信息及運行參數(shù)，并對柴油發(fā)電機的輸出電源的質(zhì)量進行監(jiān)測。在動力機房室，通過發(fā)電機廠家提供的專用智能通訊接口，及相應的通訊協(xié)議，軟件工程師將根據(jù)此編集驅(qū)動，導入系統(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)驗公式可以較為準確地預估水下爆炸沖擊波峰值壓力。因此,在探究網(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為當?shù)芈曀?ε為體積變化率。

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

為探究水下爆炸數(shù)值仿真中網(wǎng)格尺寸和Q2對沖擊波峰值壓力的影響,對78 g TNT 水下自由場爆炸過程進行數(shù)值仿真計算。當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結果如圖2～圖6 所示。由圖2 和圖3 可知,當λ較小時,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),當網(wǎng)格密度較大時,過小的Q2反而導致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隨λ的變化關系,如圖7 所示。整體來看,隨著λ的增大和Q2的減小,EP逐漸降低。但以Q2=0.10為例,隨著λ的增大,EP先減小后增大,這是因為網(wǎng)格密度較大時,過小的Q2會加大偽振蕩,致使Pm數(shù)值解與經(jīng)驗的相對誤差增大,因此在數(shù)值計算中不能一味地減小網(wǎng)格尺寸和一次項系數(shù)。同時,當λ=1 時,EP很大,即使調(diào)整Q2也未能使EP滿足工程精度。因此在水下爆炸數(shù)值計算中應首先確定λ,同時調(diào)整Q2方能得到精度較高的計算結果。

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

3 計算誤差預估模型

為了便于應用,實現(xiàn)在預定精度下快速確定網(wǎng)格尺寸和一次項系數(shù),需構建出關于λ和Q2的誤差EP預估模型?？紤]到實際工程需要,誤差應控制在20%以內(nèi)。同時,由前面分析可知,較大的網(wǎng)格密度以及過小的一次項系數(shù)均可能導致計算誤差增大。因此,將誤差預估模型的變量區(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λ的關系曲線如圖9 所示,k為斜率,b為截距?？梢钥吹?各曲線呈近似平行的線性關系,因此可構造線性表達式為

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

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

表6 lgEP 關于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ù)進行線性擬合,結果如表7 所示,可以看出數(shù)據(jù)具有較高擬合精度。將式(9)代入式(8),整理得

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

從式(10)可以看到,誤差EP是關于的函數(shù)。利用圖8 中的數(shù)據(jù)繪制EP與的關系曲線如圖11 所示?？梢钥吹?EP與近似呈線性關系。因此,為了進一步提高誤差預估模型的擬合精度,對圖11 中數(shù)據(jù)重新進行線性擬合,其中,擬合模型的復合相關系數(shù)R=0.998,決定系數(shù)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 所示。可以看到,預估誤差與實際計算誤差相近,因此誤差預估模型對于不同裝藥量和裝藥形狀也具有很高的適用性,有助于建立水下爆炸模型對整體計算精度和網(wǎng)格數(shù)量進行綜合分析,為預定精度的仿真模型設計提供依據(jù)。

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

4 結論

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

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

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

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