水介質初始參數設置對水下爆炸載荷的影響*

鄭永輝，魏繼鋒

（北京理工大學爆炸科學與技術國家重點實驗室，北京 100081）

隨著數值計算方法的進步和計算能力的提升，數值仿真技術為解決水下爆炸機理分析、結構響應機制規律揭示等科學技術難題提供了有力的手段。在水下爆炸數值仿真研究中，研究人員已對水介質狀態方程、人工黏性系數、網格密度等影響計算精度的因素開展了較為細致的探討。

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

除直接設置水介質初始參數外，也可通過施加重力、采用程序提供的特定功能等方式實現流場壓力初始化。這些方式簡化了操作，但程序仍需提供水介質的初始參數使其具有一定的壓力。以LSDYNA 提供的關鍵字INITIAL_EOS_ALE 為例，其說明文檔指出，程序將采用一種迭代方法根據用戶提供的壓力計算介質的初始內能和比容。簡而言之，該類方法其本質上也是對水介質的初始參數進行設置。

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

1 水介質狀態方程

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

式中：、、、、、和均為狀態方程系數，e為質量內能增量，、、、、、和均為狀態方程系數，、、、、、和均為狀態方程系數。

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

表1 水介質狀態方程參數Table 1 EOS parameters of water

2 初始參數設置

2.1 設置方式及狀態參數變化

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

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

式中：= 273 K。當0 ≤≤ 0.025 時，隨的增大而逐漸增大（見圖1）。在實際深水環境中，隨著水深的增大，水溫逐漸減小；在通過加壓構造的模擬深水環境中，水溫與環境溫度基本相同，即溫度可認為保持不變。因此，從溫度變化趨勢的角度來看，方式Ⅱ與實際情況有一定的差距。

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

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

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

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

2.2 流體可壓縮性

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

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

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

3 初始參數設置對水下爆炸載荷的影響

3.1 數值仿真模型

建立了一維球對稱水下爆炸數值仿真模型，并應用LS-DYNA 數值仿真軟件進行計算。其中，炸藥采用MAT_HIGH_EXPLOSIVE_BURN 材料模型以及JWL 狀態方程，參數取自文獻[17]。水介質的狀態方程采用式（4），參數在表1 中列出。

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

3.2 沖擊波載荷特性

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

表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 種設置方式在不同水深下獲得的沖擊波超壓峰值、沖量和能流密度，其中，和分別由下式給出；

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

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

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

表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 種設置方式下氣泡的脈動曲線，如圖5 所示，相應的氣泡最大半徑和第1 次脈動周期列于表4。可以看出，除方式Ⅰ外，其余兩條曲線基本重合；方式Ⅲ、方式Ⅱ、方式Ⅰ對應的和依次增大。顯然，該規律與水介質壓縮比的變化一致。由于氣泡脈動與水介質的慣性密切相關，當體積不變時，水介質密度的增大導致其慣性增大，在密度變化較小時慣性增幅有限，因此，3 種方式下氣泡脈動曲線的差別并不明顯。

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

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

4 對比分析

4.1 沖擊波載荷對比分析

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

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

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

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

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

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

圖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 氣泡脈動參數對比分析

Cole基于不可壓縮理論模型推導了和與水深和裝藥量的關系，Swift 等基于試驗數據給出了計算氣泡脈動參數的經驗公式：

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

4.3 設置方式合理性探討

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

5 結 論

分析了4 種水介質初始參數設置方式對應的熱力學過程，并采用有限元數值仿真方法，深入探討了初始參數設置方式對水下爆炸載荷特性的影響，并與相關試驗與理論研究成果進行了對比分析，主要結論如下。

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

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

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