## Abstract

In this paper, a multiscale moving contact line (MMCL) theory is presented and employed to simulate liquid droplet spreading and capillary motion. The proposed MMCL theory combines a coarse-grained adhesive contact model with a fluid interface membrane theory, so that it can couple molecular scale adhesive interaction and surface tension with hydrodynamics of microscale flow. By doing so, the intermolecular force, the van der Waals or double layer force, separates and levitates the liquid droplet from the supporting solid substrate, which avoids the shear stress singularity caused by the no-slip condition in conventional hydrodynamics theory of moving contact line. Thus, the MMCL allows the difference of the surface energies and surface stresses to drive droplet spreading naturally. To validate the proposed MMCL theory, we have employed it to simulate droplet spreading over various elastic substrates. The numerical simulation results obtained by using MMCL are in good agreement with the molecular dynamics results reported in the literature.

## 1. Introduction

In fluid mechanics, the hydrodynamics theory of moving contact line (MCL) has been studied for more than four decades, and it is still a fascinating subject that attracts many research interests. This is because of its potential applications in solving a broad range of scientific and engineering problems such as droplet spreading, capillary motion, cell motility, wetting and wet friction, surfactant assembly, and many other colloidal physics and chemistry phenomena. The standard macroscale hydrodynamics of MCL theory usually employs the so-called no-slip condition as part of the interface boundary condition. Even though this approach has been very successful in solving many fluid mechanics problems, when it comes to the MCL problem, the no-slip condition adopted at the liquid–solid boundary will lead to singularity in shear stress distribution at the vicinity of the MCL front, which poses a serious challenge in solving the MCL problem (e.g. [1–6]).

The no-slip boundary condition referred here is a boundary condition that imposes a zero relative velocity condition at the fluid–solid interface, and it may lead to a non-integrable singular shear stress at the MCL front of the triple interface. The unbounded shear stress implies infinite dissipation, and hence it removes the ability of the MCL theory to solve any practical problems. Because of its profound potential in applications and its intellectual appeal, this challenge has attracted much attention from the field of applied mathematics and fluid mechanics, bioengineering and chemical engineering, as well as computational materials and computational physics.

In fact, the MCL problem is not a purely fluid mechanics or applied mathematical problem, but a chemomechanical problem at multiscale near the interface. In recent years, many attempts have been proposed to solve this problem by considering molecular interaction at the fluid–solid interface (e.g. [7–10]). So far, the atomistic or molecular force-based continuum MCL formulation is involved with complex multiscale boundary conditions, and it is still a work-in-progress in solving actual dynamic wetting problems. Despite much research effort over so many years, we are still looking for a simple, viable and predictive moving contact line solution at continuum scale such that it can match and predict experimental measurement and solve engineering problems. In a recent paper [11], the present authors proposed a multiscale dynamic wetting model (MDWM) that uses a coarse-grained adhesive contact model (CGCM) developed by Sauer & Li [12–14] to avoid the singularity problem resulting from the no-slip boundary condition. However, the model proposed in [11] has several shortcomings. First, the computational cost of the MDWM is very large for practical three-dimensional simulations, because it makes use of a double volume integral to calculate the contact/interaction force between the droplet and the substrate. Second, in the MDWM, the interface adhesion and the interface moving contact line formulations are somewhat disjointed. Last, MDWM does not correctly account for the interface dynamic or inertial effects. In this work, we systematically derive and formulate a rigorous multiscale moving contact line (MMCL) theory that seamlessly couples the MCL theory with the adhesive contact theory by using a fluid interface model, which is an analogue of the Gurtin–Murdoch theory of elastic interface [15–18]. A key feature or component of MMCL is that we are able to convert the double volume integral for interface adhesive traction into a double surface integral that not only reduces the computational cost significantly, but also provide the solution for calculating the interface traction jump across the interface of the triple phase system.

The paper is organized into five sections. In §2, we discussed the basic idea of MMCL and how it works; in §3, the dynamic governing equations of MMCL and its finite-element formulation are presented; in §4, several numerical examples are presented to demonstrate the capability of the MMCL model, and finally we close the presentation in §5 with a few remarks.

## 2. The basic idea of multiscale moving contact line and how it works

