Oedometric consolidation test

Hello,

(5) depends where and how you want to use it

(6) An alternative is to use facets only for the cylinderical (non-flat) surface and used boxes for bottom and top (actually you can also use boxes for the cylinder itself). You can then easily load the top box as a single body (instead of dealing with a bunch of facets). In that case, you would have to modify engines to reflect Box shape

> if O.time < 3600 ( say the time for the experiment is 1 hour)

It is easy to estimate "simulation speed" from O.dt and number of time steps per real second. I am not sure that simulating one hour is realistic..

cheers
Jan

thanks Jan :)

For (5) will have to explore while running the simulation

      (6) that brings me to a new question. How to use boxes to create a cylindrical geometry? Then pick out the top box and define it as the loading plate.

(6) you can use this or similar function to create a cylindrical surface. Top and bottom box is easy. Have a look at [1] for engines to use
###
from math import sin,cos

def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):
   """creates cylinder made of boxes. Axis is parallel with z+"""
   center = Vector3(center)
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles]
   ret = []
   for i,pt in enumerate(pts):
      l = pi*radius/nSegments
      es = (.5*thick,l,.5*height)
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      ret.append(box(pt,es,ori,**kw))
   return ret

surf = cylSurf((1,2,3),4,5,nSegments=12,thick=.1)
O.bodies.append(surf)
###

cheers
Jan

[1] https://github.com/yade/trunk/blob/master/examples/simple-scene/simple-scene.py

Thanks Jan but I did not fully understand the code you have written here. Only grabbed the overview. Could you please comment on each step above for a clearer understanding and implementation.

Upon implementing the above code, I am curious as to how we can achieve a circular box at the top and bottom of the cylinder. I am searching for the commands in the documentation but cannot find any.

###
def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):
   """creates cylinder made of boxes. Axis is parallel with z+
   center ... center of bottom disc
   radius ... radius of cylinder
   height ... height of cylinder
   nSegments ... number of boxes forming cylinder
   thick ... thickness of the boxes
   **kw ... keywords passed to box function (e.g. material, color ...)
   """
   # just converts (a,b,c) to Vector3(a,b,c) to enable "center + ..."
   center = Vector3(center)
   # nSegments angles to define individual boxes, from 0 to 2*pi
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   # centers of boxes. Vector 3(radius,0,0) rotated by angle 'a'
   # z coordinate is .5*height
   # center + shift = center of corresponding box
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles] # centers of plates
   ret = []
   for i,pt in enumerate(pts): # for each box center create corresponding box
      # "length" of the box along circumference, cylinder circumference / nSegments
      l = pi*radius/nSegments
      # extents, half dimensions along x,y,z axis
      es = (.5*thick,l,.5*height)
      # by default box is axis-aligned, we need to rotate it along z axis (0,0,1) by angle i*2*pi/nSegments
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      # have a look at utils.box documentation. pt=center of box, es=extents (half dimension), ori=orientation
      ret.append(box(pt,es,ori,**kw))
   return ret
###

you can use print to see the values if you want to understand what is going on. You can also write your function from the scratch.
cheers
Jan

Alright Jan. Thanks, now I have some idea about this :)

But I still cannot figure out the way to create circular top and bottom boxes (plate shapes) from the documentation

Boxes are just boxes, cuboids. The point is that you do not need circular top and bottom, the enlarged "square" box is just fine, the overhangs do not matter.
cheers
Jan

Yes, Jan it works fine so far. I tried the following:
upper_plate=O.bodies.append(box((1,2,3),(4,4,0.0002)))
lower_plate=O.bodies.append(box((1,2,8),(4,4,0.0002)))

The set-up looks fine. My only confusion is whether due to its square geometry it would move into the cylinder upon loading or not? I mean it should behave similar to a round plate in function (as you are suggesting that the overhangs do not matter)

I will try that part along with the rest tomorrow :)

> whether due to its square geometry it would move into the cylinder upon loading or not?

it would move as it was circular, because there would be no box-box interaction (you would not have Ig2_Box_Box.. in engines, because it is not (yet) implemented in Yade).

cheers
Jan

Okay, I will be on that part soon :)

firstly, I am trying to fill in the created cylinder with spheres. I did the following (full script):
# """creates cylinder made of boxes. Axis is parallel with z+""
from math import sin,cos

def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):

   center = Vector3(center)
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles]
   ret = []
   for i,pt in enumerate(pts):
      l = pi*radius/nSegments
      es = (.5*thick,l,.5*height)
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      ret.append(box(pt,es,ori,**kw))
   return ret

surf = cylSurf((1,2,3),4,5,nSegments=12,thick=.1)
O.bodies.append(surf)

