Applying confinement (rigid boundary) to a cylindrical sample

Hello,
you can use walls or facets or boxes. See e.g. [1] for static boundary or [2] for active boundary applying lateral confining pressure.
cheers
Jan

###
from yade import pack, geom
pred = pack.inCylinder((0,0,0),(0,0,0.25),0.10)
sp = pack.randomDensePack(pred,radius=0.01,spheresInCell=300)
O.bodies.append(sp)
facets = geom.facetCylinder((0,0,0.125),0.1,0.25,wallMask=2|4)
O.bodies.append(facets)
###

[1] https://yade-dem.org/doc/yade.geom.html#yade.geom.facetCylinder
[2] https://github.com/yade/trunk/blob/master/examples/concrete/triax.py

Thanks Jan.
Option [1] leads to the movement of spheres out of the boundary (I tried the script). I do not want the spheres to move out of the boundary. Therefore, I would have to go for application of lateral confining pressure. That you say is depicted in [2].

# facets division
nw = 24,
nh = 15,

# 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))
O.bodies.append(facets)
mass = O.bodies[0].state.mass
for f in facets:
 f.state.mass = mass
f.state.blockedDOFs = 'XYZz'

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

As I can understand from reading this:

Using [1] creates a default set of facets and combines them to form a cylindrical boundary of the prescribed dimensions. Isn't it?
On the other hand in [2], the user is creating facets (using vectors) and then combining them followed by application of pre-defined stresses(confining pressure). My question here is can this operation be performed like in [1]? that is using a simpler way to create facets and then applying stresses.

For e.g something like (just an idea to keep it simple):

#create facet sets and combine them to form a cylindrical shaped boundary around the sample.

facets = geom.facetCylinder((0,0,0.125),0.1,0.25,wallMask=2|4)
O.bodies.append(facets)

# for all facets block restrain the desired motion
facets.state.blockedDOFs='XYZz'

#apply prestress to each one
O.forces.addF(f.id,preStress*a*n)

In short, I am looking for a simpler way to perform the application of laterally confining pressures.

> I do not want the spheres to move out of the boundary

in both cases, you will have to use different engines, namely Bo1_Facet, Ig2_Facet_Sphere, see [2]

> Using [1] creates a default set of facets and combines them to form a cylindrical boundary of the prescribed dimensions. Isn't it?

yes

>On the other hand in [2], the user is creating facets (using vectors) and then combining them followed by application of pre defined stresses(confining pressure). My question here is can this operation be performed like in [1]? that is using a simpler way to create facets and then applying stresses.

yes, you can use facetCylinder and then apply force on resulting facets in the same way as in [2]

> In short, I am looking for a simpler way to perform the application of laterally confining pressures.

Do you really want active confining pressure? Wouldn't static boundary be enough? At least for beginning, I would try the static boundary..

cheers
Jan

Yes mate it works with the inclusion of those engines. Foolishly, I completely forgot about those for the new facet bodies.
Thanks :)

Yes, I am first trying the static boundary case ( where the material of the facetCylinder itself holds the spheres intact). Thanks to you, that appears to work successfully now.

As you are suggesting, for the active confining pressure case we can directly use the facetCylinder as one body and apply the desired stress, isn't it? Then should I try the following in addition to [1]:

 facets.state.blockedDOFs='XYZz' --> this command seems okay

O.forces.addF(f.id,preStress*a*n) --> I am confused about the application of forces using the facetcylinder as a whole. The (f.id,Prestress*a*n) used here is when we apply forces through each facet. In order to apply force from the facetcylinder as a whole (such that each facet does the same) can I use: O.forces.addF(facets.id,0.20) where, for e.g., 0.20 represents the magnitude of force(stress) applied.

Would be trying the rest on Monday when I return to the system ;)

Hi,

> O.forces.addF(f.id,preStress*a*n) --> I am confused about the application of forces using the facetcylinder as a whole. The (f.id,Prestress*a*n) used here is when we apply forces through each facet.

