钢筋混凝土框架结构动力超自由度单元

Dynamic Super-Degree-of-Freedom Element of Reinforced Concrete Frame Structures

摘要: 钢筋混凝土框架结构的超自由度单元被拓展到钢筋混凝土框架结构的动力分析中，拓展的单元称之为动力超自由度单元。根据ABAQUS提供的UEL用户程序接口，动力超自由度单元被应用到ABAQUS动力分析中。通过一个三层三跨的钢筋混凝土框架结构在小震、中震和强震三种地震波作用下振动台试验结果的对比分析，本文开发的用户程序和动力超自由度单元的可靠性得到验证。结果表明程序能够很好地描述地震导致的钢筋混凝土框架结构强度和刚度的退化，说明本文开发的程序和动力超自由度单元可以广泛地用于钢筋混凝土框架结构的动力分析。

Abstract: The super-degree-of-freedom element, which is an element of beam-column joint of reinforced concrete frame structures, is developed to be suitable for dynamic analysis of reinforced concrete frame structures. The element is called dynamic-super-degree-of-freedom element. According to UEL interface provided by ABAQUS platform, the element is applied to ABAQUS dynamic analysis. By means of comparison with the results of a shaking table test of 3-story-3-bay reinforced concrete frame building under minor-moderate, moderate-severe and severe ground motions, reliabilities of the user subroutine and the element were validated. The results show that the subroutine can well reflect strength and stiffness degradation of reinforced concrete frame structures produced by earthquake loading. The subroutine and the element are widely applicable to refined dynamic simulation of concrete frame structures.

1. 引言

钢筋混凝土结构的梁柱节点是制约结构非线性分析走向实际应用的关键因素，自然也是国内外土木工程数值仿真研究者一直关注的研究热点 [1] [2] [3] [4] [5] 。自文 [1] 作者将其节点的有限元分析程序置于OpenSees平台后，极大地改善了混凝土结构非线性分析方法走向实际应用在此之前的窘迫状态，同时也吸引了更多研究者对于梁柱节点的研究兴趣。继文 [1] 后，文 [2] 也将其节点单元开发到OpenSees平台之上。根据节点变形和受力机理，文 [1] 和 [2] 将节点从功能上区分为节点核心区剪切块的剪切变形、梁柱与核心区交界面的剪切变形和梁柱受力钢筋与核心区混凝土的粘结滑移三大功能组成部分，并分别用三组共13根一维弹簧来描述该三大组成部分的材料模型。由于该类梁柱节点单元的组成方式不同于一般单元，因此，该类梁柱节点单元被称为“超级单元”。对于“超级单元”核心区剪切块材料模型怎样建立，不同的学者有不同的主张。因此，超级单元中核心区的剪切材料模型怎样建立，对于超级单元的使用者来说是一个难点。为了避免核心区材料模型建立的困难，文 [6] [7] [8] 开发了以平面4节点和平面8节点钢筋混凝土单元代替核心区剪切弹簧的钢筋混凝土框架梁柱节点单元，并通过ABAQUS平台二次开发了用户应用程序，实现了通过ABAQUS可以进行钢筋混凝土框架结构的非线性滞回分析。本文在文 [8] 的钢筋混凝土梁柱节点超自由度单元基础上，将超自由度单元发展到钢筋混凝土框架结构动力分析，并在ABAQUS平台上开发了用户程序。适用于动力分析的超自由度单元，不妨称为动力超自由度单元。

2. 超自由度单元简述

超自由度单元的力学特性，由8根滑移弹簧、4根剪切弹簧和4节点平面单元共同构成。超自由度单元由4个外节点和4个内节点组成，其中外节点的自由度为3个，与普通梁单元自由度一致，见图1(a)；而内节点自由度与普通平面应力单元的自由度一致，内节点自由度见图1(b)。

根据超自由度单元的自由度构成原则，该单元有4个节点，单元节点的位置与4个外节点一致，而每个单元节点的自由度组成则是由外节点自由度与内节点自由度混合而成。具体地说，单元的总自由度由公式(1)确定：

$U={\left\{\begin{array}{cccc}{D}_1& D_2& D_3& D_4\end{array}\right\}}^T$ (1)

