Hi All,
I am a new user.
I implemented the following simulation [in.packing]:
Box with periodic boundaries filled with 9261 particles with size and velocities read from data file I created [data.packing].
I am using a linear spring dashpot contact model. Quantities are scaled in such a way that are dimensionless: M=10-9 kg, L=10-3m, T=10-6s.
I acquire the dump files in Matlab and, for each particle, I compute the overlaps with other particles and the elastic forces (output from LIGGGHTS) acting on it.
The problem is: many contacts are missing (see pdf).
If I take the output from the first time step, compute the contacts in Matlab and check the forces in the output file of LIGGGHTS, the two results are in total agreement. If I go further with the simulation for 400 time steps, 400000 time steps, I have the problem I explained. While the time goes on more and more contacts are missing.
I was not confident about the setup of the neighbor search and contact model, so I played a lot with parameters like the skin value etc.
Finally, I use the same setup to run another simulation, in this case I create almost 4000 particles in LIGGGHTS, with exactly the same parameters (kn, timestep...), my comparison using Matlab gives positive check.
My questions are:
Why this contacts are missing?
Can I increment the accuracy of the output (more digits per variable on the output file)? I tried
dump_modify 1 format
but something was wrong. Any suggestions?
Thank you!
Giuseppe
| Attachment | Size |
|---|---|
 in.not_read.txt | 1.86 KB |
| data.packing.txt | 1.42 MB |
| in.packing.txt | 1.24 KB |
 matlaboutput.pdf | 51.9 KB |
giuraso | Fri, 01/16/2015 - 16:14
particle radius are between 0
particle radius are between 0.5 and 1.5, given by a uniform distribution.
JoshuaP | Fri, 01/16/2015 - 16:30
hi
you use Matlab for postprocess, so you get LIGGGHTS running on Windows?
I would prefer the commands from LIGGGHTS to calculate the number of contacts for each atom:
'compute ID group-ID contact/atom'
for ID u can use 'all'
and then you can put it to the dump file as 'c_group-ID'
btw doesnt it take long to let matlab calculate all distances between the atoms? I mean for 4000 particles you have ~16.000.000 possible contacts. Because I'm also trying to find a way to calculate the overlapped volume.
Regards
Joshua
giuraso | Sat, 01/17/2015 - 12:13
I am using Matlab (on Windows
I am using Matlab (on Windows), and LIGGGHTS is running under LINUX, on another computer.
Yes, Matlab calculates all the distances, it does not take too long, even for 9000 particles. For one output file just one or two minutes are required. But my Matlab function is still a draft... I am not taking periodic boundaries into account. I didn't go on with the postprocessing implementation because of this big problem of missing contacts.
However, implementing a small neighbor search in Matlab, that reduce the computation of the distances, it shouldn't be a big issue.
Joshua, thanks for your advice. Feel free to send me a message if you think we can exchange some kind of information.
Cheers,
Giuseppe
richti83 | Sat, 01/17/2015 - 18:20
Some Suggestions
Why won't you use compute_pair_gran_local and output it as local dump file. I think you are interested in the stress of each contact, so I would dump id, pos, force, contact area:
compute fc all pair/gran/local id pos force contactArea
dump forcechain all local 1000 post/fc*.dump 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] #x1 y1 z1 x2 y2 z2 id1 id2 periodic fx fy fz area
Than you can calculate the distance from pos2-pos1, evaluate stress by F/A or do some other things. And don't need to find the contacts in extern program again.
To increase precission of output you can change the output format with dump modify this way (this is an example for dump forcechain command, for dump costum the number of columns and datatypes is different as you know):
dump_modify forcechain format "%.12e %.12e %.12e %.12e %.12e %.12e %1.0f %1.0f %1.0f %.12e %.12e %.12e %.12e" #x1 y1 z1 x2 y2 z2 id1 id2 periodic fx fy fz Area
This dumps 12 digits after the dot in scintific notation.
Have a nice weekend,
Christian.
Btw: it is not a big problem to compile liggghts for windows using mpich for windows and cygwin or visual studio
giuraso | Mon, 01/19/2015 - 11:57
Thanks Christian for your
Thanks Christian for your suggestion. It’s really helpful and actually it’s one of the next step I am interested in.
Additional info (my input file build a box with solid fraction of 0.66 – In the case of particle generated by Liggghts, the solid fraction is 0.27)
-“And don't need to find the contacts in extern program again.”
The problem is that using an external computation to verify my simulation in Liggghts, something is wrong. By looking at this calculation, I realized that contacts are missing and this is verified only when I use an input file to insert particles.
Using Joshua’s advice and have as output the number of contacts for each particle:
compute 1 all contact/atom
In dump file there is also the number of contacts for each particle. If I compute the overlap from the particle position (using Matlab), many contacts are missing.
I am also computing PE and KE. KE is fine but Potential Energy is always 0.
compute 2 all pe/atom
compute 3 all ke/atom
And then I would like to compute the total PE as follows:
compute pp all reduce sum c_2
I tried also:
compute 4 all pair/gran/local contactArea
compute compE all reduce sum (c_4*1.0e5) #kn = 1.0e5
variable totalPE equal c_compE
but this is not working
richti83 | Mon, 01/19/2015 - 14:10
external program
I think it is dangerous to export/import data.
If I compute the overlap from the particle position (using Matlab), many contacts are missing
What precission can you handle (eg. 32bit windows vs 64bit dumpfile ?) Some overlaps will be very small.
I would 1st dump pair/gran/local, import the pairs (column 7 and 8) and search this atoms in dump custom file, calculate the overlap, when overlap is zero print id1 id2 position1 position2 and try to find the bug manualy)
You can try to use my extended compute/pair/gran with output of radi and radj and distance calculated internal in liggghts, than you can avoid using the dump custom file.
(discussed here: http://www.cfdem.com/comment/14698#comment-14698)
Sorry I've no Idea about the Energy, maybe you need a special compute or it is simply not supported for pair gran.
Best,
Christian.