Nanoscale sheared droplet: volume-of-fluid, phase-field and no-slip molecular dynamics

Abstract

The motion of the three-phase contact line between two immiscible fluids and a solid surface arises in a variety of wetting phenomena and technological applications. One challenge in continuum theory is the effective representation of molecular motion close to the contact line. Here, we characterize the molecular processes of the moving contact line to assess the accuracy of two different continuum two-phase models. Specifically, molecular dynamics simulations of a two-dimensional droplet between two moving plates are used to create reference data for different capillary numbers and contact angles. We use a simple-point-charge/extended water model. This model provides a very small slip and a more realistic representation of the molecular physics than Lennard-Jones models. The Cahn–Hilliard phase-field model and the volume-of-fluid model are calibrated against the drop displacement from molecular dynamics reference data. It is shown that the calibrated continuum models can accurately capture droplet displacement and droplet break-up for different capillary numbers and contact angles. However, we also observe differences between continuum and atomistic simulations in describing the transient and unsteady droplet behaviour, in particular, close to dynamical wetting transitions. The molecular dynamics of the sheared droplet provide insight into the line friction experienced by the advancing and receding contact lines. The presented results will serve as a stepping stone towards developing accurate continuum models for nanoscale hydrodynamics.

1. Introduction

The motion of a two-phase interface over a solid surface has turned out to be a challenging problem for continuum fluid mechanics (CFM) models. This is most evident for models that assume a sharp interface between the phases and the no-slip velocity condition at the solid surface. Under these classical assumptions, one ends up naturally with an immobile line at which both fluid phases meet with the solid, the so-called contact line. Clearly, a fixed contact line is in contradiction with, for example, our observations of a spreading drop on a surface, or of a liquid imbibition into a porous medium. This theoretical issue was identified already half a century ago by Huh & Scriven (1971) as a stress singularity at the contact line. Motivated by the importance of the moving contact line in applications such as printing (Kumar 2015), ${\rm CO}_2$ storage and water management in fuel cells (Singh et al. 2019), a number of approaches have been suggested for overcoming the stress singularity in Navier–Stokes based solvers. The extensive work on the topic (Bonn et al. 2009; Snoeijer & Andreotti 2013; Sui, Ding & Spelt 2014) suggests that there exists no single continuum approach for describing the moving contact line. Instead, different models are suitable depending on the application and the representation of interfacial physics, and each model comes with its own sets of empirical parameters. Certain guidelines are required to determine these parameters.

Despite the vast amount of theoretical developments that have been presented (Voinov 1976; Cox 1986; Shikhmurzaev 1994, 1997; Kalliadasis, Bielarz & Homsy 2000; Pismen & Pomeau 2000; Eggers 2004; Flitton & King 2004; Snoeijer 2006; Wilson et al. 2006; Pismen & Eggers 2008; Snoeijer & Eggers 2010; Chan, Snoeijer & Eggers 2012; Nold et al. 2018; Chan et al. 2020), a theoretical consensus is yet to be reached. To a large extent, progress is hindered by a lack of understanding of the nanoscale physics of wetting phenomena (Afkhami, Gambaryan-Roisman & Pismen 2020; Afkhami 2022). For example, a fundamental question relates to the nanoscopic contact angle. The typical choice until now has been a constant equilibrium angle (Kronbichler & Kreiss 2017). However, recent experimental (Deng et al. 2016) and numerical (Fernández-Toledano, Blake & De Coninck 2021) evidence suggests a contact line velocity-dependent dynamic contact angle even at nanoscale. The development of more accurate measurement techniques is ongoing (Thormann et al. 2008; Eriksson et al. 2019); however, as of yet, we lack a complete insight into interface shape and velocity field near moving contact lines in nanoscale. Atomistic simulations have provided significant insight into the molecular physics at this scale (Gentner, Ogonowski & De Coninck 2003; Smith et al. 2016, 2018; Perumanath et al. 2019; Lācis et al. 2020), although most systems used so far have been based on idealized force models between the liquid and the substrate.

In practice, there are few common approaches to model numerically moving contact lines in standard continuum methods. For sharp interface models, such as volume-of-fluid (VOF) and level-set (LS), the movement of the contact line is typically allowed by an explicit Navier-slip condition (Navier 1823; Spelt 2005) or by an implicit numerical slip at the contact line (Renardy, Renardy & Li 2001; Afkhami, Zaleski & Bussmann 2009). For diffuse interface models, such as the Cahn–Hilliard phase-field (PF) model (Jacqmin 2000), the contact line advances through diffusion even in the no-slip scenario. The lattice Boltzmann method (LBM) can be leveraged to solve different types of flow problems. The standard application of the LBM typically uses Shan & Chen (1993) or Gunstensen et al. (1991) models. It is argued that in the Shan–Chen model, the contact line moves through phase change (Kamali et al. 2011), similar to the LBM with thermodynamically consistent potentials (Briant, Wagner & Yeomans 2004). The LBM has also been used to solve PF equations (Briant & Yeomans 2004); consequently, the motion of the contact line then occurs through diffusion. In the current work, we focus particularly on two of these models, namely, geometric VOF and Cahn–Hilliard PF models.

For the VOF method, different components to model the contact line (static, dynamic angle, hysteresis window, Navier slip, etc.) have been proposed over time. A recent comparison between different options can be found in Legendre & Maglio (2015). They concluded that the models incorporating dynamic contact angle better represent the experiments of a spreading drop. For the PF model, there are guidelines to select the model parameters through calibration with experiments (Yue & Feng 2011). A part of the calibration is choosing the diffuse interface thickness in a manner to satisfy the so-called sharp interface limit (Yue, Zhou & Feng 2010; Magaletti et al. 2013; Xu, Di & Yu 2018). With this approach, a good agreement between PF simulations and capillary rise experiments has been demonstrated.

In general, however, the connection between the selected CFM models and molecular reality is not clear. It is not known how to choose the model parameters for an accurate representation of nanoscopic physics that determines the moving contact line speed. To address this, comparisons between molecular dynamics (MD) simulations and the chosen CFM models have been carried out (Qian, Wang & Sheng 2003; Barclay & Lukes 2019; Mohand et al. 2019). The typical MD work considers Lennard-Jones types of surfaces with large slip lengths. Only more recently have water MD simulations been carried out by some of us (Johansson & Hess 2018; Lācis et al. 2020) over surfaces with negligible slip.

In this work, we generate benchmark data from MD simulations of wetting over no-slip surfaces. These conditions can be reproduced employing simple-point-charge/extended (SPC/E) water on a smooth silica-like substrate (Johansson, Carlson & Hess 2015). We choose a forced wetting set-up (Blake et al. 2015) – instead of a capillary-driven one (Villanueva & Amberg 2006) – due to the more versatile control of the wetting process. In particular, we choose a sheared droplet configuration (see figure 1a), which is a well-studied canonical problem (Jacqmin 2004; Sbragaglia, Sugiyama & Biferale 2008; Gao & Lu 2013; Wang & McCarthy 2013) allowing simultaneous access to both receding and advancing contact lines. Furthermore, depending on the wall velocity $U_w$, the system either (i) stabilizes at a steady state (if $U_w < U_{w,c}$) or (ii) exhibits a non-trivial unsteady behaviour (if $U_w > U_{w,c}$). Here, $U_{w,c}$ is the critical wall velocity describing the boundary between (i) and (ii). Hence the sheared droplet configuration provides rich interfacial dynamics that are challenging to capture with CFM models.

Figure 1. (a) Dimensions and properties of the sheared-droplet configuration considered in the current work. (b) Close-up of the molecular system near the moving contact lines. Overlying solid mesh illustrates the binning boundaries for the collection of flow data from MD simulations.

By assuming a negligible slippage, the space of input parameters is essentially reduced to wall velocities ($U_w$) and equilibrium contact angles between water and silica ($\theta _0$). We then adopt a two-step approach. In the first step, we calibrate the continuum simulations against MD for a given pair $(U_w,\theta _0)$. This yields the necessary PF and VOF parameters that best reproduce the steady droplet displacement measured from MD. In the second step, we fix the calibrated PF and VOF parameters and assess the predictive capability of the CFM models for different $U_w$ by characterizing the interface shape and the drop displacement both in the stable and unstable regimes. This is an extension of the work presented by Lācis et al. (2020), who evaluated CFM model performance in matching the MD results for a single steady wall velocity ($U_w < U_{w,c}$) and a single equilibrium contact angle.

The MD of the sheared droplet reported herein provide insight into the friction experienced by the advancing and receding contact lines. We demonstrate asymmetric features of advancing and receding lines, and report evidence of large-scale temporal ‘stick-slip’ like oscillations. These observations not only enhance our physical understanding of moving contact lines, but also aid the development needed to increase the accuracy of continuum models.

The paper is organized as follows. In § 2, we describe the flow configuration and CFM models, and demonstrate the effect of the unknown parameters. The reference MD simulations of the sheared droplet system are described in § 3. We calibrate the CFM models against MD in § 4. Predictions from CFM are evaluated against MD in § 5. Then, in § 6, we provide insights into the molecular physics of the sheared droplet system. Following that, in § 7, limitations of CFM models, fluid slippage, friction and potential future modelling directions are discussed. We conclude the paper in § 8. In Appendices A–F, important physical and technical details are provided.

2. Flow configuration and continuum models

We consider a two-dimensional system that is periodic in the streamwise direction and bounded in the vertical direction by two parallel horizontal plates located at $y=0$ and $y=H$. A liquid drop of density $\rho _\ell$ and viscosity $\mu _\ell$ is sandwiched between the plates such that its maximum width is $W$. The drop is surrounded by water vapour with density $\rho _v$ and viscosity $\mu _v$. The surface tension between the phases is constant and denoted by $\sigma$. The numerical values of the geometry and fluid properties are reported in figure 1(a).

We study the response of the droplet to two parameters: (i) $U_w$, the constant velocity of the upper and lower walls moving in opposite directions, and (ii) the equilibrium contact angle $\theta _0$. The former is discussed interchangeably in its non-dimensional form, using the capillary number

