Navier-Stokes方程式での移流項 (u・▽)u について



移流方程式

移流方程式

移流方程式

空間微分の離散化

風上差分(Upwind differencing)

移流方程式の時間微分を前進オイラー法で離散化すると,

upwind.eq1.gif

1次元の場合を考える.グリッドupwind.eq2.gif(座標upwind.eq3.gif)について,

upwind.eq4.gif

ここで,upwind.eq5.gif.このupwind.eq6.gifを差分を使って求める.

風上差分ではその名の通り,風上側の値を使って差分を行う1次精度の離散化法である. 1次元の場合,upwind.eq7.gifでは座標軸上の右から左に風が吹いているので,風上はi+1, upwind.eq8.gifでは左から右に風が吹いているので,風上はi-1となる.

upwind.eq9.gif

風上差分はCFL(Courant-Friedreichs-Lewy)条件を満たしていれば安定である. CFL条件によるタイムステップ幅の制限は,

upwind.eq10.gif

タイムステップ幅は通常これよりもさらに小さい数を用いた方がよい. 係数:CFL数(CFL number)をupwind.eq11.gifとすると,

upwind.eq12.gif

ここで,upwind.eq13.gif.よく使われるのはupwind.eq14.gifである.

振動と数値拡散

風上差分の離散式(upwind.eq8.gifの場合)を位置upwind.eq3.gifでテーラー展開してみる.

upwind.eq15.gif

upwind.eq16.gifは空間2階微分であり,これは拡散を表している. 拡散に関する項は,例えば,流体のナビエ・ストークス方程式にもある(upwind.eq17.gif). しかし,これはそれらの項とは異なり,離散化によって生まれた拡散であり,これを数値拡散(numerical diffusion)と呼ぶ. 下の矩形波の移流を見ると風上差分において数値拡散が発生しているのがよく分かる.

さて,上式の微分を下で述べる中心差分で離散化してみる.

upwind.eq18.gif

風上差分は中心差分に拡散項を加えたものであることが分かる. 中心差分では非常に大きな振動が発生する. その振動を拡散項により抑えているのが風上差分である. この振動を抑えるという考え方は重要である. より高次の関数を使って補間することで振動を抑えるのが,Lax-WendroffやQUICK,QUICKEST,河村・桑原スキームなどで, 2階,3階の微分を使って抑えるのがENOやWENOである. これらに対して,元の形状を保持するという考え方に基づき,異なるアプローチをとるのがCIP法,RCIP法などである. ちなみにRCIP法では数値拡散が発生してしまうので,これを抑えるためのSTAA手法では逆に数値拡散項を引くことで拡散を抑えている.

  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    風上差分では下の中心差分のように振動が全くない代わりに,形状がなまっている. これは数値拡散と呼ばれる現象によるためである.
    rectan_upwind.jpg
    矩形波の風上差分による移流

中心差分(Central differencing)

流れの速度central.eq1.gifに関係なく,グリッド点を中心とした周囲のグリッドの値の平均で空間微分を近似するのが中心差分である. つまり,1次元の場合,以下となる.

central.eq2.gif

前進オイラー法と組み合わせると,

central.eq3.gif
  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    風上差分ではt=5まで移流させた結果のみ示したが,中心差分では非常に振動が大きくなるため, t=0.5, 1.0, 5.0の3つのグラフを示した.流れの後方に多くの振動が現れ成長していく様子が分かる.
    rectan_central1.jpg
    矩形波の中心差分による移流(t=0.5)
    rectan_central2.jpg
    矩形波の中心差分による移流(t=1.0)
    rectan_central3.jpg
    矩形波の中心差分による移流(t=5.0)

Lax-Wendroff法

Lax-Wendroff法(または,ライス(Leith)法)*3では線形補間の代わりに2次多項式を使った補間で近似を行う.2次多項式の係数を求めるために,グリッドiとその両隣(i-1,i+1)の3つの値Lax-Wendroff.eq1.gifを用いる.

2次多項式は,

Lax-Wendroff.eq2.gif

ここで各係数は,

Lax-Wendroff.eq3.gif

Lax-Wendroff.eq4.gifとすると移流方程式の近似更新式が得られる.

Lax-Wendroff.eq5.gif

風上差分と同様に中心差分に拡散項を加えたものになっている.拡散項の係数が異なるのが風上差分との違いである.

  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_lax-wendroff.jpg
    矩形波のLax-Wendroff法による移流

QUICK(Quadratic Upstream Interpolation for Convective Kinematics)

QUICK*4quick.eq1.gifの内挿にquick.eq2.gifquick.eq3.gifに加えて,風上側のもう一つのquick.eq4.gifを用いる方法である.

quick.eq5.gif

これを移流方程式の離散化に用いると,

quick.eq6.gif
  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_quick.jpg
    矩形波のQUICKによる移流

QUICKEST(QUICK with estimated streaming terms)

QUICKEST(QUICK with estimated streaming terms)はLax-Wendroffの3つ(quickest.eq1.gif,quickest.eq2.gif,quickest.eq3.gif)に加えて,QUICKと同様に風上側のもう一つのquickest.eq4.gifを追加して近似する手法である(quickest.eq5.gifならi-2,quickest.eq6.gifならi+2).

quickest.eq7.gif
  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_quickest.jpg
    矩形波のQUICKESTによる移流

河村・桑原スキーム(Kawamura-Kuwabara scheme)

