Research ArticleCHEMICAL PHYSICS

A novel physical mechanism of liquid flow slippage on a solid surface

See allHide authors and affiliations

Science Advances  27 Mar 2020:
Vol. 6, no. 13, eaaz0504
DOI: 10.1126/sciadv.aaz0504

Abstract

Viscous liquids often exhibit flow slippage on solid walls. The occurrence of flow slippage has a large impact on the liquid transport and the resulting energy dissipation, which are crucial for many applications. It is natural to expect that slippage takes place to reduce the dissipation. However, (i) how the density fluctuation is affected by the presence of the wall and (ii) how slippage takes place through forming a gas layer remained elusive. Here, we report possible answers to these fundamental questions: (i) Density fluctuation is intrinsically enhanced near the wall even in a quiescent state irrespective of the property of wall, and (ii) it is the density dependence of the viscosity that destabilizes the system toward gas-layer formation under shear flow. Our scenario of shear-induced gas-phase formation provides a natural physical explanation for wall slippage of liquid flow, covering the slip length ranging from a microscopic (nanometers) to macroscopic (micrometers) scale.

INTRODUCTION

The no-slip boundary condition of liquid flow is one of the most fundamental assumptions in fluid dynamics. However, there is no rigorous physical foundation for this condition. Therefore, its validation and violation have been discussed since the early days of the development of fluid mechanics (18). Recently, notable progress in experimental technology has unambiguously shown that the condition does not generally hold on a microscopic scale for a high enough shear rate. The widely used measure that characterizes the degree of the violation is the so-called slip length b, which is defined as the length at which vwall = b(∂v/∂y)wall, where v is the flow velocity, y axis is set to the direction normal to the wall, and the subscript “wall” means the value at the wall surface. For example, the slip length of water, which is the most typical and essential liquid and has thus been studied by many researchers, usually takes a value ranging from a few nanometers to 10 nm but sometimes reaches a few micrometers (2, 913).

There have been a large number of studies on the physical mechanism behind such flow slippage. For example, the nature of fluid slip has been studied from a microscopic aspect by molecular dynamics (MD) simulations, which take into account microscopic interactions between the fluid and solid atoms (14, 15). These studies have shown that a depletion layer of a few angstrom is formed between the fluid and the wall when the contact angle has a large value (e.g., a hydrophobic wall for water). Although a consistent scaling relationship has been established among the slip length, contact angle, and depletion layer, the physical mechanism behind flow slippage has remained elusive. For example, it is not clear how simple shear flow that is not coupled to the volume change can induce density fluctuations toward the formation of a gas layer. The shear rate used in MD simulations is typically 109 to 1012 s−1, which is higher than the one in experiments by many orders of magnitude. Furthermore, the slip length estimated in this way is far smaller than the typical experimental values, and thus, it is rather difficult to explain a slip length exceeding a few tens of nanometer. Although this approach is very transparent and attractive, thus, the relationship between the slippage phenomenon seen in microscopic MD simulations and that in macroscopic experiments is not so clear yet.

There is another type of theoretical approach, which assumed a priori that localized low-viscous gas layers (1620) or nanobubbles (2123) are formed on the wall. Once we accept this assumption, we may explain the apparent violation of the no-slip boundary condition, including the case of a large slip length. This scenario is rationalized by assuming the presence of nanobubbles on the wall (8, 20, 24). The existence of nanobubbles has been shown for a quiescent state by recent experimental studies [see (8) for review]. Consistent with this scenario, there has been experimental (25, 26) and numerical (27) evidence that dissolved gases in liquids play a crucial role in flow slippage of simple liquids; for example, it was shown experimentally that the fluid slip sensitively depends on the amount and kind of a dissolved gas (25, 26). We note that under ordinary experimental conditions, gas dissolution is inevitable. If the liquid is oversaturated with dissolved gas, gas bubbles will be spontaneously nucleated even in a quiescent condition if we wait for a sufficiently long time. That is, we may regard such a liquid to be metastable against gas formation.

The above scenario is attractive, and the surface enrichment of gas certainly takes place for hydrophobic surfaces. However, there have been debates on whether it is due to nanobubble formation or just the gas enrichment on the wall (15, 27, 28). Furthermore, there has so far been no satisfactory physical understanding of the impact of shear on nanobubble or gas-layer formation. For example, it is mostly unknown whether nanobubble formation under shear flow is the same as that in a quiescent state. The existence of stable nanobubbles was even paradoxical since a simple classical estimate suggests that they should dissolve quickly because of the large internal Laplace pressure coming from the large gas-liquid interface curvature (8). The stability of nanobubbles in a quiescent state has recently been explained by the pinning of their three-phase contact line. This stabilization mechanism suggests a possible fundamental difference between in a quiescent state and under shear. For example, it is not clear how the presence of a tiny amount of nanobubbles can lead to the formation of a rather thick gas layer on the wall. In addition, the formation of nanobubbles requires unique surface properties such as roughness and chemical inhomogeneity. This fact implies that the nucleation of the gas phase in bulk may suffer from a substantial barrier and is practically hard to take place. Then, the critical question is how shear flow affects the nucleation kinetics. This question is also related to another fundamental question of what determines the critical shear rate. The scenario based on nanobubbles seems to be difficult to explain the fact that flow slippage takes place even for walls of θ ≤ 90°.

Here, it is worth noting that slippage phenomena similar to those of simple liquids under flow were observed also for much more viscous fluids such as lubricants (29, 30), for which direct experimental evidence for shear-induced cavitation (or bubble formation) has been shown. The results of these two types of systems, simple nonviscous liquids and highly viscous liquids, apparently look quite different; at the same time, however, there is a commonality that both phenomena involve gas-phase formation. Thus, it is natural to consider that there may be a common physical mechanism behind the two phenomena. This motivates us to explain the flow slippage phenomena based on shear-induced cavitation or gas-layer formation on solid surfaces in a unified manner.

