Findind Pxy for calculating Viscosity in a couette flow

LAMMPS LAMMPS Mailing List Mirror
1

I want to calculate viscosity in a couette flow, for which I require Pxy.
From what I read in the manual Lammps does that for us. But its my understanding
that it calculates that for an instant of time. Hence I require to average over time.

I am using the ave/correlate to average it:

variable pxy equal pxy
.
fix 6 all ave/correlate 1 50000 50000 v_pxy type auto file pres_fb.dat
.
thermo_style custom step temp etotal c_mobile_thermo v_pxy

but I am getting the following error:

ERROR: Thermo keyword in variable requires thermo to use/init press

when checked in manual it says:
"Thermo keyword in variable requires thermo to use/init press
You are using a thermo keyword in a variable that requires pressure to be calculated, but your thermo
output does not use it. Add it to your thermo output"

I am unable to figure out what is it that I am doing wrong.

Also find at the end the complete code.

Thank you.

BR/Joseph

# 3-d LJ flow simulation

units lj
dimension 3
boundary p s p

atom_style atomic
neighbor 0.3 bin
neigh_modify delay 1

# create geometry

region fluid block 0 11.95 1.5 26.07 0.0 7.22 units box
region btm_wall block 0 11.95 0 1.5 0.0 7.22 units box
region top_wall block 0 11.95 26.07 28.0 0.0 7.22 units box
region simulation_box union 3 btm_wall fluid top_wall units box

create_box 3 simulation_box

lattice fcc 0.81 orient x 1 1 -2 orient y 1 1 1 orient z 1 -1 0
create_atoms 1 region fluid units box
lattice fcc 0.81 orient x -1 1 0 orient y 1 1 1 orient z 1 1 -2
create_atoms 2 region btm_wall units box
lattice fcc 0.81 orient x -1 1 0 orient y 1 1 1 orient z 1 1 -2
create_atoms 3 region top_wall units box

mass 1 1.0
mass 2 1.0
mass 3 1.0

# LJ potentials

pair_style lj/cut 2.2
pair_coeff * * 1.0 1.0
pair_coeff 1 2*3 0.1 1.0
pair_modify shift yes

# define groups

region 1 block INF INF INF 1.5 INF INF units box
group lower region 1
region 2 block INF INF 26.07 INF INF INF units box
group upper region 2
group boundary union lower upper
group flow subtract all boundary

set group lower type 2
set group upper type 3

# initial velocities

#compute mobile flow temp
#velocity flow create 1.10 482748 temp mobile units box
#fix 1 all nve

# thermostating

compute mobile_thermo all temp/partial 0 0 1
fix 2 all langevin 1.10 1.10 1.0 699483
fix 1 all nve
fix_modify 2 temp mobile_thermo
thermo_modify temp mobile_thermo

# Couette flow

velocity lower set 0.0 0.0 0.0 units box
velocity upper set 1.0 0.0 0.0 units box
fix 3 boundary setforce 0.0 0.0 0.0

variable pxy equal pxy

fix 5 flow ave/spatial 1 500000 500000 x center 11.95 z center 7.22 y 1.5 1.02375 vx region fluid file vel_fb.profile units box
fix 6 all ave/correlate 1 50000 50000 v_pxy type auto file pres_fb.dat

# Run

timestep 0.005
thermo 5000

thermo_style custom step temp etotal c_mobile_thermo v_pxy

run 1000000

2

Dear Joseph,

To average over time, you'd better use fix ave/time.

But beware that the definition of pxy is only valid for a bulk system.
Therefore pxy is ill-defined in a confined system as the one you are
considering. In fact, you can simply measure the shear stress as the
shear force on the walls divided by the walls surface, something like
(I haven't check it):

fix 3upper upper setforce 0.0 0.0 0.0
fix 3lower lower setforce 0.0 0.0 0.0
variable shearupper equal -f_3upper[1]/lx/lz
variable shearlower equal f_3lower[1]/lx/lz
fix 6 all ave/time 10 10000 100000 v_shearupper v_shearlower file
name_of_the_file

Of course pressure and stresses always fluctuate a lot in MD
simulations, so you'll need to average over very long times.

Best,
Laurent

3

Thank you for your suggestions. I shall try it out.

The reason I was using Pxy as the paper I am trying to reproduce(by P A Thompson & S M Troian)
says that viscosity was calculated by viscosity = Pxy/Shear rate. Shear rate being calculated from the slope of the velocity profile.

Also on a different note I wanted to confirm if the damp parameter in 'fix langevin' is reciprocal of the friction constant(Gamma) given in the langevin equations mentioned in Schneider & Stoll. In the paper they say that 1/Gamma must be much much smaller than Tch-equilibrium chain length and much much greater than Tc-characteristic time of the dynamics. I am not sure how to calculate Tch or Tc. Is Tc the same as characteristic time of L-J potential. If yes then Tc comes as 1 for LJ units, hence the criteria given is not being satisfied.

BR/Joseph

4

Also on a different note I wanted to confirm if the damp parameter in
'fix langevin' is reciprocal of the friction constant(Gamma) given in
the langevin equations mentioned in Schneider & Stoll. In the paper they
say that 1/Gamma must be much much smaller than Tch-equilibrium chain
length and much much greater than Tc-characteristic time of the
dynamics. I am not sure how to calculate Tch or Tc. Is Tc the same as
characteristic time of L-J potential. If yes then Tc comes as 1 for LJ
units, hence the criteria given is not being satisfied.

I don't know without reading their paper. But the fix langevin doc
page clearly states what it is. And you have the source code
for src/fix_langevin.cpp, so you can check out the details.

Failing that, you can simply run with different Tdamp params and
monitor the Temp and see what works for your system. It's
a fudge factor to start with.

Steve

5

thanxx

I was referring to the paper of (Schneider) Schneider and Stoll, Phys Rev B, 17, 1302 (1978) that has been mentioned in the Lammps manual "Apply a Langevin thermostat as described in (Schneider)" under description for fix Langevin command.

BR/Joseph