Home | About | Journals | Submit | Contact Us | Français |

**|**Proc Math Phys Eng Sci**|**PMC4528669

Formats

Article sections

- Abstract
- 1. Introduction
- 2. The basic idea of multiscale moving contact line and how it works
- 3. Dynamic interface moving contact line theory
- 4. Numerical examples
- 5. Conclusion
- References

Authors

Related links

Proc Math Phys Eng Sci. 2015 July 8; 471(2179): 20150224.

PMCID: PMC4528669

Department of Civil and Environmental Engineering, University of California, Berkeley, CA 94720, USA

e-mail: ude.yelekreb@nafoahs

Received 2015 April 4; Accepted 2015 May 28.

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

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.

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.

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.

$$\begin{array}{rl}{\mathit{\Gamma}}_{\mathrm{LS}}& ={\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})\cup {\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{S}),\phantom{\rule{1em}{0ex}}{\mathit{\Gamma}}_{\mathrm{GL}}={\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{G})\cup {\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{L})\\ \mathrm{and}\phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}{\mathit{\Gamma}}_{\mathrm{GS}}& ={\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{G})\cup {\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{S}),\end{array}\}$$

2.1

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 *Ω*_{α}, *α*=*G*,*L*,*S* and *Ω*_{G}=*Γ*_{GL}(G)*Γ*_{GS}(G), *Ω*_{L}=*Γ*_{GL}(L)*Γ*_{LS}(L) and *Ω*_{G}=*Γ*_{LS}(S)*Γ*_{GS}(S).

The MMCL is a hybrid continuum model. The Lagrange description of continuum mechanics is used to model the dynamic deformation of each phase,

$$\mathrm{\nabla}\cdot {\mathit{\sigma}}_{k}+{\rho}_{k}{\mathbf{\text{b}}}_{k}={\rho}_{k}{\ddot{\mathbf{\text{u}}}}_{k},\phantom{\rule{1em}{0ex}}k=G,L\hspace{0.17em}\mathrm{and}\hspace{0.17em}S,$$

2.2

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

$$\mathit{\sigma}={\kappa}_{\mathrm{L}}(\mathrm{ln}J)\mathbf{\text{I}}+{\mu}_{\mathrm{L}}(\mathrm{\nabla}\otimes \mathbf{\text{v}}+{(\mathrm{\nabla}\otimes \mathbf{\text{v}})}^{\mathrm{T}}),$$

2.3

where *κ*_{L} and *μ*_{L} are, respectively, the bulk modulus and viscosity of the liquid phase; ** v** is the velocity of the liquid phase,

The solid phase is modelled as a St. Venant-Kirchhoff material, and its constitutive equation is given as follows:

$$\mathbf{\text{S}}={\mathrm{\lambda}}_{\mathrm{S}}\hspace{0.17em}\mathrm{tr}(\mathbf{\text{E}})\mathbf{\text{I}}+2{\mu}_{\mathrm{S}}\mathbf{\text{E}},$$

2.4

where **S** is the second Piola–Kirchhoff stress; $\mathbf{\text{E}}=\frac{1}{2}({\mathbf{\text{F}}}^{\mathrm{T}}\mathbf{\text{F}}-\mathbf{\text{I}})$ is the Green–Lagrangian strain; **F** is the deformation gradient and λ_{S}, *μ*_{S} are the Lame constants of the solid substrate.

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

$${\mathit{\Pi}}^{\mathrm{ac}}={\int}_{{\mathit{\Omega}}_{k}}{\int}_{{\mathit{\Omega}}_{l}}{\beta}_{k}{\beta}_{l}\varphi (r)\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{l}\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{k},\phantom{\rule{1em}{0ex}}r=|{\mathbf{\text{x}}}_{k}-{\mathbf{\text{x}}}_{\ell}|,$$

2.5

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

$$\varphi (r)=\epsilon [{\left(\frac{{\sigma}_{0}}{r}\right)}^{12}-2{\left(\frac{{\sigma}_{0}}{r}\right)}^{6}],$$

2.6

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

$$\begin{array}{rl}\delta {\mathit{\Pi}}^{\mathrm{ac}}& ={\int}_{{\mathit{\Omega}}_{k}}{\int}_{{\mathit{\Omega}}_{l}}{\beta}_{k}{\beta}_{\ell}(\frac{\mathrm{\partial}\varphi (r)}{\mathrm{\partial}{\mathbf{\text{x}}}_{k}}\cdot \delta {\mathbf{\text{u}}}_{k}+\frac{\mathrm{\partial}\varphi (r)}{\mathrm{\partial}{\mathbf{\text{x}}}_{\ell}}\cdot \delta {\mathbf{\text{u}}}_{\ell})\mathrm{d}{\mathit{\Omega}}_{\ell}\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{k}\\ & =-{\int}_{{\mathit{\Omega}}_{k}}{\beta}_{k}{\mathbf{\text{b}}}_{k}\cdot \delta {\mathbf{\text{u}}}_{k}\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{k}-{\int}_{{\mathit{\Omega}}_{l}}{\beta}_{l}{\mathbf{\text{b}}}_{\ell}\cdot \delta {\mathbf{\text{u}}}_{\ell}\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{\ell}\hspace{0.17em},\end{array}$$

2.7

where we define the adhesive body force in each phase as

$${\mathbf{\text{b}}}_{k}({\mathbf{\text{x}}}_{k}):=-\frac{\mathrm{\partial}{\mathit{\Phi}}_{\ell}}{\mathrm{\partial}{\mathbf{\text{x}}}_{k}},\phantom{\rule{1em}{0ex}}{\mathit{\Phi}}_{\ell}:={\int}_{{\mathit{\Omega}}_{l}}{\beta}_{\ell}\varphi (r)\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{l}$$

