A Geometric Diffuse-Interface Method for Droplet Spreading

TL;DR

This paper introduces a geometric diffuse-interface method for modeling droplet spreading, providing new insights into thin-film dynamics through numerical simulations and comparisons with existing models.

Contribution

It presents a novel regularization of the thin-film equation based on geometric gradient flows, offering a new approach to analyze droplet spreading behavior.

Findings

New perspectives on large-scale droplet-spreading dynamics

Numerical simulations comparing different thin-film models

Insights into the solution properties of the proposed regularization

Abstract

This paper exploits the theory of geometric gradient flows to introduce an alternative regularization of the thin-film equation. The solution properties of this regularization are investigated via a sequence of numerical simulations whose results lead to new perspectives on thin-film behavior. The new perspectives in large-scale droplet-spreading dynamics are elucidated by comparing numerical-simulation results for the solution properties of the current model with corresponding known properties of three different alternative models. The three specific comparisons in solution behavior are made with the slip model, the precursor-film method and the diffuse-interface model.

Tables

Table 1. Table 1: Summary of the various small length scales used in the different regularisation methods

Method	Lengthscale
Navier Slip Model	Slip length
Attractive / Repulsive Potential	Precursor-film thickness
Diffuse-Interface Method	Diffuse-interface thickness
The present method	Length scale of smoothing kernel

Equations

$$\frac{\partial h}{\partial t}=-\frac{\partial }{\partial x}\left(h^n\frac{{\partial }^{3}h}{\partial x^3}\right),t>0,x\in (-\infty ,\infty ),$$

$$h(x,t=0)=h_0(x),h_0(x)\ge 0,x\in (-\infty ,\infty ),$$

$$\frac{\mathrm{d}M}{\mathrm{d}t}=0,M={\int }_{-\infty }^{\infty }h(x,t)\mathrm{d}x.$$

$$h(x,t)=t^{a}f(x/t^a),a=\frac{1}{n+4}.$$

$$f^n{f}^{\prime \prime \prime }=\frac{nf}{n+4}.$$

$$\frac{\partial h}{\partial t}=-\frac{\partial }{\partial x}\left[\left(h^3+\lambda h^2\right)\frac{{\partial }^{3}h}{\partial x^3}\right],t>0,x\in (-\infty ,\infty ),$$

$$u(x,y=0,t)={\lambda }_{\ast }{\left(\frac{\partial u}{\partial y}\right)}_{y=0},$$

$${\left(\frac{\partial h}{\partial x}\right)}_{x=x_{\mathrm{m}}}^3\sim \left(\frac{\mathrm{d}{x}_{\mathrm{m}}}{\mathrm{d}t}\right)\ln \left[\frac{x_{\mathrm{m}}}{\lambda b}{\left(\frac{\mathrm{d}{x}_{\mathrm{m}}}{\mathrm{d}t}\right)}^{1/3}\right],$$

$$\left(\frac{\mathrm{d}{x}_{\mathrm{m}}}{\mathrm{d}t}\right)t^{6/7}=\text{Const.}+\text{Logarithmic corrections},$$

$$\overline{h}(x,t)={\int }_{-\infty }^{\infty }K(x-y;\alpha )h(y,t)\mathrm{d}y:=K\ast h,$$

$$E=\frac{1}{2}{\int }_{-\infty }^{\infty }|{\partial }_{x}h{|}^2\mathrm{d}x.$$

$$\frac{\partial h}{\partial t}=\frac{\partial }{\partial x}\left[h\mu (h)\frac{\partial }{\partial x}\frac{\delta E}{\delta h}\right],\mu (h)=h^2.$$

$$\frac{\mathrm{d}E}{\mathrm{d}t}=-{\int }_{-\infty }^{\infty }h\mu (h)|{\partial }_{xx}h{|}^2\mathrm{d}x\le 0.$$

$$\bar{E}=\frac{1}{2}{\int }_{-\infty }^{\infty }|{\partial }_x\overline{h}{|}^2\mathrm{d}x,$$

$$\frac{\delta \bar{E}}{\delta h}=K\ast \frac{\delta E}{\delta \overline{h}}=-K\ast {\partial }_{xx}\overline{h},$$

$${\partial }_{t}h=-{\partial }_x\left[h\overline{\mu }(h,\overline{h}){\partial }_{x}K\ast {\partial }_{xx}\overline{h}\right],$$

$$⟨\zeta \diamond \nu ,u⟩:=⟨\zeta ,-{\pounds }_u\nu ⟩,$$

$$a_t=-{\pounds }_{{\left(\mu (a)\diamond \frac{\delta E}{\delta a}\right)}^{♯}}a,$$

$$⟨\zeta ,a_t⟩=⟨\frac{\delta E}{\delta a},\delta a⟩,\text{ with }\delta a=-{\pounds }_{{\left(a\diamond \zeta \right)}^{♯}}\mu (a).$$

$$a(x,t)=\sum_{i=1}^{N}\alpha_{i}(t)\delta(x-q_{i}(t))\,,\qquad\text{ with }\qquad\dot{q}_{i}=\left(\mu(a)\diamond\frac{\delta E}{\delta a}\right)\bigg{|}_{x=q_{i}}\,.$$

$$\dot{q}_{i}(t)=\Big{[}\bar{\mu}(h,\bar{h})\partial_{x}\big{(}\bar{h}-\bar{\bar{h}}\big{)}\Big{]}_{x=q_{i}(t)}\,,$$

