増分型非線形解法の導入

前項で示した弾塑性体の境界値問題を弱形式定式化することにより得られた基

礎式(1.568)は, 単純引張など変位, 応力などの分布が明らかになっている場合を

除き解析的に解を求めることはできない. そのため後述するような方法により有 限要素離散化を行い, 近似解を求める.

一般に弱形式であらわされた偏微分方程式の境界値問題を有限要素離散化する と, 以下のようになる.

find u_h ∈ W_h such that

w_h·(Q(u_h)−F) = 0 ∀w_h ∈ Wh (1.570) ただし, W_h は許容関数全体の空間 W を有限要素離散化したもので, W_h ⊂ W で ある. Q,F はそれぞれ内力ベクトル,外力ベクトルと呼ばれる.

式(1.570)が成立するための必要十分条件は,

Q(u_n) =F (1.571)

であり,式(1.571)は平衡方程式と呼ばれている.

外力ベクトル F に対して {F^k} を

0 =F⁰ <F¹ <F² <· · ·<Fⁿ=F (1.572) のように定める. 増分型非線形解法では,

Q(u^k) =F^k (1.573)

なる u^k を u^k−1 をもとに順次求めていく. ただし F⁰ = 0 に対応する解u⁰ は明 らかに 0である.

a. forward-Euler 法

u^k−1が得られているとき,F^k−F^k−1が十分小さければu^kは以下の方程式を解 くことにより近似できる.

