Extract normal force

Asked by NGANDU KALALA Gauthier on 2020-11-14

Hi,

I'm working on a tumbling ball mill simulation and I am interested on the extraction of the normal contact forces between the balls and the mill wall;
I use i.phys.normalForce.dot(i.geom.normal), the length of the simulation is 4098, for a step time of 2.44140e-4s and a total simulation time of 1 second. However by extracting the normal forces, I expected to have a force vector of the same length as the number of simulation, but I receive instead a force vector with a length less than the number of iterations ( i.e. 2362). I would like to know why I don't get the length of the vector equal to the length of the simulation?

How can i do to have the length of the vector normal force, with the same length of the total number of iteration?

Can someone help me?

Question information

Language:
English Edit question
Status:
Answered
For:
Yade Edit question
Assignee:
No assignee Edit question
Last query:
2020-11-24
Last reply:
2020-11-24

This question was reopened

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

Hello,

> I am interested on the extraction of the normal contact forces

What is "the extraction"? printing, saving, ...?
What is "the normal contact forces"? individual force vectors? magnitudes? sum of the vectors? sum of magnitudes? ... ?

> I use i.phys.normalForce.dot(i.geom.normal)

what is "i"?

> the length of the simulation is 4098

"length" is number of time steps?

> by extracting the normal forces

what is "extracting"?
what is "the normal forces"? vectors? magnitudes? ... ?

> I expected to have a force vector of the same length as the number of simulation

what is "the number of simulation"?
what is "a force vector"? I would expect "a force vector" to be a "standard Euclidean vector"..

> I would like to know why I don't get the length of the vector equal to the length of the simulation?
> How can i do to have the length of the vector normal force, with the same length of the total number of iteration?

What is "the vector"?
What is "the length of the vector"?
What is "the length of the simulation" and why would it equal to "the length of the vector"?
What is "the vector normal force"?

From the question, we can hardly even guess, please read [1] and provide more information.

> Can someone help me?

Yes, but first please help us. E.g. by providing a MWE [1], which would (should) have answered most of the "what is ..." questions..

cheers
Jan

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

Hi,

Thanks for your response, and you can find here after some response at your questions;

What is "the extraction"? printing, saving, ...?
The term extraction used here means , Saving force
What is "the normal contact forces"? individual force vectors? magnitudes? sum of the vectors? sum of magnitudes? ... ?
I call normal contact force, the magnitudes of the force of contact between ball and mill in the normal direction, saving at each step time; i need to save this force of contact at each time step. To know the evolution of this in time.
what is "the number of simulation"?
the number of simulation is the number of iteration here is 4098
what is "a force vector"? I would expect "a force vector" to be a "standard Euclidean vector"..
I call force vector, or normal force, in my question here, the vector contenning the value of magnitudes force in normal direction saving at each step time.

 the force at each step time. if the length of simulation, or the vector time length is 4098, i expect to have the saving vector force with the same length.

Briefly i want to save the evolution of interaction force in time; reason why i insert a loop for j in O.time in my script; But i can't have the result i need.
I don't know if there is an other way, to save the contact force at each step time, and to plot his evolution in time;

I' m trying to do it in my mwe here after:

# Material
rho = 7640.0 # [kg/m^3] density (of the balls)
kn = 1.0e5 # [N/m] equivalent normal stiffness
kt = kn # [N/m] equivalent tangential stiffness
en = 0.3 # [-] normal coefficient of restitution (sphere/sphere)
mu = 0.75 # [-] friction coefficient
cn = 0.75
ct = cn

# Geometry/Filling
length = 0.4 # [m] axial length of the mill
rMean = 0.772/2.0 # [m] liner profile mean radius
rBall = 0.01 # [m] ball radius
fRatio = 1.0 # [%] filling ratio
sphRadFuzz = 0.8

# Simulation

tEnd = 1.0 # [s] total simulation time
O.dt = 2.44140e-4 # [s] fixed time step
tRec = 1.0e-2 # [s] recording interval

#material of cylinder
cylMat = O.materials.append(ViscElMat(kn=2*kn,ks=2*kt,cn=2*cn,
                                             cs=2*ct,frictionAngle=.1))
#material of spheres
ballMat = O.materials.append(ViscElMat(kn=2*kn,ks=2*kt,cn=2*cn,
                                             cs=2*ct,frictionAngle=.1,density=rho))

ids=O.bodies.append(geom.facetCylinder((0,0,0),0.8,0.4,segmentsNumber=64, wallMask=7, angleRange=None,
                   radiusTopInner=-1, radiusBottomInner=-1,material=cylMat))

# define domains for initial cloud of red and blue spheres
packHt=.8*rMean # size of the area
bboxes=[(Vector3(-.5*length,-.5*packHt,-.5*packHt),Vector3(.5*length,0,.5*packHt)),(Vector3(-.5*length,-.5,-.5*packHt),Vector3(.5*length,.5*packHt,.5*packHt))]
colors=(1,0,0),(0,0,1)
for i in (0,1): # red and blue spheres
 sp=pack.SpherePack();
 bb=bboxes[i];
 vol=(bb[1][0]-bb[0][0])*(bb[1][1]-bb[0][1])*(bb[1][2]-bb[0][2])
 sp.makeCloud(bb[0],bb[1],rBall,sphRadFuzz,int(.25*vol/((4./3)*pi*rBall**3)),False)
 O.bodies.append([utils.sphere(s[0],s[1],material=ballMat,color=colors[i]) for s in sp])

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
    InteractionLoop(

        [Ig2_Sphere_Sphere_ScGeom(),
         Ig2_Facet_Sphere_ScGeom(),
         Ig2_Wall_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys()],
        [Law2_ScGeom_MindlinPhys_Mindlin()]
    ),

   #GravityEngine(gravity=(0,-3e4,0)), # gravity artificially high, to make it faster going ;-)
   RotationEngine(ids=ids,angularVelocity=1000,rotateAroundZero=True,rotationAxis=(0,0,1)),
   NewtonIntegrator(damping=0,gravity=[0,0,-9.81]),
   PyRunner(command='exportForces()',iterPeriod=1)
]
def exportForces():

  # get data

for j in O.time:
  for i in O.interactions:
    if isinstance(O.bodies[i.id1].shape,Sphere) \
        and isinstance(O.bodies[i.id2].shape,Facet):

  totalNormalForce = i.phys.normalForce

  # Save data

  out=open(os.getcwd()+'/output/NormalForceMagnitude','a')
  out.write('%s\n'%str(NormalForceMagnitude))
  out.close()
###
O.run(int(math.ceil(tEnd/O.dt))+1)

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

Hello,

> I call normal contact force, the magnitudes of the force of contact between ball and mill in the normal direction
> > the force at each step time
> i want to save the evolution of interaction force in time

still this is very confusing.
In the simulation, you have (possibly) many balls and many ball-mill interactions and many forces.
What is "(the) force"?
- force of one specific interaction?
- sum of force vectors?
- sum of force magnitudes?
- ... ?

> I call force vector, or normal force, in my question here, the vector contenning the value of magnitudes force in normal direction saving at each step time.

I strongly recommend to use a different name for the "vector" in this case not to be confused with "ordinary" force vector.
Maybe array? field?

> i insert a loop for j in O.time in my script; But i can't have the result i need.

it results in an error, no surprise you do not have the result.
Fix errors first.

> I don't know if there is an other way, to save the contact force at each step time

There are many ways.
To choose a suitable one, first of all you need to know (or explain it properly in the case you want help from others) what you want. Which does not seem to be the case (see above).

> ... mwe ...

W = working, which is not the case here. Please, fix the indentation to be valid Python code.
Moreover, O.time is a number and you cannot iterate it "for j in O.time"

> totalNormalForce = i.phys.normalForce

**guessing** it should be inside the sphere-facet condition:
- why the equal sign? Shouldn't the value of totalNormalForce be incremented by i.phys.normalForce instead of rewriting?
- you write you want normal force magnitudes. here you assign/sum normal force vector. What exactly do you want?
- you can save the data "manually" as you do here. But you should clear the output file first (the code appends data to whatever already is in the file)
- instead of "manual" saving, Yade "standardly" uses plot module for this purposes [1]

cheers
Jan

[1] https://yade-dem.org/doc/user.html#tracking-variables

Hi, Thank you, This some response,

still this is very confusing.
In the simulation, you have (possibly) many balls and many ball-mill interactions and many forces.
What is "(the) force"?
- force of one specific interaction?
- sum of force vectors?
- sum of force magnitudes?

Exactly you are right. In the simulation, several interactions take place between balls and mill, because of the presence of several balls. That's why I want to save at each iteration (or step time), the sum of forces magnitudes.

