PeriTriaxController of different materials

Hi,
If dynCell=True then the cell deformation is dynamic. I mean, inertial, like a mass-spring system.
Depending on stiffness of particles, pseudo-mass of the cell, and damping, it will oscillate more or less.
Related parameters in PeriTriaxController: growDamping, mass, maxStrainRate.

Thanks Bruno,

I don't quite understand what you mean by oscillations. I don't think the pressure is going to go down no matter what the shock is. Should I change it to dynCell=False？

James

Thanks Bruno,

I don't quite understand what you mean by oscillations. I don't think the pressure is going to go down no matter what the shock is. Should I change it to dynCell=False？

James

> I don't think the pressure is going to go down

Well... no need to think here. Check output and you will see (converging) oscillations around the goal value.
Please fix errors if you need more insight.

####### ERROR1
Yade [1]: plot.plot()
---------------------------------------------------------------------------
IndexError Traceback (most recent call last)
/usr/bin/yadedaily in <module>()
----> 1 plot.plot()

/usr/lib/x86_64-linux-gnu/yadedaily/py/yade/plot.py in plot(noShow, subPlots)
    593 global currLineRefs
    594 figs=set([l.line.axes.get_figure() for l in currLineRefs])
--> 595 if not hasattr(list(figs)[0],'show') and not noShow:
    596 import warnings
    597 warnings.warn('plot.plot not showing figure (matplotlib using headless backend?)')

IndexError: list index out of range

####### ERROR2

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
/usr/bin/yadedaily in <module>()

/usr/bin/yadedaily in compactionFinished()
     43 print('aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa')
     44 O.cell.trsf=Matrix3.Identity
---> 45 for i in n:
     46 O.bodies[i].material = O.materials[s]
     47 triax.goal=(sigmaIso,sigmaIso,-1000000)

NameError: global name 'n' is not defined

Also note that saving whole plot data to text file each time you add one point in addPlotData is a sure way to waste cpu time and kill your HD in the long run.
Bruno

Thanks Bruno,

I made a small mistake and corrected my code.

# encoding: utf-8

sigmaIso=-25000

from yade import pack, qt, plot

O.periodic=True

spheres=O.materials.append(FrictMat(young=64e8,poisson=0.12,density=2650,frictionAngle=0.0005))
s=O.materials.append(FrictMat(young=64e8,poisson=0.12,density=2650,frictionAngle=0.24))

'''sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.1,rRelFuzz=.3,periodic=True)
n=sp.toSimulation(material=spheres)'''

psdSizes=[0.15,.18,.2,.21,.25,.27,.30,.31,.34,.38]
psdCumm=[0.025,0.11,0.27,0.45,0.73,0.88,0.95,0.98,0.99,1]

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1.5,1.5,1.5),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,periodic=True)
n=sp.toSimulation(material=spheres)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
 PeriTriaxController(label='triax',
  goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
  dynCell=True,maxStrainRate=(10,10,10),
  maxUnbalanced=.1,relStressTol=1e-3,
  doneHook='compactionFinished()'
    ),
 NewtonIntegrator(damping=.1),
 PyRunner(command='addPlotData()',iterPeriod=1000),
]
O.dt=.5*PWaveTimeStep()

def addPlotData():
 plot.addData(sxx=triax.stress[0],syy=triax.stress[1],szz=triax.stress[2],
        exx=triax.strain[0],eyy=triax.strain[1],ezz=triax.strain[2],
 q=triax.stress[2]-0.5*(triax.stress[0]+triax.stress[1]),
 p=(triax.stress[0]+triax.stress[1]+triax.stress[2])/3,
  )
 print('szz=',triax.stress[2],'ezz=',triax.strain[2])
 #plot.saveDataTxt('qp1.txt',vars=('q','p','szz','ezz'))

def compactionFinished():
 print('aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa')
 O.cell.trsf=Matrix3.Identity
 for i in n:
         O.bodies[i].material = O.materials[s]
 triax.goal=(sigmaIso,sigmaIso,-1000000)
 triax.stressMask=7
 triax.maxStrainRate=(0,0,.000001)
 triax.doneHook='triaxFinished()'
 triax.maxUnbalanced=0.000001

def triaxFinished():
 print 'Finished'
 O.pause()

There is still no plot despite the addPlotData.
Anyway, don't you see the oscillations in the printed numbers?

Bruno,

I saw the oscillations.But the oscillation is purely in the first isotropic compression and does not cause a reduction in pressure in the second compression.

#############################
# encoding: utf-8

sigmaIso=-25000

from yade import pack, qt, plot

O.periodic=True

spheres=O.materials.append(FrictMat(young=64e9,poisson=0.12,density=2650,frictionAngle=0.231))

length=3

psdSizes=[0.15,.18,.2,.21,.25,.27,.30,.31,.34,.38]
psdCumm=[0.025,0.11,0.27,0.45,0.73,0.88,0.95,0.98,0.99,1]

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(length,length,length),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,periodic=True)
n=sp.toSimulation(material=spheres)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
        [Ig2_Sphere_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_FrictPhys()],
        [Law2_ScGeom_FrictPhys_CundallStrack()]
    ),
 PeriTriaxController(label='triax',
  goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
  dynCell=True,maxStrainRate=(.1,.1,.1),
  maxUnbalanced=.01,relStressTol=1e-3,
  doneHook='compactionFinished()'
    ),
 NewtonIntegrator(damping=.1),
 PyRunner(command='addPlotData()',iterPeriod=10000),
]
O.dt=.5*PWaveTimeStep()

