非圧縮性粘性流れの数値計算プログラムにおける時間積分精度を調査する

概要

非圧縮性粘性流れの数値計算プログラムにおける，時間積分精度を確認するための一つの方法を示す．

非圧縮性粘性流れに限らず，空間と時間を離散化する問題で時間精度を確認するには，空間離散化に起因する誤差を極力小さくする必要がある．

時間積分法として，1次の前進Euler法，Sanderse-Korenによる3次のRunge-Kutta法¹，Aithal-FerranteによるFastRK3²³を用いる．

支配方程式および数値計算法

流れに非圧縮性を仮定すると，質量および運動量は，それぞれ式(1), (2)で表される．

\begin{align}
    \nabla \cdot \boldsymbol{u} &= 0\qquad (1)\\
    \frac{\partial \boldsymbol{u}}{\partial t} + \left( \boldsymbol{u}\cdot \nabla \right)\boldsymbol{u} &= -\frac{1}{\rho}\nabla p + \nu \nabla^2 \boldsymbol{u}\qquad (2)
\end{align}

ここで，$\boldsymbol{u}$は速度，$p$は圧力，$\rho$は密度，$\nu$は動粘度である．

$\boldsymbol{u}$と$p$の時間変化を計算するために，部分段階法を適用する．部分段階法では，まず式(2)中の圧力勾配の項を無視して速度の時間発展を計算した後，圧力勾配の影響による速度の時間変化を計算する．

まず，式(2)の時間微分項を前進Euler法によって離散化する．上付き添字$n$は時刻を表している．圧力は時刻$n+1$の値を，速度は時刻$n$の値を用いると，次式が得られる．

    \frac{\boldsymbol{u}^{n+1} - \boldsymbol{u}^n}{\Delta t} + \left( \boldsymbol{u}^n\cdot \nabla \right)\boldsymbol{u}^n = -\frac{1}{\rho}\nabla p^{n+1} + \nu \nabla^2 \boldsymbol{u}^n \qquad(3)

ここで，$p^{n+1}$は直ちには判らないので，圧力を無視して式(3)を計算する．

    \frac{\boldsymbol{u}^{*} - \boldsymbol{u}^n}{\Delta t} + \left( \boldsymbol{u}^n\cdot \nabla \right)\boldsymbol{u}^n = \nu \nabla^2 \boldsymbol{u}^n \qquad(4)

圧力の影響を無視しているので，得られる速度は$\boldsymbol{u}^{n+1} $ではない．便宜上求めるこの速度$\boldsymbol{u}^{*}$を，仮速度や中間速度とよぶ．式(3)と式(4)の差を取ると，圧力勾配項の影響による速度の時間変化を記述する式が得られる．

    \frac{\boldsymbol{u}^{n+1} - \boldsymbol{u}^*}{\Delta t}  = -\frac{1}{\rho}\nabla p^{n+1} \qquad(5)

速度の時間発展を，式(4),(5)に分けることが，部分段階法の名前の由来である．

$p^{n+1}$を決定するために，式(5)の発散を取る．

    \frac{1}{\Delta t}(\nabla\cdot\boldsymbol{u}^{n+1} - \nabla\cdot\boldsymbol{u}^*)  = -\frac{1}{\rho}\nabla^2 p^{n+1} \qquad(6)

$\boldsymbol{u}^*$は圧力の影響が入っていないので連続の式は満たさないが，$\boldsymbol{u}^{n+1}$は連続の式を満たさなければならない．すなわち，$\nabla\cdot\boldsymbol{u}^{n+1}=0$である必要がある．そのためには，式(6)において

    \frac{1}{\Delta t}\nabla\cdot\boldsymbol{u}^* - \frac{1}{\rho}\nabla^2 p^{n+1} =0

でなければならない．速度が関係する項を右辺に，圧力を左辺に置くと，これは速度と圧力の関係を表すPoisson方程式となる．

     \nabla^2 p^{n+1}=\frac{\rho}{\Delta t}\nabla\cdot\boldsymbol{u}^*\qquad(7)

これを圧力Poisson方程式とよぶ．

部分段階法では，$\boldsymbol{u}^n$が既知の時，次のような手順で流れ場の時間発展を計算する．

