##はじめに
　今回は、第一章の最後としてより実践的な計算法について考えます。ここで考えるのは、スカラー移流方程式に対するTVD法といわれるものです。ここでは、線形のスカラー移流方程式を考えることにします。

##シリーズ構成
####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法による数値計算法

##今回扱う内容
　今回は、TVD法の中でも、HartenやYeeらによるnon-MUSCL型とvan LeerによるMUSCL型の2つを使って数値計算を行っていきます。これまでの計算方法よりも難易度は上がりますが、圧縮性の流れの計算において重要な基礎となる部分になります。

##non-MUSCL型TVD法(c>0の場合)
　第3回で考えた数値流束を用いたLax-Wendroff法を再掲します。

\tilde{f}^{n}_{j+1/2}=\frac{1}{2}\Big\{(1-\nu)f^{n}_{j+1}+(1+\nu)f^{n}_{j}\Big\}\tag{1}

これを変形して、(1次精度風上法の数値流束)+(おつりの項)と変形してみると、

\tilde{f}^{n}_{j+1/2}=f^{n}_{j}+\frac{1}{2}(1-\nu)(f^{n}_{j+1}-f^{n}_{j})\tag{2}

このおつりの項を修正流束といいます。
　Lax-Wendroff法が数値的な振動を起こすのは、物理量が急激に変化するところでこの修正流束のせいでした。そこで、物理量が急激に変化するところでは修正流束が0、そのほかのところでは1となるような係数$B_{j+1/2}$を導入すると、

\tilde{f}^{n}_{j+1/2}=f^{n}_{j}+\frac{1}{2}B_{j+1/2}(1-\nu)(f^{n}_{j+1}-f^{n}_{j})\tag{3}

　係数$B_{j+1/2}$を流速の概念から見ると、「1次精度ではすべての流束が左からきていて、2次精度で状況に応じて量を加減しながら少しづつ右からの流束を入れ込んでる」と考えることができます。この係数$B_{j+1/2}$を流束制限係数といいます。また、式(3)のような手法をHartenの流束修正法といいます。また、この係数$B_{j+1/2}$をTVDと呼ばれる条件を満たすように設定する手法をTVD法といいます。

　ここからは、数値振動が起きないような係数$B_{j+1/2}$の制約を考えていきます。式(3)を、数値流束を用いた基本的な差分式

q^{n+1}_{j}=q^{n}_{j}-\frac{\Delta{t}}{\Delta{x}}\Big(\tilde{f}^{n}_{j+1/2}-\tilde{f}^{n}_{j-1/2}\Big)\tag{4}

に代入すると、

\frac{q^{n+1}_{j}-q^{n}_{j}}{q^{n}_{j-1}-q^{n}_{j}}=\nu\Big\{1-\frac{1}{2}(1-\nu)B_{j-1/2}\Big\}+\frac{1}{2}\nu(1-\nu)\frac{B_{j+1/2}}{r_j}\tag{5}

となります。ここで$r_j$は、

r_j\equiv\frac{q^n_j-q^n_{j-1}}{q^n_{j+1}-q^n_j}

と定義できます。
　さて、数値振動が起きないように下の図に示すように$q^{n+1}_j$が$q^n_j$と$q^n_j-1$の間をとるように制限を加えてみましょう。
　　　　　　　　　　　
したがって、

0\leqq\frac{q^{n+1}_{j}-q^{n}_{j}}{q^{n}_{j-1}-q^{n}_{j}}\leqq1\tag{6}

式(5)の右辺を代入すると、

0\leqq\nu\Big\{1-\frac{1}{2}(1-\nu)B_{j-1/2}\Big\}+\frac{1}{2}\nu(1-\nu)\frac{B_{j+1/2}}{r_j}\leqq1

\Leftrightarrow-\frac{2}{\nu}\leqq B_{j-1/2}-\frac{B_{j+1/2}}{r_j}\leqq\frac{2}{1-\nu}\tag{7}

CFL条件より、$0\leqq\nu\leqq1$であることから、式(7)が変形できて、

-2\leqq B_{j-1/2}-\frac{B_{j+1/2}}{r_j}\leqq2\tag{8}

