2007 2

学位論文題目

Title

界面追跡法に基づく二相流数値予測手法に関する研究

氏名

Author

林, 公祐

専攻分野

Degree

博士（工学）

学位授与の日付

Date of Degree

2007-03-25

資源タイプ

Resource Type

Thesis or Dissertation / 学位論文

報告番号

Report Number

甲4024

URL

http://www.lib.kobe-u.ac.jp/handle_kernel/D1004024

※当コンテンツは神戸大学の学術成果です。無断複製・不正使用等を禁じます。

著作権法で認められている範囲内で、適切にご利用ください。

Create Date: 2016-10-08

界面追跡法に基づく

二相流数値予測手法に関する研究

2007

年

2

月

神戸大学大学院自然科学研究科

林 公祐

目次

第1章 序論 1 1.1 背景 . . . 1 1.2 二相流数値計算法 . . . 2 1.3 固定格子に基づく界面追跡法 . . . 4 1.3.1 体積追跡法 . . . 4 1.3.2 レベルセット法. . . 5 1.3.3 フェーズフィールドモデル . . . 10 1.3.4 フロント・トラッキング法 . . . 12 1.4 代表的な体積追跡法の性能評価 . . . 15 1.4.1 DA法 . . . 15 1.4.2 FCT . . . 16 1.4.3 PLIC型 . . . 18 1.4.3.1 FLAIR . . . 18 1.4.3.2 MARS . . . 19 1.4.4 CIP法. . . 21 1.4.5 既存手法の性能評価 . . . 22 1.5 体積追跡法における表面張力評価 . . . 26 1.5.1 高さ関数 . . . 26 1.5.2 CSFモデル . . . 27 1.5.3 CSSモデル . . . 29 1.5.4 最小二乗法を応用したモデル . . . 30 1.6 本研究の目的 . . . 31 1.7 本論文の構成 . . . 32 第1章の参考文献 . . . 34 第2章 界面追跡法NSSの開発 42

2.1 緒論 . . . 42 2.2 基礎方程式. . . 43 2.2.1 単相流の基礎方程式および二相界面における跳躍条件. . . 43 2.2.2 一流体近似による二相流の基礎方程式 . . . 44 2.2.3 一流体近似に基づく基礎方程式と跳躍条件の整合性 . . . 45 2.2.4 仮定・構成式 . . . 46 2.3 基礎方程式の数値解法 . . . 47 2.3.1 格子および変数配置 . . . 47 2.3.2 改良SOLA法 . . . 48 2.3.3 境界条件 . . . 50 2.3.4 運動方程式の離散化 . . . 52

2.4 Non-uniform subcell scheme (NSS)の開発. . . 56

2.4.1 体積率移流方程式 . . . 56 2.4.2 サブセルに基づく移流方程式の解法 . . . 57 2.4.3 体積率移流方程式の時間積分 . . . 61 2.5 表面張力評価法の提案 . . . 64 2.5.1 適切なカラー関数の関数形の考察 . . . 64 2.5.2 局所レベルセット関数による表面張力評価 . . . 69 2.5.3 秩序変数による表面張力評価 . . . 72 2.6 結論 . . . 74 第2章の参考文献 . . . 76 第3章 開発手法の性能評価 79 3.1 緒論 . . . 79 3.2 表面張力評価精度の評価 . . . 80 3.2.1 回転体に対する曲率評価 . . . 80 3.2.1.1 局所レベルセット関数算出時のパラメーター . . . 80 3.2.1.2 秩序変数算出時のパラメーター . . . 81 3.2.1.3 曲率ベクトル計算結果 . . . 81 3.2.2 静止液中単一液滴計算 . . . 83 3.3 基本特性の評価 . . . 87 3.3.1 体積保存およびシャープな界面の維持 . . . 87 3.3.1.1 ダクト内静止水中空気泡 . . . 87 3.3.1.2 水面への水滴落下 . . . 90 3.3.2 複数の界面が存在する系への適用性 . . . 90

3.4 気泡および液滴運動の予測 . . . 94 3.4.1 無限静止液中単一流体粒子 . . . 94 3.4.1.1 無限静止水中単一空気泡 . . . 94 3.4.1.2 無限静止液中単一空気泡 . . . 96 3.4.1.3 無限静止液中単一液滴 . . . 100 3.4.2 一様せん断流中気泡に働く揚力 . . . 102 3.4.2.1 揚力係数計算法 . . . 102 3.4.2.2 計算結果 . . . 103 3.5 低空間分解能計算の定性的妥当性 . . . 104 3.5.1 静止水中単一上昇気泡 . . . 106 3.5.2 せん断流中気泡に働く揚力 . . . 106 3.5.3 障害物周りを通過する気泡 . . . 109 3.6 結論 . . . 114 第3章の参考文献 . . . 115 第4章 曲線座標系への拡張 117 4.1 緒論 . . . 117 4.2 ３次元円柱座標系用NSS . . . 118 4.2.1 基礎方程式 . . . 118 4.2.2 基礎方程式の変数変換 . . . 119 4.2.3 改良SOLA法 . . . 121

4.2.4 Non-uniform subcell scheme . . . 124

4.2.5 レベルセット関数 . . . 127 4.3 ３次元一般曲線座標系用NSS . . . 128 4.3.1 一般曲線座標 . . . 128 4.3.2 境界適合座標 . . . 129 4.3.3 体積要素および面積要素ベクトル . . . 130 4.3.4 基礎方程式 . . . 132 4.3.5 流れ場の解法 . . . 133 4.3.5.1 T-CUP法 . . . 133 4.3.6 変数配置および近傍セルの表記法 . . . 133 4.3.7 圧力ポアッソン方程式と速度修正 . . . 135 4.3.7.1 (A)速度補間 . . . 136 4.3.7.2 (B)速度修正量補間 . . . 136 4.3.8 CIP法. . . 137

4.3.8.1 方向分割型CIP法. . . 137

4.3.9 Non-uniform subcell scheme . . . 140

4.3.10 レベルセット関数 . . . 143 4.3.11 表面張力評価 . . . 145 4.3.12 格子生成 . . . 146 4.3.13 ブロック分割法. . . 147 4.4 結論 . . . 148 第4章の参考文献 . . . 149 第5章 曲線座標系用NSSの性能評価 151 5.1 緒論 . . . 151 5.2 円管内液滴に対する壁効果に関する実験. . . 152 5.2.1 既存の修正係数モデル . . . 152 5.2.2 実験装置および実験方法 . . . 156 5.2.2.1 実験装置 . . . 156 5.2.2.2 液体の準備. . . 156 5.2.2.3 撮影の準備. . . 158 5.2.2.4 画像処理による液滴径および液滴速度の算出 . . . 159 5.2.3 実験結果 . . . 159 5.3 円管内静止液中単一液滴運動の数値予測. . . 163 5.3.1 計算条件 . . . 163 5.3.2 計算結果 . . . 164 5.3.3 数値計算を利用した抗力係数モデル作成の試み . . . 169 5.4 2 x 2ロッドバンドル内気泡の数値計算 . . . 172 5.4.1 計算条件 . . . 172 5.4.2 計算結果 . . . 173 5.5 円管内層流中空気泡の数値計算 . . . 179 5.6 結論 . . . 183 第5章の参考文献 . . . 185 第6章 結論 187 付録A ブレント法 191 A.1 ブレント法による非線形方程式の解法 . . . 191 付録Aの参考文献 . . . 192

付録B Thermo-CUP法 193 B.1 基礎方程式のフェイズ分割 . . . 193 B.1.1 移流フェイズ . . . 194 B.1.2 拡散フェイズ . . . 194 B.1.3 音響フェイズ . . . 196 B.2 方程式のまとめ . . . 197 付録Bの参考文献 . . . 197 付録C 圧力ポアッソン方程式の離散化 199 C.1 圧力ポアッソン方程式の離散化 . . . 199 C.2 速度修正 . . . 204 付録D 円柱座標系の基礎方程式 205 D.1 直交曲線座標系における勾配と発散 . . . 205 D.2 直交曲線座標系の基礎方程式 . . . 206 D.3 円柱座標系の基礎方程式 . . . 207 付録E 任意セル形状へのNSSの拡張 209 E.1 NSSによる任意形状セル内の界面位置算出法 . . . 209 付録Eの参考文献 . . . 213 付録F 乱流モデルの導入 214 F.1 フィルタ操作を施した一流体近似に基づく基礎方程式 . . . 214 F.2 サブグリッドスケールモデルの組み込み. . . 216 付録Fの参考文献 . . . 217

図一覧

