サブセルに基づく移流方程式の解法

式(2.4.3)左辺第２項の評価を注意深く行わねば，数値誤差が生じてしまう．例えば，

数値粘性を含む風上差分を適用した場合，体積率が時間経過とともに拡散する，すなわち 界面のぼやけが生じる．また，αの上下限値(0< α <1)を考慮していないため体積誤差 を生じる．数値拡散や体積誤差が生じると物理的に正しい解が得られないだけでなく，計 算を破綻させてしまう場合もある．本研究では，非一様なサブセルを一時的に界面セル内 に導入して流体体積移流量を計算する手法を開発する．本手法は，サブセルを導入するこ とで複雑な幾何計算を避けるため，容易に３次元化できる．また，詳細は2.5節で述べる が，本手法の界面再構築過程で得られる距離関数を利用することにより，高精度表面張力 評価が可能となる．

式(2.4.3)をセル(i, j,k)で離散化すると次式を得る．

αⁿ⁺¹_i,j,kΩ_i,j,k= αⁿ_i,j,kΩ_i,j,k−





(αu1∆t∆x2∆x3)_i+1/2,j,k−(αu1∆t∆x2∆x3)_i−1/2,j,k +(αu2∆t∆x1∆x3)_i,j+1/2,k−(αu2∆t∆x1∆x3)_i,j−1/2,k +(αu3∆t∆x1∆x2)_i,j,k+1/2−(αu3∆t∆x1∆x2)_i,j,k−1/2



 (2.4.4) ここで，Ωi,j,k =(∆x1)i(∆x2)j(∆x3)kである．以下，簡単のためx1 方向の移流についてのみ 説明する．この場合，上式は，

αⁿ⁺¹_i,j,kΩi,j,k= αⁿ_i,j,kΩi,j,k−h

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

(2.4.5) となる．上式右辺第２項大括弧内の２つの項は，各々セル表面(i+1/2, j,k), (i−1/2, j,k) を通じて隣接セルから流れ込むあるいは隣接セルに流れ出す流体体積量を表わしてい る．以下，セル (i, j,k) から隣接セルに流体が流れ出す条件，すなわち，(u₁)_i+1/2,j,k >

0, (u1)i−1/2,j,k < 0 の場合について説明する．この条件を図 2.4.1 に示す．図中，体積

(u1)i+1/2,j,k∆t∆x2∆x3, (u1)i−1/2,j,k∆t∆x2∆x3 に含まれる流体が隣接セルに移る．これら 被輸送体積中に含まれるαを，界面セル内に一時的に配置するサブセルを利用して算出す る．セル(i, j,k)内で，界面を平面と仮定する．平面の方程式は次式で与えられる．

n·xint +φint = 0 (2.4.6)

ここで，φ_int はセル中心‐界面間の符号付距離，x_int はセル中心を始点，界面上の点を終 点とする距離ベクトルである．nは界面法線であり，次式で計算する．

n= ∇φ

|∇φ| or ∇α

|∇α| (2.4.7)

φは局所レベルセット関数である．φについては2.5節で詳述する．Tomiyamaらによれ ば，法線の計算にφ,αのいずれを用いても，法線方向の算出精度に大きな差はない[11]．

|(u₁)i+1/2, j, k∆t|

(u₁)i+1/2, j, k

(u₁)i−1/2, j, k

|(u₁)i−1/2, j, k∆t|

∆x₁

∆x₂

∆x₃

x₁ x₂ x₃

interface i, j, k

|(u₁)i+1/2, j, k∆t|∆x₂∆x₃

|(u₁)_{i−1/2, j, k}∆t|∆x₂ ∆x₃

n

x_int

図2.4.1 界面セル

|(u₁)i+1/2, j, k∆t|

|(u₁)I−1/2, j, k∆t|

x₁ x₂ x₃

∆x3sub

∆x₂^sub non-transferred

volume N_C

transferred volumes N_1m, N_1p

∆x₁^sub

図2.4.2 サブセル配置

