復(fù)雜建筑結(jié)構(gòu)環(huán)境下沖擊波超壓的快速估算方法

陳鑫， 高軒能， 付詩(shī)琦

(1. 華僑大學(xué) 土木工程學(xué)院， 福建 廈門 361021； 2. 福建農(nóng)林大學(xué) 交通與土木工程學(xué)院， 福建 福州 350002)

爆炸荷載的破壞力巨大，快速評(píng)估爆炸襲擊或意外爆炸對(duì)建筑結(jié)構(gòu)及人員安全的影響具有重要的意義.以往的研究已經(jīng)總結(jié)了無(wú)約束環(huán)境下的爆炸沖擊波超壓的經(jīng)驗(yàn)公式，但當(dāng)在城市街區(qū)中或大型結(jié)構(gòu)內(nèi)等復(fù)雜三維環(huán)境下發(fā)生爆炸時(shí)，爆炸點(diǎn)高度、地面剛度、地面成坑情況、建筑表面形狀、墻面和窗戶的剛度和破損程度等諸多不確定因素都使爆炸荷載變得極為復(fù)雜，經(jīng)驗(yàn)公式無(wú)法適用.因此，必須尋求一種在復(fù)雜建筑結(jié)構(gòu)環(huán)境下，爆炸沖擊波超壓荷載分布的快速計(jì)算方法.

楊亞?wèn)|等[1]以結(jié)構(gòu)壁面作為鏡面，把反射波按鏡像爆炸源的入射波計(jì)算，并疊加各壁面的鏡像入射波，快速獲得了長(zhǎng)方體房間內(nèi)的內(nèi)爆炸沖擊波.但這一計(jì)算方法必須將壁面假設(shè)為剛性，且適用于較為簡(jiǎn)單的幾何空間區(qū)域.另外，由于進(jìn)行足尺模型爆炸試驗(yàn)的成本昂貴，而縮尺模型試驗(yàn)受到材料的高應(yīng)變率效應(yīng)的影響，存在不完全相似關(guān)系，與實(shí)際情況存在偏差[2].因此，采用數(shù)值方法建立流固耦合模型模擬，是目前分析復(fù)雜環(huán)境內(nèi)超壓荷載分布的主要手段.

為了獲得較高的計(jì)算精度，流固耦合模型需要大量精細(xì)的任意拉格朗日歐拉(ALE)網(wǎng)格[3]，計(jì)算成本高.LS-DYNA有限元軟件提供的ALE與LBE相結(jié)合的方法[4-5]，以及AUTODYN軟件提供的重映射(REMAP)技術(shù)[6]，都是將爆炸荷載分成簡(jiǎn)單自由場(chǎng)傳播和流固耦合兩階段計(jì)算，節(jié)約了計(jì)算成本.這兩種方法適用于爆炸源周圍空曠無(wú)干擾，需要流固耦合分析的目標(biāo)區(qū)域小的情況，若要對(duì)較大的耦合場(chǎng)進(jìn)行分析時(shí)，依然需要大量的ALE網(wǎng)格.鑒于降低計(jì)算成本和提高計(jì)算精度的矛盾，本文提出采用粗網(wǎng)格進(jìn)行低成本的流固耦合模擬，初步獲得爆炸荷載場(chǎng)，再利用基于經(jīng)驗(yàn)公式擬合的修正公式進(jìn)行修正，獲得精度較高的超壓估算方法.

圖1 自由場(chǎng)空爆模型Fig.1 Model of free-field explosion

1 有限元模型與材料參數(shù)

估算方法的核心在于擬合出超壓修正公式，為確定超壓修正公式的適用范圍，利用LS-DYNA有限元軟件，建立自由場(chǎng)空爆ALE模型，根據(jù)對(duì)稱性建立1/8模型.模型的計(jì)算空間尺寸為4 m×4 m×15 m(長(zhǎng)×寬×高)的空氣域，如圖1所示.利用“體積分?jǐn)?shù)”語(yǔ)句在空氣域中分割一個(gè)球形的區(qū)域，用以填充炸藥材料.在模型的3個(gè)對(duì)稱面上施加對(duì)稱約束，其余表面均采用透射邊界以模擬無(wú)限空氣域.

炸藥選用*MAT_HIGH_EXPLOSIVE_BURN材料模型，并用Jones-Wilkins-Lee狀態(tài)方程描述爆炸物的特性，表達(dá)式為

(1)

式(1)中:P為壓力;V為體積；A,B,R1,R2,ω為狀態(tài)方程參數(shù)，通過(guò)試驗(yàn)確定參數(shù)的取值；E0為炸藥的單位體積初始內(nèi)能.炸藥的材料參數(shù)[7],如表1所示.表1中：v為爆速；PCJ為爆壓.炸藥密度ρ取1 630 kg·m-3.

