静的弾性問題の有限要素法解析アルゴリズム

概要
基礎理論
1.応力とひずみおよび平衡方程式
2.降伏条件式
3.構成式（応力－ひずみ関係式）
 有限要素法
1.有限要素法の概要
2.仮想仕事の原理式と変分原理
3.平面ひずみ弾性有限要素法定式化

FEMの基礎方程式
1.
平衡方程式
  x  xy  zx
 x  y  z  Gx  0 
  xy  y  yz

 G y  0 　

x

y

z

 
 yz  z
zx

 Gz  0 　

x

y

z

2.　ひずみ－変位関係式
u  xy  1  xy  1  u  v 

 x  x 
2
2  y x 


v 
1
1  v w 

 y 




    yz
yz

y
2
2
 z y 



w 
1
1  w u 


 z
 zx   zx    
z 
2
2  x z 
3.　応力－ひずみ関係式
(構成式)　 ij  Dijkl kl
4.　境界条件式　 力の境界条件　ti  Pi on　St  変位の境界条件　ui  Vi on　Su 応力とひずみおよび
平衡方程式
物体にはたらく力と応力
応力の定義
応力ベクトル：
 P dP
Tn  lim

A0  A dA
垂直応力：
 N dN
 n  lim

A0  A dA
せん断応力：
 Q dQ
 n  lim

A0  A dA
応力ベクトル
Tx   xx e x   xy e y   xz e z

Ty   yx e x   yy e y   yz e z

Tz   zx e x   zy e y   zz e z
テンソル標記で
Ti 
 e
ij j
j
  ij e j
応力テンソル
 xx

 ij     yx
 zx
 x

  xy
 zx
 xy  xz 

 yy  yz 
 zy  zz 
 xy  zx 

 y  yz 
 yz  z 
モーメントの釣合
＝応力テンソルの対称性
 xy   yx   xy
 ij   ji
二次元応力行列と主応力
Tx   x  yx  nx 
nx 
  
    

Ty   xy  y  n y 
n y 
 x  
 
 xy
 yx  nx  0
  

 y    n y  0
 x 
 yx
0
 xy
 y 
 2  I1  I 2  0  (  1 )(   2 )  0  主応力1 , 2
三次元応力の座標変換
Tnx   x l   yx m   zx n

Tny  xy l   y m   zy n

Tnz  xz l   yz m   z n
Tni 

ji n j
  ji n j
j
 n  Tnxl  Tny m  Tnz n
 n 2  Tn 2   n 2
三次元応力行列と主応力
今考えている面を主応力  がはたらく主応力面とすると
Tni   ji n j   ni    ij n j
( ij    ij ) n j  0
上式が nj=0 以外の解をもつためには
det ( ij    ij )  0
上式を展開すると
 3  I1 2  I 2  I 3  0
ここで J1，J2，J3 を応力の不変量という．
上式の3実根を 1 ，2，3 とすれば
   1    2    3   0  主応力  1 , 2, 3
三次元応力の不変量
 I1   x   y   x   ii  3 m

1

2
2
2
 ii  jj   ij ij 
I
















 2
x y
y z
z x
xy
yz
zx
2

 I 3   x y z  2 xy yz zx   x yz 2   y zx 2   z xy 2  det  ij 

あるいは主応力を用いて表すと
 I1   1   2   3

 I 2   1 2   2 3   3 1
I    
 3
1 2 3
平均垂直応力と偏差応力
平均垂直応力＝静水応力 ⇒塑性変形に無関係
m 
 x  y  z
3
1
1
  ii  J1
3
3
偏差応力 ⇒塑性変形を引き起こす
 x '   x   m
 '    
y
m
 y
,

 z '   z   m
 xy ,  yz ,  zx

 1 '   1   m

 2 '   2   m   ij '   ij   m  ij
 '    
 3
3
m
偏差応力の不変量
 J1   x ' y ' z '   1 ' 2 ' 3 '  0

 J 2  1  x '2  y '2  z '2  2( xy 2   yz 2   zx 2 )
2


1
2
2
2




'

'


'

'


'

'



'


'


'

