Sample Scripts from GB Books
GB005:  1次元/2次元過渡問題

ここでは熱伝導方程式に時間依存の項が含まれるケース（過渡問題、transient problems）について考察します。時間依存型の問題の場合にはGB004で扱った方程式

に対し時間微分項が加わる形となります。

ここに ρ は物質の密度を、cp は比熱（specific heat capacity）を表すパラメータです。数式(2)は体積要素中におけるエネルギー収支を表現したものと言えます。なお、この数式において T は絶対温度を、t は時間を表す変数である点に注意してください。
以下においてはまず1次元の問題を扱います。すなわち数式モデル上、空間変数は x のみであるわけですが、FlexPDE上は2次元の問題として表現されている点に注意してください。

本来は1次元の熱伝導なのでy軸方向の幅は0として考えるべきかも知れませんが、FlexPDEの適性を考慮し、ここでは右図のような2次元の広がりを持った鉄棒をモデル化することにします。鉄棒の温度を T(x, t) で表すことにしたとき、その初期状態 T(x, 0) は300度Kであるとします。この鉄棒の左端に t = 0 以降熱を加え、400度Kに保ったとしたとき、鉄棒内の温度分布が時間と共にどう変化するかをFlexPDEで計算してみます。

実は鉄棒が無限に長い場合には解析解が知られていて

という数式で T が表現されます。ここに T_i は初期状態における温度、T₀ は x = 0 における温度、λ は熱伝導率を表すパラメータです。また erf は誤差関数

を意味し、FlexPDE中にも用意されています。ここではFlexPDEによる数値解がこの理論値とどの程度の一致を見せるかについても検証してみます。

まずタイトルを設定します。
  TITLE
    'Bar with a Temperature Step'    { transient01a.pde }

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

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

偏微分方程式の定義に先立ち、パラメータ類を定義します。これらは境界の定義、及び境界条件の設定に際して使用されます。解析対象のドメイン上に熱源がないため上記数式中における h は0、また鉄の電導率は82に、ρ *cp の値は 7.87e3*449 に設定されています。一方、temp_ex は厳密解を表す数式です。
  DEFINITIONS
    Lx = 1.0  Ly = 0.1
    heat = 0                  { Heat source }
    k = 82  rcp = 7.87e3*449  { Iron parameters }
    tempi = 300  temp0 = 400  { Initial and boundary temperatures }
    fluxd_x = -k*dx(temp)  fluxd_y = -k*dy(temp)
    fluxd = vector(fluxd_x, fluxd_y)  fluxd_m = magnitude(fluxd)
    temp_ex = (tempi - temp0)*erf[x/(2*sqrt(k*t/rcp))] + temp0

ここで時間依存型問題の場合に必要となる初期状態の規定を行っておきます。
  INITIAL VALUES
    temp = tempi

次に方程式を上記(2)の形で定義します。DEFINITIONSセクション中におけるfluxdの定義と合せると温度 temp に関する2階の偏微分方程式が導かれることがわかります。
  EQUATIONS
    div(fluxd) + rcp*dt(temp) = heat

今度は境界の形状と境界条件を設定します。鉄棒の左端についてはValue文で定数値を指定しますが、その他の境界についてはNatural文にて 断熱型の境界を設定しています。
  BOUNDARIES
    Region 1
      Start(0,0)
        Natural(temp) = 0  Line to (Lx, 0) to (Lx, Ly) to (0, Ly)
        Value(temp) = temp0  Line to Close

これも時間依存型問題に固有の指定項ですが、TIME セクションを用いて解析対象の時間帯を指定します。
  TIME
    from 0 to 3000

最後に出力すべき情報を指定します。ここでは鉄棒の中心軸上における温度 temp の値を、その理論値 temp_ex と合せてelevationプロットの形でプロットするわけですが、それを t = 100 から 3000 まで100秒単位に出力すべく、for ... by ... to クローズが使用されている点にご注意ください。
  PLOTS
    for t = 100 by 100 to 3000
    Elevation(temp, temp_ex) from (0, Ly/2) to (Lx, Ly/2)

  END

上記スクリプトをそのままFlexPDEに実行させるとすべてのプロットが一気に生成され、t = 3000 に対応したプロットが表示された状態で処理が終了します。最初からすべてのプロットを確認するには「File」メニュー：「View File」と操作し、該当する.pg5ファイルを選択してください。「View」メニュー：「Next」あるいは「Back」と操作することにより手動で前後に進めることができます。一方、「View」メニュー：「Movie」と操作すれば画面切替えは一定時間間隔で自動的に行われるようになります。

t = 0 における温度分布には不連続性が含まれており（x = 0 においては400度、その点以外では300度）、その扱いには難しい面があるのですが、FlexPDEはそれを無難に処理し、厳密解と良く一致したプロットを出力してきています。ここでは t = 1800 と t = 3000 のときのプロットを掲載しておくことにします。

(1) Elevation(temp, temp_ex) from (0, Ly/2) to (Lx, Ly/2) at t = 1800
t = 1800 の時点におけるelevationプロットです。temp と temp_ex のグラフが異なる色で2本描かれているわけですが、実質的に重なってしまっているため1本の曲線のように見えます。プロット下部に示されている積分値についても良い一致が見られます。

(2) Elevation(temp, temp_ex) from (0, Ly/2) to (Lx, Ly/2) at t = 3000
tが3000に近づくにつれ、鉄棒右端における温度分布には理論値との乖離が見られるようになります。これは数値演算に伴う誤差に起因するものではありません。右端における境界条件の違いによるものです。厳密解は鉄棒が無限遠まで続いているとしたときに成り立つものである点に注意してください。FlexPDEが計算したのは有限長の鉄棒中における熱伝導です。