calculate the mass of the particles

Asked by Estefany Carmona Alvarez on 2020-02-18

I have generated a code for the simulation of the sand free fall.

This simulation consists of dropping the spheres that simulate the sand, so that they are deposited inside two cylindrical vessels. As a result of this simulation I need to obtain the weight or mass of the particles that remain within each cylinder.

In the Yade documentation that I have consulted, I saw that the state.mass function fulfills this need, but since I do not have enough domain in this programming language, I needed to know if anyone has any example of how to use this code or if they have more information ...

I appreciate your collaboration

The code I am using is as follows:
postscript: the code comments are in Spanish (Latin)

from yade import pack, plot

#Altura de caida del material
hc=0

#EQUIPO

#coordenadas centro
x=0
y=0
z=0

#coordenadas extension

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion
O.bodies.append(geom.facetBox((x,y,z),(f,a,(h+hc)),wallMask=31,wire=True))

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

      #O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Radio
radius=0.05

#Altura de la probeta
height=0.1

#Probeta 1
#Coordenadas centro
xp1=0.1
yp1=0
zp1=1

O.bodies.append(geom.facetCylinder((xp1,yp1,(-zp1-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0.933,0.090,.894), wire=False))

#Probeta 2
#Coordenadas centro
xp2=0.1
yp2=0
zp2=1

O.bodies.append(geom.facetCylinder((-xp2,yp2,(-zp2-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0,0,1), wire=False))

#coordenadas perimetro de particulas
#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

#porosidad
p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.008
t2=0.015
t3=0.030
t4=0.085
t5=0.200
t6=0.475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
angulo=radians(36) #radianes
#densidad del material
densidad=1418 #(kg/m3)
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
porosidad=0.9
#Algoritmo YADE para el material
snd=O.materials.append(FrictMat(young=young,frictionAngle=radians(angulo),density=densidad,
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
    Bo1_Facet_Aabb(),
    Bo1_Cylinder_Aabb(),
    Bo1_Wall_Aabb()
 ]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 #GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner()
]
O.dt=PWaveTimeStep()

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-02-24
Last reply:
2020-02-24
Jan Stránský (honzik) said : #1

Hello,

I have tried your code, but the spheres are larger than the cylinders and therefore no spheres are deposited inside any cylinder and their mass is zero.
Please update the code such that it results in something meaningful.

Otherwise the task is easy, something like:
###
cyl1 = pack.inCylinder(...) # [1], predicate for a cylinder
mass = 0.0
for b in O.bodies:
   if cyl1(b.state.pos): # b is inside cyl1
      mass += b.state.mass
###

or a one-liner :-)
mass = sum(b.state.mass for b in O.bodies if cyl1(b.state.pos))

cheers
Jan

[1] https://yade-dem.org/doc/yade.pack.html#yade._packPredicates.inCylinder

1. In this code is the actual size of the particles, which if they are deposited in the cylinders ...

from yade import pack, plot

#Altura de caida del material
hc=0

#EQUIPO

#coordenadas centro
x=0
y=0
z=0

#coordenadas extension

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion
O.bodies.append(geom.facetBox((x,y,z),(f,a,(h+hc)),wallMask=31,wire=True))

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

      #O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Radio
radius=0.05

#Altura de la probeta
height=0.1

#Probeta 1
#Coordenadas centro
xp1=0.1
yp1=0
zp1=1

O.bodies.append(geom.facetCylinder((xp1,yp1,(-zp1-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0.933,0.090,.894), wire=False))

#Probeta 2
#Coordenadas centro
xp2=0.1
yp2=0
zp2=1

O.bodies.append(geom.facetCylinder((-xp2,yp2,(-zp2-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0,0,1), wire=False))

#coordenadas perimetro de particulas
#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

#porosidad
p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.008
t2=0.015
t3=0.030
t4=0.085
t5=0.200
t6=0.475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
angulo=radians(36) #radianes
#densidad del material
densidad=1418 #(kg/m3)
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
porosidad=0.9
#Algoritmo YADE para el material
snd=O.materials.append(FrictMat(young=young,frictionAngle=radians(angulo),density=densidad,
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
    Bo1_Facet_Aabb(),
    Bo1_Cylinder_Aabb(),
    Bo1_Wall_Aabb()
 ]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 #GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner()
]
O.dt=PWaveTimeStep()

2. Could you please help me with more information about the function you indicated to calculate the weight of the particles that remain inside the cylinder ...

###
cyl1 = pack.inCylinder(...) # [1], predicate for a cylinder
mass = 0.0
for b in O.bodies:
   if cyl1(b.state.pos): # b is inside cyl1
      mass += b.state.mass
###

or a one-liner :-)
mass = sum(b.state.mass for b in O.bodies if cyl1(b.state.pos))

I am new to Yade and in the available documents there is not much information about it.

Thanks for your help Jan !!

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

> I am new to Yade and in the available documents there is not much information about it.

Yade uses Python as an user interface. In python, you have access to almost any data in the simulation (tutorial [1] and user's manual [2] is full of this).
So you have access to positions of all particles, then it is easy to check if they are in a cylinder or not and sum mass of those which are. There are several options e.g.
###
def isBodyInCylinder(body,center,radius):
   d = body.state.pos - center # vector between cylinder center and particle center
   d[2] = 0 # only xy components are relevant
   distance = d.norm() # horizontal distance
   return distance < radius # particle is inside, if distance < raidus
def masses():
   mass1 = mass2 = 0.0 # intially zero
   center1 = Vector3(xp1,yp1,0) # center of 1st cylinder
   center2 = Vector3(-xp2,yp2,0) # center of 2nd cylinder
   for b in O.bodies: # test all particles
      if not isinstance(b.shape,Sphere): continue # skip non-spherical
      if isBodyInCylinder(b,center1,radius): # if a particle is inside 1st cylinder ...
         mass1 += b.state.mass # ... increase mass1 by the particle mass
      if isBodyInCylinder(b,center2,radius): # if a particle is inside 2nd cylinder ...
         mass2 += b.state.mass # ... increase mass2 by the particle mass
   return mass1,mass2 # return computed values
###

cheers
Jan

Thanks Jan for your help!

> I included the code you indicated, but now I don't know where I can check the mass of the particles or if it is necessary to include some other function to be able to visualize this value.

from yade import pack, plot

#Altura de caida del material
hc=0

#EQUIPO

#coordenadas centro
x=0
y=0
z=0

#coordenadas extension

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion
O.bodies.append(geom.facetBox((x,y,z),(f,a,(h+hc)),wallMask=31,wire=True))

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

      #O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Radio
radius=0.05

#Altura de la probeta
height=0.1

#Probeta 1
#Coordenadas centro
xp1=0.1
yp1=0
zp1=1

O.bodies.append(geom.facetCylinder((xp1,yp1,(-zp1-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0.933,0.090,.894), wire=False))

#Probeta 2
#Coordenadas centro
xp2=0.1
yp2=0
zp2=1

O.bodies.append(geom.facetCylinder((-xp2,yp2,(-zp2-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0,0,1), wire=False))

#coordenadas perimetro de particulas
#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

#porosidad
p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.00008
t2=0.00015
t3=0.00030
t4=0.00085
t5=0.00200
t6=0.00475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
angulo=radians(36) #radianes
#densidad del material
densidad=1418 #(kg/m3)
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
porosidad=0.9
#Algoritmo YADE para el material
snd=O.materials.append(FrictMat(young=young,frictionAngle=radians(angulo),density=densidad,
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
    Bo1_Facet_Aabb(),
    Bo1_Cylinder_Aabb(),
    Bo1_Wall_Aabb()
 ]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 #GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner()
]
O.dt=PWaveTimeStep()

def isBodyInCylinder(body,center,radius):

 d=body.state.pos-center
 d[2]=0
 distance=d.norm()
 return distance<radius

def masses():
 mass1=mass2=0.0
 center1=(xp1,yp1,0)
 center2=(-xp2,yp2,0)
 for b in O.bodies:
  if not isinstance(b.shape,Sphere): continue
  if isBodyInCylinder(b,center1,radius):
   mass1+=b.state.mass
  if isBodyInCylinder(b,center2,radius):
   mass2+=b.state.mass
  return mass1,mass2

> In the previous answer, you mentioned some supporting documents. You could attach them to me to consult them. They would be helpful to understand some functions that I need and for which there is not much information on the web.

thanks!!!

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

> but now I don't know where I can check the mass of the particles or if it is necessary to include some other function to be able to visualize this value.

run your simulation until a point when it is considered finished. Stop it. Call
masses()
in yade interactive terminal. The returned value are the masses in the two cylinders. You can print it or do whatever you want:
##
mass1,mass2 = masses()
print(mass1,mass2)
totalMass = mass1 + mass2
...
##

> In the previous answer, you mentioned some supporting documents. You could attach them to me to consult them. They would be helpful to understand some functions that I need

sure, just forgot to add it last time, see below

> for which there is not much information on the web.

I don't understand this..
"for b in O.bodies:" is surely documented as well as that b.state.mass is the mass of the body and b.state.pos its center.
Definitely you can easily find information on the web how to do a sum or how to check if a point is inside a cylinder..
Of course Yade documentation does not describe exactly your scenario.
BUT, Yade is open source and you are free (and welcome actually) to add such example or improve the documentation (or its part) if you find it not good! :-)

cheers
Jan

[1] https://yade-dem.org/doc/tutorial.html
[2] https://yade-dem.org/doc/user.html

Jan !! Help...

It seemed an easy function to use, but despite defining the function that stops the simulation in relation to the reduction of force, the simulation continues without showing me the results.

I included the following code in the simulation

####
O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
    Bo1_Facet_Aabb(),
    Bo1_Cylinder_Aabb(),
    Bo1_Wall_Aabb()
 ]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 #GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner(command="checkUnbalanced()"),
 PyRunner(command='isBodyInCylinder()'),
 PyRunner(command='masses()')
]
O.dt=PWaveTimeStep()

def checkUnbalanced():
 if unbalancedForce<0.05: O.pause()

def isBodyInCylinder(body,center,radius):

 d=body.state.pos-center
 d[2]=0
 distance=d.norm()
 return distance<radius

def masses():
 mass1=mass2=0.0
 center1=(xp1,yp1,0)
 center2=(-xp2,yp2,0)
 for b in O.bodies:
  print(mass1,mass2)
  if not isinstance(b.shape,Sphere): continue
  if isBodyInCylinder(b,center1,radius):
   mass1+=b.state.mass
  if isBodyInCylinder(b,center2,radius):
   mass2+=b.state.mass
  return mass1,mass2

###
According to what I understand, this function causes the simulation to stop when the equilibrium force of the particles is less than 0.05, but although the particles visually no longer interact, the simulation continues and in the script there is no result of mass.

Thanks for the documents, they were very necessary ... I hope that with the development of this project I can contribute something very interesting to the Yade community ...

It seems easy to use this software, but it has become an excellent professional challenge ...

thanks Jan!

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

> PyRunner(command="checkUnbalanced()"),
> PyRunner(command='isBodyInCylinder()'),
> PyRunner(command='masses()')

use exatly one of periods (iterPeriod [3], virtPeriod [4] or realPeriod [5]) for PyRunner, otherwise it is not run.

cheers
Jan

[3] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PyRunner.iterPeriod
[4] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PyRunner.virtPeriod
[5] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.PyRunner.realPeriod

Jan!
With the commands you indicated, I can now see that the simulation stops for the indicated condition.

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
    Bo1_Facet_Aabb(),
    Bo1_Cylinder_Aabb(),
    Bo1_Wall_Aabb()
 ]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 #GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner(command='checkUnbalanced()',realPeriod=1),
 PyRunner(command='isBodyInCylinder()',realPeriod=1),
 PyRunner(command='masses()',realPeriod=1)
]
O.dt=PWaveTimeStep()