Here, we propose a novel physical mechanism of flow slippage based on the shear-induced liquid-to-gas transition in a liquid metastable against gas formation, taking into account the boundary conditions, including the friendliness (wettability) of the gas phase to the wall. To this end, we numerically study the mechanism of gas generation on a flat solid wall in a viscous metastable liquid, which mimics a liquid oversaturated with a gas. We assume that a metastable liquid obeys the mass conservation and compressible Navier-Stokes equations with thermal fluctuations and satisfies the no-slip boundary condition on a flat solid wall. We note that, in this thermal system, the thermally activated nucleation of gas bubbles should occur for a metastable liquid even in a quiescent state if we wait for a sufficiently long time. This incubation time can be enormously long, and thus, nucleation practically does not happen besides in the form of nanobubbles. This fact implies that there should be some mechanism of shear-induced reduction of the nucleation barrier. We investigate how the nucleation barrier for gas-bubble formation is reduced by shear flow and how shear-induced gas-bubble formation is affected by the presence of the solid wall and its surface wettability. Our study demonstrates that the amplitude of density fluctuation increases near the wall, and the application of shear further enhances the fluctuation and helps the formation of a gas layer on the wall. This indicates that the presence of nanobubbles is a sufficient condition to have flow slippage for low-viscosity liquids, such as water, but may not be a necessary condition. For very viscous liquids, for example, nanobubbles may be unnecessary, and spontaneous shear-induced slippage can take place even at a low shear rate.

RESULTS AND DISCUSSIONS

Here, we consider a van der Waals liquid, whose viscosity η weakly depends on the density ρ, near a gas-liquid transition as a model, and use a general coarse-grained model of compressible liquid in the framework of fluctuating hydrodynamics (31) (see Materials and Methods on the details). Unless otherwise stated, we set the homogeneous initial normalized density of the fluid as ψ0 = 0.480, which corresponds to a thermodynamically metastable state (see the Supplementary Materials for the van der Waals phase diagram). We introduce a solid wall, which has a contact angle θ with a liquid, i.e., the angle within the body of the liquid formed at the gas-liquid-solid interface, and apply simple shear, whose shear rate is γ̇.

Enhancement of density fluctuations near walls in a quiescent condition (at γ̇=0)

It is well known that the average density of ψ0(y) near the wall increases or decreases compared to the average density in bulk, depending on the contact angle θ. Here, the average density ψ0(y) is defined asψ0(y)ψ(x,y,t)x,t(1)where 〈⋯〉x and 〈⋯〉t are the averaging operations in terms of the coordinate x (along the wall) and the time t, respectively. Such behavior can be explained by the standard density functional theory (32). We show in the inset of Fig. 1A the spatial profile of the average density of ψ0(y) for various values of the contact angle θ. We can see that the density profile near the wall is crucially dependent on the contact angle θ, as expected. The liquid density near the wall increases (decreases) for θ ≤ 90° (θ ≥ 90°).

Fig. 1 Density fluctuation near the wall at a quiescent condition and under shear.

(A) The amplitude of density fluctuation near the boundary wall at γ̇=0. The red dashed line represents the thermal equilibrium value theoretically estimated for a bulk system without boundary walls. The inset shows the average density profile near the boundary wall. (B) Dependence of the normalized amplitude of density fluctuation on the distance to the boundary wall in the one-dimensional model B system. The results are for the case of the contact angle θ = 90°. The dashed curves are the theoretical prediction, Eq. 3. (C) Contact-angle dependence of the amplitude of density fluctuation near the boundary wall for γ̇=0.01. (D) Shear rate dependence of the amplitude of density fluctuation near the boundary wall for θ = 90°. Data are obtained in the metastable liquid state under steady shear.

Besides the well-known dependence of the density profile near the wall on θ mentioned above, we find a nontrivial influence of the wall on the spatial profile of the amplitude of density fluctuation, δψ(y) in a quiescent condition, i.e., for the shear rate γ̇=0. Figure 1A shows the dependence of the spatial profile of the amplitude of density fluctuation δψ(y) on the contact angle θ. Here, the density fluctuation δψ(y) are defined asδψ(y)ψ(x,y,t)2xψ(x,y,t)x2t(2)

We can see that the amplitude of density fluctuation at the wall is larger than that in bulk for any contact angle. As shown in Fig. 1, this enhancement can be seen even when θ = 90°, i.e., for the case of no preferential wetting. This behavior is unexpected since there is no fluid-wall coupling in the free-energy functional for θ = 90°. This implies that the enhancement should arise solely from the fluid-solid boundary condition itself.

Here, we explain the physical mechanism behind this wall-induced enhancement of density fluctuation. The equilibrium amplitude of thermal density fluctuation is determined by the competition between the strength of thermal fluctuation and the compressive resilience against the fluctuation. Now, we note that the boundary condition sets the influx and outflux through the boundary wall to be zero. Thus, the boundary condition restricts the direction of the compressive resilience only along the wall, leading to the anisotropy in the compressibility. This anisotropy effectively decreases the compressive resilience, resulting in the enhancement of density fluctuation near the wall than in bulk.

Since this enhancement mechanism is just a consequence of the constraint that the density flux is zero through the wall, we can consider this problem in the framework of model B (33, 34), where diffusion is the only transport mechanism, without losing the generality. By solving the one-dimensional problem, we can derive the density fluctuation near the wall analytically for θ = 90° as follows (see the Supplementary Materials for the derivation)δψ(y)=1+exp [yΛ]·δψbulk(3)where Λ is the decay length of density fluctuation and δψbulk is the equilibrium value of thermal fluctuation in bulk. Therefore, the amplitude of density fluctuation should be theoretically 2 times (i.e., by about 40%) larger at the wall surface than in bulk for θ = 90°. This prediction is approximately consistent with our simulation results (see Fig. 1B). The amplitude of density fluctuation near the wall is larger for a larger contact angle, as shown in Fig. 1A. Here, we note that the decay length Λ is proportional to the root of the compressibility for θ = 90° (see the Supplementary Materialss for the theoretical expression of Λ). We also note that this effect is enhanced by the application of shear, as shown in Fig. 1 (C and D).