# create the lower plate
lower_plate=O.bodies.append(box((1,2,3),(4,4,0.0002)))

#pack the body (named "surf" with spherical particles)

sp=pack.randomDensePack(surf,radius=0.01,spheresInCell=300)
O.bodies.append(sp)
yade.qt.Controller()

This gives the following error:
File "/usr/bin/yade", line 182, in runScript
    execfile(script,globals())
  File "box-cylinder.py", line 19, in <module>
    sp=pack.randomDensePack(surf,radius=0.01,spheresInCell=300)
  File "/usr/lib/x86_64-linux-gnu/yade/py/yade/pack.py", line 418, in randomDensePack
    if not dim: dim=predicate.dim()
AttributeError: 'list' object has no attribute 'dim'

Is there another method to pack the body with a set of spheres?

use pack.inCylinder(...) as a predicate for randomDensePack, not surf (list of boxes)
Jan

I tried
pred=pack.inCylinder((1,2,3),(1,2,8),4)
sp=pack.randomDensePack(pred,radius=0.01,spheresInCell=300)

and the system is overloaded and hung. Cannot see the simulation :)
Is there another way of packing other than the randomDense. Like where we can pack particles with a certain porosity.

very roughly: you have volume with dimensions 8x8x5, trying to fill there particles with 0.02 diameter, i.e. 400x400x250=40M particles.. probably that is the reason of "the system is overloaded and hung".

> Is there another way of packing other than the randomDense. Like where we can pack particles with a certain porosity.

randomDensePack is a quick predefined function. You can have a look at the source code and create your own function based on your needs (e.g. frictionAngle of used material during compaction has influence on resulting porosity). The idea is simple, use makeCloud to create a loose packing and then somehow compact it.
There are also geometrical approaches to create a packing, but I am not very familiar with them. If you are interested in this topic, please open a new question.

cheers
Jan

Okay Jan, lets keep the porosity topic for later.

Now, I managed to get the set up ready. the following is the script:

from yade import pack,plot

(1) create cylinder made of boxes.

from math import sin,cos

def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):

   center = Vector3(center)
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles]
   ret = []
   for i,pt in enumerate(pts):
      l = pi*radius/nSegments
      es = (.5*thick,l,.5*height)
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      ret.append(box(pt,es,ori,**kw))
   return ret

surf = cylSurf((1,2,3),4,5,nSegments=12,thick=.1)
O.bodies.append(surf)

(2) create a bottom plate
lower_plateID=O.bodies.append(box((1,2,3),(4,4,0.0002)))

(3) Fill in the cylinder with spjheres
pred=pack.inCylinder((1,2,3),(1,2,8),4)
sp=pack.randomDensePack(pred,radius=0.25)
O.bodies.append(sp)
yade.qt.Controller()

(4) create an upper (top plate)
upper_plateID=O.bodies.append(box((1,2,8),(4,4,0.0002)))

(5) Load the top plate with a certain constant force
O.forces.addF(upper_plateID,(0,0,-0.02),permanent=True)

(6) Apply Engines for interactions
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),NewtonIntegrator(gravity=[0,0,-9.81],damping=0.1)]

((7) time step
O.dt= 0.3*PWaveTimeStep()

At this point when I run the simulation, it is blowing up (exploding). (Please try the run on your system to cross-check)
Is there some error in my process -implementation so far :) ?

a) use
surf = cylSurf( ..., wire=True)
to have transparent boxes to see what is hapenning inside. In this case, the shape is strange and spheres are overlapping, causing explosion

b) I would use
sp=pack.randomDensePack(..., spheresInCell=300)
I have almost no experience without spheresInCell, and that I have is a bad experience :-)

c) make the boxes (except the one that is loaded) fixed
for b in surf: b.state.blockedDOFs = 'xyzXYZ'
O.bodies[lower_plateID].state.blockedDOFs = 'xyzXYZ'

Jan

Thanks Jan, that was right on the current target ;)
Coming to the next part of the problem which involves the plotting of the data, I wish to track the displacement of the upper plate w.r.t the time (the square-root of time to be precise). For this purpose I have come up with the following:

adding to the engines: PyRunner(iterPeriod=100,command="addPlotData()") followed by:

import math
     def addPlotData():
           Dz=plate.state.displ()[2]
           t= O.time
           plot. addData(Tx=math.sqrt(t), Dz=Dz)
     plot.plots= { 'Tx':('Dz')}
plot.plot()

I can see the plot but upon running the simulation the data is not updated on the plot. What am I doing wrong.. :) ?

correction: I have used, Dz=upper_plateID.state.displ()[2]

Hello,
please post a complete code, the segment you gave should be a syntax error..
thanks
Jan

