Liquid slip over gas nanofilms

Copyright and reuse: The Warwick Research Archive Portal (WRAP) makes this work by researchers of the University of Warwick available open access under the following conditions. Copyright © and all moral rights to the version of the paper presented here belong to the individual author(s) and/or other copyright owners. To the extent reasonable and practicable the material made available in WRAP has been checked for eligibility before being made available.


I. INTRODUCTION
A fundamental understanding of the flow of fluids over solid surfaces is important in many emerging technologies. Examples include self-cleaning surfaces, de-icing coatings for aerospace control surfaces, and drag-reducing or antifouling marine coatings. The classical no-slip fluid boundary condition is not necessarily applicable at micro-and nanoscales and the velocity slip at a surface is required in order to predict the overall flow behavior. The Navier boundary condition can describe the finite slip velocity at a planar wall and is given by V wall = b | wall , where | wall = dV /dy is the gradient of the velocity profile normal to the wall (in the y direction) and b is the slip length (which is the normal distance extrapolated into the wall at which the velocity would be zero). The slip length is a phenomenological parameter that determines the characteristics of the surface-fluid friction.
Several experimental [1,2] and theoretical studies [3,4] have suggested that the presence of surface nanobubbles and gas nanofilms [see Fig. 1(a)] may be responsible for large slip at liquid-solid interfaces. Vinogradova [5] has proposed the gas-cushion model (GCM) to calculate the apparent slip length b at a liquid-solid interface in the presence of a lubricating gas film [see Fig. 1(b)]. The slip length in this model is given by where the viscosity contrast C μ = μ l /μ g is the ratio of the liquid to the gas viscosities and Y g is the gas film thickness. This simple model follows from the assumption of a one-dimensional (1D) steady parallel flow, and continuity of stress and velocity (no slip) at both the solid-gas and gas-liquid interfaces. The GCM also assumes that there is a nonzero mass flow rate in the gas film, which may Schematic of (a) a surface nanobubble at the liquid-solid interface, (b) a zoomed in view showing the liquid-gas-solid section, and (c) the proposed velocity profile in the liquid-gas-solid section, with gas slip velocities indicated at both the solid-gas and gas-liquid interfaces. Note that Y g represents the local nanobubble thickness, which is modeled in (b) as a constant thickness film of height Y g . The effect of recirculation inside the bubble will be discussed in Sec. II. not be the case. A generalization of the GCM to a superhydrophobic surface was therefore proposed by Nizkaya et al. [6] to account for the zero mass flow rate in the entrapped gas phase, which may have an impact on the overall slip [7]. Theoretical predictions using Eq. (1) seem to agree well with experiments when the characteristic length scale of the gas film (i.e., Y g ) is much larger than the mean free path of the lubricating gas molecules. For instance, Nizkaya et al. [8] used atomic force microscopy to confirm the GCM theoretical prediction by measuring the local slip length in the gas regions. In their experiments the gas film thickness Y g = 1.9 μm is much larger than the gas mean free path λ (e.g., at standard atmospheric conditions, air molecules have λ = 68 nm [9]), which is given by where R s is the specific gas constant and T and ρ g are the temperature and mass density of the gas, respectively. The degree of gas rarefaction can be expressed by the Knudsen number Kn = λ/Y g , which in their experiment is 0.04, indicating that the gas is in the early slip regime. The GCM is commonly used to provide slip conditions at liquid-vapor interfaces to calculate the effective slip length of flows over surfaces consisting of alternating hydrodynamic slip boundary conditions [10,11]. We refer the reader to other recent articles [12,13] for various definitions of slip length. In the limiting case where Y g λ, de Gennes [14] estimated the slip length by equating the liquid viscous stress to the thermal stress in a gas such that where v x = √ k B T /2πm is the tangential velocity component of a gas molecule colliding with an interface, k B is the Boltzmann constant, and m is the mass of a gas molecule. 084003-2 By using Eq. (2) and expressing v x in terms of R s , Eq. (3) can be recast in a form similar to the GCM equation (1), i.e., The major limitation of the de Gennes model is that it is only applicable in the free-molecular regime (i.e., Kn > 10), so it does not provide a full description of the slip length in all Knudsen regimes [15]: the de Gennes model does not retrieve the predictions of the GCM at lower Kn [15]. While slip length predictions from the GCM are perfectly valid in the continuum regime (Kn < 0.001) and perhaps in the early slip regime (0.001 < Kn < 0.1), it is still not clear how the slip length varies when the gas is in the transition regime (0.1 < Kn < 10). In this paper we propose a modification to the GCM to include thermodynamic nonequilibrium gas effects in any gas film that interfaces with the liquid-solid regions in a shear flow. We term this model the rarefied-gas-cushion model (r-GCM). The r-GCM adopts the same basic assumptions (i.e., 1D parallel flow and continuity of stress) as the GCM, but relaxes the condition that the velocity should also be continuous. Instead, velocity slips at the two interfaces are modeled by the classical Maxwell slip boundary condition. Figure 1(c) shows a schematic of the steady-state velocity profile in one dimension throughout a liquid-gas-solid system with infinite dimension in the direction parallel to the surface, assuming that velocity-slip effects (due to rarefaction) occur in the gas layer. We assume that there is slip at both the solid-gas and gas-liquid interfaces and that there are two different strain rates in the liquid and the gas due to their contrasting viscosities. At steady state, the shear stress is constant and equal in the gas and liquid layers. We adopt a simple velocity slip model, so this means that our GCM is strictly only valid for small Kn (<0.1); however, in Sec. III we investigate its performance beyond this regime to test its qualitative predictive capability.

