数理解析研究所講究録 1271 巻 2002 年 112-124

112
レイリー・ テイラー不安定性による液滴の分布村上洋一 (Youichi Murakami)
シーディーアダプコ 飯田和雄 (Kazuo Iida)

阪府大工

1

はじめに

自由境界をもつ液膜が絡む問題は工学や自然界に数多く存在する. 例えば, 胃酸から胃壁を守るためにある種の液体が粘膜から分泌され膜を形成している. これがうまく出ないと胃潰瘍になってしまうそうだ. このことを参考にしたのかどうかは知らないが, タービ

ン翼を熱から保護するために薄い空気の層を壁近傍で流すことが行われている. まぶたがスムーズに開閉できるのも眼球に薄い涙の膜があるためである. これが不十分だとドライアイになるのだろうか. このように摩擦の低減に利用されることもある. この講究録が書ける (読める) のも液膜の働きかと思うと, 何だかありがたい気持ちになる. さまざまな
地学, 生物物理学や化学の応用例については, Oron et al. [1] にまとめてある. 種々の問
題に液膜は関連しているので, 流体物理学の一分野として液膜の運動をとりあげることは,各分野の基盤的知識を高めることにつながるはすである. 液膜の一様な状態がどのような状況で維持され, 不安定性によりどのような変化が生じるのか. これらを明らかにすることが主な課題となる.
ここで考える非常に薄い液膜では壁面での粘性の効果が膜全体に及ぶため, 非粘性を仮定できない. そのため, この力学は保存系ではなく散逸系である. 保存系の揚合は, 1 次
元的に局在する孤立波はよく観察されるが, 2 次元的に局在する孤立波はあまり観察されない. 一方, 散逸系の場合は 2 次元的に局在した構造が観察されることが多い. 斜面を下る流れにおける馬蹄形の波 [2], レイリー. テイラー不安定性における液滴 [3] や散逸の強い場合のファラデー波における 「指」 に似た構造の波 [4] などがその例である. このように, 本質的に 3 次元の流れ場とその 2 次元表面変形を取り扱う必要がある点が散逸波動の1 つの特徴と考えられる.
通常, 非圧縮のナヴイエ. ストークス方程式が基礎となりそれにさまざまな効果を付加し適切な境界条件を与えれば, 現象を記述することができる. しかし, このような流体の基礎方程式を表面変形を考慮して数値的に解くことは原理的には可能であるが, 非常に時間がかかる. 例えば, 臨界値を求めるような問題では 2 次元計算でも莫大な時間がかかる[5]. L たがって, 流体層が薄いという点を考慮して基礎方程式を簡略化することが適切である. この近似をすることにより表面変形 $h=h(x, y, t)$ のみの方程式を導出し, それを解析的およひ数値的に解くニとが考えられる. 一般に方程式に非線形かつ減衰を表す高階微を非常分が含まれるため, 数値計算が困難になる. 時間について陽的に解くと時間刻みに小さくとる必要があり, 現実的ではない. 陰解法を用いた 2 次元の数値シミュレーショ
ンはごく最近になって行われるようになった. 例えば, Sharma and Khanna[6] は, ファ$\Delta t$

ン・デル. ワールスカを考慮した場合について $30\cross 30$ およひ 60$\cross 60$の分割数の数値シ

113
ミュレーションをギアの方法を用いて行ってぃる. Oron[7] はマランゴニ効果のある揚合について $51\cross 51$ およひ $81\cross 81$ の分割数を用いてシミュレーションしてぃる. このように, あまり大きなサイズの数値シミュレーションは行われてぃない. 局在構造の空間分布を議論するような揚合はより大きなシステムサイズの系を扱うことが望ましい.
この研究では, 近似によって得られている 2 次元の長波長方程式を比較的大きなシステ

