20 50 4 20 50 12 4
8 4
12 8
8 4
8
図2.5.1 曲率ベクトル評価制度検証用回転体
図2.5.2 κnの真値
図2.5.3 体積率αを評価法Aに適用して得たκn
16, 19.5, 21.6とした．各形状の模式図を図2.5.1に，κnの真値を図2.5.2に示す．αを評 価法Aに適用して得たκnを図2.5.3に，評価法Bに適用して得たκnを図2.5.4に示す．
評価法Aを用いた場合，曲率ベクトルの大きさκのみならず，その方向nにも大きな誤 差が生じている．評価法Bは，評価法Aに比べて計算精度がよい．本結果から，カラー 関数にαを用いる場合，より良い表面張力評価を行うためには，計算負荷の大きい評価法 Bを用いなければならないことが確認できる．
ι ‑戸
o . . . . . . '
r
， .. ' .
図2.5.4 体積率αを評価法Bに適用して得たκn
図2.5.5 レベルセット関数φを評価法Aに適用して得たκn
図2.5.6 レベルセット関数φを評価法Bに適用して得たκn
一方，レベルセット法[33]では，次式，
|∇φ|=1 (2.5.14)
を満足する符号付距離関数(レベルセット関数)φをC とすることで，αよりも精度よく κnを評価できる．また，式(2.5.14)より式(2.5.6)右辺大括弧内第１項はゼロとなる．し たがって式(2.5.6)は式(2.5.5)と等価になる．このため，図2.5.5, 2.5.6に示すように，φ を用いた場合，評価A, Bのいずれを用いてもほぼ同じ良好な結果が得られる．
以上の検討結果から，κを計算するためのカラー関数には，表面張力が作用する界面領 域(３セル幅程度)内で，
|∇C|= const. (2.5.15)
を満足する関数が適切であるといえる．
2.5.2 局所レベルセット関数による表面張力評価
前節の考察より，良好に表面張力を評価するには，αではなくレベルセット関数のよう
に|∇C|= const.を満足するカラー関数を使用する必要があることがわかった．本節では，
体積率αから界面近傍のセルにのみレベルセット関数φを生成し，そのφを用いて表面 張力評価を行う方法を提案する．この界面近傍にのみ生成するφを局所レベルセット関 数と呼ぶ．以下にφの算出法を説明する．
界面セル(i, j,k)(0 < αi,j,k <1)のφを，2.4節のNSSアルゴリズムにより求まるφint で 与える．
φi,j,k =φint (if 0< αi,j,k <1) (2.5.16)
φint に基づいて，界面セル周辺のα = 1およびα= 0の計算セルのφを求める．図 2.5.7 に示す，界面セル(i, j,k)とその隣接セル(i, j,k+1)(αi,j,k+1 = 0)を考える．この場合，セ ル(i, j,k+1)のφは，セル(i, j,k)のφint に基づいて次式で与える．
φi,j,k+1 = h
ri,j,k+1−
ri,j,k+φintni,j,ki
·ni,j,k (2.5.17)
ここで，rはセル中心の位置ベクトルである．一般には，セル(i, j,k+1)の周囲にはセル (i, j,k)以外にも界面セルがあるため，それら界面セルのφint も考慮する必要がある．そこ で，次の重み付き平均によりφi,j,k+1 を算出する．
φi,j,k+1 = P
l,m,n∈NS×NS×NS Wl,m,nφl,m,n
P
l,m,n∈NS×NS×NS Wl,m,n (2.5.18)
ここで，φl,m,n は，式(2.5.17)により界面セル(l,m,n)からセル (i, j,k+1)に対して計算
されたレベルセット関数である．NS は重み付き平均を行う範囲内にあるセル数である．
Wl,m,n は重み関数であり，次式で与える．
Wl,m,n =( ri,j,k+1−rl,m,n−p if 0< αl,m,n < 1
0 otherwise (2.5.19)
ここで，pは経験的に与えるパラメーターである．
以上の方法で求めた局所レベルセット関数φを用いて，表面張力Fσを次式で評価する．
Fσ= σκ(φ)∇Hε(φ)
|∇φ| (2.5.20)
φintni,j,k O
ri,j,k ri,j,k+1
|(ri,j,k+1–rint).ni,j,k|
ni,j,k φi,j,k+1=(ri,j,k+1–rint).ni,j,k
(i,j,k+1)
(i,j,k)
interface
fluid 1 fluid 2
ri,j,k+φintni,j,k
図2.5.7 局所レベルセット関数φの算出
ここで，Hε は平滑化されたヘビサイド関数であり，次式で与えられる．
Hε(φ)=
0 if φ <−ε
1 2
1+φ/ε+(1/π) sin (πφ/ε)
if |φ| ≤ε
1 if φ > ε
(2.5.21) εはHεが記述する界面の幅を規定するパラメーターで，1.5∆x程度にとる．
より高精度な局所レベルセット関数φ，すなわち距離関数としての性質(|∇φ|= 1)をよ り厳密に満足するφが必要とされる場合には，上述の手順に加えて以下に述べる再初期化 を行う．
次の方程式の定常解は距離関数としての性質(|∇φ|=1)を満足する[33]．
∂φ
∂τ = signε(φ0) (1− |∇φ|) (2.5.22)
ここで，τは擬似的な時間，φ0 は再初期化開始前のφの値，signε は次式で定義される平 滑なsign関数である．
signε(φ0)= φ0 q
φ20+ ∆x2
(2.5.23)
式(2.5.22)を次の形に変形する．
∂φ
∂τ +w· ∇φ=signε(φ0) (2.5.24)
ここで，
w=signε(φ0) ∇φ
|∇φ| (2.5.25)
式(2.5.24)左辺には移流方程式の解法が適用できる．本研究では，時間微分項に２次精度
ルンゲ・クッタ法，移流項にENO(essentially non-oscillatory)スキーム[34]を用いる[35]． ２次精度ルンゲ・クッタ法では，まず中間段階m+1/2のφを次式で求める．
φm+1/2 =φm+ ∆τL (φm) (2.5.26)
ここで，L(φm)=signε(φ0)(1− |∇φm|)とおいた．L(φ)中の|∇φ|は後述するENOスキーム により計算する．次に，上式で得たφm+1/2よりL(φm+1/2)を計算し，次ステップの値φm+1 を次式で計算する．
φm+1 =φm+ ∆τ 2
hL (φm)+L(φm+1/2) i
(2.5.27) L(φ)中の|∇φ|= [(∂φ/∂x1)2+(∂φ/∂x2)2 +(∂φ/∂x3)2]1/2 をENOスキームに基づいて計算 する．例えば∂φ/∂x1 は次式で与えられる．
∂φ
∂x1
!
=
∂φ
∂x1
+
if
∂φ
∂x1
+
signε(φ0)<0 and
∂φ
∂x1
−
signε(φ0)< −
∂φ
∂x1
+
signε(φ0)
∂φ
∂x1
−
if
∂φ
∂x1
−
signε(φ0)>0 and
∂φ
∂x1
+
signε(φ0)> −
∂φ
∂x1
−
signε(φ0)
0 if
∂φ
∂x1
−
signε(φ0)<0 and
∂φ
∂x1
+
signε(φ0)> 0
(2.5.28) (∂φ/∂x1)+, (∂φ/∂x1)− は以下の諸式で定義した．
∂φ
∂x1
!−
= φi,j,k−φi−1,j,k
(∆x1)i−1/2 + ∆x1 2 m
Dx1φi,j,k,Dx1φi−1,j,k
(2.5.29)
∂φ
∂x1
!+
= φi+1,j,k−φi,j,k (∆x1)i+1/2 + ∆x1
2 m
Dx1φi,j,k,Dx1φi+1,j,k
(2.5.30) ここで，Dx1 は差分演算子であり，これをφに作用させて得られる差分式を以下に示す．
Dx1φi,j,k = 1 (∆x1)i
"
φi+1,j,k−φi,j,k
(∆x1)i+1/2 − φi,j,k−φi−1,j,k (∆x1)i−1/2
#
Dx1φi−1,j,k = 1 (∆x1)i−1
"
φi,j,k−φi−1,j,k
(∆x1)i−1/2 − φi−1,j,k−φi−2,j,k
(∆x1)i−3/2
#
Dx1φi+1,j,k = 1 (∆x1)i+1
"
φi+2,j,k−φi+1,j,k
(∆x1)i+3/2 − φi+1,j,k−φi,j,k
(∆x1)i+1/2
#
また，関数m(a,b)は次式で定義した．
m (a,b)=
( a if |a| ≤ |b|
b otherwise (2.5.31)
x2, x3方向も同様に計算する．
上述の再初期化過程において，界面セルのφはφintに固定する．こうすることで，レベ ルセット法における再初期化の場合に生じる，φの記述する界面が動いてしまうという問 題を解決できる．
2.5.3 秩序変数による表面張力評価
本節では，フェーズフィールドモデル[36–38]において界面を記述するために用いられ る秩序変数をカラー関数として利用する方法を提案する．2.5.2節では，界面再構築時に 界面セルの距離関数が求まるという NSSの特徴を活かして，局所レベルセット関数によ る表面張力評価法を提案したが，秩序変数はカーン‐ヒリヤード方程式 [39]を差分法で 解けば得られるため，NSS以外の体積追跡法にも応用できる．
まず，秩序変数の関数形が，表面張力評価に適しているか否かを調べる．フェーズ フィールドモデルでは，界面は秩序変数ϕにより数セル幅の領域として記述される．ϕの 時間変化は次のCahn-Hilliard方程式に従う．
∂ϕ
∂t =∇ ·Γ(ϕ)∇η (2.5.32)
ここで，Γはオンサガー係数，ηは化学ポテンシャルである．ϕはバルク流体中で最小値 ϕminあるいは最大値ϕmax を，界面領域ではその間の値をとり，相界面を記述する．ϕmin
およびϕmaxは，平衡条件，
η(ϕmin)=η(ϕmax) (2.5.33)
を満足するように決定される．系の自由エネルギー汎関数F [40]を，
F = Z
Ω
f (ϕ)+ Kε
2 |∇ϕ|2
dΩ (2.5.34)
と定義すると，化学ポテンシャルηは，
η= ∂f
∂ϕ −Kε∇ · ∇ϕ (2.5.35)
と与えられる．ただし，f はバルク自由エネルギー，Kε は界面厚さを規定するパラメー ターである．式(2.5.35)を式(2.5.32)に代入して次式を得る．
∂ϕ
∂t =∇ ·Γ(ϕ)∇
"
∂f
∂ϕ −Kε∇ · ∇ϕ
#
(2.5.36)
-30 -20 -10 0 10 20 30 0.30
0.35 0.40
x
ϕ
K
ε= 0.05 K
ε= 0.10 K
ε= 0.20 K
ε= 0.40
ϕ
maxϕ
min図2.5.8 秩序変数ϕの関数形
簡単のため，式(2.5.36)に支配される ϕの平衡解の特性を１次元の場合について説明す る．シャープな界面(x < 0でϕ(x) = ϕmin, x > 0でϕ(x)= ϕmax)を記述するϕ(x)を初期 値として式(2.5.36)を解くと，その平衡解は，図2.5.8に示すtanh関数となる．また，多 次元の場合，界面の厚みはKε1/2 に比例し，いたる所で一定である．この秩序変数分布は，
精度よく曲率評価を行うためにカラー関数に求められる性質(|∇C|=const.)を界面中央近 傍で近似的に満足する．そこで，体積率αをϕの初期値としてCahn-Hilliard方程式を解 き，その時間進行解をC とする方法を提案する．本手法では，瞬時のαから擬似的な時 間進行によりϕを生成する．よって，式(2.5.36)における時間tを，擬似的な時間τに置 き換える．ファンデアワールス流体モデル[38]に基づき，f を次式で与える．
f = ϕ
"
RT ln ϕ
1−bϕ −ϕa
#
(2.5.37) また，Γ(ϕ)= Γ0ϕ(Γ0 は定数)とおくと，式(2.5.36)は，
∂ϕ
∂τ =∇ ·Γ0ϕ ("
RT
ϕ(1−bϕ)2 −2a
#
∇ϕ− ∇Kε∇ · ∇ϕ )
(2.5.38) となる．ただし，Rは気体定数，T は温度，a, bは各々粒子間の長距離引力作用，短距離 反発作用の強さを表わす定数である．本手法では，式(2.5.38)中のパラメーターおよび定 数は，所定の界面厚さを実現する目的のためだけに使用し，物理的な意味は付与しない．
Kεを変更することで異なる界面厚さが実現される様子を図2.5.8に示す．
曲率評価の手順を以下に示す．
(1) α の値をϕmin から ϕmax の間に補間して与え，ϕ の初期値とする．本研究では，
ϕ= (ϕmax−ϕmin)α+ϕminと与える．
(2) 式(2.5.38)を解き，所定の界面厚さまで拡散したϕを得る．本研究では，時間差分 に１次のオイラー時間進行法，空間差分に４次精度中心差分を用いて式(2.5.38)を 離散化する．
(3) ϕをカラー関数とし，式(2.5.7)あるいは式(2.5.10)を用いて曲率κを計算する．
x1方向の４次精度の補間および中心差分は以下の諸式で与えられる．
ϕi = −ϕi−3/2+9ϕi−1/2+9ϕi+1/2−ϕi+3/2
16 (2.5.39)
∂ϕ
∂x1
!
i
= ϕi−3/2−27ϕi−1/2+27ϕi+1/2−ϕi+3/2
24∆x1
(2.5.40)
∂2ϕ
∂x21
i
= ϕi−3/2−ϕi−1/2−ϕi+1/2+ϕi+3/2
2 (∆x1)2 (2.5.41)
x2, x3方向についても同様である．
曲率を精度よく評価するために複数のセルに拡散させたϕを式(2.5.2)中の∇Ci,j,kの評 価に用いると，広範囲に表面張力が作用してしまう．そこで，曲率評価のために求めたϕ とは別に，よりシャープな界面を記述するϕを生成し，それを用いて∇Ci,j,k の評価を行 うことで表面張力の作用する領域を制限することができる．
ここで提案した二種の表面張力評価法の精度検証は次章で行う．
2.6 結論
本章では，基本特性(体積保存・シャープな界面の維持・高精度表面張力評価)を備える 界面追跡法の開発を目的として，サブセルに基づく界面追跡法NSS(non-uniform subcell
scheme)を開発した．NSSの特徴は，(a)計算法が単純な既存手法の多くは計算セル内の
界面の傾斜を無視するのに対し，NSSは界面の傾斜を考慮している点，(b)界面の傾斜を 考慮する既存手法の多くは複雑な幾何計算を要するが，NSSは計算セル内に一時的に配 置するサブセルを利用した簡単なアルゴリズムに基づくため，容易に３次元化できる点に ある．
体積追跡法では，CSF(continuum surface force)モデルを用いて表面張力評価を行う際，
体積率をカラー関数とすると精度が低いことが課題であった．そこで，CSFモデルのカ ラー関数として適切な関数形について考察した．その結果，以下の結論を得た．
(1) 勾配の大きさが界面近傍で一定である関数がカラー関数として望ましい．
(1)の考察に基づいて，カラー関数として望ましい性質を有する関数を採用した，以下 の二種の表面張力評価法を提案した．(a)NSSのアルゴリズムを利用して界面近傍に局所 レベルセット関数を生成し，カラー関数として利用する方法．(b)カーン‐ヒリヤード方 程式を解いて秩序変数を生成し，カラー関数として利用する方法．