FrontISTR に実装されている定式化を十分に理解し, 解きたい問題に対してソースコードを自由にカスタマイズ ( 要素タイプを追加, 材料の種類を追加, ユーザサブルーチンを追加 ) できるようになること を最終目標とします 第 3 回, 第 7 回, 第 10 回の研究会では,FrontIST

FrontISTR による弾性弾性解析 ( 直交異方弾性体 ) 東京大学新領域創成科学研究科人間環境学専攻橋本学 014 年 7 月 30 日第 11 回 FrontISTR 研究会 < 機能 例題 定式化 プログラム解説編 弾性解析 ( 直交異方弾性体を中心に ) >

FrontISTR に実装されている定式化を十分に理解し, 解きたい問題に対してソースコードを自由にカスタマイズ ( 要素タイプを追加, 材料の種類を追加, ユーザサブルーチンを追加 ) できるようになること を最終目標とします 第 3 回, 第 7 回, 第 10 回の研究会では,FrontISTR に実装されている弾性解析 ( 等方弾性体 ) の定式化, ソースコードの関連するサブルーチンについて紹介しました 第 3 回 FrontISTR 研究会プログラミング編,013/5/ 開催 第 7 回 FrontISTR 研究会産業応用事例, 有限変形定式化, ユーザーの声への対応編, 013/1/3 開催 第 10 回 FrontISTR 研究会有限変形定式化と実装,Ver.4.3 公開編,014//1 開催 今回は,FrontISTR に実装されている直交異方弾性体に焦点を当てます

線形弾性体 微小変形 ( 微小変位 ) 有限変形 ( 有限変位 ) 微小ひずみ 微小ひずみ 大ひずみ 弾塑性体粘弾性体線形弾性体粘弾性体弾塑性体超弾性体 直交異方性がある場合を扱います ( 講演では, 微小変形理論の場合を説明します ) 有限変形大ひずみ 1 0 0 T 0 0 T E= { ( u ) + ( u ) + ( u ) ( u ) } t 0S= f( t 0E, t 0E t 0E,...) ひずみ変位こう配の 次項がある応力ひずみの 次以上の項がある t t t t t 0 3

目次 前半 解析機能/ サンプル例題 1. 解析機能とユーザマニュアル該当箇所. メッシュファイルと解析制御ファイルの設定方法 注意点 3. サンプル例題 ( 内圧を受ける血管を模擬した円管モデル ) 後半 定式化/ プログラム 4. 直交異方弾性体の有限要素法定式化 5. プログラム説明 4

目次 前半 解析機能/ サンプル例題 1. 解析機能とユーザマニュアル該当箇所. メッシュファイルと解析制御ファイルの設定方法 注意点 3. サンプル例題 ( 内圧を受ける血管を模擬した円管モデル ) 後半 定式化/ プログラム 4. 直交異方弾性体の有限要素法定式化 5. プログラム説明 5

等方性 (isotropy) 弾性定数や線膨張係数があらゆる方向で等しい 異方性 (anisotopy) 弾性定数や線膨張係数が方向によって異なる - 特に, 直交する三つの軸の方向で異なる性質が直交異方性 (orthotropy) である 3 e z O z e y e x y x x e 3 e e 1 1 e, e, e 1 3 : 材料で定義される直交基底ベクトル 直交異方性材料の代表的な例 - 木のように繊維方向がある材料 - ある方向に補強材を入れた材料 - 血管壁 ( 径方向, 周方向, 長さ方向 ) など 6

応力とひずみ 3 e z O z e y e x y x x e 3 e e 1 1 e 1, e, e 3 : 材料で定義される直交基底ベクトル 応力 σ = σ e e + σ e e + σ e e xx x x xy x y zx x z + σ e e + σ e e + σ e e xy y x yy y y yz y z + σ e e + σ e e + σ e e zx z x yz z y zz z z = σ e e ひずみ ij i j ε = ε e e + ε e e + ε e e xx x x xy x y zx x z + ε e e + ε e e + ε e e xy y x yy y y yz y z + ε e e + ε e e + ε e e zx z x yz z y zz z z = ε e e ij i j (1.1) (1.) 直交異方弾性体の構成方程式 σ = C ε ij ijkl kl 弾性定数 (1.3) 7

等方弾性体の構成方程式 弾性定数があらゆる方向で等しい 1 ν ν 0 0 0 E E E ν 1 ν 0 0 0 εxx σ E E E xx ε yy 1 σyy ν ν 0 0 0 ε zz E E E σ zz = ε 1 xy 0 0 0 0 0 σxy ε µ yz σyz 1 ε 0 0 0 0 0 zx σ zx µ 1 0 0 0 0 0 ε µ σ ひずみ S Complianceに相当 応力 σxx λ + µ λ λ 0 0 0 εxx σ yy 0 0 0 εyy λ λ + µ λ σ zz λ λ λ + µ 0 0 0 ε zz = σxy 0 0 0 µ 0 0 εxy σ 0 0 0 0 µ 0 yz εyz σ 0 0 0 0 0 µ zx ε zx σ 応力 D = S 1 D マトリックス Stiffness に相当 ε ひずみ (1.4) (1.5) E ν : Young 率 [Pa] : Poisson 比 [-] λ= Eν E µ = G= (1 + ν )(1 ν ) (1 + ν ) : Lame 定数 [Pa] 8

直交異方弾性体の構成方程式 (1) 弾性定数が直交する三つの軸の方向で異なる 1 ν1 ν13 0 0 0 E1 E1 E 1 ν1 1 ν3 0 0 0 ε 11 E E E σ 11 ε ν31 ν3 1 σ 0 0 0 ε 33 E3 E3 E 3 σ 33 = ε 1 1 σ1 0 0 0 0 0 ε G 3 1 σ 3 ε 1 31 σ 0 0 0 0 0 31 G 3 ε 1 σ 0 0 0 0 0 G 31 ひずみ 応力 E1, E, E3,,,,, : Young 率 [Pa] ν ν ν ν ν ν 1 13 1 3 31 3 S Compliance に相当 : Poisson 比 [-] (1.6) G1, G3, G 31 : 横弾性定数 [Pa] ν1 ν1 = E1 E ν 3 ν3 = E E3 ν 31 ν13 = E E 3 1 9

