数理解析研究所講究録 1271 巻 2002 年 248-257248
2 層円環形液体シートの安定性解析と薄膜近似
京大・情報佐野雅之 (Masayuki Sano)

京大・情報船越満明 (Mitsuaki Funakoehi)

Department of Applied Analysis and Complex Dynamical Systems
Graduate School of Informatics, Kyoto Univercity

1Iri.troduction.

液体シート の安定性の問題は掖体スプレーの形成や、液体シートによる遮蔽などの観点から重要であり埋論的には今までにさまざまな場合が調べられてきた [1]。液体が薄いシート状になって流このような液体シートの形としては平面 2 次元状シート、 中空の円環形シートなどが考えられる。

1 層シートにおいては長波長で 2 つの不安定モードがある事が知られている [2][3]。 一つは界面が同じ方向に変位する対称 (sinuous) モードであり、 もう一つは界面が逆方向に変位する反対称(varicoee) モードである。
線形解析の結果として反対称モードの不安定成長率は対称モードの不安定成長率よりも大きい事が知られている。 1 層平面 2 次元シートにおいては不安定は界面における Kelvin-Helmholtz 不安

定であり、真空中では不安定は生じないが、 円環形シートにおける不安定は界面における表面張力のみで生じる不安定であって真空中でも発生する [6][7]。

1 層液体シートの崩壊過程の詳細を知るためには非線形の効果を取り入れた時間発展を調べる必要がある。 しかし液体シートは移動する自由境界を含む系であるのでポテンシャル理論を用いた直接計算は計算量が多いという面で困難である。 したがって何らかの近似を用いて系をモデル化し、より計算しやすいようなモデル方程式を導く事が有効である。そのようなモデルとして境界面を渦層として表す離散渦層モデル [5] とシートの厚さ方向に平均化したナビエ・ストークスを用いる薄膜近似モデルがある [61[7]。 これらのモデル化によって 1 層液体シートの崩壊過程が計算されており、 1 層円環形シートの場合には中空部が閉じて崩壊する場合やシートの厚さが 0 になって崩壊する場合があることが示されている。
ここでは 2 層の中空円環形液体シートを考える。 2, 層の液体シートでは対称モード、 反対称モードの他にもう一つのモードが現れる。以下では 2 層円環形液体シートの線形解析を行ない待間的不安定の成長率の振舞いを調べる。 また非線形の効果を取り入れた解析を行なう事を日標として 2層円環形液体シートに薄膜近似を適用しモデル方程式を導き、その妥当性をポテンシャル理論との比較によって検証する事とする。

2

ポテンシャル理論による 2 層円環形液体シートの安定性解析

21

基礎方程式と分散関係の導出

ここでは 2 層円環形液体シートの時間的安定性をポテンシャル理論によって調べる。 図 1 のよ軸とする円筒座標系で考える。以下では液体シートおよひ周囲の流体を構成する流体はすべて非粘性の非圧縮流体と仮定する。 内側の円環シートを構成
外側の円環シートを構成する液体を流体 2 とし、 図のように 2 層に重なった円する液体を流体
うに円環の軸を$z$

軸、 円環の半径方向を$r$$1_{\text{、}}$

環形シートを円環の内側の領域、 内側の液体シート、外側の液体シート、 円環の外側の領域の領域

249

図 1: 2 層円環形液体シート[こ分け、 それぞれの領域での速度ポテンシャ)、 圧力、 密度を $\phi_{:}(r, z, t),p_{i}(r, z, t),$ $\rho_{i}(r, z, t)$ とおく。 ここで $i=0\ldots 3$ がそれぞれ円環の内側の領域、 内側の液体シート、外側の液体シート、 円環の外側の領域を表す。 流体 1 と円環の内側の領域との界面が $r=R_{1}$ ( ,t)(内側界面)、 流体 1 と流体$z$

2 との界面が $r=R_{2}$ ( , t)(中心界面)、 流体 2 と円環の外側の領域との界面が $r=R_{3}$ ( , t)(外側界面) 与えられるとし、それぞれの界面での表面張力係数がでであるとする。流体$z$$z$$\sigma_{01},$ $\sigma_{12},$の液体の $r,$$z$方向の速度をそれぞれ$v_{r,z}^{(i)}$$1_{\text{、}}2$$\sigma_{23}$とする。安定性を調べるべき定常流はシートの内部では$z$

軸方向に一定速度で一様な流れであり、 $v_{z}^{(1)}=V_{1},v_{z}^{(2)}=V_{2}$ とするとこの定常流に対する速度ポテ
ンシャルは$\phi_{i}=V_{\dot{l}}z$$(i=1,2)$ ,で与えられる。 また定常流では境界面は $R_{i}=R_{i0}$(1)$(i=0,3)$$\phi_{i}=0$$(i=1,2,3)$ であるとする。ただし R 一ま定数である。 また定常流での各領域内での圧力は一定であり、 表面張力による境界面での法線応力の差を考慮して$p_{3}=0$,$p_{2}= \frac{\sigma_{23}}{R_{30}}$,$p_{1}= \frac{\sigma_{23}}{R_{30}}\neq\frac{\sigma_{1_{\mathrm{A}}^{\circ}}}{R_{20}}$$(_{\backslash }2)$$p_{0}= \frac{\sigma_{23}}{R_{30}}+\frac{\sigma_{12}}{R_{20}}+\frac{\sigma_{01}}{R_{10}}$とする。 この定常流の微小撹乱にたいする線形安定性を調べるために、軸対称で$z$

軸方向への正弦

波的な微小擾乱として次の形を仮定する。

(3)$R_{i}(z, t)=R_{i0}+\eta_{i}\exp(i\omega t-ikz))$

