空間分布をもつ時間発展問題の例としての拡散(熱伝導)問題

【注: この行のように【】でかこまれた段落は脚注のようなものです。】

空間分布をもつ量の時間発展を考える必要がある問題がある。 そのなかで比較的簡単な例として、ここでは、 「量 q の流れが、量 q の空間微分(傾き、勾配)に比例する」 という問題を考えてみよう。 量を物質の濃度としたとき、これは拡散問題である。 量を温度としたとき、これは熱伝導問題である。

【温度の流れというのは無理があるが、エネルギーが「温度×熱容量」 (熱容量は定数)と書けるとし、エネルギーの流れを考えるのである。 なお、温度はふつう T と書くが、Fortran では T と t の区別が つかないので、変数名として T は避けることにし、 ここでは温度も q としておく。】

式の上では偏微分方程式というものになるが、こわがらないでほしい。

話を簡単にするために、q は、y 方向、z 方向にはどこまで行っても 変わらないとしよう。そうすると、x 方向の流れと傾きだけを考えればよい。 これが、空間1次元の拡散問題である。 建物の壁を通しての熱伝導を考えるとよいだろう。

ある面の単位断面積あたり、単位時間あたりに q が通過する量を、 q のフラックス密度(flux density)という。 ここでは、xの正の方向に向かうフラックス密度を F と書いておく。 Fが、q の空間微分に比例するということは、

F = - κ ∂q/∂x ... (1)

のように書くことができる。ここで、「∂」は偏微分記号である。 今の例では、qは空間座標xと時間tの関数q(x,t)だと考えている。 この場合、 tの値を固定してxだけを変化させた場合のxによる微分を「xによる偏微分」といい、 ∂q/∂xのように書く。 xの値を固定してtだけを変化させた場合のtによる微分を「tによる偏微分」といい、 ∂q/∂tのように書く。 (もっと一般的には、数量は3次元の空間座標と時刻の関数だと考えることが多い。 ただし、 とくに気象の場合は鉛直座標として高さ、気圧などなん通りもの量が使われるので、 偏微分を使う場合は、数量が何の関数であるのかを明示してから話を始める 必要がある。)

q の大きいほうから小さいほうへqが動くのがふつうなので、κが正になる ように、負の符号をつけた。κを拡散係数という。 SI(国際)単位を採用する場合には、 時間の単位は秒、長さの単位は m なので、 拡散係数の単位は m2 / s である。 対象によっては、これと違う時間や長さの単位を使ったほうが計算につごうが よい場合もある。その場合も、 拡散係数の単位は「(長さの単位)2/(時間の単位)」である。

【熱伝導の場合は、温度Tの勾配に対してエネルギーフラックス密度Feを考えて
Fe = - k ∂T/∂x ... (2)
の k を熱伝導係数(熱伝導率)ということが多い。 体積あたりのエネルギーが 比熱容量(質量あたりの熱容量、比熱ともいう)c×密度ρ×温度Tであると考え、 温度のフラックスをFTとすれば、
Fe = c ρ FT ... (3)
(1)式を温度について書けば
FT = - κ ∂T/∂x ... (4)
比較すると、熱伝導係数kはcρκに等しいことがわかる。】

【また、流体の運動に関する動粘性係数は、運動量の拡散係数に相当する。】

ところで、Fだけがqの時間変化の原因であるとすると、qの時間変化は どのように書くことができるだろうか。 幅 Δx、Δy、Δz の箱をとって考えてみる。 x方向の両端の座標を、x, x+Δxとする。 箱の中の q の総量は、およそ q Δx Δy Δz である。 箱の壁を通しての、時間Δtの間の q の出入りは、およそ { - F(x+Δx) + F(x) } Δy Δz Δt である。

したがって、

{q(t+Δt) - q(t)} Δx Δy Δz = { - F(x+Δx) + F(x) } Δy Δz Δt

{q(t+Δt) - q(t)}/Δt = -{F(x+Δx) - F(x)}/Δx

Δt → 0, Δx → 0 とすれば、

∂q/∂t = -∂F/∂x ... (5)

となる。

(1)と(5)を組み合わせて F を消すと

∂q/∂t = κ ∂2q/∂x2 ... (6)

となる。ただし、κが定数であると仮定した。κがxによって、あるいは qの値によって変わる場合は、

∂q/∂t = ∂/∂x { κ ∂q/∂x } ... (7)

と書かなければならない。ここでは(6)でよいとしよう。

この問題の答えを具体的に求めるためには、 t = 0 での各空間点でのqの値(初期条件)のほかに、 t > 0 での領域の両端の q に関するなんらかの条件(境界条件)が必要である。 その一例としては、常に、

という定数で与えるという固定境界条件がある。 境界条件の与えかたは、対象とする問題によって使い分けなければならない。

微分方程式(6)を差分近似してシミュレーションを行なう。 「前方差分」その他の差分法の用語に関する教材は [別ページ]参照。

拡散問題は減衰問題に似ているので、時間に関して中央差分を使うのは 危険である。ここでは前方差分を使う場合を例に示す。 つまり、dq/dtを

