- Email: [email protected]

PII: DOI: Reference:

S0045-7825(18)30549-8 https://doi.org/10.1016/j.cma.2018.10.045 CMA 12150

To appear in:

Comput. Methods Appl. Mech. Engrg.

Received date : 1 January 2018 Revised date : 30 October 2018 Accepted date : 31 October 2018 Please cite this article as: Y. Liu, C. Peco and J. Dolbow, A fully coupled mixed finite element method for surfactants spreading on thin liquid films, Computer Methods in Applied Mechanics and Engineering (2018), https://doi.org/10.1016/j.cma.2018.10.045 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

*Manuscript Click here to download Manuscript: LiuEtAl-version3.pdf

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Click here to view linked Referenc

A fully coupled mixed finite element method for surfactants spreading on thin liquid films Yingjie Liua,b , Christian Pecoa,b , John Dolbowa,b,c,d a

Information and Initiative Center, Duke University, Durham, NC 27708, USA Department of Civil and Environmental Engineering, Duke University, Durham, NC 27708, USA c Department of Mathematics, Duke University, Durham, NC 27708, USA d Department of Mechanical Engineering and Material Science, Duke University, Durham, NC 27708, USA b

Abstract A model for the spreading of surfactants over thin liquid films on rough surfaces is presented, along with a novel finite element based discretization. An existing framework for surfactants spreading is first augmented by allowing for the influence of a non-trivial surface roughness on the film height and surfactant concentration. In a standard approach, the fourth order equation for the film height is then split into two second order systems. However, distinct from previous approaches, the proposed method also approximates the surface and depth-averaged velocity fields as independent variables, giving rise to a five-field system of coupled nonlinear equations. Benchmark calculations indicate that this approach allows for a ∆t ∼ h2 scaling without loss of convergence in the Newton algorithm. Consistent with analytical estimates, simulations of surfactant drops spreading over thin liquid films on smooth substrates exhibit a t1/4 temporal scaling for the surfactant leading edge, while constant surfactant sources show faster evolutions and exhibit a t1/2 scaling. Finally, simulations of surfactants spreading over extremely thin films demonstrate that fingering instabilities can be triggered by perturbations to either the film height or the substrate roughness. Keywords: surfactant spreading; lubrication; thin films; mixed method; fingering instabilities 1. Introduction Understanding how surfactants spread is important to a wide range of problems spanning industrial and biomedical applications [1–6]. In this work, we focus on surfactants spreading on thin liquid films. Due to the complexity of the underlying physics and the difficulties of obtaining detailed experimental measurements, a great deal of research on this problem has relied on computational science. In the present study, we generalize a well established lubrication theory for this class of problems to account for the roughness of the underlying surface, and develop a novel mixed finite element method to discretize the governing equations. The accuracy, robustness, and convergence behavior of the method is demonstrated with a range of benchmark studies and numerical investigations. The proposed method successfully captures the fingering instability of surfactants spreading over extremely thin liquid films, and also reveals the sensitivity of the spreading process to the surface roughness of the underlying substrate. When a drop of surfactant is deposited on a thin liquid film, the non-uniform interfacial concentration causes a surface tension gradient that gives rise to Marangoni forces [7–9]. These stresses drive the surfactant to spread rapidly towards the uncontaminated liquid area. For this class of problems, Jensen [10] proposed a lubrication model to capture the evolution of the surfactant concentration and the height of the liquid film. For a reduced order model of the theory, Jensen [10] also provided a similarity solution. Over the Email address: [email protected] (John Dolbow)

Preprint submitted to Computer Methods in Applied Mechanics and Engineering

October 30, 2018

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

past decade, versions of this model have been examined with the aid of numerical simulations in both one and two-dimensional settings [11–19]. Moreover, Karapetsas et al. [20] established a two dimensional model (one horizontal dimension and one vertical dimension) to incorporate the effect of vertical velocity, in which the spreading process of surfactant, the effects of surfactant at the moving contact line, and sorption kinetics above and below the critical micelle concentration are investigated. These studies have provided considerable insight into the complicated physics behind the spreading of surfactants over thin liquid films. The fully coupled lubrication model developed by Jensen [10] presents a number of challenges to robust numerical simulation. In particular, we note the fourth order equation governing the film height and the nonlinear coupling with the surfactant concentration on the surface. The flow of surfactant from a localized initial state often gives rise to sharp gradients in the fields that can be difficult to capture [21–23]. Many of the previous results were obtained using high order methods, e.g., high order Runge-Kutta (RK) schemes, high order finite volume methods (FVM) and high order finite difference methods (FDM) [12, 13, 16, 17, 24]. Because the evolution equation for the film height is a degenerate parabolic equation of fourth order, numerical approximations to the film height can become negative unless special care is taken to prevent such numerical artifacts [25–27]. Considerable effort has gone into the construction of positivity-preserving solutions in the continuous setting [28–31] and the development of positivity-preserving numerical schemes for fourth order degenerate parabolic equations [11, 32, 34, 37]. In the present study, we propose a mixed finite element solver in which the primary fields are approximated by relatively low order basis functions. Our approach to ensuring positivity is based on previous efforts for degenerate fourth order parabolic equations, e.g., thin-film systems [32, 33, 37], lubrication type equations [11, 34, 35, 38] and degenerate Cahn-Hilliard equations [39–41]. Our approach is distinguished by the fact that it treats the fully coupled problem between the surfactant concentration and the film height, and that the velocity fields are treated as independent variables. To our knowledge, a fully coupled finite element method for the lubrication theory has yet to be developed. Perhaps the most advanced model to date is the work of Barrett et al. [11], who relaxed the strong coupling between the surfactant concentration and the surface deformation in order to reduce some of the nonlinearity inherent to the model. A significant challenge for fully coupled methods is the issue of convergence in the nonlinear system of equations at any given time step. Our mixed method treats the spreading velocities as independent variables, and solves for all fields simultaneously in a monolithic manner. Our approach reduces the complexity of the original system, without relaxing any of the nonlinear coupling. In this manuscript, we demonstrate that these choices give rise to a much greater degree of robustness, as evidenced by the ability to employ time steps that are significantly larger than those of some previous schemes developed for this class of problems. When the underlying liquid film is extremely thin, the surfactant spreading is sensitive to the physical process near the contact line, and it can become unstable when the spreading is externally driven. This instability of the contact line can cause fingering phenomena that have been extensively reported [42, 43]. A variety of different mechanisms have been proposed to explain this fingering instability [44–46]. Warner et al. [47] carefully explored the underlying mechanisms and reproduced the fingering pattern with the aid of numerical simulation. The fingers are characterized by high gradients in the primary fields, and capturing them requires the use of either highly refined uniform grids or spatial adaptivity. In this work, the use of relatively low order finite element functions makes the adaptation of spatial adaptivity relatively straightforward. Fingering instabilities represent a bifurcation from a uniform state. In a numerical setting, such bifurcations are often triggered by perturbing the problem in some manner. For example, Warner et al. [47] successfully simulated fingering by introducing a small perturbation to the initial height of the liquid film. In this work, we also examine the manner in which perturbations to the film height can enable this sort of bifurcation. Such considerations led us to hypothesize that the roughness of the underlying surface might also be sufficient to break the symmetry in these types of problems. To investigate this hypothesis, we modified the 2

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Fig. 1: Schematic illustration of a surfactant spreading over a thin liquid film. A drop of surfactant is introduced on top of a thin liquid film resting on a solid substrate with roughness f . The surface tension gradient drives the surfactant to spread over the film while the Marangoni force deforms the liquid film. Here, c denotes the surfactant concentration, and h is the height of the liquid film. The fields v and u represent the surface velocity and depth-averaged velocity, respectively. The size of the indicated roughness relative to the actual film heights studied is not shown to scale. In the simulations presented in this work, roughness functions of maximum magnitude | f |max ≈ 0.01h are examined.

original lubrication model by introducing the surface roughness. This has enabled us to perform simulations of surfactants spreading over rough surfaces and demonstrate that the fingering patterns can be sensitive to their roughness. In a similar manner, Sellier and Panda [50] explored the possibility of controlling the free surface of thin liquid films via modification of the underlying substrate. Their model is distinct from the one we consider here in that it did not involve surfactants and any form of Marangoni flow. Moreover, Sellier and Panda [50] assumed the magnitude of the surface topology to be comparable to the thickness of the liquid film. In the present study, the magnitude of the surface roughness is assumed to be much smaller than the film height, such that the main assumptions underlying lubrication theory still hold. This paper is organized as follows: Section 2 develops the fully coupled lubrication model with the additional consideration of substrate roughness. Section 3 describes the mixed finite element discretization and examines the issue of stability. Representative numerical results, along with a comparison against those obtained in literature, are presented in Section 4. Section 5 is devoted to the study of fingering instabilities and the influence of surface roughness. Finally, some concluding remarks and future lines of research are provided in Section 6. 2. Formulation The presence of surfactant lowers the surface tension of the underlying liquid. When a drop of surfactant is deposited on the surface a liquid film, the variation in the surfactant concentration induces a surface tension gradient which causes the fluid to flow. A schematic illustration of the problem is provided in Figure 1. Variations of this problem have been modeled by Jensen [10] and subsequently investigated by other researchers [12–14, 16, 47]. To our knowledge, however, existing models have assumed that the underlying liquid film is coated on a smooth substrate. As such, the influence of surface roughness has yet to be considered. In this section, a comprehensive model is proposed to characterize a surfactant spreading over a thin liquid film which is deposited on a rough solid substrate. In our derivation, we assume that the magnitude of the surface roughness is small compared to the film height, such that the main assumptions underlying the lubrication theory still hold. 2.1. Model Derivation Let U = (u, v, w) represent the velocity profile within the underlying fluid, and V = (u, v) be the projection of U on the x − y plane. We define P(x, y, z, t) and h(x, y, t) as the pressure and height of the liquid 3

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

film, respectively. Assuming that the spreading velocity is sufficiently small, the flow can be regarded as Stokesian. Stokes flow is characterized by the equation −∇P + µ∆U + ρgk = 0,

(1)