この微小擾乱に対応する速度ポテンシャルの解を求める。速度ポテンシャルは各領域でラブラス方
程式$\Delta\phi_{i}=0$を満たしている。 円筒内部では$rarrow \mathrm{O}$(4)$(i=0\ldots 3)$

で有限、 円筒外部では r\rightarrow \propto 》で有限であるという条件のもとで微小擾乱に対する速度ポテンシャルの解は次のように書ける。

ここで$K_{0},$$I_{0}$$\phi_{3}(r, z, t)$$=$$A_{3}I\mathrm{f}\mathrm{o}(kr)\mathrm{e}\mathrm{x}_{\mathrm{t}}\mathrm{p}i(kz-\omega t)$(5)$\phi_{2}(r, z,t)$$=$$V_{2}z+(A_{2}K_{0}(kr)+B_{2}I_{0}(kr))\exp i(kz-\omega t)$(6)$\phi_{1}(r, z,t)$$=$$V_{1}^{-}z+(A_{1}IC_{0}(kr)+B_{1}I_{0}(kr))\exp i(kz-\omega t)$(7)$\phi \mathrm{o}(r, z, t)$$=$$B_{0}I_{0}(kr)\exp i(kz-\omega t)$.

.

はそれぞれ 0 次の第 1 種、第 2 種の変形ペツセル関数であり、ら決まる定数である。

(8)$A_{1}.,$$B$: は境界条件か250
運動学的境界条件は.

at$\frac{\partial\phi_{i}}{\partial z}\frac{\partial R_{i}}{\partial z}=i\partial\phi\partial r$$\frac{\partial R_{i}}{\partial t}[perp]$$r=R_{i}(z,t\grave{)}$(.9)$i=1,2,3$であり、 また力学的境界条件は境界面での表面張力を考慮して境界面での圧力差 \Delta 乃が

ュp:

$=$

$\sigma:-1,i(-\frac{7z\partial^{2}K\tau}{(1+(^{\partial}\neq_{z})^{2})^{(3/2)}}.+.\frac{1}{R\sqrt{1+(^{\epsilon_{\partial z}}-\mathrm{A}_{\mathrm{i}})^{2}}})$

at r=I ら (z, ) $:=1,2,3$

(10)

$t$

と表せるということを表す。ポテンシャル流であるので各領域での圧力はベルヌーイの定理から(11)$p:=- \frac{\rho}{2}.\cdot|\nabla\phi:|^{2}-\rho:\frac{\partial\phi}{\partial t}.\cdot$と書ける。

境界条件 (9) (10) に速度ポテンシャル (5)(6) $(7)(8)$ およひ境界の形 (3) を代入し
$A_{2},B_{2},A_{3},\eta_{1},$$\eta_{2},\eta_{3}$$B_{0,}.A_{1}.B_{1}$,が微小量として線形化するとこれら 9 個の変数に対する次のような連立方程式

を得る。
$C$

(

$B_{0},$ $A_{1},$ $B_{1},$ $A_{2},$ $B_{2},$ $A_{3},$ $\eta_{1}$, 砲, 男)T $=0$(12)

