水下爆炸沖擊波數值模擬的網格尺寸確定方法

第一作者張社榮男，教授，1960年生

通信作者王高輝男，講師，1986年生

郵箱：wanggaohui@whu.edu.cn

水下爆炸沖擊波數值模擬的網格尺寸確定方法

張社榮1，李宏璧1，王高輝1,2，孔源1

(1.天津大學水利工程仿真與安全國家重點實驗室，天津300072; 2.武漢大學水資源與水電工程科學國家重點實驗室，武漢430072)

摘要：網格尺寸對水下爆炸沖擊波傳播過程的數值模擬精度有較大影響，且不同當量下可接受的網格尺寸也有較大區別。通過研究不同比例爆炸距離處水下爆炸沖擊波峰值、比沖量對有限元網格尺寸的依賴性，以“炸藥半徑與網格尺寸之比”作為網格尺寸劃分依據，提出了一種對不同當量均具有較強適用性的網格尺寸確定方法；對數值結果的誤差分析表明，采用炸藥半徑的1/3作為數值模型的網格尺寸，其計算精度可滿足工程要求；同時通過對數值計算誤差進行非線性擬合，得到了不同網格尺寸下沖擊波峰值壓力、比沖量的誤差估計式，可為研究人員根據可接受的誤差范圍確定任意當量下合適的網格尺寸或估計數值計算的整體誤差提供參考。

關鍵詞：網格尺寸效應；峰值壓力；比沖量；炸藥量；比例爆炸距離

基金項目：國家創新研究群體科學基金(51021004)；國家自然科學基金(51379141)；天津市應用基礎與前沿技術研究計劃(13JCYBJC19400)；水利工程仿真與安全國家重點實驗室開放基金資助

收稿日期：2014-01-08修改稿收到日期：2014-04-24

中圖分類號:TV312文獻標志碼：A

A method to determine mesh size in numerical simulation of shock wave of underwater explosion

ZHANGShe-rong1,LIHong-bi1,WANGGao-hui1,2,KONGYuan1(1. State Key Laboratory of Hydraulic Engineering Simulation and Safety, Tianjin University, Tianjin 300072, China;2. State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China)

Abstract:The numerical simulation accuracy of underwater explosion shock wave propagation depends on mesh size, and the acceptable mesh size varies with charge weight. Through analyzing the mesh size effect on different scaled distances of shock wave peak pressure and impulse, taking the ratio of explosive radius to grid size as a criterion to determine mesh size, a mesh size-determining method being suitable for various charge weights was proposed. The error analysis of numerical results showed that if taking one third of the explosive radius as the mesh size, the numerical accuracy can meet the requirements of engineering; in addition, through nonlinear fitting of numerical simulation errors, the error estimation formulas for peak pressure and impulse of shock wave under different mesh sizes are derived. These error estimation formulas provided a reference for researchers to select an appropriate mesh size according to the charge weight and the acceptable error range or to estimate the accuracy of numerical results.

Key words:mesh size effect; peak pressure; impulse; charge weight; scaled distance

對水下爆炸沖擊波傳播過程進行數值模擬時，數值仿真的精度對有限元網格尺寸具有很大的依賴性[1]。通常為了提高計算精度，網格必須劃分的足夠小，且人工粘性的引入也要求計算網格劃分的尺度不能太大，否則計算會忽略沖擊波陣面信息，降低沖擊波陣面峰值壓力，從而降低計算的精度。當采用較小的精確網格尺寸時，單元數量增加，計算的最小時間步長也隨之減小，降低計算效率并增加對計算機硬件水平的要求。因此當模型尺寸較大、爆炸工況較為復雜時，受計算機硬件、軟件以及計算時間的限制，不得不采用較大的網格。

