Sample Script

〔10-1〕
3D バイメタルの熱応力・変形解析、及び温度分布解析 【3d_bimetal.pde】

1.概要

より大きな鉄板(メッシュ図の青色)に接合された、小さなアルミニウム角棒(黄色)の問題を考えます。この結合体に対して、底面の温度を100℃で固定し、側面と上面は対流によって冷却します。3Dの温度分布と、関連した変形、及び応力を解析します。 結合体の全ての表面は、拘束条件を与えていません。 温度分布によって自由に変形します。移動と回転をキャンセルするための積分拘束は使用しません。 3Dでの積分拘束は非常に高くつきますので、積分拘束は2Dサンプルで扱います。 その代わり、このケースでは FlexPDE が見つけた解で、作図の前に移動と回転の平均値を除去します。

2.メッシュ図

3d_bimetal-メッシュ

3.解析結果

図中の {Fig.A},{Fig.B},{Fig.C} は、4項で対応するスクリプトを示します。

3D bimetal 変形図

3d bimetal 温度分布図

3d bimetal 応力

4.スクリプト

下記のスクリプトをマウスでコピーし、FlexPDE エディット・ウィンドウに貼り付けて実行する際には、日本語のコメントを除去して下さい。 そのままですと、コンパイル・エラーが発生する場合があります。

{ 3D_BIMETAL.PDE
 3Dバイメタル
  This problem considers a small block of aluminum bonded to a larger block
  of iron.  The assembly is held at a fixed temperature at the bottom, and
  is convectively cooled on the sides and top.  We solve for the 3D temperature
  distribution, and the associated deformation and stress.

  All faces of the assembly are unconstrained, allowing it to grow as the
  temperature distribution demands.  We do not use an integral constraint
  to cancel translation and rotation, as we have done in 2D samples,
  because in 3D this is very expensive.  Instead, we let FlexPDE find a solution,
  and then remove the mean translation and rotation before plotting.
}

title 'Bimetal Part'

coordinates
    cartesian3

select
    painted             { show color-filled contours 等高線図を色分けで塗りつぶします。}

variables
    Tp                  { temperature  温度 }
    U                   { X displacement  X方向変位 }
    V                   { Y displacement  Y方向変位 }
    W                   { Z displacement  Z方向変位 }

