Energy tracking functions

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 ---> https://we.tl/t-0e9ZcKzmhp

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

#parameters

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

#materials
#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'))

#bodies

    # 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)))

#engines

O.engines = [
    ForceResetter(),
    InsertionSortCollider([Bo1_Sphere_Aabb(), Bo1_Facet_Aabb()]),
    InteractionLoop(
        [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
        Ip2.es = 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():
    plot.addData(
    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(),
    **O.energy
    )
   # 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'}
plot.plot()

O.saveTmp()