陈 楠 濮嘉铭 刘 龙
上海海事大学 上海 201306
摘 要：结构材料非线性行为分析是结构稳定性研究的重点。基于向量式有限元三维梁单元的基本理论，引入材料弹塑性模型和单元进入塑性的屈曲准则，从而实现对空间框架结构的弹塑性行为分析。利用Matlab 编制相应的数值计算程序，以两个典型框架为例，并与文献的计算结果进行对比，验证了理论推导和程序的准确性。结果表明，向量式有限元理论能够准确地应用于分析框架结构的弹塑性行为，较为真实地模拟出框架结构受载荷时会产生塑性铰的现象。
关键词：结构稳定性；向量式有限元；塑性铰
中图分类号：TU313 文献标识码：A 文章编号：1001-0785（2020）20-0099-05
0 引言
随着起重载荷、起升高度等参数的不断提高，起重机械钢结构朝着高耸、柔长、轻量化方向发展，结构稳定性问题更加突出，造成严重后果。因此对起重机结构进行非线性稳定性研究，分析臂架失稳失效具体过程和获取准确的稳定极限承载力对臂架的结构设计十分重要，理论分析、有限元法等为分析结构的弹塑性问题提供了有力工具，刘永华[1] 通过弹性屈曲分析和考虑大变形的弹塑性非线性分析考察结构的非线性行为，利用Ansys 软件进行了仿真验证；高健利[2] 推导了高阶非线性空间桁架单元，基于OpenSEES 软件开发了考虑剪切变形影响的几何非线性梁柱单元。孟丽霞[3] 推导出计及二阶效应的惯性矩二次变化变截面Timoshenko梁单元转角位移方程，提出了一种可用于分析变截面Timoshenko 梁单元几何非线性大位移问题的计算方法。王金平[4] 提出了基于共旋坐标法的几何非线性子结构方法，通过对微分方程的求解并结合工程中的结构失稳判断准则得到臂架结构的失稳载荷。赵二飞[5] 从桁架臂稳定性计算过程中的长度系数和稳定系数出发，利用梁柱理论建立臂架力学模型和稳定方程，得到了交叉式和点对点式桁架臂弦杆的屈曲规律，并提出了桁架臂弦杆等效长度系数和拟合稳定曲线作为桁架臂局部稳定性的计算依据。胡燕东等[6] 从精确有限元理论出发，推导了桁架结构的临界稳定表达式，改进了有限元软件分析超静定结构稳定性的方法。
传统弹塑性分析中需要先作预应力分析，不断更新结构形态，计算效率较低且易遇到求解不收敛的问题。向量式有限元(vector form intrinsic finite element，简称：VFIFE) 是美国普渡大学丁承先[1] 教授结合向量力学和有限元概念提出的一种新型力学计算方法，能够有效求解结构的弹塑性、屈曲失稳、碰撞、倒塌等复杂力学行为。利用此优势，在非线性梁理论研究中，黄正[8] 等运用向量式有限元梁单元理论，推导了计算构件内力的非线性公式，并考虑弯扭耦合变形；曲激婷[9] 等基于向量式有限元建立半刚接钢框架结构模型，对结构的竖向倒塌进行了动力学响应分析；陈冲[10] 利用向量式有限元梁理论，分析受到外载冲击后，结构从局部失效到整体倒塌的全过程力学特性。
本文基于向量式有限元三维梁单元理论，编制Matlab 计算程序，进行空间框架的弹塑性研究。通过计算单层单跨空间框架和六层空间框架两个案例，对比弹性分析和弹塑性分析，展示了结构材料非线性变形过程。
1 向量式有限元基本原理
向量式有限元基于点值描述、途径单元和逆向运动三个基本概念。首先将结构的空间形态由有限数量的点值描述，构成节点和单元；通过牛顿运动定律控制节点运动，利用单元逆向运动获取纯变形以计算内力，再通过坐标转换来确定方向；将结构在一时间段内的运动变形划分为多个时间微段的途径单元，每个途径单元都会更新局部坐标系，实现结构位置形状的步步更新。向量式有限元的计算流程如下：
1) 确定结构的材料、截面性质、阻尼、计算步长及总时长等参数。
2) 将结构划分为节点和梁杆单元。
3) 通过中央差分控制方程计算节点的运动。
4) 由逆向运动获得单元的纯变形，并确定内力的大小和方向。
5) 更新控制方程中的外力和内力。
6) 判断是否达到计算终止时间。若否，则重复3)、
4) 和5)；若是，则输出分析结果。
1.1 控制方程及求解
向量式有限元研究对象是描述结构的一系列空间质点，故其控制方程是一组由牛顿运动定律得到的点运动公式，并运用显式的中央差分公式进行求解。对单个质点的平动和转动微分方程：
式中：m、I 分别为质点的质量、质量惯性矩；x、θ 分别为质点的线位移、角位移；F、M 分别为作用在质点上的力、力矩。
将上述运动控制方程转化为中央差分形式：
式中：h 为时间增量的步长，可设其为常数；n 为步数。根据胡克定律，受力后的无阻尼结构将维持持续振动的状态，因此需在方程中加入虚拟的计算阻尼力，上述差分公式则可变为
1.2 梁单元内力计算
向量式有限元中梁单元分为平面梁单元和空间梁单元。平面梁单元的节点位移有3 个分量，分别是两个平移量和一个转动量，单元内力除了轴向力，还有一组剪力和弯矩。对于空间梁单元，其节点分量包括3 个平移量和3 个转动量，经过简化后，单元内力包括轴力、扭矩和两组剪力和弯矩。另外，平面梁单元的局部坐标系中其中一个主轴方向（垂直于运动平面）不作变化，而空间梁单元中局部坐标系三个主轴方向均会发生改变。
单元内力大小与纯变形有关，向量式有限元利用逆向运动概念来计算纯变形，建立局部坐标系并根据虚功原理得到内力大小和方向，再通过正向运动（坐标转换）让内力方向回到原来位置。如图1a 所示，杆件单元1a—2a 经运动变形到达1-2，令其作一个虚拟的逆向运动，包括平移（-u1）和转动（-θ），则可得到纯变形为两者的长度变化Δ。
(a)
(b)
图1 逆向运动及变形计算
对于梁单元，其纯变形除了线位移还包括角位移。如图1b 所示，为梁单元由ta 运动至t 时的形态，u1 为节点1 的位移。单元1-2 对1a-2a 的转角θ 为刚体转角，节点转角θ1、θ2 分别为β1a、β2a 与β1、β2 的差量。节点转角θ1、θ2 中，一部分为单元的刚体转角，另一部分为单元自身几何变形产生的弯角。可得到梁单元两节点的弯角即纯变形φ 1、φ 2 为
式中：E 为弹性模量，l 为单元长度，Δ 为单元变形，A 为单元截面积。
相同的处理方法，可以得到梁单元的转角和另一组弯角，对应计算得到单元的扭矩和另一组剪力和弯矩。
1.3 弹塑性模型
结构的破坏中，一般同时伴随着几何和材料非线性问题，只有同时考虑几何大变形和材料弹塑性才能比较真实的反应结构屈曲破坏过程。
为了分析方便，材料采用如图2 所示的双线性强化弹塑性模型，在应力达到屈服应力σy 以前为完全弹性，应力—应变曲线的斜率为弹性模量E，达到屈服应力进入塑性，斜率为切线模量Et，一般可取Et=0.03E。
图2 双线性强化弹塑性模型
1.4 屈服准则
材料随着载荷增加，由弹性阶段进入塑性阶段，结构的力学行为也因此发生显著变化。运用向量式有限元进行弹塑性分析时，需要施加一个屈服准则，用以判断单元是否进入塑性变形阶段。屈服准则与单元的应力或者应变相关联，一般用应力作为界定条件在计算中更为方便。该屈服准则对应的单元应力条件，就是单元弹性变形阶段应力的上界，也是塑性变形阶段应力的起始。
主要考虑轴力N 和弯矩M 的单元的屈服条件
式中：MY 和NY 分别为单元内开始出现屈服时的弯矩和轴力。
2 数值算例
2.1 单层单跨空间框架
单层单跨空间框架的几何尺寸和载荷施加情况如图3 所示[2]。构件弹性模量和剪切模量分别为E=2×105MPa 和G=7.929 3×104 MPa， 屈服强度σy=234.43MPa。框架的梁和柱的截面分别为W18×60 和W10×60，长度均为L=3.657 6 m，框架的底端固定，柱顶作用有四个不同大小的竖向载荷，其中一个柱顶作用有一水平载荷，外载荷均为为等比例施加。
图3 单层单跨空间框架
运用向量式有限元对空间框架进行分析时，需要选择空间梁单元建模计算。将框架的每根构件划分为4 个单元，分析得到载荷F 和侧向位移Δ 的关系曲线如图5、图6 所示，并将向量式有限元分析结果与参考文献2 进行对比。
如图4a 所示，不考虑材料非线性而进行弹性分析时，结构的整体侧移和外载荷呈一阶关系，结构的弹性稳定极限载荷远远大于参考文献的极限载荷值。采用只考虑轴力的双线性强化弹塑性材料模型时，结构进入塑性阶段，载荷—位移曲线斜率逐渐减小然后趋于定值，其极限载荷也远大于参考文献。
在向量式有限元中引入综合考虑轴力和弯矩的塑性承载力函数（屈服准则），对构件的强轴和弱轴方向分别给定塑性弯矩承载力进行计算，并用CRC 切线模量考虑残余应力的影响。分析结果如图4b 所示，由曲线对比可知，向量式有限元的计算结果与参考文献很接近，两者的极限载荷值误差在5% 以内。
(a)
(b)
图4 单层单跨空间框架弹性及弹塑性变形曲线
2.2 六层空间框架
六层空间框架如图5 所示[2]，材料的弹性模量为E=2.068 5×105 MPa，屈服强度为σy=250 MPa，泊松比为ν=0.3。框架每层的均布重力载荷大小为9.6 kN/m2，可等效为作用在每层柱顶的集中载荷。沿Y 方向作用有风载荷，可等效为作用在框架前立面的所有梁柱的连接点的集中载荷，大小为53.376 kN。
采用向量式有限元空间梁单元对此复杂空间框架进行精细塑性铰高等分析计算。综合考虑计算精度和效率，每根梁柱均划分为3 个单元。分析得到框架顶部A 点X、Y 轴方向的位移与外载荷的关系曲线，并与参考文献2的载荷—位移曲线对比，如图6 所示。
图5 六层空间框架
(a)
(b)
图6 六层空间框架弹性及弹塑性变形曲线
不考虑材料非线性的影响，首先对框架进行了弹性分析，如图6a 所示，为保证一定计算效率，受到向量式有限元时间步长和划分单元数目的影响，位移响应存在误差，但其仍在合理范围内。通过塑性铰模型考虑材料非线性，对此框架进行高等分析计算，结果如图6b所示，与参考文献对比可以看出，两者进入塑性的载荷情况基本一致，两曲线的大部分误差在10%～20%之间，向量式有限元计算只需约2 min，而一般有限元计算需要约10 min。分析结果表明：基于向量式有限元的精细塑性铰高等分析具有其较好的有效性和计算效率。
3 结论
本文基于向量式有限元梁单元理论，应用Matlab软件编制程序进行框架结构的弹塑性行为分析，不需要作特殊处理即可有效分析大变形、大变位、弹塑性等复杂力学问题。主要得出以下结论：
1) 考虑几何和材料非线性，向量式有限元能够有效求解得到空间框架结构的弹塑性变形，从而为工程上结构极限分析提供参考。
2) 对于载荷- 位移曲线单纯的弹性分析不能够准确描述结构的变形，而弹塑性模型和屈曲准则的引入，能够很好地解决这一问题。
3) 本文计算结果表明：向量式有限元直接通过对所有节点和单元的数值计算，能够判断结构弹塑性变形，并能有效获取结构的屈曲行为。
本文工作还有待深入，可以进一步提高向量式有限元法计算精度和效率；结合可靠的试验数据，对计算结果进行修正；引入向量式有限元板单元、固体单元等，进行精细化结构非线性分析。
参考文献
[1] 刘永华. 空间钢框架高等分析方法研究[D]. 哈尔滨：哈尔滨工业大学,2007.
[2] 高健利. 基于OpenSees 的钢结构高等分析二次开发[D].郑州：郑州大学,2016.
[3] 孟丽霞. 起重机变截面复杂梁杆系统稳定性与非线性大位移研究[D]. 哈尔滨工业大学,2013.
[4] 王金平. 起重机臂架结构几何非线性稳定性分析[D].大连: 大连理工大学,2015.
[5] 赵二飞. 起重机桁架臂局部稳定性研究[ D ] . 长春:吉林大学,2018.
[6] 胡燕东, 兰朋, 陆念力. 高耸超静定桁架结构的稳定性分析[J]. 建筑机械化, 2010, 31(9):50-55.
[7] 丁承先, 段元锋, 吴东岳. 向量式结构力学[M]. 北京：科学出版社,2012.
[8] 黄正,刘石,杨毅,等.基于非线性梁理论的有限质点法[J].计算力学学报,2019,36(5):610-617.
[9] 曲激婷, 李志彩. 半刚接钢框架结构抗竖向连续倒塌动力响应分析[J]. 大连理工大学学报,2019,59(4):409-416.
[10] 陈冲. 基于向量式有限元的索杆结构精细化分析和倒塌破坏研究[D]. 杭州: 浙江大学, 2017.
