PeriTriaxController of different materials

Asked by jamespaul on 2019-02-05

Hi,

I do isotropic compression at first, and then I do one-dimensional compression.I get stress and strain in real time.

I found that if I used the default material, everything would be fine. In a one-dimensional compression, the pressure slowly increases.But when I use my own definition of the material, in one-dimensional compression, the pressure goes down to zero and then goes up.I have not made any changes to the code except for the addition of material=spheres when generating the particles.

Why is it that just changing the material creates something that can't happen in reality? Do I need to adjust the maxUnbalanced and relStressTol for isotropic compression?

Please tell me exactly where I made a mistake.Run my code if you can.

Thanks,

James

#####################################
# 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.0005))
s=O.materials.append(FrictMat(young=64e9,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)
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=.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,.0001)
   triax.doneHook='triaxFinished()'
   triax.maxUnbalanced=0.000001

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

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2019-02-06
Last reply:
2019-02-07

This question was reopened

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.

jamespaul (jamespauljames) said : #2

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

jamespaul (jamespauljames) said : #4

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

jamespaul (jamespauljames) said : #7

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?

jamespaul (jamespauljames) said : #9

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?

jamespaul (jamespauljames) said : #11

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? :)

jamespaul (jamespauljames) said : #13

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()

Can you help with this problem?

Provide an answer of your own, or ask jamespaul for more information if necessary.

To post a message you must log in.