王超申，康道輝，游世輝，謝純凱，曽憲任

(1.湘潭大學(xué) 土木工程與力學(xué)學(xué)院,湘潭 411105; 2.棗莊學(xué)院 機(jī)電工程學(xué)院,棗莊 277160)

1 引 言

近年來，如何提高地下防護(hù)工程在防御戰(zhàn)時抵抗各種武器殺傷破壞的能力是國內(nèi)外學(xué)者一直研究的問題[1]。李硯召等[2]通過平面裝藥加載試驗,證明了分配層分層結(jié)構(gòu)對核爆炸防護(hù)效果優(yōu)于單層結(jié)構(gòu)。文獻(xiàn)[3，4]都以空殼復(fù)合材料為分配層進(jìn)行集團(tuán)裝藥化爆試驗，證明這類材料既對爆炸波有衰減彌散作用，也具備承受多次打擊的能力。趙躍堂等[5]模擬發(fā)現(xiàn)，分配層厚度越大，對沖擊波的削弱作用越明顯，并且支撐結(jié)構(gòu)頂部荷載越均勻。葉中豹[6]提出新型三維Kelvin-Maxwell雙粘性系數(shù)的等效本構(gòu)模型，并用于模擬分配層中的空殼復(fù)合材料，通過試驗和模擬結(jié)果對比，證明了這種本構(gòu)模型的可靠性。然而，以上研究的各層材料之間的界面是水平的，波的入射角是固定的，但根據(jù)波動理論，應(yīng)力波的幾何擴(kuò)散效應(yīng)會影響波的衰減彌散[7]，利用這種效應(yīng)，可改變分配層的幾何形狀來擴(kuò)散沖擊荷載。其次，分配層厚度在一定的情況下，各層材料的順序還有待優(yōu)化。

光滑粒子流體動力學(xué)法SPH(Smoothed Particle Hydrodynamics)可用于解決爆炸沖擊的問題，但存在精度與效率低和邊界難以處理的缺陷。而FEM法具有計算效率高和邊界條件設(shè)置方便的優(yōu)點(diǎn)，可以很好地彌補(bǔ)SPH法的不足[8]。SPH-FEM耦合法是在LS-DYNA的框架內(nèi)建立FEM單元和SPH粒子并耦合的數(shù)值分析方法，SPH法中不同材料粒子的相對運(yùn)動可形成界面滑移，光滑連續(xù)地處理物質(zhì)交界面，能很好地模擬爆炸沖擊引起的成坑和物質(zhì)飛散現(xiàn)象。

典型的成層式防護(hù)工程的結(jié)構(gòu)主要有偽裝土層、遮彈層、分配層和支撐結(jié)構(gòu)4部分。本文采取有限元網(wǎng)格與光滑粒子相耦合的方法，重點(diǎn)研究不同排列方式的馬蹄形分配層結(jié)構(gòu)對爆炸波衰減彌散的影響，旨在尋找更優(yōu)良的分配層結(jié)構(gòu)來提高地下防護(hù)工程的抗爆防震性能。

2 SPH-FEM耦合模擬方法和模型建立

2.1 SPH-FEM耦合數(shù)值模擬方法

SPH法的核心是核函數(shù)近似法。對于任意函數(shù)f(x),都可以采用積分形式

(1)

式中Ω是包含x積分域，δ是狄拉克函數(shù)。采用光滑函數(shù)W(x-x′,h)取代狄拉克δ(x-x′)函數(shù)，f(x)的積分表達(dá)式為

(2)

式中W(x-x′,h)稱為光滑函數(shù)或者核函數(shù)kernel function，h定義為核函數(shù)的光滑長度。用粒子體積ΔV對dx′進(jìn)行替換，粒子的質(zhì)量mj為

(3)

粒子在i處的表達(dá)式為

(4)