2.8

and

$${\mathbf{\text{b}}}_{\ell}({\mathbf{\text{x}}}_{k}):=-\frac{\mathrm{\partial}{\mathit{\Phi}}_{k}}{\mathrm{\partial}{\mathbf{\text{x}}}_{\ell}},\phantom{\rule{1em}{0ex}}{\mathit{\Phi}}_{k}:={\int}_{{\mathit{\Omega}}_{k}}{\beta}_{k}\varphi (r)\hspace{0.17em}\mathrm{d}{\mathit{\Omega}}_{k},$$

2.9

where *Φ*_{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.

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

$$\mathrm{d}{\mathbf{F}}_{k}=\{{\beta}_{k}{\beta}_{\ell}({\mathbf{n}}_{l}\otimes {\mathbf{s}}_{kl})\cdot {\mathbf{n}}_{k}\psi (s)\}\hspace{0.17em}\mathrm{d}{a}_{k}\hspace{0.17em}\mathrm{d}{a}_{\ell},$$

2.10

where **n**_{k} and **n**_{} are the unit surface out-normal at points **x**_{k} and **x**_{}, respectively, *k*,=*G*,*L*,*S*, *k*≠, and

$$\psi (s)=\frac{1}{{s}^{3}}{\int}_{s}^{\mathrm{\infty}}\varphi (r){r}^{2}\hspace{0.17em}\mathrm{d}r,\phantom{\rule{1em}{0ex}}\mathrm{with}\hspace{0.17em}s=|{\mathbf{\text{s}}}_{k\ell}|,$$

2.11

According to Jagota & Argento [23], the only requirement for such conversion is that the interaction potential *ϕ*(*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,

$$\psi (s)=\frac{2}{3}\epsilon [\frac{1}{6}{\left(\frac{{\sigma}_{0}}{s}\right)}^{12}-{\left(\frac{{\sigma}_{0}}{s}\right)}^{6}].$$

2.12

Similar to equation (2.10), the adhesive contact force applied on an infinitesimal surface element d*a*_{}*Ω*_{} due to the presence of an infinitesimal surface element surface d*a*_{k}*Ω*_{k} is expressed as

$$\mathrm{d}{\mathbf{F}}_{\ell}=\{{\beta}_{\ell}{\beta}_{k}({\mathbf{n}}_{k}\otimes {\mathbf{s}}_{\ell k})\cdot {\mathbf{n}}_{\ell}\psi (s)\}\hspace{0.17em}\mathrm{d}{a}_{\ell}\hspace{0.17em}\mathrm{d}{a}_{k}.$$

2.13

It may be noted that d**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. $\sum _{i=1}^{{N}_{k}}\mathrm{d}{\mathbf{F}}_{k}=\sum _{i=1}^{{N}_{l}}\mathrm{d}{\mathbf{F}}_{\ell}$. *N*_{k} and *N*_{} are the total number of infinitesimal surface elements on *Ω*_{k} and *Ω*_{}, respectively.

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:

$$\frac{\mathrm{d}{\mathbf{F}}_{k}}{\mathrm{d}{a}_{k}}=\{{\beta}_{k}{\beta}_{\ell}({\mathbf{n}}_{\ell}\otimes {\mathbf{s}}_{k\ell})\psi (s)\hspace{0.17em}\mathrm{d}{a}_{\ell}\}\cdot {\mathbf{n}}_{k},$$

2.14

which indicates that the quantity in the curly bracket of the equation above can be viewed as an infinitesimal surface stress contribution from d*a*_{}*Ω*_{} at **x**_{k}*Ω*_{k}. Thus, if we perform the integration over the surface *Ω*_{}, we may obtain the adhesive surface stress

$${\mathit{\sigma}}_{k}^{\mathrm{adh}}={\int}_{\mathrm{\partial}{\mathit{\Omega}}_{\ell}}{\beta}_{k}{\beta}_{\ell}({\mathbf{n}}_{\ell}\otimes {\mathbf{s}}_{kl})\psi (s)\hspace{0.17em}\mathrm{d}{a}_{\ell},$$

2.15

which represents the surface stress tensor at point **x**_{k}*Ω*_{k}, due to the interaction from the phase . The corresponding surface traction at point **x**_{k}*Ω*_{k} can then be expressed as

$${\mathbf{\text{t}}}_{k}^{\mathrm{adh}}={\mathit{\sigma}}_{k}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{k},$$

2.16

Similarly, we can derive the surface stress tensor and traction at point **x**_{}*Γ*_{k}() as

$${\mathit{\sigma}}_{\ell}^{\mathrm{adh}}={\int}_{\mathrm{\partial}{\mathit{\Omega}}_{k}}{\beta}_{\ell}{\beta}_{k}({\mathbf{n}}_{k}\otimes {\mathbf{s}}_{\ell k})\psi (s)\hspace{0.17em}\mathrm{d}{a}_{k},\phantom{\rule{1em}{0ex}}{\mathbf{t}}_{\ell}^{\mathrm{adh}}={\mathit{\sigma}}_{\ell}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\ell}.$$

2.17

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

$$\begin{array}{rl}{\mathit{\sigma}}_{k}^{\mathrm{adh}}& =-\{2\pi a{\beta}_{k}{\beta}_{\ell}{\int}_{a}^{+\mathrm{\infty}}\psi (s)s\hspace{0.17em}\mathrm{d}s\}{\mathbf{e}}_{z}\otimes {\mathbf{e}}_{z}\\ & =-\left\{\frac{\pi \epsilon {\beta}_{k}{\beta}_{\ell}{a}^{3}}{45}[{\left(\frac{{\sigma}_{0}}{a}\right)}^{12}-15{\left(\frac{{\sigma}_{0}}{a}\right)}^{6}]\right\}{\mathbf{e}}_{z}\otimes {\mathbf{e}}_{z},\end{array}$$

2.18

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

$${\mathit{\sigma}}_{k}^{\mathrm{adh}}={\sigma}_{kxx}^{\mathrm{adh}}{\mathbf{e}}_{x}\otimes {\mathbf{e}}_{x}+{\sigma}_{kyy}^{\mathrm{adh}}{\mathbf{e}}_{y}\otimes {\mathbf{e}}_{y}+{\sigma}_{kzz}^{\mathrm{adh}}{\mathbf{e}}_{z}\otimes {\mathbf{e}}_{z},$$

2.19

where

$$\begin{array}{rl}{\sigma}_{kxx}^{\mathrm{adh}}& ={\sigma}_{kyy}^{\mathrm{adh}}=\pi {\beta}_{k}{\beta}_{l}\epsilon {\sigma}_{0}^{6}\{\frac{4(5{a}^{4}+14{a}^{2}{R}^{2}+5{R}^{4}){R}^{3}}{135{(a-R)}^{8}{(a+R)}^{8}}{\sigma}_{0}^{6}+\frac{({a}^{2}+{R}^{2})R}{3{a}^{2}{({a}^{2}-{R}^{2})}^{2}}-\frac{1}{12{a}^{3}}\mathrm{log}\left[{\left(\frac{a-R}{a+R}\right)}^{2}\right]\}\end{array}$$

2.20

and

$$\begin{array}{rl}{\sigma}_{kzz}^{\mathrm{adh}}& =\pi {\beta}_{k}{\beta}_{l}\epsilon {\sigma}_{0}^{6}\{-\frac{4{R}^{3}(55{a}^{6}+207{a}^{4}{R}^{2}+117{a}^{2}{R}^{4}+5{R}^{6})}{135{(a-R)}^{9}{(a+R)}^{9}}{\sigma}_{0}^{6}\\ & \phantom{\rule{1em}{0ex}}+\frac{2({a}^{4}+4{a}^{2}{R}^{2}-{R}^{4})R}{3{a}^{2}{(a-R)}^{3}{(a+R)}^{3}}+\frac{1}{6{a}^{3}}\mathrm{log}\left[\frac{{(a-R)}^{2}}{{(a+R)}^{2}}\right]\}.\end{array}$$

2.21

In the equations above, *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.

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.

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

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:

$${\mathbf{\text{f}}}^{D,k}:={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k}+{\mathbf{\text{t}}}_{k},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{k},\phantom{\rule{1em}{0ex}}k=G,L,S,$$

3.1

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

$${\mathrm{\nabla}}_{\mathrm{s}}:=\mathrm{\nabla}-\mathbf{\text{n}}(\mathbf{\text{n}}\cdot \mathrm{\nabla}),$$

3.2

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 **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

$$\begin{array}{rl}{\mathbf{\text{f}}}^{\mathrm{D},\mathrm{G}}& ={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{G}}+{\mathbf{\text{t}}}_{\mathrm{G}},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{\mathrm{G},}\\ {\mathbf{\text{f}}}^{\mathrm{D},\mathrm{L}}& ={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{L}}+{\mathbf{\text{t}}}_{\mathrm{L}},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{\mathrm{L}}\\ \mathrm{and}\phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}{\mathbf{\text{f}}}^{\mathrm{D},\mathrm{S}}& ={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{S}}+{\mathbf{\text{t}}}_{\mathrm{S}},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{\mathrm{S}},\end{array}$$

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

