ABSTRACT

Two fiber-beam-element concrete models are presented. The models are comprised of a tension/compression envelope and unload/reload path. The tension envelope contains an elastic range and a linear soft zone. The compression envelope also includes a linear elastic range followed by the formulation of Eurocode No. EN 1992 and KENT-PARK model, whose peak stress is the concrete cylinder compressive strength for the unconfined concrete and is calculated by the KENT-PARK model for the confined concrete. The unload/reload paths of OpenSees CONCRETE01 and CONCRETE02 are used. According to UMAT interface provided by ABAQUS platform, two user subroutines of concrete, F1-CON and F2-CON, are developed for fiber beam element. By means of comparison with the results of a full scale shaking table test, reliabilities of the user subroutines and modeling methods were validated. The results show that the subroutines can well reflect strength and stiffness degradation of reinforced concrete columns produced by earthquake loading. The subroutines are widely applicable to refined simulation of concrete columns or other flexural members widely.

Keywords:Concrete Structures, Earthquake Loading, Material Model, Numerical Analysis, ABAQUS

全尺寸钢筋混凝土桥墩柱ABAQUS动力分析

方自虎，李向鹏

深圳大学土木工程学院，广东 深圳

收稿日期：2017年10月21日；录用日期：2017年11月5日；发布日期：2017年11月13日

摘 要

本文提出两个纤维梁单元混凝土模型。模型由拉压骨架曲线和加载卸载路径组成，其中，受拉骨架曲线由弹性和线性软化两段组成；受压骨架曲线由线性弹性和随后EN 1992版欧盟混凝土规范建议公式结合KENT-PARK模型组成，其混凝土峰值压应力，对于非约束混凝土采用圆柱体抗压强度，对于约束混凝土采用KENT-PARK模型计算。加载卸载路径分别采用OpenSees的CONCRETE01和CONCRETE02模型路径。依据ABAQUS平台的用户材料接口UMAT，开发了两个适用于纤维梁单元的混凝土模型用户子程序，F1-CON和F2-CON。通过对比全尺寸振动台试验各项结果，本文开发用户子程序与模拟方法的可靠性得到验证，结果表明用户子程序能很好地反映受地震荷载作用的钢筋混凝土柱强度与刚度的退化，模型可以广泛用于精细化模拟钢筋混凝土柱及其它弯曲构件。

关键词 :混凝土结构，地震荷载，材料模型，数值分析，ABAQUS

Copyright © 2017 by authors and Hans Publishers Inc.

This work is licensed under the Creative Commons Attribution International License (CC BY).

http://creativecommons.org/licenses/by/4.0/

1. 引言

钢筋混凝土结构动力非线性性能的研究越来越受到土木工程研究者的重视，自然试验研究的成果也越广泛见诸学术刊物 [1] [2] [3] 。除钢筋混凝土结构的动力试验越来越被研究者重视外，钢筋混凝土结构的动力非线性数值分析也在学术界渐渐受到关注，一些国际竞赛也油然而生，特别是由美国太平洋地震工程研究中心(PEER)于2010年在加州大学进行了一次全尺寸钢筋混凝土桥墩柱动力试验，并将试验相关资料公开放置于网络进行全球征集数值分析结果以来，各种钢筋混凝土结构动力分析的理论研究与材料模型研究也越来越受到学者们的青睐 [4] [5] 。由于通用商业软件，如：ABAQUS、ANSYS、SAP等，具有成熟高效的非线性算法和功能强大的前后处理能力，结构非线性分析的研究者更愿意在这类通用平台上二次开发自己的材料模型或单元来进行结构的非线性分析 [5] 。本文也是在ABAQUS通用平台上开发出用于纤维单元的混凝土材料模型以及钢筋模型，并用这些开发的模型对钢筋混凝土桥墩柱进行非线性动力分析。

2. 混凝土模型简述

本文混凝土受压受拉骨架曲线如图1所示，其中受压曲线分成4段，受拉曲线分成2段。