{q(x, t+Δt) - q(x, t)}/Δt

で近似する。

空間に関しては、x方向に順次積分するわけではないし、 xの正・負の両方向を公平に扱うのが望ましいので、中央差分をとる。 1階差分をすなおにとると(ここで時刻はtで固定するので関数のtへの依存性は省略して書く)、

{q(x+Δx) - q(x-Δx)}/2Δx

だが、ここではわざと、幅Δx / 2をもとにした中央差分

{q(x+Δx/2) - q(x-Δx/2)}/Δx = QX(x) [と臨時に名づける]

を使うことにする。もう1回同様な差分をとると、

{QX(x+Δx/2) - QX(x-Δx/2)}Δx = {q(x+Δx) - 2 q(x) + q(x-Δx)}/(Δx)2

となって、実際にΔxの整数倍でない位置の値を扱う必要はなくなる。

以上の考え方をプログラムにしたものを用意した。
Unix用(改行コードLFのみ)Windows用(改行コードCRLF)
時間発展シミュレーションFortran kakusan.f kakusan.f
C kakusan.c kakusan.c
Pascal kakusan.p kakusan.p

プログラムに書きこんである定数の値は
       kappa   =   1.0e-6 拡散係数、単位例 m2 / s
       deltat  =  25.0    時間間隔         s
       deltax  =   0.01   空間間隔         m
       nt      = 100      時間ステップ数
である。 ただし、1.0e-6 (Fortranの倍精度の場合は1.0d-6と表記する)は 1.0×10-6を表わす。
境界条件は
       q0      =  20.0
       q1      =  30.0
とした。

また、初期条件は、領域全体でq1と同じ値とした。 (ただし、x = 0の壁では境界条件を優先することにした。) たとえば、へやの外の気温が 30℃で、中の気温をある時刻から 20℃に たもつことができたと仮定したとき、 厚さ10 cmの中の温度分布を求める問題と考えてもよいだろう。 また、実際に求めたいのは温度分布よりも 壁を通じての熱伝導による熱フラックス密度(時間あたり、面積あたりに通過 するエネルギー量)かもしれないが、このプログラムでは直接にそれを求めて はいない。

計算結果は図でみたほうがよいのだが、 ここでは、数字の形で書き出して、それを目で見ることにする。 ファイルに redirect して、lessなどで見ることを勧める。 そのため、ある時刻の各格子点での q の値を横1行にならべて書くようにした。
Fortranの場合を例に、プログラムの構成をもう少し説明する。 このプログラムでは、「配列変数」を使っている。 Fortranには、配列宣言専用のDIMENSION文もあるが、ここでは型宣言文で 兼ねる方法を示す。
      REAL*8 q(0:nx)
という宣言は、倍精度実数型の変数 q(0), q(1), ... q(nx-1), q(nx) を まとめて確保するものと考えればよい。かっこの中の式を、添え字という ことがある。配列宣言では、かっこの中で添え字のとりうる値の下限、 上限を宣言する。下限は省略することができ(Fortranでは省略するのが もともとの形)、省略すると1とみなされる。(他の言語では違うので注意。 Pascalでは下限は省略できない。Cでは下限は0に固定されている。) なお、ここで使うnxは、上で
      PARAMETER (nx=10)
として値を決めた定数である。nxを変数にすると、Fortranコンパイラは、 配列のための記憶領域をどれだけとったらよいのかわからないので、 コンパイルすることができない。直接
      REAL*8 q(0:10)
と書いても、機械語に翻訳される内容は同じなのだが、差分の区間の数を 変えたくなったとき、別の意味で出てきた10という定数をまちがえて変えて しまわないように、名前をつけておいたのである。
まず、初期条件の数値を配列qoldに与えて、その数値を書き出す。
それから、ktの値を1つずつ変えながら時間ステップを進めていく。 まず、両端以外のすべての格子点について、 qoldの値を参照しながら、新しいqの値を計算する。 両端の格子点については境界条件に基づく値を与える。 結果を書き出す。 それから、qの値をqoldにコピーしておく。 ここまでが各時間ステップの作業である。

プログラムの実行方法(FortranまたはCの場合)は次のようになる。 ただし、kakusan.txtは計算結果を書き出すファイル名の例である。 (Windows上では「./a.out」のところは 「./a」とすればよいはず。) 例文の左端の「%」は、OS (シェル)の出すプロンプト文字を代表させたもので、 たとえば「bash$」かもしれない。 ユーザーが打ちこむ文字には含まれない。

% g77 kakusan.f
% ./a.out > kakusan.txt
または
% gcc kakusan.c -lm
% ./a.out > kakusan.txt

実行できたら、 プログラムを編集し、 拡散係数kappa, 空間間隔deltax, 時間間隔deltat, 時間ステップ数ntなどを 書き変えて動かしてみてほしい。 (deltatをある程度大きくすると結果が拡散問題に対応しなくなると思う。) 境界条件を変えてみるのもよいだろう。 結果を比較したい場合は、出力ファイル名をうまくつけて区別してほしい。

