鄭永輝，魏繼鋒

（北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點實驗室，北京 100081）

隨著數(shù)值計算方法的進(jìn)步和計算能力的提升，數(shù)值仿真技術(shù)為解決水下爆炸機(jī)理分析、結(jié)構(gòu)響應(yīng)機(jī)制規(guī)律揭示等科學(xué)技術(shù)難題提供了有力的手段。在水下爆炸數(shù)值仿真研究中，研究人員已對水介質(zhì)狀態(tài)方程、人工黏性系數(shù)、網(wǎng)格密度等影響計算精度的因素開展了較為細(xì)致的探討。

在進(jìn)行水下爆炸問題計算時，需設(shè)定水介質(zhì)的初始參數(shù)以設(shè)置流場初始壓力。已有研究中，大多采用以下兩種方式：一是認(rèn)為水介質(zhì)密度變化很小，可忽略不計，僅通過改變內(nèi)能設(shè)置流場初始壓力；二是認(rèn)為流體靜壓力的增大導(dǎo)致水的密度增大，僅通過調(diào)整比容來設(shè)置初始壓力。當(dāng)采用前一種方式時，沖擊波超壓峰值隨水深的增大略有減小，后一種方式下的變化趨勢則相反?？梢钥闯觯跏紖?shù)的設(shè)置方式會帶來水下爆炸載荷計算結(jié)果的差異，而對此的細(xì)致分析尚未見報道。

除直接設(shè)置水介質(zhì)初始參數(shù)外，也可通過施加重力、采用程序提供的特定功能等方式實現(xiàn)流場壓力初始化。這些方式簡化了操作，但程序仍需提供水介質(zhì)的初始參數(shù)使其具有一定的壓力。以LSDYNA 提供的關(guān)鍵字INITIAL_EOS_ALE 為例，其說明文檔指出，程序?qū)⒉捎靡环N迭代方法根據(jù)用戶提供的壓力計算介質(zhì)的初始內(nèi)能和比容。簡而言之，該類方法其本質(zhì)上也是對水介質(zhì)的初始參數(shù)進(jìn)行設(shè)置。

水介質(zhì)初始參數(shù)設(shè)置本質(zhì)上是對水介質(zhì)參數(shù)隨水深變化的規(guī)律的反映。本文中將首先從常用的狀態(tài)方程中選定符合參考狀態(tài)時參數(shù)的水介質(zhì)狀態(tài)方程；然后從熱力學(xué)角度分析現(xiàn)有兩種設(shè)置方式對應(yīng)的熱力學(xué)過程，并給出第3 種參數(shù)設(shè)置方式，同時以關(guān)鍵字INITIAL_EOS_ALE 為例，對程序給出的初始化結(jié)果進(jìn)行分析；隨后采用LS-DYNA 程序進(jìn)行一維球形裝藥水下爆炸仿真研究，細(xì)致分析前三種初始參數(shù)設(shè)置方式對水下爆炸載荷特性的影響，并與已有研究成果進(jìn)行對比；最后確定適當(dāng)?shù)某跏紖?shù)設(shè)置方式。

1 水介質(zhì)狀態(tài)方程

數(shù)值仿真軟件中常見的水介質(zhì)狀態(tài)方程為Mie-Grüneisen 狀態(tài)方程和Polynomial 狀態(tài)方程。Mie-Grüneisen 狀態(tài)方程以沖擊絕熱線作為參考線，最終形式為：

式中：、、、、、和均為狀態(tài)方程系數(shù)，e為質(zhì)量內(nèi)能增量，、、、、、和均為狀態(tài)方程系數(shù)，、、、、、和均為狀態(tài)方程系數(shù)。

與Mie-Grüneisen 狀態(tài)方程類似，式（2）～（3）也存在與實際壓力不符的問題。式（4）在=時，可滿足= 0，e= 0 以及=。因此，宜選擇Polynomial 狀態(tài)方程式（4）作為水介質(zhì)狀態(tài)方程，具體參數(shù)見表1。

