第 3 章 の参考文献
4.1 緒論
4.3.5 流れ場の解法
本研究では，Yabe らが提案した C-CUP(CIP-combined, unified procedure) 法 [22] の 従属変数，速度・圧力・内部エネルギーを，速度・圧力・温度に変更した姫野らの T-CUP(thermo-CUP)法[23]に基づいて式(4.3.22), (4.3.23)を解く．ただし，本研究で扱う 二相系は非圧縮かつ等温に限るため，計算過程はオリジナルの T-CUP法よりも簡略化さ れる．以下にその計算過程を示す．なお，実際に開発した計算コードではエネルギー方程 式も加味しており，非圧縮，等温などの仮定は条件分岐によって取り入れている．エネル ギー方程式も含んだT-CUP法については付録を参照されたい．
4.3.5.1 T-CUP法
T-CUP法では，フェイズ分割の概念に基づいて，連続の式および運動方程式を移流フェ
イズ，拡散フェイズ，音響フェイズに分割する．
移流フェイズでは，次の速度ベクトルuの移流方程式を解く．
u{∗}−un
∆t +(u· ∇u)n =0 (4.3.24)
上式はCIP法により直接解くため，有限体積的には扱わない．
拡散フェイズでは，次式を解いて，分子粘性および表面張力による運動量変化を計算 する．
Z
Ω
ρ{∗}u{∗∗}−u{∗}
∆t dΩ = Z
S
(τ−Tσ){∗}·dS (4.3.25)
音響フェイズでは，次の圧力ポアッソン方程式を解き，n+1時刻の圧力を得る．
Z
S
∆t
ρ{∗∗}∇pn+1 ·dS= Z
S
u{∗∗}+ ∆t g
·dS (4.3.26)
次式により速度場をu{∗∗}から発散なしの条件を満足するun+1 に更新する．
Z
Ω
ρ{∗∗}un+1−u{∗∗}
∆t dΩ =− Z
S
pn+1I·dS+ Z
Ω
ρ{∗∗}gdΩ (4.3.27)
4.3.6 変数配置および近傍セルの表記法
変数配置にはコロケート配置[24]を用いる(図4.3.5)．すなわち，速度ベクトルのデカ ルト成分u1,u2,u3，圧力 p，体積率α，レベルセット関数φをセル中心に配置するが，連 続の式を満足するように反復修正される速度成分は各セル表面上に配置する．セル表面に
x1 x3 x2
U1
ξ1 ξ2 ξ3
U3
U2
u1, u2, u3, p, α, φ
図4.3.5 コロケート配置
配置される速度は曲線座標ξi の反変速度Uiである．以下，ξi 座標の成分を大文字で記述 する．Ui は反変成分の変換則に従う．
Ui = ∂ξi
∂xjuj (4.3.28)
ただし，セルの面積ベクトルdSは式(4.3.12)で定義されているため，dSを利用して，
JUi
±ξi = (u·dS)±ξi (4.3.29)
と計算し，(JUi)+ξi を配列に記憶する．Ui は常にJUi の形で基礎方程式中に現れるため，
この方が都合がよい[25]．
ある量がセルのどの位置に定義されているのかは，下付添字により区別する．セル中心 の量には大文字のローマ字，セル表面の量には小文字のローマ字を付す(図4.3.6)．具体 的には，セルC から見て座標ξ1, ξ2, ξ3 が増える方向に位置するセル表面を各々下付添字 e, n, t (east, north, topより)，減る方向に位置するセル表面を下付添字w, s, b (west, south,
bottomより)で表わす．セルCの各隣接セルは，セルC と共有するセル表面の記号の大
文字で表わす．例えば，セル表面eを共有する隣接セルはEである．また，Cとはセル表 面を共有しないが，Cの２つの隣接セルとセル表面を共有するセルが，C の斜め方向に存 在する．このようなセルには，Cの２つの隣接セルを示す記号を両方付す．例えば，Eと Nの両方とセル表面を共有するセルを ENと書く．
以上の記法を用いると，例えば，セル表面eにおける圧力 pのξ1 方向の偏導関数の離 散式は次式で与えられる．
∂p
∂ξ1
!
e
= pE −pC
∆ξ1 = pE− pC
∵ ∆ξ1 =1
x1 x3 x2
C
E W
e w
ξ1 ξ2 ξ3
図4.3.6 隣接セル
4.3.7 圧力ポアッソン方程式と速度修正
音響フェイズにおける圧力ポアッソン方程式(4.3.26)を離散化すると次式となる．
−CCpn+1C +X C pn+1
neighbor =X
f
h
u{∗∗}+ ∆t g
·dS i
f (4.3.30)
左辺のCは，圧力勾配ベクトルの面積分を離散化する際に現れる係数で，面積ベクトルの 成分，ξi のxi についての偏導関数，密度ρ，時間刻み幅∆tを含む．左辺第２項のP
は，
セルC のすべての近傍セル(neighbor)に対して和をとる．圧力勾配ベクトルの離散式と 係数C の詳細は付録を参照されたい．ポアッソン方程式は標準的な反復解法で解ける．
本研究ではSOR法を用いて解いた．
ポアッソン方程式の解 pn+1 により，セル表面の反変速度 JUi を修正する．例えば，セ ル表面eにおける速度(JU1)eの修正式は次式となる．
JU1
n+1
e =
JU1 {∗∗}
e −
"
∆t
ρ ∇pn+1+ ∆t g
!
·dS
#
e
(4.3.31) 次に，セル中心Cに定義された速度のデカルト成分u1, u2, u3 の修正法について説明す る．セル表面の反変速度とは異なり，セル中心の速度に対しては連続の式を満足するとい う拘束条件が適用されておらず，その修正法には任意性がある．セル中心速度の修正法に は，以下の２つが考えられる．
4.3.7.1 (A)速度補間
連続の式を満足するセル表面上の反変速度成分をセル中心に補間し，その後デカルト座 標成分に変換する．すなわち，速度修正を次式により行う．
(ui)n+1 = 1 JC
∂xi
∂ξj
JUj n+1
C (4.3.32)
各成分を書き下すと，以下のようになる．
(u1)n+1 = 1 JC
"
∂x1
∂ξ1
JU1 n+1
C + ∂x1
∂ξ2
JU2 n+1
C + ∂x1
∂ξ3
JU3 n+1
C
#
(4.3.33)
(u2)n+1 = 1 JC
"
∂x2
∂ξ1
JU1 n+1
C + ∂x2
∂ξ2
JU2 n+1
C + ∂x2
∂ξ3
JU3 n+1
C
#
(4.3.34)
(u3)n+1 = 1 JC
"
∂x3
∂ξ1
JU1 n+1
C + ∂x3
∂ξ2
JU2 n+1
C + ∂x3
∂ξ3
JU3 n+1
C
#
(4.3.35) ここで，
JUi
n+1
C =
JUi
+ξin+1 +
JUi
−ξin+1
2 (4.3.36)
4.3.7.2 (B)速度修正量補間
反変速度の速度修正量をセル中心に補間した後，デカルト成分に変換してセル中心速度 に加える．各セル表面における反変速度の修正量を次式により計算する．
∆JUi
±ξi =−
"
∆t
ρ ∇pn+1+ ∆t g
!
·dS
#
±ξi
(4.3.37) セル中心Cに修正量を補間する．
∆UCi =
∆JUi +ξi
+
∆JUi
−ξi 2JC
(4.3.38) 修正量∆UCi をデカルト成分に変換する．
∆ui = ∂xi
∂ξj∆Uj (4.3.39)
速度修正を次式で行う．
un+1i =u{∗∗}i + ∆ui (4.3.40)
前者 (A) の場合，連続の式を満足する反変速度をセル中心に補間するため，結果と して得られるデカルト成分 ui は連続の式を大きく破らないと期待できる．しかし，梶 島[25]は，以下の問題点を指摘している．更新後のセル表面の反変速度Ui は，例えば U1e
u{∗∗}C , u{∗∗}E
のように，隣接セルの速度から補間して求めた中間量に，修正量∆U1 を 加えたものである．u{∗∗}E は運動方程式を解いて求めるが，その際Cから離れた E の隣接 セルの量を用いて計算されているため，最終的に求められた更新後のセル表面の反変速度 U1e には，遠方の計算点の情報が入り込んでいる．よって，前者の方法では補間を重ねた 結果解がぼやける．このことは，例えば水中気泡などを計算する場合に，気泡に働く浮力 や界面に働く表面張力が本来作用しない遠方のバルク流体中に作用してしまうことを意味 する．一方，後者(B)の場合，セル表面で求めた速度修正量そのものをセル中心に補間し て速度修正を行うため，前者(A)の方法のような問題は生じない．よって，本研究では後 者(B)の速度修正法を採用する．