Velocity issue

Submitted by robin on Tue, 02/10/2015 - 13:48

Hello everybody,

I'm working on a very simple test-case (see input file below). Two particles, one on the top of the other with a normal spring between them (pair_style gran model/hooke/stifness) are compacted in one direction (fix comp all deform 1 z trate -1).

In the dump file, I get the expected positions and forces (components z and fz), but the velocity in the z direction (vz) is 0 for both particles at all time steps.

What am I doing wrong?

Thanks for your help.

Robin

##########################
#Input file:
atom_style granular
atom_modify map array
boundary f f f
newton off
echo both

communicate single vel yes
units si
region reg block -0.5 0.5 -0.5 0.5 -0.5 0.5 units box
create_box 1 reg
neighbor 0.5 bin
neigh_modify delay 0

#Insertion of the particles
create_atoms 1 single 0 0 0.25
create_atoms 1 single 0 0 -0.25
set region reg diameter 0.5

#Material properties required for new pair styles
fix m4 all property/global coefficientFriction peratomtypepair 1 0.3
fix m6 all property/global kn peratomtypepair 1 10
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 * *
timestep 0.00001

#output settings, include total thermal energy
compute 1 all ke/atom
thermo_style custom step atoms ke vol
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

fix comp all deform 1 z trate -1
dump dmp all custom 2000 post/dump*.packing id type x y z vx vy vz fx fy fz radius
run 100000

Forums:

aaigner's picture

aaigner | Tue, 02/17/2015 - 17:55

There is no integration at all?! If you don't integrate the particles, the velocities are not updated, so add a
fix integr all nve/sphere

Cheers
Andi

robin | Tue, 02/17/2015 - 19:38

OK, thanks to you, I can now compute velocities.

However, I think I still don't understand what's going on. Indeed if I set :
-> kn peratomtypepair 1 0
Velocities go back to zero again.
Second problem, even if I set for example :
-> gamman_abs peratomtypepair 1 10
Forces also drop to zero if kn is set to zero.

What am I misunderstanding? Is a non-zero kn value necessary for the computations?

Thanks!
Robin

JoshuaP | Wed, 02/18/2015 - 10:33

if kn is 0 the contact force is always calculated also as 0, so no force is transmitted at the contact.
Contact force is always stiffness times overlap. I'm interested in what case do you wanna simulate a 0 stiffness?

regards

robin | Wed, 02/18/2015 - 15:44

1) Contact force and velocity
The contact force being zero doesn't imply that the velocity should be zero. As far as I understand, without interaction forces, the velocity should be directly set by the "fix deform".

2) Kn and normal force
In the hooke/stifness model you have two contributions in the normal force : the elastic and the damping part. By setting kn=0 (elastic part) gamman_abs=10 (damping part) you should get some interaction forces.

3) Why bother?
In order to implement my own contact laws, I need to understand well enough how Liggghts works !

Cheers

JoshuaP | Wed, 02/18/2015 - 16:06

I think the fix deform changes the volume of the simulation box and doesnt deform the particles. If you want to see some interaction forces in your script you should set the gravity by (fix gravity).
Actually you are right there should be some interaction forces by the damping part.
regards

robin | Fri, 04/03/2015 - 14:42

Part of the answer is in the documentation of the fix deform command :
" 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. "

However it looks like there is no way to take into account this deformation to update the velocites. Indeed, if the keyword "remap v" is used, the velocities is not updated as I would want it to, but only if a periodic boundary is crossed:
" If remap is set to v, then any atom in the fix group that crosses a periodic boundary will have a delta added to its velocity "

In both cases, the velocites remains zero and thus the forces if only damping forces are set.

Cheers