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