Particle Deposition structure (help)

Submitted by mahdi_roozbeh on Wed, 11/23/2011 - 05:42

Hi Friends:

I almost read the documentation of the LIGGGHTS and I'm somehow getting familiar with the commands. I would be thankful if you can help me in fix deposit command.

I followed the pour_dev example to create my own input script, but the output are unacceptable.

What I'm going to do is to pour 1000 particles in cylindrical container (we have boundary interaction with particles) ; and in every pre-specified step, I want to calculate the void ratio (porosity) of deposited particles.

My question is that Can I put fix deposit command in pour_dev example ? If I can, what parts should be changed to make this input script rational ?

TX

ckloss's picture

ckloss | Thu, 11/24/2011 - 08:22

Hi,

not sure if I got your question? What do you mean with "void ratio of deposited particles"?
Yes, you should be able to use fix deposit with any script

>>If I can, what parts should be changed
That's something you will have to find out

Have fun, Christoph

mahdi_roozbeh | Thu, 11/24/2011 - 10:27

Thanks Christoph

I have used pour/legacy command ( or it can be pour deposit).
void ratio = porosity = empty regions of the container filled by particles

I have some questions and I would thankful if you can answer them.

1)The problem is that when I' using for example 1000 particles to be inserted in the defined region, I'm seeing just a fraction of those particles in Paraview not all of them. I dont know what is the problem and also I have used different radius particles, while the Glyph of the Paraview doesnt scale them properly (particles are shown bigger than they are). How can I scale them properly?

2) Is there any command that I can find the porosity in a specified height of the container ( from the bottom of the container until the specified height)?

3) Based on the DEM, is there any overlaps between particles in pour commands at the end of implementation of the script?

TX

raguelmoon's picture

raguelmoon | Mon, 11/28/2011 - 12:31

Hi TX,
You dont need to scale if you are using new version of pizza available at this site. After converting dump file into vtk, open them in paraview and do this in Glyph:
1. in "Glyph Type" select "Sphere"
2. in "Sphere", set radius to 1, don't change anything
3. set "Scale Mode" to Scalar
4. change "Set Scale Factor" value to 1 (first check Edit option right to it)
5. Uncheck both options at bottom "Mask points" and Random Mode"
6. Now go to the top of Glyph and set "Scalars" value as radius

You can see all particles in correct dimensions.
Cheers,
Ram

Ram

Maryam | Mon, 06/18/2012 - 02:20

So, is it possible to output the porosity of a packing or the total grain volume? or it should be calculated by post-processing the dump file?

Maryam

ckloss's picture

ckloss | Mon, 06/18/2012 - 19:11

There are multiple ways to do it:

(a) write your own post-processing for that
(b) use the CFDEMcoupling from this website to evaluate the proposity based on a CFD mesh
(c) use a third-aprty tool (like voro++) for that - just search the DEM discussion boards for "voro++"

Christoph

JoshuaP | Fri, 01/16/2015 - 10:34

Hi,

I'm trying also to compute the void ratio and this was the only topic I found. I tried it with (c) the voronoi volume, so I calculate the volume of each particle divided by his voronoi volume. That takes a lot of computation time and in the end I dont get my assumed value. So I tried (a) and calculate the volumes of all particles and sum them. Then I divide it by the simulation volume. This is fine for low stress states, where not many overlaps occur.
So for a good estimation of the void ratio I need the overlap, maybe through the potential energy, but this is not just calculated by the overlap. Is there anyway to calculate the overlapvolume? Or maybe direktly the void area?

Thanks
Joshua

NTT1508 | Sun, 01/18/2015 - 03:44

Hi Joshua,

Can we calculate overlap based on contact area and position of dump/gran/local? Contact area is calculated with reference to overlap volume of particles I think.

Regards,

Nathan,

richti83's picture

richti83 | Sun, 01/18/2015 - 11:19

contact area is computed this way:

radi = atom->radius[i];
radj = atom->radius[j];
vectorSubtract3D(atom->x[i],atom->x[j],del);
rsq = vectorMag3DSquared(del);
r = sqrt(rsq);
contactArea = - M_PI/4 * ( (r-radi-radj)*(r+radi-radj)*(r-radi+radj)*(r+radi+radj) )/rsq;
array[ipair][n++] = contactArea;

I see no way to eleminate radi and radj. But I extended compute_pair_gran local to dump ri, rj, |del| and overlap.
see this commit for the changes (or just pull from richti83 repository)
https://github.com/richti83/LIGGGHTS-WITH-BONDS/commit/dd9ee0523ca053b16...