Here, we emphasize that this enhancement of density fluctuation near the wall plays an essential role in shear-induced destabilization of a liquid against the gas-phase formation. This effect may be crucial to understand why the flow slippage takes place near the wall even for θ ≤ 90°.

Nontrivial effects of shear flow on density fluctuation

Now, we consider the impact of this enhancement of density fluctuation near the wall on shear-induced cavitation. For a metastable liquid, we should eventually observe inhomogeneous nucleation of gas bubbles via a thermal activation process if we wait for a sufficiently long time. According to the classical nucleation theory, the probability of the gas-phase nucleation pN is estimated as pN ∝ exp [− ΔG*/kBT], where ΔG* is the height of the free-energy barrier for nucleation. Considering the above boundary effect, which increases density fluctuation near the wall, it is natural to expect that the probability of inhomogeneous nucleation on the solid surface should be larger than the above bulk value. Here, we emphasize that this conclusion should hold even when θ ≤ 90°.

Keeping this in mind, we consider the effect of shear flow on density fluctuation. Usually, it is believed that shear flow is not coupled to density fluctuation simply because simple shear deformation does not cause any volume deformation. However, it was shown recently that this is limited to the case when the viscosity does not depend on the density (31, 35). Since the viscosity of a liquid generally depends on its density, shear flow can enhance density fluctuation through modification of the effective compressibility along compression and extension directions of shear deformation (31). When the imposed shear rate γ̇ is larger than the critical shear rate γ̇c=(η/P)T1 (with P being the pressure), the sign of the effective compressibility becomes negative, and the liquid intrinsically becomes unstable against gas-phase formation. For γ̇<γ̇c, on the other hand, the system remains metastable against cavitation. However, the amplitude of density fluctuation must increase even in such a case. We show in Fig. 1 (C and D) the contact-angle dependences of density fluctuation under shear flow in a metastable liquid and its shear rate dependence for the case of θ = 90°, respectively. We can see that shear flow further amplifies density fluctuation near the wall and makes the decay length of the amplitude toward bulk longer, which is due to the increase in the effective compressibility by shearing. This amplification of density fluctuation by shear flow means that the degree of supersaturation in the metastable liquid increases, and thus, the system effectively approaches its spinodal (instability) point. This leads to the reduction of the nucleation barrier ΔG*. Thus, inhomogeneous nucleation is accelerated by this flow-induced instability mechanism, and gas-bubble formation should occur on the wall much easier under shear flow than in a quiescent state. We expect that the application of shear can induce the nucleation of a much larger amount of the gas phase than nanobubbles formed in a quiescent state with the help of unique surface properties of the wall.

For a shear rate of γ̇ larger than the critical shear rate γ̇c, the sign of the effective compressibility becomes negative. In this case, we expect spontaneous liquid-to-gas phase transition, i.e., spontaneous formation of gas bubbles. The spinodal-type transition should lead to the exponential growth of density fluctuation in the early stage. We indeed see such exponential growth for γ̇<γ̇c, as shown in Fig. 2A.

Fig. 2 Bubble formation on the wall induced by shear flow.

(A) Temporal growth of density fluctuation on the boundary wall for θ = 90°. Data are obtained by averaging more than eight datasets simulated with different thermal noises. The solid curves represent the fitting results by using Eq. 4. (B) Amplification factor ω(γ̇) of density fluctuation. (C) Dependence of the characteristic shear rates γ̇ on the contact angle θ. The purple squares are obtained by the relation ω(γ̇c)=0 in (B). The purple dashed line represents the critical shear rate for the onset of instability of bulk liquid obtained from simulations in a system without walls (i.e., under a periodic boundary condition). We note that the theoretical values of the critical shear rate γ̇c predicted by Furukawa and Tanaka (31) and Steinberg et al. (35) are 0.0093 and 0.226, respectively. The upper white region corresponds to the unstable region under shear with walls. The orange and green regions correspond to two types of metastability, characterized by intermittent and steady behaviors, respectively (see below for the estimation of the black dashed line). (D) Time evolution of the slip length for an intermittent case (γ̇=0.010 and θ = 27.6°). Here, the time axis is t = tt0, where t0 is arbitrarily chosen in the time trajectory, in which the slip length is fluctuating. Thus, this figure shows only a part of the fluctuating trajectory. The state at t = t0 in the figure is the gas-liquid coexisting state, but later, the gas phase is torn by the Taylor mechanism and eventually dissolved into the liquid phase. The insets are the corresponding snapshots of the density distribution. The brighter regions have lower density.

Here, we note that in this simulation, thermal fluctuation is switched on at t = 0. Thus, the time development of density fluctuation should be a superposition of the thermal equilibration process and the exponential growth process of density fluctuation. The former process can be obtained by the time evolution of the fluctuation δψwall(γ̇=0;t) for γ̇=0 (see the lowest curve in Fig. 2A). Now, we denote the time evolution of density fluctuation on the wall under shear flow as δψwall(γ̇;t). For the spinodal decomposition–type instability, then, δψwall(γ̇;t) can be written asδψwall(γ̇;t)=exp [ω(γ̇)t]·δψwall(γ̇=0;t)(4)where ω(γ̇) is the amplification factor. The solid curves in Fig. 2A are the results of fitting of Eq. 4 to the data, using ω(γ̇) as the only fitting parameter. We conclude from the excellent fitting results that the spinodal-type instability can well describe the time evolution of density fluctuation in the initial stage. We also show the contact-angle dependence of ω(γ̇) obtained by the fitting in Fig. 2B. According to the theory of flow-induced instability for bulk systems (31), one can show that the amplification factor should satisfy the relation ω(γ̇)γ̇/γ̇c1, which we can confirm in Fig. 2B. Then, the critical shear rate γ̇c can be obtained from the relation ω(γ̇c)=0. The dependence of the critical shear rate estimated in this way is shown for various contact angles in Fig. 2C. We can see that the critical shear rate γ̇c monotonically decreases with an increase in the contact angle. Therefore, spontaneous bubble formation in a liquid should occur at a lower shear rate for a solid wall more friendly (or more wettable) to the gas phase (e.g., a hydrophobic wall in the case of water), as expected.

