- Email: [email protected]
Variational approach to evolutionary free boundary problems
Nonlinear Analysis 71 (2009) e1547–e1552
Contents lists available at ScienceDirect
Nonlinear Analysis journal homepage: www.elsevier.com/locate/na
Variational approach to evolutionary free boundary problems Seiro Omata, Masaki Kazama ∗ , Hideaki Nakagawa Kanazawa University, College of Science and Engineering, Kakuma-machi, 920-1192 Kanazawa, Japan
article
info
MSC: 35L70 35R35 47J30 58E50 74F10
abstract The analysis of evolutionary nonlocal problems represented by volume–constrained motion is discussed. In the physical sense, these problems can describe the motion of droplets on a solid plane. The variational approach called the discrete Morse flow combined with a fluid-dynamics model appears to be a suitable tool. © 2009 Elsevier Ltd. All rights reserved.
Keywords: Partial differential equation of hyperbolic type Free boundary Volume constraint Numerical solution Coupled model Variational method
1. Introduction The aim of this work is to develop a simple model for the motion of droplets on surfaces, such as a water droplet sliding on a tilted plane. There exist various models describing this kind of motion. One group of models relies on the lubrication approximation of de Gennes and considers droplets overcoming shear exerted by the solid surface without changing their shape (see, e.g., [1] or [2]). Another group of models is based on Navier–Stokes equations, where the treatment of the interface so as to express surface tension is the main task (see [3]). One example may be the recently developed phasefield method adopting Navier–Stokes equations [4]. In such methods it is usually difficult to model the dynamical contact angle without prescribing it explicitly and the interface in many of them becomes smeared. Our model addresses these points in the sense that the contact angle is determined implicitly from physical parameters of the model such as surface tensions or fluid velocity, which seems to be more natural. Moreover, the surface of the drop is given exactly as a graph. The aspect of changing contact angle is important in many applications, including micro- and nano-fluidics, development of surfaces with special wetting properties, effective heat transfer, spraying of pesticides or paints, printing, etc. A method for changing the qualities of the surface using nanomachines has been reported recently, which makes it possible to regulate the motion of a drop with high accuracy. 2. Model The idea of the model is to divide the droplet into two parts: one is a film representing its surface and governing surface tension and the other is the fluid inside this film.
∗
Corresponding author. E-mail addresses: [email protected] (S. Omata), [email protected] (M. Kazama), [email protected] (H. Nakagawa). 0362-546X/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.na.2009.01.231
e1548
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
The film is used to express surface tension and contact angle and should preserve volume since it prescribes the domain where incompressible fluid moves. Moreover, the droplet moves on an obstacle so a free boundary appears. The fluid inside the film behaves according to Navier–Stokes equations. We consider only the cases when the contact angle θ is smaller than 90◦ . Then we can describe the surface as a scalar function u : (0, T ) × Ω → R, where (0, T ) is the time interval and Ω is the domain where the motion is considered. The plane, on which the drop rests, corresponds to 0-level set of the function u. The boundary of the set {u > 0} ≡ {(t , x) ∈ (0, T ) × Ω : u(t , x) > 0} is the free boundary. First, we consider the motion of the constrained film adopting for it the following equations (see [14]):
ρ utt (t , x) = γ ∆u(t , x) + λ(t ) for (t , x) ∈ (0, T ) × Ω ∩ {u > 0}, γ ρ |∇ u|2 − u2t = Q 2 on (0, T ) × Ω ∩ ∂{u > 0},
(1)
u(0, x) = u0 (x),
(3)
2
2
ut (0, x) = v0 (x) in Ω .
(2)
Here ρ, γ and Q are given positive constants. The parameter Q describes the adhesive or surface tension properties of the materials. The value of u is set to zero outside of {u > 0}. Finally, λ is a Lagrange multiplier originating in the requirement of volume preservation
Z Ω
u(t , x)dx = V
∀t ∈ [0, T ],
and is defined by
λ=
1
Z
V
ρ utt u + γ |∇ u|2 dx.
Ω
This equation as it stands can model motion of soap bubbles attached to surfaces [5]. It can be derived by the following consideration (see also [6]). First, we suppose that the potential energy for the shape of the film is described by the formula
Z γ 2
Ω
|∇ u|2 + Q 2 χu>0 dx,
where χu>0 is the characteristic function of the set {u > 0}. We define the action integral of the film as J (u) =
T
Z
Z γ 2
Ω
0
|∇ u|2 + Q 2 χu>0 −
ρ 2
u2t χu>0 dxdt .
(4)
Then we can derive the film equation (1) by calculating the first variation using volume-preserving perturbations dJ
u+δζ R 1+(δ/V ) ζ
dδ
=0
δ=0
for any ζ ∈ C0∞ ((0, T ) × Ω ∩ {u > 0}). Next, let η ∈ C0∞ ((0, T ) × Ω , R × Rm ), z = (t , x) ∈ (0, T ) × Ω , define τδ (z ) = z + δη and choose Kδ in such a way that the perturbation Kδ u(τδ−1 (z )) preserves volume. By carrying out the inner variation
dJ (Kδ u(τδ−1 (z ))) dδ
= 0, δ=0
we get the free boundary condition (2). Now we consider the coupled model including fluid motion inside the film. We treat only the two-dimensional case for the time being. In this case, the region filled with fluid is Ωu (t ) = {(x, y)| u(t , x) > 0, y ∈ (0, u(t , x))}. We introduce the following equations as the model for the interaction between the fluid and the film.
ρχu>0 utt = γ uxx − Q 2 χ0 (u) + p|y=u − σ gu + χu>0 λ in (0, T ) × Ω Z 1 with λ = ρ uutt + γ u2x + Q 2 χ0 (u)u + σ gu2 − up|y=u dx V
(5)
Ω
σ (vt + (v · ∇)v) = −∇ p + µ∆v in
[
{ t } × Ωu ( t )
t ∈(0,T )
(6)
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
∆p = −σ ∇ · [(v · ∇)v]
[
in
{ t } × Ωu ( t )
e1549
(7)
t ∈(0,T )
∂ p = µˆy · ∆v y=0 ∂ y y =0
v(x, 0, t ) = 0,
(8)
∂(v · τ) ut = 0, v · n|y=u = p ∂ n y=u 1 + u2x ∂ p σ = − p γ uxx − Q 2 χε (u) + p|y=u − σ gu + λ 2 ∂ n y=u ρ 1 + ux σ −p 2(v · xˆ )|y=u uxt + (v · xˆ )2 |y=u uxx + µn · ∆v|y=u . 2 1 + ux
(9)
(10)
We shall explain the meaning of new symbols appearing in the above identities. Function χ (u) ∈ C 2 (R) is a smoothing of χu>0 satisfying
χε (s) =
0, 1,
s≤0 ε≤s
and χε0 (s) ≤ C /ε for s ∈ (0, ε). The regularization is necessary from the viewpoint of numerical computation but is also thought to express the real physical situation in a better way than the sharp contact angle. p(x, y, t ) and v(x, y, t ) are pressure distribution and velocity distribution of the fluid, respectively.pσ denotes the density of fluid, µ is viscosity coefficient and g is a symbol for gravitational acceleration. Moreover, n = 1/ 1 + u2x (−ux , 1) is the √ unit outer normal vector to the surface u, and τ = 1/ 1 + ux (1, ux ) is its unit tangent vector, while xˆ , yˆ denote unit vectors in the direction of axes, i.e., (1, 0) and (0, 1), respectively. By Ωu we mean the time-dependent region filled with the fluid, surrounded by the obstacle and the function u. Eq. (5) is the equation for the film mentioned earlier but written as only one condition. By a formal calculation of the limit as → 0 one can derive the free boundary condition (2) and show the equivalence of the formulations. We have added into the equation a gravity term and outer force term coming from the effect of fluid flow. Particularly, the term p|y=u describes the pressure force exerted by the fluid on the film. Eqs. (6) and (7) are incompressible viscous fluid equations with boundary conditions (8)–(10). We impose no-slip condition on the rigid boundary y = 0, and free slip condition on the film boundary y = u.
3. Variational method Our next goal is to solve the system (5)–(10) numerically. Regarding the equation for the film, we use the discrete Morse flow approach introduced first in [7]. The discrete Morse flow is a variational method that solves time-dependent problems by discretizing time and defining a sequence of minimization problems approximating the original problem. The corresponding minimizers are then interpolated with respect to time and discretization parameter is sent to zero (see also [8,9]). The advantage of this method over other approaches lies mainly in its constructivity and in the absence of restrictive assumptions, such as convexity. We shall explain the idea on the example of the wave equation [10]. We consider the following problem: utt (t , x) = ∆u(t , x) u( t , x ) = 0
in QT = (0, T ) × Ω ,
(11)
on (0, T ) × ∂ Ω ,
u(0, x) = u0 (x),
(12)
ut (0, x) = v0 (x) in Ω .
(13)
First, we fix a natural number N > 0, determine the time step h = T /N and put u1 (x) = u0 (x) + hv0 (x). Function u0 corresponds to the approximate solution at time level t = 0, while function u1 is the approximate solution at time level t = h. We define the approximate solution un on further time levels t = nh for n = 2, 3, . . . , N, to be a minimizer of the following functional in H01 (Ω ): Jn (u) =
|u − 2un−1 + un−2 |2
Z Ω
2h2
dx +
1 2
Z Ω
|∇ u|2 dx.
(14)
We observe that the second term of the functional is lower-semicontinuous with respect to sequentially weak convergence in H 1 (Ω ) and the first term is continuous in L2 (Ω ). The existence of minimizers then follows immediately from the fact that the functionals are bounded from below. This is a crucial advantage over the continuous version of this functional, the Lagrangian of the type (4). As the next step, we define the approximate solutions u¯ h and uh through interpolation of the minimizers {un }Nn=0 in time:
e1550
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
u¯ h (t , x) =
u0 (x), un (x),
t =0 t ∈ ((n − 1)h, nh], n = 1, . . . , N ,
(15)
u0 (x), t = 0 t − (n − 1)h nh − t un (x) + un−1 (x), t ∈ ((n − 1)h, nh], n ≤ N . h h Since un is a minimizer of Jn , the first variation of Jn at un vanishes. Thus, for any ϕ ∈ H01 (Ω ) we have
(
u ( t , x) = h
T
Z
Z
uht (t ) − uht (t − h)
Ω
h
h
ϕ + ∇ u¯ h ∇ϕ dxdt = 0 ∀ϕ ∈ L2 (0, T ; H01 (Ω )).
(16)
Now, we take the time step to zero. To be able to do so, some estimate on the approximate solutions is needed. We obtain the following energy estimate:
kuht (t )kL2 (Ω ) + k∇ u¯ h (t )kL2 (Ω ) ≤ CE for a.e. t ∈ (0, T ),
(17)
1
where constant CE depends on H -norms of the initial data but is independent of h. Thanks to estimate (17), we can extract a subsequence converging weakly to a function u ∈ H 1 (QT ) and pass to the limit in (16) as h → 0+ to conclude that T
Z
Z Ω
0
(−ut ϕt + ∇ u∇ϕ)dxdt −
Z
v0 ϕ(0, x)dx = 0 ∀ϕ ∈ C0∞ ([0, T ) × Ω ).
Ω
Moreover, it can be shown that u satisfies boundary condition (12) and remaining initial condition (13) in the sense of traces. Therefore, we have proved by the discrete Morse flow method that there exists a weak solution u ∈ H 1 (QT ) to problem (11)–(13). The method of discrete Morse flow can be naturally applied to nonlocal problems and extends even to free-boundary problems. The advantage of our approach regarding globally constrained problems lies in the fact that a semi-discretization of time allows us to use results from elliptic theory. Moreover, the variational principle enables us to deal with the constraint by incorporating it in the set of admissible functions. In other words, we minimize a functional corresponding to (14) in the set
K = {u ∈ H01 (Ω );
Z udx = V },
(18)
Ω
instead of minimizing in H01 (Ω ). In this way, we avoid the direct treatment of the nonlocal term. Specifically, to solve problem (1)–(3), [5] suggests the following version of discrete Morse flow method: Minimize Jn (u) :=
Z Ω
ρ
|u − 2un−1 + un−2 |2 2h2
χu>0 dx +
γ 2
Z
2
Z
|∇ u| dx + Ω
Ω
Q 2 χ (u)dx
(19)
in the function set
n
KV = u ∈ H (Ω ); u = 0 on ∂ Ω , 1
Z
o
Ω
uχu>0 dx = V .
One of the important points in the construction of approximate solution of the type (15) to free boundary constrained problems is to show the continuity of the corresponding minimizers and thus justify that the sets {uh > 0} are open (see [5] for the proof in the above case). The discrete Morse flow method has been applied to the analysis of hyperbolic volume-constrained evolutionary problems. We state here a theorem from [11] dealing with a volume-constrained hyperbolic equation. Theorem 1. Let us consider the following hyperbolic problem: utt (t , x) = ∆u(t , x) + λ(u),
λ=
1 V
Z Ω
(utt u + |∇ u|2 )dx,
(20)
u(t , x) = g (t , x) on (0, T ) × ∂ Ω ,
(21)
u(0, x) = u0 (x),
(22)
ut (0, x) = v0 (x) in Ω .
Let T > 0 and Ω be a bounded domain in Rm with Lipschitz continuous boundary ∂ Ω . We assume that g ∈ L2 (∂ Ω ) but put here g = 0 for u0 , v0 belonging to H 1 (Ω ) satisfy the compatibility conditions u0 (x) = g (x), v0 (x) = 0 for R simplicity. Further R x ∈ ∂ Ω and Ω u0 dx = V , Ω v0 dx = 0. Then there is a weak solution u ∈ H 1 (0, T ; L2 (Ω )) ∩ L∞ (0, T ; H01 (Ω )) satisfying the R condition of constant volume, u(0) = u0 and the following identity for all test functions φ ∈ C0∞ ([0, T ) × Ω ) with Φ = Ω φ dx T
Z
Z
0
=
Ω
1 V
(−ut φt + ∇ u∇φ)dxdt − T Z
Z 0
Z Ω
(−ut (uΦ )t + |∇ u| Φ )dxdt − 2
Ω
v0 φ(0)dx 1 V
Z Ω
u0 v0 Φ (0)dx.
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
e1551
We cite one more result from [12], related to the case investigated here, i.e., a volume-constrained free boundary hyperbolic equation in one space dimension. Theorem 2. Let Ω ⊂ R1 be a domain and T > 0 a given final time. Assume that f (t , x, u) is continuous in the variable u and satisfies |f (t , x, u)| ≤ Cf |u|+ Γ (t , x) for some small constant Cf and a nonnegative Γ ∈ L2 ((0, T )× Ω ). Further, let u0 and v0 be functions from the space H01 (Ω ) fulfilling the compatibility conditions mentioned in the previous theorem. Then there exists a weak solution to the volume-preserving hyperbolic problem (5) with an obstacle defined as the plane {u = 0} in the following sense. A function u ∈ H 1 (0, T ; L2 (Ω )) ∩ L∞ (0, T ; H01 (Ω )) is called a weak solution if u(0, x) = u0 , if u = 0 outside {u > 0} and if ∞ ∞ the following R identity holds for all φ(t , x) ∈ C0 (([0, T ) × Ω ) ∩ {u > 0}) and an arbitrary fixed u˜ ∈ C0 (([0, T ] × Ω ) ∩ {u > 0}) satisfying Ω u˜ (t , x)dx = V : T
Z 0
=
Z Ω
1
(−ρ ut φt + γ ∇ u∇φ + f (u)φ)dxdt − Z
V
0
T
Z Ω
where Φ (t ) denotes
Z Ω
ρv0 φ(0, x)dx
(−ρ ut (˜uΦ )t + (γ ∇ u∇ u˜ + u˜ f (u))Φ )dxdt −
R
Ω
1
Z
V
Ω
ρ u˜ 0 v0 Φ (0)dx,
(23)
φ(t , x)dx and u˜ 0 = u˜ (0, x).
In the above theorem we have set f (u) = Q 2 χ0 (u) − p|y=u + σ gu. 4. Numerical results In order to solve numerically Eqs. (5)–(10), we have developed a numerical algorithm combining the discrete Morse flow method with MPS (Moving Particle Semi-implicit) method (see [13]). Particle methods are convenient for problems with moving boundaries since they do not use any mesh. The whole system is solved by alternately solving the model for the fluid and the model for the film. To be more precise, at each time level, we solve (6)–(10) by MPS method. After that, by using the obtained pressure distribution, we are able to solve (5) by discrete Morse flow method. We performed numerical computations for a droplet hanging from a tilted ceiling using this method. Eqs. (5)–(10) were solved with the following values of parameters: ρ = 1., γ = 1/2.2, Q 2 = 0.5/γ , σ = 5. and µ = 0.25. The plane was inclined at an angle π /8. The equation for the film was discretized in time by standard finite element method. We took an interval of length L = 1.5 and partitioned it into 800 elements of length Dx = L/800. The discrete Morse flow uses a time step Dx/10 and the minimization is achieved through the steepest descent method. The discretization for the fluid uses 952 particles. The initial shape of the drop is taken close to the stationary shape on a horizontal plane. Pictures show the results of numerical computation. We observed that the droplet moves while oscillating. One can notice the presence of rotating flow inside the droplet by tracing one fixed particle (the marked particle in the figures). This agrees with observations from real experiments.
5. Conclusion We have developed a model for motion of small droplets on rigid surfaces, where the surface tension and, thus, contact angle play an important role. The resulting numerical algorithm couples a minimization method for the surface of the drop
e1552
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
and a particle method for the fluid inside. Numerical results show qualitative agreement with observed facts. However, the model also needs quantitative comparison, which is our future goal. Acknowledgements The work of the first author was partly supported by the Grant-in-Aid for Scientific Research (B 18340047), Japan Society for the Promotion of Science (JSPS). The work of the second and third authors was supported by JASSO (Japan Student Services Organization). References [1] F. Brochard, Motions of droplets on solid surfaces induced by chemical or thermal gradients, Langmuir 5 (1989) 432–438. [2] R.S. Subramanian, N. Moumen, J.B. McLaughlin, Motion of a drop on a solid surface due to a wettability gradient, Langmuir 21 (2005) 11844–11849. [3] D.B. Kothe, Perspective on Eulerian finite volume methods for incompressible interfacial flows, in: Free surface flows (Udine, 1997), in: CISM Courses and Lectures, vol. 391, Springer, Vienna, 1998, pp. 267–331. [4] N. Takada, J. Matsumoto, S. Matsumoto, N. Ichikawa, Application of a phase-field method to the numerical analysis of motions of a two-phase fluid with high density ratio on a solid surface, J. Comput. Sci. Technol. 2 (2) (2008) 318–329 (special issue on Computational methods for multiphase flows). [5] T. Yamazaki, S. Omata, K. Švadlenka, K. Ohara, Construction of approximate solution to a hyperbolic free boundary problem with volume constraint and its numerical computation, Adv. Math. Sci. Appl. 16 (2006) 57–67. [6] K. Švadlenka, S. Omata, Mathematical analysis of a constrained parabolic free boundary problem describing droplet motion on a surface, Indiana Univ. Math. J. (in press). [7] N. Kikuchi, An approach to the construction of Morse flows for variational functionals, in: J.-M. Coron, J.-M. Ghidaglia, F. Hélein (Eds.), Nematics – Mathematical and Physical Aspects, in: NATO Adv. Sci. Inst. Ser. C: Math. Phys. Sci., Kluwer Acad. Publ., Dordrecht, Boston, London, 1991, pp. 195–198. [8] T. Nagasawa, S. Omata, Discrete Morse semiflows of a functional with free boundary, Adv. Math. Sci. Appl. 2 (1) (1993) 147–187. [9] S. Omata, A numerical method based on the discrete Morse semiflow related to parabolic and hyperbolic equation, Nonlinear Anal. 30 (4) (1997) 2181–2187. [10] A. Tachikawa, A variational approach to constructing weak solutions of semilinear hyperbolic systems, Adv. Math. Sci. Appl. 4 (1994) 93–103. [11] K. Švadlenka, S. Omata, Mathematical modelling of surface vibration with volume constraint and its analysis, Nonlinear Anal. 69 (2008) 3202–3212. [12] E. Ginder, K. Švadlenka, Mathematical analysis of a constrained parabolic free boundary problem describing droplet motion on a surface, 2008. Preprint. [13] S. Koshizuka, H. Tamako, Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dyn. 4 (1) (1995) 29–46. [14] K. Kikuchi, S. Omata, A free boundary problem for a one dimensional hyperbolic equation, Adv. Math. Sci. Appl. 9 (1999) 775–786.
Contents lists available at ScienceDirect
Nonlinear Analysis journal homepage: www.elsevier.com/locate/na
Variational approach to evolutionary free boundary problems Seiro Omata, Masaki Kazama ∗ , Hideaki Nakagawa Kanazawa University, College of Science and Engineering, Kakuma-machi, 920-1192 Kanazawa, Japan
article
info
MSC: 35L70 35R35 47J30 58E50 74F10
abstract The analysis of evolutionary nonlocal problems represented by volume–constrained motion is discussed. In the physical sense, these problems can describe the motion of droplets on a solid plane. The variational approach called the discrete Morse flow combined with a fluid-dynamics model appears to be a suitable tool. © 2009 Elsevier Ltd. All rights reserved.
Keywords: Partial differential equation of hyperbolic type Free boundary Volume constraint Numerical solution Coupled model Variational method
1. Introduction The aim of this work is to develop a simple model for the motion of droplets on surfaces, such as a water droplet sliding on a tilted plane. There exist various models describing this kind of motion. One group of models relies on the lubrication approximation of de Gennes and considers droplets overcoming shear exerted by the solid surface without changing their shape (see, e.g., [1] or [2]). Another group of models is based on Navier–Stokes equations, where the treatment of the interface so as to express surface tension is the main task (see [3]). One example may be the recently developed phasefield method adopting Navier–Stokes equations [4]. In such methods it is usually difficult to model the dynamical contact angle without prescribing it explicitly and the interface in many of them becomes smeared. Our model addresses these points in the sense that the contact angle is determined implicitly from physical parameters of the model such as surface tensions or fluid velocity, which seems to be more natural. Moreover, the surface of the drop is given exactly as a graph. The aspect of changing contact angle is important in many applications, including micro- and nano-fluidics, development of surfaces with special wetting properties, effective heat transfer, spraying of pesticides or paints, printing, etc. A method for changing the qualities of the surface using nanomachines has been reported recently, which makes it possible to regulate the motion of a drop with high accuracy. 2. Model The idea of the model is to divide the droplet into two parts: one is a film representing its surface and governing surface tension and the other is the fluid inside this film.
∗
Corresponding author. E-mail addresses: [email protected] (S. Omata), [email protected] (M. Kazama), [email protected] (H. Nakagawa). 0362-546X/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.na.2009.01.231
e1548
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
The film is used to express surface tension and contact angle and should preserve volume since it prescribes the domain where incompressible fluid moves. Moreover, the droplet moves on an obstacle so a free boundary appears. The fluid inside the film behaves according to Navier–Stokes equations. We consider only the cases when the contact angle θ is smaller than 90◦ . Then we can describe the surface as a scalar function u : (0, T ) × Ω → R, where (0, T ) is the time interval and Ω is the domain where the motion is considered. The plane, on which the drop rests, corresponds to 0-level set of the function u. The boundary of the set {u > 0} ≡ {(t , x) ∈ (0, T ) × Ω : u(t , x) > 0} is the free boundary. First, we consider the motion of the constrained film adopting for it the following equations (see [14]):
ρ utt (t , x) = γ ∆u(t , x) + λ(t ) for (t , x) ∈ (0, T ) × Ω ∩ {u > 0}, γ ρ |∇ u|2 − u2t = Q 2 on (0, T ) × Ω ∩ ∂{u > 0},
(1)
u(0, x) = u0 (x),
(3)
2
2
ut (0, x) = v0 (x) in Ω .
(2)
Here ρ, γ and Q are given positive constants. The parameter Q describes the adhesive or surface tension properties of the materials. The value of u is set to zero outside of {u > 0}. Finally, λ is a Lagrange multiplier originating in the requirement of volume preservation
Z Ω
u(t , x)dx = V
∀t ∈ [0, T ],
and is defined by
λ=
1
Z
V
ρ utt u + γ |∇ u|2 dx.
Ω
This equation as it stands can model motion of soap bubbles attached to surfaces [5]. It can be derived by the following consideration (see also [6]). First, we suppose that the potential energy for the shape of the film is described by the formula
Z γ 2
Ω
|∇ u|2 + Q 2 χu>0 dx,
where χu>0 is the characteristic function of the set {u > 0}. We define the action integral of the film as J (u) =
T
Z
Z γ 2
Ω
0
|∇ u|2 + Q 2 χu>0 −
ρ 2
u2t χu>0 dxdt .
(4)
Then we can derive the film equation (1) by calculating the first variation using volume-preserving perturbations dJ
u+δζ R 1+(δ/V ) ζ
dδ
=0
δ=0
for any ζ ∈ C0∞ ((0, T ) × Ω ∩ {u > 0}). Next, let η ∈ C0∞ ((0, T ) × Ω , R × Rm ), z = (t , x) ∈ (0, T ) × Ω , define τδ (z ) = z + δη and choose Kδ in such a way that the perturbation Kδ u(τδ−1 (z )) preserves volume. By carrying out the inner variation
dJ (Kδ u(τδ−1 (z ))) dδ
= 0, δ=0
we get the free boundary condition (2). Now we consider the coupled model including fluid motion inside the film. We treat only the two-dimensional case for the time being. In this case, the region filled with fluid is Ωu (t ) = {(x, y)| u(t , x) > 0, y ∈ (0, u(t , x))}. We introduce the following equations as the model for the interaction between the fluid and the film.
ρχu>0 utt = γ uxx − Q 2 χ0 (u) + p|y=u − σ gu + χu>0 λ in (0, T ) × Ω Z 1 with λ = ρ uutt + γ u2x + Q 2 χ0 (u)u + σ gu2 − up|y=u dx V
(5)
Ω
σ (vt + (v · ∇)v) = −∇ p + µ∆v in
[
{ t } × Ωu ( t )
t ∈(0,T )
(6)
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
∆p = −σ ∇ · [(v · ∇)v]
[
in
{ t } × Ωu ( t )
e1549
(7)
t ∈(0,T )
∂ p = µˆy · ∆v y=0 ∂ y y =0
v(x, 0, t ) = 0,
(8)
∂(v · τ) ut = 0, v · n|y=u = p ∂ n y=u 1 + u2x ∂ p σ = − p γ uxx − Q 2 χε (u) + p|y=u − σ gu + λ 2 ∂ n y=u ρ 1 + ux σ −p 2(v · xˆ )|y=u uxt + (v · xˆ )2 |y=u uxx + µn · ∆v|y=u . 2 1 + ux
(9)
(10)
We shall explain the meaning of new symbols appearing in the above identities. Function χ (u) ∈ C 2 (R) is a smoothing of χu>0 satisfying
χε (s) =
0, 1,
s≤0 ε≤s
and χε0 (s) ≤ C /ε for s ∈ (0, ε). The regularization is necessary from the viewpoint of numerical computation but is also thought to express the real physical situation in a better way than the sharp contact angle. p(x, y, t ) and v(x, y, t ) are pressure distribution and velocity distribution of the fluid, respectively.pσ denotes the density of fluid, µ is viscosity coefficient and g is a symbol for gravitational acceleration. Moreover, n = 1/ 1 + u2x (−ux , 1) is the √ unit outer normal vector to the surface u, and τ = 1/ 1 + ux (1, ux ) is its unit tangent vector, while xˆ , yˆ denote unit vectors in the direction of axes, i.e., (1, 0) and (0, 1), respectively. By Ωu we mean the time-dependent region filled with the fluid, surrounded by the obstacle and the function u. Eq. (5) is the equation for the film mentioned earlier but written as only one condition. By a formal calculation of the limit as → 0 one can derive the free boundary condition (2) and show the equivalence of the formulations. We have added into the equation a gravity term and outer force term coming from the effect of fluid flow. Particularly, the term p|y=u describes the pressure force exerted by the fluid on the film. Eqs. (6) and (7) are incompressible viscous fluid equations with boundary conditions (8)–(10). We impose no-slip condition on the rigid boundary y = 0, and free slip condition on the film boundary y = u.
3. Variational method Our next goal is to solve the system (5)–(10) numerically. Regarding the equation for the film, we use the discrete Morse flow approach introduced first in [7]. The discrete Morse flow is a variational method that solves time-dependent problems by discretizing time and defining a sequence of minimization problems approximating the original problem. The corresponding minimizers are then interpolated with respect to time and discretization parameter is sent to zero (see also [8,9]). The advantage of this method over other approaches lies mainly in its constructivity and in the absence of restrictive assumptions, such as convexity. We shall explain the idea on the example of the wave equation [10]. We consider the following problem: utt (t , x) = ∆u(t , x) u( t , x ) = 0
in QT = (0, T ) × Ω ,
(11)
on (0, T ) × ∂ Ω ,
u(0, x) = u0 (x),
(12)
ut (0, x) = v0 (x) in Ω .
(13)
First, we fix a natural number N > 0, determine the time step h = T /N and put u1 (x) = u0 (x) + hv0 (x). Function u0 corresponds to the approximate solution at time level t = 0, while function u1 is the approximate solution at time level t = h. We define the approximate solution un on further time levels t = nh for n = 2, 3, . . . , N, to be a minimizer of the following functional in H01 (Ω ): Jn (u) =
|u − 2un−1 + un−2 |2
Z Ω
2h2
dx +
1 2
Z Ω
|∇ u|2 dx.
(14)
We observe that the second term of the functional is lower-semicontinuous with respect to sequentially weak convergence in H 1 (Ω ) and the first term is continuous in L2 (Ω ). The existence of minimizers then follows immediately from the fact that the functionals are bounded from below. This is a crucial advantage over the continuous version of this functional, the Lagrangian of the type (4). As the next step, we define the approximate solutions u¯ h and uh through interpolation of the minimizers {un }Nn=0 in time:
e1550
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
u¯ h (t , x) =
u0 (x), un (x),
t =0 t ∈ ((n − 1)h, nh], n = 1, . . . , N ,
(15)
u0 (x), t = 0 t − (n − 1)h nh − t un (x) + un−1 (x), t ∈ ((n − 1)h, nh], n ≤ N . h h Since un is a minimizer of Jn , the first variation of Jn at un vanishes. Thus, for any ϕ ∈ H01 (Ω ) we have
(
u ( t , x) = h
T
Z
Z
uht (t ) − uht (t − h)
Ω
h
h
ϕ + ∇ u¯ h ∇ϕ dxdt = 0 ∀ϕ ∈ L2 (0, T ; H01 (Ω )).
(16)
Now, we take the time step to zero. To be able to do so, some estimate on the approximate solutions is needed. We obtain the following energy estimate:
kuht (t )kL2 (Ω ) + k∇ u¯ h (t )kL2 (Ω ) ≤ CE for a.e. t ∈ (0, T ),
(17)
1
where constant CE depends on H -norms of the initial data but is independent of h. Thanks to estimate (17), we can extract a subsequence converging weakly to a function u ∈ H 1 (QT ) and pass to the limit in (16) as h → 0+ to conclude that T
Z
Z Ω
0
(−ut ϕt + ∇ u∇ϕ)dxdt −
Z
v0 ϕ(0, x)dx = 0 ∀ϕ ∈ C0∞ ([0, T ) × Ω ).
Ω
Moreover, it can be shown that u satisfies boundary condition (12) and remaining initial condition (13) in the sense of traces. Therefore, we have proved by the discrete Morse flow method that there exists a weak solution u ∈ H 1 (QT ) to problem (11)–(13). The method of discrete Morse flow can be naturally applied to nonlocal problems and extends even to free-boundary problems. The advantage of our approach regarding globally constrained problems lies in the fact that a semi-discretization of time allows us to use results from elliptic theory. Moreover, the variational principle enables us to deal with the constraint by incorporating it in the set of admissible functions. In other words, we minimize a functional corresponding to (14) in the set
K = {u ∈ H01 (Ω );
Z udx = V },
(18)
Ω
instead of minimizing in H01 (Ω ). In this way, we avoid the direct treatment of the nonlocal term. Specifically, to solve problem (1)–(3), [5] suggests the following version of discrete Morse flow method: Minimize Jn (u) :=
Z Ω
ρ
|u − 2un−1 + un−2 |2 2h2
χu>0 dx +
γ 2
Z
2
Z
|∇ u| dx + Ω
Ω
Q 2 χ (u)dx
(19)
in the function set
n
KV = u ∈ H (Ω ); u = 0 on ∂ Ω , 1
Z
o
Ω
uχu>0 dx = V .
One of the important points in the construction of approximate solution of the type (15) to free boundary constrained problems is to show the continuity of the corresponding minimizers and thus justify that the sets {uh > 0} are open (see [5] for the proof in the above case). The discrete Morse flow method has been applied to the analysis of hyperbolic volume-constrained evolutionary problems. We state here a theorem from [11] dealing with a volume-constrained hyperbolic equation. Theorem 1. Let us consider the following hyperbolic problem: utt (t , x) = ∆u(t , x) + λ(u),
λ=
1 V
Z Ω
(utt u + |∇ u|2 )dx,
(20)
u(t , x) = g (t , x) on (0, T ) × ∂ Ω ,
(21)
u(0, x) = u0 (x),
(22)
ut (0, x) = v0 (x) in Ω .
Let T > 0 and Ω be a bounded domain in Rm with Lipschitz continuous boundary ∂ Ω . We assume that g ∈ L2 (∂ Ω ) but put here g = 0 for u0 , v0 belonging to H 1 (Ω ) satisfy the compatibility conditions u0 (x) = g (x), v0 (x) = 0 for R simplicity. Further R x ∈ ∂ Ω and Ω u0 dx = V , Ω v0 dx = 0. Then there is a weak solution u ∈ H 1 (0, T ; L2 (Ω )) ∩ L∞ (0, T ; H01 (Ω )) satisfying the R condition of constant volume, u(0) = u0 and the following identity for all test functions φ ∈ C0∞ ([0, T ) × Ω ) with Φ = Ω φ dx T
Z
Z
0
=
Ω
1 V
(−ut φt + ∇ u∇φ)dxdt − T Z
Z 0
Z Ω
(−ut (uΦ )t + |∇ u| Φ )dxdt − 2
Ω
v0 φ(0)dx 1 V
Z Ω
u0 v0 Φ (0)dx.
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
e1551
We cite one more result from [12], related to the case investigated here, i.e., a volume-constrained free boundary hyperbolic equation in one space dimension. Theorem 2. Let Ω ⊂ R1 be a domain and T > 0 a given final time. Assume that f (t , x, u) is continuous in the variable u and satisfies |f (t , x, u)| ≤ Cf |u|+ Γ (t , x) for some small constant Cf and a nonnegative Γ ∈ L2 ((0, T )× Ω ). Further, let u0 and v0 be functions from the space H01 (Ω ) fulfilling the compatibility conditions mentioned in the previous theorem. Then there exists a weak solution to the volume-preserving hyperbolic problem (5) with an obstacle defined as the plane {u = 0} in the following sense. A function u ∈ H 1 (0, T ; L2 (Ω )) ∩ L∞ (0, T ; H01 (Ω )) is called a weak solution if u(0, x) = u0 , if u = 0 outside {u > 0} and if ∞ ∞ the following R identity holds for all φ(t , x) ∈ C0 (([0, T ) × Ω ) ∩ {u > 0}) and an arbitrary fixed u˜ ∈ C0 (([0, T ] × Ω ) ∩ {u > 0}) satisfying Ω u˜ (t , x)dx = V : T
Z 0
=
Z Ω
1
(−ρ ut φt + γ ∇ u∇φ + f (u)φ)dxdt − Z
V
0
T
Z Ω
where Φ (t ) denotes
Z Ω
ρv0 φ(0, x)dx
(−ρ ut (˜uΦ )t + (γ ∇ u∇ u˜ + u˜ f (u))Φ )dxdt −
R
Ω
1
Z
V
Ω
ρ u˜ 0 v0 Φ (0)dx,
(23)
φ(t , x)dx and u˜ 0 = u˜ (0, x).
In the above theorem we have set f (u) = Q 2 χ0 (u) − p|y=u + σ gu. 4. Numerical results In order to solve numerically Eqs. (5)–(10), we have developed a numerical algorithm combining the discrete Morse flow method with MPS (Moving Particle Semi-implicit) method (see [13]). Particle methods are convenient for problems with moving boundaries since they do not use any mesh. The whole system is solved by alternately solving the model for the fluid and the model for the film. To be more precise, at each time level, we solve (6)–(10) by MPS method. After that, by using the obtained pressure distribution, we are able to solve (5) by discrete Morse flow method. We performed numerical computations for a droplet hanging from a tilted ceiling using this method. Eqs. (5)–(10) were solved with the following values of parameters: ρ = 1., γ = 1/2.2, Q 2 = 0.5/γ , σ = 5. and µ = 0.25. The plane was inclined at an angle π /8. The equation for the film was discretized in time by standard finite element method. We took an interval of length L = 1.5 and partitioned it into 800 elements of length Dx = L/800. The discrete Morse flow uses a time step Dx/10 and the minimization is achieved through the steepest descent method. The discretization for the fluid uses 952 particles. The initial shape of the drop is taken close to the stationary shape on a horizontal plane. Pictures show the results of numerical computation. We observed that the droplet moves while oscillating. One can notice the presence of rotating flow inside the droplet by tracing one fixed particle (the marked particle in the figures). This agrees with observations from real experiments.
5. Conclusion We have developed a model for motion of small droplets on rigid surfaces, where the surface tension and, thus, contact angle play an important role. The resulting numerical algorithm couples a minimization method for the surface of the drop
e1552
S. Omata et al. / Nonlinear Analysis 71 (2009) e1547–e1552
and a particle method for the fluid inside. Numerical results show qualitative agreement with observed facts. However, the model also needs quantitative comparison, which is our future goal. Acknowledgements The work of the first author was partly supported by the Grant-in-Aid for Scientific Research (B 18340047), Japan Society for the Promotion of Science (JSPS). The work of the second and third authors was supported by JASSO (Japan Student Services Organization). References [1] F. Brochard, Motions of droplets on solid surfaces induced by chemical or thermal gradients, Langmuir 5 (1989) 432–438. [2] R.S. Subramanian, N. Moumen, J.B. McLaughlin, Motion of a drop on a solid surface due to a wettability gradient, Langmuir 21 (2005) 11844–11849. [3] D.B. Kothe, Perspective on Eulerian finite volume methods for incompressible interfacial flows, in: Free surface flows (Udine, 1997), in: CISM Courses and Lectures, vol. 391, Springer, Vienna, 1998, pp. 267–331. [4] N. Takada, J. Matsumoto, S. Matsumoto, N. Ichikawa, Application of a phase-field method to the numerical analysis of motions of a two-phase fluid with high density ratio on a solid surface, J. Comput. Sci. Technol. 2 (2) (2008) 318–329 (special issue on Computational methods for multiphase flows). [5] T. Yamazaki, S. Omata, K. Švadlenka, K. Ohara, Construction of approximate solution to a hyperbolic free boundary problem with volume constraint and its numerical computation, Adv. Math. Sci. Appl. 16 (2006) 57–67. [6] K. Švadlenka, S. Omata, Mathematical analysis of a constrained parabolic free boundary problem describing droplet motion on a surface, Indiana Univ. Math. J. (in press). [7] N. Kikuchi, An approach to the construction of Morse flows for variational functionals, in: J.-M. Coron, J.-M. Ghidaglia, F. Hélein (Eds.), Nematics – Mathematical and Physical Aspects, in: NATO Adv. Sci. Inst. Ser. C: Math. Phys. Sci., Kluwer Acad. Publ., Dordrecht, Boston, London, 1991, pp. 195–198. [8] T. Nagasawa, S. Omata, Discrete Morse semiflows of a functional with free boundary, Adv. Math. Sci. Appl. 2 (1) (1993) 147–187. [9] S. Omata, A numerical method based on the discrete Morse semiflow related to parabolic and hyperbolic equation, Nonlinear Anal. 30 (4) (1997) 2181–2187. [10] A. Tachikawa, A variational approach to constructing weak solutions of semilinear hyperbolic systems, Adv. Math. Sci. Appl. 4 (1994) 93–103. [11] K. Švadlenka, S. Omata, Mathematical modelling of surface vibration with volume constraint and its analysis, Nonlinear Anal. 69 (2008) 3202–3212. [12] E. Ginder, K. Švadlenka, Mathematical analysis of a constrained parabolic free boundary problem describing droplet motion on a surface, 2008. Preprint. [13] S. Koshizuka, H. Tamako, Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dyn. 4 (1) (1995) 29–46. [14] K. Kikuchi, S. Omata, A free boundary problem for a one dimensional hyperbolic equation, Adv. Math. Sci. Appl. 9 (1999) 775–786.