三维纤维梁单元增量非线性有限元分析

胡郑州，吴明儿

（同济大学 建筑工程系，上海200092）

根据式（3），得：

杆系结构在极限荷载、强震等作用下，表现高度非线性，若按照线性理论分析，很难满足工程设计要求.Zeris［1］和Spacone［2］提出的纤维梁单元，并成功解决钢筋混凝土柱的材料非线性问题；文献［3］基于柔度法并考虑了双向弯曲与变轴力作用下对空间杆系结构进行了非弹性分析，其分析结果与试验结果相当吻合；文献［4］采用纤维梁单元分析了钢－混凝土组合结构在地震作用下的受力机理和破坏规律；文献［5］采用纤维梁单元在火灾作用下模拟整体结构连续性倒塌过程、机理及破坏规律；文献［6］采用纤维梁单元分析了在考虑方钢管局部屈曲情况下薄壁方钢管混凝土梁柱结构在单调和循环荷载下的性能；文献［7－8］采用考虑几何和材料非线性的纤维梁单元对桥墩在轴力和弯矩作用下进行了二阶分析，亦得出了满意结果.可见纤维梁单元是分析杆系结构高度非线性的最好单元之一.

本文首先根据连续性介质力学和虚位移原理推导了通用的考虑大位移增量非线性有限元增量更新拉格朗日（Updated Lagrangian，UL）列式，同时修正大位移增量矩阵.据此理论推到了三维纤维梁大位移增量有限元UL列式.以上理论编制了分析程序，通过对几个算例分析，证明该单元的精确性、通用性.

1 增量UL非线性列式

基于连续性介质力学和虚位移原理［9］，参考现实时刻t构形坐标系定义得：

式中：kL表示小位移刚度矩阵；kg表示为几何刚度矩阵；kU表示为大位移增量刚度矩阵；Δue表示节点增量位移；N表示单元局部坐标插值函数矩阵；t＋Δttq表示单元节点参考t时刻构形在t＋Δt时刻施加的荷载向量；t＋Δttρ表示单元参考t时刻构形在t＋Δt时刻的密度；t＋Δttf表示单元参考t时刻构形在t＋Δt时刻的体荷载向量；t＋Δt tt表示单元参考t时刻构形在t＋Δt时刻的面荷载向量；ttσ表示t时刻的应力向量；tdv表示t时刻单元积分体积；tda表示t时刻单元积分面积；tV、tA分别表示为t时刻单元的体积、面积.在基于上一时刻构形时下面公式推导过程省略左上角t＋Δt和左下角t，其中kL、kg、kU表达式为

式（2）—（4）中：BL表示线性应变—位移矩阵；BN表示非线性应变—位移矩阵；DT表示材料本构关系矩阵；G表示对形函数偏导矩阵；P表示应力矩阵.式（4）仍然是不对称矩阵，进行下面简化，

于是可使kU矩阵对称化，从而使整个刚度矩阵是对称的，这对单元切线刚度矩阵元素在计算机内存中存储是有利的.Bathe等［9］在推导UL 时将kU项忽略掉，本文保留该项对高度非线性的影响，利用在增量迭代步上一步迭代步的节点增量位移收敛解去求解大位移增量矩阵.

2 三维纤维梁单元刚度矩阵推导

2.1 小位移刚度矩阵

三维纤维梁单元的位移函数采用Euler－Bernoulli梁的位移函数，局部坐标如图1所示，截面上任意一点的轴向应变为

式中：ε11为截面上任意一点的轴向应变；ε0为截面形心轴向应变；y为纤维在局部坐标y－z上y轴坐标；z为纤维在局部坐标y－z上z轴坐标；v″为局部坐标系沿y方向位移两阶偏导；w″为局部坐标系沿z方向位移两阶偏导.

图1 三维纤维梁单元Fig.1 Three dimensional fiber beam element

利用虚位移原理便可推导截面切线刚度矩阵ks［1－2］：

式中：Etan为纤维弹性模量；A为截面面积，dA截面积分面积变量，其中L表达式如下：

得出ks表达式如下：

式中：n为纤维截面中纤维总数；Ei，tan表示第i纤维弹性模量；Ai表示第i个纤维面积；yi，zi定义见图1所示.

根据式（2）得到单元的小位移刚度矩阵：