1.3.1 DA法の概念図 . . . 6 1.3.2 PLICの概念図 . . . 6 1.3.3 レベルセット関数 . . . 7 1.3.4 再初期化の効果 . . . 8 1.3.5 パーティクルレベルセット法 . . . 10 1.3.6 粒子位置情報によるレベルセット関数の修正 . . . 10 1.3.7 フロント・トラッキング法における界面表現 . . . 12 1.3.8 フロント・トラッキング法における表面張力評価 . . . 13 1.4.1 DA法におけるセルの記法 . . . 16 1.4.2 界面形状の分類 . . . 18 1.4.3 サブケース. . . 19 1.4.4 サブケース判別ダイアグラム . . . 19 1.4.5 傾斜線分による界面の近似 . . . 20 1.4.6 計算領域 . . . 23 1.4.7 円形状輸送結果(d/∆x = 20, u∆t/∆x = 0.10) . . . . 24 1.4.8 円形状輸送結果(d/∆x = 10, u∆t/∆x = 0.25) . . . . 24 1.4.9 球形状輸送結果(d/∆x = 20, u∆t/∆x = 0.10) . . . . 25 1.4.10 球形状輸送結果(d/∆x = 20, u∆t/∆x = 0.25) . . . . 25 1.5.1 高さ関数H . . . 26 1.5.2 高精度高さ関数計算用ステンシル . . . 28 1.5.3 界面を近似した多角形 . . . 31 2.2.1 界面に隔てられた不浸透の二流体 . . . 44 2.2.2 界面を含み，各相側に突出した計算セル. . . 46 2.3.1 計算格子 . . . 48 2.3.2 スタガード配置 . . . 48

2.3.3 境界セル . . . 51 2.3.4 CIP法の概念 . . . 53 2.4.1 界面セル . . . 58 2.4.2 サブセル配置 . . . 58 2.4.3 距離ベクトル d(v) . . . 59 2.4.4 界面がセルから切り取る体積 . . . 60 2.4.5 再構築された界面および被輸送体積 . . . 60 2.4.6 EI(Eulerian implicit) . . . 63 2.4.7 LE(Lagrangian explicit) . . . 64 2.5.1 曲率ベクトル評価制度検証用回転体 . . . 67 2.5.2 κnの真値 . . . 67 2.5.3 体積率αを評価法Aに適用して得たκn . . . 67 2.5.4 体積率αを評価法Bに適用して得たκn . . . 68 2.5.5 レベルセット関数φを評価法Aに適用して得たκn . . . 68 2.5.6 レベルセット関数φを評価法Bに適用して得たκn . . . 68 2.5.7 局所レベルセット関数φの算出 . . . 70 2.5.8 秩序変数ϕの関数形 . . . 73 3.2.1 曲率ベクトル評価精度検証用回転体 . . . 80 3.2.2 κnの真値 . . . 80 3.2.3 局所レベルセット関数φを評価法Aに適用して得たκn . . . 82 3.2.4 局所レベルセット関数φを評価法Bに適用して得たκn . . . 82 3.2.5 再初期化を施したφを評価法Aに適用して得たκn . . . 82 3.2.6 再初期化を施したφを評価法Bに適用して得たκn . . . 83 3.2.7 秩序変数ϕを評価法Aに適用して得たκn . . . 83 3.2.8 秩序変数ϕを評価法Bに適用して得たκn . . . 83 3.2.9 体積率αを評価法Aに適用して得たκn . . . 84 3.2.10 体積率αを評価法Bに適用して得たκn . . . 84 3.2.11 擬似流れ評価用計算体系 . . . 85 3.2.12 体積率αを表面張力評価に使用して得た液滴周囲の流れ場 . . . 85 3.2.13 局所レベルセット関数φを表面張力評価に使用して得た液滴周囲の流れ場 86 3.2.14 秩序変数ϕを表面張力評価に使用して得た液滴周囲の流れ場 . . . 86 3.3.1 ダクト内気泡計算用計算体系 . . . 88 3.3.2 ダクト内静止水中空気泡の形状(a)および気泡周りの瞬時の速度場(b) . . 89 3.3.3 気泡周りの瞬時の圧力場(a)および気泡端部の体積率分布(b) . . . 89

3.3.4 気泡体積θの時間変化 . . . 90 3.3.5 水滴計算用計算体系 . . . 91 3.3.6 水面への水滴落下 . . . 91 3.3.7 液滴周囲の体積率分布 . . . 92 3.3.8 気泡プルームの計算 . . . 93 3.4.1 無限静止水中空気泡の終端上昇速度 . . . 95 3.4.2 無限静止水中単一空気泡計算用計算体系. . . 95 3.4.3 無限静止液中単一空気泡計算用計算領域. . . 97 3.4.4 無限静止液中単一空気泡の上昇速度 . . . 98 3.4.5 気泡形状および気泡周りの流線(いずれも左半分はHnat [7]らの実験写真) 99 3.4.6 無限静止液中単一液滴の液滴レイノルズ数(実験結果はMyintらによる[8])101 3.4.7 液滴形状の計算結果と実験結果[8] . . . 101 3.4.8 一様せん断流中空気泡計算用計算体系 . . . 102 3.4.9 一様せん断流中空気泡の横方向運動 . . . 104 3.4.10 一様せん断流中空気泡の揚力係数 . . . 105 3.5.1 静止水中単一気泡 . . . 107 3.5.2 気泡終端上昇速度の格子依存性 . . . 107 3.5.3 一様せん断流中単一気泡 . . . 108 3.5.4 せん断流中を上昇する気泡の横方向運動. . . 109 3.5.5 揚力係数CL の予測値 . . . 110 3.5.6 計算領域 . . . 110 3.5.7 せん断流中気泡の運動 . . . 112 3.5.8 気泡周りの体積率分布 . . . 113 4.2.1 スタガード配置 . . . 120 4.2.2 界面セルに準備したサブセル . . . 125 4.2.3 距離ベクトル d(v) . . . 126 4.3.1 曲線座標 . . . 129 4.3.2 物理空間(a)と計算空間(b) . . . 130 4.3.3 物理空間における体積要素 . . . 131 4.3.4 面積要素ベクトルdS . . . 132 4.3.5 コロケート配置 . . . 134 4.3.6 隣接セル . . . 135 4.3.7 微係数対応. . . 138 4.3.8 計算空間における界面セル(a)とサブセル配置(b) . . . 141

4.3.9 距離ベクトルξi_(v) . . . 143 4.3.10 ブロック分割法 . . . 147 4.3.11 仮想セル . . . 148 5.2.1 固体粒子に対する壁効果 . . . 153 5.2.2 固体粒子に対する修正係数 . . . 154 5.2.3 流体粒子に対する壁効果 . . . 155 5.2.4 気泡に対する修正係数 . . . 155 5.2.5 実験装置の模式図 . . . 157 5.2.6 液滴速度の時間変化 . . . 158 5.2.7 画像の屈折補正(左：原画像，右：屈折補正後) . . . 160 5.2.8 画素番号の記法 . . . 161 5.2.9 液滴体積計算 . . . 161 5.2.10 円管内液中単一液滴の終端上昇速度 . . . 162 5.2.11 液滴レイノルズ数 . . . 162 5.2.12 修正係数の実験値 . . . 163 5.2.13 Habermanの修正係数モデル . . . 164 5.3.1 円管内静止液中液適用計算体系 . . . 165 5.3.2 界面近傍の体積率分布 . . . 165 5.3.3 液滴形状(C1:log M = −3.0, KF96-30) . . . 166 5.3.4 液滴形状(C2:log M = −3.0, KF96-100) . . . 167 5.3.5 液滴形状(C3:log M = −4.8, KF96-30) . . . 167 5.3.6 液滴内部および周囲の速度場 . . . 168 5.3.7 修正係数の比較 . . . 168 5.3.8 ２次元円柱座標系用NSSによる計算結果 . . . 169 5.3.9 液滴レイノルズ数の数値予測結果と作成した抗力係数モデル . . . 170 5.3.10 提案した抗力係数モデル式による液滴終端速度 . . . 171 5.3.11 提案した抗力係数モデル式に基づく修正係数 . . . 172 5.4.1 2 × 2ロッドバンドル . . . 174 5.4.2 計算格子の例 . . . 174 5.4.3 気泡界面周囲の体積率分布(d=10 mm) . . . 175 5.4.4 気泡終端上昇速度 . . . 176 5.4.5 気泡形状の比較 . . . 177 5.4.6 蒸気泡の予測形状 . . . 178 5.5.1 円管内気泡計算用計算格子 . . . 179

