Physics | Hydrodynamics » Molecular Hydrodynamics of the Moving Contact Line in Two-Phase Immiscible Flows

COMMUNICATIONS IN COMPUTATIONAL PHYSICS Vol. 1, No 1, pp 1-52 Commun. Comput Phys February 2006 REVIEW ARTICLE Molecular Hydrodynamics of the Moving Contact Line in Two-Phase Immiscible Flows Tiezheng Qian1,∗ , Xiao-Ping Wang1 and Ping Sheng2 1 Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China. 2 Department of Physics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China. Received 1 June 2005; Accepted (in revised version) 17 September 2005 Communicated by Weinan E Abstract. The no-slip boundary condition, ie, zero fluid velocity relative to the solid at the fluid-solid interface, has been very successful in describing many macroscopic flows. A problem of principle arises when the no-slip boundary condition is used to model the hydrodynamics of immiscible-fluid displacement in the vicinity of the moving contact line, where the interface separating two immiscible

fluids intersects the solid wall. Decades ago it was already known that the moving contact line is incompatible with the no-slip boundary condition, since the latter would imply infinite dissipation due to a non-integrable singularity in the stress near the contact line. In this paper we first present an introductory review of the problem We then present a detailed review of our recent results on the contact-line motion in immiscible two-phase flow, from molecular dynamics (MD) simulations to continuum hydrodynamics calculations. Through extensive MD studies and detailed analysis, we have uncovered the slip boundary condition governing the moving contact line, denoted the generalized Navier boundary condition. We have used this discovery to formulate a continuum hydrodynamic model whose predictions are in remarkable quantitative agreement with the MD simulation results down to the molecular scale. These results serve to affirm the validity of the generalized Navier boundary condition,

as well as to open up the possibility of continuum hydrodynamic calculations of immiscible flows that are physically meaningful at the molecular level. Key words: Moving contact line; slip boundary condition; molecular dynamics; continuum hydrodynamics. Correspondence to: Tiezheng Qian, Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China. Email: maqian@usthk ∗ http://www.global-scicom/ 1 c °2006 Global-Science Press 2 1 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 Introduction The no-slip boundary condition, i.e, zero relative tangential velocity between the fluid and solid at the interface, serves as a cornerstone in continuum hydrodynamics [3]. Although fluid slipping is generally detected in molecular dynamics (MD) simulations for microscopically small systems at high flow rate [2, 8, 39, 41], the no-slip boundary condition still works well for macroscopic flows at low flow rate.

This is due to the Navier boundary condition which actually accounts for the slip observed in MD simulations [2, 8, 39, 41]. This boundary condition, proposed by Navier in 1823 [30], introduces a slip length ls and assumes that the amount of slip is proportional to the shear rate in the fluid at the solid surface: £ ¤ vslip = −ls n · (∇v) + (∇v)T , slip is the slip velocity at the surface, measured relative to the (moving) wall, [(∇v)+ where v ¤ (∇v)T is the tensor of the rate of strain, and n denotes the outward surface normal (directed out of the fluid). According to the Navier boundary condition, the slip length is defined as the distance from the fluid-solid interface to where the linearly extrapolated tangential velocity vanishes (see Fig. 1) Typically, ls ranges from a few angstroms to a few nanometers [2,8,39,41] For a flow of characteristic length R and velocity U , the slip velocity is on the order of U ls /R. This explains why the Navier boundary condition is

practically indistinguishable from the no-slip boundary condition in macroscopic flows where v slip /U ∼ ls /R 0. It has been well known that the no-slip boundary condition is not applicable to the moving contact line (MCL) where the fluid-fluid interface intersects the solid wall [11, 12, 20] (see Fig. 2 for both the static and moving contact lines). The problem may be simply stated as follows In the two-phase immiscible flow where one fluid displaces another fluid, the contact line appears to “slip” at the solid surface, in direct violation of the no-slip boundary condition. Furthermore, the viscous stress diverges at the contact line if the no-slip boundary condition is applied everywhere along the solid wall. This stress divergence is best illustrated in the reference frame where the fluid-fluid interface is time-independent while the wall moves with velocity U (see Fig. 2b) As the fluid velocity has to change from U at the wall (as required by the no-slip boundary

condition) to zero at the fluid-fluid interface (which is static), the viscous stress varies as ηU/x, where η is the viscosity and x is the distance along the wall away from the contact line. Obviously, this stress diverges as x 0 because the distance over which the fluid velocity changes from U to zero tends to vanish as the contact line is approached. In particular, this stress divergence is non-integrable (the integral of 1/x yields ln x), thus implying infinite viscous dissipation. Over the years there have been numerous models and proposals aiming to resolve the incompatibility of the no-slip boundary condition with the MCL. For example, there have been the kinetic adsorption/desorption model by Blake [4], the fluid slip models by Hocking [19], by Huh and Mason [21], and by Zhou and Sheng [43], and the Cahn-Hilliard-van der Waals models by Seppecher [38], by Jacqmin [24], and by Chen et al. [7] In the kinetic adsorption/desorption model by Blake [4], the role of molecular

processes was investigated. A deviation of the dynamic contact angle from the static contact angle was shown to be responsible for the fluid/fluid displacement at the MCL. In the slip model by Hocking [19], the effect of a slip coefficient (slip Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 v slip ls 3 fluid solid Figure 1: Slip length introduced in the Navier boundary condition, defined as the distance from the fluid-solid interface to where the linearly extrapolated tangential velocity vanishes. length) on the flow in the neighborhood of the contact line was examined. Two slip models were used by Huh and Mason [21]: Model I for classical slippage assumes the slip velocity is proportional to the shear stress exerted on the solid; Model II for local slippage assumes that over a (small) distance the liquid slips freely where fluid stress vanishes, but thereafter the liquid/solid bonding has been completed and the no-slip boundary condition is applied. In the

slip model by Zhou and Sheng [43], the incompressible Navier-Stokes equation was solved using a prescribed tangential velocity profile as the boundary condition, which exponentially interpolates between the complete-slip at the contact line and the no-slip far from the contact line. The Cahn-Hilliard-van der Waals models by Seppecher [38], by Jacqmin [24], and by Chen et al. [7] suggested a resolution when perfect no-slip is retained With the fluid-fluid interface modeled to be diffuse, the contact line can thus move relative to the solid wall through diffusion rather than convection. All the above models are at least mathematically valid because the divergence of stress has been avoided, either by introducing some molecular process to drive the contact line [4], or by allowing slip to occur [19, 21, 43], or by modeling a diffuse fluid-fluid interface [7,24,38]. Pismen and Pomeau have presented a rational analysis of the hydrodynamic phase field (diffuse interface) model based on the

lubrication approximation [31]. The most usual (and natural) approach to resolve the stress divergence has been to allow slip at the solid wall close to the contact line. In fact, molecular dynamics (MD) studies have clearly demonstrated the existence of fluid slipping in the molecular-scale vicinity of the MCL [28, 40]. However, the exact rule that governs this relative slip has eluded numerous prior attempts. As a matter of fact, none of the existing models has proved successful by quantitatively accounting for the contact-line slip velocity profile observed in MD simulations. In a hybrid approach by Hadjiconstantinou [15], the MD slip profile was used as the boundary condition for finite-element continuum calculation. The continuum results so obtained match the corresponding MD results, therefore demonstrating the feasibility of hybrid algorithm [16, 36]. But the problem concerning the boundary condition governing the contact-line motion was still left unsolved. In this paper we

first present an introductory review of the problem, including (1) the 4 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 (a) (b) fluid 1 γ contact line γ1 fluid 2 θs γ 2 solid wall fluid 1 γ1 γ θd fluid 2 x solid wall γ2 U Figure 2: (a) Static contact line. A fluid-fluid interface is formed between two immiscible fluids and intersects the solid wall at a contact line. The static contact angle θs is determined by the Young’s equation: γ cos θs + γ2 − γ1 = 0, where γ1 , γ2 , and γ are the three interfacial tensions at the three interfaces (two fluid-solid and one fluid-fluid). (b) Moving contact line When one fluid displaces another immiscible fluid, the contact line is moving relative to the solid wall. (Here fluid 1 displaces fluid 2) Due to the contact-line motion, the dynamic contact angle θd deviates from the static contact angle θs . We will show that this deviation is primarily responsible for the near-complete slip at the

contact line. origin of the stress singularity, (2) the ad hoc introduction of the slip boundary condition, (3) the MD evidence of fluid slipping, (4) the gap between the existing MD results and a continuum hydrodynamic description, and (5) a preliminary account on how to bridge the MD results and a continuum description. We then present a detailed review of our recent results on the contact-line motion in immiscible two-phase flow, from MD simulations to continuum hydrodynamics calculations [34]. In our MD simulations, we consider two immiscible dense fluids of identical density and viscosity, with the temperature controlled above the gas-liquid critical point. (Similar results would be obtained if the temperature is below the critical point) For fluid-solid interactions, we choose the wall density and interaction parameters to make sure (1) there is no epitaxial locking of fluid layer(s) to the solid wall, (2) a finite amount of slip is allowed in the single-phase flow for each of

the two immiscible fluids, (3) the fluid-fluid interface makes a finite microscopic contact angle with the solid surface. Through extensive MD simulations and detailed analysis, we have uncovered the boundary condition governing the fluid slipping in the presence of a MCL. With the help of this discovery, we have formulated a continuum hydrodynamic model of two-phase immiscible flow. Numerical solutions have been obtained through an explicit finite-difference scheme. A comparison of the MD and continuum results shows that velocity and fluid-fluid interface profiles from the MD simulations can be quantitatively reproduced by the continuum model. The paper is organized as follows. In Section 2 we review a few known facts in mathematics and physics concerning the contact line motion. Together, they point out the right direction to approach and elucidate the problem. In fact, they almost tell us what is expected for a boundary condition governing the slip at the MCL and in its vicinity. In

Section 3 we outline the main results from MD simulations to continuum hydrodynamic modeling. In Section 4 we present Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 5 θ=φ r PHASE B PHASE A θ θ=π θ=0 U Figure 3: The Hua-Scriven model [20] of MCL. the first part of the MD results. From the slip velocity and the tangential wall force measured at the fluid-solid interface, a slip boundary condition is deduced. In Section 5 we formulate a continuum hydrodynamic model of two-phase immiscible flow. The continuum differential expression for the tangential stress at the solid surface is derived, from which the generalized Navier boundary condition (GNBC) is obtained from the slip boundary condition deduced in Section 4. In Section 6 we show a systematic comparison of the MD and continuum results The validity of the GNBC and the continuum model is demonstrated. In Section 7 we present the second part of the MD results, concerning the tangential force balance in a

boundary layer at the fluid-solid interface and the decomposition of the tangential stress in the fluid-fluid interfacial region. In Section 8 we establish the correspondence between the stress components measured in MD and those defined in the continuum hydrodynamics. Based on this correspondence, the continuum GNBC is obtained in an integrated form from the MD results in Sections 4 and 7. This may be regarded as a direct MD verification of the continuum GNBC. It also justifies the use of the Cahn-Hilliard hydrodynamic formulation of two-phase flow, from which the continuum form of GNBC, as verified by the MD results, is derived. The paper is concluded in Section 9 2 2.1 Stress and slip: A brief review Non-integrable stress singularity: the Huh-Scriven model Hua and Scriven [20] considered a simple model for the immiscible two-phase flow near a MCL. A flat solid surface is moving with steady velocity U in its own plane and a flat fluid-fluid interface, formed between two immiscible

phases A and B, intersects the solid surface at a contact line (see Fig. 3) The contact line is taken as the origin of a polar coordinate system (r, θ), in which the contact angle is φ. The two-dimensional flow of Newtonian and incompressible fluids is dominated by viscous stress for r ¿ η/ρU , where η is the viscosity and ρ the mass density. In the viscous flow approximation, the Navier-Stokes equation is linearized and the steady flow is solved from the biharmonic equation for the stream function ψ(r, θ): ∇2 (∇2 ψ) = 0. 6 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 The boundary conditions to be imposed at θ = 0, π (solid surface) and θ = φ (fluid-fluid interface) are: (i) vanishing normal component of velocity at the solid surface and fluid-fluid interface, (ii) continuity of tangential velocity at the fluid-fluid interface, (iii) continuity of tangential viscous stress at the fluid-fluid interface, (iv) no slip at the solid surface. Here

conditions (i)–(iii) are well justified while (iv) is usually taken as a postulate of continuum hydrodynamics. The no-slip boundary condition can be justified a posteriori in macroscopic flow by checking the correctness of its consequences. In the present model, however, it leads to physically incorrect results for stress, in the microscopic vicinity of the MCL. The similarity solution of the biharmonic equation is in the form of ψ(r, θ) = r(a sin θ + b cos θ + cθ sin θ + dθ cos θ), in which the eight constants (4 for phase A and 4 for phase B) are determined by the eight boundary conditions in (i)–(iv). What Hua and Scriven found is that the shear stress and pressure fields vary as 1/r and hence increase in magnitude without limit as the contact line r = 0 is approached. As a consequence, the total tangential force exerted on the solid surface, which is an integral of the tangential stress along the surface, is logarithmically infinite. That indicates a singularity in

viscous dissipation, which is physically unacceptable. Obviously, the Hua-Scriven model is defective in the immediate vicinity of the MCL. This is also seen from the normal stress difference across the fluid-fluid interface, which varies as 1/r as well. According to the Laplace’s equation, the interface curvature should increase rapidly as the contact line (r = 0) is approached, in order to balance the normal stress difference by curvature force. This is clearly inconsistent with the assumption of a flat fluid-fluid interface Nevertheless, the flow field solved from the Hua-Scriven model may approximately describe the asymptotic region at a large distance from the contact line (where the viscous flow approximation is still valid). In that region, the no-slip boundary condition is considered valid and the fluid-fluid interface is almost flat, due to the reduced stress. 2.2 Introducing the slip boundary condition The deficiency of the Hua-Scriven model results from the

incompatibility of the no-slip boundary condition with a MCL: no slip means vr = ±U at the solid surface (θ = 0, π) where r > 0 while at r = 0 the MCL requires a perfect slip. That is, the no-slip boundary condition leads to a velocity discontinuity at the MCL. In order to remove the stress singularity at the MCL, while retaining the Newtonian behavior of stress, the continuity of velocity field must be restored. For this purpose, a slip profile can be introduced to continuously interpolate between the complete slip at the MCL and the no-slip boundary condition that must hold at regions far away. Dussan V. [10] considered a plate of infinite extent either inserted into or withdrawn from a semi-infinite domain of fluid at a constant velocity U0 (see Fig. 4) The contact line is taken as the origin of a polar coordinate system (r, φ), in which the apparent contact angle is α at r ∞. The equation of motion is the Navier-Stokes equation with the incompressibility condition ∇ ·

u = 0. The boundary conditions at the solid surface φ = 0 and the free surface {(r, φ(r)}|0 ≤ r < ∞} are: (i) the kinematic boundary conditions u · φ̂ = 0 at φ = 0 and Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 k 7 Free surface α φ(φ) r r (φ) φ U0 Fluid φ=0 Figure 4: A plate is inserted into a semi-infinite domain of fluid at a constant velocity U0 . The angle between the plate and the free surface at r ∞ is α. u · n = 0 at φ = φ(r), where n is an outward unit vector normal to the free surface, (ii) the dynamic boundary condition t · T · n = 0 at φ = φ(r), where T is the stress tensor and t is a unit vector tangent to the free surface, i.e, t · n = 0, (iii) the Laplace condition n · T · n = σκ at φ = φ(r), where σ is the surface tension and κ the interface curvature, (iv) the slip boundary condition u = U (r)r̂ at φ = 0, where U (r) is a prescribed function, and (v) the critical static contact angle φ(0) = φs

. The slip boundary condition must continuously interpolate between the complete slip at the MCL (r = 0) and the no-slip boundary condition at r ∞: lim U (r) = 0, r0 lim U (r) = U0 . r∞ However, the form of U (r) in the intermediate region, where U varies from 0 to U0 , is unknown. Three different slip boundary conditions were used for U (r), in order to assess the sensitivity of the overall flow field to the form of the slip boundary condition. They are U1 = r/Ls (r/Ls )2 U0 , U2 = U0 , 1 + r/Ls 1 + (r/Ls )2 U1 = (r/Ls )1/2 U0 . 1 + (r/Ls )1/2 where Ls is a length scale called the slip length (not the one defined in the Navier boundary condition). It was found that on the slip length scale the flow fields are quite different whereas on the meniscus length scale, i.e, the length scale on which almost all fluid-mechanical measurements are performed, all the flow fields are virtually the same That is, identical macroscopic flow behaviors are expected from different slipping

models. 2.3 Slip observed in molecular dynamics simulations According to the conclusion in [10], only in the microscopic slip region (of length scale ∼ Ls ) can different slip models be distinguished. This makes the experimental verification of a particular 8 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 U θapp Phase A Phase B Phase A z U x near complete slip near complete slip small partial slip small partial slip Figure 5: The xz projection of the simulated system. Thick solid lines represent solid walls and dashed lines represent fluid-fluid interfaces. The arrows indicate fluid velocity close to the solid walls, from which the variation of the amount of slip is clearly seen. slip model very difficult because experiments usually probe distances (∼ 1µm) much larger than Ls . Naturally, computer experiments become very useful in elucidating the problem [1] Non-equilibrium MD simulations were carried out to investigate the fluid motion in the

vicinity of the MCL, in both the Poiseuille-flow [28] and Couette-flow [40] geometries. In these MD simulations, interactions between fluid molecules were modeled with the Lennard-Jones potential, modified to segregate immiscible fluids. The confining walls were constructed with a molecular structure. Wall-fluid interactions were also modeled with the Lennard-Jones potential, with energy and length scales different from those of the fluid-fluid interactions In the simulations performed in the Couette geometry [40], two immiscible fluids were confined between two planar walls parallel to the xy plane and a shear flow was induced by moving the two solid walls in ±x directions at constant speed U (see Fig. 5) Steady-state velocity fields were obtained from the time average of fluid molecular velocities in small bins. There was clean evidence for a slip region in the vicinity of the MCL: appreciable slip occurs within a length scale ∼ 10σ, where σ is the length scale in the

fluid-fluid Lennard-Jones potential, and at the MCL the slip is near-complete. Far from the interfaces, viscous damping makes the flow insensitive to the fast variations near the interfaces, and hence uniform shear flow was observed, with a negligible amount of slip. Therefore, MD simulations provide evidence for a cutoff, below which the no-slip boundary condition breaks down, introduced phenomenologically in various slip models to remove the stress singularity. However, the exact boundary condition that governs the observed slip was left unresolved. In particular, a breakdown of local hydrodynamics in the molecular scale slip region was suggested [40], considering the extreme velocity variations there. The Navier slip model assumes that the amount of slip is proportional to the tangential component of the stress tensor, Pxz , in the fluid at the solid surface. In Ref [40], the microscopic value of Pxz was directly measured A comparison to the slip profile is roughly consistent with

the Navier slip model. However, a large discrepancy was found between the microscopic Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 9 L S/1 fluid 1 γ12 γs1 v θ r fluid 2 S/2 γs2 S Figure 6: A displacement of fluid 2 by fluid 1 over distance L in a cylindrical capillary of radius r. value of Pxz and the shear rate ∂vx /∂z. According to the authors, this discrepancy suggests a breakdown of local hydrodynamics. On the contrary, we will show in this paper: (i) The Navier slip model is the correct model describing the fluid slipping near the MCL. (ii) The tangential stress in the Navier model is not merely viscous. (iii) The interfacial tension plays an important role in governing the contact-line slip in immiscible two-phase flows. (iv) There is no breakdown of local hydrodynamics. 2.4 2.41 Fluid/fluid displacement driven by unbalanced Young force Unbalanced Young force For a cylindrical capillary of radius r, if the steady displacement is

sufficiently slow, then the pressure drop across the moving fluid-fluid interface is given by ∆p = 2γ12 cos θ/r, where γ12 is the interfacial tension, and θ is an appropriate dynamic contact angle. At equilibrium, the pressure drop is given by ∆p0 = 2γ12 cos θ0 /r, where θ0 is the equilibrium contact angle. In general, θ and θ0 differ. Physically, ∆p0 measures the change of the interfacial free energy at the fluid-solid interface when the fluid-fluid interface moves relative to the wall, and ∆p measures the external work done to the system. Therefore, the difference ∆p − ∆p0 is a measure of the dissipation due to the presence of the MCL. Let’s consider a displacement of fluid 2 by fluid 1 over distance L (see Fig. 6) According to the Young’s equation for the equilibrium contact angle, the force πr2 ∆p0 equals 2πrγ12 cos θ0 = 2πr(γS1 − γS2 ), where γS1 and γS2 are the interfacial tensions for the S/1 and S/2 interfaces, respectively. Thus the

change of the interfacial free energy at the fluid-solid interface is given by πr2 ∆p0 L = 2πrL(γS1 − γS2 ). The external work done to the system is simply πr2 ∆pL It follows that the dissipation due to the presence of the MCL is given by πr2 (∆p − ∆p0 )L = 2πrLγ12 (cos θ − cos θ0 ), where γ12 (cos θ − cos θ0 ) = γ12 cos θ + γS2 − γS1 is the unbalanced Young force [12]. In order to find a relation between the displacement velocity v and the unbalanced Young force FY , two classes of models have been proposed to describe the contact line motion: a) An Eyring approach, based on the molecular adsorption/desorption processes at the contact 10 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 ξ (a) fluid 1 v θ fluid 2 (b) U K+ K− W ~FY λ2 λ λ x Figure 7: (a) A molecular picture of the three-phase zone. (b) Shifted potential profile line [4]. b) A hydrodynamic approach, assuming that dissipation is dominated by viscous

shear flow inside the wedge [12]. Viscous flows in wedges were calculated by Hua and Scriven [20] For wedges of small (apparent) contact angle, a lubrication approximation can be used to simplify the calculations [12]. As discussed in Sections 21 and 22, a (molecular scale) cutoff has to be introduced to remove the logarithmic singularity in viscous dissipation. On the contrary, the Eyring approach assumes that the molecular dissipation at the tip is dominant. A quantitative theory is briefly reviewed below. 2.42 Blake’s kinetic model The role of interfacial tension was investigated in a kinetic model by Blake and Haynes [4]. Consider a fluid-fluid interface in contact with a flat solid surface at a line of three-phase contact (see Fig. 7a) Viewed on a molecular scale, the three-phase line is actually a fluctuating threephase zone, where adsorbed molecules of one fluid (at the solid surface) interchange with those of the other, either by migration at the solid surface or through

the contiguous bulk phases. In equilibrium the net rate of exchange will be zero. For a three-phase zone moving relative to the solid wall (Fig. 7a), the net displacement, driven by the unbalanced Young force FY = γ12 cos θ + γS2 − γS1 , is due to a nonzero net rate of exchange. Let ξ be the interfacial thickness, σ be the area of an adsorption site, and λ be the hopping distance of molecules. The force per adsorption site is approximately σFY /ξ and the energy shift over distance λ is approximately λσFY /ξ ∼ FY λ2 (see Fig. 7b) This energy shift leads to two different rates K+ and K− : ¶¸ · µ 1 1 2 , K± = k exp − W ∓ FY λ kB T 2 for forward and backward hopping events, respectively. Here W is an activation energy ¡ ¢ for hopping. It follows that the velocity of the MCL is v = λ(K+ − K− ) ∝ sinh FY λ2 /2kB T For very small FY , v ∝ FY . Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 z=z0 11 z z=0 x Figure 8: A boundary

layer of fluid at the fluid-solid interface, responsible for viscous momentum transport between fluid and solid. Circles represent fluid molecules and solid squares represent wall molecules Blake’s kinetic model shows that fluid slipping can be induced by the unbalanced Young force at the contact line. Therefore, it emphasizes the role of interfacial tension, though not in a hydrodynamic formulation. On the contrary, the authors of Ref [40] considered the viscous shear stress as the only driving force. In fact a large discrepancy was found between the shear rate and the microscopic value of the tangential stress (the driving force according to the Navier slip model). In this paper we will show that in the two-phase interfacial region, such a discrepancy is exactly caused by the neglect of the non-viscous contribution from interfacial tension. 2.5 From the Navier boundary condition to the generalized Navier boundary condition: a preliminary discussion Here we give a preliminary

account on the main finding reported in Ref. [34], the GNBC Based on the results reviewed in Sections 2.1, 22, 23, and 24, we try to show: a) The Navier boundary condition may actually work for the single-phase slip region near the MCL. b) In the two-phase contact-line region, the GNBC is a natural extension of the Navier boundary condition, with the fluid-fluid interfacial tension taken into account. 2.51 Navier boundary condition and slip length £ ¤ The validity of the Navier boundary condition vslip = −ls n · (∇v) + (∇v)T has been well established by many MD studies for single-phase fluids [2, 8, 39, 41]. This boundary condition is a constitutive equation that locally relates the amount of slip to the shear rate at the solid surface, though in most of the reported simulations, the hydrodynamics involves no spatial variation along the solid surface. Physically, the Navier boundary condition is local in nature simply because intermolecular wall-fluid interactions are

short-ranged. The slip length ls is a phenomenological parameter that measures the local viscous coupling between fluid and solid. Fig 8 illustrates the viscous momentum transport between fluid and solid through a boundary layer of fluid. The thickness of this boundary layer, z0 , must be of molecular scale, within the range of wall-fluid interactions. Now we show that the slip length ls is defined based on a linear law for tangential wall force and the Newton’s law for shear stress. Hydrodynamic motion of fluid at the solid surface is measured by the slip velocity vxslip in the x direction, defined relative to the (moving) wall. When vxslip is present, a tangential wall 12 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 force Gw x is exerted on the boundary-layer fluid, defined as the rate of tangential (x) momentum transport per unit wall area, from the wall to the (boundary-layer) fluid. Physically, this force represents a time-averaged effect of wall-fluid

interactions, to be incorporated into a hydrodyw slip , where β is namic slip boundary condition. The linear law for Gw x is expressed by Gx = −βv called the slip coefficient and the minus sign means the fluid-solid coupling is viscous (frictional). w For the boundary layer of molecular R z0thickness, the tangential wall force Gx must be balanced f by the tangential fluid force Gx = 0 dz (∂x σxx + ∂z σzx ), where the integral is over the z-span of the boundary layer (see Fig. 8), σxx(zx) are the xx(zx) components of the fluid stress tensor, and ∂x σxx + ∂z σzx is the fluid force density in the x direction. From the equation of tangential f force balance Gw x + Gx = 0, we obtain Z z0 Z z0 dzσxx (x, z) + σzx (x, z0 ) = βv slip (x). dz (∂x σxx + ∂z σzx ) = ∂x 0 0 This equation should be regarded as a microscopic expression for the Navier slip model: the amount of slip is proportional to the tangential fluid force at the fluid-solid interface. When the normal

stress σxx varies slowly in the x direction and the tangential stress σzx is caused by shear viscosity only, the above equation becomes η∂z vx (x, z0 ) = βv slip (x), where η is the viscosity and ∂z vx (x, z0 ) is the shear rate “at the solid surface”. For flow fields whose characteristic length scale is much larger than z0 , the boundary layer thickness, we replace z0 by z = 0 and recover the Navier boundary condition for continuum hydrodynamics, in which the slip length is given by ls = η/β. To summarize, the hydrodynamic viscous coupling between fluid and solid is actually measured by the slip coefficient β. The Navier boundary condition, in which the slip length is introduced, is due to the Newton’s law for viscous shear stress. For very weak viscous coupling between fluid and solid, β 0, and thus ls ∞, making ∂z vx 0: the fluid-solid interface becomes a free surface. 2.52 Single-phase region A number of MD studies have shown that for single-phase flows,

the Navier boundary condition is valid in describing the fluid slipping at solid surface [2, 8, 39, 41]. Therefore, we expect that it can also describe the partial slip observed in the single-phase region near the MCL [40]. However, according to the authors of Ref. [40], the Navier boundary condition failed even in the single-phase region: the velocity gradient ∂vx /∂z was not proportional to the amount of slip. This discrepancy was attributed to a breakdown of local hydrodynamics, considering the very large velocity variations observed near the MCL. Here we present a heuristic discussion, to explain why such a discrepancy is expected even if the Navier boundary condition actually works for the single-phase region. First, the success of the hybrid approaches outlined below strongly indicates that local hydrodynamics doesn’t break down in the slip region. In Ref [40], an apparent contact angle θapp was defined at half the distance between two solid surfaces (see Fig. 5) For θapp

< 135◦ , fluid-fluid interfaces could be approximated by planes. The Navier-Stokes equation was solved for the simplified geometry. For this purpose, the tangential velocity along fluid-fluid interface Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 13 z θ app =π/2 no−slip condition Navier−Stokes equation Navier condition vx O U x Figure 9: Two-dimensional corner flow caused by one rigid plane sliding over another, with constant inclination θapp = 90◦ . In the reference frame of the vertical plane, the horizontal plane is moving to the left with constant speed U . was set to be zero (according to the simulation results) and the slip velocity ∆V (r) (measured relative to the moving wall) was specified as a function of distance r from the MCL: ∆V (r) = ±U exp[−r ln 2/S], with + and − for the lower and upper walls, respectively. This ∆V (r), proposed according to the MD slip profile, continuously interpolates between the complete slip

at the MCL (∆V (r) = ±U at r = 0) and the no-slip boundary condition (∆V (r) 0 at r ∞). With a proper value chosen for S (≈ 18σ), the MD flow fields were reproduced by solving the Navier-Stokes equation with the above boundary conditions. Recently, an improved hybrid approach has been used to reproduce the MD simulation results for contact-line motion in a Poiseuille geometry [15]. To further test the validity of continuum modeling, we have solved the Navier-Stokes equation for a corner flow [33]. Consider one rigid plane sliding steadily over another, with constant inclination θapp (see Fig. 9) The fluid is Newtonian and incompressible The no-slip boundary condition is applied at the vertical plane and the Navier boundary condition applied at the horizontal plane. The kinematic condition of vanishing normal component of velocity at the solid surface requires vx = 0 at the vertical plane and vz = 0 at the horizontal plane, and hence v = 0 at the intersection O, which is

taken as the origin of the coordinate system (x, z). This corner-flow model resembles the continuum model used in the hybrid approach above, and produces similar flow fields for quantitative comparison with MD results. In particular, the corner flow exhibits a slip profile very close to ∆V (r): the slip velocity vx + U shows a linear decrease over a length scale ∼ ls , the slip length in the Navier condition, followed by a more gradual decrease. (Note that vx + U = U at O implies complete slip) From the hybrid approach with the prescribed slip velocity ∆V (r) to the corner-flow model with the Navier boundary condition, they indicate that the single-phase region near the MCL can be modeled by the Navier-Stokes equation for an incompressible Newtonian fluid, supplemented by appropriate boundary conditions. Then a simple question arises: Given a continuum hydrodynamic model that uses the Navier boundary condition and approximately reproduces 14 Qian, Wang and Sheng / Commun.

Comput Phys, 1 (2006), pp 1-52 the MD flow fields in the single-phase region near the MCL, why do the MD simulation results show a large discrepancy between the velocity gradient ∂vx /∂z and the amount of slip? The answer lies in the fast velocity variation in the vicinity of the MCL, where the flow is dominated by viscous stress. With the characteristic velocity scale set by U and the characteristic length scale set by ls , the normal stress σxx must show as well a fast variation along the solid surface: ∂x σxx ∼ ηU/ls2 . According to Section 251, the microscopic value of the tangential Rz fluid force Gfx is given by ∂x 0 0 dzσxx (x, z) + σzx (x, z0 ). This expression is necessary because Gfx is distributed in a boundary layer of molecular thickness z0 . Obviously, to represent Gfx by σzx ≈ η∂vx /∂z at z = z0 alone, the normal stress σxx near the solid surface has to vary slowly in the x direction. This is not the case in the vicinity of the MCL if the Navier

slip model works there, as implied by the continuum corner-flow calculation which yields ∂x σxx ∼ ηU/ls2 near the intersection. (It is reasonable to expect that if the Navier slip model is valid, then the normal stress measured in MD simulations should in general from continuum R z0 agree with that 2 model calculation with the Navier condition.) Given ∂x 0 dzσxx ∼ ηU z0 /ls according to the corner-flow model, and that z0 and ls are both ∼ 1σ, it is obvious that considering only η∂vx /∂z at z = z0 would lead to an appreciable underestimate of Gfx ∼ ηU/ls in the slip region. In short, the large discrepancy between the velocity gradient ∂vx /∂z and the amount of slip, observed in the single-phase slip region near the MCL, cannot be used to exclude the microscopic Navier slip model and the hydrodynamic Navier slip boundary condition, for if the Navier model is valid, then the tangential viscous stress η∂vx /∂z, measured at some level away from the solid

surface, is not enough for a complete evaluation of the tangential fluid force Gfx . 2.53 Two-phase region slip The linear law for tangential wall force, Gw x = −βvx , describes the hydrodynamic viscous coupling between fluid and solid. Assume that in the two-phase region, the two fluids interact slip with the wall independently. Then the tangential wall force becomes Gw at the x = −βvx contact line, with the slip coefficient given by β = (β1 ρ1 + β2 ρ2 )/(ρ1 + ρ2 ), where β1 and β2 are the slip coefficients for the two single-phase regions separated by the fluid-fluid interface, and ρ1 and ρ2 are the local densities of the two fluids in the contact-line region. Obviously, β varies between β1 and β2 across the fluid-fluid interface. f The equation of tangential force balance Gw x + Gx = 0 must hold as well in the two-phase region. Therefore the Navier slip model is still of the form Z z0 ∂x dzσxx (x, z) + σzx (x, z0 ) = βv slip (x). 0 Nevertheless, this does

not lead to the Navier boundary condition η∂z vx = βv slip anymore because in the two-phase region, the tangential stress is not contributed by shear viscosity only: there is a non-viscous component in σzx , caused by the fluid-fluid interfacial tension. Put in an ideal form of decomposition, the tangential stress σzx (x, z) at level z can be expressed as Y v (x, z), (x, z) + σzx σzx (x, z) = σzx Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 15 v is the viscous component due to shear viscosity and σ Y (the tangential Young stress) is where σzx zx the non-viscous component, which is narrowly distributed in the two-phase interfacial region and R R Y related to the interfacial tension through int dxσzx (x, z) = γ cos θ(z). Here int dx denotes the integration across the fluid-fluid interface along the x direction, γ is the interfacial tension, and θ(z) is the angle between the interface and the x direction at level z. Physically, the existence of a

fluid-fluid interface causes an anisotropy in the stress tensor in the two-phase interfacial region. The interfacial tension is an integrated measure of that stress anisotropy. When expressed in a coordinate system different from the principal system, the stress tensor is not diagonalized and a Y appears. In the presence of shear flow, while the shear viscosity leads to the viscous nonzero σzx v in σ , the non-viscous component σ Y is still present. A detailed discussion on component σzx zx zx the stress decomposition will be given in Section 7.3 0 is balanced by the gradient Consider an equilibrium state in which the tangential stress σzx 0 : of the normal stress σxx Z z0 ∂x 0 0 0 dzσxx (x, z) + σzx (x, z0 ) = 0, where the superscript 0 denotes equilibrium quantities. Here the equilibrium tangential stress 0 is narrowly distributed in the two-phase interfacial region and related to the interfacial σzx R 0 (x, z) = γ cos θ 0 (z). (There is no viscous stress in

equilibrium, and tension through int dxσzx 0 is caused by the interfacial tension only.) Subtracting the equation of equilibrium hence σzx force balance from the expression of Navier slip model, we obtain Z z0 0 0 dz[σxx (x, z) − σxx (x, z)] + [σzx (x, z0 ) − σzx (x, z0 )] ∂x 0 Z z0 0 v Y 0 dz[σxx (x, z) − σxx (x, z)] + [σzx (x, z0 ) + σzx (x, z0 ) − σzx (x, z0 )] = ∂x 0 = βv slip (x). This equation will be the focus of our continuum deduction from molecular hydrodynamics. In fact it leads to the GNBC which governs the fluid slipping everywhere, from the two-phase contact-line region to the single-phase regions away from the MCL. This will be elaborated in Sections 4, 5, 7, and 8. To summarize, our preliminary analysis shows that compared to the single-phase region where the tangential stress is due to shear viscosity only, the two-phase region has the tangential Young stress as the additional component. Naturally, the Navier boundary condition, which

considers the tangential viscous stress only, needs to be generalized to include the tangential Young stress. 3 Statement of results We have carried out MD simulations for immiscible two-phase flows in a Couette geometry (see Fig. 10) [34] The two immiscible fluids were modeled by using the Lennard-Jones (LJ) potentials for the interactions between fluid molecules. The solid walls were modeled by crystalline plates 16 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 V H z V y x Figure 10: Schematic of the immiscible Couette flow. Two immiscible fluids (empty and solid circles) are confined between two solid walls (solid squares) that are parallel to the xy plane and separated by a distance H. The Couette flow is generated by moving the top and bottom walls at a speed V along ±x, respectively. More technical details will be presented in Sections 4 and 7. The purpose of carrying out MD simulations is threefold: (1) To uncover the boundary condition governing

the MCL, denoted the GNBC; (2) To fix the material parameters (e.g viscosity, interfacial tension, etc) in our hydrodynamic model; (3) To produce velocity and interface profiles for comparison with the continuum hydrodynamic solutions. Our main finding is the GNBC: v Y βvxslip = σ̃zx = σzx + σ̃zx . Here β is the slip coefficient, vxslip is the tangential slip velocity measured relative to the (moving) wall, and σ̃zx is the hydrodynamic tangential stress, given by the sum of the viscous stress v and the uncompensated Young stress σ̃ Y . The validity of the GNBC has been verified by a σzx zx detailed analysis of MD data (Sections 4, 7, and 8) plus a comparison of the MD and continuum results (Section 6). Compared to the conventional Navier boundary condition which includes the tangential viscous stress only, the GNBC captures the uncompensated Young stress as the additional component. Together, the tangential viscous stress and the uncompensated Young stress give rise to the

near-complete slip at the MCL. The uncompensated Young stress arises from the deviation of the fluid-fluid interface from its static configuration and is narrowly distributed in the fluid-fluid interfacial region. Obviously, far away from the MCL, the uncompensated Young v . stress vanishes and the GNBC becomes the usual Navier Boundary condition βvxslip = σzx We have incorporated the GNBC into the Cahn-Hilliard (CH) hydrodynamic formulation of two-phase flow [7, 24] to obtain a continuum hydrodynamic model [34]. The continuum model may be h briefly describedi as follows. The CH free energy functional [5] is of the form R F [φ] = dr K (∇φ)2 /2 + f (φ) , where φ(r) is the composition field defined by φ(r) = (ρ2 − ρ1 )/(ρ2 + ρ1 ), with ρ1 and ρ2 being the local number densities of the two fluid species, f (φ) = −rφ2 /2 + uφ4 /4, and K, r, u are parameters which can be determined in MD simulations by Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52

17 p √ 2 measuring the interface width ξ = K/r, pthe interfacial tension γ = 2 2r ξ/3u, and the two homogeneous equilibrium phases φ± = ± r/u (= ±1 in our case). The two coupled equations of motion are the CH convection-diffusion equation for φ and the Navier-Stokes equation (with the addition of the capillary force density): ∂φ + v · ∇φ = M ∇2 µ, ∂t ¸ · ∂v + (v · ∇) v = −∇p + ∇ · σ v + µ∇φ + mρgext , ρm ∂t (3.1) (3.2) together with the incompressibility condition ∇·v = 0. Here M is the phenomenological mobility coefficient, ρm is the fluid mass density, p is the pressure, σ v is the Newtonian viscous stress tensor, µ∇φ is the capillary force density with µ = δF/δφ being the chemical potential, and mρgext is the external body force density (for Poiseuille flows). The boundary conditions at the solid surface are vn = 0, ∂n µ = 0 (n denotes the outward surface normal), a relaxational equation for φ at the solid surface:

∂φ + v · ∇φ = −ΓL(φ), ∂t (3.3) and the continuum form of the GNBC: βvxslip = −η∂n vx + L(φ)∂x φ, (3.4) Here Γ is a (positive) phenomenological parameter, L(φ) = K∂n φ + ∂γwf (φ)/∂φ with γwf (φ) being the fluid-solid interfacial free energy density, β is the slip coefficient, vxslip is the (tangential) slip velocity relative to the wall, η is the viscosity, and L(φ)∂x φ is the uncompensated Young stress. Compared to the Navier boundary condition, the additional component captured by the Y . Its differential expression is given GNBC in Eq. (34) is the uncompensated Young stress σ̃zx by Y σ̃zx = L(φ)∂x φ = [K∂n φ + ∂γwf (φ)/∂φ]∂x φ (3.5) at z = 0. Here the z coordinate is for the lower fluid-solid interface where ∂n = −∂z (same below), with the understanding that the same physics holds at the upper interface. It can be shown that the integral of this uncompensated Young stress along x across the fluid-fluid

interface yields Z int dx[L(φ)∂x φ](x, 0) = γ cos θdsurf + ∆γwf , (3.6) R where int dx denotes the integration along x across the fluid-fluid interface, γ is the fluid-fluid surface, and ∆γwf is the interfacial tension, θdsurf is the dynamic contact angle at the solid R change of γwf (φ) across the fluid-fluid interface, i.e, ∆γwf ≡ int dx∂x γwf (φ) The Young’s equation relates ∆γwf to the static contact angle θssurf at the solid surface: γ cos θssurf + ∆γwf = 0, (3.7) 18 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 which is obtained as a phenomenological expression for the tangential force balance at the contact line in equilibrium. Substituting Eq (37) into Eq (36) yields Z dx[L(φ)∂x φ](x, 0) = γ(cos θdsurf − cos θssurf ), (3.8) int implying that the uncompensated Young stress arises from the deviation of the fluid-fluid interface from its static configuration. Equations (35), (36), (37), and (38) will be

derived in Section 5. In essence, our results show that in the vicinity of the MCL, the tangential viscous stress −η∂n vx as postulated by the usual Navier boundary condition can not account for the contact-line slip profile without taking into account the uncompensated Young stress. This is seen from both the MD data and the predictions of the continuum model. Besides the external conditions such as the shear speed V and the wall separation H, there are nine parameters in our continuum model, including ρm , η, β, ξ, γ, |φ± |, M , Γ, and θssurf . The values of ρm , η, β, ξ, γ, |φ± |, and θssurf were directly obtained from MD simulations. The values of the two phenomenological parameters M and Γ were fixed by an optimized comparison with one MD flow field. The same set of parameters (corresponding to the same local properties in a series of MD simulations) has been used to produce velocity fields and fluid-fluid interface shapes for comparison with the MD results

obtained for different external conditions (V , H, and flow geometry). The overall agreement is excellent in all cases, demonstrating the validity of the GNBC and the hydrodynamic model. The CH hydrodynamic formulation of immiscible two-phase flow has been successfully used to construct a continuum model. We would like to emphasize that while the phase-field formulation does provide a convenient way of modeling that is familiar to us, it should not be conceived as the unique way. After all, what we need is to incorporate our key finding, the GNBC, into a continuum formulation of immiscible two-phase flow. The GNBC itself is simply a fact found in MD simulations, independent of any continuum formulation. 4 4.1 Molecular dynamics I Geometry and interactions MD simulations have been carried out for two-phase immiscible flows in Couette geometry (see Figs. 10 and 11) [34] Two immiscible fluids were confined between two planar solid walls parallel to the xy plane, with the fluid-solid

boundaries defined by z = 0, H. The Couette flow was generated by moving the top and bottom walls at a constant speed V in the ±x directions, respectively. Periodic boundary conditions were imposed along the x and y directions Interaction between fluid molecules separated by a distance r was modeled by a modified LJ potential i h Uf f = 4² (σ/r)12 − δf f (σ/r)6 , where ² is the energy scale, σ is the range scale, with δf f = 1 for like molecules and δf f = −1 for molecules of different species. (The negative δf f was used to ensure immiscibility) The Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 19 average number density for the fluids was set at ρ = 0.81σ −3 The temperature was controlled at 2.8²/kB (This high temperature was used to reduce the near-surface layering induced by the solid wall.) Each wall was constructed by two [001] planes of an fcc lattice, with each wall molecule attached to a lattice site by a harmonic spring. The mass of

the wall molecule was set equal to that of the fluid molecule m. The number density of the wall was set at ρw = 186σ −3 The wall-fluid interaction was modeled by another LJ potential h i Uwf = 4²wf (σwf /r)12 − δwf (σwf /r)6 , with the energy and range parameters given by ²wf = 1.16² and σwf = 104σ, and δwf for specifying the wetting property of the fluid. There is no locked layer of fluid molecules at the solid surface. We have considered two cases In the symmetric case, the two fluids have the identical wall-fluid interactions with δwf = 1. Consequently, the static contact angle is 90◦ and the fluid-fluid interface is flat, parallel to the yz plane. In the asymmetric case, the two fluids have different wall-fluid interactions, with δwf = 1 for one and δwf = 0.7 for the other As a consequence, the static contact angle is 64◦ and the fluid-fluid interface ispcurved in the xz plane. In most of our simulations, the shearing speed V was on the order of 01 ²/m, the

sample dimension along y was 6.8σ, the wall separation along z varied from H = 68σ to 68σ, and the sample dimension along x was set to be long enough so that the uniform single-fluid shear flow was recovered far away from the MCL. Steady-state quantities were obtained from time average p over 105 τ or longer where τ is the atomic time scale mσ 2 /². Throughout the remainder of this paper, all physical quantities are given in terms of the LJ reduced units (defined in terms of ², σ, and m). 4.2 Boundary-layer tangential wall force We denote the region within z0 = 0.85σ from the wall the boundary layer (BL) (see Fig 12) It must be thin enough to ensure sufficient precision for measuring the slip velocity at the solid surface, but also thick enough to fully account for the tangential wall-fluid interaction force, which is of a finite range. The wall force can be singled out by separating the force on each fluid molecule into wall-fluid and fluid-fluid components. The fluid

molecules in the BL, being close to the solid wall, can experience a strong periodic modulation in interaction with the wall. This lateral inhomogeneity is generally referred to as the “roughness” of the wall potential [41]. When coupled with kinetic collisions with the wall molecules, there arises a nonzero tangential wall force density gxw that is sharply peaked at z ≈ z0 /2 and vanishes beyond z ≈ z0 (see Fig. R z0 13).w From w w the force density gx , we define the tangential wall force per unit area as Gx (x) = 0 dzgx (x, z), which is the total tangential wall force accumulated across the BL. The short saturation range of the tangential wall force may be understood as follows. Out of the BL, each fluid molecule can interact with many wall molecules on a nearly equal basis. Thus the modulation amplitude (the roughness) of the wall potential would clearly decrease with increasing distance from the wall. That’s why the tangential wall force tends to saturate at z ≈ z0 ,

which is still within the cutoff distance of wall-fluid interaction. On the contrary, the normal wall force arises from the direct wall-fluid interaction, independent of whether the wall 20 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 (a) fluid 1 fluid 2 fluid 1 symmetric (b) V fluid 1 fluid 1 fluid 2 fluid 1 asymmetric V fluid 2 fluid 1 z y x Figure 11: Schematic of simulation geometry. (a) Static configurations in the symmetric and asymmetric cases. Here fluid 2 is sandwiched between fluid 1 due to the periodic boundary condition along the x direction. In the symmetric case, the static contact angle is 90◦ and the fluid-fluid interface is flat, parallel to the yz plane. In the asymmetric case, the static contact angle is not 90◦ and the fluid-fluid interface is curved in the xz plane. (b) Dynamic configuration in the symmetric case potential is rough or not. Consequently, the normal wall force saturates much slower than the tangential

component. We have measured the slip velocity and the tangential wall force in the BL. Spatial resolution along the x direction was achieved by evenly dividing the BL into bins, each ∆x = 0.85σ or 0.425σ The slip velocity vxslip was obtained as the time average of fluid molecules’ velocities inside the BL, measured relative to the moving wall; the tangential wall force Gw x was obtained from the time average of the total tangential wall force experienced by the fluid molecules in the BL, divided by the bin area in the xy plane. As reference quantities, we also measured Gw0 x in the slip w static (V = 0) configuration. Fig 14 shows vx and Gx measured in the dynamic configuration and Gw0 x measured in the static configuration. It is seen that in the absence of hydrodynamic motion (V = 0), the static tangential wall force Gw0 x is not identically zero everywhere. Instead, it has some fine features in the contact-line region (a few σ’s) (see Fig. 14b) This nonzero static component

in the tangential wall force arises from the microscopic organization of fluid molecules in the contact-line region. The static component is also present in Gw x measured in the dynamic configuration, as shown by Fig. 14b To see the effects arising purely from the hydrodynamic motion of the fluids, we subtract Gxw0 from Gw x through the relation w w0 G̃w x = Gx − Gx , w where G̃w x is the hydrodynamic part in Gx . In the notations below, the over tilde will denote the difference between that quantity and its static part. We find the hydrodynamic tangential Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 21 z z=z0 x Figure 12: Boundary layer at the lower fluid-solid interface. The empty and solid circles indicate the instantaneous molecular positions of the two fluids projected onto the xz plane. The solid squares denote the wall molecules. The dashed line indicates the level of z = z0 slip wall force per unit area, G̃w x , is proportional to the local slip

velocity vx : slip G̃w x (x) = −βvx (x), (4.1) where the proportionality constant β is the slip coefficient. In Fig 15, G̃w x is plotted as a function slip measured in the BL. The lines represent and v of vxslip . The symbols represent the values of G̃w x x slip slip w the values of G̃x calculated from −βvx using vx measured in the BL and β = (β1 ρ1 + β2 ρ2 )/(ρ1 + ρ2 ), with β1,2 the slip coefficients for the two fluid species and ρ1,2 the molecular densities of the two fluid species measured in the BL. Independent measurements determined √ √ √ β1 = β2 = 1.2 ²m/σ 3 for the symmetric case, β1 = 12 ²m/σ 3 and β2 = 0532 ²m/σ 3 for the asymmetric case. (The dependence of β on β1,2 and ρ1,2 assumes the two fluids interact with the wall independently. The desired expression is obtained by expressing G̃w x as the weighted slip1 w2 = −β v slip2 and noting v slip1 ≈ v slip2 to within 10%). and G̃ average of G̃w1 = −β v x x x x 2 1 x x 4.3

Tangential fluid force In either the static equilibrium state (where G̃w x = 0) or the dynamic steady state (where G̃w = 6 0), local force balance necessarily requires the stress tangential to the fluid-solid interface x to be the same on the two sides. Therefore, the hydrodynamic tangential fluid force per unit area, G̃fx , must be proportional to the slip velocity vxslip : G̃fx (x) = βvxslip (x), (4.2) such that G̃fx (x) + G̃w x (x) = 0 according to Eq. (41) (The MD evidence for this force balance will be presented in Section 7.) Physically, G̃fx is the hydrodynamic force along x exerted on a 22 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 8 −1 gxw(z) / Gxw ( σ ) 6 4 2 0 0 0.2 0.4 0.6 z /σ 0.8 1 Figure 13: Subdividing the BL into thin sections, we plot the reduced tangential wall force density as a function of distance z away from the boundary. The solid lines are averaged gxw (z) in thin sections at different x, normalized by the

corresponding total wall force per unit area The dashed line is a smooth Gaussian fit. It is seen that gxw (z) is a function sharply peaked at z ≈ z0 /2 Note that R z0 dz[gxw (x, z)/Gw x (x)] = 1. 0 BL fluid element by the surrounding fluids, and may be expressed as G̃fx (x) = Z z0 dz[∂x σ̃xx (x, z) + ∂z σ̃zx (x, z)] Z z0 dzσ̃xx (x, z), = σ̃zx (x, z0 ) + ∂x 0 (4.3) 0 using the fact that σ̃zx (x, 0) = 0. (More strictly, σ̃zx (x, 0− ) = 0 because there is no fluid below (0) 0 z = 0, hence no momentum transport across z = 0.) Here σ̃xx(zx) = σxx(zx) − σxx(zx) , with σxx (0) being the normal component and σzx the tangential component of the fluid stress tensor in the dynamic (static) configuration. 4.4 Sharp boundary limit The form of G̃fx in Eq. (43) is derived from the fact that the tangential wall force is distributed in a BL of finite thickness (Figs. 12 and 13) Now we take the sharp boundary limit by letting the tangential wall force

strictly concentrate at z = 0: g̃xw (x, z) = G̃w x (x)δ(z) with the same G̃w (x) per unit area. As shown in Fig. 13, the tangential wall force density is a sharply peaked x function. By taking the sharp boundary limit the normalized peaked function is replaced by Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 23 0.3 (a) vxslip 0.2 0.1 0 0.1 (b) Gxw(0) 0.0 −0.1 −0.2 −0.3 −30 −20 −10 0 x/σ 10 20 30 Figure 14: Slip velocity and tangential wall force (in reduced units) measured in the BL at the lower fluid-solid interface. (a) The slip velocity vxslip = vx + V is plotted as a function of x The solid line p denotes the dynamic symmetric case with V = 0.25 ²/m and H = 136σ; the dashed line denotes the p dynamic asymmetric case with V = 0.2 ²/m and H = 136σ The slip at the contact line (x ≈ 0) is near-complete, i.e, vxslip ≈ V (b) The tangential wall force is plotted as a function of x The solid line w denotes Gw x in the dynamic

symmetric case; the dashed line denotes Gx in the dynamic asymmetric case. w0 The dotted line denotes Gw0 x in the static symmetric case; the dot-dashed line denotes Gx in the static w0 asymmetric case. Note that Gx vanishes out of the contact-line region δ(z). Rewriting G̃fx in Eq (43) as G̃fx (x) = Z z0 dz[∂x σ̃xx (x, z) + ∂z σ̃zx (x, z)] Z z0 + dz[∂x σ̃xx (x, z) + ∂z σ̃zx (x, z)], = σ̃zx (x, 0 ) + 0− 0+ we obtain G̃fx (x) = σ̃zx (x, 0+ ) = βvxslip (x), (4.4) because local force balance requires ∂x σ̃xx + ∂z σ̃zx = 0 above z = 0+ . Therefore, in the sharp boundary limit σ̃zx varies from σ̃zx (x, 0− ) = 0 to σ̃zx (x, 0+ ) = G̃fx (x) at z = 0 such that (∇ · σ̃) · x̂ = G̃fx (x)δ(z), in balance with the tangential wall force density g̃xw (x, z) = G̃w x (x)δ(z). Equation (44) may serve as a boundary condition in hydrodynamic calculation if a continuum (differential) form of σ̃zx (x, 0+ ) is given. This will be accomplished

in Section 52 24 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 30 20 ~ 1 / Gxw 10 0 −10 −20 −30 −20 −10 0 10 20 1 / vxslip slip Figure 15: 1/G̃w (in reduced units). Symbols are MD data measured x plotted as a function of 1/vx p in the BL at different x locations. The circles denote the symmetric case with V = 025 ²/m and p H = 13.6σ; the squares denote the asymmetric case with V = 02 ²/m and H = 136σ The solid and dashed lines are calculated from Eq. (41) for the symmetric and asymmetric cases, respectively, as described in the text. The lower-right data segment corresponds to the lower BL, whereas the upper-left −1 segment corresponds to the upper BL. The slopes of the two dotted lines are given by β1,2 , which are proportional to the slip length. 5 Continuum hydrodynamic model For decades numerous models have been proposed to resolve the boundary condition problem for the contact-line motion [4,7,19,21,24,38,43], but so far none has

proved successful by reproducing the slip velocity profiles observed in MD simulations [28,40]. In particular, based on the extreme (velocity) variations in the slip region, a breakdown of local hydrodynamic description in the immediate vicinity of the MCL has been suggested [40]. The main purpose of this paper is to present a continuum hydrodynamic model that is capable of reproducing MD results in the molecular-scale vicinity of the MCL [34]. For this purpose, we have derived a differential form for Eq. (44) (the continuum GNBC, Eq (34)) using the Cahn-Hilliard hydrodynamic formulation of two-phase flow [7, 24]. Our model consists of the convection-diffusion equation in the fluid-fluid interfacial region (Eq. (31)), the Navier-Stokes equation for momentum transport (Eq. (32)), the relaxational equation for the composition at the solid surface (Eq. (33)), and the GNBC (Eq (34)) 5.1 Cahn-Hilliard free energy functional The CH free energy was proposed to phenomenologically describe

an interface between two coexisting phases [5]. In terms of the composition order parameter φ = (ρ2 − ρ1 )/(ρ2 + ρ1 ), the Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 25 CH free energy functional reads F = Z dr · ¸ 1 K (∇φ)2 + f (φ) , 2 p with f (φ) = − 21 rφ2 + 41 uφ4 . Two thermally stable phases are given by φ± = ± r/u at which ∂f /∂φ = 0. An interface can be formed between the phases of φ+ and φ− in coexistence 5.11 Chemical potential The chemical potential µ is defined by µ= δF = −K∇2 φ − rφ + uφ3 , δφ from which the diffusion current J = −M ∇µ is obtained with M being the mobility coefficient. The convection-diffusion equation (Eq. (31)) comes from the continuity equation ∂φ Dφ = + v · ∇φ = −∇ · J. Dt ∂t 5.12 Interfacial tension A few important physical quantities can be derived from the CH free energy. We first derive the interfacial tension γ for the interface formed between φ+

and φ− . In equilibrium, the spatial variation of φ is determined by the condition that µ(r) is constant, i.e, −K∂z2 φ − rφ + uφ3 = constant. Here the interface is assumed to be in the xy plane with the interface normal along the z direction and the constant equals to zero because limz±∞ φ = φ± and limz±∞ µ = 0. The interfacial profile is solved to be z φ0 (z) = φ+ tanh √ , 2ξ p with ξ = K/r being a characteristic length along the interface normal. The first integral is 1 − K (∂z φ)2 + f (φ) = C, 2 where the integral constant C equals f (φ± ). It follows the interfacial free energy per unit area, i.e, the interfacial tension, is given by ¸ Z · Z 1 2 γ = dz K (∂z φ) + f (φ) − f (φ± ) = dzK (∂z φ)2 . 2 Using the interfacial profile φ0 (z), we obtain √ √ Z 2 2Kφ2± Kφ2± 2 2r2 ξ −4 γ= √ dz̄ cosh z̄ = = . 3ξ 3u 2ξ 26 5.13 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 Capillary force and Young

stress We now turn to the forces arising from the interface. Consider a virtual displacement u(r) and the corresponding variation in φ, δφ(r) = −u(r) · ∇φ. The change of the free energy due to this δφ is i  h  2 1 ¸ Z · Z   ∂ K (∇φ) 2 ∂f (φ) δφ + dr δ (∂j φ) δF = dr   ∂φ ∂(∂j φ) Z Z = dr [µδφ] + ds [K∂n φδφ] Z Z £ Y ¤ = − dr [g · u] + ds σni ui , where g = µ∇φ is the capillary force density in the Navier-Stokes equation (Eq. (32)), and Y = −K∂ φ∂ φ is the tangential Young stress (the i direction is along the fluid-solid interface, σni n i i ⊥ n). The body force g(r) = µ∇φ can be reduced to the familiar curvature force in the sharp interface limit [6]. The unit vector normal to the level sets of constant φ is given by m = ∇φ/|∇φ| and µ∇φ = [−K∇2 φ − rφ + uφ3 ]|∇φ|m 2 φ − rφ + uφ3 ]|∇φ|m = −K∇2t φ|∇φ|m + [−K∂m where ∇t and ∂m denote the differentiations

tangential and normal to the interface respectively. For gently curved interfaces, the order parameter φ along the interface normal can be approx2 φ − rφ + uφ3 ≈ 0. Hence, imated by the one-dimensional stationary solution φ0 , i.e, −K∂m 2 µ∇φ ≈ −K∇t φ|∇φ|m, from which we obtain the desired relation µ∇φ ≈ K|∇φ|2 κm ≈ γκδ(lm )m, R R where κ = −∇2t φ/|∇φ| is the curvature and γ ≈ dlm K (∇φ)2 ≈ dlm K (∂m φ)2 is the interfacial tension, with lm being the coordinate along the interface normal and the interface located at lm = 0. For gently curved interfaces, K∂n φ ≈ K∂m φ cos θsurf , where n is the outward (solid) surface normal, m the (fluid-fluid) interface normal, and θsurf the angle at which the interface intersects Y = K∂ φ∂ φ at the solid surface (n · m = cos θsurf ). For the Rtangential Young stress σzx n x Y along x across the interface equals z = 0 where n = −z and i = x, the integral int dxσzx R R R

R to int dxK∂n φ∂x φ = ( int dφK∂m φ) cos θsurf , where int dφK∂m φ = int dlm K (∂m φ)2 = γ. Hence, Z int Y dxσzx = γ cos θsurf , (5.1) where θsurf may be the dynamic contact angle at the solid surface θdsurf or the static contact R Y is the tangential force per unit length at the contact line (aligned angle θssurf . This int dxσzx along y), exerted by the fluid-fluid interface of tension γ, which intersects the solid wall at the contact angle θsurf . So it equals to γ cos θsurf Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 5.14 27 Young’s equation The Young’s equation for the static contact angle θssurf Rcan be derived as well. Consider the interfacial free energy at the fluid-solid interface, Fwf = dsγwf (φ). Minimizing the total free energy F + Fwf with respect to φ at the solid surface yields · ¸ ∂γwf (φ) = 0, (5.2) K∂n φ + ∂φ φeq from which an equation of local tangential force balance £ Y¤ £ Y ¤ 0

+ ∂x γwf (φeq ) = 0, σ̃zx φeq = σzx + ∂x γwf (φ) φeq = σzx (5.3) Y = σ Y +∂ γ is obtained at z = 0. Here σ̃zx x wf is the uncompensated Young stress (first introduced zx 0 denotes the static Young stress in Eq. (35)), φeq is the equilibrium composition field, and σzx Y σzx (φeq ). Integrating Eq (53) along x across the interface leads to the Young’s equation R R 0 and ∆γ γ cos θssurf + ∆γwf = 0 (Eq. (37)), where γ cos θssurf = int dxσzx wf ≡ int dx∂x γwf (φ) is the change of fluid-solid interfacial free energy per unit area across the fluid-fluid interface. A microscopic picture for the Young’s equation as an (integrated) equation of tangential force balance will be elaborated in Section 7.21 5.2 Two boundary conditions Equations (5.2) and (53) are boundary conditions for the equilibrium state In the dynamic Y + ∂ γ (φ) = L(φ)∂ φ steady state, however, neither K∂n φ + ∂γwf (φ)/∂φ = L(φ) nor σzx x wf x vanishes. In

fact, the nonzero L(φ) is responsible for the relaxation of φ at the solid surface while the nonzero L(φ)∂x φ is necessary to a slip boundary condition that is able to account for the near-complete slip at the MCL. The convection-diffusion equation (Eq. (31)) is fourth-order in space Consequently, besides the usual impermeability condition ∂n µ = 0, one more boundary condition is needed. The dynamics of φ at the solid surface is plausibly assumed to be relaxational, governed by the first-order extension of Eq. (52) More explicitly, when the system is driven away from the equilibrium, both ∂φ/∂t + v · ∇φ and L(φ) become nonzero, and they are related to each other by a linear relation ∂φ + v · ∇φ ∝ L(φ). ∂t This leads to Eq. (33) with Γ introduced as a phenomenological parameter The GNBC (Eq. (34)) is obtained by substituting 0 (x, 0) = σ v (x, 0) + σ Y (x, 0) − σ 0 (x, 0) σ̃zx (x, 0+ ) = σzx (x, 0) − σzx zx zx zx v (x, 0) + σ Y (x, 0) + ∂

γ = σzx x wf zx v (x, 0) + σ̃ Y . = σzx zx (5.4) into Eq. (44) Here the hydrodynamic tangential stress σ̃zx is decomposed into a viscous v and a non-viscous component σ̃ Y . The viscous component is simply given by component σzx zx 28 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 v = η∂ v ; the non-viscous component is the uncompensated Young stress σ̃ Y , given by σzx z x zx Y = σ Y + ∂ γ (φ) (Eq. (35)) According to Eq (53), this uncompensated Young stress σ̃zx x wf zx Y along vanishes in the equilibrium state. But in a dynamic configuration, from the integral of σ̃zx x across the fluid-fluid interface (Eqs. (36), (37), and (38)) Z Y dxσ̃zx = γ cos θdsurf + ∆γwf = γ(cos θdsurf − cos θssurf ), int there is always a non-viscous contribution to the total tangential stress σ̃zx as long as the fluidfluid interface deviates from its static configuration. In Section 6 we will show that the GNBC, with the uncompensated Young

stress included, can account for the slip velocity profiles in the vicinity of the MCL, especially the near-complete slip at the contact line. In Sections 72 and 73 we will present more MD evidence supporting the GNBC. A “derivation” of the GNBC, based on the tangential force balance (Section 72) and the tangential stress decomposition (Section 7.3), will be given in Section 8 5.3 Dimensionless equations Dimensionlessp equations suitable for p numerical computation are obtained as follows. We scale φ by |φ± | = r/u, length by ξ = K/r, velocity by the wall speed V , time by ξ/V , and pressure/stress by ηV /ξ. In dimensionless forms, the convection-diffusion equation is ∂φ + v · ∇φ = Ld ∇2 (−∇2 φ − φ + φ3 ), ∂t the Navier-Stokes equation is ¸ · ∂v + (v · ∇) v = −∇p + ∇2 v + B(−∇2 φ − φ + φ3 )∇φ, R ∂t the relaxational equation for φ at the solid surface is " # √ ∂φ 2 surf + vx ∂x φ = −Vs ∂n φ − cos θs sγ

(φ) , ∂t 3 and the GNBC is [Ls (φ)]−1 vxslip √ # 2 cos θssurf sγ (φ) ∂x φ − ∂n vx . = B ∂n φ − 3 " (5.5) (5.6) (5.7) (5.8) Here sγ (φ) = (π/2) cos(πφ/2) is from the fluid-solid interfacial free energy γwf (φ) = (∆γwf /2) sin(πφ/2), which denotes a smooth interpolation between ±∆γwf /2. Five dimensionless parameters appear in the above equations. They are (1) Ld = M r/V ξ, √which is the ratio of a diffusion length M r/V to ξ, (2) R = ρV ξ/η, (3) B = r2 ξ/uηV = 3γ/2 2ηV , which is inversely proportional to the capillary number Ca = ηV /γ, (4) Vs = KΓ/V , and (5) Ls (φ) = η/β(φ)ξ, which is the ratio of the slip length ls (φ) = η/β(φ) to ξ, where β(φ) = (1 − φ)β1 /2 + (1 + φ)β2 /2. A numerical algorithm based on a fixed uniform mesh has been presented in Ref. [34] Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 6 29 Comparison of MD and continuum results To demonstrate the validity of

our continuum model, we have obtained numerical solutions that can be directly compared to the MD results for flow field and fluid-fluid interface shape. We have carried out the MD-continuum comparison in such a way that virtually no adjustable parameter is involved in the continuum calculations. This is achieved as follows There are totally nine material parameters in our continuum model. They are ρm , η, β, ξ, γ, |φ± |, M , Γ, and θssurf . (Note (1) For the asymmetric case, two unequal slip coefficients β1 and β2 are involved in β; (2) The three parameters ξ, γ, and |φ± | are equivalent to the three parameters K, r, and u in the CH free energy density; (3) θssurf is for ∆γwf = −γ cos θssurf .) Among the nine parameters, seven are directly obtainable (measurable) in MD simulations. They are ρm , η, β1,2 , ξ, γ, |φ± |, and θssurf . (The fluid mass density ρm is set in MD simulations, the viscosity η and the slip coefficients β1,2 can be measured in

suitable single-fluid MD simulations, the interfacial thickness ξ can be obtained by measuring the interfacial profile φ = (ρ2 − ρ1 )/(ρ2 + ρ1 ) in MD simulations, the interfacial tension γ can be obtained by measuring an integral of the pressure/stress anisotropy in the interfacial region [26], |φ± | = 1 means the total immiscibility of the two fluids, and the static contact angle θssurf is directly measurable.) The two phenomenological parameters M and Γ have been introduced to describe the composition dynamics in the interfacial region. Their values are fixed by an optimized MD-continuum comparison That is, one MD flow field is best matched by varying the continuum flow field with respect to the values of M and Γ. Once all the parameter values are obtained (7 measured in MD simulations and 2 fixed by one MD-continuum comparison), our continuum hydrodynamic model can yield predictions that can be readily compared to the results from a series of MD simulations with

different external conditions (V , H, and flow geometry). The overall agreement is excellent in all cases, thus demonstrating the validity of the GNBC and the hydrodynamic model. We emphasize that the MD-continuum agreement has been achieved both in the molecular-scale vicinity of the contact line and far way from the contact line. This opens up the possibility of not only continuum simulations of nano- and microfluidics involving immiscible components, but also macroscopic immiscible flow calculations that are physically meaningful at the molecular level. (Molecular-scale details may be resolved through the iterative grid redistribution method without significantly compromising computation efficiency, see [35, 37]). 6.1 6.11 Immiscible Couette flow Two symmetric cases In Figs. 16 and 17 we show the MD and continuum velocity fields for two symmetric cases of immiscible Couette flow. In MD simulations, these two cases have the same local properties (fluid density, temperature,

fluid-fluid interaction, wall-fluid interaction, etc) but different external conditions (H and V ). Correspondingly, the continuum results are obtained using the same set of nine material parameters ρm , η, β (= β1 = β2 ), ξ, γ, |φ± |, M , Γ, and θssurf . 30 6.12 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 Two asymmetric cases In Figs. 18 and 19 we show the MD and continuum velocity fields for two asymmetric cases of immiscible Couette flow. In MD simulations, these two cases have the same local properties (fluid density, temperature, fluid-fluid interaction, wall-fluid interaction, etc) but different external conditions (H and V ). Correspondingly, the continuum results are obtained using the same set of ten material parameters ρm , η, β1 , β2 , ξ, γ, |φ± |, M , Γ, and θssurf . In particular, among these parameters, β2 and θssurf are measured in MD simulations while all the others directly come from the symmetric cases. Therefore, the

comparison here is without adjustable parameters 6.13 From near-complete slip to uniform shear flow From Figs. 16, 17, 18, and 19, we see that at the MCL, the slip is near-complete, ie, vx ≈ 0 and |vxslip | ≈ V , while far away from the contact line, the flow field is not perturbed by the fluid-fluid interface and the single-fluid uniform shear flow is recovered. The slip amount in the uniform shear flow is 2ls V /(H + 2ls ), vanishing in the limit of H À ls . Here we encounter an intriguing question: In a mesoscopic or macroscopic system, what is the slip profile which consistently interpolates between the near-complete slip at the MCL and the no-slip boundary condition that must hold at regions far away? Large-scale MD and continuum simulations have been carried out to answer this question [35]. 6.14 Steady-state fluid-fluid interface In Fig. 20 we show the MD and continuum fluid-fluid interface profiles for one symmetric and one asymmetric cases whose velocity fields are

shown in Figs. 16 and 18 6.2 Immiscible Poiseuille flow In order to further verify that the continuum model is local and the parameter values are local properties, hence applicable to different flow geometries, we have carried out MD simulations and continuum calculations for immiscible Poiseuille flows. We find that the continuum model with the same set of parameters is capable of reproducing the MD results for velocity field and fluid-fluid interface profile, shown in Fig. 21 Similar to what we have observed in Couette flows, here at the MCL the slip is near-complete, i.e, vx ≈ 0 and |vxslip | ≈ V , while far away from the contact line, the flow field is not perturbed by the fluid-fluid interface and the single-fluid unidirectional Poiseuille flow is recovered. In particular, the slip amount in the unidirectional Poiseuille flow vanishes in the limit of H À ls . We emphasize that the overall agreement is excellent in all cases (from Fig. 16 to 21), therefore the validity of

the GNBC and the hydrodynamic model is well affirmed. 6.3 Flow in narrow channels It is generally believed that continuum hydrodynamic predictions tend to deviate more from the “true” MD results as the channel is further narrowed [42]. This tendency has indeed been Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 31 0.05 vx 0.00 −0.05 −0.10 −0.15 vz −0.20 0.05 0.00 −0.05 −0.10 −20 −10 0 x/ σ 10 20 Figure 16: Comparison of the MD (symbols) and continuum (lines) velocity profiles (vx (x) and vz (x) at different z levels) for a symmetric case of immiscible Couette flow (V = 0.25(²/m)1/2 and H = 136σ) The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425σ (circles and solid lines), 2125σ (squares and dashed lines), 3825σ (diamonds and dotted line), and 5.525σ (triangles and dot-dashed lines) The MD velocity profiles were measured by dividing the fluid space into 16 layers along

z, each of thickness H/16 = 0.85σ 0.00 vx −0.01 −0.02 −0.03 −0.04 0.02 vz 0.01 0.00 −0.01 −0.02 −30 −20 −10 0 x/ σ 10 20 30 Figure 17: Comparison of the MD (symbols) and continuum (lines) velocity profiles (vx (x) and vz (x) at different z levels) for a symmetric case of immiscible Couette flow (V = 0.05(²/m)1/2 and H = 272σ) The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.85σ (circles and solid lines), 425σ (squares and dashed lines), 765σ (diamonds and dotted line), and 11.05σ (triangles and dot-dashed lines) The MD velocity profiles were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 1.7σ 32 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 0.2 vx 0.1 0.0 −0.1 −0.2 vz 0.05 0.00 −0.05 −20 −10 0 10 x/ σ 20 Figure 18: Comparison of the MD (symbols) and continuum (lines) velocity profiles (vx (x) and vz (x) at

different z levels) for an asymmetric case of immiscible Couette flow (V = 0.2(²/m)1/2 and H = 136σ), shown at z = 0.425σ (circles and solid lines), 2975σ (squares and long-dashed lines), 5525σ (diamonds and dotted line), 8.075σ (up-triangles and dot-dashed lines), 10625σ (down-triangles and dashed lines), 13.175σ (left-triangles and solid lines) Although the solid lines are used to denote two different z levels, for each solid line, whether it should be compared to circles or left-triangles is self-evident (same for the next figure). The MD velocity profiles were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85σ 0.10 vx 0.05 0.00 −0.05 −0.10 0.04 vz 0.02 0.00 −0.02 −0.04 −40 −30 −20 −10 0 x/ σ 10 20 30 40 Figure 19: Comparison of the MD (symbols) and continuum (lines) velocity profiles (vx (x) and vz (x) at different z levels) for an asymmetric case of immiscible Couette flow (V = 0.1(²/m)1/2 and H =

272σ), shown at z = 0.85σ (circles and solid lines), 595σ (squares and long-dashed lines), 1105σ (diamonds and dotted line), 16.15σ (up-triangles and dot-dashed lines), 2125σ (down-triangles and dashed lines), 26.35σ (left-triangles and solid lines) The MD velocity profiles were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 1.7σ Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 33 14 12 10 z/σ 8 6 4 2 0 −4 −2 0 2 x/ σ 4 Figure 20: Comparison of the MD (symbols) and continuum (lines) fluid-fluid interface profiles, defined by ρ1 = ρ2 (φ = 0). The circles and dotted line denote the symmetric immiscible Couette flow with V = 0.25(²/m)1/2 and H = 136σ; the squares and dashed line denote the asymmetric immiscible Couette flow with V = 0.2(²/m)1/2 and H = 136σ The MD profiles were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85σ z/ σ 12 8 4 0 −4 −3 0.2

−2 x/σ −1 0 1 vx 0.1 0.0 −0.1 −0.2 −0.3 −20 −10 0 x/σ 10 20 Figure 21: Comparison of the MD (symbols) and continuum (lines) results for an asymmetric case of immiscible Poiseuille flow. An external force mgext = 005²/σ is applied on each fluid molecule in the x direction, and the two walls, separated by H = 13.6σ, move at a constant speed V = 051(²/m)1/2 in the −x direction in order to maintain a time-independent steady-state interface. (a) Fluid-fluid interface profiles, defined by ρ1 = ρ2 (φ = 0). (b) vx (x) at different z levels The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425σ (circles and solid line), 2125σ (squares and dashed line), 3.825σ (diamonds and dotted line), and 5525σ (triangles and dot-dashed line). The MD profiles were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85σ 34 Qian, Wang and Sheng / Commun. Comput Phys, 1

(2006), pp 1-52 0.00 vx −0.05 −0.10 −0.15 vz 0.05 0.00 −0.05 −15 −10 −5 0 x/ σ 5 10 15 Figure 22: Comparison of the MD (symbols) and continuum (lines) velocity profiles (vx (x) and vz (x) at different z levels) for a symmetric case of immiscible Couette flow (V = 0.25(²/m)1/2 and H = 68σ) The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425σ (circles and solid lines), 1.275σ (squares and dashed lines), 2125σ (diamonds and dotted line), and 2975σ (triangles and dot-dashed lines). The MD velocity profiles were measured by dividing the fluid space into 8 layers along z, each of thickness H/8 = 0.85σ observed but the deviation is not serious for H as small as 6.8σ, as shown in Fig 22 This deviation is presumably due to the short-range molecular layering induced by the rigid wall [23]. As the channel becomes narrower, the layered part of the fluids occupies a relatively larger space, thus making the

MD-continuum comparison less satisfactory. 6.4 Temperature effects Most of our MD results have been obtained by setting the temperature at 2.8²/kB , above the liquid-gas coexistence region. Such a high temperature was used to reduce the fluid layering at the solid wall [23]. Similar two-fluid simulations have also been performed for temperatures ranging from 1.2²/kB to 30²/kB We find that the MD results can always be reproduced by our continuum model, with material parameters (e.g viscosity, interfacial tension, and slip length) varying with the temperature. In Fig 23 we show the MD velocity profiles obtained at the temperatures 1.4²/kB and 28²/kB It can be seen that they are qualitatively very close to each other. The quantitative difference is due to the different material parameters at different temperatures. Finally we list in Table 1 the parameter values in the continuum hydrodynamic model, used for the MD-continuum comparison at T = 2.8²/kB Qian, Wang and Sheng /

Commun. Comput Phys, 1 (2006), pp 1-52 35 Table 1: Parameter values used in the continuum hydrodynamic calculations for the MD-continuum comparison at T = 2.8²/kB ρm ≈ 0.81m/σ 3 ls1 = η/β1 ≈ 1.3σ ξ ≈ 0.33σ √ M ≈ 0.023σ 4 / m² 7 √ η ≈ 1.95 ²m/σ 2 ls2 = η/β2 ≈ 1.3σ or 33σ γ ≈ 5.5²/σ 2 √ Γ ≈ 0.66σ/ m² |φ± | = 1 cos θssurf = 0 or ≈ 0.38 Molecular dynamics II We have formulated a continuum hydrodynamic model based on the CH free energy and the GNBC. The solutions of the model equations agree with the MD results remarkably well This indicates that our model captures the right physics, and hence more MD evidences can be obtained to support the continuum GNBC (Eq. (34)) which takes into account the uncompensated Young stress. This necessarily requires a reliable measurement of fluid stress near the solid surface, plus a decomposition of the tangential stress into two components, one being viscous and the other interfacial, as expressed

by Eq. (54) 7.1 Measurement of fluid stress Irving and Kirkwood [22] have shown that in the hydrodynamic equation of momentum transport, the stress tensor (flux of momentum) may be expressed in terms of the molecular variables as σ(r, t) = σ K (r, t) + σ U (r, t), where σ K is the kinetic contribution to the stress tensor, given by * + · ¸· ¸ X pi pi σ K (r, t) = − mi − V(r, t) − V(r, t) δ(xi − r) , mi mi i and σ U is the contribution of intermolecular forces to the stress tensor, given by + * 1 XX σ U (r, t) = − (xi − xj )Fij δ(xi − r) . 2 i j6=i Here mi , pi , and xi are respectively the mass, momentum, and position of molecule i, V(r, t) is the local average velocity, Fij is the force on molecule i due to molecule j, and h· · ·i means taking the average according to a normalized phase-space probability distribution function. The Irving-Kirkwood expression has been widely used for stress measurement in MD simulations. However, as pointed out by the

authors themselves [22], the above expression for σ U represents only the leading term in an asymptotic expansion, accurate when the interaction range 36 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 vx 0 (a) T=1.4 −0.1 −0.2 vx 0 (b) T=2.8 −0.1 −0.2 −30 −20 −10 0 x/ σ 10 20 30 Figure 23: Comparison of the MD (symbols) and continuum (lines) velocity profiles (vx (x) at different z levels) for two symmetric cases of immiscible Couette flow at different temperatures. (a) The case of T = 1.4²/kB , V = 025(²/m)1/2 , and H = 136σ, with a weak wall-fluid interaction and a high ratio of ρw to ρ (compared to case b here). (b) The case of T = 28²/kB , V = 025(²/m)1/2 , and H = 136σ, with ²wf = 1.16², σwf = 104σ, δwf = 1, and ρw /ρ = 23 The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425σ (circles), 2125σ (squares), 3825σ (diamonds), and 5.525σ (triangles) The MD

velocity profiles were measured by dividing the fluid space into 16 layers along z, each of thickness H/16 = 0.85σ The slip length for (a) is larger than that for (b) Therefore, far away from the contact line, the slip amount in (a) (≈ 0.1(²/m)1/2 ) is larger than that in (b) (≈ 0.05(²/m)1/2 ) is small compared to the range of hydrodynamic variation. As a consequence, this leading-order expression for σ U is not accurate enough near a fluid-fluid or a fluid-solid interface. Unfortunately, this point has not been taken seriously For the MCL problem, a knowledge of the stress distributions at both the fluid-fluid and the fluid-solid interfaces is of fundamental importance to a correct understanding of the underlying physical mechanism. Therefore, a reliable stress measurement method is imperative. To have spatial resolution along the x and z directions, the sampling region was evenly divided into bins, each ∆x = 0.425σ by ∆z = 085σ in size The stress components σxx and

σzx were obtained from the time averages of the kinetic momentum transfer plus the fluidfluid interaction forces across the fixed-x and z bin surfaces. More precisely, we have directly measured the x component of the fluid-fluid interaction forces acting across the x(z) bin surfaces, in order to obtain the xx(zx) component of σ U . For example, in measuring σU zx at a designated z-oriented bin surface, we recorded all the pairs of fluid molecules interacting across that surface. Here “acting/interacting across” means that the line connecting a pair of molecules intersects the bin surface (the so-called Irving-Kirkwood convention [22]). For those pairs, we then computed Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 37 1 A B 2 C 3 D 4 z 5 z bin surface normal along z x Figure 24: Schematic illustration for the measurement of the zx component of σ U . The horizontal solid lines (separated by short vertical lines) represent bin surfaces with surface

normal along the z direction. Circles denote fluid molecules. The dashed lines connect pairs of interacting molecules Here the bin surfaces and the molecules are projected onto the xz plane. Molecules that appear to be close to each other may not be in the interaction range if their distance along y is too large. A pair of interacting molecules may act across more than one bin surface. Here the (1,3) pair acts across the surfaces A and C while the (1,5) pair acts across the surfaces B and D. At each bin surface the stress measurement must run over all the pairs acting across that surface. For surface D, there are three pairs of interacting molecules (1,5), (2,4), and (2,5) that contribute to the zx component of σ U . σU zx at the given bin surface from σU zx = 1 X Fijx , δsz (i,j) where δsz is the area of z-oriented bin surface, (i, j) indicate all available pairs of fluid molecules interacting across the bin surface, with molecule i being “inside of ẑδsz ” and molecule

j being “outside of ẑδsz ” (molecule i is below molecule j), and Fijx is the x component of the force on molecule i due to molecule j. For a schematic illustration see Fig 24 7.2 Boundary-layer tangential force balance From the data of stress measurement, we now present the MD evidence for the BL tangential force balance, first introduced in Section 4.3 for obtaining Eq (42) from Eq (41) 7.21 Static tangential force balance We start from the tangential force balance in the static configuration (V = 0). As first pointed out in Section 4.2, the static tangential wall force Gw0 x shows molecular-scale features in the 38 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 contact-line region, due to the microscopic organization of fluid molecules there. Then according to the local force balance, static fluid stress must vary in such a way that the total force density f0 vanishes. An integrated form of the static tangential force balance is given by Gw0 x (x)+Gx

(x) = R z0 f0 w0 0, where Gw0 x (x) = 0 dzgx (x, z) and the static tangential fluid force Gx (x) is of the form Z z0 Z z0 0 0 0 0 dzσxx (x, z). (7.1) dz[∂x σxx (x, z) + ∂z σzx (x, z)] = σzx (x, z0 ) + ∂x Gfx0 (x) = 0 0 σxx 0 0 σzx Here and are the xx and zx components of fluid static configuration, both R z0stress0 in the w0 0 measured as reference quantities. In Fig 25 we show R 0 dzσxx , σzx (z0R), Gfx0 , and R z0 Gx 0 (which 0 is the same as in Fig. 14b) In the symmetric case, int dxσzx (x, z0 ) int dx∂x 0 dzσxx (x, z), R Rz and int dx 0 0 dzgxw0 (x, z) all vanish because θssurf = 90◦ . For the asymmetric case, Gw0 x (x) + f0 Gx (x) = 0 means Z z0 Z z0 0 0 dzσxx (x, z) + dzgxw0 (x, z) = 0, (7.2) σzx (x, z0 ) + ∂x 0 0 £ Y ¤ which corresponds to the continuum equation (Eq. (53)) σzx + ∂x γwf φeq = 0 at the solid 0 (x, z ) corresponds to the continuum σ Y (φ ) at the solid surface while surface. Here Rσzx 0 zx eq R z0 z 0 w0 0 + ∂x 0 dzσxx 0

dzgx corresponds to the continuum ∂x γwf (φeq ) at the solid surface. The Young’s equation (Eq. (37)) γ cos θssurf + ∆γwf = 0 is then obtained through integration, using Z Z Y 0 dxσzx = γ cos θssurf (7.3) dxσzx (x, z0 ) = int int and Z dx∂x int Z z0 0 0 dzσxx + Z dx int Z 0 z0 dzgxw0 = Z dx∂x γwf = ∆γwf . (7.4) int Here γ cos θssurf and ∆γwf are the two tangential forces per unit length along the contact line (along y), the former due to the tilt of the fluid-fluid interface (θssurf 6= 90◦ ) while the latter due to the different wall-fluid interactions for the two fluid species. In fact, equations (73) and (74) are the microscopic definitions for the two continuum quantities γ cos θssurf and ∆γwf in the Young’s equation, whose validity is based on the microscopic tangential force balance expressed in Eq. (72) 7.22 Dynamic tangential force balance f An integrated form of the dynamic tangential force balance is given by Gw x

(x) + Gx (x) = 0, R z0 f w where Gw x (x) = 0 dzgx (x, z) and the dynamic tangential fluid force Gx (x) is of the form Z z0 Z z0 Gfx (x) = dz[∂x σxx (x, z) + ∂z σzx (x, z)] = σzx (x, z0 ) + ∂x dzσxx (x, z). (7.5) 0 0 Here σxx and σzxRare the xx and zx components of fluid stress in the dynamic configuration. In z Fig. 26 we show 0 0 dzσxx , σzx (z0 ), and Gw x . From a comparison between Figs 25 and 26, we Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 0.5 0.2 (a) 0.4 (b) 0.1 0.3 0.2 0.0 0.1 −0.1 0.0 −0.1 −6 −4 −2 0 2 4 6 1.5 −0.2 −6 −4 −2 0 2 4 6 1.5 (c) 1.0 (d) 1.0 0.5 0.0 0.5 −0.5 −1.0 −1.5 39 0.0 −6 −4 −2 0 2 4 6 −6 −4 −2 0 2 4 6 Rz 0 0 , σzx (z0 ), Gfx0 , and Gw0 Figure 25: Profiles of 0 0 dzσxx x for the lower BL. The horizontal axes are x/σ R z0 0 2 We show 0 dzσxx (²/σ ) in (a) for the symmetric case and in (c) for the asymmetric case. For clarity, 0 0 = 0 far from

the interface, and in has been vertically displaced such that in the symmetric case, σxx σxx 0 0 3 the asymmetric case, σxx = 0 at the center of the interface. The profiles of σzx (z0 ), Gfx0 , and Gw0 x (²/σ ) 0 are plotted in (b) for the symmetric case and in (d) for the asymmetric case. The squares denote σzx (z0 ), w0 f0 the diamonds denote Gfx0 , and the solid triangles denote Gw0 x . Both (b) and (d) show Gx (x) = −Gx (x) R z0 w can see that the dynamic 0 dzσxx , σzx (z0 ), and Gx indeed show features seen from R z0 quantities w0 0 0 the static quantities 0 dzσxx , σzx (z0 ), and Gx , respectively. This is particularly evident in the asymmetric case where the static quantities vary more appreciably. The reason to take the static quantities as reference quantities is now clear: a hydrodynamic quantity must be obtained from the corresponding dynamic quantity by subtracting its static part, formally expressed as Q̃ = [Q]dynamic − [Q]static , where the over tilde

denotes the hydrodynamic quantity (first introduced for G̃w x in Section 4.2) In Fig. 27 we show the MD evidence for the BL hydrodynamic tangential force balance, which f is expressed as G̃w x (x) + G̃x (x) = 0. This equation is necessary for Eqs (41) and (42) to hold simultaneously. In summary, to verify the static/dynamic tangential force balance, we need to (1) identify the w(0) is distributed; (2) measure the normal and tangential BL where the tangential wall force Gx (0) (0) components of stress σxx and σzx according to the original definition of stress; (3) calculate the f (0) tangential fluid force Gx according to Eq. (71)/(75) 40 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 −9.5 0.8 (a) (b) 0.6 −10.0 0.4 0.2 −10.5 −11.0 0.0 −0.2 −20 −10 0 10 20 −7 −20 −10 0 10 20 1.8 (c) −8 (d) 1.2 −9 0.6 −10 0.0 −11 −20 −10 0 10 20 −20 −10 0 10 20 Rz BL. The horizontal axes are x/σ We Figure 26:

Profiles of 0 0 dzσxx , σzx (z0 ), and Gw x for the lower p R z0 2 show 0 dzσxx (²/σ ) in (a) for the symmetric case (V = 0.25 ²/m and H = 136σ) and in (c) for the p 3 asymmetric case (V = 0.2 ²/m and H = 136σ) The profiles of σzx (z0 ) and Gw x (²/σ ) are plotted in (b) for the symmetric case and in (d) for the asymmetric case. The squares denote σzx (z0 ) and the diamonds denote Gw contact line, Gw x . Note that in the vicinity of the x + σzx (z0 ) 6= 0. The importance of R z0 the x-gradient of the z-integrated normal stress ∂x 0 dzσxx is therefore evident. 0.00 ~w ~f Gx and −Gx −0.05 −0.10 −0.15 −0.20 −0.25 −20 −10 0 x/σ 10 20 Figure 27: Hydrodynamic force balance for the lower BL. The circles represent the symmetric case (V = p p 0.25 ²/m and H = 136σ); the squares represent the asymmetric case (V = 02 ²/m and H = 136σ) f w f The empty symbols denote G̃w x ; the solid symbols denote −G̃x . It is seen that G̃x (x) = −G̃x (x)

Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 7.3 41 Tangential Young stress Here we present the MD evidence for the decomposition of the tangential stress. In a dynamic v = η(∂ v + ∂ v ) configuration, away from the interfacial region the tangential viscous stress σzx z x x z is the only component in the (single-fluid) tangential stress σzx . But in the (two-fluid) interfacial v and a nonregion, the tangential stress σzx can be decomposed into a viscous component σzx Y : viscous component σzx v Y σzx = σzx + σzx , (7.6) v is still η(∂ v + ∂ v ) and σ Y is the tangential Young stress, satisfying where σzx z x x z zx Z Y Σd ≡ dxσzx (x, z) = γ cos θd (z). (7.7) int Here θd (z) is the dynamic interfacial angle at level z. Equation (76) is essential to obtaining v vanishes and σ Eq. (54) for σ̃zx (x, 0) In a static configuration, the viscous stress σzx zx becomes 0 σzx , satisfying Z Σs ≡ int 0 dxσzx (x, z) = γ cos θs (z),

(7.8) Y and σ 0 are where θs (z) is the static interfacial angle at level z. Fig 28 shows that both σzx zx nonzero in the interfacial region only. The inset to Fig 28 shows the evidence for Eqs (77) Y and σ 0 as the dynamic and static Young stresses. and (7.8), which identify σzx zx Equations (7.7) and (78) can be derived from the mechanical definition for the interfacial tension γ [26]: Z £ ¤ γ= dlm P⊥ (lm ) − Pk (lm ) , int i.e, γ is the integral (along the interface normal across the interface) of the difference between the normal and parallel components of the pressure, where lm is the coordinate along the interface normal m, and P⊥ and Pk are the pressure-tensor components normal and parallel to the interface, respectively. (Note that far away from the interface the pressure is isotropic and P⊥ = Pk ). When the interfacial angle θd (or θs ) is 90◦ , the interface normal m = x and the non-viscous stress tensor in the interfacial region is diagonal in the xyz

coordinate system:   −P⊥ 0 0 −Pk 0  = −P⊥ I + (P⊥ − Pk )(I − mm), σ non−viscous =  0 0 0 −Pk where I is the identity matrix. According to this expression, when the interfacial angle θd (or θs ) Y (or σ 0 ) arises from the interfacial stress deviates from 90◦ (see Fig. 29), the Young stress σzx zx anisotropy as the off-diagonal zx component of the microscopic stress tensor: Y (0) σzx = z · σ non−viscous · x = (P⊥ − Pk ) [z · (I − mm) · x] = −(P⊥ − Pk )mz mx = (P⊥ − Pk ) cos θd(s) sin θd(s) , 42 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 4 Σ d,s tangential Young stress 1.5 1.0 3 2 1 0 0 0.5 1 2 3 γ cos θd,s 4 0.0 −4 −2 0 x/σ 2 4 6 0 Figure 28: Dynamic and static Young stresses at z = z0 . The solid circles denote σzx in the static p Y symmetric case; the empty circles denote σzx in the dynamic symmetric case (V = 0.25 ²/m and 0 Y H = 13.6σ); the solid squares denote

σzx in the static asymmetric case; the empty squares denote σzx in p Y the dynamic asymmetric case (V = 0.2 ²/m and H = 136σ) Here σzx was obtained by subtracting the viscous component η(∂z vx + ∂x vz ) from the total tangential stress σzx . Inset: Σd,s plotted as a function of γ cos θd,s at different z levels. Here θd,s was measured from the time-averaged interfacial profiles The symbols have the same correspondence as in the main figure. The data indicate Σd,s = γ cos θd,s where mz = − cos θd(s) and mx = sin θd(s) . It follows that Σd(s) ≡ Z int Y (0) dxσzx (x, z) = = Z Zint int dx(P⊥ − Pk ) cos θd(s) sin θd(s) dlm (P⊥ − Pk ) cos θd(s) = γ cos θd(s) , where θd(s) is treated as a constant along x and dx sin θd(s) = dlm . 8 The generalized Navier boundary condition We want to “derive” the continuum GNBC using the MD results in Sections 4 and 7. For this purpose we first need to establish the correspondence between the stress

components measured in MD and those defined in the continuum hydrodynamics. This correspondence is essential to obtaining the microscopic dynamic contact angle θdsurf , which is defined in the continuum hydrodynamics (see Eq. (36) and (38)) but not directly measurable in MD simulations (because of the diffuse BL). Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 43 fluid−fluid interface dlm dx z m θd(s) x Y (0) Figure 29: Schematic illustration for the origin of the tangential Young stress σzx component of the microscopic stress tensor. 8.1 as an off-diagonal MD-continuum correspondence It has been verified that for a BL of finite thickness, the GNBC is given by Gfx0 (x) βvxslip (x) = Gfx (x) Z z− 0 £ ¤ £ ¤ ∂ 0 0 = dz σxx (x, z) − σxx (x, z) + σzx (x, z0 ) − σzx (x, z0 ) , ∂x 0 (8.1) in which only MD measurable quantities are involved. Now we interpret these MD-measured quantities in terms of the various continuum variables in the

hydrodynamic model. In so doing it is essential to note the following. (1) σxx can be decomposed into a molecular component HD . Meanwhile, σ 0 can be composed into and a hydrodynamic component: σxx = Txx + σxx xx 0 = T HS the same molecular component and a hydrostatic component: σxx xx + σxx . Physically, 0 Txx is the normal stress σxx for the case of a flat, static fluid-fluid interface. Such an interface exists in the symmetric case (θssurf = 90◦ ) for any value of H. It also exists in the asymmetric case (θssurf 6= 90◦ ) for H ∞ (with vanishing curvature ∼ 1/H) In either case, the HS vanishes according to the Laplace’s interface has zero curvature and the hydrostatic stress σxx 0 HS equation: σxx Txx as σxx 0. (To be more precise, γ cos θssurf and ∆γwf should be defined in Eqs (73) and (74) when the fluid-fluid interface has zero curvature However, in the asymmetric case where ∆γwf (= −γ cos θssurf ) is nonzero, there is a static interfacial

curvature HS , which should be subtracted from σ 0 in the ≈ 2 cos θssurf /H. This results in a hydrostatic σxx xx 0 replaced by T . Meanleft-hand side of Eq (74) That is, ∆γwf should be obtained with σxx xx while, due to the static interfacial curvature, the static interfacial angle θs (z0 ) determined by R 0 (x, z ) is a bit different from the true θ surf . In fact, cos θ (z ) deviates cos θs (z0 ) = γ −1 int dxσzx s s 0 0 from cos θssurf by the BL-integrated curvature ≈ 2z0 cos θssurf /H.) The molecular component Txx exists even if there is no hydrodynamic fluid motion or fluid-fluid interfacial curvature. On HD arises from the hydrodynamic fluid motion the contrary, the hydrodynamic component σxx HD becomes σ HS , which comes from the and interfacial curvature. In the static configuration, σxx xx 44 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 interfacial curvature. (2) σzx (x, z0 ) can be decomposed into a viscous component R plusY a

Young component : σzx (x, z0 ) = v Y v σzx (x, z0 ) + σzx (x, z0 ) with σzx = η(∂z vx + ∂x vz ) and int dxσzx (x, z0 ) = γ cos θd (z0 ) (Eqs. (76) and (7.7)) R 0 0 (x, z ) is the static Young stress: i.e, (3) σzx 0 int dxσzx (x, z0 ) = γ cos θs (z0 ) (Eq. (78)) Using the above relations, we integrate Eq. (81) along x across the fluid-fluid interface and obtain ¸ Z ·Z z 0 Z HD slip v dzσxx (x, z) + dxβvx (x) = ∆ dxσzx (x, z0 ) + γ cos θd (z0 ) 0 int int ·Z z 0 ¸ HS −∆ dzσxx (x, z) − γ cos θs (z0 ), (8.2) 0 where ∆ hR z0 0 HD(HS) dzσxx i ∆ HD(HS) is the change of the z-integrated σxx ·Z z0 0 HD(HS) dzσxx ¸ ≡ Z dx∂x int Z 0 z0 across the interface: HD(HS) dzσxx . According to the Laplace’s equation, the hydrostatic stress is directly related to the static curvature κs : HS −∆σxx = γκs , and the z-integrated curvature −∆ Z z0 0 R z0 0 HS dzσxx dzκs equals to cos θs (z0 ) − cos θssurf .

Hence, =γ Z z0 0 h i dzκs = γ cos θs (z0 ) − cos θssurf . (8.3) Substituting Eq. (83) into Eq (82) yields Z int dxβvxslip (x) = ∆ Z 0 z0 HD dzσxx (x, z) + Z int v dxσzx (x, z0 ) + γ cos θd (z0 ) − γ cos θssurf . (8.4) In order to interpret Eq. (84) in the continuum hydrodynamic formulation with a sharp BL, it is essential to note the following: (1) The sum of the first three terms on the right-hand side of Eq. (84) is the net fluid force along x exerted on the three fluid sides of a BL fluid element in the interfacial region; (2) The last Rterm in the right-hand side of Eq. (84), −γ cos θssurf , is the net wall force along x, ∆γwf = int dx∂x γwf , which arises from the wall-fluid interfacial free energy jump across the fluid-fluid interface, in accordance with the Young’s equation ∆γwf + γ cos θssurf = 0. Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 8.2 45 Extrapolated dynamic contact angle Now we take the

sharp boundary limit to relate the net fluid force Z Z z0 HD v dzσxx (x, z) + dxσzx (x, z0 ) + γ cos θd (z0 ) ∆ int 0 in Eq. (84) to the tangential stresses (viscous and non-viscous) at the solid surface The purpose of doing so is to obtain the surface contact angle θdsurf through extrapolation. Note that θdsurf is not directly measurable in MD simulations due to the diffuse BL. Only the extrapolated θdsurf can be compared to the contact angle in continuum calculations. In Section 4.4, we take the sharp boundary limit by assuming a tangential wall force conw w centrated at z = 0: g̃xw (x, z) = G̃w x (x)δ(z). While g̃x becomes a δ function, G̃x (x) per unit area remains the same. Using the equation of local force balance ∂x σ̃xx + ∂z σ̃zx = 0 above z = 0+ , we obtain σ̃zx (x, 0+ ) = G̃fx (x) as the tangential stress at the solid surface (Eq. (44)) The extrapolation here follows this spirit. We turn to the Stokes equation in the BL: v v −∂x p + ∂x σxx +

∂z σzx + µ∂x φ = 0, (8.5) obtained from the x-component of Eq. (32) by dropping the inertial term and the external force term. Integrating Eq (85) along z across the BL and then along x across the fluid-fluid interface, we obtain ½Z z0 ¾ Z v v ∆ dz [−p(x, z) + σxx (x, z)] + dxσzx (x, z0 ) + γ cos θd (z0 ) 0 int Z v dxσzx = (x, 0) + γ cos θdsurf . (8.6) int Two relations have been used in obtaining Eq. (86): (1) The capillary force density in the sharp interface limit [6] is given by µ∂x φ γκδ(x − xint ), where κ is the interfacial curvature and xint the location of the interface along x (see Section 5.13) (2) The z-integrated curvature gives Z z0 dzκ = cos θd (z0 ) − cos θdsurf . 0 The local force balance along x is expressed by REq. (85) R z0 Accordingly, the tangential force balance region ( int dx 0 dz) is expressed by Eq. (86), where ª R z0 for the BL fluidsv in the integration dz [−p(x, z) + σ (x, z)] is the net fluid force on the left

and right (∓x-oriented) surfaces, ∆ xx 0 R R v v int dxσzx (x, z0 ) + γ cos θd (z0 ) is the fluid force on the z = z0 surface, and int dxσzx (x, 0) + γ cos θdsurf is the tangential fluid force at the solid surface. 46 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 (a) z=z0 (b) nonzero tangential stress A zero z=0 tangential stress wall force distributed inside B D z=z0 C z=0 +− z=0 wall force concentrated at z=0 Figure 30: (a) For the real tangential wall force continuously distributed between z = 0 and z = z0 , the tangential stress in the fluid is continuous, vanishing at z = 0. (b) In the sharp boundary limit, the tangential wall force is considered to be concentrated at z = 0. Accordingly, the net fluid force on the fluids bound by A, B, C, and D (at z = 0+ ) is considered to be zero. In other words, the net fluid force on the three surfaces A, B, and C is fully transmitted to the tangential stress at the surface D at z = 0+ (Eq. (86)) So

there arises an abrupt change of the tangential stress from 0 at z = 0− to G̃fx at z = 0+ , by which the concentrated wall force is balanced. v with σ HD : Substituting Eq. (86) into Eq (84) and identifying the normal stress −p + σxx xx v HD −p + σxx = σxx , we obtain Z int dxβvxslip (x) = Z int v dxσzx (x, 0) + γ cos θdsurf − γ cos θssurf , (8.7) which is identical to the integral of the continuum GNBC along x across the fluid-fluid interface (Eqs. (34), (35), (36), and (38)) By doing so, we have assumed the tangential wall force is concentrated at z = 0. (This leads to Eqs (85) and (86) between z = 0+ and z = z0 ) In essence, the above extrapolation is to obtain the tangential stresses (viscous and non-viscous) at the solid surface when the limit of g̃xw (x, z) = G̃w x (x)δ(z) is taken, as illustrated in Fig. 30 For w G̃x (x) distributed in the diffuse BL, there is actually no tangential stress at the solid surface. Only in the sharp boundary limit

does a nonzero tangential stress appear at z = 0+ , equal to the net fluid force accumulated from z = 0− to z = z0 in the diffuse BL. R It is worth emphasizing that the right-hand side of Eq. (87) is int dxσ̃zx (x, 0+ ), where σ̃zx (x, 0+ ) is the extrapolated tangential stress in Eq. (44) Its continuum expression is given by Eq. (54) Equation (87) concludes our “derivation” of (an integrated form of) the continuum GNBC from the MD results presented in Sections 4 and 7. Equation (8.6) constitutes theR basis for obtaining Rthe dynamic contact angle θdsurf from MD z v ) = z0 dzσ HD is a sharp drop across the fluiddata. The dominant behavior of 0 0 dz (−p + σxx xx 0 fluid interface. This stress drop implies a large curvature in the BL, £ ¤ R z which pullsR zthe extrapolated 0 θdsurf closer to θssurf . We show the BL-integrated normal stress 0 0 dzσ̃xx = 0 0 dz σxx − σxx in Fig. 31, where the large stress drop across the fluid-fluid interface is clearly seen (In

the 0 and T , σ̃ asymmetric case, due to the small difference between σxx xx xx Ris not precisely the z HD HD HS HS is on the hydrodynamic stress σxx . In fact, σ̃xx = σxx − σxx Nevertheless, ∆ 0 0 dzσxx surf order of 2γz0 cos θs /H, much smaller than the magnitude of the stress drop shown in Fig. 31) Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 47 BL−integrated normal stress −8.5 −9.0 −9.5 −10.0 −10.5 −11.0 −20 −10 0 x/ σ 10 20 £ ¤ Rz Rz 0 (x, z) plotted as a function of x. The circles denote Figure 31: 0 0 dzσ̃xx (x, z) = 0 0 dz σxx (x, z) − σxx p the symmetric case (V = 0.25 ²/m and H = 136σ) and the squares denote the asymmetric case p 0 0 (V = 0.2 ²/m and H = 136σ) For clarity, σxx was vertically displaced such that σxx = 0 far from the 0 interface in the symmetric case, and for the asymmetric case, σxx = 0 at the center of the interface (same as in Fig. 25) Rz In the partial-slip region at the

vicinity of the MCL, 0 0 dzσ̃xx shows a fast variation along x as well. This means that the BL tangential force balance cannot be established unless the gradient of the BL-integrated normal stress is taken into account (see Eqs. (43), (71), and (7.5)) It is worth pointing out that the stress variation depicted in Fig 31 is also produced by the continuum hydrodynamic calculations, in semi-quantitative agreement with the MD data. In Fig. 32 we show the continuum profiles of layer-integrated pressure, obtained for a symmetric case of Couette flow. It is readily seen that the magnitude of the pressure change across the fluid-fluid interface decays away from the solid surface quickly. This conforms to the interface profile shown in Fig. 20: the interfacial curvature quickly decreases with the increasing distance from the solid surface. 8.3 The importance of the uncompensated Young stress According to Eq. (76), the hydrodynamic tangential stress σ̃zx (x, z0 ) can be decomposed into v

and the non-viscous component σ̃ Y : the viscous component σzx zx v Y σ̃zx (x, z0 ) = σzx (z) + σ̃zx (x, z0 ). v (x, z ) = In Fig. 33 we show that away from the interfacial region the tangential viscous stress σzx 0 Y = σ η(∂z vx + ∂x vz )(x, z0 ) is the only nonzero component, but in the interfacial region σ̃zx zx − v − σ 0 = σ Y − σ 0 is dominant, thereby accounting for the failure of the Navier boundary σzx zx zx zx condition to describe the contact line motion. Therefore away from the MCL region the Navier Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 layer−integrated pressure 48 0.6 0.4 0.2 0.0 −0.2 −20 −10 0 x/ σ 10 20 p R z+z /2 Figure 32: z−z00/2 dzp(x, z) obtained for the symmetric case of V = 0.25 ²/m and H = 136σ, plotted as a function of x. The profiles are symmetric about the center plane z = H/2, hence only the lower half is shown at z = 0.425σ (solid line), 2125σ (dashed line), 3825σ (dotted line),

and 5525σ (dotdashed lines) The solid line denotes the BL-integrated pressure, in semi-quantitative agreement with Rz the 0 0 dzσ̃xx (x, z) profile (circles) shown in Fig. 31 (Note that in Fig 31 the MCL is placed at x ≈ 0 while here it is shifted to x ≈ −3.9σ, same as in Figs 16, 20, and 23b) boundary condition is valid, but in the interfacial region it clearly fails to describe the contact line motion. 0 in the BL. The dominant behavior is We have measured the z-integrated σ̃xx = σxx − σxx a sharp drop across the interface, as shown in Fig. 31 for both the symmetric and asymmetric cases. As already pointed out in Section 82, this stress drop means a large curvature in the BL, which pulls the extrapolated θdsurf closer to θssurf . The value of θdsurf obtained through extrapolation is 88 ± 0.5◦ for the symmetric case and 63 ± 05◦ for the asymmetric case at the lower boundary, and 64.5 ± 05◦ at the upper boundary These values are noted to be very close to

θssurf . Yet the small difference between the dynamic and static contact angles is essential in accounting for the near-complete slip at the MCL. v In essence, our results show that in the vicinity of the MCL, the tangential viscous stress σzx as postulated by the Navier boundary condition can not give rise to the near-complete MCL slip Y in combination with the gradient of without taking into account the tangential Young stress σzx the (BL-integrated) normal stress σxx . For the static configuration, the Young stress is balanced by the normal stress gradient, leading to the Young’s equation. It is only for a MCL that there is a component of the Young stress which is no longer balanced by the normal stress gradient, and this uncompensated Young stress is precisely the additional component captured by the GNBC but missed by the traditional Navier boundary condition. Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 49 −3 tangential stress (εσ ) 0.6 0.4

0.2 0.0 −10 −5 0 5 10 15 x/σ Figure 33: Two components of the hydrodynamic tangential stress at z = z0 , plotted as a function Y of x. The symbols connected by dashed lines denote σ̃zx ; the solid and dotted lines represent the viscous component. Here the symmetric case is represented by circles and solid line; the asymmetric case represented by squares and dotted line. In the contact line region the non-viscous component is almost one order of magnitude larger than the viscous component. The difference between the two components, however, diminishes towards the boundary, z = 0, due to the large interfacial pressure drop (implying a large curvature) in the BL, thereby pulling θdsurf closer to θssurf . 9 Concluding remarks In summary, we have found the slip boundary condition, i.e, the GNBC, for the contact-line motion. Based on this finding, we have formulated a hydrodynamic model that is capable of reproducing the MD results of slip profile, including the

near-complete slip at the MCL. It should be noted, however, that the present continuum formulation can not calculate fluctuation effects that are important in MD simulations. Based on the results and experiences obtained in the present study of the MCL, we have been working on the following. (1) Large scale MD simulations on two-phase immiscible flows show that associated with the MCL, there is a very large 1/x partial-slip region, where x denotes the distance from the contact line. This power-law partial-slip region has been reproduced in large-scale adaptive continuum calculations [35, 37] based on the local, continuum hydrodynamic formulation presented in this paper. MD simulations and continuum solutions both indicate the existence of a universal slip profile in the Stokes-flow regime, well described by v slip (x)/V = 1/(1 + x/als ), where v slip is the slip velocity, V the speed of moving wall, ls the slip length, and a is a numerical constant ∼ 1 [35]. A large 1/x partial slip

region is significant, because the outer cutoff length scale directly determines the integrated effects, such as the total steady-state dissipation. While in the past the 1/x stress variation away from the MCL has been known [20], to our knowledge the fact that the partial slip also exhibits the same spatial dependence has not been previously 50 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 seen, even though the validity of the Navier boundary condition at high shear stress has been verified [2, 8, 39, 41]. (2) By explicitly taking into account the long-range wall-fluid interactions, our hydrodynamic model for two-phase immiscible flow at the solid surface can be used to investigate the dry spreading of a pure, nonvolatile liquid, attracted towards the solid by long-range van der Waals forces [12]. The precursor film, driven by the gradient of disjoining pressure due to the long-range force [9], was observed decades ago [13, 17, 25]. Nevertheless, the theoretical

analysis remains to be difficult, because of a wide separation of different length scales involved [12, 14, 18, 32]. Application of our model to this multi-scale problem is currently underway. (3) Decades ago it became well known that the driven cavity flow is incompatible with the noslip boundary condition, since the latter would lead to stress singularity and infinite dissipation (known as corner-flow singularity) [3, 29]. While MD studies [27] have clearly demonstrated relative fluid-wall slipping near the corner intersection, the exact rule that governs this relative slip has been left unresolved. Based on the simulation technique developed in the present study, we have verified the validity of the Navier boundary condition in governing the fluid slipping in driven cavity flows. We have used this discovery to formulate a continuum hydrodynamics whose predictions are in remarkable quantitative agreement with the MD simulation results at the molecular level [33]. Acknowledgments The

authors would like to thank Weinan E and Weiqing Ren for helpful discussions. This work was partially supported by the grants DAG03/04.SC21 and RGC-CERG 604803 References [1] M. Allen and D Tildesley, Computer Simulation of Liquids, Clarendon, New York, 1987 [2] J-L. Barrat and L Bocquet, Large slip effect at a nonwetting fluid-solid interface, Phys Rev Lett, 82 (1999), p. 4671; Influence of wetting properties on hydrodynamic boundary conditions at a fluid/solid interface, Faraday Discuss., 112 (1999), p 119 [3] G.K Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, 1991 [4] T. D Blake and J M Haynes, Kinetics of liquid/liquid displacement, J Colloid Interf Sci, 30 (1969), p. 421; T D Blake, Dynamic contact angles and wetting kinetics, in: J C Berg (Ed), Wettability, Marcel Dekker Inc., 1993, p 251 [5] J. W Cahn and J E Hilliard, Free energy of a nonuniform system I Interfacial free energy, J Chem. Phys, 28 (1958), p 258 [6] R. Chella and J Vinals,

Mixing of a two-phase fluid by cavity flow, Phys Rev E, 53 (1996), p 3832 [7] H. Y Chen, D Jasnow and J Vinals, Interface and contact line motion in a two phase fluid under shear flow, Phys. Rev Lett, 85 (2000), p 1686 [8] M. Cieplak, J Koplik and J R Banavar, Boundary conditions at a fluid-solid interface, Phys Rev Lett., 86 (2001), p 803 [9] B. V Derjaguin, Theory of the capillary condensation and other capillary phenomena taking into account the disjoining effect of long-chain molecular liquid films (in Russian), Zh. Fiz Khim, 14 (1940), p. 137 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 51 [10] E. B Dussan V, The moving contact line: the slip boundary condition, J Fluid Mech, 77 (1976), p. 665 [11] E. B Dussan, V, On the spreading of liquids on solid surfaces: static and dynamic contact lines, Ann. Rev Fluid Mech, 11 (1979), p 371 [12] P. G de Gennes, Wetting: statics and dynamics, Rev Mod Phys, 57 (1985), p 827 [13] H. Ghiradella, W Radigan and H L Frisch,

Electrical-resistivity changes in spreading liquid-films, J. Colloid Interf Sci, 51 (1975), p 522 [14] K. B Glasner, Spreading of droplets under the influence of intermolecular forces, Phys Fluid, 15 (2003), p. 1837 [15] N. G Hadjiconstantinou, Combining atomistic and continuum simulations of contact-line motion, Phys. Rev E, 59 (1999), p 2475 [16] N. G Hadjiconstantinou, Hybrid atomistic-continuum formulations and the moving contact-line problem, J. Comput Phys, 154 (1999), p 245 [17] W. B Hardy, The spreading of fluids on glass, Philos Mag, 38 (1919), p 49 [18] H. Hervet and P G de Gennes, Physique des surfaces et des interfaces Dynamique du mouillage: films precurseurs sur solide, C. R Acad Sci Paris II, 299 (1984), p 499 [19] L. M Hocking, A moving fluid interface Part 2 The removal of the force singularity by a slip flow, J. Fluid Mech, 79 (1977), p 209 [20] C. Hua and L E Scriven, Hydrodynamic model of steady movement of a solid/liquid/fluid contact line, J. Colloid Interf Sci,

35 (1971), p 85 [21] C. Huh and S G Mason, The steady movement of a liquid meniscus in a capillary tube, J Fluid Mech., 81 (1977), p 401 [22] J. H Irving and J G Kirkwood, The statistical mechanical theory of transport processes IV The equations of hydrodynamics, J. Chem Phys, 18 (1950), p 817 For the validity of the IrvingKirkwood expression, see the paragraph following equation (515) in this paper [23] J. N Israelachivili, Intermolecular and Surface Forces (2nd ed), Academic Press, London, 1992 [24] D. Jacqmin, Calculation of two-phase Navier-Stokes flows using phase-field modeling, J Comput Phys. 155 (1999), p 96; D Jacqmin, Contact-line dynamics of a diffuse fluid interface, J Fluid Mech., 402 (2000), p 57 [25] H. P Kavehpour, B Ovryn and G H McKinley, Microscopic and macroscopic structure of the precursor layer in spreading viscous drops, Phys. Rev Lett, 91 (2003), p 196104 [26] J. G Kirkwood and F P Buff, The statistical mechanical theory of surface tension, J Chem Phys, 17

(1949), p. 338 [27] J. Koplik and J R Banavar, Corner flow in the sliding plate problem, Phys Fluids, 7 (1995), p. 3118 [28] J. Koplik, J R Banavar and J F Willemsen, Molecular dynamics of Poiseuille flow and moving contact lines, Phys. Rev Lett, 60 (1988), p 1282; J Koplik, J R Banavar and J F Willemsen, Molecular dynamics of fluid flow at solid surfaces, Phys. Fluids A, 1 (1989), p 781 [29] H. K Moffatt, Viscous and resistive eddies near a sharp corner, J Fluid Mech, 18 (1964), p 1 [30] C. L M H Navier, Memoire sur les lois du mouvement des fluides, Memoires de l’Academie Royale des Sciences de l’Institut de France, 6 (1823), 389-440. [31] L. M Pismen and Y Pomeau, Disjoining potential and spreading of thin liquid layers in the diffuseinterface model coupled to hydrodynamics, Phys Rev E, 62 (2000), p 2480; L M Pismen, Mesoscopic hydrodynamics of contact line motion, Colloid Surface A, 206 (2002), p 11 [32] L. M Pismen, B Y Rubinstein and I Bazhlekov, Spreading of a wetting film

under the action of van der Waals forces, Phys. Fluid, 12 (2000), p 480 [33] T. Z Qian and X P Wang, Driven cavity flow: from molecular dynamics to continuum hydrody- 52 Qian, Wang and Sheng / Commun. Comput Phys, 1 (2006), pp 1-52 namics, SIAM Multiscale Model. Simul, 3 (2005), p 749 [34] T. Z Qian, X P Wang and P Sheng, Molecular scale contact line hydrodynamics of immiscible flows, Phys. Rev E, 68 (2003), p 016306 [35] T. Z Qian, X P Wang and P Sheng, Power-law slip profile of the moving contact line in two-phase immiscible flows, Phys. Rev Lett, 93 (2004), p 094501 [36] W. Ren and Weinan E, Heterogeneous multiscale method for the modeling of complex fluids and micro-fluidics, J. Comput Phys, 204 (2005), p 1 [37] W. Ren and X P Wang, An iterative grid redistribution method for singular problems in multiple dimensions, J. Comput Phys, 159 (2000), p 246 [38] P. Seppecher, Moving contact lines in the Cahn-Hilliard theory, Int J Engng Sci, 34 (1996), p 977 [39] P. A Thompson and M

O Robbins, Shear flow near solids: epitaxial order and flow boundary conditions, Phys. Rev A, 41 (1990), p 6830 [40] P. A Thompson and M O Robbins, Simulations of contact-line motion: slip and the dynamic contact angle, Phys. Rev Lett, 63 (1989), p 766; P A Thompson, W B Brinckerhoff and M O Robbins, Microscopic studies of static and dynamic contact angles, J. Adhesion Sci Tech, 7 (1993), p. 535 [41] P. A Thompson and S M Troian, A general boundary condition for liquid flow at solid surfaces, Nature, 389 (1997), p. 360 [42] K. P Travis and K E Gubbins, Poiseuille flow of Lennard-Jones fluids in narrow slit pores, J Chem. Phys, 112 (2000), p 1984 [43] M. Y Zhou and P Sheng, Dynamics of immiscible-fluid displacement in a capillary tube, Phys Rev. Lett, 64 (1990), p 882; P Sheng and M Y Zhou, Immiscible-fluid displacement: Contact-line dynamics and the velocity-dependent capillary pressure, Phys. Rev A, 45 (1992), p 5694