$$\begin{array}{rl}{\mathbf{\text{f}}}_{\mathrm{GL}}^{\mathrm{D}}& :={\mathbf{\text{f}}}^{\mathrm{D},\mathrm{G}}+{\mathbf{\text{f}}}^{\mathrm{G},\mathrm{L}}={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{GL}}+({\mathbf{\text{t}}}_{\mathrm{G}}+{\mathbf{\text{t}}}_{\mathrm{L}}),\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{\mathrm{GL}},\end{array}$$

3.3

$$\begin{array}{rl}{\mathbf{\text{f}}}_{\mathrm{LS}}^{\mathrm{D}}& :={\mathbf{\text{f}}}^{\mathrm{D},\mathrm{L}}+{\mathbf{\text{f}}}^{\mathrm{D},\mathrm{S}}={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{LS}}+({\mathbf{\text{t}}}_{\mathrm{L}}+{\mathbf{\text{t}}}_{\mathrm{S}}),\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{\mathrm{LS}}\end{array}$$

3.4

$$\begin{array}{rl}\mathrm{and}\phantom{\rule{2em}{0ex}}{\mathbf{\text{f}}}_{\mathrm{GS}}^{\mathrm{D}}& :={\mathbf{\text{f}}}^{\mathrm{D},\mathrm{G}}+{\mathbf{\text{f}}}^{\mathrm{D},\mathrm{S}}={\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{GS}}+({\mathbf{\text{t}}}_{\mathrm{G}}+{\mathbf{\text{t}}}_{\mathrm{S}}),\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{\mathrm{GS}},\end{array}$$

3.5

where

$${\mathit{\varsigma}}_{\mathrm{GL}}:={\mathit{\varsigma}}_{\mathrm{G}}+{\mathit{\varsigma}}_{\mathrm{L}};\phantom{\rule{1em}{0ex}}{\mathit{\varsigma}}_{\mathrm{LS}}:={\mathit{\varsigma}}_{\mathrm{L}}+{\mathit{\varsigma}}_{\mathrm{S}}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{\mathit{\varsigma}}_{\mathrm{GS}}:={\mathit{\varsigma}}_{\mathrm{G}}+{\mathit{\varsigma}}_{\mathrm{S}}.$$

