What is the force balance equation in the BL?

In the BL, the wall-fluid interaction means the existence of a dynamic, tangential wall force density g̃x w such that the force balance equation is given by (“•s̃)• x̂1 g̃x w50 inside the BL.

How did the authors increase the wall thickness in their simulations?

of increasing the wall thickness in their simulations from two layers of wall molecules to four layers and to infinite layers ~by using the continuum approximation beyond the four layers!.

What is the rate of work done to the system in the steady state?

In the steady state, the external work done to the system is fully dissipated in the system through convection-diffusion of the composition, slipping at surface, and shear viscosity in the bulk.

What is the udsur f value for a moving CL?

It is only for a moving CL 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 NBC.

How is the rate of work calculated?

The rate of that work is positive, given by the integrated local force times the wall velocity, i.e., *dxuG̃x wuV5*dxbuvx slipuV per unit length.

What is the rate of dissipation due to the moving CL?

For the symmetric case, the resulting heat generation rate due to the CL is thus bV2WsL ~for one wall!, where L is the length of the CL and Ws defines the width of the CL region:Ws5 1 VE ~ uvxslipu2v0slip!dx .

What was the dw f for specifying the wetting property of the fluid?

The wall-fluid interaction was also modeled by a LJ potential Uw f , with energy and range parameters ew f51.16e and sw f51.04s , and a dw f for specifying the wetting property of the fluid.

What is the way to verify that the boundary conditions and the parameter values are local properties?

To further verify that the boundary conditions and the parameter values are local properties and hence applicable to flows with different macroscopic conditions, the authors have varied the wall speed V, the system size H, and the flow geometry to check that the same set of parameters plus the GNBC are valid for reproducing ~a!

What is the first order approximation of the CH free energy functional?

As a first-order approximation, the authors formulate a hydrodynamic model based on the GNBC and the CH free energy functional @21# that has been successful in the calculations of fluid-fluid interfacial phenomena:F@f#5E drF12 K~¹f!21 f ~f!G , ~6! where f5(r22r1)/(r21r1), f (f)52 1 2 rf 21 14 uf 4, and K, r, u are the parameters that can be directly obtained from MD simulations through the interface profile thickness j 5AK/r @22#, the interfacial tension g52A2r2j/3u , and the two homogeneous equilibrium phases given by the condition of ] f /]f50, yielding f656Ar/u (561 in their case!.

(Open Access) Molecular scale contact line hydrodynamics of immiscible flows. (2003) | Tiezheng Qian

Q: How was the kinetic momentum transfer measured?

Gxw 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; sxx(nx) was obtained from the time averages of the kinetic momentum transfer plus the fluidfluid interaction forces across the constant-x(z) bin surfaces,6-2and vx(z) was measured as the time-averaged velocity component~s!

Molecular scale contact line hydrodynamics of immiscible ﬂows

Tiezheng Qian

Department of Physics and Institute of Nano Science and Technology, The Hong Kong University of Science and Technology,

Clear Water Bay, Kowloon, Hong Kong, China

Xiao-Ping Wang

Department of Mathematics, The Hong Kong University of Science and Technology, Clear Water Bay, Kowloon, Hong Kong, China

Ping Sheng

Department of Physics and Institute of Nano Science and Technology, The Hong Kong University of Science and Technology,

Clear Water Bay, Kowloon, Hong Kong, China

共Received 23 October 2002; revised manuscript received 2 May 2003; published 17 July 2003兲

From extensive molecular dynamics simulations on immiscible two-phase ﬂows, we ﬁnd the relative slip-

ping between the ﬂuids and the solid wall everywhere to follow the generalized Navier boundary condition, in

which the amount of slipping is proportional to the sum of tangential viscous stress and the uncompensated

Young stress. The latter arises from the deviation of the ﬂuid-ﬂuid interface from its static conﬁguration. We

give a continuum formulation of the immiscible ﬂow hydrodynamics, comprising the generalized Navier

