混凝土爆破损伤的SPH-FEM耦合法数值模拟*

王志亮,毕程程,李鸿儒

(合肥工业大学土木与水利工程学院，安徽 合肥 230009)

爆炸过程涉及炸药爆炸能量的高速释放、冲击波在介质中的快速传播、混凝土的损伤破裂乃至飞溅等方面的问题，目前的理论研究水平在准确描述爆炸过程方面尚存在不足。实验研究由于成本高、实施周期长以及很难全面地测定相关参量，最终往往只能得到混凝土在爆炸破坏作用下的终态响应结果以及获取部分监测点的时程信息,而数值模拟可以很好地再现爆炸响应的中间过程，并且能有效地描述混凝土在爆炸作用下的真实动态响应情况。相比于传统的FEM数值算法，如Lagrange、Euler以及ALE算法，光滑粒子流体动力学(smoothed particle hydrodynamics, SPH)方法[1]集合了无网格法、拉格朗日法和粒子流法的优势，能够避免Lagrange算法中的网格畸变问题、Euler算法中的重分和输运计算问题以及ALE算法兼具双重网格属性带来的计算耗时较长等问题。SPH法所具有的无网格属性可以避免网格扭曲变形造成的精度损失，拉格朗日属性可以得到较为准确的模拟结果，粒子属性能够较好地展现爆炸冲击所产生的成坑、物质飞溅等破坏现象。但是由于算法差异，SPH法相对于传统的FEM数值算法计算效率较低，且在边界处理上也存在着一定的问题。

针对SPH法存在的问题，诸多学者提出了SPH-FEM耦合法，即建模时将爆炸近区大变形区域采用SPH建模，而在爆炸远区的小变形区域采用FEM建模。这样既能避免大变形区域内的网格畸变，又可以减少SPH的计算域，从而大大提高计算效率。胡英国等[2]采用SPH-FEM耦合法对深孔梯段爆破的动力效应进行了数值模拟，模型在爆破近区采用SPH粒子，远区采用有限单元，得出该耦合算法能够较为准确地模拟爆破近区的岩体运动和损伤特征；王维国等[3]模拟了炸药在砂土地基中的触地爆炸，分别比较了纯SPH模型和SPH-FEM耦合模型对计算结果的影响，得出纯SPH模型计算耗费的时间约是SPH-FEM耦合模型的3倍，而计算精度基本保持一致，模拟得到的爆炸成坑效果与经验值相一致；崔溦等[4]采用SPH-FEM耦合数值法，研究了土中钢筋混凝土箱涵的爆炸动力响应问题，其将SPH法用于模拟爆破近域土体，而FEM法用于模拟远场土体和钢筋混凝土箱涵；Lu等[5]研究地下结构在浅埋炸药爆炸作用下的动态响应时，分别针对二维和三维的SPH-FEM模型进行数值比较，得出二维模型计算的结果更为精确且波形稳定；Koneshwaran等[6]在分析地表炸药爆炸对地下隧道的破坏效应时，对炸药及其附近的土体采用SPH粒子，并讨论了常见的几种SPH-FEM耦合方法；VUYST等[7]结合平板对撞的例子，讨论了SPH法、FEM法以及SPH-FEM耦合算法在计算效率和精度上的差异，得出SPH-FEM耦合算法效果最佳，并在刚形体跌落入水以及子弹击穿靶板例子上得到良好验证；杨刚等[8]模拟炸药在混凝中爆炸时，分别采用纯Euler算法、纯SPH算法、Euler-Lagrange耦合算法以及SPH-Lagrange耦合算法，指出用纯Euler算法或纯SPH算法描述混凝土时，计算得到的压力时程曲线在衰减段会有较大的数值波动。相比较而言，Euler-Lagrange耦合算法以及SPH-Lagrange耦合算法模拟效果较好。

Holmquist-Johnson-Cook(HJC)本构模型[9]是针对混凝土材料提出的一种率相关的损伤本构模型，能有效地计算混凝土在高应变率下的大变形问题。由于该模型能够较好地描述混凝土在高速碰撞、侵彻以及爆炸作用下的力学行为，适用于拉格朗日和欧拉算法，现已被LS-DYNA程序引入,并被广泛运用[10-15]，纪冲等[16]、梁超[17]还将SPH-FEM耦合法结合HJC本构模型运用到弹体侵彻混凝土上，而目前基于SPH-FEM耦合法模拟混凝土爆破的公开发表的成果很少，结合爆破问题对HJC本构中各参数敏感度进行分析的研究更少。本文中针对C30混凝土中的爆破成坑等问题，开展较深入的研究，力图得出一些具有参考价值的结论。

