大型柔性空间结构热_非线性动力学耦合有限元分析

CN11-2223/NJTsinghuaUniv(Sci&Tech),2008,Vol.48,No.2w29

http://qhxbw.chinajournal.net.cn

大型柔性空间结构热-非线性动力学耦合有限元分析

段　进,　薛明德,　向志海

(清华大学工程力学系,北京100084)

摘　要:该文从三维固体力学基本方程和虚功原理出发,考虑几何非线性和截面翘曲的影响,发展了一种热-动力学耦合分析薄壁杆单元,并采用Newmark方法和Newton-Raphson迭代联合求解。其大转动问题采用Rodrigues公式考虑,其截面温度分解为平均温度和摄动温度。该文将瞬态温度场分析和几何非线性很好地统一起来,大大降低了计算规模,为复杂柔性空间结构的热-动力学分析提供了一套求解方案。通过与文献中的数值算例比较,验证了该方法的可靠性。利用该方法分析哈勃望远镜的太阳翼,解释了其扭转振动和动力屈曲等问题。

关键词:空间结构;热-动力学分析;大转动;几何非线

性;翘曲中图分类号:V414.1文章编号:1000-0054(2008)02-0276-04

文献标识码:A

Coupledthermalandnonlineardynamicsfiniteelementanalysisofflexiblelarge

spacestructures

DUANJin,XUEMingde,XIANGZhihai

(DepartmentofEngineeringMechanics,TsinghuaUniversity,

Beijing100084,China)Abstract:Aconsistentbeamelementwasdevelopedforacoupledthermallyanddynamicsfiniteelementanalysisfromsolidmechanicstheoryandthebasicprinciplesofvirtualwork.Thegeometricnonlinearitiesandsectionalwarpingaretakenintoaccount.LargerotationisconsideredbythestructuralelementusingRodriguesformula.Besides,thenodaltemperaturevariablesaredecomposedintoboththeaverageandtheperturbationtemperatures.ThenonlineargoverningequationissolvedusingtheNewmarkmethodandNewton-Raphsoniteration.Atransienttemperatureanalysisisconsolidatedwithageometricnonlineardisplacementanalysisbythebeamelementforanalyzingthecouplednonlinearthermalanddynamicsresponseoflargeflexiblespacestructures.Thevalidityofthemethodisverifiedbycomparingwithwell-knownnumericalsimulationsintheliterature.AthermalanddynamicsanalysisoftheHSTsolararrayispresentedtodemonstratetheperformanceofthescheme.

Keywords:spacestructures;thermallydynamicanalysis;large

rotation;geometricnonlinearity;warping

坏。Thornton和Kim对该太阳翼进行热诱发振动理论分析,较好地解释了该结构在进出地球阴影时,由

[2]

于热流条件的突变容易诱发弯曲振动。Xue和Ding采用有限元方法,提出了一种考虑辐射的闭口薄壁温度杆单元,轴向温度采用两节点线性插值,截面内温度分解成平均温度和摄动温度,并构造周向插值函数使

[3]

各谐温度解耦。Xue、Duan等将该温度单元推广到开口薄壁截面,并导出了小变形情况下的热-动力学耦合方程,利用该单元对哈勃望远镜太阳翼进行热诱发振动分析,其结果显示:进出地球阴影时,太阳翼会发生明显的弯曲和扭转振动。

文[1-4]均以小变形假设为前提,而航天器的附件大多为柔性构件,应该考虑几何非线性的影响。本文考虑大转动的影响,采用更新的Lagrange格式,推导Euler-Bernoulli梁几何非线性热-动力学方程,并采用Newmark方法和Newton-Raphson迭代联合求解。梁单元温度场不仅包含传统的平均温度,还包括截面温差[2-3]。梁单元的大转动分为侧向转动和绕轴线的转动,前者采用Rodrigues公式计算,后者是平面转动。

