空间网格结构几何非线性有限元分析方法的研究

收稿日期:2001-01-20;修改稿收到日期:2001-07-17.作者简介:董石麟(1932-),男,教授,中国工程院院士.
第19卷第3期2002年8月
　计算力学学报　
Chinese Journal of Computational Mechanics
V ol .19,N o .3A ug ust 2002
文章编号:1007-4708(2002)03-0365-04
空间网格结构几何非线性有限元分析方法的研究
董石麟1,　张志宏1　,李元齐2
(1.浙江大学土木系空间结构研究中心,浙江杭州310027;2.同济大学土木系上海200092)
摘　要:在对T .L /U .L .法和C .O ra n 梁柱单元有限元法进行系统研究比较的基础上,推导了结合以上两种理论的几何非线性有限元列式。

C.O ran 梁柱单元在分析轴力为主的结构时具有非常高的效率,但在分析纯弯曲问题时却存在很大困难。

本文采用V C++6.0编制了面向对象的非线性有限元程序,对其进行了考题验算分析,得出了有用的结论。

关键词:几何非线性有限元;修正的拉格朗日描述;梁柱单元中图分类号:T U 311 文献标识码:A
1　引　言
空间网格结构的非线性有限元静力分析主要有两种方法,一是基于梁柱理论有限元法
[1-4]
;二是
基于T L/UL 描述的有限变形理论有限元方法[5-8]。

方法一在分析轴力为主的结构时具有非常
高的效率,但在分析纯弯曲问题时却存在很大困难。

本文推导了基于U.L.描述的梁单元切线刚度矩阵,单元内力求解分别采用两种方法,使得本文方法可以高效求解以轴力为主的结构,而且在分析以弯矩为主的结构仍能保持内力的精确求解。

最后,用VC ++ 6.0编制了采用本文方法的空间结构非线性有限元程序。

2　基于T .L /U .L .描述的几何　非线性梁单元
2.1　全量型变分原理
材料类型为超弹性材料,同各文献一样采用全量型变分原理: ∏=
∫
　0V
S ij (k ) ij 0
d V -∫
　0
V
　0F i
0u i 0
d V -∫
　0
S
　0T i
　0u i 0
d S =0(1)
式(1)是在Lag rang e 描述下的即以初始构形为参
考构形,若采用于Updated Lagrange 描述下则以上各量应基于t 时刻构形[5-7]即:
∏=∫　t
V t S ij
(k ) t ij d t V -∫
　t
V
　t F i
t u i t
d V -∫　t
S
t
T i
　t u i
t
d S =0(2)
其中,
∏为系统能量泛函;t
S
ij (k )
为基于t 时刻构形
的第二Piola-Kirchhoff 应力;t ij 为基于t 时刻构形
的Green 应变;t F i 、t T i
分别为t 时刻单位体积中的
体力及t 时刻单位边界上的面力;t V 、t
S 分别代表t 时刻的体积和边界面。

考察t + t 时刻的构形,由式(1)可以推导出T.L.描述的单元切/割线刚度矩阵有限元列式,由式(2)可推导出U.L.描述下单元切线刚度矩阵,具体推导请参考文献[5-8]。

2.2　坐标变换和内力计算
(1)坐标变换J.Argy ris
[9]
提出的描述大转动的转换方法,
该方法非常有效,本文采用了该方法。

此外还有文献[1,3,4]节点定向矩阵的方法以及文献[6]中的坐标变换方法等。

(2)内力计算
内力的计算是影响计算结果正确性的主要因素,内力的求解将直接影响残差的大小,而残差大小决定了近似解逼近真解的程度。

本文采用C.
Oran 梁柱理论计算内力,这一点与文献[5,6]不同,对于纯弯曲结构或以弯矩为主的结构,本文程序采用T.L.法中的单元割线刚度矩阵来求解。

用Saafan 关于杆件挠曲对轴向变形的影响以及Liveslay 和Chandler 推导的稳定函数,C.Oran 运用传统的梁柱理论得出梁单元在随动局部坐标系下的杆端力与位移的关系[3,4],其中单元轴力的计算公式为
N =E A (
u
L 0
-C b 2-C b 3)(3)
式中,E 、L 0、A 分别为材料的弹性模量、单元初始长度、单元截面面积;N 为单元轴力(压为正);u 、C b 2和C b 3分别为单元轴向缩短量、单元由弯曲变形引起的轴向变形。

梁柱理论有限元列式的具体推导可参考文献[1,3,4]。

