基于节点坐标变量的桁架大位移非线性有限元法

文章编号：16732049(2014)02013805

收稿日期：20131117

摘要：为求解桁架大位移问题, 提出了一种基于节点坐标变量的非线性有限元法——以杆端节点坐标向量为显函变量写出单元杆端力向量表达式,由单元杆端力向量装配结构非线性平衡方程。求解时,首先根据矩阵微分理论求出单元杆端力向量关于杆端节点坐标向量的导数矩阵,由该导数矩阵装配结构非线性微分平衡方程;然后按照Newton切线法原理建立等效线性逼近方程,引入边界约束条件得到结构节点坐标的迭代公式。研究结果表明：该方法稳定性好、精度高、收敛速度快且简单易用, 为求解桁架大位移问题提供了一种有效方法。

关键词：非线性分析;Newton切线法;大位移;有限元法;平面扁桁架；节点

中图分类号：TU311.4 文献标志码：A

Nonlinear Finite Element Method Solving Large Displacement Problems of Trusses Subjected to Node Coordinate Variable

LIU Shutang

Abstract: In order to solve the large displacement problems of trusses, a nonlinear finite element method was proposed, which used the node coordinate variables as unknowns. An expression of the member end force vector was written in terms of members end coordinates, by which the global nonlinear equilibrium equation was fabricated. In solution, firstly, based on matrix differential theory, element derivative matrix was obtained with respect to the member end force vector, by which the global differential nonlinear equilibrium equation was fabricated. Secondly, based on Newton tangent method theory, the global equivalent linear matrix equation was established, then the structures boundary restraints was introduced, the iterative formulae for structure node coordinates were obtained. The research results show that the present method has good stability, high precision, quick convergence and easiness in use, and is very efficient for solving the large displacement problems of trusses.

Key words: nonlinear analysis; Newton tangent method; large displacement; finite element method; planar flat truss; node

0 引 言

在杆系结构的设计分析中,结构通常假设为桁架(节点铰接结构),杆件轴线为直线,杆件材料为线弹性状态。当结构位移较小时,可忽略位移对结构平衡方程的影响,并通常以结构零态构形建立结构平衡方程。当结构位移较大时(如高耸、大跨及柔性结构受荷载作用后[HJ]其位移通常已达到与结构尺度同一量级的情况),位移对结构平衡方程的影响不能忽略,需要基于荷载态构形建立结构非线性平衡方程,且需要采用非线性求解方法来求解。

桁架大位移问题非线性分析涉及以下3个方面的问题:①用何变量作为结构非线性平衡方程的未知量;②采用何种方法建立结构非线性平衡方程;③利用何种方法求解结构非线性平衡方程。

结构非线性平衡方程的未知量常用的有2种:节点位移和节点坐标。传统上,人们通常取节点位移为未知量建立结构非线性平衡方程，但是近年来一些研究则选取节点坐标作为未知量［13］。Liu等［4］研究也发现,以节点坐标为未知量建立的结构非线性平衡方程列式简单,使求解方法也变得简单。

结构非线性平衡方程的建立方法常用的有2种:直接建立方法［5］和势能最小化方法［13］。直接建立方法是根据杆件内力和外荷载直接建立总节点力向量平衡方程,该方法概念明确，应用简单。势能最小化方法是：先建立结构总势能泛函,将非线性几何方程和应力应变关系代入总势能泛函中,经离散化处理得到关于未知量的泛函表达式,再根据总势能最小化条件得到关于未知量的非线性平衡方程。势能最小化方法也常忽略一些无法处理的高阶项。

结构非线性平衡方程的求解方法可分为2种主要类型:直接迭代法和微分线性逼近法［6］。直接迭代法的特点是:不涉及求导过程,只是利用原非线性方程建立迭代公式,它适用于不能或很难求导的非线性方程。该方法适用性较强,但是收敛速度慢。微分线性逼近法的特点是:求出结构非线性方程关于未知量的导数矩阵(Jacobi矩阵),据此建立原方程的线性逼近方程,得到关于未知量的迭代公式，进而求解。目前建立桁架结构非线性平衡方程时习惯于以节点位移作为未知量,很难实现关于节点位移未知量的求导过程,并且也没有相关的研究报道。