II. RAREFIED GAS-CUSHION MODEL
We start by defining the boundary conditions for the gas layer [as illustrated in Fig. 1(c)]. The slip velocity magnitudes at the solid-gas (sg) and gas-liquid (gl) interfaces can be expressed using Maxwell boundary conditions [16], assuming an isothermal system: where σ sg and σ gl are the tangential momentum accommodation coefficients and τ xy is the shear stress (which is constant in the y direction). For simplicity, we take σ sg ≈ 1 and σ gl ≈ 1 as experimentally typical values [17]. As illustrated in Fig. 1(c), the apparent slip length can be extrapolated from the strain rate in the liquid (the dotted line in the figure), i.e., where V l is the velocity of the liquid at the gas-liquid interface: By substituting Eqs. (5), (6), and (8) into Eq. (7), we obtain the apparent slip length prediction for the r-GCM: Inspection of this equation reveals that in the continuum limit (Kn → 0) the GCM is retrieved and in the free molecular limit (Kn 1) the de Gennes equation is retrieved.
084003-3 TABLE I. NEMD simulation data for the four cases (A-D) with different gas Knudsen numbers as shown in the second column. The gas layer thicknesses and bulk densities are shown in the third and fourth columns, respectively. The fifth and sixth columns show the gas bulk viscosities extracted from our NEMD simulations and the NIST database [19], respectively. The seventh and eighth columns show the gas bulk pressures from our NEMD simulations and the NIST database [19], respectively. The rest of this paper focuses on using molecular simulations to validate the theoretical slip length obtained using Eq. (9) for infinite gas nanofilms. However, for completeness we include here a note on extending our model to account for practical nanofilms, which are of finite size (e.g., a surface nanobubble, or entrapped gas in a nanotextured surface). For these cases, the slip length prediction in Eq. (9) needs to be recast to account for a backflow in the gas from consideration that there is no leakage of gas from the sides, which results in a net zero mass flux across the film, as suggested by Nizkaya et al. [18]. With this assumption, the derivation in the Appendix shows that Eq. (9) should be modified to In the continuum limit (Kn → 0), Eq. (10) retrieves the 1 4 prefactor in the local slip length, as found by Nizkaya et al. [18], which means that the slip is 4 times smaller in the finite nanofilm case. Equation (10) shows that for increasing Knudsen numbers, the reduction in local slip at gas-liquid interfaces increases to around 2 times in the free-molecular regime. Verification of this relation using molecular dynamics simulations of flows over nanobubbles would require huge computational resources and is beyond the scope of the present work.

