工学系大学院生のブログ

2人の大学院生による雑記ブログ

第1-3回 非定常拡散方程式(ディリクレ条件)[python]

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

前回は具体的な5個の点を用いて離散化し、下記の行列に変形しました。

$$\left[\begin{array}{r}u_0^{n+1} \\ u_1 ^{n+1} \\u_2 ^{n+1} \\u_3 ^{n+1} \\u_4 ^{n+1} \end{array}\right]=\left[\begin{array}{ccccc}1-2c&c&0&0&0 \\ c&1-2c&c&0&0 \\ 0&c&1-2c&c&0 \\ 0&0&c&1-2c&c \\ 0&0&0&c&1-2c \end{array}\right] \left[\begin{array}{r}u_0^{n} \\ u_1 ^{n} \\u_2 ^{n} \\u_3 ^{n} \\u_4 ^{n} \end{array}\right] $$

ただし、\(c=D\frac{\Delta t}{\Delta x^2}\)


今回考える境界条件はディリクレ条件と呼ばれる条件です。これは次の式のように境界の値を直接与える方法です。

$$\begin{eqnarray}\begin{cases}{}u(0,t)=u_0(t)&\\u(L,t)=u_L(t)&\end{cases}\end{eqnarray}$$

5個に離散化した場合に適用すると次のように書けます。最初と最後の行が少し変わっていることに注意してください。

$$\left[\begin{array}{r}u_0^{n+1} \\ u_1 ^{n+1} \\u_2 ^{n+1} \\u_3 ^{n+1} \\u_4 ^{n+1} \end{array}\right]=\left[\begin{array}{ccccc}1&0&0&0&0 \\ c&1-2c&c&0&0 \\ 0&c&1-2c&c&0 \\ 0&0&c&1-2c&c \\ 0&0&0&0&1 \end{array}\right] \left[\begin{array}{r}u_0^{n} \\ u_1 ^{n} \\u_2 ^{n} \\u_3 ^{n} \\u_L ^{n} \end{array}\right] $$

ここで上記の式を次のように置きます。

$$\bf u^{\rm n+1}=\bf Au^{\rm n}$$

$$\bf A= \left[\begin{array}{ccccc}1&0&0&0&0 \\ c&1-2c&c&0&0 \\ 0&c&1-2c&c&0 \\ 0&0&c&1-2c&c \\ 0&0&0&0&1 \end{array}\right] $$

この\(\bf A\)は常に一定であるため、実際に時間ステップを進める方法は\(\bf A\)と、わかっている \(\bf u^{\rm n}\) の内積をとることで行えます。

さて、これらを用いて次のような問題を解きたいと思います。

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

今回は一次元ですので、条件を以下のように設定します。

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

とる点の数を決めているので、\(\Delta x\)の計算は次のようになります。

$$\Delta x =L/(x方向の空間点の数-1)$$

そのため切りのいい数字に1を足しています。また初期条件は点(0.5, 100)を頂点とし、長さL=1.0の両端が0となる二次関数としました。

結果はこの通りです。

解析結果

熱伝導により温度が徐々に下がっていることがわかります。
対象とする時間を増やせば温度が均等に0になるところまで見えると思います。ただし、今回の計算方法では時間を増やすなら時間方向にとる点も増やさないと発散してしまいます。


境界条件にノイマン条件を課したものを第1-4回で行っています。

最後に使用したコードを張っておきます。

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 : x方向にとる点の数
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)
    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()
このエントリーをはてなブックマークに追加

コメント

  • 設定を変更したので、日本標準時になったかと思います。
    ご指摘ありがとうございます。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です

トラックバックURL: https://teru-hide.com/cfd-python-unsteady-diffusion-dirichlet/trackback/
にほんブログ村 英語ブログへ
にほんブログ村