To start with, we consider the following triple phase system of gaseous (G), liquid (L) and solid (S) as shown in figure 1. Each interface consists of two surface layers, e.g.
*Γ*_{LS}(L) denotes the liquid surface of the interface *Γ*_{LS} while *Γ*_{LS}(S) denotes the solid surface of the interface *Γ*_{LS}. Similarly, *Γ*_{GL}(G) denotes the gaseous surface of the interface *Γ*_{GL} while *Γ*_{GL}(L) denotes the liquid surface of the interface *Γ*_{LS}; and *Γ*_{GS}(G) denotes the gaseous surface of the interface *Γ*_{GS} while *Γ*_{GS}(S) denotes the solid surface of the interface *Γ*_{GS}. For the triple phase system, the total boundary of each phase is denoted as ∂*Ω*_{α}, *α*=*G*,*L*,*S* and ∂*Ω*_{G}=*Γ*_{GL}(G)∪*Γ*_{GS}(G), ∂*Ω*_{L}=*Γ*_{GL}(L)∪*Γ*_{LS}(L) and ∂*Ω*_{G}=*Γ*_{LS}(S)∪*Γ*_{GS}(S).

### (a) Governing equations of the bulk continuum

The MMCL is a hybrid continuum model. The Lagrange description of continuum mechanics is used to model the dynamic deformation of each phase,
*k* denotes the phase; *σ*_{k} is the Cauchy stress; *ρ*_{k} is the mass density per unit volume; **b**_{k} is the body force per unit mass and **u**_{k} is the displacement field.

Since the stress in the gas phase is negligible, we only consider the equations of motion in phase *k*=*L* and *S*. In this work, the liquid phase is modelled as a compressible Newtonian fluid, with the constitutive relation given as
*κ*_{L} and *μ*_{L} are, respectively, the bulk modulus and viscosity of the liquid phase; ** v** is the velocity of the liquid phase,

*J*=

*ρ*

_{0}/

*ρ*, and

*ρ*

_{0},

*ρ*are the density of the liquid in the reference and current configuration, respectively.

The solid phase is modelled as a St. Venant-Kirchhoff material, and its constitutive equation is given as follows:
**S** is the second Piola–Kirchhoff stress; **F** is the deformation gradient and λ_{S}, *μ*_{S} are the Lame constants of the solid substrate.

### (b) Interaction between two different continuum phases

#### (i) Body–body interaction

The key of the MMCL theory is how it treats the interaction and interface condition between two distinct continuum phases. To do so, we employ a general adhesive contact model that stems initially from Bradley's van der Waals force model [19] and DLVO (Derjaguin and Landau, Verwey and Overbeek) theory [20–22], which supplements the microscale or mesoscale adhesive force condition at the interface.

If we choose the interphase potential as an average between the Lennard-Jones potential (the van der Waals force) and the zeta potential (the double layer electrokinetic potential), we can quantify it when a given interface is specified.

For simplicity, we denote this interatomic potential as *ϕ*(|**x**_{kℓ}|), where **x**_{kℓ}=**x**_{ℓ}−**x**_{k},*k*,ℓ=*G*,*L*,*S*, *k*≠ℓ and **x**_{k}∈*Ω*_{k},**x**_{ℓ}∈*Ω*_{ℓ}. Then the total interphase adhesive-contact potential energy *Π*^{ac} for two interacting phase is
*ϕ*(*r*) is the interatomic potential, and in this work, the 12-6 type Lennard-Jones potential is used as a model to demonstrate the effect of the van der Waals force and double-layer force potentials
*ε* is the potential well (in the unit of energy) and *σ*_{0} is the atomic equilibrium distance.

By simply taking the first variation of the interaction potential energy, we can obtain
*Φ*_{k} and *Φ*_{ℓ} are the homogenized macro interphase interaction potentials between phases *k* and *l*, *k*,ℓ=*G*,*L*,*S*. By doing so, we have expressed the interbody intermolecular forces in the form of a time-dependent body force distribution. For detailed computation implementation of CGCM, readers may consult [12–14].

