結(jié)構(gòu)爆炸毀傷的浸沒(méi)多介質(zhì)有限體積物質(zhì)點(diǎn)法1)

倪銳晨 孫梓賢 李家盛 張雄

(清華大學(xué)航天航空學(xué)院,北京 100084)

引言

結(jié)構(gòu)爆炸毀傷現(xiàn)象廣泛存在于國(guó)防科技、公共安全和國(guó)民生產(chǎn)等領(lǐng)域中,例如結(jié)構(gòu)抗爆設(shè)計(jì)、爆炸事故分析、建筑物定向爆破等.炸藥爆炸是典型的多介質(zhì)流動(dòng)問(wèn)題,涉及爆炸產(chǎn)物和傳播介質(zhì)兩種物質(zhì).爆炸發(fā)生后,伴隨著劇烈的化學(xué)反應(yīng)和能量轉(zhuǎn)化,爆炸產(chǎn)物急劇膨脹與周?chē)膫鞑ソ橘|(zhì)相互作用,推動(dòng)介質(zhì)產(chǎn)生強(qiáng)烈的壓縮進(jìn)而形成沖擊波向外傳播.沖擊波還會(huì)進(jìn)一步與固體結(jié)構(gòu)相互作用,使結(jié)構(gòu)產(chǎn)生變形和毀傷.實(shí)驗(yàn)研究能夠反映真實(shí)的物理過(guò)程,但可重復(fù)性差、成本高,通常只能得到一些經(jīng)驗(yàn)結(jié)論或經(jīng)驗(yàn)公式.由于實(shí)際工程結(jié)構(gòu)復(fù)雜的幾何外形和非線性的材料本構(gòu)關(guān)系,根據(jù)特定實(shí)驗(yàn)得到的經(jīng)驗(yàn)結(jié)論不具有普遍性.因此,針對(duì)這類(lèi)極端工程問(wèn)題開(kāi)發(fā)有效且高效的流固耦合算法一直都是熱門(mén)的研究領(lǐng)域.

根據(jù)物質(zhì)界面的處理方式,多介質(zhì)流體求解器主要分為兩類(lèi):精確界面方法(sharp interface method,SIM)和彌散界面方法(diffused interface method,DIM)[1].SIM 需要跟蹤和重構(gòu)物質(zhì)界面,并在界面處引入額外的控制方程處理密度和內(nèi)能等熱力學(xué)量的強(qiáng)間斷,包括流場(chǎng)體積分?jǐn)?shù)(volume of fluid,VOF)方法[2-5]和基于水平集函數(shù)的虛擬流體方法(ghost fluid method,GFM)[6-8].DIM 將物質(zhì)界面處理成彌散區(qū)域,通過(guò)體積分?jǐn)?shù)加權(quán)計(jì)算相應(yīng)的熱力學(xué)量[9-12].與SIM 相比,DIM 全場(chǎng)采用統(tǒng)一的控制方程和數(shù)值方法,能夠自然地處理新生的物質(zhì)界面.因此,本文采用DIM 進(jìn)行多介質(zhì)流動(dòng)現(xiàn)象的模擬.

結(jié)構(gòu)毀傷模擬涉及固體結(jié)構(gòu)的極端變形和破壞破碎.拉格朗日方法受限于嚴(yán)重畸變的網(wǎng)格,而歐拉方法則難以記錄材料變形歷史相關(guān)的內(nèi)變量,都不能很好地模擬建筑物的結(jié)構(gòu)毀傷現(xiàn)象.Harlow[13-14]結(jié)合拉格朗日方法和歐拉方法的思想,提出了質(zhì)點(diǎn)網(wǎng)格(particle-in-cell,PIC)方法,避免了上述網(wǎng)格畸變和內(nèi)變量記錄的數(shù)值困難.Sulsky等[15-16]將FLIP PIC方法引入到固體力學(xué)領(lǐng)域,提出了物質(zhì)點(diǎn)法(material point method,MPM)[17-18].近30年來(lái),物質(zhì)點(diǎn)法在極端變形問(wèn)題的模擬中表現(xiàn)優(yōu)異,例如超高速碰撞[19-22]、靶板侵徹[23-24]、斷裂裂紋演化[25-26]、流體流動(dòng)[27-29]以及地質(zhì)滑坡失效[30-32]等.

浸沒(méi)邊界法(immersed boundary method,IBM)在流固耦合模擬中拋棄了傳統(tǒng)的貼體網(wǎng)格,在涉及固體結(jié)構(gòu)大變形及復(fù)雜幾何外形的問(wèn)題模擬中表現(xiàn)優(yōu)異[33].Mittal等[34]根據(jù)流固界面處邊界條件的施加方式將浸沒(méi)邊界法分為兩類(lèi):連續(xù)力浸沒(méi)邊界法和離散力浸沒(méi)邊界法.其中連續(xù)力浸沒(méi)邊界法[35-40]采用罰函數(shù)法在描述流固界面的拉格朗日質(zhì)點(diǎn)上計(jì)算流固耦合作用力,并通過(guò)離散近似的狄拉克 δ-函數(shù)將其彌散到流固界面附近的歐拉網(wǎng)格格心作為控制方程的源項(xiàng).離散力浸沒(méi)邊界法在流固界面處對(duì)流場(chǎng)施加速度邊界條件、對(duì)固體結(jié)構(gòu)施加力邊界條件,根據(jù)速度邊界條件的施加方式可以進(jìn)一步分為:虛擬網(wǎng)格方法[41-45]和網(wǎng)格切割方法[46-48].

