To anyone who has experience using the granular materials package,
I am trying to duplicate a granular materials simulation that has previously
been done using some custom software. This software has become outdated and not
easy to maintain, so the hope was to use LAMMPS for future work.
A basic simulation consists of 4000 particles randomly distributed inside of a
cylinder. Once these particles are settled at the bottom, the floor is raised
ever so slightly. The resulting stress on the particles is analyzed.
The issue I am running into is during the settling of particles. I can get them
to do one of two things. 1) The particles are pushed outside of the cylinder's
wall. I am guessing this is because the particle-particle force is too large.
2) The particles fall down and gather at the bottom, only to have half of them
pushed back upward. What really does not make sense is that in this situation,
the energy in the system increases. The energy gradually decreases as the
particles settle, but when this explosion happens, you can see the output "ke"
value increase.
I created my script from the in.pour example. That one runs fine. Perhaps the
problem is going from the default "lj" units to "si".
If it helps, I can post my script. I can also post videos of what is happening.
For reference, I am using LAMMPS version "11Jan10". Particle density 1220
kg/m^3. Time step = 0.0001 s. I use gran/hertz/history. Particle diameter =
0.1.
The first thing I would do is try to make the examples/in.pour
work exactly as it does now, but in SI units. That is possible,
but not trivial. You have to get every parameter correct in
its new units. If you can do that, then you should be able
to start adding features one at a time, like moving the wall, etc
and make sure it does what you expect. Even if I were doing
the simulation in LJ units, that's what I would do. Granular
particles are nearly rigid, so if you mess things up, you generate
tremendous forces.
A follow-up question - I tried searching the mailing list but could not find an
answer. How does one go about converting from LJ to SI? The documentation
talks about making use of a sigma and epsilon which are properties of the
materials being simulated. Unfortunately I am a CS major and not ME, MSE, etc.
and have no clue on where to look for this stuff. Also, are the values found in
update.cpp any use for this conversion?
Or is the typical result of this conversion a simulation that blows up, much
like I had before. Again, any tips are greatly appreciated.
I have some notes where I've done it myself for a couple
simple systems (real <-> LJ <-> metal, SI would be similar),
which I can send you. I should probably make them part
of the LAMMPS doc pages.
It really just involves carefully thinking about each input/output
quantitity and doing the units conversion in a consistent way.
Unfortunately, one mistake and you're hosed.
It will be next week before I have time to post it.
Does anyone know if LAMMPS can output per particle contact forces. I see that it
can output a force vector per particle. The end goal is to use these per particle
contact forces to calculate the stress on each particle. I see that LAMMPS
outputs a stress tensor, but those values may not be what I am looking for. If I
can just output the per particle contact forces, I have other code that does the
calculations that I need. So if it is possible to output these forces, that would
be the ideal situation for me right now.
To clarify "contact forces", I mean the forces between particles that are in
contact with each other. So if particles A, B, and C were all in contact with
each other, I would like LAMMPS to output something like:
A B force(AB)
A C force(AC)
B C force(BC)
From that information I can use reliable code to do everything else I need. The
results of "compute stress/atom" may not be what is desired. Unfortunately, I
am a CS person and I am not sure that the equation provided for "compute
stress/atom" is what I need. It very well could be though. But the values that
it outputs do not seem to be what I am looking for... perhaps its due to being
in LJ units and not SI like the other data I have.
To clarify "contact forces", I mean the forces between particles that are in
contact with each other. So if particles A, B, and C were all in contact with
each other, I would like LAMMPS to output something like:
with soft interaction potentials, the concept of "contact" is ill defined.
when are two particles far enough apart that they are not in contact?
A B force(AB)
A C force(AC)
B C force(BC)
this seems horrible inefficient to me. do you have an idea how many
lines of text this would be generating? if you want that kind of output,
you have to modify the source.
>From that information I can use reliable code to do everything else I need. The
results of "compute stress/atom" may not be what is desired. Unfortunately, I
am a CS person and I am not sure that the equation provided for "compute
sorry, but whether you are a CS person or not, doesn't change the physics.
you either know what you are doing or not.
stress/atom" is what I need. It very well could be though. But the values that
it outputs do not seem to be what I am looking for... perhaps its due to being
in LJ units and not SI like the other data I have.
if you are not even sure about units, then you won't get any improvements
from outputting the forces, since they will be in lj units as well
similar to the
stress tensor. note that LAMMPS *does* support running in si units.
Compute pair//local is close to what you need. Its calculations can be output
with dump local. However, contact forces are usually for granular
systems, which
means a pair of particles can have normal and tangential forces. The granular
pair styles do not ouput this info to compute pair/local (just simpler
pair styles do),
because it involves contact history and various complicated factors. So you
cannot use this command with the granular pair styles. It would have to
be extended to do so.