Shear force measurement of the hydrodynamic wall position in molecular dynamics

Walls impose boundary conditions (BCs) to the flow of liquids. The commonly used no-slip BC fails to describe flows in nanofluidic systems¹ and needs to be replaced by the more general partial slip BC.² This BC can be expressed in terms of stress, as initially done by Navier:^3,4 the viscous shear stress in the liquid at the wall, $η∂zv|z=zs$ (with η being the shear viscosity, z being the normal to the wall, z_s being the wall position, and $v$ being the velocity parallel to the wall) is equal to an interfacial friction stress $τw$ proportional to the slip velocity $vs$⁠, i.e., the jump of parallel velocity at the interface,

defining the (fluid) friction coefficient λ. The partial slip BC can also be written as a kinematic relation on the velocity field at the interface,

where the slip length b is uniquely related to the friction coefficient λ for a fluid with a given viscosity: b = η/λ.

Regardless of its form, the partial slip BC involves two independent parameters: the slip length b (or equivalently the friction coefficient λ) and the hydrodynamic wall position (HWP) z_s, where the BC applies. To accurately predict flows in nanofluidic systems, both parameters must be known. Molecular dynamics (MD) simulations can provide such information^5–23 and have been used to explore the molecular mechanisms of liquid-solid slip.^24–36 A liquid confined between parallel flat walls, with periodic BCs in the lateral directions, is commonly used. To measure both b and z_s, one can simulate two types of flow in the same system, typically a Couette and a Poiseuille flow, or simulate a Poiseuille flow with two system heights; see Ref. 37 and references therein. However, these measurements are usually delicate, and the measured b and z_s are very sensitive to the fits of the flow profile, which are affected by thermal fluctuations. Alternatively, Bocquet and Barrat⁵ derived Green-Kubo formulas for the friction coefficient and for the position of the wall; however, applying these formulas in finite-size MD simulations is delicate.^9,15–17

Few studies have attempted to measure the hydrodynamic position of the wall. For generic Lennard-Jones (LJ) fluids and walls of different wettabilities and corrugation, the reported shift Δ between the wall surface and the HWP varied between ∼1.1 and 2.5σ, with σ being the molecular diameter.^5–7,17 In practice, most studies assume a given position for the hydrodynamic wall, ranging typically from the physical wall position up to the position of the liquid first absorption layer^{2,11,18,21,30,38–40} or at the Gibbs dividing plane (GDP),^10,12,32,41 which leaves the simpler task to determine only one parameter, slip length or friction coefficient.

In this article, we will show that the HWP can be accurately extracted from a single Poiseuille flow simulation by measuring the shear stress on the wall. We will also show that the common practice of applying a gravity-like force per particle to generate a Poiseuille flow in MD simulations imposes by construction that the HWP identifies with the GDP. Finally, we will compare the HWP obtained for a pressure-driven flow generated by a fluid piston and the GDP measured at equilibrium.

We will start by showing the relation between interfacial shear stress and HWP in a Poiseuille flow. To that aim, we would like to emphasize that the partial slip BC effectively takes into account all the phenomena occurring in the molecular vicinity of the interface and provides a BC for the flow far from the interface, where the liquid is described by its bulk properties. As a consequence, in the partial slip BC, the velocity and the shear rate at the interface need to be obtained from the extrapolated bulk velocity field (i.e., using the bulk liquid viscosity), regardless of the true velocity field at the interface, where viscosity may locally change.^42,43 This is the rule we will follow in the derivation below.

We consider a Poiseuille flow induced by a constant force density f in a liquid confined between two parallel walls located at a vertical position z = ±H/2; see Fig. 1. The force density can be due to a pressure gradient, f = (−∇p), or to a gravity-like field g, f = ρg, with ρ being the bulk liquid density. The walls impose a partial slip BC (with a slip length b) applying at a distance Δ from the physical walls (defining the hydrodynamic height h = H − 2Δ). The velocity profile is obtained by integrating the Stokes equation, $∂z2v=−f/η$⁠, using the symmetry of the profile at z = 0, and the partial slip BC at the bottom wall, $v$(−h/2) = b ∂_z$v$|_z=−h/2,

The shear stress applied by the liquid to the wall is $τw$ = η ∂_z$v$|_z=−h/2 = f h/2. For a given force density f, the hydrodynamic height can therefore be measured via the interfacial shear stress,

Using this relation, we now would like to discuss a common approach used in MD simulations to impose a Poiseuille flow, hereafter referred to as “gravity-like flow,” where one applies a force per particle f_i = f/n_bulk (with n_bulk being the bulk number density) to liquid particles. In that case, the total force applied to the liquid, F = Nf_i = Nf/n_bulk, with N being the number of liquid particles, is equal to the total shear force between the liquid and the two confining walls, F = 2$Sτw$ = Sfh, with S being the wall area. From this equality, it results that

Introducing the number density profile n(z), Eq. (5) can be rewritten as

