Ncount bin sensitivity

Hello,

I am new to LAMMPS, but I have generated a Couette flow and the profile is almost linear (attached), however the profile is slightly non-linear near x=0.1 and 0.9. I noticed that the number of counts (Ncount) is approximately half at x=0.1 and 0.9. Is there any reason why there should be a discrepancy at the ends of the profile? These are the lines I am using to generate the profile:

run 100000

compute cc1 flow chunk/atom bin/1d y center 0.05 units reduced
fix 8 flow ave/chunk 2 50000 1100000 cc1 vx vz norm sample file vel.profile

run 1000000

Is there any way to make the number of counts the same in all bins? Or is there another reason for the difference - e.g. bin sensitivity? (I think I have enough timesteps)

Kind regards,

Andrew

vx_vs_y.png

Hello,

I am new to LAMMPS, but I have generated a Couette flow and the profile is
almost linear (attached), however the profile is slightly non-linear near
x=0.1 and 0.9. I noticed that the number of counts (Ncount) is approximately
half at x=0.1 and 0.9. Is there any reason why there should be a discrepancy
at the ends of the profile? These are the lines I am using to generate the

my guess would be, that this is likely due to the ordering and
excluded volume near the walls. fix ave/chunk averages over points,
i.e. the position of the center of your particles. just make some
snapshots of your system and check manually.

axel.

Also, is the total count for all the bins equal

to the total # of particles you expect to be counting?

From one snapshot in a dump file, you should

be able to manually count the particles in the

first/last bin and verify what LAMMPS is counting for the same.

I.e. use fix ave/time with no averaging.

Steve

Hi Steve,

Yes I just checked that the number of particles is 1002 as defined by the number of counts in each of the 20 bins across the profile and yes, this equals the total number of particles that were originally created. To prove it, this is the snapshot of the polydisperse system:

https://www.dropbox.com/s/sw6jvvccrtepq0l/atoms_2.png?dl=0

I have uploaded my profile in comparison with the literature:

https://www.dropbox.com/s/p1zonm686lxzzoq/vx_vs_y_2.png?dl=0

There seems to be some kind of crystallization near the walls for the LAMMPS simulation, although the central flow region seems to be fairly linear. The daCruz case also contains no friction - presumably this means that the particles have no tangential force on them and if there is also no tangential damping, they cannot rotate. Could a tightly packed wall encourage crystallisation if the particles cannot rotate free of this zone - like raking the particles through the flow? I noticed that the literature ignores the zones close to the wall and has possibly rescaled the profile to show only the linear zone, but I’m not sure if that is the reason. Do you know what could cause this problem? Presumably increased polydispersity would cause more slippage at the wall and less polydispersity would cause more crystallisation?

But why would re-ordering the wall particles influence the velocity profile? - but it’s true that the wall particles are added first in order, then the flow particles are randomly placed between them. Should the wall particles be added with random id numbers?

Also, why is there an excluded volume near the walls? - the velocity profile is being computed for all atoms in the simulation.

My input script looks like this:

https://www.dropbox.com/s/6zltghadqx1kq7w/input.txt?dl=0

Kind regards,

Andrew

I’d encourage you to match the de Cruz set-up exactly if you can. From your visualization with very nearly hexagonal packings in many places, I wouldn’t be surprised by your velocity profile. As Axel states, its probably due to ordering near your walls. To answer some of your other points:

There seems to be some kind of crystallization near the walls for the LAMMPS simulation, although the central flow region seems to be fairly linear. The daCruz case also contains no friction - presumably this means that the particles have no tangential force on them and if there is also no tangential damping, they cannot rotate. Could a tightly packed wall encourage crystallisation if the particles cannot rotate free of this zone - like raking the particles through the flow?

You do not require friction to have crystallization. However, friciton on the particle DOES affect many things in these systems from anisotropy in ordering, stress, instabilities, etc. I suggest matching these between comparisons as best you can.

I noticed that the literature ignores the zones close to the wall and has possibly rescaled the profile to show only the linear zone, but I’m not sure if that is the reason. Do you know what could cause this problem?

Its typically ignored for exactly the reason you are seeing it - they are interested in ‘statistically homogeneous simple shear’. Your walls are rigid and ordered, why would you expect confined particles to not have similar ordering near them. Most of these ideas are just geometry - how many ways can additional circles be arranged in a space if there are additional particles that cannot be moved - compared with a case where there are no immobile circles.

Presumably increased polydispersity would cause more slippage at the wall and less polydispersity would cause more crystallisation?

The idea that polydispersity decreases crystallization is generally accepted in the DEM community. Slippage is a bit more vague in my opinion.

But why would re-ordering the wall particles influence the velocity profile? - but it’s true that the wall particles are added first in order, then the flow particles are randomly placed between them. Should the wall particles be added with random id numbers?

The only thing affecting your velocity profile is that you essentially have a non-yielding solid layer near the walls that moves at roughly the same velocity - That is what your velocity profile tells me. From your input script, re-ordering your commands should not change a thing - this is VERY likely NOT an averaging/implementation error. You’re welcome to confirm this. Try seeing what your velocity profile looks like if you shear with a triclinic box instead, if you still require more convincing.

Also, why is there an excluded volume near the walls? - the velocity profile is being computed for all atoms in the simulation.

If you still have any questions on what is meant by ‘excluded volume,’ please do some searching - if your dealing with DEM (relatively hard particles) you will benefit from understanding what this means.

As a side point, the primary benefit of shearing with walls is that you can easily control the imposed normal force without manipulating the spaces in between relatively hard particles throughout the domain as is done with a barostat. You don’t need to use them if that’s not your concern. Gravity can also be a concern depending on the problem.

HTH,

Hi Eric,

Thanks a lot for your help. I am new to DEM, so sorry if I sound a bit confused, you definitely made things clearer. I read another version of the paper, and they say that they ignored the first 5 layers from the wall and chose H to limit the wall effects.

I also need to measure the pressure on the upper wall. Do you know how to measure the pressure on the upper wall? This is the snapshot, with the upper blue wall moving at velocity 1 in direction x.

https://www.dropbox.com/s/sw6jvvccrtepq0l/atoms_2.png?dl=0

So far I have this:

Pressure

compute 11 all temp
compute 12 all pressure 11

Pressure Profile Output (upper is the upper wall)

compute cc2 upper chunk/atom bin/1d x center 0.05 units reduced
fix 13 upper ave/chunk 50 8000 451500 cc2 fy norm sample file force.profile

However, the file force.profile contains zeros, I’m assuming because I set my force to zero using:

fix 8 boundary setforce 0.0 0.0 0.0

Is there a fix for pressure to output pressure.profile instead of force.profile? Or do I do something with the compute with id 12?

This is the rest of the script:

https://www.dropbox.com/s/xtla5nazgxsdhba/input2.txt?dl=0

Sorry for my naive questions.

Kind regards,

Andrew

You are correct in that your forces on your wall particles are being zeroed out by your fix setforce command. There are a few alternatives of course - outputting the forces between all particles with compute pair/local and selecting only those which interact with the walls - making use of Newton’s 3rd law. You could also make use of compute stress/atom and do some appropriate averaging to get the pressure in the bulk - probably more appropriate if you are trying to match a stress/strain rate plot.

There is also a fix store/force command which

can store forces on particles before they are

modified by other fixes. They you could

dump these “original” forces on the wall

particles with dump snapshots. See the

doc page for the command.

Steve