2020 年 1 月 24 日

2019

(

)

推進薬タンク内スロッシングに対する 磁場による抑制効果の数値解析

2020 ^年 1 ^月 24 ^日

首都大学東京大学院

システムデザイン研究科 システムデザイン専攻 航空宇宙システム工学域 博士前期課程

学修番号 18863635 氏名 古市 侑太郎

指導教員 田川 俊夫 准教授

目次

第1章 緒言 1

1.1 ^研究背景 . . . 1

1.2 ^研究目的 . . . 2

第2章 計算手法 3 2.1 ^{解析モデル・計算格子}. . . 3

2.1.1 ^{解析モデル} . . . 3

2.1.2 ^計算格子. . . 4

2.2 ^{二相流のモデル化} . . . 4

2.2.1 ^{各相の取扱い} . . . 4

2.2.2 ^{表面張力の算出} . . . 5

2.3 ^{界面捕獲法} . . . 7

2.3.1 Level Set^法 . . . 7

2.3.2 VOF (Volume of Fluid)^法 . . . 8

2.3.3 CLSVOF (Coupled Level Set and VOF method)^法 . . . 15

2.4 ^{支配方程式} . . . 16

2.4.1 ^連続の式. . . 16

2.4.2 Navier-Stokes^方程式 . . . 16

2.4.3 ^{移流方程式} . . . 17

2.4.4 Biot-Savart^の法則 . . . 17

2.4.5 磁場に対するガウスの法則 . . . 19

2.5 ^境界条件 . . . 20

2.6 IB(Immersed Boundary)^法 . . . 20

2.7 ^{無次元支配方程式} . . . 21

2.8 ^{圧力修正計算法} . . . 22

2.9 ^{計算パラメータ} . . . 22

第₃章 結果及び考察 ₂₄ 3.1 ^{妥当性検証} . . . 24

3.1.1 ^{液滴の振動} . . . 24

3.1.2 ^{気泡の上昇} . . . 25

3.2 磁場によるスロッシング抑制効果の解析 . . . 28

3.3 ^{今後の課題と方針} . . . 31

第₄章 結言 ₃₂

Appendix A 体積補正 33

Appendix B 無次元化 ₃₅

B.1 Navier-Stokes^方程式 . . . 36

B.2 ^{物性値の定義} . . . 38

B.3 連続の式，磁場に対するガウスの法則 . . . 38

B.4 ^{移流方程式} . . . 38

B.5 ^{無次元支配方程式} . . . 39

Appendix C Multigrid法 40 C.1 二相流における圧力ポアソン方程式 . . . 40

C.2 圧力ポアソン方程式の離散化と係数行列の設定 . . . 41

C.3 幾何的マルチグリッド法の導出. . . 42

C.4 2^{段グリッド法} . . . 43

記号表

Nomenclature

A(t) 各時間における無次元の界面近傍領域の体積 [-]

b_c コイル中心での磁束密度 [T]

b 修正後の磁束密度ベクトル [T]

b_Biot Biot-Savartで求めた磁束密度ベクトル [T]

B 修正後の無次元磁束密度ベクトル [-]

BBiot Biot-Savartで求めた無次元磁束密度ベクトル [-]

C VOF^関数 [-]

C_err ^{体積誤差率} [-]

D ^{無次元コイル直径} [-]

eg 重力ベクトル [-]

fmag 磁化力ベクトル [N/m³]

Fmag 無次元磁化力ベクトル [-]

f_{sur f} ^{表面法線力ベクトル} [N/m³]

F_{sur f} 無次元表面法線力ベクトル [-]

g ^{重力加速度の大きさ} [m/s²]

Ga ^{ガリレイ数} [-]

H ^{無次元コイル高さ} [-]

H_α^scaled(ϕ) Density-scaled Heaviside^関数 [-]

i_c ^{コイル電流値} [A]

l ^代表長さ(^{タンク直径}) [m]

L ^{無次元代表長さ}(^{タンク直径}) [-]

M 磁場強度を表す無次元数 [-]

p ^圧力 [Pa]

P ^{無次元圧力} [-]

r ^{コイル半径} [m]

r コイル位置ベクトル [m]

R 無次元コイル位置ベクトル [-]

t ^時間 [s]

u ^{流速ベクトル} [m/s]

U ^{無次元流速ベクトル} [-]

U_B 壁面における無次元流速ベクトル [-]