(2.1)\begin{equation} {Ca} = \frac{2\mu_\ell U_w}{\sigma}. \end{equation}

The capillary number corresponding to critical wall velocity $U_{w,c}$ is ${Ca}_c$.

For small ${Ca}$, the liquid slips on the solid, resulting in a steady deformed droplet shape as shown in figure 2(a). Above a certain critical capillary number ${Ca}_c$, the liquid is deformed to such a degree that its interface to the surrounding vapour breaks, leaving behind multiple disconnected liquid droplets. The second control parameter, $\theta _0$, is a measure of the surface's affinity to water. In hydrophilic conditions ($\theta _0<90^\circ$), the affinity is strong, resulting in a larger droplet deformation compared to hydrophobic conditions ($\theta _0>90^\circ$). Under dynamic conditions, the contact line can be different from the equilibrium one. For advancing contact lines (marked ‘A’ in figure 2a), the liquid displaces the vapour, while for receding contact lines (marked ‘R’ in figure 2a), the vapour displaces the liquid. Both these processes are largely determined by molecular interactions between the water and the substrate as depicted in figure 1(b).

Figure 2. (a) Drop displacement measurement at the left ($\Delta x_l$) and right ($\Delta x_r$) sides at selected time instances. The interface angle is defined with respect to the horizontal line and is sampled as a function of the vertical coordinate $\theta (y)$. Note that $\theta (y)$ is measured only for ${Ca} < {Ca}_{c}$, while $\Delta x$ is measured for all $U_w$. (b) Variations of drop displacement $\Delta x = ( \Delta x_l + \Delta x_r )/2$ over time in VOF simulations for $\lambda = 0.47, 0.66, 0.94, 1.87, 3.74$ nm (increasing with the arrow). The equilibrium contact angle is $\theta _0 = 95^\circ$, and the capillary number is ${Ca} = 0.20 < {Ca}_{c}$.

The drop deformation induced by the moving walls is characterized with measures defined in figure 2(a). As a global measure, we introduce drop displacement of the left ($\Delta x_l$) and right ($\Delta x_r$) two-phase interface, respectively. For more detailed characterization, we also evaluate the interface angle $\theta (y)$ for steady configurations. It is obtained from the slope of each linear segment on the interface.

In the continuum setting, the flow and pressure fields $(\boldsymbol {{u}}, p)$ are obtained by solving the incompressible Navier–Stokes equations in the domain containing both phases. In two dimensions, the equations are

(2.2)$$\begin{gather} \rho ( x,y ) \left[ \frac{\partial \boldsymbol{{u}}}{\partial t} + \left(\boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{{u}} \right] ={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla} \boldsymbol{\cdot} \left[\mu ( x,y ) \left\{ \boldsymbol{\nabla} \boldsymbol{{u}} + (\boldsymbol{\nabla} \boldsymbol{{u}} )^{\rm T} \right\}\right] + \boldsymbol{{f}}_\sigma, \end{gather}$$

(2.3)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{{u}} = 0. \end{gather}$$

Here, $\boldsymbol {{f}}_\sigma$ is the surface tension force, and $\rho ( x,y )$ and $\mu ( x,y )$ are spatially dependent density and viscosity, respectively. The functions $\rho ( x,y )$ and $\mu ( x,y )$ take liquid and vapour values in the region occupied by each phase, and undergo a sharp transition at the boundary between the phases. In this transition region, the volume force $\boldsymbol {{f}}_\sigma$ is applied to model the force induced by the surface tension.

Zero wall-normal velocity, $u_y = 0$, is imposed on the moving impermeable walls. For the tangential velocity component, we impose a Navier-slip boundary condition

(2.4)\begin{equation} u_x = u_w + \ell_s\,\frac{\partial u_x}{\partial n}, \end{equation}

where $\ell _s$ is the slip length, $u_w$ is the wall velocity ($U_w$ at the top wall and $-U_w$ at the bottom wall), and $n$ is the wall-normal coordinate. Although it is expected that $\ell _s = 0$ for the chosen liquid–solid combination, the implementation of the CFM models allows for tests with non-zero value.

Equations (2.2) and (2.3) with corresponding boundary conditions are shared between the Cahn–Hilliard PF model and the geometric VOF model. However, the way surface tension and evolution of the two-phase interface are treated differs.

2.1. Geometric volume-of-fluid model

In this model, a phase variable $C(x,y)$ is defined as 0 in the gas and 1 in the liquid. The phase variable thus represents the liquid volume fraction in each cell. It satisfies the convection equation

(2.5)\begin{equation} \frac{\partial C}{\partial t} ={-} \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} C. \end{equation}

The surface tension force is applied by the continuous surface force (CSF) method (Brackbill, Kothe & Zemach 1992):

(2.6)\begin{equation} \boldsymbol{{f}}_\sigma ={-}\sigma \kappa \boldsymbol{\nabla} C, \end{equation}

where $\kappa$ is the curvature of the interface.

As a boundary condition for the phase variable, we impose a dynamic contact angle $\theta _{num}$. The contact angle is the angle between the interface ($\hat {n}_i$) and wall ($\hat {n}$) normals, so

(2.7)\begin{equation} \hat{n}_i\boldsymbol{\cdot} \hat{n} = \cos(\theta_{num}), \quad \text{where} \quad \hat{n}_i ={-}\frac{\boldsymbol{\nabla} C}{|\boldsymbol{\nabla} C|}, \end{equation}

and the sign might vary depending on convention. The relationship for $\theta _{num}$ is inspired by Cox's theory, as described and evaluated by Legendre & Maglio (2015). It takes the form

(2.8)\begin{equation} G(\theta_{num}) = G(\theta_0) + {Ca}_{cl}\ln\left(\frac{\varDelta/2}{\lambda} \right), \end{equation}

where $G$ is a monotonically increasing function, $\lambda$ is a microscopic cut-off length scale, and $\varDelta$ is the wall-normal cell height. We confirm the grid convergence reported by Legendre & Maglio (2015) in Appendix A.1. The capillary number ${Ca}_{cl}$ is based on the velocity of the contact line with respect to the wall. It is estimated by interpolating the velocity field linearly at the first grid cell. For the most hydrophobic and hydrophilic configurations, only the constant wall velocity was used for increased robustness.

The geometric VOF model is solved with an open-source finite volume code, called the PArallel, Robust, Interface Simulator (PARIS) (Aniszewski et al. 2021). The full set of equations (2.2), (2.3), (2.5) and (2.6) together with the boundary conditions ((2.4), (2.7) and periodic inlet/outlet) are discretized on a regular cuboid grid with a staggered spatial representation. The cell spacing is constant in all directions. More details about the numerical method are given in Appendix A.1. In order to compute the curvature near the wall for very hydrophobic droplets ($\theta _0 = 127^\circ$), we use a high-order scheme to approximate derivatives (Appendix A.2 provides further details).

The only unknown free parameter for VOF simulations is the length scale $\lambda$. To test the effect of $\lambda$, we select $\theta _0 = 95^\circ$ and ${Ca} = 0.20 < {Ca}_{c}$. The evolution of the drop displacement $\Delta x$ in time for $\lambda$ values $0.47$ nm, $0.66$ nm, $0.94$ nm, $1.87$ nm, and $3.74$ nm is shown in figure 2(b). Due to symmetry, $\Delta x_l = \Delta x_r = \Delta x$ in the CFM simulations, which we have verified numerically. We observe that larger $\lambda$ correspond to smaller steady $\Delta x$. As $\lambda$ is increased, the contact angle at the surface (2.8) and thus the interface shape (2.7) near the wall are modified. The interface curvature is modified in such a way that the surface tension force across the interface opposes the friction force from the wall, which leads to a smaller $\Delta x$. The $\Delta x$ for the initial time agrees between all $\lambda$ values. In this short period, the contact line does not slip with respect to the substrate, and the slope of $\Delta x (t)$ is equivalent to the wall separation velocity ($2 U_w$).

2.2. Cahn–Hilliard PF model

We choose a model of a binary mixture with a classical fourth-order polynomial potential $\varPsi$ (see (B1) in Appendix B) (Jacqmin 2000; Carlson 2012). To describe the evolution of both phases, the PF model uses a phase variable $C(x,y)$ ranging from $1$ in the liquid to $-1$ in the vapour. At the interface, the function exhibits a smooth transition. The variable $C$ is governed by a convection–diffusion equation

(2.9)\begin{equation} \frac{\partial C}{\partial t} = \boldsymbol{\nabla} \boldsymbol{\cdot} \left[ M \boldsymbol{\nabla} \phi \right] - \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} C . \end{equation}

Here, $M$ is the PF mobility and $\phi$ is the chemical potential. The latter is defined as

(2.10)\begin{equation} \phi = \frac{2\sqrt{2}}{3}\,\frac{\sigma}{\epsilon}\,\varPsi'(C ) - \frac{2\sqrt{2}}{3}\,\sigma\epsilon \nabla^2 C. \end{equation}

The chemical potential (2.10) contains the surface tension ($\sigma$) and the interface thickness ($\epsilon$). The derivation of (2.9) and (2.10) is standard and can be found in Carlson (2012) and Jacqmin (2000). The surface tension $\sigma$ is a physical parameter and is set to a value corresponding to the water liquid–vapour interface. The constants $M$ and $\epsilon$, on other hand, are typically treated as numerical parameters. Having solved for function $C$, the surface tension force in (2.2) is inserted as

(2.11)\begin{equation} \boldsymbol{{f}}_\sigma = \phi\, \boldsymbol{\nabla} C. \end{equation}

The Cahn–Hilliard equation (2.9) is a fourth-order partial differential equation, thus two boundary conditions are needed on solid walls. The lowest-order boundary condition is

