王 强，常 凯，侯康康，路 炯，刘 琳

(沈阳建筑大学土木工程学院，辽宁沈阳 110168)

0 引 言

在强震作用下，复杂高层建筑结构会进入塑性状态，因此进行结构的动力弹塑性时程分析对于评估结构的抗震安全性具有重要意义。整体结构分析时，采用梁单元能够以较少的计算量合理描述梁柱构件的力学性态。梁单元模型通常有塑性铰和纤维模型2类，纤维梁单元模型直接从材料单轴应力-应变本构关系和截面构造配筋出发，相对于塑性铰模型，其能够更准确地描述梁柱构件的轴力-双向弯曲耦合效应，在工程界已得到了广泛认可[1-2]。

通用有限元软件ABAQUS拥有强大的非线性分析能力，在分析大型复杂结构时得到了广泛应用。然而目前该软件缺乏用于三维纤维梁单元的混凝土本构模型，无法采用梁单元描述钢筋混凝土梁柱杆件的非线性力学行为，因此采用该软件进行建筑结构的动力弹塑性时程分析存在困难。为此，已有学者通过ABAQUS的用户子程序接口，基于Scott-Kent-Park的混凝土本构模型对用于三维梁单元的混凝土单轴本构模型进行二次开发[3-4]。

2008年汶川地震后，中国修订了《混凝上结构设计规范》(GB 50010—2010，以下简称《规范》)，给出了考虑损伤的混凝土本构关系。采用该本构关系进行结构的弹塑性分析对于中国工程界具有重要的现实意义。目前，已有部分学者基于该本构模型进行了二次开发，如胡晓斌等[5]、齐虎等[6]，但前者并未给出受拉段的滞回规则且模拟结果与试验结果差距较大，而后者缺少实例数据对比。因此，本文对《规范》附录C混凝土本构关系进行了深入研究，首先补充了受拉加卸载滞回规则与考虑箍筋对混凝土的约束作用，并修正了受压损伤恢复的问题，随后利用ABAQUS用户子程序接口UMAT，VUMAT进行二次开发，编写了用于显、隐式分析的梁单元混凝土单轴本构模型子程序，用于钢筋混凝土梁柱构件在多维受力状态下的弹塑性分析，满足结构动力弹塑性分析的需求。

1 纤维梁单元

基于材料单轴本构关系的纤维模型是将梁柱构件沿纵向划分为若干子段，再沿构件横截面划分成纤维束。ABAQUS中的B31纤维梁单元构型如图1所示。每个纤维(以截面积分点表征)只考虑它的轴向本构关系，且可定义不同的本构关系。构件横截面变形符合平截面假定。对截面纤维的当前状态积分就可以得到截面的双向抗弯刚度、双向抵抗矩以及轴力，进而沿构件长进行积分，可以得到精确的构件单元刚度矩阵。纤维模型可以自然、简单地描述构件的双向弯曲-轴力耦合效应[7-8]。

2 材料的本构模型

纤维模型的准确度在很大程度上依赖其使用的材料本构模型与材料参数的合理性。为准确描述梁柱构件中钢筋与混凝土的材料性能，本文采用2个共节点的梁单元来分别模拟构件中的纵向钢筋与混凝土材料，且纵向钢筋通过等面积、等位置原则代换成箱形截面单元，从而分别定义混凝土和钢筋的本构模型与材料参数。

2.1 钢筋的本构模型

大部分学者在模拟钢筋混凝土构件受力性能时，钢筋本构通常采用ABAQUS自带本构，但该本构对于包兴格效应考虑不足。因此，本文钢筋本构采用自编单轴本构，可较为充分地考虑钢材的包兴格效应以及屈服、硬化与软化特性，其中本构的骨架曲线采用Esmaeily-Xiao模型[9]，加卸载曲线采用Légeron模型[10]，具体图形如图2所示。图2中，σ为应力，ε为应变。该本构模型的准确性验证和模型的参数取值方法可参考文献[11]。

2.2 混凝土的本构模型

2.2.1 受压区骨架曲线与滞回规则

《规范》给出了素混凝土的受压区骨架曲线及滞回规则，如图3所示，其中，fc,r为混凝土抗压强度代表值，εc,r为混凝土受压峰值应力对应的应变，σun，εun分别为受压混凝土骨架曲线开始卸载时的应力和应变，Er为受压混凝土卸载/再加载的变形模量，εz为受压混凝土卸载至零应力点时的残余应变。骨架曲线具体公式参见文献[12]。

(1)受压区骨架曲线上的应力-应变关系

受压区骨架曲线上的应力-应变关系采用混凝土单轴受压损伤演化参数dc表征，即

σ=(1-dc)E0ε

(1)

(2)

(3)

(4)

(5)

式中：E0为素混凝土的弹性模量；αc为单轴受压应力-应变曲线下降段参数。

