钢筋混凝土结构改进型分离式数值模型*

2022-07-11冯晓伟肖桂仲

爆炸与冲击 2022年6期

曾繁，冯晓伟，黄超，徐权，肖桂仲,4，田荣

（1. 中国工程物理研究院高性能数值模拟软件中心，北京 100088；2. 北京应用物理与计算数学研究所，北京 100088；3. 中国工程物理研究院总体工程研究所，四川绵阳 621999；4. 南京理工大学理学院，江苏南京 210094）

钢筋混凝土作为结构材料广泛应用于土木工程中，其力学性能非常复杂，与混凝土基体、钢筋及它们之间相互作用等力学行为紧密相关。当钢筋混凝土结构承受爆炸荷载时，钢筋与混凝土间的粘结应力上升，当超过一定的限制值后，钢筋与混凝土间发生滑移。尽管界面的粘结应力-滑移关系建立在微观尺度层面，却显著影响结构在爆炸荷载下的力学性能。

在钢筋混凝土结构的有限元分析方法中，围绕钢筋的数值建模，常用有限元模型有3 种：(1) 弥散式模型；(2) 嵌入式模型；(3) 分离式模型。

弥散式钢筋模型常用于由均匀分布的钢筋组成的墙体、楼板等表面型结构。该模型将钢筋混凝土简化为均质材料，易于使用。但是对于梁柱结构中箍筋和纵筋等复杂的钢筋几何问题，弥散式模型由于离散精度问题而不被采纳。

在嵌入式模型中，钢筋混凝土构件横截面分为混凝土部分和钢筋部分，依据平截面假定计算截面应变，再根据钢筋和混凝土应力-应变关系和力平衡条件得到单元刚度矩阵，其中钢筋强化效果等效为嵌入在混凝土单元中的附加轴向分量。嵌入式模型使得钢筋的运用完全独立于其几何形状，且不增加数值模型的自由度。该模型假定嵌入的钢筋与混凝土完美粘结，即钢筋节点位移与混凝土节点位移一致。

在分离式模型中，分别对钢筋和混凝土精细建模，混凝土采用实体单元离散，钢筋采用一维杆或者梁单元离散。在数值精度方面，该模型优于嵌入式模型，且通过大量的钢筋混凝土构件爆炸实验验证，故而在构件一级的抗爆分析中常为工程师所青睐。由于用于离散钢筋的单元需布置于混凝土单元的边界，这导致混凝土的网格严重受限于钢筋几何形状。

同时，在分离式模型中通过显式构造界面单元，可以考虑界面的粘结滑移效应。例如，师燕超等在混凝土网格节点和钢筋节点间构造一维的虚拟弹簧单元，通过一维滑动接触模型考虑二者间的粘结滑移效应。Kwak 等通过混凝土网格节点与钢筋网格节点共点方式，构造了无物理尺寸的一维连接单元。De Groot 等在混凝土与钢筋之间定义二维粘结带，考虑平移、旋转等不同类型的相互作用。基于钢筋、混凝土网格节点建立的界面模型能较好地捕捉粘结滑移效应。但是界面模型给计算带来额外的开销。更重要的一点是，基于节点的界面几何建模的复杂度随着工程结构几何复杂度的增加而急剧增长，将耗费巨大的时间人力成本，这会极大地制约该类模型在钢筋混凝土结构级的工程应用。因此针对大型钢筋混凝土建筑结构，如何高效高精度地建立钢筋混凝土数值模型依旧是一个开放性问题。

根据混合物理论，钢筋混凝土的力学行为可由混凝土基体和强化的钢筋两项组分的本构关系来描述，并且无需对钢筋进行显式离散化。该理论具有不同组分间应变相容性的局限性假定，可以通过修正钢筋本构关系，对混合物理论进行修正。本文构建一种基于修正混合物理论的改进型分离式数值模型，以期给出在钢筋本构方程中如何考虑界面粘结滑移机制的简单有效的修正方法，并避免显式构造界面即可考虑钢筋尺度的局部效应；对改进型分离式数值模型进行钢筋混凝土构件级、结构级的实验验证，以探索该模型在宏观大型钢筋混凝土结构有限元分析中的实用性和可靠性。