2.2 几何刚度矩阵和大位移增量刚度矩阵

根据式（3），得：

其中，P，G 表达式如下：

其中：

其中，σ11表示截面形心出轴应力，x表示单元纤维截面形心到单元端部距离.

根据式（5）大位移增量刚度矩阵：

其中：由式（9）—（11）得到单元的局部坐标下参考t时刻构形在t＋Δt时刻构形的切线刚度矩阵，并通过单元坐标转换矩阵，将局部坐标系的切线刚度矩阵转化为整体坐标系下，最后得到U.L.列式的大位移增量方程如下：

3 混凝土和钢筋材料本构

3.1 混凝土本构关系

本文混凝土本构模型采用修正的Kent－Park模型［12］，该模型通过修改素混凝土（无约束混凝土）受压骨架曲线应变软化段斜率来考虑箍筋对混凝土的约束影响，其骨架曲线如图2所示，σc，εc为是混凝土的应力和应变；ε0为混凝土强度峰值应变；K为考虑箍筋对混凝土约束作用所引起混凝土强度增大系数；εC20为约束混凝土强度等级为C20对应的极限应变；ε′C20为素混凝土强度等级为C20 对应的极限应变；fc为混凝土轴心抗压强度，MPa.分段函数加以描述为

式中：εu为混凝土强度极限应变；fcu为混凝土极限强度，MPa；ε0＝0.002K，K＝1＋ρsvfyv／fc，ρsv为从箍筋外边缘计算核心混凝土的体积配箍率；fyv为箍筋屈服强度.其中：Z为混凝土应变软化段斜率；b′为从箍筋外边缘计算的核心混凝土宽度；sv为从箍筋中心算起的箍筋间距.

3.2 钢筋本构关系

本文采用的钢筋本构模型最初是Menegotto等［13］提出的，后经Filippou等［12］修正已考虑等向强化影响的本构模型（图3）.Minegotto等所建议的模型如下，

图2 约束和素混凝土本构关系Fig.2 Constitutive relationship of confined concrete and unconfined

式中，R0，a1，a2均为参数.

Filippou等［12］建议将线性的屈服渐近线进行应力平移，平移大小取决塑性应变的最大值，表达式如下：

式中：σst为钢筋拉应力；fy为钢筋屈服应力；εmax为钢筋最大应变；εmin为钢筋最小应变；εy为钢筋屈服应变.

图3 钢筋本构关系Fig.3 Constitutive relationship of steel

4 大转动计算方法

本文采用Yang［14］更新单元几何坐标的大转动的研究成果，将每步增量位移分解为刚体位移ur和自然变形un，

t时刻构形的轴tx，ty，tz将转化为t＋Δt时刻的构形t＋Δtx，t＋Δty，t＋Δtz，如图4所示，自然变形un参考t＋Δt时刻的构形的表达式可以写成如下，

式中，Δu可以表示为

式中：t＋ΔtL表示t＋Δt时刻构形的单元长度；tL表示在t时刻构形的单元长度.单元节点i转动位移表示为

式中：φ和ni表示如下：

其中，αi、βi、γi表示单元i端截面形心坐标系（αi，βi，γi）向量.

图4 单元刚体位移和自然变形Fig.4 Rigid body displacements and natural deformations

5 算例分析与讨论

5.1 钢筋混凝土桥墩柱push－over分析

钢筋混凝土桥墩柱，墩高2 265mm，矩形截面，沿横桥向宽800mm，顺桥向宽400mm，将桥墩支承的上部结构的重量等效为墩顶的集中力P，大小为393.2kN.本文采用一个纤维梁单元模拟此桥墩柱，每个单元截面由32块核心混凝土纤维、24块保护层混凝土纤维和18根钢筋纤维组成，如图5 所示.保护层和核心混凝土纤维的本构模型采用本文的3.1节的模型，材料特性参数［15］包括混凝土峰值压应力fc、峰值压应变ε0、极限压应力fcu、极限压应变εcu，其值见表1.钢筋纤维采用3.2节的本构模型，初始弹性模量Es0＝200GPa，屈服强度fy＝357 MPa，屈服后刚度硬化系数b＝0.01，式（16）—（17）参数取值为，R0＝20，a1＝18.5，a3＝0.01，a4＝7.

