|Home | About | Journals | Submit | Contact Us | Français|
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 , 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  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.
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,
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; 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  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 ϕ(|xk|), where xk=x−xk,k,=G,L,S, k≠ and xkΩk,xΩ. Then the total interphase adhesive-contact potential energy Πac for two interacting phase is
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
where we define the adhesive body force in each phase as
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 xkΩk and xΩ so that sk:=x−xk=:−sk. Following the approach outlined by Jagota & Argento , 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Ωk due to the presence of an infinitesimal surface element surface daΩ is expressed as
where nk and n are the unit surface out-normal at points xk and x, respectively, k,=G,L,S, k≠, and
According to Jagota & Argento , 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Ω due to the presence of an infinitesimal surface element surface dakΩk is expressed as
It may be noted that dFk≠−dF, 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. . Nk and N are the total number of infinitesimal surface elements on Ωk and Ω, respectively.
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Ω at xkΩk. Thus, if we perform the integration over the surface Ω, we may obtain the adhesive surface stress
which represents the surface stress tensor at point xkΩk, due to the interaction from the phase . The corresponding surface traction at point xkΩk can then be expressed as
Similarly, we can derive the surface stress tensor and traction at point xΓk() as
Note that if one of the two interacting phases is rigid, we may be able to analytically integrate equation (2.15) to obtain the adhesive surface stress tensor for a given interatomic potential ϕ(r). For instance, for the 12-6 Lennard-Jones potential given by equation (2.6), if the phase l is a rigid infinite plate as shown in figure 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 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.
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.
In this section we discuss the dynamic interface moving contact line theory.
where fD,k is the surface D'Alembert force density, ςk is the surface stress, tk is the adhesive surface traction vector and 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
which allows us to introduce the following interface D'Alembert force densities:
In general, on a specific interface Γk,k,=G,L,S, and k≠, we have nk=−n. Subsequently, we may denote
where is the magnitude of the interface normal velocity (the superscript ‘intf’ means interface), and is the average mass density on the interface (k).
Therefore, the dynamic equations of interface MCL read as
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 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
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; 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  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.
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. [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
where C is the right Cauchy–Green tensor
The mean curvature is related to the divergence of the surface unit out-normal (e.g. ) as
where is the derivative of the deformation gradient tensor, and it is defined as
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:
in which k and are two adjacent phases, and Γk is the interface.
In (3.27), the interface traction jump [t]=(tk−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. ). 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 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
In the above equation, both traction forces tk and t are unknown. In order to solve tk 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.
We call this procedure the unilateral interface traction approximation.
If we adopt the no-slip boundary condition, the relative slip velocity vanishes, i.e. vk−v=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 , we can obtain the traction on the boundary of the -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
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
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),
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 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
Integration by parts yields
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. , 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 .
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×109Pa, μ=0.6×10−3Pas, and the surface tension between the droplet and the atmosphere γGL=7.28×10−2Nm−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 dt=2.0×10−14s. Time history of the ellipse changing to perfect sphere process are shown in figure 6. We plot the evolution of the lengths of longest and shortest semi-principal axes of the ellipsoidal droplet and the pressure change at the centre of the droplet in figure 7.
Since the relaxation process is a dynamic process, so it is expected that there are some oscillations in the length of the semi-principal axes at the beginning of the simulation. But eventually, both the principal axes and the pressure at the centre of the droplet will reach to the stable values, owing to the damping effect of the droplet viscosity. To check the validity of the final state, we would like to confirm whether the following Young–Laplace equation for a spherical droplet is satisfied or not:
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
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 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 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 8b. The material properties are given as follows. For the droplet (water): we choose the density ρd=1.0×103kgm−3, viscosity μ=0.6×10−3Pas and bulk modulus κ=2.2×109Pa; for the substrate (copper): we have the density ρs=8.94×103kgm−3, Young's modulus E=1.2×1011Pa, and Poisson's ratio ν=0.34.
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−10m, εds=8.4725×10−21J, which are obtained based on the Lorentz-Berthelot rule:
In this simulation, the time step is Δt=5×10−15s and the total steps is nsteps=200000. Reduced units are used throughout the simulation, with the unit of time t0=1.0×10−12s, length l0=1.0×10−10m and mass m0=1.0×10−27kg.
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 . 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. . The only unknown parameter for the interface layer is γL, which can be obtained based on the target contact angle . 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. . One can see that the dynamic contact angle of the MMCL are in general agreement with the result obtained by molecular dynamics simulations . 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.
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.
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 . 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 . 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.
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−16s, and the total simulation time steps is nsteps=1000000. From t=0 to 0.025ns, the rigid sphere is gradually pushed down at the speed of Δz=2.0×10−14m per time step (figure 17a–c). In time duration of t=0.025ns to t=0.05s, the sphere tip is pulled up with the same speed (figure 17d–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 17g–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 . 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.