サ ブ セ ル を 利 用 し て φint を 求 め る ．図 2.4.2 に 示 す よ う に ，被 輸 送 体 積 (u₁)i+1/2,j,k∆t∆x2∆x3, (u₁)i−1/2,j,k∆t∆x2∆x3 お よ び 残 留 体 積 に 各 々 N1p, N1m, NC 個 のサブセルを準備する．図中∆x₁^sub,∆x₂^sub,∆x₃^sub はx₁,x₂,x₃ 方向のサブセル幅である．す べてのサブセル体積Ω^sub = ∆x₁^sub∆x^sub₂ ∆x₃^subを加えるとセル体積Ωi,j,kに一致する．

NC

X

l=1

Ω_l^sub+

N_1p

X

m=1

Ω^sub_m +

N1m

X

n=1

Ω_n^sub = Ωi,j,k (2.4.8)

サブセルに含まれる流体量のセル体積Ω_i,j,kに対する割合をサブセル体積率α^sub と定義す

， ， ， ， -~ーーーーーーーーーー

φintn n

d₍₇₎

condition for d_(v) d_(v)^. n < d_(v+1)^. n cell center

subcell

interface

d₍₈₎

d₍₄₎ d₍₃₎

d₍₁₎ d₍₅₎

d₍₆₎

d₍₂₎

図2.4.3 距離ベクトルd(v)

る．定義より，α^sub の値がとりうる範囲は0< α^sub <Ω^sub/Ωi,j,kである．α^sub はφintの関 数である．α^sub(φ_int)を次の区分一次関数で与える．

α^sub_m (φint)=max









min



 Ω^sub_m

4φintn−P₈

v=5d(v)

·n

Ωi,j,k

P₄

v=1d(v)−P₈

v=5d(v)

·n , Ω_m^sub

Ω_i,j,k



,0









 (2.4.9)

なお，上式の代わりに他の関数形を用いることも可能である．d_(v) はセル中心からサブセ ルmの頂点vへの距離ベクトルである(図2.4.3)．次式を満足するように頂点番号を付す．

d_(v)·n≤ d_(v+1)· n (2.4.10)

図2.4.4に示す，界面をその１つの面とし，かつ法線nが内向きである多面体の体積は，

αⁿ_i,j,kΩi,j,kに一致しなければならない．したがって，φ_int は次式を満足しなければならない．

f (φint)=

NC

X

l=1

α_l^sub(φint)+

N_1p

X

m=1

α_m^sub(φint)+

N1m

X

n=1

α^sub_n (φint)−αⁿ_i,j,k = 0 (2.4.11) 上式を反復法により解いてφ_int を得る．体積率α, α^sub の定義より，f (φ_int)のとりうる値

x₁ x₂ x₃

interface

i, j, k n

x_int

n

α_i,j,kΩ_i,j,k

(1−α_i,j,k)Ω_i,j,k

図2.4.4 界面がセルから切り取る体積

transferred fluid volume

|(αu)_i−1/2,j,k∆t∆x₂∆x₃|=

(Σ^N1mn=1α^subn)Ωi,j,k

x₁ x₂ x₃

transferred fluid volume

|(αu)_i+1/2,j,k∆t∆x₂∆x₃|=

(Σ^N1pm=1α^subm)Ωi,j,k

n

φ_int

図2.4.5 再構築された界面および被輸送体積

は−αⁿ_i,j,k < f (φint)< 1−αⁿ_i,j,k の範囲にある．f (φint) = −αⁿ_i,j,k, 1−αⁿ_i,j,k に対応するφint は nに応じて定まり，各々φint の最小値・最大値となる．f (φint)の具体的な関数形は未知で あるため，ニュートン法のように f (φ_int)の導関数を必要とする反復解法は適用できない．

一方，ブレント法[28, 29]では導関数は不要である．また，ブレント法では解の探索範囲 [φmin, φmax]を指定する必要があるが，上述の通りφint の最大値・最小値は既知である．以 上より，ブレント法は式(2.4.11)の解法として適当であると考えられる．本研究ではブレ ント法を用いて上式を解く．ブレント法により方程式を解く手順については付録を参照さ れたい．

式(2.4.11)の解 φint により，式(2.4.5)右辺大括弧内の２項は各々以下の諸式で与えら