(2.12)\begin{equation} - \mu_f \epsilon \left( \frac{\partial C}{\partial t} + \boldsymbol{{u}} \boldsymbol{\cdot} \boldsymbol{\nabla} C \right) = \frac{2\sqrt{2}}{3}\,\sigma\epsilon\,\boldsymbol{\nabla} C \boldsymbol{\cdot} \hat{n} - \sigma \cos ( \theta_0 )\,g'( C ), \end{equation}

where $\mu _f$ is contact line friction, and $g$ is switch function (B2). Equation (2.12) is also known as the non-equilibrium wetting condition and requires the equilibrium contact angle $\theta _0$ as an input. Setting a non-zero $\mu _f$ yields a dynamic contact angle different from $\theta _0$ (Jacqmin 2000; Qian et al. 2003; Carlson, Do-Quang & Amberg 2011). The second boundary condition on solid walls is $\boldsymbol{\nabla} \phi \boldsymbol {\cdot } \hat {n} = 0$, which states that there is no diffusive flux through the walls.

The full system of fluid and PF equations (2.2), (2.3), (2.9) and (2.10), together with the boundary conditions, is discretized and solved using open-source code FreeFEM (Hecht 2012). Adaptive mesh is used near the two-phase interface to capture the variation of the function $C$. More details about the PF equations and simulations can be found in Appendix B.

The PF model has three unknown free parameters: the PF mobility $M$, the contact line friction $\mu _f$, and the interface thickness $\epsilon$. To provide an intuition about each parameter, we vary them one at a time. As before, $\theta _0 = 95^\circ$ and ${Ca} = 0.20$. Figure 3(a) shows $\Delta x$ for $\mu _f = 0$, $\epsilon = 0.7$ nm and $M = (3.5; 7.0; 10.5; 14.0; 17.5) \times 10^{-16}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$. We observe that the final steady $\Delta x$ value is reduced as $M$ is increased, i.e. smaller deformation with increased diffusion. In figure 3(b), the displacement for $\epsilon = 0.7$ nm, $M = 1.75 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$ and $\mu _f = (0.0; 0.5; 1.0; 1.5; 2.0)\,\mu _\ell$ is shown. Here, $\Delta x$ increases with $\mu _f$, therefore the contact line friction has an opposite effect compared to mobility. This competition has been previously clearly showcased by Yue & Feng (2011). Finally, the evolution of $\Delta x$ in time for $\epsilon = (0.18;0.35;0.70;1.40;2.80 )$ nm is shown in figure 3(c) with $\mu _f = 0$ and $M = 1.08 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$. It can be observed that – in contrast to the $M$ and $\mu _f$ variations – notable differences are present for initial times. As $\epsilon$ is decreased (in the direction against the arrow), the steady $\Delta x$ value is converging. This is a signature of the sharp interface limit (Yue et al. 2010; Xu et al. 2018). For enriched understanding, in Appendix F we report streamlines from PF near a receding contact line for similar parameter variations.

Figure 3. Drop displacement $\Delta x$ over time in PF simulations, varying PF mobility (a), contact line friction (b), and interface thickness (c). In (a), $\mu _f = 0$, $\epsilon = 0.7$ nm and PF mobility is $M = (3.5; 7.0; 10.5; 14.0; 17.5) \times 10^{-16}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$. In (b), $\epsilon = 0.7$ nm, $M = 1.75 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$ and contact line friction is $\mu _f = (0.0; 0.5; 1.0; 1.5; 2.0)\,\mu _\ell$. In (c), $\mu _f = 0$, $M = 1.08 \times 10^{-15}\ {\rm m}^4\ ({\rm N}\ {\rm s})^{-1}$ and interface thickness is $\epsilon = (0.18;0.35;0.70;1.40;2.80 )$ nm. Varying parameters are increasing along the black arrows in all panels. In all simulations, the equilibrium contact angle is $\theta _0 = 95^\circ$ and the wall velocity is $U_w = 6.67\ {\rm m}\ {\rm s}^{-1}$.

3. Molecular dynamics simulations of the sheared droplet

We describe the polar molecules of the water droplet using the SPC/E model. This is the simplest model allowing hydrogen bonds with the solid substrate. It also offers an accurate description of water bulk and interfacial properties, and retains a relatively low computational cost. The bounding walls are formed as mono-layers of ${\rm SiO}_2$ quadrupoles (figure 4a) that are restrained into a hexagonal lattice. A quasi-two-dimensional system with depth $4.68$ nm (figure 4c) is constructed. Albeit the composition of walls is structurally unrealistic, this simple surrogate configuration allows emulating the two fundamental electrostatic interactions characteristic to hydrophilic substrates (Johansson & Hess 2018). The first is the hydrogen bond between water and silica (dotted line in figure 4b). The second interaction is the adsorption of water molecules on the substrate. The adsorption occurs due to the attraction between water oxygen and silicon atoms (dashed line in figure 4b). The strong electrostatic interaction is responsible for a very small hydrodynamic slip at the wall. Note that the current MD configuration does not include any chemical reactions that would occur at a real crystalline or amorphous silica surface.

Figure 4. (a) Molecular geometry of silica quadrupoles. (b) Water molecules adsorbing (bonds sketched with dashed lines) and forming hydrogen bond (sketched with dotted line) with silica substrate. (c) Top view of the contact line region, showing the lattice of silica quadrupoles and water up to $\approx$1.5 nm above the lower periodic boundary (blue line in b).

The strength of the water–substrate interaction can be tuned by adjusting the charge distribution in ${\rm SiO}_2$ (figure 4a). Different interaction strengths will result in different equilibrium contact angles $\theta _0$. We simulate the system via atomistic MD in the NVT ensemble, using well-established force fields and thermostats. Details regarding the physical and numerical simulation set-up can be found in Appendix C.1.

3.1. Equilibration runs

First, we measure the equilibrium properties of the water drop between two static plates and generate thermodynamically consistent initial conditions. This is done through so-called equilibration runs (see Appendix C.2). During the run time of the MD simulation, we collect flow data in regular $0.2\ {\rm nm} \times 0.2\ {\rm nm}$ bins for a time interval 12.5 ps. This yields instantaneous density $\rho ^i(x,y,t)$ and flow velocity $u^i_x(x,y,t)$, $u^i_y(x,y,t)$ data as functions of space and time. After the initial transient, $\rho ^i(x,y,t)$ is averaged over time to reduce noise in the liquid–vapour interface shape. The interface shape is extracted based on this averaged density, i.e. $\rho (x,y)$. The interface position is determined by seeking the location where the liquid density transitions from bulk density to very small vapour density. Examples of density distributions are shown in figures 5(e, f). More details on interface extraction are provided in C.3.

Figure 5. Liquid water density variation near the bottom wall for equilibrium angles (a) $\theta _0 = 127^\circ$, (b) $95^\circ$, (c) $69^\circ$, and (d) $38^\circ$. Dashed black lines illustrate the boundaries of the bins. The bin filled with red shows the selected location of the solid wall, while the bin filled with green shows the first reliable interface measurement. Water density distributions along the green bin for $\theta _0 = 95^\circ$ configuration are shown (e) over the full span of $x$ coordinates, and (f) near the left liquid–vapour interface.

To determine the equilibrium contact angle, we extrapolate the interface angle towards a hydrodynamic wall position. We assume that the wall position is at the centre of the bin coloured in red in figures 5(a–d). The $q$ values are tuned to yield $\theta _0 = 127^\circ$, $95^\circ$, $69^\circ$ and $38^\circ$. This allows us to investigate hydrophobic, neutral and hydrophilic wetting conditions. In parallel to $\theta _0$ extraction, we also identify the first reliable bin for interface shape measurement. This bin is shown in green in figures 5(a–d). The interface angle computed from points closer to the wall exhibits extreme deviation from continuum description (figure 20d). Consequently, for comparisons with CFM predictions (presented later), we always extract the interface shape from MD simulations neglecting the unreliable data points.

In previous investigations, density variations have been observed for Lennard-Jones (L-J) liquids near solid walls and two-phase interfaces (Bugel, Galliéro & Caltagirone 2011; Stephan et al. 2018). We determine the extent of $\rho$ oscillations near the wall in our MD system. The water liquid density distribution along the height of the channel is computed as

(3.1)\begin{equation} \rho_y (y) = \int\limits_{x_l}^{x_r} \rho (x,y)\,{{\rm d}\kern0.06em x}. \end{equation}

Here, the boundaries for integration $x_l$ and $x_r$ are selected for each $\theta _0$ to fall within the liquid phase for all $y$ coordinates. The close-up of $\rho _y (y)$ near the bottom wall is shown in figures 5(a–d). Oscillations in $\rho$ have smaller amplitude and occur over smaller distances than observed typically in L-J systems. The small layering is an outcome of the combination of the SPC/E water model and the ${\rm SiO}_2$ surface model. For the same value of $\theta _0$, the layering would propagate further into the liquid phase if an L-J surface were used instead of ${\rm SiO}_2$.

3.2. Dynamic configuration

Having obtained the initial state from the equilibrium simulations, we turn to the dynamic configuration. Since the configuration is symmetric in a continuum sense, the final $\Delta x$ is obtained as an average between the left ($\Delta x_l$) and right ($\Delta x_r$) interfaces. We determine the drop displacement using the first reliable bins near the top and bottom walls (figure 5, green bins). Since there are no interface data near the hydrodynamic wall, the final steady drop displacement is obtained by extrapolating the interface shape. We use polynomial extrapolation as detailed in Appendix C.4.

Figure 6(a) shows $\Delta x (t)$ for $\theta _0 = 95^\circ$ for different capillary numbers, starting at ${Ca} = 0.05$ and then increased incrementally by $\Delta {Ca} = 0.05$. For all ${Ca}$ up to and including ${Ca} = 0.25$, we observed a stable configuration. The obtained steady drop displacements for $\theta _0 = 95^\circ$ are summarized in the second row of table 1. The table also reports displacement for the other equilibrium contact angles. Different ${Ca}$ values are gradually tested until at least three simulations in a stable regime are gathered. For all $\theta _0$, as expected, we observe that $\Delta x$ increases with ${Ca}$.