について数値的に解いた結果 (こつぃて報告する. 液膜につぃて長波長方程式や近似方程式のニレーションの研究は多い. しがし, 1 次元の揚合につぃ
ムサイズ$(128\cross 128)$$\backslash \grave{\nearrow}_{\backslash }^{\backslash }\backslash$ては直接数値シミュレーションとの比較もなされてぃるが, 重要となる 2 次元の場合については実験との比較があまりないようである. そこで, 長波長方程式の有効性を検討できるように室内実験のある場合を取り上げる. 幸いなことに, 薄い液膜に対して重力と表面張力のみが働くレイリー. ティラー不安定性に関する実験が, Fermigier et al.[3] にょって行われている. また, この系のポテンシャルは非常に簡単な形をしてぃるので, 定常解を解析的に求めることができるをいう利点もある.
彼らは, 直径の容器にの薄さのシリコンオイル層を作った. また, シリ
コンオイルは粘性率の大きなものを用い, 以下のことを見出した. (i) 一般に 1 点に局在した撹乱は同心円状に広がりながら発達し, 最終的には 6 回対称性を持っように軸対称な釣鐘状ピーク (液滴) が分布する. 1 これを 6 角形パターンと呼ぶ. (ii) 線状に撹乱を与え$56\mathrm{m}\mathrm{m}$$0.2\mathrm{m}\mathrm{m}$初期はその形状を保ちながら成長し回りに平行な線状のピークが広がってぃく. 線
状のピークが不安定になり局在した釣鐘状ピークが並ひ, 最終的には 6 角形パターンになる. (iii) 釣鐘状ピークの間隔は線形安定性理論の最大増幅率に対応する波数とほぼ一致す
ると,(iv) 時間が十分経過すると, ピークは雫となって落下する.
この現象を理論的に説明するために, Fermigier et al.[3] では, この長波長方程式を導き, その枠組みで振幅方程式を導出し非線形増幅率の観点から 6 角形パターンが生じることを議論している. Hammond[8], Yiantsios and Higgins[9] およひ Oron and Rosenau[10]る.

は 1 次元に限定した長波長方程式をもとに時間発展と定常解につぃて考察してぃる.

実験

結果からも明らかなように表面変形を 1 方向に一様と仮定することはできない. そこで,この研究では 2 次元の長波長方程式を数値的に解いた. Schwartz[ll] が, 差分法にょり予
備的な数値シミュレーションを行ってぃるが, 最終状態の 1 例を考えてぃるに過ぎず, その詳細については述べていない. 数値計算結果の詳細につぃては, 既に飯田と村上 [12] お
よひ村上と飯田 [13] で行い, 定性的およひ一部分は定量的に実験と一致することを確認している. また, 釣鐘状ピークの形状, 時間発展, 安定性, 融合条件等につぃてはそこで報告した. この研究では特に, 現れる釣鐘状ピークとその空間分布の特性を明らかにする.
まず, 長波長方程式の導出方法およひその簡単な性質を述べる. 数値シミュレーション
によって得られたパターンを, 実験との比較を交えながら与える. さらに, 釣鐘状ピークについて今まで得られている結果についてまとめる. 以上は, 村上と飯田 [13] で得られた結果をまとめたものである. 釣鐘状ピークの分布の特性をボロノイ分割 (ドローネ三角分
割) を用いて調べる. 最後にこの研究のまとめと今後の課題につぃて説明する.

1

タイトルでは液滴としたが,$\sim-\vee\vee$\mbox{\boldmath $\tau$}.IR 釣鐘状ピークと呼ぷことにする.

114
2

長波長方程式

2.1

問題の定式化

図 1: Schematic representation of the fluid film bounded by a$\mathrm{w}\mathrm{d}\mathrm{l}$and air.

固体表面に粘性流体を薄く塗った状態で生じるレイリー. テイラー不安定性をここでは扱う. 図 1 に示すような座標系をとり, T 向きを 座標の正方向としている. 密度 , 粘の層をなしている. 重力 $g=gz$ が T 方向に働き, 粘性流体
性率 の粘性流体が厚さと空気の界面では表面張力係数 に比例する表面張力が作用する. 空気の運動は一切考えず, 表面では応力自由の条件を課す. 厚さは位置 $r=(x,y)$ の関数 :は変形がないときの厚さである. なは表面変形であり,で表される. ここで,$z$$\rho$$h_{0}$$\eta$$\gamma$$h(r,t)=h_{\mathrm{O}}+\zeta(\mathrm{r}, t)$$h_{\mathrm{O}}$$\zeta(\tau, t)$

お図 1 において水平方向の速度成分が描かれているが, これは表面変形によって生じる弱い流れを示し, 最初から水平方向の流れを仮定しているものではない.

3 次元非圧縮粘性流体が満たす連続の式およひナヴイエストークス方程式は以下のように与えられる.
$\nabla\cdot u+\frac{\partial w}{\partial z}=0$,$\rho \mathrm{t}\frac{\partial u}{\partial t}+(u\cdot\nabla+w\frac{\partial}{\partial z})u\}=-\nabla p+\eta\nabla^{2}u$$\rho \mathrm{t}\frac{\partial w}{\partial t}+$