def checkUnbalanced():
 if unbalancedForce()<0.05: O.pause()

def isBodyInCylinder(body,center,radius):

 d=body.state.pos-center
 d[2]=0
 distance=d.norm()
 return distance<radius

def masses():
 mass1=mass2=0.0
 center1=(xp1,yp1,zp1)
 center2=(-xp2,yp2,zp2)
 for b in O.bodies:

  if not isinstance(b.shape,Sphere): continue
  if isBodyInCylinder(b,center1,radius):
   mass1+=b.state.mass
  if isBodyInCylinder(b,center2,radius):
   mass2+=b.state.mass
  return mass1,mass2
  print(mass1,mass2)
>>>
1. But when executing the isBodyInCylinder function the script shows the following alert:

###
TypeError: isBodyInCylinder () takes exactly 3 arguments (0 give)
###

I understand that the arguments requested are the body (b), the center (center1 or center2) and the radius. I understand that these arguments are defined, so I don't know what the error might be.

>>>
2. If it is ok to add a PyRunner for each command created
###
PyRunner (command = 'checkUnbalanced ()', realPeriod = 1),
PyRunner (command = 'isBodyInCylinder ()', realPeriod = 1),
PyRunner (command = 'masses ()', realPeriod = 1)
###

>>>>>
3. Is there any reference value for the real Period or how can I ensure that the period indicated for the function is sufficient.