$$\overline{\mu }(h,\overline{h})=\frac{3}{2}h\overline{h}-\frac{1}{2}{h}^2.$$

$$K\ast f=(1-{\alpha }^2{\partial }_{xx}{)}^{-1}f=\frac{1}{2\alpha }{\int }_{-\infty }^{\infty }{\mathrm{e}}^{-|x-y|/\alpha }f(y)\mathrm{d}y,$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\int }_{-\infty }^{\infty }\overline{h}(x,t)\mathrm{d}x=0.$$

$$h=\overline{h}-{\alpha }^2{\partial }_{xx}\overline{h}.$$

$${\int }_{-\infty }^{\infty }h(x,t)\mathrm{d}x={\int }_{-\infty }^{\infty }\overline{h}(x,t)\mathrm{d}x.$$

$$\partial_{t}\overline{h}=-K*\partial_{x}\bigg{\{}\left[(1-\alpha^{2}\partial_{xx})\overline{h}\right]\mu(h,\overline{h})\partial_{x}\big{(}K*\partial_{xx}\overline{h}\big{)}\bigg{\}}.$$

$$\frac{\bm{\overline{h}}^{n+1}-\bm{\overline{h}}^{n}}{\Delta t}=-K\mathcal{D}_{1}\bigg{\{}\left[(1-\alpha^{2}\mathcal{D}_{2})\bm{\overline{h}}^{n}\right]\bullet(\bm{\overline{h}}^{n})\bullet(\bm{\overline{h}}^{n})\bullet(\mathcal{D}_{1}K\mathcal{D}_{2}\bm{\overline{h}}^{n+1})\bigg{\}},$$

$$h_0(x)=1+ϵ\cos (kx),k=n(2\pi /L),n\in \{1,2,\cdots \}$$

$$\frac{\partial {\overline{h}}_1}{\partial t}=-K\ast {\partial }_x({\partial }_{x}K\ast {\partial }_{xx}{\overline{h}}_1)=-{\partial }_{xxxx}(K\ast (K\ast {\overline{h}}_1)),$$

$$\sigma (k)=-\frac{k^4}{(1+{\alpha }^2{k}^2{)}^2}$$

Full text

I Introduction

This work is concerned with the thin-film equation

[TABLE]

and its modifications. Boundary conditions are chosen as $|x|\rightarrow\infty$ so that Equation (1b) conserves mass:

[TABLE]

The particular value $n=3$ is physically relevant, as then Equation (1b) is a model for the free surface of a viscous thin-film flow. Indeed, equation (1b) with $n=3$ amounts to the Navier–Stokes equations for a thin-film flow, in the limit of lubrication theory. This derivation is developed in Reference oron1997long . A sketch of the physical scenario is given in Figure 1. The case $n=3$ is the subject of the present article. In particular, we revisit the problem of droplet spreading, that is, we seek to model the time evolution of an initial profile $h_{0}(x)$ with compact support.

For a general value of $n$, the spreading of droplets governed by Equation (1b) admits the similarity solution bertozzi1996lubrication

[TABLE]

Substituting this trial solution into Equation (1b) yields the ordinary differential equation

[TABLE]

For $n<3$ Equation (4) possesses smooth solutions with compact support. In this case, equation (4) with initial conditions $f(0)=1$, $f^{\prime}(0)=0$, and $f^{\prime\prime}(0)=-\mu<0$, can be solved by adjusting $\mu$ until $f=f^{\prime}=0$ at some $\eta=\eta_{0}>0$, corresponding to the outermost extent of the droplet. Thus, the position $x_{0}$ of the microscopic contact line where the free surface $h(x,t)$ touches down to zero is described by $x_{0}=\eta_{0}t^{1/(n+4)}$. Unfortunately, this description breaks down precisely for the physically relevant value of $n=3$, at which $f(\eta)$ degenerates into a Dirac delta function centred at $\eta=0$, and the droplet does not spread.

Physically, the breakdown of Equation (1b) as a model of droplet spreading for $n=3$ is due to a small but not ignorable effect which occurs in the vicinity of the microscopic contact line. Namely, the modelling assumptions which enable the passage from the Navier–Stokes equations for a thin film to the simplified free-surface evolution equation (1b) assume there is no relative motion between the film and the underlying substrate. (This is the no-slip condition.) However, the no-slip condition is not consistent with a moving contact line. Many different approaches have been proposed in the literature to restore the missing physics. Three of these approaches will be summarised below. These approaches all exhibit the same qualitative behaviour. However, they each have certain drawbacks. The short summaries of the three main approaches given below will provide the context in which our own proposal for healing the contact-line singularity will be introduced.

Slip-length modelling: A common approach in the modelling literature to resolving the paradox of the moving contact-line is to modify Equation (1b) (with $n=3$) as follows,

[TABLE]

where $\lambda$ is a positive (dimensionless) constant related to the slip length. Equation (5) can be derived using lubrication theory, starting from the Navier–Stokes equations. Instead of imposing the no-slip condition on the velocity component $u(x,y=0,t)$ tangent to the substrate, one instead imposes the condition

[TABLE]

where $\lambda_{*}$ is the dimensional slip length. Working in the limit of lubrication theory, and using appropriate non-dimensionalization, one obtains Equation (5) in this manner.

