第 26 回数値流体力学シンポジウム
講演番号 C02-3

気泡の合体の効果を考慮したキャビテーション流れのための気泡モデル
Bubble Model Considering Bubble Coalescence Effect for Cavitating Flow
○ 田村善昭, 東洋大, 〒350-8585 川越市鯨井 2100, E-mail: tamura@toyo.jp
鶴見伸夫, 東洋大院, 〒350-8585 川越市鯨井 2100, E-mail: tsurumi@cse.eng.toyo.ac.jp
松本洋一郎，東大，〒113-8656 東京都文京区本郷 7-3-1, E-mail: ymats@mech.t.u-tokyo.ac.jp
Yoshiaki Tamura, Toyo University, Kujirai 2100, Kawagoe, Saitama 350-8585 Japan
Nobuo Tsurumi, Toyo University, Kujirai 2100, Kawagoe, Saitama 350-8585 Japan
Yoichiro Matsumoto, University of Tokyo, Hongo 7-3-1, Bunkyo-ku, Tokyo 113-8656 Japan
In this paper, an improved cavitation model for cavitating flow simulations is proposed. The present model based on
the previously proposed bubble model in which the Rayleigh - Plesset equation (R-P) is solved to obtain the volumetricmotion of a bubble. In the previous model, cavitation bubbles never coalesce while the present model includes theeffect of coalescence into R-P. As a result, developed cavitations, especially sheet cavitations, are expected to staylonger than in the previous model. Flows around Clark-Y 11.7% and NACA0015 are simulated to show the validity ofthe present model, compared with the previous results and the corresponding experiments.

１．はじめに	 キャビテーション流れの数値解析は気液両相の空間・時間スケールが異なるため何らかのモデル化をする必要がある．これまでに提案されている解析手法には主に 3 つある．1 つは VOF 法に基づくもの(1)(2)でキャビティが計算格子よりも大きいことが前提であり，シートキャビテーションのように大きな気泡を解くのに適している．この対極にあるのが均質流モデル(3)(4)と呼ばれるもので，気泡は計算格子よりも十分小さく，格子内では均質に分布しているとする方法である．3 つ目はちょうどこれらの中間で，気泡は計算格子よりも小さいが流れに影響する程度の大きさがあり，周囲の圧力に応じて体積変化をする他，液相との間ですべり速度をもつ気泡モデル(5)である．気泡モデルは，気泡の運動を記述できるため，キャビテーション気泡の移流や集積を計算できる，気泡が収縮する際の圧力を推定できるなどの利点がある一方で，計算が不安定になりやす
い，気相の体積率（ボイド率）が上がって気泡同士の干渉や合体が起こるようになるとモデル化が難しい等の欠点もある。
著者らは，気泡モデルに基づくキャビテーション流れの解析について，基本的な解法(6)を構築し，高ボイド率でも計算できるモデル(7)(8)
を考案してきた．ここでは文献(8)に加え，ボイド率が上がると起こるであろう合体の効果を取り入れる手法を新たに考案したの€で，これについて報告する．後で述べるように，合体の効果は第一
義的には気泡の成長速度に影響を及ぼすので，気泡の半径運動を記述する Rayleigh-Plesset の式に合体の効果を含めるようにした．この手法を Clark-Y11.4%，NACA0015 の二種類の２次元翼型まわりのキャビテーション流れに適用し，実験結果および従来の方法(6)で計算したものと比較し，本提案手法の有効性について検討する．

・ 液相は気相に比べて十分に大きな熱容量を持ち，気液界面における熱伝達による温度変化はない．
これらを元に，以下の基礎方程式が導かれる．
1) 気泡にかかる力の釣り合い式(1)
FIi + FAi + FPi + FDi + FLi + FHi = 0
ここで，FIi は気泡の慣性力，FHi は履歴力だが無視する．FAi は付加慣性力で，

€€

⎧⎛
3
∂ rG3 ρ L uGi ⎞⎪ ∂ rG ρ L uGi4⎟πβ ⎨⎜+ uGj⎟⎜∂t∂x j3 €€⎪⎠
⎩⎝