受压曲线第一段从坐标原点到压应力峰值的40%为结束点，为直线段，反应混凝土的弹性特性。第二段是应变从最大弹性应变 ${\epsilon }_{ce}$ 经压应力峰值对应的峰值应变 ${\epsilon }_{co}$ 再到极限压应变 ${\epsilon }_{cu}$ ，为二次曲线段，反应混凝土的非线性特性。第二段的应力应变关系，采用EN1992欧盟混凝土规范中建议的应力应变关系公式，即下式(1)。

$\frac{\sigma }{f_c}=\frac{k\eta -{\eta }^2}{1+\left(k-2\right)\eta }$ (1)

式中， $\eta =\epsilon /{\epsilon }_{co}$ ， ${\epsilon }_{co}$ 为混凝土压应力峰值对应的应变。可以取实测值，如果没有实测值，取欧盟规

图1. 混凝土模型受压受拉骨架曲线

范建议值，即 $0.7f_c^{0.31}$ 的千分之一，最大值不超过0.0028。而 ${\epsilon }_{cu}$ ，普通混凝土取0.0035，高强混凝土取0.0028。假定 ${\epsilon }_{ce}/{\epsilon }_{co}$ 的比用 ${\eta }_{0.4}$ 表示，则式(1)中的系数k可以由式(2)计算求得：

$k=\frac{2}{3{\eta }_{0.4}}+\frac{{\eta }_{0.4}-0.8}{0.6}$ (2)

混凝土弹性模量如果没有实测值，就由欧盟混凝土规范公式计算，即 $22{\left(0.1{f^{\prime }}_c+0.8\right)}^{0.3}$ (单位为MPa)， ${f^{\prime }}_c$ 为素混凝土棱柱体抗压强度。

第四段是混凝土强度的残存段，来源于修正的Kent-Park混凝土模型 [6] ，该模型认为，当混凝土应变大于最大混凝土压应变 ${\epsilon }_{cm}$ 后，混凝土强度保持为 $\alpha f_c$ ， ${\epsilon }_{cm}$ 的计算公式如式(3)所示。

${\epsilon }_{cm}=0.004+0.9{\rho }_s\frac{f_{yh}}{300}$ (3)

式中， ${\rho }_s$ 为约束混凝土的箍筋体积配筋率， $f_{yh}$ 为箍筋屈服强度。第三段取直线连接第二段与第四段即可。而混凝土峰值应力 $f_c$ 的取值，非约束混凝土取棱柱体抗压强度；对于约束混凝土，根据修正的Kent-Park混凝土模型 [6] ，由式(4)计算得到：

$f_c=\left(1+\frac{{\rho }_s{f}_{yh}}{{f^{\prime }}_c}\right){f^{\prime }}_c$ (4)

本文在ABAQUS中开发的F1-CON和F2-CON模型，在没有实测弹性模量和峰值应变时，输入参数为3或者4个，前3个参数分别对应图1： $f_c$ ， $\alpha$ 和 ${\epsilon }_{cm}$ ，其中F1-CON模型不计混凝土拉应力的影响。第4个参数对应图1的最大混凝土拉应变 ${\epsilon }_{tm}$ ，缺省值为0.01。而 $f_t$ 也用欧盟混凝土规范公式计算，即 $f_t=0.3{\left({f^{\prime }}_c\right)}^{2/3}$ 。如果有实测弹性模量和峰值应变时，输入参数采用6个，前4个同前述一样，后面2个参数分别是弹性模量 $E_c$ 和峰值应变 ${\epsilon }_{co}$ 。两种模型的卸载和重加载路径分别与OpenSees中CONCRETE01和CONCRETE02模型一致，此处不再赘述。

文中计算ABAQUS自带混凝土两种模型所需要的参数时，混凝土拉压应力应变关系均使用图1所示的骨架曲线，这样无论用ABAQUS自带混凝土模型还是本文开发模型，进行数据对比时，就不会因混凝土骨架曲线不同而导致没有可比性。本文的钢筋模型采用作者开发的钢筋混凝土结构的钢筋滞回模型 [7] ，在ABAQUS上开发材料模型的接口程序同文 [8] 。

