減圧環境における上昇気泡の数値解析

Synopsis:

In order to give appropriate assessment of designed systems for secondary refining processes, the predictive methods with computational fluid dynamics are strongly desired. In the present study, a numerical method, called ‘CIP-LSM’ (CIP based Level Set & MARS), was developed to simulate incompressible and compressible two-phase flows including surface tension and gravity effect. With the developed code, the dynamic behavior of bubbles rising in a liquid pool under reduced-pressure condition were tried to be simulated. As to the bubble shape, rising velocity and volume expansion, good agreement was obtained between the computed results and corresponding experiments.

1. 研究背景と目的

転炉や電気炉のような酸化精練炉から出鋼される溶鋼が鋳造されるまでの間に，温度と組成の均一化や介在物の除去および形態制御を行う過程は二次精錬と呼ばれ，高純度鋼製造のための標準工程となっている。

中でもRuhrstahl-Hausen真空脱ガス法(RH)は，高純度鋼の大量生産における二次精錬を代表する手法である。RHでは，真空槽に接続された2本の管を取鍋内の溶鋼中に浸したうえで，片方の管(上昇管)にアルゴンガスを吹込む。すると，比重差によるエアリフトポンプ効果のため，溶鋼は真空槽内へ導かれ，脱炭，脱窒素，脱水素などの処理を連続的かつ高効率に行うことができる。真空処理を経た溶鋼は，他方の管(排出管)を通って取鍋へ戻され，再び真空槽との間を循環する。

近年，鋼材の高純度化が進む中，更に高効率な介在物除去装置を開発し運用するためには，その設計段階から，二次精錬設備内部における熱流動特性を定量的に予測する技術が求められる。しかしながら，減圧環境に置かれた溶鋼の内部を運動する気泡の挙動を直接的に可視化するのは難しく，速度場や温度場などを計測するのも容易でない。

このような技術課題に対する実験的な接近法として，水を模擬液体とした気液二相流の観察実験が精力的に行われてきたが，それに加えて，蒸気圧の小さなシリコーン油を模擬液体とし，減圧容器内の油槽に気体を吹込み，鉛直方向の大きな圧力比に伴って膨張しながら浮上する気泡の挙動を観察した研究例¹⁾も報告されている。

その一方で，理論と実験を補完する手段として，気液界面を直接的に解像する数値解析法の確立も強く期待されている。実際，1990年代後半以降，自由表面流の分野でも，経験的な構成方程式やモデルをできるだけ排除した数値解析法が精力的に提案され^2,3,4)，界面張力と濡れ性を考慮して，大変形しながら移動する気液界面を精度良く予測できる界面捕獲法の技術が確立されつつある。

このような背景から，Kumagaiらは，減圧容器内のシリコーン油槽を膨張しながら上昇する単一気泡を詳細に観察する実験に取り組んでいる⁵⁾。一連の実験では，数値解析法の検証に役立てられることを視野に入れ，気泡挙動の高い再現性を得るべく，容器内部の雰囲気圧制御や，気泡放出時の装置動作が適切に管理されている。本論文では，これと対応して，近年開発された自由表面流解析手法であるCIP-LSMを用い，減圧液相内単一上昇気泡挙動の再現を試みた結果を報告する。対応する実験結果との比較に基づき，流れ場を詳しく記述するとともに，同解法の適用性と課題について議論する。

2. 数値解析法

2・1 界面追跡法と界面捕獲法

一般に，自由表面流数値解析の難しさは，液面形状が先決的に与えられないところにある。解法は，液面を移動格子やマーカー粒子を用いてラグランジュ的に追いかける界面追跡法と，固定格子に定義された識別関数の等高面として液面をオイラー的に捕獲する界面捕獲法に大別される。移動格子法は液面位置を陽に獲得できるため，界面張力の正確な評価が要求される解析に適用され，静止液体中を浮力の効果で上昇する気泡の数値解析にその成功例を見る。先駆的な例としては，気泡の上昇速度と抵抗係数を軸対称粘性計算により整理した研究⁶⁾を挙げることができる。これに対し，界面捕獲法は，液面形状の再現に劣るとされてきたが，界面張力を体積力として表現するCSFモデル⁷⁾と，識別関数として距離関数を用いるLevel Set法⁴⁾の出現により，その欠点は大きく改善された。また，界面捕獲法の一種であるVOF法は，体積保存に優れ，液面の合体や分裂も容易に表現できる反面，Fig.1に模式的に示すように，数値粘性に起因して液面が散逸する欠点を抱えていたが，PLIC-VOF法⁸⁾をはじめとして，液面を横切る密度の跳びを1格子幅で捕獲する手法が提案されている。

