第 25 回数値流体力学シンポジウム
C07-3

水中放電により誘起される強非定常なセカンドキャビテーションの数値解析
Numerical Simulation of Strongly Unsteady Secondary Cavitation Induced
by Underwater Electric Discharge○小板丈敏, 東北大院, 宮城県仙台市青葉区荒巻字青葉 6-3, E-mail: kotia@iswi.cir.tohoku.ac.jp
孫明宇，東北大学際センター, 宮城県仙台市青葉区荒巻字青葉 6-3, E-mail: sun@cir.tohoku.ac.jp
Taketoshi Koita, Tohoku University, Aramaki aza Aoba 6-3, Aoba-ku, Sendai
Mingyu Sun, CIR, Tohoku University, Aramaki aza Aoba 6-3, Aoba-ku, Sendai
Strongly unsteady secondary cavitation induced by underwater high electric discharge voltage in a rectangular tube andin a container were studied numerically. Numerical simulation was carried out with a two-fluid model and a modified
Merkle's cavitation model. It is considered that the behavior of the secondary cavitation is related to the definition of
the sound speed for the two-phase media. Therefore, two numerical models based on different definitions of soundspeed for the cavitation region are tested. One is the classic homogeneous sound speed, and the other introduces amodification, in which the liquid sound speed is adopted if the vapor volume fraction is below certain threshold value.
These numerical results are compared with experimental result. It was found that numerical solution depends on thedefinition of the two-phase sound speed. The new defined sound speed with this modification does give a betteragreement.
Euler 座標にマッピングをし，次のステップでの初期値とした[7]．
ここで，セカンドキャビテーションに見られるように，セルより小さい気相の流体界面構造の近似計算を行うために，subgrid モデ
ル法を用いて[8]，セル内の気相，液相の体積，運動量及びエネルギー変化の相互作用を計算した．気相と液相の状態量には，それぞれ，理想状態方程式と Tait 型の状態方程式を用いた．また，キャビテーションの解析を行うため，相変化を考慮した．相変化モデルには Merkle’s キャビテーションモデルを二流体モデルに適合し，一流体モデルでの圧力を二流体の平均圧力とした[6]．総質量の変化がないと仮定し，ml を液相の質量とすると，二相流のキャビテーションモデルは以下に示される．

１．緒言
矩形管内部で水中放電により水中爆発が生じた際，水中爆発によって誘起された気泡の上の水は加速され，液体ジェットが形成される．そして，気泡と底壁面との間でセカンドキャビテーションが発生することが報告されている[1，2]．この液体ジェットは医学の分野でジェットメスの開発に応用され[3]，また，火山学において，実験室での火山噴火の挙動の基礎研究として応用されている[4]．ここで，このセカンドキャビテーションは管壁面でエロージョンを引き起こすと考えられる．また，容器内部でも水中爆発直後にキャビテーションの発生が観測されている[5]．よって，ジェットメスの開発や基礎研究に対して，セカンドキャビテーションの挙動の解明が重要であると言える．従って，本研究では矩形管内部および容器で発生するセカンドキャビテーションの挙動の数値解析を行う．ここで，キャビテーション現象の解析に対して，代表的に翼周りに発生する周期的なキャビテーション流れの解析およびキャビテーションモデルの開発が行われてきた．Y.
ZhuとM. Sun[6]は翼周りのキャビテーション流れ場で翼の揚力係数は音速伝播の時間のスケールにおいて緩和時間を伴って劇的に変化することを解析結果より示している．よって，数値解析において，キャビテーション現象は二相媒体の音速の定義に依存すると考えられる．また，本研究で解析を行うセカンドキャビテーシ
ョンは周期的に発生するのではなく，水中爆発の発生直後に生じる強非定常現象である．従って，このキャビテーションの数値解析において，二相媒体の音速の定義が重要であると考えられる．
以上より，本研究では矩形管内部および容器で発生するセカンドキャビテーションを解析対象とし，二相媒体の音速の定義の変化させて解析を行う．そして，得られたキャビテーションの体積率の解析値と実験値を比較をし，二相媒体の音速の定義がキャビテーションの発生に与える影響を評価する．

