损伤梁单元等效节点荷载的计算及应用

郝效强，李 杰

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

建立精细化的结构分析模型是准确模拟结构非线性反应的基础.由于问题的复杂性，在有限元分析中关于非线性问题在很多方面往往采用简化处理的方式，如单元的刚度、质量和等效节点荷载等，其中，单元刚度矩阵和等效节点荷载直接影响结构非线性分析的数值稳定性和计算效率.针对梁单元非线性刚度矩阵已经提出了很多模型，如集中塑性 铰 模 型［1］、多 分 量 模 型［2-3］、分 段 变 刚 度 模型［4-5］、纤维模型［6］等.而关于等效节点荷载的研究则相对较少.在通常基于刚度法的有限元结构分析中，等效节点荷载只与非节点荷载的形式、分布及插值形函数等有关［7-8］.由于在分析过程中很难动态调整插值形函数，因此，等效节点荷载随着非节点荷载呈线性变化.这一近似处理虽能满足线性分析的精度要求，但在非线性分析中，由于梁单元端部会出现损伤、破坏，非线性变化的单元抗力与线性变化的等效节点荷载之间的矛盾变得突出起来，使节点变形与实际情况存在较大误差，也会增加平衡迭代计算次数，甚至引起计算不收敛的严重问题.另一方面，基于有限元柔度法建立的梁单元模型可以直接获得单元截面柔度，有条件根据单元损伤状态更新等效节点荷载.但是，由于单元柔度需要通过迭代确定，使非线性分析中更新等效节点荷载的计算量较大，同时与常规有限元分析流程不兼容.因此，虽然针对基于柔度法的梁柱单元模型已有大量的研究［9-12］，而关于等效节点荷载计算的改进则鲜有涉及.文献［6］提出了基于柔度法计算等效节点荷载，但只能在加载分析前进行，无法考虑单元损伤的影响.

笔者通过引入截面割线柔度，提出损伤梁单元等效节点荷载的计算方法，使等效节点荷载的变化与单元状态保持一致，一定程度上解决了前述计算量大以及有限元分析流程的兼容性问题.将这一方法引入自主开发的钢筋混凝土结构分析、设计软件ADARCS （advanced design and analysis of reinforced concrete structures）之中，进行了一系列结构分析研究，结果表明，上述改进确能提高结构非线性求解的效率和收敛性.

1 基于刚度法形成梁单元等效节点荷载的缺点

如图1所示，采用刚度法分析梁弯曲问题时，对于节点2的Hermite单元，非节点横向荷载（q（x），P1，P2，…）和 弯 矩（M1，M2，…）的 等 效 节 点 荷 载（Pe）可表示为［7］，

图1 承受横向荷载作用的等截面梁Fig.1 Uniform beam with transverse load

式中：N为形函数；ξ为单元自然坐标，ξt和ξr分别是单元内横向集中力和弯矩作用点的自然坐标；L为单元长度.

对于常见的均匀分布荷载，等效节点荷载为

从式（2）可以看出，在结构分析过程中，若不改变形函数，则单元的等效节点荷载不变.这样，在结构非线性分析过程中，等效节点荷载形式在非线性分析之前已经确定，并与节点荷载一起形成荷载向量，然后按照一定的增量方式加载.因此，在非线性分析中，等效节点荷载将与外加荷载呈线性规律变化，不会根据单元损伤状态进行调整.

事实上，当单元进入强非线性状态时，在梁单元端部将出现一定的刚度退化区域，甚至出现“塑性铰”，此时，按式（2）计算的等效节点荷载将大于实际的等效节点荷载.为了平衡结构内力，必然增加非线性平衡迭代次数，在结构承载力退化阶段，甚至会引起非线性迭代分析的不收敛.

2 基于柔度法计算损伤梁单元的等效节点荷载

众所周知，建立单元节点力与变形关系主要有刚度法和柔度法.在小变形情况下，梁单元任意截面的内力可以通过平衡关系准确地由杆端力表示，而且通过迭代可以平衡单元非协调变形.显然，基于柔度法确定的梁单元非线性状态能同时满足变形协调和平衡条件，而基于刚度法确定的单元状态则只能满足变形协调条件.因此，基于柔度法的非线性梁单元模型更适合精细化的结构分析.

