tuse Youngs' modulus for scale-invariant contact stiffness - 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 ab4f66ef7e057437d51be8d867b0a1eb4ffe6f21
(DIR) parent 7c54905221253a4bee645c39829c788eb3121eac
(HTM) Author: Anders Damsgaard Christensen <adc@geo.au.dk>
Date: Tue, 16 Aug 2016 11:52:51 -0700
use Youngs' modulus for scale-invariant contact stiffness
Diffstat:
M python/sphere.py | 33 +++++++++++++++++++++++++++++--
M src/contactmodels.cuh | 51 +++++++++++++++++++++++++------
M src/datatypes.h | 1 +
M src/device.cu | 2 ++
M src/file_io.cpp | 2 ++
M src/sphere.cpp | 7 +++++++
M src/version.h | 2 +-
7 files changed, 85 insertions(+), 13 deletions(-)
---
(DIR) diff --git a/python/sphere.py b/python/sphere.py
t@@ -24,7 +24,7 @@ numpy.seterr(all='warn', over='raise')
# Sphere version number. This field should correspond to the value in
# `../src/version.h`.
-VERSION = 2.12
+VERSION = 2.13
# Transparency on plot legends
legend_alpha = 0.5
t@@ -161,6 +161,11 @@ class sim:
# rotations. UNUSED
self.k_r = numpy.zeros(1, dtype=numpy.float64)
+ # Young's modulus for contact stiffness [Pa]. This value is used
+ # instead of the Hookean stiffnesses (k_n, k_t) when self.E is larger
+ # than 0.0.
+ self.E = numpy.zeros(1, dtype=numpy.float64)
+
# The viscosity normal to the contact [N/(m/s)]
self.gamma_n = numpy.zeros(1, dtype=numpy.float64)
t@@ -524,6 +529,9 @@ class sim:
elif self.k_r != other.k_r:
print('k_r')
return 31
+ elif self.E != other.E:
+ print('E')
+ return 31.5
elif self.gamma_n != other.gamma_n:
print('gamma_n')
return 32
t@@ -1042,6 +1050,10 @@ class sim:
self.k_n = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.k_t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.k_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ if self.version >= 2.13:
+ self.E = numpy.fromfile(fh, dtype=numpy.float64, count=1)
+ else:
+ self.E = numpy.zeros(1, dtype=numpy.float64)
self.gamma_n = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.gamma_t = numpy.fromfile(fh, dtype=numpy.float64, count=1)
self.gamma_r = numpy.fromfile(fh, dtype=numpy.float64, count=1)
t@@ -1363,6 +1375,7 @@ class sim:
fh.write(self.k_n.astype(numpy.float64))
fh.write(self.k_t.astype(numpy.float64))
fh.write(self.k_r.astype(numpy.float64))
+ fh.write(self.E.astype(numpy.float64))
fh.write(self.gamma_n.astype(numpy.float64))
fh.write(self.gamma_t.astype(numpy.float64))
fh.write(self.gamma_r.astype(numpy.float64))
t@@ -3239,7 +3252,11 @@ class sim:
r_min = numpy.min(self.radius)
m_min = self.rho * 4.0/3.0*numpy.pi*r_min**3
- k_max = numpy.max([self.k_n[:], self.k_t[:]])
+
+ if self.E > 0.001:
+ k_max = numpy.max(numpy.pi/2.0*self.E*self.radius)
+ else:
+ k_max = numpy.max([self.k_n[:], self.k_t[:]])
# Radjaii et al 2011
self.time_dt[0] = epsilon/(numpy.sqrt(k_max/m_min))
t@@ -3766,6 +3783,18 @@ class sim:
'''
self.k_t[0] = k_t
+ def setYoungsModulus(self, E):
+ '''
+ Set the elastic Young's modulus (`E`) for the contact model. This
+ parameter is used over normal stiffness (`k_n`) and tangential
+ stiffness (`k_t`) when its value is greater than zero. Using this
+ parameter produces size-invariant behavior.
+
+ :param E: The elastic modulus [Pa]
+ :type E: float
+ '''
+ self.E[0] = E
+
def setDampingNormal(self, gamma, over_damping=False):
'''
Set the dampening coefficient (gamma) in the normal direction of the
(DIR) diff --git a/src/contactmodels.cuh b/src/contactmodels.cuh
t@@ -53,8 +53,18 @@ __device__ Float contactLinear_wall(
const Float3 vel_t = vel - n * (dot(n, vel));
const Float vel_t_length = length(vel_t);
+ // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
+ // 2013) based on Young's modulus if params.E > 0.0.
+ // Use the calculated stiffness for normal and tangential components.
+ Float k_n = devC_params.k_n;
+ Float k_t = devC_params.k_t;
+ if (devC_params.E > .001) {
+ k_n = 3.141592654/2.0 * devC_params.E * (radius_a + radius_b)/2.;
+ k_t = k_n;
+ }
+
// Normal force component: Elastic - viscous damping
- Float3 f_n = fmax(0.0, -devC_params.k_n*delta
+ Float3 f_n = fmax(0.0, -k_n*delta
- devC_params.gamma_wn*vel_n) * n;
const Float f_n_length = length(f_n); // Save length for later use
t@@ -168,8 +178,18 @@ __device__ void contactLinearViscous(
Float3 vel_t_ab = vel_ab - (n_ab * dot(vel_ab, n_ab));
Float vel_t_ab_length = length(vel_t_ab);
+ // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
+ // 2013) based on Young's modulus if params.E > 0.0.
+ // Use the calculated stiffness for normal and tangential components.
+ Float k_n = devC_params.k_n;
+ Float k_t = devC_params.k_t;
+ if (devC_params.E > .001) {
+ k_n = 3.141592654/2.0 * devC_params.E * (radius_a + radius_b)/2.;
+ k_t = k_n;
+ }
+
// Normal force component: Elastic - viscous damping
- f_n = fmax(0.0, -devC_params.k_n * delta_ab
+ f_n = fmax(0.0, -k_n * delta_ab
- devC_params.gamma_n * vel_n_ab) * n_ab;
// Make sure the viscous damping doesn't exceed the elastic component,
t@@ -304,7 +324,7 @@ __device__ void contactLinear(
Float3 delta_t;
// Normal force component: Elastic - viscous damping
- f_n = fmax(0.0, -devC_params.k_n*delta + devC_params.gamma_n * vel_n) * n;
+ f_n = fmax(0.0, -k_n*delta + devC_params.gamma_n * vel_n) * n;
Float f_n_length = length(f_n);
// Store energy dissipated in normal viscous component
t@@ -328,7 +348,7 @@ __device__ void contactLinear(
delta_t = delta_t0 + vel_t * devC_dt;
// Tangential force: Visco-Elastic, before limitation criterion
- Float3 f_t_elast = -devC_params.k_t * delta_t;
+ Float3 f_t_elast = -k_t * delta_t;
Float3 f_t_visc = -devC_params.gamma_t * vel_t;
f_t = f_t_elast + f_t_visc;
Float f_t_length = length(f_t);
t@@ -356,7 +376,7 @@ __device__ void contactLinear(
// In the sliding friction case, the tangential spring is adjusted
// to a length consistent with Coulombs (dynamic) condition (Luding
// 2008)
- delta_t = -1.0/devC_params.k_t
+ delta_t = -1.0/k_t
* (devC_params.mu_d * length(f_n-f_c) * t
+ devC_params.gamma_t * vel_t);
t@@ -463,9 +483,19 @@ __device__ void contactHertz(
// New tangential displacement vector
Float3 delta_t;
+ // Use scale-invariant contact model (Ergenzinger et al 2010, Obermayr et al
+ // 2013) based on Young's modulus if params.E > 0.0.
+ // Use the calculated stiffness for normal and tangential components.
+ Float k_n = devC_params.k_n;
+ Float k_t = devC_params.k_t;
+ if (devC_params.E > .001) {
+ k_n = 3.141592654/2.0 * devC_params.E * (radius_a + radius_b)/2.;
+ k_t = k_n;
+ }
+
// Normal force component
- f_n = (-devC_params.k_n * powf(delta_ab, 3.0f/2.0f)
- -devC_params.gamma_n * powf(delta_ab, 1.0f/4.0f) * vel_n_ab)
+ f_n = (-k_n * powf(delta_ab, 3.0f/2.0f) -devC_params.gamma_n *
+ powf(delta_ab, 1.0f/4.0f) * vel_n_ab)
* n_ab;
// Store energy dissipated in normal viscous component
t@@ -494,7 +524,7 @@ __device__ void contactHertz(
if (delta_t0_length > 0.f || vel_t_ab_length > 0.f) {
// Shear force: Visco-Elastic, limited by Coulomb friction
- Float3 f_t_elast = -devC_params.k_t * powf(delta_ab, 1.0f/2.0f) * delta_t0;
+ Float3 f_t_elast = -k_t * powf(delta_ab, 1.0f/2.0f) * delta_t0;
Float3 f_t_visc = -devC_params.gamma_t * powf(delta_ab, 1.0f/4.0f) * vel_t_ab;
Float f_t_limit;
t@@ -523,12 +553,13 @@ __device__ void contactHertz(
// length which is consistent with Coulomb's equation
// (Hinrichsen and Wolf, 2004)
delta_t = (f_t + devC_params.gamma_t * powf(delta_ab, 1.0f/4.0f) * vel_t_ab)
- / (-devC_params.k_t * powf(delta_ab, 1.0f/2.0f));
+ / (-k_t * powf(delta_ab, 1.0f/2.0f));
// Shear friction heat production rate:
// The energy lost from the tangential spring is dissipated as heat
//*es_dot += -dot(vel_t_ab, f_t);
- *es_dot += length(delta_t0 - delta_t) * devC_params.k_t / devC_dt; // Seen in EsyS-Particle
+ *es_dot += length(delta_t0 - delta_t) * k_t / devC_dt; // Seen in
+ EsyS-Particle
//*es_dot += fabs(dot(delta_t0 - delta_t, f_t)) / devC_dt;
} else { // Static case
(DIR) diff --git a/src/datatypes.h b/src/datatypes.h
t@@ -75,6 +75,7 @@ struct Params {
Float k_n; // Normal stiffness
Float k_t; // Tangential stiffness
Float k_r; // Rotational stiffness
+ Float E; // Young's modulus
Float gamma_n; // Normal viscosity
Float gamma_t; // Tangential viscosity
Float gamma_r; // Rotational viscosity
(DIR) diff --git a/src/device.cu b/src/device.cu
t@@ -233,6 +233,8 @@ __global__ void checkConstantValues(int* dev_equal,
*dev_equal = 28; // Not ok
if (dev_params->nb0 != devC_params.nb0)
*dev_equal = 29; // Not ok
+ if (dev_params->E != devC_params.E)
+ *dev_equal = 30; // Not ok
}
__global__ void checkParticlePositions(
(DIR) diff --git a/src/file_io.cpp b/src/file_io.cpp
t@@ -172,6 +172,7 @@ void DEM::readbin(const char *target)
ifs.read(as_bytes(params.k_n), sizeof(params.k_n));
ifs.read(as_bytes(params.k_t), sizeof(params.k_t));
ifs.read(as_bytes(params.k_r), sizeof(params.k_r));
+ ifs.read(as_bytes(params.E), sizeof(params.E));
ifs.read(as_bytes(params.gamma_n), sizeof(params.gamma_n));
ifs.read(as_bytes(params.gamma_t), sizeof(params.gamma_t));
ifs.read(as_bytes(params.gamma_r), sizeof(params.gamma_r));
t@@ -497,6 +498,7 @@ void DEM::writebin(const char *target)
ofs.write(as_bytes(params.k_n), sizeof(params.k_n));
ofs.write(as_bytes(params.k_t), sizeof(params.k_t));
ofs.write(as_bytes(params.k_r), sizeof(params.k_r));
+ ofs.write(as_bytes(params.E), sizeof(params.E));
ofs.write(as_bytes(params.gamma_n), sizeof(params.gamma_n));
ofs.write(as_bytes(params.gamma_t), sizeof(params.gamma_t));
ofs.write(as_bytes(params.gamma_r), sizeof(params.gamma_r));
(DIR) diff --git a/src/sphere.cpp b/src/sphere.cpp
t@@ -295,6 +295,11 @@ void DEM::checkValues(void)
exit(1);
}
+ if (params.E < 0.0) {
+ cerr << "Error: E = " << params.E << " N/m^3" << endl;
+ exit(1);
+ }
+
if (params.rho <= 0.0) {
cerr << "Error: rho = " << params.rho << " kg/m3" << endl;
exit(1);
t@@ -827,6 +832,8 @@ void DEM::forcechains(const std::string format, const int threedim,
// Normal force on contact
f_n = -params.k_n * delta_n;
+ // TODO: Use Young's modulus if it is greater than 0
+
if (f_n < lower_cutoff)
continue; // skip the rest of this iteration
(DIR) diff --git a/src/version.h b/src/version.h
t@@ -2,6 +2,6 @@
#define VERSION_H_
// Define source code version
-const double VERSION = 2.12;
+const double VERSION = 2.13;
#endif