表1 炸藥的材料參數(shù)Tab.1 Material parameters of explosive

空氣選用*MAT_NULL材料模型，方程選用*EOS_LINEAR_POLYNOMIAL，即

P=C0+C1μ+C2μ2+C3μ3+(C4+C5μ+C6μ2)Ea.

(2)

式(2)中：μ=ρa(bǔ)/ρa(bǔ)0-1，ρa(bǔ)0為空氣的參考質(zhì)量密度，ρa(bǔ)為空氣密度，取1.290 kg·m-3；單位體積內(nèi)能Ea取2.5×105J·m-3；實(shí)常數(shù)C0=C1=C2=C3=C6=0，C4=C5=0.4.

2 超壓誤差和網(wǎng)格尺寸經(jīng)濟(jì)性分析

利用節(jié)1中自由場(chǎng)空爆模型對(duì)超壓模擬誤差進(jìn)行分析，并進(jìn)一步定義網(wǎng)格尺寸的粗細(xì)程度.為了便于分析，引入網(wǎng)格密度λ,即

λ=r/b.

(3)

2.1 沖擊波超壓的誤差分析

取b=0.05 m，W=100 kg進(jìn)行計(jì)算,易得炸藥半徑r=0.245 m，網(wǎng)格密度λ=4.90.將文中模擬得到的爆炸沖擊波超壓峰值ΔPM與文獻(xiàn)[8-11]中的入射超壓經(jīng)驗(yàn)公式計(jì)算值進(jìn)行對(duì)比，結(jié)果如表2所示.表2中：ΔPH為根據(jù)Henrych等[8]的試驗(yàn)公式計(jì)算出的超壓值，有

(4)

(5)

ΔPY為楊濤春等[10]統(tǒng)計(jì)18組入射波超壓公式獲得的回歸公式計(jì)算值；ΔPB為根據(jù)GB 6722—2014《爆破安全規(guī)程》[11]公式計(jì)算的超壓值.

由表2可知：在4組經(jīng)驗(yàn)公式中，ΔPH的數(shù)值最小，且與ΔPM最為接近,ΔPW最為保守.

表2 文中模擬值與經(jīng)驗(yàn)公式計(jì)算值的對(duì)比Tab.2 Comparison of overpressure values by numerical simulation and empirical formula

圖2 比例距離與超壓模擬誤差曲線(λ=4.90)Fig.2 Curves of scale distance and overpressure simulation error (λ=4.90)

圖3 網(wǎng)格密度對(duì)超壓模擬誤差的影響Fig.3 Influence of mesh density on overpressure simulation error

圖4 炸藥質(zhì)量對(duì)超壓模擬誤差的影響Fig.4 Influence of explosive weight on overpressure simulation error

對(duì)于相同的網(wǎng)格尺寸，炸藥質(zhì)量W的變化會(huì)改變網(wǎng)格密度λ，從而影響超壓誤差δ.分析炸藥質(zhì)量W對(duì)δ的影響，結(jié)果如圖4所示.由圖4可知：不同炸藥質(zhì)量下，相同網(wǎng)格密度的超壓誤差曲線基本重合，說(shuō)明只要保證相同的網(wǎng)格密度，可以忽略炸藥質(zhì)量對(duì)超壓誤差的影響.

2.2 網(wǎng)格尺寸經(jīng)濟(jì)性分析

圖5 網(wǎng)格密度限值與比例距離的擬合曲線Fig.5 Fitting curve of mesh density limit and scale distance

(6)

由此，定義λ<λlim的網(wǎng)格尺寸為粗網(wǎng)格；反之，則為細(xì)網(wǎng)格.采用細(xì)網(wǎng)格對(duì)3D環(huán)境進(jìn)行模擬，所需要的單元數(shù)量十分巨大.以裝藥量860 kg的“飛毛腿B”導(dǎo)彈為例，殺傷半徑達(dá)到150 m.若建立尺寸為300 m×300 m×150 m(長(zhǎng)×寬×高)的爆炸區(qū)域，采用由λlim計(jì)算的0.2 m細(xì)網(wǎng)格，至少需要16.875億個(gè)單元，還未計(jì)入城市建筑模型的單元數(shù)量.由于細(xì)網(wǎng)格的計(jì)算成本極高，在計(jì)算資源有限的情況下，可以采用粗網(wǎng)格模擬計(jì)算，并通過(guò)修正獲得更精確的沖擊波超壓.

3 沖擊波超壓的快速估算

(7)