I'will provide my mwe soon

This is my MWE;

As I said before, I want to save for each iteration, the sum of force (amplitude) of interactions between balls and mill;
In this MWE, I have rather calculated the forces when there is interaction; But I wish to calculate it rather for each iteration; this will allow me to represent the evolution of its forces in time;

from yade import pack,utils, plot,ymport,qt
#######################################################################################

fr = 0.1;rho=7640.0;length = 4.0;rMean = 0.772/2.0;sphRadFuzz = 0.8;rBall = 0.04 ;

MillMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=2*fr))
ballMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=fr,density=rho))

#######################################################################################

# Create Mill from facets
Mill=O.bodies.append(geom.facetCylinder((0.0,0.0,0.5),8.0,height=4.0,segmentsNumber=64,wallMask=7,material=MillMat))
#Walls=O.bodies.append([utils.wall(position=(0,0,0),axis=2,sense=0,material=MillMat),
# utils.wall(position=(0,length,0),axis=2,sense=0,material=MillMat)
# ])
# define domains for initial cloud of red and blue spheres
packHt=.8*rMean # size of the area
bboxes=[(Vector3(-.5*length,-.5*packHt,-.5*packHt),Vector3(.5*length,0,.5*packHt)),(Vector3(-.5*length,-.5,-.5*packHt),Vector3(.5*length,.5*packHt,.5*packHt))]
colors=(1,0,0),(0,0,1)
for i in (0,1):
 sp=pack.SpherePack();
 bb=bboxes[i];
 vol=(bb[1][0]-bb[0][0])*(bb[1][1]-bb[0][1])*(bb[1][2]-bb[0][2])
 sp.makeCloud(bb[0],bb[1],rBall,sphRadFuzz,int(.25*vol/((4./3)*pi*rBall**3)),False)
 O.bodies.append([utils.sphere(s[0],s[1],material=ballMat,color=colors[i]) for s in sp])

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(

   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),

    NewtonIntegrator(gravity=(0,-9.81,0),damping=0.0),

    RotationEngine(ids=Mill,rotationAxis=[0,0,-1],rotateAroundZero=True,
                 angularVelocity=1),

   ]

O.dt=.5*PWaveTimeStep()

def saveForces():

   for i in O.interactions:

    if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
       MagnitudeNormalForce = sum((i.phys.normalForce.norm()))
    else :
       MagnitudeNormalForce = 0

   # out=open(os.getcwd()+'/Resultat/MagnitudenormalForce','a')
   # out.write('%s\n'%str(MagnitudenormalForce))
   # out.close()

O.run(100000)
#O.save('Mymodel2.yade')

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

You defined saveForces, but it is not used anywhere.
You have to use saveForces somehow, most likely form a PyRunner (as you already had in the previous codes).
cheers
Jan

Hi Jan,

This is my updated MWE, as i tell you before, i need the simulation show the plot on the screen, and update while the simulation runs in time, and save the result in txt format.

rom yade import pack,utils, plot,ymport,qt
#######################################################################################

fr = 0.1;rho=7640.0;length = 4.0;rMean = 0.772/2.0;sphRadFuzz = 0.8;rBall = 0.04 ;

MillMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=2*fr))
ballMat = O.materials.append(ViscElMat(kn=2*1.0e5,ks=2*1.0e5,frictionAngle=fr,density=rho))

#######################################################################################

# Create Mill from facets
Mill=O.bodies.append(geom.facetCylinder((0.0,0.0,0.5),8.0,height=4.0,segmentsNumber=64,wallMask=7,material=MillMat))
#Walls=O.bodies.append([utils.wall(position=(0,0,0),axis=2,sense=0,material=MillMat),
# utils.wall(position=(0,length,0),axis=2,sense=0,material=MillMat)
# ])
# define domains for initial cloud of red and blue spheres
packHt=.8*rMean # size of the area
bboxes=[(Vector3(-.5*length,-.5*packHt,-.5*packHt),Vector3(.5*length,0,.5*packHt)),(Vector3(-.5*length,-.5,-.5*packHt),Vector3(.5*length,.5*packHt,.5*packHt))]
colors=(1,0,0),(0,0,1)
for i in (0,1):
 sp=pack.SpherePack();
 bb=bboxes[i];
 vol=(bb[1][0]-bb[0][0])*(bb[1][1]-bb[0][1])*(bb[1][2]-bb[0][2])
 sp.makeCloud(bb[0],bb[1],rBall,sphRadFuzz,int(.25*vol/((4./3)*pi*rBall**3)),False)
 O.bodies.append([utils.sphere(s[0],s[1],material=ballMat,color=colors[i]) for s in sp])