Fig. 1.

 Dissipation of VOF due to numerical viscosity.

2・2 CIP-LSM

研究背景で述べたように，二次精錬の真空処理が行われる溶鋼の液面近傍では，通常の大気圧下とは異なり，深さ方向に圧力が数倍のオーダーで変化するため，気相の体積変化を無視できず，非圧縮流の仮定は不適切となる。流動の速さは各相の音速と比べてかなり遅く，圧縮性流体解法を適用するのも難しい。伝熱現象を考慮するためには，温度を従属変数とせず直接的に解くのが望ましく，熱流束をより適切に評価するためには，界面近傍における物性値の平滑化を避け，人工的な界面厚を排除したい。

様々な環境に置かれた自由表面流を解析対象とできるよう開発された数値解法(CIP-LSM)は，流体解法として，CIP-CUP法⁹⁾の一種で，温度を独立変数とするTCUP(Thermo CIP-CUP)法¹⁰⁾を，界面捕獲法としてPLIC-VOF法の一種であるMARS法³⁾を採用しており，上述の要求に適合できる可能性がある。

離散化手法の詳細は，文献¹¹⁾に詳述している。本報では解法を構成するアルゴリズムのうち，気液の状態方程式を適切に考慮した流体解法と，界面捕獲法について，付録の章で説明する。

3. 解析結果

前章および付録に述べたる数値解析手法を用い，熊谷らの実験結果を参照し，大気圧環境ならびに減圧環境において浮力の効果で上昇する単一気泡の再現を試みた結果を述べる。

3・1 実験概要

まず，解析対象とした実験の概要を説明する。

実験に供した透明アクリル製の減圧容器は，Fig.2に示すように，内径100 mm，高さ200 mmの円筒形で，底面から155 mmまでシリコーン油(密度ρ=935 kg/m³，界面張力係数σ=20.1 mN/m)を満たした。容器底部の中心には空気を導入するノズルが配置されている。ノズル上方には，透明樹脂製で直径40 mmの半球形状をしたダンピングカップが伏せられた状態で設置されており，ノズルを通過した空気はカップ内で合体し，所定体積の単一気泡を得られるようになっている。カップの頂部は，減圧容器の中心軸上にあり，液面から100 mmの深さに位置している。カップは，その中心を貫く水平な軸に支持されており，遠隔駆動によりカップを軸周りに反転させると，気泡を放出できる。

Fig. 2.

 Experimental set-up⁵⁾. Test vessel and damping cup.

一方，減圧容器の上蓋にはバッファ槽を経て真空ポンプに繋がる配管が接続されており，液面位置での背景圧力を1.0 atm(大気開放)から0.013 atmの範囲で保持できる。

放出された単一気泡は，浮力によりシリコーン油の中を上昇し，その後，液面に達して破裂する。この間の変化は高速度カメラにより撮影され，画像から気泡形状と鉛直方向位置が判別される。

実験は，背景圧P₀，初期気泡体積V_i，液体動粘度νを様々に変えて行われたが，本報では，実験結果のうち特徴的ものとしてTable 1に掲げる3例を選び，その結果を数値的に模擬した。

Table 1.  Conditions of test cases.

	P0 (atm)	Vi (cm3)	ν (cSt)
Case 1	1.0	5.0	1.0
Case 2	0.013	1.0	10
Case 3	1.0	1.0	1.0

3・2 計算条件