ここで係数行列 $C$ の 0 でない係数は$C_{11}=C_{22}=-kK_{1}(kR_{30})$$C_{27}=C_{38}=i(.\omega-V_{2}k)$,$C_{47}=i(\omega-V_{1}k)$,$C_{17}=C_{69}=$ 劫,$C_{23}=kI_{1}(kR_{30})$$C_{32}=C_{44}=-kK_{1}(kR_{20})$ ,,$C_{54}=-kK_{1}(kR_{10})$ ,$C_{33}=C_{45}=-kI_{1}(kR\underline,0)$$C_{55}=C_{66}=kI_{1}(kR_{10})$$C_{59}=i(\omega-\mathrm{r}_{1}^{r}.k)$,$C_{71}=-\dot{i}\rho_{3}K_{0}(R_{30})\omega$$C_{73}=\dot{r,}\rho_{2}I_{0}(Rs\mathrm{o})(\omega-V_{2}k)$,$C_{8_{\vee}^{\circ}}=-i\rho_{2}K_{0}(R_{20})(\omega-V_{2}k)$$C_{84}=-i\rho_{1}K_{0}(R_{20})(\omega-V_{1}k)$$C_{94}=-i\rho_{1}I\mathrm{f}_{0}(Rl.0)(\omega-V_{1}k)$$C_{77}=- \sigma_{23}(k^{2}-\frac{1}{R_{30}^{2}})$

,

,

,

,

,

,

$C_{\acute{7}2}=i\rho_{2}K_{0}(R_{30})(\omega-V_{\sim},k)$,$C_{83}=-i\rho_{2}I_{0}(R_{20})(\omega-V_{2}k.)$$C\mathrm{s}\mathrm{s}=-i\rho_{1}I_{0}(R_{20})(\omega-V_{1}k)$$C_{96}=-i\rho_{1}I_{0}(R_{10})(\omega-V_{1}k)$$C_{88}=- \sigma_{12}(k^{2}-\frac{1}{R_{20}^{2}})$

,

,

,

.

,

$C\mathrm{o}\mathrm{s}=i\rho_{0}I_{0}(R_{10})\omega$,$C_{99}=- \sigma_{01}(k^{2}-\frac{1}{R_{10}^{2}})$で与えられる。 係数行列 $C$ の行列式が 0 という条件から分散関係は$\det C=f(\omega, k;R_{10}, R_{20}, R_{30,01}\sigma,\sigma_{12},\sigma_{23},\rho 0,\rho_{1},\rho_{2},\rho_{3}, V_{1}, V_{2})=0$となる。 これはることとし、$k$$\omega$

に対する実係数 6 次式になっていることが分かる。 ここでは時間的不安定を考え

を実数、$\omega$を複素数として分散関係を $\omega=\omega(k)$ の形で解く。 このとき、 固有モードの位相速度は &(\mbox{\boldmath $\omega$}(k))/k、時間不安定の成長率はある$\omega(k)$(13)に対応したモードの振幅比$\eta:/\eta_{j}$${\rm Im}(\omega(k))$

で与えられる。 また分散関係の解でを計算する事によって固有モードの形状を知ることが251

図 2: 波数22$k$に対する固有値の実部及び虚部 (不安定成長率)

層$:(\mathrm{a})1$$\acute{(}\mathrm{b}$

) 層
$2$

2 層円環形液体シートの固有モードの解析

以下の 2 層円環形液体シートの固有モードの解析では周囲の流体の効果を無視し $\rho_{3}---\rho_{0}=0$ とおく。 また 2 層のシートに速度差があると Kelvin-Hebnholtz 型の不安定が現れるが、 ここでは表面張力不安定に重点を置くこととし 2 層のシートの速度差がない場合 (以下の解析では物理量の基準と で,’$V_{1}=\underline{V_{2})\text{を}考}-$$\text{し}$$v_{\mathrm{r}\mathrm{e}\mathrm{f}}=\sqrt{\frac{2}{\rho_{1}}\underline{\sigma}_{R_{10}}\Omega\llcorner}$$p_{\mathrm{r}\mathrm{e}\mathrm{f}}=\mathrm{A}_{\mathrm{L}}\sigma$$t_{\mathrm{r}\mathrm{e}\mathrm{f}}=\sqrt{}^{\mathrm{B}\mathrm{L}\lrcorner\alpha}2\sigma R\rho \mathrm{l}$”える。 また、をとり、 これ$1$を用いてその他の物理量を無次元化してぃる。
代表的な場合として $R_{20}=1,$ $b_{10}=b_{20}=0.2,$ $\sigma_{01}=\sigma_{12}=\sigma_{23}---1,$ $\rho_{2}=\rho_{1}=1$ を考え、これを refer.ence state と呼ぶことにする。 この場合の安定性解析の結果を図$b_{10}=R_{20}-R_{10}$まず図 2 は$\omega$,$b_{20}=R_{30}-R_{20}$は定常流での各$\mathrm{i}$の厚さを表すパラメータである。