水下爆炸數值模擬中網格尺寸效應研究的核心就在于如何針對具體的爆炸問題確定合適的網格尺寸，以在確保計算精度的同時盡可能減少網格數量提高計算效率。目前，已有部分學者提出或研究過爆炸波模擬中的網格尺寸效應。如Krauthammer等[2]研究了爆炸沖擊荷載作用下鋼筋混凝土結構動力響應數值模擬中的網格尺寸、重力以及靜力荷載等的影響，研究結果表明有限元模型的網格尺寸對于爆炸荷載作用下鋼筋混凝土結構的變形和應力結果有較大影響；Shi等[1]研究了不同比例爆炸距離下空中爆炸沖擊波特征參量——沖擊波峰值壓力、波前到達時間、沖擊波沖量等對網格的依賴性，并提出了一個減小網格效應所造成數值誤差的方法；Zukas等[3]探討了網格尺寸大小、網格長寬比、網格形狀等對數值模擬精度的影響；Jayasinghe等[4]采用顯示非線性有限元程序，利用250mm的均勻網格研究了飽和砂土中樁基礎的動態響應；Luccioni等[5]討論了利用流體力學軟件模擬預測爆炸荷載時的網格尺寸效應，認為100mm的網格尺寸就可以較為精確的模擬爆炸荷載的傳播規律，而較粗的網格尺寸則僅僅可以用來定性的模擬爆炸荷載在城市復雜環境中的傳播規律；Wang等[6]采用100 mm的網格對比研究了水下和空中爆炸作用下混凝土重力壩的動力響應規律；余曉菲等[7]采用MSC.DYTRAN研究了50mm、70mm等網格尺寸下水下爆炸沖擊波峰值的計算誤差，并分析了爆炸荷載作用下加筋圓柱殼的動態響應；龔順風等[8]采用10 mm級的網格建立了鋼筋混凝土柱的數值模型，并闡述了空中近爆作用下鋼筋混凝土柱的損傷機理。可見，不同研究者在研究中使用的網格尺寸有較大差別。這一方面是由于不同問題的數值模擬可接受的誤差范圍有較大區別，另外則是因為爆炸問題的網格尺寸效應與爆炸所發生的環境介質(水體、空氣、巖土體等)有一定關系，且與炸藥量有較大關系。已有文獻對網格尺寸效應的研究往往針對某一特定當量，故提出的網格尺寸在炸藥當量改變時其使用有較大的局限性。

為提出一種對不同當量炸藥水下爆炸問題均具有較強適用性的網格尺寸確定方法，本文基于非線性顯示動力分析程序AUTODYN平臺，建立了不同炸藥量下的水下爆炸數值模型，對比分析了不同網格尺寸和不同比例爆炸距離處沖擊波峰值壓力和比沖量對有限元網格尺寸的依賴性，提出了一種對不同炸藥當量均具有較強適用性的網格尺寸確定方法，數值計算結果的誤差分析表明該方法合理可行，并通過對不同當量下的數值計算誤差進行非線性擬合給出了水下爆炸沖擊波峰值壓力、比沖量計算誤差與網格尺寸的關系曲線及誤差估計式。

1材料模型及狀態方程

1.1　炸藥

在爆轟性能的考察和計算中，炸藥爆轟產物的狀態方程必不可少，在爆轟高溫高壓條件下，知道某時刻產物的組成和各組分的熱力學參數，即可按照某種混合法則建立爆轟產物的狀態方程。高能炸藥材料采用JWL狀態方程[9-11]描述了爆轟壓力P和每單位體積內能E及相對體積V的關系：

(1)

式中：P為爆轟壓力；V為爆轟產物體積和炸藥初始體積之比；e為炸藥的初始內能；A、B、R1、R2、ω為特征參數，各參數取值如下：A=3.738×1011Pa，B=3.747×1011Pa，R1=4.15，R2=0.9，ω=0.35，D=6.93×103m/s，e=5.999×109J/m3，PCJ=2.10×1010Pa，ρ0=1.63×103kg/m3。

1.2　水

炸藥在水介質中爆炸時，在裝藥內形成高溫高壓的爆轟產物，其壓力遠遠大于周圍水介質的靜壓力，在水介質中會產生水中沖擊波和氣泡脈動現象。水介質采用多項式狀態方程[12-13]進行描述，其形式根據水所處狀態的不同而定。

當水為壓縮狀態(μ>0)時，其狀態方程為：

(2)