If someone can tell me how to calculate the common volume of two overlapping spheres I can add this too.

When my formulas are correct maybe Christoph can add these few lines to LIGGGHTS-PUBLIC master.

Best,
Christian.

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

NTT1508 | Sun, 01/18/2015 - 13:08

Hi Christian,

It is more convenient to find out overlap volume with respect to radius of Ri, Rj (that you have extracted) and distance (d) between them (based on there positions thanks to compute/pair/gran/local). Refer to some basic mathematical sources (ex wolfram...) which is quite easy to prove, the overlap volume between two particles can be calculated as:

V= PI*(Ri+Rj-d)^2*[d^2+2d*(Ri+Rj)-3*(Ri-Rj)^2]*/(12*d)

If Ri+Rj=d, V return to zero.

In terms of overlap area, we can get radius of overlap disk (circle) and then it becomes the same way as distance based solution.

Thanks,

Regards,

Nathan,

richti83's picture

richti83 | Sun, 01/18/2015 - 15:23

Ok, I understand.

Now as we have the overlap volume, we could calculate the "real" void fraction by subracting all overlap-volumes from particle volume ?!
I see two ways:
1.) a new property/atom which holds the adjusted volume (V_Sphere-Sum(V_overlap)), I'm not sure if this is multicore save when the contacts are on different domains
2.) do it in ParaView even in the forcechain reader or with a python programmable filter, in this case it does not costs any simulation time and there are no multicore problems, but it is more difficult to get the contribution of particle-particle and particle-wall contacts. (maybe the filter needs two inputs, holds an array with id->volume and substract for each line of the reader-output the overlaping volume...)
I would prefer the 2nd choice because it is more flexible and can be used in matlab or other postprozessing tools too. I will cook up a prototype.

cheer,
Christian.

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

NTT1508 | Wed, 01/21/2015 - 09:45

Hi Christian,

Yes I think second one is more convenient and helpful in the future use.

Anyway, what is the new syntax now for compute pair/gran/local on your modified version to extract Ri,j and distance while ramaining pos,id,contactA,force and torque.

Thanks,

Regards,

richti83's picture

richti83 | Wed, 01/21/2015 - 11:14

Hi Nathan,

The compute command is the same, I extended "contactArea" by four values.
To dump add c_fwc[13] c_fwc[14] c_fwc[15] c_fwc[16] c_fwc[17] to your dump local command.
At the moment I'm struggeling with spheres which overlap more than one other sphere because than we need to calculate the overlapping value of three spheres.
I tried to proof my code with CAD Data, see this comparisation for the results:
http://www.richtisoft.de/transfer/overlapping_spheres.pdf
As you can see at the moment the corrected volume calculated by my filter is to low, I think I need to decide if the overlap is allready substracted by a prior contact.
Testssript to generate these arrangements:
http://www.richtisoft.de/transfer/in_void.lmp
Statefile for ParaView with programmable filter included:
http://www.richtisoft.de/transfer/voidratio3.pvsm
Needed reader code (not officall yet because may be changed in near future), use build/libliggghts_forcechainreader.so when you won't like to recompile.
http://www.richtisoft.de/transfer/forcechain_reader.tar.gz

Let me know when you find my mistake :-)

Best,
Christian

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

NTT1508 | Wed, 01/21/2015 - 23:33

Hi Christian,

Thank for your kind and hard work,

I think the problem of your code is possible due to common volume of 3 or more spheres that is subtracted more than one time. For example if we have (1), (2) or (3) particles, the overlap volume of (1) and (2) might be partially in common with the overlap volume of [(2),(3)] and [(1),(3)]. So a part of overlap volume is subtracted more than one time when we calculate overlap volume of a pair separately. This probably causes your result smaller than CAD calculation.

We can build a library for several cases from (2) to ( 6) spheres sharing common volume and recall them for each suitable situation but I think it is really time consuming. There are also other algorithms to deal with overlap volume of objects but definitely more complicated. In PFC, there's also a command to estimate overlap volume of particles If I remember correctly.

Thank Christian very much,

Best regards,

Nathan,

JoshuaP | Tue, 03/03/2015 - 11:49

hey,

I found a way to calculate void ratio without calculating the overlap volume. void ratio is equal to particle density divided by the bulk density minus 1.
particle density is known, total mass can be calculated and is always the same even after compression, bulk volume is known.


compute rad all property/atom radius
variable perpartvolume atom c_rad*c_rad*c_rad*PI*4/3
compute partvolume all reduce sum v_perpartvolume
variable bulkdensity equal ${rho}*c_partvolume/${simvolume}
variable voidratio equal ${rho}/${bulkdensity}-1

