yade-users team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Bruno Chareyre said:     (by the date of Mon, 05 Jul 2010 18:00:06 +0200)

> >
> >> (If an interaction is broken, then it is a different issue, i.e. that
> >> dissipated energy of broken interactions should be saved somehow, since
> >> otherwise we lose the number).
> >>     
> >
> > the calculation of plasticDissipation is incremental. If the contact
> > breaks, then nothing is lost, you only stop incrementing it.
> >   
> Good point. Something is missed each time a contact is lost, since this 
> contact dissipated something on [t,t+dt] even if it is lost at time t+dt.
> This can be negligible in dense packing (triaxial) but significant in 
> bouncing spheres.
> It doesn't explain the excessive growth of plasticE though...
> 
> I'm wondering if the problem is not just the time integration scheme : 
> if a bouncing is described in 2-3 steps, the approximation of 
> derivatives is horrible, and it might well create energy artificialy. 
> Perhaps Cundall introduced global damping for a reason...

is it enough if we decrease timestep 1000 times?

I made a comparison with jobs=1 and with the same scene (so the
curves should be all identical). They are not.

See attached plot. And also the script (Chiara: it is without
my shearDisp modification, it is using the correct formula).


on unrelated note:

A more technical question, since I am learning to use plot module.
Can I set a plot title? - a big text above the plot box.

-- 
Janek Kozicki                               http://janek.kozicki.pl/  |

Attachment: energy_BasicLaw_comparison.png
Description: PNG image

#!/usr/local/bin/yade-trunk -j 1
# -*- 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.saveTmp()

O.loadTmp()
O.materials[0].frictionAngle=radians(0)
O.run(10000,True)
O.materials[0].frictionAngle=radians(25)
O.run(30000,True)
plot.plot()
plot.splitData()

O.loadTmp()
O.materials[0].frictionAngle=radians(0)
dry.traceEnergy=True
O.run(10000,True); 
O.materials[0].frictionAngle=radians(25)
O.run(30000,True)
plot.plot()
plot.splitData()

O.loadTmp()
O.materials[0].frictionAngle=radians(0)
dry.traceEnergy=True
dry.neverErase=True
O.run(10000,True); 
O.materials[0].frictionAngle=radians(25)
O.run(30000,True)
plot.plot()
plot.splitData()

O.loadTmp()
O.materials[0].frictionAngle=radians(0)
dry.traceEnergy=True
dry.neverErase=True
O.dt=.001*utils.PWaveTimeStep()
O.run(1000000,True)
O.materials[0].frictionAngle=radians(25)
O.run(3000000,True)
plot.plot()
plot.splitData()