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

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

郵箱：wanggaohui@whu.edu.cn

水下爆炸沖擊波數(shù)值模擬的網(wǎng)格尺寸確定方法

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

(1.天津大學(xué)水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室，天津300072; 2.武漢大學(xué)水資源與水電工程科學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，武漢430072)

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

關(guān)鍵詞：網(wǎng)格尺寸效應(yīng)；峰值壓力；比沖量；炸藥量；比例爆炸距離

基金項(xiàng)目：國(guó)家創(chuàng)新研究群體科學(xué)基金(51021004)；國(guó)家自然科學(xué)基金(51379141)；天津市應(yīng)用基礎(chǔ)與前沿技術(shù)研究計(jì)劃(13JCYBJC19400)；水利工程仿真與安全國(guó)家重點(diǎn)實(shí)驗(yàn)室開放基金資助

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

中圖分類號(hào):TV312文獻(xiàn)標(biāo)志碼：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

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

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

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

1材料模型及狀態(tài)方程

1.1　炸藥

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

(1)

式中：P為爆轟壓力；V為爆轟產(chǎn)物體積和炸藥初始體積之比；e為炸藥的初始內(nèi)能；A、B、R1、R2、ω為特征參數(shù)，各參數(shù)取值如下：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　水

炸藥在水介質(zhì)中爆炸時(shí)，在裝藥內(nèi)形成高溫高壓的爆轟產(chǎn)物，其壓力遠(yuǎn)遠(yuǎn)大于周圍水介質(zhì)的靜壓力，在水介質(zhì)中會(huì)產(chǎn)生水中沖擊波和氣泡脈動(dòng)現(xiàn)象。水介質(zhì)采用多項(xiàng)式狀態(tài)方程[12-13]進(jìn)行描述，其形式根據(jù)水所處狀態(tài)的不同而定。

當(dāng)水為壓縮狀態(tài)(μ>0)時(shí)，其狀態(tài)方程為：

(2)

當(dāng)水為膨脹狀態(tài)(μ<0)時(shí)，其狀態(tài)方程為：

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

(3)

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

(4)

式中，ρ和h分別是水的密度和深度，g和p0分別是重力加速度和大氣壓強(qiáng)。各參數(shù)取值如下：ρ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沖擊波經(jīng)驗(yàn)公式

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

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；θ為沖擊波的指數(shù)衰減時(shí)間常數(shù)(s)，它是沖擊波壓力峰值Pm衰減到Pm/e所需的時(shí)間，e=2.718 28為自然對(duì)數(shù)的底；Ish為沖擊波比沖量，N·s/m2；W為TNT炸藥質(zhì)量，kg；R為測(cè)點(diǎn)到爆心的距離，m；R0為炸藥初始半徑，m；td為波前到達(dá)時(shí)間，s；tp為沖擊波正壓作用時(shí)間，s；P0為爆心處流體靜水壓強(qiáng)，Pa；Patm為大氣壓，Pa；c為水中聲速，m/s；H0為爆心的初始深度，m。

3網(wǎng)格尺寸確定方法及其驗(yàn)證

3.1　數(shù)值計(jì)算模型

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

Z=R/W1/3

(17)

式中，R為測(cè)點(diǎn)到爆心的距離，m；W為TNT炸藥質(zhì)量，kg。

文中所用數(shù)值計(jì)算模型如圖1所示。

圖1　自由場(chǎng)水下爆炸計(jì)算模型 Fig.1 Numerical model of free-field explosions in water

計(jì)算模型中起爆點(diǎn)位于炸藥的中心，水和炸藥均采用多物質(zhì)Euler算法(Euler算法適合于描述液體和氣體的行為，爆轟產(chǎn)物和水可在網(wǎng)格間流動(dòng)，故大變形或者有流動(dòng)的情形不會(huì)導(dǎo)致網(wǎng)格畸變)。為了消除人工截?cái)噙吔绲挠绊?，在模型截?cái)噙吔缣幨┘油干溥吔鐥l件(Transmit邊界)，使得沖擊波在人工截?cái)噙吔缣師o反射。

3.2　網(wǎng)格尺寸確定方法的提出

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

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

圖2　不同當(dāng)量及不同網(wǎng)格尺寸下沖擊波峰值壓力對(duì)比 Fig.2 Comparison of shock wave peak pressures of different mesh sizes and charges

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

(18)

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

在3.3、3.4小節(jié)中通過數(shù)值解與經(jīng)驗(yàn)值對(duì)比、數(shù)值解誤差分析兩個(gè)角度對(duì)采用“炸藥半徑與網(wǎng)格尺寸之比”作為網(wǎng)格尺寸確定方法的適用性及合理性進(jìn)行了驗(yàn)證分析。

3.3　水下爆炸沖擊波數(shù)值模擬網(wǎng)格尺寸確定及其驗(yàn)證

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

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

圖3　沖擊波峰值壓力數(shù)值解與經(jīng)驗(yàn)值對(duì)比 Fig.3 Peak pressure comparison of numerical and empirical results

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

圖4　沖擊波比沖量數(shù)值解與經(jīng)驗(yàn)值對(duì)比 Fig.4 Shock wave impulse comparison of numerical and Empirical results

3.4　誤差分析

為進(jìn)一步給出依據(jù)λ值確定有限元模型網(wǎng)格尺寸時(shí)的數(shù)值計(jì)算誤差，探究當(dāng)λ=3(網(wǎng)格尺寸為炸藥半徑的1/3)時(shí)，數(shù)值計(jì)算誤差能否滿足工程要求，并給出任意λ值下數(shù)值計(jì)算的誤差估計(jì)式，特針對(duì)沖擊波峰值及比沖量的數(shù)值計(jì)算結(jié)果進(jìn)行了誤差分析。需要指出的是，因?yàn)槿狈Ρㄔ囼?yàn)數(shù)據(jù)，本文的誤差分析是以已有沖擊波經(jīng)驗(yàn)公式(式(5)~(16))的計(jì)算值為準(zhǔn)，分析數(shù)值計(jì)算結(jié)果偏離經(jīng)驗(yàn)公式計(jì)算值的程度。