5.5.2 気泡体積θの時間変化 . . . 180 5.5.3 気泡形状と気泡周りの速度場の時間変化. . . 181 5.5.4 気泡重心の軌跡 . . . 182 5.5.5 壁面への気泡の衝突による渦の形成 . . . 182 E.1.1 三角形セルに準備された三角形サブセル. . . 210 E.1.2 距離ベクトル dv . . . 211 E.1.3 与えられた３角形セルの体積率αに対して算出した距離関数φint . . . 211 E.1.4 五角形セルに準備された三角形サブセル. . . 212 E.1.5 与えられた五角形セルの体積率αに対して算出した距離関数φint . . . 212 F.1.1 界面を含む検査体積 . . . 215

表一覧

1.2.1 二相流計算手法[10] . . . 3 1.3.1 界面追跡法の特性比較[44] . . . 14 3.3.1 水‐空気系の物性値 . . . 88 3.4.1 液‐液系の無次元数 . . . 100 3.5.1 水および蒸気の物性値(圧力7MPa,温度559K) . . . 106 3.5.2 せん断流中気泡計算条件 . . . 108 5.2.1 静止液中液滴の実験条件 . . . 158 5.3.1 円管内液中液滴の実験条件 . . . 170 5.4.1 水‐空気系の物性値 . . . 173 5.4.2 水‐蒸気の物性値(圧力7MPa,温度559K) . . . 173

主な使用記号

[英文字] a, b, R, T ファンデアワールス流体モデルのパラメーター C 平滑化カラー関数 C0 カラー関数 CD 抗力係数 CL 揚力係数 d 球等価直径[m] D 管直径[m] dv セル中心からサブセル頂点vへの距離ベクトル dS 面積要素[m2] dS 面積要素ベクトル[m2] dΩ 体積要素[m3] E アスペクト比 Eo エトベス数 f 修正係数 f 体積力ベクトル[N/kg] F 自由エネルギー汎関数 Fσ 表面張力ベクトル g 共変計量テンソルの行列式 g 重力加速度ベクトル g1, g2, g3 重力加速度ベクトルのデカルト成分 g1, g2, g3 共変基底ベクトル g1, g2, g3 反変基底ベクトル gi j 共変計量テンソル gi j 反変計量テンソル h セル幅[m]

Hε 平滑化ヘビサイド関数 I 単位テンソル J 座標変換のヤコビアン(= |∂xi/∂ξj|) M モルトン数 n 界面法線ベクトル N1, N2, N3 界面法線の曲線座標ξi における共変成分 Kε 界面の厚みを規定するパラメーター p 圧力[Pa] qr, qθ, qz 速度ベクトルの円柱座標成分(= rur, uθ, uz) r, θ, z 円柱座標 r 位置ベクトル Re レイノルズ数 sign_ε 平滑化sign関数 t 時間[s] T 応力テンソル[Pa s] Tσ キャピラリーテンソル[J/m3] u 速度ベクトル[m/s] ur, uθ, uz 速度ベクトルの円柱座標物理成分[m/s] U1_{, U}2_{, U}3 _{曲線座標ξ}i _{における反変速度成分} VT 終端速度[m/s] w 再初期化方程式の移流速度(= sign_ε(φ0)(∇φ/|∇φ|)) W 重み関数 x1, x2, x3, x1, x2, x3 デカルト座標 xint セル中心から界面への距離ベクトル [ギリシャ文字] α セル平均体積率 αsub _{サブセル体積率} Γ オンサガー係数 δS デルタ関数 δε 平滑化デルタ関数 ∆t 時間刻み幅[s] ∆x1, ∆x2, ∆x3 物理空間におけるセル幅[m] ∆ξ1_{, ∆ξ}2_{, ∆ξ}3 _{計算空間におけるセル幅}

ε ヘビサイド関数の平滑化幅を規定するパラメーター η 化学ポテンシャル µ 粘度[Pa s] µ∗ _粘度比 κ 曲率[1/m] λ 直径比(= d/D) ξ1_{, ξ}2_{, ξ}3 _{一般曲線座標} ξi (v) セル中心からサブセル頂点vへの距離ベクトル ρ 密度[kg/m3_] σ 表面張力[N/m] τ 時間tとは無関係な擬似時間 φ 局所レベルセット関数，レベルセット関数 φ0 再初期化開始前のレベルセット関数 φint 界面セルの局所レベルセット関数 Φ 速度修正量のポテンシャル ϕ 秩序変数 χ 相定義関数 Ω セル体積，被積分体積[m3] Ωsub _{サブセル体積} [添字] D 液滴 G 気相 int 界面 L 液相 n 相n p 相 p P 粒子，流体粒子

第

1

章

序論

1.1

背景

工業機器に見られる流れは二相流であることが少なくない．二相流が関連する工業機器 を安全かつ高効率に設計するためには，二相流の理解が必須である．単相流の場合には， 支配方程式を解く際，系の内部および外部境界における境界条件を考慮する [1]．一方， 二相流では相界面で成立する条件，すなわち跳躍条件[2]も考慮する必要がある．しかも， 界面の位置および形状は一般に非定常である．このことは，二相流の理解を困難にする主 因である．例えば，１未満の低粒子レイノルズ数(慣性力と粘性力の比)で流体中を運動す る単一固体球や単一球形流体粒子のように，非常に単純な系に対しては理論解が得られて いる[3, 4]が，実際の工業機器における粒子レイノルズ数は１よりもはるかに大きく（水 中を0.2 m/sで上昇する直径8 mmの気泡のようなありふれた条件でも気泡レイノルズ数 は1600程度になる[5]），かつ粒子数も膨大になる．そのような状況においては，理論解 を得ることは期待できない．このため，1910年頃のボイラー水循環系の安全設計を目的 とした研究に端を発する気液二相流研究は，主として実験的方法で行われた[6]．原子力 発電が登場してからは，原子力プラントの熱水力設計に高度の安全性が要求されるため， プラント内気液二相流に関する膨大な実験データが収集された[7]．一方，基礎方程式を 数値的に解く数値流体力学[8]が発展し，二相流数値予測のために種々の計算モデルが考 案された[7, 9, 10]．初期に考案された均質流モデル[9]やドリフトフラックスモデル[11] は，装置の大域的な物理量の評価に対して有力である．さらに，より厳密な基礎方程式に 基づく二流体モデル[9]が提示されたことで，より多くの情報が得られるようになった． 最近では，計算機性能の向上によって，個々の気泡を直接追跡する界面追跡法[2, 12]や 気泡追跡法[2]による詳細計算が可能になってきた．

1.2

二相流数値計算法

個々の気泡の挙動を追跡するのではなく，気相の存在割合を表わすボイド率や速度，圧 力，温度などの平均量を従属変数として定式化する方法を，平均化二相流モデルという． 平均化二相流モデルは，混合流モデルおよび二流体モデルに分類される．混合流モデルは さらに均質流モデル，スリップ比モデル[13]，ドリフトフラックスモデルに分類される． 均質流モデルでは，気相と液相が同一速度で流れ，かつ均質な混合流であるとみなす．こ の仮定のもとで導かれた基礎方程式は単相流のものと同様な形になるが，密度などの変数 は各相の存在割合によって重み付けされた形で定義される．均質流モデルは，相間速度差 が小さな高液流速の気泡流などには適用可能であるが，相間速度差が大きい流れには適用 できない．スリップ比モデルは，二相の速度比を考慮して均質流モデルを拡張した方法で ある．スリップ比モデルでは，例えば液相が静止している気泡流ではスリップ比が無限大 となることや，その物理的意味が曖昧なことが問題であった．そこで，相間相対速度を導 入したドリフトフラックスモデルが Zuberによって提案された．これらの混合流モデル では，基礎方程式を閉じるために実験相関式が必要である． 混合流モデルでは，均質の仮定やドリフトフラックスパラメーターと呼ばれる実験定数 の導入により，基礎方程式の数を単相流の場合と同数にしている．一方，多流体モデルで は，各相に保存則を適用して，各相で成立する基礎方程式を導く．相間の保存量輸送は， 基礎方程式に含まれる相間相互作用の項が担う．本モデルは，多くの大胆な仮定に基づく 混合流モデルに比べてより厳密である．現在では，３次元多流体モデルによる二相流計算 が広く行われている．ただし，相間相互作用の評価に必要となる関係式の数は，混合流モ デルよりも多い．また，基礎方程式が数学的に不適切であることが指摘されている[14]． 気泡追跡法[15–21]とは，個々の気泡をラグランジュ的に追跡し，液相流れはモデル あるいは基礎方程式の解として与える方法である．液相流れをモデルで与える方法を one-way，基礎方程式を解いて求める方法をtwo-way気泡追跡法と呼ぶ．各気泡をラグラ ンジュ的に追跡するという手法の性格上，その適用先は気泡流やスラグ流などの流動様式 に限定されている．３次元気泡追跡法は，Tomiyama らによって初めて提案された[18]． ´