この式は、

0\leqq B_{j-1/2}\leqq 2 \tag{9}

0\leqq\frac{B_{j+1/2}}{r_j}\leqq2\tag{10}

を共に満たせば成立することがわかります。式(9)と(10)の満たす範囲は下の図の塗りつぶされた部分になります。
　　　　　　　　　
Lax-Wendroff法では、$B_{j+1/2}=1$であるので、$r<1/2$の範囲内では許容範囲外となり数値振動が起こります。このように、図の許容範囲内で流束制限関数が用いられるように、例えば、以下の関数を適用する方法があります。

f(r)=
\left \{
\begin{array}{l}
0　(r<0) \\
r　(0\leqq r\leqq1) \tag{11}\\
1　(r>1)
\end{array}
\right.

##TVD条件
　TVD条件とは、数値振動を抑えるための定量的な条件です。
　1次元の線形スカラー移流方程式について、

U=\int \left|\frac{du}{dx}\right| dx \tag{12}

この量は波の振幅の総和と等しくなります。前節で考えた格子点ごとの変化量に注目して、その総和を以下のように定義します。

TV(u^n)\equiv\displaystyle\sum_j|q^n_{j+1}-q^n_{j}| \tag{13}

このTVをTotal Variationといいます。TVが時間とともに増大していかない条件を考える(数値振動していかないという意味です)と、

TV(u^{n+1})\leqq TV(u^n) \tag{14}

この条件をTVD条件といいます。

##non-MUSCL型TVD法(cの条件なし)
　1次精度風上法を修正して、

\tilde{f}^{n}_{j+1/2}=\frac{1}{2}\Big\{\Big(\bar f^{n}_{j+1}+\bar f^{n}_{j}\Big)-|c|\Big(q^n_{j+1}-q^n_{j}\Big)\Big\}\tag{15}

ここで、
$$\bar f=f+g\tag{16}$$
$$\bar c=c+\gamma\tag{17}$$
$$g_j=minmod(\sigma(c_{j+1/2})\Delta_{j+1/2},(c_{j-1/2})\Delta_{j-1/2})\tag{18}$$
$$\sigma(z)=\frac{1}{2}\Big[\psi(z)-\frac{\Delta{t}}{\Delta{x}}z^2\Big]\tag{19}$$

\left \{
\begin{array}{l}
|z|　(|z|\geqq\delta) \\
\tag{20}
\frac{z^2+\delta^2}{2\delta}　(|z|\leqq\delta) \\
\end{array}
\right.

\left \{
\begin{array}{l}
\frac{g_{j+1}-g_j}{\Delta_{j+1/2}}　(\Delta_{j+1/2}\neq0) \\
\tag{21}
0　(\Delta_{j+1/2}=0) \\
\end{array}
\right.

のように修正することから、式(15)の流束に「-」がついています。ここで、minmod関数というのは以下で定義されます。
$$minmod(x,y)=\Big(sgn(0.5,x)+sgn(0.5,y)\Big)min(|x|,|y|)\tag{22}$$
ここで、sgn関数は符号だけを取り出す関数です。minmod関数の意味は簡単に言うと、「符号が違うと0、同じなら小さいほうを返す関数」といえます。

　また、この方法の発想としては、
「不連続や急こう配のところでは1次精度で計算する。その時は変化があったとしても傾きの小さいほうを採用して計算する。」
という考え方といえます。

##Harten-Yeeのnon-MUSCL型TVD法
　HartenとYeeは前節の修正方法をプログラムで書きやすいように単純化しました。

\tilde{f}^{n}_{j+1/2}=\frac{1}{2}\Big\{\Big(f^{n}_{j+1}+f^{n}_{j}\Big)+\phi^n_{j+1/2}\Big\}\tag{23}

$$\phi^n_{j+1/2}=\sigma(c_{j+1/2})(g_j+g_{j+1})-\psi(c_{j+1/2}+r_{j+1/2})(q^n_{j+1}-q^n_j)\tag{24}$$
$$\sigma(z)=\frac{1}{2}\Big[\psi(z)-\frac{\Delta{t}}{\Delta{x}}z^2\Big]\tag{25}$$
$$minmod(\Delta_{j+1/2},\Delta_{j-1/2})\tag{26}$$

\left \{
\begin{array}{l}
\frac{\sigma(c_{j+1/2})(g_{j+1}-g_j)}{\Delta_{j+1/2}}　(\Delta_{j+1/2}\neq0) \\
\tag{27}
0　(\Delta_{j+1/2}=0) \\
\end{array}
\right.

これをHarten-Yeeの風上型TVD法という。
　では、これを用いて移流方程式を数値計算しましょう。

import numpy as np
import matplotlib.pyplot as plt

def init(q1, q2, dx, jmax, sp=0.5):
    x = np.linspace(0, dx * (jmax-1), jmax)
    q = np.array([(float(q1) if i < sp*jmax else float(q2)) for i in range(jmax)])
    return (x, q)

def minmod(x, y):
    sgn = np.sign(x)
    return sgn * max(min(abs(x), sgn*y), 0.0)

def do_computing2(x, q, c, dt, dx, nmax, ff, interval = 2):
    plt.figure(figsize=(7,7), dpi=100) # グラフのサイズ
    plt.rcParams["font.size"] = 22 # グラフの文字サイズ
    
    # 初期分布
    plt.plot(x, q, marker='o', lw=2, label='n=0') 
    
    delta = np.zeros(jmax)
    g = np.zeros(jmax)
    flux = np.zeros(jmax)
    
    for n in range(1, nmax + 1):
        qold = q.copy()
        
        for j in range(0, jmax - 1):
            delta[j] = qold[j + 1] - qold[j]
        
        for j in range(1, jmax - 1):
            g[j] = minmod(delta[j], delta[j-1])
        
        for j in range(0, jmax - 1):
            flux[j] = ff(qold, c, delta, g, dt, dx, j)
        
        for j in range(1, jmax - 1):
            q[j] = qold[j] - dt / dx * (flux[j] - flux[j-1])
            
        # 各ステップの可視化
        if n % interval == 0:
            plt.plot(x, q, marker='o', lw=2, label=f'n={n}')

    # グラフの後処理
    plt.grid(color='black', linestyle='dashed', linewidth=0.5)
    plt.xlabel('x')
    plt.ylabel('q')
    plt.legend()
    plt.show()

def TVD(q, c, delta, g, dt, dx, j):
    sigma = 0.5 * (np.abs(c) - dt/dx * c**2)
    gamma = sigma * (g[j+1] - g[j]) * delta[j]/(delta[j]**2 + 1e-12)
    phi = sigma * (g[j] + g[j+1]) - abs(c + gamma) * delta[j]

    ur = q[j+1]
    ul = q[j]
    fr = c * ur
    fl = c * ul
    return 0.5 * (fr + fl + phi)

c = 1
dt = 0.05
dx = 0.1

jmax = 21
nmax = 20
q1 = 1
q2 = 0
x, q = init(q1, q2, dx, jmax, sp=0.05)

do_computing2(x, q, c, dt, dx, nmax, TVD, interval=4)

実行すると、
　　　　　　

##MUSCL型TVD法
　MUSCL型TVD法は、van LeerのMUSCL内挿という方法を利用します。すでに流束や物理量の内挿の仕方については考えてきましたので、ここからは、非線形問題を考えていきます。
　MUSCL法は、本来は有限体積法に基づく方法であるので、セル平均を基本として計算が行われることになります。各セルを代表する物理量の分布を以下に示します。一次精度の場合、セル内で物理量が一定であると仮定します。2次精度の場合は、各セル内での物理量が線形であると仮定します。これを考えると、3次精度にするには、セル内の物理量を2次関数で仮定するれ場良いことがわかると思います。

　1次精度セル分布
　　　　　　　

　2次精度セル分布
　　　　　　　
　さて、セル内の物理量は一般的に、

q(x)=q_j+\frac{1}{\Delta{x}}(x-x_j)\delta_jq+\frac{3\kappa}{2(\Delta{x})}\Big[(x-x_j)^2-\frac{(\Delta{x})^3}{12}\Big]\delta_jq　(x_{j-1/2}\leqq x \leqq x_{j+1/2})\tag{28}

ただし、

q_j=\frac{1}{\Delta{x}}\int_{j-1/2}^{j+1/2} q(x) dx\tag{29}

と定義します。
　セル境界での物理量は、

(q_L)_{j+1/2}=q_j+\frac{1}{4}(1-\kappa)(q_j-q_{j-1})+\frac{1}{4}(1+\kappa)(q_{j+1}-q_j)\tag{30}

(q_R)_{j+1/2}=q_j-\frac{1}{4}(1+\kappa)(q_{j+1}-q_j)-\frac{1}{4}(1-\kappa)(q_{j+2}-q_{j+1})\tag{31}

となります。添え字の$L$は「左から」を意味し、$R$は「右から」を意味します。いずれも第1項が1次精度を表していて、第2項を加えることにより高精度化できます。プリグラムで書きやすいように書き直すと、

(q_L)_{j+1/2}=q_j+\frac{\epsilon}{4}\Big[(1-\kappa)(q_j-q_{j-1})+(1+\kappa)(q_{j+1}-q_j)\Big]\tag{30'}

(q_R)_{j+1/2}=q_j-\frac{\epsilon}{4}\Big[(1+\kappa)(q_{j+1}-q_j)-(1-\kappa)(q_{j+2}-q_{j+1})\Big]\tag{31'}

$\epsilon=0$のとき1次精度、$\epsilon=1$のときは高次精度を表します。$\epsilon=1$かつ$\kappa=-1$とすると、2次精度の完全風上法を表します。$\kappa=0$のときは風下側が1点加わる風上法で、$\kappa=1/3$のときは、3次精度となります。
　このように、物理量の内挿、外挿の係数を選ぶことにより、1次～3次までの精度をもつ方法が定義できます。今後の変形が簡単になるように、デルタ形式を導入したいと思います。

(\Delta_{+})_{j}=q_{j+1}-q_j,　(\Delta_{-})_{j}=q_{j}-q_{j-1} \tag{32}

式(26')(27')をデルタ形式で書き直しておきましょう。

(q_L)_{j+1/2}=q_j+\frac{1}{4}(1-\kappa)(\Delta_{-})_{j}+\frac{1}{4}(1+\kappa)(\Delta_{+})_{j}\tag{33}

(q_R)_{j+1/2}=q_j-\frac{1}{4}(1+\kappa)(\Delta_{-})_{j+1}-\frac{1}{4}(1-\kappa)(\Delta_{+})_{j+1}\tag{34}

TVD法にするために制御関数を入れましょう。

(q_L)_{j+1/2}=q_j+\frac{1}{4}(1-\kappa)(\bar\Delta_{-})_{j}+\frac{1}{4}(1+\kappa)(\bar\Delta_{+})_{j}\tag{35}

(q_R)_{j+1/2}=q_j-\frac{1}{4}(1+\kappa)(\bar\Delta_{-})_{j+1}-\frac{1}{4}(1-\kappa)(\bar \Delta_{+})_{j+1}\tag{36}

ここで、

\bar \Delta_{+}=minmod(\Delta_{+},b\Delta_{-}) \\
\bar \Delta_{-}=minmod(\Delta_{-},b\Delta_{+}) \\
b=\frac{3-\kappa}{1-\kappa}

と定義しています。
　いくつかの$\kappa$をイメージして常識を見てみると、「変化量の～倍と比較して小さいほうを利用している」ことがわかると思います。

##まとめ
　本記事では、TVD法と呼ばれる実践的な計算法を紹介しました。Hartenの流束修正法はシステム方程式に対する離散化の基礎となっている概念でありとても重要です。
　第1回～第6回では、スカラー移流方程式の数値計算法について考えました。次回からは、引き続きスカラー移流方程式を扱いますが、その中でも時間積分の方法に注目していきたいと思います。

　長々と書いてきたので、どこかで間違いがあるかもしれません。そのときはお知らせくださいませ。。

##参考文献
藤井孝蔵、立川智章、Pythonで学ぶ流体力学の数値計算法、2020/10/20