yade-users team mailing list archive

[Rahul R <question701028@xxxxxxxxxxxxxxxxxxxxx>]

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

Hi, I'm new in Yade. I wanted to do a triaxial test on a cylindrical specimen. I found a code on github and modified the code for my specimen (the original code creates a spherical pack and do triaxial test but I created a cylindrical sample  of desired porosity and imported the coordinates to this code). I ran the code but the cylindrical membrane that was supposed to be created around the spherical pack is created away from the pack. Can somebody help me?

Here is the link to the actual code: http://bazaar.launchpad.net/~yade-pkg/yade/git-trunk/view/head:/examples/concrete/triax.py
Here is the code:

# -*- encoding=utf-8 -*-
from __future__ import print_function
################################################################################
#
# Triaxial test. Axial strain rate is prescribed and transverse prestress.
# Test is possible on prism or cylinder
# An independent c++ engine may be created from this script in the future.
#
################################################################################
from builtins import range
from yade import ymport, plot
from yade import pack, plot
import os

# default parameters or from table
readParamsFromTable(noTableOk=True,
	# material parameters
	young = 10e7,
	poisson = .3,
	frictionAngle = 0.15,
	sigmaT = 1.5e6,
	
	# prestress
	preStress = -3e6,
	# axial strain rate
	strainRate = -100,

	# assamlby parameters
	rParticle = 0.0005,
	width = 5,
	height = 9.85,
	bcCoeff = 5,

	# facets division
	nw = 48,
	nh = 30,

	# output specifications
	fileName = 'test',
	exportDir = '/tmp',
	runGnuplot = False,
	runInGui = True,
)
from yade.params.table import *

# materials
sphereMat = O.materials.append(CohFrictMat(young=young,poisson=poisson,frictionAngle=frictionAngle,alphaKr=0.25,alphaKtw=0,etaRoll=0.005,etaTwist=0,normalCohesion=5e6,shearCohesion=5e6,momentRotationLaw=True,density=2530))

frictMat = O.materials.append(FrictMat(
	young=young,poisson=poisson,frictionAngle=frictionAngle
))

# spheres
sp= ymport.textExt('DensepackRD70-80',format='x_y_z_r',shift= Vector3(0,0,0.5*height), scale=1.0, material=sphereMat)
spheres = O.bodies.append(sp)

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

# 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)
top_limit = 0
top_id = 0
for s in top:
	if s.state.pos[2]>=top_limit:
		top_limit = s.state.pos[2]
		top_id = s.id
bot_limit = height
bot_id = 0
for s in bot:
	if s.state.pos[2]<=bot_limit:
		bot_limit = s.state.pos[2]
		bot_id = s.id
# facets
facets = []
rCyl2 = .5*width / cos(pi/float(nw))
for r in range(nw):
	for h in range(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))

O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
	f.state.mass = mass
	f.state.blockedDOFs = 'XYZz'

# 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)
	e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) / (height-rParticle*2*bcCoeff)
	plot.addData(
		i = O.iter,
		s = s,
		e = e,
	)

# 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=5e-3):
	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=100),
	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

# run
O.run()


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