直交異方弾性材料の構成方程式 () 弾性定数が直交する三つの軸の方向で異なる σ E1(1 ν3ν3) E1( ν31ν3+ ν1) E1( ν1ν3+ ν31) 0 0 0 ε 11 11 E1( ν31ν3+ ν1) E(1 ν13ν31) E( ν1ν31+ ν3) σ 0 0 0 ε σ 33 ε33 = E1( ν1ν3+ ν31) E( ν1ν31+ ν3) E3(1 ν1ν1) σ 0 0 0 1 ε 1 σ 3 0 0 0 G1 0 0 ε 3 σ 31 0 0 0 0 G3 0 ε 31 0 0 0 0 0 G 31 σ 応力 E1, E, E3,,,,, D = S = 1 ν ν ν ν ν ν ν ν ν : Young 率 [Pa] ν ν ν ν ν ν 1 13 1 3 31 3 D マトリックス Stiffness に相当 1 1 3 3 31 13 1 3 13 : Poisson 比 [-] 1 ε ひずみ (1.7) (1.8) G1, G3, G : 横弾性定数 [Pa] 31 ν1 ν1 = E1 E ν 3 ν3 = E E3 ν 31 ν13 = E3 E1 10

直交異方弾性体の引張変形の例 E =.00 10 Pa, E = E = 1.00 10 Pa 7 7 1 3 6 1 = 13= 3= 0.30, G1 = G1 = G1 = 7.69 10 Pa ν ν ν 0.3 m 1/8 モデル 0.05 m p =.0 MPa 3 1 fixed 1 3 p =.0 MPa 1 3 p =.0 MPa 3 p =.0 MPa 1 11

FrontISTR の解析機能を確認するため,FrontISTR のユーザマニュアル ( ファイル名 FrontISTR_user_manual_Ver35.pdf ) の該当箇所を見ます FrontISTR ソースコード FrontISTR_V43_p1.tar.gz を解凍すると, ディレクトリ FrontISTR_V43 ができます FrontISTR のユーザマニュアルはディレクトリ FrontISTR_V43/ doc 内にあります FrontISTR のユーザマニュアルの 118 ページ ~10 ページに直交異方弾性体の記述があります 現在の FrontISTR のバージョンでは,!SOLUTION, TYPE=NLSTATIC のときのみ直交異方弾性体に対応しています 次にリリースされる修正版は,!SOLUTION, TYPE=STATIC でも直交異方弾性体に対応する予定です 1

FrontISTR ユーザマニュアルより (1) FrontISTR のユーザマニュアルの 10 ページ 材料定数を定義する座標系での弾性定数を入力 ν1 ν1 = E1 E ν 3 ν3 = E E3 ν 31 ν13 = E3 E1 13

FrontISTR ユーザマニュアルより () FrontISTR のユーザマニュアルの 118 ページ ~119 ページ 材料定数を定義する座標系の直交基底ベクトルを入力 セクション情報と座標系情報を結びつける I x 3 x 3 e e 3 e 1 ca e1 = ca 1 ca cb e3 = ca cb x 1 e = e e 3 1 14

目次 前半 解析機能/ サンプル例題 1. 解析機能とユーザマニュアル該当箇所. メッシュファイルと解析制御ファイルの設定方法 注意点 3. サンプル例題 ( 内圧を受ける血管を模擬した円管モデル ) 後半 定式化/ プログラム 4. 直交異方弾性体の有限要素法定式化 5. プログラム説明 15

!NODE 983, 3.51159307E-0,.0117903E+00, 0.00000000E+00 984, 3.5370485E-0,.038855E+00, 0.00000000E+00!ELEMENT, TYPE=361 160, 983, 984, 1177, 1175, 1377, 1469, 1470, 1379 159, 97, 973, 984, 983, 1375, 1468, 1469, 1377!EGROUP, EGRP=E0000160 160,!SECTION, TYPE=SOLID, EGRP=E0000160, MATERIAL=M01!EGROUP, EGRP=E0000159 159,!SECTION, TYPE=SOLID, EGRP=E0000159, MATERIAL=M01!END 16

!VERSION 3!SOLUTION, TYPE=NLSTATIC 3. 3.0.8.6.4 a!material, NAME=M01!ELASTIC, TYPE=ORTHOTROPIC 0.1, 0.1, 1.0, 0.049, 0.049, 0.049, 0.048, 0.34, 0.34..0 b c 1.8-1. -1.0-0.8-0.6-0.4-0. 0.0 0.!ORIENTATION, DEFINITION=COORDINATES, NAME=O0000160.6E-0, 3.0E+00,.5E-0, -9.8E-01,.0E+00,.5E-0, 1.8E-0,.0E+00,.5E-0,!SECTION, SECNUM=1, ORIENTATION=O0000160!ORIENTATION, DEFINITION=COORDINATES, NAME=O0000159 7.9E-0, 3.0E+00,.5E-0, -9.5E-01,.0E+00,.5E-0, 5.3E-0,.0E+00,.5E-0,!SECTION, SECNUM=, ORIENTATION=O0000159!END ca e1 = ca ca cb e3 = ca cb e = e e 3 1 17

目次 前半 解析機能/ サンプル例題 1. 解析機能とユーザマニュアル該当箇所. メッシュファイルと解析制御ファイルの設定方法 注意点 3. サンプル例題 ( 内圧を受ける血管を模擬した円管モデル ) 後半 定式化/ プログラム 4. 直交異方弾性体の有限要素法定式化 5. プログラム説明 18

( ) 本モデルおよび計算結果は, 帝京大学ジョイントプログラムセンター田沼唯士教授との共同研究の成果です 節点数 :947,583 要素数 :864,000 微小変形, 線形弾性体六面体 1 次要素 ( 非適合要素 ) 1/4 対称モデル 管の外径 :.3mm 管の内径 :.0mm y 40.0 mm (800 分割 ) z x 内圧 ( 圧力差 ) p = 0.01351335 MPa Copyright (c) 014 Teikyo University & The University of Tokyo 19

1 0.04 mm 0.015 mm 0.073 mm 0.058 mm y 周方向長さ方向 周方向 周方向 長さ方向 z x Copyright (c) 014 Teikyo University & The University of Tokyo 0

5 3 8 6 7 z e e y z e x y x e e 3 4 e 1 1 1 3 軸 : 径方向軸 : 周方向軸 : 長さ方向 1 3 E E 1 E G G G 1 3 1 3 31 方向の繊維によって剛性が大きくなる場合, 下記のようにパラメータを設定 = 1.0MPa = 0.1MPa = 0.1MPa = 0.33557MPa = 0.04766MPa = 0.33557MPa ν ν ν 1 13 3 = 0.49 = 0.49 = 0.049 : 方向に引張る場合の 方向の縮みを意味する 1 ν ν ν E = ν = 0.049 1 1 E1 E = ν = 0.049 3 31 13 E1 E = ν = 0.049 3 3 3 E : 方向に引張る場合, 方向には縮みにくい Copyright (c) 014 Teikyo University & The University of Tokyo 1 1

径方向の繊維によって剛性が大きくなる場合 ε 11 1/1.0 0.049/0.1 0.049/0.1 0 0 0 σ 11 ε 0.49/1.0 1/0.1 0.049/0.1 0 0 0 σ ε 33 0.49/1.0 0.049/0.1 1/0.1 0 0 0 σ 33 = ε 1 0 0 0 1/0.33557 0 0 σ 1 ε 3 0 0 0 0 1/0.04766 0 σ 3 ε 0 0 0 0 0 1/0.33557 31 σ 31 周方向の繊維によって剛性が大きくなる場合 ε 11 1/0.1 0.49/1.0 0.049/0.1 0 0 0 σ 11 ε 0.049/0.1 1/1.0 0.049/0.1 0 0 0 σ ε 33 0.049/0.1 0.49/1.0 1/0.1 0 0 0 σ 33 = ε 1 0 0 0 1/0.33557 0 0 σ 1 ε 3 0 0 0 0 1/0.33557 0 σ 3 ε 0 0 0 0 0 1/0.04766 31 σ 31 長さ方向の繊維によって剛性が大きくなる場合 ε 11 1/0.1 0.049/0.1 0.49/1.0 0 0 0 σ 11 ε 0.049/0.1 1/0.1 0.49/1.0 0 0 0 σ ε 33 0.049/0.1 0.049/0.1 1/1.0 0 0 0 σ 33 = ε 1 0 0 0 1/0.04766 0 0 σ 1 ε 3 0 0 0 0 1/0.33557 0 σ 3 ε 0 0 0 0 0 1/0.33557 31 σ 31 Copyright (c) 014 Teikyo University & The University of Tokyo

計算結果 ( 変位の大きさ ) u 0.366 0.644 [mm] Reference configuration y z x 変位の大きさの分布および最大値 最小値は Abaqus の結果と一致する Copyright (c) 014 Teikyo University & The University of Tokyo 3

計算結果 (Mises 応力 ) von Mises stress 0.009193 0.190 [MPa] Reference configuration y z x Mises 応力の分布および最大値 最小値は Abaqus の結果と一致する Copyright (c) 014 Teikyo University & The University of Tokyo 4

目次 前半 解析機能/ サンプル例題 1. 解析機能とユーザマニュアル該当箇所. メッシュファイルと解析制御ファイルの設定方法 注意点 3. サンプル例題 ( 内圧を受ける血管を模擬した円管モデル ) 後半 定式化/ プログラム 4. 直交異方弾性体の有限要素法定式化 5. プログラム説明 5

e z z e y e x ex, ey, ez y x 5 1 3 e e 3 e 1 g 3 6 ζ 8-node solid element 4 8 g g 1 : 全体座標系の直交基底ベクトル x x = x g =, g =, g ξ ξ ξ 1 1 3 3 : 共変基底ベクトル e, e, e 1 3 1 : 材料定数を定義する座標系の直交基底ベクトル η ξ 7 3 応力 σ = σ e e + σ e e + σ e e xx x x xy x y zx x z + σ e e + σ e e + σ e e xy y x yy y y yz y z + σ e e + σ e e + σ e e zx z x yz z y zz z z = σ e e ij i j ひずみ ε = ε e e + ε e e + ε e e xx x x xy x y zx x z + ε e e + ε e e + ε e e xy y x yy y y yz y z + ε e e + ε e e + ε e e zx z x yz z y zz z z = ε e e ij i j 直交異方弾性材料の計算に必要なパラメータは σ = C ε ij ijkl kl e, e, e 1 3 線形弾性定数 直交基底ベクトル (4.1) (4.) (4.3) 6

ひずみベクトルの変換 (1) ε = e ε e ij i j = ( e e )( e e ) ε i m j n mn = ( e e )( e e ) ε + ( e e )( e e ) ε + ( e e )( e e ) ε i x j x xx i x j y xy i x j z xz + ( e e )( e e ) ε + ( e e )( e e ) ε + ( e e )( e e ) ε i y j x yx i y j y yy i y j z yz + ( e e )( e e ) ε + ( e e )( e e ) ε + ( e e )( e e ) ε i z j x zx i z j y zy i = ( e e )( e e ) ε + ( e e )( e e ) ε + ( e e )( e e ) ε 1 + { ( e i ex )( e j ey ) + ( e i ey )( e j ex ) } εxy 1 + { ( e i ey )( e j ez ) + ( e i ez )( e j ey ) } εyz 1 + { ( e i ez )( e j ex ) + ( e i ex )( e j ez ) } εzx (4.4) z j z zz i x j x xx i y j y yy i z j z zz 7

