yade-users team mailing list archive

[Gianni Pellegrini <question701259@xxxxxxxxxxxxxxxxxxxxx>]

Question #701259 on Yade changed:
https://answers.launchpad.net/yade/+question/701259

    Status: Answered => Open

Gianni Pellegrini is still having a problem:
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 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

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.