表1 水介質(zhì)狀態(tài)方程參數(shù)Table 1 EOS parameters of water

2 初始參數(shù)設(shè)置

2.1 設(shè)置方式及狀態(tài)參數(shù)變化

根據(jù)水介質(zhì)狀態(tài)方程，設(shè)置初始壓力時可改變參量和e?，F(xiàn)今大多數(shù)仿真研究是通過改變e（方式Ⅰ）來進(jìn)行流場初始壓力設(shè)置。從熱學(xué)角度看，方式Ⅰ中不同水深下的水介質(zhì)參數(shù)按等容過程（d=0）變化，此時外界做功為零，流場的內(nèi)能增量全部由外界傳導(dǎo)的熱量提供（d= d）。這表明采用方式Ⅰ時，流體靜壓力的增大單純源于外界傳熱，與重力等外力做功等因素?zé)o關(guān)，這與實際深水環(huán)境以及加壓模擬深水環(huán)境嚴(yán)重不符。

少數(shù)研究人員通過調(diào)整（方式Ⅱ）來進(jìn)行初始壓力設(shè)置。此時，內(nèi)能變化為零（d= 0），即等內(nèi)能形式，外界對流體做功與外界向流體傳導(dǎo)的熱量之和為零（d= d-d= 0），水介質(zhì)被壓縮后向外放出熱量。為分析溫度變化趨勢，引入含溫度的內(nèi)能表達(dá)式：

式中：= 273 K。當(dāng)0 ≤≤ 0.025 時，隨的增大而逐漸增大（見圖1）。在實際深水環(huán)境中，隨著水深的增大，水溫逐漸減?。辉谕ㄟ^加壓構(gòu)造的模擬深水環(huán)境中，水溫與環(huán)境溫度基本相同，即溫度可認(rèn)為保持不變。因此，從溫度變化趨勢的角度來看，方式Ⅱ與實際情況有一定的差距。

熱力學(xué)研究結(jié)果表明，溫度到處相同是重力場的熱平衡條件。據(jù)此，本文中提出第3 種初始參數(shù)設(shè)置方式，即假設(shè)水介質(zhì)參數(shù)隨水深按照等溫過程（d= 0，方式Ⅲ）變化。對于該方式，Δ= 0，將式（5）代入式（4），可得等溫過程對應(yīng)的狀態(tài)方程。顯然，對于實際深水環(huán)境，方式Ⅲ只是一種粗糙的近似，不過理論上優(yōu)于方式Ⅰ和方式Ⅱ；對于模擬深水爆炸試驗，近似程度較高。

除上述設(shè)置方式外，數(shù)值仿真軟件往往也具備根據(jù)重力和水深信息或者設(shè)定的流體靜壓力計算水介質(zhì)初始參數(shù)的功能，只是用戶往往無法直接獲知其遵循的規(guī)律?，F(xiàn)以關(guān)鍵字INITIAL_EOS_ALE 為例，通過讀取不同水深下的流場壓力初始化結(jié)果來初步了解LS-DYNA 程序中和e隨水深的變化規(guī)律（方式Ⅳ），結(jié)果如圖1 所示，可以看出，方式Ⅳ對應(yīng)的熱力學(xué)過程接近方式Ⅱ。

圖1 中給出了4 種設(shè)置方式在不同水深條件時的壓縮比、內(nèi)能增量e和溫度變化Δ的計算結(jié)果。方式Ⅱ和方式Ⅳ的計算結(jié)果幾乎重合，二者均與方式Ⅲ較接近；方式Ⅰ的計算結(jié)果與其他3 種相去甚遠(yuǎn)。采用方式Ⅱ、Ⅲ、Ⅳ時，隨水深增加而增大；采用方式Ⅰ時則始終為零；在同一水深處，的值從小到大依次為：方式Ⅰ、方式Ⅳ、方式Ⅱ、方式Ⅲ。采用方式Ⅰ和方式Ⅳ時，e隨水深增加而增大，采用方式Ⅱ時e始終為零，采用方式Ⅲ時e隨水深的增加而減??；在同一水深處，e的值從小到大依次為：方式Ⅲ、方式Ⅱ、方式Ⅳ、方式Ⅰ。采用方式Ⅰ、方式Ⅱ和方式Ⅳ時，Δ隨水深增加而增大；采用方式Ⅲ時Δ始終為零；在同一水深處，Δ的值從小到大依次為：方式Ⅲ、方式Ⅱ、方式Ⅳ、方式Ⅰ。