SPH-FEM耦合算法多采用固定耦合算法，設(shè)置點(diǎn)面固連接觸，根據(jù)線動量守恒和角動量守恒調(diào)整粒子與單元節(jié)點(diǎn)的速度和位置，保證兩者間的位移協(xié)調(diào)，并且將粒子區(qū)的力學(xué)信息傳遞到有限元網(wǎng)格。耦合界面處的粒子與單元節(jié)點(diǎn)一一對應(yīng)，粒子與其對應(yīng)的單元節(jié)點(diǎn)具有相同的加速度，其大小是由粒子與單元節(jié)點(diǎn)的質(zhì)量以及作用于兩者的力共同確定[9,10]。在本文的模擬中，近場SPH區(qū)域和遠(yuǎn)場FEM區(qū)域是相同的巖石介質(zhì)，粒子與其對應(yīng)的單元節(jié)點(diǎn)將在整個計算過程中保持完全一致，如圖1所示。

2.2 材料模型

本文采用的減震吸能的材料有砂土和泡沫混凝土。砂土是含水細(xì)砂，含水率控制在3%左右。泡沫混凝土模型使用可壓扁材料泡沫模型*MAT_CRUSHABLE_FOAM。泡沫混凝土相關(guān)參數(shù)如下[11]，密度為0.799 g/cm3，泊松比為0.1，楊氏模量為342.2 MPa，抗拉強(qiáng)度為0.2 MPa，阻尼系數(shù)為 0.1。巖石材料和混凝土材料均采用適應(yīng)高壓、高應(yīng)變率的模型*MAT_JOHNSON_HOLMQUIST_CONCRETE(HJC),HJC材料模型在爆炸荷載下能夠很好地反映巖石等脆性材料損傷效應(yīng)的動態(tài)響應(yīng)，其等效屈服強(qiáng)度是壓力、應(yīng)變率和損傷的函數(shù)，巖石材料[12]、鋼筋混凝土[13]和混凝土[14]的相關(guān)參數(shù)列入表1。本文使用Wang[15]推薦的細(xì)砂，其本構(gòu)模型采用*SOIL_AND_FOAM，主要的參數(shù)值如下，密度為1800 kg/m3，剪切模量為 63.8 MPa，初始體積模量為134.0 MPa,內(nèi)摩擦角為35°，表觀粘聚力近似取 5 kPa，含水率為3.4%。另外，為增加分配層的抗爆性和增大材料間的不匹配度，在分配層中加入鋼鐵。

圖1 SPH與FEM的耦合

Fig.1 Coupling of SPH to FEM

2.3 幾何模型

在本文數(shù)值模擬中，考慮到巖體觸地爆炸的問題具有對稱特點(diǎn)，故采用1/4模型。采用SPH-FEM耦合法建模，先在LS-DYNA中建立TNT炸藥-巖體-錨桿-襯砌的有限元模型，數(shù)值模型尺寸為470 cm×400 cm×840 cm，地下硐室采用馬蹄形結(jié)構(gòu)，為典型的成層式防護(hù)工程。地下巷道的最外層是100 mm的噴射混凝土層，向內(nèi)是 300 mm 的分配層，再將分配層分為6層，如圖2所示，這 6層 分配層每層的厚度均為50 mm。耦合模型中鋼筋材料采用BEAM161四節(jié)點(diǎn)梁單元，其余材料均為SOLID164六面體單元。整個模型共劃分為496960個有限單元和30762個SPH粒子，其中炸藥粒子為216個，爆炸近區(qū)巖體粒子為30546個。圖3為Lagrange單元和SPH粒子的分布情況。

SPH粒子與有限單元面采用自動點(diǎn)面接觸，對稱面施加全局約束以及SPH對稱邊界，為描述半無限空間，其余面定義無反射邊界條件，并設(shè)置吸收剪切波和膨脹波，防止爆炸波在邊界處發(fā)生反射，并對時間步長和人工體積粘性系數(shù)作相應(yīng)調(diào)整。為了更好地觀測鋼筋混凝土破壞情況，設(shè)置了*MAT_ADD_EROSION失效準(zhǔn)則,以最大拉應(yīng)力為破壞準(zhǔn)則。

圖2 分配層結(jié)構(gòu)

