# yade-users team mailing list archive

## Re: Elastic energy

• From: Janek Kozicki <janek_listy@xxxxx>
• Date: Mon, 5 Jul 2010 14:44:29 +0200

```Bruno Chareyre said:     (by the date of Mon, 05 Jul 2010 13:35:42 +0200)

> FYI, energy is conservated almost exactly in a periodic triaxial test
> (taking into account also the work input from boundaries). I just checked.
> It seems errors come only in dynamic bouncing spheres, I hope you'll
> find the problem folk!

dynamic you say...

So - what exactly happens when the contact breaks? The plastic energy
is changing into kinetic energy??? Look at the elastic energy, it's
almost zero, how is it possible that kinetic energy is increasing?

I made some change to plots, see attachment.

Slip_how is what we are expected to calculate.

But Ek_Eel is more interesting, you can see here, that there is no
way that elastic energy could contribute to kinetic energy *that much*.

How plastic energy can turn back into kinetic energy??? Because where
else that increase of kinetic energy would be coming from?

To check if current calculation of plasticDissipation is only correct
when it *increases*, but skips the case when it decreases, I made
another plot - the negative_Ek_only_decreasing. It is calculated from
remaining global energy (kinetic and elastic), but only when that
global energy *decreases*, thus the case when plasticDissipation
increases.

Look at how strikingly similar negative_Ek_only_decreasing is to
plasticDissipation ! And they are both calculated in a totally
different way.

There is a weird factor of *2.0 included. Why, again? They look so
similar, that I don't think it is a coincidence. I guess that the
only difference between two curves is coming from the same thing that
caused me to add the lines:

if(O.materials.frictionAngle == 0):
return 0

to the python code. Comment that out and you will see an integral of
energy instabilities (summing only when Ek+Eel decreases).

To conclude, I suppose that this formula:

plasticDissipation +=
((1/currentContactPhysics->ks)*(trialForce-shearForce))//plastic disp.
.dot(shearForce);//active force

is correct only when plasticDissipation increases. But totally
ignores the case when it decreases.

--
Janek Kozicki                               http://janek.kozicki.pl/  |
```
```#!/home/janek/bin/yy --threads=1
# -*- 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',traceEnergy=True)]
),
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 *

prevE=28571.0
sum_negative_E=0.0

def only_decreasing():
if(O.materials.frictionAngle == 0):
return 0
global prevE
global sum_negative_E
if(prevE>utils.kineticEnergy()+dry.elasticEnergy()):
sum_negative_E += (prevE - utils.kineticEnergy() - dry.elasticEnergy())*2.0
prevE=utils.kineticEnergy()+dry.elasticEnergy()
return sum_negative_E

plot.plots={'t':('Ek','Eel','Ek_Eel','Slip','Etot','Slip_how','negative_Ek_only_decreasing'),'t_':('Slip')}
def monitoring():
Etot=utils.kineticEnergy()+dry.elasticEnergy()+dry.plasticDissipation(),Ek_Eel=utils.kineticEnergy()+dry.elasticEnergy(),negative_Ek_only_decreasing=only_decreasing())
O.run(10000,True);
print "\n**** Now friction angle is set from 0 to 25 degrees. ****\n";
O.run(30000,True); plot.plot()

```

Attachment: scr_3224.png
Description: PNG image