Typically, one works with $\lambda\ll 1$, as the effect of slip is a small but not ignorable. Then, Equation (5) can be solved via the method of matched asymptotic expansions eggers2004characteristic . In this approach, one distinguishes between inner and outer solutions, separated by a lengthscale $x_{\mathrm{m}}(t)$ which demarcates the regions of validity of the different solutions. The geometry of this setup is shown in Figure 2.

In the outer region, with $|x|\ll x_{\mathrm{m}}(t)$, the effect of slip is ignorable, and one can work with $\lambda=0$ in Equation (5). The free-surface profile in this outer region can therefore be well approximated by the similarity solution (4) with $n=3$. By not continuing this solution past $|x|=x_{\mathrm{m}}$, the singularity that occurs in the similarity solution is avoided.

In the inner region, we identify the microscopic contact line $x=x_{0}(t)$ where the free-surface height touches down to zero in a smooth fashion, $h(x_{0},t)=0$ and $h_{x}(x_{0},t)=0$. As such, the inner region corresponds to $|x-x_{0}|\ll\lambda$. In this region, one therefore solves Equation (5) without the $h^{3}$ term. It is the disappearance of the $h^{3}$ term in this limit which enables the smooth touchdown of the solution at $x=x_{0}(t)$. Finally, the inner and outer solutions are matched at the scale $x_{\mathrm{m}}$. As such, the position $x_{\mathrm{m}}$ is interpreted as the macroscopic contact line, and $\tan\theta_{\mathrm{m}}=h_{x}(x_{\mathrm{m}},t)$ has the interpretation of the macroscopic contact angle. Correspondingly, the position $x_{0}$ has the interpretation of the microscopic contact line, meaning that the microscopic contact angle $\tan\theta_{0}=h_{x}(x_{0},t)$ is zero. It can be noted that this description deals with perfect wetting, such that the droplet spreads indefinitely, and the macroscopic contact angle never reaches a constant value.

Up to prefactors, the matched-asymptotic expansion procedure produced above gives the following relation for the macroscopic contact angle:

[TABLE]

where $b$ is a constant. (Equation (7) was also discovered via a different approach, in Voinov Voinov1977 .) Upon identifying the outer solution on the left-hand side of Equation (7) with the similarity solution (3) $h=t^{-a}h(x/t^{a})$ with $a=1/7$, Equation (7) reduces to

[TABLE]

hence, the leading-order behaviour of Equation (7) is given by $x_{\mathrm{m}}\sim t^{1/7}$, which is the experimentally validated Tanner’s Law oron1997long , valid for a droplet spreading on a perfectly hydrophobic substrate (‘complete wetting’). As such, the leading-order behaviour of the Navier slip model is consistent with experimental findings.

Although the Navier slip model (6) alleviates the singularity in the free-surface height $h(x,t)$ at the microscopic contact line, the higher spatial derivatives of $h(x,t)$ remain singular there. This means the capillary pressure $P=-h_{xx}$ is singular at the microscopic contact line. Although the resulting singularity is rather mild, it does correspond to infinite pressure, which is undesirable in a physical model (although the integral of the pressure, or the force, does remain finite). Also, the prescription of the slip model (6), while convenient from the modelling point of view, does not have an a priori theoretical basis, beyond the obvious connection to atomic-scale fluid-substrate interactions. These two drawbacks have motivated the development of other models of droplet spreading.

Precursor-film modelling: Here, a potential function is prescribed which governs the molecular interactions between the fluid particles and the substrate. Typically, the chosen model potential function has the form of a Van der Waals potential with attractive and repulsive components oron1997long . This model is motivated by a solid understanding of the fluid-substrate interaction. In the lubrication limit, the result of this modelling step is again a single equation for the droplet profile $h(x,t)$; instead of $h_{t}=-\partial_{x}(h^{3}\partial_{x}P)$ with $P=-h_{xx}$ one has $P=-h_{xx}+\Phi(h)$, where $\Phi$ is fluid-substrate interaction potential o2002theory . This description gives a solution for $h(x,t)$ which allows for droplet spreading. Furthermore, the model allows for the possibility of equilibration of the spreading process whereby the film assumes a constant shape characterized by an equilibrium contact angle, thereby modelling partial wetting. Crucially, the value of the equilibrium contact angle is expressed in terms of the coefficients in the Van der Waals potential; hence, the model is fully characterized by basic physical considerations. A further advantage of the model is the droplet-spreading solution maintains a finite stress at the contact line – and yet the leading-order behaviour of the precursor-film and Navier-slip models agree, in a theory of matched asymptotic expansions sui2014numerical . However, a drawback of the model is the requirement for the droplet solution to possess a precursor film – an ultra-thin but non-vanishing film that stretches out indefinitely beyond the droplet core, whose thickness is set by a balance between the attractive and repulsive components in the potential $\Phi$. It can be readily understood that this feature is undesirable in certain applications – for instance, in drop deposition on a substrate, where the substrate should properly be assumed to be initially free from contamination by any fluid film (no matter how thin).

Diffuse-interface modelling: We summarise the diffuse-interface model in the context of the full Navier–Stokes equations, from which the thin-film equation (1b) emerges as a limit. In the full Navier–Stokes equations, the motion of a droplet over a wall can be modelled using such a method. In addition to the velocity and pressure fields, an auxiliary order-parameter variable $C$ (the ‘phase field’) is introduced, which tracks the phases, such that $C=1$ inside the droplet, and $C=0$ outside, in the surrounding phase. There is a smooth transition between these two extreme values across a finite width – hence, a diffuse interface. The phase field evolves according to its own evolution equation, which is typically taken as the Cahn–Hilliard equation ding2007diffuse , thereby resulting in a mathematically consistent framework. By introducing the diffuse interface in this manner (rather than having a sharp interface at the contact line), the model effectively produces slip, through the diffusive fluxes. Hence, the stress singularity at the moving contact line is removed even when a no-slip velocity boundary condition is imposed ding2007wetting .

