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.
_______________________________________________
Mailing list: https://launchpad.net/~yade-users
Post to : yade-users@lists.launchpad.net
Unsubscribe : https://launchpad.net/~yade-users
More help : https://help.launchpad.net/ListHelp
