弹塑性力学的非线性有限元

有较高的二次收敛速率， 切向刚度矩阵m+1[K](i-1)在每个迭代步中都要重新计算和分解 对于理想弹塑性材料或应变软化材料，切向刚度矩阵病态，s困难
P u
改进的Newton-Raphson法
使用第n个(n<m+1)加载步时计算所得的切向刚度矩阵n[K]替代切向刚度 矩阵m+1[K](i-1)。
准Newton法
（1）是N-R法和改进的N-R法之间的一个折衷方法。 （2）使用低秩矩阵去更新刚度矩阵m+1[K](i-1)的逆矩阵。Broyden–Fletcher-
Goldfarb-Shanno(BFGS)方法就是其中的一种。 （3）准Newton法的收敛速率介于线性收敛和二次收敛之间。 （4）可适用于应变强化、应变软化或理想塑性等分析。可以考虑卸载。
p u
改进N-R法的特点 （1）比 N-R法减少了刚度矩阵的计算和分解。 （2）是线性收敛，通常比N-R法收敛得慢，如在分析应变软化材料时， 收敛将会特别地慢。 （3）刚度矩阵可能变成奇异矩阵或病态矩阵的问题仍然存在。 （4）如果出现卸载，应力状态从塑性状态卸载到弹性状态，这个算法 可能得不到一个收敛结果，除非一旦卸载出现，刚度矩阵重新计算。
本构方程
（1）增量本构关系，是无穷小应力增量与应变增量的关系。
（2）加载步中的荷载增量是有限值，应力和应变增量也为有限值。
（3）必须对增量本构关系在加载步内积分，确定有限应变增量ij 与有限应力增量ij的关系
m1
m1
ij
dij
C ep ijkl
d
kl
m
m
其中
C ep ijkl
切线模量为
C ep ijkl
力边界S上的面力是 m1 X i mX i X i
位移边界Su上的位移是
u m1 i
m
ui
ui
求m+1增量加载步物体内引起的位移增量u i和应力增量ij
用增量表示的求解方程
平衡方程
mij,i +ij,i+mXj+Xj=0
几何方程
1
1
mij +ij = 2 (mui,j+muj,i) + 2 (ui,j+uj,i)
高斯点的应变增量 m+1{}(i)= [B] m+1{u}(i) {}=m+1{}(i) m{}
分三种情况讨论由高斯点的{}计算实际产生的{} 假设产生的是弹性变形，应力增量是 {e}=[D]e{}
（1）第m加载步末，高斯点的应力状态位于屈服面内，f (m{}，m)<0， 在第m+1步加载后，f (m{}+{e}，m)<0，应力状态仍位于屈服面 内， 这说明该高斯点实际就处于弹性状态。那么，实际应力增量就是 {}={e}
其中
P
{u}(i) m1 {P}(i1) (m1) {F} 0
{u} m1{u}(i1)
{u}(i)=m+1{u} m+1{u}(i-1)
m1{P}(i1) m1 {P( m1{u}(i1) )} [B]T m1 {}i1 dQ Q
P
{u} m1{u}(i1)
[B]T
Q
{}
{}
{} {u}
dQ
[B]T [D]ep
Q
[B]dQ m1 k i1
m 1{u}(i 1)
N-R法的迭代格式可写为
m+1[K](i-1) {u}(i) = m+1{F} m+1{P}i1
m+1{u}(i)=m+1{u}(i-1) +{u}i
m+1{u}(0)=m{u}， m+1[K](0)=m[K]， m+1{P}(0)=m{P} N-R算法的特点
（2）对于每一步荷载增量，求解物体内应力和应变的增量， （3）“积分”（累计）得到最终的应力和应变。
具体求解步骤
第m增量加载步末 体积力是mXi， 力边界S上的面力是，位移边界Su上的位移是 物体内的位移mui、应变mij和应力mij均已求出
在第m+1增量加载步，给定增量的荷载和位移
体积力是
m+1Xi= mXi +Xi
关键的问题是如何选定[K] 的值 增量法的主要优缺点：
（1）优点：它适用于各种类型的非线性状态， （2）缺点：计算效率不高。当要求足够精度的解答，必须采用足够小的 荷载增量，这导致较大的计算量。
P u
增量迭代法
思路： 将加载步内的荷载增量适当加大，在每个加载步内进行若干次迭代
目的： 为了提高计算效率
d kl
d
p kl
d
1 h
f ij
d ij
dij
Cijkl dkl
d f kl
d
f ij
Cijkl dkl
h
f ij
Cijkl
f kl
d ij
Hale Waihona Puke Cep ijkld
kl
C ep ijkl
Cijkl
Cijmn
f mn
f pq
C pqkl
h
f ij
Cijkl
f kl
与弹性刚度张量具有相同的对称性
S
Q
控制方程改写成
[B]T {}dQ = {F}+ m{F} m{P} Q
不平衡节点力
(m{}) = m{F} m{P(m{})} 最终的增量求解方程
[K]{u}== {F}+(m{u})
增量法
每一个荷载增量步，解一次线性方程组，其实质就是以分段线性来近似 非线性。
将求解的位移增量、应变增量和应力增量叠加到上一步的位移、应变和 应力中去
Cijkl
Cijmn
f mn
f pq
C pqkl
h
f ij
C ijkl
f kl
ij
C
ep ijkl
kl
其中代表第m+1加载步中间的某一状态。 Euler近似方法：是使用增量加载前切线模量
C C ep m ep
ijkl
ijkl
弹塑性切线刚度的推导
dij
C
ijkl
d
e kl
Cijkl
ij已精确满足m加载步末的所有方程和边界条件则可以从中消去它们得到只有增量的一组方程在求解过程中它们不一定能精确满足方程和边界条件将它们保留下来可减小累积误差增量有限元格式对于m1加载步末假想它发生虚位移uijijm1加载步末的节点力平衡控制方程改写成每一个荷载增量步解一次线性方程组其实质就是以分段线性来近似非线性
本构积分
必要性 第m+1加载步第i次迭代计算已完成，进行下一次迭代(i+1)，需要得
到第i次迭代结束时应力m+1{}(i)，应变m+1{}(i)及内变量m+1{}(i)。 原则
根据当前加载步内（从m加载步末到m+1加载步第i次迭代）求出的位 移增量{u}，确定应变增量{}、应力增量{}和内变量增量{}，然 后与m加载步末的对应量相加 方法： （1）确定应力所处的状态，即是处于塑性加载、弹性、还是卸载 （2）使用弹塑性本构关系通过数值积分计算相应的应力增量{}
方法 有限元控制方程用位移表示为 (m+1{u}) = m+1{P(m+1{u}) m+1{F} = 0 m+1{u}= m{u}+{u}， m+1{}= m{}+{} 问题归结为求m+1{u}，使得不平衡力为零，即(m+1{u})=0。
取上一步末的位移作为下一步的初值，m+1{u}(0) = m{u}， 不平衡力一般不为零，(m+1{u}(0))0， 通过迭代方法逐步消除不平衡力使之为零。 在有限元分析中广泛应用三种Newton迭代法。
m+1[K](0) = m[K] = n[K] 改进的N-R法的迭代方案可表示为
n[K]{u}(i) = m+1{F} m+1{P}i1
m+1{u}(i)=m+1{u}(i-1) +{u}i
m+1{u}(0)=m{u}， m+1{P}(0)=m{P}
(i=1,2,…)
若在每一个步的每一次迭代中，始终采用初始弹性矩阵[K]0进行计算， 称为初应力法。
{}={e} 如果f (m{}+{e}，m)>0，则高斯点继续处于塑性加载状态，中取r=0 来计算应力增量。
m1
m r
m1
{} [D]ep{d}
m
m
[ D]e {d}
[D]ep {d}
m r
m1
{} r[D]e{}
[D]ep d
m r
（3）第m加载步末，高斯点处于塑性状态，f (m{}，m)=0。 在第m+1步加载后，f (m{}+{e}，m)<0，处于卸载状态或中性变
载状态，其应力增量是
塑性力学问题的有限元方法
增量方法 增量迭代方法
塑性力学问题的提法
增量理论： 应力全量和应变全量不单值对应，取决于历史，且是非线性
边值问题的特点： 物体内的应力应变分布取决于：
（1）边界上荷载和位移的最终值（2）达到最终值的路径 求解方法：
（1）通常将荷载分成若干个增量，从初始状态开始，沿着给定的 加载路径，将荷载增量逐步施加
，
有限元控制方程
，
[B]T m1 {}dQ [N ]T m1 {X }d [N ]T m1 {X }dQ
Q
S
Q
m+1{P}= m+1{F}
m+1 加载步末的节点力平衡
m1{P} [B]T m1{}dQ Q
m+1{F} [N ]T m1 {X }d [N ]T m1 {X }dQ
（2）几何关系的矩阵形式 {}=[B]{u}e
{}={x，y，z，xy，yz，zx}T （3）本构方程式的矩阵形式为
｛｝= [D]ep{}
{}={x，y，z，x，y，yz，zx}T
f
f
x
f y
，
f f f
z xy yz
T
f
zx
D ep
De
De f f T De h f T De ，f
边界条件
ni( mij +ij) = m X j X j mui +ui = m ui ui
两点说明：
在力边界S上 在位移边界Su上
（1）除本构方程外上述其他方程及边界条件都是线性的。 （2）若mui、mij和mij已精确满足m加载步末的所有方程和边界条件
则可以从中消去它们，得到只有增量的一组方程
Newton-Raphson迭代方法
假设在（m+1）步，第（i-1）次迭代近似值m+1{u}(i-1)已经得到。 使用Taylor级数将Ψ(m+1{u})在m+1{u}(i-1)处展开，并忽略高次项
( m1{u}(i1) )
( m1{u} m1 {u}(i1) ) 0
{u} m1{u}(i1)
在求解过程中，它们不一定能精确满足方程和边界条件
将它们保留下来，可减小累积误差
增量有限元格式
对于m+1加载步末，假想它发生虚位移(ui)，虚应变是 (ij ))
m1 ij ij dQ＝ m1 X i ui d m1 X i ui dQ
Q
S
Q
（1）将单元内的位移增量表示成节点位移增量的插值形式 {u}=[N]{u}e
这个方法为常规弹塑性非线性联立方程组的求解提供了一个安全有效 的方法。这是目前使用得最好的算法。
位移准则
收敛准则
{u}(i)
2 D
m1{u}(i) m {u} 2
力准则
m1{F} m1 {P}(i)
2 F
m1{F} m {P} 2
能量准则
{u}(i)T ( m1{F} m1 {P}(i) ) E {u}(i)T ( m1{F} m {P})
（2）第m加载步末，高斯点处于弹性状态，f (m{}，m)<0， 但是f(m{}+{e}，m)>0，即在第m+1步中间，高斯点进入塑性状态 因此，存在一个比例因子r，使得 f(m{}+r{e}，m)=0 应变增量因而被分成两部分： 第一部分是纯弹性变形 r{} 第二部分将产生塑性变形。 (1r){}