definitions
    long = 1
    wide = 0.3
    high = 1
    tabx = 0.2
    taby = 0.4

    K                  { thermal conductivity  熱伝達率 }
    E                  { Youngs modulus  ヤング率 }
    alpha              { expansion coefficient  熱膨張係数 }
    nu                 { Poisson's Ratio  ポアソン比 }

    Q = 0              { Thermal source  熱源 }
    Ta = 0.            { define the ambient thermal sink temperature  周囲の熱シンク温度 }

    { define the constitutive relations  構成関係の定義 }
    G = E/((1+nu)*(1-2*nu))
    C11 = G*(1-nu)
    C12 = G*nu
    C13 = G*nu
    C22 = G*(1-nu)
    C23 = G*nu
    C33 = G*(1-nu)
    C44 = G*(1-2*nu)/2
    b = G*alpha*(1+nu)

    { Strains  歪み }
    ex = dx(U)
    ey = dy(V)
    ez = dz(W)
    gxy = dy(U) + dx(V)
    gyz = dz(V) + dy(W)
    gzx = dx(W) + dz(U)

    { Stresses  応力 }
    Sx  =  C11*ex + C12*ey + C13*ez - b*Tp
    Sy  =  C12*ex + C22*ey + C23*ez - b*Tp
    Sz  =  C13*ex + C23*ey + C33*ez - b*Tp
    Txy =  C44*gxy
    Tyz =  C44*gyz
    Tzx =  C44*gzx

    { find mean translation and rotation  移動と回転の平均値 }
    Vol = Integral(1)
    Tx = integral(U)/Vol                    { X-motion }
    Ty = integral(V)/Vol                    { Y-motion }
    Tz = integral(W)/Vol                    { Z-motion }
    Rz = 0.5*integral(dx(V) - dy(U))/Vol    { Z-rotation }
    Rx = 0.5*integral(dy(W) - dz(V))/Vol    { X-rotation }
    Ry = 0.5*integral(dz(U) - dx(W))/Vol    { Y-rotation }

    { displacements with translation and rotation removed  移動と回転を除去した変位量 }
    { This is necessary only if all boundaries are free  全ての境界条件がフリーの時だけに必要です }
    Up = U - Tx + Rz*y - Ry*z
    Vp = V - Ty + Rx*z - Rz*x
    Wp = W - Tz + Ry*x - Rx*y

    { scaling factors for displacement plots 変形図を作図するための縮小/拡大率 }
    Mx = 0.2*globalmax(magnitude(y,z))/globalmax(magnitude(Vp,Wp))
    My = 0.2*globalmax(magnitude(x,z))/globalmax(magnitude(Up,Wp))
    Mz = 0.2*globalmax(magnitude(x,y))/globalmax(magnitude(Up,Vp))
    Mt = 0.4*globalmax(magnitude(x,y,z))/globalmax(magnitude(Up,Vp,Wp))

initial values
    Tp = 5.
    U = 1.e-5
    V = 1.e-5
    W = 1.e-5

equations
    Tp: div(k*grad(Tp)) + Q = 0.         { the heat equation  熱方程式 }
    U:  dx(Sx) + dy(Txy) + dz(Tzx) = 0   { the U-displacement equation  X方向変位の方程式 }
    V:  dx(Txy) + dy(Sy) + dz(Tyz) = 0   { the V-displacement equation }
    W:  dx(Tzx) + dy(Tyz) + dz(Sz) = 0   { the W-displacement equation }

extrusion z = 0,long

boundaries
    surface 1 value(Tp)=100             { fixed temp bottom  底部の固定温度 }
    surface 2 natural(Tp)=0.01*(Ta-Tp)  { poor convective cooling top  上面の対流冷却 }

    Region 1    { Iron  鉄 }
       K = 0.11
       E = 20e11
       nu  =0.28
       alpha =  1.7e-6
       start(0,0)
         natural(Tp) = 0.1*(Ta-Tp)   { better convective cooling on vertical sides 垂直側面の対流冷却 }
         line to (wide,0)
           to (wide,(high-taby)/2)
           to (wide+tabx,(high-taby)/2)
           to (wide+tabx,(high+taby)/2)
           to (wide,(high+taby)/2)
           to (wide,high)
           to (0,high)
           to close

     Region 2   { Aluminum アルミニウム }
       K = 0.5
       E = 6e11
       nu  =0.25
       alpha =  2*(2.6e-6)              { Exaggerate expansion 誇張した熱膨張係数 }
       start(wide,(high-taby)/2)
         line to (wide+tabx,(high-taby)/2)
           to (wide+tabx,(high+taby)/2)
           to (wide,(high+taby)/2)
           to close

monitors
    contour(Tp) on y=high/2  as "Temperature"
    contour(Up) on y=high/2 as "X-displacement"
    contour(Vp) on x=4*wide/5 as "Y-displacement"
    contour(Wp) on y=high/2 as "Z-displacement"
    grid(x+My*Up,z+My*Wp) on y=high/2 as "XZ Shape"
    grid(y+Mx*Vp,z+Mx*Wp) on x=wide/2 as "YZ Shape"
    grid(x+Mz*Up,y+Mz*Vp) on z=long/4 as "XY Shape"
    grid(x+Mt*Up,y+Mt*Vp,z+Mt*Wp)  as "Shape"

plots
	cdf(Tp)
    contour(Tp) on y=high/2  as "XZ Temperature"
    contour(Up) on y=high/2 as "X-displacement"
    contour(Vp) on x=4*wide/5 as "Y-displacement"
    contour(Wp) on y=high/2 as "Z-displacement"
    grid(x+My*Up,z+My*Wp) on y=high/2 as "XZ Shape"
    grid(y+Mx*Vp,z+Mx*Wp) on x=4*wide/5 as "YZ Shape"
    grid(x+Mz*Up,y+Mz*Vp) on z=long/4 as "XY Shape"
    grid(x+Mt*Up,y+Mt*Vp,z+Mt*Wp)  as "Shape"		{ Fig. A を出力 }
    contour(Sx) on y=high/2 as "X-stress"		{ Fig. B を出力 }
    contour(Sy) on y=high/2 as "Y-stress"
    contour(Sz) on y=high/2 as "Z-stress"		{ Fig. C を出力 }
    contour(Txy) on y=high/2 as "XY Shear stress"
    contour(Tyz) on y=high/2 as "YZ Shear stress"
    contour(Tzx) on y=high/2 as "ZX Shear stress"

end