3. 桥墩柱试验简介

该试验是PEER于加州大学举办的一个钢筋混凝土桥墩柱动力响应试验测试竞赛 [1] [2] [3] (如图2照片所示)，桥墩柱的截面配筋如图3。

图2. 试验照片

图3. 柱截面尺寸、配筋和单元分割

图3中，纵筋为直径36 mm，箍筋为2根直径16 mm，间距15 3mm，箍筋外的混凝土保护层厚度51 mm，具体材料力学性能参数参见文 [1] [2] [3] 。试验时柱顶配有一个质心位于柱顶的混凝土质量块，重2322 kN，即质量236.7吨。该质量块的转动惯量 [1] 为： $J_Z=0.9923\times {10}^9\text{ton}\cdot {\text{mm}}^2/G$ 。

4. 模态分析

首先进行模态分析，确定系统的固有频率。试验时由于与振动台激励方向垂直的方向上，试件位移受到约束，因此，可以认为系统在作平面振动。

在ABAQUS平台上系统的模态分析，整个柱子长7.315米，平分成15个节点，每2个节点间采用线性平面铁摩辛柯梁单元B21。考虑到柱子箍筋对混凝土的约束作用，同时ABAQUS对圆形截面平面梁单元的截面积分最大积分点数只能是9，也即，从轴线到梁的受拉或受压边缘的辛普森积分点只有5点，这样的积分点数很难准确描述如此巨大截面的混凝土非线性特性。因此，本文分析时，每2个节点之间采用3个B21单元来描述柱子截面(如图3所示)。其中最外围的圆环截面表示柱子混凝土保护层截面，环的厚度51 mm，ABAQUS梁单元截面特性选择SECTION = PIPE。图3中间的环，取以钢筋为中心位置的圆环，其厚度为67.5 mm，ABAQUS截面选择同保护层一样。这两个圆环的积分点都选择11点，这样相当于将环分成20等分。

由于外围的环是非约束混凝土，同时，其混凝土的弹性模量和峰值应变均为已知，因此，ABAQUS的输入命令如下：

*MATERIAL, NAME = F1-CONe

*Depvar

10,

*User Material, constants = 6

40.9,.2,0.004,0.01,22900.,0.0026

上述命令中第5行的参数分别为：混凝土峰值应力，残存应力比值，最大压应变，弹性模量和峰值应变。

而第二层环的类似上述第5行则为：

45.2,.2,0.008,0.01,22900.,0.0028

其中，45.2是按照式(4)计算得到，0.0028是0.0026按照同比例放大但最大值为0.0028的结果。

最内层是圆柱体，ABAQUS截面为CIRC，直径982 mm。柱顶的质量块用一个MASS单元，附加转动惯量的ROTARYI单元，进行模态分析和后续动力分析。

ABAQUS模态分析的前2个频率结果为：8.7914 rad/s和42.271 rad/s。后续动力分析时，模型参数与模态分析的一致。

5. 估计系统Rayleigh阻尼

ABAQUS通用数值分析平台处理阻尼，是在输入材料特性时，采用下面的命令来考虑阻尼的：

*MATERIAL，NAME=……

*DAMPING，ALPHA=……，BETA=……

上面的命令行中，ALPHA就是要输入质量影响系数，BETA就是要输入刚度影响系数。

但是，如果材料是用户开发的，即命令行为如下形式时：

*USER MATERIAL，CONSTANTS=……

则上面阻尼系数，只能输入ALPHA而不能输入BETA，即对于用户开发的材料模型，ABAQUS无法考虑BETA阻尼。

当系统的阻尼比为 $\xi$ ，可由下式(5)来估计Rayleigh阻尼的α和β值。

$\alpha =\frac{2{\omega }_1{\omega }_2\xi }{{\omega }_2+{\omega }_1}$ (5a)

