Energy tracking functions

Asked by Stefano Simone on 2020-04-23

Good morning,
I'm new to YADE and I'm using it for a part of my master thesis.

I want to model a granular damper. In the case study the phases are the following
1) The spheres are inserted in a parallelepiped box (approx. 2cmx1cmx0.5cm)
2) The vibration of the all system makes the box shear cyclically.

(code at the end)
IMG to schematize the system --->

The idea is to build before all the objects (box as facets and sphere pack), then to compress the pack of spheres by imposing a motion on the cover and finally to excite the system with HarmonicRotationEngine.

I use Hertz-Mindlin contact law (code below).

I'm facing some "physical" problems:
1) I want to get the total energy given to the system, so I save the Torque aruond z and the rotation to integrate over an half cycle. Is there another way to get this total energy directly from YADE? Am I taking the wrong data or interpreting them bad?

2) From Law2_ScGeom_MindlinPhys_Mindlin I take the terms of dissipation but what I noticed is the loss rate (dissipated energy per cycle) is higher than the stored energy per cycle estimated before. I think there must be some problem with how I compute the stored energy or within the simulation.

3) I want to relate the confinement pressure (F_closing in code) with the dissipation in this system, so I expect that at low pressure friction is not effective and at very high pressure is not effective as well, with a sort of optimum in the middle. What I see is that I can increase the pressure on the pack more and more (up to complete compenetration of the sphere layers) and the dissipation by friction only increases to not physical values? What is the reason of this? Am I missing something in controlling?

4) Are there other ways to achieve the my result? If yes, are there some examples to referr to?

Thank you very much for the helping and for the care you all have in this community. I really hope that you can help me with this problems that make me stuck. I hope as well that I did respect all the rules for making questions on this forum.

Thank you again and regards.

Stefano Simone

The code is the following:

from yade import plot,utils
from yade import pack
import math
import random
import numpy


f = 250
period = 1/f
wrad = 2*pi*f
ang_ampl = 0.04
e_n = 0.63
e_s = 0.63
frictionAngle = 0.197
num_per = 1
g = 9.81
F_closing= 15000

#O.materials.append(FrictMat(young = 210e9, poisson = 0.35, frictionAngle = 0.7, density = 7800, label = 'spmat'))
O.materials.append(FrictMat(young = 70e9, poisson = 0.22, frictionAngle = 0.197, density = 2500 , label = 'spmat')) #sphere of glass
O.materials.append(FrictMat(young = 210e9, poisson = 0.35, frictionAngle = 0, density = 0, label = 'wmat'))


    # sphere packing creation
pred = pack.inAlignedBox((-0.010,-0.005,0), (.010,0.005,0.005))
O.bodies.append(pack.regularHexa(pred, radius=0.00075, gap=0.00001, material = 'spmat'))

    # random distribution of radii
Dev = 0.0001
media = 0.00075
val_max = media + Dev
val_min = media - Dev

radii = numpy.linspace(val_min, val_max, num=50, endpoint=True, retstep=False)

for b in O.bodies:
    if isinstance(b.shape,Sphere):
        b.shape.radius = random.choice(radii)

    # total volume and number of the spheres
V_spheres = 0
number = 0

for b in O.bodies:
    if isinstance(b.shape, Sphere):
        V_spheres = V_spheres + 4/3*pi*b.shape.radius**3
        number = number + 1

    # Container

Down = O.bodies.append(geom.facetBox((0,0,0),(0.02,0.005,0)))

Right = O.bodies.append(geom.facetBox((0.01,-0.005,0.0025), (0,0.005*2,0.0025)))

Left = O.bodies.append(geom.facetBox((-0.01,-0.005,0.0025), (0,0.005*2,0.0025)))

Bottom = O.bodies.append(geom.facetBox((0,-0.0050,0.005), (0.01,0,0.005)))