1 改进型分离式数值模型

参照《混凝土结构设计规范（GB50010—2010）》，本文中以6 层钢筋混凝土框架结构为例详细阐述改进型分离式数值模型的离散步骤：(1) 混凝土结构几何建模；(2) 钢筋骨架线几何建模；(3) 采用实体单元对框架结构进行网格离散；(4) 根据钢筋骨架线几何位置及钢筋特征尺寸信息，对包含钢筋骨架线的网格定义成钢筋混凝土网格，未包含钢筋骨架线的网格定义成素混凝土网格（这种离散方式巧妙地避免了传统分离式建模中混凝土网格节点与钢筋网格节点的复杂拓扑约束关系，极大缩短建模时间）；(5) 定义两类网格的材料本构模型，素混凝土网格采用经典Holmquist-Johnson-Cook （HJC）模型，钢筋混凝土网格采用本文中提出的基于修正混合物理论的钢筋混凝土复合材料本构模型。

图1 为网格尺寸为25 mm 六面体单元的有限元模型，网格数为1.08×10。图1(a)为整体框架结构有限元模型，1～3 层为混凝土部分，4～6 层为钢筋骨架线，其中图1(b)为局部放大图。图1(c)展示了构件级混凝土网格和钢筋混凝土网格。对于大型混凝土结构而言，如采用实体单元—结构单元的传统分离式数值模型建模，以25 mm 的网格为例，需处理千万级钢筋混凝土节点与混凝土节点的约束关系，前处理在网格剖分过程中将耗费大量的时间人力成本，不利于工程实用化。基于改进型分离式数值模型，可快速实现工程级有限元模拟，且网格分辨率达到钢筋尺寸。

图1 钢筋混凝土框架结构分离式有限元模型Fig. 1 Improved discrete finite element model of reinforced concrete frame structures

1.1 修正型混合物理论模型

根据经典混合物理论，每个组分变形满足应变协调假定，即：

式中：ε和 ε（ i=1,2···）分别为应变和第个组分的应变张量。应变协调假定在复合材料理论中具有一定的局限性。对于由混凝土基体和钢筋组成的复合材料而言，当钢筋和混凝土间的界面发生变形或者滑移时，界面滑移将导致复合材料应变场不连续，并影响混凝土与钢筋间的应力传递，从而导致钢筋中应力下降。通过考虑界面变形或者滑移，应变协调方程（式(1)）可写成：

式中：下标f 和m 分别代表钢筋和混凝土，ε为界面变形或者滑移的应变张量。不用显式构造界面，而将界面的粘结滑移效应引入到钢筋应力应变关系中。等效钢筋弹塑性本构关系可表达为：

基于修正混合物理论，兼顾混凝土基体和具有粘结滑移效应的钢筋力学行为，复合材料本构方程可表达为：

1.2钢筋等效模型

Belarbi 等根据实验结果获取了嵌入在混凝土中钢筋的实际屈服强度满足如下关系式：

式中：ρ为有效配筋率，为混凝土在开裂应变为8×10时的拉伸断裂应力。

在粘结滑移模型中，钢筋混凝土粘结破坏有3 种形式：

(1)混凝土劈裂破坏，当由混凝土保护层或者箍筋提供的环向约束较弱时，钢筋周边环向拉伸应力将造成混凝土劈裂失效。

(2)拉拔式破坏，当环向约束较强时，避免了混凝土劈裂破坏模式；当界面粘结强度较小时候，界面剪应力容易达到粘结强度，从而诱发拉拔式破坏模式。

(3)钢筋断裂失效，当界面粘结长度较大、粘结强度较高时，钢筋应力优先达到最大拉伸应力，从而诱发钢筋断裂失效。

图2钢筋尺度模型Fig.2Schematicofthemodelatsteelbarscale