実験に於ける円筒側壁の影響は小さいと仮定し，計算では，Fig.3に示すように，幅×奥行×高さ=L×L×2L(L=100 mm)の直方体容器を計算領域とした。初期気泡体積が比較的大きな場合であるCase 1では，Fig.4(a)に示すように，水平断面内に51×51点，高さ方向に201点の格子を配置した。他方，初期気泡体積が比較的小さな場合であるCase 2およびCase 3では，Fig.4(b)に示すように，鏡面対称性を仮定して，水平断面の1/4領域(L/2×L/2)だけを解き，その面内に41×41点，高さ方向に201点の格子を配置した。いずれの場合も，格子は気泡周辺で側壁近傍より密になるようになっている。なお，最小格子幅はCase 1で0.50 mm，Case 2およびCase 3では0.20 mmの程度である。

Fig. 3.

 Computational domain.

Fig. 4.

 Computational grid on horizontal plane.

作動流体は，非圧縮流体であるシリコーン油と理想気体の状態方程式に従う空気を用いた。初期条件として，ダンピングカップ球面の内側に接するように，所定体積を持つ気相領域を与えた。温度と圧力をT₀=288.15 KとP₀でそれぞれ一様に与え，対応する粘性係数や界面張力係数などの物性値を算出した。境界条件として，容器側壁を断熱の粘着壁，対称面(Case 2およびCase 3)を断熱の滑り壁とした。重力加速度は9.81 m/s²であり，計算開始と共に体積力として鉛直下向きに加えられ，下壁の存在を表現した境界条件とともに，鉛直方向の圧力(静水圧)分布が流れ場に形成される。静水圧分布と気液密度差が干渉した結果として浮力の効果が表現される。

3・3 大気圧・大気泡(Case 1)

まず，大気圧条件における特徴的な単一上昇気泡の例として，容器を減圧せずにP₀=1.0 atmとし，気泡初期体積を5.0 cm³，動粘度を10 cStとしたCase 1を模擬した結果を述べる。

計算結果で得られた界面形状を識別関数H_Vのゼロ等高面により可視化し，対応する実験結果と共にFig.5に示す。実験で得られた画像からは，ダンピングカップからの気泡放出時刻を定義して同定するのは困難で，且つ，放出直後の気泡形状もカップの動きに影響されると報告されていることから，今回は，容器底面から155 mmに位置するシリコーン油の液面位置に気泡が到達する瞬間を基準時刻t=0.00 secと定義し，前後の相対時刻を揃えて気泡形状を比較した。図中，実験と同じ容器側面方向に加え，斜め上方からの可視化結果も示している。

Fig. 5.

 Dynamic deformation if rising bubble. Case 1:P₀ = 1.0 atm, V_i = 5.0 cm³, 10 cSt: Big bubble in ambient pressure.

計算結果からは，カップから放出された直後の気泡が，(b)過渡的に半球状のSkirted形状を呈し，続いて，(c)上部のSpherical-cap部分と下部の環状気泡に分裂，さらに，(d)環状気泡が周方向に分裂，(e)液面に近づくと気泡下面側が窪み水平方向の視直径が大きくなり，(f)液面に到達して破裂するとともに，随伴流によって液面が隆起する様子まで，よく再現されていることが見てとれる。

以上の結果から，界面捕獲法でしばしば採用される人工的な界面厚を排除し，界面を横切る物性値の跳びを1格子幅以内で捕捉した場合であっても，界面張力，重力，慣性力の影響を受けて大きく変形する気泡の挙動を，分裂と合体も含め，安定かつ適切に模擬できる可能性を確認できた。

3・4 減圧・小気泡(Case 2)

続いて，初期気泡体積を1.0 cm³まで小さくし，容器をP₀=0.013 atmまで減圧した場合(Case 2)について同様の模擬を行った。計算結果で得られた界面形状を同様に可視化し，対応する実験結果と共にFig.6に示す。

Fig. 6.

 Dynamic deformation if rising bubble. Case 2:P₀ = 0.013 atm, V_i = 1.0 cm³, 10 cSt: Small bubble in reduced pressure.

まず気泡形状に注目すると，計算結果では，カップから放出された直後の気泡が，(b)過渡的に半球形状を呈した後，周期的な変動をする形状振動を経て，最終的に安定な(d)Spherical-cap形状に至り，その後，(e)液面に近づくと気泡下面側が窪み，(f)液面に到達して破裂する様子まで，概ね模擬されている。しかし，実験結果では気泡の分裂が認められないのに対し，計算では，気泡が半球形状からSpherical-cap形状に移行する間，僅かではあるが，最外周の尖った部分から分裂が認められる。形状振動が減衰して安定形状に至るまでの時間は，実験と比べて計算の方が長い。