1 算 法

1.1 SPH基本原理

SPH方程可表述[1]如下。

第一步：使用积分表达式对函数进行近似，任一函数f(x)的积分近似表达式为：

(1)

式中：积分域为整个问题域空间,W(x-x′,h)为光滑核函数,h为光滑长度。

第二步：粒子近似法，即用粒子函数对式(1)的积分表达式离散化，通过对离散化的粒子插值得到函数值，也就是使用光滑核函数W影响范围内的临近粒子j的运动信息求和平均代替参考粒子i的运动状态，则粒子i处的函数近似式为：

(2)

式中:mj、ρj分别为粒子j的质量和密度(j=1, 2, …,N),N为在粒子i的支持域内粒子数的总量。

对粒子i的求解可通过使用光滑核函数W支持域内的粒子j求和得到(见图1)。图中h为光滑长度，W为光滑核函数，κ是与x处光滑函数相关的常数，Ω表示支持域，S为粒子作用边界。

1.2 SPH-FEM耦合方式

SPH-FEM耦合法的关键是在最靠近有限元网格的一排 SPH 粒子如何顺利地将各物理力学信息传递给有限单元。目前比较常见的有三种方法[6]：第一种是将靠近SPH 粒子的有限单元上的节点赋予粒子属性，SPH 粒子将各信息传递给有限单元上的节点，再通过节点传递给单元；第二种是在SPH粒子与有限单元之间建立混合区，该区域粒子与单元相互重合，将粒子携带的信息过渡到单元上；第三种是通过定义节点与面之间的接触来传递信息(见图2)，常见的包括滑动接触和固连接触[3]。

本文中采用固连接触的方式将SPH粒子与有限单元联系起来，最靠近有限元网格的SPH粒子将从其他粒子传递过来的能量、动量以及运动方程等物理力学信息传递给有限单元，同时以点-面固结的方式保证两者间位移与变形的协调。接触通过在LS-DYNA中建立模型后生成的k文件中添加关键字 CONTACT_TIED_NODES_TO_SURFACE来定义[18]，并通过关键字CONTROL_CONTACT调节接触刚度等参量。

2 本构模型

2.1 混凝土材料模型

HJC损伤本构模型考虑了高静水压力、高应变率和材料的损伤效应，主要包括三个部分：强度模型、状态方程和损伤模型[9,18]。

强度模型是以特征化等效应力进行描述，可表示为：

(3)

损伤模型如图4所示，由等效塑性应变和体积应变累积而成，其演化方程为：

(4)

D1(P*+T*)D2≥EFmin

(5)

本文中涉及到的C30混凝土[19-20]基本力学参数为：密度ρ0=2 400 kg/m3，抗压强度fc=39.2 MPa，抗拉强度T=3.162 MPa，弹性模量E=33.4 GPa，泊松比ν=0.202。混凝土剪切模量G和体积模量K根据下式确定：

G=E/[2(1+ν)]

(6)

K=E/[3(1-2ν)]

(7)

求得G=13.89 GPa，K=18.68 GPa。

混凝土达到弹性极限时的静水压力Pcrush和体积应变μcrush可由下式确定：

Pcrush=fc/3

(8)

μcrush=Pcrush/K

(9)

求得Pcrush=13.07 MPa，μcrush=0.000 7。

损伤参数D1可按照公式D1=0.01/(1/6+T*)确定，其中T*=T/fc，求得D1=0.04。参数D2和Smax对模拟结果的影响很小[10,14]，因此D2按文献[9]取为1，Smax取为7。

文献[19]针对C30混凝土在考虑率型微损伤演化的条件下改进了Johnson-Cook模型，并基于相关实验给出了极限面参数A、B、N及应变率影响参数C的值，所提出的模型在混凝土动态实验中得到了良好的应用。本文中基本力学参数大多取自文献[19]，考虑到HJC各参数之间耦合紧密，故参数A、B、N和C仍取自文献[19]。状态方程参数Plock、μlock、K1、K2和K3在缺乏实验数据的情况下，大多取自文献[9]。根据文献[11]对FS的假设，认为混凝土达到极限密度ρmax时，体积应变达到最大值μmax，即：

