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

レーザー誘起液体ジェットの数値解析モデル
The Numerical Model for Laser-Induced Liquid Jet⃝

石川 大樹，東北大・航空，宮城県仙台市青葉区荒巻字青葉 6-3，E-mail: ishikawa@shock.ifs.tohoku.ac.jp
孫明宇，東北大・流体研，宮城県仙台市青葉区片平 2-1-1，E-mail: sun@fmail.ifs.tohoku.ac.jp

中川 敦寛，東北大・医，宮城県仙台市青葉区星陵町 1-1，E-mail: nakg neurosurg@yahoo.co.jp
冨永 悌二，東北大・医，宮城県仙台市青葉区星陵町 1-1，E-mail: tomi@nsg.med.tohoku.ac.jp

Daiki Ishikawa, Tohoku University, 6-3, Aramaki aza Aoba, Aoba-ku, Sendai, Miyagi 980-8578, Japan
Mingyu Sun, Tohoku University, 2-1-1, Katahira, Aoba-ku, Sendai, Miyagi 980-8577, Japan
Atsuhiro Nakagawa, Tohoku University, 1-1, Seiryo-machi, Aoba-ku, Sendai, Miyagi 980-8574, Japan
Teiji Tominaga, Tohoku University, 1-1, Seiryo-machi, Aoba-ku, Sendai, Miyagi 980-8574, Japan1.

諸言

Laser-induced liquid jet (LILJ) is a promising surgical knife because it reserves blood vessel and emits less heat.
LILJ is ejected at high speed by shock wave and expansion of a bubble caused by the pulsed laser focusing inwater. In practice, this underwater explosion brings elastic deformation of the thin and narrow tube as well. In
this study, we developed a numerical model for laser-induced liquid jet under the influence of wall elasticity andrealistic fluids with phase change, covering a wide range for temperature and pressure. Using these models, thephenomenon of laser breakdown in water which is the early stage of LILJ was investigated. In the simulation, it
was found that cavitation bubbles are formed near the acrylic resin tube wall by low pressure caused by the large
oscillation of the tube wall.
Water
Fiber
Tube
Pulse
Ho: YAG Laser

近年，従来の液体ジェットメスを改良した，レーザー誘起液体ジェット（LILJ）を手術用メスとして用いる技術が開発されている (1,2) ．液体ジェットメスの最大の特徴
は，熱損傷がなく，組織選択性が高いため，細血管や神経の温存が可能な点である．従来の液体ジェットに比べて LILJ は，微調整が可能なことや，噴出水量が微量なため，顕微鏡手術の妨げとなっていた飛沫や気泡の発生を軽減できることが明らかになっている．また，このデバイスは構造が単純なため，カテーテル内や内視鏡への導入も比較的容易にできると考えられている (1,3) ．このように，手術用メスとして LILJ を用いることは，臨床実験からその有用性が確認されており，現在，商業化が見込まれている (2) ．そのため，より安全で安定した手術を行うために，LILJ の現象を，数値解析を用いて予測・再現することが重要となっている．

Y Connector

Laser Irradiation

Laser-Induced Liquid Jet
Gas

LiquidUnderwater Explosion

Solid

Fig. 1: Schematic diagram of the laser-induced liquidjet (LILJ) system is shown. LILJ is generated by irradiating pulsed laser within tube ﬁlled with water.
本研究では，考案した弾性変形モデルと実在流体状態
方程式，水中爆発モデルを用いて，噴出装置の弾性変形と，高温高圧の相変化を考慮できる完全圧縮性固気液連成数値解析手法を構築した．また，この手法により，LILJ
およびその発生過程の解析を行っている．
本稿では，この手法と，LILJ 発生初期の現象である，
レーザーエネルギーの解放によって水から気泡が誘起され，高圧領域が伝播する現象の解析について，得られた知見と共に紹介する．

LILJ は 次 の よ う な 狭 空 間 内 の 極 限 現 象 で 生 じ る（Fig. 1）．(1) 噴出装置内に満たされた水に，レーザーを照射させることで水中爆発を起こす．(2) 水は相変化し，高温高圧の気泡が生成する．(3) 生じた衝撃波と気泡の膨張によりジェットが噴出される．LILJ の現象を数値解析により解明するためには，気泡の膨張とジェットの噴出という圧縮性流体と非圧縮性流体の特徴を併せ持つ流体を解析することと，水中爆発のような高温高圧かつ相変化を伴う現象を解析すること，また，カテーテルや内視鏡への導入に向けて，噴出装置の変形を考慮できるアルゴ
リズムの構築が必要となる．本研究グループでは，航空宇宙の分野で培った圧縮性スキームのノウハウを活かし，圧縮性流体と非圧縮性流体の特徴を併せ持つ流体の解析を行ってきた (4,5) ．また，筆者らは，薄肉構造体の弾性変形を解析可能な弾性変形モデルの構築を行い，初期気泡の存在をあらかじめ仮定するモデルを用いて LILJ の解析を行った (6,7) ．さらに，実在流体状態方程式を考案
し，水中爆発モデルを構築した．その結果，エネルギーを得た水が蒸気へと相変化し，気泡となる現象を数値解
析により再現した (8−10) ．

