{'fraction_non_alphanumeric': 0.03582591969223371, 'fraction_numerical': 0.018153402260158692, 'mean_word_length': 4.334081436357807, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Lattice Boltzmann simulations of apparent slip in hydrophobic microchannels

Various experiments have found a boundary slip in hydrophobic microchannel flows, but a consistent understanding of the results is still lacking. While Molecular Dynamics (MD) simulations cannot reach the low shear rates and large system sizes of the experiments, it is often impossible to resolve the needed details with macroscopic approaches. We model the interaction between hydrophobic channel walls and a fluid by means of a multi-phase lattice Boltzmann model. Our mesoscopic approach overcomes the limitations of MD simulations and can reach the small flow velocities of known experiments. We reproduce results from experiments at small Knudsen numbers and other simulations, namely an increase of slip with increasing liquid-solid interactions, the slip being independent of the flow velocity, and a decreasing slip with increasing bulk pressure. Within our model we develop a semi-analytic approximation of the dependence of the slip on the pressure.

During the last century it was widely assumed that the velocity of a Newtonian liquid at a surface is always identical to the velocity of the surface. However, in recent years well controlled experiments have shown a violation of the no-slip boundary condition in sub-micron sized geometries. Since then, experimental [1] and theoretical works [2], as well as computer simulations [3][4][5][6][7][8][9] have tried to improve our understanding of boundary slip. The complex behavior of a fluid close to a solid interface involves the interplay of many physical and chemical properties. These include the wettability of the solid, shear rate, pressure, surface charge, surface roughness, as well as impurities and dissolved gas. Since all those quantities have to be determined very precisely, it is not surprising that our understanding of the phenomenon is still unsatisfactory. Due to the large number of different parameters, a significant dispersion of the results can be observed for ostensibly similar systems [1], e.g. observed slip lengths vary between nanometres [10] and micrometers [11] and while some authors find a dependence of the slip on the flow velocity [12], others do not [11,13]. Most computer simulations apply Molecular Dynamics (MD) and report increasing slip with decreasing liquid density [6,7] or liquid-solid interactions [8,14], while slip decreases with increasing pressure [4]. These simulations are usually limited to some tens of thousands of particles, lengths scales of nanometres and timescales of nanoseconds. Also, shear rates are orders of magnitude higher than in any experiment [1]. We overcome these limitations using the lattice Boltzmann (LB) algorithmc EDP Sciences a powerful method for simulating fluid dynamics [15]. Rather than tracking individual atoms and molecules, the dynamics of the single-particle distribution function η of mesoscopic fluid packets is described. In contrast to MD simulations, this method is less computationally demanding and allows to simulate experimentally accessible length and time scales. Our ansatz differs from other LB approaches where slip is introduced by generalizing no-slip bounce back boundary conditions to allow specular reflections with a given probability [3] or where the viscosity is modified due to local density variations [16]. In both cases, parameters determining the properties at the boundaries are not easily mappable to experimentally available values. Our approach is based on Shan and Chen's multi-phase LB model [17]. Here, interactions between different species are modelled by mesoscopic forces between the phases. This naturally opens the way to introduce similar interactions between each fluid species and the channel walls, where the strength of the interaction is determined by the fluid densities, free coupling constants, and a wall interaction parameter which is treated in a similar manner as a local fluid density. The model allows the simulation of multi-phase flows along hydrophobic boundaries and is introduced in the following. However, in order to study the influence of hydrophobicity on the boundary slip and to demonstrate the basic properties of the model, we focus on single phase flow in this paper. Results of multi-phase simulations will be presented in a future work. A multi-phase LB system can be represented by a set of equations [18] where η α i (x, t) is the single-particle distribution function, indicating the amount of species α with velocity c i , at site x on a D-dimensional lattice of coordination number b (D3Q19 in our implementation), at time-step t. For the collision operator Ω α i we choose the Bhatnagar-Gross-Krook (BGK) form where τ α is the mean collision time for component α and determines the fluid viscosity. The system relaxes to an equilibrium distribution η α eq i which can be derived imposing restriction on the microscopic processes, such as explicit mass and momentum conservation for each species [19]. η α (x, t) ≡ i η α i (x, t) is the fluid density and u α (x, t) is the macroscopic velocity of the fluid, defined as η α (x, t)u α (x, t) ≡ i η α i (x, t)c i . Interactions between different fluid species are introduced as a mean field body force between nearest neighbors [17]: where ψ α (x, t) = (1 − e −η α (x,t)/η0 ) is the so-called effective mass with η 0 being a reference density that is set to 1 in our case [17]. gᾱ α is a force coupling constant, whose magnitude controls the strength of the interaction between component α andᾱ. The dynamical effect of the force is realized in the BGK collision operator in Eq. (2) by adding to the velocity u in the equilibrium distribution an increment δu α = τ α F αᾱ /η α . For the interaction of the fluid components with the channel walls we apply mid-grid bounce back boundary conditions [15] and assign interaction properties to the wall which are similar to those of an additional fluid species. I.e., we specify constant values for the force coupling constant gᾱ α = g wall,α and the density ηᾱ = η wall at wall boundary nodes of the lattice. This results in a purely local force as given in Eq. 3 between the flow and the boundaries. Even though one could argue that a single parameter to tune the fluid-wall interaction would be sufficient, we keep our approach as close as possible to the original idea of Shan and Chen in order to benefit from the experience obtained from other works using the original model. Furthermore, the additional parameter allows more flexibility when simulating not only a single fluid, but a multi-phase system. The fluid-wall interaction can be linked to a contact angle between fluid droplets and solid walls as it is often used to quantitatively describe hydrophobic interactions [20]. Recently, Benzi et al. have shown how to compute the contact angle within the Shan-Chen model [21]. The same authors also developed an approach to model apparent slip which is related to ours, but instead of using only local fluid-solid interactions, they add an exponential decay of the interaction with distance from the wall [22]. We simulate pressure driven flow between two infinite planes (Poiseuille flow), where pressure driven boundary conditions are implemented in a similar way as in most experiments: a fixed pressure is set at the channel inlet and an open boundary at the outlet. The outlet is realized by interpolating the particle distribution function at the end of the channel as given by leading to a linear pressure gradient. Already in 1823, Navier proposed a boundary condition where the fluid velocity at a surface is proportional to the shear rate at the surface, i.e. v z (x 0 ) = β∂v z (x)/∂x at x = x 0 [23]. Following his hypothesis the velocity in flow direction (v z ) at position x between the planes is given by where 2h is the distance between the planes, and µ the viscosity. In contrast to a no-slip formulation, the last term in Eq. 4 linearly depends on the slip length β. Since β is typically of the order of nanometers or micrometers, it can be neglected in macroscopic experiments. In order to obtain β from our data, we measure the pressure gradient ∂P /∂z at the center of the channel and the velocity profile between the two planes at a fixed position z. β is then obtained by a least square fit with Eq. 4. Our simulation parameters are as follows: the lattice size is kept constant with the channel length (z direction) being 256 sites, the distance between the plates 2h being 60 sites (x direction). We approximate infinite planes by using a 16 sites wide channel with periodic boundaries in y direction. In order to assure a fully equilibrated system we simulate for at least 40000 time steps before measuring and assured our results being independent of the discretization level by comparing to simulations of 28 and 124 sites wide channels. Each data point in the figures below corresponds to about six hours simulation time on eight IBM Power 4 1.7GHz CPUs. All units in this paper are in lattice units with the lattice constant c and timestep ∆t set to 1 if not stated otherwise.
The dependence of the slip length β on the interaction parameter g wall,α is studied for η wall =1.0 and 5.0. The bulk pressure P = ρc 2 s , where ρ is the fluid density and c s = 1/ √ 3 the speed of sound, is kept at P =0.11, while the flow velocity is set to V =0.033. As shown in Fig. 1a we vary g wall,α from 0.06 to 0.22 and find a steady increase of β for increasing g wall,α . As expected, the curve for η wall =5.0 is growing substantially faster than for η wall =1.0. The maximum available β are at about 5.2 for g wall,α =0.26 and η wall =1.0. At these strong fluidwall interactions, the force as given in Eq. 3 becomes very large and results in a large area of low fluid density close to the wall. Increasing the interaction even further results in numerical instabilities due to too steep density gradients. In order to study the dependence of the slip on other parameters, the coupling constant g wall,α is kept constant at 0.08 from now on. Fig.1b depicts the dependence of β on η wall for different bulk pressures P =0.033, 0.1, and 0.3 and fixed flow velocity V = 3.5 · 10 −3 in the system. While all three graphs grow constantly with increasing η wall , the one for P =0.033 grows the fastest demonstrating that absolute values for β are higher for lower pressure.
We have measured the magnitude of the boundary slip over a very wide range of flow velocities V from 1 · 10 −4 to 3 · 10 −2 for wall interactions η wall =0.0, 0.5, 1.0, and 2.0. V is measured at the center of the channel and given on a logarithmic scale in Fig. 2. For η wall =0.0 we do not find any boundary slip confirming that our method properly reproduces no slip behavior in the interaction free case. With increasing wall interactions, we achieve an increase of the magnitude of β to up to ≃1.1 for η wall =2.0. We are not able to find any velocity dependence of β, but find constant slip for fixed fluid-wall interactions, which is consistent with many experiments [13,24]. The fluctuations of the data for very low flow velocities are due to numerical uncertainties of the fit at very low curvature of the parabolic velocity profile. For V > 0.01 we find a slight deviation of β from the constant measurements. This is due to a small variation of the bulk pressure from P =0.097 for V = 1·10 −4 to P =0.106 for V = 0.03 that cannot easily be avoided for technical reasons. We have checked for a few data points that β stays constant if P can be kept at exactly fixed values, too. The slip length being independent of the flow velocity is consistent with many experiments and computer simulations, like the MD simulations of Cottin-Bizonne et al. [25] and the experiments of Cheng et al. [13] and Baudry et al. [26]. We speculate that an increase of β with increasing flow velocity as measured by some experiments [12] is due to surface roughness of the channel boundaries or other nonlinear effects. Since our model is not able to treat roughness on an atomic scale, we do not expect to conform with those results. MD simulations which find a non-constant value for β operate at very high shear rates which are orders of magnitude higher than what can be obtained by our approach [9]. Computing the exact slip in dependence of the interaction parameters from first principle analytically is a very hard or even impossible task since our interaction as given in Eq. 3 modifies the equilibrium distribution in the BGK operator. Therefore, we present a semianalytic approximation which utilizes the common two-layer model. Here, it is assumed that a thin fluid layer with thickness δ and different viscosity as the bulk fluid exists near the channel walls. As calculated by various authors [2], within this model the slip length can be computed as β = (µ bulk /µ 1 − 1)δ, where µ bulk is the viscosity of the bulk fluid, and µ 1 the viscosity close to the wall. Since the dynamic viscosity is given by the kinematic viscosity times the fluid density, µ = ρν = ρ(2τ α − 1)/6 [15], we write β = (ρ bulk /ρ 1 − 1)δ. ρ bulk can be measured in the channel's center and ρ 1 at the first lattice site next to the wall. Fig. 3 shows the dependence of ρ 1 on g wall,α for η wall =1.0, 5.0, P =0.11, and V =0.033. Since ρ 1 cannot easily be computed analytically, we postulate an interaction term that depends on the bulk density and the fluid-wall interaction as well as a free fit parameter k, and fit ρ 1 with an exponential function ρ 1 = ρ bulk (x, t) exp(−I). With only a single value for k we are able to utilize this equation to fit ρ 1 for all our simulation parameters. k is found to be 8.35 for our data. The lines in Fig. 3 illustrate the good quality of our approximation. A similar approach is applied to model the thickness of the layer at the wall which strongly depends on the fluid-wall interaction and bulk density. Here, we set δ = exp(I). As a result, β can be estimated by β = (exp(I)− 1) exp(I). The semi-analytic approximation is used to fit the dependence of the slip length β on the bulk pressure P . Fig. 4 shows the simulation data (symbols) and the approximation (lines) for wall interactions η wall =0.5, 1.0, and 2.0. The bulk pressure is varied from 0.03 to 0.33. We find a decrease of β with increasing pressure P . An increase of η wall leads to an increasing slope of the curves and to higher absolute values for β. Furthermore, we find a decrease of the slip with increasing bulk pressure. These results qualitatively agree with MD simulations [4,8]. Even with a single value for the fit parameter k, the semi-analytic description of β agrees very well for low fluid-wall interactions. For strong interactions (η wall =2.0), the fit qualitatively reproduces the behavior of the slip length. Higher order terms in the exponential ansatz for δ are needed for a better agreement.
To demonstrate that our approach is able to achieve experimentally available length and time scales, we scale our simulations to the experimental setup of Tretheway and Meinhart [11]. They use a 30µm high and 300µm wide microchannel with typical flow velocities of V = 10 −2 mm/s. For water, they measure a slip length of 0.92µm. The Reynolds number Re = 2hV /ν in their experiment is ≃0.3. To reproduce the observed slip, we set g wall,α =0.16 and η wall = 1.0 (see Fig. 1a). We are able to cover a wide range of flow velocities, i.e. for the setup given above, velocities can range from as low as 1 × 10 −4 to as high as 0.05 corresponding to Re between 0.038 and 19. The Knudsen number is given by Kn = ν/(c s 2h) which corresponds to 4.8 × 10 −3 for the simulations presented here. At these low Kn, the hydrodynamic approach is well valid. However, it has been shown that the LB method can be applied for Kn much larger than 1 if one uses modified boundary conditions [27]. Our mesoscopic force is expected to be able to properly describe fluid-wall interactions in such systems as well.
In conclusion, we presented a new approach to investigate boundary slip in hydrophobic microchannels by means of a multi-phase LB model. In contrast to MD simulations, our model is able to reach the length and time scales of typical experiments and is applicable for a wide range of realistic flow velocities. We qualitatively reproduced the dependence of slip on the hydrophobicity of the channel walls and found constant slip for varying flow velocities. The decrease of the slip with increasing pressure can be approximated by a semi-analytic approach. Our results are consistent with MD simulations [4,8,25] and experiments [11].