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

第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 u M && (4.141)

式中，Ku u C 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 V

V =-=?δψεσδR (4.151)

式中, R 为载荷列阵；*δ为虚位移；*

ε为虚应变

用应变的增量形式d d =εB δ代入上式，消去*

δ项，可以得到非线性问题的一般平衡方程式为

()d 0T V

V =-=?ψδ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 V

V

V V =+??ψB σB σ (4.155)

不考虑初应变和初应力的影响，得

d d d ==σD εDB δ

并且

d d L =B B

这样可得

d d d d T L V

V =+?ψB σK δ (4.156)

这里

0d T L V

V ==+?K B DB K K (4.157)

式中0K 为通常的小位移的线性刚度矩阵。L K 矩阵则是由于大位移引起，它可以写成

()00d T T T L L L L L V

V =++?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 V

V σ=?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 y

y y ε??????=

+++???????? ??

?

?????+??+??+??=

222)()()(21z w z v z u z w z ε

z

w y w z v y v z u y u y w z v xy ????+????+????+??+??=

γ (4.161)

z

w

x w z v x v z u x u z u x w yz ????+

????+????+??+??=

γ y

w x w y v x v y u x u x v y u zx ????+????+????+??+??=

γ 对于微小位移情况，可以略去二次以上的偏导数项，得到小变形时的应变公式。在有限变形中，假设应变仍为小量。应变和位移之间的关系为：

0L =+εεε (4.162)

式中0ε为线性应变部分。对于非线性部分，可以写成：

0000011

22000T x T y x T z L y T T z y T T z z

x T T

y x ??????

???????

?==???????????????

?θθθθεθCθθθθθθθθ (4.163) 式中

[][

][

]T

x T

y T

z u v w x x x u v w y y y u v w z z

z

???=??????=??????=???θθθ (4.164) 式中C 为96?矩阵。

根据θ的定义，可以将θ表示成任意一点位移的函数，引入形函数N 后，可以得到

e =θGδ (4.165)

对于(4.163)式进行微分，得

11

d d d d 22

L =+=εCθC θC θ (4.166)

因此，

d e e L L ==εCGδB δ (4.167)

B 矩阵为

0L =+B B B (4.168)

这样得到

0d T L V

V +=?K K B DB (4.169)

另外，有

d d d d d

e T T T L V

V

V V σ==??K δB σG C σ (4.170)

利用矩阵C 和列阵θ的性质，得到

d d d x xy zx T

e yx y yz zx zy z στττστττσ????

==??????

I I I C σI I I θMG δI I I (4.171)

式中I 为三阶单位矩阵，M 是99?的六个应力分量组成的矩阵。因此几何刚度矩阵为

d T V

V σ=?k G MG (4.172)

故此，非线性三维单元的切线刚度矩阵为

0T L σ=++k k k k (4.173)

作为特例，可以直接写出三角形单元的上述有关表达式。由三角形单元的位移模式

e e

i j m u N N N v ??

??===??????

f N δI I

I δ (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 j

j i i

e v u v u v u }{}{=δ，式中的i i i c b a ,,等由

结点坐标确定，?为三角形单元的面积。

根据式(4.164)

T

x y u v u v x

x y

y ????????==???

???????

??θθθ 把式(4.174)代入上式得

00000010002000i j m i j

m e i j

m i j m b b b b b b c c c c c c ????

???=

???

????

?θδ (4.175) 由(4.165)式可以知道

0000001000200

i j

m i j m i j

m i j

m b b b b b b c c c c c c ????

???=

????????

G (4.176) 根据定义，由式(4.163)确定的平面问题的C 矩阵为

0000

0T

x T y T T y x u v

x x u u y y u v u v y

y

x

x ??????????????

????==???????????

?????????????????

θC θθθ (4.177) 这样可以得到C 的显式为

10

02i i j j m m

i i j j m m

i i j j m m i i j j m m i i j j m m i i j j m m

i i j j m m

i 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 由线性问题给出，即

000010

002i j m i j m i i

j

j

m

m b b b c c c c b c b c b ?????

=?????

?

B (4.180) 至此，0B 、L B 和G 都是常数矩阵，只与单元结点坐标和结点位移有关。因此，线性

刚度矩阵为：

000T t ?=k B DB (4.181)

式中t 为单元厚度。初始刚度矩阵为：

00()T T T L L L L L t ?=++k B DB B DB B DB (4.182)

几何刚度矩阵为：

T t σ?=k G MG (4.183)

式中

00000000x xy x xy xy y xy y στσττστσ??

???

?

=????????

M (4.184) 因此，对于几何非线性问题，平面三角形单元的切线刚度矩阵可以由式(4.173)求出。

10.2 材料非线性问题的有限元法

材料非线性是指材料的本构方程是非线性的。一般主要分为两类: 一类是非线性弹性问题，如橡胶、塑料、岩石等，在加载时应力应变关系等性质呈现非线性的物理现象，卸载时可逆。另一类是指材料的弹塑性问题。材料超过屈服极限后呈现出非线性。

在机械结构分析中，常见的本构方程主要有线弹性和非线性弹性模型，其特点是应力仅应变的函数，加卸载规律相同。公式如下：

kl ijkl ij C εσ= (4.185)

其中，对于线弹性材料为常数，对于非线性弹性材料，ijkl C 是kl ε的函数。

此外，还有超弹性模型、次弹性模型、弹塑性模型等。对于弹塑性模型，可以认为弹塑材料发生塑性变形时，其总应变可以分解为两部分：

p ij e

ij ij εεε+= (4.186)

即总应变为弹性应变和塑性应变之和。加载时遵循一定规律，如Prandtl-Reuss 方程，而卸

载时为弹性。对于应力足够大时的金属、土壤、岩石等材料，都有此类特征。

在非线性动力有限元分析时，具有黏性特性的模型十分重要。对于黏弹性模型，一般包括松驰（指突加应变作用下应力逐渐减少）和蠕变（指突加应力作用下应变会逐渐增加）。典型的松驰模型是如下Maxwell 模型：

φσσ

ε

+=E && (4.187) 蠕变模型如Voigt-Kelvin 模型为：

ε

φεσ&+=E (4.188) 而弹/粘塑性模型是指材料的塑性变形与时间有关，本构方程中出现非齐次的时间微分。

下面以弹/粘塑性材料的变形分析为例，说明非线性动力有限元的基本步骤。即弹/粘塑性材料的变形过程中要考虑时间效应，材料开始屈服后，塑性流动、应力和应变均与时间有关。

10.2.1 弹/粘塑性问题的基本表达式

假设总应变ε分离成弹性应变e ε和粘塑性应变vp ε，即

e vp =+εεε&&& (4.189)

其中) (&表示对时间的求导。总应力率取决于弹性应变率，有

e =σDε&& (4.190)

式中，D 是弹性矩阵。粘塑性性质为

0),(0=-F F vp εσ (4.191)

式中，0F 是单向屈服应力，是硬化参数的函数。

设粘塑性应变率仅取决于当前的应力，如下式所示：

}

{)(}{σφγε??>

<=Q

F vp & (4.192) 式中，),(vp Q Q εσ=是塑性势，γ是流动参数，><)(x φ项对于0>x 是正单调增量函数，即

?

?

?≤>>=<0 00

)()(x x x x 当当φφ (4.193) 在这里只讨论F Q ≡的情况。函数φ的两种常用形式为

1)()(

-=-F F F M e

F φ (4.194)

或

N

F F F F )(

)(0

0-=φ (4.195) 式中，N M ,为常数。

10.2.2 粘塑性应变增量和应力增量

对于式(4.192)所表示的应变定律，定义时间间隔n n n t t t -=?+1出现的应变增量为

,vp n ?ε，有

, ,,1[(1) ]vp n n vp n vp n t ?ε?ΘεΘε+=-+&& (4.196)

其中Θ为常数， 01Θ≤≤。

上式中的,1vp n +ε&可用如下近似公式表达

,1,vp n vp n n n ?+=+εεH σ&& (4.197)

式中, vp n n

???

=

????εH σ&. n H 取决于应力水平，

可导出显式公式，具体内容参见有关文献[欧文，塑性有限元]。

对于应力增量n ?σ，有

,()n en vp n ????==-σD εD εε (4.198)

用位移增量表示时，有

,()n n n vp n n t ???=-σD B δε& (4.199)

10.2.3 弹/粘塑性平衡方程

在任何瞬时n t 的平衡方程应满足

d 0T

n n n V +=?B σR (4.200)

在时间增量中，平衡方程由增量形式给出，即

d 0T n

n

n

V ??+=?B σR

(4.201)

在时步n t ?中出现的位移增量能为

1

n T n n V ??-=δK (4.202)

,d T

n n n vp n n V t V ???=+?B D εR & (4.203)

这里，,T n K 是切线刚度矩阵，使用如下形式，即

,d T

T n n n n V =?K B D B (4.204)

把位移增量n ?δ代回式(4.199)，可以得到应力增量n ?σ. 再有

1n n n ?+=+σσσ (4.205) 1n n n ?+=+δδδ (4.206)

且有

1,vp n n n n ???-=-εB δD σ (4.207)

,1,,vp n vp n vp n ?ε?ε?ε+=+ (4.208)

应力增量的计算是基于增量平衡方程(4.201)的线性化形式，累积所有这样的应力增量得到总应力1n +σ是不正确的，并不会真实地满足平衡方程(4.200)。为此，可以采用计算残余平衡力的方法进行迭代求解，即

111d 0T n n n n V +++=+≠?ψB σR (4.209)

对于几何非线性问题，1n +B 由位移1n +δ计算得出。然后，残余平衡力1n +ψ迭加到下一个时间步的载荷增量上。

在弹/粘塑性分析中，时间步长的选择十分重要，限于篇幅，这里不再加以讨论。

10.3 材料非线性问题的动力有限元法

在材料非线性问题的动力有限元分析中，首先考虑材料的非线性本构关系，此外还要考虑是否包含大位移、大应变等因素，在实际工程分析中，需要对具体不同问题具体处理，目前没有统一的方法。在这里只限于考虑小应变条件下的材料非线性问题，且对有限元动平衡方程进行时间上的数值积分，以实现动力分析。

根据虚功原理，可以导出结构材料在n t 时刻的动平衡方程。平衡方程与材料本身的性质无关。对于每一个结点i 应满足相应的平衡方程：

0n n n n n

i Bi Ii Di Ti -++-=P f f f f (4.211)

其中n

i P 为内阻力：

d d n nT n nT n

n n i i i vevp i i V

V

V d V ==??P B σB D B (4.212)

体力一致力为(n

b 为作用的体力矢量)

d n

nT Bi i n V

V =?f N b (4.213)

式中，n

i N 为n t 时刻的形函数。 惯性力为(n

ρ为密度)

d n nT n n n Ii i i i V

V ρ=?f N N d && (4.214)

阻尼力为(n

C 为阻尼系数阵)：

d n nT n n n Di i i i V

V =?f N C N d & (4.215)

边界面力的一致力为(n

s 为面力矢量)：

d n

nT n Ti i Γ

Γ=?f N s (4.216)

式中, Γ为受面力作用的边界与单元边界重合的部分。

利用Gauss-Legendre 乘积方法，依据单元形函数，可以对上述各式进行数值计算。 设材料的粘弹性本构方程vp ve e εεεε++=中各部分如下： 1) 弹性应变

σεD

e 1

=

(4.217) 其中D 为弹性矩阵，设为常数矩阵。

2) 粘弹性应变增量