def addPlotData():
 plot.addData(sxx=triax.stress[0],syy=triax.stress[1],szz=triax.stress[2],
         exx=triax.strain[0],eyy=triax.strain[1],ezz=triax.strain[2],
  q=triax.stress[2]-0.5*(triax.stress[0]+triax.stress[1]),
  p=(triax.stress[0]+triax.stress[1]+triax.stress[2])/3,
  unbalanced=unbalancedForce(),i=O.iter,
  #k=(triax.stress[0]*0.5+triax.stress[1]*0.5)/triax.stress[2]
  )
 print('szz=',triax.stress[2],'ezz=',triax.strain[2])
 #plot.saveDataTxt('qp1.txt',vars=('q','p','szz','ezz'))

plot.plots={'i':('unbalanced',),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz'),}
plot.plot()

def compactionFinished():
 print('aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa')
 O.cell.trsf=Matrix3.Identity
 s=O.materials.append(FrictMat(young=64e9,poisson=0.12,density=2650,frictionAngle=0.24))
 for i in n:
         O.bodies[i].material = O.materials[s]
 triax.goal=(sigmaIso,sigmaIso,-0.005)
 dynCell=False,
 triax.stressMask=3
 triax.maxStrainRate=(0,0,.001)
 triax.doneHook='triaxFinished()'
 triax.maxUnbalanced=0.000001

def triaxFinished():
 print 'Finished'
 O.pause()

I don't understand your question, sorry.
Do you mean that I should try this 3rd version of the script to discover something new?

Yeah.

In the second operation, a strange phenomenon occurred, that is, as the one-dimensional compression, the pressure decreased.

Running the program can take 10 minutes or more. You can adjust some non-critical parameters to speed up the simulation time.

The problem has been bothering me for a long time. I wish you could help me.

Thank you

> Running the program can take 10 minutes or more. You can adjust some non-critical parameters

Would you mind adjusting them for me? :)

Hi Bruno,

Please run my program. It takes about two minutes.Thank you so much!

1.I have defined a material called 'sphere'. First, it does isotropic compression. The measured data are very correct.

2.Then,using 'doneHook', it does one-dimensional compression. So before compression, I applied the 'changeMaterial' function to change the friction.Am I correct in changing the material?

3.The strange thing is that in one-dimensional compression, the pressure becomes zero.If I use the default material, everything is fine, that is, the pressure increases slowly with one-dimensional compression.I used my own material. Why is this anti-science?

James

##########################
# encoding: utf-8
from yade import pack, qt, plot

sigmaIso=-25000

O.periodic=True

sphere=O.materials.append(FrictMat(young=64e9,poisson=0.12,density=2650,frictionAngle=0.231,label='sphere'))
#s=O.materials.append(FrictMat(young=64e9,poisson=0.12,density=2650,frictionAngle=0.24))

multiple=.01
length=2

psdSizes=[0.15*multiple,.18*multiple,.2*multiple,.21*multiple,.25*multiple,.27*multiple,.30*multiple,.31*multiple,.34*multiple,.38*multiple]
psdCumm=[0.025,0.11,0.27,0.45,0.73,0.88,0.95,0.98,0.99,1]

sp=pack.SpherePack()
sp.makeCloud((0,0,0),(length*multiple,length*multiple,length*multiple),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,periodic=True)
n=sp.toSimulation()

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 PeriTriaxController(label='triax',
  goal=(sigmaIso,sigmaIso,sigmaIso),stressMask=7,
  dynCell=True,maxStrainRate=(5,5,5),
  maxUnbalanced=.1,relStressTol=1e-3,
  doneHook='compactionFinished()'
 ),
 NewtonIntegrator(damping=.1),
 PyRunner(command='addPlotData1()',iterPeriod=10000),
]
O.dt=.5*PWaveTimeStep()
qt.Controller()

def addPlotData1():
 plot.addData(sxx=-triax.stress[0],syy=-triax.stress[1],szz=-triax.stress[2],
         exx=-triax.strain[0],eyy=-triax.strain[1],ezz=-triax.strain[2],
  unbalanced=unbalancedForce(),i=O.iter,
  )
 print('szz=',-triax.stress[2],'ezz=',-triax.strain[2],'Force=',unbalancedForce())

plot.plots={'i':('unbalanced',),'i ':('sxx','syy','szz'),' i':('exx','eyy','ezz')}
plot.plot()

def changeMaterial():
    s=O.materials.append(FrictMat(young=64e9,poisson=0.12,density=2650,frictionAngle=0.24))
 for i in n:
         O.bodies[i].material.frictionAngle=0.24
         #O.bodies[i].material = O.materials[s]

def compactionFinished():
    print('Now, one-dimensional compression!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!')
 O.engines += [PyRunner(command='changeMaterial()',initRun=1,nDo=1)];
    O.cell.trsf=Matrix3.Identity
   triax.goal=(-1000000,-1000000,-1000000)
 triax.dynCell=True
 triax.relStressTol=1e-3
 triax.stressMask=7
 triax.maxStrainRate=(0,0,0.00001)
 triax.doneHook='triaxFinished()'
 triax.maxUnbalanced=0.1

def triaxFinished():
    print 'Finished'
    O.pause()