式(7)中：D1，D2和D3均為無(wú)量綱的擬合系數(shù).

3.1 修正系數(shù)式1

圖6 網(wǎng)格密度與擬合系數(shù)曲線Fig.6 Curves of mesh density and fitting coefficient

(8)

式(8)中:D1=2.34;D2=1.475;D3=1.265.

提取不同網(wǎng)格密度下的D1，D2和D3，得到網(wǎng)格密度與擬合系數(shù)曲線，如圖6所示.根據(jù)各數(shù)據(jù)點(diǎn)的分布規(guī)律，將D1，D2和D3擬合為與網(wǎng)格密度λ相關(guān)的函數(shù)，即

(9)

(10)

(11)

圖7 修正前、后超壓的對(duì)比Fig.7 Comparison of overpressure before and after correction

利用式(10),(11)或式(8)的修正系數(shù)，可以快速修正數(shù)值模擬的ΔPM，修正后的沖擊波超壓ΔPR表達(dá)式為

ΔPR=η×ΔPM.

(12)

以λ=0.65的算例為例，提取11個(gè)測(cè)點(diǎn)，將修正后的ΔPR與ΔPH進(jìn)行對(duì)比.修正前、后超壓的對(duì)比，如圖7所示.由圖7可知：修正后的超壓值與經(jīng)驗(yàn)公式計(jì)算值ΔPH吻合良好；分析11個(gè)測(cè)點(diǎn)的數(shù)據(jù)，按式(8)或式(10)修正后，ΔPR相對(duì)于ΔPH的最大誤差分別為4.2%和1.8%，誤差較小.

3.2 修正系數(shù)式2

由節(jié)2.1的分析可知，根據(jù)Henrych的經(jīng)驗(yàn)公式(式(4))計(jì)算出的超壓峰值偏小，存在低估爆炸沖擊波荷載的風(fēng)險(xiǎn)；而吳彥捷等[9]綜合21組經(jīng)驗(yàn)公式的超壓平均值，擬合獲得的超壓計(jì)算公式(式(5))更具有參考意義，計(jì)算出的超壓ΔPW也更為安全.因此，以ΔPW作為修正目標(biāo)，再次進(jìn)行模擬超壓ΔPM的修正擬合.與式(10)，(11)的擬合方法相同，先計(jì)算并擬合出不同網(wǎng)格密度λ下的D1，D2和D3，再將D1，D2和D3擬合為與λ相關(guān)的函數(shù)，最終獲得修正系數(shù)式為

(13)

(14)

3.3 理想剛性墻面反射超壓的估算與驗(yàn)證

在獲得修正系數(shù)η之后，即可按照式(12)修正低成本粗網(wǎng)格模型的模擬超壓，獲得較為精確的估算超壓值.為了檢驗(yàn)估算超壓的準(zhǔn)確性，模擬計(jì)算剛性墻面的反射超壓，并與理論值進(jìn)行對(duì)比.對(duì)于無(wú)限大的理想剛性墻面，李翼祺等[19]給出了正反射超壓ΔP2與入射超壓ΔP1之間的關(guān)系式,即

(15)

式(15)中:P0為標(biāo)準(zhǔn)大氣壓，取0.1 MPa.以W=100 kg為例，根據(jù)式(4)計(jì)算得到入射超壓ΔP1=0.078 4 MPa，代入式(15)，即可算得的理論上的ΔP2=0.204 MPa.

表3 理想剛性墻面反射超壓估算Tab.3 Estimation of reflection overpressure on ideal rigid wall

由表3可知：經(jīng)過(guò)超壓修正后，超壓誤差δ明顯減小，修正效果良好.利用粗網(wǎng)格(λ為0.49，1.23)模擬并修正后的ΔPR-H，可以獲得與細(xì)網(wǎng)格(λ為2.50，4.90)的ΔPM相當(dāng)?shù)哪M精度，且只需要細(xì)網(wǎng)格模型的1/125，甚至更少的單元數(shù)量，極大地降低了計(jì)算成本.

3.4 封閉空間內(nèi)爆炸的超壓估算與驗(yàn)證