在小变形条件下，根据柔度法（结构力学中的方法），图1 中q（x）的等效节点弯矩可由式（3）得到

对于以弯曲变形为主的梁单元，轴向变形和剪切变形的影响一般都很小，可以忽略不计.式（3）中的系数可以按如下公式计算：

式中：f（x）是单元的截面柔度；Mi（x）和Mj（x）分别是在i，j端的单位弯矩作用下单元的弯矩分布函数；Mq（x）是荷载q作用下单元的弯矩分布函数.式（4）中的积分可以采用高斯积分法求解，结合式（3）即可求得均布荷载作用下单元的等效节点弯矩（M1，M2）.同理，可以得到等效节点剪力（Q1，Q2）.对于承受轴向荷载的梁单元，也可以得到等效节点轴力（N1，N2）.

基于柔度法可以得到梁单元各种形式非节点荷载的等效节点荷载.在非线性分析的每一加载步，可以根据单元（编号为w，单元总数为ne）的损伤状态更新等效节点荷载（如图2所示）.但对处于非线性状态的单元，由于f（x）随受力状态发生变化，一般需要引入迭代分析确定等效节点荷载.而无论是采用增量迭代法或全量迭代法计算，所增加的计算量都是难以接受的；其次，在有限元非线性分析流程中，平衡迭代分析中需要更新单元的刚度或抗力，对应的单元状态也需要迭代分析确定，这样，单元刚度和荷载向量的形成处于不同的分析步骤中（如图2所示），迭代过程的差异将会出现同一个非线性状态对应不同单元截面柔度的情况.

图2 增量非线性分析流程图（第k 加载步）Fig.2 Flowchart of nonlinear analysis with incremental method（load step k）

上述问题使得柔度法形成等效节点荷载很难很好地融入常规结构有限元分析框架.为解决这一问题，提出了割线柔度法.如图3所示，对于某一确定的状态（如每一加载步的初始状态），引入截面割线柔度（fs（x））表示单元的截面变形和内力关系

式中，d和D分别为截面变形和内力.

图3 截面割线柔度变化示意图Fig.3 Variation of sectional flexibility

这样，更新等效节点荷载的计算过程将和线性分析一样（式（4）中的f（x）替换为fs（x）），不需要采用迭代分析.其次，在加载步长较小时，截面割线柔度fs在加载步内的变化不大，可以根据加载步起点的割线柔度计算等效节点力.显然，利用截面割线柔度可以反映梁单元的损伤状态，使单元等效节点荷载与结构损伤状态保持一致的变化趋势.同时，这一方法可以提高平衡迭代分析的效率和稳定性.

引入割线柔度计算等效节点荷载的流程如图4所示.结合图2可见，该方案合理地避免了不一致迭代分析中的矛盾，在不影响结构分析整体框架的情况下，实现了单元等效节点荷载的动态调整.

图4 割线柔度法计算等效节点荷载流程Fig.4 Calculation of equivalent nodal load based on sectional flexibility

3 应用实例

3.1 模型概况

框架模型取自某工业厂房，结构形式如图5 所示，边跨3.9m，中跨5.1m，纵向跨度4.5m，底层高6.0m，其余楼层高5.4m.1～3层为工厂车间，4和5层为办公室.C30混凝土，框架柱、梁截面尺寸分别为600mm×500mm 和700mm ×350mm，配筋如表1，2所示.

图5 框架结构模型Fig.5 Frame structure model

表1 柱单元截面配筋Tab.1 Reinforcement of columns

各楼层框架梁的设计线荷载及水平地震作用标准值如表3所示.其中，恒载包括楼板的自重及其传递过来的楼面恒载，活载来自于楼面.按照抗震设防要求，以7 度抗震设防，设计基本地震加速度为0.10g，Ⅱ类场地、第一设计地震分组，场地特征周期为0.35s，结构自振周期为0.40s.水平地震作用采用底部剪力法计算，考虑了楼板质量和楼板传递过来恒载和活载的等效质量.

