読者です 読者をやめる 読者になる 読者になる

とつとつとしてろうとせず

ひまつぶしにどうぞ。

液膜のRayliegh-Taylor不安定性に対する線形安定性解析

注意

以下の記述は私的な備忘録であり、内容の正確さに関しては一切保証いたしませんので、悪しからず。

問題設定


図のように壁面に塗布された薄い液膜の2次元流動を考える(たとえば、風呂場の天井に付着した水の膜など)。壁面表面から鉛直下向きに$z$軸を取り、壁面に沿って$x$軸を取る。最初に厚さ$h_{0}$の均一な液膜があり、そこに重力$g$によって正弦波状の表面変形$\zeta(x,t)$が生じたとすると、液膜厚さ$z=h(x,t)$は$h=h_{0}+\zeta(x,t)$と表される。液体の密度を$\rho$、粘度を$\mu$とし、液膜表面には表面張力$\gamma$が働くものとする。液体の下にある空気の運動は考えず、気液界面には応力が働かないとする(自由表面)。また、$x$方向は周期境界条件とし、側壁の影響は無視する。

基礎式の導出

まず、解析の基礎式となる液膜方程式を導く。非圧縮性流体の2次元運動は次のNavier-Stokes方程式によって表される。
\begin{align}
& \rho \paren{ \pdif{u}{t} + u \pdif{u}{x} + w\pdif{u}{z} } = - \pdif{p}{x} + \mu \paren{ \pddif{u}{x} + \pddif{u}{z} } \\
& \rho \paren{ \pdif{w}{t} + u \pdif{w}{x} + w\pdif{w}{z} } = - \pdif{p}{z} + \mu \paren{ \pddif{w}{x} + \pddif{w}{z} } + \rho g
\end{align}また、質量保存則は次式の連続の式で表される。
\begin{align}
& \pdif{u}{x} + \pdif{w}{z} = 0
\end{align}ここで、$(u,w)$はそれぞれ$(x,z)$方向の流速、$p$は圧力である。
境界条件としては壁面上ですべりなし条件$u|_{z=0} = w|_{z=0} =0$を課す。また、自由表面上($z=h$)では速度の境界条件として接線応力なし、圧力の境界条件として次のYoung-Laplaceの式を用いる。
\begin{align}
& \mu \pdif{u}{z} = 0 \\
& p = - \gamma \kappa
\end{align}ここで、$p$は大気圧を基準とするゲージ圧、$\gamma$は表面張力、$\kappa$は表面形状の曲率であり、液膜厚さが基板からの高さ(高さ関数)で与えられる場合には次式で表される。
\begin{align}
\kappa = \frac{\partial_{x}^{2} h}{\sqrt{ 1+ \paren{\partial_{x} h}^{2} }}
\end{align}ここで、$\partial_{x} = \partial/\partial x$である。
解析を簡単にするため、基礎式に対して潤滑近似を適用する。すなわち、表面変形によって生じる流速は小さく、慣性力は無視できるとする。また、液膜が非常に薄く、表面変形の周期$\lambda$を代表長さにとると、$h_{0} \ll \lambda$であるとする(長波長近似)。この結果、連続の式のオーダーの比較により、$w \ll u$が分かるので、$z$方向の流れは無視できる。同様に$\partial_{x}^{2} u \ll \partial_{z}^{2} u$であることも分かる。これより、Navier-Stokes方程式は次式のようになる。
\begin{align}
\label{eq:u}
& 0 = - \pdif{p}{x} + \mu \pddif{u}{z} \\
\label{eq:w}
& 0 = - \pdif{p}{z} + \rho g
\end{align}
式\eqref{eq:w}を$z$で積分して境界条件を考慮すると、圧力の式を得る。
\begin{align}
p = \rho g \paren{z-h} - \gamma \kappa
\end{align}また、式\eqref{eq:u}を$z$で積分して境界条件を考慮すると、$u$は次式で与えられる。
\begin{align}
u = \frac{1}{2\mu} \pdif{p}{x} z \paren{z - 2h}
\end{align}したがって、
\begin{align}
u = \frac{1}{2\mu} z \paren{z - 2h} \pdif{}{x} \paren{ - \rho g h - \gamma \kappa }
\end{align}上式を連続の式に代入して$z$について$0$から$h$まで積分すると、$w|_{z=h} = \partial_{t} h$より、
\begin{align}
& \pdif{h}{t} = - \frac{1}{3\mu} \pdif{}{x} \paren{ h^{3} \pdif{}{x} \paren{ \rho g h + \gamma \kappa } }
\end{align}これが液膜の時間発展を表す方程式である。いまの場合、$h=h_{0}+\zeta(x,t)$としているので、
\begin{align}
& \pdif{\zeta}{t} = - \frac{1}{3\mu} \pdif{}{x} \paren{ \paren{h_{0} + \zeta}^{3} \pdif{}{x} \paren{ \rho g \zeta + \gamma \kappa } }
\end{align}
解析の見通しをよくするため、上式を無次元化する。代表長さを$L$、代表時間を$T$、代表速度を$U$とし、
\[ L \sim h_{0} \ ,\ T \sim \frac{3\mu h_{0}}{\gamma} \ ,\ U \sim \frac{L}{T} = \frac{\gamma}{3\mu} \]
のようにスケールを設定し、基礎式を無次元化すると次式のようになる。
\begin{align}
\label{eq:film}
& \pdif{\zeta}{t} = - \pdif{}{x} \paren{ \paren{1 + \zeta}^{3} \pdif{}{x} \paren{ B \zeta + \kappa } }
\end{align}ここで、$B = \rho g h_{0}^{2}/\gamma$はBond数(あるいはEötvös数)である。

