Transcript 超弾性体解析で現れる剛性行列の性質とその解法に関して

超弾性体解析で現れる剛性行列の
性質とその解法に関して
鷲尾 巧 久田 俊明
東京大学
日本応用数理学会「行列・固有値の解法とその応用」研究部会
第６回研究会
本発表の概要
• Aが有限要素離散化においてどのようなプロセスを経
て構成されるか
• Aの正定値性に関する議論
• 全体行列がどのようなプロセスを経て構成されるか
• 反復解法における前処理の構成法について
• 収束性について
ほとんど話題は、不思議に思うこと、自分ではわかっていないことです!!
A B 
A 

B  C
T
ゴムと金属の違い
ヤング率 (MPa)
鉄鋼 2x106
ゴム 1～3
圧縮率 (Pa－1)
鋳鉄
ゴム
0.6x10－11
53.7x10－11
弾性変形の限界の伸び
金属 1%以下
ゴム 1000%
ゴム構造
自然状態
伸ばす
ヘルムホルツの自由エネルギー
結晶構造
エントロピーSの減少
ポテンシャルEの増加
F  E  TS
統計力学的ミクロモデルにおける非圧縮性の関与
（参照：ゴムの弾性（裳華房） 久保亮五 著）
非圧縮的伸長
高さ方向のエントロピー
伸ばす
Sz  a 2  c
１
１
1
1


 2 2 
S    a    3
 

2 

F    TS    aT  2   3
 

dF
1

 2aT    2  応力と歪みの関係式
d
 

水平方向のエントロピー
1
Sx  S y  a  c

応力

１
1
2
伸び
3
4
連続体の変形を表すテンソル
Lagrange座標
（変形前の位置を表す）
Euler座標
（現配置での位置を表す）
xX  X  uX
X
dx  FdX
dX
変形勾配テンソル
 x1 x1 x1 



X

X

X
2
3
 1
 x2 x2 x2 
F

 X1 X 2 X 3 
 x x x 
3
3
 3
 X1 X 2 X 3 
右極分解 本質的な変形を取り出す
F  RU
R: 回転行列, U:正定値対称行列
F F  UR RU  U
T
T
2
歪みテンソル
CF F
T
右Cauchy-Green変形テンソル
T
T

1
1  u u u u 
E  C  I   

2
2  X X X X 
Green-Lagrangeひずみ
圧縮
E11
引張り
E22
e2
せん断 E12
e1
E11 E12 
E

E
E
 21 22 
弾性ポテンシャルWとその変分
W  WE
歪みEを与えるために要する単位体積あたりのエネルギー
第１変分式
S11 S12 S13 
W
W 
Eij  SijEij S  S21 S22 S23
Eij
S31 S32 S33 
第２変分式
2

W
W 2
2
 W  Eij
Ekl 
 Eij
Eij Ekl
Eij
 Eij DijklEkl  Fki SijFkj
初期変位項
初期応力項
第２Piola-Kirchhoff応力テンソル
応力テンソルの物理的意味
F1dfn
F1
N
dS
参照配置
df n
n
ds
現配置
dfn  TT  nds Cauchy応力テンソル
F1dfn  ST  NdS 第2Piola - Kirchhoff応力テンソル
剛性行列の正定値性
uT Au Eij DijklEkl  Fik SijFjk
2W
Dijkl 
Eij Ekl
E Bu, F Z1u
初期変位項 : Wが凹関数ならば正定値
初期応力項 : Sが正定値ならば正定値
初期応力項の正定値性は不確定：
Sが正定値 ⇔ 主応力方向がすべって引張り
主不変量
I  TrC  1  2  3
1
II  TrC2  Tr C2  12  23  31
2
III  123
1 , 2 , 3 : Cの固有値

 
Cayley- Hamiltonの定理
C3  I C2  II C  III I  0
等方的連続体の弾性ポテンシャルは主不変量の関数として表わされる。
W  W I , II , III 
不変量から定まる剛性行列の符号
第１不変量
初期応力項：正定値
初期変位項：ゼロ
I

S

I
I
SijI  2
 2ij , Dijkl
 2 ij  0
Cij
Ckl
第２不変量
初期応力項：正定値
初期変位項：ゼロ３負３
II

