第6-3回 流束制限関数(minmod関数)[python]
この第6回を通して高次精度化を行っており,時間,空間についての高次精度化を説明してきました.
今回は流束すなわちフラックスの高次精度化を考えており,『一次精度風上差分から精度を上げたい』と考えています.
しかし,以下の事を証明したGodunovの定理と呼ばれる定理が存在します.
「\( \frac{\partial u}{\partial t} + c\frac{\partial u}{\partial x}=0\)に対して,次の形をした2次以上の精度を持つどのようなスキームも解の単調性を維持することができない」
$$u_j^{n+1}= \sum _k c_k u_{j+k}^n$$
簡単に言いますと,「全ての領域で同じスキーム(一次風上または中心差分など)を用いると2次精度以上では数値振動してしまう」ということになります.
数値振動は次のように値が離散的に分散して,値が大きくなっていくことを指します(0.6付近で振動してる).
つまり,フラックスを高次精度にしたい場合,周りの情報からどのスキームを用いるか決める非線形のスキームでなければなりません.
以下では一次風上差分とLax-Wendroffスキームと呼ばれる二次精度のスキームについて考えます.(Lax-Wendroffスキームについては他の章で導出等詳しく述べたいと思います.)
一次精度風上差分 $$ \ \ \ \tilde{f}_{j+1/2}=cu_j$$
Lax-Wendroff $$ \ \ \ \tilde{f}_{j+1/2}=\frac{c}{2}(1+\nu)u_j+\frac{c}{2}(1-\nu)u_{j+1} $$
$$\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ =c\left[ u_j+\frac{1}{2}(1-\nu)(u_{j+1}-u_j)\right]$$
ここで,\(\nu=c\frac{\Delta t}{\Delta x}\)です.
2つのスキームを比べると第2項が存在するかどうかになります.そこで次のように\(B_{j+1/2}\)を定義します.
$$\tilde{f}_{j+1/2}=c\left[ u_j+\frac{1}{2}(1-\nu)B_{j+1/2}(u_{j+1}-u_j)\right]$$
この\(B_{j+1/2}\)は流速制限関数と呼ばれ,選び方によっては一次風上差分にもLax-Wendroffスキームにもなり, \(B_{j+1/2}=0\)で一次風上, \(B_{j+1/2}=1\)で二次精度のLax-Wendroffスキームになります.この\(B\)はそれぞれのセル境界で値が違うため添え字が付けられています.
有限体積法で用いた次の形の式に代入し,変形していきます.
$$\frac{u_j^{n+1}-u_j^n}{\Delta t}=-\frac{\tilde{f}_{j+1/2}-\tilde{f}_{j-1/2} }{\Delta x}$$
$$⇔ u_j^{n+1}-u_j^n =-c\frac{\Delta t}{\Delta x} ( u_j+\frac{1}{2}(1-\nu)B_{j+1/2}(u_{j+1}-u_j) $$ $$- u_{j-1}-\frac{1}{2}(1-\nu)B_{j-1/2}(u_{j}-u_{j-1}))$$
$$⇔ u_j^{n+1}-u_j^n =-\nu(u_j^{n}-u_{j-1}^n)-\frac{1}{2}\nu(1-\nu)({B_{j+1/2}(u^n_{j+1}-u^n_j) \\ – B_{j-1/2}(u^n_{j}-u^n_{j-1}))}$$
$$⇔ \frac {u_j^{n+1}-u_j^n}{u^n_{j-1}-u_j^n} =\nu+\frac{1}{2}\nu(1-\nu){\frac{B_{j+1/2}}{r_j}- B_{j-1/2}}$$
ここで\(r_j\)は次の通りです.
$$r_j=\frac{u_{j-1}^{n}-u^n_{j}}{u^n_{j}-u^n_{j+1}}$$
数値振動を避けるための十分条件は「\(0\le(左辺)\le1\)」になります.
風上差分である後退差分では,風上から情報が伝わるため,\(u_j^{n}<u_j^{n+1}<u_{j-1}^{n}\)になります.
この式を解くと十分条件「\(0\le(左辺)\le1\)」が求まります.
$$ 0 < \frac{u_j^{n+1} – u_j^{n}}{u_{j-1}^{n} – u_j^{n}} < 1$$
左側の不等式を解いていきます.第7回で述べますが,安定的に計算を行う条件は\(0\le\nu\le1\)であり,式変形にこの条件を用います.
$$ 0 \le\nu+\frac{1}{2}\nu(1-\nu)(\frac{B_{j+1/2}}{r_j}- B_{j-1/2}) $$
$$⇔-\nu\le\frac{1}{2}\nu(1-\nu)(\frac{B_{j+1/2}}{r_j}- B_{j-1/2}) $$
$$⇔-\frac{2}{1-\nu}\le\frac{B_{j+1/2}}{r_j}- B_{j-1/2} $$
$$⇔B_{j-1/2}-\frac{B_{j+1/2}}{r_j} \le\frac{2}{1-\nu} $$
右側の不等式を解いていきます.
$$\nu+\frac{1}{2}\nu(1-\nu)(\frac{B_{j+1/2}}{r_j}- B_{j-1/2})\le1$$
$$⇔\frac{1}{2}\nu(1-\nu)(\frac{B_{j+1/2}}{r_j}- B_{j-1/2})\le1-\nu$$
$$⇔\frac{B_{j+1/2}}{r_j}- B_{j-1/2}\le\frac{2}{\nu}$$
$$⇔-\frac{2}{\nu}\le B_{j-1/2}-\frac{B_{j+1/2}}{r_j}$$
以上の結果から
$$-\frac{2}{\nu}\le B_{j-1/2}-\frac{B_{j+1/2}}{r_j}\le\frac{2}{1-\nu} $$
この条件を満たす十分条件を探します.\(0\le\nu\le1\)であるので
$$-\frac{2}{\nu}\le -2,\ \ \ \ \ 2\le\frac{2}{1-\nu} $$
よって,
$$-2\le B_{j-1/2}-\frac{B_{j+1/2}}{r_j}\le2 $$
さらに正負を考慮すると,次の2式を満たせば良いことになります.
$$0\le B_{j-1/2}\le2, \ \ \ \ \ 0\le\frac{B_{j+1/2}}{r_j}\le2 \ \ \ \ \ \ (1)$$
また流れがスムーズな領域では\(r_j\)は1に近く,この領域では二次精度を維持したいので\(B_{j+1/2}\)も1に近いと良い,つまり下記の点を通ると良いです.
$$(r_j, B_{j+1/2})=(1,1) \ \ \ \ \ \ (2)$$
この式を図示すると次のようになります.
\(r_j\le0\)で\(B_{j+1/2}=0\),\(0\le r_j\)で水色の領域内であれば条件を満たし,\((r_j, B_{j+1/2})=(1,1)\)を通ると良い.
これらの条件を満たす関数は複数考えられますが,今回はminmod関数を用います.
$$\begin{eqnarray} minmod(r) = \begin{cases}{}
0 \ \ \ \ \ (r \le 0) &\\ r \ \ \ \ \ (0 \le r \le 1) &\\ 1 \ \ \ \ \ (1 \le r) \end{cases}\end{eqnarray}$$
次回はMUSCL法とminmod関数を組み合わせる方法について行います.
*第6-5回にコードは書いてあります。
By hide
Tweet
コメントを残す