魔女の小さな冒険

魔女のちいさな探検

ゆっくりゆっくり進みます。

拡散方程式をガチで解く

 今回は本を見ながら計算した話です。文中に計算用語バリバリなので、計算アレルギーのある人は読み飛ばしてください(^-^)

 1次元問題で水蒸気濃度 c(x,t) の拡散方程式は

  dc(x,t)/dt = D d²c(x,t)/dx²

これを、フーリエ級数を使ってガチで解きます。
 初期状態では、x=0 が水面、0<x<L が乾燥空気とします。   c(x,0) = 0 (0<x<2L)

計算上、計算領域は x=[0,4L] まで拡張します。まず、対称性を良くするために水面を x=0、x=2L に想定して、ここでの水蒸気量を飽和水蒸気量 a とします。

  c(0,t) = a
  c(2L,t) = a

 一見、x=[0,2L] まで拡張すればきれいな周期関数になりそうですが、フーリエ級数展開してみたら水面で濃度c(0,0) のみを a にすることができませんでした(;_;) そこで、x=4Lまで拡張して、(0,2L)を0、(2L,4L)を2aの矩形波にします。このフーリエ級数展開なら本に載っています(^-^)b

 パラメーターをいろいろ調整すると、

  c(x,0) = a - (4a/π){sin(k_1 x) + (1/3)sin(k_3 x) + (1/5)sin(k_5 x) + ・・・}
  k_n = (π/2L)(2n+1)

が解です。
 時間発展まで考えると、前回 k の項は

  c_k exp(ikx-t/τ_k)
  τ_k = 1/(k² D)

としましたが、xの部分はexp関数からsin関数に代わったので、各項に時間の部分だけくっつけます。

  c(x,t) = a - (4a/π){sin(k_1 x)exp(-k_1² Dt) + (1/3)sin(k_3 x)exp(-k_3² Dt) + (1/5)sin(k_5 x)exp(-k_5² Dt) + ・・・}

 これを計算機に載せます。

 最初はエクセルに入れてみたのですが、手動ではおそろしく収束が悪いと判明。
 次に、何か良さそうな計算ツールを探したら、Windows用のパール(Strawberry Perl)を発見!エディター(emacs)も入手してPerl入門ゼミを見ながら、プチプチ打ち込みました。結果をエクセルで図表化しました。
 今度は、t=1[s] より t=0[s] のほうが拡散が進んでいることが発覚!ナニコレ(@o@)

 t=0 と t>0 では、

  exp(-k_n² Dt)

の部分が違うようです。この部分は、k_n² Dt>0 となる大きな n の項は 0 になることを意味します。
 t=0 の式にこの部分は含まれないので、n の収束は 1/n 程度です。n=1000 では、最初の項の 1/1000 くらいになるので、ここで計算をやめても誤差は 1/1000 くらいになります。ところが、誤差が 1/1000の 項は 1000 くらいあるので、級数が収束しているのかどうか評価できません。なのに、収束しているように見えます(^^;

 やってみてわかったことは、水面付近の水蒸気の濃度分布は時間によって変化することです。上田政文先生の葉の表面の水蒸気分布を調べた論文で、拡散が定常化したときの水蒸気分布が示されていたのですが、葉を置いてから何分後かという情報が極めて重要だと気付きました。なんとなく、定常状態があるような気がしていました。思い込みってコワイです(^^;

 計算の結果から1m³の箱で、見た感じほぼ全体に拡散する時間は3時間程度、完全に蒸発の止まる拡散時間は1.2日くらいです。時間を10、10²、10³、10⁴、10^5、・・・と計算しているので、この見積もりは桁しか合っていません。正確な数字は次回に続くヨ。

参考:
フーリエ級数展開 大石進一(1989)『フーリエ解析岩波書店 https://www.iwanami.co.jp/book/b260886.html
Perl入門ゼミ https://tutorial.perlzemi.com/ (参照日2019/04/17)
Strawberry Perl http://strawberryperl.com/ (入手したバージョン strawberry-perl-5.28.1.1-64bit-PDL)
Emacs http://cha.la.coocan.jp/doc/NTEmacs.html#sec3 (入手したバージョン emacs-26.1-i686-minimum-ime.zip)