表2 梁单元截面配筋Tab.2 Reinforcement of beams

表3 各楼层框架梁承受的均布荷载及水平地震作用标准值Tab.3 Uniformly distributed load on beamand earthquake action

为了模拟结构屈服、卸载阶段的性态，在地震作用拟静力全过程非线性分析中，水平推力采用地震作用标准值乘以一定的放大系数.

3.2 单元模型

采用分段纤维梁单元模型模拟梁、柱构件，如图6所示，AC、BD段为构件端部刚域部分，模拟框架节点对梁、柱单元的影响；CD段采用纤维模型，对CE和DF区域考虑材料（钢筋、混凝土）的非线性行为；EF段是弹性区域.如果单元弹性区域的长度（LE）为零，则该模型就是一般的纤维模型.根据相关研究［13-14］，取框架梁、柱单元非线性区域长度为2.5h0（h0为单元截面有效高度）.

图6 分段纤维梁单元模型Fig.6 Multi-zone fiber model

3.3 等效节点荷载与损伤的关系

模拟水平地震作用下的结构状态，应先进行重力荷载（恒载、活载等）作用下的静力非线性分析，然后以此为初始受力状态分析地震作用下的结构非线性.分析中，竖向荷载（恒载、活载等）不变，属典型的非比例加载分析.因此，竖向荷载等效节点荷载的计算考虑构件损伤很有必要，能提高分析的收敛性和稳定性.采用自主开发的钢筋混凝土结构分析、设计软件（ADARCS）对结构进行地震作用下的拟静力全过程非线性分析，以第一层梁B3和B4（见图5）为对象，针对恒载分别采用经典刚度法和本方法对比分析，说明随着结构非线性的发展，单元等效节点荷载的变化趋势，并分析等效节点荷载与损伤的关系.

从表4可见，在线性阶段，刚度法和本方法计算的等效节点荷载是一致的.然而，随着损伤的发展，结构进入较强的非线性发展阶段后，等效节点荷载必将显著变化.如前所述，按刚度法计算的等效节点恒载在整个非线性分析中是保持不变的，而如图7所示，按本方法计算的单元等效节点荷载则能根据构件损伤状态动态调整.对于B3梁，在非线性发展的初期，等效节点弯矩只是小幅变化；随着非线性的发展，损伤较为严重，局部截面混凝土开裂面积超过50%（如图8所示），此时严重损伤端部的等效节点弯矩大幅降低.同时，等效节点剪力也直观地反映了非节点力作用效应的动态变化.对于B4梁，虽然等效节点荷载变化规律不尽相同，但也反映了梁的损伤发展.从整个非线性发展过程看，由于损伤的存在，单元的等效节点弯矩小于线性阶段，且在损伤较为严重的梁端，退化更明显.为了说明梁端损伤发展状况，图8a，b给出了单元非线性区域（图6中的CE段和FD段）高斯积分点处的损伤情况（高斯积分点1，2，3位于单元I端，高斯积分点4，5，6位于单元J端（见图5），积分点3无损伤）.

表4 线性状态下梁的等效节点恒载Tab.4 Equivalent nodal dead load on beam in the linear state

对于不同的单元，由于截面特性、结构内力重分布等因素的影响，随着非线性的发展，等效节点荷载的变化趋势或许不同，但都是单元刚度分布变化导致荷载效应在单元内部重新分配的表现.

图7 B3和B4梁的等效节点恒载Fig.7 Equivalent nodal loads on B3and B4

图8 B3和B4梁积分点处混凝土开裂的发展Fig.8 Concrete cracking at integral points of B3and B4

3.4 计算效率和稳定性