の $k$ に対する依存性を示したものである。 実線が Im(\mbox{\boldmath $\omega$})、いる。 この図より 3 個のモードが存在し、それぞれのモードは部の大きい順に$\omega_{1},\omega_{2},\omega_{3}$とする。 ) また各

し、 k=ki。が不安定となる最大の大きい。すなわち$\omega_{2},\omega_{3}$$\omega_{1}$$k$2 に示す。 ただし

$k$

破線が

${\rm Re}(\omega)$の小さい領域で不安定である。 (虚はとなる k=k,。で純虚数から実数へ変化を与える。 さらに不安定領域ではは. より 10 倍程度$\omega_{i}$$\omega=0$$\omega_{1}$$\omega_{2},$$\omega_{3}$の表すモードの不安定は他のモードの不安定よりかなり強い。 または 0 に近づくが、$\omega_{1}$を表して$karrow \mathrm{O}$では 0 でない値に近づくということがわかる。

ここで不安定成長の機構と不安定を与える最大波数 k,。の関係を考えてみる。不安定を与える最大波数 ki。では分散関係の解は $\omega=0$ になってぃるが、 これは係数行列 $C$ の要素の中で$C_{77},$ $C_{88},$ $C_{99}=\sigma_{i-1i},$という関係が満たされる時である。 ここで$F_{1}\sim\sigma_{5^{\frac{R}{z}}\sim}^{\partial^{2}}\neq$

表面張力

(14)

$(k_{ic}^{2}- \frac{1}{R_{i}^{2}},)=0$$)1$$\sigma_{i.-1i}(k_{\tilde{i}c}^{9}-\overline{R}_{\backslash }I$.

の第 1 項は第 1 曲率にょる表面張カを表し、 変位に対して復元力として安定化に寄与する。 また第 2 項は第 2 曲率にょる$F_{2} \sim\frac{\sigma}{R_{i}}$を表し、 変位に対して変位を促進させる方向に働き不安定化に寄与する。 すな

わち不安定を与える最大波数 ki。は界面での第 1 曲率にょる表面張 ] と第 2 曲率にょる表面張 fJ が$f$

つりあうような波数である。
次に 2 層円環形液体シートの固有モードの形状を振幅比面での振幅$\eta_{2}$$\eta_{1}$と外側の界面の振幅$\eta_{3}$$\eta_{i}/\eta_{j}$を用いて調べる。 図 3 は内側の界の比及ひ内側の界面での振幅に対する中心界面での振幅の関数として図示したものである。
に対応する最も成長率が大きなモードはがすべの比を図 2 の場合と同じ reference state の場合に$\eta_{1}$$k$付近での振幅を見てみると、て同符号であり 1 層の場合の反対称モードに対応すると考えられる。 また$k\sim \mathrm{O}$$\omega_{1}$$\eta 1,\mathit{7}|2,$$\omega_{2}$$\eta \mathrm{s}$に対応する 2 番目に

252

図 3: 振幅比$\eta \mathrm{s}/\eta_{1}$及ひ$\eta_{2}/\eta_{1}$反対称モード内部モード

図 4: 2 層円環形シートでの固有モード

成長率の大きなモードは内側と外側の変位が互いに反対方向であり、それらに対して中心界面の変位が小さく 1 層における対称モードに対応すると考えられる。
に対応する最も或長率が小さな$\omega_{3}$

モードは外側と内側の界面の変位が同符号であり、それに対して中心界面の変位が逆符号で外側と内側の変位に対して比較的大きいというモードである。 これを内部モードと呼ぶことにする 図に近くなると各モードの形状はが大きくなり不安定の限界であるのものとは異なっ$4)_{\text{。}}$