粘弹性应变率可以表示为

1

2

2)1(

ηεσεE A E n ve n n ve -=& (4.218) 即

1

2

1221ηεσηεE A E E E n ve n n

ve -+=

& (4.219)

在1+n t 时刻为

1

2

1112211

ηεσηεE A E E E n ve n n ve +++-+=

& (4.220)

且有

??????+=?+=++n

n n n n n ε

εεσ

σσ11 (4.221) 采用如下插分格式

])1[(1

1+++-?=?n ve n ve n n ve t εθεθε&& (4.222)

其中，θ为插分系数，如0，0.5，1等。则粘弹性应变增量与应力增量之间的关系为

n n n

n ve n

n

n

ve t E t A E t E t σηηθθεηθ

ε????

? ???+?+

?+?=

?112

21

2

11& (4.223)

3)粘塑性应变增量

当满足一定的屈服准则时，材料产生粘塑性变形。采用类似的插分方程

])1[(1

1+++-?=?n vp n vp n n vp t εθεθε&& (4.224)

其中1

+n vp ε&用Taylor 级数来近似

n

n n vp n vp H σεε?+=+&&1 (4.225)

其中n

H 由屈服准则和粘塑性流动准则确定

n

T T n

vp

n

aa df

d a H }{)(>Φ<+??>Φ<=??=σγσε& (4.226)

采用Ducker-Prager 准则：

k J J '='+2

1α (4.227) 式中k ',α由材料粘聚力和内摩擦角定；

1J 为应力第一不变量；2J '为应力偏量的第二不变量，当取关联的流动法则时有

n n J J

f Q a )()(21σ

σασσ?'?+??=??=??=

(4.228) 这样，由式(4.224)和式(4.225)得到

n n

n n n vp n vp t H t σθεε?+?=?& (4.229)

4)总应变增量

由应变分解可知

n

vp

n ve n e n εεεε?+?+?=? (4.230) 其中总应变增量可以由结点位移得到

n n n d B ?=?ε (4.231)

式中, n B 为位移应变矩阵；n

d ?为单元结点位移增量。

由式(4.223)、式(4.229)和式 (4.230)可得应力增量为：

?

??????-??+-?=?n n vp n n ve n n n vevp n t t t E D εεηθεσ&&12/11 (4.232)

n

vevp

D 是粘弹粘塑性模量 1

11221)/1(--?

?????+?+?+=n n n n n vevp t H t E t A E D D θηηθθ (4.233)

5) 动态有限元刚度矩阵

