← Back to team overview

yade-users team mailing list archive

[Question #679405]: 2D simulation with sphere using PeriTriaxController

 

New question #679405 on Yade:
https://answers.launchpad.net/yade/+question/679405

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

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.