boundary condition, the Navier-Stokes equation, and the Cahn-Hilliard interfacial free energy. Our hydrody-

namic model yields interfacial and velocity proﬁles matching those from the molecular dynamics simulations

at the molecular-scale vicinity of the contact line. In particular, the behavior at high capillary numbers, leading

to the breakup of the ﬂuid-ﬂuid interface, is accurately predicted.

DOI: 10.1103/PhysRevE.68.016306 PACS number共s兲: 47.11.⫹j, 68.08.⫺p, 83.10.Mj, 83.10.Ff

I. INTRODUCTION

Immiscible two-phase ﬂow in the vicinity of the contact

line 共CL兲, where the ﬂuid-ﬂuid interface intersects the solid

wall, is a classical problem that falls beyond the framework

of conventional hydrodynamics 关1–12兴. In particular, mo-

lecular dynamics 共MD兲 studies have shown relative slipping

between the ﬂuids and the wall, in violation of the no-slip

boundary condition 关6,7兴. There have been numerous ad-hoc

models 关1,8,10–12兴 to address this phenomenon, but none

was able to give a quantitative account of the MD slip ve-

locity proﬁle in the molecular-scale vicinity of the CL. While

away from the moving CL the small amount of relative slip-

ping was found to follow the Navier boundary condition

共NBC兲关13兴, i.e., relative slipping proportional to the tangen-

tial viscous stress, in the molecular-scale vicinity of the CL

the NBC failed totally to account for the near-complete slip.

This failure casts doubts on the general applicability of the

NBC to immiscible ﬂows and hinders a continuum formula-

tion of the hydrodynamics in the CL region. In particular, a

共possible兲 breakdown in the hydrodynamic description for

the molecular-scale CL region has been suggested 关7兴.In

another approach 关14兴, it was shown that the MD results can

be reproduced by continuum ﬁnite element simulations, pro-

vided the slip proﬁle extracted from MD is used as input.

This work demonstrated the feasibility of the hybrid algo-

rithm, but left unresolved the problem concerning the bound-

ary condition governing the CL motion. Without a continuum

hydrodynamic formulation, it becomes difﬁcult or impos-

sible to have realistic simulations of micro- or nanoﬂuidics,

or of immiscible ﬂows in porous media where the relative

wetting characteristics, the moving CL dissipation, and be-

havior over undulating solid surfaces may have macroscopic

implications.

From MD simulations on immiscible two-phase ﬂows, we

report the ﬁnding that the generalized Navier boundary con-

dition 共GNBC兲 applies for all boundary regions, whereby the

relative slipping is proportional to the sum of tangential vis-

cous stress and the uncompensated Young stress. The latter

arises from the deviation of the ﬂuid-ﬂuid interface from its

static conﬁguration 关10兴. By combining GNBC with the

Cahn-Hilliard 共CH兲 hydrodynamic formulation of two-phase

ﬂow 关11,12兴, we obtained a consistent, continuum descrip-

tion of immiscible ﬂow with material parameters 共such as

viscosity, interfacial tension, etc兲 directly obtainable from

MD simulations. The convective-diffusive dynamics in the

vicinity of the interface and the moving CL also means the

introduction of two phenomenological dynamic parameters

whose values can be ﬁxed by comparison with one MD ﬂow

proﬁle. Once the parameter values are determined from MD

simulations, our continuum hydrodynamics can yield predic-

tions matching those from MD simulations 共for different

Couette and Poiseuille ﬂows兲. Our ﬁndings suggest the no-

slip boundary condition to be an approximation to the

GNBC, accurate for most macroscopic ﬂows but failing in

immiscible ﬂows. These results open the door to efﬁcient

simulations of nano- or microﬂuidics involving immiscible

components, as well as to macroscopic immiscible ﬂow cal-

culations, e.g., in porous media, that are physically meaning-

ful at the molecular level 关15兴. The latter is possible, for

example, by employing the adaptive method based on the

iterative grid redistribution introduced in Ref. 关15兴. This

Author to whom correspondence should be addressed. Email ad-

