数理解析研究所講究録
1522 巻 2006 年 111-119

111

水面に拘束されたバブルの運動の数値シミュレーション

Numerical simulation of motion of a bubble
restrained to water surface
自然科学研究科山崎崇史 (YAMAZAKI, TAKASHI)
Graduate School of Natural Science and Technology,
Kanazawa University
金沢大学 自然科学研究科小俣正朗 (OMATA, SEIRO)
Graduate School of Natural Science and Technology,
Kanazawa University
自然科学研究科長山雅晴 (NAGAYAMA, MASAHARU)
Graduate School of Natural Science and Technology,
Kanazawa University

金沢大学

金沢大学

Abstract. The aim of this paper is to develop the numerical method for solving. The
the hyperbolic free-boundary problem under the volume conservation
equation treated in this paper has a degenerate hyperbolic operator and non-local
term (Lagrange multiplier) coming from volume constraint. The typical phenomena
described by this equation are the motion of bubble on water surface and oil droplet
motion on a plane. We try to realize them numcricaly. It is the method called discrete
Morse flow $(DMF)$ of hyperbolic type that plays an important role in the numericalsimulation or the stage of construction of an approximate weak solution. The DMF isbased on the variational method using time-semidiscretized functional and is suitable
to treat a volume constraint condition.
$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\mathrm{r}\mathrm{a}\dot{\mathrm{i}}\mathrm{t}$1Introduction

In this paper, we treat numerically the motion of a bubble which moves on a water
surface or an oil droplet on a flat plane. There are some common features in these
phenomena. The first common feature is that the volume does not change along theirmotion. The second common feature of these phenomena is that the area where thebubble or oil droplet exist is localized. We call the set of points where the bubble
touches the water surface or where the surface of oil droplet touches the planeboundary. In this section we explain the model equation that rules the motion of the
bubble. However, the model of the bubble can be applied also to the motion of the oil
droplet because of these two common properties.
We use the graph of a scalar function to describe the shape of the bubble. The
zero level set of coincides with the water surface. The set of points where the bubble
touchcs the water surface is called free-boundary. In order to simplify the model, we
assume that the water layer is so thin that its flow or movement does not influence
the bubble. Moreover, areal density of film is expected to be constant, and the stress$fi\epsilon e-$

$u$

$u$

112

tensor density of the bubble and water surface to be homogeneous and isotropic.
We also assume that the volume of air which is surrounded by the bubble is pre, that is, the bubble movement can beserved at any time, i.e.,
described by wave equation with volume constraint. The following equation describes$\int_{\Omega}udx=M=\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$the phenomena wcll:$(x\in\Omega, 0<t<\tau)$$\chi_{u>0}u_{tt}=\Delta u-R^{2}\chi_{\epsilon}’(u)+\lambda\chi_{u>0}$

Here$\Omega$

is a bounded domain in
is used), andis a smoothing of$\mathrm{R}^{m}$$\Omega\cross(0,\tau)\subset \mathrm{R}^{m+1}$

and

.

(1.1)

