# Create packing with given void ratio at different isotropic pressures

Hello everyone,
I am Antrox and a new member of this forum.
I just started using YADE and more generally DEM, and I would like some clarifications to my doubts if possible.
I created a periodic cell and I am using a regular orthogonal packing by using a given particle radius, which I do not want to change.
By applying O.cell.velGrad=Matrix3(-rate,0,0,0,-rate,0,0,0,-rate), I am supposed to generate an isotropic pressure which should not change the void ratio of the cell. Indeed, I computed void ratio = 0.91 which should be close to the analytical solution for spheres (0.87).
1) therefore, is it correct to affirm that I could have a regular packing with any isotropic pressure I decided ?
2) then, after reaching the target isotropic pressure,, I tried to apply an ideal shearing by O.cell.velGrad=Matrix3(0,rate,rate,rate,0,rate,rate,rate,0), imaging that this deviatoric component will produce compaction without changing the isotropic component. But I obtained the opposite effect, the void ratio did not change, whilst the mean pressure did increase. Why so?
3) if so, I would like to obtain the desired void ratio at a given isotropic pressure without changing the particle radius. What is the best approach?

## Question information

Language:
English Edit question
Status:
Solved
For:
Assignee:
No assignee Edit question
Solved by:
Jan Stránský
Solved:
Last query:
 Revision history for this message Jan Stránský (honzik) said on 2022-04-11: #1

Hello,

> I am ... a new member of this forum.
> I just started using YADE and more generally DEM

welcome :-)

> I am supposed to generate an isotropic pressure which should not change the void ratio of the cell.

it depends on definition of void ratio or definition of values from which it is computed.
If you apply pressure / velocity gradient, the total volume of the cell is changing, so intuitively also void ratio would change.

> I computed void ratio = 0.91

how?
A MWE (see below) would help

> 1) therefore, is it correct to affirm that I could have a regular packing with any isotropic pressure I decided ?

yes, you can have a regular packing with any (reasonable) isotropic pressure.
Reasonable = e.g. with linear contact law (e.g. Law2_ScGeom_FrictPhys_CundallStrack) the pressure is limited by stiffness (with higher pressures resulting in degenerated collapsed "one-sphere" packing)

> 2) then, after reaching the target isotropic pressure,, I tried to apply an ideal shearing by O.cell.velGrad=Matrix3(0,rate,rate,rate,0,rate,rate,rate,0), imaging that this deviatoric component will produce compaction without changing the isotropic component. But I obtained the opposite effect, the void ratio did not change, whilst the mean pressure did increase. Why so?

ideal shearing ... compaction ... isn't it contradicting?
Void ratio does not change, because total volume does not change (because of ideal shearing).
To help you why the pressure increased, we would need a MWE (see below).

> 3) if so, I would like to obtain the desired void ratio at a given isotropic pressure without changing the particle radius. What is the best approach?

Depending on definition, but I feel like one of the values is dependent on the other.
I.e. you cannot independently vary pressure and void ratio.

Cheers
Jan

PS: also please read [1] and consider point 5 and consider providing a MWE.

 Revision history for this message Gianni Pellegrini (antrox) said on 2022-04-11 #2

Thanks Jan,
I calculated the void ratio via :
u=utils.porosity()
e=u/(1-u)

I am also checking that the overlapping volume is small enough to avoid what you mentioned above through
ds = [2*r - i.geom.penetrationDepth for i in O.interactions]
volContactss = [1/12.*pi*(4*r+d)*(2*r-d)**2 for d in ds]
volContacts = sum(volContactss)
To check that the overlapping volume is smaller than 0.1%.
But this has been removed from the MWE below because it is not interesting now.

Briefly, I am reaching my target pressure by O.cell.velGrad=Matrix3(-rate,0,0,0,-rate,0,0,0,-rate) which cannot change the void ratio of a regular ortho packing.
Then I am applying a axial load to modify the void ratio
As you point out, a pure shearing through O.cell.velGrad=Matrix3(0,rate,rate,rate,0,rate,rate,rate,0) will not produce any change of volume (as expected, my bad)
So, I am just now pushing one wall by O.cell.velGrad=Matrix3(-rate,0,0,0,0,0,0,0,0), reaching the target void ratio.
In the MWE below, I reach 0.89 void ratio at 1E5 and then the target 0.7 at 1.8e7

Now, since I am using just a linear model for the normal interaction of two bodies, my idea is I could use either 2 approaches:
1) modify the elasticity of the material (the variable young) proportionally to have the desired pressure :
factor=IsoGoal/(getStress().trace() / 3.)
O.materials[0].young =En*factor
O.interactions.clear()

2) modify the branch vector to release the pressure to the target pressure. In practical terms, it is like having a rigid link and therefore reducing the flexible length between the 2 spheres without changing their distance. Something I read in the forum, it can be simulated through
for i in O.interactions:
i.phys.unp = i.geom.penetrationDepth

Now, in the MWE below, I tested 1) .
I am not sure now about the implications of doing that. I am achieving to change to have a sample with desired pressure and void ratio. Please consider that being a MWE, I cut the final part and hence, at every loop is updating back the old stiffness.
What do you think about this approach? I can see the drawback is the change of the contact stiffness and indeed I would like to explore the approach 2) to avoid that.

