Hi everybody,
I'm currently testing the "pair_style gran model hooke/stiffness", applying a true deformation rate on pairs of particles. I'm expecting, whatever the configuration, to get the theoretical result :
force = kn x indentation
In the elas.txt script (Source) I'm putting four pairs of particle side to side, in very simillar positions (initially I wanted to test the round-off behavior).
I then plot the force versus the indentation length (negative if overlapp) in various configuration (see files Case1 and Case2.pdf).
I'm surprised to find that not all points are on the theoretical line.
In the case 1 (timestep=5e-2 , kn=4e11) for certain time steps I get no interaction force at all, even though the particles are overlapping.
In the case 2 (timestep=5e-2 , kn=4e5) not all points are on the theoretical line, but a second line appears, with same slope, but with an offset.
More generally, only looking at the Liggghts output file, with no post-processing, I can find particle overlapp without contact force. It does not look like the overlapp is "too small" to be seen, I can get a correct force for the same order of magnitude of indentation (eg in Case 1).
Playing around with the time step (see for example the commented lines with smaller values) changes the behavior. It looks like the case 2, but with a smaller offset.
Here we go for the questions :
- How can I get an overlap without interaction force?
- Why do I get this offset?
- How a small time step can mess up the force computation?
Thanks for any hint, cheers,
Robin
| Attachment | Size |
|---|---|
 Source | 1.49 KB |
 timestep=5e-2 , kn=4e11 | 3.35 KB |
| timestep=5e-2 , kn=4e5 | 3.36 KB |
ckloss | Fri, 04/17/2015 - 09:40
Hi Robin,
Hi Robin,
after having a very quick look, I guess there is a problem with the setup:
>>applying a true deformation rate on pairs of particles
Fix deform will just deform the simulation box. However, it will not impose a wall boundary condition on the particles if you have "boundary f f f" as in your case
I would suggest you use a moving mesh to do the same thing. Moving mesh is strain-controlled, you can also use a fix mesh/surface/servo if you want a strain-controlled simulation
Christoph
robin | Tue, 04/21/2015 - 14:35
Imposed strain
Hi Christoph,
Fix deform changes the box size and as as wrote in this other thread (http://www.cfdem.com/forums/velocity-issue) :
" Each time the box size or shape is changed, the remap keyword determines whether atom positions are remapped to the new box. If remap is set to x (the default), atoms in the fix group are remapped; otherwise they are not. Note that their velocities are not changed, just their positions are altered. "
In the case I give, the boundaries are far from the particles and are not supposed to interact with them. Using a moving mesh would not be a handy solution, as for bigger stack of particles I will need a remapping of the position of all particle to accelerate convergence.
Anyhow this part, moving the particle, seems to work well. The problem is that for a given indentation, the result should not depend on the way you obtain this indentation, the force computed by Liggghts is not theoretical result.
I'll try to investigate further the cases where the force is not as expected and post it with a little Octave code which exhibits the expected behavior.
Cheers, thanks fo any hints!
Robin
robin | Wed, 04/22/2015 - 09:59
Simple case
Hi everybody,
I pasted at the end of the post a 2 particle case (very similar to the first I posted) where the interaction laws doesn't seem to be respected.
In the result file for the seventh time step I get :
ITEM: ATOMS id z fz radius
1 4.999999825003373e-04 1.199973167706481e+01 5.000000000000000e-04
2 -4.999999825003342e-04 -1.199973167706481e+01 5.000000000000000e-04
So the interaction force given by Liggghts is F_liggghts = 1.199973167706481e+01 N
If I compute "by hand" what the intraction force F_theory should be :
F_theory = Kn * indentation
F_theory = Kn * (distance - 2*radius)
F_theory = 4e11*((4.999999825003373e-04+4.999999825003342e-04)-(2*5e-04))
F_theory = 1.39997314088647e+01 N
And if I plot the result I get the same kind of behavior as in my initial post... I'll be glad if any of you can spot a mistake!
Cheers
Robin
______________________________________________________
atom_style granular
atom_modify map array
boundary f f f
newton off
echo both
communicate single vel yes
units si
region reg block -4.0e-3 4.0e-3 -4.0e-3 4.0e-3 -4.0e-3 4.0e-3 units box
create_box 1 reg
neighbor 1e-3 bin
neigh_modify delay 0
timestep 5e-5
#Insertion of the particles
create_atoms 1 single 0.0e-3 0.0e-3 0.5e-3
create_atoms 1 single 0.0e-3 0.0e-3 -0.5e-3
set region reg diameter 1e-3
set region reg density 1e18
#Material properties required for new pair styles
fix m4 all property/global coefficientFriction peratomtypepair 1 0
fix m6 all property/global kn peratomtypepair 1 4e+11
fix m7 all property/global kt peratomtypepair 1 0
fix m8 all property/global gamman_abs peratomtypepair 1 0
fix m9 all property/global gammat_abs peratomtypepair 1 0
pair_style gran model hooke/stiffness absolute_damping on limitForce off tangential_damping on
pair_coeff * *
fix integr all nve/sphere
fix comp all deform 1 z trate -1e-4
dump dmp all custom 1 postUnits/coord* id z fz radius
dump_modify dmp format "%d %.15e %.15e %.15e"
robin | Fri, 04/24/2015 - 13:07
Time step issue
Hello,
What I understand so far is that in the file given by the dump command for time step n :
- the position are given for the step n
- the forces for the step n+1
Does anyone know if it's actually the case? Can you confirm it's only a "printing" issue?
Thanks for any feedback,
Robin