MOLECULAR DYNAMICS AND CONTINUUM
SIMULATIONS OF FLUID FLOWS WITH SLIP
BOUNDARY CONDITIONS

By
Anoosheh Niavaranikheiri

A DISSERTATION
Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of
DOCTOR OF PHILOSOPHY

Mechanical Engineering
2011

ABSTRACT
MOLECULAR DYNAMICS AND CONTINUUM
SIMULATIONS OF FLUID FLOWS WITH SLIP
BOUNDARY CONDITIONS
By
Anoosheh Niavaranikheiri

Microﬂuidics is a rapidly developing ﬁeld with applications ranging from molecular biology, environmental monitoring, and clinical diagnostics. Microﬂuidic systems are characterized by large surface-to-volume ratios, and, therefore, ﬂuid ﬂows are signiﬁcantly inﬂuenced
by boundary conditions. The fundamental assumption in ﬂuid mechanics is the no-slip
boundary condition, which states that the tangential ﬂuid velocity is equal to the adjacent
wall speed. Although this assumption is successful in describing ﬂuid ﬂows on macroscopic
length scales, recent experimental and numerical studies have shown that it breaks down at
microscopic scales due to the possibility of slip of the ﬂuid relative to the wall. The eﬀect of
slip is more pronounced for highly viscous liquids like polymer melts or in the region near
the moving contact line due to the large gradient in shear stress at the liquid/solid interface.
The measure of slip is the so-called slip length, which is deﬁned as a distance between the
real interface and imaginary plane where the extrapolated velocity proﬁle vanishes. The slip
length value is sensitive to several key parameters, such as surface energy, surface roughness,
ﬂuid structure, and shear rate.
In this dissertation, the slip phenomena in thin liquid ﬁlms conﬁned by either ﬂat or structured surfaces are investigated by molecular dynamics (MD) and continuum simulations. It
is found that for ﬂows of both monatomic and polymeric ﬂuids over smooth surfaces, the
slip length depends nonlinearly on shear rate at suﬃciently high rates. The laminar ﬂow

away from a curved boundary is usually described by means of the eﬀective slip length,
which is deﬁned with respect to the mean roughness height. MD simulations show that for
corrugated surfaces with wavelength larger than the size of polymer chains, the eﬀective
slip length decreases monotonically with increasing corrugation amplitude. A detailed comparison between the solution of the Navier-Stokes equation with the local rate-dependent
slip condition and results of MD simulations indicates that there is excellent agreement
between the velocity proﬁles and the eﬀective slip lengths at low shear rate and small-scale
surface roughness. It was found that the main cause of the slight discrepancy between MD
and continuum results at high shear rates is the reduction of the local slip length in the
higher pressure regions where the boundary slope becomes relatively large with respect to
the mainstream ﬂow. It was further shown that for the Stokes ﬂow with the local no-slip
boundary condition, the eﬀective slip length decreases with increasing corrugation amplitude and a ﬂow circulation is developed in suﬃciently deep grooves. Analysis of a numerical
solution of the Navier-Stokes equation with the local slip condition shows that the inertial
eﬀects promote the asymmetric vortex ﬂow formation and reduce the eﬀective slip length.

TABLE OF CONTENTS

LIST OF TABLES

vi

LIST OF FIGURES

vii

1 Introduction

1

2 Rheological study of polymer ﬂow past rough surfaces with slip boundary
conditions
2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 The details of the numerical simulations . . . . . . . . . . . . . . . . . . . .
2.2.1 Molecular dynamics model . . . . . . . . . . . . . . . . . . . . . . . .
2.2.2 Continuum method . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 MD results for ﬂat walls . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Results for periodically corrugated walls: large wavelength . . . . . . . . . .
2.4.1 MD simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4.2 Comparison between MD and continuum simulations . . . . . . . . .
2.5 Results for periodically corrugated walls: small wavelengths . . . . . . . . .
2.5.1 Comparison between MD and continuum simulations . . . . . . . . .
2.5.2 The polymer chain conﬁguration and dynamics near rough surfaces .
2.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4
4
7
7
8
10
11
11
12
14
14
15
17

3 Slip boundary conditions for shear ﬂow of polymer melts past atomically
ﬂat surfaces
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Details of molecular dynamics simulations . . . . . . . . . . . . . . . . . . .
3.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Fluid density and velocity proﬁles . . . . . . . . . . . . . . . . . . . .
3.3.2 Rate dependence of the slip length and shear viscosity . . . . . . . .
3.3.3 Analysis of the ﬂuid structure in the ﬁrst layer . . . . . . . . . . . . .
3.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31
31
33
36
36
36
39
41

4 The eﬀective slip length and vortex formation in laminar ﬂow over a rough
surface
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Details of numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3.1 Analytical solution of the Stokes equation for viscous ﬂow over a wavy
wall . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
iv

54
54
57
60
60

4.4

4.3.2 Flow over a rough surface with the local no-slip boundary condition
4.3.3 Eﬀect of the local slip on the ﬂow pattern near the rough surface .
4.3.4 Eﬀect of the Reynolds number on the eﬀective slip length . . . . . .
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 Modeling the combined eﬀect of surface roughness and shear
ﬂow of simple ﬂuids
5.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.2 The details of numerical methods . . . . . . . . . . . . . . . . .
5.2.1 Molecular dynamics simulations . . . . . . . . . . . . . .
5.2.2 Continuum simulations . . . . . . . . . . . . . . . . . . .
5.3 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.3.1 MD simulations: the intrinsic and eﬀective slip lengths .
5.3.2 Comparison between MD and continuum simulations . .
5.3.3 A detailed analysis of the ﬂow near the curved boundary
5.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.
.
.
.

60
63
65
67

rate on slip
.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

.
.
.
.
.
.
.
.
.

86
. 86
. 89
. 89
. 91
. 93
. 93
. 95
. 96
. 100

6 Slip boundary conditions for the moving contact line in molecular dynamics and continuum simulations
116
6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116
6.2 Molecular dynamics (MD) simulations . . . . . . . . . . . . . . . . . . . . . 118
6.2.1 Equilibrium case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 118
6.2.2 Shear ﬂow . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.3 Details of continuum simulations . . . . . . . . . . . . . . . . . . . . . . . . 123
6.3.1 Fixed grid method for immersed objects . . . . . . . . . . . . . . . . 123
6.3.2 Continuum Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
BIBLIOGRAPHY

145

v

LIST OF TABLES
2.1

ˆ ˆ
ˆ
Averaged x , y , and z components of the radius of gyration at equilibrium
and in the shear flow. The Rg values are reported in the bulk, near the flat

upper wall, above the peaks, and inside the grooves. The wavelength of the
lower wall is λ = 7.5σ . The size of the averaging region above the peaks
and inside the grooves is σ × 12.5σ × 1.5σ . The estimate of the error bars is
0.03σ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
2.2

ˆ ˆ
ˆ
Averaged x , y , and z components of the radius of gyration at equilibrium
and in the shear flow. The Rg values are reported in the bulk, near the flat

upper wall, above the peaks, and inside the grooves. The wavelength of the
lower is λ = 3.75σ . The dimensions of the averaging region are the same as
in the Table 2.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.1

The fluid pressure P at equilibrium, i.e., with the stationary upper wall, and
the channel height h (defined as a distance between the inner FCC planes)
as a function of the fluid density. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34

vi

LIST OF FIGURES

2.1

2.2

2.3

2.4

2.5

2.6

2.7

2.8

For interpretation of the references to color in this and all other ﬁgures,
the reader is referred to the electronic version of this dissertation. A
snapshot of the ﬂuid monomers (open circles) conﬁned between solid
walls (closed circles) obtained from the MD simulations. The atoms
of the stationary lower wall are distributed along the sinusoidal curve
with the wavelength λ and amplitude a. The ﬂat upper wall is moving
with a constant velocity U in the x direction. The eﬀective slip length
ˆ
Leﬀ is determined by the linear extrapolation of the velocity proﬁle to
ux = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

20
18

Averaged velocity proﬁles in the cell with ﬂat upper and lower walls.
The solid lines are the best linear ﬁt to the data. The vertical axes
indicate the location of the F CC lattice planes. The velocities of the
upper wall are tabulated in the inset. . . . . . . . . . . . . . . . . . .

21
19

Viscosity of the polymer melt µ/ετ σ −3 as a function of shear rate. The
dashed line with the slope −0.33 is plotted for reference. . . . . . . .

22
20

Variation of the slip length as a function of shear rate in the cell with
ﬂat upper and lower walls. The solid line is a fourth-order polynomial
ﬁt to the data given by Eq. (2.9). . . . . . . . . . . . . . . . . . . . .

21
23

Averaged density proﬁles near the stationary upper wall (a), above the
peak and in the groove of the lower wall with amplitude a = 8.16 σ and
wavelength λ = 66.60 σ (b). . . . . . . . . . . . . . . . . . . . . . . . .

24
22

Averaged velocity proﬁles for the indicated values of the corrugation
amplitude a. The vertical dashed line denotes a reference plane for
calculation of the eﬀective slip length at the corrugated lower wall.
The velocity of the ﬂat upper wall is U = 0.5 σ/τ . . . . . . . . . . . .

23
25

The eﬀective slip length as a function of wavenumber ka obtained from
the MD simulations ( ), the solution of the Stokes equation with rateindependent slip length L0 (◦) and with L0 (γτ ) given by Eq. (2.9) (×),
˙
the solution of the Navier-Stokes equation with L0 ( ). . . . . . . . .

26
24

Averaged ﬂuid density proﬁles near ﬂat upper wall (a) and corrugated
lower wall with amplitude a = 1.4 σ and wavelength λ = 7.5 σ (b). The
velocity of the upper wall is U = 0.5 σ/τ . . . . . . . . . . . . . . . . .

27
25

vii

2.9

The eﬀective slip length as a function of wavenumber ka for the indicated values of wavelength λ. Continuum results are denoted by the
dashed lines and open symbols, while the MD results are shown by
straight lines and ﬁlled symbols. The inset shows the same data for
ka 0.25. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28
26

2.10 A snapshot of four polymer chains near the lower corrugated wall for
wavelength λ = 7.5 σ and amplitude a = 1.4 σ. The velocity of the
upper wall is U = 0.5 σ/τ . . . . . . . . . . . . . . . . . . . . . . . . .

29

2.11 A snapshot of ﬁve polymer chains near the lower corrugated wall for
the wavelength λ = 3.75 σ and amplitude a = 1.0 σ. The ﬁgure shows
the side view (top) and the top view (bottom). The velocity of the
upper wall is U = 0.5 σ/τ . . . . . . . . . . . . . . . . . . . . . . . . .

30
30

3.1

3.2

3.3

3.4

3.5

3.6

A schematic representation of the steady-state shear ﬂow in the Couette cell. The upper wall is translated with a constant velocity U in
the x direction. The ﬂuid slip velocity is determined from the relation
ˆ
Vs = γLs , where γ is the slope of the velocity proﬁle and Ls is the slip
˙
˙
length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43
43

Averaged ﬂuid density proﬁles near the stationary thermal wall with
κ = 1200 ε/σ 2 and εwf /ε = 0.9. The velocities of the upper wall are
tabulated in the inset. The ﬂuid density in the center of the cell is
ρ = 0.91 σ −3 . The left vertical axis denotes the position of the liquid/solid interface at z = − 14.74 σ. . . . . . . . . . . . . . . . . . . .

44
44

Average normalized velocity proﬁles, Vx / U , for the indicated values of
the upper wall velocity U and the ﬂuid density ρ = 0.91 σ −3 . Vertical
axes denote the position of the F CC lattice planes at z = − 15.24 σ
and z = 12.03 σ. The dashed lines indicate the location of liquid/solid
interfaces at a distance 0.5 σw away from the F CC lattice planes, i.e.,
z = − 14.74 σ and z = 11.53 σ. The inset shows the same data near the
lower wall. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

45

Behavior of the ﬂuid viscosity µ/ ετ σ −3 as a function of shear rate for
the tabulated values of the polymer density. The straight dashed line
with a slope −0.37 is shown for reference. Solid curves are a guide for
the eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

46

Slip length Ls /σ as a function of the shear rate for the polymer melt
with linear chains N = 20. Solid curves are a guide for the eye. . . . .

46
47

Log-log plot of the friction coeﬃcient, k = µ/Ls , as a function of the
slip velocity for the indicated values of the ﬂuid density. The dashed
curve is the best ﬁt to Eq. (3.10) with C1 = 1.77 mσ −3 and C2 =
0.71 σ 2 τ −2 . Solid curves are a guide for the eye. . . . . . . . . . . . .

47
48

viii

3.7

Master curve for the friction coeﬃcient as a function of the slip velocity.
The same data as in Fig. 3.6. The values of the friction coeﬃcient
k ∗ [ ετ σ −4 ] and the slip velocity Vs∗ [ σ/τ ] used for normalization of
the data are shown in the inset. The dashed line y = [1 + (4x)2 ]−0.35
is plotted for reference. . . . . . . . . . . . . . . . . . . . . . . . . . .

48
49

Structure factor S(kx , ky ) in the ﬁrst ﬂuid layer for the ﬂuid density ρ = 0.91 σ −3 . The upper wall velocity is U = 0.05 σ/τ (a) and
U = 5.0 σ/τ (b). The same ﬂow conditions as in Fig. 3.3. . . . . . . .

49
50

Temperature proﬁles for the indicated values of the upper wall speed
U and the ﬂuid density ρ = 0.91 σ −3 . Vertical dashed lines denote the
position of the liquid/solid interface. The inset shows the x, y and z
ˆ ˆ
ˆ
components of the temperature proﬁle near the stationary lower wall
when the velocity of the upper wall is U = 5.0 σ/τ . . . . . . . . . . .

50
51

3.10 Log-log plot of Ls /µ as a function of the ratio T1 kB /ε [S(G1 ) ρc σ 3 ]−1
for the indicated values of the polymer density. The solid line with a
slope 0.9 is plotted for reference. . . . . . . . . . . . . . . . . . . . . .

51
52

3.8

3.9

3.11 Log-log plot of the ratio Ls /µ as a function of the variable S(0)/ [S(G1 ) ρc σ 3 ].
The values of the ﬂuid density are tabulated in the inset. The same
data as in Fig. 3.10. The straight solid line with a slope 1.15 is plotted
as a reference. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
52
4.1

4.2

4.3

4.4

Schematic diagram of the steady-state Couette ﬂow over a rough surface. The upper ﬂat wall is moving with a constant velocity U in the
x direction. The lower stationary wall is modeled as a sinusoidal wave
ˆ
with amplitude a and wavelength λ. The wavenumber ka = 2 πa/λ
varies in the range 0 ka 1.26. . . . . . . . . . . . . . . . . . . . .

68

Diagram of a bilinear element in (a) the physical coordinate system
(x, z) and (b) a transformed element in the natural coordinate system
(ξ, η). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

68
69

The eﬀective slip length as a function of wavenumber ka computed
from the solution of the Stokes equation with the no-slip boundary
condition. The primary vortex is formed at the bottom of the valley
for ka 0.79 (see the vertical arrow at ka = 0.79). The solid line is
calculated using Eq. (4.18). The inset shows the normalized velocity
proﬁles obtained from the Stokes solution for the selected values of ka.
The dashed line located at a = 0 is the reference for calculating the
eﬀective slip length. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

69
70

Shear stress along the lower wavy wall computed from the Stokes solution with no-slip boundary condition (L0 = 0) for the indicated values
∗
of wavenumber ka. The value τw used for normalization is the shear
stress at the ﬂat wall. The intersection of the curves with the dashed
line determines the location of the ﬂow separation and attachment inside the valley. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70
71

ix

4.5

The normalized pressure along the lower corrugated boundary extracted from the Stokes solution with no-slip boundary condition for
the same values of ka as in Fig. 4.4. The maximum pressure P ∗ for
ka = 1.12 and L0 = 0 is located in the bulk region at (x/λ, z/λ)
(0.07, 0.23). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71
72

Pressure contours and streamlines for the wavenumber ka = 1.12
and the slip length at the lower wall L0 = 0 (a), L0 /h = 0.03 (b),
L0 /h = 0.06 (c), and L0 /h = 0.08 (d). The pressure contours are
normalized by the maximum value P ∗ located on the left side of the
peak (x/λ, z/λ)
(0.07, 0.23) in the case of L0 = 0. The horizontal
dashed line denotes the location of the eﬀective no-slip boundary plane.
The vertical dashed line inside the valley at x/λ = 0.75 indicates the
cross-section used to compute the velocity proﬁles shown in the inset.

72
73

Normalized pressure along the lower wavy wall extracted from the
solution of the Stokes equation for the tabulated values of L0 and
ka = 1.12. The maximum pressure P ∗ for L0 = 0 is located above the
lower surface at (x/λ, z/λ) (0.07, 0.23). . . . . . . . . . . . . . . . .

76
77

Shear stress along the corrugated lower wall (ka = 1.12) for the indicated values of the local slip length L0 . The intersection of the dashed
line with the shear stress proﬁles shows the location of the ﬂow separation and attachment inside the valley. . . . . . . . . . . . . . . . . .

78
77

The eﬀective slip length, extracted from the Stokes solution as a function of L0 for ka = 0.28 (top) and ka = 1.12 (bottom). The dashed
line is computed from Eq. (4.18). For ka = 1.12 and L0 /h 0.067, the
vortex is formed in the valley (see the vertical arrow at L0 /h = 0.067).
The error bars are smaller than the symbol size. . . . . . . . . . . . .

78
79

4.10 The intrinsic slip length L0 above which there is no vortex at the
bottom of the valley (top) and the corresponding eﬀective slip length
Leﬀ (bottom) computed from the Stokes solution and plotted as a
function of wavenumber ka. The increment in the slip length used to
determine the threshold values is about the symbol size (top). The
dashed lines are quadratic ﬁts to guide the eye. . . . . . . . . . . . .

79
80

4.11 The normalized pressure (top) and shear stress (bottom) along the
lower wavy wall as a function of the Reynolds number for ka = 1.12
and no-slip boundary conditions. The value P ∗ is the maximum pressure located above the lower boundary on the left side of the peak
(x/λ, z/λ) (0.07, 0.23) for L0 = 0 and Re = 4.0. . . . . . . . . . . .

80
81

4.12 Pressure contours and streamlines extracted from the solution of the
Navier-Stokes equation for ka = 1.12 and L0 = 0 (a) and L0 /h = 0.25
(b). The value P ∗ is the maximum pressure located at (x/λ, z/λ)
(0.07, 0.23) for ka = 1.12, L0 = 0, and Re = 4.0. The horizontal
dashed line denotes the eﬀective no-slip boundary plane. The vertical
dashed line inside the valley indicates the cross-section used to compute
the velocity proﬁles shown in the inset. . . . . . . . . . . . . . . . . .

81
82

4.6

4.7

4.8

4.9

x

4.13 The normalized pressure (top) and shear stress (bottom) along the
lower wavy wall as a function of the Reynolds number for ka = 1.12 and
L0 /h = 0.25. The normalization value P ∗ is the maximum pressure
located above the surface on the left side of the peak (x/λ, z/λ)
(0.07, 0.23) for L0 = 0 and Re = 4.0. . . . . . . . . . . . . . . . . . .

83
84

4.14 The eﬀective slip length computed from the solution of the NavierStokes equation as a function of the Reynolds number for ka = 1.12
and L0 = 0 ( ), L0 /h = 0.09 ( ), L0 /h = 0.12 ( ), and L0 /h = 0.25
( ). The dashed line indicates the upper bound of the region where
a vortex is present in the groove of a rough surface. The error bars
associated with the threshold values of the Reynolds number are about
the symbol size. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

85
84

5.1

The projection of x and z coordinates of the ﬂuid monomers (open
circles) conﬁned between solid walls (ﬁlled circles). The atoms of the
lower stationary wall are uniformly distributed along the sinusoidal
curve with wavelength λ, amplitude a, and density ρw = 3.1 σ −3 . The
shear ﬂow is induced by the ﬂat upper wall with the density 0.67 σ −3
moving with a constant velocity U in the x direction. The eﬀective
ˆ
slip length Leﬀ is computed by extrapolation of the velocity proﬁle to
ux = 0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5.2

Averaged velocity proﬁles for the tabulated upper wall speeds in the
cell with ﬂat upper and lower walls. The solid lines are the best linear
ﬁts to the data. The vertical axes indicate the location of the F CC
lattice planes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

5.3

Shear viscosity in the cell with ﬂat upper and lower walls when the temperature of the Langevin thermostat is set to T = 1.1 ε/kB ( ). Fluid
viscosity as a function of the thermostat temperature 1.1 T kB /ε
1.4 for the upper wall speed U = 1.0 σ/τ ( ) and 1.1 T kB /ε 1.35
for U = 3.5 σ/τ (♦). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

5.4

The intrinsic slip length L0 /σ as a function of shear rate for atomically
ﬂat walls ( ). The blue solid curve is a ninth-order polynomial ﬁt to
the data. The eﬀective slip length Leﬀ /σ as a function of shear rate for
the ﬂow over the corrugated wall with wavenumber ka = 0.27 (◦). The
eﬀective slip length computed from the solution of the Navier-Stokes
equation with the local rate-dependent slip length (•). . . . . . . . . . 105

5.5

Temperature of the ﬁrst ﬂuid layer T (x) (in units of ε/kB ) along the
lower corrugated wall (ka = 0.27) for the indicated upper wall speeds.
The vertical dashed lines denote the location of the crest (x/λ = 0.25)
and the bottom of the valley (x/λ = 0.75). . . . . . . . . . . . . . . . 106

5.6

The normal pressure Pn (in units of ε/σ 3 ) along the lower corrugated
wall (ka = 0.27) for the selected upper wall speeds. The vertical dashed
lines indicate the location of the crest (x/λ = 0.25) and the bottom of
the valley (x/λ = 0.75). . . . . . . . . . . . . . . . . . . . . . . . . . . 107
xi

5.7

A snapshot of the ﬂuid monomers (open circles) near the lower corrugated wall (ﬁlled circles) with wavenumber ka = 0.27. The rectangular boxes denote averaging regions on the left side of the peak (I)
(x/λ
0.05), above the crest (II) (x/λ
0.25), on the right side
of the peak (III) (x/λ
0.45), and at the bottom of the valley (IV)
(x/λ
0.75). The unit vectors t and n indicate the tangential and
normal directions at the boundary. . . . . . . . . . . . . . . . . . . . 108

5.8

Fluid density proﬁles averaged inside regions I–IV (shown in Fig. 5.7)
near the corrugated lower wall with wavenumber ka = 0.27. The upper
wall speeds are U = 0 (a) and U = 6.0 σ/τ (b). . . . . . . . . . . . . . 109

5.9

Averaged tangential velocity proﬁles inside regions I–IV obtained from
MD simulations (open symbols) and the solution of the Navier-Stokes
equation (ﬁlled symbols). The location of the averaging regions near
the lower corrugated wall is shown in Fig. 5.7. The upper wall speed is
U = 1.0 σ/τ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

5.10 Averaged normal velocity proﬁles inside regions I–IV (shown in Fig. 5.7)
for the upper wall speed U = 1.0 σ/τ . The proﬁles are extracted from
MD simulations (open symbols) and the solution of the Navier-Stokes
equation (ﬁlled symbols). The horizontal dashed lines indicate the
location of the minima in density proﬁles (see text for details). . . . . 111
5.11 Tangential velocity proﬁles averaged inside regions I–IV (shown in
Fig. 5.7) near the lower corrugated wall. The upper wall speed is
U = 6.0 σ/τ . The velocity proﬁles are extracted from MD simulations
(open symbols) and the solution of the Navier-Stokes equation (ﬁlled
symbols). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.12 The intrinsic slip length along the lower corrugated wall (ka = 0.27)
for the upper wall speed U = 6.0 σ/τ . The data are obtained from
MD simulations (open symbols) and the solution of the Navier-Stokes
equation (ﬁlled symbols). The vertical dashed lines denote the position
of the crest (x/λ = 0.25) and the bottom of the valley (x/λ = 0.75). . 113
5.13 Averaged normal velocity proﬁles inside regions I–IV obtained from
MD simulations (open symbols) and the solution of the Navier-Stokes
equation (ﬁlled symbols). The upper wall speed is U = 6.0 σ/τ . The
horizontal arrows indicate the location of the minima in density proﬁles
between the ﬁrst and second ﬂuid layers in regions I and III. . . . . . 114
5.14 The intrinsic slip length L0 /σ as a function of the bulk pressure P
(in units of ε/σ 3 ) in the cell with ﬂat upper and lower walls. The
temperature of the Langevin thermostat is varied in the range 1.1
T kB /ε
1.35 for the upper wall speed U = 3.5 σ/τ (a) and 1.1
T kB /ε 1.4 for U = 1.0 σ/τ (b). . . . . . . . . . . . . . . . . . . . . 115
6.1

The schematic representation of the contact line and the shape of the
ﬂuid/ﬂuid interface for the static (left) and dynamic (right) cases. . . 130
131
xii

6.2

The schematic picture of the moving contact line problem where the
slip velocity and the friction coeﬃcient far from the contact line are
denoted by uslip and βslip and in the contact region by ucontact and
βcontact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
131

6.3

Initial arrangement of ﬂuid molecules and wall atoms in the system
that consists of two immiscible ﬂuids conﬁned between ﬂat walls. . . . 131
132

6.4

The modiﬁed Lennard-Jones potential Eq. (6.1) as a function of monomer
monomer distance for diﬀerent values of the parameter δ. . . . . . . . 132
133

6.5

A snapshot of an equilibrated system with two immiscible ﬂuids when
both walls are at rest. . . . . . . . . . . . . . . . . . . . . . . . . . . 132
133

6.6

The averaged density proﬁles for ﬂuid 1 across the channel width in
the z direction (a) and density proﬁles for ﬂuids 1 and 2 parallel to the
ˆ
walls in the x direction (b). . . . . . . . . . . . . . . . . . . . . . . . 134
ˆ
133

6.7

The normal stress distributions in the x, y , and z directions in equilibˆ ˆ
ˆ
rium case (U = 0) calculated from the Irving-Kirkwood relation. The
normal stresses are measured in units of ε/σ 3 . . . . . . . . . . . . . . 134
135

6.8

(a) The snapshot of the ﬂuid molecules near the solid wall within 5σ
from the contact line at equilibrium. The horizontal arrows near the
contact line indicate the direction of the tangential forces on each ﬂuid.
(b) The shear stress distribution τxz (x) along the wall as a function of
the distance from the contact line. . . . . . . . . . . . . . . . . . . . . 135
136

6.9

The shape of the ﬂuid/ﬂuid interfaces as a function of the wall speeds
for (a) U = 0, (b) U = 0.05 σ/τ , and (c) U = 0.1 σ/τ . . . . . . . . . . 136
137

6.10 The contact angle averaged in two ﬂuid phases is plotted as a function
of the wall speeds U . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137
138
6.11 The snapshot of the system near the contact line for U = 0.1 σ/τ (top),
the averaged velocity vectors (bottom) obtained from MD simulations. 138
139
6.12 The slip velocity of the ﬁrst ﬂuid layer for upper and lower wall speeds
U = 0.1 σ/τ . The velocity proﬁles are averaged separately in each
ﬂuid phase as shown by symbols. The dashed line is the slip velocity
computed regardless of the ﬂuid type. The black solid is a guide for
the eye. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139
140
6.13 (a) The slip velocity in the ﬁrst ﬂuid layer as a function of upper/lower
wall speeds. (b) The velocity proﬁles across the channel far from the
contact lines in the z direction as a function of the upper/lower wall
ˆ
speeds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140
141
xiii

6.14 The setup of the continuum simulations where the black solid lines
show the ﬁxed Cartesian grids, the red circles are the initial position
of the marker points, and the blue circles are the initial location of the
contact points. The boundary conditions are applied along the green
regions away from the contact point and along blue regions close to
the contact point (top). The schematic velocity proﬁle in the bottom
ﬁgure shows an enlarged view of the velocity proﬁle inside top ﬁgure’s
dashed-boxes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 142
141
6.15 The results from the Navier-Stokes solution of the moving contact line
as a function of upper/lower wall speeds U = 0.01 (a), U = 0.05 (b),
and U = 0.1 (c). The velocity vectors are shown by blue arrows,
streamlines by black lines, and ﬂuid/ﬂuid interfaces by red lines. . . . 142
143
6.16 The ﬁlled symbols represent the contact angle extracted from the NavierStokes solution as a function of upper/lower wall speeds. The open
symbols are the results from the molecular dynamics simulations. Solid
and dashed lines are a guide for the eye. . . . . . . . . . . . . . . . . 143
144
6.17 The symbols and solid lines show the slip velocities near the moving
contact line obtained from the Navier-Stokes solution. The dashed
lines are the slip velocities from molecular dynamics simulations in the
ﬁrst ﬂuid layer as a function of upper/lower wall speeds. . . . . . . . 144
145

xiv

CHAPTER 1
Introduction
Microﬂuidics is a multidisciplinary research ﬁeld that bridges between engineering, chemistry, and biology [1]. Microﬂuidic systems are characterized by high surface to volume ratios
and low Reynolds number ﬂows. These features allow precise control and manipulation of
ﬂuids, reduced sample volume usage, and fast mixing [2]. As new areas of application are
currently being explored and the systems under investigation are becoming more complex,
there is an increasing demand to understand the physics of ﬂows on the submicron scale
and to develop more eﬃcient numerical modeling techniques.
In the interfacial region of about a few molecular diameters from the solid wall, the
ﬂuid ﬂow has very diﬀerent physical properties than in the bulk region. One of the classical
assumptions in the ﬂuid mechanics textbooks is the no-slip boundary condition, which proved
to be very successful in describing ﬂows on large scales. However, recent experimental
and molecular dynamics simulation studies have shown that this condition needs revision
for microscopic scales due to the possibility of slip of the ﬂuid relative to the wall. The
breakdown of the no-slip condition, the so-called slip boundary condition, refers to a situation
where the ﬂuid tangential velocity adjacent to the wall is diﬀerent from the wall speed. The
measure of slip is the slip length, which is deﬁned as a distance from the real liquid-solid
interface to an imaginary plane where an extrapolated velocity proﬁle vanishes. There are
several key factors that inﬂuence the degree of slippage including surface roughness, surface
energy, ﬂuid structure, relative size of ﬂuid molecules with respect to the wall atoms, the
ratio of the ﬂuid and wall densities, and the rate of shear.

1

