yade-users team mailing list archive

[ytang <question693852@xxxxxxxxxxxxxxxxxxxxx>]

New question #693852 on Yade:
https://answers.launchpad.net/yade/+question/693852

Hi all,

recently, I'm trying to study energy consumption in YADE, so I just use the gravity deposition process to check the energy dissipation process.

here is the code: the code is almost the same as the example in YADE [1].
#########################################################################
from yade import pack, plot
damping = 0.4
O.bodies.append(geom.facetBox((.5,.5,.5),(.5,.5,.5),wallMask=31))
sp=pack.SpherePack()
sp.makeCloud((0,0,0),(1,1,1),rMean=.23,rRelFuzz=0.0,seed =1)
sp.toSimulation()
print(len(O.bodies))

for b in O.bodies:
	if isinstance(b.shape,Sphere):	
		positionz = b.state.pos[2]
		mass = b.state.mass
		print(positionz)
		print(mass)
		print("#######")		
print("############################################")
for b in O.bodies:
	if isinstance(b.shape,Sphere):	
		radius = b.shape.radius
		print(radius)
O.engines=[
	ForceResetter(),
	InsertionSortCollider([Bo1_Sphere_Aabb(),Bo1_Facet_Aabb()]),
	InteractionLoop(
		# handle sphere+sphere and facet+sphere collisions
		[Ig2_Sphere_Sphere_ScGeom(),Ig2_Facet_Sphere_ScGeom()],
		[Ip2_FrictMat_FrictMat_FrictPhys()],
		[Law2_ScGeom_FrictPhys_CundallStrack(traceEnergy = True, label = 'law')]
		#[Ip2_FrictMat_FrictMat_MindlinPhys( betan =0.2 , betas = 0.2, krot = 0.2 , eta = 0.5  )],
		#[Law2_ScGeom_MindlinPhys_Mindlin()]
	),
	#NewtonIntegrator(gravity=(0,0,-9.81),kinSplit=True,damping=damping),
	NewtonIntegrator(gravity=(0,0,-9.81),damping=0.4),	
	PyRunner(command='checkUnbalanced()',realPeriod=2),
	PyRunner(command='addPlotData()',iterPeriod=2000),
	#PyRunner(command='record_all_info()',iterPeriod=5000)	
]
O.dt=.5*PWaveTimeStep()
O.trackEnergy=True
def checkUnbalanced():
	if unbalancedForce()<.001:
		O.pause()
		plot.saveDataTxt('bbb.txt.bz2')
######################################################################################
######################################################################################
###################### method 1, energy 1 ##############################

"""
def addPlotData():
	# each item is given a names, by which it can be the unsed in plot.plots
	# the **O.energy converts dictionary-like O.energy to plot.addData arguments
	elastic_part = law.elasticEnergy()
	plastic_part = law.plasticDissipation()
	E_kin_translation = 0
	E_kin_rotation    = 0
	E_pot		  = 0
	E_plastic	  = 0	
	E_nonviscoDamp =	0
	E_tracker	  = dict(list(O.energy.items()))

	E_kin_translation = E_tracker['kinTrans']
	E_kin_rotation    = E_tracker['kinRot']
	E_pot             = E_tracker['gravWork']
	if('plastDissip' in E_tracker):
		E_plastic	  = E_tracker['plastDissip']
	if(damping!=0):
		E_nonviscoDamp = E_tracker['nonviscDamp']
	plot.addData( elasticEnergy = elastic_part, plasticenergy = plastic_part,
		translationenergy = E_kin_translation, rotationenergy = E_kin_rotation,
		gravitywork = E_pot, plasticdissip = E_plastic, nonviscodamp = E_nonviscoDamp,
		i=O.iter,unbalanced=unbalancedForce())
    	plot.saveDataTxt('energy.txt')
"""
#####################################################################################
################### method 2, energy 2 #####################################
def addPlotData():
	plot.addData(i=O.iter,unbalanced=unbalancedForce(),**O.energy)
plot.plots={'i':('unbalanced',None,O.energy.keys)}
plot.plot()
O.saveTmp()
################################################################################
################################################################################
I use two methods to get the data, as you can see the annotation from the code. I put the data in this folder [2].
###################
I have some basic questions.
(1) how do we calculate the gravity work? 
In my simulation, there are four particles.  As we all know, the particles should have the maximum gravity work at the height position (which means at the start position). 
I get the z position of the particles and use this equation (m*g*h) to calculate the gravity work. The result is 1022.689 for these four particles. However, the recorded data shows that gravity work changes from -461.95 to -510.83. (as for -461.95, this data is recorded at timesteps = 2000. ) I find the source code on how to calculate the gravity work [3], but I do not fully understand it.

(2) I plot the curves on the energy changes with the timesteps. it seems the energy is not conserved. I do not know whether the energy in YADE is conserved. Can anyone tell me about this?

The reason why I say the energy is not conserved: if we add all the energy at every certain timestep, the result is not zero. for instance : 
i	elastPotential	gravWork	kinRot	kinTrans	nonviscDamp
72000	0.344990744	-510.833081	8.65006E-28	2.57632E-27	291.8306236
Here the total energy is not zero.

(3) I compare the results from method one and method two at the final stage:
##################################################################
energy1	i	elasticEnergy	gravitywork	rotationenergy	translationenergy	nonviscodamp	plasticdissip	plasticenergy
	62000	0.344990744	-510.833081	1.44064E-27	4.73862E-27	291.8306236	0	145.4170567
								
								
energy2	i	elastPotential	gravWork	kinRot	kinTrans	nonviscDamp		
	72000	0.344990744	-510.833081	8.65006E-28	2.57632E-27	291.8306236		
###############################################################################								
The results for these two methods are almost the same for some of the energies. However, there is no plastic dissipation in method 2 (although plastic dissipation in method one is zero). 
In my opinion,  plastDissip can be automatically recorded when we use this command (O.trackEnergy=True) even if the value is zero. 
why the plastDissip does not being recorded here? 
what is the relationship between the plastDissip and the plasticDissipation[4]? (it seems they are not the same from method one. one is zero, the other one is 145.41. )

I'm sorry that I have several questions.
thanks!



references: [1] https://yade-dem.org/doc/tutorial-examples.html
[2]https://www.dropbox.com/sh/fd13w0yednmiz8g/AABkLkUYBiInTSVZpKZmcSYya?dl=0
[3] https://gricad-gitlab.univ-grenoble-alpes.fr/cailletr/yade/-/blob/a583fdbea2c5b63f5d26b96655a6e05e8bccf375/pkg/dem/NewtonIntegrator.cpp
[4] https://yade-dem.org/doc/yade.wrapper.html?highlight=law2_scgeom_frictphys_cundallstrack#yade.wrapper.Law2_ScGeom_FrictPhys_CundallStrack

-- 
You received this question notification because your team yade-users is
an answer contact for Yade.