爆炸荷载作用下大型箱涵安全性的SPH-FEM耦合分析

崔 溦，卢珊珊，宋慧芳，张社荣

(天津大学水利工程仿真与安全国家重点实验室，天津 300072)

地下结构物的安全性是工程界关注的重要问题．以南水北调中线工程的大型输水箱涵为例，其单孔净尺寸达到 4.4,m×4.4,m，一旦受到损坏，将直接影响到沿线20多座大中城市的用水安全．

爆炸冲击荷载作用下，地下结构物的响应与地面结构物有很大不同，其最主要的差别在于必须考虑土与结构物的相互作用．虽然修正后的单自由度分析方法可以近似考虑这种相互作用，但在模型建立、参数选择以至计算结果上都存在较大复杂性和不确定性[1]．数值分析方法可以精确模拟爆炸冲击下地下结构物的响应，但现有的研究方法，无论是有限元、有限差分还是它们之间的耦合算法，一般将爆破过程与结构物响应割裂开来，直接将爆破荷载以压力的形式施加到结构物上，其计算精度难以得到保证[2-3]．精确模拟地下结构物爆破响应的主要难点在于数值模型必须包含药包周围的大变形土体与地下结构物的复杂体型，且要在计算精度与计算效率上都能满足工程要求[4]．

笔者以 SPH法模拟药包周围土体，以拉格朗日有限元法模拟复杂结构物，通过二者之间耦合研究大型输水箱涵在爆破荷载作用下的响应，以为工程设计与维护提供一定参考．

1 SPH与FE耦合算法

在采用常规有限元方法分析爆炸引起的大变形问题时，由于网格的剧烈扭曲经常会导致计算困难，为了克服这种困难，光滑粒子流体动力学(SPH)作为一种重要的无网格法得到了较为广泛的应用[5]．

SPH法把连续的物理量用多数粒子的集合来进行插值求解，连续体的运动用有限数量的粒子运动来离散化，由于没有使用网格，因此不会发生界面变形大所引起的计算溢出问题．计算空间导数时不需要使用任何网格，而是通过一个称为“核函数”的积分核进行“核函数估值”近似，将流体动力学基本方程组转换成数值计算的 SPH方程，整个流场离散成一系列“粒子”，所有力学量(密度、速度、压力、内能等)由这些粒子负载，这些粒子可以按计算公式，亦即按流体动力学运动的规律任意地运动．在 SPH法中，每个粒子 I与在和它相距给定距离范围内(一般假定为2,h，h为光滑长度)的其他所有粒子J发生相互作用．SPH法要引入一个特殊函数W(x-x′,h)作为核函数，场函数经过核函数“光滑化”再在整个求解域上积分，便得到了场函数的核估计．对场函数 f，其核估计表示为

式中x和 x′均为位置矢量．

图 1给出了核估计的基本概念．细节问题可参见文献[6]．

SPH法在模拟复杂体型的类似于墙的结构物时，由于相对于结构尺寸时墙体厚度较小，也就需要较小颗粒和时间步长，在计算效率与精度上，该方法受到一定限制．采用SPH与FE耦合技术，在材料小变形和结构复杂区域采用有限元模拟，在爆炸近域采用SPH模拟，可以克服上述缺点，并保证计算精度．

SPH颗粒和FEM网格有2种耦合方法，一种是通过可以相对滑动的接触面连接，另一种则是固接．在本次分析中，由于SPH颗粒和FEM网格的交界面位于土介质内，颗粒与网格之间不存在相对变形而固接在一起，因此根据运动方程，其他颗粒传递到交界面颗粒上的力与有限元网格传递到交界面颗粒上的力相同．图 2给出了 SPH粒子与传统拉格朗日有限元网格之间的固结耦合关系．

图1 SPH法的核估计Fig.1 SPH approximation

图2 SPH与FEM的耦合关系Fig.2 Coupling of SPH and FEM

2 材料模型

2.1 TNT

冲击波是爆破过程的重要特征．在本次研究中，炸药状态方程用JWL(Jones-Wilkins-Lee)方程来描述．JWL状态方程可表示为

式中：v为比容，1/vρ=；e为比内能；1C、1r、2C、2r和ω为常数，且对于常规爆破来说，它们的数值可以由动力试验来确定．本次数值分析所用 TNT的材料参数如表1所示．

表1 炸药的材料参数Tab.1 Material parameters for TNT

2.2 土

土用冲击状态方程和基于 D-P准则的弹塑性模型来描述，并定义其拉伸极限．

对于冲击荷载，即使在冲击速度约 2倍于初始声速0c和冲击压力达到 100,GPa量级的情况下，Rankine-Hugoniot状态方程都可以表示为

采用 D-P准则和拉伸极限来描述土的强度特性．为了避免土的剪胀问题，采用非关联的流动法则，其塑性势函数定义为

式中Y为屈服极限．

本次数值分析采用土体的具体参数如表2所示．

表2 土的材料参数Tab.2 Material parameters for soil

