有限元-第6章

（6-8）
式中， [ B0 ] 是一般线性分析由小位移应变得到的单元几何矩阵， [ BL ] 由大 位移非线性应变引起的。 将式（6-7）引入式（6-6）得单元的平衡方程
∫∫∫[ B ]
Ve
T
{σ }dv − { f e } = 0
（6-9）
在分析非线性问题时，多采用增量列式方法，以下着重介绍这一方法。 将式（6-9）所示的平衡方程写成微分的形式
第六章 结构几何非线性分析的有限元法
§ 6-1 引言 上述各章介绍的都是关于线性弹性体的有限元分析。线性弹性体的特点 是材料的应力应变关系满足虎克定律，位移是微小的。位移与载荷呈线性关 系。 但是许多实际问题，位移与载荷不呈线性关系。一类属于材料非线性， 材料的应力与应变的关系是非线性的。另一类属于几何非线性，几何非线性 问题指的是大位移问题，但是有很多大位移问题中，结构内部的应变是微小 的，材料的应力应变关系还是线性的。对于几何非线性问题，必须依照弹性 体变形后的位移建立平衡方程，由于变形后的位移是未知的，这给处理非线 性问题带来一定的复杂性。但是结构的非线性分析可以按一系列的线性段来 求解，这样可将适用于线性问题的分析方法和结果经过修改来解决非线性问 题。本章仅对工程结构中遇到的结构大位移的几何非线性问题的有限元分析 进行讨论。 几何非线性问题的平衡方程必须建立在结构变形之后的状态上，为了描 述结构的变形需要设立一定的参考系统，一种做法是让单元的局部坐标系跟 随结构一起发生变化，由此便产生了带有流动坐标的迭代法，这种做法对于 杆系的非线性有限元分析特别显示其优越性。另一种做法是让单元的局部坐 标系始终固定在结构发生变形之前的位置，以结构变形前的原始位形作为基 本的参考位形，这种分析方法称为 Lagrange 方法。对于同一个物理问题来说， 其本身的规律并不因为选择不同的参考系统而发生变化。
§ 6-2 带有流动坐标的迭代法 带有流动坐标的迭代法是指结构在发生大位移的过程中，使各单元的局 部坐标系跟随结构一起运动，由此来描述结构的非线性。这一方法对于杆系 结构的大位移分析特别显示其优越性。尤其在杆件发生比较大的转动时，采 用这一方法比采用以下要介绍的总体 Lagrange 方法更为适宜。这是因为通过 局部坐标系的流动可以方便地描述单元的刚体运动，从而较容易地得到变形 后的单元在变形后结构中的位置。 无论是小位移问题还是大位移问题，结构在承受载荷发生变形后，必定 满足平衡方程。结构坐标系中的单元节点力是由结构坐标系中的单元刚度矩 阵与该单元的节点位移的乘积得到。结构坐标系中的单元刚度矩阵不仅取决 于单元本身的属性，还与该单元所处的方位有关。在线性小位移分析中，由 于节点位移引起单元方位的变化十分微小因而可以不计，在计算结构的节点 位移和合力时可仍然用处于变形以前方位的单元刚度矩阵。而在结构大位移 的非线性问题的分析中，单元方位的变化不能忽略，必须利用处于变形后方 位的单元刚度矩阵来计算节点合力。 首先，讨论流动坐标系中单元节点位移和节点力的计算方法。图 6-1（a）
§ 6-3 总体的 Lagrange 列式法 如果始终以结构未变形的原始位形作为参考位形进行有限元列式称为总 体的 Lagrange 列式法。采用这种列式方法，单元局部坐标系始终固定在结构 变形前的位置，单元局部坐标系与结构总体坐标系之间的转换关系可始终保 持不变，此时，按线性理论推导的单元刚度矩阵已不再适用，而需要推导出 在大位移情况下的单元刚度矩阵。可以推断，此时的刚度矩阵应是节点位移 的函数。 在非线性问题中，刚度矩阵的含义是载荷—位移曲线的斜率，这个斜率 在不同的位移时是不同的，故称为切线刚度矩阵。一般来说，在求解非线性 问题时，可以把原属于非线性的载荷—位移关系看成是由一系列的线性段组 成，于是，就希望求得单元节点力增量与单元节点位移增量之间的关系，这 种关系可通过单元切线刚度矩阵表达，由单元切线刚度矩阵组装成结构的切 线刚度矩阵，切线刚度矩阵是节点位移的函数。 下面介绍总体的 Lagrange 列式法的基本理论和切线刚度矩阵的一般表达 式。 若结构处于平衡时，单元也是平衡的，由最小势能原理 δΠ e = 0 ，有
由以上分析可以看出，进行结构的大位移分析，可以将按线性分析得到 的节点位移作为结构位移的第一次近似值，然后根据上述节点位移对单元刚 度矩阵进行修改，从而反映变位后的单元在变位后的结构中所发挥的作用。 带有流动坐标的迭代法分析结构大位移问题的基本思路是：根据已求得 的节点位移值修改单元的坐标变换矩阵 [λ ] ，从而达到修正结构坐标系中单元 刚度矩阵的目的。根据修正后的各单元刚度矩阵乃至刚度方程可以计算出节 点合力。按上述结构位移第一次近似值算出的节点合力与节点所受到的外载 荷并不相等，即此时并未达到节点的平衡，这是因为按线性分析所得到的节 点位移并不是结构的真正平衡位置，结构当然也就无法保持平衡。为了找到 结构真正的平衡位置，可以将结构的节点载荷与上述节点合力之差、各个节 点不平衡力作为一组新的外载荷施加于上述已发生变形的结构上，求得节点 位移的增量，将这个节点位移增量与原求得的节点位移相加便得到节点位移 的第一次修正值。根据节点位移的修正值再重新修改各单元的坐标变换矩阵， 并重新计算新的节点合力和节点不平衡力，继而再将这些节点不平衡力施加 于结构。重复进行上述过程一般可以使节点不平衡力减少到可以被忽略的水 平，此时的节点位移所对应的变式结构在发生大位移后的真正平衡位置，按 照这个平衡位置可以计算处结构在大位移情况的内力。 综上所述，在采用带流动坐标的迭代法时，一个典型的迭代过程应包括 下列步骤首先为外载荷 {P} 作用下的结构假定一组节点位移 {∆} 。 (1) 根据结构整体坐标下的节点位移列阵 {∆} 确定单元的位置，建立单元的 局部坐标系； (2) 计算各单元在局部坐标系中的节点位移列阵 {δ ′ e } ，形成在局部坐标下 个单元的刚度矩阵 [k ′ e ] ，并计算相应的单元节点力 { f ′ e } ；
Ve
T
{σ }dv + ∫∫∫[ B ]T d {σ }dv = d{ f e }
Ve
（6-10）
对于线弹性材料，应力应变关系为
{σ } = [ D]{ε }
（6-11）
应力增量与应变增量之间的关系是
d{σ } = [ D]d{ε } = [ D][ B ]d{δ e }
（ c）
且从式（6-8）有
中给出在结构坐标系 xy 的一个未变形的平面刚架的梁单元 ij ，由节点 i 和 j 的 坐标值算出 x0 、 y 0 、 L0 和 α 0 。单元变形后，整体坐标下的单元节点位移为
{δ e } = [u i vi θ i uj vj
θ j ]T
（6-1）
单元变形后的位置和形状如图 6-1（b）所示。 x ′y ′ 为单元的流动坐标。也即局 部坐标，它的原点位于变形后的杆端 i ， x ′ 轴沿变形后的杆端节点连线方向。 由图 6-1（a）(b)可得
δΠ e = ∫∫∫δ {ε }T {σ }dv − δ {δ e }T { f e } = 0
Ve
（6-6）
增量形式的应变—位移关系可表示为
d{ε } = [ B ]d{δ e }
（6-7a）
上式 d{δ e } 表示单元节点位移 {δ e } 的微分，根据变分与微分运算在形式上的相 似性，有
δ {ε } = [ B ]δ {δ e }
(3) 将 [k ′ e ] 和 { f ′ e } 变换到整体坐标系得到 [k e ] 和 { f e } ； (4) 集 合 各 单 元 刚 度 矩 阵 ， 形 成 结 构 刚 度 矩 阵 [ K ] = ∑ [k e ] 和 节 点 合 力
{F } = ∑ { f e } 。 [ K ] 为在当前变形位置的结构刚度矩阵。
[k L ] = ∫∫∫[ B0 ]T [ D][ B L ]dv + ∫∫∫[ B L ]T [ D][ B0 ]dv + ∫∫∫[ BL ]T [ D][ BL ]dv
∫∫∫ d ([ B ]
Ve
T
{σ }) dv − d{ f e } = 0
（ a）
由于几何矩阵 [ B ] 和应力 {σ } 都是单元节点位移的函数，因此有
d ([ B ]T {σ }) = d [ B ]T {σ } + [ B ]T d{σ }
（b ）
将式（b）引入式（a） ，则有
∫∫∫ d[ B ]
N = {δ }T {δ }
则位移收敛条件的一种形式是
N n − N n −1 ≤e Nn
式中， e 是精度要求，可根据工程的要求和问题的性质而定，一般可取 10 −6 ～
10 −2 。作为另一种更加严格的收敛判断，位移收敛条件可以要求每一个节点自
由度的位移满足
∆δ i δ i ≤ e
有时，位移收敛条件和力收敛条件同时运用。 结构所承受的载荷 {P} 可以一次计入，也可以分成 n 个载荷增量。当给定 的载荷水平达到收敛后，可以增加到下一新的载荷水平，在进行迭代找出一 个新的平衡位置。这样把迭代计算和增量计算结合起来，为了达到最后的载 荷水平，可以取许多增量步而每步用较少的迭代次数，也可以取较少的增量 步而每步用较多的迭代次数，以提高计算的精度。
(5) 计算不平衡力 {∆F } = {P} − {F } 。 (6) 求解结构方程 [ K ]{∆δ } = {∆F } ，得到节点位移增量 {∆δ } 。将 {∆δ } 加到前 次迭代中积累起来的节点位移 {δ } 中去，就给出节点位移的新近似值。 (7) 检查收敛性，如果不满足收敛性判断则返回到步骤（1） 。 以上的过程可概括为如下的迭代公式
通过坐标变换矩阵 [λ ] ，得到在结构总体坐标系中单元节点力与节点位移之间 的关系
{ f e } = [k e ]{δ e }
而
[k e ] = [λ ]T [k ′ e ][λ ]
由于单元的方位角 α 与单元节点位移有关，使得单元刚度矩阵 [k e ] 是单元 节点位移列阵的函数，即