土壌の温度の季節変化を考える。 原因は地表面での熱収支であり、 地中では熱伝導でその影響が伝わっていっているはずである。
これを1次元拡散問題[別記事参照]として考え、
差分計算してみる。今回は、次のような境界条件を設定する。
(鉛直座標にはzを使うことが多いのだが、
ここでは、1次元拡散問題の例文とそろえるため、xを使っている。)
x = 0)では、温度を、1年周期の正弦波(サインカーブ)として与える。
地表面温度の平均値q0、振幅qampと、
位相を示すものとして1年のうちで地表面温度が極大となる季節
tphaseを与えている。(tphase = 0は、
計算開始の時点を地表面温度の極大と一致させることを示す。)x = nx * deltax) では、温度勾配がないものとする。
エネルギーの流れの原因が熱伝導だけなので、これは、そこでエネルギーの
流れがないとするのと同じことである。
具体的には、1つ余分に格子点をとり、
毎ステップ、q(nx+1) = q(nx)とする。プログラム例を次のように用意した。
| Unix用(改行コードLFのみ) | Windows用(改行コードCRLF) | ||
|---|---|---|---|
| 時間発展シミュレーション | Fortran | tion.f |
tion.f |
| C | tion.c |
tion.c | |
| 作図 (gnuplot) | 時間変化グラフ | tion_time.gplt |
tion_time.gplt |
| 時間・深さ等値線 | tion_cont.gplt |
tion_cont.gplt | |
| 鉛直分布グラフ | tion_depth.gplt |
tion_depth.gplt | |
| 出力データの形式変更 (Awk) | 等値線用 | tion_txv.awk |
tion_txv.awk |
| 鉛直分布グラフ用 | tion_txv_sepa.awk |
tion_txv_sepa.awk |
物理量の単位は、計算中では、SI (m, kg, s, K)でなるべく統一した。 ただし、温度は、 今回の問題の範囲では絶対温度のケルビン(K)とセルシウス(摂氏, C)度の どちらを使ってもプログラムを変更する必要がないので、 例では(生活上よく使う)セルシウス度を想定した。 黒体放射のように絶対温度を含む式を組み合わせる場合は注意が必要である。 また、出力では、時間の単位を「年」にした。このため換算をしている。
温度拡散係数kappa (近藤 2000では「a」(122ページ))は、
「熱伝導率/(比熱容量×密度)」である。
その値は、土壌の種類とそれが水をどれだけ含んでいるかによって変わるが、
ここでは近藤(2000) 表4.2の「温度拡散係数の目安」のうち
「乾燥砂地・粘土」の値を採用して、
2.3×10-7 m2 / sとした。
計算する深さの範囲は 4 m、間隔は10 cmとした。
また、時間ステップは (1/5000)年とした。
(時間ステップを(1/1000)年とすると、この時間空間間隔とkappaの
値では差分計算がうまくいかない。これは差分方式に時間に関する前方差分を
使っているためで、差分方式を変えればこの程度でも安定して計算できるだろう。
前方差分のままで、kappaを大きくした場合、たとえば
近藤(2000)の「湿り砂地・粘土」の値
6.7×10-7 m2 / sとした場合は、
差分の時間間隔をもっと短くする必要があるかもしれない。)
全部の時間ステップの全部の格子点の値を書き出すとファイルが大きくなり、
実行時間もかかってしまうので、
深さについてはkxouti (= 2)個にひとつ、
時間についてはktouti (= 25)個にひとつを
書き出すようにした。
(等間隔でない書き出しも、プログラムがもう少し複雑になるが、
不可能ではない。)
【C言語版のプログラムは、Fortran版のプログラムを
f2cというプログラムで変換したものから
手作業で書きかえて作った。
f2cの出力は、f2cのプログラムライブラリ
(部品集)をリンクする必要があるが、必要となるライブラリの大部分は、
Fortranの入出力(WRITE文など)の仕事をCで表現したものである。
これとCの標準入出力関数(printfなど)とはぴったりとは対応しないが、
その場の意図を考えれば書きかえることができる。
書きかえた結果のtion.cはf2cのライブラリを
必要としない。Cコンパイラに付属する数学関数ライブラリは
必要とする(gccコマンドで-lmを指定すること)。】
実行方法は次のようになる。
ただし、tion.txtは計算結果を書き出すファイル名の例である。
(Windows上では「./a.out」のところは 「./a」とすればよいはず。)
例文の左端の「%」は、OS (シェル)の出すプロンプト文字を代表させたもので、
たとえば「bash$」かもしれない。
ユーザーが打ちこむ文字には含まれない。
% g77 tion.f % ./a.out > tion.txt
% gcc tion.c -lm % ./a.out > tion.txt
この部分のプログラム例は、
一般的熱伝導の例題のために
作ったkakusanシリーズの例文と基本的に同じだが、
時間・空間の数値の範囲、時間ステップ、数値の桁数などについて
細かい変更をしたものである。
gnuplot> plot "tion.txt" using 1:2, "tion.txt" using 1:3, "tion.txt" using 1:21gnuplot>」はgnuplotが出すプロンプトであり、
ユーザーが打ちこむ文字には含まれない。)
また、例として用意したgnuplot用のスクリプト(コマンド群)ファイル
tion_time.gpltを次のような形で呼び出せば、
もう少し整ったグラフができるようにしてある。gnuplot> load "tion_time.gplt"tion.txtに固定してある。
必要に応じて、スクリプトファイルを編集し
読むファイル名、深さのレベル番号、グラフの目盛りのつけかたなどを変更して
実行してみてほしい。
いろいろな深さの温度の時間変化をグラフにし、
勝手な初期条件からどのようにして1年周期の同じ値のくりかえしに近づいていくか、
また、年変化の振幅と位相が深さによってどう変わるか、見てほしい。tion.txtを読んでこの形式に変換するプログラム
tion_txv.awkを用意した。
次のようにして実行する。 % gawk -f tion_txv.awk tion.txt > tion_txv.txttion_txv.txtで固定しているが
必要に応じて変更してほしい。
"unset surface" などの指定は、
gnuplotのバージョンによっては "set nosurface"などの形に変える
必要があるかもしれない。gnuplot> load "tion_cont.gplt"tion_txv.txtから分ける方法もあるが、
ここでは、tion.txtから組みかえるプログラム
tion_txv_sepa.awkを用意した。
ただし、作図に必要な時間ステップについてファイルを作ったほうがよいので、
開始時刻、終了時刻、時間間隔をそれぞれfrom, to, intervalで指定するようにした。
ただし、書き出されるのは、時刻の値がintervalで割り切れる場合に限る。
例として、(定常的1年周期変動に近づいたところの季節変化を見たいので)
8年めの1年間について1/8 (=0.125)年ごとに書き出すことにしてみる。
(実行例の行が長いので折り返して表示されるかもしれないが1行である。)% gawk -v from=8 -v to=8.875 -v interval=0.125
-f tion_txv_sepa.awk tion.txttion_8.000.txtなどの複数のファイルができる。
それをgnuplotで次のようにして作図することができる。gnuplot> load "tion_depth.gplt"現実の地中温度の問題に近づけるには、 前に0次元モデルの例 でやったように、地表面でのエネルギーの流れの年変化を与えたほうがよいが、 それは機会があれば別の記事でとりあげたい。
温度拡散係数が一定の場合の熱伝導問題は、 差分計算によらなくても、微分方程式を直接(数式変形で)解くこともできる。 フーリエ変換で有名なフーリエ(FOURIER)が扱った問題である。 興味のある人は、深さに伴って振幅と位相がそれぞれどう変わるか、 理論と数値計算の結果が合っているか確かめてみてほしい。
また、0次元モデルとの対比のうえでは、 鉛直方向に集計した、土壌層のもっているエネルギー量の合計の季節変化を みるとよさそうである。 プログラムに少しだけ手を加え、土壌の密度と比熱容量の値を与えれば これを計算して書き出すようにすることもできる。
日本の気象庁の測候所で、昔はいくつかの深さの地温の観測をしていた。 (わたしが小学生のとき学校で見学に行って見せてもらった。) しかし今は続けられていない。農業試験場や大学などでは観測しているところがある。 中国(一般気象観測所)、タイ(農業気象観測所)、インドネシア(気候観測所)などでは 今も毎日の地温が記録されている。 多くは紙の形で記録されているが、 一部ディジタル化されたものもある。 自由に公開されたものはあまりないが、 再配布しないという条件で教材用に提供できるものはある。 興味のある人は、観測データと数値計算の結果の特徴を比較して見てほしい。
[「ローカルなエネルギー収支の0次元モデル」のページの文献]参照。