$($

$k$

てくる。 例えば$k\sim \mathrm{O}$$k_{\dot{\alpha}}$$\omega_{2}$に対応するモードは k=k。
。付近では内部モードのような振舞をしていることがわかる。 またさらに$k$が k:。を越えて大きくなると例えば$\omega_{1}$に対応するモードは $k\gg k_{1e}=\pi^{\frac{1}{10}}$では内部モードのような振舞いを示し、 またこのような領域では$\omega_{3}$に対応するモードが反対称モードの振舞いを示す。 このようなモードの入れ替わりは 1 層シートの場合でも起こることがわかっている。

次にパラメータ変化に対する不安定成長率の変化を調べる。 図 5 は 2 層の厚さの比 $b_{20}/b_{10}$ の変化 (図 5-b)、中心界面と内側界面の表面張力係数の比の変化 図

層の密度比内の変化 (図 5-d)、に対して不安定或長率の大きさが各モードでどのように変わるかを横軸に をと$\sigma_{12}/\sigma_{01}$$($$5-\mathrm{c})_{\text{、}}2$$/\rho_{1}$

$k$

253

図 5: パラメータの変化に対する不安定成長率の変化(a) 1 層、 (b) 2 層の厚さの初期値の比を変化、 (c) 中心界面と内外界面の表面張力係数の比を変化、(d) 2 層の密度の比を変化

り縦軸に成長率を表したグラフで示したものである。 また図 5-a として 1 層シートの振幅比を比較のために示している。 図 5-b からは $b_{20}/b_{10}$ を 1 より大きくしていくと反対称モードの或長率は減少するものの、極大に関しては $k=0$ 以外で成長率の極大が現れることがわかる。 これはパラメー
タの変化によって成長率の極大がは $k\neq 0$ になるということを示しており興味深いと考えられる。
で 1 層の場
を大きくすると全体として成長率が増加すること、
図 5-c からは$\sigma_{12}/\sigma_{01}arrow 0$$\sigma_{12}/\sigma_{01}$合に漸近することがわかる。 図 5-d からはの成長率が減少するのに対し、$\omega_{1}$3$\rho_{2}/\rho_{1}$$\rho_{2}/\rho_{1}$を 1 より大きくするとを小さくすると$\omega_{2}$$\omega_{3},\omega_{1}$に対応するモードに対応する対称モードの成長率が増加し、に対応する反対称モードの成長率については $k=0$ 以外で成長率の極大が現れることがわかる。

薄膜近似の 2 層円環形シートへの適用
薄膜液体シートの形状の時間発展を考える時は、 非線形効果が重要であると考えられる。 薄膜

シートでは界面が時間的空間的に変化する自由境界面であり、 直接数値計算によるシミュレーショ

ンは困難である。そのために系に近似を行って現象をうまく記述できるモデル方程式を導き出すことが重要である。 1 層円環形シー } では薄膜の厚さ方向にナビエ・ストークス方程式を積分した平均化方程式を導く手法 (薄膜近似) が有効であることが知られている [ f7]。 ここでは [8] に従つて$\backslash$$6^{1}\rfloor$

254
薄膜近似を非粘性流体からなる 2 層円環形シートに適用し非線形効果が重要な領域の振舞いを記述できるようなモデル方程式を導き、そのモデル化の正当性を線形安定性解析の結果がポテンシャル
解析の結果と合っているかを調べることで確認する。

基礎方程式は連続の式とオイラー方程式である。 周囲の流体の効果を無視し、 2 層の流体の基
礎方程式を各層ごとに円環の厚さ方向に平均化するために$R_{2}(z,t)$$r$をかけて $R_{1}(z, t)$ から$R_{2}(_{\backslash }z_{\backslash }t)$

又はから $R_{3}(z,t)$ まで積分する。