表1 计算结果比较混凝土材料特性Tab.1 Material properties of concrete

对图5的钢筋混凝土桥墩结构，采用本文建议的三维纤维梁单元并对截面的核心混凝土分成32、40、和64 块三种情况进行位移控制的增量pushover分析，得到图6所示的该桥墩的墩底反力—位移关系曲线.

由图6所示，当位移增量增加到0.438m 时，由于单元出现大量的复特征值而首先出现不收敛状况，程序自动停止求解；核心混凝土分别划分32、40和64块计算结果近似，通过试算知一般钢筋混凝土梁柱核心混凝土截面划分到30～40块基本可以满足工程需求；当顶点水平位移到0.012m，保护层混凝土已经退出工作状态，引起桥墩构件的刚度部分退化，但考虑箍筋对核心混凝土约束作用而提高构件的承载力，故在0.012～0.438m 范围内构件内混凝土、钢筋出现应力强化.钢筋混凝土桥墩柱在push－over分析中，墩顶此时发生较大的位移和转角，当达到0.438m 时结构承载力达到极限而发生推覆破坏.特别是在桥墩位移在0.012～0.438m 区间，钢筋混凝土桥墩出现强非线性特性，需要考虑大位移、大转动几何非线性.

5.2 William 肘式框架分析

如图7所示肘式框架（P 为集中力），两端嵌固，杆件为截面为矩形，截面的宽为19.13mm，高6.17 mm，弹性模量为71 018.5 MPa.William［16］对该结构进行试验研究和非线性分析，其分析结果与试验数据相当吻合.William 在进行分析时，考虑了构件的大位移、轴力对弯曲刚度影响及弯曲引起构件缩短等非线性耦合关系.本文在分析该结构时，对每个构件采用6个纤维单元，每个纤维单元截面沿高度方向划分3块纤维，沿宽度方向划分2块纤维进行模拟分析，并与William 对该结构的实验数据进行对比分析，如图8所示，表明本文计算结果在前、后屈曲阶段与William 的实验结果基本一致.

5.3 Remseth空间框架穹顶分析

如图9—10所示Remseth空间框架穹顶结构，支座嵌固，杆件的截面为矩形，其截面尺寸为0.76 m×1.22 m.截面材料特性：弹性模量为20 690 MPa，剪切模量为8 621 MPa.本文计算分析结果与Shi等［17］对该结构进行弹塑性大变形、大转动分析结果和Park等［18］对该结构大位移分析结果进行了对比，如图11所示.本文对该结构每个构件用一个纤维梁单元模拟，单元截面沿宽度划分8块纤维和沿高度划分14块纤维；Shi对该结构分析时：采用能量法和单元任意分布的塑性铰模型，并考虑弯曲—拉伸高度耦合非线性、大变形、大转动作用，用一个单元来模拟一个构件；Park在分析该结构时：构件所有截面沿宽度采用8积分点和沿高度采用14积分点，每个构件采用16个单元进行该结构的大位移、大转动前屈曲和空间框架大变形倒塌分析.图11显示，结果基本一致，表明本文建议的三维纤维梁单元可以精确分析空间框架考虑大位移、大转动的屈曲分析和倒塌分析.图11中横坐标值是根据本文计算结果每3 000个子增量步计算结果来描述的，是根据区域来划分数值而非常规等间距.

6 结语

本文以三维连续体介质力学虚功原理和虚位移原理为基础，推导了大位移增量UL 列式，此增量列式中极大程度上保留大位移增量刚度矩阵kU.然后基于UL增量列式推导出小应变、大位移、大转动三维纤维梁大位移增量非线性有限元UL 列式.该梁单元的切线刚度矩阵考虑了复合材料的非线性、大位移大转动高度几何非线性.通过对几个算例分析，得到非常有用的以下两个结论：

（1）算例表明，本文建议的三维纤维梁单元是一种比较精确的单元，单元截面可以模拟真实的非线性复合材料本构关系；能通过大荷载步增量、相对大的位移增量、每步增量步较少的子迭代步精确求解大位移、大转动问题.

（2）此外本文建议的三维纤维梁可以用于pushover分析；可精确分析结构的前屈曲和后屈服状态荷载与变形高度非线性特性；同时该单元也可用于结构倒塌分析.