對(duì)于固體結(jié)構(gòu)存在大變形乃至極端變形的流固耦合問(wèn)題,近年來(lái)基于浸沒(méi)邊界法和無(wú)網(wǎng)格法的流固耦合算法也得到了初步的應(yīng)用.Gilmanov等[44-45]將混合浸沒(méi)邊界法(hybrid immersed boundary method,HIBM)與物質(zhì)點(diǎn)法相結(jié)合,針對(duì)魚(yú)和浮游橈足類(lèi)生物在水中的游動(dòng)現(xiàn)象進(jìn)行了模擬與研究.張雄等[49]提出了基于拉格朗日乘子的連續(xù)力浸沒(méi)邊界法(continuous-forcing immersed boundary method based on Lagrangian multiplier,lg-CFIBM),并將相應(yīng)的lg-CFIBM-MPM 耦合算法應(yīng)用于單介質(zhì)流場(chǎng)中激波與結(jié)構(gòu)的相互作用及結(jié)構(gòu)動(dòng)態(tài)斷裂問(wèn)題的模擬.

本文采用基于黎曼求解器的多介質(zhì)有限體積法作為流體求解器,采用物質(zhì)點(diǎn)法作為固體求解器,并將lg-CFIBM 擴(kuò)展到多介質(zhì)流體中處理流固耦合邊界條件,建立了浸沒(méi)多介質(zhì)有限體積物質(zhì)點(diǎn)法(immersed multi-material finite volume-material point method,iMMFV-MPM)的流固耦合求解器,對(duì)近場(chǎng)爆炸下方形鋼筋混凝土靶板的失效模式、外爆載荷下建筑物的毀傷現(xiàn)象以及多腔室內(nèi)爆炸試驗(yàn)進(jìn)行了模擬,模擬結(jié)果與相關(guān)實(shí)驗(yàn)數(shù)據(jù)吻合良好,驗(yàn)證了當(dāng)前建立的流固耦合算法的有效性及精度.

1 流體求解器

1.1 流體控制方程

針對(duì)炸藥爆炸的模擬,相應(yīng)的流體控制方程為多組分可壓縮歐拉方程

通量G,H的具體表達(dá)式與F類(lèi)似.ρk為流體組分k的密度,αk為流體組分k的體積分?jǐn)?shù)

為流場(chǎng)的單位質(zhì)量總能量,ek為流體組分k的內(nèi)能,Yk=αkρk/ρ為流體組分k的質(zhì)量分?jǐn)?shù).

為了封閉上述的多組分可壓縮歐拉方程,還需要引入各組分的輸運(yùn)方程以及狀態(tài)方程

其中V=[u,v,w]為流場(chǎng)x-,y-和z-方向的速度,μ為保證各組分間壓強(qiáng)相等的罰函數(shù)系數(shù),本文取μ=∞.

1.2 有限體積法黎曼求解器

將守恒形式的微分控制方程(1)改寫(xiě)為有限體積法的積分形式

其中φ(x,y)為坡度限制器[50],本文取

然后,采用HLLC 三波近似求解器[51]計(jì)算單元表面Si-1/2,j,k上的黎曼問(wèn)題

對(duì)于非守恒形式的輸運(yùn)方程(6),采用算子分裂方法進(jìn)行計(jì)算:首先,允許各組分間的壓強(qiáng)存在差異進(jìn)行輸運(yùn)計(jì)算;然后,通過(guò)壓力松弛調(diào)整各組分的體積分?jǐn)?shù)保證各組分間的壓強(qiáng)相同.輸運(yùn)部分的控制方程為

上述3 個(gè)方程均由3 部分組成:時(shí)間偏導(dǎo)數(shù)項(xiàng)、守恒形式的通量項(xiàng)以及非守恒形式的速度散度項(xiàng).因此以體積分?jǐn)?shù)為例,結(jié)合式(10)給出的界面黎曼狀態(tài),輸運(yùn)方程的半離散形式為

對(duì)于半離散形式的式(7)和式(13),均采用具備TVD性質(zhì)的3 階龍格庫(kù)塔格式進(jìn)行時(shí)間積分.

1.3 壓力松弛

壓力松弛過(guò)程需要保證單元內(nèi)部的質(zhì)量、動(dòng)量和能量守恒,通過(guò)調(diào)整各組分的體積分?jǐn)?shù)、密度和內(nèi)能使得各組分的壓強(qiáng)相等.具體的控制方程為

通過(guò)推導(dǎo)化簡(jiǎn)[1]可以得到最終的松弛方程為

2 固體求解器

2.1 固體控制方程

更新拉格朗日格式下質(zhì)量、動(dòng)量和能量的守恒方程為

其中 ρ 表示當(dāng)前時(shí)刻的密度,u表示位移,σ 表示當(dāng)前構(gòu)型下的柯西應(yīng)力,b表示單位質(zhì)量的體力,ε 表示應(yīng)變張量,s表示單位質(zhì)量的熱源,q表示單位時(shí)間下單位面積的熱流密度.

大部分的固體材料本構(gòu)描述的是應(yīng)變率張量和應(yīng)力率張量之間的關(guān)系

其中 σ?為柯西應(yīng)力張量的焦曼率導(dǎo)數(shù).