1. 式(4)を用いて$\boldsymbol{u}^n$から中間速度$\boldsymbol{u}^*$を計算する．
2. 式(7)を解いて，圧力$p^{n+1}$を計算する
3. 式(5)を用いて$\boldsymbol{u}^{n+1}$を計算する．

時間積分の高次精度化

上述のアルゴリズムは，時間積分に1次の前進Euler法を用いているため，形式的には時間1次精度（$\Delta t$を$1/10$にしても誤差は$1/10$にしかならない）しかない．時間積分に起因する誤差を低減するのが難しく，計算結果が理論値や実験値から乖離する原因となる．また，空間に2次の中心差分，時間に1次の前進Euler法を用いるいわゆるFTCSスキームは，全ての$\Delta t$に対して不安定で，長時間計算を続けると解が発散に至る．特に粘性よりも慣性の影響が大きいような流れを計算する場合には注意が必要である．

これらの理由から，非圧縮性流れの数値計算では少なくとも2次精度以上の時間積分法が用いられる．時間積分に2次精度のAdams-Bashforth法とCrank-Nicolson法の組合せや，3次精度Runge-Kutta法を見かけることが多い．

SKRK3

Sanderse and Koren¹は，3次精度のRunge-Kutta法（以降，SKRK3と記述する）を用いた非圧縮性流れの数値計算法を示している．SKRK3でもFractional Step法が用いられるが，Runge-Kutta法の各段階で計算される中間速度の式が異なる．SKRK3では，各段階$s$において，以下の式で中間速度を計算する．

{\boldsymbol{u}^{*}}^{(s)} = \boldsymbol{u}^n+\Delta t\sum_{i=1}^s a_{s,i}F(\boldsymbol{u}^{(s-1)})

$F(\boldsymbol{u})=-\left( \boldsymbol{u}\cdot \nabla \right)\boldsymbol{u} + \nu \nabla^2 \boldsymbol{u}$である．

中間速度を計算した後，各段階での圧力$p^{(s)}$の勾配を利用して式(1)を満たす速度を計算する．

\boldsymbol{u}^{(s)} = {\boldsymbol{u}^*}^{(s)} -c_s\Delta t\frac{1}{\rho}\nabla p^{(s)} 

式(6)から式(7)を求めたのと同じ理屈を使えば，上式から圧力Poisson方程式が得られる．

\nabla^2 p^{(s)}=\frac{\rho}{c_s\Delta t}\nabla\cdot{\boldsymbol{u}^*}^{(s)}

SKRK3を用いると，次のような手順で流れ場の時間発展が計算される．

1. $\boldsymbol{u}^{(0)}:=\boldsymbol{u}^{n}$とする．
2. ${\boldsymbol{u}^{*}}^{(1)} = \boldsymbol{u}^n+\Delta t a_{1,1}F(\boldsymbol{u}^{(0)})$を計算する．
3. $\nabla^2 p^{(1)}=\frac{\rho}{c_1\Delta t}\nabla\cdot{\boldsymbol{u}^*}^{(1)}$を解いて，圧力$p^{(1)}$を得る．
4. $\boldsymbol{u}^{(1)} = {\boldsymbol{u}^*}^{(1)} -c_1\Delta t\frac{1}{\rho}\nabla p^{(1)}$を計算して第1段階の速度$\boldsymbol{u}^{(1)}$を得る．
5. ${\boldsymbol{u}^{*}}^{(2)} = \boldsymbol{u}^n+\Delta t \left[a_{2,1}F(\boldsymbol{u}^{(0)})+a_{2,2}F(\boldsymbol{u}^{(1)})\right]$を計算する．
6. $\nabla^2 p^{(2)}=\frac{\rho}{c_2\Delta t}\nabla\cdot{\boldsymbol{u}^*}^{(2)}$を解いて，圧力$p^{(2)}$を得る．
7. $\boldsymbol{u}^{(2)} = {\boldsymbol{u}^*}^{(2)} -c_2\Delta t\frac{1}{\rho}\nabla p^{(2)}$を計算して第2段階の速度$\boldsymbol{u}^{(2)}$を得る．
8. ${\boldsymbol{u}^{*}}^{(3)} = \boldsymbol{u}^n+\Delta t \left[a_{3,1}F(\boldsymbol{u}^{(0)})+a_{3,2}F(\boldsymbol{u}^{(1)})+a_{3,3}F(\boldsymbol{u}^{(2)})\right]$を計算する．
9. $\nabla^2 p^{(3)}=\frac{\rho}{c_3\Delta t}\nabla\cdot{\boldsymbol{u}^*}^{(3)}$を解いて，圧力$p^{(3)}$を得る．
10. $\boldsymbol{u}^{(3)} = {\boldsymbol{u}^*}^{(3)} -c_3\Delta t\frac{1}{\rho}\nabla p^{(3)}$を計算して第3段階の速度$\boldsymbol{u}^{(3)}$を得る．
11. $\boldsymbol{u}^{n+1}:=\boldsymbol{u}^{(3)}, p^{n+1}:=p^{(3)}$とする．

