第10章 非线性动力有限元法 (1)10.1 几何非线性问题的有限元法 (2)10.1.1 几何非线性问题的牛顿迭代法 ........................................................................... 2 10.1.2 典型单元的切线刚度矩阵 ................................................................................. 4 10.2 材料非线性问题的有限元法 (8)10.2.1 弹/粘塑性问题的基本表达式 .............................................................................. 8 10.2.2 粘塑性应变增量和应力增量 ............................................................................... 9 10.2.3 弹/粘塑性平衡方程 ............................................................................................ 10 10.3 材料非线性问题的动力有限元法 ................................................................................ 11 10.4 应用举例 (14)10.4.1 粘弹粘塑性动力有限元分析举例 ................................................................... 14 习题.. (15)第10章 非线性动力有限元法当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形, 以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。

另外, 对于多数非线性动力学问题，还需要考虑材料非线性、接触非线性等方面的影响。

非线性动力学分析求解的基本方程有如下形式0=-+P I uM (4.141) 式中，Ku uC I += 为粘性效应项，考虑阻尼、粘塑、粘弹等效应。

P 为外部激励。

对于考虑各种非线性效应的动力学问题求解，需要对动力学方程进行直接时间积分。

即非线性动力有限元分析具有如下特点：（1）问题分析过程需要考虑时间积分效应，不必做模态分析，不必提取固有频率；（2）采用直接积分方法求解非线性动力学方程，需要对时间作积分计算，因此计算量远远大于线性模态动力学方法；（3）非线性动力学分析中可以施加不同类型的载荷，包括结点力、非零位移、单元载荷；（4）在每个时间步上，进行质量、阻尼、及刚度的集成，采用完整矩阵，不涉及质量矩阵的近似；（5）可以同时考虑几何、材料和接触等多种非线性效应。

非线性动力有限元分析程序常采用隐式Hilber-Hughes-Taylor 法进行时间积分运算。

这种方法适于模拟非线性结构的动态问题，对于冲击、地震等激发的结构动态响应以及一些由于塑性或粘性阻尼造成的能量耗散，隐式算法特别有效。

隐式积分方法需要对刚度矩阵求逆计算，并通过多次迭代求解增量步平衡方程。

隐式Hilber-Hughes-Taylor 时间积分算法为无条件稳定，对时间步长没有特别的限制。

采用子空间法也可以对动力学平衡方程作时间积分运算。

子空间法是提取模态分析得到的各阶特征模态，并采用与线性模态动力学分析方法相近的分析方式进行求解。

对于带有微小非线性效应的问题，如材料小范围进行入屈服、结点转角不大的情况，子空间法效率比进接积分法要高。

此外，非线性动力有限元分析还可以采用显式动态算法，如中心差分法。

显式时间积分算法为有条件稳定，其临界稳定时间步长限制了时间步长的大小，与有限元模型最小单元尺寸、材料应力波速等有关。

显式时间积分法适于模拟高速冲击、接触等问题。

上述方法的选择需要综合考虑计算量、分析问题的规模、单元限制等多方面因素，需要丰富的有限元模拟的理论、经验和实践知识。

以下以几何非线性问题和材料非线性问题为例介绍非线性有限元法，其中粘弹粘塑性非线性材料问题的分析是典型的非线性动力有限元的求解思想。

10.1 几何非线性问题的有限元法几何非线性问题一般是指物体经历大的刚体位移和转动，但固连于物体坐标系中的应变分量仍假设为小量, 即大位移小应变情况。

10.1.1 几何非线性问题的牛顿迭代法由数值分析技术可知，求解非线性方程组的数值方法的常规方法是Newton-Raphson 法，即牛顿迭代法，这是一种近似线性化迭代求解方法。