為使式(16)定解,相應(yīng)的動(dòng)力學(xué)條件和初/邊值條件為

2.2 物質(zhì)點(diǎn)法

如圖1 所示,物質(zhì)點(diǎn)法將物質(zhì)區(qū)域 Ω 離散為一組拉格朗日質(zhì)點(diǎn),因此連續(xù)體密度可以近似為

圖1 物質(zhì)點(diǎn)法的空間離散Fig.1 Spatial discretization of MPM

其中np表示質(zhì)點(diǎn)總數(shù),mp表示質(zhì)點(diǎn)p的質(zhì)量,xp表示質(zhì)點(diǎn)p的位置坐標(biāo),δ表示狄拉克函數(shù).

物質(zhì)點(diǎn)法在每個(gè)時(shí)間步中的拉格朗日步進(jìn)行動(dòng)量方程的求解,質(zhì)點(diǎn)和網(wǎng)格保持相對(duì)位置固定,因此可以采用有限單元法中的形函數(shù)NI(x) 插值得到計(jì)算域的位移場(chǎng)

其中uI是網(wǎng)格節(jié)點(diǎn)的位移,ng是背景網(wǎng)格節(jié)點(diǎn)數(shù)量.

將式(20)和式(21)代入式(19)中,可以得到定義在背景網(wǎng)格節(jié)點(diǎn)上的離散動(dòng)量方程

分別是網(wǎng)格節(jié)點(diǎn)I處的內(nèi)力和外力,NIp=NI(xp)是質(zhì)點(diǎn)p處的形函數(shù),σp=σ(xp) 是質(zhì)點(diǎn)p的應(yīng)力,h是為了將弱形式(19)中最后一項(xiàng)面積分轉(zhuǎn)換為體積分而引入的假想邊界層厚度.

顯式物質(zhì)點(diǎn)法采用中心差分法的蛙跳格式進(jìn)行動(dòng)量方程的時(shí)間步推進(jìn)

2.3 物質(zhì)點(diǎn)法格式

在顯式物質(zhì)點(diǎn)法中,物質(zhì)點(diǎn)上的應(yīng)力狀態(tài)更新既可以在每個(gè)時(shí)間步開(kāi)始時(shí)進(jìn)行,也可以在每個(gè)時(shí)間步結(jié)束時(shí)進(jìn)行,分別稱(chēng)為USF(update-stressfirst)格式和USL(update-stress-last)格式[18].在改進(jìn)的USL(modified USL,MUSL)格式中,將更新后的質(zhì)點(diǎn)動(dòng)量映射回網(wǎng)格節(jié)點(diǎn),并采用該節(jié)點(diǎn)速度場(chǎng)進(jìn)行質(zhì)點(diǎn)應(yīng)力狀態(tài)的更新.

Ni等[52]指出拉格朗日屬性的質(zhì)點(diǎn)位置和歐拉屬性的網(wǎng)格間相互作用都會(huì)影響物質(zhì)點(diǎn)法的穩(wěn)定性,并且證明了USL 格式的本質(zhì)不穩(wěn)定性.Bardenhagen[53]發(fā)現(xiàn)MUSL 格式在模擬過(guò)程中總能量會(huì)緩慢耗散,而USF 格式的總能量則會(huì)逐漸增加.相應(yīng)的數(shù)值黏性使得MUSL 格式相比于USF 格式更穩(wěn)定.因此,本文將采用MUSL 格式的物質(zhì)點(diǎn)法.

3 流固耦合算法

3.1 網(wǎng)格單元分類(lèi)

為了施加流固邊界條件,需要先定位浸沒(méi)邊界的位置.如圖2 所示,處于固體結(jié)構(gòu)內(nèi)部的網(wǎng)格單元被分類(lèi)為type I,處于流場(chǎng)內(nèi)部的網(wǎng)格單元被分類(lèi)為type II,而被浸沒(méi)邊界穿過(guò)的網(wǎng)格單元被分類(lèi)為type III.

圖2 網(wǎng)格單元分類(lèi)Fig.2 Grid cell classification

固體結(jié)構(gòu)的運(yùn)動(dòng)由拉格朗日屬性的物質(zhì)點(diǎn)進(jìn)行跟蹤,而圖2 中紅線所示的浸沒(méi)邊界僅用于說(shuō)明,在實(shí)際模擬過(guò)程中并不會(huì)顯式重構(gòu)邊界.因?yàn)楫?dāng)固體結(jié)構(gòu)發(fā)生破碎后,通過(guò)三維空間內(nèi)的點(diǎn)云跟蹤并重構(gòu)新生的浸沒(méi)邊界過(guò)于復(fù)雜、不穩(wěn)定并且耗時(shí).可以從浸沒(méi)邊界的定義出發(fā)

其中 Ωp是質(zhì)點(diǎn)p占據(jù)的區(qū)域,通過(guò)物質(zhì)點(diǎn)的分布來(lái)確定網(wǎng)格單元的分類(lèi).具體步驟如下.

(1) 識(shí)別每個(gè)物質(zhì)點(diǎn)位于哪個(gè)網(wǎng)格單元內(nèi),并初步將網(wǎng)格單元分為兩類(lèi):物質(zhì)點(diǎn)單元和空單元.空單元即為流體網(wǎng)格單元(type II).