2.3 混凝土

冲击荷载下混凝土响应是一个复杂的非线性和率相关过程，在目前分析中，以 RHT模型应用最多，RHT模型可以分为强度模型和状态方程模型2部分[7]．对于压碎材料来说，强度模型主要包括弹性极限函数、破坏函数和残余强度函数3部分．

在硬化阶段之后，混凝土塑性应变导致的材料损伤和强度降低表示为

在该模型中，状态方程采用 P-α模型描述，其主要假定是在相同的压力和温度条件下，多孔介质材料的比内能和相同密度实心材料是一致的，该模型重点关注高应力下的混凝土性态，但同时也可以提供低应力状态下压缩过程的合理描述．混凝土模型的主要参数如表3～表5所示.

表3 混凝土的RHT模型物理参数Tab.3 Physical parameters for RHT model of concrete

表4 混凝土的RHT模型强度参数Tab.4 Strength ntensity parameters for RHT model of concrete

表5 混凝土的RHT模型状态参数Tab.5 State parameters for RHT model of concrete

2.4 钢 筋

在本次分析中，钢筋采用率相关的John-Cook弹塑性模型模拟．屈服应力表示为

在本次分析中，直接将钢筋和混凝土单元共结点连接，钢筋材料参数如表6所示．

表6 钢筋材料参数Tab.6 Parameters used for reinforcement steel bar

2.5 接触面模型

根据 Huck的试验成果，当土与结构物的接触面比较粗糙时，接触面应力超出土的最大剪应力时发生破坏，且高压力时接触面的强度特性与土体属性非常接近[8]．基于上述原因，本次分析将土与混凝土箱涵之间接触面按完全连接考虑．

2.6 边界条件

为满足辐射条件，数值模型中土为透射边界．透射边界条件允许应力波通过物理边界传输而没有反射[5]．假定边界的法向速度矢量为 un，边界处压力 p采用计算式如下：

式中：pref和 uref分别为相关压力和速度分量；I为材料阻抗．

3 数值模型

图3给出本次分析的SPH-FEM耦合模型，考虑无水工况，模型共包含钢筋混凝土箱涵、土体、炸药3种介质，并近似按二维轴对称平面应变问题分析．炸药及周边土体采用 SPH模拟，箱涵及周边土体采用FEM模拟，为了提高计算效率，近炸药及箱涵周边土体 SPH和 FEM 网格较密，距离较远则网格逐渐稀疏，FEM 网格大约 8万个，采用二维四结点平面单元，高斯积分法；SPH粒子大约2万个．

图3 SPH-FEM耦合模型Fig.3 SPH-FEM coupled model

图 4给出本次分析模型的具体尺寸及计算中监测点的位置．TNT质量为 80,kg，埋深为 d(1,m，2,m)，距箱涵最小水平距离为 h(6,m，8,m)．箱涵单孔净尺寸为 4.4,m×4.4,m，钢筋布置及其他尺寸如图4所示，钢筋保护层厚度为50,mm．

图4 模型布置Fig.4 Configuration of model

4 结果分析

4.1 爆坑形成

图 5给出不同埋药深度下土中爆坑形成过程．从图 5中可以看出，随着埋药深度的增加，爆坑形态发生明显变化，并逐渐由可见爆坑向爆腔过渡．数值模拟很好地揭示了爆破鼓包的运动过程和表面土体的层状剥离现象．在埋药深度1,m情况下，以抛掷爆破为主，形成爆坑最大直径约为 3.8,m．在埋药深度2,m情况下，以松动爆破形式为主．在地面形成比较明显的松动破坏区，主要为自由面的反射拉伸波所引起，最大空腔直径约为 4.2,m．根据Henrych[9]针对不同属性土、不同爆破方式下的爆腔经验公式计算得到爆腔直径分别为 3.67,m和4.18,m，与数值模拟结果吻合较好．

图5 不同埋药深度下爆坑形成过程Fig.5 Crater formation under different buried depths

4.2 土中爆炸冲击波传播

图 6给出了爆破后土中速度时程曲线．图中测点位置如图 4所示，其中测点 8、9位于 SPH网格内，测点 10、11位于 FEM 网格内．从图 6中可以看出，随着与药包距离增加，速度峰值逐渐减小；药包近域的压力变化具有冲击波的特征，单峰值且持续时间很短，随着与药包距离增加，速度变化开始具有应力波的特征，幅值减小且持续时间增加．通过与文献[9-10]对比可以看出，无论是测点 8、9还是测点 10、11，速度时程曲线特征与试验和理论计算结果一致，表明本文所采取的分析方法是较准确的．

图6 土中冲击速度时程Fig.6 Velocity time history in soil

4.3 混凝土损伤

图7给出不同爆破距离下箱涵损伤分布．从图7中可以看出，损伤严重区域主要出现在边墙与底板和顶板交界处，主要为冲击荷载引起的压损伤．爆炸距离对箱涵累积损伤有较大影响，随着水平距离h和埋深d的增加，累积损伤区域逐渐减小．