对于目前桁架大位移问题的分析方法,孙焕纯等［7］和许强等［8］认为“以往对桁架结构的大变形非线性分析,都是应用最小势能原理建立关于节点位移的非线性联立平衡方程,求解的工作量大,尤其对多自由度的大型复杂桁架更为突出”。为了解决这些问题, 文献［7］,［8］中提出了“交替迭代线性逼近法”,并有效应用于2个典型的桁架大位移问题,但是对于一些荷载工况也不能得到问题的收敛解。

为了求解桁架大位移问题,笔者提出一种新的非线性有限元分析方法：①以杆端节点坐标向量为显函变量写出单元杆端力向量表达式,由单元杆端力向量装配结构总节点力向量平衡方程;②根据矩阵微分理论,求出单元杆端力向量关于杆端节点坐标向量的导数矩阵,由单元导数矩阵装配结构关于总节点坐标向量的导数矩阵;③按照Newton切线法原理,建立结构等效线性逼近方程;④引入边界约束条件得到节点坐标的迭代公式。本文方法的结构非线性平衡方程列式简单且精确,求解过程也简单,收敛速度快,可有效求解桁架大位移问题。

1 结构非线性平衡方程的建立方法

设荷载态下单元的方向向量和长度分别为

νc(c)=［-I I］c

(1)

lc(c)=cTI66c

(2)

I66= I3 -I3-I3I3

(3)

式中:c为荷载态杆端节点坐标向量;νc(c)，lc(c)分别为单元的方向向量和长度;I为3阶单位向量。

荷载态下单元轴力为

N(c)=EAlc(c)-lz lz

(4)

式中:N(c)为荷载态下单元轴力;lz为零态单元长度;A为单元截面面积;E为材料弹性模量。

荷载态下单元杆端力向量f(c)为

f(c)=1 lc(c)-νc(c) νc(c)N(c)

(5)

将式(1)~(4)代入式(5),得到单元杆端力向量f(c)为

f(c)=EA(1 lz-1 cTI66c)I66c

(6)

由单元杆端力向量f(c),根据节点编号和对号入座原则,写出以总节点坐标向量为未知量的结构非线性平衡方程为

F(C)=f(c)-B=0

(7)

式中:F(C)为结构荷载态构形总节点坐标向量C的函数;B为结构荷载态总节点力向量;f(c)为由f(c)按节点编号和对号入座原则进行装配得到的向量。

2 结构非线性平衡方程的求解方法

令式(7)关于结构总节点坐标向量C的导数为FC(C),即

FC(C)=dF(C) dC=df(c) dc

(8)[HJ1.9mm]

根据矩阵微分理论,方程式(8)中df(c) dc的求导结果为

df(c) dc=EA (cTI66c)3/2I66ccTI66+EA(1 lz-1 cTI66c)I66

(9)

根据Newton切线法原理,结构非线性平衡方程的线性逼近方程式为

F(C(k))=F(C(k-1))+FC(C(k-1))• (C(k)-C(k-1))=0

(10)

根据节点坐标约束条件,C可分为非约束节点坐标Cf和约束节点坐标Cr。设有一元素为0或1的正交变换矩阵It可完成以下变换

CfCr=ItC或

C=ITtCfCr

(11)

将式(11)代入式(10),然后方程两边再左乘It,有

ItF(C(k-1))+ItFC(C(k-1))ITt•（CfCr(k)-CfCr(k-1)）=0

(12)

令P=ItF(C(k-1)),K=ItFC(C(k-1))ITt,则式(12)可简写为

P+K（CfCr(k)-

CfCr(k-1)）=0

(13)

对应Cf和Cr的分块,P和K也进行相应分块,对于约束节点坐标向量，有C(k)r-C(k-1)r=0，则式(13)可写为

PfPr+Kff KfrKrf Krr

C(k)f-C(k-1)f 0=0

(14)

按矩阵乘法展开式(14),得到节点坐标向量的迭代公式为

C(k)f=C(k-1)f-（Kff）-1Pf

(15)

3 算例分析

一平面扁桁架及其节点、杆件编号与节点荷载P见图1。扁桁架杆件截面均为圆形,杆件截面直径见表1，弹性模量E=200 GPa。扁桁架零态构形节点x,y坐标见表2，取右半部的结构进行分析。

图1 平面扁桁架（单位:cm）

Fig.1 Planar Flat Truss (Unit:cm)

采用本文方法对1.0P荷载系（图1）进行分析,节点位移计算结果见表3,杆件轴力计算结果见表4。为便于比较,文献［7］和ANSYS的计算结果也列于表3，4。通过比较可以看出,本文结果与文献［7］及ANSYS计算结果是一致的。