where µ is the viscosity of the fluid, ρ is the density, g is gravity and k = (0, 0, −1). The fluid is assumed to be incompressible, i.e., ∇ · U = 0. (2)

The height of the liquid film is assumed to be small compared to the size of the domain, i.e., = H/L 1, where H and L are representative lengths in the z direction and x − y plane, respectively. As a result of dimensional analysis, w u, w v and u xx uzz , uyy uzz ;

v xx vzz , vyy vzz ;

w xx wzz , wyy wzz ,

(3)

where subscripts denote a partial derivative with respect to the spatial variable. Therefore, Eq. (1) can be rewritten as P x = µuzz ,

(4)

Py = µvzz ,

(5)

Pz = −ρg.

(6)

The lubrication model essentially invokes these assumptions to reduce the fully three-dimensional system to a planar one [48, 49]. As such, it is convenient to invoke planar versions of the gradient ∇ and Laplacian ∆ operators. Accordingly, in the following, we use ∇ p = ∂ x i + ∂y j for the planar gradient operator and ∆ p = ∂ xx + ∂yy for the planar Laplacian operator. The balance of normal stress at the surface of the liquid film (z = h) implies P = Patm + wz + σ∇ · n; at z = h,

(7)

where Patm is the atmospheric pressure, n is the unit outward normal to the surface, and σ represents the surface tension [8]. Considering that wz is small compared to the curvature term and ∇ · n ≈ − h xx + hyy = −∆ p h under the assumption 1, (7) can be rewritten as P = Patm − σ∆ p h; at z = h.

(8)

The surface tension σ is assumed to vary with the surfactant concentration c, in such a way that the following conditions are satisfied: σm ; for c = cm , (9) σ(c) = σ0 ; for c = 0,

where cm is the saturation concentration of the surfactant, and σm and σ0 are the surface tensions of the fully contaminated and uncontaminated liquid surfaces, respectively. The precise form of the relationship between the surface tension σ and the surfactant concentration c will be provided in Section 2.2. Combining Eq. (6) and (8) yields P = ρg(h − z) + Patm − σ∆ p h.

(10)

In many models, the surface tension σ is approximated by the uncontaminated surface tension σ0 [11, 13, 47]. Such approximations neglect the coupling between the pressure and the surfactant concentration. 4

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

In the current model, no such approximation is invoked. Taking the derivative of the pressure with respect to the in-plane coordinates, and substituting it in Eqs. (4) and (5), one can replace the pressure with the second derivative of the in-plane velocity µ

∂2 V = ∇ p ρgh − σ∆ p h , 2 ∂z

(11)

where V satisfies the no-slip boundary condition at the bottom of the liquid, i.e., the rough surface of the solid substrate (z = f ): V|z= f = 0, (12) where f is a function of the coordinates (x, y), representing the roughness of the substrate. We assume that the roughness of the substrate is sufficiently small compared to the film height, such that f h. When f (x, y) = 0, this model reduces to a standard problem with a smooth substrate (z = f = 0). The shear stress at the free surface (z = h) is balanced by the gradient of the surface tension, i.e., µ

∂V = ∇ p σ. ∂z

(13)

Integrating Eq. (11) along the z direction twice and applying boundary conditions (12) and (13), the velocity field V takes the parabolic form: ! 1 1 z− f 2 V= ρg∇ p h − ∇ p σ∆ p h (z − f ) − (h − f )(z − f ) + ∇ p σ. (14) µ 2 µ

We now derive the equations governing the spreading of the surfactant at the surface of the fluid. To that end, we write v = V(x, y, h) for velocity field on the surface. From Eq. (14), this yields v = V(x, y, h) =

h− f (h − f )2 ∇pσ + ∇ p σ∆ p h − ρgh . µ 2µ

(15)

We assume that the surfactant concentration c is governed by mass conservation on the surface: ∂c + ∇ s · J = 0, ∂t

(16)

where ∇ s is the surface operator defined on the deforming surface z = h and can be written as 1 + h2 ∂ − h h ∂ − h ∂ x y x y ∂x ∂y ∂z ∂ ∂ 1 2 ∂ −h h + 1 + h − h ∇ s = (I − nn) · ∇ = x y y x ∂x ∂y ∂z 1 + h2x + h2y ∂ ∂ ∂ −h x − hy + h2x + h2y ∂x ∂y ∂z

(17)

Since h x ∼ H/L = and hy ∼ H/L = , it follows that

and

1 2 2 = 1 − (h + h ) + · · · = 1 + O 2 ; x y 2 2 1 + h x + hy

(18)

1 + h2x = 1 + O 2 ; 1 + h2y = 1 + O 2 ; h x hy = O( 2 )

(19)

Therefore, the surface gradient operator can be reasonably well approximated as ∇ s ≈ ∇ p . Most of the 5

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

existing studies of lubrication theory have adopted this approximation (see, for example, Jensen [10] or Barrett et al. [11] or Warner et al. [47]) and we invoke it here as well. Interested readers can check Stone [51] and Pereira and Kalliadasis [52] for more discussions on the surface operators. The total flux J can be decomposed into a diffusive flux Jdiff = −D∇ p c and an advective flux Jadv = cv, with D the surface diffusivity of the surfactant. Substituting these fluxes into the conservation equation (16), we obtain ∂c + ∇ p · (cv) = D∆ p c. (20) ∂t We define u as the average spreading velocity of the underlying liquid, and calculate it as 1 u(x, y) = h− f

ˆ

h f

V(x, y, z)dz =

h− f (h − f )2 ∇pσ + ∇ p σ∆ p h − ρgh . 2µ 3µ

(21)

Finally, we derive the equations governing the height h of the fluid layer. During the spreading of the surfactant, the surface tension difference induces a deformation of the underlying liquid. Consider a column of fluid with cross section S on the x − y plane. Conservation of mass implies ¨ ˆ ¨ d (h − f )dA = − [(h − f )u]dS = − ∇ p · [(h − f )u]dS . (22) dt S ∂S S Considering that ∂ f /∂t = 0, the governing equation for the film height h is ∂h + ∇ p · (h − f )u = 0. ∂t

(23)

Combining Eqs. (15), (20), (21) and (23) with the relationship between surface tension and surfactant concentration (9), we arrive at the complete system of coupled equations. Provided with suitable boundary and initial conditions, the model is well posed. A zero flux condition is applied on the boundary, i.e., ∇ p c · n = 0; ∇ p (h − f ) · n = 0; ∇ p (σ∆ p h) · n = 0 on Γ.

(24)

With an initial distribution of the surfactant concentration and liquid height, the model is complete. Intermolecular forces (van der Waals forces) are included in some models [10, 11]. Here, we confine attention to cases in which intermolecular forces are negligible compared to the gravitational and Marangoni forces. In fact, there exist a variety of physical phenomena that can influence the evolution of thin liquid films. Oron et al. [48, 49] developed a general mathematical model for lubrication approximations and discussed several particular cases which incorporate the effects of thermocapillarity, evaporation/condensation, movement of the solid substrate under the film, and etc. Interested readers are directed to Oron et al.[48, 49] for additional details. 2.2. Dimensionless Form ˆ f = H fˆ, σ = ξσ, Let x = L xˆ, y = Lˆy, h = H h, ˆ and c = cm cˆ , t = T tˆ. The model proposed above can be rewritten in the following dimensionless form: ∂ˆc + ∇ p · (ˆcvˆ ) = δ∆ p cˆ , ∂tˆ (hˆ − fˆ)2 vˆ = (hˆ − fˆ)∇ p σ ˆ+ κ∇ p σ∆ ˆ p hˆ − βhˆ , 2 h i ∂hˆ + ∇ p · (hˆ − fˆ)uˆ = 0, ∂tˆ 6

(25) (26) (27)

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

uˆ =

(hˆ − fˆ)2 hˆ − fˆ ∇pσ ˆ+ κ∇ p σ∆ ˆ p hˆ − βhˆ , 2 3

(28)

where δ = µD/(ξH), κ = H 2 /L2 and β = ρgL2 /ξ are defined as the diffusive, capillary and gravitational coefficients. Here the characteristic time is defined as T = µL2 /(ξH), and ξ is the difference in surface tension between the surfactant and the underlying liquid, ξ = σ0 − σm . Several different functional relationships between the dimensionless surface tension and the dimensionless concentration have been proposed in the literature [10, 11, 47]. Barrett et al. [11] adopted the simplest version σ(c) = 1 − c, which only suffices for modeling a single layer of surfactant. In the present study, we use a formulation that corresponds to multiple layers [13]: n γ + m (1 − cˆ ) ; for cˆ < 1, σ ˆ = (29) γ; for cˆ ≥ 1,

where γ = σm /ξ, and m is the number of surfactant layers. Unless otherwise specified, we use a cubic decay of the surface tension with concentration (n = 3). Inserting Eqs. (29) and (24) into (32) and (34), we arrive at the result that the velocity fields must satisfy vˆ · n = 0; uˆ · n = 0 on Γ.

(30)

This implies that the flow cannot penetrate the boundary. In the following sections, the hat symbol is omitted for convenience, and all variables correspond to dimensionless quantities. 3. Numerical Solver A complete model was proposed in the last section to characterize the spreading of surfactants on a thin, viscous liquid film. The surfactant concentration, film height and spreading velocities are fully coupled, presenting several challenges for a numerical solver. In the following, a fully coupled mixed finite element solver is proposed to approximate solutions to the present lubrication model. 3.1. Variational Form A rearrangement of the governing equations (25-28) brings considerable convenience in the numerical simulations. Defining the physical height as h p = h − f , and replacing h by h p + f in each of (25) - (28), we obtain: ∂c + ∇ p · (cv) = δ∆ p c, (31) ∂t h2p h2p (32) v = h p ∇ p σ + κ∇ p σ∆ p h p − βh p + κ∇ p θ f , 2 2 ∂h p + ∇ p · h p u = 0, (33) ∂t h2p h2p hp u = ∇ p σ + κ∇ p σ∆ p h p − βh p + κ∇ p θ f , (34) 2 3 3 where θ f = σ∆ p f − β f . This rearrangement indicates that the governing equations of h and h p basically have the same structure, with the exception that the velocity fields include an additional term h2p κ∇ p θ f . In the following of the present manuscript, we will directly solve for h p and then compute h = h p + f . Inserting Eq. (34) into (33) yields 2 h p h3p ∂h p + ∇ p · ∇ p σ ˆ + κ∇ p σ∆ p (h p + f ) − β(h p + f ) = 0. (35) ∂t 2 3 7

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