［1］ Zeris C A，Mahin S A.Analysis of reinforced concrete beamcolumns under uniaxial excitation［J］.Journal of Structural Engineering，ASCE，1988，114（4）：804.

［2］ Spacone E，Filippou F C，Taucer F F.Fiber beam－column model for nonlinear analysis of R／C frames： part I.Formulation［J］.Earthquake Engineering and Structural Dynamics，1996，25：711.

［3］ 陈滔，黄宗明.基于有限单元柔度法的材料与几何双重非线性空间梁柱单元［J］.计算力学学报，2006，23（5）：524.CHEN Tao，HUANG Zongming.Material and geometrically nonlinear spatial beam－column element based on the finite element flexibility method ［J］. Chinese Journal of Computational Mechanics，2006，23（5）：524.

［4］ 聂建国，陶慕轩.采用纤维梁单元分析钢－混凝土组合结构地震反应的原理［J］.建筑结构学报，2011，32（10）：1.NIE Jianguo，TAO Muxuan.Theory of seismic response analysis of steel－concrete composite structures using fiber beam elements［J］.Journal of Building Structures，2011，32（10）：1.

［5］ 李易，陆新征，叶列平，等.混凝土框架结构火灾连续倒塌数值分析模型［J］.工程力学，2012，29（4）：96.LI Yi，LU Xinzheng，YE Lieping，et al.Numerical models of fire induced progressive collapse analysis for reinforced concrete frame structures［J］.Engineering Mechanics，2012，29（4）：96.

［6］ Zubydan A H，ElSabbagh A I.Monotonic and cyclic behavior of concrete－filled steel－tube beam－columns considering local buckling effect［J］.Thin－Walled Structures，2011，49（4）：465.

［7］ Thai H T，Kim S E.Second－order inelastic analysis of cablestayed bridges［J］.Finite Elements in Analysis and Design，2012，53：48.

［8］ Zupan D，Saje M.The linearized three－dimensional beam theory of naturally curved and twisted beams：the strain vectors formulation ［J］. Computer Methods in Applied Mechanics and Engineering，2006，195（33／36）：4557.

［9］ Bathe K J，Wilson E L.Nonsap—a nonlinear structural analysis program［J］.Nuclear Engineering and Design，1974，29：266.

［10］ Riks E.An incremental approach to the solution of snapping and buckling problems［J］.International Journal of Solids and Structures，1979，15（7）：529.

［11］ Crisfield M A.An arc－length method including line searches and accelerations［J］.International Journal for Numerical Methods in Engineering，1983，19（9）：1269.

［12］ Taucer F F，Spacone E，Filippou F C.A fiber beam－column element for seismic response analysis of reinforced concrete structures［R］.Berkeley：Earthquake Engineering Research Center of University of California，Berkeley，1991.

［13］ Menegotto M，Pinto P E，Method of analysis for cyclically loaded reinforced concrete Plane frames including changes in geometry and non－elastic behavior of elements under combined normal force and bending ［C］／／ Proceedings，IABSE Symposium on Resistance and Ultimate Deformability of Structures Acted on by Well－Defined Repeated Loads.Lisbon：IABSE，1973：15－20.

［14］ YANG Yeongbin.Theory and analysis of nonlinear framed structures［M］.［S.l.］：Prentice－Hall，Inc，1994.

［15］ 禚一，李忠献.基于显式算法的纤维梁柱单元模型［J］.工程力学，2011，28（12）：39.ZHUO Yi，LI Zhongxian.Explicit algorithm－based fiber beam column element model［J］.Engineering Mechanics，2011，28（12）：39.

［16］ William F W.An approach to the non－linear behaviour of the members of a rigid jointed plane framework with finite deflections［J］.The Quarterly Journal of Mechanics ＆Applied Mathematics，1964，17（4）：451.

［17］ Shi G，Atluri S N.Elasto－plastic large deformation analysis of space－frames： a plastic－hinge and stress－based explicit derivation of tangent stiffnesses［J］.International Journal for Numerical Methods in Engineering，1988，26（3）：589.

［18］ Park M S，Lee B C.Geometrically non－linear and elastoplastic three dimensional shear flexible beam element of von－misestype hardending material［J］.International Journal for Numerical Methods in Engineering，1996，99（3）：383.