【CFD/格子法】時間差分スキーム【格子法入門 #12】

子法入門第11回は、ナビエ・ストークス方程式の離散化に関する「時間差分スキーム」について説明します。

前回は、、偏微分方程式の分類分けである「双曲型・放物型・楕円型」についてお話しました。下記からどうぞ。

リンクーーーーーーーーーーーーーーー

重要なポイントは下記の通りです。

• ２階微分方程式は、双曲型・放物型・楕円型の３つに分けられる
• 移流方程式は双曲型にあたる
• 非定常拡散方程式は放物型にあたる
• 定常拡散方程式は楕円型にあたる
• 楕円型であるポアソン方程式には、収束計算が使われる

「双曲型・放物型・楕円型」はあくまで数学的な表現のため理解は必須ではないですが、覚えておくと書籍や論文が読みやすくなるでしょう。

今回は、結構重要だけで今まで説明できていなかった時間差分スキームについて説明したいと思います。

時間差分スキーム

流体解析で使われるナビエ・ストークス方程式は下記のように書かれます。

ここで左辺の項に注目すると、時間$t$の微分が含まれることがわかります。これは非定常項と呼ばれ、時間変化を意味する項です。

現実の時間は連続的に流れるので時間微分という概念が存在しますが、PCの中ではそうはいきません。

０と１の世界ですから、滑らかな変化をはっきりさせるために「離散化」のプロセスが必要になります。

そこで今回は、時間差分のスキームについて説明していきます。

非定常問題

ナビエ・ストークス方程式だとややこしすぎるので、時間微分だけに注目できる式を代わりに使いましょう。

ここからは、下記の非定常問題について扱います。

$$ \frac{\partial \phi}{\partial t} = f (\phi) $$

この式を離散化すると、$t=t_0$からはじまり$t=t_1,t_2 …$と進むと、任意の変数$\phi$も$\phi = \phi_0,\phi_1,\phi_2 …$と進んでいきます。

計算対象が決まったので、ここから時間離散化スキームについて見ていきましょう。

オイラー法（陽解法）

まずは最もシンプルな離散化について考えましょう。

現在の時刻を$n$とすると、未来の時刻$n+1$のときの値$\phi^{n+1}$を求めたいとします。

その場合、オイラー法を使うと非定常項を下記のように離散化できます。

$$ \frac{\partial \phi}{\partial t} = \frac{\phi^{n+1} – \phi^n}{ \Delta t}$$

これは最も一般的な離散化で、「$\Delta t$進む間に、１ステップ分$\phi$が変化する」ということを意味しています。

高校数学で考えると当たり前かもしれません。しかしこれは線形補間であるため、精度を上げるならもっと良い方法がたくさんあります。今後、もっと良い方法について説明していきます。

オイラー法を使って未来の時刻$n+1$の値を予測しようとすると、下記のように書けます。

$$ \phi^{n+1} = \phi^n + \Delta t f (\phi^n) $$

オイラー法は「陽解法」に分類されます。

陽解法とは、既知の値だけで計算できる式のことです。といっても、ここでは当たり前にしか見えないと思うので、ここからは陽解法と対になる「陰解法」について紹介します。

陰解法

陰解法は陽解法と対になる手法です。陰解法で使うのは「未知の値を含む式」となります。

具体的には、下記のような式を使います。

$$ \phi^{n+1} = \phi^n + \Delta t f (\phi^{n+1}) $$

ここで陰解法と異なるのは、右辺の$f (\phi^{n+1})$だけです。

陰解法では右辺に未知数が含まれることで、単純な代入操作では解けなくなります。

一般的には、１回の時間ステップに対してさらに細かいステップを加えた収束計算が行われます。

クランク・ニコルソン法

クランク・ニコルソン法は、陽解法と陰解法の中間に当たる手法です。

クランク・ニコルソン法では下記の式が用いられます。

$$ \phi^{n+1} = \phi^n + \frac{1}{2}\Delta t f (\phi^{n+1} + \phi^n) $$

右辺の関数$f$が$f (\phi^{n+1} + \phi^n)$となっており、それを$frac{1}{2}$しています。

つまり、陽的な$\phi^n$と陰的な$\phi^{n+1}$をブレンドした手法であることがわかります。

オイラー法は１次精度であるのに対し、クランク・ニコルソン法は２次精度であることが特徴です。

アダムス・バッシュフォース法

アダムス・バッシュフォース法は、陽的な時間スキームの一種です。

オイラー法では過去の１ステップの値しか使わなかったのに対し、アダムス・バッシュフォース法では2,3…ステップの値も使うことで高精度にします。

２次精度だと下記の式となります。

$$ \phi^{n+1} = \phi^n + \Delta t [ \frac{3}{2} f (\phi^n) – \frac{1}{2} f (\phi^{n-1}) ] $$

現在の値$f(\phi^n)$と前のステップの値$f(\phi^{n-1})$というように、前ステップの情報を追加することで予測の精度を高めています。

