yade-users team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

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
we compute velocities, as already discussed in this thread
http://www.mail-archive.com/yade-users@xxxxxxxxxxxxxxxxxxx/msg01729.html).
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

_______________________________________________
>> Mailing list: https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
>> Post to     : yade-users@xxxxxxxxxxxxxxxxxxx
>> Unsubscribe : https://launchpad.net/~yade-users<https://launchpad.net/%7Eyade-users>
>> More help   : https://help.launchpad.net/ListHelp
>>
>
>

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

# Script to monitor the energy of a system
from yade import utils
#-----------------------------------------------------
# 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,label='mat'))
#-----------------------------------------------------
# geometry
#-----------------------------------------------------
# create a random packing in a box space
from yade import pack
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[0],s[1]))
# 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')
from yade import qt
qt.View()
qt.Controller()
#-----------------------------------------------------
# plot some results
#-----------------------------------------------------
from math import *
from yade import plot
plot.plots={'t':('Ek','Eel','Slip','Etot'),'t_':('Slip')}
def monitoring():
	plot.addData(Ek=utils.kineticEnergy(),Eel=dry.elasticEnergy(),Slip=dry.plasticDissipation(),t=O.time,t_=O.time,
		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.materials[0].frictionAngle=radians(25)
O.run(30000,True); plot.plot()

Attachment: Energy.png
Description: PNG image