Calculating stress-strain/force-displacement of granular solid

LAMMPS LAMMPS Mailing List Mirror
1

Dear LAMMPS-users,

I am trying to model a solid cube based on the discrete element method (DEM). To do this, I am using pair_style from the granular package. In my model, when particles are in contact they interact and results are given following the formula in lammps-docs pair_style gran/… When compressing the cube in one direction, particles move apart in the other directions due to Poisson effect. For the cube to remain a solid, these particles still have to be bonded following DEM interactions. Now, I have bonded these particles using bond_style thanks to the help of Axel.

Particles now can move away from eachother, while still being bonded. However, the forces/stresses that are calculated are now due to bond_style, which don’t include equations related to DEM.

Does anybody know if and how stresses/forces can be calculated related to pair_style gran, while particles are moving away from eachother?

I imagined that if you could assign a radius wherein the fromula in pair_style gran/… is solved, interactions could be calculated. Thought that adjusting the cutoff would be an approach. However, I cannot find anything in pair_style gran where it is possible to specify the cutoff. Could this work and is there a way to assign cutoff?

Thanks a lot for the input in advance. If there just is not a way with LAMMPS to do this, I would also gladly like to know.

Below you can find my LAMMPS-script:

#------------------------------------Initialization--------------------------------------------------------------------------------------------

atom_style hybrid sphere bond
boundary p p p

newton off
comm_modify vel yes

------------------------------------Make geometry--------------------------------------------------------------------------------------------

lattice sc 0.5
region reg block 0 10 0 10 0 10 units box
create_box 1 reg bond/types 1 extra/bond/per/atom 6
create_atoms 1 box

neighbor 1 bin
neigh_modify delay 5

#-------------------------------------Define binding criteria----------------------------------------------------------------------------------
special_bonds lj/coul 1.0 1.0 1.0 extra 6
bond_style harmonic
bond_coeff 1 80.0 0.5
create_bonds all all 1 0.5 1.5

pair_style gran/hooke/history 561e6 241e6 50 NULL 5 0
pair_coeff * *

2

When you use the pair styles granular, you do specify

a radius (or diameter) for every particle. Pairs of

particles only interact when they overlap. Those forces

will contribute to the pressure. All of this is independent

of whether 2 particles are bonded. If you don’t use special_bonds

0 * * to turn off the pairwise interaction between 2 bonded particles,

then those forces and stress contributions will still be present.

If you want 2 granular particles to “interact” when they are not

overlapping, e.g. at some longer cutoff, then you need

a different form of granular pair style, i.e. add a new one to LAMMPS.

Steve

3

Aksel,

This sound like your model is quickly getting needlessly complicated, IMO. It looks like you just have to lose some memory in your line of reasoning, and take a step back. Why use granular models to begin with?

The line of reasoning that I’ve read has been as follows:

I want to observe the Poisson effect in elastic solids by modeling them with the discrete element method (DEM), i.e. granular particles
DEM models only have repulsive conservative forces
Due to a lack of attractive forces non-confined solid masses of granular particles are fragile and undergo flow/explosion or in the best case plastic (non-affine) rearrangement with very small imposed deformations
Introducing bonds can keep the solid mass together in a non-confined set-up

However, your bond style harmonic already is a quasi-generalized version of the Hookean granular model anyway! (with no cut-off at r_0=radius of the granular particle) Why not just use that? There doesn’t seem to be a need for granular models at all.

4

Dear Steve and Eric,

Thank you for your reactions. To be honest, I am not an expert on computational methods for solid mechanics other than the finite element method (FEM). I understand that the approach (DEM) seems cumbersome and it maybe is. It is an approach of my research group to set up a simulation for bone mechanics other then FEM, since the latter solves complexer equations, is computational costly and also less easy to use for nonlinear solid mechanics.
(A part of my project is to compare linear with nonlinear solid mechanics, using a method based on simple mechanics, such as those describing the interactions of DEM)

Using bond style harmonic and nonlinear might indeed also be an approach. However, to model a true solid like bone, shear forces should also be able to be applied. DEM/granular interactions seem more applicable because of containing a tangential component.
A study by Tavarez and Plesha (2007) showing a bonded contact between elements in order to model a solid based on DEM, gave me the confirmation that the concept can be used, hence the motivation for trying it. However, it might be that another framework, maybe LIGGGHTS is more appropriate, or that DEM in general might not be the way to go. I guess it is also a part of my project to find this out.

