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

オブジェクトが熱的に定常的な状態にあるということは、任意の体積要素から出て行く熱量とその中で生成される熱量とが等しいことを意味します。 これを数学的に表現すると次のようになります。

ここに f は

によって定義され、熱流束密度（heat flux density、単位面積当りに伝達される熱量）を表します。ただし λ は熱伝導率（thermal conductivity）、T は絶対温度、h は単位体積当りに発生する熱量を意味します。この場合の熱源としては材質中を流れる電流、電磁波、化学反応等が考えられます。(1), (2)より支配方程式は次のような2階の偏微分方程式となります。

GB003で扱った電導の問題とは次のような対応関係があります。

なお、以下に示す用例においてはSI単位系（m, K, W）を使用されている点にご注意ください。

右の図は中に2本の管を含む直径 1m の円筒状のパイプの断面図です。左右の管の中にはそれぞれ50度C、80度Cの温水が流れているとしたとき、断熱性素材からなる灰色を施した部分における熱伝導の様子をFlexPDEを使って調べてみることにします。ただしパイプは0度Cの土壌中に埋設されているものとします。
以下では温度分布や熱流束密度ベクトル場のプロットと共に、熱流の出入りの帳尻が合っているかどうかも検証してみます。なお、(3)では温度 T に関する2階の方程式が規定されているわけですが、スクリプト作成時には T や t という変数名を使わないようにしてください。t は時間変数を表す予約語です。

Note: FlexPDEでは大文字、小文字は区別されません。

まずタイトルを設定します。
  TITLE
    'Two Insulated Tubes'    { hconduction01a.pde }

次に演算精度に関するセレクタをセットします。デフォルトは 0.002 なのですが、ここでは精度を多少高めに設定します。
  SELECT
    Errlim = 3e-4

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

偏微分方程式の定義に先立ち、パラメータ類を定義します。これらは境界の定義、及び境界条件の設定に際して使用されます。解析対象のドメイン上に熱源がないため上記数式中における h は0、また素材の熱伝導率は 0.03 と設定している点に注意してください。
  DEFINITIONS
    r0 = 0.1  d = 0.15  r1 = 0.5
    k = 0.03              { Thermal conductivity }
    heat = 0              { Power density }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)  { Heat flux density }
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)
    f_angle = sign(fluxd_y)*arccos(fluxd_x/fluxd_m)/PI*180

方程式を定義します。従属変数 temp の数式表現になっていませんが、DEFINITIONSセクションの情報に基づき、temp に関する2階のPDEが自動的に導出されます。
  EQUATIONS
    div(fluxd) = heat     { 2nd order PDE in temp }

次に境界の形状と境界条件を設定します。外界との境界条件、及び2本の管との境界条件、共にDirichlet型であるためValue文を使用します。
  BOUNDARIES
    Region 1
      Start 'outer' (0, -r1)  Value(temp) = 273     { Frozen soil }
        Arc(Center = 0,0) Angle = 360 to Close
      Start 'left' (-d - r0, 0)  Value(temp) = 323 
                                 { Exclude left hot water tube }
        Arc(Center = -d,0) Angle = 360 to Close
      Start 'right' (d - r0, 0)  Value(temp) = 353 
                                 { Exclude right hot water tube }
        Arc(Center = d,0) Angle = 360 to Close

最後に出力すべき情報を指定します。
  PLOTS
    Grid(x, y)
    Contour(temp)    Surface(temp)
    Vector(fluxd) norm
    Contour(f_angle)
    Contour(fluxd_x)    Contour(fluxd_y)    Contour(fluxd_m)
    Elevation(normal(fluxd)) on 'outer'
    Elevation(-normal(fluxd)) on 'left'    { Outwards from tube }
    Elevation(-normal(fluxd)) on 'right'

  END

(1) Grid(x, y)
FlexPDEによって設定されたグリッドの形状を示しています。

(2) Contour(temp)
解析対象領域（ドメイン）上での等温線図は次のようになります。

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

(4) Vector(fluxd) norm
熱流束密度 f のベクトルプロットを示したものです。等温線図とは直交したものと なっています。左の管は熱を発散すると同時に吸収もしている点にご注意ください。

(5) Contour(f_angle)
熱流束密度ベクトル f の偏角に関する等高線図をプロットしたものです。

(6) Contour(fluxd_x)
熱流束密度ベクトル f のx成分の値に関する等高線図です。

(7) Contour(fluxd_y)
熱流束密度ベクトル f のy成分の値に関する等高線図です。

(8) Contour(fluxd_m)
熱流束密度ベクトル f の絶対値に関する等高線図です。

次に境界上における熱流束密度ベクトル f の 法線成分 f_n の値をElevationコマンドによりプロットしてみます。グラフの形状もさることながら、プロット下部に出力されるその積分値に意味があります。外界に流出して行く熱量と温水チューブから流入して来る熱量とがバランスしていなくてはならないからです。

(9) Elevation(normal(fluxd)) on 'outer'
外部境界上での熱流束密度ベクトル f の外向き法線成分の値をプロットしたのが次の図です。境界上での積分値は 12.225 という値になっています。

(10) Elevation(-normal(fluxd)) on 'left'
次に左側のチューブ上で熱流束密度ベクトル f の 外向き法線成分の値をプロットしてみます。ここで言っている“外向き”というのは解析対象のドメインから見て外向きということなので、normalという演算子に対し負号が付いている点に注意してください。境界上での積分値は -2.806 という値になっています。

(11) Elevation(-normal(fluxd)) on 'right'
同様にして右側のチューブ上で熱流束密度ベクトル f の 外向き法線成分の値をプロットしてみます。こちらの場合はすべて流入となるため、関数値はすべてマイナスとなっています。境界上での積分値は -9.415 であり、(10)の -2.806 と合わせると流入分の合計値は 12.221 となります。(9)で求めた外界への流出分 12.225 と比べたときの誤差は非常に小さい（0.04%）と言えます。