In this dissertation, we investigate ﬂuid ﬂows and the slip phenomena in conﬁned systems
using molecular dynamics (MD) and continuum simulations. The molecular dynamics technique, despite its computational cost, provides a detailed information about the structure,
dynamics, and thermodynamics of complex ﬂows, especially near the solid/liquid interface.
On the other hand, the continuum modeling is computationally more aﬀordable but requires
speciﬁcation of the proper boundary conditions. One of the main themes in this dissertation
is a comparative analysis of ﬂuid ﬂows using both atomistic and continuum descriptions. In
particular, we numerically solve the Navier-Stokes equation for ﬂows in various geometries
using the slip boundary conditions obtained from MD studies. The continuum predictions
are compared with the results of MD simulations for a wide range of parameters describing
liquid-solid interfaces, i.e., ﬂuid structure (monatomic and polymeric liquids), wall structure (wall density and topological roughness), and ﬂow conditions (low/high shear rates,
low/high Reynolds numbers, and single/two phase ﬂows). The results are categorized as
follows.
In chapter II, molecular dynamics simulations are carried out to investigate the dynamic
behavior of the slip length for polymeric liquids conﬁned between either atomically ﬂat or
periodically corrugated surfaces. The MD results show that for atomically ﬂat walls the
slip length depends non-linearly on shear rate. In contrast to the description of the ﬂow
over smooth boundaries (with microscopic surface roughness) by the intrinsic slip length, it
is more appropriate to characterize the ﬂow away from macroscopically rough surfaces by
the eﬀective slip length, which is usually deﬁned with respect to the location of the mean
roughness height. For periodically corrugated surfaces, the eﬀective slip length gradually
decreases as a function of wavenumber. The polymer chain conﬁguration and dynamics are
studied near rough surfaces. For ﬂows over periodically corrugated surfaces, the solution
of the Stokes equation is compared with the results of MD simulations in a wide range of
wavelengths and amplitudes.
The eﬀects of shear rate and ﬂuid density on slip boundary conditions for polymer melt
ﬂows past ﬂat crystalline walls are investigated by molecular dynamics simulations in chapter
III. It is shown that the rate-dependent slip length exhibits a local minimum and then it
increases rapidly at higher shear rates. The friction coeﬃcient at the liquid/solid interface

2

follows a power law decay as a function of the slip velocity. The ratio of the viscosity to the
slip length (i.e., friction coeﬃcient at the liquid-solid interface) is determined by the surface
induced structure in the ﬁrst ﬂuid layer.
In chapter IV, we investigate the eﬀects of local slip boundary conditions and the Reynolds
number on the ﬂow structure near periodically corrugated surfaces and the eﬀective slip
length. It was shown that for the Stokes ﬂow with the local no-slip boundary condition, the
eﬀective slip length decreases with increasing corrugation amplitude and a ﬂow circulation
is developed in suﬃciently deep grooves. The analysis of numerical solution of the NavierStokes equation with the local slip condition shows that the inertial eﬀects promote the
asymmetric vortex ﬂow formation and reduce the eﬀective slip length.
In chapter V, the slip ﬂow of monatomic ﬂuids over a periodically corrugated surface is
studied in a wide range of shear rates using molecular dynamics and continuum simulations. The eﬀective slip length in both methods is nearly constant at low shear rates and
it gradually increases at higher shear rates. The slight discrepancy between the eﬀective
slip lengths computed from MD and continuum methods at high shear rates is explained
by careful examination of the local pressure, velocity, and density proﬁles along the wavy
surface.
The moving contact line problem is considered in the last chapter VI using molecular dynamics and continuum simulations. The MD simulations are carried out to resolve the ﬂow
ﬁelds in the vicinity of the moving contact line at the molecular scale. Then, the boundary conditions extracted from MD simulations is implemented in the continuum solution of
the Navier-Stokes equations (using ﬁnite diﬀerence/particle tracking methods) to reproduce
velocity proﬁles and the shape of the ﬂuid/ﬂuid interface.

3

CHAPTER 2
Rheological study of polymer ﬂow
past rough surfaces with slip
boundary conditions
2.1

Introduction

The dynamics of ﬂuid ﬂow in conﬁned geometries has gained renewed interest due to the
recent developments in micro- and nanoﬂuidics [3]. The investigations are motivated by
important industrial applications including lubrication, coating, and painting processes.
The ﬂow behavior at the sub-micron scale strongly depends on the boundary conditions at
the liquid/solid interface. A number of experimental studies on ﬂuid ﬂow past nonwetting
surfaces have shown that the conditions at the boundary deviate from the no-slip assumption [4]. The most popular Navier model relates the slip velocity (the relative velocity of
the ﬂuid with respect to the adjacent solid wall) and the shear rate with the proportionality
coeﬃcient, the slip length, which is determined by the linear extrapolation of the ﬂuid velocity proﬁle to zero. The magnitude of the slip length depends on several key parameters,
such as wettability [5, 6, 8, 7], surface roughness [9, 10, 11, 12, 13, 14], complex ﬂuid structure [15, 16], and shear rate [17, 18, 19]. However, the experimental determination of the
slip length as a function of these parameters is hampered by the presence of several factors
with competing eﬀects on the wall slip, e.g. surface roughness and wettability [9] or surface

4

roughness and shear rate [10].
In recent years, molecular dynamics (MD) simulations have been widely used to examine
the slip ﬂow past atomically smooth, homogeneous surfaces [20, 21, 125, 23, 24, 25, 26, 28,
27, 29]. The advantage of the MD approach is that the velocity proﬁles and shear stresses
are resolved at the molecular level. The slip length in the shear ﬂow of simple ﬂuids past
crystalline walls is a function of the wall-ﬂuid density ratio [21, 23], the relative size of wall
atoms and ﬂuid molecules [23, 28], the surface energy [21, 23, 29], and the interfacial shear
rate [23, 27, 29]. Weak wall-ﬂuid interactions and incommensurable structures of the solid
and ﬂuid phases at the interface usually lead to enhancement of slip [21, 23, 26, 29]. If the slip
length at low shear rates is about several molecular diameters then it increases with the shear
rate, and the slope of the rate dependence is greater for weaker wall-ﬂuid interactions [23, 29].
The rate dependence of the slip length in the ﬂow of polymer melts is more complicated
because of the additional length and time scales associated with the dynamics of polymer
chains at the interface and the shear thinning viscosity [125, 30, 31, 27, 32].
In the presence of surface roughness [33, 34, 35, 36, 37] or chemical patterning [38, 39,
40, 41, 42], the ﬂuid ﬂow near the solid boundary is perturbed on the length scales of the
surface heterogeneities and its description requires deﬁnitions of the eﬀective slip length
and the average location of the reference plane. Most commonly, the location of the reference plane is deﬁned as the mean height of the surface asperities, and the shear rate is
determined by averaging of the ﬂuid ﬂow over the typical length scale of the surface inhomogeneities. In general, the surface roughness is expected to reduce the eﬀective slip
length for wetting liquids [34, 35, 37]. For suﬃciently rough surfaces the no-slip boundary
condition can be achieved even if the local condition is of zero shear stress [43]. However, in special cases, when the ﬂuid is partially dewetted at the nanostructured, or the
so-called “superhydrophobic” surfaces, the slip length might be enhanced up to a few microns [44, 45, 46, 47, 48, 49, 50].
In recent MD studies on shear ﬂow of simple ﬂuids, the behavior of the eﬀective slip length
was investigated in the Couette cell with either mixed boundary conditions [40] or periodic
surface roughness [35]. In the ﬁrst study, the lower stationary wall with mixed boundary
conditions was patterned with a periodic array of stripes representing alternating regions of

5

ﬁnite slip and zero shear stress. In the other study [35], the periodically roughened surface
was modeled by introducing a sinusoidal oﬀset to the position of the wall atoms. At the wavy
wall, the local slip length is modiﬁed by the presence of curvature and becomes positiondependent along the curved boundary [51, 52]. A detailed comparison between continuum
analysis and MD simulations shows an excellent agreement between the velocity proﬁles
and eﬀective slip lengths when the characteristic length scale of substrate inhomogeneities is
larger than approximately thirty molecular diameters [40, 35]. In the case of rough surfaces,
an additional correction due to variable wall density was incorporated in the analysis [35].
The problem of applicability of the results obtained for monoatomic ﬂuids to polymer melts
is important for modeling polymer ﬂows in conﬁned geometries and design of the hybrid
continuum-atomistic algorithms.
In this chapter, the MD simulations are carried out to study the dynamic behavior of
the slip length at the interface between a polymer melt and atomically ﬂat or periodically
corrugated surfaces. The MD results for ﬂat crystalline walls conﬁrm previous ﬁndings [32]
that the slip length goes through a local minimum at low shear rates and then increases
rapidly at higher shear rates. For periodically corrugated surfaces and wetting conditions,
the eﬀective slip length decreases gradually with increasing values of the wavenumber. The
solution of the Stokes equation with either constant or rate-dependent local slip length
is compared with the MD simulations for corrugated surfaces with wavelengths ranging
from molecular dimensions to values much larger than the radius of gyration of polymer
chains. The orientation and the dynamics of linear polymer chains are signiﬁcantly aﬀected
by surface roughness when the corrugation wavelengths are comparable with the radius of
gyration.
The rest of the chapter is organized as follows. The details of molecular dynamics and
continuum simulations are described in the next section. The dynamic response of the slip
length and shear viscosity in the cell with atomically ﬂat surfaces is presented in Sec. 2.3.
The results of MD simulations for corrugated walls with large wavelengths and comparison
with continuum predictions are reported in Sec. 2.4. The slip behavior for small wavelengths
and the conformational properties of the polymer chains near the rough surfaces are analyzed
in Sec. 2.5. The summary is given in the last section.

6

2.2

The details of the numerical simulations

2.2.1

Molecular dynamics model

The computational domain consists of a polymeric ﬂuid conﬁned between two atomistic
walls. Figure 2.1 shows the MD simulation setup. Any two ﬂuid monomers within a cut-oﬀ
distance of rc = 2.5 σ interact through the truncated Lennard-Jones (LJ) potential
VLJ (r) = 4ε

σ 6
σ 12
,
−
r
r

(2.1)

where ε is the energy scale and σ is the length scale of the ﬂuid phase. The LJ potential
was also employed for the wall-ﬂuid interaction with εwf = ε and σwf = σ. The wall
atoms do not interact with each other, and the wall-ﬂuid parameters are ﬁxed throughout
the study. In addition to the LJ potential, the neighboring monomers in a polymer chain
(N = 20 beads) interact with the ﬁnite extensible nonlinear elastic (FENE) potential
1 2
r 2
VF EN E (r) = − kro ln 1 −
,
2
ro

(2.2)

where ro = 1.5 σ and k = 30 εσ −2 [53]. The MD simulations were performed at a constant
density ensemble with ρ = 0.88 σ −3 . The total number of ﬂuid monomers is Nf = 67200.
The motion of the ﬂuid monomers was weakly coupled to an external thermal reservoir [54]. To avoid a bias in the ﬂow direction, the random force and the friction term were
added to the equation of motion in the y direction [21]
ˆ
myi = −
¨
i=j

∂(VLJ + VF EN E )
− mΓyi + fi (t),
˙
∂yi

(2.3)

where Γ = 1.0 τ −1 is a friction constant that regulates the rate of heat ﬂux between the ﬂuid
and the heat bath, and fi (t) is the random force with zero mean and variance 2mΓkB T δ(t)
determined from the ﬂuctuation-dissipation theorem [55]. Temperature of the Langevin
thermostat is set to T = 1.1 ε/kB , where kB is the Boltzmann constant. The equations of
motion were integrated using the Verlet algorithm [56] with a time step ∆t = 0.005 τ , where
τ=

mσ 2 /ε is the characteristic LJ time.

The dimensions of the system in the xy plane, unless speciﬁed otherwise, were set to
Lx = 66.60 σ and Ly = 15.59 σ. The upper wall was composed of two layers of an F CC
7

lattice with density ρw = 1.94 σ −3 , which corresponds to the nearest-neighbor distance of
d = 0.9 σ between wall atoms in the (111) plane. The lower wall was constructed of two F CC
layers of atoms distributed along the sinusoidal curve with the wavelength λ and amplitude
a. For the largest wavelength λ = 66.60 σ, the density of the lower wall ρw = 1.94 σ −3 was
kept uniform along the sinusoid (by including additional rows of atoms parallel to the y axis)
ˆ
to avoid additional analysis of the eﬀective slip length due to variable wall density [35]. In
the present study, the corrugation amplitude was varied in the range 0

a/σ

12.04. In

the absence of the imposed corrugation (a = 0) the distance between the inner F CC planes
is set to Lz = 74.15 σ in the z direction. Periodic boundary conditions were imposed in the
ˆ
x and y directions.
ˆ
ˆ
The initial velocities of the ﬂuid monomers were chosen from the Maxwell-Boltzmann
probability distribution at the temperature T = 1.1 ε/kB . After an equilibration period
of about 3 × 104 τ with stationary walls, the velocity of the upper wall was gradually
increased in the x direction from zero to its ﬁnal value during the next 2 × 103 τ . Then
ˆ
the system was equilibrated for an additional period of 6 × 103 τ to reach steady-state.
Averaging time varied from 105 τ to 2 × 105 τ for large and small velocities of the upper wall
respectively. The velocity proﬁles were averaged within horizontal slices of Lx × Ly × ∆z,
where ∆z = 0.2 σ. Fluid density proﬁles near the walls were computed within slices with
thickness ∆z = 0.01 σ [29].

2.2.2

Continuum method

A solver based on the ﬁnite element method was developed for the two-dimensional steadystate and incompressible Navier-Stokes (NS) equation. The NS equation with these assumptions is reduced to
ρ (u ·

u) = − p + µ 2 u,

(2.4)

where u is the velocity vector, ρ is the ﬂuid density, and p and µ are the pressure ﬁeld and
viscosity of the ﬂuid respectively.
The incompressibility condition is satisﬁed by a divergence-free velocity ﬁeld u. In order
to avoid the decoupling of the velocity and the pressure ﬁelds in the numerical simulation

8

of the incompressible ﬂow, the penalty formulation is adopted [57]. This method replaces
the continuity equation,

· u = 0, with a perturbed equation
p
·u=− ,
Λ

(2.5)

where Λ is the penalty parameter. For most practical applications, where computation is
performed with double-precision 64 bit words, a penalty parameter between 107 and 109 is
suﬃcient to conserve the accuracy [57]. In our simulations the incompressibility constraint
was set to Λ = 107 . The pressure term in Eq. (2.5) is then substituted into the NS equation.
The equation Eq. (2.4) can be rewritten as follows:
ρ (u ·

u) = Λ (

· u) + µ 2 u,

(2.6)

where the continuity equation is no longer necessary [57]. The NS equation is integrated
with the Galerkin method using bilinear rectangular isoparametric elements [57].
Four boundary conditions must be speciﬁed for the continuum simulation. Periodic
boundary conditions are used for the inlet and outlet along the x direction. A slip boundary
ˆ
condition is applied at the upper and lower walls. In the local coordinate system (spanned
by the tangential vector t and the normal vector n), the ﬂuid velocity along the curved
boundary is calculated as
ut = L0 [(n ·

)ut + ut /R],

(2.7)

where ut is the tangential component of u = ut t + un n, L0 is the slip length at the ﬂat
liquid/solid interface, the term in the brackets is the local shear rate, and R is the local
radius of curvature [52]. The radius of curvature is positive for the concave and negative for
the convex regions. For a ﬂat surface, R → ∞, the boundary condition given by Eq. (2.7)
simply becomes the Navier slip law.
The simulation is started by applying the no-slip boundary condition as the initial guess.
Once the equations of motion are solved implicitly, the local slip velocities at the lower and
upper boundaries are updated using Eq. (2.7). In the next step, the equations of motion
are solved with the updated slip velocities used as a new boundary condition. The iterative
procedure is repeated until the solution converges to a desired accuracy. The convergence
rate of the iterative solution remains under control with the under-relaxation value about
9

0.001 for the boundary nodes. In all continuum simulations, the grids at the lower boundary
have an aspect ratio of about one. The computational cost is reduced by increasing the
aspect ratio of the grids in the bulk region.
The normalized average error value in the simulation is deﬁned as
Np

error =
i=1

| un − un+1 |
i
i
| un+1 |
i

/Np ,

(2.8)

where Np is the number of nodes in the system, un is the velocity at the node i and time
i
step n, and un+1 is the velocity in the next time step. The typical error in the converged
i
∂u

solution is less than 10−9 . At the boundaries the solution satisﬁes ut = Llocal ∂nt , where
the local slip length is Llocal = (L−1 − R−1 )−1 [52].
0

2.3

MD results for ﬂat walls

The averaged velocity proﬁles for selected values of the upper wall speed U are presented
in Fig. 2.2. The proﬁles are linear throughout the cell, except for U

6.0 σ/τ where a

slight curvature appears in the region of about 4 σ near the walls. Note that the relative
slip velocity at the upper and lower walls increases with increasing upper wall speed. The
shear rate was determined from the linear ﬁt to the velocity proﬁles across the whole width
of the channel (see Fig. 2.2). The uncertainty in the estimated value of shear rate is due
to the thermal ﬂuctuations and the slight curvature in the velocity proﬁles near the walls.
The typical error bars for the shear rate are about 2 × 10−5 τ −1 and 6 × 10−4 τ −1 for small
and large upper wall speeds respectively (not shown).
In this study, the shear stress in steady-state ﬂow was computed from the Irving-Kirkwood
relation [58]. The dynamic response of the ﬂuid viscosity with increasing shear rate is
presented in Fig. 2.3. At higher shear rates, the ﬂuid exhibits shear thinning behavior with
the slope of about −0.33. Although the power law coeﬃcient is larger than the reported
values in experimental studies [59], the results are consistent with previous MD simulations
of polymer melts for similar ﬂow conditions [60, 61, 32].
The slip length was calculated by the linear extrapolation of the ﬂuid velocity proﬁle to
zero with respect to a reference plane, which is located 0.5 σ away from the F CC lattice
10

plane [29, 32]. Figure 2.4 shows the rate dependence of the slip length in the same range
of shear rates as in Fig. 2.3. The slip length goes through a shallow minimum at low shear
rates and then increases rapidly at higher rates. The error bars are larger at low shear rates
2
because the thermal ﬂuid velocity vT = kB T/m is greater than the average ﬂow velocity.

The nonmonotonic behavior of the slip length in sheared polymer ﬁlms with atomically ﬂat
surfaces can be interpreted in terms of the friction coeﬃcient at the liquid/solid interface,
which undergoes a gradual transition from a nearly constant value to the power law decay
as a function of the slip velocity [32]. The data for the slip length shown in Fig. 2.4 are well
ﬁtted by the fourth-order polynomial
L0 (x)/σ = 16.8 − 72.0 × 10 x + 44.0 × 103 x2 − 97.3 × 104 x3 + 80.5 × 105 x4 ,

(2.9)

where x = γτ is the shear rate. The polynomial ﬁt will be used to specify the boundary
˙
conditions for the continuum solution described in the next section.

2.4

Results for periodically corrugated walls: large
wavelength

2.4.1

MD simulations

The monomer density proﬁles were computed in the averaging regions located in the grooves
and above the peaks of the corrugated lower wall with wavelength λ = Lz = 66.60 σ. The
dimensions of the averaging regions are set to 0.5 σ and Ly = 15.59 σ in the x and y directions
ˆ
ˆ
respectively. Figure 2.5 shows the monomer density proﬁles near the upper and lower walls at
equilibrium. The pronounced density oscillations are attributed to the successive layering of
the ﬂuid monomers near the walls. These oscillations decay within a few molecular diameters
from the walls to a uniform proﬁle characterized by the bulk density of ρ = 0.88 σ −3 . Note
that the height of the ﬁrst peak in the density proﬁle inside the groove is slightly larger
than its value above the crest [see Fig. 2.5(b)]. The ﬂuid monomers experience stronger net
surface potential in the groove than above the crest because of the closer spatial arrangement
of the wall atoms around the location of the ﬁrst density peak in the groove. This eﬀect

11

is ampliﬁed when the local radius of curvature at the bottom of the grooves is reduced at
smaller wavelengths (see below).
The averaged velocity proﬁles in the cell with periodically corrugated lower and ﬂat upper
walls are presented in Fig. 2.6. The ﬂuid velocity proﬁles were averaged within horizontal
slices of thickness ∆z = 0.2 σ distributed uniformly from the bottom of the grooves to the
upper wall. With increasing corrugation amplitude, the slip velocity decreases and the
proﬁles acquire a curvature near the lower wall. Note also that the relative slip velocity
at the upper wall increases with the shear rate. The linear part of the velocity proﬁles,
30

z/σ

60, was used to determine the eﬀective slip length Leﬀ at the corrugated lower

wall. The variation of the eﬀective slip length as a function of wavenumber ka = 2 πa/λ is
shown in Fig. 2.7. The slip length decreases monotonically with increasing wavenumber and
becomes negative at ka

0.7. The results of comparison between the MD and continuum

simulations for three diﬀerent cases are presented in the next section.

2.4.2

Comparison between MD and continuum simulations

In continuum simulations, the length, time and energy scales are normalized by the LJ
parameters σ, τ and ε respectively. The continuum nondimensional parameters are denoted
by the (˜) sign. The size of the two-dimensional domain is ﬁxed to 66.60×73.15 in the x and
ˆ
z directions, respectively. The following three cases were examined: the Stokes solution with
ˆ
constant slip length in Eq. (2.7), the Stokes solution with shear-rate-dependent slip length
given by Eq. (2.9), and the Navier-Stokes solution with constant slip length in Eq. (2.7). For
˜
all cases considered, the ﬂat upper wall is translated with a constant velocity U = 0.5 in the
x direction. Similarly to the MD method, the eﬀective slip length is deﬁned as a distance
ˆ
from the reference plane at a = 0 to the point where the linearly extrapolated velocity
proﬁle vanishes.
In the ﬁrst case, the ﬁnite element method was implemented to solve the Stokes equation
˜
with boundary conditions at the upper and lower walls speciﬁed by Eq. (2.7) with L0 = 14.1.
This value corresponds to the slip length L0 = 14.1±0.5 σ extracted from the MD simulations
in the cell with ﬂat walls and the velocity of the upper wall U = 0.5 σ/τ . Note that at
large amplitudes a
˜

8.16, the normal derivative of the tangential velocity ∂ut /∂n at the
12

bottom of the grooves is negative, while the x component of the slip velocity is positive
ˆ
everywhere along the corrugated lower wall. The dependence of the eﬀective slip length on
the corrugation amplitude is shown in Fig. 2.7. The continuum results agree well with the
approximate analytical solution [52] for ka

0.5 (not shown). For larger amplitudes, ka >

0.5, where the analytical solution is not valid, our results were tested to be grid independent.
There is an excellent agreement between slip lengths obtained from the MD and continuum
simulations for ka

0.3. With further increasing the amplitude, the slip length obtained

from the continuum solution overestimates its MD value. The results presented in Fig. 2.7
are consistent with the analysis performed earlier for simple ﬂuids [35], although a better
agreement between MD and continuum solutions was expected at ka > 0.5 because of the
larger system size considered in the present study.
As discussed in the previous section, the slip length for atomically ﬂat walls is ratedependent even at low shear rates (see Fig. 2.4). In the second case, we include the eﬀect
of shear rate in the analysis of the eﬀective slip length at the corrugated lower wall and
ﬂat upper wall. The Stokes equation is solved with boundary conditions given by Eq. (2.7),
where the slip length Eq. (2.9) is a function of the local shear rate at the curved and ﬂat
boundaries. The results obtained from the Stokes solution with constant and rate-dependent
slip lengths are almost indistinguishable (see Fig. 2.7). This behavior can be attributed to
a small variation of the intrinsic slip length Eq. (2.9) at low shear rates. For example, at
the largest amplitude, a = 12.04, the local shear rate at the corrugated wall is position˜
dependent and bounded by |∂ut /∂n + ut /R|

0.0035. In this range of shear rates, the

normalized value of the slip length in Eq. (2.9) varies between 14.7

˜
L0

16.6. It is

expected, however, that the eﬀect of shear rate will be noticeable at larger values of the top
˜
wall speed U .
˜
In the third case, the Navier-Stokes equation is solved with a constant slip length L0 = 14.1
in Eq. (2.7) at the ﬂat upper and corrugated lower walls. The upper estimate of the Reynolds
number based on the ﬂuid density ρ = 0.88, viscosity µ = 20.0, and the ﬂuid velocity diﬀer˜
˜
ence across the channel is Re ≈ 1.3. It was previously shown by Tuck and Kouzoubov [63]
˜
that at small ka and L0 = 0 the magnitude of the apparent slip velocity at the mean surface
increases due to ﬁnite Reynolds number eﬀects for Re

13

30. In our study, the diﬀerence

between the slip lengths extracted from the Stokes and Navier-Stokes solutions is within
the error bars (see Fig. 2.7). These results conﬁrm that the slip length is not aﬀected by the
inertia term in the Navier-Stokes equation for Re

1.3. To check how sensitive the bound-

ary conditions are to higher Reynolds number ﬂows, we have also repeated the continuum
˜
simulations for larger velocity of the upper wall U = 50, which corresponds to Re ≈ 130. For
the largest corrugation amplitude a = 12.04, the backﬂow appears inside the groove and the
˜
˜
eﬀective slip length becomes smaller than its value for U = 0.5 by about 0.7 (not shown).

2.5

Results for periodically corrugated walls: small
wavelengths

2.5.1

Comparison between MD and continuum simulations

The MD simulations described in this section were performed at corrugation wavelengths
(λ/σ = 3.75, 7.5 and 22.5) comparable with the size of a polymer coil. Periodic surface roughness of the lower wall was created by displacing the F CC wall atoms by ∆z = a sin(2 πx/λ)
in the z direction [35]. In order to reduce the computational time the system size was
ˆ
restricted to Nf = 8580 ﬂuid monomers and Lx = 22.5 σ, Ly = 12.5 σ and Lz = 35.6 σ. All
other system parameters were kept the same as in the previous section.
The representative density proﬁles near the upper and lower walls are shown in Fig. 2.8 for
the wavelength λ = 7.5 σ. The height of the ﬁrst peak in the density proﬁle is larger in the
grooves than near the ﬂat wall or above the crests of the corrugated surface. The eﬀective slip
length as a function of wavenumber ka is plotted in Fig. 2.9. For all wavelengths, the slip
length decreases monotonically with increasing values of ka. At the smallest wavelength
λ = 3.75 σ, the slip length rapidly decays to zero at ka ≈ 0.4 and weakly depends on the
corrugation amplitude at larger ka. Inspection of the local velocity proﬁles for λ = 3.75 σ
and ka

1.0 indicates that the ﬂow is stagnant inside the grooves [62].

˜
In the continuum analysis, the Stokes equation with a constant slip length L0 = 14.1 in
Eq. (2.7) is solved for the three wavelengths. The comparison between the MD results and
the solution of the Stokes equation is presented in Fig. 2.9. The error bars are larger for the

14

smallest wavelength because of the ﬁne grid resolution required near the lower boundary at
ka

0.2 (see inset in Fig. 2.9). The results shown in Fig. 2.9 conﬁrm previous ﬁndings for

simple ﬂuids [35] that the slip length obtained from the Stokes ﬂow solution overestimates
its MD value and the agreement between the two solutions becomes worse at smaller wavelengths. It is interesting to note that the curves for diﬀerent wavelengths intersect each
other at ka ≈ 0.63 in the MD model and at ka ≈ 1.02 in the continuum analysis. The same
trend was also observed in the previous study on slip ﬂow of simple ﬂuids past periodically
corrugated surfaces [35].

2.5.2

The polymer chain conﬁguration and dynamics near rough
surfaces

In this section, the properties of polymer chains are examined in the bulk and near the
corrugated boundary with wavelengths λ = 3.5 σ and λ = 7.5 σ. The radius of gyration Rg
was computed as
1
2
Rg =

N

N

(Ri − Rcm )2 ,

(2.10)

i=1

where Ri is the three-dimensional position vector of a monomer, N = 20 is the number of
monomers in the chain, and Rcm is center of mass vector deﬁned as
1
Rcm =
N

N

Ri .

(2.11)

i=1

The chain statistics were collected in four diﬀerent regions at equilibrium (U = 0) and
in the shear ﬂow induced by the upper wall moving with velocity U = 0.5 σ/τ in the x
ˆ
direction. Averaging regions were located above the peaks, in the grooves, near the ﬂat
upper wall and in the bulk (see Fig. 2.10 for an example). The dimensions of the averaging
regions above the peaks and in the grooves of the lower wall are σ × 12.5 σ × 1.5 σ, and near
the upper wall and in the bulk are 22.5 σ × 12.5 σ × 1.5 σ. Three components of the radius
of gyration were computed for polymer chains with the center of mass inside the averaging
regions.
In the bulk region, the components of the radius of gyration remain the same for both
wavelengths, indicating that the chain orientation is isotropic at equilibrium and is not
15

aﬀected by the conﬁning walls. In the steady-state ﬂow, the eﬀective slip length is suppressed by the surface roughness and, therefore, the shear rate in the bulk increases with
the corrugation amplitude. This explains why the x component of the radius of gyration
ˆ
Rgx increases slightly at larger amplitudes (see Tables 2.1 and 2.2). Near the upper wall,
the polymer chains become ﬂattened parallel to the surface and slightly stretched in the
presence of shear ﬂow. These results are consistent with the previous MD simulations of
polymer melts conﬁned between atomically ﬂat walls [65, 30, 64].
In the case of a rough surface with the wavelength λ = 7.5 σ, a polymer chain can be
accommodated inside a groove (see Fig. 2.10 for an example). With increasing corrugation
amplitude, the polymer chains inside the grooves elongate along the y direction and contract
ˆ
in the x direction (see Table 2.1). The tendency of the trapped molecules to orient parallel to
ˆ
the grooves was observed previously in MD simulations of hexadecane [34]. In the presence
of shear ﬂow, polymer chains are highly stretched in the x direction above the crests of the
ˆ
wavy wall. A snapshot of the unfolded chains during migration between neighboring valleys
is shown in Fig. 2.10. The ﬂow conditions in Fig. 2.10 correspond to a negative eﬀective slip
length Leﬀ ≈ −2 σ.
For the smallest corrugation wavelength λ = 3.75 σ, polymer chains cannot easily ﬁt in
the grooves unless highly stretched. Therefore, the y component of the radius of gyration
ˆ
is relatively large when the center of mass is located in the deep grooves (see Table 2.2).
Figure 2.11 shows a snapshot of several polymer chains in contact with the lower corrugated
wall. The chain segments are oriented parallel to the grooves and stretched above the crests
of the surface corrugation. Visual inspection of the consecutive snapshots reveals that the
chains near the corrugated wall move, on average, in the direction of shear ﬂow; however,
their tails can be trapped for a long time because of the strong net surface potential inside
the grooves. For large wavenumbers ka