當水為膨脹狀態(μ<0)時，其狀態方程為：

P=T1μ+T2μ2+B0ρ0e

(3)

式中：μ=ρρ0-1，ρ0是初始密度；A1、A2、A3、B0、B1、T1和T2是由AUTODYN材料庫直接賦值的常數；e是比內能，其定義式如下：

(4)

式中，ρ和h分別是水的密度和深度，g和p0分別是重力加速度和大氣壓強。各參數取值如下：ρ0=1×103kg/m3，A1=2.2×109Pa，A2=9.54×109Pa，A3=1.457×1010Pa，B0=0.28，B1=0.28，T1=2.2×109Pa，T2=0 Pa。

2沖擊波經驗公式

Zamyshlyaev[14-15]在其關于水下爆炸問題的權威著作《Dynamics Loading in Underwater Explosion》中系統總結了當時水下爆炸問題試驗研究及數值研究的成果，對Cole[16]的研究成果作了進一步的發展，將水下爆炸的爆炸載荷分為5個階段：指數衰減階段、倒數衰減階段、倒數衰減后段、氣泡膨脹收縮階段和脈動壓力段，并給出了每個階段相應的載荷計算公式，公式直觀明了，且具有很高的計算精度，本文使用該公式驗證和分析數值計算結果的精度。水下爆炸中，峰值超壓往往比靜水壓力大幾個數量級，因此對于淺水爆炸，靜水壓力可忽略不計，峰值超壓也被簡稱為峰值壓力。沖擊波的壓力時程曲線在某一點由0迅速攀升到峰值壓力，接著又呈指數形式衰減，水下爆炸沖擊波指數衰減階段的沖擊波壓力時程曲線、峰值壓力及比沖量依Zamyshlyaev給出的經驗公式有：

P(t)=Pme-t/θ,t<θ

(5)

(6)

(7)

(8)

(9)

(10)

(11)

(12)

(13)

(14)

P0=Patm+ρgH0

(15)

(16)

式中：Pm為沖擊波峰值壓力，Pa；P(t)為沖擊波壓力，Pa；θ為沖擊波的指數衰減時間常數(s)，它是沖擊波壓力峰值Pm衰減到Pm/e所需的時間，e=2.718 28為自然對數的底；Ish為沖擊波比沖量，N·s/m2；W為TNT炸藥質量，kg；R為測點到爆心的距離，m；R0為炸藥初始半徑，m；td為波前到達時間，s；tp為沖擊波正壓作用時間，s；P0為爆心處流體靜水壓強，Pa；Patm為大氣壓，Pa；c為水中聲速，m/s；H0為爆心的初始深度，m。

3網格尺寸確定方法及其驗證

3.1　數值計算模型

對于無限水域或自由空氣中球形裝藥的爆炸問題，究其本質可作為一維問題模擬。這是因為在自由場中，球形炸藥起爆后，沖擊波會從起爆點沿著球體半徑方向呈球狀波向外傳播，因此可以根據球對稱性利用一維中心對稱模型來模擬爆炸波的傳播。利用對稱性減少問題的維數，可顯著減少網格數量，提高計算效率，這種處理使得在計算中采用mm級的精細網格成為可能，在較大的尺度跨度內取值對于增強網格尺寸效應研究的說服力是必要的。本文利用AUTODYN中的一維楔體中心對稱單元(Wedge)模擬球形炸藥在自由場中的爆炸問題。由于Zamyshlyaev計算公式不適用于比例爆炸距離小于0.4 m/kg1/3的情況，且當比例爆炸距離大于1.5 m/kg1/3時，爆炸沖擊波對一般鋼筋混凝土結構的損傷作用已較小。故本文僅在比例爆炸距離0.4~1.6 m/kg1/3范圍內布置監測點。比例爆炸距離Z的定義式為：

Z=R/W1/3

(17)

式中，R為測點到爆心的距離，m；W為TNT炸藥質量，kg。

文中所用數值計算模型如圖1所示。

圖1　自由場水下爆炸計算模型 Fig.1 Numerical model of free-field explosions in water