(2) 搜索每個(gè)物質(zhì)點(diǎn)單元的相鄰單元.若所有相鄰單元均為物質(zhì)點(diǎn)單元,則該物質(zhì)點(diǎn)單元為固體網(wǎng)格單元(type I);否則,該物質(zhì)點(diǎn)單元為流固兩相網(wǎng)格單元(type III).

3.2 基于拉格朗日乘子的連續(xù)力浸沒(méi)邊界法

如圖3 所示,傳統(tǒng)的連續(xù)力浸沒(méi)邊界法通過(guò)罰函數(shù)法實(shí)現(xiàn),無(wú)法在每個(gè)時(shí)間步嚴(yán)格滿(mǎn)足滑移或無(wú)滑移邊界條件,因此無(wú)法捕捉與邊界速度方向密切相關(guān)的激波結(jié)構(gòu),例如馬赫反射.傳統(tǒng)的離散力浸沒(méi)邊界法在每個(gè)時(shí)間步都需要計(jì)算流固耦合界面的法向和垂直距離,計(jì)算效率較低.并且在爆轟波驅(qū)動(dòng)下結(jié)構(gòu)動(dòng)態(tài)斷裂問(wèn)題的模擬中,固體結(jié)構(gòu)會(huì)產(chǎn)生尺寸不足一個(gè)網(wǎng)格的細(xì)小碎片,還會(huì)產(chǎn)生大量的新生流固耦合界面,這些都使得流固界面難以追蹤并重構(gòu).因此,本文將Ni等[49]提出的基于拉格朗日乘子的連續(xù)力浸沒(méi)邊界法(lg-CFIBM)擴(kuò)展到多介質(zhì)流體中.該算法在每個(gè)時(shí)間步內(nèi)嚴(yán)格滿(mǎn)足速度邊界條件的同時(shí),不再需要重構(gòu)流固耦合界面,適用于結(jié)構(gòu)爆炸毀傷問(wèn)題的模擬.

圖3 各浸沒(méi)邊界法邊界條件施加示意圖Fig.3 Diagram of applying boundary conditions of different IBMs

根據(jù)式(23)和式(24)重構(gòu)得到的背景網(wǎng)格節(jié)點(diǎn)動(dòng)量和集中質(zhì)量,固體結(jié)構(gòu)在type III 單元格心處的速度場(chǎng)可以由下式給出

其中nen是每個(gè)單元的節(jié)點(diǎn)數(shù)量.通過(guò)將type III 單元格心處的流場(chǎng)速度vc設(shè)置為浸沒(méi)邊界的外插速度,浸沒(méi)邊界上的速度邊界條件就可以在每個(gè)時(shí)間步嚴(yán)格滿(mǎn)足.

式(29) 速度邊界條件的施加方式等價(jià)于在每個(gè)type III 單元格心處施加相應(yīng)的歐拉源項(xiàng)

其中k表示龍格庫(kù)塔循環(huán)的階次,.在lg-CFIBM 中,type III 單元只記錄單元表面的數(shù)值通量,因此根據(jù)全微分公式Δ(ρv)=ρΔv+vΔρ,式(30)可以改寫(xiě)為

根據(jù)動(dòng)量守恒定律,式(31)歐拉源項(xiàng)的反作用力需要施加到固體結(jié)構(gòu)上

另外,type III 單元還扮演著在物質(zhì)點(diǎn)法中計(jì)算流固耦合表面應(yīng)力的角色.通過(guò)將式(26) 中的替換為Fp,物質(zhì)點(diǎn)法中的節(jié)點(diǎn)外力公式可以化簡(jiǎn)為

4 數(shù)值驗(yàn)證

為驗(yàn)證所建立的流固耦合算法的數(shù)值精度和有效性,本文分別模擬了二維激波與氦氣泡相互作用、近場(chǎng)爆炸下方形鋼筋混凝土靶板的失效模式、外爆載荷下建筑物的毀傷現(xiàn)象以及多腔室內(nèi)爆炸試驗(yàn)4 個(gè)算例,并將模擬結(jié)果與相關(guān)實(shí)驗(yàn)結(jié)果進(jìn)行了對(duì)比.

4.1 二維激波與氦氣泡相互作用

為了驗(yàn)證所建立的多介質(zhì)有限體積黎曼求解器具備描述多介質(zhì)流動(dòng)及激波捕捉的能力,本小節(jié)模擬了經(jīng)典的標(biāo)準(zhǔn)算例“二維激波與氦氣泡相互作用問(wèn)題”[54].

激波前后的流場(chǎng)狀態(tài)應(yīng)滿(mǎn)足激波間斷關(guān)系

其中MS為入射激波的馬赫數(shù).激波前空氣區(qū)域的無(wú)量綱參數(shù)設(shè)置為

則激波后空氣區(qū)域的無(wú)量綱參數(shù)可以由式(35)和MS=1.22計(jì)算得到

氦氣泡區(qū)域相應(yīng)的無(wú)量綱參數(shù)為

空氣和氦氣的力學(xué)行為均由理想氣體狀態(tài)方程

進(jìn)行描述,相應(yīng)的氣體常數(shù)分別為 γ1=1.4(空氣)和γ2=1.67(氦氣).