UI B IB点における無次元流速ベクトル [-]

UI P IP点における無次元流速ベクトル [-]

U_{I P} UI P のX^方向成分 [-]

V_{I P} U_{I P} ^のY ^方向成分 [-]

V(t) 各時間における無次元の液相体積 [-]

Vinit 無次元の初期液相体積 [-]

x 計算領域内の位置ベクトル [m]

X 計算領域内の無次元位置ベクトル [-]

X_n,Y_n,Z_n X, Y, Z^{方向の格子数} [-]

∆X,∆Y,∆Z X, Y, Z^{方向の格子幅} [-]

Greek letters

α ^{界面遷移幅} [m]

αˆ ^{無次元界面遷移幅} [-]

γ ^表面張力 [N/m]

Γ ^{ラプラス数} [-]

δI P IP^{点までの距離} [-]

κ ^界面曲率 [1/m]

κˆ ^{無次元界面曲率} [-]

µG 気相の粘性係数 [Pa·s]

µL 液相の粘性係数 [Pa·s]

µ¯ ^{粘性係数比} [-]

µϕ 無次元粘性係数 [-]

µm 真空の透磁率 [H/m]

π ^円周率 [-]

ρG 気相の密度 [kg/m³]

ρ^{L 液相の密度} [kg/m³]

ρ¯ ^密度比 [-]

ρ_ϕ ^{無次元密度} [-]

τ ^{無次元時間} [-]

τMAX 無次元計算時間 [-]

ϕ Level Set^関数 [m]

ϕp スカラポテンシャル [-]

Φ ^無次元Level Set^関数 [-]

χ^{G 気相の質量磁化率} [m³/kg]

χL 液相の質量磁化率 [m³/kg]

χ¯ ^{質量磁化率比} [-]

χϕ 無次元質量磁化率 [-]

ψ ^{壁面からの距離関数} [-]

第

1

緒言

1.1

人間の活動領域が地球周回軌道上にまで拡大するのに伴い，ロケットや宇宙探査機のような地上と は異なる加速度環境下で液体を取り扱う機会が増加している．更に近年，将来の安価で高頻度の宇宙 輸送需要に対応するため盛んに再使用型ロケットの開発がされており、日本においても、発射場の占 有面積や既存施設の流用の観点から、垂直離着陸式の再使用型ロケットが注目されている．また，地 上から低軌道に至る打ち上げロケットでは初段の部分再使用が実現されつつある．

推力や姿勢が変化する宇宙機での動的加速度環境では，容器内の自由表面を持つ液体が外部からの 加振によって揺動する現象(^{スロッシング})が推進薬タンク内で起き，タンク内の温度と圧力や宇宙 機の重心が変動する．そのため，推進薬の位置や温度を制御することがミッションの成功性に関わっ

てくる^[1][2]．推進薬は宇宙機の重量の大半を占めており，スロッシングの発生は機体の重心を変動さ

せ姿勢制御に悪影響を及ぼす他，推進薬を加圧するためのアレッジガスがタンク底部に到達すると気 相が配管に混入し，エンジンの失火や動作不良を招く危険性がある．さらに，タンク内部が極低温の 液相と比較的高温の気相が共存する環境である場合，スロッシングによって気液間の熱交換や相変化 が促進されタンク圧力が急降下し，推進薬流量の低下やターボポンプ入り口でのキャビテーションの 発生原因となるため，エンジンの動作不良に直結する危険性もはらんでいる．例えば, H-IIA^ロケッ トの第一段ではLOX(^液体酸素)をエンジンで熱交換させガス化した後,機体側へ還流してタンク加圧 ガスとして用いる方法を採用したが,一成分二相系となるタンク圧力制御と相変化量予測の難しさが 認識されている^[3][4][5]. ^{また，最近では}JAXA/ISASが開発する再使用型ロケット実験機においても, 飛行成立性にスロッシングが少なからず影響を及ぼすことが認識されている^[6].