計算模型中起爆點位于炸藥的中心，水和炸藥均采用多物質Euler算法(Euler算法適合于描述液體和氣體的行為，爆轟產物和水可在網格間流動，故大變形或者有流動的情形不會導致網格畸變)。為了消除人工截斷邊界的影響，在模型截斷邊界處施加透射邊界條件(Transmit邊界)，使得沖擊波在人工截斷邊界處無反射。

3.2　網格尺寸確定方法的提出

為初步探究炸藥當量與網格尺寸效應之間的關系，分別采用50 mm、100 mm、200 mm的網格對炸藥量1 kg、100 kg、500 kg及1 000 kg的水下爆炸沖擊波傳播特性進行分析，并將數值計算結果與經驗公式計算值進行對比分析，對采用50 mm、100 mm、200 mm…這種確定性數值進行網格劃分所存在的問題做了簡要分析。

圖2給出了炸藥量分別為1 kg、100 kg、500 kg及1 000 kg時不同網格尺寸下水下爆炸的沖擊波峰值壓力與經驗公式的對比曲線。值得說明的是由于1 kg炸藥的半徑為52.7 mm，因此圖2(a)中沒有給出1 kg炸藥100 mm及200 mm網格下的數值計算結果。由圖2(a)~(d)對比可知，不同炸藥量的水下爆炸問題對網格尺寸大小的依耐性不同。在同一網格尺寸下，炸藥當量較大的情況下能獲得較為精確的數值結果，而炸藥當量較小時，則與經驗公式計算值偏離較大。如對于50 mm的網格，1 kg、100 kg、500 kg、1 000 kg當量時比例爆炸距離0.4 m/kg1/3處沖擊波峰值的數值解分別為53.92 MPa、145.83 MPa、157.01 MPa、160.75 MPa，與經驗公式計算值相比，計算誤差分別為69.07%、16.35%、9.93%、7.79%。

圖2　不同當量及不同網格尺寸下沖擊波峰值壓力對比 Fig.2 Comparison of shock wave peak pressures of different mesh sizes and charges

可見，水下爆炸數值模擬的網格尺寸效應與炸藥當量有較大關系，因此若采用50 mm、100 mm、200 mm…這樣的確定數值進行網格劃分，在炸藥當量不同時其使用具有較大的局限性。為便于研究炸藥當量對數值模擬網格尺寸效應的影響，本文特引入一無量綱參量λ來描述炸藥半徑與有限元網格尺寸之間的關系，并以“炸藥半徑與網格尺寸之比”作為網格尺寸劃分依據：

(18)

式中R0為藥包半徑，Le為有限元模型的網格尺寸。

在3.3、3.4小節中通過數值解與經驗值對比、數值解誤差分析兩個角度對采用“炸藥半徑與網格尺寸之比”作為網格尺寸確定方法的適用性及合理性進行了驗證分析。

3.3　水下爆炸沖擊波數值模擬網格尺寸確定及其驗證

為驗證上文提出的網格尺寸確定方法的可行性，探究其是否能適用于不同當量的情形，并在保證數值解能較為接近經驗公式計算值的前提下確定“炸藥半徑與網格尺寸之比”的合適取值。以1 kg、100 kg、500 kg及1 000 kg炸藥水下爆炸的問題為例，對該網格尺寸確定方法的適用性進行了驗證分析。

圖3分別給出了1 kg、100 kg、500 kg、1 000 kg的TNT炸藥在自由場水下爆炸的沖擊波峰值壓力與經驗公式計算值的對比曲線。從圖3可以看出，隨著比例爆炸距離的增加，沖擊波峰值壓力逐漸減小，這與圖3中利用經驗公式計算所得規律一致。此外，無論是炸藥當量較小(1kg)還是炸藥當量很大(1 000 kg)的情形，當有限元網格尺寸取為炸藥半徑的1/2倍或1倍(λ=2和λ=1)時，在比例爆炸距離較小的情況下，數值結果與經驗公式計算值偏離較大；而當有限元計算網格的尺寸取為炸藥半徑的1/3(λ=3)時，比例爆炸距離0.4~1.6 m/kg1/3范圍內數值計算所得的沖擊波峰值壓力與經驗公式計算值均較接近，繼續減小網格尺寸可使數值計算結果愈加接近經驗公式計算值，但計算量也會顯著增加。

