yade-users team mailing list archive
-
yade-users team
-
Mailing list archive
-
Message #18729
[Question #677795]: Triaxial test of flexible membrane particles flying out of the boundary
New question #677795 on Yade:
https://answers.launchpad.net/yade/+question/677795
Hello, everyone.
I am doing a triaxial test on a flexible film. I need an axial strain of about 10%. When I started running the program, the particles flew out of the cylindrical flexible membrane boundary. If anyone knows please tell me how to solve this problem. Thank you!
from yade import pack, plot
import os
# default parameters or from table
readParamsFromTable(noTableOk=True,
# type of test ['cyl','cube']
testType = 'cyl',
# material parameters
young = 30e6,
poisson = .2,
frictionAngle = 0.27,
sigmaT = 1.5e3,
epsCrackOnset = 1e-4,
relDuctility = 10,
# prestress
preStress = -1.5e5,
# axial strain rate
strainRate = -1,
# assamlby parameters
rParticle = .075e-3, #
width = 2e-3,
height = 5e-3,
bcCoeff = 5,
# facets division
nw = 24,
nh = 15,
# output specifications
fileName = 'young = 30e6,frictionAngle = 1.5',
exportDir = '/tmp',
runGnuplot = False,
runInGui = True,
)
from yade.params.table import *
assert testType in ['cyl','cube']
# materials
concMat = O.materials.append(CpmMat(
young=young,frictionAngle=frictionAngle,poisson=poisson,sigmaT=sigmaT,
epsCrackOnset=epsCrackOnset,relDuctility=relDuctility
))
frictMat = O.materials.append(FrictMat(
young=young,poisson=poisson,frictionAngle=frictionAngle
))
# spheres
pred = pack.inCylinder((0,0,0),(0,0,height),.5*width) if testType=='cyl' else pack.inAlignedBox((-.5*width,-.5*width,0),(.5*width,.5*width,height)) if testType=='cube' else None
sp=SpherePack()
sp = pack.randomDensePack(pred,spheresInCell=2000,radius=rParticle,memoizeDb='/tmp/triaxTestOnCylinder.sqlite',returnSpherePack=True)
spheres=sp.toSimulation(color=(0,1,1),material=concMat)
# bottom and top of specimen. Will have prescribed velocity
bot = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]<rParticle*bcCoeff]
top = [O.bodies[s] for s in spheres if O.bodies[s].state.pos[2]>height-rParticle*bcCoeff]
vel = strainRate*(height-rParticle*2*bcCoeff)
for s in bot:
s.shape.color = (1,0,0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0,0,-vel)
for s in top:
s.shape.color = Vector3(0,1,0)
s.state.blockedDOFs = 'xyzXYZ'
s.state.vel = (0,0,vel)
print 'Number of elements: ', len(O.bodies)
print 'Timestep',O.dt
print 'young = 30e6,frictionAngle = 0.27,preStress = -1.5e5,sigmaT = 1.5e2'
# facets
facets = []
if testType == 'cyl':
rCyl2 = .5*width / cos(pi/float(nw))
for r in xrange(nw):
for h in xrange(nh):
v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+0)/float(nh) )
v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+0)/float(nh) )
v3 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)), rCyl2*sin(2*pi*(r+1)/float(nw)), height*(h+1)/float(nh) )
v4 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)), rCyl2*sin(2*pi*(r+0)/float(nw)), height*(h+1)/float(nh) )
f1 = facet((v1,v2,v3),color=(0,0,1),material=frictMat)
f2 = facet((v1,v3,v4),color=(0,0,1),material=frictMat)
facets.extend((f1,f2))
elif testType == 'cube':
nw2 = nw/4
for r in xrange(nw2):
for h in xrange(nh):
v11 = Vector3( -.1*width + (r+0)*width/nw2, -.1*width, height*(h+0)/float(nh) )
v12 = Vector3( -.1*width + (r+1)*width/nw2, -.1*width, height*(h+0)/float(nh) )
v13 = Vector3( -.1*width + (r+1)*width/nw2, -.1*width, height*(h+1)/float(nh) )
v14 = Vector3( -.1*width + (r+0)*width/nw2, -.1*width, height*(h+1)/float(nh) )
f11 = facet((v11,v12,v13),color=(0,0,1),material=frictMat)
f12 = facet((v11,v13,v14),color=(0,0,1),material=frictMat)
v21 = Vector3( +.5*width, -.5*width + (r+0)*width/nw2, height*(h+0)/float(nh) )
v22 = Vector3( +.5*width, -.5*width + (r+1)*width/nw2, height*(h+0)/float(nh) )
v23 = Vector3( +.5*width, -.5*width + (r+1)*width/nw2, height*(h+1)/float(nh) )
v24 = Vector3( +.5*width, -.5*width + (r+0)*width/nw2, height*(h+1)/float(nh) )
f21 = facet((v21,v22,v23),color=(0,0,1),material=frictMat)
f22 = facet((v21,v23,v24),color=(0,0,1),material=frictMat)
v31 = Vector3( +.5*width - (r+0)*width/nw2, +.5*width, height*(h+0)/float(nh) )
v32 = Vector3( +.5*width - (r+1)*width/nw2, +.5*width, height*(h+0)/float(nh) )
v33 = Vector3( +.5*width - (r+1)*width/nw2, +.5*width, height*(h+1)/float(nh) )
v34 = Vector3( +.5*width - (r+0)*width/nw2, +.5*width, height*(h+1)/float(nh) )
f31 = facet((v31,v32,v33),color=(0,0,1),material=frictMat)
f32 = facet((v31,v33,v34),color=(0,0,1),material=frictMat)
v41 = Vector3( -.5*width, +.5*width - (r+0)*width/nw2, height*(h+0)/float(nh) )
v42 = Vector3( -.5*width, +.5*width - (r+1)*width/nw2, height*(h+0)/float(nh) )
v43 = Vector3( -.5*width, +.5*width - (r+1)*width/nw2, height*(h+1)/float(nh) )
v44 = Vector3( -.5*width, +.5*width - (r+0)*width/nw2, height*(h+1)/float(nh) )
f41 = facet((v41,v42,v43),color=(0,0,1),material=frictMat)
f42 = facet((v41,v43,v44),color=(0,0,1),material=frictMat)
facets.extend((f11,f12,f21,f22,f31,f32,f41,f42))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
f.state.mass = mass
f.state.blockedDOFs = 'XYZz'
#############################################
plot.plots={'eps':('sigma')} #,'sigma.50')},'t':('eps')} #'sigma.25','sigma.50','sigma.75')}
O.saveTmp('initial');
O.timingEnabled=False
#############################################
# plots
plot.plots = { 'e':('s',), }
def plotAddData():
f1 = sum(O.forces.f(b.id)[2] for b in top)
f2 = sum(O.forces.f(b.id)[2] for b in bot)
f = .5*(f2-f1)
s = f/(pi*.25*width*width) if testType=='cyl' else f/(width*width) if testType=='cube' else None
e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
plot.addData(
i = O.iter,
s = s,
e = e,
f1 = f1,
f2 = f2
)
# apply prestress to facets
def addForces():
for f in facets:
n = f.shape.normal
a = f.shape.area
O.forces.addF(f.id,preStress*a*n)
# stop condition and exit of the simulation
def stopIfDamaged(maxEps=10e-2):
extremum = max(abs(s) for s in plot.data['s'])
s = abs(plot.data['s'][-1])
e = abs(plot.data['e'][-1])
if O.iter < 1000 or s > .5*extremum and e < maxEps:
return
f = os.path.join(exportDir,fileName)
print 'gnuplot',plot.saveGnuplot(f,term='png')
if runGnuplot:
import subprocess
os.chdir(exportDir)
subprocess.Popen(['gnuplot',f+'.gnuplot']).wait()
print 'Simulation finished'
O.pause()
#sys.exit(0) # results in some threading exception
O.dt=.5*utils.PWaveTimeStep()
enlargeFactor=1.5
O.engines=[
ForceResetter(),
InsertionSortCollider([
Bo1_Sphere_Aabb(aabbEnlargeFactor=enlargeFactor,label='bo1s'),
Bo1_Facet_Aabb()
]),
InteractionLoop(
[
Ig2_Sphere_Sphere_ScGeom(interactionDetectionFactor=enlargeFactor,label='ss2d3dg'),
Ig2_Facet_Sphere_ScGeom(),
],
[
Ip2_CpmMat_CpmMat_CpmPhys(cohesiveThresholdIter=O.iter+5),
Ip2_FrictMat_CpmMat_FrictPhys(),
Ip2_FrictMat_FrictMat_FrictPhys(),
],
[
Law2_ScGeom_CpmPhys_Cpm(),
Law2_ScGeom_FrictPhys_CundallStrack(),
],
),
PyRunner(iterPeriod=1,command="addForces()"),
NewtonIntegrator(damping=.3),
CpmStateUpdater(iterPeriod=50,label='cpmStateUpdater'),
PyRunner(command='plotAddData()',iterPeriod=10),
PyRunner(iterPeriod=50,command='stopIfDamaged()'),
]
# run one step
O.step()
# reset interaction detection enlargement
bo1s.aabbEnlargeFactor=ss2d3dg.interactionDetectionFactor=1.0
# initialize auto-updated plot
if runInGui:
plot.plot()
try:
from yade import qt
renderer=qt.Renderer()
# uncomment following line to exagerate displacement
#renderer.dispScale=(100,100,100)
except:
pass
--
You received this question notification because your team yade-users is
an answer contact for Yade.