空间刚架稳定性分析的泡函数有限元法

陈 梁，陈振中(沈阳航空航天大学 航空航天工程学部(院)，沈阳 110136)

空间刚架稳定性分析的泡函数有限元法

陈 梁，陈振中
(沈阳航空航天大学 航空航天工程学部(院)，沈阳 110136)

将泡函数有限元法理论用于空间刚架的稳定性分析。为了全面地考虑各种载荷因素对几何刚度矩阵的影响，提出了变几何刚度法和失稳判据。为了解决单元坐标系方位的确定这一个性化情况很突出的问题，运用三点定位算法编制了规范的计算程序。在一杆一单元的情况下，运用这三大理论编程计算出了比较精确的临界载荷值，说明在空间刚架的稳定性计算方面，泡函数有限元法具有突出的简洁性和精确性。

空间刚架结构；稳定性分析；泡函数；变几何刚度法；三点定位法；有限单元法

在生产实践中存在大量的空间刚架结构，如塔吊、汽车吊、输变电塔、钢构桥等，在遭受横向扰动的情况下这些刚架能否做到不失稳将是考验该类产品设计成败的关键。然而这类刚架存在杆件数量多，杆件截面空间布置的个性化突出等实际情况，这就使空间复杂刚架系统的稳定性分析具有很大的挑战性。王仁辉、高轩能等[1]将泡函数引入平面刚架的稳定性计算，并获得了简洁而精确的结果；龙驭球、袁驷等[2]提出极值点失稳问题的处理方案。在前人研究成果的基础上，本文将泡函数有限元法应用于空间刚架稳定性分析。为了全面地考虑各种载荷因素对几何刚度矩阵的影响，提出了变几何刚度法和失稳判据。为了解决繁琐的杆件单元坐标系的确定问题，利用三点定位算法编制了规范的计算程序。

1 泡函数有限元法

在单一载荷(即待求的临界载荷)作用下，结构的失稳形式为分岔点失稳。采用一杆一单元，以杆件的一端为原点建立单元局部坐标系o-x′y′，x′轴的正向指向杆的另一端，在刚架所在平面内y′轴垂直于x′指向杆件的一侧。 设杆件的长度为L，横截面积为A，弹性模量为E，抗弯轴惯性矩为I。根据文献[1]，只考虑弯曲变形影响的泡函数单元的形函数矩阵为

N=[N1，N2，N3，N4，N5],N1=(2η+1)(η-1)2,N2=Lη(η-1)2,N3=-η2(2η-3),N4=Lη2(η-1),N5=αη2(η-1)2+βη3(η-1)3

(1)

式(1)中N5为泡函数，α、β为任意常数，η=x/L，为单元的自然坐标。

当杆的内力为压力P时，由式(1)所示的形函数并利用虚功原理可得到平面刚架单元的变形刚度矩阵ke和几何刚度矩阵ge，如式(2)的前两式所示。

(2)

其中c1=EI/L3，c2=EA/L，a=4α2/5-12αβ/35+2β2/35，b=(α-3β/14)L，c=(44α2-22αβ+3β2)/77。令α=3，β=-3，则a=54/5，b=51L/14，c=621/77。

设平面刚架结构的整体坐标系为o-xy。x轴与x′轴间的夹角为φ。规定从x轴到x′ 轴以逆时针方向为φ的正向。则局部坐标系与整体坐标系间的坐标转换矩阵T如式(2)的第三式所示。

对于空间刚架结构，类似地采用一杆一单元，以杆件的一端为原点建立单元局部坐标系o-x′y′z′，x′轴的正向指向杆的另一端，而y′轴和z′轴均为杆件横截面的主惯性轴。杆件长度为L，横截面积为A，弹性模量为E，剪切模量为G，横截面对y′轴和z′轴的惯性矩为Iy和Iz，横截面的极惯性矩为J。设空间刚架问题的单元节点位移列向量为

q=(d1xd1yd1zφ1xφ1yφ1zd2xd2yd2zφ2xφ2yφ2z)T

(3)

