爆炸冲击作用下三种混凝土本构模型对比研究

杜 闯，宋 帅，张江鹏

(1.河北工业大学 土木与交通学院， 天津 300401; 2.河南省特种防护材料重点实验室， 河南 洛阳 471023)

1 引言

数值模拟是结构抗爆研究领域中常用的研究方法，数值模拟的准确性很大程度上取决于混凝土动态本构模型的选取。近30年来，国内外研究人员通过试验及理论分析等方法对混凝土结构在动态荷载作用下的力学性能展开了一系列研究，并提出了许多混凝土的动态本构模型[1]。如焦楚杰等[2]在ZWT模型基础上，以分形为损伤变量，提出了分形损伤演化的高强混凝土动态损伤本构方程。刘新荣等[3]以ZWT模型为基础，考虑损伤对混凝土本构关系的影响，提出了聚丙烯纤维混凝土动态本构模型。Biani等[4]在经典Perzyna本构模型基础上，建立了混凝土的粘塑性模型。Holmquist等[5]在Ottosen模型中引入了应变率效应，提出了Holmquist-Johnson-Cook模型。Riedel W等[6-8]基于塑性损伤力学，提出了著名的RHT模型。Marlvar等[9]对LLNL模型进行修正，提出了混凝土塑性损伤力学模型K&C模型。如何更好地了解混凝土在爆炸荷载作用下的规律，对研究人员在本构的理解、适用性的认知、本构模型的选择等方面都提出了较高的要求。因此，许多学者对不同的本构模型进行理论比较，通过分析不同本构模型的特点，为本构模型的选择提供了参考意见[10-12]。如李世民[10]对常见的混凝土损伤本构模型进行综合评述，对混凝土本构模型的选取具有重要意义，但该分析仅停留在本构模型的理论方面，没有模拟试验的进一步论证。张若棋等[11]通过理论分析和侵彻模拟对HJC、RHT本构模型进行对比研究，认为HJC本构模型对拉压子午线的描述存在不足，这为混凝土抗侵彻问题本构模型的选取提供了良好建议，但未涉及在爆炸冲击工况下本构模型的研究。蒋轲[12]对HJC、RHT、K&C 3种本构模型进行阐述，通过模拟爆炸试验分别描述了3种本构模型的模拟效果，但分析仅限于对毁伤形态的对比，没有对挠度结果进行比较分析。

为了明晰不同动态本构模型在模拟爆炸时所能突出表现的不同效果，本文对3种常用混凝土本构模型，HJC、RHT、K&C从理论上进行了分析。同时采用LS-DYNA软件对已有文献中的钢筋混凝土板抗爆试验进行数值模拟，分析比较各种本构模型在描述爆炸冲击响应时的毁伤机理。研究结果为了解动态本构模型适用性，达到特定模拟效果而选取更适合的本构模型提供了参考。

2 常见的3种混凝土本构模型

目前，在数值模拟中较常用的混凝土动态本构模型主要有HJC、RHT、K&C 3种本构模型。3种本构模型从极限面、损伤、应变率、状态方程和参数标定5个方面进行对比，如表1所示。

表1 3种本构模型

3种本构模型将混凝土变形分为形状改变和体积改变，形状改变由极限面方程表示，体积改变由状态方程描述。同时，3种本构模型在极限面和体积改变的描述上又具有各自的特点，包括极限面个数，极限面的强化和软化，极限面对压力、应变率、损伤的依赖性和状态方程三阶段性的描述等。下面针对上述3种动态本构模型的优缺点进行简要评述。

HJC本构模型是目前抗爆领域中使用最为广泛的本构模型之一。在极限面的描述上，HJC本构模型仅有一个弹性极限面，且没有考虑偏应力张量的第三不变量影响，也有学者指出在强度模型中，不同损伤程度的强度在一个公式中表达，造成不同压力下，无损伤和完全损伤的强度差异较小[13]。该本构模型认为损伤受等效塑性应变和塑性体积应变的共同影响，且损伤的主要因素是等效塑性应变，考虑了压缩损伤和剪切损伤对强度的影响效应，对损伤的解释比较全面。但在应变率效应的体现中，仅通过固定不变的应变率增强系数和参考应变率来描述应变率效应，并且该模型在描述混凝土拉伸行为时没有考虑到应变率的影响。在本构模型中包含了三段多项式状态方程，考虑了混凝土材料从破碎到孔隙压实的全过程，充分体现了混凝土的力学非线性特征。HJC本构模型各参数物理意义明确，可以通过具体试验获得，国内外对该本构模型研究较多，在缺乏试验的基础上，有学者对参数的确定提出了简易方法[14]，这极大方便了研究人员对该本构模型的使用。

