静弾性解析 1. 定式化と離散化の概要 1.1 線形弾性体の定式化 Fig.1 に示される線形弾性体の境界値問題を考える. ただし, 微小変形を仮定する.Fig.1 N において,N を次元数とすると, は有界領域であり, はその境界である. ここで, d は変位境界条件が与えられる境界, t は応力境界条件が与えられる境界である. d と t の間には, d および t d の関係が成り立つとする. t Linar lastic matrial Prscribd displacmnt u 0 t t Surfac forc (traction) d b Body forc d t Fig.1 Boundary valu problm of linar lastic matrial 1 関数空間 V { v vh ( ) N, vu on d} を考えると, 線形弾性体の境界値問題は [B] のように記述される. [B] 以下を満たすような変位 u V を求めよ. ( 平衡方程式 ) b 0 in (1.a) ( 応力 -ひずみの関係式) 1
W = C: (tr ) (1.b) ( ひずみ- 変位の関係式 ) 1 u( u) (1.c) 2 ( 境界条件式 ) uu on (1.d) d t n (1.) t on t ただし, 右上添え字 はテンソルの転置であることを意味する. また, は座標 x でのナブラ, は応力テンソル, は密度, bは単位質量当たりの体積力, は微小ひずみテン ソル, I は恒等テンソル,W は弾性ポテンシャル関数,C は 4 階の弾性テンソル,Young 率 E と Poisson 比 を用いると, Lamé 定数, は E (1 )(12 ) E 2(1 ) である. 次のエネルギー最小問題 [M] を考える. [M] 以下を満たすような変位 u V を求めよ. (2.a) (2.b) ( u) ( v) v V (3.a) ( v) W( v) d t v d bv d (3.b) t 1 N 関数空間 M { u uh ( ), u0 on } を定義すると, v uuとなる変分量 um が存在する. ここで, 変分量 d ( uu) ( u ) (4) と停留条件 を考えると, 0 (5) W : d d d t u b u t : d t u d bu d (6) 0 t が得られる. よって, エネルギー最小問題 [M] は, 以下の仮想仕事の原理 [V] と等価で 2
ある. [V] 以下を満たすような変位 u V を求めよ. : d t u d bu d um (7) t ただし, 1 ( ) 2 u u (8) である. さらに, 式 (8) を式 (7) の左辺に代入すると 1 : d : ( ) 2 u u :( u) d ( u) d ( ) ud ( n) ud ( ) ud ( n) ud ( ) ud t と式変形できるため, b in (10.a) n t on (10.b) t が得られる. よって, 仮想仕事の原理 [V] と境界値問題 [B] は等価である. 以上より, 境界値問題 [B], エネルギー最小問題 [M], 仮想仕事の原理 [V] は等価である. 仮想仕事の原理 [V] は境界値問題 [B] より解の微分可能性に対する要求が弱くなるため, 弱形式と呼ばれる. それに対して, 境界値問題 [B] は強形式と呼ばれる. d (9) 1.2 線形弾性体の離散化 有限要素法による離散化式は, 式 (7) の仮想仕事の原理 [V] を d d (11) のような有限要素 によって分割することによって求まる. 以下では,Fig.2 に示されるような 3 次元六面体 8 節点要素を考える. 各要素での節点数を 8 であるため, 各要素での補間関数と写像関数は 8 ( ) ( ) u N u (12.a) 1 8 ( ) ( ) x N x (12.b) 1 3
のように書くことができる. (5) (8) (6) (7) z 2 3 1 y x (4) g 3 g 2 g 1 (1) (2) (3) g 1 x, g = 2 x, g x 3 : Covariant basis vctors 1 1 1 1 8 : Shap function ( ) ( ) ( ) ( ) N ( ) 1, 2,..., 8 1, 1, 1, 1, 1, 1, 1, 1 ( ) 1, 2,..., 81, 1, 1, 1, 1, 1, 1, 1 ( ) 1, 2,..., 81, 1, 1, 1, 1, 1, 1, 1 Fig.2 8-nod hxahdral lmnt 式 (7) の左辺は, 式 (11) より : d : d 2 2 2 xx yy zz xy yz zx σ σ σ σ σ σ xx yy zz xy yz zx d (13) 4
のように式変形できる. 式 (1.b) より σ xx xx σ yy yy σ zz zz D σ xy 2 xy σ yz 2 yz σ zx 2 zx の関係が成り立つ. ただし, D マトリックスは 2 0 0 0 2 0 0 0 2 0 0 0 D = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 (14) (15) である. また, 式 (1.c) より xx yy zz Bu 2 xy 2 yz 2 zx の関係が成り立つ. ただし, B マトリックスは ( ) N 0 0 x ( ) N 0 0 y N 0 0 ( ) z ( ) (1) (2) ( ) (8) B, B B B B B ( ) ( ) N y 0 N z N x N z 0 ( ) ( ) ( ) ( ) 0 N y N x (16) (17) であり, 列ベクトル u は 5
(1) (1) (1) (2) (2) (2) ( ) ( ) ( ) (8) (8) (8) u u x uy uz ux uy uz u x u y u z ux uy u (18) z である. さらに, 式 (12.a) より ux uy u z Nu (19) の関係が成り立つ. ただし, N マトリックスは ( ) N 0 0 0 0 N ( ) ( ) (1) (2) ( ) (8) N 0 N 0, N N N N N ( ) である. 式 (13) は, 式 (14) と式 (16) より (20) : d 2 2 2 δ ( B δ u ) D ( B u ) u B D B d u δu K u xx yy zz xy yz zx d xx yy zz D 2 xy 2 yz 2 zx d (21) のように式変形できる. ただし, K は要素剛性マトリックスであり, K B D B d (22) である. 全体剛性マトリックスを K とすると, 式 (21) は δ δ u K u u K u (23) のように書くことができる. ただし, 総節点数を N とする列ベクトル u は 1 1 1 2 2 2 I I I N N N x y z x y z x y z x y z u u u u u u u u u u u u u (24) である. 式 (18) は各要素で付けられた通し番号であるが, 式 (24) は全体で付けられた通し番号である. 式 (7) の右辺は, 式 (19) より 6
t b d d t d b d x x t u b u ( ) ( ) y y N u N u t z b z t b u N t d u N b d u x x y y t z b z のように式変形できる. ただし, f は要素外力ベクトルであり, tx bx f N t y d b y d N t z b z f である. 全体外力ベクトルをf とすると, 式 (25) は δu f δu f (27) のように書くことができる. 式 (23) と式 (27) より, 有限要素法による離散化式 Kuf (28) を得ることができる. (25) (26) 1.3 非適合要素による線形弾性体の離散化 Wilson ら (1973) が提案した要素 [1] は, 各要素で曲げ変形を表現可能な関数を式 (12.a) の補間関数に付加した要素である. このような要素は, 要素境界で関数が不連続となり, 非適合要素と呼ばれる. 以下では,Wilson ら (1973) が提案した非適合要素を用いて, 離散化式を求める. 式 (19) を次式のように置き換える. ただし, u x uy Nu Pa (29) u z 7
P P P 1 1, 1 (1) 2 (2) 2 (3) 2 P 0 0 ( ) ( ) P 0 ( ) P 0, (1) (2) (3) ( ) 0 0 P P P P P (30) であり, 列ベクトル a は曲げ変形を表現するために追加される自由度 (1) (1) (1) (2) (2) (2) (3) (3) (3) a a x ay az ax ay az ax ay a (18) z である. すると, 式 (16) は xx yy zz Bu Ga (31) 2 xy 2 yz 2 zx のように変更される. ただし, ( ) P 0 0 x ( ) P 0 0 y P 0 0 ( ) z ( ) (1) (2) (3) G, ( ) ( ) P y 0 P z P x P z 0 ( ) ( ) ( ) ( ) 0 P y P x G G G G (32) である. 式 (13) は, 式 (14) と式 (31) より 8
: d 2 2 2 ( B δu G δ a ) D ( B u G a ) δu B D B d u B D G d a δa G D B d u G D G d a xx yy zz xy yz zx d δu Kuu Kua u δa Kau K aa a xx yy zz D d 2 xy 2 yz 2 zx のように式変形できる. ただし, K uu, K ua, K au, K aa は要素剛性マトリックスであり, (33) Kuu B D B d Kua B D G d Kau G D B Kaa G D G d d (34) である. 式 (33) と式 (35) より, 有限要素ごとの離散化式 Kuu Kua u f u Kau Kaa a 0 (35) を得ることができる. 式 (35) より a K K u (36) 1 aa au であるから, 列ベクトルa を消去することができる. この式変形は, 静的縮約と呼ばれる. すると, δ u K u δu f (37) となり, 式 (28) が得られる. ただし, である. K K K K K (38.a) 1 uu ua aa au f f u (38.b) 9
要素形状が平行六面体でない場合 ( 要素内で共変基底ベクトルが一定でない場合 ), 一定応力状態でa 0となるため,Wilson ら (1973) が提案した非適合要素はパッチテストを通らない. そこで, aylorら (1976) は G = O (39) d のような条件を課すことによって, パッチテストを通すことを可能にしている [2]. FrontISR では, 微小変形弾性解析の場合, 式 (39) の条件を課した Wilson らの非適合要素を採用している. 2. 各サブルーチンの解説 微小変形静弾性解析において,3 次元六面体 8 節点要素を使用する場合を考える. D マトリックスを計算するサブルーチンは lib/physics/elasticlinar.f90 の m_elasticlinar :: calelasticmatrix() と lib/physics/calmatmatrix.f90 の m_matmatrix :: MatlMatrix(), B マトリックスを計算するサブルーチンは lib/static_lib_3dic.f90 の m_static_lib_3dic :: SF_C3D8IC(), 要素剛性マトリックス K の計算を行うサブルーチンは,analysis/static/static_ mat_ass_main.f90 の m_static_mat_ass_main :: FSR_MA_ASS_MAIN() と m_static_mat_ ass_main :: FSR_LOCAL_SF_CREAE() である. 以下では, これらのサブルーチンについて解説する. 2.1 lib/physics/elasticlinar.f90 の m_elasticlinar :: calelasticmatrix() lib/physics/elasticlinar.f90 のモジュール m_elasticlinar の概要を Fig.3 に示す. また, モジュール m_elasticlinar のメンバであるサブルーチン calelasticmatrix() の概要を Fig.4 に示す. 2.2 lib/physics/calmatmatrix.f90 の m_matmatrix :: MatlMatrix() lib/physics/ calmatmatrix.f90 のモジュール m_matmatrix の概要を Fig.5 に示す. また, モジュール m_matmatrix のメンバであるサブルーチン MatlMatrix() の概要を Fig.6 に示す. 2.3 lib/static_lib_3dic.f90 の m_static_lib_3dic :: SF_C3D8IC() 10
lib/static_lib_3dic.f90 のモジュール m_static_lib_3dic の概要を Fig.7 に示す. また, モジュール m_static_lib_3dicのメンバであるサブルーチン SF_C3D8IC() の概要を Fig.8に示す. 2.4 analysis/static/static_mat_ass_main.f90 の m_static_mat_ass_main :: FSR_MA_ASS_MAIN() と m_static_mat_ass_main :: FSR_LOCAL_SF_CREAE() analysis/static/static_mat_ass_main.f90 のモジュール m_static_mat_ass_main の概要を Fig.9 に示す. また, モジュール m_static_mat_ass_main のメンバであるサブルーチン FSR_MA_ASS_MAIN() の概要を Fig.10 に示し, サブルーチン FSR_LOCAL_ SF_CREAE() の概要を Fig.11 に示す. 参考文献 [1] Wilson, E.L., aylor, R.L., Dohrty, W.P. and Ghaboussi, J., Incompatibl displacmnt modls, Numrical and Computr Mthods in Structural Mchanics (d. Fnvs, S.. t al.), pp.43-57, (1973). [2] aylor, R.L., Brsford, P.J. and Wilson, E.L., A non-conforming lmnt for strss analysis, Intrnational Journal for Numrical Mthods in Enginring, Vol.10, pp.1211-1219, (1976). モジュール名 :m_elasticlinar 線形弾性体の D マトリックスを計算するモジュール 使用する他のモジュール mmatrial 材料物性の情報を管理するモジュール メンバ変数 整数型 kral 実数型の種別値 11
メンバ関数 サブルーチン calelasticmatrix() 3 次元問題, 平面ひずみ問題, 平面応力問題, 軸対称問題の D マトリックスを計算するサブルーチン サブルーチン calelasticmatrix_ortho() 直交異方性がある場合,3 次元問題の D マトリックスを計算するサブルーチン サブルーチン LinarElastic_Shll() シェル要素を使用する場合, 埋め込み座標系成分の D マトリックスを計算するサブルーチン Fig.3 Collaborator and Rsponsibility of m_elasticlinar サブルーチン名 :calelasticmatrix() 3 次元問題, 平面ひずみ問題, 平面応力問題, 軸対称問題の D マトリックスを 計算するサブルーチン 引数 構造体(tMatrial) matl 材料に関連するデータ 整数型 sctyp 問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) 実数型 D(:, :) Dマトリックスの成分 12
実数型 tmp ( 省略可能 ) 温度 上位 サブルーチン m_matmatrix :: MatlMatrix() サブルーチン mcrp :: iso_crp() サブルーチン m_elastoplastic :: calelastoplasticmatrix() サブルーチン mviscoelastic :: calviscolasticmatrix() 下位 サブルーチン:m_tabl :: ftch_abldata() Fig.4 Argumnts and associatd subroutins of m_elasticlinar :: calelasticmatrix() モジュール名 :m_matmatrix 各材料の D マトリックスを計算するサブルーチンを呼ぶモジュール 使用する他のモジュール mmatrial 材料物性の情報を管理するモジュール mmchgauss Gauss 積分点の情報を管理するモジュール m_elasticlinar 線形弾性体の D マトリックスを計算するモジュール mhyprelastic 13
超弾性体の 4 階の弾性テンソルC を計算するモジュール m_elastoplastic 弾塑性体の D マトリックスを計算するモジュール mviscoelastic 粘弾性体の D マトリックスを計算するモジュール mcrp クリープを考慮した剛性マトリックス K を計算するためのモジュール muelastic ユーザ定義の弾性体の D マトリックスを計算するモジュール mumat ユーザ定義の材料の D マトリックスを計算するモジュール メンバ変数 整数型 kral 実数型の種別値 メンバ関数 サブルーチン gtnlgomflag() 未使用のサブルーチン サブルーチン MatlMatrix() 各材料の D マトリックスを計算するサブルーチンを呼ぶサブルーチン サブルーチン StrssUpdat() 各材料の応力とひずみを計算するサブルーチンを呼ぶサブルーチン サブルーチン mat_c2d() 材料が超弾性体の場合,4 階の弾性テンソルを問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) に応じた D マトリックスに変換するサブルーチン サブルーチン MatlMatrix_Shll() シェル要素を使用する場合, 各材料 ( 現バージョンでは, 線形弾性体のみ ) の応力とひずみを計算するサブルーチンを呼ぶサブルーチン サブルーチン mat_c2d_shll() シェル要素を使用する場合,4 階の弾性テンソルを D マトリックスに変換するサブルーチン Fig.5 Collaborator and Rsponsibility of m_matmatrix 14
サブルーチン名 :MatlMatrix() 各材料の D マトリックスを計算するサブルーチンを呼ぶサブルーチン 引数 構造体(tGaussStatus) gauss Gauss 積分点に関連するデータ 整数型 sctyp 問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) 実数型 matrix(:, :) Dマトリックスの成分 実数型 dt 時間増分 15
実数型 cdsys(3, 3) 直交異方性がある場合に使用する座標系 実数型 tmpratur ( 省略可能 ) 温度 上位 サブルーチン m_static_lib_2d :: SF_C2() サブルーチン m_static_lib_2d :: UPDAE_C2() サブルーチン m_static_lib_2d :: UpdatS_C2() サブルーチン m_static_lib_3d :: SF_C3() サブルーチン m_static_lib_3d :: LOAD_C3 () サブルーチン m_static_lib_3d :: UPDAE_C3() サブルーチン m_static_lib_3d :: UpdatS_C3() サブルーチン m_static_lib_3dic :: SF_C3D8IC() サブルーチン m_static_lib_3dic :: UpdatS_C3D8IC() サブルーチン m_static_lib_c3d8 :: SF_C3D8Bbar() サブルーチン m_static_lib_c3d8 :: Updat_C3D8Bbar() サブルーチン m_static_lib_c3d8 :: LOAD_C3D8Bbar() 下位 サブルーチン m_matmatrix :: mat_c2d() サブルーチン muelastic :: uelasticmatrix() サブルーチン mviscoelastic :: calviscolasticmatrix() サブルーチン m_elasticlinar :: calelasticmatrix() サブルーチン m_elasticlinar :: calelasticmatrix_ortho() サブルーチン mhyprelastic :: calelasticmoonyrivlin() サブルーチン mhyprelastic :: calelasticarrudaboyc() サブルーチン m_elastoplastic :: calelastoplasticmatrix() サブルーチン mumat :: umatlmatrix() サブルーチン mcrp :: iso_crp() Fig.6 Argumnts and associatd subroutins of m_matmatrix :: MatlMatrix() 16
モジュール名 :m_static_lib_3dic 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,B マトリックスおよび要素剛性マトリックス K を計算したり,Gauss 積分点における応力とひずみを計算したりするモジュール 使用する他のモジュール hcmw HECMWのモジュール m_utilitis 補助的なサブルーチンや関数を集めたモジュール lmntinfo 要素の情報を管理するモジュール mmchgauss 17
Gauss 積分点の情報を管理するモジュール m_matmatrix 各材料の D マトリックスを計算するサブルーチンを呼ぶモジュール メンバ変数 整数型 kral 実数型の種別値 メンバ関数 サブルーチン SF_C3D8IC() 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,B マトリックスおよび要素剛性マトリックス K を計算するサブルーチン サブルーチン UpdatS_C3D8IC() 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,Gauss 積分点における応力とひずみを計算するサブルーチン Fig.7 Collaborator and Rsponsibility of m_static_lib_3dic サブルーチン名 :SF_C3D8IC() 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,B マトリックスおよび要素剛性マトリックス K を計算するサブルーチン 引数 整数型 typ 要素タイプ 整数型 nn 各要素の節点数 (nn=8) 実数型 coord(3, nn) 各要素の節点座標 構造体(tGaussStatus) gausss(:) Gaussの積分点に関連するデータ 実数型 stiff(:, :) 18
要素剛性マトリックス K 実数型 nddisp(3, nn) 各要素の節点変位 実数型 hdisp(3, 3) u 各要素の節点自由度 a 上位 サブルーチン m_fstr_updat :: fstr_updat3d() サブルーチン m_static_mat_ass_main :: FSR_LOCAL_SF_CREAE() 下位 サブルーチン:lmntInfo :: gtjacobian サブルーチン:m_MatMatrix :: MatlMatrix () サブルーチン:lmntInfo :: gtquadpoint () サブルーチン:lmntInfo :: gtglobaldriv () サブルーチン:lmntInfo :: gtwight () サブルーチン:m_utilitis :: calinvrs () Fig.8 Argumnts and associatd subroutins of m_static_lib_3dic :: SF_C3D8IC() モジュール名 :m_static_mat_ass_main 全体剛性マトリックス K を計算するモジュール 使用する他のモジュール m_fstr FrontISRの基本情報を管理するモジュール m_static_lib 静解析で必要となるモジュール群 ( 補助的なサブルーチンや関数を集めたモジュール,B マトリックスおよび要素剛性マトリックス K を計算したり,Gauss 積分点における応力とひずみを計算したりするモジュール, 線形ソルバーの情報を管理するモジュール ) を使用するモジュール mmchgauss Gauss 積分点の情報を管理するモジュール 19
メンバ変数なし メンバ関数 サブルーチン FSR_MA_ASS_MAIN() 全体剛性マトリックス K を計算するサブルーチン サブルーチン FSR_LOCAL_SF_CREAE() 各要素タイプの要素剛性マトリックス K を計算するサブルーチンを呼ぶサブルーチン Fig.9 Collaborator and Rsponsibility of m_static_mat_ass_main サブルーチン名 :FSR_MA_ASS_MAIN() 全体剛性マトリックス K を計算するサブルーチン 引数 構造体(hcmwS_matrix) hcmesh HECMWが管理するメッシュのデータ 構造体(hcmwS_local_msh) hcma HECMWが管理するマトリックスのデータ 構造体(fstr_solid) fstrsolid FrontISRの静解析用データ 上位 サブルーチン fstr_solv_dynamic_xplicit() 20
サブルーチン サブルーチン fstr_solv_dynamic_implicit() m_static_mat_ass :: fstr_mat_ass() 下位 サブルーチン:hcmw_mat_clar() サブルーチン:m_static_mat_ass_main :: fstr_local_stf_crat() サブルーチン:hcmw_mat_ass_lm() Fig.10 Argumnts and associatd subroutins of m_static_mat_ass_main :: FSR_MA_ ASS_MAIN() サブルーチン名 :FSR_LOCAL_SF_CREAE() 各要素タイプの要素剛性マトリックス K を計算するサブルーチンを呼ぶサブルーチン 引数 構造体(hcmwS_matrix) hcmesh HECMWが管理するメッシュのデータ 整数型 ndof 各節点の自由度数 整数型 ic_typ 要素タイプ 整数型 icl 要素番号 21
実数型 xx(:), yy(:), zz(:) 各要素の節点座標 構造体(tGaussStatus) gausss(:) Gaussの積分点に関連するデータ 整数型 ist 問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) 実数型 stiffnss(:, :) 要素剛性マトリックス K 上位 サブルーチン m_static_mat_ass_main :: FSR_MA_ASS_MAIN() 下位 サブルーチン m_static_lib_2d :: SF_C2() サブルーチン m_static_lib_1d :: SF_C1() サブルーチン m_static_lib_3dic :: SF_C3D8IC() サブルーチン m_static_lib_3d :: SF_C3() サブルーチン m_static_lib_shll :: SF_Shll_MIC() サブルーチン m_static_lib_bam :: SF_Bam() サブルーチン hcmw_abort() Fig.11 Argumnts and associatd subroutins of m_static_mat_ass_main :: FSR_LOCAL_ SF_CREAE 22