(10)

表1 C30混凝土HJC参数Table 1 HJC parameters of C30 concrete

2.2 炸药材料模型

炸药按均匀连续介质考虑，采用MAT_HIGH_EXPLOSIVE_BURN材料模型和JWL状态方程精确描述爆炸过程中爆炸产物的体积、压力以及能量特性，表达式为[18]：

(11)

式中：P为爆轰压力，V为相对体积，E0表示初始体积内能，A、B、R1、R2和ω为炸药常数。

炸药各参数取值分别为[21]：ρ=1 630 kg/m3，Dv=6 930 m/s，Pcj=27.0 GPa，A=374 GPa，B=7.33 GPa，R1=4.15，R2=0.95，ω=0.3，E0=7.0 GJ/m3，V0=1.0。

2.3 土体材料模型

土体采用MAT_SOIL_AND_FOAM材料模型，其理想塑性屈服函数为[18]：

φ=J2-(a0+a1σ+a2σ2)

(12)

式中：J2=σijσij/2，σij为偏应力分量；σ为平均应力；常数α0、α1和α2为无量纲曲线J2-σ二次拟合曲线常数项。

土体主要计算参数为[22]：密度为1.80 g/cm3，剪切模量为16.01 MPa，体积卸载模量为1.26 GPa，剪切屈服面参数α0=2.4×106Pa2，α1=1.360×104Pa和α2=0.123 2。

3 模拟分析

3.1 模型建立

穆朝民等[20]、李重情等[23]为研究变埋深条件下混凝土中爆炸压力和加速度的传播规律，进行了多组实验。实验设计在室外空旷地区，在地上开挖一个边长为3 m的立方体爆坑，沿深度方向布设6个测点，用于安装压力传感器和加速度传感器，然后现场浇筑混凝土直至填满实验坑，并在浇筑时预留装药孔，填装1 kg TNT，最后进行爆炸实验。本文以此实验为基础，对其进行数值模拟。文献[20,23]给出了炸药在不同埋深条件下包括装药比例距离为-0.25 m/kg1/3的空爆(炸药悬置空中距混凝土表面0.25 m)、-0.053 m/kg1/3的触地爆(炸药置于混凝土表面)、0 m/kg1/3的半埋爆(炸药半埋于混凝土中，且炸药中心处于混凝土顶面中心处)以及0.8 m/kg1/3的全埋爆(炸药全部埋于混凝土内且距混凝土表面0.8 m)等爆炸实验结果，而其中只给出了装药比例埋深为0和0.8 m/kg1/3时相应的压力和加速度时程数据。因此，为方便将模拟结果与实验结果进行对比，同时更好地展现爆炸破坏的效果，本文中选择模拟装药比例埋深为0 m/kg1/3的半埋爆情况。

物理模型如图5所示，立方体混凝土边长为300 cm，炸药选用1 kg TNT，按照等体积原则将炸药简化成尺寸为8.5 cm×8.5 cm×8.5 cm的立方体。混凝土的四周及底部与土体接触，土体厚度均取10 cm。测点1～6距炸药中心的距离依次为25.5、51.5、81.0、125.6、187.9和250.4 cm。

计算模型如图6所示，由于模型具有对称性，为节省计算时间只建了1/4模型。考虑到混凝土采用的是HJC本构模型，该模型无法模拟混凝土在爆炸荷载作用下的裂纹扩展、破碎等过程，因此通常情况下需要添加单元失效准则MAT_ADD_EROSION，使处于破坏区的单元失效以形成破坏效果。然而，单元失效准则对粒子并不起作用，粒子只会产生运动而不会凭空消失，此时如果混凝土粒子区域选择得过小，爆炸荷载传递到混凝土近区的有限单元上导致该处的单元删除，则粒子与单元之间的接触算法会失效，影响压力波的传播和单元与粒子之间位移变形的协调性。如果不添加单元失效准则，爆炸近区混凝土有限单元会发生变形，形成“杯口状”的凸起，跟实际的爆破形态并不一致；若将粒子区域取得过大，这样不仅会造成计算时间的增加，对应的存储空间也会急剧增加，一般的计算机很难满足要求。经过多次尝试，最终将该部分混凝土(近区)的尺寸设定为30 cm×30 cm×30 cm。由于粒子能够很好地展现混凝土破碎与抛掷等过程，因此本模型中无需添加单元失效准则。