III. MOLECULAR SIMULATIONS
To test our r-GCM [Eq. (9)], we perform four nonequilibrium molecular dynamics (NEMD) simulations (cases A-D) with varying Knudsen numbers ranging across the slip and transition regimes, as listed in Table I. We choose combinations of gas densities ρ g and film thicknesses Y g to give us a well-spaced Knudsen number range. All the simulations have a setup similar to the case shown in Fig. 2(a). This consists of a solid molecular bounding wall of fixed thickness ∼0.7 nm, a gas film of thickness Y g , and a liquid layer of fixed thickness ∼10 nm. The size of the computational domain in the x and z directions is L x = L z = 4.62 nm. The liquid region is decomposed into three zones: near-interface, bulk-measurement, and forcing zones. For computational efficiency we model the system as periodic in the x and z directions but nonperiodic in the y direction. The top boundary in contact with the liquid is a zero-shear-stress specular wall: molecules colliding with this boundary conserve tangential velocity, while their normal components are reversed, i.e., perfect slip. The system slip reference plane is taken in all cases as the first layer of atoms of the solid substrate, i.e., at y = 0.
Water is chosen as the contact liquid because of its relevance to marine engineering and recent advances in superhydrophobic drag-reducing surfaces. We take nitrogen as the gas, because it is the least soluble in water of all the species in air [20]. Platinum is chosen as the solid bounding wall in this work because its wetting characteristics are known both experimentally [21] and through molecular dynamics simulations [22,23]. Newton's equations of motion for all molecules in the system are integrated using the velocity-Verlet scheme, with a time step of 2 fs. For the molecular interactions, we use the (12-6) Lennard-Jones (LJ) potential and the smoothly truncated Coulomb potential for charge sites, i.e., where U (r ij ) is the total pair potential energy, r ij is the distance between atoms i and j , and σ are the LJ energy and distance parameters, respectively, r cut = 1 nm is the cutoff radius used in all our simulations for both the LJ and the Coulomb force calculations, 0 is the dielectric constant, and q i and q j are the atomic charges. Interactions between atomic sites of different molecules are derived from the Lorentz-Berthelot mixing rules [24]. In all our simulations, we use the interatomic potential parameters listed in Table II. Water is modeled using the rigid TIP4P/2005 force field [25], which consists of four sites: an oxygen atom (O), two positively charged hydrogen sites (H), and a negatively charged massless site (M). The H-O-H bond angle and the O-H bond length in a water molecule is fixed at 104.53 • and 0.095 nm, respectively, using Hamilton's quaternions. A three-site model is used for nitrogen gas [26], with two nitrogen atoms (N) and a positively charged massless site (M). The N-N bond length is fixed at 084003- 5  TABLE II. Lennard-Jones and Coulomb potential parameters for water [25], nitrogen [26], and platinum [27] molecules in our NEMD simulations.

