基于MSC.MARC的纤维梁子程序开发与验证

2018-11-09閤东东苗启松

土木工程与管理学报 2018年5期

閤东东，张纯，苗启松

(1. 北京工业大学建筑工程学院，北京 100124； 2. 北京市建筑设计研究院有限公司，北京 100045；3. 碧桂园集团北京区域公司，北京 100070)

动力弹塑性时程分析是罕遇地震作用下结构变形验算的重要方法，也是结构性能化设计的重要手段。在钢筋混凝土结构的动力弹塑性时程分析中，应用最为广泛的通用有限元软件主要有ABAQUS和MSC.MARC(以下简称MARC)，由于这两款软件并非结构工程专用分析软件，采用其进行动力弹塑性时程分析需要进行纤维梁子程序二次开发。

已有部分学者对MARC在结构动力弹塑性分析中的应用进行了研究。陆新征等[1,2]基于MARC平台开发的THUFIBER程序可以准确模拟地震作用下结构的动力响应。Han等[3]将MARC用于组合结构的动力弹塑性分析。汪训流等[4]编制的混凝土和钢筋单轴本构分别采用Légeron-Paultre模型[5]和修正的Légeron模型[6]，对钢筋混凝土柱的受力状态进行了预测，与试验结果吻合良好。聂建国等[7]也在MARC的基础上，开发了纤维梁单元模型，用于分析钢-混凝土组合结构。柳国环等[8]也在MARC的基础上，开发了可用于模拟钢材和混凝土纤维梁单元的用户子程序。Lu等[9]运用MARC成功的模拟了632 m高的上海中心大厦在强震作用下的连续倒塌过程。

为采用MARC对钢筋混凝土结构进行动力弹塑性时程分析和倒塌分析奠定基础，本文开发了基于刚度法的纤维梁单元子程序，并与2011年中国建筑学会抗震防灾分会建筑结构抗倒塌专业委员会进行了钢筋混凝土框架结构拟静力倒塌试验进行了对比，验证了程序的正确性和可靠性。

1 MARC中的纤维梁单元

纤维模型的主要原理是将杆件截面离散化为若干个纤维，并假定构件在各受力阶段，整个截面满足平截面假定，各纤维受力分别采用各自的单轴应变-应力曲线来描述，并根据相应纤维材料的单轴应变-应力关系来计算截面内力与变形的非线性关系。纤维梁单元通过截面方向和长度方向上两次积分得到单元的刚度矩阵，通过材料级别的应变-应力关系实现精细化的有限元建模。纤维梁模型包含基于刚度法和柔度法的两种模型，刚度法采用基于位移的插值函数，并忽略构件横截面剪力对构件正截面受力的影响。目前在通用有限元软件中，都采用基于刚度法的模型。

MARC提供的纤维梁子程序接口UBEAM，包含有98号可以考虑剪切变形的两节点单积分点线性Timoshenco梁单元，和52号不考虑剪切变形的两节点三积分点Euler-Bernoulli梁单元。在接口UBEAM中，只需要给出单元截面层次的刚度矩阵，即可实现粱式结构的复杂非线性行为模拟。因此，在UBEAM中，可以采用纤维梁模型通过截面方向上的积分得到截面的刚度矩阵，也可以提供特征截面的塑性铰模型。相比较而言，前者更容易实现，且更容易考虑轴力与双向弯矩的耦合作用。下面主要对基于刚度法的纤维梁单元截面内力与变形的关系进行介绍。

图1所示为矩形构件截面示意图，截面被离散为混凝土纤维和钢筋纤维。对于第i根纤维，其形心坐标为(yi,zi)，面积为Ai。首先进行截面分析，根据平截面假定，截面各点的轴向应变εi可以表达为：

εi=ε0+κyyi+κzzi

(1)

式中：ε0为形心处的轴向应变，以受拉为正；κy和κz分别为截面主轴的曲率。

图1 矩形截面纤维模型示意

以增量形式表达的混凝土纤维和钢纤维的应力-应变关系为：

Δσi=EiΔεi

(2)

式中：Ei为第i根纤维的切线模量。

截面上所有纤维的合力增量为：

(3)

(4)

式(4)可简化为：

ΔF=DTΔε

(5)

式中:ΔF为与正应力相关的截面内力增量向量；DT为截面切线刚度矩阵；Δε为截面应变增量向量。

动力时程分析时， MARC采用隐式积分算法和Newton-Raphson迭代法求解。迭代过程中，在每个增量步的开始，各单元上一步的截面应变、应力和应变增量由主程序传递到UBEAM中，在增量步的结束时，UBEAM完成截面应力和切线刚度矩阵的更新并返回主程序。在应力更新的过程中，需要借助于一些状态变量，并开辟数组空间记录所有单元的纤维关键历史应力、应变信息，并在每一步迭代收敛时进行更新。

