Sci. Rep. Kanazawa Univ.
Vol. 57, pp.1–12, 2013

Three-dimensional numerical simulation of the
motion of a droplet on an inclined plane
N GUYEN T RI C ONG∗
Graduate School of Natural Science and Technology,
Kanazawa University, Kakuma, Kanazawa, 920-1192, Japan(Received July 4, 2013 and accepted in revised form August 9, 2013)
Abstract This paper introduces a mathematical model represented by a system of equations for a droplet moving on an inclined plane. A three-dimensional simulation for thisphenomenon is obtained by solving a system of equations based on both discrete Morseflow and the smoothed particle hydrodynamics method.
Keywords.
coupled model, free boundary, volume constraint, hyperbolic problem,
droplet motion

1

Introduction

In this work, we derive a model for a moving droplet on a surface. The motion of droplet
has been described by various models, in which it is usually difficult to model the dynamical contact angle [10]. In our model, we consider a coupling of film representing thesurface of the droplet and of fluid inside the film. Physical parameters of the film modelsuch as surface tensions or velocity determine implicitly the contact angle. We assume
that the fluid motion follows the Euler equations and that the film and fluid interact via
pressure force. This leaves us with the problem of deriving a model for the motion of the
film.
Surface tensions are important factors in describing the shape of droplets. For a liquid drop, there are three types of surface tensions: solid-surface tension γSG , the liquidsurface tension γLG and the solid-liquid interfacial surface tension γSL (Figure 1). In
equilibrium, the contact angle of the droplet is determined by the properties of the liquid
and the material on which the droplet is lying [16]. In particular, the relationship between
the contact angle and the surface tensions is given by Young’s equation:

γSG − γSL = γLG cos θ .
∗ Present address: Camau Teacher Training College, 155A, Nguyen Tat Thanh Street, Camau City, Viet-nam. e-mail: nctri.cpc@gmail.com

2

Three-dimensional numerical simulation of the motion of a droplet on an inclined plane

Figure 1: The setting of the model.
In our work, we focus only on the case where the dynamic contact angle θ ∈ [0, 90◦ ).
In this setting, we describe the shape of the film by a scalar function u. In addition, since
the motion of the droplet is constrained by the surface over which it is moving, there
appears an obstacle in our model equation. Moreover, as the droplet moves and changes
its shape, we assume that the volume of the fluid is preserved.
To summarize, we develop a model for a film with positive contact angle undervolume-constraint and an obstacle. The equation of motion for the film in this paper
is obtained as the stationary point of the action integral for the system. The numericalsolution of the model equation is obtained by the method of discrete Morse flow, which
has been developed by N. Kikuchi to solve parabolic problems [5]. This method was also
used to solve hyperbolic problems [2], [9], and one of the extensions of this method wasused to solve free-boundary problems [3], [6], [10]. Solving volume-preserving problems is another extension of this method [13], [12] and we remark that this method cannaturally apply to the free boundary problems with volume constraint [14], [1].
For the fluid motion, we use the smoothed particle hydrodynamics method (SPH) tosolve the Euler equations. The system of equations is then solved by the combinationof both methods. Through this coupling approach, we are able to realize the dynamic
contact angle in an implicit way, i.e., only through the consideration of the energy ofthe system. In this sense, we believe that this method can overcome the necessity ofprescribing the contact angle evolution explicitly. This necessity causes difficulties inmost of the existing models.

2

The model equations

In this section, we derive the model equation for the motion of a droplet on an inclined
plane. Our model is described by a system of equations which consists of the governing
equation for the motion of a film, together with the equation for the motion of fluid inside
the film.

2.1

The equation of motion for the film

In this section, we examine the surface of the droplet. The equilibrium shape of thedroplet is our starting point. From our assumptions that the droplet can be described by

3

