Thin-film free boundary problems for partial wetting

Journal of Computational Physics 295 (2015) 770–778

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Thin-film free boundary problems for partial wetting Dirk Peschka Weierstrass Institute, Mohrenstr. 39, 10117 Berlin, Germany

a r t i c l e

i n f o

Article history: Received 8 October 2014 Received in revised form 23 April 2015 Accepted 24 April 2015 Available online 4 May 2015 Keywords: Thin films Free boundary problem Numerical methods

a b s t r a c t We present a novel framework to solve thin-film equations with an explicit non-zero contact angle, where the support of the solution is treated as an unknown. The algorithm uses a finite element method based on a gradient formulation of the thin-film equations coupled to an arbitrary Lagrangian–Eulerian method for the motion of the support. Features of this algorithm are its simplicity and robustness. We apply this algorithm in 1D and 2D to problems with surface tension, contact angles, and gravity. © 2015 Elsevier Inc. All rights reserved.

1. Introduction and model statement In this paper we describe the creeping motion of a viscous liquid over a flat solid substrate. If we assume that at each time t the liquid occupies the volume

(t ) = {(x, y , z) ∈ R3 : 0 < z < h(t , x, y )}, then the interface between the liquid and the surrounding air is located at z = h(t , x, y ) > 0 and the solid substrate is at z = 0. The shape of the liquid domain can be entirely parametrized by the function h(t , x, y ). In the thin-film limit of the Stokes equation, the evolution of the domain (t ) is given by solutions of

∂t h = ∇ · (m(h)∇ π ),

(1)

with given initial data h(0, x, y ) = h0 (x, y ) ≥ 0. The mobility m is a nonnegative, monotone function with m(0) = 0 and typically m(h) = hα with 0 < α ≤ 3 depending on the type of friction at the substrate [1]. The most studied case is certainly the no-slip boundary condition with α = 3, see [2] and references therein. To avoid a singular energy dissipation at a moving contact line [3,4] a slip condition with α = 2 was proposed. The case α = 1 is relevant for liquid motion on porous rough surfaces or in Hele-Shaw cells [5,6]. The evolution is driven by an energy



E (h) =

1 2

|∇ h|2 dx d y + V (h) ,

(2)

R2

which defines π = δ E /δh in (1). What makes this problem sometimes difficult to treat analytically and numerically is the fact, that the function h may vanish in certain regions and thereby the support of h defined as





ω(t ) = (x, y ) ∈ R2 : h(t , x, y ) > 0 ,

(3)

depends on time and (1) becomes a free boundary problem. Depending on the choice of V (h) there are different strategies to tackle this difficulty.

E-mail address: [email protected]. http://dx.doi.org/10.1016/j.jcp.2015.04.041 0021-9991/© 2015 Elsevier Inc. All rights reserved.

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

771

For V ≡ 0 we are in the situation where the liquid spreads over the substrate with zero contact angle. Algorithms were constructed which use a globally defined h(t , ·) and enjoy the crucial property of maintaining positivity or non-negativity [7–9]. This is based on the observation in [10] that the thin-film equation has, in addition to its energy dissipation



d D = − dt E (h) =

m(∇ π )2 dx d y ≥ 0,





also a nonlinear entropy dissipation d/dt G (h) = − (h)2 . In addition to an entropy-based discretization it might be necessary to regularize the mobility to ensure non-negativity. This class of algorithms is certainly computationally efficient, since the domain is fixed and solutions are defined globally in space. Thereby people can apply highly developed numerical methods and obtain beautiful numerical results. For a review concerning the mathematics of the contact line problem we recommend the review [11]. In order to treat finite contact angles, intermolecular potentials such as

 

V (h) =

A h2

−

B



h8

dx d y ,

(4)

R2

were introduced, e.g., [12,13] or [2] and references therein. For numerical results in 2D with heterogeneous substrates we refer to [14]. Appropriate A , B lead to a so-called precursor film, i.e., regions where 0 < h ∼ ε  1 with strictly positive h, and even ensure global positivity of h with contact angles in a certain sense. Compared to the case V ≡ 0 these methods are still quite easy to implement but usually require spatial adaptivity near the approximate contact line. This makes the algorithms potentially sluggish in higher dimensions. The big advantage of these methods is that there is no problem with changes in topology, i.e., when connectivity properties of ω change. The main reason for the success of these methods might be that one actually never has to deal with a free boundary problem at all. We pursue a different strategy. The intent of this paper is to consider the case of non-smooth potentials such as