文献[13]指出，钢筋混凝土梁柱构件中由于箍筋对核心区混凝土的约束作用，使得混凝土的受压应力-应变曲线随着约束效应的增大由明显的陡峰向平缓饱满的曲线过渡。因此在构件滞回性能模拟时有必要考虑箍筋的约束作用。本文参照文献[11]，采用约束指标λt来考虑钢筋混凝土梁柱构件中箍筋对混凝土的约束作用，对《规范》中的混凝土应力-应变曲线方程进行修正。λt计算公式如下

(6)

式中：fyv为箍筋屈服强度；ρv为混凝土梁柱构件中箍筋的体积配箍率。

此时，约束混凝土的抗压强度fcc、峰值应变εpc及曲线下降段参数αdc为

fcc=(1+0.5λt)fc,r

(7)

εpc=(1+2.5λt)εc,r

(8)

αdc=(1-1.75λt)αc

(9)

将fcc，εpc，αdc代替式(2)～(5)中相应的fc,r，εc,r，αc，则可得梁柱构件中考虑箍筋约束作用的混凝土应力-应变骨架曲线。

考虑箍筋约束作用的混凝土与非约束混凝土的应力-应变曲线对比如图4所示。由图4可以看出，当柱构件中混凝土考虑箍筋约束作用时，混凝土抗压强度、峰值应变等都有所提高，尤其是曲线下降段形状趋于平缓。

(2)受压区滞回规则的修正

《规范》已给出了重复荷载作用下受压混凝土卸载与再加载路径方程，即

σ=Er(ε-εz)

(10)

(11)

(12)

(13)

式中：εca为附加应变。

已有研究表明，对于混凝土材料，重复加卸载条件下的材料损伤是不可恢复的，即受压混凝土卸载/再加载的变形模量单调降低[13-14]。当卸载应变较大时，采用式(11)～(13)确定的受压混凝土卸载/再加载的变形模量Er，会出现随卸载应变增大而增大的现象，即混凝土受压损伤恢复。

由图5可以看出，采用《规范》公式，在卸载应变超过A点时出现了变形模量Er增大的情况。A点应变为2.87εc，即为εc/(εc+εun)=0.09εun/εc时对应的卸载应变εun。因此，本文对式(13)进行了修正，采用式(14)计算附加应变εca，即

(14)

计算结果表明，按公式(14)确定的附加应变得到的变形模量Er单调下降，能够合理表征循环加载下混凝土材料的受压损伤不可恢复特性。

2.2.2 受拉区骨架曲线与滞回规则

《规范》中给出了受拉区骨架曲线，如图6所示，其中，ft,r，εt,r分别为受拉区应力、应变，具体公式参见文献[12]，但是《规范》中并没有给出明确的受拉区滞回规则。目前王瑾等[15-16]在OpenSEES中采用Concrete02模型在模拟混凝土构件滞回性能的分析中取得了较好的结果。本文参照OpenSEES中Concrete02模型，采用原点指向规则，对《规范》混凝土本构补充混凝土受拉区滞回规则。

混凝土在拉-压与压-拉往复加载条件下，本文采用拉区应力-应变曲线随受压残余应变点迁移的滞回规则。若混凝土初始受力即为受拉，则沿受拉区骨架曲线加载。若受拉卸载，卸载曲线由受拉卸载点B指向初始受拉点A，如图7所示。随后受压加载，则沿受压区骨架曲线加载，加载至C点时卸载至受压残余应变点D点。继续受拉加载时其路径为DE(DE∥BA)，至F点进入受拉卸载段时，依据原点指向规则，其路线为指向受压残余应变点D的FD直线。受压再加载按DC路线回到受压曲线上的C点。

3 用户自定义材料子程序的实现

ABAQUS软件提供了基于隐式算法的UMAT材料子程序接口和基于显式算法的VUMAT材料子程序接口，分别应用于隐式求解器(Standard)和显式求解器(Explicit)[17]。

采用前述的混凝土与钢筋本构关系，笔者基于UMAT接口进行二次开发，编制了材料子程序，并嵌入ABAQUS/Standard模块中。由于隐式求解必须引入雅可比矩阵，而雅可比矩阵和N+1时刻的应力、应变以及状态变量是实时更新的，这导致材料非线性发展较为强烈时会出现雅可比矩阵奇异而产生隐式求解的数值不收敛问题。因此，本文又基于VUMAT接口，编制了用于显式分析的材料子程序。在进行梁柱构件乃至整体结构的非线性分析时，对于竖向加载工况下结构非线性发展较弱的工况，采用隐式求解方法；对于水平加载等结构非线性发展比较强烈的工况则采用显式分析方法。通过显、隐式结合来完成结构在不同荷载工况下的接力计算。

4 算例验证