Zunらも最近３次元 one-way気泡追跡法を開発している[21]．また，two-way気泡追跡 法において，計算セルサイズよりも大きい楕円，キャップ，テイラー気泡を扱えるのは Tomiyamaらの方法[19]のみである．two-way気泡追跡法では抗力，揚力などの相間相互 作用の相関式が，one-way気泡追跡法ではさらに液相流れに関する理論あるいは実験モデ ルが必要である． 界面追跡法とは，界面を記述する変数の時間発展方程式を解いて界面を追跡する方法の 総称である．界面を記述する方法によって以下の方法に大別される：(1)格子そのものが

表1.2.1 二相流計算手法[10] d: bubble diameter, ∆x: cell size

Spatial Method Fundamental Eqs. Applicability to

Resolution Practical Problem

low Averaging Method Averaged Field Eqs. + high

(∆x >> d) (Homogeneous, Drift-Flux, Two-Fluid Models) Constitutive Eqs.

intermediate Bubble Tracking Method Eq. of Bubble Motion + intermediate (∆x ∼ d) (One-Way, Two-Way Methods) Constitutive Eqs.

high Interface Tracking Method Navier-Stokes Eq. + low (∆x < d) (Front Tracking, Volume Tracking, Level Set) Surface Tension

high Microscopic Method Translation & Collision low (∆x << d) (Lattice Gas, Lattice Boltzmann Methods) of Pseudo Molecules

界面を表わす方法，(2)固定格子の各格子点上に定義したスカラー変数により界面を記述す る方法，(3)流れ場は固定格子上で解くが，界面を記述するために別の非構造格子を準備す る方法．代表的なものに，(1)BFC (boundary-fitted coordinate)法[22–24], ALE (arbitrary Lagrangian-Eulerian) 法[25–27], (2) 体積追跡法 (volume tracking method) [12, 28], レベ ルセット法 [29],フェーズフィールドモデル (Phase-Field Model) [30–33], (3)フロント・ トラッキング法[34, 35]がある．また，仮想粒子の並進および衝突を各格子点上で行い マクロな流れ場を予測するLBM(lattice Boltzmann method) [36]を，二相流に拡張して界 面追跡計算を行っている研究例 [37, 38]や，多数の粒子により流体運動を計算する粒子 法[39, 40]により界面追跡計算を行っている研究例もある． 平均化二相流モデルおよび気泡追跡法では実験相関式が必須であり，したがって計算精 度および手法の適用範囲は相関式に依存する．一方，界面追跡法は実験相関式を必要とし ないため，相関式が整備されていない二相流にも適用可能である．界面追跡法は，平均化 二相流モデルや気泡追跡法に比べると計算負荷が大きいが，計算機性能の発達に伴いその 適用範囲は拡大している．今後より実際的な問題，例えば円管内や原子炉燃料集合体内二 相流などに適用するには，界面追跡法の高精度化に加えて，円柱座標系や一般曲線座標系 用への拡張が必要となる．Tomiyamaがまとめた各計算法の特徴[10]を表1.2.1に示す． 様々な空間スケールを有する分散性二相流を，多流体モデル，気泡追跡法あるいは界面 追跡法単独で計算するのは困難である．例えば，界面追跡法は計算セルサイズよりも十分 に大きい気泡は扱えるが，計算セルよりも小さい気泡は扱えない．逆に，多流体モデルは 計算セルより小さい多数の気泡を扱えるが，複数の計算セルにわたる大気泡は扱えない． この課題を克服するために，多流体モデル，気泡追跡法，界面追跡法を統合する研究が行 われている．多流体モデルと体積追跡法を結合した先駆的な方法として，Tomiyamaらの 提案した(N+2)-fieldモデル[41]がある．彼らは，界面追跡法により数セル以上の大きさ

になる大気泡や自由表面を計算し，セルよりも小さい気泡を多流体モデルにより計算する ことで，従来法単独では計算できなかった大小気泡が混在する気泡プルームや気泡塔を模 擬した流れの計算を行っている．最近では，化学反応[42]や，気泡同士の合体・分裂[43] を考慮できるように，(N+2)-fieldモデルの拡張が進められている．また，Tomiyamaら は，(N+2)-fieldモデルの概念を発展させて，多流体モデル，体積追跡法およびtwo-way 気泡追跡法を統合した二相流計算手法を提案した[44, 45]．計算者は，各手法を組み合わ せることにより，様々な空間スケールを有する二相流に対応することができる．基礎方程 式はよく整理されており，各手法が取り扱う体積率の有効・無効を選択することで，各手 法単独，あるいは複数手法を組み合わせた計算を容易に実行できる．本手法により，これ まで単独の手法では計算が困難であった多くの二相流が予測可能になると期待される．こ れら統合手法では，界面追跡機能として体積追跡法が採用されている．その理由として， 体積追跡法は界面を体積率で記述するため，体積率を扱う多流体モデルおよび気泡追跡法 と結合しやすいことが考えられる．ただし，統合計算法によって信頼性の高い予測を得る ためには，各々の計算手法が高精度でなければならない．すなわち，精度よく界面運動を 予測できる体積追跡法の開発が必須である．

1.3

固定格子に基づく界面追跡法

体積追跡法に代表される固定格子に基づく界面追跡法は，格子そのものが変形してラグ ランジュ的に界面を追跡するBFC法やALE法に比べて，複数の界面を含む流れに適用し やすく，計算負荷も小さい．このため，工業的に興味のある二相流に対する適用性は，ラ グランジュ的方法よりも固定格子に基づく方法の方が高いと考えられる．本節では，代表 的な固定格子に基づく界面追跡法を説明する．

1.3.1

体積追跡法

体積追跡法の研究は，DeBar [46] やNohら[47]が1970年代に報告した先駆的方法に 始まる．現在広く知られているのは，HirtらによるVOF(volume of fluid) 法[28]であろ う．体積追跡法では，セル平均体積率 α(計算セルを流体１が占めるときα = 1,流体２が 占めるときα = 0,セル内に界面が存在する場合0 < α < 1)により界面を記述する．次の 体積率移流方程式を解き，界面を追跡する． ∂α ∂t + u · ∇α = 0 (1.3.1) ここで，t は時間，u は速度である．体積追跡法により精度よい二相流予測を行うには， 上式を精度よく解く必要がある．ただし，通常の差分法を適用すると数値拡散が生じて界

面がぼやけるという問題が生じる．そこで，各セル内での界面形状や界面位置を考慮す る方法が適用される．例えば，VOF法ではDA(donor-acceptor)法[48]を採用している． DA法では，図1.3.1に示すように，勾配のない線分で界面を近似表現する．この界面形 状と，セル表面に定義した速度成分uおよび時間刻み幅∆tに基づいて，隣接セルに輸送 される流体量を算出する．DA法では，計算の進行に伴い，体積率の数値拡散および流体 体積の誤差が生じる．そこで，図1.3.2に示すように，界面を表現する線分の傾きも考慮 することで計算精度を改善する PLIC(piecewise linear interface calculation)[49]が提案さ れた．現在では，いくつかのPLIC型解法[49–52]が提案されている．しかし，それらの 多くは複雑な幾何計算を要するため，多次元化が容易ではない．次節で代表的な体積率移 流方程式の解法の詳細と性能評価について述べる．

また，表面張力評価精度が低いことは，体積追跡法の本質的な問題である[53]．体積追 跡法では，CSF(continuum surface force)モデル[54]を用いて，表面張力を体積力として 運動方程式に組み込むことが多い．表面張力Fσは次式で与えられる． Fσ= σκnδS (1.3.2) ここで，σは表面張力，κは界面曲率，nは界面法線，δS はデルタ関数である．上式から 明らかなように，表面張力を良好に評価するには界面曲率κを精度よく求めなければなら ない．κは次式で与えられる． κ = 1 |∇C| " ∇C |∇C| · ∇ |∇C| − ∇ · ∇C # (1.3.3) Cは一定の厚みを有する界面を記述するカラー関数である．体積追跡法では，Cに体積率 αを採用する．しかし，CSFモデルとαの組み合わせでは，精度良く表面張力を評価でき ないことが指摘されている[53]．

1.3.2

レベルセット法

