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>>
#!/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()
_______________________________________________ Mailing list: https://launchpad.net/~yade-users Post to : [email protected] Unsubscribe : https://launchpad.net/~yade-users More help : https://help.launchpad.net/ListHelp