Figure 6. Time evolution of $\Delta x$ from MD with $\theta _0 = 95^\circ$ for simulations reaching steady state. In (a), ${Ca}$ increases along the black arrow, starting at ${Ca} = 0.05$ with increments of $0.05$. Windows of ${\pm }5$ nm around the steady mean value are shown for (b) ${Ca} = 0.20$ and (c) ${Ca} = 0.25$. With red arrows we indicate thermal oscillations and motion similar to stick-slip.

Table 1. Overview of MD simulations carried out in this work. For each simulation, ${Ca}$ and steady $\Delta x$ (in nm) are given. If simulation is unstable, then ‘unst.’ is reported instead of a $\Delta x$ value.

$^a$Calibration simulation.

$^{b}$Stick-slip like behaviour.

The largest stable ${Ca} = 0.25$ at $\theta _0=95^\circ$ exhibits oscillations similar to the so-called ‘stick-slip’ behaviour (Orejon, Sefiane & Shanahan 2011; Varma, Roy & Puthenveettil 2021). For a prolonged time, the contact line shows more resistance towards movement (i.e. stick). This period is followed by another in which the contact line exhibits less resistance towards movement (i.e. slip). We show a zoomed view of $\Delta x (t)$ for ${Ca} = 0.25$ in figure 6(c), where the stick-slip behaviour is identified with red arrows. Note that the partial stick-slip effect for $({Ca}, \theta _0) = (0.25,90^\circ )$ is distinct from the oscillations observed for ${Ca} = 0.20$. We show the enlarged view of $\Delta x$ versus time for ${Ca} = 0.20$ in figure 6(b), where the red arrows depict the oscillation magnitude and time scale. We observe that the oscillations at ${Ca} = 0.20$ are much smaller in magnitude and span smaller time scales compared to the stick-slip-like motion (figure 6c). More discussion on the physics behind these oscillations is provided in § 6.2. Similar oscillations with increased magnitude are observed also for $\theta _0 = 127^\circ$ at ${Ca} = 0.90$ and $\theta _0 = 38^\circ$ at ${Ca} = 0.02$. Simulations exhibiting stick-slip oscillations are indicated with a superscript $b$ in table 1.

3.3. Droplet break-up

As ${Ca}$ is increased further, $\Delta x(t)$ measured from MD does not stabilize around some finite value. Instead, $\Delta x(t)$ grows continuously. An example of $\Delta x(t)$ at $({Ca},\theta _0) = (0.30,95^\circ )$ is shown in figure 7. At $10$ ns, the drop is only moderately deformed (figure 7a). At $23$ ns (figure 7b), there are two drops at the top and bottom walls, connected by thin thread of liquid water. Then at around $23.5$ ns break-up occurs, and at $24$ ns we observe two completely separated drops (figure 7c). Note that at the break-up instant, the slope of $\Delta x(t)$ changes distinctively. This is due to the absence of the surface tension force that was resisting the displacement. For the two separate drops, there is no competition between the friction at the top and bottom contact lines. Instead, the friction at the contact lines now ensures that the two drops follow the wall velocity and separate with speed $2 U_w$. Hence the critical capillary number ${Ca}_c$ lies in between ${Ca} = 0.25$ and ${Ca} = 0.30$. In the sheared droplet configuration, the exact value of ${Ca}_{c}$ is determined by the most unstable contact line (either advancing or receding). For investigations of individual contact lines, other types of experiments should be performed, such as plunging/withdrawn plate (Eggers 2005) or hydrodynamic assist (Afkhami et al. 2018; Liu, Carvalho & Kumar 2019a; Fullana, Zaleski & Popinet 2020).

Figure 7. Drop displacement in MD with $\theta _0 = 95^\circ$ at ${Ca} = 0.30$. In insets (a–c), we show the drop shape at three selected time instances. The selected time is shown by a green circle with an arrow pointing to the inset. The green dashed line corresponds to the relative wall speed $2 U_w$.

4. Calibrating CFM models against MD

The aim of this section is to identify the free parameters in PF and VOF such that the displacement $\Delta x$ obtained from the continuum models matches the displacement obtained from MD simulations. The four calibration pairs $({Ca},\theta _0)$ are (marked with superscript $a$ in table 1) $(0.60,127^\circ )$, $(0.20,95^\circ )$, $(0.10,69^\circ )$ and $(0.015,38^\circ )$. For each $\theta _0$, we have chosen the largest steady ${Ca}$ value available from MD simulations. If the chosen steady ${Ca}$ value leads to large stick-slip like oscillations of $\Delta x$ (such as observed in figure 6c), then we select the previous (smaller) steady ${Ca}$ value. The unsteady simulations are not suitable for calibration, because $\Delta x(t)$ grows in time (figure 7).

Before calibration, we have to make sure that the same system – in terms of geometry and fluid properties – is represented in CFM and MD. The dimensions and physical properties are reported in figure 1(a). The bulk liquid density is obtained by taking an average of $\rho _y$ over the height of the channel:

(4.1)\begin{equation} \rho_\ell = \int_{y_b}^{y_t} \rho_y (y)\, {{\rm d} y} \approx 990\ \mathrm{kg}\ \mathrm{m}^{{-}3}. \end{equation}

This density is valid for all equilibrium angles. Viscosity and surface tension of liquid SPC/E water are taken from previous work (Lācis et al. 2020). The viscosity measurement of the vapour SPC/E phase from MD is impractical. Therefore, both viscosity and density of the vapour are determined from engineering tables (Appendix D). For VOF simulations, we increase the vapour density to $\rho _v^{VOF} = 9.9\ {\rm kg}\ {\rm m}^{-3}$. This ensures numerical stability of the simulations, while keeping the influence of the vapour inertia negligible.

The size of the channel in the $x$ direction (the distance between periodic boundary conditions) is matched with the distance between the left and right periodic boundary conditions in the MD system. The channel height is based on the chosen hydrodynamic wall position at the centre of the red bin (figures 5a–d), which results in $H = 29.22$ nm. This position is the subject of an investigation in itself (Herrero et al. 2019) and is discussed further in Appendix E.

The slip length (2.4) used in CFM models in principle can be related directly to MD simulations. Previous work (Huang et al. 2008) has demonstrated that the slip length of the MD system follows a quasi-universal relationship with respect to equilibrium contact angle $\theta _0$. However, for the chosen MD system (SPC/E water on the ${\rm SiO}_2$ surrogate wall), accurate slip length quantification has not yet been done. To obtain an indication about the validity of $\ell _s = 0$ for the selected MD system, we compare streamwise velocities near the wall between PF and MD (Appendix E). Through this comparison, we found that $\ell _s = 0$ holds for $\theta _0 = 38^\circ \text {--} 95^\circ$, while for $\theta _0 = 127^\circ$, the appropriate choice is $\ell _s = 0.44$ nm. These values are used in both VOF and PF as input parameters without any additional fitting; see the sixth and ninth columns of table 2. Note that $\ell _s$ in VOF is fixed and distinct from $\lambda$ – the length scale in condition (2.8), which is used for calibration.

Table 2. Results of CFM calibration against MD. For each $({Ca},\theta _0)$ pair, we report MD reference data (steady displacement $\Delta x$) in the third column. In the following four columns, we report VOF parameters (2.4) and (2.8), and resulting steady displacements. In the final five columns, PF parameters (2.4), (2.9) and (2.12), and resulting displacements, are given.

4.1. VOF calibration

Figure 8 compares $\Delta x$ obtained from VOF (red) and MD (black) for the four calibration configurations. We have adjusted the parameter $\lambda$ in (2.8) to find the best fit of $\Delta x$ as $t\rightarrow \infty$. We observe that the droplet displacement from VOF stabilizes at values that agree well with the atomistic simulations, with the most challenging calibration configuration being $\theta _0=127^\circ$ (figure 8a). The $\lambda$ values that reproduce the displacement obtained from MD are listed in the sixth column of table 2. For the hydrophobic configurations, we obtained $\lambda = 4.325$ nm for $\theta _0 = 127^\circ$ and $\lambda = 0.935$ nm for $\theta _0 = 95^\circ$. We observe that for $\theta _0 = 95^\circ$, $\lambda$ is roughly four times smaller than for $\theta _0 = 127^\circ$. Naively, $\lambda$ can be regarded as a slip-related length scale near the contact line. Larger $\lambda$ for $\theta _0 = 127^\circ$ then suggest smaller friction, in line with the findings from PF calibration (§ 4.2).

Figure 8. Time evolution of $\Delta x$ in MD, PF and VOF for all calibration configurations.

For $\theta _0 = 69^\circ$ and $38^\circ$, we found even smaller $\lambda$ values, $0.313$ nm and $0.146$ nm, respectively. Qualitatively, this continuous behaviour of decaying $\lambda$ for smaller $\theta _0$ follows the contact line friction argument presented later (§ 4.2). It is also interesting to note that the obtained $\lambda$ values from this calibration procedure are of order 0.1–4  nm. This is similar to what is used by Legendre & Maglio (2015), who set $\lambda = 1$ nm in macroscopic simulations. The mesh spacing in the VOF simulations is $\varDelta /2 = 0.457$ nm. For $\theta _0 = 38^\circ$ and $69^\circ$, we have $\lambda < \varDelta /2$ (table 2). On the other hand, for $\theta _0 = 95^\circ$ and $127^\circ$, we have $\lambda > \varDelta /2$ (table 2). Here, the sign of the logarithm in (2.8) is negative. This is typically not the case when applying Cox inspired relationships such as (2.8). Therefore, imposed $\theta _{num}$ for $\theta _0 = 95^\circ$ and $127^\circ$ should be viewed as a numerical parameter to tweak the curvature near the wall.

4.2. PF calibration

The droplet displacement $\Delta x$ produced by PF (green) and MD (black) are compared in figure 8 for the calibration configurations. The interface thickness, mobility and line friction that provide the best match with the displacement obtained from MD are listed in table 2. For PF, different combinations of parameters may provide a good fit; therefore it is prudent to have physically motivated guidelines.

