Create packing with given void ratio at different isotropic pressures

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.

[1] https://www.yade-dem.org/wiki/Howtoask

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:
from yade import pack,plot

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()
O.cell.velGrad=Matrix3(-rate,0,0,0,-rate,0,0,0,-rate)

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)
        O.cell.velGrad=Matrix3(-rate,0,0,0,0, 0, 0, 0,0)
    if e < VoidRatio:
         print('e', e)
         print('shearing done',getStress().trace() / 3)
         O.cell.velGrad=Matrix3(0,0,0,0,0, 0, 0, 0,0)
         factor=IsoGoal/(getStress().trace() / 3.)
         O.materials[0].young =En*factor
         O.interactions.clear()

O.run() # run forever

Thank you

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

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

> 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