Sanderse and Korenは，係数$a_{s,i}, c_s$を次のとおり定めている．

a = \left[
\begin{matrix}
\frac{1}{3}& 0 & 0\\
-1& 2& 0\\
0& \frac{3}{4}& \frac{1}{4}
\end{matrix}
\right]

c = \left[
\begin{matrix}
\frac{1}{3}\\
1\\
1
\end{matrix}
\right]

FastRK3

SKRK3（に限らないが）では，各段階で圧力Poisson方程式を解く必要がある．非圧縮性粘性流れの数値計算では，一般的に圧力Poisson方程式の求解が計算時間の大半を占めるため，可能な限り圧力Poisson方程式は解きたくない．

Aithal and Ferrante²は，SKRK3の第1，第2段階で用いている圧力を過去の圧力から外挿し，圧力Poisson方程式の求解を第3段階での1回に限定するFastRK3を提案している．

FastRK3を用いると，次のような手順で流れ場の時間発展が計算される．

1. 第1，第2段階における擬似的な圧力$\phi^{(1)}, \phi^{(2)}$を，それぞれ$\phi^{(1)}=(1+c_2)p^{n}-c_2p^{n-1}, \phi^{(2)}=(1+c_3)p^{n}-c_3p^{n-1}$によって外挿する．
2. $\boldsymbol{u}^{(0)}:=\boldsymbol{u}^{n}, \tilde{\boldsymbol{u}}^{(0)}:=\boldsymbol{u}^{n}$とする．
3. ${\boldsymbol{u}^{*}}^{(1)} = \boldsymbol{u}^n+\Delta t a_{1,1}F(\tilde{\boldsymbol{u}}^{(0)})$を計算する．
4. $\tilde{\boldsymbol{u}}^{(1)} = {\boldsymbol{u}^*}^{(1)} -c_1\Delta t\frac{1}{\rho}\nabla \phi^{(1)}$を計算して第1段階の近似速度$\tilde{\boldsymbol{u}}^{(1)}$を得る．
5. ${\boldsymbol{u}^{*}}^{(2)} = \boldsymbol{u}^n+\Delta t \left[a_{2,1}F(\tilde{\boldsymbol{u}}^{(0)})+a_{2,2}F(\tilde{\boldsymbol{u}}^{(1)})\right]$を計算する．
6. $\tilde{\boldsymbol{u}}^{(2)} = {\boldsymbol{u}^*}^{(2)} -c_2\Delta t\frac{1}{\rho}\nabla \phi^{(2)}$を計算して第2段階の近似速度$\tilde{\boldsymbol{u}}^{(2)}$を得る．
7. ${\boldsymbol{u}^{*}}^{(3)} = \boldsymbol{u}^n+\Delta t \left[a_{3,1}F(\tilde{\boldsymbol{u}}^{(0)})+a_{3,2}F(\tilde{\boldsymbol{u}}^{(1)})+a_{3,3}F(\tilde{\boldsymbol{u}}^{(2)})\right]$を計算する．
8. $\nabla^2 p^{n+1}=\frac{\rho}{c_3\Delta t}\nabla\cdot{\boldsymbol{u}^*}^{(3)}$を解いて，圧力$p^{n+1}$を得る．
9. $\boldsymbol{u}^{n+1} = {\boldsymbol{u}^*}^{(3)} -c_3\Delta t\frac{1}{\rho}\nabla p^{n+1}$を計算して速度$\boldsymbol{u}^{n+1}$を得る．

係数$a_{s,i}, c_s$は，SKRK3と同じ値が用いられる．

後に，擬似的な圧力の外挿式はパラメータ$\alpha, \beta$を用いた次式に修正されている³．

