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

管壁の弾性変形を考慮したレーザー誘起液体ジェットの数値解析モデルの開発
Development of a Numerical Model for Laser-Induced Liquid Jet under the Influence of Wall Elasticity⃝

石川 大樹，東北大学工学研究科，宮城県仙台市青葉区荒巻字青葉 6-3 E-mail：ishikawa@iswi.cir.tohoku.ac.jp
孫 明宇，東北大学学際科学国際高等センター，宮城県仙台市青葉区荒巻字青葉 6-3 E-mail：sun@cir.tohoku.ac.jp
中川 敦寛，東北大学医学系研究科，宮城県仙台市青葉区星陵町 2-1
冨永 悌二，東北大学医学系研究科，宮城県仙台市青葉区星陵町 2-1

Daiki Ishikawa, Tohoku University, Aramaki aza Aoba 6-3,Aoba-ku,Sendai,Miyagi
Mingyu Sun, Tohoku University, Aramaki aza Aoba 6-3,Aoba-ku,Sendai,Miyagi
Atsuhiro Nakagawa, Tohoku University, Seiryo-machi 2-1,Aoba-ku,Sendai,Miyagi
Teiji Tominaga, Tohoku University, Seiryo-machi 2-1,Aoba-ku,Sendai,Miyagi
Laser-Induced Liquid Jet (LILJ) is expected to be widely used as a surgical knife because a liquid jet
reserves blood vessel and emits less heat. LILJ is ejected at high speed by shock wave and expansion
wave caused by the pulsed laser. Strong shock and expansion wave bring elastic deformation of the thintube wall. In this study, we try to develop a numerical model considering the eﬀect of wall elasticity. In
this paper, we ﬁrst propose a simple elastic model. The eﬀect of the wall elasticity is investigated andcompared in the simulation of an oscillating high pressure channel.
1.

緒言

Wall

近年，加圧された液体をノズルから高速で放出する技術を手術用ジェットメスとして，臨床応用することが期待されている．
液体ジェットを手術用メスとして臨床応用する際の利点は，温熱効果がないことと，血管温存が可能であるという点である．
一方，欠点は，装置のサイズや重量が過大であり，また，気泡や組織が飛散してしまうという点である．

Laser-Induced Liquid Jet

BubbleLiquidAtmosphereFiberFig. 1 Laser-Induced Liquid Jet

従来の液体ジェットを用いた手術用メスの利点を生かしつつ，欠点を解消する技術として，レーザー誘起液体ジェット（LILJ）（Fig. 1）を手術用メスとして用いる技術を開発してきた (1) ．
LILJ は，液体中でパルスレーザーを発振させ，内部の液体を一部気化させた後，気泡の膨張と，それに伴い生じる圧縮波により高速で噴出される液体ジェットのことである．以上のような原理でパルス液体ジェットを発生させることで，装置を小さくすることができる．

管壁の弾性変形を考慮した数値解析モデルの構築

2.

本章ではまず，弾性変形モデルについて紹介した後，計算手法について述べる．

2.1
2.1.1

LILJ を手術用ジェットメスとして用いることは臨床現場からもその有用性及び実用性が報告されており，市場化に向けて液体ジェットの特性を詳細に評価する必要があると考えられる．
しかしながら LILJ の現象は，気泡の膨張に代表されるような圧縮性流体の現象と，低マッハ数の非圧縮性流体の現象が同時に起こっているため，既存の解析技術では，解析が困難である．
そのため，本研究グループでは，LILJ の特性を数値解析で評価可能なアルゴリズムの開発を行ってきた (2) ．

弾性変形モデル
基礎方程式の導出

本研究では，円管が内圧と外圧を受けたときに生じる半径方向の微小振動の運動方程式を立て，弾性変形モデルの基礎方程式を算出することとした．円管が，軸方向に細長く，平面ひずみ状態にあると仮定すると，半径 r 方向の変位 s(r, t) の運動方程式は，式 (1) となる (3) ．

∂ 2 s(r, t)
1=
∂t2
ρw r

LILJ 発生装置の管内では，非常に強い圧縮波や膨張波が発生しているため，薄い材料で作られた管壁は弾性変形し，LILJの特性に大きな影響を与えると考えられる．そこで本研究では，管壁の弾性変形を考慮したモデルを考案し，本研究グループで開発した気液二相の数値解析技術と統合することで，固気液三相の現象を解析できるアルゴリズムの開発を行ってきた．更に，構築したアルゴリズムを用いて，LILJ の特性評価を行っている．

()∂σr (r, t)−σθ (r, t) + σr (r, t) + r∂r(1)