4次精度の中心差分(kk.eq1.gifを使用)に 拡散項として4階微分項を追加したのが,河村・桑原スキーム*5である. これにより3次精度の風上差分を得る(1次の風上差分は2次中心差分+2階微分項).

4次精度の中心差分は以下.

kk.eq2.gif

4階微分項の差分式は,

kk.eq3.gif

(4階微分の差分式の導出は4階微分を参照).

4階微分項に係数kk.eq4.gifを掛けた上で先ほどの式に追加する.

kk.eq5.gif

kk.eq6.gifの場合,

kk.eq7.gif
  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_kk.jpg
    矩形波の河村・桑原スキームによる移流

ENO(Essentially Non-Oscillatory polynomial interpolation)

ENO(Essentially Non-Oscillatory polynomial interpolation)は風上差分を改良し, 3次多項式の形で近似することで3次精度を実現した方法である. ENOの最初のアイデアは *6 で提案され, *7, *8 で数値計算に適用され, *9 でHamilton-Jacobi(HJ)方程式へ適用された(HJ ENO). 移流方程式はHJ方程式であるので,ここからは,HJ ENOについて説明する.

HJ ENOの式を述べる前に,準備として差分式を定義しておく. eqa_phii.gifのn階差分をeqa_Dinphi.gifと表記する.とすると,n=0は,

eqa_eno1.gif

となる.ここで,iはグリッド番号(座標値eqa_xi.gif).

1階差分はグリッド間(i-1/2とi+1/2)で定義される.

eqa_eno2.gif
eqa_eno3.gif

2階差分はi-1/2とi+1/2での1階差分値を使って,iで定義される.

eqa_eno4.gif

同様に3階差分は,

eqa_eno5.gif
eqa_eno6.gif

ENOでは3次多項式によりeqa_phi.gifを近似する.

eqa_eno7.gif

ここで,eqa_Qm.gifはm次項である. eqa_phix+i.gif, eqa_phix-i.gifの式がほしいので, 上式を微分して,x=eqa_xi.gifとすると,

eqa_eno8.gif
  • 1次項eqa_Q1d.gif
    eqa_phix-.gifの場合k=i-1,eqa_phix+.gifでk=iとすると,
    eqa_eno9.gif
    eqa_eno10.gif
    eqa_Q1d.gifのみの場合が風上差分に相当する.これに,2次項,3次項を加えることで,3次精度を得る.
  • 2次項eqa_Q2d.gif
    eqa_eno11.gif
    eqa_eno12.gif
    ここで,
    eqa_eno13.gif
  • 3次項eqa_Q3d.gif
    eqa_eno14.gif
    eqa_eno15.gif
    ここで,
    eqa_eno16.gif
    eqa_eno17.gif

これらによって,eqa_phixxi.gifを求め,

eqa_eno18.gif

によりeqa_phi.gifを更新する.

  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_eno.jpg
    矩形波のENOによる移流

WENO(Weighted Essentially Non-Oscillatory polynomial interpolation)

ENOを改良して,重み付き平均をとることで5次精度を実現したのがWENOである.

3次精度のENOではeqa_phix-i.gifについて, eqa_weno1.gifの6つの内,4つの値を用いて補間を行う. そして,その組み合わせは3パターン (eqa_weno2.gifeqa_weno3.gifeqa_weno4.gif)である. 実際に展開してみと以下となる(eqa_phix-.gifでは,k=i-1).

eqa_weno5.gif

1+2+4, 1+2+5, 1+3+6, 1+3+7の4パターンになるが, 以下に示すように1+2+5と1+3+6が同じなので,組み合わせは3パターンになる. それぞれさらに展開してφiの式にする.

eqa_weno6.gif

差分をeqa_weno7.gifと表記すると上記3パターンは,

eqa_weno8.gif

ENO近似は結局この3つのパターンに集約される. これらの内どれかが状況に応じて用いられる. WENO(Weighted ENO)はその名の通り,これら3つの凸結合 *10eqa_phix.gifを近似する *11. つまり,

eqa_weno12.gif

ここで,eqa_weno13.gif

問題となるのはeqa_weno14.gifの選択である. 関数が滑らかならば,eqa_weno15.gifで5次精度が得られる *12. ただし,矩形波のように滑らかでない関数だととたんに不安定になる. *13では滑らかでない関数でも安定な重みの計算法を提案している.その方法を以下に示す.

eqa_phix.gifの滑らかさ(smoothness)を以下のように定義する.

eqa_weno16.gif

先ほどの重みの組み合わせ(eqa_weno17.gif)をsmoothnessの2乗で割ることで重みを得る.

eqa_weno18.gif

ここで,

eqa_weno19.gif

eqa_weno20.gifはゼロ割を防ぐための項である. eqa_weno21.gifeqa_phi.gifが符号付距離場であることが前提で設定されている. つまり,eqa_phix.gifがほぼ1であり,eqa_weno22.gifの場合である. 最終的にeqa_weno23.gifを満たすために,正規化を行う.

eqa_weno24.gif

これらから,

eqa_weno25.gif

によりeqa_phi.gifを更新することで5次精度の結果が得られる.

  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_weno.jpg
    矩形波のWENOによる移流

セミラグランジュ(semi-Lagrangian)法

CG分野でももっともポピュラーな方法.下図のように計算点xから-u方向にΔtだけバックトレースして,その位置での速度場u(x_bt, t)を次のステップの速度場u(x, t+Δt)とする方法.下図ではバックトレースに1次関数を用いている. この方法は大きなタイムステップ幅でも安定して計算できるのが大きな特徴である.