Fig.2 Distribution layer structure

圖3 Lagrange單元SPH粒子分布情況

Fig.3 Lagrange unit and SPH particle distribution

3 模擬結(jié)果對比分析

對分配層結(jié)構(gòu)的Lagrange單元進(jìn)行監(jiān)測，研究11個方案下各單元的位移、應(yīng)力和應(yīng)變。根據(jù)應(yīng)力波傳播的特點(diǎn)，距爆源中心越近，應(yīng)力波峰值越大，因此選取模型對稱軸線上的 7個 測點(diǎn)(A= H371761，B= H371600，C= H437961，D= H349280，E= H394561，F(xiàn)= H21440，G= H511280)作為研究對象，每個測點(diǎn)的單元在厚度上長度為 1.7 cm，表2是各測點(diǎn)應(yīng)力峰值。測點(diǎn)分布如圖4所示。

3.1 含不同材料的分配層結(jié)構(gòu)對比分析

當(dāng)分配層采用砂土?xí)r，等效應(yīng)力呈現(xiàn)折線增大的趨勢，圖5(a)中，應(yīng)力波在大約750 μs時開始從0增加，1500 μs時幅值開始降低，2250 μs時才達(dá)到峰值。而圖5(b)應(yīng)力波大約在600 μs時開始從0增加，1750 μs時幅值開始降低，3500 μs時到達(dá)峰值。從圖6可以看出，5000 μs時，鋼筋混凝土支撐結(jié)構(gòu)出現(xiàn)失效單元(淺灰色單元)，表明砂土能更大程度地削弱到達(dá)鋼筋混凝土支撐結(jié)構(gòu)的應(yīng)力波強(qiáng)度。從表2中測點(diǎn)F的應(yīng)力峰值可知，方案〈1〉和〈2〉的爆炸波應(yīng)力峰值分別衰減了99.2%和 96.4%。由于鋼鐵與砂土的波阻抗比值要大于鋼鐵與泡沫混凝土的波阻抗比值，分配層厚度一定的情況下，含水砂土對爆炸應(yīng)力波的衰減明顯好于泡沫混凝土。

表2 測點(diǎn)應(yīng)力峰值

Tab.2 Stress peak at the measuring point

方案應(yīng)力峰值/MPa遮彈層ABCDEFG〈1〉ggssss180.0180.00180.001.501.351.611.4220.58〈2〉gghhhh180.0180.00180.006.886.066.156.4722.57〈3〉gghhss180.0180.00180.006.845.990.890.7816.70〈4〉gghshs180.0180.00180.006.966.926.861.0120.89〈5〉ggshsh180.0180.00180.001.196.571.546.8120.12〈6〉gshgsh180.0180.002.116.1315.361.304.376.03〈7〉gshshg180.0180.001.386.827.011.1415.5322.12〈8〉ghsghs180.0180.006.971.8023.526.861.0816.68〈9〉shshsh180.02.166.900.776.521.176.128.46〈10〉hshshs180.06.855.356.725.245.700.787.25〈11〉shgshg180.02.216.9621.620.443.5620.194.77注:將鋼鐵材料縮寫為g,泡沫混凝土為h,砂土為s。

圖4 測點(diǎn)分布

Fig.4 Survey point distribution map

3.2 不同排列方式分配層結(jié)構(gòu)的對比分析

爆炸波在不同介質(zhì)中能量損耗和衰減速度相差很大，分配層中相鄰介質(zhì)材料的波阻抗、應(yīng)變率影響因子黏度系數(shù)和材料的屈服強(qiáng)度等相差越大，對爆炸波的吸能彌散作用越明顯[16]，并且爆炸波從硬材料入射到軟材料中對波的削弱效果更強(qiáng)。另外，分配層材料采用周期性布置時，彌散效應(yīng)會得到增強(qiáng)，本文的數(shù)值模擬也采取這種布置方式。

3.2.1 測點(diǎn)應(yīng)力對比分析