現在のスロッシング対策の一つとしてバッフル板が挙げられる(Fig. 1.1)．バッフル板はタンク内 壁に設置されるドーナツ状の板で，スロッシングを減衰させるデバイスである．しかし，先行研究に よれば，バッフル板はスロッシング時の液面の運動力学的な変動を抑制できるが,^{気液間での熱的な} 変動は必ずしも抑制できるわけではなく,バッフル板が存在した方がそうでない場合より圧力変動の 勾配が大きくなる場合があり(Fig. 1.2)^[7]，その原因として，バッフル板内縁により液相が乱流となり 温度成層状態が崩されることで，下層の冷たいサブクール層が気液界面に露出して気液間の温度交換 が促進されているということが見出された(Fig. 1.3, 1.4)^[8]．これは，現在バッフル板が搭載されて いる宇宙機にとって飛行成立性に深刻な影響を与える可能性がある．もし，バッフル板を用いずにス ロッシングを抑制できれば温度成層状態は崩されず，スロッシングに伴う動的挙動と熱的挙動を同時 に制御できる可能性があるが，物理的なデバイス以外での宇宙機の推進薬制御は未だ知見の蓄積が少 ないのが現状である．

Fig. 1.1:Baffle plates in a tank^[9] Fig. 1.2:Pressure-time history^[7]

Fig. 1.3:Schematic of thermal stratification and equilib- rium state

Fig. 1.4:Schematic of turbulence flow due to baffle plates

1.2

そこで本研究では，ロケットに多用される推進薬の一つである液体酸素が常磁性体である性質に着 目し，液体酸素のスロッシングが発生した際，磁場による抑制効果にどの程度の有用性が望めるか，

数値計算を用いて評価を試みた．

第

2

計算手法

これより本論では太字はベクトルを示すこととするが，太字のナブラ演算子のみ無次元化を施した ナブラ演算子を表す．

2.1

2.1.1 解析モデル

本研究での解析モデルをFig. 2.1に示す．無重力状態における直径1 [m]の球形の宇宙機タンクを 模した容器中で，気体酸素と液体酸素の非圧縮性気液二相流を仮定し，x方向の加振力として最大 0.40 [G] (0.40×9.8[m/s²])^{の三角波状の衝撃力}(Fig. 2.2)を仮定した．なお，本計算では熱及び濡れ 性の影響は考慮しない．計算領域は球形タンクを含有する立方体とする．また，磁場は流体計算領域 外に存在する一巻きコイルの作る静磁場を想定し，コイルの位置やコイル中心での磁場の強さを様々 に変化させ，数値計算を行った．コイルの各条件をTable 2.1^{にまとめる．}

Fig. 2.1:Computational model Fig. 2.2:Acceleration

Table 2.1:Computational conditions for coils

Case No. of computational conditions Case 0 Case 1 Case 2 Case 3 Case 4 Dimensionless height of the coils (H) [-] -0.2 0.0 0.0 0.0 Ratio of coil diameter to characteristic length (D/L) [-] No coil 1.0 2.0 2.0 2.0 Magnetic flux density at the center of coils (bc) [T] 1.0 0.1 1.0 3.0

Fig. 2.3:Staggered grid

2.1.2 計算格子

空間の離散化には，デカルト直交座標系における等間隔スタッガード格子を用いる．Figure 2.3^に 示すように，セル中心(i, j)^に圧力p,^密度 ρ, Level Set^関数ϕ^やVOF^関数C^{などのスカラー量を定} 義し，セル界面(i+1/2,j), (i, j+1/2)^に速度u, v, w^{，磁束密度ベクトル}B_x,B_y,B_z ^{などのベクトル量を} 定義する．

2.2

2.2.1 各相の取扱い

本研究では気液二相流を想定しているため，自由界面を捕獲する必要がある．界面捕獲法として は,表面張力の算出精度が良いが体積保存性が悪いLevel Set^法[10],体積保存性に優れているが表面 張力の算出精度が悪いVOF^法(Volume Of Fluid method)^[11],^{それらを組み合わせ},^{互いの欠点を補い} つつ,^{利点を活かす}CLSVOF^法(Coupled Level Set and VOF method)^{[12]などがある}. ^{本研究で取り上} げる現象は界面形状が大変形し，かつ表面張力の影響が大きい無重力環境下を想定するので，体積保 存性と表面張力の精度を両立を図ってCLSVOF法を採用した．詳細は後述する．

また，界面は本来不連続面であり，特に気液の密度比が非常に大きい場合は，数値計算で不連続面 を解くと計算が不安定になる．そこで，密度などの物性値における気液界面の不連続面に対して一定 の厚みを持った遷移領域を設ける必要がある．以下に，遷移領域を持つ気液二相の物性値（密度 ρ^， 粘度 µ^，磁化率 χ）を表現した式を示す．