4.2.1. PF calibration proposed in the literature

For the Cahn–Hilliard PF model, there is a standard calibration procedure proposed by Yue et al. (2010) and Yue & Feng (2011). The sequential steps are: (i) choosing the interface thickness $\epsilon$ that is suitable to describe the physical problem; (ii) setting the PF mobility $M$ according to the sharp interface limit (Yue et al. 2010); and (iii) calibrating the contact line friction $\mu _f$ against experiments.

The interface thickness has to be smaller than the important physical length scales in the chosen system, which for the sheared droplet configuration are the water drop height $H = 29.22$ nm and width $W \approx 38$ nm (figure 1a). Based on this, the interface thickness is set to $\epsilon = 0.7$ nm. According to the sharp interface limit (Yue et al. 2010), the criterion for choosing $M$ is

(4.2)\begin{equation} M > 1/16 \epsilon^2 / \sqrt{\mu_v \mu_\ell} . \end{equation}

Inserting the chosen values for interface thickness, vapour viscosity and liquid viscosity, we obtain $M > 3.21 \times 10^{-16}\ {\rm m}^4\ {\rm N}^{-1}\ {\rm s}^{-1}$. In the MD simulations, we do not observe any physical effect that would hint towards a large $M$ value. Therefore, for all $\theta _0$ in this section, we select $M = 3.5 \times 10^{-16}\ {\rm m}^4\ {\rm N}^{-1} \ {\rm s}^{-1}$, which is close to the lower limit.

The last step is to carry out simulations with different $\mu _f$ until a steady state is attained that matches the displacement obtained from MD. Table 3 summarizes the displacement obtained for each calibration pair $({Ca},\theta _0)$ at $\mu _f = 0, 0.1, 0.2, 0.3, 2.36,11.84$. The fourth column of table 3 shows the displacements for $\mu _f = 0$. We observe that no steady solution has been obtained for $\theta _0 = 127^\circ$. The steady $\Delta x$ obtained for $\theta _0=95^\circ$ overestimates the MD value (third column of table 3) by roughly $50\,\%$. Increasing $\mu _f$ only deteriorates the agreement. Therefore, we conclude that the matching procedure proposed by Yue & Feng (2011) is not adequate to calibrate the PF model for the chosen nanoscale configuration at $\theta _0 = 127^\circ$ and $95^\circ$.

Table 3. Calibrating PF with MD by changing $\mu _f$ (strategy proposed by Yue & Feng 2011). In the third column, we show the reference $\Delta x$ from MD. To the right, we show steady $\Delta x$ (in nm) obtained from PF for specified $\mu _f / \mu _\ell$. If no steady state exists, then we write ‘unst.’. In all simulations, $M = 3.5 \times 10^{-16}\ {\rm m}^4\ {\rm N}^{-1}\ {\rm s}^{-1}$ and $\epsilon = 0.7$ nm.

For $\theta _0 = 69^\circ$ and $38^\circ$, on the other hand, $\mu _f = 0$ results in an underestimated drop displacement (table 3) compared to MD. We observe that the required PF contact line friction for $\theta _0 = 38^\circ$ ($\mu _f = 11.84\mu _\ell$) is roughly five times larger than the one for $\theta _0 = 69^\circ$ ($\mu _f = 2.36\mu _\ell$). More hydrophilic $\theta _0$ entails a stronger affinity between the liquid and the wall. Stronger affinity, in turn, can yield larger friction, which is consistent with the obtained $\mu _f$ values.

4.2.2. Calibrating PF by adjusting mobility for $\theta _0 = 127^\circ$ and $95^\circ$

A way to reduce friction near the contact line is to increase $M$ (figure 3a) and thus allow for more diffusion. Therefore, we impose $\mu _f = 0$ and increase $M$ for $\theta _0 = 127^\circ$ and $95^\circ$ until the steady $\Delta x$ from the PF agrees with MD. The required PF mobility values are reported in table 2. We observe that for $\theta _0 = 95^\circ$, a three times larger $M$ is required compared to the $M$ set by the sharp interface limit (see the tenth column of table 2), whereas for $\theta _0 = 127^\circ$, the $M$ value has to be increased by a factor of hundred.

Note that this calibration procedure is compatible with the sharp interface limit. The condition (4.2) can be rewritten in terms of interface thickness as $\epsilon < 4\,M^{1/2}\mu _v^{1/4}\mu _\ell ^{1/4}$. As we increased $M$, the interface thickness $\epsilon = 0.7$ nm was kept constant, therefore the condition is satisfied. However, this calibration procedure produces a noticeable diffusion near the contact line (figure 23c), which is not observed in the MD simulations.

5. Predictions from PF and VOF models

We have shown that continuum systems can be tuned to match the final steady droplet displacement computed from MD simulations. It is also necessary to understand how well the CFM models capture other key features of the system, including the interface shape and the time-dependent transient behaviour of the droplet. Moreover, an important practical aspect is the accuracy of the CFM when it comes to predicting the droplet behaviour to characterize the sheared droplet system quantitatively for a range of capillary numbers. Specifically, in this section, we fix the parameters for PF and VOF to values reported in table 2.

5.1. Time evolution of drop displacement

Figure 9 shows the droplet displacement as a function of time for a fixed $\theta _0=95^\circ$ but different ${Ca}$. Figure 9(d) is identical to figure 8(b) and corresponds to the conditions for which the system was calibrated, i.e. ${Ca}=0.20$. We observe that both PF and VOF capture the transient dynamics very well overall. The PF model is slightly slower (predicts smaller $\Delta x$ at the same time instant) than VOF and MD.

Figure 9. Time evolution of $\Delta x$ in MD, PF and VOF for ${Ca} \in (0.05, 0.30)$, with equilibrium contact angle $\theta _0 = 95^\circ$. The calibrated functions $\Delta x (t)$ from PF and VOF are shown in (d) with green and red lines, respectively. The a priori measurements of $\Delta x (t)$ from MD (§ 3) are shown with a black line. This colour code is retained for all comparisons that follow.

Figures 9(a,b,c,e) show the transient behaviour of CFM models in off-calibration conditions for ${Ca} = 0.05$, $0.10$, $0.15$ and $0.25$. For moderate capillary numbers ${Ca} \geqslant 0.10$, we observe qualitatively similar behaviour as for ${Ca} = 0.20$. The CFM models predict the transient and steady $\Delta x$ values rather accurately. The PF, however, is always slightly slower compared to the VOF model. For ${Ca} = 0.05$, the steady $\Delta x$ value is slightly lower in CFM models than in MD. The PF model is slower in the transient compared to the VOF, and the agreement with MD is arguably worse.

To conclude the $\theta _0 = 95^\circ$ investigations, we consider the CFM model predictions for the unsteady configuration with ${Ca} = 0.30$ (figure 9f). The results from PF, VOF and MD are indistinguishable for the first $5$ ns, showing excellent predictive capability. For later times, VOF simulation over-predicts and PF simulation under-predicts the $\Delta x$ observed from MD. This is in line with observations in steady situations (figures 9a–e). Both VOF and PF exhibit the rapid change of slope at around $23.5$ ns, corresponding to the drop break-up (discussed in § 3).

We have carried out the same investigation for $\theta _0=38^\circ, 95^\circ$ and $127^\circ$. Qualitatively similar results to those observed in figure 9 are obtained; see supplementary figures 1–3 available at https://doi.org/10.1017/jfm.2022.219.

5.2. Interface shape

For interface shape comparisons, the data are presented as the variation of the angle along the interface $\theta ( y )$ (figure 2a). Figure 10 compares $\theta (y)$ obtained from MD, PF and VOF for $\theta _0=95^\circ$ and different ${Ca}$. The equilibrium angles at the bottom and top walls are shown with black dotted lines, while the hydrodynamic wall positions are presented with black dashed lines. The interface shapes in both calibration (${Ca}=0.20$) and in off-calibration conditions (${Ca} = 0.05$, $0.10$, $0.15$ and $0.25$) are similar between the models. We observe that both PF and VOF exhibit more pronounced differences from MD at the top wall (near the advancing contact line). At ${Ca} = 0.25$ (figure 10e), the CFM model predictions are notably different from MD data compared to the other capillary numbers. The loss of accuracy could arise from the large $\Delta x$ oscillations in MD at ${Ca} = 0.25$ (§§ 3 and 6.2) that are not captured with the CFM models.

Figure 10. Steady interface shape from MD, PF, and VOF simulations, $\theta _0 = 95^\circ$. Equilibrium angles and solid wall locations are shown with black dotted and dashed lines, respectively.

The MD data in figure 10 do not extend all the way to the wall, so the exact value of the dynamic contact angle is not known. At $\theta _0=95^\circ$, the dynamic contact angle in PF simulations is equal to the equilibrium angle since $\mu _f = 0$ (table 2). For VOF, the dynamic contact angle is different from $\theta _0$. However, the angle is larger (smaller) than the equilibrium angle for the receding (advancing) contact line. This is opposite to the understanding arising from analysis of uncompensated Young stress. The source of this effect is the value $\lambda = 0.935 \ {\rm nm} > 0.467\ {\rm nm} = \varDelta / 2$ (table 2).

We have repeated this comparison of interface shapes at different ${Ca}$ values between MD, PF and VOF for equilibrium contact angles $\theta _0 = 127^\circ$, $69^\circ$ and $38^\circ$. The results are reported in supplementary figures 4–6 and are similar to what is presented above.

5.3. Steady drop displacement

