第五章非线性有限元法

11
5.1.3 THE ARC-LENGTH METHOD 荷载增量法 [k (u )]{u} = λ ⋅ { f }
n (∆u ) = Ψ (u n ) = [k (u n )]{u n } − λ ⋅ { f }n = R n [ k ]T n
p
Softening: (∆λ ) n +1 需要一个负值
切线变刚度法
9
1. a new KT matrix has to be computed at each iteration; 2. if direct solution for Eq. is used the matrix needs to be factored at each iteration; 3. on some occasions the tangent matrix is symmetric at a solution state but unsymmetric otherwise (e.g. in some schemes for integrating large rotation parameters or non-associated plasticity). In these cases an unsymmetric solver is needed in general.
p=
σ 11 + σ 22 + σ 33
3
1 = tr (σ ) 3
4
5.1.1 Newton–Raphson method
(modified Newton methods) The most rapidly convergent process for solutions of problems.
Ψ (u ) = [ K (u )]u − R = 0
u
Ψu
n +1
5
Ψu
( )
n +1
∂Ψ (∆u )n = 0 ≈ Ψ (u ) + ∂u
n
n
解出⊿un
n

(∆u )
n T
n
= − [k ]
(
n −1 T
)
(∆u )
= − [k ]
(
n −1 T
)
Ψ (u )
n
Ψ (u n )
u n +1 = u n + ∆u n
10
5.1.2 荷载增量法的思想
引入标量荷载参量，通过调整荷载参量得到适当 的荷载。
在荷载 P 作用下得到有限元方程
[k (u )]{u} = { f }
定义参数 0≤λ≤1 λ0 = 0 λ1 = 0.1 ⋯ λ9 = 0.9 参数
[k (u )]{u} = λi ⋅ { f }
每一步荷载，迭代收敛后，进行下下一时步加载
n T n n T
原方程 [k (u )]{u} = { f }
n n n
[k ]a = -k(u )u + (λ ) (f E ) = R n
不平衡力
原方程 n (∆u )n = [k (u n )]{u n } − λ{ f }n = R n [ k ]T
14
整理归纳
[k ]a = R
n T
n T
2
5.1 Solution of non-linear algebraic equations
T u = [u1 u 2 u3 ...u n ] 方程 [k (u )]{u} = { f } 初值 u 0 = [u10 u 2 0 u3 0 ...u n 0 ]T
例如： 2 2 2 x1 + x2 + x3 − 1.0 = 0 1 2 2 初值 1 2 x1 + x2 − 4 x3 = 0 1 2 2 3 x1 − 4 x2 + x3 = 0 x1 K (u ) x2 x3 x1 1.0 与便利有关 2 x x − 4 x = 0 1 2 2 非常数 3 x − 4 x x 0 3 3 1
∂Ψ [k ] = ∂u
n
[k ]1 T
un
6
例如
k (u ){u} = R
3 u13 − u2 u1 = 2 3 3 3 u u − − 1 2 u2 = 2 3
3 2 u13 − u 2 R= 3 2 k (u ) = 3 3 − u1 − u 2 3
2
n
T
n
n
T
n
T
n 2
α1 + α 2 ∆λ + α 3 (∆λ ) = 0
n n 2
(∆λn )1, 2 =
确定根
− α2 ± α4 2α 3
n
α 4 = α 2 − 4α1α 3
2
若a4<0，虚根。求解需要用一个缩减因子
(∆λ ) = −αλ , α = ε /[ε ] 0.1 ≤ α ≤ 0.5
(
)
8
(4) 计算n+1时步的变量值
u1 u2
( n +1)
u1 = u2
(n)
∆u1 + ∆u2
(n)
(5) 判断收敛否
|| u
( n +1)
−u
(n)
||≤ ε
收敛，结束；否则，返回(1)继续计算 每次需要计算逆矩阵
∂u1
( n +1) 第n步计算结束，求n+1步的值 u1( n +1) u2 (n) 4 3 2 u − u 1 2 ( n ) −1 [ K 计算逆矩阵 3 n T ] (1) [k ]T = − u 2 4 u 3 1 2 3 4 (n) (n) 3 (n) u u u ( ) − ( ) ( (n) 1 1 2) (Ψ1 ) = -2 3 (2) 计算函数值 3 (n) (n) 4 (n) u u − u ( ) ( ) ( (Ψ ) ( n ) = − 1 2 2) -2 2 3 (n) (n) Ψ ( u ) ∆ u 1 计算增量 1 ( n ) − 1 (3) = − [k ]T (n) ∆ u 2 Ψ2 (u )
n
n
不平衡力 弧长法的基本公式
[k ]b = (f E )
定义变量
∆u n = a + b(∆λ ) n
s 2 = (u n + ∆u n ) T (u n + ∆u n ) s 2 = (u n + a + b∆λn ) T (u n + a + b∆λn )
弧长
= (u n + a ) T (u n + a ) + 2(u n + a ) T b∆λn + b T b(∆λn ) 2
(λ ) n +1 = (∆λ ) n +1 + (λ ) n
Softening range
u
土，混凝土，岩石等材 料，表现出软化现象。 混凝土开裂后的分析
12
n (∆u ) = [k (u n )]{u n } − λ{ f E }n = R n [ k ]T n
写成矢量形式 F(u , λ ) = [k (u )] u − λ ⋅ f E
3
表示为 Ψ (u ) = 0 → [ K (u )]u − R = 0 iterative method
Start K (u 0 ) = K 0
→ u 1 = [ K 0 ]−1 R to || u n +1 − u n ||< ε
→ u n +1 = [ K n ]−1 R
常用解法：牛顿法，拟牛顿法 讨论：牛顿迭代法 弧长法
n n n n n
n
0
把u和λ作为未知量，一阶Taylor展开
n
F(u
n +1
n
,λ
n +1
∂F ∂k ( u ) n = → [ k ] T ∂ u ∂ u
n n n n
∂F ∂F n ) ≈ F(u , λ ) + ∆u + ∆λn ∂u ∂λ n n n
σ 11 σ 12 σ 13 σ ij = σ σ 22 12 对称 σ 33 s11 s12 s13 sij = s s 22 12 对称 s33
σ = s + pI
1 0 0 I= 1 0 对称 1
切线刚度矩阵
∂Ψ1 4 3 u ∂u2 3 1 = ∂Ψ2 2 − u1 ∂u2 −u 4 3 u2 3
2 2
7
3 u13 − u 2 ∂Ψ1 Ψ1 = u1 - 2 = 0 3 ∂u 1 3 3 Ψ = − u1 − u 2 u - 2 = 0 [k ]T = ∂Ψ2 2 2 3
Ψ u n + (∆u ) = Ψ (u n +1 )
n
(
= u + ∆u
n
n
)
( )
n +1
∂Ψ (∆u )n = 0 ≈ Ψ (u ) + ∂u
n
n
∂Ψ ∂u
n
导数
K' (u n )
∂Ψ1 ∂Ψ1 ∂Ψ1 ∂u ∂u ..... ∂u 2 m 1 ∂Ψ2 ∂Ψ2 ∂Ψ2 ..... ∂u m = ∂u1 ∂u 2 ......... ∂Ψm ∂Ψm ..... ∂Ψm ∂u1 ∂u 2 ∂u m