本方法也是一种自适应的算法.它使结构的等效节点荷载，尤其是等效弯矩和结构的内力保持一致的变化，从而加速平衡迭代，提高分析效率.对3.1节的框架模型，图10是按照传统方法和本方法计算的基底剪力和顶层位移关系曲线.表5为分析步数和耗时对比，分析采用同一台计算机，荷载增量步等参数均取相同值.从表5 可以看出，采用传统刚度法，虽能节省每一加载步重新计算等效节点荷载的工作量，但整体计算效率反而更低.这主要是由于不变的等效节点荷载会增加平衡迭代分析的工作量，而且，结构进入强非线性状态后，在很多加载步中程序需要多次折减荷载增量，重新迭代分析以得到满足收敛标准的数值解，比较耗时.当结构进入屈服、卸载阶段，一部分梁、柱构件出现塑性铰.从图9 可以看出，传统方法由于不调整等效节点弯矩而使节点受力无法平衡，导致计算分析不收敛，无法模拟结构承载力下降严重的阶段.而本方法则能在合理控制分析误差的情况下模拟单元屈服引起的结构承载力完整卸载过程，得到结构全过程非线性发展路径，且不需要改变结构模型和加载形式.因此，本方法对结构屈服、卸载阶段的分析具有较好的自适应性.

图9 结构基底剪力与顶层位移关系曲线Fig.9 Caver of base shear force and roof displacement

表5 非线性分析步数及耗时Tab.5 Steps and time-consuming of nonlinear analysis

4 结论

（1）基于截面割线柔度法形成梁单元等效节点荷载，能反映单元的损伤情况，适合于结构非线性尤其是强非线性状态的分析.本方法不仅能提高非线性求解的稳定性、有效控制分析误差，而且不需要改变结构模型和加载形式.

（2）在强震等动力非线性分析中，对于动力作用过程中不变的荷载（如重力荷载），本文所建议的等效节点荷载计算方法完全适用.但是，若非节点荷载为动力荷载，其等效节点荷载的处理则需进一步研究.

［1］ Giberson MF. The response of nonlinear multi-story structures subjected to earthquake excitations［D］.Pasadena：Earthquake Engineering Research Laboratory，1967.

［2］ Clough R W，Benuska K L，Wilson E L.Inelastic earthquake response of tall buildings［C］∥Proceedings of the 3rd World Conference on Earthquake Engineering， Auckland and Wellington，1965：68-89.

［3］ Angelo D Ambrisi，Filip C Filippou.Modeling of cyclic shear behavior in RC members ［J］.Journal of Structural Engineering，ASCE，1999，125（10），1143.

［4］ Soleimani D，Popov E P，Bertero V V.Hysteretic behavior of reinforced concrete beam-column subassemblages［J］.ACI Journal，Proceedings，1979，76（11）：1179.

［5］ 汪梦甫，王海波，尹华伟，等.钢筋混凝土平面杆件非线性分析模型及其应用［J］.上海力学，1999，20（1）：82.WANG Mengfu，WANG Haibo，YIN Huawei，et al.Amodel for nonlinear analysis of reinforced concrete planar bar members and the application ［J］.Shanghai Journal of Mechanics，1999，20（1）：82.

［6］ 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，University of California，1991.

［7］ 王勖成.有限单元法［M］.北京：清华大学出版社，2003.WANG Xuncheng.Finite element method ［M］.Beijing：Tsinghua University Press，2003.

［8］ Crisfield MA.Non-linear finite element analysis of solids and structures（Volume 1）［M］.New York：John Wiley & Sons Inc.，1991.

［9］ Zeris C A，Mahin S A.Analysis of reinforced concrete beamcolumns under uniaxial excitation［J］.J Struct Eng，1988，114（4）：804.

［10］ Zeris C A，Mahin S A.Behavior of reinforced concrete structures subjected to biaxial excitation［J］.J Struct Eng，1991，117（9）：2657.

［11］ Spacone E，Ciampi V， Filippou F C.Mixed formulation of nonlinear beam finite element［J］.Comput Struct，1996，58：71.

［12］ Neuenhofer A，Filippou F C.Evaluation of nonlinear frame finite-element models［J］.Struct J Eng，1997，123（7）：958.

［13］ Barker A L L.Recent research in reinforced concrete and its application to design［J］.Jorunal of the Institution of Civil Engineers，1951，35（4）：262.

［14］ 朱伯龙，董振祥.钢筋混凝土非线性分析［M］.上海：同济大学出版社，1985.ZHU Bolong，DONG Zhenxiang.Nonlinear analysis of reinforced concrete［M］.Shanghai：Tongji University Press，1985.