数値流体力学8～時間陰解法～

##はじめに
　前回は、時間陽解法について考えました。今回は、それと対をなす時間陰解法について考えていきたいと思います。前回と同様、基本的な考え方を抑えた後、主な計算方法について扱っていきます。
　時間陽解法の記事で扱った置換を本記事でも引き続き使っていきます。
　もちろん、今回も扱う題材はスカラー移流方程式です。

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

##時間陰解法の考え方
　**離散化式全体の中で時間n+1の項が複数存在する計算法を時間陰解法といいます。**時間陰解法の重要な性質は、時間ステップnの任意の格子点の物理量の変化が時間ステップn+1のすべての物理量に影響を及ぼすことです。つまり、物理量が極めて速く伝播するのです。これにより、時間陰解法では一般的に、CFL条件を考慮しなくてよいことになっています。

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

##主な時間陰解法
####Crank-Nicholson法
　空間微分に2次精度の中心差分を適用した離散化方法をCrank-Nicholson法といいます。

\frac{q^{n+1}_j-q^n_j}{\Delta{t}}=\frac{1}{2}c\Big\{\frac{q^{n+1}_{j+1}-q^{n+1}_{j-1}}{2\Delta{x}}-\frac{q^{n}_{j+1}-q^{n}_{j-1}}{2\Delta{x}}\Big\}=0\tag{2}

クーラン数を$\nu$と書き換えて、時間ステップn+1の項を左辺に、時間ステップnの項を右辺にまとめると、

-\frac{1}{4}\nu q^{n+1}_{j-1}+q^{n+1}_{j}+\frac{1}{4}\nu q^{n+1}_{j+1}=\frac{1}{4}\nu q^{n}_{j-1}+q^{n}_{j}-\frac{1}{4}\nu q^{n}_{j+1}\tag{3}

これは、行列式としてまとめることができます。

