Velocity of Particles

Submitted by mhk on Tue, 11/12/2013 - 06:38

Hi everyone,
I have tried to calculate the velocity of particles based on their position, but the results show that the velocities which I calculated are exactly twice the dumped velocity. I have checked the time interval, I used the real time as follows:
"real time=timestep*No. of steps"
I was wondering if anyone has any idea why this happens?or what am I doing wrong?

Cheers,
Mahsa

rahulsoni | Fri, 11/15/2013 - 13:04

mhk

Its strange to see this post. Anyway, to crosscheck this, I simplified one of my complex in.script to get only 1-2 particles, turned off the gravity, given a constant velocity to particles and the results I observed are correct. Dump files recorded the correct velocity and it matched with the calculation that I made with the real time and also with the natural principles of physics.

I do not understand exactly which method you followed to calculate but there might be some issues with them:
1. time step should be small enough to make simulation stable. Better to 20-30 % of Rayleigh time.
2. Check the frequency (dump dmp all custom 100) of dump files being dumped. So, real time = no. of dump files x steps (for dump) x time step
3. If you directly looked into raw data then missing the other component of vector.
4. If you have used a geometry then scaling factor may have a role.

Please explain exactly how did you calculated so that I can understand the issue.

All the best..

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

mhk | Mon, 11/18/2013 - 00:45

Actually I have checked almost every thing but I could not find where my mistake is.
Regarding No.1, I use the command below:
fix ts all check/timestep/gran 1000 0.1 0.1
so, I don't get any warning message.

regarding No. 2, could you please explain more? why do we need to multiple by No. of dump files? below is the line that I have used for dumping the position of particles:
dump d10 all custom 1000 dump_position.liggghts id type type x y z

Regarding no.4 here is my related line of script:
fix cad1 all mesh/surface/stress file opentwoends.stl type 2 scale 1.0 curvature 1e-5

does the curvature effect the position components?
Could you please explain more about No. 3?I don't get your point.
Thanks again

mhk | Mon, 11/18/2013 - 03:47

to make it clear, here is dump of position for two timesteps(2600000 and 2601000) to calculate dx/dt(dt=timestep(0.000001)*1000):

ITEM: TIMESTEP
2600000
ITEM: NUMBER OF ATOMS
910
ITEM: BOX BOUNDS pp mm mm
-0.5 0.5
-1 1
-1 1
ITEM: ATOMS id type type x y z
353 1 1 -0.122143 -0.0317908 -0.42442
803 1 1 -0.0115618 -0.0150474 -0.42488

2601000
ITEM: NUMBER OF ATOMS
910
ITEM: BOX BOUNDS pp mm mm
-0.5 0.5
-1 1
-1 1
ITEM: ATOMS id type type x y z
353 1 1 -0.120459 -0.0301441 -0.425168
803 1 1 -0.00986969 -0.0127653 -0.42536

Here is dump of velocity for the same timesteps (2600000 and 2601000):

2600000
ITEM: NUMBER OF ATOMS
910
ITEM: BOX BOUNDS pp mm mm
-0.5 0.5
-1 1
-1 1
ITEM: ATOMS id type type vx vy vz
353 1 1 0.909155 0.728317 -0.673066
803 1 1 0.85725 1.11407 -0.374521
ITEM: TIMESTEP
2601000
ITEM: NUMBER OF ATOMS
910
ITEM: BOX BOUNDS pp mm mm
-0.5 0.5
-1 1
-1 1
ITEM: ATOMS id type type vx vy vz
353 1 1 0.768014 0.93035 -0.0406483
803 1 1 0.833402 1.16876 -0.086947

rahulsoni | Mon, 11/18/2013 - 09:44

Dear Mahsa

Regarding 1: Okay
Regarding 2: This is something similar I suggested as you have already done as dt=timestep(0.000001)*1000
Regarding 3: Will be cleared later
Regarding 4: Okay

As per the information you have provided the magnitude of velocity that I found based on dx/dt(dt=timestep(0.000001)*1000) is V=2.4712306428 (for atom 353).
However, the magnitude of V1 and V2 as per the dump files are 1.3453721815 and 1.2070828103, respectively. Please see the spreadsheet for calculation details

.

If you see in spreadsheet or even in general sense then you see that cosines of positions and velocities differs from each other. Meaning all have different directions (an essential part of the vector). Now, V calculated as 2.4712306428 is the average velocity by which the particle has moved from position 1 to position 2. However, the V1 in the dump file is velocity that particle possessed at position 1 and V2 is that particles possesses at position 2. Ironically, this two velocities does not have any relation with the V. Please note the dump files are written with the gap of 1000 steps. Several collisions might have happened during this period, particle could possess different values and directions of velocity. That's why I written V as the average velocity.