evaporatio n :∂ml= k E ⋅ ml ⋅ ( ~p − p st (T ))∂t(1)∂ml(2)= k C ⋅ ( M − ml ) ⋅ ( ~p − p st (T ))∂tkE，kC は相変化の定速度であり，CE = 1.0×103，CC = 1.0×103，流condensation :れ場の代表長さ L∞ を用いて，式(3),(4),(5)より求められる．
kE =CE
20.5 ρ ∞U ∞ t ∞(3)

kC =

CC
20.5 ρ ∞U ∞ t ∞(4)

L∞
U∞(5)

t∞ =

また，pst(T)は水温T [℃]での水の飽和蒸気圧を表し，以下のTetensの式を用いた．
7.5T  2+
237.3+T p st (T ) = 6.11 × 10 
２．計算手法
本研究では，，非粘性の圧縮性気液二圧力二速度二流体モデルを用いて計算した．また，非粘性圧縮性方程式は有限体積法を用いて離散化した．また，Lagrange-Remap 法を用いた半陰解法で計算を行った．この Lagrange-Remap 法では Lagrange ステップにおいて，Lagrange 座標での基礎方程式を ILU+GMRES 反復法で解き，次に Remap ステップとして Lagrange 座標で求められた保存量を

[ Pa]

(6)

セル内の平均圧力 ~p は液相と気相の圧力 p，密度ρ，体積 V に添字l，v を用いて，以下のように表される．
β ⋅ p + β l ⋅ pv
~p= v lβl + βv
βl =1

ρ l al2Vl,βv =(7)ρ v a v2(8)VvCopyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C07-3て，-140 mm < y < 15 mm は水とした．本計算では二相しか扱えないため，y > 15 mm は水蒸気とした．実験では y > 15 mm は大気であるが，y > 15 mm での空気と水蒸気の違いは気泡現象およびキャビテーションの発生にはあまり影響を及ぼさないと仮定した．
放電電圧は V = 4.5 kv，放電エネルギーへの変換率はα = 0.3 とし，放電エネルギーは Ev = 0.675 J と計算した[5]．気泡の初期条件は気泡半径 rb = 0.05 mm とし，また，上記の放電エネルギーで算出された気泡の内部圧力pv0 = 5421.2 MPa および内部温度Tv0 = 11756 Kを用いた[5]．水および 15 mm < y < 140 mm での水蒸気は圧力1013.25 hPa，温度 293.15 K とした．境界条件は壁は固体壁とし，y = 140 mm において開放境界とした．本解析の初期格子は 4×74であり，気泡と気液界面の密度変化に対し解適合格子法[10]を用いた．図 4(b)に計算開始後の気泡周りおよび水面近くの格子を示す．解適合格子法において，0 段階から 4 段階までの解適合を行い解析を行った．クーラン数は 0.2 とした．また，ε で定義された音速 a *がキャビテーションの発生および体積率に与える影響を調べるため，ε = 0.00 %，0.05 %，0.07 %，0.10 %，0.15 %，0.20 %と変化させて解析を行った．

式(2) が平衡状態に収束する際，気相の凝縮は止まり，平均圧力は~p = p st となる．よって，式(9)，式(10)が求められる．

∂mleee= − k E ⋅ f ( m l ) ⋅ ( m l − m l ) = −λ E ( m l − m l )∂t(9)∂mleee= − k C ⋅ g ( m l ) ⋅ ( m l − m l ) = −λ C ( ml − ml )∂t(10)

従って，平衡状態の液相の質量 ml は式(5),(6)を陽解法で解き，

ml = ml + (ml 0 − ml ) ⋅ exp( −λ ⋅ δt )ee(11)となる．
以上より，本研究では上記の平衡キャビテーションモデルを用いて，δ t が非常に大きいと仮定し，セル内の平均圧力が水の飽和蒸気圧 pst(T)に到達するまで反復計算を行った．
また，本研究で解析を行うキャビテーション現象は二相媒体の音速の定義に依存すると考えられ，セル内の音速の定義が重要であると言える．液相と気相を含むセルの音速には式(13)で示されるセル内で均質された Homogeneous 型の音速 a homogeneous が一般的に用いられている．ここで，矩形管および容器の大部分の領域が純水に満たされた試験部に対して，セカンドキャビテショーンの体積率はとても小さく，セル内に占める気相の割合は 0.1%以下であると考えられる[9]．よって，セル内には二相が存在するが，液相が支配的なセルが存在すると考えられる．従って，本研究ではLagrange ステップにおいて，以下のようにセルの音速 a *を定義した．