圖1 不同設(shè)置方式下的μ、eV 和ΔT 值Fig. 1 Values of μ, eV and ΔT in different setting modes

2.2 流體可壓縮性

流體的可壓縮性決定流體內(nèi)微弱擾動波的傳播速度，即流體內(nèi)的聲音傳播速度，因此常用聲速來表示流體的可壓縮性。聲速的表達(dá)式為：

圖2 中給出了4 種設(shè)置方式在不同水深條件時的聲速。從圖2 中可以看出，隨著水深的增大，4 種方式下的水介質(zhì)聲速逐漸增大，即可壓縮性均出現(xiàn)減弱，但減弱幅度不同，由小到大的排序為：方式Ⅰ、方式Ⅳ、方式Ⅱ、方式Ⅲ?？梢?，聲速的變化與壓縮比的變化情形一致。

圖2 不同設(shè)置方式下的聲速變化Fig. 2 Acoustic velocity in different modes

3 初始參數(shù)設(shè)置對水下爆炸載荷的影響

3.1 數(shù)值仿真模型

建立了一維球?qū)ΨQ水下爆炸數(shù)值仿真模型，并應(yīng)用LS-DYNA 數(shù)值仿真軟件進(jìn)行計算。其中，炸藥采用MAT_HIGH_EXPLOSIVE_BURN 材料模型以及JWL 狀態(tài)方程，參數(shù)取自文獻(xiàn)[17]。水介質(zhì)的狀態(tài)方程采用式（4），參數(shù)在表1 中列出。

炸藥選用1 kg TNT，水深范圍0 ～ 5 km，由于方式Ⅱ和方式Ⅳ對應(yīng)的初始參數(shù)極為接近，因此只對比前3 種方式，各自對應(yīng)的和e值取自圖1。為保證計算精度，網(wǎng)格尺寸為裝藥半徑的1/100；為減小邊界反射沖擊波對氣泡脈動的影響，水域半徑取237.2 m，大于水深為0 m 時聲波在一個脈動周期內(nèi)傳播距離的1/2。

3.2 沖擊波載荷特性

選取水深為5 km，爆距= 55（為裝藥半徑），計算得到了3 種方式下的沖擊波壓力-時間曲線，如圖3 所示，對應(yīng)的沖擊波到達(dá)時間、超壓峰值Δ和正壓持續(xù)時間列于表2 中?？梢钥闯?，3 條曲線并不完全重合，除方式Ⅰ外，其他2 種方式的仿真結(jié)果相差較小；方式Ⅲ、Ⅱ、Ⅰ的及值依次減小，而Δ依次增大。該規(guī)律與流體可壓縮性結(jié)果基本一致，即流體可壓縮性越弱，沖擊波到達(dá)時間越短，沖擊波超壓峰值越大，正壓持續(xù)時間越短。

表2 H =5 km 及 R =55R0 時3 種方式下的 ta 、 Δ Pm 及tcTable 2 Values of ta , Δ Pm and tc in three modes when H =5 km andR=55R0

圖3 H=5 km 及R=55R0 時3 種方式下的沖擊波壓力-時間曲線Fig. 3 Shock wave pressure-time curves in three modes when H=5 km and R=55R0

圖4 中給出了= 55時3 種設(shè)置方式在不同水深下獲得的沖擊波超壓峰值、沖量和能流密度，其中，和分別由下式給出；