RHT本构模型由HJC本构模型进一步发展而来，并成功应用于混凝土、岩石等脆性材料的爆炸数值模拟。该本构模型的屈服准则采用三个独立的极限面来表达，考虑了偏应力张量的第三不变量影响，弥补了HJC本构模型的不足。RHT 模型的弹性极限面类似于帽盖极限面，其残余强度极限面为压力的指数函数。RHT本构模型可以清晰的描述出弹性极限面与失效面之间的线性强化阶段，同时可以较好的描述混凝土材料在部分损伤和完全损伤条件下继续抵抗剪切变形的特征。此外，RHT本构模型的屈服应力引入了罗德角，能够准确的描述压缩子午线失效强度的变化。在损伤的定义中RHT本构模型仅考虑了等效塑性应变，损伤只考虑等效应力超过失效面之后的软化阶段的损伤量，没有考虑等效应力在弹性极限面与失效面之间的线性强化阶段的损伤，这是不足之处。状态方程采用p-α方程，反映了混凝土在不同孔隙压实程度下的体积改变特征。然而，RHT本构模型对断裂能的描述存在不足，在拉伸荷载作用下无法细致的表现出应变率效应带来的影响。在应用中，由于该本构模型参数较多，且大多数参数难以从试验中获得，其参数难以准确标定。

K&C本构模型同RHT本构模型相似，都考虑了3个极限面和偏应力张量的第三不变量影响，K&C 本构模型采用分段函数来表示拉压子午线之间的关系，对拉压子午线的描述更加细致。在损伤表达方面，K&C本构模型不仅考虑了剪切变形损伤并进一步划分，同时考虑拉伸、压缩不同情况的损伤，但没有考虑混凝土体积膨胀引起的损伤软化。在应变率效应的表达中采用了径向放大法，其应变率增强因子可以通过试验确定，具有一定的合理性。但该本构模型自身不包含状态方程，需要另外采用8号状态方程来定义。该本构模型的诸多参数可以由试验来确定，即使在缺乏试验的情况下，LS-DYNA 软件可以根据无侧限抗压强度自动生成相关参数，极大方便了参数的标定过程。

3 数值模拟对比分析

综上所述，以上3种动态本构模型均较为复杂，模型具有各自优点和局限性。为了进一步分析3种动态本构模型的数值模拟效果，本文对已有文献试验进行模拟对比分析。

3.1 计算模型

采用LS-DYNA有限元软件对文献[15]中的RC板爆炸试验进行数值模拟。如图1所示，试验采用临空板爆炸方式。其中被爆炸的RC板尺寸为1.25 m×1.25 m×0.05 m，板两侧被夹板约束，炸药位于板的中心正上方，距离板0.5 m，TNT装药量为0.64 kg，混凝土板配筋方式为单层双向Φ6@75。

图1 抗爆试验构件及钢筋布置图

在相同条件下，分别采用HJC、RHT、K&C3种材料模型进行数值模拟，对比不同本构模型下钢筋混凝土板的破坏形态和最大挠度。

有限元建模时，钢筋材料选择线性强化模型MAT_PLASTIC_KINEMATIC_TITLE、混凝土材料分别选择MAT_111 (HJC)、MAT_272 (RHT)、MAT_72R3(K&C)3种本构模型。如图2所示，建立1/4模型，在对称面上施加法向约束。在混凝土的一侧上下建立刚性、全约束夹板，来约束构件的垂直、转动位移，夹板使用MAT_RIGID_TITLE刚性材料。由于爆炸荷载过大可以忽略钢筋与混凝土之间的滑移影响，二者采用共节点方式建立联系。分析使用流固耦合法，钢筋混凝土板采用Lagrange网格，炸药和空气采用Euler网格，单元使用多物质ALE算法。

图2 钢筋混凝土板模型示意图

