Sample Scripts from GB Books
GB010: 3次元の熱伝導

ここで紹介するスクリプトはGunnar Backstrom氏の承諾のもと、書籍 “Simple Fields of Physics by Finite Element Analysis” に記されている多数のFlexPDE適用事例の中からその一部を紹介するものです。

PDF版 (349KB)

FlexPDEで3次元の場を扱う場合には、x-y平面に定義されたbase planeをz軸方向に押出す（extrusion）という考え方が基本になります。

1. 平衡状態

ここでは右図のような鉄の板に銅製の円柱を溶接した構造体（単位はm）を考え、その円柱部上面に1e6の熱流束を与えたときの平衡状態における熱伝導の様子をFlexPDE を用いて解析します。ただし境界条件としては次の平面図、立面図に示したようなものを想定します。

なお、支配方程式に関する背景情報については GB004 を参照ください。

平面図	立面図

正面は熱的に絶縁されており、他は一定温度300度Kに保つものとします。この条件はz軸方向にextrudeされた側壁に対しても適用される点に注意してください。	底面は熱的に絶縁状態、一方円柱部上面からは1e6の熱流束が与えられるものとします。

1.1 Problem descriptor [ 3dfields2a.pde ]

まずタイトルを設定します。
TITLE
'Steady Heat Conduction' { 3dfields2a.pde }

次に演算精度に関するセレクタをセットします。
SELECT
Errlim = 1e-3

次に座標系が3次元直交座標系であることを明示します（デフォルトはCartesian2）。
COORDINATES
Cartesian3

従属変数を定義します。
VARIABLES
temp

関連するパラメータや数式を定義します。r0 は円柱の半径、中心の位置は (L/4, L/4)、z0 は鉄板の厚さを意味しています。
DEFINITIONS
    L = 0.1 r0 = 0.01 z0 = 0.02
    k heat = 0       { Thermal conductivity and power density }
    in_fluxd = 1e6    { Input heat flux density at top }
    fluxd_x = -k*dx(temp) fluxd_y = -k*dy(temp) fluxd_z = -k*dz(temp)
    fluxd = vector(fluxd_x, fluxd_y, fluxd_z)
    fluxdm = magnitude(fluxd)

3次元の熱伝導方程式を定義します。この偏微分方程式の誘導については GB004 を参照ください。
EQUATIONS
    div(fluxd) = heat

z軸方向へのextrusionを定義します。下から順に境界面（Surface）とそれらにはさまれた層（Layer）を規定して行きます。
EXTRUSION                  { Parallel surfaces }
    Surface 'bottom' z = 0
    Layer 'iron'
    Surface 'middle' z = z0 { Interface }
    Layer 'copper'
    Surface 'top' z = L

BOUNDARIESセクションにおいてbase plane上でのリージョン構成と境界条件を定義します（上記平面図参照）。これらbase plane上での形状定義や境界条件は基本的にz軸方向にextrudeされて行きます。ただし上面と下面についてはSurface文にて境界条件を明示します。なおリージョン'domain'の定義中、レイヤ'copper'の部分についてはvoidという修飾子を用いている点に注意してください。リージョン'cylinder'内では鉄、銅双方の熱伝導率が定義されます。
BOUNDARIES
    Surface 'bottom' Natural(temp) = 0       { Insulated }
    Surface 'top' Natural(temp) = -in_fluxd { Input flux density }

    Region 'domain'                          { Full solution domain }
      Layer 'iron' k = 82                    { Fe }
      Layer 'copper' void
      Start 'outer' (0, 0)
        Natural(temp) = 0 Line to (L, 0)
        Value(temp) = 300 Line to (L, L) to (0, L) to Close

    Region 'cylinder'
      Layer 'iron' k = 82                    { Fe }
      Layer 'copper' k = 400                 { Redefine void to Cu }
        Start (L/4+r0, L/4) Arc(Center = L/4, L/4)
                            Angle = 360 to Close

最後に出力すべき情報を規定します。
PLOTS
    Grid(x, y, z)
    Contour(temp) on x = y
    Contour(temp) painted on x = y
    Contour(fluxdm) painted on x = y
    Vector(fluxd) norm on x = y
    Contour(temp) on y = 0
    Contour(temp) on y = L/4
    Contour(temp) on x = L/4
    Contour(temp) on z = z0 on 'iron'
    Contour(temp) on z = 0

END

1.2 実行結果

(1) Grid(x, y, z)
FlexPDEによって生成された3次元のメッシュ構成を示しています。

(2) Contour(temp) on x = y
対角面 x = y 上における等温線図を示したものです。銅製の円柱部は頂部が約600度Kで、下に向かってほぼリニアに温度は下がって行きます。また鉄製の基盤部両端は300度Kに保たれていることがわかります。

(3) Contour(temp) painted on x = y
プロット(2)と同一内容ですが、塗りつぶしを施した等温線図です。

(4) Contour(fluxdm) painted on x = y
熱流束密度ベクトル f の絶対値についての等高線図です。円柱頂部から 1e6 の熱流束が流れ込んでいることを反映して、円柱内での |f| はほぼ一様に 1e6 という値になっています。なお、最大値は鉄板との接合部で発生しています。