基于共旋法与稳定函数的几何非线性平面梁单元

邓继华，谭建平，谭 平，田仲初

(1. 长沙理工大学土木工程学院，长沙 410076；2. 广州大学工程抗震研究中心，广州 510405)

随着计算机软硬件技术的发展，基于非线性有限元开展结构破坏乃至失效倒塌分析已经成为可能[1 − 2]，这其中研究开发出各种高效精确的非线性单元模型和算法是基础和前提[3 − 5]。对于几何非线性梁单元，许多学者从考虑应变的高阶项推导精确的单元切线刚度矩阵、杆端力与变形后单元坐标转换矩阵的精确计算等入手，进行了卓有成效的研究[4, 6]。在各种几何非线性有限元列式中，共旋坐标(CR)法认为结构的几何非线性主要由结构大位移引起，因此在分析中每一个单元都建立一个共旋坐标系(局部坐标系)，该坐标系随单元位移发生平移和转动，该坐标系下单元的纯变形可以通过扣除单元刚体平移和转动而得到，相应可求得单元切线刚度矩阵和抗力，再通过变换的坐标转换矩阵将其转换到结构坐标系下，此过程中几何非线性效应被自动计入[7 − 9]。研究表明，基于共旋坐标法研究几何非线性，不仅具有列式简单、计算准确，以及适用于杆、梁、平面、壳、实体单元等各类型单元的优点[10 − 14]，而且能将几何与材料非线性脱藕，即材料非线性在局部坐标系内考虑，几何非线性则通过局部坐标系与结构坐标系之间的转换予以实现[15−16]。对于梁单元而言，除了大位移能引起几何非线性以外，梁-柱效应也是引起几何非线性的一个重要因素[17]，这其中采用能考虑构件轴力与弯矩相互作用的稳定函数[18 − 19]来考虑梁柱效应是一种最精确的方法[20 − 21]，实际工程分析中常采用的几何刚度矩阵实质是将稳定函数法表达的单元刚度矩阵用级数展开，取其一阶近似而得到[22]。但应指出的是，以往基于稳定函数法的研究采用的有限元列式一般为T.L或U.L 列式，笔者尚未见到将稳定函数法与共旋坐标(CR)法结合起来进行几何非线性平面梁单元的研究。因此，本文基于共旋坐标法，在局部坐标系下分离单元刚体位移和变形，采用稳定函数考虑梁柱效应；再基于结构坐标系与局部坐标系下节点力之间及节点位移之间的总量关系及微分导出的增量关系，获得平面梁单元在结构坐标系下的全量平衡方程和切线刚度矩阵；在此基础上根据带铰梁端弯矩为零的受力特征，导出了能考虑梁端带铰的单元切线刚度矩阵。最后通过数个经典算例对本文算法及程序进行了较全面和严格的检验。

1 平面梁单元切线刚度矩阵

如图1(a)、图1(b)所示分别为初始时刻和计算t时刻的平面梁单元，XY为结构坐标系，xy为单元的共旋坐标系，它具有随单元变形而转动的特点，在运动中始终以节点i为原点，以节点i到j的连线方向为x轴，y轴则由x轴按逆时钟方向旋转90°形成。

图1 变形前后平面梁单元Fig. 1 Plane beam element before and after deformation

初始时刻设单元在结构坐标系下水平及竖向的投影长度为x0和y0，在计算t时刻结构坐标系中的位移向量d=[ui viθi uj vjθj]T，在共旋坐标系中的位移向量dL=[uLi vLiθiLuLj vLjθLj]T，联立图1(a)、图1(b)，有：

平面梁单元在共旋坐标系中的节点位移为：

其中：

式中：si=sinαt,co=cosαt。

图2(a)、图2(b)所示分别为梁单元及微段的受力平衡情况。

图2 梁单元及微段Fig. 2 beam element and micro-segment

对式(11)进行与单元受拉时同样的推导，可得到与式(9)相同的表达式，只是各参数s、c、k、ω的计算式变成以下形式：

kT的具体表达式较繁琐，限于篇幅此处未列出。

联立式(19)和式(20)，可将式(17)写成：

2 梁端带铰时的单元切线刚度矩阵

3 非线性方程组求解和计算流程

非线性方程组的增量法求解主要有荷载增量法和位移增量法[23]。理论上而言，荷载增量法只能计算出荷载-位移曲线的上升段，而位移增量法能计算荷载-位移曲线的下降段，从而能获得极限荷载值与位移值；但在实际计算中发现，由于在荷载-位移曲线的顶点附近结构切线刚度接近奇异，即使采用位移增量法，非线性计算也很难收敛。鉴于此，很多文献都大力推荐弧长增量法，文献[24]对传统弧长增量法[25]研究后认为，该方法能顺利进行的前提是必须对结构的特性有深刻的了解并进行复杂的试算，显然这会给分析工作带来极大困难和不便；因此，该文献在传统弧长增量法的基础上提出一种如图3 所示的基于牛顿-拉菲逊法的投影增量法，该方法能有效克服传统弧长增量法的上述缺点，且该方法收敛速度要稍快，计算量也略小，本文采用该方法进行非线性方程组的求解。

本文非线性计算流程如下：

