Sample Scripts from GB Books
GB005:  1次元/2次元過渡問題

ここまでは本質的に1次元の問題を便宜的に2次元のモデルとして扱ってきたわけですが、今回は2次元における熱伝導の問題を正面から扱うことにします。右図中央にあるのは断面が長方形の鉄製の棒で、電磁誘導によって熱せられるものとします。周囲は熱伝導率の低いガラス繊維によって囲まれており、その外側を囲う正方形の箱は一定温度300度Kに保たれているものと仮定します。

以下に示すスクリプトの中では従属変数 temp に対し threshold = 0.1 という指定が行われている点に注意してください。これより小さな温度変動に対しては興味がないことを示すもので、計算時間の短縮に効果があります。

基本形は transient01a.pde と変わりません。
  TITLE
    'Transient Heating of Steel Bar'    { transient01d.pde }

  SELECT
    Errlim = 1e-3

計算効率を考慮し、従属変数 temp に関する0.1以下の変動には関心がないことを宣言します。
  VARIABLES
    temp(threshold = 0.1)     { Absolute tolerance }

全体のエンクロージャのサイズは 20cm x 20cm、鉄棒断面のサイズは 2cm x 2/3cm とします。また熱源のパワー密度を 1e7 に設定します。鉄棒から発せられる単位時間当りの熱量はこれに断面積をかけた値となります。
  DEFINITIONS
    L = 100e-3  L0 = 10e-3
    heat                      { Heat source }
    k  rcp                    { Heat parameters }
    tempi = 300               { Initial and boundary temperatures }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)
    power_d = 1e7                    { Power density }
    power = (2*L0)*(2*L0/3)*power_d  { Power from bar }

温度の初期値は300度Kとします。
  INITIAL VALUES
    temp = tempi

今回の場合、鉄棒部分では右辺は0とはなりません。
  EQUATIONS
    div(fluxd) + rcp*dt(temp) = heat

以下のBOUNDARIESセクションにおいて、Region 1ではドメイン全体としての形状とデフォルトの材質特性を設定、Region 'steel'では鉄棒部分の範囲と固有パラメータの設定を行います。
  BOUNDARIES
    Region 1  heat = 0  k = 0.5  rcp = 2e6  { Glass fiber, default }
      Start(-L,-L)
        Value(temp) = tempi  Line to (L, -L) to (L, L)
                                  to (-L, L) to Close
    Region 'steel'  heat = power_d  k = 45  rcp = 3.5e6
                                            { Stainless steel }
      Start(-L0, -L0/3)  Line to (L0, -L0/3) to (L0, L0/3)
                              to (-L0, L0/3) to Close

  TIME
    from 0 to 1e6

プロットを出力するタイミングは対数スケールを意識して設定します。 等温線図については鉄棒周囲のみを拡大した図も作成します。
  PLOTS
    for t = 100, 300, 1e3, 3e3, 1e4, 3e4, 1e5, 3e5, 1e6
    Contour(temp)
    Contour((temp - tempi)*rcp) Report(power*t)
    Contour(temp) zoom(-5*L0, -5*L0, 10*L0, 10*L0)

  END

ここでは出力結果の一部のみを紹介するに留めます。

(1) Contour(temp) at t = 100
t = 100 の時点では周囲のガラス繊維の温度はほとんどが300度Kのままです。

(2) Contour(temp) at t = 1000
t = 1000 の時点における等温線図です。鉄棒の温度は780度Kに達しています。

(3) Contour(temp) at t = 10000
t = 10000 の時点における等温線図です。エンクロージャに近い部分の等温線の形状は角の丸い正方形に近づいて行きます。鉄棒の温度は1300度Kまで上昇しています。

(4) Contour(temp) at t = 1e5
t = 100000 の時点における等温線図です。

(5) Contour(temp) at t = 1e6
t = 1000000 の時点における等温線図です。(4)と比べてほとんど変化が見られなくなっています。鉄棒部分の最高温度は1430度Kです。

(6) Contour(temp) zoom(-5*L0, -5*L0, 10*L0, 10*L0) at t = 1e6
最終時点における等温線図(5)の中心部分を拡大表示したものです。鉄棒に比較的近い部分で等温線の形状はほぼ円形になっている ことがわかります。

(7) Contour((temp - tempi)*rcp) Report(power*t) at t = 1000
これはドメイン全体でのエンタルピー（熱エネルギー）増分を示す図です。Report文によってプロット下部に鉄棒から発せられた熱量 (power*t) の積分値が示されていますが、t = 1000 あたりまでは良い一致が見られます。

(8) Contour((temp - tempi)*rcp) Report(power*t) at t = 1e6
これに対し t = 1000000 の時点では積分値の乖離が非常に大きくなっています。これは外部境界から熱が失われていることによるものです。境界条件として断熱（Natural(temp) = 0）という条件を設定すれば積分値は一致するようになります。