Molecule
Atom 1.098 nm. The platinum substrate is modeled using FCC rigid atoms, but we list the LJ interactions for platinum from [27] in Table II for completeness, as these are needed for the Lorentz-Berthelot mixing rules. The steps of the simulation procedure are as follows. First, the lattices of the solid, gas, and liquid regions are initialized and allowed to equilibrate at a temperature of T = 298 K using a Berendsen thermostat. High gas adsorption occurs at the platinum surface during the equilibration step, which makes it very challenging to know beforehand how many gas molecules are needed in the initial gas lattice in order to produce the target gas density in the bulk zone of the gas film. To achieve the target gas and liquid bulk densities, we introduce a second simulation step to tune the densities using the FADE insertion-deletion protocol [28], while applying a simple velocity constraint procedure [29] on the water region to prevent the gas film departing from its target thickness.
The viscosity of water for the TIP4P/2005 model [30] at the set density ρ W = 1000 kg/m 3 is μ l = 0.855 × 10 −3 Pa s, which is needed for the slip length calculations. The viscosities of nitrogen gas μ g at 298 K are calculated from our NEMD simulations; our values are found to be in good agreement with the National Institute of Standards and Technology (NIST) database [19], as shown in Table I. Figures 2(a) and 2(b) show an example of one of the test cases (case A) at the end of the second simulation step. There is clearly a gas adsorption layer at the platinum surface. The bulk gas density and bulk gas layer thickness measured at the end of the equilibrium step are close to the prescribed values shown in Table I. Note that our measurements of the bulk density (for both liquid and gas) and viscosity are in regions away from the interfaces and do not contain any adsorption layer effects. At steady state and by Newton's third law, the pressure in the gas balances the pressure in the water. Using the TIP4P/2005 water model with a cutoff of 1 nm, in conjunction with nonperiodic boundary conditions in the y direction, produces <0.5% error in the water density and <20% error in the measured pressure between our NEMD simulations and experiments [19], which could be improved at a computational cost by increasing the intermolecular potential cutoff distance.
After the system is equilibrated, the flow simulation proceeds by imposing a constant shear stress of τ xy = 0.3 MPa to the liquid molecules in the forcing zone via a uniform force F given by where Y = 3.06 nm is the forcing zone thickness and ρ n is the liquid number density in the forcing zone. This value for τ xy is chosen to be large enough to produce a reasonable signal-to-noise ratio in the particle ensemble, but well below the magnitude at which the strain varies nonlinearly with stress and also below the critical stress at which the slip length diverges. A velocity-unbiased Berendsen thermostat is applied in bins of width 0.8 nm within the liquid region and one bin in the gas region, for the entire duration of the simulation. After the system reaches a steady state (after ∼60-100 ns), which is dependent on each problem, the simulation is run for a further period (∼10-20 ns) in order to measure the gas bulk densities and pressures (see Table I) as well as the average velocity V in the bulk zone of the liquid region. The apparent 084003-6 TABLE III. Apparent slip length results for the four cases. The second, third, and fourth columns show the slip lengths calculated from the NEMD simulations, the standard GCM equation (1), and the r-GCM equation (9). The values in parentheses in the second column indicate upper and lower 95% confidence limits. The last two columns indicate the absolute percent errors with respect to our NEMD results.
where μ l is the known liquid viscosity in the bulk zone and H is the known height from the reference plane (y = 0) to the midpoint of the bulk measurement zone, as shown in Fig. 2(a). All our simulations were performed using the mdFoam+ solver [29,31], a highly parallelised molecular dynamics code that is implemented within the open-source OpenFOAM C++ libraries. These simulations were run on ARCHER, the UK's national supercomputer.