MWE:

import time
import datetime
import os

from yade import pack, plot, export
import numpy as np

#PARAMETERS
frictionAngle = 0.6
VoidRatio = 0.7
IsoGoal=-100000
poisson=0.2
R=1e-3
rate=1e-4
dimcell = 0.02
density= 1e12
En=1e9

#SETTINGS
O.periodic = True
O.cell.hSize = Matrix3(dimcell , 0, 0, 0, dimcell , 0, 0, 0, dimcell )

#ENGINE
pp=O.materials.append(FrictMat(young=En,poisson=poisson,frictionAngle=frictionAngle,density=density))
O.bodies.append(pack.regularOrtho(pack.inAlignedBox((0, 0, 0), (dimcell , dimcell , dimcell)), radius=R, gap=0, color=(0, 0, 1), material=pp))

O.engines = [
ForceResetter(
),
InsertionSortCollider([Bo1_Sphere_Aabb()]),
InteractionLoop([Ig2_Sphere_Sphere_ScGeom()], [Ip2_FrictMat_FrictMat_FrictPhys()], [Law2_ScGeom_FrictPhys_CundallStrack()]),
NewtonIntegrator(damping=.2),
PyRunner(command='Compaction()', realPeriod=1)
]

O.dt = .5 * PWaveTimeStep()

def Compaction():
u=utils.porosity()
e=u/(1-u)
print('e', e)
if np.abs(getStress().trace() / 3.) > np.abs(IsoGoal):
print('compaction done', e)
if e < VoidRatio:
print('e', e)
print('shearing done',getStress().trace() / 3)
factor=IsoGoal/(getStress().trace() / 3.)
O.materials[0].young =En*factor
O.interactions.clear()

O.run() # run forever

Thank you

 Revision history for this message Jan Stránský (honzik) said on 2022-04-11: #3

Hello,

> But this has been removed from the MWE below because it is not interesting now.

+1 for understanding the meaning of MWE (which is not a matter of course at all)

> Briefly, I am reaching my target pressure by O.cell.velGrad=Matrix3(-rate,0,0,0,-rate,0,0,0,-rate) which cannot change the void ratio of a regular ortho packing.

But isotropic pressure (applied as constant isotropic velGrad) DOES change the void ratio (because it changes total volume).
Just see the output of the script (where "e" is decreasing).
Am I missing something?

> So, I am just now pushing one wall by O.cell.velGrad=Matrix3(-rate,0,0,0,0,0,0,0,0), reaching the target void ratio.

Why not 3 walls? (related to above (?))

> 1) modify the elasticity of the material

yes, introducing another variable, you can now vary pressure AND void ratio independently and stiffness is now the dependent variable.

> O.materials[0].young =En*factor

be aware that modifying material has no effect on existing interactions (only on the newly created)

> 2) modify the branch vector ...
> ...
> I tested 1)
> ...

I like 1) more.
One reason is that I have almost no experience with 2 :-)
But more importantly I **think** like the two approaches are very similar - modifying material to introduce another variable to let void ratio and pressure be independent.

> I can see the drawback is the change of the contact stiffness and indeed I would like to explore the approach 2) to avoid that.

Yes, there are drawbacks. E.g. that you define material to get good result, but usually the material properties (stiffness is usually the most basic) are the inputs.
What exactly is meant by you?

Cheers
Jan

 Revision history for this message Gianni Pellegrini (antrox) said on 2022-04-13: #4

Hello Jan,
yes, you are right, it is better to contextualize.
I would like to simulate triaxial testing of samples with different void ratios.

the first drop of void ratio you observed is just because the spheres are floating in the air (no contact), therefore for my case, this is not a realistic scenario. The analysis starts when the void ratio is about 0.89 as analytically calculated and the mean pressure becomes different from 0.
After this point, if the packing is regular and the loading is isotropic, the void ratio should stay constant because the configuration of the spheres cannot change (if pressure is reasonable).

So far, method 1 works and I think that by using O.interactions.clear() , the new forces in the following loop are calculated correctly since the interactions are continuously computed (so there is no new or old). Maybe I am wrong but your answer solved my problem.
I will investigate method 2 with a new question in the future.
Thanks again

 Revision history for this message Gianni Pellegrini (antrox) said on 2022-04-13: #5

Thanks Jan Stránský, that solved my question.

 Revision history for this message Jan Stránský (honzik) said on 2022-04-14: #6

> After this point, if the packing is regular and the loading is isotropic, the void ratio should stay constant because the configuration of the spheres cannot change (if pressure is reasonable).

again, it depends on the definition of "constant", "cannot change" etc.
Because with increasing pressure and decreasing overall volume decreasing the void ration is NOT constant and the spheres DO change (although the effect may by negligible).

> ... method 1 works and I think that by using O.interactions.clear() , the new forces in the following loop are calculated correctly

yes, after O.interactions.clear() new interactions are created. New interactions are based on current materials.

Cheers
Jan