3.6

In general, on a specific interface *Γ*_{k},*k*,=*G*,*L*,*S*, and *k*≠, we have **n**_{k}=−**n**_{}. Subsequently, we may denote

$$\begin{array}{rl}{\mathbf{\text{t}}}_{k}+{\mathbf{\text{t}}}_{\ell}& =({\mathit{\sigma}}_{k}-{\mathit{\sigma}}_{l})\cdot {\mathbf{\text{n}}}_{k}:=[{\mathbf{\text{t}}}_{k}],\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{k\ell}\phantom{\rule{1em}{0ex}}\mathrm{or}\end{array}$$

3.7

$$\begin{array}{rl}{\mathbf{\text{t}}}_{k}+{\mathbf{\text{t}}}_{\ell}& =({\mathit{\sigma}}_{\ell}-{\mathit{\sigma}}_{k})\cdot {\mathbf{\text{n}}}_{\ell}:=[{\mathbf{\text{t}}}_{\ell}],\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{k\ell}.\end{array}$$

3.8

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

$${\mathbf{\text{f}}}_{k\ell}^{\mathrm{D}}={\overline{\rho}}_{k\ell}{v}_{k}^{\mathrm{intf}}({\mathbf{\text{v}}}_{\ell}-{\mathbf{\text{v}}}_{k})\phantom{\rule{1em}{0ex}}\mathrm{or}\phantom{\rule{1em}{0ex}}{\mathbf{\text{f}}}_{k\ell}^{\mathrm{D}}={\overline{\rho}}_{k\ell}{v}_{\ell}^{\mathrm{intf}}({\mathbf{\text{v}}}_{k}-{\mathbf{\text{v}}}_{\ell}),$$

3.9

where ${v}_{\ell}^{\mathrm{intf}}$ is the magnitude of the interface normal velocity (the superscript ‘intf’ means interface), and ${\overline{\rho}}_{k\ell}$ is the average mass density on the interface (*k*).

Therefore, the dynamic equations of interface MCL read as

$${\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k\ell}+[{\mathbf{\text{t}}}_{k}]={\overline{\rho}}_{k\ell}{v}_{k}^{\mathrm{intf}}({\mathbf{\text{v}}}_{\ell}-{\mathbf{\text{v}}}_{k}).$$

3.10

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

$$\begin{array}{rl}{\mathit{\varsigma}}_{\mathrm{G}}& ={\gamma}_{\mathrm{G}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)},\end{array}$$

3.11