IV. RESULTS AND DISCUSSION
The slip lengths predicted by our four NEMD simulations (cases A-D, with increasing gas Kn) at steady state are compared with calculations from the GCM [Eq. (1)] and our proposed r-GCM [Eq. (9)] in Table III, in the second, third, and fourth columns, respectively. The values in parentheses in the second column indicate upper and lower 95% confidence limits due to the standard error in V . The last two columns are the absolute percentage error of the theoretical slip lengths with respect to the NEMD slip length, which are calculated using where b NEMD is the apparent slip length calculated from our NEMD simulations and b theor is the theoretical apparent slip length prediction from either the r-GCM or the GCM. We observe that the error in the slip length predictions of the GCM increases with Kn in the transition regime (∼10%-50% for these four cases), while our r-GCM error remains roughly constant at 4%-8%. Figure 3 synthesizes the results presented in Table III: The normalized apparent slip lengths b/Y g are plotted against Kn for the four cases (A-D) using molecular dynamics (closed circles), the r-GCM (solid line), the GCM (dotted line), and the de Gennes model (dot-dashed line). Note that for the three theoretical calculations we use the average gas viscosity across all our cases (see the fifth column in Table I); this approximation is reasonable as we have verified that the slip lengths do not change significantly when individual gas viscosity values are used. As can be seen in Fig. 3, the r-GCM is in good agreement with our NEMD results in the transition regime and also converges to the GCM slip length in the slip and continuum limit. These molecular dynamics simulation results indicate the usefulness of our simple rarefied-gas-cushion model, as well as the important role that rarefaction (through Kn) plays in determining the slip properties of interfacial and multiphase problems at the micro-and nanoscales. Our r-GCM also agrees very well with the de Gennes model towards the free-molecular regime, with a 4% difference between the two models at Kn = 10 that decreases to only 0.5% at Kn = 100.  (1) and our proposed r-GCM equation (9). The theoretical slip lengths are calculated using a gas dynamic viscosity of 1.82 × 10 −3 Pa s, which is the average of the gas viscosities of the four cases, at 298 K. Note that error bars in the NEMD results are much smaller than the size of the symbols and so have been omitted from this graph.
All our NEMD simulations in this paper are computationally expensive, requiring around 60-100 ns of problem time and 30-60 days of computational time on 48 processors. For Knudsen numbers below the transition regime (i.e., in the continuum regime and early slip regime) simulations are much more computationally costly, since a decrease in Knudsen number requires the gas film thickness to be increased, which increases the evolution time scales. For instance, if Kn = 0.02 (corresponding to the early slip regime) and the gas thickness is 100 nm (ρ g = 40 kg/m 3 ), the problem time to reach steady state would be ∼600 ns. A similar difficulty is encountered at higher Knudsen numbers (in the upper transition regime); at Kn = 4 and for a gas thickness of around 15 nm (ρ g = 1.2 kg/m 3 ), there are only about ten gas molecules in the gas film (which opens up questions on the stability or existence of the gas layer), and in any case the problem time to reach steady state can be estimated to be around 1 μs. This has prevented us from using NEMD simulations to investigate how the r-GCM performs in the upper range of the transition regime and into the free-molecular regime.
To complete our analysis, we briefly examine the cause of the small errors (∼5%) seen between our r-GCM and the NEMD results. Figure 4 shows the steady-state velocity profiles along the y direction extracted from our NEMD simulations, with the expected steady-state velocity profiles from the GCM and the r-GCM. As expected, at low Knudsen numbers (i.e., in the slip regime) our NEMD results are in good agreement with the r-GCM, although the GCM solution is still not a bad prediction. This demonstrates why the GCM might still be a good estimator of slip length in the early parts of the slip regime. While the r-GCM results seem to agree better with our NEMD results in the transition regime than the GCM results do, we observe that the deviations are due to the nonlinear strain rate within the gas and possibly also the assumption of σ sg = σ gl ≈ 1. The nonlinear velocity profile is due to the well-known Knudsen layer phenomenon; in order to enable higher accuracy in these theoretical gas slip length models one would need in the future to incorporate more sophisticated results from kinetic theory into the model.
We now use the r-GCM equation (9) and the gas mean free path equation (2) to make predictions for a range of typical gas layer thicknesses and gas densities, without resorting to additional expensive NEMD simulations. From the discussion above, we expect around 5% error in our calculations. The literature concerning surface nanobubbles is still in its infancy [32] and there is uncertainty about 084003-8 the actual pressures and densities of the confined gas. We choose three sets of density ranges for gas films at the nanoscale (Y g < 100 nm), which may cover the operating conditions of these gas bubbles or thin films: the atmospheric density range (ρ g = 0.5-3 kg/m 3 ), a medium density range (ρ g = 3-30 kg/m 3 ), and a much-higher-density range (ρ g = 30-300 kg/m 3 ) as observed by Peng et al. [33] and very recently by Che and Theodorakis [34] in their molecular dynamics studies of nitrogen gas adsorption at a water-graphite interface. Contour plots of Knudsen number and apparent slip lengths are in Fig. 5 for these three sets of densities. Our first observation is that across these gas densities the operating conditions have rarefied-gas features [Figs. 5(a)-5(c)]. In the low-density range, the flow is predominantly in the transition regime [ Fig. 5(a)], which means that the r-GCM [ Fig. 5(d)] predicts much larger slip lengths than the standard GCM [ Fig. 5(g)]. For example, a gas film with a height 100 nm and ρ g = 1.5 kg/m 3 has a Knudsen number of ∼0.5 and produces an apparent slip length of 9.4 μm, which is twice the GCM prediction. As the gas density increases, the Knudsen number decreases; in the medium-density range, the transition regime still covers 50% of the operating conditions, so the r-GCM is still required to make accurate predictions. In higher-density conditions, the gas is less rarefied and the slip Knudsen regime is now dominant, meaning that the GCM and r-GCM results are very similar [compare a continuous thin film, such as the case of gas entrapped within a series of corrugated ridges (Fig. 6, inset), then the effective slip length b eff of the surface is a combination of the slip lengths of the solid-liquid (b sl ) and the gas-liquid (b gl ) interfaces, approximated by for b sl > 0 and b sl ,b gl > L/10, where ζ is the gas-surface coverage and L is the distance between the solid ridges [10,11]. Several experiments on smooth hydrophobic surfaces have shown that the slip length b sl does not exceed 100 nm [36][37][38][39], which means that bubble coverage plays a substantial role in the overall effective slip length of the surface, according to the form of Eq. (15). Only surfaces with bubble coverage greater than 95% (ζ = 0.95) produce effective slip lengths that are as large as those seen locally at gas-liquid interfaces; at lower bubble coverage, b eff rapidly approaches b sl , making the presence of bubbles ineffective at reducing drag. Our results in Fig. 6 however show that the effective slip length calculated using the GCM deviates from that using the r-GCM when incorporated in the mixed hydrodynamic model of Eq. (15) and this error depends on the Knudsen number. Whether entrapped gas with a very large surface coverage or thin films completely covering a surface can be generated and be stable at the nanoscale with low gas density remains an open problem and would need further investigation.

V. CONCLUSION
From fundamental results in gas kinetic theory, we have derived a theoretical model to estimate apparent slip lengths of shear-driven flows over gas nanofilms. For these nanoscale flow problems, the standard gas-cushion model and the de Gennes model are each applicable only at extreme ends of the Knudsen number scale. We used nonequilibrium molecular dynamics simulations to investigate water flows over gas nanofilms of varying thickness and Knudsen number and these results compare favorably with calculations using our rarefied-gas-cushion model, as well as demonstrating that the GCM is not accurate outside the gas continuum or slip regimes. We showed that the r-GCM is surprisingly effective over a wide range of Knudsen numbers, although validation in the upper part of the transition regime and into the free-molecular regime is computationally intractable. In such cases and to verify the flow behavior for lower gas densities and larger gas film thicknesses, we suggest performing hybrid NEMD-kinetic simulations of these nanofilms or nanobubbles (e.g., using the direct simulation Monte Carlo technique for the bulk gas region, with NEMD at the interfaces) to provide further testing of our theoretical model. The r-GCM and the contour plots presented in Sec. IV can support any future experimental investigations in this research area.  Furthermore, the stress and strain at the gas-liquid interface are zero for the pressure-driven component, so it follows that the shear in the liquid in the Couette component is equal to the shear in the liquid in the total flow, i.e., where b * [indicated in Fig. 7(c)] is the apparent slip length with the zero-mass-flow-rate constraint and b is the apparent slip length derived for infinite films [see Eq. (9)]; both can contain rarefied-gas effects. Rearranging Eq. (A4) and substituting b from Eq. (9) gives b * = C μ (2 Kn + 1) 1 + 6 Kn 4 + 12 Kn − 1 Y g . (A5) In the continuum limit, i.e., as Kn → 0 and assuming that C μ 1, Eq. (A5) retrieves the slip length over entrapped gas in cavities [18], i.e., which demonstrates that the slip is a quarter of that predicted by the GCM, i.e., b/b * ≈ 4. For the rarefied gas films presented in this work, the local slip ratio can be written as b b * ≈ 4 + 12 Kn 1 + 6 Kn .
If our slip model can be applied at high Kn, this indicates that for free molecular gas flow the slip for finite nanofilms is half that of an infinite film, as demonstrated in Fig. 8.