結構爆炸毀傷的浸沒多介質有限體積物質點法1)

倪銳晨 孫梓賢 李家盛 張雄

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

引言

結構爆炸毀傷現象廣泛存在于國防科技、公共安全和國民生產等領域中,例如結構抗爆設計、爆炸事故分析、建筑物定向爆破等.炸藥爆炸是典型的多介質流動問題,涉及爆炸產物和傳播介質兩種物質.爆炸發生后,伴隨著劇烈的化學反應和能量轉化,爆炸產物急劇膨脹與周圍的傳播介質相互作用,推動介質產生強烈的壓縮進而形成沖擊波向外傳播.沖擊波還會進一步與固體結構相互作用,使結構產生變形和毀傷.實驗研究能夠反映真實的物理過程,但可重復性差、成本高,通常只能得到一些經驗結論或經驗公式.由于實際工程結構復雜的幾何外形和非線性的材料本構關系,根據特定實驗得到的經驗結論不具有普遍性.因此,針對這類極端工程問題開發有效且高效的流固耦合算法一直都是熱門的研究領域.

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

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

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

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

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

1 流體求解器

1.1 流體控制方程

針對炸藥爆炸的模擬,相應的流體控制方程為多組分可壓縮歐拉方程

通量G,H的具體表達式與F類似.ρk為流體組分k的密度,αk為流體組分k的體積分數

為流場的單位質量總能量,ek為流體組分k的內能,Yk=αkρk/ρ為流體組分k的質量分數.

為了封閉上述的多組分可壓縮歐拉方程,還需要引入各組分的輸運方程以及狀態方程

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

1.2 有限體積法黎曼求解器

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

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

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

對于非守恒形式的輸運方程(6),采用算子分裂方法進行計算:首先,允許各組分間的壓強存在差異進行輸運計算;然后,通過壓力松弛調整各組分的體積分數保證各組分間的壓強相同.輸運部分的控制方程為

上述3 個方程均由3 部分組成:時間偏導數項、守恒形式的通量項以及非守恒形式的速度散度項.因此以體積分數為例,結合式(10)給出的界面黎曼狀態,輸運方程的半離散形式為

對于半離散形式的式(7)和式(13),均采用具備TVD性質的3 階龍格庫塔格式進行時間積分.

1.3 壓力松弛

壓力松弛過程需要保證單元內部的質量、動量和能量守恒,通過調整各組分的體積分數、密度和內能使得各組分的壓強相等.具體的控制方程為

通過推導化簡[1]可以得到最終的松弛方程為

2 固體求解器

2.1 固體控制方程

更新拉格朗日格式下質量、動量和能量的守恒方程為

其中 ρ 表示當前時刻的密度,u表示位移,σ 表示當前構型下的柯西應力,b表示單位質量的體力,ε 表示應變張量,s表示單位質量的熱源,q表示單位時間下單位面積的熱流密度.

大部分的固體材料本構描述的是應變率張量和應力率張量之間的關系

其中 σ?為柯西應力張量的焦曼率導數.

為使式(16)定解,相應的動力學條件和初/邊值條件為

2.2 物質點法

如圖1 所示,物質點法將物質區域 Ω 離散為一組拉格朗日質點,因此連續體密度可以近似為

圖1 物質點法的空間離散Fig.1 Spatial discretization of MPM

其中np表示質點總數,mp表示質點p的質量,xp表示質點p的位置坐標,δ表示狄拉克函數.

物質點法在每個時間步中的拉格朗日步進行動量方程的求解,質點和網格保持相對位置固定,因此可以采用有限單元法中的形函數NI(x) 插值得到計算域的位移場

其中uI是網格節點的位移,ng是背景網格節點數量.

將式(20)和式(21)代入式(19)中,可以得到定義在背景網格節點上的離散動量方程

分別是網格節點I處的內力和外力,NIp=NI(xp)是質點p處的形函數,σp=σ(xp) 是質點p的應力,h是為了將弱形式(19)中最后一項面積分轉換為體積分而引入的假想邊界層厚度.

顯式物質點法采用中心差分法的蛙跳格式進行動量方程的時間步推進

2.3 物質點法格式

在顯式物質點法中,物質點上的應力狀態更新既可以在每個時間步開始時進行,也可以在每個時間步結束時進行,分別稱為USF(update-stressfirst)格式和USL(update-stress-last)格式[18].在改進的USL(modified USL,MUSL)格式中,將更新后的質點動量映射回網格節點,并采用該節點速度場進行質點應力狀態的更新.

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

3 流固耦合算法

3.1 網格單元分類

為了施加流固邊界條件,需要先定位浸沒邊界的位置.如圖2 所示,處于固體結構內部的網格單元被分類為type I,處于流場內部的網格單元被分類為type II,而被浸沒邊界穿過的網格單元被分類為type III.

圖2 網格單元分類Fig.2 Grid cell classification

固體結構的運動由拉格朗日屬性的物質點進行跟蹤,而圖2 中紅線所示的浸沒邊界僅用于說明,在實際模擬過程中并不會顯式重構邊界.因為當固體結構發生破碎后,通過三維空間內的點云跟蹤并重構新生的浸沒邊界過于復雜、不穩定并且耗時.可以從浸沒邊界的定義出發

其中 Ωp是質點p占據的區域,通過物質點的分布來確定網格單元的分類.具體步驟如下.

(1) 識別每個物質點位于哪個網格單元內,并初步將網格單元分為兩類:物質點單元和空單元.空單元即為流體網格單元(type II).

(2) 搜索每個物質點單元的相鄰單元.若所有相鄰單元均為物質點單元,則該物質點單元為固體網格單元(type I);否則,該物質點單元為流固兩相網格單元(type III).

3.2 基于拉格朗日乘子的連續力浸沒邊界法

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

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

根據式(23)和式(24)重構得到的背景網格節點動量和集中質量,固體結構在type III 單元格心處的速度場可以由下式給出

其中nen是每個單元的節點數量.通過將type III 單元格心處的流場速度vc設置為浸沒邊界的外插速度,浸沒邊界上的速度邊界條件就可以在每個時間步嚴格滿足.

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

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

根據動量守恒定律,式(31)歐拉源項的反作用力需要施加到固體結構上

另外,type III 單元還扮演著在物質點法中計算流固耦合表面應力的角色.通過將式(26) 中的替換為Fp,物質點法中的節點外力公式可以化簡為

4 數值驗證

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

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

為了驗證所建立的多介質有限體積黎曼求解器具備描述多介質流動及激波捕捉的能力,本小節模擬了經典的標準算例“二維激波與氦氣泡相互作用問題”[54].

激波前后的流場狀態應滿足激波間斷關系

其中MS為入射激波的馬赫數.激波前空氣區域的無量綱參數設置為

則激波后空氣區域的無量綱參數可以由式(35)和MS=1.22計算得到

氦氣泡區域相應的無量綱參數為

空氣和氦氣的力學行為均由理想氣體狀態方程

進行描述,相應的氣體常數分別為 γ1=1.4(空氣)和γ2=1.67(氦氣).

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

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

從圖5(a)和圖6 的對比中可以看到當前的多介質有限體積算法(MMFVM)可以很好地捕捉氦氣泡表面的激波演化,包括入射激波、折射激波、反射激波、側面激波和膨脹扇區.同時,從圖5 各個時刻計算模擬的紋影圖結果和實驗照片的對比中可以看出,MMFVM能夠很好地模擬多介質的流動現象并捕捉相應的物質界面形狀.

圖5 二維激波與氦氣泡相互作用問題的模擬結果(左)和實驗結果[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 近場爆炸下方形鋼筋混凝土靶板的失效模式

如圖7 所示,Wang等[55]針對近場爆炸下方形鋼筋混凝土靶板的失效模式進行了實驗研究.本小節將分別采用iMMFV-MPM 算法和MPM 算法對該問題進行計算模擬,從而證明iMMFV-MPM 算法相比于MPM 算法的優勢.

圖7 實驗裝置及鋼筋混凝土的布筋形式[55]Fig.7 Geometry setup of experiments[55]

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

圖8 計算模擬的幾何參數設置Fig.8 Geometry set for simulation

TNT 炸藥采用JWL 狀態方程進行描述[55]

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

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

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

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

鋼筋的力學行為采用簡化Johnson-Cook 模型進行描述[55]

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

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

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

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

iMMFV-MPM 模擬的鋼筋混凝土靶板中心點的時程位移曲線如圖11 所示,而MPM 的模擬結果由于靶板中心損傷過于嚴重而無法給出相應的中心點位移曲線.Wang等[55]在實驗中測量得到靶板中心點的最大位移為 3 5 mm,而iMMFV-MPM 的模擬結果為 3 5.9 mm,與實驗結果的相對誤差為 2.6%,吻合良好.

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

4.3 外爆載荷下建筑物的毀傷現象

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

圖12 實驗裝置及計算模擬的幾何建模Fig.12 Setup of experiment and simulation

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

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

圖13 流場壓強云圖與建筑物損傷云圖Fig.13 Pressure contour in fluid domain and damage contour in building

圖14 流場體積分數截面云圖與建筑物壓力云圖Fig.14 Volume of fraction contour on the plane of z=300mm and pressure contour in building

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

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

4.4 多腔室內爆炸試驗

郭志昆等[59]對扁平箱形密閉結構內爆炸問題進行了實驗研究,具體的幾何參數設置如圖16(a)所示.本文建立了如圖16(b)所示的幾何模型(未顯示頂板)進行計算模擬.

圖16 實驗和模擬的幾何參數設置Fig.16 Geometry setup of experiment and simualtion

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

主爆室外側墻中心點處的壓強時程曲線如圖17 所示,爆轟波作用在側墻中心點處的超壓峰值為0.938 MPa,實驗中相應測點處測得的超壓峰值為0.987 MPa,相對誤差為4.96%,與實驗結果吻合良好.幾個典型時刻中截面上的流場壓強云圖和體積分數云圖如圖18 所示.方形炸藥包爆炸后,爆炸產物急劇膨脹并推動周圍空氣產生爆轟波如圖18(a)所示,在t=1.5 ms 時,爆轟波與主爆室的側墻相互作用產生反射波,形成圖17 中的第1 個超壓峰值(peak 1).t=3.5 ms 時,頂板與底板的反射波在中截面處相交形成二次反射波如圖18(b)所示,該二次反射波與主爆室的側墻相互作用,形成圖17 中的第2 個超壓峰值(peak 2).相鄰側墻上的爆轟反射波傳播到外側墻中心點處形成第3 個超壓峰值(peak 3);頂板與底板的二次反射波在相鄰側墻上反射后傳播到外側墻中心點處形成第4 個超壓峰值(peak 4);正對側墻上的爆轟反射波傳播到外側墻中心點處形成第5 個超壓峰值(peak 5).

圖17 主爆室側墻中心點處的壓強曲線Fig.17 Pressure curve at the center of sidewall in explosion room

圖18 各時刻中截面上的流場壓強云圖和體積分數云圖Fig.18 Pressure contour and volume of fraction contour of middle plane at different instants

從圖18 各時刻的壓強云圖和壓強標尺中可以發現,主爆室的超壓峰值在1 MPa 左右,而側爆室的超壓峰值驟降為0.05 MPa～0.1 MPa,對角爆室的超壓峰值則為0.02 MPa 左右,這與實驗觀測到的各爆室的超壓峰值比例一致.從各時刻的體積分數云圖中可以發現,爆炸產物基本都殘留在主爆室內,主/側爆室之間的內墻開窗方向上形成對外的高壓射流,而在4 個角區方向上形成“蘑菇傘”狀的物質界面結構.

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

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

5 結論

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

(1) 相比于傳統的采用純MPM 模擬流固耦合問題,建立的iMMFV-MPM 方法在模擬炸藥爆炸產生爆轟波的過程中可以得到更銳利的波陣面,從而在后續爆轟波與固體結構相互作用以及固體結構毀傷現象的模擬中得到高精度的結果.

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

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