Search tips
Search criteria 


Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
Proc Math Phys Eng Sci. 2015 July 8; 471(2179): 20150224.
PMCID: PMC4528669

On multiscale moving contact line theory


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.

Keywords: adhesive contact, droplet spreading, surface tension

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. [16]).

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. [710]). 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 [1214] 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 [1518]. 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.


where Γ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 [partial differential]Ωα, α=G,L,S and [partial differential]ΩG=ΓGL(G)[union or logical sum]ΓGS(G), [partial differential]ΩL=ΓGL(L)[union or logical sum]ΓLS(L) and [partial differential]ΩG=ΓLS(S)[union or logical sum]ΓGS(S).

Figure 1.
The triple phase system of gaseous (G), liquid (L) and solid (S). (Online version in colour.)

(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,


where k denotes the phase; σk is the Cauchy stress; ρk is the mass density per unit volume; bk is the body force per unit mass and uk 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


where κ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:


where S is the second Piola–Kirchhoff stress; E=12(FTFI) is the Green–Lagrangian strain; 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 [2022], 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 ϕ(|xk[ell]|), where xk[ell]=x[ell]xk,k,[ell]=G,L,S, k[ell] and xk[set membership]Ωk,x[ell][set membership]Ω[ell]. Then the total interphase adhesive-contact potential energy Πac for two interacting phase is

Πac=ΩkΩlβkβlϕ(r)dΩldΩk,r=| xkx|,

where ϕ(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


where ε 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

δΠac=ΩkΩlβkβ(ϕ(r)xkδuk+ϕ(r)xδu)dΩdΩk=Ωkβkbkδ ukdΩkΩlβl bδudΩ,

where we define the adhesive body force in each phase as




where Φk and Φ[ell] are the homogenized macro interphase interaction potentials between phases k and l, k,[ell]=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 [1214].

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[ell]=Γk[ell](k)[union or logical sum]Γk[ell]([ell]) between the phases k and l and choose two arbitrary points xk[set membership][partial differential]Ωk and x[ell][set membership][partial differential]Ω[ell] so that sk[ell]:=x[ell]xk=:−s[ell]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 dak[subset or is implied by][partial differential]Ωk due to the presence of an infinitesimal surface element surface da[ell][subset or is implied by][partial differential]Ω[ell] is expressed as


where nk and n[ell] are the unit surface out-normal at points xk and x[ell], respectively, k,[ell]=G,L,S, k[ell], and


According to Jagota & Argento [23], the only requirement for such conversion is that the interaction potential ϕ(r) decays faster than 1/r3, 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,


Similar to equation (2.10), the adhesive contact force applied on an infinitesimal surface element da[ell][subset or is implied by][partial differential]Ω[ell] due to the presence of an infinitesimal surface element surface dak[subset or is implied by][partial differential]Ωk is expressed as


It may be noted that dFk≠−dF[ell], but the total force applied on phase k due to the presence of phase [ell] is equal to that of the force applied on phase [ell] by the phase k, i.e. i=1NkdFk=i=1NldF. Nk and N[ell] are the total number of infinitesimal surface elements on [partial differential]Ωk and [partial differential]Ω[ell], respectively.

Figure 2.
Derjaguin approximation: a schematic to convert the body–body interaction to the surface-surface interaction. (Online version in colour.)

Furthermore, one can rewrite equation (2.10) as follows:


which indicates that the quantity in the curly bracket of the equation above can be viewed as an infinitesimal surface stress contribution from da[ell][subset or is implied by][partial differential]Ω[ell] at xk[set membership][partial differential]Ωk. Thus, if we perform the integration over the surface [partial differential]Ω[ell], we may obtain the adhesive surface stress


which represents the surface stress tensor at point xk[set membership][partial differential]Ωk, due to the interaction from the phase [ell]. The corresponding surface traction at point xk[set membership][partial differential]Ωk can then be expressed as

tkadh=σkadh nk,

Similarly, we can derive the surface stress tensor and traction at point x[ell][set membership]Γk[ell]([ell]) 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 3a, then one can obtain an explicit expression for surface stress tensor as


where a is the shortest distance of the particle xk to the plate, ez is the unit vector along the z-axis. Alternatively, if the phase [ell] is a rigid sphere with radius R as shown in figure 3b, the surface stress tensor can be found as a diagonal tensor






In the equations above, a is the shortest distance of the surface particle xk to the centre of the rigid sphere. ex, ey and ez 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.

Figure 3.
The surface–surface adhesion integrations: (a) a rigid half plane with an arbitrary body and (b) a rigid sphere with an arbitrary body. (Online version in colour.)

(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 4b 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 4b). 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 4a), 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.

Figure 4.
Comparison of shear stress (σ12) distributions of the front of the MCL: (a) MCL and (b) MMCL. (Online version in colour.)

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. [1518]), we consider the following surface equations of motion for the surface of all three phases:

fD,k:=sςk+ tk,xΩk,k=G,L,S,

where fD,k is the surface D'Alembert force density, ςk is the surface stress, tk is the adhesive surface traction vector and [nabla]S is the surface gradient operator that is defined as


in which 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 tk with the coarse-grain adhesive molecular force. In fact, the adhesive traction force is not the only contribution to the surface traction tk,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

fD,G=sςG+ tG,xΩG,fD,L=sςL+ tL,xΩLandfD,S=sςS+ tS,xΩS,

which allows us to introduce the following interface D'Alembert force densities:






In general, on a specific interface Γk[ell],k,[ell]=G,L,S, and k[ell], we have nk=−n[ell]. Subsequently, we may denote



Moreover, Gurtin et al. [16,17] have shown that the interface D'Alembert force can be written as

fkD=ρ¯kvkintf( vvk)orfkD=ρ¯kvintf(vk v),

where vintf is the magnitude of the interface normal velocity (the superscript ‘intf’ means interface), and ρ¯k is the average mass density on the interface (k[ell]).

Therefore, the dynamic equations of interface MCL read as

sςk+[tk]=ρ¯kvkintf(v vk).

Without considering surface diffusion and friction, the following surface constitutive relations are chosen




where u is the three-dimensional surface displacement field; μL is the surface viscosity. Note that Is(2) denotes the unit tensor on a smooth surface or two-dimensional manifold, which is defined as


where I is the unit second-order tensor in a three-dimensional Euclidean space, and P is the projection tensor defined as


where n is the unit out-normal of the surface at the point of interest. It may be noted that P is a symmetric tensor. ds is the surface or interface rate of deformation that is defined as

ds:=PSym(v) PT=PSym(v)P,

where 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


where κ 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; [nabla]s is the surface gradient operator defined in equation (3.2); and the operator [multiply sign in circle] 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:


in which the surface elastic tensor is related to surface tension γS as well,


where λS and μS are the surface Lame constants. Hence, the surface constitutive equation for the solid phase becomes


Subsequently, one can readily derive the interface constitutive relation, for instance,


where the surface strain ϵs is determined by projecting the bulk strain onto the surface, i.e.


where ϵ is the Eulerian–Almansi strain in the bulk, which is defined as


where F is the deformation gradient of the bulk continuum.

Remark 3.1 —

In general, the surface adhesive stress tensor


thus [tk]≠0. Consider a simple case of gaseous and liquid interface. Let


The interface MCL equation yields


If we assume that the no-slip boundary condition may hold in this case, and it then leads to


which is the classical Young–Laplace equation of gas–liquid interface, and κ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. [2426]). 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


where C is the right Cauchy–Green tensor


The mean curvature is related to the divergence of the surface unit out-normal (e.g. [27]) as


By substituting equation (3.22) into equation (3.24), the following equation can be obtained:


where G is the derivative of the deformation gradient tensor, and it is defined as