ひずみベクトルの変換 () 1 1 1 ( e 1 ex)( e 1 ex) ( e 1 ey)( e 1 ey) ( e 1 ez)( e 1 ez) {( e 1 ex)( e 1 ey) + ( e 1 ey)( e 1 ex) } {( e 1 ey)( e 1 ez) + ( e 1 ez)( e 1 ey) } {( e 1 ez)( e 1 ex) + ( e 1 ex)( e 1 ez) } ε 11 εxx ε 1 1 1 ( e ex)( e ex) ( e ey)( e ey) ( e ez)( e ez) {( e ex)( e ey) + ( e ey)( e ex) } {( e ey)( e ez) + ( e ez)( e ey) } {( e ez)( e ex) + ( e ex)( e ez) } εyy ε 33 = 1 1 1 ε zz ε ( e 3 ex)( e 3 ex) ( e 3 ey)( e 3 ey) ( e 3 ez)( e 3 e z) {( e 3 ex)( e 3 ey) + ( e 3 ey)( e 3 ex) } {( e 3 ey)( e 3 ez) + ( e 3 ez)( e 3 ey) } {( e 3 ez)( e 3 ex) + ( e 3 ex)( e 3 ez) } 1 εxy ε 3 ( e 1 ex)( e ex) ( e 1 ey)( e ey) ( e 1 ez)( e ez) ( e 1 ex)( e ey) + ( e 1 ey)( e ex) ( e 1 ey)( e ez) + ( e 1 ez)( e ey) ( e 1 ez)( e ex) + ( e 1 ex)( e ez) ε yz ε 31 ( e ex)( e 3 ex) ( e ey)( e 3 ey) ( e ez)( e 3 ez) ( e ex)( e 3 ey) + ( e ey)( e 3 ex) ( e ey)( e 3 ez) + ( e ez)( e 3 ey) ( e ez)( e 3 ex) + ( e ex)( e3 e z) ε zx ( 3 x)( 1 x) ( 3 y)( 1 y) ( 3 z)( 1 z) ( 3 x)( 1 y) + ( 3 y)( 1 x) ( 3 y)( 1 z) + ( 3 z)( 1 y) ( 3 z)( 1 x) + ( 3 x)( 1 z) e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e e l1 m1 n1 l1m1 m1n1 nl ε 1 1 xx ε l m n lm mn nl εyy l3 m3 n3 l3m3 m3n3 n3l ε 3 zz = ll 1 m1m n1n l1m + ml 1 m1n + n1m nl 1 + l1n ε xy (4.5) ll3 mm3 nn3 lm3+ ml3 mn3+ nm3 nl3+ ln 3 εyz l3l1 m3m1 n3n1 l3m1 + m3l1 m3n1 + n3m1 n3l1 + l3n 1 ε zx R ε ε li = ( e i ex) mi = ( e i ey) ni = ( e i ez) (4.6) 8

応力ベクトルの変換 (1) σ = e σ e ij i j = ( e e )( e e ) σ i m j n mn = ( e e )( e e ) σ + ( e e )( e e ) σ + ( e e )( e e ) σ i x j x xx i x j y xy i x j z xz + ( e e )( e e ) σ + ( e e )( e e ) σ + ( e e )( e e ) σ i y j x yx i y j y yy i y j z yz + ( e e )( e e ) σ + ( e e )( e e ) σ + ( e e )( e e ) σ i z j x zx i z j y zy i = ( e e )( e e ) σ + ( e e )( e e ) σ + ( e e )( e e ) σ {( e e )( e e ) ( e e )( e e )} {( e e )( e e ) ( e e )( e e )} {( e e )( e e ) ( e e )( e e )} z j z zz i x j x xx i y j y yy i z j z zz + + σ i x j y i y j x xy + + σ i y j z i z j y yz + + σ i z j x i x j z zx (4.7) 9

応力ベクトルの変換 () σ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 11 ( e ex)( e ex) ( e ey)( e ey) ( e ez)( e ez) ( e ex)( e ey) + ( e ey)( e ex) ( e ey)( e ez) + ( e ez)( e ey) ( e ez)( e ex) + ( e ex)( e ez) σxx σ ( e e x) ( e ex) ( e ey)( e ey) ( e ez)( e ez) ( e ex)( e ey) + ( e ey)( e ex) ( e ey)( e ez) + ( e ez)( e ey) ( e ez)( e ex) + ( e ex)( e ez) σyy σ ( e 33 3 ex)( e 3 ex) ( e 3 ey)( e 3 ey) ( e 3 ez)( e 3 ez) ( e 3 ex)( e 3 ey) + ( e 3 ey)( e 3 ex) ( e 3 ey)( e 3 ez) + ( e 3 ez)( e 3 ey) ( e 3 ez)( e 3 ex) + ( e 3 ex)( e 3 ez) σ = zz σ 1 ( e 1 ex)( e ex) ( e 1 ey)( e ey) ( e 1 ez)( e ez) ( e 1 ex)( e ey) + ( e 1 ey)( e ex) ( e 1 ey)( e ez) + ( e 1 ez)( e ey) ( e 1 ez)( e ex) + ( e 1 ex)( e ez) σxy σ 3 ( e ex)( e 3 ex) ( e ey)( e 3 ey) ( e ez)( e 3 ez) ( e ex)( e 3 ey) + ( e ey)( e 3 ex) ( e ey)( e 3 ez) + ( e ez)( e 3 ey) ( e ez)( e 3 ex) + ( e ex)( e 3 ez) σyz σ 31 ( e 3 ex)( e 1 ex) ( e 3 y)( 1 y) ( 3 z)( 1 z) ( 3 x)( 1 y) ( 3 y)( 1 x) ( 3 y)( 1 z) ( 3 z)( 1 y) ( 3 z)( 1 x) ( 3 x)( 1 z) σ e e e e e e e e e e e + e e e e e e e e + e e e e e e e e + e e e e zx l1 m1 n1 l1m1 m1n1 nl σ 1 1 xx σ l m n lm mn nl σyy l3 m3 n3 l3m3 m3n3 n3l σ 3 zz = ll 1 mm 1 n1n l1m + ml 1 mn 1 + n1m nl 1 + l1n σ xy ll3 mm3 nn3 lm3 + ml3 mn3 + nm3 nl3+ ln3 σ (4.8) yz l3l1 m3m1 n3n1 l3m1 + m3l1 m3n1 + n3m1 n3l1 + l3n 1 σ zx R σ σ li = ( e i ex) mi = ( e i ey) ni = ( e i ez) (4.9) 30

