# yade-users team mailing list archive

## Re: Elastic energy

```Hi guys,

I am experiencing again a bit of troubles with total energy conservation.
For a sphere-sphere impact test everything seems fine. Now I have tried
another simple example: few spheres with an input initial velocity in a box.
I attach the script, if you run it you will obtain the plot in the attached
figure. If there is no friction, then total energy is more or less constant
as shown in the first part of the plot (not exactly constant due to the way
However, as long as I include some friction, both total and kinetic energy
decreases and no slip occurs since the function plasticDissipation() of
Law2_ScGeom_FrictPhys_Basic contact law returns 0. I have already checked on
paper and discussed with Bruno and the code which computes
plasticDissipation in ScGeom is formally correct.

So, where is the energy going if no dissipation occurs?

Thanks for any help or suggestions.

Chiara

_______________________________________________
>>
>
>
```
```#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

# Script to monitor the energy of a system
#-----------------------------------------------------
# material
#-----------------------------------------------------
young=600.0e6
poisson=0.6
density=2.60e3
O.materials.append(FrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle,label='mat'))
#-----------------------------------------------------
# geometry
#-----------------------------------------------------
# create a random packing in a box space
sp=pack.SpherePack()
mn=Vector3(0,0,0)
mx=Vector3(1.0,1.0,1.0)
ntot=10
R=0.1
std=0.0
sp.makeCloud(minCorner=mn,maxCorner=mx,rMean=R,rRelFuzz=std,num=ntot,periodic=False,porosity=0.85)
for s in sp:
O.bodies.append(utils.sphere(s,s))
# create some boundaries
wires=True
O.bodies.append(utils.box(center=[-0.1,0.5,0.5],extents=[.0025,1.0,1.0],dynamic=False,wire=wires,material='mat'))
O.bodies.append(utils.box(center=[1.1,0.5,0.5],extents=[.0025,1.0,1.0],dynamic=False,wire=wires,material='mat'))
O.bodies.append(utils.box(center=[1.1/2,-0.1,0.5],extents=[1.0,0.0025,1.0],dynamic=False,wire=wires,material='mat'))
O.bodies.append(utils.box(center=[1.1/2,1.1,0.5],extents=[1.0,0.0025,1.0],dynamic=False,wire=wires,material='mat'))
O.bodies.append(utils.box(center=[0.5,0.5,-0.1],extents=[1.0,1.0,0.0025],dynamic=False,wire=wires,material='mat'))
O.bodies.append(utils.box(center=[0.5,0.5,1.1],extents=[1.0,1.0,0.0025],dynamic=False,wire=wires,material='mat'))
#-----------------------------------------------------
# initial condition (velocities)
#-----------------------------------------------------
v=7.0
for b in O.bodies:
if b.id%2 == 0:
b.state.vel=Vector3(-v,v,-v) # assign an initial velocity
else:
b.state.vel=Vector3(v,-v,v) # assign an initial velocity
#-----------------------------------------------------
# list of engines
#-----------------------------------------------------
O.engines=[
ForceResetter(),
BoundDispatcher([Bo1_Sphere_Aabb(),Bo1_Box_Aabb()]),
InsertionSortCollider(),
InteractionDispatchers(
[Ig2_Sphere_Sphere_ScGeom(),Ig2_Box_Sphere_ScGeom()],
[Ip2_FrictMat_FrictMat_FrictPhys()],
[Law2_ScGeom_FrictPhys_Basic(label='dry')]
),
NewtonIntegrator(damping=0.0), # *** NO DAMPING ***
PeriodicPythonRunner(iterPeriod=10,command='monitoring()')
]
#-----------------------------------------------------
# time step
#-----------------------------------------------------
O.dt=.1*utils.PWaveTimeStep()
O.saveTmp('init')
qt.View()
qt.Controller()
#-----------------------------------------------------
# plot some results
#-----------------------------------------------------
from math import *
plot.plots={'t':('Ek','Eel','Slip','Etot'),'t_':('Slip')}
def monitoring():
Etot=utils.kineticEnergy()+dry.elasticEnergy()+dry.plasticDissipation())
O.run(10000,True);
print "\n**** Now friction angle is set from 0 to 25 degrees. ****\n";
O.run(30000,True); plot.plot()

```

Attachment: Energy.png
Description: PNG image