2016年2月2日
第25回 FrontISTR研究会
＜事例・サービスの紹介／
FrontISTRの構造要素 (シェル要素／梁要素)の解説＞
FrontISTRの梁要素／シェル要素の解説
東京大学
新領域創成科学研究科
人間環境学専攻
橋本 学
『
FrontISTRに実装されている定式化を十分に理解し，
解きたい問題に対して
ソースコードを自由にカスタマイズ
(要素タイプを追加，材料の種類を追加，ユーザサブルーチン
を追加
)
できるようになる
こと』
を最終目標とします
今回は，
FrontISTRに実装されている構造要素
(Bernoulli-Euler梁要素／MITCシェル要素) を解説します
また，シェル要素／梁要素とソリッド要素を混ぜた
解析メッシュを使用する場合の計算方法について解説します
2FrontISTRのメッシュファイル
ヘッダー
節点の番号
I と
節点の座標値
(x
I, y
I, z
I)
要素の番号
e と
要素内の
節点同士の
つながり
(コネクティビティ)
節点グループ
セクションデータ
材料データ
セクションの設定によって
要素グループに材料データが
与えられる
境界条件や集中荷重(外力)を与える
ファイルの終わり
梁要素: TYPE=BEAM
シェル要素: TYPE=SHELL
梁要素：TYPE=611
シェル要素：
TYPE=741
TYPE=743
3(注意) 解析メッシュにソリッド要素，梁要素，シェル要素が混在する場合，
節点自由度数を
3に揃えるため，*ではなく**の要素タイプ番号を使用してください
4要素種類
要素タイプ番号
節点数
節点自由度数
説明
線要素
111
2
3
2節点リンク要素
112
3
3
3節点リンク要素
平面要素
231
3
3
三角形1次要素
232
6
3
三角形2次要素
241
4
3
四角形1次要素
242
8
3
四角形
2次要素 (Serendipity族)
ソリッド要素
301
2
3
2節点トラス要素
341
4
3
四面体1次要素
342
10
3
四面体2次要素
351
6
3
プリズム1次要素
352
15
3
プリズム
2次要素
361
8
3
六面体
1次要素
362
20
3
六面体2次要素 (Serendipity族)
インターフェース要素
541
4×2
3
四角形面1次要素
542
8×2
3
四角形面2次要素
梁要素
611
*
2
6
2節点梁要素 (Bernoulli-Euler梁)
641
**
2×2
3
2節点梁要素 (Bernoulli-Euler梁)
シェル要素
731
*
3
6
三角形1次要素 (MITC3シェル)
761
**
3×2
3
三角形1次要素 (MITC3シェル)
741
*
4
6
四角形1次要素 (MITC4シェル)
781
**
4×2
3
四角形
1次要素 (MITC4シェル)
743
9
6
四角形
2次要素 (MITC9シェル)
5 611
3次元ソリッド要素
3○○
シェル要素
7○○
梁要素
6○○
要素形状
要素形状
要素形状
6
静解析
動解析
荷重増分解析不要
荷重増分解析要
直接時間積分
固有値解析
周波数応答解析
微小変形理論
(線形弾性体のみ)
微小変形理論
有限変形理論
微小変形理論
(線形弾性体のみ)
有限変形理論
!SOLUTION,TYPE=DYNAMIC
!DYNAMIC下のidx_respを1に設定
!SOLUTION,TYPE=EIGEN
!SOLUTION,TYPE=DYNAMIC
!DYNAMIC下のidx_respを2に設定
!SOLUTION,TYPE=NLSTATIC
!SOLUTION,TYPE=STATIC
!DYNAMIC,TYPE=LINEAR
!DYNAMIC,TYPE=NONLINEAR
材料定義において
INIFINITEを追記
例：!ELASTIC, INFINITE
構造解析
熱伝導解析
定常解析
非定常解析
!SOLUTION,TYPE=HEAT
!HEATの下に記述無または0.0と記述
ひずみの
2次以上の項がある
微小変形理論
(微小変位)
微小ひずみ
線形弾性体
弾塑性体
粘弾性体
クリープ材
有限変形理論
(有限変位)
微小ひずみ
線形弾性体
弾塑性体
粘弾性体
クリープ材
大ひずみ
弾塑性体
超弾性体
粘弾性体
クリープ材
7
0 0 T 0 0 T
01
(
) (
)
(
) (
)
2
t
t
t
t
tE
u
u
u
u
0tS
f
(
0tE
,
0tE
0tE
, ...)
変位こう配の
2次項がある
有限変形理論
大ひずみ
ひずみ
応力
!ELASTIC
!ELASTICと!PLASTIC
!HYPERELASTIC
!ELASTICと!VISCOELASTIC
!ELASTIC
!ELASTICと!PLASTIC
!ELASTICと!CREEP
!ELASTICと!VISCOELASTIC
!ELASTICと!CREEP
青字はFrontISTRで解析可能な材料
8
x
u
t
b
有界領域
[m
N]
境界
[m
N - 1]
物質点の位置ベクトル
[m]
変位
[m]
トラクション
[Pa]
t
時刻
[s]
ナブラ
[1/m]
単位質量当たりの体積力
[N/kg]
次元
(3次元：N = 3)
N
密度
[kg/m
3]
0 0 T 0 0 T
01
(
) (
)
(
) (
)
2
t
t
t
t
tE
u
u
u
u
時刻
t の物理量
基準となる時刻が時刻
0の意味
, ,
A B a
b
,
a b
,
a b
スカラー
ベクトル
2階のテンソル
:
A B
ij ijA B =
i ia b
a b =
n
外向き単位法線ベクトル
[-]
9
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
２．シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
10
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・
理論
・ 現バージョンのプログラム
・ 解析事例
２．シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
(a) 時刻 t = 0に梁の中立軸に垂直な直線は梁が変形する間
も直線であるが，変形した中立軸に垂直である
(せん断変形は考慮されない)
Bernoulli-Eulerの仮定
(b) 梁の断面は変形しない
微小ひずみの仮定
M
M
t = 0
11(2) (2) (2) (2) (2) (2) x y z x y z u u u
(1) (1) (1) (1) (1) (1) x y z x y zu
u
u
(1)
(2)
Bernoulli-Euler梁要素
1ˆe
ref 1 2 ref 1ˆ
ˆ
ˆ |
e
e
e
| e
e
3 1 2ˆ
ˆ
ˆ
e
e
e
x e y e z e x y z 1 ˆx 3 ˆx 2 ˆx鷲津久一郎
, 宮本博, 山田嘉昭, 山本善之, 川井忠彦共編, ''有限要素法ハンドブック I 基礎編'' (1981), p.218-220.
中立軸方向の一様引張・一様圧縮
曲げ
ねじり
方向変位：
の
1次式を仮定
1ˆx
ˆx
1方向変位：
の
3次式を仮定
2ˆx
ˆx
1方向変位：
の
3次式を仮定
3ˆx
ˆx
1方向回転：
の
1次式を仮定
1ˆx
ˆx
1微小変形を仮定
2 3 1ˆ
ˆ
ˆ
u
x
3 2 1ˆ
ˆ
ˆ
u
x
方向回転：
3ˆx
方向回転：
2ˆx
並進自由度と
回転自由度
12(2) (2) (2) (2) (2) (2) x y z x y z u u u
(1) (1) (1) (1) (1) (1) x y z x y zu
u
u
(1)
(2)
Bernoulli-Euler梁要素
1ˆe
ref 1 2 ref 1ˆ
ˆ
ˆ |
e
e
e
| e
e
3 1 2ˆ
ˆ
ˆ
e
e
e
x e y e z e x y z 1 ˆx 3 ˆx 2 ˆx 1 1 1 1 2 2 2 2 3 3 3 3ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
x y z x x y z y x y z zu
u
u
u
u
u
e e
e e
e e
e e
e e
e e
e e
e e
e e
鷲津久一郎
, 宮本博, 山田嘉昭, 山本善之, 川井忠彦共編, ''有限要素法ハンドブック I 基礎編'' (1981), p.218-220.
eT
局所座標成分
全体座標成分
ˆ
ˆ
ˆ
ˆ
ˆ
(
ˆ
)
(
ˆ
)
ij i j ij i j e e i j i k kl l jT
ik klT
jl
e
e
e
e
e
e
e e
e e
ˆ
ˆ
ˆ
(
ˆ
)
i i i i e i i j j ij ju
u
u
T u
u
e
e
e u
e e
13変位ベクトル
応力テンソル
並進自由度と
回転自由度
1 1 1 1 2 2 2 2 3 3 3 3ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
ˆ
x y z x x y z y x y z z
e e
e e
e e
e e
e e
e e
e e
e e
e e
eT
局所座標成分
全体座標成分
2節点Bernoulli-Euler梁要素 (3/4)
14 3 3 3 3 3 2 3 2 2 2 2 2 3 2 3 2 1 1 2 2 2 2 2 2 3 3 3 3 2 2 3 3 0 0 0 0 0 0 0 0 0 0 12 6 12 6 0 0 0 0 0 0 0 0 12 6 12 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 4 6 2 0 0 0 0 0 0 0 0 6 4 6 2 0 0 0 0 0 0 0 0 ˆ 0 0 0 0 0 0 0 0 0 0 12 0 e EA EA l l EI EI EI EI l l l l EI EI EI EI l l l l GJ GJ l l EI EI EI EI l l l l EI EI EI EI l l l l EA EA l l EI l K 3 3 3 2 3 2 2 2 2 2 3 2 3 2 1 1 2 2 2 2 2 3 3 3 3 6 12 6 0 0 0 0 0 0 0 12 6 12 6 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 2 6 4 0 0 0 0 0 0 0 0 6 2 6 4 0 0 0 0 0 0 0 0 EI EI EI l l l EI EI EI EI l l l l GJ GJ l l EI EI EI EI l l l l EI EI EI EI l l l l (1) 1 (1) 2 (1) 3 (1) 1 (1) 2 (1) 3 (2) 1 (2) 2 (2) 3 (2) 1 (2) 2 (2) 3 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ e u u u u u u
uˆ
ˆ ˆ
e e
eK u
f
局所座標系
: Young率
: 梁の長さ
：断面
2次モーメント
: ねじり定数
E
l
I
2,
I
3J
2節点Bernoulli-Euler梁要素 (4/4)
15 T Tˆ
ˆ
e e e e e e e e e e e e e e e
T
O
O O
T
O
O O
T
O O O
O
T
O O
O
T
O O
O
T
O O
K
u
f
O
O
T
O
O
O
T
O
O
O T O
O
O
O T
O
O
O T
O
O O T
eK
f
e全体座標系の要素剛性マトリックス
を作成後，
を全体剛性マトリックスに加える
(1) (1) (1) (1) (1) (1) (2) (2) (2) (2) (2) (2) x y z x y z e x y z x y z u u u u u u
u eK
eK
(a) 時刻 t = 0にシェルの中立面に垂直な直線はシェルが変形する
間も直線であるが，変形した中立面に垂直である必要はない
(せん断変形は考慮される)
Reissner-Mindlin板
の仮定
(b) シェルのディレクタは変形しない
微小ひずみの仮定
M
M
t = 0
161 g 2 g 3 g 1 3 1 2 (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (3) 1 V (3) 2 V (3) 3 V (4) 1 V (4) 2 V (4) 3 V 3 4 1 e 2 e 3 e 1 x 2 x 3 x 2
FrontISTRでは微小変形を仮定したMITCシェル要素が実装されている
1 g 2 g 3 g 1 3 3 1 (3) 1 V (3) 2 V (3) 3 V (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V 2 1 e 2 e 3 e 1 x 2 x 3 x 2 MITC4シェル要素
1 2 (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (3) 1 V (3) 2 V (3) 3 V (4) 1 V (4) 2 V (4) 3 V 3 4 1 e 2 e 3 e 1 x 2 x 3 x 1 g 2 g 3 g 1 3 2 5 6 8 7 9 (6) 1 V (6) 2 V (6) 3 V (5) 1 V (5) 2 V (5) 3 V (7) 1 V (7) 2 V (7) 3 V (8) 1 V (8) 2 V (8) 3 V (5) 1 V (5) 2 V (5) 3 VMITC9シェル要素
MITC3シェル要素
Bathe, K.J., Finite Element
Procedures, Prentice Hall, (1995).
(a) Isoparametric
degenerated
shell element
(b) The
covariant components
measured in the convected
system are used as the Green-Lagrange strain
components.
(c) The transverse shear strain components are evaluated
using
the values interpolated at sampling points
.
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 2 x x y y z z
u
u
u
u
e
e
e
V
V
Fig. MITC4 shell element
(α) (α) (α) 1 2 3 V , V , V
: 節点での単位直交ベクトル
: ディレクタベクトル
a: 厚さ
位置ベクトル
変位ベクトル
(微小変形を仮定)
1 2 3= x x x g , g , g: 共変基底ベクトル
節点あたり
5自由度
(α) 3 V 0 0 (α) 3 3 0 3| g V | g 1 g 2 g 3g
1 2 (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (3) 1 V (3) 2 V (3) 3 V (4) 1 V (4) 2 V (4) 3 V 3 4 x e y e z e x y z 4 ( ) ( ) ( ) 3 1( , )
2
a
N
x
x
V
4 ( ) ( ) ( ) 0 ( ) 3 3 1 4 ( ) ( ) ( ) 0 ( ) 3 1( , )
(
)
2
( , )
(
)
2
a
N
a
N
u
u
V
V
u
V
0 4 ( ) ( ) ( ) 0 ( ) 3 3 1( , )
(
)
2
a
N
u
x
x
u
V
V
Dvorkin, E.N. and Bathe, K.J., “A Continuum Mechanics Based Four-node Shell Element for General Non-linear Analysis,”
Engineering Computations, Vol.1, pp.77-88, (1984).
AS AS AS 31 31 31, A 31,C AS AS AS 23 23 23, D 23, B
1
1
(1
)
(1
)
2
2
1
1
(1
)
(1
)
2
2
* 31 DI 23 DI 31 31 23 23(
)
(
)
V
dV
V
dV
31 A C 23 D B( ) (1
)
( ) (1
)
(1
) ( )
(1
) ( )
1 g 2 g 1 2 4 3 1 g 2 g 1 2 4 C 3 A D BDvorkin, E.N. and Bathe, K.J., “A Continuum Mechanics Based Four-node Shell Element for General Non-linear Analysis,”
Engineering Computations, Vol.1, pp.77-88, (1984).
0 0 0 0 ij i j ij
i jg
g
g
g
Tensorial components
:
C
t1
:
2
VdV
SdS
V
dV
u t
u
b
20t
:
V
dV
S
dS
V
dV
u t
u
b
*0
より，仮想仕事の原理
AS DI DI 31 31 31 31 AS DI DI 23 23 23 231
1
(1
)
(0, 1, 0)
(1
)
(0, 1, 0)
2
2
1
1
(1
)
( 1, 0, 0)
(1
)
(1, 0, 0)
2
2
Mixed Interpolation of Tensorial Components
Dvorkin, E.N. and Bathe, K.J., “A Continuum Mechanics Based Four-node Shell Element for General Non-linear Analysis,”
Engineering Computations, Vol.1, pp.77-88, (1984).
1 g 2 g 1 2 4 3 1 g 2 g 1 2 4 C 3 A D B 21
1 g 2 g 1 2 1 2 4 3
Dvorkin, E.N. and Bathe, K.J., “A Continuum Mechanics Based Four-node Shell Element for General Non-linear Analysis,”
Engineering Computations, Vol.1, pp.77-88, (1984).
1 g 2 g 1 2 1 2 4 3
せん断ひずみ成分
1g
2 g 1 2 1 2 4 C 3 A D B 22( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 2 x x y y z z
u
u
u
u
e
e
e
V
V
Fig. MITC4 shell element
(α) (α) (α) 1 2 3 V , V , V
: 節点での単位直交ベクトル
: ディレクタベクトル
a: 厚さ
位置ベクトル
変位ベクトル
(微小変形を仮定)
: 共変基底ベクトル
節点あたり
5自由度
(α) 3 V 0 0 (α) 3 3 0 3| g V | g 1g
2 g 3g
3 1 (3) 1 V (3) 2 V (3) 3 V (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V 2 x e y e z e x y zLee, P.S. and Bathe, K.J., “Development of MITC Isotropic Triangular Shell Finite Elements, Computers & Structures, Vol.82, pp.945-962, (2004). 1 2 3= x x x g , g , g 3 ( ) ( ) ( ) 3 1
( , )
2
a
N
x
x
V
3 ( ) ( ) ( ) 0 ( ) 3 3 1 3 ( ) ( ) ( ) 0 ( ) 3 1( , )
(
)
2
( , )
(
)
2
a
N
a
N
u
u
V
V
u
V
0 3 ( ) ( ) ( ) 0 ( ) 3 3 1( , )
(
)
2
a
N
u
x
x
u
V
V
23Lee, P.S. and Bathe, K.J., “Development of MITC Isotropic Triangular Shell Finite Elements, Computers & Structures, Vol.82, pp.945-962, (2004). B C 1
g
2 g 1 2 4 Aせん断ひずみ成分
24( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) 1 1 2 2 x x y y z z
u
u
u
u
e
e
e
V
V
Fig. MITC4 shell element
(α) (α) (α) 1 2 3 V , V , V
: 節点での単位直交ベクトル
: ディレクタベクトル
a: 厚さ
位置ベクトル
変位ベクトル
(微小変形を仮定)
: 共変基底ベクトル
節点あたり
5自由度
(α) 3 V 0 0 (α) 3 3 0 3| g V | g x e y e z e x y z 1 2 (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (3) 1 V (3) 2 V (3) 3 V (4) 1 V (4) 2 V (4) 3 V 3 4 1 g 2 g 3 g 5 6 8 7 9 (6) 1 V (6) 2 V (6) 3 V (5) 1 V (5) 2 V (5) 3 V (7) 1 V (7) 2 V (7) 3 V (8) 1 V (8) 2 V (8) 3 V (5) 1 V (5) 2 V (5) 3 V 1 2 3= x x x g , g , g 9 ( ) ( ) ( ) 3 1( , )
2
a
N
x
x
V
9 ( ) ( ) ( ) 0 ( ) 3 3 1 9 ( ) ( ) ( ) 0 ( ) 3 1( , )
(
)
2
( , )
(
)
2
a
N
a
N
u
u
V
V
u
V
0 9 ( ) ( ) ( ) 0 ( ) 3 3 1( , )
(
)
2
a
N
u
x
x
u
V
V
Bucalem, M.L. and Bathe, K.J., “Higher-order MITC general shell element,” International Journal for Numerical Methods in
Engineering, Vol.36, pp.3729-3754, (1993).
Bucalem, M.L. and Bathe, K.J., “Higher-order MITC general shell element,” International Journal for Numerical Methods in
Engineering, Vol.36, pp.3729-3754, (1993). 1
g
2 g 1 2 1 5 2 4 7 3 8 9 6 1g
2 g 1 2 1 5 2 4 7 3 8 9 6 1g
2 g 1 2 1 5 2 4 7 3 8 9 6せん断ひずみ成分
260 0 0 0
ˆ
ˆ
ˆ
ˆ
ˆ
ijkl i j k l ijkl i j k lC
C
C
g
g
g
g
e
e
e
e
1 2 3 ˆ ˆ ˆ e , e , e: 単位直交ベクトル
1 2 3= x x x g , g , g: 共変基底ベクトル
弾性テンソル
k : せん断補正係数
E
: Young率
: Poisson比
3 2 3 3 1 2 3 1 3 2 3 ˆ ˆ = , ˆ , ˆ ˆ ˆ ˆ | | | g g e e e = e e e | g g e 1 g 2 g 3 g 1 2 3 1 2 (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (3) 1 V (3) 2 V (3) 3 V (4) 1 V (4) 2 V (4) 3 V 3 4 1 ˆe 2 ˆe 3 ˆe (α) (α) (α) 1 2 3 V , V , V: 節点での単位直交ベクトル
各節点において
(α) (α) 3 1 (α) (α) (α) 3 3 2 1 2 3 3 1 ˆ ˆ , , ˆ | | e e V e V V = V V e e弾性テンソル成分の変換
1111 1122 1112 1123 1131 2211 2222 2212 2223 2231 1211 1222 1212 1223 1231 2 2311 2322 2312 2323 2331 3111 3122 3112 3123 3131 1 0 0 0 ˆ ˆ ˆ ˆ ˆ 1 0 0 0 ˆ ˆ ˆ ˆ ˆ 1 0 0 0 0 ˆ ˆ ˆ ˆ ˆ 2 1 ˆ ˆ ˆ ˆ ˆ 1 0 0 0 ˆ ˆ ˆ ˆ ˆ C C C C C C C C C C E C C C C C C C C C C k C C C C C 0 2 1 0 0 0 0 2 k 0 0 0 0 0 0 0 0 0 0 0 0(
:
) :
ˆ
(
ˆ
)(
ˆ
)(
ˆ
)(
ˆ
)
ijkl i j k l mnop i j k l m n o pC
C
C
g
g
g
g
e
g
e
g
e
g
e
g
27(A) 隣接要素の節点で同じディレクタベクトルを使用する場合
(B) 隣接要素の節点で別のディレクタベクトルを使用する場合
(1) 5自由度として計算する場合
(ディレクタ回りの回転自由度を計算しない場合)
(2) 6自由度として計算する場合
(ディレクタ回りの回転自由度をdrilling DOFとして
計算する場合
)
FrontISTRでは，(2) と (B) の定式化を使用しています
28シェル要素の節点あたりの自由度数を
6として，六つ目の自由度 (
drilling
DOF
と呼ばれる三つ目の回転自由度
) を考慮 (Hughes and Brezzi 1989)
( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) x x y y z z x x y y z zu
u
u
u
e
e
e
e
e
e
t ** 31 DI 23 DI 31 31 23 23 0 0 3 31
:
(
)
(
)
2
1
1
:
:
2
2
2
V S V V V VdV
dS
dV
dV
dV
dV
u t
u
b
V
e
V
e
drilling DOFを考慮した場合のエネルギー
( /
0.0001 1.0)
( ) 1 2 ( ) 3 ( ) 0 ( ) 3 3 1 ( ) 1 2 ( ) 3 ( ) 0 ( ) 3 1( ,
)
(
)
2
( ,
)
(
)
2
n na
N
a
N
u
u
V
V
u
V
0 0 T
1
(
)
2
u
u
e
e
ijke
i
e
je
kHughes, T.J.R. and Brezzi, F., Computer Methods in Applied Mechanics and Engineering, Vol.72, pp.105-121, (1989).
Nguyen-Van, Hieu and Mai-Duy, Nam and Tran-Cong, Thanh, CMES: Computer Modeling in Engineering and Sciences,
t 0 0 3 3
1
1
:
:
:
2
2
V V S VdV
V
V
dV
dS
dV
e
e
u t
u
b
**0
より，仮想仕事の原理
AS DI DI 31 31 31 31 AS DI DI 23 23 23 231
1
(1
)
(0, 1, 0)
(1
)
(0, 1, 0)
2
2
1
1
(1
)
( 1, 0, 0)
(1
)
(1, 0, 0)
2
2
( /
0.0001 1.0)
300 3
V
0 3V
0 3V
0V
3 0 3V
0 3V
0 3V
(A) 隣接要素の節点で
同じディレクタベクトル
を使用する場合
(B) 隣接要素の節点で
別のディレクタベクトルを
使用する場合
3132
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・ 理論
・
現バージョンのプログラム
・ 解析事例
２．シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
33
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・ 理論
・ 現バージョンのプログラム
・
解析事例
２．シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
シェル要素とソリッド要素の計算結果の比較
集中荷重
1mN
集中荷重
1mN
薄膜構造の円筒
半径
100mm
高さ
200mm
Pinched cylinder model
(Bucalem and Bathe 1993)
メッシュの規模による計算精度の違い
Bucalem, M.L. and Bathe, K.J., International Journal for Numerical Methods in Engineering,
Vol.36, pp.3729-3754, (1993).
34「奥田洋司
, 早田浩平, 橋本学, 上島豊, “クラウドCAEシステムを用いた効率的な有限要素モデリング,”
42
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
２．
シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
43
要素混在問題における
3×3BCSR形式での格納
1 g 2 g 3 g 1 3 (3) (1) (3) 1 V (3) 2 V (3) 3 V (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (2) 2 1 g 2 g 3 g 1 3 (1) (2) (1) 1 V (1) 2 V (1) 3 V (2) 1 V (2) 2 V (2) 3 V (3) 1 V (3) 2 V (3) 3 V (4) 1 V (4) 2 V (4) 3 V (3) (4) 2 ( ) ( ) ( ) x y zu
u
u
( ) ( ) ( ) ( ) 1 ( ) 2 x y zu
u
u
( ) ( ) ( ) ( ) ( ) ( ) x y z x y zu
u
u
Hughes, T.J.R. and Brezzi, F., Computer Methods in Applied Mechanics and Engineering, Vol.72, pp.105-121, 1989.
Nguyen-Van, Hieu and Mai-Duy, Nam and Tran-Cong, Thanh, CMES: Computer Modeling in Engineering and Sciences, Vol.49, No.2, pp.81-110, 2009.
11 12 1 3 3 3 3 3 3 21 22 2 3 3 3 3 3 3 1 2 3 3 3 3 3 3 N N N N NN
K
K
K
K
K
K
K
K
K
K
drilling DOF
[7][8] ( ) ( ) ( ) x y z
( ) ( ) ( ) x y zu
u
u
3 DOFs at each node
: Number of original and dummy nodes
N
Dummy node
5 DOFs
at each node
6 DOFs at each node
3 DOFs at each node
Stiffness matrix (sparse matrix)
FrontISTRの既存の要素データの利用 (1/2)
44
611梁要素
6○○
シェル要素
7○○
3次元ソリッド要素
3○○
ソリッド要素と梁要素／シェル要
素を混在させる場合，
節点あたり
3自由度
の
データ構造に統一して，
剛性マトリックスを作成し，
線形ソルバーに渡すようにする
梁要素のデータ構造が
節点あたり
3自由度を持つように，
四面体要素
341のデータ構造を
利用した要素タイプ
641
を
導入する
節点あたりの自由度
・
3次元ソリッド要素：3
・ シェル要素：
6
・ 梁要素：
6
FrontISTRの既存の要素データの利用 (2/2)
45
611梁要素
6○○
シェル要素
7○○
3次元ソリッド要素
3○○
ソリッド要素と梁要素／シェル要素
を混在させる場合，
節点あたり
3自由度
の
データ構造に統一して，
剛性マトリックスを作成し，
線形ソルバーに渡すようにする
MITC4シェル要素のデータ構造が
節点あたり
3自由度を持つように，
六面体要素
361のデータ構造を
利用した要素タイプ
761
を
導入する
節点あたりの自由度
・
3次元ソリッド要素：3
・ シェル要素：
6
・ 梁要素：
6
要素タイプ
641 (1)
46
(2) 1 (2) 2 (2) 3 x x x (1) 1 (1) 2 (1) 3 x x x (1)
(2)
(2)
(1)
(3)
(4)
(1) 1 (1) 2 (1) 3 x x x (1) 1 (1) 2 (1) 3 x x x (2) 1 (2) 2 (2) 3 x x x (2) 1 (2) 2 (2) 3 x x x FrontISTRのメッシュファイルでは，節点座標を二つ用意する
(メッシュを表示する際に都合が良いためであるが，
計算には使用しないので
(3) と (4) の節点座標はどのような値でも良い)
計算したい梁要素
611
入力データとして用意するのは
四面体要素
431のデータ構造
要素タイプ
641 (2)
47
FrontISTRのプログラム内では，(3) と (4) の節点変位は梁の回転自由度となる
(2) 1 (2) 2 (2) 3 (2) 1 (2) 2 (2) 3 u u u (1) 1 (1) 2 (1) 3 (1) 1 (1) 2 (1) 3 u u u (1)
(2)
計算したい梁要素
611
(1) 1 (1) 2 (1) 3 u u u (2)
(1)
(3)
(4)
(2) 1 (2) 2 (2) 3 u u u (2) 1 (2) 2 (2) 3 (1) 1 (1) 2 (1) 3 FrontISTRのプログラム内では
四面体要素
431のデータ構造
要素タイプ
641 (3)
48
1,1 1,2 1,11 1,12 2,1 2,2 2,11 2,12 11,1 11,2 11,11 11,12 12,1 12,2 12,11 12,12K
K
K
K
K
K
K
K
K
K
K
K
K
K
K
K
K =
剛性マトリックス
(梁要素611)
1,1 1,2 1,11 1,12 2,1 2,2 2,11 2,12 11,1 11,2 11,11 11,12 12,1 12,2 12,11 12,12K
K
K
K
K
K
K
K
K
K
K
K
K
K
K
K
K =
FrontISTRのプログラム内では，
梁要素の剛性マトリックス成分をテーブルに従って変更する
(2) 1 (2) 2 (2) 3 (2) 1 (2) 2 (2) 3u
u
u
(1) 1 (1) 2 (1) 3 (1) 1 (1) 2 (1) 3u
u
u
(1) 1 (1) 2 (1) 3u
u
u
(2) 1 (2) 2 (2) 3u
u
u
(1) 1 (1) 2 (1) 3
(2) 1 (2) 2 (2) 3
1 2 3 4 5 6 7 8 9 10 11 12 1 2 3 4 5 6 7 8 9 10 11 12テーブル
4,8 7,5K
K
例えば，
FrontISTRのプログラム内で計算する
剛性マトリックス
49
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
２．
シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
50
目次
１．
FrontISTRにおけるBernoulli-Euler梁要素／MITCシェル
要素を用いた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
２．
シェル要素／梁要素とソリッド要素を混ぜた解析
・ 理論
・ 現バージョンのプログラム
・ 解析事例
梁の曲げ解析
(解析モデル)
51
T
xx ,max= 60 MPa
T
xx ,min= －60 MPa
x
y
L
H
M
T
xx,maxT
xx,minYoung’s modulus: E = 200,000 MPa
Poisson’s ratio: ν = 0
Geometry: L = 100 mm, H = 10 mm
O
梁の曲げ解析
(メッシュ)
52
611 411 111 711 211 811 511 311 11 601 701 801 501 301 1 11 40 5 10 5 703 603 20 5 307 503 803 303 605 705 805 105 205 305 405 5 505 507 707 807 309 509 809 109 409 609 607 709 9 209 1 11Case A : 梁1次要素
Case B : ソリッド1次要素 (非適合要素)
Case C : 梁1次要素とソリッド1次要素の混合
梁の曲げ解析
(梁1次要素とソリッド1次要素の混合)
53
x x z y y y x z z xu
u
y
z
u
u
z
u
u
y
中立軸上の節点
5の変位を
，回転自由度を
とすると，
同じ断面上にある節点
105，205，305，405，505，605，705，805の
変位
を
以下の式によって拘束する
(ただし，微小変形を仮定)
x
y
z
40 5 10 5 703 603 20 5 307 503 803 303 605 705 805 105 205 305 405 5 505 507 707 807 309 509 809 109 409 609 607 709 9 209 1 11梁の中立軸と
x 軸を一致させる
( ,
u u u
x y,
z)
( ,
x y,
z)
( ,
u u u
x y,
z
)
梁の曲げ解析
(梁1次要素とソリッド1次要素の混合)
54
105 5 5 205 5 305 5 5 405 5 5 505 5 5 605 5 5 705 5 805 5 52
2
2
2
2
2
x x z x x x x z x x z x x z x x z x x x x zH
u
u
u
u
H
u
u
H
u
u
H
u
u
H
u
u
u
u
H
u
u
梁の曲げ解析
(計算結果：Case AとCase B)
55
0.0E+00 5.0E-02 1.0E-01 1.5E-01 2.0E-01 2.5E-01 3.0E-01 3.5E-010.0E+00 2.0E+01 4.0E+01 6.0E+01 8.0E+01 1.0E+02
uy [mm] x [mm] Exact Case A
変形の拡大率 ×100
0.0
0.3 [mm]
Case Bの変形形状
Displacement
0.0E+00 5.0E-02 1.0E-01 1.5E-01 2.0E-01 2.5E-01 3.0E-01 3.5E-010.0E+00 2.0E+01 4.0E+01 6.0E+01 8.0E+01 1.0E+02
uy
[mm]
x [mm] Exact
梁の曲げ解析
(計算結果：Case C)
56
変形の拡大率 ×100
0.0
0.3 [mm]
Case Cの変形形状
Displacement
0.0E+00 5.0E-02 1.0E-01 1.5E-01 2.0E-01 2.5E-01 3.0E-01 3.5E-010.0E+00 2.0E+01 4.0E+01 6.0E+01 8.0E+01 1.0E+02
uy
[mm]
x [mm] Exact
