Applying confinement (rigid boundary) to a cylindrical sample

Asked by Swapnil on 2017-09-06

How can we apply a confining cylindrical boundary to a cylindrical
arrangement (packing) of the spheres ( in order to constrain the motion of
spheres outwards when a compressive load is applied) ? I have been searching but can only find information
about 'Wall'.

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2017-09-06
Last reply:
2019-01-09
Jan Stránský (honzik) said : #1

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

Swapnil (swapspace) said : #2

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.

Jan Stránský (honzik) said : #3

> 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

Swapnil (swapspace) said : #4

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.

Swapnil (swapspace) said : #5

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

Jan Stránský (honzik) said : #6

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

Swapnil (swapspace) said : #7

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

Jan Stránský (honzik) said : #8

Yes, it should be fine
Jan

Swapnil (swapspace) said : #9

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

Jan Stránský (honzik) said : #10

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

Swapnil (swapspace) said : #11

Understood this Jan, thanks a ton :)

trucgiao91 (giao) said : #12

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

Jan Stránský (honzik) said : #13

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

Can you help with this problem?

Provide an answer of your own, or ask Swapnil for more information if necessary.

To post a message you must log in.