(in the following, to simplify the notation,is the characteristic function of the set $\{u>0\}$$\Omega_{\tau}=$$\chi_{u>0}$$\chi_{\epsilon}\in C^{2}(\mathrm{R})$$\chi$

satisfying,$s\leq 0$

1,$\epsilon\leq s$$0$$\chi_{\epsilon}(s)=\{$. The termwith interpolating in $0<s<\epsilon$ in such a way that I
describes the net adhesive force. It is due to this rcstriction force that the generationof the new surface or movement of free-boundary are obstructed. The specificity ofequation (1.1) lies in the coefficienton the left-hand side. Bccause of this coefficient, non-negativity of the solution is guaranteed. Function , which appears in
the last term of (1.1) is a Lagrange multiplier originating in the volume preservation
condition. The explicit form of the Lagrange multiplier A is obtained as follows:$\chi_{\epsilon}’(s)|\leq C/\epsilon$$R^{2}\chi_{\epsilon}’(u)$$\chi_{u>0}$$\lambda$(1.2)$\lambda=\frac{1}{M}\int_{\Omega}[uu_{\ell t}\chi_{u>0}+|\nabla u|^{2}+R^{2}u\chi_{e}’(u)]dx$

by formal calculation with volume constraint under the assumption that A does notdepend on space variables. The integral representation of Lagrange multiplier makes
the problem more difficult. However, we can calculate an approximate solution to (1.1)
by use of a time-semidiscretized functional which is called the discre Morse flow ofhyperbolic type (see [7]).
$te$2Minimizing method for the bubble problem

Like in [8], we introduce another approximation problem to (1.1). Here, we givethe volume constraint in the admissible space for finding a minimizer of a discretizedfunctional corresponding to the Lagrangian.

Problem 2.1. Let be a bounded domain inof the follouring functional:$\Omega$$\mathrm{R}^{m}$. For $n=2,3,$$\ldots$, find minimizer$u_{n}$$J_{n}(u):= \int_{\Omega}\frac{|u-2u_{n-1}+u_{n-2}|^{2}}{2h^{2}}dx+\frac{1}{2}\int_{\Omega}|\nabla u|^{2}dx+\int_{\Omega}R^{2}\chi_{e}(u)dx$

in the function set$\mathcal{K}_{M}:=\{u\in W^{1,2}(\Omega, \mathrm{R});u=0$

on$\partial\Omega,$$\int_{\Omega}u\chi_{u>0}dx=1\}$

.

,

(2.1)

113

is to be dewith $u_{1}=u_{0}+hv_{0}$ are given and the sequencermined inductively. Moreover, by use of these minimizers, construct an approximate
weak solution to (1.1).

Functions$u_{0},$$\{u_{n}\}$$u_{1}\in \mathcal{K}_{M}$$te$Let us use test function
and take the first variation of$(u+ \delta\zeta)/\int_{\Omega}(u+\delta\zeta)\chi_{u+\delta\zeta>0}dx,$$\zeta\in C_{0}^{\infty}(\Omega\cap\{u>0\})$$\sqrt:n$$\int_{\Omega}(\frac{u-2u_{n-1}+u_{n-2}}{h^{2}}\zeta+\nabla u\nabla\zeta+R^{2}\chi_{\epsilon}’(u)\zeta)dx=\int_{\Omega}(\lambda_{n}dx$$\forall_{\zeta}\in C_{0}^{\infty}(\Omega\cap\{u>0\})$$\int_{\Omega}\nabla u\nabla\zeta dx=0$

(2.2)

,

$\forall_{\zeta}\in C_{0}^{\infty}(\Omega\cap\{u\leq 0\}^{\mathrm{o}})$

.

(2.3)

Here,
$\mathrm{A}_{n}=\int_{\Omega}(\frac{u-2u_{m-1}+u_{m-2}}{h^{2}}u+|\nabla u|^{2}+R^{2}\chi_{e}’(u)u)dx$

is the Lagrange multiplier coming from the volume constraint. From the second identity, we can conclude thatoutside the set $\{u>0\}$ .
$u\equiv 0$Interpolation in time and approximate solution

3and introduce
In this section, we carry out interpolation in time of minimizers
the approximate weak solution. First we state the definition of weak solution.
$\{u_{n}\}$Definition 3.1. We call$u$a weak solution to (1.1), if$u$satisfies the following:$\int_{0}^{\tau}\int_{\Omega}(-\mathrm{u}_{t}\zeta_{t}+\nabla u\nabla\zeta+R^{2}\chi_{\epsilon}’(u)\zeta)dxdt-\int_{\Omega}v_{0}\zeta(x, \mathrm{O})dx=\int_{0}^{\tau}\int_{\Omega}\lambda\zeta dxdt$$\forall(\in C_{0}^{\infty}(\Omega \mathrm{x}[0,\tau)\cap\{u>0\})$$\mathrm{u}\equiv 0$and$u(\mathrm{O})=u_{0}$

(3.1)

,

outside of $\{u>0\}$(3.2)

in the sense of traces.

Now, we consider the approximate solutions. We define$\overline{u}^{h}$and$u^{h}$

on$\Omega \mathrm{x}(0, \infty)$by$\overline{u}^{h}(x, t)=u_{n}(x)$,$u^{h}(x, t)= \frac{t-(n-1)h}{h}u_{n}(x)+\frac{nh-t}{h}u_{n-1}(x)$$\overline{\lambda}^{h}(t)=\lambda_{n}$,,for $(x, t)\in\Omega\cross((n-1)h,nh],$ $n\in N$ . We can construct the approximate weak solution
to the bubble problem in terms of and .
$\overline{u}^{h}$$u^{h}$

114

Definition 3.2 (Approximate solution). Let
defined as above. If the following conditions$\{u_{n}\}\subset \mathcal{K}_{M}$and let$\overline{u}^{h}$and$\int_{h}^{\tau}\int_{\Omega}(\frac{u_{t}^{h}(t)-u_{t}^{h}(t-h)}{h}\zeta+\nabla\overline{u}^{h}\nabla\zeta+R^{2}\chi_{e}’(u^{h})()dxdt=\int_{h}^{\tau}\int_{\Omega}\overline{\lambda}^{h}\zeta dx$$\forall_{\zeta\in C_{0}^{\infty}}(\Omega\cross[0, \tau)\cap\{u^{h}>0\})$in$u^{h}\equiv 0$$u^{h}$

be,,$\Omega \mathrm{x}(h,\tau)\backslash \{u^{h}>0\}$

and the initial conditions $u^{h}(0)=u_{0},$ $u^{h}(h)=u_{0}+hv_{0}$ are satisfied, then we callapproximate solutions to the bubble problem.
and$\overline{u}^{h}$$u^{h}$

If one can pass to the limit as, then the above approximate solutions are
expected to converge to the solution of $(3.1)-(3.2)$ . We expect that a good regularity ofminimizers
should imply that the limit of agrees with the Lagrange multiplierof (1.1). By now, we eould not get any result concerning the convergence of approximate
solutions. However, we can still carry out numerical computations using a minimizing
method.
$harrow \mathrm{O}$$\overline{\lambda}^{h}$$\{u_{n}\}$4$\lambda$

Numerical method

Here we present the numerical method and experimental results. We apply a finite
element method with minimizing algorithm and find minimizer of the approximate
functional
defined above via steepest descent method for a fixed time step .
The time step and diameter of each finite element are chosen small enough related tothe approximation parameter 6. However, the functional (2.1) does not correspond tothe case when two or more bubbles exist. In that case, we should detect each bubble
and find minimizer for each bubble. After detection of each bubble, the basic courseof our minimizing algorithm is the following:$J_{n}(u)$

$n$

$h$

1.

$u_{0},$$u_{1}(=u_{0}+hv_{0})\in \mathcal{K}_{M}$2. for $n=1,2,$.

$\mathrm{a}$$\mathrm{b}$$\ldots,$$N$

is an initial value, determine$u_{n+1}$

by the following procedure$v^{1}=u_{n}$. for $k=1,2,$($\mathrm{i}$$K_{n}$.

$\mathrm{i}\mathrm{i}$.

,is made so large that a certain error estimator becomes small enough)as retrieval direction
search the minimizer
in the(steepest descentmethod)$\ldots$$K_{n}$$p^{k}=\nabla_{u}J_{n}(v^{k})$$\hat{v}^{*+1}$$\mathrm{d}\mathrm{i}\mathrm{r}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}-p^{k}$$\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{h}\mathrm{o}\mathrm{d}+\mathrm{b}\mathrm{i}\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$$\mathrm{i}\mathrm{i}\mathrm{i}.\hat{v}^{k+1}=\max\{v^{k+1},0\sim\}$$\mathrm{i}\mathrm{v}$. projection onto the volume-constraint plane:$v^{k+1}=P(\hat{v}^{k+1})$$\mathrm{c}$.

$u_{\mathrm{n}+1}=v^{K_{\hslash}+1}$

115

By executing tasks a.
, we get an approximate minimizer for a fixed time step .
, i.e.
In the following simulations, we use equation with a damping term$n$$-\mathrm{c}.$$\gamma u_{t}$$\chi_{u>0}u_{u}+\gamma u_{t}=\Delta u-R^{2}\chi_{\epsilon}’(u)-\lambda\chi_{u>0}$

We choose the parameters as follows: $h=5\cross 10^{-3},$

4.1$\epsilon=0.1,$.

$\gamma=0.5$.

Neumann condition

The first example uses Neumann boundary condition and $R^{2}=0.35$ . An initial
velocity is imparted to the bubble. In this case, it approaches the boundary of .
After touching the boundary, the bubble moves along the boundary. The more thebubble leans against the boundary, the smaller the area of the surface of the bubble
becomes. The bubble stops and keeps the smallest surface when reaching the corner of.
$\Omega$$\Omega$(a)$t=0.\mathrm{O}$(d) $t=2.65$(b) $t=1.75$(c) $t=2.2$(e) $t=6.35$(f) $t=8.0$

Figure 1: Neumann boundary condition. The bubble which has an initial velocity
chosen in such away that it moves diagonally, moves toward the boundary. Aftertouching the boundary, the bubble moves along it and finds the corncr of the domain.

116

4.2

Collision of two bubbles

The second example treats a collision of two bubbles with the same volume. Afterthe collision, the bubbles merge and keep the smallest surface near the center of initial
position of two bubbles.

$(\mathrm{a})t=0.\mathrm{O}$$(\mathrm{b})t=0.8$$(\mathrm{c})t=0.9$(d) $t=1.1$(e) $t=1.2$(f) $t=1.5$

Figure 2: collision of two bubbles. Two bubbles separated initially merge into eachother.

4.3

Division of oil droplet

In the third example, the motion of oil drople on the flat plane on which there is aagent is considered.
nonuniform distribution of
The contact angle between the oil droplet and the plane is prescribed by the affinityof the material of the water surface to the oil droplet surface. Therefore, if there isa nonuniform distribution of surface-active agent, we can observe that the oil droplet
moves on the plane driven by the force of non-uniformity of the contact angle. The
free-boundary condition$\mathrm{s}\mathrm{u}\mathrm{r}\mathrm{f}\mathrm{a}\mathrm{c}\triangleright \mathrm{a}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}$$|\nabla u|^{2}-(u_{t})^{2}=R^{2}$on$\partial\{v>0\}$

,

(4.1)

which is obtained by passing to the limit asin (1.1) tells us that we can rulethe contact angle of the oil droplet and the plane by defining the coefficient . The
oil droplet lies down if
distribution becomes small and the droplet moves gradually$\mathit{6}arrow 0$$R$$R^{2}$

117(a) $t=0.0$(b) $t=0.85$(c) $t=1.0$(d) $t=1.15$(e) $t=1.2$(f) $t=1.5$

Figure 3: Division of bubble (case 2).
aiming at the part of$\Omega$with smaller$R^{2}$. We set

0.35
3.5 otherwise.

$(x_{1}x_{2}(x_{1}^{2}-x_{2}^{2})\leq 0)$$R^{2}=\{$

4.4

Combination with reaction-diffusion equation

In [11] it is reported that an oil droplet on the glass plane immersed in surface-active
agent solution immediately moves according to the gradient of the concentration of thesurface-active agent. We are trying to realize the motion of the oil droplet numericallyby applying the model of the bubble (1.1) to this phenomenon and combiming it withreaction-diffusion equation which shows the time development of the concentration of
surface-aetive agent. Let be the concentration of surface-active agent, then our model
equation of is as follows:
$v$

$v$

$R^{2}=( \frac{Q}{T})^{2}+\frac{\gamma_{w}(v)-\gamma_{b}}{T}$,$v_{t}=\kappa\Delta v+\mu(v_{0}-v)v\chi_{u\leq 0}-\nu v\chi_{u>0}$$\gamma_{w}=\frac{\gamma_{0}}{1+\alpha v^{\beta}}$

Here
and

$T$

$\gamma_{b}$,.

is the tension of the oil droplet surface,is the tension of glass plane surface
expresses the tension of the region where both the oil droplet and glass plane$\gamma_{w}$118surface have come in contact, $Q>0$ is adhesion between the oil and the surfaceof plane. In the numerical experiment we set a $=1,$ $\beta=2$ and impose Neumannboundary condition. We choose the uniform distribution as the initial condition for .
We use the explicit method to solve this equation. It is observed that the oil droplet
put near the boundary starts moving aiming at the area where the concentration of
the surface-active agent is high. However, the droplet slows down and stops before
reaching the center of the domain, although it should keep moving as the experimentshows. We think that it is nec\’esary to improve the model, e.g. by considering theinternal energy of the oil in addition to the surface energy.
$v$$(\mathrm{a})t=0.5$$(\mathrm{b})t=1.5$$(\mathrm{c})t=2.5$(d) $t=3.0$(e) $t=4.0$(f) $t=4.5$

Figure 4: Combination with the reaction-diffusion system.

5

Conclusions

We model the motion of a bubble on a water surface by a free-boundary problem with volume conservation constraint and we have established a numerical method
for its solution. The bubble moves on the water surfaee while changing its shape
but preserving its volume. The model equation becomes free-boundary problem ofdegenerate hyperbolic type which is difficult to treat. We have introduced the time.
semidiscretized variational method to solve this problem and it gives good numericalresults. This model can also be applied to the motion of oil on the bottom of water or to
problems related to the phenomenon of a water-drop dripping from ceiling. Therefore,

119

this work has many applications and is significant for the future studying of hyperbolic
free-boundary problems.

References[1] H. W. Alt- L. A. Caffarelli, “Existence andfor a minimum prvblem
with free boundary“, J. Reine Angew. Math., 325 (1981), 105-144.
$regula7^{\cdot}ity$[2] H. Berestycki L. A. Cafafrelli L. Nirenberg, “Uniform estimates for regularization of ffoe boundary problems“, in Analysis and Partial Differential Equation,
Marcel Dekker, New York, 1990.
--[3] M. Giaquinta, “Multiple integmb in the calculus of variations and nonlinear elliptic systems”, Princeton University Press, 1983.
[4] H. Imai K. Kikuchi K. Nakane S. Omata- T. Tachxikawa,
nume ricalapproach to the asymptotic behavior of solutions of a one-dimensional hyperbolic
ffee $bounda\eta$ prvblem“, JJIAM 18 No.1(2001), 43-58.
--$‘ {}^{t}A$-[5] K. Kikuchi-S. Omata, $u_{Affee}$ boundary problem for a one dimensional $h_{\Psi}erbolic$equation”, Adv. Math. Sci. Appl. 9 No.2 (1999), 775-786.

[6] O. Ladyzhenskaya- N. Uraltseva, ttLinear and Quasilinear Elliptic Equations”,
Academic Press, New York and London, 1968.

[7] T. Nagasawa-K. Nakane- S. Omata, “NumeComputations for motion ofvortioes govemed by a Hyperbolic Ginzburg-Landau System”, Nonlinear Anal. 51(2002) No.1 Ser A: Theory Methods, 67-77.
$r\dot{\tau}cal$[8] S. Omata,Numerical treatment of film motion utth ffee boundary“, Adv. Math.
Sci. Appl., 14, (2004), 129-137.
$‘ {}^{t}A$[9] T. Yamazaki S. Omata- H. Yoshiuchi K. Ohara,
motion on water
surface”, Gakuto Intern. Ser. Math. Sci. Appl. 23 (2005), 209-216.
--$‘ {}^{t}Bubble$[10] T. Yamazaki-S. Omata-K. Svadlcnka-K. Ohara, “Constru ction of approximate solution to a hyperbolic ffee boundary prvblem utth volume constraint and its
numerical computation“, Intem. Ser. Math. Sci. Appl. 23 (2005), 209-216. Adv.
Math. Sci. Appl. 16 No.1 (2006), 57-67.

[11] Y. Sumino-N. Magome-T. Hamada-K. Yoshikawa, “Self-Running
Emergence of Regular Motion from Nonequilibrium Noise”, Phys. Rev. Lett., 94,
068301, (2005).
$\mathrm{D}\mathrm{r}\mathrm{o}\mathrm{p}\mathrm{l}\mathrm{e}\mathrm{t}+$