##はじめに
数値流体解析(CFD)の基礎シリーズの第2回です。前回は、スカラー移流方程式について数値計算をしてみました。今回はその続きです。よかったら、前回の記事も見てみてください。
何か間違いなどありましたら気軽にコメントしてください！
##シリーズ構成
####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法による数値計算法
##今回扱う内容
今回は、輸送速度の符号がわからない場合のスカラー移流方程式の数値計算法について考えていきます。基本的な考え方を扱った後、具体的にプログラムに落とし込んでいきます。
##基本的な考え方
輸送方程式の輸送速度があらかじめわからないとき、以下の分岐をプログラムに組み込む必要があります。
\left \{
\begin{array}{l}
q^{n+1}_{j}=q^{n}_{j}-c\Delta{t}\frac{q^{n}_{j}-q^{n}_{j-1}}{\Delta{x}} (c>0) \\
q^{n+1}_{j}=q^{n}_{j}-c\Delta{t}\frac{q^{n}_{j+1}-q^{n}_{j}}{\Delta{x}} (c<0) \\
\end{array}
\right.
これをプログラムに落とし込むにあたって、If分の分岐を考えると思います。しかし、If文を使ってしまうと計算の過程でコンピュータに負担がかかってしまうのです。なので、以下のような評価式を使うことで解決します。ここでは、一次精度風上法を用いています。
q^{n+1}_{j}=q^{n}_{j}-\Delta{t}\Big(\frac{c+|c|}{2}\frac{q^{n}_{j}-q^{n}_{j-1}}{\Delta{x}}+\frac{c-|c|}{2}\frac{q^{n}_{j+1}-q^{n}_{j}}{\Delta{x}}\Big)
##プログラムへの実装
import numpy as np
import matplotlib.pyplot as plt
#初期値の設定、空間初期分布の設定
c = -1
dt = 0.05
dx = 0.1
jmax = 21
nmax = 6
x = np.linspace(0, dx*(jmax-1), jmax)
q = np.zeros(jmax)
for j in range(jmax):
if (j < jmax/2):
q[j] = 0
else:
q[j] = 1
def do_computing(x, q, c, dt, namax, dx):
plt.figure(figsize=(7,7), dpi=100)
plt.rcParams["font.size"] = 22
#初期分布の可視化
plt.plot(x, q, marker='o', lw=2, label='n=0')
for n in range(1, nmax+1):
qold = q.copy()
for j in range(1, jmax-1):
c2 = (c + abs(c)) / 2
c3 = (c - abs(c)) / 2
q[j] = qold[j] - dt * (c2 * (qold[j] - qold[j-1]) / dx + c3 * (qold[j+1] - qold[j]) / dx)
#各ステップの可視化
if n % 2 == 0:
plt.plot(x, q, marker='o', lw=2, label=f'n={n}')
#グラフの後処理
plt.grid(color="black", linestyle='dashed', linewidth='0.5')
plt.xlim([0,2.0])
plt.ylim([0,1.2])
plt.xlabel('x')
plt.ylabel('q')
plt.legend()
plt.show()
do_computing(x, q, c, dt, nmax, dx)
プログラム結果は、
しっかり計算結果が表れていることがわかりますね。なお、一次精度風上法における数値計算は[前回の記事](https://qiita.com/vardia/items/3bc2b21e75a4b7557897)で詳しく解説しています。
##まとめ
今回は、輸送速度の符号がわからないときの移流方程式の数値計算法を考えました。If文を使うのではなくて、絶対値を使った処理をすることでコンピュータに負担をかけずに計算を行うことができます。
##参考文献
藤井孝蔵、立川智章、Pythonで学ぶ流体力学の数値計