Yade

Final velocity

Asked by Rodolfo França de Lima on 2017-10-30

I am in doubt as to the plot of a clump. I want to plot the speed of an object falling freely in different output positions. However in different positions the final speed is the same, as follows the scripts:

from yade import pack, plot, qt, utils
import sys
raio=0.0015
global angulo
angulo=0

materialsoja=CohFrictMat(density=788,young=10.9e6,poisson=0.3,alphaKr=0.05,frictionAngle=radians(31.8),label="milho")
idsoja=O.materials.append(materialsoja)

O.bodies.appendClumped([
#camada1
sphere([0,-2*raio,10*raio],raio,material=idsoja),
sphere([0,-1*raio,9.95*raio],raio,material=idsoja),
sphere([0,2*raio,10*raio],raio,material=idsoja),
sphere([0,1*raio,9.95*raio],raio,material=idsoja),
sphere([0,0,9.95*raio],raio,material=idsoja),

#camada2
sphere([0,-1.95*raio,9*raio],raio,material=idsoja),
sphere([0,1.95*raio,9*raio],raio,material=idsoja),
sphere([0,-0.95*raio,9*raio],raio,material=idsoja),
sphere([0,0.95*raio,9*raio],raio,material=idsoja),
sphere([0,0,9*raio],raio,material=idsoja),

#camada3
sphere([0,-1.9*raio,8*raio],raio,material=idsoja),
sphere([0,1.9*raio,8*raio],raio,material=idsoja),
sphere([0,-0.9*raio,8*raio],raio,material=idsoja),
sphere([0,0.9*raio,8*raio],raio,material=idsoja),
sphere([0,0,8*raio],raio,material=idsoja),

#camada4
sphere([0,-1.85*raio,7*raio],raio,material=idsoja),
sphere([0,1.85*raio,7*raio],raio,material=idsoja),
sphere([0,-0.85*raio,7*raio],raio,material=idsoja),
sphere([0,0.85*raio,7*raio],raio,material=idsoja),
sphere([0,0,7*raio],raio,material=idsoja),

#camada5
sphere([0,-1.8*raio,6*raio],raio,material=idsoja),
sphere([0,1.8*raio,6*raio],raio,material=idsoja),
sphere([0,-0.8*raio,6*raio],raio,material=idsoja),
sphere([0,0.8*raio,6*raio],raio,material=idsoja),
sphere([0,0,6*raio],raio,material=idsoja),

#camada6
sphere([0,-1.75*raio,5*raio],raio,material=idsoja),
sphere([0,1.75*raio,5*raio],raio,material=idsoja),
sphere([0,-0.75*raio,5*raio],raio,material=idsoja),
sphere([0,0.75*raio,5*raio],raio,material=idsoja),
sphere([0,0,5*raio],raio,material=idsoja),

#camada7
sphere([0,-1.7*raio,4*raio],raio,material=idsoja),
sphere([0,1.7*raio,4*raio],raio,material=idsoja),
sphere([0,-0.7*raio,4*raio],raio,material=idsoja),
sphere([0,0.7*raio,4*raio],raio,material=idsoja),
sphere([0,0,4*raio],raio,material=idsoja),

#camada8
sphere([0,-1.65*raio,3*raio],raio,material=idsoja),
sphere([0,1.65*raio,3*raio],raio,material=idsoja),
sphere([0,-0.7*raio,3*raio],raio,material=idsoja),
sphere([0,0.7*raio,3*raio],raio,material=idsoja),
sphere([0,0,3*raio],raio,material=idsoja),

#camada9
sphere([0,-1.6*raio,2*raio],raio,material=idsoja),
sphere([0,1.6*raio,2*raio],raio,material=idsoja),
sphere([0,-0.7*raio,2*raio],raio,material=idsoja),
sphere([0,0.7*raio,2*raio],raio,material=idsoja),
sphere([0,0,2*raio],raio,material=idsoja),

#camada10
sphere([0,-1.55*raio,raio],raio,material=idsoja),
sphere([0,1.55*raio,raio],raio,material=idsoja),
sphere([0,-0.7*raio,raio],raio,material=idsoja),
sphere([0,0.7*raio,raio],raio,material=idsoja),
sphere([0,0,raio],raio,material=idsoja)
])

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   DragEngine(ids=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40],Cd=0.6, Rho=1.225),
   NewtonIntegrator(damping=0.05,gravity=(0,0,-9.81),label="NI"),
   PyRunner(command='para()',iterPeriod=10000),
   PyRunner(iterPeriod=1000,command="plotAddData()"),

]

def para():
if O.time>=7:
O.pause()

O.dt=utils.PWaveTimeStep()

O.saveTmp()

#O.run()

def plotAddData():
   v = O.bodies[0].state.vel
   vz = v[2]
   plot.addData(
      iter = O.iter,
      time = O.time,
      vz = vz,
   )
