異材接合体の特異性のオーダの計算

異材接合体に関する
特異性のオーダの計算
倉橋貴彦
応力のつり合い方程式（x-y座標）
  x  xy

0

y
 x

  xy   y  0
 x
y

x 
2
y
2
, y 
式(1)
2
x
2
,  xy
Χ：エアリーの応力関数
（式(2)により表される応力は，
式(1)を満たしていることが分かる。
これにより，式(1)を表す応力は式(2)に示すエアリーの応力関数
（スカラー変数）Χにより与えられることが分かる。）
2

xy
式(2)
応力のつり合い方程式（r-θ座標）
  r 1  r  r   
0
 r  r  
r

 1     r  2 r  0
式(3)
 r 

r
1  1  2 
1 2
1 2 1 2
r 
 2
,    2 2 ,  r  2

2
2
r r r 
r r
r r
r 
式(4)
フックの法則（r-θ座標）
1

 r  E  r   

1

   r 


 
E


1
21   




 r
 r G r
E

式(5)
ひずみの定義式（r-θ座標）
ur

 r  r

ur 1 u




 r
r r 


1 ur u u


 r 
r



r
r

式(6)
ひずみの定義式からur,uθを除去すると・・・
 2
r
2

1  2 r
r
2

2
（極座標系の適合条件式）
2  1  1  2 r
1 



 2 r
r r
r r
r r r 
式(7)
式(4)【エアリーの応力関数により表した応力（つり合い条件式）】を
式(5)【フックの法則】に代入し，
得られた式を式(7)【適合条件式（ひずみ変位関係式）】に代入すると
以下に示す式(8)【重調和方程式 (Biharmonic equation)】が得られる．
2
2 
 
1

1


  0


 r 2 r r r 2  2 


式(8)
4
 2   1  2  2  3
4 2
 2
 2 2 
 4
4
3
r r
r
r r r r
r  2

2 3
r r
3
2

2
4
r r 
2
 r ,   r p 1F  
2
2

1 4
r 
4
4
0
重調和方程式の答えであるΧは，
以下の式により表すことができる．
（重調和方程式に代入することで，解となることを確認できる．）
 r p 1a sin( p  1)  b cos(p  1)  c sin( p  1)  d cos(p  1) 
 r ,   r p 1F  
 r p 1a sin( p  1)  b cos(p  1)  c sin( p  1)  d cos(p  1) 
応力の式にΧを代入すると・・・
1  1  2 
p 1 
F ( )    1F ( )
 r (  rr ) 
 2

r
2
r r r 
  (   ) 
1 2
r 2 r 2
1 2
 r p 1   1F ( )
1 2
 r (  r )  2 2 
 r p 1F ( )
r r
r 
 ij  r p 1
 ij  r p 1
ui  r p
Bogyの特性方程式
この方程式のpが特性根
(特異性のオーダ λ=1-p)
A  2 B  C  2 D  2 E  F  0
2
2
K  p, x   sin 2  px   p 2 sin 2 x
Singular point
A  4 K  p, 1 K  p,  2 
Material 1
E1, ν1
B  2 p 2 sin 2 1 K  p,  2   2 p 2 sin 2  2 K  p, 1 
 
D  2 p sin  sin  p   sin
C  4 p 2 p 2  1 sin 2 1 sin 2  2  K  p, 1   2 
θ1
O
2
θ2
Material 2
E2, ν2
2
1
2
E   D  K  p,  2   K  p, 1 
θ2はマイナスを入れる。
（図の場合θ=-90°）
F  K  p, 1   2 
Dundursパラメータ α、β

2
1  1   2  1
1  1   2  1
, 
1  1   2  1
1  1   2  1

2
 2 sin 2  p1 
Bogyの論文と
角度の取り方が異なるため符号が
逆になっていることに注意
2
1
i 
3  4 i

 i   3  i
1 
i