如圖4 所示,整個(gè)計(jì)算域[0,325]×[-44.5,44.5]采用單元尺寸為h=0.25 的均勻網(wǎng)格進(jìn)行空間離散.x=325處為入流邊界條件,并以式(37)描述相應(yīng)的入流流場(chǎng)狀態(tài);x=0處為出流邊界條件;其余邊界均為對(duì)稱(chēng)邊界條件.直徑D=50的氦氣泡的圓心位于(x=175,y=0) 處,而入射激波的初始間斷面位于x=225處.各個(gè)時(shí)刻的計(jì)算模擬結(jié)果與實(shí)驗(yàn)數(shù)據(jù)的對(duì)比如圖5 所示.文獻(xiàn)[54]給出了氦氣泡表面的激波演化示意圖,如圖6 所示.

圖4 二維激波與氦氣泡相互作用問(wèn)題描述Fig.4 Diagram of 2D shock-helium bubble problem

從圖5(a)和圖6 的對(duì)比中可以看到當(dāng)前的多介質(zhì)有限體積算法(MMFVM)可以很好地捕捉氦氣泡表面的激波演化,包括入射激波、折射激波、反射激波、側(cè)面激波和膨脹扇區(qū).同時(shí),從圖5 各個(gè)時(shí)刻計(jì)算模擬的紋影圖結(jié)果和實(shí)驗(yàn)照片的對(duì)比中可以看出,MMFVM能夠很好地模擬多介質(zhì)的流動(dòng)現(xiàn)象并捕捉相應(yīng)的物質(zhì)界面形狀.

圖5 二維激波與氦氣泡相互作用問(wèn)題的模擬結(jié)果(左)和實(shí)驗(yàn)結(jié)果[54](右)Fig.5 Numerical results(left) and experiment data[54](right) of 2D shock-helium bubble problem

圖6 氦氣泡表面的雙反射-折射示意圖[54]Fig.6 Schematic for twin regular reflection-refraction[54]

4.2 近場(chǎng)爆炸下方形鋼筋混凝土靶板的失效模式

如圖7 所示,Wang等[55]針對(duì)近場(chǎng)爆炸下方形鋼筋混凝土靶板的失效模式進(jìn)行了實(shí)驗(yàn)研究.本小節(jié)將分別采用iMMFV-MPM 算法和MPM 算法對(duì)該問(wèn)題進(jìn)行計(jì)算模擬,從而證明iMMFV-MPM 算法相比于MPM 算法的優(yōu)勢(shì).

圖7 實(shí)驗(yàn)裝置及鋼筋混凝土的布筋形式[55]Fig.7 Geometry setup of experiments[55]

如圖8 所示,根據(jù)問(wèn)題的對(duì)稱(chēng)性,采用 1 /4 模型進(jìn)行計(jì)算模擬.計(jì)算域尺寸L=W=H=600mm,采用單元尺寸h=8 mm 的均勻網(wǎng)格進(jìn)行離散.尺寸為550mm×500mm×40mm 的靶板距離爆炸中心R=400mm,采用體積為2 mm×2 mm×2 mm 的物質(zhì)點(diǎn)進(jìn)行離散.質(zhì)量為m=0.46 kg 的球形炸藥設(shè)置于計(jì)算域原點(diǎn)處.x=0mm,y=0mm和z=0mm 處為對(duì)稱(chēng)邊界條件,其余邊界處均為出流邊界條件.

圖8 計(jì)算模擬的幾何參數(shù)設(shè)置Fig.8 Geometry set for simulation

TNT 炸藥采用JWL 狀態(tài)方程進(jìn)行描述[55]

其中V=ρ/ρ0為相對(duì)體積,炸藥的初始密度為ρ0=1.63 mg/mm3,E為單位初始體積的內(nèi)能,相應(yīng)的初始值為E0=7000mJ/mm3,其余參數(shù)設(shè)置為A=371.2 GPa,B=3.21 GPa,R1=4.15,R2=0.95,w=0.3.

混凝土是水泥作為膠凝材料,砂石作為集料,通過(guò)與一定比例的水?dāng)嚢瓒傻牟牧?目前常用的混凝土材料本構(gòu)為Holmquist-Johnson-Cook(HJC)本構(gòu)[56],其將混凝土大變形過(guò)程中的孔洞坍縮所產(chǎn)生的效應(yīng)引入了強(qiáng)度模型中,能夠很好地描述混凝土材料在高應(yīng)變率、高應(yīng)力下的大變形和破壞行為.在本算例中,混凝土的密度為 ρs=2.44 g/cm3,楊氏模量為E=35.66 GPa,泊松比為 ν=0.2,其余參數(shù)如表1和表2 所示.

表1 鋼筋混凝土靶板算例中HJC 強(qiáng)度模型參數(shù)[56]Table 1 The material constants of HJC strength model in reinforced concrete slab simulation[56]

表2 鋼筋混凝土靶板算例中HJC 狀態(tài)方程參數(shù)[56]Table 2 The material constants of HJC EOS model in reinforced concrete slab simulation[56]

鋼筋的力學(xué)行為采用簡(jiǎn)化Johnson-Cook 模型進(jìn)行描述[55]

