$\partial u/\partial t = D \cdot \partial^2 u/\partial x^2$ をCrank-Nicolson法で安定に求解。初期条件・境界条件を自由に設定してアニメーションで確認。
熱工学(熱伝導解析):エンジンシリンダーやブレーキディスク内部の温度分布の時間変化を予測します。材料の熱応力や耐久性評価の基礎データとしてCAE解析で広く利用されています。
化学工学・材料工学(物質拡散):半導体製造プロセスにおけるドーパント(不純物)の拡散、金属表面への浸炭処理による炭素濃度の浸透深さの予測などに応用されます。
環境工学:地下水や土壌中における汚染物質の拡散・移流をモデル化し、汚染範囲の予測や浄化計画の立案に役立てられます。
金融工学:オプション価格の評価モデルであるブラック・ショールズ方程式は、拡散方程式の一種(対数収益率の拡散)として定式化されます。このツールで学ぶ数値解法はその基礎となります。
このツールで遊び始めるときに、つまずきやすいポイントをいくつか挙げておくよ。まず「拡散係数 D の値の現実感覚」。教科書では D=1.0 で計算されることが多いけど、実世界の値は桁が大きく異なる。例えば、空気中の水蒸気の拡散係数は約 2.5e-5 m²/s、鋼材の熱拡散率は約 1e-5 m²/s 程度だ。シミュレーターで D=1 にすると、あっという間に広がって見えるのは、この「スケール感」を無視しているから。実問題に当てはめる時は、領域の長さ L や時間スケールと合わせて、無次元数であるフーリエ数 $F_o = D \Delta t / (\Delta x)^2$ を意識しよう。Crank-Nicolson法は無条件安定だけど、精度のためには $F_o$ を大きすぎない値(例えば10以下)に抑えるのがコツだね。
次に「Neumann境界条件=0の本当の意味」。断熱や物質の出入りなしを表すのは正解だけど、これは「勾配が0」、つまり端っこで分布のカーブが平らになる条件だ。だから、初期分布が端で値を持っていると、そこから全く動かなくなる。例えば、金属棒の中心だけを加熱し、両端を断熱にしたら、熱は端にたどり着いても外に出られず、最終的には全体が一様な温度に落ち着く。Dirichlet(値固定)と混同しないようにしよう。
最後に「数値的拡散と振動」の落とし穴。このソルバーは優秀だけど、初期条件に急峻な変化(ステップ関数など)があると、計算上でわずかに「にじんで」見えたり、ごく初期の時間ステップで値が微妙に振動することがある。これは離散化に伴う避けられない現象で、メッシュを細かくする(Δxを小さくする)ことで軽減できる。実務では「計算結果が滑らかすぎるな?」と思ったら、メッシュ依存性をチェックする癖をつけよう。
銅パイプ(D=1.17×10⁻⁴m²/s)長さL=1.0m、格子数Ngrid=100の場合、Crank-Nicolson法でdx=0.01m、フーリエ数Fo=0.45のとき安定求解が得られる。初期条件u(x,0)=100sin(πx)で開始すると、t=10s(無次元時間Fo≈4.7)で最大値u_max≈8.2°Cに減少。このとき領域全体のエネルギーE(t)は初期値から約92%低下する