1)根据式(2)由上一荷载步末单元i、j节点在结构坐标系下的总位移向量d求出共旋坐标系下的节点位移向量dL；

2)由式(13)中最后一行单元轴力N的计算式得到N，再根据轴力N是受拉还是受压分别基于式(10)或式(12)计算参数s、c、k、ω，进而根据式(9)得到单元弯矩Mi和Mj；

3)根据式(13)可得到单元在共旋坐标系下的杆端力向量为f；

4)由式(4)和式(16)分别计算矩阵T和t；

5)由式(19)计算Kn，再由式(22)得到单元在结构坐标系下的切线刚度矩阵K；

6)由式(16)得到单元在结构坐标系下的节点抗力向量F；

7)对所有单元重复步骤1)～步骤6)，基于常规的“对号入座”原则叠加形成结构总刚矩阵和总抗力矩阵；

8)由总抗力矩阵与总荷载矩阵之差形成不平衡力矩阵作为非线性平衡方程的右端项，解方程得到增量位移矩阵 ∆d，与前述总位移向量d叠加形成新的总位移向量d；

9)判断是否收敛，如是，转至下一个荷载步；如否，进行下一轮迭代计算。

图3 投影增量法示意Fig. 3 Projection increment method

4 算例分析

图4 Lee’s 框架及荷载 /cm Fig. 4 Lee’s frame geometry and loading

图5 Lee’s 框架Fig. 5 Lee’s frame

采用3 个经典算例对本文算法及程序从正确性、计算效率及精度方面进行较全面和严格的检验。算例1. 如图4 所示为Lee’s 框架[26]，材料及几何参数如下：线弹性模量E=7060.8 kN/m2，梁和柱具有同样的截面，且截面面积为6 cm2，惯性矩为2 cm4，对结构进行几何非线性分析，分析中柱均分为8 个单元，荷载作用点左边梁均分为2 个单元，右边梁均分为6 个单元，图5(a)所示为框架的初始构型以及不同变形时刻的构型，图5(b)、图5(c)所示分别为荷载作用点处荷载-水平位移及荷载-竖向位移曲线，可看出与相关文献[27]的计算结果很吻合(该文献未提供荷载-水平位移曲线)。算例2 . 如图6 所示悬臂梁，集中荷载P=35 kN作用在自由端，梁跨为L=100 m、弯曲刚度为EI=3.5×104kN·m2、梁端水平位移为u、垂直位移为w、转角为θ，对本算例按以下两种方法求解：

图6 端部承受集中荷载的悬臂梁Fig. 6 Cantilever subjected to concentrated load at free end

方法1. 采用本文算法研制的程序进行计算；

方法2. 按文献[11]和文献[28](文献[28]就是将文献[11]的算法从平面梁延伸到空间梁)的方法进行计算，即共旋坐标系下的杆端力就由节点位移(去除了刚体位移)与小位移线弹性刚度矩阵相乘得到。

为使计算结果具有可比性，两个程序采用的方程组求解方法、数据精度及收敛精度标准等都完全相同，方法2 中杆端力也是基于全量计算，详见文献[11]及文献[28]。

两种方法得到的计算结果及比较见表1 和表2，解析解见文献[29]。可以看出，无论是从非线性计算能力、精度还是稳定性方面，方法1 均优于方法2(方法2 在划分为1 个单元及4 个荷载步情况下还出现错误的收敛解)；还可以看出，由于方法1 与方法2 中杆端力的计算均是基于全量计算，所以在单元数划分固定的前提下，计算精度与荷载步数无关。

表1 悬臂梁自由端的位移Table 1 Displacement of cantilever at free end

表2 计算误差 /(%)Table 2 Calculation error

算例3. 如图7(a)所示，一对拉力2P作用于铰接方棱形框架的对角点，各项几何、材料参数及位移变量详见图7，EI为杆件的抗弯刚度。

图7 对角点受拉力作用时的铰接方棱形框架Fig. 7 Diamond-shaped frame under a pare of opposite concentrated tensions in opposite angles

为验证本文所导出的带铰梁单元模型，基于荷载及结构对称的特征，取图7(b)所示的1/2 结构进行计算，按每根杆件分成1、2、3 和4 个单元(中间铰左、右两个相邻单元类型选用本文推导的带铰单元)，在荷载参数PL2/EI=10 时所得的计算结果如表3 所示，其中解析解见文献[29]。从表3可看出，计算结果收敛非常好，如按工程精度误差不超过5%的要求，则每根杆件划分成3 个单元就已经足够。

表3 铰接方棱形框架位移Table 3 Displacement of diamond-shaped frame

5 结论

本文在局部坐标系下采用稳定函数考虑P-δ效应，再基于结构坐标系与局部坐标系下节点力之间及节点位移之间的总量关系及微分导出的增量关系，获得平面梁单元在结构坐标系下的切线刚度矩阵，同时还获得单元杆端力的全量算法；在此基础上，根据带铰梁端弯矩为零的受力特征，导出了能考虑梁端带铰的单元切线刚度矩阵表达式。多个算例表明本文算法是正确的，相对于已有文献算法，本文算法在非线性计算能力、精度以及稳定性方面均占有一定优势，可用于平面杆系结构的几何非线性分析。