Upper = O.bodies.append(geom.facetBox((0,0.005,0.005), (0.014,0,0.005)))

Cover = O.bodies.append(geom.facetBox((0,0,0.005),(0.02,0.005,0)))


O.engines = [
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
        [Ig2_Sphere_Sphere_ScGeom(), Ig2_Facet_Sphere_ScGeom()],
        [Ip2_FrictMat_FrictMat_MindlinPhys(frictAngle = 0.05, en=e_n , es=e_s, label='Ip2')],
        [Law2_ScGeom_MindlinPhys_Mindlin(includeAdhesion=False, calcEnergy=True, label='Mindlin')]
    NewtonIntegrator(gravity=(0, 0, -g), damping=0.1),
    TranslationEngine(dead = True, translationAxis=[0, 0, 1], velocity=-0.002, ids=Cover, label='pres'),
    HarmonicRotationEngine(dead=True, f =f, A = ang_ampl*pi, fi = 0*pi/2 , ids=Right, rotateAroundZero=True, rotationAxis=(0, 0, 1), zeroPoint=(0.01, -0.005, 0), label = 'rotR'),
    HarmonicRotationEngine(dead=True, f =f, A = ang_ampl*pi, fi = 0*pi/2 , ids=Left, rotateAroundZero=True, rotationAxis=(0, 0, 1), zeroPoint=(-0.01, -0.0050, 0), label = 'rotL'),
    #HarmonicMotionEngine(dead=True, f =[f,0,0], fi = [-pi/2*0, 0, 0],A = [sin(ang_ampl*pi),0,0] , ids=Upper, label = 'rotU'),
    PyRunner(dead = False, command = 'startPressure()', iterPeriod = 100, label = 'phase'),
    PyRunner(dead = False, command = 'plotData()', iterPeriod = 1000, label = 'data')

O.dt = PWaveTimeStep()
O.trackEnergy = True

def startPressure():
        if O.iter > 250000 and unbalancedForce() < 0.05:
            pres.dead = False
            pres.velocity = -0.05
            phase.command = 'stopPressure()'
            phase.iterPeriod = 100

def stopPressure():
    top_sphere = max([b.state.pos[2] for b in O.bodies if isinstance(b.shape, Sphere)])
    if (O.bodies[Cover[0]].state.pos[2] - top_sphere)< 0.0005:
        pres.velocity = -0.05
    if abs((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(1))*1 >= F_closing:
    #if O.bodies[Cover[0]].state.pos[2] <=0.5:
        Ip2.frictAngle = frictionAngle
        O.materials[0].frictionAngle = frictionAngle
        O.materials[1].frictionAngle = frictionAngle
        Ip2.en = 0.63 = 0.63
        print(((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(1))*1)
        pres.dead = True
        O.bodies[Cover[0]].state.blockedDOFs = 'xyzXYZ'
        O.bodies[Cover[0]].state.vel = (0, 0, 0)
        O.bodies[Cover[1]].state.blockedDOFs = 'xyzXYZ'
        O.bodies[Cover[1]].state.vel = (0,0,0)
        phase.command = 'Rotation_control()'

def Rotation_control():
    if unbalancedForce() < 0.02 and O.time > 0.4 :
        rotR.dead = False
        rotR.zeroPoint[0] = O.bodies[Right[0]].state.pos[0]
        rotL.dead = False
        #rotU.dead = False
        data.dead = False

def plotData():
    t = O.time,
    Edf = Mindlin.frictionDissipation,
    Eds = Mindlin.shearDampDissip,
    Edn = Mindlin.normDampDissip,
    Etot = Mindlin.normElastEnergy() + Mindlin.shearEnergy,
    Ed = Mindlin.frictionDissipation + Mindlin.shearDampDissip + Mindlin.normDampDissip,
    v_av = numpy.average([b.state.vel.norm() for b in O.bodies if isinstance(b.shape, Sphere)]),
    TorqueR = O.forces.t(O.bodies[Right[0]].id)[2] + O.forces.t(O.bodies[Right[1]].id)[2],
    TorqueL = O.forces.t(O.bodies[Left[0]].id)[2] + O.forces.t(O.bodies[Left[1]].id)[2],
    ForceUP = -O.forces.f(O.bodies[Cover[0]].id)[2] - O.forces.f(O.bodies[Cover[1]].id)[2],
    v_UP = O.bodies[Cover[0]].state.vel[2],
    Angle = O.bodies[Right[0]].state.ori[2],
    AngVel = O.bodies[Right[0]].state.angVel[2],
    Pressure = ((O.forces.f(O.bodies[Cover[0]].id)[2] + O.forces.f(O.bodies[Cover[1]].id)[2])/(0.0002))*1e-6,
    Unbalanced = unbalancedForce(),
   # if O.iter > 1000000:
   # plot.saveDataTxt('Side.dat', vars=('t', 'Edf', 'Eds', 'Edn', 'TorqueL','TorqueR','ForceUP', 'v_UP','AngVel','Angle','Unbalanced','Pressure'))

plot.plots = {'t':('Ed','kinetic','Etot'), 't ':'TorqueR', 't ':'Angle', 't ': 'TorqueL', 't ':'v_UP'}


Question information

English Edit question
Yade Edit question
No assignee Edit question
Solved by:
Jérôme Duriez
Last query:
Last reply:
Jérôme Duriez (jduriez) said : #1


> I hope as well that I did respect all the rules for making questions on this forum.

Not really ;-) See maybe In short

- please one Launchpad question for one question

- you would increase your chances getting help / solving yourself your problems minimalizing the shown script. Commented commands show this was not done.. (I have to regret you're not the first one to do so, and won't be the last..)

- generally speaking, you're also more prone to getting answers about meanings / usages / existence of specific YADE functions. Less probable for general discussion about a research workflow for a given project

As for 1), there is for instance automatic energy tracking in YADE through (which you actually have in your script ?). Doc is for instance at, with several usage examples elsewhere, eg at

Stefano Simone (stefanosimone) said : #2

Thank you for the answer!

I'm sorry, next time I'll follow the rules. Thank you also for the useful tips for asking issue.

In this theard I will limitate my questions to the energy tracking.

1) I know that in the terms of ** there is kinetic energy, gravity Work and nonviscous damping. If I save the does it compute the sum of these three?