Nguyen Tri Conga scalar function, we express the evolution of the film by a function:u : Ω × (0, T ) −→ R,
where (0, T ) is the time interval and Ω ⊂ R2 is the domain where the motion is considered. The boundary of ∂ Ω is assumed to be Lipschitz on which a Dirichlet condition
is prescribed and the film is assumed not to go under the plane. Moreover, the volume
enclosed under the film is assumed to be constant, that isΩuχu>0 dx = V > 0.

The surface energy of the stable droplet can thus be written as2
E = γg 1 + |∇u| χu>0 dx + γs χu>0 dx,
Ω

(2.1)

Ω

where γg = γLG and γs = γSL − γSG .
If we assume that the minimizer exists and is smooth, we obtain the following result(see [11]).
Lemma 1. Let the minimizer of (2.1) be smooth in {u > 0}. Then Young’s equation

γs = −γg cos θholds on ∂ {u > 0}.
We rewrite equation (2.1) as follows:2
E = γg 1 + |∇u| dx + (γg + γs )χu>0 dx − γg |Ω|.
Ω

(2.2)

Ω

Now suppose that the gradient of u remains small (i.e, the deformation of the film is verysmall), then by Taylor expansion we have1
1 + |∇u|2  1 + |∇u|2 .
2Thus, we can write the approximation of the surface energy (2.2) in the formγg|∇u|2 dx + R2 χu>0 dx
Ẽ =
Ω 2
Ω

(2.3)

where R2 = γg + γs . We want to know the motion of the film when an outer force f (u)with potential F(u) is applied on it. The potential energy due to the outer force f (u) isΩF(u)χu>0 dx.

4Three-dimensional numerical simulation of the motion of a droplet on an inclined plane

On the other hand, the kinetic energy of the vertical movement of the film is given by
 

σ 2
ut χu>0 dx,
Ω2

where σ is the area density of the surface. Therefore, the Lagrangian for the film can be
written as
 

γg
σ 2
ut χu>0 − |∇u|2 − R2 χε (u) + F(u)χu>0 dx,
L(u,t) =2
Ω 2
where χu>0 in equation (2.3) is replaced by a smoothing function χε ∈ C2 (R) satisfyingχε (s) =if s ≥ ε ,if s ≤ 010

and |χε (s)| ≤ C/ε for s ∈ (0, ε ). The purpose of smoothing is to avoid the presence ofdelta function in equation [11].
The action integral can be defined byJ(u) = T

L(u,t)dt,
0

and the problem is to find the stationary point of the functional J in a suitable function
space satisfying the given volume constraint. The problem is thus stated as follows.
Problem 1. Find the stationary state u of the functional T 
γg
σ 2ut χu>0 − |∇u|2 − R2 χε (u) + F(u)χu>0 dxdt
J(u) =0Ω2

2

in the function space
K = u ∈ H 1 (ΩT ); u|∂ Ω = 0, uχu>0 = V .
ΩLet u be a stationary point of J. We select an arbitrary function ϕ ∈ C0∞ ((0, T ) × Ω ∩{u > 0}), and denoteΦ(t) =

Ωϕ (t, x)dx.

uε = (u + εϕ )V.
V + εΦ

Since u is a stationary point, we havedJ(uε )= 0.
d ε ε =0Nguyen Tri Cong

5

We arrive at the relation
 T0Ω∩{u>0}− γg Δu − f (u) + σ utt + R2 χε (u) − λ ϕ dxdt = 0.

Then we can derive the film equation:

χu>0 σ utt = γg Δu + f χu>0 − R2 χε (u) + λ χu>0 ,
where

1
λ=
V


Ω

(2.4)

γg |∇u|2 − F(u)χu>0 + R2 uχε (u) + σ utt uχu>0 dx.

The derivation of this equation, in particular the presence of characteristic function in thecoefficients, follows the reasoning in [15].

2.2

The model of fluid