O.engines=[
   ForceResetter(),
   InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb(),Bo1_Wall_Aabb()]),
   InteractionLoop(

   [Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom(),Ig2_Wall_Sphere_ScGeom()],
      [Ip2_FrictMat_FrictMat_FrictPhys()],
      [Law2_ScGeom_FrictPhys_CundallStrack()]
   ),

    NewtonIntegrator(gravity=(0,-9.81,0),damping=0.0),

    RotationEngine(ids=Mill,rotationAxis=[0,0,-1],rotateAroundZero=True,
                 angularVelocity=1),
    PyRunner(command='addPlotData()',iterPeriod=100),
   ]

O.dt=.5*PWaveTimeStep()

# collect history of data which will be plotted
def addPlotData():

  for i in O.interactions:

    if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
       Force = sum((O.forces.f(i).norm()))
    else :
       Force = 0

  plot.addData(i=O.iter,Force=sum(O.forces.f(i).norm()))

plot.plots={'i':('Force')}

# show the plot on the screen, and update while the simulation runs
plot.plot()

O.run(100000)
#O.save('Mymodel2.yade')

but when i run the simulation i get this error;

ArgumentError: Python argument types in
    ForceContainer.f(ForceContainer, Interaction)
did not match C++ signature:
    f(yade::pyForceContainer {lvalue}, long id, bool sync=False)
In [1]: Erreur de segmentation (core dumped)

Can i have some correction in my code please?

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

code 1:
> def exportForces():
> totalNormalForce = i.phys.normalForce

code 2:
> def saveForces():
> MagnitudeNormalForce = sum((i.phys.normalForce.norm()))

code 3:
> def addPlotData():
> Force = sum((O.forces.f(i).norm()))

Every update, you completely changed the namings. Yes, it is just naming and here it is clear, but next time please try to stay consistent.
There is no reason for such renaming and then I would prefer not to rename the stuff to prevent confusion.

Every update, you changed the force computation method, without any comment. Next time please consider giving some comment about the reasons.
Ok, from the discussion it came out that the magnitude should be force.norm() instead of force vector itself, so i.phys.normalForce -> i.phys.normalForce.norm() is reasonable.
But why the sum (inside for loop)?
And why i.phys.normalForce -> O.forces.f ?(!?!?)

> This is my updated MWE
> rom yade import pack,utils, plot,ymport,qt

W=working, starting a script with "rom yade..." surely is not :-)

> for i in O.interactions:
> ... O.forces.f(i)

O.forces.f needs body ID, see documentation [1]. In your code you pass "i", which is Interaction... (this is what the error says)

O.forces.f returns resulting force on the body. It has no relation to individual interaction forces.
(O.forces.f is sum of all interaction forces on the body, but with no backwards influence)

> sum((O.forces.f(i).norm()))

force.norm() is one number.
What you expect from sum(number)?
try it in python. e.g. sum(2)

> for i in O.interactions:
> if ...:
> Force = sum((O.forces.f(i).norm()))
> else :
> Force = 0

this way, you loop over all interactions just to assign Force to be the force from the very last interaction or 0...

A possible solution:
###
def addPlotData():
    Force = 0.0 # initial zero value
    for i in O.interactions:
        if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet):
            Force += i.phys.normalForce.norm() # increment the Force
    # or even a one liner
    # Force = sum(i.phys.normalForce.norm() for i in O.interactions if isinstance(O.bodies[i.id1].shape,Facet) or isinstance(O.bodies[i.id2].shape,Facet))
    plot.addData(i=O.iter,Force=Force)
###

You can save the data with plot.saveDataTxt [2], either inside addPlotData, or in some other PyRunner with lower frequency, or after the simulation is finished, or ... whenever you need.

cheers
Jan

[1] https://yade-dem.org/doc/yade.wrapper.html#yade.wrapper.ForceContainer.f
[2] https://yade-dem.org/doc/yade.plot.html#yade.plot.saveDataTxt

Can you help with this problem?

Provide an answer of your own, or ask NGANDU KALALA Gauthier for more information if necessary.

To post a message you must log in.