Ei
, (i  1,2)
2(1  i )
・・・Plane strain
・・・Plane stress
Fortranによる各特性根pに対する特性方程式の値を出力するプログラム
implicit double precision ( a-h, o-z )
c
c Unit of E1, E2, Pa
c Angle theta1, theta2,radian
c
open(11,file='output.dat')
c
E1 = 210.d0 * (10.d0**9)
E2 = 72.d0 * (10.d0**9)
c
anu1 = 0.30d0
anu2 = 0.35d0
C
pi = atan(1.d0)*4.d0
theta1 = 90.d0 * pi / 180.d0
theta2 = -90.d0 * pi / 180.d0
c
c===== Computation of Dundurs Parameter
c
gamma = (E2/(2*(1.d0+2*anu2))) / (E1/(2*(1.d0+2*anu1)))
c
akk1 = ( 3.d0 - anu1 ) / ( 1.d0 + anu1 )
akk2 = ( 3.d0 - anu2 ) / ( 1.d0 + anu2 )
c
alpha = ( gamma * ( akk1 + 1.d0 ) - ( akk2 + 1.d0) )
& / ( gamma * ( akk1 + 1.d0 ) + ( akk2 + 1.d0) )
c
beta = ( gamma * ( akk1 - 1.d0 ) - ( akk2 - 1.d0) )
& / ( gamma * ( akk1 + 1.d0 ) + ( akk2 + 1.d0) )
c
write(*,*) 'alpha=',alpha
write(*,*) 'beta=',beta
c
write(*,*) 'alpha*(alpha-2*beta)', alpha*(alpha-2.d0*beta)
c
c----c
dp = 0.0001
imax = 1/dp
c
c===== Computation of characteristic equation
c
c
do i = 1,imax
c
p = i * dp
c
ak1 = sin(p*theta1)*sin(p*theta1)
& - p*p*sin(theta1)*sin(theta1)
ak2 = sin(p*theta2)*sin(p*theta2)
& - p*p*sin(theta2)*sin(theta2)
ak3 = sin(p*(theta1+theta2))*sin(p*(theta1+theta2))
& - p*p*sin(theta1+theta2)*sin(theta1+theta2)
ak4 = sin(p*(theta1-theta2))*sin(p*(theta1-theta2))
& - p*p*sin(theta1-theta2)*sin(theta1-theta2)
c
aa = 4.d0 * ak1 * ak2
bb = 2.d0 * p * p * sin(theta1)*sin(theta1) * ak2
& + 2.d0 * p * p * sin(theta2)*sin(theta2) * ak1
cc = 4.d0*p*p*(p*p-1.d0)*sin(theta1)
& *sin(theta1)*sin(theta2)*sin(theta2)
& + ak3
dd = 2.d0*p*p*
& ( sin(theta1)*sin(theta1) * sin(p*theta2)*sin(p*theta2)
& - sin(theta2)*sin(theta2) * sin(p*theta1)*sin(p*theta1) )
ee = - dd + ak2 - ak1
ff = ak4
c
afun = aa*(beta**2) + 2.d0*bb*alpha*beta + cc*(alpha**2)
&
+ 2.d0*dd*beta + 2.d0*ee*alpha + ff
c
write(11,100) p, afun
100 format(2f15.10)
c
end do
c
stop
end
Bogyの特性方程式の計算結果
特性根の刻みを0.0001として0～1まで計算した結果
拡大
p=0.9287の時
特性方程式の値は
概ね0である．
E1 = 210GPa
E2 = 72GPa
ν1 = 0.30
ν2 = 0.35
θ1=90 deg.
θ2 = -90 deg
特異性のオーダλ
λ=1-p=1-0.9287=0.0713
数値解法により求める場合は、係数を入れた特性方程式に対して，
二分法等による反復法により特性根pを求める。また、虚部が現れる
問題もあるため、その際は別途、特性方程式の計算プログラムを構
築する。

Download Report