Rheological properties of thin liquid films by molecular dynamics simulations

j,m~d lqon-l~mui~ l~id ELSEVIER

J. Non-Newtonian Fluid Mech., 69 (1997) 169-193

Rheological properties of thin liquid films by molecular dynamics simulations A. Jabbarzadeh *, J.D. Atkinson, R.I. Tanner Department of Mechanical and Mechatronic Engineering, The University of Sydney, Sydney, N.S. I4/. 2006, Australia Received 13 May 1996; accepted 16 September 1996

Abstract

In this paper we present the molecular dynamics simulations of thin fluids films sheared in Couette flow geometry between two structured plane walls. An NVT ensemble of atoms was chosen and simulation conducted in isothermal conditions. To keep the temperature at the required level a Gaussian thermostat was employed. This method was shown to be superior to the simple velocity rescaling method, especially at high shear rates. The Gaussian thermostat method gave results for viscosity in good agreement with the results of other researchers who used the reservoir method. The results for density and velocity profiles were obtained for a wide range of simulation parameters. The effects of shear rate and wall-fluid interaction strength were investigated in detail over a wide range of parameters. The material functions and normal stress differences were also obtained and the effects of shear rate and wall strength parameter on these properties were studied. The effect of film thickness on the viscosity was investigated and was compared with what we found for bulk fluid using the SLLOD algorithm, the existence of a non-Newtonian region with shear-thinning effect is found and examined for various films. The results suggest an increase in viscosity for thinner films in the Newtonian regime, though this is valid only for a limited range of wall-fluid interaction strength. A decrease in viscosity was also observed when the attraction force of the wall was increased. © 1997 Elsevier Science B.V.

Keywords: Molecular dynamics; Thin films; NEMD; Shear flow; Wall slip; Rheology; Gaussian thermostat

1. Introduction

The behaviour of the flow in the vicinity of an impermeable wall has been a subject of interest for many years. The understanding of this behaviour has great importance in lubrication, wetting and coating problems. Israelachvili and coworkers [1-3] and Chan and Horn [4] have experimentally investigated the structure and dynamics of lubrication films using the SFA (surface force apparatus). Computer * Corresponding author. 0377-0257/97/$17.00 © 1997 Elsevier Science B.V. All rights reserved. PII S 0 3 7 7 - 0 2 5 7 ( 9 6 ) 0 1 5 2 0 - 0

170

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

simulations have also been employed as a virtual microscopic laboratory to examine the microstructure and transport properties of the fluid confined between solid walls separated by only a few molecular diameters. Molecular dynamics simulation has been widely used for this purpose both for equilibrium and nonequilibrium conditions. These simulations reveal the inhomogeneous nature of the fluid in proximity to the walls [5-8]. They also examine the non-slip boundary condition which is generally accepted for fluid mechanics problems [9,10]. These simulations reveal the layered structure of the fluid atoms in the vicinity of the solid boundaries. The results also suggest that the fluid properties in channels where wall separation is comparable to the molecular diameter are different from those of the bulk fluid [2,4]. One of the problems in the moving boundary shear simulations is how to remove the heat from the system. In most of the simulations this is done by simply rescaling the peculiar momenta of the particles. This however is not effective at high shear rates where the amount of generated heat is much more than can be handled by a simple rescaling scheme. Another possibility is to thermostat the particles on the walls to let the heat flow to them from the fluid particles. Liem et al. [11] have employed this method but have not been able to increase the shear rate beyond a certain limit because of the penetration of the fluid particles into the walls. The use of the Gaussian thermostat, which is the subject of our study in this paper, is another option. However, in using this thermostat we should not make any assumption as to the form of the velocity profile which is normally assumed to be linear from hydrodynamics theory. This linear streaming velocity profile assumption is not valid for the inhomogeneous system that we study. This assumption also remains highly dubious at high Reynolds numbers [12]. In this paper we will examine the fluid dynamic and static properties for a fluid consisting of Lennard-Jones spheres confined between two solid walls and undergoing a Couette shear flow. A N V T ensemble of atoms where temperature (T), volume (V) and the number of atoms (N) are kept constant during the course of the simulation has been considered. The structural and dynamic properties of the sheared fluid have been analysed, and compared with previous results. We also put stronger emphasis on the rheological properties of the fluid under investigation. In this regard the material functions of the fluid are examined, and normal stress differences are also obtained for a number of different geometries. The existence of a non-Newtonian regime was found and the effect of the wall properties and other simulation parameters on the material functions is established. The simulations have been done in such a way that the average density of the system remains constant for all the geometries. The effect of the system size has also been considered to make sure that different geometries are comparable.

2. Simulation details 2.1. Simulation in the presence o f the walls

In this situation the fluid is confined between two solid walls. These walls are parallel to the X Y plane and Couette flow is generated by moving two walls in opposite directions along the X axis with the same velocity so that shear is applied in the X Z plane. The usual periodic

boundary conditions have been applied in the X and Y directions. Interaction between the fluid atoms is governed by a 6-12 Lennard-Jones potential described by

A. Jabbarzadeh et a l . / J . Non-Newtonian Fluid Mech. 69 (1997) 169-193

~bu(r) = 4e