在爆炸荷载钢筋混凝土构件毁伤破坏过程中，通常表现为钢筋屈服、混凝土剥落等失效形式。即钢筋屈服时，界面剪应力满足上升段本构关系，根据图2 中力平衡关系给出滑移位移计算公式：

钢筋等效弹性模量可进一步由式(5)、(6)和(7)计算。

钢筋等效模型中，除了修正屈服强度、弹性模量等参数，同时需对塑性硬化模量进行修正。等效塑性模量为：

其中钢筋应力可由图2 中力平衡关系和界面最大剪应力给出，δ=δ。

因此，通过式(5)、(6)、(8)，可获得考虑界面粘结滑移的钢筋等效应力-应变关系，其具体表达式为：

图3 弹性模量缩放因子 /Ef 与粘结区长度l 的关系Fig. 3 Relationship between the scaling factor of the elastic modulus ( /Ef ) and the length of the bonding zone (l) under different physical parameters

图4 塑性硬化模量缩放因子/Hf 与粘结区长度l 的关系Fig. 4 Relationship between the scaling factor of the plastic hardening modulus ( /Hf ) and the length of the bonding zone (l) under different physical parameters

从图3 和图4 中可以得出：钢筋等效弹性、塑性硬化模量与粘结区长度呈正相关性，且均收敛到稳定值。塑性模量的收敛速度大于弹性模量的收敛速度。钢筋直径越粗、屈服强度越高、混凝土强度越高，则钢筋等效模量越大，规律趋势与Dehestani 等的结果一致。当粘结区长度较小时，钢筋混凝土粘结容易发生拉拔式破坏，界面变形受钢筋直径、屈服强度、混凝土强度等因素影响较大。此时考虑界面粘结滑移的等效模量对钢筋直径、屈服强度、混凝土强度等因素较为敏感。当粘结区长度较大时，容易诱发钢筋断裂的破坏模式，其等效模量主要与钢筋模型相关。

1.3 混凝土HJC 本构模型

混凝土材料在爆炸荷载作用下将处于大应变、高围压及高应变率的状态。在众多混凝土动态本构模型中，HJC 本构模型因其简明合理的描述和计算程序的适应性，在混凝土爆炸问题中获得了广泛应用。因此，复合材料中混凝土基体部分采用HJC 本构模型，其应力-应变关系为：

式中： σ为混凝土应力，为混凝土抗压强度，为内聚力参数，为损伤变量，为压力硬化系数，为压力，为压力硬化指数，为应变率参数。

1.4 钢筋混凝土复合材料本构模型

结合钢筋等效模型及混凝土HJC 模型，复合材料本构方程（式(4)）可表达为：

2 层次化验证与确认

2.1 构件级精度验证

改进型数值模型的功能模块已在冲击波结构毁伤程序中实现，该程序是一款大规模流固耦合并行程序，能够高效高精度地模拟爆炸冲击类问题，其中流场分析模块是基于JASMIN 框架研发，结构动力学响应分析模块是基于JAUMIN 框架研发，数值模拟结果采用大规模可视分析平台TeraVAP 进行后处理分析。为验证数值模型在钢筋混凝土结构爆炸毁伤方面的计算精度以及粘结滑移模型的有效性，选取了轴压作用下钢筋混凝土柱及板爆炸实验。

2.1.1 轴压作用柱爆炸实验

采用18.2 kg 岩石乳化炸药（等效11.1 kg TNT），爆高为1.5 m。钢筋混凝土柱长2500 mm，截面为边长200 mm 的正方形。主筋选用4 根直径为20 mm 的HRB335 级钢筋，对称布置，配筋率为3.14%，钢筋抗拉屈服强度为466.7 MPa，极限强度为601.3 MPa，弹性模量为200 GPa，塑性硬化模量为1.12 GPa。箍筋选用直径为8 mm HRB235 级钢筋，间距为150 mm，配筋率为0.34%，钢筋抗拉屈服强度为483.5 MPa，极限强度为582.4 MPa，弹性模量为200 GPa，塑性硬化模量为0.94 GPa。混凝土为C40 级商品混凝土，保护层厚度为20 mm。爆炸实验前在柱端部施加轴向压力12 MPa，爆炸过程中柱构件下表面通过LVDT 位移计测量位移时程，位移计布置于柱跨中位置。

