\underline{円盤領域での拡散方程式解  1 拡散方程式  1}

2 次元円盤領域での拡散型方程式の解析解

竹広真一

平成 20 年 4 月 13 日

この文章では, フーリエ級数–多項式関数展開を利用した 2 次元円盤領域での流体計算のための SPMODEL ライブラリ (spml) のモジュール eq module のテストのために, 2 次元円盤領域での拡散型方程式のさまざまな境界条件に対する解を (できるだけ) 解析的に求めることを行う.

1 拡散方程式

半径 ao の円盤内の 2 次元領域での拡散問題を考える. 支配方程式は 2 次元拡散方程式

\frac{∂ψ}{∂t} = κ∇^{2}ψ  (1)

ϕ は拡散する物理量であり, 例えば温度, 組成濃度, 磁場などである. x, y はそれぞれ壁に平行および垂直方向の座標, κ は拡散率 (正数) である. ∇^{2} は 3 次元ラプラシアンであり, 2 次元極座標系 (r, ϕ) での表現は

∇^{2} = \frac{1}{r} \frac{∂}{∂r} (r\frac{∂}{∂r} ) + \frac{1}{r^{2}} \frac{∂^{2}}{∂ϕ^{2}} .  (2)

一般的な境界条件は, 外側境界円周での物理量とその r 微分の線形結合で与えられる (混合型境界条件).

α_{1} \frac{∂ψ}{∂r} + α_{0}ψ = 0, at r = a,  (3)

(4)

\underline{円盤領域での拡散方程式解  1 拡散方程式  2}

ここで α_{i} は定数係数であり, 以下の条件を満たしているものとする.

α_{1} · α_{0} ≤ 0,  (5)

ただし α_{1} と α_{0} は同時には 0 にはならないものとする.

1.1 フーリエ級数展開

方位角方向にはフーリエ級数展開して解くことができる. 以下, その 1 波数成分のみ取りだして扱う. 時間に関しても指数型の解を仮定し,

ψ(x, y, t) = ψ̃_{m}(r)e^{imϕ+σt}  (6)

これを支配方程式に代入すると,

σψ̃_{m} = κ [\frac{1}{r} \frac{d}{dr} (r\frac{dϕ̃_{m}}{dr} ) − \frac{m^{2}}{r^{2}} ϕ̃_{m}] ,  (7)

境界条件は,

α_{1} \frac{dψ̃_{m}}{dr} + α_{0}ψ̃_{m} = 0 at r = a.  (8)

1.2 最大成長率の見積もり

成長率の最大値を領域積分から見積もることができる. (2) rψ̃_{m}^{∗} をかけて r 方向に全領域積分を行い, 部分積分を繰り返すと,

σ ∫_{0}^{a} r|ψ̃_{m}|^{2}dr = ∫_{0}^{a} drrψ̃_{m}^{∗}κ [\frac{1}{r} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) − \frac{m^{2}}{r^{2}} ψ̃_{m}]
= κ ∫_{0}^{a} drψ̃_{m}^{∗} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) − κm^{2} ∫_{0}^{a} \frac{1}{r}|ψ̃_{m}|^{2}dr
= κ [rψ̃_{m}^{∗} \frac{dψ̃_{m}}{dr} ]_{0}^{a} − κ ∫_{0}^{a} r |\frac{dψ̃_{m}}{dr} |^{2} dr − κm^{2} ∫_{0}^{a} \frac{1}{r}|ψ̃_{m}|^{2}dr

ここで, 境界条件を用いて右辺第 1 項目を見積もると

[rψ̃_{m}^{∗} \frac{dψ̃_{m}}{dr} ]_{0}^{a} = \{ \begin{matrix}−\frac{α_{1}}{α_{0}} a^{2} |\frac{dψ̃_{m}}{dr} |_{r=a}^{2} & (α_{0} ≠ 0) \\ −\frac{α_{0}}{α_{1}} a^{2}|ψ̃_{m}|_{r=a}^{2} & (α_{1} ≠ 0) \end{matrix} \} ≤ 0.

eq ḋiffusion.tex  平成 20 年 4 月 13 日 (竹広真一)

\underline{円盤領域での拡散方程式解  1 拡散方程式  3}

したがって,

σ ∫_{0}^{a} r|ψ̃_{m}|^{2}dr ≤ −κ ∫_{0}^{a} r |\frac{dψ̃_{m}}{dr} |^{2} dr − κm^{2} ∫_{0}^{a} \frac{1}{r}|ψ̃_{m}|^{2}dr < −κm^{2} ∫_{0}^{a} |ψ̃_{m}|^{2}dr,
σ < \frac{−κm^{2} ∫_{0}^{a} \frac{1}{r}|ψ̃_{m}|^{2}dr}{∫_{0}^{a} r|ψ̃_{m}|^{2}dr}.

