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

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

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

0 引言

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

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

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

1.1 數(shù)值模型

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

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

1.2 數(shù)值模型

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

式中:V為相對(duì)體積;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)時(shí),其狀態(tài)方程為

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

為了便于應(yīng)用,實(shí)現(xiàn)在預(yù)定精度下快速確定網(wǎng)格尺寸和一次項(xiàng)系數(shù),需構(gòu)建出關(guān)于λ和Q2的誤差EP預(yù)估模型?？紤]到實(shí)際工程需要,誤差應(yīng)控制在20%以內(nèi)。同時(shí),由前面分析可知,較大的網(wǎng)格密度以及過小的一次項(xiàng)系數(shù)均可能導(dǎo)致計(jì)算誤差增大。因此,將誤差預(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和λ取對(duì)數(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)對(duì)圖9 中的數(shù)據(jù)進(jìn)行擬合得到擬合系數(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 截距-對(duì)數(shù)粘性系數(shù)線性擬合Fig.10 Linear fitting of intercept-logarithm viscosity coefficient

利用式(9)對(duì)圖10 中數(shù)據(jù)進(jìn)行線性擬合,結(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)系。因此,為了進(jìn)一步提高誤差預(yù)估模型的擬合精度,對(duì)圖11 中數(shù)據(jù)重新進(jìn)行線性擬合,其中,擬合模型的復(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àn)證誤差預(yù)估模型具有較高的普適性,需要對(duì)不同藥量和裝藥形狀的炸藥進(jìn)行仿真計(jì)算,同時(shí)結(jié)合研究背景,對(duì)0.2、5、500、1 500 和5 000 kg TNT 柱形裝藥(長(zhǎng)徑比為1)和球形裝藥水下爆炸進(jìn)行計(jì)算。在距裝藥中心0.327～2.8 m/kg1/3比例爆距間設(shè)置16 個(gè)觀測(cè)點(diǎn),取Q2=0.06,λ=6。不同工況下的峰值壓力平均誤差及預(yù)估誤差如表8 所示?？梢钥吹?預(yù)估誤差與實(shí)際計(jì)算誤差相近,因此誤差預(yù)估模型對(duì)于不同裝藥量和裝藥形狀也具有很高的適用性,有助于建立水下爆炸模型對(duì)整體計(jì)算精度和網(wǎng)格數(shù)量進(jìn)行綜合分析,為預(yù)定精度的仿真模型設(shè)計(jì)提供依據(jù)。

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

4 結(jié)論

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

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

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

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