D マトリックスの変換 R ε R l1 m1 n1 l1m1 m1n1 nl 1 1 l m n lm mn nl l3 m3 n3 l3m3 m3n3 n3l3 = l1l mm 1 n1n l1m + ml 1 m1n + n1m nl 1 + l1n ll3 mm3 nn3 lm3+ ml3 mn3 + nm3 nl3+ ln3 l3l1 m3m1 n3n1 l3m1 + m3l1 m3n1 + n3m1 n3l1 + l3n 1 l1 l l3 ll 1 ll3 l3l 1 m1 m m3 mm 1 mm3 m3m1 n n n n n n n n n = l1m1 lm l3m3 l1m + ml 1 lm3+ ml3 l3m1 + m3l 1 m1n1 mn m3n3 m1n + n1m mn3 + nm3 m3n1 + n3m1 nl 1 1 nl n3l3 nl 1 + l1n nl3 + ln3 n3l1 + l3n 1 T 1 3 1 3 3 1 ε (4.10) R σ l1 m1 n1 l1m1 m1n1 nl 1 1 l m n lm mn nl l3 m3 n3 l3m3 m3n3 n3l3 = l1l m1m n1n l1m + ml 1 m1n + n1m nl 1 + l1n ll3 mm3 nn3 lm3+ ml3 mn3 + nm3 nl3+ ln3 l3l1 m3m1 n3n1 l3m1 + m3l1 m3n1 + n3m1 n3l1 + l3n 1 (4.11) (4.11) R = R 1 T σ ε (4.1) σ = D ε (4.13) R σ = D R ε σ R R 1 1 σ σ σ ε σ = = 1 σ T ε σ = ε R D R R D R T ε R D R ε ε ε D= R D R ε ε (4.14) ε (4.15) 剛性マトリックス ( 非適合要素の場合 ) の求め方は, 第 3 回 FrontISTR 研究会資料 静弾性解析 の1.3と同じです [1] Bathe, K.J., Finite Element Procedures, Prentice Hall, 1996. 31

目次 前半 解析機能/ サンプル例題 1. 解析機能とユーザマニュアル該当箇所. メッシュファイルと解析制御ファイルの設定方法 注意点 3. サンプル例題 ( 内圧を受ける血管を模擬した円管モデル ) 後半 定式化/ プログラム 4. 直交異方弾性体の有限要素法定式化 5. プログラム説明 3

FrontISTR_V43_p1.tar.gz を解凍します ディレクトリ src の下がソースファイル群です FrontISTR Ver.3.5 のメインプログラムです 四つのディレクトリ main, common, analysis, lib があります 33

データの読み込み関係のプログラム 動解析用プログラム伝熱解析用プログラム 静解析用プログラム 有限要素の幾何情報を計算するプログラム B マトリックスの計算で使用 材料情報を計算するプログラム D マトリックスの計算で使用 34

[main/fistr_main.f90] PROGRAM fstr_main メインプログラム hecmw_init() hecmw_get_mesh() [main/fistr_main.f90]fstr_init() 変数初期化 入力データ読み込み hecmw_nullify_matrix() hecmw_nullify_result_data() [main/fistr_main.f90] fstr_init_file() hecmw_mat_con() [main/fistr_main.f90] fstr_condition() hecmw_ctrl_get_control_file() [main/fistr_main.f90] fstr_linear_static_analysis() 線形静解析用のルーチンへ [analysis/static/fstr_solve_linear.f90] m_fstr_linear::fstr_solve_linear() [analysis/static/static_mat_ass.f90] m_static_mat_ass::fstr_mat_ass() 全体剛性マトリックスの作成 [analysis/static/static_mat_ass_main.f90] m_static_mat_ass_main::fstr_mat_ass_main() hecmw_mat_clear() [analysis/static/static_mat_ass_main.f90] m_static_mat_ass_main::fstr_local_stf_create() 要素剛性マトリックスの計算 [analysis/static/static_lib_3dic.f90] m_static_lib_3dic::stf_c3d8ic() 3 次元六面体 1 次要素 [lib/element/element.f90] elementinfo::getquadpoint() Gaussの積分点数 [lib/element/element.f90] elementinfo::getglobalderiv() 形状関数の微分値 [lib/physics/calmatmatrix.f90] m_matmatrix::matlmatrix() Dマトリックス [lib/physics/elasticlinear.f90] m_elasticlinear::calelasticmatrix_ortho() 直交異方弾性体の場合 hecmw_mat_ass_elem() 要素剛性マトリックスをassemble [analysis/static/fstr_ass_load.f90] m_fstr_ass_load::fstr_ass_load() 外力ベクトルの計算 [analysis/static/fstr_addbc.f90] m_fstr_addbc::fstr_addbc() 境界条件の処理 hecmw_allreduce_r1() [lib/solve_lineq.f90] m_solve_lineq::solve_lineq() 線形ソルバーによる求解 hecmw_solve_33() hecmw_update_3_r() [analysis/static/fstr_update.f90] m_fstr_update::fstr_update3d() [lib/static_lib_3dic.f90] m_static_lib_3dic::updatest_c3d8ic() 応力の計算 [lib/static_lib_3dic.f90] m_static_lib_3dic::stf_c3d8ic() 要素剛性マトリックスの計算(3 次元六面体 1 次要素の場合 ) [analysis/static/static_output.f90] m_static_output:: fstr_static_output() 結果の出力 [analysis/static/static_make_result.f90] m_static_make_result::fstr_write_static_result() [main/fistr_main.f90] fstr_main::fstr_finalize() 変数の削除 hecmw_finalize() 35 [ ディレクトリ / ファイル名 ] モジュール名 :: サブルーチン名 () を意味しています