模型主要在LS-DYNA软件中建立，结合关键字的添加与修改，大体可分为三步：(1) 在LS-DYNA中建立FEM模型；(2)采用solid center建立SPH模型；(3)在LS-DYNA中求解，用LS-PREPOST处理计算结果，详见图7。在边界处理上，土体外侧及底面设定透射边界，用来模拟该方向的无限区域；土体内侧与混凝土相接触的部分定义面面接触，并考虑摩擦力的影响；对称面处设置SPH对称边界和FEM对称边界，用以约束单元和粒子在对称面处的法向运动；混凝土(近区)粒子与混凝土(远区)有限单元之间采用固连接触的方式定义耦合算法，即通过添加关键字 CONTACT_TIED_NODES_TO_SURFACE来定义。

3.2 结果分析

图8展示了1 kg TNT炸药半埋爆炸时爆坑的形成过程。在爆炸初期，随炸药能量的释放，由于冲击波压力远大于混凝土的动态抗压强度，炸药爆炸产生的高压气体迅速破碎了炮孔周围的混凝土。同时，应力波在混凝土中引起了径向和环向拉应力，导致大量裂隙出现，并向自由面方向扩展。自由面处的混凝土在反射拉应力作用下被拉裂、发生片落，并在爆生气体的作用下产生抛掷。随着爆生气体膨胀和楔入，径向和环向裂隙进一步扩大，抛掷现象加剧，爆破漏斗逐步形成。最后，炸药能量耗散殆尽，其中一部分通过自由面直接扩散到空气中，另一部分能量则被用于破碎、抛掷混凝土以及形成爆炸震动波，爆破漏斗形态渐趋稳定。

为研究爆炸波在混凝土中的传播特性，同时与实测数据作对比，采用1 kg TNT炸药，则测点1～6距炸药中心的比例距离Z依次为0.255、0.515、0.810、1.256、1.879和2.504 m/kg1/3。考虑到各测点曲线峰值相差较大，图9～10仅显示了中间4个测点相关的压力时程曲线和加速度时程曲线。

由图9～10可以看出，炸药起爆后应力波迅速传播，当达到测点后，压力和加速度在很短的时间内迅速达到峰值，然后逐渐衰减，最后趋近于零。不同测点处的压力和加速度有所差异，即随着测点与药包比例距离的增加，测点处的峰值压力和峰值加速度逐渐减小，正压作用时间有所延长。图9中各测点的压力在经历过正压段后出现持续时间较长且峰值较小的负压段，这可能是由于混凝土试样较小，爆炸波传播到四周边界后经多次反射作用造成的[24]；图10中各测点加速度呈现单峰且峰后出现连续的小幅度数值振荡，这也与实际结果较为吻合[20]。

表2为不同测点处实测数据[20,23]与计算结果的对比，可以看出计算得到的峰值压力与对应的实测数据吻合度较高，误差基本在±12%以内，而峰值加速度相对而言误差偏大，但大多也都在±20%以内。

表2 不同测点处计算结果与实测结果对比Table 2 Comparison between calculated results and measured results at different test points

3.3 HJC参数敏感度分析

模型参数的敏感度分析是指研究一个或多个不确定性参数的变化对计算结果产生的影响，即模型对某个参数或某组参数变化反应的敏感程度。通过敏感度分析，既可以识别对模型计算结果起决定性作用的参数，又可较好地评价参数偏差对计算结果产生的影响，这在本构参数确定与反演等方面有着广泛的应用。

考虑到混凝土采用的HJC本构参数涉及到21个之多，加上状态方程参数大多沿用文献[9]，因此想要进一步提高模拟的精度，还需对各参数的敏感度进一步分析。依照每次分析所涉及到的参数个数，可将敏感度分析分为全局敏感度分析和局部敏感度分析，全局敏感度分析需要检验多个参数的变化对计算结果总的影响，同时评价每个参数及各参数之间对计算结果的影响；局部敏感度分析只需每次检验单个参数对计算结果的影响。考虑到HJC本构模型各参数之间耦合紧密，若采用全局敏感度分析，一方面难以确定参数变化所遵循的标准和范围，另一方面要考虑到每个参数及各参数之间对计算结果的影响，计算量太大。因此，本文中采用局部敏感度分析法[10,13-14]，即研究HJC模型中单个参数变化对计算结果的影响。