距離関数により界面を記述するレベルセット法の研究は，Osher & Sethian によって始

められた[55, 56]．位置xにおける界面からの距離 Lに，正負の符号を付けた以下の符号 付距離関数φを定義する(図1.3.3)． φ (x, t) = ±L (1.3.4) 界面はφ = 0の等値面として記述される．φの移流方程式は次式で与えられる． ∂φ ∂t + u · ∇φ = 0 (1.3.5)

u

u

∆

t

i

i+1

transferred volume

interface represented

by line segments

図1.3.1 DA法の概念図

u

u

∆

t

i

i+1

interface represented

by line segments

図1.3.2 PLICの概念図 固定格子上で連続の式および運動方程式を解き，得られた速度場で上式を解いて次時刻の 界面位置を求める．移流項には ENO(essentially non-oscillatory) スキーム[57]が適用さ れることが多い． 各格子点がどちらの相に属するかはφの正負から判断されるが，密度や粘度などの物性 は数値的安定化のために界面近傍で平滑化される．同様に，表面張力も界面近傍の有限幅

φ

(t) = 0

on interface

φ

(t+

∆

t) = 0

on interface

φ

(t) = +L’’

φ

(t) =

−

L’

L’

L’’

図1.3.3 レベルセット関数 領域に連続的に作用する．レベルセット法における表面張力は次式で評価される． Fσ= σκnδε = σκ∇Hε (1.3.6) ここで，δεおよびHε は各々平滑化されたデルタ関数およびヘビサイド関数である．この 表面張力評価は体積追跡法におけるCSFモデルと一致する[58]．ただし，φによる曲率 評価精度は，体積率αによる評価に比べて著しく良いことが知られている．第2章で述べ るが，この理由は φの関数形が曲率を評価するカラー関数としてαよりも適切なためで ある． 移流方程式を解いてφを移流させると，本来距離関数が有する性質（空間勾配の大きさ が１）が失われる(図1.3.4(a))．Sussmanら[29]は，次の再初期化方程式を解くことによ り，本来の距離関数としての性質を回復させる方法を提案した． ∂φ ∂τ = signε(φ0) (1 − |∇φ|) (1.3.7) ここで，τは再初期化の時間，sign_ε は平滑化された sign関数，φ0 は再初期化開始時の φである．上式を用いた計算例を図1.3.4(b)に示す．再初期化をしない図1.3.4(a)と比べ て，距離関数としての性質が時間進行後も維持されていることが確認できる．しかし，再 初期化を実行する場合，以下の２つの問題が生じる．まず，再初期化方程式の時間進行は 計算負荷を増大させる．次に，再初期化に伴って界面位置が動いてしまい，流体体積が保 存しない．体積を保存させながら再初期化を行うために，いくつかの方法が提案されてい る．Changら[58]は，通常の再初期化方程式に加えて，初期体積に戻すために次式を解 く方法を提案した． ∂φ ∂τ + (θ (t = 0) − θ (τ)) (−P + κ) |∇φ| = 0 (1.3.8) ここで，θ(t = 0)は計算開始時の流体体積，θ(τ)は時間 τにおける流体体積，Pおよびκ は適当な値をとるパラメーターである．Meierらの方法[59]は，Changらの方法を簡略

0 20 40 60 -10 0 10 20 x φ cycle=0 cycle=100 cycle=200 0 20 40 60 -10 0 10 20 x φ cycle=0 cycle=100 cycle=200 0 cycle 100 cycle 200 cycle 0 cycle 100 cycle 200 cycle

without reinitialization with reinitialization

(a) (b) 図1.3.4 再初期化の効果 化したものと解釈できる． φ = φ0+ S−1[θ (t = 0) − θ (τ)] (1.3.9) ここで，S は界面積である．これらの方法では，体積誤差がどこで生じたのかは考慮せず に全体の形状を修正するため，移流方程式の予測した界面形状とは異なった形状になって しまう． 一方，Sussmanら[60]は局所的に体積を保存させる方法を提案した．各計算セル内に

含まれる流体体積は，平滑化ヘビサイド関数Hε(φ)により表される． Z Ω Hε(φ0) dΩ 再初期化の時間進行中に上式で定義された体積が変化しない条件， ∂ ∂τ Z Ω Hε(φ0) dΩ = 0 (1.3.10) を，再初期化方程式に束縛条件として組み込む． ∂φ ∂τ = signε(φ0) (1 − |∇φ|) + λH 0 ε(φ0) |∇φ0| (1.3.11) ここで，λはラグランジュ未定乗数である．本手法は各計算セルの流体体積の保存を保証 するので，再初期化が大域的界面形状に及ぼす影響は小さいと考えられる．しかしなが ら，通常の再初期化のみの場合と比べて，体積誤差は約半分に減少する程度の改善しか見 られない．Takahiraら[61]は，Sussmanらの方法によるφの修正量を，経験的に決めた 定数を掛けて補正する方法を試みているが，適切な定数の値を決定するために何度も計算 を繰り返す必要がある． Wheeler [62]は，通常の再初期化方程式において，界面からセル幅程度の距離h内にあ る計算セルに対しては，条件， sign_ε(φ0) = 0 (if φ0 < h) (1.3.12) を課すのみという単純な方法を提案している．Sussmanの方法とWheelerの方法では，再 初期化中に体積変動がないように工夫しているが，移流方程式の解が含む体積誤差は時間 進行とともに蓄積していく． 上述の体積誤差対策は再初期化方程式に基づいている．一方，再初期化の代わりに， 流れに追従する粒子を利用して体積誤差を抑えるパーティクルレベルセット法が最近 Enrightらにより提案された[63]．粒子は流れ場に影響を及ぼさず，レベルセット関数の 補正にのみ用いられる．粒子は図1.3.5のように界面近傍に配置される．このとき，いろ いろな大きさの粒子を準備する．また，粒子の重なりを許すことで，粒子群による界面表 現精度を向上させる．移流前にφ > 0の領域にある粒子を正粒子，φ < 0の領域にある粒 子を負粒子とラベリングする．φおよび全粒子を移流した後，φと粒子の符合が一致しな い場合，φ(界面)の移流に誤差が生じたと判断する．φの修正は近傍粒子の情報に基づい て行われる(図1.3.6)．粒子位置の時間変化は３次精度TVDスキームおよびルンゲ・クッ タ法を用いて高精度に解かれる．この方法によれば再初期化方程式を修正する従来法に比 べて良好に体積を保存できることが報告されているが，アルゴリズムの複雑化，粒子数の 増加に伴う計算時間の増加が課題である．

r

_p

r

_p

r

p

r

p

interface

p

₁

p

₂

p

₃

p

4 図1.3.5 パーティクルレベルセット法

r

_p

r

_p

r

_p

φ

< 0

φ

> 0

escaped negative particle

correct interface

interface after

reinitialization

positive particles

図1.3.6 粒子位置情報によるレベルセット関数の修正

1.3.3

フェーズフィールドモデル

二相系の時間発展は，ギンツブルグ‐ランダウ型自由エネルギー汎関数F [64]， F = Z Ω  f (c) + Kε 2 |∇c| 2  dΩ (1.3.13) が最小になる方向に進む．ここで，cは濃度，f は自由エネルギー，Kεは界面厚さを既定 するパラメーターである．上式右辺第２項は界面自由エネルギーに相当する．濃度に関す る保存則， ∂c ∂t = ∇ · Dcµ (1.3.14) における化学ポテンシャルµの関数形はF の変分を調べることで得られる． ∂c ∂t = ∇ · " Dc∇ ∂ f ∂c + Kε∇ · ∇c !# (1.3.15)