Thanks Jan!!

> TypeError: isBodyInCylinder () takes exactly 3 arguments (0 give)

Dear Estefany,
Please keep in mind point 5 in [1].
Although Jan may keep replying the questions because he is super-nice, it's now clearly not about computing total mass.
It may be beneficial for you to step back maybe, and study python functionalities and error messages a bit, before diving again into Yade.
It's very confusing for someone to not distinguish what's python and what's Yade, and I smell it could be your case. :)
Regards
Bruno

Side note: consider indenting with tabs or four spaces (at least), most text editors can do that for you. It's much more readable.

[1] https://www.yade-dem.org/wiki/Howtoask

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

Hello,

1.
as the error says, isBodyInCylinder () takes exactly 3 arguments, but you gave 0 arguments..
3 arguments = 3 values inside "calling" brackets.

solution:
isBodyInCylinder is just an auxiliary function for masses() function, just do not call it from PyRunner (delete that PyRunner from the O.engines)

2.
I guess yes for this case. If to use iterPeriod, virtPeriod or realPeriod depends always on specific context.

3.
see above, it is always different. If a specific value is sufficient, it is easy to test.

cheers
Jan

Jan !!! you cannot imagine how grateful I am with your help, it is very valuable and important to me.
If there is any way I can contribute to your help, please let me know.