2) Considering that in my system the mass of spheres is very very small (~10^-6 kg) it should result Etot =~ nonviscDamp. What is the meaning of this term?

3) Is it correct to compare this directy with friction dissipation term from HM Law2 (Law2_ScGeom_MindlinPhys_Mindlin.frictionDissipation) or the other computed dissipated energies?

Thank you in advance.


Best Jérôme Duriez (jduriez) said : #3

1) The name of is pretty well chosen since it sums everything what's inside What's exactly inside depends on your simulation, but can be easily checked eg with

2) nonviscDamp is the work of (non-physical) Cundall's damping forces, added by NewtonIntegrator.damping. As these forces are non-physical, it does not have much of a meaning, except telling how much one is cheating for his/her DEM simulations.

3) On the other side, I guess Law2_ScGeom_MindlinPhys_Mindlin.frictionDissipation is the plastic work of tangential forces at contact (at least, that's what is Law2_ScGeom_FrictPhys_CundallStrack.plasticDissipation. Doc of the other one might help).
Friction-limited tangential forces do exist in nature, and this one has a sound meaning.

PS: a more meaningful title to the thread might help (updating question titles is welcome, updating initial questions once someone answered never is, for your information), now we fortunately got a little more specific ;-)

Stefano Simone (stefanosimone) said : #4

Ok, thank you for the answer and the patience. I updated the title.

Stefano Simone (stefanosimone) said : #5

Thanks Jérôme Duriez, that solved my question.