圖7是方案〈6〉～〈11〉測點(diǎn)G的等效應(yīng)力曲線?？梢姡瑳_擊波大約在800 μs時到達(dá)鋼筋混凝土支撐結(jié)構(gòu)，隨著時間的推移，沖擊波在分配層結(jié)構(gòu)內(nèi)不斷發(fā)生反射和透射，使得鋼筋混凝土支撐結(jié)構(gòu)的應(yīng)力一直在增大。當(dāng)分層材料相同并且各種材料厚度相同時，方案〈9〉和〈10〉都只將砂土和泡沫混凝土采用周期性布置的方式，方案〈9〉中每個周期內(nèi)材料波阻抗值是遞增的，而方案〈10〉每個周期內(nèi)材料波阻抗值是遞減的，應(yīng)力波峰值分別衰減了95.4%和96.0%。

圖5 砂土和泡沫混凝土測點(diǎn)F有效應(yīng)力時程曲線

Fig.5F-point effective stress time-course curve of sand and foam concrete

方案〈9〉和方案〈11〉的泡沫混凝土和砂土排列順序一致，方案〈9〉中無鋼鐵，其應(yīng)力波峰值衰減了95.3%；方案〈11〉分配層從上往下，第三層和第六層為鋼鐵，應(yīng)力波衰減了97.4%，削波能力較好。分配層中加入鋼鐵后，更大程度削弱了應(yīng)力波的傳播，這是由于鋼鐵受反射波和大變形的作用，使得其也隨多孔材料一起變形，吸收了一部分能量，因此局部的破壞作用加大，而鋼筋混凝土支撐結(jié)構(gòu)受到的破壞減小。

圖6 砂土與泡沫混凝土分配層等效應(yīng)力云圖

Fig.6 Equivalent stress cloud diagram of sand and foam concrete distribution layer

圖7 方案〈6〉～〈11〉測點(diǎn)G的等效應(yīng)力曲線

Fig.7 Scheme 〈6〉～〈11〉 equivalent stress curve of pointG

方案〈6〉和方案〈8〉分配層從上往下，第一層和第三層為鋼鐵，其他四層分別按波阻抗值遞增和遞減來排列，應(yīng)力波峰值分別衰減了96.7%和90.9%。分層結(jié)構(gòu)中材料按波阻抗比遞增的順序排列時，會增強(qiáng)應(yīng)力波強(qiáng)度，因此，本節(jié)模擬只考慮按波阻抗比順序排列減小的情況。方案〈8〉中一個周期的材料順序是按波阻抗值遞減和波阻抗比減小的方式來布置的，應(yīng)力波衰減了 90.9%；方案〈11〉一個周期的材料順序是采用相鄰界面材料波阻抗值遞增和波阻抗比減小的方式來布置，應(yīng)力波衰減97.4%,略優(yōu)于方案〈8〉。很明顯，分配結(jié)構(gòu)中各層材料順序按波阻抗值遞增和波阻抗比減小、周期性布置和加入鋼鐵是較優(yōu)異的方案。

3.2.2 能量對比分析

圖8是六種方案中鋼筋混凝土支撐結(jié)構(gòu)的總能量變化曲線。可以看出，方案〈6〉和〈9〉～〈11〉在5000 μs時鋼筋混凝土支撐結(jié)構(gòu)的總能量大約是方案〈7〉和〈8〉的1/10。其中，方案〈11〉的支撐結(jié)構(gòu)總能量最低，爆炸開始后，動能上升速度加快，內(nèi)能上升速度較緩，5000 μs時，支撐結(jié)構(gòu)動能占總能量的75%；而方案〈7〉的支撐結(jié)構(gòu)總能量最高，5000 μs時，其動能占總能量的60%。對比兩種方案可知，動能是總能量最重要的組成部分，如圖9所示。方案〈7〉的支撐結(jié)構(gòu)動能為39636 J，內(nèi)能為28455 J；方案〈11〉的支撐結(jié)構(gòu)動能為5578 J，較方案〈7〉約低86%，內(nèi)能為1838 J，較方案〈7〉約低94%。方案〈8〉每一周期的材料是按波阻抗遞減的方式來布置的，總能量約為70450 J，對比分析可知，此方案吸能效果較差。從圖8能量曲線的斜率來看，5000 μs后，方案〈7〉和〈8〉中鋼筋混凝土結(jié)構(gòu)的總能量增加的速率要比方案〈11〉大，而方案〈11〉能量的增長明顯變緩。