With your contribution I corrected the code and this was the result:

from yade import pack, plot

#Altura de caida del material
hc=0

#EQUIPO

#coordenadas centro
x=0
y=0
z=0

#coordenadas extension

f=0.45
a=0.50
h=1

#estrutura del equipo de pluviacion
O.bodies.append(geom.facetBox((x,y,z),(f,a,(h+hc)),wallMask=31,wire=True))

#Mallas
#O.bodies.append(utils.wall((0,0,(.2+hc)),axis=2,color=(1,0.945,0.078),sense=-1))

#(0,0,(.2+hc)),axis=2,color(1,0.945,0.078),sense=0))

      #O.bodies.append(utils.wall((0,0,.2),axis=2,sense=0))

#PROBETA

#Radio
radius=0.05

#Altura de la probeta
height=0.1

#Probeta 1
#Coordenadas centro
xp1=0.1
yp1=0
zp1=1

O.bodies.append(geom.facetCylinder((xp1,yp1,(-zp1-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0.933,0.090,.894), wire=False))

#Probeta 2
#Coordenadas centro
xp2=0.1
yp2=0
zp2=1

O.bodies.append(geom.facetCylinder((-xp2,yp2,(-zp2-hc+(height/2))), radius=radius, height=height, fixed=True, wallMask=6, color=(0,0,1), wire=False))

#coordenadas perimetro de particulas
#esquina inferior
xi=0.225
yi=0.225
zi=0.2

#esquina superior
xs=0.225
ys=0.25
zs=0.5

#porosidad
p=0.903

#tamices
#los diametrso y los porcentajes acumulados siempre se deben ubicar de menor a mayor
#por cada t debe haber un p

#Diametro de la abertura de tamices psdSizes
t1=0.00008
t2=0.00015
t3=0.00030
t4=0.00085
t5=0.00200
t6=0.00475

#Porcentajes retenidos del material psdCumm
p1=0
p2=0.002
p3=0.111
p4=0.889
p5=0.9994
p6=1

#Material
#angulo de friccion
angulo=radians(36) #radianes
#densidad del material
densidad=1418 #(kg/m3)
#Modulo de Young
young=(1e8) #Pa
#Poisson
poisson=0.5
#contenidod de vacios
porosidad=0.9
#Algoritmo YADE para el material
snd=O.materials.append(FrictMat(young=young,frictionAngle=radians(angulo),density=densidad,
#label="SndG"
))

#Algoritmo YADE para el cuerpo de las esferas

sp=pack.SpherePack()

sp.makeCloud((xi,-yi,(zi+hc)), (-xs,ys,(zs+hc)), porosity=porosidad, psdSizes=[t1,t2,t3,t4,t5,t6], psdCumm=[p1,p2,p3,p4,p5,p6], distributeMass=False)

sp.toSimulation(material=snd,color=(0.882,0.866,0.839))

O.engines=[
 ForceResetter(),
 InsertionSortCollider([Bo1_Sphere_Aabb(),
    Bo1_Facet_Aabb(),
    Bo1_Cylinder_Aabb(),
    Bo1_Wall_Aabb()
 ]),
 InteractionLoop(
  # handle sphere+sphere and facet+sphere collisions
  [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
  [Ip2_FrictMat_FrictMat_FrictPhys()],
  [Law2_ScGeom_FrictPhys_CundallStrack()]
 ),
 #GravityEngine(gravity=(0,0,-9.81)),
 NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),
 PyRunner(command='checkUnbalanced()',realPeriod=1),
 PyRunner(command='masses()',realPeriod=1)
]
O.dt=PWaveTimeStep()

def checkUnbalanced():
 if unbalancedForce()<0.05: O.pause()

def isBodyInCylinder(body,center,radius):

 d=body.state.pos-center
 d[2]=0
 distance=d.norm()
 return distance<radius

def masses():
 mass1=mass2=0
 center1=(xp1,yp1,zp1)
 center2=(-xp2,yp2,zp2)
 for b in O.bodies:

  if not isinstance(b.shape,Sphere): continue
  if isBodyInCylinder(b,center1,radius):
   mass1+=b.state.mass

  if isBodyInCylinder(b,center2,radius):
   mass2+=b.state.mass
   return mass1, mass2
mass1,mass2=masses()
print"masa de cada cilindro (g):",(mass1, mass2)

totalMass=mass1+mass2
print"masa total de los cilindros (g):",(totalMass)

masa=sum([b.state.mass for b in O.bodies])
print"masa del material de simulacion (g):",(masa)

Now I only have a small question, the value that the script shows me before visualizing the simulation, corresponds to the value of the mass of the particles that remain inside each cylinder? is that when compared with the tax experiment they are much inferior.

Thanks again for your contributions!

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

> the value ... corresponds to the value of the mass of the particles that remain inside each cylinder?

yes

> is that when compared with the tax experiment they are much inferior.

yes, it is a simulation, not an experiment :-)

cheers
Jan

Chareyre (bruno-chareyre-9) said : #13

Is there O.run() somewhere, or you are speaking of mass at iteration 0?
B

Le lun. 24 févr. 2020 19:48, Jan Stránský <
<email address hidden>> a écrit :

> Question #688837 on Yade changed:
> https://answers.launchpad.net/yade/+question/688837
>
> Jan Stránský posted a new comment:
> > the value ... corresponds to the value of the mass of the particles
> that remain inside each cylinder?
>
> yes
>
> > is that when compared with the tax experiment they are much inferior.
>
> yes, it is a simulation, not an experiment :-)
>
> cheers
> Jan
>
> --
> You received this question notification because your team yade-users is
> an answer contact for Yade.
>
> _______________________________________________
> Mailing list: https://launchpad.net/~yade-users
> Post to : <email address hidden>
> Unsubscribe : https://launchpad.net/~yade-users
> More help : https://help.launchpad.net/ListHelp
>

Can you help with this problem?

Provide an answer of your own, or ask Estefany Carmona Alvarez for more information if necessary.

To post a message you must log in.