If V1 and V2 would have been in the same direction then youer methodology of calculation might work (might because we can't be sure that there are no collisions in between).

Vectors must be handled properly with their direction as well. Therefore in the point no. 2 of my first post I echoed about the vectors and their components.

So, the simulation is not doing anything wrong rather the approach you have adopted to compare is not correct!!!

Feel free to ask if any doubt.

All the best :)

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

mhk | Tue, 11/19/2013 - 02:02

Thanks again for your expanded explanation. Actually, I am not saying that the simulation is doing something wrong. Simply, I don't know where I am doing wrong. I have not checked it for total velocity, I did only in x direction. In the link below you can find my calculations.
https://docs.google.com/spreadsheet/ccc?key=0AqOEyhOKsKvodEZfWGxhSTFmRlo...

I have done similar calculations for all my particles. I also plot the velocity distribution. they are exactly twice. I mean the velocity based on dx/dt is twice mean instantaneous velocities. Basically, the displacement in x direction gives the average velocity so the mean of instantaneous velocities in x direction should be equal to average velocity. Don't you think so? As it's exactly twice, I guess there might be something wrong somewhere. What do you think?

rahulsoni | Tue, 11/19/2013 - 04:40

Mahsa, sorry for incorrectly getting your query and taking into the wrong direction.

Anyway the spreadsheet you have linked is not shared with me. My email is

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

rahulsoni | Tue, 11/19/2013 - 05:46

Yes, I can access it now. Thanks

I checked for all three directions separately and they are coming twice.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

mhk | Tue, 11/19/2013 - 10:09

that is exactly my problem! I am thinking of making a mistake in my input script, but I have no idea how or where. you said you tried it for your case and it was correct. would you mind to check it for me again?if yours works fine, may I ask for a copy of your input script (I mean the related lines)? actually, I have copied the related lines of my script above and you said they are OK, so I don't know what is wrong.

rahulsoni | Tue, 11/19/2013 - 13:27

Okay Mahsa, here it is. I willl explain first what I did. Actually, I have a hexagonal shape of drum then released two particles from some height. Gravity was off and particles was given the velocity 1 m/s. I simplified the simulation to get the clear picture.

dump1000.mixer
ITEM: TIMESTEP
1000
ITEM: NUMBER OF ATOMS
2
ITEM: BOX BOUNDS ff ff ff
-0.08 0.08
-0.17 0.17
-0.17 0.17
ITEM: ATOMS id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
1 1 1 0.0347623 -0.0564282 0.069001 0 0 0 0 0 -1 0 0 0 0 0 0 0.0025
2 2 2 0.0347623 0.0485718 0.069001 0 0 0 0 0 -1 0 0 0 0 0 0 0.0025

dump2000.mixer
ITEM: TIMESTEP
2000
ITEM: NUMBER OF ATOMS
2
ITEM: BOX BOUNDS ff ff ff
-0.08 0.08
-0.17 0.17
-0.17 0.17
ITEM: ATOMS id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
1 1 1 0.0347623 -0.0564282 0.068001 0 0 0 0 0 -1 0 0 0 0 0 0 0.0025
2 2 2 0.0347623 0.0485718 0.068001 0 0 0 0 0 -1 0 0 0 0 0 0 0.0025

Motion is only there in z-direction.
dz/dt=(0.068001-0.069001)/(1000x0.000001)=-1 m/s
(Vz_1+Vz_2)/2=(-1+(-1))/2=-1 m/s