よって成長率の最大値は負であることがわかる.

1.3 一般解

支配方程式を変形すると

\frac{1}{r} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) + [\frac{−σ}{κ} − \frac{m^{2}}{r^{2}} ] ψ̃_{m} = 0.

さらに動径座標を

r̃ = \sqrt{\frac{−σ}{κ}} r, r = \sqrt{\frac{κ}{−σ}} r̃,  (9)

と変形すると,

\frac{1}{r̃} \frac{d}{dr̃} (r̃\frac{dψ̃_{m}}{dr̃} ) + [1 − \frac{m^{2}}{r̃^{2}} ] ψ̃_{m},

この微分方程式はベッセルの微分方程式であり, その原点で特異性のない一般解は m 次の第 1 種ベッセル函数 J_{m}(r̃) で表される.

ψ̃_{m} = A_{1}J_{m}(r̃) = A_{1}J_{m} (\sqrt{\frac{−σ}{κ}} r)  (10)

1.4 外側ディリクレ型境界条件の場合

外側での境界条件がディリクレ型の

ψ̃_{m} = 0 at r = a.  (11)

場合を考える. これを一般解 (10) に代入すると

J_{l} (aα) = 0.  (12)

ただし α = \sqrt{\frac{−σ}{κ}} とおいた. ベッセル函数の零点から動径波数 α と成長率 σ = −κα^{2} が定められる. 係数 A_{1} は初期条件から定まる.

eq ḋiffusion.tex  平成 20 年 4 月 13 日 (竹広真一)

\underline{円盤領域での拡散方程式解  1 拡散方程式  4}

1.5 外側ノイマン型境界条件の場合

外側での境界条件がノイマン型の

\frac{dψ̃_{m}}{dr} = 0 at y = a,  (13)

の場合を考える. これを一般解 (10) に代入して整理すると

αJ′_{m} (aα) = 0.  (14)

ただし α = \sqrt{\frac{−σ}{κ}} とおいた. 係数行列の行列式が 0 であることから, 成長率を定める分散関係に相当する式が得られる. ベッセル函数の微分の零点から動径波数 α と成長率 σ = −κα^{2} が定められる. 係数 A_{1} は初期条件から定められる^{1}.

1.6 領域外部の境界領域でラプラス方程式に従う場合

有限の熱伝導率を持つ物質で円盤内外領域が満たされている場合の熱拡散問題, あるいはポロイダル磁場の拡散の問題の場合には円盤領域の拡散方程式を次のよう

J′_{n}(x) ＝\frac{n}{x}J_{n}(x) − J_{n+1}(x)

1微分計算にはベッセル関数の漸化式を用いればよい.

eq ḋiffusion.tex  平成 20 年 4 月 13 日 (竹広真一)

\underline{円盤領域での拡散方程式解  2 非圧縮流体の粘性拡散解  5}

な混合型の境界条件の下で解くことになる^{2}.

\frac{dψ̃_{m}}{dr} |_{a} + \frac{m}{a}ψ̃_{m}|_{a} = 0,  (15)

これを一般解 (10) に代入して整理すると,

J′_{m} (aα) + \frac{m}{aα}J_{m} (aα) = 0.  (16)

ここでベッセル関数の性質

xJ′_{n}(x) ＝ xJ_{n−1}(x) − nJ_{n}(x)

を適用すると,

J_{l−1} (aα) = 0.  (17)

ベッセル函数の零点から動径波数 α と成長率 σ = −να^{2} が定められる.

2 非圧縮流体の粘性拡散解

半径 a の 2 次元円盤領域での非圧縮流体の 2 次元運動を考える. 非圧縮速度場は流線関数を用いて次のように表すことができる.

u = e_{z} × ∇ψ,  (18)
u_{r} = −\frac{1}{r} \frac{∂ψ}{∂ϕ} , u_{ϕ} = \frac{∂ψ}{∂r} .  (19)
∇^{2}ψ = 0,

\frac{1}{r} \frac{d}{dr} (r \frac{dψ̃_{m}}{dr} ) − \frac{m^{2}}{r^{2}} ψ̃_{m} = 0

ψ̃_{m}|_{r=a} = C_{1}r^{−ma}, \frac{dψ̃_{m}}{dr} |_{a} = −mC_{1}r^{−(m+1)a}

\frac{dψ̃_{m}}{dr} |_{a} + \frac{m}{a}ψ̃_{m}|_{a} = 0,

2領域外部の支配方程式は方位角方向に関してフーリエ変換するとこの方程式の独立な 2 つの解は r^{−m}, r^{m} である. ∞ で発散しない解を求めると

ψ̃_{m} = C_{1}r^{−m} (r ≥ a) .

境界 r = a で値と 1 階微分が連続であることからこれらより C_{1} を消去すると r = a での境界条件が得られる.

