yade-users team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

Hello all,

sorry to bother you again about energy, but I recently had some problems
with that and so now I had to go back to the simplest problem.

I am attaching the script, just run it to see the plot (attached as well,
anyway). The test could not be simpler: elastic impact in the normal
direction. No mechanical energy is present, only elastic potential energy
and kinetic energy. I am simply using the linear spring model. I plot the
results from functions utils.kineticEnergy() and elasticEnergy() of
Law2_ScGeom_FrictPhys_Basic contact law. It seems that elastic energy is
half of what it should be. Could you help me to see where I am wrong?

Thank you.

Chiara

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

# Test of elastic potential energy (normal direction only hereby)
from yade import utils
import math
#---------------------------------
# geometry
r1,r2=1e-2,1e-2 # radii 
p1,p2=[0,0,0],[r1+r2,0,0] # center positions 
#---------------------------------
# material
young=600.0e6
poisson=0.6 
density=2.60e3 
frictionAngle=radians(0) 
O.materials.append(FrictMat(young=young,poisson=poisson,density=density,frictionAngle=frictionAngle))
O.bodies.append(utils.sphere(p1,r1,dynamic=False,wire=True)) # body id=0 
O.bodies.append(utils.sphere(p2,r2,dynamic=True,wire=True)) # body id=1
O.bodies[1].state.vel=Vector3(-0.005,0,0) # set 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='contact')]
	),
	NewtonIntegrator(damping=0.0),
	PeriodicPythonRunner(iterPeriod=1,command='energy()')
]
#---------------------------------
# time step
O.dt=.01*utils.PWaveTimeStep()
O.saveTmp('init')
#---------------------------------
from yade import qt
qt.View()
qt.Controller()
#---------------------------------
# plot some results
from math import *
from yade import plot
plot.plots={'t':('Eel','Ek','Etot')}
def energy():
	plot.addData(Eel=contact.elasticEnergy(), # potential elastic energy
		     Ek=utils.kineticEnergy(), # kinetic energy
		     Etot=contact.elasticEnergy()+utils.kineticEnergy(), # total energy
		     t=O.time)
O.run(1000,True); plot.plot()

Attachment: Plot.png
Description: PNG image