結果の作図

数値モデルの計算結果は、大量の数値になることが多い。 結果が正しく数値と読める形で出力されているかどうかという意味では 数字を目で見て確かめる必要もある。 (数値と読めない形になってしまう例だが、 Fortranでは、計算の結果の桁数が多すぎて決められた書式におさまらない ときは、「*****」のような形で表示される。)

多くの数値を目で見て確認するために、図にしてみることがよく行なわれる。 今の問題は、空間1次元だが、時間を合わせると2次元であり、 計算結果は2次元時空間に分布するデータとして見ることもできる。 [別ページ]の図は、 上記とほぼ同じだが、数量の桁が大きく違う(ほか細かい違いもある)例文の計算結果の 等値線表示(横軸が時間、縦軸が空間)である。 作図したソフトウェアは、 1993-1996年度に増田が都立大学で作った グラフィックス(GKS)の教材である。

gnuplotによる作図

プログラム例を次のように用意した。

Unix用(改行コードLFのみ)Windows用(改行コードCRLF)
作図 (gnuplot)時間変化グラフ kakusan_time.gplt kakusan_time.gplt
時間・空間(x)等値線 kakusan_cont.gplt kakusan_cont.gplt
空間(x)分布グラフ kakusan_x.gplt kakusan_x.gplt
出力データの形式変更 (Awk)等値線用 kakusan_txv.awk kakusan_txv.awk
空間(x)分布グラフ用 kakusan_txv_sepa.awk kakusan_txv_sepa.awk
計算プログラムkakusan.{c,f,p}の出力形式は、gnuplotで いろいろな深さの温度の時間発展のグラフを作るのに便利なようにしてある。 gnuplotの中で
gnuplot> plot "kakusan.txt" using 1:2, "kakusan.txt" using 1:5
のようにすれば、書き出されたうち、の1(=2-1)番め、4(5-1)番めの場所の q (温度)の時間変化が同じ図に表示されるはずである。 (「gnuplot>」はgnuplotが出すプロンプトであり、 ユーザーが打ちこむ文字には含まれない。) また、例として用意したgnuplot用のスクリプト(コマンド群)ファイル kakusan_time.gpltを次のような形で呼び出せば、 多数の位置の温度変化をいっしょにグラフにできるようにしてある。
gnuplot> load "kakusan_time.gplt"
このスクリプトファイル例では、読むファイルの名前をkakusan.txtに固定してある。 複数のシミュレーションを比較する場合は、 スクリプトファイルを編集して使ってほしい。

時間(t)と空間(x)の2次元の等値線表示をgnuplotを使ってするためには、 数値データをならべかえておく必要がある。 1行に3つ、この場合で言えば、時間、x、温度の値があり、 同じ時刻の値を示す行がつながったあと、次の時刻の値の前に空行がある。 さきのkakusan.txtを読んでこの形式に変換するプログラム kakusan_txv.awkを用意した。 次のようにして実行する。
% gawk -f kakusan_txv.awk kakusan.txt > kakusan_txv.txt
そして、gnuplotの中で次のように実行すると等値線図が作られる。 スクリプト中で読むファイルの名前は便宜上kakusan_txv.txtで固定しているが 必要に応じて変更してほしい。 "unset surface" などの指定は、 gnuplotのバージョンによっては "set nosurface"などの形に変える 必要があるかもしれない。
gnuplot> load "kakusan_cont.gplt"

位置(x)に対する温度分布のグラフをgnuplotを使って作るためにも、 数値データをならべかえておく必要がある。 この場合も、1行に3つ、時間、x、温度の値がある形だが、 時刻ごとにファイルを作る。 上のkakusan_txv.txtから分ける方法もあるが、 ここでは、kakusan.txtから組みかえるプログラム kakusan_txv_sepa.awkを用意した。 ただし、作図に必要な時間ステップについてファイルを作ったほうがよいので、 時間間隔をintervalで指定するようにした。 書き出されるのは、時刻の値がintervalで割り切れる場合に限る。 例として、tの間隔500ごとに書き出すことにしてみる。
% gawk -v interval=500 -f kakusan_txv_sepa.awk kakusan.txt
その結果、kakusan_500.0.txtなどの複数のファイルができる。 それをgnuplotで次のようにして作図することができる。
gnuplot> load "kakusan_x.gplt"
"set style data" の指定は、 gnuplotのバージョンによっては "set data style"の形に変える 必要があるかもしれない。

今回の教材用のAwkプログラムは、 ここで作成したkakusan.txtの 時間ステップ、数値の桁数、空間座標値の記録のしかた などに対応するように作ったものであり、汎用的ではない。 他のデータに適用する際にはあちこち変更が必要になると思う。 実際、地中熱伝導の例に適用する場合は 基本的に同じプログラムを細かく調整した。

(1999-07-16), 2008-06-15, 2008-06-18, 2008-06-24
増田 耕一 (MASUDA Kooiti)