\begin{align}
\phi^{(1)}=(1+\beta+c_1\alpha)\phi^{n-\beta}-(\beta+c_1\alpha)\phi^{n-1-\beta}\\
\phi^{(2)}=(1+\beta+c_2\alpha)\phi^{n-\beta}-(\beta+c_2\alpha)\phi^{n-1-\beta}
\end{align}

また，手順8.で計算されている$p^{n+1}$も$\phi^{n+1-\beta}$に置き換えられ，$p^{n+1}$は必要に応じて$\phi^{n-1-\beta}, \phi^{n-\beta}, \phi^{n+1-\beta}$を用いた補間によって計算される．

空間高次精度化

時間積分精度の記事ではあるが，空間離散化の高次精度化についても言及する必要がある．時間積分の実質的な（プログラムを作成し，それを実行することで得られる）精度を確認するには，ある$\Delta t_1$を用いて得られた誤差と，$\Delta t_2=\Delta t_1/10$を用いて得られた誤差の関係を見る．その際，空間離散化に起因する誤差が混在していると，時間積分に起因する誤差が，理論通り減少しなくなってしまう．
その例は，著者の過去の記事からも確認できる．リンク先の図では，クーラン数が0.1以下では4次精度の時間積分法を使っても誤差が減少しない．これはプログラムミスではなく，時間積分に起因する誤差が空間離散化に起因する誤差よりも小さくなったためである．

このような問題を回避するため，空間離散化に起因する誤差を十分小さくする必要がある．理論的には2次精度の中心差分を用いたとしても，離散点間隔$\Delta x$を十分に小さく取れば誤差は小さくなる．しかし，それには非常に多くの格子点が必要となり，計算負荷が高くなる．そこで，高次の空間離散化手法を用いることにして，格子点をさほど増やさずに誤差を十分に小さくすることを狙う．

支配方程式はスタガード格子上で離散化する．式(2)の移流項は，式(1)を考慮して発散型$\nabla\cdot(\boldsymbol{u}\otimes\boldsymbol{u})$に変形する．空間離散化には，2次，6次，8次，10次精度の中心差分を用い，移流項の離散化には，森西が示した適切な離散化形式⁴を用いる．

これまでベクトル表記していた量を，成分で表示する．すなわち，$\boldsymbol{u}=(u_1, u_2, u_3)$と表記し，軸方向も同様に$(x_1, x_2, x_3)$と書く．

離散演算子として，差分演算子と補間演算子を次のように定義する⁵．

\left. \frac{\delta_n f}{\delta_n x_1}\right|_{x_1, x_2, x_3}=\frac{f(x_1+n\Delta x_1/2, x_2, x_3)-f(x_1-n\Delta x_1/2,x_2,x_3)}{n\Delta x_1}

\left. \overline{f}^{nx_1}\right|_{x_1, x_2, x_3}=\frac{f(x_1+n\Delta x_1/2, x_2, x_3)+f(x_1-n\Delta x_1/2,x_2,x_3)}{2}

差分演算子は添字$1,2,3$について総和規約に従うが，補間演算子は総和規約に従わず，同項にある総和規約に従う添字に伴って変化する．

上式は2次精度の差分と補間である．高次精度の差分と補間は，異なる$n$を用いた2次精度の差分と補間を組み合わせて表現する．