In this work, we propose a regularisation of Equation (1b) in the spirit of the diffuse-interface method. As such, we propose a modified version of Equation (1b) which depends not only on the free-surface height $h(x,t)$, but also on a diffuse free-surface height,

[TABLE]

where $K(s;\alpha)$ is a smoothing kernel which smooths out small-scale features on a lengthscale $\alpha$ or less. The regularisation is not ad-hoc, instead, it is introduced in the context of a rigorous gradient-energy theory in Section II. The proposed regularisation method has some advantages over the other methods considered herein, in particular:

• •

In contrast to the slip model, the proposed regularisation produces a continuous pressure profile everywhere. We demonstrate below that our new theoretical model amounts to imposing the usual fluid-mechanical interfacial conditions on the filtered free-surface height $\overline{h}$, rather than on $h$. In a context where the free surface $h$ comes into contact with a substrate which has some atomic level of roughness, this makes physical sense, and the use of $\overline{h}$ reflects our uncertain knowledge of the precise location where the droplet, the surrounding medium, and the substrate all come into contact. The parameter $\alpha$ can therefore be viewed as expressing this uncertainty.

• •

In contrast to the precursor-film method, we do not require a precursor film of small-but-finite thickness to extend to infinity. Our model effectively has a precursor film whose thickness falls off to zero at large distances from the droplet core, the falloff scale is exactly $\alpha$.

• •

Although motivated by the diffuse-interface concept, our model does not dispense with the classical description of a sharp interface separating the fluid phases. As such, the sharp interface is still contained in our model; although it is no longer a dynamic variable, and it is recovered from what is effectively a diffuse interface via a deconvolution operation (cf. Equation (8)).

By now one may have recognised that each of the above canonical methods involves a small length scale, as summarised in Table 1. Therefore, it can be seen that each of the methods is concerned effectively with parametrisation missing the small-scale physics, to produce the same large-scale droplet-spreading dynamics in each case. Our model also fits into this framework, as shown in Table 1. We will further illustrate this and other features of our model by using numerical simulations in Section III below.

II Geometric Gradient-Flow structure

The starting-point of the theoretical analysis is to notice that under the dynamics of Equation (1b) the following free-energy functional decays with time:

[TABLE]

Indeed, since $\delta E/\delta h=-\partial_{xx}h$, Equation (1b) can be re-written as

[TABLE]

By multiplying both sides of Equation (10) by $\delta E/\delta h$ and integrating from $x=-\infty$ to $x=\infty$, one obtains

[TABLE]

The proposal for the regularized version of Equation (9) is

[TABLE]

where $\overline{h}$ is the filtered free-surface height given by Equation (8). Then,

[TABLE]

so that equation (10) becomes

[TABLE]

where we have defined $\overline{\mu}(h,\overline{h})$ so that the mobility depends in general on both $h$ and $\overline{h}$.

Equation (10) belongs to a class of geometric gradient-flow equations which appeared in HP2 ; HP3 ; HPT1 . These gradient flows were first inspired by Darcy’s Law for highly viscous flows, which establishes a proportionality relation between the velocity and the force experienced by the fluid. This construction generally applies to arbitrary tensor fields on the configuration manifold $M$ and it involves concepts of Geometric Mechanics, such as Lie derivatives and momentum maps. More specifically, if $V$ denotes the space of tensor fields and $V^{*}$ its dual, one usually considers the duality pairing $\langle\cdot,\cdot\rangle:V^{*}\times V\to\mathbb{R}$ given by the standard $L^{2}-$pairing. This pairing can be used to define a momentum map $\diamond:V^{*}\times V\to\mathfrak{X}(M)^{*}$ whose target space $\mathfrak{X}(M)^{*}$ is identified with the space of one form-densities on $M$, that is $\mathfrak{X}(M)^{*}=\Lambda^{1}(M)\otimes\operatorname{Den}(M)$. In practice, upon defining by $\mathfrak{X}(M)$ the space of vector fields on $M$, this momentum map is defined as

[TABLE]

for any $u\in\mathfrak{X}(M)$ and any couple $(\nu,\zeta)\in V\times V^{*}$. Here, $\pounds_{u}$ denotes the Lie derivative with respect to $u$ and we use the $L^{2}-$pairing on both sides of the equality. Further, we assume that $M$ is a Riemmanniann manifold so that one can define the musical isomorphism $(\,\cdot\,)^{\flat}:\mathfrak{X}\to\mathfrak{X}^{*}$ (flat) and its inverse $(\,\cdot\,)^{\sharp}:\mathfrak{X}^{*}\to\mathfrak{X}$ (sharp). In terms of these operations, a geometric gradient-flow on $V$ is given by an equation of motion of the type

[TABLE]

where $\mu:V\to V^{*}$ is a generalized mobility and $E=E(a)$ is the energy functional, whose functional derivative is denoted by $\delta E/\delta a$. The geometric equation (15) has the following variational formulation HP3 ; HPT1 , which unfolds its gradient-flow nature: for an arbitrary $\zeta\in V^{*}$, one writes