Now we consider the fluid motion inside the film. From the assumption, the domain of
fluid flow at time t is given as:Ω f (t) = {(x1 , x2 , z) ∈ R3 ; (x1 , x2 ) ∈ Ω, z ∈ (0, u(t, x1 , x2 ))}(2.5)

In this domain, we propose that the motion of the fluid follows the following equations:
Conservation of massDρ+ ρ ∇ · v = 0 in ∪t∈(0,T ) Ω f (t) × {t},
Dt

(2.6)

Conservation of momentumDv
1= − ∇P + g in ∪t∈(0,T ) Ω f (t) × {t},
Dt
ρ

(2.7)

where v is the velocity, g is the gravitation field vector and P is the pressure . In general
P is a function of ρ and the thermal energy, but in this paper the pressure is taken as afunction of ρ , namely
P = c2 (ρ − ρ0 ),where c, ρ0 are given.

2.3

The model of droplet motion

In this section, we consider the influence of fluid flow inside the film on the motion of
the droplet. In this case, the outer force f is the pressure force pushing the film from theinside. The pressure force per unit area is written as Pn, where1(−ux1 , −ux2 , 1)n= 
1 + |∇u|2

6

Three-dimensional numerical simulation of the motion of a droplet on an inclined planeis the unit outer normal vector of the surface. Therefore, P(x, u,t) is the net force which
is applied to the film [4]. Thus, the equation (2.4) becomes

χu>0 σ utt = γg Δu + Pχu>0 − R2 χε (u) + λ χu>0 ,
where

λ=

1
V


Ω

(2.8)

γg |∇u|2 − uP + R2 uχε (u) + σ utt uχu>0 dx.

For the fluid flow, we impose v(x, u,t) = (0, 0, 0) on the plane z = 0, and v(x, u,t) =(0, 0, ut ) on the film z = u(t, x1 , x2 ).
In summary, a model of droplet motion is given as

χu>0 σ utt = γg Δu + Pχu>0 − R2 χε (u) + λ χu>0Dρ= −ρ ∇ · v
Dt
1Dv= − ∇P + g
Dtρ
P = c2 (ρ − ρ0 )in Ω × (0, T ),(2.9)

in ∪t∈(0,T ) Ω f (t) × {t},(2.10)

in ∪t∈(0,T ) Ω f (t) × {t},(2.11)

in ∪t∈(0,T ) Ω f (t) × {t},(2.12)

v|z=0 = 0, v|z=u (x, u,t) = (0, 0, ut ).

(2.13)Taking the characteristic length L, characteristic time t∗ , density ρ0 and velocity v∗ , weobtain the following nondimensional form of the equations:

χu>0 utt = ΓΔu + ΣPχu>0 − Πχε (u) + λ χu>0Dρ= −ρ ∇ · v
Dt
1Dv= − ∇P + κ gDtρ
2
P = cs (ρ − 1)in Ω × (0, T ),(2.14)

in ∪t∈(0,T ) Ω f (t) × {t},(2.15)

in ∪t∈(0,T ) Ω f (t) × {t},(2.16)

in ∪t∈(0,T ) Ω f (t) × {t},(2.17)

v|z=0 = 0, v|z=u (x, u,t) = (0, 0, ut ),
where
Γ=

3

(2.18)

γgt∗2
R2t∗2
t∗2
c
ρ0 L
, cs = .
,
Σ
=
,
Π
=
,
κ
=
2
2
σL
σ
σL
L
v∗

The numerical method

In this section, we describe the numerical algorithm for solving equations (2.14)–(2.18).
These equations are solved by combining the discrete Morse flow method with the SPH
method. More specifically, we solve (2.15)–(2.18) by the SPH method (see [8], [7]), andwe construct an approximate solution to the equation of film motion using the discrete
Morse flow method.

7

Nguyen Tri Cong