ρ= ρG+(ρL−ρG)H_α(ϕ) (2.1)

µ= µ^G+(µ^L− µ^G)H_α(ϕ) (2.2)

χ= χG+(χL− χG)H_α(ϕ) (2.3)

ここで，添え字L^は液相，G^{は気相を表し，}ϕ^{は界面からの}Level Set^{関数である．また，}H_α(ϕ) は遷移領域を持つ近似Heaviside^{関数であり，}αを界面厚さとする．本計算ではα=1.75∆x^とした．

Fig. 2.4:Smoothed Heaviside function and Density-scaled Heaviside function

H_α(ϕ)= 



0 (ϕ < −α)

1 2

[

1+ ^ϕ_α+ ¹_πsin (πϕ

α

)] (−α≤ ϕ≤ α)

1 (α < ϕ)

(2.4)

ただし，本研究では気液の密度差が大きいため，表面張力の取り扱いには Yokoi^{氏が提唱する} Density-scaled balanced CSF^{モデル[13]}を用いて計算の安定化を図った．H_α^scaled(ϕ)^はDensity-scaled Heaviside^{関数を表す．}

H_α^scaled(ϕ)=



0 (ϕ <−α)

1 2

[1

2 + ^ϕ_α + ₂^ϕ_α²2 − ₄¹_π2

{ cos

(2πϕ α

)−1

}+ ^α+ϕ_απ sin (2πϕ

α

)] (−α≤ ϕ≤ α)

1 (α < ϕ)

(2.5)

また，Fig. 2.4^にHeaviside^関数とDensity-scaled Heaviside^{関数を示す．}Density-scaled Heaviside 関数は通常のHeaviside関数より液相側に寄った形をしている．これにより気相側に掛かる表面張力 を弱めることで，高密度比の場合に界面付近に生じる非物理的な速度を抑制することができる．

2.2.2 表面張力の算出

二相流計算で重要なのが表面張力の算出である．表面張力は本来，二相界面において単位長さあた りに働く接線力であるが，温度が均一である場合には接線方向の力が打ち消しあい，結果として界面 の局所平均曲率に応じた法線方向の見かけ上の力だけが残る．界面における単位長さ当たりに働く法 線方向の力の大きさ fs[N/m²]は，次のように表される(^{下付き添え字}s^はsurface force (^面積力)^を 表す)^．

fs =γκ (2.6)

ここで，γ ^{は界面張力}[N/m],κ^{は局所平均曲率}[1/m]を表す．本来，界面は物理的に不連続である ので，界面の両側において流体の圧力や物性値（密度，粘度など）に跳びが生じる．そこで，Blackbill らの提案したCSF (Continuum Surface Force)^{モデル[14]を用いる．}CSFモデルでは，界面を一定幅 の遷移領域を持つ連続的な領域と捉える(Fig. 2.4^参照)．表面張力によりもたらされる応力は遷移領 域全体に対して分布的に作用するものであると考え，遷移領域を横切った時の法線力の積分値が，本 来の不連続面での圧力ジャンプγκに等しくなるように設定される．ただし，遷移領域の厚さは曲率

半径より充分小さいことが要求される．遷移領域の局所部分での法線力の次元は[N/m³]^{であり，遷} 移域に作用する法線力はLevel Set^{関数を変数とする}Heaviside^{関数，または}Heaviside^{関数の傾きを}

表すDiracのデルタ関数を用いて次式で与えられる．

fv =γκδ(ϕ)∇ϕ=γκ∇H_α(ϕ) (2.7)

