数値計算のテクニックの紹介
倉橋貴彦
1. ペナルティ法の応用（接触問題に対するペナルティ法の利用）
※部材1が少しずつ降下して、下の複合部材に接触して、埋め込まれていく問題
荷重P
節点1
荷重P
領域（１）
部材(1)
変形
節点1
降下
節点2
節点3
部材(1)
節点2
節点3
部材(2)
節点4
節点4
接触
節点1
部材(2)
節点2
節点3
節点4
領域（２）
x
節点5
u5  0
x
部材(3)
節点5
u5  0
荷重P
部材(1)
部材(2)
x
部材(3)
節点5
u5  0
部材(3)
重み付き残差法の計算式に立ち返って考えてみると、

xb
u x r x dx 
*
xa

ua*

xb
xa
ub*

d 
dux  
u x   EA
dx  0
dx 
dx 
*

 K(e)  ua 
*
    ua
K(e)  ub 
 K (e)

 K(e)
ub*

 f a 


f
 b 
K (i ) 
E(i ) A(i )
L(i )
E(i ) (i)部材のヤング率
A(i ) (i)部材の断面積
L(i ) (i)部材の長さ
荷重P
節点1
制約条件式
部材(1)
節点2
節点3
接触
部材(2)
節点4
x
節点5
u5  0
部材(3)
 u2  u3   E
 *
*
*
u

u

E
2
3

α：ペナルティ係数
①節点2，3における変位量の差は接触しているため，
数値計算の結果として，何倍しても限りなくゼロに近くなるようにする
必要がある。
②重み付き残差法の定式化に伴い、重み関数を節点2，3における
変位量の差と同様の形で定義する。
重ね合わせ後の重み付き残差方程式
領域① 部材(1)
(1)



(e) (1) 
 
x( e ) b
x( e ) a



 K(1)
u(*e) x *(e) r x dx   u1* u2* 

 K(1)



 K(1)   u1 
*
* P 
    u1 u2    0
K(1)  u2 
0 
領域② 部材(2)+(3)
(3)



( e )  ( 2) 
 
x( e ) b
x( e ) a
 K ( 2)


u(*e) x *(e) r x dx   u3* u4* u5*  K( 2)


 0


制約条件式


E*E  u2*  u3*  u2  u3   0

 K ( 2)
K( 2)  K(3)
 K(3)
0  u3 
0
 
*
*
* 
 K(3)  u4   u3 u4 u5  0   0

R 
K(3)  
u5 
 5


各式の総和は重み付き残差法の考えよりゼロになるため、以下のように書くことができる。
(1)



(e) (1) 
 
x( e ) b
x( e ) a
u(*e)
x 
(3)
*
(e) r



x dx 

 ( e )  ( 2) 
 
x( e ) b
x( e ) a

u(*e) x *(e) r x dx   E*E  0


よって、全体系の方程式は、以下のように書くことができる。
u
*
1
u2* u3* u4*
  K (1)

  K (1)

u5*   0

 0

 0


 K (1)
K (1)  

0
0
0

K ( 2)  
 K ( 2)
0
0
0
 K ( 2)
K ( 2)  K (3)
 K (3)
0   u1   P  


0  u2   0  
    

0 u3    0    0

 K (3)  u4   0  
   
K (3)  u5   R5  

重み関数の値は任意であるため、それ以外の部分が零ベクトルになることから、
以下の式が得られる。
 K (1)

 K (1)
 0

 0
 0

 K (1)
K (1)  

0
0
0

K ( 2)  
 K ( 2)
0
0
0
 K ( 2)
K ( 2)  K (3)
 K (3)
0   u1   P 

0  u2   0 
   

0 u3    0 

 K (3)  u4   0 
   
K (3)  u5   R5 
固定条件を考慮することで、最終的に解く方程式は以下のようになる。
 K (1)

 K (1)
 0

 0
 K (1)
K (1)  

0
0

K ( 2)  
 K ( 2)
0
  u1   P 
   
0
 u2    0 
  
 K ( 2)  
u3   0 


K ( 2)  K (3)  
0

u4 
 
例題ペナルティ関数法による接触部における変位の連続性を考慮した変形解析
この例題では，接触点が接合されているモデルと
同じ答えが出る必要がある。
荷重P = 1.0N
節点1
部材(1)
節点2
節点3
接触
x
K (i ) 
 u1  1 

  
u2 3   2
 u  3 
 4   
部材(2)
節点4
節点5
 1 1 0   u1  1
  


 1 2  1 u23   0
  
 0  1 2  
 u4  0