[TABLE]

It is easy to see that Equation (15) follows from the above by integrating by parts and using the definition (14).

Equation (10) belongs to the class of geometric gradient-flow equations (15). This may be shown as follows. Let $V=\operatorname{Den}(\mathbb{R})$, so that $\pounds_{u}a=\partial_{x}(ua)$ and $\zeta\diamond a=a\partial_{x}\zeta$. Then, the sharp operator becomes trivial and equation (15) reduces to (10) for $a=h{\rm d}x$. Notice that in the case of equation (13), we have extended the notion of generalized mobility such that $\mu:V\times V\to V^{*}$. In this case, equation (13) is associated to the regularized version (12) of the energy functional (9). Interestingly enough, the latter belongs to the Burbea-Rao class BR1 ; BR2 of information norms on probability densities. More specifically, the energy functional (9) identifies the norm associated to a $2^{\rm nd}-$order entropy metric in the Burbea-Rao class. We shall leave this connection to information geometry as a direction for future studies.

An important property of geometric-gradient flows of the type (15) is that, when the generalized mobility and the functional derivative $\delta E/\delta a$ are sufficiently smooth, equation (15) admits singular solutions of the type HP3 ; HPT1

[TABLE]

The dynamics of the weights $\alpha_{i}(t)$ can be found on a case-by-case basis by direct substitution. In the case of equation (13), the existence of these solutions depends on the specific expression of $\bar{\mu}(h,\bar{h})$. If this is smooth enough after replacing the singular solution ansatz, then one easily verifies that the weights $\alpha_{i}$ are all constant and

[TABLE]

where we have denoted $\bar{\bar{h}}=K*(K*h)$. We make three remarks about equation (16):

$\bar{\mu}(h,\bar{h})$ is a function(al) of both ${h}(x,t)=\sum_{j}\alpha_{j}\delta(x-q_{j}(t))$ and $\bar{h}(x,t)=\sum_{k}\alpha_{k}K(x-q_{k}(t))$; 2. 2.

$(\bar{h}-\bar{\bar{h}})$ is a function of $(x-q_{j}(t))$ and constant weights $\alpha_{j}$, summed over indices $j$; 3. 3.

After the functional dependence of the mobility $\bar{\mu}$ has been specified, the $x$-dependences of $(\bar{h}-\bar{\bar{h}})(x-q_{j}(t))$ and $\bar{\mu}(h,\bar{h})$ are both evaluated at $x=q_{i}(t)$ to produce a closed dynamical system for the positions $q_{i}(t)$ for each $i=1,2,\dots,N$.

The singular solutions (16) exist, provided $\delta E/\delta h$ is a smooth functional derivative, which holds for our previous energy functional (12). However, the singular solutions also require a smooth generalized mobility. Indeed, the above notation $\bar{\mu}(h,\bar{h})$ is suggestive of an extra smoothing possibly occurring in the mobility function(al). For example, given a mobility function $\mu(h,\bar{h})$, a smooth mobility may be introduced by writing $\bar{\mu}=K*\mu(h,\bar{h})$. In certain cases, previous work has shown HP1 ; HP2 ; HP3 that the singular solutions of geometric gradient-flow equations emerge spontaneously from arbitrary smooth initial conditions and this behavior was exploited to model self-aggregation and alignment of particles with anisotropic interactions HP3 ; HPT2 ; HPT4 ; HONT4 .

In the present work, we shall take an alternative route. Instead of studying singular solution dynamics, we exploit the construction of geometric gradient-flows to introduce a new regularization of the thin-film equation. Then, we shall focus on studying the properties of this regularization without inserting any extra smoothing into the generalized mobility $\bar{\mu}(h,\bar{h})$. While this case precludes the existence of the singular solutions (16), it still leads to new perspectives on thin-film behavior. Here, we choose the functional form

[TABLE]

Equation (13) with the mobility (17) also has a physical basis: it can be obtained from the Navier–Stokes equations and lubrication theory by imposing the interfacial conditions at the smoothened free surface $y=\overline{h}(x,t)$, rather than on the sharp free surface $y=h(x,t)$. In more detail, in the lubrication limit of the Navier–Stokes equations we apply the condition of vanishing shear stress and the matching condition for the droplet pressure to equal the capillary pressure at $y=\overline{h}(x,t)$. It can be noted also that Equations (13) and (17) reduce to the thin-film equation when the convolutions are dropped.

III Numerical Simulations

In this section we explore the solutions of Equation (13) using numerical simulations. For definiteness, we take the filter $K$ to be the inverse of the Helmholtz operator,

[TABLE]

for all continuous, integrable functions on the real line. An advantage of working with the Helmholtz kernel is that it confers on $\overline{h}(x,t)$ the following property:

Theorem 3.1: Under suitable boundary conditions, the integral of the diffuse free-surface height $\overline{h}$ is conserved,

[TABLE]

Proof: Starting with Equation (13), it can be seen that the integral of the bare free-surface height $h(x,t)$ is conserved (subject to appropriate boundary conditions as $|x|\rightarrow\infty$), since the equation for $h_{t}$ is written in conservative form. However, we for the Helmholtz kernel, we have $\overline{h}=(1-\alpha^{2}\partial_{xx})^{-1}h$, hence

[TABLE]

We integrate Equation (19) over the whole real line. Assuming $\partial_{x}\overline{h}\rightarrow 0$ as $|x|\rightarrow\infty$, we obtain