δ(ϕ)= { 1

2α

{ 1+cos

(πϕ α

)} (|ϕ| ≤ α)

0 (|ϕ| > α) (2.8)

ここで，下付き添え字のv ^はvolume force (^体積力)^{を表す．式}(2.7)を界面遷移領域にわたって積 分すると，

f_s =

∫ _∞

−∞fv·dx=

∫ _∞

−∞γκδ(ϕ)∇ϕ·dx=γκ

∫ _∞

−∞δ(ϕ)dϕ=γκ (2.9)

となり，式(2.6)に一致することがわかる．ただし，Dirac^{のデルタ関数の性質}(^{全域で積分すると} 1)を用いていること及び積分区間にわたってκが定数とみなしていることに注意する．

さて，式 (2.7)のうち，どちらかを用いてプログラムを実装する必要がある．δ(ϕ)^{を用いた式が}

Standard CSF^モデル，H_α(ϕ)^{を用いた式は}Balanced CSF^{モデルである．}Balanced CSF^{モデルは離} 散化式が圧力項と同様の形になるので，圧力項の離散化誤差を打ち消しやすいという性質がある．ま た，H_α(ϕ)の代わりに前項で述べたDensity-scaled Heaviside^関数H_α^scaled(ϕ)^{を用いて離散化したも}

のがDensity-scaled balanced CSFモデルであり，本研究ではこちらを採用した．

また，界面曲率κを計算する必要がある．界面の局所平均曲率κは次式に従って界面法線方向単位 ベクトルnの発散を取ることで得られる．

κ =−∇ ·n=−∇ · ( ∇ϕ

|∇ϕ|

)

=−∇²ϕ (∵ |∇ϕ| =1) (2.10)

ここで，|∇ϕ| =1^はLevel Set^{関数の性質である}(^{次項で述べる})^．今，ϕ^{の下付き添え字}x^をx^に よる偏微分を表すとして，二次元デカルト座標での曲率(^式(2.11))及び三次元デカルト座標での曲率 (^式(2.12))^{を以下に示す．}

κ= 2ϕxyϕxϕy−ϕx xϕ²y−ϕyyϕ²x

(ϕ²x+ϕ²y

)³₂ (2.11)

κ= 2(

ϕxϕyϕxy+ϕyϕzϕyz+ϕzϕxϕz x

)−ϕx x

(ϕ²y+ϕ²z

)−ϕyy

(ϕ²z +ϕ²x

)−ϕzz

(ϕ²x+ϕ²y

) (ϕ²x+ϕ²y +ϕ²z

)³₂ (2.12)

ただし，この式は式(2.10)^の分母|∇ϕ|も微分していることに注意する(^つまり，|∇ϕ| =1^を用い ていない)^．|∇ϕ| =1^{を用いて簡易的に}

κ =−∇²ϕ=−(ϕx x+ϕyy+ϕzz) (2.13)

としても良いが，式(2.11), (2.12)と比較して精度が劣る可能性がある．

(a)Without re-initialization (b)With re-initialization Fig. 2.5:The effect of re-initialization

2.3

本研究ではCLSVOF法を界面捕獲法として用いることを前項で述べたので，本項ではその詳細に ついて述べる．CLSVOF^法はLevel Set^法とVOF法を併用したものなので，まずLevel Set^法とVOF 法単体についてそれぞれ説明したのち，CLSVOF^{法について解説する．}

2.3.1 Level Set法

Level Set^関数ϕを界面からの符号付き垂直方向距離を値として持つ距離関数として定義し，界面

を一定の厚みを持った遷移領域として表現し固定格子上でオイラー的に捕獲する．つまり，界面上で 0の値を取り，液相で正，気相で負の距離の値を持つ関数として定義される（もしくは液相が負，気 相が正と定義しても良い）．なお，本計算では液相が正，気相が負となるようにLevel Set^{関数を定義} した．このように定義されたLevel Set関数は，距離関数としての以下の性質をもつ．

|∇ϕ|=1 (2.14)

また，Level Set関数の移流方程式，つまり界面の移流方程式は次式で表される．

∂ϕ

∂t +(u· ∇)ϕ=0 (2.15)

相変化による界面の更新を考えないことにすれば，界面上に存在し続ける流体粒子がもつLevel Set^{関数の値は}0のまま不変である．したがって，移流方程式(2.15)^によってLevel Set^{関数を更新} することで新たなゼロ等高面へと移流する．しかしながら，式(2.15)^によってLevel Set^{関数を移流} すると，距離関数としての性質を失われてしまう（Fig. 2.5b^[15]）．そこで，再初期化という操作を行 うことで距離関数の性質を回復させる（次項で説明する）．

Level Set法の再初期化

移流方程式に従ってLevel Set関数を移流させると，距離関数としての性質を失い，さらに計算を 続けると界面厚さが不均一になり，やがて計算が発散してしまう恐れがある．そこで、各タイムス テップで式(2.15)を計算した後に再初期化という操作を施すことで、距離関数の性質を保つように計 算を進めることができる．具体的には、以下の式を用いて再初期化の手続きを行う．

∂ϕ

∂t˜ =SLS M(1− |∇ϕ|) (2.16)

S_{S L M} = ϕold

√ϕ²_old+ϵ² (2.17)

式(2.16)^{はそれぞれの格子で}|∇ϕ| = 1となる（距離関数の性質）ようにLevel Set^{関数にそのエ}

ラー値を加えていくことを意味している．また，式(2.17)^の添え字oldはその値が再初期化前のもの であることを意味する．これらの式で記述される再初期化の操作を疑似時間t˜^{の進行による反復計算} で収束させることで、距離関数の性質を回復させることができる（Fig. 2.5a^{[15]）．このとき、界面近} 傍では因子S_{LS M} ^は0に近づくので再初期化を施しても界面位置は動かずに距離関数の性質を回復 することができる．ただし，界面位置が動かないのは解析的な話であり，離散化して数値的に解くと 実際には界面が移動してしまうので，再初期化の反復計算回数を多くしすぎると実現象にそぐわない 計算結果となることに注意されたい．なお，ϵ は界面近傍でのゼロ割りを防ぐための値であり，界面 遷移幅α程度の値を設定する．実際に計算する場合、式(2.16)を移流方程式のように変形し，対流 項を風上差分法により評価する．

∂ϕ

∂t˜ +(W_{LS M}· ∇)ϕ)= S_{LS M} (2.18)