$\beta =\frac{2\xi }{{\omega }_2+{\omega }_1}$ (5b)

系统的阻尼比，一般对于钢筋混凝土结构，其变化范围大约在0.02~0.05。因此，按照式(5a)可以估计Rayleigh阻尼质量影响系数从0.291到0.728之间变化。

6. 数值分析结果

该系统受6组工况地震波作用，分别标记为：EQ1到EQ6，每组工况的地震波前面约20秒的白噪音波，最后大约也是20秒以上的白噪音波，其目的应该是消除本组地震波下结构的振动对下一组地震波引起结构振动的影响。从文 [3] 可知，EQ1工况下，柱子纵向钢筋没有屈服，而EQ2工况时，钢筋刚达到或接近屈服，EQ3工况时，钢筋才完全进入屈服阶段。因此，要分析该系统的非线性动力响应，至少要连续分析EQ1到EQ3，因为，EQ1虽然钢筋没有进入屈服阶段，但混凝土无疑已经有非线性表现，其非线性应力应变历史路径对后续的分析都存在影响。

本文以下分析都是针对EQ1到EQ3的，为了节省分析时间，又因重点考虑EQ3的分析，因此，分析的3个工况的时间分别取成3个地震波的前60 s、70 s和50 s。这样取波的考虑如下：虽然去掉了EQ1工况的尾部白噪音波，考虑到EQ1工况时柱子钢筋处于弹性阶段，而EQ2首部的20 s白噪音波仍然存在，EQ1工况的结构振动对EQ2的影响，对分析而言不是主要的，而EQ3工况分析的重点是20 s之后与50 s之前的时间段。

本文首先选择系统的Rayleigh阻尼质量影响系数 $\alpha$ 为0.291进行分析。

图4是本文使用混凝土F1-CON和文 [7] 钢筋模型，EQ3工况柱顶位移ABAQUS分析结果与试验结果的对比。

从图4可以看出，无论最大位移值还是波形，本文分析结果都能很好地反映试验结果。这个结果表明，本文提出的钢筋和混凝土材料模型在动力分析中有着不俗的表现，而且，选取的Rayleigh阻尼的质量影响系数为0.291，也可以基本确定是正确的。也许仅有位移结果还不能充分说明模型的可靠性，下面再看看其它物理量的分析结果与试验结果的对比情况。

从图5的柱顶加速度和图6的柱底弯矩对比图可以看出，本文计算结果与试验结果基本都是吻合的，而图7的柱底剪力与试验结果相差很大。但是，如果仔细对比可以看出，柱底剪力在数值上可以说有很高的吻合度，欠缺的只是波形上的不吻合。

图4. EQ3柱顶位移本文结果与试验结果对比

图5. EQ3柱顶加速度分析结果与试验结果对比

图6. EQ3柱底弯矩分析结果与试验结果对比

图7. EQ3柱底剪力分析结果与试验结果对比

不仅EQ3工况的各项物理指标本文分析值与试验值相当吻合，EQ1的位移结果本文分析结果也是非常吻合的，其结果见图8。

为了明确本文分析结果的精度如何，图9和图10分别给出了2个其他研究者的分析成果，以供比较本文结果。

图8. EQ1工况的柱顶位移分析结果与试验结果对比

图9. BIANCHI分析的EQ1柱顶位移

图10. 文 [4] 分析的EQ1柱顶位移

图9是EQ1工况的柱顶位移直接从文 [9] 截取的图片，文 [9] 是直接引用2011年于希腊召开的结构动力与地震工程计算方法第三次国际会议意大利研究者BIANCHI的研究成果，由于没有看到BIANCHI的原文，仅从引用者提供的图形可知，图中显示的是EQ1工况的10秒到40秒的柱顶位移分析结果。对比本文10秒到40秒的分析结果，可以看出，从EQ1工况的18秒到28秒，这是EQ1工况中主要激励阶段。这阶段本文分析结果，无论是波形还是数值，与试验值可以说是极度吻合，而BIANCHI的分析结果较之试验结果，无疑是偏大的。

