（1. 湖南大学 土木工程学院，湖南 长沙 410082；2. 长沙理工大学 土木与建筑工程学院， 湖南 长沙 410076） 摘 要：针对已有的钢筋混凝土梁单元非线性分析模型采用较多的假定和近似从而导致计算量增大及计算精度下降的问题，基于共旋坐标法建立了考虑材料和几何非线性的任意截面钢筋混凝土梁的数值分析模型.首先利用虚功原理计算共旋坐标系下完全粘结钢筋混凝土梁考虑材料非线性的切线刚度矩阵，再通过结构坐标系与共旋坐标系下节点力之间及节点位移之间的总量关系及微分导出的增量关系，最终获得钢筋混凝土梁在结构坐标系中考虑几何与材料双重非线性的切线刚度矩阵.算例结果表明，本文算法可减少计算量、不累积误差、精度高.
关键词：钢筋混凝土； 梁单元；共旋坐标法；双非线性；微分法；切线刚度矩阵
中图分类号：TU323.3 文献标识码：ACorotational Procedure for the Binonlinear Analysis
of Reinforced Concrete Beam Element
DENG Jihua1，2， SHAO Xudong1
（1. College of Civil Engineering， Hunan Univ， Changsha， Hunan 410082， China； 2. School of Civil
Engineering and Architecture， Changsha Univ of Science and Technology， Changsha， Hunan 410076， China） Abstract：Multiple assumptions and approximations in nonlinearity analysis models of existing reinforced concrete beam element result in low calculation efficiency and low calculation accuracy. In this paper， based on corotational procedure， a numerical model for a given section considering material and geometrical nonlinear analysis of reinforced concrete beam element was developed. Firstly， by means of virtual work， a tangent stiffness matrix for the material nonlinearity of perfectlybonded reinforced concrete beam element was derived in corotational coordinate system. Then， by building total and incremental relationships derived from differential equations of nodal displacements and forces between global coordinate system and corotational coordinate system， respectively， tangent stiffness in global coordinate system reinforced concrete beam element was developed by considering geometric and material nonlinearity. A comparison between the results in this paper and those from existing references has demonstrated that the algorithm developed is highly efficient and accurate with many advantages， such as noncumulative calculation errors and reduction in computation.
Key words：reinforced concrete； beam element；corotational procedure；binonlinear；variational method；tangent stiffness matrix
--------------------------------------------------------------------------------
在目前已有的钢筋混凝土杆系结构考虑几何与材料双重非线性的各种分析模型中，分层梁模型由于具有不受截面形状限制、不同梁层可采用不同材质、能模拟已开裂的梁结构等优点而应用较广＼[1-3＼].但目前的分层梁模型存在以下缺点：1）在分层梁模型中，将钢筋划分成截面的一层，认为截面每一层的应变沿梁轴向均匀分布.因此，为保证计算精度，单元长度需划分得很短，这在非线性计算中是非常不利的；2）将钢筋沿单元轴向理想化为平行于梁轴线的直线段，而在实际钢筋混凝土梁中，由于受力或构造的需要，钢筋并不总是完全平行于梁轴线；3）材料非线性分析中，一般通过假定截面形心处的应变和曲率得到各层的应变，由材料本构关系得到应力，再由截面平衡得到计算内力，将其与实际内力比较以确定单元刚度，这一过程往往需要多次反复迭代，甚至由于混凝土和钢筋的应力
--------------------------------------------------------------------------------
应变曲线为分段曲线而不收敛＼[4＼]；4）几何非线性一般通过应变计算中计入高阶项，基于总体拉格朗日法（T.L法）或修正拉格朗日法（U.L法）来考虑，这不仅导致单元切线刚度矩阵异常复杂，还具有T.L法和U.L法本身所固有的缺点，如T.L法在非线性程度较高时所得计算结果精度很差＼[5＼]，以及U.L法为保证计算精度，需将荷载步取得较小，导致计算量显著增加和误差累积的问题[6].
为解决上述问题，须对现有非线性算法进行改进.在几何非线性分析方法中，共旋坐标法相对于T.L法和U.L法而言，具有列式简单、力学概念明确以及计算精度高的特点，因而成为研究的热点＼[7-8＼].本文在上述文献的基础上，首先利用共旋坐标系下应变与扣除刚体位移后的变形呈线性关系的特点计算出完全粘结钢筋混凝土梁考虑材料非线性的切线刚度矩阵，再基于静力平衡通过微分获得钢筋混凝土梁在结构坐标系中考虑几何与材料双非线性的切线刚度矩阵，多个算例表明本文算法是正确的.
湖南大学学报（自然科学版） 2013年
第8期 邓继华等：钢筋混凝土梁单元双非线性分析的共旋坐标法
1 钢筋混凝土梁元切线刚度矩阵
图1所示为钢筋混凝土平面梁元，基于实际情况，可设钢筋在梁单元内为直线.
图2示意了在结构坐标系XY中初始时刻和计算时刻t钢筋混凝土梁元的几何参数及位移的即时变量；设钢筋混凝土梁元的共旋坐标系为xy，该坐标系是随单元变形而转动的，它始终以节点i为原点，以节点i到j的连线方向为x轴，由x轴逆时针转90°为y轴.
图1 钢筋混凝土梁元
Fig.1 Reinforced concrete beam element
图2 变形前后梁单元
Fig.2 Beam element before and after deformation
1.1 共旋坐标系下混凝土梁元切线刚度矩阵
设初始时刻单元节点在结构坐标系里的坐标为（xi，yi）和（xj，yj），在计算时刻t结构坐标系中的位移向量为d=uiviθiujvj θjT，在共旋坐标系中的位移向量为dl=u′i （1）
混凝土梁单元在共旋坐标系中的节点位移为：
（2）
式中：0l和tl分别为变形前、后梁单元长度.
设单元任一截面形心轴处的应变及截面曲率为ε0，φ，在共旋坐标系下单元的应变
--------------------------------------------------------------------------------
位移关系只需考虑线性项，有：
（3）
式中：B0为线性应变矩阵.
采用沿梁高分层的方法，对计算截面任意梁层k，基于平截面假定，由ε0和φ可求出该梁层的应变值εk为：
（4）
式中：yk为梁层k到形心轴的距离.
对式（4）微分，有：
（5）
通过混凝土的应力
--------------------------------------------------------------------------------
应变关系由εk可求得梁层k的应力σk和切线模量ETk，截面力P=（ncmc）由截面平衡条件求得为：
（6）
式中： Ak为梁层k的面积；n为截面分层总数.
对式（6）微分，有：
（7）
将式（5）代入式（7）并写成矩阵形式有：
（8）
（9）
由虚功原理可建立单元的平衡方程组有：
（10）
式中：fc=fxicfyicmicfxjcfyjcmjcT为共旋坐标系下混凝土单元等效节点力列阵.
将式（9）代入式（10）可得：
（11）
对式（11）微分，并考虑式（8）有：
（12）
式中：kTc为共旋坐标系下考虑材料非线性的混凝土元切线刚度矩阵.
1.2 共旋坐标系下钢筋单元切线刚度矩阵
参考图1所示，设钢筋上距离左节点所在截面为x的任意一点A在竖向距离中性轴为d（x），且x与d（x）的值在变形过程中始终不变.对于钢筋只需考虑轴向拉压应变，与前面混凝土应变计算类似，在共旋坐标系下钢筋的应变计算也只考虑线性项，由应变旋转公式＼[9＼]知钢筋的应变为：
（13）
式中：εc，A为混凝土在A点的应变.
对式（13）微分，并联立式（3）有：
（14）
式中：B0，s为钢筋的应变矩阵值，只需将钢筋具体位置的值代入B0就很容易求得.
通过钢筋的应力
--------------------------------------------------------------------------------
应变关系由εs可求得该点的应力σs及切线模量ETs，钢筋的轴向力Ps为：
（15）
式中：As为钢筋的面积.
对式（15）微分，并联立式（14）有：
（16）
由虚功原理可建立共旋坐标系下钢筋单元的平衡方程有：
（17）
式中：fs=fxisfyismisfxjsfyjsmjsT为共旋坐标系下钢筋元等效节点力列阵.
将式（14）代入式（17）可得：
（18）
对式（18）微分，并联立式（16）有：
（19）
式中：kTs=cos 3θ∫0l0BT0，sAsETsB0，sdx即为共旋坐标系下考虑材料非线性的钢筋元切线刚度矩阵.
1.3 结构坐标系下钢筋混凝土梁切线刚度矩阵
从式（2）知，由于u′i，v′i，v′j恒为0，对式（2）的后3项微分，不难得到共旋坐标系下位移微分δdl用结构坐标系下位移微分δd表达的形式：
（20）
设钢筋混凝土梁单元在结构坐标系下的节点力向量F=FxiFyiMiFxjFyjMjT，由静力平衡可知：
（21）
式中：t为坐标转换矩阵.
将式（21）两边微分可得：
（22）
对矩阵t微分，为表述方便，将δtT（fc+fs）改写成：
（23）
为得到tT（δfc+δfs）用δd表达的形式，联立（12）（19）（20）有：
（24）
联立式（22）（23）和（24），可得到结构坐标系下钢筋混凝土梁单元考虑几何与材料双重非线性的单元切线刚度矩阵为：
（25）
2 非线性分析流程
1）根据上一计算时刻单元i，j节点在结构坐标系下的总位移向量d，由式（2）求出共旋坐标系中的位移向量dl，由式（11）和（12）求出混凝土梁元在共旋坐标系下的切线刚度矩阵kTc及等效杆端力fc；再基于式（18）和（19）求出钢筋梁元在共旋坐标系下的切线刚度矩阵kTs及等效杆端力fs.
2）通过式（25）得到钢筋混凝土单元在结构坐标系下的切线刚度矩阵KT，基于式（21）得到结构坐标系下的等效杆端力F.
3）重复1）至2）的步骤，生成结构切线刚度矩阵∑KT和等效杆端力合力∑F.
4）计算不平衡力ΔR=L-∑F，其中L为到计算时刻t施加的总外荷载的等效节点力.
5）求解结构方程∑K·Δd=ΔR，得到节点位移增量Δd，将其叠加到总位移向量d中.
6）收敛条件判断，如收敛，转到t+Δt时刻计算，如不收敛，返回1）进行下一次迭代.
3 材料的应力
--------------------------------------------------------------------------------
应变关系
受压区混凝土采用的应力
--------------------------------------------------------------------------------
应变关系为：
（26）
式中：fc为混凝土轴心抗压强度；ε0为与fc对应的应变，且有ε0=0.002；εu为极限压应变，且有εu=0.003 8.
受拉区混凝土应力
--------------------------------------------------------------------------------
应变关系为：
普通钢筋受拉和受压时都采用理想弹塑性的应力
--------------------------------------------------------------------------------
应变关系.
4 算例分析
例1 如图3所示的肘式框架，两端嵌固，由William完成其试验分析工作，由于该结构的几何非线性十分突出，同时又具有相应的试验结果，因而成为众多研究者，如Wood， Papadrakakis， Yang， Chan， Meek等人，检验各自所建立理论的有效性与准确性的标准算例.杆件为一十分细长矩形等截面直杆，截面宽19.1 mm，高6.71 mm，弹性模量为71 000 MPa.分析时将每根杆件划分成10个单元，采用位移增量法求解，跨中截面位移与荷载的关系如图4所示，可看出与William的试验结果是比较吻合的.
图3 William肘式框架（单位：mm）
Fig.3 William’s toggle frame（Unit：mm）
挠度Δ/mm
图4 荷载
--------------------------------------------------------------------------------
挠度曲线
Fig.4 Loaddeflection curves
例2 如图5所示均布荷载作用下的悬链线无铰拱模型＼[3＼]，拱的跨度为L=4 m，矢跨比为f/L=0.2，拱轴系数m=2.24.横截面为5.5 cm（宽）×20 cm（高）的等截面矩形，截面顶和底各配有5根φ4的钢筋，钢筋形心到截面上下缘的距离取为5 mm，材料性质见文献＼[10＼].
图5 模型拱的节点划分
Fig.5 Node partition of arch model
表1列出了试验结果、本文及文献＼[1＼]的计算结果，图6示出了本文计算的4，5，6号关键节点在线性、几何非线性、材料非线性及几何与材料双非线性下的荷载
--------------------------------------------------------------------------------
挠度曲线.从表1可看出，本文计算结果是正确的，由图6可看出，本模型几何及材料非线性效应均比较明显，必须考虑两者的耦合作用.
表1 试验与计算极限荷载值
Tab.1 Loadcarrying capacity of test and calculation
kN/m
试验值＼[3＼] 材料非线性 双非线性
本文 文献＼[1＼] 本文 文献＼[1＼]
63.42 90.86 90.63 67.67 60.84
挠度/mm
（a） 4号节点
挠度/mm
（b） 5号节点
挠度/mm
（c） 6号节点
图6 荷载
--------------------------------------------------------------------------------
挠度曲线
Fig.6 Loaddeflection curves
5 结 论
本文基于共旋坐标系下应变与位移的线性关系，利用虚功原理直接由截面切线刚度矩阵通过积分导出单元考虑材料非线性的切线刚度矩阵，再通过结构坐标系与共旋坐标系下节点力之间及节点位移之间的总量关系及微分导出的增量关系，获得钢筋混凝土梁在结构坐标系中考虑几何与材料双重非线性的切线刚度矩阵，避免了已有文献须反复迭代求解单元切线刚度矩阵的缺点，考虑了应变沿梁轴向的变化；同时，从以上推导过程也知钢筋方向与梁轴线方向可斜交，算例表明本文算法具有较高精度，适用于钢筋混凝土结构的几何与材料非线性分析.诚然，在混凝土开裂及接近破坏时，混凝土与钢筋完全粘结的假定不再成立，这是以后研究中应仔细考虑的问题.
参考文献
[1] 刘磊，张光卿，袁长卿. 杆系结构的高精度非线性分析＼[J＼].土木工程学报，2006 ，39（1）：25-28.
LIU Lei， ZHANG Guangqing， YUAN Changqing. The high precision nonlinear analysis for frame structures＼[J＼]. China Civil Engineering Journal， 2006，39（1）：25-28. （In Chinese）
[2] 奉龙成，罗小华. 钢筋混凝土拱面内极限承载能力的非线性分析＼[J＼].土木工程学报，2002，35（3）：20-24.
FENG Longcheng， LUO Xiaohua. Nonlinear analysis for ultimate bearing capacity of arch bridge＼[J＼].China Civil Engineering Journal，2002，35（3）：20-24. （In Chinese）
[3] 陈克济. 钢筋混凝土拱桥面内极限承载力的非线性分析＼[J＼]. 桥梁建设， 1983（1）： 24-36.
CHEN Keji. Nonlinear analysis on inplane ultimate loadcarrying capacity of reinforced concrete arch bridge＼[J＼]. Bridge Construction， 1983（1）：24-36.（In Chinese）
[4] 曾永革，李传习. 钢筋混凝土梁单元截面刚度求解方法＼[J＼]. 铁道科学与工程学报， 2009，6（1）：62-67.
ZENG Yongge， LI Chuanxi. Calculating methods for section rigidity of reinforced concrete beam element＼[J＼]. Journal of Railway Science and Engineering， 2009，6（1）：62-67. （In Chinese）
[5] 吕和祥， 朱菊芬， 马莉颖. 大转动梁的几何非线性分析讨论＼[J＼]. 计算结构力学及其应用， 1995，12（4）：485-490.
LV Hexiang， ZHU Jufen， MA Liying. Discussion of analyzing of geometric nonlinear beams with large rotations＼[J＼]. Chinese Journal of Computational Mechanics， 1995， 12（4）： 485-490. （In Chinese）
[6] 邓继华，邵旭东. 基于U.L列式的带刚臂平面梁元非线性分析＼[J＼]. 湖南大学学报：自然科学版， 2012，39（5）：8-12.
DENG Jihua， SHAO Xudong. Nonlinear analysis of plane beam element with rigid arms based on U.L formulation＼[J＼]. Journal of Hunan University： Natural Sciences， 2012，39（5）： 8-12. （In Chinese）
[7] CRISFIELD M A， MOITA G F. A corotational formulation for 2D continua including incompatible modes＼[J＼]. International Journal for Numerical Methods in Engineering， 1996， 39： 2619-2633.
[8] 蔡松柏，沈蒲生.大转动平面梁有限元分析的共旋坐标法＼[J＼].工程力学，2006，23（S1）： 69-72.
CAI Songbai， SHEN Pusheng. Corotational procedure for finite element analysis of plane beam element of large rotational displacement＼[J＼].Engineering Mechanics， 2006，23（S1）： 69-72. （In Chinese）
[9] 康清梁. 钢筋混凝土有限元分析＼[M＼]. 北京： 中国水利水电出版社，1995.
KANG Qingliang. Reinforced concrete finite element analysis ＼[M＼]. Beijing： China Water & Power Press， 1995. （In Chinese）
[10] 刘磊. 大跨度混凝土桥梁的双非线性分析＼[D＼]. 北京：北方交通大学土木建筑工程学院，2000： 70-71.
LIU Lei. The coupling nonlinear analysis of long span concrete bridges＼[D＼]. Beijing： College of Civil Engineering and Architecture， Northern Jiaotong University， 2000： 70-71. （In Chinese）