\left[
\begin{array}{ccc}
1 & \frac{\nu}{4} & \cdots & \cdots & 0 \\
-\frac{\nu}{4} & \ddots & & & \vdots \\
\vdots & &\ddots & &\vdots \\
\vdots &  &  &\ddots & \frac{\nu}{4}\\
0 & \cdots &\cdots & -\frac{\nu}{4} & 1
\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{4}

ここで、RHSはRight Head Sideの略で時間ステップnで評価される式(4)の右辺をまとめたものです。知りたい情報は左辺の第2項ですから、左辺第1項の行列を反転させなければいけないのです。このように、時間陰解法は演算量が増加する問題があります。
　当然、陰解法の性質である、時間ステップnの任意の格子点の物理量の変化が時間ステップn+1のすべての物理量に影響を及ぼしていることも式(4)からわかります。

　von Meumann解析でも安定性について確認しておきましょう。式(2)から、

g=\frac{1-\frac{1}{2}i\nu sin\theta}{1+\frac{1}{2}i\nu sin\theta}\tag{5}

つまり、

\left \{
\begin{array}{l}
|g|=1 \\
\tag{6}
\phi=-tan^{-1}\frac{\nu sin\theta}{1+\frac{1}{4}\nu^2 sin^2\theta}
\end{array}
\right.

$|g|=1$であることから、クーラン数がどのような値であっても中立安定であることがわかります。しかし、Crank-Nicholson法にも弱点が2つあります。まずは前述した計算量の問題、それから、$\nu$が大きくなると、位相誤差の観点から精度が悪くなってしまいます。
　Crank-Nicholson法を使って移流方程式を数値計算していきましょう。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

jmax = 21
nmax = 100
dx = 0.1
dt = 0.05

c = 0.1

qq = np.zeros([jmax, nmax])

#初期条件の設定
for j in range(0,jmax):
    if (j < jmax/2-1):
        qq[j,:] = 1.0
    
    else:
        qq[j,:] = 0.0

#輸送速度の行列を作成    
coe = np.ones([jmax, nmax])
for i in range(jmax):
    coe[i,:] = c

#クーラン数をもつ行列の作成
q = np.zeros([jmax, nmax])
nu = np.ones([jmax, nmax])
for i in range(jmax):
    nu[i,:] = dt * coe[i,:] / dx

#クランクニコルソン法の適用
for n in range(nmax-1):
    #連立方程式の係数行列Aの作成
    A = np.zeros([jmax-2,jmax-2])
    for i in range(jmax-2):
        A[i,i] = 4/nu[i,n]
        if i >=1 :
            A[i-1,i] = 1
        if i <= jmax-4 :
            A[i+1,i] = -1
    
    #係数行列bの作成
    b = np.zeros([jmax-2])
    for i in range(jmax-2):
        b[i] = qq[i,n] + 4/nu[i,n] * qq[i+1,n] - qq[i+2,n]
    b[0] += qq[0,n+1]
    b[jmax-3] += qq[-1,n+1]
    
    #Aq=bをqについて解く
    q = np.linalg.solve(A,b)
    for i in range(jmax-2):
        qq[i+1,n+1] = q[i]
        
fig = plt.figure()

ims = [] 
for i in range(nmax):
    T = list(qq[:,i])
    x = np.linspace(0,dx*(jmax-1),jmax)
    if i % int(nmax*0.02) == 0:
        ##im=plt.plot(x,T, '-', color='red',markersize=10, linewidth = 2, aa=True)
        ##anim.append(im)
        
        if i % 20 == 0:
            im=plt.plot(x,T,marker='o', lw=2, label=f'n={i}')

##ani = animation.ArtistAnimation(fig, ims)
plt.xlabel('x')
plt.ylabel('q')
plt.legend()
plt.show()

実行すると、
　　　　　　　
位相誤差に起因する数値振動がみられると思います。この解析結果からもわかる通り、位相誤差による振動はすべて波の後ろで起こっていることもわかります。アニメーションにするとこれがより分かりやすいです。

####近似LDU分解
　まずは非線形方程式について考えていきます。流速関数については時間方向にテイラー展開ができて、

f^{n+1}=f^n+\left.\frac{\partial{f}}{\partial{q}}\right|^n(\Delta{q^n})+\left.\frac{\partial^2{f}}{\partial{q^2}}\right|^n(\Delta{q^n})^2+･･･\tag{6}

ここでは時間方向と書いていますが、式(6)の展開は$\Delta{q}$でやっています。保存型のスカラー移流方程式に対してオイラー陰解法を適用した形は以下のように近似できます。

\Delta q^n_j=-\Delta{t}\Big(\frac{\partial{f}}{\partial{x}}\Big)^{n+1}_j
=-\Delta{t}\frac{\partial}{\partial{x}}\Big[f^n+\left.\frac{\partial{f}}{\partial{q}}\right|^n(\Delta{q^n})+\left.\frac{\partial^2{f}}{\partial{q^2}}\right|^n(\Delta{q^n})^2+･･･\Big] \tag{7}

時間($\Delta{q}$)の2次精度をとると、

\bigg(1+\Delta{t}\frac{\partial}{\partial{x}}\Big(\frac{\partial{f}}{\partial{q}}\Big)\bigg)^n_j \Delta{q^n_j}=-\Delta{t}\left.\frac{\partial{f}}{\partial{x}} \right|^n_j\tag{8}

さらに、$\partial{f}/\partial{q}=c(const)$のとき、

\Big(1+\Delta{t}\frac{\partial{c^+}}{\partial{x}}+\Delta{t}\frac{\partial{c^-}}{\partial{x}}\Big)=-\Delta{t}R^n(q_j)\tag{9}

ここで、式(9)の左辺を近似化すると、

\Big(1+\Delta{t}\frac{\partial{c^+}}{\partial{x}}+\Delta{t}\frac{\partial{c^-}}{\partial{x}}\Big)\approx \Big(1-\frac{\Delta{t}}{\Delta{x}}c^-_j+\Delta \delta^b_xc^+\Big)\Big(1+\frac{\Delta{t}}{\Delta{x}}|c| \Big)^{-1}\tag{10}

この右辺を3段階のスイープで計算していくことになります。

【第1スイープ】

\Big(1+\frac{\Delta{t}}{\Delta{x}}|c|_j\Big)^n_j\Delta{q^*_j}=-\Delta{t}R^n(q_j)+\Big(\frac{\Delta{t}}{\Delta{x}}|c|_j\Big)^n_j\Delta{q^*_{j-1}}\tag{11}

【第2スイープ】

\Delta{q^{**}_j}=\Big(1+\frac{\Delta{t}}{\Delta{x}}|c|_j\Big)^n_j\Delta{q^*_j}\tag{12}

【第3スイープ】

\Big(1+\frac{\Delta{t}}{\Delta{x}}|c|_j\Big)\Delta{q^n_j}=\Delta{q^{**}_j}-\Big(\frac{\Delta{t}}{\Delta{x}}c^-_{j+1}\Big)\Delta{q^*_{j+1}}\tag{13}

このLDU分解を移流方程式に適用してみましょう。

import numpy as np
import matplotlib.pyplot as plt

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

def UPWIND1(alf, q, c, dt, dx, jmax):
    for j in range(0, jmax-1):
        ur, ul = q[j+1], q[j]
        fr, fl = c * ur, c * ul
        alf[j] = 0.5 * (fr + fl - abs(c) * (ur-ul))

def do_computing_LDU(x, q, c, dt, dx, nmax, ff, interval=2, xlim=None):
    plt.figure(figsize=(7,7), dpi=100)
    plt.rcParams["font.size"] = 22

    plt.plot(x, q, marker='o', lw=2, label='t=0')

    alf = np.zeros(jmax)
    dq = np.zeros(jmax)
    for  n in range(1, nmax-1):
        qold = q.copy()

        #近似LUD分解
        c_a = abs(c)
        c_p = 0.5 * (c + c_a)
        c_n = 0.5 * (c - c_a)
        nu_a = c_a * dt/dx
        nu_p = c_p * dt/dx
        nu_n = c_n * dt/dx

        ff(alf, qold, c, dt, dx, jmax)
        R = np.append(0.0, np.diff(alf) / dx)

        ##第1スイープ
        for j in range(1, jmax-1):
            dq[j] = (-dt * R[j] + nu_p * dq[j-1]) / (1 + nu_a)

        ##第2スイープ
        for j in range(jmax-2, 0, 1):
            dq[j] = dq[j] - nu_n * dq[j+1] / (1 + nu_a)
        
        ##第3スイープ
        for j in range(1, jmax-1):
            q[j] = qold[j] + dq[j]

        #各ステップの可視化
        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')
    if xlim is not None:
        plt.xlim(xlim)
    plt.legend()
    plt.show()

c = 1
dt = 0.8
dx = 0.1

jmax = 70
nmax = 5

xs = -1

q1 = 1
q2 = 0
x, q = init(q1, q2, xs, dx, jmax)
do_computing_LDU(x,q, c, dt, dx, nmax, UPWIND1, interval=1, xlim=[-1,5])

時間刻み幅$\Delta{t}$を色々変えると、
①$\Delta{t}=0.05,\Delta{x}=0.1$
　　　　　　　

②$\Delta{t}=0.1,\Delta{x}=0.1$
　　　　　　　

③$\Delta{t}=0.2,\Delta{x}=0.1$
　　　　　　　
④$\Delta{t}=0.8,\Delta{x}=0.1$
　　　　　　　
どのケースでも時間積分が安定して行われていることがわかります、一方、クーラン数が$2.0$程度までであればそこまで変わらないものの、④のケースでは、かなり解が鈍ってしまっていることがわかりますね。これは、時間陰解法による時間的な離散化誤差の蓄積が原因となります。

##まとめ
　今回は、時間陰解法について扱いました。代表的なCrank-Nicholson法やLDU分解を用いることによってその性質について議論していきました。
　次回からは、2階偏微分方程式について数値計算を扱っていきたいと思います。

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