0.5, the magnitude of the negative eﬀective

slip length is approximately equal to the sum of the corrugation amplitude and Rgz of the
polymer chains above the crests of the wavy wall.

16

2.6

Conclusions

In this chapter the eﬀects of the shear rate and surface roughness on slip ﬂow of a polymer melt was studied using molecular dynamics and continuum simulations. The linear
part of the velocity proﬁles in the steady-state ﬂow was used to calculate the eﬀective slip
length and shear rate. For atomically ﬂat walls, the slip length passes through a shallow
minimum at low shear rates and then increases rapidly at higher shear rates. In the case
of periodic surface heterogeneities with the wavelength larger than the radius of gyration,
the eﬀective slip length decays monotonically with increasing the corrugation amplitude.
For small wavenumbers, the eﬀective slip length obtained from the solution of the Stokes
equation with constant and shear-dependent local slip length is in a good agreement with
its values computed from the MD simulations, in accordance with the previous analysis for
simple ﬂuids [35]. At low Reynolds numbers, the inertial eﬀects on slip boundary conditions
are negligible. When the corrugation wavelengths are comparable to the radius of gyration,
polymer chains stretch in the direction of shear ﬂow above the crests of the surface corrugation, while the chains located in the grooves elongate perpendicular to the ﬂow. In this
regime the continuum approach fails to describe accurately the rapid decay of the eﬀective
slip length with increasing wavenumber.

17

U

Solid Wall

ux(z)

Leff

Lz

λ

z
y

a
x

Figure 2.1. For interpretation of the references to color in this and all other ﬁgures, the
reader is referred to the electronic version of this dissertation. A snapshot of the ﬂuid
monomers (open circles) conﬁned between solid walls (closed circles) obtained from the MD
simulations. The atoms of the stationary lower wall are distributed along the sinusoidal
curve with the wavelength λ and amplitude a. The ﬂat upper wall is moving with a constant
velocity U in the x direction. The eﬀective slip length Leﬀ is determined by the linear
ˆ
extrapolation of the velocity proﬁle to ux = 0.

18

8
U = 1.0 σ/τ
U = 2.0 σ/τ
U = 4.0 σ/τ
U = 6.0 σ/τ
U = 8.0 σ/τ

7

< ux> τ / σ

6
5
4
3
2
1
0

10

20

30

40

z /σ

50

60

70

Figure 2.2. Averaged velocity proﬁles in the cell with ﬂat upper and lower walls. The solid
lines are the best linear ﬁt to the data. The vertical axes indicate the location of the F CC
lattice planes. The velocities of the upper wall are tabulated in the inset.

19

1.4
Lx = 66.60 σ
Lz = 74.15 σ

log10 ( µ /ε τ σ −3)

1.3
1.2
1.1
1.0

-0.33

ρw = 1.94 σ -3
εwf = 1.0 ε

0.9
0.8
-3.5

σwf = 1.0 σ

-3.0

-2.5

. -2.0
log10(γτ)

-1.5

-1.0

Figure 2.3. Viscosity of the polymer melt µ/ετ σ −3 as a function of shear rate. The dashed
line with the slope −0.33 is plotted for reference.

20

30
25

L s /σ

20

ρw = 1.94 σ -3
εwf = 1.0 ε
σwf = 1.0 σ

15
10
Lx = 66.60 σ
Lz = 74.15 σ

5
0
-3.5

-3.0

-2.5

. -2.0
log10(γτ)

-1.5

-1.0

Figure 2.4. Variation of the slip length as a function of shear rate in the cell with ﬂat upper
and lower walls. The solid line is a fourth-order polynomial ﬁt to the data given by Eq. (2.9).

21

4

68

(a)

70

72

74

4

Flat upper wall

ρσ3

2

2

0

0

4

(b)

4

Inside the groove
Above the peak

2

0

2

-6

-4

z /σ

10

12

0

Figure 2.5. Averaged density proﬁles near the stationary upper wall (a), above the peak
and in the groove of the lower wall with amplitude a = 8.16 σ and wavelength λ = 66.60 σ
(b).

22

0.5
Flat walls
a = 3.52 σ
a = 7.24 σ
a = 12.04 σ

< ux> τ / σ

0.4
0.3
0.2

U = 0.5 σ /τ
λ = 66.60 σ
Lz = 74.15 σ

0.1
0.0
-10

0

10

20

30

z /σ

40

50

60

70

Figure 2.6. Averaged velocity proﬁles for the indicated values of the corrugation amplitude
a. The vertical dashed line denotes a reference plane for calculation of the eﬀective slip
length at the corrugated lower wall. The velocity of the ﬂat upper wall is U = 0.5 σ/τ .

23

16
MD simulation

12

Stokes L0= constant
.
Stokes L0( γ τ )
Navier-Stokes

Leff /σ

8
4
0
-4
0.0

U = 0.5 σ /τ
λ = 66.60 σ
Lz = 74.15 σ

0.2

0.4

0.6

0.8

1.0

1.2

ka
Figure 2.7. The eﬀective slip length as a function of wavenumber ka obtained from the MD
simulations ( ), the solution of the Stokes equation with rate-independent slip length L0
(◦) and with L0 (γτ ) given by Eq. (2.9) (×), the solution of the Navier-Stokes equation with
˙
L0 ( ).

24

28

4

30

(a)

32

34

36

4

Flat upper wall

ρσ3

2

2

0
6

0
6
Inside the groove
Above the peak

(b)
4

4

2
0

2
0

2

4

z /σ

6

0
8

Figure 2.8. Averaged ﬂuid density proﬁles near ﬂat upper wall (a) and corrugated lower
wall with amplitude a = 1.4 σ and wavelength λ = 7.5 σ (b). The velocity of the upper wall
is U = 0.5 σ/τ .

25

12

Leff /σ

12

Leff /σ

8

8
4

4

0.0

0.1

0.2

ka

0
-4
-8
0.0

λ = 3.75 σ
λ = 7.5 σ
λ = 22.5 σ

0.4

0.8

1.2

1.6

ka
Figure 2.9. The eﬀective slip length as a function of wavenumber ka for the indicated values
of wavelength λ. Continuum results are denoted by the dashed lines and open symbols,
while the MD results are shown by straight lines and ﬁlled symbols. The inset shows the
same data for ka 0.25.

26

λ = 7.5 σ
Bulk

Equilibrium
Shear ﬂow
Upper Equilibrium
wall
Shear ﬂow
Peak Equilibrium
Shear ﬂow
Groove Equilibrium
Shear ﬂow

a = 0.2 σ
Rgx
Rgy
1.18
1.18
1.70
1.11
1.34
1.37
1.63
1.29
1.35
1.36
1.65
1.28
1.31
1.37
1.50
1.34

Rgz
1.18
1.04
0.66
0.63
0.69
0.67
0.63
0.62

a = 0.6 σ
Rgx
Rgy
1.18
1.18
1.76
1.10
1.35
1.36
1.64
1.29
1.40
1.31
1.85
1.21
1.10
1.47
1.19
1.44

Rgz
1.18
1.02
0.66
0.63
0.76
0.75
0.61
0.60

a = 1.4 σ
Rgx
Rgy
1.18
1.18
1.79
1.10
1.35
1.35
1.67
1.29
1.46
1.27
2.17
1.10
0.77
1.69
0.75
1.89

Rgz
1.18
1.01
0.66
0.63
0.84
0.93
0.61
0.59

27

Table 2.1. Averaged x, y , and z components of the radius of gyration at equilibrium and in the shear ﬂow. The Rg values are reported in
ˆ ˆ
ˆ
the bulk, near the ﬂat upper wall, above the peaks, and inside the grooves. The wavelength of the lower wall is λ = 7.5 σ. The size of the
averaging region inside the grooves and above the peaks is σ × 12.5 σ × 1.5 σ. The estimate of the error bars is ±0.03 σ.

27

λ = 3.75 σ
Bulk

Equilibrium
Shear ﬂow
Upper Equilibrium
wall
Shear ﬂow
Peak Equilibrium
Shear ﬂow
Groove Equilibrium
Shear ﬂow

a = 0.07 σ
Rgx
Rgy
1.18
1.18
1.70
1.12
1.33
1.36
1.60
1.32
1.38
1.33
1.61
1.27
1.33
1.36
1.60
1.29

Rgz
1.18
1.03
0.66
0.63
0.67
0.65
0.64
0.61

a = 0.2 σ
Rgx
Rgy
1.18
1.18
1.76
1.10
1.36
1.34
1.64
1.29
1.42
1.30
1.66
1.26
1.25
1.42
1.60
1.33

Rgz
1.18
1.02
0.66
0.63
0.71
0.69
0.61
0.59

a = 1.0 σ
Rgx
Rgy
1.18
1.18
1.78
1.10
1.37
1.34
1.65
1.27
1.41
1.25
1.78
1.17
0.48
2.64
0.47
2.70

Rgz
1.18
1.01
0.65
0.63
0.92
0.84
0.55
0.54

28

Table 2.2. Averaged x, y , and z components of the radius of gyration at equilibrium and in the shear ﬂow. The Rg values are reported in
ˆ ˆ
ˆ
the bulk, near the ﬂat upper wall, above the peaks, and inside the grooves. The wavelength of the lower wall is λ = 3.75 σ. The dimensions
of the averaging region are the same as in Table 2.1.

28

12.5 σ

σ
1.5 σ

z y
x

Figure 2.10. A snapshot of four polymer chains near the lower corrugated wall for wavelength
λ = 7.5 σ and amplitude a = 1.4 σ. The velocity of the upper wall is U = 0.5 σ/τ .

29

z

σ

x

1.5 σ

y
x

σ

12.5 σ
Figure 2.11. A snapshot of ﬁve polymer chains near the lower corrugated wall for the
wavelength λ = 3.75 σ and amplitude a = 1.0 σ. The ﬁgure shows the side view (top) and
the top view (bottom). The velocity of the upper wall is U = 0.5 σ/τ .

30

CHAPTER 3
Slip boundary conditions for shear
ﬂow of polymer melts past atomically
ﬂat surfaces
3.1

Introduction

The ﬂuid ﬂow in microﬂuidic channels can be signiﬁcantly inﬂuenced by liquid slip at the
solid boundary because of the large surface to volume ratio [67]. The Navier model for
the partial slip boundary conditions relates the ﬂuid slip velocity to the tangential viscous
stress at the wall via the friction coeﬃcient, which is deﬁned as the ratio of the shear
viscosity to the slip length. Geometrically, the slip length is determined as a distance from
the solid wall where the linearly extrapolated ﬂuid velocity proﬁle vanishes. Experimental
studies have demonstrated that for atomically smooth surfaces the slip length is relatively
large for non-wetting liquid/solid interfaces [5, 15, 117], polymeric ﬂuids [68, 16] and high
shear rates [69, 17, 19]. By contrast, surface roughness even on the submicron length scale
strongly reduces the degree of slip for both polymeric [12] and Newtonian [10, 14, 13] ﬂuids.
However, experimentally it is often diﬃcult to isolate and, consequently, to analyze the
eﬀects of wettability, surface roughness and shear rate on slip boundary conditions because
of the small scales involved.
During the past two decades, a number of molecular dynamics (MD) simulations [20,

31

21, 121, 25, 24, 26, 27, 28, 40, 29] have been performed to investigate the dependence of
the slip length on the structure of dense ﬂuids near atomically ﬂat surfaces. For simple
ﬂuids, it has been well established that the slip length correlates inversely with the wallﬂuid interaction energy and the surface induced order in the ﬁrst ﬂuid layer [21, 121, 29].
The strong wall-ﬂuid coupling and commensurate structures of the liquid/solid interface
lead to the no-slip boundary condition for monoatomic ﬂuids [21, 28]. The interfacial slip
is enhanced in the ﬂow of polymer melts due to the higher shear viscosity and reduced
molecular ordering near the ﬂat wall [125, 30, 31, 27]. The magnitude of the slip length in
the ﬂow of polycarbonate blends past atomically smooth nickel surfaces is determined by the
competition between the apparent slip of the adsorbed polymer layer and its thickness [70].
The slip boundary conditions for the polymer melt ﬂow over the brush can be controlled
by changing the density of the end-grafted polymer layer [71]. Finally, the degree of slip at
the interface between immiscible polymers is larger for longer chain molecules and greater
repulsion between unlike species [72].
Until ten years ago, the slip length in the Navier model was assumed to be independent of
shear rate. In the seminal MD study by Thompson and Troian [23] on shear ﬂow of simple
ﬂuids past atomically ﬂat rigid walls, it was found that the slip length increases nonlinearly
with the shear rate. Upon increasing the surface energy, a gradual transition from the power
law behavior [23] to the linear rate dependence was recently reported for weak wall-ﬂuid
interactions and incommensurate structures of the liquid/solid interface [29]. In a wide
range of shear rates and surface energies, the slip length in simple ﬂuids is determined by a
function of the in-plane structure factor, contact density and temperature of the ﬁrst ﬂuid
layer near the wall [29]. The rate-dependent slip length at the interface between a polymer
melt composed of short linear chains and a crystalline substrate is also well described by the
power law function used to ﬁt the data for simple ﬂuids [27, 23]. Whether these boundary
conditions hold for higher molecular weight polymers or rough surfaces is an open question
that has to be addressed for the correct interpretation of the experimental results and
modeling ﬂuid ﬂows in microﬂuidic channels.
In this chapter, the eﬀects of shear rate and ﬂuid density on slip boundary conditions for
the ﬂow of a polymer melt past a crystalline wall are investigated by molecular dynamics

32

simulations. We will show that the rate-dependent slip length reaches a local minimum
and then increases rapidly at higher shear rates. In the shear thinning regime, the friction
coeﬃcient at the liquid/solid interface follows a power law decay as a function of the slip
velocity. In a wide range of shear rates and polymer densities, the ratio of the slip length
to the shear viscosity is determined by the surface induced structure in the ﬁrst ﬂuid layer.
These results suggest that the correlation between the ﬂuid structure near the ﬂat wall and
the slip length is determined by a universal combination of parameters for both polymer
melts and simple ﬂuids [29].
This chapter is organized as follows. The molecular dynamics simulation model and
parameter values are described in the next section. The results for the rate dependence of
the slip length and the analysis of the ﬂuid structure near the solid wall are presented in
Section 5.3. The conclusions are given in the last section.

3.2

Details of molecular dynamics simulations

The polymer melt is conﬁned between two atomically ﬂat walls and is subject to planar
shear in the x direction (see Fig. 3.1 for a schematic representation). The ﬂuid monomers
ˆ
interact through the pairwise Lennard-Jones (LJ) potential
VLJ (r) = 4 ε

σ 12
σ 6
−
,
r
r

(3.1)

with a cutoﬀ distance rc = 2.5 σ, where ε and σ represent the energy and length scales of
the ﬂuid phase. The total number of monomers is Nf = 6000. The wall-ﬂuid interaction
parameters of the LJ potential are set to εwf = 0.9 ε and σwf = σ throughout the study.
The polymer melt is modeled as a collection of bead-spring chains of N = 20 LJ monomers
interacting through the FENE (ﬁnitely extensible nonlinear elastic) potential [59]
k 2
2
VF EN E (r) = − ro ln[1 − r2 /ro ],
2

(3.2)

with Kremer and Grest parameters k = 30 εσ −2 and ro = 1.5 σ [53]. This choice of parameters
prevents polymer chains from unphysical crossing or bond breaking [54].
The system was coupled to an external heat bath through a Langevin thermostat [54]
with a random force and a damping term with a small value of the friction coeﬃcient
33

Γ = 1.0 τ −1 [73]. The thermostat was only applied to the equation of motion for a ﬂuid
monomer along the y axis to avoid a bias in the shear ﬂow direction [21]. The equations of
ˆ
motion for the x, y and z components are given by
ˆ ˆ
ˆ
(3.3)

i=j

∂Vij
,
∂xi

(3.4)

i=j

∂Vij
+ fi ,
∂yi

(3.5)

i=j

∂Vij
,
∂zi

m¨i = −
x
m¨i + mΓyi = −
y
˙
m¨i = −
z

where fi is a random force with zero mean and variance fi (0)fj (t) = 2mkB T Γδ(t)δij , the
temperature is T = 1.1 ε/kB and kB is the Boltzmann constant. The equations of motion
for ﬂuid monomers and wall atoms were integrated using the ﬁfth-order gear-predictor algorithm [56] with a time step

t = 0.002 τ , where τ =

mσ 2 /ε is the characteristic time of

the LJ potential.
Each wall was constructed out of 576 LJ atoms, which formed two (111) layers of the
face-centered cubic (F CC) lattice with the density ρw = 1.40 σ −3 .

The corresponding

nearest-neighbor distance between equilibrium positions of the wall atoms in the xy plane is
d = 1.0 σ. The dimensions of the Couette cell (xyz) were held ﬁxed to 20.86 σ × 12.04 σ × h,
where h is the channel height (see Table 3.1). The MD simulations were performed at a
constant density ensemble. The ﬂuid density is deﬁned as a ratio of the total number of
monomers to the volume accessible to the ﬂuid, i.e., 0.5 σw away from the F CC lattice planes
in the z direction. The range of ﬂuid densities considered in this study is 0.86
ˆ

ρ σ3

1.02.

The pressures of the system at equilibrium are given in Table 3.1 for each value of the ﬂuid
density. Periodic boundary conditions were employed along the x and y directions parallel
ˆ
ˆ
to the conﬁning walls.
ρ [ σ −3 ]
P [ ε σ −3 ]
h [σ]

0.86 0.88 0.91 0.94 1.00
0.0
0.5
1.0
2.0
4.0
28.93 28.18 27.27 26.32 24.85

1.02
5.0
24.36

Table 3.1. The ﬂuid pressure P at equilibrium, i.e., with the stationary upper wall, and the
channel height h (deﬁned as a distance between the inner F CC planes) as a function of the
ﬂuid density.

34

The wall atoms were attached to their equilibrium F CC lattice positions by the harmonic
spring with the potential Vsp = 1 κ r2 . The spring stiﬀness κ = 1200 ε/σ 2 is large enough so
2
that the mean-square displacement of the wall atoms is less than the Lindemann criterion
for melting, i.e., δu2 /d 2

0.023 (see Ref. [74]). In a recent paper on shear ﬂow of

simple ﬂuids [37], it was shown that the slope of the rate-dependent slip length remains
approximately the same for the thermal walls with κ = 1200 ε/σ 2 and rigid F CC walls. The
friction term and the random force were also added to the x, y and z components of the
ˆ ˆ
ˆ
wall atoms equations of motion
(3.6)

i=j

∂Vij
∂Vsp
−
+ fi ,
∂xi
∂xi

(3.7)

i=j

∂Vij
∂Vsp
−
+ fi ,
∂yi
∂yi

(3.8)

i=j

∂Vij
∂Vsp
−
+ fi ,
∂zi
∂zi

mw xi + mw Γxi = −
¨
˙
mw yi + mw Γyi = −
¨
˙
mw zi + mw Γzi = −
¨
˙

where the mass of the wall atoms is mw = 10 m and the summation is performed only over
the ﬂuid monomers within the cutoﬀ distance rc = 2.5 σ.
The simulations began with an equilibration period of about 5×104 τ with both walls being
at rest. Then, the steady-state shear ﬂow was generated by translating the upper wall with
a constant velocity U in the x direction, while the lower wall always remained stationary
ˆ
(see Fig. 3.1). The ﬂuid velocity proﬁles were computed by averaging the instantaneous
monomer velocities in the x direction within bins of thickness ∆z = 0.2 σ for a time interval
ˆ
up to 4×105 τ at the lowest shear rates. The density proﬁles were averaged for a time period
5 × 104 τ within bins of thickness of only ∆z = 0.01 σ to accurately resolve the ﬂuid structure
near the walls [29]. The maximum value of the Reynolds number is Re ≈ 10.5, which was
estimated from the maximum ﬂuid velocity diﬀerence across the channel, the shear viscosity,
and the channel height. The small values of the Reynolds number correspond to laminar
ﬂow conditions even at the highest shear rates considered in this study.

35

3.3
3.3.1

Results
Fluid density and velocity proﬁles

The ﬂuid structure in the direction perpendicular to a ﬂat surface usually consists of several
distinct layers on the molecular scale [75]. Examples of the averaged monomer density
proﬁles near the lower thermal wall are presented in Fig. 3.2 for three values of the upper
wall velocity. The highest peak is associated with a layer of ﬂuid monomers (each connected
to a polymer chain) in contact with the wall atoms. The contact density ρc corresponds
to the maximum value of the peak. The density oscillations gradually decay away from
the wall to a uniform proﬁle with bulk density ρ = 0.91 σ −3 . The amplitude of the density
oscillations is reduced at higher values of the slip velocity.
Figure 3.3 shows the averaged velocity proﬁles for the same values of the upper wall speed
as in Fig. 3.2 and the ﬂuid density ρ = 0.91 σ −3 . For a given U , the relative slip velocity is
the same at the upper and lower walls. Surprisingly, the scaled slip velocity is smaller at
the intermediate speed of the upper wall U = 1.0 σ/τ . Because of the thermal ﬂuctuations,
longer averaging times (up to 4 × 105 τ ) were required to accurately resolve the ﬂuid velocity
proﬁle for U = 0.05 σ/τ . A small curvature in the velocity proﬁle is evident for the largest
velocity of the upper wall U = 5.0 σ/τ . The deviation from the linearity might be related to
the heating up of the ﬂuid near the walls (see discussion below). The slope of the velocity
proﬁle (used for calculation of the slip length and the shear rate) was computed from the
central part of the velocity proﬁle excluding the data from the region of about 2 σ near the
walls.
We comment that at small shear rates, γτ
˙

0.01, and higher ﬂuid densities, ρ

1.04 σ −3 ,

the ﬂuid velocity proﬁles acquired a pronounced curvature within a distance of about 4 σ
from the walls. The linearity of the ﬂuid velocity proﬁles is restored at higher shear rates.
This eﬀect will be examined in a separate study.

3.3.2

Rate dependence of the slip length and shear viscosity

In the MD study by Thompson and Troian [23] on shear ﬂow of Newtonian liquids past
atomically ﬂat surfaces, the slip length was found to increase nonlinearly with the shear
36

rate. The data for diﬀerent wall densities and weak wall-ﬂuid interactions were ﬁtted well
by the power law function
Ls (γ) = Lo (1 − γ/γc )−0.5 ,
˙
˙ ˙
s

(3.9)

where γ˙c and Lo are the parameters used for normalization. The rate-dependent slip length
s
in the Poiseuille ﬂow of simple ﬂuids was also well described by the function [Eq. (3.9)] for
incommensurate structures of the liquid/solid interface and weak surface energy [29].
The nonlinear behavior of the slip length can be interpreted in terms of the velocitydependent shear stress at the wall [76]. By deﬁnition, the ﬂuid slip velocity Vs is related to
the slip length as Vs = γLs . On the other hand, the shear stress through any plane parallel
˙
to the walls is the same in the steady-state shear ﬂow, and it is equal to the ﬂuid viscosity
multiplied by the shear rate, i.e., σxz = µγ. At the liquid/solid interface the viscous shear
˙
stress becomes σxz = kVs , where k is the friction coeﬃcient per unit area. The functional
form given by Eq. (3.9) implies that the friction coeﬃcient, k = µ/Ls , can be expressed as
a function of the slip velocity as follows:
k(Vs ) = C1 ( C2 + Vs2 − Vs ),

(3.10)

where C1 = µ/2γc (Lo )2 and C2 = (2Lo γc )2 , and µ is the shear-independent viscosity [23,
˙
s
s˙
27].
In the previous MD study [27], it was shown that the slip length increases monotonically
with shear rate for polymer melts with short linear chains N

16. The data for the slip

length were well ﬁtted by the power law function given by Eq. (3.9) for shear rates exceeding
γ
˙

5×10−3 τ −1 . Although the ﬂuid viscosity exhibited a slight shear thinning behavior for

chain lengths N

8, no minimum in the rate dependence of the slip length was observed.

In the present study, the shear thinning eﬀect is more pronounced because the polymer melt
consists of longer chains and the simulations are performed at higher ﬂuid densities.
Figure 3.4 shows the rate behavior of the ﬂuid viscosity, which was computed from the
Kirkwood relation for the shear stress [59] in steady-state shear ﬂow. In the range of
accessible shear rates, the transition from the Newtonian regime to the shear thinning is
more evident at lower ﬂuid densities. At higher shear rates, the bulk viscosity approximately
follows a power law decay with the exponent −0.37 (see Fig. 3.4). A similar shear thinning
37

behavior of polymer melts was reported in the previous MD studies [60, 61]. The error bars
are larger at lower shear rates because of an increasing role of thermal ﬂuctuations [77].
The dynamic response of the slip length with increasing shear rate is shown in Fig. 3.5 for
the indicated values of the polymer density. The rate dependence of the slip length exhibits
a distinct local minimum, which is followed by a rapid increase at higher shear rates. The
minimum is shifted to lower shear rates and becomes more shallow at higher ﬂuid densities.
The slip length tends to saturate to a constant value at the lowest shear rates for the ﬂuid
density ρ = 0.86 σ −3 . Typical error bars are included for three points corresponding to the
ﬂuid velocity proﬁles shown in Fig. 3.3. In the range of shear rates from the local minimum
up to the highest γ, the data for the slip length (as well as the ratio Ls /µ) cannot be ﬁtted
˙
by the power law function given by Eq. (3.9). For each value of the ﬂuid density, the highest
shear rate corresponds to the maximum sustainable shear stress at the interface [29].
The data for the rate-dependent slip length and ﬂuid viscosity are replotted in Fig. 3.6
in terms of the friction coeﬃcient at the liquid/solid interface and the slip velocity. The
constant value of the friction coeﬃcient at small slip velocities indicates that the ratio of
the slip length and viscosity is rate-independent. At higher ﬂuid densities (or pressures),
the friction coeﬃcient is larger, and the transition to the power law regime is shifted to
smaller slip velocities. The transition point approximately corresponds to the location of
the minimum in the shear rate dependence of the slip length shown in Fig. 3.5. The data for
the lowest ﬂuid density ρ = 0.86 σ −3 can be fairly well ﬁtted by Eq. (3.10) [ see Fig. 3.6 ], but
the ﬁt becomes worse at higher ﬂuid densities (not shown). The simulation results indicate
that the nonmonotonic rate behavior of the slip length for polymer melts can be ascribed
to the velocity dependence of the friction coeﬃcient, which undergoes a transition from a
nearly constant value to a decreasing function of the slip velocity.
Furthermore, the data for diﬀerent ﬂuid densities can be collapsed onto a single master
curve by scaling both axes (see Fig. 3.7). The values of the normalization parameters for
the friction coeﬃcient and the slip velocity are listed in the inset. The decay at large
slip velocities approximately follows the power law function with the exponent −0.7 (see
Fig. 3.7). The onset of the power law decay depends on the polymer density and might occur
when the time scale of the in-plane diﬀusion of the ﬂuid monomers over the lattice spacing

38

d is comparable to d/Vs . We also note that the data for the viscosity and slip length at the
interface between a polymer melt (with linear chains N = 20 and density ρ = 0.88 σ −3 ) and
a rigid F CC wall (with density ρw = 1.94 σ −3 and εwf = ε) can be normalized to follow the
master curve shown in Fig. 3.7 [93].

3.3.3

Analysis of the ﬂuid structure in the ﬁrst layer

The surface induced order in the adjacent ﬂuid layer is described by the structure factor at
the reciprocal lattice vectors of the crystal wall [21]. The structure factor is deﬁned as
S(k) =

1
N

e

i k·rj 2

,

(3.11)

j

where rj = (xj , yj ) is the position vector of the j-th monomer and the summation is performed over N monomers within the ﬁrst layer. The friction force between a layer of
adsorbed monomers and a solid substrate is proportional to the peak value of the in-plane
structure factor estimated at the ﬁrst reciprocal lattice vector [78, 79]. Previous MD studies
of simple ﬂuids have also reported a strong correlation between the slip length (which is
inversely proportional to the friction coeﬃcient at the liquid/solid interface) and the surface
induced order in the ﬁrst ﬂuid layer [21, 121, 29].
Figure 3.8 shows the averaged structure factor in the ﬁrst ﬂuid layer moving with a
constant velocity against the stationary lower wall. The circular ridge at the wavevector |k| ≈ 2π/σ is attributed to the short-range ordering of the ﬂuid monomers. The periodic
potential of the crystalline substrate induces several sharp peaks in the ﬂuid structure factor, which are reduced at higher slip velocities. The ﬁrst reciprocal lattice vector in the
direction of shear ﬂow is G1 = (7.23 σ −1 , 0). Although the nearest-neighbor distance between equilibrium lattice sites in the xy plane is equal to the diameter of the LJ monomers,
the surface induced order in the ﬁrst ﬂuid layer is frustrated by the topological constraints
associated with packing of polymer chains near the surface.
In the recent MD study on shear ﬂow of simple ﬂuids past rigid F CC walls [29], it was
shown that the slip length scales according to
Ls ∼ [ T1 /S(G1 ) ρc ] 1.44 ,

39

(3.12)

where G1 is the the ﬁrst reciprocal lattice vector, T1 and ρc are the temperature and the
contact density of the ﬁrst ﬂuid layer respectively. This relation was found to hold in a wide
range of shear rates and for weak wall-ﬂuid interactions 0.3

εwf / ε

1.1. It is expected,

however, that the scaling relation given by Eq. (3.12) is not valid at higher surface energies
because the ﬁrst ﬂuid layer is locked to the solid wall and, as a consequence, the slip length
is negative [21]. Note that the surface induced peaks in the structure factor are much higher
for polymers (see Fig. 3.8) than for simple ﬂuid [29]. Since the viscosity is independent of
shear rate for simple ﬂuids, the scaling relation Eq. (3.12) also predicts the dependence of
the friction coeﬃcient on microscopic parameters of the adjacent ﬂuid layer.
We computed the parameters entering Eq. (3.12) in the ﬁrst ﬂuid layer for the polymer
melt. The ﬂuid temperature was estimated from the kinetic energy
m
kB T =
3N

N