ここで，Dcは拡散係数である．上式は保存系時間依存型ギンツブルグ‐ランダウ方程式 あるいはカーン‐ヒリヤード方程式 [65]と呼ばれる．cの空間分布が界面形状を記述す る．本手法では，自由エネルギーが最小となる方向へ向かって自発的に界面構造が決まる ため，体積追跡法において必要とされるような，界面形状の再構築などは不要である．マ クロな界面の移流計算を行うために，カーン‐ヒリヤード方程式に移流項を追加する．高 田ら[32, 33]は，従属変数を相の秩序状態を表わす秩序変数ϕに置き換え，移流項を付加 した次式を用いて界面追跡計算を行っている． ∂ϕ ∂t + ∇ · uϕ = ∇ · " Γ (ϕ) ∇ ∂ f ∂ϕ + Kε∇ · ∇ϕ !# (1.3.16) ここで，Γはオンサガー係数[66]である．速度および圧力場は，連続の式および表面張力 を考慮したナビエ‐ストークス方程式を解くことにより求める． ∇ · u = 0 (1.3.17) ∂ρu ∂t + ∇ · ρuu = −∇p + ∇ · τ + ρKW∇ (∇ · ∇ρ) (1.3.18) ここで，pは熱力学的圧力，τは粘性応力テンソル，KW は界面自由エネルギーと界面密 度勾配に関係するパラメーターである．表面張力項は，系の自由エネルギー汎関数と質量 保存則から導かれる．カーン‐ヒリヤード方程式とナビエ‐ストークス方程式を組み合わ せてマクロな界面を追跡する本手法をフェーズフィールドモデルと呼ぶ．なお，本モデル では Kεにより規定される界面厚さを有するが，界面厚さをゼロとする極限操作により得 られる境界条件は，界面を不連続面とした場合に導かれる跳躍条件と等価である[67]． カーン‐ヒリヤード方程式と，ナビエ‐ストークス方程式の圧力および表面張力項は， 自由エネルギーが最小になるように系が進展するという概念から導かれており，熱力学的 に一貫性がある．これは，体積追跡法,レベルセット法,フロント・トラッキング法などの 方法にはない長所である．Jametら[68]は，この長所を活かして，表面張力評価誤差に起 因する擬似流れを消去できる方法を提案している． フェーズフィールドモデルは，他の界面追跡法に比べて新しく，実際の計算例はまだ少 ないものの，相変化 [69, 70]や濡れ性[71]を考慮した計算が実施されるなど，界面追跡 法としての有用性が確認されてきている．一方，本来メゾスケールの二相混合・分離など を扱うカーン‐ヒリヤード方程式を，マクロな流れを予測するナビエ‐ストークス方程式 と結合する点には疑問が残る．カーン‐ヒリヤード方程式右辺最終項のϕのラプラシア ンは界面の曲がり具合を表わす．この項にはさらに２階偏微分がかかっているが，これは 界面の曲がりを平滑化する（界面自由エネルギーを減少させる），すなわち表面張力とし て働く．しかし，ナビエ‐ストークス方程式にも表面張力は含まれているため，表面張力 が２重カウントされている可能性がある．また，多くの場合，ナビエ‐ストークス方程式

x

_p

x

₁

x

₂

x

₃

u

_p

x

₁

x

₂

x

_p

u

_p (a) (b) 図1.3.7 フロント・トラッキング法における界面表現 の表面張力項を適当な形に変形して，その一部を圧力項に吸収させる．この場合，圧力ポ アッソン方程式を解いて得られる圧力はもはや熱力学的圧力ではないため，熱的効果を考 慮する二相流に適用するには不便である．フェーズフィールドモデルは，従来の界面追跡 法では計算不可能であった核生成などの現象を予測できる可能性を有するものの，そのモ デルの妥当性を注意深く検討する必要がある．

1.3.4

フロント・トラッキング法

体積追跡法やレベルセット法はオイラー的に，BFC法やALE法はラグランジュ的に界 面を追跡する．一方，フロント・トラッキング法[34]は，流れ場をオイラー的に解き，界 面追跡をラグランジュ的に行う．界面は，２次元の場合線要素，３次元の場合三角形面要 素で構成される非構造格子により表現される(図1.3.7)．フロント・トラッキング法では， この要素をラグランジュ的に追跡するため界面位置が明確である．要素端点の位置 xp を 次式を解いて更新する． d xp d t = up (1.3.19) ここで，下付添字 pは要素端点の量を意味する．フロント・トラッキング法では，気泡全 体の正味の表面張力がゼロという事実を満足するように表面張力が評価される[72]．２次 元の場合，線要素eに働く表面張力 Fσは，弧長sを曲線のパラメーター，t を曲線の接

e

Legendre

polynomial for e

e

n

t

t

t

t x n

t(s+

∆

s)

s+

∆

s

s

t(s)

(a) (b) 図1.3.8 フロント・トラッキング法における表面張力評価 線として， Fσ= Z s σκnds = σ [t (s + ∆s) − t (s)] (1.3.20) と与えられる(図1.3.8)．３次元の場合， Fσ= Z S σκndS = σ Z s t × nds (1.3.21) ２次元の場合，e自身および隣接する線要素の端点の座標からルジャンドル多項式を生成 して t を求める．３次元の場合には，面要素端点の座標値から tを，面要素端点を通る２ 次放物面から nを求める．各要素端点あるいは隣接する面要素と共有する辺上の表面張 力は，各要素に対して求めた表面張力を平均することで求める．こうすると，気泡表面全 体の正味の表面張力はゼロとなる．表面張力は重み付き平均によって滑らかに近傍格子点 に配分される．上述のように界面追跡，表面張力評価において望ましい特徴を有する一方 で，以下の短所を有する．(1)計算の進行とともに，要素の疎密を防ぐための要素の再構 成が必要になる．その際，各要素間の隣接関係も更新する必要がある．膜沸騰や液滴衝突 などに見られるような界面トポロジーの変化に対しては，これは特に重要な課題である． このような隣接関係を考慮せずに解く工夫も最近なされている[73]が，体積追跡法のよ うにセル内の界面形状を推定する必要が生じる．(2)界面近傍の流れ場計算用固定格子上 に，各格子点がどちらの相に属するかを定義する相定義関数χ(x, t)を準備するが，計算の 安定化のためにχを平滑化する．この平滑関数を得るために，次のポアッソン方程式を解 く必要がある． ∇ · ∇χ = ∇ · Z S nδ (x − xint) dS (1.3.22)

表1.3.1 界面追跡法の特性比較[44] A: Excellent, B: Good, C: Poor

Requirements Volume tracking Front tracking Level Set BFC

(a) Volume conservation C A C A

(b) Interface Sharpness B B B A

(c) Surface tension force C A A A

(d) Density ratio A C A C

(e) Multiple Interface A A C C

(3)二相の密度比が大きいと計算が不安定になる．このため，初期のフロント・トラッキ ング法の適用例は密度比10程度の二相流に限れられていた．ただし，この問題は最近解 決されており，高密度比の計算例も見られるようになっている[74, 75]．(4)原理的に流体 体積の保存を保証していない．体積誤差の原因として，要素を移動する速度は周囲固定格 子に定義された速度からの補間であるが，この速度は発散なしの条件を満たさないこと， および要素の疎密を調整する操作が体積を変えてしまうことが挙げられる．Tryggvason ら[72]は，空間分解能が十分あれば気泡を直径の100倍移動した後の体積誤差は1, 2%程 度であるが，分解能が十分でない場合許容できないほどの誤差が生じると述べている．彼 らは，体積保存を保証するために数ステップ毎に要素端点位置を修正しているが，その場 合界面形状が変化してしまう． 表1.3.1にTomiyamaらがまとめた界面追跡法の特性比較[44]を示す．体積追跡法は， 複数の界面が混在する二相流に対する適用性が他の界面追跡法に比べて高い．また，水‐ 空気系などの高密度比二相流に適用できる．一方，流体体積保存，シャープな界面の維 持，高精度表面張力については，他の方法に劣ると考えられている．これらの欠点を改善 すれば，体積追跡法は非常に有用な方法になると考えられる． また，工業的に重要な二相流は円管内や複雑流路内に見られることが多いが，そのよう な二相流に適用するためにはデカルト座標系用に開発された界面追跡法を拡張する必要が ある．近年，非構造格子に拡張した体積追跡法が提案されている[76–82]．非構造格子に 基づく方法では，計算セル内で再構成した界面形状をラグランジュ的に移流するが，その 場合一般に体積が保存されないことが課題である．また，ラグランジュ的体積率移流法は アルゴリズムが複雑であり，既存手法の多くは２次元計算用である[76, 77, 79–82]．３次 元計算も行われているが，その移流アルゴリズムは明確でない[78]．一方姫野らは，最近， レベルセット法を一般座標系用に拡張し[83]，さらに体積保存性を改善するために体積率 を導入している[84]．彼らの方法は構造格子を用いているため，アルゴリズムは非構造格 子の場合ほど複雑化しない．しかしながら，彼らの体積率移流方法および界面近傍のレベ ルセット関数の推定方法では，界面位置が正しく考慮されていないことが課題である．

1.4

代表的な体積追跡法の性能評価

