yade-users team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

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[0].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
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',traceEnergy=True)]
	),
	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

prevE=28571.0
sum_negative_E=0.0

def only_decreasing():
	if(O.materials[0].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():
	plot.addData(Ek=utils.kineticEnergy(),Eel=dry.elasticEnergy(),Slip=dry.plasticDissipation(),Slip_how=28571-utils.kineticEnergy()-dry.elasticEnergy(),t=O.time,t_=O.time,
		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.materials[0].frictionAngle=radians(25)
O.run(30000,True); plot.plot()

Attachment: scr_3224.png
Description: PNG image