Triaxial test cylindrical membrane created away from the pack

Hello,

> I'm new in Yade.

welcome :-)

> I found a code on github ...
> http://bazaar.launchpad.net ...

we are using GitLab currently as a main repository platform [1]

> the cylindrical membrane ... is created away from the pack.

> width = 5,
> height = 9.85
> ...
> sp= ymport.textExt(...)
> spheres = O.bodies.append(sp)

How did you create the packing? What is its dimensions and position?
I guess it simply does not fit together with defined cylinder dimensions and position..

Shift and/or resize the packing and/or the cylinder.

Cheers
Jan

[1] https://gitlab.com/yade-dev/trunk

Hi Jan,
Thank you for your reply.
I tried scaling the packing, but that doesn't seem to eliminate the issue.
As for the first question, I created the specimen by creating a cylindrical
pack and depositing it in a cylindrical facet box of the desired diameter
(5 units), and doing compression and extension cycles manually till the
desired porosity range and height of the specimen was reached. I've shared
the coordinates file in the mail. Here is the code I used to create the
sample:

from yade import pack

#From PSD 1 of PSD curves. Density used is 2530kg/m3
psdSizes,psdCumm=[0.0694, 0.07508, 0.08895, 0.11967, 0.16797, 0.20143,
0.26611, 0.30037, 0.34524, 0.39682, 0.414, 0.44251, 0.46447, 0.49646,
0.52426, 0.60258, 0.70105, 0.73584, 0.80093, 0.89857, 0.99597, 1.23104,
1.64618, 1.99808, 2.57658, 3.44547, 4.69182],[0.00039, 0.01361, 0.02507,
0.04094, 0.08854, 0.12644, 0.20049, 0.23928, 0.29833, 0.38472, 0.41469,
0.45788, 0.49931, 0.55573, 0.59892, 0.69853, 0.77434, 0.80078, 0.83604,
0.87924, 0.90833, 0.94358, 0.96915, 0.98325, 0.99207, 0.99647, 1.]

#Create a cylinder pack
sp=pack.SpherePack();
sp.makeCloud((25,25,55),(29,29,110),psdSizes=psdSizes,psdCumm=psdCumm,distributeMass=True,num=35000);
#This works with f=0.15.
cyl=pack.inCylinder((27,27,55),(27,27,110),radius=2) #Creating a predicate
to fill the pack.
sp1=pack.filterSpherePack(cyl,sp,True) #Create a cylindrical pack using the
predicate "cyl" and the psd sizes and %fines provided in "makeCloud".
sp1.toSimulation();

#Creating a cylindrical box.
O.bodies.append(geom.facetCylinder((27,27,81),radius= 2.5,height=
60,fixed=True));

#Creating and identifying interactions btw particles, walls and facets.
O.engines=[
ForceResetter(),
InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
InteractionLoop(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_CundallStrack()]
              ),
NewtonIntegrator(gravity=(0,0,-9.81),damping=0.5),
#PyRunner(command='compaction()',iterPeriod=50)
]
O.dt=.5*PWaveTimeStep()

On Tue, 22 Mar 2022 at 16:41, Jan Stránský <
<email address hidden>> wrote:

> Your question #701028 on Yade changed:
> https://answers.launchpad.net/yade/+question/701028
>
> Status: Open => Answered
>
> Jan Stránský proposed the following answer:
> Hello,
>
> > I'm new in Yade.
>
> welcome :-)
>
> > I found a code on github ...
> > http://bazaar.launchpad.net ...
>
> we are using GitLab currently as a main repository platform [1]
>
> > the cylindrical membrane ... is created away from the pack.
>
> > width = 5,
> > height = 9.85
> > ...
> > sp= ymport.textExt(...)
> > spheres = O.bodies.append(sp)
>
> How did you create the packing? What is its dimensions and position?
> I guess it simply does not fit together with defined cylinder dimensions
> and position..
>
> Shift and/or resize the packing and/or the cylinder.
>
> Cheers
> Jan
>
> [1] https://gitlab.com/yade-dev/trunk
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
> https://answers.launchpad.net/yade/+question/701028/+confirm?answer_id=0
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
> https://answers.launchpad.net/yade/+question/701028
>
> You received this question notification because you asked the question.
>

> I tried scaling the packing, but that doesn't seem to eliminate the issue.

if the issue is that the cylinder is not positioned correctly, then the solution is to position it correctly :-)

> I've shared the coordinates file in the mail.

launchpad questions does not work well with mail.. you can put it here directly,

> Here is the code I used to create the sample:

what are the stop conditions? what number of iterations to run it?

Cheers
Jan

Hi Jan,
Sorry for the delayed reply. Your solution of positioning the cylinder pack
correctly worked, but the membrane is not applying confining stress to the
pack. It shows some error in the terminal saying:

<FATAL ERROR> ThreadRunner:35 void yade::ThreadRunner::run(): Exception
occured:
PyRunner error.

