yade-dev team mailing list archive

[chiara modenese <c.modenese@xxxxxxxxx>]

Hi all (Bruno?),


I spent a bit of time on energy considerations. I remember Bruno saying that
in the TT energy is conserved (could you post the example you tried,
please?). I post a simple script which plots energy terms and total energy
as well. I also post the graph I obtained (your output if you run the
script). I am using default parameters, and NO damping at all. The total
energy seems to increase in the triaxial phase. Why?
I checked in the code the way we calculate the external work done by the
walls and it looks correct. Look at the unbalanced force in the graph. My
conclusion is that as long as the system is dynamic, energy appears to be
not conserved. If I am not wrong, this could be a problem for people running
simulation in a dynamic state rather than in quasi-static equilibrium. Still
unclear why to me.


If you have a bit of time just have a look, the script could not be simpler.



Cheers. Chiara

Attachment: Triax_energy.png
Description: PNG image

#!/usr/local/bin/yade-trunk -x
# -*- coding: utf-8 -*-
# -*- encoding=utf-8 -*-

# Script to test energy conservation in the 'default' Triaxial Test
TriaxialTest(dampingForce=0.0,dampingMomentum=0.0).load() # NO Cundall'S DAMPING


contactLaw=O.engines[3].lawDispatcher.functors[0] # get contact law
contactLaw.traceEnergy=True # calculate energy terms
ttcomp=O.engines[5] # get TriaxialCompressionEngine class
O.engines=O.engines+[PeriodicPythonRunner(iterPeriod=10,command='energy_Basic()')] # add PeriodicEngine


def energy_Basic(): 
	plot.addData(
		workWalls=ttcomp.externalWork, # external work
		Ek=utils.kineticEnergy(), # kinetic energy
		EnplusEs=contactLaw.elasticEnergy(), # elastic potential energy (both normal and shear)
		Df=contactLaw.plasticDissipation(), # sliding dissipated energy
		Etot=ttcomp.externalWork+utils.kineticEnergy()+contactLaw.elasticEnergy()+contactLaw.plasticDissipation(), # total energy
		unbalanced=ttcomp.UnbalancedForce, # Unbalanced force
		t=O.time)


from yade import plot
plot.plots={'t':('workWalls','Ek','EnplusEs','Df',('Etot','g--'),'|||',('unbalanced'))}
plot.plot(); O.run(30000, True);