Nucleation of gas bubbles accelerated by shear flow

As discussed above, even without shear flow, gas bubbles will be formed by the thermal activation process for a metastable state if we wait for a sufficiently long time. So, we need to consider nucleation-growth–type gas formation in addition to spinodal decomposition–type one. We show in Fig. 3A a typical gas-bubble generation process observed in a metastable state. We can see that the system remains spatially homogeneous in the metastable state until gas-bubble generation happens (see left panel). At a particular moment, a gas-bubble nucleus appears and rapidly grows in its size. When a gas bubble is nucleated, it often appears only on one of the two walls in our simulations since the full growth of a bubble takes place before the next bubble nucleation event takes place in our simulation box size. We define this birth time of the first gas bubble as the “incubation time.” The dependence of the incubation time on the shear rate and the contact angle is shown in Fig. 3B. We can see that the incubation time for bubble formation on the walls decreases with an increase in the contact angle and/or the shear rate.

Fig. 3 Bubble formation on the walls in a metastable state under shear flow.

(A) A bubble formation process for γ̇=0.063 and θ = 135°: the metastable liquid (t = 1000) (left), immediately after bubble formation (t = 70,500) (middle), and a long time after the bubble formation (t = 82,500) (right). (B) Incubation time versus contact angle θ. Each point is obtained by averaging more than eight datasets, which are simulated with different sets of thermal noises. The black dashed curve represents the spinodal line, which is obtained from the purple square points in (C). Although in the spinodal region the system is thermodynamically unstable, there is apparently a finite incubation time in our criteria to detect the gas phase generation time. (C) Shear rate dependence of the incubation time for various contact angles. The dashed curves are guides to the eye. Data are obtained by using the same procedures discussed in the main text.

Here, we also show in Fig. 3C how the incubation time depends on the shear rate of γ̇ to consider the relevance of our scenario to a realistic situation. In Fig. 3C, the vertical axis is the same as Fig. 3B, but the horizontal axis is the shear rate. For large contact angles, we can see that the incubation time tinc decays monotonically with an increase in γ̇ in this weak shear regime: The decay rate α(θ) increases with a decrease in the contact angle θ. We can at least infer that α(θ) steeply increases when θ approaches 90°. Therefore, the incubation time should drastically decrease with an increase in γ̇ in a very small shear rate regime. Furthermore, for large θ, gas bubbles are to be formed even at a quiescent condition if we wait for a long enough time. This fact indicates that our mechanism should be valid even for realistic shear rates used in experiments, implying that nucleation-growth–type gas formation is practically much more important than spinodal decomposition–type one. We will discuss the situation relevant for experiments later.

Behavior of the gas phase under shear flow

If we stop shearing immediately after the nucleation of the gas phase and wait for a sufficient time to reach the thermodynamic equilibrium, the gas phase takes a stable equilibrium shape on the wall, depending on the contact angle. Under a steady shear flow, however, the structure of the gas phase is distorted from its equilibrium shape by the Kelvin-Helmholtz (KH) instability mechanism, which is the hydrodynamic instability occurring when two fluid layers with different densities move with the different horizontal speeds (36).

We investigate the temporal change of the shape of the gas phase under shear flow by performing the following simulations. We prepare the thermal equilibrium state of the gas phase coexisting with the liquid phase on the wall as the initial state. Then, we apply shear flow at t = 0 to investigate the temporal change of the gas-phase shape induced by shear flow.

Time evolution of the gas-phase shape under shear. First, we show an example of the time evolution of the gas-phase shape under shear flow in Fig. 4A for γ̇=0.00010 and θ = 27. 6°. In this case, the gas bubble peels off from the wall and goes into the bulk liquid without disintegration of the droplet. For a sufficiently high shear rate, on the other hand, the gas phase peels off from the wall and splits into several pieces, which eventually dissolve into the bulk liquid and disappear (see Fig. 2D).

Fig. 4 Time evolution of a gas-phase droplet under shear flow.

(A) The results are for γ̇=0.00010 and θ = 27.6°. The fluctuation in the background is due to thermal noise (see movie S1). (B) The results are for γ̇=0.00016 and θ = 90.0° (see movie S2). (C) The results are for γ̇=0.0040 and θ = 90.0° (see movie S3). (D) The results are for γ̇=0.0016 and θ = 150.8° (see movie S4). (E) The results are for γ̇=0.0040 and θ = 150.8° (see movie S5). (F) Dependence of the critical shear rate for detachment of the gas phase on the wall on the contact angle. The green filled circles indicate the shear rate above which the gas phase is peeled off from the wall while keeping a spherical shape. The purple squares indicate the shear rate above which the SD of the wall coverage fluctuation exceeds 0.03. The critical shear rate is not sensitive to this choice of the threshold coverage.

Next, we show the case for the contact angle θ = 90° and γ̇=0.0016 in Fig. 4B. We can see that the shape of the gas-phase droplet changes from the initial hemispherical shape to the pancake-like elongated shape. This shape change is induced by the shear stress imposed on the gas-phase droplet. In this case, the gas phase contacts with the wall over a larger area under shear flow than in a quiescent state. Moreover, the gas phase keeps the pancake-like shape rather steadily. For a sufficiently strong shear rate, on the other hand, the gas-phase droplet is strongly deformed by the Taylor mechanism and destabilized by the KH instability mechanism, which leads to the transient dissolution and disappearance of the droplet after a sufficiently long time (see Fig. 4C).

Last, we show the case for a large contact angle, θ = 150. 8° and γ̇=0.0016, in Fig. 4D. In this case, when the shear flow is applied, the gas phase takes a layered shape. When the shear rate is sufficiently high, however, the layer structure is destabilized by the KH instability, as shown in Fig. 4E. We confirm that the fluctuation increases as the shear rate increases.