1.0P荷载系节点y方向位移和杆件轴力迭代[CM)][HJ]

表1 扁桁架杆件截面直径

Tab.1 Diameters of Member Section for Flat Truss

杆件编号 1,2,4,8,9 3,5,7,10 6

]截面直径/cm 3.8 2.0 5.2

表2 扁桁架零态构形节点x,y坐标

Tab.2 Node Coordinates in x，y Directions for Flat Truss

cm

节点编号 x坐标 y坐标

1 0.000 000 00 50.000 000 00

2 1 000.000 000 00 0.000 000 00

3 300.166 574 15 44.997 224 54

4 299.833 425 85 35.002 435 47

5 700.166 574 15 24.997 224 54

6 699.833 425 85 15.002 775 41

曲线分别如图2，3所示。从图2，3可以看出,迭代仅4次就满足了节点坐标容差ε=10-7的要求,文献［7］中迭代20次达到相同结果,本文方法收敛速度较快。

由图2，3还可以看出,迭代收敛前的曲线为较长一段水平线,说明迭代收敛解是稳定和可靠的。

对2.5P荷载系进行分析,通过对比可以看出,本文结果与文献［7］及ANSYS计算结果基本是一致的,只是在小数点后3位后有一些差别，这是由于方法不同产生的误差。

2.5P荷载系节点y方向位移和杆件轴力迭代变化曲线分别如图4，5所示。迭代8次就满足了节点坐标容差ε=10-7的要求,文献［7］中迭代20次达到相同结果,本文方法收敛速度较快。由图4，5可以看出,迭代收敛前的曲线为较长一段水平线,说明迭代收敛解是稳定和可靠的。

