第6-2回　MUSCL法（有限体積法）[python]

今回は空間高次精度化にあたり、MUSCL法（ Monotone Upwind Scheme for Conservation Laws ）を用います。

この方法の考え方をまず述べます。

有限体積法では検査体積をとり、そのセル内で一定の値を取るとして考えています。

MUSCL法では、セル内で一次関数や二次関数の分布があると考え、高次精度化します。

上の図のようにセル内の値を考えると、左から一次精度、二次精度、三次精度で表すことができます。

セル内の物理量の分布を求めるため、物理量\(u\)の分布関数を\(x_{j-1/2} \le x \le x_{j+1/2}\)の範囲でテーラー展開します。

$$u(x)=u(x_j)+(x-x_j)u^{\prime}(x_j)+\frac{1}{2}(x-x_j)^2u ^{\prime \prime} (x_j)$$ $$+ \frac{1}{6}(x-x_j)^3u ^{\prime \prime \prime} (x_j)+O(x^4) \ \ \ \ \ \ \ \dots ①$$

MUSCL法では、セルでの平均値をセルの値とします。つまり \(x_{j-1/2} \le x \le x_{j+1/2}\) での平均をセルの値とします。

$$u_j=\frac{1}{\Delta x} \displaystyle\int _ {x_{j-1/2}} ^{x_{j+1/2}} u(x)dx$$

$$=\frac{1}{\Delta x} [u(x_j)x+ \frac{1}{2}(x-x_j)^2u^{\prime}(x_j)+\frac{1}{6}(x-x_j)^3u ^{\prime \prime} (x_j) $$ $$+ \frac{1}{24}(x-x_j)^4u^{\prime \prime \prime }(x_j)+O(x^5) ] ^{x_{j+1/2}}_ {x_{j-1/2}} $$

$$=\frac{1}{\Delta x} \left(\Delta x u(x_j)+\frac{1}{6} \left(\frac{\Delta x}{2} \right)^3 u^{ \prime \prime }(x_j) \times 2+O(\Delta x^5)\right)$$

$$⇔u_j=u(x_j)+\frac{\Delta x^2}{24}u^{ \prime \prime }(x_j)+O(\Delta x^4)$$

さらに①から\(u(x)\)の一階微分、 二階微分 は次の通りになります。

$$u^{\prime}(x)= u^{\prime}(x_j)+(x-x_j)u ^{\prime \prime} (x_j)$$ $$+ \frac{1}{2}(x-x_j)^2u ^{\prime \prime \prime} (x_j)+O(\Delta x^3) $$

$$u^{\prime \prime}(x)=u ^{\prime \prime} (x_j)+ (x-x_j)u ^{\prime \prime \prime} (x_j)+O(\Delta x^2) $$

また、セルでの平均値\(u^{\prime}_j\), \(u^{\prime \prime}_j\) は積分を行うことで、次のようになります。

$$u^{\prime}_j= \frac{1}{\Delta x}\int ^{x_{j+1/2}}_ {x_{j-1/2}} u^{\prime}(x)dx = u^{\prime}(x_j)+\frac{\Delta x^2 }{24}u ^{\prime \prime \prime} (x_j)+O(\Delta x^4) $$ $$= u^{\prime}(x_j)+O(\Delta x^2)$$

$$ u^{\prime \prime}_j= \frac{1}{\Delta x}\int ^{x_{j+1/2}}_ {x_{j-1/2}} u^{\prime \prime}(x)dx =u^{\prime \prime}(x_j)+O(\Delta x^2)$$

最後に、これらの結果を①に代入します。これでセル内の分布が定義できます。

$$u(x)=u(x_j)+(x-x_j)u^{\prime}(x_j)+\frac{1}{2}(x-x_j)^2u ^{\prime \prime} (x_j)$$ $$+ \frac{1}{6}(x-x_j)^3u ^{\prime \prime \prime} (x_j)+O(\Delta x^4) $$

$$u(x)=(u_j-\frac{\Delta x^2}{24}u^{\prime \prime}(x_j))+(x-x_j)(u^{\prime}(x_j))+\frac{1}{2}(x-x_j)^2(u ^{\prime \prime}_j ) $$ $$+ \frac{1}{6}(x-x_j)^3u ^{\prime \prime \prime} (x_j)+O(\Delta x^4) $$

$$u(x)=u_j+(x-x_j)u^{\prime}_j +\frac{1}{2}\left [ (x-x_j)^2-\frac{\Delta x^2}{12}\right]u^{\prime \prime}_j +O(\Delta x^3) $$