eq ḋiffusion.tex  平成 20 年 4 月 13 日 (竹広真一)

\underline{円盤領域での拡散方程式解  2 非圧縮流体の粘性拡散解  6}

ストークス近似の下でこれらポテンシャルの従う方程式は

\frac{∂}{∂t} ∇^{2}ψ = ν∇^{2}∇^{2}ψ.  (20)

運動学的境界条件は, 外側円周にて法線方向の速度成分が 0 となることである. 流線関数で表現すれば,

ψ = 0, at r = a  (21)

力学的境界条件が自由すべりの場合は

r\frac{∂}{∂r} (\frac{1}{r} \frac{∂ψ}{∂r} ) = 0, at r = a.  (22)

粘着条件の場合には,

\frac{∂ψ}{∂r} = 0, at r = a.

2.1 フーリエ級数展開

方位角方向にはフーリエ級数で展開して解くことができる. 以下, その 1 波数成分のみ取りだして扱う. 時間に関しても指数型の解を仮定し,

ψ(x, y, t) = ψ̃_{m}(r)e^{imϕ+σt},  (23)

これを支配方程式に代入すると,

σ [\frac{1}{r} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) − \frac{m^{2}}{r^{2}} ] ψ̃_{m} = ν [\frac{1}{r} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) − \frac{m^{2}}{r^{2}} ] [\frac{1}{r} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) − \frac{m^{2}}{r^{2}} ] ψ̃_{m}

境界条件は自由すべりの場合, 力学的境界条件が

ψ̃_{m} = r\frac{d}{dr} (\frac{1}{r} \frac{dψ̃_{m}}{dr} ) = 0,  at r = a.  (24)

粘着条件の場合には,

ψ̃_{m} = \frac{dψ̃_{m}}{dr} = 0,  at r = a.  (25)

2.2 最大成長率の見積もり

境界条件が粘着条件の場合には, 最大成長率が負となることを支配方程式を全積分することで求められる. トロイダルポテンシャルに関しては, 拡散方程式の両端ディリクレ境界条件の場合と同じであるから最大成長率が負であることが分かる.

ポロイダルポテンシャルについては, 演算子 D_{l} を

D_{l} ≡ \frac{1}{r} \frac{d}{dr} (r\frac{d}{dr} ) − \frac{m^{2}}{r^{2}}

と定義すると, 任意の 2 つの関数 f(r), g(r) に対して

∫_{0}^{a} f(r)D_{l}g(r)rdr = ∫_{0}^{a} f [\frac{1}{r} \frac{d}{dr} (r\frac{d}{dr} ) − \frac{m^{2}}{r^{2}} ] grdr
= ∫_{0}^{a} f \frac{d}{dr} (r\frac{dg}{dr} ) dr − ∫_{0}^{a} m^{2} \frac{fg}{r} dr
= [fr\frac{dg}{dr} ]_{0}^{a} − ∫_{0}^{a} \frac{df}{dr} (r\frac{dg}{dr} ) dr − ∫_{0}^{a} m^{2} \frac{fg}{r} dr
= [fr\frac{dg}{dr} ]_{0}^{a} − [r\frac{df}{dr} g]_{0}^{a} + ∫_{0}^{a} g \frac{d}{dr} (r\frac{df}{dr} ) − ∫_{0}^{a} m^{2} \frac{fg}{r} dr
= [fr\frac{dg}{dr} ]_{0}^{a} − [r\frac{df}{dr} g]_{0}^{a} + ∫_{0}^{a} g(r)D_{l}f(r)rdr

f = ψ_{m}^{∗}(r), g = D_{l}ψ_{m}(r) と選ぶと境界条件 (25) を用いて

∫_{0}^{a} ψ_{m}^{∗}(r)D_{l}D_{l}ψ_{m}(r)rdr = [ψ_{m}^{∗}r\frac{dD_{l}ψ_{m}}{dr} ]_{0}^{a} − [r\frac{dψ_{m}^{∗}}{dr} D_{l}ψ_{m}]_{0}^{a} + ∫_{0}^{a} D_{l}ψ_{m}^{∗}D_{l}ψ_{m}rdr
= ∫_{0}^{a} |D_{l}ψ_{m}|^{2}rdr

流線関数の式 (2.1) に ψ_{m}^{∗}r をかけて積分すると

σ ∫_{0}^{a} ψ̃_{m}^{∗}D_{l}ψ̃_{m}rdr = ν ∫_{0}^{a} ψ̃_{m}^{∗}D_{l}D_{l}ψ̃_{m}rdr,

左辺は拡散方程式の右辺の式変形と同様に計算できて,

