1. Introduction
Wetting is defined as the collective set of phenomena that occur when a solid is exposed to a fluid and it comprises a number of interfacial effects. These include the interplay of capillary forces, which are net forces on the fluid, especially in the three-phase contact region, resulting from solid–fluid and fluid–fluid interactions as well as external (e.g. gravitational) forces. In the dynamic case, viscous forces additionally influence the wetting behaviour. Wetting is important in several industrial areas (de Gennes, Brochard-Wyart & Quere Reference de Gennes, Brochard-Wyart and Quere2004; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009), such as in the chemical industry for the spreading of paints (Schoff Reference Schoff1992; Prajapati & Arya Reference Prajapati and Arya2024), in soil science for the study of the penetration of liquids into porous rocks (Melrose Reference Melrose1965; Garfi et al. Reference Garfi, John, Rücker, Lin, Spurin, Berg and Krevor2022) or in the construction industry for waterproofing of concrete (McCarter Reference McCarter1998; Muhammad et al. Reference Muhammad, Keyvanfar, Abd. Majid, Shafaghat and Mirza2015). It is also ubiquitous in life sciences, e.g. the rise of sap in plants, the adhesion of parasites on wet surfaces or even the wetting of the eye aided by special proteins (Holly & Lemp Reference Holly and Lemp1971; de Gennes et al. Reference de Gennes, Brochard-Wyart and Quere2004). The static (equilibrium) contact angle
$\theta$
of sessile droplets is an important measure of the static wetting behaviour and has been the subject of several theoretical (Graham-Eagle & Pennell Reference Graham-Eagle and Pennell2000; Vafaei & Podowski Reference Vafaei and Podowski2005; Pereira & Kalliadasis Reference Pereira and Kalliadasis2012; Becker et al. Reference Becker, Urbassek, Horsch and Hasse2014; Jasper & Anand Reference Jasper and Anand2019; Yatsyshin & Kalliadasis Reference Yatsyshin and Kalliadasis2021) and experimental (Vafaei et al. Reference Vafaei, Borca-Tasciuc, Podowski, Purkayastha, Ramanath and Ajayan2006; Ghasemi & Ward Reference Ghasemi and Ward2010; Sharp, Farmer & Kelly Reference Sharp, Farmer and Kelly2011; Schuster, Schvezov & Rosenberger Reference Schuster, Schvezov and Rosenberger2015) studies. In the case of total wetting,
$\theta ={0}{^{\circ }}$
, and for total dewetting,
$\theta = {180}{^{\circ }}$
, with partial wetting between these limits.
Analogously to the equilibrium case, the dynamic wetting behaviour, i.e. when a droplet placed on a solid spreads or is driven by an external force (e.g. gravity), can be characterised by dynamic contact angles (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009). Importantly, depending on the direction of the movement of the three-phase contact region, an advancing or receding dynamic contact angle,
$\Theta _{a}$
or
$\Theta _{r}$
, which generally differ, is observed as visualised in figure 1(b). Dynamic contact angles depend on velocity or, equivalently, on the driving force causing the movement of the contact region (
$f_x$
in figure 1) (Voinov Reference Voinov1976; Cox Reference Cox1986). The difference between advancing and receding contact angles is commonly attributed to chemical or spatial inhomogeneities (roughness) of the solid surface and sometimes termed dynamic contact angle hysteresis in the literature (Butt et al. Reference Butt2022).
The dynamics close to the three-phase contact region (in macroscopic studies often called contact line) is at the centre of research in many works (Koplik, Banavar & Willemsen Reference Koplik, Banavar and Willemsen1988; Decker & Garoff Reference Decker and Garoff1997; Jacqmin Reference Jacqmin2000; Narhe, Beysens & Nikolayev Reference Narhe, Beysens and Nikolayev2004; Lee & Müller-Plathe Reference Lee and Müller-Plathe2022). This can be attributed to the inherent challenges of modelling the time evolution of contact lines, where continuum fluid dynamics fails due to a stress singularity. Particularly, the widely used no-slip boundary condition requires the contact line to be at rest, which is obviously not correct for spreading droplets (Huh & Scriven Reference Huh and Scriven1971). Since the dynamic behaviour of the contact line is strongly affected by molecular interactions, it needs to be studied at the microscopic scale of individual molecules (i.e. a few Å to nm).
Snapshot of a droplet moving parallel to the solid–fluid interface due to an external force
$f_x$
in a system with dimensions
$L_x$
and
$L_y$
. (a) Atomistic model used for equilibrium and non-equilibrium molecular dynamic (NEMD) simulations with individual solid (red) and fluid (green) particles. (b) Density profile from hydrodynamic DFT with molecular layering, surface roughness as well as advancing and receding dynamic contact angles,
$\Theta _{a}$
and
$\Theta _{r}$
, respectively.
Several routes have been taken to address this challenge. Microscopic information can be used as an input for macroscopic models. This procedure is followed, for example, by sharp interface models, where interfaces are taken to be infinitely thin. These models describe a slip length based on microscopic considerations, which is then included in continuum fluid dynamics models to avoid the stress singularity (Hocking Reference Hocking1983, Reference Hocking1992; Eggers Reference Eggers2004). A second route is through (atomistic) molecular simulations, where equilibrium molecular dynamics and non-equilibrium molecular dynamics (NEMD) simulations provide insights into a variety of microscopic phenomena related to wetting (Koplik et al. Reference Koplik, Banavar and Willemsen1988; Hong, Ha & Balachandar Reference Hong, Ha and Balachandar2009; Yuan & Zhao Reference Yuan and Zhao2013; Li et al. Reference Li, Yan, Fichthorn and Yu2018; Li & Zhang Reference Li and Zhang2019; Lee & Müller-Plathe Reference Lee and Müller-Plathe2022). MD provides detailed results, but is limited to small length and time scales. A further challenge for these models is that MD provides quantities with statistical uncertainties, whereas macroscopic (continuum) models are deterministic. This effectively limits models to using averaged or fitted parameters, which are determined from MD, in macroscopic models instead of density and velocity fields (Qian, Wang & Sheng Reference Qian, Wang and Sheng2006; Nakamura et al. Reference Nakamura, Carlson, Amberg and Shiomi2013; Diewald et al. Reference Diewald, Lautenschlaeger, Stephan, Langenbach, Kuhn, Seckler, Bungartz, Hasse and Müller2020). A third approach are diffuse interface models, which consider the interfaces to be of finite thickness. This physical thickness is typically of the order of a few molecular diameters (Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009), which renders numerical resolution of macroscopic problems computationally expensive (Shokrpour Roudbari & van Brummelen Reference Shokrpour Roudbari and van Brummelen2019). In many studies, however, these diffuse interface models make no attempt to resemble the physical interface; rather, the diffuse interface is a numerical mean to distinguish and track the interface between fluid phases. These models avoid the stress singularity at the contact line and provide evolution equations for density or composition and sometimes also for velocity fields. Successful examples are models based on the coupled Cahn–Hilliard/Navier–Stokes (Cahn & Hilliard Reference Cahn and Hilliard1958; Novick-Cohen Reference Novick-Cohen2008; Demont, Stoter & van Brummelen Reference Demont, Stoter and van Brummelen2023), Navier–Stokes–Korteweg equations (Dunn & Serrin Reference Dunn and Serrin1985; Heida & Málek Reference Heida and Málek2010; Diehl et al. Reference Diehl, Kremser, Kröner and Rohde2016) or dynamic density functional theory (DDFT) (Marconi & Tarazona Reference Marconi and Tarazona1999, Reference Marconi and Tarazona2000).
The Cahn–Hilliard equation was initially derived for the study of phase separation in a two-component system, where a phase-field variable or order parameter is employed to differentiate between the two immiscible phases. The coupled Cahn–Hilliard/Navier–Stokes equations include the gradient of the order parameter in the momentum balance to describe the time evolution of the velocity and the order parameter for incompressible fluids. While the Cahn–Hilliard/Navier–Stokes equations have been applied to the study of contact region dynamics in incompressible systems (Jacqmin Reference Jacqmin2000; Qian et al. Reference Qian, Wang and Sheng2006; Yue, Zhou & Feng Reference Yue, Zhou and Feng2010; Yue & Feng Reference Yue and Feng2011), the approach is not applicable to compressible fluids.
For compressible two-phase systems, the Navier–Stokes–Korteweg equations can be employed (Korteweg Reference Korteweg1901; Dunn & Serrin Reference Dunn and Serrin1985; Heida & Málek Reference Heida and Málek2010; Diehl et al. Reference Diehl, Kremser, Kröner and Rohde2016). The so-called Korteweg tensor, a constitutive relation which includes density gradients and accounts for capillarity effects, is incorporated into the momentum balance of the Navier–Stokes equations in addition to viscous stresses. Microscopic investigations of dynamic contact angles using the Navier–Stokes–Korteweg model yield accurate results when compared with MD simulations (Diewald et al. Reference Diewald, Heier, Lautenschläger, Horsch, Kuhn, Langenbach, Hasse and Müller2019, Reference Diewald, Lautenschlaeger, Stephan, Langenbach, Kuhn, Seckler, Bungartz, Hasse and Müller2020). The contact angle in the Navier–Stokes–Korteweg model, however, is prescribed by an input parameter chosen such that the contact angle agrees with the MD simulations. Additionally, to obtain the correct vapour–liquid surface tension, a second parameter,
$\kappa$
, needs to be determined, e.g. from MD simulations. The Helmholtz energy functional, which is used in Cahn–Hilliard/Navier–Stokes and Navier–Stokes–Korteweg to model phase behaviour and capillary forces, only contains information of the immediate surroundings (by means of density gradients), and molecular layering at the solid–fluid interface can therefore not be resolved (Sauer & Gross Reference Sauer and Gross2017; Diewald et al. Reference Diewald, Lautenschlaeger, Stephan, Langenbach, Kuhn, Seckler, Bungartz, Hasse and Müller2020).
Dynamic density functional theory (DDFT) is an alternative diffuse interface model. It can be viewed as a dynamic extension of classical density functional theory (DFT), which in general employs non-local Helmholtz energy functionals, to non-equilibrium situations (Michael te Vrugt & Wittkowski Reference Michael te Vrugt and Wittkowski2020). It provides time evolution equations for the density and, in some extensions, also for the momentum density. It was initially derived by Marconi & Tarazona (Reference Marconi and Tarazona1999, Reference Marconi and Tarazona2000) based on the Langevin equation. While stochastic and deterministic versions of DDFT can be obtained (Archer & Rauscher Reference Archer and Rauscher2004), this paper focuses on the deterministic version. Alternative routes for the derivation of DDFT are based on integrating the Smoluchowski equation (see e.g. Archer & Evans Reference Archer and Evans2004) or using the projection operator formalism (te Vrugt & Wittkowski Reference te Vrugt and Wittkowski2019), where the phase-space probability density is projected onto the density (see e.g. Yoshimori Reference Yoshimori2005; Español & Löwen Reference Español and Löwen2009). A theoretical existence proof based on the Liouville equation was provided by Chan & Finken (Reference Chan and Finken2005). Archer (Reference Archer2009) derived the first DDFT for atomic/molecular fluids including a momentum balance, whereas Goddard, Nold & Kalliadasis (Reference Goddard, Nold and Kalliadasis2013) extended DDFT towards mixtures of colloids. DDFT was applied to spinodal decomposition (Archer & Evans Reference Archer and Evans2004), to phase separation of colloids confined in a cavity (Archer Reference Archer2005), to determine the van Hove correlation function (Hopkins et al. Reference Hopkins, Fortini, Archer and Schmidt2010; Stopper et al. Reference Stopper, Marolt, Roth and Hansen-Goos2015a ,Reference Stopper, Roth and Hansen-Goos b , Reference Stopper, Roth and Hansen-Goos2016) and to study the phase transition of colloidal systems under shear influence (Stopper & Roth Reference Stopper and Roth2018). In addition, DDFTs including inertia (Archer Reference Archer2006), hydrodynamic interactions between colloids (Goddard et al. Reference Goddard, Nold, Savva, Pavliotis and Kalliadasis2012a ,Reference Goddard, Nold, Savva, Yatsyshin and Kalliadasisb; Durán-Olivencia et al. Reference Durán-Olivencia, Goddard and Kalliadasis2016; Goddard, Nold & Kalliadasis Reference Goddard, Nold and Kalliadasis2016) or fluctuations (Archer & Evans Reference Archer and Evans2004; Donev & Vanden-Eijnden Reference Donev and Vanden-Eijnden2014; Durán-Olivencia et al. Reference Durán-Olivencia, Yatsyshin, Goddard and Kalliadasis2017; Russo et al. Reference Russo, Durán-Olivencia, Yatsyshin and Kalliadasis2020, Reference Russo, Perez, Durán-Olivencia, Yatsyshin, Carrillo and Kalliadasis2021) were developed.
In contrast to equilibrium DFT, DDFT is not an exact theory, since the so-called adiabatic approximation is employed. The non-equilibrium two-body spatial correlation function is assumed to be the same as in equilibrium (for the same density profile) (Marconi & Tarazona Reference Marconi and Tarazona1999, Reference Marconi and Tarazona2000; Archer Reference Archer2009; Michael te Vrugt & Wittkowski Reference Michael te Vrugt and Wittkowski2020). Consequently, the equilibrium Helmholtz energy functional can be applied to the non-equilibrium system (Stierle & Gross Reference Stierle and Gross2021). A drawback of this approximation is that history-dependent superadiabatic forces (e.g. memory effects) can not be described by DDFT. An exact version of DDFT which includes superadiabatic forces is given by power functional theory (Schmidt & Brader Reference Schmidt and Brader2013; Brader & Schmidt Reference Brader and Schmidt2015; de las Heras & Schmidt Reference de las Heras and Schmidt2018; Schmidt Reference Schmidt2018).
Stierle & Gross (Reference Stierle and Gross2021) proposed hydrodynamic DFT based on a variational principle (Herivel Reference Herivel1955; Serrin Reference Serrin1959). It describes the time evolution of density and momentum including inertial effects, and viscous dissipation for pure substances and mixtures. The Cauchy pressure tensor is modelled assuming a Newtonian fluid with a spatially varying viscosity coefficient. A DFT term with the functional derivative of the Helmholtz energy functional is included in the momentum balance and non-local Helmholtz energy functionals can be employed. In a recent work, Nold et al. (Reference Nold, Goddard, Sibley and Kalliadasis2024) employ a similar approach to study the fluid close to the contact line. The no-slip boundary condition is assumed between the solid and the fluid. The evolution from a non-equilibrium contact angle to equilibrium is investigated, thereby analysing the influence of compression and shear on the slip length and contact line friction. The DFT term, which appears in hydrodynamic DFT, captures molecular forces, such as capillary forces, in a predictive manner. The predictive power is not limited to vapour–liquid interfaces, but also captures adsorption and wetting effects. While not restricted to a certain Helmholtz energy functional, it was combined with a Helmholtz energy functional based on the PC-SAFT (Gross & Sadowski Reference Gross and Sadowski2000, Reference Gross and Sadowski2001, Reference Gross and Sadowski2002; Gross et al. Reference Gross, Spuhl, Tumakaka and Sadowski2003) equation of state, which shows good agreement with experimental data in a wide range of equilibrium applications (Gross Reference Gross2009; Rehner & Gross Reference Rehner and Gross2018; Sauer & Gross Reference Sauer and Gross2019; Rehner, Bursik & Gross Reference Rehner, Bursik and Gross2021; Nitzke et al. Reference Nitzke, Stierle, Stephan, Pfitzner, Gross and Vrabec2023; Bursik, Eller & Gross Reference Bursik, Eller and Gross2024a ; Stierle et al. Reference Stierle, Bauer, Thiele, Bursik, Rehner and Gross2024). Specifically, density profiles in confined systems (Sauer & Gross Reference Sauer and Gross2017; Bursik et al. Reference Bursik, Eller and Gross2024b ), adsorption isotherms (Sauer & Gross Reference Sauer and Gross2019; Stierle et al. Reference Stierle, Bauer, Thiele, Bursik, Rehner and Gross2024) and contact angles of sessile droplets including molecular layering effects (Sauer et al. Reference Sauer, Terzis, Theiss, Weigand and Gross2018) are accurately predicted. Other Helmholtz energy functionals, such as that describing density gradient theory and leading to the Navier–Stokes–Korteweg equation (with a less accurate description of solid–fluid interfaces), can also be employed. An important feature of hydrodynamic DFT is that in a bulk phase, i.e. far enough away from interfaces, the model is equal to the isothermal (continuum) Navier–Stokes equation (Stierle & Gross Reference Stierle and Gross2021).
For these reasons, we propose hydrodynamic DFT as a suitable candidate for a unified model applicable to the description of the dynamic behaviour of fluids – and especially wetting phenomena – from the microscopic to the macroscopic scales. However, to date, hydrodynamic DFT has only been applied in a qualitative study of coalescence phenomena in one dimension (Stierle & Gross Reference Stierle and Gross2021). The contribution of this paper is, thus, to conduct a (i) higher-dimensional, (ii) quantitative, (iii) wetting-related investigation of the predictive capabilities of hydrodynamic DFT to assess its potential for modelling dynamic wetting behaviour. In particular, we apply hydrodynamic DFT to predict dynamic contact angles of two-dimensional microscopic sessile droplets driven by an external force. We use the already mentioned Helmholtz energy functional based on the PC-SAFT equation of state to model methane as an exemplary fluid. We obtain local values for the shear viscosity from a recently developed generalised entropy scaling model (Bursik et al. Reference Bursik, Eller and Gross2024b ). We analyse the importance of this local transport coefficient model by contrasting the results obtained from a model for shear viscosities in bulk phases, as it is typically employed in continuum approaches. All results are assessed by comparison to NEMD simulations for a Lennard–Jones fluid using PC-SAFT parameters for methane. Because this study focuses on the dynamic behaviour of droplets, with emphasis on the viscosity model, we choose molecular models for which reasonable agreement for (static) equilibrium properties can be expected. The PC-SAFT model has been shown to provide rather satisfactory predictions of the equilibrium properties of the Lennard–Jones fluid (Sauer & Gross Reference Sauer and Gross2017). We consider a Lennard–Jones fluid using methane parameters to generate illustrative results that are more readily interpretable.
This paper is structured as follows. We present the balance equations of hydrodynamic DFT, the Helmholtz energy functional based on PC-SAFT and the shear viscosity model in § 2. Numerical details on the discretisation, a description of NEMD simulations and the procedure for determining contact angles from density profiles are given in § 3. In § 4, we present and discuss results for density and velocity profiles including advancing and receding contact angles from hydrodynamic DFT, using the generalised entropy scaling model and a viscosity model for bulk phases, and compare them with results from NEMD simulations.
2. Hydrodynamic density functional theory
First, we introduce classical (equilibrium) DFT, followed by a short presentation of the hydrodynamic DFT model. The local shear viscosity is modelled using generalised entropy scaling (Bursik et al. Reference Bursik, Eller and Gross2024b ), which is summarised subsequently.
2.1. Classical (equilibrium) density functional theory
Classical DFT is an exact theory which describes inhomogeneous fluids in equilibrium using a grand canonical density functional
$\Omega$
at constant chemical potential
$\mu$
, volume
$V$
and temperature
$T$
according to
\begin{equation} \Omega \left [\rho (\boldsymbol{r})\right ]=F\left [\rho (\boldsymbol{r})\right ]-\int \rho (\boldsymbol{r}) \left (\mu -V^{{\textit{ext}}}(\boldsymbol{r})\right )\, \mathrm{d}\boldsymbol{r} \end{equation}
where
$F$
is the Helmholtz energy functional and
$\rho$
is the density profile which depends on the spatial coordinate
$\boldsymbol{r}$
. The square brackets denote that
$F$
and
$\Omega$
are functionals of the (number) density
$\rho$
. In this work, we split the external potential
$V^{{\textit{ext}}}$
into two contributions according to
\begin{equation} V^{{\textit{ext}}}= V^{{\textit{ext},g}} + V^{{\textit{ext},\textit{sf}}}, \end{equation}
where
$V^{{\textit{ext},g}}$
captures the effect of external potentials, such as gravity, on the fluid. The quantity
$V^{{\textit{ext},\textit{sf}}}$
captures the interactions of solid interaction sites with the fluid and therewith the effect of confinement. From the perspective of a fluid, the solid acts as a (solid–fluid) external potential
$V^{{\textit{ext},\textit{sf}}}$
.
At equilibrium, the grand canonical functional from (2.1) has a minimum and becomes the grand potential. Mathematically, this leads to the Euler–Lagrange equation
\begin{equation} \frac {\delta F[\rho ]}{\delta \rho (\boldsymbol{r})}-\mu +V^{{\textit{ext}}}(\boldsymbol{r})=0, \qquad \qquad \end{equation}
which includes the functional derivative
${\delta F[\rho ]}/{\delta \rho (\boldsymbol{r})}$
and allows to calculate the equilibrium density profile. In the context of (2.3), we use
$\rho (\boldsymbol{r})$
to denote the equilibrium density, while in latter chapters,
$\rho (\boldsymbol{r})$
is used for non-equilibrium density fields.
In this work, the Helmholtz energy functional is based on the PC-SAFT model, i.e. for homogeneous phases, the PC-SAFT equation of state is obtained. In PC-SAFT, fluid molecules are assumed to be chains of
$m_{{f}}$
tangentially bonded spheres, which are called segments. Here,
$m_{{f}}=1$
is used for the molecules in this work, but an extension to non-spherical molecules is readily available (Sauer & Gross Reference Sauer and Gross2017). For this choice of parameters, the PC-SAFT model reduces to a model for spherical molecules, which provides accurate results for the Lennard–Jones fluid. Each segment has a diameter parameter
$\sigma _{{\textit{ff}}}$
and an energy interaction parameter
$\varepsilon _{{\textit{ff}}}$
. The PC-SAFT parameters, which were previously fitted to vapour–liquid equilibrium data (Gross & Sadowski Reference Gross and Sadowski2001), are taken from the literature and given in § 3. The Helmholtz energy functional, similarly to the PC-SAFT equation of state, comprises additive contributions. For the molecules studied here (i.e. with
$m_{{f}}=1$
), the ideal gas, hard-sphere and dispersion contributions are employed as
$$F\left [\rho (\boldsymbol{r})\right ]=F_{{ig}}\left [\rho (\boldsymbol{r})\right ]+F_{{\textit{hs}}}\left [\rho (\boldsymbol{r})\right ] +F_{{disp}}\left [\rho (\boldsymbol{r})\right ] $$
$$ \equiv F_{{ig}}\left [\rho (\boldsymbol{r})\right ] + F_{{\textit{res}}}\left [\rho (\boldsymbol{r})\right ], $$
defining the residual Helmholtz energy
$F_{{\textit{res}}}[\rho (\boldsymbol{r})]$
.
The ideal gas contribution can be derived exactly from statistical mechanics. For the hard-sphere contribution, which models each segment as an impenetrable sphere, fundamental measure theory (Roth et al. Reference Roth, Evans, Lang and Kahl2002; Yu & Wu Reference Yu and Wu2002b ) is employed. The dispersion contribution describes the attractive, non-polar interactions of chain molecules modelled here by a functional developed by Sauer & Gross (Reference Sauer and Gross2017). This model can be extended to chain molecules (Tripathi & Chapman Reference Tripathi and Chapman2005a ,Reference Tripathi and Chapman b ), associating molecules (Yu & Wu Reference Yu and Wu2002a ; Sauer & Gross Reference Sauer and Gross2017) and polar molecules (Gross Reference Gross2005; Gross & Vrabec Reference Gross and Vrabec2006; Vrabec & Gross Reference Vrabec and Gross2008; Sauer & Gross Reference Sauer and Gross2017).
The interactions between solid and fluid are accounted for in the (solid–fluid) external potential
$V^{{\textit{ext},\textit{sf}}}$
. If the solid consists of atomistic Lennard–Jones interaction sites, the external potential due to dispersive interactions between the solid and the fluid can be calculated according to
\begin{equation} V^{{\textit{ext},\textit{sf}}}(\boldsymbol{r})=\sum _{\zeta =1}^{N_{{s}}} 4\varepsilon _{\zeta {f}} \left (\left (\frac {\sigma _{\zeta {f}}}{|\boldsymbol{r}_\zeta -\boldsymbol{r}|}\right )^{12}-\left (\frac {\sigma _{\zeta {f}}}{|\boldsymbol{r}_\zeta -\boldsymbol{r}|}\right )^6\right ) \end{equation}
with the index of the Lennard–Jones interaction sites
$\zeta$
and their position
$\boldsymbol{r}_\zeta$
as well as the total number of solid interaction sites
$N_{{s}}$
. In this work, the solid and the fluid consist of only one species respectively, such that the solid–fluid interaction parameters are equal for each solid–fluid pair (
$\sigma _{\zeta {f}}=\sigma _{{\textit{sf}}}$
and
$\varepsilon _{\zeta {f}}=\varepsilon _{{\textit{sf}}}$
). Rather than defining
$\sigma _{{\textit{sf}}}$
and
$\varepsilon _{{\textit{sf}}}$
, we prefer to specify solid interaction site parameters
$\sigma _{{ss}}$
and
$\varepsilon _{ss}$
assuming the Berthelot–Lorentz combining rules as
\begin{equation} \begin{aligned} \sigma _{{\textit{sf}}} & = \frac {\left (\sigma _{{ss}} + \sigma _{{\textit{ff}}}\right )}{2}, \\ \varepsilon _{{\textit{sf}}} & = \sqrt {\varepsilon _{{ss}} \varepsilon _{{\textit{ff}}}}. \end{aligned} \end{equation}
Knowing the positions and parameters of the solid interaction sites, the solid–fluid external potential
$V^{{\textit{ext},\textit{sf}}}$
, as given in (2.5), can be evaluated and used in equilibrium and hydrodynamic DFT. To obtain a two-dimensional external potential
$V^{{\textit{ext},\textit{sf}}}(x,y)$
, the three-dimensional external potential is first calculated for the depth of a single unit cell of the solid in
$z$
-direction (cf. figure 1) and then free-energy averaged in the same direction (Eller & Gross Reference Eller and Gross2021). The resulting external potential varies strongly perpendicular to the solid–fluid interface due to decaying solid–fluid interactions. Minor variations occur also in the direction parallel to the solid–fluid interface, which resemble the molecular roughness of the solid. More pronounced roughnesses are investigated by removing solid interaction sites at the interface to create stronger variations of the external potential parallel to the solid. Other external fields, such as gravity, are included by means of the gravitational external potential
$V^{{\textit{ext},g}}$
. In this work, an external force is applied in hydrodynamic DFT simulations to induce the movement of liquid droplets. An external force is equivalent to a negative gradient in the (gravitational) external potential.
DFT calculations are carried out using FeO
$_{{s}}$
(Rehner & Bauer Reference Rehner and Bauer2021; Rehner, Bauer & Gross Reference Rehner, Bauer and Gross2023), a framework for equations of state and DFT. The system is discretised in a two-dimensional Cartesian grid with 256 grid points and periodic boundary conditions in each direction. Initial density profiles for sessile droplets can be constructed from equilibrated droplets in a surrounding vapour phase or taken from MD results; both approaches provide the same final density profile within numerical accuracy. The same number of molecules as in the MD simulations projected to the two-dimensional system (i.e. the same average density) is chosen.
Since in the grand-canonical (
$\mu ,V,T$
) ensemble, different (correct) solutions can be obtained from the Euler–Lagrange equation (2.3), the number of molecules in the system is kept constant using a mathematical reformulation proposed by Rehner & Gross (Reference Rehner and Gross2018). Essentially, a minimisation with constraints is performed using Lagrange multipliers, during which the chemical potential is adjusted to achieve the desired ensemble-averaged number of molecules in the system (see Appendix A). Picard iterations and an Anderson mixing scheme (Anderson Reference Anderson1965, Reference Anderson2019) are used to obtain the equilibrium density profile from (2.3). The convolutions which appear in the Helmholtz energy functionals and in (2.11), are implemented using fast Fourier transforms (Stierle et al. Reference Stierle, Sauer, Eller, Theiss, Rehner, Ackermann and Gross2020; de Morais Sermoud et al. Reference de Morais Sermoud, de Freitas Gonçalves, Barreto, Franco, Tavares and Castier2024).
2.2. Model equations
We present the hydrodynamic DFT model following Stierle & Gross (Reference Stierle and Gross2021), originally developed for pure fluids by Archer (Reference Archer2009). Even though it is applicable to mixtures, we solely provide equations for pure substances to improve clarity. The equations comprise a mass balance (continuity equation) and momentum balance according to (Stierle & Gross Reference Stierle and Gross2021)
$$ \qquad\qquad\qquad \frac {\partial \left (M\rho \right )}{\partial t} + {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! (M\rho \boldsymbol{v}) = 0, $$
$$ \frac {\partial (M\rho \boldsymbol{v})}{\partial t} + {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! \left (M\rho \boldsymbol{v} \boldsymbol{v}^\intercal \right ) = - \rho {\boldsymbol\nabla} \left ( \frac {\delta F}{\delta \rho } + V^{{\textit{ext}}} \right ) - {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! \boldsymbol{\tau }, $$
where
$M$
is the molecular mass and
$\boldsymbol{v}$
is the velocity vector. The Helmholtz energy functional
$F$
captures intermolecular fluid–fluid interactions and is described in § 2.1. The DFT term
$\rho {\boldsymbol\nabla} ( {\delta F}/{\delta \rho } )$
includes the functional derivative
${\delta F}/{\delta \rho }$
and captures the influence of fluid–fluid interfaces on the momentum, whereas solid–fluid interactions, e.g. in a confined system, and gravitational external potentials are included using the external potential
$ V^{{\textit{ext}}}$
according to (2.2). In bulk phases, i.e. sufficiently far from solid–fluid or fluid–fluid interfaces, the solid–fluid external potential
$\rho {\boldsymbol\nabla} V^{{\textit{ext},\textit{sf}}}$
is zero, the gravitational external potential term
$\rho {\boldsymbol\nabla} V^{{\textit{ext},g}}$
becomes an external force
$\rho f^{{\textit{ext},g}}$
and the DFT term
$\rho {\boldsymbol\nabla} ( {\delta F}/{\delta \rho } )$
simplifies to the pressure gradient
${\boldsymbol\nabla} p$
using the Gibbs–Duhem equation. This results in the momentum balance for homogeneous fluids known from the isothermal Navier–Stokes equations.
The viscous stresses in (2.7b ) are modelled assuming a Newtonian fluid according to
\begin{equation} \boldsymbol{\tau } = -\chi \left ( {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! \boldsymbol{v} \right ) \mathbb{I} - \eta \left ( {\boldsymbol\nabla} \boldsymbol{v} + ({\boldsymbol\nabla} \boldsymbol{v})^\intercal - \frac {2}{3} \left ( {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! \boldsymbol{v} \right ) \mathbb{I} \right ), \end{equation}
where
$\chi$
is the volume viscosity and
$\eta$
is the shear viscosity. The first term in (2.8) describes the dilatation, while the second, traceless term describes the viscous shear contribution (Stierle & Gross Reference Stierle and Gross2021). In this work, the volume viscosity
$\chi$
is neglected, as often done for liquid systems, and more detail on the shear viscosity
$\eta$
is given in § 2.3.
2.3. Modelling the shear viscosity
For hydrodynamic DFT, i.e. in inhomogeneous systems, local values for the shear viscosity can be determined by applying entropy scaling locally following our previous work (Bursik et al. Reference Bursik, Eller and Gross2024b
), which was based on the homogeneous approach of Lötgering-Lin et al. (Reference Lötgering-Lin, Fischer, Hopp and Gross2018). Instead of a residual entropy as in the homogeneous case, local values for the residual entropy density
$\tilde {s}_{{\textit{res}}}(\boldsymbol{r})$
, an entropy per unit volume, can be calculated. From the residual Helmholtz energy density
$f_{{\textit{res}}}(\boldsymbol{r})$
, which is related to the Helmholtz energy by
$F_{{\textit{res}}} = \int f_{{\textit{res}}}(\boldsymbol{r}) \,\mathrm{d} \boldsymbol{\it r}$
, the residual entropy density is derived according to
\begin{equation} \tilde {s}_{{\textit{res}}}(\boldsymbol{r}) = -\left (\frac {\partial f_{{\textit{res}}}(\boldsymbol{r})}{\partial T}\right )_{\rho ,V}. \end{equation}
We define a dimensionless residual entropy profile per molecule
$s_{{\textit{res}}}^\#(\boldsymbol{r})$
as
\begin{equation} s_{{\textit{res}}}^{\#}(\boldsymbol{r}) = \frac {\tilde {s}_{{\textit{res}}}(\boldsymbol{r})}{\bar {\rho }^{{ES}}(\boldsymbol{r})k_{{B}}}, \end{equation}
where a weighted density, an average density that takes non-local effects into account, is used as
\begin{equation} \bar {\rho }^{{ES}}(\boldsymbol{r}) = \frac {3}{4\pi \left (\psi \varepsilon _{{\textit{sf}}}d\right )^3} \int \rho (\boldsymbol{r}^{\prime}) \Theta (\psi \varepsilon _{{\textit{sf}}} d - |\boldsymbol{r}-\boldsymbol{r}'|) \,\mathrm{d} \boldsymbol{r}^{\prime} ,\end{equation}
with the temperature-dependent, effective hard-sphere diameter (Gross & Sadowski Reference Gross and Sadowski2001)
$d(T) = \sigma _{{\textit{ff}}} \bigl (1-0.12\exp (-3{\varepsilon _{{\textit{ff}}}}/{k_{B}T}) \bigr )$
and the Heaviside step function
$\Theta$
. An adjustable parameter
$\psi$
is required to quantitatively capture the influence of solid–fluid interactions on the viscosity in the close vicinity of the interface by varying the convolution radius
$\psi \varepsilon _{{\textit{sf}}} d$
. It was shown that once adjusted for a given system (by a single fluid-phase MD simulation), the parameter
$\psi$
is transferable to different temperatures, densities, shear rates and solid–fluid interactions described by
$\varepsilon _{{\textit{sf}}}$
, which influences the wetting behaviour (Bursik et al. Reference Bursik, Eller and Gross2024b
). In this work,
$\psi =1$
is adjusted to the velocity profile of liquid-phase flow in the same geometry from a single steady-state NEMD simulation and no wetting information enters the parameter adjustment (see Appendix C).
A third-order polynomial ansatz function (Lötgering-Lin et al. Reference Lötgering-Lin, Fischer, Hopp and Gross2018; Bursik et al. Reference Bursik, Eller and Gross2024b ) is then evaluated locally
\begin{equation} \ln \left (\frac {\eta ^{{ES}}(\boldsymbol{r})}{\eta _{{ref}}(\boldsymbol{r})}\right )=A + B s_{{\textit{res}}}^\#(\boldsymbol{r}) + C \left (s_{{\textit{res}}}^\#(\boldsymbol{r})\right )^2 + D \left (s_{{\textit{res}}}^\#(\boldsymbol{r})\right )^3, \end{equation}
where
$A$
,
$B$
,
$C$
and
$D$
are substance specific parameters that were previously fitted to experimental data of pure substances in homogeneous phases by Lötgering-Lin et al. (Reference Lötgering-Lin, Fischer, Hopp and Gross2018). The reference viscosity is defined by
\begin{equation} \eta _{{ref}}(\boldsymbol{r}) = \eta _{\textit{CE}} + \bar {\rho }_{{s}}(\boldsymbol{r})\sigma _{{ss}}^3\eta _{{s},\infty }. \end{equation}
The Chapman–Enskog viscosity
$\eta _{\textit{CE}}$
(Hirschefelder, Curtiss & Bird Reference Hirschefelder, Curtiss and Bird1954) is used as the reference for calculating the dimensionless viscosity in homogeneous systems and inhomogeneous systems without solid–fluid interfaces. It is calculated as
\begin{equation} \eta _{\textit{CE}} = \frac {5}{16}\frac {\sqrt {M k_{{B}}T/(N_{{A}}\pi )}}{\sigma ^2 \Omega ^{(2,2)\#}}, \end{equation}
where an empirical correlation (Neufeld, Janzen & Aziz Reference Neufeld, Janzen and Aziz2003) is used for the dimensionless collision integral
$\Omega ^{(2,2)\#}$
. In the vicinity and within the solid, the reference viscosity needs to be adjusted, since the Chapman–Enskog viscosity is valid only for dilute vapour phases. The second contribution in (2.13) accounts for the low mobility of fluid molecules within the solid and the influence of the solid–fluid interface on viscosity (Bursik et al. Reference Bursik, Eller and Gross2024b
). The weighted density of solid interaction sites
$\bar {\rho }_{{s}}$
is determined from
\begin{equation} \bar {\rho }_{{s}} =\int \rho _{{s}} (\boldsymbol{r}' ) \Theta (R-|\boldsymbol{r}-\boldsymbol{r}'| ) \,\mathrm{d} \boldsymbol{r}', \end{equation}
where
$\rho _{{s}}$
is the density of solid interaction sites and
$\eta _{{s},\infty }$
is the hypothetical viscosity of the fluid inside the solid far away from the solid–fluid interface. Since the latter is very difficult to obtain and, in general, will be much larger than the fluid viscosity, it is here set to a large positive value.
3. Case study and numerical methods
3.1. Set-up of the contact angle study
An exemplary snapshot of a moving droplet is provided in figure 1. The individual Lennard–Jones particles from NEMD are shown in figure 1(a), where the solid consists of Lennard–Jones particles in a lattice. The density profile of the fluid from hydrodynamic DFT is given in the lower part. An external force in the
$x$
-direction
$f_x$
is applied to induce the movement of the droplet and the surrounding fluid in both approaches. This is analogous to the gravitational force causing the movement of droplets on a vertical wall in macroscopic systems. Advancing and receding contact angles
$\Theta _{a}$
and
$\Theta _{r}$
were reported for NEMD simulations of droplets in similar systems (Hong et al. Reference Hong, Ha and Balachandar2009) and can be determined from the density profiles of hydrodynamic DFT. The equilibrium MD and NEMD simulations are carried out in three dimensions, where the average density does not change in the
$z$
-direction (a side view is provided in figure 1
a). Hydrodynamic DFT is simulated in a two-dimensional system with lengths
$L_x$
and
$L_y$
.
Since at the microscopic scale capillary forces have a much stronger influence on the droplet movement compared with the macroscopic scale, the external (body) force must be large compared with macroscopic driving forces such as gravity (see Appendix B). Here, external forces ranging from
$f_x={0.056}\,{\mathrm{pN}}$
to
$f_x={0.224}\,{\mathrm{pN}}$
per particle are chosen, which for methane, corresponds to an acceleration of approximately
$2.1\times 10^{12}$
and
$8.4\times 10^{12}\,{\mathrm{m\,s^{-2}}}$
, respectively. The dynamic contact angles take on values that can be observed for larger (macroscopic) droplets in earth’s gravity field. All parameters used for the generation of results are listed in table 1. The temperature is set to
$T={120.02}\,{\mathrm{K}}$
for all simulations in this work, which is well below the critical temperature, and a pressure of
$p={0.191}\,{\mathrm{MPa}}$
is obtained from the PC-SAFT model.
PC-SAFT, generalised entropy scaling and bulk viscosity model parameters for methane.
3.2.
Implementation of hydrodynamic DFT in DuMu
$^{\textit{x}}$
In the following, the numerical model for solving (2.7) is briefly described. For the discretisation of temporal and spatial differential operators, and for solving nonlinear and linear systems of equations, we rely on the open-source software package DuMu
$^{\textit{x}}$
(Koch et al. Reference Koch2021). DuMu
$^{\textit{x}}$
is based on the numerics framework DUNE (Bastian et al. Reference Bastian2021), and is a simulator for flow and transport processes, adaptable to various multiphysics problems. For the implementation of hydrodynamic DFT, its modular architecture allows DuMu
$^{\textit{x}}$
to be coupled with advanced thermodynamic models, such as DFT based on PC-SAFT, to predict molecular interactions that influence the fluid behaviour.
DFT calculations are conducted using FeO
$_{{s}}$
(Rehner & Bauer Reference Rehner and Bauer2021; Rehner et al. Reference Rehner, Bauer and Gross2023), i.e. the evaluation of the Helmholtz energy functional
$F$
and the functional derivative
${\delta F}/{\delta \rho }$
as well as generalised entropy scaling for the viscosity calculation. To couple DuMu
$^{\textit{x}}$
and FeO
$_{{s}}$
, which is written in the Rust programming language, we build a 
interface for the required FeO
$_{{s}}$
functions. This interface can then be dynamically linked and accessed directly from DuMu
$^{\textit{x}}$
.
For the discretisation of coupled partial differential equations as in (2.7), it is well known that collocated finite-volume methods or standard finite-element schemes lead to numerical instabilities. In this work, a staggered-grid finite-volume scheme (Harlow & Welch Reference Harlow and Welch1965) is applied, where densities are defined at element centres, while velocity components are placed on element faces, around which dual control volumes are constructed. Such staggering of degrees of freedom naturally leads to a stable scheme. Further details and a compact notation of the discrete equations can be found from Schneider et al. (Reference Schneider, Weishaupt, Gläser, Boon and Helmig2020). The DFT term
$\rho {\boldsymbol\nabla} ( {\delta F}/{\delta \rho } + V^{{\textit{ext}}})$
is discretised with a central difference approximation related to the dual control volumes, leading to its evaluation on element centres.
When using an implicit time discretisation scheme, such as an implicit Euler method, the convolutions appearing in the Helmholtz energy functional
$F$
lead to a non-local stencil for this DFT term and consequently to dense matrices when considering the related functional derivatives within nonlinear solvers. As a result, this term is treated semi-implicitly, i.e. the derivatives are not directly accounted for within the nonlinear solver (in this case, Newton’s method). This approach maintains the sparsity pattern of the matrices.
3.3. Non-equilibrium molecular dynamics simulations
Hydrodynamic DFT is assessed by comparing results for dynamic contact angles to NEMD simulations. Simulations are conducted using the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS, stable release 2 Aug 2023) (Thompson et al. Reference Thompson2022) in a configuration analogous to that employed in our previous development of generalised entropy scaling (Bursik et al. Reference Bursik, Eller and Gross2024b
). The equations of motion in the canonical ensemble (constant
$N,V,T$
) are integrated by means of a velocity-Verlet method. The time step is set to
$\Delta t^*=0.005t^*$
with
$t^*= {t}/{\sigma _{{\textit{ff}}} \sqrt {M / \varepsilon _{{\textit{ff}}}}}$
, where the asterisk
$*$
denotes dimensionless quantities. For the conversion from dimensionless numbers to real units, the PC-SAFT parameters of the fluid,
$\sigma _{{\textit{ff}}}$
and
$\varepsilon _{{\textit{ff}}}$
, and the molecular mass
$M$
are employed. Effectively, we compare results for methane from hydrodynamic DFT using the PC-SAFT model to results from NEMD for the Lennard–Jones fluid, where the PC-SAFT parameters are used to convert dimensionless results to real units.
Three Nose–Hoover chains (Tuckerman et al. Reference Tuckerman, Alejandre, López-Rendón, Jochim and Martyna2006) with a damping time constant
$t_{{D}}^*=100\Delta t^*$
are employed as a thermostat for the temperature in the system. The usage of a global thermostat in an inhomogeneous non-equilibrium system can lead to temperature variations due to large local shear stresses. For the system studied here, the temperature is found to be constant throughout the majority of the droplet, with temperature deviations of up to 10 % in the contact regions as discussed in Appendix F. The solid is modelled as a frozen rigid wall.
As described previously and visualised in figure 1, NEMD simulations are performed in three dimensions, whereas a two-dimensional approach is employed for hydrodynamic DFT. A sessile cylinder is simulated with both methods. The same length
$L_x$
and the same height
$L_y$
are chosen for NEMD and hydrodynamic DFT, respectively. The depth of the third-dimension in NEMD is chosen as
$L_z=L_x/3$
to balance computational time and statistics of the results. Periodic boundary conditions are applied in all directions for NEMD in accordance with hydrodynamic DFT. A solid block is introduced at the bottom of the system, which consists of Lennard–Jones sites at their equilibrium distance in a body-centred cubic (bcc) lattice with dimensions
$L_x \times 5\sigma _{{\textit{ff}}}\times L_z$
. Due to the periodic boundary conditions, fluid at the upper end of the system interacts with the lower end of the solid. The Lennard–Jones parameters for solid–solid and fluid–fluid interactions are set to
$\sigma _{{\textit{ff}}}=\sigma _{{ss}}$
and
$\varepsilon _{{\textit{ff}}}=\varepsilon _{{ss}}$
, which, for the solid–fluid interactions, leads to
$\sigma _{{\textit{sf}}}=\sigma _{{\textit{ff}}}=\sigma _{{ss}}$
using the Lorentz–Berthelot combining rules. However, the dimensionless solid–fluid energy interaction parameter
${\varepsilon _{{\textit{sf}}}^*} = {\varepsilon _{{\textit{sf}}}}/{\varepsilon _{{\textit{ff}}}}$
is varied as noted in the respective paragraphs or figures, which allows for the study of different wetting behaviours (see e.g. Becker et al. Reference Becker, Urbassek, Horsch and Hasse2014). The fluid–fluid and solid–fluid interactions are cutoff and shifted at
$r_{{c,ff}} = 5 \sigma _{{\textit{ff}}}$
and at
$r_{{c,\textit{sf}}} = 4 \sigma _{{\textit{ff}}}$
, respectively.
For the simulation of static contact angles, equilibrium MD simulations are performed. Following Becker et al. (Reference Becker, Urbassek, Horsch and Hasse2014), fluid particles were initially arranged in a cuboid on top of the solid. The initial density in the cuboid and in the gas need to be chosen such that a stable droplet is formed. These densities were taken from Vrabec et al. (Reference Vrabec, Kedia, Fuchs and Hasse2006), who studied droplets in a gas phase (without solid) and Lennard–Jones interactions with a cutoff at
$r_{{c,ff}} = 2.5 \sigma _{{\textit{ff}}}$
. While not representing exactly the same system as in this work, these densities proved to be sufficient as an initial estimate for obtaining stable sessile droplets. A conjugate gradient method and an equilibration of
$2.5\times 10^{5}$
steps were used to obtain an equilibrated sessile droplet. Production runs of
$4.75\times 10^{6}$
steps were carried out where atom positions were written out every 1000 steps to determine equilibrium/static density profiles and contact angles. Accumulation of net momentum due to numerical inaccuracies is avoided by removing the momentum in the
$x$
- and
$z$
-directions of the centre of mass every time step.
For dynamic contact angles, the equilibrium system is used as an initial condition in the NEMD simulations and an external driving force is applied to induce the movement of the droplet. The external driving force parallel to the solid–fluid interface
$f_x$
is added to the
$x$
-component of the force vector for each individual fluid atom in each time step. Therefore, the
$x$
-component of the particle velocities is not included in the calculation of the fluid temperature needed for thermostatting. Particle positions and velocities are written to a file every 100 time steps. All errorbars provided in this work denote the
$95\,\%$
confidence interval.
3.4. Contact angles from density profiles
The static and dynamic contact angles are determined from density profiles for both, hydrodynamic DFT and NEMD. For diffuse interface models, the exact location of the vapour–liquid interface of the droplet is not unambiguously defined. One possible approach is the Gaussian convolution method, which can be employed to determine the vapour–liquid interface from atomic coordinates (Willard & Chandler Reference Willard and Chandler2010). Here, we follow a procedure similar to Sauer et al. (Reference Sauer, Terzis, Theiss, Weigand and Gross2018) and Heier et al. (Reference Heier, Stephan, Diewald, Müller, Langenbach and Hasse2021), which is summarised in figure 2. We also note that Sauer et al. (Reference Sauer, Terzis, Theiss, Weigand and Gross2018) performed a finite size study for static contact angles and confirmed convergence towards the contact angle determined from Young’s equation. First, a dividing density
$\rho _{{\textit{iso}}}={11.44}\,\mathrm{kmol \,m^{-3}}$
is chosen. From this, iso-density points are obtained, which represent the vapour–liquid interface of the droplet (red and blue points in figure 2). Half-circles are fitted to it for each half of the droplet (red and blue lines). This is necessary as the density profiles of moving droplets are, in general, not axis-symmetric. The contact angle can be calculated from the tangent to each half of the circle (red and blue dashed lines) at the solid–fluid interface (green dash-dotted line). The molecular layering at the solid–liquid interface leads to oscillations in the density. As noted in the literature, an improvement of the fit can be obtained by disregarding the region very close to the solid–fluid interface (Sauer et al. Reference Sauer, Terzis, Theiss, Weigand and Gross2018; Heier et al. Reference Heier, Stephan, Diewald, Müller, Langenbach and Hasse2021). In addition, in the present geometry, the specific location of the solid–fluid interface is not uniquely defined due to the molecular roughness of the solid. Minor differences in the contact angle value were observed by varying either the location of the solid–fluid interface, disregarding a certain region close the interface or changing the dividing density
$\rho _{{\textit{iso}}}$
. Because this study is concerned with the comparison of hydrodynamic DFT and NEMD, we apply the same procedure for hydrodynamic DFT and NEMD results, and thereby ensure comparability between both modelling frameworks.
Visualisation of the methodology for determining contact angles from density profiles. The red and blue iso-density points represent the vapour–liquid interface of the droplet, the red and blue solid lines show the circular fit, the red and blue dashed lines are the tangents to the circle at the solid–fluid interface. The latter is shown as a green dashed line.
4. Results and discussion
In this section, we present results from equilibrium and hydrodynamic DFT for sessile droplets on a solid and validate the results by comparison with equilibrium and non-equilibrium MD simulations, respectively.
4.1. Equilibrium droplets
Equilibrium droplets from DFT and equilibrium MD are used as initial profiles for the dynamic simulations in hydrodynamic DFT and NEMD, respectively. To achieve good agreement of the dynamic simulations, it is essential that the equilibrium densities agree. Figure 3 visualises the two-dimensional density profiles of the droplets from equilibrium DFT (figure 3
a) and equilibrium MD (figure 3
b) for a solid–fluid energy interaction parameter
$\varepsilon _{{\textit{sf}}}^*=0.5$
. We note that the DFT calculations are performed in two dimensions, whereas for MD, a three-dimensional simulation was performed. In the MD simulations, the simulation box dimensions and the number of methane molecules were chosen to ensure that a stable or moving droplet adopts the shape of a cylindrical cap (entirely comparable to the geometry considered in the hydrodynamic DFT). Results from DFT are predictions without adjustable parameters. Sessile droplets with a finite contact angle are obtained in both cases and two important phenomena are observed. First, inside the droplet, alternating layers of varying density are observed along the normal direction of the solid–fluid interface. This molecular layering (as reported e.g. by Becker et al. Reference Becker, Urbassek, Horsch and Hasse2014; Sauer et al. Reference Sauer, Terzis, Theiss, Weigand and Gross2018; Lee & Müller-Plathe Reference Lee and Müller-Plathe2022) extends for several layers into the lower part of the droplet. Second, the density profile exhibits a wave-like inhomogeneity in the direction parallel to the solid–fluid interface, which is caused by the molecular roughness of the solid. This is true even outside of the droplets, which can be explained by the adsorption of particles from the vapour phase.
Density profile of equilibrium droplet (
$f_x=0$
) from (a) DFT and (b) equilibrium MD at
$T={120.02}\,{\mathrm{K}}$
with
$\varepsilon _{{\textit{sf}}}^*=0.5$
.
The droplets are nearly symmetrical and the small asymmetry is caused by the solid roughness. The mean contact angles from DFT are in good agreement with equilibrium MD with
${101.4}{^{\circ }}$
and
${104.4}{^{\circ }}$
, respectively. A similar minor underestimation of contact angles from DFT compared with molecular simulations was previously reported (Sauer et al. Reference Sauer, Terzis, Theiss, Weigand and Gross2018). Several reasons may contribute to it: The PC-SAFT model (here with methane parameters) used in DFT does not exactly reproduce the properties of the Lennard–Jones fluid, which is employed in MD (Sauer & Gross Reference Sauer and Gross2017). In addition, the different treatment of solid–fluid interactions (explicit pair interactions versus external potential) and the different dimensionalities of the simulations (two- versus three-dimensional) may add to the deviation. Overall, we find very satisfactory agreement of DFT compared with MD calculations, both, regarding density profiles and contact angles, thus providing a suitable basis for the investigation of dynamic phenomena.
4.2. Importance of local viscosity model
This section assesses the role of a local viscosity model. The goal is to ascertain whether incorporating molecular details in the viscosity model by using generalised entropy scaling is essential for accurately modelling wetting phenomena or if continuum models are sufficient. Therefore, we compare generalised entropy scaling to a different viscosity model, which represents a typical continuum fluid dynamics approach and is based on viscosities from bulk phases (thus called ‘bulk’ model).
Viscosity profiles from (a) local entropy scaling model and (b) bulk viscosity model at
$f_x={0.112}\,{\mathrm{pN}}$
per particle and
$T={120.02}\,{\mathrm{K}}$
after
$t={1000}\,{\mathrm{ps}}$
.
The bulk viscosity model uses viscosities from bulk phases depending on the local density
$\rho (\boldsymbol{r})$
according to
\begin{equation} \eta ^{{bulk}}({r}) = \begin{cases} \eta _{l}, & \rho ({r})\gt 0.5(\rho _{l}+\rho _{{v}})+\rho _{{v}},\\ \eta _{{v}}, & \rho ({r})\lt = 0.5(\rho _{l}+\rho _{{v}})+\rho _{{v}},\\ \eta _{{s}}, & \rho ({r})\lt \rho _{{cut}},\\ \end{cases} \end{equation}
where
$\rho _{l}$
and
$\rho _{v}$
are densities from a vapour–liquid equilibrium at
$T={120.02}\,{\mathrm{K}}$
, and the corresponding bulk viscosities
$\eta _{l}$
and
$\eta _{v}$
are determined using entropy scaling for homogeneous phases. The procedure is as follows: first, the density of the liquid and the vapour phase are determined from the phase equilibrium of pure methane at
$T={120.02}\,{\mathrm{K}}$
. The residual entropy is then determined from the PC-SAFT model as the partial derivative of the Helmholtz energy with respect to the temperature. Using the appropriate dimensionless form of the entropy, the value for
$\ln (\eta ^{{ES}}/\eta _{{ref}})$
is determined from the ansatz function in (2.12), where the parameters
$A$
,
$B$
,
$C$
and
$D$
were previously adjusted to experimental data (Lötgering-Lin et al. Reference Lötgering-Lin, Fischer, Hopp and Gross2018). The Chapman–Enskog reference viscosity
$\eta _{\textit{CE}}$
is calculated according to (2.14), which finally allows to calculate the viscosity
$\eta _{l}$
and
$\eta _{v}$
, respectively. A very large viscosity
$\eta _{{s}}$
is chosen for the liquid within the solid domain, which is defined by a cutoff density
$\rho _{{cut}}$
, to ensure comparability with the generalised entropy scaling approach, but the impact of the viscosity in this region on velocity profiles is negligible. All parameters are provided in table 1.
The resulting viscosity profile is compared with that from generalised entropy scaling in figure 4. The viscosity from entropy scaling (figure 4 a) exhibits a more elaborate behaviour with strong oscillations perpendicular but also parallel to the solid–fluid interface. This is expected since the viscosity from entropy scaling is a nonlinear and non-local function of the density. In contrast, the viscosity from the bulk model (figure 4 b) assumes only two values in the vapour–liquid domain and a different constant value in the solid domain. While the entropy scaling model provides a continuous transition between liquid-like and vapour-like viscosity values, discontinuities are found at the transition between the different regions for the bulk viscosity model.
Velocities inside the droplet from hydrodynamic DFT using (a) the generalised entropy scaling viscosity model and (c) the bulk viscosity model as well as from (b) NEMD with the solid as the frame of reference for
$f_x={0.112}\,{\mathrm{pN}}$
per particle and
$T={120.02}\,{\mathrm{K}}$
averaged over at least
${700}\,{\mathrm{ps}}$
after a steady state is reached. Arrows denote the direction and magnitude of the velocity, whereas the colours correspond to the magnitude of the velocity.
To assess the accuracy of the viscosity models, we compare velocity profiles determined from hydrodynamic DFT using the two viscosity models and validate the results by NEMD simulations as provided in figure 5 for the medium force
$f_x={0.112}\,{\mathrm{pN}}$
, where the droplet moves into the positive
$x$
-direction (to the right). The velocities are averaged over at least 700 ps in hydrodynamic DFT and 10 000 ps in NEMD after a steady state is reached. For the NEMD results, in addition to averaging velocities in the third dimension, the centre of mass of the droplet is determined for each configuration entering the calculation of the velocity profile. Using coordinates relative to this centre of mass, substantial time-averaging of the velocity is performed to obtain a reasonable signal-to-noise ratio. For the initial acceleration of the droplet, the velocity results have considerable uncertainty and we do not report velocities for the acceleration period of the droplet. Noise-free velocity profiles are a significant advantage of hydrodynamic DFT over NEMD. For all three cases, the largest velocity relative to the solid is found at the top of the droplet and the smallest velocity at the solid–fluid interface. Good quantitative agreement is observed between hydrodynamic DFT with generalised entropy scaling (figure 5
a) and NEMD (figure 5
b). For the bulk viscosity model (figure 5
c), significantly lower velocities are found throughout the droplet, where especially the velocity close to the solid–fluid interface in the centre of the droplet is much lower than in NEMD and hydrodynamic DFT with generalised entropy scaling. The substantially different behaviour of the viscosity in the first molecular layers at the solid–fluid interface for the bulk and entropy scaling viscosity models (cf. figure 4) explains these results. This emphasises that molecular details, which are captured in the entropy scaling viscosity model, have a severe effect on the microscopic wetting behaviour.
Distance travelled by the centre of mass of the droplet from hydrodynamic DFT (HDFT) with entropy scaling viscosity model (blue line) and bulk viscosity model (light blue line) as well as from NEMD (red crosses) for
$f_x={0.112}\,{\mathrm{pN}}$
per particle and
$T={120.02}\,{\mathrm{K}}$
.
Furthermore, the influence of the viscosity models on the velocity of the entire droplet is analysed by plotting the distance travelled by the centre of mass of the droplet versus the simulation time (figure 6). In the initial phase of the simulation until approximately 100 ps, the distance covered by the droplet increases with increasing slope, i.e. the droplet moves with increasing velocity for all three models. After this initial phase, the distance travelled by the droplet increases almost linearly, which corresponds to a constant velocity
$v^{{avg}}_{{drop}}$
of the droplet and shows that a steady state is reached. At the very beginning of the simulation (until approximately 50 ps), the velocity of the droplet is small and results from all models agree. After approximately 50 ps, the results from the bulk viscosity model increasingly deviate from NEMD results, whereas good agreement between hydrodynamic DFT with generalised entropy scaling and NEMD is obtained.
Steady-state velocity of the centre of mass of the moving droplet for different external forces (per particle) from hydrodynamic DFT (HDFT) with entropy scaling viscosity model (blue circles) and bulk viscosity model (light blue triangles) as well as from NEMD (red crosses) at
$T={120.02}\,{\mathrm{K}}$
.
The steady-state velocities of the droplets
$v^{{avg}}_{{drop}}$
are summarised in figure 7 for different forces. For all models, the steady-state velocities increase with increasing force. For the low and medium forces (
$f_x={0.056}\,{\mathrm{pN}}$
and
$f_x={0.112}\,{\mathrm{pN}}$
), the velocity of the droplet predicted from hydrodynamic DFT with generalised entropy scaling accurately represents the results from NEMD. A small underestimation is observed for the largest force
$f_x={0.224}\,{\mathrm{pN}}$
, noting however that the statistical uncertainty of the velocity from NEMD is also large compared with the lower forces. In contrast, the bulk viscosity model substantially underestimates the steady-state velocity for all forces, where the underestimation is largest for the largest force
$f_x={0.224}\,{\mathrm{pN}}$
. This increasing influence of the viscosity at larger forces can be explained by the presence of larger velocities, which in turn lead to larger viscous dissipation.
Compared with macroscopic systems, our study investigates high velocity-gradients that are caused by large external forces needed to overcome the surface tension on the atomistic scale. The transition to macroscopic systems could best be investigated using a size-scaling study, where an extrapolation to macroscopically large droplets is performed. We note however that our study does not reveal any indication about nonlinear shear–strain relations. Our results are, in that sense, also applicable to low deformation tensors and low stress tensors. We conclude that using generalised entropy scaling provides more accurate results than the simpler model based on bulk viscosities. This is likely due to the importance of molecular details close to the solid–fluid interface, which are captured by entropy scaling. Thus, in the remainder of this work, we present results only for hydrodynamic DFT using the generalised entropy scaling model and we will exclude the bulk viscosity model.
4.3. Droplet movement: slip and rolling motion
In this section, the velocity profiles of droplets (some of which were already presented in figure 5) are discussed in more detail to evaluate whether hydrodynamic DFT with generalised entropy scaling captures mechanisms of wetting. Studies of microscopic models, such as diffuse interface models or MD simulations, usually report slip, i.e. a non-vanishing velocity at the solid–fluid interface in the contact region (Yue et al. Reference Yue, Zhou and Feng2010; Yue & Feng Reference Yue and Feng2011; Nakamura et al. Reference Nakamura, Carlson, Amberg and Shiomi2013; Blake et al. Reference Blake, Fernandez-Toledano and Doyen2015; Fernández-Toledano et al. Reference Fernández-Toledano, Blake and Limat2019; Lee & Müller-Plathe Reference Lee and Müller-Plathe2022). Furthermore, it was reported from MD simulations (Li et al. Reference Li, Yan, Fichthorn and Yu2018; Fernández-Toledano et al. Reference Fernández-Toledano, Blake and Limat2019; Li & Zhang Reference Li and Zhang2019) that microscopic droplets move along the solid–fluid interface by a ‘rolling’ mechanism and that this mechanism requires the presence of slip within the contact region.
Velocity relative to centre of mass velocity of the droplet from (a) hydrodynamic DFT and (b) NEMD for
$f_x={0.112}\,{\mathrm{pN}}$
per particle and
$T={120.02}\,{\mathrm{K}}$
averaged over at least
${700}\,{\mathrm{ps}}$
after a steady state is reached. Arrows denote the direction and magnitude of the flow, whereas the colours correspond to the
$y$
-component of the velocity.
The profiles of velocity with the solid as the reference determined from hydrodynamic DFT and NEMD (figure 5
a,b) reveal that the velocity is close to zero at the solid–fluid interface in the centre of the droplet. However, it is non-zero in the solid–liquid–vapour contact regions. This is an important observation as it demonstrates that slip is found in hydrodynamic DFT and in NEMD. In addition, figure 8 provides the velocity vectors relative to the velocity of the droplet’s centre of mass from hydrodynamic DFT and NEMD. The
$y$
-component of these vectors is visualised by colours, where positive and negative values correspond to flow away from and towards the solid, respectively. The translational motion of the droplet parallel to the solid–fluid interface is superimposed by a clearly visible circular motion as also reported in the literature (Li et al. Reference Li, Yan, Fichthorn and Yu2018; Fernández-Toledano et al. Reference Fernández-Toledano, Blake and Limat2019). According to figure 8, the maximum and minimum values of the relative velocity in the
$y$
-direction are approximately
${25}$
and
${- 40}\,{\mathrm{m\,s^{-1}}}$
, respectively. Comparing these values to the maximum values of the absolute velocity in figure 5, which are in the range of
${100}\,{\mathrm{m\,s^{-1}}}$
, demonstrates that the circular motion is quantitatively significant. Furthermore, analysing the vorticity of the flow field supports the finding that the droplet exhibits a significant rolling motion (see Appendix E). The presented results demonstrate that hydrodynamic DFT predicts the complex mechanisms of droplet motion at the molecular scale, particularly slip in the contact region and the circular motion of droplets, in good agreement with our NEMD simulations. We emphasise that noise-free velocity profiles are obtained from hydrodynamic DFT, whereas in NEMD, substantial averaging is required to generate meaningful velocity profiles. In addition, the good agreement regarding droplet motion further supports the validity of the generalised entropy scaling viscosity model.
4.4. Advancing and receding contact angles
The appearance of different dynamic contact angles is commonly attributed to inhomogeneities of the solid (Voinov Reference Voinov1976; Cox Reference Cox1986; Decker & Garoff Reference Decker and Garoff1997; Jamet et al. Reference Jamet, Lebaigue, Coutris and Delhaye2001; de Gennes et al. Reference de Gennes, Brochard-Wyart and Quere2004; Bonn et al. Reference Bonn, Eggers, Indekeu, Meunier and Rolley2009; Butt et al. Reference Butt2022) and they were found to depend on the velocity of the contact region (Voinov Reference Voinov1976; Cox Reference Cox1986). In the literature (Butt et al. Reference Butt2022), the term dynamic contact angle hysteresis is sometimes used for the difference between advancing and receding contact angle. Since in a strict thermodynamic sense, the term hysteresis is disputable in this context, we will speak of the difference between advancing and receding contact angles instead. The microscopic mechanisms behind this effect are studied in several works (Ren & E Reference Ren and W.2007; Bertrand, Blake & Coninck Reference Bertrand, Blake and Coninck2009; Blake et al. Reference Blake, Fernandez-Toledano and Doyen2015; Lukyanov & Likhtman Reference Lukyanov and Likhtman2016; Lukyanov & Pryer Reference Lukyanov and Pryer2017), where a particularly convincing explanation is based on a solid–fluid friction force (Lukyanov & Likhtman Reference Lukyanov and Likhtman2016). This friction force occurs mainly in the contact region and depends on its velocity, since it is caused by dynamic interactions of fluid molecules with the solid interaction sites.
From a molecular dynamics point of view, the solid–fluid friction in the contact region can be explained by a net force counteracting the flow of fluid particles. Due to the net flow in the positive
$x$
-direction (parallel to the solid surface), fluid particles are more likely to enter regions where they experience a repulsive force from the solid atoms acting in the direction opposite to the flow. In hydrodynamic DFT, the equivalent mechanism is captured by the external potential (cf. § 2.1). The gradient of the external potential multiplied by the density enters the momentum balance in (2.7b
). Close to the location of solid atoms, the gradient of the external potential in the
$x$
-direction is negative and, thus, effectively reduces the momentum in the
$x$
-direction, which in turn results in the solid–fluid friction force. Its velocity dependence results from the fact that higher velocities in the
$x$
-direction lead to larger densities in the regions where the gradient of the external potential reduces the momentum in the
$x$
-direction. Thus, the product of density and external potential gradient increases with velocity, which leads to a larger reduction of the momentum in the
$x$
-direction. A visual representation of the gradient of the external potential is provided in Appendix D. The molecular-kinetic theory (MKT), as opposed to our work, assumes a certain mechanism causing this solid–fluid friction force (Blake & Haynes Reference Blake and Haynes1969; Blake Reference Blake, Wettability and Berg1993): the fluid molecules perform random displacement (jumps) of average distance, which occur with a certain temperature-dependent frequency. The solid surface is thought of as a series of potential energy wells, at which the molecules remain for a short period between these displacements.
Density profiles of droplets moving along the solid–fluid interface with external force
$f_x={0.112}\,{\mathrm{pN}}$
per particle and
$\varepsilon _{{\textit{sf}}}^*=0.5$
from hydrodynamic DFT (HDFT) and NEMD at different simulation times
$t$
. NEMD results show statistical noise, whereas hydrodynamic DFT provides deterministic density profiles.
We compare results from hydrodynamic DFT and NEMD for the transition between the equilibrium and a dynamic steady state to analyse if differences between the advancing and receding contact angles are observed. Figure 9 compares density profiles from hydrodynamic DFT and NEMD at different times after the start of the simulation for an external force of
$f_x={0.112}\,{\mathrm{pN}}$
. The results from NEMD show significant statistical noise, since at each time, a limited number of particles located at discrete positions enter the calculation of the density profiles (compare figure 1
a). The density profiles from NEMD are smoothed following Hong et al. (Reference Hong, Ha and Balachandar2009), and all profiles are relocated such that the centre of mass is positioned at
$L_x/2$
. Starting from the equilibrium density profile (see figure 3), the droplet from hydrodynamic DFT begins to deform after 20 ps, as shown in figure 9(a). With increasing simulation time, the deformations become more pronounced (see figures 9
c, 9
e and 9
g). The molecular layering and adsorption from the gas are observed similarly to the equilibrium case. Advancing contact angles increase while receding contact angles decrease compared with the static contact angle. The results from NEMD (figures 9
b, 9
d, 9
f and 9
h) exhibit very similar behaviour, even though the analysis is obscured by fluctuations of the density profile. The previous explanation for the difference between advancing and receding contact angles, and its velocity dependence based on a solid-fluid friction force also applies to the transition behaviour in figure 9. After 100 ps, the change in droplet shape is small (not shown), which suggests that a steady state is approached, where the external force is in balance with the solid–fluid friction and viscous forces.
The occurrence of a solid–fluid friction force and, thus, of differences between the advancing and receding contact angle are commonly attributed to inhomogeneities of the solid surface. Our results demonstrate that these differences appear even for inhomogeneities at the molecular scale, which is consistent with findings of other studies (Bertrand et al. Reference Bertrand, Blake and Coninck2009; Blake et al. Reference Blake, Fernandez-Toledano and Doyen2015; Lukyanov & Likhtman Reference Lukyanov and Likhtman2016; Fernández-Toledano et al. Reference Fernández-Toledano, Blake and Limat2019). The molecular roughness of the solid introduced by the solid atoms that are arranged in a perfect lattice structure is sufficient to cause different advancing and receding contact angles. The presented results emphasise that hydrodynamic DFT is capable of predicting differences between advancing and receding contact angles at the molecular scale. Furthermore, the dynamic behaviour along the path from an initial equilibrium profile to the dynamic steady state is accurately described by hydrodynamic DFT. This includes the interplay between external forces, velocity-dependent solid–fluid friction and viscous forces in the contact region. These findings indicate that generalised entropy scaling provides meaningful values for the viscosity.
4.5. Influence of external force
For steady states, larger external forces lead to an increase of contact region velocity, which in turn affects microscopic dynamic contact angles as discussed previously and in the literature (Ren & E Reference Ren and W.2007; Bertrand et al. Reference Bertrand, Blake and Coninck2009; Blake et al. Reference Blake, Fernandez-Toledano and Doyen2015; Lukyanov & Likhtman Reference Lukyanov and Likhtman2016; Lukyanov & Pryer Reference Lukyanov and Pryer2017; Fernández-Toledano et al. Reference Fernández-Toledano, Blake and Limat2019). The general finding is that advancing contact angles increase and receding contact angles decrease compared with the equilibrium contact angle.
In figure 10, results for the density profiles of moving droplets from hydrodynamic DFT using the generalised entropy scaling viscosity model and NEMD simulations are depicted for different external forces (results for the bulk viscosity model are provided in the supplementary material). Density profiles were averaged over at least 300 ps after a steady state was reached, for which the centre of mass of the droplets in each snapshot was shifted to position
$L_x/2$
. This postprocessing is required for the NEMD results to average out fluctuations, whereas hydrodynamic DFT provides noise-free results and averaging is performed solely for better comparability. Due to this averaging, inhomogeneities on the solid–fluid interface are not visible in contrast to the equilibrium droplets (see figure 3) and the snapshots of the accelerating droplets (see figure 9).
Density profiles of droplets moving along the solid–fluid interface with different external forces
$f_x$
from hydrodynamic DFT (HDFT) and NEMD at
$T={120.02}\,{\mathrm{K}}$
with
$\varepsilon _{{\textit{sf}}}^*=0.5$
averaged over
${700}\,{\mathrm{ps}}$
after a steady state is reached.
For the lowest external force given in figures 10(a) and 10(b), a difference between advancing and receding contact angles is observed for both hydrodynamic DFT and NEMD. In particular, advancing contact angles are
${109.3}{^{\circ }}$
and
${116.5}{^{\circ }}$
, while receding contact angles are found as
${82.2}{^{\circ }}$
and
${87.6}{^{\circ }}$
from hydrodynamic DFT and NEMD, respectively. The small underestimation of contact angles, which was identified for the equilibrium droplet, transfers to the dynamic case. Furthermore, the viscosity coefficients obtained from generalised entropy scaling and the inherent approximations of hydrodynamic DFT (e.g. adiabatic approximation) might introduce additional error. Nevertheless, these results demonstrate good agreement between hydrodynamic DFT and NEMD for
$f_x={0.056}\,{\mathrm{pN}}$
, qualitatively in terms of the droplet shape given by the density profiles and quantitatively regarding the advancing and receding contact angles.
For an increasing external force (see figure 10 c–f), the deformation of the droplets and, consequently, the difference between advancing and receding contact angle become more pronounced. The density profiles show that the droplets flatten and their width increases. The adsorption at the solid–vapour interface also becomes stronger with increasing external force. Remarkably, this subtle effect is captured by hydrodynamic DFT.
Figure 11 summarises the contact angles in the system with
$\varepsilon _{{\textit{sf}}}^*=0.5$
for the different external forces. In addition, the relation between contact angles and the velocity of the contact region is visualised and discussed in the supplementary material. Interestingly, in the system studied here, while the advancing contact angle increases from the equilibrium case to the lowest external force
$f_x={0.056}\,{\mathrm{pN}}$
, it does not increase with further increases in the force, but rather approaches values similar to the equilibrium contact angle. In contrast, the receding dynamic contact angle monotonically decreases with increasing external forces. The dependence of the dynamic contact angle and droplet shape on the external force is captured by the hydrodynamic DFT in good agreement with NEMD simulations with a small systematic underestimation of the contact angles in all results. We expect that this underestimation could be further reduced by improving the agreement in the equilibrium case.
Summary of advancing and receding contact angles from hydrodynamic DFT (blue points) and NEMD (red crosses) for different external forces (per particle) at
$T={120.02}\,{\mathrm{K}}$
and with
$\varepsilon _{{\textit{sf}}}^*=0.5$
.
4.6. Influence of wetting strength
The wetting strength significantly affects the contact angles of droplets. It is modelled in NEMD and hydrodynamic DFT by altering the solid–fluid energy interaction parameter
$\varepsilon _{{\textit{sf}}}^*$
. A large value for
$\varepsilon _{{\textit{sf}}}^*$
corresponds to strong wetting and small contact angles. In the following, results for stronger wetting, modelled with
${\varepsilon _{{\textit{sf}}}^*}=0.7$
, are presented in comparison to the results provided previously, where
${\varepsilon _{{\textit{sf}}}^*}=0.5$
was used. The goal is to analyse if the hydrodynamic DFT model captures the influence of the wetting strength on dynamic contact angles and droplet velocities.
Density profiles with
$\varepsilon _{{\textit{sf}}}^*=0.7$
for the equilibrium and dynamic case are shown in figure 12. We remind that for equilibrium and hydrodynamic DFT a two-dimensional representation is employed, whereas (NE)MD is simulated in three dimensions. The equilibrium density profile from DFT (figure 12
a) confirms that increasing the solid–fluid energy interaction parameter to
${\varepsilon _{{\textit{sf}}}^*}=0.7$
leads to a flatter equilibrium droplet with a lower mean contact angle of
${62.6}{^{\circ }}$
compared with the results with
${\varepsilon _{{\textit{sf}}}^*}=0.5$
, where a mean contact angle of
${101.4}{^{\circ }}$
was determined (cf. figure 3
a). In addition, the molecular layering at the solid–fluid interface is more pronounced, which is caused by the stronger solid–fluid interactions. The density profile and contact angle from DFT is in good agreement with results from equilibrium MD (figure 12b
). Analogously to the results for
${\varepsilon _{{\textit{sf}}}^*}=0.5$
, the mean contact angle from DFT is slightly smaller than from MD (
${62.6}{^{\circ }}$
versus
${67.1}{^{\circ }}$
).
Figure 12(c) visualises the density profile from hydrodynamic DFT for the medium force (
$f_x={0.112}\,{\mathrm{pN}}$
). The droplet deforms strongly and the advancing contact angle differs from the receding one. This difference
$\Theta _{{a}}-\Theta _{{r}}$
is larger for the case with stronger wetting (
${44.3}{^{\circ }}$
for
${\varepsilon _{{\textit{sf}}}^*}=0.7$
versus
${34.5}{^{\circ }}$
for
${\varepsilon _{{\textit{sf}}}^*}=0.5$
). This might be due to a larger solid–fluid friction at higher wetting strengths, as a result of stronger solid–fluid interactions in the contact region. At the largest force (figure 12
e), the droplet elongates and the tendency to separate into two smaller droplets can be observed. This behaviour was not observed in the previous results (cf. figure 10
e) and meaningful dynamic contact angles cannot be obtained by adjusting a half-circle to the density profiles. This shape deformation, as predicted from hydrodynamic DFT, agrees well with the results from NEMD simulations, as provided in figures 12(d) and 12(f).
Density profiles of droplets moving along the solid–fluid interface with different external forces
$f_x$
from (a) equilibrium and hydrodynamic DFT and (b) equilibrium and non-equilibrium MD at
$T={120.02}\,{\mathrm{K}}$
with
$\varepsilon _{{\textit{sf}}}^*=0.7$
averaged over
${700}\,{\mathrm{ps}}$
after a steady state is reached.
Summary of advancing and receding contact angles for varying solid–fluid interaction parameter
$\varepsilon _{{\textit{sf}}}^*$
determined for different external forces (per particle) from hydrodynamic DFT (HDFT) with entropy scaling viscosity model (circles) and from NEMD (crosses) at
$T={120.02}\,{\mathrm{K}}$
.
Contact angles determined from these density profiles are summarised and contrasted to previous results for
${\varepsilon _{{\textit{sf}}}^*}=0.5$
in figure 13. For both wetting strengths, similar qualitative trends are observed from hydrodynamic DFT and NEMD: the receding contact angles decrease in comparison to the equilibrium contact angles, whereas the advancing contact angles first increase and then approach values similar to the equilibrium contact angles. The small underestimation of the contact angles from DFT and hydrodynamic DFT, which appears in the equilibrium case and transfers to the dynamic case, is observed for both wetting strengths. For the larger wetting strength, this effect vanishes at larger forces. Density profiles and contact angles from hydrodynamic DFT and NEMD are in good quantitative agreement for all forces investigated for both wetting strengths.
Steady-state velocity of the centre of mass of the moving droplet for varying solid–fluid interaction parameter
$\varepsilon _{{\textit{sf}}}^*$
determined for different external forces (per particle) from hydrodynamic DFT (HDFT) with the entropy scaling viscosity model (circles) and from NEMD (crosses) at
$T={120.02}\,{\mathrm{K}}$
.
According to figure 14, the steady-state velocity of droplets for stronger wetting (
${\varepsilon _{{\textit{sf}}}^*}=0.7$
) is much smaller compared with weaker wetting (
${\varepsilon _{{\textit{sf}}}^*}=0.5$
). Furthermore, the steady-state velocity increases less strongly between the second largest and largest forces. This is consistent with the more pronounced difference between advancing and receding contact angle for the stronger wetting case, since both can be explained by a larger solid–fluid friction. At the largest force, the solid–fluid friction becomes very large leading to a small increase in steady-state velocity and the tendency to separate into smaller droplets. The NEMD results confirm the accuracy of steady-state velocities from hydrodynamic DFT.
These results show that hydrodynamic DFT accurately captures the influence of the wetting strength on the wetting behaviour. Importantly, this also illustrates the transferability of the adjustable parameter in the generalised entropy scaling approach.
4.7. Influence of molecular roughness of the solid
In addition to the wetting strength, the molecular roughness of the solid is a major influence factor for the wetting behaviour, as it determines the energetic landscape of the solid–fluid interface, which causes solid–fluid friction (Lukyanov & Likhtman Reference Lukyanov and Likhtman2016). We employ the term molecular roughness to emphasise that this roughness is on the molecular scale and, thus, not directly comparable to the macroscopic roughness of a solid. We use the height difference
$h$
between the top and bottom layer of solid atoms exposed to the fluid as a measure for the molecular roughness of the solid, as depicted in figure 15. By removing solid atoms (visualised in blue) from the lattice, this roughness parameter is increased from
$h=0.5l_{{cell}}$
(all previous results) to
$h=1.5l_{{cell}}$
(results presented in the following), where
$l_{{cell}}$
is the length of a unit cell of the bcc lattice.
Visualisation of the methodology for modelling different molecular roughnesses of the solid. The smoother solid (previous results) contains red and blue atoms; the increased roughness is obtained by removing the blue atoms from the solid. The height difference
$h$
between the top and lowest layers which are in contact with the fluid is determined as multiples of the length of a unit cell
$l_{{cell}}$
.
Density profiles of droplets moving along the solid–fluid interface for an increased molecular roughness of the solid (
$h=1.5l_{{cell}}$
) with different external forces
$f_x$
from (a) equilibrium and hydrodynamic DFT as well as (b) equilibrium and non-equilibrium MD at
$T={120.02}\,{\mathrm{K}}$
with
$\varepsilon _{{\textit{sf}}}^*=0.5$
averaged over
${700}\,{\mathrm{ps}}$
after a steady state is reached.
Density profiles for the rougher surface (
$h=1.5l_{{cell}}$
) are presented in figure 16. In the equilibrium density profiles (figures 16
a and 16
b), the increased roughness is clearly visible in the density profiles. While the droplet from DFT is slightly less high, its shape and the determined contact angles agree sufficiently well with results from equilibrium MD. The difference in the shape of the droplets can be explained by the averaging of the external potential, which is required due to the different dimensionalities of the models (two-dimensional hydrodynamic DFT versus three-dimensional MD) and has a stronger influence on the results for an increased solid roughness. The contact angles from DFT and equilibrium MD are slightly smaller compared with the previous results with
${98.0}{^{\circ }}$
and
${96.1}{^{\circ }}$
for
$h=1.5l_{{cell}}$
versus
${101.4}{^{\circ }}$
and
${104.4}{^{\circ }}$
for
$h=0.5l_{{cell}}$
. On the macroscopic scale, the influence of solid roughness on the contact angle can be correlated using the Wenzel equation, while on the microscopic scale, its validity remains to be determined.
In the dynamic density profiles (figure 16
c–f), the solid roughness is not explicitly visible due to the averaging of density profiles over time, but it leads to a larger gas adsorption, which can be observed in all cases. For the medium force (
$f_x={0.112}\,{\mathrm{pN}}$
), the density profiles from hydrodynamic DFT (figure 16
c) and NEMD (figure 16
d) are in very good agreement. While for the strongest force (figures 16
e and 16
f), the agreement with NEMD is still sufficient, the receding contact region is smeared out more strongly in NEMD than in hydrodynamic DFT. As in the equilibrium case, this can be explained by the different dimensionalities of the models.
Summary of advancing and receding contact angles for varying solid roughness
$h$
determined for different external forces (per particle) from hydrodynamic DFT (HDFT) with entropy scaling viscosity model (circles) and from NEMD (crosses) at
$T={120.02}\,{\mathrm{K}}$
.
Contact angles for the different solid roughness are compared in figure 17. Similar values for the contact angles are obtained for the different roughness throughout all forces. For the cases studied here (
$h=0.5l_{{cell}}$
and
$h=1.5l_{{cell}}$
), the solid roughness apparently does not significantly affect the contact angles. Hydrodynamic DFT captures this effect in agreement with NEMD simulations.
A different behaviour is found for the steady-state velocity, as depicted in figure 18. Smaller steady-state velocities are found for the increased solid roughness for all forces, where results from hydrodynamic DFT and NEMD are in very good agreement. The reduced steady-state velocity can be explained by increasing attractive solid–fluid interactions acting in the direction opposite of the droplet motion as a result of the changed geometry when increasing the solid roughness. This increases the solid–fluid friction and leads to smaller steady-state velocities. The hydrodynamic DFT model correctly describes the influence of the molecular roughness of the solid on contact angles and steady-state velocities. Notably, no information on different molecular roughnesses of the solid entered the adjustment of the parameter for the generalised entropy scaling approach.
Steady-state velocity of the centre of mass of the moving droplet for varying solid roughness
$h$
determined for different external forces (per particle) from hydrodynamic DFT (HDFT) with entropy scaling viscosity model (circles) and from NEMD (crosses) at
$T={120.02}\,{\mathrm{K}}$
.
5. Conclusion
In this work, we have conducted a two-dimensional and quantitative investigation of the predictive capabilities of hydrodynamic DFT with respect to wetting at the microscopic scale by studying droplets moving along the solid–fluid interface. The model captures the influence of fluid–fluid interfaces on the dynamics by using a DFT term and the influence of the solid by an external potential. Generalised entropy scaling, an approach which incorporates microscopic detail, is employed for the viscosity, where one transferable parameter is adjusted to a single NEMD simulation of liquid-phase Poiseuille flow, i.e. no wetting information enters the model. We have evaluated the hydrodynamic DFT model through comparison to results from NEMD simulations considering velocity profiles and steady-state velocities, as well as density profiles and resulting dynamic contact angles.
The study demonstrated good quantitative agreement between hydrodynamic DFT and NEMD simulations for a wide range of wetting phenomena and yields three major findings: First, the microscopic generalised entropy scaling model for the local viscosity using a single transferable parameter, which does not include wetting data, ensures good quantitative agreement with NEMD results and more accurate results than a bulk viscosity model representing a typical continuum approach. Second, wetting phenomena at the microscopic scale, including differences between advancing and receding contact angles, the transition from equilibrium to steady state and the rolling motion of droplets, can be predicted with hydrodynamic DFT. Third, hydrodynamic DFT combined with generalised entropy scaling is transferable and captures the influence of different external forces, wetting strengths and molecular roughness of the solid on the wetting behaviour.
We conclude that hydrodynamic DFT, a unified molecular and continuum mechanics approach, is capable of predicting wetting phenomena at the microscopic scale, while reducing to the Navier–Stokes equations far away from interfaces. It is therefore, a suitable candidate for the study of wetting on multiple scales, for example, with application to mixtures in (multiphase) transport in porous media.
Funding
This research is funded in part by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – Project Number 327154368 – SFB 1313. We thank the German Research Foundation (DFG) for supporting this work by funding EXC 2075/1-390740016 under Germany’s Excellence Strategy as well as the Center for Digitalization and Technology Research of the Armed Forces of Germany (dtec.bw) through the project Macro/Micro-Simulation of Phase Decomposition in the Transcritical Regime (MaST); dtec.bw is funded by the European Union–NextGenerationEU. We acknowledge support by the Stuttgart Center for Simulation Science (SC SimTech). The authors acknowledge support by the state of Baden-Württemberg through bwHPC.
Supplementary material
Supplementary material is available at https://doi.org/10.1017/jfm.2025.10285.
Declaration of interests
The authors have no conflict to disclose.
Data availability statement
The data that support the findings of this study are openly available in the data repository of the University of Stuttgart (DaRUS) at https://doi.org/10.18419/DARUS-4528.
Author contributions
B.B.: Conceptualisation (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (equal); Software (equal); Validation (lead); Visualisation (lead); Writing – original draft (lead). R.S.: Conceptualisation (equal); Data curation (supporting); Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Visualisation (equal); Writing – original draft (equal). H.O.: Methodology (equal); Software (equal); Writing – original draft (supporting). M.S.: Funding acquisition (equal); Methodology (equal); Software (equal); Writing – original draft (supporting). G.B.: Methodology (equal); Software (equal); Writing – original draft (supporting). J.G.: Conceptualisation (equal); Funding acquisition (equal); Methodology (equal); Project administration (lead); Supervision (lead); Validation (supporting).
Appendix A. DFT with constraints
In this work, MD simulations use the canonical (
$N,V,T$
) ensemble, whereas DFT is formulated in terms of the grand canonical (
$\mu , V,T$
) ensemble. In DFT, the number of molecules
$N$
is not constant and different (but correct) solutions to the Euler–Lagrange equation can be obtained. Thus, algorithms keeping the number of molecules constant were proposed. Following Rehner & Gross (Reference Rehner and Gross2018), using the Lagrange multiplier
$\lambda$
, the unconstrained minimisation can be written as
\begin{equation} \mathcal{L}([\rho (\boldsymbol{r})], \lambda )=F+\int \rho (\boldsymbol{r}) V^{{\textit{ext}}}(\boldsymbol{r})\, \mathrm{d}{r} + \lambda \left (N-\int \rho (\boldsymbol{r})\, \mathrm{d} \boldsymbol{r}\right ) \stackrel {!}{=} \mathrm{min}. \end{equation}
After rearrangement, this provides an additional equation to the Euler–Lagrange equation. The density profiles can be determined from
$$ \rho =z e^{-\beta \left (\frac {\delta F_{{\textit{res}_{\!}}}}{\delta \rho }+V^{{\textit{ext}}}\right )}, $$
$$ z =\frac {N}{\int e^{-\beta \left (\frac {\delta F_{{\textit{res}_{\!}}}}{\delta \rho }+V^{{\textit{ext}}}\right )}\, \mathrm{dr}} ,$$
with the ensemble-averaged number of molecules
$N$
and the inverse thermodynamic temperature
$\beta ={1}/{k_{B} T}$
. Note that this approach can be used for mixtures where the ensemble-averaged number of molecules of each species is constant and with some adaption, the total ensemble-averaged number of molecules can be kept constant. We note that with these equations, the system is not canonical; it is a mathematical modification to obtain the solution with the desired ensemble-averaged number of molecules in a grand-canonical ensemble (Rehner & Gross Reference Rehner and Gross2018).
Appendix B. Estimating the external force
An external (body) force is employed in this work to induce the movement of droplets. In NEMD, the external force is added to the
$x$
-component of the force vector for each individual particle. In hydrodynamic DFT, the external force
$f_x$
is added to the gradient of the external potential in the momentum balance as an acceleration
$\rho f_x$
according to
\begin{equation} \frac {\partial (M\rho \boldsymbol{v})}{\partial t} + {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! \left (M\rho \boldsymbol{v} \boldsymbol{v}^\intercal \right ) = - \rho {\boldsymbol\nabla} \left ( \frac {\delta F}{\delta \rho } + V^{{\textit{ext},\textit{sf}}} \right ) + \rho f_x - {\boldsymbol\nabla} \!\boldsymbol{\cdot }\! \boldsymbol{\tau }, \end{equation}
where
$f_x$
has units of a force (per molecule). If divided by the molecular mass,
$f_x/M$
is equivalent to the earth’s gravitational acceleration
$g$
.
The movement of droplets is governed by the relation of the external driving force to capillary and viscous forces, cf. the momentum balance in (B1). Compared with macroscopic droplets, capillary forces have a much stronger influence on the droplet movement on the microscopic scale. The reason is that the external force is a body force, which acts on the whole volume of the droplet, while the capillary force acts mostly in the contact region. Since on the microscopic scale, the relation of the contact region to the droplet volume is much larger than at the macroscopic scale, much larger external (body) forces are required to obtain a similar motion of the droplet. The magnitude of the force required to overcome the capillary forces and to cause the droplet to move can be estimated by assuming that the capillary forces act on the contact region only and by neglecting viscous forces. Then, the influence of the capillary forces scales with
$c/V$
, where
$c$
is the length of the contact region (or contact line). For both spherical and cylindrical droplets,
$c/V \sim R^{-2}$
, where
$R$
is the radius of the respective droplet. Comparing a macroscopic droplet with
$R={1}\,{\mathrm{mm}}$
and gravitational acceleration of
$g\approx {10}\,{\mathrm{m\,s^{-2}}}$
, then a droplet with
$R={1}\,\mathrm{nm}$
needs to experience an acceleration of
$f_x/M\approx 10\times 10^{12}\,{\mathrm{m\,s^{-2}}}$
for an equivalent motion. This is roughly the same order of magnitude as the values used for methane used in this study, which are between
$2.1\times 10^{12}\,{\mathrm{m\,s^{-2}}}$
and
$8.4\times 10^{12}\,{\mathrm{m\,s^{-2}}}$
.
Appendix C. Adjustment of entropy scaling parameter
The parameter
$\psi$
scales the convolution radius for the weighted density
$\bar {\rho }^{{ES}}$
in (2.11) to quantitatively capture the influence of solid–fluid interactions on the viscosity close to the solid–fluid interface. As shown in previous work (Bursik et al. Reference Bursik, Eller and Gross2024b
), one liquid-phase reference NEMD simulation is sufficient to obtain an accurate and transferable parameter estimate. In this work, we employ a liquid-phase NEMD simulation in the same geometry used for contact angle simulations (cf. figure 1) at
$T={120.02}\,{\mathrm{K}}$
,
$f_x={0.112}\,{\mathrm{pN}}$
and
$\varepsilon_{\textit{sf}}^*=0.5$
. Since the velocity profiles are not strongly sensitive to the exact numerical value of
$\psi$
, we content ourselves with a rough estimate instead of a rigorous optimisation.
Figure 19 provides the velocity parallel to the solid–fluid interface
$v_x$
over the height of the system. Results are shown for hydrodynamic DFT with different values for
$\psi$
and for NEMD. We note that not all values tested are shown here for clarity. If
$\psi =3$
, hydrodynamic DFT underestimates the velocity profile compared with NEMD results, which suggests that the influence of solid–fluid interactions on the viscosity are overestimated. In contrast, using
$\psi =1$
, the agreement of the velocity profiles is considered to be sufficient and
$\psi =1$
is used throughout this work.
Appendix D. Mechanism for appearance of dynamic contact angles
In § 4.4, we argue that in hydrodynamic DFT, the mechanism behind the appearance of advancing and receding contact angles and their speed dependence is captured by the gradient of the external potential. Figure 20 shows the gradient of the (dimensionless) external potential in the
$x$
-direction. Due to the molecular roughness of the solid, the gradient can become negative which corresponds to a reduction of the momentum in the
$x$
-direction at this position. This is analogous to the potential wells and barriers described by molecular-kinetic theory (Blake & Haynes Reference Blake and Haynes1969; Blake Reference Blake, Wettability and Berg1993).
Velocity
$v_x$
profiles from hydrodynamic DFT using different values for the parameter
$\psi$
and from NEMD for
$f_x={0.112}\,{\mathrm{pN}}$
,
$T={120.02}\,{\mathrm{K}}$
and
$\varepsilon_{\textit{sf}}^*=0.5$
.
Gradient of the dimensionless external potential
${\partial ( \beta V^{{\textit{ext}}})}/{\partial x}$
in the
$x$
-direction as used in hydrodynamic DFT.
Appendix E. Quantifying the degree of rolling motion
In § 4.3, it is determined that the droplets adhere to a rolling motion. This was supported by the relative velocity in the
$y$
-direction, which is larger than
$1/3$
of the maximum absolute velocity. To further quantify the degree of rolling motion, figure 21 shows the length of the vorticity vector
$w$
of the two-dimensional flow field according to
\begin{equation} w = \left ( \frac {\partial v_y}{\partial x} - \frac {\partial v_x}{\partial y} \right ) ,\end{equation}
where negative values correspond to a clockwise rotation. The largest negative values are obtained in the advancing and receding part of the droplet close to the vapour–liquid interface and are approximately
${-0.04}\,{\mathrm{ps}^{-1}}$
. This would correspond to a full rotation in approximately 25 ps; since the droplet for the presented case requires approximately 300 ps to travel through the entire system (200 Å), the rolling motion is significant.
Length of vorticity vector
$w = ( {\partial v_y}/{\partial x} - {\partial v_x}/{\partial y} )$
of the droplet determined from hydrodynamic DFT.
Appendix F. Temperature in NEMD simulations
In the NEMD simulations, a global Nose–Hoover thermostat is applied to keep the temperature constant at
$T={120.02}\,{\mathrm{K}}$
. In a non-equilibrium simulation, shear effects can lead to local temperatures which deviate from the desired average temperature. Figure 22 shows the temperature from NEMD for the medium external force (
$f_x={0.112}\,{\mathrm{pN}}$
). The temperature within the droplet is close to the desired temperature (
$T={120.02}\,{\mathrm{K}}$
). The temperature is determined using relative velocities of molecules, which are obtained by subtracting the mean velocity in each bin from the absolute velocity of the molecules. In the advancing contact region, the temperature is slightly above 130 K, whereas in the receding contact region, the temperature is slightly below 110 K. Thus, local temperature variations are observed using the global Nose–Hoover thermostat, within approximately 10 % from the desired value.
Temperature
$T$
of the droplet for
$f_x={0.112}\,{\mathrm{pN}}$
determined from NEMD simulations and averaged over more than 10 000 ps.
Open access