We note that Eq. (35) is a fourth order degenerate parabolic equation in the physical film height h p whose solution cannot be readily approximated using H 1 functions. Here we follow the standard approach and employ the mixed formulation [11, 32–35, 37, 39–41, 53–55]. We begin by defining the pressure θ as an independent variable θ = σ(c)∆ p (h p + f ) − β(h p + f ).

(36)

Recalling Eq. (10) and ignoring the constant Patm , σ(c)∆ p h corresponds to the pressure at the liquid surface. Distinct from the mixed formulation developed by Barrett et al. [11], in this approach θ is a function of both the surfactant concentration c and the film height h. Compared to Eq. (10), in this case −θ corresponds to the pressure for z = 0. When the underlying liquid is deposited on a smooth substrate (z = f = 0), −θ represents the pressure at the bottom of the underlying liquid. When f , 0, although θ has no clear physical meaning, −θ could be used to approximate the pressure at the bottom surface since f h. In the derivation that follows, we use the following compact form for integrals over the volume: ˆ (37) hu, vi = uvdΩ. Ω

Multiplying the governing equations by wc ∈ C = H 1 (Ω), wh ∈ H = H 1 (Ω), wθ ∈ Θ = H 1 (Ω), integrating by parts and applying the boundary conditions, we obtain E D E

c D w , c˙ − ∇ p wc , cv + δ ∇ p wc , ∇ p c = 0, (38) D E D E wh , h˙ p − ∇ p wh , h p u = 0, E D D E * + D E θ θ θ ∂σ w , θ + ∇ p w , σ∇ p h p + f + w , ∇ p c · ∇ p h p + f + β wθ , h p + f = 0. ∂c

(39) (40)

We treat the velocity fields v and u as independent variables, and note that no derivatives of v and u are involved in the weak form. Multiplying Eqs. (32) and (34) by wvi ∈ V = L2 (Ω) and wui ∈ U = L2 (Ω) (i = 1, 2) and integrating, we obtain the following weak forms for the velocity fields: + D E * E 1 D dσ v v (41) wi , vi − wi , h p c,i − κ wvi , h2p θ,i = 0, dc 2 + D E 1* E 1 D dσ u u wi , ui − wi , h p c,i − κ wui , h2p θ,i = 0. (42) 2 dc 3

3.2. Discretization and Stabilization Scheme

The core variables of the mixed form lubrication system are the surfactant concentration c, the physical underlying film height h p and the pressure θ, which are approximated by piecewise polynomial finite element functions. In particular, we use standard bilinear quadrilateral elements. We denote the finite dimensional approximations for the surfactant concentration, the underlying film height, and the pressure as cA , hAp and θA : cA (x, y, t) =

Nd X k

Nk (x, y)ck (t); hAp (x, y, t) =

Nd X k

Nk (x, y)hk (t); θ A (x, y) =

Nd X

Nk (x, y)θk (t),

(43)

k

where Nk (x, y) are the nodal shape functions, ck , hk and θk are nodal values, and Nd represents the number of nodes in the mesh. 8

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The components of the velocity fields v and u are approximated by piecewise constant functions: viA (x, y, t) =

Ne X

Pk (x, y)vik (t); uiA (x, y) =

k

Ne X

Pk (x, y)uik (t); (i = 1, 2),

(44)

k

where viA and uiA represent the finite dimensional approximations for the velocity fields, and Pk (x, y) are the piecewise constant shape functions. vki and uki are elemental values for dimension i, and Ne represents number of elements. We note that Eqs. (31) and (33) are essentially advection-diffusion and advection equations, respectively. In general, numerical approximations to advection-diffusion systems can lose stability when the ratio of advection to diffusion is sufficiently large. The advection/diffusion balance is characterized by the P´eclet number. In this work, we examine flows with relatively small P´eclet numbers. Nevertheless, in what follows, we examine the use of standard stabilization techniques in the mixed solver. Shortly after the surfactant is deposited on the liquid film, a shock is induced due to the jump in the surface tension gradient. As a result, the spreading velocity can become very large near the edge of the surfactant drop. Our simulations indicate that the peak magnitude of the velocity field can be as large as |v|max ≈ 1. Considering that the diffusive coefficient δ is usually ∼ 10−5 , the element P´eclet number for a mesh size he = 0.1 is approximately |v|he αe = ≈ 104 . (45) 2δ In the present study, we examine the use of the standard Streamline Upwind Petrov Galerkin (SUPG) scheme [56, 57] to stabilize the system. The operator for the concentration c can be written as L = Ladv + Ldiff = v · ∇ p + ∇ p · v − δ∆ p . The SUPG scheme used in the present study is given by ˆ ˆ c τ (Ladv (w ))(˙c + L(c)) = τ (v · ∇ p wc )(˙c + v · ∇ p c + c∇ p · v − δ∆ p c)dΩ, Ω0

(46)

(47)

Ω0

P he ψ where Ω0 = Ne e Ωe , the τ = 2|v| and ψ = coth αe − 1/αe . Because the concentration c and both components of the velocity field vi are approximated by piecewise constant shape functions, ∆ p cA = 0 and ∇ p · vA = 0 within an element. As a result, the SUPG formulation for (31) is given by ˆ D E D E D E D E c A A c A c A A c A w , c˙ + τ v · ∇ p w , c˙ − ∇ p w , c v + δ ∇ p w , ∇ p c + τ (vA · ∇ p wc )(vA · ∇ p cA )dΩ = 0. (48) Ω0

The time independent terms of the bilinear form are BcSUPG (wc , cA )

ˆ D E D E c A A c A = − ∇ p w , c v + δ ∇ p w , ∇ p c + τ (vA · ∇ p wc )(vA · ∇ p cA )dΩ, Ω0

9

(49)

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

satisfying BcSUPG (wc , wc )

ˆ D E D E c c A c c = − ∇pw , w v + δ ∇pw , ∇pw + τ |vA · ∇ p wc |2 dΩ, Ω0 2 ˆ ˆ A v c c2 c = δ |∇ p w | dΩ + κ A · ∇ p w dΩ, |v |

(50) (51)

Ω0

Ω

where κc = |vA |he ψ/2. The stability is ensured even for the advection dominated case. Instabilities are likewise anticipated for the advection Eq. (33) governing the film height. In this case, the complete form with SUPG stabilization is ˆ D E D E D E h ˙A A h ˙A h A A w , h p + τ u · ∇ p w , h p − ∇ p w , h p u + τ (uA · ∇ p wh )(uA · ∇ p hAp )dΩ = 0 (52) Ω0

The time independent terms of the bilinear form are ˆ D E BhSUPG (wh , hAp ) = − ∇ p wh , hAp uA + τ (uA · ∇ p wh )(uA · ∇ p hAp )dΩ,

(53)

Ω0

satisfying BhSUPG (wh , wh )

ˆ ˆ D E h h A A h A h h = − ∇ p w , w u + τ (u · ∇ p w )(u · ∇ p w )dΩ = κ Ω0

Ω0

A 2 u h A · ∇ p w dΩ, |u |

(54)

where κh = |uA |he ψ/2 and ensures stability. The discrete velocity fields vA and uA are coupled to the discrete concentration cA and height hAp fields. Note that, in the dot products with the weight function gradient ∇ p wc and ∇ p wh , we lag the velocity fields by one time step to improve the convergence behavior. In what follows, the lagged fields are indicated by v˜ A , u˜ A . 3.3. Solver Description

We define the residual vector R = Rc , Rv1 , Rv2 , Rh , Ru1 , Ru2 , Rθ , whose components are provided as follows: ˆ D E D E D E D E c c A c A c A c A A c A c R = w , c˙ +τ v˜ · ∇ p w , c˙ − ∇ p w , c v +δ ∇ p w , ∇ p c +τ v˜ A · ∇ p wc vA · ∇ p cA dΩ, (55) Ω

+ D E * 1 v A2 v A dσ A v v A c − κ wi , h p θ,i , R = wi , vi − wi , h p dc ,i 2 ˆ D E D E D E h h ˙A h ˜A h ˙A h A A h R = w , hp + τ u · ∇pw , hp − ∇pw , hp u + τ u˜ A · ∇ p wh uA · ∇ p hAp dΩ, Ω

(56) (57)

+ E 1* D 1 2 u A u A dσ A R = wi , ui − wi , h p c,i − κ wui , hAp θ,iA , (58) 2 dc 3 D E D E * + D E θ θ A θ A θ dσ A A R = w , θ + ∇ p w , σ∇ p (h p + f ) + w , ∇ p c · ∇ p h p + f + β wθ , hAp + f . (59) dc We define the solution vector d = cA , v1A , v2A , hAp , u1A , u2A , θ A . The variational form of the problem can be u

10

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

˙ = 0 for any w ∈ W. Here C A , H A stated as: find d in W = C A ×V1A ×V2A ×H A ×U1A ×U2A ×ΘA , such that R(d, d) A 1 A A 2 and Θ ⊂ H , while Vi and Ui ⊂ L . The resulting nonlinear system of equations is implemented within the Multiphysics Object-Oriented Simulation Environment (MOOSE) [58]. An iterative Newton solver is applied. In each iteration, the full Jacobian matrix is invoked as the preconditioner. The second-order accurate, implicit scheme BDF2 is used in this work for temporal discretization. To ˙ = 0 as describe the method in the current context, we begin by rewriting R(d, d) d˙ = F (d, t) ; t ∈ [0, T ].

(60)

We partition the temporal domain into discrete time steps [0, t1 , t2 , . . . , T ] and denote the solution vector at tn as dn . At time tn+2 , assume that the vectors d0 , d1 , . . . , dn+1 are known. The solution vector at time tn+2 is then given by dn+2 =

w2n+1 (1 + wn+1 )2 1 + wn+1 dn+1 − dn + ∆tn+2 F(tn+2 , dn+2 ) 1 + 2wn+1 1 + 2wn+1 1 + 2wn+1