2.1 手法の概要
本研究の手法は Lagrange-Remap 法による (11) ．その概要を下記に箇条書きで記す．弾性変形モデル，水中爆発モデルについては 2. 2 章，2. 4 章でそれぞれ説明する．

1

c 2012 by JSFM
Copyright ⃝

完全圧縮性固気液三相の数値解析手法
気液二相の流体解析には Lagrange-Remap 法を用いた
気液二圧力二速度二流体モデルを用いる (6,11) ．この手法と弾性変形モデルを統合することで，固気液三相の解析モデルとする．本章では，初めに手法の概要を述べ，その後に，弾性変形モデルと実在流体状態方程式，水中爆発モデルについて説明する．

2.

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

なお，水・蒸気の状態方程式は後述の実在流体状態方程式（2. 3 章）を使用する．

Lagrange 座標系における質量保存則および運動量保存則は，二相を平均化した物理量を用いて，微分形式で式 (7), (8) のように表せる．

(0) 水中爆発モデルより初期条件を設定する．

Dρ̄+ ρ̄∇ū = 0DtDū+ ∇p̄ = 0ρ̄Dt(1) Euler 座標系における格子を Lagrange 座標系の粒子の初期形状と仮定し，二相の物理量を平均化する．
(2) Lagrange 座標系の保存則および弾性変形モデルの基礎方程式を有限体積法で離散化し，粒子界面および弾性壁の圧力，速度を陰解法で計算する．

(7)
(8)

ここで，状態変化が等エントロピー的に行われると仮定
し，δ p̄ = ā2 δ ρ̄ を用いる．これらの式を連立させ，有限体積法を用いて離散化すると，式 (9), (10) のように記述できる．

(3) Lagrange 座標系の粒子の保存量を計算する．
(4) 粒子の体積変化，運動量変化，エネルギー変化を計算し，二相に分配する．

(6) Subgrid モデル (13) を適用する．

{}∑
p̄n+1 − p̄ni=− β (I¯in+1 )2(ū∗n+1 · nSij )
mi i
∆t
{
}
∑
− (1 − β) (I¯in )2(ū∗n · nSij ) (9)(7) (6) で求まった Lagrange 座標系の保存量を Euler 座
標系に再びマッピングする．

mi(5) キャビテーション・沸騰モデル(12)を適用する．

(1)-(7) を繰り返し計算する．
Lagrange 座標系における保存則は，流体に働く体積力を考慮しないとき，物質微分 D/Dt を用いて式 (1) のよ
うに表せる．
∫ID
U dΩ =
F · ndS(1)
Dt Ω(t)S(t)ここで，Ω(t) は検査体積，S(t) は検査体積の境界面，nは検査体積に対して外向きの単位法線ベクトルである．
また，非粘性流体の場合，密度 ρ，圧力 p，速度ベクト
ル u，比内部エネルギー e とすると，U と F は以下のように表せる．




ρ
0



, F = 
ρu
U =
(2)
 −pI 1−puρ (e + u · u)2
ただし，I は単位テンソルである．
Lagrange 座標系における保存則を解くための行列計
算に用いる，気液二相の平均値は以下のように定義する．
ρ̄，p̄，ā，ū はそれぞれ，密度，圧力，音速，流速の平均値である．ただし添え字の g と l はそれぞれ気相と液相
を表す．

ρ̄ =