ここで,

( .
$u$

そ

$(u, w),$ $u=(u, v)$ , $\nabla=$ ($w \frac{\partial}{\partial z}$$\partial/\partial x$

,

)

(1)

,

$w\}=-\nabla p+\eta\nabla^{2}w+\rho g$$/\partial y$

(2)

.

(3)

) としてあり, 水平速度成分と垂直速度成は圧力ある.
次に境界条件について考える. 上面の固体壁 $z=0$ では, 次の粘性境界条件を課す.

分をわけている. ここで,$p$$u(x,y,0,t)=0$ ,$w(x,y,0,t)=0$ .

(4)下面の表面 $z=h(r,t)$ では, 空気の粘性率が液体に比べて非常に小さいという仮定のも

01 接線 9”i1116 条件$\eta\frac{\partial u}{\partial z}=0$

(5)

115
と法線応力のつりあいの条件(6)$P_{a}-p= \gamma\nabla\cdot(\frac{\nabla h}{\sqrt{1+(\nabla h)}})$

を課す. ここで,$P_{a}$

は大気圧である. さらに, 運動学的条件として,(7)$\frac{\partial h}{\partial t}+u\cdot\nabla h=w$が必要である. なお, 薄層が広がっている水平方向 ( $x,$を適用し, 側壁の影響を無視した.

2.2

$y$

方向) については周期境界条件

澗滑近似による発展方程式

低レイノルズ数の流れ (ストークス近似) と長波長不安定 (
, ここで
は波長
の厚さ,
は不安定モードの波長) を仮定する. 連続の式 (1) 上り, 、垂直速度成分 $w$ は$h\text{。}\ll\lambda_{M}$$h_{0}$$\lambda_{M}$水平速度成分 $u,$ と比べて非常に小さいことがわかる. このことから, 式 (2) と (3) において $w$ の項を落とす. さらに, 表面の勾配が小さいという仮定から法線応力のっりあいの$v$

条件において表面張力の項を簡単にする.
このような近似のもとで, 圧力, 水平速度をもとめ, 表面での運動学的条件に代入する
ことで, 以下の発展方程式を得る.

.

$\frac{\partial\zeta}{\partial t}+\frac{1}{3\eta}\nabla\cdot[(h_{0}+\zeta)^{3}\nabla(\rho g\zeta+\gamma\nabla^{2}\zeta)]=0$

ここで,

次に,

$h=h_{0}+\zeta(\mathrm{r}, t)$

(8)

である.

この式を無次元化するため代表長さを$h_{0}$

として, 式 (8) を書き直すと, 次のよ

うになる.
$\frac{\partial\zeta}{\partial t}+[\nabla\cdot(1+\zeta)^{3}\nabla(B\zeta+\nabla^{2}\zeta)]=0$$T= \frac{3h_{0}}{\gamma}$

ここで,

2.3

$B$

,

$B= \frac{pgh_{0}^{2}}{\gamma}$

.

,

(9)

(10)

はボンド数である.

長波長方程式の性質$\zeta\propto\exp(\sigma t+ikx)$と仮定し, それを式 (9) の線形項のみに代入しすると, 以 T の関係式が得られる.
$\sigma=-(k^{2}-\frac{B}{2})^{2}+\frac{B}{4}$

.

(11)

長波長の不安定が生じることがわかり, 方程式を導く際に表面勾配が小さいと仮定したことと整合性がある. また,
の最大値
のときの
は,
$\sigma$$\sigma_{\max}$$k_{M}=\sqrt{\frac{B}{2}}$$\sigma m\text{ }=\frac{B}{4}$

(12)

116
長波長方程式 (9) の全エネルギーの時間発展について考える.

