![[PukiWiki]](image/pukiwiki.png "[PukiWiki]")
|
Navier-Stokes方程式での移流項 (u・▽)u について
移流方程式†移流方程式  空間微分の離散化†風上差分(Upwind differencing)†移流方程式の時間微分を前進オイラー法で離散化すると，  1次元の場合を考える．グリッド  ここで， 風上差分ではその名の通り，風上側の値を使って差分を行う1次精度の離散化法である．
1次元の場合，  風上差分はCFL(Courant-Friedreichs-Lewy)条件を満たしていれば安定である． CFL条件によるタイムステップ幅の制限は，  タイムステップ幅は通常これよりもさらに小さい数を用いた方がよい．
係数:CFL数(CFL number)を  ここで， 振動と数値拡散†風上差分の離散式( 
さて，上式の微分を下で述べる中心差分で離散化してみる．  風上差分は中心差分に拡散項を加えたものであることが分かる． 中心差分では非常に大きな振動が発生する． その振動を拡散項により抑えているのが風上差分である． この振動を抑えるという考え方は重要である． より高次の関数を使って補間することで振動を抑えるのが，Lax-WendroffやQUICK,QUICKEST,河村・桑原スキームなどで， 2階，3階の微分を使って抑えるのがENOやWENOである． これらに対して，元の形状を保持するという考え方に基づき，異なるアプローチをとるのがCIP法,RCIP法などである． ちなみにRCIP法では数値拡散が発生してしまうので，これを抑えるためのSTAA手法では逆に数値拡散項を引くことで拡散を抑えている．
中心差分(Central differencing)†流れの速度  前進オイラー法と組み合わせると， 
Lax-Wendroff法†Lax-Wendroff法(または，ライス(Leith)法)*3では線形補間の代わりに2次多項式を使った補間で近似を行う．2次多項式の係数を求めるために，グリッドiとその両隣(i-1,i+1)の3つの値 2次多項式は，  ここで各係数は， 
 風上差分と同様に中心差分に拡散項を加えたものになっている．拡散項の係数が異なるのが風上差分との違いである．
QUICK(Quadratic Upstream Interpolation for Convective Kinematics)†QUICK*4は  これを移流方程式の離散化に用いると， 
QUICKEST(QUICK with estimated streaming terms)†QUICKEST(QUICK with estimated streaming terms)はLax-Wendroffの3つ( 
河村・桑原スキーム(Kawamura-Kuwabara scheme)†4次精度の中心差分( 4次精度の中心差分は以下．  4階微分項の差分式は，  (4階微分の差分式の導出は4階微分を参照)． 4階微分項に係数 
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の式を述べる前に，準備として差分式を定義しておく．
 となる．ここで，iはグリッド番号(座標値 1階差分はグリッド間(i-1/2とi+1/2)で定義される．   2階差分はi-1/2とi+1/2での1階差分値を使って，iで定義される．  同様に3階差分は，   ENOでは3次多項式により  ここで， 
これらによって，  により
WENO(Weighted Essentially Non-Oscillatory polynomial interpolation)†ENOを改良して，重み付き平均をとることで5次精度を実現したのがWENOである． 3次精度のENOでは  1+2+4, 1+2+5, 1+3+6, 1+3+7の4パターンになるが， 以下に示すように1+2+5と1+3+6が同じなので，組み合わせは3パターンになる． それぞれさらに展開してφiの式にする．  差分を  ENO近似は結局この3つのパターンに集約される．
これらの内どれかが状況に応じて用いられる．
WENO(Weighted ENO)はその名の通り，これら3つの凸結合
*10
で  ここで， 問題となるのは
 先ほどの重みの組み合わせ(  ここで， 
 これらから，  により
セミラグランジュ(semi-Lagrangian)法†CG分野でももっともポピュラーな方法．下図のように計算点xから-u方向にΔtだけバックトレースして，その位置での速度場u(x_bt, t)を次のステップの速度場u(x, t+Δt)とする方法．下図ではバックトレースに1次関数を用いている． この方法は大きなタイムステップ幅でも安定して計算できるのが大きな特徴である．  セミラグランジュ法ではバックトレースした位置での値を周囲の値を使って補間しなければならない． 一番簡単なのは線形補間を用いる方法である．ただし線形補間は1次精度であり，上記の風上差分と同じぐらい拡散が発生する． より高精度な補間法を用いることで拡散を抑えることができる．例えば，上記のENOやWENO，下記で示す方法などである．
BFECC(Back and Forth Error Compensation and Correction)†BFECCはSemi-Lagrangian法においてバックトレースする際に，フォワードトレースにより誤差を求め， その誤差を修正することで高精度に移流を行う方法である． BFECCの最初のアイデアは *14 で提案され， 2005年にCG分野で移流方程式へ適用された *15, *16． Semi-Lagrangian法のバックトレースによる  Lは直線を使った1次精度のバックトレース，uは速度場を表す．  ここで， 
BFECCの式を上の図を使って説明する．  となる．上式から，誤差eは以下のように求められる．  より正確な 
MacCormack†CIP法(Constrained Interpolation Profile scheme)†高次多項式を用いた補間法では補間式を得るためにより多くのグリッドを必要とする．例えば，ENO,WENOだと自身も含めて6グリッド(ENOで実際に用いられるのは4グリッド)の値が必要．風上差分だと2グリッドであったことを考えると3倍になっている．これは当然のことで，高次多項式の係数を得るためには多くの関数値を必要とするからである．このように用いるグリッドが増えると，特に境界付近での処理が難しくなってくる．これに対して，風上差分と同様に2グリッドのみを用いつつ，勾配値を利用することで急激に変化する関数でも対応したのがCIP法(Constrained Interpolation Profile scheme)*17*18*19である*20． CIP法の基本的なアイデアは，移流式をxで微分したとき，その導関数  グリッド上の4拘束条件  この手法は非常に簡単に見えるが，その効果は非常に良好である．特に矩形波のように急激に変化する特徴を持つ関数を移流させた場合でも，その形状を維持したまま移流させることができる．
RCIP法(Rational Constrained Interpolation Profile scheme)†有理関数を使ったCIP法． CIP法の欠点は関数を移流させるとき，若干のオーバーシュートがでることである． このオーバーシュートは移流が進んでも大きくはならないため問題がないように思えるが， 水と空気の二相流など，計算したい物理量の差が非常に大きい場合に， 小さなオーバーシュートが大きな問題となる． RCIP法 *21, *22 はこの問題を解決する． RCIP法はCIP法の3次関数の代わりに有理関数  を用いる．
ここで，  ここで，  RCIP法以外にも，Tangent変換を用いることで オーバーシュートを抑える方法も提案されている *23． Tangent変換を用いた手法では物体界面がシャープに保たれるという特徴を持つ． これは数値流体計算上は非常に有益であるが， グラフィックスで用いた場合，グリッド解像度によっては 液体表面があまり滑らかにならず，保存性も保証されない． RCIP法による密度関数の移流では，界面を表す0から1への変化部分が拡散し， 関数の傾きがなだらかになる． この滑らかな界面はCGでの利用においては美しいレンダリング結果を生むが， シミュレーションでは数値的な不安定性を引き起こす． 密度関数の体積移動によるSTAA手法 *24 を使うことで， 界面拡散を抑制することができる． また，拡散をある程度の幅で制御する方法 *25 もある．
CIP-CSL法†CIP法では格子iのプロファイル CIP-CSL4†CIP-CSL4法*26，では，CIP法の4拘束条件に関数  つまり，  の5拘束条件を満足させるために以下の4次多項式をプロファイルに用いるのが CIP-CSL4法である．  ここで， 5拘束条件  CIP法では算出した係数
 ここで，  である． CIP-CSL2†CIP-CSL2法*27,*28,*29は，
値  $D$の微分  これは1次元の移流方程式と同じであり，
CIP法での 格子点i内での 
 CIP法における空間微分値gは 
 また，  これらの条件により係数 
  CIP-CSL3†CIP-CSL3法*30,*31,*32は，グリッドセル[ CIP-CSL3法による移流項の解法では，まず，nステップでの値  である． 移流速度場の方向により以下の補間を使い分ける．  そして，  左側要素  
 ここで，  となる．
また，積分値  ここで，  これらにより算出した  である． slope
体積保存を確かめるために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秒ごとの履歴を示している．
保存系セミラグランジュ法(Conservative Semi-Lagrangian scheme)†保存系セミラグランジュ法*36は，セミラグランジュ法の重みを正規化することで 体積保存を実現した移流法である． スカラー量  密度  移流方程式×  積の微分の法則(  ここで， グリッド中心座標を  ここで， 本来，完全に質量が保存されるならばどのグリッドにおいても，
最終的に正規化した重みを用いて値を更新する．  実装†
体積保存を確かめるためにSin速度場で移流させた結果を示す． #include(): Included already: Sin速度場 - 設定
時間微分の離散化†前進オイラー法†空間微分  これは前進オイラー法と呼ばれている．これに対して，レベルセット法などでよく用いられるのがTVDルンゲクッタ(Total-Variation Diminishing Runge-Kutta:以下TVD RK)である．ルンゲクッタ法は時間方向の離散化法としてもっともポピュラーな方法であり，次数を上げていけばどんどん精度も上がる(高次(8次精度など)のルンゲクッタ(RK)は衛星の軌道計算などにも用いられているらしい)．実際には前進オイラー法は1次精度のRKと同じである．TVD RKはTVDを満たすRKである．また，計算時間はかかるが安定した数値計算が可能な陰解法や予測子・修正子法などもある． ルンゲクッタ法(Runge-Kutta scheme)†ルンゲクッタ法(以下RK)は 移流方程式の時間微分以外の項を  ここで， RK2(改良オイラー法)†2段2次のルンゲクッタは改良オイラー法(modified Eular method)とも呼ばれている．  RK3†3段3次のルンゲクッタの式は以下．  RK4†4段4次のルンゲクッタの式は以下．  例:矩形波の移流†矩形波(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の場合を示す．
ルンゲクッタの次数が上がるにつれてより大きなタイムステップでも安定して計算できるようになる．下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¬e{Shu1988:C.-W. Shu and S. Osher, "Efficient implementation of essentially non-oscillatory shock capturing schemes", J. Comput. Phys. 77, pp.439-471, 1988.};である(TVDについては下参照)． 移流方程式の時間微分以外の項を  ここで， これに対して，TVD RKは以下である．  ここで，  TVD RK2†2次精度のTVD RKはRKにおける修正オイラー法(RK2)に対応する．
 ?から，  ?から，  ここで，通常のRKの  よって，  TVD RK2では
TVD RK3†
 ?から，  ?から，  ここで，通常のRKの 
 よって，  TVD RK3では3回の前進オイラー法の組み合わせ(
TVD RK4†
 各係数の導出は長くなりそうなので省略(¬e{Shu1988};参照}．  TVD,TVB†TVD(Total-Variation Diminishing)¬e{Harten1983:A. Harten, "High resolution schemes for hyperbolic conservation laws", J. Comput. Phys. 49, pp.357-393, 1983.};は非線形方程式の収束条件のひとつである． nステップ目でのTV(Total-Variation : 全変動)は以下のように定義される．  これを離散化すると，  となる．そして，  を満たすとき，TV安定であるといい，その数値計算手法はTVDと呼ばれる． また，  を満たすとき，TVB(Total-Variation Bounded)であるという．
ここで アダムス・バシュフォース法(Adams-Bashforth scheme)†
 完全陰解法(full implicit scheme)†オイラーやルンゲクッタなどでは，求めたいステップでの値(例えば  時間微分以外に  一つの式に未知数が2つ(  という行列方程式を得る．そしてこれを逆行列計算により解く．  数値流体解析ではこの行列 クランク・ニコルソン法(Crank-Nicolson scheme)†クランク・ニコルソン法は陰解法の一種であり，時間微分以外の項に  クランク・ニコルソン法は2次精度であり，陰解法であるので行列方程式を解く必要がある． 予測子・修正子法(predictor-corrector scheme)†アダムス・バシュフォース法などの陽解法で予測子を求め，陰的な方法で予測子を修正することで |
でアフィン結合(
)，かつ，
(座標
)について，
．この
を差分を使って求める．
では座標軸上の右から左に風が吹いているので，風上はi+1，
では左から右に風が吹いているので，風上はi-1となる．
とすると，
．よく使われるのは
である．
は空間2階微分であり，これは拡散を表している．
拡散に関する項は，例えば，流体のナビエ・ストークス方程式にもある(
)．
しかし，これはそれらの項とは異なり，離散化によって生まれた拡散であり，これを数値拡散(numerical diffusion)と呼ぶ．
下の矩形波の移流を見ると風上差分において数値拡散が発生しているのがよく分かる．
に関係なく，グリッド点を中心とした周囲のグリッドの値の平均で空間微分を近似するのが中心差分である．
つまり，1次元の場合，以下となる．
を用いる．
とすると移流方程式の近似更新式が得られる．
の内挿に
と
に加えて，風上側のもう一つの
を用いる方法である．
,
,
)に加えて，QUICKと同様に風上側のもう一つの
を追加して近似する手法である(
ならi-2，
ならi+2)．
を使用)に
拡散項として4階微分項を追加したのが，河村・桑原スキーム
を掛けた上で先ほどの式に追加する．
の場合，
のn階差分を
と表記する．とすると，n=0は，
)．
を近似する．
はm次項である．
, 
の式がほしいので，
上式を微分して，x=
の場合k=i-1，
でk=iとすると，
を求め，
について，
の6つの内，4つの値を用いて補間を行う．
そして，その組み合わせは3パターン
(
と
と
)である．
実際に展開してみと以下となる(
では，k=i-1)．
と表記すると上記3パターンは，
を近似する
．
の選択である．
関数が滑らかならば，
で5次精度が得られる
)をsmoothnessの2乗で割ることで重みを得る．
はゼロ割を防ぐための項である．
は
が符号付距離場であることが前提で設定されている．
つまり，
の場合である．
最終的に
を満たすために，正規化を行う．
の更新を以下の式で表す．
)．
このとき，B点とC点の間の
とすると，
は，
からeを引いた場
より得られる．
をフォワードトレース
を計算
もまた，
と同様に移流するということである(ただし，
とする)．2メッシュ[i, i+1]間のプロファイルは，以下の3次補間関数で表される．
から，
はパラメータで0か1の値をとる．
のときはCIP法と同じであり，
のとき，
この関数は単調性と凹凸性を保存する有理関数となる．
であり, 
である．
上式の各係数は以下である．
の積分値
をさらに加える．
,
,
,
,
はグリッドi内での多項式補間パラメータである．
,
,
,
,
と
を更新していたが，
CIP-CSL4法ではさらに積分値
の更新が必要となる．
とその積分値
を用いる．
積分値
に対する移流方程式は，
と定義すると
と置き換えることで
CIP法と同様に計算ができる．
を以下のように定義する．
から補間点までの積分値を表し，
3次関数で補間する．
で置き換えられるので，
[
]のプロファイルは，
,
は以下のように求められる．
の時間発展は，
に挟まれた部分の積分値となる(下図)．
, 
]間の積分値とグリッド中心
における勾配(グリッド端
と区別するためにslopeと呼び 
で表す)を用いてCIP補間を行う手法である．
から，
中間値
を求める．セミラグランジュ法による中間値の算出式(1次元)は，
では，グリッドの端点
である．
その他の点での補間にはグリッド間の値を積分した積分値
とグリッド中心の傾き(slope) 
を用いる．
,
,
を算出する．
である．
また，
についても，
間に境界
の流束(flux)である．
は，Hyman
の移流方程式
を用いた質量保存式
)より以下の式が導かれる．
とすると，
は保存量として扱える．
とすると，セミラグランジュ法など移流法は基本的には以下のように重み付き和で表すことができる．
は重みで，
．
となるべきであるが，
全グリッドで移流した後に 
を調べると，
や 
が起こりうる．これを 
と なるように修正する．
の場合
とすると正規化した重みは以下となる．
を計算し，重み 
を計算
を算出．
を計算．
を算出
に関する離散化法について多く述べてきた一方，時間微分
に関しては以下の単純な1次精度の離散化しか述べてない．
から複数の段階を経てより高精度な近似を行う．
段数を上げることでより高次の精度を得られる．
ただし，段数=次数となるのは4段4次までで，それ以降は5段4次や8段6次など段数よりも低い次数になる．
そのため4次のルンゲクッタ(RK4)がもっともよく使われる．
ちなみに，前進オイラー法は1段1次のルンゲクッタと同じである．
として，一般的なRKは，
である．
上付きの
はRKのステップごとの中間値を示している．
として，一般的なRKは，
である．
上付の
はRKのステップごとの中間値を示している．
で，
を用い，
さらに，
とすると，
と
についての2回の前進オイラー法の組み合わせで成り立っている．
この事実を使った実装を以下に示す．
を求める．
と
で，
を用い，
とすると，
について)で成り立っている．
TVD RK2と同様に，この事実を使った実装を以下に示す．
を求める．
を算出する．
を近似する．
で，
は
のみに依存する定数である．
の計算に，
だけでなくさらに前の値(
)を使うことで高精度にする．
例えば2次精度では以下となる(Adams-Bashforthの2点公式)．
)のために
それ以前の値(
など)だけを空間微分項に用いていた．
このため，前のステップでの値が既知であれば，単に式にその値を代入するだけで計算が可能であった．
これに対し，
)含まれており，陽解法のように単純に代入だけでは解けない．
そのため，すべての
に関する式を連立させて解く(連立方程式)．
流体シミュレーションではほとんどの場合，線形連立方程式(線形システム)となる．
つまり，上記の式をすべての
が巨大になってしまい，計算コストが大きくなる．
例えば，3次元空間を
のグリッドで分割した場合，
の行列となってしまう．
ただし，
と
の平均値を用いる．
を得る方法．
最終的に陰解法を用いるが，陽解法による予測子を初期値を用いることで，少ない反復で解に達し，計算時間を減らすことができる．
eqa_phii.gif
eqa_advec_equation1.gif