back_trace.jpg

セミラグランジュ法ではバックトレースした位置での値を周囲の値を使って補間しなければならない. 一番簡単なのは線形補間を用いる方法である.ただし線形補間は1次精度であり,上記の風上差分と同じぐらい拡散が発生する. より高精度な補間法を用いることで拡散を抑えることができる.例えば,上記のENOやWENO,下記で示す方法などである.

  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_sl.jpg
    矩形波のセミラグランジュ法(線形補間)による移流

BFECC(Back and Forth Error Compensation and Correction)

BFECCはSemi-Lagrangian法においてバックトレースする際に,フォワードトレースにより誤差を求め, その誤差を修正することで高精度に移流を行う方法である. BFECCの最初のアイデアは *14 で提案され, 2005年にCG分野で移流方程式へ適用された *15, *16

Semi-Lagrangian法のバックトレースによるeqa_phi.gifの更新を以下の式で表す.

eqa_bfecc1.gif

Lは直線を使った1次精度のバックトレース,uは速度場を表す.
BFECCの場合,上式は以下のようになる.

eqa_bfecc2.gif

ここで,

eqa_bfecc3.gif


  • BFECCの考え方
    BFECCの概念図

BFECCの式を上の図を使って説明する.
まず,A点を計算点xとする. もし,特徴曲線に完璧に沿ってバックトレースしたならば,バックトレース点はBになるはずである. だが,実際には直線でトレースしているのでC点になっている(eqa_bfecc4.gif). このとき,B点とC点の間のeqa_phi.gifの値の違い(誤差)をeとする.
C点からさらにフォワードトレースした位置Dは,A点からの誤差2eを含んでいる(誤差なしならA点に戻る). D点での関数値をeqa_barphi.gifとすると,

eqa_bfecc5.gif

となる.上式から,誤差eは以下のように求められる.

eqa_bfecc6.gif

より正確なeqa_phin+1.gifは,eqa_phin.gifからeを引いた場eqa_bfecc7.gifより得られる.

eqa_bfecc8.gif
  • BFECCの実装
    1. eqa_phin.gifをバックトレース
      eqa_bfecc9.gif
    2. eqa_phia.gifをフォワードトレース
      eqa_bfecc10.gif
    3. 誤差eを引いたeqa_phiaa.gifを計算
      eqa_bfecc11.gif
    4. eqa_phiaa.gifをバックトレースして,eqa_phin+1.gifを算出
      eqa_bfecc12.gif
  • 例:矩形波の移流

    矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_bfecc.jpg
    矩形波のBFECCによる移流

MacCormack

参考文献1*1, 参考文献2*2

CIP法(Constrained Interpolation Profile scheme)

高次多項式を用いた補間法では補間式を得るためにより多くのグリッドを必要とする.例えば,ENO,WENOだと自身も含めて6グリッド(ENOで実際に用いられるのは4グリッド)の値が必要.風上差分だと2グリッドであったことを考えると3倍になっている.これは当然のことで,高次多項式の係数を得るためには多くの関数値を必要とするからである.このように用いるグリッドが増えると,特に境界付近での処理が難しくなってくる.これに対して,風上差分と同様に2グリッドのみを用いつつ,勾配値を利用することで急激に変化する関数でも対応したのがCIP法(Constrained Interpolation Profile scheme)*17*18*19である*20

CIP法の基本的なアイデアは,移流式をxで微分したとき,その導関数eqa_cip1.gifもまた,eqa_phi_x.gifと同様に移流するということである(ただし,eqa_cip2.gifとする).2メッシュ[i, i+1]間のプロファイルは,以下の3次補間関数で表される.

eqa_cip3.gif

グリッド上の4拘束条件eqa_cip4.gifから,

eqa_cip5.gif

この手法は非常に簡単に見えるが,その効果は非常に良好である.特に矩形波のように急激に変化する特徴を持つ関数を移流させた場合でも,その形状を維持したまま移流させることができる.

  • 矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_cip.jpg
    矩形波のCIPによる移流

RCIP法(Rational Constrained Interpolation Profile scheme)

有理関数を使ったCIP法. CIP法の欠点は関数を移流させるとき,若干のオーバーシュートがでることである. このオーバーシュートは移流が進んでも大きくはならないため問題がないように思えるが, 水と空気の二相流など,計算したい物理量の差が非常に大きい場合に, 小さなオーバーシュートが大きな問題となる. RCIP法 *21, *22 はこの問題を解決する. RCIP法はCIP法の3次関数の代わりに有理関数

eqa_rcip1.gif

を用いる. ここで,eqa_rcip2.gifはパラメータで0か1の値をとる. eqa_rcip3.gifのときはCIP法と同じであり,eqa_rcip4.gifのとき, この関数は単調性と凹凸性を保存する有理関数となる. eqa_rcip5.gifであり, eqa_rcip6.gifである. 上式の各係数は以下である.

eqa_rcip7.gif

ここで,

eqa_rcip8.gif

RCIP法以外にも,Tangent変換を用いることで オーバーシュートを抑える方法も提案されている *23. Tangent変換を用いた手法では物体界面がシャープに保たれるという特徴を持つ. これは数値流体計算上は非常に有益であるが, グラフィックスで用いた場合,グリッド解像度によっては 液体表面があまり滑らかにならず,保存性も保証されない.

RCIP法による密度関数の移流では,界面を表す0から1への変化部分が拡散し, 関数の傾きがなだらかになる. この滑らかな界面はCGでの利用においては美しいレンダリング結果を生むが, シミュレーションでは数値的な不安定性を引き起こす. 密度関数の体積移動によるSTAA手法 *24 を使うことで, 界面拡散を抑制することができる. また,拡散をある程度の幅で制御する方法 *25 もある.

  • 矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_tcip.jpg
    矩形波のTangent変換CIP法による移流
    rectan_rcip.jpg
    矩形波のRCIPによる移流

CIP-CSL法

CIP法では格子iのプロファイルeqa_phii.gifの決定のために, 4拘束条件から3次多項式の4係数を求めたが, この条件では値とその勾配値は保存されるが体積は保存されない. 保存系CIP移流法を用いることでこれを解決できる. 保存系でのCIP法の実装として,CIP-CSL2法やCIP-CSL4法, CIP-CSL3法などが提案されている.

CIP-CSL4

CIP-CSL4法*26,では,CIP法の4拘束条件に関数eqa_phi.gifの積分値eqa_rhoi-12.gifをさらに加える.

eqa_cipcsl4_1.gif

つまり,

eqa_cipcsl4_2.gif

の5拘束条件を満足させるために以下の4次多項式をプロファイルに用いるのが CIP-CSL4法である.

eqa_cipcsl4_3.gif

ここで,eqa_ai.gif,eqa_bi.gif,eqa_ci.gif,eqa_sdi.gif,eqa_ei.gifはグリッドi内での多項式補間パラメータである.

5拘束条件eqa_phii.gif,eqa_gi.gif,eqa_phii-1.gif,eqa_gi-1.gif,eqa_rhoi-12.gifから,

eqa_cipcsl4_4.gif

CIP法では算出した係数eqa_ai.gif,eqa_bi.gif,eqa_ci.gif,eqa_sdi.gif,eqa_ei.gifより, eqa_phi_x.gifeqa_g_x.gifを更新していたが, CIP-CSL4法ではさらに積分値eqa_rho.gifの更新が必要となる.

eqa_rho.gifの時間発展は,

eqa_cipcsl4_5.gif

ここで,

eqa_cipcsl4_6.gif

である.

CIP-CSL2

CIP-CSL2法*27,*28,*29は, 値eqa_phi_x.gifとその積分値eqa_cipcsl2_1.gifを用いる. 積分値eqa_D.gifに対する移流方程式は,

eqa_cipcsl2_2.gif

$D$の微分eqa_cipcsl2_3.gifと定義すると

eqa_cipcsl2_4.gif

これは1次元の移流方程式と同じであり, CIP法でのeqa_cipcsl2_5.gifと置き換えることで CIP法と同様に計算ができる.

格子点i内でのeqa_Di.gifを以下のように定義する.

eqa_cipcsl2_6.gif

eqa_Di.gifeqa_xi.gifから補間点までの積分値を表し, 3次関数で補間する.

eqa_cipcsl2_7.gif

CIP法における空間微分値gはeqa_phi.gifで置き換えられるので, [eqa_xi.gif, eqa_xi+1.gif]のプロファイルは,

eqa_cipcsl2_8.gif

eqa_D.gifの定義より,

eqa_cipcsl2_9.gif

また,

eqa_cipcsl2_10.gif

これらの条件により係数eqa_ai.gif,eqa_bi.gifは以下のように求められる.

eqa_cipcsl2_11.gif

eqa_rho.gifの時間発展は,eqa_xi.gif,eqa_xi+1.gifをバックトレースした位置eqa_cipcsl2_12.gifに挟まれた部分の積分値となる(下図).

eqa_cipcsl2_13.gif
csl2_int.gif

CIP-CSL3

CIP-CSL3法*30,*31,*32は,グリッドセル[eqa_xi-1.gif, eqa_xi.gif]間の積分値とグリッド中心eqa_xi-12.gifにおける勾配(グリッド端eqa_xi.gifにおける勾配(gradient) eqa_gi.gifと区別するためにslopeと呼び eqa_sdi.gif で表す)を用いてCIP補間を行う手法である.

CIP-CSL3法による移流項の解法では,まず,nステップでの値eqa_phin.gifから, 中間値eqa_tilde_phi.gifを求める.セミラグランジュ法による中間値の算出式(1次元)は,

eqa_cipcsl3_1.gif

である. 移流速度場の方向により以下の補間を使い分ける.

eqa_cipcsl3_2.gif

そして,

eqa_cipcsl3_3.gif

左側要素eqa_phiiL.gifでは,グリッドの端点eqa_xi.gif,eqa_xi-1.gifにおいて,上の式より,eqa_cipcsl3_4.gifである. その他の点での補間にはグリッド間の値を積分した積分値eqa_rhoi-12n.gifとグリッド中心の傾き(slope) eqa_di-12n.gifを用いる.

eqa_cipcsl3_5.gif
eqa_cipcsl3_6.gif

eqa_phiin.gif,eqa_phii-1n.gif,eqa_rhoi-12n.gif,eqa_di-12n.gifから,3次補間の係数eqa_cipcsl3_7.gifを算出する.

eqa_cipcsl3_8.gif

ここで,eqa_cipcsl3_9.gifである. また,eqa_phiiR.gifについても,

eqa_cipcsl3_10.gif

となる. また,積分値eqa_rhoi-12n.gifは次式で更新する.

eqa_cipcsl3_11.gif

ここで,eqa_gi.gifeqa_delta_t.gif間に境界eqa_xi.gifを通るeqa_rho.gifの流束(flux)である.

eqa_cipcsl3_12.gif

これらにより算出したeqa_tilde_phi.gifより,n+1ステップでの更新された値は,

eqa_cipcsl3_13.gif

である.

slope eqa_di-12.gifは,Hyman*33,CW*34,UNO*35などの近似手法でeqa_phin.gifより計算する. この3つの近似法の詳細な比較はXiaoらの論文に述べられている.

  • 矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_cip-csl4.jpg
    矩形波のCIP-CSL4法による移流
    rectan_cip-csl2.jpg
    矩形波のCIP-CSL2法による移流

体積保存を確かめるためにSin速度場で移流させた結果を示す.

矩形波(0.25 < x < 0.75で1,それ以外で0)をSin速度場(u=sin(πx/5))でt=5まで移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.青い線は速度場,灰色の細い線は1秒ごとの履歴を示している.

sin_rectan_cip.jpg
矩形波をSin速度場でCIP法により移流
sin_rectan_cip-csl4.jpg
矩形波をSin速度場でCIP-CSL4法により移流
sin_rectan_cip-csl2.jpg
矩形波をSin速度場でCIP-CSL2法により移流

保存系セミラグランジュ法(Conservative Semi-Lagrangian scheme)

保存系セミラグランジュ法*36は,セミラグランジュ法の重みを正規化することで 体積保存を実現した移流法である.

スカラー量 eqa_phi.gif の移流方程式

eqa_csl1.gif

密度 eqa_rho.gif を用いた質量保存式

eqa_csl2.gif

移流方程式×eqa_rho.gif + 質量保存式×eqa_phi.gifで,

eqa_csl3.gif

積の微分の法則( eqa_csl4.gif )より以下の式が導かれる.

eqa_csl5.gif

ここで,eqa_csl6.gif とすると,eqa_hat_phi.gif は保存量として扱える.

グリッド中心座標をeqa_Vxi.gif とすると,セミラグランジュ法など移流法は基本的には以下のように重み付き和で表すことができる.

eqa_csl7.gif

ここで,eqa_wij.gif は重みで,eqa_csl8.gif

本来,完全に質量が保存されるならばどのグリッドにおいても,eqa_csl9.png となるべきであるが, 全グリッドで移流した後に eqa_sigma_j.gif を調べると,eqa_csl10.gifeqa_csl11.gif が起こりうる.これを eqa_csl12.gifと なるように修正する.

  • eqa_csl13.gifの場合
    eqa_csl13.gifの場合は簡単で,単純に重みをeqa_sigma_j.gifで割る.
    eqa_csl14.gif
  • eqa_csl11.gifの場合
    eqa_csl11.gifの場合は,足りない物理量を足す必要がある. そのため,フォワードトレースした結果を追加する. フォワードトレースしたときの重みを eqa_fij.gif とすると正規化した重みは以下となる.
    eqa_csl15.gif

最終的に正規化した重みを用いて値を更新する.

eqa_csl16.gif

実装

  1. バックトレースで eqa_csl17.gif を計算し,重み eqa_wij.gif を変数に格納しておく.
  2. 各グリッドで eqa_csl18.gif を計算
  3. eqa_csl11.gif ならば,eqa_fij.gif をフォワードトレースで求め,正規化した重み eqa_csl19.gif を算出.eqa_csl13.gif ならば,eqa_csl20.gif を計算.
  4. 正規化した重みで eqa_phin+1.gif を算出
  • 矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.灰色の細い線は1秒ごとの履歴を示している.

    rectan_csl.jpg
    矩形波の保存系セミラグランジュ法による移流

体積保存を確かめるためにSin速度場で移流させた結果を示す.

矩形波(0.25 < x < 0.75で1,それ以外で0)をSin速度場(u=sin(πx/5))でt=5まで移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=0.1*Δx/uとした.青い線は速度場,灰色の細い線は1秒ごとの履歴を示している.

sin_rectan_sl.jpg
矩形波をSin速度場でセミラグランジュ法(線形補間)により移流
sin_rectan_csl.jpg
矩形波をSin速度場で保存系セミラグランジュ法により移流

時間微分の離散化

前進オイラー法

空間微分eular.eq1.gifに関する離散化法について多く述べてきた一方,時間微分eular.eq2.gifに関しては以下の単純な1次精度の離散化しか述べてない.

eular.eq3.gif

これは前進オイラー法と呼ばれている.これに対して,レベルセット法などでよく用いられるのがTVDルンゲクッタ(Total-Variation Diminishing Runge-Kutta:以下TVD RK)である.ルンゲクッタ法は時間方向の離散化法としてもっともポピュラーな方法であり,次数を上げていけばどんどん精度も上がる(高次(8次精度など)のルンゲクッタ(RK)は衛星の軌道計算などにも用いられているらしい).実際には前進オイラー法は1次精度のRKと同じである.TVD RKはTVDを満たすRKである.また,計算時間はかかるが安定した数値計算が可能な陰解法や予測子・修正子法などもある.

ルンゲクッタ法(Runge-Kutta scheme)