線形安定性解析

得られた基礎式\eqref{eq:film}に対して線形安定性解析を行う。いま、液膜の表面変形$\zeta$が次式のようなFourierモードで表されるとする。
\begin{align}
\label{eq:fourier}
\zeta(x,t) = \varepsilon \exp\paren{\omega t + i k x}
\end{align}ここで、$\varepsilon$は微小な振幅、$\omega$は成長速度、$k$は波数である。上式を式\eqref{eq:film}に代入し、$\varepsilon$について1次の項までを残して線形化すると、$\kappa \approx \partial_{x}^{2}\zeta$より、$O\left(\varepsilon\right)$で次式を得る。
\begin{align}
& \pdif{\zeta}{t} = - \paren{ B \pddif{\zeta}{x} + \frac{\partial^{4}\zeta}{\partial x^{4}} }
\end{align}上式に式\eqref{eq:fourier}を代入することで次式の線形分散関係式を得る。
\begin{align}
& \omega = - k^{2}\paren{ k^{2} - B }
\end{align}上式は平方完成によって次のように変形できる。
\begin{align}
\label{eq:omega}
& \omega = - \paren{ k^{2} - \frac{B}{2} }^{2} + \frac{B^{2}}{4}
\end{align}これより、$\omega = \omega(k)$のグラフは$k > 0$の右半平面で上に凸の曲線であることがわかる。最も成長しやすいモードは最大成長速度の条件$\displaystyle \dif{\omega}{k} = 0$より、$\displaystyle k_{\max} = \sqrt{\frac{B}{2}}$のモードであり、そのときの最大値は$\displaystyle \omega_{\max} = \frac{B^{2}}{4}$である(このことは\eqref{eq:omega}の関数形からもすぐに分かる)。$\displaystyle k = \frac{2\pi}{\lambda}$より、無次元波長$\lambda$について表せば最小波長は$\displaystyle \lambda_{\min} = 2\pi \sqrt{\frac{2}{B}}$で与えられる。

図に$B$の値を変えたときの$\omega$のグラフを示す。$\omega > 0$のとき界面は不安定であるので、この系は$B>0$のとき無条件で不安定であることがわかる。$B$は重力(浮力)と表面張力の比であり、$B$を大きくしていったときに$\omega > 0$の領域が大きくなることから、重力は界面を不安定に、表面張力は界面を安定にする効果があることがわかる。

おわりに

線形安定性解析の練習としてこの問題が取り上げられることが多いが、自分でやってみたことがなかったので計算してみた。流体力学に限らず、色んな物理系に対して安定性解析をしたり、数値計算の方法をまとめた教科書とかないのかね。力学系の話とかは具体例から戻る形で勉強しないと個人的に分かりづらくてね……。