(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:


sςk+(tk t)=ρ¯vkinft(v vk);xΓk

andσ+ρb=ρu¨and t=σn,xΩt,

in which k and [ell] are two adjacent phases, and Γk[ell] is the interface.

In (3.27), the interface traction jump [t]=(tkt[ell]) across the interface Γk[ell] 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 tk and t[ell] 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 tk on traction boundary. Based on equation (3.28), we can find tk as

tk=ρ¯vkintf(v vk)sςk t.

In the above equation, both traction forces tk and t[ell] are unknown. In order to solve tk by using an iterative solver, we employ the following interface traction approximation to evaluate t[ell], so that it only depends on the interface adhesive force, i.e.


We call this procedure the unilateral interface traction approximation.

If we adopt the no-slip boundary condition, the relative slip velocity vanishes, i.e. vkv[ell]=0. Then from equation (3.28), we find that the surface traction on the boundary of the k-phase becomes


Similarly, when we are solving the displacement and velocity fields for the bulk phase [ell], we can obtain the traction on the boundary of the [ell]-phase by first applying the unilateral adhesive traction approximation on the traction on the kth phase boundary so that


and with the help of the no-slip boundary condition, we can obtain


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

sςk+[tk]=ρ¯kvkinf(vk v),whereςkl=γkIs,xΓk,k=GL, LS, GS,

Then, based on the no-slip condition and the unilateral adhesive approximation, we can find the traction force on the liquid surface [partial differential]ΩL=ΓLS(L)[union or logical sum]ΓGL(L),




One may note that at the MCL front, three interphases merge together, and the traction on the gaseous phase should be added to the MCL as well. Therefore, we also consider


where σSadhnG is the solid surface traction force due to gaseous phase interaction.

We now consider the following surface finite-element interpolation displacement and virtual displacement fields:


where NSnode is the total number of surface element nodes, and it may be different in each surface Γkl; NIs(x) are the surface finite-element shape functions, dIs and wIs are surface nodal displacement and the surface virtual displacement at surface node Is, Is=1,2,…,Nnode. 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 Mapc(Is), to connect the two, i.e. Mapc(Is)=I or vice versa.

The weak form of the surface traction can be written as

Is=1NGLwIs{ΓGL(L)NIs(x)sςGLds+ΓGL(L)NIs(x)σGadh nLds},

Is=1NLSwIs{ΓLS(L)NIs(x)sςLSds+ΓLS(L)NIs(x)σSadh nLds}

andIs=1NGSwIs{ΓGS(G)NIs(x)sςGSds+ΓGS(G)NIs(x)σSadh nGds}.

Integration by parts yields

Is=1NGL{NIs(x)ςGLqt|ΓLS(L)+ΓGL(L)NIsx:ςGLds+ΓGL(L)NIs(x)σGadhnLds},Is=1NLS{NIs(x)ςLSq1|ΓLS(L)+ΓLS(L)NIsx:ςLSds+ΓLS(L)NIs(x)σSadh nLds}andIs=1NGS{NIs(x)ςLS(q1)|ΓGS(G)+ΓGS(G)NIsx:ςLSds+ΓLS(L)NIs(x)σSadhnLds},

where vectors q1 and qt 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. xLmcl=ΓGLΓLSΓGS, we have NIs(x)=1, and the external resultant acting on the MCL is




as shown in figure 5. If we only consider the surface tension contribution to surface stress, we then have


Thus for the MCL shown in figure 5, the external force acting on the MCL is



Under the equilibrium condition, which is implied in the no-slip condition, the right-hand side of first equation of (3.41) vanishes, and one can obtain,


which is exactly the classical Young's equation [30].

Figure 5.
MCL front. (Online version in colour.)

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×109 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 dt=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.

Figure 6.
Time history of the simulation of a three-dimensional ellipsoidal droplet embedded in atmosphere, driven by the surface tension effect. (a) t=0.004 ns, (b) t=0.04 ns, (c) t=0.08 ns, (d) t= 0.2 ns, (e) t=0.28 ns ...
Figure 7.
Evolution history of the ellipsoidal droplet: (a) length of principal axes and (b) pressure at the centre of the droplet. (Online version in colour.)

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:


where Δ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


and we end up with an error percentage of 0.54% for the pressure change compared to the value Δ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 8a. 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 8b. The material properties are given as follows. For the droplet (water): we choose the density ρd=1.0×103 kg m−3, viscosity μ=0.6×10−3 Pa s and bulk modulus κ=2.2×109 Pa; for the substrate (copper): we have the density ρs=8.94×103 kg m−3, Young's modulus E=1.2×1011 Pa, and Poisson's ratio ν=0.34.

Figure 8.
Computational model of the droplet wetting. (a) The MMCL computational model and (b) FEM mesh of the MMCL model. (Online version in colour.)

The parameters for the coarse-grained contact model are provided as the droplet atomic density βd=3.33×1028/m3, the substrate atomic density βs=8.47×1028/m3. 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: σds=0.5(σd+σs),εds=εdεs

In this simulation, the time step is Δt=5×10−15 s and the total steps is nsteps=200 000. Reduced units are used throughout the simulation, with the unit of time t0=1.0×10−12 s, length l0=1.0×10−10 m and mass m0=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 9ab) and then gradually spreads over it (figure 9ch).

Figure 9.
Time sequence of the droplet spreading on an elastic substrate. The dots (yellow online) at the fringe denote the current positions of the MCLs. (a) t=0.010 ns, (b) t=0.020 ns, (c) t=0.030 ns, (d) t=0.040 ns, (e) t=0.070 ns, ...

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 γGL=7.28×102Nm1. The surface material parameters of the substrate are γ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.