Thus, for a high shear rate regime, the gas-phase structure is strongly perturbed and detached from the wall by shear flow. For the realization of a large slip length, the presence of a sufficient amount of the gas phase on the wall is necessary. For a small contact angle, the detachment of the gas phase takes place easily, even for weak shear. We summarize the critical shear rate for detachment or large temporal fluctuations in Fig. 4F. The black dashed curve separating the steady and intermittent regions in Fig. 3C is drawn on the basis of the critical shear rate plotted in Fig. 4F. In the steady state, there is a well-defined stable slip length, whereas in the intermittent state, a slip length is not well defined and fluctuating with time as indicated in Fig. 2D.

Theoretical estimation of the critical shear rate of the KH instability. Here, we consider a mechanism by which a two-dimensional gas layer formed on a wall is destabilized under shear flow. This mechanism is known as the KH instability. We consider a situation that the gas and liquid phases coexist with the flat interface parallel to the wall located at y = 0. We also assume that the gas and liquid layers exist in the upper half plane and lower half plane, respectively. Then, we apply a simple shear flow, whose spatially averaged velocity field is given by vxγ̇y, to the system. Now, we study the response of the interface to an infinitesimally small perturbation of a vertical displacement in the form of η ∝ exp [ikx(xct)].

This problem has already been discussed in several seminal papers (3740). Here, we follow the discussion in (39) and estimate the critical shear rate γ̇KH under the assumption of a small Reynolds number ργ̇/ηkx2<1, which is satisfied in our simulations. The vertical perturbation is amplified by the advection of disturbance vorticity with the growth rate of ωgrowthO(ργ̇2/kx2η). On the other hand, the surface tension σLG tends to suppress the vertical perturbation with the characteristic decay rate of ωdecayOLGkx/η). These two competing effects are balanced when ωgrowth ∼ ωdecay. Thus, we obtain the critical shear rate for the KH instability as followsγ̇KHkx3σLGρ(5)

Now, we calculate the actual value of the critical shear rate by substituting the parameters used in our simulations. The surface tension is estimated as σLG = 0.0029 by using the standard density functional theory (41). By setting the characteristic perturbation wave number as kx = 2π/Lx, the critical shear rate is estimated as γ̇KH=0.00089. We note that this critical shear rate γ̇KH satisfies the small Reynolds number condition ργ̇/ηkx2<1. The above theoretical prediction almost coincides with the value γ̇=0.0022 obtained by numerical simulations for θ = 180° (see Fig. 4F).

Stability of gas bubbles formed on walls in a metastable state

Destabilization of gas bubbles and the resulting intermittent behavior. Destabilization of a gas phase by shear takes place not only for a flat gas layer but also for gas bubbles. For example, gas bubbles are deformed and split into small ones by shear. Thus, we consider the stability of gas bubbles generated in a metastable state under shear flow. If the generation and growth of a gas bubble occur in a quiescent condition, the bubble should be stable and exist as it is. This is simply because the final state is in thermal equilibrium. Under a steady shear flow, however, after the formation of gas bubbles, they are repeatedly broken up and coalesce unless the shear rate is very low (γ̇0.004). Gas bubbles have the stable form of a wetting layer homogeneously covering the wall under shear flow only when the shear rate is low enough (low γ̇), and the gas phase is wettable to the wall (large θ). We refer this situation to the steady state of the metastable region. Otherwise (i.e., large γ̇ or small θ), gas bubbles nucleated detach from the wall, dissolve, and disappear under shear flow. In this case, thus, we see successive disappearance and emergence of gas bubbles near the boundary walls, and thus, we refer this situation to the intermittent state of the metastable region. We show in Fig. 2C the state diagram, where the state is first divided into the unstable and metastable regions, and the latter is further divided into the intermittent and steady states.

Origin of the intermittent behavior. To understand the mechanism of the intermittent behavior, we study the dynamical behavior of the gas phase under steady shear flow. We show an example of the time evolution of the density field for γ̇=0.0010 and θ = 26. 7° in Fig. 2D. We can see there that a gas bubble formed on the wall is squeezed, torn, and lastly dissolved. Two candidate mechanisms may be responsible for this behavior: One is the Rayleigh-Taylor mechanism (34, 42), which leads to the emergence of tearing shear stress due to the difference in the viscosity between the liquid and gas, and the other is the advection due to the advective term in the Navier-Stokes equation. It is well known that the advective term stirs the system; however, this nonlinear term may not be so crucial because the Reynolds number is very small in our simulations. Once the tearing shear stress splits the gas bubble to pieces smaller than the critical gas-bubble nucleus size, they dissolve, and accordingly, the system goes back to a homogeneous state. Detachment of a gas bubble from the wall, which more easily happens for the case of lower gas wettability to the wall, also enhances their dissolution. It is the origin of the intermittent behavior. Here, we note that, in the intermittent state, the slip length b also fluctuates with time, as shown in Fig. 2D. When the wettability of the gas phase to the wall is high enough and the shear rate is low enough, the stabilizing force coming from the solid- and liquid-gas interfacial energy wins over the destabilizing force due to Taylor mechanism, leading to the emergence of a steady state. Thus, we can conclude that the formation of gas bubbles and their survival under shear is both easier under a lower shear rate and for a larger contact angle. In this steady state, the slip length b does not fluctuate largely with time.

Slip length

Next, we estimate the slip length for shear-induced flow slippage described above.

Slip in the steady state. In Fig. 5A, we show the y dependence of the shear velocity vx(y) = 〈vx(x, y, t0)〉x and the effective shear viscosity η(y) = 〈η(x, y, t0)〉x at a certain moment t = t0 when the steady state is reached in the metastable region after gas-bubble generation and growth. The decrease in the averaged effective shear viscosity near the wall indicates the nonlinear shear velocity profile vx(y) near the wall. We can see the apparent slip behavior by looking at the deviation of the profile from the red dashed line in Fig. 5A. We also show the dependence of the time-averaged apparent slip length on the shear rate and the contact angle in Fig. 5B. The slip length depends on both the shear rate and the contact angle, although its dependence is rather weak. Here, we take the time averaging operation only for the state with gas bubbles near the wall. We note that the way of this time average has a considerable effect on the estimation of slippage for the intermittent state. Although the slip length itself rather smoothly depends on the shear rate and the contact angle, the amplitude of the slip length fluctuations shows a distinct transition across the border between the steady and intermittent states (see Fig. 2C).