i.e., the HWP identifies with the GDP, corresponding to a partitioning of space between a region filled with a homogeneous liquid and another one without any liquid (see Fig. 3). Therefore, the choice made in previous work^10,12,32,41 to fix the HWP at the GDP can be justified based on hydrodynamics arguments. Additionally, because parallel flows do not affect significantly density profiles perpendicular to the walls,⁴⁴ Eq. (5) indicates that the HWP for a gravity-like flow can in fact be measured from the GDP in equilibrium simulations as we will do in the following.

To test the theoretical predictions presented above, we performed MD simulations of a liquid confined between two parallel walls (see Fig. 1) using the LAMMPS package.⁴⁵ Liquid-liquid and liquid-solid interactions were modeled with a Lennard-Jones (LJ) potential, $V(r)=4εij[(σij/r)12−(σij/r)6]$⁠, with r being the interparticle distance, ε_ij and σ_ij being the interaction energy and size, where i and j can be L for liquid particles or S for solid ones. In the following, we will use reduced units based on particle mass m, and liquid-liquid interaction energy ε = ε_LL and size σ = σ_LL. In particular, the unit of time is $τ=σm/ε$⁠. The liquid-solid interaction energy ε_LS was varied between 0.3 and 0.6ε, while keeping σ_LS = σ. The potential was truncated at 2.5σ. For the walls, we used three atomic layers of a frozen face centered cubic crystal exhibiting a (001) face to the liquid, with an interparticle distance corresponding to mechanical equilibrium, d = 2^1/6σ. We used periodic BCs along the lateral x and y directions, and all the measurements were taken in a region with a lateral size L ≡ L_x = L_y ≈ 19σ with 5206 liquid particles and 1728 solid ones. The temperature was set to T = 0.83ε/k_B by applying a Nosé-Hoover thermostat to liquid particles, only along the y and z degrees of freedom perpendicular to the flow, and with a damping time of 0.5τ. Equivalent results were obtained for different damping times and using a Berendsen thermostat. The pressure was set to 0.094ε/σ³ by using the top wall as a piston during an equilibration stage and fixing it at its equilibrium position during production. The resulting system physical height H, defined as the distance between the first inmost layers of the walls, varied between 21 and 22σ for different ε_LS values. The equations of motion were integrated using the velocity-Verlet algorithm, with a time step of 0.005τ.

To generate a pressure-driven flow, we used a fluid piston^44,46,47 [see Fig. 1(c)]: we increased the box size along the x direction, from L_x ≈ 19σ to L_x ≈ 43σ (using 11 713 liquid particles and 3888 solid ones), we applied along the x direction a force per particle to liquid particles in a thin slab of length $lxpiston≈8σ$ (the fluid piston region), and we measured the wall shear stress τ_w and the bulk pressure gradient f = (−∇p) in a measurement region of length $lxmeasure≈19σ$ far from the fluid piston. We did not observe any significant difference in the results for a bigger region between the fluid piston and the measurements region. We computed the hydrodynamic height h using Eq. (4), and the corresponding hydrodynamic shift Δ = (H − h)/2. From the fit of the Poiseuille flow profile in the bulk region with Eq. (3), we could also extract the slip length b and the viscosity η. We also computed the distance between the GDP and the physical wall, denoted Δ′, using Eq. (5) in equilibrium simulations (with no external force applied to the system). Note that we also measured the position of the GDP in the fluid piston simulations, and the results were identical to the equilibrium ones.

Finally, to compare the fluid piston results with another more common approach, we performed independent Couette flow simulations on the system illustrated in Fig. 1(a), with lateral size L ≡ L_x = L_y ≈ 19σ and 5206 liquid particles. We used the hydrodynamic height h measured in the fluid piston simulations; we measured the viscosity as the ratio between the wall shear stress and the bulk shear rate, and the slip length from a fit of the bulk velocity profile. In all the simulations, the equilibration stage lasted 2 · 10⁵ time steps, and the production lasted 10⁷ time steps.

Both for Poiseuille and Couette flows, we simulated a number of different forcing (pressure gradient, f ∈ [2.7; 4.7] · 10⁻³ε/σ⁴, or shear velocity, U ∈ [0.1; 0.5]σ/τ). Five independent simulations were run for each value of the force density. All the results shown in this paper (Δ, Δ′, b, and η) were obtained for a given ε_LS by averaging the results which belonged to the linear response regime, i.e., the forcing range in which the measured quantity remained constant. Correspondingly, the maximum shear rates produced were 0.033 and 0.048τ⁻¹ for the fluid piston and the Couette simulations, respectively. These shear rates are below the shear thinning regime of the LJ fluid, around 0.07τ⁻¹; see Ref. 48. The given error bars correspond to the statistical error within 95% of confidence level.