[TABLE]

Hence, since the integral of $h(x,t)$ is conserved, it follows that the integral of $\overline{h}$ is conserved also, and the result follows.

For the purpose of numerical simulations, we further solve Equation (13) on a truncated domain $x\in(-L,L)$ with periodic boundary conditions; this mimics an infinite domain for sufficiently large $L$. The meaning of the Helmholtz kernel (18) in the context of periodic boundary conditions is explained below.

Methodology

Rather than solving Equation (13) directly with the Helmholtz kernel (18), we solve the evolution for the smoothened free-surface height $\overline{h}(x,t)$:

[TABLE]

Therefore, we view $\overline{h}$ as the dynamical variable to be evolved in time. This is more appropriate than working with $h$ as the dynamical variable, as $\overline{h}$ is smoother; hence the numerical method is more stable than would otherwise be the case.

The numerical method used herein is a semi-implicit finite-difference scheme, based on the already-validated method developed elsewhere in a different context by Ó Náraigh and Thiffeault naraigh2010nonlinear . We provide a brief description of this method (and accompanying validations) in what follows. A more detailed development of the numerical methodology (and a comparison with alternative approaches, for instance, the particle method carrillo2019blob ) will be the subject of future work.

As such, we discretize $\overline{h}(x,t)$ on a uniform grid in space and time, with $i\in\{0,1,\cdots,N\}$ labelling the discrete spatial grid points and $n$ labelling the discrete temporal grid points. The spatial grid has a grid spacing $\Delta x$, hence, the spatial grid points are located at $x_{i}=i\Delta x-L$, with $\Delta x=2L/N$. Each partial derivative is discretized using centred finite differences. Hence, $\bm{\overline{h}}^{n}=(\overline{h}(x_{0},t^{n}),\cdots,\overline{h}(x_{N},t^{n}))^{T}$ is a column vector, and the corresponding centred difference operators (with periodic boundary conditions) are $N\times N$ square matrices, denoted here (in an obvious notation) as $\mathcal{D}_{1}$, and $\mathcal{D}_{2}$. Hence, the discretized convolution operator $K$ is itself a matrix, $K=(\mathbb{I}_{N\times N}-\alpha^{2}\mathcal{D}_{2})^{-1}$. In this way, we discretize Equation (20) in the temporal domain as follows:

[TABLE]

where the $\bullet$ denotes pointwise multiplication of vectors, and $K\mathcal{D}_{1}$ etc. denote standard matrix products. Equation (21) can be re-arranged as $M\bm{\overline{h}}^{n+1}=\bm{\overline{h}}^{n}$, where $M$ is an $N\times N$ (invertible) square matrix. Thus, the numerical method is semi-implicit, and $\bm{\overline{h}}^{n+1}$ is extracted from $\bm{\overline{h}}^{n}$ by a matrix inversion at each timestep. The semi-implicit treatment stabilizes the numerical method and allows for a larger timestep than would otherwise be the case naraigh2010nonlinear (the corresponding explicit method involves a fourth-order diffusion operator, which places severe constraints on the timestep for numerical stability).

Validation

We use the model initial condition

[TABLE]

to validate the numerical method. Here, $\epsilon\ll 1$ is a small positive parameter. Physically, this corresponds to a flat interface which is perturbed by a sinusoidal disturbance. This is realistic in the context of either Equation (1b) or (20), as in such a scenario, the initial condition (22) corresponds to damped capillary waves oron1997long . As such, we substitute Equation (22) into Equation (20) and expand the solution $\overline{h}(x,t)=1+\epsilon\overline{h}_{1}(x,t)+\cdots$ in powers of $\epsilon$, keeping only leading-order terms in $\epsilon$. The result is

[TABLE]

which is a linear partial differential equation. We substitute the normal-mode solution $\overline{h}_{1}(x,t)=\mathrm{e}^{\mathrm{i}kx+\sigma t}$ into Equation (23) to produce the dispersion relation

[TABLE]

Motivated by the exact solution encoded in Equation (24), we substitute the initial condition (22) into the full nonlinear numerical partial differential equation (21) and examine the resulting time evolution from the numerical simulation, for a range of values of the wavenumber $k$. For each considered value of $k$, we monitor the disturbance $\Delta(t)=\|\overline{h}(x,t)-1\|_{\infty}$. The result of the numerical simulations is fitted to an exponential decay law $\Delta(t)\propto\mathrm{e}^{-s_{k}t}$, where $s_{k}$ is a fitting parameter, different for each wavenumber $k$. The values of $s_{k}$ are tabulated and the results shown in Figure 3. It can be seen that the dispersion relation thus generated for the numerical results agrees exactly with the analytical dispersion relation in Equation (24), thereby confirming the correctness of our numerical methods.

Results

We solve Equation (21) with the initial condition

[TABLE]

with $y_{0}=0.5$ and $h_{0}=3$. The effect of different initial conditions has been investigated. Specifically, we have also looked at Gaussian initial conditions and a piecewise-defined initial condition, with $\overline{h}(x,t=0)=(1/2)h_{0}(y_{0}^{2}-x^{2})$ inside $|x|<y_{0}$ and $\overline{h}(x,t=0)=0$ outside. We thereby confirm that the following results are robust with respect to the choice of initial conditions. We have also carefully tested the results for numerical convergence: a convergence study is provided in Appendix A.

Sample numerical results are now shown in Figure 4. In Figure 4(a) we demonstrate spacetime plot of the diffuse free surface height $\overline{h}(x,t)$ as it evolves in space and time. The lateral extent of the region where $\overline{h}$ is significantly different from zero (i.e. the droplet) spreads over time. Figure 4(b) shows a snapshot of the free-surface profile at $t=50$. For the purposes of investigation of the numerical results, the macroscopic contact line $x_{\mathrm{m}}(t)$ is defined operationally. As such, $x_{\mathrm{m}}(t)$ is taken to be the realization of the maximum

[TABLE]

The corresponding tangent line is also shown in the figure.

The time evolution of the contact line $x_{\mathrm{m}}(t)$ is shown in Figure 5. It can be seen that $x_{\mathrm{m}}(t)$ behaves as a power law at late times, with $x_{\mathrm{m}}(t)\sim t^{p}$ and $p\approx 0.135$, which is obtained by least-squares fitting. This is very close to the theoretical value $p=1/7$ given by Tanner’s Law.

We next examine the structure of the numerical solution. Figure 6 shows a space-time plot of the solution, this time in similarity variables, with $t^{-1/7}\overline{h}$ plotted on the $z$-axis, using the scaled spatial variable $\eta=x/t^{1/7}$ (the third dimension along the $z$-axis is shown via a contourplot). As such, after transient effects have died away, when viewed on the scale of the computational domain, the solution of the regularized model (20) relaxes to a self-similar functional form.

To understand the results in Figure 6, a plot of $f_{\alpha}(\eta,t):=t^{-1/7}\overline{h}(\eta t^{1/7},t)$, with $\eta t^{1/7}=x$ is shown in Figure 7, for $t=50$. The numerical solution is compared with a solution of the ordinary differential equation $f^{2}f^{\prime\prime\prime}=\eta f/7$, which is the (singular) similarity equation for the un-regularized dynamics (1b). The ordinary differential equation is seeded with the initial condition $f^{\prime}(0)=0$; the additional required initial conditions on $f(0)$ and $f^{\prime\prime}(0)$ are fed in from the numerical solution of the partial differential equation; specifically, $f(0)=f_{\alpha}(0,t=50)$, and $f^{\prime\prime}(0)=f_{\alpha}(0,t=50)$. It can be seen that the profiles of $f_{\alpha}(\eta,t)$ and $f(\eta)$ agree for $|\eta|\ll 1$. Once the macroscopic contact line at $\eta\approx 1$ is reached, the singular nature of the solution of the un-regularised problem becomes apparent, and $f(\eta)$ diverges. One may take this as an equivalent definition of the mascroscopic contact line, i.e. equivalent to the operational definition (26). In contrast, it is precisely in this region where the smooth nature of the solution of the regularized problem begins to appear, and $f_{\alpha}(\eta,t)$ tends to zero as $|\eta|\rightarrow\infty$.

To understand the far-field structure of the diffuse free-surface height, we plot in Figure 8 the numerical value of $\overline{h}(x,t=50)$, on a semilogarithmic scale.

The tail of the profile shows a clear exponential decay $\overline{h}(x,t)\sim\mathrm{e}^{-|x-x_{\mathrm{m}}|/\alpha}$. Indeed, it is clear from Figure 7 and 8 that the late-time solution of the regularised model (20) with the smoothing kernel (18) is a patchwork of two distinct types:

[TABLE]

with the solution judiciously matching between the two extremes such that $x_{\mathrm{m}}(t)\sim t^{1/7}$, in agreement with Tanner’s Law.

Discussion

The numerical solutions use $\overline{h}$ as a dynamical variable, the variable $h=(1-\alpha^{2}\partial_{x}^{2})\overline{h}$ therefore plays a passive role. A snapshot of $h$ for the above numerical parameters is shown in Figure 9. From the figure, it is apparent that $h=0$ for $|x|\gg x_{\mathrm{m}}(t)$, this is consistent with $\overline{h}(x,t)\sim A(t)\mathrm{e}^{-|x-x_{\mathrm{m}}(t)|/\alpha}$ for $|x|\gg x_{\mathrm{m}}(t)$, where $A(t)$ is a time-dependent prefactor. We have deliberately retained the numerical gridpoints in Figure 9: this shows an apparent jump discontinuity in $h(x,t)$ near the position of the macroscopic contact line.

The presence of a jump discontinuity in $h(x,t)$ is consistent with the moving contact line. The validity of this claim is demonstrated by the following theorem:

Theorem 3.2: Consider the solution $h(x,t)$ to Equation (13). Regarding the spatial variable $x$, if $h(x,t)$ is in the function class $C^{0}(-\infty,\infty)$ and piecewise differentiable on $(-\infty,\infty)$, then there is no moving-contact-line solution to Equation (13).

Proof: We suppose that there is a moving contact-line solution $h(x,t)=\phi(x,x_{0})H(x-x_{0})H(x+x_{0})$ to Equation (13), where $x_{0}(t)$ is the microscopic contact line, i.e. the minimum positive value of $x$ for which $h(x_{0},t)=0$. Here also, $H(s)$ denotes the Heaviside function, and $\phi(x,x_{0})$ is a differentiable function on the interval $(-x_{0},x_{0})$, which by the assumption of the moving-contact-line solution satisfies

[TABLE]

We apply this solution to Equation (13) and integrate from $x=-\infty$ to $x=\infty$. We obtain

[TABLE]