W_{S L M} =S_{S L M} ∇ϕ

|∇ϕ| (2.19)

式(2.18)を用いる理由としては，界面に近い側を上流として一方向に情報を伝達していくことで，

効率良く計算ができるからである．

2.3.2 VOF (Volume of Fluid)法

VOF^法はFig. 2.6^{に示すような},計算セル内に占める体積充填率をVOF^関数C^{を用いて定義し},

以下の移流方程式を解くことにより界面を求める界面捕獲法の一種である．

∂C

∂τ +∇· (CU) −C(∇·U)=0 (2.20)

充填率は0^〜1^{で定義され},充填率の値によって物性値を考慮する. ^なお,^充填率が0^〜1^の間では2 つの流体の物性値を充填率で加重平均して考える. 移流方程式の移流計算は有限体積的に行われるた め,体積を完全に保存しながら計算を進めることが可能となる. ^しかし, VOF^{関数は界面上でしか勾} 配を持たないので,法線ベクトルと界面曲率が不連続となり,表面張力を精度を精度良く計算するのは 難しいという欠点を持つ. Fig. 2.7^のような2次元の場合を考える．実際の計算では式(2.20)^を各軸 方向に次元分割を行い、図の丸1,^丸2のように順番に計算を行う．

Fig. 2.6:Schematic of VOF function

Fig. 2.7:Schematic of two-dimensional advection of VOF function

(2.20)^式中の第2項に対して有限体積法を適用し, セルを横切る流束として項を評価する. ^また,

(2.20)式は以下のように離散化される.

C_i^∗_,j =C_iⁿ_,j− F_xⁿ_,i+1/2,j−F_xⁿ_,i−1/2,j

∆x +∆tC_iⁿ_,j

uⁿ_i+1/2,j−u_iⁿ_−1/2,j

∆x (2.21)

C_iⁿ⁺¹_,j =C_i^∗_,j− F_y,i,ⁿ _j+1/2−F_y,i,ⁿ _j−1/2

∆y +∆tC_i^∗_,j

v_i,ⁿ_j+1/2−v_i,ⁿ_j−1/2

∆y (2.22)

ここで, F_x, F_y ^{はそれぞれΔ}t^の間にx ^面とy 面のセルを横切る流束であり, ^{以下のように計算} する．

F_x,i+1/2,j =−

∫ y_i,j+1/2 yi,j−1/2

∫ x_i+1/2,j−u_i+1/2,j∆t x_i+1/2,j

χi,jdxdy (2.23)

F_y,i,j+1/2 =−

∫ _{yi,j+1/2−vi,j+1/2∆t}

yi,j+1/2

∫ _xi+1/2,j

xi−1/2,j

χi,jdxdy (2.24)

Fig. 2.8:Schematic of one-dimensional surface

χは補間関数と呼ばれる関数であり，VOF関数と以下の関係を満たす．

Ci,j = 1

∆x∆y

∫ _yi,j+1/2

yi,j−1/2

∫ _xi+1/2,j

xi−1/2,j