体積追跡法は，複数の界面が混在する二相流に対する適用性および高密度比二相流への 適用性が他の界面追跡法に比べて高い．また，界面の記述に体積率を用いるため，多流体 モデルや気泡追跡法と結合する統合計算手法の界面追跡機能として組み込みやすい．この ため，体積追跡法は各種工業機器内二相流に対して幅広い応用が期待できる．しかしなが ら，良好な二相流予測を行うために必須である，流体体積保存，シャープな界面の維持お よび良好な表面張力評価に課題がある．本節では，これまでに提案された代表的な体積追 跡法の概要を述べるとともに，簡単な移流計算の結果から，それら既存手法の特性および 課題を明らかにする．本研究では，それら既存手法の課題を克服した界面追跡法の開発を 目指す．なお，本節では体積率移流方程式の解法に議論を限定する．表面張力評価につい ては1.5節で述べる．

1.4.1

DA

法

体積率移流方程式を保存形に変形した上で離散化すると次式となる． αn+1 i, j,k = α n i, j,k+ 1 Ωi, j,k h

(αu1)i+1/2, j,k− (αu1)i−1/2, j,k i ∆t∆x2∆x3 + 1 Ωi, j,k h (αu2)i, j+1/2,k− (αu2)i, j−1/2,k i ∆t∆x1∆x3 + 1 Ωi, j,k h (αu3)i, j,k+1/2− (αu3)i, j,k−1/2 i ∆t∆x1∆x2 (1.4.1) ここで，Ωはセル体積，u1, u2, u3は各々 x1, x2, x3方向の速度成分である．下付添字i, j, k は x1, x2, x3 方向のセル番号で，変数配置はスタガード配置とする．∆x1, ∆x2, ∆x3 はセル 幅である． 簡単のため，セル表面(i + 1/2, j, k)を通じたx1 方向の移流のみ考える．DA法では，被 輸送体積∆Ωi+1/2, j,k を次式で与える．

∆Ωi+1/2, j,k = (αu1)i+1/2, j,k∆t∆x2∆x3

= minhαIAD|u1∆t| ∆x2∆x3+ ΩS U P, αIDΩi, j,k i _(1.4.2) ここで， ΩS U P = max h (αU P− αIAD) |u1∆t| ∆x2∆x3− (αU P− αID) Ωi, j,k, 0 i (1.4.3) αU P= max (αID, αIDU) (1.4.4)

u

∆

t

ID

IDU

IAD

IA

u

fluid1

fluid2

x

₁

x

