第3-3回　定常移流拡散方程式[python]

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

上記の定常移流拡散方程式を離散化し、整理し、次のような行列にしました。

$$\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] $$

ただし、 \(r=\frac{c\Delta x}{D}\) です。

今回は以下のように条件を設定します。
x方向の長さ : \(L=1.0\)
x方向に取る点の数 : \(101\)
移流項の係数 : \(c=-10, 0, 10\)
拡散項の係数 : \(D=1\)
境界条件 : \(u(0)=0, u(L)=1\)

第3-1回で説明したペクレ数\(Pe\)による差を見るために、cを変えて行います。

数値解析結果は次の通りです。

\(c=-10, Pe<0\)のとき

\(c=0, Pe=0\)のとき

\(c=10, Pe>0\)のとき

解析解通り結果がでました。

pythonのコードを最後に張っておきます。初期値cを変更すれば、上記の3つの結果を出せます。

by hide

import numpy
import matplotlib.pyplot as pyplot

# -----------------------------
# -- initial value           --
# -----------------------------
x_L = 1.0             # x_L : x方向の長さ
x_point = 101         # x_point : ｘ方向の空間分割数
c = 0.0               # c : 移流項の係数
D = 1.0               # D : 拡散項の係数
u_0 = 0.0             # u_0,u_L : 境界条件
u_L = 1.0

# -----------------------------
# -- create matrix           --
# -----------------------------
""" dx, rを計算 """
dx = x_L/(x_point-1)
r = c*dx/D

""" プロット用の位置x """
x = []
for i in range(x_point):
    x.append(i*dx)


"""  下対角成分 a_element 
     対角成分   b_element 
     上対角成分 c_element 
     右辺       b_element
"""
a_element = []
b_element = [1.0]
c_element = [0.0]
d_element = [u_0]

for i in range(x_point-2):
    a_element.append(1.0+r)
    b_element.append(-2.0-r)
    c_element.append(1.0)
    d_element.append(0.0)

a_element.append(0.0)
b_element.append(1.0)
d_element.append(u_L)

# -----------------------------
# -- function                --
# -----------------------------

def TDMA(a, b, c, d):
    """
    a,cはn-1個の要素
    b,dはn個の要素
    """
    e = [ c[0]/b[0] ]
    f = [ d[0]/b[0] ]

    for i in range(1, len(d)-1):
        e.append(c[i]/(b[i]-a[i-1]*e[i-1]))
        f.append((d[i]-a[i-1]*f[i-1])/(b[i]-a[0]*e[0]))

    u = [(d[len(d)-1]-a[len(d)-2]*f[len(d)-2]) /
         (b[len(d)-1]-a[len(d)-2]*e[len(d)-2])]

    for i in range(len(d)-2, -1, -1):
        u.append(f[i]-e[i]*u[len(d)-2-i])

    u.reverse()

    return u

# -----------------------------
# -- main                    --
# -----------------------------

u = TDMA(a_element, b_element, c_element, d_element)

pyplot.plot(x, u, 'bo')
pyplot.xlabel("length L [m]")
pyplot.ylabel("u")


pyplot.show()