plot.plots = {'time':'vz'}
plot.plot()

and now modificad

from yade import pack, plot, qt, utils
import sys
raio=0.0015
global angulo
angulo=0

materialsoja=CohFrictMat(density=788,young=10.9e6,poisson=0.3,alphaKr=0.05,frictionAngle=radians(31.8),label="milho")
idsoja=O.materials.append(materialsoja)

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
   InteractionLoop(
      [Ig2_Sphere_Sphere_L3Geom(),Ig2_Facet_Sphere_L3Geom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_L3Geom_FrictPhys_ElPerfPl()]
   ),
   DragEngine(ids=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,
25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40],Cd=0.6, Rho=1.225),
   NewtonIntegrator(damping=0.05,gravity=(-9.81,0,0),label="NI"),
   PyRunner(command='para()',iterPeriod=10000),
   PyRunner(iterPeriod=1000,command="plotAddData()"),

]

def para():
if O.time>=7:
O.pause()

O.dt=utils.PWaveTimeStep()

O.saveTmp()

#O.run()

def plotAddData():
   v = O.bodies[0].state.vel
   vz = v[0]
   plot.addData(
      iter = O.iter,
      time = O.time,
      vz = vz,
   )
plot.plots = {'time':'vz'}
plot.plot()

I believe that Yade is not calculating the projection area of the whole clump, but of only one sphere. Is it possible to make the Yade calculate velocities considering the different potions of the falling particle?

Thank you very much.

Question information

Language:: English Edit question

Status:: Solved

For:: Yade Edit question

Assignee:: No assignee Edit question

Solved by:: Rodolfo França de Lima

Solved:: 2017-11-01

Last query:: 2017-11-01

Last reply:: 2017-11-01

Link existing bug

Revision history for this message

Jan Stránský (honzik) said on 2017-10-31:

Hello,

first of all, thanks for the code you are using.
However, I am sorry, but I did not understand your problem (or more specifically, the results makes perfect sense to me). Please try to elaborate a bit more the problem.

> I want to plot the speed of an object falling freely in different output positions.

what does "different output positions" mean? Just that the gravity is applied perpendicularly in the two cases, or something else?

> I believe that Yade is not calculating the projection area of the whole clump, but of only one sphere.

I am sorry, here I did not get at all the part with projection area..

> Is it possible to make the Yade calculate velocities considering the different potions of the falling particle?

If I got it correctly, then no, see below

> However in different positions the final speed is the same, as follows the scripts:

If I got it correctly, the final speeds SHOULD be the same :-) you apply the same gravity (just rotated) to the same object, so its motion is governed by the same equation
acceleration = -9.81
so I see no way how the final velocities could differ..

thanks
Jan

Revision history for this message

Rodolfo França de Lima (rodolfolima) said on 2017-10-31:

I will try to explain it better. I am using the dragengine in YADE, so that I’m interested in the terminal velocity of the body. The terminal velocity depends on the mass, the drag coefficient and the projected area. With simple spheres, Yade calculates it right, but my problem is with clumps. I have a clump composed of some spheres, if it is falling sideways, the terminal velocity will be smaller than when falling pointing to the ground. It happens because of 2 factors, because the projected area “like the shadow”, will be smaller, and also because the drag coefficient would change. Even if I keep the same drag coefficient, the terminal velocity will be different because of the change in the projected area. When I simulate the clump in YADE, it does not happens. I suspect that it is because YADE just computes the drag force for each independent sphere, and not the clump. Is there a way to simulate the drag force at the clump level, and not just for the spheres as if they were not part of a clump?

Revision history for this message

Jan Stránský (honzik) said on 2017-10-31:

Hi Rodolfo,
thanks for more detailed info.
I am sorry, in the code I totally missed the DragEngine, it is very important related to your question, next time please emphasize it :-)
In this case, you are right: the current implementation [1] only considers individual spheres, not clumps.

> I believe that Yade is not calculating the projection area of the whole clump, but of only one sphere.

yes :-)

> Is it possible to make the Yade calculate velocities considering the different potions of the falling particle?
> Is there a way to simulate the drag force at the clump level, and not just for the spheres as if they were not part of a clump?

not with current implementation of DragEngine. But I think it would not be difficult to write your own DragEngine for clumps in Python, we can help you start if you want. So yes, it is possible :-)

cheers
Jan

[1] https://github.com/yade/trunk/blob/master/pkg/common/ForceEngine.cpp#L45

Revision history for this message

Bruno Chareyre (bruno-chareyre) said on 2017-10-31:

It seems a question is "does DragEngine support clumps?".
I don't know (I guess not), but you can probably check the code to see how the equations are.

There is also maybe an inconsistency between what you expect and what your script does:
>I suspect that it is because YADE just computes the drag force for each independent sphere, and not the clump.