れる(図2.4.5)．

(αu₁)_i+1/2,j,k∆t∆x₂∆x₃ =

N1p

X

m=1

α^sub_m (φ_int)Ω_i,j,k (2.4.12)

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

N_1m

X

n=1

α_n^sub(φint)Ωi,j,k (2.4.13) これらの値を式(2.4.5)に代入して，αⁿ_i,j,kをαⁿ⁺¹_i,j,kに更新する．

以上に示した，サブセルに基づく界面再構築・移流手法を NSS(non-uniform subcell scheme)と称する．

2.4.3 体積率移流方程式の時間積分

３次元の場合，単純な計算法で，かつ体積誤差を小さく抑えられる方向分割法を用い る[6]．２次元計算では，時間積分の計算負荷が多少増しても，全体の計算負荷に対する 影響は小さい．そこで，２次元の場合は，方向分割法に比べてやや複雑であるが，非常 に良好な体積保存を実現できるEI-LE(Eulerian implicit - Lagrangian explicit) [30]を採用 する．

３次元方向分割法を以下に説明する．体積率移流方程式，

∂α

∂t +u· ∇α= 0 を保存形に変形すると，

∂α

∂t +∇ ·αu= α∇ ·u (2.4.14)

となる．連続の式(∇ ·u = 0)によれば上式右辺は消えるが，消さずに計算する．上式右 辺は発散修正項 [6]と呼ばれており，以下に述べるように非常に重要な役割がある．式

(2.4.14)を各軸方向に分割する．

α^∗−αⁿ

∆t + ∂αⁿu₁

∂x1

=αⁿ∂u₁

∂x1

(2.4.15) α^∗∗−α^∗

∆t + ∂α^∗u2

∂x₂ =α^∗∂u2

∂x₂ (2.4.16)

αⁿ⁺¹−α^∗∗

∆t + ∂α^∗∗u3

∂x3

= α^∗∗∂u3

∂x3

(2.4.17) これら３式を順次計算して１時間ステップの移流が完了する．このとき，速度場uが連続 の式，∇ ·u = ∂u1/∂x1 +∂u2/∂x2+∂u3/∂x3 = 0を満足していても，一般に∂u1/∂x1 , 0,

∂u₂/∂x₂ , 0, ∂u₃/∂x₃ , 0であるから，１次元に分割した式(2.4.15)の左辺のみを解いた 段階で体積率のオーバーシュートやアンダーシュートが発生する．そこで，発散修正項に より，0< α < 1の範囲に納まるように修正がなされる．式(2.4.15)-(2.4.17)を離散化し

た式を以下に示す．

α^∗_i,j,kΩi,j,k=αⁿ_i,j,kn

Ωi,j,k+h

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

∆t∆x2∆x3

o

−h

(αⁿu1∆t)_i+1/2,j,k−(αⁿu1∆t)_i−1/2,j,ki

∆x2∆x3

(2.4.18) α^∗∗_i,j,kΩi,j,k=α^∗_i,j,kn

Ωi,j,k+h

(u₂)_i,j+1/2,k−(u₂)_i,j−1/2,ki

∆t∆x1∆x3

o

−h

(α^∗u2∆t)_i,j+1/2,k−(α^∗u2∆t)_i,j−1/2,ki

∆x1∆x3

(2.4.19) αⁿ⁺¹_i,j,kΩi,j,k =α^∗∗_i,j,kn

Ωi,j,k+h

(u₃)_i,j,k+1/2−(u₃)_i,j,k−1/2i

∆t∆x1∆x2

o

−h

(α^∗∗u3∆t)_i,j,k+1/2−(α^∗∗u3∆t)_i,j,k−1/2i

∆x1∆x2

(2.4.20)

EI-LEをNSSと組み合わせる場合には，界面再構築アルゴリズムに若干の変更を伴う．

計算手順を以下に説明する．

x1 方向の移流には EI(Eulerian implicit)を適用する．EI とは，発散修正項中のαを陰 的に扱う移流方法である．すなわち，

α^∗−αⁿ

∆t + ∂αⁿu1