If $h\in C^{0}(-\infty,\infty)$ in the spatial variable and if $h$ is piecewise differentiable also in the spatial variable, then the integral on the right-hand side can be broken up into different parts and evaluated to give zero, giving $\Phi(x_{0})(\mathrm{d}x_{0}/\mathrm{d}t)=0$, hence $\mathrm{d}x_{0}/\mathrm{d}t=0$, hence no moving contact line.

Therefore, we conclude that the moving contact line in the numerical simulations corresponds to a sharp free-surface height $h(x,t)$ which is piecewise differentiable, but with jump discontinuities. As such, the sharp free-surface height $h(x,t)$ satisfies Equation (13) in a weak sense, with $h\in C^{-1}(-\infty,\infty)$ (i.e. $h$ possesses a finite number of jump discontinuities). Using the convolution (18), we can conclude that $\overline{h}\in C^{1}(-\infty,\infty)$. We furthermore look at the convolution of Equation (13) (i.e. Equation (20)) – this is the equation satisfied by $\overline{h}$. By counting derivatives on both sides, Equation (20) is consistent with $\overline{h}\in C^{1}(-\infty,\infty)$, indicating that $\overline{h}$ satisfies Equation (20) in a strong sense. We emphasize that although this discussion is consistent with the very precise numerical simulations carried out herein, the rigorous proof that $\overline{h}\in C^{1}(-\infty,\infty)$ has not yet been established. This proof will be the subject of future work.

IV Discussion and Conclusions

In summary, we have introduced a regularsed thin-film equation which describes contact-line motion. The method does not rely on a slip length or a precursor film. The method is inspired by the diffuse-interface concept and it involves a smoothened or diffuse free-surface profile $\overline{h}(x,t)$. However, the method still contains a sharp interface, which can be obtained via deconvolution. The method reproduces Tanner’s law for droplet spreading. Based on the numerical results and on counting the derivatives in the regularized thin-film equation, the diffuse profile $\overline{h}$ gives rise to a strong solution of the thin-film equation. However, this should be checked rigorously using theoretical methods (e.g. along the lines of Reference bernis1990higher ). The model in its present guise can also be used as a description for spreading over heterogeneious surfaces, by introducing a spatial dependence into the lengthscale $\alpha$, e.g. letting $K*f=\{1-\partial_{x}(\alpha^{2}(x)\partial_{x}]\}^{-1}f$ in Equation (18).

The model as formulated currently allows only for indefinite droplet spreading, corresponding to spreading on a hydrophobic surface. In order to allow for arrested spreading (and hence, a static contact angle), the model will require the introduction of extra physics, for instance, by adding a body-force potential of the Van der Waals type to Equation (13). Equally, the model may be extended beyond the limit of lubrication theory, by combining the theoretical arguments in Section II with the general level-set formulation of two-phase flow. In this way, it is hoped that the present relatively simple model can serve as a template for geometric diffuse-interface methods for general two-phase flows.

Acknowledgements

This work has been produced as part of ongoing work within the ThermaSMART network. The ThermaSMART network has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska–Curie grant agreement No. 778104. DH was partially supported by EPSRC Standard [Grant number EP/N023781/1], entitled, “Variational principles for stochastic parameterisations in geophysical fluid dynamics”. CT acknowledges support from the Alexander von Humboldt Foundation (Humboldt Research Fellowship for Experienced Researchers) as well as from the German Federal Ministry for Education and Research. LON further acknowledges the UCD Research Sabbatical Leave for Faculty scheme, as well as helpful discussions with Vakhtang Putkaradze.

Appendix A Convergence analysis

In this Appendix, we look at the convergence of the numerical method for the base case considered in Section III with $\alpha=0.05$, $L=2\pi$, $N=500$, and $\Delta t=10^{-2}$. We show the effect of varying the number of grid points $N$ and the timestep $\Delta t$. In this way, we demonstrate that the numerical results shown in Section III are converged. As such, the structure of $\overline{h}(x,t)$ is shown in Figure 10(a) for $\Delta t=10^{-2}$, $t=100$, and various values of $N$. There is no visible change in the structure of $\overline{h}(x,t)$ when $N$ is varied between $250$ and $1000$. Similarly, the position of the macroscopic contact line $x_{\mathrm{m}}(t)$ is plotted in Figure 10(b) for the various values of $N$ between $250$ and $1000$. There is little or no difference between the different plots of $x_{\mathrm{m}}(t)$ versus $t$ for the various values of $N$.

In Figure 11 we further show the time evolution of $x_{\mathrm{m}}(t)$ for fixed $N=500$ and various values of $\Delta t$. Again, there is little or no difference between the different plots of $x_{\mathrm{m}}(t)$ showing that the numerical results presented in the main paper are converged.

Finally, it can be noted that the convergence of the numerical method is rather senstive to the choice of mobility. For instance, using $\overline{\mu}=h^{2}$ rather than Equation (17) leads to non-convergent results. The choice $\overline{\mu}=h^{2}$ corresponds to the Navier–Stokes equations in the lubrication limit with a regularized pressure $p=-K*\partial_{xx}(K*h)$ but the application of the no-stress boundary condition $\partial u/\partial y=0$ on $y=h$ rather than on $y=\overline{h}$. The non-convergence of the numerical results in this instance underlines the importance of using the diffuse-interface $\overline{h}$ consistently in the formulation of the interfacial stress conditions; it also underlines the importance of choosing a mobility function with a physical rationale.