V (h) = S

χω dx d y ,

(5)

R2

which measures the size of the support (wetted region) using the characteristic function



χω of the support ω

1 (x, y ) ∈ ω, χω = 0 otherwise, with the negative spreading coefficient S ≥ 0. This is interesting because for S = 0 we cover the case V ≡ 0 and for S > 0 we use a V which is a limit of the intermolecular potential in (4) as was shown in [15]. For the case V ≡ 0 short-time existence and uniqueness of classical solutions was proven in [16] also in the context of a moving support, which is in contrast to the conventional approach using the concept of weak solutions from [10]. Existence of weak solutions for partial wetting was investigated in [17,18] and later existence and uniqueness of classical solutions was shown in [19]. The contact line singularity for an explicit contact line has been investigated in [20] using formal asymptotics and in [21] rigorously. This paper is in the spirit of the aforementioned analytical and asymptotical results. We solve (1) as a free boundary problem, where the support ω in (3) is one of the unknowns. Our strategy is to solve a variational formulation to obtain ∂t h and then consistently decompose the time-derivative into convective time-derivative and into the velocity of the corresponding comoving frame. Since we have explicit control over the support it is the authors opinion that positivity is an issue of lesser importance compared to algorithms for globally (in space) defined solutions. Nevertheless, all the techniques that have been developed for global solutions apply here the same. However, even with such an enhanced method, the inherent disadvantage of our method remains that it cannot handle topological changes of the support. Still we want to convince the reader that, with a little knowledge of finite elements, such a formulation can be superior to the classical approaches. A similar but technically more involved derivation of a gradient formulation for bilayer films was given in [22]. From a numerical standpoint our algorithm might appear similar to the one in [23] but without surfactants. The main difference is that we evaluate local conditions in a gradient-type formulation. So all boundary conditions follow naturally from the variational formulation. In contrast to [23] this locality allows us to generalize the algorithm to higher dimensions. 2. Numerical algorithm 2.1. Weak formulation We assume h is a globally continuous function with continuous first spatial derivatives inside ω , and the support of h moves with finite speed. Then the main idea of the algorithm is to solve (1) only inside ω with additional boundary conditions and kinematic conditions at the boundary γ = ∂ ω . The kinematic condition is implied by the fact that h(t , ·) ≡ 0 on γ . In order to solve (1) on ω we construct a variation formulation of the thin-film equation. Therefore we integrate (1) by parts and use the fact that the flux j = −m∇ π is zero on γ . Then we get for all test functions φ that



ω



∂t h φ + m(h)∇ π · ∇φ dx d y = 0.

(6a)

772

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

π in this weak formulation, let us formally compute the derivative of  1 2 |∇ h| + S χω dx d y .

In order to see how to define

 

E (h) =

2 R2

Using Reynolds’ transport theorem with h˙ = ∂t h and x˙ · n the normal component of the velocity of



d dt

∇ h · ∇ h˙ dx d y +

E= ω

 

γ



∇ h · ∇ h˙ dx d y +

= ω



γ we get



S + 12 |∇ h|2 x˙ · n d

 |∇ h|−1 S + 12 |∇ h|2 h˙ d .

γ

The last line is obtained by using the outer normal n = −∇ h/|∇ h| and then writing the normal velocity of the boundary ˙ This is a simple result of h ≡ 0 at the boundary γ , which by taking a time-derivative becomes x˙ · n as a function of h.

h˙ + x˙ · ∇ h = 0,

on γ .

(6b)

π = δ E /δh in a weak sense on ω means π , ϕ = diffE h [ϕ ] and thereby    πϕ dx d y = ∇ h · ∇ ϕ dx d y + |∇ h|−1 S + 12 |∇ h|2 ϕ d .

Finally, setting



ω

ω

(6c)

γ

We seek (h˙ , π ) by solving (6a) together with (6c) as our weak formulation. We discuss (6b) separately in the next section. All boundary conditions in (6a), (6c) are natural boundary conditions, time derivatives h˙ are in the Eulerian reference frame. Let us examine the weak formulation from an energetic point of view. Using the flux as above we have h˙ + ∇ · j = 0 and then one can write the energy dissipated by the viscous flow as

 D=



m(h) ∇ π

2

 dx d y =

ω

ω

| j |2 dx d y . m(h)