()

()

(2)

⎫⎛ ∂ r 3 ρ u∂ rG3 ρ L u Li ⎞⎪G L Li⎟⎬−⎜+ u Lj⎟⎜∂t∂x j⎠⎪⎝
⎭

()

()

βは球形気泡なので 1/2 とする．添字の L は液相，G は気相，i, jは x,y 方向を表す．u は速度成分である． FPi は周囲流体の加速に€よる力で，
FPi = −⎛ ∂u∂u Li ⎞4
⎟
πrG3 ρ L ⎜⎜ Li + u Lj
3∂x j ⎟⎠⎝ ∂t(3)€FDi は抗力， FLi は揚力で，1
 
FDi = πrG2 ρ L C D uG − u L ( uGi − u Li )€2
€ F = 1 πr 2 ρ C u − u ε ω ( u − u ) / ω
LiG L L GL i, j,k Lk GiLiL2(9)ω L は渦度ベクトルである．抗力係数 と揚力係数(10)は，24
CD =1+ 0.15Re 0.687
bubRebub€ 2r ρ u − uRe bub = G L G L(€２．モデル化と計算手法	 
まず，文献(6)-(8)に共通するモデル化のための仮定は以下の通りである．
・ 液相は非圧縮性，気相は圧縮性とする．
・ 気相の密度及び運動量は液相に比べて十分に小さい．
・ 気相は全て気泡で，個々の気泡は球形を保ち，合体や分裂は起こらない．
・ 気液間の相変化による質量変化は液相の持つ質量に比べて十分小さい．
・ 気泡内は，不凝縮ガス及び蒸気で満たされているとし，内部で圧力は一定とする．気泡内の不凝縮ガスは気液界面での蒸発と凝縮を模擬するようにその圧力変化をモデル化する．
€FAi =(4))(5)µL⎛ ω r ⎞0.25
C L = 0.59⎜⎜  L G ⎟⎟⎝ uG − u L ⎠で与えた．
2) 気泡数密度の輸送方程式€1()∂nG ∂ nG uGj+=0∂t∂x jnG は気泡数密度である．
3) 気泡の体積運動2€
D 2 rG 3 ⎛ DrG ⎞ pB − pL 1   2rG+=+ uG − u L⎜⎟ρ
4
Dt 2 2 ⎝ Dt ⎠

(6)

(7)

Copyright © 2012 by JSFM€

第 26 回数値流体力学シンポジウム
講演番号 C02-3
また，(7c-2)式ではすべり速度の項も落としているが，これは，気泡径とは別に，ボイド率が 10%を超えた場合にすべり速度を 0 と置く処理をしているためである(8)．
さらに，本研究では，ボイド率に応じて気泡の合体の効果を考慮する．まず，(13)式より，
4(14)fG = πrG3 nG3であるが，気泡数密度が変わらないとすると，両辺の微分をとって，
f˙G = 4 πrG2 nG r˙G˙f˙ = 4 πr 2 n r˙˙ + 8πr n r˙ 2€

これは，いわゆるRayleigh-Plessetの式である．
Dは実質微分を表す．
pB は気泡内圧力で1 DrG
2T	  	  	  	  	  	  	  	  	  	  	  (8)pB = pv + pG −− 4µLrG DtrGと表される．T は表面張力，pV は蒸気圧で共に一定とする．pG は不凝縮ガスの圧力で，気泡は膨張時に等温変化，収縮時に断熱変化に€ 近い変化をすることから，
pG rG3 = const. (DrG / Dt > 0)(9)pG rG3κ = const. (DrG / Dt < 0)とモデル化した．ここでκ は比熱比で 1.4 としている．
4) 圧力方程式
1 ∂pLDr(10)€+ ∇ i f L u Li + ∇ i fG uGi − 4 πrG2 nG G = 0Dt
c 2 ∂t €この式は擬似圧縮性法の考え方に基づき，液相体積率の保存式や気泡数密度の式等を用いて導出している(6)．	  	 € 5) 気泡流の運動量保存式(∂( ρ L f L u Li ) ∂ ρ L f L u Li u Lj+J∂t∂x jG)(11)