本文将文[2-3]的温度场单元应用于几何非线性动力学分析,将薄壁杆的瞬态温度场和几何非线性统一起来,采用统一的梁单元模型求解其非线性热-动力学响应,以期降低计算规模,分析大型柔性空间结构的耦合热诱发振动问题。

[1][1]

1　当前构型下梁的几何非线性热-动力学

方程

本文采用的单元模型是2节点Euler-Bernoulli梁单元,满足如下假设条件:1)刚周边假定;2)大转动,小应变;3)各向同性线弹性材料。

在几何非线性分析中,梁单元的位移可分解为刚体位移和弹性变形2部分迭加。前者为几何非线性分析提供参考构型,后者为其提供每一个迭代步所需的预应力值。本文采用更新的Lagrange格式,以t时刻的梁构型(见图1)为参考构型,建立梁单元的平衡方程,以求解t+󰀁t时刻的响应。由虚功原

收稿日期:2007-03-14

基金项目:国家“八六三”高技术项目(863-2-2-1-9)作者简介:段进(1978—),男(汉),湖南,博士研究生。通讯联系人:薛明德,教授,

E-mail:xuemd@mail.tsinghua.edu.cn开口和闭口薄壁杆件大量应用于航天结构的柔性

附件中,比如太阳翼和天线等,在进出地球阴影时,由于热流条件的突变,其温度场会发生变化,从而容易诱发该类结构的热-动力学问题。最有名的如1990年,哈勃太空望远镜发射升空后,其太阳翼发生明显的弯曲和扭转振动,绕地球几个周期后,结构发生屈曲破段　进,等:　大型柔性空间结构热-非线性动力学耦合有限元分析

277

[3]

理得

∫󰀁

0

l

[　　t󰀂xx t!xx+

t+󰀁t

t+󰀁t

　　txy

󰀂 t∀xy+

t+󰀁t

󰀂 ∀]dydzdx=(1)R.

t+󰀁tt+󰀁t

其中:　　t󰀂、󰀂和　　t󰀂xz表示t+󰀁t时刻的轴向应力(包括热应力)和剪应力;t!xx、t∀xy和t∀xz表示由t至t+󰀁t时刻的增量应变(包括线性和非线性部分);t+󰀁t

R表示外载(包括惯性力)在该时间步内的虚功。式(1)积分得

tte*~t+󰀁t󰀁󰀁e*

Mta+(tKL+tKNL)ta=

t+󰀁tt+󰀁tt

(2)　　tPe+　　tPT-tP.

*

+󰀁t󰀁~表示梁单元集中质量阵;t　　󰀁e

其中:Mta表示t+󰀁t时刻梁节点的加速度;tKL表示梁单元线性刚度阵;te*

tKNL表示梁单元几何非线性刚度阵;ta表示梁节

t+󰀁t

点的增量位移;　　tPe表示t+󰀁t时刻梁节点所受的

t+󰀁t

外载;　　tPT表示t+󰀁t时刻梁节点的等效温度载

t

荷,详见文[3];tP表示t时刻的节点内力;符号*表示扭转中心。

t

At+󰀁t

　　txztxz

t+󰀁txx　　txy

可描述截面内任意点的翘曲值)。对其作坐标转换

时,忽略翘曲自由度的转换,其转换矩阵可表示为

t-0R

t-0RO

t-1

.(5)0R=t-0R

t-O0R

1t-000

其中:0R为3×3矩阵,它表示坐标系(x,y,z)与ttt

(x,y,z)的转换关系;O为零矩阵。

梁单元局部坐标通过刚体位移从0时刻构型转到t时刻构型,其过程可分为侧向转动和轴向转动2

000

步:1)坐标轴x、y和z绕空间某一轴线(沿该轴

tttt

~和t~线的单位矢量为e∃)转动∃角至x、yz,如图2a[6]

所示,该转动按Rodrigues公式计算,转动矩阵记tdttt为-R(∃),其关键在于求解空间角矢量∃;2)~y和ttttt~z绕坐标轴x转动%角至y和z,该过程是平面转