[(~)12- (~)61-q~shift,

171

~( ~ -( 0") 61

qSshirt(r) = 4e__---a_'" L \ r c /

-

"

In order to reduce the computational effort this potential is truncated at r~ = 2.5o. To preserve the continuity of the potential it is shifted at r = rc so that it goes smoothly to zero as it approaches re. In Eq. (1) a and ~: are L e n n a r d - J o n e s length and energy parameters and r is the distance between pairs of atoms. The walls comprise atoms of a [100] plane of a face-centred cubic lattice. Each atom on the wall is attached to a stiff spring on its lattice position. The density of the wall atoms is Pw = 1; however, this can be adjusted by packing them at different distances from each other. The dimensions of the wall are 8,~f20- in the X direction and 5xfl20- in the Y direction, where 0- is the molecular diameter of the fluid particles, so that it can accommodate 80 atoms of the same diameter. For larger systems the wall dimensions are 14x~0" and 7xf20" in the X and Y directions respectively with 196 atoms on each wall. The confined fluid has a density of p -- 0.555 and a temperature of T = 1.2 where the usual reduced L e n n a r d - J o n e s units have been used as shown in Table 1. We have to examine the fluid properties for different wall separations while keeping the above mentioned fluid average density. Depending on the wall to wall distance, different numbers of fluid atoms are put in the simulation box to fulfil this requirement. The number of wall atoms has been the same for most simulations except for those with wall separations of 4, 7 and 10 molecular diameters and larger where the number is increased to ensure that the size of the wall is comparable to the wall separation. In Table 2 the relevant data have been given for different wall separations. The same potential has been used for the interaction of wall-fluid atoms with the wall length parameter 0"w= 0". Depending on the wall-fluid interaction strength, different values for the wall energy parameter ew have been considered. H o o k e a n springs attach the wall atoms to their lattice sites with a potential of 1

(2)

~Sspring --- ~ H R 2,

Table 1 L e n n a r d - J o n e s reduced units; cr, m and e are fluid molecular diameter, mass and energy parameter respectively Quantity

Unit

Length Mass Energy Shear rate Time Velocity Density Stress Viscosity Normal stress coefficient

a m t:

(g/m~T2)t,2 (mG2/g)I,'2

(e/m)l 2 a-3 a/a 3 (me ) l /2 / a 2 m/a

172

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

Table 2 Simulation parameters for various simulated systems Wall separation (Z)

Density (p)

Fluid atoms (Nf)

Wall atoms (Nw)

Temp. (T)

2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.5 9.0 10.0 15.0 20.0

0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555 0.555

218 272 133 156 435 200 222 244 266 289 761 333 377 400 1088 1632 2176

392 392 160 160 392 160 160 160 160 160 392 160 160 160 392 392 392

1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2 1.2

where H is the spring stiffness and R the distance between the wall atoms order to keep the mean average displacement of wall atoms at a low level and state, a stiff spring with H = 6000 e o-- 2 has been chosen. An N V T ensemble has been used for the molecular dynamics (MD) temperature, volume and the number of atoms are kept constant during simulation. The equations of motions for the fluid particles are i', =

t~,

b , = F~ -

f(p,

- nxUx,,),

and their sites. In to mimic the solid simulation where the course of the

(3)

m

where p~ and ri are atomic laboratory momenta and positions and Fg represents the interatomic forces. Ux, i is the average streaming velocity at the position of particle i and is evaluated separately in the simulation as the time goes on. This lets the velocity profile develop without interference from a linear-profile-biased thermostat, nx is the unit vector in the direction of the X axis and ( is the thermostatting multiplier which keeps the kinetic temperature fixed and is described as

y F¢,- y. ~',Px,P:, ( =

(4)

All p in Eq. (4) are peculiar (thermal) momenta; Px, and P~i are the X and Z components of the peculiar momentum of atom i and ~i = (SU/SZ):~ is the local shear rate determined at the Z position of particle i. Here U is the time average streaming velocity with the average being taken

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

173

over all previous time. We compute the streaming velocity by a method which will be described in Section 3.1. Use of this Gaussian thermostat [12,13] makes no assumption as to the form of the velocity profile. The streaming velocity profile is not known in advance so we let the velocity profile develop during the course of the simulation. However we assume there is no streaming velocity in Y and Z directions; this must be so by symmetry. Besides this thermostatting method, fluid particle velocity components in the Z and Y directions have been rescaled as a standard method for removing heat. Through the rescaling method alone is not effective at high shear rates, Liem et al. [11] have used the rescaling method to remove the heat from the boundary particles. Their method however is only effective for limited shear rates up to 2 in reduced units. Attempts to apply their method beyond this limit lead to fluid particles penetrating beyond the wall. This can be avoided by increasing the spring stiffness but this makes the wall particles too stiff and unable to conduct the heat out from the fluid. The Gaussian thermostat makes it possible to simulate the problem in isothermal conditions for higher shear rates. The equations of motion were integrated by a fifth-order Gear predictor-corrector algorithm. The time step used for the simulation was 0.00253 where v (m~72/c) 1/2. To make sure that the chosen time step was small enough, we performed simulations of the same length at smaller time steps up to 0.0005v but no significant change was observed in the results for material functions, particularly for the viscosity. Each simulation was performed for 253 for equilibrium purposes. It was established that a stable density profile developed and that a steady level of energy was achieved to confirm the equilibrium state. Then the program was run for 3003 and time averages of the relevant properties were computed during the course of the simulation. For low shear rates the simulations were performed for longer times of 6003. In order to increase the efficiency of the program and reduce the computing time the Verlet neighbour list method was employed with the list updated every 50 time steps. =

2.2. S L L O D

algorithm

In order to compare the results with the properties of bulk fluid a nonequilibrium molecular dynamics simulation was conducted using the SLLOD algorithm. The algorithm was developed by Evans and Morriss [12,13] and uses a sliding brick periodic boundary condition method to apply the shear. A Gaussian thermostat keeps the temperature constant. For this algorithm simulations were performed for a system of 1372 atoms with 0.0043 time steps.

3. Results and discussion

The simulations were performed and results were obtained for a number of time average fluid properties. The properties examined were density and velocity profiles, first and second normal stress differences and material functions including fluid viscosity (r/) and first and second normal stress coefficients ~bl and ~'2 [14].

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

174

3. I. Density and velocity profiles To determine density we divided the region between the walls into a sufficient number of slices. The number of slices was adjusted so that the width of each slice was 0.1 a. We calculated the time average density for each slice. Similarly we added the X components of the velocities for particles in each slice and then averaged over the number of particles in that slice at each time step. Finally velocity profiles were obtained by computing the time average velocity for each slice during the course of simulation. Thus average density for each slice can be described as

P .... =

/

i=, Ns[Z,(t)l Av

/

'

(5)

where p .... is the time average density for each slice, AV is the volume of the slice and Ns[Zi(t)] is a counting operator, which is described by

{Ns[Zi(t)] = 1 if ( s - 1) AZ < Zi < s AZ, Ns[Z~(t)] = 0 otherwise;

(6)

it simply counts the number of particles which happen to be in that particular slab at time t. Likewise the time average velocity in the X direction for each slice will be as follows where u~x is the X component of particle i's velocity and the angle brackets denote the usual time average:

(7)

3. I. 1. Effect of the wall separation on density and velocity profiles In order to examine the effect of the film thickness, simulations were done for different wall distances ranging from 2.00- to 20o-. All the simulations were performed for the same shear rate of 0.2(~/m0-2) ~/2 and wall strength parameter ew = 1.0~. Velocity profiles were normalised by the wall velocity for that particular wall separation so that a velocity of 1.0 is equal to the wall velocity. The results are shown in Fig. 1. It can be seen from Fig. 1 that the density of the fluid particles in the vicinity of the walls is much higher than at other places. The sharp jumps in the density profiles indicate the layered structure of the fluid particles at that point. This layering effect is more pronounced at smaller separations but the density of these layers decreases as we increase the fluid film thickness. Although the strength of the layering diminishes at larger wall separations one can observe two layers of fluid particles in the proximity of each wall for all systems except for those where only a few fluid layers can be accommodated in the space between the walls, as we see for Z = 2.5-4.00- where only 2 - 3 layers can be formed in the slit. Increasing the wall separation, one can observe that the density of the fluid in the middle of the slit approaches about 0.555 which is the average density of the confined fluid. Some snapshots taken from the simulation box can give us a better idea of the layering effect. Fig. 2 shows these snapshots for different wall separations. Snapshots in Fig. 2 illustrate the layering effect. For a

A. Jabbarzadeh et al.//J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

Z=2.5~ ~:~=1.0

8

o

Z=3.5~ ¢. =1.0

Z=3.0¢ ¢. =1.0

"

Z=4.0~ ¢=1,0

Z=4.5o

175

~3

Z=5.0d ¢=1.0

~ =l,O

'°i

~oo

D

Z=6.5~ £ =l,()

Z=6.0o C =l.O

Z=5.5c~ ¢~ =1.() ,0

°

~'°°~

~

~oo ,o

i. z (o) Z=7.0~ ¢. = [ 111

Z (o1 Z=8.5o %=1.0

Z=7.5o %=1.0

,o

~oo C3 ,o

Z=9.0o £ =l.()

Z (o)

) (o)

Z=IO.O~ ~: =l.(I

z=20.O~ ~ =l.O

> D ,o

~o

:~ (o}

(o)

;

Fig. 1. Velocity and density profiles for a number of different wall separations. All profiles are for the same shear rate of 0.2(~/m~2) 12 and wall strength parameter e~ = 1.0e,. The effect of film thickness on the velocity profiles is more pronounced for Z up to 4 molecular diameters.

176

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

wall separation of 2o- there is only one layer at the middle of the slit. The n u m b e r of layers for Z = 4o- is three with stronger layering in the proximity of the walls. For Z = 7o- one can see the layering near the walls but it is weaker than in the two other thinner films. Note a depletion zone in the vicinity of each wall. F r o m the density profiles in Fig. 1 it can be seen that the depth of the depletion zone is about 0.7o-. This means that fluid particles cannot effectively come closer than 0.7o" to the walls or, equivalently, that the effective wall position is 0.7o- further in than that shown in Fig. 2. The distance between the walls also greatly affects the flow velocity in the X direction. In Fig. 1 for the velocity profile the solid line represents the perfect linear velocity which is expected from hydrodynamics theory. The velocity profile is approximately linear for Z = 20o.. As we decrease the film thickness, velocity profiles begin to deviate from linearity. This trend gets stronger as films become thinner and eventually velocity profiles become very distorted for films thinner than 4o-. Qualitatively the results for density and velocity profiles are in agreement with results of Refs. [5,15]. Our results are not comparable quantitatively with the result of these references because their fluid was of different density and temperature. 3.1.2. Density and velocity profiles and shear rate

The way the shear rate affects the density and velocity profiles could be of interest in the study of slip and non-slip boundary conditions. To examine this issue we performed simulations for a range of shear rates between 0.05 and 3.5(e/mo.2) 1/2 for films of 4o., 7o- and 10o- thickness. The results are illustrated in Fig. 3. It can be seen from Fig. 3 that the shear rate strongly affects the velocity profiles. For all three films density profiles remain more or less the same; there is only a slight change in the width of the layers and the density of the fluid in the middle part of the

(a) (b)

(c) Fig. 2. Snapshots of simulated system projected in the XZ plane for: (a) Z = 2.0or, (b) Z = 4.0or; (c) Z = 7.0a. Smaller dots represent the centre of fluid atoms and larger ones, at the top and bottom, the wall atoms.

A. Jabbarzadeh et al./ J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

177

7

?--0.8

0 "~

0.0

5

-0.5

>

d,~'- ...-/,

-~

?=,0 ?=2.0

----

Linear Velocity

~0.5

2"~..,J .,¢5.~5Jo-

?=3.0

.,,4[.~-~

° .,~-

5

o

3

r~ ¢~ -1.0 0

> -1.5

-1

o

1

Z 1.o

. ,

. . . .

~

. . . .

r

. . . .

i

. . . .

i

. . . .

i

. . . .

i .

Z=7~

0

~

~

o.o

~

0 .2

-o.5

/J.~r •~ .~.o,

0

'f

J~'*" .

~

~

+

~

,~

/~

-2.0 ~

? =0.45

•

:~--065 "?=1.5

tm

"[=2.0

~,=3.o .

-3

,

.

.

.

-2

.

.

.

.

.

0

-1

'

'

'

1

'

'

'

1

~2

B

t

'

, m ~ 6 ~ O J

2

0

3

Z

......

,7,=,0

.........

. • 0.5 O O

~

0,0

5-0.5

.g,

8 _10

--~---T =0.8 +7=1.0 ---o--7=15 ~ 7 =2.0 ? =3.0

~..

>

r0~

a ,~

-1 5

-2.0 -4

-2

0

2

4

Z Fig. 3. D e n s i t y and velocity profiles for various shear rates for wall separations o f 4~r, 7 a and 10a. All graphs are for a wall strength parameter o f ew = 1.0e.

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

178

20 f

. . . .

i

. . . .

i

. . . .

i

15

~I0

•

Z=4o"

~,

Z-7O

//

"6

a5

0

1

2

3

Shear Rate ('~)

Fig. 4. Degree of slip against the shear rate for various wall separations; for all three cases the wall separations ew = 1.

slit. Increasing the shear rate makes layers slightly thinner and their density profile steeper close to the walls, and leaves the density in the middle region slightly higher. This can be explained by looking at the velocity profiles. These profiles are third-order polynomials fitted through the average velocity of slices as described in Section 3.1. The solid line again represents the perfect linear velocity described by hydrodynamics theory. For lower shear rates the velocity profile remains either linear or slightly curved along the perfect linear line. As we increase the shear rate the profiles begin to deviate from linearity and their slopes decrease. If we further increase the shear rate the slope of the profiles eventually decreases to a point where the velocity is zero compared to that of the wall. The onset of this deviation and loss of slope depends on the wall separation. Also the degree of shear that the system can tolerate and still maintain the flow depends on the film thickness or, in other words, the wall separation. For a film of 10o thickness the velocity profile remains almost linear up to a shear rate of 0.8(g,/mo2) 1/2. At shear rates of 1.O(e/m02) l/2 and higher the profile is curved with higher slope in the vicinity of the wall and lower slope at the middle of the slit. The onset of this change for Z = 70 is at a higher shear rate of 1.0(6/mo2) 1/2 and for Z = 40 it happens at 7 = 1.0(,o,,/mo2)1/2. The average slope of the velocity curve can be regarded as a measure of the degree of slip or non-slip condition. A slope lower than the slope of the perfect linear solid line is an indication of slip and the lower this slope the higher is the degree of slip. We define the degree of slip D~ by Ds = '-- - 1,

(8)

where ~'s is the slope of the average linear velocity fitted through velocity points. A positive value for Ds represents a slip condition and zero or negative values indicate a non-slip condition. The absolute value of Ds is an indication of how strong the slip or non-slip condition is for the flow. In Fig. 4 we have plotted Ds against the shear rate. It can be seen from Fig. 4 that the non-slip boundary condition exists at only small shear rates for all three wall separations. As we increase

A. Jabbarzadeh et al./'J. Non-Newtonian Huid Mech. 69 (1997) 169-193

1.0

179

lO

0.5

8 ".,~ 0.0 ~-0.5 -10 J/.I

--0.5 e~ =2.0 e~=3.o

.

~,o

/111 -- ,. I[/^l

~-15

,9,

Illll

" ~ -2.0

;

E: =9.0

--

>

LinearVelocity

2

-2.5 -30 -2

-1

0

0

Z 1.0 , , q . . . .

I

.

.

.

.

.

.

.

.

I

.

.

.

.

.

.

.

1 1 '

>-o~:

10

~ ~,

~

'g

.5

• ~ =1.0 -1.5 L

rk

I~I

• -2,0

~

*

~ --4.0

~; =5.0

.

- • g~=2.0

4

~~

20

-25 -30

-3

-2

-1

0

1

2

Z

1~o

'

'

I

.

.

.

.

r ....

0.5

i ....

I ,

,

y

"~'~

0.0

~

~0

8

-0.5

8 ~

-2.0 -2.5 -3.O ~

'

-2

0

2

4

Z Fig. 5. Density and velocity profiles for various wall-fluid interaction strength parameters for wall separations of 46, 7a and 10G. In all cases the shear rate was 0.2(t'/mo-2) t'2.

180

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

the shear rate the slip boundary condition begins to appear and to grow rapidly at higher shear rates. Fig. 4 confirms our discussion regarding the onset o f slip for different film thicknesses, based on the analysis of the velocity profiles. It can be seen that the onset of the slip begins at higher shear rates for the thinner films. The a m o u n t of slip at the same shear rate is always less for thinner films. In the velocity profiles of Fig. 3 we can see that increasing the shear rate leads to a point where the state of pure slip is observed over the whole slit and a slug-flow condition is stabilised. We call the shear at which this state happens the non-flow shear. For Z = 100. non-flow shear is about 2.0(elm0.2) ~/2, for Z = 70- it is about 3.0(elm0.2) 1/2, while for Z = 4.00. non-flow shear is much higher than 3.0(elm0.2) ~/2. So we can conclude that increasing the shear rate enhances the slip boundary effect, while reducing the wall separation reduces the slip and makes it possible to maintain the simple ideal flow at higher shear rates.

3.1.3. The effect o f wall-fluid interaction strength The importance of the wall strength parameter (ew) is revealed when we study its effect on the density and velocity profiles. Again the results for three different film thicknesses of 100., 70- and 4o- are presented in Fig. 5 for various values o f ew ranging from 0.5e to 5.0e. For all the runs the shear rate was 0.2 (elm0.2) 1/2. F r o m Fig. 5 it can be seen that increasing ew from 0.5e to 5.0e results in some interesting consequences. The effect on density profiles in all three systems is quite pronounced. As we increase ew the density profiles become sharper and the fluid layers shift closer towards the walls. The other interesting phenomenon which has not been reported before is that increasing ew further leads to a secondary depletion zone just after the first fluid layer near the wall. This is an indication that the first layer has solidified and stuck to the wall. This may result in a situation where even more layers become solidified. As mentioned we think this phenomenon happens because the first layer is stuck to the wall and has a solid state in the form of a crystal layer. This can be better imagined if one looks at the snapshots of the simulation box taken for a fluid slab of 100. thickness for two values of e~v= 1.0e and 10.0e. These snapshots are presented in Fig. 6. F r o m Fig. 6 it can be seen that for a moderate

(a)

(b)

Fig. 6. Snapshots taken from the simulation box in the XZ plane for a film of lO~rthickness under Couette shear flow for two values of the wall-fluid interaction strength parameter: (a) ew= 1.0e; (b) e~v= 10.0e. For both cases the shear rate w a s 0.2(c,/m~r2) I/2.

A. Jabbarzadeh et a l . / J . Non-Newtonian Fluid Mech, 69 (1997) 169-193

181

1.0

0.5

•

Z=4c

A

Z=7o

•

Z=IOo

¢/)

•

0.o

o a -0.5

-1.0

. . . .

i 1

. . . .

i 2

. . . .

[ 3

. . . .

i 4

. . . .

i 5

. . . .

i

. . . .

6

Wall Strength (Ew) Fig. 7. Degree of slip against the wall strength parameter ew for various wall separations of 4or, 7or and 10a. For all the simulations shear rate was 0.2(~/mo-2) 1:2.

interaction between the wall and fluid particles at e.w= e there is a layer of fluid near the wall and the rest of the fluid in the middle region has a fairly homogeneous structure. However increasing the wall strength tenfold leads to solidification of the first layer near the walls in a crystallised structure copied from the structure of the wall atoms. The second layer is also very distinct but it is not in a solid state. These two concentrated layers leave the middle part of the fluid strongly depleted and with a low average density. The way velocity profiles are affected by changing ew is highly correlated with the density profiles and this correlation is obvious in the profiles of Fig. 5. In the secondary depletion zone the average velocity data are not statistically reliable. These points are not displayed in the graph. The solid lines for the velocity profiles are the third-order polynomial fits of the velocities. In the cases where the secondary depletion zone existed, the curves are fitted only through the velocities in the middle region. Velocity profiles also confirm that the first layer sticks to the wall when e.w reaches about 4.0e for a film thickness of 10~r. This happens for film thicknesses of 7o- and 4a at slightly lower values of ew. While velocity profiles in the middle region of the slit are almost linear at low and moderate values of ew up to 1.0e, increasing it further makes the velocity profiles become curved, as happens when we increased the shear. However there is an important difference in this case. The average slope of the velocity curves is higher than the slope of the perfect linear velocity, which means that increasing wall-fluid interaction reduces the slip boundary condition to the point of no-slip, as can be seen. We used Eq. (8) to evaluated the degree of slip at different wall strengths. In some cases, because the velocity curves at some high values of ~ were distorted, the average slopes of the velocity curves were evaluated in the middle region. Dependence of the degree of slip D~ is illustrated in Fig. 7. It can be seen that for most cases there is a non-slip condition for all three wall separations, corresponding to negative values of Ds. The slip condition is present only at low values of ew. As we increase e.w,the degree of slip decreases further which is an indication that the no-slip condition provides a locking state between the fluid particles and the wall through the solidification of the fluid particles.

182

A . J a b b a r z a d e h et a l . / J .

Non-Newtonian

Fluid M e c h . 69 (1997) 1 6 9 - 1 9 3

Nevertheless, D, appears to reach a limiting value, depending on the wall separation, for wall strengths greater than 2-3e.

4. Stress tensor and material functions

4.1. Stress tensor Stress tensor components were found for microscopic system of particles by the Irving-Kirkwood [16] method. According to this method the contribution of each particle to the stress tensor is in two parts; a configuration part and kinetic part. This can be written as (7~/3 ~- - -

milgi~li[3 q-

. . ~ . rii=F(i/l t

(9)

.

,l>!

The first sum in the right hand side of Eq. (9) denotes the kinetic contribution where mi is the atomic mass and ~ and fl are coordination system axes which for a Cartesian system can be simply substituted by X, Y or Z, and ui~ and u~/~ are the peculiar velocity components of particle i in the c~ and fl directions. The second sum represents the configurational or potential contribution where r~/~ is the c~ component of the distance vector between particles i and j and F~//~ is the fl component of the force exerted on particle i by particle j. We have to exclude the mean flow velocity when we consider the laboratory velocity component of a particle in the flow direction. Then for shear stress we can rewrite Eq. (9) as ~Tzx - -

mibliz(Uix -- U,-j) +

-~

r~/:F,-/x

(10)

i j> i

where U.,.j is the average flow velocity at the position of particle i. The angle brackets denote the time average. As a test of accuracy of Eq. (10) for calculation of shear stress, we computed the time average of the force in the X direction applied to the wall particles by fluid particles during the simulations. Then we calculated the shear stress by dividing this force by the area of the walls. This shear stress is given by Eq. (11): Nw NF

~rxz.w = Z Z r
(11)

j

Table 3 shows the shear stresses computed by two different methods and it can be seen that there is a good agreement between them except at shear rate 1. Table 3 Shear stresses computed from Eq. (10) and (11) for a slit of 7cr thickness

a,__ a ......

0.05

0.1

0.45

1.0

2.0

3.5

0.0543 0.0546

0.1006 0.1007

0.4704 0.4607

0.8791 0.9865

0.5171 0.5451

0.1803 0.1812

A. Jabbarzadeh et al.//J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

183

We can see the same agreement in Ref. [5]. We adopted Eq. (10) for the calculation of the shear stress for the rest of the simulations. 4.2. Viscosity and material functions The shear rate is slightly different at different distances from the wall because the gradient of the velocity is not constant along the Z axis; therefore the local viscosity will be slightly different. More useful is the average over the whole width of the slit. Thus we use the following constitutive equation to find the viscosity r / - o.=" ?;

(12)

For a rheologist, normal stress differences are interesting and important properties to measure and examine. The first and second normal stress differences are N: = a,..,. - o-=_.,

N2 = o-:: - o->.:..

(13)

We use the normal sign convention of mechanics [14], so that a positive sign represents tension. N o r m a l stress differences for Newtonian fluids are zero but non-Newtonian fluids exhibit normal stress difference effects. We can also define the normal stress coefficients as o-xx - - o-~~1

5,2 -- S,

o-:: - - o-yy

~b2 --

72

,

(14)

where ~,: and ~'2 are first and second normal stress coefficients, sometimes referred to as the viscometric functions [17]. For Newtonian fluids the normal stress coefficients are zero whereas for polymeric fluids we expect a positive value for ~,~ and a negative value for ~'2 with a ratio of l~,/~'2l of about 10 [17]. In Section 3.1 we examined the effect of many parameters on the structural properties of the sheared fluid. In this section we will examine the material properties of the fluid. We seek a correlation between what we found previously for density and velocity profiles and what we find now for viscosity and normal stress differences. 4.2.1. The effect o f the system size To ensure that the results from different system sizes can be compared we check whether system size affects fluid properties. We characterise the system size by the n u m b e r of particles in the system. In our case this is the number of fluid particles plus the number of wall particles. We performed the simulation for various films with different thicknesses and analysed the results for density and velocity profiles and also viscosity. In Table 4 we present the results obtained for viscosity for three different system sizes at three wall separations of 2o., 4o- and 7o-. It can be clearly seen that the results for viscosity change only very slightly. We have the largest difference of about 4.90% for Z = 2o- when we increase the system size from 249 atoms to 895 atoms. We performed several simulations for the same system size of 249 atoms for Z = 2o- and calculated the standard deviation in the computed viscosity. This standard deviation was 2.87%. We see that the difference in the viscosity for different system sizes is only slightly more than the statistical error in this case. The system-size-induced difference in the case of Z = 4o- is even less

184

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

Table 4 Computed viscosity for film thicknesses of 2or, 4or and 7ct obtained at three different systems sizes Wall separation (0)

Viscosity

Total no. of atoms

2

249 610 895

0.5963 0.6061 0.6255

4

338 827 1215

1.6113 1.5794 1.5841

7

471 1154 1695

1.0746 1.0718 1.0800

(mg)l/2/a 2

and is limited to 1.72% when we increase the number of total atoms 4 times from 338 to 1215. We think this difference can be considered as negligibly small, of the same order as the statistical fluctuations due to the finite run length. We also considered velocity and density profiles to ensure that there is no big difference in the results obtained for different system sizes. Fig. 8 shows the velocity and density profiles for a wall separation of 4o-. It can be easily seen that the density profiles are identical for systems simulated with 178 and 435 fluid atoms respectively. However the velocity profiles show a slight difference which could, however, be attributed to statistical fluctuations. Nevertheless we think that in the range of accuracy that we sought the results for all system sizes are almost the same. We conclude that the results for different system sizes can be compared qualitatively and quantitatively.

1.0

,

,mmmmlm • 0'5I ,~.

• •

N f = 1 7 8 & N =160 Nf=435 & N,~=392

3.5 3.0

ill mill

0.0

4.0

2.5

L) 0 " ~ -0.5

>,

2.0 ~

>

1,5 -1.0

1.0 -1.5

0.5 -2.0 -2

::::::J

. - - 0.0 -1

0

z

(o)

1

Fig. 8. Velocity and density profiles for a 4a thick film at different system sizes.

A. Jabbarzadeh et al./'J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

-

185

Shear Viscosity Of Confined Fluid Shear Viscosity Of Bulk Fluid

-

"~ 3 O .u) 2 > 1

-i

'

'

'

'

51

,

,

,

,

101

,

,

,

,

151

,

,

,

, 210

,

,

,

,

Wall Separation Z (~) Fig. 9. Effect of film thickness on the viscosity of the sheared fluid. The horizontal line shows the viscosity of the bulk fluid obtained from the NEMD algorithm of SLLOD.

4.2.2. The effect of film thickness on viscosity F r o m what we saw regarding the effect of the film thickness on the velocity and density profiles, we expect to see a p r o f o u n d effect on the material functions. F r o m velocity profiles obtained from Eq. (9) the elements of the stress tensor were obtained. Then from Eq. (12) we obtained the viscosity. Simulations were conducted for a broad range of different film thicknesses ranging from 20- to 200.. Also, in order to examine the bulk fluid viscosity, we ran the S L L O D program for the same fluid at the same density and temperature conditions. We used Eq. (12) to compute the viscosity for bulk fluid. The results have been plotted in Fig. 9, from which it can be seen that for films 2 . 5 - 4 molecular diameters thick the fluid viscosity is several time the viscosity of the bulk fluid. The reason for this dramatic increase is not k n o w n exactly. However some believe the layering and confinement-induced glassy transition [9] are responsible for this event. In this range of wall separation we observe an interesting p h e n o m e n o n . As we decrease the wall separation, viscosity increases to the point where the wall separation is an exact integral n u m b e r of molecular diameter and then viscosity falls. For smaller wall separations this fluctuation in the viscosity is more pronounced. It seems that higher viscosities result when fluid layers have difficulty in passing each other. For example if we have a look at the snapshots in Fig. 2 we can see that for a wall separation of 2a there is only one layer formed in the space between the walls. Since there is only one layer competition between layers does not exist and there is plenty of r o o m for sliding between this layer and the fixed wall layers. This is the main reason for the sharp drop in the viscosity in this case. If one increase the wall separation to 2.50., it can be seen from the density profiles in Fig. 1 that two layers form between the walls. These layers can only just fit between the walls and shear flow requires a high interaction between the layers passing each other which results in higher viscosity. For a wall separation of 30- there are also two layers between the walls but they are not as close together and it is easier for the layers to pass each other and hence lower viscosity is observed in this case. The film viscosity approaches the bulk fluid viscosity as we increase the film thickness. F o r films thicker than 10o~ the viscosity is about equal to the bulk fluid viscosity.

186

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193 1.0

. . . .

~

. . . .

i

. . . .

i

. . . .

0.5

r

. . . .

• .1.

Z=4m Z=7a,

•

Z=lOc,

i

¢ =£ Ew=

g

{:~=~ B u l k Fluid

• 0.0

~.--,o0.5 v {3') 0 -1.0

-1.5

-2.0

.... -2.5 -1.5

i

....

-1.0

i

....

-0.5

i

,

,

,

0.0.

,

i

. . . .

i

,

1.0

0.5

Iog(, ) Fig. 10. Logarithmic plot of film viscosity against the shear rate for three films of different thickness. The relationship between shear and viscosity is also displayed for the bulk fluid.

4.2.3. Viscosity and shear rate The viscosity was obtained for a number of shear rates between 0.05 and lO.O(e/m~r2) 1/2 at the same wall-fluid interaction strength parameter of ew = 1.0e. The simulations were repeated for three wall to wall distances of 10a, 7o- and 4o-. We also computed the shear viscosity from the SLLOD algorithm for the bulk fluid over the shear range O.05-7.0(e/m62) 1/2. The results for all four systems are displayed in Fig. 10. From Fig. 10 it can be seen that for all three film thicknesses we have two regions, one where the fluid is Newtonian and the other where it is non-Newtonian. The non-Newtonian regime exhibits a shear-thinning effect and viscosity obeys a power law. Fig. 11 shows the power-law region for the same systems. Solid lines are linear fits of the l o g - l o g plot of viscosities. 1.0

. . . .

4

. . . .

n

. . . .

0.5

I

. . . .

•

n

Z:lOo,

. . . .

I;w=¢

O.C

'~0.5

o -1.0

-1.5

-2.0 0.0

0.2

0.4

0.6

08

1.0

log(y) Fig. 11. Viscosity vs. Shear rate in the power-law region.

A. Jabbarzadeh et a l . / J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

187

Table 5 Power-law region coefficients for three different wall separations presented in Fig. 11.

z (0)

A

4.0 7.0 10.0

1.816 0.943 0.563

- 1.943 - 2.129 -2.111

F r o m Fig. 1 1 it can be seen that the viscosity obeys a power law of the form q = A/~.

(15)

Values of A and ~ for three different wall separations are given in Table 5. The value for A in Eq. (15) drops as we increase the film thickness, while ~ is about the same for wall separations of 4o- and 10(7. For Z = 7o- the slope is a little higher. An interesting fact is that for the bulk fluid shear thinning begins at much higher shear rates of about 4-5(e/mo.2) ~/2. For lower shear rates the behaviour of bulk fluid remain Newtonian. For the non-Newtonian behaviour observed in the presence of the walls it seems that the onset of shear thinning depends on the wall separation. We can see that shear thinning begins at higher shear rates for films 4~r and 7o- thick compared to the thicker film of 10o.. Considering the values of c( in Table 5 one can see that the slope o f the viscosity lines in the p o w e r - l a w region is higher than what we usually expect for non-Newtonian fluids (~ > - 1). This implies that the shear stress decreases as we increase the shear rate. This decrease in the shear stress can be seen in Fig. 12. We can see the same behaviour in the shear stress for bulk fluid at higher shear rates. Heyes [18] has performed simulations for a L e n n a r d - J o n e s fluid using the S L L O D algorithm and a profile-biased (assuming linear streaming velocity profile) Gaussian thermostat. He has 1,0 Am 0.5

Z = 4 C , ~:.= E Z = 7 o , (:w = • Z=lOo, ~ = ¢

•

...--.. 0.0

v O')-o.s 0 -1.0

-1.5 !

-2.0

, -1.5

+

~

~

r -I.0

,

,

+

,

I

,

,

,

,

-0.5

I 0.0

,

,

,

,

I 0.5

,

~

~

, 1.0

log(y) Fig. 12. Shear stress against shear rate for various wall separations and bulk fluid. For all simulations in the presence of the walls t~, = 1. Bulk fluid results are obtained from the SLLOD algorithm.

188

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193 20

•

Z=I Oo

•

Z=7o

1.5

o 0

0.5

0.0

I

2

a

i

ew

i

a

I

,

Fig. 13. Dependence of viscosity o n ew for three wall to wall distance of 4.0or, 7.0~r a n d 10.0o.

observed a sharp decrease in the shear viscosity proportional to :~-1.5. He also observed a decrease in shear stress beyond a certain shear rate and attributed these effects to the 'string phases' developed at these shear rates. This behaviour can be an indication o f instability but we do not know whether this conclusion has any physical validity, or whether it is an artifact of the particular model we have adopted. In particular, it is hard to see why the point of transition to the power-law region does not approach that o f the SLLOD (bulk fluid) simulation as film thickness Z increases, but instead moves in the opposite direction.

4.2.4. The effect of wall-fluid interaction strength on viscosity The effect of wall strength parameter ~ has been considered in some previous works [8,19]. We think a more detailed insight into the effect of this parameter is important. Fig. 13 shows how the viscosity of the sheared film changes for three wall to wall distances of 4o., 7o- and 10oover a range 0.01e-5.0e for ew. It can be seen that the viscosity increases for 0.01~ < ~ < 1.0e. for all three film thicknesses. There is a turning point at e.w= 1.0e where all three films reach their peak viscosities. Increasing the strength of the interaction further results in a dramatic drop in viscosity for all three cases. This decrease in viscosity can be compared to the shear-thinning effect that was observed before. We will call this effect wall-induced thinning. The pattern of this however is not the same for different wall separations. The solid lines in Fig. 13 are the third-order polynomial fits of the viscosities in the range 0.01 < ew < 5.0~. There is a point at which all the curves meet. A T e.w= 3~, regardless of the film thickness, the viscosities become almost identical. This point may be an important point since viscosity loses its dependence on the film thickness. In Section 3.1 we came to the conclusion that thinner films exhibit higher viscosities. It seems that for ~ > 3~ this conclusion is no longer valid. After this point we get higher viscosities for Z = 10o- than for Z = 4o-. This viscosity drop effect may be attributed to the solidification of the fluid particles close to the walls which leaves the central region of the slit with lower density. The viscosity decreases until most particles solidify in crystal structures close to the walls and then it remains at about

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

189

the same value. For a thin film of 4~r this happens at lower values for ~ because only two layers have to be solidified near each wall. For thicker films the number of layers increases and because of the short range forces involved in the interactions it becomes very difficult for the wall to solidify all the fluid particles. As we saw before in Fig. 6 for a film 10o- thick even for ~,, = 10.0e only two layers are solidified at each wall.

4.3. Normal stress differences 4.3.1. Shear rate and normal stress differences We examined the normal stress differences for three typical films at shear rates between 0.05 and lO(e/m~r2) ~/2. Our findings are shown in Fig. 14. The results suggest that in all films studied, for shear rates lower than about 0.2(g/m~r2) 1/2, both NI and N 2 remain at low values close to zero. This again confirms the existence of a Newtonian regime at lower shear rates in the same region as we saw before for the viscosity. The first normal stress difference (NO is positive for all three films. It can also be seen that thinner films exhibit higher magnitudes for both stress differences. For shear rates greater than l.O(~/m~r2) n/2 for Z = 4a, N~ is several times its value for Z = 10o-. We find that the second normal stress difference N2 is always negative. We also find that N2 is almost equal in magnitude to N~ with the opposite sign. This may be an artifact of the periodic boundary conditions used. 4.3.2. Wall strength and normal stress differences Fig. 15 shows the effect of the wall strength parameter ew on the normal stress differences. Although the normal stress components in the stress tensor can have any value according to the principles of continuum mechanics, it is expected that Nj be positive [14,20]. It is interesting to see that increasing the wall strength parameter changes the normal stress differences in such a way that N] shifts to negative values and N2 becomes positive. This sort of behaviour only can be expected in strongly oriented structures as in liquid crystals [20]. It is also seen from Fig. 15 2.0

.

.

.

.

i

.

.

.

.

i

..

..

..

..

it

..

..

..

.

.

ii

,.

.

.

,

,,-J

.

.....

Z ¢9 (~)1.5

II) '4"-°B

a

1.0

(/) U)

~0.5

E o Z

~.

0.0

N~, Z=4cf N 1 , Z=7~ N~ , Z = 1 0 , ~

•

LL -0.5

,

i

i

I 2

,

i

,

i

I 4

. . . .

~ 6

. . . .

i

8

.

.

.

.

10

Shear Rate ~, Fig. 14. First normal stress difference against shear rate.

190

A. Jabbarzadeh et a l , / J. Non-Newtonian Fluid Mech. 69 (1997) 1 6 9 - 1 9 3

2

. . . .

i

. . . .

i

'

' N~, Z=4.0o N~, Z=7.0o N~, Z = I O . O ~

¢- 0 ~

o

~.6

o Z

LL

i

,

,

l

l

l

l

,

l

2

l

l 4

l

l

,

l

,

,

,

,

6

gw

Fig. 15. Dependence of first normal stress difference on ew for various films. that N~ decreases and N2 increases as e~v increases. Negative N1 can be an indication that the fluid under shear tends to shrink inwards and we have to exert extra force to keep the walls from moving towards each other. This can be imagined knowing that at high values of e~vmost of the fluid particles are attracted to the walls, leaving the central region, which can be considered as a microscopic vacuum, at low density. Thus the fluid particles tend to shrink inwards and fill this vacuum. However, the simulation box and hence the wall separation are forced to be constant, so is the number of particles. Therefore this vacuum can not be filled and there is no way out but to exert a force pulling the walls inwards which results in a negative value for N~. This phenomenon thus appears to follow the assumptions of the simulation and may not occur in practice. 4.4. Shear rate and first and second normal stress coejficients Some rheologists prefer to use normal stress coefficients ~'1 and ~'2 instead of NI and N2 to study the normal stress differences. We used Eq. (14) to compute the time average values for Oj and ~'2. Fig. 16 shows the logarithmic plot of ~,1 against shear rate. The solid lines are the third-order polynomial fits to the data. The shape of the curves for thicker films of 7o- and 10ois close to linear. Nevertheless, it becomes more strongly curved for a thinner film of thickness 4o-, corresponding to the maximum in N~. The sign of g'l remains positive for the whole range of shear rates and film thicknesses displayed in Fig. 16.

5. The effectiveness of the Gaussian thermostat

In the case of microscopically confined fluid which undergoes shear flow, as shown before, the flow velocity is not known and to remove the generated heat some researchers [15] rescale the velocities only in the nonflow directions (Y and Z). Although this method might work well at low shear rates we will show that it is not effective at higher shear rates. One of the important

A. Jabbarzadeh et a l . / J . Non-Newtonian Fluid Mech. 69 (1997) 169-193 2.0

'

,

.

.

.

.

i

.

.

.

.

i

191

r

15 •

\'~

•

Z = 7 o ,

•

Z = 4(;

10

E.=E

,

~w = E

I

,

05

o -0,5

-10

-15

h

-20

,

,

,

-2

I

h

,

,

,

-I

P

r

,

,

h

0

,

,

,

I

Iog(~,) Fig. 16. First normal stress coefficient against shear rate ['or three different wall separations.

features of this work has been the use of a Gaussian thermostat which was introduced by Eq. (4). In this section we will show how effective this method is in comparison to the mentioned velocity rescaling method. In order to make a judgement between the two methods we did the simulation for a film of 4o- thickness at a constant temperature of 1.2(e/kB) over a range of shear rates between 0.05 and 5.0(~:/mo-2) '/2. We performed the simulations once by using a Gaussian thermostat and velocity rescaling and once more by using only the rescaling method. We obtained time averages for temperature and viscosity and plotted them against the shear rate logarithm. The results are shown in Fig. 17. The smaller graph in Fig. 17 is the same graph zoomed to show the results in the range 0.05-1.0 more clearly. The results reveal that at lower shear rates up to 0.2 both methods work excellently in terms of fixing the temperature. The results for the viscosity for both methods remain almost identical. At shear rates between 0.2

14

12

0

-1

E

'~'~-

....

~)

....

,~6' °

"e

6

-2

I--

- - - -o- - - - T e m p e r a t u r e - - - -,- - - - V i s c o s i t y

With

With

Temperature Viscosity

i

0 -2

i

,

Only With

With

i

I

I

Only

Velocity

Velocity Gaussian

Gaussian

.

.

i

log('"f )

Rescaling

Rescaling Thermostat

Thermostat

,

I

-3

P

,

,

h

'

-4

o

Fig. 17. Temperature and viscosity against shear rate for rescaling and thermostatting methods. The thickness is 4a.

192

A. Jabbarzadeh et al./J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

and 0.65 the temperature begins to deviate from its constant value for the rescaling method, although this deviation is still low. For this method in this region the maximum percentage deviation is about 10% at shear rate of 0.65. In the same range of shear rate, there is only a very slight deviation in temperature for the thermostat method, at most 0.11%. Increasing the shear rate further in the range 0.65-3.0 has catastrophic consequences for the rescaling method. It is in this range that the thermostatting method demonstrates its superiority over the rescaling method. In this range there is a sudden j u m p in the temperature for the rescaling method and rescaling can not cope with the heat generated by shear flow. Eventually the temperature rises to about 4 times the required temperature when the shear rate reaches 5.0. At a shear rate of 5.0 the deviation from fixed temperature for the thermostatting method is only about 1% which can be considered as excellent. The poor performance of the rescaling method at higher shear rates is perhaps attributable to the fact that the velocity in the direction of the flow is not rescaled. However it is possible to do a full velocity rescaling provided the flow velocity is subtracted from the component of the velocity in flow direction. This possibility has not been examined in this work. Since viscosities are almost identical in the region where both methods still keep the temperature fixed, we can conclude that our method works well and is reliable in terms of using the thermostat. Somers and Davis [15] have computed the viscosity for a film of 4o- thickness for the same fluid density and temperature. Their computed value is about 1.57-1.59(me.)~/2/o- 2. In that work they employed the reservoir method [21] to extract the heat generated by shear flow and to keep the temperature of the system constant. Our result for viscosity in Table 4 is in excellent agreement with their results. This agreement is another confirmation of the effectiveness of our method. Our method, in comparison to the reservoir method, has the advantage that we do not need to the extra number of atoms for the reservoir, therefore reducing the required number of atoms for the system and increasing the computational efficiency of the simulation. We can also conclude that the results at higher shear rates are still reliable for the thermostatting method making it applicable to a larger range of shear rates.

6. Summary In this paper we presented the molecular dynamics simulations of a thin fluid film sheared in Couette flow geometry between two structured plane walls. An N V T ensemble was chosen and simulation conducted in isothermal conditions. To keep the temperature at the required level a Gaussian thermostat was employed. The superiority of this method was demonstrated over the simple velocity rescaling method, especially at high shear rates. It was also demonstrated that the final results for this method are compatible with the previous results of some other researchers. The results for density and velocity profiles were obtained for a wide range of simulation parameters. The effects of shear rate and wall-fluid interaction strength were investigated in detail for a wide range of parameters. The material functions were also obtained and the effects of shear rate and wall strength parameter on these properties were studied. The effect of film thickness on the viscosity was investigated and compared with what we found for bulk fluid from the SLLOD algorithm. The results show an increase in viscosity for thinner films. This is in agreement with previous work, but we observed that this was valid only for a specific range of ew. The existence of a non-Newtonian region with shear thinning is also found

A. Jabbarzadeh et al./ J. Non-Newtonian Fluid Mech. 69 (1997) 169-193

193

a n d e x a m i n e d for v a r i o u s films. H o w e v e r we c o u l d n o t rule out the possibility o f instability because o f the high slope o f the v i s c o s i t y - s h e a r rate curve in the n o n - N e w t o n i a n region. A viscosity r e d u c t i o n by increasing the wall strength was also observed. W e called this phen o m e n o n the wall-induced-thinning effect and a t t r i b u t e d it to the solidification o f fluid particles at the walls. N o r m a l stress differences were f o u n d for sheared fluid in the presence o f the walls. It was s h o w n that the n o r m a l stress differences are highly d e p e n d e n t on the w a l l - f l u i d i n t e r a c t i o n strength and that they are higher for t h i n n e r films.

Acknowledgements W e gratefully a c k n o w l e d g e the s u p p o r t o f an A u s t r a l i a n R e s e a r c h S c h e m e for this work. W e also wish to t h a n k P r o f e s s o r D.J. E v a n s for p r o v i d i n g the S L L O D algorithm; P r o f e s s o r N. P h a n - T h i e n for helpful and c o n s t r u c t i v e discussions a n d Dr. Peter Daivis for a critical reading o f the m a n u s c r i p t .

References [1] J.N. Israelachvili, P.M. McGuiggan and A.M. Homola, Science, 240 (1988) 190. [2] J.N. Israelachvili, Colloid Polym. Sci., 264 (1986) 1060 1065. [3] J.N. Israelachvili, J. Colloid Interface Sci., 44 (1985) 263. [4] D.Y. Chan and R.G. Horn, J. Chem. Phys., 83 (1985) 5311. [5] 1. Bitsanis, J.J. Magda, M. Tirrell and H.T. Davis, J. Chem. Phys., 87 (1987) 1733. [6] I. Bitsanis, S.A. Somers, H.T. Davis and M. Tirrell, J. Chem. Phys., 93 (1990) 3427. [7] I. Bitsanis, T.K. Vanderlick, M. Tirrell and H.T. Davis, J. Chem. Phys., 89 (1988) 3152. [8] E. Manias, G. Hadziioannou, I. Bitsanis and G.T. Brinke, Europhys. Lett., 24 (1993) 99--104. [9] P.A. Thompson and G.S. Grest, Phys. Rev. Lett., 68 (1992) 3448. [10] P.A. Thompson and M.O. Robbins, Phys. Rev. A, 41 (1990) 6830. [11] S.Y. Liem, D. Brown and H.R. Clarke, Phys. Rev. A, 45 (1992) 3706. [12] R. Edberg, D.J. Evans, and G.P. Morriss, J. Chem. Phys., 84 (1986) 6933. [13] D.J. Evans and G.P. Morriss, Statistical Mechanics of Nonequilibrium Liquids, Academic Press, 1990. [14] R.I. Tanner, Engineering Rhelogy, Oxford Science, 1985. [15] S.A. Somers and H.T. Davis, J. Chem. Phys., 96 (1992) 5389. [16] J.H. Irving and J.G. Kirkwood, J. Chem. Phys., 18 (1950) 817. [17] R.B. Bird, R.C. Armstrong and O. Hassager, Dynamics of Polymeric Liquids, Vol. 1, Wiley-lnterscience, 1987. [18] D.M. Heyes, J. Chem. Soc., Faraday Trans. 2, 82 (1986) 1365-1383. [19] U. Heinbuch and J. Fischer, Phys. Rev. A, 40 (1989) 1144. [20] H.A. Barnes, J.F. Hutton and K. Walters, An Introduction to Rheology, Elsevier, 1989. [21] T. Ashurst and W.G. Hoover, Phys. Rev. A, 11 (1975) 658.