(61)

where wn+1 = ∆tn+2 /∆tn+1 , ∆tn+2 = tn+2 − tn+1 and ∆tn+1 = tn+1 − tn . In our experience, it is advantageous to employ an adaptive time-stepper for the current system of equations. Accordingly, a temporal adaptivity scheme is employed in the present study. The scheme employs an adaptivity fraction φT as a parameter, where 0 ≤ φT < 1. The adaptive procedure is given by: 1. Advance the simulation from tn to tn+1 with the time step ∆tn+1 = tn+1 − tn . 2. If the nonlinear solve converges, accept the solution and store the elapsed time T n+1 of the current t t step. Calculate the ratio rn+1 of elapsed time T n+1 to the time step ∆tn+1 : rn+1 = T n+1 /∆tn+1 . t t t , ∆t • if rn+1 ≥ rnt and rn+1 ≥ rn−1 n+2 = (1 − φT ) ∆tn+1 .

• otherwise ∆tn+2 = (1 + φT ) ∆tn+1 .

3. If the nonlinear solve fails to converge, halve the time step ∆tn+1 = ∆tn+1 /2 and repeat the calculation. The gradients of the surfactant concentration and liquid film height are large near the spreading front but much smaller over the remainder of the domain. The error near the spreading front dominates the system. Therefore, the mesh is adaptively refined and coarsened during the surfactant spreading process, using a self-similar adaptivity algorithm. In the present study, a self-similar isotropic adaptivity algorithm is adopted based on an error indicator using the jump in field gradients across element edges. We employ two custom markers to guide adaptivity: the refining fraction φR and the coarsening fraction φC . The adaptive procedure is given by: 1. For each element, the jump in the gradient of the element variables is calculated across the element edges and recorded as the element error e . 2. Element errors are sorted and summed to a global error g . 3. The set of elements with the largest errors and whose element errors sum to g φR are refined. 4. The set of elements with the smallest errors and whose element errors sum to g φC are coarsened. Additional details are provided in the MOOSE wiki on adaptivity. 3.4. Positivity-Preserving Numerical Scheme Equation (35) is a degenerate parabolic equation of fourth order in the film height h p . Unless special care is taken in the solution scheme, the numerical approximation for h p can become negative, which is clearly not physical [25–27]. 11

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

The simulations presented in Section 4 do not exhibit negative film heights over the entirety of the spatial and temporal domains considered. Nevertheless, for the sake of completeness, here we propose some modifications that can be incorporated within our method to improve robustness for problems in which conditions could give rise to these types of numerical artifacts. Accordingly, we follow Gr¨un and Rumpf [32] and Zhornitskaya and Bertozzi [34] to develop a numerical scheme that ensures the positivity of the discrete approximations to the solution. We partition the domain Ω into a collection of finite elements comprising the approximated domain Ωh , and use J to denote the node set of Ωh and {p j } j∈J the coordinates of these nodes. Following the work of Zhornitskaya and Bertozzi [34], we define a subspace T h of piecewise bilinear functions and introduce the interpolation operator πh : H 1 (Ω) → T h such that for χ ∈ H 1 (Ω), (πh χ)(p j ) = χ(p j ). We first define the discrete semi-inner product ˆ hχ1 , χ2 iπ =

πh (χ1 (x), χ2 (x))dΩ.

(62)

Ω

Then we arrive at a new discrete variational form, given by: E D E D Eπ D wc , c˙ A − ∇ p wc , cA vA + δ ∇ p wc , ∇ p cA = 0,

(63)

D Eπ D E wh , h˙ Ap − ∇ p wh , hAp uA = 0, D E D E * + D E θ A π θ A θ ∂σ A A w ,θ + ∇ p w , σ∇ p h p + f + w , ∇ p c · ∇ p h p + f + β wθ , hAp + f = 0. ∂c

(64) (65)

Substituting (34) into (64) yields κ D Eπ 1 2 3 ∇ p wh , hAp ∇ p σ − ∇ p wh , hAp ∇ p θ A = 0. wh , h˙ Ap − 2 3

(66) 3

The essential idea of Gr¨un and Rumpf [32] is to replace the nonlinear scalar variable hAp with an element-wise diffusivity matrix with entries given by certain harmonic integral means involving nodal values of the discrete approximation to the solution. We denote the matrix corresponding to the scalar function n hAp with power n as Mn hAp . In the following, we briefly introduce the method to construct this matrix. Additional details can be found in [32]. Let E ∈ Ωh be an arbitrary element in the domain. We assume that there is a point (x0 , y0 ) ∈ R2 , a vector (α1 , α2 ) ∈ R2 and an orthogonal matrix A such that ˆ E = x0 + AE,

(67)

where Eˆ is a convex hull of (0, α1 e1 , α2 e2 ). Here e1 and e2 are the unit vectors of the x and y axes, respectively. In fact, Eˆ is a reference element whose two edges intersect at 0 and form a right angle. Then the element-wise matrix Mn hAp for element E is given by ˆ n hAp | ˆ A−1 , Mn hAp |E = AM E

(68)

ˆ n hAp | ˆ is the diffusivity matrix of the reference element E. ˆ Before providing the explicit form of where M E ˆ n , we introduce the scalar function fν proposed in [32]: M −1 if a , b, (a − b) G0ν,n (a) − G0ν,n (b) fν,n (a, b) := (69) 1/G00 (a) if a = b, ν,n

12

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

A An where Gν,n is an entropy function satisfying 1/G00 ν,n h p = h p and depends on the power n and a small artificial parameter ν > 0. Explicit forms of Gν,n are provided in [32]. Gr¨un and Rumpf [32] proved that discrete approximations to solutions that are constructed using this strategy can preserve positivity if the initial fields satisfy certain conditions. ˆ n hAp | ˆ are given by The entries of the matrix M E " fν,n (h0p , h1p ) A ˆ Mn h p |Eˆ = 0

# 0 , fν,n (h0p , h2p )

(70)

with h0p := hAp (x0 , y0 ), h1p := hAp (x0 + α1 Ae1 , y0 ) and h2p := hAp (x0 , y0 + α2 Ae2 ). Inserting (70) into (68) completes the construction of the matrix Mn hAp . Bringing everything together, the discrete evolution equations for our positivity-preserving scheme are given by: D Eπ D E D E wc , c˙ A − ∇ p wc , cA vA + δ ∇ p wc , ∇ p cA = 0, (71)

with

E κD E D Eπ 1 D wh , h˙ Ap − ∇ p wh , M2 hAp ∇ p σ − ∇ p wh , M3 hAp ∇ p θ A = 0, 2 3 + D E D Eπ D E * A A θ A θ A θ ∂σ w ,θ + ∇ p w , σ∇ p h p + f + w , ∇ p c · ∇ p h p + f + β wθ , hAp + f = 0, ∂c * + D E 1 2 v A π v A dσ A c,i − κ wvi , hAp θ,i = 0. wi , vi − wi , h p dc 2

(72)

(73)

(74)

To complete the formulation, we also incorporate SUPG stabilization as described in Section 3.2 within (71) and (72). 4. Numerical Results We have presented a fully coupled model for simulating the spreading of surfactants over thin liquid films, and developed a robust solver for the coupled nonlinear system of equations. Throughout this section, we consider a series of problems to evaluate the accuracy and robustness of the proposed model and numerical method. Where appropriate, we make comparisons to analytical solutions, experimental observations, and simulations using other methods. 4.1. Similarity Solution When a surfactant is spreading on a thin liquid film laying on a smooth substrate, experiments show that the leading edge of the spreading front follows a R s ∼ t1/4 law. The term involving the surface tension gradient, which represents the Marangoni forces, is dominant in the spreading process. Setting the diffusive, capillary, and gravitational coefficients along with the roughness function to zero ( i.e., δ = β = κ = f (x, y) = 0) yields ∂c + ∇ p · (cv) = 0, v = h∇ p σ(c), (75) ∂t ∂h h + ∇ p · (hu) = 0, u = ∇ p σ(c). (76) ∂t 2 Jensen [10] proposed a similarity solution for the above simplified lubrication system that exhibits a t1/4 scaling on the spreading front. To test the proposed mixed form finite element solver, we perform a

13

2 1.5 h

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

No Stabilization SUPG

1 0.5 0

0

1

2

3

r

4

5

6

Fig. 2: Distribution of the film height at t = 20 along a radial line extending from the origin.

simulation for this simplified lubrication system. We use the values γ = 1.2724 and m = 1 in Eq. (29). A drop of surfactant with radius 1 is introduced at the center of a [−8, 8]×[−8, 8] film with uniform height, q 1 h(x, y, 0) = 1, c(x, y, 0) = (1 − tanh(10(r − 1))) , where r = x2 + y2 . (77) 2

The simulation is performed using a 320×320 mesh of four-node quadrilateral elements and within the temporal range t = (0, 50). Figure 2 illustrates the distribution of the film height at t = 20. We report the results only for the relevant part of the domain (i.e., r = [0, 6]), as the surfactant front has not yet covered a large portion of the liquid surface. As expected form the pure advection nature of Eqs. (75) and (76), the absence of stabilization lead to oscillations in the solution. These oscillations disappear when a SUPG scheme is incorporated, showing the adequacy of this stabilization scheme. While the stability of the results without the SUPG formulation cannot be guaranteed, the differences between the simulation results with and without stabilization appear to be minor. Furthermore, the oscillations are relatively small compared to the magnitude of the approximations, and they do not grow monotonically with time. Additionally, although the initial velocities can be large, they decay quickly as the surfactant spreads. Considering that the reduced model given in Eqs. (75) and (76) represents the most challenging case to deal with, as there is effectively no diffusion, the instability concern for the fully coupled system is not as daunting. In consequence, for the remainder of the numerical examples in this section, no SUPG stabilization is employed in the simulations. With the previous insights about the numerical behavior of the system, we solve the coupled system (75)-(76). To make a direct comparison between our numerical simulations and the similarity solution, we now follow Jensen [10] and assume that σ(c) = 1 − c. The temporal scaling of the surfactant leading edge is reported in Figure 3. A good match´is observed between the numerical simulation and the similarity solution, i.e., R s = (16Qt/π)1/4 where Q = Ω cdΩ is the total mass of the injected surfactant. In other words, upon invoking similar conditions and assumptions, our numerical simulations reproduce the t1/4 scaling in the surfactant front predicted by the similarity solution. 4.2. Allowable Time Step Size and Scaling with Mesh Refinement Given the strong nonlinearity of the system of equations describing the spreading of the surfactant, the size of the time step required to achieve convergence can become very restrictive. In this work, in 14

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Fig. 3: Variation of the surfactant leading edge with time. The leading edge of spreading scales as t1/4 when the amount of injected surfactant is fixed.

