数理解析研究所講究録1081 巻 1999 年 84-99

84

泡の流れのシミュレーション

名古屋大学大学院工学研究科 奥薗透

1(Tohru$\mathrm{o}\mathrm{k}\mathrm{u}\mathrm{Z}\circ \mathrm{n}\circ$)

はじめに多数の泡が集まってできた柔らかい物質をフォームと呼ぶ。 ビールの泡やシェービングフォームを想像してもらえれば良い。フォームを形成するひとつひとつの泡はセル状の形をなし、それらを隔てる液体二重膜はネットワーク状のパターンを形成する。 このようなパターンはセルパターンと呼ばれ、 ある程度乱れた (多分散の) 構造をとることが多い。
フォームはそのまま放置すればゆっくりとパターンの粗大化が起こり、セルの平均サイズは時間とともに大きくなる。 このようなパターン成長を主に律速しているのは、 セル内部の気体の (膜を通しての) 拡散過程である。
シェービングクリームを手にとると液体のようには流れず形を保持していることからわかるように、 フォームは小さな変形に対しては弾性体のように振舞う。 -方 . フォームは容易に変形可能であり、少し大きな力を加え続けるとあたかも流体のように “流れる” このような力学的性質は塑性と呼ばれ、 その流動現象は塑性流動と呼ばれる。通常、 このような塑性流動に特徴
的は時間スケールは上述の拡散過程の時間スケールに比べて非常に短く、液$(\mathrm{f}_{\circ \mathrm{a}}\mathrm{m})$$\text{、}$。

体膜中の粘性流体の運動が重要なエネルギーの散逸過程となる。.

以上のような物理過程はセルパターンという複雑な幾何学的構造の下で起こる。 したがって、パターン全体の時間発展、統計量、 あるいはマクロな応答を知ることは非常に複雑な多体問題となる。我々はこれに対し簡単な描像の下にモデル化を行い、計算機シミュレーションを行う。そうすることによって、今まで見えなかったことが明らかになり、現象への理解を深めることができる。

2モデル

本稿では簡単のため 2 次元のモデルに限定する。界面がネットワ一ク状につながったパターン (セルパターン) の時間発展を記述するための変数の選び方

はいろいろ考えられる。我々は、非常に単純に、 ネットワ一クの頂点の座標とそれらのつながりの情報のみに着目する。即ち、系の時間発展は頂点の運85
動方程式と頂点の衝突過程に対応するトポロジカルプロセスとで記述される。 ひとつの頂点は常に 3 つの界面によって共有される点であると仮定すれば (界面エネルギーにより 4 つ以上の界面によって共有される点は不安定である)、 トポロジカルプロセスは図 1 に示すような二種類の素過程 (Tl および T2 プロセスと呼ばれる) からなる。実際のシミュレーションでは、 ある辺の長さがある –定の小さな距離以下になるとそのまわりのネットワークのつながりの状況によって Tl あるいは T2 プロセスが起きる。
では、頂点の運動方程式はどのように決まるのか。以下でそれを考察する。
$\text{、}$.

2.1

成長過程の動力学 [1]

今、 滑らかな界面を考え、 界面の局所的な運動が曲率駆動型であるとする。と即ち、界面の位置 での法線方向の界面の速度を v(a)、平均曲率をすると、$\kappa(a)$

$a$

(1)

$v(a)=L\mathcal{K}(a)$$\mathrm{c}$

こで、 $L$ は運動学係数 (正の定数) である。 この式は形式的に次のように書き換えることができる。

$-$$\iota$.

(2)$\frac{\delta \mathcal{R}}{\delta v(a)}+\frac{\delta \mathcal{H}}{\delta q(a)}=0^{\cdot}\mathfrak{k}$$\mathcal{R}\overline{=}\frac{\sigma}{2L}\int \mathrm{d}av^{2}(a)$(3)$\mathcal{H}\equiv\sigma\int \mathrm{d}a$は界面張力、 $q(a)$ は位置 での界面の法線方向への変位である。
また、は散逸関数、は界面の自由エネルギーである。 したがって、 (2)は摩擦力とポテンシャルによる力の釣合の式である。
セルパターンの場合に頂点の運動方程式を得るため、頂点と頂点を結ぶと界面を直線で近似する。 この場合ダイナミカルな変数は頂点の位置
速度である。 ここで、 $i=1,2,$
これらの変数.
を用いて、頂点の運動方程式は (2) (3) に対応して、
ここで、

$a$

$\sigma$$\mathcal{R}$$\mathcal{H}$$\{r_{i}\}$$\mathrm{r}\cdot\cdot \text{は頂点_{の番}号^{を}表す。}$$\{v_{i}\}$(4)$\frac{\partial R}{\partial v_{i}}+\frac{\partial H}{\partial r_{i}}=0$$R \equiv\frac{1}{2}\sum_{i,j}v_{i}\cdot Dij$

. $vj$$H \equiv\sigma\sum_{\rangle\langle ij}|r_{ij}|$(5)と書ける。 ここで、はすべての辺 \langle ゆについて和をとることを意味する。 また、 テンソルの具体的な表式については文献 [1]を参照していただきたい。 (4) は頂点に働く力の釣合の式であると考えることができる。
$r_{ij}\equiv r_{i}-r_{j}-\backslash \text{、}\Sigma_{\langle ij\rangle}$$D_{ij}$86