根据文献结论[16]，在自由空气场的爆炸模拟中比例距离大于0.42 m/kg1/3时，炸药和空气可采用尺寸为25 mm的网格。在该工况下炸药与靶板的比例距离为0.591 m/kg1/3，因此对炸药和空气采用10 mm网格尺寸可以满足计算精度要求。对混凝土网格分别取15 mm、10 mm、5 mm进行计算，对比RC板几何中心处的最大挠度、最大压力。结果如表2、表3所示，得到的挠度、压力结果差异不大。但随着网格尺寸的减小，对混凝土损伤失效裂纹的扩展描述更细致，为了清晰呈现构件的裂纹扩展现象，混凝土采用5 mm网格尺寸。

表2 不同网格尺寸下的最大挠度

表3 不同网格尺寸下的最大压力

3.2 模型参数

本文中的基本材料参数与原试验材料参数一致，其他参数基于前人总结并结合参数意义进行调整[17-18]。在失效准则关键词的定义中，结合不同本构模型的模拟效果，以模拟结果最接近试验的原则进行确定，其中HJC以最大抗拉强度不超过4.2 MPa，RHT和K&C以最大拉应变不超过0.003为侵蚀标准。

根据文献[14]给出的方法，可计算出混凝土的压实体积应变、压碎体积应变等参数。抗压强度和抗拉强度采用原试验材料参数值。由于试验中的混凝土强度接近35 MPa混凝土，其他参数可引用Riedel[6]提出的35 MPa混凝土参数取值。将以上参数代入到3种本构模型中。K&C模型的参数较多，本文仅对泊松比,密度、单轴抗拉强度等基本参数进行定义。普通混凝土的泊松比在0.14～0.2，计算中一般取0.2。剪切强度增强因子与应变率之间的关系参考LS-DYNA手册[19]中混凝土相关参数来确定。以此来保证3种本构模型在参数上统一、自洽。主要参数取值见表4、表5、表6。

表4 HJC本构模型主要参数

表5 RHT本构模型主要参数

表6 K&C本构模型主要参数

3.3 加载方式验证

试验中的爆炸可以近似看作自由空气场中的爆炸。当入射冲击波随空气传递到靶板表面时，经障碍物表面反射后形成反射波，反射波与入射波相互作用，在靶板表面产生巨大的压力，该压力峰值称为反射超压峰值。通过数值模拟获取RC板迎爆面的入射波超压和反射波超压，以验证加载方式的可靠性。

采用与上述工况相同的装药量，模拟自由空气场中的爆炸。提取比例距离为0.591 m/kg1/3处的压力作为RC板表面中心处的入射波压力。取临近RC板中心处的空气压力作为反射波压力。如图3所示，爆炸产生的冲击波到达板中心，使该处压力提升至最大，随着能量的不断损耗，压力值迅速下降至0附近，并且逐渐趋于稳定，曲线走势符合爆炸超压的基本规律。

图3 超压时程曲线

在相同比例距离下，根据Henrych给出的超压公式[20]得到入射波超压峰值为1.99 MPa，模拟值为2.28 MPa，基本一致。根据美军防护设计手册[21]给出的经验预估曲线，该条件下的反射峰值超压提高6.7倍，模拟值为10.9 MPa，提高了4.8倍，误差约为28.4%。可以看出，这种由炸药状态方程施加爆炸载荷的模拟结果具有一定可靠性。

3.4 破坏形态对比

炸药爆炸的瞬间，爆炸产物压缩周围空气而形成冲击波[22]，随着冲击波传播距离的增加，波阵面不断扩大，当延伸到RC板表面时，在板的迎爆面形成较强的压缩波，对板造成超压破坏。当板中应力波传递到下表面时，在背爆面发生反射形成较强的拉伸波，使RC板的背爆面主要产生受拉破坏。数值模拟中将板的迎爆面和背爆面分别与试验做比较。在模拟中发现，当计算时间超过3 ms时，靶板单元的加速度曲线在0附近发生微小扰动，爆炸响应基本结束，RC板的破坏程度不再发生明显变化，因此模拟采用计算到3 ms时的结果作为破坏的最终结果。

在破坏形态的模拟中，对裂纹开展现象的表示方法通常有侵蚀失效法和损伤失效法。侵蚀失效法是通过添加失效关键字对变形或受力较大的单元进行删除来得到裂纹，这种方法得到的结果主要基于本构模型对状态方程和屈服面的表达。损伤失效法是通过损伤度来体现裂纹扩展等破坏现象，其结果与本构模型中损伤演化方程有直接联系。下面分别用2种方法来比较板的破坏，分析本构模型中不同方程的影响效果。

3.4.1侵蚀失效法结果对比

迎爆面破坏形态如图4所示[15]。试验板的迎爆面在中央呈现出一条较粗的弯曲裂缝。板中心出现剥落破坏，四周伴有明显的环形裂纹和少量斜裂缝。

在模拟结果中，HJC本构模型的模拟结果没能体现出中部较粗的横向裂缝。在中心处出现少量层裂脱落，在靶心的周围，出现了环形裂纹，4个角部开展斜裂缝。这些特点与试验板的迎爆面一致。RHT本构模型复现了试验板的中央横向裂缝，在中心处，出现损伤坑，在靶板的角部，没有出现斜裂缝。K&C本构模型没有出现中央横向裂缝，在中心处，模拟出损伤坑，在靶板的角部及周围没有体现出斜裂缝和环形裂纹。3种本构模型的屈服面中，HJC模型的偏平面呈圆形，由于不考虑J3的影响，其失效曲线形状不变，表现出各向同性特征。在对称荷载中，其破坏形态也体现出一定对称性，对于环形裂纹的描述较好。

图4 迎爆面破坏形态图

背爆面破坏形态如图5所示[15]。试验板的背爆面破坏主要表现为板中心处混凝土的剥落与崩塌，其剥落区半径为120 mm[15]。板的中央出现若干弯曲裂缝，四周呈现由中央向外扩展的径向裂缝及细小环形裂纹。在接近角部，有少量斜裂缝产生。HJC本构模型模拟出了试验板中心处的剥落破坏，其中剥落圆半径为107 mm，误差为10.8%。在靶板中央，HJC本构模拟出弯曲裂缝，在靶板的四周也描述出了径向裂缝和环形裂纹。在接近角部，HJC本构体现出向角部延伸的斜裂缝。RHT本构模型对中心处的剥落破坏也有良好的体现，其中剥落圆半径为70 mm，误差为41.7%。RHT本构在靶板中央模拟出了十字交叉裂缝，在四周也出现了少量环形裂纹。K&C本构模型也模拟出了试验板中心剥落区的破坏效果，其中剥落半径为155 mm，误差为29.2%。可以看出，K&C本构在靶板中央模拟出了横向裂缝，在角部出现较为密集的斜裂缝，然而对环形裂纹的开展没有描述。在屈服面的表达中，K&C和RHT模型定义拉伸和压缩子午线不同，在相同静水压力下表现出材料抗压强度大于抗拉强度。在被爆面的破坏中主要受拉伸波影响，其静水压力为负值，K&C和RHT模型的屈服范围远小于静水压力为正值的情况。因此2种本构模型背爆面的破坏程度远大于迎爆面。

图5 背爆面破坏形态图

3.4.2损伤结果对比

为了对比3种本构模型损伤表达能力的不同，利用损伤度来对RC板破坏情况进行对比。其中损伤度由0-1表示破坏程度逐渐升高。

迎爆面损伤云图如图6。HJC模型的损伤度通过等效塑性应变和塑性体积应变的累积来定义。其结果在中央呈现出由弯曲变形导致的横向损伤带，对应于中央弯曲破坏。中心呈现出近似圆形的损伤区，对应RC板中心处破坏。由于整体损伤值小于0.3，表示材料发生轻微的塑性变形，不足以体现混凝土的层裂、成坑现象。RHT模型通过等效塑性应变来定义损伤度。其结果在中央呈现出横向弯曲损伤，损伤值为0.5，表示为裂缝产生。RC板中心处的损伤度为1，即材料完全破碎。RHT模型在四周也模拟出放射状损伤裂缝。整体损伤特征与试验接近。K&C模型将损伤定义为受拉损伤、受压损伤和体积变形损伤。目前LS_DYNA软件中没有直接给出K&C模型损伤度的提取方式，一般通过塑性应变值来表达损伤，其变化范围为0～2。可以看出，在RC板中心处塑性应变大于0.6，表示为材料的层裂破坏。在四周出现大范围的塑性破坏区，其塑性应变值大于1，对应于板的环形破碎裂纹的产生。在四个角部塑性应变值接近1.8，表明RC板在角部产生严重的破裂现象。整体损伤特征与试验一致。相比RHT模型，K&C本构模型在损伤的描述上考虑三向等拉时的体积变形损伤，这使强静水压力对靶板单元造成的破坏得到良好的体现。同K&C本构模型相似，HJC模型也考虑体积应变带来的损伤，但K&C模型采用表控状态方程，可以更准确描述压力与体积之间的关系，对体积损伤的描述也更完备。由于HJC模型在不同损伤程度下都只在一个极限面中表达材料的屈服，而K&C模型RHT模型有3个极限面，在材料损伤时极限面的变化范围大于HJC模型。因此在相同条件下HJC板的损伤程度要远小于K&C模型和RHT模型。