Figure 11 shows the steady displacement as function of $Ca$ for different $\theta _0$. As expected, the MD, PF and VOF points collapse for the calibration configurations (marked with vertical arrows). For the $\theta _0 = 127^\circ$ configuration, the PF and VOF models overestimate the steady displacement for large ${Ca}$. A possible reason behind the discrepancy between PF, VOF and MD could be the slip length $\ell _s$. It has been determined for ${Ca} = 0.60$ and kept constant for runs with different ${Ca}$. In general, the slip length can increase for larger wall shear stress (Thompson & Troian 1997). This, consequently, would reduce $\Delta x$ and possibly move the PF and VOF predictions closer to the MD results. For $\theta _0 = 38^\circ$, the agreement between PF, VOF and MD diverges as we increase ${Ca} > 0.02$. A reason for the disagreement could be the increased contact line friction (§ 6.1) at receding contact line, from which a liquid film is formed (§ 5.5 and figure 14). This asymmetry is not taken into account in the CFM simulations.

Figure 11. Steady $\Delta x$ as a function of ${Ca}$ for all $\theta _0$.

5.4. Critical capillary number

Figure 12 shows the critical capillary number, ${Ca}_c$, as a function of $\theta _0$ from the three methods. To determine ${Ca}_c$ for a given $\theta _0$, we take the mean of the largest steady (${Ca}_s$) and smallest unsteady (${Ca}_u$) capillary numbers that were simulated, i.e.

(5.1)\begin{equation} {Ca}_c = \frac{{Ca}_u + {Ca}_s}{2} \pm \frac{{Ca}_u - {Ca}_s}{2}. \end{equation}

The uncertainty is determined as half of the difference between these capillary numbers. To reduce $\Delta {Ca}_c$ in CFM, we sample the ${Ca}$ space with smaller intervals. It can be observed from figure 12 that ${Ca}_c$ increases for larger $\theta _0$, which has also been reported in previous numerical (Sbragaglia et al. 2008) and analytical (Hocking 2001; Eggers 2004) investigations.

Figure 12. Critical capillary number ${Ca}_c$ as a function of $\theta _0$.

We did not manage to run any VOF simulations for $\theta _0 = 127^\circ$ and ${Ca} \geqslant 0.80$. Due to a large $\lambda$ (table 2), the dynamic angle for ${Ca} \geqslant 0.80$ falls outside the physically admissible range $(0^\circ, 180^\circ )$. Consequently, in figure 12 we compare only PF predictions with the MD data at $\theta _0 = 127^\circ$. The PF model slightly under-predicts the ${Ca}_c$ value. This discrepancy could again be due to the uncertainty in the slip length $\ell _s$.

Figure 12 shows that predictions of ${Ca}_c$ for $95^\circ$ and $69^\circ$ agree with MD results within the accuracy bounds of MD. For $\theta _0=38^\circ$, however, discrepancies are observed, where ${Ca}_c$ values computed from both CFM models are larger compared to the MD results. Note that the differences at $38^\circ$ are enhanced due to the logarithmic scale of the ${Ca}_c$ axis in figure 12.

5.5. Above the critical capillary number

In this subsection, we investigate the accuracy of the CFM models for predicting the unsteady breakage of the droplet. This constitutes the most challenging test of the CFM models in off-calibration conditions. Two configurations are selected with ${Ca} > {Ca}_c$, namely, $\theta _0 = 95^{\circ }$ at ${Ca} = 0.30$, and $\theta _0 = 38^{\circ }$ at ${Ca} = 0.05$.

Figure 13 shows the drop shape at four time instances for $({Ca},\theta _0)=(0.30,95^{\circ })$. We observe that the three models provide similar deformed states at $t \approx 5$ ns and $t \approx 10$ ns; see the top two rows in figure 13. The time instance just before and right after the break-up is shown in the third and fourth rows of figure 13. The thread connecting the lower and upper drops is very thin in the MD and VOF simulations, compared to PF. In addition, the thread in VOF is comparably long and exhibits grid-to-grid like oscillations. There are also pronounced differences in the neck of each satellite drop – the region in which the thread transitions to the drop shape. The PF simulations show the thickest neck (figure 13g), followed by MD with a slightly thinner neck (figure 13c) and VOF with a very small drop neck (figure 13k). The time instants at which the break-up occurs are remarkably similar: $t^{MD}_b = 23.41$ ns for MD, $t^{PF}_b = 24.37$ ns for PF, and $t^{VF}_b = 23.23$ ns for VOF. For a complete time-dependent animation of MD, PF and VOF simulations side by side, see supplementary movie 1.

Figure 13. Snapshots of interface shape evolution over time from MD (a–d), PF (e–h) and VOF (i–l), with equilibrium contact angle $\theta _0 = 95^\circ$ and capillary number ${Ca} = 0.30$.

We turn our attention to $\theta _0 = 38^{\circ }$ and ${Ca} = 0.05$. To avoid a premature break-up of the liquid bridge in the PF simulation, we use $\epsilon = 0.35$ nm for this particular simulation since in this unsteady example, a third length scale – width of the liquid bridge in the deformed state – becomes important. Droplet shapes at five time instances from MD, PF and VOF are shown in figure 14. We observe that initially (top two rows in figure 14) the drop is deforming only slightly. After a certain time (figures 14c,h,m), a liquid film is generated at the receding contact lines of the drop. Very similar drop shapes are observed in MD, PF and VOF at different time instances. The VOF simulation takes roughly $5$ ns less to form a film similar to one observed in MD. For PF, one has to wait for around $15$ ns more than in MD. As time progresses, the liquid film becomes longer (figures 14d,i,n), and the liquid bridge becomes thinner. We also note that the tip of the film forms a thicker drop-like region. Finally, the liquid film from the periodic image merges with the liquid bridge (figures 14e,j,o). For MD and VOF simulations, the coalescence occurs only near one of the walls. For PF, on the other hand, a symmetric configuration is obtained with fully wetted top and bottom walls. Animations of the simulations can be found in supplementary movie 2.

Figure 14. Snapshots of interface shape evolution over time from MD (a–e), PF ( f–j), and VOF (k–o), with equilibrium contact angle $\theta _0 = 38^\circ$ and capillary number ${Ca} = 0.05$. For PF, we have $\epsilon = 0.35$ nm.

6. Molecular physics of the sheared droplet

In this section, we present molecular phenomena of the sheared droplet that are particularly challenging to deal with in a continuum model. First, we present contact line friction measurements directly from MD, and asymmetry between advancing and receding lines for hydrophilic and hydrophobic configurations. Second, we discuss the nature of the stick-slip like oscillations. Finally, we show the inevitable three-dimensional nature of the drop break-up.

6.1. Contact line friction measurements from MD

To extract the $\mu _f$ from MD, we use the formula proposed by Yue & Feng (2011):

(6.1)\begin{equation} \left[\frac{\sqrt{2}}{3}\,\frac{\mu_f}{\mu_\ell}\sin\theta\right]{Ca} = \cos\theta_0 - \cos\theta . \end{equation}

Here, $\theta$ is the dynamic contact angle at the wall. This expression is an approximation of the wetting boundary condition (2.12) in case of no-slip and small capillary number. Note that there are alternative approaches to determine contact line friction, for example, based on equilibrium simulations, as proposed by Fernández-Toledano, Blake & De Coninck (2019) and Fernández-Toledano et al. (2020). For this work, however, we have determined that a non-equilibrium approach based on fitting (6.1) to MD data is the most efficient approach.

To determine the $\mu _f$ from (6.1), the dynamic contact angle $\theta$ has to be determined first for each ${Ca}$ value. We extract the dynamic contact angle above the first reliable bin (green bins in figures 5a–d). To reduce the noise, we use polynomial interpolation of the interface shape to read the dynamic contact angle at the chosen location. For consistency, we also re-evaluate $\theta _0$ at the same distance from the wall. The obtained dynamic contact angles for all $\theta _0$ values are gathered in figure 15. For each $\theta$, we also display error bars, obtained by adding $\pm$2 to the polynomial order of the interpolation. As expected, when ${Ca}$ increases, so does the deviation of $\theta$ from $\theta _0$.

Figure 15. Measured dynamic contact angles $\theta$ for different ${Ca}$ along with uncertainty intervals for advancing (red squares) and receding (blue rhombus) contact lines. The dashed lines are the best least-squares fits of (6.1). The equilibrium contact angles are (a) $\theta _0 = 127^\circ$, (b) $95^\circ$, (c) $69^\circ$, and (d) $38^\circ$.

We use a least-squares fit to match (6.1) to $\theta$ measurements. The fit provides $\mu _f/\mu _\ell$ values and error intervals. We fit the measurements taken at the top-left/bottom-right and top-right/bottom-left contact lines separately. This allows the observation of an asymmetric line friction. The obtained best-fit lines for all equilibrium contact lines are shown in figure 15. The contact line frictions are listed in table 4 along with previously reported $\mu _f$ from calibration of PF simulations.

Table 4. Comparison of contact line friction used in PF simulations and values obtained directly from MD results. For PF, the same contact line friction is used for advancing and receding contact lines, while in MD, the contact line friction can be different.

We observe that larger line friction parameters are measured for smaller $\theta _0$. This observation is consistent with the molecular kinetic theory (MKT). It states that the line friction scales with the work of adhesion needed to desorb water molecules from the substrate, which in turn increases as the surface becomes more hydrophilic (Blake & Haynes 1969). Interestingly, we observe a difference between advancing and receding line friction for $\theta _0 = 127^\circ$ and $38^\circ$. Future dedicated work is required to investigate this asymmetry in depth, in particular, to determine if the asymmetry depends on the position of the hydrodynamic wall. The observation of asymmetry is not specific to (6.1). A linearized MKT model – where line friction is defined in a slightly different way – would not change the conclusions. The line frictions obtained by fitting expression (6.1) directly to MD data are larger than the ones obtained through calibration of PF against MD (table 4). This fact seems to entail some missing physical effects in the PF model. Ideally – in a potentially more advanced PF model – one would employ $\mu _f$ obtained from MD directly in the PF boundary condition (2.12).

6.2. Stick-slip like oscillations