Let us focus on the discrete Morse flow. First, we fix a large number N > 0, which
determines the time step h = T /N, and consider the approximate shapes of the film un
at time levels tn = nh, (n = 0, 1, 2, .., N). The shape u0 is given by the initial conditionu(0, x) and u1 can be approximated using u0 and the initial velocity as u1 = u0 + v0 h, here
v0 = ut (0, x). The approximate solution un on further time levels t = nh for n = 2, 3, .., N
is defined as the minimizer of the following functional
Jn (u) =Γ|u − 2un−1 + un−2 |2
χu>0 + |∇u|2 + Πχε (u) − ΣF(u)χu>0 dx
2
2h
2

Ω

in the admissible set

(3.1)




K = u ∈ H01 (Ω); uχu>0 = V .
Ω

Here F(u) denotes a primitive of P(u) with respect to u.
The existence of a minimizer un , n = 2, 3, .., N is shown in [11]. Moreover, the continuity of the corresponding minimizers is also shown and thus it is justified that the sets{un > 0} are open (see [14]). Now, we choose test functions with support in the set{un > 0} and compute the first variation of Jn . We find thatΩun − 2un−1 + un−2
ϕ + Γ∇un ∇ϕ − ΣP(un )ϕ + Πχε (un )ϕ dx =h2Ω

λn ϕ dx∀ϕ ∈ C0∞ (Ω ∩ {un > 0}),Ω∇un ∇ϕ dx = 0∀ϕ ∈ C0∞ (Ω ∩ {un ≤ 0}c ),

where
1
λn =
V



un − 2un−1 + un−2
un + Γ|∇un |2 − ΣP(un )un + Πχε (un )un dx.
h2Ω

Next, we define the approximate solutions uh and uh through interpolation of the minimizers {un }Nn=0 in time (Figure 2) as follows:uh (t, x) =u (t, x) =hhu0 (x),t =0

un (x),t ∈ ((n − 1)h, nh], n = 1, 2, . . . , N

u0 (x), t = 0
t−(n−1)hun (x) + nh−thh un−1 (x),λ (t) = λn ,t ∈ ((n − 1)h, nh], n = 1, 2, . . . , Nt ∈ ((n − 1)h, nh], n = 1, 2, . . . , N

(3.2)

(3.3)
(3.4)

8

Three-dimensional numerical simulation of the motion of a droplet on an inclined plane

Figure 2: Interpolation of minimizers
Definition 1. Functions uh and uh , determined in (3.2) and (3.3) from a sequence {un } ⊂K are called approximate solutions to (2.14), if the following conditions
 T
h

Ω

 

T
uth (t) − uth (t − h)h
ϕ + Γ∇uh ∇ϕ − ΣP(uh )ϕ + Πχε (uh )ϕ dxdt =λ ϕ dxdt,hh Ω∀ϕ ∈ C0∞ ([0, T ) × Ω ∩ {uh > 0}),uh ≡ 0 in (0, T ) × Ω \ {uh > 0},hand the initial conditions uh (0) = u0 , and uh (h) = u0 + v0 h are satisfied. Here, λ isdefined as
h

λ =


Ω