3　C .Oran 梁柱理论有限元方法与　有限变形理论有限元方法的比较
C .Oran 梁柱单元在计算轴力为零或以弯矩为主的结构或构件时非常困难,而在分析以轴力为主的结构时基于有限变形理论的有限元的效率不如C.Oran 梁柱理论有限元,例如算例2要达到相同的精度,梁柱理论有限元只需要最少的单元数为二,而有限变形理论的有限元则需要十个单元[5]。

下面从两个方面分析其原因。

3.1　理论假设与数值模型
划分单元之后的结构即成为一个数值模型,几何非线性分析主要考虑大位移的影响,如果单元两节点的相对位移较大,用有限变形理论有限元,无论是T .L .法还是U .L .法,割线刚度矩阵还是切线刚度矩阵,位移的高阶项均已省略,故随着节点相对位移的增大,误差会增大。

如果单元数目增加一倍。

则左右节点相对位移的值要减少一半即 /2(如图1),位移高阶项的影响衰减很快,相应的误差也会小很多,即弱化了结构刚体位移(大位移)的影响,这里的刚体位移指32节点 单元的刚体位移,即不引起变形的位移,而对1-2节点,作为一个单元来讲 为引起变形的位移,故在每一增量步就更接近于其假设每一增量步仍为小变形小应变的假设,因此精度也越高。

此外,U.L.法只能推导出单元切线刚度矩阵,如没有单元内力的精确求解,误差会越积越大,最终出现应力漂移。

其二,基于有限变形理论的有限元如文献[6]K.J.Bathe 梁元的坐标变换是针对直梁单元,如果梁的变形为曲线则该坐标变换是不精确的,同一结构若划分的单元越多,则梁的变形就越接近小变形,则该坐标变换引起的误差也将越小。

其三,从数学的角度来讲,划分的单元越多,位移函数空间就越完备,因此数值模型就越接近实际结构。

理论假设与数值模型不符是有限变形理论有限元较少的单元不能实现高精度的原因。

虽然如果
单元较少可以通过每一增量步取得很小来逼近小
变形小应变,对此作者也对算例2作过分析,精度提高的效果并不明显。

3.2　理论假设与实际结构
梁柱理论来源于T im oshenko 所提出的稳定函数理论。

如果轴力为零,半解析C .Oran 梁柱理论应用于悬臂梁的内力轴力计算误差很大,因此在计算悬臂梁的纯弯曲时导致其收敛非常困难。

因此,梁柱理论的缺陷是理论上的,主要表现在:
首先,C .Oran 梁柱单元计算内力的公式及由其得到的增量切线刚度矩阵,均是在平面压弯梁的解析解的基础上迭加而得到的,并不是空间梁真实的变形曲线,而空间梁真实变形曲线的推导在弹性力学中还没有得到完善解决。

其次, C.Or an 梁柱单元在计算轴力时(式(3)),不能精确考虑弯矩引起的梁弦长的缩短,从而在纯弯曲悬臂梁中出现了轴力。

可以通过迭代或数值积分来求梁变形后弧长长度,前提是必须首先得到梁的解析变形曲线。

也可以对式(3)中轴力的计算另行推导,在求梁弧长时考虑四阶或更高阶项,具体推导可参考文献[10]。

再者, C.Oran 梁柱单元在推导刚度矩阵时忽略了剪切、扭转及弯剪、剪扭耦合应变能,该刚度矩阵是不完全的。

文献[8]指出梁柱单元刚度矩阵只是有限变形理论有限元刚度矩阵在引入一些假定之后得到的显式,因此其适用范围要小一些。

4　算例分析
算例1　悬臂梁
图3(a)为由图1悬臂梁纯弯曲采用本文方法
20个单元的分析数据形成的曲线,该结果与文献[1]非常吻合,而梁柱理论有限元对该问题完全失效。

图3(b)、(c)、(d)为梁柱理论有限元与Ansys5.5分析结果的对比,可以发现一个非常有趣的现象,即梁柱理论有限元当轴力很小时误差也很小。

从以上分析可以看出,对该模型而言梁柱理论分析结果与Ansy s5.5基本吻合,但收敛性随着偏心距的增大而越来越差,在纯弯情况下则非常困难。

366
计
算力学学报
　第19卷　
算例2　威廉平面刚架
威廉平面刚架是几何非线性分析中的经典算例,其基本尺寸如图4(a)所示。

许多文献[1,5,10]对此结构做过分析。