2.2

流動現象のモデル [2, 3, 4]

気体相 (分散相) の体積分率が 1 に近いフォームを考える。 このような系の流動現象においては、気体を非粘性、非圧縮とすれば、液体二重膜中の流体の運動による粘性散逸が重要なエネルギー散逸過程となる。特に Plateau borderと呼ばれる頂点近傍の領域での散逸が重要であり、 この領域が二重膜に沿って速さ $U$ で動くときのエネルギー散逸率はキャピラリー数が小さい場合には、に比例することが流体力学的計算からわかる [5]。
液体二重膜中の粘性散逸は頂点の運動に対して余分の摩擦力を生じる。 したがって、頂点の運動方程式は次のように変更されるであろう。
$U^{5/3}$(6)$f_{i}^{D}+f^{P}i=- \sigma\sum_{j}(i)rij/|r_{i}j|$は (4) の第 1 項を表し、は新たに加わった粘性散逸による摩擦力である。和は頂点 とつながっている 3 つの頂点 についてとる。
第 1 節でも述べたように、上で議論した流体運動の時間スケールは、拡散によってパターンが粗大化する時間スケールに比べて非常に短い。 したがって、泡の流動現象に興味があるならば、 その時間スケールではパターンの粗大化は起こらないとしてもよいであろう。そこで我々は、パターンの粗大化に関する摩擦力を考慮するかわりに、 パターンの粗大化が起こらない、
即ち各セルの面積が時間とともに変化しないという拘束を課す。 したがって、頂点の運動方程式は次のように書ける。
ここで、$f_{i}^{D}$$\Sigma_{j}(i)$$f_{i}^{P}$

$i$

$j$

$f_{i}^{D}$$f_{i}^{P}+ \sum_{\alpha}\lambda\frac{\partial}{\partial r_{i}}A\alpha=-\sigma\sum_{j}(i)rij/|r_{i}j|$

ここで、 A。はセル$\alpha$の面積、$\lambda_{\alpha}$(7)は Lagrange の未定乗数で、 条件(8)$\underline{\mathrm{d}}A_{\alpha}=0$$\mathrm{d}t$から決定される。
次に、の表式を得るため、 各々の Plateau border での粘性散逸が独を立に起こると仮定し、系全体での粘性散逸率$f_{i}^{P}$

$Q$

$Q$

$=$

$\sum_{i}\sum_{j}Q_{ij}(i)$$Q_{ij}$$\equiv$$v_{i}^{(j)}\cdot f_{ij}^{P}$

(9)

(10)

87と書く。 ここで、

$f_{ij}^{P}$

は頂点$i$に働く辺$\langle ij.

\rangle$.

$:\text{に沿う摩擦力_{で}あ^{る}}$ 。

また、

(11)

$v_{i}^{(j)}\equiv v_{i^{-}}u_{ij}$はは頂点 から遠く離れた辺 \langle 切上の流体の速度である。
であり、旧 i 力が伸び縮みするときのエネルギー散逸率であり、それに関係した摩擦である。流体力学的な計算の結果より力が$Q_{ij}$$i$$u_{ij}$$\langle$$f_{ij}^{P}$(12)$Q_{ij}.\propto|v_{i}\hat{r}_{i}j|(j).5/3$であるから、$\hat{r}_{ij}$方向の摩擦力はその方向の速度の 2/3 乗に比例する:(13)$f_{ij}^{P}\propto|v_{i}\hat{r}_{i}(j).|2/3j\hat{r}_{ij}$ここで、$\hat{r}_{ij}\equiv r_{ij}/|r_{ij}|_{0}$はまだ決定されておらず任意性が残っている。 これを解決.
するため次の 2 つの方法を考える。
ところで、$u_{ij}$$\omega$

方法 I

平均流 $u(r)$ (例えば shear flow 等) を与え、
$u_{ij}=u(r_{i})$(14)$Q=. \sum_{i}[v_{i}-u(r.i)]\cdot fiP$(15)とする。 このとき、 (10) はとなり、$f_{i}^{P}= \sum_{j}\backslash \prime f_{ij}^{P}$(16)の速度の平均流からのずれに対する摩擦力である。

は、頂点

$i$

方法 II

対称性 $u_{ij}=u_{ji\text{、}}f_{ij}^{P}=-f_{j}^{P}i.\text{を仮定すると、}$ $Q$ は$=$$\sum_{i}\sum_{j}(i)(vi-u_{ij})\cdot f^{P}ij$

$Q-$

$=$

$\sum_{i}v_{i}\cdot f_{i}^{P}$(17)

88と書ける。 したがって、頂点
だしはおよび$i$に働く摩擦力は (16) の形で与えられる。 たの関数で上の対称性を満たすようなものである。上の対称性に加え、頂点の運動方程式 (7) がガリレイ変換に対して不変である$f_{ij}^{P}$$v_{i}$$v_{j}$

