Sample Scripts from GB Books
GB006:  非線形熱伝導

今度は薄い鉄製のフォイルを使って同様の熱伝導について考察します。右図のように今回は円形の穴は設けていませんが、薄い形状であるため、表面、裏面からの熱損失を無視できず、その影響を含んだ方程式を設定する必要があります。なお熱伝導率 λ は温度に依存しないものとします。

方程式上は h の項が0でなくなります。今の場合 h < 0 となり、単位体積当りの熱損失を表現することになります。厳密なことを言えば、熱損失は表面と裏面から起るため、3次元オブジェクトとしてのモデル化を行い、z軸方向の法線微分について考えるべきなのでしょうが、ここでは2次元のモデルで代用することにします。すなわちフォイルは十分に薄く、z軸方向の温度差はないものとします。換言するなら、放射は表面から起っているのではなく、体積全体から起っているものと考えるわけです。

単位面積当りに放射される熱エネルギーを f とすると、

で与えられます。ここに ε は放射率（emissivity）、σ はStefan-Boltzmann定数を意味します。フォイルは単位時間当りこの率で熱を失うわけですが、逆に周囲（300度Kと仮定）からの赤外線放射による熱吸収という側面も考慮する必要があります。従って以下のスクリプト中では

という設定を行っています。2を乗じているのは表面と裏面双方を配慮してのことです。dはフォイルの厚さを意味し、それで割ることによって単位体積当りの熱損失が表現されます。

基本形は htransfer01a.pde と変わりません。
  TITLE
    'Foil Emitting IR Radiation'    { htransfer01b.pde }

  SELECT
    Errlim = 1e-4

  VARIABLES
    temp                { Temperature }

放射率 ε としては0.3を、Stefan-Boltzmann定数 σ としては5.67e-8という値を使用します。また偏微分方程式の左辺と右辺の値を別個に比較するために、equ1, equ2というパラメータを定義しています。
  DEFINITIONS
    Lx = 0.6  Ly = 0.4  d = 1e-3
    heat = -2*0.3*5.67e-8*(temp^4 - 300^4)/d
                        { Heat source per unit volume }
    k = 45              { Thermal conductivity of steel }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)
    equ1 = div(fluxd)
    equ2 = heat

従属変数 temp の初期値を設定します。非線形の問題の場合には適切な初期値を設定することが望ましいわけですが、指定が省略された場合には0が仮定されます。
  INITIAL VALUES
    temp = 600

  EQUATIONS
    div(fluxd) = heat

  BOUNDARIES
    Region 1
      Start(0,0)
        Natural(temp) = 0   Line to (Lx, 0)  { Insulated }
        Value(temp) = 300   Line to (Lx, Ly)
        Natural(temp) = 0   Line to (0, Ly)  { Insulated }
        Value(temp) = 1000  Line to Close

  PLOTS
    Grid(x, y)
    Contour(temp)  Surface(temp)
    Elevation(temp) from (0, Ly/2) to (Lx, Ly/2)
    Elevation(heat) from (0, Ly/2) to (Lx, Ly/2)
    Vector(fluxd) norm  Contour(fluxd_m)
    Contour(fluxd_x)  Contour(fluxd_y)
    Elevation(equ1, equ2) from (0, Ly/2) to (Lx, Ly/2)

  END

(1) Grid(x, y)
メッシュ再構成が2度ほど行われた結果、左端の領域でのメッシュ密度が高くなっています。

(2) Contour(temp)
等温線をプロットしたものです。曲面の形状は次のsurfaceプロットの方がわかりやすいでしょう。

(3) Surface(temp)
当然ですが熱損失があるため、温度が高い部分は左側境界に近い部分に限られます。

(4) Elevation(temp) from (0, Ly/2) to (Lx, Ly/2)
y = 0.2 という中心線に沿って温度の値をプロットしたものです。

(5) Elevation(heat) from (0, Ly/2) to (Lx, Ly/2)
y = 0.2 という中心線に沿って熱損失の大きさをプロットしたものです。

(6) Vector(fluxd) norm
熱流束密度 f のベクトルプロットを示したものです。normを指定しているため矢印の長さは均一ですが、ベクトル長は色で識別できます。

(7) Contour(fluxd_m)
熱流束密度 f の絶対値に関する等高線図です。温度勾配のきつい左端の近傍で大きな値となっています。

(8) Contour(fluxd_x)
熱流束密度 f のx成分の値 fx に関する等高線図です。 水平方向の流れであるため実質的に(7)の図と変わりません。

(9) Contour(fluxd_y)
熱流束密度 f のy成分の値 fy に関する等高線図です。 水平方向の流れであるためほぼ一様に0に近い値になっています。

(10) Elevation(equ1, equ2) from (0, Ly/2) to (Lx, Ly/2)
偏微分方程式の左辺 div(fluxd) の値（equ1）と右辺 heat の値（equ2）を y = 0.2 の線し沿ってelevationプロットの形で示したものです。概ね良い一致を示していると言えます。左辺のequ1は2階の導関数に基づくものであるため、 区間別に定義された定数値として表現されています。