潜水艇タイタン:深海における気泡の爆縮シミュレーション
By M.Sato
概要
先日、大西洋沖において潜水艇タイタンが深海において遭難する事故が発生したことは記憶に新しい。
事故原因についての詳細は不明であるが、深海での水圧に耐えられずに潜水艇は圧潰したのではないかと言われている。
そこで、今回は流体力学的観点から、深海において気泡の爆縮が発生するとどのようになるのかをシミュレーションを行った。
なお、計算に際してはいくつかの仮定、モデル化を行っているが、これらの妥当性については今後、機会を見て検討して行きたいと思う。
想定した流れは、深海3000mにおける半径1mの気泡の爆縮に伴う流れである。
本記事で使用した解析データについては、下記よりお問い合わせください。
👉 お問い合わせ
気泡の挙動を表す方程式
無限領域内に置かれた球(気泡)周囲の液体の流れを考える（対象とするのは、球内部ではなく球周囲の液体部分である）。
液体は非圧縮性で層流とする。気泡は原点に置かれているとする。気泡周囲の流れは半径方向のみとする。
気泡モデル
運動方程式
気泡周囲の液体部について運動量保存式を球座標系で記述すると、 速度成分は$r$方向のみであるから、以下のように簡略化される。
$$ \rho_l \left(\frac{\partial u}{\partial t} + u\frac{\partial u}{\partial r} \right)= - \frac{\partial p}{\partial r} + \mu_l \left[ \frac{1}{r^2}\frac{\partial}{\partial r} \left( r^2 \frac{\partial u}{\partial r} \right) -\frac{2u}{r^2} \right] \tag{1} $$
ここで、$t$は時間、$r$は気泡中心からの半径方向距離、$u$は半径方向速度、$\rho_l$は液体密度、$\mu_l$は液体粘度、$p$は圧力である。
連続方程式
非圧縮性流体では質量が保存することから、以下の関係式が成り立つ。
$$ u r^2 = \dot{R} R^2 = \frac{d R}{dt}R^2 \tag{2} $$
$R(t)$は気泡の半径である。変形すると
$$ u(r,t) = \frac{R^2 \dot{R}}{r^2} \tag{3} $$
非回転流れであることから速度ポテンシャルが存在する。 半径方向速度は速度ポテンシャル$\Phi$を使い、$u=\frac{\partial \Phi}{\partial r}$と表されることから
$$ \Phi = \frac{-R^2 \dot{R}}{r^2} \tag{4} $$
(1)式右辺の粘性項は(3)式を使うと相殺されてゼロになる。従って最終的に(1)式は以下のようになる。
$$ \rho_l \left[ \frac{\partial}{\partial r} \left(\frac{\partial\Phi}{\partial t}\right)+ \frac{\partial}{\partial r} \left(\frac{u^2}{2}\right) \right] +\frac{\partial p}{\partial r} =0 \tag{5} $$
この式を$r$について積分すると、非定常ベルヌイ式が得られる。$F(t)$は積分定数である。
$$ \rho_l \left[ \frac{\partial\Phi}{\partial t}+ \frac{u^2}{2} \right] +p =F(t) \tag{6} $$
$r\rightarrow \infty$の場合、$\Phi=0,p=p_\infty$であることから $F(t)=p_\infty$となる。 従って
$$ p=p_\infty - \rho_l \left( \frac{\partial\Phi}{\partial t}+\frac{u^2}{2} \right) \tag{7} $$
(4)式を代入すると
$$ p(r,t)=p_\infty - \rho_l \left[ -\frac{R^2\ddot{R}}{r}-\frac{2R\dot{R}^2}{r}+\frac{\dot{R}^2R^4}{2r^4} \right] \tag{8} $$
(8)式において$r=R(t)$とすると、以下のRayleigh 方程式が得られる。
$$ R\ddot{R}+\frac{3}{2}\dot{R}^2=\frac{p_a-p_\infty}{\rho_l} \tag{9} $$
ここで $p_a$は、液体-気泡界面における液体側圧力である。
Rayleigh-Plesset方程式
気泡表面での液体側応力は次式で表される。
$$ \sigma_{rr} =-p_a + 2\mu\frac{\partial u_r}{\partial r}-\frac{2}{3}\mu\left( \frac{\partial u_r}{\partial r}+\frac{2u_r}{r} \right) $$ $$ =-p_a+\frac{4\mu}{3} \left( \frac{\partial u_r}{\partial r} - \frac{u}{r} \right) \tag{10} $$
ここで$\partial u_r/\partial r$は次式で表される。
$$ \frac{\partial u_r}{\partial r}=\frac{\partial}{\partial r} \left(\frac{R^2\dot{R}}{r^2} \right)= -2\frac{R^2\dot{R}}{r^3} \tag{11} $$
$r=R$とすると、以下の応力表現式が得られる。
$$ \sigma_{rr}(r=R)=-p_a - 4\mu\frac{\dot{R}}{R} \tag{12} $$
気泡の表面においてラプラスの式が成り立つと仮定すると
$$ p_B+\sigma_{rr}=\frac{2S}{R} \tag{13} $$
ここで $S$は気液界面での表面張力である。
液体-気泡界面における気体側圧力を $p_B$とすると、界面における力の釣り合いから次式が成り立つ。
$$ p_a = p_B - 4\mu_l \frac{\dot{R}}{R} - \frac{2S}{R} \tag{14} $$
気泡圧力$p_B$を使ってRayleigh方程式を整理すると、以下のRayleigh-plesset方程式が得られる。
$$ \frac{p_B-p_\infty}{\rho_l} = R\ddot{R}+\frac{3}{2}\dot{R}^2+\frac{4\mu_l}{\rho_l}\frac{\dot{R}}{R}+\frac{2S}{\rho_l R} \tag{15} $$
気泡が極めて短い時間で爆縮する場合、気泡内外での熱の移動は起こらないとみなしてもよい。 気泡内部では均一の圧力になっていると仮定すると、気泡半径$R$と気泡圧力$p_B$には以下の等エントロピー関係式が成り立つ。
$$ \frac{p_B}{p_{B0}}=\left(\frac{V_0}{V}\right)^\gamma=\left(\frac{R_0}{R}\right)^{3\gamma} \tag{16} $$
ここで、$V, R$はそれぞれ、気泡体積、気泡半径を表しており、添字ゼロは初期値を表している。$\gamma$は比熱比である。
これらの式を組み合わせると最終的に$R$に関する二階常微分方程式が得られる。
$$ R\ddot{R}+\frac{3}{2}\dot{R}^2+\frac{4\mu_l}{\rho_l}\frac{\dot{R}}{R}+\frac{2S}{\rho_l R}-\frac{1}{\rho_l}\left[ p_{B0}\left( \frac{R_0}{R} \right)^{3\gamma} - p_\infty \right]=0 \tag{17} $$
液相内の圧力の空間分布は次式で表される。
$$ p(r,t)=p_\infty\left[ 1-\frac{R}{r} \right] +p_B\frac{R}{r} -\frac{2}{r}(2\mu\dot{R}+\sigma) +\rho_l\frac{R\dot{R}^2}{2r}\left[1-\frac{R^3}{r^3} \right] \tag{18} $$
juliaを使った気泡挙動の計算
連立一階常微分方程式
Rayleigh-plesset方程式は、二階常微分方程式となっているため連立一階常微分方程式に変換して計算を行った。 以下に連立微分方程式を示す。
$$ \frac{d R}{d t} =U \tag{19} $$
$$ \frac{d U}{d t} = \frac{1}{R}\left[ \frac{1}{\rho_l}\left[ p_{B0}\left( \frac{R_0}{R} \right)^{3\gamma} - p_\infty \right] -\frac{3}{2}U^2-\frac{4\mu_l}{\rho_l}\frac{U}{R}-\frac{2S}{\rho_l R} \right] \tag{20} $$
ここで、
- $R$:気泡半径
- $U$:気泡界面の速度（外向きを正）
- $\rho_l$:液体密度(=1000[kg/m^3])
- $\mu_l$:液体粘度(=1e-3[Pa.s])
- $S$:表面張力係数(=0.07[N/m])
- $\gamma$:比熱比（=1.4)
- $p_{B0}$:初期気泡圧力（＝大気圧＝1e5[Pa]）
- $p_\infty$:気泡周囲の液相圧力（=3e7[Pa]）
- $R_0$:初期気泡半径（=1[m])
解析プログラム
常微分方程式を計算するために juliaを使用した。
式中に含まれるパラメータは、深度3000mにおける、半径1[m]の球形気泡を想定した。
なお気泡内部は人間が滞在できる環境を想定し20℃、大気圧の空気とした。具体的な数値は以下のプログラムリストに記述している。
解析は時刻ゼロにおいて水深3000[m]に配置した球形気泡が爆縮する様子を計算する。
using ParameterizedFunctions, DifferentialEquations
using Plots
using Printf
# 常微分方程式の定義
eqns = @ode_def bubble begin
dR = V
dV = ((patm*(Rini/R)^(3.0*gamma)-pinf)/rhoL-1.5*V^2-4.0*nuL*V/R-2.0*sigma/(rhoL*R))/R
end patm pinf Rini rhoL nuL gamma sigma
patm = 1e5 # [Pa]
pinf = 3e7 # [Pa]
rho0 = 1.2 # [kg/m3]
Tini = 293.0 # [K]
Rini = 1.0e-0 # [m]
rhoL = 1000.0 # [kg/m3]
nuL = 1e-6 # [m2/s]
gamma = 1.4 # [-]
sigma = 0.07 # [N/m]
# 問題定義
u0 = [Rini; 0.0] # 初期 R, V
tspan = (0.0, 1e-2) # 時間積分期間
p = [ patm, pinf, Rini, rhoL, nuL, gamma, sigma ]
prob = ODEProblem(eqns,u0,tspan,p)
# 計算
@time sol = solve(prob,Tsit5(),reltol=1e-8,abstol=1e-8)
(nrow,ncol)=size(sol)
times=zeros(Float64,0)
radis=zeros(Float64,0)
velos=zeros(Float64,0)
press=zeros(Float64,0)
temps=zeros(Float64,0)
@printf("j : time radi velo pres rho temp\n")
for j=1:ncol
time = sol.t[j]
radi = sol.u[j][1]
velo = sol.u[j][2]
pres = patm*(Rini/radi)^(3.0*gamma)
rho = rho0*(Rini/radi)^3.0
temp = Tini*(Rini/radi)^(3.0*(gamma-1.0))
append!(times,time)
append!(radis,radi)
append!(velos,velo)
append!(press,pres)
append!(temps,temp)
@printf("j=%d : %16.8e %16.8e %16.8e %16.8e %16.8e %16.8e\n",
j,time,radi,velo,pres,rho,temp)
end
p1 = plot(times, radis, title="history of radius[m]",
legend=:bottomleft, linewidth=1, linestyle=:solid,
xlabel="time[sec]", ylabel="r[m]", label="radius" )
savefig("radi.png")
p2 = plot(times, velos, title="history of velocity[m/s]",
legend=:bottomright, linewidth=1, linestyle=:solid,
xlabel="time[sec]", ylabel="v[m/s]", label="velocity" )
savefig("velo.png")
p3 = plot(times, press, title="history of pressure[Pa]",
legend=:topleft, linewidth=1, linestyle=:solid,
xlabel="time[sec]", ylabel="p[Pa]", label="pressure" )
savefig("pres.png")
p4 = plot(times, temps, title="history of temperature[K]",
legend=:topright, linewidth=1, linestyle=:solid,
xlabel="time[sec]", ylabel="T[K]", label="temperature" )
savefig("temp.png")
plt = plot(p1, p2, p3, p4; layout=(4,1), size=(640, 960) )
savefig("all.png")
解析結果
計算結果によれば気泡半径は周期的に増減を繰り返し、またそれに付随して、気泡界面速度、圧力、温度も周期的に増減を繰り返していることがわかった。
今回の計算では気泡が崩壊、移動する影響は考慮していないため、長期間に渡ってこのような振動する解になっているが、実際の現象は最初の振動ピーク時刻の手前において気泡は破壊するものと考えられる。
気泡の半径$R$
以下に時刻 0から0.01[秒]までの結果図を示す。
気泡半径の経時変化
気泡界面の速度
気泡半径が増減することから、気泡界面速度は周期的に符号が反転している。ピークとなる流速絶対値は$O(10^4)$となっており、実際にはこのようなピークは発生しないと考えられる。
気泡界面速度の経時変化
気泡の圧力
気泡の圧力は気泡半径が最小となる時刻において最大となっている。 ただし圧力のピーク値は$O(10^{12})$[Pa]となっており非現実的な値となっている。
気泡圧力の経時変化
気泡の温度
時刻$t=5\times10^{-3}[sec]$における圧力の空間分布を示す。
この図は気泡半径$R(=0.402[m])$から$r=10[m]$の範囲の圧力分布を示している。 気泡表面付近において圧力が急激に上昇していることがわかる。
圧力の半径方向分布
fluentを使ったシミュレーション
Rayleigh-plesset方程式では、液相を非圧縮性流体と仮定している。
実際には、液相内の流速が音速に比較して相対的に無視できない程度に大きい場合には、圧縮性流体（密度変化考慮）として現象を調べることが必要である。
ここでは、fluentの圧縮性液体の機能を使い解析を行ってみた。
ただし気泡界面での速度は前記juliaの計算結果を利用した。
解析モデル
解析領域
軸対称モデルを作成した。 水面下3000[m]の位置に気泡があり、海底は水深4000[m]であると仮定した。気泡の半径は1[m]とした。
モデル形状
検査点位置
メッシュ分割
全体領域のメッシュ分割を示す。
メッシュ分割図
物性値
fluentではテイト(Tait)の状態方程式を使って圧縮性液体をモデル化している。
これは、等温条件下における密度と圧力の非線形関係を記述するものである。 テイトの式を単純化すると、次のように表すことができる。
$$ \left( \frac{\rho}{\rho_0} \right)^n=\frac{K}{K_0} $$
ここで
$$ K=K_0 + n\Delta p $$
$$ \Delta p=p-p_0 $$
- $p_0$：参照液体圧力（絶対値）
- $\rho_0$：参照液体密度（$p_0$における密度）
- $K_0$：参照体積弾性係数（$p_0$における体積弾性係数）
- $n$：密度指数
- $p$：液体圧力（絶対値）
- $\rho$：圧力$p$における液体密度
- $K$：圧力$p$における体積弾性係数
音速 $c$は次式で表される。
$$ c=\sqrt{\frac{K}{\rho}} $$
ここでは、以下の値を設定した。
液体密度
境界条件
今回の計算では気泡内部のガスの計算は行わない。 また、気泡境界面の移動はないと仮定した（気泡半径1[m]）。
作動条件
流体密度の計算で使用する圧力は絶対圧である。
従って静水圧を考慮するように基準密度をゼロと設定した。
作動条件
気泡界面速度の時間パターン
Rayleigh-plesset方程式の計算で得られた気泡界面の移動速度から境界値を推定した。
気泡挙動の計算結果によれば、時刻$t=5\times10^{-3}[sec]$において、気泡半径 $R=0.4[m]$、気泡界面速度$\dot{R}=500[m/s]$となっている。
従って、$t=5\times10^{-3}[sec]$において最大値となるような三角形パルス状の半径方向流入速度を仮定した。
パルス入力
流出入境界
気泡表面では速度を時間パターンデータとして指定した（ttabファイル使用）。
側面境界では、圧力が静水圧分布となっていると仮定した。
また、気泡から伝わると想定される圧力波が到達した場合には、圧力波がそのまま透過できるように無反射境界とした。
流出入境界
水面
水面はスリップ境界と近似した。
パラメータ設定
計算で使用するパラメータを以下のように設定した。
static_presは、側面境界における静水圧を設定するために使用している。
diff_presは静水圧からの変動圧力量を求めている(ポスト処理で使用)。
パラメータ
解析結果
解析は非定常で実施することになるが、予め気泡界面速度をゼロと設定して初期圧力分布を求めた。
以下に全体場における圧力分布を示す。 気泡位置は水深3000[m]であり、$p=3\times10^7[Pa]$程度となっている。
初期圧力
その後、時間刻みを $\Delta t=10^{-5}[sec]$ として非定常解析を実施した。
以下に解析結果を示す。
流速ベクトル
時刻 $t=10^{-2}[sec]$における流速ベクトル図を示す。 矢印長さが同一になるように表示している。
流速ベクトル
流速強度コンター
気泡中心から100[m]程度の範囲を対象として流速強度のコンターを示す。 気泡中心から同心円状に変動が伝わっていることが分かる。
流速強度コンター
圧力コンター
気泡中心から100[m]程度の範囲を対象として圧力のコンターを示す。 圧力レンジを調整して圧力変動が分かるようにしている。
圧力コンター
圧力差分コンター
気泡中心から100[m]程度の範囲を対象として圧力差分(diff_pres)のコンターを示す。 これは計算された圧力値と静水圧の差分を表している。
圧力差分コンター
音速コンター
気泡中心から100[m]程度の範囲を対象として音速のコンターを示す。 水中の音速は大体 1500[m/s]程度であり、水深が深い（圧力が大きい）ほど音速が大きくなっている。
音速コンター
検査点における圧力変動
気泡中心から10[m]、100[m]の位置に検査点を配置し、そこでの圧力の経時変化を調べてみた。
- probe-Aでは、時刻t=0.006～0.02[sec]付近で大きい圧力変動が発生している
- probe-Bでは、時刻t=0.07～0.09[sec]付近で微小な圧力変動が見られ、減衰した圧力波が伝播している
ことがわかる。
圧力経時変化
結論
潜水艇タイタンの事故で発生したと思われる爆縮現象の解析を行ってみた。
不自然な仮定も含まれていると思いますが今後の調査、研究の参考になれば幸いです。
弊社の解析事例
弊社の流体解析事例については、下記のリンクからご覧ください．