σ ∫_{0}^{a} ψ̃_{m}^{∗}D_{l}ψ̃_{m}rdr = −σ ∫_{0}^{a} [r |\frac{dψ̃_{m}}{dr} |^{2} + m^{2} ∫_{0}^{a} \frac{1}{r}|ψ̃_{m}|^{2}] dr

\underline{円盤領域での拡散方程式解  2 非圧縮流体の粘性拡散解  8}

であるから,

−σ ∫_{0}^{a} [r^{2} |\frac{dψ̃_{m}}{dr} |^{2} + \frac{m^{2}}{r} |ψ̃_{m}|^{2}] dr = ν ∫_{0}^{a} |D_{l}ψ_{m}|^{2}rdr,

よって

σ = −ν \frac{∫_{0}^{a} |D_{l}ψ_{m}|^{2}rdr}{∫_{0}^{a} [r |\frac{dψ̃_{m}}{dr} |^{2} + \frac{m^{2}}{r} |ψ̃_{m}|^{2}] dr}.

2.3 一般解

拡散方程式の場合と同様に支配方程式を変形すると

\{\frac{1}{r} \frac{d}{dr} (r\frac{dψ̃_{m}}{dr} ) + [\frac{−σ}{κ} − \frac{m^{2}}{r^{2}} ]\} [\frac{1}{r} \frac{d}{dr} (r\frac{d}{dr} ) − \frac{m^{2}}{r^{2}} ] ψ̃_{m} = 0.

成長率 σ < 0 に対してのこの一般解は

ψ̃_{m} = B_{1}r^{m} + B_{2}J_{m} (\sqrt{\frac{−σ}{κ}} r)  (26)

流線関数に作用する 2 つのカッコでくくられた演算子は交換可能であることに注意されたい.

2.4 外側粘着条件の場合の流線関数

この場合の境界条件は (24) より

ψ̃_{m} = \frac{dψ̃_{m}}{dr} = 0, at r = a.

一般解 (26) をこれに代入して整理すると

( \begin{matrix}a^{l} & J_{m} (aα) \\ ma^{m−1} & αJ′_{m} (aα) \end{matrix} ) ( \begin{matrix}B_{1} \\ B_{2} \end{matrix} ) = 0,  (27)

左辺の係数行列式が 0 なる条件が固有値 α を定める分散関係となる.

| \begin{matrix}a^{m} & J_{m} (aα) \\ ma^{m−1} & αJ′_{m} (aα) \end{matrix} | = 0.  (28)
eq ḋiffusion.tex  平成 20 年 4 月 13 日 (竹広真一)

\underline{円盤領域での拡散方程式解  2 非圧縮流体の粘性拡散解  9}

この分散関係から動径波数 α と成長率 σ = −να^{2} が定められる. 係数 B_{1}, B_{2} の組合せは, 求まった α を用いて

( \begin{matrix}B_{1} \\ B_{2} \end{matrix} ) = ( \begin{matrix}J_{m}(aα) \\ −a^{m} \end{matrix} )  (29)

と与えられる.

2.5 上側自由すべり条件の場合の流線関数

この場合の境界条件は (24) より

ψ̃_{m} = r\frac{d}{dr} (\frac{1}{r} \frac{∂ψ̃_{m}}{∂r} ) = 0, at r = a

2 つめの境界条件を書き換えると

r\frac{d}{dr} (\frac{1}{r} \frac{∂ψ̃_{m}}{∂r} ) = \frac{d^{2}ψ̃_{m}}{dr^{2}} − \frac{1}{r} \frac{dψ̃_{m}}{dr} = 0. at r = a

これに一般解 (26) を代入して整理すると

( \begin{matrix}a^{m} & J_{m} (aα) \\ m(m − 2)a^{m−2} & α^{2}J′′_{m} (aα) − \frac{α}{a}J′_{m} (aα) \end{matrix} ) ( \begin{matrix}B_{1} \\ B_{2} \end{matrix} ) = 0,  (30)

左辺の係数行列式が 0 なる条件が固有値 α を定める分散関係となる. ^{3}

| \begin{matrix}a^{m} & J_{m} (aα) \\ m(m − 2)a^{m−2} & α^{2}J′′_{m} (aα) − \frac{α}{a}J′_{m} (aα) \end{matrix} | = 0.  (31)

この分散関係から動径波数 α と成長率 σ = −να^{2} が定められる. 係数 B_{1}, B_{2} の組合せは, 求まった α を用いて

( \begin{matrix}B_{1} \\ B_{2} \end{matrix} ) = ( \begin{matrix}J_{m} (aα) \\ −a^{m} \end{matrix} )  (32)

で定められる.

J′′ = −\frac{1}{r}J′ − [1 − \frac{m^{2}}{r^{2}} ] J,

3ベッセル函数の 2 階微分はベッセルの満たす微分方程式からと計算される.

eq ḋiffusion.tex  平成 20 年 4 月 13 日 (竹広真一)