1= −∇i pL + ∇ k ( µ∇ k u Li ) + ∇ i ( µ∇ k u Lk )3
ここで，気相の密度は小さいとして無視している．また， f L は液相の体積率， µ は混相流の実効粘性係数で，⎧⎛ µ + 5/ 2 µ ⎞⎫€G(12)µ = ⎨1+ fG ⎜ L⎟⎬ ≅ (1+ fG ) µ L€⎝ µ L + µG ⎠⎭⎩€で与える．
f はボイド率である．
G6) 体積率の拘束条件
4€(13)
fG + f L = 1 = πrG3 nG + f L3€
以上が解くべき式となる．
次に，高ボイド率でも計算ができるようにするため，以下の操作を行う．まず，(7)式を，移流を計算する
2 式と気泡の体積運動そ€のものを計算する 2 式の計以下の 4 式に分ける．
∂rG∂r(7a)+ uGj G = 0∂t∂x j∂r˙G∂r˙(7b)+ uGj G = 0∂t∂x j2€∂ 2 rG 3 ⎛ ∂rG ⎞ pB − pL 1   2+ ⎜ ⎟ =+ uG − u L2
ρL2 ⎝ ∂t ⎠4∂t∂rG€= r˙G∂t
この時，(8)，(9)式はそれぞれ以下のように変形した．
2Tr˙
pB = pv + pG −− 4µL GrGrG€
pG rG3 = const. (r˙G > 0)

rG

€

€

€

€

G G G

G G Gより，€r f˙r˙G = G G
3 fG
€(15)
rG ˙f˙G 2rG f˙G 2
˙r˙G =−3 fG
9 fG2と書ける．これを(7c)式に代入すると，すべり速度を無視して，
2rG2 1 ∂ 2 fG rG2 ⎛ 1 ∂fG ⎞ pB − pL(16)− ⎜⎟ =€3 fG ∂t 2 18 ⎝ fG ∂t ⎠ρLを得る．この式を簡単に考察すると，Rayleigh-Plesset の式と係数は異なるものの，形はほぼ同じであり，ボイド率の変化は気泡内外の圧力差（式の右辺）によることが分かる．ただし，左辺に気泡半径€の二乗が掛かっていることから，同じボイド率でも気泡径が大きいほど変化が小さい，言い換えると圧力差の影響が小さくなることが分かる．そこでここでは合体のモデル化として，合体すると気泡径が大きくなることから，その分，気泡の成長（や収縮）が（合体しない場合に比して）遅くなると考えた．具体的には，(7c-1)式および(7c-2)式の圧力差の項に，気泡半径の二乗に反比例する係数を掛ける．しかし，そのときの気泡径がいくらであるかは，合体の度合いによる．そこで，最も簡便な仮定として，初期ボイド率で気泡数
密度 nG0 ，ボイド率１では気泡数密度 0（つまりすべての気泡が合体）とし，その間の気泡数はボイド率で線形補間するとすれば，1− fG(17)nG = nG0
1− fG 0となる．ボイド率を固定すると，1 rG2 ∝ nG2 3 であるので，これから，(7c-1)及び(7c-2)式の圧力差の項を，

€

(7c)

2pB − pL
ρLp − p ⎛ 1− fG ⎞ 3⇒€ B L ⎜⎟ρ L ⎝ 1− fG0 ⎠と置き換えた．これにより，真に合体を扱った訳ではないが，合体により気泡の成長（や収縮）が遅れる現象を模擬できると考える．

(7d)€３．計算結果
まず，主な計算条件を以下に示す．
Table 1 Computational conditions
Clark-Y11.4%
NACA0015
No. of Grid Points
565×82501×60
Angle of Attack
2°8°
Chord Length10cm15cm
Flow Velocity
10m/s
8m/s
Reynolds Number
6.0×105
1.2×106
Time Step
1×10-7s
1.875×10-7s
6No. of Iterations
1×102×106
Initial Bubble Radius
10µm
Initial Void Fraction
0.10%
Vapor Pressure
2.3×103Pa
Surface Tension
7.2×10-2N/m
計算格子は C 型で，翼面上にそれぞれ 363 点および 301 点とっている（図１）．計算法の詳細は文献(6)-(8)と同じであるためここでは略すが，(1)式以下を一般座標変換し，対流項は固有値の最大値