(13)$B= \int(\frac{1}{2}(\nabla\zeta)^{2}-\frac{B}{2}\zeta^{2})dV$と定義し, 式 (9) を考慮すると,$\frac{\partial E}{\partial t}=-\int(1+\zeta)^{3}|\nabla(B\zeta+\nabla^{2}\zeta)|^{2}W<0$(14)となる. エネルギー $E$ は単調減少することを示しており, 時間に依存するアトラクターはありえないことがわかる.

次のような特別な場合を考えると, 式 (9) の定常解の条件を満たす.
$B\zeta+\nabla^{2}\zeta=\alpha$

.

(15)

は任意定数である. このヘルムホルツの方程式の解は自動的に定常解の条件を満たす. 軸対称解は次のように 0 次のベッセル関数で表すことができる.
ここで,$\alpha$$\zeta=AJ_{0}(\sqrt{B}r)+\zeta_{0}$

.

(16)

で最小値をとるとすると,は任意定数である.
時間発展のシミュレーション結果において釣$AJ\mathrm{o}(r_{\min})+\zeta_{0}>-1$ を満たす必要がある.
鐘状のピークが現れるが, $0<r<r_{m-n}$ の範囲で定常軸対称解で近似できることがわかる. 興味深いことに, この円領域の体積を求めると,ここで,$r=\prime_{m}:_{\hslash}=\sqrt{B}r_{1}$$A,$ $\zeta 0=\alpha/B$(17)$2 \pi\int\zeta rdr=2\pi r_{m1n}^{2}.\zeta 0$のように振幅 $A$ に依存しない. したがって,$\zeta_{0}$を保ったまま $A$ が変化したとしてもこの円領域での流量は保存したままである.

3

数値シミュレーション結果

今回用いたパラメータが表 1 に示されている. Fermigier et al.[3] が用いたものを採用となり, 長波長近似している. ここで \lambda M=2\pi 酩/\sim \approx 13.2mm より,の仮定は満たされている. また, $B=0.018$ である.
$h_{\mathrm{o}}/\lambda_{M}\approx 0.0152$$\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{c}\mathrm{a}1\mathrm{v}\mathrm{d}\mathrm{u}\mathrm{e}0.0002$$\mathrm{d}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{s}\mathrm{i}\mathrm{o}\mathrm{n}m$$\frac{\mathrm{p}\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{r}}{\mathrm{t}\mathrm{y}\mathrm{p}\mathrm{i}\mathrm{c}\mathrm{a}1\mathrm{t}\mathrm{h}\mathrm{i}\mathrm{c}\mathrm{k}\mathrm{n}\mathrm{a}\mathrm{e}\mathrm{s}}$h。

visicosity yy
density
kinematic viscosity$\rho$

surface tension$\gamma$

$\nu$

1.0
970
0.001031
0.021$kg/m\cdot s$$kg/m^{3}$$m^{2}/s$$N/m$

表 1: Definitions and values of parameters.

117
数値計算方法については, 空間についてはフーリエガラーキン法を, 時間発展につぃては前進オイラー法を用いた. また, 最小の表面厚さが初期厚みに対して $\zeta=(-1.0+10^{-4})h_{0}$にまで達したとき (液膜の厚みが 0 に近づくとき) に壁に付着するとみなし, 計算を止めた. 空間領域 $[Lx2\pi/k_{M}, L_{\mathrm{Y}}2\pi/k_{M}]$ , $(Lx=L_{\mathrm{Y}}=8)$ とし, 最大増幅率に対応する波長$\lambda_{M}=2\pi/k_{M}$ の長さを16 分割した. $dt$ は振幅の成長とともに変化させた.
実験と関連して, 初期条件としては,方向に伸ひたロールがトタンの屋根のように並んだもの (ロール撹乱) ,点に局在したガウス分布 (1 点局在撹乱) , (iii) 線状に小さな丘を与えたもの (線状撹乱) , (iv) 2 本の直交する線状に小さな丘を与えたもの (直交線状撹乱) 以上の 4 っの揚合には計算領域全域に微小なランダムを加える. (v) ランダ$h_{\mathit{0}}$$(\mathrm{i})1$$(\mathrm{i}\mathrm{i})1$

ムのみ (ランダム撹乱) の合計 5 っを用いた.

計算スキームのチェックとして, 時間発展の初期段階で振幅の成長が線形最大増幅率に一致することとエネルギーが単調減少することを確認した.

3.1

表面のパターン

