鄭永輝，魏繼鋒

（北京理工大學(xué)爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081）

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

在進(jìn)行水下爆炸問(wèn)題計(jì)算時(shí)，需設(shè)定水介質(zhì)的初始參數(shù)以設(shè)置流場(chǎng)初始?jí)毫ΑＲ延醒芯恐?，大多采用以下兩種方式：一是認(rèn)為水介質(zhì)密度變化很小，可忽略不計(jì)，僅通過(guò)改變內(nèi)能設(shè)置流場(chǎng)初始?jí)毫Γ欢钦J(rèn)為流體靜壓力的增大導(dǎo)致水的密度增大，僅通過(guò)調(diào)整比容來(lái)設(shè)置初始?jí)毫?。?dāng)采用前一種方式時(shí)，沖擊波超壓峰值隨水深的增大略有減小，后一種方式下的變化趨勢(shì)則相反?？梢钥闯?，初始參數(shù)的設(shè)置方式會(huì)帶來(lái)水下爆炸載荷計(jì)算結(jié)果的差異，而對(duì)此的細(xì)致分析尚未見(jiàn)報(bào)道。

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

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

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

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

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

與Mie-Grüneisen 狀態(tài)方程類似，式（2）～（3）也存在與實(shí)際壓力不符的問(wèn)題。式（4）在=時(shí)，可滿足= 0，e= 0 以及=。因此，宜選擇Polynomial 狀態(tài)方程式（4）作為水介質(zhì)狀態(tài)方程，具體參數(shù)見(jiàn)表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è)置初始?jí)毫r(shí)可改變參量和e。現(xiàn)今大多數(shù)仿真研究是通過(guò)改變e（方式Ⅰ）來(lái)進(jìn)行流場(chǎng)初始?jí)毫υO(shè)置。從熱學(xué)角度看，方式Ⅰ中不同水深下的水介質(zhì)參數(shù)按等容過(guò)程（d=0）變化，此時(shí)外界做功為零，流場(chǎng)的內(nèi)能增量全部由外界傳導(dǎo)的熱量提供（d= d）。這表明采用方式Ⅰ時(shí)，流體靜壓力的增大單純?cè)从谕饨鐐鳠?，與重力等外力做功等因素?zé)o關(guān)，這與實(shí)際深水環(huán)境以及加壓模擬深水環(huán)境嚴(yán)重不符。

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

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

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

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

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

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

2.2 流體可壓縮性

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

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

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

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

3.1 數(shù)值仿真模型

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

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

3.2 沖擊波載荷特性

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

表2 H =5 km 及 R =55R0 時(shí)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 時(shí)3 種方式下的沖擊波壓力-時(shí)間曲線Fig. 3 Shock wave pressure-time curves in three modes when H=5 km and R=55R0

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

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

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

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

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

3.3 氣泡脈動(dòng)特性

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

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

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

4 對(duì)比分析

4.1 沖擊波載荷對(duì)比分析

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

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

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

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

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

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

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

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

4.2 氣泡脈動(dòng)參數(shù)對(duì)比分析

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

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

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

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

5 結(jié) 論

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

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

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

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