为验证本文所采用混凝土本构的合理性，分别对文献[18]中钢筋混凝土柱试件边柱A及文献[19]中钢筋混凝土柱试件TP74，TP77进行数值模拟。

4.1 试件边柱A

本文首先模拟了文献[18]中轴压-单轴受弯钢筋混凝土柱试件边柱A在低周往复荷载下的滞回性能。

4.1.1 试验概况

试件构造和配筋如图8所示，加载图示见图9，加载规则如图10所示，材料参数取值见表1。

4.1.2 结果分析

试验与模拟结果对比如图11所示。可以看出，当不考虑钢筋混凝土柱中箍筋对混凝土的约束作用时，本文采用的混凝土本构模型已能够较好地模拟混凝土柱的滞回性能，滞回曲线饱满程度与试验结果接近，但承载力略低。当考虑混凝土柱中箍筋对混凝土的约束作用时，计算所得承载能力有所提高，模拟结果与试验结果更为吻合，表明本文所采用的混凝土本构模型较为合理。

表1 材料参数Tab.1 Material Parameters

图12为在考虑箍筋约束情况下修正附加应变前后计算得出的滞回曲线对比。可以看出，2条曲线基本重合，证明本文对公式(13)的修正是可行的。

图13为试件边柱A在加载幅值10，25，45 mm时的单圈滞回曲线对比。由图13可知：加载幅值为

10 mm时，滞回环峰值点与实测曲线峰值点较为接近，相对误差达到2%，但计算曲线的单圈耗能能力(单圈滞回环包含的面积)略高于实测曲线的耗能能力；加载幅值达到25 mm时，计算曲线的滞回环峰值点略高于实测曲线，计算曲线的耗能能力略高于实测曲线的耗能能力，相对误差低于16%。当加载幅值达到45 mm时，滞回环峰值点较实测曲线峰值点高约20%，且计算曲线较“瘦”，其单圈耗能能力略低于实测曲线，相对误差约为4%。这可能与模拟分析未考虑试件的剪切非线性、轴压力波动、黏结滑移等因素有关。此外，计算曲线能够反映实测曲线的捏缩现象。

图14为计算所得试件边柱A底部钢筋、混凝土单元截面中的角部纤维应力-应变曲线。由图14可知，在单轴往复加载过程中，混凝土纤维受压已进入骨架曲线下降段，钢筋纤维受拉已进入骨架曲线软化段，且应力-应变曲线符合本文所构建的混凝土本构模型的滞回规则。

4.2 试件TP74，TP77

本文对文献[19]中的试件TP74和TP77进行模拟分析。构件单元及截面划分与试件边柱A相同。试件及试验概况详见文献[19]。

4.2.1 单轴加载试件TP74

图15为TP74实测滞回曲线和采用本文构建的材料本构计算所得的滞回曲线对比，考虑箍筋约束和不考虑箍筋约束2种情况。可以看出，计算所得滞回曲线饱满程度均与试验结果较为接近，且考虑箍筋约束后计算所得承载力有所提高，与实测曲线更为吻合，证明了本文所构建的混凝土本构模型的合理性。

图16为在考虑箍筋约束情况下计算得出的附加应变修正前后滞回曲线对比。通过对比发现，修正附加应变前后最大承载力一致，但修正后的承载力并无明显退化现象，且滞回曲线与图15(a)的实测滞回曲线更为吻合，这表明了本文修正附加应变的正确性。

4.2.2 双轴加载试件TP77

图17为试验与数值模拟所得滞回曲线，由图17可知，计算所得的承载力、强度退化程度、滞回曲线饱满程度均与试验结果较为接近。

此外，在试验曲线与计算曲线中，在保持方向1位移恒定时，在方向2加载均会使方向1的承载力出现明显的下降趋势，反之亦然。此现象表明，本文提出的模型能够合理地反映梁柱构件的双向弯曲耦合性能。

5 结 语

(1)本文采用原点指向模型与拉区应力-应变曲线随受压残余应变点迁移法，补充完善了《规范》混凝土本构模型的混凝土受拉段加卸载、拉压受力状态转换的滞回规则。

(2)本文采用基于体积配箍率建立的约束指标λt考虑钢筋混凝土梁柱构件中箍筋对混凝土的约束作用，对《规范》混凝土受压应力-应变骨架曲线方程进行修正，有效提高了钢筋混凝土梁柱构件滞回性能模拟结果的准确性。

(3)本文对《规范》混凝土本构模型中的附加应变算法进行了修正，所得的受压混凝土卸载/再加载变形模量Er能够合理地考虑循环加载下混凝土受压损伤不可恢复特性。

(4)本文构建的混凝土本构模型及相应的用户材料子程序能够细致地描述梁柱构件受力过程中截面的内力变化情况，准确地模拟钢筋混凝土柱在低周往复荷载作用下的滞回性能，能进一步用于大型复杂结构的弹塑性反应分析。