If someone finds a mistake in this please let me know.

regards
Joshua

NTT1508 | Wed, 03/04/2015 - 00:30

Hi Joshua,

Thank for your interest. Theoretically, it is a good way to calculate void ratio through bulk density.

But I cannt understand why void ratio=[(particle density)/(bulk density)-1].? In my understanding, Void Ratio=Vv/Vp=(Vcell-Vp)/Vp=Vcell/Vp-1=[(m_fluid/m_par)-1][particle density/bulk density] - 1=(w+1)(particle density/bulk density)-1. here w can be water content in Geo_mechanics. Void fraction, on the other hand, = 1-Vvoid/Vcell.
total mass cannot be constant because when overlap occurs, total volume of particles decreased and results in reduction of total mass. Regardless overlap, it is simple case.

Anyway, Christian has created a powerful technique on Paraview with reference to Monte Carlo concept to calculate overlap volume. I think that technique is reliable.

Regards,

richti83's picture

richti83 | Wed, 03/04/2015 - 09:04

total volume of particles decreased and results in reduction of total mass
I think this is not true. The mass of two overlapping spheres should be equal to the two spheres without overlap. In reality the material density at contact point will increase localy under compression (volume is reduced, mass is constant -> rho increases). Think of a accident with a car, the mass is the same, but the body will be smaller.
That's why I think you can not mix mass and volume in this case, because the mass of two totaly overlapping particles will remain 2*mass. But the volume will be 4/3Pi*r^3.

In the above case partvolume would be 2*Vpart , simvolume would be r^3 (asuming your variable simvolume represents the max extends of the bounding box) and voidratio=simvolume/(2*partvolume)-1 = (1/(2*4/3 Pi)) -1= const. (r^3 is cutted in denominator and numerator)

When it would be so easy, I would not use the time consuming monte carlo aproach :-).

Best,
Christian.

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

NTT1508 | Wed, 03/04/2015 - 14:02

Hi,

It is interesting to hear from Christian. You are right, in many practical cases, the density of material will be increased in the contact region. Volume decreases but mass is conserved. In DEM, though we accept overlap that is used for void ratio estimation (volume concern), however the program still considers total mass=mass of part 1 + mass of part 2 ( in case particle size is different). Note that Overlap concept is used for soft particles by Cundall, 1979. In many practice cases, hard material is nearly not overlapped. Also in soft material, overlaps is common. However, I believe in many real cases of soft material, the mass at contact region is not conserved due to heterogeneous of material.

Thank for sharing,

Regards,

JoshuaP | Wed, 03/04/2015 - 14:42

yes I forgot that particle density changes by overlapping of two particles. I'm now trying to install your paraview plugin but I didnt get cmake running yet. Is it possible to install it without cmake?
regards

richti83's picture

richti83 | Wed, 03/04/2015 - 15:27

just use the precompiled version in build folder. (you need to download the coresponding paraview binary from kitware first)

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

JoshuaP | Thu, 03/05/2015 - 11:13

hi Christian,

many thanks for these plugins. Do you have maybe somewhere a description how to use it for plotting force chains and voidratio? I searched for it already in the forum and on your blog but can't find sth.. I mean not how to install, that worked fine with your description from your blog, but how to visualize it in paraview.

many thanks
Joshua

richti83's picture

richti83 | Thu, 03/05/2015 - 13:19

Have a look at the 2nd picture here:
http://techblog.richtisoft.de/?p=41
first append a Calculator filter, name result array name Z and input in the equitation field coordsZ (or select from scalars drop down).
To the calculator append 1..N Threshold filters, set Scalars to Z and range to zLo-zHi for every region you want to extract.
To every threshold appaed the 000_VoidFraction filter you have downloaded from my blog [1].
To create a chart you need to adjust the output array of the filter to Type vtkFloatArray and outputarray.SetValue(0,void), tile the view, select sphreadsheed, choose Showing VoidFractionFilter1 .. N, Attribute row data, select the first line. Than filters->PlotSelectionOverTime -> wait some time and you will get a voidfraction over time plot of the selected region. You can repead this for every regions and show all graphs in one line chart view.

[1] http://www.richtisoft.de/transfer/voidfraction_custom_filter.xml

I'm not an associate of DCS GmbH and not a core developer of LIGGGHTS®
ResearchGate | Contact

saman | Thu, 05/23/2019 - 16:38

Hi everyone
After all this discussion is now the best way to calculate the ratio void? I'm really tired of being here so much looking.
(ratio void =volume of voids /volume of solids)