2D simulation with sphere using PeriTriaxController

Asked by nie jiayan

Hello, everyone. Now I do a 2D sphere biaxial compression simulation in Yade using PeriTriaxController. But one problem exists, though I use a low loading compression rate in y direction, the confine stress in x direction can not confined to a certain pressure. In the simulation, p.state.blockedDOFs = 'zXY' is setted. I do not figure out the reason, could someone give some suggestions? The following is my code. Thanks in advance.

from yade import pack, qt, plot
import time

############################################
### DEFINING VARIABLES AND MATERIALS ###
############################################

ts = time.time()
pressure = -1e5
size=0.24

## create materials for spheres for initial consolidation
O.materials.append(FrictMat(young=6.0e8,poisson=0.8,frictionAngle=.0))

# setup periodic boundary, insert the packing

sp = pack.SpherePack()
sp.makeCloud(minCorner=(0,0,.05),maxCorner=(size,size,.05),rMean=.005,rRelFuzz=.4,num=400,periodic=True,seed=1)
sp.toSimulation()
O.cell.hSize = Matrix3(size,0,0, 0,size,0, 0,0,.1)

O.dt=.5*PWaveTimeStep()

for p in O.bodies:
        p.state.blockedDOFs = 'zXY'
        p.state.mass = 2650 * 0.1 * pi * p.shape.radius**2 # 0.1 = thickness of cylindrical particle
        inertia = 0.5 * p.state.mass * p.shape.radius**2 # in 2D simulation, the inertia is important, should be calculate carefully.
        p.state.inertia = (.5*inertia,.5*inertia,inertia)

print(len(O.bodies))
############################
### DEFINING ENGINES ###
############################

triax=PeriTriaxController(
  # specify target values and whether they are strains or stresses
  goal=(pressure,pressure,0),stressMask=3,
  # type of servo-control
  dynCell=True,maxStrainRate=(0.5,0.5,0),
  # wait until the unbalanced force goes below this value
  maxUnbalanced=0.001,relStressTol=1e-3,
                # call this function when goal is reached and the packing is stable
  doneHook='consolidationFinished()'
 )

newton=NewtonIntegrator(damping=.1)

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb()]),
 InteractionLoop(
  [Ig2_Sphere_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 triax,
 newton
]

def addPlotData():

   pp = utils.porosity() #this is the porosity of the cell.
   ee = pp / (1-pp) #this is the void ratio of the 3D cell.
   evv=-triax.strain[0]-triax.strain[1]
   ezv=-triax.strain[2]
   szv=-triax.stress[2]
   sigmadv=-triax.stress[1]+triax.stress[0]
   sigmam=triax.stress[1]+triax.stress[0]
   sigma0v=-triax.stress[0]
   sigma1v=-triax.stress[1]
   ssdv=-triax.strain[1]
   srdv=-2.0*sigmadv/sigmam

   plot.addData(unbalanced=unbalancedForce(),
    ev=evv,
                qp=srdv,
                ssd=ssdv,
                ez=ezv,
                sz=szv,
                sigma0=sigma0v,
                sigma1=sigma1v
   )

# define what to plot//since they have the same x (i), the latter i should have a space in it 'i '
plot.plots={'ssd':('unbalanced',),'ssd ':('qp',),' ssd':('ev',),'ssd ':('ez',),' ssd':('sz',),' ssd':('sigma0',)
   #each 'i' should have different number or location of space. Or the latter will overwrite the former one.
}
# show the plot
plot.plot()

yade.qt.Controller(), yade.qt.View()

# output the necessary data which are useful to postprocess
def outputData():
   fout = open(filename, 'a')
   pp = utils.porosity() #this is the porosity of the cell.
   ee = pp / (1-pp) #this is the void ratio of the 3D cell.
   evv=-triax.strain[0]-triax.strain[1]
   sigmadv=-triax.stress[1]+triax.stress[0]
   sigmam=triax.stress[1]+triax.stress[0]
   ssdv=-triax.strain[1]
   srdv=-2.0*sigmadv/sigmam
   fout.write(str(ssdv)+' '+str(srdv)+' '+str(evv)+' '+str(ee)+' '+str(unbalancedForce())+'\n')
   fout.close()

def consolidationFinished():
        triax.goal=(pressure,-.20,0)
        triax.stressMask=1
        triax.doneHook='triaxFinished()'
        O.engines += [PyRunner(command='addPlotData()',iterPeriod=100)]
        # set the current cell configuration to be the reference one
 setContactFriction(0.5)
        O.cell.trsf=Matrix3.Identity
        O.cell.velGrad=Matrix3.Zero
        for p in O.bodies:
            p.state.vel = Vector3.Zero
            p.state.angVel = Vector3.Zero
            p.state.refPos = p.state.pos
            p.state.refOri = p.state.ori
 print 'consolidationFinished'
        print('time: '+str((time.time()-ts)/60.)+'min')

def triaxFinished():
 print 'Finished'
        print('time: '+str((time.time()-ts)/60.)+'min')
        O.save('shear_dense_packing_100kPa.yade.gz')
 O.pause()

Question information

Language:
English Edit question
Status:
Expired
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Launchpad Janitor (janitor) said :
#1

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