アダムス・バッシュフォース法は、時間方向にステンシルが長いため、１ステップで急激な変化が起こっても異常な値を取りにくい手法です。そのため、計算が安定化しやすいです。

ただし欠点として、$f(\phi^{n-1})$という変数を追加で保持する必要があり、メモリ使用量が急激に増加するため、大規模流体解析では避けられがちです。

ルンゲクッタ法

ルンゲクッタ法も陽的な手法です。

ルンゲクッタ法は、陽的に多段の計算、つまり１時間ステップに対して何ステップ家に分けた計算を行うことで精度を向上させます。

オイラー法は１段１次のルンゲクッタ法とも言えます。

下記に４段４次のルンゲクッタ法を示します。

$$ \phi^1 = \phi^n + \Delta f (\phi_n) $$

$$ \phi^2 = \phi^n + \Delta f (\frac{\phi_n + \phi^1}{2}) $$

$$ \phi^3 = \phi^n + \Delta f (\frac{\phi_n + \phi^2}{2}) $$

$$ \phi^4 = \phi^n + \Delta f (\phi_3) $$

$$ \phi^{n+1} = \frac{1}{6} ( \phi^1 + 2 \phi^2 + 2 \phi^3 + \phi^4 )$$

元の値$fai^n$に対して$\phi^1 , \phi^2 , \phi^3 , \phi^4 $と補正を入れながら進めていくことで、それらを統合した$\phi^{n+1}$を計算します。

ルンゲクッタ法は、メモリ使用量が少なく、手軽に精度が上げられるので、流体解析ではよく利用されています。

TVDルンゲクッタ法

TVDルンゲクッタ法は、TVD条件と呼ばれる安定条件に基づいたルンゲクッタ法です。

TVD条件を用いることで、非常に安定した解析を行うことができます。

まずはTVDについて説明していきましょう。

TVD(Total Valiation Diminishing：全変動減少)

TVDは単調性を維持するためにスキームをコントロールする手法です。

単調性とは、時間とともに変動が大きくならないという条件です。

TVDではTV（全変動）を利用します。TVは下記の式で表されます。

$$ TV(u^n) = \Sigma_i | \phi^n_j – \phi^n_i|$$

そして、下記の式を満たすとき、TV安定、つまりTVDであるということになります。

$$ TV(u\phi^{n+1}) \leq TV(\phi^n) $$

上記の式は、時間ステップが進行してもTVが大きくならないことを意味しています。

TVDルンゲクッタ

一般的なルンゲクッタ式は下記のように表せます。

$$ \phi^{(i)} = \phi^{(0)} + \Delta \sum^{i-1}_{k=0} c_{ik} f (\phi^{(k)}) $$

ここで、$i=1,..m$であり、$\phi^{(0)} = \phi^n$、$\phi^{(m)} = \phi^{n+1}$です。

これに対して、TVDがルンゲクッタ式は少し違います。

$$ \phi^{(i)} = \sum^{i-1}_{k=0} { \alpha_{ik} \phi^{k} + \beta_{ik} f (\phi^{(k)}) }$$

ここで、$\alpha_{ik} \geq 0$かつ$\sum^{i-1}_{k=0} alpha_{ik} = 1$であり、$\beta_{ik} = c_{ik} – \sum^{i-1}_{l=k+1} c_{lk} \alpha_{il}$です。これらの条件を元に$\alpha$と$\beta$を決定します。

今回は$\alpha$と$beta$に関する計算は省略します。

TVD３段ルンゲクッタ

ここでは、TVDにより安定させた３段のルンゲクッタを紹介します。

式は下記のようになります。

$$ \phi^{(1)} = \phi^{(0)} + \Delta t f (\phi^{(0)}) $$

$$ \phi^{(2)} = \frac{3}{4} \phi^{(0)} + \frac{1}{4} (\phi^{(1)} + \Delta t f (\phi^{(1)}) ) $$

$$ \phi^{(3)} = \frac{1}{3} \phi^{(0)} + \frac{2}{3} (\phi^{(2)} + \Delta t f (\phi^{(2)}) ) $$

１回の時間ステップにて上記計算を行うことで、精度が高くかつ安定した計算を行うことができます。

おわりに

今回は時間差分スキームについて説明しました。

格子法では主に移流項の空間差分スキームが注目されますが、時間差分スキームも当然精度に影響を与えます。

重要なポイントは下記の通りです。

• 一番簡単で一般的なのは、オイラー法（陽解法）である
• 陰解法を使うと、収束計算が必要になる
• クランクニコルソン法は、陽解法と陰解法の中間に当たる
• アダムス・バッシュフォース法は、前ステップの情報を増やして計算を安定させる
• ルンゲクッタ法は１タイムステップに対して多段に分けて計算する
• TVDルンゲクッタ法は、TV条件と呼ばれる安定条件に基づいたルンゲクッタ法である

様々な種類があり、どれも一長一短です。

CFD計算で高精度を出したいときによく利用されるのは、ルンゲクッタ法でしょう。

ルンゲクッタ法はCFD以外の他の分野でも重要な考え方なので、ぜひ基本だけでも理解しておきましょう。