COMMAND: 'stopIfDamaged()'

ERROR:
max() arg is an empty sequence

STACK TRACE:
Traceback (most recent call last):

  File "<string>", line 1, in <module>

  File "Triax.py", line 119, in stopIfDamaged
    extremum = max(abs(s) for s in plot.data['s'])

ValueError: max() arg is an empty sequence

Why does this happen? This is the updated 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(-27,-27,-51), scale=1.0, material=sphereMat)
spheres = O.bodies.append(sp)

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

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

> What are the stop conditions? what number of iteration run it?
The second code I shared is just for gravity deposition. It does not
contain the code for applying compression and extensions to the pack. I did
compression and extension cycles manually in the terminal after gravity
deposition. Each time a cycle is completed, I check if the porosity of the
pack reaches a constant value. If it does, I repeat the cycle again until
the porosity reaches the desired value. So there are no stop conditions in
the code.

On Tue, 22 Mar 2022 at 22:05, Jan Stránský <
<email address hidden>> wrote:

> Your question #701028 on Yade changed:
> https://answers.launchpad.net/yade/+question/701028
>
> Status: Open => Needs information
>
> Jan Stránský requested more information:
> > I tried scaling the packing, but that doesn't seem to eliminate the
> issue.
>
> if the issue is that the cylinder is not positioned correctly, then the
> solution is to position it correctly :-)
>
> > I've shared the coordinates file in the mail.
>
> launchpad questions does not work well with mail.. you can put it here
> directly,
>
> > Here is the code I used to create the sample:
>
> what are the stop conditions? what number of iterations to run it?
>
> Cheers
> Jan
>
> --
> To answer this request for more information, you can either reply to
> this email or enter your reply at the following page:
> https://answers.launchpad.net/yade/+question/701028
>
> You received this question notification because you asked the question.
>

> Why does this happen?
> ValueError: max() arg is an empty sequence

this means, that the argument of max function is empty. Try in Python: max([]).
max is used in this line:
xtremum = max(abs(s) for s in plot.data['s'])
meaning that at the time you call it, plot.data['s'] is empty.
You can have some checking beforehand or try-except construction to prevent the error.

Cheers
Jan

Hi Jan,
Sorry for the very delayed response. I might have solved that error. I
think the iteration interval that I gave for the command 'plotAddData' was
higher than the period given for other commands that use plotAddData. But
now another error occurred:

<FATAL ERROR> ThreadRunner:35 void yade::ThreadRunner::run(): Exception
occured:
PyRunner error.

COMMAND: 'plotAddData()'

ERROR:
list index out of range

STACK TRACE:
Traceback (most recent call last):

  File "<string>", line 1, in <module>

  File "Triax.py", line 134, in plotAddData
    e = (top[0].state.displ()[2] - bot[0].state.displ()[2]) /
(height-rParticle*2*bcCoeff)

IndexError: list index out of range

Why does this error occur? I don't understand meaning of some of these
errors. This is the updated 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(-27,-27,-51), scale=1.0, material=sphereMat)
spheres = O.bodies.append(sp)

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

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

#create top & bot facet plate
facets3 = []
nw=45
rCyl2 = (.6*width+2*rParticle) / cos(pi/float(nw))
for r in range(nw):
if r%2==0:
v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),
rCyl2*sin(2*pi*(r+0)/float(nw)), height )
v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),
rCyl2*sin(2*pi*(r+1)/float(nw)), height )
v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)),
rCyl2*sin(2*pi*(r+2)/float(nw)), height )
v4 = Vector3( 0, 0, height )
f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat)
f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat)
facets3.extend((f1,f2))
topcap = O.bodies.append(facets3)
facets3 = []
for r in range(nw):
if r%2==0:
v1 = Vector3( rCyl2*cos(2*pi*(r+0)/float(nw)),
rCyl2*sin(2*pi*(r+0)/float(nw)), 0 )
v2 = Vector3( rCyl2*cos(2*pi*(r+1)/float(nw)),
rCyl2*sin(2*pi*(r+1)/float(nw)), 0 )
v3 = Vector3( rCyl2*cos(2*pi*(r+2)/float(nw)),
rCyl2*sin(2*pi*(r+2)/float(nw)), 0 )
v4 = Vector3( 0, 0, 0 )
f1 = facet((v1,v2,v4),color=(0,0,0),material=frictMat)
f2 = facet((v2,v3,v4),color=(0,0,0),material=frictMat)
facets3.extend((f1,f2))
botcap = O.bodies.append(facets3)

# C.g). define top & bot wall id
for s in topcap:
top_id = s
bot_id = 0
for s in botcap:
bot_id = s

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

# run
O.run()

> You can have some checking beforehand or try-except construction to
prevent the error.
Also can you please explain this soln that you gave for the previous reply.

Thank you,
Rahul.

Maybe your bot or top list variables are empty, making it impossible to ask for either bot[0] or top[0]

Thanks Jan Stránský, that solved my question.