空间刚架动态失稳全过程分析

第18卷第 6期 工 程 力 学 Vol.18 No.6 2001年 12 月 ENGINEERING MECHANICS Dec. 2001 ——————————————— 收稿日期：2000-02-28；修订日期：2000-04-27 作者简介：王策(1968)，男，博士后，从事大跨度钢结构振动与稳定性研究 文章编号：1000-4750(2001)06-068-08 空间刚架动态失稳全过程分析 王 策，刘西拉 (清华大学土木系，北京 100084) 摘 要：采用静动力非线性方法研究结构的后屈曲问题。基于连续介质力学对物质运动的描 述方法，采用增量拉格朗日公式推导空间梁单元的切线刚度矩阵，包括节点大位移、大转角 产生的几何非线性，及弹塑性本构关系。动力分析采用 Newmark 数值积分，结合 Newton- Raphson平衡迭代。采用弧长法计算结构的屈曲前非线性反应，该方法不用卸荷跟踪荷载位移 反应全过程，而在临界点继续缓慢加荷，用动力时程方法进行空间刚架结构的失稳破坏全过 程仿真。 关键词：空间刚架；后屈曲；时程分析；失稳全过程仿真 中图分类号：TU318.2 文献标识码：A 1 概述 单层球面网壳为缺陷敏感性结构，在大雪、强风和强烈地震作用下，结构可能丧失整 体稳定性。跨度 93.5m 的布加勒斯特穹顶网壳结构，于 1963 年在一场持续长时间的大雪 后突然倒塌，这一事故使工程师们认识到网壳结构屈曲问题的重要性。目前关于球面网壳 的研究主要集中在结构静力稳定性及静力后屈曲分析，采用 Risk[1]提出的弧长法研究结构 跳跃屈曲，即在数值分析中通过对结构施加某种约束，求得结构荷载位移曲线的下降段。 而在实际结构中，外荷载往往不能自动卸荷，并且结构跳跃失稳过程为动态的，因此用静 力方法计算得到结构荷载位移曲线的下降段是虚假的。如果结构的后屈曲路径不稳定，荷 载随位移增加而下降，结构不能继续吸收外力功，而且要释放应变能，转化为驱动结构失 稳的动能。 Meek[2]用动力时程方法分析弹性刚架结构的跳跃屈曲；Abedi[3]将结构后屈曲释放的应 变能转化为结构失稳的初速度，求解结构的跳跃屈曲过程。本文采用球面弧长法计算结构 的屈曲前非线性反应，当结构随荷载增加达到临界点后，与一般后屈曲计算方法不同，不 用卸荷跟踪荷载位移反应全过程，而在临界点继续缓慢施加荷载，用动力时程方法模拟空 间刚架结构真实的跳跃失稳全过程。 空间刚架动态失稳全过程分析 69 空间梁单元的几何非线性切线刚度矩阵推导的主要方法有梁柱理论和有限元法。 Oran[4]采用梁柱理论，切线刚度矩阵达为超越函数的形式，并提出“节点方向矩阵”的 概念，对大转角进行修正。Meek[5]在线性有限元刚度矩阵中增加几何刚度矩阵，采用 Oran 的“节点方向矩阵”。文献[6]采用三维等参六面体单元推导曲梁的切线刚度矩阵，每段梁 采用 3个结点，但增加了结构整体分析计算量。 笔者基于连续介质力学对物质运动的描述方法，采用更新的拉格朗日公式，运用虚功 原理和等参有限元，推导空间梁单元大位移、大转角，小应变的切线刚度矩阵。由于节点 转角不满足交换律，不能进行矢量叠加，采用欧拉角对节点大转角进行计算。在材料非线 性分析中，根据Mises屈服准则和 Prandtl-Reuss塑性流动法则，推导材料弹塑性本构关系。 非线性时程分析的关键问题是数值积分的稳定性，本文采用 Newmark数值积分方法，在每 个时间步长内结合 Newton-Raphson平衡迭代，以保证计算精度和数值积分的稳定性。 2 空间梁单元非线性计算公式 按不同的参考位置，有限元法分 为全拉格朗日方法(Total Lagrangian) 和 修 正 拉 格 朗 日 方 法 (Updated Lagrangian)，修正拉格朗日方法的计 算效率更高。空间刚架的节点经历大 位移、大转角，但杆件的应变较小， 因此是大位移、大转角、小应变分析。 空间梁单元在笛卡尔坐标系中的运动 如图 1所示。 为求物体在离散时间点 0， tD ， 2 tD ，3 tD ⋯，上的平衡位置，假设已 知从时间 0 到时间 t 所包含的全部时 间步长下的静态变量和运动学变量， 通过某种递推方式计算 tt D+ 时刻物 体的平衡位置和全部状态变量。在拉格朗日增量分析方法中，用虚位移原理表示物体在时 间 tt D+ 的平衡, RVe ttttijtt Vtt ij tt DD DD D dt ++++ + =ò d (1) 其中 ijtt tD+ 为 Cauchy应力张量的笛卡尔分量， ijtt eD+ 为无穷小应变张量的笛卡尔分量，即 )( 2 1 )( 2 1 i tt j j tt i i tt j j tt i ijtt x u x u x u x u e DDDDD ¶ ¶d ¶ ¶d ¶ ¶ ¶ ¶ dd +++++ +º+= (2) 外力虚功 Rtt D+ 为 SufVufR ttsiS S i tt iV B i tttt tt tt tt dd DDDD dd D D D ++++ òò + + + += (3) 图 1 笛卡尔坐标系中梁的运动 70 工 程 力 学 B i tt fD+ 和 Si tt fD+ 分别为外部作用的体积力和表面力向量的分量， iud 是虚位移向量的第 i 个 分量。时间 t 的位移 i tu ，时间 tt D+ 的位移 utt D+ ， =i 1,2,3；因此得到 i tt ii tt i t ii t uxx uxx DD ++ += += 0 0 ( =i 1,2,3) (4) 从时间 t 到 tt D+ 的位移增量为 i t i tt i uuu -= +D ( =i 1,2,3) (5) 根据连续介质力学，式(1)可写成利用第二 Piola-Kirchhoff应力和 Green-Lagrange应变 张量表达的虚功表达式， RVS tttij tt tt ij tt t V DDD ed +++ =ò d (6) 式中第二 Piola-Kirchhoff 应力张量 ijtt t SD+ 和 Green-Lagrange 应变张量 ijtt t eD+ ，为从位形 t 至 位形 tt D+ ，但参照位形 t ，它们分别表示为下式 rj t ttsr tt si t tttt t ij tt t xxS ,, D D DD D t r r + + ++ + = (7) )( 2 1 ,,,, jktiktijtjitij tt t uuuu ++= + ded D (8) 应力分解为 ijtij t ij tt t SS += + tD (9) 式中 ijtt 为 Cauchy应力张量， ijt S 为第二 Piola-Kirchhoff应力张量增量参照位形 t ，考虑应 变增量 ijtt t eD+ ，有下式 ijtij tt t ee D =+ (10) ijtijtijt e he += (11) )( 2 1 ,, ijtjitijt uue += (12) jktiktijt uu ,,2 1 =h (13) 增量应力—应变本构关系采用下式 rstijrstijt CS e= (14) 方程(6)写为 òòò -=+ + V t ijtij ttt V t ijtij t V t ijtrstijrst ttt VeRVVC ddd dthdtede D (15) 上式为位移增量 iu 的非线性方程，不能直接求解，须进行线性化处理。 rstijrstijt ijtijt eCS e = =ded (16) 将方程(15)线性化为 空间刚架动态失稳全过程分析 71 òòò -=+ + V t ijtij ttt V t ijtij t V t ijtrstijrst ttt VeRVVeeC ddd dthdtd D (17) 方程(17)表示为位移增量的线性方程，是等参有限元的基本方程。 推导有限元控制方程的基本步骤为，选择插值函数，并在连续介质力学方程中利用这 些函数选择单元坐标和位移的插值，通过对每个节点位移依次引用虚位移原理，即得到有 限元控制方程。在等参有限元中，节点坐标和位移用相同的插值函数 å = = N i k iki xhx 1 00 ， å = = N i k i t ki t xhx 1 ， å = ++ = N i k i tt ki tt xhx 1 DD ( i =1,2,3) (18) å = = N i k iki uhu 1 00 ， å = = N i k i t ki t uhu 1 ( i =1,2,3) (19) k i t x 表示节点 k ， i 方向， t 时刻的坐标， kit u 表示节点 k ， i 方向， t 时刻的位移， kh 为 节点 k 的插值函数。N 是单元节点数。用式(18)、(19)计算式(17)中的积分，由单元分析可 以得到 FRuKK tt tt NL t tL t t -=+ +D)( (20) VBCBK tL t tV t T L t tL t t t dò= (21) VBBK tNL t tV tT NL t tNL t t t dtò= (22) VBF t V tT L t t t t t dtˆò= (23) 式中 Ltt K 为线性应变增量刚度矩阵; NLtt K 为非线性应变增量矩阵； Ftt 为与时间 t 的单元应 力相应的节点力向量； Rtt D+ 为外荷载向量； Ltt B 为线性应变位移变换矩阵； NLtt B 为非线性 应变位移变换矩阵； Ct 为增量应力应变材料特征矩阵。 tt 和 tˆt 为 t 时刻 Cauchy应力矩阵 和向量。构造适当的插值函数，计算 Ltt B 矩阵和 NLtt B 矩阵，通过积分求出 Ltt K 和 NLtt K ，得 到单元刚度矩阵。对节点大转角进行计算，将局部坐标系下单元刚度矩阵转化为整体坐标 系下单元刚度矩阵。 上述分析梁单元位移和应变—位移关系式的计算，这些运动学关系式的计算是直接 的，并且由于推导基于一般的虚功原理，从而这些关系式准确地描述了单元的几何非线性， 对物质运动描述能精确地表达大位移、大转角和大应变。实际结构采用普通碳素钢，钢材 的屈服强度较低。保证结构在任意荷载作用下为弹性是不可能的，也是不必要的。分析结 构的稳定性不仅要考虑几何非线性，还要考虑材料本构关系的物理非线性。本文采用双线 性等向硬化模型，具体推导见文献[7]。 3 算例分析 3.1 平面刚架 平面刚架由两根梁组成，将半跨结构划分为 5个梁单元，截面参数 EA=8.253×106N， EI=2.66×107 2mmN × 。分别假设结构为线弹性和理想弹塑性，屈服强度 ys =235N/mm2。 采用球面弧长法计算结构弹性和弹塑性的荷载位移全过程(图 2)。 72 工 程 力 学 该结构具有较高的几何非线性，结构弹性和弹塑性静力极限荷载分别为 P =155N， P =130N。弹性结构在超越极限点后承载力略微下降，然后进入稳定的后屈曲路径。弹塑 性结构在超越极限点后承载力大幅下降，释放大量的应变能后进入稳定的后屈曲路径。当 结构顶点位移刚过静力临界点，继续缓慢增加荷载，加荷速度 0.1N/s，进行非线性时程计 算。取节点集中质量，质量密度 r =0.0078kg/cm3，结构基本周期 =fT 0.233s，时间步长 tD =0.05s。采用瑞利阻尼，阻尼比z =0.02，阻尼系数a =0.13， b =0.00227。结构通过承 载力临界点后，外荷载继续缓慢增加，结构自然进入后屈曲跳跃失稳过程，由于阻尼作用， 结构动力反应迅速衰减，达到新的稳定平衡位置，完成跳跃屈曲全过程。结构跳跃后的弹 性和弹塑性变形由于结构失稳后释放应变能而都大于静力计算结果(图 2)。结构通过极限点 后的弹性及弹塑性荷载位移全过程见图 3。 图 2 双杆结构荷载位移全过程 (a) 弹性跳跃屈曲 (b) 弹塑性跳跃屈曲 图 3 跳跃屈曲荷载位移全过程 3.2 六角星型穹顶 六角星型穹顶边界铰接，截面参数 E=3.03×105N/cm2，G=1.096×105N/cm2，A=3.17cm2， 2I =2.377cm 4, 3I =0.295cm 4, J=0.918cm4。分别设结构为弹性和理想弹塑性，屈服强度 ys =500N/mm 2。采用球面弧长法计算结构弹性和弹塑性荷载位移全过程(图 4)。 结构顶点位移刚过静力临界点，继续缓慢增加荷载，加荷速度 0.056N/s。取节点集中 质量，质量密度 r =0.0051kg/cm3。结构基本周期 =fT 0.198s，时间步长 tD =0.025s。阻尼比 空间刚架动态失稳全过程分析 73 z =0.02，阻尼系数a =0.108， b =0.00368。结构通过极限点的荷载位移曲线见图 5。通过 极限点后结构的动位移持续增长，产生逃逸运动，结构整体翻转跳跃失稳。梁由压弯构件 转化为拉弯构件，由于构件的弹性卸荷作用，结构在新的振动平衡位置刚度矩阵恢复为正 定，结构围绕新的振动平衡位置继续衰减振动。 图 4 顶点荷载位移全过程 (a) 弹性跳跃屈曲 (b) 弹塑性跳跃屈曲 图 5 跳跃屈曲荷载位移曲线 3.3 60米跨度单层穹顶 Kaiwet型单层球面网壳，跨度 60m，矢高 6m，杆件采用f 194×10的 Q235钢管。边 界条件为固支。设材料为理想弹塑性，屈服强度 ys =235N/mm2。采用球面弧长法求出结构 弹性和弹塑性屈曲临界点，过临界点继续缓慢增加满跨均布荷载，加荷速度 9.8N/m2s。节 点集中质量,质量分布与静力临界荷载相同，时间步长 tD =0.02 s。阻尼比z =0.02，阻尼系 数a =0.286， b =0.00139。结构过极限点后弹性及弹塑性荷载位移全过程曲线见图 6。 74 工 程 力 学 (a) 弹性跳跃屈曲 (b) 弹塑性跳跃屈曲 图 6 顶点跳跃屈曲荷载位移全过程 结构失稳破坏模态对称，沿 6个径向肋出现局部凹陷变形，凹陷变形逐渐加深向周围 扩展，最终产生大面积的凹陷变形。考虑杆件的弹塑性，结构的大批杆件进入塑性，结构 没有稳定的后屈曲路径，结点位移持续增长，结构一直处于失稳运动状态，直至完全失稳 倒塌破坏，由正向壳变为反向凹陷壳，结构最终失稳破坏模态如图 7。 弹性失稳模态 弹塑性失稳模态 图 7 结构跳跃失稳模态 4 结论 通过静动力相结合的非线性时程方法研究结构的后屈曲问题，不用人为卸荷跟踪虚假 的荷载位移反应全过程，而在结构跨越临界点后继续缓慢加荷，使结构自然按实际情况进 入后屈曲路径，完成跳跃屈曲的动态过程。该方法不仅可以仿真空间刚架结构的失稳破坏 全过程，并从动力学角度和能量守恒原理进一步揭示结构的失稳机理，求出结构实际的承 空间刚架动态失稳全过程分析 75 载力和真实的安全储备。 参考文献： [1] E Risk. The application of Newton’s method to the problem of elastic stability[J]. J. Appl. Mech. 1972, 39: 1060-1066. [2] J L Meek, Qiang Xue. A study on the instability problem for 3D frame[J]. Computer Methods Appl. Mech. Engrg. 1998, 158: 235-254. [3] Karim Abedi, Gerard A. R. Parke. Dynamic propagation of local instability in single-layer braced domes[C]. In: Proceedings of Asia-Pacific Conference on Shell and Spatial Structures. T. Lan eds. Beijing. 1996. 556-563. [4] C Oran. Tangent stiffness in space frames[J]. ASCE, ST. Div. 1973, 99(ST6): 987-1001. [5] J L Meek, S Loganathan. Large displacement analysis of space frames structures[J]. Computer Methods Appl. Mech. Engrg. 1989, 72: 57-75. [6] K S Surana. Geometrically non-linear formulation for three dimensional curved beam element with large rotations[J]. International journal for numerical methods in engineering, 1989, 28: 43-73. [7] 王策. 单层球面网壳动力稳定性[D]. 哈尔滨建筑大学, 1997. DYNAMIC INSTABILITY ANALYSIS OF SPATIAL FRAME WANG Ce , LIU Xi-la (Tsinghua University, Beijing 100084) Abstract: In the present paper the non-linear static and dynamic post-buckling analysis of spatial frames is conducted. Based on continuum mechanics principles, the updated Lagrangian formulation is employed to develop a three-dimensional beam element which includes joint large displacements and large rotations. In material nonlinear analysis, the Mises yield criterion and Prandtl-Reuss flow rule are adopted to describe the elastic-plastic constitutive relation of the material. The Newmark integration combined with Newton-Raphon equilibrium iteration is used to solve the structural nonlinear vibration equation. The pre-buckling analysis employs the constant arc length method. The new method does not need to reduce the load in order to trace the complete load displacement curve. The applied load is increased slowly after the static critical point and the snap through process is simulated by the step by step integration method. Key words: spatial frame; post buckling; step by step integration; instability process simulation