基于冲击波结构毁伤程序分别建立考虑钢筋与混凝土间粘结滑移的改进型分离式数值模型和不考虑粘结滑移效应的分离式数值模型，网格尺寸均为10 mm。素混凝土网格采用HJC 混凝土本构模型，模型参数为：初始密度为2.4 g/cm，剪切模量为14.4 GPa，泊松比为0.20，=50.16 MPa,=0.28，=1.84，=0.84，=0.007。对于钢筋混凝土网格，在考虑界面粘结滑移的数值模型中，参照Dehestani 等的研究，界面粘结区长度=300 mm，通过纵筋、箍筋等效材料参数（表1），及混凝土参数，再结合式(12)，给出钢筋混凝土网格中复合材料模型参数。在不考虑粘结滑移数值模型中，纵筋、箍筋材料参数为原始参数。

表1 钢筋等效材料参数Table 1 Equivalent material parameters of rebar

对于爆炸荷载的模拟，根据药量大小和目标距离，采用ConWep 爆炸荷载模型，通过结合反射超压和入射超压，并把入射角θ 的影响考虑进来，计算得到目标上压力，直接作用在结构上，计算模型为：

在近场爆炸实验中，通过预备性爆炸实验，标定ConWep 模型有效性。图5 为爆心距离1.5 m 位置处反射超压时程曲线的ConWep 模拟结果与实验结果对比情况。从图5 可以看出，通过修正ConWep 模型中药量参数，获取与实验测量结果相当的超压时程曲线模拟结果，因此该模型可作为有效的构件荷载输入。

图5 反射超压时程曲线对比Fig. 5 Comparison of the reflected overpressure-time history curves

图6 为柱跨中竖向位移及破坏形貌的模拟结果与实验结果的对比。从对比结果可以看出：考虑钢筋与混凝土间粘结滑移效应时，数值模拟的峰值位移和残余位移与实验值误差分别为6.28%和8.54%，试验件损伤破坏形貌与模拟的损伤破坏形貌吻合较好。相比于不考虑粘结滑移效应有限元模型，改进型分离式模型的模拟结果与实验结果更接近，这说明在爆炸荷载作用下钢筋混凝土结构动力学响应有限元模拟中，应该考虑钢筋与混凝土间粘结滑移。

图6 柱竖向位移时程曲线及破坏形貌对比Fig. 6 Comparison of the vertical displacement-time history curves and failure morphology of columns

2.1.2 板爆炸实验

实验中药量为10 kg TNT，爆高2.0 m。钢筋混凝土板几何尺寸为1200 mm×500 mm×50 mm，钢筋型号为直径8 mm 的HRB235 级钢筋，采用双层双向配筋，主方向间距100 mm，次方向按构造配筋间距200 mm，保护层厚度10 mm，配筋率为0.96%，钢筋参数详见2.1.1 节。混凝土为C60 级商品混凝土，保护层厚度为10 mm。

类似地，基于冲击波结构毁伤程序分别建立考虑粘结滑移的改进型分离式数值模型和不考虑粘结滑移效应的分离式数值模型，网格尺寸均为10 mm。素混凝土单元采用HJC 混凝土本构模型，混凝土强度为39.2 MPa，模型参数参照2.1.1 节内容。钢筋等效材料参数为：等效强度为420.3 MPa，等效弹性模量为170.1 GPa，等效塑性硬化模量657.0 MPa，钢筋混凝土单元中钢筋体积分数为50%。

图7 为C30 板跨中竖向位移的模拟结果与实验结果的对比。从对比结果可以看出：考虑粘结滑移效应时，数值模拟的残余位移与实验结果更为接近。进一步说明在爆炸荷载作用下钢筋混凝土结构动力学响应有限元模拟中，钢筋与混凝土间粘结滑移的必要性。