為進(jìn)一步驗(yàn)證文中估算方法的有效性和準(zhǔn)確性，對(duì)試驗(yàn)?zāi)Ｐ瓦M(jìn)行超壓估算.文獻(xiàn)[20]中進(jìn)行的封閉空間內(nèi)爆炸試驗(yàn)示意圖，如圖8所示.試驗(yàn)?zāi)Ｐ筒捎?50 mm厚的鋼筋混凝土制成，內(nèi)壁布置了5個(gè)超壓測(cè)點(diǎn)(1#～5#).試驗(yàn)分別進(jìn)行多個(gè)炸藥量(W分別為75,150,200,300,500 g)的內(nèi)爆炸，起爆點(diǎn)均設(shè)在模型的幾何中心處.根據(jù)對(duì)稱性，建立1/8模型，空氣域尺寸(長(zhǎng)×寬×高)為0.8 m×0.8 m×1.6 m，以b=0.025 m進(jìn)行網(wǎng)格劃分，同時(shí)，按照?qǐng)D8布置剛性墻面，建模方式同節(jié)3.3.封閉空間內(nèi)爆炸數(shù)值模型示意圖，如圖9所示.

圖8 封閉空間內(nèi)爆炸試驗(yàn)示意圖(單位：mm) 圖9 封閉空間內(nèi)爆炸數(shù)值模型示意圖 Fig.8 Schematic of explosion test Fig.9 Schematic of explosion numerical in closed space (unit： mm) model in closed space

(a) 第一峰值(W=75 g) (b) 最大峰值(2#測(cè)點(diǎn))圖10 超壓實(shí)測(cè)值與估算值的對(duì)比Fig.10 Comparison of overpressure test values and estimated values

取W=75 g(λ=0.89)進(jìn)行模擬，并提取各測(cè)點(diǎn)超壓時(shí)程曲線的第一次峰值(ΔPM)，再分別利用式(10)，(11)和式(13)，(14)，計(jì)算修正后的超壓ΔPR-H和ΔPR-W，結(jié)果如圖10(a)所示.圖10(a)中：ΔPtest為實(shí)測(cè)超壓.由圖10(a)可知：修正結(jié)果明顯降低了超壓誤差，各測(cè)點(diǎn)的修正超壓與實(shí)測(cè)超壓吻合較好;僅1#測(cè)點(diǎn)和3#測(cè)點(diǎn)的ΔPR-H偏小，正如節(jié)3.2所述，ΔPR-H存在低估超壓的風(fēng)險(xiǎn).將1#，3#測(cè)點(diǎn)的ΔPM按照式(13)，(14)進(jìn)行修正后，兩個(gè)測(cè)點(diǎn)的ΔPR-W與ΔPtest的誤差僅為11.6%和17.2%，與細(xì)網(wǎng)格的模擬精度相當(dāng)，說(shuō)明通過(guò)兩組修正公式的互補(bǔ)可以較好地控制誤差.

封閉空間內(nèi)爆炸的超壓曲線具有多峰值的特點(diǎn)，2#測(cè)點(diǎn)實(shí)測(cè)超壓時(shí)程曲線的最大峰值均大于第一次峰值.因此，進(jìn)一步模擬并提取不同W下，2#測(cè)點(diǎn)的超壓最大峰值，將修正超壓與實(shí)測(cè)超壓進(jìn)行對(duì)比，如圖10(b)所示.由圖10(b)可知：除W=75 g的工況外，其他彈藥量下的ΔPtest均介于ΔPR-H和ΔPR-W之間，且誤差基本控制在20%以下，達(dá)到細(xì)網(wǎng)格的模擬精度要求.

綜上所述，利用文中估算方法獲得的修正超壓與理論值和試驗(yàn)實(shí)測(cè)值吻合良好.限于篇幅，文中僅驗(yàn)證了剛性壁面環(huán)境下的模擬超壓值.對(duì)于考慮壁面剛度和壁面破壞等更為復(fù)雜情況下的超壓場(chǎng)，修正公式的準(zhǔn)確性還有待進(jìn)一步驗(yàn)證.將實(shí)際壁面假設(shè)為剛性壁面進(jìn)行模擬，由于未考慮壁面的耦合吸能，會(huì)使反射超壓偏大，計(jì)算獲得的超壓場(chǎng)偏保守，因此，可以將其作為實(shí)際壁面環(huán)境下超壓場(chǎng)的上限值.

4 結(jié)論

為了確定復(fù)雜環(huán)境下的爆炸超壓場(chǎng)，提出利用粗網(wǎng)格建立低成本模型，再對(duì)模擬超壓進(jìn)行修正，以獲得高精度超壓的快速估算方法,得到以下3點(diǎn)結(jié)論.

3) 利用擬合的沖擊波超壓修正系數(shù)公式，能夠有效地對(duì)網(wǎng)格密度小于λlim的粗網(wǎng)格模型進(jìn)行超壓修正，獲得較為精確的估算超壓.通過(guò)估算超壓與理論值、試驗(yàn)實(shí)測(cè)值的對(duì)比，驗(yàn)證提出的超壓估算方法是可行的，且可以大幅度降低計(jì)算成本，為復(fù)雜環(huán)境下的超壓荷載計(jì)算提供參考.