图7 不同爆破距离下混凝土损伤分布Fig.7 Distribution of damage in concrete at different detonation distances

4.4 箱涵响应

图8给出不同爆炸距离下箱涵变形情况．

图8 混凝土变形时程曲线Fig.8 Concrete deformation time history curves

从图 8中可以看出，不同测点变形时程差别明显．从测点1、2、3的水平向变形来看，不同爆炸距离下变形规律基本一致，边墙中间位置水平变形最大，最大量达到 12.5,mm，出现在最近爆炸距离(D＝1,m、h＝6,m)情况，残余变形量为不可回复的塑性变形，而边墙与顶板底板交界处变形较小．从测点 1、4、6、7的竖向变形来看，不同爆炸距离下变形规律存在一定差别，爆炸深度越大，箱涵整体将呈现一定向上移动趋势．

图9给出不同爆炸距离下箱涵应力变化情况．从图 9中可以看出，各测点应力时程差别明显．从测点

1、2、4的水平向应力来看，变化规律基本一致．在冲击峰值之后，存在较长时间的摆动现象，主要为结构自身响应引起．水平向应力最大值为-33,MPa，出现在最近爆炸距离(D＝1,m、h＝6,m)情况．边墙中间位置的应力时程与边墙和顶板交界处的应力时程明显不同，存在明显的二次加载现象，主要为箱涵周边土体流动变形导致压应力增加引起[11]．从测点 1、2、4的竖向应力来看，变化规律存在较大差别．测点 1存在较为明显的冲击峰值，2、4点则不明显．不同爆炸距离下测点2的应力时程差别较大，拉压应力的变化与埋深密切相关．

图9 混凝土应力时程曲线Fig.9 Concrete stress time histories curves

5 结 论

(1) 采用 SPH 法模拟爆破近域土体大变形，FE法模拟远域土体及箱涵结构，爆破后40,ms计算时间约为 14,h，具有较好的计算效率，计算模型较好地揭示了土中爆坑形成和冲击波传播过程，计算精度满足要求．

(2) 不同爆炸位置下箱涵结构爆破响应存在较大差异，尤其是边墙中间位置，受土体变形流动影响，存在较为明显二次加载现象，且爆破位置对其拉压变化影响明显．

［1］ Wdidlinger P，Hinman E. Analysis of underground protective structures[J].J Struct Eng，1987，114(7)：1658-1673.

［2］ Lu Yong，Wang Zhongqi，Chong Karen. A comparative study of buried structures in soil subjected to blast load using 2D and 3D numerical simulations[J].Soil Dynamics and Earthquake Engineering，2005，25(4)：275-288.

［3］ Wang Zhongqi，Lu Yong，Hao Hong. A full coupled numerical analysis approach for buried structures subjected to subsurface blast[J].Computersand Structures，2005，83(4/5)：339-356.

［4］ 杜义欣，刘晶波，伍 俊，等. 常规爆炸下地下结构的冲击震动环境[J]. 清华大学学报：自然科学版，2006，46(3)：322-326.

Du Yixin，Liu Jingbo，Wu Jun，et al. Blast shock and vibration of underground structures with conventional weapon[J].Journal of Tsinghua University：Science and Technology，2006，46(3)：322-326(in Chinese).

［5］ 王 吉，王肖钧，卞 梁. 光滑粒子法与有限元的耦合算法及其在冲击动力学中的应用[J]. 爆炸与冲击，2007，27(6)：522-528.

Wang Ji，Wang Xiaojun，Bian Liang. Linking of smoothed particle hydrodynamic method to standard finite element method and its application in impact dynamics [J].Explosion and Shock Waves，2007，27(6)：522-528(in Chinese).

［6］ Liu G R，Liu M B.Smoothed Particle Hydrodynamics：A Meshfree Particle Method[M]. London：World Scientific，2003.

［7］ 王 政，倪玉山，曹菊珍. 冲击载荷下混凝土动态力学性能研究进展[J]. 爆炸与冲击，2005，25(6)：519-524.

Wang Zheng，Ni Yushan，Cao Juzhen. Recent advances of dynamic mechanical behavior of concrete under impact loading[J].Explosion and Shock Waves，2005，25(6)：519-524(in Chinese).

［8］ Huck P J，Saxena S K. Response of soil-concrete interface at high pressure[C]//Proceedings of the Tenth International Conference on Soil Mechanics and Foundation Engineering. Stockholm，Sweden，1981：141-144.

［9］ Henrych J.The Dynamics of Explosion and Its Use[M].Amsterdam：Elsevier Scientific Publishing Company，1979.

［10］ Wang Z，Hao H，Lu Y. A three-phase soil model for simulating stress wave propagation due to blast loading[J].Int J Numer Anal Geomech，2004，28(1)：33-56.

［11］ James T，Baylot P E. Effect of soil flow changes on structures loads[J].J Struct Eng，2000，126(12)：1434-1441.