以上の初期条件において計算を行った結果が次の図 2 に示されてぃる. 最終的な状態は,釣鐘状の軸対称ピークが空間分布し, 1 次元性は保たれないことがゎがる. 表面変形の発達過程において釣鐘状ピークが成長して形成される. 十分に発達してぃない隣接したピー
クは互いに引き合い融合し, 1 っの釣鐘状ピークとして成長することがあった. なお, すべてのパターンは最終的に破断 (表面が固体壁に付着する) し, その破断時間の直前付近のパターンが示されてぃる.
図は微小なロールにさらに微小なランダムを加えた場合の時間発展が, 示されて$(T\approx 1.7\cross 10^{5})$$2(\mathrm{a})$ロールが発達しかなり大きくなった後に, 一様な方向に不安定となり釣鐘状ピーク
がほぼ等間隔で形成された. 隣の列とはピークの並ひ方の位相が異なり, 6 回対称性に極めて近い分布を示すことがわかる.
図は Fermigier et al.[3] が行った実験と同様の初期パターンを再現したものである.
彼らの実験では, 初期パターンで単独の微小なピークを, 液層表面にある小さな塵が原因
でできたものとしており, それと同じ状態をガウス分布にょり再現してぃる. 時間発展が進むと局所的な撹乱から同心円状に振幅が広がりながら成長してぃく. 中心の振幅の大き
い円環状の部分から分裂しながら軸対称ピークが形成されてぃった. 最初の円環の部分ではピークが 6 つ形成されているのは興味深い. なお, 破断時間はおよそ $T=2.45\cross 10^{5}$で, これは実験の約 2 時間とおおよそ一致する.
いる.

$2(\mathrm{b})$図は,実験で一本の長い針金で作った撹乱に似せた初期条件の時間発展を示してぃ
定性的なパターンは実験とよく似てぃることがゎがる.
図は, 実験で 2 本の長い針金を直角に交差させて作った撹乱に似せた初期条件の$2(\mathrm{c})$る.

$2(\mathrm{d})$時間発展を示している. この場合も定性的なパターンは実験とよく似ており, 軸対称ピー

クの正方形格子が他の場合と異なり形成されてぃる. この場合はピークの振幅が他の揚合よりも小さいという特徴がある. このように振幅のピークは初期条件にょらずに決定されるものではなく,全体の配置のような他のピークとどのように隣接してぃるかに密接に関係しているようである.
図$2(\mathrm{e})$は,微小なランダムからの時間発展を示してぃる. それほどきれいに並んでぃる

118
わけではなく,6 回対称性に近いかどうか検討する余地がある.

図 3 には, 最終状態における $B\zeta+\Delta\zeta$ の分布が示されている. この値が定数であるということは定常であるということを示している. したがって, ほぼ一定値となっているところは定常で止まっていると考えられ, 変化しているところが非定常性を示している. 図で見られる円領域は釣鐘状ピークの位置と一致しており, 軸対称定常解で近似できることが示唆される. 初期条件に依らず円領域は各々同じ大きさで重なっているところはないことがわかる. この点について後に述べる.

119.

$.\cdot..\cdot\cdot..\cdot\cdot.\cdot\ddot{\mathrm{a}}^{1}..\cdot.\cdot\cdot.\cdot..\cdot.\cdot.\cdot.\cdot..\cdot..\cdot..\cdot...\cdot.\cdot..\cdot..\cdot.\ddot{\dot{\ovalbox{\tt\small REJECT}}}..\cdot.\cdot.\cdot.\cdot..\cdot.\cdot.\cdot.\cdot.\cdot..\cdot.\cdot...\cdot.\cdot..\cdot.\cdot.\cdot..\cdot..\cdot..\cdot.\cdot.\cdot\cdot\cdot.\cdot.\cdot.\cdot.\cdot.\cdot.\cdot..\ovalbox{\tt\small REJECT}.\cdot.\cdot.\cdot.\cdot..\cdot.\cdot.\cdot..\cdot.\cdot.\cdot\cdot.\cdot.\cdot.\cdot.\cdot\cdot.\cdot...\cdot.\cdot.\cdot...\cdot..\cdot..\cdot..\cdot\cdot\cdot.\cdot.\cdot.\cdot.\cdot.\cdot..\cdot..\cdot...\cdot.\cdot...\cdot..\cdot\cdot....\cdot..\cdot.\cdot\cdot$$(\mathrm{a})\mathrm{R}\mathrm{o}11,$$T=1.68\mathrm{x}10^{5}$$(\mathrm{c})\mathrm{L}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k},$.