式中： $D_i$ ( $i=1-4$ )为第i节点的自由度。以 $D_1$ 为例，具体的自由度由公式(2)组成，其它节点依次类推。

$D_1={\left\{\begin{array}{ccccc}{\Delta }_1& {\Delta }_2& {\Delta }_3& u_1& v_1\end{array}\right\}}^T$ (2)

Figure 1. (a) Node and degree of freedom in the element; (b) degree of freedom in the internal node

图1. (a) 单元的节点与自由度；(b) 内节点自由度

式中： ${\Delta }_1$ 、 ${\Delta }_2$ 和 ${\Delta }_3$ 为外节点1的自由度，对应单元节点1的水平、竖向和转角位移；而 $u_1$ 和 $v_1$ 为内节点1的自由度，对应内节点1的水平与竖向位移。因单元的节点与外节点一致，因此，内节点位移对于单元的节点1而言，没有实质位移与之对应。在ABAQUS具体实现时，其命令如下：

*USER ELEMENT, NODES=4, TYPE=U4, PROPERTIES=36, ……

1, 2, 6, 14, 15

也就是说，上述命令中的自由度1、2、6与实质位移 ${\Delta }_1$ 、 ${\Delta }_2$ 和 ${\Delta }_3$ 对应，即：自由度1和2与节点的水平和竖向位移对应，自由度6与转动位移对应；而上述命令中的自由度14和15不与本节点的实质位移对应，它们实际上是内节点的水平与竖向位移 $u_1$ 和 $v_1$ 。

3. 单元动力分析的质量矩阵

从文 [1] [2] 单元在OpenSees中的源程序来看，均采用零质量矩阵的形式，即单元给定的质量全部为0。这在OpenSees系统内应该是可以用零质量矩阵进行动力分析的，但在ABAQUS系统内是不行的，因为ABAQUS系统无论静力分析还是动力分析，每个增量步系统首先要检测单元的质量矩阵，如果单元的质量矩阵为0矩阵，则系统停止运算。超自由度单元是基于4节点平面单元组合钢筋与节点核心区混凝土黏结滑移特性而成，同时考虑与梁柱单元的连接。由于平面梁柱单元每个节点有水平、竖向的惯性力和转动惯性力，因此，超自由度单元的每个节点也应该包含该三个方向的惯性力。但超自由度单元的每个节点的5个自由度中的后2个自由度(见式(2))，对于该节点而言是没有实质物理量与之对应的，因此，对应的2个自由度质量应该为0。基于这些因素，本文假定，超自由度单元的总质量矩阵为：

$M=\left[\begin{array}{cccc}{M}_1& & & 0\\ & M_2& & \\ & & M_3& \\ 0& & & M_4\end{array}\right]$ (3)

式中： $M_i$ ( $i=1-4$ )为第i节点的质量矩阵，它由下面的式(4)组成：

$M_i=\left[\begin{array}{ccccc}m& & & & 0\\ & m& & & \\ & & I& & \\ & & & \alpha & \\ 0& & & & \alpha \end{array}\right]$ (4)

式中：m取单元总质量的四分之一；I取绕单元质心的总转动惯量的八分之一； $\alpha$ 取0或者取很小的值(如：0.001m)。

4. ABAQUS用户程序设计

ABAQUS用户程序UEL，就是为每一个计算步提供单元所需的AMATRX与RHS。其中，AMATRX在静力分析时主要是当前运算步的单元刚度矩阵，而在动力分析时，则根据系统的LFLAGS数组的指令，可以是刚度矩阵、质量和阻尼矩阵。而RHS则是单元的节点抗力向量，动力分析时要考虑惯性力以及上一计算步的节点抗力等。

根据动力分析过程的需要，ABAQUS平台控制过程会根据LFLAGS(3)的控制参数确定需要计算的单元参量，具体如下：LFLAGS(3) = 1，定义RHS和当前步刚度矩阵AMATRX；LFLAGS(3) = 2，定义当前步刚度矩阵AMATRX；LFLAGS(3) = 3，定义当前步阻尼矩阵AMATRX；LFLAGS(3) = 4，定义当前步质量矩阵AMATRX；LFLAGS(3) = 5，定义RHS；LFLAGS(3) = 6，定义RHS和当前步质量矩阵AMATRX。因此，将超自由度单元设计成可以用于钢筋混凝土框架结构的动力分析，就是要根据LFLAGS(3)上述数值定义相应的AMATRX和RHS。具体公式如下：