ルンゲクッタ法(以下RK)はrk.eq1.gifから複数の段階を経てより高精度な近似を行う. 段数を上げることでより高次の精度を得られる. ただし,段数=次数となるのは4段4次までで,それ以降は5段4次や8段6次など段数よりも低い次数になる. そのため4次のルンゲクッタ(RK4)がもっともよく使われる. ちなみに,前進オイラー法は1段1次のルンゲクッタと同じである.

移流方程式の時間微分以外の項をrk.eq2.gifとして,一般的なRKは,

rk.eq3.gif

ここで,rk.eq4.gifである. 上付きのrk.eq5.gifはRKのステップごとの中間値を示している.

RK2(改良オイラー法)

2段2次のルンゲクッタは改良オイラー法(modified Eular method)とも呼ばれている.

rk.eq6.gif

RK3

3段3次のルンゲクッタの式は以下.

rk.eq7.gif

RK4

4段4次のルンゲクッタの式は以下.

rk.eq8.gif

例:矩形波の移流

矩形波(0.25 < x < 0.75で1,それ以外で0)を一定速度(u=0.75)で移流させた結果を示す.グリッド解像度N=128, シミュレーション空間の大きさL=5.0,グリッド幅Δx=L/N,タイムステップ幅Δt=(1.0/CFL)*Δx/uとして,CFL=7, 1.1, 0.535の場合を示す.

rectan_rk_cfl7.jpg
前進オイラー,RK2, RK3での移流結果.CFL=7でΔt=0.00744.
rectan_rk_cfl11.jpg
RK2,RK3,RK4での移流結果.CFL=1.1でΔt=0.0473.
rectan_rk_cfl0535.jpg
RK3,RK4での移流結果.CFL=0.535でΔt=0.0972.

ルンゲクッタの次数が上がるにつれてより大きなタイムステップでも安定して計算できるようになる.下2つのグラフではCFL=1.1, 0.535と中途半端な数字になっているが,これはRK2,RK3が振動を始める限界を探って設定したためである.ちなみにRK4の限界はCFL=0.49ぐらいなので,RK3とRK4の安定性は大きくは変わらない.

TVDルンゲクッタ(Total-Variation Diminishing Runge-Kutta)

ルンゲクッタ法などの解法は1次元,一定幅グリッドではTV安定であるが, 多次元,可変幅グリッドではどうか?という問題がある. これに対して,TVDを満たすルンゲクッタがTVD RK(1)である(TVDについては下参照).

移流方程式の時間微分以外の項をtvdrk.eq1.gifとして,一般的なRKは,

tvdrk.eq2.gif

ここで,tvdrk.eq3.gifである. 上付のtvdrk.eq4.gifはRKのステップごとの中間値を示している.

これに対して,TVD RKは以下である.

tvdrk.eq5.gif

ここで,

tvdrk.eq6.gif

TVD RK2

2次精度のTVD RKはRKにおける修正オイラー法(RK2)に対応する.

tvdrk.eq7.gifで,

tvdrk.eq8.gif

,ら,

tvdrk.eq9.gif

△ら,

tvdrk.eq10.gif

ここで,通常のRKのtvdrk.eq11.gifを用い, さらに,tvdrk.eq12.gifとすると,

tvdrk.eq13.gif

よって,

tvdrk.eq14.gif

TVD RK2ではtvdrk.eq15.giftvdrk.eq16.gifについての2回の前進オイラー法の組み合わせで成り立っている. この事実を使った実装を以下に示す.

  1. 前進オイラー法によりtvdrk.eq16.gifを求める.
    tvdrk.eq17.gif
  1. 前進オイラー法によりtvdrk.eq16.gifからtvdrk.eq18.gifを求める.
    tvdrk.eq19.gif
  1. tvdrk.eq20.giftvdrk.eq18.gifの平均をとることで,最終的な値を得る.
    tvdrk.eq21.gif

TVD RK3

tvdrk.eq22.gifで,

tvdrk.eq23.gif

,ら,

tvdrk.eq24.gif

△ら,

tvdrk.eq25.gif

ここで,通常のRKのtvdrk.eq26.gifを用い,

tvdrk.eq27.gif

tvdrk.eq28.gifとすると,

tvdrk.eq29.gif

よって,

tvdrk.eq30.gif

TVD RK3では3回の前進オイラー法の組み合わせ(tvdrk.eq15.giftvdrk.eq16.giftvdrk.eq31.gifについて)で成り立っている. TVD RK2と同様に,この事実を使った実装を以下に示す.

  1. TVD RK2と同じく2回の前進オイラー法によりtvdrk.eq32.gifを求める.
    tvdrk.eq33.gif
  1. 以下の凸結合によりtvdrk.eq31.gif
    tvdrk.eq34.gif
  1. tvdrk.eq31.gifから前進オイラー法でtvdrk.eq35.gifを算出する.
    tvdrk.eq36.gif
  1. 最終的に,以下の凸結合によりtvdrk.eq37.gifを近似する.
    tvdrk.eq38.gif

TVD RK4

tvdrk.eq39.gifで,

tvdrk.eq40.gif

各係数の導出は長くなりそうなので省略((1)参照}.

tvdrk.eq41.gif

TVD,TVB

TVD(Total-Variation Diminishing)(2)は非線形方程式の収束条件のひとつである. nステップ目でのTV(Total-Variation : 全変動)は以下のように定義される.

tvdrk.eq42.gif

これを離散化すると,

tvdrk.eq43.gif

となる.そして,

tvdrk.eq44.gif

を満たすとき,TV安定であるといい,その数値計算手法はTVDと呼ばれる. また,