uth (t) − uth (t − h) h
u + Γ|∇uh |2 − ΣP(uh )uh + Πχε (uh )uh dx.
hIn order to get a minimizer un , n = 2, 3, .., N of functional Jn (u), we use the following
minimizing algorithm:
1. Given the initial condition u0 and v0 , set u1 = u0 + hv0 .
2. For n = 1, 2, .., N, determine un+1 as follows:(a) a1 = un(b) For k = 1, 2, .., Kn (maximum number of iterations) repeat:i. compute the gradient pk = u Jn (ak ),ii. search for the minimizer ak+1 of Jn in the direction −pk ,ak+1 , 0)iii. set ak+1 = max(ak+1 ),
iv. project ak+1 onto the volume-constraint hyperplane: ak+1 = Proj(v. if convergence criterion is fulfilled, leave the cycle.
(c) un+1 = aKn
In this algorithm, Jn (ak ) is approximated using finite element method for space discretization and minimizers are determined by the steepest descent method, combined with abisection method (step ii).

Nguyen Tri Cong

9

The whole system is solved by combining the discrete Morse flow with the SPH
method. At each time level t = nh, we have the shape of the film un , the positions and veMm M
locities of fluid particles {xmn }m=1 , {vn }m=1 , from which we can find the new shape un+1MmM
of the film, the new positions and velocities {xmn+1 }m=1 , {vn+1 }m=1 of the fluid particles
as follows:
1. Predict the shape of the film u∗ using the discrete Morse flow method without apressure force.
mm∗
2. Determine position xm
n+1 , velocity vn+1 , and pressure Pn+1 in the region below u ,using the SPH method.

3. Determine the new shape un+1 of the film, using the discrete Morse flow methodwith the pressure force.

4

Numerical results

4.1

Qualitative comparison with experiment

We use the above procedure to simulate the motion of a droplet under an inclined planewith angle α = 20◦ (Figure 3).

Figure 3: Motion of a droplet under an inclined plane (experiment).
The domain Ω = (0, 1) × (0, 0.6) is divided into 80 × 48 squares with Δx = 1/80 andeach square is divided into two triangle elements. In addition, the parameters of equation(2.14) are given as√Γ = 1, Σ = 0.5, Π = 1.65, ε = 3.2Δx, h = 4 × 10−4 , κ = 1, cs = 5
and the fluid inside the drop is represented by 1451 particles.

10

Three-dimensional numerical simulation of the motion of a droplet on an inclined plane

By observing the simulation results (Figure 4), it can be seen that the shape of droplet
oscillates and the volume is preserved within the discretization error while the droplet
moves. In addition, there are no particles moving out of the film during the motion. In
other words, the number of particles representing the fluid is controlled well by the film
during the motion. Judging from the obtained shapes of the droplet these results show aqualitative agreement with observations from the real experiments.
The above parameters were chosen so that they are close to the real values. However,a precise comparison with experimental data remains a future task.

t = 0.

t = 0.28

t = 0.56
t = 1.12
Figure 4: A droplet suspended under an inclined plane (simulation). Blue dots represent
the film, green points represent the fluid inside the film and black dots represent the
inclined plane.

4.2

Comparison with analytical result

Here we compute the stationary shape of a droplet with a given volume V using ourcoupled model and compare the result with the corresponding analytical solution. The
analytical solution is given by the minimizer u of the functionalJ(u) =

Ω1
1|∇u|2 + ρ gu2 χu>0 + γ χu>0 dx.
22Assuming radial symmetry u(x) = u(|x|) = u(r), it is possible to compute the minimizer
analytically in the form
√√I0 (r ρ g)
2γ I0 (α ρ g)
1−,
u(r) =√√
ρ g I1 (α ρ g)
I1 (α ρ g)

Nguyen Tri Cong

where α satisfies

11

√

α α

√
ρ g I0 (α ρ g)V ρg
−1 = √ ,√2 I1 (α ρ g)
2π 2γ

and I0 , I1 are modified Bessel functions of the first kind.
For the parameters ρ = 1, g = 9.8, γ = 0.5,V = 0.00687 we obtain α = 0.2073.

Figure 5: Cross-section of the analytical and numerical solutions.
The same parameters are used to calculate the steady-state solution by the coupledmodel. The initial condition is taken as a spherical cap with the volume V and the computations are continued until the change in the shape of the droplet is sufficiently small.
The cross-section of analytical and numerical solutions is plotted in Figure 5. The obtained radius of stationary droplet is αnum = 0.2209, which is slightly bigger than theanalytical value. This is probably caused by the ε -smoothing of contact angle in our numerical model. Nevertheless, it is necessary to carry out detailed convergence analysis inorder to fully validate the model.

5Conclusions

We derive a simple model for the motion of a droplet moving on the plane. The film
motion is a hyperbolic free boundary problem with volume-constraint, which was solved
using the discrete Morse flow method. The fluid motion equations, on the other hand,were solved using the SPH method. Moreover, we showed a three-dimensional simulation of a droplet moving under an inclined plane. Comparison of numerical results withexperimental data and analytical solution was performed.

Acknowledgement
The author expresses his gratitude to Professor Seiro Omata, for many ideas and advice
about the problem.

12

Three-dimensional numerical simulation of the motion of a droplet on an inclined plane

References[1] E. G INDER , K. ŠVADLENKA , A variational approach to a constrained hyperbolic free boundary
problem, Nonlinear Analysis 71, 2009, pp. e1527-e1537.
[2] K. H OSHINO , N. K IKUCHI , On a construction of weak solutions to linear hyperbolic partial differential systems with the higher integrable gradients, J. Math. Sci. 93, 1999, pp. 636-652.
[3] H. I MAI , K. K IKUCHI , K. NAKANE , S. O MATA , T. TACHIKAWA , A numerical approach to theasymptotic behavior of solutions of a one-dimensional hyperbolic free boundary problem, JJIAM 18(1), 2001, pp. 43-58.
[4] M. K AZAMA , S. O MATA , Modelling and computation of fluid-membrane interaction, Nonlinear
Analysis 71, 2009, pp. e1553-e1559.
[5] N. K IKUCHI , An approach to the construction of Morse flows for variational functionals, in ”Nematics- Mathematical and Physical Aspects”, ed.: J.-M. Coron, J.-M. Ghidaglia, F. Hélein, NATO Adv. Sci.
Inst. Ser. C: Math. Phys. Sci. 332, Kluwer Acad. Publ., Dodrecht-Boston-London, 1991, pp. 775-786.
[6] K. K IKUCHI , S. O MATA , A free boundary problem for a one dimensional hyperbolic equation, Adv.
Math. Sci. Appl. 9 (2), 1999, pp. 775-786.
[7] G.R. L IU , M.B. L IU , Smooth Particle Hydrodynamics a Meshfree Particle Method, World Scientific
Publishing, London, 2003.
[8] J.J. M ONAGHAN , Simulating free surface flows with SPH, Journal of Computational physics 110,
1994, pp. 399-406.
[9] S. O MATA , A numerical method based on the discrete Morse semiflow related to parabolic and hyperbolic equation, Nonlinear Analysis 30 (4), 2003, pp. 2181-2187.
[10] S. O MATA , M. K AZAMA , H. NAKAGAWA , Variational approach to evolutionary free boundary
problems, Nonlinear Analysis 71, 2009, pp. e1547-e1552.
[11] K. Š VADLENKA , Mathematical analysis and numerical computation of volume-constrained evolutionary problems, involving free boundaries, PhD. thesis, Kanazawa University, Japan, 2008.
[12] K. ŠVADLENKA , A new numerical model for propagation of tsunami waves, Kybernetika 43 (6), 2007,
pp. 893-902.
[13] K. ŠVADLENKA , S. O MATA , Construction of solutions to heat-type problems with volume constraint
via the discrete Morse flow, Funks. Ekvac. 50, 2007, pp. 261-285.
[14] T. YAMAZAKI , S. O MATA , K. ŠVADLENKA , K. O HARA , Construction of approximate solution to a
hyperbolic free boundary problem with volume constraint and its numerical compution, Advances inMathematical Sciences and Applications 16, 2006, pp. 57-67.
[15] H. YOSHIUCHI , S. O MATA , K. ŠVADLENKA , K. O HARA Numerical solution of film vibration withobstacle, Advances in Mathematical Sciences and Applications 16, 2006, pp. 33-43.
[16] T. YOUNG , An essay on the cohesion of fluids, Philos. Trans. R. Soc. London 95, 1805, pp. 65-87.