MUSCL法の肝は、このように分布を持たせ境界付近の値を算出できることです。この値を用いてフラックスを決めることで精度が上がります。

続きまして、セル境界付近での値を求めます。

今までの式ではセル\(j\)について行いました。この式の右端、すなわち境界\(j+1/2\)の左側の\(u_L\)を求めます。境界に対して右左で名前を付けるので、気を付けてください。

セル境界での物理量は\(u^L_{j+1/2}=u^j (x_{j+1/2}),
u^R_{j+1/2}=u^{j+1}(x_{j-1/2}) \)と定義し代入し、計算を行います。

$$u^L_{j+1/2}=u_j+\frac{\Delta x}{2}u_j^{\prime}+
\frac{\Delta x^2}{12}u_j^{\prime \prime} +O(\Delta x^3)$$

$$u^R_{j+1/2}=u_{j+1}-\frac{\Delta x}{2}u_{j+1}^{\prime}+
\frac{\Delta x^2}{12}u_{j+1}^{\prime \prime} +O(\Delta x^3)$$

\(k\)を用いて、この式を一般化します。このページの最後に書きましたが、この\(k\)を変えることで、一次精度から三次精度まで自由に変えることができます。

$$u^L_{j+1/2}=u_j+\frac{\Delta x}{2}u_j^{\prime}+ \frac{k\Delta x^2}{4}u_j^{\prime \prime} $$

$$u^R_{j+1/2}=u_{j+1}-\frac{\Delta x}{2}u_{j+1}^{\prime}+ \frac{k\Delta x^2}{4}u_{j+1}^{\prime \prime} $$

差分近似を用いて\(u^{\prime}_j\)と\(u^{\prime \prime}_j\)を表します。

$$ u^{\prime}_j=\frac{u_{j+1}-u_{j-1}}{2\Delta x} $$ $$ u^{\prime \prime}_j =\frac{ u_{j+1} -2 u_{j} + u_{j-1} }{\Delta x^2}$$

よって

$$u^L_{j+1/2}=u_j+\frac{1}{4}( u_{j+1}-u_{j-1} )+ \frac{k}{4}( u_{j+1}-2u_j+u_{j-1} ) $$ $$= u_j+\frac{1}{4}(({ u_{j+1}-u_j)+(u_j-u_{j-1}) + k( u_{j+1}-u_j)-k(u_j-u_{j-1} )}) $$ $$ = u_j+\frac{1-k}{4}(u_j-u_{j-1})+ \frac{1+k}{4}(u_{j+1}-u_{j}) $$

$$ u^L_{j+1/2} = u_j+\frac{1-k}{4}\Delta ^-_j+ \frac{1+k}{4} \Delta ^+_j $$

ここで、\( \Delta ^+_j =u_{j+1}-u_{j}, \ \Delta ^-_j =u_{j}-u_{j-1} \)と置きました。

同様に \(u^R_{j+1/2} \)は

$$ u^R_{j+1/2} = u_{j+1}-\frac{1-k}{4}\Delta ^+_{j+1}- \frac{1+k}{4} \Delta ^-_{j+1} $$

最後にこれらの値を用いて、フラックスを定義します。

第5-2回で用いたFVSでは次の式のように、一次風上を用いていたので\(A^+\)に\(A^+_j\)、 \(A^-\)に\(A^-_{j+1}\) を適用しました。

$$ \tilde{f}_{j+1/2}=F_j^++F_{j+1}^- =A_j^+Q_j+A_{j+1}^-Q_{j+1} $$

考え方は同じですがより境界に近い値を用います。つまり、\(A^+\)に\((A^+_{L})_{j+1/2}\)、\(A^-\)に\((A^-_R)_{j+1/2}\)を適用します。

$$ \tilde{f}_{j+1/2}=(F_L^+)_{j+1/2}+(F_R^-)_{j+1/2} = (A^+_{L}Q_L)_{j+1/2} + (A^+_{L}Q_R)_{j+1/2} $$

これで、MUSCL法の適用はできました。

このMUSCL法は基本的に精度を上げるために導入しているので\(k\)の値は三次精度となる\(1/3\)にしますが、\(k\)の値によってどのようになるか、以下に書いておきます。

\(k=1/3\) : 三次精度
\(k=0\) : 二次精度風上差分
\(k=1\): 一次精度の代数平均、境界での不連続面なし
\(k=-1\):片側のみを用いた外挿

次回は流束制限関数というものについて行います。

*第6-5回にコードは書いてあります。

by hide