t-at

动,如图2b所示,转动矩阵记为R(%),其关键在于

t

求解%。

图1　开口薄壁梁单元运动图

　　上述方程通过3次坐标转换并组集可得结构总体方程:1)扭转中心转换到形心,转换矩阵记

[3]

为#;2)t时刻构型转换到0时刻构型,转换矩阵

t

记为0-R(详见第2节);3)0时刻单元局部坐标转换到总体坐标。结构总体的动力学方程为

t+󰀁t󰀁󰀁tt+󰀁ttMa+K󰀁a=Q-F.(3)

t+󰀁t

其中:󰀁a表示形心增量位移;Q表示外载(含等

t

效温度载荷);F表示结构内力。

采用Newmark方法和Newton-Raphson迭代

[5]

可将式(3)改写为t

K+42M󰀁a(k)=t+󰀁tQ-t+󰀁tF(k-1)-󰀁t

4t+󰀁t(k-1)t4t󰀁t󰀁󰀁

a-a)-a-aM2(.(4)󰀁t󰀁t

t+󰀁t(0)tt+󰀁t(0)t

其中:k表示迭代步,且F=F,a=a.

须指出的是,薄壁杆单元的温度分解为平均温

tt+󰀁t

度和摄动温度,反映在方程(4)的K和Q中,具体计算参考文[3]。

图2　梁单元大转动

t-初始构型到当前构型的转换矩阵0R可表示为

t-t-att-dt

R(%)R(∃).(6)0R=

　　如图2a所示,梁单元初始时刻的轴方向单位矢

0

量记为ex,是已知的,t时刻的轴方向单位矢量记t

为ex,它由梁节点原始坐标和节点位移确定,也是

0tt

已知的。由ex和ex可确定空间角矢量∃,即:

2　大转动情况下梁节点位移的转换关系

对于本文的开口薄壁梁单元,每个节点共有7个自由度,3个位移自由度,3个转角自由度,以及截面翘曲自由度(该自由度与截面翘曲函数的乘积278

0

t

t

0

t

清华大学学报(自然科学版)2008,48(2)

ex×ex

∃=[arccos(ex󰀁ex)]0t.(7)

󰀂ex×ex󰀂

[6]tdt

　　按Rodrigues公式,侧向转动矩阵-R(∃)可

表示如下:

t

t-dtTsin∃t

[R(∃)]=I+tS(∃)+

∃t

1-cos∃t2

t2S(∃).(8)∃

t

其中:I为3×3的单位矩阵;S(∃)为3×3反对称矩[6]

阵。

轴向转动是平面转动,Bathe和Bolourchi在文[7]中已给出详细表达式。

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]一致。

3　数值算例

下面先给出2个典型的几何非线性动力学算例,通过与已有文献结果比较,验证本文方法的可靠性,然后利用本文方法分析哈勃望远镜的太阳翼,并解释1990年著名的哈勃太阳翼失效事件。3.1　固支梁受集中力作用

两端固支梁中部受突加的集中力作用,杆长为L,横截面为矩形,宽和高分别为b和h,突加载荷P=2848N,几何尺寸和材料参数如图3所示。梁中部挠度w与时间t的关系曲线如图4所示,其中非线性解对应左侧坐标,线性解对应右侧坐标。

图4显示:1)线性分析的振幅远大于非线性分析振幅,其准静态变形(图中未显示)也有同样规律,这是因为该工况下非线性刚度远大于线性刚度;2)本文结果与文[8-9]的结果吻合较好,这在一定程度上证明了本文方法的可靠性。另外,由于非线性振动频率明显高于线性振动频率,计算时所取的时间步长也不相同,本算例中非线性分析时间步长为25&s,线性分析为200&s。

图5　简支曲梁模型图

3.3　哈勃太阳翼热-动力学分析

哈勃望远镜太阳翼模型如图7所示,太阳毯采用梁单元构成的网来模拟。太阳热流始终垂直表面入射,数值为1.350kW/m。太阳翼几何尺寸和材料