u5  0
よって、以下の式の答えも、上のモデルと同じ解が
得られる必要がある。
部材(3)
E(i ) A(i )
L(i )
1
0
1

 1 1    
 0  1

0
1
 0
0   u1  1
   
0  u2  0
  
 1 u3  0

2  u4  0
 u1  1 
   
u2  2
  
u3  2
u4  3
理想とする答え
K(1)  K(2)  K(3)  1N/m
では、どの程度のαの値を設定する必要があるか？
例題
α=1の場合
荷重P = 1.0N
節点1
部材(1)
節点2
節点3
接触
部材(2)
節点4
x
節点5
K (i ) 
u5  0
部材(3)
E(i ) A(i )
L(i )
K(1)  K(2)  K(3)  1N/m
0   u1  1
 1 1 0

   
 1 2  1 0  u2   0
 0  1 2  1 u3  0


0  1 2  u4  0
 0
 u1  4.00 精度悪い
  

u2  3.00
 

u
2
.
00
3
  

u4  1.00
α=10の場合
1
0
0   u1  1
1

   

1
11

10
0

 u2   0
 0  10 11  1 u3  0


0
 1 2  u4  0
 0
 u1  3.10
  

u2  2.10
 

u
2
.
00
3
  

u4  1.00
α=100の場合
1
0
0   u1  1
1

   

1
101

100
0

 u2   0
 0  100 101  1 u3  0


0
1
2  u4  0
 0
 u1   3.01
  

u2   2.01
 

u3  2.00
u4  1.00 精度良い
αの値は小さい値を設定すると、違った答えが得られる
ことがある。係数行列の数値にも依るがなるべく大きい
値を設定する必要がある。
2. 周期構造に対する数値解析（均質化法への導入）
1
2
3
4
5
(1)
L1
(2)
L2
(3)
L1
(4)
L2
荷重P
1
(1)
L1+L2
(2)
L1+L2
2
3
荷重P
周期構造の場合、ひと固まりの構造を１つの要素として取り扱うことで、最終的に得られる
有限要素方程式の本数を減らすことができる。
⇒では、どのように一つの要素として取り扱いができるようにするか？
例題
巨視的な弾性定数K
1
2
3
K(1)  1N/m K( 2)  2 N/m
P  1N
(1)
個々の要素の弾性定数K
(2)
K (e) 
E(e) A(e)
L(e)
巨視的弾性定数は
結果として、
巨視的変位は
荷重P
まず、巨視的変位を考える。
（周期構造の一つの塊としての変位）
~ P
U ~
K
~ 1 2 2
K

1 2 3
~ P 1 3
U  ~    1.5m
K 2 2
3
また、有限要素法により解くと、
~
P
P
U  u1  u2 

K1 K 2
巨視的変位を分解すると、
右のように書くことができる。
巨視的弾性定数
~ P
P
1
K1K 2
K ~


P
P
1
1
K1  K 2
U


K1 K 2
K1 K 2
 1 1 0   u1   R1 

   

1
3

2

 u2    0 
  
 0  2 2  
u3   1 
同じ答えが
得られる。
 3 2 u2  0

    
 2 2  u3  1
u2   1m 
   

u
1
.
5
m
3

  
例題巨視的弾性定数を考慮した有限要素解析
1
2
3
4
5
(1)
L1
(2)
L2
(3)
L1
(4)
L2
荷重P
K (1)  1N/m K ( 2)  2 N/m
K (3)  1N/m K ( 4)  2 N/m
P  1N
1
(1)
L1+L2
(2)
L1+L2
2
3
荷重P
~
2
K (1)  N/m
3
~
2
K ( 2)  N/m
3
P  1N
通常の有限要素法による計算
均質化法により考えた場合
 2
 3
 2

 3

 0

2
3
4
3
2

3


0 
U1   R1 

2    
  U 2    0 
3    
2  U 3   1 
3 
0
0   u1   R1 
 1 1 0

   

1
3

2
0
0

 u2   0 
 0  2 3  1 0  u3    0 

   
0
0

1
3

2

 u4   0 
0
0
0  2 2  u5   1 

固定条件の式を除くと・・・
固定条件の式を除くと・・・
 4
 3
 2

 3
2
  U  0
3 2 
2 U 3  1

3 
0  u2  0
 3 2 0

   

2
3

1
0

 u3   0
 0  1 3  2 u4  0

   
0
0

2
2

 u5  1
結果として変位は・・・
結果として変位は・・・
U 2  1.5m 
 

U
3
.
0
m

 3 
対応する節点において、
同じ答えが得られる。
u2  1.0m 
  

u3  1.5m 
 

u
2
.
5
m
4
  

u5  3.0m 

Download Report