K. Lykov Blog

LAMMPS: How to Compute Fluid Viscosity

The theory is described in “Poiseuille flow to measure the viscosity of particle model fluids” by J. A. Backer et al. Below I describe how to use this approach in LAMMPS.

(1) Run LAMMPS with the following script

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
boundary p p p

units     lj
atom_style    atomic

lattice custom 3.0 a1 1.0 0.0 0.0 a2 0.0 1.0 0.0 a3 0.0 0.0 1.0 &
      basis 0.5 0.0 0.0 basis 0.0 0.5 0.0 basis 0.0 0.0 0.5

region box block -7.0 7.0  -7.0 7.0  -14.0 14.0

region left  block -7.0 7.0  -7.0 7.0 -14.0 0.0
region right block -7.0 7.0  -7.0 7.0 0.0 14.0

# Uncomment it if you don't use restart file
create_box  1 box
create_atoms    1 box
mass        1 1.0

neighbor    0.3 bin
neigh_modify    delay 0 every 4 check no

#******************DPD******************
#to store velocities by ghost atoms
#communicate single vel yes - for old versions of lammps
comm_style brick
comm_modify vel yes

# T cutoff seed
pair_style  dpd 0.1 1.0 34387
# atom_type atom_type a gamma=sigma^2/2 cutoff(optional)
# where a is Fc coefficent.
pair_coeff  1 1 25.0 45.0 1.0

thermo          500
timestep 0.01

fix 1 all nve
fix 2 all addforce -0.055 0.0 0.0 region left
fix 3 all addforce 0.055 0.0 0.0 region right
fix 4 all ave/spatial 50 1000 50000 z center 0.5 vx file vel-visc.txt

run 100000

(2) Open vel-visc and copy in a separate document data for one time step.

(3) Open gnuplot, type:

1
gnuplot> plot "visc_vel.txt" using 2:4

The result should be something like:

(4) From analytical solution for the problem, it is known that . Where , p - is numeric density(3.0 in our case, determined by custom lattice), g is driving force (0.055), n - dynamic viscosity. In order to find alpha we will use gnuplot’s fit command. As you might see on the Figure above, there are 2 parabolas. I pick the left one, so the analytical solution look like . Then type

1
2
3
gnuplot> f(x)=a*(x*14 + x*x)
gnuplot> fit f(x) 'visc_vel.txt' using 2:4 via a
gnuplot> plot "visc_vel.txt" using 2:4, f(x)

The result should be , thus viscosity n=2.68 in DPD units. The plot with velocities from simulation and with the fitting plot should look like that: