钢筋混凝土梁冲击试验数值模拟研究

姜 华，贺拴海，王君杰

钢筋混凝土结构被广泛用于土木、海洋和防护工程，由于结构冲击实验比较复杂且在很多情况下不经济，开展此类结构的冲击数值模拟对研究结构的动态响应具有重要意义。近10年来，伴随着计算机技术的快速发展，采用三维有限元模拟钢筋混凝土结构在冲击荷载下的响应已引起了越来越多研究者的兴趣。

回顾采用三维有限元技术数值模拟钢筋混凝土结构受冲击作用时，发现采用的混凝土本构模型比较简单以致无法反映混凝土材料复杂的力学行为。例如，Von-Mises(VM)和Drucker-Prager(DP)模型用到了巨型钢筋混凝土梁受落重的冲击模拟［1－2］，受压侧为双线性联合受拉侧拉断的模型用于数值模拟钢筋混凝土梁的非线性响应［3］。然而，这些模型无法反映实验中观察到的由材料损伤引起的尤其在峰值强度之后所在区域出现的材料刚度退化现象。此外，K＆C模型用到了受横向冲击荷载作用的混凝土柱上，这个模型无法预测材料的非弹性体积膨胀现象［4］。

基于这些背景，为了验证采用高级混凝土本构模型数值模拟的精度，研究有关混凝土材料的几个重要行为对冲击响应的影响以及在复杂结构中推广运用高级混凝土模型，本文采用LS－DYNA程序对已有钢筋混凝土梁的冲击实验进行了数值模拟，其中弹塑性损伤帽盖模型［5－6］用于混凝土材料，该模型可以反映实验中观察得到的混凝土材料的主要宏观行为，而弹塑性模型用于钢筋材料模拟。

1 试验介绍

日本国防研究所Fujikake等［7］开展的钢筋混凝土梁冲击试验用于验证弹塑性损伤帽盖本构模型和讨论相应的数值模拟技术。钢筋混凝土试验梁截面尺寸为1 700 mm×250 mm×150 mm，梁体为双筋梁、上下层均采用2根直径为16 mm的钢筋、钢筋截面积为397 mm2、屈服强度为426 MPa，箍筋采用直径10 mm钢筋，布置间距为75 mm，屈服强度为295 MPa。根据上面的截面形式，构件达到弯曲抗力时的剪力为91.1 kN，而设计抗剪能力232 kN，因此梁体发生弯曲破坏。

混凝土的级配骨料最大尺寸为10 mm，试验在试件浇注并养护70天后进行，测得的单轴抗压强度为42.0 MPa。试验中浇注4片相同尺寸和配筋的梁，供4个工况使用。

图1 钢筋混凝土梁落锤冲击试验装置Fig.1 Drop hammer impact test setup of reinforced concrete beam

冲击试验采用图1所示的自由落锤方式，落锤总重400 kg，冲头采用直径90 mm的半球形刚体，冲头和冲击体之间放置冲击力传感器以测量冲击梁体时产生的冲击力;梁体放置在支撑间距为1 400 mm的特定支撑装置上，梁体挠度时程采用激光测试仪测量，测量时在梁体贴橡胶薄片。试验时，系统对力和挠度数据的采样频率为100 kHz。试验时分共分为4个工况(表1)，每个工况试验前落锤离梁高度分别为0.15 m、0.3 m、0.6 m和1.2 m，以使梁出现不同的损伤和破坏状态。

表1 钢筋混凝土梁落锤冲击试验工况列表Tab.1 The impact test of reinforced concrete beam

2 弹塑性损伤帽盖模型

由于钢筋混凝土结构在冲击荷载下的力学行为复杂和混凝土材料本身的复合特性，已提出了多种混凝土本构关系模型。瞬态软件LS-DYNA、AUTODYN和ABQUS材料库中的混凝土本构关系模型按塑性变形计算方法总体上可以分为三大类:第一类采用关联流动准则计算塑性应变增量，塑性体积应变增量基于增量流动准则得到，因此可以反映低静水下出现的由剪切力引起的塑性体积膨胀现象，代表模型有Mohr-Coulomb、DP、弹塑性帽盖［6］和弹塑性损伤帽盖模型［5－6］;第二类采用Prandtl-Reuss流动准则 (VM准则用作塑性势函数)计算塑性应变增量，塑性体积应变由状态方程(EOS)［5］得到，塑性体应变增量独立与增量流动准则，换句话时偏应力和静水力响应不耦合，这种情况下无法反映膨胀行为，代表模型有 HJC［8］、伪涨量［6］、K＆C［6，9］、RHT［10］本构等，这些模型被广泛用于模拟高速冲击作用下的混凝土;最后一类采用非关联准则计算塑性应变增量，其中屈服函数不同于塑性势函数，这样膨胀行为得到了有效控制，代表模型有塑性损伤模型［11］。实际上如果在第一类的本构屈服面上外加帽盖，通过强化准则收缩帽盖同样可以控制膨胀现象。

对于承受低速冲击的混凝土来说，混凝土材料的重要行为如剪缩，峰值强度前后的膨胀［12］，峰值强度前的强化，峰值强度后的应变软化和刚度退化需得到反映，LS-DYNA中的弹塑性损伤帽盖模型可以较好反映这些特性以及不可恢复的变形，因此被用于模拟混凝土。

弹塑性损伤帽盖模型在LS-DYNA材料库中称为*MAT_SCHWER_MURRY_CAP MODEL模型，在弹塑性帽盖模型基础上扩展得到，由SCHWER和 MURRY提出，该模型对塑性流动和损伤积累采取分开处理方式，其中塑性流动，受剪应力控制，会导致永久变形但不会出现弹性模量退化，而损伤导致弹性模量和强度的不断退化，引入帽盖是为了模拟材料中空隙破坏引起的塑性体积变化，因此该模型在这里称为弹塑性损伤帽盖模型。

模型失效面由剪切失效面和帽盖面组合成，两种失效面之间采用光滑过渡，剪切失效面为三不变量模型，帽盖面由两部分函数组成，其值为1或者为椭圆。

由于理想塑性响应为混凝土材料在高静水压下的特征，并不是混凝土在拉伸和低静水区域的特征，因此需要模拟拉伸和低静水压缩区的应变软化。模型分别采用拉伸和压缩损伤描述拉伸和压缩区的应变软化和弹性模量退化，其中损伤认为由微裂缝的增长和合并导致。

在延性区，采用联合各向同性损伤和塑性变形的弹塑性损伤，而在脆性区，采用弹性损伤忽略了塑性。

损伤定义基于连续损伤力学框架，损伤和损伤能释放率为热力学共轭力，满足Clausius-Duhem耗散不等式。因此定义了损伤面用于确定损伤加载是否发生，当能量模(为损伤能释放率的函数)，超过损伤阀值时，损伤开始和积累，这里可以称为损伤准则［13］，损伤变量的演化根据正交法则。损伤处理方式非常类似于基于热力学框架的经典塑性力学，其中应力张量和塑性应变张量热力学共轭，屈服面定义为应力张量或者不变量的函数，塑性应变的演化根据塑性流动法则。有关该模型的详细介绍和参数确定可参考文献［14］。

3 数值模拟冲击试验

数值模拟时先对图1所示的冲击试验装置建立几何模型、划分网格并赋予相应的材料属性，然后施加边界条件，接触条件和输出控制参数等。

混凝土材料采用弹塑性损伤帽盖模型，28天圆柱体抗压强度取30.0 MPa(因原试验没有测试该强度，这里为经计算比较后的推测值)，模型强度参数采用C30混凝土，参数取值见表2。

表2 弹塑性损伤帽盖模型材料参数Tab.2 Parameters of elastoplastic damage cap model

混凝土单元类型采用8节点六面体单元，积分方式采用单点积分，单元网格最小尺寸为10 mm，最大尺寸为12.5 mm，全梁共划分55 300个体单元;纵向钢筋和箍筋材料采用弹塑性随动强化模型［14］，弹性模量和强化模量分别为200 GPa和3 GPa，单元类型采用杆单元，共计1 690个，钢筋的网格尺度和混凝土单元一致，以使钢筋单元和混凝土单元共节点;冲头采用刚性材料，单元采用4节点四面体单元，网格最小尺度为2 mm(位于冲头最前端与梁体接触部分)，单元数共计25 944个，冲头最前端球面部分和图1中的实际冲头尺寸完全一致，因冲头以外部分不影响接触碰撞采用图2中类似弹体的几何体代替落锤，两者满足质量相等。冲头和梁体之间采用自动面对面接触，支座装置和梁体直接设定自动点对面接触(支座装置材料采用刚性材料并赋予竖向约束)，两者的初始接触距离设为0以反映支座装置对梁体的竖向约束作用。数值模拟时为了节省计算机时，冲头和梁顶面的初始距离设为1mm，将表1给出的冲击速度赋给冲头并施加重力加速度，接触力和跨中挠度输出频率采用和试验数据相同的采样频率为100 kHz，输出响应的时间间隔为10－5s。

图2 钢筋混凝土梁冲击试验计算模型Fig.2 Numerical model of impact test of reinforced concrete beam

图3 至图6为采用弹塑性损伤帽盖模型模拟得到的撞击力和跨中挠度与试验值对比，其中撞击力时程均出现多次振荡。当梁的运动速度大于碰撞体时，撞击力减小出现下降段，之后落锤和梁速度继续减小，当梁体速度慢于落锤速度时，撞击力力开始出现上升段，上述过程持续进行直至出现落锤反弹和分离，不再有接触、撞击力变为0。

从撞击力对比可见计算得到的撞击力峰值和持续时间与实验情况均吻合较好，且随着冲击速度的增大，模拟结果和试验值之间的误差减小。从梁体跨中挠度对比可知，模拟结果和试验值在峰值挠度之前吻合非常好，下降段主要由梁体弹性变形释放产生，预测值和实测值之间的误差随落锤速度增大而减小。

图7(a)为各工况试验后梁体的破坏情况，其中工况1冲击速度较小，梁体出现了多条由梁底部往上开展的竖向裂缝，主要由梁体的弯曲破坏引起;工况2同样出现了多条竖向裂缝，但离跨中较远方向竖直裂缝慢慢倾斜;工况3出现了斜裂缝，由剪力和弯矩共同作用下的主拉应力引起，梁体出现了斜剪破坏;工况4除了出现斜裂缝之外，裂缝宽度变宽，在靠近跨中附近裂缝宽度最大，并且直接与落锤碰撞部位的混凝土出现了脱落，梁体出现了较大的不可恢复挠度。

图7 试验后的梁破坏形态对比Fig.7 Comparsion damage models of four beams after test

图7 (b)为数值模拟得到的各工况下梁体拉损伤分布图，因拉损伤分布和混凝土裂缝分布状况直接联系，拉损伤云图和实验观测得到的裂缝分布总体上吻合较好。

4 混凝土本构模型讨论

弹塑性损伤帽盖模型(原模型)经简化后可以退化为工程上常用的几种模型和多种新模型。对原模型的简化包括子午线函数、偏平面角隅函数、帽盖面和损伤软化四个方面，其中压缩子午线可以简化为线性函数，角隅函数可以由Rubin模型退化为直线和圆弧形式，通过向静水轴远端平移帽盖面可以达到取消帽盖面目的，取消软化段后可变为帽盖1模型，此时材料不存在刚度退化现象，若在此基础上将Rubin角隅函数退化为圆弧则该模型变为帽盖2模型。表3给出了各种退化模型和原有模型之间的关系，如前面提及的DP模型对应于表中DP(无软化)模型，由原模型子午线，角隅函数作简化后并取消帽盖面和软化段得到。考虑到结构出现较严重损伤情况对工程具有更重要的指导意义，下面各模型的比较和讨论均基于工况4展开。

表3 弹塑性损伤帽盖模型各退化模型对比Tab.3 Comparisons of the elastoplastic damage cap model and its degenerate models

4.1 应变软化段影响

混凝土材料峰值强度之后的应变软化行为包括材料强度退化和刚度退化，材料此时出现宏观裂缝，并不断开展导致不稳定断裂，最终达到破坏状态而丧失凝聚力。这一阶段力学行为最为复杂，已有的理论包括弹塑性理论，损伤力学和断裂力学均用到了描述这些现象。这里将软化段基于自由能释放形式的原模型和取消软化段的简化模型(对应表3中的帽盖1)作对比以研究应变软化段对模拟结果的影响，取消软化段可通过将软化参数设为0实现。图8所示为两种模型模拟得到的撞击力和挠度与试验值对比:简化模型对应的撞击力持续时间比原模型对应持时缩短，撞击力峰值稍大;两者模拟得到的挠度结果和试验值在0.01 s时刻以内吻合较好，主要由于两种本构关系曲线在峰值点以前完全重合，在该时刻以后简化模型得到的挠度小于原模型相应值。其它简化模型取消应变软化段和有软化的结果对比与上述现象类似，可见材料层面的应变软化行为最终会影响到结构的宏观响应，描述材料的应变软化现象在混凝土本构中是必要的，后面各工况比较时均考虑应变软化行为。

4.2 帽盖面影响

帽盖面用于描述静水压力导致的塑性体积屈服现象，塑性体积应变计算时采用关联流动准则。当帽盖屈服面与静水轴初始交点位置参数X0取较大值时可以退化为弹塑性损伤无帽盖模型，在强度准则方面意味着不考虑静水力导致的塑性体积屈服，在材料变形方面意味着忽略混凝土材料内部的空隙压实破坏导致的塑性体积变形。图9所示为原模型和取消帽盖后的模型(对应表3中的弹塑性损伤模型)模拟结果对比，可知取消帽盖后出现撞击力时程围绕试验值较大幅度振荡，最大撞击力幅值减小但撞击力持时基本不变，因此撞击力和试验值间的偏差增大;挠度对比可知取消帽盖后的模型得到的挠度比原有模型值和试验值稍大。可见，塑性体积应变描述在混凝土变形中占据了重要地位，直接影响到结构的宏观响应如撞击力，后面比较和应用中均保留帽盖面描述塑性体积变化。

4.3 偏平面形状影响

破坏面的偏平面形状受角隅函数控制，模型A由原模型的角隅函数作VM简化得到。由图10可知得到的撞击力峰值稍大于原模型值和试验值，并且振荡幅度增大，与试验值之间的吻合度略差于原模型;挠度对比可知在0.012 5 s以前几乎无影响，之后幅值稍小于原模型对应值。对比结果说明当VM角隅函数用于偏平面时，撞击力计算值与试验值之间的误差变大，但对跨中挠度影响较小。

模型B由原模型的角隅函数作MC简化得到，其撞击力和挠度、时程和原模型对应值大致相同(图11)。联合上面的简化结果可见，在保持原模型其它行为不变的情况下，偏平面角隅函数对梁撞击响应的影响较小，并且MC角隅函数得到的结果要比VM角隅函数得到的结果要好。

4.5 子午面形状影响

这里的DP模型和工程上常用的DP模型有所区别在于前者考虑了应变软化现象，可由原模型子午线作线性简化，偏平面采用VM角隅函数得到。由图12给出的对比可知，模型DP模拟得到的撞击力时程峰值和振荡幅度均变大，且荷载持续时间缩短，得到的挠度曲线在0.0125s以前几乎无影响，之后幅值小于原模型对应值和试验值。因为模型DP同模型A一样也对角隅函数作了圆弧简化，故得到的挠度继承了模型A的缺点。

类似DP模型，这里讨论的模型MC同样考虑了应变软化，由原模型子午线和偏平面图像均做线性简化后得到。从图13得到的结果对比可知，模型MC得到的撞击力和挠度和模型DP类似，撞击力时程持续时间稍大，挠度稍大于试验值和原模型对应值。然而该模型相比模型DP而言，在计算梁体挠度方面与实验值之间的误差稍小。

4.6 混凝土划分技术影响

考虑钢筋划分混凝土网格方法有两种:一种方法在钢筋单元和混凝土单元之间采取共节点方式，这种方法很难用于复杂结构如桥墩和箱梁.。第二种方法采用 *CONSTRAINTED LARANGE IN SOLID［6］技术将钢筋单元耦合到混凝土单元，即梁单元所在从组耦合于混凝土所在主组，有关该技术详细的描述可参见参考文献［6］。两种划分方式的混凝土本构均采用原模型，且混凝土单元网格大小相同，采用耦合方法建模时不合并具有相同坐标的钢筋节点和混凝土节点，取而代之的是将钢筋单元组成的集合绑定于混凝土单元组成的集合中。两种划分方法得到的结果和实验值比较见图14，共节点方法得到的响应要优于耦合方法得到的响应，然而，考虑到采用公节点方法很难用于复杂结构的混凝土单元划分而耦合方法非常方便，如果忽略两者间的较小差别，耦合方法可以作为共节点方法的替代物。

图14 钢筋和混凝土建模方式对撞击力和挠度影响对比Fig.14 The influence of finite model of steel and concrete to collision force and midspan deflection

5 结论

为了验证弹塑性损伤帽盖模型的精度，展开了钢筋混凝土梁冲击试验数值模拟，预测得到的冲击力、梁体跨中挠度以及梁体试验后的裂缝类型和试验情况均吻合较好。混凝土材料的损伤行为、帽盖面、子午平面和偏平面形状对梁冲击响应影响如下:

总的说来，退化模型得到的预测结果不如原模型合理;损伤、帽盖面和子午线面形状对梁撞击响应的影响较大;偏平面形状影响不大，MC角隅函数优于VM角隅函数，当较难获得混凝土剪切状态下的强度数据时，可以用作Rubin缩放函数的替代物;共节点方法得到的结果要优于耦合方法得到的结果，如果忽略小的差别，后者可以用于复杂结构混凝土网格划分。

［1］Bhatti A Q，Kishi N，Mikami H，et al.Elasto-plastic impact response analysis of shear-failure-type RC beams with shear rebars［J］.Journal of Material and Design，2009，30(3):502－510.

［2］ Bhatti A Q，Kishi N，Konno H，et al.Elasto-plastic dynamic response analysis of prototype RC girder under falling-weight impact loading considering mesh size effect［J］.Structure and Infrastructure Engineering，2012，8(9):817 －827.

［3］Kishi N，Bhatti A Q.An equivalent fracture energy concept for nonlinear dynamic response analysis of prototupe RC girders subjected to falling-weight impact loading［J］.International Journal of Impact Engineering，2010，37(1):103－113.

［4］Thilakarathna H M I，Thambiratnam D P，Dhanasekar M，et al.Numerical simulation of axially loaded concrete columns under transverse impact and vulnerability assessment［J］.International Journal of Impact Engineering，2010，37(11):1100－1112.

［5］Murray Y D，Lewis B A.Numerical simulation of damage in concrete.Technical report submitted to the Defense Nuclear Agency by APTEK［R］.DNA－TR－95－190，1995.

［6］ LS-DYNA KEYWORF USER'S MANUAL.Version971［M］.Livermore Software Technology Corporation，2007.

［7］Fujikake K，Li B.Impact response of reinforced concrete beam and its analytical evaluation［J］.Journal of Structure，ASCE，2009，135(8):938－950.

［8］Holmquist T J，Johnson G R，Cook W H.A computational constitutive for concrete subjected to large strains，high strain rates and high pressure［C］.Proceeding of the Fourteenth International Symposium on Ballistics，1993:591 －600.

［9］ Malvar L J，Crawford J E，et al.A plasticity concrete material model for dyna3d［J］.Int.J.Impact Engineering，1997，19(9 －10):847 －873.

［10］ Riedel W，Thoma K，Hiermaier S，et al.Penetration of reinforced concrete by BETA2B2500 numerical analysis using a new macroscopic concrete model for hydrocodes［C］.9th International Symposium，Interaction of the Effects of Munitions with Structures，Berlin-Strausberg，IBMAC，1999:315－322.

［11］ ABQUS User Manual， version 6.4. ABQUS Inc.，Pawtucket，Rhode Island，2004.

［12］ Kupfer H，Hilsdorf H K，Rusch H.Behavior of concrete under biaxial stresses［J］.ACI Journal，1969，66(8):545－666.

［13］ Simo J，Ju J.On strain-based continuum damage models I，formulation［J］， International Journal of Solids and Structures，1987，23(7):821 －840.

［14］姜 华，贺拴海，王君杰.混凝土弹塑性损伤帽盖模型参数确定研究［J］.振动与冲击，2012，31(15):132 －139.

［15］ Cowper G R，Symonds P S.Strain hardning and strain rate effects in the impact loading of cantilever beams［R］.Brown University，Technical Report No.28，1957.