K^k(u^k−u^k−1 =F^k−F^k−1 (1.574) K^k = ∂Q

∂u

u=u^k−1 (1.575)

Q K

^k

u

^k

u

^k−1

F

^k−1

F

^k

図 1.7: 近似解u^k

K は接線剛性マトリックスと呼ばれている.

これをもとに, forward-Euler法では, まずstep 数 n を設定し,外力 F から外力 ベクトルの増分 ΔF = ¹

nF を求め, 下式にしたがって変位増分 Δu^k を求める.

K^kΔu^k =ΔF (1.576)

これにより u^k は以下のように求められる.

u^k=u^k−1+Δu^k (1.577)

= k

i=1

Δuⁱ (1.578)

このアルゴリズムの概念図を図1.8に示す.

この図からわかるように, step 数を十分大きく設定しないと,誤差が蓄積してい く. ただし弾塑性解析では, 1step 内で降伏あるいは除荷する積分点の数はなるべ く少ない方が好ましいため step 数を大きく設定しなければならず, そのため従来 から弾塑性解析においては forward-Euler 法はよく用いられている.

u Q

u

³

u

²

u

¹

Δu

³

Δu

²

Δu

¹

ΔF

ΔF ΔF

K

¹

K

²

K

³

図 1.8: forward-Euler法の概念図 b. Newton-Raphson 法

これに対し Newton-Raphson 法では反復計算を行うことにより誤差の蓄積をな くす. 具体的にはu^k−1 が得られている時,次のような手順で u^k を求める(図 1.9 参照).まず,

u

^k₀

= u

^k−1

u F

K

^k₂

F

^k

F

^k−1

Q

^k₁

Q

^k₂

Q

^k

K

^k₁

Δu

^k₁

Δu

^k₂

u

^k

u

^k₂

u

^k₁

F

^k

− Q

^k₀

F

^k

− Q

^k₁

F

_k

− Q

^k₂

図 1.9: Newton-Raphson法 概念図

u^k₀ =u^k−1 (1.579)

Q^k₀ =Q(u^k₀) (1.580)

K^k₁ = ∂Q

∂u

u=u^k₀ (1.581)

とする. K は接線剛性マトリックスである. ここで,

K^k₁Δu^k₁ =F^k−Q^k₀ (1.582) u^k₁ =u^k₀ + Δu^k₁ (1.583) の様に u^k₁ を定めると, u^k₀ よりも u^k₁ の方が u^k に対するよい近似となる.以下, u^k_i−1 が得られた時

Q^k_i−1 =Q(u^k_i−1) (1.584)

K^k_i = ∂Q

∂u

u=u^k_i−1 (1.585)

K^k_i Δu^k_i =F^k−Q^k_i−1 (1.586) u^k_i =u^k_i−1+ Δu^k_i (1.587) の様に u^k_i を求めることを繰り返す. u^k_i−1 よりもu^k_i の方がu^k に対するよい近似 となっており, Q に極値がなくF^k に対する解が存在する,F^k−Q^k₀ が大きす ぎない, などの条件が整っていれば u^k_i は u^k に,F^k−Q^k_i は 0にそれぞれ収束 する. コンピュータで計算を行なう場合は, F^k−Q^k_i= 0 となることはないの で, 適当な収束判定で計算を打ち切る. なお, F^k は,解析の step 数をn とすると,

F^k = k

nF (1.588)

とすることが多い.

前述のように,弾塑性解析では, 降伏,除荷の判定を行なうため, step数を大き設 定する必要がある. そのため反復を行なわなくてもほぼ正確に近い値が得られて いると考えらるため, ここでは反復を省略して解析を進めることにする.

具体的には u^k−1 が得られている時,下式に従い u^k を求める.

K^kΔu^k =F^k−Q(u^k−1) (1.589)

u^k=u^k−1+Δu^k (1.590)

このアルゴリズムの概念図を図1.10に示す. なお比較のため前述の forward-Euler 法による結果をあわせて示す. 図中で がついたものが forward-Euler 法のもので ある. forward-Euler 法とほぼ同様のアルゴリズムであるが, 内力と外力の不平衡 力を考慮しているため, forward-Euler 法よりも誤差の蓄積が軽減できる.

u

³

u

²

u

¹

= u

¹

Δu

³

Δu

²

Δu

¹

ΔF

ΔF ΔF

K

¹

K

²

K

³

F

¹

F

²

F

³

Q(u

¹

)

Q(u

³

)

Q(u

²

)

Δu

¹

Δu

²

Δu

³

u

³

u

²

F

²

− Q(u

¹

)

F

³

− Q(u

²

) K

²

K

³

図 1.10: 反復しない forward-Euler 法の概念図 c. 解析例による手法の比較

第 ?? 章において説明するモデルを例にとり, forward-Euler 法と反復しない

newton-Raphson法により解析を行い両者の結果を比較する. ただし解析モデルの

肉厚は均一で,ソリッド要素を用いる. Prager 則による移動硬化則を用い,n 乗則 でカーブフィットを行った. 材料物性は Young率 19592 kgf/mm², Poisson 比 0.3, 降伏応力 25 kgf/mm², n = 0.2, C₂ = 0.001 とする.

図1.11 は, forward-Euler 法と反復しない Newton-Raphson 法により, 幾何学 的非線形性を考慮し, エルボに閉→開→閉の強制変位入力を与える 1 サイクル

を 1200 step に分割し, 2 サイクル解析を行った結果の比較である. 反復しない

Newton-Raphson 法によるload-stroke 線図はほぼ一定の所をまわっているのに対 して, forward-Euler 法による load-stroke 線図はドリフトしている. 図1.12 は同 じ条件でエルボを開く側から解析を行った結果である. 同様にドリフトしている ことがわかる. これに対して図1.13に示すように幾何学的非線形性を考慮しない 場合にはドリフトはみられなかった.

-6 -4 -2 0 2 4 6 8 10 12 14

-40 -30 -20 -10 0 10 20 30 40

load (ton)

stroke (mm) forward-Euler Newton-Raphson without iteration

図 1.11: load-stroke線図 (幾何学的非線形)

-4 -2 0 2 4 6 8 10 12 14

-40 -30 -20 -10 0 10 20 30 40

load (ton)

stroke (mm) forward-Euler Newton-Raphson without iteration

図 1.12: load-stroke 線図(幾何学的非線形 逆サイクル)

-6 -5 -4 -3 -2 -1 0 1 2 3 4

-40 -30 -20 -10 0 10 20 30 40

load (ton)

stroke (mm)

forward-Euler Newton-Raphson without iterarion

図 1.13: load-stroke 線図 (幾何学的線形)

forward-Euler 法も, 反復しない Newton-Raphson 法も, step 数を増加させれば 正確に近づくことが期待できる. そこで1サイクルを10倍の 12000 stepに分割し て解析を行った. forward-Euler法の結果を図1.14に,反復しないNewton-Raphson 法の結果を図1.15 に示す. 図1.16に 12000 stepで解析した際の forward-Euler 法 と, 反復しない Newton-Raphson 法の結果の比較を示す. 図1.15からわかるよう に, 反復しない Newton-Raphson 法を用いた場合, 1 サイクルあたり 1200 step の

結果も 12000 step の結果もほぼ一致している. これより 1 サイクルあたり 1200

step の分割で十分であることがわかる. 図1.14からわかるように, forward-Euler 法を用いた場合でも step 数を増加させればドリフトは軽減される. しかしながら 図1.16からわかるように, いぜんとして forward-Euler 法の結果にはドリフトが 含まれており, 1 サイクルあたり 12000 stepでは不十分である.

以上の結果から, 本研究では, 反復しない Newton-Raphson 法を用いることに する.

-6 -4 -2 0 2 4 6 8 10 12 14

-40 -30 -20 -10 0 10 20 30 40

load (ton)

stroke (mm) 1200 step

12000 step

図 1.14: load-stroke線図 (forward-Euler 法)

-6 -5 -4 -3 -2 -1 0 1 2 3 4 5

-40 -30 -20 -10 0 10 20 30 40

load (ton)

stroke (mm) 1200 step

12000 step

図 1.15: load-stroke線図 (反復しない Newton-Raphson 法)

-6 -4 -2 0 2 4 6

-40 -30 -20 -10 0 10 20 30 40

load (ton)

stroke (mm) forward-Euler Newton-Raphson without iteration

図 1.16: load-stroke 線図(幾何学的非線形 12000step)

1.1.13 updated Lagrange 法による内力ベクトルと剛性マトリッ

ドキュメント内 Report98.dvi (ページ 65-77)