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

ここで紹介するスクリプトはGunnar Backstrom氏の承諾のもと、書籍 “Simple Fields of Physics by Finite Element Analysis” に記されている多数のFlexPDE適用事例 の中からその一部を紹介するものです。

PDF版 (428KB)

熱伝導方程式でも例えば

  

(1)

のように線形で同次型(右辺が0)の場合には、重ね合せの原理が成り立ちます。すなわち T1, T2 が方程式(1)の解であるとしたとき、T = T1 + T2 も(1)の解となります。
しかし熱伝導率 λ が温度 T に依存するような場合には方程式は T に関し線形とは言えなくなるため、重ね合せの原理は適用できなくなります。このため解析的手法で解を求めることはほぼ絶望的となります。しかしFlexPDEの場合にはこのような非線形の問題であっても取扱いが可能です。

1. 温度依存の熱伝導率

右の図は中央部に円形の穴のあいた酸化マグネシウム(MgO)製の板を示しています。この材料の熱伝導率は T に反比例することが知られています。左端は2300度Kに、右端は300度Kに保たれているとしたときの温度分布をFlexPDEを用いて解析してみます。ただし、その他の境界は熱的に絶縁されているものとします。また、板は十分に厚く、赤外線放射による表面からの熱のロスは無視しうるものとします。

1.1 Problem descriptor [ htransfer01a.pde ]

まずタイトルを設定します。
  TITLE
    'MgO Plate with a Hole, k(T)'    { htransfer01a.pde }

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


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


偏微分方程式の定義に先立ち、パラメータ類を定義します。熱伝導率 k が温度 temp の関数となっているため、方程式が非線形となる点にご注意ください。
  DEFINITIONS
    Lx = 0.6  Ly = 0.4  r0 = 0.1
    heat = 0            { Heat source }
    k = 15000/temp      { Thermal conductivity k(T) }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)


時間依存型問題の場合、INITIAL VALUESセクションは初期条件の設定のために使用されますが、定常型問題の場合には解探索の初期値という意味合いを持ちます。非線形問題の場合にはシステムをうまく解に誘導する意味において、あるいは計算時間を短縮する意味において初期値の設定は有用と言えます。ここでは左右の境界値の中間を取って T(x, y) = 1300 という初期状態を設定します。
  INITIAL VALUES
   
temp = 1300

次に方程式を

の形で設定します。f 熱流束密度(heat flux density)を表すベクトルで

で定義されます。
  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) = 2300  Line to Close
      Start(Lx/2 - r0, Ly/2)                 { Exclude hole }
        Natural(temp) = 0
          Arc(Center = Lx/2, Ly/2) Angle = 360 to Close


最後に出力すべき情報を規定します。
  PLOTS
   
Grid(x, y)
    Contour(temp)
    Contour(k) painted
    Vector(fluxd) norm
    Contour(fluxd_x)  Contour(fluxd_y)  Contour(fluxd_m)
    Elevation(temp) from (0, 0) to (Lx, 0)

  END

1.2 実行結果

(1) Grid(x, y)
FlexPDEによって生成されたメッシュ構成を示しています。メッシュ再構成は1回だけで済んでいます。

(2) Contour(temp)
等温線のプロットです。熱抵抗率(thermal resistivity, 1/λ)が大きなプレート左端、及び穴によってくびれている部分で等温線の密度が高くなっています。

(3) Contour(k) painted
熱伝導率 λ (スクリプト上は k)の値の等高線図です。温度の高い左側では小さな値、温度の低い右側では大きな値となっています。

(4) Vector(fluxd) norm
熱流束密度 f のベクトルプロットを示したものです。等温線図とは直交したものと なっています。

(5) Contour(fluxd_x)
熱流束密度 f のx成分 fx の値に関する等高線図です。等温線図(2)が左右非対称であるのに対し、fx のプロットは左右対称となっています。これが非対称であればエネルギー保存則に反することになります。

(6) Contour(fluxd_y)
熱流束密度 f のy成分 fy の値に関する等高線図です。

(7) Contour(fluxd_m)
熱流束密度 f の絶対値に関する等高線図です。円形の穴の上端と下端部で値が最大となっています。

(8) Elevation(temp) from (0, 0) to (Lx, 0)
境界下端に沿って温度がどう変化しているかをグラフ化したものです。

次へ

page_top_icon