はじめに
今までは、スカラー移流方程式を扱ってきました。今回からは、二階の偏微分方程式を扱っていきます。今回はその中でも基本的な、1次元拡散方程式について扱っていきます。
シリーズ構成
1. スカラー移流方程式
1.1 輸送速度が正の場合
1.2 輸送速度の符号が不明の時の線形問題
1.3 数値流束を利用した方法
1.4 輸送速度が未知量であるとき(Burgers方程式)
1.5 多次元への拡張(2次元スカラー移流方程式)
1.6 実践的な計算法(TVD方程式)
2. 移流方程式の時間積分法
2.1 時間陽解法について
2.2 時間陰解法について(Crank-Nicholson法、近似LDU分解)
3. 2階偏微分方程式
3.1 1次元熱伝導方程式の数値計算法
3.2 楕円方程式の数値計算法
4. 圧縮性流れの数値計算法
4.1 オイラー方程式の数値計算法
4.2 MacCormack法による数値計算法
4.3 TVD法による数値計算法
陽的解法
1次元拡散方程式は以下の式で表せます。
\frac{\partial{q}}{\partial{t}}-\alpha\frac{\partial^2{q}}{\partial{x^2}}=0 \tag{1}
ここで、$\alpha$は拡散係数です。時間微分項にオイラーの陽的差分を、空間微分項に中心差分を適用すると、
\frac{q^{n+1}-q^{n}}{\Delta{t}}=\frac{q_{j+1}-2q_{j}-q_{j-1}}{\Delta{x^2}}\tag{2}
これは、次のような簡単な前進スキームになります。
q^{n+1}-q^{n}=q^{n}+\mu(q_{j+1}-2q_{j}-q_{j-1})\tag{3}
ここで、
$$\mu=\alpha\frac{\Delta{t}}{\Delta{x^2}}\tag{4}$$
である。移流方程式と同様に、von Nuemannの安定性解析を行うことができます。
$$g=1+\mu(e^{-i\theta}-2+e^{i\theta})\tag{5}$$
または、
$$g=1-2\mu(1-cos{\theta})\tag{6}$$
これから、位相誤差は生じないことがわかります。$|g|<1$であるためには、
$$-1\leqq1-2\mu(1-cos{\theta})\leqq\tag{7}1$$
$$\Leftrightarrow 0\leqq 2\mu sin^2\frac{\theta}{2} \leqq 1\tag{8}$$
任意の$\theta$について安定であるためには、
$$0 \leqq \mu \leqq \frac{1}{2} \tag{9}$$
では、実際に1次元拡散方程式を数値計算してみましょう。今回は、初期条件・境界条件を以下のように定めます。今回は、拡散現象の中でも温度拡散について扱っていこうと思います。
\left \{
\begin{array}{l}
境界条件:q(0,t)=q(1,t)=0, \\
初期条件:q(x)=sin(\pi x)
\end{array}
\right.
初期温度$sin(\pi x)$の棒の両端を$0℃$で冷やした時の熱拡散を計算することになります。
コードは以下の通りです。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def init(jmax):
x = np.linspace(0, 1, jmax)
q = np.sin(np.pi * x)
return (x, q)
def do_computing(x, q, a, dt, dx, nmax):
fig=plt.figure(figsize=(7,7), dpi=100)
plt.rcParams["font.size"] = 22
for n in range(1, nmax+1):
qold = q.copy()
for j in range(1, jmax-1):
dq = a * dt * (qold[j+1] - 2.0 * qold[j] + qold[j-1]) / (dx ** 2)
q[j] = qold[j] + dq
q[0] = 0
q[-1] = 0
im = plt.plot(x, q, color='red', lw=2)
ims.append(im)
ani = animation.ArtistAnimation(fig,ims)
plt.grid(color='black', linestyle='dashed', linewidth=0.5)
plt.xlim([0,1])
plt.xlabel('x')
plt.ylabel('q')
plt.show()
jmax = 21
nmax = 500
a = 1
dt = 0.01
dx = 0.2
x, q = init(jmax)
ims = []
do_computing(x, q, a, dt, dx, nmax)
実行すると、熱拡散の様子が見られます。アニメーションの時間が長いので一部だけ掲載したいと思います。
陰的解法(Crank-Nicholson法)
陽的解法では、時間が1次精度となってしまいます。そこで、時間微分項にも中心差分を適用することにより時間も2次精度にすることを考えましょう。そこで、第8回で扱ったCrank-Nicholson法を1次元拡散方程式に適用しましょう。すると、
\frac{q^{n+1}-q^{n}}{\Delta{t}}=\frac{\alpha}{2}\Bigl\{\frac{q^{n+1}_{j+1}-2q^{n+1}_{j}+q^{n+1}_{j-1}}{\Delta{x^2}}-\frac{q^{n}_{j+1}-2q^{n}_{j}+q^{n}_{j-1}}{\Delta{x^2}}\Bigr\}\tag{10}
式(10)を変形をすると、
(\frac{2}{\nu}+2)q^{n+1}_{j}-(q^{n+1}_{j+1}+q^{n+1}_{j-1})=(\frac{2}{\nu}-2)q^{n}_{j}+(q^{n}_{j+1}+q^{n}_{j-1})\tag{4}
もちろん、これを解くときは3項方程式の形をとります。
\left[
\begin{array}{ccc}
\frac{2}{\nu}+2 & 1 & \cdots & \cdots & 0 \\
1 & \ddots & & & \vdots \\
\vdots & &\ddots & &\vdots \\
\vdots & & &\ddots & 1\\
0 & \cdots &\cdots & 1 & \frac{2}{\nu}+2
\end{array}
\right]
\left[
\begin{array}{ccc}
\vdots \\
q^{n+1}_{j-1} \\
q^{n+1}_{j} \\
q^{n+1}_{j+1} \\
\vdots
\end{array}
\right]
=
\left[
\begin{array}{ccc}
\\
\\
\\
RHS \\
\\
\\
\\
\end{array}
\right]
\tag{11}
RHSはRight Side Handの略で式(11)の右辺の係数ベクトルを表します。
von Neumannの安定性解析もしておきましょう。
g=\frac{1-\mu(1-cos\theta)}{1+\mu(1-cos\theta)}\tag{12}
$\theta=0$のとき$g=1$となります。もし。$\mu$が0よりも大きい場合、$g$が$\theta$とともに減少して、$\pi$のときに最小値をとります。したがって、安定限界は、
\frac{1-2\mu}{1+2\mu}>-1\tag{13}
ですが、これはすべての$\mu$に対して成り立ちます。つまり、無条件で安定ということになります。ですが、計算が安定に進むということであり、その計算結果が正しいということを保証するものではありません。
では、実際に計算してみましょう。
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation
Nx = 100
Nt = 5000
delta_x = 1/(Nx-1)
delta_t = 0.05
count_T = []
count_x = []
ims = []
alpha = 0.1
TT = np.zeros([Nx, Nt])
#初期条件の設定
for i in range(Nx):
TT[:,0] = np.sin(np.pi*i)
#境界条件の設定
for i in range(Nt):
TT[0, i] = 0
TT[-1, i] = 0
#熱伝達率の行列を作成
coe = np.ones([Nx, Nt])
for i in range(Nx):
coe[i,:] = alpha
#安定性を決める数 alpha*dt/dx^2 をもつ行列の作成
q = np.zeros([Nx, Nt])
vT = np.ones([Nx, Nt])
for i in range(Nx):
vT[i,:] = delta_t * coe[i,:] / (delta_x ** 2)
fig = plt.figure(figsize=(7,7), dpi=100)
plt.rcParams["font.size"] = 18
plt.rcParams['xtick.direction'] = 'in'
plt.rcParams['ytick.direction'] = 'in'
plt.rcParams['legend.frameon'] = True
plt.rcParams['legend.edgecolor'] = 'black'
count_x.append(0)
count_T.append(0)
plt.plot(count_x, count_T, c='black', marker='o', lw=2, label='x=0.8[mm]')
#クランクニコルソン法の適用
for t in range(Nt-1):
#連立方程式の係数行列Aの作成
A = np.zeros([Nx-2,Nx-2])
for i in range(Nx-2):
A[i,i] = 2/vT[i,t] +2
if i >=1 :
A[i-1,i] = -1
if i <= Nx-4 :
A[i+1,i] = -1
#係数行列bの作成
b = np.zeros([Nx-2])
for i in range(Nx-2):
b[i] = TT[i,t]+(2/vT[i+1,t]-2) * TT[i+1,t] + TT[i+2,t] + q[i+1,t]
b[0] += TT[0,t+1]
b[Nx-3] += TT[-1,t+1]
#Aq=bをqについて解く
T = np.linalg.solve(A,b)
for i in range(Nx-2):
TT[i+1,t+1] = T[i]
count_T.append(T[i])
count_x.append(i)
im = plt.plot(count_x, count_T, c='black', marker='o', lw=2)
ims.append(im)
ani = animation.ArtistAnimation(fig,ims)
plt.gca().spines['right'].set_visible(False)
plt.gca().spines['top'].set_visible(False)
plt.xlabel('x s')
plt.ylabel('T ℃')
plt.legend()
plt.show()
まとめ
今回は、二階偏微分方程式の数値計算初回として、1次元拡散方程式を扱いました。次回は、ポアソン方程式で有名な楕円方程式を扱います。
参考文献
藤井孝蔵、立川智章、Pythonで学ぶ流体力学の数値計算法、2020/10/20