$\int_{R_{i}}^{R}\cdot.+1\frac{\partial}{\partial r}(rv_{r}^{(\dot{l})})dr+\int_{R-}^{R}\cdot.+1\frac{\partial v_{z}^{(i)}}{\partial z}rdr=0$$I_{R}^{R_{\mathrm{i}\star 1}} \dot{.}\frac{\partial v_{z}^{(i)}}{\partial t}rdr+\int_{R_{i}}^{R}\cdot.+1\frac{\partial v_{z}^{(i)2}}{\partial z}rdr+\int_{R}^{R_{+1}}.\frac{\partial}{\partial r}(rv_{z}^{(1)}.v_{r}^{(i)})dr=-\int_{R\mathrm{a}}^{R_{j+1}}\frac{1}{\rho_{i}}\ovalbox{\tt\small REJECT}^{i)}dr$$\int_{R}^{\alpha_{+1}}.\frac{\partial v_{r}^{(\dot{l})}}{\partial t}rdr+\int_{R}^{R\mathrm{s}+1}\dot{.}\frac{\partial}{\partial \mathrm{r}}(rv_{r}^{(1)2}.)dr+\int_{R}^{R_{+1}}$$\cdot\frac{\partial}{\partial z}(v_{z}^{(\dot{l})}v_{r}^{(\cdot)}.)rdr=-\int_{h}^{R+1}\frac{1}{\rho}.\cdot\frac{\partial p^{(\dot{l})}}{\partial r}rdr$

.

(15)

(.16)
(17)

境界条件は運動学的境界条件と界面での圧力の差と表面張力の関係を与える力学的境界条件である。
実際に平均化された方程式を導くためには厚さ方向の各物理量の関数を決定する必要がある。そのために「シートの厚さがその上に生じる擾乱の波長に比べて小さい」という仮定をおく。 このことを物理量に関する仮定として次のように表す。 これが薄膜近似の考え方である。

.
.

$v_{z}$は$r$

方向 (シートの厚さ方向) に一定である。

$v^{(\dot{l})},.(r, z, t)=v_{z}^{(\dot{l})}(z,t)$$v_{r}^{(\dot{l})},p_{i}$の $r$ 方向の分布形は一次関数$v_{r}^{()} \dot’(r, z, t)=..\frac{v_{r}^{(1)}(R_{+1})+v_{r}^{(\dot{l})}(R_{+1}\dot{)}}{2}.+\frac{v_{r}^{(i)}(R_{i+1})-v_{r}^{(1)}(R_{t})}{R_{i+1}-R_{1}}.\cdot(r-\cdot\frac{R_{+1}+R_{i}}{2})$$p^{(:)}(r, \approx, t)=\frac{p^{(i)}(R_{t\neq 1})+p^{(i)}(R_{+1})}{2}.+\cdot\frac{p^{(\dot{\mathrm{s}})}(R_{+1})-p^{(1)}(R)}{R_{+1}-R_{i}}.\cdot.(_{\backslash }r-\cdot\frac{R_{+1}+R}{2}.\cdot)$

円環の外側およひ内側での圧力は境界面ての境界条件より表面張力による圧力差から決まる。 1 層の場合にはこの値のみで圧力分布の関数形を決める事ができる。 しかし 2 層の場合には圧力分布の関数系を決めるためには中心界面での圧力差だけでなく、 中心界面での圧力の値を決める事が必要である。 したがって中心界面の圧力を表す $P(z, t)$ を導入する (図 6)。 これを中心界面での内側の層での圧力と外側の層での圧力の平均値とみなす。 すなわち
$P(z,t)= \frac{1}{2}(p^{(1)}(R_{2})+p^{(2)}(R_{2}))$

$P$

(18)

をもちいると各界面での圧力は$p^{(1)}(R_{1})=p_{0}-\Delta p_{1}$,$p^{(1)}(R_{2})=P+ \frac{\Delta p_{2}}{2}$$(19\grave{)}$$p^{(2)}(R_{2})’=P- \frac{\Delta p_{2}}{2}$,$p^{(2)}(R_{3})=p_{3}+ \frac{\Delta p_{2}}{2}$(20)基礎方程式を簡単にするために物理量 $f(f=v_{z}^{(1)}., v_{r}^{(1)}.,p^{(1)}.)$ に対し厚さ方向に平均化された量を定義する。
.

$\overline{f}=.\cdot\frac{\int_{\alpha^{+1}}^{R}frdr}{\int_{R_{5}}^{R}+1rdr}$

(21)

255

図 6: 薄膜近似における圧力分布関数

薄膜近似の仮定の下での平均量を求めると$\overline{v_{z}}^{(i\rangle}(z, t)$$=$$\overline{v_{r}}^{(i)}(z, t)$$=$$\overline{p}^{(i)}(z, t)$$=$となる。 これらの平均量$v_{7,}^{(i)}(z, t)$(22)$\frac{(v_{r}^{(i)}(R_{\dot{\iota}+1})+v_{r}^{(i)}(R_{i}))}{2}+\frac{1}{6}\frac{R_{i+1}-R_{i}}{R_{i+1}+R_{i}}(v_{r}^{(i)}(R_{i+1})-v_{r}^{(\dot{|})}(R_{i}))$(23)$\frac{(p^{(i)}(R_{i+1})+p^{(i)}(R_{i}))}{2}+\frac{1}{6}\frac{R_{i+1}-R_{i}}{R_{i+1}\cdot+R_{i}}(p^{(i)}(R_{i+1})-p^{(i)}(R:))$(24)$\mathit{5}^{(i)},\overline{v_{r}}^{(i)},\overline{p}^{(i)}$.

$\Delta p^{(i\mathrm{I}}$及び$=p^{(i)}(R^{i+1})-p^{(i\mathrm{I}}(R^{i}.)$$\Delta v_{r}^{(i)}=v_{r}^{(i)}(R_{\grave{J}}^{1+1}.-v_{r}^{(i)}(R^{i}.)$を用いて基礎方程式を書き直すことにより、おける各層の圧力$v^{(1)}z’ v^{(1)}z’\overline{v^{(1)}r},\overline{v^{(_{\sim}^{i})}r.},$ $\Delta v^{(1)}r’\Delta v^{(2)}r’ F_{\mathrm{t}}1,$$R_{3}\mathrm{D}$$\underline{|},,2,$$p^{(2)}(R_{3}),p^{(2)}(R_{2}),p^{(1)}(R_{2}),p^{(1)}(R_{1})$で書かれた 2 層円環$\sqrt[\backslash ]{}\backslash$と界面に–} のモデル方程式$\cdot$が得られる。 これは空間 1 次元の時間発展方程式であり、一般的な方法で解ける方程式であることに注意する。

この方程式を導くに当たって薄膜近似の仮定をし、 中心界面での圧力 $P(z, t)$ という変数を導入した。 このようにして導かれた方程式の正当性を調べるためにこの方程式について時間的不安定の安定性解析を行いポテンシャル解析による結果との比較を行う事にする。 変数は$\mathrm{I}J_{r}(\overline{(1)}z, t),$

る波数

$k$

$\overline{v_{r}^{(}.}(z, t)|2),$.

$\Delta v^{(1)},,.(z, t),\Delta v_{r}^{(2)}(z, t),R_{1}(z, t),R_{2}(z, t),R_{3}(z, t),P(z, t)’$$v_{z_{\backslash }}^{(1)_{(}}z,$

$+.$

),

$v_{\approx}^{(2)}(z.t)$,てあり、 定常流に対すの正弦波的な微小擾乱として$R_{i}(z, t)=R_{i0}+\eta_{i}.\exp i(\omega t-kz)$$v_{z}^{(i)}(z, t)=A_{z}^{(1)}.\exp i(\omega t-kz)$,$\overline{v_{r}^{(i)}}(z,$,$t1,=A_{f}^{(i)}\exp i(\omega t-kz),$$\Delta v_{r}^{(i)}(z, t)=A_{\Delta v}^{(i)}\exp i(\omega t-kz)$$P(z, t)= \frac{\sigma_{23}}{R_{3\mathit{0}}}+\frac{1}{2}\frac{\sigma_{12}}{R_{20}}+A_{P}\exp i(\omega t-kz)$

(25)

(26)
(27)
(28)

256
を仮定し、 これを基礎方程式に代入し擾乱の係数$\eta_{\ovalbox{\tt\small REJECT}},AS^{)},$),$A\ovalbox{\tt\small REJECT}$$A\ovalbox{\tt\small REJECT},$$A_{P}$

が微小量として線形化し係数行列の行列式が 0 である条件を求めると
d 杭 A $=f$ ($\omega,$$k;R_{10},$ $R_{20},$ $R_{30},$ $\sigma_{01},$ $\sigma_{12},$ $\sigma_{23},$がであり、 これが分散関係になる. これは$\omega$$\rho_{1}$, 内,$V_{1},$ $V_{2}$) $=0$

についての実係数の 6 次式である。 これを(29)$\omega$について解くことで固有モードの位相速度 &(\mbox{\boldmath $\omega$})/k、時間不安定の成長率を求める事ができる。 こ
れらをポテンシャル理論によるものと比較し 2 層系に対する薄膜近似の妥当性を (微小擾乱の範囲${\rm Im}(\omega)$内) で検証する。

ence state$\mathrm{k}$${\rm Re}(\varpi)$$\mathbb{R}\alpha)$

図 7: ポテンシャル理論と薄膜近似による固有値の比較:

係数$\text{、}$state (内外の層で表面張カ厚さの初期値、密度、 の比がすべて 1 の場合) (b) 密度の比が O.ml の場合$(\mathrm{a})\mathrm{r}\mathrm{e}\mathrm{f}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}$図 7 はポテンシャル理論及ひ薄膜近似によって分散関係の解を求め、不安定成長率 (固有値の虚部) 及ひ位相速度を波数で割ったもの (固有値の実部) を横軸を波数$k$にとってグラフにしたものである。 これから分かるように不安定成長率はパラメータによらず非常によく一致する。具体的には位相速度は$k$