in.script [Find .stl file at https://docs.google.com/file/d/0B56LTL-18kBaRmZMNnNPQk1rdjA/edit]
# Moving mesh example

atom_style granular
boundary f f f
newton off

communicate single vel yes
units si

#region
region reg block -0.08 0.08 -0.17 0.17 -0.17 0.17 units box
create_box 3 reg

neighbor 0.0025 bin
neigh_modify delay 0

#Material properties
fix m1 all property/global youngsModulus peratomtype 7.e10 7.e10 3.e9
fix m2 all property/global poissonsRatio peratomtype 0.22 0.22 0.35
fix m3 all property/global coefficientRestitution peratomtypepair 3 &
0.67 0.67 0.67 &
0.67 0.67 0.67 &
0.67 0.67 0.67
fix m4 all property/global coefficientFriction peratomtypepair 3 &
0.95 0.95 0.8 &
0.95 0.95 0.8 &
0.8 0.8 0.5
fix m5 all property/global characteristicVelocity scalar 0.5.

#pair style
pair_style gran/hooke/history #Hooke without cohesion
pair_coeff * *

timestep 0.000001

fix nveint all nve/sphere
fix grav all gravity 0.0 vector 0.0 0.0 -1.0

#box walls
fix boxwalls_x1 all wall/gran/hooke/history primitive type 3 xplane -0.075
fix boxwalls_x2 all wall/gran/hooke/history primitive type 3 xplane +0.075
fix boxwalls_y1 all wall/gran/hooke/history primitive type 3 yplane -0.165
fix boxwalls_y2 all wall/gran/hooke/history primitive type 3 yplane +0.165
fix boxwalls_z1 all wall/gran/hooke/history primitive type 3 zplane -0.165
fix boxwalls_z2 all wall/gran/hooke/history primitive type 3 zplane +0.165

#import mesh from cad:
group cad_group_1 region reg
fix cad1 cad_group_1 mesh/surface file Drum_Hexagonal.stl type 3 move 0. 0. 0. &
scale 0.001 rotate axis 0. 1. 0. angle 0. rotate axis 0. 0. 1. angle 0.

fix splitter_wall all wall/gran/hooke/history primitive type 1 yplane 0.0

#particle insertion
region bc_1 block -0.068 0.068 -0.10 -0.005 0.06 0.07 units box
group nve_group_1 region reg
fix ins_1 nve_group_1 pour/legacy 1 1 1 vol 0.5 10 diam 0.005 0.005 dens 2700 2700 vel 0. 0. 0. 0. -1.0. region bc_1

#particle insertion
region bc_2 block -0.068 0.068 0.005 0.10 0.06 0.07 units box
group nve_group_2 region reg
fix ins_2 nve_group_2 pour/legacy 1 2 1 vol 0.5 10 diam 0.005 0.005 dens 2700 2700 vel 0. 0. 0. 0. -1.0. region bc_2

#thermo settings
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol
thermo 100
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

#make a dump of particles and the stl file
dump dmp all custom 100 post/dump*.mixer id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstl cad_group_1 mesh/stl 100 post/dump*.stl cad1

#restart command
restart 100000 tmp.restart

#use the imported mesh (drum) as granular wall
fix drum_mixer all wall/gran/hooke/history mesh n_meshes 1 meshes cad1
run 1 upto

#run for particle insertion
run 3000000 upto
#stopping particle insertion
unfix ins_1
run 1
unfix ins_2
run 1

#removing drum_splitter
unfix splitter_wall
run 1

#now lift the particles up
fix rotate_drum all move/mesh mesh cad1 rotate origin 0. 0. 0. axis 1. 0. 0. period 14.0.
run 100000000

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

rahulsoni | Tue, 11/19/2013 - 13:32

Mahsa, it will be really helpful if you can help me briefing about your simulation. What exactly is being done. I there would be catch in between. Have you checked the similar with other set of steps (other than 2600000 and 2601000)?

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

mhk | Wed, 11/20/2013 - 05:59

>> Have you checked the similar with other set of steps (other than 2600000 and 2601000)?
Yes I have checked for different steps and different frequencies of dumping.
I also checked your input script and compared with mine. there are only minor differences e.g. I have used "hooke" 's model for forces, I have used "neighbor 0.01 bin", "scale 1.0" for imported mesh, I have not used box walls I only used the imported mesh as wall, I used periodic boundary conditions in x direction.

rahulsoni | Wed, 11/20/2013 - 08:36

Mahsa, we both know that hooke's model, neighbor, scale and boxwalls are irrelevant with this problem. Can you please tell me Mahsa what actually is your simulation. I am suspecting that the catch may be among the energy distribution while simulating. 1000 steps is a big difference between two snapshots, therefore it is good to know whats happening inside.

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

rahulsoni | Fri, 11/22/2013 - 07:06

Sorry Mahsa, for delayed reply. But I am trying all possible combinations but could not yet debug the issue. Tried to think whether it is happening because of rotation, or the direction of gravity and similar ones. But, none of them of worked. Surprisingly, the simulation appears to be unstable for me. I simplified it, made the gravity only in z-direction, started dumping the file from very beginning.

Case I: Gravity was on according to my calculation the particle its drop should have touched the bottom of drum after 400000 steps approx. But the particle did it after 200000 only.

Case II: Gravity was turned off after 10000 steps because I wanted to allot a constant velocity of 0.1 m/s. Everything goes fine but when I looked into the position dump files it shown the constant velocity as 0.1962 m/s. Also, the velocity calculated by dx/dt was 0.392 m/s.

I followed the following script:
(mine version of stl file can be found at https://drive.google.com/file/d/0B56LTL-18kBaREtvWFdpRGFXZDA/edit?usp=sh... )

# Ball mill simulation

atom_style granular
boundary p m m
newton off

communicate single vel yes
units si

#region
region reg block -0.5 0.5 -1 1 -1 1 units box

create_box 2 reg

#neighbor 0.01 bin
neigh_modify delay 0

#Material properties required for new pair styles

fix m1 all property/global youngsModulus peratomtype 1.e7 2.e11
fix m2 all property/global poissonsRatio peratomtype 0.27 0.3
fix m3 all property/global coefficientRestitution peratomtypepair 2 0.7 0.7 0.7 0.7
fix m4 all property/global coefficientFriction peratomtypepair 2 0.3 0.3 0.3 0.3
fix m5 all property/global coefficientRollingFriction peratomtypepair 2 0.1 0.1 0.1 0.1

#New pair style
pair_style gran/hertz/history rolling_friction cdt #Hertzian with rolling friction and without cohesion
pair_coeff * *

timestep 0.000001

fix 1 all nve/sphere
#fix 2 all gravity 9.81 vector 0.0 0.0 -1.0

#import mesh from cad:
fix cad1 all mesh/surface/stress file opentwoends.stl type 2 scale 0.00667 curvature 1e-5

#use the imported mesh as granular wall with rolling Friction
fix bucket_wall all wall/gran/hertz/history mesh n_meshes 1 meshes cad1 rolling_friction cdt

#distributions for insertion/ Partikelverteilung
fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 7850 radius constant 0.025
fix pts2 all particletemplate/sphere 1 atom_type 1 density constant 7850 radius constant 0.005
fix pdd1 all particledistribution/discrete 1000 2 pts1 1 pts2 0.0

#region and insertion
region bc block -0.3 -0.1 -0.1 0.1 -0.2 0.3 units box
group nve_group region reg
fix ins nve_group insert/pack seed 100 distributiontemplate pdd1 insert_every once overlapcheck yes particles_in_region 1 region bc ntry_mc 10000

#apply nve integration to all particles that are inserted as single particles
fix integr nve_group nve/sphere

#fix ts all check/timestep/gran 1000 0.1 0.1

#thermo settings
compute 1 all erotate/sphere
thermo_style custom step atoms ke c_1 vol dt
thermo 1000
thermo_modify lost ignore norm no
compute_modify thermo_temp dynamic yes

#contact details
compute fc all wall/gran/local
dump forcedump all local 10000 dump_force.liggghts c_fc[1] c_fc[2] c_fc[3] c_fc[4] c_fc[5] c_fc[6] c_fc[7] c_fc[8] c_fc[9] c_fc[10] c_fc[11] c_fc[12] c_fc[13] c_fc[14] c_fc[15] c_fc[16] c_fc[17]

dump d0 all custom 1000 post/dump_*.liggghts id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
dump dumpstl all stl 1000 post/dump*.stl
dump d10 all custom 1000 dump_position.liggghts id type type x y z
dump d20 all custom 1000 dump_velocity.liggghts id type type vx vy vz

#insert the first particles so that dump is not empty
run 1
#unfix ts
unfix ins

fix 2 all gravity 9.81 vector 0.0 0.0 -1.0
run 10000
unfix 2 #awarding a constant velocity ~0.1m/s

#run with particle insertion
run 480000 upto
#unfix ins

#run to let particle settle
run 160000

#unfix 2
#fix 3 all gravity 9.81 vector 0.1745 0 -0.985

#moving mesh, rotate
fix movecad2 all move/mesh mesh cad1 rotate origin 0. 0. 0. axis 1. 0. 0. period 1.5
run 2600000 upto

run 3000000 upto

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

rahulsoni | Fri, 11/22/2013 - 07:54

Mahsa, after so much brainstorming for last 2-3 days. I have reached to very very silly and simple solution. Feeling ashamed on myself.

Go back to your in.script and find that you have given following two lines:
fix 1 all nve/sphere
fix integr nve_group nve/sphere

That means it is integrating or calculating each and everything two times. Please make a note that a precedence of calculation is acceleration followed by velocity followed by position. Now, whats happening here is the velocities that are dumped in the velocity file is correct. Now, the dx in dx/dt has got doubled.

In general:
dx = ( v1 + v2 ) / 2 x dt
So, in general dx is dx but in your case it is done for the extra second time that means dx was calculated then again one more dx was added resulting in 2dx. And therefore, the velocities are coming doubled.

PS: nve integration command can be placed anywhere (no need to put after particle insertion). Go and modify the script and get the correct results.

Anyway, thanks Mahsa, this problem has taught me a big thing. Any silly mistake in simulation can be disaster and sometimes you never that it has happened..

:):):)

--
Regards

Rahul Kumar Soni
Scientist, CSIR - IMMT, India

mhk | Thu, 11/28/2013 - 02:23

Thank you for following up. I really appreciate it. I didn't know if we write it twice it will update twice! :)

mhk | Thu, 11/28/2013 - 03:37

I was wondering if you also know how I can get the dissipative force for collision of particles from LIGGGHTS? I need to know how much energy is lost in a period of time.
Thanks again for your help.