$T=1.66\mathrm{x}10^{5}$$(\mathrm{b})\mathrm{G}\mathrm{a}\mathrm{u}\mathrm{s}\mathrm{s},$.

$T=2.46\mathrm{x}10^{5}$.

(d)Crosspeak, $T=1.32\mathrm{x}10^{5}$ .

(e) Random, $T=1.74\mathrm{x}10^{5}$ .
図 2: Final patterns of the surface deformation$\zeta$.

120$(\mathrm{a})\mathrm{R}\mathrm{o}11,$$T=1.68\mathrm{x}10^{5}$$(\mathrm{c})\mathrm{L}\mathrm{i}\mathrm{n}\mathrm{e}\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k},$.

$T=1.66\mathrm{x}10^{5}$$(\mathrm{b})\mathrm{G}\mathrm{a}\mathrm{u}\mathrm{a}\mathrm{e},$.

$(\mathrm{e})\mathrm{R}\mathrm{a}\mathrm{n}\mathrm{d}\mathrm{o}\mathrm{m},$$T=2.46\mathrm{x}10^{5}$$(\mathrm{d})\mathrm{C}\mathrm{r}\mathrm{o}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k},$$T=1.74\mathrm{x}10^{5}$

Ik 3: Final statae of the term.

$T=1.32\mathrm{x}10^{5}$.
$B\zeta+\nabla\zeta$

.

.

121
4

