工学系大学院生のブログ

2人の大学院生による雑記ブログ

第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

このエントリーをはてなブックマークに追加

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です