采用本文方法(两单元)进行了分析。

该结果与文献[1,5,10]的结果是完全吻合的。

而有限变形理论有限元如文献[5,8]则是采用10个单元。

5　结　论
C.Oran梁柱单元在求解纯弯曲问题时非常困难,而经典的T.L./U.L.法在求解以轴力为主的结构时效率不如梁柱理论有限元高。

反过来,求解以弯矩为主的结构时,虽然梁柱单元误差也可以很小,但是其计算效率将随轴力与弯矩比值的减小而降低,甚至不收敛。

因此,在不了解结构杆件内力分布的情况下采用传统的梁柱单元进行分析,存在一定的问题,有待进一步的研究。

参考文献(References):
[1]　李元齐.大跨度拱支网壳稳定性研究[D].上海:同济
大学,1998.(L i Y ua nqi.Stability r esearch on lar ge-
span a rch-suppor ted ret iculated shell structures[D].
Shanghai:T o ngji U niv ersit y,1998.(in Chinese)) [2]　沈士钊,陈　昕.网壳结构稳定性[M].北京:科学出
版社,1999.(Shen Shizhao,Chen Xin.T he S tability
of L attice S hells[M].Beijing:Science Publishing,
1999.(in Chinese))
[3]　Cena p Or an and A SCE M.T angent st iffness in plane
fr ames,ASCE[J].J our nal of the Str uctur al
367
　第3期董石麟,等:空间网格结构几何非线性有限元分析方法的研究
D iv is ion,1973,6:973-985.
[4]　Cenap O ran and ASCE M.T ang ent stif fness in space
fr ames[J].A SCE.J our nal of the Structur al
D iv is ion,1973,6:987-1001.
[5]　朱忠义.球面组合网壳结构的几何非线性分析[D].
杭州:浙江大学,1995.(Zhu Z ho ngy i.A nalysis of
geo metr ically no n-linear stability o f co mposite lattice
spher ical shell[D].Zhejiang U niv ersity,1995.(in
Chinese))
[6]　Bat he K J and Said Bo lo ur chi.L ar ge displacement
analy sis o f t hree-dimensio nal beam str uctur es[J].
I nt.J.f or N um.M eth.in Engin,1979,14:961-986.
[7]　董永涛,张耀春.板壳结构非线性屈曲分析的修正拉
格朗日法[J].土木工程学报,1997,30(2):34-41.
(D ong Y o ng tao,Zhang Yao chun.U L fo rm u-lat ion
fo r no nlinear buckling analysis o f plates and shells
[J].J.Civ il Eng ineer ing,1997,30(2):34-41.(in
Chinese))
[8]　陈政清,曾庆元,颜全胜.空间杆系结构大挠度问题内
力分析的U L列式法[J].土木工程学报,1992,25
(5):34-44.(Chen Zhengqing,Z eng Q ing yuan,Y an
Q uansheng.A U L for mulat ion fo r inter nal for ce
analysis o f spatial fr ame str uctur es w ith lar g e
displacement[J].J.Civil Engineer ing,1992,25(5):
34-44.(in Chinese))
[9]　Ar g yr is J.A n ex cur sion into lar ge ro tatio ns[J].
Comp.M eth.A pp l.M ech.E ng,1982,32:85-155. [10]Saafan S A.N onlinear behav io ur of st ructural plane
fr ames[J].J ournal of the S tructur al D ivision,
A SCE,1963,89(4):557-579.
Research on geometrical nonlinar finite element
method of spatial reticulated structures
Dong Shilin1,　Zhang Zhihong1,　Li Yuanqi2
(1.Depar tment o f Civil Eng ineering,Z hejiang U niver sity,Zhejiang.Hang zho u310027,China;
2.Depar tment o f Civil Eng ineering,T ongji U niver sity,Shang hai200094,China)
Abstract:Based on systemic study o n the geometrical nonlinear finite elem ent theory(C.Oran beam-co lum n theor y and T.L.m ethod/U.L.method)o f spatial structures,A g eo metr ical nonlinear finite elem ent method and som e discussion is presented in this paper,a new phenom non that C.Or an beam-co lum n elem ent can not solv e cantiliever only bearing monment is fir stly put fo rw ard and discussed.A nonlinear finite elem ent prog ram has also been dev elo ped by ado pting VC++6.0language.
Key words:g eo metr ical no nlinear finite element method;updated Lagrange descriptio n;beam-column element
368计算力学学报　第19卷　。