Another way to write the flux is j = h v¯ with v¯ = −h−1 m∇ π being the average velocity of the fluid. Near a moving contact line the average velocity v¯ would have to approach a nonzero constant vector. Then near the contact line the integrand in the dissipation is approximately | j |2 /m ≈ h−α (h v¯ )2 = | v¯ |2 h2−α as (x, y ) → (¯x, y¯ ) ∈ γ . Using the assumption that the moving contact line (for simplicity in 1D) has a small Young contact angle θ , we get a dissipation

x¯ +ε ...

2 −α | j |2 D≈ | v¯ |2 tan(θ) |x − x¯ |2−α dx + dx, x¯

x¯ +ε

m

which for α = 3 has a logarithmic singularity in the first integral and thereby the dissipation becomes infinite. As long as d the driving force dt E is finite this implies that for α ≥ 3 the contact line velocity v¯ needs to vanish, i.e., the contact line does not move. One approach to “cure” this problem is introducing a precursor using (4), which provides a cut-off for the integration. For the sharp-interface model presented in this paper, however, the most obvious cure is to modify the mobility and use mixed mobilities of the form m(h) = h3 + bh2 . The parameter b is called the slip-length [3,4,24]. Note, for mobilities with α ≥ 3 the discretization length x provides a cut-off for the integration. This leads to a contact line velocity converging to zero either logarithmically like 1/| log x| for α = 3 or algebraically for α > 3. For the interpretation of E recall that the surface energy of a liquid film is

E˜ =

 



γlg 1 + |∇ h|2 + γls − γsg dx d y , ω

where γlg , γls , γsg are the surface tensions of liquid–gas, liquid–solid, solid–gas surfaces respectively. In situations where the slope |∇ h| is small (by definition it is zero outside ω ) we get by Taylor expansion

E˜ ≈





γlg 1 + ω

1 |∇ h|2 2





1 |∇ h|2 2

+ γls − γsg dx d y = γlg

+



γlg +γls −γsg γlg



χω dx d y ,

R2

so that we get E ≈ γlg−1 E˜ if we identify S = (γlg + γls − γsg )/γlg . Note that up-to the normalization γlg−1 the constant S is the negative spreading coefficient. The sign S > 0 corresponds to partial wetting, whereas S ≤ 0 corresponds to complete

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

773

wetting. To see what the weak formulation with different signs of S implies formally, we assume h is smooth enough and integrate the first integrand on the right side of (6c) by parts and get





πϕ dx d y = − ω



 |∇ h|−1 S − 12 |∇ h|2 ϕ d ,

ϕ h dx d y + ω

γ

√

which implies π = −h in ω and a contact angle |∇ h| = 2S on γ . These are also the expressions for π and the boundary condition, that one would use in (1) from a PDE viewpoint. For S < 0 the construction above does not produce a meaningful equilibrium condition at γ . A mathematical problem that occurs for negative S is that the functional E is not lower-semicontinuous and thereby it is open if classical variational methods can even give existence of solutions for the corresponding gradient system. For the reasons given above in this paper we concentrate on the situations where S > 0 and we have either 0 < α < 3 or use a mixed mobility with a positive slip-length. 2.2. Kinematic condition Now assume h(t , x, y ) is already known and its time-derivative h˙ was computed by solving (6a), (6c). What remains is to determine h(t + τ , x, y ) and ω(t + τ ). Setting h(t + τ , ·) = h(t , ·) + τ h˙ makes no sense since h(t + τ , ·) and h(t , ·) are defined on different domains. Instead we propose an arbitrary Lagrangian–Eulerian (ALE) method based on (6b) and define a ψ : ω(t ) → R2 mapping points from ω(t ) to ω(t + τ ). The normal component of the velocity ψ˙ ≡ x˙ is fixed by (6b), but one is free to choose its tangential component and the interpolation of the boundary values to the interior. When the tangential component is fixed, then one can use (6b) to determine ψ˙ uniquely, e.g., fixing ψ˙ · t = 0 produces ψ˙ · n = h˙ /|∇ h|. One might ˙ x). then think of ψ˙ generating an infinitesimal transformation ψ(x) = x + τ ψ( If one moves the domain according to ψ , then one also needs to transform the Eulerian time-derivative into the consistent ALE time-derivative. In ALE coordinates we have H (t , x) = h(t , ψ(x)), then

˙ = h˙ + ψ˙ · ∇ h, H

(7)

follows for all x ∈ ω(t ). For simplicity assume τ  1 so that to leading order for spatial derivatives we have ∇ψ ≈ I2 , ∇ H ≈ ∇ h. The considerations before imply H˙ = 0 at x± (t ) so that H has constant support. If h˙ is known, then one can use ˙ H˙ ) is one-to-one. ∇ h to decompose it into ψ˙ and H˙ . Once the interpolation is defined, this decomposition h˙ → (ψ, To simplify things, we outline the construction in the 1D case where h(t , x) = h(t , x, y ) for all y and give a more detailed description of the method in Appendix A. The case is treated separately, because we can define the interpolation explicitly. Assume we have ω(t ) = (x− (t ), x+ (t )) being a single interval. Using (6b) we solve for x˙ ± and get x˙ ± = −h˙ /∂x h at x± , which in turn allows to define ψ˙ by linear interpolation as

˙ x) = (˙x+ − x˙ − )ξ(x) + x˙ − , ψ(

(8)

where ξ(x) = (x − x− )/(x+ − x− ). As explained, we also need to transform the time-derivatives accordingly. We get





˙ (x) = h˙ (x) + (˙x+ − x˙ − )ξ(x) + x˙ − ∂x h(x). H Now consider the general case for 2-dimensional h(t , x, y ). While the 1-dimensional case admits a generic choice of the interpolation, an optimal choice in the higher-dimensional situation depends much on the behavior of the solution that one is interested in. Let us first consider an interpolation, where ψ˙ has no tangential component on γ . One possible choice is

˙ x, y ) = α (x, y )∇ h(t , x, y ), ψ( where α (x, y ) is a yet-to-be-determined function. Using (6b) we get the boundary data α = −h˙ /|∇ h|2 on γ . Since ψ˙ is uniquely defined by α and we have its boundary data, we can, for example, interpolate ψ˙ by defining α as harmonic −α = 0 with the Dirichlet boundary condition above. In general, the choice of the interpolation and the tangential direction of ψ˙ is ambiguous and driven by computational simplicity and robustness. In one example of the next section we show, how the tangential part of ψ˙ can be modified to obtain a more robust deformation. 2.3. Algorithm summarized Let us summarize the algorithm, given the solution h(t , ·) at an arbitrary time t with support 1. First compute the time derivative h˙ ∈ V by solving the weak formulation



ω





h˙ φ + m(h)∇ π · ∇φ dx d y = 0,

 ω





πϕ − τ ∇ h˙ · ∇ ϕ dx d y =



 ∇ h · ∇ ϕ dx d y +

ω

γ

 |∇ h|−1 S + 12 |∇ h|2 ϕ d ,

ω(t ).

774

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

Fig. 1. (Left) Surface energy driven droplet evolving towards a stationary state at t = 10−5 , t = 10−3 , t = 10 as indicated by the arrow. (Right) Surface energy and normal gravity driving a droplet towards a stationary state at t = (0, 10−3 , 0.01, 0.1, 1)/2 as indicated by the arrow.

for all ϕ , φ ∈ V from a suitable discrete function space V . Note this will produce the proper contact angle, since the boundary term from the derivative δ E /δh is included here. Compared to the definition of π in (6c), we replaced ∇ h · ∇ ϕ in the second equation by ∇(h + τ h˙ ) · ∇ ϕ to obtain a stabilized semi-implicit time-discretization. 2. Based on the knowledge of the normal part of the velocity x˙ and after a possible modification of the tangential part one defines an interpolation ψ˙ to the inside. In one dimension this can be done using (8) as outlined before. The choice for the tangential part and for the interpolation depends very much on the problem at hand. For instance, if one wants to find traveling wave solutions (in higher dimensions) then one could choose the tangential part such that ψ˙ might be constant in space. 3. When the interpolation ψ˙ is chosen, the time-derivative of h and the one of H in the ALE frame transform according to (7). Finally, one can move the vertices of the finite element domain decomposition by x → x + τ ψ˙ and the nodal ˙ and obtains the new solution at time t + τ . This corresponds to setting ψ(x) = x + τ ψ( ˙ x) and values by h → h + τ H ˙ (x). The implementation details are explained in Appendix A. H (t + τ , x) = H (t , x) + τ H The principles used in the first step of the algorithm are rather standard, but to our knowledge were not applied to a free-boundary problem of thin-film type. In particular it was not used in combination with our decomposition of the ˙ in the corresponding comoving frame. In the work time-derivative into the transport velocity ψ˙ and the time derivative H by Karapetsas et al. [25] the velocity/position of the contact-line is determined from the global mass conservation constraint. However, since this is only one scalar constraint such an approach only works in one spatial dimension with a single contact point. Karapetsas et al. argued that the evaluation of the usual kinematic conditions x˙ = (m/h)∇h is not feasible, since it contains the product of the singular term ∇h and the degenerate term m/h as one approaches the boundary. It is certainly unclear how to design a robust algorithm which takes this condition literally. In this sense our approach has the advantage that we only work with h˙ , ∇ h, both of which are neither degenerate nor singular at the contact line. In addition we use local conditions, which make it possible to extend the algorithm to higher dimensions. To underline this claim we present some examples in 1D and 2D. 3. Examples in 1D Here we consider an extension of the energy considered before, namely



V (h) =





S + f (h; x) dx d y ,

(9)

ω

where f may explicitly depend on x in a parametric way. Thereby we can include effects such as partial wetting and gravity in x (tangential) and z (normal) direction. Without the tangential component we get stationary solutions, whereas with a tangential component of gravity and sufficiently small volumes we get stable traveling wave solutions. The numerical solution of the algorithm above with m(h) = h2 , f ≡ 0, S = 1, L = 1 with initial data h0 (x) = L /2 −|x − L /2| is shown in the left panel of Fig. 1. Note that we choose the mobility with this exponent to circumvent the certainly interesting discussion on the contact √ line singularity, which would prevent any motion of the contact √ line. The energy above implies the contact angle is |∂x√ h| = 2 at x± and the stationary solution is h(x) = a(1 − η2 )/ 2 where η = (x − 1/2)/a, x ∈ (1/2 − a, 1/2 + a), and a2 = 18/16. In the left panel of Fig. 1 the exact solution is shown using a few points overlaying the numerical solution with full lines. Already at t = 10 there is no visible change in the solution and the comparison with the exact stationary solution is excellent. The simulation shown here was performed using n = 100 elements and a nonuniform time-step size. Even though the initial data provided here is not continuously differentiable at x = 1/2, the numerical algorithm had no problems finding the stationary state on coarse or on fine meshes with n = 100 or n = 1000 elements and big time-step sizes such as τ = 103 . In this sense the algorithm might be also suited just to find stationary solutions of this type of PDE problems. Some details of the implementation can be found in Appendix A. The author encourages study of the MATLAB code, which can be downloaded from [26]. Gravity in tangential direction can be included in the problem by using f (h; x) = − g 1 xh, where g 1 encodes the strength of gravitational forces. Here we use g 1 = 1. Then the total energy of the system can be decreased by moving the droplet in +x direction. The corresponding PDE solution is shown in Fig. 2 for the same initial data as before but

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

775

Fig. 2. Convergence to a traveling wave solution under the action of gravitational forces in tangential direction for small volumes L = 1 at t = (0, 10−3 , 0.1, 1, 2, 3, 4) (left) and bigger volumes L = 5 and t = (0, 0.1, 1, 2, 3, 4) (right).

Fig. 3. Evolution towards stationary droplet at t = 0, 3, 12 (left) to (right) where height is indicated by the shading.

with L = 1 (left panel) and L = 5 (right panel). For large times h converges to a traveling wave solution of the form h(t , x) = F (x − v¯ t ) for some compactly supported F . The traveling wave solution F for the bigger droplet volume is more asymmetric. In the third simple example we add normal gravitational forces, which can be achieved by using for instance f (h; x) = g 2 h2 , where g 2 encodes the strength of normal gravitation forces. Here we use g 2 = 102 . Then the total energy of the system can be decreased by making the droplet slightly flat, as one can see by comparing the left and the right panel in Fig. 1. The expressions for the gravitational energy follow from the well-known (specific) gravitational energy of a point-mass e grav (x, y ) = ρ gz , where in the tilted coordinate system z = −x sin α + y cos α . Then the contribution of gravity to the energy of the thin film is

 h

 e grav (x, y ) d y dx =

ω 0

with g 1 = ρ g sin α , g 2 =





g 2 h2 − g 1 hx dx,

ω 1 2

ρ g cos α .

4. Examples in 2D Similar to the 1D example we use the mobility m(h) = h2 and the extended potential (9) with

f (h; x) = − g 1 hx + g 2 h2 , where we use f again to implement gravitation in the tangential and normal direction using g 1 and g 2 respectively. In the first example in Fig. 3 we use an initial flower shaped domain, where the boundary in polar coordinates is at r (φ) = 5(1 + 12 sin 3φ). The initial height solves −h0 = 1 with homogeneous Dirichlet boundary conditions. Here we use S = 2, g 1 = g 2 = 0 which gives rise to a contact angle |∇ h| = 2 at γ . As expected, the evolution depicted in Fig. 3 is headed toward the stationary droplet state, i.e., a paraboloid which is cut-off at h = 0 and the proper contact angle and volume. Degradation of mesh quality during the evolution was handled using remeshing and interpolating the previous solution in the new mesh by linear interpolation. In this case this was the most severe source of loss of mass (< 2%). This can be avoided if one chooses an interpolation which maintains a good mesh quality for a longer time or by improving the interpolation of course. The next example shows how we can maintain mesh quality in certain situations somewhat longer. This example uses an initial disc of radius r = 5 and a paraboloid for h0 (x, y ). As a driving force we use gravity with g 1 = 0.2 and g 2 = 1. As in the 1D case this results in convergence to a traveling wave solution, i.e., a droplet moving in +x direction via h(t , x, y ) = F (x − t v , y ). Clearly Fig. 4 shows that the motion of the support is mainly a translation with constant velocity v in x direction. However, in the simulation in the middle panel of this figure we only use the normal component of ψ˙ to move the domain.

776

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

Fig. 4. Droplet under the action of gravity at t = 0 (left) and at t = 14 (middle, right). Comparison between choice of the tangential component (middle) none (right) optimal constant velocity according to (10).

Fig. 5. Motion of droplets due to gravitational force g 2 = 1 and different g 1 at increasing times (left to right).

In order to find out what this implies let us for a moment assume that the domain is still a ball of radius r, which is transported using ψ˙ = ( v ex · n)n. Then for points on the boundary we get in a comoving frame centered at (0, 0)

ψ˙ −

  v 0

=v

x r2

  x y

 

−

v 0

=

v r2



− y2 xy



,

so there is the clear tendency to transport points from the front of the droplet to its back. This is exactly what happens in the middle panel of Fig. 4. In order to prevent this effect, we choose the tangential part of ψ˙ more carefully. In particular we want this method to work with traveling wave solutions. We do this by finding ( v x , v y ) ∈ R2 , so that

 

2 ( v x , v y ) − x˙ · n d ,

(10)

γ

is minimal. The remaining part we choose so that ψ˙ = ( v x , v y ) + α ∇ h where the boundary values for α are chosen accordingly and α is harmonic inside ω . The resulting deformation at t = 14 is shown in the right panel of Fig. 4. This result is obtained without using any remeshing. With stronger gravity g 1 the droplet shape again becomes more asymmetric in time and real droplets would eventually decay into smaller satellite droplets after some finite time. This can be seen in Fig. 5, where we combine the optimized transport with remeshing every 50 time-steps. For the largest value g 1 = 0.8 the final solution h is close to a pinch-off. In Fig. 6 we show the same investigation performed with the mixed mobility m(h) = 13 h3 + 0.1h2 , where the slip length b = 0.1 is included in order to avoid a slowdown of the contact line as x → 0. We basically observe the same instability mechanism, i.e., at sufficiently large tangential gravitational forces g 1 satellite droplets pinch off. Compared to the case α = 2 this happens already at smaller values of g 1 . Some aspects of droplet shapes on inclined substates have been investigated both experimentally and theoretically [27,28]. In particular there exists a corner singularity at the rear end of the droplet, which is observed only in a certain parameter regime. Otherwise droplets are stable small volumes or develop slender rivulets for large volumes. Thus, a systematic study of possible droplet shapes for different values of g 1 , g 2 is certainly an interesting question for future studies. 5. Conclusion We presented a novel algorithm with which one can solve the thin-film free boundary problem with an explicit contact line efficiently in 1D and 2D. Therefore we used a gradient-type finite element method, together with a practical method to enforce the kinematic condition. The method is very robust in the sense that it works on fine but also on very coarse

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

Fig. 6. Motion of droplets due to gravitational force as in Fig. 5 but with m(h) =

777

1 3 h 3

+ 0.1h2 .

meshes and moderate time-steps. In two spatial dimensions mesh quality can be an issue. We addressed this issue by showing that sometimes one can use the ambiguity in the tangential component of ψ˙ to obtain a mesh update, by which the mesh quality does not degrade too quickly. The feasibility of the approach was underlined by providing examples of droplets driven by surface tension and gravity. Acknowledgements I am thankful for many fruitful discussions with Robert Huth, with whom I started to discuss different ways to tackle this problem. The extension of this work to liquid bilayers is joint work with Sebastian Jachalski and Georgy Kitavtsev, who I wish to thank as well. The financial support by the Research Center Matheon through project C 10 (-2014) and the Einstein Center for Mathematics in Berlin through project D-OT 1 (2014-) is kindly acknowledged. Appendix A. Detailed explanation of FE method and ALE decomposition in 1D

˙ from a finite In order to compute approximate weak solutions of (6a), (6c) we choose representations of h, h˙ and H dimensional  −1 function space as follows. First we decompose the domain ω ⊂ R into a union of mutually disjoint intervals ω = nk= 1 I k . Then we define V = { v ∈ C (ω ) : v affine on each I k } the space of piecewise affine continuous functions. Note that V has no restriction due to essential boundary conditions. Letting I k = (xk , xk+1 ) and with x− = x1 and x+ = xn we have a decomposition of ω = (x− , x+ ) into n − 1 intervals. The space V has the dimension dim V = n and we can construct basis functions φi (x) by defining φi (x j ) = δi j where the φi are piecewise linear on each I k . Now we can write solutions as n n linear combinations of φi with real coefficients u i , p i so that h˙ (x) = i =1 u i φi (x), π (x) = i =1 p i φi (x), and thereby the real n × n matrices to construct the weak formulation are



Mi j =



Sij =

φi (x)φ j (x) dx, ω

φi (x)φ j (x) dx,

S˜ kij =

ω







m h(t k , x) φi (x)φ j (x) dx,

ω

where ω is based on the values x− , x+ at time t k . Since these so-called hat functions φi are only supported on the set I max(i −1,1) ∪ I min(i ,n−1) , these matrices are all sparse. This gives us the discrete problem Ax = b, where x = (u, p ) ,



A=

M

−τ S

S˜ k M



 

,

b=

0 b˜

,



where b˜ i = ω ∇ h(x) · ∇φi (x) dx + δin |∇ h|−1 ( S + 12 |∇ h|2 )x+ + δi1 |∇ h|−1 ( S + 12 |∇ h|2 )x− . Solving this we get a piecewise linear   ˙ At t = t k we represent h(t k , x) = n hi φi (x) and H˙ (x) = n U i φi (x) in V with the difference that representation for h. i =1 i =1 ˙ H˙ ∈ V but h is piecewise constant it is impossible to by definition h1 = hn = 0 and U 1 = U n = 0. Since we assumed h,  interpret (6b) pointwise we do instead is to project h on a linear function g ∈ V using g , φ

 = − h, φ for  inside ω . What  − 1 all φ ∈ V . Using g = i g i φi , h = i h i φi then this can be computed by solving g = M Ch where C i j = ω φi (x)φ j (x) dx. Having computed this g and using ξi = (xi − x1 )/(xn − x1 ) we construct two matrices

P i j = δ1 j (1 − ξi ) g i + δnj ξi g i ,

I i j = δi j − δ1i δ1 j − δni δnj ,

X i j = δ1 j (1 − ξi ) + δnj ξi ,

and then have ( I − P )U = u, where U = (˙x− , U 2 , . . . , U n−1 , x˙ + ) and u = (u 1 , . . . , un ) . Now we can solve U = ( I − P )−1 u and then update the nodal values and the mesh according to

H(t k+1 ) = H(t k ) + τ I U, xk+1 = xk + τ X U,



where H i (t k ) = h t k , xi (t k ) , x = (x1 , . . . , xn ) , H(t k ) = 0, H 2 (t k ), . . . , H n−1 (t k ), 0 .

778

D. Peschka / Journal of Computational Physics 295 (2015) 770–778

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]