(8)’

(9)’
pG rG3κ = const. (r˙G < 0)€
このまま解けば従来法と同じであるが，文献(8)では，ボイド率と気泡の膨張／収縮により，(7c)式を以下の式に置き換えるとした．
2∂ 2 r 3 ⎛ ∂r ⎞ p − p 1   2rG € 2G + ⎜ G ⎟ = B L + uG − u L , pB ≤ pL or fG < 0.3
ρL2 ⎝ ∂t ⎠
4∂t(7c-1)
23 ⎛ ∂rG ⎞ pB − pL(7c-2), pB > pL and 0.3 ≤ fG ≤ 0.9⎜ ⎟ =2 ⎝ ∂t ⎠ρL∂rG4(7c-3)= 0, rG = 3 0.9 πnG , pB > pL and fG > 0.9∂t3
これらの根拠は文献(8)に譲るが，膨張時にボイド率が上がった場合，気泡半径の２回微分項を落とす，さらに上限値を超えた場合，ボイド率が上がらないように気泡径を固定するというものである．
2Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 C02-3
で人工粘性項を作る Lax-Friedrich 的な風上差分，その他は中心差分で離散化し，(10)式を各時間ステップで収束させるため内部反復を行うのに合わせて全体の時間進行も陰的に行っている．

Fig. 5 Pressure distribution around Clark-Y 11.4% (σ=0.8)(a) Clark-Y 11.4%(b) NACA0015
Fig. 1 Close-up view of grids
0.6

0.5

0.4

experimentpreviousref. 8PresentCL0.3

Fig. 6 Pressure distributions of s=0.7 (blue), 0.8 (red) and 1.0 (black)
図３にキャビテーション数 0.4 のある時刻でのキャビティ形状とキャビティ後端付近の速度ベクトルを示す．キャビティはボイド
率 30%以上が白，それ以下がグレースケールで表されている．翼の上面にシート状のキャビティが長く伸びていること，キャビティの後端に逆流が生じていることが分かる．図４には翼近傍と翼まわりの圧力分布を示す．キャビティのところではほぼ圧力が一定になっていることが分かる．比較のため，キャビテーション数0.8 の結果を図５に示す．
画像からキャビティの存在は確認できな
いが，右のグラフから分かるように，上面の流れが加速するあたりで圧力が頭打ちになっており，実際この辺りでボイド率 10%弱の領域が確認できる．
実験との差異は主にキャビテーション数が0.6から0.8あたりである．そこで，キャビテーション数が 0.7，0.8 および 1.0 の場合の翼まわりの圧力分布を示したのが図６である．ここでは，単にキャビテーションが起こる領域で圧力が上がっているだけで，キャビティが伸びて上面の揚力を増やすといった現象は見られていない．実験を正しいとすれば，キャビティの成長を促進する何か別のメカニズムが必要であるように思われる．

0.2
0.4 0.6 0.811.2 1.4 1.6 1.82σFig. 2 Lift coefficient vs. cavitation index for Clark-Y 11.4%
まず，図２に Clark-Y 11.4%の場合の，キャビテーション数に対する平均の揚力係数を示す．実験は文献(11)，previous は文献(6)以前の著者らの初期の結果を示す．流れ場は非定常であるので，揚力係数の履歴を見て適当なところを平均した値である．実験と差はあるものの，これまでで最もよい結果を得ている．

0.90.8Fig. 3 Developed cavity on Clark-Y 11.4% (σ=0.4)

0.7

CL0.6exp.
Previousref. 7
Present

0.5
0.4
11.5

22.53σFig. 7 Lift coefficient vs. cavitation index for NACA0015
Fig. 4 Pressure distribution around Clark-Y 11.4% (σ=0.4)
次に，NACA0015 まわりの流れについて示す．まず，図７はキ3Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 C02-3
ャビテーション数に対する揚力係数のグラフである．これも実験(12)との一致はよくないが，これまででは最もよい結果を示している．

