yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #19289
[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.