图6 迎爆面损伤云图

背爆面损伤云图如图7。HJC模型没能体现出损伤现象。RHT模型模拟出十字交叉损伤裂缝，在RC板中心处也出现圆形的损伤破碎区。在四周及角部呈现出斜裂纹的开展以及放射状损伤裂缝。整体损伤特征与试验基本一致。K&C模型的损伤区集中在靶板中部，在角部也出现放射状损伤裂缝。HJC模型通过一个“最大拉伸压力”来定义拉伸极限，使该本构不能很好的表现出拉伸损伤裂纹的扩展现象。与RHT模型相比，K&C模型多考虑了体积拉伸损伤，在损伤破坏的程度上要高于RHT模型。

图7 背爆面损伤云图

对比3种本构模型的模拟效果，HJC本构模型可以结合失效准则模拟出靶板的弯曲裂缝、环形裂纹、角部斜裂纹以及中心剥落圆破坏，其中对剥落圆的描述最接近试验。但该本构不适合以损伤度来描述具体的破坏特征。RHT本构模型结合失效准则可以模拟出RC板的弯曲裂缝、环形裂纹以及中心剥落圆。相比之下RHT本构模型更适合以损伤度来体现裂纹扩展等破坏现象。K&C本构模型结合失效准则可以复现RC板的弯曲裂缝，中心剥落圆以及斜裂缝的扩展，但在损伤表示方面对裂纹扩展的描述不及RHT模型。综合考虑3种本构模型的模拟效果，HJC本构模型符合的破坏指标最多，对RC板破坏形态的整体描述直观，模拟效果最好，RHT和K&C本构模型也可以一定程度上反映RC板的破坏特征。

3.5 最大挠度对比

试验中，RC板在爆炸荷载作用下发生塑性破坏，板中央区域出现挠度变形。在数值模拟中，通过比较构件中心的竖向节点位移来描述RC板的最大挠度变形。

图8是3种本构模型下RC板背爆面中心处的加速度时程曲线。从3种本构模型的模拟结果看出，在同一时刻，应力波传递到RC板下表面使该处单元的加速度迅速提高，随后在极短的时间内单元加速度降低，在2.5 ms后3种本构的模拟结果均表现出加速度曲线逐渐平稳并趋近于零的现象，此时爆炸响应基本结束。3种本构模型的加速度曲线走势一致，与实际情况相符。

图8 板跨中加速度时程曲线

图9表示3种本构模型的中心挠度变化。试验中测得板的最大挠度为19.0 mm[15]，模拟计算得出HJC本构下的最大挠度为22.0 mm，误差13.6%，而RHT本构和K&C本构下的挠度呈线性趋势，其最大挠度值均超过40 mm。结合加速度曲线分析，当爆炸引起的动力响应结束时，RC板受到自身惯性作用仍产生部分位移。在计算时间内，HJC本构下的RC板首先达到最大挠度，RHT和K&C本构下的RC板仍处于惯性阶段。

通过对比3种本构模型下的挠度最大值可以看出，HJC本构模型对中央挠度的描述最接近试验，RHT本构模型和K&C本构模型的模拟效果误差较大。

图9 板跨中挠度时程曲线

4 结论

基于对3种混凝土动态本构模型模拟爆炸冲击响应的分析及数值模拟结果的对比，得出以下结论：

1) 在本构模型的理论方面。HJC本构模型仅采用一个弹性极限面来描述材料的屈服准则，RHT和K&C本构模型采用3个极限面来描述混凝土的屈服，更加全面。在损伤上，K&C本构模型更加细致，RHT本构模型比较粗糙。在应变率上，3种本构模型各自表示不同。在状态方程上，K&C本构模型不包括状态方程，需要另外定义，HJC和RHT本构模型包括状态方程。参数标定上，K&C本构模型最容易。

2) 在破坏形态的描述方面。3种本构模型都能够在一定程度上描述试验现象。整体上看，HJC本构模型效果最佳，RHT本构模型最适合通过损伤度来描述裂纹扩展现象。

3) 在最大挠度的描述方面。HJC本构模型对挠度的模拟程度相对接近试验。