tadded new methods for other energy balance components - sphere - GPU-based 3D discrete element method algorithm with optional fluid coupling
(HTM) git clone git://src.adamsgaard.dk/sphere
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) LICENSE
---
(DIR) commit 22e628495585c3145e362c918c9e62dab291a147
(DIR) parent d951e69722dbdda4e924831e6ef544c468ad533e
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 20 Jun 2014 13:10:50 +0200
added new methods for other energy balance components
Diffstat:
M python/sphere.py | 67 ++++++++++++++++++++++++++++---
M tests/contactmodel.py | 50 ++++++++++++++++++++++++++++++-
2 files changed, 111 insertions(+), 6 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -2686,9 +2686,20 @@ class sim:
'''
return self.rho[0]*self.volume(idx)
+ def momentOfInertia(self, idx):
+ '''
+ Returns the moment of inertia of a particle.
+
+ :param idx: Particle index
+ :type idx: int
+ :returns: The moment of inertia [kg*m^2]
+ :return type: float
+ '''
+ return 2.0/5.0*self.mass(idx)*self.radius[idx]**2
+
def kineticEnergy(self, idx):
'''
- Returns the kinetic energy for a particle.
+ Returns the (linear) kinetic energy for a particle.
:param idx: Particle index
:type idx: int
t@@ -2700,7 +2711,7 @@ class sim:
def totalKineticEnergy(self):
'''
- Returns the total kinetic energy for all particles.
+ Returns the total linear kinetic energy for all particles.
:returns: The kinetic energy of all particles [J]
'''
t@@ -2709,9 +2720,32 @@ class sim:
esum += self.kineticEnergy(i)
return esum
+ def rotationalEnergy(self, idx):
+ '''
+ Returns the rotational energy for a particle.
+
+ :param idx: Particle index
+ :type idx: int
+ :returns: The rotational kinetic energy of the particle [J]
+ :return type: float
+ '''
+ return 0.5*self.momentOfInertia(idx) \
+ *numpy.sqrt(numpy.dot(self.angvel[idx,:], self.angvel[idx,:]))**2
+
+ def totalRotationalEnergy(self):
+ '''
+ Returns the total rotational kinetic energy for all particles.
+
+ :returns: The rotational energy of all particles [J]
+ '''
+ esum = 0.0
+ for i in range(self.np):
+ esum += self.rotationalEnergy(i)
+ return esum
+
def viscousNormalEnergy(self, idx):
'''
- Returns the viscous absorbed energy for a particle in the
+ Returns the viscous dissipated energy for a particle in the
normal component of its contacts.
:param idx: Particle index
t@@ -2723,17 +2757,40 @@ class sim:
def totalViscousNormalEnergy(self):
'''
- Returns the total viscous absorbed energy for all particles
+ Returns the total viscous dissipated energy for all particles
(normal component).
:returns: The normal viscous energy of all particles [J]
-
+ :return type: float
'''
esum = 0.0
for i in range(self.np):
esum += self.viscousNormalEnergy(i)
return esum
+ def frictionalEnergy(self, idx):
+ '''
+ Returns the frictional dissipated energy for a particle.
+
+ :param idx: Particle index
+ :type idx: int
+ :returns: The frictional energy lost of the particle [J]
+ :return type: float
+ '''
+ return self.es[idx]
+
+ def totalFrictionalEnergy(self):
+ '''
+ Returns the total frictional dissipated energy for all particles.
+
+ :returns: The total frictional energy lost of all particles [J]
+ :return type: float
+ '''
+ esum = 0.0
+ for i in range(self.np):
+ esum += self.frictionalEnergy(i)
+ return esum
+
def energy(self, method):
'''
Calculates the sum of the energy components of all particles.
(DIR) diff --git a/tests/contactmodel.py b/tests/contactmodel.py
t@@ -133,4 +133,52 @@ pytestutils.compareFloats(orig.totalKineticEnergy(),
"Viscoelastic normal collision (4/4):", tolerance=0.05)
-#orig.cleanup()
+
+## Oblique elastic collisions
+
+# Normal impact, low angle, no slip
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+#orig.radius[:] = [1.0, 2.0]
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.vel[0,2] = 1
+orig.vel[0,0] = 1
+orig.mu_s[0] = 1e9 # no slip
+orig.mu_d[0] = 1e9 # no slip
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+print(after.es)
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalRotationalEnergy(),\
+ "Oblique normal collision (1/4):", tolerance=0.05)
+
+# Normal impact, low angle, slip
+orig = sphere.sim(np=2, sid='contactmodeltest')
+after = sphere.sim(np=2, sid='contactmodeltest')
+sphere.cleanup(orig)
+#orig.radius[:] = [1.0, 2.0]
+orig.radius[:] = [1.0, 1.0]
+orig.x[0,:] = [5.0, 5.0, 2.0]
+orig.x[1,:] = [5.0, 5.0, 4.05]
+orig.vel[0,2] = 1
+orig.vel[0,0] = 1
+orig.mu_s[0] = 0.0
+orig.mu_d[0] = 0.0
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+print(after.es)
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalRotationalEnergy(),\
+ "Oblique normal collision (2/4):", tolerance=0.05)
+