Calculation of pressure in Cylinder by fix mesh/gran/stressanalysis

Submitted by Silias on Thu, 06/14/2012 - 11:49

Hi everybody,

I'd like to calculate the hydrostatic pressure in a cylinder filled with granular particles on its cylindric wall.
I think the command

fix xxx all mesh/gran/stressanalysis xxx.stl ...

is not an appropriate measure to calculate these forces, because by outputting the forces on the mesh with f_xxx[1], f_xxx[2], f_xxx[3] I get the total forces in these directions. That means I get the sum of positive and negative forces and the result will be somewhere near zero...

How could I proceed instead?

Thanks in advance,
kind regards,

Sebastian

Forums:

moritzhoefert | Thu, 06/14/2012 - 17:41

Hi Sebastian,

f_xxx[i] return the resulting forces. If you want to analyze local quantities you have to look at the forces exerted on individual triangles that the wall is composed of. Probably you need a time average but I don't know out of the top of my head how to do that.

fix meshBottom all mesh/gran/stressanalysis ../bottom.stl 2 0.001 0. 0. 191. 0. 0. 0.
fix walls all wall/gran/hooke/history 1 0 mesh/gran 4 meshPole meshPaddle1 meshVessel meshBottom
dump dumpstl4 all mesh/gran/VTK ${dumpMod} post/externals*.vtk stress meshBottom

I'd appreciate if you kept me posted about your progress on that because I might have to do something similar in the future.

Cheers,
Moritz

zamir | Thu, 06/14/2012 - 19:59

Have a look at my blog http://coolsimulations.wordpress.com and scroll down to the post titled "The Aster Template." In it, you will find a python function that grabs a vtk output file from Liggghts and returns a 2d list with coordinates and pressures that you can further manipulate. You can loop the function to grab all vtk files in a folder so you can get a time history of pressures, and you can also average them to get an equivalent confining pressure vs. time.

What is your application? I'm curious...

Zamir

Silias | Fri, 06/15/2012 - 12:43

Hi Zamir,
thanks for the hint with the vtk-files. Since I don't know python well, I'd like to find another possibility to evaluate the pressure on the cylindric walls.
Do you think it's possible to evaluate these files directly in Paraview to get for example a time history of the pressure of different triangles?

The application is an uniaxial compression experiment I'd like to simulate.

Cheers,
Sebastian

ckloss's picture

ckloss | Sat, 06/16/2012 - 11:17

Hi Sebastian,

you can also output the per-particle forces that the wall exerts, see fix wall/gran doc for details. If you post-process that slice-wise, you should be able to reproduce the Janssen curve
Cheers, Christoph

Silias | Mon, 06/18/2012 - 17:50

Hi Christoph,
thank you very much - the Janssen curve is also exactly what I want to reproduce!

Sorry, but I don't understand what you mean with "slice-wise"-postprocessing?
I reckon, that I have to use this command

fix cylwall all wall/gran/hertz/history 1 0 zcylinder 0.05 1 store_force yes

to store the per-particle-forces in a fix property/atom with id force_(ID).
To get the effective force of each atom in the xy-plane I use this command:
variable xyforce atom (f_force_cylwall[1]^2+f_force_cylwall[2]^2)^0.5

My plan would have been to spatially average these effective per-particle force-vectors within diffrent slices in the z-direction...
is that the right method?

My problem here is, that I dont get peratom values when I dump xyforce... xyforce is a scalar or lets say every atom has the same value... Do you have a clue, how that comes? Because dumping f_force_cylwall[1] results in diffrent values for each particle...

Thanks in advance,
kind regards,

Sebastian

ckloss's picture

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

Hi Sebastian,

>>My plan would have been to spatially average these effective per-particle force-vectors within diffrent slices in the z-direction...
>>is that the right method?
That's what I meant

>>xyforce is a scalar or lets say every atom has the same value... Do you have a clue, how that comes?
Actually not... does it work at all to output per-atom values?

Cheers, Christoph

Silias | Tue, 06/19/2012 - 10:53

Yes, I can dump xyforce and f_force_cylwall[1]. Dumping v_xyforce every particle has the same value while dumping f_force_cylwall[1] every particle has it's own value...

Something else, which is funny, is with

variable a atom f_force_cylwall[1]

dumping v_a results in the same values for every particle while dumping f_force_cylwall[1] results in different values for all particles...

I use version 1.5.3

Cheers,
Sebastian

ckloss_ | Sat, 06/23/2012 - 17:45

Hi Sebastian,

difficult to say if this is a problem in your script or an internal issue. version 2.0RC does not have the stressanalysis feature yet, but this is about to come in weeks - then you can try if 2.0 solves your problem. If your problem persists, I will have a closer look.

In the meantime, you could dump directly (w/o the variable command) and do the evaluation as post-procesing

Cheers, Christoph