圖3　沖擊波峰值壓力數值解與經驗值對比 Fig.3 Peak pressure comparison of numerical and empirical results

圖4給出了1 kg、100 kg、500 kg、1 000 kg炸藥水下爆炸時比例爆炸距離0.4~1.6 m/kg1/3處沖擊波比沖量數值計算值與經驗公式計算值的對比曲線。從圖4可以看出，在比例爆炸距離0.4~1.6 m/kg1/3范圍內，隨著比例爆炸距離的增加，數值計算所得沖擊波比沖量減小，這與經驗公式預測值所示規律一致。當有限元網格取為1/3倍炸藥半徑時，數值計算結果已非常接近經驗公式預測結果。由圖3、圖4可看出，當炸藥半徑與網格尺寸之比λ=3時，無論小當量1 kg、中等當量100 kg，還是大當量500 kg、1 000 kg的水下爆炸問題，沖擊波峰值壓力及比沖量數值計算結果均與經驗公式計算值比較接近，這說明①本文采用的根據“炸藥半徑與網格尺寸之比”確定有限元網格尺寸的方法具有較強的適用性，對于不同當量，采用相同的λ值劃分有限元模型的網格，均能獲得精度較為相近的數值計算結果；②采用炸藥半徑的1/3作為網格尺寸時，不同當量下數值計算所得沖擊波峰值壓力和比沖量均與經驗公式計算值差別不大。

圖4　沖擊波比沖量數值解與經驗值對比 Fig.4 Shock wave impulse comparison of numerical and Empirical results

3.4　誤差分析

為進一步給出依據λ值確定有限元模型網格尺寸時的數值計算誤差，探究當λ=3(網格尺寸為炸藥半徑的1/3)時，數值計算誤差能否滿足工程要求，并給出任意λ值下數值計算的誤差估計式，特針對沖擊波峰值及比沖量的數值計算結果進行了誤差分析。需要指出的是，因為缺乏爆炸試驗數據，本文的誤差分析是以已有沖擊波經驗公式(式(5)~(16))的計算值為準，分析數值計算結果偏離經驗公式計算值的程度。

圖5以1 kg、100 kg、500 kg、1 000 kg炸藥水下爆炸為例，給出了比例爆炸距離0.4、0.8、1.2、1.6 m/kg1/3處沖擊波峰值壓力計算誤差與λ的關系曲線。對比圖5(a)~(d)可以發現：①沖擊波峰值計算誤差隨著λ的增大而減小，也即網格尺寸越精細，數值計算誤差越小；②采用相同的λ值進行網格劃分，不同當量下數值計算誤差基本一致，說明依據λ值確定有限元網格尺寸的方法對不同當量炸藥的爆炸問題均具有較強的適用性，且能獲得精度大致相當的數值結果；③比例爆炸距離較大處峰值壓力的計算誤差大于比例爆炸距離較小處，這是由于沖擊波從近場傳播至遠場過程中由于不精細的網格而過快衰減誤差不斷積累所致；而且在λ較小的情況下(見曲線的彎曲段)，比例爆炸距離較大的監測點對應的誤差曲線斜率明顯小于比例爆炸距離較小的監測點，說明比例爆炸距離越遠，沖擊波峰值對網格的依賴性越小。因此，數值計算中可在近場區域采用較密網格而遠場區域采用較大網格——近場區使用較密網格可避免沖擊波在近場區過快衰減進而影響遠場區的計算精度，利用遠場區計算結果對網格依賴性較小的特點適當擴大網格以減小整體計算量——這正是爆炸問題數值模擬中常用的網格劃分策略。

圖5　沖擊波峰值壓力計算誤差分析 Fig.5 Calculation error analysis of shock wave peak pressure

圖6以1 kg、100 kg、500 kg、1 000 kg炸藥水下爆炸為例，給出了比例爆炸距離0.4、0.8、1.2、1.6 m/kg1/3處沖擊波比沖量計算誤差與λ的關系曲線。對比圖5和圖6可以發現：不同比例爆炸距離處，沖擊波峰值壓力的計算誤差有一定區別，而比沖量的計算誤差沒有明顯區別。在比沖量計算誤差與λ關系曲線的彎曲段，不同比例爆炸距離所對應曲線的斜率幾乎一樣，說明不同比例爆炸距離處比沖量的網格依賴性無顯著區別。結合圖5和圖6可看出，當網格尺寸定為炸藥半徑的1/3時，沖擊波峰值壓力的計算誤差約為35%，而比沖量的計算誤差可控制在15%左右，這種計算精度可滿足工程要求(實際中還可根據上文提到的在近場區域適當加密網格的辦法提高峰值壓力的計算精度)，說明采用炸藥半徑的1/3作為有限元網格尺寸在工程中是合理可行的。

前文已對不同網格尺寸下比例爆炸距離處0.4、0.8、1.2、1.6 m/kg1/3處沖擊波比沖量計算誤差與λ的關系進行了分析，得到了計算誤差隨比例爆炸距離變化而變化的一般規律。但是在具體數值計算中，研究者往往更關注于計算的整體誤差而不是若干特定點的計算誤差。由于不同比例爆炸距離處沖擊波峰值壓力的計算誤差有一定區別，為了評價某一λ值下數值計算的整體誤差，本文以各比例爆炸距離處沖擊波峰值計算誤差的平均值作為相應λ值下數值計算的整體誤差，比沖量的計算誤差也做了同樣處理。基于這種處理，采用非線性擬合方式獲得了爆炸沖擊波峰值壓力及比沖量整體計算誤差與λ的關系曲線及相應的誤差估計式。圖中擬合曲線所用數據點為1 kg、50 kg、100 kg、250 kg、500 kg、1 000 kg炸藥水下爆炸不同網格劃分方式下的整體計算誤差。

圖6　比沖量計算誤差分析 Fig.6 Calculation error analysis of shock waveimpulse

圖7　沖擊波峰值壓力及比沖量計算誤差擬合曲線 Fig.7 Fitting curves of numerical calculation errors about shock wave pressure and impulse

圖7(a)給出的是沖擊波峰值計算整體誤差Ep的擬合曲線，相應的誤差估計公式為：

Ep=69.699 46·λ-0.651 51

(19)

圖7(b)給出了沖擊波比沖量計算整體誤差Ei的擬合曲線，相應的誤差估計公式為：

Ei=31.610 78·λ-0.641 76

(20)

采用非線性擬合方式獲得的沖擊波峰值及比沖量整體計算誤差擬合曲線所得的相關系數的平方分別為0.988 53和0.995 43，均在0.95以上，說明曲線擬合精度和可信度較高。利用這些誤差估計式可方便的估算特定網格尺寸和炸藥當量下數值模擬的精度，這可為研究者根據可接受的誤差限值快速地確定所需要的網格尺寸或根據具體的有限元模型估計整體計算誤差提供參考。

4結論

本文基于非線性顯示動力分析程序AUTODYN平臺，建立了不同炸藥量、不同網格尺度下的水下爆炸數值模型，對不同比例爆炸距離處爆炸沖擊波峰值、比沖量的網格尺寸效應進行了敏度分析和計算結果的誤差分析，主要得到以下結論：

(1)網格尺寸對水下爆炸沖擊波傳播過程的數值模擬精度有較大影響。在一定的數值計算精度要求下，炸藥當量越大，可使用的網格尺寸越大；

(2)本文以“炸藥半徑與網格尺寸之比”作為網格尺寸劃分依據，提出了一種對不同當量均具有較強適用性的網格尺寸確定方法。對水下爆炸的數值模擬，當網格尺寸定為炸藥半徑的1/3時，計算精度可滿足工程要求；對數值結果的誤差分析也表明，當不同數值模型的“炸藥半徑與網格尺寸之比”相等時，不同當量情況下可得到數值精度相當的結果；