You can apply force to the cylinder as a whole, but then it would not apply confining pressure, but tried to move as a rigid body in the direction of the force, so very different result than confining boundary..

> In order to apply force from the facetcylinder as a whole (such that each facet does the same) can I use: O.forces.addF(facets.id,0.20) where, for e.g., 0.20 represents the magnitude of force(stress) applied.

no
- facets.id would give error as facets is a list of bodies, not body.
- O.forces.addF accepts vector, not only magnitude

Jan

Alright, got it :)

It implies that after I create a facet cylinder using the command:
 facets = geom.facetCylinder((0,0,0.125),0.1,0.25,wallMask=2|4)
O.bodies.append(facets)

I need to extract each facet that forms the cylinder and then apply stress to it. Isn't it?

I guess I have to perform the following.
preStress= -2e4
for f in facets:
f.state.blockedDOFs = 'XYZz'
def addForces():
 for f in facets:
        n = f.shape.normal
        a = f.shape.area
        O.forces.addF(f.id,preStress*a*n)
----
and then it should be fine .Please correct me if the above deduction is wrong :)

Yes, it should be fine
Jan

Thanks Jan. I tried running the script. It does not appear to be working the way it should. Instead, I see a few spheres and a different interaction behavior. Can you try running the script to figure out what is being done wrong here :)

from yade import pack,plot
pred=pack.inCylinder((0,0,0),(0,0,0.25),0.10)
sp=pack.randomDensePack(pred,radius=0.01,spheresInCell=300)
O.bodies.append(sp)
facets = geom.facetCylinder((0,0,0.125),0.1,0.25,wallMask=2|4)
O.bodies.append(facets)
yade.qt.Controller()
idSteel=O.materials.append(FrictMat(young=210e9,poisson=.25,frictionAngle=.8,label="steel"))
hammer=sphere((0,0,0.30),0.04,material=idSteel)
hammerID=O.bodies.append(hammer)
hammer.state.blockedDOFs='z'
calm()
hammer.state.vel=Vector3(0,0,-0.5)
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),NewtonIntegrator(),PyRunner(iterPeriod=100,command="addPlotData()"),PyRunner(iterPeriod=100,command="addForces()")]
preStress= -2e4
for f in facets:
  f.state.blockedDOFs = 'XYZz'
def addForces():
  for f in facets:
    n = f.shape.normal
    a = f.shape.area
    O.forces.addF(f.id,preStress*a*n,permanent=True)
def addPlotData():
  Dz=hammer.state.displ()[2]
  Fz=O.forces.f(hammerID)[2]
  plot.addData(i=O.iter,Fz=Fz,Dz=Dz)
plot.plots={'Dz':('Fz')}
plot.plot()
O.dt=0.5*PWaveTimeStep()
O.saveTmp()
O.run()

Hello,

in this case, try to find the problem step by step. I.e. do not run directly (remove O.run() form the script) and have a look what is hapenning. You would see perfect state at the beginning, but after one step, all facets disappears. In a more complicated configuration, such information could lead to faster solution.

In this case, the solution is easy. Facets has zero mass and inertia by default, leading to NaNs (and then disappearance) because of accel=force/mass

##
...
for f in facets:
 f.state.blockedDOFs = 'XYZz'
 f.state.mass = 1 # possible a different value
 f.state.inertia = (1,1,1) # possible a different value
...
##

cheers
Jan

Understood this Jan, thanks a ton :)

Hi, I am also new to DEM, is it possible if I can just block particles from the outer layer not to move in the x and y direction?
regards,
G

Hello,

> I am also new to DEM

welcome :-)

> is it possible if I can just block particles from the outer layer not to move in the x and y direction?

b.state.blockedDOFs = "xy"

please next time open a new question and do not reopen long time inactive questions

cheers
Jan