Figure 2 presents the measured shifts between the wall surface and the HWP, Δ from fluid piston simulations using Eq. (4) and Δ′ from equilibrium simulations using Eq. (5), i.e., using the position of the GDP. Note that, even though Δ′ was measured from equilibrium simulations, we will still refer to it as the gravity-like flow hydrodynamic shift because Eq. (5) was derived analytically for a gravity-driven flow. In this figure, one can see that Δ′ is comparable to Δ although they slightly differ for large ε_LS. In general, the hydrodynamic shifts for a pressure-driven and for a gravity-like flow do not have to be identical: indeed, while the pressure gradient and the gravity-like force identify in the bulk, they will generally differ in the molecular vicinity of the interface where the fluid becomes heterogeneous. In particular, the gravity-like force distribution will follow that of the density and can result in a different effective BC for the bulk flow. In Fig. 2, one can also see that, for both flows, the distance between the wall surface and the HWP decreases with the wall wettability controlled by the ε_LS parameter (a higher ε_LS corresponds to a more wettable system). We can understand this result in terms of GDP (Fig. 3): for a more wettable system, the peaks of the density profile close to the wall are more pronounced, so the area of an effective liquid with constant density is bigger, hence a larger value of h and a smaller value of Δ. We also varied ε_LL and found that its impact was much smaller than that of ε_LS. This can be rationalized based on the GDP. Indeed, changing ε_LL should (at the first order) rescale the whole liquid density profile so that the GDP should not change.

The results we obtained here for the HWP are significantly different from the ones obtained by other studies^5–7,17 (with Δ between ∼1.1 and 2.5σ) and from the sometimes used assumption^11,18 that Δ = 0σ. For the less wettable walls, the common assumption^2,30,38–40 that Δ ∼ 1σ agrees with our results. However, for higher ε_LS values, this assumption is generally not valid, especially in systems with significantly small slip lengths like the one discussed in this paper.

From the pressure-driven flow, in addition to the determination of the HWP, one can also measure the system transport coefficients. Table I reports the η and b results for a set of ε_LS parameters. Because the shear viscosity is a property of the bulk liquid and in all simulations the temperature and pressure are set constant, there is no effect of the wall wettability on η, and its value is comparable to the one obtained from Couette simulations. As rationalized in previous work,⁴⁹ the slip length decreases with ε_LS: for a less wettable system, the fluid friction coefficient is smaller which implies a higher value of b = η/λ. If we compare these results with those obtained from Couette simulations, we can see that by means of the pressure-driven flow method, we obtained equivalent b and η measurements with the same order of magnitude in the error precision.

We have shown that the position where the hydrodynamic BC imposed by walls should be applied can be efficiently determined by measuring the wall shear stress in MD simulations of a Poiseuille flow. As a consequence, we have shown that for gravity-driven flows, the HWP is only controlled by the static density profile of the fluid close to the wall and identifies with the GDP, which can be measured from equilibrium simulations. Accordingly, the HWP could be estimated from previous work where the equilibrium structure of liquid-solid interfaces was modeled.⁵⁰ We then simulated a LJ fluid confined between two parallel frozen walls and measured the HWP by using a fluid piston to generate a pressure-driven flow. The pressure-driven flow hydrodynamic wall was comparable (although not identical) to the GDP. We investigated the effect of wetting by varying the liquid-solid interaction energy. We found that the hydrodynamic BC applies in the liquid, at a distance Δ from the wall surface varying from ∼1σ (with σ the atomic diameter) on nonwetting walls to a fraction of σ on wetting walls. The decrease in Δ for increasing wetting can be rationalized in terms of GDP, which is shifted toward the solid when the adsorption of the fluid increases on more wetting surfaces. The measured values of Δ are generally lower than previous measures, which ranged between 1.1 and 2.5σ, but they correspond approximately to the standard assumption made in MD studies of liquid-solid slip that Δ ∼ 1σ. Finally, we have shown that, in addition to the HWP, the Poiseuille flow simulation also provides an accurate estimate of the slip length and fluid viscosity by comparing the measured values with those obtained from independent Couette flow simulations.

Overall, we have presented a simple, fast, and accurate method to fully characterize the transport properties of a confined fluid by measuring the viscosity, slip length, and HWP in a single Poiseuille flow simulation. Note that this method is not limited to the simple slab geometry considered here. For instance, it can easily be extended to cylindrical channels, where the wall shear stress is τ_w = fr/2, with f being the pressure gradient and r being the hydrodynamic radius of the pore. The method should also apply to mixtures—for a gravity-like flow, one can show that Eqs. (5) and (6) apply when replacing the number of liquid particles and the number density by the total mass of the liquid and the mass density, respectively—and to thermalized walls. We plan to investigate the latter case in the future. An analogous approach could also be applied to characterize the effective wall position for interfacial heat transfer. We hope these results will help improve the characterization of the hydrodynamic BC by MD simulations in systems of interest for nanofluidic applications.

This work was supported by the ANR, Project No. ANR-16-CE06-0004-01 NECtAR. T.O. and Y.Y. are supported by JSPS KAKENHI Grant Nos. JP18K03929 and JP18K03978, Japan, respectively. Y.Y. is also supported by JST CREST Grant No. JPMJCR18I1, Japan. L.J. is supported by the Institut Universitaire de France. L.J. benefited from a JSPS international fellowship for research in Japan.