次に，計算結果から気泡重心を算出し，容器底面からの位置と鉛直方向速度を描けば，Fig.7を得る。同図に打点された実験結果と比較すると，安定形状に至った後の終端速度について，良い一致が認めることができる。

Fig. 7.

 Rising velocity against to bubble position. Case 2:P₀ = 0.013 atm, V_i = 1.0 cm³, 10 cSt.

続いて，減圧した場合(Case 2)の計算で得られた気泡の体積について，初期気泡体積で除した膨張比を定義し，容器底面からの測った気泡重心の位置に対して描くとFig.8を得る。図中，減圧せずにP₀=1.0 atmとした場合(Case 3)の計算結果も示されているが，膨張比に明らかな違いが認められる。液面からダンピングカップ頂部までの浴深は100 mmであるから，初期気泡が感じるシリコーン油の静水圧は約0.009 atm(917 Pa)である。減圧したCase 2では，背景圧P₀と静水圧が同じオーダーとなっているので，気泡が放出されて液面に到達する間に，気泡周囲の圧力比は約0.022 atmから0.013 atmまで，約0.4倍にまで減少する。このため，気泡は上昇につれて大きく膨張することが理解できる。また，気泡内部の空気が理想気体の等エントロピー変化をすると仮定し，静水圧比から膨張比を算出すると点線のようになる。このことから，計算により模擬された気泡は，断熱的な変化を経たと推定される。

Fig. 8.

 Volume Expansion against to bubble position. Case 2:P₀ = 0.013 atm, V_i = 1.0 cm³, 10 cSt. Case 3: P₀ = 1.0 atm, V_i = 1.0 cm³, 10 cSt.

以上の結果から，減圧環境で周囲静圧が大きく分布する場合についても，気体の体積変化を考慮しつつ，気泡形状や終端上昇速度について，概ね妥当な数値解を得られる可能性が示された。その一方，気液界面の曲率半径が格子幅と同等以下になる気泡の最外周部から，実験では観察されない気泡分裂が認められた。これは，今回採用した界面捕獲法では，格子幅より小さいスケールの気泡，液滴，気膜，液膜を正しく捕獲できないことに起因すると考えられる。

気液間の密度比が約5000倍にも達するような場合についても，人工的な界面厚を排除した安定な計算が行えたことは，今後，液面を挟む熱移動と物質移動のサブグリッドモデルを考案するうえで重要であることを強調したい。

4. 結言

本論文では，様々な環境における自由表面流予測手法として開発が進められているCIP-LSMを用い，大気圧および減圧環境における単一上昇気泡挙動の数値的模擬を試みた。

最終的に気泡が Spherical-cap形状を呈するような場合について，実験で観察された気泡形状と終端速度，および，減圧環境での体積膨張を，数値解析でも概ね再現することができた。また，気液の密度比が5000倍を超える減圧環境においても，界面厚を排除した計算は安定に進行した。

格子幅より小さな寸法の液面形状が正しく捕獲されず，実現象と異なる気泡分裂が計算で発生するのが認められた。格子幅や時間刻みの適切な設定，界面捕獲法ならびに界面張力モデルの改良に関しては，更なる検討が必要であることが認識された。

記号

ρ:密度 [kg/m³]

p:圧力 [Pa]

T:温度 [K]

e:内部エネルギ [J/kg]

$\overrightarrow{u}$ :速度 ［m/s］

$\overrightarrow{g}$ :重力加速度 ［m/s²］

$\overrightarrow{q}$ :熱流束 ［W/m²］

C_S:音速(∂p/∂ρ)_S [m/s]

C_p:定圧比熱 [J/(kg・K)]

B:体膨張係数 [kg/(m³・K)]

μ:粘性係数 [Pa・s]

k:熱伝導係数 [W/(m・K)]

σ:界面張力 [N/m]

ξ，η，ζ:計算空間の座標 [m]

付録

CIP-LSMを構成するアルゴリズムのうち，流体解法と界面捕獲法の概要を，以下で詳しく説明する。

A1 支配方程式

