2D clump in the compression test which the mean stress is constant

Hello,

please read [1].

Thanks for sharing the code. However, it is of little information for us, because:
1) it loads files we do not have access to (Packing2D_area.txt)
2) (more importantly) it loads previously saved simulation, with all its bodies, engines... anything not present here. So we have no idea what's going on - how the sample is loade (using periodicity? walls? ...), how many clumps, what shapes, what laws......

So please try to make your script minimal - deleting unnecessary commented lines (improves readability), using "just a few" clumps (then no need of external files, runs quickly) etc. - and working (we can copy-paste it and test it or at least see all information).
By creating a true MWE, you often find the problem yourself.

thanks
Jan

[1] https://yade-dem.org/wiki/Howtoask

Hi，
The first script is clump consolidation；

from yade import pack,plot,qt,export,utils
O.periodic=True
O.materials.append(FrictMat(young=6.e8,poisson=.8,frictionAngle=radians(26),density=2650))

sp = pack.SpherePack()
size = .01
sp.makeCloud(minCorner=(0,0,.0020),maxCorner=(size,size,.0020),rMean=.0002,rRelFuzz=.4,num=400,periodic=True,seed=1)
sp.toSimulation()
O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,.004)
print len(O.bodies)
for p in O.bodies:
   p.state.blockedDOFs = 'zXY'
   p.state.mass = 2650 * 0.004 * pi * p.shape.radius**2 # 0.1 = thickness of cylindrical particle
   inertia = 0.5 * p.state.mass * p.shape.radius**2
   p.state.inertia = (.5*inertia,.5*inertia,inertia)

O.engines = [
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),
   GlobalStiffnessTimeStepper(),
   NewtonIntegrator(damping=0.7),
   PeriTriaxController(
      dynCell=True,
      goal=(-1.e2,-1.e2,0),
      stressMask=3,
      relStressTol=.001,
      maxUnbalanced=.001,
      maxStrainRate=(.3,.3,.0),
      doneHook='Finished()',
      label='biax'
   ),
     PyRunner(command='plotAddData()',iterPeriod=100),
]

def plotAddData():
 plot.addData(
  iter=O.iter,iter_=O.iter,
  sxx=biax.stress[0],syy=biax.stress[1],szz=biax.stress[2],
  exx=O.cell.size[0],eyy=O.cell.size[1],ezz=O.cell.size[2],
  Z=avgNumInteractions(),
  Zm=avgNumInteractions(skipFree=True),
                voidratio=utils.voidratio2D(0.004),
  unbalanced=utils.unbalancedForce(),
  t=O.time,
  gWork=O.energy['gravWork'],
  Ep=O.energy['elastPotential'],
  Edamp=O.energy['nonviscDamp'],
  Ediss=O.energy['plastDissip'],
  Ekin=utils.kineticEnergy(),
        Etot=O.energy.total(),**O.energy

 )
        plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio'))
        plot.saveDataTxt('energyFile',vars=('t','Etot','unbalanced','gWork','Edamp','Ekin'))

O.trackEnergy=True

def Finished():
 O.save('isotropicState.xml')
 print 'Finished'
 print utils.voidratio2D(0.004)
 O.pause()

yade.qt.Controller(), yade.qt.View()
O.run()
plot.plot(subPlots=True)

The second script is the compression test:

from yade import pack,plot,qt,export
from math import fabs
import numpy as np
O.load('isotropicState.xml')
O.trackEnergy=True
setContactFriction(radians(26.6))

for p in O.bodies:
   p.state.blockedDOFs='zXY'

def servo():
 tol=0.01
 b=1
 sr_x=-.1
 sr_y=0.
 sr_z=0.
        epsilonmax=10
 p=-100
 sx_curr=utils.getStress()[0,0]
 sy_curr=utils.getStress()[1,1]
        sy_req=2*p-sx_curr
        gy=2*epsilonmax/sy_req
 sry_curr=O.cell.velGrad[1,1]
 srz_curr=O.cell.velGrad[2,2]
 if fabs(sy_curr-sy_req)>=tol:
                sr_y=-gy*(sy_req-sy_curr)
 else:
  sr_y=sry_curr
 O.cell.velGrad=Matrix3(sr_x,0,0,0,sr_y,0,0,0,sr_z)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
     [Ip2_FrictMat_FrictMat_FrictPhys()],
     [Law2_ScGeom_FrictPhys_CundallStrack()]),
 PyRunner(command='fabric()',iterPeriod=10000),
 NewtonIntegrator(damping=0.2),
        GlobalStiffnessTimeStepper(),
 PyRunner(command='servo()', iterPeriod=1),
 PyRunner(command='plotAddData()',iterPeriod=100),
]

plot.live=True
plot.plots={'iter':('sxx','syy','szz'),'iter_':('exx','eyy','ezz'), ' iter':('voidratio')
}
def plotAddData():
 plot.addData(
  iter=O.iter,iter_=O.iter,
  sxx=utils.getStress()[0,0],
  syy=utils.getStress()[1,1],
  szz=utils.getStress()[2,2],
  exx=O.cell.size[0],
  eyy=O.cell.size[1],
  ezz=O.cell.size[2],
  Z=avgNumInteractions(),
  Zm=avgNumInteractions(skipFree=True),
                voidratio=utils.voidratio2D(0.004),
  unbalanced=utils.unbalancedForce(),
  t=O.time,
  gWork=O.energy['gravWork'],
  Edamp=O.energy['nonviscDamp'],
  Ekin=utils.kineticEnergy(),
        Etot=O.energy.total(),**O.energy

 )
 plot.saveDataTxt('macroFile',vars=('t','exx','eyy','ezz','sxx','syy','szz','Z','Zm','voidratio'))
 plot.saveDataTxt('energyFile',vars=('t','Etot','unbalanced','gWork','Edamp','Ekin'))

O.run()

Hello,

thanks for improvements, but please really try to make the script minimal, not only from code point of view, but the situation itself..
I stopped the initial compaction after 1 mio iterations running 20 minutes (and still the target unbalanced force order of magnitude far away)..

To illustrate and reproduce the problem, I think creating a few clumps "by hand" should be enough..
Or (for example) the problem might be in the servo() function and would be there running just spheres.. By simplifying the problem, often you localize the problem much more precisely and then also get the answer, solution or advices better and faster.
So please make the example (code, model itself) minimal. Help us to help you :-)

thanks
Jan

This question was expired because it remained in the 'Needs information' state without activity for the last 15 days.