第 2 章 の参考文献
3.2 表面張力評価精度の評価
3.2.1 回転体に対する曲率評価
提案した二種の手法による界面曲率評価の妥当性を検証するため，図3.2.1に示す球，
回転楕円体，冠球形状の曲率ベクトルκnを計算した．これら回転体形状に対するκnの真
値を図3.2.2に示す[1]．各手法が含むパラメーターの設定値を以下に説明する．
3.2.1.1 局所レベルセット関数算出時のパラメーター
次式で与えられる界面セル(i, j,k)の局所レベルセット関数φi,j,k，
φi,j,k =φint (if 0< αi,j,k <1) (3.2.1)
に基づいて，界面セル近傍のセルのφを次式で与える．
φi,j,k = P
l,m,n∈NS×NS×NS Wl,m,nφl,m,n
P
l,m,n∈NS×NS×NS Wl,m,n (3.2.2)
.••.
ι ‑戸
o . . . . . . .
~.〆
ここで，NS は平均をとる領域のセル数，Wl,m,n は次式で定義される重み関数である．
Wl,m,n =
( ri,j,k− rl,m,n−p if 0< αl,m,n <1
0 otherwise (3.2.3)
rはセル中心の位置ベクトル，αは体積率である．以下の計算では，NS =7, p= 16を用 いた．再初期化を実施する場合，∆τ=0.375∆x，再初期化回数を4とした．
3.2.1.2 秩序変数算出時のパラメーター
ファンデアワールス流体モデルを適用したカーン‐ヒリヤード方程式[2]，
∂ϕ
∂τ =∇ ·Γ0ϕ ("
RT
ϕ(1−bϕ)2 −2a
#
∇ϕ− ∇Kε∇ · ∇ϕ )
(3.2.4) において，各パラメーターを次のように設定した．R = 1, T = 0.293, a = 1, b = 1． この場合，秩序変数の最大値，最小値は，ϕmax = 0.405, ϕmin = 0.265 となる．ϕ が平 衡解となるまで上式を時間進行させると計算コストが高くなる．短い時間進行によ る時間進行解によっても，精度よく曲率を評価できることが重要である．そこで，本 計算では時間進行回数を 10 とした．ϕ が平衡解となるのに要する時間進行回数を少 なくするため，安定に計算を進められる範囲で拡散係数を大きくとる．本計算では，
Kε =0.4, Γ0 =32, ∆τ=0.01,∆x1 = ∆x2 = ∆x3 = 1.0 (Γ0∆τ/(∆x1)2 = 0.32)とした．
3.2.1.3 曲率ベクトル計算結果
以下，曲率評価に次の式(3.2.5)を用いるものを評価法A，式(3.2.6)を用いるものを評 価法Bと称する．
κ=−∇ · ∇C
|∇C|
!
(3.2.5) κ= 1
|∇C|
"
∇C
|∇C|· ∇ |∇C| − ∇ · ∇C
#
(3.2.6) ここで，Cはカラー関数であり，φあるいはϕを代入する．
φを評価法A, B に適用して得た回転体の曲率ベクトルκnを図3.2.3, 3.2.4，再初期化 を施したφを評価法A, Bに適用して得たκnを図3.2.5,3.2.6，ϕを評価法A, Bに適用し て得たκnを図3.2.7, 3.2.8に示す．また，Cとして体積率 αを用いて得たκnを図3.2.9, 3.2.10に示す．
図から，φ, ϕを用いれば，A, Bどちらの評価法を用いても，各回転体形状に対して，α を用いるよりも誤差の小さい曲率ベクトルが得られることが確認できる．これは，φとϕ が，良好な表面張力評価を行うためにカラー関数に求められる特性，
|∇C|= const. (3.2.7)
図3.2.3 局所レベルセット関数φを評価法Aに適用して得たκn
図3.2.4 局所レベルセット関数φを評価法Bに適用して得たκn
図3.2.5 再初期化を施したφを評価法Aに適用して得たκn
を有しているためである．図3.2.3, 3.2.4と図3.2.7, 3.2.8を比べると，φを用いた方が，ϕ を用いるよりも曲率の絶対値が大きくなる箇所の精度が高いことがわかる．これは，カー ン‐ヒリヤード方程式(3.2.4)の右辺が曲率の大きい部分の回転体形状をやや平滑化して しまうためと考えられる．本研究では，局所レベルセット関数を利用した方法により，
NSSに高精度表面張力評価性能を付与する．ただし，ϕを利用した曲率評価法は，αを直 接用いるよりも大幅に精度を改善でき，また通常の差分法で式(3.2.4)を解くだけでよい ためNSS以外の体積追跡法にも容易に組み込めるという長所を有している．
図3.2.6 再初期化を施したφを評価法Bに適用して得たκn
図3.2.7 秩序変数ϕを評価法Aに適用して得たκn
図3.2.8 秩序変数ϕを評価法Bに適用して得たκn
3.2.2 静止液中単一液滴計算
体積率αを直接カラー関数として用いて，垂直ダクト内に満たされた静止液中を上昇 する単一液滴を計算した．計算条件を以下に示す．エトベス数 Eo = 0.22，モルトン数 M =1×10−3，密度比ρD/ρL =0.79，粘度比µD/µL =0.36．下付添字Lは周囲流体，Dは
図3.2.9 体積率αを評価法Aに適用して得たκn
図3.2.10 体積率αを評価法Bに適用して得たκn
液滴を意味する．エトベス数およびモルトン数は以下の諸式で定義した．
Eo= g (ρL−ρD) d2
σ (3.2.8)
M = gµ4L(ρL−ρD)
ρ2Lσ3 (3.2.9)
ここで，gは重力加速度，dは液滴球等価直径，σは表面張力である．エトベス数は液滴 が受ける浮力と表面張力の比を表わす．本計算では，表面張力評価誤差が顕著に現れるよ うに，浮力に比べて表面張力が大きい条件（Eoの小さい条件）に設定している．また，液 滴直径dあたりのセル数d∗(= d/∆x)は12とした．計算体系は，x1, x2, x3 軸方向にそれ ぞれ3d, 3d, 8dの直方体とし，セル幅は∆x1 = ∆x2 = ∆x3 で一定とした(図3.2.11)．計算 領域上部から液滴終端上昇速度の大きさで一様に周囲流体を流入させ，液滴が計算領域内 に長時間とどまるようにした．計算領域下部には自由流出条件を，側壁には−x3方向に動 く移動壁条件を課して無限静止液中液滴の条件を模擬した．図3.2.12に計算結果を示す．
αをカラー関数として用いると表面張力誤差が大きいため，液滴表面近傍に非物理的な渦 構造が生じている．
ぶ明、
ぷ(す ふ
x1 x2 x3
drop
top boundary: uniform inflow bottom boundary:
continuous out flow side boundaries:
moving boundary domain:
3d x 3d x 8d cell size:
∆x1=∆x2=∆x3 spatial resolution:
d*=12 g
図3.2.11 擬似流れ評価用計算体系
図3.2.12 体積率αを表面張力評価に使用して得た液滴周囲の流れ場
次に，φおよびϕを利用した表面張力評価法を用いて同様の計算を行った．その結果を
図3.2.13,図3.2.14に示す．これらの図から，表面張力評価法を改善することによって，
物理的により妥当な流れ場が得られることを確認できる．
図3.2.13 局所レベルセット関数φを表面張力評価に使用して得た液滴周囲の流れ場
図3.2.14 秩序変数ϕを表面張力評価に使用して得た液滴周囲の流れ場