一般に，固定格子を用いた自由表面流の解析手法は，均質二相流に関する質量，運動量，ならびに，内部エネルギの式，   

$$\frac{D\rho }{Dt}=-\rho \nabla \cdot \overrightarrow{u}$$(1)

  

$$\rho \frac{D\overrightarrow{u}}{Dt}=-\nabla p+\nabla \left(T_v+T_{\sigma }\right)+\rho \overrightarrow{g}$$(2)

  

$$\rho \frac{De}{Dt}=-p\nabla \cdot \overrightarrow{u}+\dot{\Theta }$$(3)

  

$$T_v=-\left(2\mu /3\right)\left(\nabla \cdot \overrightarrow{u}\right)I+\mu \left(\nabla \overrightarrow{u}+\nabla {\overrightarrow{u}}^T\right)$$(4)

  

$$T_{\sigma }=\sigma {\delta }_S\left(I-{\overrightarrow{n}}_S{\overrightarrow{n}}_S\right)$$(5)

  

$$\dot{\Theta }=\left(T_v+T_{\sigma }\right)\nabla \cdot \overrightarrow{u}-\nabla \cdot \overrightarrow{q}+\dot{L}$$(6)

を解く流体解法の部分と，識別関数H_Sの移流方程式，   

$$\frac{D{H}_S}{Dt}=\frac{\partial H_S}{\partial t}+\left(\stackrel{\rightharpoonup }{u}\cdot \nabla \right)H_S=0$$(7)

  

$$\begin{array}{c}\text{where}{H}_S=0.5：その点が液相に属する場合\\ H_S=0：その点が界面に属する場合\\ H_S=-0.5：その点が気相に属する場合\end{array}$$(8)

を解いて移動境界の位置と形状を特定する界面捕獲法の部分とで構成される。上記のT_vは粘性応力テンソル，T_σは界面張力テンソルであり，${\overrightarrow{n}}_S$は界面上で液相側へ向けた単位法線ベクトルである。なお，本報では，気液間の相変化を考慮していない。

A2 流体解法

局所熱平衡を仮定して熱力学関係式を用いると，式(1)~(3)は，(u, v, w, T, p)を状態量とする方程式系   

$$\frac{D\overrightarrow{u}}{Dt}=-\frac{\nabla p}{\rho }+\frac{1}{\rho }\nabla \left(T_v+T_{\sigma }\right)+\overrightarrow{g}$$(9)

  

$$\frac{DT}{Dt}=T\frac{B}{\rho C_p}\frac{Dp}{Dt}+\frac{\dot{\Theta }}{\rho C_p}$$(10)

  

$$\frac{1}{\rho C_S^2}\frac{Dp}{Dt}=-\nabla \cdot \overrightarrow{u}+B\frac{\dot{\Theta }}{\rho C_p}$$(11)

に変換できる。ここに現れた物性値は，   

$$B=-\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial T}\right)}_p：体膨張係数$$(12)

  

$$\frac{1}{\rho C_S^2}=\frac{1}{\rho }{\left(\frac{\partial \rho }{\partial p}\right)}_T-\frac{T}{\rho C_p}{B}^2：音速 \left(C_s\right)$$(13)

であり，いずれも気液各相に対する状態方程式，   

$$\rho =\left(0.5+H_S\right)\cdot {\rho }_{Liq}\left(T,p\right)+\left(0.5-H_S\right)\cdot {\rho }_{Gas}\left(T,p\right)$$(14)

を決めれば直ちに(T, p)の関数として導出される。

任意の状態方程式に従う作動流体の場合，式(10)を考慮した，温度と圧力の連成問題となるが，TCUP法では，温度と圧力についてSIMPLE法に類似した外部反復を行うことで，式(9)~(11)を同時に満足する解を探索することができるように設計されている。

なお，作動流体が非圧縮の場合はC_s=∞かつB=0であるから，TCUP法のアルゴリズムは非圧縮流体解法のHSMAC法と一致する。

A3 MARS法による体積率の移流

式(7)の離散化にあたっては，各格子点に対応した検査体積Ω毎にH_Sを体積積分し，これに含まれる各相の体積率H_Vを，   