1 2
2
3
3 1
1
2
3

2

  1    2     2     2  6( 2   2   2 )
x
y
y
z
z
x
xy
yz
zx

6

1
   1   2 2   2   3 2   3   1 2
6

 J 3   1 ' 2 ' 3 '









二次元x方向応力の平衡方程式
（＝釣合方程式）
  yx 


 x 
dx dy   x dy   yx 
dy dx   yx dx  Fx dxdy  0
 x 
x
y




  x   yx

 Fx  0
x
y
三次元応力の平衡方程式
（＝釣合方程式）
   x   yx   zx
  x   y   z  Fx  0

   xy   y   zy


 Fy  0

y
z
 x

  yz   z
xz


 Fz  0

y
z
 x

 ji , j  Fi  0
ひずみ（微少ひずみ）の定義
垂直ひずみ
（＝垂直微少ひずみ）
 yy 
u y  B   u y  A
dy
 u y  A
u y  A 
dy  u y  A
y

dy
 uy

 y
y

 ux
 xx   x   x

 uy

 y
 yy 
y


 uy
 z
 zz 
y

せん断ひずみ
（＝微少せん断ひずみ）
1
1
2
2
1  u z  B   u z  A u y  D   u y  A 
 


2
dy
dz

 yz       tan   tan  
u y
 u  u z dy  u

 z
u

dz

u
z
y
y
1
y
z



2
dy
dz





1   uz  u y  1
 

   yz
2  y  z  2

 ux  u y 1
 xy   y   x  2  xy

 u y  uz 1


  yz
 yz 
z y 2


 uz  ux 1



  zx
 zx
x z 2