2 混凝土和钢筋材料应变-应力关系

2.1 混凝土材料本构模型

混凝土采用了两种单轴本构模型，分别命名为Concrete01，Concrete02。Concrete01模型采用了OpenSees中的Concrete02模型。Concrete02为GB 50010-2010《混凝土结构设计规范》中提供的模型[10]。两种混凝土材料模型的循环加载曲线如图2所示。

图2 Concrete01与Concrete02混凝土材料本构模型

2.2 钢筋材料本构模型

钢筋单轴本构共包括Steel01和Steel02两种模型，其中Steel01为运用广泛的双线性随动强化模型。Steel02采用OpenSees中的修正Menegotto-Pinto模型。

3 纤维梁子程序测试

3.1 试验概况

2011年，中国建筑学会抗震防灾分会建筑结构抗倒塌专业委员会进行了钢筋混凝土框架结构拟静力倒塌试验，其中对一组钢筋混凝土框架柱进行了滞回试验(图3)，并公布了其试验数据和材料参数[11]。试验过程中，水平加载位置位于距柱底750 mm处，轴压千斤顶球铰距柱底1030 mm。底层中柱C的配筋如图4所示。其中混凝土立方体强度为30.1 MPa。Φ10钢筋屈服强度和极限强度分别为481和745 MPa，弹性模量为265433 MPa，屈服应变为0.002，伸长率为23.6%；Φ8钢筋屈服强度和极限强度分别为582和855 MPa，弹性模量为289850 MPa，屈服应变为0.002，伸长率为28.8%。

图3 钢筋混凝土框架柱滞回试验

图4 底层中柱C配筋/mm

3.2 C柱测试结果

将柱单元沿长度方向上划分为5个单元，柱截面纤维划分示意图如图1所示，其中混凝土部分一共划分了8×8=64个纤维，钢筋根据实际钢筋根数划分了8根纤维。采用上述混凝土和钢筋材料单轴本构模型，基于UBEAM子程序接口，分别编制了52号Euler-Bernoulli梁和98号Timoshenko梁单元的截面本构关系。采用MARC对C柱进行时程分析，持时与拟静力加载试验时间相同。对采用Euler-Bernoulli梁和Timoshenco梁单元，及不同混凝土和钢筋单轴本构时的结果进行了归纳总结，结果定义如下：E和T分别代表52号和98号梁的计算结果，Concrete01Steel01代表混凝土和钢筋分别采用Concrete01和Steel01单轴本构模型时的结果，其余工况依此类推。

中柱C的配筋结果详细数据可见网站。其中位移加载时程和轴压时程分别如图5a，5b所示，加载持时265.427 s。

图5 中柱C往复加载方案

图6所示为52号Euler-Bernoulli梁数值模拟结果与试验结果对比，从图中可以看出，数值模拟结果与试验结果较为一致。分别从图6a，6b可以看出，采用相同的混凝土和不同的钢筋单轴本构模型时，得到的计算结果较为接近，说明混凝土单轴本构对计算结果不会产生严重影响。图7所示为98号Timoshenko梁数值模拟结果与试验结果对比，可见数值模拟结果与试验结果吻合较好，仔细对比图6，7，可以发现当计算参数相同时，52号梁和98号梁得到的计算结果略有区别。数值模拟结果与试验结果较为一致，说明计算结果具有一定的可靠性。

图6 采用Euler-Bernoulli梁时中柱C滞回曲线

图7 采用Timoshenko梁时中柱C滞回曲线

为详细和直观地考察子程序的正确性，利用UBEAM接口输出了混凝土和钢筋纤维应变、应力以及切线刚度时程曲线进行分析。图8所示为52号Euler-Bernoulli梁沿长度方向第一个积分点处，角部第一根混凝土纤维滞回曲线，从图中可以看出，混凝土材料滞回曲线与预期的本构模型一致，当钢筋本构模型不同时，得到的混凝土滞回曲线略有区别。图9所示为角部第一根钢筋纤维应变-应力关系曲线，钢筋材料滞回曲线与采用的模型一致，对比图8可以看出，相对而言，因混凝土材料模型不同引起的钢筋滞回曲线差别较因钢筋材料模型不同引起的混凝土滞回曲线差别大。

图8 Euler-Bernoulli梁角部第一根混凝土纤维滞回曲线

图9 Euler-Bernoulli梁长度方向上第一个积分点滞回曲线

图10输出了沿长度方向第一个积分点处，角部第一根混凝土纤维的切线刚度时程曲线，从图中可以看出不同钢筋材料本构得到的时程曲线较为一致，采用Concrete01模型时局部情况下出现了负刚度，采用Concrete02模型时，可以看出随着构件进入屈服，混凝土材料刚度呈下降趋势。图11所示为角部第一根钢筋纤维的切线刚度时程曲线，从图中可以看出，初始刚度约为2.8×105N/mm2。