$$\frac{H_V}{J}=\underset{\Omega }{\int }{H}_{S}dV=\underset{k-1/2}{\overset{k+1/2}{\int }}\underset{j-1/2}{\overset{j+1/2}{\int }}\underset{i-1/2}{\overset{i+1/2}{\int }}{H}_S\frac{d\xi d\eta d\zeta }{J}$$(15)

と定義して，独立な状態量とする。即ち，H_Vは[−0.5, 0.5]を定義域とするボイド率である。続いて，Fig.9に描くように，一般座標ξの一定面でΩを切り取った断面を考え，この断面の面積分率H_ξを，   

$$\frac{H_{\xi }}{J}=\underset{\delta \xi \to 0}{\lim }\left\{\frac{1}{\delta \xi }\underset{k-1/2}{\overset{k+1/2}{\int }}\underset{j-1/2}{\overset{j+1/2}{\int }}\underset{\xi -\delta \xi /2}{\overset{\xi +\delta \xi /2}{\int }}{H}_S\frac{d\xi d\eta d\zeta }{J}\right\}$$(16)

Fig. 9.

 Re-construction of Interface; Tangential plane of zero level set and the interpolation of H_ξ in ξ -direction based on MARS¹¹⁾.

により定義する。そのうえで，MARS法に則り，|ξ−i|≤1/2の区間におけるH_ξの分布を，近似的に，   

$$H_{\xi }=\max \left\{-0.5,\min \left[0.5,{H^{\prime }}_{\xi }\left(\xi -{\xi }_S\right)\right]\right\}$$(17)

と区分的一次関数で内挿補間し，方向別に再構成する。式中の界面勾配H'_ξは，気液界面の単位法線ベクトル${\overrightarrow{n}}_S$から計算される。すると，Fig.9(a)に示すように，有限時間Δtの間にΩの表面を通過するH_Vの流束を，区分的一次関数の台形積分により算出できる。

体積率に加えて界面勾配も考慮することで，移流計算に伴うH_Vの非物理的な数値散逸を回避しつつ，気液界面を横切る密度の跳びを1格子幅で捕獲できるという点が重要である。

A4 Level Set法による界面形状捕捉

上述の界面勾配H'_ξを精度良く算出するために，形状捕捉に優れるLevel Set法を援用する。即ち，Fig.10に示すように，界面を基準とした距離関数φを生成し，${\overrightarrow{n}}_S$および界面曲率κを，それぞれ，   

$${\overrightarrow{n}}_S=\frac{\nabla \varphi }{\left|\nabla \varphi \right|}$$(18)

  

$$\frac{\kappa }{J}=\underset{\Omega }{\int }\nabla \cdot {\overrightarrow{n}}_{S}dV=\underset{\partial \Omega }{\int }{\overrightarrow{n}}_S\cdot d\overrightarrow{S}=\sum _{\hat{\xi }=\xi ,\eta ,\zeta }\left[\frac{\nabla \varphi }{\left|\nabla \varphi \right|}\cdot \frac{\nabla \hat{\xi }}{J}\right]$$(19)

Fig. 10.

 Normal vector ${\overrightarrow{n}}_S$ of interface calculated from the distribution of φ¹¹⁾.

に従って獲得する。φを効率的に生成(再初期化)するためには，独立変数H_Vのゼロ等高面として表現される界面近傍において，従属変数であるφに関するHamilton-Jacobi型方程式，   

