$\partial u/\partial t = D \cdot \partial^2 u/\partial x^2$ をCrank-Nicolson法で安定に求解。初期条件・境界条件を自由に設定してアニメーションで確認。
支配方程式は1次元拡散方程式(または熱伝導方程式)です。物理量$u$(温度、濃度など)の時間発展を記述します。
$$\frac{\partial u}{\partial t}= D \frac{\partial^2 u}{\partial x^2}$$$u(x, t)$: 物理量(温度[K]、濃度[mol/m³]など), $t$: 時間[s], $x$: 位置[m], $D$: 拡散係数[m²/s](物質拡散係数や熱拡散率)。この式は「時間変化(左辺)は、空間的な分布のカーブの大きさ(右辺の2階微分)に比例する」ことを意味します。
本シミュレーターで用いているCrank-Nicolson法の離散化式です。時間積分に陰解法と陽解法の考え方を組み合わせています。
$$\frac{u_i^{n+1}-u_i^n}{\Delta t}=\frac{D}{2}\left(\delta^2_x u_i^{n+1}+\delta^2_x u_i^n\right)$$$u_i^n$: $n$時ステップ、$i$格子点での値, $\Delta t$: 時間刻み, $\delta^2_x$: 空間2階中心差分オペレータ。右辺の$n$と$n+1$の項を平均している点が特徴です。この式を整理すると、各時間ステップで解くべき方程式は三重対角行列の連立一次方程式となります。
熱工学(熱伝導解析):エンジンシリンダーやブレーキディスク内部の温度分布の時間変化を予測します。材料の熱応力や耐久性評価の基礎データとして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を小さくする)ことで軽減できる。実務では「計算結果が滑らかすぎるな?」と思ったら、メッシュ依存性をチェックする癖をつけよう。
この1次元拡散の考え方は、実はめちゃくちゃ多くの分野の基礎になってるんだ。まず「電池工学」。リチウムイオン電池の充放電時、電極内部でのリチウムイオンの拡散が性能を左右する。このシミュレーターで濃度 u をリチウム濃度に、D を固体内拡散係数に見立てれば、充電速度が速すぎると電極表面で濃度が飽和して劣化する現象(プラッシング)のイメージがつかめる。
次に「構造力学の世界への発展」。拡散方程式の解き方をマスターすると、実は梁のたわみの微分方程式や、一次元波動方程式の数値解法への架け橋になる。支配方程式が似ているからね。例えば、梁のたわみの基本式は $EI \frac{d^4 w}{dx^4} = q(x)$ で、これも差分法で解く。連立一次方程式を解くアルゴリズム(Thomas算法)は全く同じ土台なんだ。
さらに「画像処理」にも応用されている。画像の「ガウシアンフィルタ」や「アニスィオトロピックディフュージョン」は、拡散方程式を2次元の画像上で解くことで、ノイズ除去やエッジ保存平滑化を実現している。このツールで一次元のガウス分布が時間とともに滑らかに広がる様子を見ることは、画像がぼやけていくプロセスの本質的な理解に直結するんだ。
このツールに慣れたら、次のステップとして「2次元への拡張」を考えてみよう。1次元ではThomas算法でサクッと解けた連立方程式が、2次元ではより大きな疎行列になる。これをどう効率的に解くかが次の関門で、「反復法(例えばSOR法、共役勾配法)」の出番だ。まずは正方形領域で、両端を冷やした金属板の温度分布がどう時間変化するか、をイメージしてみるといい。
数学的な背景を深めたいなら、「偏微分方程式の分類(放物型、楕円型、双曲型)」を学ぼう。今回の拡散方程式は「放物型」の代表格。これに対して、定常状態を表すラプラス方程式(楕円型)や、振動を表す波動方程式(双曲型)との違いを理解すると、CAEで出会うあらゆる物理現象が体系化されて見えてくる。このツールのCrank-Nicolson法の式をよく見ると、時間項をどう扱うかが鍵なんだ。
最後に、実務でCAEツールを使うなら、「無次元化」の概念は必須だ。先ほど少し触れたフーリエ数もそう。方程式を無次元化すると、パラメータの数が減り、現象の本質が浮き彫りになる。例えば、拡散係数D、長さL、時間tがバラバラにある問題も、無次元時間 $t^* = D t / L^2$ を使って整理すれば、たった一つの曲線で現象を表現できるようになる。これは実験データとシミュレーションを比較する時にも超重要なスキルだから、今のうちに感覚を掴んでおこう。