对2.6P荷载系进行分析,对于2.6P荷载系,文献［7］方法和ANSYS均不能求解。本文2.6P荷载系下节点位移和杆件轴力迭代曲线如图6所示。迭代20次满足了节点坐标容差ε=10-7的要[CM(22]求[KG-*7]。[KG*3]由图6可以看出,迭代收敛前的曲线为较长一

表3 节点位移

Tab.3 Node Displacements

cm

1.0P荷载系 2.5P荷载系 2.6P荷载系

节点位移 本文方法 文献［7］方法 ANSYS 本文方法 文献［7］方法 ANSYS 本文方法

v1 -1.046 800 -1.047 400 -1.047 -1.620 90 -1.619 -1.619 -103.900 0

u3 -0.018 554 -0.018 563 -0.019 -0.149 12 -0.150 -0.150 -0.739 2

v3 -0.920 740 -0.921 225 -0.921 -4.583 30 -4.607 -4.606 -77.410 0

u4 -0.013 704 -0.013 711 -0.014 -0.180 38 -0.182 -0.182 0.255 5

v4 -0.920 910 -0.921 390 -0.921 -4.582 30 -4.605 -4.605 -77.390 0

u5 -0.020 650 -0.020 660 -0.021 -0.218 67 -0.220 -0.220 -0.837 6

v5 -0.502 240 -0.150 250 -0.503 -3.973 10 -3.998 -3.998 -35.840 0

u6 -0.004 703 -0.004 705 -0.005 -0.155 15 -0.156 -0.156 0.243 5

v6 -0.502 780 -0.503 030 -0.503 -3.975 00 -4.000 -4.000 -35.820 0

注:ui,vi分别为节点i在x,y方向的位移。

表4 杆件轴力

Tab.4 Axial Forces of Members

N

1.0P荷载系 2.5P荷载系 2.6P荷载系

杆件轴力 本文方法 文献［7］方法 ANSYS 本文方法 文献［7］方法 ANSYS 本文方法

N1 -155 84.000 -15 584.300 -15 584.502 -64 293.400 -64 485.182 -64 484.939 -5 568.84

N2 -150 75.000 -15 075.100 -15 074.934 -13 266.800 -13 071.182 -13 071.001 77 937.50

N3 77.098 77.102 77.104 359.310 358.438 358.449 -1 301.68

N4 -128 98.000 -12 898.300 -12 898.483 -56 338.500 -56 521.954 -56 521.950 -9 621.14

N5 -2 704.200 -2 704.200 -2 704.233 -8 006.830 -8 015.200 -8 015.164 4 140.54

N6 -123 78.000 -12 378.300 -12 378.148 -4 942.960 -4 735.709 -4 735.696 73 808.40

N7 77.213 77.217 77.220 351.755 350.715 350.723 -1 304.25

N8 -156 28.000 -15 627.900 -15 628.101 -64 740.200 -64 935.253 -64 935.106 -5 420.73

N9 -150 80.000 -15 080.400 -15 080.215 -12 940.600 -12 741.762 -12 741.534 77 929.50

N10 -2 696.400 -2 696.400 -2 696.384 -8 321.230 -8 332.917 -8 332.953 4 105.19

注:N1~N10分别为1~10号杆件的轴力。

图2 1.0P荷载系下节点y方向位移迭代曲线

Fig.2 Iteration Curves for Node ydirection Displacements Under 1.0P Load System

图3 1.0P荷载系下杆件轴力迭代曲线

Fig.3 Iteration Curves for Axial Forces of Members Under 1.0P Load System

图4 2.5P荷载系下节点y方向位移迭代曲线

Fig.4 Iteration Curves for Node ydirection Displacements Under 2.5P Load System

图5 2.5P荷载系下杆件轴力迭代曲线

Fig.5 Iteration Curves for Axial Forces of Members Under 2.5P Load System

图6 2.6P荷载系下节点y方向位移迭代曲线

Fig.6 Iteration Curves for Node ydirection Displacement Under 2.6P Load System

图7 2.6P荷载系下结构平衡状态

Fig.7 Equilibrium State of Structure Under 2.6P Load System

段水平线,说明迭代收敛解是稳定和可靠的。

2.6P荷载系下结构平衡状态如图7所示。对大于2.6P荷载系也进行了分析,均能够获得稳定收敛解(分析结果略)。

4 结 语

(1)以节点坐标为未知量列出桁架大结构的非线性平衡方程,建立方程过程没有任何近似,所得到的非线性平衡方程是精确的列式。

(2)基于本文的桁架结构非线性平衡方程列式,可有效利用微分线性逼近这种精确方法进行求解,所形成的求解方程很简单,易于应用。

(3)大位移平面扁桁架算例研究表明,本文方法对于结构非常接近失稳的稳定平衡状态能够获得精[LL][LL]确解,且具有很高的稳定性和很快的收敛速度。

(4)本文方法稳定性好、精度高、收敛速度快且简单易用,为求解桁架大位移问题提供了一种有效方法,具有重要理论意义和工程应用价值。

参考文献：

References:

［1］OHKUBO S,WATADA Y,FUJIWAKI T.Nonlinear Analysis of Truss by Energy Minimization［J］.Computers & Structures,1987,27(1):129145.

［2］ KAZBERUK A,MIEDZIALOWSKI C Z,TRIBILLO R.Finite Element Discretization by Minimization of Elastic Strain Energy Method［J］.Finite Elements in Analysis and Design,1999,32(2):6370.

［3］ GRECOA M,GESUALDOA F A R,VENTURINIB W S.Nonlinear Positional Formulation for Space Truss Analysis［J］.Finite Elements in Analysis and Design,2007,42(12):10791086.

［4］ LIU S T,LONG Q L.A New Method Tracing Loaddeflection Equilibrium Path of a Doubly Nonlinear Truss［J］.Applied Mechanics and Materials,2012,166169:6872.

［5］ 沈世钊,徐崇宝,赵 臣,等.悬索结构设计［M］.2版.北京:中国建筑工业出版社,2006.

 SHEN Shizhao,XU Chongbao,ZHAO Chen,et al.Design of Cable Structures［M］.2nd ed.Beijing:China Architecture & Building Press,2006.

［6］ 熊洪允,曾绍标,毛云英.应用数学基础［M］.天津:天津大学出版社,1994.

 XIONG Hongyun,ZENG Shaobiao,MAO Yunying.Fundamentals of Applied Mathematics［M］.Tianjin:Tianjin University Press,1994.

［7］ 孙焕纯,许 强,龙武智.桁架结构几何大变形分析的精确方法［J］.应用力学学报,2009,26(1):4550.

 SUN Huanchun,XU Qiang,LONG Wuzhi.Accurate Algorithm for Geometrically Large Deflection Analysis of Truss Structures［J］.Chinese Journal of Applied Mechanics,2009,26(1):4550.

［8］ 许 强,陈 庆,孙焕纯.大跨度桁架几何大变形结构分析的一种数值方法［J］.土木工程学报,2009,42(1):1622.

 XU Qiang,CHEN Qing,SUN Huanchun.An Algorithm for Analysis of Large Span Truss Structures with Large Deformation［J］.China Civil Engineering Journal,2009,42(1):1622.