J. Eggers, Toward a description of contact line motion at higher capillary numbers, Phys. Fluids 16 (9) (2004) 3491–3494. A. Oron, S.H. Davis, S.G. Bankoff, Long-scale evolution of thin liquid films, Rev. Mod. Phys. 69 (3) (1997) 931. C. Huh, L. Scriven, Hydrodynamic model of steady movement of a solid/liquid/fluid contact line, J. Colloid Interface Sci. 35 (1) (1971) 85–101. L. Hocking, Sliding and spreading of thin two-dimensional drops, Q. J. Mech. Appl. Math. 34 (1) (1981) 37–55. P. Neogi, C.A. Miller, Spreading kinetics of a drop on a rough solid surface, J. Colloid Interface Sci. 92 (2) (1983) 338–349. P. Constantin, T.F. Dupont, R.E. Goldstein, L.P. Kadanoff, M.J. Shelley, S.-M. Zhou, Droplet breakup in a model of the Hele-Shaw cell, Phys. Rev. E 47 (6) (1993) 4169. L. Zhornitskaya, A.L. Bertozzi, Positivity-preserving numerical schemes for lubrication-type equations, SIAM J. Numer. Anal. 37 (2) (1999) 523–555. G. Grün, M. Rumpf, Nonnegativity preserving convergent schemes for the thin film equation, Numer. Math. 87 (1) (2000) 113–152. J.A. Diez, L. Kondic, Computing three-dimensional thin film flows including contact lines, J. Comput. Phys. 183 (1) (2002) 274–306. F. Bernis, A. Friedman, Higher order nonlinear degenerate parabolic equations, J. Differ. Equ. 83 (1) (1990) 179–206. A.L. Bertozzi, The mathematics of moving contact lines in thin liquid films, Not. Am. Math. Soc. 45 (1998) 689–697. E. Ruckenstein, R.K. Jain, Spontaneous rupture of thin liquid films, J. Chem. Soc. Faraday Trans. 2, Mol. Chem. Phys. 70 (1974) 132–147. G.F. Teletzke, H. Ted Davis, L.E. Scriven, How liquids spread on solids, Chem. Eng. Commun. 55 (1–6) (1987) 41–82. L.W. Schwartz, R.R. Eley, Simulation of droplet motion on low-energy and heterogeneous surfaces, J. Colloid Interface Sci. 202 (1) (1998) 173–188. S. Jachalski, R. Huth, G. Kitavtsev, D. Peschka, B. Wagner, Stationary solutions of liquid two-layer thin-film models, SIAM J. Appl. Math. 73 (3) (2013) 1183–1202. L. Giacomelli, H. Knüpfer, A free boundary problem of fourth order: classical solutions in weighted Hölder spaces, Commun. Partial Differ. Equ. 35 (11) (2010) 2059–2091. F. Otto, Lubrication approximation with prescribed nonzero contact angle, Commun. Partial Differ. Equ. 23 (11–12) (1998) 63–103. M. Bertsch, L. Giacomelli, G. Karali, Thin-film equations with partial wetting energy: existence of weak solutions, Phys. D, Nonlinear Phenom. 209 (1) (2005) 17–27. H. Knüpfer, Well-posedness for the Navier slip thin-film equation in the case of partial wetting, Commun. Pure Appl. Math. 64 (9) (2011) 1263–1296. J. Flitton, J. King, Surface-tension-driven dewetting of Newtonian and power-law fluids, J. Eng. Math. 50 (2–3) (2004) 241–266. L. Giacomelli, M.V. Gnann, F. Otto, Regularity of source-type solutions to the thin-film equation with zero contact angle and mobility exponent between 3/2 and 3, Eur. J. Appl. Math. 24 (05) (2013) 735–760. R. Huth, S. Jachalski, G. Kitavtsev, D. Peschka, Gradient flow perspective on thin-film bilayer flows, J. Eng. Math. (2014) 1–19. G. Karapetsas, R.V. Craster, O.K. Matar, On surfactant-enhanced spreading and superspreading of liquid drops on solid surfaces, J. Fluid Mech. 670 (2011) 5–37. E. Lauga, M. Brenner, H. Stone, Microfluidics: the no-slip boundary condition, in: Springer Handbook of Experimental Fluid Mechanics, Springer, 2007, pp. 1219–1240. G. Karapetsas, R.V. Craster, O.K. Matar, Surfactant-driven dynamics of liquid lenses, Phys. Fluids 23 (12) (2011) 122106. D. Peschka, A MATLAB algorithm for thin-film free boundary problems in one spatial dimension, https://github.com/dpeschka/thinfilm-freeboundary, 2015. T. Podgorski, J.-M. Flesselles, L. Limat, Corners, cusps, and pearls in running drops, Phys. Rev. Lett. 87 (3) (2001) 036102. J. Snoeijer, N. Le Grand-Piteira, L. Limat, H. Stone, J. Eggers, Cornered drops and rivulets, Phys. Fluids 19 (4) (2007) 042104.