2

图3　固支梁模型图

图6　简支曲梁挠度参数时间响应图

图4　固支梁中部挠度随时间变化图

参数见文[1,3]。左梁和右梁均是开口薄壁截面,由

于展开机构的原因,无法保证开口方向的对称性,本文假设左梁开口朝Y方向,右梁由Y向Z轴偏转(段　进,等:　大型柔性空间结构热-非线性动力学耦合有限元分析

279

角。当(=0°时,太阳翼是左右对称结构,不会出现扭转现象,该情况等效于文[1]的闭口薄壁模型;当(≠0°时,由于左右开口薄壁杆的温度分布不同,其热-动力学响应也不同,结构会发生扭转振动:1)当(=10°时,左梁和右梁端部的挠度时间曲线如图8

所示,左右梁端部的挠度偏差约为0.2m;2)当(=

20°时,挠度时间曲线如图9所示,非线性分析结果显示哈勃望远镜的太阳翼此时已经发生动力屈曲,该结论跟实际观测的结果一致。

[1]

图7　哈勃太阳翼模型图

方法分析哈勃望远镜的太阳毯,得出了与实际观测一致的结论。

参考文献　(References)

[1]

ThorntonEA,KimYK.Thermallyinducedbendingvibrationsofaflexiblerolled-upsolararray[J].JSpacecraftRockets,1993,30(4):438-448.

[2]DingY,XueMD,KimJK.Thermo-structuralanalysisof

spacestructuresusingFouriertubeelements[J].ComputMech,2005,36(4):289-297.

[3]XUEMingde,DUANJin,XIANGZhihai.Thermally-inducedbending-torsioncouplingvibrationoflargescalespacestructures[J].ComputationalMechanics,2007,40(4):707-723.

[4]程乐锦,薛明德.大型空间结构热-动力学耦合有限元分析

[J].清华大学学报(自然科学版),2004,44(5):681-684,688.

CHENGLejin,XUEMingde.Coupledthermal-dynamicFEManalysisoflargescalespacestructures[J].JTsinghuaUniv(Sci&Tech),2004,44(5):681-684,688.(inChinese)

[5]BatheKJ.FiniteElementProcedures[M].NJ:

Prentice-Hall,1996.

[6]CrisfieldMA.Non-LinearFiniteElementAnalysisofSolids

andStructures(AdvancedTopics)[M].Chichester:Wiley,1997.

[7]BatheKJ,BolourchiS.Largedisplacementanalysisof

three-dimensionalbeamstructures[J].IntJNumerMethEng,1979,14:961-986.

[8]YangTL,SaigalS.Asimpleelementforstaticanddynamic

responseofbeamswithmaterialandgeometricnonlinearities[J].IntJNumerMethEng,1984,20:851-867.

[9]MondkarDP,PowellGH.Finiteelementanalysisof

non-linearstaticanddynamicresponse[J].IntJNumerMethEng,1977,11:499-520.

[10]HumphreysJS.Ondynamicsnapbucklingofshallowarches

[J].AIAAJ,1966,4:878-886.

[11]BatheKJ,RammE,WilsonEL.Finiteelement

formulationsforlargedeformationdynamicanalysis[J].IntJNumerMethEng,1975,9:353-386.

[12]NoorAK,KnightF.Nonlineardynamicanalysisofcurved

beams[J].ComputMethApplMechEng,1979,23:225-251.图8　(=10°梁端部挠度图

图9　(=20°梁端部挠度图

4　结　论

本文考虑大转动的影响,导出薄壁杆结构几何

非线性热-动力学有限元方程。对于厚壁杆或者实心杆,可忽略截面翘曲或者选取特殊的翘曲函数,因此本文方法对于普通梁杆结构有一定的通用性。本文将薄壁构件的瞬态温度场计算与几何非线性动力学分析统一起来,大大降低了计算规模,提高了计算效率,可用来方便地分析大型柔性空间结构。用本文