对于非线性方程0)(=x ψ，具有一阶导数，在n x 点作一阶泰勒级数展开，它在n x 点的线性近似为d ()()()()d n n n x x x x xψψψ=+- (4.142) 因此，非线性方程0)(=x ψ在n x 附近似为线性方程：d ()()()0d n n n x x x xψψ+-= (4.143) 当d () 0d n xψ≠时，由上式求得n 步的修正项 1d ()/()d n n n X x xψ∆ψ+=- (4.144) Newton-Raphson 方法的迭代公式为11++∆+=n n n X x X (4.145)在几何非线性有限元法中，结构的刚度矩阵与其几何位置有关，平衡方程由变形后的位形描述，因此，结构的刚度矩阵是几何变形的函数。

设变形为δ, 结构的平衡方程式()0-=K δδR (4.146)为一个非线性方程组。

记非线性方程()0K =-=ψδδR (4.147)用Newton-Raphson 方法求()0=ψδ的根时，迭代公式分别为11n n n ++=+δδδ (4.148)其中,1n ∆+δ满足下式1()T n n n ∆+=-K δR K δδ (4.149)式中, T n K 称为切线刚度矩阵，表达式为d ()()d T n n =ψδK δ(4.150) 在每一个迭代步中，通过求解切线刚度矩阵T n K ，进而用1n ∆+δ进行迭代求解，称为Newton-Raphson 方法，又称切线刚度法。

牛顿法的收敛性是好的。

但是某些非线性问题中，使用牛顿法迭代时，若T n K 出现奇异或病态，则对T K 的求逆出现困难。

关于这一点也可以采用其它修正办法，如引入阻尼因子。

对于已经建立的有限元方程，设ψ表示内为和外力矢量的总和，有***d 0T T T VV =-=⎰δψεσδR (4.151)式中, R 为载荷列阵；*δ为虚位移；*ε为虚应变用应变的增量形式d d =εB δ代入上式，消去*δ项，可以得到非线性问题的一般平衡方程式为()d 0T VV =-=⎰ψδB σR (4.152)该式不论位移或应变的大小与否均成立。

在有限变形中，应变和位移之间的关系是非线性的，即B 矩阵是δ的非线性函数。

但是，近似地可将进行如下分解：0L =+B B B (4.153)式中, 0B 为线性应变分析的部分； L B 为由非线性变形引起的，与δ有关。

假定应力应变关系为线弹性，于是有00()=-+σD εεσ (4.154)式中 ][D 为材料的弹性矩阵； }{0ε为初应变列阵；}{0σ为初应力列阵对于式(4.152)的非线性平衡方程式，可用Newton-Raphson 方法进行迭代求解。

对该式微分，有d d d d d T T VVV V =+⎰⎰ψB σB σ (4.155)不考虑初应变和初应力的影响，得d d d ==σD εDB δ并且d d L =B B这样可得d d d d T L VV =+⎰ψB σK δ (4.156)这里0d T L VV ==+⎰K B DB K K (4.157)式中0K 为通常的小位移的线性刚度矩阵。

L K 矩阵则是由于大位移引起，它可以写成()00d T T T L L L L L VV =++⎰K B DB B DB B DB (4.158)式(4.156)又可记成：()0d d d L T σ=++=ψK K K δK δ (4.159)式中d d T L VV σ=⎰K B σ (4.160)式中，σK 是关于应力水平的对称矩阵，称之为初应力矩阵或几何刚度矩阵。

因此，用Newton-Raphson 方法迭代求解几何非线性问题的步骤为： (1) 用线弹性解作为1δ，即一次近似； (2) 通过定义1()δB 求出1σ，求出1ψ； (3) 确定切线刚度矩阵1T K ；(4)211/T ∆=-δψK ， 212∆=+δδδ；(5) 重复上述迭代步骤，直至n ψ足够小。

在这里，没有考虑载荷R 可能由于变形而发生的变化，即在这里假设了载荷不因变形而改变其大小和方向，否则是非保守力作用下的大变形问题，在此不做讨论。

10.1.2 典型单元的切线刚度矩阵求解具体的几何非线性问题时，必须计算单元的切线刚度矩阵。

对于一般空间问题，无论位移和应变大小，都可以利用应变的基本定义写出位移和应变的关系式。

用变形前的坐标),,(z y x 做为自变量，可以用位移w v u ,,定义如下大变形问题的应变分量表达式⎥⎦⎤⎢⎣⎡∂∂+∂∂+∂∂+∂∂=222)()()(21x w x v x u x u x ε2221()()()2y v u v w y yy y ε⎡⎤∂∂∂∂=+++⎢⎥∂∂∂∂⎣⎦⎥⎦⎤⎢⎣⎡∂∂+∂∂+∂∂+∂∂=222)()()(21z w z v z u z w z εzwy w z v y v z u y u y w z v xy ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ (4.161) zw x w z v x v z u x u z u x w yz ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ ywx w y v x v y u x u x v y u zx ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ 对于微小位移情况，可以略去二次以上的偏导数项，得到小变形时的应变公式。

在有限变形中，假设应变仍为小量。

应变和位移之间的关系为：0L =+εεε (4.162)式中0ε为线性应变部分。

对于非线性部分，可以写成：000001122000T x T y x T z L y T T z y T T z zx T Ty x ⎡⎤⎢⎥⎢⎥⎧⎫⎢⎥⎪⎪⎢⎥==⎨⎬⎢⎥⎪⎪⎢⎥⎩⎭⎢⎥⎢⎥⎣⎦θθθθεθCθθθθθθθθ (4.163) 式中[][][]Tx Ty Tz u v w x x x u v w y y y u v w z zz∂∂∂=∂∂∂∂∂∂=∂∂∂∂∂∂=∂∂∂θθθ (4.164) 式中C 为96⨯矩阵。