The engine does what you ask: "DragEngine(ids=[0,1,2,3,4,5,6,7,8,9..." means that you declare the clump members to the engine. I would think you just need to declare the clump.
The numbering method is also unclear. You have 51 bodies after first section (50 members + 1 clump) and 40 DragEngine ids. Why?

Bruno

Revision history for this message

Rodolfo França de Lima (rodolfolima) said on 2017-10-31:

Thank you for the help, I will look into the implementation and maybe implement a custom dragengine.

Revision history for this message

Jan Stránský (honzik) said on 2017-11-01:

Hi Rodolfo,

based on personal communication, please find attached an example.
Drag coefficient can be also implemented as tensor if needed.
Before relying on this engine, please properly test it.(!!!)

You have to define area yourself, in local coordinate system (with respect to principal axes of inertia). But it also can be automated.
If performance would be a problem, C++ implementation is also possible.

cheers
Jan

####
from yade import plot

class PyDragEngine:
  def __init__(self,rho,cd,area,id):
    """
    Python implementation of
    rho ... mass density of fluid
    cd ... drag coefficient
    area ... "area tensor" in local coordinate system
    id ... body id to apply force on

    TODO cd could also be a tensor (as area)
    TODO testing
    """
    self.rho = float(rho)
    self.cd = float(cd)
    self.area = Matrix3(area)
    self.id = int(id)
  def __call__(self):
    b = O.bodies[self.id]
    v = b.state.vel.norm() # velocity magnitude
    if abs(v) < 1e-6: # do nothing for very low velocity
      return
    n = b.state.vel / v # velocity direction as unit vector
    r = b.state.ori.toRotationMatrix() # rotation matrix
    a = r*self.area*r.transpose() # rotation of area tensor according to rot of body. Maybe it is rT*a*r, please check
    a = n.dot(a*n) # projection of "a" to direction of n.
    f = .5*self.rho*self.cd*a*pow(v,2) # drag force magnitude
    f = -f*n # drag force vector
    O.forces.addF(self.id,f) # add it to the body

cid,sids = O.bodies.appendClumped(( # testing clump, ax=2, ay=3, az=6
  sphere((0,0,0),.5),
  sphere((1,0,0),.5),
  sphere((2,0,0),.5),
  sphere((0,1,0),.5),
  sphere((1,1,0),.5),
  sphere((2,1,0),.5),
))
clump = O.bodies[cid]
O.dt = 0
clump.state.ori = Quaternion((1,0,0),.5*pi)*Quaternion((0,0,1),.5*pi) # to test with non-unit ori, then ax=3, ay=6, az=2

a = Matrix3(3,0,0, 0,6,0, 0,0,2) # global area tensor, diag(ax,ay,az)
r = clump.state.ori.toRotationMatrix()
a = r.transpose()*a*r # local area tensor # ?? r*a*rT ??
pyDragEngine = PyDragEngine(rho=1.225,cd=.6,area=a,id=cid)

O.engines=[
  ForceResetter(),
  InsertionSortCollider([Bo1_Sphere_Aabb()]),
  InteractionLoop(
    [Ig2_Sphere_Sphere_ScGeom()],
    [Ip2_FrictMat_FrictMat_FrictPhys()],
    [Law2_ScGeom_FrictPhys_CundallStrack()]
  ),
  PyRunner(iterPeriod=1,command="pyDragEngine()"), # before NowtonIntegrator (!)
  NewtonIntegrator(damping=0.,label='newton'),
  PyRunner(iterPeriod=1,command="plotAddData()"),
]

def plotAddData():
plot.addData(t=O.time,v=clump.state.vel.norm())

# test for different gravity directions
O.dt = 3e0
newton.gravity = Vector3(1,0,0) # ax=3
O.run(100,True)
calm()
newton.gravity = Vector3(0,1,0) # ay=6
O.run(100,True)
calm()
newton.gravity = Vector3(0,0,1) # az=2
O.run(100,True)
calm()
newton.gravity = Vector3(1,1,1).normalized()
O.run(100,True)

plot.plots = {'t':'v'}
plot.plot()
####

Hi Rodolfo,

based on personal communication, please find attached an example.
Drag coefficient can be also implemented as tensor if needed.
Before relying on this engine, please properly test it.(!!!)

cheers
Jan

####
from yade import plot

class PyDragEngine:
  def __init__(self,rho,cd,area,id):
    """
    Python implementation of 
    rho ... mass density of fluid
    cd ... drag coefficient
    area ... "area tensor" in local coordinate system
    id ... body id to apply force on

def plotAddData():
  plot.addData(t=O.time,v=clump.state.vel.norm())

plot.plots = {'t':'v'}
plot.plot()
####

Revision history for this message

Rodolfo França de Lima (rodolfolima) said on 2017-11-01:

Thanks a lot for the help, this will be very useful. I will share any improvement I eventually make, like area calculation.

To post a message you must log in.

Ask a question

Edit question

Yade

Final velocity

Question information

Related bugs

Related FAQ:

Subscribers