S
II
II
SijII  2
 2Iij  Cij  , Dijkl
 2 ij  4ij kl  ik jl 
Cij
Ckl
第３不変量 初期応力項：正定値 初期変位項：ゼロ３負３
III

S

III
III
SijIII  2
 III C1 ij , Dijkl
 2 ij  4III C1 ij C1 kl  C1 ik C1
Cij
Ckl
 
        
第２、３不変量の剛性行列（初期応力項＋初期変位項）の符号判定は困難
jl
ゴム状物質の弾性ポテンシャル例
Mooney-Rivlin体
W  c1 I  3  c2 II  3

III  1

非圧縮性制約条件付き
あるいは
1
W  c1 I  3  c2 II  3   III 2
2
ペナルティ方的なポテンシャル
など
第２不変量の凹関数性を保証できないため、Newton-Raphson解法の収束性を
一般に保証できない。
ポテンシャルWのLandscape(想像）
境界および制約条件を満たす解の集合
初期状態
心筋の弾性ポテンシャルの例
C
~
Mazhari and McCulloch,
W  e 1 , Q  b E , b  0 Usyk,
Journal of Elasticity 61,2000
2
Q
ij
e3
e2
2
ij
ij
~
e3
e1
自然基底
~
e2
sheet
~
e1
fiber
異方性に従った基底
  2Q

Q ~

Q

Q
~ W
Sij  ~  W ~ , Dijkl  W  ~ ~  ~ ~ 
 Eij Ekl Eij Ekl 
Eij
Eij


初期応力項： 歪みに依存
初期変位項： 正定値
Q
 2Q
~
~  2bij Eij , ~ ~  2bij ik jl
Eij
Eij Ekl
異方性連続体モデルの不安定性




異方性連続体モデルは、初期変位項が正定値となるように設計される。
しかし、安定性議論においては初期応力項の影響が無視されているようである。
異方性モデルは、変形過程の初期において、解けなくなるケースが多々生じる。
安定な第１不変量を少し加えることにより、不安定性を回避できる。
非線形
（２次項あり）
線形
変位
u
変位勾配テンソル
u
Z
 FI
X
凹関数
弾性ポテンシャル
歪みテンソル


1
E  Z  ZT  ZT Z
2
W E
 歪みテンソルが変位の２次関数となっていることが問題。
 しかし、回転不変性などを考慮しつつ弾性ポテンシャルを変位勾配Z
の凹関数として直接的に定義することは困難。
体積弾性の取扱について
1 2
W W0  
2
  J , J  det F 体積剛性
 
未定乗数の導入
W0


0
Sij 

 Sij  
Eij
Eij
Eij
第１変分式 （制約条件式も付加）
W  SijEij     
第２変分式 （制約条件式も付加）
 2W  Eij DijklEkl  Sij 2 Eij  Eij


  
Ekl  
Eij
Ekl
2W0
2
Dijkl 

Eij Ekl
Eij Ekl
不定値か正定値悪条件か？
未定乗数導入時


 W  Eij DijklEkl  Sij Eij  Eij
  
Ekl  
Eij
Ekl
2
2
A0 BT 
A 

B  C 
剛性行列
直接的変分

不定値

2 W0 1 2 2
 W  Eij
Ekl  Sij 2 Eij
Eij Ekl
2
剛性行列
A  A0  BT C1B
正定値 しかし悪条件
固有値分布 (1)
1

A B  A 0 A 0  A BT 
A 
   

1  
B  C  B S  0  S   0 S 
T
S  C  BA B
1
T
#positive eigenvalues = dim (A)
#negative eigenvalues =dim(C)
a
b
c
0
d
固有値分布 (2)
A BT  u  ru 

    
B C    r 
 A BT  u  ru 

    
 B C     r 
固有値は右半面に含まれる
T

 u 
A
B
* *
u λ 
 
 B C   


 u*Au  *C  2 1 Im(*BT p)
| u*BT p |
ブロックタイプ前処理行列
A BT  A 0 A1 0  A BT 
A 
   

1  
B  C  B S  0  S   0 S 
PLU
QA 0  QA1 0  QA BT 





B QB   0  QB  0 QB 
1
A
QA  A QB  S or QB  C  BQ B
良いQB の構築は簡単ではない。
T
ブロックタイプ前処理の効果
T

M
(I  M)N 
1
1
T U PLU AT U  
T
 N(I  M) N(2I  M)N 
M  QA1/2AQA1/2, N  QB1/2BQA1/2
M  I as QA  A
 T
1/2
1 T 1/2
1 T
NN

Q
BQ
B
Q

I
as
Q

BQ
B
A
B
B
AB

Q1/2
QA1/2BT 
A
TU  

1/2
QB 
 0
正定値ブロックタイプ前処理
A BT  A 0 A1 0  A BT 
A 
   

1  
B
S
B  C     0  S   0 S 
QA 0  QA1 0  QA BT 
~
PLU  




B QB   0 QB  0 I 
T


M
(I

M)N
~ 1
1
T U PLU AT U  
T
N(I  M)  N(2I  M)N 
I 0 
0  I


MINRESなどの対称行列向き（しかし非正定値）反復法が適用可能
ＩＬＵ分解とSchur complement
M  L  DD1 (D  U)
d
u
L, D およびUを足し合わせ Q=L+D+Uとおく。
Set Q  A
for i  1,, n
for k  1,,i 1 only for i, k  P
for j  k  1,, n only for k , j  P
if i, j  P then
l
Q
d
u
l
Q－l u/d on P
Q
on PC
qij  qij  qik  qkk 1  qkj
ILU分解の安定性
対角行列Dの符号が完全分解時の符号に
従うことが望ましい
不完全なSchur complementを計算
より実践的な前処理（ILU前処理）
•変数分離のみで純代数的に構成できる。
•フィルインの設定に変数分離情報を利用する。
~T 
0  DA1 0  DA  UA B
DA  LA
U
P  ~




1 
B
D

L
0
D

U
0

D
L
B
B 

B
B
B 


k
j
(bi,bk,bj)=(block(i),block(k),block(j))
k
i
(bi,bk,bj)
Allowed fill-ins
(2,1,1) , (1,1,2)
Level 0 or 1
(2,1,2)
All fill-ins
(1,1,1), (2,2,2)
No fill-ins
収束性の傾向
問題 : ２次元Mooney体 (4/3c(MINI)-要素)
収束判定 : 相対残差L2ノルム＜1.e-8
反復法 ： Full-GMRES
h
1/8
1/16
1/32
1/64
不定値前処理
41
58
113
258
正定値前処理
70
127
258
541
• 反復数はともに1/hにおおよそ比例
•正定値前処理はおおよそ２倍の反復数
正定値前処理での固有値分布
収束性の理論的評価 (対称行列の場合)
CG
ek
e0
A
A
a
 min
pk , pk ( 0) 1
max
i 1,...,n
i  [ a, b] ( a  0) 
0
pk (i )
ek
e0
MINRES
rk
 min max pk (i )
r0 pk , pk (0)1 i 1,...,n
b
A
A
 b / a 1 

 2 

 b / a 1 
a
k
a=O(h2), b=O(1)
 #Iterations=O(1/h)
b
c
d
0
c=O(h2), a,b,d=O(1)
#Iterations=O(1/h)
 ad / bc 1 
rk

i  [ a, b] [ c, d ] ( b  0  c) 
 2
 ad / bc  1 
r0


k /2
収束性の関連性
 A BT  I 0 A BT  I 0  I 0  A BT 








0

I

B
C
B

C
B

C


 
 0 1I 0 1I 

k
 A B  I 0 I 0  A B  I 0 




 


B
C
B

C

 0 1I0 1I 
 0 1I
T
k
T
I 0   A 1BT 



0 1I  1B C 
 A 1BT 


 1B C 
A BT 


B

C


k
I 0 


0 1I
I 0 


0 1I
1
1
収束性を関連付けられるか？
数値実験では、右の方が２倍速い。
まとめ
• 変位部剛性行列の安定性について
初期応力項が不安定性を引き起こす可能性あり
非圧縮性制約条件が大変形に対する安定性に貢献？
Eを基準としない妥当な弾性ポテンシャルの定義法はあるか？
• 線形化方程式の解法に関して
行列のブロック構造および符号を考慮して前処理行列を構築す
ることが重要。
収束性の理論的評価の残された問題。