が $k\sim 5$ のあたりまでは非常に良く一致している。 ここでは $R_{20}=1,$ $b_{10}=0.2$ よ

り $k\sim 5$ は波長が薄膜の厚さ程度を意味しており、 これは薄膜近似の仮定とあっている。

4

結論
ここでは 2 層円環形液体シートの安定性について考えた。 1 層液体シートで生じる対称モード、

反対称モードについては 2 層円環シートにおいても生じ全体として同様の変形をシートにもたらす。その不安定成長率の比は 1 層円環形液体シートでそうであったように反対称モードの成長率が対称モードの成長率より大きくなっている。 反対称モードは 1 層の反対称モードが 2 つある場合、

257
対称モードは中心界面の表面張力が効果を及ぼしていない場合に対応すると考えられるのでこれは自然な結果と考えられる。

2 層シートでは 1 層シートでは現れなかったモードとして中心界面がより内外の界面に比べてより変動する内部モードが現れる。 しかしこの不安定成長率は対称モード、反対称モードと比較して弱い。

不安定を与える最大の波数 k:。に注目すると、 この値の付近ではモードの形状はとはかなり異なっている。不安定を与える最大の波数$k$$k\sim \mathrm{O}$のあたりはある界面の変位を増加させる方向に働く

第一曲率による表面張力と、 界面の変位を減少させる方向に働く第二曲率がつりあうような波数であり、対応するモードは対応する界面の変位が主に変位しているようなモードになっていると考えられる。 したがって固有モードの形状はから k\sim k:。にかけて変化していくと考えられる。
中心界面における圧力を変数に導入する事で薄膜近似を 2 層円環形シートに適用し、近似方程式$k\sim \mathrm{O}$を得た。波長がシートの厚さ程度の波長より大きいという近似の範囲内て、位相速度と不安定成長
率がポテンシャル理論の結果と一致する事を確認した。
パラメータ変化に対する不安定成長率の変化の振舞いや、実際のシートの崩壊過程を調べるためには非線形効果を取り入れた数値シミュレーションが必要であると考えられる。薄膜近似による方
程式は 1 層の場合と同様 2 層の場合においても数値シミュレーションに有効であると考えられる。

参考文献[1] W. A. Sirignano and C. Mehring, “Review of theory of distortion and distintegration ofliquid streams” Prog. Energ. Comb. Sci. (2000), 609 and references therein.

[2] H. B. Squire, “Investigation of the instability of moving liquid film”Brit.J.Appl.Mech. 4(1953),167.
[3] G. D. Crapper, N. Dombrowski and G. A. D. Pyott, “Kelvin-Helmholtz wave growth on
cylindrical sheets” J.Fluid Mech. 68 (1975) ,497.

[4] J. H. Dumbleton and J. J. Hermans, “Capillary Stability of aHollow Inviscid Cylinder”
Phys Fluids 13 (1970), 12.
[5] R. H. Rangel and W. A. Sirignano, Phys$.Fl\mathrm{u}ids\mathrm{A}3(1991),2392$ .
[6] T. Yoshinaga and K. Kotani, Lect RIMS, Kyoto Univ $.(2000),28$ [in Japanese][7] C. Mehring and W. A. Sirignano, “Axisymmetric capillary waves on thin annular liquid
sheets” Phys Fluids 12(2000),1417.
[8] J. I. Ramos, “Annular Liquid Jets: Formulation and Steady State Analysis” ZAMM
72(1992),565.

[9] J. M. Kendall, “Experiment on annular liquid jet instability and on the formation of liquid
shells”Phys. Fluids 29(1986),2086.