第3-1回　定常移流拡散方程式（離散化）[python]

定常移流拡散方程式とはなにか？

はじめに

$$c\frac{\partial u(x)}{\partial x}=D\frac{\partial^2 u }{\partial x^2}$$

今回解く方程式は上に示します、定常移流拡散方程式になります。
（c,Dは定数 ）

左辺が移流項と呼ばれるxの一階微分項、右辺が拡散項と呼ばれるxの二階微分項になります。

解析解を求める

まずはこの方程式の解析解を求めてみたいと思います。

境界条件は\(u(0)=u_0, u(L)=u_L\)とします。\(A, B\)を定数として次のように変形します。\((c\neq0\))

$$\frac{c}{D}\frac{\partial u(x)}{\partial x}=\frac{\partial^2 u }{\partial x^2} ⇔ \frac{c}{D}u =\frac{\partial u }{\partial x}+A $$

両辺に\(exp(-\frac{c}{D}x)\)を掛けます。

$$⇔ (\frac{c}{D} \cdot exp(-\frac{c}{D}x))u – exp(-\frac{c}{D}x)
\cdot \frac{\partial u }{\partial x}=Aexp(-\frac{c}{D}x)$$

$$⇔ \frac{\partial}{\partial x} (-exp(-\frac{c}{D}x) \cdot u)=A (exp(-\frac{c}{D}x) $$

$$⇔ -exp(-\frac{c}{D}x) \cdot u= \frac{D}{C} A (exp(-\frac{c}{D}x)+B $$

$$⇔ u= -\frac{D}{c} A – Bexp(\frac{c}{D}x) $$

境界条件より

$$\begin{eqnarray} \begin{cases} u_0= -\frac{D}{C} A – B & \\
u_L= -\frac{D}{c} A – Bexp(\frac{c}{D}L) & \end{cases} \end{eqnarray}$$

連立方程式を解くと、下記のように\(A,B,u\)が求まります。

$$A=-\frac{c}{D}\frac{u_{0}exp(\frac{c}{D}L)-u_{L}}{exp(\frac{c}{D}L)}-1, B= \frac{u_{0}-u_{L}}{exp(\frac{c}{D}L)-1} $$

$$u= \frac{ exp(\frac{c}{D}L)- exp(\frac{c}{D}x) }{ exp(\frac{c}{D}L) -1}u_{0}+ \frac{ exp(\frac{c}{D}x)-1 }{ exp(\frac{c}{D}L) -1}u_{L} $$

また、\(c=0\)のときは次のように求めることができます。

$$0= \frac{\partial^2 u }{\partial x^2} ⇔ u=Ax+B$$

境界条件より、 下記のように\(A,B,u\)が求まります

$$A=\frac{u_L-u_0}{L}, B=u_0$$

$$u= \frac{u_L-u_0}{L}x+u_0$$

以上のことからペクレ数\(Pe=\frac{c}{D}L\)によって、解が変化することがわかります。

数値解析の為の準備

解析解は求まったので、数値解析に向け、左辺に対して後退差分,右辺に対して二階微分の中心差分を取り式変形していきます。

$$c\frac{u_i-u_{i-1}}{\Delta x}=D\frac{u_{i+1}-2u_i+u_{i-1}}{\Delta x^2}$$

ここで,\(r=\frac{c\Delta x}{D}\)と置き、式変形を行います。

$$r\left(u_i-u_{i-1}\right)=u_{i+1}-2u_i+u_{i-1}$$

$$(1+r)u_{i-1}-(2+r)u_i+u_{i+1}=0$$

この式を行列で書くと次のようになります。

$$\left[\begin{array}{cccccc}1&0&&&& \\ 1+r&-2-r&1+r&&& \\ &\ddots& \ddots & \ddots && \\&& \ddots & \ddots & \ddots &\\ &&&1+r&-2-r&1+r \\ &&&&0&1 \end{array}\right] \bf u= \left[\begin{array}{c} u_{0,bd} \\ 0 \\ \vdots \\ \vdots \\0 \\ u_{n-1,bd} \end{array}\right] $$

ただし，最初と最後の行はディリクレ条件で境界の値を与えているため値が変わっています。

この行列を解くことで解を求めることができます。

以前紹介したガウスの消去法でもこの行列は解くことができますが，今回は三重対角行列です．しかし\(0\)の部分に対して計算するのは無駄なので，三重対角行列に特化したThomasのアルゴリズムを次回に紹介したいと思います。

by hide