(σ > ε )aa* =  homogeneousal(σ ≤ ε )ahomogeneous =(12)(a)(b)
Fig.1 Mesh of initial grid, Maximum refinement level = 4(a) Overall initial grid, (b) Initial grid zoomed view of near the bubble andwater surface

Vl + Vv
 VlV + v 2 2ρvCv  ρl Clρ al =３．１．２ 解析結果および考察
図 2，図 3 にε = 0.05 %で定義された音速 a *を用いて解析を行った矩形管内部の水蒸気の体積率分布，および圧力分布を示す．放電により水中爆発が生じた際，爆轟波が発生し，水中を伝播する．
図 3 の圧力分布の時間履歴より，この爆轟波は底壁面で反射され，気泡と底壁面の間で負圧の領域が生じることが分かった．また，
図 2 の体積率分布の時間履歴において，気泡と底壁面の間でキャビテーションの発生が観測された．従って，このセカンドキャビテーションは底壁面から反射された爆轟波により生じる負圧によって誘起されたと考えられる．図４に各εの条件で得られたキャビテーション領域の界面速度の時間履歴の解析結果および実験値を示す．図４より，界面速度の時間履歴はε に依存し，Homogeneous型の音速 a homogeneous を用いた場合より，本研究で定義した音速 a *を使用して得られた解析結果の方が実験値と近い傾向を示すことが分かった．また，各条件εで解析された -80 mm < y < -20 mm での水蒸気体積率を 10 mm 間隔で平均化した体積率の最大値の時間履歴を図 5 に示す．図 6 に図 5 における時間 t = 4.0 ms までで得られた各解析条件εでの水蒸気の最大体積率を赤丸で示し，実験で同条件で得られた最大体積率[9]を青線で示す．図 5，6 より，キャビテーションの水蒸気体積率はε に依存し，ε が減少するにつれて，最大体積率が増加することが明らかになった．また，ε = 0.0 %の場合，すなわち，a *= a homogeneous であるとき，キャビテーションの体積率は実験値より高くなことが分かった．以上より，強非定常なセカンドキャビテーションの解析に対して，式(12)で定義したセルの音速 a *を用いた数値解析モデルは Homogeneous 型の音Vl V ρ  l 2  ρ l Cl 
σはセル内に占める水蒸気の体積率であり，εは水の音速を用いる際のセル内の水蒸気の体積率の限界値である．ここで，式(12)において，ε が 0 であるとき，a *は a homogeneous となる．
３．キャビテーション解析および a *の検証
矩形管内部で水中放電により生じるキャビテーションの発生は底壁面で反射した爆轟波により誘起される負圧領域で生じると考えられる[1]．また，容器内部で水中爆発直後に発生するキャビテーションは爆轟波背後に生じる負圧により誘起されると考えられる[5]．よって，本研究では矩形管および容器で発生する要因が異なる強非定常なキャビテーションを解析対象する．また，本研究で定義した音速a *およびHomogeneous 型の音速a homogeneous を用いた解析結果を実験結果と比較し，音速の定義がキャビテーションの発生に与える影響の検証を行う．
３．１ 矩形管内部で生じるセカンドキャビテーションの解析
３．１．１ 解析条件
図 1(a)に計算開始後の格子を示す．気泡は水蒸気とし，(x, y) = (0,0)の位置に設置した．本解析では水面からの深さ H = 15 mm において水中爆発が発生した際に生じる現象を解析対象とした．よっ2Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C07-3
速 a homogeneous を用いた解析より，有効であると考えられる．

Fig.5 Histories of difference maximum volume fraction obtained bynumerical model with sound speed definition(ε )Fig.2 Time variation of volume fraction in a tube at ε = 0.05 % fordischarge energy Ev = 0.675 J and explosion depth H = 15 mm withphase change

Fig.6 Dependence of maximum volume fraction on sound speeddefinition(ε )３．２ 容器内部で生じるセカンドキャビテーションの解析
３．２．１ 解析条件
図 7(a)に計算開始後の格子を示す．-150 mm < y < 150 mm は水とし，y > 150 mm は水蒸気とした． (x, y) = (0, 0)の位置に水蒸気の気泡を設置した．気泡の初期条件は気泡半径 rb = 0.05 mm とした．ここで，本解析では放電電圧は V = 4.5 kv において，放電エネルギーへの変換率をα = 0.2 およびα = 0.4 と変化させて解析を行った．よって，気泡の内部圧力および温度の初期条件は，それぞれの変換率αに対して，α = 0.2 のとき pv0 = 3358.6 MPa，Tv0 = 7283
K，および，α = 0.4 では pv0 = 7483.9 MPa，Tv0 = 16230 K と算出し，解析を行った[5]．
水およびy > 150 mm の水蒸気は圧力1013.25 hPa，温度 293.15 K とした．境界条件は壁は固体壁とし，y > 165 mm において開放境界とした．計算開始後の格子数は 17752 であり，0段階から 3 段階までの解適合格子法[10]を用いて解析をした．図7(b)に計算開始後の気泡周りの格子を示す．また，ε = 0.00 %，0.05 %，0.10 %，0.15 %，0.20 %と変化させ解析を行った．クーラ
ン数は 0.1 とした．

Fig.3 Time variation of pressure in a tube at ε = 0.05 % for dischargeenergy Ev = 0.675 J and explosion depth H = 15 mm with phase change

３．２. ２ 解析結果および考察
図 8，図 9 にα = 0.4 およびε = 0.05 %で定義された音速 a *を用いて解析を行った容器内部の t = 0.04 ms，0.08 ms における水蒸気の体積率分布，および圧力分布を示す．ここで，奥行き 5 mm の容Fig.4 Histories of cavitation interface velocity for y direction at x = 0.0mm obtained by numerical model with sound speed definition(ε )3Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C07-3(a) t = 0.04ms(b) t = 0.08m
Fig.8 Time variation of volume fraction in a container at ε = 0.05 % fordischarge energy Ev = 0.90 J with phase change(a)(b)
Fig.7 Mesh of initial grid, Maximum refinement level = 3(a) Overall initial grid, (b) Initial grid zoomed view of near the bubble
器内部で放電電圧 V = 4.5 kv で水中放電をした際，t = 0.04 ms 以降で-50 mm < x < 50 mm，-50 mm < y < 50 mm の範囲において水中爆発により膨張する気泡の周りにキャビテーションが発生することが報告されている[5]．しかし，図 8 の水蒸気の体積率分布図で見られるように，数値解析において，キャビテーションの発生は確認できなかった．ここで，α = 0.2，0.4 において，x = 0.0 mm における y 軸方向の全ε 解析条件に対する t = 0.04 ms，0.08 ms での圧力分布をぞれぞれ図 10，11 に示す．また，水温 293.15 K (T = 20 ℃)での水の飽和蒸気圧は式(6)より，pst(T) = 2336.56 Pa であり，この値を図 10，11 に赤線で示す．図 10，11 より，キャビテーションが生じる-50 mm < y < 50 mm において，全てのε の解析条件で圧力分布は水の飽和蒸気圧より高く，キャビテーションが発生していないことが分かった．ここで，報告された実験結果[5]では，容器はステンレス製であり，また，幅 10 mm の 2 枚のアクリル製の可視化窓で挟まれている．石川と孫ら[11]の壁面の弾性変形を考慮したレーザー誘起液体ジェットの数値解析モデルの検証結果において，高圧状態により壁面が弾性変形をする際，壁面の圧力は非定常的に負圧になることが報告されている．また，水中爆発により生じた衝撃波および膨張波はアクリル製の可視化窓に干渉し，透過波および反射波が生じると考えられる．よって，以上より，この透過波および反射波の影響および壁面の負圧によって，水中爆発で誘起された気泡の周囲の圧力が負圧になり，セカンドキャビテーションが発生すると考えられる．従って，容器内部で水中放電により生じる強非定常なセカンドキャビテーションの解析を行うためには可視化窓の透過波と反射波を考慮した 3 次元数値解析を行う必要があると言える．