根据文献[3]中所示平面刚架刚度矩阵和空间刚架刚度矩阵的继承与发展关系，并结合上述平面刚架泡函数有限元法的特点，得空间刚架单元的变形刚度矩阵ke和几何刚度矩阵ge，如式(4)所示。

(4)

其中P为单元轴向压力，a=EA/L，b=GJ/L，c=EIy/L3，d=EIz/L3，a1=4α2/5-12αβ/35+2β2/35，b1=(α-3β/14)L，c1=(44α2-22αβ+3β2)/77。令α=3，β=-3，则a1=54/5，b1=51L/14，c1=621/77。

设空间刚架结构的整体坐标系为o-xyz，nx′、ny′、nz′分别为各单元坐标轴的方向余弦，则单元坐标系与整体坐标系间的坐标转换矩阵

(5)

根据以上的稳定性理论并结合有限元法基本步骤，归纳出如下分岔点失稳问题的求解步骤。

(1)单位载荷作用下杆件内力的计算

在研究点(将要施加待求载荷F)处施加与F同方向的单位载荷，根据有限元法的基本流程算得各杆件的内力(P)，不必考虑结构自重的影响。

(2)单元刚度矩阵ke和ge的计算与组装

各单元都要计算ke，而只有当P<0时才要计算ge(式(2)中的第2式以-P代替P)，接着将它们分别组装成整体变形刚度矩阵K和整体几何刚度矩阵G。

(3)考虑约束条件和求解广义特征值问题

剔除受约束的自由度，即从K和G中删除该自由度所对应的行和列，从而得到空间刚架结构的变形方程为(K-FG)q=0，其具有非零解q的条件就是

|K-FG|=0

(6)

可利用逆迭代法来求解此广义特征值问题，得到最小特征值和相应的特征向量。最小特征值即为临界荷载Fcr，该特征向量即为屈曲模态，将此特征向量归一化后再结合节点坐标数据，可绘制屈曲模态实物图。

(4)评估泡函数对计算结果精度的影响。

2 变几何刚度法及失稳判据

根据文献[2]可知：和分岔点失稳问题处理法相比，在计算各单元的几何刚度矩阵时，极值点失稳问题处理法也只考虑了待求载荷(F)，而忽略了其它载荷(如结构自重、水平集中载荷、水平风载荷等)的影响，不同点在于式(6)的右端为其它载荷的等效节点力向量(Fp)，被称为结构的屈曲测试向量。此时的空间刚架结构的变形方程为

(K-FG)q=Fp

(7)

结构的屈曲条件为q中至少有一个元素是无穷大量。为了增强可操作性，提出了如下失稳判据：

|qi|max≥h/100,i=1…n

(8)

其中h表示刚架高度，n表示总体线位移自由度数。满足此判据的最小F值即为Fcr。

因为都要求q的系数矩阵行列式为零，所以上述两种方法求得的临界载荷值也相差不大。由于具有相同且固定的G矩阵，可将这两种方法归为一类，即定几何刚度法。

结构承受多种载荷作用时，宜采用极值点失稳问题处理法。当待求载荷和其它载荷对结构的变形具有相当的影响时，仍然采用处理极值点失稳问题的定几何刚度法将会得出错误的计算结果，给结构设计带来巨大的风险。为了克服这个缺点，处理极值点失稳问题的变几何刚度法被提出。该法的变形方程为

(K-G)q=Fp

(9)

其中G矩阵包含了全部载荷的影响。变几何刚度法和迭代法紧密结合，其具体实施步骤如下：

(1)采用分岔点失稳问题处理法计算出F的参考值F0，同时也生成K矩阵和Fp向量。

(2)采用变几何刚度法计算出F的临界值Fcr

取F0/2作为F的迭代初始值，上限值为3F0/2，步长为单位力值。考虑所有载荷的影响，计算出各杆件的内力(P)，按式(4)计算单元几何刚度矩阵，并组装成总体几何刚度矩阵G。式(9)中剔除受约束的自由度后，可解得总体节点位移向量q。这个迭代过程将循环地进行下去，直至式(8)所示的判据得到满足，此时所得的F值即为Fcr。

3 三点定位法