$$\begin{array}{rl}{\mathit{\varsigma}}_{\mathrm{L}}& ={\gamma}_{\mathrm{L}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+{\mathrm{\nabla}}_{\mathrm{s}}{\gamma}_{\mathrm{L}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+{\mu}_{\mathrm{S}}{\mathbf{\text{d}}}_{\mathrm{s}}\end{array}$$

3.12

$$\begin{array}{rl}\mathrm{and}\phantom{\rule{2em}{0ex}}{\mathit{\varsigma}}_{\mathrm{S}}& ={\gamma}_{\mathrm{S}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+\frac{\mathrm{\partial}{\mathit{\Gamma}}_{\mathrm{S}}}{\mathrm{\partial}{\mathit{\u03f5}}_{\mathrm{s}}}+{\gamma}_{\mathrm{S}}{\mathrm{\nabla}}_{\mathrm{s}}\otimes \mathbf{\text{u}},\end{array}$$

3.13

where **u** is the three-dimensional surface displacement field; *μ*_{L} is the surface viscosity. Note that ${\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}$ denotes the unit tensor on a smooth surface or two-dimensional manifold, which is defined as

$${\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}:=\mathbf{\text{P}}\mathbf{\text{I}}=\mathbf{\text{P}},$$

3.14

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

$$\mathbf{\text{P}}:=\mathbf{\text{I}}-\mathbf{\text{n}}\otimes \mathbf{\text{n}},$$

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. **d**_{s} is the surface or interface rate of deformation that is defined as

$${\mathbf{\text{d}}}_{\mathrm{s}}:=\mathbf{\text{P}}\mathrm{Sym}(\mathrm{\nabla}\otimes \mathbf{\text{v}}){\mathbf{\text{P}}}^{\mathrm{T}}=\mathbf{\text{P}}\mathrm{Sym}(\mathrm{\nabla}\otimes \mathbf{\text{v}})\mathbf{\text{P}},$$

3.15

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

$${\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}={\mathrm{\nabla}}_{\mathrm{s}}\cdot \mathbf{\text{P}}=2\kappa \mathbf{\text{n}},$$

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; _{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:

$${\mathit{\Gamma}}_{\mathrm{S}}=\frac{1}{2}{\u03f5}_{ij}^{\mathrm{s}}{C}_{ijk\ell}^{\mathrm{S}}{\u03f5}_{k\ell}^{\mathrm{s}},\phantom{\rule{1em}{0ex}}i,j,k,\ell =1,2,$$

3.16

in which the surface elastic tensor is related to surface tension *γ*_{S} as well,

$${C}_{ijk\ell}^{\mathrm{S}}=({\mathrm{\lambda}}_{\mathrm{S}}+{\gamma}_{\mathrm{S}}){\delta}_{ij}{\delta}_{k\ell}+{\mu}_{\mathrm{S}}({\delta}_{ik}{\delta}_{j\ell}+{\delta}_{i\ell}{\delta}_{jk}),\phantom{\rule{1em}{0ex}}i,j,k,\ell =1,2,$$

3.17

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

$${\mathit{\varsigma}}_{\mathrm{S}}={\gamma}_{\mathrm{S}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+2({\mu}_{\mathrm{S}}-{\gamma}_{\mathrm{S}}){\mathit{\u03f5}}_{\mathrm{s}}+({\mathrm{\lambda}}_{\mathrm{S}}+{\gamma}_{\mathrm{S}})\hspace{0.17em}\mathrm{tr}({\mathit{\u03f5}}_{\mathrm{s}}){\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+{\gamma}_{\mathrm{S}}{\mathrm{\nabla}}_{\mathrm{s}}\otimes \mathbf{\text{u}}.$$

3.18

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

$$\begin{array}{rl}{\mathit{\varsigma}}_{\mathrm{LS}}& ={\gamma}_{\mathrm{LS}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+{\mathrm{\nabla}}_{\mathrm{s}}{\gamma}_{\mathrm{L}}{\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+2({\mu}_{\mathrm{S}}-{\gamma}_{\mathrm{S}}){\mathit{\u03f5}}_{\mathrm{s}}+({\mathrm{\lambda}}_{\mathrm{S}}+{\gamma}_{\mathrm{S}})\hspace{0.17em}\mathrm{tr}({\mathit{\u03f5}}_{\mathrm{s}}){\mathbf{\text{I}}}_{\mathrm{s}}^{(2)}+{\gamma}_{\mathrm{S}}{\mathrm{\nabla}}_{\mathrm{s}}\otimes \mathbf{\text{u}}+{\mu}_{\mathrm{L}}{\mathbf{\text{d}}}_{\mathrm{s}},\end{array}$$

3.19

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

$${\mathit{\u03f5}}_{\mathrm{s}}:=\mathbf{\text{P}}\cdot \mathit{\u03f5}\cdot \mathbf{\text{P}},$$

3.20

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

$$\mathit{\u03f5}=\frac{1}{2}(\mathbf{\text{I}}-\mathbf{\text{c}})\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}\mathbf{\text{c}}:={\mathbf{\text{F}}}^{-\mathrm{T}}\cdot {\mathbf{\text{F}}}^{-1},$$

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

In general, the surface adhesive stress tensor

$${\mathit{\sigma}}_{k}^{\mathrm{adh}}=({\beta}_{k}{\beta}_{\ell}{\int}_{\mathrm{\partial}{\mathit{\Omega}}_{\ell}}{\mathbf{\text{n}}}_{\ell}\otimes {\mathbf{\text{s}}}_{kell}\psi (s)\hspace{0.17em}\mathrm{d}{a}_{\ell})\ne {\mathit{\sigma}}_{\ell}^{\mathrm{adh}}=({\beta}_{k}{\beta}_{\ell}{\int}_{\mathrm{\partial}{\mathit{\Omega}}_{k}}{\mathbf{\text{n}}}_{k}\otimes {\mathbf{\text{s}}}_{\ell k}\psi (s)\hspace{0.17em}\mathrm{d}{a}_{k}),$$

thus [**t**_{k}]≠0. Consider a simple case of gaseous and liquid interface. Let

$${\mathit{\sigma}}_{\mathrm{G}}=-{p}_{\mathrm{G}}{\mathbf{\text{n}}}_{\mathrm{G}}\otimes {\mathbf{\text{n}}}_{\mathrm{G}},\phantom{\rule{1em}{0ex}}{\mathit{\sigma}}_{\mathrm{L}}=-{p}_{\mathrm{L}}{\mathbf{\text{n}}}_{\mathrm{L}}\otimes {\mathbf{\text{n}}}_{\mathrm{L}}\to [{\mathbf{\text{t}}}_{\mathrm{L}}]=({p}_{\mathrm{G}}-{p}_{\mathrm{L}}){\mathbf{\text{n}}}_{\mathrm{L}}\leftarrow {\mathbf{\text{n}}}_{\mathrm{G}}=-{\mathbf{\text{n}}}_{\mathrm{L}}.$$

The interface MCL equation yields

$${\mathrm{\nabla}}_{\mathrm{s}}\cdot ({\gamma}_{\mathrm{GL}}{\mathbf{\text{I}}}_{\mathrm{S}}^{(2)})+({p}_{\mathrm{G}}-{p}_{\mathrm{L}}){\mathbf{\text{n}}}_{\mathrm{L}}=\overline{\rho}{v}_{\mathrm{L}}^{\mathrm{intf}}({\mathbf{\text{v}}}_{\mathrm{L}}-{\mathbf{\text{v}}}_{\mathrm{G}}).$$

3.21

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

$${p}_{\mathrm{L}}-{p}_{\mathrm{G}}=2{\gamma}_{\mathrm{GL}}{\kappa}_{\mathrm{GL}},$$

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

By using the Lagrangian description, we can obtain the relationship between the unit out-normal vector in the current and reference configurations as

$$\mathbf{n}=\frac{{\mathbf{F}}^{-\mathrm{T}}\mathbf{N}}{\sqrt{\mathbf{N}\cdot {\mathbf{C}}^{-1}\mathbf{N}}},$$

3.22

where **C** is the right Cauchy–Green tensor

$$\mathbf{C}={\mathbf{F}}^{\mathrm{T}}\mathbf{F}.$$

3.23

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

$$2\kappa =-\mathrm{\nabla}\cdot \mathbf{\text{n}}.$$

3.24

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

$$\begin{array}{rl}2\kappa & =\frac{{\mathbf{F}}^{\mathrm{T}}\mathcal{G}\vdots ({\mathbf{C}}^{-1}\mathbf{N}\otimes {\mathbf{C}}^{-1}\mathbf{N}\otimes {\mathbf{C}}^{-1}\mathbf{N})-{\mathrm{\nabla}}_{\mathbf{X}}\mathbf{N}:({\mathbf{C}}^{-1}\mathbf{N}\otimes {\mathbf{C}}^{-1}\mathbf{N})}{{(\mathbf{N}\cdot {\mathbf{C}}^{-1}\mathbf{N})}^{3/2}}\\ & \phantom{\rule{1em}{0ex}}-({\mathbf{n}}^{\mathrm{T}}\mathcal{G}):{\mathbf{C}}^{-1}+\frac{1}{\sqrt{\mathbf{N}\cdot {\mathbf{C}}^{-1}\mathbf{N}}}{\mathbf{C}}^{-1}:{\mathrm{\nabla}}_{\mathbf{X}}\mathbf{N},\end{array}$$

3.25

where $\mathcal{G}$ is the derivative of the deformation gradient tensor, and it is defined as

$$\mathcal{G}:=\frac{{\mathrm{\partial}}^{2}\mathbf{x}}{\mathrm{\partial}\mathbf{X}\otimes \mathrm{\partial}\mathbf{X}}.$$

3.26

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.

For a monolithic solution, we must solve the following coupled different equations simultaneously:

$$\begin{array}{rl}\mathrm{\nabla}\cdot {\mathit{\sigma}}_{k}+{\rho}_{k}{\mathbf{\text{b}}}_{k}& ={\rho}_{k}{\ddot{\mathbf{\text{u}}}}_{k}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{\mathbf{\text{t}}}_{k}={\mathit{\sigma}}_{k}{\mathbf{\text{n}}}_{k},\hspace{0.17em}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{tk};\end{array}$$

3.27

$$\begin{array}{rl}{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k\ell}+({\mathbf{\text{t}}}_{k}-{\mathbf{\text{t}}}_{\ell})& =\overline{\rho}{v}_{k}^{\mathrm{inft}}({\mathbf{\text{v}}}_{\ell}-{\mathbf{\text{v}}}_{k});\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{k\ell}\end{array}$$

3.28

$$\begin{array}{rl}\mathrm{and}\phantom{\rule{2em}{0ex}}\mathrm{\nabla}\cdot {\mathit{\sigma}}_{\ell}+{\rho}_{\ell}{\mathbf{\text{b}}}_{\ell}& ={\rho}_{\ell}{\ddot{\mathbf{\text{u}}}}_{\ell}\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{\mathbf{\text{t}}}_{\ell}={\mathit{\sigma}}_{\ell}{\mathbf{\text{n}}}_{\ell},\hspace{0.17em}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{t\ell},\end{array}$$

3.29

in which *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.

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

$${\mathbf{\text{t}}}_{k}=\overline{\rho}{v}_{k}^{\mathrm{intf}}({\mathbf{\text{v}}}_{\ell}-{\mathbf{\text{v}}}_{k})-{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k\ell}-{\mathbf{\text{t}}}_{\ell}.$$

3.30

In the above equation, both traction forces **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.

$${\mathbf{\text{t}}}_{\ell}\approx {\mathbf{\text{t}}}_{\ell}^{\mathrm{adh}}={\mathit{\sigma}}_{\ell}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\ell}=-{\mathit{\sigma}}_{\ell}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{k}.$$

We call this procedure *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

$${\mathbf{\text{t}}}_{k}=-{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k\ell}+{\mathit{\sigma}}_{\ell}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{k},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{tk}.$$

3.31

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

$${\mathbf{\text{t}}}_{k}\approx {\mathbf{\text{t}}}_{k}^{\mathrm{adh}}={\mathit{\sigma}}_{k}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{k}=-{\mathit{\sigma}}_{k}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\ell},$$

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

$${\mathbf{\text{t}}}_{\ell}=-{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k\ell}+{\mathit{\sigma}}_{k}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\ell},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in \mathrm{\partial}{\mathit{\Omega}}_{t\ell}.$$

3.32

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

$${\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{k\ell}+[{\mathbf{\text{t}}}_{k}]={\overline{\rho}}_{k\ell}{v}_{k}^{\mathrm{inf}}({\mathbf{\text{v}}}_{k}-{\mathbf{\text{v}}}_{\ell}),\phantom{\rule{1em}{0ex}}\mathrm{where}\hspace{0.17em}{\mathit{\varsigma}}_{kl}={\gamma}_{k\ell}{\mathbf{\text{I}}}_{\mathrm{s}},\hspace{0.17em}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{k\ell},\hspace{0.17em}k\ell =\text{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 *Ω*_{L}=*Γ*_{LS}(L)*Γ*_{GL}(L),

$${\mathbf{\text{t}}}_{\mathrm{L}}=-{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\varsigma}_{\mathrm{GL}}+{\mathit{\sigma}}_{\mathrm{G}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{L})$$

3.33

and

$${\mathbf{\text{t}}}_{\mathrm{L}}=-{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\varsigma}_{\mathrm{LS}}+{\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L}).$$

3.34

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

$${\mathbf{\text{t}}}_{\mathrm{G}}=-{\mathrm{\nabla}}_{\mathrm{s}}\cdot {\varsigma}_{\mathrm{GS}}+{\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{G}},\phantom{\rule{1em}{0ex}}\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{G}),$$

3.35

where ${\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{G}}$ 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:

$${\mathbf{u}}_{\mathrm{s}}(\mathbf{\text{x}},t)=\sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{Snode}}}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathbf{d}}_{{I}_{\mathrm{s}}}(t)\phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{\mathbf{w}}_{\mathrm{s}}(\mathbf{\text{x}},t)=\sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{Snode}}}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathbf{w}}_{{I}_{\mathrm{s}}}(t),$$