\begin{align}
\frac{\partial^\text{6th} f}{\partial x_1}&
=\frac{150}{128}\frac{\delta_1 f}{\delta_1 x_1}
-\frac{25}{128}\frac{\delta_3 f}{\delta_3 x_1}
+\frac{3}{128}\frac{\delta_5 f}{\delta_5 x_1}\\
\frac{\partial^\text{8th} f}{\partial x_1}&
=\frac{1225}{1024}\frac{\delta_1 f}{\delta_1 x_1}
-\frac{245}{1024}\frac{\delta_3 f}{\delta_3 x_1}
+\frac{49}{1024}\frac{\delta_5 f}{\delta_5 x_1}
-\frac{5}{1024}\frac{\delta_7 f}{\delta_7 x_1}\\
\frac{\partial^\text{10th} f}{\partial x_1}&
=\frac{39690}{32768}\frac{\delta_1 f}{\delta_1 x_1}
-\frac{8820}{32768}\frac{\delta_3 f}{\delta_3 x_1}
+\frac{2268}{32768}\frac{\delta_5 f}{\delta_5 x_1}
-\frac{405}{32768}\frac{\delta_7 f}{\delta_7 x_1}
+\frac{35}{32768}\frac{\delta_9 f}{\delta_9 x_1}\\
\overset{\text{6th}}{\overline{f}^{x_1}}&
=\frac{150}{128}\overline{f}^{1x_1}
-\frac{25}{128}\overline{f}^{3x_1}
+\frac{3}{128}\overline{f}^{5x_1}\\
\overset{\text{8th}}{\overline{f}^{x_1}}&
=\frac{1225}{1024}\overline{f}^{1x_1}
-\frac{245}{1024}\overline{f}^{3x_1}
+\frac{49}{1024}\overline{f}^{5x_1}
-\frac{5}{1024}\overline{f}^{7x_1}\\
\overset{\text{10th}}{\overline{f}^{x_1}}&
=\frac{39690}{32768}\overline{f}^{1x_1}
-\frac{8820}{32768}\overline{f}^{3x_1}
+\frac{2268}{32768}\overline{f}^{5x_1}
-\frac{405}{32768}\overline{f}^{7x_1}
+\frac{35}{32768}\overline{f}^{9x_1}
\end{align}

移流項は，2次精度の補間と差分を用いると，次のように離散化することが適切であることが示されている⁶．

\frac{\partial u_j u_i}{\partial x_j} = \frac{\delta_1}{\delta_1 x_j}\left(\overline{u_j}^{1x_i}\overline{u_i}^{1x_j}\right)

森西は，移流項の高次精度化において，補間と差分をそれぞれ独立に高次精度化するのではなく，移流速度を高次精度補間して求め，差分の各項でその移流速度を使うことで適切に高次精度化できることを示している⁴⁵．
つまり，移流速度を

\begin{align}
\overset{\text{6th}}{\overline{U_j}^{x_i}}&:=
\frac{150}{128}\overline{u_j}^{1x_i}
-\frac{25}{128}\overline{u_j}^{3x_i}
+\frac{3}{128}\overline{u_j}^{5x_i}\\
\overset{\text{8th}}{\overline{U_j}^{x_i}}&:=
\frac{1225}{1024}\overline{u_j}^{1x_i}
-\frac{245}{1024}\overline{u_j}^{3x_i}
+\frac{49}{1024}\overline{u_j}^{5x_i}
-\frac{5}{1024}\overline{u_j}^{7x_i}\\
\overset{\text{10th}}{\overline{U_j}^{x_i}}&:=
\frac{39690}{32768}\overline{u_j}^{1x_i}
-\frac{8820}{32768}\overline{u_j}^{3x_i}
+\frac{2268}{32768}\overline{u_j}^{5x_i}
-\frac{405}{32768}\overline{u_j}^{7x_i}
+\frac{35}{32768}\overline{u_j}^{9x_i}
\end{align}

とおき，移流速度を用いて移流項を下記のように離散化する．

