hp120250 「京」以外 HPCI 一般利用 HPCI other than K General Use
2 つの自由表面を有する薄い液膜挙動の解析
Numerical analysis of thin liquid film with two free surfaces
白崎 実 Minoru SHIRAZAKI
横浜国立大学 大学院環境情報研究院
Graduate School of Environment and Information Sciences, Yokohama National University 要旨 薄い液膜は我々の身のまわりの自然の中でもよく見られる他，熱交換器や液体膜電極などの産業 機器にも応用されている．しかし，薄い液膜の挙動を把握することは，実験的にも数値解析にお いても一般に難しい．とりわけ，シャボン玉に見られるような，その両側に自由表面を有する薄 い液膜は特に難しいと言える．本研究では，流入する気体によって変形するシャボン玉の定量的 な解析を目的とする．直交格子系においてLevel Set 法を用いて，主に表面張力がシャボン玉の挙 動に与える影響について調べた．その結果，表面張力の大きさがシャボン玉の運動に大きな影響 を与えることが分かった． キーワード：計算流体力学，自由表面流れ，非圧縮性流れ，混相流，薄い液膜 Abstract
Thin liquid film can be found often in nature and is applied in a lot of industrial devices such as the heat exchanger, the liquid membrane electrode and so on. It is generally difficult to acquire the knowledge of the behavior of thin liquid film by experimental and numerical analysis, especially when it has free surface on both sides like soap bubble. This study aims at qualitative analysis of soap bubble deformation induced by gas injection. The Level Set Method in the Cartesian Grid is used for investigating the influence of surface tension on the behavior of soap bubble. The numerical results show that surface tension have a major effect on the motion of soap bubble.
Keywords: Computational fluid dynamics, Free surface flow, Incompressible flow, Multiphase flow, Thin liquid film
© 2019 Research Organization for Information Science and Technology All rights reserved. Received: 16 January 2018
Accepted: 19 April 2019 Available online: 22 April 2019
1. 研究の背景と目的 自由表面を有する泡は，飲料の充填時や排水処理工程など産業界において様々な場面で見られ る．製造ラインで発生した泡は噴きこぼれによる充填効率の低下や，製品の不良化率の上昇，配 管ポンプの摩耗促進の原因となることが知られている[1]．このため，旧来より機能，用途に応じ た消泡技術が開発されてきている[2]．近年，より効率性，汎用性を高めた消泡装置，消泡剤の研 究が化学メーカーを中心に活発に行われているが，消泡過程における液膜，特に自由表面を有す る液膜の挙動に着目した研究はないようである．その理由の一つとしては，自由表面を有する液 膜の性質が挙げられる．液膜は，絶対量としては小さな変形であっても，膜の厚さと比較すれば 大変形となることから非線形性が強く現れる．更に，そこへ液体の粘性による散逸効果や界面の 乱れによる非一様な表面張力が加わるため，一般に液膜は複雑な挙動を示す[3]．液膜には，片側 だけが自由表面となっているものもあるが，両面が自由表面である液膜ではこれら複雑な挙動は 特に顕著に表れる． このような両面に2 つの自由表面を有する液膜の力学的挙動を把握することは，工学的な面か らも，また理学的側面からも意義深いと言える．2 つの自由表面を有する液膜のうち，もっとも 単純で身近なものとして図1 のようなシャボン玉が挙げられる．我々はシャボン玉を具体例とし て，流れと薄い液膜とが連成する現象についての研究を行ってきた．シャボン玉は，吹き込まれ る空気によって液膜が大きく変形をしながら成長し，それにともなって膜厚は薄くなっていく． シャボン玉の成長や変形は身近な現象でありながらも非常に複雑であり，その液膜の詳細な挙動 や液膜内部の流動を捉えることは実験的にも数値計算を行う上でも難しいと考えられる．シャボ ン玉に関する研究は，映像素材としての期待からコンピュータ・グラフィックスの分野でなされ ている[4]ものの，シャボン玉の成長過程における力学的挙動について気液二相を考慮した研究は 見られない．そこで，本研究課題では流入する気体によって成長，変形するシャボン玉について， 数値計算によりその力学的な挙動について2 次元，3 次元での解析を行い，各パラメータがシャ ボン玉の成長や挙動，そして破断に与える影響について調べるとともに，今後，さらなる大規模 詳細解析へと発展させるための基礎的な確認と検証を行うことを目的とする． 図1 実際のシャボン玉（上）と内部の流速分布の計算結果（下）
2. 計算モデル
気液二相流解析では，直交格子上で界面形状を表現することが出来るVOF 法や Level Set 法[5] が広く用いられている．シャボン玉の挙動においては表面張力による影響が大きいことが予想さ れるため，今回の計算では界面の法線ベクトルや曲率を精度よく算出することが可能なLevel Set 法を採用した．計算手法は以下の通りである．気相と液相とを統一的に扱い，支配方程式である 非圧縮性流れの連続の式とNavier-Stokes 方程式を連立させ，流体の流速と圧力を算出した後，気 相と液相の界面をとらえるためにLevel Set 関数の移流方程式を解く．表面張力の評価には Level Set 関数から求まるデルタ関数と密度および平均密度を用い，気液界面に仮想界面幅を設定して 計算を行った．一般に，時間が経過するにつれLevel Set 関数は距離関数としての性質を満たさな くなり，質量保存を維持しなくなる傾向がある．そのため，Level Set 関数の再初期化と体積補正 [6]を行う．この際，Level Set 関数の再構築に関しては Sussman の再初期化式[7]を用いた．計算 領域はStaggered 格子によって離散化し，非圧縮性流体の解法として Fractional Step 法を用いた． 移流項には5 次精度 WENO 法，粘性項には 2 次精度中心差分法，圧力 Poisson 方程式の解法には Bi-CGStab 法，時間積分には 3 次精度の TVD Runge-Kutta 法を用いた． 具体的なモデルとしては，図2 に示すように，ストローから吹き込まれた気体により薄い液膜 であるシャボン玉が膨らむ現象を想定し，2 次元および 3 次元解析を行った．2 次元モデルでは 図2 の y 方向中央位置での xz 断面を計算対象としている．境界条件は，流入口に Poiseuille 流れ を与え，その反対側壁面では流出条件を他の壁面では滑りなし条件を，ストロー表面では各方向 共に滑りなし条件を課し，初期の液膜厚さを 0.2 cm に設定した．これ以降は，流入口に与える Poiseuille 流れの最大速度を流入速度と呼ぶことにする．シャボン玉の初期形状は半球状であると し，ストロー壁面との接触面で粘着し続け膨らんでいく．計算領域として，直交等間隔格子で 2 次元計算では各方向最大600 分割，3 次元計算では各方向最大 400 分割した計算を行った．初期 状態の液膜は厚さ方向に2 次元計算では 30 格子，3 次元計算では 20 格子に分割されている．気 相には20 ℃の空気（密度 1.2 kg/m3，粘性係数1.8×10-5 Pa･s）を，液相には 20 ℃の 20 %グリセ リン水溶液（密度1045.9 kg/m3，粘性係数2.051×10-5 Pa･s) [8] に界面活性剤を混ぜ，表面張力を 落としたシャボン玉溶液を想定した．なお，液膜の界面活性剤濃度は常に一定とし，図 14 の計 算を除き重力は考慮していない． 3. 並列計算の方法と効果（性能） 計算には京都大学学術情報メディアセンターのスーパーコンピュータCRAY XE6 を利用した． 2 次元計算に対しては OpenMP による 1 ノード内（最大 32 コア）での並列化を行い，32 コアを 用いた場合におよそ 14 倍の速度向上を得た．また，3 次元計算では領域分割にもとづいた MPI による並列化を行い，64 ノードでノード当り 4 コアの合計 256 コアを用いた計算では，64 ノー ドでノード当り1 コアの合計 64 コアを用いた計算との比で 3.16 倍（理論値 4 倍）の速度向上を 得た．これらの結果を図3 および図 4 に示す．
図2 計算モデル 4. 検証計算 開発したコードの妥当性を確認するための検証を行った．シャボン玉は液膜によって形成され ており，シャボン玉の液膜挙動には表面張力による影響が大きいと考えられることから，2 次元 液膜流れの計算と，表面張力の効果が強く働く3 次元単一気泡上昇の計算について，既往の計算 結果との比較を行った． 4.1 2 次元液膜流れ 図5 に示すように，計算領域左側に高密度流体，右側に低密度流体が配置されている状態を初 期の状態とし，重力によって高密度流体が液膜となり流れていく計算を行い，Bonometti らの計 算結果[9]と比較した．格子分割数は x，y 各方向に 3328×640 とし，計算条件は Bonometti らが行 ρ/ρ = 10-3 Reynolds 数 Re = 5.0×104 図3 2 次元計算における速度向上 図4 3 次元計算における速度向上
をH，横方向長さ L = 13.5H，初期の自由表面位置 x0 = 3.25 H として計算を行った．ここで，ρ は 流体の密度であり，Reynolds 数は H L H H H H L H gAtH H Re= , At =
,
=
(1) で定義されている．g は重力加速度，μは粘性係数であり，下付き添え字H と L はそれぞれ高密 度流体，低密度流体の値であることを表わす．なお，At は Atwood 数である． 計算領域の縦方向長さ H で無次元化した高密度流体先端位置 * / H H x x Hと低密度流体先端位 置 * / L L x x Hの時間変化のグラフを図 6 に示す．横軸が無次元時間t*t/ H gAt/ ，縦軸が流体 先端位置を表しており，紫の実線で示したものが開発したコードでの計算結果，黒の実線と点線 で示したものがBonometti らの計算結果である． * H x は計算開始直後に多少の差は見られるが，そ の後は概ね一致している． * L x は良好に一致しており，これらより概ね妥当な結果が得られたと言 える． 4.2 3 次元単一気泡の上昇 3 次元単一気泡の上昇の計算を行い，Balcazar らの計算結果[10]と比較した．図 7 に示すように 計算領域は6cm×6cm×8cm であり，高密度流体中にある低密度流体（気泡）が浮力により浮上す る現象である．低密度流体は直径d = 1cm の球形であり，その中心は下壁から 1cm の高さに位置 している．格子分割数は各軸方向に150×150×200 としており，このとき気泡直径は 25 格子分の 大きさとなっている．壁に対する流速の境界条件として全壁にno-slip 条件を課し，圧力の境界条 件として全壁に法線方向圧力勾配ゼロを課した．二流体の密度比，粘性係数比はいずれも100 で あり，(2)式で定義される Eötvös 数は Eo = 97.1，Morton 数は M = 0.971 として計算を行った．こ れらの条件も全てBalcazar らの計算条件と合わせている． 4 2 2 3 H H g gd Eo= , M = (2) 図5 2 次元薄膜流れの計算条件 図6 各流体の先端位置 𝑥 𝑥 Heavy fluid Light fluid 𝐿 𝐻 𝑔 𝑥 𝑥 𝑦ここで，∆ρ は二流体の密度差H ，L σは表面張力係数である．横軸を無次元時間t*t/ d g/ とし，縦軸を(3)式で定義される Reynolds 数としたグラフを図 8 に示す．ただし，それぞれの文 献[9, 10]における定義に合わせるため，ここでの無次元時間と Reynolds 数は 4.1 節での定義とは 異なることに注意されたい． H H wd Re=
(3) このReynolds 数は気泡上昇速度 w を変数としている．また，横軸を無次元時間t ，縦軸を無次元* 化した気泡重心のz 座標z*z d/ としたグラフを図9 に示す．いずれの図においても，紫の実線 で示したものが開発したコードによる数値計算の結果，黒の線で示したものがBalcazar らの計算 結果である．黒の線によるM1～M4 は格子解像度が異なる結果であり，M1 が最も粗い格子，M4 が最も細かい格子である．本計算はM3 相当の格子解像度での計算としている．これらの結果に ついても良好な一致を示しており，開発したコードによる計算は妥当であると考えられる． 図7 3 次元単一気泡の上昇の計算条件 図8 気泡上昇速度に対応した Reynolds 数の時間変化5. 計算結果 開発した計算コードを用いて，気体が吹き込まれることにより成長するシャボン玉の挙動につ いての解析を行った．これにより得られた結果をまとめると次のようになる． シャボン玉の変形挙動は，予想されるように，吹き込まれた空気によるシャボン玉内部の流れ と液膜内部の流動，液膜の変形とその伝播，表面張力の影響が相互に連成した非常に複雑なもの となる．例えば，図 10 に示す 2 次元モデル計算結果のように，シャボン玉の形状は，ストロー から流入する空気の流れとシャボン玉に働く表面張力の効果により，液膜は空気の流入方向への 伸びと，その垂直方向への伸びを交互に繰り返して「振動」するように変形する．さらに，空気 の流入によってシャボン玉が大きくなるに伴って，液膜の厚さは薄くなるが，液膜がある限界以 上に薄くなって不連続になると液膜は破れてしぼんでしまう．ここでは，2 次元モデルの場合は 液膜が途中でちぎれてしまうこと，3 次元では液膜に穴があくことを「破断」と定義する． 図11 のようにシャボン玉のサイズ X，Y を定義し，初期形状である半球時のサイズ X0，Y0か らの変位の時間的な変化を図12 に示す．ここで，液膜の表面張力係数σを0.02N/m，0.03N/m， ストローからの空気の流入速度Uinを0.1m/s，0.25m/s と変化させている．条件によってグラフの 右端の時刻が異なっているが，このグラフの右端の時刻が各シャボン玉の破断する時刻に対応し ている．図12 から，シャボン玉成長時の振 動の「周期」は表面張力が大きいほど短く なり，流入流速が小さいと両方向に細かく 振動しながらシャボン玉は成長することが 分かる．流入流速が大きいとシャボン玉の サイズの変化も大きく，表面張力による変 形の抑制も相対的に大きくなる傾向が見ら れる．この傾向は，3 次元計算でも同様で ある．シャボン玉が破断する時刻や破断す る位置は，表面張力をはじめとする各パラ メータの値に依存する．3 次元モデルにお いて，x 方向に吹き込まれた空気によって， シャボン玉はx方向とyz方向に振動しながら 成長し，破断に至る計算結果の例を図13 に 示す．さらに，空間解像度が低いと早い段階 で破断することも確認している．シャボン玉 が成長した際の液膜の厚さを考えると今回 の格子分割数では空間解像度が十分である とは言えないため，今後，より高い解像度で の詳細な解析が必要であり，これによってシ ャボン玉はより安定する可能性がある． 図10 シャボン玉の成長 （2 次元モデル，流入流速 0.25m/s， 表面張力係数0.02N/m（左），0.03N/m（右））
図11 シャボン玉の x, y 各方向のサイズ
図12 シャボン玉成長時の変位の時間変化（2 次元モデル，x 方向（左），y 方向（右））
また，今回の計算モデルの場合，ストロー先端部（固体部分）と液膜の接触部の取り扱いは， この部分でLevel Set 関数をどのように処理するかに帰着する．現在はこの部分で Level Set 関数 が移流しないという処理をしているが，このストローと液膜の接触部分については，シャボン玉 の安定性と密接な関係があると考えられ，今後の課題でも述べるように，シャボン液の適切な供 給を含むより現実的なモデルを導入する必要がある．最後に，図 14 は重力がシャボン玉の成長 に与える影響を調べたものである．重力により液膜が下方向に流動することで，液膜厚さの不均 一が顕著になり，破断が早まる傾向が見られた． 図14 シャボン玉形状と流速分布（2 次元モデル，重力なし（左），重力あり（右）） 6. まとめと今後の課題 流入する気体によって成長，変形するシャボン玉の力学的な挙動について，CFD による 2 次元 モデル，3 次元モデルでの解析を行い，各物理パラメータが成長過程におけるシャボン玉の成長 や挙動，そして破断に与える影響について調べた．その結果，気体の流入速度や表面張力などの パラメータがシャボン玉の変形挙動に与える影響について次のような知見を得た． 流入する空気によって成長するシャボン玉は，振動しながら大きくなるが，表面張力の影響 を強く受ける．振動の「周期」は表面張力が大きいほど短くなり，また流入流速に依存して その挙動が変化する シャボン玉が「破断」する時刻や位置は，表面張力をはじめとする各パラメータに依存する 重力により液膜厚さの不均一が顕著になり，破断が早まる傾向が見られる そして，単純に思えるこの現象も非常に複雑な側面を持っていることが確認でき，さらに詳細な 解析が必要であることも分かった． 今後，対応および検討すべき課題をいくつかを挙げると次のようになる．まず，空間解像度を 高めるために，より大規模，高解像度な計算の実施が挙げられる．単に大規模にするだけではな く局所的な細分化をはじめとした効率的な詳細解析技術の導入が必要になると思われる．また， 計算モデルの見直しとして，実際のシャボン玉では，膨らむことに合わせてシャボン玉溶液はス
トロー内側から適切に供給され続けていると考えられ，今回，これが計算において考慮出来てい ないことがシャボン玉の挙動が不安定となっている原因の一つと考えられる．また，本研究課題 では表面張力係数を一定としているが，界面活性剤の濃度分布がシャボン玉の安定性に寄与して いる可能性もある．ストローと液膜の接触部分の境界条件の改良や接触角の考慮なども含め，よ り現実的なモデルに近づけていくことが必要である． 参考文献 [1] 石井淑夫, 『泡のエンジニアリング』, テクノシステムズ (2005) [2] 菅原隆, 『泡コントロールと消泡・脱泡事例集』, 技術情報協会, (2007) [3] 吉永隆夫, 数理解析研究所講究録, 1271 巻, pp. 135-144 (2002)
[4] X. Wei et al., Proc. 2003 ACM SIGGRAPH/Eurographics SCA, pp. 75-85 (2003) [5] S.J. Osher and J.A. Sethian, J. Comp. Phys., Vol. 79, pp. 12-49 (1988)
[6] K. Tsubogo, 応用力学論文集, Vol. 6, pp. 201-208 (2003) [7] M. Sussman et al., J. Comput. Phys., Vol. 114, pp. 146-159 (1994)
[8] (社)日本機械学会, 『技術資料―流体の熱物性値集』, (社)日本機械学会, pp. 477-478 (1983) [9] T. Bonometti et al., Journal of Fluid Mechanics, Vol. 616, pp. 445-475 (2008)