tvdrk.eq45.gif

を満たすとき,TVB(Total-Variation Bounded)であるという. ここでtvdrk.eq46.giftvdrk.eq47.gifのみに依存する定数である.

アダムス・バシュフォース法(Adams-Bashforth scheme)

adams-bashforth.eq1.gifの計算に,adams-bashforth.eq2.gifだけでなくさらに前の値(adams-bashforth.eq3.gif)を使うことで高精度にする. 例えば2次精度では以下となる(Adams-Bashforthの2点公式).

adams-bashforth.eq4.gif

完全陰解法(full implicit scheme)

オイラーやルンゲクッタなどでは,求めたいステップでの値(例えばfull-implicit.eq1.gif)のために それ以前の値(full-implicit.eq2.gifなど)だけを空間微分項に用いていた. このため,前のステップでの値が既知であれば,単に式にその値を代入するだけで計算が可能であった. これに対し,full-implicit.eq1.gifを時間微分項以外に用いる方法を陰解法と呼ぶ.

full-implicit.eq3.gif

時間微分以外にfull-implicit.eq2.gifがない場合を完全陰解法と呼ぶ. 移流方程式の完全陰解法による差分式(前進オイラー+風上差分)は,

full-implicit.eq4.gif

一つの式に未知数が2つ(full-implicit.eq5.gif)含まれており,陽解法のように単純に代入だけでは解けない. そのため,すべてのfull-implicit.eq6.gifに関する式を連立させて解く(連立方程式). 流体シミュレーションではほとんどの場合,線形連立方程式(線形システム)となる. つまり,上記の式をすべてのfull-implicit.eq6.gifに関して立てることで,

full-implicit.eq7.gif

という行列方程式を得る.そしてこれを逆行列計算により解く.

full-implicit.eq8.gif

数値流体解析ではこの行列full-implicit.eq9.gifが巨大になってしまい,計算コストが大きくなる. 例えば,3次元空間をfull-implicit.eq10.gifのグリッドで分割した場合,full-implicit.eq11.gifの行列となってしまう. ただし,full-implicit.eq9.gifは正定対称疎行列であることが多く,前処理付共役勾配法など効率的な解法を使うことができる. それでも,やはり陽解法に比べると圧倒的に計算時間が必要である. 一方で陰解法は数値的に安定しており,大きなタイムステップでも発散することなしに解ける(CFL条件にも左右されない)ため,CG用アプリケーションなどではよく用いられている.

クランク・ニコルソン法(Crank-Nicolson scheme)

クランク・ニコルソン法は陰解法の一種であり,時間微分以外の項にcrank-nicolson.eq1.gifcrank-nicolson.eq2.gifの平均値を用いる.

crank-nicolson.eq3.gif

クランク・ニコルソン法は2次精度であり,陰解法であるので行列方程式を解く必要がある.

予測子・修正子法(predictor-corrector scheme)

アダムス・バシュフォース法などの陽解法で予測子を求め,陰的な方法で予測子を修正することでpredictor-corrector.eq1.gifを得る方法. 最終的に陰解法を用いるが,陽解法による予測子を初期値を用いることで,少ない反復で解に達し,計算時間を減らすことができる.