在任一n t 时刻，可以证明，该时刻的切线刚度矩阵为

?=V

n n

vevp nT n T dV B D B K (4.234)

3. 显式时间积分法

任一n t 时刻的动平衡方程可以写成如下矩阵形式：

n n n n f P d C d M =++&&& (4.235)

式中M 为总质量矩阵，C 为总阻尼矩阵，n

P 为内力的总矢量，n

f 为作用的体力、面力的

一致结点力矢量。

为简便起见，用中心差分公式对上述动平衡方程进行离散化。加速度、速度分别为：

}2{)

(11

2

-+-?≈

n n n n d d d t d && (4.236) }{2111-+-?≈n n n d d t

d & (4.237) 通过差分处理，在1+n t 时刻的位移可用n n t t ,1-时刻的位移显式地表过出来。

根据上述思路，可以对这类非线性材料进行动力有限元计算，得到相应的动应力、动应

变响应等结果。

10.4 应用举例

10.4.1 粘弹粘塑性动力有限元分析举例

编制粘弹粘塑性动力有限元程序，对受到振动压实物料的位移及应力应变响应特性进行了如下研究。

假设物料是均匀的，形成粘弹粘塑性材料结构，如图4-8所示为该结构某断面的平面四边形单元网格划分。在上边界的某一点受到周期性集中载荷的作用，载荷具有不对称性。两侧边界无约束。与基础层结合处的底边结点均受约束。假设材料满足等向强化的Drucker-Prager屈服条件。材料参数如表4-3。

表4-3 材料参数取值

Poision 比弹性

模量

弹性

模量

粘滞

系数

粘滞

系数

粘聚力内摩

擦角

强化

系数

密度流动

性系

数

μ1E

(N/mm2)

2

E

(N/mm2)

1

η

(sec-1)

2

η

(sec-1)

c

(N/mm2)

φκρ

(kg/mm3)

γ

0.3 6.8804 2.164 7.5e-5 1.9e-7 7.07 62.73 1.0e-3 3.0e-6 1.0e-2

X

Y

图4-8 受集中载荷的物料层的网格划分

设载荷)(t

P为

??

??

?≤≥=0)(;

00)(;2sin )(0t P t P vt b t P 当当π (4-174) 其中，60Hz v =，5

0 1.010N b =-?。取时间步长为dt =3.333X10-6sec 。在不同的时刻，材

料断面的应力场分布不同。某结点垂直方向(y 方向)上的应力时间历程如图4-9（a ）所示。由于作用力的不对称性，应力应变关系呈现近似的不对称滞回特性，如图4-9(b)所示。

(a) (b)

图4-9 材料断面某结点垂直方向的应力历程及其应力应变关系

习题