Figure 10.
The dynamic contact angle evolution with different substrates. (Online version in colour.)
Figure 11.
Time evolution of the droplet contact radius for the three different contact angles. (Online version in colour.)

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.

Figure 12.
Droplet configurations for the three different meshes at t=52 ps. (a) mesh 1, (b) mesh 2 and (c) mesh 3. (Online version in colour.)
Figure 13.
Time evolution of the droplet contact radius for three different meshes. (Online version in colour.)

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=1000r=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 t1/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~t1/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.

Figure 14.
Evolution of the contact radius for the large droplet with radius R= 5 μm. (a) Contact radius versus the simulation time (b) log–log plot of the contact radius evolution curve. t0=1 μs; R0=1 μm. ...
Figure 15.
Final configuration of the large droplet spreading. (Online version in colour.)

(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.

Figure 16.
The model for capillary motion on a spherical tip. (Online version in colour.)

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 17ac). In time duration of t=0.025 ns to t=0.05 s, the sphere tip is pulled up with the same speed (figure 17df). 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 17gj). 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.

Figure 17.
Time sequence of the capillary motion on a spherical tip. (a) t=0.005 ns, (b) t=0.015 ns, (c) t=0.025 ns, (d) t= 0.030 ns, (e) t=0.040 ns, (f) t=0.050 ns, (g) t= 0.055 ns, (h) t=0.060 ns, ...

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: 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.


