tviscous energy in tangential component added to ev - 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 d5d09a59fa14e9a6b46f1a5712d782e76258c352
(DIR) parent 5dfe62b3e9c41649433d64737c1c27459b567940
(HTM) Author: Anders Damsgaard <anders.damsgaard@geo.au.dk>
Date: Fri, 20 Jun 2014 13:21:09 +0200
viscous energy in tangential component added to ev
Diffstat:
M python/sphere.py | 16 +++++++---------
M src/contactmodels.cuh | 5 +++--
M tests/contactmodel.py | 28 ++++++++++++++++++++++++++++
3 files changed, 38 insertions(+), 11 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -2743,29 +2743,27 @@ class sim:
esum += self.rotationalEnergy(i)
return esum
- def viscousNormalEnergy(self, idx):
+ def viscousEnergy(self, idx):
'''
- Returns the viscous dissipated energy for a particle in the
- normal component of its contacts.
+ Returns the viscous dissipated energy for a particle.
:param idx: Particle index
:type idx: int
- :returns: The normal viscous energy of the particle [J]
+ :returns: The energy lost by the particle by viscous dissipation [J]
:return type: float
'''
return self.ev[idx]
- def totalViscousNormalEnergy(self):
+ def totalViscousEnergy(self):
'''
- Returns the total viscous dissipated energy for all particles
- (normal component).
+ Returns the total viscous dissipated energy for all particles.
- :returns: The normal viscous energy of all particles [J]
+ :returns: The normal viscous energy lost by all particles [J]
:return type: float
'''
esum = 0.0
for i in range(self.np):
- esum += self.viscousNormalEnergy(i)
+ esum += self.viscousEnergy(i)
return esum
def frictionalEnergy(self, idx):
(DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
t@@ -385,8 +385,6 @@ __device__ void contactLinear(Float3* F, Float3* T,
// watt = N*m/s = N*s/m * m/s * m/s * s / s
*ev_dot += 0.5 * devC_params.gamma_n * vel_n * vel_n;
-
-
// Add max. capillary force
f_c = -devC_params.kappa * sqrtf(radius_a * radius_b) * n;
t@@ -450,6 +448,9 @@ __device__ void contactLinear(Float3* F, Float3* T,
* (devC_params.mu_d * length(f_n-f_c) * t
+ devC_params.gamma_t * vel_t);
+ // Store energy dissipated in tangential viscous component
+ *ev_dot += 0.5 * devC_params.gamma_t * vel_t * vel_t;
+
// Shear friction heat production rate:
// The energy lost from the tangential spring is dissipated as heat
//*es_dot += -dot(vel_t_ab, f_t);
(DIR) diff --git a/tests/contactmodel.py b/tests/contactmodel.py
t@@ -186,3 +186,31 @@ pytestutils.test((after.angvel[:,1] < 0.0).all(),
"Oblique normal collision (4/5):")
pytestutils.test(after.totalFrictionalEnergy() > 0.0,
"Oblique normal collision (5/5):")
+
+# 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.3
+orig.mu_d[0] = 0.3
+orig.gamma_t[0] = 1.0e6
+orig.defineWorldBoundaries(L=[10,10,10])
+orig.initTemporal(total = 0.1, file_dt = 0.01)
+
+orig.run(verbose=False)
+after.readlast(verbose=False)
+pytestutils.compareFloats(orig.totalKineticEnergy(),
+ after.totalKineticEnergy()
+ + after.totalRotationalEnergy()
+ + after.totalFrictionalEnergy(),
+ "Oblique normal collision (6/5):", tolerance=0.05)
+pytestutils.test((after.angvel[:,1] < 0.0).all(),
+ "Oblique normal collision (7/5):")
+pytestutils.test(after.totalFrictionalEnergy() > 0.0,
+ "Oblique normal collision (8/5):")