3.36

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

$$\begin{array}{rl}& \sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{GL}}}{\mathbf{\text{w}}}_{{I}_{\mathrm{s}}}\{-{\int}_{{\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{GL}}\hspace{0.17em}\mathrm{d}s+{\int}_{{\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\sigma}}_{\mathrm{G}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}}\hspace{0.17em}\mathrm{d}s\},\end{array}$$

3.37

$$\begin{array}{rl}& \sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{LS}}}{\mathbf{\text{w}}}_{{I}_{\mathrm{s}}}\cdot \{-{\int}_{{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{LS}}\hspace{0.17em}\mathrm{d}s+{\int}_{{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}}\hspace{0.17em}\mathrm{d}s\}\end{array}$$

3.38

$$\begin{array}{rl}\mathrm{and}\phantom{\rule{2em}{0ex}}& \sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{GS}}}{\mathbf{\text{w}}}_{{I}_{\mathrm{s}}}\cdot \{-{\int}_{{\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{G})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathrm{\nabla}}_{\mathrm{s}}\cdot {\mathit{\varsigma}}_{\mathrm{GS}}\hspace{0.17em}\mathrm{d}s+{\int}_{{\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{G})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{G}}\hspace{0.17em}\mathrm{d}s\}.\end{array}$$

3.39

Integration by parts yields

$$\begin{array}{rl}& \sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{GL}}}\{-{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\varsigma}}_{\mathrm{GL}}\cdot {\mathbf{\text{q}}}_{t}{|}_{\mathrm{\partial}{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}+{\int}_{{\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{L})}\frac{\mathrm{\partial}{N}_{{I}_{\mathrm{s}}}}{\mathrm{\partial}\mathbf{\text{x}}}:{\mathit{\varsigma}}_{\mathrm{GL}}\hspace{0.17em}\mathrm{d}s+{\int}_{{\mathit{\Gamma}}_{\mathrm{GL}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\sigma}}_{\mathrm{G}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}}\hspace{0.17em}\mathrm{d}s\},\\ & \sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{LS}}}\{-{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\varsigma}}_{\mathrm{LS}}\cdot {\mathbf{\text{q}}}_{1}{|}_{\mathrm{\partial}{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}+{\int}_{{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}\frac{\mathrm{\partial}{N}_{{I}_{\mathrm{s}}}}{\mathrm{\partial}\mathbf{\text{x}}}:{\mathit{\varsigma}}_{\mathrm{LS}}\hspace{0.17em}\mathrm{d}s+{\int}_{{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}}\hspace{0.17em}\mathrm{d}s\}\\ \mathrm{and}\phantom{\rule{2em}{0ex}}& \sum _{{I}_{\mathrm{s}}=1}^{{N}_{\mathrm{GS}}}\{-{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\varsigma}}_{\mathrm{LS}}\cdot (-{\mathbf{\text{q}}}_{1}){|}_{\mathrm{\partial}{\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{G})}+{\int}_{{\mathit{\Gamma}}_{\mathrm{GS}}(\mathrm{G})}\frac{\mathrm{\partial}{N}_{{I}_{\mathrm{s}}}}{\mathrm{\partial}\mathbf{\text{x}}}:{\mathit{\varsigma}}_{\mathrm{LS}}\hspace{0.17em}\mathrm{d}s+{\int}_{{\mathit{\Gamma}}_{\mathrm{LS}}(\mathrm{L})}{N}_{{I}_{\mathrm{s}}}(\mathbf{\text{x}}){\mathit{\sigma}}_{\mathrm{S}}^{\mathrm{adh}}\cdot {\mathbf{\text{n}}}_{\mathrm{L}}\hspace{0.17em}\mathrm{d}s\},\end{array}$$

where vectors **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. $\mathrm{\forall}\hspace{0.17em}\mathbf{\text{x}}\in {\mathcal{L}}_{\mathrm{mcl}}={\mathit{\Gamma}}_{\mathrm{GL}}\cap {\mathit{\Gamma}}_{\mathrm{LS}}\cap {\mathit{\Gamma}}_{\mathrm{GS}}$, we have *N*_{Is}(**x**)=1, and the external resultant acting on the MCL is

$${\mathbf{F}}_{\mathrm{L}}^{\mathrm{mcl}}=-{\mathit{\varsigma}}_{\mathrm{GL}}\cdot {\mathbf{\text{t}}}_{\mathrm{GL}}-{\mathit{\varsigma}}_{\mathrm{GS}}\cdot {\mathbf{\text{t}}}_{\mathrm{GS}}-{\mathit{\varsigma}}_{\mathrm{LS}}\cdot {\mathbf{\text{t}}}_{\mathrm{LS}}$$

where

$${\mathbf{\text{t}}}_{\mathrm{GL}}={\mathbf{\text{q}}}_{t}={\mathbf{\text{q}}}_{1}\mathrm{cos}\theta -{\mathbf{\text{q}}}_{2}\mathrm{sin}\theta \phantom{\rule{1em}{0ex}}\mathrm{and}\phantom{\rule{1em}{0ex}}{\mathbf{\text{t}}}_{\mathrm{GS}}=-{\mathbf{\text{q}}}_{1},\phantom{\rule{1em}{0ex}}{\mathbf{\text{t}}}_{\mathrm{LS}}={\mathbf{\text{q}}}_{1},$$

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