当LFLAGS(3) = 1时，需要计算：

$AMATRX=\left(1+\alpha \right)K+\frac{M}{\beta {\left(\Delta t\right)}^2}$ (5)

和

$RHS=F_J\left(t+\Delta t\right)-\left(1+\alpha \right)F\left(t+\Delta t\right)+\alpha F\left(t\right)$ (6)

式(5)和(6)中的 $\alpha$ 和 $\beta$ 是ABAQUS动力算法的参数， $\Delta t$ 为时间增量，( $t+\Delta t$ )为当前时间步，t为前一时间步，K为刚度矩阵，M为质量矩阵，F为单元的节点抗力， $F_J$ 为惯性力。

当LFLAGS(3) = 2时：AMATRX = K。

当LFLAGS(3) = 3时：AMATRX为阻尼矩阵，本文方法不考虑单元的阻尼。

当LFLAGS(3) = 4时：AMATRX = M。

当LFLAGS(3) = 5时：RHS在式(6)基础上在加上 $\alpha$ F( $t-\Delta t$ )。

当LFLAGS(3) = 6时：RHS = F，AMRTRX = M。

根据式(4)的假定，LFLAGS(3) = 1和5，计算RHS时，式(4)中 $\alpha$ 取0；而LFLAGS(3) = 4和6时，AMATRX实际是质量矩阵，为了质量矩阵非奇异， $\alpha$ 应该取很小的值，本文取0.001m。

5. 验证

5.1. 试验简介

该框架试验是BRACCI等人1992年在纽约州立大学进行的一系列框架振动台综合试验 [9] [10] ，按照实际结构1/3缩小比例建造的两榀三层三跨钢筋混凝土框架模型，测试在地震下的结构响应。振动台输入的地震波数据，在该大学开发的钢筋混凝土结构非线性动力分析有限元软件IDARC的examples中。该试验一共测试了三个不同地震强度下的结构响应，地震波加速度峰值分别为0.05 G、0.2 G及0.3 G (G为重力加速度)，代表了小震、中震和强震情况。

试验模型宽3048 mm，两榀框架柱之间间距1828.8 mm，一榀框架的柱距也是1828.8 mm，每层柱净高1066.8 mm，每层梁高度152.4 mm，梁宽76.2 mm，板厚50.8 mm，柱子为方形，边长101.6 mm。混凝土平均强度为26.9 Mpa，弹性模量为26,890 Mpa，峰值应力对应的应变为0.0021。边柱梁端上部配2根4号钢筋，底部配2根5号钢筋。中柱梁端上部配3根4号钢筋，底部配2根4号钢筋。柱子配4根4号钢筋。4号钢筋直径5.7 mm，屈服强度468.84 Mpa，弹性模量214,000 Mpa。5号钢筋直径6.4 mm，屈服强度262.01 Mpa，弹性模量同4号。

试验模型为了模拟实际结构的动力特性，还根据三分之一的比例关系在模型上设置了配重质量。模型的每楼层质量包括：纵向与横向的梁4.23 kN，板21.48 kN，混凝土配重条53.376 kN，铅砖33.94 kN，以及木支撑与电缆合计4.89 kN，平均分配到8个柱顶的质量约为1.47 ton。

5.2. ABAQUS分析模型

本文选择一榀框架建立模型，框架柱和梁的分析单元选B21梁单元，梁柱节点选本文的动力超自由度单元，楼层的配重选mass单元，每层mass单元放置于每层柱的顶部。其中，柱子单元截面是矩形，梁单元截面是T型，T型翼缘宽度选600 mm (相当于2榀框架间距的三分之一)。每层柱子和每跨梁都等距分8个单元，扣除分析模型自身的梁与板质量后，mass单元的质量约为1.4 ton。钢筋采用rebar的形式，钢筋采用文 [11] 滞回模型，4号钢筋的3个输入参数为：2.14E5, 468.84, 0.0011；5号钢筋的输入参数为：2.14E5, 262.01, 0.003。混凝土模型采用文 [12] 的F1-CON和F2-CON，混凝土模型的6个输入参数为：26.9, 0.2, 0.017, 0.01, 26890, 0.0021。

