第1-4回　非定常拡散方程式（ノイマン条件）[python]

今回は境界条件が少し違う場合を行います。今回の条件は下記になります。

$$\left(\frac{\partial u}{\partial x}\right)_{wall} = \alpha$$

つまり境界における勾配を\(\alpha\)にするというわけです。
この勾配\(\alpha\) を用いると、次のような問題を解くことに相当します。

“長さ1.0 mの物体があり、その中心部で100 Kとなる分布を持っています。両端で断熱壁とします。その時の温度分布は時間とともにどのようになっていくでしょうか。”

断熱壁とは温度の出入りがない壁のことです。
つまり断熱壁では温度勾配が常に生まれません。そのため、境界条件は次の通りになります。

$$\left(\frac{\partial u}{\partial x}\right)_{wall} = 0$$

このままでは計算できないので、近似を行います。壁面における変化率が0になればよいので次のように書けます。

$$(\Delta u)_{wall}=0$$

$$ u_1-u_0=0 ⇔ u_0=u_1$$

同様に反対の壁面でも計算できます。

ただし、この場合の近似は壁面の位置が位置0と位置1の中間にあると近似していることを頭に入れておいてください。

条件を次のように設定し、計算を行います。

時間 : \(t=0.1\)
時間方向にとる点の数 ： \(501\)
x方向の長さ : \(L=1.0\)
x方向にとる点の数 ： \(xpoint=51\)
拡散係数 : \(D=1.0\)
初期条件：\(u(x,0)=400x(1-x)\)
境界条件：\(\left(\frac{\partial u}{\partial x}\right)_{wall}=0\)

コード上において、前回(第1-3回)と違うところは下記の2点のみです。

$$ u_0^n=u_1^n$$ $$u_{ xpoint -1}^n=u_{ xpoint -2}^n $$

したがって最後に上記のuを上書きすればよいわけです。（65,66行目）

結果になります。

両端が断熱壁ですので、中心の温度が拡散し、壁面付近の温度が上昇して一定温度に収束していく様子が確認できました。

by hide

import numpy
from scipy.sparse import dia_matrix, identity
import matplotlib.pyplot as pyplot
import matplotlib.animation as animation

# -----------------------------
# -- initial value           --
# -----------------------------
time = 0.1                     # time : 時間
time_point = 501               # time_point : 時間方向にとる点の数
x_L = 1.0                      # x_L : x方向の長さ
x_point = 51                   # x_point : ｘ方向にとる点の数
D = 1.0                        # D : 拡散係数

dt = time / (time_point-1)
dx = x_L / (x_point-1)
c = D*dt / (dx**2)

# -----------------------------
# -- create matrix           --
# -----------------------------
"""
I.C.   u(x,0) = 400x(1-x)
B.C.   u[0] = 0
       u[n-1] = 0
"""
""" x """
x = []
for i in range(x_point):
    x.append(dx*i)

""" u """
u = []
for i in range(x_point):
    u.append(400*x[i]*(1-x[i]))

""" matrix A """
data = numpy.array(
    [numpy.ones(x_point), -2.0*numpy.ones(x_point), numpy.ones(x_point)])

data[2][1] = 0
data[0][x_point-2] = 0
data[1][0] = 0
data[1][x_point-1] = 0

offsets = numpy.array([-1, 0, 1])
B = dia_matrix((data, offsets), shape=(x_point, x_point))
E = identity(x_point)
A = E+c*B

# -----------------------------
# -- preparation             --
# -----------------------------
fig = pyplot.figure()
img = pyplot.plot(x, u)
ims = []
ims.append(img)

# -----------------------------
# -- main                    --
# -----------------------------
for i in range(time_point):
    u = A.dot(u)

    u[0] = u[1]
    u[x_point-1] = u[x_point-2]

    img = pyplot.plot(x, u)
    ims.append(img)

pyplot.xlabel("length L [m]")
pyplot.ylabel("temperature [K]")

ani = animation.ArtistAnimation(fig, ims, interval=10)

pyplot.show()