图7 板竖向位移时程曲线对比Fig. 7 Comparison of the vertical displacement-time history curves of the slab

研究有限元网格尺寸对爆炸荷载下板变形的影响，选取10、5、2.5 mm 等3 种网格尺寸，其中对应网格数分别为30000、240000、1920000。图8 为3 种网格尺寸的数值模拟结果：数值模拟结果与网格尺寸的大小相关，如10 mm 与5 mm 的模拟结果略有差异；随着网格尺寸逐渐减小，跨中位移时程曲线差异减少并逐步收敛（5 mm 和2.5 mm 模拟结果）；进一步的细化网格尺寸可有限地提高计算精度，但是计算成本成倍增加，因此10 mm 的网格尺寸能够兼顾计算精度及计算成本。

图8 网格尺寸对板跨中竖向位移时程曲线的影响Fig. 8 Effect of mesh size on the vertical displacement-time history curves of the slab

2.2 结构级精度验证

通过结构级爆炸实验，进一步验证改进型分离式数值模型计算精度。Woodson 等和Baylot 等开展了1/4 缩尺的两层框架结构外场爆炸实验。为了避免爆炸冲击波与框架结构强耦合相互作用的干扰，即爆炸冲击波在结构内外表面复杂入射、反射、绕射等现象，选取了文献[26-27]中含砖墙的框架结构爆炸实验用于验证模型。该框架结构宽3.2 m，深1.53 m，高2 m。框架结构由6 根柱支撑，结构中柱截面尺寸为89 mm×89 mm，柱高为2 m，背面中柱截面尺寸为76 mm×76 mm，其余的柱截面尺寸为152 mm×152 mm。底层和顶层楼板厚度均为41 mm，层高0.91 m。立柱高于顶层楼板0.178 m。边梁位于中柱和边柱间，截面尺寸为81 mm×54 mm。四块填充墙均由13×15 块砖墙堆砌而成，砖墙尺寸为宽99 mm，厚48 mm，高48 mm。实验中采用7.1 kg C4 半球形炸药（1 kg C4 炸药等效为1.19 kg TNT），爆点位置离地面高度为305 mm，距离底层中柱1.07 m。爆炸前，第2 层中柱顶部施加1.7 t 重物，模拟2.1 MPa 预应力荷载。图9 为改进型分离式数值模型，网格尺寸为5 mm，总计814 万单元，其中暴露在外部分为钢筋混凝土网格。图9 中插图为爆炸前实验设置。素混凝土网格采用HJC 混凝土本构模型，其基本参数参照2.1.1 节。钢筋混凝土网格采用复合材料模型，参数由混凝土强度和钢筋基本参数给出。在炸药起爆前，在目标柱顶部施加2.1 MPa 预应力分析，当结构所有节点速度低于0.1 m/s 时，便认为结构达到新的静力平衡。随后开展爆炸模拟，总模拟时间为0.3 s。

图9 改进型分离式数值有限元模型及实验（单位：m）Fig. 9 Improved discrete finite element model and experiment (Unit: m)

图10 为框架结构破坏形貌的模拟结果与实验结果的定性对比，其中图10(a)为爆炸荷载作用刚结束时框架结构损伤云图，图10(b)为结构最终状态损伤破坏云图，图10(c)为钢筋骨架的塑性区云图，图10(d)为实验结果。模拟结果表明：(1) 迎爆面砖墙均失效；(2) 低层中柱损伤区域最早出现在柱上下两端，随后发展至柱跨中背爆面处，由此可见该柱失效模式为剪切为主，弯曲为辅的耦合失效（图10(a)和10(b)）；(3) 低层中柱钢筋骨架的塑性区域主要集中在柱两端部及跨中背面位置；(4) 低层中柱失效后，失去支撑的边梁通过其抗弯承载力提供框架的抗倒塌能力，在“梁机制”作用下框架结构未发生垮塌。