[˙ i − v(ri )]2 ,
r

(3.13)

i=1

˙
where ri is the instantaneous velocity of the ﬂuid monomer and v(ri ) is the local average ﬂow
velocity. Averaged temperature proﬁles are shown in Fig. 3.9 for the same ﬂow conditions
as in Figs. 3.2 and 3.3. At low shear rates γ
˙

0.01 τ −1 , the ﬂuid temperature is equal to its

equilibrium value of T = 1.1 ε/kB determined by the Langevin thermostat. As in the case of
simple ﬂuids [29, 37], with further increase of the shear rate, the polymer melt heats up and
the temperature proﬁle becomes non-uniform across the channel. The ﬂuid temperature is
higher near the walls because of the large slip velocity, which is comparable to the thermal
2
velocity, vT = kB T/m, at high shear rates. The non-uniform temperature distribution might

be related to the slight curvature in the velocity proﬁle for the top wall speed U = 5.0 σ/τ
shown in Fig. 3.3. We can also notice that at high shear rates, the temperature in the y
ˆ
direction, in which the Langevin thermostat is applied, is slightly smaller than its value in
the x and z directions (see inset of Fig. 3.9). This diﬀerence implies that the kinetic energy
ˆ
ˆ
in the y direction dissipates faster than the energy transfer from the other directions. In
ˆ
principle, it is important to investigate whether the choice of the thermostat, e.g. dissipative
particle dynamics thermostat, has an eﬀect on the velocity and temperature proﬁles at high
shear rates [80]. In the present study, the temperature in the ﬁrst ﬂuid layer was computed

40

from
T1 =

z1
z0

T (z)ρ(z)dz

z1
z0

ρ(z)dz,

(3.14)

where the limits of integration (z0 = −14.4 σ and z1 = −13.8 σ) are given by the width of
the ﬁrst peak in the density proﬁle (e.g. see Fig. 3.2).
The ratio of the slip length to viscosity is plotted in Fig. 3.10 as a function of T1 /[S(G1 ) ρc ]
for diﬀerent ﬂuid densities and shear rates. Except for higher densities at low shear rates, the
data collapse onto a single curve, which approximately follows the power law function with
a slope of 0.9 (see Fig. 3.10). Although the slope is diﬀerent from the value 1.44 reported for
simple ﬂuids [29], it is remarkable that the same combination of parameters determines the
boundary conditions for both polymer melts and simple ﬂuids. The simulation results also
show that for each value of the ﬂuid density, the number of ﬂuid monomers in the ﬁrst layer
N = S(0) decreases by about 4% at the highest shear rates. The same data for the slip
length and viscosity are replotted as a function of the variable S(0)/ [S(G1 ) ρc ] in Fig. 3.11.
The collapse of the data is consistent with the previous results for the rate-independent
slip length at the interface between simple ﬂuids and crystalline walls [21]. The correlation
between the friction coeﬃcient and the ﬂuid structure in the adjacent layer presented in
Figs. 3.10 and 3.11 suggests a possible existence of a general functional relation between the
slip length, shear rate and polymer density.

3.4

Conclusions

In this chapter, the rate dependence of the slip length in the shear ﬂow of polymer melts
past atomically smooth, thermal surfaces was studied using molecular dynamics simulations.
For weak wall-ﬂuid interactions, the slip length passes through a minimum as a function of
shear rate and then increases rapidly at higher shear rates. The nonlinear rate dependence
of the slip length was analyzed in terms of the dynamic friction at the liquid/solid interface.
In a wide range of ﬂuid densities, the friction coeﬃcient (the ratio of the shear viscosity
and the slip length) undergoes a universal transition from a constant value to the power law
decay as a function of the slip velocity. The simulation results for polymer melts conﬁrm the
previous ﬁndings for simple ﬂuids [29] that the surface induced structure and the contact

41

density of the adjacent ﬂuid layer are crucial factors for determining the value of the friction
coeﬃcient.

42

U
z

h

y

Vx (z)

x

Polymer melt

Ls

Vs

Thermal wall

Figure 3.1. A schematic representation of the steady-state shear ﬂow in the Couette cell.
The upper wall is translated with a constant velocity U in the x direction. The ﬂuid slip
ˆ
velocity is determined from the relation Vs = γLs , where γ is the slope of the velocity proﬁle
˙
˙
and Ls is the slip length.

43

ρσ3

3

ρ = 0.91 σ−3
U = 0.05 σ /τ
U = 1.0 σ /τ
U = 5.0 σ /τ

2

1

0

-14

-13

z /σ

-12

-11

-10

Figure 3.2. Averaged ﬂuid density proﬁles near the stationary thermal wall with
κ = 1200 ε/σ 2 and εwf /ε = 0.9. The velocities of the upper wall are tabulated in the inset. The ﬂuid density in the center of the cell is ρ = 0.91 σ −3 . The left vertical axis denotes
the position of the liquid/solid interface at z = − 14.74 σ.

44

1

V /U
x

0.8

0.6

V /U
x

0.4

z /σ

0.2
-14

-12

-10

-8

ρ = 0.91 σ−3
U = 0.05 σ /τ
U = 1.0 σ /τ
U = 5.0 σ /τ

0.4

0.2

0

-12

-8

-4

z /σ

0

4

8

12

Figure 3.3. Average normalized velocity proﬁles, Vx / U , for the indicated values of the
upper wall velocity U and the ﬂuid density ρ = 0.91 σ −3 . Vertical axes denote the position
of the F CC lattice planes at z = − 15.24 σ and z = 12.03 σ. The dashed lines indicate the
location of liquid/solid interfaces at a distance 0.5 σw away from the F CC lattice planes,
i.e., z = − 14.74 σ and z = 11.53 σ. The inset shows the same data near the lower wall.

45

ρ = 0.86 σ−3
ρ = 0.88 σ−3
ρ = 0.91 σ−3

µ / ετσ

−3

100

ρ = 0.94 σ−3
ρ = 1.00 σ−3
ρ = 1.02 σ−3

10

0.01
.

0.001

γτ

0.1

Figure 3.4. Behavior of the ﬂuid viscosity µ/ ετ σ −3 as a function of shear rate for the
tabulated values of the polymer density. The straight dashed line with a slope −0.37 is
shown for reference. Solid curves are a guide for the eye.

46

20

ρ = 0.86 σ−3
ρ = 0.88 σ−3
ρ = 0.91 σ−3

Ls / σ

16

ρ = 0.94 σ−3
ρ = 1.00 σ−3
ρ = 1.02 σ−3

12

8

4

0

0.01
.

0.001

γτ

0.1

Figure 3.5. Slip length Ls /σ as a function of the shear rate for the polymer melt with linear
chains N = 20. Solid curves are a guide for the eye.

47

µ / Ls [ετσ−4 ]

ρ = 0.86 σ−3
ρ = 0.88 σ−3
ρ = 0.91 σ−3
ρ = 0.94 σ−3

10

1

ρ = 1.00 σ−3
ρ = 1.02 σ−3
0.01

0.1

Vs τ /σ

1

Figure 3.6. Log-log plot of the friction coeﬃcient, k = µ/Ls , as a function of the slip velocity
for the indicated values of the ﬂuid density. The dashed curve is the best ﬁt to Eq. (3.10)
with C1 = 1.77 mσ −3 and C2 = 0.71 σ 2 τ −2 . Solid curves are a guide for the eye.

48

Vs*
1.00
0.89
0.70
0.45
0.10
0.03

1

k / k*

ρ = 0.86 σ−3
ρ = 0.88 σ−3
ρ = 0.91 σ−3
ρ = 0.94 σ−3
ρ = 1.00 σ−3
ρ = 1.02 σ−3

0.1

0.01

0.1

1

Vs

10

k*
1.41
1.72
2.30
3.59
13.45
35.64

100

/ V*
s

Figure 3.7. Master curve for the friction coeﬃcient as a function of the slip velocity. The
same data as in Fig. 3.6. The values of the friction coeﬃcient k ∗ [ ετ σ −4 ] and the slip
velocity Vs∗ [ σ/τ ] used for normalization of the data are shown in the inset. The dashed
line y = [1 + (4x)2 ]−0.35 is plotted for reference.

49

Figure 3.8. Structure factor S(kx , ky ) in the ﬁrst ﬂuid layer for the ﬂuid density ρ = 0.91 σ −3 .
The upper wall velocity is U = 0.05 σ/τ (a) and U = 5.0 σ/τ (b). The same ﬂow conditions
as in Fig. 3.3.

50

2.4

T kB /ε

1.8

T kB /ε

2

1.6

U = 5.0 σ /τ

1.6

Tx
Ty
Tz

1.4

z /σ

1.2
-14

-12

-10

1.2

ρ = 0.91 σ−3
U = 0.05 σ /τ

0.8

-12

-8

-4

z /σ

U = 1.0 σ /τ
U = 5.0 σ /τ
0

4

8

12

Figure 3.9. Temperature proﬁles for the indicated values of the upper wall speed U and
the ﬂuid density ρ = 0.91 σ −3 . Vertical dashed lines denote the position of the liquid/solid
interface. The inset shows the x, y and z components of the temperature proﬁle near the
ˆ ˆ
ˆ
stationary lower wall when the velocity of the upper wall is U = 5.0 σ/τ .

51

10

Ls / µ [σ 4 / ετ]

Slope = 0.9
1

ρ = 0.86 σ−3
ρ = 0.88 σ−3
ρ = 0.91 σ−3
ρ = 0.94 σ−3
ρ = 1.00 σ−3
ρ = 1.02 σ−3

0.1

0.01

0.01

0.1
3 −1

0.3

T1 kB /ε [ S(G1)ρcσ ]

Figure 3.10. Log-log plot of Ls /µ as a function of the ratio T1 kB /ε [S(G1 ) ρc σ 3 ]−1 for
the indicated values of the polymer density. The solid line with a slope 0.9 is plotted for
reference.

52

10

Ls / µ [σ 4 / ετ]

Slope = 1.15
1

ρ = 0.86 σ−3
ρ = 0.88 σ−3
ρ = 0.91 σ−3
ρ = 0.94 σ−3
ρ = 1.00 σ−3
ρ = 1.02 σ−3

0.1

0.01

2

10

3 −1

50

S(0) [ S(G1)ρcσ ]

Figure 3.11. Log-log plot of the ratio Ls /µ as a function of the variable S(0)/ [S(G1 ) ρc σ 3 ].
The values of the ﬂuid density are tabulated in the inset. The same data as in Fig. 3.10.
The straight solid line with a slope 1.15 is plotted as a reference.

53

CHAPTER 4
The eﬀective slip length and vortex
formation in laminar ﬂow over a
rough surface
4.1

Introduction

An accurate ﬂow prediction in micro-channels is important for the optimal design and
fabrication of microﬂuidic devices whose applications range from medicine to biotechnology [81, 82]. The boundary conditions and the surface topology are major factors aﬀecting
the ﬂow pattern near the solid boundary and the formation of recirculation zones. The
ﬂow separation at rough surfaces can modify the wall shear stress distribution or initiate
instability towards turbulence. In microﬂuidic channels, the vortex ﬂow enhances the mixing eﬃciency [83, 84] and promotes convective heat transfer [85, 86, 87]. In cardiovascular
systems, the separation region at the entrance of branching vessels may trap lipid particles
which could lead to arterial diseases [88, 89, 90]. In the present study we examine the role
of slip boundary condition in determining the ﬂow properties near rough surfaces including
the separation phenomena, and distribution of pressure and shear stress along the surface.
Although the validity of the no-slip boundary condition is well accepted at the macroscopic level, recent experiments [5, 16, 9, 19, 91, 14] and molecular dynamics (MD) simulations [21, 23, 25, 92, 27, 29] reported the existence of a boundary slip in microﬂows. The

54

model ﬁrst proposed by Navier relates the slip velocity to the rate of shear via the proportionality coeﬃcient, the so-called slip length. The MD simulations are particularly suitable
for examining the inﬂuence of molecular parameters on the microscopic slip length at the
liquid/solid interface. The advantage of the MD simulations is that a detailed ﬂow analysis
can be performed at the molecular level while the explicit speciﬁcation of the boundary
conditions is not required. In contrast to description of the ﬂow near boundary by means of
microscopic slip length, it is convenient to characterize the ﬂow over macroscopically rough
surfaces by the eﬀective slip length, which is deﬁned as a distance from the level of the
mean height of the surface roughness to the point where linearly extrapolated bulk velocity
proﬁle vanishes. Recent MD studies have demonstrated that the eﬀective slip length in a
ﬂow of simple ﬂuids [35] and polymer melts [93] over a wavy surface agrees well with hydrodynamic predictions [51, 52] when the corrugation wavelength is larger than approximately
thirty molecular diameters.
The inﬂuence of surface roughness on ﬂuid ﬂow with either local no-slip or zero shear stress
(i.e. perfect slip) boundary conditions has been extensively studied in the past decades [43,
62, 94, 95, 63, 96, 97], see also a review section in [35]. Analytical calculations have shown
that in a shear ﬂow over a corrugated surface with microscopic no-slip or zero shear stress
conditions, the eﬀective boundary slip is insigniﬁcant macroscopically [43, 94]. The eﬀective
no-slip boundary plane is located at the intermediate position between crests and valleys of
the rough surface when the no-slip condition is imposed along the solid boundary [62, 63,
98, 14, 36]. For an arbitrary surface roughness with small amplitudes, the slip coeﬃcient
in the Navier model is proportional to the average amplitude of the wall roughness and
depends on the position of the origin of the coordinate system [95]. Applying the no-slip
boundary condition along the wavy surface, Tuck and Kouzoubov [63] have demonstrated
that the eﬀective slip length is inversely proportional to the corrugation wavelength and
quadratically proportional to the amplitude of the surface roughness. However, the series
expansion method used in [63] fails at large wavenumbers, ka

0.5, when a backﬂow

appears inside the grooves of the substrate. The eﬀective slip length for a ﬂow above the
surface with deep corrugations only weakly depends on the depth of the grooves [62, 97, 93].
Despite considerable analytical eﬀorts, the relation between the vortex ﬂow structure in

55

deep grooves and the eﬀective slip length has not yet been systematically investigated.
The laminar ﬂow separation at the corrugated surface with the local no-slip boundary
conditions depends on the depth of the grooves and the Reynolds number [99, 100, 101, 102,
103, 104, 105, 106]. In a creeping ﬂow over a sinusoidal surface, the ﬂow circulation appears
in suﬃciently deep grooves and, as the corrugation amplitude increases, the vortex grows
and remains symmetric [104, 107, 108]. With increasing Reynolds number, the vortex ﬂow
forms even in shallow grooves, the circulation region expands, and the center of vorticity is
displaced upstream [99, 102, 103, 105]. In the limit of small-scale surface roughness and for
no-slip boundary conditions, the apparent slip velocity at the mean surface becomes more
negative as the Reynolds number increases [63]. A noticeable change in the eﬀective slip
length was also observed at Re

100 for laminar ﬂow over deep grooves when the local slip

length is comparable to the corrugation amplitude [93]. However, the inﬂuence of the local
slip condition at the curved boundary on the vortex ﬂow formation has not been considered
at ﬁnite Reynolds numbers.
This chapter is focused on investigation of the eﬀects of local slip boundary conditions and
the Reynolds number on the ﬂow structure near periodically corrugated surfaces and the
eﬀective slip length. We will show that for the Stokes ﬂow with the local no-slip boundary
condition, the eﬀective slip length decreases with increasing corrugation amplitude and
the ﬂow circulation develops in suﬃciently deep grooves. In the presence of the local slip
boundary condition along the rough surface, the eﬀective slip length increases and the size of
the vortex is reduced but its structure remains symmetrical. The analysis of the numerical
solution of the Navier-Stokes equation with the local slip condition shows that the inertial
eﬀects promote the asymmetric vortex ﬂow formation and reduce the eﬀective slip length.
This chapter is organized as follows. The details of a continuum model and the implementation of the local slip boundary conditions are described in the next section. The
analytical results for the Stokes ﬂow over a wavy surface by Panzer et al. [52] are summarized
in Sec. 4.3.1. The analysis of the eﬀective slip length and the ﬂow structure is presented in
Sec. 4.3.2 for the no-slip case and in Sec. 4.3.3 for the ﬁnite microscopic slip. The eﬀect of
Reynolds number on the eﬀective slip ﬂow over a periodically corrugated surface is studied
in Sec. 4.3.4. A brief summary is given in the last section.

56

4.2

Details of numerical simulations

The two-dimensional incompressible and steady Navier-Stokes (NS) equation is solved using
the ﬁnite element method. The computational setup consists of a viscous ﬂuid conﬁned
between an upper ﬂat wall and a lower sinusoidal wall (see Fig. 4.1). The corrugation
wavelength of the lower wall is set to λ and is equal to the system size in the x direction.
ˆ
The upper wall is located at h = λ above the reference line at a = 0, which is deﬁned as
the mean height of the surface roughness. The continuity and NS equations are
· u = 0,
ρ (u ·

u) = − p + µ 2 u,

(4.1)
(4.2)

ˆ
where u = u ˆ + v j is the velocity vector in the cartesian coordinate system, p is the pressure
i
ﬁeld, ρ is the ﬂuid density, and µ is the Newtonian viscosity.
The penalty formulation is employed to avoid decoupling between the pressure and velocity ﬁelds [109]. In this method, the continuity equation is replaced with a perturbed
equation
p
·u=− ,
Λ

(4.3)

where Λ is the penalty parameter, which ensures the incompressibility condition. Thus, the
modiﬁed momentum equations in the x and z directions are
ˆ
ˆ
ρ (u ·

u) = Λ (

· u) + µ 2 u,

(4.4)

ρ (u ·

v) = Λ (

· u) + µ 2 v.

(4.5)

The advantage of the penalty formulation is the elimination of pressure and the continuity
equation. The penalty parameter Λ must be large enough so that compressibility errors are
minimal. The upper bound of Λ is determined from the condition that the viscous eﬀects
are not destroyed by the machine precision [109, 110]. The penalty parameter Λ should be
chosen according to the rule
Λ = c max{µ, µRe},

(4.6)

where Re is the Reynolds number and the constant c is recommended to be about 107 for
computations with double-precision 64 bit words [110].
57

The Galerkin formulation of Eq. (4.4) and Eq. (4.5) can be explicitly rewritten as
∂Nj
∂Nj
+ vi uj
¯
∂x
∂z

+

∂Nj
∂Nj
+ vi vj
¯
ρNi ui vj
¯
∂x
∂z
Ω

+

Ω

ρNi ui uj
¯

∂Nj
∂Ni ∂Nj
uj +
v dΩ +
∂x ∂x
∂z j
Ω
∂Ni ∂Nj ∂Ni ∂Nj
µ
+
uj dΩ = RHSx , (4.7)
∂x ∂x
∂z ∂z
Ω
Λ

∂Nj
∂Ni ∂Nj
uj +
v dΩ +
∂z ∂x
∂z j
Ω
∂Ni ∂Nj ∂Ni ∂Nj
vj dΩ = RHSz , (4.8)
µ
+
∂x ∂x
∂z ∂z
Ω
Λ

where Ni is the weight function, Nj , uj , vj are the node shape function and the velocities
in each element, ui and vi are the lagged velocities, and the right hand side (RHSx , RHSz )
¯
¯
terms include the boundary velocities.
In our simulation, the bilinear quadrilateral elements (i, j = 1, 2, 3, 4) with non-orthogonal
edges are transformed to straight-sided orthogonal elements by introducing the natural
coordinates ξ = ξ(x, z) and η = η(x, z). The shape functions Ni in the natural coordinate
system are deﬁned as
Ni =

(1 + ξi ξ)(1 + ηi η)
4

i = 1, .., 4

(4.9)

where ξi and ηi are the corner points of each element (see Fig. 4.2).
In the next step, Eq. (4.7) and Eq. (4.8) are integrated numerically using four-point
Gaussian quadrature [57]. The ﬁnal system of equations is constructed as follows:
ρ[K1 ] + Λ[K2 ] + µ[K3 ]

u
v

=

RHSx
,
RHSz

(4.10)

where the terms RHSx and RHSz contain the velocities at the boundary nodes.
The boundary conditions must be speciﬁed at the inlet, outlet and upper and lower walls
of the Couette cell. The periodic boundary conditions are imposed at inlet and outlet along
the x direction. A ﬁnite slip is allowed along the lower wall while the boundary condition at
ˆ
the upper wall is always no-slip. In the local coordinate system (spanned by the tangential
t and normal n vectors), the ﬂuid velocity along the lower wavy wall is computed from
ut = L0 [(n ·

)ut + ut /R(x)],

(4.11)

where ut is the tangential component of u = ut t + un n, L0 is the intrinsic (or microscopic)
slip length at the ﬂat surface, and R(x) is the local radius of curvature [52]. The radius of
58

curvature is positive for concave and negative for convex regions. The Navier slip condition
for a ﬂat wall is recovered from Eq. (4.11) when R(x) → ∞. The eﬀective slip length Leﬀ at
the corrugated lower wall is obtained by extrapolating the linear part of the velocity proﬁle
(0.45

z/h

0.9) to zero velocity with respect to the reference line a = 0.

The simulation begins by setting the no-slip boundary condition at the upper and lower
walls as an initial guess. Once Eq. (4.10) is solved, the ﬂuid velocities at the lower boundary
are updated using Eq. (4.11). This iteration is repeated until the solution is converged to a
desired accuracy. The convergence rate of the solution remains under control by using the
under-relaxation value 0.001 for the boundary nodes. The results presented in this chapter
are obtained with the grid resolution 150 × 150 in the x and z directions, respectively. In
ˆ
ˆ
order to check the accuracy of the results, several sets of simulations were also carried out
with a ﬁner grid 180 × 180. The maximum relative error of the eﬀective slip length due
to the grid size is Leﬀ /h = 0.003. The converged solution of the Navier-Stokes equation
satisﬁes the following boundary condition:
ut = Ls (x)

∂ut
,
∂n

1
1
1
=
−
,
Ls (x)
L0 R(x)

(4.12)

where Ls (x) is the local slip length in the presence of surface curvature [52].
The accuracy of the numerical solution is checked by the normalized average error, which
is deﬁned as

Np

error =
i=1

| un − un+1 |
i
i
| un+1 |
i

/Np ,

(4.13)

where Np is the total number of computational nodes, un is the velocity at the node i and
i
time step n, and un+1 is the velocity in the next time step. The typical value of the error
i
in the converged solution is less than 10−9 . Throughout the study, the results are presented
in the non-dimensional form. The length scale, shear rate, shear stress, and velocity are
∗
˙
normalized by h, γ ∗ , τw , and U ∗ , respectively, where γ ∗ is the shear rate in the case of
˙
∗
˙
˙
no-slip boundary condition at the ﬂat upper and lower walls, and τw = µγ ∗ and U ∗ = hγ ∗ .

59

4.3
4.3.1

Results
Analytical solution of the Stokes equation for viscous ﬂow
over a wavy wall

The eﬀect of small periodic surface roughness on the eﬀective slip length has been previously
investigated for pressure-driven ﬂows in a Couette geometry [51, 52]. The analytical solution
of the Stokes equation with boundary conditions Eq. (4.12) at the wavy wall with amplitude
a and wavelength λ was obtained for two limiting cases and small ka. For L0 /λ

1 the

eﬀective slip length is given by
lim Leﬀ = L0 − ka2 ω◦ (ka),

kL0 →0

while in the limit of L0 /λ

(4.14)

1 it reduces to
1
k 3 a2 −1
lim Leﬀ =
+
,
L0 ω∞ (ka)
kL0 →∞

(4.15)

where the functions ω◦ (ka) and ω∞ (ka) are deﬁned as
ω◦ (ka) =

1 − 1/4(ka)2 + 19/64(ka)4 + O[(ka)6 ]
,
1 + (ka)2 − 1/2(ka)4 + O[(ka)6 ]

(4.16)

ω∞ (ka) =

1 − 5/4(ka)2 + 61/64(ka)4 + O[(ka)6 ]
.
1 + (ka)2 − 1/2(ka)4 + O[(ka)6 ]

(4.17)

and

An approximate analytical expression for the eﬀective slip length that interpolates between
the two bounds Eq. (4.14) and Eq. (4.15) is given by
L0 ω∞ (ka) − ka2 ω0 (ka)/(1 + 2kL0 )
Leﬀ =
,
1 + k 3 a2 L0
with the range of applicability ka

(4.18)

0.5. For larger wavenumbers ka > 0.5, the function

ω∞ (ka) overestimates the numerical solution and the interpolated formula Eq. (4.18) does
not apply [52].

4.3.2

Flow over a rough surface with the local no-slip boundary
condition

In this section, the Stokes equation with the local no-slip condition at the upper and lower
walls is solved numerically to study the eﬀect of corrugation amplitude on the eﬀective
60

slip length. The velocity proﬁles, averaged over the period of corrugation λ, are plotted
in the inset of Fig. 4.3 for several values of wavenumber ka = 2 πa/λ. As ka increases, the
normalized velocity proﬁles remain linear in the bulk region and become curved near the
lower corrugated wall. The linear part of the velocity proﬁles is used to compute the eﬀective
slip length, which is plotted as a function of wavenumber in Fig. 4.3. With increasing
corrugation amplitude of the lower wall, the eﬀective slip length decays monotonically and
becomes negative, indicating that the eﬀective no-slip boundary is shifted into the ﬂuid
domain. For ka

1, the numerical results agree well with the analytical solution Eq. (4.18)

denoted by the solid line in Fig. 4.3. The deviation from the analytical solution becomes
signiﬁcant at larger wavenumbers where the streamlines extracted from the Stokes solution
indicate the presence of backﬂow at the bottom of the valley.
In order to investigate the ﬂow behavior above the sinusoidal surface, the shear stress and
pressure along the lower wall were computed from the solution of the Stokes equation. In
the presence of surface curvature the wall shear stress τw has two components
τw = µ

∂ut
+ ut /R(x)
,
∂n
w

(4.19)

where ∂ut /∂n is the normal derivative of the tangential velocity ut , and R(x) is the local
radius of curvature. In the case of no-slip boundary condition (ut = 0), the local shear
∂u

stress at the wall is reduced to τw = µ ∂nt . The normalized shear stress along the lower
corrugated wall is plotted in Fig. 4.4 for diﬀerent corrugation amplitudes. The maximum
value of the shear stress is located at the peak of the surface corrugation (x/λ = 0.25)
and it increases with increasing amplitude, which is consistent with the results of previous
analytical studies of a laminar ﬂow over a wavy wall [111, 100]. The ﬂuid tangential velocity
near the boundary is proportional to the wall shear stress shown in Fig. 4.4. Therefore, the
tangential velocity is also maximum above the peak and, as the ﬂow moves downstream, it
decelerates and the velocity becomes zero inside the valley at suﬃciently large amplitudes.
For ka

0.79, the shear stress proﬁles intersect the dashed line (τw = 0) at two points and a

clockwise ﬂow circulation develops inside the valley. As the corrugation amplitude increases,
the intersection points move away from each other and the ﬂow recirculation region becomes
larger. These results are in agreement with previous estimates of the critical wavenumber

61

ka ≈ 0.77 for the onset of ﬂow separation in suﬃciently thick ﬁlms [107, 108].
The pressure along the lower wavy wall is plotted in Fig. 4.5 for the same amplitudes as
in Fig. 4.4. The value P ∗ used for normalization is the maximum pressure, which is located
above the wavy surface with ka = 1.12 on the left side of the peak. For each wavenumber,
the pressure along the surface is maximum on the left side of the peak, where the surface
faces the mainstream ﬂow. The surface pressure reaches its minimum value on the right
side of the peak (see Fig. 4.5). As the ﬂow moves further downstream into the valley, it
encounters an adverse pressure gradient, which becomes larger as ka increases. At large
wavenumbers ka

0.83, the ﬂow near the surface cannot overcome the combined resistance

of the viscous forces and the adverse pressure gradient, and separates from the surface at
the point where τw = 0.
The pressure contours and streamlines near the corrugated surface with wavenumber
ka = 1.12 are depicted in Fig. 4.6(a). The pressure contours indicate the presence of an
adverse pressure gradient in the region 0.3

x/λ

0.6 on the right side of the peak (see

also Fig. 4.5). The streamlines illustrate the ﬂow separation inside the valley at x/λ

0.52.

After the separation point the ﬂow near the wall reverses direction and moves against the
mainstream. The local velocity proﬁle inside the valley is shown in the inset of Fig. 4.6(a).
As the ﬂow cross-section decreases towards the peak, the ﬂow attaches to the surface at
the point where the shear stress is zero again. Note also that for ka

0.83 in Fig. 4.5, the

pressure proﬁles along the lower surface exhibit a slight drop inside the valley where the
vortex is present. After the attachment point, due to the periodic boundary condition in
the x direction, the surface pressure increases up to its maximum value, which is located on
ˆ
the left side of the peak. The nonlinearity due to inertia is absent in the Stokes ﬂow and
the vortex in the valley remains symmetrical.
We also comment that a secondary vortex is formed inside a deeper cavity (ka ≈ 2.33),
which is counter-rotating with respect to the primary vortex (not shown). Previous analytical study of a creeping ﬂow over a wavy wall demonstrated that a secondary vortex appears
at ka

2.28 when the channel width is larger than the corrugation wavelength [108]. In

the limiting case, when the cavity consists of two parallel walls, an inﬁnite sequence of
counter-rotating vortices appears between the walls [112, 113].

62

4.3.3

Eﬀect of the local slip on the ﬂow pattern near the rough
surface

Next, we present the results of the numerical solution of the Stokes equation with the
local slip condition at the lower wavy wall while the boundary condition at the upper
wall remains no-slip. The pressure contours and the streamlines are plotted in Fig. 4.6 for
several values of the intrinsic slip length L0 and the wavenumber ka = 1.12. The streamline
patterns indicate that with increasing slip length L0 , the size of vortex inside the valley is
progressively reduced and the vortex eventually disappears. Similar to the analysis in the
previous section, the pressure and shear stress along the lower wall are computed in the
presence of the local boundary slip.
The normalized pressure along the lower wavy wall is plotted in Fig. 4.7 as a function of
the slip length L0 . Similar to the case of no-slip boundary condition, the proﬁles exhibit a
maximum and a minimum in pressure on the right and left sides on the peak, respectively.
In a wide range of L0 , the diﬀerence between the locations and the magnitudes of the
extrema is barely noticeable. As the ﬂow moves down a slope (0.3

x/λ

0.5 region

in Fig. 4.7), the pressure gradient along the surface becomes positive and its magnitude
decreases with increasing values of L0 . The separation and attachment points are located
at the intersection of the shear stress proﬁles with the horizontal line (τw = 0) shown in
Fig. 4.8. In comparison to the no-slip case, the smaller combined eﬀect of the adverse
pressure gradient and the wall shear stress causes a shift of the separation point deeper into
the valley (e.g. see Fig. 4.6). As the slip length L0 increases, the separation and attachment
points move closer to each other, the vortex becomes smaller and eventually vanishes. We
also note that at the separation point both components of the shear stress ∂ut /∂n and
ut /R(x) become zero in agreement with the MRS criterion for the ﬂow separation of the
boundary layer at a moving substrate [114].
The inﬂuence of the local slip on the shape of the velocity proﬁle and the normal derivative
of the tangential velocity at the bottom of the valley is illustrated in the insets of Fig. 4.6.
The local slip length Ls (x) in Eq. (4.12) is a function of the radius of curvature R(x) and
the intrinsic slip length L0 . For all cases considered in Fig. 4.6, the condition L0 < |R(x)|

63

holds and the local slip length Ls remains positive everywhere along the lower wall, which
means that the slip velocity ut and the normal derivative ∂ut /∂n carry the same sign. It is
expected, however, that for the opposite condition, L0 > |R(x)|, the values ut and ∂ut /∂n
would have diﬀerent signs at the bottom of the valley, and the corresponding velocity proﬁle
will be qualitatively similar to the proﬁle shown in the inset of Fig. 4.6(d) but shifted by a
negative slip velocity (see next section).
The eﬀective slip length computed from the numerical solution of the Stokes equation and
Eq. (4.18) is plotted in Fig. 4.9 as a function of the intrinsic slip length L0 for wavenumbers
ka = 0.28 (a/h = 0.04) and ka = 1.12 (a/h = 0.18). For small values of L0 , the eﬀective
slip length approaches a negative value previously reported in Fig. 4.3 for the no-slip case.
As L0 increases, Leﬀ grows monotonically and appears to saturate to a constant value. The
transition of the eﬀective slip length from a growing function of L0 to a nearly constant
value occurs at larger L0 when ka decreases. For the small wavenumber ka = 0.28, the
eﬀective slip length computed from Eq. (4.18) is in a good agreement with the numerical
solution of the Stokes equation. Visual inspection of the streamline patterns indicates
that there is no backﬂow at any L0 . In the saturation regime, L0 → ∞, the wall shear
stress becomes zero everywhere along the lower boundary, the streamlines near the lower
wall follow the boundary curvature, and the eﬀective slip length in Eq. (4.15) approaches
Leﬀ /λ

1/ [2π(ka)2 ] − 9/8π. For the large wavenumber ka = 1.12, the analytical results

Eq. (4.18) overestimate Leﬀ computed from the numerical solution at L0 /h
ﬂow circulation is developed in the valley at L0 /h

0.02 and the

0.067. The vortex vanishes at the

bottom of the valley at suﬃciently large values of L0 (denoted by the vertical arrow in
Fig. 4.9), and the ﬂow streamlines are deformed to follow the boundary curvature [e.g. see
Fig. 4.6(d)]. The results for the intrinsic slip length, which determines the threshold for
the onset of the ﬂow circulation at the bottom of the groove, are summarized in Fig. 4.10.
For the wavenumbers examined in this study, ka

1.26, the numerical simulations indicate

that if the ﬂow circulation is present in the valley then the eﬀective slip length is negative
and Leﬀ increases with decreasing vortex size.

64

4.3.4

Eﬀect of the Reynolds number on the eﬀective slip length

The analysis of the Stokes equation discussed in the previous section demonstrated that
increasingly large local slip at the lower wavy wall eliminates the ﬂow circulation in the
valley and leads to a larger eﬀective slip. In this section, the inﬂuence of the inertia term in
the Navier-Stokes equation on the ﬂow pattern and the eﬀective slip length is investigated.
For the shear ﬂow with slip condition at the lower corrugated wall L0 , the Reynolds number
is deﬁned as

ρ U ∗ h (1 + Leﬀ /h)
Re =
,
µ

(4.20)

where ρ is the ﬂuid density, U ∗ is the upper wall velocity, and h (1 + Leﬀ /h) is the distance
between the upper ﬂat wall and the eﬀective no-slip boundary plane. In the case of no-slip
boundary condition at the lower ﬂat wall, Leﬀ is zero and the standard deﬁnition of the
Reynolds number is recovered, i.e., Re = ρ U ∗ h/µ.
The pressure and shear stress along the lower wavy wall are plotted in Fig. 4.11 for the
selected values of the Reynolds number and no-slip boundary conditions. As Re increases,
the adverse pressure gradient along the right side of the corrugation peak (0.3
becomes larger and the pressure drop inside the valley (0.5

x/λ

x/λ

0.5)

0.9) increases. For

each value of the Reynolds number in Fig. 4.11(b), similar to the Stokes ﬂow case, the shear
stress above the peak is maximum and, as the ﬂow moves downstream along the right side
of the peak, it decelerates and eventually separates from the surface when τw = 0. With
increasing upper wall velocity, the shear stress above the peak increases and causes the
ﬂow to decelerate faster along the right side of the peak. The separation and attachment
points (determined from the condition τw = 0) move further apart from each other and the
circulation region inside the valley expands [see Fig. 4.11(b)].
The eﬀect of the inertia term in the Navier-Stokes equation on the shape of the wall shear
stress and pressure proﬁles can be seen in Fig. 4.11. The shear stress proﬁles above the
peak are not symmetric with respect to x/λ = 0.25, which means that the ﬂow decelerates
faster on the right side of the corrugation peak than it accelerates on the left side. Also,
the average adverse pressure gradient on the left side of the peak is larger than its value
on the right side of the peak (0.3

x/λ

0.5). By increasing the upper wall velocity,

65

the separation point moves further upstream than the attachment point downstream. The
formation of asymmetric vortex ﬂow at ﬁnite Reynolds numbers is consistent with previous
ﬁndings for a ﬂow in an undulated tube [103, 99, 102, 105]. The pressure contours and
streamlines extracted from the NS equation with the no-slip boundary condition are plotted
in Fig. 4.12(a) for ka = 1.12 and Re = 79. The ﬂow streamlines in the valley indicate an
asymmetric clockwise circulation, which is larger than the ﬂow circulation region shown in
Fig. 4.6(a) for the Stokes case.
In the presence of the local slip condition along the lower corrugated wall, the size of the
vortex becomes smaller while the ﬂow structure remains asymmetric (see Fig. 4.12). The
decrease of the vortex size is similar to the case of the Stokes ﬂow shown in Fig. 4.6 and
can also be described in terms of the pressure and shear stress along the lower wall. In
Fig. 4.13 the pressure and shear stress proﬁles are plotted for the same values of the upper
wall velocity as in Fig. 4.11 but with the slip boundary condition (L0 /h = 0.25) along the
lower wall. Note that the Reynolds numbers in Fig. 4.13 are slightly larger than the values
reported in Fig. 4.11 for the same U ∗ because of the larger eﬀective slip length entering
the deﬁnition of the Reynolds number [see Eq. (4.20)]. For each value of the upper wall
velocity, the adverse pressure gradient and the wall shear stress on the right side of the
peak are smaller than in the case of the no-slip boundary condition, and, as a result, the
vortex either becomes smaller or vanishes (see Figs. 4.11 and 4.13). As the Reynolds number
increases, however, the vortex forms and expands asymmetrically to ﬁll the bottom of the
groove. The inset of the Fig. 4.12(b) demonstrates that the slip velocity at the bottom of the
valley is negative while its normal derivative ∂ut /∂n is positive, in contrast to the velocity
proﬁles shown in Fig. 4.6 for the same ka = 1.12 and smaller slip lengths L0 /h

0.08.

The eﬀective slip length is plotted in Fig. 4.14 as a function of the Reynolds number for
the selected values of the intrinsic slip length L0 and ka = 1.12. With increasing Reynolds
number, the ﬂow streamlines move away from the lower boundary and straighten out, the
slope of the normalized velocity proﬁles increases and the eﬀective no-slip boundary plane
is shifted into the ﬂuid domain. For L0 /h

0.067, the circulation is always present in

the valley and the ﬂow streamlines in the bulk of the ﬂuid do not penetrate deep into the
valley [e.g. see Fig. 4.12(b)]. For larger slip lengths L0 /h > 0.067, the ﬂow streamlines show
66

that there is no backﬂow at low Re, and the vortex is formed at the bottom of the valley
only at suﬃciently large Reynolds numbers indicated by the dashed line in Fig. 4.14. The
numerical results obtained from the solution of the Navier-Stokes equation demonstrate that
the growth or decay of the vortex as a function of the Reynolds number or the intrinsic slip
length is accompanied by the decrease or increase of the eﬀective slip length.

4.4

Conclusions

In this chapter the eﬀects of local slip boundary condition and the Reynolds number on
the ﬂow structure near sinusoidally corrugated surfaces and the eﬀective slip length were
investigated numerically by solving the Stokes and Navier-Stokes equations. The eﬀective
slip length was deﬁned with respect to the mean height of the surface roughness by extrapolating the linear part of the velocity proﬁle averaged over the corrugation period. In the
case of the Stokes ﬂow with the local no-slip boundary condition, the eﬀective slip length
decreases with increasing corrugation amplitude and the vortex ﬂow develops in the groove
of the rough surface for ka

0.79. In the presence of the local slip boundary condition

along the wavy wall, the eﬀective slip length increases and the size of the recirculation zone
is reduced. The vortex vanishes at suﬃciently large values of the intrinsic slip length. The
analysis of the pressure and wall shear stress computed from the Navier-Stokes equation
shows that the asymmetric vortex ﬂow develops in the groove due to the inertia term even
when the local slip boundary condition is applied. The eﬀective slip length decreases with
increasing Reynolds number. The numerical simulations suggest that the variation of the
vortex size as a function of either the Reynolds number or the intrinsic slip length correlates
with the magnitude of the eﬀective slip length.

67

U

h

u (z )

z
x

Leff

Solid wall

a

λ
Figure 4.1. Schematic diagram of the steady-state Couette ﬂow over a rough surface. The
upper ﬂat wall is moving with a constant velocity U in the x direction. The lower stationary
ˆ
wall is modeled as a sinusoidal wave with amplitude a and wavelength λ. The wavenumber
ka = 2 πa/λ varies in the range 0 ka 1.26.

68

(a)

(b)

2
1

η

2 (−1,1)

1(1,1)
ξ

3
4

z
x

3 (−1,−1)

4 (1,−1)

Transformation

Figure 4.2. Diagram of a bilinear element in (a) the physical coordinate system (x, z) and
(b) a transformed element in the natural coordinate system (ξ, η).

69

Panzer et al. Eq. (18)
Stokes solution

0.00

-0.10
-0.15

1.0
0.8
0.6
0.4

-0.20
-0.25
0.0

z/h

Leff / h

-0.05
ka = 0
.ka = 0.07
.ka = 1.12

0.2
*

u /U
0.0
-0.2 0.0 0.2 0.4 0.6 0.8 1.0
0.2

0.4

0.6

0.8

1.0

1.2

ka
Figure 4.3. The eﬀective slip length as a function of wavenumber ka computed from the
solution of the Stokes equation with the no-slip boundary condition. The primary vortex is
formed at the bottom of the valley for ka 0.79 (see the vertical arrow at ka = 0.79). The
solid line is calculated using Eq. (4.18). The inset shows the normalized velocity proﬁles
obtained from the Stokes solution for the selected values of ka. The dashed line located at
a = 0 is the reference for calculating the eﬀective slip length.

70

3.0

ka = 0.1
ka = 0.28
ka = 0.7
ka = 0.83
ka = 1.12

2.5

τw / τw

2.0
*

1.5
1.0
0.5
0.0
0.0

0.25

0.5

x/λ

0.75

1.0

Figure 4.4. Shear stress along the lower wavy wall computed from the Stokes solution with
no-slip boundary condition (L0 = 0) for the indicated values of wavenumber ka. The value
∗
τw used for normalization is the shear stress at the ﬂat wall. The intersection of the curves
with the dashed line determines the location of the ﬂow separation and attachment inside
the valley.

71

0.9
0.6

L0 / h = 0.0

P/ P

*

0.3
0.0
-0.3

ka = 0.1
ka = 0.28
ka = 0.7
ka = 0.83
ka = 1.12

-0.6
-0.9
-1.2
0.0

0.25

0.5

x/λ

0.75

1.0

Figure 4.5. The normalized pressure along the lower corrugated boundary extracted from
the Stokes solution with no-slip boundary condition for the same values of ka as in Fig. 4.4.
The maximum pressure P ∗ for ka = 1.12 and L0 = 0 is located in the bulk region at
(x/λ, z/λ) (0.07, 0.23).

72

1.0

( a)
-0.05

P/P *

z/h

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1

z /h

-0.1

0.5
-0.15
-0.002

0

0 u/U * 0.002

L 0 /h = 0
0

0.5

x/λ

1.0

Figure 4.6. Pressure contours and streamlines for the wavenumber ka = 1.12 and the slip
length at the lower wall L0 = 0 (a), L0 /h = 0.03 (b), L0 /h = 0.06 (c), and L0 /h = 0.08
(d). The pressure contours are normalized by the maximum value P ∗ located on the left
side of the peak (x/λ, z/λ) (0.07, 0.23) in the case of L0 = 0. The horizontal dashed line
denotes the location of the eﬀective no-slip boundary plane. The vertical dashed line inside
the valley at x/λ = 0.75 indicates the cross-section used to compute the velocity proﬁles
shown in the inset.

73

1.0

(b)
-0.05

z/ h

z /h

-0.1

0.5

-0.15
-0.002

0

0 u/U * 0.002

L 0 /h = 0.03
0

0.5

x/λ

Figure 4.6. continued.

74

1.0

1.0

(c)
-0.05

z/h

z /h

-0.1

0.5

-0.15
-0.002

0

0 u/U * 0.002

L 0 /h = 0.06
0

0.5

x/λ

Figure 4.6. continued.

75

1.0

1.0

(d)
-0.05

z/ h

z /h

-0.1

0.5

-0.15
-0.002 u/U * 0

0

0.002

L 0 /h = 0.08
0

0.5

x/λ

Figure 4.6. continued.

76

1.0

0.9
0.6

ka = 1.12

P/ P

*

0.3
0.0
-0.3

L0 / h = 0.0
L0 / h = 0.02
L0 / h = 0.04
L0 / h = 0.06
L0 / h = 0.08

-0.6
-0.9
-1.2
0.0

0.25

0.5

x/λ

0.75

1.0

Figure 4.7. Normalized pressure along the lower wavy wall extracted from the solution of
the Stokes equation for the tabulated values of L0 and ka = 1.12. The maximum pressure
P ∗ for L0 = 0 is located above the lower surface at (x/λ, z/λ) (0.07, 0.23).

77

3.0

L0 / h = 0.0
L0 / h = 0.02
L0 / h = 0.04
L0 / h = 0.06
L0 / h = 0.08

2.5

τw / τw

2.0
*

1.5
1.0
0.5
0.0
0.0

0.25

0.5

x/λ

0.75

1.0

Figure 4.8. Shear stress along the corrugated lower wall (ka = 1.12) for the indicated values
of the local slip length L0 . The intersection of the dashed line with the shear stress proﬁles
shows the location of the ﬂow separation and attachment inside the valley.

78

1.8
1.2

ka = 0.28
Panzer et al. Eq. (18)

Leff /h

0.6
0.0

0.1

0.1

1

10

100

ka = 1.12
Panzer et al. Eq. (18)

0.0
-0.1
0.01

0.1

L0 / h

1

10

Figure 4.9. The eﬀective slip length, extracted from the Stokes solution as a function of L0
for ka = 0.28 (top) and ka = 1.12 (bottom). The dashed line is computed from Eq. (4.18).
For ka = 1.12 and L0 /h 0.067, the vortex is formed in the valley (see the vertical arrow
at L0 /h = 0.067). The error bars are smaller than the symbol size.

79

0.12

L0 / h

0.09

(a)

no backflow

0.06
0.03

vortex flow

0.00

Leff / h

-0.07
-0.08

(b)

-0.09
0.8

0.9

1.0

ka

1.1

1.2

1.3

Figure 4.10. The intrinsic slip length L0 above which there is no vortex at the bottom of
the valley (top) and the corresponding eﬀective slip length Leﬀ (bottom) computed from
the Stokes solution and plotted as a function of wavenumber ka. The increment in the slip
length used to determine the threshold values is about the symbol size (top). The dashed
lines are quadratic ﬁts to guide the eye.

80

P/ P

*

30

(a)

L0 / h = 0.0

15
0
-15

τw / τw

80
*

(b)

Re = 39
Re = 79
Re = 118

60
40
20
0
0.0

0.25

0.5

x/λ

0.75

1.0

Figure 4.11. The normalized pressure (top) and shear stress (bottom) along the lower wavy
wall as a function of the Reynolds number for ka = 1.12 and no-slip boundary conditions.
The value P ∗ is the maximum pressure located above the lower boundary on the left side
of the peak (x/λ, z/λ) (0.07, 0.23) for L0 = 0 and Re = 4.0.

81

1.0

( a)

0

z/h

P/P*

z /h

-0.06

20
17
15
10
5
0
-5
-10
-15
-17
-20

-0.12

0.5

-0.006

0

-0.003

0 u/U *

Re= 79
L 0 /h = 0
0

0.5

x/λ

1.0

Figure 4.12. Pressure contours and streamlines extracted from the solution of the NavierStokes equation for ka = 1.12 and L0 = 0 (a) and L0 /h = 0.25 (b). The value P ∗ is the
maximum pressure located at (x/λ, z/λ) (0.07, 0.23) for ka = 1.12, L0 = 0, and Re = 4.0.
The horizontal dashed line denotes the eﬀective no-slip boundary plane. The vertical dashed
line inside the valley indicates the cross-section used to compute the velocity proﬁles shown
in the inset.

82

1.0
0

(b)

z/h

z /h

-0.06
-0.12

0.5

-0.006

0

-0.003

0 u/U *

Re= 85
L 0 /h = 0.25
0

0.5

x/λ

Figure 4.12. continued.

83

1.0

30

P/ P

*

(a)

L0 / h = 0.25

15
0
-15

τw / τw

20

(b)

Re = 43
Re = 85
Re = 127

*

10

0
0.0

0.25

0.5

x/λ

0.75

1.0

Figure 4.13. The normalized pressure (top) and shear stress (bottom) along the lower
wavy wall as a function of the Reynolds number for ka = 1.12 and L0 /h = 0.25. The
normalization value P ∗ is the maximum pressure located above the surface on the left side
of the peak (x/λ, z/λ) (0.07, 0.23) for L0 = 0 and Re = 4.0.

84

-0.02
-0.03

L0 / h = 0.0
L0 / h = 0.09
L0 / h = 0.12
L0 / h = 0.25

ka = 1.12

Leff / h

-0.04
-0.05
-0.06
-0.07
-0.08
-0.12
-0.13
0

20

40

60

80

100

120

Reynolds number

Figure 4.14. The eﬀective slip length computed from the solution of the Navier-Stokes
equation as a function of the Reynolds number for ka = 1.12 and L0 = 0 ( ), L0 /h = 0.09
( ), L0 /h = 0.12 ( ), and L0 /h = 0.25 ( ). The dashed line indicates the upper bound
of the region where a vortex is present in the groove of a rough surface. The error bars
associated with the threshold values of the Reynolds number are about the symbol size.

85

CHAPTER 5
Modeling the combined eﬀect of
surface roughness and shear rate on
slip ﬂow of simple ﬂuids
5.1

Introduction

An adequate description of ﬂuid ﬂow over a solid boundary with surface roughness on multiple length scales often requires resolution of ﬁne microscopic details of the ﬂow structure
while retaining peculiarities of a macroscopic picture [115]. The most popular and practically important example is the problem of transport through heterogeneous porous media,
which is an inherently multiscale system [116]. Modeling ﬂuid ﬂow at multiple scales is a
nontrivial problem itself; however, it becomes much more diﬃcult if one also needs to incorporate slip boundary conditions into equations of motion. Typically, a local slip appears
if the wall-ﬂuid interaction is suﬃciently low, but the value of slip length (an extrapolated
distance of the velocity proﬁle to zero) generally depends not only on the surface energy
and structure [5, 117, 10, 13, 14] but on shear rate as well [17, 10, 118, 19, 119]. Therefore,
a correct form of the Navier-Stokes equation can be obtained only if the combined eﬀect
of surface roughness and shear rate on the slip length is taken properly into account. It is
intuitively clear that the surface roughness and shear rate have an opposite impact on the
slip length. A precise evaluation of relative contributions of these two factors is the subject

86

of the present study.
The boundary conditions for the ﬂow of monatomic ﬂuids over atomically ﬂat surfaces can
be described in terms of the intrinsic slip length, which is deﬁned as the ratio of the shear
viscosity to the friction coeﬃcient at the liquid/solid interface. The friction is determined by
the strength of the interaction potential and molecular-scale corrugations of the crystalline
wall. Recent MD studies [120, 20, 21, 23, 25, 121, 28, 40, 29, 37] have reported values
of the intrinsic slip length up to about a hundred molecular diameters for weak wall-ﬂuid
interaction energies. It is well established that the formation of commensurate structures of
liquid and solid phases at the interface leads to stick boundary conditions [21, 28]. The slip
length decreases with increasing bulk pressure in the ﬂow of nonwetting simple ﬂuids over a
smooth substrate [25]. The MD simulations have also demonstrated that, depending on the
wall-ﬂuid interaction energy and the ratio of wall and ﬂuid densities, the slip length is either
relatively small (below a couple of molecular diameters) and shear-rate-independent [21] or
larger than several molecular diameters and increases with shear rate [23, 29]. Although
several functional forms for the rate-dependent intrinsic slip length have been suggested
for monatomic ﬂuids [23, 122, 29, 37], as well as for polymer melts [27, 32, 93, 123, 124],
the rate-dependent boundary conditions were not used for modeling slip ﬂows in complex
geometries.
In contrast to the description of the ﬂow over smooth boundaries (with microscopic surface
roughness) by the intrinsic slip length, it is more appropriate to characterize the ﬂow away
from macroscopically rough surfaces by the eﬀective slip length, which is usually deﬁned
with respect to the location of the mean roughness height. Recently, we investigated the
behavior of the eﬀective slip length in laminar ﬂow over periodically corrugated surfaces for
simple ﬂuids [35] and polymer melts [93]. Both MD and continuum simulations have shown
that the eﬀective slip length decreases with increasing amplitude of the surface roughness. In
the continuum analysis, the local (rate-independent) slip length is modiﬁed by the presence
of boundary curvature [51, 52] and, therefore, becomes position-dependent along the wavy
wall. A direct comparison between MD and continuum simulations at low shear rates
demonstrated that there is an excellent agreement between the eﬀective slip lengths for
large wavelengths (above approximately 30 molecular diameters) and small wavenumbers

87

ka

0.3 [35, 93]. In the present chapter, the analysis of the ﬂow over a periodically

corrugated surface is extended to higher shear rates and to the situation where the local
slip length at the curved boundary is a function of both the radius of curvature and shear
rate.
In the last decade, a number of studies have focused on the development of the
hybrid continuum-atomistic methods to simulate ﬂuid ﬂows near smooth solid boundaries [125, 126], rough walls [127], superhydrophobic surfaces [128], corner ﬂows [129], and
ﬂows near a moving contact line [130, 131]. In multiscale algorithms, the time and length
scales accessible in molecular simulations are considerably extended by employing the continuum method in the bulk region and applying the MD technique only near the interfaces.
The two methods are coupled via constrained dynamics in an overlap region by matching
mass, momentum and energy ﬂuxes or velocity ﬁelds. The full MD simulations are then
performed to validate the hybrid scheme. The drawback of the hybrid methods, however,
is the cumbersome procedure for coupling of the two methods in the overlap region. The
approach pursued in the present chapter is diﬀerent from that in the hybrid schemes. We
ﬁrst compute by MD simulations the intrinsic slip length as a function of shear rate in
the ﬂow over an atomically smooth surface. The Navier-Stokes equation is then solved numerically (in the whole computational domain) for the ﬂow over a rough surface with the
rate-dependent intrinsic slip length used as a local boundary condition. We found that the
full MD simulations recover the continuum results at low and intermediate shear rates and
small-scale surface roughness.
In this chapter, the dynamic behavior of the eﬀective slip length is investigated in steady
shear ﬂow over a periodically corrugated surface using molecular dynamics and continuum
simulations. Both methods show that the eﬀective slip length is nearly constant at low
shear rates and it gradually increases at higher shear rates. The slight discrepancy between
the eﬀective slip lengths obtained from MD and continuum simulations at high shear rates
is analyzed by examination of the local velocity and density proﬁles, as well as pressure and
temperature distributions along the wavy surface. It is also shown by MD simulations that
the singular behavior of the intrinsic slip length (in a ﬂow over a smooth wall) at high shear
rates is suppressed by the surface roughness.

88

The rest of this chapter is organized as follows. The details of molecular dynamics and
continuum simulations are described in the next section. The shear rate dependence of the
intrinsic and eﬀective slip lengths and the results of the comparative analysis are presented
in Sec. 5.3. The results are brieﬂy summarized in the last section.

5.2
5.2.1

The details of numerical methods
Molecular dynamics simulations

The system geometry and the deﬁnition of the eﬀective slip length are illustrated in Fig. 5.1.
The total number of ﬂuid monomers in the cell is ﬁxed to Nf = 16 104. The ﬂuid monomers
interact through the truncated pairwise Lennard-Jones (LJ) potential
VLJ (r) = 4ε

σ 12
σ 6
−
,
r
r

(5.1)

where ε is the energy scale and σ is the length scale of the ﬂuid phase. The cutoﬀ radius
for the LJ potential is set to rc = 2.5 σ for both ﬂuid/ﬂuid and wall/ﬂuid interactions. The
size of the wall atoms is the same as ﬂuid monomers. Wall atoms are ﬁxed at the lattice
sites and do not interact with each other.
The motion of the ﬂuid monomers was weakly coupled to a thermal reservoir via the
Langevin thermostat [54]. To avoid bias in the shear ﬂow direction, the random force and
friction terms were added to the equations of motion in the y direction, perpendicular to
ˆ
the plane of shear [21]. The equations of motion for ﬂuid monomers in all three directions
are as follows:
(5.2)

i=j

∂Vij
,
∂xi

(5.3)

i=j

∂Vij
+ fi ,
∂yi

(5.4)

i=j

∂Vij
,
∂zi

m¨i = −
x
m¨i + mΓyi = −
y
˙
m¨i = −
z

where Γ = 1.0 τ −1 is a friction coeﬃcient that regulates the rate of heat ﬂux between the ﬂuid
and the heat bath, and fi (t) is the random force with zero mean and variance 2mΓkB T δ(t)

89

determined from the ﬂuctuation-dissipation relation. Unless otherwise speciﬁed, the temperature of the Langevin thermostat is set to T = 1.1 ε/kB , where kB is the Boltzmann
constant. The equations of motion were integrated using the Verlet algorithm [56] with a
time step ∆t = 0.002 τ , where τ =

mσ 2 /ε is the characteristic time of the LJ potential.

The intrinsic slip length was computed at the lower ﬂat wall composed out of two layers
of the F CC lattice with density ρw = 3.1 σ −3 and nearest-neighbor distance d = 0.77 σ between atoms in the xy plane. The partial slip is permitted along the lower boundary due
to relatively low wall-ﬂuid interaction energy εwf = 0.5 ε. The no-slip boundary condition
is enforced at the ﬂat upper wall with a lower density 0.67 σ −3 by allowing the formation of
commensurate structures of liquid and solid phases at the interface. The nearest-neighbor
distance between wall atoms is 1.28 σ in the (111) plane of the F CC lattice. The wall-ﬂuid
interaction energy ε ensures the no-slip condition at the upper wall holds even at high shear
rates. The cell dimensions are set to Lx = 42.34 σ and Ly = 10.0 σ parallel to the conﬁning
walls. The MD simulations were performed at a constant ﬂuid density (ρ = 0.81 σ −3 ) ensemble, i.e., the relative distance between the walls was always ﬁxed to Lz = 41.41 σ in the
z direction. Periodic boundary conditions were applied along the x and y directions.
ˆ
ˆ
ˆ
The eﬀect of surface roughness on slip ﬂow was investigated in the cell with a lower
sinusoidal wall with wavelength λ = Lx and amplitude a = 1.82 σ (see Fig. 5.1). Special
care was taken to make the lower corrugated wall with a uniform density ρw = 3.1 σ −3 (by
including additional rows of atoms along the y direction) so that the lattice spacing along
ˆ
the sinusoidal curve remained d = 0.77 σ.
The ﬂuid monomer velocities were initialized using the Maxwell-Boltzmann probability
distribution at the temperature T = 1.1 ε/kB . After an equilibration period of about 2 ×
104 τ , the velocity and density proﬁles across the channel were averaged for about 2 × 105 τ
within bins of Lx × Ly × ∆z, where ∆z = 0.01 σ. The location and dimensions of the
averaging regions for the computation of the local velocity and density proﬁles will be
speciﬁed in Sec. 5.3.3.

90

5.2.2

Continuum simulations

The two-dimensional, steady-state and incompressible Navier-Stokes (NS) equation is solved
numerically using the ﬁnite element method. The equation of motion based on these assumptions is written as follows:
ρu ·

u = − p + µ 2 u,

(5.5)

ˆ
where u = u ˆ + v j is the velocity vector, ρ is the density of the ﬂuid, and p and µ are the
i
pressure ﬁeld and ﬂuid viscosity, respectively.
The penalty formulation was implemented for the incompressible solution to avoid decoupling of the velocity and pressure terms [57]. In the penalty formulation, the continuity
equation,

· u = 0, is replaced with a perturbed equation
p
·u=− ,
Λ

(5.6)

where Λ is the penalty parameter, which enforces the incompressibility condition. For most
practical situations, where computation is performed with double-precision 64 bit words,
a penalty parameter (Λ) between 107 and 109 is suﬃcient to conserve the accuracy [57].
After substitution of the pressure term in Eq. (5.6) into Eq. (5.5), the modiﬁed momentum
equation can be rewritten as follows:
ρu ·

· u) + µ 2 u.

u=Λ (

(5.7)

The Galerkin method with rectangular bilinear isoparametric elements is used to integrate
the Navier-Stokes equation [57, 132].
The boundary conditions are applied at the inlet, outlet, and upper and lower walls.
The periodic boundary conditions are implemented along the x direction for the inlet and
ˆ
outlet. The boundary condition at the upper wall is always no-slip. A partial slip boundary
condition is applied at the lower wall. In addition to the global coordinate system (ˆ, z ), a
x ˆ
local coordinate system (t, n) is deﬁned along the lower curved boundary with unit vectors
t and n representing the tangential and normal directions, respectively. The tangential
component of the ﬂuid velocity in the local coordinate system can be computed as
ut = L0 [(n ·

) ut + ut /R],
91

(5.8)

where R is the local radius of curvature and L0 is the intrinsic slip length at the ﬂat
liquid/solid interface [52]. The radius of curvature R is positive and negative for concave
and convex regions, respectively. The second-order forward-diﬀerencing scheme was used to
compute accurately the normal derivative of the tangential velocity (n ·

) ut at the lower

boundary.
At the beginning of the simulation procedure, the boundary conditions at the upper and
lower walls are set to no-slip. Then, the equations of motion are solved implicitly and the
tangential velocities at the lower boundary are updated according to Eq. (5.8). In the next
step, the updated velocities at the lower boundary are used as a new boundary condition and
the equations of motion are solved again. The iterative process is repeated until a desired
accuracy is achieved. The convergence rate of the solution was controlled by applying an
under-relaxation factor of 0.001 for the boundary nodes. For all results reported in the
present chapter, the continuum simulations were performed with 150 × 150 grid resolution.
However, the accuracy of the results was also checked by performing a set of simulations
using a ﬁner grid 180 × 180. The maximum error in the eﬀective slip length due to the grid
resolution is Leﬀ /Lz = 0.003. The same error bars were reported in the previous study of
laminar ﬂow over a periodically corrugated surface with either no-slip or (rate-independent)
slip boundary conditions [132].
The averaged error of the solution is deﬁned as
Np

error =
i=1

| un − un+1 |
i
i
| un+1 |
i

/Np ,

(5.9)

where Np is the number of nodes in the computational domain, un is the magnitude of
i
the velocity at node i and time step n, and un+1 is the magnitude of the velocity in
i
the next time step. The solution is converged when error
∂u

10−9 and the tangential

velocities along the lower boundary satisfy ut = Llocal ∂nt , where the local slip length is
Llocal = (L−1 − R−1 )−1 [52].
0

92

5.3
5.3.1

Results
MD simulations: the intrinsic and eﬀective slip lengths

We ﬁrst consider ﬂow in the cell with ﬂat upper and lower walls. The averaged velocity
proﬁles are plotted in Fig. 5.2 for the selected upper wall speeds. The proﬁles are linear
throughout the channel except for a slight curvature near the lower wall at U = 5.0 σ/τ .
Note that with increasing upper wall speed, the slip velocity at the lower wall increases
monotonically up to about 2 σ/τ at U = 5.0 σ/τ . The ﬂow velocity at the top boundary
remains equal to the upper wall speed due to the commensurability of liquid and solid
structures at the interface and suﬃciently high wall-ﬂuid interaction energy. The shear rate
was estimated from the linear slope of the velocity proﬁles across the entire cell.
The ﬂuid viscosity computed from the Kirkwood relation [58] is plotted in Fig. 5.3 as a
function of shear rate in the cell with ﬂat walls. In agreement with the results of previous MD
studies on shear ﬂow of simple ﬂuids [23, 29, 37], the shear viscosity µ = 2.15 ± 0.15 ετ σ −3
remains rate-independent up to γτ
˙

0.072, which corresponds to the upper wall speed

U = 5.375 σ/τ . It is also known that at high shear rates the ﬂuid temperature proﬁles become
nonuniform across the cell and the heating up is larger near the interfaces [29, 37, 32]. To
ensure that the ﬂuid viscosity is not aﬀected by the heating up, we performed two additional
sets of simulations at the upper wall speeds U = 1.0 σ/τ and U = 3.5 σ/τ and increased the
temperature of the Langevin thermostat with an increment of 0.05 ε/kB . In the range
of the thermostat temperatures and error bars reported in Fig. 5.3, the shear viscosity
is independent of temperature for both upper wall speeds. In Sec. 5.3.3, the results of
simulations performed at higher temperatures will be used to estimate the eﬀect of pressure
on the slip length.
The intrinsic slip length L0 in the ﬂow over a ﬂat lower wall was determined from the
linear extrapolation of velocity proﬁles to zero with respect to the reference plane located
0.5 σ above the F CC lattice plane. The variation of the intrinsic slip length as a function
of shear rate is presented in Fig. 5.4. At low shear rates γτ
˙

0.005, the slip length is nearly

constant. The leftmost data point reported in Fig. 5.4 is γτ
˙

1.7 × 10−4 at the upper wall

speed U = 0.01 σ/τ . At higher shear rates, the intrinsic slip length ﬁrst gradually increases
93

and then grows rapidly as it approaches γτ
˙

0.072 when U = 5.375 σ/τ . The behavior of

the slip length in the ﬂow over a ﬂat wall can not be well described by the power law function
suggested in Ref. [23]. Instead, a ninth-order polynomial function was used to ﬁt the data
(see the blue curve in Fig. 5.4). In the continuum analysis presented in the next section, the
polynomial function will be used as a local (rate-dependent) boundary condition, Eq. (5.8),
for the ﬂow over a periodically corrugated wall.
When the upper wall speed becomes greater than U = 5.375 σ/τ , the rate-dependence
of the intrinsic slip length passes through a turning point at the maximum shear rate
γτ
˙

0.072, and the slip length increases (L0

Lz ) while the shear rate decreases, leading

to a negative slope of the rate-dependent curve (not shown). In this regime the slip velocity
2
at the lower ﬂat wall is greater than the ﬂuid thermal velocity vT = kB T /m. The behavior

of the intrinsic slip length at very large slip velocities was not considered further in the
present study.
Next, the inﬂuence of surface roughness on the eﬀective slip length is investigated in
the ﬂow over the lower corrugated wall with wavenumber ka = 2πa/λ = 0.27. The velocity
proﬁles, averaged over the period of corrugation λ = Lx , are linear in the bulk of the ﬁlm and
curved near the lower wall within about 2 a from the bottom of the valley (e.g., see [93, 132]).
Therefore, the eﬀective slip length and shear rate were extracted from a linear ﬁt to the
velocity proﬁles inside the bulk region 10

z/σ

30. The eﬀective slip length as a

function of shear rate is also plotted in Fig. 5.4. Both the magnitude and slope of the ratedependence of the eﬀective slip length are reduced in comparison with the results for the ﬂat
wall. Following the rate-independent regime at γτ
˙

0.01, the eﬀective slip length increases

monotonically up to a maximum value Leﬀ (0.14 τ −1 )
At higher shear rates γτ
˙

16.6 σ at U = 8.0 σ/τ .

0.14, the eﬀective slip length starts to decrease with further

increase in shear rate (see Fig. 5.4). We found, however, that the ﬂuid temperature (density)
becomes anomalously high (low) on the right side of the corrugation peak. For example,
the ﬂuid temperature in the region 0.5

x/λ

0.9 varies up to T kB /ε

1.6 − 1.85 at

U = 12.0 σ/τ (the rightmost point in Fig. 5.4). Also, the amplitude of the characteristic
density oscillations normal to the curved boundary at x/λ

0.43 is reduced by about 40%

at U = 12.0 σ/τ in comparison with the ﬂuid layering at equilibrium conditions (not shown).

94

In the next section, the comparison between the eﬀective slip lengths estimated from MD
and continuum simulations is studied at shear rates γτ
˙

5.3.2

0.14.

Comparison between MD and continuum simulations

In this section, we describe the behavior of the rate-dependent eﬀective slip length computed
from the solution of the Navier-Stokes equation. All variables in continuum simulations are
normalized by the LJ length σ, time τ , and energy ε scales. The system is designed to mimic
the MD setup. All nondimensional parameters in continuum simulations are marked by the
˜
˜
(∼) sign. The cell dimensions are set to Lx = 42.34 and Lz = 40.41 in the x and z directions,
ˆ
ˆ
respectively. Similar to the MD setup, the wavenumber of the lower corrugated wall is ﬁxed
to ka = 0.27. This value was chosen based on the previous analysis of the eﬀective slip
length in laminar ﬂow of simple [35] and polymeric [93] ﬂuids over a rough surface. It was
shown [35] that at low shear rates there is an excellent agreement between the eﬀective
slip lengths extracted from MD and continuum simulations when a local (rate-independent)
intrinsic slip length is used as a boundary condition along the wavy wall with λ
ka

30 σ and

0.3.

It is important to emphasize that in continuum simulations the intrinsic slip length in
Eq. (5.8) is a function of the total shear rate at the curved boundary. The total shear rate is
computed as a sum of the normal derivative of the tangential velocity ∂ut /∂n and the ratio
of the slip velocity to the radius of curvature ut /R. The eﬀective slip length is extracted
from a linear part of the velocity proﬁles in the bulk region 10

z
˜

30.

Figure 5.4 shows the eﬀective slip length extracted from the continuum solution with the
local rate-dependent intrinsic slip length at the lower corrugated boundary. The upper wall
speed is varied in the range 0.025

˜
U

6. As expected, the eﬀective slip length obtained

from the continuum simulations agrees well with the MD results at low shear rates γτ
˙
where L0 is nearly constant. With increasing shear rate up to γτ
˙

0.01

0.1, the eﬀective slip

length increases monotonically. There is an excellent agreement between the two methods
at intermediate shear rates 0.01

γτ
˙

0.04. At higher shear rates, 0.04

γτ
˙

0.1, the

continuum predictions slightly overestimate the MD results. The maximum discrepancy in
the eﬀective slip length is about 2.5 σ at γτ
˙
95

0.1, which corresponds to the upper wall

˜
˜
speed U = 6.0. For U > 6.0, the local shear rate at the crest of the lower boundary,
[(n ·

) ut + ut /R], becomes larger than the highest shear rate, γτ
˙

0.072, achieved in the

MD simulations for the intrinsic slip length (see Fig. 5.4), and, therefore, the local boundary
condition, Eq. (5.8), cannot be applied.

5.3.3

A detailed analysis of the ﬂow near the curved boundary

In order to determine what factors cause the discrepancy between the eﬀective slip lengths
obtained from MD and continuum simulations at high shear rates, we performed a detailed
analysis of the ﬂuid temperature, pressure, density and velocity proﬁles near the lower
boundary. The ﬂuid temperature in the ﬁrst ﬂuid layer along the lower wall is plotted in
Fig. 5.5 for the indicated upper wall speeds. The temperature was computed by subtracting the mean ﬂow in the x and z directions and averaging the kinetic energy in all three
ˆ
ˆ
directions. As expected, at low shear rates γτ
˙

0.019 (U

1.0 σ/τ ), the ﬂuid temperature

remains equal to the value T = 1.1 ε/kB set in the Langevin thermostat. With increasing
upper wall speed, the temperature increases and eventually becomes nonuniformly distributed along the curved boundary (see Fig. 5.5). Somewhat surprisingly, we ﬁnd that at high
upper wall speeds U

5.0 σ/τ the ﬂuid temperature is lowest above the crest even though

the slip velocity in this region is expected to be higher than in the other locations.
The normal pressure distribution along the lower corrugated wall is shown in Fig. 5.6 for
the selected upper wall speeds. Each data point represents the ratio of the averaged normal
force on the wall atoms from the ﬂuid monomers (within the cutoﬀ radius) to the surface
area 0.77 σ × 10.0 σ in the ty plane. The small variation of the pressure at equilibrium is due
to the boundary curvature, i.e., the number of ﬂuid monomers within the cutoﬀ distance
from the wall atoms is on average slightly larger above the peak than inside the valley.
With increasing upper wall speed, the normal pressure increases in the valley (x/λ
and x/λ

0.12

0.54) while it becomes smaller than the equilibrium pressure on the right slope

of the peak (see Fig. 5.6). The normal pressure distribution extracted from the solution of
the NS equation (not shown) is qualitatively similar to the proﬁles presented in Fig. 5.6.
However, in the region x/λ

0.54 the MD data overestimate the continuum predictions

for the pressure distribution at U

3.0 σ/τ due to the increase in ﬂuid temperature (see
96

Fig. 5.5). It is clear from Fig. 5.6 that at higher upper wall speeds there is a considerable
deviation from the equilibrium pressure, which can aﬀect the local intrinsic slip length along
the curved boundary (see discussion below). We ﬁnally comment that the normal pressure
on the ﬂat lower wall (ka = 0) at equilibrium is equal to P
with the pressure on the locally ﬂat surface at x/λ
increases up to P

2.36 ε/σ 3 (which agrees well

0.5 for U = 0 in Fig. 5.6) and it

2.52 ε/σ 3 at the highest upper wall speed U = 5.375 σ/τ .

A snapshot of the ﬂow near the lower corrugated wall is shown in Fig. 5.7, where the
rectangular boxes illustrate the location of the averaging regions along the surface. The
dimensions of the boxes are d = 0.77 σ, 5.0 σ, and Ly = 10.0 σ along the curved boundary
(xz plane), normal to the curved boundary (xz plane), and in the y direction, respectively.
ˆ
The velocity and density proﬁles discussed below were averaged inside the boxes within bins
of thickness 0.01 σ in the direction normal to the surface.
In Figure 5.8 the averaged density proﬁles are plotted in regions I–IV for U = 0 and 6.0 σ/τ .
In both cases, a pronounced ﬂuid layering is developed near the lower wall. In the stationary
case, Fig. 5.8(a), the density proﬁles are almost the same in all four regions. Similar to the
results of previous MD studies on slip ﬂow over smooth surfaces [29, 37], when the upper
wall speed increases, the contact density (magnitude of the ﬁrst peak) is reduced due to a
ﬁnite slip velocity along the curved boundary [see Fig. 5.8(b)]. Among regions I–IV, it is
expected that the local slip velocity is highest above the crest (region II) and lowest at the
bottom of the valley (region IV), and, therefore, it is not surprising that the contact density
in the region IV is larger than in the region II [see Fig. 5.8(b)]. However, this correlation
does not hold in regions I and III, where it is expected that at ﬁnite Re the slip velocity in
the region I would be larger than in the region III; but, as shown in Fig. 5.8(b), the contact
density in the region I is slightly larger than in the region III. Instead, we correlate this
behavior with the nonuniform pressure distribution along the curved boundary where the
local normal pressure in the region I is higher than in the region III (see Fig. 5.6). Note also
that the location of the second and third ﬂuid layers in regions II and III is slightly shifted
away from the surface and the magnitude of the peaks is reduced due to the lower pressure
in the region 0.2

x/λ

0.5 at U = 6.0 σ/τ in Fig. 5.6.

The local tangential velocity proﬁles normal to the curved boundary are presented in

97

Fig. 5.9 for the upper wall speed U = 1.0 σ/τ . Note that the proﬁles are almost linear away
from the surface except for a slight curvature in the MD velocity proﬁles near the wall.
The slip velocity above the crest (region II) is larger than in the other regions and ut at
the bottom of the valley (region IV) is the lowest, as expected. The tangential velocity on
the left side of the peak (region I) is larger than on the right side (region III) due to the
inertial eﬀects [132]. The tangential velocity ut and its normal derivative ∂ut /∂n extracted
from MD simulations in regions I, III, and IV agree rather well with the continuum results.
Despite a slight overestimation of the continuum velocity with respect to the MD results
above the crest (region II), the eﬀective slip length is essentially the same in both methods
for the upper wall speed U = 1.0 σ/τ (see the data at γτ
˙

0.019 in Fig. 5.4).

Figure 5.10 shows local proﬁles of the velocity component normal to the curved boundary
for the upper wall speed U = 1.0 σ/τ . Within statistical uncertainty, the normal velocities
at the boundary are zero in both MD and continuum simulations. In regions I and IV the
normal velocity is negative, while in regions II and III the velocity is positive away from the
boundary. The negative and positive signs of the normal velocity are expected in regions I
and III because of the positive and negative slopes of the boundary with respect to the shear
ﬂow direction. Although the boundary slope at the crest and at the bottom of the valley is
parallel to the upper wall, the local symmetry of the ﬂow with respect to x/λ = 0.25 and
0.75 is broken at ﬁnite Reynolds numbers (Re ≈ 15.6 for U = 1.0 σ/τ in Fig. 5.10), and, as
a result, the normal velocity in regions II and IV is not zero. Furthermore, the MD velocity
proﬁles exhibit an unexpected oscillatory pattern away from the surface. These oscillations
correlate well with the locations of the minima among the ﬁrst three ﬂuid layers [e.g., see
Fig. 5.8(a)]. Visual inspection of the individual molecule trajectories indicates that the ﬂuid
monomers mostly undergo thermal motion within the moving layers but occasionally jump
back and forward between the layers. The net ﬂux between the layers, however, is not zero
and its sign is determined by the direction of the ﬂow near the boundary.
The tangential velocity proﬁles averaged in regions I–IV at a higher upper wall speed
U = 6.0 σ/τ are shown in Fig. 5.11. The ﬂow pattern near the lower boundary is similar
to that described in Fig. 5.9 for U = 1.0 σ/τ , but the relative diﬀerence between the MD
and continuum velocities is larger in each region which results in a maximum discrepancy

98

between the eﬀective slip lengths of about 2.5 σ at γτ
˙

0.1 in Fig. 5.4. We tried to

estimate the local intrinsic slip length from the MD velocity proﬁles in regions I–IV and
to compare L0 with the continuum predictions. However, the normal derivative of the
tangential velocity cannot be computed accurately because of the curvature in the velocity
proﬁles and ambiguity in determining the range for the best linear ﬁt. The resulting error
bars associated with the linear extrapolation of the MD velocity proﬁles in Fig. 5.11 are
greater than the maximum diﬀerence between MD and continuum results for Leﬀ at high
shear rates in Fig. 5.4. Instead, a more reliable data for the local intrinsic slip length (see
Fig. 5.12) are obtained by computing the friction coeﬃcient (kf = µ/L0 ), which is the ratio
of the wall shear stress to the slip velocity. The agreement between the MD and continuum
results is quite good on the right side of the peak in Fig. 5.12. The intrinsic slip length
computed from MD simulations is about 3 − 4 σ smaller than the continuum L0 on the left
side of the peak due to higher ﬂuid pressure in that region (see also discussion below).
The normal velocity proﬁles averaged in regions I–IV along the lower boundary are displayed in Fig. 5.13 for the upper wall speed U = 6.0 σ/τ . In each region, the amplitude of
oscillations in the MD velocity proﬁles is increased with respect to the mean ﬂow predicted
by the continuum analysis. Similar to the correlation between density and velocity proﬁles
observed in Fig. 5.10 for U = 1.0 σ/τ , the locations of the maxima in the normal velocity
proﬁles in Fig. 5.13 correspond to the minima in the density proﬁles shown in Fig. 5.8(b).
At the higher upper wall speed, U = 6.0 σ/τ , however, the location of the ﬁrst peak in the
velocity proﬁle above the crest (region II) and on the right side of the peak (region III) is
slightly displaced away from the boundary due to a small shift in the position of the second
and third ﬂuid layers in regions II and III [see Fig. 5.8(b)].
The eﬀect of pressure variation along the curved boundary is not included in the ratedependent intrinsic slip length, Eq. (5.8), used in continuum simulations. It was shown in
Fig. 5.6 that at higher upper wall speeds, U

3.0 σ/τ , the ﬂuid pressure along the lower

corrugated wall deviates signiﬁcantly from the equilibrium pressure. Next, we estimate the
inﬂuence of pressure on the intrinsic slip length in shear ﬂow over a ﬂat wall by increasing the
temperature of the Langevin thermostat. The MD simulations described in Sec. 5.3.1 were
repeated at the upper wall speeds U = 1.0 σ/τ and 3.5 σ/τ to examine L0 at low (γτ ≈ 0.018)
˙
99

and high (γτ ≈ 0.055) shear rates. The results are presented in Fig. 5.14 for the indicated
˙
temperatures of the Langevin thermostat. As the bulk pressure increases, the intrinsic slip
length is reduced by about 3 σ at U = 3.5 σ/τ and about σ at U = 1.0 σ/τ . Note, however,
that under these ﬂow conditions the shear viscosity remains nearly constant (see Fig. 5.3),
and the variation of the friction coeﬃcient (kf = µ/L0 ) is determined only by the intrinsic
slip length.
The maximum increase in pressure and the diﬀerence in L0 reported in Fig. 5.14(a) correlate well with the variation of pressure in Fig. 5.6 and the intrinsic slip length in Fig. 5.12
on the left side of the peak at U = 6.0 σ/τ . We conclude, therefore, that the discrepancy
between the eﬀective slip lengths computed from MD and continuum simulations at high
shear rates (shown in Fig. 5.4) is caused by the increase in ﬂuid pressure on the left side of
the peak (where the solid boundary faces the mainstream ﬂow) and, as a result, the local
intrinsic slip length is reduced. Thus, a more accurate continuum description of the slip ﬂow
over a curved boundary can be achieved by including the eﬀect of pressure and temperature
on the rate-dependent intrinsic slip length.

5.4

Conclusions

In this chapter, molecular dynamics simulations were used to determine the intrinsic slip
length as a function of shear rate in steady ﬂow over an atomically ﬂat surface. It was found,
in agreement with previous MD studies, that the intrinsic slip length is nearly constant (rateindependent) at low shear rates and it grows rapidly as the shear rate approaches a critical
value. In the presence of periodic surface roughness, the magnitude of the eﬀective slip
length is signiﬁcantly reduced and its dependence on shear rate becomes less pronounced.
For the same surface roughness, the Navier-Stokes equation was solved numerically with
the rate-dependent boundary conditions speciﬁed by the intrinsic slip length obtained from
MD simulations. The continuum simulations reproduce accurately the MD results for the
eﬀective slip length at low and intermediate shear rates. The small diﬀerence between the
eﬀective slip lengths at high shear rates was analyzed by performing a detailed analysis of
the local velocity ﬁelds and pressure distribution along the curved boundary. We found

100

that the main cause of the discrepancy between MD and continuum results at high shear
rates is the reduction of the local intrinsic slip length in the region of higher pressure where
the boundary slope becomes relatively large with respect to the mainstream ﬂow. These
ﬁndings suggest that the continuum analysis of the ﬂow over a rough surface at high shear
rates should take into account the eﬀect of pressure on the rate-dependent intrinsic slip
length. One way to include pressure in the numerical analysis is to perform an additional
set of MD simulations (for ﬂat walls) at diﬀerent values of the constant normal load and
vary the upper wall speed. Thus, the intrinsic slip length as a function pressure and shear
rate can be obtained from interpolation of the MD data and then used as a local boundary
condition in the continuum simulation.

101

U

Solid Wall

ux( z )

Leff

z
y

Lz

a
x

λ

Figure 5.1. The projection of x and z coordinates of the ﬂuid monomers (open circles)
conﬁned between solid walls (ﬁlled circles). The atoms of the lower stationary wall are
uniformly distributed along the sinusoidal curve with wavelength λ, amplitude a, and density
ρw = 3.1 σ −3 . The shear ﬂow is induced by the ﬂat upper wall with the density 0.67 σ −3
moving with a constant velocity U in the x direction. The eﬀective slip length Leﬀ is
ˆ
computed by extrapolation of the velocity proﬁle to ux = 0.

102

7
6

< ux> τ / σ

5
4

U = 0.1 σ/τ
U = 0.5 σ/τ
U = 1.0 σ/τ
U = 3.0 σ/τ
U = 5.0 σ/τ

3
2
1
0
10

20

z /σ

30

40

Figure 5.2. Averaged velocity proﬁles for the tabulated upper wall speeds in the cell with
ﬂat upper and lower walls. The solid lines are the best linear ﬁts to the data. The vertical
axes indicate the location of the F CC lattice planes.

103

3.2

µ / ε τ σ −3

2.8

T = 1.1 ε / kB, 0.01 ≤ U τ / σ ≤ 5.375
1.1 ≤ T kB/ε ≤ 1.4, U = 1.0 σ / τ
1.1 ≤ T kB/ε ≤ 1.35, U = 3.5 σ / τ

2.4
2.0
1.6

0.00 0.01 0.02 0.03 .0.04 0.05 0.06 0.07

γτ

Figure 5.3. Shear viscosity in the cell with ﬂat upper and lower walls when the temperature
of the Langevin thermostat is set to T = 1.1 ε/kB ( ). Fluid viscosity as a function of the
thermostat temperature 1.1 T kB /ε 1.4 for the upper wall speed U = 1.0 σ/τ ( ) and
1.1 T kB /ε 1.35 for U = 3.5 σ/τ (♦).

104

Leff /σ

30

MD lower flat wall
MD rough wall ka = 0.27
Navier-Stokes ka = 0.27

20

10
0.001

.0.01
γτ

0.1

Figure 5.4. The intrinsic slip length L0 /σ as a function of shear rate for atomically ﬂat
walls ( ). The blue solid curve is a ninth-order polynomial ﬁt to the data. The eﬀective
slip length Leﬀ /σ as a function of shear rate for the ﬂow over the corrugated wall with
wavenumber ka = 0.27 (◦). The eﬀective slip length computed from the solution of the
Navier-Stokes equation with the local rate-dependent slip length (•).

105

Stationary
U = 1.0 σ/τ
U = 3.0 σ/τ
U = 5.0 σ/τ
U = 6.0 σ/τ

1.5

T kB / ε

1.4
1.3

ka = 0.27

1.2
1.1
1.0
0.00

0.25

0.50

x/λ

0.75

1.00

Figure 5.5. Temperature of the ﬁrst ﬂuid layer T (x) (in units of ε/kB ) along the lower
corrugated wall (ka = 0.27) for the indicated upper wall speeds. The vertical dashed lines
denote the location of the crest (x/λ = 0.25) and the bottom of the valley (x/λ = 0.75).

106

4.5
Stationary
U = 1.0 σ/τ
U = 3.0 σ/τ
U = 5.0 σ/τ
U = 6.0 σ/τ

4.0

P n σ 3/ ε

3.5
3.0
2.5
2.0

ka = 0.27

1.5
0.00

0.25

0.50

x/λ

0.75

1.00

Figure 5.6. The normal pressure Pn (in units of ε/σ 3 ) along the lower corrugated wall
(ka = 0.27) for the selected upper wall speeds. The vertical dashed lines indicate the
location of the crest (x/λ = 0.25) and the bottom of the valley (x/λ = 0.75).

107

z
y

x
I

II
III

IV

n
t

0

0.25

0.5

0.75

1.0

x/λ
Figure 5.7. A snapshot of the ﬂuid monomers (open circles) near the lower corrugated wall
(ﬁlled circles) with wavenumber ka = 0.27. The rectangular boxes denote averaging regions
on the left side of the peak (I) (x/λ 0.05), above the crest (II) (x/λ 0.25), on the right
side of the peak (III) (x/λ 0.45), and at the bottom of the valley (IV) (x/λ 0.75). The
unit vectors t and n indicate the tangential and normal directions at the boundary.

108

(a)

3

Region I
Region II
Region III
Region IV

U=0

2

ρσ3

1
0

1

3

2
(b)

3

4

3

4

U = 6.0 σ/τ

2
1
0

1

2

n /σ

Figure 5.8. Fluid density proﬁles averaged inside regions I–IV (shown in Fig. 5.7) near the
corrugated lower wall with wavenumber ka = 0.27. The upper wall speeds are U = 0 (a)
and U = 6.0 σ/τ (b).

109

2.5
Region I

n /σ

2.0

Region II
Region III

1.5

Region IV

1.0

0.5

MD
NS
MD
NS
MD
NS
MD
NS

U = 1.0 σ/τ

0.1

0.2

< ut> τ / σ

0.3

Figure 5.9. Averaged tangential velocity proﬁles inside regions I–IV obtained from MD
simulations (open symbols) and the solution of the Navier-Stokes equation (ﬁlled symbols).
The location of the averaging regions near the lower corrugated wall is shown in Fig. 5.7.
The upper wall speed is U = 1.0 σ/τ .

110

2.5
Region I

n /σ

2.0

Region II
Region III

1.5

1.0

0.5
-0.06

Region IV

MD
NS
MD
NS
MD
NS
MD
NS

U = 1.0 σ/τ
-0.04

-0.02

< un > τ / σ

0.00

0.02

Figure 5.10. Averaged normal velocity proﬁles inside regions I–IV (shown in Fig. 5.7) for
the upper wall speed U = 1.0 σ/τ . The proﬁles are extracted from MD simulations (open
symbols) and the solution of the Navier-Stokes equation (ﬁlled symbols). The horizontal
dashed lines indicate the location of the minima in density proﬁles (see text for details).

111

2.5
Region I

n /σ

2.0

Region II
Region III

1.5
Region IV

1.0

0.5
0.0

MD
NS
MD
NS
MD
NS
MD
NS

U = 6.0 σ/τ

0.5

1.0

1.5

< ut> τ / σ

2.0

2.5

Figure 5.11. Tangential velocity proﬁles averaged inside regions I–IV (shown in Fig. 5.7)
near the lower corrugated wall. The upper wall speed is U = 6.0 σ/τ . The velocity proﬁles
are extracted from MD simulations (open symbols) and the solution of the Navier-Stokes
equation (ﬁlled symbols).

112

34
32

MD rough wall ka = 0.27
Navier-Stokes ka = 0.27

L 0 /σ

30
28
26
24
22

U =6.0 σ/τ

20
0.00

0.25

0.50
x/λ

0.75

1.00

Figure 5.12. The intrinsic slip length along the lower corrugated wall (ka = 0.27) for the
upper wall speed U = 6.0 σ/τ . The data are obtained from MD simulations (open symbols)
and the solution of the Navier-Stokes equation (ﬁlled symbols). The vertical dashed lines
denote the position of the crest (x/λ = 0.25) and the bottom of the valley (x/λ = 0.75).

113

2.5
Region I

n /σ

2.0

Region II
Region III

1.5
Region IV

MD
NS
MD
NS
MD
NS
MD
NS

1.0

U = 6.0 σ/τ
0.5
-0.4

-0.3

-0.2

-0.1

< un> τ / σ

0.0

0.1

Figure 5.13. Averaged normal velocity proﬁles inside regions I–IV obtained from MD simulations (open symbols) and the solution of the Navier-Stokes equation (ﬁlled symbols). The
upper wall speed is U = 6.0 σ/τ . The horizontal arrows indicate the location of the minima
in density proﬁles between the ﬁrst and second ﬂuid layers in regions I and III.

114

L 0 /σ

24

U = 3.5 σ/τ

22
20
18

(a)

L 0 /σ

17

U =1.0 σ/τ

16
15
14

(b)

2.5

3.0

P σ 3/ ε

3.5

Figure 5.14. The intrinsic slip length L0 /σ as a function of the bulk pressure P (in units
of ε/σ 3 ) in the cell with ﬂat upper and lower walls. The temperature of the Langevin
thermostat is varied in the range 1.1 T kB /ε 1.35 for the upper wall speed U = 3.5 σ/τ
(a) and 1.1 T kB /ε 1.4 for U = 1.0 σ/τ (b).

115

CHAPTER 6
Slip boundary conditions for the
moving contact line in molecular
dynamics and continuum simulations
6.1

Introduction

The problem of contact line motion, which involves ﬂow of two immiscible viscous ﬂuids
in contact with a solid substrate, has captivated a lot of attention because of its fundamental diﬃculties and important applications including droplet spreading, spray painting,
lubrication, and coating [134, 135].
In wetting problems, the line of contact of the ﬂuid/ﬂuid interface with the solid wall
is called the contact line. In the absence of the ﬂuid ﬂow relative to the solid wall, the
ﬂuid/ﬂuid interface forms an angle with the solid wall, called the static contact angle θs ,
which is determined by the balance of the surface tensions between three media. In the
presence of ﬂow, the shape of the ﬂuid/ﬂuid interface near the wall becomes deformed and
it is described by the dynamic contact angle θd . The deformation of the interface depends on
the ﬂow ﬁelds near the contact line, surface tensions, and energy dissipation at the contact
line.
The analytical calculation of ﬂow ﬁelds and motion of the contact line with no-slip boundary condition along the solid surface leads to a stress singularity and inﬁnite viscous dissipa-

116

tion near the contact line [136, 137]. Due to multi-scale nature of the moving contact line
problem and breakdown of the local hydrodynamics near the contact point, the problem has
inspired many theoretical/numerical studies over the last 30 years [62, 137, 138, 141, 144].
On the other hand, experimental studies are limited to measuring a macroscopic contact
angle, known also as an apparent contact angle θa , which is accessible at a distance R (order
of microns) from the surface. When the shape of the ﬂuid/ﬂuid interface is deformed on a
molecular distance from the wall, the static contact angle θs and dynamic contact angle θd
are distinguishable while θa depends on the value of R.
Most of analytical and numerical solutions of the moving contact line problem were performed based on the fact that the far ﬁeld solution of the problem is almost insensitive
to the exact dynamics of the contact line close to wall [138]. One of the possibilities to
remove the singularity at the contact line is to replace the no-slip boundary condition with
slip BC via the Navier model. Some of the challenges in applying the slip models near the
ﬂuid/solid interface include a correct speciﬁcation of the region near the contact line where
the boundary slip must be applied, the magnitude of the slip length, and the shape of the
ﬂuid/ﬂuid interface.
The molecular dynamics (MD) simulations have been only recently used to study the
moving contact line problem [20, 130, 131, 139, 140, 141, 142, 143, 144, 145]. One of the
advantages of the MD method is that it oﬀers a detailed information about microscopic
ﬂow structure near the contact line without any assumptions about boundary conditions or
constitutive relations. Thus, an accurate physical description of the ﬂow ﬁelds and motion of
the contact line can be obtained from MD simulations, which in turn might help to develop
models for characterizing the ﬂow behavior in continuum modeling.
In this chapter, the moving contact line problem is studied using a combined approach of
molecular dynamics and continuum simulations. First, MD simulations are carried out to
resolve the ﬂow ﬁelds in the vicinity of the moving contact line on the molecular scale and
to establish a relation between the interfacial shear stresses and slip length as a function of
the distance from the contact line. This is done by extracting the velocity proﬁles and the
shear and normal stress tensors at the ﬂuid/solid interface. Second, in the presence of shear
ﬂow, the shape of the ﬂuid/ﬂuid interface is accurately resolved using MD simulations. In

117

the next step, the boundary conditions near the moving contact line extracted from MD
simulations are implemented as input conditions in the continuum solution of the NavierStokes equation (using the ﬁnite diﬀerence/particle tracking methods) to reproduce velocity
proﬁles and the shape of the ﬂuid/ﬂuid interface.
The rest of this chapter is organized as follows. The details of MD simulations, equilibration procedure, careful evaluation of normal stresses near the solid substrate, averaging
velocity proﬁles are presented in the Sec. 6.2. The continuum model and comparative analysis are described in Sec. 6.3. A brief summary is given in the last section.

6.2
6.2.1

Molecular dynamics (MD) simulations
Equilibrium case

The system consists of two immiscible ﬂuids conﬁned between ﬂat walls. The initial conﬁguration of the atoms is presented in Fig. 6.3. Each of two solid walls is modeled as two layers
of F CC lattice in [111] direction. The number of ﬂuid monomers in ﬂuid 1 and 2 are ﬁxed
to Nf1 = 18900 and Nf2 = 18900, respectively. The total number of wall atoms is set to
Nw = 10600. The MD simulations were performed at a constant ﬂuid density (ρ = 0.81 σ −3 )
ensemble, while the wall density is set to ρw = 1.86 σ −3 to allow a partial slip ﬂow along the
walls in each single ﬂuid phase. The thermostat temperature is set to T = 1.1 ε/kB , where
kB is the Boltzmann constant. The system size in the x and y directions is Lx = 241.87 σ
ˆ
ˆ
and Ly = 7.90 σ. The relative distance between the walls is ﬁxed to Lz = 25.41 σ in the z
ˆ
direction. Periodic boundary conditions were applied along the x and y directions.
ˆ
ˆ
The ﬂuid monomers interact through the truncated pairwise Lennard-Jones (LJ) potential
VLJ (r) = 4ε

σ 12
σ 6
−δ
,
r
r

(6.1)

where ε is the energy scale and σ is the length scale of the ﬂuid phase. The cutoﬀ radius for
the LJ potential is set to rc = 2.5 σ for the ﬂuid/ﬂuid and wall/ﬂuid interactions. For the
contact line problem, the LJ potential is modiﬁed by adding a parameter δ in the attractive
part of the potential. For ﬂuid 1/ﬂuid 1, ﬂuid 2/ﬂuid 2, and wall/ﬂuid interactions δ = 1.
However, for immiscible ﬂuids (ﬂuid 1/ﬂuid 2), δ = −1 and the LJ potential only consists
118

of the repulsive term. In Fig. 6.4 the LJ potential is plotted as a function of the distance
between two monomers for diﬀerent values of δ. It is clear from Fig. 6.4 that for all values
of r/σ when δ = −1 the potential is repulsive. Throughout the simulations the interaction
energies between the ﬂuids and walls are the same (εf1w = εf2w ).
For a system with stationary walls, after an equilibration period of about 2 × 104 τ ,
the shapshot of ﬂuid molecules of both phases and the shape of the ﬂuid/ﬂuid interface
are presented in Fig. 6.5. The ﬂuid/ﬂuid interface is almost a straight line and the static
contact angle is θs ≈ 90. This value is in agreement with the angle predicted by the Young’s
equation γf1w + γcosθs = γf2w , when γf1w = γf2w .
The ﬂuid properties in equilibrium case are studied by computing the density and the
normal stress distributions inside the ﬂuid domain. First, the ﬂuid density is averaged in
bins measuring 50 σ × 7.90 σ × 0.01 σ in the x, y , and z directions in ﬂuid 1 away from the
ˆ ˆ
ˆ
contact line, respectively. The density proﬁle along the z direction is shown in Fig. 6.6(a).
ˆ
In Fig. 6.6(b) the ﬂuid density is averaged in bins inside the bulk region with dimensions
0.1 σ × 7.90 σ × 5 σ in the x, y , and z directions, respectively. The density proﬁles are
ˆ ˆ
ˆ
plotted along the channel in the x direction. In both cases the density in the bulk region is
ˆ
0.81 σ −3 . In Fig. 6.6(a) pronounced oscillations in the density proﬁle are due to the layering
of the ﬂuid molecules near the walls. The ﬂuid densities are reduced in Fig. 6.6(b) for both
ﬂuids as they get closer to the ﬂuid/ﬂuid interface region. The location of the contact line
is estimated from the intersection of two density proﬁles. The width of the interface is
approximately 5 σ.
In addition to the ﬂuid densities, calculating the normal stresses provides an insight into
the ﬂuid characteristic near the ﬂuid/ﬂuid or ﬂuid/solid interfaces. Estimating the stresses
correctly is a diﬃcult task even in the equilibrium case. Equilibrium simulations can provide
an opportunity to check whether our numerical algorithm for the stress calculation is correct
and, as a result, the surface tension is computed accurately. As discussed in previous
chapters, the Irving-Kirkwood relation is used to compute normal and shear stresses [58].
However, in our previous studies the stresses were only estimated in the bulk region to
compute the ﬂuid viscosity [32, 93, 146]. In this chapter, the algorithm for computing
stresses near the interfaces is modiﬁed. If we assume that any two monomers are connected

119

through an imaginary line, while the stress is calculated from the combination of kinetic and
potential energies between the two monomers, the stress associated with each pair should be
partitioned in proportion to the fraction of the line segment between them that lies in each
spatial bin [147]. In this study the imaginary line between any two monomers is divided
into 100 segments regardless of the distance between the monomers.
The results for τxx , τyy , and τzz normal stresses in the x, y , and z directions, are plotted
ˆ ˆ
ˆ
in Fig. 6.7. The legends show the values of normal stresses in units of ε/σ 3 . As it is
required by the conservation of the momentum, the component of pressure tensor normal to
any interface is constant across that interface. Therefore, the stress τxx is constant across
the ﬂuid/ﬂuid interface (ˆ direction) and τzz is constant across the ﬂuid/wall interface (ˆ
x
z
direction). This also indicates that the modiﬁed stress algorithm is accurate and it can be
used for further analysis, in particular, the evaluation of stresses in a shear ﬂow. The next
step is to calculate the surface tension of the ﬂuid/ﬂuid interface
γ=

(τ⊥ − τ ) dx,

(6.2)

where τ⊥ and τ are the normal stresses perpendicular and parallel to the interface, respectively. The value of the surface tension calculated from Eq. (6.2) is about γ = 3.7 ± 0.2εσ −2 ,
which agrees well with previous studies for similar ﬂow conditions [141, 144].
Fig. 6.8(a) shows an enlarged view of an equilibrated system near the contact line. The
depletion region of ﬂuid monomers near the ﬂuid/ﬂuid interface approximately deﬁnes the
location of the contact line. This is also consistent with the density reduction in each ﬂuid
phase close to the interface in Fig. 6.6(b). The arrows indicate the direction of the tangential
forces Fx from the solid wall applied on each ﬂuid phase at the ﬂuid/solid interface near the
contact line. The distribution of the tangential force per unit area is deﬁned as the shear
stress τxz (x). In Fig. 6.8(b), far away from the contact line region, the shear stress vanishes
as expected in the equilibrium case. However, the net tangential force on the left and right
sides of the contact line becomes negative and positive, respectively. The ﬁnite value of the
tangential force is due to the repulsive interaction between ﬂuid monomers that belong to
diﬀerent ﬂuid phases.

120

6.2.2

Shear ﬂow

In this section, the eﬀects of shear ﬂow on the moving contact line, the contact angle, and
velocity proﬁles are presented. The shear ﬂow is induced by translating the upper and lower
walls in opposite directions while the center of mass of the ﬂuid phases remains in the center
of the channel. The conﬁguration of monomers for the selected upper/lower wall speeds is
plotted in Fig. 6.9. As the wall velocity U increases, the ﬂuid/ﬂuid interface remains straight
and the contact angle, measured in the advancing ﬂuid, becomes larger. The values of the
contact angle as a function of the ﬂuid velocity are shown in Fig. 6.10. For the velocities
considered in the present study, the contact angle increases linearly from θs = 90◦ when
U = 0 to θd = 115◦ when U = 0.1 σ/τ . However, at higher velocities, U

0.12 σ/τ , the

ﬂuid/ﬂuid interface undergoes an unsteady motion and breaks up, while the contact angle
varies signiﬁcantly (not shown).
For U = 0.1 σ/τ , the conﬁguration of ﬂuid monomers near the contact line region is
shown in Fig. 6.11 (top). The ﬂow behavior in this system is studied by averaging the
instantaneous velocities of ﬂuid monomers in small spatial bins over a long period of time.
The average velocity vectors are estimated inside bins measuring 1.0 σ × 7.90 σ × 0.9 σ in
the x, y , and z directions, respectively. An example of the averaged velocity vectors are
ˆ ˆ
ˆ
plotted in Fig. 6.11 (bottom). The direction of vectors indicates the presence of circulation
zones, which are centered in the bulk region of each ﬂuid phase while rotating in a clockwise
direction. The ﬂuid ﬂow cannot penetrate the ﬂuid/ﬂuid interface, and therefore, it is forced
to adjust its direction and move parallel to the interface.
Inspired by the detailed resolution of the velocity proﬁles inside the ﬂuid domain, the
slip velocities can also be calculated near the ﬂuid/solid interface. Fig. 6.12 shows the
slip velocity of the ﬁrst ﬂuid layer averaged separately in each ﬂuid phase. Due to weak
wall/ﬂuid interaction energy and incommensurable wall/ﬂuid densities, a ﬁnite value of the
slip velocity is observed far away from the contact line. In each ﬂuid phase (indicated by
symbols in Fig. 6.12), the slip velocity increases near the contact line and approaches the
wall speed. However, the the combined eﬀect of slip ﬂows from two ﬂuids (shown by the
dashed line) results in the peak of the slip velocity which is smaller than the wall speed. A

121

slightly smaller value of the peak in the slip velocity can be explained by the fact that the
contact line region is not a line with an inﬁnitesimal thickness but an area expanding about
5 σ in each ﬂuid (see the density proﬁles in Fig. 6.6(b) and stress distributions in Fig. 6.8(b)).
Moreover, the diﬀusive motion of the contact line can be another factor attributing to the
reduced slip velocity at the contact line.
In Fig. 6.13(a) the slip velocity at the solid/liquid interface is plotted for the selected values
of the walls speed U = 0.01σ/τ , U = 0.05σ/τ , and U = 0.1σ/τ . The slip velocity increases
with the wall speed. An important question that needs to be answered for the continuum
modeling is what should be the appropriate boundary condition near the contact line and
whether the Navier slip law is valid in the region near the contact line. To investigate this
issue, the slip length and the friction coeﬃcient away from the contact line in a single ﬂuid
phase must be calculated. In Fig. 6.13(b), the velocity proﬁles in the bulk region away from
the ﬂuid/ﬂuid interface are shown as a function of the wall speed. The velocity proﬁles
are linear, and, therefore, the slip length is extracted from the extrapolation of a linear ﬁt
to the proﬁles. The intrinsic slip length for the ﬂow of a single phase is estimated to be
L0 = 2.2 ± 0.5 σ. The ﬂuid viscosity was calculated from the ratio of the bulk shear stress
and shear rate away from the contact line τxz /γ = µ = 2.0 ± 0.2ετ σ −3 . As a result, the
˙
friction coeﬃcient at the interface between a single phase ﬂuid and the solid wall is the ratio
of the shear viscosity to the slip length, i.e., βslip = µ/L0 = 0.9. These information will be
used in the next section as a boundary condition in continuum modeling.

122

6.3

Details of continuum simulations

In moving boundary problems the interface treatment including shape, location, and velocity is highly complicated because the information on the interface need to be computed
along with the ﬁeld equations. Numerous techniques exist for tracking moving interfaces.
These techniques can be classiﬁed under three main categories a) moving grid (Lagrangian)
methods, b) Eulerian methods, and c) combined Lagrangian-Eulerian methods. All these
approaches have their relative merits and none of them is superior for all applications.
In this section, the problem of the contact line motion is modeled using a LagrangianEulerian algorithm. The Navier-Stokes (N-S) equations are solved with the ﬁnite diﬀerence
method on a ﬁxed Eulerian grid combined with tracking Lagrangian marker points assigned
on the ﬂuid/ﬂuid interface. The N-S equations provide the description of the pressure and
velocity distributions in the ﬂuid domain [149]. The marker points determine the time
evolution of the ﬂuid/ﬂuid interface [150]. Meanwhile, the molecular dynamics simulations
provide correct boundary conditions near the contact line (i.e., the slip length, shear stresses,
and friction coeﬃcients).

6.3.1

Fixed grid method for immersed objects

This method operates on a ﬁxed Cartesian mesh (Eulerian part) while the ﬂuid/ﬂuid interfaces move through the mesh (Lagrangian part). The advantage of the ﬁxed grid approach
is that the grid topology remains simple while large distortions of the interface take place.
In the immersed boundary method, originally used by Peskin [148] to model heart ﬂows,
the interface can be handled via a set of marker points.
Interface
Consider a ﬂuid domain, which is covered with a ﬁxed Cartesian grid. The interface between
Fluid 1 and Fluid 2 represented by a curve C. The grid spacing is h. The immersed boundary
denoted by C(t) (a curve in 2D case) is represented by k markers or interfacial points with
coordinate xk (s) with k = 1, 2, ..., K. The markers are uniformly distributed along C(t) at
a fraction of the grid spacing, 0.5 < ds < 1.5h. The ﬂuid/ﬂuid interface is parameterized as

123

a function of the arclength s by ﬁtting a local cubic polynomial to every four neighboring
points. Once the interface position is known, the normal vector and curvature are evaluated
and the material properties in each ﬂuid are assigned. The normal vector in 2D is given by
nx = −

ys
2 + y2
xs
s

and

ny =

xs
,
2 + y2
xs
s

(6.3)

where subscript s denotes d/ds. The curvature is then obtained by taking the divergent of
the normal vector:
κ=

· n.

(6.4)

Knowing the interface location with respect to the grid system, the material properties
are assigned in each ﬂuid phase based on a Heaviside step function:
β = β1 + (β2 − β1 )H(x − xk ),

(6.5)

where β is the material property such as ρ density and µ shear viscosity. The subscripts
1 and 2 denote Fluid 1 and Fluid 2, respectively. The function H(x − xk ) is the discrete
Heaviside step function deﬁned as:

 dim 1

(xm − (xm )k ) 1
π(xm − (xm )k )


1+
+ sin

2
d
π
d
m=1
H(x − xk ) =
 1



 0

if |x − xk |

d

,
if x − xk > +d
if x − xk < −d
(6.6)

where dim is the spatial dimension, d = 2h with h is the grid spacing, x is the grid coordinate,
and xk is the interfacial marker coordinate. The Heaviside step function allows a smooth
distribution of the material properties over a transition zone d. The transition zone is twice
the grid spacing d = 2h. Therefore, the interface has the ﬁnite thickness d. Once the
material properties are assigned to each ﬂuid, the ﬂuid ﬂow equations can be solved.
Immersed boundary treatment
The immersed boundary technique facilitates the communication between the moving markers, which represent the interface and the ﬁxed grid. In particular, since the location of the
marker points does not coincide with the grid points, the velocity ﬁeld stored at the cellcenter of each grid is interpolated to obtain the velocity of the interface points, and the
124

interface forces acting on the marker points is distributed to the nearby grid points using a
discrete Delta function. The Delta function is given by:
 dim

π(xm − (xm )k )
1

1 + cos
δ(x − xk ) =
2d
d
 m=1

0

if |x − xk |

d

.

(6.7)

otherwise

The delta function is the derivative of the Heaviside function and typically spreads over
the region of 4h.
Estimation of the force F : The force exerted by the ﬂuid/ﬂuid interface on the ﬂow is
added to the momentum equation. The integral of the source terms is expressed as:
F =

C(t)

γκˆ δ(x − xk )ds,
n

(6.8)

where the integral is taken over C(t), γ is the surface tension, κ is the curvature, n is the
ˆ
normal unit vector, δ(x − xk ) is the ’Dirac’ Delta function, x gives the location of the grid
points, xk gives the position of the markers as a function of arclength s. In the discretized
form, the force is estimated by:
γκk nk δ(x − xk )∆sk ,
ˆ

Fp =

(6.9)

k

The force at the point P , Fp , is evaluated based on the summation of all interfacial forces
of the marker points located inside a circle of radius 2h weighted by the Delta function.
Estimation of the interface velocity Uk : The velocity of a marker point or interface velocity denoted by Uk should satisfy the continuity condition. So Uk is:
Uk =

Ω

u δ(x − xk )dx,

(6.10)

where the integral is over the entire ﬂuid domain Ω, u is the ﬂuid velocity vector, and
δ(x − xk ) is the Delta function. In the discretized form, the velocity of the marker point is:
uij δ(x − xk )h2 ,

Uk =

(6.11)

ij

where h is the grid spacing, i and j are the grid indices and u is the ﬂuid velocity. The
velocity at the k th marker, Uk , is the sum of the ﬂuid velocity of the grid points located
inside a circle of radius 2h, weighted by the Delta function.
125

The interface velocity Uk is by deﬁnition the time derivative of the marker point position:
Uk =

∂x
.
∂t

(6.12)

Therefore, the interface is advected as:
xn+1 = xn +
k
k

n
t(Uk ).

(6.13)

The dynamics of two immiscible ﬂuids in a channel is governed by the incompressible
Navier-Stokes equations given below
· u = 0,

∂(ρu)
+
∂t

(6.14)

· (uρu) = − p + µ 2 u + F.

(6.15)

In these equations u is the ﬂuid velocity vector, p the pressure, F the body force, µ the
dynamic viscosity of ﬂuid, ρ the density of ﬂuid, and t time. The properties are constant in
each ﬂuid.
The Navier-Stokes equations are solved using the projection method on ﬁxed Cartesian
grids. The projection method is divided into three steps:
Step1: Solve the momentum equations without pressure. Evaluation of the intermediate
velocities ux and uz for a two dimensional problem:
˜
˜
˜
u − un
3(u ·
+
∆t

u)n − (u ·
2

u)n−1

=

˜
ν 2u n − ν 2u
+ F,
2

(6.16)

where n is the time step level and F is the source term.
The pressure equation is derived by assuming that the velocity satisﬁes the continuity
equation at n + 1 time step level.
Step 2: The following equation is solved for the pressure
· u n+1 = 0 ⇒

·(

1
1 n+1
p
)=
ρ
∆t

˜
· u.

(6.17)

Step 3: Combining the velocity and pressure information: Finally, the correction step is
done as follows
˜
u n+1 = u − ∆t(

126

1 n+1
p
).
ρ

(6.18)

6.3.2

Continuum Results

All variables in continuum simulations are normalized by the LJ length σ, time τ , and
energy ε scales. The cell dimensions are set to Lx = 100 and Lz = 25 in the x and z direcˆ
ˆ
tions, respectively. The system setup is designed to mimic the MD simulations with only
one ﬂuid/ﬂuid interface for simplicity. The schematic picture of the system is shown in
Fig. 6.14 (top).
The upper and lower walls are moving in the opposite directions with velocity U . At
the inlet and outlet it is assumed that the ﬂow enters the system parallel to the walls
uz |inlet = 0 and ∂ux |inlet = 0. In Fig. 6.14 (top) the red circles represent the initial location
∂x
of the Lagrangian marker points. The total number of the marker points is set to Np = 30 to
ensure an accurate resolution of the ﬂuid/ﬂuid interface with respect to the number of grids
in the z direction (Nz = 25). The motion of these marker points is governed by Eq. (6.13).
ˆ
Far away from the contact line, along the green line in the Fig. 6.14, the slip velocity uslip is
calculated using the Navier slip condition with βslip = 0.9 and µ = 2.0 extracted from the
MD simulations. At the contact lines indicated by the blue circles in the schematic picture,
the velocity ucontact is estimated from the following equation
ucontact = γ(cos(θs ) − cos(θ))/βCL ,

(6.19)

where βCL is the friction coeﬃcient at the contact line, γ = 3.7 the surface tension, and
θs = 90◦ the static contact angle. The value of βCL = 10.5 is selected by trial and error so
that the contact angle θ = θd in the converged solution stays within an error bar of about
5 degrees from the dynamic contact angle predicated by the MD simulations. Applying
the correct boundary condition for the rest of the domain (blue lines of in Fig. 6.14 (top))
is challenging. In this region neither the far-ﬁeld Navier condition nor the contact point
dynamics can be applied. Therefore, motivated by the shape of the slip velocity proﬁles
from the MD simulations (see Fig. 6.13(a)), a Delta function similar to Eq. (6.7) with the
maximum value of ucontact and a minimum value of uslip is used. The schematic shape of
velocity distribution is shown in Fig. 6.14 (bottom). The values of the ucontact and uslip
vary during the numerical iterations based on the temporary contact angle and the local
shear rate until the solution converges.
127

The velocity vectors and the streamlines are shown in Fig. 6.15 for U = 0.01, U = 0.05,
and U = 0.1. The streamline directions indicate that the ﬂuid can not penetrate into the
contact region, and, therefore, a circulation zone is formed in each side of the contact line,
which is shown by a red line. The dynamic contact angle as a function of the wall speed is
plotted in Fig. 6.16 and the results are compared with MD simulations. The dynamic contact
angle in the advancing ﬂuid phase increases with increasing wall speed. As discussed earlier
in Sec. 6.2 and shown in the Fig. 6.16, the value of the friction coeﬃcient βCL was chosen
to keep the contact angle extracted from the Navier-Stokes solution within an error bar of
about 5◦ from the MD results.
The slip velocity proﬁles from the Navier-Stokes solution are plotted in Fig. 6.17 by symbols and compared with MD results (dashed lines) for the selected wall speeds. The continuum results successfully reproduce the slip velocity proﬁles away from the contact line.
The results show qualitative agreement between the shape of the proﬁles near the contact
line. Note that one of the diﬀerences between the MD and continuum results is that the slip
velocity proﬁle from N-S solution is narrower near the contact region. This discrepancy can
be in part attributed to an ambiguity in deﬁning the width of the Delta function Eq. (6.7)
in the contact region shown in Fig. 6.14 (bottom). It is not clear at present whether the
region near the contact line should be about 3 − 5 σ (as deduced from the density proﬁles
and shear stress distribution) or even about 10 σ (the width of slip velocity proﬁles).

6.4

Conclusions

The problem of moving contact line between two immiscible ﬂuids on a smooth surface is
studied using molecular dynamics (MD) and continuum simulations. The ﬂow boundary
conditions near the moving contact line extracted from MD simulations and applied in the
continuum solution of the Navier-Stokes equation to reproduce velocity proﬁles and the
shape of the ﬂuid/ﬂuid interface. The algorithm based on the Irving-Kirkwood relation
for stress calculation provides an accurate description of shear stresses near ﬂuid/ﬂuid and
ﬂuid/solid interfaces. With increasing wall speed, the contact angle increases and the slip
velocity near the contact line is signiﬁcantly larger than in shear ﬂow of a single ﬂuid

128

phase. By implementing the slip boundary condition extracted from MD simulations, the
continuum method successfully reproduces the ﬂow pattern and velocity proﬁles.

129

γ
Fluid 1

Fluid 2

γ1 θ s

γ2
γ

Fluid 1

γ1 θ d

Fluid 2

γ2

U
Figure 6.1. The schematic representation of the contact line and the shape of the ﬂuid/ﬂuid
interface for the static (left) and dynamic (right) cases.
130

ucontact,βcontact
uslip ,βslip

U
Figure 6.2. The schematic picture of the moving contact line problem where the slip velocity
and the friction coeﬃcient far from the contact line are denoted by uslip and βslip and in
the contact region by ucontact and βcontact .

Fluid 1

Fluid 2

Fluid 1

Figure 6.3. Initial arrangement of ﬂuid molecules and wall atoms in the system that consists
of two immiscible ﬂuids conﬁned between ﬂat walls.

131

1

δ = -1

VLJ/ ε

0.5

δ=0
0

-0.5

δ=1
-1
1

1.5

r/ σ

2

Figure 6.4.
The modiﬁed Lennard-Jones potential Eq. (6.1) as a function of
monomer/monomer distance for diﬀerent values of the parameter δ.

Figure 6.5. A snapshot of an equilibrated system with two immiscible ﬂuids when both
walls are at rest.

132

(a)
4

ρσ

3

3
2
1
0

5

10

15

20

120

180

z/σ

25

(b)
0.8
0.7

ρσ3

0.6
0.5
0.4
0.3
0.2
0.1
0

60

x/σ

240

Figure 6.6. The averaged density proﬁles for ﬂuid 1 across the channel width in the z
ˆ
direction (a) and density proﬁles for ﬂuids 1 and 2 parallel to the walls in the x direction
ˆ
(b).
133

τ xx ,τ yy

τ zz

10.16
8.76
7.37
5.97
4.58
3.18
1.79
0.39
2.86
2.68
2.50
2.31
2.13
1.95
1.76
1.58

τ xx
τ yy
τ zz

Figure 6.7. The normal stress distributions in the x, y , and z directions in equilibrium case
ˆ ˆ
ˆ
(U = 0) calculated from the Irving-Kirkwood relation. The normal stresses are measured
in units of ε/σ 3 .

134

(b)

0.2
0.15

τxz σ3/ ε

0.1
0.05
0

-0.05
-0.1
-0.15
-0.2

40

50

x /σ

60

Figure 6.8. (a) The snapshot of the ﬂuid molecules near the solid wall within 5σ from
the contact line at equilibrium. The horizontal arrows near the contact line indicate the
direction of the tangential forces on each ﬂuid. (b) The shear stress distribution τxz (x)
along the wall as a function of the distance from the contact line.
135

Figure 6.9. The shape of the ﬂuid/ﬂuid interfaces as a function of the wall speeds for (a)
U = 0, (b) U = 0.05 σ/τ , and (c) U = 0.1 σ/τ .

136

Uτ/σ

0.1

Molecular Dynamics

0.05

0
90

95

100

105

θd

110

115

Figure 6.10. The contact angle averaged in two ﬂuid phases is plotted as a function of the
wall speeds U .

137

Figure 6.11. The snapshot of the system near the contact line for U = 0.1 σ/τ (top), the
averaged velocity vectors (bottom) obtained from MD simulations.

138

0.1

uslipτ / σ

Fluid 1
Fluid 2

0.05

0

0

50

x /σ

100

Figure 6.12. The slip velocity of the ﬁrst ﬂuid layer for upper and lower wall speeds U =
0.1 σ/τ . The velocity proﬁles are averaged separately in each ﬂuid phase as shown by
symbols. The dashed line is the slip velocity computed regardless of the ﬂuid type. The
black solid is a guide for the eye.

139

(a)
0.1

uslipτ / σ

U = 0.01 σ/τ
U = 0.05 σ/τ
U = 0.1 σ/τ

0.05

0
100

150

x /σ

200

(b)
0.1
U = 0.01 σ/τ
U = 0.05 σ/τ
U = 0.1 σ/τ

u(z) τ / σ

0.05

0

-0.05

-0.1

0

5

10

15

z /σ

20

25

Figure 6.13. (a) The slip velocity in the ﬁrst ﬂuid layer as a function of upper/lower wall
speeds. (b) The velocity proﬁles across the channel far from the contact lines in the z
ˆ
direction as a function of the upper/lower wall speeds.
140

U

U
z
x

u contact

u slip

u slip

Figure 6.14. The setup of the continuum simulations where the black solid lines show the
ﬁxed Cartesian grids, the red circles are the initial position of the marker points, and the
blue circles are the initial location of the contact points. The boundary conditions are
applied along the green regions away from the contact point and along blue regions close
to the contact point (top). The schematic velocity proﬁle in the bottom ﬁgure shows an
enlarged view of the velocity proﬁle inside top ﬁgure’s dashed-boxes.

141

(a)

U = 0.01

z

x

(b)

U = 0.05

(c)

U = 0.1

Figure 6.15. The results from the Navier-Stokes solution of the moving contact line as
a function of upper/lower wall speeds U = 0.01 (a), U = 0.05 (b), and U = 0.1 (c).
The velocity vectors are shown by blue arrows, streamlines by black lines, and ﬂuid/ﬂuid
interfaces by red lines.

142

Uτ/σ

0.1

Molecular Dynamics
Navier-Stokes

0.05

0
90

95

100

105

θd

110

115

Figure 6.16. The ﬁlled symbols represent the contact angle extracted from the Navier-Stokes
solution as a function of upper/lower wall speeds. The open symbols are the results from
the molecular dynamics simulations. Solid and dashed lines are a guide for the eye.

143

0.1

U = 0.01

uslipτ / σ

U = 0.05
U = 0.1

0.05

0

0

25

50

x /σ

75

100

Figure 6.17. The symbols and solid lines show the slip velocities near the moving contact
line obtained from the Navier-Stokes solution. The dashed lines are the slip velocities from
molecular dynamics simulations in the ﬁrst ﬂuid layer as a function of upper/lower wall
speeds.

144

BIBLIOGRAPHY

145

[1] G. Karniadakis, A. Beskok, and N. Aluru, Microﬂows and nanoﬂows: Fundamentals
and simulation (Springer, New York, 2005).
[2] H. A. Stone, A. D. Stroock, and A. Ajdari, “Engineering ﬂows in small devices: Microﬂuidics toward a lab-on-a-chip,” Annu. Rev. Fluid Mech. 36, 381 (2004).
[3] A. A. Darhuber and S. M. Troian, “Principles of microﬂuidic actuation by modulation
of surface stresses,” Annu. Rev. Fluid Mech. 37, 425 (2005).
[4] C. Neto, D. R. Evans, E. Bonaccurso, H. J. Butt, and V. S. J. Craig, “Boundary slip
in Newtonian liquids: A review of experimental studies,” Rep. Prog. Phys. 68, 2859
(2005).
[5] N. V. Churaev, V. D. Sobolev, and A. N. Somov, “Slippage of liquids over lyophobic
solid-surfaces,” J. Colloid Interface Sci. 97, 574 (1984).
´
[6] C. Cottin-Bizonne, S. Jurine, J. Baudry, J. Crassous, F. Restagno, and E. Charlaix,
“Nanorheology: An investigation of the boundary condition at hydrophobic and hydrophilic interfaces,” Eur. Phys. J. E 9, 47 (2002).
[7] L. Joly, C. Ybert, and L. Bocquet, “Probing the nanohydrodynamics at liquid-solid
interfaces using thermal motion,” Phys. Rev. Lett. 96, 046101 (2006).
[8] T. Schmatko, H. Hervet, and L. Leger, “Friction and slip at simple ﬂuid-solid interfaces:
The roles of the molecular shape and the solid-liquid interaction,” Phys. Rev. Lett. 94,
244501 (2005).
[9] R. Pit, H. Hervet, and L. Leger, “Direct experimental evidence of slip in hexadecane:
Solid interfaces,” Phys. Rev. Lett. 85, 980 (2000).
[10] Y. Zhu and S. Granick, “Limits of the hydrodynamic no-slip boundary condition,”
Phys. Rev. Lett. 88, 106102 (2002).
[11] E. Bonaccurso, H. J. Butt, and V. S. J. Craig, “Surface roughness and hydrodynamic
boundary slip of a newtonian ﬂuid in a completely wetting system,” Phys. Rev. Lett.
90, 144501 (2003).
[12] J. Sanchez-Reyes and L. A. Archer, “Interfacial slip violations in polymer solutions:
Role of microscale surface roughness,” Langmuir 19, 3304 (2003).
[13] T. Schmatko, H. Hervet, and L. Leger, “Eﬀect of nanometric-scale roughness on slip
at the wall of simple ﬂuids,” Langmuir 22, 6843 (2006).
[14] O. I. Vinogradova and G. E. Yakubov, “Surface roughness and hydrodynamic boundary
conditions,” Phys. Rev. E 73, 045302(R) (2006).
[15] K. B. Migler, H. Hervet, and L. Leger, “Slip transition of a polymer melt under shear
stress,” Phys. Rev. Lett. 70, 287 (1993).

146

[16] R. G. Horn, O. I. Vinogradova, M. E. Mackay, and N. Phan-Thien, “Hydrodynamic
slippage inferred from thin ﬁlm drainage measurements in a solution of nonadsorbing
polymer,” J. Chem. Phys. 112, 6424 (2000).
[17] Y. Zhu and S. Granick, “Rate-dependent slip of Newtonian liquid at smooth surfaces,”
Phys. Rev. Lett. 87, 096105 (2001).
[18] V. S. J. Craig, C. Neto, and D. R. M. Williams, “Shear-dependent boundary Slip in an
aqueous Newtonian liquid,” Phys. Rev. Lett. 87, 054504 (2001).
[19] C. H. Choi, K. J. A. Westin, and K. S. Breuer, “Apparent slip ﬂows in hydrophilic and
hydrophobic microchannels,” Phys. Fluids 15, 2897 (2003).
[20] J. Koplik, J. R. Banavar, and J. F. Willemsen, “Molecular dynamics of ﬂuid ﬂow at
solid interfaces,” Phys. Fluids A 1, 781 (1989).
[21] P. A. Thompson and M. O. Robbins, “Shear ﬂow near solids: Epitaxial order and ﬂow
boundary conditions,” Phys. Rev. A 41, 6830 (1990).
[22] P. A. Thompson, M. O. Robbins, and G. S. Grest, “Structure and shear response in
nanometer thick ﬁlms,” Israel Journal of Chemistry 35, 93 (1995).
[23] P. A. Thompson and S. M. Troian, “A general boundary condition for liquid ﬂow at
solid surfaces,” Nature (London) 389, 360 (1997).
[24] A. Jabbarzadeh, J. D. Atkinson, and R. I. Tanner, “Wall slip in the molecular dynamics
simulation of thin ﬁlms of hexadecane,” J. Chem. Phys. 110, 2612 (1999).
[25] J.-L. Barrat and L. Bocquet, “Large slip eﬀect at a nonwetting ﬂuid-solid interface,”
Phys. Rev. Lett. 82, 4671 (1999).
[26] M. Cieplak, J. Koplik, J. R. Banavar, “Boundary conditions at a ﬂuid-solid interface,”
Phys. Rev. Lett. 86, 803 (2001).
[27] N. V. Priezjev and S. M. Troian, “Molecular origin and dynamic behavior of slip in
sheared polymer ﬁlms,” Phys. Rev. Lett. 92, 018302 (2004).
[28] T. M. Galea and P. Attard, “Molecular dynamics study of the eﬀect of atomic roughness
on the slip length at the ﬂuid-solid boundary during shear ﬂow,” Langmuir 20, 3477
(2004).
[29] N. V. Priezjev, “Rate-dependent slip boundary conditions for simple ﬂuids,” Phys. Rev.
E 75, 051605 (2007).
[30] R. Khare, J. J. de Pablo, and A. Yethiraj, “Rheology of conﬁned polymer melt,” Macromolecules 29, 7910 (1996).
[31] A. Koike and M. Yoneya, “Chain length eﬀects on frictional behavior of conﬁned ultrathin ﬁlms of linear alkanes under shear,” J. Phys. Chem. B 102, 3669 (1998).

147

[32] A. Niavarani and N. V. Priezjev, “Slip boundary conditions for shear ﬂow of polymer
melts past atomically ﬂat surfaces,” Phys. Rev. E 77, 041606 (2008).
[33] J. P. Gao, W. D. Luedtke, and U. Landman, “Structures, solvation forces and shear of
molecular ﬁlms in a rough nano-conﬁnement,” Tribol. Lett. 9, 3 (2000).
[34] A. Jabbarzadeh, J. D. Atkinson, and R. I. Tanner, “Eﬀect of the wall roughness on slip
and rheological properties of hexadecane in molecular dynamics simulation of Couette
shear ﬂow between two sinusoidal walls,” Phys. Rev. E 61, 690 (2000).
[35] N. V. Priezjev and S. M. Troian, “Inﬂuence of periodic wall roughness on the slip
behaviour at liquid/solid interfaces: Molecular-scale simulations versus continuum predictions,” J. Fluid Mech. 554, 25 (2006).
[36] C. Kunert and J. Harting, “Roughness induced boundary slip in microchannel ﬂows,”
Phys. Rev. Lett. 99, 176001 (2007).
[37] N. V. Priezjev, “Eﬀect of surface roughness on rate-dependent slip in simple ﬂuids,” J.
Chem. Phys. 127, 144708 (2007).
[38] E. Lauga and H. A. Stone, “Eﬀective slip in pressure-driven Stokes ﬂow,” J. Fluid
Mech. 489, 55 (2003).
[39] S. C. Hendy, M. Jasperse, and J. Burnell, “Eﬀect of patterned slip on micro- and
nanoﬂuidic ﬂows,” Phys. Rev. E 72, 016303 (2005).
[40] N. V. Priezjev, A. A. Darhuber, and S. M. Troian, “Slip behavior in liquid ﬁlms on
surfaces of patterned wettability: Comparison between continuum and molecular dynamics simulations,” Phys. Rev. E 71, 041608 (2005).
[41] T. Qian, X. P. Wang, and P. Sheng, “Hydrodynamic slip boundary condition at chemically patterned surfaces: A continuum deduction from molecular dynamics,” Phys.
Rev. E 72, 022501 (2005).
[42] S. C. Hendy and N. J. Lund, “Eﬀective slip boundary conditions for ﬂows over nanoscale
chemical heterogeneities,” Phys. Rev. E 76, 066313 (2007).
[43] S. Richardson, “On the no-slip boundary condition,” J. Fluid Mech. 59, 707 (1973).
´
[44] C. Cottin-Bizonne, J.-L. Barrat, L. Bocquet, and E. Charlaix, “Low-friction ﬂows of
liquid at nanopatterned interfaces,” Nature Mat. 2, 237 (2003).
´
[45] C. Cottin-Bizonne, C. Barentin, E. Charlaix, L. Bocquet, and J.-L. Barrat, “Dynamics of simple liquids at heterogeneous surfaces: Molecular-dynamics simulations and
hydrodynamic description,” Eur. Phys. J. E 15, 427 (2004).
[46] M. Sbragaglia, R. Benzi, L. Biferale, S. Succi, and F. Toschi, “Surface roughnesshydrophobicity coupling in microchannel and nanochannel ﬂows,” Phys. Rev. Lett. 97,
204503 (2006).

148

[47] M. Sbragaglia and A. Prosperetti, “A note on the eﬀective slip properties for microchannel ﬂows with ultrahydrophobic surfaces,” Phys. Fluids 19, 043603 (2007).
[48] J. Ou, B. Perot, and J. P. Rothstein, “Laminar drag reduction in microchannels using
ultrahydrophobic surfaces,” Phys. Fluids 16, 4635 (2004).
[49] C. H. Choi and C. J. Kim, “Large slip of aqueous liquid ﬂow over a nanoengineered
superhydrophobic surface,” Phys. Rev. Lett. 96, 066001 (2006).
[50] P. Joseph, C. Cottin-Bizonne, J.-M. Benoit, C. Ybert, C. Journet, P. Tabeling, and
L. Bocquet, “Slippage of water past superhydrophobic carbon nanotube forests in microchannels,” Phys. Rev. Lett. 97, 156104 (2006).
[51] D. Einzel, P. Panzer, and M. Liu, “Boundary condition for ﬂuid ﬂow: Curved or rough
surfaces,” Phys. Rev. Lett. 64, 2269 (1990).
[52] P. Panzer, M. Liu, and D. Einzel, “The eﬀects of boundary curvature on hydrodynamic
ﬂuid ﬂow: Calculation of slip lengths,” Int. J. Mod. Phys. B 6, 3251 (1992)
[53] K. Kremer and G. S. Grest, “Dynamics of entangled linear polymer melts: A molecular
dynamics simulation,” J. Chem. Phys. 92, 5057 (1990).
[54] G. S. Grest and K. Kremer, “Molecular-dynamics simulation for polymers in the presence of a heat bath,” Phys. Rev. A 33, 3628 (1986).
[55] J. P. Boon and S. Yip, Molecular Hydrodynamics (McGraw-Hill, New York, 1980).
[56] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon, Oxford,
1987).
[57] J. C. Heinrich and D. W. Pepper, Intermediate Finite Element Method: Fluid Flow
and Heat Transfer Applications (Taylor and Francis, Philadelphia, 1999).
[58] J. H. Irving and J. G. Kirkwood, “The statistical mechanical theory of transport
processes. IV. The equations of hydrodynamics,” J. Chem. Phys. 18, 817 (1950).
[59] R. B. Bird, C. F. Curtiss, R. C. Armstrong, and O. Hassager, Dynamics of Polymeric
Liquids 2nd ed. (Wiley, New York, 1987).
[60] Z. Xu, J. J. de Pablo, and S. Kim, “Transport properties of polymer melts from nonequilibrium molecular dynamics,” J. Chem. Phys. 102, 5836 (1995).
[61] J. T. Bosko, B. D. Todd, and R. J. Sadus, “Viscoelastic properties of denrimers in the
melt from nonequilibrium molecular dynamics,” J. Chem. Phys. 121, 12050 (2004).
[62] L. M. Hocking, “A moving ﬂuid interface on a rough surface,” J. Fluid Mech. 76, 801
(1976).
[63] E. O. Tuck and A. Kouzoubov, “A laminar roughness boundary condition,” J. Fluid
Mech. 300, 59 (1995).

149

[64] T. Aoyagi, J. Takimoto, and M. Doi, “Molecular study of polymer melt conﬁned between walls,” J. Chem. Phys. 115, 552 (2001).
[65] I. Bitsanis and G. Hadziioannou, “Molecular-dynamics simulations of the structure and
dynamics of conﬁned polymer melts,” J. Chem. Phys. 92, 3827 (1990).
[66] S. J. Plimpton, J. Comp. Phys. 117, 1 (1995); see also URL http://lammps.sandia.gov.
[67] L. Bocquet and J.-L. Barrat, “Flow boundary conditions from nano- to micro-scales,”
Soft Matter 3, 685 (2007).
[68] V. Mhetar and L. A. Archer, “Slip in entangled polymer solutions,” Macromolecules
31, 6639 (1998).
[69] K. M. Awati, Y. Park, E. Weisser, and M. E. Mackay, “Wall slip and shear stresses
of polymer melts at high shear rates without pressure and viscous heating eﬀects,” J.
Non-Newton. Fluid Mech. 89, 117 (2000).
[70] X. Zhou, D. Andrienko, L. Delle Site, and K. Kremer, “Flow boundary conditions for
chain-end adsorbing polymer blends,” J. Chem. Phys. 123, 104904 (2005).
[71] C. Pastorino, K. Binder, T. Kreer, and M. Muller, “Static and dynamic properties of
the interface between a polymer brush and a melt of identical chains,” J. Chem. Phys.
124, 064902 (2006).
[72] S. Barsky and M. O. Robbins, “Molecular dynamics study of slip at the interface
between immiscible polymers,” Phys. Rev. E 63, 021801 (2001).
[73] M. Tsige and G. S. Grest, “Molecular dynamics simulation of solventpolymer interdiffusion: Fickian diﬀusion,” J. Chem. Phys. 120, 2989 (2004).
[74] J. L. Barrat and J. P. Hansen, Basic concepts for simple and complex liquids (Cambridge
University Press, Cambridge, 2003).
[75] J. N. Israelachvili, Intermolecular and Surface Forces, 2nd ed. (Academic Press, San
Diego, 1992).
[76] M. O. Robbins, “Fluid dynamics: All stressed out,” Nature (London) 389, 331 (1997).
[77] Z. Xu, R. Khare, J. J. de Pablo, and S. Kim, “On the calculation of transport properties
of polymer melts from nonequilibrium molecular dynamics,” J. Chem. Phys. 106, 8285
(1997).
[78] E. D. Smith, M. O. Robbins, and M. Cieplak, “Friction on adsorbed monolayers,” Phys.
Rev. B 54, 8252 (1996).
[79] M. S. Tomassone, J. B. Sokoloﬀ, A. Widom, J. Krim, “Dominance of phonon friction
for a xenon ﬁlm an a silver (111) surface,” Phys. Rev. Lett. 79, 4798 (1997).
[80] C. Pastorino, T. Kreer, M. Muller, and K. Binder, “Comparison of dissipative particle
dynamics and Langevin thermostats for out-of-equilibrium simulations of polymeric
systems,” Phys. Rev. E 76, 026706 (2007).
150

[81] D. J. Beebe, G. A. Mensing, and G. M. Walker, “Physics and applications of microﬂuidics in biology,” Annu. Rev. Biomed. Eng. 4, 261 (2002).
[82] T. M. Squires and S. R. Quake, “Microﬂuidics: Fluid physics at the nanoliter scale,”
Rev. Mod. Phys. 77, 977 (2005).
[83] V. Mengeaud, J. Josserand, and H. H. Girault, “Mixing processes in a zigzag microchannel: Finite element simulations and optical study,” Anal. Chem. 74, 4279 (2002).
[84] A. Gigras and S. Pushpavanam, “Early induction of secondary vortices for micromixing
enhancement,” Microﬂuid. Nanoﬂuid. 5, 89 (2008).
[85] H. M. Metwally and R. M. Manglik, “Enhanced heat transfer due to curvature-induced
lateral vortices in laminar ﬂows in sinusoidal corrugated-plate channels,” Int. J. Heat
Mass Transfer 47, 2283 (2004).
[86] P. E. Geyer, N. R. Rosaguti, D. F. Fletcher, and B. S. Haynes, “Thermohydraulics
of square-section microchannels following a serpentine path,” Microﬂuid. Nanoﬂuid. 2,
195 (2006).
[87] N. R. Rosaguti, D. F. Fletcher, and B. S. Haynes, “Low-Reynolds number heat transfer
enhancement in sinusoidal channels,” Chem. Eng. Sci. 62, 694 (2007).
[88] J. H. Forrester and D. F. Young, “Flow through a converging-diverging tube and its
implications in occlusive vascular disease. 1. Theoretical development,” J. Biomech. 3,
297 (1970).
[89] D. M. Wootton and D. N. Ku, “Fluid mechanics of vascular systems, diseases, and
thrombosis,” Annu. Rev. Biomed. Eng. 1, 299 (1999).
[90] S. A. Berger and L. D. Jou, “Flows in stenotic vessels,” Annu. Rev. Fluid Mech. 32,
347 (2000).
[91] C. Cottin-Bizonne, B. Cross, A. Steinberger, and E. Charlaix, “Boundary slip on
smooth hydrophobic surfaces: Intrinsic eﬀects and possible artifacts,” Phys. Rev. Lett.
94, 056102 (2005).
[92] V. P. Sokhan, D. Nicholson, and N. Quirke, “Fluid ﬂow in nanopores: An examination
of hydrodynamic boundary conditions,” J. Chem. Phys. 115, 3878 (2001).
[93] A. Niavarani and N. V. Priezjev, “Rheological study of polymer ﬂow past rough surfaces
with slip boundary conditions,” J. Chem. Phys. 129, 144902 (2008).
[94] K. M. Jansons, “Determination of the macroscopic (partial) slip boundary condition for
a viscous ﬂow over a randomly rough surface with a perfect slip microscopic boundary
condition,” Phys. Fluids 31, 15 (1988).
[95] M. J. Miksis and S. H. Davis, “Slip over rough and coated surfaces,” J. Fluid Mech.
273, 125 (1994).

151

[96] K. Sarkar and A. Prosperetti, “Eﬀective boundary conditions for Stokes ﬂow over a
rough surface,” J. Fluid Mech. 316, 223 (1996).
[97] C. Y. Wang, “Flow over a surface with parallel grooves,” Phys. Fluids 15, 1114 (2003).
[98] N. Lecoq, R. Anthore, B. Cichocki, P. Szymczak, and F. Feuillebois, “Drag force on a
sphere moving towards a corrugated wall,” J. Fluid Mech. 513, 247 (2004).
[99] J. C. F. Chow and K. Soda, “Laminar ﬂow in tubes with constriction,” Phys. Fluids
15, 1700 (1972).
[100] G. L. Bordner, “Nonlinear analysis of laminar boundary layer ﬂow over a periodic
wavy surface,” Phys. Fluids 21, 1471 (1978).
[101] S. Taneda, “Visualization of separating Stokes ﬂows,” J. Phys. Soc. Jpn. 46, 1935
(1979).
[102] I. Sobey, “On ﬂow through furrowed channels. Part 1. Calculated ﬂow patterns,” J.
Fluid Mech. 96, 1 (1980).
[103] T. Nishimura, Y. Ohori, and Y. Kawamura, “Flow characteristics in a channel with
symmetric wavy wall for steady ﬂow,” J. Chem. Eng. Japan 17, 466 (1984).
[104] S. Tsangaris and D. Potamitis, “On laminar small Reynolds-number ﬂow over wavy
walls,” Acta Mechanica 61, 109 (1986).
[105] G. Leneweit and D. Auerbach, “Detachment phenomena in low Reynolds number
ﬂows through sinusoidally constricted tubes,” J. Fluid Mech. 387, 129 (1999).
[106] H. Zhou, R. J. Martinuzzi, R. E. Khayat, A. G. Straatman, and E. Abu-Ramadan,
“Inﬂuence of wall shape on vortex formation in modulated channel ﬂow,” Phys. Fluids
15, 3114 (2003).
[107] C. Pozrikidis, “Creeping ﬂow in two-dimensional channels,” J. Fluid Mech. 180, 495
(1987).
[108] M. Scholle, A. Wierschem, and N. Aksel, “Creeping ﬁlms with vortices over strongly
undulated bottoms,” Acta Mechanica 168, 167 (2004).
[109] J. C. Heinrich and R. S. Marshall, “Viscous incompressible ﬂow by a penalty function
ﬁnite element method,” Computers and Fluids 9, 73 (1981).
[110] T. J. R. Hughes, W. K. Liu, and A. Brooks, “Finite element analysis of incompressible
viscous ﬂows by the penalty function formulation,” J. Comp. Phys. 30, 1 (1979).
[111] T. B. Benjamin, “Shearing ﬂow over a wavy boundary,” J. Fluid Mech. 6, 161 (1959).
[112] F. Pan and A. Acrivos, “Steady ﬂows in rectangular cavities,” J. Fluid Mech. 28, 643
(1967).
[113] H. K. Moﬀatt, “Viscous and resistive eddies near a sharp corner,” J. Fluid Mech. 18,
1 (1964).
152

[114] H. Schlichting and K. Gersten, Boundary-Layer Theory (Springer, Berlin, 2000).
[115] P. Koumoutsakos, “Multiscale ﬂow simulations using particles,” Annu. Rev. Fluid
Mech. 37, 457 (2005).
[116] D. A. Nield and A. Bejan, Convection in Porous Media (Springer, New York, 2006).
[117] J. Baudry, E. Charlaix, A. Tonck, and D. Mazuyer, “Experimental evidence for a large
slip eﬀect at a nonwetting ﬂuid-solid interface,” Langmuir 17, 5232 (2001).
[118] Y. Zhu and S. Granick, “No-slip boundary condition switches to partial slip when
ﬂuid contains surfactant,” Langmuir 18, 10058 (2002).
[119] U. Ulmanella and C.-M. Ho, “Molecular eﬀects on boundary condition in micro/nanoliquid ﬂows,” Phys. Fluids 20, 101512 (2008).
[120] U. Heinbuch and J. Fischer, “Liquid ﬂow in pores: Slip, no-slip, or multilayer sticking,”
Phys. Rev. A 40, 1144 (1989).
[121] J.-L. Barrat and L. Bocquet, “Inﬂuence of wetting properties on hydrodynamic boundary conditions at a ﬂuid/solid interface,” Faraday Discuss. 112, 119 (1999).
[122] S. C. Yang and L. B. Fang, “Eﬀect of surface roughness on slip ﬂows in hydrophobic
and hydrophylic microcannels by molecular dynamics simulation,” Molecular Simulation 31, 971 (2005).
[123] A. Martini, H. Y. Hsu, N. A. Patankar, and S. Lichter, “Slip at high shear rates,”
Phys. Rev. Lett. 100, 206001 (2008).
[124] N. V. Priezjev, “Shear rate threshold for the boundary slip in dense polymer ﬁlms,”
Phys. Rev. E 80, 031608 (2009).
[125] S. T. O’Connell and P. A. Thompson, “Molecular dynamics-continuum hybrid computations: A tool for studying complex ﬂuid ﬂows,” Phys. Rev. E 52, R5792 (1995).
[126] E. G. Flekkøy, G. Wagner, and J. Feder, “Hybrid model for combined particle and
continuum dynamics,” Europhys. Lett. 52, 271 (2000).
[127] X. B. Nie, S. Y. Chen, W. E, and M. O. Robbins, “A continuum and molecular
dynamics hybrid method for micro- and nano-ﬂuid ﬂow,” J. Fluid Mech. 500, 55 (2004).
[128] Q. Li and G.-W. He, “An atomistic-continuum hybrid simulation of ﬂuid ﬂows over
superhydrophobic surfaces,” Biomicroﬂuidics 3, 022409 (2009).
[129] X. B. Nie, S. Y. Chen, and M. O. Robbins, “Hybrid continuum-atomistic simulation
of singular corner ﬂow,” Phys. Fluids 16, 3579 (2004).
[130] N. G. Hadjiconstantinou, “Hybrid atomistic-continuum formulations and the moving
contact-line problem,” J. Comput. Phys. 154, 245 (1999).
[131] W. Q. Ren and W. E, “Heterogeneous multiscale method for the modeling of complex
ﬂuids and microﬂuidics,” J. Comput. Phys. 204, 1 (2005).
153

[132] A. Niavarani and N. V. Priezjev, “The eﬀective slip length and vortex formation in
laminar ﬂow over a rough surface,” Phys. Fluids 21, 052105 (2009).
[133] W. Shyy, H.S. Udaykumar, M. M. Rao, and R. W. Smith, Computational Fluid Dynamics with Moving Boundaries (Dover Publications, Inc., New York, 1996).
[134] G. de Gennes, “Wetting: Statics and dynamics,” Rev. Mod. Phys. 57, 827 (1985).
[135] A. L. Yarin, “Drop impact dynamics: Splashing, spreading, receding, bouncing,”
Annu. Rev. Fluid Mech. 38, 159 (2006).
[136] C. Huh and L. E. Scriven, “Hydrodynamic model of steady movement of a
solid/liquid/ﬂuid contact line,” J. Colloid Interface Sci. 35, 85 (1971).
[137] E. B. Dussan V., “On the spreading of liquids on solid surfaces: Static and dynamic
contact lines,” Annu. Rev. Fluid Mech. 11, 371 (1979).
[138] R. G. Cox, “The dynamics of the spreading of liquids on a solid surface. Part 1.
Viscous ﬂow,” J. Fluid Mech. 168, 169 (1986).
[139] J. Koplik, J. R. Banavar, and J. F. Willemsen, “Molecular dynamics of Poiseuille ﬂow
and moving contact lines,” Phys. Rev. Lett. 60, 1282 (1988).
[140] P. A. Thompson and M. O. Robbins, “Simulations of contact-line motion: Slip and
the dynamic contact angle,” Phys. Rev. Lett. 63, 766 (1989).
[141] P. A. Thompson, W. B. Brinckerhoﬀ and M. O. Robbins, “Microscopic studies of
static and dynamic contact angles,” J. Adhesion Science and Technology 7, 535 (1993).
[142] N. G. Hadjicostantinou, “Combining atomistic and continuum simulations of contactline motion,” Phys. Rev. E 59, 2475 (1999).
[143] T. Qian, X. P. Wang, and P. Sheng, “Molecular scale contact line hydrodynamics of
immiscible ﬂows,” Phys. Rev. E 68, 016306 (2003).
[144] W. Ren and W. E, “Boundary conditions for the moving contact line problem,” Phys.
Fluids 19, 022101 (2007).
[145] T. Ito, Y. Hirata, and Y. Kukita, “Fluid epitaxialization eﬀect on velocity dependence
of dynamic contact angle on molecular scale,” J. Chem. Phys. 132, 054702 (2010).
[146] A. Niavarani and N. V. Priezjev, “Modeling the combined eﬀect of surface roughness
and shear rate on slip ﬂow of simple ﬂuids,” Phys. Rev. E 81, 011606 (2010).
[147] C. Denniston and M. O. Robbins, “Mapping molecular models to continuum theories
for partially miscible ﬂuids,” Phys. Rev. E 69, 021505 (2004).
[148] C. S. Peskin, “Numerical analysis of blood ﬂow in the heart,” J. Comput. Phys. 25,
220 (1977).

154

[149] M. Griebel, T. Dornseifer, and T. Neunhoeﬀer, Numerical simulation in ﬂuid dynamics: A practical introduction (Society for industrial and applied mathematics, Philadelphia, 1998).
[150] W. Shyy, M. Francois, H. Udaykumar, N. N’dri, and R. Tran-son-tay, “Moving boundaries in micro-scale bioﬂuid dynamics,” Appl. Mech. Rev. 54, 405 (2001).

155