清华大学学报(自然科学版)J T singh ua Un iv (Sci &Tech ),2008年第48卷第2期
2008,V o l.48,N o.2w 29
http://qhx bw.chinajo urnal.net.cn
大型柔性空间结构热-非线性动力学耦合有限元分析
段 进, 薛明德, 向志海
(清华大学工程力学系,北京100084)
收稿日期:2007-03-14
基金项目:国家“八六三”高技术项目(863-2-2-1-9)作者简介:段进(1978—),男(汉),湖南,博士研究生。通讯联系人:薛明德,教授,
E-mail:xuemd@mail.tsingh ua.edu.cn
摘 要:该文从三维固体力学基本方程和虚功原理出发,考虑几何非线性和截面翘曲的影响,发展了一种热-动力学耦合分析薄壁杆单元,并采用N ew mar k 方法和N ew ton -Raphson 迭代联合求解。其大转动问题采用Rodr igues 公式考虑,其截面温度分解为平均温度和摄动温度。该文将瞬态温度场分析和几何非线性很好地统一起来,大大降低了计算规模,为复杂柔性空间结构的热-动力学分析提供了一套求解方案。通过与文献中的数值算例比较,验证了该方法的可靠性。利用该方法分析哈勃望远镜的太阳翼,解释了其扭转振动和动力屈曲等问题。
关键词:空间结构;热-动力学分析;大转动;几何非线
性;翘曲中图分类号:V 414.1文献标识码:A
文章编号:1000-0054(2008)02-0276-04
Coupled thermal and nonlinear dynamics finite element analysis of flexible large
space structures
DUAN Jin ,XUE Mingde ,XIA NG Zhihai
(Department of Engineering Mechanics ,Tsinghua University ,
Beij ing 100084,China )Abstract :A consis tent beam elem ent was developed for a coup led ther mally and dynamics finite element an alys is from s olid m echan ics theory and th e basic principles of virtual w ork.Th e geom etric nonlin earities and s ectional w arping are tak en into accou nt.Large rotation is con sidered by the structural elem ent u sing Rodrigu es formula.Besides,the nodal temperature variables are d ecom pos ed into both the average and th e perturbation tem peratures.T he nonlin ear gover ning equation is s olved usin g the New mark meth od an d New ton-Raphs on iteration.A trans ient tem perature analys is is cons olidated w ith a geometric nonlin ear dis placement analys is b y the beam element for analyzin g the coupled n on linear th ermal and dynamics respon se of large flexible s pace structur es.T he validity of the meth od is verified by comparing w ith w ell-kn ow n num erical sim ulation s in the literature.A thermal and dynamics analysis of the HST solar array is pres ented to demons trate the performance of the schem e.
Key words :s pace s tructures;th erm ally dynamic analysis ;large
rotation ;geom etric nonlinearity;warping
开口和闭口薄壁杆件大量应用于航天结构的柔性
附件中,比如太阳翼和天线等,在进出地球阴影时,由于热流条件的突变,其温度场会发生变化,从而容易诱发该类结构的热-动力学问题。最有名的如1990年,哈勃太空望远镜发射升空后,其太阳翼发生明显的弯曲和扭转振动,绕地球几个周期后,结构发生屈曲破
坏[1]。Thornton 和Kim [1]
对该太阳翼进行热诱发振动理论分析,较好地解释了该结构在进出地球阴影时,由
于热流条件的突变容易诱发弯曲振动。Xue 和Ding [2]
采用有限元方法,提出了一种考虑辐射的闭口薄壁温度杆单元,轴向温度采用两节点线性插值,截面内温度分解成平均温度和摄动温度,并构造周向插值函数使
各谐温度解耦。Xue 、Duan 等[3]
将该温度单元推广到开口薄壁截面,并导出了小变形情况下的热-动力学耦合方程,利用该单元对哈勃望远镜太阳翼进行热诱发振动分析,其结果显示:进出地球阴影时,太阳翼会发生明显的弯曲和扭转振动。
文[1-4]均以小变形假设为前提,而航天器的附件大多为柔性构件,应该考虑几何非线性的影响。本文考虑大转动的影响,采用更新的Lagr ange 格式,推导Euler-Bernoulli 梁几何非线性热-动力学方程,并采用New mark 方法和New ton -Raphson 迭代联合求解。梁单元温度场不仅包含传统的平均温度,还包括截面温差[2-3]。梁单元的大转动分为侧向转动和绕轴线的转动,前者采用Ro drigues 公式计算,后者是平面转动。
本文将文[2-3]的温度场单元应用于几何非线性动力学分析,将薄壁杆的瞬态温度场和几何非线性统一起来,采用统一的梁单元模型求解其非线性热-动力学响应,以期降低计算规模,分析大型柔性空间结构的耦合热诱发振动问题。
1 当前构型下梁的几何非线性热-动力学
方程
本文采用的单元模型是2节点Euler -Bernoulli 梁单元,满足如下假设条件:1)刚周边假定;2)大转动,小应变;3)各向同性线弹性材料。
在几何非线性分析中,梁单元的位移可分解为刚体位移和弹性变形2部分迭加。前者为几何非线性分析提供参考构型,后者为其提供每一个迭代步所需的预应力值。本文采用更新的Lagrange 格式,以t 时刻的梁构型(见图1)为参考构型,建立梁单元的平衡方程,以求解t + t 时刻的响应。由虚功原
理得
∫l
A [t + t
t xx t !x x +
t + t
t xy
t ∀x y +
t + t
t xz t ∀xz
]d y d z d x =t + t
R .
(1)其中:t + t t xx 、t + t t x y
和t + t
t xz 表示t + t 时刻的轴向应力(包括热应力)和剪应力;t !x x 、t ∀x y 和t ∀x z 表示由t 至t + t 时刻的增量应变(包括线性和非线性部分);t + t
R 表示外载(包括惯性力)在该时间步内的虚功。式(1)积分得
M ~t + t t a e *
+(t t K L +t t K N L )t a e *=
t + t t P e +t + t t P T -t
t P .
(2)其中:M ~表示梁单元集中质量阵;t + t t a e
*
表示t + t 时刻梁节点的加速度;t
t K L 表示梁单元线性刚度阵;t t K N L 表示梁单元几何非线性刚度阵;t a e *
表示梁节
点的增量位移;t + t
t P e 表示t + t 时刻梁节点所受的
外载;t + t
t P T 表示t + t 时刻梁节点的等效温度载
荷,详见文[3];t
t P 表示t 时刻的节点内力;符号*表示扭转中心。
图1 开口薄壁梁单元运动图
上述方程通过3次坐标转换并组集可得结构总体方程:1)扭转中心转换到形心,转换矩阵记
为#[3]
;2)t 时刻构型转换到0时刻构型,转换矩阵
记为t
0R -(详见第2节);3)0时刻单元局部坐标转换到总体坐标。结构总体的动力学方程为
M t + t a +t K a =t + t Q -t F .(3)
其中: a 表示形心增量位移;t + t
Q 表示外载(含等
效温度载荷);t
F 表示结构内力。
采用New mark 方法和New ton -Raphson 迭代
可将式(3)改写为
[5]
t
K +4 t
2M a (k )=t + t Q -t + t F (k -1)-M 4 t
(t + t a (k -1)-t a )-4 t t a -t a
.(4)其中:k 表示迭代步,且t + t F (0)=t F ,t + t a (0)=t
a .
须指出的是,薄壁杆单元的温度分解为平均温
度和摄动温度,反映在方程(4)的t K 和t + t
Q 中,具体计算参考文[3]。
2 大转动情况下梁节点位移的转换关系
对于本文的开口薄壁梁单元,每个节点共有7个自由度,3个位移自由度,3个转角自由度,以及截面翘曲自由度(该自由度与截面翘曲函数的乘积
可描述截面内任意点的翘曲值[3]
)。对其作坐标转换
时,忽略翘曲自由度的转换,其转换矩阵可表示为
t 0R -=t 0R
-t 0R -O
1
t 0R
-O t 0R
-1
.(5)其中:t 0R -为3×3矩阵,它表示坐标系(0x ,0y ,0
z )与(t x ,t y ,t
z )的转换关系;O 为零矩阵。
梁单元局部坐标通过刚体位移从0时刻构型转到t 时刻构型,其过程可分为侧向转动和轴向转动2
步:1)坐标轴0x 、0y 和0
z 绕空间某一轴线(沿该轴
线的单位矢量为t e ∃)转动t ∃角至t x 、t
y ~和t z ~,如图2a 所示,该转动按Rodrigues 公式[6]
计算,转动矩阵记为t R -d (t ∃),其关键在于求解空间角矢量t ∃;2)t y ~和t z ~绕坐标轴t x 转动t %角至t y 和t z ,该过程是平面转
动,如图2b 所示,转动矩阵记为t R -a (t
%),其关键在于
求解t
%。
图2 梁单元大转动
初始构型到当前构型的转换矩阵t 0R -可表示为
t 0R -=
t R -a (t %)t R -d (t
∃).(6) 如图2a 所示,梁单元初始时刻的轴方向单位矢
量记为0
e x ,是已知的,t 时刻的轴方向单位矢量记为t
e x ,它由梁节点原始坐标和节点位移确定,也是
已知的。由0e x 和t e x 可确定空间角矢量t
∃,即:
277
段 进,等: 大型柔性空间结构热-非线性动力学耦合有限元分析
t
∃=[arccos(0
e x t
e x )]0
e x ×t
e x
0e x ×t e x
.(7)
按Rodr ig ues 公式[6],侧向转动矩阵t R -d (t
∃)可
表示如下:
[t R -d (t ∃)]T =I +sin t
∃t ∃S (t
∃)+
1-cos t
∃t ∃
2S (t ∃)2
.(8)其中:I 为3×3的单位矩阵;S (t
∃)为3×3反对称矩阵[6]
。
轴向转动是平面转动,Bathe 和Bo lourchi 在文[7]中已给出详细表达式。
3 数值算例
下面先给出2个典型的几何非线性动力学算例,通过与已有文献结果比较,验证本文方法的可靠性,然后利用本文方法分析哈勃望远镜的太阳翼,并解释1990年著名的哈勃太阳翼失效事件。3.1 固支梁受集中力作用
两端固支梁中部受突加的集中力作用,杆长为L ,横截面为矩形,宽和高分别为b 和h ,突加载荷P =2848N ,几何尺寸和材料参数如图3所示。梁中部挠度w 与时间t 的关系曲线如图4所示,其中非线性解对应左侧坐标,线性解对应右侧坐标。
图4显示:1)线性分析的振幅远大于非线性分析振幅,其准静态变形(图中未显示)也有同样规律,这是因为该工况下非线性刚度远大于线性刚度;2)本文结果与文[8-9]的结果吻合较好,这在一定程度上证明了本文方法的可靠性。另外,由于非线性振动频率明显高于线性振动频率,计算时所取的时间步长也不相同,本算例中非线性分析时间步长为25&s ,线性分析为200&s 。
图3 固支梁模型图
图4 固支梁中部挠度随时间变化图
3.2 简支曲梁受均布力作用
两端简支的圆弧曲梁受突加的分布法向载荷作用,截面为矩形,宽和高分别为b 和h ,几何尺寸和材料参数见图5。这是一个经典的动力屈曲问题,文[10-12]都曾经对它做过分析。该曲梁的挠度时间响应如图6所示(量纲1参数 和∋的具体意义见文[11]),当曲梁的载荷参数增加到一定数值时,结构会发生动力屈曲,其现象为:1)挠度参数 短时间内大幅度增加;2)平衡位置突然发生显著变化。图6显示,当载荷参数(其定义见图6)P =0.195时,曲梁绕平衡位置 =0.25附近振动,振幅约为0.5;当P =0.205时,平衡位置突变到 =1.6,振幅激增至3.2,结构显然已经屈曲。因此,该曲梁的屈曲点应该在P =0.195~0.205,略低于文[10]的计算结果P =0.210~0.215,略高于文[11]的P =0.190~0.20。当P =0.205时,本文的计算结果与文[12]的结果吻合,当P =0.25时,本文结果与文[10-11]的结果吻合。另外,观察图6
可知:当载荷参数P 达到屈曲载荷后,随着P 的增加,结构平衡位置和振幅均不会明显增加,但屈曲时间点会明显提前。该结论与文[10]一致。
图5 简支曲梁模型图
3.3
哈勃太阳翼热-动力学分析
哈勃望远镜太阳翼模型如图7所示,太阳毯采用梁单元构成的网来模拟。太阳热流始终垂直表面入射,数值为1.350kW/m 2
。太阳翼几何尺寸和材料
图6 简支曲梁挠度参数时间响应图
参数见文[1,3]。左梁和右梁均是开口薄壁截面,由
于展开机构的原因,无法保证开口方向的对称性,本文假设左梁开口朝Y 方向,右梁由Y 向Z 轴偏转(
278
清华大学学报(自然科学版)2008,48(2)
角。当(=0°时,太阳翼是左右对称结构,不会出现扭转现象,该情况等效于文[1]的闭口薄壁模型;当(≠0°时,由于左右开口薄壁杆的温度分布不同,其热-动力学响应也不同,结构会发生扭转振动:1)当(=10°时,左梁和右梁端部的挠度时间曲线如图8
所示,左右梁端部的挠度偏差约为0.2m ;2)当(=
20°时,挠度时间曲线如图9所示,非线性分析结果显示哈勃望远镜的太阳翼此时已经发生动力屈曲,该结论跟实际观测的结果[1]
一致。
图7
哈勃太阳翼模型图
图8 (=10°
梁端部挠度图
图9 (=20°梁端部挠度图
4 结 论
本文考虑大转动的影响,导出薄壁杆结构几何
非线性热-动力学有限元方程。对于厚壁杆或者实心杆,可忽略截面翘曲或者选取特殊的翘曲函数,因此本文方法对于普通梁杆结构有一定的通用性。本文将薄壁构件的瞬态温度场计算与几何非线性动力学分析统一起来,大大降低了计算规模,提高了计算效率,可用来方便地分析大型柔性空间结构。用本文
方法分析哈勃望远镜的太阳毯,得出了与实际观测一致的结论。
参考文献 (References )
[1]
T hornton E A,Kim Y K.Th ermally indu ced bending vibr ations of a flexible rolled-up solar array [J ].J S p acecr af t R oc kets ,1993,30(4):438-448.
[2]Ding Y,Xue M D,Kim J K.Th ermo-s tructural analysis of
s pace s tru ctures usin g Fou rier tube elements [J].Comp ut M e ch ,2005,36(4):2-297.
[3]XUE M in gde,DU AN J in,XIANG Zh ihai.T hermally-in duced bending-torsion coupling vibration of large scale s pace structures [J ].Comp utational M echanics ,2007,40(4):707-723.
[4]程乐锦,薛明德.大型空间结构热-动力学耦合有限元分析
[J ].清华大学学报(自然科学版),2004,44(5):681-684,688.
CHE NG Lejin,XUE M in gde.Coupled ther mal-dynamic FEM analysis of large s cale space s tru ctures [J ].J Tsinghua Univ (S ci &T ech ),2004,44(5):681-684,688.(in Ch ines e)
[5]Bathe K J.Finite Element Procedures [M ].NJ:
Prentice-Hall,1996.
[6]Cr isfield M A.Non-Linear Finite Element Analysis of Solids
and Structu res (Advanced T opics )[M ].Chiches ter:Wiley,1997.
[7]Bathe K J,Bolourchi S.Large dis placemen t analysis of
thr ee-d imens ional beam structu res [J ].I nt J N umer M eth E ng ,1979,14:961-986.
[8]Yang T L,S aigal S.A s imple element for s tatic and dynamic
respons e of b eams w ith material and geometric n onlinearities [J ].I nt J N umer M eth Eng ,1984,20:851-867.
[9]M on dkar D P,Pow ell G H.Fin ite elem ent analysis of
non -linear s tatic and d ynam ic respons e [J ].Int J N umer M e th Eng ,1977,11:499-520.
[10]Hu mphreys J S.On dynamic snap bucklin g of sh allow arches
[J ].A I AA J ,1966,4:878-886.
[11]Bathe K J,Ramm E,W ilson E L.Finite element
formulations for large deform ation dynamic analysis [J ].I nt J Nume r M eth E ng ,1975,9:353-386.
[12]Noor A K,Knight F.Nonlinear d ynam ic analysis of cu rved
beams [J ].Comput M e th A p p l M ech Eng ,1979,23:225-251.
279
段 进,等: 大型柔性空间结构热-非线性动力学耦合有限元分析