はじめに
前回までは、1次元の移流方程式について考えてきました。今回は、多次元への拡張として2次元の移流方程式について数値計算を行っていきます。今回は、Gauss Hill問題と呼ばれる、ガウス分布が一定速度で移動する様子をアニメーションで描画してみようと思います。
シリーズ構成
1. スカラー移流方程式
1.1 輸送速度が正の場合
1.2 輸送速度の符号が不明の時の線形問題
1.3 数値流束を利用した方法
1.4 輸送速度が未知量であるとき(Burgers方程式)
1.5 多次元への拡張(2次元スカラー移流方程式)
1.6 実践的な計算法(TVD方程式)
2. 移流方程式の時間積分法
2.1 時間陽解法について
2.2 時間陰解法について(Crank-Nicholson法、近似LU分解)
3. 2階偏微分方程式
3.1 1次元熱伝導方程式の数値計算法
3.2 楕円方程式の数値計算法
4. 圧縮性流れの数値計算法
4.1 オイラー方程式の数値計算法
4.2 MacCormack法による数値計算法
4.3 TVD法による数値計算法
支配方程式
2次元スカラー移流方程式は、
\frac{\partial{q}}{\partial{t}}+c\frac{\partial{q}}{\partial{x}}+d\frac{\partial{q}}{\partial{y}}=0
ここで、$c$および$d$が正の実数であるときを考えます。一次精度風上法を利用すると、
q^{n+1}_{j,k}=q^{n}_{j,k}-c\Delta{t}\frac{q^{n}_{j,k}-q^{n}_{j-1,k}}{\Delta{x}}-d\Delta{t}\frac{q^{n}_{j,k}-q^{n}_{j,k-1}}{\Delta{y}}
と離散化できます、今回は$c=d=1.0$、$\Delta{x}=\Delta{y}=1.0[m]$として計算を行いましょう。
では、実際にプログラムを実行してみましょう。
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
import matplotlib.animation as animation
c = 1.0
d = 1.0
dt = 0.05
dx = 0.1
dy = 0.1
jmax = 41
kmax = 41
nmax = 50
ims = []
def init(dx, dy, jmax, kmax):
x = np.linspace(0, dx * (jmax-1), jmax)
y = np.linspace(0, dy * (kmax-1), kmax)
q = init_q(x, y, jmax, kmax)
return (x, y, q)
def init_q(x, y, jmax, kmax):
def gausfunc(x, y, mu, V):
detV = np.linalg.det(V)
invV = np.linalg.inv(V)
A = 1.0 / (2*np.pi*np.sqrt(detV))
dX = (np.array([x, y]) - mu[:, None, None]).transpose(1, 2, 0)
return A * np.exp(-0.5 * dX[:,:, None] @ invV[None, None] @ dX[:,:,:, None])
mu = np.array([dx * jmax / 4, dy * kmax / 4])
V = np.array([[dx * jmax / 40, 0], [0, dy * kmax / 40]])
X, Y = np.meshgrid(x, y)
q = gausfunc(X, Y, mu, V)[:, :, 0, 0]
return q
def do_computing(x, y, q, dt, dx, dy, nmax, interval = 2):
X, Y = np.meshgrid(x, y)
fig = plt.figure(figsize=(10, 7), dpi=100) # グラフのサイズ
plt.rcParams["font.size"] = 10 # グラフの文字サイズ
# 初期分布の可視化
ax1 = fig.add_subplot(1, 1, 1, projection='3d')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
for n in range(1, nmax + 1):
qold = q.copy()
for j in range(1, jmax - 1):
for k in range(1, kmax - 1):
式(2)
q[j][k] = qold[j][k] - dt / dx * (
c * (qold[j][k] - qold[j-1][k]) + d * (qold[j][k] - qold[j][k-1])
)
im = ax1.plot_wireframe(X,Y,q, color = "black", rstride = 1, cstride = 1, linewidth = 0.5)
ims.append([im])
ani = animation.ArtistAnimation(fig, ims)
plt.show()
x, y, q = init(dx, dy, jmax, kmax)
do_computing(x, y, q, dt, dx, dy, nmax)
実行すると、
ガウス分布が移動していく様子が見られました。
参考文献
藤井孝蔵、立川智章、Pythonで学ぶ流体力学の数値計算法、2020/10/20