₂ 図1.4.1 DA法におけるセルの記法 下付添字IDは流体の輸送元のセル，IDU は輸送元のさらに上流側のセル，IADは輸送が 行われるセル表面を意味する．これらの記法を図1.4.1に示した．IAは流体を受け取るセ ルを意味する．αIDA は，界面の傾きと隣接セルの体積率の値によって決める． IAD = ( ID if |n1| < |n2| or |n1| < |n3| IA otherwise (1.4.5) ここで，n1, n2, n3 は界面法線nのx1, x2, x3 方向の成分である．上式による判定の後，次 式で再度判定を行う．

IAD = IA (if αIA = 0 or αIDU = 0) (1.4.6)

また，セル IDが界面セルでない場合は IDA = IAとする．これらのセル設定を式(1.4.3) に適用して被輸送体積を求める． 式(1.4.2)は，IDA = IDの場合風上差分，IDA = IAの場合風下差分に相当する．１次 風上差分は安定性が高く，また解の単調性を維持できる．しかし，数値拡散が生じる．一 方，風下差分は不安定であるが，数値拡散の係数が負であるため界面をシャープに維持し やすい．DA法は，状況に応じてこれら風上差分と風下差分を切り替えて安定な計算を行 う差分法であるが，スキーム切り替えの判定を界面形状に基づいて行うのが特徴である． 各軸方向に分割した体積率移流方程式を順次解くという方法でDA法を多次元化した方 法は，オリジナルのDA法と区別して改良 DA法[85]と呼ばれる．しかし，現在では方 向分割による多次元化は他の方法でもよく行われている．以下ではオリジナルのDA法を DA-a，改良DA法をDA-bと呼ぶ．

1.4.2

FCT

Rudman [86]は，FCT(flux-corrected transport) [87, 88]を体積率移流方程式に応用した． FCTの概略を以下に述べる．簡単のためx1 方向の移流のみ考える．まず，数値拡散を生

じるが，安定かつ解の単調性を維持する低次の解法により移流項を評価し，αの推定値を 得る． α∗_i = αn i + F_i+1/2L − F_i−1/2L (∆x1)i (1.4.7) ここで，FL は低次の解法で評価した流束である．単調性を維持する解法により FL を計 算しているため，α∗_{は0 < α}∗ _{< 1の範囲に収まる．逆拡散を生じる不安定な高次の解法} により評価した流束をFH _{とする．これらの流束から非拡散流束F}A _{を定義する．} FA= FH − FL (1.4.8) 次式によりαを更新する． αn+1 i = α ∗ i +

Ci+1/2F_i+1/2A − Ci−1/2F_i−1/2A

(∆x1)i (1.4.9) ここで，C (0 < C < 1)は制限関数である．解の単調性を維持するように制限関数を決定 するため，α∗_{の最大値α}max_{および最小値α}min_{を調べる．} αmax i = max  αn i−1, α n i, α n i+1, α ∗ i−1, α ∗ i, α ∗ i+1  (1.4.10) ここで，以下の量を定義する． αmin i = min  αn i−1, α n i, α n i+1, α ∗ i−1, α ∗ i, α ∗ i+1  (1.4.11) Q_i+1/2I = αmax_i − α∗_i
(∆x1)i ∆t (1.4.12) Q_i+1/2D =  α∗_i − αmin i  (∆x1)i ∆t (1.4.13)

P_i+1/2I = max0, F_i−1/2A − min0, FA_i+1/2 (1.4.14)

P_i+1/2D = max  0, F_i+1/2A  − min0, FA_i−1/2  (1.4.15) RI_i+1/2= min  0, QI_i+1/2/PI_i+1/2  (1.4.16)

RD_i+1/2= min0, QD_i+1/2/PD_i+1/2 (1.4.17)

制限関数Cを次式で与える． Ci+1/2= min  RI_i+1/2, RD_i+1/2  (1.4.18) Rudman は，体積率移流方程式を各軸方向に分割し，上述の１次元 FCT を各方向 に適用するという方法で多次元化している．また，FCT ではなく TVD(total variation diminishing)を用いる方法もある．

i

j

i+1

x

₁

x

₂

case1

case4

case7

case2

case5

case8

case3

case6

case9

図1.4.2 界面形状の分類

1.4.3

PLIC

型

1.4.3.1 FLAIR

Ashgrizらが提案したFLAIR(flux line-segment model for advection and interface

recon-struction) [50] は，界面セルとその隣接セル間における界面の連続性を考慮して界面形状 を体積率から再構成する．２次元の場合，界面セルおよび隣接セルの状態は図1.4.2に示 す状態のいずれかになる．どの状態にあるのかは，セルiとセルi + 1の体積率αi, αi+1か ら判別できる．線分を表わす式は， y = ax1+ b (1.4.19) 上式をセルiについて積分するとαi，セルi + 1について積分するとαi+1に一致するとい う条件から，勾配aと切片bが定まる．セル表面を通じて輸送される領域u∆t∆x2 に含ま れる体積は解析的に計算できる．ケース9の界面形状は，さらに図1.4.3に示す４つのサ ブケースに分類される．これらサブケースのいずれに該当するかは，図1.4.4に示すサブ ケース判別ダイアグラムを用いて，体積率のみから判別できる．また，ケース1, 2, 7, 8 もさらに４つのサブケースに分類されるが，それらも別の判別ダイアグラムを用意し，αi, αi+1から判別する． FLAIRのプログラムでは，αi, αi+1の値に応じて適当なケースに分類する判別処理がほ とんどの部分を占める．このアルゴリズムをそのまま３次元に拡張することは困難であ る[82]．人見ら[89]は，体積率輸送方程式を各軸方向に分割し，それぞれの輸送に２次 元FLIARアルゴリズムを適用している．

x

₁

x

₂

a

b

c

d

図1.4.3 サブケース 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 αi+1 αi αi+1+[αi+1(1-αi)] 1/2 =0.5 3αi-αi+1=2 αi-3αi+1=0 αi=αi+1 αi-[αi+1(1-αi)] 1/2 =0.5 a b c d 図1.4.4 サブケース判別ダイアグラム 1.4.3.2 MARS 界面を傾いた線分で表現する(図1.4.5)．すなわち， y = ax1+ 1 2 (1.4.20) ただし，セルを∆x1 = ∆x2 = 1に規格化，切片を1/2に設定している．MARS [52]では， FLAIRとは異なり，隣接セル間の線分の連続性を考慮せず界面セル毎に線分を定義する (図1.4.5)．線分の傾きaを界面法線nから求める．切片を1/2としているため，線分の

O

X

_a

y(x

₁

)=ax

₁

+b

x

₁

x

₂

y(X

_a

)

∆

x

₁

∆

y

a=

∆

y/

∆

x

₁

n

∆

x

₁

= 1

∆

x

₂

= 1

図1.4.5 傾斜線分による界面の近似 下部の面積を求めると， Z _1/2 −1/2 ydx1 = 1 2 (1.4.21) となる．一般にセルの体積率は1/2ではないから，正しい界面位置，すなわち， A (X0) = Z _X 0+1/2 X0−1/2 ydx1 = αi, j (1.4.22) となるX0 を求めなければならない．0 < a < 1の場合，任意の Xに対するA(X)の値は次 式で与えられる． A (X) =                                    0 if − Xb ≤ X a 2[X + Xb] 2 _{if − X} b < X ≤ −Xa aX + 1 2 if − Xa < X < Xa 1 − a 2[Xb− X] 2 _{if X} a < X ≤ Xb 1 if Xb ≤ X (1.4.23) 上式を任意の勾配値に一般化すると次式となる．

Rsign = sign (a) (1.4.24)

Xa =

1 2 

Xb = 1 2  Rsigna−1+ 1  (1.4.26) X0 = Rsign 2 h min  1, Rsigna i₋₁h 2αi, j− 1 i (1.4.27) X0 =                                Xb−       − 2 a   αi, j− 1 2  1 + Rsign _         1/2 if Xa ≤ X0 Rsign h min  1, Rsigna i−1 αi, j− 1 2    if − Xa < X0 < Xa −Xb+        2 a   αi, j− 1 2  1 − Rsign _         1/2 if X0 ≤ −Xa (1.4.28) 上式でX0 を求めた後，隣接セルへの被輸送体積を次式で算出する． ∆Ωi+1/2, j = Z _X 0+1/2 X0+1/2−u1∆t ydx1 (1.4.29)

1.4.4

CIP

法

移流方程式の高精度な解法として知られるCIP(cubic interpolated propagation)法[90– 92]を用いて体積率移流方程式を解く．CIP法ではαの移流方程式， ∂α ∂t + u · ∇α = 0 を，αの補間多項式に基づいて解く．その際，αの勾配ベクトル g = ∇α = (g1, g2, g3)も 移流方程式を満足させるように解くことで，精度よい解を得る．移流方程式の解は次式で 与えられる． α (x, t + dt) = α (x − udt, t) (1.4.30) したがって，時刻tにおけるαの空間分布の位置 x − u∆tにおける値を，t + ∆t における 位置 xの値とすればよい．αの空間分布は次の３次多項式で与える． α = C000+ C100δ1+ C010δ2+ C001δ3 + C110δ1δ2+ C011δ2δ3+ C101δ1δ3 + C021δ22δ3+ C102δ1δ23+ C201δ21δ3 + C120δ1δ2₂+ C210δ2₁δ2+ C012δ2δ2₃ + C200δ21+ C020δ22+ C002δ23 + C300δ3₁+ C030δ3₂+ C003δ3₃ + C111δ1δ2δ3 (1.4.31)

ここで，

δ1 = −u1∆t, δ2 = −u2∆t, δ3 = −u3∆t (1.4.32)

式(1.4.31)の係数は2.3節で与える． 界面をシャープに維持するための工夫として，デジタイザー[93]が提案されている． この方法では，体積率に変数変換， H = tan " CHπ α − 1 2 !# (1.4.33) を施した変数 Hを移流する．体積率に戻す場合は次の逆変換を施す． α = 1 CHπ arctan H + 1 2 (1.4.34) CH は定数で，１よりも若干小さい値に設定する． 式(1.4.31)のCIP補間によると，1 < α, α < 0となる場合がある．その場合はCIP補 間を線形補間に変更すれば解の単調性を維持できる．しかしながら，体積率移流方程式を 解く際にこの補間変更を適用すると，体積誤差が大きくなることが指摘されている [94]． そこで，以下の計算では補間変更は行わずに，移流後に0 < α < 1となるよう強制的に値 を修正した．

1.4.5

既存手法の性能評価

図1.4.6(a)に示すように，矩形領域をその対角線に沿って直線移動する円形状を上述の 既存手法により計算した．流れ場は一様かつ定常とした．よって，直線移動後も体積率移 流方程式の真の解は円形状を保つ．本計算に用いた体積追跡法は，DA-a, DA-b, FLAIR,

MARS, FCT, CIP 法, デジタイザーを適用した CIP 法 (以下，DCIP 法) である．なお，

DA-b, FLAIR, MARS, FCTは体積率移流方程式を方向分割して解いている．一定の移動

速度は20.5_{uで，円形状の直径dに対して4.24倍の距離だけ移動させた．計算セル幅∆x} および時間刻み幅∆tは一定で，クーラン数u∆t/∆xが0.1 となるように設定した．また， dに対して割り当てたセル数d/∆xは20である．4.24d移動後(600ステップ後)の形状を 図1.4.7に示す．図中 EV，EAは体積誤差および形状誤差で，各々以下の諸式で定義した． EV = |θ (t = 0) − θ (tend)| θ (t = 0) (1.4.35) EA = X ∀i, j  αi, j(t = 0) − αi, j(tend) h 1 − αi, j(t = 0) i (1.4.36)

d

d

4.24d

5.2d

u

2

0.5

_u

(a) (b) 図1.4.6 計算領域 ここで，θは円形状の体積，tendは計算終了時刻である．また，d/∆x = 10, u∆t/∆x = 0.25 に設定して同様の計算を行った．計算結果を図1.4.8に示す．どちらの条件でも，各手法 の予測結果の特徴は同じである．まず，DA-aは体積誤差，数値拡散ともに非常に大きく， もはや円形状を維持できていない．体積率移流方程式を方向分割して解くDA-b, FLAIR, MARS, FCTは体積を非常に良好に維持した．PLIC型の解法は形状誤差が小さく，特に MARSは形状誤差がすべての手法で最も小さかった．一方，DA-b, FCTによる予測形状 は多角形に近い形状となった．CIP法は数値拡散が大きいのに対して，DCIP法は界面を シャープに維持できているが，多角形に近い形状になった．

DA-a, DA-b, FCT, CIP法, DCIP法の長所のひとつは，３次元化が容易なことである．

これらの方法を３次元に拡張し，球の直線移動計算を行った(図1.4.6(b))．球の直径dに 対して割り当てたセル数d/∆xは20で，u∆t/∆x = 0.1, 0.25の２条件とした．移動距離 は5.2d である．計算結果を図1.4.9, 1.4.10に示す．DA-aは特にクーラン数が高い条件で 球形からは程遠い形状となった．DA-b, FCTは多面体に近い形状となった．一方CIP法 の予測形状は球形に近い滑らかな界面を維持した．ただし，形状誤差が大きいことに現れ ているように，体積率の数値拡散が大きい．デジタイザーを適用すると数値拡散は改善さ れたものの，界面が滑らかではなくなった．また，方向分割を用いたDA-b, FCTは良好 に体積を保存した．これらの特徴は２次元の場合と同様である． 以上の結果より，高精度二相流予測には，体積を保存し，かつシャープな界面を維持で きるPLIC型解法が望ましいと言える．しかしながら，PLIC型解法は，例えばFLAIRの ように，計算アルゴリズムが２次元計算に限定されており，３次元化が容易でない場合が 多い．よって，３次元化が容易なPLIC型解法を開発する必要がある．

DA-a

CIP CIP with digitizer

DA-b FLAIR MARS

FCT E_V= 10.0 % E_A= 11.3 % E_V= 0.0 % E_A= 6.89 % E_V= 0.0 % E_A= 1.22 % E_V= 0.0 % E_A= 1.41 % E_V= 0.641 % E_A= 9.19 % E_V= 4.06 % E_A= 3.62 % E_V= 0.0 % E_A= 2.54 % 図1.4.7 円形状輸送結果(d/∆x = 20, u∆t/∆x = 0.10) DA-a

CIP CIP with digitizer

DA-b FLAIR MARS

FCT E_V= 20.4 % E_A= 45.0 % E_V= 0.0 % E_A= 8.62 % E_V= 0.0 % E_A= 4.46 % E_V= 0.0 % E_A= 1.10 % E_V= 1.58 % E_A= 13.2 % E_V= 13.1 % E_A= 8.34 % E_V= 0.0 % E_A= 2.27 % 図1.4.8 円形状輸送結果(d/∆x = 10, u∆t/∆x = 0.25)