[main/fistr_main.f90] PROGRAM fstr_main メインプログラム hecmw_init() hecmw_get_mesh() [main/fistr_main.f90] fstr_init() 変数初期化 入力データ読み込み hecmw_nullify_matrix () hecmw_nullify_result_data () [main/fistr_main.f90] fstr_init_file() hecmw_mat_con () [main/fistr_main.f90] fstr_condition() hecmw_ctrl_get_control_file () [main/fistr_main.f90] fstr_nonlinear_static_analysis() 非線形静解析用のルーチンへ [analysis/static/fstr_solve_nlgeom.f90] m_fstr_solve_nlgeom::fstr_solve_nlgeom() 荷重増分のループ [analysis/static/fstr_solve_nonlinear.f90] m_fstr_nonlinearmethod::fstr_newton(). loading step (substep) 1. [analysis/static/fstr_solve_nonlinear.f90] m_fstr_nonlinearmethod::fstr_newton(). loading step (substep)... [analysis/static/static_output.f90] m_static_output::fstr_static_output() 計算結果の出力 [analysis/static/static_make_result.f90] m_static_make_result::fstr_write_static_result() [main/fistr_main.f90] fstr_finalize() 変数の削除 hecmw_finalize () [ ディレクトリ / ファイル名 ] モジュール名 :: サブルーチン名 () を意味しています 36

[analysis/static/fstr_solve_nonlinear.f90] m_fstr_nonlinearmethod::fstr_newton() Newton-Raphson 反復 hecmw_allreduce_i1() [analysis/static/fstr_ass_load.f90] m_fstr_ass_load::fstr_ass_load() 外力ベクトルの計算 [analysis/static/fstr_stfiffmatrix.f90] m_fstr_stiffmatrix::fstr_stiffmatrix() 全体接線剛性マトリックスの計算 hecmw_mat_clear() [lib/static_lib_c3d8.f90] m_static_lib_c3d8::stf_c3d8bbar() 3 次元六面体 1 次要素 [lib/element/element.f90] elementinfo::getquadpoint() Gaussの積分点数 [lib/element/element.f90] elementinfo::getglobalderiv() 形状関数の微分値 [lib/physics/calmatmatrix.f90] m_matmatrix::matlmatrix() Dマトリックス [lib/physics/elasticlinear.f90] m_elasticlinear::calelasticmatrix_ortho() 直交異方弾性体の場合 hecmw_mat_ass_elem() 要素接線剛性マトリックスをassemble [analysis/static/fstr_addbc.f90] m_fstr_addbc::fstr_addbc() 境界条件の処理 [lib/solve_lineq.f90] m_solve_lineq::solve_lineq() 線形ソルバーによる求解 Newton-Raphson iteration (iter) 1 hecmw_update_3_r() [analysis/static/fstr_update.f90] m_fstr_update::fstr_updatenewton() 応力 内力ベクトルの計算 [lib/static_lib_c3d8.f90] m_static_lib_c3d8::update_c3d8bbar() 3 次元六面体 1 次要素 [lib/element/element.f90] elementinfo::getquadpoint() Gaussの積分点数 [lib/element/element.f90] elementinfo::getglobalderiv() 形状関数の微分値 [lib/physics/calmatmatrix.f90] m_matmatrix::matlmatrix() Dマトリックス [lib/physics/elasticlinear.f90] m_elasticlinear::calelasticmatrix_ortho() 直交異方弾性体の場合 hecmw_update_3_r() [analysis/static/fstr_residual.f90] m_fstr_residual::fstr_update_ndforce() 残差ベクトルの計算 hecmw_allreduce_r1() [analysis/static/fstr_stfiffmatrix.f90] m_fstr_stiffmatrix::fstr_stiffmatrix().. [analysis/static/fstr_residual.f90] m_fstr_residual::fstr_update_ndforce() hecmw_allreduce_r1() [analysis/static/fstr_stfiffmatrix.f90] m_fstr_stiffmatrix::fstr_stiffmatrix().. [ ディレクトリ / ファイル名 ] モジュール名 :: サブルーチン名 () を意味しています Newton-Raphson iteration (iter) 37

D マトリックスの計算で使用するサブルーチン [lib/utilities/utilities.f90] m_utilities::transformation() [lib/physics/elasticlinear.f90] m_elasticlinear::calelasticmatrix_ortho() [lib/physics/calmatmatrix.f90] m_matmatrix::matlmatrix() B マトリックスの計算で使用するサブルーチン [lib/element/element.f90] elementinfo::getglobalderiv() 剛性マトリックスの計算で使用するサブルーチン [analysis/static/static_lib_3dic.f90] m_static_lib_3dic::stf_c3d8ic() 全体剛性マトリックスを計算するサブルーチン [analysis/static/static_mat_ass_main.f90] m_static_mat_ass_main::fstr_local_stf_create() 38

モジュール名 :m_utilities 補助的なサブルーチンや関数を集めたモジュール 使用する他のモジュール hecmw HECMW のモジュール メンバ変数 整数型 kreal 実数型の種別値 実数型 PI=3.1415965358979 円周率 メンバ関数 サブルーチン memget() 使用されるメモリを記録するサブルーチン サブルーチン append_intname() ファイル名の後ろに整数を入れるサブルーチン サブルーチン insert_intarray() 整数を整数配列へ入れるサブルーチン サブルーチン tensor_eigen3() 3 3 の対称マトリックスの固有値を計算するサブルーチン サブルーチン eigen3() 3 3 の対称マトリックスの固有値と固有ベクトルを計算するサブルーチン 関数 Determinant() 3 3 のマトリックスの determinant を計算するサブルーチン サブルーチン fstr_chk_alloc() メモリオーバーフローをチェックするサブルーチン サブルーチン calinverse() マトリックスの逆を計算するサブルーチン サブルーチン cross_product() ベクトルの外積を計算するサブルーチン サブルーチン transformation() 座標変換マトリックスを作成するサブルーチン 39

引数 実数型 jacob(3, 3) 単位基底ベクトルの方向余弦 実数型 tm(6, 6) 座標変換マトリックスの成分 上位なし 下位なし サブルーチン名 :transformation() 座標変換マトリックスを作成するサブルーチン li = ( e i ex) = jacob( i,1) mi = ( e i ey) = jacob( i, ) ni = ( e i ez) = jacob( i, 3) R ε tm(1,1) tm(1, ) tm(1, 3) tm(1, 4) tm(1, 5) tm(1, 6) tm(,1) tm(, ) tm(, 3) tm(, 4) tm(, 5) tm(, 6) tm(3,1) tm(3, ) tm(3, 3) tm(3, 4) tm(3, 5) tm(3, 6) = tm(4,1) tm(4, ) tm(4, 3) tm(4, 4) tm(4, 5) tm(4, 6) tm(5,1) tm(5, ) tm(5, 3) tm(5, 4) tm(5, 5) tm(5, 6) tm(6,1) tm(6, ) tm(6, 3) tm(6, 4) tm(6, 5) tm(6, 6) 40