iMMFV-MPM 算法和MPM 算法模擬得到的幾個(gè)典型時(shí)刻的壓力云圖如圖9 所示.iMMFV-MPM能夠很好地捕捉爆轟波的傳播和反射現(xiàn)象,保持尖銳的激波面.并且在t=0.18ms 時(shí)刻,鋼筋混凝土靶板背面出現(xiàn)拉伸失效區(qū)域,從而引發(fā)層裂現(xiàn)象.然而在MPM 的模擬過(guò)程中,由于跨網(wǎng)格噪聲和數(shù)值斷裂,激波面無(wú)法維持完整的球面,并且在激波面附近產(chǎn)生大量速度過(guò)高的質(zhì)點(diǎn),導(dǎo)致t=0.08 ms 時(shí)刻爆轟波就與鋼筋混凝土靶板產(chǎn)生相互作用.這些散亂的高速質(zhì)點(diǎn)撞擊到鋼筋混凝土靶板上,產(chǎn)生大量零散的、非物理的高應(yīng)力區(qū),導(dǎo)致無(wú)法觀測(cè)到相應(yīng)的層裂現(xiàn)象,如圖9(d)所示.

圖9 iMMFV-MPM和MPM 的壓力云圖結(jié)果Fig.9 Pressure contour results of the iMMFV-MPM and the pure MPM

鋼筋混凝土靶板上的失效云圖及失效模式如圖10所示,并且在后處理顯示中將鋼筋網(wǎng)格人為移動(dòng)到靶板表面用于定位裂紋和層裂區(qū)域的位置.從圖10(a)中可以看出,靶板呈現(xiàn)中等程度的損傷并且在背面出現(xiàn)明顯的層裂剝落.層裂剝落區(qū)域的平均半徑為 1 20mm,層裂深度為靶板厚度的一半,正好暴露出鋼筋網(wǎng)絡(luò)[55].靶板背面存在一系列從層裂區(qū)域出發(fā)的主裂紋,其中水平和豎直方向上的3 條主裂紋均沿著鋼筋網(wǎng)絡(luò)分布.從圖10(a)和圖10(b)的對(duì)比中可以發(fā)現(xiàn),iMMFV-MPM 模擬捕捉到的失效模式與實(shí)驗(yàn)結(jié)果吻合良好.相應(yīng)的層裂剝落區(qū)域平均半徑約為 1 12 mm,略小于實(shí)驗(yàn)結(jié)果,并且在夾支方向上(橫向)剝落區(qū)域的尺寸略大于非夾支方向(縱向),與實(shí)驗(yàn)結(jié)果一致.同時(shí)從y=25 mm 的截面失效云圖結(jié)果中可以發(fā)現(xiàn),層裂現(xiàn)象終止于靶板厚度方向的中部.然而,在MPM 模擬結(jié)果中,整個(gè)鋼筋混凝土靶板幾乎完全失效,無(wú)法觀測(cè)到相應(yīng)的層裂現(xiàn)象,如圖10(c)所示.

圖10 方形鋼筋混凝土靶板的失效模式Fig.10 Damage modes of the reinforced concrete slab

iMMFV-MPM 模擬的鋼筋混凝土靶板中心點(diǎn)的時(shí)程位移曲線如圖11 所示,而MPM 的模擬結(jié)果由于靶板中心損傷過(guò)于嚴(yán)重而無(wú)法給出相應(yīng)的中心點(diǎn)位移曲線.Wang等[55]在實(shí)驗(yàn)中測(cè)量得到靶板中心點(diǎn)的最大位移為 3 5 mm,而iMMFV-MPM 的模擬結(jié)果為 3 5.9 mm,與實(shí)驗(yàn)結(jié)果的相對(duì)誤差為 2.6%,吻合良好.

圖11 iMMFV-MPM 模擬的鋼筋混凝土靶板中心點(diǎn)位移時(shí)程曲線Fig.11 Central deflection curve of the reinforced concrete slab by iMMFV-MPM

4.3 外爆載荷下建筑物的毀傷現(xiàn)象

如圖12(a)所示,Baylot等[57-58]對(duì)外爆載荷作用下建筑物的毀傷現(xiàn)象進(jìn)行了實(shí)驗(yàn)研究.基于問(wèn)題的對(duì)稱(chēng)性,本文建立了如圖12(b)所示的幾何模型,具體的幾何尺寸及鋼筋的布筋方式參照文獻(xiàn)[57].

圖12 實(shí)驗(yàn)裝置及計(jì)算模擬的幾何建模Fig.12 Setup of experiment and simulation

實(shí)驗(yàn)中采用 7.1 kg 的C4 炸藥,等價(jià)于8.449 kg的TNT 炸藥,TNT 炸藥的相關(guān)參數(shù)設(shè)置與式(40)相同.混凝土的力學(xué)行為采用HJC 模型[58]進(jìn)行描述,密度 ρ=2.068 g/cm3,楊氏模量E=28.7 GPa,泊松比 ν=0.19,抗壓強(qiáng)度=42 MPa,其余的HJC 模型參數(shù)與表1和表2 相同.鋼筋的力學(xué)行為采用線性強(qiáng)化彈塑性模型進(jìn)行描述,密度 ρ=7.5 g/cm3,楊氏模量E=200GPa,泊松比 ν=0.3,屈服強(qiáng)度 σy=449 MPa,極限強(qiáng)度 σb=513 MPa,失效應(yīng)變 εf=0.18.