dress: sheng@ust.hk

PHYSICAL REVIEW E 68, 016306 共2003兲

method has demonstrated the capability of resolving, at the

same time, both the global behavior of a partial differential

equation solution with coarse mesh and a strong singularity

in a localized region with a reﬁned local mesh of over 10

ratio to the coarse mesh.

II. MOLECULAR DYNAMICS SIMULATIONS

The MD simulations were performed for both static and

dynamic conﬁgurations in the Couette and Poiseuille ﬂows.

The two immiscible ﬂuids were conﬁned between two paral-

lel walls separated along the z direction, with the ﬂuid-solid

boundaries deﬁned by z⫽ 0,H 共see Fig. 1 for the Couette

geometry兲. Interaction between the ﬂuid molecules was mod-

eled by a modiﬁed Lennard-Jones 共LJ兲 potential U

⫽ 4

⑀

关

(

␴

/r)

⫺

␦

(

␴

/r)

兴

, where r is the distance between

the molecules,

⑀

and

␴

are the energy scale and the range of

interaction, respectively, and

␦

⫽ 1 for like molecules and

␦

⫽⫺1 for molecules of different species. Each of the two

walls was constructed by two 共or more兲关001兴 planes of an

fcc lattice 共see Appendix A兲, with each wall molecule at-

tached to the lattice site by a harmonic spring. The mean-

squared displacement of wall molecules was controlled to

obey the Lindemann criterion. The wall-ﬂuid interaction was

also modeled by a LJ potential U

, with energy and range

parameters

⑀

⫽ 1.16

⑀

and

␴

⫽ 1.04

␴

, and a

␦

for speci-

fying the wetting property of the ﬂuid. Both U

and U

were cut off at r

⫽ 2.5

␴

. The mass of the wall molecule was

set equal to that of the ﬂuid molecule m, and the average

number densities for the ﬂuids and wall were set at

␳

⫽ 0.81/

␴

and

␳

⫽ 1.86/

␴

, respectively. The temperature

was controlled at 2.8

⑀

, where k

is Boltzmann’s con-

stant. Moving the top and bottom walls at a constant speed V

in the ⫾ x directions, respectively, induced the Couette ﬂow

关7兴. Applying a body force mg

ext

to each ﬂuid molecule in

the x direction induced the Poiseuille ﬂow 关6兴. Periodic

boundary conditions were imposed on the x and y boundaries

of the sample. Most of our MD simulations were carried out

on samples consisting 6144 atoms for each ﬂuid and 2880

atoms for each wall. The sample is 163.5

␴

by 6.8

␴

along the

x and y, respectively, and H⫽ 13.6

␴

. Our MD results repre-

sent time averages over 20–40 million time steps. For tech-

nical details of our MD simulations, we followed those de-

scribed in Ref. 关16兴.

Two different cases were considered in our simulations.

The symmetric case refers to identical wall-ﬂuid interactions

for the two ﬂuids 共both

␦

⫽ 1), which leads to a ﬂat static

interface in the yz plane with a 90° contact angle. The asym-

metric case refers to different wall-ﬂuid interactions, with

␦

⫽ 1 for one and

␦

⫽ 0.7 for the other. The resulting

static interface is a circular arc with a 64° contact angle. We

measured six quantities in the Couette-ﬂow steady states of

V⫽ 0.25(

⑀

/m)

1/2

, H⫽ 13.6

␴

for the symmetric case and V

⫽ 0.2(

⑀

/m)

1/2

, H⫽13.6

␴

for the asymmetric case:

slip

, the

slip velocity relative to the moving wall; G

, the tangential

force per unit area exerted by the wall; the

␴

compo-

nents of the ﬂuid stress tensor (n denotes the outward surface

normal兲, and

We denote the region within 0.85

␴

⫽ z

of the wall the

boundary layer 共BL兲. It must be thin enough to render sufﬁ-

cient precision for measuring

slip

, while thick enough to

fully account for the tangential wall-ﬂuid interaction force,

due to the ﬁnite range of the LJ interaction. Thus, it is not

possible to do MD measurements strictly at the ﬂuid-solid

boundary, not only because of poor statistics, but also be-

cause of this intrinsic limitation. The wall force can be iden-

tiﬁed by separating the force on each ﬂuid molecule into

wall-ﬂuid and ﬂuid-ﬂuid components. For 0⬍ z⭐z

the ﬂuid

molecules can detect the atomic structure of the wall. When

coupled with kinetic collisions with the wall molecules, there

arises a nonzero tangential wall force that varies along the z

direction and saturates at z⯝z

. G

is the saturated total

tangential wall force per unit wall area 共Fig. 2兲. In Appendix

A we give account of our MD results on both the tangential

and normal components of the wall force, plus the effect共s兲

of increasing the wall thickness in our simulations from two

layers of wall molecules to four layers and to inﬁnite layers

共by using the continuum approximation beyond the four lay-

ers兲.

Spatial resolution along the x and z directions was

achieved by evenly dividing the sampling region into bins,

each ⌬x⫽ 0.425

␴

by ⌬z⫽ 0.85

␴

in size.

slip

was obtained

as the time average of ﬂuid molecules’ velocities inside the

BL, measured with respect to the moving wall; G

was ob-

tained from the time average of the total tangential wall force

experienced by the ﬂuid molecules in the BL, divided by the

bin area in the xy plane;

␴

xx(nx)

was obtained from the time

averages of the kinetic momentum transfer plus the ﬂuid-

ﬂuid interaction forces across the constant-x(z) bin surfaces,

FIG. 1. 共Color兲 Segments of the MD simulation sample for the

immiscible Couette ﬂows. The colored dots indicate the instanta-

neous molecular positions of the two ﬂuids projected onto the xz

plane. The black 共gray兲 circles denote the wall molecules. The up-

per panel illustrates the symmetric case; the lower panel illustrates

the asymmetric case. The red circles and the blue squares represent

the time-averaged interface proﬁles, deﬁned by

␳

⫽

␳

(

␾

⫽ 0), for

the two cases. The black solid lines are the interface proﬁles calcu-

lated from the continuum hydrodynamic model with the GNBC.

QIAN, WANG, AND SHENG PHYSICAL REVIEW E 68, 016306 共2003兲

016306-2

and

x(z)

was measured as the time-averaged velocity com-

ponent共s兲 within each bin. For the contribution of intermo-

lecular forces to the stress, we have directly measured the

ﬂuid-ﬂuid interaction forces across bin surfaces instead of

using the Irving-Kirkwood expression 关17兴, whose validity

was noted to be not justiﬁed at the ﬂuid-ﬂuid or the wall-

ﬂuid interface 关see the paragraph following Eq. 共5.15兲 in the

above reference兴. In Appendix B we give some details on our

MD stress measurements. As reference quantities, we also

measured G

␴

in the static (V⫽ 0) conﬁguration.

In addition, we measured in both the static and dynamic

conﬁgurations the average molecular densities

␳

and

␳

for

the two ﬂuid species in each bin to determine the interface

proﬁle. The shear viscosity

␩

⫽ 1.95

冑

⑀

␴

and the interfa-

cial tension

␥

⫽ 5.5

⑀

␴

were also determined.

We have also measured the interface and velocity proﬁles

for the Poiseuille ﬂow in the asymmetric case, as well as for

the Couette ﬂows with different V and H in the symmetric

case.

III. GENERALIZED NAVIER BOUNDARY CONDITION

In the presence of a ﬂuid-ﬂuid interface, the static ﬂuid

stress tensor

␴

reﬂects the static Young stress 共surface ten-

sion兲 as well as those stresses arising from wall-ﬂuid inter-

action. This is the case in spite of the fact that in all the MD

ﬂuid stress measurements only the ﬂuid-ﬂuid interaction was

counted 共see Appendix B 2兲. The reason is that because the

MD measurements were carried out either in the static equi-

librium state or in the dynamic steady state, local force bal-

ance necessarily requires the ﬂuid stress components to fully

reﬂect the inﬂuence of the wall-ﬂuid interaction. For the con-

sideration of moving CL, we will be concerned with the part

of the ﬂuid stress tensor, which is purely dynamic in origin,

i.e., arising purely from the hydrodynamic motion of the

ﬂuid 共and the CL兲. In the notations below, the over tilde

denotes the difference between that quantity and its static

part. Thus, if

␴

is the total stress, we will be concerned only

with the hydrodynamic part, denoted by

␴

⫽

␴

⫺

␴

. We note

that in the absence of body forces, the momentum equation

in bulk ﬂuid is given by

␳

关

⳵

t⫹ (v•“)v

兴

⫽ “•

␴

. In the

BL, the wall-ﬂuid interaction means the existence of a dy-

namic, tangential wall force density g

such that the force

balance equation is given by (“•

␴

)•x

⫹ g

⫽ 0 inside the

BL. The tangential wall force density g

, shown explicitly

in the inset to Fig. 2, is a function sharply peaked at z

⬇z

/2. Here we note that the boundary layer thickness is

extremely small (z

⫽ 0.85

␴

), hence the inertial effect may

be neglected (m

␳

␩

⬍ 0.1). MD evidence for an inte-

grated form of the steady-state force balance is shown in Fig.

3. The total tangential force exerted by the wall on the ﬂuid

is given by G

⫽

兰

dzg

per unit wall area. In steady state,

this wall force is necessarily balanced by the tangential ﬂuid

force G

⫽

兰

dz(

⳵

␴

⫹

⳵

␴

) 共inset to Fig. 3兲关18兴. Here

⳵

x,z,n

means taking partial derivative with respect to x, z,or

surface normal.

We now present evidences to show that everywhere on the

boundaries, relative slipping is proportional to G

关the

GNBC, see also Eq. 共3兲 below兴:

⫽

␤

slip

, 共1兲

where

␤

is the slip coefﬁcient and G

can be written as

⫽

⳵

冕

␴

共

兲

⫺

␴

共

兲

, 共2兲

where we have used the fact that

␴

(0)⫽ 0. 关More strictly,

␴

⫺

)⫽ 0 because there is no ﬂuid below z⫽ 0 and hence

no momentum transport across z⫽ 0] Here the z coordinate

is for the lower ﬂuid-solid boundary 共same below兲, with the

understanding that the same physics holds at the upper

boundary, and

⳵

⫽⫺

⳵

for the lower boundary.

Force balance means that at steady state, the frictional

force exerted by the solid wall on the moving 共slipping兲 ﬂuid

is fully accounted for in G

. Thus, the GNBC 共or NBC兲 can

be expressed in either G

or G

, but not both. In Fig. 3 we

show the measured MD data for the symmetric and asym-

metric cases in the Couette geometry. The symbols represent

the values of G

measured in the BL. The solid lines repre-

sent the values of G

calculated from

␤

slip

by using

␤

⫽

␤

⫽

␤

⫽ 1.2

冑

⑀

␴

for the symmetric case and

␤

⫽ 1.2

冑

⑀

␴

␤

⫽ 0.532

冑

⑀

␴

for the asymmetric case

away from the CL region 共straight line segments in Fig. 3兲,

and

␤

⫽ (

␤

␳

⫹

␤

␳

)/(

␳

⫹

␳

) in the CL region 关19兴, with

slip

and

␳

1,2

obtained from MD simulations. It is seen that

FIG. 2. By subdividing the boundary layer into thin sections, we

plot the accumulated wall force per unit wall area as a function of

distance z away from the boundary. Here G

(z) is deﬁned by

(z)⫽

兰

⬘

), where g

is the density of tangential wall

force. For different x positions, the absolute value of the saturating

total wall force is different. However, when normalized by the cor-

responding saturated total wall force per unit area at each x, all

points fall on a universal curve, nearly independent of x. It is seen

that at z⫽ z

the wall force has reached its saturation value. Inset:

Tangential wall force density plotted as a function of distance z

away from the boundary. The solid lines are averaged g

in thin

sections at different x, normalized by the corresponding saturated

total wall force per unit area. The dashed line is a smooth Gaussian

ﬁt. It is seen that g

is a function sharply peaked at z⬇z

/2. In the

sharp boundary limit this peaked wall force density is approximated

by G

␦

(z).

MOLECULAR SCALE CONTACT LINE HYDRODYNAMICS... PHYSICAL REVIEW E 68, 016306 共2003兲

016306-3

for the lower boundary 共upper right panel兲, the MD data

agree well with the predictions of Eq. 共1兲. For the upper

boundary 共lower left panel兲 the straight line segments also

agree well with Eq. 共1兲. However, there is some discrepancy

in the interfacial region of the upper boundary that seems to

arise from a ‘‘shear thinning’’ effect of decreasing

␤

at very

large tangential stresses 关13兴.

The fact that the wall force density is distributed inside a

thin BL and vanishes beyond the BL necessitates the form of

as deﬁned by Eq. 共2兲. However, it is intuitively obvious

that the ﬂuids would experience almost the identical physical

effect共s兲 from a wall force density G

␦

(z), concentrated

strictly at the ﬂuid-solid boundary with the same total wall

force per unit area. In the inset to Fig. 2, it is shown that the

MD-measured wall force density is a sharply peaked func-

tion. The sharp boundary limit involves the approximation of

replacing this peaked function by

␦

(z). The replacement of a

diffuse boundary by a sharp boundary can considerably sim-

plify the form of the GNBC, because local force balance

along x then requires

⳵

␴

⫹

⳵

␴

⫽ 0 away from the

boundary z⫽ 0. Integration of this relation from 0

⫹

to z

yields

⳵

冕

␴

共

兲

⫹

␴

共

兲

⫺

␴

共

⫹

兲

⫽ 0

FIG. 3. 共Color兲

␤

V/G

plotted as a function of V/

slip

. Sym-

bols are MD data measured in the BL at different x locations, where

the red circles denote the symmetric case and the blue squares de-

note the asymmetric case. The solid lines were calculated from Eq.

共1兲 with values of

␤

1,2

and the expression of

␤

given in the text.

The statistical errors of the MD data are about the size of the sym-

bols. The upper-right data segment corresponds to the lower bound-

ary, whereas the lower-left segment corresponds to the upper

boundary. The slopes of the two dashed lines are given by

␤

1,2

⫺ 1

Inset: G

plotted as a function of G

, measured in the two BL’s at

different values of x. The symbols have the same correspondence as

in the main ﬁgure. The data are seen to lie on a straight line with a

slope of ⫺ 1, indicating G

⫹ G

⫽ 0.

FIG. 4. 共Color兲 Two components of the dynamic tangential

stress at z⫽ z

, plotted as a function of x. The dashed lines denote

␴

; solid lines represent the viscous component. Here red indicates

the symmetric case and blue indicates the asymmetric case. In the

CL region the nonviscous component is 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

␪

closer to

␪

. Inset: ⌺

d,s

plotted as a

function of

␥

cos

␪

d,s

at different values of z. Here ⌺

⫽⫺

兰

dx(

␴

⫺

␴

), ⌺

⫽⫺

兰

␴

, and

␪

d,s

was measured from

the time-averaged interfacial proﬁles 共Fig. 1兲. The red circles denote

the symmetric case, the blue squares denote the asymmetric case,

the solid blue squares denote the asymmetric static case, and the

single solid red circle at the origin denotes the symmetric static

case. The data are seen to follow a straight 共dashed兲 line with slope

1, indicating ⌺

d,s

⫽

␥

cos

␪

d,s

FIG. 5. 共Color兲 S⫽

兰

␴

(z)dz⫽

兰

关

␴

(z)⫺

␴

(z)

兴

dz plot-

ted as a function of x. Here red circles denote the symmetric case

and blue squares denote the asymmetric case. For clarity,

␴

was

vertically displaced such that

␴

⫽ 0 far from the interface in the

symmetric case, and for the asymmetric case,

␴

⫽ 0 at the center

of the interface.

QIAN, WANG, AND SHENG PHYSICAL REVIEW E 68, 016306 共2003兲

016306-4

and as a consequence 关by comparing with Eq. 共2兲兴 G

⫽⫺

␴

⫹

). Therefore,

␴

changes from

␴

⫺

)⫽ 0to

␴

⫹

)⫽ G

at z⫽ 0, leading to

共

“•

␴

兲

•x

⫽ G

␦

共

兲

Comparing with the diffuse boundary, where (“•

␴

)•x

⫹ g

⫽ 0, we see that the form of the equation remains the same,

but the BL is now from 0

⫺

to 0

⫹

, instead of from 0

⫹

to z

as in the diffuse case. Thus, GNBC 共1兲 becomes

⫺

␴

共

兲

⫽

␤

slip

in the sharp boundary limit.

The tangential stress

␴

can be decomposed into a vis-

cous component and a non-viscous component:

␴

共

兲

⫽

␴

共

兲

⫹

␴

共

兲

In Fig. 4 we show that away from the interfacial region the

tangential viscous stress

␴

(z)⫽

␩

(

⳵

⫹

⳵

)(z) is the

only nonzero component, but in the interfacial region

␴

⫽

␴

⫺

␴

⫺

␴

⫽

␴

⫺

␴

is dominant, thereby account-

ing for the failure of NBC to describe the CL motion. There-

fore, away from the CL region the NBC is valid, but in the

interfacial region the NBC clearly fails to describe the CL

motion. We wish to clarify the origin of

␴

and

␴

as the

dynamic and static Young stresses, respectively, so that

␴

⫽

␴

⫺

␴

is the uncompensated Young stress. As shown in

the inset to Fig. 4, the integrals 共across the interface兲 of

␴

(⫽

␴

⫺

␴

, calculated by subtracting the viscous compo-

nent

␩

(

⳵

⫹

⳵

) from the total tangential stress

␴

) and

␴

are equal to

␥

cos

␪

and

␥

cos

␪

, respectively, at differ-

ent values of z, i.e.,

⫺

冕

int

␴

共

兲

⫽

␥

cos

␪

共

兲

and

⫺

冕

int

␴

共

兲

⫽

␥

cos

␪

共

兲

where

␪

(z) and

␪

(z) are, respectively, the dynamic and the

static interfacial angles at z 关20兴. Here

兰

int

dx denotes the

integration across the ﬂuid-ﬂuid interface along x. These re-

sults clearly show the origin of the extra tangential stress in

the interfacial region to be the interfacial 共uncompensated兲

Young stress. Thus, the GNBC is given by

␤

slip

⫽⫺

␴

共

兲

⫽⫺

关

␩

⳵

兴共

兲

⫺

␴

共

兲

. 共3兲

Here only one component of the viscous stress is nonzero,

due to

⫽ 0 at the boundary, and ⫺

␴

(0) is the uncom-

pensated Young stress, satisfying

⫺

冕

int

␴

共

兲

dx⫽

␥

共

cos

␪

surf

⫺ cos

␪

surf

兲

with

␪

d(s)

surf

being a microscopic dynamic 共static兲 contact angle

at the ﬂuid-solid boundary. The fact that

␴

(0)⬇0 away

from the CL shows that the GNBC implies NBC for single

phase ﬂows.

Due to the diffuse nature of the BL in the MD simula-

tions, the contact angle

␪

d(s)

surf

cannot be directly measured.

Nevertheless, they are obtainable through extrapolation by

using the integrated interfacial curvature within the BL. That

is, in the sharp boundary limit the force balance in the ﬂuids

is expressed by

⳵

␴

⫹

⳵

␴

⫽ 0. Integration in z across the

BL gives

⳵

冕

␴

共

兲

⫺

␴

共

兲

⫹

␴

共

兲

⫺

␴

共

兲

⫹

␴

共

兲

⫽ 0.

共4兲

Integration 关of Eq. 共4兲 along x] across the ﬂuid-ﬂuid inter-

face then yields

⌬

冋

冕

␴

共

兲

册

⫺

冕

int

␴

共

兲

⫹

冕

int

␴

共

兲

⫹

␥

⫺

␥

⫽ 0, 共5兲

where ⌬

关

兰

␴

(z)

兴

is the change of the z-integrated

␴

across the interface, K

and K

denote the dynamic and the

static z-integrated interfacial curvatures:

⫽ cos

␪

共

兲

⫺ cos

␪

surf

and

⫽ cos

␪

共

兲

⫺ cos

␪

surf

Here ⌬

关

兰

␴

(z)

兴

␴

␪

), and

␪

) are ob-

tainable from MD simulations, K

⯝⫾ 2z

cos

␪

surf

/H for the

circular static interfaces, while

␴

(0)⫽

␩

关

⳵

兴

(0) may be

obtained by extrapolating from the values of tangential vis-

cous stress at z⫽ z

,2z

, and 3z

. Therefore, the micro-

scopic dynamic contact angle

␪

surf

can be obtained from Eq.

共5兲. In Appendix B 3 we give a more detailed account of the

relationship between the MD measured stresses and the

stress components in the continuum hydrodynamics. The

above extrapolation is based on this correspondence.

We have measured the z-integrated

␴

⫽

␴

⫺

␴

in the

BL. The dominant behavior is a sharp drop across the inter-

face, as shown in Fig. 5 for both the symmetric and asym-

metric cases. The value of

␪

surf

obtained is 88°⫾ 0.5° for the

symmetric case and 63°⫾ 0.5° for the asymmetric case at the

lower boundary, and 64.5°⫾ 0.5° at the upper boundary.

These values are noted to be very close to

␪

surf

. Yet the

small difference between the dynamic and static 共micro-

scopic兲 contact angles is essential in accounting for the near-

complete slip in the CL region.

In essence, our results show that in the vicinity of the CL,

the tangential viscous stress ⫺

␴

as postulated by the NBC

cannot give rise to the near-complete CL slip without taking

into account the tangential Young stress ⫺

␴

in combina-

MOLECULAR SCALE CONTACT LINE HYDRODYNAMICS... PHYSICAL REVIEW E 68, 016306 共2003兲

016306-5

Molecular scale contact line hydrodynamics of immiscible flows.

Figures

Citations

Wetting and Spreading

Heterogeneous multiscale methods: A review

The heterogeneous multiscale method

A phase field model for the mixture of two incompressible fluids and its approximation by a Fourier-spectral method

A multi-phase SPH method for macroscopic and mesoscopic flows

References

Computer Simulation of Liquids

Related Papers (5)

Hydrodynamic Model of Steady Movement of a Solid / Liquid / Fluid Contact Line

The dynamics of the spreading of liquids on a solid surface. Part 1. Viscous flow

Wetting: statics and dynamics

Wetting and Spreading

On the Spreading of Liquids on Solid Surfaces: Static and Dynamic Contact Lines

Frequently Asked Questions (12)

Q1. What are the contributions mentioned in the paper "Molecular scale contact line hydrodynamics of immiscible flows" ?

Q2. What is the force balance equation in the BL?

Q3. How did the authors increase the wall thickness in their simulations?

Q4. What is the rate of work done to the system in the steady state?

Q5. How was the kinetic momentum transfer measured?

Q6. What is the udsur f value for a moving CL?

Q7. How is the rate of work calculated?

Q8. What is the rate of dissipation due to the moving CL?

Q9. What was the dw f for specifying the wetting property of the fluid?

Q10. What is the way to verify that the boundary conditions and the parameter values are local properties?

Q11. What is the first order approximation of the CH free energy functional?

Q12. What is the energy scale of the interaction between the two molecules?