(3)本文以各比例爆炸距離處計算誤差的平均值作為該網格尺寸下數值計算的整體誤差，通過非線性擬合得到了沖擊波峰值及比沖量計算誤差的誤差估計式，可為研究人員根據可接受的誤差范圍確定任意當量下合適的網格尺寸或根據計算模型網格尺寸估計數值計算的整體誤差提供參考。

參考文獻

[1]Yanchao S, Zhongxian L I, Hong H. Mesh size effect in numerical simulation of blast wave propagation and interaction with structures[J]. Transactions of Tianjin University,2008, 14(6)：396-402.

[2]Krauthammer T, Otani R K. Mesh, gravity and load effects on finite element simulations of blast loaded reinforced concrete structures [J]. Computers & structures,1997, 63(6): 1113-1120.

[3]Zukas J A, Scheffler D R. Practical aspects of numerical simulations of dynamic events: effects of meshing[J]. International Journal of Impact Engineering,2000, 24(9): 925-945.

[4]Jayasinghe L B, Thambiratnam D P, Perera N, et al. Computer simulation of underground blast response of pile in saturated soil [J]. Computers & Structures, 2013, 120: 86-95.

[5]Luccioni B, Ambrosini D, Danesi R. Blast load assessment using hydrocodes [J]. Engineering Structures,2006, 28(12): 1736-1744.

[6]Wang G, Zhang S, Kong Y, et al. A comparative study on the dynamic response of concrete gravity dams subjected to underwater and air explosions[J]. Journal of Performance of Constructed Facilities, 2014.

[7]余曉菲，劉土光，張濤，等. 水下爆炸載荷作用下加筋圓柱殼的響應分析[J]. 振動與沖擊, 2006, 25(5):106-115.

YU Xiao-fei, LIU Shi-guang, ZHANG Tao,et al. Response analysis of cylindrical shell subjected to underwater explosion shock waves [J]. Journal of Vibration and Shock,2006, 25(5):106-115.

[8]龔順風，夏謙，金偉良. 近爆作用下鋼筋混凝土柱的損傷機理研究[J]. 浙江大學學報(工學版), 2011(8): 1405-1410.

GONG Shun-feng, XIA Qian, JIN Wei-liang. Investigation on damage mechanism of RC column under close-in explosion [J]. Journal of Zhejiang University (Engineering Science),2011(8): 1405-1410.

[9]John B. Numerical modelling of shock wave and pressure pulse generation by underwater explosions[J]. NASA,1998(19990007753).

[10]Larcher M. Simulation of the effects of an air blast wave[J]. JRC Technical Note, Pubsy JRC41337, Ispra, 2007.

[11]石磊，杜修力，樊鑫. 爆炸沖擊波數值計算網格劃分方法研究[J]. 北京工業大學學報, 2010(11): 1465-1470.

SHI Lei, DU Xiu-li, FAN Xin. A study on the mesh generation method for numerical simulation of blast wave [J]. Journal of Beijing University of Technology,2010(11):1465-1470.

[12]Barras G, Souli M, Aquelet N, et al. Numerical simulation of underwater explosions using an ALE method. The pulsating bubble phenomena[J]. Ocean Engineering, 2012, 41: 53-66.

[13]胡毅亭，賈憲振，饒國寧，等. 水下爆炸沖擊波和氣泡脈動的數值模擬研究[J]. 艦船科學技術，2009(2): 134-140.

HU Yi-ting, JIA Xian-zhen, RAO Guo-ning, et al. Numerical study of underwater explosion shock wave and bubble pulse [J]. Ship Science and Technology,2009(2): 134-140.

[14]劉科種，徐更光，辛春亮，等. AUTODYN水下爆炸數值模擬研究[J]. 爆破, 2009(3): 18-21.

LIU Ke-zhong, XU Geng-guang, XIN Chun-liang,et al. Research on numerical simlnation in underwater explosion by Autodyn [J].Blasting,2009(3): 18-21.

[15]Zamyshlyaev B V, Yakovlev Y S. Dynamic loads in underwater explosions[M]. Sudostroenie, Leningrad, 1973.

[16]Cole R H. Underwater explosions [M]. New York: Dover Publications, 1965.