ここで ij は
工学的せん断ひずみ
ひずみテンソル
(ひずみ－変位の関係式）
1   ui  u j  1
 ij  

 ui , j  u j ,i 

2   x j  xi  2
 xx

     xy
 zx
 xy
 yy
 yz
 
x
 zx  
 1
 yz     xy
2

 zz  1
  zx
 2
1
 xy
2
y
1
 yz
2
1 
 zx 
2
1 
 yz 
2 
z 

体積ひずみと偏差ひずみ
体積ひずみ
V 

ii
  ii   xx   yy   zz   x   y   z
i
偏差ひずみ
 '   '    1     1 
x
xx
V
x
V
 xx
3
3

 yy '   y '   yy  1  V   y  1  V
1
3
3
  ij '   ij   ij  V

3

1
1
 zz '   z '   zz   V   z   V
3
3

 xy '   xy ,  yz '   yz ,  zx '   zx
降伏条件
単軸応力状態の降伏条件
多軸応力状態の降伏条件 (1)
多軸応力状態の降伏条件 (2)
主せん断応力と最大せん断応力
 1 ,  2 ,  3 : 主応力
 1 ,  2 ,  3 : 主せん断応力
  1   
3
1 2 2

1

 2   3   1
2

1

 3  2  1   2
 max :最大せん断応力
 max  max  1 , 2 , 3 
1   3
 1   2   3   max 
2
Trescaの降伏条件 (1864)
最大せん断応力 max が材料固有の臨界値CT
（せん断降伏応力）に
達したとき降伏する．
 max  max  1 , 2 , 3 
1
1
1

 max   1   2 ,  2   3 ,  3   1   CT  kT
2
2
2

あるいは 1   2   3 のとき
 max 
1   3
 CT  kT
2
ここでkT : せん断降伏応力
Trescaの降伏条件における
臨界値の決定
単軸引張試験の降伏応
力を  Y とすると，
 1   Y , 2   3  0
1  3  Y
 CT 

2
2
純粋せん断の降伏状態
にあるとき
 1  kT ,  2  0,  3  kT
1   3
 CT 
 kT 
Y
2
2
 kT
Misesの降伏条件 (1913)
材料中のせん断ひずみ
エネルギー＝偏差応力
の2次の
不変量 J 2 が材料固有の臨界値CM に達したとき降伏．


1
1
J 2   ij ' ij '   1 '2  2 '2  3 '2
2
2
1
  1   2 2   2   3 2   3   1 2  CM  k M 2
6
ここでk M : はせん断降伏応力


Misesの降伏条件における
臨界値の決定
単軸引張試験の降伏応
力を  Y とすると，
 1   Y , 2   3  0
 CM 
Y 2
3
純粋せん断の降伏状態
にあるとき
 1  k M ,  2  0,  3  k M
 CM  k M 2
 kM 
Y
3
降伏曲面・降伏曲線
主応力空間における降伏曲面
π平面上の降伏曲線
降伏条件式の実験的検証
応力－ひずみ関係式
＝構成式
弾性体の構成式 (1)
（一般化されたフックの法則）
1
1 e  xy
e
  x   y   z ,  xy   xy 
E
2
2G
1
1 e  yz
e
e
 y   y   z   x ,  yz   yz 
E
2
2G
1
1 e  zx
e
e
 z   z   x   y ,  zx   zx 
E
2
2G
ここでE はヤング率，
G は横弾性係数，
 はポアソン比
 xe
である．テンソル標記
では
1 

1

e
 ij 
 ij   kk ij 
 ij   kk ij
E
E
2G
E
弾性体の構成式 (2)
（一般化されたフックの法則）
あるいはその逆関係と
して









E
e
e
(1  ) xe   ye   ze ,  xy  2G xy
 G xy
(1  )(1  2 )
E
e
e
y 
(1  ) ye   ze   xe ,  yz  2G yz
 G yz
(1  )(1  2 )
E
e
e
z 
(1  ) ze   xe   ye ,  zx  2G zx
 G zx
(1  )(1  2 )
x 
テンソル標記では
 ij  2G   ik  jl   il jk  
1
2

1  2
e
 ij kl  kle  Dijkl
 kle

弾性体の構成式 (3)
（一般化されたフックの法則）
またフックの法則を偏
差応力を用いて表すと
 xy
1


m ,
  xy 
2G
E
2
2G
 y ' 1  2
 yz
1
e
e
y 

 m ,  yz   yz 
2G
E
2
2G

 ' 1  2
1
 ze  z 
 m ,  zxe   zx  zx
2G
E
2
2G
テンソル標記では
 x ' 1  2
 xe
 ije

 ij ' 1  2
2G

E
e
 xy
 m ij
Reussの構成式
”塑性ひずみ増分
d ijp の方向は偏差応力 ij '
の方向に一致する”と
仮定した塑性構成式
p
p
p
d xp d y d zp d xy d yz d zxp





 d
 x '  y '  z '  xy
 yz
 zx
テンソル標記すると
d ijp   ij ' d
上式は塑性体積一定の
条件を満足している．
d xp  d yp  d zp   x ' y ' z 'd  0
剛塑性体の構成式
（Levy-Misesの式）
上式を変形し，一般応
力成分で表すと
2
1
1

p
  x   y   z d , d xy  d xy   xy d
3
2
2

2
1
1

p
p
d y   y   z   x d , d yz  d yz   yz d
3
2
2

2
1
1

p
p


d z   z   x   y d , d zx  d zx   zx d
3
2
2

d xp
弾塑性体の構成式
（Prandtle-Reussの式）
全ひずみ増分d ij  弾性ひずみ増分d ije  塑性ひずみ増分d ijp
d xy d xy
d x ' 1  2
d x 

d m   x ' d , d xy 

  xy d
2G
E
2
2G
d y ' 1  2
d yz d yz
d y 

d m   y ' d , d yz 

  yz d
2G
E
2
2G
d
d
d z ' 1  2
d z 

d m   z ' d , d zx  zx  zx   zx d
2G
E
2
2G
テンソル標記すると
d ij ' 1  2
e
p
d ij  d ij  d ij 

d m ij   ij ' d
2G
E
相当応力と相当塑性ひずみ増分
降伏応力に
降伏応力の程度を単軸
多軸応力状態における
ーゼス
数値  を相当応力と呼び，ミ
換算して評価できる関
る
材では次式のようにな



1
3
 1   2 2   2   3 2   3   1 2
 1 '2  2 '2  3 '2 
2
2
1
 x   y 2   y   z 2   z   x 2  6  xy 2   yz 2   zx 2

2





るとき　塑性仕事増分dW p に関して次式が成立す
dW p   ij d ij p   ij ' d ij p   d p
れる
といい，次式で定義さ
d p を相当塑性ひずみ増分
d p 

2
p

d
x
3 
  d   d 
p 2
p 2
2
y
z


1
p
d xy
2
  d   d  
p 2
p 2
2
yz
zx
二次元平面ひずみ
弾性有限要素法
有限要素法とは


FEM＝Finite Element Method
解析対象物体(連続体)を有限個の要素に分割し，各要
素について剛性方程式を構成し，それらを全要素につい
て重ね合わせる
固体力学解析用有限要素法


弾塑性有限要素法
・弾性有限要素法（静的陽解法）
・微少変形弾塑性有限要素法
（静的陽解法・静的陰解法）
・大変形弾塑性有限要素法
（静的陽解法・静的陰解法・動的陽解法）
剛塑性有限要素法
（静的陰解法）
弾性FEM定式化の流れ
(1) 釣合方程式
(3) 変分原理
ガウスの発散定理
ポテンシャル
停留の原理
(5) 形状関数
(2) 仮想仕事の原理式
(4) 構成方程式
離散化
(7) 有限要素方程式
(6) ひずみ－変位関係式
弾性FEMの基礎方程式
＝弾性境界値問題
1.
平衡方程式
  x  xy  zx
 x  y  z  Gx  0 
  xy  y  yz

 G y  0 　

x

y

z

 
 yz  z
zx

 Gz  0 　

x

y

z

3.　応力－ひずみ関係式
(構成式)　e
 ij  Dijkl
 kl

 ij 
E
E
 ij 
 kk ij
1 
1  1  2 
Ue
1 e
, U e  Dijkl
 ij  kl
 ij
2
2.　ひずみー変位関係式
u



 x x

v

 y 
y


w


 z
z

1
1  u v 




  
xy
xy

2
2  y x 


1
1  v w 




    yz
yz
2
2
 z y 


1
1  w u 
 zx   zx    
2
2  x z 

4.　境界条件式　 力の境界条件　ti  Pi on　St  変位の境界条件　ui  Vi on　Su 弾性FEM定式化の流れ
(1) 釣合方程式
(3) 変分原理
ガウスの発散定理
ポテンシャル
停留の原理
(5) 形状関数
(2) 仮想仕事の原理式
(4) 構成方程式
離散化
(7) 有限要素方程式
(6) ひずみ－変位関係式
仮想仕事の原理式
静的可容応力：平衡方程式と力学的境界条件を満足する応力
動的可容変位：ひずみ－変位関係式と幾何学的境界条件を満足する変位
仮想変位：動的可容変位の変分
静的可容応力と仮想変位に対して次式が成り立つ．
 (
ji , j

 Gi )ui dS  (ti  Pi )ui dS  0
V
St
上式にガウスの発散定理を適用すると次の仮想仕事の原理式を得る

V
ji ij dV

 G u dV  Pu dS　i
V
i
i
i
St
可容応力と仮想変位によってなされる内部仕事が外部仕事に等しいことを表す．
変分原理
仮想仕事の原理式は弾性体の全ポテンシャルエネルギΦの第一変分が零である
ことを表しているポテンシャルエネルギ停留の原理に置き換えることができる．



  0 ,   U e dV  Pi ui dS  Gi ui dV
V
St
V
今，真の変位をui，それからわずかに異なる任意の可容変位をui+uiとすると，
 2U e
ij kl dV
V  ij  kl

1
ui  ui   ui  
2
ひずみエネルギ関数Ueが正値2次形式の場合，上式右辺第2項は正であるから
 ui  ui    ui 
となり，真の変位に対するポテンシャルエネルギは最小値をとる．
弾性FEM定式化の流れ
(1) 釣合方程式
(3) 変分原理
ガウスの発散定理
ポテンシャル
停留の原理
(5) 形状関数
(2) 仮想仕事の原理式
(4) 構成方程式
離散化
(7) 有限要素方程式
(6) ひずみ－変位関係式
2次元平面ひずみ変形状態の
ひずみと応力
11 12
   21  22
0
 0
0
0

0
0 
 11  12
   21  22 0 
0  33 
 0
 33
E
 11   22 
(11   22 ) 
(1  )(1  2 )

平面ひずみ変形状態における
応力－ひずみ関係式

1

 x 
(1  v) E  
 

 y  
  (1  )(1  2 ) 1 
 z
 0

または
   D 

1 
1
0

0 
 
 x 
0  y 


1  2  2 xy 
2(1  ) 
弾性FEM定式化の流れ
(1) 釣合方程式
(3) 変分原理
ガウスの発散定理
ポテンシャル
停留の原理
(5) 形状関数
(2) 仮想仕事の原理式
(4) 構成方程式
離散化
(7) 有限要素方程式
(6) ひずみ－変位関係式
三角形3節点要素と形状関数
形状関数 N  x, y  の性質
i  N  x, y  は節点 で 1 ，それ以外
の2つの節点で0 の値をとる．
ii  N  x, y  は線形の関数である．
形状関数の具体形
u ( x, y )  N1u1  N 2u2  N 3u3  N u
v( x, y )  N1v1  N 2v2  N 3v3  N v
あるいはマトリックスの形で
u   N1
 
v   0
または
u  N d
0
N2
0
N3
N1
0
N2
0
 u1 
v 
 1
0  u2 
 

N 3   v2 
u3 
 
 v3 
形状関数の計算
1
N1 
( x2 y3  x3 y2 )  ( y2  y3 ) x  ( x3  x2 ) y
2
1
N2 
( x3 y1  x1 y3 )  ( y3  y1 ) x  ( x1  x3 ) y
2
1
N3 
( x1 y2  x2 y1 )  ( y1  y2 ) x  ( x2  x1 ) y
2
ただし，
 x , y  は節点  におけるx, y の座標値であり，
1 x1
2  det 1 x2

1 x3
y1 
y2   2  (要素の面積)

y3 
弾性FEM定式化の流れ
(1) 釣合方程式
(3) 変分原理
ガウスの発散定理
ポテンシャル
停留の原理
(5) 形状関数
(2) 仮想仕事の原理式
(4) 構成方程式
離散化
(7) 有限要素方程式
(6) ひずみ－変位関係式
ひずみ－変位マトリックス (1)
（Bマトリックス）
x 
y 
2 xy
N
u N1
N

u1  2 u2  3 u3 
x x
x
x
N
v N1
N

v1  2 v2  3 v3 
y y
y
y




 N
u
x
 N
v
y
u v
  xy  
y x
N
N
N
N
N
N
 1 u1  1 v1  2 u2  2 v2  3 u3  3 v3
y
x
y
x
y
x



N 
 N
u

v 


x
 y

ひずみ－変位マトリックス (2)
（Bマトリックス）
マトリックスの形式で書くと
 N1

  x    x   x

   

 y    y   0
2    
 xy   xy   N1

 y
0
N1
y
N1
x
N 2
x
0
N 2
y
0
N 2
y
N 2
x
N 3
x
0
N 3
y
u 
 1 
0  v1
 
 u 
N 3  2
 
y  v2 
N 3  u 
 3 
x  v
 
3
ひずみ－変位マトリックス (3)
（Bマトリックス）
N 
1
(a b x  C y )
2
の形になっているから，そのｘおよびｙに関する勾配は
N
N
1
1

b ,

c
x
2
y
2
ただし
b1  y2  y3 , b2  y3  y1 , b3  y1  y2
c1  x3  x2 , c2  x1  x3 , c3  x2  x1
と書けるので，これをひずみ－変位マトリックスに代入すると
ひずみ－変位マトリックス (4)
（Bマトリックス）
したがってひずみ－変位関係式は
  x    x  b1 0 b2

   

 y     y    0 c1 0
2    c b c
 xy   xy   1 1 2
または
さらに
   Bd
   D    DB d
0
b3
c2
0
b2
c3
 u1 
v 
0 1
u2 

c3  
 v2
b3   
u3 
 
 v3 
弾性FEM定式化の流れ
(1) 釣合方程式
(3) 変分原理
ガウスの発散定理
ポテンシャル
停留の原理
(5) 形状関数
(2) 仮想仕事の原理式
(4) 構成方程式
離散化
(7) 有限要素方程式
(6) ひずみ－変位関係式
離散化(要素剛性方程式) (1)
三角形3節点要素について，仮想仕事の原理式の左辺（内部仕事）は
   dV
    

ij
V
V
ij
11

 x

V


11
22 22
 y
 T  dV
V
 2 1212 dV
 x 
 
2 xy   y dV
 
 xy 
離散化(要素剛性方程式) (2)
仮想仕事の原理式の右辺（外部仕事）は


bi  ui dV  ti  ui dS
V

S

 (b1  u1  b2  u2 ) dV  (t1  u1  t 2  u2 ) dS
V
S
bx 
 u y    dV 
by 
  u
  u
  u bdV   u tdS



V
x
T
V
T
S
V
x
t x 
 u y    dS
t y 
離散化(要素剛性方程式) (3)
ここで以下の関係式がる
   D   DBd
 u  N  d
    B d
よって三角形3節点要素に関する仮想仕事の原理式は

V

T
T


d  B  D B d dV
e

V



d T [ N ]T bdV 
e
S


d T N T tdS
e
離散化(要素剛性方程式) (4)
ここで仮想変位は定数であり，積分の外に出してもよいので
 d T 


Ve
BT DBdV d  
Ve

[ N ]T bdV 
任意の仮想変位に対して上式が成立するためには [


Ve
BT DBdV d  
Ve
Se
N T tdS   0

] 内は常に0

[ N ]T bdV 
Se
N T tdS
これが解くべき剛性方程式である．左辺の積分内のマトリックスを

Ve
BT DBdV  BT DB  K e 
とおくとことにする．△ は三角形要素の面積である．
離散化(要素剛性方程式) (5)
仮想仕事の原理式の右辺第１項の物体力の項は

 N1
0

N2
T
N  bdV  e 
e
0
V
V

 N3

0

ただし物体力は要素内で一定と仮定
0
bx 
 
by 
N1 


0  bx 
 bx 
  dV   

N 2 by 
3 by 

bx 
0
 

N3 
by 
離散化(要素剛性方程式) (6)
右辺第１項表面力の項は，
例えば面2-3に右図のように
表面力が分布しているなら
形状関数マトリックスを
次のように書き直して
l

0
 u  0 0 1  L
 
l
 v  0 0
0
1

L
l
L
0
 u1 
 v 
 1
0   u 
2
l   v2 


L   u 
3


 v3 
離散化(要素剛性方程式) (7)
これより表面力の項は次式のようになる

0 
 0
 0
0 
0
 l

0
1 
0 
 
 L
 t
t x 

L
x
T
l
 0
  dS   

N  tdS 
1

2 t y 
Se
Se 
L  t y 
 l

t x 
0
 L

 

t y 
l 
 0


L 

ただし表面力は面2-3上で等分布荷重とした．
離散化(要素剛性方程式) (8)
最終的に要素剛性方程式は次式のように書き換えられる
K d f 
e
K 
e
e
: 要素剛性マトリックス
 u1 
v 
 1
u2 
d     : 節点変位
 v2 
u3 
 
 v3 
f e 
 f x1 
f 
 y1 
 f x2 
   : 節点力
 f y2 
 f x3 
 
 f y3 
全体剛性方程式
下図に示すような複数要素からなる系の全体系に関する仮想仕事の原理
のマトリックス表示は
 d 
T




f 0

K d  
e
e
e
e
これより全体系に関する剛性方程式は次のように得られる
K d   f 
K    K e 
e
f  
e

f
 
e
弾性有限要素法解析の流れ
領域の要素分割，境界条件の設定
Pre-Processor
要素剛性マトリックスの計算
全体剛性マトリックスの計算
等価節点力，変位拘束の導入
FEM Analysis
連立一次方程式を解き節点変位を求める
節点変位からひずみ，応力の計算
結果の出力，可視化
Pre-Processor

Download Report