圖4 R=55R0 時不同水深下3 種方式所對應(yīng)的沖擊波載荷參數(shù)Fig. 4 Values of shock wave load in three modes when R=55R0

式中：為沖擊波能流密度，ρ、分別為未擾動流體的密度和聲速，為便于比較，統(tǒng)一取ρ= ρ=1 000 kg/m，== 1 483.2 m/s。

表3 中給出了沖擊波載荷參數(shù)在5 km 水深處相對于0 m 的變化幅度?？梢钥闯?，隨著水深的增大，采用方式Ⅰ時Δ逐漸減小，其他2 種方式則逐漸增大，且變化幅度略大于方式Ⅰ。對于和，3 種方式給出的仿真結(jié)果均隨水深的增大而減小，方式Ⅱ和Ⅲ的變化幅度略小于方式Ⅰ。

表3 沖擊波載荷在5 km 水深處相對于0 m 的變化幅度Table 3 Changing amplitudes of shock wave load when depth changes from 0 m to 5 km

3.3 氣泡脈動特性

分析=5 km 時3 種設(shè)置方式下氣泡的脈動曲線，如圖5 所示，相應(yīng)的氣泡最大半徑和第1 次脈動周期列于表4?？梢钥闯觯绞舰裢?，其余兩條曲線基本重合；方式Ⅲ、方式Ⅱ、方式Ⅰ對應(yīng)的和依次增大。顯然，該規(guī)律與水介質(zhì)壓縮比的變化一致。由于氣泡脈動與水介質(zhì)的慣性密切相關(guān)，當(dāng)體積不變時，水介質(zhì)密度的增大導(dǎo)致其慣性增大，在密度變化較小時慣性增幅有限，因此，3 種方式下氣泡脈動曲線的差別并不明顯。

圖5 H=5 km 時3 種方式下氣泡的半徑-時間曲線Fig. 5 Bubble radius-time curves in three modes when H=5 km

表4 H =5 km 時3 種方式下對應(yīng)的 R max 和TTable 4 Values of R max and T in three modes whenH=5 km

4 對比分析

4.1 沖擊波載荷對比分析

自20 世紀(jì)60 年代起，研究人員開展了較多的水深對沖擊波載荷特性影響的研究。Baum 等在模擬深水爆炸容器中進(jìn)行了深水爆炸試驗，當(dāng)1 ＜/＜ 7 時，Δ在0 ～ 4 km 的模擬水深范圍內(nèi)基本保持不變。Vanzant 等給出了14.75 g Pentolite 炸藥在更大爆距范圍內(nèi)的試驗結(jié)果，本文對其試驗數(shù)據(jù)進(jìn)行了線性擬合，結(jié)果表明，Δ隨水深的增大而略有增大，且在較大爆距處的變化更明顯，如圖6 所示。鐘帥等開展了0～150 m 模擬水深水下爆炸試驗，指出Δ也隨水深增大而略有增加。Slifko、Xiao 等也進(jìn)行了類似研究，不過試驗在大洋中進(jìn)行，爆源和測點并不位于同一深度，而且爆距較大，與本文分析工況不符。

圖6 Vanzant 等[19]的試驗結(jié)果及線性擬合曲線Fig. 6 Test results from Vanzant et al[19] and linear fitting results

理論上，Baum 等從相似性和量綱理論出發(fā)，在淺水經(jīng)驗公式的基礎(chǔ)上，給出了考慮流體靜壓力時Δ的計算式：

式中：和α 為常數(shù)，當(dāng)6 ≤/＜ 12 時，=44.1 MPa，α = 1.5；當(dāng)12 ≤/＜ 240 時，=52.4 MPa，α = 1.13；= 304.9 MPa 為Tait 狀態(tài)方程中的常數(shù)；為幾何指數(shù)，對于球面沖擊波，= 3。

式（12）表明，Δ隨水深的增大而緩慢增大；同時隨著爆距的增大，α 逐漸減小，流體靜壓力對Δ的影響逐漸增大。