Fig. 5 Apparent slip of fluid under shear.

(A) An example of the y dependence of the averaged shear velocity vx(y) ≡ 〈vx(x, y, t0)x〉 after gas-bubble growth in the steady state for γ̇=0.0010 and θ = 135°. The green and purple symbols are for the initial densities of ψ0 = 0.46 and ψ0 = 0.48, respectively. The small averaged viscosities of the gas phase formed on the wall lead to finite slip lengths indicated by the double-headed arrows. The insets show snapshots of the density distribution for ψ0 = 0.46 and 0.48. (B) Dependence of the time-averaged apparent slip length on the contact angle θ ≥ 90.0° for the shear rate γ̇ and the initial densities of ψ0 = 0.46, 0.47, and 0.48. The dashed arrows indicate the theoretically predicted slip lengths for θ = 180.0°, which are estimated from the thickness of the gas layers on the basis of the thermodynamic lever’s rule and the balance of the shear stress. (C) Time dependence of the slip length for spinodal states. The results are shown for (γ̇,θ)=(0.0398,27.6), (0.0398, 90.0°), and (0.0100, 169.4°). (D) Snapshots of density distribution for spinodal states. These figures correspond to the results in (C).

Slip in the unstable regime. Here, we note that in the unstable region, the flow slippage also takes place, but the slip length strongly fluctuates. We investigated the time dependence of the slip length and density distribution for the spinodal region. We select the shear rate and contact angle as (γ̇,θ)=(0.0398,27.6), (0.0398,90.0°), and (0.0100,169.4°). These parameter sets correspond to spinodal (i.e., unstable) cases (see Fig. 2C). The slip lengths and the snapshots for the parameter sets are given in Fig. 5 (C and D, respectively). We can confirm that the slip lengths have large temporal fluctuations, and the density distributions are much disturbed.

System size effect on the slip length. Until now, we have shown the results of the slip length for the system size of Lx = Ly = 128. However, we find that when we change the system size while keeping the box shape, the apparent slip length is proportional to the size of the system. The reason is that the ratio of the total volume of generated gas bubbles localized near the wall to that of the system is determined by the lever rule in a nonequilibrium steady state. Therefore, even if the metastability of the liquid is weak, the slip length can become large for a system of large size. It may explain some experimental reports of large slip lengths (∼a few micrometers). The validity of our mechanism may be checked experimentally by examining the system size dependence of the slip length. It should be noted that the slip length should also depend on the amount of dissolved gas. This may explain a broad distribution of the slip lengths reported even for the same liquid-substrate pairs.

Comparison of physical parameters between simulation and real systems

Here, we make a connection of our simulation results to experimental observation.

Spinodal decomposition–type gas-phase formation. As in many other simulations, there is a large gap in the time scale between the simulations and experiments. We estimate the actual physical values of the viscosity coefficient ηc, unit length l, and unit time τ used in our simulation by substituting typical parameters for liquids like water. In our simulation, we impose the following constraints: v01/3/l=0.1, ηc=7lm kB T/v0, and τ = (ηcv0)/(5kBT) between numerical parameters. Here, ηc is the viscosity, m is the molecular mass, v0 is the molecular volume, and T is the temperature. By using these relations and setting T = 300 K, m = 3.0 × 10−26 kg, and v0 = 3.0 × 10−29 m3, we obtainηc8.1 × 103 Pa·s,l3.1×109 m,τ1.2 × 1011s(6)

In particular, the viscosity coefficient η and the critical shear rate γ̇c=(η/P)T1 for the liquid-side thermal equilibrium density ψ = 0.495 are given byη(ψ=0.495)1.5×102 Pa·s,γ̇c2.1×109 s1(7)

On the other hand, it is well known that the experimental value of the viscosity coefficient for pure water at ambient condition is about 1 mPa·s. The value of γ̇c of pure water can be estimated from the pressure dependence of the viscosity, which was measured by many researchers (4345): γ̇c1012 s−1 at ambient condition.

There are not many experimental studies on the viscosity of water saturated with gas, but there are several studies on the aqueous solution of carbon dioxide (4648). The value of γ̇c of this system estimated from the pressure dependence of the viscosity at relatively high pressure range (from 10 to 100 MPa) is about 1012 s−1 (46, 48). For a low pressure range (from 0.7 to 5 MPa), on the other hand, the value of γ̇c is estimated to be about 1010 s−1 (47), although large experimental error prevents accurate estimation. Anyway, γ̇c is much lower at lower pressure than at higher pressure.

The time unit in our simulations is ∼12 ps, and the unit of the shear rate is thus ∼8 × 1010 s−1. So, the shear rate in our simulations is quite high compared to experiments that reported flow slippage, which is comparable to MD simulations [see, e.g., (15)]. In any case, these estimations indicate that the spinodal decomposition–type gas-phase formation is not expected for low viscosity liquids such as water in a realistic experimental condition.

However, this is not necessarily the case for very viscous liquids near the glass transition (31, 49). Now, we discuss the value of the unit time τ for highly viscous liquids. In simulations, we investigated the instability mechanism by solving the Navier-Stokes equation under a low Reynolds number condition. However, it is expected that the same conclusion should be obtained even under the Stokes approximation because the inertia term in the Navier-Stokes equation does not play an essential role in our instability mechanism. Under the Stokes approximation, it is unnecessary to assume the relation ηc=7lm kB T/v0, and we can treat the viscosity coefficient as an independent parameter. This means that we can increase the unit time τ for highly viscous liquids. Table 1 lists the unit time τ for a phenyl ether oligomer (5P4E) and a metallic glass (Vitreloy-1) near their glass transition temperatures, which are estimated from the experimental data. The large unit time τ implies that the spinodal-type instability may be relatively easily induced even at a rather low shear rate, which can be realized experimentally, for highly viscous liquids. This may provide an interesting opportunity to check our prediction directly. Shear-induced cavitation phenomena were experimentally studied for lubricant oligomers (29, 30). Surface effects on such instability is a fascinating problem for future experimental investigation. Flow slip in these viscous liquids should lead to a drastic reduction of energy dissipation in their hydrodynamic transport, which may have a considerable impact on various applications.

