はじめに

　今回から、移流方程式の時間積分に注目していきたいと思います。まずは、時間陽解法について本記事で整理したいと思います。

シリーズ構成

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法による数値計算法

今回扱う内容

　時間積分法として、時間陽解法の考え方とその主な手法について整理していきます。

時間陽解法の考え方

　時間陽解法については今までの記事でも扱ってきました。数値流束を用いた移流方程式の離散化式は以下の式(1)のようでした。

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{1}

この空間微分項を$R_j$と書くことで、

q^{n+1}_{j}=q^{n}_{j}-\Delta{t}R_j(q^n) \tag{2}

この式のように、$R_j$がすべて時間ステップnで評価され、離散化式全体で時間ステップn+1の項が一つしかない時間積分法を時間陽解法といいます。

主な時間陽解法

　これまでも見てきた通り、基本的には時間陽解法にはCFL条件があるので、その条件に引っかからないように、多段階の時間陽解法が開発され、主流となっています。

2段階時間陽解法

\left \{
\begin{array}{l}
q^*_j=q^n_j-\Delta{t}R_j(q^n) \\
\tag{3} \\
q^{n+1}_{j}=q^{n}_{j}-\frac{1}{2}\Delta{t}[R_j(q^n)+R_j(q^*)]
\end{array}
\right.

2段階ルンゲクッタ法

　実用的によく使われているのはルンゲクッタ法である。ルンゲクッタ法は、空間と時間を別々に考えて、状微分方程式で使用される方法を適用することができるのです。

\left \{
\begin{array}{l}
q^{(1)}_j=q^n_j-\frac{1}{2}\Delta{t}R_j(q^n) \\
\tag{4} \\
q^{n+1}_{j}=q^{n}_{j}-\Delta{t}R_j(q^{(1)})
\end{array}
\right.

多段階ルンゲクッタ法

　ルンゲクッタ法は2段階に限らず、多段階の方法を用いることによって、時間積分の安定領域を広げることができます。以下には、4段階のルンゲクッタ法を示しておきましょう。

\left \{
\begin{array}{1}
q^{(1)}_j=q^n_j \\ \\
q^{(2)}_j=q^n_j-\frac{1}{2}\Delta{t}R_j(q^{(1)}) \\ \\
q^{(3)}_j=q^n_j-\frac{1}{2}\Delta{t}R_j(q^{(2)}) \tag{5}\\ \\
q^{(4)}_j=q^n_j-\Delta{t}R_j(q^{(3)}) \\ \\
q^{n+1}_j=q^n_j+\frac{1}{6}\Delta{t}[R_j(q^{(1)})+2R_j(q^{(2)})+2R_j(q^{(3)})+R_j(q^{(4)})]
\end{array}
\right.

多段階になると、安定性が向上しますが、計算負荷がその分増えてしまうことは言うまでもありません。

2段階Lax-Wendroff法

　第1回から扱ってきたLax-Wendroff法は、スカラー移流方程式では多少の数値振動が起こる程度で問題ないのですが、圧縮性オイラー方程式や圧縮性ナビエ-ストークス方程式を解くには計算負荷が多くなり、面倒です。そこで、Lax-Wendroff法にも、2段階に分けることで演算量を減らす方法がとられています。

\left \{
\begin{array}{1}
q^{n+1/2}=\frac{1}{2}(q^{n}_{j+1}+q^{n}_{j})-\frac{1}{2}c\Delta{t}\frac{q^{n}_{j+1}-q^{n}_{j}}{\Delta{x}} \\
\tag{6} \\
q^{n+1}_{j}=q^{n}_{j}-c\Delta{t}\frac{q^{n+1/2}_{j+1/2}-q^{n+1/2}_{j-1/2}}{\Delta{x}}
\end{array}
\right.

2段階のLax-Wendroff法のアイデアは、下の図を見てみるとわかると思います。
　　　　　　　　　
1段階目で時間ステップ$n$での3点を使って時間ステップ$n+1/2$の2点を作り、2段階目で時間ステップ$n+1/2$の2点を使って時間ステップ$n+1$の点を作るといった感じです。

MacCormack法

2段階Lax-Wendroff法をさらに簡単にしたのがMacCormack法です。

\left \{
\begin{array}{1}
\bar q_j=q^n_j-c\Delta{t}\frac{q^{n}_{j+1}-q^{n}_{j}}{\Delta{x}} \\
\tag{7} \\
q^{n+1}_j=\frac{1}{2}(q^{n}_{j}+\bar q_{j})-c\Delta{t}\frac{\bar q_j-\bar q_{j-1}}{\Delta{x}}
\end{array}
\right.

pythonで計算してみましょう。

import numpy as np
import matplotlib.pyplot as plt

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

def do_computing_mc(x, q, c, dt, dx, nmax, 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')

    for n in range(1, nmax + 1):
        qbar = q.copy()
        
        # MacCormack法
        for j in range(0, jmax-1):
            qbar[j] = q[j] - dt * c * (q[j+1] - q[j]) / dx
        
        for j in range(1, jmax-1):
            q[j] = 0.5 * ((q[j] + qbar[j]) - dt * c * (qbar[j] - qbar[j-1]) / dx)
        
        # 各ステップの可視化
        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.xlim([0, 2])
    plt.ylim([0,1.2])
    plt.xlabel('x')
    plt.ylabel('q')
    plt.legend()
    plt.show()

c = 1
dt = 0.05
dx = 0.1

jmax = 21
nmax = 6
q1 = 1
q2 = 0
q, x = init(q1, q2, dx, jmax)

do_computing_mc(x, q, c, dt, dx, nmax, interval=2)

実行すると、当然Lax-Wendroff法と同じ結果が得られるはずです。
　　　　　　　　

まとめ

　今回は時間陽解法の基本的な考え方と、その主な計算方法について考えました。多段階にすることにより、CFL条件など安定性の条件をうまく避けるような方法が多かったかと思います。一方で、多段階にすることによって計算負荷がどんどん大きくなってしまうジレンマがあるのも忘れてはいけません。
　次回は、時間陰解法について考えていきたいと思います。

参考文献

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