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

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

PDF版 (428KB)

熱伝導方程式でも例えば

(1)

のように線形で同次型（右辺が0）の場合には、重ね合せの原理が成り立ちます。すなわち T₁, T₂が方程式(1)の解であるとしたとき、T = T₁ + T₂ も(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