# encoding: utf-8
from woo.dem import *
import woo.core
import woo.dem
import woo.pyderived
import woo.triangulated
import math
from minieigen import *
import woo.models
[docs]class TriaxTest(woo.core.Preprocessor,woo.pyderived.PyWooObject):
r'''Preprocessor for triaxial test with rigid boundary. The test is run in 2 stages:
* **compaction** where random loose packing is compressed to attain :math:`\sigma_{\rm iso}` (:obj:`sigIso`) in all directions. The compaction finishes when the stress level is sufficiently close to :obj:`sigIso` and unbalanced force drops below :obj:`maxUnbalanced`.
* **Triaxial compression**: displacement-controlled compression along the ``z``-axis, with strain rate increasing until :obj:`maxRates` is reached; the test finished when axial strain attains :obj:`stopStrain`. During this phase, lateral (:math:`x` and :math:`y`) stresses are maintained at :math:`\sigma_{\rm iso}`, as much as possible.
.. youtube:: qWZBCQbS6x4
'''
_classTraits=None
_PAT=woo.pyderived.PyAttrTrait # less typing
def postLoad(self,I):
if self.preCooked:
print('Applying pre-cooked configuration "%s".'%self.preCooked)
if self.preCooked=='Spheres in cylinder':
self.generator=woo.dem.PsdSphereGenerator(psdPts=[(.01,0),(.04,1.)])
self.shape='cylinder'
elif self.preCooked=='Capsules in cylinder':
self.generator=woo.dem.PsdCapsuleGenerator(psdPts=[(.01,0),(.04,1.)],shaftRadiusRatio=(.6,1.2))
self.shape='cylinder'
elif self.preCooked=='Ellipsoids in box':
self.generator=woo.dem.PsdEllipsoidGenerator(psdPts=[(.01,0),(.04,1.)],axisRatio2=(.5,.9),axisRatio3=(.3,.7))
self.shape='box'
elif self.preCooked=='Sphere clumps in box':
self.generator=woo.dem.PsdClumpGenerator(psdPts=[(.01,0),(.04,1.)],clumps=[
SphereClumpGeom(centers=[(0,0,-.15),(0,0,-.05),(0,0,.05)],radii=(.05,.08,.05),scaleProb=[(0,1.)]),
SphereClumpGeom(centers=[(.05,0,0) ,(0,.05,0) ,(0,0,.05)],radii=(.05,.05,.05),scaleProb=[(0,.5)]),
])
self.shape='box'
else: raise RuntimeError('Unknown precooked configuration "%s"'%self.preCooked)
self.preCooked=''
_attrTraits=[
_PAT(str,'preCooked','',noDump=True,noGui=False,startGroup='Predefined config',choice=['','Spheres in cylinder','Capsules in cylinder','Ellipsoids in box','Sphere clumps in box'],triggerPostLoad=True,doc='Apply pre-cooked configuration (i.e. change other parameters); this option is not saved.'),
# noGui would make startGroup being ignored
_PAT(float,'sigIso',-500e3,unit='kPa',startGroup='General',doc='Confining stress (isotropic during compaction)'),
_PAT(Vector3,'maxRates',(2e-1,2e-1,1.),'Maximum strain rate during the compaction phase (for all directions), during the triaxial phase in axial sense, and during the triaxial phase in radial sense(s).'),
_PAT(float,'stopStrain',-.3,unit=r'%',doc='Goal value of axial deformation in the triaxial phase'),
_PAT(bool,'planeStrain',False,doc='(For demonstration purposes only:) during the triaxial phase, prescribe zero displacement along :math:`x` and stress-control only :math:`y`-axis.'),
_PAT(str,'shape','cell',choice=('cell','box','cylinder'),doc='Shape of the volume being compressed; *cell* is rectangular periodic cell, *box* is rectangular :obj:`~woo.dem.Wall`-delimited box, *cylinder* is triangulated cylinder aligned with the :math:`z`-axis'),
_PAT(Vector3,'iniSize',(.3,.3,.6),unit='m',doc='Initial size of the volume; when :obj:`shape` is ``cylinder``, the second (:math:`y`) dimension is ignored.'),
_PAT(woo.dem.ParticleGenerator,'generator',woo.dem.PsdCapsuleGenerator(psdPts=[(.01,0),(.04,1.)],shaftRadiusRatio=(.6,1.2)),
# woo.dem.PsdEllipsoidGenerator(psdPts=[(.01,0),(.04,1.)],axisRatio2=(.5,.9),axisRatio3=(.3,.7))
# woo.dem.PsdSphereGenerator(psdPts=[(.01,0),(.04,1.)])
'Particle generator; partices are then randomly placed in the volume.'),
_PAT(woo.models.ContactModelSelector,'model',woo.models.ContactModelSelector(name='linear',damping=.7,numMat=(1,2),matDesc=['particles','boundary'],mats=[woo.dem.FrictMat(density=1e7,young=1e8)]),doc='Select contact model. The first material is for particles; the second, optional, material, is for the boundary (the first material is used if there is no second one).'),
_PAT(str,'reportFmt',"/tmp/{tid}.xhtml",startGroup="Outputs",doc="Report output format; :obj:`Scene.tags <woo.core.Scene.tags>` can be used."),
_PAT(str,'saveFmt',"/tmp/{tid}-{stage}.bin.gz",'''Savefile format; keys are :obj:`Scene.tags <woo.core.Scene.tags>`; additionally ``{stage}`` will be replaced by
* ``init`` for stress-free but compact cloud,
* ``iso`` after isotropic compaction,
* ``backup-011234`` for regular backups, see :obj:`backupSaveTime`, 'done' at the very end.
'''),
# _PAT(int,'backupSaveTime',1800,doc='How often to save backup of the simulation (0 or negative to disable)'),
_PAT(float,'dtSafety',.7,startGroup='Tunables',doc='See :obj:`woo.core.Scene.dtSafety`.'),
_PAT(float,'maxUnbalanced',.1,'Maximum unbalanced force at the end of compaction'),
_PAT(int,'cylDiv',40,hideIf='self.shape!="cylinder"',doc='Number of segments to approximate the cylinder with.'),
_PAT(float,'massFactor',.2,'Multiply real mass of particles by this number to obtain the :obj:`woo.dem.WeirdTriaxControl.mass` control parameter'),
_PAT(float,'rateStep',.01,'Increase strain rate by this relative amount at the beginning of the triaxial phase, until the value given in :obj:`maxRates` is reached.'),
]
[docs] def __new__(klass,**kw):
self=super().__new__(klass)
self.wooPyInit(TriaxTest,woo.core.Preprocessor,**kw)
return self
def __init__(self,**kw):
woo.core.Preprocessor.__init__(self)
self.wooPyInit(TriaxTest,woo.core.Preprocessor,**kw)
[docs] def __call__(self):
# preprocessor builds the simulation when called
return prepareTriax(self)
[docs]def prepareTriax(pre):
import woo
woo.master.usesApi=10102
S=woo.core.Scene(fields=[DemField()])
S.pre=pre.deepcopy()
S.lab.partMat=pre.model.mats[0]
wallMat=S.lab.partMat if len(pre.model.mats)==1 else pre.model.mats[1]
# initially no friciton
S.lab.partMat.tanPhi=0.
partMask=0b0001
wallMask=0b0011
loneMask=0b0010
S.dem.loneMask=loneMask
margin=.1
factoryKw=dict(maxMass=-1,maxNum=-1,generator=pre.generator,massRate=0,maxAttempts=5000,materials=[S.lab.partMat],shooter=None,mask=partMask,collideExisting=False)
if pre.shape=='cell':
factory=woo.dem.BoxInlet(box=((0,0,0),pre.iniSize),**factoryKw)
S.lab.relVol=1.
margin=0. # no margins, use the full volume
elif pre.shape=='box':
mid=(margin+.5)*pre.iniSize;
for ax in (0,1,2):
for sense in (-1,1):
# make copy of mid
pos=Vector3(mid); pos[ax]-=sense*.5*pre.iniSize[ax]
S.dem.par.add(woo.utils.wall(pos,sense=sense,axis=ax,mat=wallMat,mask=wallMask),nodes=True)
factory=woo.dem.BoxInlet(box=(margin*pre.iniSize,(1+margin)*pre.iniSize),**factoryKw)
S.lab.relVol=AlignedBox3((0,0,0),pre.iniSize).volume()/AlignedBox3((0,0,0),(1+2*margin)*pre.iniSize).volume()
elif pre.shape=='cylinder':
d,h=pre.iniSize[0],pre.iniSize[2]
bot=Vector3((.5+margin)*d,(.5+margin)*d,margin*h)
top=bot+(0,0,h)
factory=woo.dem.CylinderInlet(node=woo.core.Node(pos=bot,ori=Quaternion((0,1,0),-math.pi/2.)),radius=.5*d,height=h,**factoryKw)
# axDiv to make sure we don't have facet larger than half of the cell (collider limitation)
S.dem.par.add(woo.triangulated.cylinder(bot,top,radius=.5*d,div=pre.cylDiv,axDiv=3,capA=True,capB=True,wallCaps=True,mat=wallMat,mask=wallMask),nodes=True)
S.lab.relVol=math.pi*d*h/AlignedBox3((0,0,0),(1+2*margin)*pre.iniSize).volume()
if isinstance(pre.generator,woo.dem.PsdEllipsoidGenerator): raise NotImplementedError('It is not possible to combine ellipsoids with cylindrical boundary, since Ellipsoid+Facet collisions are not (yet) supported.')
else: raise RuntimeError('Triax.shape must be one of "cell", "box" or "cylinder".')
# generate as many particles as possible
# list(woo.system.childClasses(woo.dem.BoundFunctor))
S.engines=[factory]
S.one()
print('Number of nodes:',len(S.dem.nodes))
S.periodic=True
S.cell.setBox((1+2*margin)*pre.iniSize)
S.dtSafety=pre.dtSafety
S.engines=[
woo.dem.WeirdTriaxControl(goal=(pre.sigIso,pre.sigIso,pre.sigIso),maxStrainRate=(pre.maxRates[0],pre.maxRates[0],pre.maxRates[0]),relVol=S.lab.relVol,stressMask=0b01111,globUpdate=1,maxUnbalanced=pre.maxUnbalanced,mass=pre.massFactor*sum([n.dem.mass for n in S.dem.nodes]),doneHook='import woo.pre.triax; woo.pre.triax.compactionDone(S)',label='triax',absStressTol=1e4,relStressTol=1e-2),
woo.core.PyRunner(20,'import woo.pre.cylTriax; woo.pre.triax.addPlotData_checkProgress(S)')
]+woo.utils.defaultEngines(model=pre.model,dynDtPeriod=100)
S.lab.stage='compact'
S.lab._setWritable('stage')
try:
import woo.gl
S.gl(woo.gl.Renderer(dispScale=(5,5,2),rotScale=0,cell=True if pre.shape=='cell' else False),woo.gl.Gl1_DemField(colorBy='radius'),woo.gl.Gl1_CPhys(),woo.gl.Gl1_Wall(div=3))
except ImportError: pass
return S
[docs]def addPlotData_checkProgress(S):
assert S.lab.stage in ('compact','triax')
import woo
t=S.lab.triax
sxx,syy,szz=t.stress.diagonal()
dotE=S.cell.gradV.diagonal()
dotEMax=t.maxStrainRate
# net volume
vol=S.cell.volume*S.lab.relVol
# current radial stress
srr=.5*(sxx+syy)
# mean stress
p=t.stress.diagonal().sum()/3.
# deviatoric stress
q=szz-.5*(sxx+syy)
qDivP=(q/p if p!=0 else float('nan'))
if S.lab.stage=='compact':
## t.strain is log(l/l0) for all components
exx,eyy,ezz=t.strain
err=.5*(exx+eyy)
# volumetric strain is not defined directly, and it is not needed either
eVol=float('nan')
else:
# triaxial phase:
# only axial strain (ezz) and volumetric strain (eVol) are known
#
# set the initial volume, if not yet done
if not hasattr(S.lab,'netVol0'): S.lab.netVol0=S.cell.volume*S.lab.relVol
# axial strain is known; xy components irrelevant (inactive)
ezz=t.strain[2]
# current volume / initial volume
eVol=math.log(vol/S.lab.netVol0)
# radial strain
err=.5*(eVol-ezz)
# undefined
exx=eyy=float('nan')
# deviatoric strain
eDev=ezz-(1/3.)*(2*err+ezz) # FIXME: is this correct?!
S.plot.addData(unbalanced=woo.utils.unbalancedForce(),i=S.step,
sxx=sxx,syy=syy,srr=.5*(sxx+syy),szz=szz,
exx=exx,eyy=eyy,err=err,ezz=ezz,
dotE=dotE,dotErr=.5*(dotE[0]+dotE[1]),
dotEMax=dotEMax,
dotEMax_z_neg=-dotEMax[2],
eDev=eDev,eVol=eVol,
vol=vol,
p=p,q=q,qDivP=qDivP,
isTriax=(1 if S.lab.stage=='triax' else 0), # to be able to filter data
# parTanPhi=S.lab.partMat.tanPhi,
#memTanPhi=S.lab.memMat.tanPhi,
#suppTanPhi=S.lab.suppMat.tanPhi
# save all available energy data
#Etot=O.energy.total()#,**O.energy
)
if not S.plot.plots:
S.plot.plots={
'i':('unbalanced',None,'vol'),'i ':(('sxx','syy','szz') if S.pre.planeStrain else ('srr','szz')),' i':(('exx','eyy','ezz','eVol') if S.pre.planeStrain else ('err','ezz','eVol')),'i ':('dotE_z','dotEMax_z'),
'eDev':(('qDivP','g-'),None,('eVol','r-')),'p':('q',),
# energy plot
#' i ':(O.energy.keys,None,'Etot'),
}
S.plot.xylabels={'i ':('step','Stress [Pa]',),' i':('step','Strains [-]','Strains [-]')}
S.plot.labels={
'sxx':r'$\sigma_{xx}$','syy':r'$\sigma_{yy}$','szz':r'$\sigma_{zz}$','srr':r'$\sigma_{rr}$','surfLoad':r'$\sigma_{\rm hydro}$',
'exx':r'$\varepsilon_{xx}$','eyy':r'$\varepsilon_{yy}$','ezz':r'$\varepsilon_{zz}$','err':r'$\varepsilon_{rr}$','eVol':r'$\varepsilon_{v}$','vol':'volume','eDev':r'$\varepsilon_d$','qDivP':'$q/p$','p':'$p$','q':'$q$',
'dotE_x':r'$\dot\varepsilon_{xx}$','dotE_y':r'$\dot\varepsilon_{yy}$','dotE_z':r'$\dot\varepsilon_{zz}$','dotE_rr':r'$\dot\varepsilon_{rr}$','dotEMax_z':r'$\dot\varepsilon_{zz}^{\rm max}$','dotEMax_z_neg':r'$-\dot\varepsilon_{zz}^{\rm max}$'
}
## adjust rate in the triaxial stage
if S.lab.stage=='triax':
t.maxStrainRate[2]=min(t.maxStrainRate[2]+S.pre.rateStep*S.pre.maxRates[1],S.pre.maxRates[1])
[docs]def compactionDone(S):
# if S.lab.compactMemoize:
print('Compaction done at step',S.step)
import woo
t=S.lab.triax
# set the current cell configuration to be the reference one
S.cell.trsf=Matrix3.Identity
S.cell.refHSize=S.cell.hSize
S.cell.nextGradV=Matrix3.Zero # avoid spurious strain from the last gradV value
## S.lab.leapfrog.damping=.7 # increase damping to a larger value
if not S.pre.planeStrain:
t.stressMask=0b0011 # z is strain-controlled, xy are stress-controlled
t.goal=(S.pre.sigIso,S.pre.sigIso,S.pre.stopStrain)
else:
t.stressMask=0b0001 # z is strain-controlled, y has zero strain, x is stress-controlled
t.goal=(S.pre.sigIso,0,S.pre.stopStrain)
# start with small rate, increases later
t.maxStrainRate=(S.pre.maxRates[2],S.pre.maxRates[2],S.pre.rateStep*S.pre.maxRates[1])
t.maxUnbalanced=10 # don't care about that
t.doneHook='import woo.pre.triax; woo.pre.triax.triaxDone(S)'
# recover friction angle
S.lab.partMat.tanPhi=S.pre.model.mats[0].tanPhi
# force update of contact parameters
S.lab.contactLoop.updatePhys='once'
try:
import woo.gl
woo.gl.Gl1_DemField.updateRefPos=True
except ImportError: pass
S.lab.stage='triax'
if S.pre.saveFmt:
out=S.pre.saveFmt.format(stage='compact',S=S,**(dict(S.tags)))
print('Saving to',out)
S.save(out)
[docs]def triaxDone(S):
print('Triaxial done at step',S.step)
if S.pre.saveFmt:
out=S.pre.saveFmt.format(stage='done',S=S,**(dict(S.tags)))
print('Saving to',out)
S.save(out)
S.stop()
import woo.utils
(repName,figs)=woo.utils.htmlReport(S,S.pre.reportFmt,'Triaxial test with rigid boundary',afterHead='',figures=[(None,f) for f in S.plot.plot(noShow=True,subPlots=False)],svgEmbed=True,show=True)
woo.batch.writeResults(S,defaultDb='triax.hdf5',series=S.plot.data,postHooks=[woo.pre.cylTriax.plotBatchResults],simulationName='triax',report='file://'+repName)