To differentiate between fluctuations caused by molecular-scale motion and large stick-slip like oscillations, we define three reference length scales. The characteristic length scale of interface fluctuations far from contact lines can be estimated by balancing the thermal energy $k_B T$ with the energy due to surface tension $\sigma l_{th}^2$, giving $l_{th}=\sqrt {k_B T/\sigma }\simeq 0.27$ nm. The other two scales are the van der Waals diameter of SPC/E water molecules $\sigma _{SPC/E} \simeq 0.32$ nm, and the hexagonal lattice spacing of the substrate $d_{hex}=0.45$ nm. We then compute the root-mean-square for fluctuations of the contact line displacement ${\rm RMSF}=(\langle \Delta x^2\rangle -\langle {\Delta x}\rangle ^2)^{1/2}$ for each of the four contact lines, both at equilibrium and under shear conditions. For $\theta _0 = 95^\circ$, we show the obtained RMSF values as a function of ${Ca}$ in figure 16(a).

Figure 16. (a) Root-mean-square fluctuations (RMSF) of all four contact lines against capillary number for $\theta _0=95^\circ$. The error bars show a conservative estimate of the standard error. The spacing in the hexagonal silica lattice (blue dashed line), the van der Waals (vdW) diameter of the SPC/E water molecules (red dot-dashed line), and the length scale of thermal fluctuations (green dotted line) are shown. (b) Drop displacement over time for $\theta _0=95^\circ$ and ${Ca}=0.25$. The mean ($\Delta x$), left ($\Delta x_l$) and right ($\Delta x_r$) displacements are shown. A close-up on a time interval showing pinning/depinning is shown in the inset. The green dotted lines show the slope corresponding to wall velocity: at the beginning of the simulation, both advancing and receding contact lines stick to the wall and thus move apart from each other with velocity $2U_w$; during stick-slip events, the receding contact line sticks to the wall while the advancing line maintains a steady motion, thus the displacement matches the velocity of a single wall $U_w$.

We observe that for ${Ca}<0.25$, the RMSF is comparable with $d_{hex}$, hinting that the observed fluctuations for moderate capillary numbers are caused by the same process producing fluctuations at the equilibrium, that is, the local thermal-induced pinning (de-pinning) of the contact line on (from) adsorption sites close to the average contact line position. It is worth noticing that the scale of these fluctuations is larger than that expected on the interface far from contact lines. This observation is confirmed by examining the RMSF across the whole interface at equilibrium (see supplementary figures 7–10). Note that for some contact angles, the fluctuations in equilibrium simulations are smaller compared to $d_{hex}$.

We examine more closely the MD simulation with $\theta _0=95^\circ$ at ${Ca}=0.25$ (figure 6a), which shows a much larger contact line RMSF, incompatible with lattice spacing driven fluctuations. The speed of the drop displacement during the stick-slip like motion is much smaller than $U_w$ (figure 16b). We conclude that the contact line does not entirely pin when resisting motion and only partially slips when complying with it. Moreover, for most of the time evolution, $\Delta x_l$ and $\Delta x_r$ are synchronized in an oscillatory motion. However, there are a few time intervals (between 17 and 20 ns, and between 45 and 47 ns) where indeed complete stick-slip occurs. In these intervals, the advancing and receding speeds match the magnitude of wall velocity; see the inset of figure 16(b).

It appears that local (pinning/depinning) and global (oscillations) processes co-exist. Pinning/de-pinning can be explained by the fact that ${Ca}$ is close to the critical capillary number ${Ca}_c$. The physical interpretation of the global oscillations is more delicate. In our modelling approach, we have assumed implicitly that the contact line motion is over-damped. This means that there should be a (possibly nonlinear) direct relation between force and speed, in which the proportionality constant is the contact line friction. This may not be true for large wall velocities. In such a case, the effects of the neighbouring flow inertia come into play. In this scenario, the coupling between positions and forces is more complex. The stochastic forcing produced by thermal fluctuations of the microscopic contact angle is no longer completely dissipated. Instead, these may excite oscillation modes of the whole interface. Stick-slip has been theorized for flat surfaces and homogeneous fluids under some flow conditions (Hocking 2001; Eggers 2005; Varma et al. 2021). However, to the best of the authors’ knowledge, it has not been observed directly.

The selected standard CFM models are not capable of describing stick-slip like oscillations. The Navier–Stokes equations, underlying the PF and VOF models, are inherently deterministic, where all the thermal oscillations are averaged out. To model stick-slip like processes in a continuum model, a possible approach could be to introduce random fluctuations on the imposed $\theta _0$. The distribution of the contact angle oscillations has been identified previously (Smith et al. 2016). Fluctuations of the contact angle would in turn induce oscillations in $\Delta x (t)$.

6.3. Three-dimensional nature of drop break-up

In figures 13(d,h), we have observed that the interface shapes obtained from MD and PF display more diffuse regions at the tip of the broken thread. It is tempting to conclude that PF captures the MD behaviour correctly. However, the MD snapshots in figure 13 have been averaged in the spanwise $z$ direction, while PF simulation is a purely two-dimensional (2-D) result.

To obtain a better physical understanding of the exact break-up process, we investigate the molecular field of the MD system exactly at the break-up. In figure 17, we show water molecule locations for $\theta _0 = 95^\circ$ at ${Ca} = 0.4$ shortly before the break-up. Observing the drop from the side (figure 17a), it seems that the thread is not interrupted. In reality, the thin thread develops three-dimensional (3-D) holes (figure 17b) before disconnecting completely. Averaging in the $z$ direction results in a lower density (than the bulk), thus giving the impression of a diffuse interface (figure 13e).

Figure 17. Detailed views of the molecular system upon breakage for $\theta _0=95^\circ$ and ${Ca}=0.4$, observed at $t=9.45$ ns: (a) overview from the side, and (b) close-up on the neck region. Two transparent periodic images are added at the sides. The close-up (b) is obtained by positioning the camera orthogonal to the thin thread in (a).

It is also worth noticing that the formation of these 3-D threads occurs very quickly. The time that it takes for the threads to form is around $50$ ps before the actual break-up. The break-up itself occurs after a few nanoseconds from the instant the 2-D neck starts to form. We can thus infer that for the chosen combination of a molecular system and domain size, there exists a time scale separation between 2-D and 3-D breakage dynamics. This in turn suggests that the selected MD system dimensions are reasonable to produce 2-D results and can be compared directly with 2-D CFM simulations.

7. Discussion

Based on the data presented in the previous sections, we discuss the accuracy, limitations and future development directions of CFM models when it comes to modelling the molecular physics in a nanoscale channel.

7.1. Challenges of the chosen CFM models

For $\theta _0 = 127^\circ$, we observed differences in the steady drop displacement for ${Ca}>0.6$; see figure 11(a). In general, very high contact angles are challenging to model using the particular CFM parameters that we have chosen for calibration. For VOF, it is important to consider condition (2.8) as a numerical model that allows adjustment of the interface curvature at the contact line. For $\theta _0 = 127^\circ$, this was achieved by increasing $\lambda$ to values larger than the grid size, i.e. $\lambda > \varDelta /2$ (table 2). This leads to a negative logarithm in condition (2.8). This means that for the bottom-left receding contact line (figure 2a), where ${Ca}_{cl}$ is negative, we have $G(\theta _{num}) > G(\theta _0)$ and thus $\theta _{num} > \theta _0$. In other words, the imposed numerical contact angle is increased such that the curvature at the receding contact line results in a sufficiently large positive forcing in the $x$ direction on the fluid. The imposed numerical contact angle therefore does not approximate the true dynamical contact angle of the system. For PF, the mobility parameter $M$ is calibrated (§ 4.2.2) to increase contact line velocity for high $\theta _0$. Here, mobility is considered as a numerical tuning parameter to match the droplet displacement and does not correspond to the actual molecular diffusion at the interface (see Appendix F and figures 23a–c).

Very low contact angles impose other challenges for CFM. For $\theta _0 = 38^\circ$, we observed that steady drop displacement diverged (figure 11d) between the three models. Indeed, the prediction of a critical capillary number showed deviations between the methods (figure 12). In addition, the time that it takes for the drop to evolve in different shapes (figure 14) is different in all simulation methods. This is despite the fact that the $\lambda$ obtained through VOF calibration (§ 4.1) is reasonable ($\lambda < \varDelta / 2$) and the standard PF calibration procedure (§ 4.2.1) works well. By investigating MD data directly, asymmetric contact line friction was observed (table 4). For $\theta _0 = 38^\circ$, we observed significantly larger $\mu _f$ for the receding contact line compared to that for the advancing contact line. In the CFM models, on the other hand, the contact line properties ($\mu _f$ for PF and $\lambda$ for VOF) were the same for both advancing and receding sides. Since the receding contact line becomes unstable first (figure 14), this is the likely reason for the CFM model inaccuracy.

Finally, there are limitations of the chosen CFM at very low and high ${Ca}$ values. At high ${Ca}$ values (close to ${Ca}_c$), we observed enhanced oscillations of $\Delta x (t)$ (figures 6(c) and 16, superscript $b$ in table 1). More detailed analyses of these oscillations were presented in § 6.2 for $\theta _0 = 95^\circ$. The CFM model does not include the intrinsic oscillations present in the molecular reality. Consequently, the agreement between CFM and MD in drop time evolution (figure 9(e) and supplementary figure 6) is degraded. This could also be a potential source of increased discrepancy between interface shapes from CFM and MD (figure 10e). Also, for low ${Ca}$ values (for example, ${Ca} = 0.05$ in figure 9a), we have observed a relatively large difference between CFM models and MD results. The underlying cause for this inaccuracy is the relatively small drop displacement. As seen in figure 6(a), the MD oscillations have roughly the same magnitude for all ${Ca}$ values. Consequently, the signal-to-noise ratio in MD is much smaller for smaller ${Ca}$ values, and the large oscillations can give the impression of larger inaccuracy in CFM models.

7.2. Fluid slippage and contact line friction