圖8 方案〈6〉～〈11〉支撐結(jié)構(gòu)總能量曲線

Fig.8 Scheme 〈6〉～〈11〉 total energy curve of the support structure

3.3 拱頂位移對比分析

表3是11種方案測點(diǎn)G的Y方向位移。在爆炸載荷作用下，鋼筋混凝土支撐結(jié)構(gòu)離爆心最近的部位豎向位移會最大，巷道拱頂?shù)呢Q向位移與支撐結(jié)構(gòu)頂部的位移很接近，因此，可以通過控制點(diǎn)G的豎向位移來減少爆炸波對支撐結(jié)構(gòu)的破壞。11種方案中，除了方案〈7〉和〈11〉位移峰值偏大，方案〈6〉和〈8〉測點(diǎn)G的位移峰值最小以外，其他方案的位移峰值比較接近，約為最大位移的60%，可以說明這兩個方案中的支撐結(jié)構(gòu)受到來自分配層的能量較小，而方案〈6〉對支撐結(jié)構(gòu)的保護(hù)更佳。分配層的變形不斷調(diào)低應(yīng)力峰值，減緩并減弱了爆炸波對支撐結(jié)構(gòu)的破壞，對于遭受爆炸荷載的防護(hù)工程來講，巷道拱頂?shù)膿隙戎?沉降量)是判定防護(hù)結(jié)構(gòu)能否抵抗二次爆炸沖擊的關(guān)鍵指標(biāo)之一。因此，提高鋼筋混凝土支撐結(jié)構(gòu)的抗壓彎破壞能力的同時,還需科學(xué)地優(yōu)化分配層結(jié)構(gòu)。

圖9 方案〈7〉和方案〈11〉支撐結(jié)構(gòu)動能和內(nèi)能曲線

Fig.9 Scheme 〈7〉 and 〈11〉 kinetic energy and internal energy curve of support structure

表3 測點(diǎn)G豎向位移峰值

Tab.3 Peak displacement ofGpoint vertical displacement

方案位移峰值/cmABCDEFG〈1〉ggssss2.212.191.220.840.580.210.03〈2〉gghhhh1.941.911.611.040.530.440.15〈3〉gghhss1.951.921.631.071.200.450.05〈4〉gghshs1.931.911.610.680.620.440.14〈5〉ggshsh1.911.901.191.150.560.390.14〈6〉gshgsh2.681.400.780.570.230.120.04〈7〉gshshg3.511.971.710.580.060.060.06〈8〉ghsghs2.441.690.600.690.480.340.10〈9〉shshsh2.191.390.820.470.150.100.04〈10〉hshshs1.460.870.590.230.160.120.04〈11〉shgshg1.900.860.180.060.040.030.03

3.4 誤差分析

在強(qiáng)烈爆炸作用下，爆炸近區(qū)的動態(tài)響應(yīng)非常復(fù)雜，現(xiàn)有的連續(xù)介質(zhì)力學(xué)和線彈性力學(xué)等理論基礎(chǔ)不足以準(zhǔn)確描述爆炸近區(qū)的巖體變形特點(diǎn)，我國對于常規(guī)武器在掩體的爆炸效應(yīng)的研究大多是在前蘇聯(lián)的試驗成果的基礎(chǔ)上進(jìn)行的。如謝苗基(Щемякин Е И)和麥德維捷夫(Медведев Н С),通過大量試驗得出沖擊波壓力和距離的關(guān)系[17]。

δm=p2/γα

(6)

(7)

