Sample Scripts from GB Books
GB004:  2次元熱伝導

今度は 断面が半円形のステンレス製鋼材（長さは十分に長いとします）に電流を流し一様に加熱します。周囲は右図のように一辺が0.2mの花崗岩で囲まれているとし、その外側は氷と水によって0度Cに保たれている状態を考えます。このモデルの場合には偏微分方程式中に熱源の項 h(x, y) を含めて考える必要があります。単位体積当り花崗岩中では h = 0、鋼材中では h = 1e6 であるとします。また熱伝導率は花崗岩が k = 3.5、ステンレス鋼が k = 45 という値を仮定します。
スクリプトのベースとしては hconduction01a.pde を用いますが、鋼材部分をexcludeするのではなく別個のリージョンとして扱う点が異なります。

基本形は hconduction01a.pde と大きく変わりませんが、ドメイン全体を規定するリージョン1の他に鋼材部分を規定するリージョン'steel'が定義されている点に注意してください。リージョン1中で設定したパラメータ値はドメイン全体に適用されますが、リージョン'steel'中でもパラメータ値の指定を行っているため、該当リージョン中では後者の値が使用されることになります。
  TITLE
    'Heated Semi-Circular Rod'    { hconduction01b.pde }

  SELECT
    Errlim = 3e-4

  VARIABLES
    temp                  { Temperature }

  DEFINITIONS
    Lx = 0.1  Ly = 0.1  r0 = 0.05
    k                     { Thermal conductivity }
    heat                  { Power per unit volume }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)  { Heat flux density }
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)

  EQUATIONS
    div(fluxd) = heat     { 2nd order PDE in temp }

  BOUNDARIES
    Region 1  k = 3.5  heat = 0             { Granite, default }
      Start 'outer' (-Lx, -Ly)  Value(temp) = 273
        Line to (Lx, -Ly) to (Lx, Ly) to (-Lx, Ly) to Close
    Region 'steel'  k = 45  heat = 1.0e6    { Steel rod }
      Start 'rod' (-r0, 0)
        Line to (r0, 0)
        Arc to (0, r0) to Close

  PLOTS
    Grid(x, y)
    Contour(temp)    Surface(temp)
    Contour(temp) painted on 'steel'
    Vector(fluxd) norm
    Contour(fluxd_m) painted
    Elevation(normal(fluxd)) on 'rod'
    Elevation(normal(fluxd)) on 'outer'
      Report(PI*r0^2/2*1e6) as 'Heating power'

  END

(1) Grid(x, y)
図面下部に“Grid#3”と表示されていることからわかるように、所定の精度を得るために2回のグリッド再構成が行われています。

(2) Contour(temp)
ドメイン全体としての等温線図は次のようになります。鋼材部分のみのより詳しい等温線図については(4)を参照ください。

(3) Surface(temp)
温度 T(x, y) の曲面の形状をプロットしたものです。

(4) Contour(temp) painted on 'steel'
鋼材部分のみのより詳しい等温線図です。

(5) Vector(fluxd) norm
熱流束密度 f のベクトルプロットを示したものです。温度の最高点から周囲に流れが広がっていく様子が示されています。

(6) Contour(fluxd_m) painted
熱流束密度ベクトル f の絶対値に関する等高線図です。鋼材の角の部分にピークがありそうなことはこの図からも見て取れますが、次のelevationプロットの方がわかりやすいでしょう。

(7) Elevation(normal(fluxd)) on 'rod'
鋼材外周上での熱流束密度ベクトル f の外向き法線成分の値をプロットしたのが次の図です。 角の部分にピークが現れています。また境界上での積分値は 3934 という値になっています。

(8) Elevation(normal(fluxd)) on 'outer' Report(PI*r0^2/2*1e6) as 'Heating power'
外部境界上での熱流束密度ベクトル f の外向き法線成分の値をプロットしたのが次の図です。境界上での積分値は 3920 という値になっています。(7)で求めた値とは多少乖離が見られます。これは鋼材の角の部分で関数値がピークとなっているため、積分計算の精度が落ちているためです。ちなみにこのelevationプロットでは

という値をレポート出力させています。これは半円全体で生み出される単位時間当りの熱量の合計値であり、その値 3927 が上記積分の理論値となる点にご注意ください。