图10是文 [4] 分析EQ1工况的柱顶位移，该图也是从文 [4] 中直接截取的，图上半部是其分析全图，下半部为全图中的一个局部放大。从局部放大图形中可以看出，其分析结果在波形上有超前现象，同时最大值小于试验值。本文认为本文的分析结果要优于图10的结果。

再看其它Rayleigh阻尼质量影响系数 $\alpha$ 的对比结果。图3.14给出了 $\alpha$ 为0.18时的分析结果与试验结果的对比。

从图11可以看出，无论是正向位移最大值还是反向位移最大值，质量阻尼影响系数为0.18时的分析结果都是大于试验值的，具体试验值分别为61.8 mm和−55.6 mm，而分析值分别是66.4 mm和−63.1 mm。同时，分析波形在55秒后明显不符合试验振动几乎停滞状态，这种现象正说明是因阻尼偏小引起的。

EQ1工况代表非屈服阶段，EQ3则代表屈服阶段，从上述分析可以看出，本文材料模型无论是结构在屈服阶段还是位于未屈服阶段，都具有很高的分析精度，说明本文开发的材料模型在动力分析中是十分有效的。

7. 与ABAQUS弥散裂缝模型对比

再看看本文材料模型与ABAQUS自带的弥散裂缝混凝土模型和塑性硬化钢筋模型分析结果的对比情况。图12给出了EQ3工况的柱顶位移对比。

从图12可以看出本文材料模型是优于ABAQUS自带的弥散裂缝混凝土模型的。如果再结合运算时间看，本文模型的优势则更加明显。同样的EQ3工况50 s的分析时间，本文模型所用的WALLCLOCK TIME是2536秒，而ABAQUS自带模型的时间是26,006秒，是本文模型分析时间的10倍以上。

8. 与ABAQUS损伤塑性模型对比

本文同样分析了ABAQUS自带的损伤塑性混凝土模型和塑性硬化钢筋模型的三种工况，然而，该模

图11. 质量阻尼系数0.18时EQ1工况柱顶位移

图12. EQ3工况柱顶位移本文与弥散裂缝模型结对比

图13. EQ3工况柱顶位移损伤塑性混凝土模型分析结果

型在EQ2的第62秒后，运算因无法继续而终止。图13给出EQ2工况柱顶位移的该模型分析结果与试验结果的对比。

从图13可以看出，ABAQUS自带损伤塑性混凝土模型的分析结果已经出现漂移，说明已经无法继续分析，因此，EQ3工况也就无法再分析了，自然其可靠性与本文开发的混凝土与钢筋模型是无法比拟的。

综上所述，可知本文开发的混凝土模型与钢筋模型，在动力分析中具有工程所需要的足够精度，说明本文开发的材料模型是有效的。

9. 结论

本文将EN1992欧盟混凝土规范的混凝土压应力应变关系曲线结合Kent-Park混凝土模型 [6] ，形成混凝土的压应力应变关系骨架曲线，拉应力应变关系骨架曲线由两段直线组成，并结合OpenSees中CONCRETE01和CONCRETE02模型的卸载和重加载路径，形成了自己的混凝土纤维单元材料模型F1-CON和F2-CON，同时在通用商业平台ABAQUS上开发出用户定义UMAT。通过对钢筋混凝土桥墩柱在屈曲与未屈曲阶段的动力分析，经与试验数据比较以及对比其他研究者的数值分析成果，验证了本文混凝土材料模型的有效性，说明本文开发的材料模型是有效的。

基金信息

国家自然科学基金(编号：51378313 & 51578336)。

文章引用

方自虎,李向鹏. 全尺寸钢筋混凝土桥墩柱ABAQUS动力分析
Dynamic Analysis of a Full-Scale Bridge Reinforced Concrete Column Using ABAQUS[J]. 土木工程, 2017, 06(06): 564-575. http://dx.doi.org/10.12677/HJCE.2017.66068

参考文献 (References)