Table 1 Simulation unit time τ for phenyl ether oligomer (5P4E) and metallic glass (Vitreloy-1).

The pressure dependence of viscosity for 5P4E and Vitreloy-1 are calculated by the empirical formula in (29) and (54), respectively. We assume the particle volume v0 ∼ 10−29m3.

View this table:

Nucleation-growth–type gas-phase formation. For nucleation-growth–type gas-phase formation, the story is quite different, since nanobubbles may exist even in a quiescent state without shear. For example, it is widely known that a liquid with dissolved gas, which is used in most of real experiments, tends to form nanobubbles on the wall if it is wettable to the gas phase (e.g., a hydrophobic wall for water) in a quiescent state (8). This fact seems to suggest that the gas layer may be easily formed by shear, which leads to the flow slippage [see, e.g., (8, 20, 24)].

The metastability should lead to the formation of such nanobubbles if we wait for a long enough time. Thus, the situation looks obvious, and no further consideration is necessary to understand the flow slippage. However, this is not the case. We stress that the existence of stable nanobubbles itself is not so obviously justified since a simple classical estimate suggests that they should dissolve quickly due to the large internal Laplace pressure coming from the large gas-liquid interface curvature (8). The stability of nanobubbles in a quiescent state has recently been explained by the pinning of their three-phase contact line. This stabilization mechanism suggests a possible fundamental difference between in a quiescent state and under shear. For example, it is not clear at all how the presence of a tiny amount of nanobubbles can lead to the formation of a rather thick gas layer on the wall. Thus, for example, we do not know the amount of the gas phase to be formed under shear. Furthermore, we should note that the formation of nanobubbles requires unique surface properties such as surface roughness and chemical inhomogeneity. This fact means that the actual free-energy barrier for gas bubble nucleation must be extremely high, and the nucleation does not happen in a realistic time scale besides nanobubbles. Even in such a situation, the shear flow should play a crucial role in gas-phase formation and the resulting slippage. Our scenario provides a simple physical picture on the impact of such shear effects on the formation of gas bubbles or gas layers.

As we show in Fig. 3 (B and C), there is a steep reduction of the incubation time of gas-bubble nucleation as a function of the shear rate, which is particularly the case for θ ∼ 90°. If we use the parameters for water, the time unit in this plot approximately corresponds to τ ∼ 10−12 s, and the unit of the shear rate is ∼109 s−1 for Fig. 3 (B and C). Thus, for example, for θ ∼ 140°, we estimate the incubation time of spontaneous formation of gas bubbles on the walls to be about 0.01 to 0.1 s, although this estimation is very rough. For lower θ, this incubation time should become even much longer and may exceed the experimental time scale. The crucial point is that our scenario can explain the occurrence of the flow slippage even for θ < 90°, despite that the gas phase is less wettable for the wall in this condition. According to Fig. 3 (B and C), the dependence of the incubation time on the shear rate is very steep, particularly for a smaller contact angle. Thus, we expect that the shear-induced cavitation takes place even for a very small shear rate used in experiments within a reasonable observation period. It means that there may practically be an onset shear rate even for a metastable state [see, e.g., (26)]. Since the amount of a gas phase under shear is not constrained by the pinning of the three-phase contact line, it should simply follow the lever rule. Thus, it can even explain a macroscopic slip length if a sample size is large. Thus, we may say that the most natural scenario of flow slippage in a low-viscosity liquid may be nucleation-growth–type shear-induced cavitation. Here, we note that the density dependence of the viscosity η or the functional form of η(ρ), affects the strength of the shear effects on the gas-phase formation, as shown in the Supplementary Materials.

Furthermore, the state of the gas phase nucleated under shear flow is not necessarily in a steady state in the form of gas layers on the wall. The interaction of the gas phase formed on the wall with the shear flow can lead to complex behavior, as shown in Figs. 2 (C and D) and 4. The intermittent behavior should be observed in a broad region of the contact angle shear rate state diagram, which may be a cause of considerable fluctuations of the slip length observed in experiments (12, 13).

CONCLUSIONS

To summarize, we propose a novel physical mechanism of flow slippage originating from intrinsic enhancement of density fluctuation near solid surfaces due to the boundary condition and the coupling between shear stress and density through the density dependence of the viscosity η(ρ). Our theory provides essential information on how we can control the flow slippage, including the case of a metastable liquid saturated by gas, which we mainly discuss in this article. However, the relevance of our scenario to real experiments should be carefully examined. We note that the information on the amount of dissolved gas is crucial for such a comparison. The comparable shear rate required for the formation of slippage between our coarse-grained simulations and MD simulations indicates that, without dissolved gas, the experimentally unrealistic and exceptionally high shear rate is necessary for shear-induced gas-layer formation. We have also shown that this mechanism should also be essential for a viscous liquid near its glass transition point since the density dependence of the viscosity is very strong there (see the preceding section). In such a case, even spinodal decomposition–type shear-induced cavitation can take place preferentially on solid walls.

Our mechanism can naturally explain experimental reports of rather large slip lengths, which have so far been difficult to explain theoretically. For example, the shear rate–dependent slip lengths reported in (11) can be rather naturally understood by our mechanism, including complex, nonideal behaviors. It is known that the large-length slip experimental results largely fluctuate with time and depend on the spatial position of their measurements (50). Such features have been quite difficult to be explained by previous models. According to our results, the shape of gas bubbles near the wall is neither uniform nor steady in the metastable intermittent state and fluctuates spatiotemporally, even for a flat homogeneous solid wall. In the above, we have discussed the gas generation process on a completely flat wall. However, actual solid walls must have surface roughness and chemical inhomogeneity. Impact of such inhomogeneities on flow slippage phenomena and their application to the control of slippage (5053) remains an interesting problem for the future investigation.