1. Huh C, Scriven LE. 1971. Hydrodynamic model of steady movement of a solid/liquid/fluid contact line. J. Colloid Interface Sci. 35, 85–101. (doi:10.1016/0021-9797(71)90188-3)
2. Dussan V EB, Davis SH. 1974. The moving contact line: the slip boundary condition. J. Fluid Mech. 65, 71–95. (doi:10.1017/S0022112074001261)
3. Dussan V EB. 1976. The moving contact line: the slip boundary condition. J. Fluid Mech. 77, 665–684. (doi:10.1017/S0022112076002838)
4. Shikhmurzaev YD. 2006. Singularities at the moving contact line. Mathematical, physical and computational aspects. Physica D 217, 121–133. (doi:10.1016/j.physd.2006.03.003)
5. Qian T, Wang X-P, Sheng P. 2006. A variational approach to moving contact line hydrodynamics. J. Fluid Mech. 564, 333–360. (doi:10.1017/S0022112006001935)
6. Wang C, Zhang LT. 2015. Numerical modeling of gas-liquid-solid interactions: gas–liquid free surfaces interacting with deformable solids. Comput. Methods Appl. Mech. Eng. 286, 123–146. (doi:10.1016/j.cma.2014.12.011)
7. Qian T, Wang X-P, Sheng P. 2003. Molecular scale contact line hydrodynamics of immiscible flows. Phys. Rev. E 68, 016306 (doi:10.1103/PhysRevE.68.016306) [PubMed]
8. Qian T, Wang X-P, Sheng P. 2006. Molecular hydrodynamics of the moving contact line in two-phase immiscible flows. Commun. Comput. Phys. 1, 1–52.
9. Wang X, Peng X, Lee D. 2003. Dynamic wetting and stress singularity on contact line. Science China 46, 407–417. (doi:10.1360/02ye0407)
10. Hadjiconstantinou NG. 1999. Hybrid atomisticcontinuum formulations and the moving contact-line problem. Comput. Phys. 154, 245–265. (doi:10.1006/jcph.1999.6302)
11. Minaki H, Li S. 2014. Multiscale modeling and simulation of dynamic wetting. Comput. Methods Appl. Mech. Eng. 273, 274–302.
12. Sauer RA, Li S. 2007. A contact mechanics model for quasi-continua. Int. J. Numer. Methods Eng. 71, 931–962. (doi:10.1002/nme.1970)
13. Sauer RA, Li S. 2007. An atomic interaction-based continuum model for adhesive contact mechanics. Finite Elements Anal. Des. 43, 384–396. (doi:10.1016/j.finel.2006.11.009)
14. Sauer RA, Li S. 2007. An atomiscically enriched continuum model for nanoscale contact mechanics and its application to contact scaling. J. Nanosci. Nanotechnol. 8, 1–17. [PubMed]
15. Gurtin ME, Murdoch AI. 1975. A continuum theory of elastic material surfaces. Arch. Rational Mech. Anal. 57, 291–323. (doi:10.1007/BF00261375)
16. Gurtin ME, Struther A. 1990. Multiphase thermomechanics with interfacial structure 3. evolving phase boundaries in the presence of bulk deformation. Arch. Rational Mech.Anal. 112, 97–160. (doi:10.1007/BF00375667)
17. Gurtin ME. 1991. Dynamics of evolving phase boundaries in deformable continua. Le Matematiche XLVI, Fasc. I, 195–201.
18. Gurtin ME, Weissmuller J, Larche F. 1998. A general theory of curved deformable interfaces in solids at equilibrium. Philos. Mag. A 78, 1093–1109. (doi:10.1080/01418619808239977)
19. Bradley RS. 1932. The cohesive force between solid surfaces and the surface energy of solids. Philos. Mag. Ser. 7, 853–862.
20. Derjaguin B, Landau L. 1941. Theory of the stability of strongly charged lyophobic sols and of the adhesion of strongly charged particles in solutions of electrolytes. Acta Physico Chemica URSS 14, 633–662.
21. Verwey EJW, Overbeek JThG. 1948. Theory of the stability of lyophobic colloids. Amsterdam, The Netherlands: Elsevier Publishing Company.
22. Derjaguin BV. 1934. Analysis of friction and adhesion, IV. The theory of the adhesion of small particles. Kolloid Z. 69, 155–164. (doi:10.1007/BF01433225)
23. Jagota A, Argento C. 1997. An intersurface stress tensor. J. Colloid Interface Sci. 191, 326–336. (doi:10.1006/jcis.1997.4933) [PubMed]
24. Jacqmin D. 2000. Contact-line dynamics of a diffuse fluid interface. J. Fluid Mech. 402, 57–88. (doi:10.1017/S0022112099006874)
25. Yue P, Zhou C, Feng JJ. 2010. Sharp-interface limit of the Cahn–Hilliard model for moving contact lines. J. Fluid Mech. 645, 279–294. (doi:10.1017/S0022112009992679)
26. Wang XP, Wang YG. 2007. The sharp interface limit of a phase field model for moving contact line problem. Methods Appl. Anal. 14, 287–294. (doi:10.4310/MAA.2007.v14.n3.a6)
27. Peters JMH. 2001. Total curvature of surfaces (via the divergence of the normal. Int. J. Math. Educ. Sci. Technol. 32, 795–810. (doi:10.1080/00207390110053766)
28. Kuttler U, Wall WA. 2008. Fixed-point fluid–structure interaction solvers with dynamic relaxation. Comput. Mech. 43, 61–72. (doi:10.1007/s00466-008-0255-5)
29. Fan H, Li S. In preparation Computational formulations of multiscale moving contact line and simulation of droplet spreading.
30. Young T. 1805. An essay on the cohesion of fluids. Phil. Trans. R. Soc. 95, 65–87. (doi:10.1098/rstl.1805.0005)
31. Choi J, Cho M, Kim W. 2010. Multiscale analysis of nanoscale thin film considering surface effects: thermomechanical properties. J. Mech. Mater. Struct. 5, 161–183. (doi:10.2140/jomms.2010.5.161)
32. Blake TD, Clarke A, de Coninck J, de Ruijter M, Voue M. 1999. Droplet spreading: a microscopic approach. Colloids Surf. A Physicochem. Eng. Aspects 149, 123–130. (doi:10.1016/S0927-7757(98)00602-5)
33. Bird JC, Mandre S, Stone HA. 2008. Short-time dynamics of partial wetting. Phys. Rev. Lett. 100, 234501 (doi:10.1103/PhysRevLett.100.234501) [PubMed]
34. Tanner LH. 1979. The spreading of silicone oil drops on horizontal surfaces. J. Phys. D Appl. Phys. 12, 1473 (doi:10.1088/0022-3727/12/9/009)
35. Rafai S, Sarker D, Bergeron V, Meunier J, Bonn D. 2002. Superspreading: aqueous surfactant drops spreading on hydrophobic surfaces. Langmuir 18, 10 486–10 488. (doi:10.1021/la020271i)
36. Fan H, Li S. 2015. Modeling universal dynamics of cell spreading on elastic substrates. Biomech. Model. Mechanobiol. (BMMB) 14, 1–16. (doi:10.1007/s10237-015-0673-1) [PubMed]
37. Fan H, Li S. 2014. Modeling microtubule cytoskeleton via an active liquid crystal elastomer model. Comput. Mater. Sci. 96 Part B, 559–566.

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society