图1为由角钢搭建而成的空间刚架结构[2]，确定各杆件单元坐标系o-x′y′z′的方位是一个繁琐而重要的问题。图2所示结构为取自此类空间刚架的一部分。

图1 铁路控制塔计算模型图

图2 单元坐标系示意图

设在整体坐标系中A点坐标为(x1,y1,z1)，C点坐标为(x2,y2,z2)，B点坐标为(x3,y3,z3)。根据下面的一系列公式，可由A、C和B三点的整体坐标确定斜杆AC的单元坐标系。类似地，可用三点定位法确定斜杆DB和其它各斜杆的单元坐标系。

(10)

(11)

4 算例

图1中x-y平面为水平面，z方向为铅垂方向。它具有28个连接点和96根杆件，其中24根立柱采用L50×5型角钢(1#)，24根横杆和48根斜杆均采用L56×5型角钢(2#)。两种角钢截面的参数如表1所示。采用Q345型钢材，其弹性模量为210 GPa，泊松比为0.3，密度为7 800 kg/m3。根据一杆一单元的网格划分法可知，此刚架结构具有28个节点和96个单元，25-28号节点固定。

(1)在1节点的-z轴方向上施加集中载荷F且不计结构自重时，采用处理分岔点失稳问题的定几何刚度法求F的临界值Fcr1。

(2)在1节点的-z轴方向上施加集中载荷F且不计结构自重时，采用处理极值点失稳问题的定几何刚度法求F的临界值Fcr2。

(3)在1节点的-z轴方向上施加集中载荷F和刚架背面(外法线方向为-x轴方向)承受水平风压载荷且考虑结构自重时，采用处理极值点失稳问题的变几何刚度法求F的临界值Fcr3(单位长度风载荷p=-2ρbv2，风速v=40 m/s，空气密度ρ=1.29 kg/m3，b为角钢断面宽度)。

角钢截面参数如表1所示，其中S表示角钢截面号，q表示单位长度的质量，b为角钢断面宽度，A表示截面面积，J、Iy和Iz分别表示截面的扭转惯性矩和两个主轴惯性矩。节点坐标数据如表2所示，其中N表示节点号，x,y和z表示节点的3个坐标值。

表1 角钢截面参数

表2 刚架的节点坐标

对于24根立柱，有8种单元坐标系；对于24根横杆，有4种单元坐标系。根据几何关系可以直接写出这12种单元坐标系的方向余弦矩阵Λ。而对于48根斜杆，将三点定位算法编成规范的计算程序，可快速地获得全部斜杆的Λ矩阵。

(1)分岔点失稳(定几何刚度法)

将上述泡函数有限元法和三点定位法理论编制成MATLAB程序，按前述分岔点失稳问题的求解步骤，最终求得此空间刚架结构的临界载荷Fcr1=951.87 kN。

当采用普通有限元法即不考虑刚度矩阵(式(4))中的泡函数项时，适当修改前述的MATLAB程序后求得此结构的临界载荷值为2 282.1 kN(一杆一单元)和952.57 kN(一杆二单元)。两类算法计算结果的详细比较如表3所示。由此可知：一杆一单元的泡函数有限元法可达到一杆二单元的普通有限元法的精度且计算时间更短；一杆一单元的普通有限元法的计算精度很差。

表3 两类定几何刚度分岔点失稳有限元算法的对比

(2)极值点失稳(定几何刚度法)

根据式(7)可知K、G和Fp均不变，只需计算一次。采用迭代法求得临界载荷Fcr2=951.88 kN。由此可知处理分岔点失稳问题的定几何刚度法和处理极值点失稳问题的定几何刚度法所求得的临界载荷值相差不大，两种方法没有本质区别。

(3)极值点失稳(变几何刚度法)

根据式(9)可知K和Fp均不变，只需计算一次，但G矩阵随F的变化而变化。采用迭代法求得临界载荷Fcr3=939.38 kN，明显小于Fcr1和Fcr2。

由此可知，采用变几何刚度法来处理极值点失稳问题具有安全性更高的优点。

5 结论

空间刚架泡函数有限元法的突出优点是既简洁又准确。通过实际算例，如下的三大创新点得到了验证。

(1)采用泡函数法，在保证计算精度的前提下大幅降低了问题的自由度数；

(2)提出变几何刚度法和失稳判据，全面考虑了各种载荷因素对几何刚度矩阵的影响，优化了极值点失稳理论；

(3)采用三点定位法，简化了单元坐标系的定位问题。

而目前的各种商用有限元分析软件都不具有上述功能，因此本文所论述的空间刚架稳定性分析的泡函数有限元法具有极大的实用价值。

[1] 高轩能,王仁辉.框架结构稳定分析的泡函数有限元法[J].华侨大学学报(自然科学版),2006,27(2):155-158.

[2] 龙驭球,包世华,匡文起,等.结构力学Ⅱ(第2版)[M].北京:高等教育出版社,2006.

[3] DARYL L.LOGAN.有限元方法基础教程(第五版)[M].张荣华,王蓝婧,李继荣,等译,北京:电子工业出版社,2014.

[4] 刘文国,任文敏,张维.用复合有限条-气泡函数法研究加肋壳体稳定性[J].清华大学学报(自然科学版),1999,39(11):61-64.

[5] 夏拥军,陆念力.梁杆结构稳定性分析的高精度Euler-Bernoulli梁单元[J].沈阳建筑大学学报(自然科学版),2006,22(3):362-366.

[6] 刘慧娟,罗永峰,杨绿峰,等.单层网壳结构稳定性分析的随机缺陷模态叠加法[J].同济大学学报(自然科学版),2012,40(9):1294-1299.

[7] 周 飞,宛树旗.单层网壳结构稳定性分析[J].国外建材科技,2007,28(3):119-122.

[8] 彭细荣,杨庆生,孙卓.有限单元法及其应用[M].北京:清华大学出版社,2012.

[9] 徐斌,高跃飞,余龙.MATLAB有限元结构动力学分析与工程应用[M].北京:清华大学出版社,2009.

[10]黄志新,刘成柱.ANSYS Workbench 14.0超级学习手册[M].北京:人民邮电出版社,2014.

[11]龙述尧,蒯形成,刘腾喜.计算力学[M].长沙:湖南大学出版社,2007.

[12]P I KATTAN.MATLAB有限元分析与应用[M].韩来彬译,北京:清华大学出版社,2006.

[13]CHAJES A.Principles of Structural Stability Theory[M].New York:Prentice-Hall Inc,1974.

[14]WEAVER W,JOHNSTON P R.Finite element for ftructural analysis[M].New York:Prentice-Hall Inc,1984.

[15]CAO JIWEI,LIU MINGFANG,CHEN SHAOCHUN.Mixed Finite Element Formats of any Order Based on Bubble Functions for Stationary Stokes Problem[J].Journal of Math,2016,31(1):87-95.

Abubble-functionfiniteelementmethodforthestabilityanalysisofspaceframestructures

CHEN Liang，CHEN Zhen-zhong
(Faculty of Aerospace Engineering,Shenyang Aerospace University,Shenyang 110136,China)

Stability of space frame structures was analyzed by a Bubble-Function Finite Element Method(BFFEM)in this paper.Both a variable geometric stiffness method and a buckling criterion were presented with the consideration of various load influences on the geometric stiffness matrix.A formal program was written by a Three Points Method to define different element coordinate systems.For one bar one element,a very precise critical load for the space frame was calculated using the above three theories,which indicated that the BFFEM was simple and accurate for analyzing the stability of space frame structures.

space frame structures;stability analysis;bubble-function;variable geometric stiffness;three points method;finite element method

2017-10-08

陈 梁(1988-)，男，湖南长沙人，硕士研究生，主要研究方向：结构动力学分析与有限元法应用，E-mail:chenliang0102@163.com；陈振中(1963-)，男，辽宁沈阳人，教授，博士，主要研究方向：结构疲劳断裂与可靠性分析，E-mail:zhenzhong_chen@hotmail.com。

2095-1248(2017)06-0040-06

O302;TU328

A

10.3969/j.issn.2095-1248.2017.06.007

吴萍 英文审校：赵欢)