第10章(非线性有限元)分解

公式号、图号等第十章 非线性动力有限元法当机械结构受到较大的外载荷，或受到持续时间较短的冲击载荷作用时，结构会产生过大的变形, 以至于必须考虑结构几何大变形对结构整体刚度及固有频率的影响，即所谓的几何非线性影响。

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

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

P 为外部激励。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4.6.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 可能由于变形而发生的变化，即在这里假设了载荷不因变形而改变其大小和方向，否则是非保守力作用下的大变形问题，在此不做讨论。

4.6.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)zwx w z v x v z u x u z u x w yz ∂∂∂∂+∂∂∂∂+∂∂∂∂+∂∂+∂∂=γ yw x 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⨯矩阵。

根据θ的定义，可以将θ表示成任意一点位移的函数，引入形函数N 后，可以得到e =θGδ (4.165)对于(4.163)式进行微分，得11d d d d 22L =+=εCθC θC θ (4.166)因此，d e e L L ==εCGδB δ (4.167)B 矩阵为0L =+B B B (4.168)这样得到0d T L VV +=⎰K K B DB (4.169)另外，有d d d d de T T T L VVV V σ==⎰⎰K δB σG C σ (4.170)利用矩阵C 和列阵θ的性质，得到d d d x xy zx Te yx y yz zx zy z στττστττσ⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦I I I C σI I I θMG δI I I (4.171)式中I 为三阶单位矩阵，M 是99⨯的六个应力分量组成的矩阵。

因此几何刚度矩阵为d T VV σ=⎰k G MG (4.172)故此，非线性三维单元的切线刚度矩阵为0T L σ=++k k k k (4.173)作为特例，可以直接写出三角形单元的上述有关表达式。

由三角形单元的位移模式e ei j m u N N N v ⎧⎫⎡⎤===⎨⎬⎣⎦⎩⎭f N δI II δ (4.174)其中，∆++=2/)(y c x b a N i i i i ，∆++=2/)(y c x b a N j j j j ，∆++=2/)(y c x b a N m m m m ，T m m jj i ie v u v u v u }{}{=δ，式中的i i i c b a ,,等由结点坐标确定，∆为三角形单元的面积。

根据式(4.164)Tx y u v u v xx yy ⎧⎫⎧⎫∂∂∂∂==⎨⎬⎨⎬∂∂∂∂⎩⎭⎩⎭θθθ 把式(4.174)代入上式得00000010002000i j m i jm e i jm i j m b b b b b b c c c c c c ∆⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦θδ (4.175) 由(4.165)式可以知道00000010002000i jm i j m i jm i jm b b b b b b c c c c c c ∆⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦G (4.176) 根据定义，由式(4.163)确定的平面问题的C 矩阵为00000Tx T y T T y x uv x x uu y y u v u v yyxx ⎡⎤∂∂⎢⎥∂∂⎢⎥⎡⎤⎢⎥∂∂⎢⎥==⎢⎥⎢⎥∂∂⎢⎥⎢⎥⎣⎦⎢⎥∂∂∂∂⎢⎥∂∂∂∂⎢⎥⎣⎦θC θθθ (4.177) 这样可以得到C 的显式为1002i i j j m m i i j j m mi i j j m m i i j j m m i i j j m m i i j j m mi i j j m mi i j j m m b u b u b u b v b v b v c u c u c u c v c v c v c u c u c u c v c v c v b u b u b u b v b v b v ∆⎡⎤++++⎢⎥=++++⎢⎥⎢⎥++++++++⎣⎦C (4.178)故L =B CG (4.179)而0B 由线性问题给出，即000010002i j m i j m i ijjmm b b b c c c c b c b c b ∆⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦B (4.180) 至此，0B 、L B 和G 都是常数矩阵，只与单元结点坐标和结点位移有关。