(a) t = 0.04ms(b) t = 0.08m
Fig.9 Time variation of pressure in a container at ε = 0.05 % for dischargeenergy Ev = 0.90 J with phase change

Fig.10 Relationship of pressure and y axis at α = 0.4 and α = 0.2 fordifferent ε at t = 0.04 ms

４．結言
矩形管および容器内部で水中放電で誘起される非強定常なセカ
ンドキャビテーションの解析において，Merkle’s キャビテーシ
ョンモデルを二流体モデルに適合した相変化モデルを用いた．
そして，本研究では Homogeneous 型の音速 a homogeneous および，新たに定義した二相媒体の音速 a *を用いて解析を行った．矩形管内部で生じるセカンドキャビテーションの解析において，キャビテーションは底壁面で反射した爆轟波により誘起される負圧領域で生じることが分かった．また，本研究で定義した音速 a *を使用した解析モデルは Homogeneous 型の音速 a homogeneous を用いた解析より，キャビテーション発生の解析に対して有効であることが確認できた．容器内部で生じるセカンドキャビテーションの解析において，キャビテーションの発生は確認できなかった．ここで，容器内部で水中爆発が生じた際，衝撃波および膨張波が壁面に干渉して発生する透過波および反射波の影響により，容器内部が負圧になり，キャビテーションが発生すると考えられる．よって，容器内部で水中放電で発生する強非定常なセカンドキャビテーションの解析には，3 次元数値解析を行う必要があると言える．

Fig.11 Relationship of pressure and y axis at α = 0.4 and α = 0.2 fordifferent ε at t = 0.08 ms
4

Copyright © 2011 by JSFM

第 25 回数値流体力学シンポジウム
C07-3
参考文献
(1) Koita, T., Zhu, Y. and Sun, M., “Investigation of bubble collapseand water jet induced by underwater explosion in a rectangular tube,” 28thInternational Symposium on Shock Waves, (2011).
(2) Koita, T., Zhu, Y. and Sun, M., “The Behavior of Bubble and Water
Jet Induced by Underwater Explosion in a Rectangular Tube,” 49th AIAA
Aerospace Sciences Meeting, (2011).
(3) 中川敦寛, 隈部俊宏, 小川欣一, 平野孝幸, 金森政之, 斎藤竜太, 渡辺みか, 橋本時忠, 中野徹, 亀井尚, 上之原宏司, 高山和喜, 冨永悌二, “パルス Ho:YAG レーザーを用いた微小ジェットの原理開発と臨床応用：東北大学における医工産学連携の取り組み”, 日本レーザー医学会誌, Vol. 30, No 2 (2009), pp. 119-125.
(4) Kedrinskiy, V., “Hydrodynamic aspects of explosive eruptions ofvolcanoes: simulation problems,” Journal of Shock Waves,(2009),
Volume 18.
(5) 小板丈敏, 孫明宇, "水中点爆発気泡現象の初期条件の設定に関する研究," 平成 22 年度衝撃波シンポジウム講演論文集, (2010),
pp.191-194.
(6) Zhu, Y. and Sun, M., “Numerical Simulation of Cavitation withTwo-Fluid Method,” AFI/TFI-2010, (2010).
(7) 矢田和之, 市東素明, 孫明宇, "非構造格子を用いたall-speed流れ解析," 第24回数値流体力学シンポジウム, C3-3(2010).
(8) Sun, M., “Direct simulation of gas-liquid two-phase flows with
subgrid particles, ” 24th Computational Fluid Dynamics, (2010), A8-2.
(9) Koita, T., Hayashi, K. and Sun, M., “Secondary Cavitation Induced
by Underwater Electric Discharge in a Tube,” 8th International
Conference on Fluid Dynamics, (2011).
(10) Sun, M. and Takayama, T. , “Conservative smoothingon anadaptive ouaelrilateral grid,” J. Com. Phys, 150 (1999), 143-180.
(11) 石川大樹, 孫明宇, 中川敦寛, 冨永悌二, "管壁の弾性変形を考慮したレーザー誘起液体ジェットの数値解析モデルの開発," 第25回数値流体力学シンポジウム, A06-1(2011).

5Copyright © 2011 by JSFM