測點(diǎn)依然是從11種方案的地下硐室中心線上選擇，應(yīng)力峰值誤差分析選擇距爆心 1.30 m,1.45 m,1.55 m,1.65 m,1.75 m,1.85 m和1.95 m處的測點(diǎn)。根據(jù)模擬的數(shù)據(jù)顯示，距爆心1.80 m處，應(yīng)力峰值變小的速率加快，1.95 m時急劇變小，與經(jīng)驗值的誤差變大,列入表5。這是由于材料具有一定的粘性，隨著應(yīng)力波傳播距離的增加，不斷消耗能量的同時還隨著波面面積的增加造成能量的擴(kuò)散而衰減，沖擊波逐漸變?yōu)閴嚎s波。經(jīng)計算及對比分析，在爆炸近區(qū)，測點(diǎn)的數(shù)值模擬結(jié)果與計算值之間誤差在2.1%～3.5%之間，考慮到參數(shù)設(shè)置以及網(wǎng)格劃分不是很準(zhǔn)確，兩者計算結(jié)果較為吻合。

表4 經(jīng)驗公式參數(shù)

Tab.4 Empirical formula parameters

裝藥初始密度/g·cm-3爆速/m·s-1巖石初始密度/g·cm-3聲傳播速度/m·s-1應(yīng)力衰減指數(shù)1.60175302.83223501.95

表5 數(shù)值模擬結(jié)果與經(jīng)驗公式對比

Tab.5 Comparison of numerical simulation results with empirical formula

項目沖擊波壓力/MPa1.30m1.45m1.55m1.65m1.75m1.85m1.95m數(shù)值模擬161.53129.60107.66100.8885.5066.6634.56經(jīng)驗公式157.07126.94111.4798.6787.9878.9471.24誤差/%2.82.13.52.22.915.651.5

4 結(jié) 論

(1)從含不同材料的分配層結(jié)構(gòu)對比分析結(jié)果發(fā)現(xiàn)，低含水砂土對波的衰減和彌散作用要優(yōu)于泡沫混凝土，可見砂土的波阻抗值較小。

(2)采用有利于削弱應(yīng)力波峰值的周期性布置分配層材料方式模擬。從方案〈9〉和〈10〉模擬的結(jié)果來看，每一周期材料按波阻抗值增加順序排列要優(yōu)于其按波阻抗值減小削波能力。而對比方案〈6〉和〈8〉，泡沫材料和砂土組合順序不變，加入鋼鐵后，每一周期材料按波阻抗值遞增順序排列要更優(yōu)于其按波阻抗值遞減的削波能力?？梢姡峙鋵咏Y(jié)構(gòu)的軟材料中加入硬材料后，更能削弱應(yīng)力波的強(qiáng)度。

(3)大波阻抗的鋼鐵與低波阻抗值的砂土緊密貼合極大增強(qiáng)了分配層的削波能力，沒加鋼鐵與加鋼鐵的方案即〈9〉和〈11〉對比，方案〈11〉的鋼筋混凝土支撐結(jié)構(gòu)G處的應(yīng)力峰值比方案〈9〉約小44%，表明加入鋼鐵的分配層吸能效果更佳。與方案〈8〉相比，方案〈11〉采用周期材料按相鄰界面材料波阻抗值遞增和波阻抗比減小的方式來布置更優(yōu)。

(4)通過對11種方案鋼筋混凝土支撐結(jié)構(gòu)的總能量和豎向位移進(jìn)行對比分析，方案〈11〉中支撐結(jié)構(gòu)的總能量和位移響應(yīng)最小，并且方案〈11〉的分配層結(jié)構(gòu)主要是削弱了應(yīng)力波的動能。因此，方案〈11〉的分配層結(jié)構(gòu)抗震吸能效果最佳，對支撐結(jié)構(gòu)的保護(hù)作用最大。

(5)將數(shù)值模擬得到的不同測點(diǎn)的應(yīng)力峰值與前蘇聯(lián)經(jīng)驗公式的計算值對比發(fā)現(xiàn)，在距爆心較近位置的應(yīng)力峰值計算結(jié)果和模擬結(jié)果比較吻合。