ここで，ρw ，σθ (r, t)，σr (r, t) はそれぞれ，円管の密度，円管円周方向の応力，円管半径方向の応力である．ただし，応力 σθ (r, t)，σr (r, t) は変位 s(r, t) を用いて式 (2),(3) のように表わせる．

)(
s(r, t)′ ∂s(r, t)′+ν σθ (r, t) = Er∂r()

 σr (r, t) = E ′ ∂s(r, t) + ν ′ s(r, t)∂rr

本稿では，本研究で構築した弾性変形モデルについて紹介し，その後，モデルの検証について報告する．

1(2)
(3)

c 2011 by JSFM
Copyright ⃝第 25 回数値流体力学シンポジウム
A06-1

Atmosphere (Outside)
Elastic wall

Liquid or Gas (Inside)phase Lphase R

Fig. 2 Elastic wall modelFig. 3 Calculation of pressure and velocity at interfaces ofcontrol volumes

なお，平面ひずみ状態を仮定しているので E ′ と ν ′ の値はヤ
ング率 E とポアソン比 ν を用いて式 (4),(5) のようになる．

 E′ =(1 − ν)E(1 + ν)(1 − 2ν) ν′ = ν1−νここで，i はセル番号，p∗ ,u∗ はそれぞれ，検査体積界面の圧力，二次元直交座標系における検査体積界面の速度ベクトルを表す．また，Ii ，mi ，n，Sij はそれぞれ，音響インピーダン
ス，i 番目の検査体積内の質量， 検査体積界面に対する単位法
線ベクトル，検査体積の境界線の長さである．

(4)(5)

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

σr (rout , t) = −pout

2.2.2 検査体積界面の圧力と速度の算出
次に，音響ソルバーを用いて，リーマン問題を近似的に解く
ことにより，検査体積界面の圧力と速度を求める (4) ．Fig. 3 のように，壁面の垂直方向を x 軸とし，検査体積界面を挟んで左側の，密度，圧力，x 軸方向速度，音速の絶対値，音響インピーダンスをそれぞれ ρL ,pL ,uL ,aL ,IL ，検査体積界面を挟んで右側の，密度，圧力，x 軸方向速度，音速の絶対値，音響インピーダンスをそれぞれ ρR ,pR ,uR ,aR ,IR とすると，検査体積界面の圧力 p∗ は，検査体積に対する法線方向速度 u∗n を用いて，式 (14),(15) のようになる．

(7)となる．
以上の条件のもと式 (1) を rin から rout まで積分すると，
R を円管の内半径，d を円管の肉厚であるとして，運動方程式は式 (8) のようになる．なお，計算の過程で，薄肉円管の仮定を用いた．

d2 s(t)E′1=−s(t) +(pin − pout )dt2
ρw R 2ρw d(8)

ここで，k1 = E ′ /(ρw R2 )，k2 = 1/(ρw d) とおくと，式 (8) は式 (9) と書くことができる．

d2 s(t)= −k1 s(t) + k2 (pin − pout )dt2= −k1 sn+1 + k2 (p∗n+1 − p0 )∆tn+1sn+1 = sn + {(1 − β)un}∆tw + βuw

(9)

(10)
(11)

2.2.1 支配方程式
ラグランジュ座標系における二次元非粘性流体に対する支配方程式は有限体積法により，式 (12),(13) のようになる．
∑Dpi=−(u∗ · n + uw )Sij
Dtmi
1 ∑ ∗Dui=−p nSij
Dtmi

(16)
(17)

なお，流れの最大マッハ数 M に対して，M ≥ 1 のときは，
β1 = β2 = 1，M < 1 のときは，β1 = 1/β2 ,β2 = M (2 − M )である (6) ．
壁面における圧力は，u∗n の代わりに，x 軸方向の壁面の速度 uw を用いて，式 (14) または，式 (15) より求めた．なお，図中の破線が壁面を表す．
以上，式 (10)〜(17) を連立させ，ラグランジュ座標におけるセル図心の圧力と速度，壁面の法線方向速度を陰解法を用いて求めた．
なお，陰解法は ILU+GMRES 反復法で解き，気液二相の解析には二流体モデルを用いた．

計算手法

Ii2(15)

IL uL + IR uR + (pL − pR )β1IL + IR
IL pR + IR pL + IL IR (uL − uR )β2∗p =IL + IRu∗n =なお，β は時間精度を決定する定数であり，本研究では β = 1である後退オイラー法を用いた．

2.2(14)p∗ = pR + IR (u∗n − uR )壁面以外の検査体積界面の値を求める際には，全速度化した音響ソルバーを用い，MUSCL 法によって空間二次精度とした (5) ．
このとき検査体積界面の圧力と法線方向速度は，式 (16),(17)のようになる．

2.1.2 基礎方程式の離散化
本研究では，Fig. 2 のような，円管の中心軸方向の影響を考慮しない，一次元弾性変形モデルを仮定し，式 (9) を離散化した．円管の外圧を p0 ，壁面の圧力を p∗ ，壁面法線方向速度をuw とすると，式 (8) は式 (10),(11) のように離散化できる．
un+1− unwwp∗ = pL − IL (u∗n − uL )2.2.3

ラグランジュ座標系からオイラー座標系への保存量のマッピング(12)

ラグランジュ座標系で解いた圧力と速度から，保存量を求め，オイラー座標系にマッピングする．Fig. 4 にその概念図を空間nn一次元にして示した．Qni−2 ,Qi−1 ,Qi は，ラグランジュ座標系(13)2c 2011 by JSFM
Copyright ⃝第 25 回数値流体力学シンポジウム
A06-1で計算する前の保存量を表す．また，Q′i−2 ,Q′i−1 ,Q′i は，ラグラn+1n+1ンジュ座標系で計算した後の保存量，Qn+1はオイi−2 ,Qi−1 ,Qi
ラー座標にマッピングした後の保存量を表す．更に， Ωi−2 ,Ωi−1n+1は内部セルの検査体積，Ωnは，n ステップ目と n + 1 ス
i ,Ωi
テップ目の壁面に接するセルの検査体積を表す．
cw
Elasti

実線の斜線で囲まれた部分は検査体積 i から壁面変位によって移動した保存量である．破線の斜線で囲まれた部分は，内部セルの界面において，x 軸の小さい方の検査体積から大きい方の検査体積に，流出した保存量を表している．各検査体積におけるマッピング後の保存量は，有限体積法により求めることができる．なお，流出した保存量は，壁面に対しては検査体積重心の物理量を，壁面以外の界面に対しては検査体積重心の物理量から勾配によって決定した．

allFig. 5 Volume change due to elastic wall displacement

Water
1.013×107 Pa
293 KΩn+1 = Ωn + ∆Ω

1.02 mm

また，式 (18) のように，n + 1 ステップ目の壁面に接するセルの検査体積 Ωn+1 は，Ωn と壁面変位 sn+1 によって増減した体積 ∆Ω から求めることができ，仮想の検査体積としてn + 1 ステップ目の計算に用いることとした（Fig. 5）．この手法により，格子を移動させることなく移動境界問題を解くことが可能となる．

(18)
52.7 mm

一方，本研究で開発した弾性変形モデルでは，∆Ω が負の値を持ち，かつ Ωn < |∆Ω| のような，格子サイズより大きい変形は取り扱うことができないことが分かる．

Fig. 6 Initial grid and condition

Initial states in the Eulerian frame

3.
Elastic wall

Grid interface

Grid interface

Grid interface

3.1

Grid interface

Elastic wall

3.2
Particle interface

Particle interface

Particle interface

Elastic wall

Grid interface

Grid interface

時間・空間収束性構築した弾性変形モデルの妥当性を検証するために，壁面変位の時間収束性，空間収束性について調べた．
まず，時間収束性を確かめるために，LEV4 の格子を用いて解析を行った．クーラン数を 0.40,0.20,0.10,0.05,0.025,0.0125と変化させたときの壁面変位の時間収束性を観察した．Fig. 7には，クーラン数が 0.40,0.10,0.025 のときの壁面外向き変位の時間履歴を結果を示す．
次に，空間収束性を確かめるために，クーラン数 0.025 のもと解析を行った．格子を LEV0 から LEV4 と変化させたときの壁面変位の空間収束性を観察した．Fig. 8 には，LEV0,LEV2,LEV4のときの壁面外向き変位の時間履歴を結果を示す．
なお，LEV4 の格子を用い，クーラン数が 0.00625 の値で解析した結果を収束解として両者の図に実線で示した．また，値はチャネル中央上壁部の値を用いた．

Final states in the Eulerian frame after remapping

Grid interface

計算条件構築したモデルを検証するために，Fig. 6 のような短軸が1.02 mm，長軸が 52.7 mm の単純チャネルを用いた．初期条
件は，チャネル内部が 1.013 × 107 Pa，293 K の水でチャネル内が満たさた高圧状態とする．境界条件はチャネルの短辺を剛体壁の境界条件とし，チャネル長辺に弾性モデルを適用させた．
初期格子は Fig. 6 の通りである．この格子の各辺を 2 分割したものから，16 分割した格子まで計 5 種類の格子を使用した．格子が粗いものから順に，LEV0,LEV1,…LEV4 とする．
また，式 (10) における定数 k1 ，k2 の値には，ステンレス鋼の値を用いた．なお，ステンレス鋼のヤング率，ポアソン比，密度はそれぞれ，200 GPa，3.0，7800 kg/m3 ，管の肉厚 d は0.20 mm とした．

Intermediate states in the Lagrangian frame

Particle interface

弾性変形モデルの検証

本章では，構築した弾性変形モデルについて行った検証について説明する．

Grid interface

Fig. 4 Mapping from Lagrangian coordinate to Euleriancoordinate

3c 2011 by JSFM
Copyright ⃝第 25 回数値流体力学シンポジウム
A06-1

どちらの場合も，時間が経つにつれて，壁面の変位が一定の値に収束していることが分かる．また，クーラン数が小さくなるか，格子が細かくなると，壁面の振動周期は小さくなり，やがて，一定の値になることが分かる．以上のことから，時間と空間の収束性を確認することができた．
また，クーラン数が大きい場合や，格子が粗い場合は，振動が収束する時間が短く，壁面変位の振幅が小さくなることも確かめられた．この理由は，時間刻みが大きい場合，高周波の波を捉えられないからであると考察できる．
それぞれの値において，クーラン数，格子幅，水の音速から，時間刻みを算出し，各時間刻みにおける壁面の振動周期を表したものが，Fig. 9 である．なお，水の音速は 1500 m/s とした．
Fig. 9 から，時間と空間の収束性が確認できる．また，振動周期は音波がチャネルを短軸方向に往復する周期のおおよそ半分になっていることも確認できる．このことから，チャネルの長軸方向の中心軸に対して対称的に，波が往復運動し，それに起因した周期で壁面が変位していることが分かる．
以上のことは Fig. 10 からも確認できる．Fig. 10 には壁面のチャネル外向きの変位，壁面の圧力，壁面に接するセルの図
心の圧力，チャネル中心の圧力を示した．なお，格子は LEV4を用い，クーラン数は 0.025 とした．壁面の圧力とチャネル中心の圧力は，逆位相になっていることが確認できる．また，壁面の変位は壁面圧力の変動と連成しながら振動していることも分かる．壁面変位の収束原因は，管壁が大気圧に対して仕事をするためであると考えられる．
1.00e-007

1.00e-006
9.50e-007
9.00e-007

Period of vibration [s]8.50e-007

7.50e-007
7.00e-007

6.00e-007
5.50e-007
5.00e-007
0.00e+000

1.00e-009

2.00e-009

3.00e-009∆t [s]

4.00e-009

5.00e-009

6.00e-009

Fig. 9 Convergence of the period of vibration
1.00e-007

1.15e+007
Outward wall displacement
Wall pressure
Pressure at the centroid of the cell on wall
Pressure at the center of channel

9.00e-008

1.10e+007

7.00e-008
1.05e+007
6.00e-008
1.00e+007

5.00e-008

Pressure [Pa]Outward wall displacement [m]8.00e-008

4.00e-008
9.50e+006
3.00e-008
2.00e-008
9.00e+006
1.00e-008
0.00e+000
0.00e+000

1.00e-006

2.00e-006

3.00e-006

4.00e-006

8.50e+006
5.00e-006

Time [s]Fig. 10 Time history of displacement and pressure

4.

8.00e-008
Outward displacement [m]8.00e-007

6.50e-007

Exact=0.40=0.10=0.025

9.00e-008

Time convergence
Grid convergence

結言

本稿ではまず，本研究で開発した，管壁の弾性変形を考慮した数値解析モデルを提案した．更に，高圧単純チャネルの管壁の変位と圧力を調べ，検証を行った．その結果，時間と空間収束性により弾性変形モデルの妥当性が確かめられた．

7.00e-008
6.00e-008
5.00e-008
4.00e-008
3.00e-008
2.00e-008

参考文献

1.00e-008
0.00e+000
0.00e+000

2.00e-006

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

8.00e-006

Fig. 7 Time convergence of wall displacement

1.00e-007

Exact
LEV 0LEV 2LEV 4

9.00e-008(2) Muhd Hilmi Bin Shapien, Mingyu Sun, “Numerical
Simulation of Laser-Induced Liquid Jet with Phase
Change”, Symposium on shock wave in Japan (2011).

Outward displacement [m]8.00e-008
7.00e-008
6.00e-008(3) 高橋幸伯, 町田進, 角洋一, 基礎材料力学, 三訂第 4 刷(2009), pp.154-160, (株) 培風館.

5.00e-008
4.00e-008
3.00e-008(4) Eleuterio F. Toro, Riemann Solvers and Numerical
Methods for Fluids Dynamics, Third Edition (1999),
pp.293-304, Springer-Verlag.

2.00e-008
1.00e-008
0.00e+000
0.00e+000

2.00e-006

4.00e-006
Time [s]6.00e-006

8.00e-006(5) 矢田和之, 市東素明, 孫明宇, “非構造格子を用いた allspeed 流れ解析”, 第 24 回数値流体力学シンポジウム,
2010.

Fig. 8 Grid convergence of wall displacement(6) Meng-Sing Liou, “A sequel to AUSM, Part II: AUSM+-up for all speeds”, Journal of Computational Physics,
Vol. 214 (2006), pp.137-170, 2006.
4c 2011 by JSFM
Copyright ⃝