DA 法

体積率移流方程式を保存形に変形した上で離散化すると次式となる．

αⁿ⁺¹_i,j,k =αⁿ_i,j,k+ 1 Ω_i,j,k

h(αu₁)_i+1/2,j,k−(αu₁)_i−1/2,j,ki

∆t∆x₂∆x₃ + 1

Ω_i,j,k

h(αu2)_i,j+1/2,k−(αu2)_i,j−1/2,k

i∆t∆x1∆x3

+ 1 Ωi,j,k

h(αu3)_i,j,k+1/2−(αu3)_i,j,k−1/2

i∆t∆x1∆x2

(1.4.1)

ここで，Ωはセル体積，u1,u2,u3は各々 x1,x2,x3方向の速度成分である．下付添字i, j,k は x₁,x₂,x₃ 方向のセル番号で，変数配置はスタガード配置とする．∆x₁,∆x₂,∆x₃ はセル 幅である．

簡単のため，セル表面(i+1/2, j,k)を通じたx1 方向の移流のみ考える．DA法では，被

輸送体積∆Ω_i+1/2,j,k を次式で与える．

∆Ωi+1/2,j,k =(αu1)_i+1/2,j,k∆t∆x2∆x3

=min

hαIAD|u1∆t|∆x2∆x3+ ΩS U P, αIDΩi,j,k

i (1.4.2)

ここで，

Ω_{S U P} = maxh

(α_{U P}−α_IAD)|u₁∆t|∆x₂∆x₃−(α_{U P}−α_ID)Ω_i,j,k,0i

(1.4.3)

αU P= max (αID, αIDU) (1.4.4)

u∆t

ID

IDU IAD IA

u fluid1

fluid2

x₁ x₂

図1.4.1 DA法におけるセルの記法

下付添字IDは流体の輸送元のセル，IDU は輸送元のさらに上流側のセル，IADは輸送が 行われるセル表面を意味する．これらの記法を図1.4.1に示した．IAは流体を受け取るセ ルを意味する．αIDA は，界面の傾きと隣接セルの体積率の値によって決める．

IAD=

( ID if |n1|<|n2|or |n1|< |n3|

IA otherwise (1.4.5)

ここで，n1,n2,n3 は界面法線nのx1,x2,x3 方向の成分である．上式による判定の後，次 式で再度判定を行う．

IAD= IA (if αIA =0 or αIDU =0) (1.4.6)

また，セル IDが界面セルでない場合は IDA= IAとする．これらのセル設定を式(1.4.3) に適用して被輸送体積を求める．

式(1.4.2)は，IDA= IDの場合風上差分，IDA= IAの場合風下差分に相当する．１次

風上差分は安定性が高く，また解の単調性を維持できる．しかし，数値拡散が生じる．一 方，風下差分は不安定であるが，数値拡散の係数が負であるため界面をシャープに維持し やすい．DA法は，状況に応じてこれら風上差分と風下差分を切り替えて安定な計算を行 う差分法であるが，スキーム切り替えの判定を界面形状に基づいて行うのが特徴である．

各軸方向に分割した体積率移流方程式を順次解くという方法でDA法を多次元化した方 法は，オリジナルのDA法と区別して改良 DA法[85]と呼ばれる．しかし，現在では方 向分割による多次元化は他の方法でもよく行われている．以下ではオリジナルのDA法を DA-a，改良DA法をDA-bと呼ぶ．

1.4.2 FCT

Rudman [86]は，FCT(flux-corrected transport) [87, 88]を体積率移流方程式に応用した．

FCTの概略を以下に述べる．簡単のためx₁ 方向の移流のみ考える．まず，数値拡散を生

じるが，安定かつ解の単調性を維持する低次の解法により移流項を評価し，αの推定値を 得る．

α^∗_i = αⁿ_i + F_i+1/2^L −F_i−1/2^L

(∆x₁)_i (1.4.7)

ここで，F^L は低次の解法で評価した流束である．単調性を維持する解法により F^L を計 算しているため，α^∗は0 < α^∗ <1の範囲に収まる．逆拡散を生じる不安定な高次の解法 により評価した流束をF^H とする．これらの流束から非拡散流束F^A を定義する．

F^A= F^H −F^L (1.4.8)

次式によりαを更新する．

αⁿ⁺¹_i =α^∗_i + Ci+1/2F_i+1/2^A −Ci−1/2F_i−1/2^A

(∆x₁)_i (1.4.9)

ここで，C (0<C < 1)は制限関数である．解の単調性を維持するように制限関数を決定 するため，α^∗の最大値α^maxおよび最小値α^minを調べる．

α^max_i =max

αⁿ_i−1, αⁿ_i, αⁿ_i+1, α^∗_i−1, α^∗_i, α^∗_i+1

(1.4.10) ここで，以下の量を定義する．

α^min_i = min

αⁿ_i−1, αⁿ_i, αⁿ_i+1, α^∗_i−1, α^∗_i, α^∗_i+1

(1.4.11) Q_i+1/2^I = α^max_i −α^∗_i (∆x₁)_i

∆t (1.4.12)

Q_i+1/2^D =

α^∗_i −α^min_i (∆x1)_i

∆t (1.4.13)

P_i+1/2^I =max

0,F_i−1/2^A

−min

0,F^A_i+1/2

(1.4.14) P_i+1/2^D =max

0,F_i+1/2^A

−min

0,F^A_i−1/2

(1.4.15) R^I_i+1/2 = min

0,Q^I_i+1/2/P^I_i+1/2

(1.4.16) R^D_i+1/2 = min

0,Q^D_i+1/2/P^D_i+1/2

(1.4.17) 制限関数Cを次式で与える．

Ci+1/2 = min

R^I_i+1/2,R^D_i+1/2

(1.4.18)

Rudman は，体積率移流方程式を各軸方向に分割し，上述の１次元 FCT を各方向

に適用するという方法で多次元化している．また，FCT ではなく TVD(total variation diminishing)を用いる方法もある．

i j

i+1

x₁ x₂

case1

case4

case7

case2

case5

case8

case3

case6

case9

図1.4.2 界面形状の分類

1.4.3 PLIC 型