*1 A. Selle, R. Fedkiw, B. Kim, Y. Liu, J. Rossignac, "An unconditionally stable MacCormack method", J. Sci. Comput. 35(2), pp.350-371, 2008.
*2 R. MacCormack, "The effect of viscosity in hypervelocity impact cratering", AIAA Hypervelocity Impact Conference, 1969.
*3 P. Lax and B.Wendroff, "Systems of conservation laws", Commun. Pure Appl Math. 13, pp.217-237, 1960.
*4 B. P. Leonard, "A Stable and Accurate Convective Modelling Procedure Based on Quadratice Upstream Interpolation", Comput. Meths. Appl. Mech. Eng. 19, pp.59-98, 1979.
*5 T. Kawamura and K. Kuwahara, "Computation of High Reynolds Number Flow Around a Circular Cylinder with Surface Roughness", AIAA Paper 84-0340, 1984.
*6 A. Harten, B. Engquist, S. Osher and S. Chakravarthy, "Uniformly high-order accurate essentially non-oscillatory schemes III", J. Comput. Phys. 71, pp.231-303, 1987.
*7 C.-W. Shu and S. Osher, "Efficient implementation of essentially non-oscillatory shock capturing schemes", J. Comput. Phys. 77, pp.439-471, 1988.
*8 C.-W. Shu and S. Osher, "Efficient implementation of essentially non-oscillatory shock capturing schemes II", J. Comput. Phys. 83, pp.32-78, 1989.
*9 S. Osher and J. Sethian, "Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations", J. Comput. Phys. 79, pp.12-49, 1988.
*10 凸結合とは一次結合eqa_weno9.gifでアフィン結合(eqa_weno10.gif),かつ,eqa_weno11.gif
*11 X. D. Liu, S. Osher, and T. Chan, "Weighted essentially non-oscillatory schemes", J. Comput. Phys. 115(1), pp. 200-212, 1994.
*12 G. Jiang and C.-W. Shu, "Efficient implementation of weighted ENO schemes", J. Comput. Phys. 126, pp.202-228, 1996.
*13 G. Jiang and D. Peng, "Weighted ENO Schemes for Hamilton Jacobi Equations", SIAM J. Sci. Comput. 21, pp.2126-2143, 2000.
*14 T. Dupont, Y. Liu, "Back and forth error compensation and correction methods for removing errors induced by uneven gradients of the level set function", J. Comput. Phys. 190(1) (2003) pp.311-324.
*15 B.-M. Kim, Y. Liu, I. Llamas, J. Rossignac, "Using BFECC for fluid simulation", Eurographics Workshop on Natural Phenomena 2005, 2005.
*16 B.-M. Kim, Y. Liu, I. Llamas, J. Rossignac, "Advections with significantly reduced dissipation and diffusion", IEEE Trans. Vis. Comput. Graph. 13(1), pp.135-144, 2007.
*17 H. Takewaki and T. Yabe, "The cubic-interpolated pseudo particle (CIP) method: application to nonlinear and multi-dimensional hyperbolic equations", J. Comput. Phys. 70, pp.355-372, 1987.
*18 T. Yabe and E. Takei, "A New Higher-Order Godunov Method for General Hyerbolic Equations", Journal of the Physical Society of Japan, 57, pp.2598-2601, 1988.
*19 T. Yabe, T. Ishikawa, Y. Kadota and F. Ikeda, "A Multidimensional Cubic-Interpolated Pseudoparticle (CIP) Method without Time Splitting Technique for Hyperbolic Equations", Journal of The Physical Society of Japan, 59, pp.2301-2304, 1990.
*20 提案された当初はCubic Interpolated Pseudo-particle Schemeの略だったが,その後さまざまな改良手法が提案され,3次(Cubic)多項式以外も用いるようになったので,Constrained Interpolation Profile schemeに変更された.
*21 F. Xiao, T. Yabe, and T. Ito, "Constructing oscillation preventing scheme for advection equation by rational function", Comput. Phys. Commun. 93, pp.1-12, 1996.
*22 F. Xiao, T. Yabe, G. Nizam and T. Ito, "Constructing a multi-dimensional oscillation preventing scheme for the advection equation by a rational function", Comput. Phys. Commun. 94, pp.103-118, 1996.
*23 T. Yabe and F. Xiao "Description fo Complex and Sharp Interface during Shock Wave Interaction with Liquid Drop", Journal of the Physical Society of Japan, 62, pp.2537-2540, 1993.
*24 池端昭夫, 肖鋒, "保存型自由表面捕獲スキームと固液気3相流への適用", 日本機械学会2002年度年次大会, 2002
*25 M. Fujisawa and K. T. Miura, "Animation of ice melting phenomenon based on thermodynamics with thermal radiation", Proc. GRAPHITE2007, pp.249-256, 2007.
*26 R. Tanaka, T. Nakamura, T. Yabe, "Constructing exactly conservative scheme in a non-conservative form", Comput. Phys. Commun. 126(3), pp.232-243, 2000.
*27 T. Yabe, R. Tanaka, T. Nakamura, F. Xiao, "An exactly conservative semi-Lagrangian scheme (CIP-CSL) in one dimension", Mon. Weather Rev. 129, pp.332-344, 2001.
*28 T. Nakamura, R. Tanaka, T. Yabe and K. Takizawa, "Exactly Conservative Semi-Lagrangian Scheme for Multi-dimensional Hyperbolic Equations with Directional Splitting Technique", J. Comput. Phys. 174, pp.171-207, 2001.
*29 K. Takizawa, T. Yabe, T. Nakamura, "Multi-dimensional semi-Lagrangian scheme that guarantees exact conservation", Comput. Phys. Commun. 148(2), pp.137-159, 2002.
*30 F. Xiao, T. Yabe, "Completely conservative and oscillationless semi-Lagrangian schemes for advection transportation", J. Comput. Phys. 170(2), pp.498-522, 2001.
*31 F. Xiao, "Unified formulation for compressible and incompressible flows by using multi-integrated moments I: one-dimensional inviscid compressible flow", J. Comput. Phys. 195, pp.629-654, 2004
*32 F. Xiao, R. Akoh and S. Ii, "Unified formulation for compressible and incompressible flows by using multi-integrated moments II: Multi-dimensional version for compressible and incompressible flows", J. Comput. Phys. 213, pp.31-56, 2006.
*33 J. H. Hyman, "Accurate monotonicity preserving cubic interpolations", SIAM Journal on Scientific Computing, 4, pp.645-654, 1983.
*34 A. Harten and S. Osher, "Uniformly high order accurate non-oscillatory schemes. I", SIAM Journal on Numerical Analysis, 24, pp.279-309, 1987.
*35 P. Collela and P. R. Woodward, "The piecewise parabolic method (PPM) for gas-dynamical simulations", J. Comput. Phys. 54, pp.174-201, 1984
*36 M. Lentine, J. T. Gretarsson and R. Fedkiw, "An unconditionally stable fully conservative semi-Lagrangian method", J. Comput. Phys. 230(8), pp.2857-2879, 2011
(1)  C.-W. Shu and S. Osher, "Efficient implementation of essentially non-oscillatory shock capturing schemes", J. Comput. Phys. 77, pp.439-471, 1988.
(2)  A. Harten, "High resolution schemes for hyperbolic conservation laws", J. Comput. Phys. 49, pp.357-393, 1983.

添付ファイル: fileeqa_phii.gif 1008件 [詳細] fileeqa_advec_equation1.gif 1242件 [詳細]

トップ   編集 凍結 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2011-11-01 (火) 17:59:47 (3199d)