拡散方程式は、濃度の高い場所から低い場所へ物質が広がる現象を記述する基本的な偏微分方程式です。
式だけを見ると突然現れた数学的な関係に感じられますが、その背景には質量保存とフィックの法則という明確な物理原理があります。
微小な領域へ入る物質量と外へ出る物質量の差を考え、内部に蓄積される量と結び付けることで連続の式が得られます。
さらに、物質流束が濃度勾配に比例するというフィックの第一法則を代入すれば、標準的な拡散方程式を導出できます。
拡散方程式の導出では、微小体積に対する質量収支、流束の定義、濃度勾配との関係を順番に整理することが重要です。
本記事では、一次元の直感的な導出から三次元のベクトル表現、発生項や移流を含む一般式、円筒座標や球座標への拡張まで詳しく解説します。
拡散方程式は質量保存則とフィックの法則から導かれます
それではまず、拡散方程式を導く基本的な考え方について解説していきます。
微小区間における物質量の収支を考える
一次元の棒や細い管の中を物質が拡散する状況を考えます。
位置xからx+Δxまでの微小区間を取り出し、断面積をAとします。
濃度をc(x,t)とすれば、微小区間に含まれる物質量は濃度と体積の積です。
微小区間の物質量=c(x,t)AΔx
濃度が時間とともに変化する場合、単位時間当たりの蓄積量は時間微分で表せます。
蓄積速度=AΔx∂c/∂t
一方、物質は左側の面から流入し、右側の面から流出します。
単位面積および単位時間当たりに通過する物質量を流束Jと呼びます。
左面からの流入量はAJ(x,t)、右面からの流出量はAJ(x+Δx,t)です。
内部で化学反応や発生がない場合、蓄積量は流入量から流出量を引いた値に等しくなります。
AΔx∂c/∂t=AJ(x,t)-AJ(x+Δx,t)
断面積Aを消去し、両辺をΔxで割ります。
∂c/∂t=-{J(x+Δx,t)-J(x,t)}/Δx
Δxをゼロへ近づけると、右辺の差分商は流束の空間微分になります。
∂c/∂t=-∂J/∂x
この式が一次元の連続の式です。
流束が位置によって減少する場所では、流入量が流出量を上回るため濃度が増加します。
反対に流束が位置とともに増加する場所では、流出量が多くなり濃度が低下します。
フィックの第一法則で流束を表す
連続の式だけでは、濃度cと流束Jという二つの未知量が含まれています。
そこで、拡散流束と濃度分布を結び付ける構成式が必要です。
通常の拡散では、フィックの第一法則を用います。
J=-D∂c/∂x
Dは拡散係数です。
負号は、物質が濃度の高い方向から低い方向へ流れることを示します。
濃度が右方向へ増加している場合、濃度勾配は正です。
このとき流束は負となり、物質は左方向へ移動します。
濃度勾配が急であるほど流束の絶対値は大きくなります。
| 量 | 意味 | 代表的な単位 |
|---|---|---|
| c | 単位体積当たりの物質量 | mol/m³ |
| J | 単位面積および単位時間当たりの物質量 | mol/m²s |
| D | 拡散の進みやすさ | m²/s |
| ∂c/∂x | 位置に対する濃度変化 | mol/m⁴ |
単位を確認すると、Dと濃度勾配の積は流束と同じ単位になります。
この次元の整合性は、式の符号や係数を検算する際にも役立ちます。
連続の式へフィックの法則を代入する
連続の式にフィックの第一法則を代入します。
∂c/∂t=-∂/∂x(-D∂c/∂x)
Dが位置や濃度に依存しない一定値なら、微分の外へ出せます。
∂c/∂t=D∂²c/∂x²
これが一次元の拡散方程式です。
フィックの第二法則と呼ばれることもあります。
左辺は濃度の時間変化、右辺は濃度分布の空間的な曲率を表します。
質量保存による連続の式と、濃度勾配に比例する流束を表すフィックの第一法則を組み合わせることで、フィックの第二法則である拡散方程式が得られます。
Dが一定でない場合は、単純にDを微分の外へ出せません。
∂c/∂t=∂/∂x{D(x,c)∂c/∂x}
この形では拡散係数の空間変化や濃度依存性も考慮できます。
Dがcに依存すれば非線形偏微分方程式となり、解析解を求める難易度が高まります。
連続の式を微分形と積分形から理解する方法
続いては、質量保存を表す連続の式をより詳しく確認していきます。
有限な検査体積に対する積分形を考える
質量保存則は、任意の検査体積に含まれる物質量の変化として表せます。
検査体積をV、その境界面をSとします。
内部の総物質量は濃度cを体積全体で積分した量です。
総物質量=∫V c dV
境界面を外向きに通過する流束をJとすると、単位時間当たりの正味流出量は面積分で表せます。
正味流出量=∫S J・n dS
nは境界面に対する外向き単位法線ベクトルです。
内部で発生や消費がなければ、総物質量の増加率は正味流出量のマイナスに等しくなります。
d/dt∫V c dV=-∫S J・n dS
この式が連続の式の積分形です。
検査体積の形や大きさにかかわらず成立するため、有限体積法の基礎にもなっています。
発散定理を用いて微分形へ変換する
境界面上の流束積分を体積積分へ変換するため、発散定理を使います。
∫S J・n dS=∫V ∇・J dV
検査体積が時間とともに動かず、濃度が十分に滑らかなら時間微分を積分の内側へ移せます。
∫V ∂c/∂t dV=-∫V ∇・J dV
両辺を一つの体積積分へまとめます。
∫V{∂c/∂t+∇・J}dV=0
この関係が任意の検査体積Vに対して成立するには、積分される関数自体がゼロでなければなりません。
∂c/∂t+∇・J=0
これが三次元の連続の式です。
∇・Jは流束の発散であり、微小領域から外へ出ていく流れの強さを表します。
発散が正なら周囲へ流出しているため濃度が減少し、発散が負なら流入が優勢なため濃度が増加します。
発生項や反応項を追加する
微小領域の内部で物質が生成または消費される場合、連続の式へ発生項Rを追加します。
∂c/∂t+∇・J=R
Rが正なら内部で物質が生成され、負なら反応などによって消費されます。
フィックの法則J=-D∇cを代入すると、反応拡散方程式が得られます。
∂c/∂t=∇・(D∇c)+R
Dが一定なら次の形です。
∂c/∂t=D∇²c+R
例えば一次反応によって物質が消費されるなら、R=-kcと置けます。
∂c/∂t=D∇²c-kc
この方程式は、拡散による空間移動と反応による濃度減少を同時に表します。
触媒粒子内部の反応、薬物の分解、生態系の個体群変化などに応用される形です。
発生項は質量保存則を破るものではなく、検査体積内部で別の物質やエネルギーから対象成分へ変換された量を収支へ加えたものです。
三次元空間における拡散方程式の導出
続いては、一次元の議論を三次元へ拡張する方法を確認していきます。
ベクトル形式のフィックの法則を使う
三次元空間では、濃度はx、y、zの各方向へ変化します。
濃度勾配はベクトルとして表されます。
∇c=(∂c/∂x,∂c/∂y,∂c/∂z)
流束Jも各方向の成分を持つベクトルです。
J=(Jx,Jy,Jz)
等方的な媒質では、どの方向にも同じ拡散係数Dを用いることができます。
ベクトル形式のフィックの第一法則は次のとおりです。
J=-D∇c
流束は濃度が最も急速に低下する方向を向きます。
山の斜面に例えるなら、濃度勾配は最も急な上り方向を示し、拡散流束はその反対の下り方向を示す関係です。
流束の発散を連続の式へ代入する
三次元の連続の式は∂c/∂t+∇・J=0です。
フィックの第一法則を代入します。
∂c/∂t=∇・(D∇c)
Dが空間的に一定なら、ラプラシアンを使って簡潔に表せます。
∂c/∂t=D∇²c
直交座標系におけるラプラシアンは、各方向の二階偏微分の和です。
∇²c=∂²c/∂x²+∂²c/∂y²+∂²c/∂z²
したがって三次元の拡散方程式は次の形になります。
∂c/∂t=D(∂²c/∂x²+∂²c/∂y²+∂²c/∂z²)
各方向の曲率が濃度の時間変化へ寄与します。
ある点が周囲より高濃度なら複数方向へ物質が流出し、濃度が低下します。
異方性拡散では拡散係数をテンソルで表す
結晶、繊維材料、生体組織などでは、方向によって拡散しやすさが異なることがあります。
このような媒質は異方性を持つといいます。
拡散係数を一つの数値ではなく、行列状の拡散テンソルで表します。
J=-Dテンソル∇c
対角成分は各座標方向の拡散の強さを示し、非対角成分は一方向の濃度勾配が別方向の流束へ与える影響を示します。
連続の式へ代入すると次の形です。
∂c/∂t=∇・(Dテンソル∇c)
繊維方向には速く、繊維を横切る方向には遅く拡散する材料などを表現できます。
座標軸を拡散の主方向へ合わせれば、拡散テンソルを対角化できる場合があります。
その場合、方向ごとに異なる拡散係数を持つ比較的分かりやすい式になります。
| 媒質 | 拡散特性 | 係数の表現 |
|---|---|---|
| 均一な液体 | 方向による差が小さい | スカラーD |
| 層状材料 | 層に平行な方向と垂直方向で異なる | 方向別の係数 |
| 繊維組織 | 繊維方向へ拡散しやすい | 拡散テンソル |
移流や反応を含む一般的な輸送方程式の導出
続いては、流体の移動や化学反応を含む拡張形を確認していきます。
移流による物質輸送を流束へ追加する
流体が速度vで移動している場合、物質は拡散だけでなく流れに乗って運ばれます。
この輸送を移流または対流と呼びます。
移流による流束は濃度と速度の積で表せます。
J移流=cv
拡散流束と移流流束を合わせた全流束は次のとおりです。
J全体=cv-D∇c
これを連続の式へ代入します。
∂c/∂t+∇・(cv)=∇・(D∇c)+R
非圧縮性流体で∇・v=0、Dが一定なら次の形へ整理できます。
∂c/∂t+v・∇c=D∇²c+R
左辺第二項が移流による濃度変化です。
河川中の汚染物質、大気中の微粒子、配管内の溶質などの解析に用いられます。
反応速度式を発生項へ組み込む
対象成分が化学反応によって生成または消費される場合、反応速度をRへ設定します。
一次反応なら濃度に比例する形です。
R=-kc
二次反応なら濃度の二乗や複数成分の濃度積が現れます。
R=-kc²
またはR=-kcA cB
複数成分が関与する場合は、各成分について連続の式を立て、化学量論係数に応じた反応項を与えます。
拡散速度が反応速度より十分に速ければ濃度はほぼ均一になります。
反対に反応が非常に速い場合は、表面付近で成分が消費され、内部まで届きにくくなるでしょう。
この競合を評価する無次元数としてダムケラー数やチーレ数などが使われます。
無次元化によって支配的な現象を判断する
代表長さをL、代表濃度をC、代表速度をUとして変数を無次元化します。
移流拡散方程式では、移流と拡散の比を表すペクレ数が現れます。
Pe=UL/D
Peが小さい場合は拡散が支配的です。
Peが大きい場合は移流が支配的となり、流れ方向へ物質が運ばれます。
反応と拡散の比を表す無次元数も導入できます。
代表的な反応拡散比=kL²/D
この値が大きいほど、拡散で内部へ到達する前に反応が進みやすくなります。
無次元化は単位を消すためだけの操作ではなく、移流、拡散、反応のうち、どの現象が系の挙動を支配しているかを判断する手段です。
実験条件や装置の大きさを変更する際にも、無次元数をそろえることで相似な挙動を予測できます。
座標系や条件によって変わる拡散方程式
続いては、平面以外の形状で使われる拡散方程式を確認していきます。
円筒座標では半径方向の面積変化を考慮する
円柱や管の内部で軸対称な拡散が起こる場合、濃度は半径rと時間tに依存します。
半径が大きくなるほど、物質が通過する円筒面の面積も増えるため、直交座標と同じ二階微分だけでは表せません。
軸方向および角度方向の変化がない場合の式は次のとおりです。
∂c/∂t=D・1/r・∂/∂r(r∂c/∂r)
右辺を展開すると、二階微分に加えて1/rを含む一階微分が現れます。
∂c/∂t=D(∂²c/∂r²+1/r・∂c/∂r)
中心r=0では1/rが現れるため特異に見えますが、対称性から中心の濃度勾配をゼロとします。
数値計算では中心点専用の差分式を使用することが一般的です。
球座標では球面積の増加を反映する
球状粒子や液滴の内部で中心対称な拡散を考える場合、濃度は半径rだけに依存します。
球面の面積は半径の二乗に比例して増えるため、拡散方程式は次の形になります。
∂c/∂t=D・1/r²・∂/∂r(r²∂c/∂r)
展開形は次のとおりです。
∂c/∂t=D(∂²c/∂r²+2/r・∂c/∂r)
触媒粒子内部への反応物の拡散、薬剤粒子からの放出、球状食品の加熱などに用いられます。
中心では対称条件として∂c/∂r=0を設定します。
表面では濃度指定、流束指定、外部との物質移動を表すロビン条件などを与えます。
拡散係数が変化する場合は積の微分を残す
Dが位置、温度、濃度によって変わる場合、D∇²cという単純な形は一般に成立しません。
正しい保存形は次のとおりです。
∂c/∂t=∇・(D∇c)
一次元で展開すると、Dの空間微分を含む項が現れます。
∂c/∂t=D∂²c/∂x²+∂D/∂x・∂c/∂x
Dが濃度に依存する場合は、連鎖律によってさらに非線形な項が生じます。
保存形のまま数値離散化すると、界面を横切る流束を適切に扱いやすくなります。
異なる材料が接する界面では、濃度そのものの連続性だけでなく流束の連続性が重要です。
拡散方程式を導出するときの注意点
続いては、導出で間違えやすい点を確認していきます。
流入と流出の符号を統一する
連続の式で最も起こりやすい間違いは、流束の向きと符号の混乱です。
流束Jの正方向を右向きと決めた場合、左面から正のJが入る量は流入としてプラスです。
右面から正のJが出る量は流出としてマイナスになります。
そのため、蓄積量は左面のJから右面のJを引いた形です。
最終的に∂c/∂t=-∂J/∂xとなることを確認しましょう。
さらにJ=-D∂c/∂xを代入すると負号が二つ重なり、拡散方程式の右辺は正のDになります。
濃度の定義と体積変化を確認する
濃度を単位体積当たりの物質量として扱う場合、媒質の密度や体積が一定であることを暗黙に仮定する場合があります。
圧縮性流体や膨潤する材料では、単純な濃度だけでなく密度や空隙率の変化を含める必要があります。
多孔質材料では、単位全体体積当たりの濃度か、空隙体積当たりの濃度かによって式が変わります。
空隙率をεとすれば、蓄積項が∂(εc)/∂tとなる場合もあります。
変数の定義を曖昧にすると、実験データとの比較で係数が合わなくなるため注意が必要です。
フィックの法則が適用できる範囲を理解する
フィックの法則は、多くの通常拡散で有効な近似です。
しかし、非常に短い時間や微小なスケール、複雑な媒質では、流束が濃度勾配へ瞬時に応答しない場合があります。
また、荷電粒子では濃度勾配だけでなく電場の影響も受けます。
温度勾配による熱拡散、圧力勾配、化学ポテンシャル勾配などが駆動力となるケースもあります。
その場合はネルンスト・プランク式や一般化された輸送方程式を用いる必要があるでしょう。
拡散方程式を使う前に、流束が本当にフィックの法則で表せるか、媒質が均一か、拡散係数を一定と見なせるかを確認することが重要です。
まとめ
拡散方程式は、微小領域における質量保存則と、濃度勾配に比例する拡散流束を組み合わせることで導出できます。
まず流入量、流出量、蓄積量の収支から連続の式を求め、次にフィックの第一法則を代入する流れです。
拡散係数が一定なら、一次元では∂c/∂t=D∂²c/∂x²となり、三次元では∂c/∂t=D∇²cとなります。
内部で反応や発生が起こる場合は発生項を加え、流体による移動がある場合は移流流束を追加します。
円筒や球形の領域では、面積が半径とともに変化するため、座標系に対応したラプラシアンを使用しなければなりません。
質量収支、流束、構成式という三段階を意識すれば、条件が変わっても拡散方程式を体系的に導けるでしょう。