MUSCLスキーム

偏微分方程式の研究において、MUSCL スキームは有限体積法であり、解がショック、不連続性、または大きな勾配を示す場合でも、特定のシステムに対して高精度の数値解を提供できます。MUSCL はMonotonic Upstream-centered Scheme for Conservation Laws (van Leer, 1979) の略で、この用語はBram van Leer (van Leer, 1979)の独創的な論文で導入されました。この論文で、彼は最初の高次全変動減少(TVD) スキームを構築し、2 次空間精度を実現しました。

このアイデアは、ゴドゥノフのスキームの区分定数近似を、前のタイムステップから取得したセル平均状態から導出された再構成状態に置き換えるというものです。各セルについて、傾斜制限された再構成された左および右状態が取得され、セル境界 (エッジ) でのフラックスを計算するために使用されます。これらのフラックスは、次に、リーマン ソルバーへの入力として使用でき、その後、ソリューションが平均化され、ソリューションを時間的に進めるために使用されます。あるいは、フラックスは、基本的にルサノフのようなスキームである リーマン ソルバーを使用しないスキームで使用することもできます。

MUSCL方式の基礎を、次の単純な1次スカラー1Dシステムを検討することで考察します。このシステムでは、正の方向に波が伝播すると仮定します。

${\displaystyle u_{t}+F_{x}\left(u\right)=0.\,}$

ここで、 は状態変数を表し、 はフラックス変数を表します。 ${\displaystyle u}$${\displaystyle F}$

ゴドゥノフの基本スキームは、各セルに対して区分定数近似を使用し、セル中心を としてインデックス付けした上記の問題の1次風上離散化を実現します。半離散スキームは次のように定義できます。 ${\displaystyle i}$

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F\left(u_{i}\right)-F\left(u_{i-1}\right)\right]=0.}$

この基本スキームでは、衝撃や急激な不連続性は不鮮明になりやすいため、これらを処理できません。この効果の例は、右に伝播するステップ波を伴う 1D 移流方程式を示す、右向きの図に示されています。シミュレーションは 200 セルのメッシュで実行され、4 次ルンゲ・クッタ時間積分器 (RK4) が使用されました。

不連続性のより高い解像度を提供するために、ゴドゥノフのスキームは各セルの区分線形近似を使用するように拡張することができ、その結果、空間的に2次の精度を持つ中心差分スキームが得られる。区分線形近似は、

${\displaystyle u\left(x\right)=u_{i}+{\frac {\left(x-x_{i}\right)}{\left(x_{i+1}-x_{i}\right)}}\left(u_{i+1}-u_{i}\right)\qquad \forall x\in (x_{i},x_{i+1}].}$

したがって、セル端でのフラックスを評価すると、次の半離散スキームが得られます。

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F\left(u_{i+1/2}\right)-F\left(u_{i-1/2}\right)\right]=0,}$

ここで、およびはセルエッジ変数の区分近似値、すなわち、 ${\displaystyle u_{i+1/2}}$${\displaystyle u_{i-1/2}}$

${\displaystyle u_{i+1/2}=0.5\left(u_{i}+u_{i+1}\right),}$

${\displaystyle u_{i-1/2}=0.5\left(u_{i-1}+u_{i}\right).}$

上記の 2 次スキームは滑らかな解に対してより高い精度を提供しますが、全変動減少(TVD) スキームではないため、不連続性や衝撃が存在する場合に解に誤った振動が生じます。この影響の例は、右に伝播するステップ波を伴う 1D 移流方程式を示す、反対側の図に示されています。この精度の低下は、ゴドゥノフの定理により予測されるものです。シミュレーションは 200 セルのメッシュで実行され、時間積分には RK4 が使用されました。 ${\displaystyle \,u_{t}+u_{x}=0}$

MUSCLベースの数値スキームは、傾斜制限された左と右の外挿状態を使用して、各セルに線形区分近似を使用するという考え方を拡張します。これにより、次の高解像度のTVD離散化スキームが得られます。

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F\left(u_{i+1/2}^{*}\right)-F\left(u_{i-1/2}^{*}\right)\right]=0.}$

これをもっと簡潔に書くと、

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}+{\frac {1}{\Delta x_{i}}}\left[F_{i+1/2}^{*}-F_{i-1/2}^{*}\right]=0.}$

数値フラックスは、連続フラックス関数の 1 次近似と 2 次近似の非線形組み合わせに対応します。 ${\displaystyle F_{i\pm 1/2}^{*}}$

記号およびは、（限定された外挿されたセルエッジ変数の）スキーム依存関数を表す。すなわち、 ${\displaystyle u_{i+1/2}^{*}}$${\displaystyle u_{i-1/2}^{*}}$

${\displaystyle u_{i+1/2}^{*}=u_{i+1/2}^{*}\left(u_{i+1/2}^{L},u_{i+1/2}^{R}\right),u_{i-1/2}^{*}=u_{i-1/2}^{*}\left(u_{i-1/2}^{L},u_{i-1/2}^{R}\right),}$

ここで、風下勾配を使用すると、

${\displaystyle u_{i+1/2}^{L}=u_{i}+0.5\phi \left(r_{i}\right)\left(u_{i+1}-u_{i}\right),u_{i+1/2}^{R}=u_{i+1}-0.5\phi \left(r_{i+1}\right)\left(u_{i+2}-u_{i+1}\right),}$

${\displaystyle u_{i-1/2}^{L}=u_{i-1}+0.5\phi \left(r_{i-1}\right)\left(u_{i}-u_{i-1}\right),u_{i-1/2}^{R}=u_{i}-0.5\phi \left(r_{i}\right)\left(u_{i+1}-u_{i}\right),}$

そして

${\displaystyle r_{i}={\frac {u_{i}-u_{i-1}}{u_{i+1}-u_{i}}}.}$

関数は、区分近似の傾きを制限するリミッタ関数であり、解が TVD であることを保証し、それによって不連続点や衝撃の周囲で発生する可能性のある誤った振動を回避します (フラックス リミッタセクションを参照)。リミッタは のときはゼロに等しく、 のときは 1 に等しくなります。したがって、TVD 離散化の精度は、局所的極値では 1 次まで低下しますが、領域の滑らかな部分では 2 次まで低下する傾向があります。 ${\displaystyle \phi \left(r_{i}\right)}$${\displaystyle r\leq 0}$${\displaystyle r=1}$

このアルゴリズムは簡単に実装できます。Kurganovと Tadmor のスキーム(下記参照)などの適切なスキームを選択すると、標準的な数値積分手法を使用してソリューションを進めることができます。 ${\displaystyle F_{i+1/2}^{*}}$

Kurganov と Tadmor (KT) の中心スキーム(Kurganov と Tadmor、2000)の前身は、 Nessyahu と Tadmor (NT) のスタガード中心スキーム(Nessyahu と Tadmor、1990) です。これは、MUSCL 再構成を使用する、リーマン ソルバーを必要としない 2 次高解像度スキームです。これは完全に離散的な方法で、実装が簡単で、スカラー問題とベクトル問題に使用でき、高次再構成が追加された Rusanov フラックス (ローカル Lax-Friedrichs フラックスとも呼ばれる) として考えることができます。このアルゴリズムは中心差分に基づいており、高勾配現象を示すシステムを記述する PDE の解を取得するために使用される場合、リーマン型ソルバーと同等のパフォーマンスを発揮します。

KT 方式は NT 方式を拡張したもので、元の NT 方式よりも数値粘性が低くなっています。また、完全離散方式または半離散方式として実装できるという利点もあります。ここでは半離散方式について検討します。

計算は以下の通りです。

${\displaystyle F_{i-{\frac {1}{2}}}^{*}={\frac {1}{2}}\left\{\left[F\left(u_{i-{\frac {1}{2}}}^{R}\right)+F\left(u_{i-{\frac {1}{2}}}^{L}\right)\right]-a_{i-{\frac {1}{2}}}\left[u_{i-{\frac {1}{2}}}^{R}-u_{i-{\frac {1}{2}}}^{L}\right]\right\}.}$

${\displaystyle F_{i+{\frac {1}{2}}}^{*}={\frac {1}{2}}\left\{\left[F\left(u_{i+{\frac {1}{2}}}^{R}\right)+F\left(u_{i+{\frac {1}{2}}}^{L}\right)\right]-a_{i+{\frac {1}{2}}}\left[u_{i+{\frac {1}{2}}}^{R}-u_{i+{\frac {1}{2}}}^{L}\right]\right\}.}$

ここで、局所伝播速度は、セル上のヤコビアン固有値の絶対値の最大値であり、 ${\displaystyle a_{i\pm {\frac {1}{2}}}\ }$${\displaystyle F\left(u\left(x,t\right)\right)}$${\displaystyle {i},{i\pm 1}}$

${\displaystyle a_{i+{\frac {1}{2}}}\left(t\right)=\max \left[\rho \left({\frac {\partial F\left(u_{i+1/2}^{L}\left(t\right)\right)}{\partial u}}\right),\rho \left({\frac {\partial F\left(u_{i+1/2}^{R}\left(t\right)\right)}{\partial u}}\right),\right]}$

スペクトル半径を表す。 ${\displaystyle \rho \left({\frac {\partial F\left(u\left(t\right)\right)}{\partial u}}\right)\ }$${\displaystyle {\frac {\partial F\left(u\left(t\right)\right)}{\partial u}}.}$

これらのCFL関連の速度を超える特性情報は必要ありません。

上記のフラックス計算は、 Lax-Friedrichs フラックスと呼ばれることがよくあります(ただし、このようなフラックス表現は Lax (1954) ではなく Rusanov (1961) に登場していることは注目に値します)。

高解像度スキームの使用の有効性の例が、右に伝播するステップ波を伴う 1D 移流方程式 を示す反対の図に示されています。シミュレーションは、Superbee リミッターを備えた Kurganov および Tadmor 中心スキームを使用して 200 セルのメッシュで実行され、時間積分には RK-4 が使用されました。このシミュレーション結果は、上記の 1 次風上差分および 2 次中心差分の結果と非常によく対照的です。このスキームは、方程式のセットに適用した場合にも良好な結果をもたらします。このスキームをオイラー方程式に適用した結果については、以下を参照してください。ただし、たとえば Superbee リミッターは、一部の滑らかな波に対して非現実的なシャープニングを引き起こす可能性があるため、適切なリミッターを選択する際には注意が必要です。 ${\displaystyle u_{t}+u_{x}=0\ }$

このスキームは、拡散項が存在する場合には、それを容易に含めることができる。例えば、上記の1Dスカラー問題を拡散項を含むように拡張すると、次のようになる。

${\displaystyle u_{t}+F_{x}\left(u\right)=Q_{x}\left(u,u_{x}\right),}$

これに対してクルガノフとタドモアは次のような中心差分近似を提案している。

${\displaystyle {\frac {\mathrm {d} u_{i}}{\mathrm {d} t}}=-{\frac {1}{\Delta x_{i}}}\left[F_{i+{\frac {1}{2}}}^{*}-F_{i-{\frac {1}{2}}}^{*}\right]+{\frac {1}{\Delta x_{i}}}\left[P_{i+{\frac {1}{2}}}-P_{i-{\frac {1}{2}}}\right].}$

どこ、

${\displaystyle P_{i+{\frac {1}{2}}}={\frac {1}{2}}\left[Q\left(u_{i},{\frac {u_{i+1}-u_{i}}{\Delta x_{i}}}\right)+Q\left(u_{i+1},{\frac {u_{i+1}-u_{i}}{\Delta x_{i}}}\right)\right],}$

${\displaystyle P_{i-{\frac {1}{2}}}={\frac {1}{2}}\left[Q\left(u_{i-1},{\frac {u_{i}-u_{i-1}}{\Delta x_{i-1}}}\right)+Q\left(u_{i},{\frac {u_{i}-u_{i-1}}{\Delta x_{i-1}}}\right).\right]}$

アルゴリズムの詳細 (完全バージョンと半離散バージョン) とその導出については、オリジナルの論文 (Kurganov と Tadmor、2000) に記載されています。また、1D と 2D の例も多数掲載されています。Nessyahu と Tadmor (1990) による以前の関連論文にも追加情報が記載されています。

注:このスキームは、もともとは線形外挿 に基づく 2 次スキームとして Kurganov と Tadmor によって発表されました。その後の論文 (Kurganov と Levy、2000) では、このスキームが 3 次スキームの基礎にもなり得ることが示されています。放物線再構成 (3 次) を使用したこのスキームの 1D 移流の例とオイラー方程式の例は、以下の放物線再構成とオイラー方程式のセクションに示されています。

線形外挿の考え方を高次の再構成に拡張することが可能であり、その例が反対側の図に示されています。ただし、この場合、左と右の状態は、2 次の風上バイアス差分方程式の補間によって推定されます。これにより、空間的に 3 次の精度を持つ放物線再構成スキームが得られます。

我々はケルマニ（Kermani, et al., 2003）のアプローチに従い、3次の風上バイアススキームを提示する。ここでもシンボルとがスキーム依存関数（限定再構成セルエッジ変数）を表す。しかし、この場合、それらは放物線再構成された状態、すなわち、 ${\displaystyle u_{i+{\frac {1}{2}}}^{*}}$${\displaystyle u_{i-{\frac {1}{2}}}^{*}}$

${\displaystyle u_{i+{\frac {1}{2}}}^{*}=f\left(u_{i+{\frac {1}{2}}}^{L},u_{i+{\frac {1}{2}}}^{R}\right),\quad u_{i-{\frac {1}{2}}}^{*}=f\left(u_{i-{\frac {1}{2}}}^{L},u_{i-{\frac {1}{2}}}^{R}\right),}$

そして

${\displaystyle u_{i+{\frac {1}{2}}}^{L}=u_{i}+{\frac {\phi \left(r_{i}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i-{\frac {1}{2}}}+\left(1+\kappa \right)\delta u_{i+{\frac {1}{2}}}\right],}$

${\displaystyle u_{i+{\frac {1}{2}}}^{R}=u_{i+1}-{\frac {\phi \left(r_{i+1}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i+{\frac {3}{2}}}+\left(1+\kappa \right)\delta u_{i+{\frac {1}{2}}}\right],}$

${\displaystyle u_{i-{\frac {1}{2}}}^{L}=u_{i-1}+{\frac {\phi \left(r_{i-1}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i-{\frac {3}{2}}}+\left(1+\kappa \right)\delta u_{i-{\frac {1}{2}}}\right],}$

${\displaystyle u_{i-{\frac {1}{2}}}^{R}=u_{i}-{\frac {\phi \left(r_{i}\right)}{4}}\left[\left(1-\kappa \right)\delta u_{i+{\frac {1}{2}}}+\left(1+\kappa \right)\delta u_{i-{\frac {1}{2}}}\right].}$

ここで= 1/3 であり、 ${\displaystyle \kappa \ }$

${\displaystyle \delta u_{i+{\frac {1}{2}}}=\left(u_{i+1}-u_{i}\right),\quad \delta u_{i-{\frac {1}{2}}}=\left(u_{i}-u_{i-1}\right),}$

${\displaystyle \delta u_{i+{\frac {3}{2}}}=\left(u_{i+2}-u_{i+1}\right),\quad \delta u_{i-{\frac {3}{2}}}=\left(u_{i-1}-u_{i-2}\right),}$

リミッター関数は上記と同じです。 ${\displaystyle \phi \left(r\right)\ }$

放物線再構成は実装が簡単で、上記の線形外挿の代わりに Kurganov および Tadmor 方式で使用できます。これは、KT 方式の空間解を 3 次まで上げる効果があります。これは、オイラー方程式を解くときにうまく機能します (以下を参照)。この空間次数の増加は、滑らかな解の場合、2 次方式よりも一定の利点がありますが、衝撃の場合はより散逸的です。反対側の図を、線形外挿と Superbee リミッタを使用した KT アルゴリズムを使用して得られた上記の解と比較してください。このシミュレーションは、同じ KT アルゴリズムを使用して放物線再構成で 200 セルのメッシュで実行されました。時間積分は RK-4 によって行われ、スプリアス振動を回避するために、van Albada リミッタの代替形式 が使用されました。 ${\displaystyle \phi _{va}(r)={\frac {2r}{1+r^{2}}}\ }$

簡単のため、熱伝達と体積力のない1次元の場合を考える。したがって、保存ベクトル形式では、一般的なオイラー方程式は次のように簡略化される。

${\displaystyle {\frac {\partial \mathbf {U} }{\partial t}}+{\frac {\partial \mathbf {F} }{\partial x}}=0,}$

どこ

${\displaystyle \mathbf {U} ={\begin{pmatrix}\rho \\\rho u\\E\end{pmatrix}}\qquad \mathbf {F} ={\begin{pmatrix}\rho u\\p+\rho u^{2}\\u(E+p)\end{pmatrix}},\qquad }$

ここで は状態のベクトルであり、 はフラックスのベクトルです。 ${\displaystyle {\mbox{U}}}$${\displaystyle {\mbox{F}}}$

上記の式は、質量、運動量、エネルギーの保存を表しています。したがって、3 つの式と 4 つの未知数(密度) (流体速度) (圧力) (総エネルギー)があります。総エネルギーは次のように表されます。 ${\displaystyle \rho }$${\displaystyle u}$${\displaystyle p}$${\displaystyle E}$

${\displaystyle E=\rho e+{\frac {1}{2}}\rho u^{2},}$

ここで、特定の内部エネルギーを表します。 ${\displaystyle e\ }$

システムを閉じるには状態方程式が必要です。私たちの目的に適したものは

${\displaystyle p=\rho \left(\gamma -1\right)e,}$

ここで、は流体の 比熱の比に等しくなります。${\displaystyle \gamma \ }$${\displaystyle \left[c_{p}/c_{v}\right]}$

上の簡単な1Dの例で示したように、各状態変数の左と右の外挿状態を取得することで、次に進めます。密度については、次の式が得られます。

${\displaystyle \rho _{i+{\frac {1}{2}}}^{*}=\rho _{i+{\frac {1}{2}}}^{*}\left(\rho _{i+{\frac {1}{2}}}^{L},\rho _{i+{\frac {1}{2}}}^{R}\right),\quad \rho _{i-{\frac {1}{2}}}^{*}=\rho _{i-{\frac {1}{2}}}^{*}\left(\rho _{i-{\frac {1}{2}}}^{L},\rho _{i-{\frac {1}{2}}}^{R}\right),}$

どこ

${\displaystyle \rho _{i+{\frac {1}{2}}}^{L}=\rho _{i}+0.5\phi \left(r_{i}\right)\left(\rho _{i}-\rho _{i-1}\right),\quad \rho _{i+{\frac {1}{2}}}^{R}=\rho _{i+1}-0.5\phi \left(r_{i+1}\right)\left(\rho _{i+1}-\rho _{i}\right),}$

${\displaystyle \rho _{i-{\frac {1}{2}}}^{L}=\rho _{i-1}+0.5\phi \left(r_{i-1}\right)\left(\rho _{i}-\rho _{i-1}\right),\quad \rho _{i-{\frac {1}{2}}}^{R}=\rho _{i}-0.5\phi \left(r_{i}\right)\left(\rho _{i+1}-\rho _{i}\right).}$

同様に、運動量、全エネルギーについても、速度は運動量から計算され、圧力は状態方程式から計算されます。 ${\displaystyle \rho u}$${\displaystyle E}$${\displaystyle u}$${\displaystyle p}$

限定された外挿状態が得られたら、これらの値を使用してエッジフラックスを構築します。エッジフラックスがわかれば、半離散スキームを構築できます。

${\displaystyle {\frac {\mathrm {d} \mathbf {U} _{i}}{\mathrm {d} t}}=-{\frac {1}{\Delta x_{i}}}\left[\mathbf {F} _{i+{\frac {1}{2}}}^{*}-\mathbf {F} _{i-{\frac {1}{2}}}^{*}\right].}$

標準的な数値手法を使用して積分することで、解を求めることができます。

上記は、MUSCL 方式の基本的な考え方を示しています。ただし、オイラー方程式の実用的な解を得るには、関数を定義するために適切な方式 (上記の KT 方式など) も選択する必要があります。 ${\displaystyle \mathbf {F} _{i\pm {\frac {1}{2}}}^{*}}$

右の図は、上記の高解像度の Kurganov と Tadmor の中央スキーム (KT) と線形外挿および Ospre リミッタを使用した、GA Sod の衝撃波管問題 (Sod、1978) の 2 次解を示しています。これは、オイラー方程式を解くための MUSCL アプローチの有効性を明確に示しています。シミュレーションは、KT アルゴリズムとOspre リミッタを使用するように調整された Matlab コード (Wesseling、2001) を使用して、200 セルのメッシュで実行されました。時間積分は、4 次 SHK (RK-4 と同等の性能) 積分器によって実行されました。次の初期条件 ( SI単位) が使用されました。

• 左圧力 = 100000 [Pa];
• 圧力右= 10000 [Pa];
• 密度左 = 1.0 [kg/m3];
• 密度右 = 0.125 [kg/m3];
• 長さ = 20 [m];
• 左速度 = 0 [m/s];
• 速度右 = 0 [m/s];
• 持続時間 =0.01 [秒];
• ラムダ = 0.001069 (Δt/Δx)。

右の図は、上記の高解像度の Kurganov および Tadmor 中心スキーム (KT) を使用し、放物線再構成とファン アルバダ リミッタを使用した、GA Sod の衝撃波管問題 (Sod、1978) の 3 次解を示しています。これもまた、オイラー方程式を解くための MUSCL アプローチの有効性を示しています。シミュレーションは、放物線外挿とファン アルバダリミッタを使用した KT アルゴリズムを使用するように適合された Matlab コード (Wesseling、2001) を使用して、200 セルのメッシュで実行されました。疑似振動を回避するために、ファン アルバダ リミッタの代替形式が使用されました。時間積分は、4 次 SHK 積分器によって実行されました。同じ初期条件が使用されました。 ${\displaystyle \phi _{va}(r)={\frac {2r}{1+r^{2}}}\ }$

オイラー方程式を高い精度で解くための様々な高解像度スキームが開発されている。そのようなスキームの例としては、

• オッシャー計画、そして
• Liou -Steffen AUSM（移流上流分割法）スキーム。

これらの方法やその他の方法の詳細については、以下の参考文献を参照してください。Kurganov と Tadmor の中央スキームのオープン ソース実装については、以下の外部リンクを参照してください。

• 有限体積法
• フラックスリミッター
• ゴドゥノフの定理
• 高解像度スキーム
• ライン法
• セルゲイ・K・ゴドゥノフ
• 総変動は減少している
• ソッドショックチューブ

• Kermani, MJ, Gerber, AG, Stockie, JM (2003)、「Roeのスキームを使用した熱力学的水分予測」、イラン航空宇宙学会第4回会議、アミール・カビール工科大学、テヘラン、イラン、1月27～29日。[1]
• Kurganov, AlexanderとEitan Tadmor (2000)、非線形保存則と対流拡散方程式のための新しい高解像度中心スキーム、J. Comput. Phys.、160、241–282。[2]
• Kurganov, Alexander and Doron Levy (2000)、「保存則と対流拡散方程式のための3次半離散中心スキーム」、SIAM J. Sci. Comput.、22、1461–1488。[3]
• Lax, PD (1954). 非線形双曲方程式の弱解とその数値計算、Comm. Pure Appl. Math.、VII、pp159–193。
• Leveque, RJ (2002).双曲型問題に対する有限体積法、ケンブリッジ大学出版局。
• van Leer, B. (1979)「究極の保守的差分スキームに向けて、V. Godunov法の2次続編」、J. Com. Phys .、32、101–136。
• Nessyahu, H. およびE. Tadmor (1990)、双曲型保存則に対する非振動中心差分、J. Comput. Phys.、87、408–463。[4]
• Rusanov, VV (1961). 非定常衝撃波と障害物の交差の計算、J. Comput. Math. Phys. USSR、1、pp267–279。
• Sod, GA (1978)「収束する円筒状衝撃波の数値的研究」流体力学誌、83、785–794。
• Toro, EF (1999)、「Riemann Solvers and Numerical Methods for Fluid Dynamics」、Springer-Verlag。
• Wesseling, Pieter (2001)、「計算流体力学の原理」、Springer-Verlag。

• Hirsch, C. (1990)、「内部流れと外部流れの数値計算」、第 2 巻、Wiley。
• Laney, Culbert B. (1998)、「計算ガスダイナミクス」、ケンブリッジ大学出版局。
• Tannehill, John C.、他 (1997)、「計算流体力学と熱伝達」、第 2 版、Taylor and Francis。

• GEES – Fortranで書かれた、Kurganov と Tadmor の中心スキームを使用してオイラー方程式を解くオープンソース コード(作成者: Arno Mayrhofer)