∂x₁ =α^∗∂u1

∂x₁ (2.4.21)

セル体積Ω_i,j について積分した上式を離散化する．

α^∗_i,j∆x1 =αⁿ_i,j∆x1−h

(αⁿu1∆t)_i+1/2,j−(αⁿu1∆t)_i−1/2,ji

+α^∗_i,j∆th

(u1)_i+1/2,j−(u1)_i−1/2,j i

(2.4.22) 整理すると，

α^∗_i,j =n

Ωi,j+ ∆th

(u1)_i+1/2,j −(u1)_i−1/2,j i∆x2

o₋₁ nαⁿ_i,jΩi,j−h

(αⁿu1∆t)_i+1/2,j−(αⁿu1∆t)_i−1/2i

∆x2

o (2.4.23)

となる．上式の幾何学的意味を，図 2.4.6(a)に示す流れ場を例に説明する．本図では，

u_i+1/2,j < 0, u_i+1/2,j >0であり，隣接セル(i+1, j), (i−1, j)の両方からセル(i, j)に流体が

流入する．対象とする流体は非圧縮であるが，計算の都合上流体はセル(i, j)内に納まる ように圧縮される(図2.4.6(b))．圧縮される流体量と圧縮係数は各々，

αⁿ_i,jΩ_i,j−h

(αⁿu₁∆t)_i+1/2,j−(αⁿu₁∆t)_i−1/2,ji

∆x₂ nΩi,j + ∆th

(u₁)_i+1/2,j−(u₁)_i−1/2,ji

∆x2

o₋₁

である．したがって，α^∗_i,j は0< α^∗_i,j < 1の範囲に納まる．

次に，x2 方向の移流を行うが，x2 方向にはLE(Lagrangian explicit)を適用する．x1 方 向の移流計算時に圧縮された流体体積は，この移流によって膨張して本来の体積に戻る．

図2.4.7のように，x₂方向の速度u₂について，(u₂)_i,j+1/2 >0, (u₂)_i,j−1/2 < 0とする．体積

u_i−1/2,j u_i+1/2,j

i−1 i i+1

j

∆x₁−[(u₁)_i+1/2,j−(u₁)_i−1/2,j]∆^t

|(u₁)_i−1/2,j∆t| |(u₁)_i+1/2,j∆t|

i

i−1 i+1

j

|(u₁)_i−1/2,j∆^{t| |(u}1)_i+1/2,j∆^t|

(a)

(b)

図2.4.6 EI(Eulerian implicit)

の移流量を求めるが，セル(i, j)内で界面再構築を行うのではなく，膨張後に流体が占め

る体積，[∆x2+(u2)i,j+1/2∆t−(u2)i,j−1/2∆t]∆x1 内で行う．再構築完了後に，Ωi,j 外に配置

したサブセル内に含まれる体積は隣接セルに移流される．ただし，再構築する際の制約条 件は次式である．

f (φ_int)

=

NC

X

l=1

α_l^sub(φint)+

N_2p

X

m=1

α^sub_m (φint)+

N_2m

X

n=1

α^sub_n (φint)− α^∗_i,j Ωi,j

nΩi,j+h

(u2)_i,j+1/2 −(u2)_i,j−1/2

i∆t∆x1

o

= 0

(2.4.24) また，x2 方向移流時には法線の成分を変更する必要がある．まず，式(2.4.7)で法線を計 算する．求めた法線のx1 方向成分n1に次の量を掛ける．

1 Ωi,j

nΩ_i,j +h

(u₂)_i,j+1/2−(u₂)_i,j−1/2i

∆t∆x₁o

その後，大きさが１のベクトル，すなわち単位ベクトルになるようにn1,n2を規格化する．

|(u₂)_i,j−1/2∆t|

i j−1

j+1

j

|(u₂)_i,j+1/2∆t|

(u₂)_i,j+1/2

(u₂)_i,j−1/2

図2.4.7 LE(Lagrangian explicit)

本節で開発したNSSの体積保存性およびシャープな界面の維持に関する性能評価は次 章で行う．

ドキュメント内 2007 2 (ページ 73-80)