圖5以1 kg、100 kg、500 kg、1 000 kg炸藥水下爆炸為例，給出了比例爆炸距離0.4、0.8、1.2、1.6 m/kg1/3處沖擊波峰值壓力計(jì)算誤差與λ的關(guān)系曲線。對(duì)比圖5(a)~(d)可以發(fā)現(xiàn)：①?zèng)_擊波峰值計(jì)算誤差隨著λ的增大而減小，也即網(wǎng)格尺寸越精細(xì)，數(shù)值計(jì)算誤差越?。虎诓捎孟嗤摩酥颠M(jìn)行網(wǎng)格劃分，不同當(dāng)量下數(shù)值計(jì)算誤差基本一致，說明依據(jù)λ值確定有限元網(wǎng)格尺寸的方法對(duì)不同當(dāng)量炸藥的爆炸問題均具有較強(qiáng)的適用性，且能獲得精度大致相當(dāng)?shù)臄?shù)值結(jié)果；③比例爆炸距離較大處峰值壓力的計(jì)算誤差大于比例爆炸距離較小處，這是由于沖擊波從近場(chǎng)傳播至遠(yuǎn)場(chǎng)過程中由于不精細(xì)的網(wǎng)格而過快衰減誤差不斷積累所致；而且在λ較小的情況下(見曲線的彎曲段)，比例爆炸距離較大的監(jiān)測(cè)點(diǎn)對(duì)應(yīng)的誤差曲線斜率明顯小于比例爆炸距離較小的監(jiān)測(cè)點(diǎn)，說明比例爆炸距離越遠(yuǎn)，沖擊波峰值對(duì)網(wǎng)格的依賴性越小。因此，數(shù)值計(jì)算中可在近場(chǎng)區(qū)域采用較密網(wǎng)格而遠(yuǎn)場(chǎng)區(qū)域采用較大網(wǎng)格——近場(chǎng)區(qū)使用較密網(wǎng)格可避免沖擊波在近場(chǎng)區(qū)過快衰減進(jìn)而影響遠(yuǎn)場(chǎng)區(qū)的計(jì)算精度，利用遠(yuǎn)場(chǎng)區(qū)計(jì)算結(jié)果對(duì)網(wǎng)格依賴性較小的特點(diǎn)適當(dāng)擴(kuò)大網(wǎng)格以減小整體計(jì)算量——這正是爆炸問題數(shù)值模擬中常用的網(wǎng)格劃分策略。

圖5　沖擊波峰值壓力計(jì)算誤差分析 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處沖擊波比沖量計(jì)算誤差與λ的關(guān)系曲線。對(duì)比圖5和圖6可以發(fā)現(xiàn)：不同比例爆炸距離處，沖擊波峰值壓力的計(jì)算誤差有一定區(qū)別，而比沖量的計(jì)算誤差沒有明顯區(qū)別。在比沖量計(jì)算誤差與λ關(guān)系曲線的彎曲段，不同比例爆炸距離所對(duì)應(yīng)曲線的斜率幾乎一樣，說明不同比例爆炸距離處比沖量的網(wǎng)格依賴性無顯著區(qū)別。結(jié)合圖5和圖6可看出，當(dāng)網(wǎng)格尺寸定為炸藥半徑的1/3時(shí)，沖擊波峰值壓力的計(jì)算誤差約為35%，而比沖量的計(jì)算誤差可控制在15%左右，這種計(jì)算精度可滿足工程要求(實(shí)際中還可根據(jù)上文提到的在近場(chǎng)區(qū)域適當(dāng)加密網(wǎng)格的辦法提高峰值壓力的計(jì)算精度)，說明采用炸藥半徑的1/3作為有限元網(wǎng)格尺寸在工程中是合理可行的。

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

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

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

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

Ep=69.699 46·λ-0.651 51

(19)

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

Ei=31.610 78·λ-0.641 76

(20)

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

4結(jié)論

本文基于非線性顯示動(dòng)力分析程序AUTODYN平臺(tái)，建立了不同炸藥量、不同網(wǎng)格尺度下的水下爆炸數(shù)值模型，對(duì)不同比例爆炸距離處爆炸沖擊波峰值、比沖量的網(wǎng)格尺寸效應(yīng)進(jìn)行了敏度分析和計(jì)算結(jié)果的誤差分析，主要得到以下結(jié)論：

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

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

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

參考文獻(xiàn)

[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]余曉菲，劉土光，張濤，等. 水下爆炸載荷作用下加筋圓柱殼的響應(yīng)分析[J]. 振動(dòng)與沖擊, 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]龔順風(fēng)，夏謙，金偉良. 近爆作用下鋼筋混凝土柱的損傷機(jī)理研究[J]. 浙江大學(xué)學(xué)報(bào)(工學(xué)版), 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]石磊，杜修力，樊鑫. 爆炸沖擊波數(shù)值計(jì)算網(wǎng)格劃分方法研究[J]. 北京工業(yè)大學(xué)學(xué)報(bào), 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]胡毅亭，賈憲振，饒國(guó)寧，等. 水下爆炸沖擊波和氣泡脈動(dòng)的數(shù)值模擬研究[J]. 艦船科學(xué)技術(shù)，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水下爆炸數(shù)值模擬研究[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.