ρg Ωg + ρl ΩlΩg + Ωlρg a2gρl a2l
pg +
pl
Ωl
Ωg
p̄ =
ρg a2gρl a2l+ΩgΩlv
uΩg + Ωl)ā = u
u (
Ωg
Ωl
t
ρ̄
+
ρg a2gρl a2lū =ρg ug Ωg + ρl ul Ωlρg Ωg + ρl Ωl

(3)

∑
ūn+1
− ūnii=− β(p̄∗n+1 nSij )∆t∑− (1 − β)(p̄∗n nSij )ここで，添え字 i, n はそれぞれ，セル番号，時間ステッ
プであり，p∗ , u∗ はそれぞれ，検査体積界面の圧力と検査体積界面の速度である．弾性壁界面においては，界面
法線方向速度 u∗ · n の代わりに，壁面法線方向速度 uw
を用いて弾性変形モデルと連成させる．なお，n ステッ
プ目の検査体積界面の圧力と速度は，MUSCL 法 を用いて空間二次精度とする．また，I, m, n, Sij はそれぞれ，
音響インピーダンス ρa，検査体積の質量， 検査体積界面に対する単位法線ベクトル，検査体積 i と j の境界表
面積である．β は時間精度を決める定数であり，本研究では後退オイラー法を採用する．また，音響インピーダ
ンス I¯in+1 は，I¯in で近似する．
検査体積界面の圧力と速度は，近似リーマン解より求める (14) ．界面を挟んで左側の検査体積を L，右側を Rの添え字で表すと，検査体積界面の圧力 p∗ は，検査体積
界面の法線方向速度 u∗n （u∗ · n または uw ） を用いて，式 (11), (12) のように表せる．

p∗ = pL − IL (u∗n − uL )∗p

= pR + IR (u∗n − uR )

(11)
(12)

内部セルの界面は，低マッハ数流れに対応するため，全速度修正を考慮した近似リーマン解を用いる．このとき検査体積界面の法線方向速度 u∗n と圧力 p∗ は，式 (13),(14) のようになる．

IL uL + IR uR + (pL − pR )β1
IL + IR
IL pR + IR pL + IL IR (uL − uR )β2p∗ =IL + IR

u∗n =
(4)

(10)

(13)
(14)

低マッハ数流れに対応するため，流れの最大マッハ数 M
に対して，M ≥ 1 のときは，β1 = β2 = 1，M < 1 のときは，β2 = M (2 − M ), β1 = β2−1 に設定する (15) ．一方，壁面の圧力は式 (11) または式 (12) より求める．これら
界面の式を式 (9), (10) に代入し，粒子界面および弾性壁の圧力，速度を陰解法（不完全 LU 分解，GMRES 反復法 (16) ）で計算する (11) ．

(5)

(6)

2

c 2012 by JSFM
Copyright ⃝第 26 回数値流体力学シンポジウム
講演番号 C04-5

行列計算後，Lagrange 座標系で求めた値を用いて，Lagrange 座標系での粒子の保存量を計算する．次に，粒子の体積変化，運動量変化，エネルギー変化を計算した後，保存則を満たすように保存量を二相に分配する．そして，キャビテーション・沸騰モデルおよび Subgrid モデルを適用させる．キャビテーション・沸騰モデルには，Merkle
らの非平衡モデルを二流体に適応させたもの (12) を用いる．また格子内の相互作用を，Subgrid モデルを用いて計算する (13) ．
最後に，求まった保存量を，二相それぞれ，Euler 座
標系にマッピングする (6) ．粒子の速度から計算した流束を用いて，有限体積法により Euler 座標系での保存量を求めている．

2.2 弾性変形モデル
本節では，薄肉円管の半径方向の運動方程式から構築した，薄肉構造体の弾性変形モデルについて紹介する (6) ．
円管が，長軸方向に細長く，平面ひずみ状態にあると仮定すると，長軸方向の影響を考慮しないとき，半径 r 方向の変位 s(r, t) の運動方程式は，式 (15) となる (17) ．
()∂ 2 s(r, t)1∂σr (r, t)
−σ
(r,
=
t)
+
σ
(r,
t)
+
r
(15)
θ
r
∂t2
ρw r
∂r
ここで，ρw , σθ (r, t), σr (r, t) はそれぞれ，円管の密度，
円周方向の応力，半径方向の応力である．ただし，応力 σθ (r, t), σr (r, t) は変位 s(r, t) を用いて下記のように表すことができる．
)(
s(r, t)′ ∂s(r, t)′ σθ (r, t) = E
+ν
(16)

r
∂r
(
)

∂s(r, t)s(r, t) σr (r, t) = E ′
+ ν′
(17)
∂r
r
ここで，平面ひずみ状態を仮定しているので E ′ と ν ′ の値はヤング率 E とポアソン比 ν を用いて式 (18), (19)のようになる．

 E′ =(1 − ν)E(1 + ν)(1 − 2ν) ν′ = ν
1−ν

(19)

(21)

円管を薄肉と仮定し，式 (15) を rin から rout まで積分
すると，d を円管の肉厚として，運動方程式は式 (22) のようになる．

d2 s(t)
1E′=−
2 s(t) + ρ d (pin − pout )2dtρw rinw(22)

2
ここで，k1 = E ′ /(ρw rin)，k2 = 1/(ρw d) とおくと，
式 (22) は式 (23) と表すことができる．

d2 s(t)= −k1 s(t) + k2 (pin − pout )dt2

ここで，このモデルの内壁面の圧力と壁面法線方向速度において，p∗ は流体の弾性壁界面の圧力，uw は流体の壁面法線方向速度となる．このように壁面圧力と壁面速度を流体の基礎方程式と連成させて，行列計算により陰的に解く．
n + 1 ステップ目の壁面に接するセルの検査体積 Ω n+1
は，n ステップ目の検査体積 Ω n と，壁面変位 s の増減から求めることができる．このように，仮想の検査体積を用いて n + 1 ステップ目の計算をすることで，格子を移動させること無く，移動境界問題を解くことが可能となる．

2.3 水・蒸気の実在状態方程式
LILJ は高温高圧下の水中爆発により生じるため，広範囲の圧力と温度において相変化を計算可能な，水と蒸気の実在状態方程式が必要となる．しかしながら，過去の圧縮性流体の研究で用いられている簡易的な状態方程式や線形化された内部エネルギーの式 (5,18,19) では，レーザー照射直後の現象のような高温高圧かつ相変化を伴う現象を正確に解析することに不十分である．また，PROPATH
や各種データベース (20,21) には計算の適用範囲が存在することや，反復計算を含む数値解析では計算効率が下がるという問題がある．
そこで，水の実在状態方程式には，Tait 式を，フィッティングを用いて高温域まで拡張させた，修正 Tait 式を考案した (8−10) ．Tait 式は，等温過程においては高圧の計算に適しており，従来から衝撃波の解析などに広く用いられている (19,22) ．構築した修正 Tait 式を式 (26) に示す．なお，psat ，ρlsat はそれぞれ飽和蒸気圧，水の飽和密度である．また，n = 7.15 である．
· 水の実在状態方程式{(p(ρ, T ) = Kp (T )(18)

また，境界条件は，円管の内圧と外圧をそれぞれ pin , pout ,円管の内半径と外半径をそれぞれ rin , rout とすると，
式 (20), (21) のようになる．
{
σr (rin , t) = −pin(20)

σr (rout , t) = −pout

次に，内壁面での圧力を p∗ ，壁面法線方向速度を uwとし，式 (23) を式 (24), (25) のように離散化する．
 n+1
 uw − unw= −k1 sn+1 + k2 (p∗n+1 − pout ) (24)
∆t
 n+1
s
= sn + {βun+1+ (1 − β)unw }∆t
(25)
w

ρ
ρlsat (T ))n }− 1 + psat (T )(26)

修正 Tait 式は，飽和曲線周りに式の展開を行うこと
で，相変化の解析を可能にしてある．また， Kp は温度の関数であり，常温常圧下での密度を基準にデータベース(20) とフィッティングを行うことで求める（Tab. 1）．以上より，修正 Tait 式は，飽和曲線周りと常温常圧下での密度で，正確な計算が期待できる水の状態方程式となっている．Fig. 2 (a) に修正 Tait 式を示す．三重点から臨界点近傍で，飽和曲線と常温常圧下での密度以外の値も，データベースと良い一致を示していることがわかる．
一方，蒸気の実在状態方程式には，谷下によって提案されたビリアル型状態方程式 (23) を用いる（Fig. 2 (b)）．

· 蒸気の実在状態方程式

 p(ρ, T ) = ρRT (1 + ρB (T ))
)−2.7
(
)−8.4
(
T
T

− 437
 B (T ) = −0.668
100
100

(27)
(28)

(23)
3

c 2012 by JSFM
Copyright ⃝第 26 回数値流体力学シンポジウム
講演番号 C04-5

Tab. 1: Fitting coeﬃcients ai and corresponding correlation coeﬃcient R for water and steam in temperature range[273.16K-647.096K]
f (T )

a0

a1

a2

a3

a4

a5

a6

R

ln(psat )−7.045364 × 10 7.037878 × 10−1 −2.716884 × 10−3 6.064245 × 10−6 −7.928070 × 10−9 5.647571 × 10−12 −1.692521 × 10−15

1.000

ρlsat

−3.726770 × 103

6.675586 × 10 −3.900921 × 10−1 1.217122 × 10−3 −2.146343 × 10−6 2.018627 × 10−9 −7.918682 × 10−130.993

ln(ρgsat )−1.989760 × 10 −2.089313 × 10−1 2.821146 × 10−3 −1.168679 × 10−5 2.365039 × 10−8 −2.390330 × 10−11 9.685766 × 10−15

0.987

elsat

1.232021 × 10

egsat

−5.006823 × 107 7.730674 × 105

Kp

−10

0.999

−4.691108 × 103

1.497363 × 10 −2.647511 × 10−2 2.459928 × 10−5 −9.396962 × 10−90.988

−1.318904 × 1010 1.648058 × 108−8.279949 × 105

2.210894 × 103

2.670858 × 10−3 −8.888164 × 10−7

0.988

fl1

−7.246856 × 10

52

−6.713922 × 10

fl2

−6.017620 × 107 7.259707 × 105−3.651604 × 103

fg−2.213051 × 106 2.195612 × 104−3.329504 × 10

64

1.083304 × 10

0.7

Pressure [GPa]

0.6
0.5
0.4
0.3

2.429734 × 10

2

−1

−8.280355 × 10

−3.328391−3−1.560331 × 10

−6−4.055142 × 10

0.992

−1.488568 × 10−2 1.206004 × 10−5 −4.082380 × 10−90.984

−9.091751 × 10 1.987885 × 10−1 −2.399263 × 10−4 1.499730 × 10−7 −3.731245 × 10−11

1.000

3.925001 × 10

900
850

1.6
1.4
1.2
1.0
0.8

Databaseρ[kgm-3] Saturation
999.975 (Fitted)
950
900
850
800
750
700
650
600
Saturation

0.4
0.2
0.0

300

400

500

600

Temperature [K](a) Water
3.0

Analysis
ρ[kgm-3]

Databaseρ[kgm-3]
0.5976
6
20
40
60
80
100
Saturation

Analysis
ρ[kgm-3]
0.5976(Fitted)
6
20

Databaseρ[kgm-3]
2.9 0.5976
6
2.8 204060
2.7 80100
2.6 Saturation

Specific Internal Energy [MJkg-1]

Pressure [MPa]

10

80
60

40

8

40
60
80
100

2.5

6

20

2.4

4
Saturation(Fitted)

2
0

300

400

600

Saturation(Fitted)

2.3

6
500

Analysis
ρ[kgm-3]
600
650
700
750
800
850
900
(Fitted)950
999.975(Fitted)

0.6(a) Water

12−1.566144 × 10

1.8

950

18

14

6.444526 × 10

9.818312

800
750
0.1
700
650
600
0.0
Saturation
300 400 500 600
Temperature [K] (Fitted)

16−62.208302

999.975(Fitted)

0.2

−3

−9

Analysis
ρ[kgm-3]

Databaseρ[kgm-3]
999.975
950
900
850
800
750
700
650
600
Saturation

1.565437 × 10

Specific Internal Energy [MJkg-1]

6

2.2

0.5976

300

400

500

600

Temperature [K]Temperature [K](b) Steam(b) SteamFig. 2: Comparison between analytical equations ofstates (lines) and database(20) (dots) of p-T diagrams
for water and steam

Fig. 3: Comparison between analytical expressions(lines) and database(20) (dots) of e-T diagrams for waterand steam4

c 2012 by JSFM
Copyright ⃝第 26 回数値流体力学シンポジウム
講演番号 C04-5

Temperature [K]また，内部エネルギーの式もデータベース (20) とのフィッティングにより構築する．水と蒸気の比内部エネル
ギーの式は，式 (29) のような飽和曲線周りの展開式から
求める．
( )
∂e
∆v
(29)
e(v, T ) = esat (T ) +∂v T
ここで，esat (T ), v, ∆v はそれぞれ，水・蒸気の飽和比
内部エネルギー，比体積，比体積の変化量である．熱力学より，右辺第二項における勾配は，式 (30) のように表せる (24) ．
( )
(
)
∂e
∂p
= −p + T
(30)
∂v T
∂T v

e(ρ, T ) = elsat (T )+(ρ− ρlsat (T )(31)

· 蒸気の比内部エネルギーe(ρ, T ) = egsat (T ) + (ρ − ρgsat (T ))fg (T )6∑

ai T i0.2 0.4 0.6 0.8 1.0
Energy density [GJm-3]

0.6
Pressure [GPa]

0.4
0.2
0.0
0.0

0.2 0.4 0.6 0.8 1.0
Energy density [GJm-3]

(b)

(32)

ここで，水の比内部エネルギーの式における fla ，flb は，
単純な形式でフィットすることは難しい (8) ．そこで，水の比内部エネルギーの式を改めて e(ρ1,2 , T ) = elsat (T ) +fl1,2 (ρ1,2 , T ) と定義し，データベースとのフィッティング
により fl1 ，fl2 を求めた後，fla ，flb を導出する手法を採っている．なお，ρ1 ，ρ2 は，常温常圧下での密度である ρ1 = 999.975 kgm−3 と， ρ2 = 900 kgm−3 とする．
水の実在状態方程式，並びに水と蒸気の内部エネルギーの式に用いる飽和曲線と，それぞれの式を構成する温度の関数は，すべて式 (33) のような温度に関する 6 次の多項式を用いている．

f (T ) =300

Fig. 4: Analytical results of initial conditions after
point energy release in water are shown. Energy density Ea /Ω0 is needed to be given: (a) Temperature T1
calculated from Eq. (26), (31), (34), (35); (b) Pressure
p1 at T1 , ρ1 calculated from Eq. (26). These values canbe initial conditions for numerical methods in ﬂuids.

)fla (T )+ (ρ−2 − ρlsat (T )−2 )flb (T )

400

(a)

· 水の比内部エネルギー
n−2

500

0.0

比体積が密度の逆数であることを用いて，水に対しては式 (26) を，蒸気に対しては式 (27), (28) を式 (30) に代入し，水と蒸気の飽和密度 ρlsat , ρgsat を基準に，それぞれ密度で積分する．その結果，飽和曲線 esat （水は elsat ，
蒸気は egsat ）からの比内部エネルギーの変化量が，水と蒸気でそれぞれ求まる．求まった変化量を，式 (29) に代入することで，比内部エネルギーの式を導出する．この後，温度 T の関数 f (T ) をデータベース (20) とフィットする．導出した比内部エネルギーの式を以下に列挙すると共に，Fig. 3 に示す．

n−2

600

(33)

i=0

フィッティング係数 ai を，三重点から臨界点までの相関係数 R と共に Tab. 1 に示す．これらすべての関数において高い相関が得られている．

2.4 水中爆発モデル
本節では，考案した水中爆発モデルについて紹介する (8−10) ．水中爆発の数値解析は，高エネルギーの水を流体解析の初期条件として与え，流れの計算の中で相変化を起こし，気泡を生成させることで行う．このモデル
は，レーザー照射初期の水中爆発の過程を次のように仮定している．まず，レーザーのエネルギーが微小体積 Ω0で水に吸収され，その領域が飽和状態になったとする．次に，その領域が等温過程を経て初期の密度 ρ0 と等しくなったとする．これは，化学反応が無視できるとき，水の質量は一定と考えられるからである．このとき，質量保存則とエネルギー保存則はそれぞれ式 (34), (35) となる．
{
ρ0 Ω0 = ρ1 Ω1(34)

ρ0 Ω0 e0 (ρ0 , T0 ) + Ea = ρ1 Ω1 e1 (ρ1 , T1 )(35)

なお，初期の温度，密度，比内部エネルギー，微小体積をそれぞれ，T0 , ρ0 , e0 , Ω0 ，圧縮後を，T1 , ρ1 = ρ0 , e1 , Ω1 ，
水がレーザーから吸収したエネルギーを Ea とする．エ
ネルギー密度 Ea /Ω0 をパラメータとして与え，式 (26),(31), (34), (35) を連立させることによって，未知数が T1のみの式となる．したがって，反復法を用いて T1 を求めることができる．求めた T1 と，既知である ρ1 を用いることで，水のすべての物性値を求めることができる．
Fig. 4 に，このモデルで計算した，エネルギー密度に対する T1 と p1 の値を示す．ただし，T0 = 293.15 K，
p0 = 0.101325 MPa とする．これらの物性値は，流体解
析の初期条件となる．このモデルは，エネルギー密度をパラメータと定義しているため，数値解析の次元を任意に選ぶことができる．
解析条件
レーザーにより気泡が誘起される LILJ 発生初期の現象について，高温高圧の相変化を考慮した，完全圧縮性固気液二次元解析を行う．LILJ 発生装置は，径 1.0 mm の細管の中央に，径 0.4 mm のファイバーが配置され，ファイバー先端から出射端までの距離（standoﬀ distance）が10.0 mm の，神経膠腫用デバイスを考える．

3.

5c 2012 by JSFM
Copyright ⃝第 26 回数値流体力学シンポジウム
講演番号 C04-5

Stainless Steel

Acrylic Resin

Young’s Modulus[GPa]

200

3.2

Poisson’s Ratio [-]

3.0

3.5

7800

1200

Density

[kgm

−3

]

A Formed Bubble

初期条件は次のように設定する．細管中を T0 =293.15 K，p0 = 0.101325 MPa の水で満たす．そこへ，
ファイバー先端より，レーザーエネルギー 130 mJ (2) を水に与える．水は，ファイバー先端面に沿った径 0.4 mmの半円形領域内（奥行 0.2 mm）で，レーザーエネルギー
を 0.1(25)（13 mJ）瞬時に吸収するものと仮定する．このとき，この範囲の水のエネルギー密度は 1.0345 GJm−3であるから，水中爆発モデルより (Fig. 4)，温度 T1 =598.86 K，圧力 p1 = 0.57384 GPa が，高温高圧水の初期条件として得られる．
管壁の境界は弾性壁とする．弾性壁の物性値には，カテーテルや内視鏡を想定し，ステンレス鋼とアクリル樹
脂（ポリメチルメタクリレート）の値を用いる (Tab. 2)．
なお，ファイバーの境界は剛体壁とする．

4.

解析結果・考察

Liquid Jet

Fiber

Acrylic Resin Tube
1

エネルギーを得た水は液体ジェット放出方向にそのエネルギーを解放し，水中爆発を起こす．そのとき，爆発中
心部では圧力が下がり，ファイバー先端の中央から，初期気泡が生成し始める．その後，気泡は膨張し，高圧領域は液体ジェット放出方向へと伝播する．
ステンレス管とアクリル管を比較すると，剛性が小さ
いアクリル管の方が管壁は大きく振動し，高圧領域の伝播が遅れることがわかる．また，アクリル管壁面では，初期気泡とは別に，管の変形により生じた低圧領域にキャビテーションが生じることがわかる．一方，剛性が大きいステンレス管の場合，管壁近傍のキャビテーションは確認されない．したがって，剛性の小さい材料をカテー
テルや内視鏡に用いた場合，エネルギー照射直後に，その壁面近傍でキャビテーションが生じる可能性が示唆される．

-0

Liquid Jet
Cavitation
0.0

Fig. 5 に LILJ 発生初期における，初期気泡生成部付
近の二次元解析結果を，管の材質がステンレス鋼のものとアクリル樹脂のものを並べて示す（t = 1.7 µs）．図には，ボイド率（二相のうち気相が占める割合）の分布と，圧力（実線），速度ベクトル（矢印）を示す．なお，圧力はゲージ圧である．

5.

t=1.7 μs

Stainless Steel Tube

Tab. 2: Physical properties of the device

1.0

Void Fraction [-] 0.001
Pressure [MPa] 0.0

2.0
0.01
5.0

10.0

3.0
[mm]

0.1

1.0

15.0

20.0

Fig. 5: Numerical results of laser breakdown in water which is the early stage of LILJ at t = 1.7µs areshown. Cavitation bubbles are formed near the acrylic
resin tube wall by low pressure caused by the large oscillation of the tube wall.

参考文献(1) T. Hirano, M. Komatsu, T. Saeki, H. Uenohara, A.
Takahashi, K. Takayama and T. Yoshimoto, “Enhancement of Fibrinolytics with a Laser-Induced
Liquid Jet”, Lasers. Surg. Med., Vol. 29 (2001),
pp. 360-368.

(2) 中川敦寛, 隈部俊宏, 小川欣一, 平野孝幸, 金森政之,
斎藤竜太, 渡辺みか, 橋本時忠, 中野徹, 亀井尚, 上之
原宏司, 高山和喜, 冨永悌二, “パルス Ho: YAG レー
ザーを用いた微小ジェットの原理開発と臨床応用：
東北大学における医工産学連携の取り組み”, 日レ医誌, Vol. 30, No. 2 (2009), pp. 119-125.

結言

本研究では，LILJ 噴出装置の弾性変形と，高温高圧の相変化を考慮できる完全圧縮性固気液連成数値解析手法を構築した．LILJ 発生初期の現象である，レーザーエネ
ルギーの解放によって水から気泡が誘起され，高圧領域が伝播する現象について解析を行った．その結果，アクリル樹脂を用いた LILJ 発生装置では，レーザーエネル
ギー照射直後，水中爆発と共に管は大きく振動し，低圧となった壁面近傍でキャビテーションが生じることが確認された．

(3) T. Ohki, A. Nakagawa, T. Hirano, T. Hashimoto,V. Menezes, H. Jokura, H. Uenohara, Y. Soto, T.
Saito, R. Shirane and T. Tominaga, “Experimental Application of Pulsed Ho: YAG Laser-Induced
Liquid Jet as a Novel Device for Rigid Neuroendoscope”, Lasers. Surg. Med., Vol. 34 (2004), pp. 227234.
6

c 2012 by JSFM
Copyright ⃝第 26 回数値流体力学シンポジウム
講演番号 C04-5(4) M. Sun, “Numerical and Experimental Studies ofShock Wave Interaction with Bodies”, Ph. D. Thesis, Tohoku University (1998), pp. 1-163.

(19) A. Chernishev, N. Petrov and A. Schmidt, “Numerical Investigation of Processes Accompanying
Energy Release in Water Near the Free Surface”,
Proc. of the 28th ISSW, Vol. 2 (2012), pp. 3-8.

(5) 菊池崇将, 孫明宇, “レーザー誘起液体ジェットの数
値計算”, 衝撃波シンポジウム講演論文集, Vol. 2007
(2007), pp. 463-464.

(20) W. Wagner, Description of the Software (Dynamic
Link Library) for the Calculation of Thermodynamic and Transport Properties of Water (IAPWS95), Lehrstuhl fuer Thermodynamik, Ruhr University Bochum, Germany (2004).

(6) 石川大樹, 孫明宇, 中川敦寛, 冨永悌二, “管壁の微小弾性変形を考慮したレーザー誘起パルス液体ジェットの数値解析モデルの開発及び液体ジェットの特性評
価”, 日本機械学会論文集（B 編）, Vol. 78, No. 793(2012), pp. 1467-1482.

(21) A Program Package for Thermo-Physical Properties of Fluid by PROPATH Group, Research Institute for Information Technology Kyushu University.

(7) 石川大樹, 孫明宇, 中川敦寛, 冨永悌二, “管壁の弾性変形を考慮したレーザー誘起液体ジェットの数値解
析モデルの開発”, 数値流体力学シンポジウム講演論
文集, Vol. 25 (2011), A06-1, pp. 1-4.

(22) G. K. Batchelor, An Introduction to Fluid Dynamics (1967), pp. 55-57, Cambridge University Press.

(8) 石川大樹, 孫明宇, “相変化の数値解析における水・蒸気の実在状態方程式を用いた水中爆発の初期状態
設定モデル”, 混相流研究の進展（年会講演会論文精選集）, Vol. 26, No. 5 (2013), (投稿中).

(23) 中嶋宏, 近藤博美, “水蒸気の第 2 ビリアル係数に関する一考察”, 長崎大学工学部研究報告, Vol. 7 (1976),
pp. 21-25.
(24) 日本機械学会編, JSME テキストシリーズ 熱力学,
初版第 6 刷 (2007), pp. 89-97, 丸善.

(9) 石川大樹, 孫明宇, “水中爆発の数値解析における水・蒸気の実在状態方程式を用いた初期状態の設定モデル”, 日本混相流学会年会講演会講演論文集, Vol. 2012(2012), pp. 330-331.

(25) 小板丈敏, 孫明宇, “水中点爆発気泡現象の初期条件の設定に関する研究”, 衝撃波シンポジウム講演論文
集, Vol. 2010 (2011), pp. 191-194.

(10) D. Ishikawa and M. Sun, “The Numerical Setup
of Conditions for Point Energy Release in Water”, Proc. of the Ninth International Conferenceon Flow Dynamics (2012), pp. 552-553.
(11) 矢田和之, 市東素明, 孫明宇, “非構造格子を用いた
all-speed 流れ解析”, 数値流体力学シンポジウム講
演論文集, Vol. 24 (2010), C3-3, pp. 1-6.
(12) Y. Zhu and M. Sun, “Numerical Simulation of Cavitation with Two-Fluid Method”, Proc. of the 10th
International Symposium on Advanced Fluid Information and Transdisciplinary Fluid Integration(AFI/TFI-2010), PS1-7 (2010).
(13) 孫明宇, “Subgrid 粒子を含む気液二相流の直接シミュ
レーション”, 数値流体力学シンポジウム講演論文集,
Vol. 24 (2010), A8-2, pp. 1-6.
(14) E. F. Toro, Riemann Solvers and Numerical Methods for Fluids Dynamics, Third Edition (1999),
pp. 293-304, Springer-Verlag.
(15) M. S. Liou, “A Sequel to AUSM, Part II: AUSM+-up for All Speeds”, J. of Computational Physics,
Vol. 214 (2006), pp. 137-170.
(16) Y. Saad and M. H. Schultz, “GMRES: A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems”, SIAM J. of Sci. Stat.
Comput., Vol. 7 (1986), pp. 856-869.
(17) 高橋幸伯, 町田進, 角洋一, 基礎材料力学, 三訂第 4 刷(2009), pp. 154-160, 培風館.
(18) A. H. Koop, “Numerical Simulation of UnsteadyThree-Dimensional Sheet Cavitation”, Ph. D. Thesis, University of Twente (2008), pp. 39-51.
7c 2012 by JSFM
Copyright ⃝