$$\begin{array}{rl}{\mathit{\varsigma}}_{\mathrm{GL}}& ={\gamma}_{\mathrm{GL}}{\mathbf{\text{I}}}_{\mathrm{GL}}^{\mathrm{S}}={\gamma}_{\mathrm{GL}}({\mathbf{\text{q}}}_{3}\otimes {\mathbf{\text{q}}}_{3}+{\mathbf{\text{q}}}_{t}\otimes {\mathbf{\text{q}}}_{t}),\\ {\mathit{\varsigma}}_{\mathrm{GS}}& ={\gamma}_{\mathrm{GS}}{\mathbf{\text{I}}}_{\mathrm{GS}}^{S}={\gamma}_{\mathrm{GS}}({\mathbf{\text{q}}}_{3}\otimes {\mathbf{\text{q}}}_{3}+{\mathbf{\text{q}}}_{1}\otimes {\mathbf{\text{q}}}_{1})\\ \mathrm{and}\phantom{\rule{2em}{0ex}}\phantom{\rule{1em}{0ex}}{\mathit{\varsigma}}_{\mathrm{LS}}& ={\gamma}_{\mathrm{LS}}{\mathbf{\text{I}}}_{\mathrm{LS}}^{S}={\gamma}_{\mathrm{LS}}({\mathbf{\text{q}}}_{3}\otimes {\mathbf{\text{q}}}_{3}+{\mathbf{\text{q}}}_{1}\otimes {\mathbf{\text{q}}}_{1}),\end{array}$$

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

$$\begin{array}{rl}{\mathbf{F}}_{\mathrm{L}}^{\mathrm{mcl}}& ={\gamma}_{\mathrm{GS}}{\mathbf{\text{q}}}_{1}-{\gamma}_{\mathrm{LS}}{\mathbf{\text{q}}}_{1}-{\gamma}_{\mathrm{GL}}({\mathbf{\text{q}}}_{1}\mathrm{cos}\theta -{\mathbf{\text{q}}}_{2}\mathrm{sin}\theta )\end{array}$$

3.40

$$\begin{array}{rl}& =\left[\begin{array}{c}{\gamma}_{\mathrm{GS}}-{\gamma}_{\mathrm{LS}}-{\gamma}_{\mathrm{GL}}\mathrm{cos}\theta \\ {\gamma}_{\mathrm{GL}}\mathrm{sin}\theta \end{array}\right]\left[\begin{array}{c}{\mathbf{\text{q}}}_{1}\\ {\mathbf{\text{q}}}_{2}\end{array}\right].\end{array}$$

3.41

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,

$${\gamma}_{\mathrm{GS}}={\gamma}_{\mathrm{LS}}+{\gamma}_{\mathrm{GL}}\mathrm{cos}\theta ,$$

3.42

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

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

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}Pas, and the surface tension between the droplet and the atmosphere *γ*_{GL}=7.28×10^{−2}Nm^{−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*=5nm, 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.

Time history of the simulation of a three-dimensional ellipsoidal droplet embedded in atmosphere, driven by the surface tension effect. (*a*) *t*=0.004ns, (*b*) *t*=0.04ns, (*c*) *t*=0.08ns, (*d*) *t*= 0.2ns, (*e*) *t*=0.28ns **...**

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:

$$\mathrm{\Delta}P=\frac{2\gamma}{R},$$

4.1

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.62MPa, the final radius of the spherical droplet *R*=5.71nm. The surface tension *γ*=*γ*_{LG}=72.75×mNm^{−1}. Note that

$$\frac{2\gamma}{R}=\frac{2\times 72.75\times {10}^{-3}}{5.71\times {10}^{-9}}\hspace{0.17em}\mathrm{Pa}=25.48\times {10}^{6}\hspace{0.17em}\mathrm{Pa}=25.48\hspace{0.17em}\mathrm{MPa},$$

4.2

and we end up with an error percentage of 0.54% for the pressure change compared to the value Δ*P*=25.62MPa, 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.

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*=5nm, and the substrate is a short cylinder shape platform with the dimension *H*×*R*=4.5×20nm, 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 5nm 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}kgm^{−3}, viscosity *μ*=0.6×10^{−3}Pas and bulk modulus *κ*=2.2×10^{9}Pa; for the substrate (copper): we have the density *ρ*^{s}=8.94×10^{3}kgm^{−3}, Young's modulus *E*=1.2×10^{11}Pa, and Poisson's ratio *ν*=0.34.

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×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: ${\sigma}^{\mathrm{ds}}=0.5({\sigma}^{d}+{\sigma}^{s}),{\epsilon}^{\mathrm{ds}}=\sqrt{{\epsilon}^{d}{\epsilon}^{s}}$

In this simulation, the time step is Δ*t*=5×10^{−15}s and the total steps is *n*steps=200000. 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*).

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.010ns, (*b*) *t*=0.020ns, (*c*) *t*=0.030ns, (*d*) *t*=0.040ns, (*e*) *t*=0.070ns, **...**

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 ${\gamma}_{\mathrm{GL}}=7.28\times {10}^{-2}\hspace{0.17em}\mathrm{N}\hspace{0.17em}{\mathrm{m}}^{-1}$. The surface material parameters of the substrate are *γ*_{S}=1.0398Nm^{−1}, surface Lame parameters λ^{S}=15.6Nm^{−1}, *μ*^{S}=−8.6Nm^{−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*=15ps, which indicates that the early-time dynamics are independent of the wettability.

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 (15657 nodes, 13892 elements) and mesh 3 (21253 nodes, 19300 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*=52ps. 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.

Droplet configurations for the three different meshes at *t*=52ps. (*a*) mesh 1, (*b*) mesh 2 and (*c*) mesh 3. (Online version in colour.)

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*=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.

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. *t*_{0}=1μs; *R*_{0}=1μm. **...**

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*=5nm. The thickness of the water film is *H*=4.5nm. 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=1000000. From *t*=0 to 0.025ns, 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.025ns to *t*=0.05s, 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.05ns, 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.

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.

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.

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.

The authors declare no competing interests in this work.

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, 10486–10488. (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**

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |