Finite volume method in partial differential equations
偏微分方程式 の研究において 、 MUSCL法は 有限体積法 であり、 与えられたシステムに対して、解がショック、不連続性、または大きな勾配を示す場合でも、非常に正確な数値解を与えることができます。MUSCLは、 Monotonic Upstream-centered Scheme for Conservation Laws (van Leer, 1979)の略で、この用語は Bram van Leer (van Leer, 1979)による独創的な論文で導入されました 。この論文で、彼は初めて 高次 の全変動減少 (TVD)法を構築し、 2次の空間精度を実現しました。
そのアイデアは、ゴドゥノフ法の 区分定数近似を 、前の時間ステップで得られたセル平均状態から導出される再構成状態に置き換えることです。各セルについて、傾きが制限された再構成された左状態と右状態が得られ、セル境界（エッジ）におけるフラックスを計算するために使用されます。これらのフラックスは、 リーマンソルバー への入力として使用され、その後、解が平均化されて、解を時間的に進めるために使用されます。あるいは、フラックスは、基本的にルサノフ法のようなリーマンソルバーフリー 法で使用できます 。
線形再構成 1次元の移流方程式 。ステップ波が右方向に伝播する。解析解と、一次の風上空間離散化スキームに基づくシミュレーションを示す。 u t + u x = 0 {\displaystyle u_{t}+u_{x}=0} MUSCL法の基礎を、次のような単純な1次スカラー1次元システムで考察する。このシステムでは、正の方向に伝播する波動があると仮定する。
u t + F x ( u ) = 0. {\displaystyle u_{t}+F_{x}\left(u\right)=0.\,} ここで、 は 状態変数を表し、は フラックス 変数を表します 。 u {\displaystyle u} F {\displaystyle F}
ゴドゥノフの基本スキームは、各セルに対して区分定数近似を使用し、上記の問題を1次の風上離散化して、セル中心を とインデックス付けする 。半離散スキームは次のように定義できる。 i {\displaystyle i}
d u i d t + 1 Δ x i [ F ( u i ) − F ( u i − 1 ) ] = 0. {\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.} この基本スキームでは、衝撃波や鋭い不連続面はぼやけてしまう傾向があるため、これらを扱うことができません。この影響の一例が、右方向に伝播するステップ波を伴う1次元移流方程式を示す右図に示されています。シミュレーションは200セルのメッシュで実行され、4次 ルンゲ・クッタ 時間積分器（RK4）が使用されました。
不連続性のより高い分解能を実現するために、ゴドゥノフの手法を各セルの区分線形近似を用いるように拡張することができ、これにより空間的に 2次の精度を持つ 中心差分 法 が得られる。区分線形近似は以下から得られる。
u ( x ) = u i + ( x − x i ) ( x i + 1 − x i ) ( u i + 1 − u i ) ∀ x ∈ ( x i , x i + 1 ] . {\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}].} したがって、セル端でのフラックスを評価すると、次のような半離散的なスキームが得られます。
d u i d t + 1 Δ x i [ F ( u i + 1 / 2 ) − F ( u i − 1 / 2 ) ] = 0 , {\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,} 1次元の移流方程式 。ステップ波が右方向に伝播する。解析解と、2次中心差分空間離散化スキームに基づくシミュレーションを示す。 u t + u x = 0 {\displaystyle u_{t}+u_{x}=0} ここで 、 およびはセルエッジ変数の区分近似値、 すなわち 、 u i + 1 / 2 {\displaystyle u_{i+1/2}} u i − 1 / 2 {\displaystyle u_{i-1/2}}
u i + 1 / 2 = 0.5 ( u i + u i + 1 ) , {\displaystyle u_{i+1/2}=0.5\left(u_{i}+u_{i+1}\right),} u i − 1 / 2 = 0.5 ( u i − 1 + u i ) . {\displaystyle u_{i-1/2}=0.5\left(u_{i-1}+u_{i}\right).} 上記の2次スキームは滑らかな解に対してより高い精度を提供しますが、 全変動減衰 （TVD）スキームではないため、不連続性や衝撃が存在する場合、解に不要な振動が生じます。この影響の例は、右方向に伝播するステップ波を伴う1次元移流方程式を示す右図に示されています。この精度の低下は、 ゴドゥノフの定理 により予測されるものです 。シミュレーションは200セルのメッシュで実行され、時間積分にはRK4を使用しました。 u t + u x = 0 {\displaystyle \,u_{t}+u_{x}=0}
MUSCL タイプの左と右の状態の線形外挿の例。 MUSCLベースの数値計算スキームは、傾き制限の ある左右の外挿状態を用いることで、各セルに線形区分近似を用いるという考え方を拡張する 。これにより、以下の高解像度TVD離散化スキームが得られる。
d u i d t + 1 Δ x i [ F ( u i + 1 / 2 ∗ ) − F ( u i − 1 / 2 ∗ ) ] = 0. {\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.} これをもっと簡潔に書くと、
d u i d t + 1 Δ x i [ F i + 1 / 2 ∗ − F i − 1 / 2 ∗ ] = 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 次近似の非線形組み合わせに対応します。 F i ± 1 / 2 ∗ {\displaystyle F_{i\pm 1/2}^{*}}
記号 およびは、 （限定された外挿されたセルエッジ変数の）スキーム依存関数を表す。 すなわち 、 u i + 1 / 2 ∗ {\displaystyle u_{i+1/2}^{*}} u i − 1 / 2 ∗ {\displaystyle u_{i-1/2}^{*}}
u i + 1 / 2 ∗ = u i + 1 / 2 ∗ ( u i + 1 / 2 L , u i + 1 / 2 R ) , u i − 1 / 2 ∗ = u i − 1 / 2 ∗ ( u i − 1 / 2 L , u i − 1 / 2 R ) , {\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),} ここで、風下斜面を使用すると、
u i + 1 / 2 L = u i + 0.5 ϕ ( r i ) ( u i + 1 − u i ) , u i + 1 / 2 R = u i + 1 − 0.5 ϕ ( r i + 1 ) ( u i + 2 − u i + 1 ) , {\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),} u i − 1 / 2 L = u i − 1 + 0.5 ϕ ( r i − 1 ) ( u i − u i − 1 ) , u i − 1 / 2 R = u i − 0.5 ϕ ( r i ) ( u i + 1 − u i ) , {\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),} そして
r i = u i − u i − 1 u i + 1 − u i . {\displaystyle r_{i}={\frac {u_{i}-u_{i-1}}{u_{i+1}-u_{i}}}.} 関数は 、区分近似の傾きを制限するリミッタ関数であり、解がTVDとなることを保証します。これにより、不連続点や衝撃波の周囲で発生する可能性のある不要な振動を回避します（ フラックスリミッタの セクションを参照）。リミッタは の場合にはゼロ 、 の場合には1です 。したがって、TVD離散化の精度は局所的極値では1次まで低下しますが、領域の滑らかな部分では2次まで低下する傾向があります。 ϕ ( r i ) {\displaystyle \phi \left(r_{i}\right)} r ≤ 0 {\displaystyle r\leq 0} r = 1 {\displaystyle r=1}
このアルゴリズムは実装が簡単です。Kurganov -Tadmor法 （下記参照）などの適切なスキームを選択すれば 、標準的な数値積分手法を用いて解を求めることができます。 F i + 1 / 2 ∗ {\displaystyle F_{i+1/2}^{*}}
クルガノフとタドモールの中央計画 KurganovとTadmor (KT) の中心スキーム （KurganovとTadmor、2000） の前身は、 NessyahuとTadmor (NT)のスタガード 中心スキーム （NessyahuとTadmor、1990）です。これは、MUSCL再構成を使用する、リーマンソルバーを必要とせず、2次の 高解像度スキーム です。これは完全に離散的な方法で実装が簡単で、 スカラー問題 と ベクトル問題に使用でき、高次の再構成が補足されたルサノフフラックス（局所Lax-Friedrichsフラックスとも呼ばれる）と見なすことができます。このアルゴリズムは 中心差分 に基づいており、 高勾配現象を示すシステムを記述するPDEの解を得るために使用した場合、リーマン型ソルバーに匹敵するパフォーマンスを発揮します。
KT法はNT法を拡張したもので、元のNT法よりも数値粘性が低くなっています。また、完全離散法 または 半離散 法のいずれでも実装できるという利点もあります 。ここでは半離散法について考察します。
計算は以下の通りです。
F i − 1 2 ∗ = 1 2 { [ F ( u i − 1 2 R ) + F ( u i − 1 2 L ) ] − a i − 1 2 [ u i − 1 2 R − u i − 1 2 L ] } . {\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\}.} F i + 1 2 ∗ = 1 2 { [ F ( u i + 1 2 R ) + F ( u i + 1 2 L ) ] − a i + 1 2 [ u i + 1 2 R − u i + 1 2 L ] } . {\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\}.} 1次元の移流方程式 。ステップ波が右方向に伝播する。解析解と、SuperBeeリミッターを用いたKurganov-Tadmor中心スキームに基づくシミュレーションを示す。 u t + u x = 0 {\displaystyle u_{t}+u_{x}=0} ここで、 局所伝播速度は 、 次式で与えられる セル上 のヤコビアン固有値の絶対値の最大値である。 a i ± 1 2 {\displaystyle a_{i\pm {\frac {1}{2}}}\ } F ( u ( x , t ) ) {\displaystyle F\left(u\left(x,t\right)\right)} i , i ± 1 {\displaystyle {i},{i\pm 1}}
a i + 1 2 ( t ) = max [ ρ ( ∂ F ( u i + 1 / 2 L ( t ) ) ∂ u ) , ρ ( ∂ F ( u i + 1 / 2 R ( t ) ) ∂ u ) , ] {\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]} スペクトル半径 を 表す 。 ρ ( ∂ F ( u ( t ) ) ∂ u ) {\displaystyle \rho \left({\frac {\partial F\left(u\left(t\right)\right)}{\partial u}}\right)\ } ∂ F ( u ( t ) ) ∂ u . {\displaystyle {\frac {\partial F\left(u\left(t\right)\right)}{\partial u}}.}
これらの CFL 関連の速度を超える特性情報は必要ありません。
上記のフラックス計算は、通常、 Lax-Friedrichs フラックス と呼ばれます(ただし、このようなフラックス表現は Lax (1954) ではなく Rusanov (1961) に登場していることは注目に値します)。
高解像度スキームを使用する有効性の例が、 右に伝播するステップ波を伴う 1D 移流方程式を示した反対の図に示されています。シミュレーションは、Kurganov および Tadmor の中心スキームを Superbee リミッター とともに使用し、時間積分に RK-4 を使用し、200 セルのメッシュで実行されました。このシミュレーション結果は、上に示した 1 次風上差分および 2 次中心差分の結果と非常によく対照的です。このスキームは、方程式のセットに適用した場合にも良好な結果をもたらします。このスキームをオイラー方程式に適用した結果については、以下を参照してください。ただし、たとえば Superbee リミッターは一部の滑らかな波に対して非現実的な鋭さを引き起こす可能性があるため、適切なリミッターを選択する際には注意が必要です。 u t + u x = 0 {\displaystyle u_{t}+u_{x}=0\ }
このスキームは、拡散項が存在する場合には容易に組み込むことができる。例えば、上記の1次元スカラー問題を拡散項を含むように拡張すると、以下のようになる。
u t + F x ( u ) = Q x ( u , u x ) , {\displaystyle u_{t}+F_{x}\left(u\right)=Q_{x}\left(u,u_{x}\right),} これに対してクルガノフとタドモアは次のような中心差分近似を提案している。
d u i d t = − 1 Δ x i [ F i + 1 2 ∗ − F i − 1 2 ∗ ] + 1 Δ x i [ P i + 1 2 − P i − 1 2 ] . {\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].} どこ、
P i + 1 2 = 1 2 [ Q ( u i , u i + 1 − u i Δ x i ) + Q ( u i + 1 , u i + 1 − u i Δ x i ) ] , {\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],} P i − 1 2 = 1 2 [ Q ( u i − 1 , u i − u i − 1 Δ x i − 1 ) + Q ( u i , u i − u i − 1 Δ x i − 1 ) . ] {\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 and Tadmor, 2000）に記載されており、1次元および2次元の例も多数掲載されています。また、Nessyahu and Tadmor (1990) による関連論文にも追加情報が掲載されています。
注：このスキームは、もともとKurganovとTadmorによって 線形外挿 に基づく2次スキームとして提示されました 。後の論文（Kurganov and Levy, 2000）では、このスキームが3次スキームの基礎にもなり得ることが示されています。このスキームの1次元移流の例と、放物面再構成（3次）を用いたオイラー方程式の例は、以下の 放物面再構成 と オイラー方程式の セクションに示されています。
区分放物線再構成 MUSCL 型状態放物線再構成の例。 線形外挿の考え方を高次の再構成に拡張することが可能であり、その例を右図に示します。ただし、この場合、左と右の状態は、2次の風上バイアス差分方程式の補間によって推定されます。これにより、空間的に3次の精度を持つ放物線状の再構成スキームが得られます。
我々はKermani (Kermani, et al., 2003)のアプローチに従い、3次の風上バイアススキームを提示する。ここでも記号 と は スキーム依存関数（限定的に再構成されたセル端変数）を表す。ただし、この場合、それらは放物線的に再構成された状態、すなわち に基づいて いる 。 u i + 1 2 ∗ {\displaystyle u_{i+{\frac {1}{2}}}^{*}} u i − 1 2 ∗ {\displaystyle u_{i-{\frac {1}{2}}}^{*}}
u i + 1 2 ∗ = f ( u i + 1 2 L , u i + 1 2 R ) , u i − 1 2 ∗ = f ( u i − 1 2 L , u i − 1 2 R ) , {\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),} そして
u i + 1 2 L = u i + ϕ ( r i ) 4 [ ( 1 − κ ) δ u i − 1 2 + ( 1 + κ ) δ u i + 1 2 ] , {\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],} u i + 1 2 R = u i + 1 − ϕ ( r i + 1 ) 4 [ ( 1 − κ ) δ u i + 3 2 + ( 1 + κ ) δ u i + 1 2 ] , {\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],} u i − 1 2 L = u i − 1 + ϕ ( r i − 1 ) 4 [ ( 1 − κ ) δ u i − 3 2 + ( 1 + κ ) δ u i − 1 2 ] , {\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],} u i − 1 2 R = u i − ϕ ( r i ) 4 [ ( 1 − κ ) δ u i + 1 2 + ( 1 + κ ) δ u i − 1 2 ] . {\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次元の移流方程式 。ステップ波が右方向に伝播する。解析解と、放物線再構成とファン・アルバダ・リミッターを用いたクルガノフ・タドモア中心スキームに基づくシミュレーション結果を示す。 u t + u x = 0 {\displaystyle u_{t}+u_{x}=0} ここで = 1/3 であり、 κ {\displaystyle \kappa \ }
δ u i + 1 2 = ( u i + 1 − u i ) , δ u i − 1 2 = ( u i − u i − 1 ) , {\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),} δ u i + 3 2 = ( u i + 2 − u i + 1 ) , δ u i − 3 2 = ( u i − 1 − u i − 2 ) , {\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),} リミッター関数は 上記と同じです。 ϕ ( r ) {\displaystyle \phi \left(r\right)\ }
放物面再構成は実装が簡単で、上記に示した線形外挿の代わりに Kurganov および Tadmor 方式で使用できます。これは、KT 方式の空間解を 3 次まで上げる効果があります。これは、オイラー方程式を解くときに優れたパフォーマンスを発揮します (以下を参照)。この空間次数の増加は、滑らかな解の場合に 2 次方式よりも有利ですが、衝撃の場合はより消散的です。反対側の図と、線形外挿と Superbee リミッタを使用した KT アルゴリズムを使用して得られた上記の解を比較してください。このシミュレーションは、同じ KT アルゴリズムを使用して放物面再構成を行い、200 セルのメッシュで実行されました。時間積分は RK-4 によって行われ、 スプリアス振動を回避するために、van Albada リミッタの代替形式 が使用されました。 ϕ v a ( r ) = 2 r 1 + r 2 {\displaystyle \phi _{va}(r)={\frac {2r}{1+r^{2}}}\ }
例: 1次元オイラー方程式 簡単のため、熱伝達と体積力を考慮しない1次元の場合を考える。したがって、ベクトル保存形式では、一般 オイラー方程式は 次のように簡約される。
∂ U ∂ t + ∂ F ∂ x = 0 , {\displaystyle {\frac {\partial \mathbf {U} }{\partial t}}+{\frac {\partial \mathbf {F} }{\partial x}}=0,} どこ
U = ( ρ ρ u E ) F = ( ρ u p + ρ u 2 u ( E + p ) ) , {\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 } ここで は 状態のベクトルであり、 はフラックスのベクトルです。 U {\displaystyle {\mbox{U}}} F {\displaystyle {\mbox{F}}}
上記の式は、 質量 、 運動量 、 エネルギー の保存則を表しています。したがって、3つの式と4つの未知数（密度 、流速、 圧力、全エネルギー） があります 。全エネルギーは、 ρ {\displaystyle \rho } u {\displaystyle u} p {\displaystyle p} E {\displaystyle E}
E = ρ e + 1 2 ρ u 2 , {\displaystyle E=\rho e+{\frac {1}{2}}\rho u^{2},} ここで、 特定の内部エネルギーを表します。 e {\displaystyle e\ }
系を閉じるには 状態方程式 が必要です。私たちの目的に適した状態方程式は
p = ρ ( γ − 1 ) e , {\displaystyle p=\rho \left(\gamma -1\right)e,} ここで、 は 流体の 比熱の比に等しくなります。 γ {\displaystyle \gamma \ } [ c p / c v ] {\displaystyle \left[c_{p}/c_{v}\right]}
上の簡単な1次元の例で示したように、各状態変数について左と右の外挿状態を求めることで、密度を求めることができます。
ρ i + 1 2 ∗ = ρ i + 1 2 ∗ ( ρ i + 1 2 L , ρ i + 1 2 R ) , ρ i − 1 2 ∗ = ρ i − 1 2 ∗ ( ρ i − 1 2 L , ρ i − 1 2 R ) , {\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),} どこ
ρ i + 1 2 L = ρ i + 0.5 ϕ ( r i ) ( ρ i + 1 − ρ i ) , ρ i + 1 2 R = ρ i + 1 − 0.5 ϕ ( r i + 1 ) ( ρ i + 2 − ρ i + 1 ) , {\displaystyle \rho _{i+{\frac {1}{2}}}^{L}=\rho _{i}+0.5\phi \left(r_{i}\right)\left(\rho _{i+1}-\rho _{i}\right),\quad \rho _{i+{\frac {1}{2}}}^{R}=\rho _{i+1}-0.5\phi \left(r_{i+1}\right)\left(\rho _{i+2}-\rho _{i+1}\right),} ρ i − 1 2 L = ρ i − 1 + 0.5 ϕ ( r i − 1 ) ( ρ i − ρ i − 1 ) , ρ i − 1 2 R = ρ i − 0.5 ϕ ( r i ) ( ρ i + 1 − ρ i ) . {\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).} 同様に、運動量 、全エネルギーについても、 速度 は 運動量から計算され、圧力 は 状態方程式から計算されます。 ρ u {\displaystyle \rho u} E {\displaystyle E} u {\displaystyle u} p {\displaystyle p}
限定された外挿状態が得られたら、これらの値を用いてエッジフラックスを構築します。エッジフラックスが既知であれ ば、半離散スキームを構築できます 。
d U i d t = − 1 Δ x i [ F i + 1 2 ∗ − F i − 1 2 ∗ ] . {\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法など）も選択する必要があります 。 F i ± 1 2 ∗ {\displaystyle \mathbf {F} _{i\pm {\frac {1}{2}}}^{*}}
GA Sodの「衝撃波管」問題に基づくオイラー方程式の高解像度シミュレーション。線形外挿とOspreリミッタを用いたKuganov-Tadmor中央スキームに基づく解析解とシミュレーション（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)。 GA Sodの「衝撃波管」問題に基づくオイラー方程式の高解像度シミュレーション（SI単位系）。KurganovとTadmorの中央法に基づき、放物線再構成とファン・アルバダ・リミッターを適用した解析解とシミュレーション（3次）解を示します。 右の図は、上記の高解像度の Kurganov および Tadmor 中心スキーム (KT) を使用し、放物線再構成とファン アルバダ リミッタを使用した、GA Sod の 衝撃波管 問題 (Sod、1978) の 3 次解を示しています。これもまた、オイラー方程式を解くための MUSCL アプローチの有効性を示しています。シミュレーションは、放物線外挿とファン アルバダ リミッタ を使用した KT アルゴリズムを使用するように適合させた Matlab コード (Wesseling、2001) を使用して、200 セルのメッシュで実行されました 。 不要な振動を避けるため、ファン アルバダ リミッタの代替形式 を使用しました。時間積分は 4 次の SHK 積分器によって実行されました。同じ初期条件が使用されました。 ϕ v a ( r ) = 2 r 1 + r 2 {\displaystyle \phi _{va}(r)={\frac {2r}{1+r^{2}}}\ }
オイラー方程式を高精度に解く様々な高解像度法が開発されている。そのような法の例として、
オッシャー 計画 、そして Liou -Steffen AUSM （移流上流分割法）スキーム。 これらの手法やその他の手法に関する詳細は、以下の参考文献をご覧ください。KurganovとTadmorの中央スキームのオープンソース実装は、以下の外部リンクからご覧いただけます。
参照
参考文献
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). 双曲型問題のための有限体積法 , Cambridge University Press. van Leer, B. (1979)「究極の保存的差分スキームに向けて、V. ゴドゥノフ法の二次続編」 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 で書かれた、クルガノフとタドモアの中央スキームを使用してオイラー方程式を解くオープンソースコード （作者：アルノ・マイヤーホファー）
![{\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.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/373bbe067f44f6824dde758f2ac0fb99f9056298)
![{\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}].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cc785f74a7dc970ff741ab6468983b943d862490)
![{\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,}](https://wikimedia.org/api/rest_v1/media/math/render/svg/3e679a445e736fc2d97ad91fbeefb3e2228eed14)
![{\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.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e4f0432d7d1f2313d6101629fcee401728e4a8b1)
![{\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.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bbcd1b79823bd121c764772d1acf9a77b22eca69)
![{\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\}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/797a23caf1cdab28cdb411ebe1022fa0dba4b535)
![{\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\}}.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b72f9a973efd1827ab3f0a304fef575ecedf0a80)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/e8c3a65255e7c37e4efb8c43541b1c62028c98a3)
![{\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].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/31422d5b11311e551a113a34d9be14d793c8c4c6)
![{\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],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7f52e9f99931a16d0de065a0fbc2f73f7d6510f6)
![{\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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1282687e8bc1b0d75d5b9e86b779a92e0c0b08a4)
![{\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],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/bfec74720862bfcee4e9304b77ff6861b56423a1)
![{\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],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/05ff2d71636cbcb7ecf52c9c6258600936856f7a)
![{\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],}](https://wikimedia.org/api/rest_v1/media/math/render/svg/88d5efa7ed3fb256c51e6b51f471311973900ef4)
![{\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].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/989a21adc92174ed95e3e9ce9b5bfa380e91f42f)
![{\displaystyle \left[c_{p}/c_{v}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/41e20eb99b7ea7d835080c9d18cce0d4c6bf3218)
![{\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].}](https://wikimedia.org/api/rest_v1/media/math/render/svg/deb9011df41d0199f1809f2a59847eb9c523fdda)