使用する他のモジュール [lib/physics/material.f90] mmaterial 材料物性の情報を管理するモジュール メンバ変数 整数型 kreal 実数型の種別値 モジュール名 :m_elasticlinear 線形弾性体の D マトリックスを計算するモジュール メンバ関数 サブルーチン calelasticmatrix() 3 次元問題, 平面ひずみ問題, 平面応力問題, 軸対称問題の D マトリックスを計算するサブルーチン サブルーチン calelasticmatrix_ortho() 直交異方性がある場合,3 次元問題の D マトリックスを計算するサブルーチン サブルーチン LinearElastic_Shell() シェル要素を使用する場合, 埋め込み座標系成分の D マトリックスを計算するサブルーチン 41

サブルーチン名 :calelasticmatrix_ortho() 直交異方性がある場合,3 次元問題の D マトリックスを計算するサブルーチン 引数 構造体 (tmaterial) matl 材料に関連するデータ 整数型 secttype 問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) 実数型 bij(3, 3) li = ( e i ex) = bij( i,1) mi = ( e i ey) = bij( i, ) ni = ( e i ez) = bij( i, 3) 実数型 DMAT(:, :) D マトリックスの成分 実数型 temp ( 省略可能 ) 温度 上位 サブルーチン [lib/physics/calmatmatrix.f90] m_matmatrix::matlmatrix() 下位 サブルーチン [lib/utilities/ttalbe.f90] m_table::fetch_tabledata() サブルーチン [lib/utilities/utilities.f90] m_utilities::transformation() 4

モジュール名 :m_matmatrix 各材料の D マトリックスを計算するサブルーチンを呼ぶモジュール 使用する他のモジュール [lib/physics/material.f90] mmaterial 材料物性の情報を管理するモジュール [lib/physics/mechgauss.f90] mmechgauss Gauss 積分点の情報を管理するモジュール [lib/physics/elasticlinear.f90] m_elasticlinear 線形弾性体の D マトリックスを計算するモジュール [lib/physics/hyperelastic.f90] mhyperelastic 超弾性体の 4 階の弾性テンソルを計算するモジュール [lib/physics/elastoplastic.f90] m_elastoplastic 弾塑性体の D マトリックスを計算するモジュール [lib/physics/viscoelastic.f90] mviscoelastic 粘弾性体の D マトリックスを計算するモジュール [lib/physics/creep.f90] mcreep クリープを考慮した剛性マトリックスを計算するためのモジュール [lib/user/uelastic.f90] muelastic ユーザ定義の弾性体の D マトリックスを計算するモジュール [lib/user/umat.f90] mumat ユーザ定義の材料の D マトリックスを計算するモジュール メンバ変数 整数型 kreal 実数型の種別値 メンバ関数 サブルーチン getnlgeomflag() 未使用のサブルーチン サブルーチン MatlMatrix() 各材料の D マトリックスを計算するサブルーチンを呼ぶサブルーチン サブルーチン StressUpdate() 各材料の応力とひずみを計算するサブルーチンを呼ぶサブルーチン サブルーチン mat_cd() 材料が超弾性体の場合,4 階の弾性テンソルを問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) に応じた D マトリックスに変換するサブルーチン サブルーチン MatlMatrix_Shell() シェル要素を使用する場合, 各材料 ( 現バージョンでは, 線形弾性体のみ ) の応力とひずみを計算するサブルーチンを呼ぶサブルーチン サブルーチン mat_cd_shell() シェル要素を使用する場合,4 階の弾性テンソルを D マトリックスに変換するサブルーチン 43

サブルーチン名 :MatlMatrix() 各材料の D マトリックスを計算するサブルーチンを呼ぶサブルーチン 引数 構造体 (tgaussstatus) gauss Gauss 積分点に関連するデータ 整数型 secttype 問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) 実数型 matrix(:, :) D マトリックスの成分 実数型 dt 時間増分 実数型 cdsys(3, 3) 直交異方性がある場合に使用する座標系 実数型 temperature 省略可能温度 上位 サブルーチン [lib/static_lib_d.f90] m_static_lib_d::stf_c() サブルーチン [lib/static_lib_d.f90] m_static_lib_d::update_c() サブルーチン [lib/static_lib_d.f90] m_static_lib_d::updatest_c() サブルーチン [lib/static_lib_3d.f90] m_static_lib_3d::stf_c3() サブルーチン [lib/static_lib_3d.f90] m_static_lib_3d::tload_c3 () サブルーチン [lib/static_lib_3d.f90] m_static_lib_3d::update_c3() サブルーチン [lib/static_lib_3d.f90] m_static_lib_3d::updatest_c3() サブルーチン [lib/static_lib_3dic.f90] m_static_lib_3dic::stf_c3d8ic() サブルーチン [lib/static_lib_3dic.f90] m_static_lib_3dic::updatest_c3d8ic() サブルーチン [lib/static_lib_3dc3d8.f90] m_static_lib_c3d8::stf_c3d8bbar() サブルーチン [lib/static_lib_3dc3d8.f90] m_static_lib_c3d8::update_c3d8bbar() サブルーチン [lib/static_lib_3dc3d8.f90] m_static_lib_c3d8::tload_c3d8bbar() 下位 サブルーチン [lib/physics/calmatmatrix.f90] m_matmatrix::mat_cd() サブルーチン [lib/user/uelastic.f90] muelastic::uelasticmatrix() サブルーチン [lib/physics/viscoelastic.f90] mviscoelastic::calviscoelasticmatrix() サブルーチン [lib/physics/elasticlinear.f90] m_elasticlinear::calelasticmatrix() サブルーチン [lib/physics/elasticlinear.f90] m_elasticlinear::calelasticmatrix_ortho() サブルーチン [lib/physics/hyperelastic.f90] mhyperelastic::calelasticmooneyrivlin() サブルーチン [lib/physics/hyperelastic.f90] mhyperelastic::calelasticarrudaboyce() サブルーチン [lib/physics/elastoplastoc.f90] m_elastoplastic::calelastoplasticmatrix() サブルーチン [lib/user/umat.f90] mumat::umatlmatrix() サブルーチン [lib/user/creep.f90] mcreep::iso_creep() 44