Best,
Aksel

Tavarez, Federico A., and Michael E. Plesha. “Discrete element method for modelling solid and particulate materials.” International Journal for Numerical Methods in Engineering 70.4 (2007): 379-404.

5

Great! That definitely gives more perspective. You might also turn off the dynamic friction coefficient all together, just using the tangential spring and damper models; the coulumb model could be a source of non-elastic behavior if you hit your large slip threshold. But again I don’t know exactly what you have in mind.

LIGGGHTS is based on the LAMMPS framework, so you probably won’t see much advantage to using that framework.

To your original question, Steve already gave a good answer about the output pressure (read over the compute doc) - and you can use just compute pressure in this case. I would still prefer a granular bond model (where there is a single normal spring and not a messy superposition), but you are using what’s available. Good luck.

6

Dear Eric,
Again thanks for your reply and tips. My supervisor and me discussed, that for the moment, we will focus on stress-strain relation in compression (z-) direction following pair style granular.Because I am not completely sure what will give me the stress in this direction for the entire system, I let the simulation both the compute per atom stress by the stress/atom style command and compute a “-pzz”.

  1. I was wondering what the difference is between the “c_p[3]” component of the per atom stress and “-pzz”. They give exactly the same number, but c_p[3] is a factor 10^3 greater then -pzz.

I also computed the total system pressure with “variable press …”.

I have to verify my model with a analytical solution, for which I am using a basic constitutive relation (isotropic linear elastic, stress strain relation).
2. Just to be sure, what output would you suggest to use for verifying a stress-strain relation, c_p[3], -pzz or the total-system pressure?

  1. In pair style granular doc, it is clearly explained how to get the normal (Kn) and tangential elastic constants (Kt) from a Youngs modulus (E) and Poisson’s ratio (v) for the Hertz-Midlin interaction. Could you also calculate Kn and Kt in a similar way for the hookean interaction? What would be the relation then between Kn or Kt and E and v?

Best,
Aksel

1 Like
7

Hi Aksel,

My responses below (note that I don’t want to get wrapped up in someone else’s research, so unless you have additional LAMMPS related Q’s, which most of these are not, this will probably be my last answer)

  1. From your script, these are related since you used compute reduce sum - and should give the same value if they are normalized correctly. I suggest reading up on the compute stress/atom and pressure documentation - and search past postings if need be. Looking at the dimensions of the two quantities would also give an answer. There has been a lot of chatter on this subject on the mail list, so you might also look there. You may find that stress/atom (without compute reduce and proper normalization) may be more useful as you continue your studies.

  2. Either give the same if normalized correctly - ideally you would make use of the whole stress or pressure tensor rather than a single component. If you don’t know what that is do some reading in any continuum mechanics book - I used Malvern’s book back in the day, though there are probably better/more accessible ones out there.

  3. Note that Hookean contacts are technically impossible for linear elastic spheres. However, they are usually good enough to give physically accurate/reasonable results in granular flows. There are some methods for estimating Hookean model constants depending on what you want to keep constant, e.g. characteristic overlap or energy loss from a collision. That’s for you to decide.

I would not focus on this point however. The discrete elements in your simulation should not be confused with real spheres, they are a coarse-grained elements that you hope will act like an FEA. As such you will also find that the Poisson ratio for these spheres is unlikely to be the same as the macroscale Poisson’s ratio that you obtain from your simulations. Though the relationship between pair force coefficients and observed constitutive relationship may be of interest.
HTH.

1 Like
8

Dear Eric,
Again thanks for your reply and tips. My supervisor and me discussed, that
for the moment, we will focus on stress-strain relation in compression (z-)
direction following pair style granular.
Because I am not completely sure what will give me the stress in this
direction for the entire system, I let the simulation both the compute per
atom stress by the stress/atom style command and compute a "-pzz".

1. I was wondering what the difference is between the "c_p[3]" component
of the per atom stress and "-pzz". They give exactly the same number,
but c_p[3] is a factor 10^3 greater then -pzz.

​pressure != virial stress. the scaling factor is the inverse volume. is
your box by any chance 10x10x10?​ :wink:

I also computed the total system pressure with "variable press ...".

​you should be able to infer the answer to your question 1 from this
formula directly. total pressure is the average of the diagonal elements of
the pressure tensor, that leaves the division by the volume...

axel.​

1 Like