從前述研究中可以看到，隨著水深的增大，采用方式Ⅱ和方式Ⅲ時Δ會逐漸增大。從圖7 中可以看到，5 km 水深條件下當(dāng)7≤≤ 65時，隨著爆距的增大，水深對Δ的影響也在逐漸增大。圖8為= 5 km 時3 種方式下Δ計算結(jié)果相對于式（12）的誤差。綜合來看，采用方式Ⅱ和方式Ⅲ獲得的計算結(jié)果與試驗值及理論值均較為吻合，方式Ⅲ的相對誤差更小。

圖7 ΔPm 在5 km 水深處相對于在0 m 處的變化幅度Fig. 7 Changing amplitudes of ΔPm when depth changes from 0 m to 5 km

圖8 H=5 km 時ΔPm 相對于計算式(8)的誤差Fig. 8 Relative errors of ΔPm between simulation results and calculation formulas when H=5 km

4.2 氣泡脈動參數(shù)對比分析

Cole基于不可壓縮理論模型推導(dǎo)了和與水深和裝藥量的關(guān)系，Swift 等基于試驗數(shù)據(jù)給出了計算氣泡脈動參數(shù)的經(jīng)驗公式：

圖9 氣泡脈動參數(shù)相對于經(jīng)驗值的誤差Fig. 9 Relative errors of Rmax and T between simulation and empirical formula

4.3 設(shè)置方式合理性探討

從熱力學(xué)角度看，2.1 節(jié)中已指出，采用方式Ⅰ時，流體靜壓力的升高單純由外界傳熱引起，這與真實情形不符；采用方式Ⅱ和方式Ⅳ時，溫度隨水深的增大而略有增大，與實際情況有一定差距。采用方式Ⅲ時，水介質(zhì)溫度不隨水深變化，與模擬深水試驗條件吻合，與實際深水環(huán)境的吻合性較方式Ⅰ、Ⅱ和Ⅳ要好。從與已有成果的對比分析看，方式Ⅱ、Ⅲ的結(jié)果與已有成果吻合較好；相較而言，采用方式Ⅲ時載荷值的相對誤差最小。綜上分析認(rèn)為，初始參數(shù)設(shè)置宜選用方式Ⅲ，方式Ⅱ和方式Ⅳ次之，在討論水深影響時不建議選用方式Ⅰ。

5 結(jié) 論

分析了4 種水介質(zhì)初始參數(shù)設(shè)置方式對應(yīng)的熱力學(xué)過程，并采用有限元數(shù)值仿真方法，深入探討了初始參數(shù)設(shè)置方式對水下爆炸載荷特性的影響，并與相關(guān)試驗與理論研究成果進(jìn)行了對比分析，主要結(jié)論如下。

（1）根據(jù)參考狀態(tài)下水介質(zhì)參數(shù)，選擇式（4）形式的狀態(tài)方程作為水介質(zhì)狀態(tài)方程；其他常用狀態(tài)方程均忽略了參考壓力項，無法滿足參考狀態(tài)下的水介質(zhì)參數(shù)。

（2）僅改變水介質(zhì)內(nèi)能增量e，流體壓力源于外界傳熱，與實際環(huán)境不符；僅改變壓縮比時，水介質(zhì)溫度隨水深的增大而略有增加；LS-DYNA 中INITIAL_EOS_ALE 關(guān)鍵字給出的初始參數(shù)計算結(jié)果與等內(nèi)能假定非常接近；與上述3 種方式相比，本文提出的等溫假定更符合真實深水環(huán)境以及模擬深水試驗條件。

（3）按等內(nèi)能假定和等溫假定設(shè)置初始壓力時，水下爆炸載荷特性結(jié)果接近，與相關(guān)試驗與理論研究成果吻合良好；綜合分析實際深水環(huán)境以及模擬深水試驗條件，在數(shù)值仿真計算的初始壓力設(shè)置時，宜選用本文提出的水介質(zhì)參數(shù)等溫變化假定。