χi,jdxdy (2.25)

つまり，定性的に説明すれば，補間関数をセル一つ分面積分すると，VOF^{関数に等しくなる．こ} の補間関数の構築方法には種々の方法があるが，本研究ではTHINC / WLIC^{法[11]を用いる．次項}

でTNINC^法及びWLIC^{法について述べる．}

THINC (Tangent of Hyperbola for INterface Capturing)法

THINC^法では,^補間関数 χ^をtanh^{を用いて表現する}. tanh^{を用いることで},^{非常に数値拡散を小}

さくすることができ,体積誤差を少なくすることができる. ^また, THINC^{法では界面を}1^{次元方向の} み考えるため,^以降ではx^{方向のみの}1次元問題として考えるものとする.

x^方向速度u^がu ≥ 0^のとき, Fig. 2.8^のような1^{次元の界面を考えると}, THINC^{法における補間} 関数は次のように表される.

χi= 1 2

[

1+αctanh

{β(x−x_i−1/2

∆x −x˜i

)}]

(2.26) ここで，αcは界面の法線方向により決定する．

αc = {

1 (C_i+1 ≥C_i−1)

−1 (Ci+1 <Ci−1) (2.27)

また，β^は界面のSmoothingパラメータであり，本研究では β=3.5としている．これで未知数は

˜

xi のみとなり，これを求めていく．

一次元の場合でも式(2.25)と同様に考えて，補間関数とVOF^{関数の間には} C_i = 1

∆x

∫ _xi+1/2,j

x_i−1/2,j

χidx (2.28)

という関係が成り立つので，これを用いると，

Ci = 1

∆x

∫ _xi+1/2,j

xi−1/2,j

χidx (2.29)

= 1

∆x

∫ _xi+1/2,j

xi−1/2,j

1 2

[

1+αctanh

{β(x−x_i−1/2

∆x −x˜i

)}]

dx (2.30)

= 1 2∆x

[

x+ αc∆x β ln

{ cosh

(β(x−x_i−1/2

∆x

)−x˜_i

)}]^xi+1/2,j x_i−1/2,j

(2.31)

= 1 2

[ 1+ αc

β ln

{cosh(β(1−x˜i)) cosh(βx˜_i)

}]

(2.32)

= 1 2

( 1+ αc

β lna₁

) (

a₁≡ cosh(β(1−x˜_i)) cosh(βx˜_i) ^と置く

)

(2.33) 従って，VOF^関数C_i は以下のように変形できる．

a₁ =exp ( β

αc(2C_i−1) )

(2.34) さらに，式(2.33)^のa₁^{の定義式について}

a₁ = cosh(β(1−x˜_i))

cosh(βx˜_i) (2.35)

= exp(β(1−x˜i))+exp(−β(1−x˜i))

exp(βx˜_i)+exp(−βx˜_i) (2.36)

= a₃²+a²₂

a3(a²₂+1) (2.37)

と変形できる．ただし，

a2≡exp(βx˜i) (2.38)

a₃≡exp(β) (2.39)

と置いた．式(2.37)^をa2 について整理して，a2を定義式に戻したうえでx˜i について整理しなお せば，

a₂²= a₃²−a₁a₃

a1a3−1 (2.40)

⇔x˜i = 1 2β ln

(a²₃−a1a3

a₁a₃−1 )

(2.41) となる．a1とa3は式(2.34), (2.39)から求めることができるので，これで未知数x˜iが求められる．

従って，(i+1/2,j)のセル界面を横切る流束F_x,i+1/2^{は次のように求まる．}

F_x,i+1/2=−

∫ _{xi+1/2−ui+1/2∆t}

x_i+1/2

χidx (2.42)

=−

∫ _{xi+1/2−ui+1/2∆t}

xi+1/2

1 2 [

1+αctanh

{β(x−x_i−1/2

∆x −x˜i

)}]

dx (2.43)

=−1 2 [

x+ αc∆x β ln

{ cosh

(β(x−xi−1/2

∆x

)−x˜_i

)}]^xi+1/2−ui+1/2∆t x_i+1/2

(2.44)

以上がx^方向速度u^がu ≥0^{の時の場合であった}. u <0の時も同様に考えることができ,^まとめた ものを以下に示す．

Fx,i+1/2=−

∫ _{xi+1/2−ui+1/2∆t}

xi+1/2

χiu pdx (2.45)

=−1 2

[

x+ αc∆x β ln

{ cosh

(β(x−x_{iu p−1/2}

∆x

)−x˜iu p

)}]^{xi+1/2−ui+1/2∆t}

xi+1/2

(2.46)

iup= {

i (u_i+1/2 ≥0)

i+1 (ui+1/2 <0) (2.47)

αc = {

1 (C_{iu p+1} ≥C_{iu p−1})

−1 (Ciu p+1 <Ciu p−1) (2.48)

˜

x_{iu p} = 1 2βln

(a²₃−a1a3

a₁a₃−1 )

(2.49)

a₁ ≡ cosh(β(1−x˜_i))

cosh(βx˜i) (2.50)

a3 ≡exp(β) (2.51)

WLIC (Weighted Line Interface Calculation)法

VOF法はセル中の液相の体積充填率を移流方程式に従って移流させるものであるが，そこで一つ 問題点が生じる．それは「界面の向きがわからない」ということである．VOF^{関数は体積充填率と} いうスカラー量であり，スタッガード格子ではセルの中央一点に定義されている．つまり，あるセル の中でVOF^{関数の値が}0.7であったとして，それだけでは自由表面がどのように分布しているか(^ど の方向に偏っているのか)を推定するのは容易ではない．もしかしたらセルの右側に0.7^{液相が寄っ} ているかもしれないし，セルの上側に0.3^・下側に0.4だけ液相が寄っているかもしれない．どちら の場合であってもVOF^関数値は0.7として換算される．このような問題点があるため，VOF^関数の 移流を考えるとき，周囲のVOF関数を基に界面の再構築を行う(界面の近似形状を求める)^必要があ る．この再構築の方法は幾つかあるが，代表的なものを以下に紹介する(Fig. 2.9^参照)^．

• SLIC (Simple Line Interface Calculation)

格子に対して平行線を引いた単純な補間関数を用いて計算を行う. コーディングが比較的容易 であるが,^{計算精度が悪い．}

• PLIC (Piecewise Linear Interface Calculation)

界面の勾配を考えた補間関数を用いて計算を行う. 様々な場合分けが必要なので，SLIC^に比 べ計算精度が高いがコーディングが複雑．

• WLIC (Weighted Line Interface Calculation)

界面を格子と平行な直線で表し,法線方向により重み付けをする. SLIC^{のようにシンプルであ}

るが, PLIC^と同様,界面方向を考慮しているので精度は高め.

Fig. 2.9:Schematic of SLIC, PLIC and WLIC method

本研究では，コーディングが容易でかつ比較的精度が高いWLIC法を採用した．以下にその詳細 を述べる．

WLIC法では界面を格子に平行な直線で表現をする. ^その際,補間関数を界面の単位法線ベクトル を用いて重み付けをする．

χi,j =ωxχx,i,j+ωyχy,i,j (2.52)

ωx = |n_x|

|nx|+|ny|, ωy = |n_y|

|nx|+|ny| (2.53)

ωx +ωy =1 (2.54)

単位法線ベクトルはVOF関数の勾配から計算すると精度が悪くなるので, 2^{次元計算の場合},^近傍 8点を用いて計算を行う．

n_x,i,j = (Ci+1,j+1+2Ci+1,j+Ci+1,j−1) − (Ci−1,j+1+2Ci−1,j+Ci−1,j−1)

8∆x (2.55)

ny,i,j = (Ci+1,j+1+2Ci,j+1+Ci−1,j+1) − (Ci+1,j−1+2Ci,j−1+Ci−1,j−1)

8∆y (2.56)

なお, 3^{次元計算の場合は近傍}26^{点を用いることで},^{以下のようになる．}

n_x,i,j = 1 32∆x

[{4C_i+1,j,k+2(C_i+1,j,k+1+C_i+1,j,k−1+C_i+1,j+1,k+C_i+1,j−1,k) +C_i+1,j+1,k+1+C_{i+1,j−1,k+1}+C_{i+1,j+1,k−1}+C_{i+1,j−1,k−1}}

−{

4Ci−1,j,k+2(Ci−1,j,k+1+Ci−1,j,k−1+Ci−1,j+1,k+Ci−1,j−1,k) + Ci−1,j+1,k+1+Ci−1,j−1,k+1+Ci−1,j+1,k−1+Ci−1,j−1,k−1

}]

(2.57)