order to make a fair comparison with previous results in the literature, we employ a standard NewtonRaphson method with a consistent tangent matrix to linearize the system of equations. Newton’s method can fail to converge for a variety of reasons, most commonly due to an initial guess that is too far from the basin of attraction. Employing a smaller time step can help alleviating this issue, but the problem tends to be exacerbated as the mesh spacing decreases. This issue appears to be particularly acute for previous simulations of surfactant spreading, forcing excessively small time steps in order to obtain convergence. Using a finite volume scheme, Wong [13] reported a time step of ∆t = 10−8 , adopting a fix grid size of ∆x = 0.01 and limiting the number of Newton’s iterations to 10. Moreover, when as many as 2000 iterations were allowed, the maximum time step for their approach only increased slightly, to 10−6 . Momoniat et al. [14] adopted a finite difference solver and used an explicit scheme to discretize the temporal derivative. The maximum time step of this explicit method is limited by the CFL condition ∆t < 1/2h2e . Barrett et al. [11] reported that their finite element based solver could use a uniform time step of ∆t = 10−3 for a mesh size of 0.0625. However, we note that their findings apply to a system with reduced nonlinearities in the governing equations, effectively trading model accuracy for computational robustness. Compared to the aforementioned approaches, the fully coupled model and mixed finite element solver described in this work allow for much larger time steps. To demonstrate this, we consider the problem of an initially circular drop of surfactant with radius 1 introduced on a [-5,5] × [-5,5] liquid surface, given by the following initial conditions: q m h(x, y, 0) = 1, c(x, y, 0) = (1 − tanh (100(r − 1))) , where r = x2 + y2 . (78) 2

To make a fair comparison, the material parameters used here are identical to those in Wong [13], i.e., γ = 0.0, m = 2, κ = 0.015, β = 18.059, δ = 3 × 10−5 , f (x, y) = 0.0, and we also follow the suggestion of that work to rewrite the pressure as P ≈ ρg(h − z) + Patm − σ0 ∆ p h. As a result, the pressure θ simplifies to θ≈

σ0 ∆ p h − βh. ξ

(79)

We set the limit on the maximum number of Newton’s iterations to 10, and explore the largest time step allowing for convergence throughout the simulation. This is repeated for a sequence of simulations using progressively smaller mesh sizes. A comparison of the maximum allowable initial time step for various 15

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

Table 1: Maximum Allowable Time Step for Various Simulation Methods

mesh size he 0.1 0.05 0.025 0.01

finite volume [13] 10−3 2.5 × 10−4 5 × 10−7 10−8

compact finite element 5.0 × 10−4 1.0 × 10−4 3.75 × 10−5 1.0 × 10−5

proposed method 1.25 × 10−2 3.75 × 10−3 1.0 × 10−3 1.25 × 10−4

methods is given in Table 1. The method proposed in this manuscript allows for much larger time steps than the finite volume method suggested by Wong [13]. For a mesh size he of 0.01, the maximum initial time step ∆t for the proposed finite element method is approximately 10−4 , in contrast to 10−8 for the finite volume based method. In general, the method proposed in this manuscript roughly allows for a scaling of ∆t ∼ h2e in the time step. This significant difference in the maximum allowable time step size can be attributed in part to the fact that the proposed finite element method is based on a mixed form of the governing equations. By treating the dimensionless pressure θ as an independent variable, the mixed form circumvents the higher order spatial derivative operators, often associated with convergence issues. Moreover, by treating the surface and average velocities as independent variables, we have reduced the highly nonlinear system to a mixed system with reduced nonlinearity. These strategies effectively yield a significant improvement in the robustness of the solver, and allow for the use of much larger time steps. We complete this study by investigating the performance of a “compact” version of a finite element method, similar to that described in Barrett et al. [11]. In this compact formulation, there are only three independent variables: the surfactant concentration c, the liquid film height h, and the pressure θ. The compact form is given by: * + E D E

c dσ 1 D c w , c˙ − ∇ p w , ch ∇ p c − κ ∇ p wc , ch2 ∇ p θ + δ ∇ p wc , ∇ p c = 0, (80) dc 2 + D E 1* E 1 D h ˙ h 2 dσ w ,h − (81) ∇ p c − κ ∇ p wh , h3 ∇ p θ = 0, ∇pw , h 2 dc 3 + D E D E * D E θ θ θ dσ w , θ + ∇ p w , σ∇ p h + w , ∇ p c · ∇ p h + β wθ , h = 0. (82) dc The maximum allowable initial time step of this compact form finite element method is also listed in Table 1. In this case, for the smallest mesh size examined (he = 0.01), the maximum time step size is still an order of magnitude smaller than the proposed method. This suggests the counter-intuitive conclusion that approximating the surface and depth-averaged velocities separately, and effectively increasing the number of fields from three to five, actually improves the convergence behavior of the nonlinear system. This is a key finding of our work. 4.3. Surfactant Droplet Spreading on a Thin Liquid Film We now apply the mixed finite element method to the fully coupled model described in Section 2. We simulate the spreading of an initially circular drop of surfactant with radius 1, placed at the center of a liquid film with dimensions [-8,8] × [-8,8]. The substrate is assumed to be smooth, so the roughness function f (x, y) = 0. The model parameters are taken to be: γ = 1.2724, m = 1, κ = 0.015, β = 28.325, δ = 3 × 10−5 . The initial mesh size he is 0.1, and the mesh is adaptively refined and coarsened during the simulation. To begin, we examine this problem with the following initial conditions for the film height and surfactant

16

2

2

1.5

1.5

1

1

c

h

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.5 0

0.5

0

1

2

r

3

4

0

5

0

1

(a)

2

r

3

4

5

(b)

Fig. 4: Spreading of a single surfactant drop on a uniform liquid film: (a) height h. (b) concentration c. The curves represent the height and concentration fields over a sequence of time at intervals of ∆T = 0.5.

concentration: h(x, y, 0) = 1,

c(x, y, 0) =

1 (1 − tanh(4(r − 1))) , 2

where r =

q

x 2 + y2 .

(83)

As the surfactant spreads over the surface, the fluid height and concentration evolve with time as shown in Figure 4. The presence of the surfactant results in a significant deformation of the surface, giving as a result a considerable variation in the height of the layer. The liquid height near the center of the domain is smaller due to the compression exerted by the Marangoni force. Due to mass conservation, the compression of the fluid near the center of the domain is compensated by increases of the film height elsewhere. With time, the surfactant spreads outward and covers more surface, decreasing its peak magnitude. When a fixed amount of surfactant is deposited on a thin liquid film, both experiments and simulations indicate that the leading edge of the surfactant scales as t1/4 [10, 13, 16]. However, this temporal scaling changes when a constant source of surfactant concentration is provided and, to our knowledge, this effect has yet to be explored for thin films. We now consider this case by examining the spreading of a surfactant fed by a constant source at the center of a circular domain of liquid film with an outer radius of 100 and initial height of 10. The initial conditions are taken to be: q 1.01 for r <= 1.5, where r = h(x, y, 0) = 10, c(x, y, 0) = x 2 + y2 . (84) 0.01 otherwise, To approximate a constant source of surfactant at the center of the domain, the concentration field is held fixed at c = 1 over an inner circle with radius 1. The temporal evolution of the underlying liquid height and the surfactant concentration are shown in Figure 5. We observe that, during the spreading process, the Marangoni force acts to deform the underlying liquid film. Figure 6 illustrates the variation in the leading edge of the surfactant front with time. In this case, our simulation results indicate that the surfactant front spreads much faster, and that the leading edge scales as R s ∼ t1/2 . Jensen [10] studied the spreading of a surfactant front on a rectangular domain and also 17

10.4

2

10.2

1.5

10

1

c

h

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.5

9.8

9.6

10

20

r

30

40

0

50

(a)

10

20

r

30

40

50

(b)

Fig. 5: Spreading of surfactant from a constant source on a uniform fluid film: (a) height h. (b) concentration c. The curves correspond to distinct dimensionless times t = 1, 2, 4, 8, 17, 40, 61, 119, 241 and 557.

reported a t1/2 scaling law. We note that the construction of the problem is very different here. In addition to the geometric differences, Jensen [10] assumed a concentration source that gives rise to a surfactant concentration at the film center that increases with time. By contrast, in this work, the concentration at the film center remains fixed with time. We contend that this condition is a much more realistic representation of what can readily be studied experimentally [59]. Finally, we note that when a constant surfactant source is introduced over a deep fluid layer, the leading edge follows a t3/4 law [60]. In the current lubrication model, the drag force induced by the (relatively close) solid substrate effectively retards the spreading of the front. 5. Fingering Instability When a drop of surfactant is spreading over an extremely thin liquid film, the behavior of the surfactant drop is very sensitive to physical perturbations near the contact line. Although the contact line is usually stable, a depression in the liquid height will be created behind the contact line if the spreading is driven by an external force. In these cases, the particularly small height of the liquid layer can give rise to fingering instabilities. Fingering instabilities have been observed for two different cases: surfactants spreading on a dry solid substrate [61] and surfactants spreading on a prewetted substrate [42, 43]. The model described in Section 2 is best suited to the latter case, when the liquid film is extremely thin. Figure 7 shows a sequence of images from an experiment reported by Afsar-Siddiqui et al. [43]. These images correspond to the spreading of a uniform drop deposited at the center of a prewetted substrate. As the surfactant spreads, a thinned region first develops near the edge of the drop. At subsequent times, several fingers in the fluid layer start to grow in this thin zone. A similar problem was simulated by Warner et al. [47] using a finite difference method. Fingering instabilities were simulated for surfactant fronts spreading over rectangular domains, using structured Cartesian grids. To our knowledge, the development of fingering instabilities from initially circular drops has not been simulated before. Here, we attempt to elucidate some of the fundamental differences induced by the particular geometry. 18

100

Rs

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

10

Simulation Data Fitted Curve Rs ~ t1/2 1

100

101

102 t

Fig. 6: Scaling of the surfactant front with time for a constant surfactant source applied to the center of a circular domain.

(a) t = 0.067s;

(b) t = 0.14s

(c) t = 0.31s

(d) t = 2s

Fig. 7: Deposition of a surfactant drop on a prewetted solid substrate [43].

5.1. Height Perturbation Induced Instability The origins of fingering instabilities for surfactants spreading over thin films have been examined extensively [44, 61–63]. As the process is inherently unstable, it is sensitive to imperfections. The spreading of the surfactant induces a thinned zone near the edge of the drop. A localized increase in the film height anywhere within the thinned zone leads to a corresponding increase in the local surface velocity and the rate of surfactant transport. The development of a low film thickness zone with relatively high surfactant concentration further exacerbates this behavior (i.e., the height decreases in the shallower regions and increases in the deeper zones). Eventually, deeper regions spread faster and fingers appear. We begin by considering a problem in which instabilities are triggered via perturbation to the initial height of the liquid film, as in Warner et al. [47]. Since the liquid film is assumed to be extremely thin, we neglect gravitational effects, and set β = 0.0. As a result, the pressure θ simplifies to θ = σ(c)∆ p h.

(85)

While the effect of van der Waals forces could be incorporated into the expression for the pressure θ, experiments suggest that fingering occurs in films with thicknesses much greater than those where intermolecular 19

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

forces would be expected to be significant [45, 46, 64, 65]. We consider the problem of a line of surfactant introduced at the left end of a [0,2π]×[0,2π] domain. The model parameters are set to: γ = 1.2724, m = 1, κ = δ = 1 × 10−4 , η = 0.0. Compared to the problem described in Section 4.3, we note that the parameter κ is considerably smaller, consistent with the assumption of a very thin liquid film. To isolate the effect of height perturbation in triggering a fingering instability, we set f (x, y) = 0 everywhere. The initial condition is decomposed as h(t = 0) = h0 + h0 , where h0 is a uniform field and h0 represents the perturbation: N X h0 (x, y, 0) = (1 − x2 + b)H(1 − x) + bH(x − 1), h0 (x, y, 0) = exp −B(x − 1)2 Ak cos(λk y),

(86)

k=1

with B = 5, b = 0.05. Here, H(x) = 0.5(1 + tanh(K x)) represents a regularized Heaviside function. In the following, K = 20 unless specified otherwise. The constants λk and Ak denote the wavenumbers and magnitudes, respectively. We use four (N = 4) cosine functions to perturb the film height in the vertical direction, with {A1 , A2 , A3 , A4 } = {0.01, 0.01, 0.01, 0.005} and {λ1 , λ2 , λ3 , λ4 } = {2, 5, 7, 20}. The initial concentration of the surfactant is set to c(x, y, 0) = H(1 − x). Warner et al. [47] decomposed a lubrication model into a base state and a perturbed state, and simulated each of the two states with the initial conditions h0 and h0 , respectively. They subsequently superimposed the two fields to obtain the final result. In other words, the perturbed fingering solution was superimposed on top of the stable base solution. In the present study, we simulate the evolution of all fields in a single, unified state with the initial condition h(t = 0) = h0 + h0 . Figure 8 shows the height of the simulated surfactant at t = 100. Figure 8(a) shows the finite element results when the mesh size is fixed to π/100, while Figure 8(b) shows the results when the mesh size is π/50 with one level of adaptive refinement and coarsening. The refining and coarsening fractions are set to φR = 0.7 and φC = 0.1, respectively. The initial time step is 0.01 and the time adaptivity fraction is set to φT = 0.1. Note that the smallest element in each case is identical. One can observe that the results of the two cases are nearly indistinguishable, indicating that the initial perturbation to the film height drives the instability, as opposed to numerical artifacts such as roundoff errors. The fingering patterns are also similar to those reported by Wong [13]. Of particular interest is the formation of the thinned zone near the edge of the drop. As the front propagates, several fingers grow from the edge of the drop and appear in this thinned region, in accordance with experimental observations. We now consider a simulation that closely resembles the experiments described in Afsar-Siddiqui et al. [43]. A drop of surfactant with radius r0 = 1 is deposited on a circular domain with radius R0 = 4.5. The initial condition is taken to be: h(x, y, 0) = (1 − r2 + b)H(1 − r) + bH(r − 1) + A(ϕ) exp −B(r2 − 1) , (87)

where A(ϕ) is uniformly distributed between [−0.01, 0.01], B = 5, and b = 0.05. In this expression, r is the distance to the center of the drop and ϕ = arctan(y/x) is the angle. The initial surfactant concentration is assumed to be a circular drop with c(x, y, 0) = H(1 − r). The domain is discretized with 23,798 quadrilateral elements. Adaptive refinement with φR = 0.7 and φC = 0.0 is employed in the simulation to ensure the accuracy of the result, while maintaining efficiency. The initial time step is 0.05 and the time adaptivity fraction is set to φT = 0.1. Figure 9 illustrates the evolution of the film height with time. We note that the initial height is close to a smooth circle with a small perturbation. As time increases and the surfactant spreads over the domain, a thinned zone appears at the edge of the drop. The size of this thinned zone gradually increases with time. Small protrusions start to form on the drop edge and gradually elongate into fingers. Comparing these results to Figure 7, we find that the simulation matches the final fingering pattern reasonably well. 20

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Fig. 8: Mesh sensitivity test for fingering instability: (a) uniform mesh with mesh size he = π/100. (b) Adaptive mesh with background mesh size he = π/50 and one level of refinement.

Zooming into the region close to the drop edge near the bottom right part of Figure 9(d), we obtain a clearer representation of the fingers (c.f. Figure 9(e)). The results illustrate how small perturbations elongate into fingers. Moreover, some fingers branch into two sub-fingers at their tips. This tip-split phenomenon implies a secondary instability which can be observed in Figure 7(d). Hamraoui et al. [66] reported experimental observations of this secondary instability for the spreading of soluble surfactants. They also reported that fingers continued to split and could eventually grow into tree-like structures under certain conditions. This secondary instability of soluble surfactants was numerically reproduced by the subsequent work of Warner and Craster [67] and Edmonstone et al. [68]. In our simulations, we have found it advantageous to employ mesh adaptivitiy as a means to maintain local refinement in the vicinity of the tip of these fingers, where branching events tend to occur. 5.2. Roughness Induced Instability In the previous example, we simulated fingering instabilities triggered by perturbations to the initial liquid layer height. Jensen and Naire [69] emphasized that the behavior of a spreading surfactant is sensitive to the physical process near the contact line. Since the roughness of the substrate can directly affect the position of the contact line, it is reasonable to expect that this type of perturbation can also trigger fingering instabilities. The model developed in this work allows for such an investigation, which we describe here. A surfactant drop is introduced at the left end of a [0,2π]×[0,2π] domain. The roughness function f is set to N X f (x, y) = exp −B(x − 1)2 Ak (cos(λk y) + 1) , (88) k=1

where the constants B, λk , and Ak are taken to be identical to those used in Eq. (86). Periodic boundary conditions are imposed on the top and bottom surfaces of the computational domain. The initial liquid height and surfactant concentration are assumed to be: h(x, y, 0) = (1 − x2 + b)H(1 − x) + bH(x − 1) + A exp −B(x − 1)2 , c(x, y, 0) = H(1 − x), (89)

where A = 0.035. Therefore, the physical height h p of the liquid film at the beginning of the simulation is 21

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a) t = 0

(b) t = 10

(c) t = 50

(d) t = 100

Single Finger

Branching Finger

(e) Fingers near the drop edge: a zoom-in image of Figure (d) Fig. 9: Formation of the fingering pattern. The height field h is shown at various times in the top panels, and a zoom of the bottom right quadrant at t = 100 is shown in the bottom panel. Contour levels for all panels are identical to those shown in the bottom.

given by N X h p (x, y, 0) = h(x, y, 0) − f (x, y) = (1 − x2 + b)H(1 − x) + bH(x − 1) − exp −B(x − 1)2 Ak cos(λk y). (90) k=1

Another case with a smooth substrate but perturbed initial height is constructed as reference. In order to minimize the differences between the rough substrate case and this reference case, the initial physical height ˜ y, 0) = h p (x, y, 0), and c˜ (x, y, 0) = c(x, y, 0). No and concentration field of the reference case are set to h(x, mesh or time adaptivity is used for these two simulations. The patterns at t = 100 for the two cases are illustrated in Figure 10. The simulations confirm that the simulations indicate that a rough substrate also gives rise to a fingering pattern. Comparing Figure 10(a) to 10(b), we find that the two patterns look qualitatively similar. Comparing the fingering patterns induced by the perturbed height and the substrate roughness carefully, clear differences emerge. Figure 10(c) plots the film height h along the vertical coordinate y direction at four different locations for the case of the rough substrate (c.f. Figure 10(a)) and for the case of the smooth substrate (c.f. Figure 10(b)). The variation in the height h is fairly similar inside and outside of the drop (Figure 11 (c)-(1) and Figure 11 (c)-(4), respectively). At the thinned zone (c.f. Figure 11 (c)-(3)), the fingers over the rough substrate are much thicker (and longer) than those of the reference case. Our simulations indicate that the 22

h

h

z=0

f

1.0 1.5 2.0 2.5

(a) 0.25

z=0

f=0

1.0 1.5 2.0 2.5

(b)

Rough Substrate & Flat Initial Film Height Smooth Substrate & Perturbed Initial Film Height

0.2

0.4

0.10.45

0.3

0.050.3

0.2

h

h

h

0.150.6

00.15 0 0

1 0

1

2 2

3 y 3 y

4 4

0.1

5 5

6

0

6

0

1

2

(1) x = 1.0 0.2

0.08

0.15

0.06

0.1

0.04

0.05 0

3 y

4

5

6

(2) x = 1.5

h

h

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

0.02 0

1

2

3 y

4

5

0

6

(3) x = 2.0

0

1

2

3 y

4

5

6

(4) x = 2.5

(c) Fig. 10: Fingering instability comparison at t = 100: (a) rough substrate and flat initial film height. (b) smooth substrate and perturbed initial film height. (c) comparison of h along various slices parallel to the y-axis.

fingers on the rough substrate grow faster than those on the smooth one, which might be attributed to the fact that the substrate roughness can affect the contact line directly. 23

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

To investigate the sensitivity of surfactant spreading to the roughness of the substrate, we decompose the roughness function in the previous example into four independent cases. In each case, the roughness function is perturbed by trigonometric functions with a single mode. The four cases employ the same initial height and concentration fields given in (89). The substrate roughness is given by f (x, y) = A˜ exp −B(x − 1)2 (cos(λy) + 1) , (91) P where A˜ = 4k Ak , and λ = λk in each of the four k cases (k = 1, 2, 3, 4). Figure 11 compares the height h of the liquid films in each of the four cases at t = 100. The results clearly indicate that the wavelength of the perturbation in the surface roughness has a significant impact on the spreading of the surfactant and the deformation of the liquid film. h

(a) λ = 2

(b) λ = 5

(c) λ = 7

(d) λ = 20

Fig. 11: Height of the liquid films (i.e. h) at t = 100 when the films are deposited on the substrates with different roughness patterns. The roughness function shares the same distribution along x direction but has different perturbations along y direction. The perturbation is added using the trigonometry functions with different wave lengths.

In the previous case, the perturbation to the surface roughness is only introduced along the y coordinate and the ensuing instability is triggered along the coordinate direction, as one would intuitively expect. In the following, a more general case is investigated. As before, a circular drop with radius 1 is deposited on a prewetted circular substrate with radius 4.5. However, perturbations are now introduced to the substrate roughness in both the radial and circumferential directions. The roughness function and initial height fields are set to N 1X f (x, y) = Ak (cos(λk r) + 1) + A(ϕ) exp −B (r − 1)2 − δ, (92) 2 k=1 1˜ c(x, y, 0) = H(1 − r), h(x, y, 0) = (1 − r2 + b)H(1 − r) + bH(r − 1) + A, 2

(93)

˜ where A(ϕ) is uniformly distributed between [−0.01, 0.01], B = 5, δ = 0.005, b = 0.05. The constants A, Ak and λk are taken to be identical to those used in the last example. The roughness function is illustrated in Figure 12(a). Since the substrate surface is perturbed along both the radial r and circumferential ϕ directions, we construct a nearly-structured mesh to minimize the effect of hoop round-off error on the radial spreading. One quarter of the mesh is illustrated in Figure 12(b). The radius of the inner circle is 0.8. The circles are partitioned into 160 divisions while the straight lines connecting inner and outer circles are partitioned into 80 divisions. The mesh inside the inner circle is unstructured while the mesh outside is structured. The size of the initial surfactant droplet is larger than that of the inner circle, such that the influence of numerical artifacts in the circumferential direction should be minimal. The refining and coarsening fractions are set to φR = 0.7 and φC = 0.1, respectively. The time step is set to be 0.05 and the time adaptivity fraction φT = 0.0. 24

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a)

(b)

Fig. 12: (a) Roughness function; (b) A quarter of the finite element mesh.

Figure 13 shows the height field at different times during the spreading process. Initially, the liquid height is smooth in both the r and ϕ directions. Shortly after, the surfactant begins to spread over the prewetted substrate, a thin zone starts to form. In contrast to the case in Figure 9, we observe several thick and thin rings due to the perturbation in the radial direction. At later times, the thin zone at the drop edge keeps expanding, and several fingers appear. Finally, short fingers evolve into long fingers (c.f. Figure 13(d)). A zoom of the upper left quadrant at t = 120 is shown in Figure 13(e). Both single fingers and secondary instability induced tip-splits are observed. One can observe an interesting multi-ring pattern in Figure 13(d). The variations in the roughness function and the film height along the x-axis at t = 120 are illustrated in Figure 14. These results indicate that the rings with relatively thick films correspond to locations with larger obstacles on the substrate surface, while the rings with thinner films correspond to regions with smaller obstacles. Larger obstacles usually reduce the spreading velocity and transport rate of the surfactant. Therefore, the surfactant accumulates at locations where larger obstacles exist, and a thick zone is observed. 6. Conclusions We developed a novel model and discretization method for the simulation of surfactants spreading on thin liquid films deposited on rough substrates (i.e., an extended lubrication theory). To our knowledge, this is the first time that the effect of substrate roughness on surfactant spreading has been explicitly incorporated into a model of this type. Our main conclusions are summarized as follows: (1) Compared to previous efforts using finite difference, finite volume, and high order Runge-Kutta methods, the proposed finite element method has more flexibility and can simulate surfactant spreading over complex geometries. Moreover, the proposed method has better convergence behavior and allows for much larger time steps as the mesh size decreases. Compared to a previous finite element based method [11], the present formulation addresses the effect of surfactant concentration on fluid pressure, and captures the full coupling between surfactant concentration and film height. (2) The present model reproduces the similarity solution associated with a simplified model of spreading. When a fixed amount of surfactant is introduced on a thin liquid film, the spreading of the leading edge scales 25

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

(a) t = 0

(b) t = 10

(c) t = 50

(d) t = 120

(e) Zoom of the upper-left quadrant near the drop edge at t = 120 Fig. 13: Formation of the fingering pattern over a circular domain with surface roughness in the radial and circumferential directions. A zoom of the upper right quadrant at the final time is shown in (e).

with time like R s ∼ t1/4 . When this initial condition is replaced by a constant surfactant source, spreading proceeds much more quickly and the leading edge follows R s ∼ t1/2 . (3) The present model exhibits fingering instabilities when the surfactant is spreading on a prewetted substrate. Fingering instabilities are observed due to both perturbations in the initial liquid height as well as the substrate roughness. The roughness of the substrate has a significant influence on the spreading behavior, since it can directly affect the contact line between the liquid and the solid. The final fingering pattern is sensitive to the details of the substrate surface. In summary, the present study proposed a comprehensive model for surfactant spreading. In contrast to previous models, the model considers both the interaction between the surfactant and the underlying liquid film as well as the effect of substrate roughness. The proposed mixed finite element solver appears to deliver better flexibility and convergence behavior. The fact that the surfactant spreading exhibits a sensitivity to the surface roughness suggests that surfactant dynamics might be used to extract information about surfaces whose roughness can be difficult to characterize through other means. The model and solver developed in this work are sufficiently robust to generate forward solutions as part of an inverse formulation. This is an area for future research.

26

100

roughness function: f film height: h physical film: hp = h - f

10-1 height

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

10-2

10-3 0

1

2

x

3

4

Fig. 14: Film height distribution along y = 0.

Acknowledgments Yingjie Liu, Christian Peco, and John E. Dolbow are grateful for the support of NSF grant CMMI1537306, to Duke University. References [1] Leenaars, A. F. M., Huethorst, J. A. M., Van Oekel, J. J., 1990. Marangoni drying: A new extremely clean drying process. Langmuir 6 (11), 1701-1703. [2] Grotberg, J. B., 1994. Pulmonary flow and transport phenomena. Annual Review of Fluid Mechanics 26 (1), 529–571. [3] Wit, A. D., Gallez, D., Christov, C. I., 1994. Nonlinear evolution equations for thin liquid films with insoluble surfactants. Physics of Fluids (1994-present) 6 (10), 3256–3266. [4] Braun, R. J., Snow, S. A., Pernisz, U. C., 1999. Gravitational drainage of a tangentially-immobile thick film. Journal of Colloid and Interface Science 219 (2), 225–240. [5] Matar, O. K., Craster, R. V., 2001. Models for Marangoni drying. Physics of Fluids (1994-present) 13 (7), 1869–1883. [6] Peco, C., Chen, W., Liu, Y., Bandi, M., Dolbow, J., Fried, E. 2017. Influence of surface tension in the surfactant-driven fracture of closely-packed particulate monolayers. Soft Matter 13, 5832-5841. [7] Edwards, D., Brenner, H., 1991. Interfacial transport processes and rheology. Butterworth-Heinemann, Boston. [8] Dettmer, W., Peric, D., 2006. A computational framework for free surface fluid flows accounting for surface tension. Computer Methods in Applied Mechanics and Engineering 195 (23), 3038-3071. [9] Rasthofer, U., Henke, F., Wall, W.A., Gravemeier, V., 2011. An extended residual-based variational multiscale method for two-phase flow including surface tension. Computer Methods in Applied Mechanics and Engineering 200 (21), 1866-1876. 27

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[10] Jensen, O. E., Grotberg, J. B., 1992. Insoluble surfactant spreading on a thin viscous film: shock evolution and film rupture. Journal of Fluid Mechanics 240, 259–288. [11] Barrett, J. W., Garcke, H., N¨urnberg, R., 2004. Finite element approximation of surfactant spreading on a thin film. SIAM Journal on Numerical Analysis 41 (4), 1427–1464. [12] Peterson, E. R., Shearer, M., 2011. Radial Spreading of a surfactant on a thin liquid film. Applied Mathematical Research Express 2011 (1), 1–22. [13] Wong, J., 2011. Simulations of surfactant spreading. Master thesis, Harvey Mudd. [14] Momoniat, E., Rashidi, M. M., Herbst, R. S., 2013. Numerical investigation of thin film spreading driven by surfactant using upwind schemes. Mathematical Problems in Engineering 2013, e325132. [15] Klinteberg, L., Lindbo, D., Tornberg, A.-K.,2014. An explicit Eulerian method for multiphase flow with contact line dynamics and insoluble surfactant. Computers & Fluids 101, 50-63. [16] Swanson, E. R., Strickland, S. L., Shearer, M., Daniels, K. E., 2014. Surfactant spreading on a thin liquid film: reconciling models and experiments. Journal of Engineering Mathematics 94 (1), 63–79. [17] Dieter-Kissling, K., Marschall, H., Bothe,D., 2015. Direct Numerical Simulation of droplet formation processes under the influence of soluble surfactant mixtures. Computers & Fluids 113, 93-105. [18] Barrett, J., Garcke, H., N¨urnberg, R., 2015. Stable finite element approximations of two-phase flow with soluble surfactant. Journal of Computational Physics 297 (15), 530-564. [19] Yang, X., Ju, L., 2017. Linear and unconditionally energy stable schemes for the binary fluid-surfactant phase field model. Computer Methods in Applied Mechanics and Engineering 318, 1005-1029. [20] Karapetsas, G., Richard, V.C., Matar, O.K., 2011. Surfactant-driven dynamics of liquid lenses. Physics of Fluids 23 (12), 122106. [21] Yang, H., Przekwas, A., 1992. A comparative study of advanced shock-capturing schemes applied to burgers’ equation. Journal of Computational Physics 102 (1), 139-159. [22] Jaisankar, S., Raghurama Rao, S., 2009. A central rankine-hugoniot solver for hyperbolic conservation laws. Journal of Computational Physics 228 (3), 770-798. [23] Kou, J., Sun, S., Wang, X., 2015. Efficient numerical methods for simulating surface tension of multicomponent mixtures with the gradient theory of fluid interfaces. Computer Methods in Applied Mechanics and Engineering 292, 92-106. [24] Khatri, S., Tornberg, A.K., 2011. A numerical method for two phase flows with insoluble surfactants. Computers & Fluids, 49, 150-165. [25] Bertozzi, A.L., Brenner, M.P., Dupont, T.F., Kadanoff, L.P., 1994. Singularities and similarities in interface flow. in Trends and Perspectives in Applied Mathematics, L. Sirovich, ed., Applied Mathematical Science. 100, Springer-Verlag, New York, 155-208. [26] Bertozzi, A.L., 1996. Symmetric singularity formation in lubrication-type equations for interface motion. SIAM Journal on Applied Mathematics, 56, 681-714. [27] Almgren, R., Bertozzi, A.L. and Brenner, M.P., 1996. Stable and unstable singularities in the unforced Hele-Shaw cell. Physics of Fluids, 8, 1356-1370. 28

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[28] Bernis, F. and Friedman, A., 1990. Higher order nonlinear degenerate parabolic equations. Journal of Differential Equations, 83, 179-206. [29] Bertozzi, A. L., Pugh, M., 1996. The lubrication approximation for thin viscous films: Regularity and long time behavior of weak solutions. Communications of Pure Applied Mathematics, 49, 85-123. [30] Beretta, E., Berstch, M., Dal Passo, R., 1995. Nonnegative solutions of a fourth order nonlinear degenerate parabolic equation. Archive of Rational Mechanics and Analysis 129, 175-200. [31] Chugunova M., Pugh, M.C., Taranets, R.M. 2010. Nonnegative Solutions for a Long-Wave Unstable Thin Film Equation with Convection. SIAM Journal on Numerical Analysis, 42, 1826-1853. [32] Gr¨un, G., Rumpf, M., 2000. Nonnegativity preserving convergent schemes for the thin film equation. SIAM Journal on Numerical Analysis, 87, 113-152. [33] Gr¨un, G., Rumpf, M., 2001. Simulation of singularities and instabilities in thin film flow, European Journal of Applied Mathematics, 12, 293-320. [34] Zhornitskaya, L., Bertozzi, A.L., 2000. Positivity-preserving numerical schemes for lubrication type equations. SIAM Journal on Numerical Analysis, 37, 523-555. [35] Gr¨un, G., 2003. On the convergence of entropy consistent schemes for lubrication type equations in multiple space dimensions, Mathematical Computation, 72, 1251-1279. [36] Barrett, J. W., N¨urnberg, R., 2004. Convergence of a finite element approximation of surfactant spreading on a thin film in the presence of van der Waals forces, IMA Journal of Numerical Analysis, 24, 323-363 [37] Barrett, J. W., Blowey, J.F., Garcke, H., 1998. Finite element approximation of a fourth order nonlinear degenerate parabolic equation, Numerische Mathematik, 80, 525-556. [38] Barrett, J. W., N¨urnberg, R., Warner, M. R. E., 2006. Finite element approximation of soluble surfactant spreading on a thin film. SIAM Journal on Numerical Analysis 44 (3), 1218–1247. [39] Barrett, J. W., Blowey, J. F., Garcke, H., 1999. Finite element approximation of the Cahn-Hilliard equation with degenerate mobility, SIAM Journal of Numerical Analysis, 37, 286-318. [40] Barrett J. W., Blowey, J. F., 2001. Finite element approximation of a degenerate Allen-Cahn/CahnHilliard system, SIAM Journal of Numerical Analysis, 39, 1598-1624. [41] Barrett, J. W., Blowey, J. F., Garcke, H., 2001. On fully practical finite element approximations of degenerate Cahn-Hilliard systems, M2AN Mathematical Modeling and Numerical Analysis, 35, 713748. [42] Afsar-Siddiqui, A. B., Luckham, P. F., Matar, O. K., 2003. Unstable spreading of aqueous anionic surfactant solutions on liquid films. Part 2. highly soluble surfactant. Langmuir 19 (3), 703–708. [43] Afsar-Siddiqui, A. B., Luckham, P. F., Matar, O. K., 2003. Unstable spreading of aqueous anionic surfactant solutions on liquid films. Part 1. sparingly soluble surfactant. Langmuir 19 (3), 696–702. [44] Troian, S. M., Herbolzheimer, E., Safran, S. A., 1990. Model for the fingering instability of spreading surfactant drops. Physical Review Letters 65 (3), 333–336.

29

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[45] Cachile, M., Cazabat, A., 1999. Spontaneous spreading of surfactant solutions on hydrophilic surfaces:??? cnem in ethylene and diethylene glycol. Langmuir 15 (4), 1515–1521. [46] Cachile, M., Cazabat, A., Bardon, S., Valignat, M., Vandenbrouck, F., 1999. Spontaneous spreading of surfactant solutions on hydrophilic surfaces. Colloids and Surfaces A 159 (1), 47–56. [47] Warner, M. R. E., Craster, R. V., Matar, O. K., 2004. Fingering phenomena associated with insoluble surfactant spreading on thin liquid films. Journal of Fluid Mechanics 510, 169–200. [48] Oron, A., Bankoff, S.G., Davis, S.H., 1996. Thermal singularities in film rupture. Physics of Fluids 8, 3433-3435. [49] Oron, A., Davis, S.H., Bankoff, S.G., 1997. Long-scale evolution of thin liquid films. Review of Modern Physics 69, 931. [50] Sellier, M., Panda, S., 2010. Beating capillarity in thin film flows. International Journal for Numerical Methods in Fluids 63, 431–448. [51] Stone, H., 1990. A simple derivation of the time-dependent convective-diffusion equation for surfactant transport along a deforming interface. Physics of Fluids A: Fluid Dynamics, 2 (1), 111–112. [52] Pereira, A., Kalliadasis, S., 2008. On the transport equation for an interfacial quantity. The European Physical Journal Applied Physics 44 (22) 211-214. [53] Johnson, C., Pitk¨aranta, J., 1982. Some mixed finite element methods related to reduced integration, Mathematics of Computation 38, 375-400. [54] Wazwaz, A.-M., 1995. On the solution of the fourth order parabolic equation by the decomposition method, International Journal of Computational Mathematics 57, 213-217. [55] Li, J., 2005. Mixed methods for fourth-order elliptic and parabolic problems using radial basis functions, Advances in Computational Mathematics, 23, 21-30. [56] Brooks, A. N., Hughes, T. J., 1982. Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations. Computer Methods in Applied Mechanics and Engineering 32 (1), 199–259. [57] Hughes, T., Franca, C., Mallet, M, 1987. A new finite element formulation for computational fluid dynamics: VI. Convergence analysis of the generalized SUPG formulation for linear time-dependent multidimensional advective-diffusive systems. Computer Methods in Applied Mechanics and Engineering 63 (1), 97–122. [58] Gaston, G., Newman, C., Hansen, G., and Lebrun-Grandi´e, D., 2009. MOOSE: A parallel computational framework for coupled systems of nonlinear equations. Nuclear Engineering Design 239 (10), 1768–1778. [59] Bandi, M. M., Tallinen, T., Nurnberg, L., 2011. Shock-driven jamming and periodic fracture of particulate rafts. Europhysics Letters 96, 36008. [60] Jensen, O. E., 1995. The spreading of insoluble surfactant at the free surface of a deep fluid layer. Journal of Fluid Mechanics 293, 349–378. [61] Starov, V. M., Kosvintsev, S. R., Velarde, M. G., 2000. Spreading of surfactant solutions over hydrophobic substrates. Journal of Colloid and Interface Science 227 (1), 185–190. 30

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 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65

[62] Homsy, G. M., 1987. Viscous fingering in porous media. Annual Review of Fluid Mechanics 19 (1), 271–311. [63] Borgas, M. S., Grotberg, J. B., 1988. Monolayer flow on a thin film. Journal of Fluid Mechanics 193, 151–170. [64] Cachile, M., Schneemilch, M., Hamraoui, A., Cazaba, A.M., 2002. Films driven by surface tension gradients.m. Advanced Colloid Interface Science 96 (1), 59–74. [65] Matar, O., Craster, R., 2009. Dynamics of surfactant-assisted spreading. Soft Matter 5 (1), 3801-3809. [66] Hamraoui, A., Cachile, M., Poulard, C., Cazabat, A., 2004. Fingering phenomena during spreading of surfactant solutions. Colloids and Surfaces A 250 (1), 215-221. [67] Warner, M.R.E. Craster, R, 2004. Fingering phenomena created by a soluble surfactant deposition on a thin liquid film. Physics of Fluids 16 (1), 2933. [68] Edmonstone, B., Craster, R., Matar, O., 2006. Surfactant-induced fingering phenomena beyond the critical micelle concentration. Journal of Fluid Mechanics 564 (1), 105-138. [69] Jensen, O. E., Naire, S., 2006. The spreading and stability of a surfactant-laden drop on a prewetted substrate. Journal of Fluid Mechanics 554, 5–24.

31