In this work, the nanoscale molecular system has a negligible hydrodynamic slip. It was observed that the slip length $\ell _s$ exhibits only small variations in $\theta _0$. Below bulk liquid, the slip length had to be $0.44$ nm for $\theta _0 = 127^\circ$, and $0$ nm for all other $\theta _0$. As discussed in § 7.1, for low-friction configurations, $\theta _0>90^\circ$, it was necessary to adjust $\lambda$ and $M$ to further enhance contact line movement relative to the wall. This serves as a hint that the friction near the contact line is much smaller compared to what is modelled through bulk slip length $\ell _s$. In general, however, other physical mechanisms that are intrinsic to the contact line are expected to co-exist with molecular slippage.

We have extracted contact line friction directly from MD (table 4). For hydrophilic substrates, a particularly high friction value was obtained. This is consistent with the formation of a microscopic water film for a relatively small capillary number (i.e. $\theta _{rec}\sim 0$ for ${Ca}\ll 1$). On the other hand, the interpretation for the almost vanishing receding friction on the $\theta _0 = 127^\circ$ surface is less obvious. To the best of the authors’ knowledge, nowhere before has an asymmetric behaviour been reported. It is tempting to explain the asymmetry by stating that hydrophilic surfaces are easier to wet rather than de-wet, and vice versa for hydrophobic surfaces. This puts the classical view – that hydrophilic substrates are high-friction surfaces and hydrophobic substrates are low-friction surfaces – under doubt. Indeed, it is not clear whether the value of contact line friction can be predicted from $\theta _0$ alone (Liu, Yu & Wang 2019b; Wang 2019). Reasoning with the frame of mind of KMT instead, line friction asymmetry hints at some complex physics modulating adsorption/desorption of molecules at the contact line. Fluid/surface interface energy alone is not sufficient to describe asymmetry between adsorption and desorption. More in-depth examination of sub-continuum fluid displacement near contact lines will be required to arrive at a physical understanding of this process. When the asymmetry is understood, it is straightforward to impose different $\mu _f$ values in PF and different $\lambda$ values in VOF for advancing and receding contact lines.

7.3. Potential modelling directions

It has been proposed previously that a better way to model the moving contact line is to use the so-called generalized Navier boundary condition (GNBC), as first hinted by Blake (1993). This approach was later evaluated against MD simulations (Qian et al. 2003; Qian, Wang & Sheng 2006; Mohand et al. 2019), and good agreement has been found. However, recently Lācis et al. (2020) have tried to match GNBC with MD simulations exhibiting negligible slip but were not successful in demonstrating any advantage of GNBC against no-slip and Navier-slip boundary conditions.

In this work, we also use substrate with negligible hydrodynamic slip ($\ell _s = 0.44$ nm for $\theta _0 = 127^\circ$, and $\ell _s = 0.0$ nm for $\theta _0 = 38^\circ \text {--}95^\circ$). Previously (Lācis et al. 2020), for $\theta _0 = 95^\circ$ we extracted $\ell _s = 0.17$ nm, which is close to what we have in the current work. However, the most important overlap between our current and previous (Lācis et al. 2020) work is that the slip length imposed at the solid wall is constant over the surface. This is the main reason why the GNBC for the selected system (Lācis et al. 2020) did not exhibit any advantage. Very small slip length corresponds to very large friction at the contact line. Consequently, the addition of uncompensated Young's stress does not lead to significantly modified flow near the contact line.

Nevertheless, there are no solid arguments as to why the effective friction exactly at the contact line should be the same as below the bulk liquid. It could very well be that the effective slippage (friction) exactly at the contact line must be prescribed larger (smaller) compared to below the bulk liquid. This has the potential for improving both PF and VOF ability to match the MD results. An alternative approach to prescribing larger friction in VOF simulations would be to use so-called staggered slip or negative slip (Hartmann et al. 2021). Another effect, which was not considered in this work, is so-called disjoining pressure (Pismen & Pomeau 2000). Including it would also allow for direct modification of contact line motion. Furthermore, more detailed studies of local rheological effects (including the orientation parameter of water) could provide insight of detailed mechanisms governing film formation and contact line friction. An accurate description of Navier-slip related friction near the contact line could improve results attainable also using the GNBC condition. This is because the friction parameter is an important input in the GNBC. An infinite friction parameter renders GNBC ineffective, whereas gradually reducing the friction parameter (increasing slip) would amplify the GNBC effect on the velocity near the contact line.

Another aspect is that the MD results hint towards asymmetric properties of advancing and receding contact lines. Consequently, contact line friction and possibly local slip length could be different. This enlarges the parameter space enormously. Fundamental investigations into the possible cause of asymmetry between the adsorption and desorption process would be welcome. These could potentially shed more light on what input should be given to CFM models to match the molecular reality. Alternatively, some kind of hybrid methods that allow the matching between MD and Navier–Stokes solvers (Hadjiconstantinou 1999; Zhang, Borg & Reese 2017; Borg et al. 2018) could be used. These approaches would alleviate the need to understand the asymmetric properties of the contact line and provide direct coupling between MD and CFM solvers.

Finally, there are other CFM models available that could be benchmarked against the molecular data produced through this work. For example, there exist different PF models, such as van der Waals (Laurila et al. 2012) or Cahn–Allen (Eggleston, McFadden & Voorhees 2001). The level-set model (Tornberg & Engquist 2000) or lattice Boltzman method (Chen et al. 2014) are other potential candidates for the simulation of multiphase flows; for instance, Latva-Kokko & Rothman (2007) showed how a no-slip LBM model is able to capture automatically the speed-dependent dynamic contact angle and the interface-local slip length. The number of freely adjustable parameters differs between all alternatives. However, the issue of the appropriate velocity boundary condition will be shared between all of models based on single continuum description. A good recent classification of multiphase models can be found in work by Soligo, Roccon & Soldati (2021). It must be recognized that the hybrid models (Zhang et al. 2017; Liu et al. 2021) alleviate the need to understand the fundamental mechanisms near the moving contact line, and provide a way to couple MD and CFM directly. This is another alternative that should be considered for efficient simulations of multiphase systems. If the stick-slip like oscillations of the contact line are important, then some other means of continuum modelling can be pursued. For example, the fluctuating hydrodynamic interfaces model proposed by Flekkøy & Rothman (1996) and Smith et al. (2016) could be evaluated as a suitable choice. Finally, we mention the possibility of regularizing the contact line stress singularity via Brinkman's model for porous surfaces (Devauchelle, Josserand & Zaleski 2007).

Despite the missing physical mechanisms in PF and VOF, we have demonstrated that sufficiently accurate predictions of interface shape, drop displacement and critical capillary number can be obtained. The only prerequisite is that the PF and VOF simulations have to be calibrated with the MD once for given $\theta _0$.

8. Conclusions

We have calibrated a standard Cahn–Hilliard PF model, as well as a standard geometric volume-of-fluid model with Cox-like wetting condition, against MD simulations of water over a no-slip substrate. The no-slip behaviour in the MD system is an outcome of electrostatic bonds between polar water molecules and polar wall molecules. Two parameters (mobility $M$ and contact line friction $\mu _f$) were adjusted in PF simulations, and one parameter (Cox cut-off length scale $\lambda$) in VOF simulations. Four different equilibrium contact angles ($\theta _0 = 127^\circ$, $95^\circ$, $69^\circ$ and $38^\circ$) were investigated. For each $\theta _0$, the largest steady stick-slip free simulation was selected for calibration. The PF and VOF models were calibrated to match the steady $\Delta x$ – single scalar macroscopic measurement – observed in MD. Using the calibrated parameters, a series of simulations was carried out for other solid wall velocities. We demonstrated that CFM simulations can sufficiently accurately predict the drop displacement without any additional adjustments. The critical capillary number predictions from CFM models also displayed good agreement with MD data.

In addition, we have showcased predictions for two unsteady sheared droplet configurations of $\theta _0 = 95^\circ$ and $\theta _0 = 38^\circ$. The CFM models predicted all qualitative features of the MD simulations. For $\theta _0 = 95^\circ$, drop displacement and break-up in half were predicted. For $\theta _0 = 38^\circ$, thin film deposition and coalescence with the periodic image were predicted. Despite the quantitative differences in the time instances of the similar shapes, the CFM predictions exhibited good agreement with MD results.

We identified molecular physics that to the best of the authors’ knowledge have not been previously reported. We extracted line friction directly from MD and compared it with the PF calibration results. The MD results showed the same trend as obtained with PF. In addition, the resulting contact line friction $\mu _f$ values were asymmetric between advancing and receding contact lines for $\theta _0 = 127^\circ$ and $38^\circ$.

Finally, we have discussed the variations of PF and VOF parameters for matching the MD results for all $\theta _0$ values. We identified that the currently chosen CFM models seem to be lacking a way to describe reduced friction near the contact line for $\theta _0 \geqslant 95^\circ$. A possible future direction to remedy this shortcoming would be to introduce larger slippage (lower friction) locally near the contact line. In a continuum setting, this could correspond to having a spatially varying slip. In addition, we anticipate that the asymmetric behaviour of the advancing and receding contact lines is the source of the inaccuracies in PF and VOF when predicting the $\theta _0 = 38^\circ$ results.

We have demonstrated that by calibrating the CFM once for each $\theta _0$ by targeting a single global measure $\Delta x$, it is possible to obtain many accurate predictions of interface shape $\Delta x$ as a function of ${Ca}$, and also prediction of ${Ca}_c$. This is despite the fact that the selected MD configuration exhibits practically negligible slippage. The selected CFM models seem to perform very well for close to neutral and slightly hydrophobic configurations ($\theta _0 = 95^\circ$ and $69^\circ$). On the other hand, more deviations were observed for hydrophobic ($\theta _0 = 127^\circ$) and hydrophilic ($\theta _0 = 38^\circ$) configurations. These accuracy limits must be kept in mind if these CFM models are applied in similar conditions.

The results of this study continue to enrich our understanding of connections between continuum mechanics simulations and molecular reality. We believe that this work provides important insights into PF and VOF models, associated open questions, and the required calibration procedures. Properly calibrated, both PF and VOF can serve as useful tools for investigations of technological applications.