Hi everybody,
Could you tell me how to dump the contact forces between granular particles.
Requesting your advice.
regards,
Samantaray
There is no current instrumentation in the code to do this.
The closest thing is compute pair/local which generates
info on pairwise forces and can produce output for the
dump local command. However, it only works for
total pairwise force, and the granular pair styles
would have to be augmented to have a single() routine
that compute pair/local would call. And some option
would need to be added to distinguish between normal
and tangential forces.
You could check if the LIGGGHTS add-on package
has instrumented the code for this. They have
added many features to the granular capability
in LAMMPS. It is linked to from the LAMMPS WWW
page under the Related Tools heading.
Steve
I just added single() functions to the granular
pair styles which compute the normal force
between pairs of particles, as well as additional
values for the tangential force components.
These can all be accessed by additions
to the compute pair/local command, so
that they can be dumped via the dump local command.
See the granular pair style and compute pair/local
doc pages for details.
This is in a 1Dec11 patch.
Please try it out and tell me if the values are what
you expect. Doing the right thing with the shear
history values is a bit tricky, so it's possible it's
still not quite right.
Steve
Thank you Steve for the function to dump tangential forces which was difficult to find out.
I also wish to get the contact normal forces. I think it will be easier.
But I dont know till now how to dump these contact normal forces for the granular style.
Best regards,
Samantaray
See the doc page for pair_gran.html. The "force"
output of the single() func, e.g. as accessed by
compute pair/local, is just the normal forces.
So that you can also get components fx,fy,fz
in the usual manner (multiply by delx,dely,delz).
You couldn't do this if the "force" were not along
the axis between the 2 particles.
The extra components p1,p2,p3 are the tangential
forces. So if you add things together (e.g. via
a per-atom variable) you can also get the total
force if you like.
Steve
Hi everybody,
I have recently written the following code for the strain controlled shear of granular material.
I am facing the following problems.
- I could not visualize the result in pizza by using the triclinic box.
- I dont think the result is correct from the position of the particles in the dump file.
Could you please help me.
Best regards,
Samantaray.
###########Code starts here #####################
units cgs
atom_style sphere
atom_modify map array
boundary f p p
newton off
lattice sc 0.011
region reg prism -0.04 0.2 -0.04 0.2 -0.06 0.06 0 0 0.02 units box
create_box 2 reg
create_atoms 2 region reg
set group all diameter 0.01 density 2.0
communicate single vel yes
neighbor 0.001 bin
neigh_modify delay 0
pair_style gran/hertz/history 20546.0 NULL 15660.46 NULL 0.5 1 #Hertzian with cohesion
pair_coeff * *
fix xwalls all wall/gran 20546.0 NULL 15660.46 NULL 0.5 1 xplane -0.04 0.2
timestep 0.000002
fix gravi all gravity 981 vector 1.0 0.0 0.0
variable dump_every equal 100000
fix integ all nve/sphere
compute str all stress/atom ke pair
compute p all reduce sum c_str[1] c_str[2] c_str[3]
variable presse equal -(c_p[1]+c_p[2]+c_p[3])/(3*vol)
compute forces all pair/local dist fx fy fz p1 p2 p3 p4
thermo_style custom step atoms temp etotal ke vol press v_presse
thermo ${dump_every}
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes
run 1
run 50000 upto
fix 2 all deform 1 yz erate 0.01 remap v units box
dump dmp2 all custom ${dump_every} post/dump.fab_mark id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
run 50000
######### end of code ############################
I visualize triclinic geometries all the time with Pizza.py
so it can do it. I would first verify that the coords
in the dump file are what you expect and that the
bounding box it prints is also correct. You should
be able to do this for a simple small system.
Steve
Hi Steve,
Thank you for the reply and best wishes for the new year.
pizza is now working fine for the triclinic box.
I have following doubts regarding the code for shearing.
The following command created the triclinic box for shearing in horizontal y direction.
boundary f p p
region reg prism -0.04 0.2 -0.04 0.2 -0.06 0.06 0 0 0.0 units box
fix xwalls all wall/gran 20546.0 NULL 15660.46 NULL 0.5 1 xplane -0.04 0.2
I have used following deform command for shearing
fix 2 all deform 1 yz erate 0.4 remap v units box
I dont see the shearing results in paraview.
Do I need to use move command for shearing the top plate.
fix F_movetop all move linear 0. 0.4 0. units box
How will I assign the topwall fix in the above move command.
Do I need to group the particles near the top wall for shearing?
How will I do that as I am using boundary f p p. So there are no walls in horizontal y and vertical z directions due to periodic boundary conditions.
Requests your help in running the shear code that I have been trying for some months.
Best regards,
Samantaray
Hi Steve,
Thank you for the reply and best wishes for the new year.
pizza is now working fine for the triclinic box.
I have following doubts regarding the code for shearing.1.
The following command created the triclinic box for shearing in horizontal y
direction.
boundary f p p
region reg prism -0.04 0.2 -0.04 0.2 -0.06 0.06 0 0 0.0 units box
fix xwalls all wall/gran 20546.0 NULL 15660.46 NULL 0.5 1 xplane
-0.04 0.2
none of these commands will create a box, let alone a
non-orthogonal one. where is the question?
2.
I have used following deform command for shearing
fix 2 all deform 1 yz erate 0.4 remap v units box
what is the question here?
3.
I dont see the shearing results in paraview.
Do I need to use move command for shearing the top plate.
there are two ways to induce shearing. 1: you tilt
the entire system (and the atoms with it), or 2: you
move atoms explicitly and don't change the box.
fix F_movetop all move linear 0. 0.4 0. units box
How will I assign the topwall fix in the above move command.
what has the wall to do with this?
4.
Do I need to group the particles near the top wall for shearing?
why do you think you might need this?
also, a group is just like a label. it doesn't
do anything to the system.
How will I do that as I am using boundary f p p. So there are no walls in
what do boundaries have to do with defining groups?
horizontal y and vertical z directions due to periodic boundary conditions.
so what? you still have a finite number of atoms.
Requests your help in running the shear code that I have been trying for
some months.
from reading your questions, it seems to me that your problem
is that you are trying to do something where you are lacking
lot of *fundamental* knowledge in doing MD simulations and
also seem to be following the (far too popular) approach of
just randomly guessing stuff until lit looks like you think it
should. this is leading to certain doom. there are just too many
subtleties that you cannot know unless you spend the time
and effort to learn them properly. the months that you mention
that you spent on trying to get this calculation going would
have been better spent on learning basics and experimenting
with simpler systems.
it is a good thing that you have (finally) started asking questions,
but a mailing list is no replacement for a classroom or a textbook
or a tutor/adviser. my suggestion is to spend a few weeks on
learning MD basics, e.g. going over the (simple) examples provided
with LAMMPS and trying to recreate simulations and their results(!)
described in an MD textbook with LAMMPS before you get back
to setting up a shear study. you'll see that suddenly things will
be much easier. also, i suggest you spend some time researching
through the mailing list archive. you are not the first person having
problems to set this kind of calculation up, and there may be lots
of valuable information in the previous discussions on the subject.
good luck,
axel.
Thanks Alex for the reply.
there are two ways to induce shearing. 1: you tilt
the entire system (and the atoms with it), or 2: you
move atoms explicitly and don’t change the box.
you tilt the entire system (and the atoms with it)
for this I used following commands
i. to make triclinic box
boundary f p p
region reg prism -0.04 0.2 -0.04 0.2 -0.06 0.06 0 0 0.0 units box
create_box 2 reg
create_atoms 2 region reg
ii. to tilt the entire system I used the following command
fix 2 all deform 1 yz erate 0.4 remap v units box
The full code is in the previous thread.
I need to know the missing steps as i dont see any shearing in paraview?
Best regards,
Samantaray
I don't use ParaView, and I imagine your problem
has to do with the viz, not with LAMMPS. If you
look at the box bounds printed in the dump file
you should be able to tell if you are indeed shearing
the box.
Steve