\begin{align}
\frac{\partial u_j u_i}{\partial x_j}&\overset{\text{6th}}{\approx}
\frac{150}{128}\left(\frac{\delta_1 \overset{\text{6th}}{\overline{U_j}^{x_i}}\overline{u_i}^{1x_j}}{\delta_1 x_j}\right)
-\frac{25}{128}\left(\frac{\delta_3 \overset{\text{6th}}{\overline{U_j}^{x_i}}\overline{u_i}^{3x_j}}{\delta_3 x_j}\right)
+\frac{3}{128}\left(\frac{\delta_5 \overset{\text{6th}}{\overline{U_j}^{x_i}}\overline{u_i}^{5x_j}}{\delta_5 x_j}\right)\\
\frac{\partial u_j u_i}{\partial x_j}&\overset{\text{8th}}{\approx}
\frac{1225}{1024}\left(\frac{\delta_1 \overset{\text{8th}}{\overline{U_j}^{x_i}}\overline{u_i}^{1x_j}}{\delta_1 x_j}\right)
-\frac{245}{1024}\left(\frac{\delta_3 \overset{\text{8th}}{\overline{U_j}^{x_i}}\overline{u_i}^{3x_j}}{\delta_3 x_j}\right)
+\frac{49}{1024}\left(\frac{\delta_5 \overset{\text{8th}}{\overline{U_j}^{x_i}}\overline{u_i}^{5x_j}}{\delta_5 x_j}\right)
-\frac{5}{1024}\left(\frac{\delta_7 \overset{\text{8th}}{\overline{U_j}^{x_i}}\overline{u_i}^{7x_j}}{\delta_7 x_j}\right)\\
\frac{\partial u_j u_i}{\partial x_j}&\overset{\text{10th}}{\approx}
\frac{39690}{32768}\left(\frac{\delta_1 \overset{\text{10th}}{\overline{U_j}^{x_i}}\overline{u_i}^{1x_j}}{\delta_1 x_j}\right)
-\frac{8820}{32768}\left(\frac{\delta_3 \overset{\text{10th}}{\overline{U_j}^{x_i}}\overline{u_i}^{3x_j}}{\delta_3 x_j}\right)
+\frac{2268}{32768}\left(\frac{\delta_5 \overset{\text{10th}}{\overline{U_j}^{x_i}}\overline{u_i}^{5x_j}}{\delta_5 x_j}\right)
-\frac{405}{32768}\left(\frac{\delta_7 \overset{\text{10th}}{\overline{U_j}^{x_i}}\overline{u_i}^{7x_j}}{\delta_7 x_j}\right)
+\frac{35}{32768}\left(\frac{\delta_9 \overset{\text{10th}}{\overline{U_j}^{x_i}}\overline{u_i}^{9x_j}}{\delta_9 x_j}\right)
\end{align}

拡散項$\frac{\partial^2 f}{\partial x_j \partial x_j}$は，1階微分の1階微分$\frac{\partial}{\partial x_j}\frac{\partial f}{\partial x_j}$として表現し，それぞれの微分を高次精度差分で計算する．

空間離散化についても，実質的な精度を確認する．$x\in[0, 2\pi]$の領域において$\frac{\partial uu}{\partial x}$を計算し，解析的に計算できる値との平均二乗平方誤差を計算する．$u=\sin(x)$とし，周期境界条件を課した．

下図は，格子点数を16, 32, 64, 128, 256, 512, 1024として計算したときの誤差の変化である．比較のために，傾き2, 6, 8, 10の線をそれぞれ点線，破線，1点鎖線，2点鎖線で併記している．格子点数を増加させると誤差は精度次数どおりに減少しており，作成したプログラムの正しさが示されている．

空間を2次精度で離散化すると，空間離散化に起因する誤差の影響を無視できるほど小さくすることは，現実的な格子点数ではできそうにない．波数1の正弦波に対して，6次精度では1024点，8次精度では256点，10次精度では128点使って解像すれば，空間離散化に起因する誤差を無視できるほど小さくできる．

時間積分の実質的な精度の調査

テスト問題

時間積分の実質的な精度を調査するために，Taylor-Green渦を計算対象とする．Taylor-Green渦はNavier-Stokes方程式の厳密解の一つであり，回転方向の異なる渦が$x_1,x_2$方向に周期的に存在する．流れ場が$x_1,x_2$方向にそれぞれ長さ$2\pi$で周期性を有している場合を考えると，Taylor-Green渦を表す速度はそれぞれ次式で表される．

\begin{align}
u_1 &= e^{-2\nu t}\sin{x_1} \cos{x_2}\\
u_2 &= -e^{-2\nu t}\cos{x_1}\sin{x_2}
\end{align}

$x_1 \times x_2 \in [0, 2\pi]\times [0, 2\pi]$を$128\times 128$の格子に分割し，$t=0\text{ s}$ から$1\text{ s}$まで計算した．計算時間間隔$\Delta t$には，$1, 0.1, 0.01, 0.001$の4通りを用い，$t=1\text{ s}$における$x_1$方向速度$u_1$の平均二乗平方誤差を調べた．動粘度$\nu$は$1/200$とした．$\Delta t=1$では安定条件は満足されないが，この条件では時間積分回数は1回だけであるため，直ちに発散するとは考えられない．

空間離散化には10次精度の補間と差分を用いるが，比較のために2次精度の補間と差分も実装する．圧力Poisson方程式は，FFTを利用して解いた．FFTの計算には，fortran-langコミュニティによって倍精度実数に拡張，パッケージ化されたFFTPACKを利用した．

時間積分には1次の前進Euler法，SKRK3，FastRK3（$\alpha=\beta=\frac{1}{2}$）を用いた．FastRK3では，過去2ステップの擬似的な圧力$\phi^{n-\beta}, \phi^{n-1-\beta}$が必要になるため，最初の2ステップはSKRK3を用いた．

実行結果

下図は$t=1$における$u_1$の平均二乗平方誤差である．空間離散化には2次精度の補間と差分を用いている．全ての時間積分法において，計算時間間隔$\Delta t$の現象と共に誤差も減少している．しかし，その減少の傾きは，図中に併記した傾き1の直線よりも緩やかであり，1次精度すら出ていない．空間離散化に起因する誤差が大きな割合を占めており，$\Delta t$を小さくしても平均二乗平方誤差が減らないためである．

空間離散化に10次精度の補間と差分を用いた場合の結果を下図に示す．空間離散化に起因する誤差が無視できる程度に小さくなったことから，時間積分に起因する誤差の変化を確認できる．

✕で示されたEuler法は，傾き1で誤差が減少しており，1次精度であることが確認できる．▲で示したSKRK3，○で示したFastRK3共に3次精度を達成できている．

試みに，FastRK3で要求される擬似的な圧力の初期値$\phi^{1-\beta}, \phi^{0-\beta}$を0として，全てFastRK3で計算した結果が●である．$\Delta t=1$の結果は1次精度のEuler法よりも4桁ほど誤差が大きく，初期値を0とした影響が非常に大きいことがわかる．また，$\Delta t$を1より小さくしていっても3次精度は達成できず，実質的な精度は2次精度に落ちる．これは1ステップと2ステップ目の計算結果の不正確さが影響していると考えられる．

まとめ

非圧縮性粘性流れの数値計算プログラムにおける時間積分精度を確認するための一つのやり方を示した．

空間離散化に起因する誤差が大きいと，時間積分の実質的な精度を確認できなくなる．そのため，空間離散化には極力高次精度の方法を用い，予め誤差が無視できるほど小さくなる格子点数を決定する．

高次精度の空間離散化手法を用いると，一般的には計算に用いる点が多くなる．壁面などの境界があると，境界近くで空間離散化精度を落とさざるを得なくなる．そういった問題を回避するために，周期境界条件を課すことができ，また厳密解を記述できるTaylor-Green渦をテスト問題として用いた．

Euler法，Sanderse-Korenによる3次のRunge-Kutta法，Aithal-FerranteによるFastRK3を実装してテストしたところ，いずれの方法も形式的な精度を実質的に達成できていることが確認できた．他の時間積分法を用いても，同様に精度を確認できると考えている．

FastRK3は，1ステップあたりに圧力Poisson方程式を1回解くだけでよく，実質的に3次精度を達成できる．ただし，擬似的な圧力の初期値を適当に与えると時間積分精度が落ちることがわかった．プログラムの実装は多少煩雑にはなるが，1,2ステップ目には自動出発可能な3次精度の時間積分法を用いる必要がある．また，3次精度を達成できるとはいっても，誤差の値自体はSanderse-Korenによる3次のRunge-Kutta法よりも7桁ほど大きくなる．

1. Sanderse, B. and Koren, B., Accuracy analysis of explicit Runge–Kutta methods applied to the incompressible Navier–Stokes equations, Journal of Computational Physics, Vol.231(8), (2012), pp.3041-3063. ↩ ↩²

2. Aithal, A.B. and Ferrante, A., A fast pressure-correction method for incompressible flows over curved walls, Journal of Computational Physics, Vol.421, (2020), 1096933. ↩ ↩²

3. Aithal, A.B., Tipirneni, M. and Ferrante, A., Temporal accuracy of FastRK3, Journal of Computational Physics, Vol.475, (2023), 111853. ↩ ↩²

4. 森西, 非圧縮性流体解析における差分スキームの保存特性（第2報, スタガードおよびコロケート格子系の差分スキーム）, 機論B, 62(604), (1996), pp.4098-4105. ↩ ↩²

5. 森西, 非圧縮性流体解析における差分スキームの保存特性（第1報, 解析的要求事項, 離散オペレータの定義, レギュラ格子系の差分スキーム）, 機論B, 62(604), (1996), pp.4090–4097. ↩ ↩²

6. 梶島, 対流項の差分形式とその保存性, 機論B, 60(574), (1994), pp. 2058-2063. ↩