首先，对上述ABAQUS分析模型进行频率分析，以便确定分析模型的混凝土材料模型和Reighlay阻尼。混凝土模型选F1-CON时，得到系统的前3个频率分别为：1.7391 Hz，5.6272 Hz和9.7105 Hz。选F2-CON时，得到系统前3个频率分别为：2.2027 Hz，6.4789 Hz和10.036 Hz。由文 [9] 可知，试验测得框架的第一频率是1.76 Hz，因此，F1-CON混凝土模型更适合本分析模型。

系统的阻尼比变化范围大约在0.02~0.05，根据上述频率可以估计Rayleigh阻尼质量影响系数约在0.33到0.83之间变化，本文选择0.33。

5.3. 数值分析结果

图2是地震波峰值为0.05 G情况下框架结构顶层与基础相对位移ABAQUS分析结果与试验结果的对比。

从图2可以看出，无论最大位移值还是波形，本文分析结果都能很好地反映试验结果。这个结果表明，动力超自由度单元在动力分析中也有着不俗的表现，而且，选取的Rayleigh阻尼的质量影响系数为0.33，也可以基本确定是正确的。

图3为框架底层柱的剪力ABAQUS分析结果与试验的对比，可以看出，本文计算结果与试验结果基本都是吻合的。

图2和图3中，从大约9 s到13 s之间顶层相对位移与底层剪力的ABAQUS分析结果与试验结果在波形上都存在不吻合的情况，说明ABAQUS分析模型的刚度变化在该时间段内与实际试验模型的刚度变化不一致，而在13 s后，分析结果的波形再次与试验结果一致，说明此时分析模型的刚度再次能够准确表达实际模型的刚度。图2和图3的结果反映了小震情况下动力超自由度单元的动力特性是比较可靠的。

图4和图5为地震波峰值为0.2 G的中震情况下，顶层相对位移和底层柱剪力的ABAQUS分析结果

Figure 2. Third floor displacement time history for 0.05 G

图2. 0.05 G第三层位移时程

Figure 3. First floor shear force time history for 0.05 G

图3. 0.05 G第一层剪力时程

Figure 4. Third floor displacement time history for 0.2 G

图4. 0.2 G第三层位移时程

Figure 5. First floor shear force time history for 0.2 G

图5. 0.2 G第一层剪力时程

Figure 6. Third floor displacement time history for 0.3 G

图6. 0.3 G第三层位移时程

Figure 7. First floor shear force time history for 0.3 G

图7. 0.3 G第一层剪力时程

与试验结果的对比。

从图4和图5可以看出，除从大约18 s到25 s分析结果与试验结果相差较大外，其它时间段内分析结果与试验结果无论是波形还是数值都基本吻合。这说明动力超自由度单元在中震情况下分析精度也是可以接受的。

图6和图7为地震波峰值为0.3 G的强震情况下，顶层相对位移和底层柱剪力的ABAQUS分析结果与试验结果的对比。

从图6和图7可以看出，除从大约6 s到9 s的时间段外，无论在波形还是数值上，分析结果与试验结果也基本是吻合的。而6 s到9 s可能正是模型的刚度变化时，分析模型的刚度变化与实际试验模型刚度没有同步实现的表现，而在9 s后，两者的刚度趋于一致。因此，可以说动力超自由度单元在强震情况下的分析精度也是令人满意的。

9. 结论

本文将钢筋混凝土梁柱节点的超自由度单元拓展到用于钢筋混凝土框架结构的动力分析中，同时在通用商业平台ABAQUS上开发出用户定义的UEL。通过对一个两榀三层三跨的钢筋混凝土框架结构在地震波加速度峰值为0.05 G、0.2 G和0.3 G分别对应小震、中震和强震情况的动力分析，经与试验数据比较，验证了动力超自由度单元模型的有效性，说明本文开发的动力超自由度单元在小震、中震和强震情况下分析钢筋混凝土框架结构的动力响应都是有效的。

基金项目

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