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

最初のケースと基本的には同一の問題設定ですが、ここでは円柱頂部に定量の熱流を与え始めた時点からの時間変化をFlexPDEを使って考察してみます。対角面 x = y 上における等温線図を0.1秒後から1000秒後までの範囲でプロットしてみます。十分な時間が経過した後の等温線図は平衡状態のそれに漸近して行くはずです。なお、時間変化の解析には数分の計算時間を要します。

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

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

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

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

関連するパラメータや数式を定義します。
  DEFINITIONS
    L = 0.1  r0 = 0.01  z0 = 0.02
    k  rcp  heat = 0                { Thermal parameters }
    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)
    f_angle = sign(fluxd_y)*arccos(fluxd_x/fluxdm)/PI*180

初期状態を定義します。当初温度は一律300度Kの状態にあるとします。
  INITIAL VALUES
    temp = 300

時間依存型の熱伝導方程式は

で与えられますが、背景情報については GB005 を参照ください。なお、ρ（密度）、c_p（比熱）の値は材質ごとに異なる点に注意する必要があります。
  EQUATIONS
    div(fluxd) + rcp*dt(temp) = heat

z軸方向へのextrusionを定義します。
  EXTRUSION                         { Parallel surfaces }
    Surface 'bottom' z = 0
    Layer 'iron'
    Surface 'middle' z = z0         { Interface }
    Layer 'copper'
    Surface 'top' z = L

BOUNDARIESセクションにおいてbase plane上でのリージョン構成と境界条件を定義します。
  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 rcp = 7.87e3*449     { 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 rcp = 7.87e3*449     { Fe }
      Layer 'copper' k = 400 rcp = 8.96e3*385  { Redefine void to Cu }
      Start (L/4+r0, L/4) Arc(Center = L/4, L/4) Angle = 360 to Close

時間依存型問題に固有の指定項ですが、TIME セクションを用いて解析対象の時間帯を指定します。
  TIME
    0 to 1e3

最後に出力すべき情報を指定します。ここでは対角面 x = y 上における等温線図を0.1秒後から1000秒後までの範囲でプロットします。
  PLOTS
    for t = 0.1, 5, 10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 100,
            200, 500, 1000
    Contour(temp) painted on x = y

  END

上記スクリプトをそのままFlexPDEに実行させるとすべてのプロットが一気に生成され、t = 1000 に対応したプロットが表示された状態で処理が終了します。最初からすべてのプロットを確認するには「File」メニュー：「View File」と操作し、該当する.pg5ファイルを選択してください。「View」メニュー：「Next」あるいは「Back」と操作することにより手動で前後に進めることができます。一方、「View」メニュー：「Movie」と操作すれば画面切替えは一定時間間隔で自動的に行われるようになります。

合計17枚のプロットが生成されますが、ここでは t = 0.1, 60, 1000 の3ケースについてのみ結果を示すことにします。

(1) Contour(temp) painted on x = y at t = 0.1
対角面 x = y 上での等温線図です。まだ時間が0.1秒しか経過していないため、円柱頂部を除けば温度は一律300度Kです。

(2) Contour(temp) painted on x = y at t = 60
t = 60における等温線図です。円柱頂部の温度はまだ530度Kであり、平衡状態には達していません。

(3) Contour(temp) painted on x = y at t = 1000
t = 1000における等温線図です。平衡状態の場合のプロット(3)と対比してみてください。