$$\{\begin{array}{c}\frac{\partial \varphi }{\partial \tau }+\left(S_{LSM}{\overrightarrow{n}}_S\right)\cdot \nabla \varphi =S_{LSM}\text{and}{S}_{LSM}=\mathrm{sgn}\left(\varphi \right)\\ \left(\text{if}\left|H_V\right|=0.50\right)\end{array}$$(20)

  

$$\begin{array}{c}\varphi =H_V\cdot \left(\left|N_{\xi }\right|+\left|N_{\eta }\right|+\left|N_{\zeta }\right|\right)\text{and}{S}_{LSM}=0\\ \left(\text{if}\left|H_V\right|<0.50:界面近傍\right)\\ 但し,N_{\hat{\xi }}\equiv \left(x_{\hat{\xi }},y_{\hat{\xi }},z_{\hat{\xi }}\right)\cdot {\overrightarrow{n}}_S\text{for}\hat{\xi }=\xi ,\eta ,\zeta \end{array}$$(21)

を擬似時間進行法により反復して収束させれば良い。その際，式(21)のような埋込境界条件を課せば，H_Vのゼロ等高面とφのゼロ等高面を良好に一致させることができる。

このように，Level Set法とMARS法を相補的に用いることで，形状捕捉に優れたLevel Set法の長所を損なうことなく，気液両相の体積保存を良好に満足することができる。

A5 物性値と界面張力

自由表面流の数値解法では，計算を安定化する便宜として，界面近傍でρなどの物性値を人工的に平滑化する操作が施されることがある。これに対し，今回の提案では，各格子点における物性値を，例えば密度について，   

$$\rho =\left(0.5+H_V\right)\cdot {\rho }_{Liq}\left(T,p\right)+\left(0.5-H_V\right)\cdot {\rho }_{Gas}\left(T,p\right)$$(22)

のように，体積率H_Vの加重平均で算出するよう変更した。また，その他の物性値Xについても同様に，   

$$X=\left(0.5+H_V\right)\cdot X_{Liq}\left(T,p\right)+\left(0.5-H_V\right)\cdot X_{Gas}\left(T,p\right)$$(23)

の形で与え，界面を横切る物性値の跳びを1格子幅以内で捕獲されている。流体解析の計算で必要となる各物性値について，Xの形をまとめ，Table 2に掲げる。

Table 2.  Physical properties X in Eq.(23).

Physical properties	X
ρ: Density	ρ
Cp: Specific heat at constant p	ρCp
CS: Sound velocity	(ρCs2)–1
B: Volume expansion coefficient	B
μ: Viscosity coefficient	μ–1
K: Heat conductivity	k–1

他方，界面張力の局所合力は，曲率κを用いて，   

$$\nabla \cdot T_{\sigma }=-\sigma \kappa \frac{d{H}_S}{d\varphi }{\overrightarrow{n}}_S=-\sigma \kappa \nabla H_S\cong -\sigma \kappa \nabla H_V$$(24)

と近似できるので，   

$$-\nabla \cdot {\left(\frac{\nabla \cdot T_{\sigma }}{\rho }\right)}_{\left(m-1\right)}^{\left(n+\lambda \right)}\delta t^{\left(n\right)}=\nabla \cdot {\left(\sigma \kappa \frac{\nabla H_V}{\rho }\right)}_{\left(m-1\right)}^{\left(n+\lambda \right)}\delta t^{\left(n\right)}$$(25)

のように，圧力項と同じく，反変加速度として運動量式に加えた。

文献

• 1)   T.  Tatsuoka and  K.  Ito: CAMP-ISIJ, 7 (1994), 1100.
• 2)   T.  Yabe and  E.  Takei: J. Phys. Soc. Jpn., 57 (1988), 2598.
• 3)   T.  Kunugi: Computational Fluid Dynamics Journal, 9 (2000), No.1.
• 4)   M.  Sussman,  P.  Smereka and  S.  Osher: J. Comput. Phys., 114 (1994), 146.
• 5)   T.  Kumagai,  Y.  Kishimoto,  O.  Nakamura,  T.  Tagawa,  T.  Himeno and  S.  Shimazaki: Tetsu-to-Hagané, 101 (2015), 93.
• 6)   G.  Ryskin and  L.G.  Leal: J. Fluid Mech., 148 (1984), 19.
• 7)   J.U.  Brackbill,  D.B.  Kothe and  C.  Zemach: J. Comput. Phys., 100 (1992), 335.
• 8)   D.L.  Youngs: Numerical Methods for Fluids Dynamics, ed. by K.W.Morton and M.J.Baines, Academic Press, Massachusetts, (1982).
• 9)   T.  Yabe and  E.  Takei: J. Phys. Soc. Jpn., 57 (1988), 2598.
• 10)   T.  Himeno and  T.  Watanabe: Trans. JSME Ser. B, 69 (2003), No.678, 226.
• 11)   T.  Himeno,  H.  Negishi,  S.  Nonaka,  C.  Inoue,  T.  Watanabe and  S.  Uzawa: Trans. JSME Ser. B, 76 (2010), No.765, 778.