B μ 左 σ 砧 : 体膨張係数 ; 粘性係数 : 熱伝導係数 : 界面張力 : 接触角 ξ, η,ζ : 計算空間の座標 : ヤコビアン : ナブラ演算子 2. 数値解析法 様々な加速度環境における自由表面流の数値解析 779 [ { 3 ( 皿 K ] [Pa si [Wf ( m K ] [

Transcription

1 778 日本機械学会論文集 (B 編 76 巻 765 号 (2010 5> 論文 No, * 様々な加速度環境における自由表面流の数値解析 (CIP 法,MARS 法,Level Set 法を協調した解法の改良 * 姫野武洋 1 *, 根岸秀世 2, 野中聡 * 井上智博 1 * 1 * 1, 渡辺紀徳, 鵜沢聖治 Numerical A 皿 alysis of Free Surface Flows under Various Conditionsin Acceleration (lmprevement of CIP LSM : CIP Based Level Set & MARS Takehiro HIMENO * 4,Hideyo NEGISHI,SatoshiNONAI (A, Chihiro INOUE,Toshinori WATANABE and SeijiUZAWA * 3 4Department of Aeronautics and Astronautics,The University f Tokyo, 7 3 1Hongo,Bunky ku,tokyo,1 工 Japah With the pr gress of human activities in space, 出 e occasion to handle liquidsin n uniform accelerati }n >r low gravity isnow growing.ollthe launchvehicles with liquid propulsion system, the dynamic acceleration during itspewered ascent or ballisticflightmakes itvery dimcult t control the position of propellantsin. the tanks For the establishment f the technology for the management of liquid propel!ant illspace vehicles,a numerical method, ca }1ed IP LSM (CIP basedlevel Set& MARS, was developed 亡 o simulate three dirnerlsionalfree surface Hows under varieus gravity conditions,which has been appl 弖 ed to clarify the dyna 而 c behavior of 1{quid propellant in the tanks of launch vchicles. Key 砺 o 鷓 : 1. 緒君 近年, 宇宙輸送系の推進機関や軌道上構造物の熱管 嚠幾器など, 地上とは異なる加速度環境で液体を利用 する場面が増えつつある. このような環境では, 貯蔵 容器内の液体を望ましい位置に保持し, 思い通りに外 部へ搬送するなど, 自由表面流を管理することが非常 に難しい. 特に, 地上から地球周回軌道まで機体姿勢 やエンジン推力を変えながら飛行する液体ロケットや, 比重差を利用した流体の駆動と安定化を期待できない 微小重力環境で長期間に亘り作動する宇宙機について, 開発コストと運用リスクを低減させるためには, その 設計段階からタンク内部における推進薬の挙動を適切に予測することが重要である. これら推進薬管 ge (t cn に関する技術課題を解決するためには, 地上では再現 が困難な 漸 に関する知見の獲得と蓄積が不可 欠であり, 理論と実験を補完する手段として数値解析手法の確立が期待されている. このような動機付けから, 筆者らは様々な加速度環 境において自由表面流の予測手法を確立すべ * * L 原稿受付 2eo9 年 12 月 9 口. Space Propulsion,Free Surface Flows,CIP LSM, Level Set Method, MARS く,CP 113 S556 東京都文 正員, 東京大学大学 ee. コ 1 学系研究科 ( 愚京区本郷 * 2 ( 独 宇宙航空研究開発機構情報計ー算工学センタ つくば市千現 211. ms ( 独 宇宙航空研究闥発機構, 宇宙科学研究所 相模原市由野台 3 1 1, E muil : uero.t.u tokyo.acjp 法,MARS 法,Level Set 法を協調的に組み合わせ, 体積保存と形状捕捉の両方に優れる数値解法として,. CIP LSM ( CIP base[1 LeyelSet& MARS (η を提案し改 良を続けている. 並行して, 落下塔や加振機を用いて 実現される微小重力環境や動的加速度環境における気 躙液界面観察実騨を行い, 対応する数値解との比較により, 手法の妥当性検証に取り組んできた. 本報では 自由表面流蠏去を構成する流体解法と界面追 跡法のそれぞれについて, 改良されたアルゴリズムを提案 するとともに様々な加速農環境における解槲列として, ダム崩壊捍題砕波を伴うスロッシング現象界面張力駆動流をそれぞれ数直的に模擬した結果を議論する. ρ ρ 7 喃 9 需 9 : 密度 : 圧力 : 温度 ; 内部工ネルギ : 速度 : 重力加速度 : 熱流束 記 号 : 音速 (Op! P s : 定圧比熱 kghn3 Pa ] [K ] [JAcg] [ mfs misz W /m2 [ mvs ] Jlakg K 68

2 B μ 左 σ 砧 : 体膨張係数 ; 粘性係数 : 熱伝導係数 : 界面張力 : 接触角 ξ, η,ζ : 計算空間の座標 : ヤコビアン : ナブラ演算子 2. 数値解析法 様々な加速度環境における自由表面流の数値解析 779 [ { 3 ( 皿 K ] [Pa si [Wf ( m K ] [N!tn ] [ ] [ m ] [ ] [11m ] 宇宙機の推進薬タンク内部では, 常温の気体と極低 温の液体がしばしば共存し, 熱交換に伴う流体粒子の 体積変化を無視できないため, 非圧縮流の仮定が不適 切な場合も少なくない, 灘動の速さは各相の音速と比 べてかなり遅く, 圧縮性流体解法を適用するのも難し い. 伝熱現象や相変化までも考慮するならば, 温度を 従属変数とせず直接的に解くのが望ましく, 熱流東を より適切に評価するためには, 界面近傍における物性値の平滑化を避け, 入工的な界面厚を極力排除したい. また, 高精度な界面追跡のためには, 数値的散逸の小 さな計算スキームで対流項を評価すべ きである. この ような方針に従って再設計された LSM の考え方を, 以下で詳しく説明する. 2.1 支 M 方程式般に, 固定格子を用いた自 肉表面流の解析手法は, 均質二相流に関する質量, 運 動量, ならびに, 内部エネルギの式, 2tL= ρ.ti (1 1t 加 P16i = Vp + (Tv + T + 廖 2 ρ 聖罵 ρ, 五 + e (3 Dt Tp = (2 μ /3( の 1 + μ( 71 露 + 爵,( 4 T σ = o δs (1 晦恥, (5 6 卩 = ( Tv + T. 譚 7 + 五, (6 を解く流体解法の部分と, 識別関数 Hs の移流方程式, H, whe [e H =: [g 0.5 : その点が液相に属する場合 魚 = 0 : その点が界面に属する場合 (8 Hs = O.5 : その点が気相に属する場合 CIP LSM の場合, 流体解法として,CIPCUP zaclo の種であり, 温度を独立変数とする皿 σthermo CP CUP 法 Cll を, 界面追跡法として PLI(>VOF ( 法 12 } の種である MARS 法を採用している 22 流体撰法以下に, 改良された TCUP 法の導 出と, 解法手順を説明する, 局所熱平衡を仮定して相 律をはじめとする熱力学関係式を用いると, 斑 1 (3 1 ま, 淋量く } =t ( u, v, w,t,p ( 1e, を状態量とする方程式系 器 里 聖み T L 聖 + 亙 (Tv ξ D 〆 Dt ρ, ら 璽嵩, 露 + β 亙 ρc 量 D ρc に変換できる, ここに現れた物性値は, β 糊 1 = 圭勉 LB 鰐 Pap. pc, 29, 3e, 驪縣数 : 音速 (Cs であり, いずれも気液各相に対する状態方程式, P = ( Hs P g ( T P + (05 s 臨 ( ア P 49, 2 (9 (10 (11 ( 12 (13 (14 を決めれば直ちに (T,p q 関数として導出される. n.1 離散化と変数配置 まず, 時間方向の離散 化にっ n いて, 第 n 時劉から有限時間掛の変化を考える. 式 (9 GD 左辺の流体微分項を時間項と移流項に分解 せず,CIP 法に基づいてセミラグランジュ し, 右辺の非移流項を (nta 時刻で評価すれば, 莞譖 髣学 諾語 κ1 的に評価 帆圭ず ( 15 (16 馴砦馴副 岬茜 暖ゾ (m を得る. なお,CimLt Nioolscm 法を採用するならば 1 ln,euler 陰解法の場合 A 1 と約束する. 次に空間方向につ いて, 図 1 に示すような検査体積 Ω を定義し, その重心を三次元般曲線座標系で, η 1 を解いて移動境界の位置と形状を特定する界面追跡法の部分とで構成される. 上記の Tv は粘性応力テンソ ル,T σ は界面張力テンソル であり, 馬は界面上で液 相側へ向けた単位法線ベクトルである. なお, 本報では, 気液問の相変化を考慮していない, 2 ( i, Jl + η! 〆へ \ + ξ Fig.1 ( 加舷 ol vo e Ω md 蜘 cs 69

3 780 様々な加速度環境における自由表面流の数値解析 (ξ,η,ζ} = (t,i k に対応する格子点と約束した構造格 子系を採用する. そのうえで, 式 15 ( 17 を Ω につい て体積積分し, 非移流項を有限体積的に離散化する. その際, 変数配置はコロケート型とし,ζ} に加えて, CIP 法の計算に必要な φ の空問勾配値 e,,,, 襲等も格子点に配置する. 必要に応じ, Ω の表面における 速度 ( + ノ 2μ を適当な補間により算出し, 検査体積の表 面 Ω に関する反変速度を, 例えば ξ 方向について, _ ( u. J _ 峠 のように与えて補助的に用いる. λ22 鯉法手順 の数値解 2 ( 網 1 UP 法において, 第 n +1 時劾 を求める手順は,CIP 法により流体微 分項を推定する移流段階と, 式 (15 ( 17 の両辺を外部 (18 反復によって釣り合わせる非移流段階とに期 される. 以下では, 表記の煩雑化を避けるため, 非移流項に作 用する空間微分の演算子 を離散化せずに記し, 階の手順を詳しく説明する. 両段 1: 移流段階 CIP 法では, 注 R する格子点と風辷 側隣接 7 点の情報を参照し,3 次多項式で内挿補間さ れた各状態量の空間分布翻叡 ξ,η,ζ を, = 娜 ζ 勢盤 (i ue ; δt( の, ノー vgr, bl,k η 曜靜 (n δ 〆り (19 に従っ て平行移動して, 各状態量の値を更新する. 同 時に, この次導関数を利用して pξ.d,,, ζ 等の勾配 値も更新される. 既報 feで ISi, 空間分布酵沸を三次元 A 型 C1P 法 (CIP A ( i3 を用いて推定し, 移流方向を, σ 謬 = σ g (i, ノ,k のように, 当該格子点の反変速度で 与えたが, 今回は, 空間分布を C 型有理関数 CIP 法 (RC [P C ( 14XI5 で推定し, 移流方向を, 瑠〇 5 協 σ 碑 δ, ノー哩傍 ω. ん δ 凸 + 万叫 9 (20 のように 2 段階 RungeKutta 法で与えるよう変更して, 移流段階の計算における数値粘性を低減させた. 時間刻みは, 計算領域に共通で, 流速に関する CFL 条件, 姻 u 釧唱 1 嚠 } 盈畦 に従って決定し, 今回,eop O.3 とした. % (21 II: 非移流段階非移流段階は, 式 (15 (17 右辺 の粘性項と伝熱頃を陰的に評価し, 速度と温度をそれ ぞれ修正する拡散段階と, それに続いて圧力項を陰的 に評価して圧力を修正する音響段階から成る. 以下, 未知の状態量 o ( +i に対して, 拡散段階と音響段階を m 回だけ反復して得られる ( 外部反復 近似解を鰍 1 り と表し,(in Z 時刻における状態 tsq { + A の第脚近似を, 1 = 1 謝 λ u 歛玄 + (1 Z 2 n (22 と与え, 非移流項の評価に用いる. また, 初期値を 1 = 銘广 ( n } とし, 初回反復を m = 1 と定義する. Il t : 拡散段階式 15 と (1 ( Dに含まれる未知量 (23 > を, 既知の第 m 1 近魄 謡 1 と修正量 δ 露 に分離して, Q 個 } 黜 + 嬬, a4 と表し, 時間項と拡散項の評価に用いれば, 拡散段階 での速度と温度の修正量,ser. を求める離散方程式 *+ 雛滋 ; (25 { 叢 : 湖 : 燈 li 饗惘 :: 彗 矗 を得る. 両式の左辺第 2 項を評価する際 怯 :1: 速度修正量 の発散. 碗島を無視すればそれぞれ, ホ l ( tsi]v {m μ 2 (Sii ; (27 ( δ 量 {; re 2 δ 異 : (28 と近似できる. すると式 25 と (26 ぽ δ 瑞 } の各成分 が互いに連成しない,4 本のスカラー方程式と表され るので, 適当な反復法を用いてそれぞれを解けば良い ( 速度と温度の内部反復 }. なお, 拡散段階では圧力の 修正を行わず, 帥 = 0 ( 29 (m とする. そのうえで, 拡散段階の第剏近似解を, 露广倉縞 + 礁こ (30 と記憶し, 拡散項を, Tv ρ :1 雫繧 チ 胤 に従っ ll 2 : 音響段階音響段階では, 拡散段階での第 m 近似解蝋に対する修正量 } 癬嘉 ; D を導入し, (33 ( + = 鉱 + 礁鵯 } 1 と表したうえで, 時間項と圧力項の評価に罵いる. 式 (15 両辺の発散をとり, 式 (17 に代入すると, 音響段階における圧力の修正量に関するスカラー方程式として, 謁留嬰 ニー 雫 匿 : 彗 鵜 1 δ :1: 31 轟 て再評価して, 音響段階の計算へ継承する t 醴習 * タ, 蕗講島転晦瑠 鉱 (32 70

4 v { 聖 RIIS(. 咢 : 翼 様々な加速度環境における自由表面流の数値解析 781 Tv 1 ρ :1} が得られるので, 反復法を適用してこれを解く ( 圧力の 内部反復. また, 式 (34 の右辺を第 m 5RHS 反復の残差 ( とし, 後の収束判定黄 38 に用いる. 圧力の修 正量に対応して, 速度と温度の修正量にっいても, 鴫 筑 黜 畷詈 :111 鰡 : と与える. 音響段階を経て得られた 副 (35 踏需 D 第 ntl 時刻における状態量の第 m 近似解が, 匪 1L + 輪鷯壇認 と与えられる. 鷓二 }; + } 簸謡 } 麺謝 + を用いれば 111: 収束判定最後に, 式 (34 で定義される残差 6 (37 を計算領域内の全ての検査体積 Ω について体積積分し, Σ ( SRIIs,! < e Σ C ω 1 (38 Ω Ω を条件として収束判定を行う. 判定式 {38 が満足されなければ, 1 を 1 だけ増分し, 寅 22} 以降の処理を繰り返す 逆に満足されれば, 囲 瓦 Q( 諞鰍丁 (39 とし, 第 1 近似解を漸近解として採用する. また, 移 流段階で必要となる空間勾配値の変化についても, 対 応する状態量の非移流段階変化について, 圃 鰯圭匝鮎 :: (4 のように. 中心差分で評価して更新する. そのうえで, 第 η 時間刻みの計算を終了し,n を 1 だけ増分して式 (19 から始まる第 n +1 時間刻みの計算に進む. M3 反復解法 項で述べたように, 式 25, (2( りおよび (34 は, 対応する Ω の内部で空間積分され, 三次元般曲線座標系で記述される計算空間で有限体 積的に離散化される. 式 34 を例にとれば, 諸 浄 毒 伽 黒 [ 讐可 : 42 とする. 魂 42 のように, 圧力項に由来する加速度が, Ω の表面で評価されることと貫性を保つべく, 式 (35 の計算においても, 例えば ξ 方向について, 詠萼 : 話響 平 + o のように, 各面を通過する体積流量の増減をまず計算 したうえで, 線形補間を用いて格子点における反変速 度の変化を算出し, 計算空問から物理空間への座標変換を行って, 対応する速度変化珂嵩 1 を求める. こうして導かれた未知数 ( λδ に関する 5 組の連立次方程式のそれぞれに対し, 本手法では不完全 LU 分解の前処理を伴う BiCG STAB 法を適用した. 即ち, 非移流段階を通じた外部反復を 1 回試行する問に, 拡 散段階で速度と温度にっ ついて 1 回の内部反復を行う. いて 4 回, 音響段階で圧力に 作動流体が非圧縮の場合は Cs = かつ B O であるから,TCUP 法は MAC 法に致ずる. 方, 任意 の状態方程式に従う作動流体の揚合, 式 16 のように 潟度と圧力の連成問題となるが, 本手法は, 温度と圧 力につ いて SIMPLE 法に類似した外部反復を行うこと で, 式 15 (17 を同時に満足する解を探索することが できるように設計されている. また, 既報とは異な り, 拡散項を陰的に評価することで, 拡散項に崗来する時間刻みの制限を受けないよう改良されている, 2.3 界面追跡法密度 ρ など, 流体解法で必要と なる各種物性値を, 状態方程式 14 をはじめとする構 成方程式に従っ て算出するためには, 棺属性の翻擱 数魚の移流方程式 (7 を離散化して数値的に解かねぱ ならない. 本研究が対象とする, 宇宙機の推進薬タン ク内における熱流動問題では, 液体推進薬の重心位置や排出に伴う残液量のより適切な推算が重要であり, 界面追跡法には, 良好な体積保存牲に舶え, 界面位置および形状の明瞭な捕捉が要求される. 2ふ 1MARS 法による体積率移流このような観点から CIP LSM では, 或 8 のような二値関数の有限体積的な移流に優れた MARS 法を界面追跡法とし て採用した. 式 (7 の離散化にあたり, まず, 図 2 に示 すように, 各格子点に対応する Ω 毎に Hs を体積積分し, 剛幅ドー / 1 ご 凝 蘿鑼廴 轡欝 1ξ1 1 IHgF (. i/ I Fig2 Re oc 霞 E 血 1c 虚 on of 畦コげ註 α } ;TangentialPlane of zer( levelset and the 洫卿 oladon of 瑞 hl e 臘 i nbased n MARS (6. 71

5 782 様々な加速度環境における自由表面流の数値解析 ここに含まれる各根の体積率 Hv を, 争ぼ腓,:lll (M 馬 :1: 響 と定義して, 独立な状態量とする. 即ち,Hv は [ 0.5, 0,5] を定義域とするボイド率である. 続いて, 般座標 ξ の定面上で Ω によって切り取られる領域を考え, そこでの面 as 率 Hl を, 4 ζ} 45 与 撮 :L Ilil :ll: 瑞 孥 に従って定義する そのうえで,MARS 法に従い, 1ξ _ il 1,, 2 の区間における He の分布を, He = max { 05. min [0 5,H } 4 9 至 } ( 紛 と区分的次関数で内挿補問し, 方向別に再溝成する. ここに現れた界面勾配瑞 } よ界面の単位法線ベクト ル恥を次項で述べる距離関数 φ を用いて評価し, H ; 2v / 鯛匹 P c whe e 丼 ξ ( ξ ン ξ ξ i. f ξ ξ η ζ と算出できる. また,1 ξ_ i1 112 の区間で He を積分 すると Hv に致する条件より, 切片塗は, ド鱒 { 去 歯帰頁画 : ぼ団 i ( 國 / 1 幽 i = ξ :if ξ 11 ( 叫 / i ( 4 Wiere ξ. fi F (50. max 1, min 1, Hl と定まる,Hv の定義域を, 通常の [O, 1 亅ではなく, [ O.5,0.5 ] として対称牲を利用したことで, 晶を与える式の記述を 2 パターンのみに簡単化できた. 流体解法の計算で必要となるのは,(ntA 時刻の値 Hsn + n であるから, 界面追跡の時間刻みを, rktr } r AOt cn + {1 Z bl n 1 ( 51 とし, 移流段階と非移流段階の間に体積率を更新する. 具体的には, まず, 式 4 のように再構成された Hg φ< φ = φ Stencil 諭 Contro1Voiume 碗 馬 〆 凝菖 Fig3 No a1 r 舜 ε of in 鰒漁 c 紐 α 戯 d 倉 dis 軸 dbu 雛 on o 埀 φ (Of. の分布を, 図 2( a に示す如く, 区間 9. [(i± 1/2,(i± II2 (, 盟声 3, 雌 (52 で台形積分すれば,Ω の表面を通じた ffs の流束 (G ζ.,,,,e fj を次のように算出できる., ズ } 53罪, H9 僻議孚吻, 掣掛 : 1 糊 1: 1, H : 撃團 1: ls, 爿割 : 岬將 鰐,,, (55 (56 (57 黄 53,(55 および 57 右辺の掃引領域に含まれる流体 粒子の重複数え上げを避けるべく, 流束を方向別に順 次推定したうえで,Hv の時間変化を, + a II 彦 }= (58 學 個 :[ 争 r:} 岬 5 ff 個 1: 糊 :鬥 1:} と与える. これは, 倒 1: [ (ntl 時刻において Ω の内部に含ま れる流体粒子が,( rr,.a 1 時刻に存在していた領域を考 え, セミラグランジュ的な移流操作を行うことを意 味する このように, 計算領域全体を過不足なく勘定 することにより,Hu の数値的散逸を本質的に回避し ながら, 体積保存を高度に満足することができる. 232 Level Set 法による界面形状捕捉 MARS 法における界面勾配の算出や, 後述する界面張力の評 価では, 界面の単位法線ベクトル恥を高精度に捕捉 することが極めて重要となる.CIP LSM では, 形状捕 捉について Leve1 set 法の考え方を採用し, 図 3 に示す ように, 界面を基準とした距離関数 φ を生成したうえ で, 法線方向 n5 = φ tv 卜 fis と曲率 κ を, それぞ糺 f 9v f d 押 具 [ 嵩 剄, (59 (60 72

6 により獲得する. φ を効率的に生成 ( 再初期化 するた めには, 式 (58 により独立変数である He + 2} を得た後, そのゼロ等高面で表わされる界面近傍において, 従属 変数である φ に関する恥血髄 Jaoobi 型方程式, (SLS 窪 { 酬ピ Q i f s v φ (iflh.1 e.50 様々な加速度環境における自由表面流の数値解析 783 韻パ即 (φ (6 り 11 N 喇 1} and S = S (iflh i, ( 50 : 界面近傍 whe e N ξ ( Xe Y ξ 4 fi, f 9 ξ η ζ を擬似時間進行法により反復して収束させれば良い. その際, 式 62 のような埋込境界条件を課せば 飾の ゼロ等高面と φ のゼロ等高面を良好に致させること ができる. (62 具体的には, まず, 相属性を判別する因子 SLsw を, 31 翻,= sgn ( の ( n + r= sgn ( FJ ( 綱 (63 で与え,φ を CIP 法により移流して, 反復の初期値 φ を, φ 欝 λ L φ 墨 λ 1 > (i ue δthi, ノー 7 さ } δtfi }, ト嶝 と定める. なお, 初回反復を Z = 1 と定義する, 次に, 式 61 左辺の計算に CIP 法を適用して, + 滝 } ipf,9 = 嬬嬉 (i iimdtit, ノ Vanf tiτ, レ馬嬲 δτ δth ( Σ 憲 レ δ τ (65 とする.CIP 法で必要となる勾配値, 移流方向, および, 局所擬似時問刻みは, それぞれ, } (φ,φ i,φ 9 閤富 (Nl, Nn, Ng 嚮 い 1 脚 蝠凱 脚姻 ノ乙跚 i,if,si1,ln! }.w 1 6r = c φ で与えればよい. 今回,c 广 α3 とした. (6 (67 (68 他亢体積率についての, O.4ggg の条件を満障驚 1 たす検査体積 Ω には, 気液界面が含まれると判別し, ここに属する格子点では, 式 65 に依らず, il; 1 n =H ; A ( 1,, 齟 3 踟 λ ニ o 國 (69 剛瑠 なる例外処理を施し, 境界条件として埋め込む. 以上, 式 {65 (69 の処理を,1 を増分しつつ所定 の上限 lmax まで反復すれば界面を基準とした距離関数の性質 1 φ1 = 1 を備えた第 近似解が得られ φ η + 1, =φ + A } 鵠 ( 70 と採用して, 界面形状の捕捉に用いる. 2A 物性値と界面張力既報では, 計箪を安定 化する便宜として, 従来の Level 默法と同様に, 界面 近傍で ρ などの物性値を人工的に平滑化する操作を施していた (D 伽 h 鰍 eme 血 od :DIM. T 自 blel P 跚 α1 rttes Physical p 脚 es. 加 Eq. ( 72. ρ :1ensdy ρ CI, : SPea7c heatat cc zstrzrzt ρ ρら 1 Cs :SbU,wivel ity ( ρ c 畫..1 _..1. 呶照騨蠏黶璽燃 ζ_. 丑..._. μ : 恥梛妙 o 噛蹴 μ しかし, 数個分の格子点幅に跨る 1 界面厚 の物理 酎慮味は不明確である. そのため, 例えば微粒化現象 のように, 慣性力と粘性力および界面張力が相互に影 響を及ぼすような自由表面流れを, 界面追跡法により 数値的に模擬する場合 ( 10, 界面厚に含まれる領域の流れを物理的にどう解釈すべきかはっきりしない, また, 界面近傍の乱流や熱貫流など, 格子幅に比してミクロ スケールあるいはメゾスケールの熱流動現象をマクロ スケールの状態量で記述される数理モデルとして考慮 する場合, 界面厚の存在が当該モデル み入れる妨げになると考えられる. を流体解法に組 これに対し, 今回の提案では, 各格子点における物 性値を, 例えば密度について, ρ ; (O.5 + 正 ly PLtg(T,p + (05 H ダ Pdes( ヱ T p (71 のように, 体積率 Hr の加重平均で算出するよう変更. した. また, その他の物性値 X についても同様に, 丿 r = (O.5 + Hp, X 君 e ( T,p + (05 fff. X (las(t,p ( 72 の形で与え, 界面を横切る物性値の跳びを 1 格子幅以 内で捕捉しようと試みた (Sh 御舳 ce Method :SIM. なお, 流体解析の計算で必要となる各物性値について, X の形をまとめ, 表 1 に掲げる. 他方, 界面張力の局所合力は, 曲率 κ を用いて,. T.. 侃些万 = 7eVH 話 φ, cr axvhtr (73 と近似できるので, 音響段階の式 (34 において, { 劉 :: 陛 : 鴛 74 のように, 圧力項と同じく, 反変加速度として運動量 式に加えた. 既報 では, 界面張力の局所合力を界面 厚の領域と重なる Ω の中心に加えていた. 方, 界面 厚を排除した今回の場合, 界面張力を Ω の表面で加え た方が, より安定な計算の進行を確認できた. なお, 拡散段階の式 (25 右辺における界面張力項についても音響段階と同様に評価し, 貫性を保っている. 界面張力に関する境界条件である固体壁面上の濡れ 性は, 便宜的に接触角命の概念を導入して表現した. 図 3 に示すように, 般座標 ξ の等箇面として壁面形状が与えられる場合, 壁面の単位法線ベクトル g と,. 73

7 784 様々な加速度環境における自由表面流の数値解析接触線に直交する単位接線ベクトル 6tは, 9. 塁, 6. φ ( 馬 φ e (75 nl 4 1 φ (9n φ 引と表されるので, 壁面上での鵡を式 59 によらず, 海., = 彦 c s ( θ + 弓 sin ( θ に致させれば良い. なお, 本手法のようにセル中心 法に基づく格子系を採用すると, 壁面上に格子点が配 置されないので, 数値解析における粘着条件と濡れ進行の矛盾は直接的には表れない. Level set 法の体積保存性を向上させる試みとしては, 筆者らの手法を含む VOF 法との組み合わせ { 邸 % 他, 多数のマーカー粒子を用いる工夫など, 様々な手法 が提案されている. それらの中にあっ て, 今回提案す る手法は, 流体解法と界面追跡法の両方に CIP 法を協 調的に用いた点が特徴である 即ち, 流体解法では, 界面での不連続を伴う密度を状態量とせず, 数値的な 拡散と散逸誤差の小さな CIP 法を適用して流体微分項 を評価したうえで, 非移流項の離散化を有限体積的に 行うことで, 界面近傍における物性値の平滑化をでき る限り排除している また, 界面追跡法でも, 界面か らの距離関数を生成する手順に CIP 法を適用することで,MARS 法と Level Set 法のそれぞれを, 三次元 般曲線座標系に拡張する困難を克服し, 体積保存と形状捕捉の両方に優れた方法を構築している. 3. 検証計算 本章では, 今回改良された CIP LSM により, 静的 あるいは動的な加速度を受けて気液問の比重差に駆動 される流れと, ほぼ無重量の環境で界面張力と濡れ性 に駆動される流れにっいて, それぞれ数値的な模擬を 試み, 手法の妥当性を評価する. なお, 以下に示す計 算では, 非移流項の評価に関し 1 尸 12 とした. また, 収束判定式 β8 について,Sc 10 4 と設定したが, い 黙 乱蓴 F 嬉 4 The 翩 oon 五 9 田洫 on ofdam b ak 血 9 P1 blen て皿 k.l ; 0 46 m w 柚 abl ずれの計算でも, 第 3 近似鰍 F3 がこれを満足した 3.1 静的加速度環境まず, 時間的に定の重力 が加わる環境に置かれた, 有限な大きさの水槽内部に おけるダム崩壊問題を数値的に模擬した結果を述べる 計算条件比較対象とする実験結果遡と同じく 図 4 に示すように, 幅 高さ 奥行叫五 4L L (L O.146 皿 の水槽を計算領域とし, 底面の左壁から L の位置に hx2h L (mo.024m の障害物を配置した 実験では水槽の上部は開放されてい るが 計算では密 閉容器とし, 対称性を仮定して奥行方向に 2 の部分 1 を行った. 作動流体 だけを模擬した三次元計 ewcase は, 非圧縮流体である水と理想気体の状態方程式に従 う空気を用いた. 初期条件として, 左壁に接する幅 高さ L 鉱の静止した直方体形状の水柱を配置し, 温 度と圧力を K IDl3x 10fPa と与え, 対応する粘 性係数や界面張力係数などの物性値を算出した. 境界 条件として, 対称面を断熱の滑り壁, その他を断熱の 粘着壁とした. 幅 高さ 奥行の方向に Xl2 個の格子点を等間隔に配置した. なお, 重力加速度は ( 廴玉 Os 02G see ( c 0.3e see 粉曜 1 プ 1 零 昌 鯨馬 i 門 8 耄可 哺鱒甲 d 0.40 蹴 ( e 0.5 sec 矗 O.IO see 02 s 030 s q40 seo O.SO sec Case1 ;Comp 晦 tiom with RCII} C and SIM Fig.5 ViSualiZedliquid 面 m of {lam breakir ]g lrrbblern 丶舳 ab 星 ockgt : 9.81m /s2.l = O, 146 m 74

8 様々な加速度環境における自由表面流の数値解析 785 ( の O.UO sec (e 0.50 scc Case 2 :CQmpu 倣 i with CIP A +IC andd 皿版 F 嬉 6 Dam hrealcing prob 匡創 m solved with a block 9,81 m!s2 であり, 水柱の崩壊開始を基準時刻广 α00s と定義した 計算結果計算結果で得られた界面形状を Hv のゼロ等高面により可視化し, 対応する実験結果と共に, 時刻 0.10 sec から0,50 sec まで並べて図 5 に示す. 水槽底面を走る水際線が障害物と衝突して跳躍した 液相の波頭形状 (c, 波頭が右壁に鋼達した後に観察さ れる液相の架橋形状 ( の, その後, 架橋により閉じ込め られた気相部分の特徴的な輪郭形状など, 実験で観 察された界面形状の変化が, 計算でもよく再現できて いることが見て取れる. 方, 界面厚を採用し, 流体微分項の評価に補間チ ェ ック付 A 型 CIP 法 (CIP A +IC を適用した場合の計算 結果 (Case2 を同様に可視化して図 6 に示す, 図 5 の 結果と比較すると, 補間チェックに伴う過大な数値粘 性により, 波頭が右壁に到達する高さが大きく失われ ており, 実験と比較して, 架橋形状 (d や気室の輪郭形状 ( e についても再現性が悪化していることが分かる. これらの結果から, 計算格子幅より小さく分裂した 液滴や気泡を解像できないものの, 数値糊生を低減し, 界面厚を排除した本手法によっ て, 加速度と比重差に 駆動される自由表面流を, 安定かつより適切に模擬で きる可能性を確認できたと言える. 3.2 動的加速度環境次に, 時間的に変化する加 速度を受ける加振容器内で生じる液体揺動 ( スロッシ ング と砕波現象にっいて, 改良された CIP LSM と物 体適合格子を用いて数値的に模擬し, 対応する実験と比較した結果を述べる 実験概要宇宙航空研究開発機構 (JAXA では, 再使用型宇宙輸送機 (RLV の実現へ 組みのつ 向けた取り として, 図 7 のような, 液体酸素と液体水 素 (LOXILHi を推進薬とする圧力供給方式の小型エ ジンを搭載した n ケット実験機を用い, ン 比較的短秒時 の垂直離着陸飛行を反復する試験 eo} を実施している. 連の試験では, 実験機が空中でエンジンを停止して 放物飛行をする運用も検討されたが, その秒時に機体 が横風突風を受けると, 推進薬タンク内で大振幅の液 体揺動が生じ, 推進薬の供給に重大な支障を生じるこ Fig,7 Reusablerocket VehicLeTest ( RVT, CDn.duσted by JA 卿 ISAS 働 へ皿剛い1 ノ ぞ } ご丶 じ驫. 芦 ド読諺鬻づごi 厂 9 i 駿.卩.,, :: 二二 : ヅ 峯 :2 毒瀞 i 申 z 工驫繍 _ { L._ 篇 _ 甲レ o,20 陥噛 呷. ; 臨 晦藷 _ ξ o,20 照 ρ 戸 1 二尸 F {g.8 Comaputed domain and g 直 d syste 皿 fbr lhe modei ta 血 k with 50 x 90 x 36 s 耄 encis 覆 蓍 1 撃鬢 10.0 奪羣ゼ 鱒 竃 G.2 窟 s o.30 歸 o. 窃 o 癖 β 贓 輜鸞 驚ゆ O 諾 霤緜毫欝 4 O 5 軍 藁澄 } 穫 o 登 蔭馨 O_OO O,26 0,50 05 鰄 飆藤潰臨呵諭 壽 鵡謠 o 滑 蝋 鰍 扇 黨羃やづ螽曩 鵡瓣篝 蚋翫 き o 昌 7s F 碆 9DiSp eme 航 and a elem 宜 ons rmvided 可 mechanical exciterp. とが懸念された, この現象を事前に評価するべ 耋 く, 筆 者らは, 図 8 の右図に示すような, 実機 LHz タンクと ほぼ相似形状のアクリル製模型タンク ( 直径比 O.20 m :O.64 m を用いた実験を行った. 具体的には, 模型 タンクを加振装置に搭載し, 図 9 に示すように, 水平 噌 75

9 786 様 7 な加速度環境における自由表面流の数値解析 D.oo sec (b 0.24sec ( c O.48 s 已唱 E 綱 m ρ 0.n sec O.00 soe 024 s (c 鰰 8 s (d α72 sec Case 3 ;Co π 甲 u 加鹹 on Fig.10 Sl( 曲 g inthemode1 寝 mk wilhout and s1(sh devioes と垂直方向の加速度が同程度になる環境を獲得すると ともに, 内部の液体 ( 常温常圧の水 が揺動する様子を 高速度カメラで撮影した. いま, 界面張力と粘性力の 影響を無視し, 慣性力と加速度力で作られる Froudle 数 の致に基づく相似則を仮定すれば, 模型タンク内部 で観察された液体揺動を, 実機タンク内部の現象に換 7 賜算できる. 実験の詳細は, 参考文献を参照されたい. 3.Z2 計算条件模型タンクとほぼ同形状の計箟領域を図 8 の左図に示す. 初期条件で, 底から O.080m の高さに液面を与え, 温度 K および圧力 1.0 atm を持っ 理想気体の空気と非圧縮液体の水を作動流体と したうえで, 実験と同じく, 図 9 に示した加速度履歴 を入力した. なお計算では, 流れ場の鏡面対称性を仮定し, タンクの半周だけを解いている, 3.23 計算結果計算結果 (Case 3 を実験と比較 して時系列に並べ, 液面を可視化して図 10 に示す. 計算でも, 横加速度の印加開始から 024s 後の波頭形 状 (b,0.48sec 後に液体がタンク頂部へ達する様子 ( c, および, それに続く砕波現象と気相巻込み ( のが明瞭 } こ 再現されている. 前節で述べたダム崩壊問題と同様, 計算の格子幅より細かく分裂した気泡や飛沫につ 比較はできないが, いて 液体揺動の振幅や砕波のタイミン グなどに注目すると, 計算は実験結果とよく致して いる, また, 液体が複雑に分裂合体しているにも関わ らず, 計算領域内における液相体積の増減は, 初期体 積に対し 0.5% 以内に収まっており, 本手法が体積保存性に優れていることを実証できた. 以上のことから, 良 CIP LSM 界面を 1 格子輻以内で捕捉する改 と物体適合格子を用い, 動的に変化する加 速度環境で, 三次元的に大変形し大移動する界面挙動を, かなり正確に模擬できる可能性を確認できた. 33 微小婁力環境最後に, 微小重力環境に置か れた円筒容器の内部で, 界面張力と濡れ性に駆動され て変形する界面について, 数値的な模擬を行って対応する落下塔実験結果との比較を行った, 3 ふ 1 実験概要実験では, 界面張力に駆動され る自由表面流の非定常現象につ いて基礎的な知見を獲 得するとともに, 初期条件と境界条件が明確に課され た容器内部の流れ場を観察することで, 数値解法の検 証に適したデータの蓄積を目的としている. 落下塔実 験の詳細は, 筆者らの既鞭鯛 に述べられている. 落下 塔設備 ( 東京大学に設置 の高さは約 10m で, 落下筐体に搭載された供試体は, 約 1 秒間に亘り 10 3G 以下の 微小重力環竟を享受できる. 本報では, 高さ 260, 内径 D 80 mm の寸法を 持つアクリル製円筒容器を観察対象とし, 液体のエタ ノールを円筒容積の半分だけ封入した場合を参照する. 落下筐体に固定された D カメラにより観察された 界面変形を, 筐体分離から時系列に並べて図 ll 上側に 示す. 重力が消失して界面張力が卓越するため, 固体 壁面上で濡れが進行し, その接触線変位が界面上を界面張力波として伝播していく様子が見て取れる 計算条件計算 (Case4 でも, 空気とエタノ ールを気液の組み合わせとし, それぞれを理想気体と 非圧縮流体に仮定のうえ, 両相に温度 K, 圧力 LOa 跏の状態における物性値を与えた, 円筒容器の直 径 D も 80 としたが, 醗 H にっいては, 端面の 影響を無視して 120 に縮めた. 陥の子午面内に 76

10 様々な加速度環境における自出表面流の数値解析 787 面とし, 流れ場に加わる加速度は無重力とした. 気液界面におけるエタノールの蒸発は考慮されていない 計算結果計算結果を時系列に並べ, 図 11 (ao.oo 弔,03 s ( a :O.00 sec (b0 40 0! 孚 3sec 混 xp 翻 m 轡 厂 Φ :e.4va sec ] Case 4 ;Compatatien (c sec (c : 瓦 80se Figll Defbrmationofliquid surface ina Cylin 面 cal vessel E 世 ianol and poly acrylate resi1 0 ; 80 mm [ 揚蟇 唇 8 貰 巳轄 籌籥 鴛巴の 鳴 Ω ミ un 議 O0ひ O 0 ロ 鰓 0 鍵舳 ーサ 呼 σ 工 2 伽 o 3 綴が侮 伍班 Fig.12 Cap 韮 y ddven 且 ows.observed defommdion は 点 を配置した, of 虹 q 由 d sur 臨 e i 皿 ac ヅ1 垣 cal vesse1 and polyecrylate resin,d = 80 る ;E 廿正 mol 周方向には 5deg おきに 3 点の格子点 式 (76 に与える接触角品の値は固気液各相の物性や 壁面粗さ, 濡れのヒステリシスなどに依存すると考え られるが, 今回は厳密な議論を避け, 実験結果のうち, 接触線が前進する進行濡れの過程にのみ注目して, & = ea (77 のように, 前進接触角 ea を定数で与えた. 壁面境界では, 速度につ いて粘着条件を課し, 温度 と圧力は近接格子点の値を 0 次外挿した, また, アク リル平面上でエタノール液滴は拡張濡れの形態を呈す るので, OA を Odeg と与えた. 封入する液体の量は円 筒容積の半分として, 中心軸と直交する平面を初期界 下側に示す. また, 計算された接触線変位を円筒壁で 追跡し, 実験データと重ねて示せば図 12 を得る. これ らから, 固体壁面上での進行濡れと, 界面張力波の伝播について, 計算結果と観察された実験結果の間に良 好な致を認めることができる. 以上より, 計算に組 み込まれる境界条件と濡れ性モデルにつ いて, 拡張濡 れの形態を呈するような固液を組み合わせた場含, 実験に基づいて式 77 の砧を Odcg と与えるとともに, 粘着壁条件を併用することで, 接触線の変位を良好に 再現できる可能性が確認された. また, 界面厚を排除 した今回の計算でも, 界面張力波の伝播とそれに伴う 界面変形を安定かつ精度良く模擬できると評価された. 4. 結言 本論文では, 様々な加速度環境における流体管理技 術の確立へ 向け, 体積保存と形状捕提の両方に優れた 自由表面流予測手法を確立すべく,CIP 法,MARS 法, Level Set 法を協調的に組み合わせた CIP LSM の大幅 な改良を提案するとともに, その計算手順を詳細に説 明した. 本手法では, 従来の LevelSet 法で採用されて いた界面厚を排除し, 捕捉するとともに,MARS 作にもセミラグランジュ 物性値の跳びを 1 格子幅以内で 法における体積率の移流操 的な離散化を考案したこと で, 体積率の散逸を本質的に回避しながら, より明瞭に界面を追跡できるよう :t 夫されている. 提案した手法を用い, 静的あるいは動的な加速度を 受けて気液間の比重差に駆動される流れと, 微小重力 環境で界面張力と濡れ性に駆動される流れにっいて, それぞれ数値的な模凝を試みた. 動的な液面変形に注 目して計算結果と実験データを比較したところ, 何れの場合でも, 両者の問に良妊な致を確認できた. 濡れ性モデルの高度化, 乱流モデルや相変化モデル の組み込みに加え, 移動境界熱流動現象につ の妥当性を評価する更なる検証など, 真服すべ いて手法 き課題 は少なくないが, 様々な環境に置かれた流れ場の解明に, 本手法が役立てられることを期待したい. 謝 本研究の部は, 科学研究費補助金 (No の助成を受けて遂行された, また, 数値計算コ 辞 ードの 構築にあたっては, 京セラ株式会社の小野智之氏ならびに宮本康治氏から, 数多くの有益な議論を頂戴した. ここに記して筆者らの深い謝意を表する. 77

11 The JapanSociety Society of Mechanical Engineers 788 rkxiigbnmaescij }L-rkSeS4 g Mft 'mb th CD X{Heewt x tu (1O '[ Yabe and P. Y. Waqg, Uijficd rntmerical pmc for oot[4rressible and ineompmssiblefuiq X Pdys. Sbc.Jqpa,B (1 Rcynol(lsW. C. and H M Sta"edoq, Liquid ptepellant VoL60K1991Xpp.210S bchavier at lowand eero g, 11aeIlyuamicBehauibl of/[ligztitis (l1. T. }liinenoand 1'.Watanahe, 'lhermo-ftuid Managetr!ende in Mom'ng Cbntainets, H. N. Almamscm, tedxacas4 su106 under I.ow-gradtyConditionsi 1stRepo4 TC(JI'mchod for (1966>,pp theanalysisofthermo-fluid,pherk}merva, Tnzns.JSME Sin:a (2 B. N. Antarmid V. S,NuetiotAntarl Fi(ndo,rentaty ofluw Vel.69,No.678 (2003, pp (injapariesej. 6avityNtiidIlynamics and Hb`ti1}tmsf14 (1993 CRC (12 D. L. YOmlgs,Ntunend Adletim[Mis.t2m Ilhmlsllynamics, Prcss,JSBN O ]3-5. editedbyk, W, Morton and M. J.BainesL,(198ny Acadernic (3> T. \abe and E. TtikrL A r}ew higher-order Godnov method I'tesEgMassachuseUs. forgenerej1ryperbolic equanian$j Ph}ay,Sbe.Jqpwz,Vel.57 (l3t, Yabq T.Ishikaw and Y. Kadeta,J IVo4stS,c.J/ua (1988Xpp.es9&2601. Vol.5g(19eq,pp.2301H23os. (4 T.Kurmgi,MARSfiN7rruudtiphasecalculatio4Cmpctiind (14 T, Acki,Cewnpaatio,mlnuid byncllnics.j:vol 4,(1995, lhid tlyncvnks Jtntmcd.,Vbl.9, No,1(2(XKI}. (5 M. ssgman,p. Smerelca, and S.Osiher, pp A 1evelset appma[in (ISF. Xiao.T. Yabe arkl T. ito,convnt PiwF.Cinnmma, Vol. lirroorrrpiiting solutions to incompressible two-phasc fiow,j 93,(l996Xpp.1. Cbmptct.P}eLs., Vol.114(1994, pp.t (6 T.HirTierro and 1'.Waianabe, 'Ihettpo-Fluid (16 inouc,c.,watarvabe,t.,hirtierk, T.and Uzawa S.,Efirect ancl Manqgernent Mochanistn of iniector internalflow or Liquid Sheet under l.ffw-grrrvityoenditions : pt Report, FneSurihce IlyTLurnics and Atomizrtian: Effectef injecti(m Veloclty Flows DrivenbySurfaceFer[x}aTrans.ISuaber.B, VoL68, Profile,7'tunstXSma &inbl Vbl.76, No.765(2010 (in Ne.6g7(2003,pp OnJapanese Japa;nse (7 T. IIimeno,S. Nonaka, Y. Naruo, Y. Inataniand T. {17 S. P. van der Pij1, A. Scgal,and C. Vuilt A Mass Wut`fNumesical Analysisof Slosliingand Wave Conserving Level Set Method for Mode}ing of BreakipginaSmallVesselbyCM-I.SM", JSME Multi-PhaseFlows, lnt,x IVitmer,Meth,FVuicts, Vol. 47, 1ittemaioual,Jbtarnae SleriesBL Vol.47, No.4eOe4, (2e05, pp pp.709r715. (18D. Enrighg R Fedkiw, j.ferziger,and I.Mitchell, A (8 T.H{rneno,T.WatflnabeandRiniai,`ENumaricalAnalysis of Fteebim FkJws Drivenby hmial 'rensien Hybrid Panicl Leyel Set Method foc Improving EEfogt", Interface Capturing, Jeurnat ofcemputational Pnptsics, JbiLualof thejtpan 5'ecity of Micng`ii,tty opiecutinag Vol, 183,(2002, pp, Vol.23,No.2{2006>pp.99-10S. (19KosihizLika, S..Taman H., arid Oka Y, Convnttaticomt (9 T.HimengT.Watanabe,S,Nonakti,Y.Narucr,Y.Inataniand blindlbuamicsj VoL ag(1995, pp.29=i6. H Acki, Numerica!and ExperimeTtalIirvestigatiori (20 Y, lrmtat}i. Y. Naruo atidk. Yoncmetq J imrtvi cnd Sloshingin Rockst Tanks with Dafrrping Deviws AL4A Rechety A:4A,VoL3& No.1 (2ootX pps6-42 Poper2(M1555K2eOD -78- NII-Electronic Library