One may note that the above adhesive contact formulation is essentially a body–body interaction. As shown in equation (2.5), the body–body approach requires a double-layer integral over the two deformable bodies, which needs extensive computation effort, especially for a macroscale three-dimensional computational model. Fortunately, the main effect of molecular adhesion is concentrated around the vicinity of the interface, and one can further convert the double volume integral to a double surface integral.

#### (ii) Surface–surface interaction

As shown in figure 2, we first consider the interface *Γ*_{kℓ}=*Γ*_{kℓ}(*k*)∪*Γ*_{kℓ}(ℓ) between the phases *k* and *l* and choose two arbitrary points **x**_{k}∈∂*Ω*_{k} and **x**_{ℓ}∈∂*Ω*_{ℓ} so that **s**_{kℓ}:=**x**_{ℓ}−**x**_{k}=:−**s**_{ℓk}. Following the approach outlined by Jagota & Argento [23], we can convert the body–body interaction forces directly to surface–surface interaction forces. The adhesive contact force applied on an infinitesimal surface element d*a*_{k}⊂∂*Ω*_{k} due to the presence of an infinitesimal surface element surface d*a*_{ℓ}⊂∂*Ω*_{ℓ} is expressed as
**n**_{k} and **n**_{ℓ} are the unit surface out-normal at points **x**_{k} and **x**_{ℓ}, respectively, *k*,ℓ=*G*,*L*,*S*, *k*≠ℓ, and
*ϕ*(*r*) decays faster than 1/*r*^{3}, so that the integral of the above equation is finite when *r* approaches infinity. For the interatomic potential *ϕ*(*r*) given by equation (2.6), this condition is satisfied, and in fact one can perform the analytical integration and obtain,
*a*_{ℓ}⊂∂*Ω*_{ℓ} due to the presence of an infinitesimal surface element surface d*a*_{k}⊂∂*Ω*_{k} is expressed as
**F**_{k}≠−d**F**_{ℓ}, but the total force applied on phase *k* due to the presence of phase ℓ is equal to that of the force applied on phase ℓ by the phase *k*, i.e. *N*_{k} and *N*_{ℓ} are the total number of infinitesimal surface elements on ∂*Ω*_{k} and ∂*Ω*_{ℓ}, respectively.

Furthermore, one can rewrite equation (2.10) as follows:
*a*_{ℓ}⊂∂*Ω*_{ℓ} at **x**_{k}∈∂*Ω*_{k}. Thus, if we perform the integration over the surface ∂*Ω*_{ℓ}, we may obtain the adhesive surface stress
**x**_{k}∈∂*Ω*_{k}, due to the interaction from the phase ℓ. The corresponding surface traction at point **x**_{k}∈∂*Ω*_{k} can then be expressed as
**x**_{ℓ}∈*Γ*_{kℓ}(ℓ) as

Note that if one of the two interacting phases is rigid, we may be able to analytically integrate equation (2.15) to obtain the adhesive surface stress tensor for a given interatomic potential *ϕ*(*r*). For instance, for the 12-6 Lennard-Jones potential given by equation (2.6), if the phase *l* is a rigid infinite plate as shown in figure 3*a*, then one can obtain an explicit expression for surface stress tensor as
*a* is the shortest distance of the particle **x**_{k} to the plate, **e**_{z} is the unit vector along the *z*-axis. Alternatively, if the phase ℓ is a rigid sphere with radius *R* as shown in figure 3*b*, the surface stress tensor can be found as a diagonal tensor
*a* is the shortest distance of the surface particle **x**_{k} to the centre of the rigid sphere. **e**_{x}, **e**_{y} and **e**_{z} are defined to be the unit vectors along the axes shown in the figure. In general, to calculate the adhesive surface traction, one needs to closely trace the positions of points on the surface and evaluate the current surface unit out-normal for each surface element, and then perform numerical integration.

### (c) Comparison study between conventional moving contact line and the multiscale moving contact line

To understand how the proposed MMCL method works, we have carried out a comparison study by using the conventional MCL method and the MMCL method proposed in the work to simulate the adhesive contact between two continuum objects (shown in figure 4). We would like to mention that the numerical example in figure 4*b* is carried out by using the surface–surface integral-based contact scheme, together with the interface effects. The main advantage to combine the coarse-grained contact model with the MCL line theory is that the interface adhesive-repulsive force can levitate the liquid droplet above the solid substrate (figure 4*b*). By doing so, it completely eliminates singularity problem of the conventional MCL theory, while retaining the surface energy description in dynamics wetting modelling. From a purely mechanics perspective, in the conventional MCL theory hydrodynamics simulation, a singular shear stress arises due to the fact that the initial contact line front between the liquid phase and solid phase forms a crack-shaped cleavage (figure 4*a*), and the abrupt change of surface tangent direction will cause stress concentration. On the other hand, if one can levitate the liquid droplet over the solid substrate, and it will separate the liquid phase and solid phase. This creates a gap between the solid surface and liquid surface and makes the contact line front a traction-free surface. A direct consequence of this approach is that the mathematical idealization induced pathology is eliminated, which is why the proposed method works.

## 3. Dynamic interface moving contact line theory

In this section we discuss the dynamic interface moving contact line theory.

### (a) Interface dynamic equations of the moving contact line

Inspired by the Gurtin–Murdoch theory of elastic interface (e.g. [15–18]), we consider the following surface equations of motion for the surface of all three phases:
**f**^{D,k} is the surface D'Alembert force density, *ς*_{k} is the surface stress, **t**_{k} is the adhesive surface traction vector and ∇_{S} is the surface gradient operator that is defined as
**n** is the unit out-normal of the surface.

In order to couple the surface equations of motion with small-scale adhesive force, we must link the adhesive traction **t**_{k} with the coarse-grain adhesive molecular force. In fact, the adhesive traction force is not the only contribution to the surface traction **t**_{k},*k*=*G*,*L*,*S*. Many other factors are to be considered for the total surface traction, i.e. surface tension, friction, surface strain gradient, surface elasticity gradient, diffusive and other dissipative forces such as surface viscosity and velocity gradients. In the following, an interface MCL theory is proposed to calculate the interface traction at the boundary of each phase. To start with, we first consider the surface D'Alembert force density, and according to equation (3.1), we have
*Γ*_{kℓ},*k*,ℓ=*G*,*L*,*S*, and *k*≠ℓ, we have **n**_{k}=−**n**_{ℓ}. Subsequently, we may denote
*et al.* [16,17] have shown that the interface D'Alembert force can be written as
*k*ℓ).

Therefore, the dynamic equations of interface MCL read as

Without considering surface diffusion and friction, the following surface constitutive relations are chosen
**u** is the three-dimensional surface displacement field; *μ*_{L} is the surface viscosity. Note that **I** is the unit second-order tensor in a three-dimensional Euclidean space, and **P** is the projection tensor defined as
**n** is the unit out-normal of the surface at the point of interest. It may be noted that **P** is a symmetric tensor. **d**_{s} is the surface or interface rate of deformation that is defined as
**v** is the three-dimensional velocity field. The last equality holds because the projection tensor **P** is a symmetric tensor.

In some part of the text, in order to emphasize the material properties of the manifold, we write it as **I**^{(k)}_{s}, *k*=*G*,*L*,*S* or LS, GS and GL, etc., in a manner that is self-evident. One may note that
*κ* is the average curvature of the surface.

In equations (3.11)–(3.13), *γ*_{G},*γ*_{L} and *γ*_{S} are the surface tensions in different phases, *Γ*_{S} is the solid surface strain energy, *ϵ*_{s} is the surface strain tensor; ∇_{s} is the surface gradient operator defined in equation (3.2); and the operator ⊗ is the standard notation for tensor product in tensor algebra. For the case of infinitesimal deformation, Gurtin & Murdoch [15] proposed the following quadratic form of the surface strain energy:
*γ*_{S} as well,
_{S} and *μ*_{S} are the surface Lame constants. Hence, the surface constitutive equation for the solid phase becomes
*ϵ*_{s} is determined by projecting the bulk strain onto the surface, i.e.
** ϵ** is the Eulerian–Almansi strain in the bulk, which is defined as

**F**is the deformation gradient of the bulk continuum.

### Remark 3.1

In general, the surface adhesive stress tensor
**t**_{k}]≠0. Consider a simple case of gaseous and liquid interface. Let
*κ*_{LS} is the mean curvature of the gas–liquid interface. This example shows that the Young–Laplace equation is consistent with the no-slip condition. Moreover, it reveals that there are two factors that contribute to the elimination of shear stress singularity: the Pauli-repulsion gap and the fluid membrane formulation. The Pauli-repulsion gap will allow general slip condition, and the Gurtin-type fluid membrane formulation may be viewed as a sharp interface model as the limit of phase-field approach to the MCL (e.g. [24–26]). The Cahn–Hilliard-type phase field model of MCLs is effective in removing singularity of shear stress distribution, but requiring the solution of an interphase diffusion zone. Comparing with the other phase field model of MCL, the MMCL liquid surface membrane theory may be more elegant in mathematics and more relevant in colloidal surface physics and chemistry.

### Remark 3.2

By using the Lagrangian description, we can obtain the relationship between the unit out-normal vector in the current and reference configurations as
**C** is the right Cauchy–Green tensor

### (b) The solution of the interface fluid membrane problem

We are now considering the general solution of the proposed MMCL theory. The MCL hydrodynamics is in fact a fluid–structure interaction problem. It has two standard solutions: (i) monolithic solution and (ii) iterative solution.

#### (i) Monolithic solution

For a monolithic solution, we must solve the following coupled different equations simultaneously:
*k* and ℓ are two adjacent phases, and *Γ*_{kℓ} is the interface.

In (3.27), the interface traction jump [**t**]=(**t**_{k}−**t**_{ℓ}) across the interface *Γ*_{kℓ} is solely dependent on the intermolecular adhesion force of two bulk phases derived in the previous section.

#### (ii) Iterative solution

Alternatively, we can use the iterative solution approach, which is often used in the solution of fluid–structure interaction problems (e.g. [28]). That is, we solve equations (3.27) and (3.29) alternately or iteratively. In this case, in addition to intermolecular adhesion, the surface traction **t**_{k} and **t**_{ℓ} will depend on other factors as well. For instance, when we solve the displacement and velocity fields in the *k* phase, we need the traction **t**_{k} on traction boundary. Based on equation (3.28), we can find **t**_{k} as
**t**_{k} and **t**_{ℓ} are unknown. In order to solve **t**_{k} by using an iterative solver, we employ the following interface traction approximation to evaluate **t**_{ℓ}, so that it only depends on the interface adhesive force, i.e.
*the unilateral interface traction approximation*.

If we adopt the no-slip boundary condition, the relative slip velocity vanishes, i.e. **v**_{k}−**v**_{ℓ}=0. Then from equation (3.28), we find that the surface traction on the boundary of the *k*-phase becomes
*k*th phase boundary so that

### Remark 3.3

It is both interesting and surprising that we have found from equations (3.30)–(3.32) that the no-slip condition may not be the ‘culprit’ of the shear stress singularity problem after all. In fact, it is the multiscale character of the moving contact interface that prevents the macroscale MCL theory from yielding the correct solution.

In the following, we shall show how to incorporate the effects of the interface equations of motion (3.28) and use them to calculate the surface traction forces. For illustration purpose, we consider a simple surface stress expression that only contains the surface tension contribution while neglecting all other contributions such as surface elasticity, viscosity, friction, etc. In this simplified case, we have
*Ω*_{L}=*Γ*_{LS}(L)∪*Γ*_{GL}(L),

We now consider the following surface finite-element interpolation displacement and virtual displacement fields:
*N*_{Snode} is the total number of surface element nodes, and it may be different in each surface *Γ*_{kl}; *N*_{Is}(**x**) are the surface finite-element shape functions, **d**_{Is} and **w**_{Is} are surface nodal displacement and the surface virtual displacement at surface node *I*_{s}, *I*_{s}=1,2,…,*N*_{node}. Note that the surface FEM nodes are a subset of the bulk FEM nodes, and the virtual displacements of the surface nodes are not independent from the virtual displacements of the bulk nodes. There is a connectivity map, say Map_{c}(*I*_{s}), to connect the two, i.e. Map_{c}(*I*_{s})=*I* or vice versa.

The weak form of the surface traction can be written as
**q**_{1} and **q**_{t} are defined in figure 5. The first terms of three surface forces are the external force acting on the MCL. If we use the Lagrangian-type finite-element method, the MCL can always be treated as part of the element boundary [11,29]. Therefore, for a given material point at the MCL, i.e. *N*_{Is}(**x**)=1, and the external resultant acting on the MCL is

## 4. Numerical examples

In this section, a few numerical examples that were conducted by using the MMCL method are reported to validate the proposed theory.

### (a) Validation of the surface tension effect

In the first example, we consider an ellipsoidal liquid droplet suspended in the air. In this case, the only driving force for the evolution of the system is the surface tension on the interface of the liquid and the atmosphere.

In the simulation, the surface tensions are applied based on equation (3.37). Constitutive equation of the droplet is given by equation (2.3). The material parameters used in the simulation are *κ*=2.2×10^{9} Pa, *μ*=0.6×10^{−3} Pa s, and the surface tension between the droplet and the atmosphere *γ*_{GL}=7.28×10^{−2} N m^{−1}. A total of 4341 particles is used in the discretization of the droplet of ellipsoidal shape. The radius of the first semi-axis of the ellipsoidal is *a*=5 nm, and the initial aspect ratio of three semi-principal axes are *a*:*b*:*c*=1:1.5:1. The simulation time step is chosen as d*t*=2.0×10^{−14} s. Time history of the ellipse changing to perfect sphere process are shown in figure 6. We plot the evolution of the lengths of longest and shortest semi-principal axes of the ellipsoidal droplet and the pressure change at the centre of the droplet in figure 7.

Since the relaxation process is a dynamic process, so it is expected that there are some oscillations in the length of the semi-principal axes at the beginning of the simulation. But eventually, both the principal axes and the pressure at the centre of the droplet will reach to the stable values, owing to the damping effect of the droplet viscosity. To check the validity of the final state, we would like to confirm whether the following Young–Laplace equation for a spherical droplet is satisfied or not:
*P* is the pressure change, *γ* is the surface tension of the liquid and *R* is the final radius of the droplet. For the final configuration of the droplet, the pressure change Δ*P*=25.62 MPa, the final radius of the spherical droplet *R*=5.71 nm. The surface tension *γ*=*γ*_{LG}=72.75×mN m^{−1}. Note that
*P*=25.62 MPa, which implies that the Young–Laplace equation is satisfied.

In fact, based on the energetic argument, the droplet tends to stay with the shape that has the lowest energy state. Thus, minimizing the surface area of an ellipsoidal shall make the ellipsoidal become a sphere. The simulation result confirms that the implementation of surface tension force agrees with the physics.

### (b) Liquid droplet spreading

In this example, simulations of the droplet spreading by using the MMCL are presented. The droplet is modelled as a perfect sphere with a radius *r*=5 nm, and the substrate is a short cylinder shape platform with the dimension *H*×*R*=4.5×20 nm, as shown in figure 8*a*. The liquid phase is the water droplet, and the solid substrate is treated as a single crystal copper. Constitutive equations of the liquid droplet and solid substrate are provided by equations (2.3) and (2.4). Initially, the droplet is placed at 5 nm above the top of the substrate. In finite-element discretization, 4341 nodes and 4000 elements are used to discretize the droplet, and 4036 nodes and 2916 elements are used to discretize the substrate, as shown in figure 8*b*. The material properties are given as follows. For the droplet (water): we choose the density *ρ*^{d}=1.0×10^{3} kg m^{−3}, viscosity *μ*=0.6×10^{−3} Pa s and bulk modulus *κ*=2.2×10^{9} Pa; for the substrate (copper): we have the density *ρ*^{s}=8.94×10^{3} kg m^{−3}, Young's modulus *E*=1.2×10^{11} Pa, and Poisson's ratio *ν*=0.34.

The parameters for the coarse-grained contact model are provided as the droplet atomic density *β*^{d}=3.33×10^{28}/m^{3}, the substrate atomic density *β*^{s}=8.47×10^{28}/m^{3}. The parameters for the Lennard-Jones potential for the intermolecular interaction between the droplet and the substrate are chosen as *σ*^{ds}=3.05×10^{−10} m, *ε*^{ds}=8.4725×10^{−21} J, which are obtained based on the Lorentz-Berthelot rule:

In this simulation, the time step is Δ*t*=5×10^{−15} s and the total steps is *n*steps=200 000. Reduced units are used throughout the simulation, with the unit of time *t*_{0}=1.0×10^{−12} s, length *l*_{0}=1.0×10^{−10} m and mass m_{0}=1.0×10^{−27} kg.

A time sequence of a dynamic droplet spreading is shown in figure 9. One can see that the droplet first adheres to the substrate (figure 9*a*–*b*) and then gradually spreads over it (figure 9*c*–*h*).

To validate the MMCL model, a comparison study of the dynamic contact angle evolutions of during droplet spreadings with MD simulations is performed. The surface tension between water and the atmosphere is chosen as *γ*_{S}=1.0398 N m^{−1}, surface Lame parameters λ^{S}=15.6 N m^{−1}, *μ*^{S}=−8.6 N m^{−1}, which are obtained from Choi *et al*. [31]. The only unknown parameter for the interface layer is *γ*_{L}, which can be obtained based on the target contact angle [11]. In this simulation, three sets of target contact angle are chosen, to illustrate the capability of the MMCL in modelling hydrophobic and hydrophilic behaviours in wetting phenomena. The dynamic contact angle is defined as the average of all tilted angles of the elements on the MCL. The time history of the dynamic angles obtained by using MMCL on different elastic substrates is displayed in figure 10, in which the molecular dynamics simulation results are taken from Blake *et al*. [32]. One can see that the dynamic contact angle of the MMCL are in general agreement with the result obtained by molecular dynamics simulations [32]. Meanwhile, the time dependence of the droplet contact radius *r* for the three different cases is shown in figure 11. One can see that the data fall into three different curves, corresponding to the three different contact angles. The curves separate after *t*=15 ps, which indicates that the early-time dynamics are independent of the wettability.

To test convergence of the method, three different meshes of the droplet model are considered: mesh 1 (4341 nodes, 4000 elements), mesh 2 (15 657 nodes, 13 892 elements) and mesh 3 (21 253 nodes, 19 300 elements). To reduce the computational cost, the mesh size of the substrate is fixed. In the convergence study, we set the equilibrium contact angle to be *θ*=50°. Figure 12 shows the droplet configurations for the three different meshes at time *t*=52 ps. It can be seen that the three configurations are almost the same. To provide a qualitative comparison during the entire spreading process, the evolution of the droplet contact radius for the droplet with respect to three different meshes versus time is plotted in figure 13. It can be seen that the three curves are almost overlapped, showing the convergence of the proposed model to some extent. We would like to mention that the detection for the locations of the MCLs is very important for the application of the surface tensions. To limit the length of the paper, the numerical detection algorithm is not included in this paper.

Moreover, in this section, a preliminary study of the droplet spreading by using the coarse-grained Lennard-Jones potential is reported, to demonstrate the capability of the proposed MMCL theory. Details of the coarse-graining techniques can be found in [14]. In this case, the coarse-graining parameter is set to be *η*=2.215, which indicates that the radius of the droplet studied here is *R*=1000*r*=5 μm. The material parameters of the droplet and the substrate are the same as the previous case, except that the equilibrium contact angle is set to be 0°. The contact radius of the droplet as a function of time and the log–log scale of the function are plotted in figure 14. One can see that at early stage the contact radius approximately grows proportionally to *t*^{1/2}, due to the capillary wave generation [33]. The linear fitting of the early stage contact radius shows the validity of the model to some degree. In the final stage, the contact radius is believed to be related with time as *r*∼*t*^{1/10} (not fitted in the figure), resulting from a balance between the surface tension and the viscous force close to the MCL [34,35]. Figure 15 shows the final configuration of the droplet. Detailed analysis of the large droplet spreading will be reported in a separate paper, in order to keep the presentation to an appropriate length.

### (c) Capillary motion on a spherical tip

In this example, in an attempt to capture the capillary motion of liquid climbing along solid surface, a rigid solid sphere is first slowly pushed into an infinite large water film, and then it is gently pulled back to the original position. In the end, the rigid sphere is held still at that position. The numerical model is shown in figure 16. The radius of the spherical tip is *r*=5 nm. The thickness of the water film is *H*=4.5 nm. Notice that the water film is rested on an infinite large substrate underneath. A total of 5812 nodes and 2904 elements are used to discretize the water film. For the computation, no mesh is needed for the spherical tip and the solid substrate, given the fact that they are both treated as rigid solids, i.e. the elastic deformation is neglected. For visualization purposes, a mesh for the spherical tip is included in the postprocess.

The interaction between the rigid sphere and the water film and the interplay between the water film and the infinite substrate are described by using the 12-6 Lennard-Jones potential. The corresponding surface stress tensors for the interaction of the water film with the rigid sphere and the infinite plate are given in equations (2.20), (2.21) and (2.18).

The constitutive equation and material properties of the water film are the same as the previous example. The simulation time step is Δ*t*=1.0×10^{−16} s, and the total simulation time steps is *n* steps=1 000 000. From *t*=0 to 0.025 ns, the rigid sphere is gradually pushed down at the speed of Δ*z*=2.0×10^{−14} m per time step (figure 17*a*–*c*). In time duration of *t*=0.025 ns to *t*=0.05 s, the sphere tip is pulled up with the same speed (figure 17*d*–*f*). One can find that during the push and pull process the water film motion follows the outer front surface of the rigid sphere tip and the substrate below the water film acts as a glue that fix the vertical displacement of the water film. After *t*=0.05 ns, the sphere is held still at that position until the end of the simulation (figure 17*g*–*j*). It can be observed that even after the sphere tip position is fixed, water still keeps climbing up along the sphere surface until an equilibrium state is reached, and this is the capillary motion.

## 5. Conclusion

In this work, we have developed an MMCL theory that couples a molecular adhesive contact theory with an interface hydrodynamics theory of MCLs, which is an extension of the Gurtin–Murdoch theory of elastic solid interface membrane to fluid interface membrane, to solve MCL problems. The MMCL theory has three novel technical features: (i) the MCL is established on an interface fluid membrane; (ii) an adoption of the Derjaguin approximation to calculate interface adhesive traction jump; and (iii) the proposal of unilateral adhesive traction approximation in the iterative solution of interface traction.

The reported validation and simulation results clearly show that the proposed method can be used to calculate general MCL problems by using continuum finite-element computation with adequate accuracy and minimum computational resources. So far we have used MMCL method to compute several scientific and engineering problems [36,37]. The largest computation is involved with a liquid crystal elastomer droplet with the radius of 5 μm and with 20 μs spreading time [36]. This computation had been performed on a desktop computer with a series computation code. This indicates that indeed MMCL method has potential to be applied in solving macroscale problems of colloidal mechanics and physics.

The most surprising finding of this work is that we have found that the no-slip boundary condition is not necessarily always causing singularity problem in shear stress distribution at the MCL front. If we view the fluid membrane model proposed here as a limit of a phase-field formulation, the no-slip condition may not cause shear stress singularity while simplifying the calculations. The real problem is that the conventional MCL theory stems from a macroscale phenomenological mechanics modelling, and it lacks the length scale to differentiate the coupling and interplay between the macroscale flow and the molecular adhesion in order to provide viable solution to describe the MCL motion.

On the other hand, the essence of the MMCL theory is to couple the molecular scale adhesive theory and surface tension formulation with the macroscale hydrodynamics theory to form an integrated multiscale description that is able to predict MCL evolution.

## Data accessibility

The data of this paper are available in the following repository site by free download: https://drive.google.com/open?id=0B2NUvmWMHW4XNHZVTFFsMVFteXc&authuser=0. The interested readers who have no internet access may also write to the corresponding author S.L. to request data.

## Author' contributions

S.L. proposed and formulated the theory and formulation; S.L. and H.F. formulated the computational algorithm; H.F. developed computer code; H.F. and S.L. performed numerical simulations; S.L. and H.F. analysed the computing data and simulation results, and S.L. and H.F. write the paper.

## Competing interests

The authors declare no competing interests in this work.

## Funding statement

H.F. was partially supported by a graduate fellowship from China Scholar Council (CSC). This support is greatly appreciated.

- Received April 4, 2015.
- Accepted May 28, 2015.

- © 2015 The Author(s) Published by the Royal Society. All rights reserved.