以0.515 m/kg1/3处的测点作为研究对象，采用控制变量法，将21个参数在表1的基础上分别单独作±40%和±20%的变化，将每次计算的结果(包括峰值压力和峰值加速度)与表1中初始参数计算结果的比值作为研究参数敏感度的指标[10,13]。通过大量的数值计算，得到C30混凝土HJC模型不同参数变化率下引起的峰值压力和峰值加速度的变化率，如图11～12所示。

可以看出，对峰值压力(Pm)影响比较大的参数有9个：ρ、fc、A、B、G、Pcrush、μcrush、Plock和μlock；对峰值加速度(am)影响比较大的参数有7个：ρ、fc、A、B、μcrush、Plock和μlock。他们在参数变化区间内峰值变化率绝对值均≥10%，最高可达50%，而其他参数对结果的影响较小(≤5%)。

为定量评价各参数对计算结果稳定性影响的敏感度，以0.515 m/kg1/3处的测点峰值的变化率M(x)作为目标函数，在某个参数变化率x=xi(i=1，2，3，4，5，对应的xi=-0.4，-0.2，0，0.2，0.4)处，目标函数M(x)对该参数变化率xi的敏感度Si可表示为[13]：

Si=∂M(x)/∂xi

(13)

式中：敏感度Si表示M(x)在x=xi处的导数值，Si的绝对值越大，表示M(x)对xi越敏感。

图11～12中各曲线大多呈线性或者抛物线变化，因此目标函数M(x)可按照多项式进行拟合，则不同参数的目标函数M(x)与参数变化率x的关系可表示为：

M(x)=a0+a1x+a2x2+…+anxn

(14)

将式(14)代入式(13)，可得：

(15)

参数变化率x区间内包含5个参数，即x=-0.4，-0.2，0，0.2，0.4，可取式(15)中绝对值的最大值来表示该参数的敏感度S：

S=|Si|max

(16)

式中：Si为不同参数变化率xi处的敏感度，i=1，2，3，4，5。

根据式(14)～(16)，将图11～12中各曲线数据利用Matlab软件按照多项式进行拟合，拟合精度设在95%以上，之后求导并代入参数变化率xi，最终得到本模型各参数的敏感度S，见表3。

表3 参数敏感度Table 3 Parameter sensitivity

表3定量给出了本模型对计算结果(峰值压力Pm和峰值加速度am)的敏感度，可见在HJC本构所有参数中，敏感度S>0.4的参数有：ρ、fc、A、B、G、Pcrush、μcrush、Plock、μlock和N，这同图11～12中峰值变化率大于或等于10%所对应的参数基本上能保持一致。

相比于侵彻类与碰撞类问题的模拟[10, 12-14]，对爆炸类模拟结果影响较大的HJC参数有很多，且对于不同的研究对象影响的程度也不同，其中除了极限面参数A和B之外，还涉及到状态方程中的三个参数μcrush、Plock和μlock，为此在模拟爆破类问题时，这些敏感度较大的参数需要结合实验结果来确定，而对于敏感度较小的参数可采用本文中所给的值，以减少重复实验确定参数的繁琐过程。

4 结 论

(1)采用SPH-FEM耦合法模拟炸药在混凝土中爆炸，能直观地展示了混凝土爆坑形态发展的全过程，固连接触的方式很好地解决了SPH法在处理边界问题上的缺陷，同时SPH-FEM耦合法计算的效率也较高。(2)给出了C30混凝土HJC模型各参数的确定过程，并基于确定的参数使用SPH-FEM耦合法对混凝土爆破成坑现象进行了模拟，计算结果与实测数据比较吻合，故所提出的HJC本构参数比较合理。(3)HJC本构参数对峰值压力和峰值加速度均有较大影响的参数主要有ρ、fc、A、B、μcrush、Plock和μlock，故在模拟爆破类问题时，这些敏感度较大的参数需要结合专门实验加以确定，而对于敏感度较小的参数可采用本文建议值。