ことを要請すれば(18)$u_{ij}=-\wedge(v_{i}2+v_{j})$が得られる。

3

シミ$\dot{I}$レーションめ結果

2.1 節のモデルによるシミュレーションから得られたパターンを初期状態とし、 Lees-Edwards の周期境界条件によって単純ずり流を与え、 シミュレ一ションを行った。以下では前節の方法 II によるシミュレーションの結果を列挙する。 シミュレーションの詳細は文献 [3] を参照されたい。
図2

初期状態 (a) および変形中 (b) の系のパターンである。

図 3 および 4

単位面積あたりの界面エネルギー $E(t)$ およびずり応力の時間変化を示す。横軸は何れもずり速度 の逆数でスケールした時間即ちずり歪みである。微小な変形領域では弾性的振舞いをしているが、大変形領域では非常に特異な振舞いをしていることがわかる。

$\tau_{xy}(t)$$\dot{\gamma}$$\ovalbox{\tt\small REJECT}$$(\dot{\gamma}t>1)$

図 5 図 4 の–部門拡大したものである。図中に描かれている縦棒はトポロ
ジカルプロセス (T1) が起こった時刻を示している。 Tl プロセスが起こることによって応力が解放されていることがわかる。 ここには示されていないがエネルギー $E(t)$ についてもほば同様の振舞いをする。
!

図 6 セルの流れの様子を図示したものである。右側の棒グラフにはセルの速度の成分が表示してある。微小変形領域ではきれいな shear flow が形成されているが大変形領域でエネルギーを放出しているときには乱れ$x$$(\mathrm{a})_{\text{、}}$た速度場が形成されている (b)。大変形領域においてもエネルギーを蓄積しているときには速度場は乱れていない (c)。

このように大変形領域においては、乱れた速度場を伴ったエネルギーの解放過程が Tl プロセズによって誘起され、このような過程は間欠的に繰り返される。我々はこのエネルギ一解放の過程を泡の系における “ なだれ (avalanche)”現象と名付けた。

89をそのなだれのサイズ1 回のなだれで解放された界面エネルギーとしたとき、 なだれのサイズ分布 $P(s)$ を両対数プロットしたものである。
$P(s)$ は巾分布していることがわかる。
図7$s$

図8

時系列 $E(t)$ に対するパワースペクトルの両対数プロットである。 これも巾則に従っていることがわかる。

図 9 Tl プロセスが起こつた場所を示したものである。 この図では 1 回のrun で起こったすべての Tl プロセスの位置を重ね書きしている。 これを見る限り、 ほぼ–様な場所で起こっている。 しかしながら、 いくつかのなだれ. について起こった Tl プロセスを取り出してみると図 10 のようになる。同–のなだれで起こった Tl プロセスは回じシンボルで示してある。 これを見ると、 1 回のなだれで起こる Tl プロセスの位置には相関があることがわかる。

Acknowledgments . 本研究は中部大学川崎恭治教授との共同研究です。

References$,[1]$T. Okuzono and K. Kawasaki, Trends$\mathrm{i}.\mathrm{n}$Stat. Phys. 1, 65 (.1994).

[2] T. Okuzono, K. Kawasaki, and T. Nagai, J. Rheol. NY 37, 571 (1993).
[3] T. Okuzono and K. Kawasaki, Phys. Rev.

$\mathrm{E}\bm{5}1$,1246 (1995).

[4] K. Kawasaki and T. Okuzono, in Dynamical Systems and Chaos, Vol.
2: Physics,. Y. Aizawa, S. Saito and K. Shiraiwa (World Scientific,Singapore, 1995)$\mathrm{e}\mathrm{d}\mathrm{s}$[5] L. W. Schwartz and H. M. Princen, J. Colloid Interface Sci. 118. 201(1987).

90

Tl
(a)

$\mathrm{T}2$(b)$)$

閃

1

91

(a)

(b)
$|^{\iota}\underline{\neg\neq}\mathrm{Z}$92.
$\gamma$$1^{\iota}\underline{\neq}\neg 3$

$t$

93

$\vee\wedgearrow$ト.
$\gamma$$l^{\backslash }\underline{\backslash \neq}\neg 4$

$t$

94

$\bigwedge_{\vee}rightarrow$ト.
$\gamma$$t$$\dot{|}^{\mathrm{t}}\mathit{4}\supset 5$

95$\overline{\underline{\mathrm{o}}}$可プ$\Phi$

粛–,..

鵡

–

$\overline{*0}$–96$\wedge\eta$$\vee \mathrm{h}$$\underline{\mathrm{O}\frac{\mathrm{o}}{\mathfrak{W}}}$$-\mathrm{O}$$-\delta$$-\angle\pm$$\log_{10^{S}}$$)$

コ\neq

7

$–\angle$

97$\vee\wedge 3$$\Phi^{\mathrm{H}}$ロ$\overline{\infty}$$\underline{\mathrm{O}}$$\log_{10}\omega$

図$8’$

98

Locations of Tl$>$X
図$\uparrow$$\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{c}.\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{e}\mathrm{s}$99

Locations of Tl during an avalanche,$\approx$

X
図$|O$