モジュール名 :m_static_lib_3dic 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,B マトリックスおよび要素剛性マトリックスを計算したり,Gauss 積分点における応力とひずみを計算したりするモジュール 使用する他のモジュール hecmw HECMW のモジュール [lib/utilities/utilities.f90] m_utilities 補助的なサブルーチンや関数を集めたモジュール [lib/element/element.f90]elementinfo 要素の情報を管理するモジュール [lib/physics/mechgauss.f90] mmechgauss Gauss 積分点の情報を管理するモジュール [lib/m_common_struct.f90] m_common_struct 有限要素解析における共通データを定義するモジュール [lib/physics/calmatmatrix.f90] m_matmatrix 各材料の D マトリックスを計算するサブルーチンを呼ぶモジュール [lib/m_fstr.f90] m_fstr FrontISTR における共通データを定義するモジュール メンバ変数 整数型 kint 整数型の種別値 実数型 kreal 実数型の種別値 メンバ関数 サブルーチン STF_C3D8IC() 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,B マトリックスおよび要素剛性マトリックスを計算するサブルーチン サブルーチン UpdateST_C3D8IC() 3 次元六面体 8 節点要素 ( 非適合要素 ) の場合,Gauss 積分点における応力とひずみを計算するサブルーチン 45

サブルーチン名 :STF_C3D8IC() 3 次元六面体 8 節点要素 (B-bar 要素 ) の場合,B マトリックスおよび要素剛性マトリックスを計算するサブルーチン 引数 整数型 etype 要素タイプ 整数型 nn 各要素の節点数 (nn=8) 実数型 ecoord(3, nn) 各要素の節点座標 構造体 (tgaussstatus) gausses(:) Gauss の積分点に関連するデータ 実数型 stiff(:, :) 要素剛性マトリックス 実数型 tincr 時間増分 実数型 coords(3, 3) 材料の局所座標系を定義するのに必要な変数 実数型 mddisp(3, :) 省略可能節点変位 実数型 ehdisp(3, 3) 省略可能曲げモードのための節点自由度 上位 サブルーチン [analysis/static/fstr_stiffmatrix.f90] m_fstr_stiffmatrix::fstr_stiffmatrix() 下位 サブルーチン [lib/element/element.f90] elementinfo::getjacobian() サブルーチン [lib/m_common_struct.f90] m_common_struct::set_localcoordsys() サブルーチン [lib/element/element.f90] elementinfo::getshapefunc() サブルーチン [lib/physics/calmatmatrix.f90] m_matmatrix::matlmatrix() サブルーチン [lib/element/element.f90] elementinfo::getquadpoint() サブルーチン [lib/element/element.f90] elementinfo::getglobalderiv() 関数 [lib/element/element.f90] elementinfo::getweight() サブルーチン [lib/utilities/utilities.f90] m_utilities::calinverse() サブルーチン [lib/element/element.f90] elementinfo::getshapefunc() サブルーチン [lib/utilities/ttable.f90] m_table::fetch_tabledata() 46

モジュール名 :m_static_mat_ass_main 全体剛性マトリックスを計算するモジュール 使用する他のモジュール [lib/m_fstr.f90] m_fstr FrontISTR における共通データを定義するモジュール [lib/static_lib.f90] m_static_lib 静解析で使用するモジュールを定義するモジュール [lib/physics/mechgauss.f90] mmechgauss Gauss 積分点の情報を管理するモジュール メンバ関数 サブルーチン FSTR_MAT_ASS_MAIN() 全体剛性マトリックスを作成するサブルーチン サブルーチン FSTR_LOCAL_STF_CREATE() 要素剛性マトリックスを計算するサブルーチンを呼ぶサブルーチン 47

サブルーチン名 :FSTR_LOCAL_STF_CREATE() 要素剛性マトリックスを計算するサブルーチンを呼ぶサブルーチン 引数 構造体 (hecmwst_matrix) hecmesh HECMW が管理するメッシュのデータ 整数型 ndof 各節点の自由度数 整数型 ic_type 要素タイプ 整数型 icel 要素番号 実数型 xx(:), yy(:), zz(:) 各要素の節点座標 構造体 (tgaussstatus) gausses(:) Gauss の積分点に関連するデータ 整数型 iset 問題の種類 (3 次元問題 / 平面ひずみ / 平面応力問題 / 軸対称問題 ) 実数型 stiffness(:, :) 要素剛性マトリックス 上位 サブルーチン m_static_mat_ass_main:: FSTR_MAT_ASS_MAIN() 下位 サブルーチン m_static_lib_d :: STF_C() サブルーチン m_static_lib_1d :: STF_C1() サブルーチン m_static_lib_3dic :: STF_C3D8IC() サブルーチン m_static_lib_3d :: STF_C3() サブルーチン m_static_lib_shell:: STF_Shell_MITC() サブルーチン m_static_lib_beam:: STF_Beam() サブルーチン hecmw_abort() 48

次回の < 機能 例題 定式化 プログラム解説編 > では, 弾塑性解析 や 熱応力解析 を行う予定です 次回は有限変形理論を扱いますプログラム解説では, 第 10 回 FrontISTR 研究会で説明できなかった荷重増分ステップや Newton-Raphson 反復に対応するサブルーチンも説明します 定式化の理解を深められる例題を増やしていく予定です 今年度の最終目標は, FrontISTR のカスタマイズ (Element/Material 追加およびユーザサブルーチン使用 ) の解説です 49

変位 u によって材料で定義された座標が変化する場合の計算方法については? 主軸が設定できると便利である LS-DYNA や Abaqus では要素番号と座標系のテーブルを入力することが可能である 繊維方向の応力成分を出力したい 設定した座標系を Revocap_PrePost で見られると良いのだが 50