1次元拡散方程式ソルバー 戻る
数値解析ツール

1次元拡散方程式ソルバー

$\partial u/\partial t = D \cdot \partial^2 u/\partial x^2$ をCrank-Nicolson法で安定に求解。初期条件・境界条件を自由に設定してアニメーションで確認。

パラメータ設定
拡散係数 D (×10⁻⁵ m²/s)
×10⁻⁵
領域長 L (m)
m
格子数 N
初期条件
境界条件
計算結果
0.00s
経過時間 t
0.000
フーリエ数 Fo
1.000
最大値 u_max
1.000
エネルギー E(t)
1.00
領域長 L
差分
スナップショット
エネルギー保存推移
理論・主要公式
$$\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)$$ 無条件安定・2次精度。Thomas算法(三重対角行列)で各時刻を高速求解。拡散時間スケール:$t^* = L^2/D$

1次元拡散方程式ソルバーとは

🙋
拡散方程式って、教科書で見たけど、実際に何が計算できるんですか?
🎓
大まかに言うと、何かが「広がっていく」様子を計算するんだ。例えば、コップにインクを一滴垂らすと、じわーっと広がるよね。あの現象を数値的に再現できるのがこのシミュレーターだよ。上の「初期条件」で「ガウス分布」を選んで、「拡散係数 D」のスライダーを動かしてみると、広がる速さがどう変わるかすぐにわかるよ。
🙋
え、そうなんですか!「境界条件」って何ですか?「Dirichlet」とか「Neumann」って難しそう…。
🎓
実務では、領域の端っこをどう扱うかが特に重要なんだ。「Dirichlet」は端の値を固定する設定。例えば、金属棒の両端を氷で常に冷やし続ける(温度0℃に固定)イメージだ。「Neumann」は端での出入りを指定する。断熱状態(熱が出入りしない)を再現できるよ。シミュレーターで切り替えて、分布の広がり方がどう変わるか確かめてみて。
🙋
「Crank-Nicolson法」って聞きますが、他の計算方法と何が違うんですか?
🎓
CAEの実務で安定性と精度のバランスが取れた優秀な手法なんだ。概略として説明すると、前の時刻と今の時刻の計算を「半々」に混ぜてるんだ。これで、時間ステップ$\Delta t$を大きく取っても計算が発散しにくく(無条件安定)、しかも精度が高い。このツールの裏側では、この方法で作られた方程式を「Thomas算法」という超速アルゴリズムで解いているんだよ。

よくある質問

境界条件は「ディリクレ条件(固定値)」と「ノイマン条件(勾配固定)」の2種類から選択できます。左端と右端で独立に設定可能で、例えば左端を温度100℃に固定し、右端を断熱(勾配0)にするといった使い方ができます。数値入力後「条件を更新」ボタンを押してください。
Crank-Nicolson法は無条件安定ですが、拡散係数Dが大きすぎる、空間刻みdxに対して時間刻みdtが大きすぎる、または初期条件に急峻な変化があると、数値振動が発生することがあります。まずdtを小さく(例:0.01→0.001)してみてください。それでも発散する場合は、Dの値を1桁下げて試してください。
支配方程式が同じ形なので、物質の拡散(例:インクが水中に広がる)、生物の個体群分散、半導体中のキャリア拡散、地下水汚染物質の移流なし拡散などに応用できます。物理量uの単位と拡散係数Dの値を対象現象に合わせて設定すれば、同様にシミュレーション可能です。
現時点ではグラフ表示とアニメーション再生のみの対応です。ただし、画面上のグラフは右クリックで画像保存が可能です。数値データが必要な場合は、初期条件やパラメータをメモし、ご自身でExcel等で再計算するか、今後の機能追加をお待ちください。

実世界での応用

熱工学(熱伝導解析):エンジンシリンダーやブレーキディスク内部の温度分布の時間変化を予測します。材料の熱応力や耐久性評価の基礎データとして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. 拡散係数D(m²/s)と領域長L(m)を入力。例えば鋼材の熱拡散率D=1.2×10⁻⁵m²/s、L=0.1mを設定
  2. 格子分割数Ngridを指定(推奨50~200)。Courant数Fo=Ddt/dx²≤0.5を満たす時間刻みdtが自動計算される
  3. 初期条件(例:左端100°C固定、右端0°C固定)と境界条件を設定後、シミュレーション開始。アニメーション表示で温度分布の時間発展を確認

具体的な計算例

銅パイプ(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%低下する

実務での注意点