Hi Jan. Sorry for being late. I was not at the system for all these days.

Here is the current code:

from yade import pack,plot
import math,time
from math import sin,cos

def cylSurf(center,radius,height,nSegments=12,thick=0,**kw):
   """creates cylinder made of boxes. Axis is parallel with z+"""
   center = Vector3(center)
   angles = [i*2*pi/nSegments for i in range(nSegments)]
   pts = [center + Vector3(radius*cos(a),radius*sin(a),.5*height) for a in angles]
   ret = []
   for i,pt in enumerate(pts):
      l = pi*radius/nSegments
      es = (.5*thick,l,.5*height)
      ori = Quaternion((0,0,1),i*2*pi/nSegments)
      ret.append(box(pt,es,ori,**kw))
   return ret

surf = cylSurf((1,2,3),4,5,nSegments=12,thick=.1,wire=True)
O.bodies.append(surf)
for b in surf: b.state.blockedDOFs='xyzXYZ'
pred=pack.inCylinder((1,2,3),(1,2,8),4)
sp=pack.randomDensePack(pred,radius=0.25,spheresInCell=350)
O.bodies.append(sp)
yade.qt.Controller()
lower_plateID=O.bodies.append(box((1,2,3),(4,4,0.0002)))
O.bodies[lower_plateID].state.blockedDOFs = 'xyzXYZ'
upper_plateID=O.bodies.append(box((1,2,8),(4,4,0.0002)))
O.forces.addF(upper_plateID,(0,0,-0.02),permanent=True)
O.engines=[ForceResetter(),InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),InteractionLoop([Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],[Ip2_FrictMat_FrictMat_FrictPhys()],[Law2_ScGeom_FrictPhys_CundallStrack()]),PyRunner(iterPeriod=100,command="addPlotData()"),NewtonIntegrator(gravity=[0,0,-9.81],damping=0.1)]
O.dt= 0.3*PWaveTimeStep()

def addPlotData():
  Dz=upper_plate.state.displ()[2]
  t= O.time
  plot.addData(Tx=math.sqrt(t), Dz=Dz)

plot.plots= { 'Tx':('Dz')}
plot.plot()

Once again, take care of errors

NameError: global name 'upper_plate' is not defined

cheers
Jan

Jan, I think I got my mistake on that one ;)

Should have defined Upper_plate and then saved its ID separately.

Yes, exactly..

The last part of this task before moving to the next:

the box at the top (upper_plate) would have a default mass. Right? How can we know that value?

thereafter, I can alter it using the command: upper_plate.state.mass= ....( some suitable value).......

> the box at the top (upper_plate) would have a default mass. Right?

depending on "default". It has mass equal to its material density times its volume

> How can we know that value?

as you wrote:
m = upper_plate.state.mass
upper_plate.state.mass = whatever

cheers
Jan

Yeah, the reason I am asking is because it becomes unstable after a certain period of load application. Therefore, I am playing with the mass to find a suitable value to make the loading stable :)

It is one option. Another (perhaps "more physical") is to play with time step.
Jan

Sorry, in the case of a plate, probably playing with mass is ok also from "physical point of view" :-)
good luck
Jan

Or I guess doing both is an even better idea. Will try slowing down the simulation speed :)

I must point out (sorry to enter this discussion late) that you seem to make your life complex for no good, for the sole reason that you want a circular cross-section.
Proof: after 29 messages it is still not solved ;)
More precisely:
- pack.inCylinder(...) + randomDensePack will give extremely large porosity along the curved wall, which implies another equilibration stage, and still a highly heterogeneous micro-structure in the end.
- you have no automatic servo-control for the applied stress, so you replace it with an applied force on top plate, which is not the best choice in terms of stability (then you'll need to have artificial damping, never a good thing).
- calculating the radial stress is less easy (unless you are only interested in the sample-scale average)
- you cannot use the fluid modules (be it 1 or 2-phase flow) because they assume a square box.

On the other hand, all the above problems disappear if you use TriaxialStressController as in many examples. You can even find example scripts for oedometer in both 1 and 2-phase flow.
Since there is nothing that distinguish a square cross section from a circular cross section in an oedometer (from a mechanical point of view), I would seriously consider throwing away the idea of a "realistic" mold shape unless there is _very_ strong reason for it. "Being like in experiment" is not a good reason IMO.

Bruno

Hi Bruno,

absolutely no problems with entering the discussion late. After all such discussions with learned people like you all here would make me understand and learn more :)

 Yes, I reckon you are right with trying the square box shape. I am yet to get the idea on the Triaxial stress controller though.

Mechanics-wise currently I am trying to get in the working with wet-granular media and hence have posted a separate question regarding the same :)

I'll come back to what you suggested here for sure..