模擬得到的幾個(gè)典型時(shí)刻的流場(chǎng)壓強(qiáng)云圖和建筑物損傷云圖如圖13 所示,z=300mm 平面上的流場(chǎng)體積分?jǐn)?shù)云圖和建筑物的壓強(qiáng)云圖如圖14 所示.炸藥爆炸形成球形的爆轟波,在t=0.5 ms 時(shí)爆轟波與下層中心柱相互作用形成高壓的反射波如圖13(a)所示,下層中心柱迎爆面上的混凝土材料在高壓爆轟波的作用下孔隙被逐漸壓實(shí)、損傷度逐漸累積如圖13(b)所示.流場(chǎng)經(jīng)過(guò)立柱和中板還會(huì)出現(xiàn)明顯的繞流現(xiàn)象,如圖13(b)和圖14(a)所示.在高壓反射波的驅(qū)動(dòng)下,中心立柱附近的爆炸產(chǎn)物出現(xiàn)了明顯的回卷現(xiàn)象,如圖14 的體積分?jǐn)?shù)云圖所示.由于中板下方的流場(chǎng)壓強(qiáng)明顯高于上方的流場(chǎng)壓強(qiáng),中板會(huì)向上隆起,表現(xiàn)為中板上表面主要分布拉應(yīng)力,如圖14所示.在本實(shí)驗(yàn)中,Baylot等[57-58]將鋼筋配置在混凝土內(nèi)部接近表面處,因此計(jì)算模擬結(jié)果中混凝土結(jié)構(gòu)上沒(méi)有出現(xiàn)層裂現(xiàn)象.而下層中心柱背爆面中部的損傷度累積則是由大變形的鋼筋擠壓混凝土引起的.下層中心柱彎曲后上下兩端連接處受到強(qiáng)剪切作用,導(dǎo)致出現(xiàn)強(qiáng)烈的毀傷現(xiàn)象如圖13(d)所示.

圖13 流場(chǎng)壓強(qiáng)云圖與建筑物損傷云圖Fig.13 Pressure contour in fluid domain and damage contour in building

圖14 流場(chǎng)體積分?jǐn)?shù)截面云圖與建筑物壓力云圖Fig.14 Volume of fraction contour on the plane of z=300mm and pressure contour in building

爆炸結(jié)束后建筑物毀傷的實(shí)驗(yàn)結(jié)果和計(jì)算模擬結(jié)果如圖15 所示.計(jì)算模擬的建筑物毀傷現(xiàn)象與實(shí)驗(yàn)結(jié)果吻合良好,下層中心柱的上下兩端連接處均出現(xiàn)了劇烈的損傷,中板的上表面都能觀測(cè)到相應(yīng)的裂紋分布.從側(cè)面損傷度云圖和實(shí)驗(yàn)側(cè)面照片的對(duì)比中可以看出,中板均存在明顯的隆起,并且中部存在一條貫穿裂紋.

圖15 最終時(shí)刻的建筑物損傷云圖Fig.15 Damage contour in building at last

4.4 多腔室內(nèi)爆炸試驗(yàn)

郭志昆等[59]對(duì)扁平箱形密閉結(jié)構(gòu)內(nèi)爆炸問(wèn)題進(jìn)行了實(shí)驗(yàn)研究,具體的幾何參數(shù)設(shè)置如圖16(a)所示.本文建立了如圖16(b)所示的幾何模型(未顯示頂板)進(jìn)行計(jì)算模擬.

圖16 實(shí)驗(yàn)和模擬的幾何參數(shù)設(shè)置Fig.16 Geometry setup of experiment and simualtion

實(shí)驗(yàn)中采用 2 20g 的TNT 炸藥,TNT 炸藥的相關(guān)參數(shù)設(shè)置與式(40)相同.實(shí)驗(yàn)采用的混凝土設(shè)計(jì)強(qiáng)度等級(jí)為C30,相關(guān)HJC 模型的材料參數(shù)來(lái)自文獻(xiàn)[60],密度 ρ=2.4 g/cm3,泊松比 ν=0.202,楊氏模量E=33.4 GPa,抗壓強(qiáng)度=39.2 MPa,最大靜水拉力T=3.162 MPa,其余參數(shù)與表1和表2 相同.實(shí)驗(yàn)采用的鋼筋為HPB235 級(jí)鋼筋,其力學(xué)行為采用理想彈塑性模型進(jìn)行描述,密度 ρ=7.8 g/cm3,楊氏模量E=220GPa,泊松比 ν=0.3,屈服強(qiáng)度 σy=235 MPa.

主爆室外側(cè)墻中心點(diǎn)處的壓強(qiáng)時(shí)程曲線如圖17 所示,爆轟波作用在側(cè)墻中心點(diǎn)處的超壓峰值為0.938 MPa,實(shí)驗(yàn)中相應(yīng)測(cè)點(diǎn)處測(cè)得的超壓峰值為0.987 MPa,相對(duì)誤差為4.96%,與實(shí)驗(yàn)結(jié)果吻合良好.幾個(gè)典型時(shí)刻中截面上的流場(chǎng)壓強(qiáng)云圖和體積分?jǐn)?shù)云圖如圖18 所示.方形炸藥包爆炸后,爆炸產(chǎn)物急劇膨脹并推動(dòng)周?chē)諝猱a(chǎn)生爆轟波如圖18(a)所示,在t=1.5 ms 時(shí),爆轟波與主爆室的側(cè)墻相互作用產(chǎn)生反射波,形成圖17 中的第1 個(gè)超壓峰值(peak 1).t=3.5 ms 時(shí),頂板與底板的反射波在中截面處相交形成二次反射波如圖18(b)所示,該二次反射波與主爆室的側(cè)墻相互作用,形成圖17 中的第2 個(gè)超壓峰值(peak 2).相鄰側(cè)墻上的爆轟反射波傳播到外側(cè)墻中心點(diǎn)處形成第3 個(gè)超壓峰值(peak 3);頂板與底板的二次反射波在相鄰側(cè)墻上反射后傳播到外側(cè)墻中心點(diǎn)處形成第4 個(gè)超壓峰值(peak 4);正對(duì)側(cè)墻上的爆轟反射波傳播到外側(cè)墻中心點(diǎn)處形成第5 個(gè)超壓峰值(peak 5).

圖17 主爆室側(cè)墻中心點(diǎn)處的壓強(qiáng)曲線Fig.17 Pressure curve at the center of sidewall in explosion room

圖18 各時(shí)刻中截面上的流場(chǎng)壓強(qiáng)云圖和體積分?jǐn)?shù)云圖Fig.18 Pressure contour and volume of fraction contour of middle plane at different instants

從圖18 各時(shí)刻的壓強(qiáng)云圖和壓強(qiáng)標(biāo)尺中可以發(fā)現(xiàn),主爆室的超壓峰值在1 MPa 左右,而側(cè)爆室的超壓峰值驟降為0.05 MPa～0.1 MPa,對(duì)角爆室的超壓峰值則為0.02 MPa 左右,這與實(shí)驗(yàn)觀測(cè)到的各爆室的超壓峰值比例一致.從各時(shí)刻的體積分?jǐn)?shù)云圖中可以發(fā)現(xiàn),爆炸產(chǎn)物基本都?xì)埩粼谥鞅覂?nèi),主/側(cè)爆室之間的內(nèi)墻開(kāi)窗方向上形成對(duì)外的高壓射流,而在4 個(gè)角區(qū)方向上形成“蘑菇傘”狀的物質(zhì)界面結(jié)構(gòu).

最終時(shí)刻頂板和主爆室側(cè)墻的破壞情況如圖19所示.從圖19(b)中可以看出,由于非爆室的超壓峰值驟降為主爆室的1/10甚至更低,破壞都集中在主爆室,而非爆室墻體仍處于彈性工作狀態(tài),與實(shí)驗(yàn)觀測(cè)結(jié)果一致.主爆室的頂板中心區(qū)域發(fā)生劇烈的貫穿破壞,其位置和大小與實(shí)驗(yàn)照片基本一致;并且大量的非貫穿裂紋從中心貫穿區(qū)域邊緣向四周延伸.從圖19(b)最右側(cè)的側(cè)墻破壞情況中可以發(fā)現(xiàn),結(jié)構(gòu)內(nèi)轉(zhuǎn)角處是薄弱部位,出現(xiàn)了明顯的通長(zhǎng)裂紋,這與實(shí)驗(yàn)觀測(cè)結(jié)果一致.

圖19 頂板和側(cè)墻的破壞情況Fig.19 Roof and sidewall fragmentation

5 結(jié)論

本文針對(duì)建筑物爆炸毀傷的流固耦合問(wèn)題,通過(guò)基于拉格朗日乘子的連續(xù)力浸沒(méi)邊界法(lg-CFIBM)將多介質(zhì)有限體積法(MMFVM)和物質(zhì)點(diǎn)法(MPM)相結(jié)合,建立了浸沒(méi)多介質(zhì)有限體積物質(zhì)點(diǎn)耦合算法(iMMFV-MPM).首先通過(guò)二維氦氣泡與激波相互作用的典型算例驗(yàn)證了上述算法具備求解多介質(zhì)流動(dòng)及激波捕捉的能力.在此基礎(chǔ)上,對(duì)近場(chǎng)爆炸下方形鋼筋混凝土靶板的失效模式、外爆載荷下建筑物的毀傷現(xiàn)象以及多腔室內(nèi)爆炸試驗(yàn)進(jìn)行了模擬,得到以下結(jié)論.

(1) 相比于傳統(tǒng)的采用純MPM 模擬流固耦合問(wèn)題,建立的iMMFV-MPM 方法在模擬炸藥爆炸產(chǎn)生爆轟波的過(guò)程中可以得到更銳利的波陣面,從而在后續(xù)爆轟波與固體結(jié)構(gòu)相互作用以及固體結(jié)構(gòu)毀傷現(xiàn)象的模擬中得到高精度的結(jié)果.

(2) 將基于單介質(zhì)流體的lg-CFIBM 流固耦合算法擴(kuò)展到多介質(zhì)流體中,能夠較好地模擬近場(chǎng)爆炸下爆炸產(chǎn)物與固體結(jié)構(gòu)間的相互作用,例如大樓外爆炸算例中,中心立柱附近爆炸產(chǎn)物的繞流和回卷現(xiàn)象.

(3) 采用iMMFV-MPM 方法首先對(duì)建筑物中最常見(jiàn)的鋼筋混凝土結(jié)構(gòu)進(jìn)行了近場(chǎng)爆炸下失效模式的模擬,能夠很好地捕捉到層裂現(xiàn)象;此外,根據(jù)爆炸載荷常見(jiàn)的兩種形式分別模擬了“外爆炸載荷下建筑物的毀傷現(xiàn)象”和“多腔室內(nèi)爆炸試驗(yàn)”,給出了爆炸產(chǎn)物和空氣的流場(chǎng)演化以及結(jié)構(gòu)的變形和毀傷特性分析,均與相應(yīng)的實(shí)驗(yàn)結(jié)果吻合良好.