軸対称定常解の性質
図 4 に示すように, 釣鐘状ピークは 0 次のベッセル関数で与えられる定常軸対称解の最. の大値と最小値の間でほぼ一致していることがわかる. 以前にも述べたが, $0<r<r_{m1n}$間での体積はちょうと 0 になる. また, この領域では B\mbox{\boldmath $\zeta$}+\Delta (が一定の値を取ってぃる.

Comparison between an axisymmetric peak and the steady solution. (a) xdirection, $nx=5,$ $ny=3,$ $\zeta_{0}=1.262$ .
-direction, $nx=5,$ $ny=5,$ $\zeta 0=1.465$ .
図 4:$(\mathrm{b})y$

時間発展が進み表面変形が進むと, 釣鐘状ピークを持っようになってくる. 釣鐘状ピー
クを特徴付ける 2 っのパラメータ $A$ とは時間とともにゅっくりと変化してぃる. ピー$\zeta_{0}$クの境界のところでは非定常性があるため, ピークはほぼ定常解に近似されっつ成長を続.
ける. このことは軸対称定常解が安定な構造であることを示唆してぃる.
との値は初期条件に依存しており, 特に一意的に決定されるわけではないことも示されてぃるが,
ここでは省略する.
定常軸対称解となる 0 次のベッセル関数を最大値がら最小値 (O<r<rm 一 まで切り$A$$\zeta_{0}$

取った部分領域での定常解の安定性を調べた. この場合はピークの構造よりも大きい長波
長の撹乱を除外して考えることになる. っまり構造の局所安定性を調べるこ ‘とに相当する.
数値的には, 円領域の外部を強制的に最小値に固定することでこの状態を実現した. このような状態を作った後に, 安定であることを確かめるために, 質量を保存する微小な撹乱を加え, もとの形に戻るかどうかを時間発展を行うことで調べた. 初期の軸対称定常解にもどることを $A=0.5\sim 2.3$ の範囲で確認した. このことは $A$ の値に依らす, また明らがにには依す, この軸対称定常解が安定であることを示してぃる. したがって, 安定性か$\zeta 0$

ら振幅 $A$ が決定されると考えるのは困難である.

計算結果から, 釣鐘状ピークが隣り合うものと融合し振幅を伸ばすことが確認できたが,その条件について考える. 詳細は省略するが, あるピークをベッセル関数で近似したとき,

2 つのピークの領域 $(0<r<r_{m\dot{*}n})$ が重なったときに融合が生じることがわかった. これは破断前の状態でをプロットした際, 定常な円領域が重なってぃないことと関$B\zeta+\nabla\zeta$

係している.

122
5

釣鐘状ピークの分布

初期条件がランダムな場合の釣鐘状ピークの空間分布 [こつ1 て考察する. これは人為的な撹乱を加えない最も自然な状況と考えられる. この\nearrow Д拭璽} , 6 角形\nearrow Д拭
ンに似ているようだがきれいなパターンではない. 波数空間で さまざまなモードが同円状に分布している. そこで, 実空藺で特徴付けることを考える. 釣鐘状ピークはすべて同じとみなして, その中心の位置の空間分布に着目する. 周期的に点力\leq 分布して y‘な\mbox{\boldmath $\nu$}‘とここでは,$\mathrm{I}$$r\llcorner\backslash \backslash$$|\mathrm{h}$きはある点に対してどの点が隣接する点であるかが不明確である. そこで,($\text{ト^{}\backslash }\backslash$

ドローネ分割
ドローネ分割と田 ボロ[14] を用いて隣接するピークを定義した.
分割の定義は以下のよう}こなる.
分割と相補的なものである. ボロ

凸一*三角分割)$\text{ノイ}$$\text{ノ・イ}$

平面上に 個の点の集合 $P=\{p_{1},p_{2}\cdots,p_{n}\}$ があるとする. 各々の : に対し鯆蟲舛垢. これらの領て, 他の $p_{j},(j\neq:)$ エり p- の方が近い点の領域 R(p
$p$

$n$

域 $R(p_{1}),$ $R(p_{2}),$$\cdots,$$R(p_{n})$とをボロノイ図と呼ぶ.
点とボロノイ辺と呼ぶ.

$p_{j}$

は平面を領域に分割することになる. この分割のこを母点, ボロノイ図の頂点と辺をそれぞれボロノイ$P$
さらに, ボロノイ辺を共有する 2 つの母点を直線線分で結ぶと, 点集合 の三角形分割力 S得られる. これをドローネ分割七呼ぶ. 図 にその例が示されている. このようにして隣接ピークを決め,ビークを結ぶ直線の長さをピーク間の距離とした.

図 5: ドローネ分割の例

このようにして隣接するピークの個数のヒストグラムを作り, 規格化した結果が表 2 にまとめてある.

表 2: Probabilitiae of the numbers of adjacent peaks (or points).

123のピークの結果によると, 隣接数は 6 が一番多いが, 5 や 7 も各々の 20%
らいある. このように 6 角形パターンと呼ぶには問題があろう. また, サンプル数が 54 と少なかったため, 他の隣接数が 0 になった可能性は高い. 上り大きなシステムで同様の計算をする必要がある. 注目すべきことは平均隣接数がちょうど 6 になることである. ボロノ
イ分割 (ドローネ分割と双対関係にある) の 3 辺が 1 点で常に合うときは一般的に 6 にな
表$2(\mathrm{a})$$\langle$のることがオイラーの関係式を用いて示すことができる. これは一般的な結果で表一様乱数の揚合にも成立する. したがって, 平均が 6 だから 6 角形パターンと断定するのはよくないだろう. 比較のため, 一様乱数についての Tanemura[15] の結果が表に示されている. この場合は隣接数 6 のところで最大確率をとるが, 分布はなだらかに広がり
非対称である. したがって, レイリー. テイラー不安定性による釣鐘状ピークの分布は単純な乱数の分布とは異なる. これは何も不思議なことではなく, 乱数の揚合は 2 点はいくらでも近づくことができるため, 空間分布は釣鐘状ピークよりも点の散らばり具合が大きくなるからである. ピークの分布を再現するモデルについては, 最後の節で述べる.
を求め, 最大増幅率に相当する
このようにして定義した隣接ピークの距離の平均波長との比をとると,という結果になった. Fermigier et al.[3]が実験で求めた比率, (\lambda p/\lambda M)。xp\approx 1.19 と非常に近い値をとることがわかった. 彼らの$2(\mathrm{b})$$5(\mathrm{b})$$(\lambda_{\mathrm{p}})$$(\lambda_{M})$$(\lambda_{p}/\lambda_{M})_{cal}\approx 1.20$

実験は局在した撹乱がいくつか重なったものでわれわれの初期条件と対応するものではな
い. なお, 最大増幅率に対応する波長で完全な 6 角形パターンができたとすると, この値は 2/ 3\approx L 祐になる. このようにほぼ最大増幅率に対応する波長で決定される間隔で分布する. 最初にピークの種ができる間隔が最大増幅率を与える波長で決定されることを$\sqrt$

反映しているためと思われる. 後にピークの吸収があるのでやや間隔が広くなったのであ
ろう.

6

おわりに
非常に薄いレイリー. テイラー不安定性でどのような表面変形が生じるを潤滑近似による長波長方程式を直接数値シミュレーションおよひ解析をすることで調べた. 破断する直前は, 軸対称定常解で近似できる釣鐘状ピーク (液滴) が重ならないように空間分布する.
特に, 今回はランダムな初期条件から生じた釣鐘状ピークの空間分布について調べた.

ピークの空間分布の統計を特徴付けるために,

ドローネ分割 (ボロ\acute {分割) を用いた.
これにより隣接数およひピーク間距離を合理的に定義できるようになった. また, 一様乱数との結果を比較し, それとは一致しないことを指摘した. ピークは固有の大きさ (0 次のベツセル関数で定義できる) を持っており, 重なると融合する. このような状態を表現するために排除面積を持っランダムな点の分布を考えればよいだろう. 1 円玉が重ならないようにランダムに並べたものに相当すると考えられる. 平べったい箱の中に 1 円玉を適$\text{ノ}$当に並べて上手に横に揺するとこのような分布になるのではないかとナイーブに考えていこれを 「1 円玉モデル」 と呼ぶことにする. 数値計算で得られた単位面積当
たりピークの面積と 1 円玉モデルのそれとを一致させて, ドローネ分割による統計量を比る. そこで,較する. これを現在行っている. 発展方程式で得られた乱れたパターンがランダムといえるのか ? 何らかの相互作用のためにランダムとは異なるのか ? といった点を解明してぃき

124
他の方程式系においても空間的に乱れた定常解が見出されている [16]. この場合は単純に液滴の融合ルールのようなものはないが, 同様の考え方で調べることができないか検討中である. 特徴的な構造を持つ系において, その乱れた空間分布を構造間の相互作用にもとづいて導出することは興味深い課題と考えている.

7

参考文献

[1] A. Oron, S. H. Davis and S. G. Bankoff (1997), Rev. Mod. Phys. 69, 931.
[2] J. J. Liu, B. Schneider and J. P. Gollub (1995), Phys. Fluids 7, 55.
[3] M. Fermigier, L. Limat, J. E. Wesffreid, P. Boudinet and C. Quilliet (1992), J. Fluid
Mech. 236, 349.
[4] 0. Lioubashevski, H. Arbell and J. Fineberg (1996), Phys. Rev. Lell. 76, 3959.
[5] Y. Murakami and M. Chikano (2001), Phys. Fluids 13,65.
[6] A. Sharma and R. Khanna (1998), Phys. Rev. Lett. 85, 3463.
[7] A. Oron (2000), Phys. Fluids 12, 1633.
[8] P. S. Hammond (1983), J. Iluid Mech. 137, 363.
[9] S. G. Yiantsios and B. G. Higgins (1989), Phys. Fluids Al 1484.
[10]A. Oron and P. Rosenau (1992), J. Phys. $(f\}nnoe)2131$ .
[11] L. W. Schwartz (1999), Advances in $C\alpha ting$ and Drying of Thin Films, 3rd $Eum\mu an$.
Sy\mbox{\boldmath $\tau$}nl 真)sium,
数理研研究集会講究録, 1209, 223.
,[12] 飯田和雄, 村上洋一
[13] 村上洋–, 飯田和雄 (2001), 数理研研究集会講究録, 12$1,83.
[14] Xuehou Tan, 平田富夫 (2 1), 計算幾何学入門 (森北出版).
$Cc\omega t_{\dot{l}}ng$$\mathrm{E}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{n}\mathrm{g}\mathrm{e}\mathrm{n},105$$(2\mathrm{N}1)$[15] M. Tanemura (2001), Statistical $D$-sttibutims of Poisson Vomnoi Cells in Two
and Three Dimensions Research Memorandum No. 796 The Institute of Statistical
Mathematics.
[16] 村上洋– (2001), 形の科学会誌 16, 135.