Fig. 8 Time sequence of pressure field for NACA0015 in every
1.875×10-2s (σ=1.0)
NACA0015 まわりの流れは非定常性が強いので，図８に 0.01875sおきの時系列の圧力分布を示す．ここに示した図やその間の変化を見ると，渦を放出した後，キャビティが伸びたり縮んだりを繰り返し，大きくキャビティが縮むと流れが剥離して大きな渦を放出しているように見える．渦放出の周期は 6Hz 程度であるが，条件は異なるものの，実験(12)では約 20Hz と異なっている．
４．おわりに	 気泡モデルに基づくキャビテーション流れ解析について，ボイド率が高くなったときに気泡の合体の影響を考慮したモデルを作り，２次元翼まわりの流れ場に適用した．同種のモデルの中では揚力係数等の比較でこれまでで最もよい結果を得たが，まだ実験とは隔たりがある．今後の課題としては，まだ考慮していないものとして，気泡の初生（今は初めから気泡核が一様に存在するとしている）と，気泡成長時の溶存気体の析出が挙げられる．これ
らを考慮することで結果がどう変わるかを調べたい．それとは別に，本手法（の従来法）は既に３次元に拡張し，実問題に適用している(13)ので，本モデルをこれらに適用していきたい．
参考文献	 
(1) Kunz, R. F., et al., “A Preconditioned Navier-Storkes Method for
Two Phase Flows with Application to Cavitation
Prediction, ”Computers and Fluids, 29 (2000) pp. 849-875.
(2) Senocak, I. and Shyy, W., “A Pressure-Based Method for Turbulent
Cavitating Flow Computations,” Journal of Computational Physics,
176 (2002) pp. 363-383.
(3) Dumont, N., et al., “Numerical Simulation of Cavitating Flows inDiesel Injectors by a Homogeneous Equilibrium ModelingApproach,” 4th International Symposium on Cavitaion (2001).
(4) Reboud, J. L. and Delannoy, Y., “Two-Phase Flow Modeling ofUnsteady Cavitation,” 2nd International Symposium on Cavitation(1994) pp. 39-44.
(5) Schnerr, G. H. and Sauer, J., “Physical and numerical modeling ofunsteady caviation dynamics,” 4th international Conference onMultiphase Flow (2001).
(6) Tamura, Y. and Matsumoto, Y., “Improvement of Bubble Modelfor Cavitating Flow Simulations,” Journal of Hydrodynamics, Ser.
B, 21, 1 (2009) pp. 41-46.
(7) Tsurumi, N., Tamura, Y. and Matsumoto, Y., “Improvement ofBubble Model in High Void Fraction for Cavitating FlowSimulations,” Journal of Computational Science and Technology, 6,3 (2012) pp. 113-128.
(8) Tamura, Y., Tsurumi, N. and Matsumoto, Y., “FurtherImprovement of Bubble Model for Cavitating Flow Simulations,”
8th International Symposium on Cavitation (2012).
(9) Schiller, L., Z. Ver. Dtsch. Ing., Vol. 77 (1993) pp. 318.
(10) Sridhar, G. and Katz, J., “Drag and Lift Forces on Microscopic
4Copyright © 2012 by JSFM

第 26 回数値流体力学シンポジウム
講演番号 C02-3
Bubbles Entrained Vortex,” Phys. Fluids, 7(2) (1995) pp.389-399.
(11) Knapp, R.T. et al., Cavitation, Engineering Societies Monographs,
McGraw-Hill (1970) pp267-320.
(12) Sakoda, M. et al., “Mechanism of Cloud Cavitation Generation on a2-D Hydrofoil,” International Symposium on Cavitation (2001).
(13) Fukaya, M., Tamura, Y. and Matsumoto, Y., “Prediction ofCavitation Intensity and Erosion Area in Centrifugal Pump byUsing Cavitating Flow Simulation with Bubble Flow Model,”Journal of Fluid Science and Technology, 5, 2 (2010) pp. 305-316.

5Copyright © 2012 by JSFM