In this work, we focus on the slippage of a simple viscous fluid, in which the flow-induced instability stems from the density dependence of the viscosity. It has been shown that complex fluids also exhibit similar flow-induced instability in bulk (54), which resembles the shear-induced instability (5557) of dynamically asymmetric complex fluids, which is a mixture of large (slow) and small (fast) components (58). The shear-induced instability in such complex fluids and simple viscoelastic liquids can be understood in a unified manner by using the mapping of c ↔ ρ and Π ↔ p (31), where c is the concentration and Π is the osmotic pressure in the complex fluids. However, rather a few studies have been conducted on the shear-induced instability of complex fluids influenced by a boundary wall (5967). In this case, the liquid-rich dilute phase is to be formed near the walls for such a mixture to reduce viscous dissipation, as the gas phase is formed for a simple liquid. We believe that the instability mechanism described here may be applied commonly to shear-induced instability of viscoelastic materials and complex fluids under the influence of walls.

MATERIALS AND METHODS

Simulation method

Here, we consider a simple viscous liquid by using the following general model of fluid dynamics (31). The time evolution of density ρ(r, t) is given by the continuity equationρt=i(ρvi)(8)where vi is the fluid velocity. The time evolution of fluid velocity is given by the momentum conservation equation, i.e., the Navier-Stokes equationρ(t+vjj)vi=jσijV+jσijRjΠij(9)where σijV, σijR, and Πij are the viscous stress, thermal noise stress, and pressure tensors, respectively. The viscous stress tensor σijV is given byσijV=η(ρ)(ivj+jvi2dδijmvm)+ζ(ρ)δijmvm(10)

Here, we assume that the ρ dependence of the shear viscosity η(ρ) and the bulk viscosity ζ(ρ) as η(ρ) = ηd · [1 + 6(ρv0/m) + 36(ρv0/m)2]/7 and ζ(ρ) = 8η(ρ)/15, respectively. To see the effect of the functional form of the density dependence of η(ρ), we also use another functional form (Supplementary Materials). The thermal noise stress σijR satisfies the fluctuation-dissipation relationσijR(r,t)σmnR(r,t)=2kBTδd(rr)δ(tt)×[η(ρ)(δimδjn+δinδjm)+[ζ(ρ)2dη(ρ)]δijδmn](11)

We use the Ginzburg-Landau–type free energy functional for the liquid, including the fluid-wall interaction=dV[f(ρ)+C2ρ2]+walldSϕS(ρ)(12)

We use the van der Waals formula f(ρ) = − εv0(ρ/m)2 − (ρ/m)kBT ln (mv0 − 1) as the bulk free-energy density f(ρ), which allows us to describe the gas-liquid phase transition properly. Here, m, v0, and ε are the molecular mass, excluded volume, and strength of attractive interaction, respectively. We use the van der Waals formula as a simple model equation that can describe the transition between the gas and liquid phases. We set the fluid-wall surface interaction as ϕS(ρ) = kρ, where k is a constant. The pressure tensor Πij is obtained by the functional derivative of ℱ. We use the normalized density, unit length, and unit time as ψ(r, t) = ρ(r, t)v0/m, l = (m2C/v0kBT)1/2, and τ = ηcv0/(5kBT), respectively. The velocity is scaled by l/τ. Furthermore, we set v01/3/l=0.1,ηcv0/(7lkBTm)=1 and 1/ε˜=kBT/ε=0.28. These scale factors and parameters are the same as those in (31).

We consider two-dimensional liquids (the spatial dimension d = 2) in this work, but the result should be essentially the same for d = 3. The system has a square shape whose horizontal and vertical lengths are Lx = Ly = 128. The system is sandwiched by two parallel flat solid walls, and the flow direction is the x direction. We impose a periodic boundary condition to x direction. The contact angle θ between the gas phase, the liquid phase, and the solid wall is analytically calculated from the density functional theory (32, 41). The velocity of the top and bottom walls are given by vupper=(γ̇Ly/2,0)t and vlower=(γ̇Ly/2,0)t, respectively, where γ̇ is the shear rate. We impose the no-slip boundary condition at the wall surfaces. Numerical simulations are conducted on a staggered lattice. The contribution from the pressure tensor (or chemical potential) to the time evolution equations near the boundary walls is derived by discretizing the free-energy functional and using the Onsager’s principle (68).

When we study the stability of gas bubbles formed on the wall under shear, we calculate the time evolution of the system while switching off the thermal noise after the gas bubble generation under thermal fluctuation. We detect the emergence of gas bubbles from a metastable liquid after a sufficiently long time, which is defined as t = 0. The incubation time of gas bubbles is defined as follows. First, we define the “gas phase” as lattices having a density less than ψ = 0.247. We calculate the time evolution of the system until the total area fraction of the gas phase exceeds 3%. If the area fraction exceeds 3%, then we go back in time and define the moment that corresponds to the gas area fraction equal to 0.03% as “gas bubble generation time” t = 0. We note that our results are not sensitive to the selection of the above thresholds, 3 and 0.03%.

SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/13/eaaz0504/DC1

Supplementary Text

Fig. S1. The van der Waals equation of state used in the main text.

Fig. S2. Normalized density dependence of the viscosity coefficient.

Fig. S3. Shear rate dependence of the incubation time for viscosity functions.

Movie S1. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0001_theta = 27.6).

Movie S2. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0016_theta = 90.0).

Movie S3. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0040_theta = 90.0).

Movie S4. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0016_theta = 150.8).

Movie S5. Time evolution of a gas-phase droplet under shear flow (ShearRate = 0.0040_theta = 150.8).

Reference (69)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

REFERENCES AND NOTES

Acknowledgments: We thank A. Furukawa for valuable discussions and suggestions. Funding: This work was partially supported by Grants-in-Aid for Specially Promoted Research (JP25000002) and Scientific Research (A) (JP18H03675) from the Japan Society of the Promotion of Science (JSPS) and by the Mitsubishi Foundation for support. Y.K. acknowledges the JSPS for financial support (JP15J00018). Author contributions: H.T. conceived and supervised the project. Y.K. performed the simulations. Y.K. and H.T. discussed the results and wrote the manuscript. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article