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

Presentation on theme: "超弾性体解析で現れる剛性行列の性質とその解法に関して"— Presentation transcript:

1 超弾性体解析で現れる剛性行列の性質とその解法に関して
鷲尾　巧　　久田　俊明 東京大学 日本応用数理学会「行列・固有値の解法とその応用」研究部会 第６回研究会

2 本発表の概要 Aが有限要素離散化においてどのようなプロセスを経て構成されるか Aの正定値性に関する議論
全体行列がどのようなプロセスを経て構成されるか 反復解法における前処理の構成法について 収束性について ほとんど話題は、不思議に思うこと、自分ではわかっていないことです!!

3 ゴムと金属の違い 鉄鋼 2x106 ゴム 1～3 鋳鉄 0.6x10－11 ゴム 53.7x10－11 金属 1%以下 ゴム 1000%
ヤング率　(MPa) 圧縮率　(Pa－1) 弾性変形の限界の伸び 鉄鋼 2x106 ゴム 1～3 鋳鉄 0.6x10－11 ゴム 53.7x10－11 金属 1%以下 ゴム 1000% ゴム構造 自然状態 伸ばす ヘルムホルツの自由エネルギー エントロピーSの減少 結晶構造 ポテンシャルEの増加

4 統計力学的ミクロモデルにおける非圧縮性の関与 （参照：ゴムの弾性（裳華房） 久保亮五 著）
統計力学的ミクロモデルにおける非圧縮性の関与 （参照：ゴムの弾性（裳華房）　久保亮五　著） 非圧縮的伸長 高さ方向のエントロピー 伸ばす １ 水平方向のエントロピー １ １ 応力と歪みの関係式

5 連続体の変形を表すテンソル 変形勾配テンソル 右極分解 本質的な変形を取り出す Euler座標 （現配置での位置を表す）
Lagrange座標 （変形前の位置を表す） 右極分解　本質的な変形を取り出す R: 回転行列, U:正定値対称行列

6 歪みテンソル 右Cauchy-Green変形テンソル Green-Lagrangeひずみ 圧縮 引張り せん断

7 弾性ポテンシャルWとその変分 第１変分式 第２変分式 歪みEを与えるために要する単位体積あたりのエネルギー
第２Piola-Kirchhoff応力テンソル 第２変分式 初期変位項 初期応力項

8 応力テンソルの物理的意味 現配置 参照配置

9 剛性行列の正定値性 初期変位項 : Wが凹関数ならば正定値 初期応力項 : Sが正定値ならば正定値 初期応力項の正定値性は不確定：

10 主不変量 等方的連続体の弾性ポテンシャルは主不変量の関数として表わされる。

11 不変量から定まる剛性行列の符号 第１不変量 初期応力項：正定値 初期変位項：ゼロ 第２不変量 初期応力項：正定値 初期変位項：ゼロ３負３
第１不変量　　初期応力項：正定値　　初期変位項：ゼロ 第２不変量　　初期応力項：正定値　　初期変位項：ゼロ３負３ 第３不変量　　初期応力項：正定値　　初期変位項：ゼロ３負３ 第２、３不変量の剛性行列（初期応力項＋初期変位項）の符号判定は困難

12 ゴム状物質の弾性ポテンシャル例 Mooney-Rivlin体 あるいは など 非圧縮性制約条件付き ペナルティ方的なポテンシャル
第２不変量の凹関数性を保証できないため、Newton-Raphson解法の収束性を 一般に保証できない。

13 ポテンシャルWのLandscape(想像）
境界および制約条件を満たす解の集合 初期状態

14 心筋の弾性ポテンシャルの例 sheet fiber
Usyk, Mazhari and McCulloch, Journal of Elasticity 61,2000 sheet fiber 自然基底 異方性に従った基底 初期応力項：　歪みに依存　　　初期変位項：　正定値

15 異方性連続体モデルの不安定性 異方性連続体モデルは、初期変位項が正定値となるように設計される。
　異方性連続体モデルは、初期変位項が正定値となるように設計される。 　しかし、安定性議論においては初期応力項の影響が無視されているようである。 　異方性モデルは、変形過程の初期において、解けなくなるケースが多々生じる。 　安定な第１不変量を少し加えることにより、不安定性を回避できる。 非線形 （２次項あり） 凹関数 線形 変位 変位勾配テンソル 歪みテンソル 弾性ポテンシャル 　歪みテンソルが変位の２次関数となっていることが問題。 　しかし、回転不変性などを考慮しつつ弾性ポテンシャルを変位勾配Zの凹関数として直接的に定義することは困難。

16 体積弾性の取扱について 未定乗数の導入 第１変分式　（制約条件式も付加） 第２変分式　（制約条件式も付加）

17 不定値か正定値悪条件か？ 未定乗数導入時 不定値 剛性行列 直接的変分 正定値　しかし悪条件 剛性行列

18 固有値分布 (1) #positive eigenvalues = dim (A)
#negative eigenvalues =dim(C) a b c d

19 固有値分布 (2) 固有値は右半面に含まれる

20 ブロックタイプ前処理行列 良いQB の構築は簡単ではない。

21 ブロックタイプ前処理の効果

22 正定値ブロックタイプ前処理 MINRESなどの対称行列向き（しかし非正定値）反復法が適用可能

23 ＩＬＵ分解とSchur complement
d u l Q L, D およびUを足し合わせ　Q=L+D+Uとおく。 d u l Q－l u/d on P Q on PC ILU分解の安定性 対角行列Dの符号が完全分解時の符号に従うことが望ましい 不完全なSchur complementを計算

24 より実践的な前処理（ILU前処理） 変数分離のみで純代数的に構成できる。 フィルインの設定に変数分離情報を利用する。 j
k j (bi,bk,bj)=(block(i),block(k),block(j)) (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

25 収束性の傾向 h 1/8 1/16 1/32 1/64 不定値前処理 41 58 113 258 正定値前処理 70 127 541
問題 : ２次元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 541 反復数はともに1/hにおおよそ比例 正定値前処理はおおよそ２倍の反復数

26 正定値前処理での固有値分布

27 収束性の理論的評価 (対称行列の場合) CG MINRES a=O(h2), b=O(1)  #Iterations=O(1/h)
a=O(h2), b=O(1)  #Iterations=O(1/h) c d a b MINRES c=O(h2), a,b,d=O(1) 　#Iterations=O(1/h)

28 収束性の関連性 収束性を関連付けられるか？ 数値実験では、右の方が２倍速い。

29 まとめ 変位部剛性行列の安定性について 初期応力項が不安定性を引き起こす可能性あり 線形化方程式の解法に関して
　　初期応力項が不安定性を引き起こす可能性あり 　　非圧縮性制約条件が大変形に対する安定性に貢献？ 　　Eを基準としない妥当な弾性ポテンシャルの定義法はあるか？ 線形化方程式の解法に関して 　行列のブロック構造および符号を考慮して前処理行列を構築することが重要。 　収束性の理論的評価の残された問題。