tadd hertz-mendlin contact model, still untested - granular - granular dynamics simulation
(HTM) git clone git://src.adamsgaard.dk/granular
(DIR) Log
(DIR) Files
(DIR) Refs
(DIR) README
(DIR) LICENSE
---
(DIR) commit 35334244861f0274779e5f5df8936931e8649c49
(DIR) parent 1530b4bb345ad9e65cdf43acda89c5b5ca48ded6
(HTM) Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Wed, 21 Apr 2021 21:55:20 +0200
add hertz-mendlin contact model, still untested
Diffstat:
M arrays.c | 15 +++++++++++++--
M arrays.h | 1 +
M grain.c | 35 ++++++++++++++++++++-----------
M grain.h | 4 +++-
M grains.c | 96 ++++++++++++++++++++++++++-----
M grains.h | 1 +
M granular.5 | 4 ++++
M packing.c | 28 ++++++++++++++--------------
M simulation.c | 22 +++++++++++++---------
9 files changed, 155 insertions(+), 51 deletions(-)
---
(DIR) diff --git a/arrays.c b/arrays.c
t@@ -290,8 +290,7 @@ normalize(const double *in, const int n)
double *out;
check_magnitude(__func__, 1, n);
- if (!(out = calloc(n, sizeof(double))))
- err(1, "%s: out calloc", __func__);
+ out = empty(n);
copy_values(in, out, n);
max_val = max(out, n);
t@@ -339,3 +338,15 @@ dot(const double *a, const double *b, const int n)
return out;
}
+
+double *
+cross(const double a[3], const double b[3])
+{
+ double *out = empty(3);
+
+ out[0] = a[1]*b[2] - a[2]*b[1];
+ out[1] = a[2]*b[0] - a[0]*b[2];
+ out[2] = a[0]*b[1] - a[1]*b[0];
+
+ return out;
+}
(DIR) diff --git a/arrays.h b/arrays.h
t@@ -54,5 +54,6 @@ double * normalize(const double *in, const int n);
double euclidean_norm(const double *a, const int n);
double euclidean_distance(const double *a, const double *b, const int n);
double dot(const double *a, const double *b, const int n);
+double * cross(const double a[3], const double b[3]);
#endif
(DIR) diff --git a/grain.c b/grain.c
t@@ -23,7 +23,7 @@ grain_defaults(struct grain *g)
{
int i;
- g->radius = 0.5;
+ g->diameter = 1.0;
for (i = 0; i < 3; i++) {
g->pos[i]
= g->vel[i]
t@@ -46,6 +46,8 @@ grain_defaults(struct grain *g)
g->youngs_modulus = 1e9;
g->poissons_ratio = 0.2;
g->friction_coeff = 0.4;
+ g->damping_n = 0.0;
+ g->damping_t = 0.0;
g->tensile_strength = 1e8;
g->shear_strength = 1e8;
g->fracture_toughness = 1e8;
t@@ -86,7 +88,7 @@ print_padded_nd_uint(FILE *stream, const size_t *arr)
void
grain_print(FILE *stream, const struct grain *g)
{
- fprintf(stream, "%.*g\t", FLOATPREC, g->radius * 2.0);
+ fprintf(stream, "%.*g\t", FLOATPREC, g->diameter);
print_padded_nd_double(stream, g->pos);
print_padded_nd_double(stream, g->vel);
print_padded_nd_double(stream, g->acc);
t@@ -105,6 +107,8 @@ grain_print(FILE *stream, const struct grain *g)
fprintf(stream, "%.*g\t", FLOATPREC, g->youngs_modulus);
fprintf(stream, "%.*g\t", FLOATPREC, g->poissons_ratio);
fprintf(stream, "%.*g\t", FLOATPREC, g->friction_coeff);
+ fprintf(stream, "%.*g\t", FLOATPREC, g->damping_n);
+ fprintf(stream, "%.*g\t", FLOATPREC, g->damping_t);
fprintf(stream, "%.*g\t", FLOATPREC, g->tensile_strength);
fprintf(stream, "%.*g\t", FLOATPREC, g->shear_strength);
fprintf(stream, "%.*g\t", FLOATPREC, g->fracture_toughness);
t@@ -118,13 +122,12 @@ grain_print(FILE *stream, const struct grain *g)
struct grain *
grain_read(char *line)
{
- double diameter;
struct grain *g;
if (!(g = malloc(sizeof(struct grain))))
err(1, "%s: grain malloc", __func__);
- if (sscanf(line, "%lg\t" /* 2.0 * radius */
+ if (sscanf(line, "%lg\t" /* diameter */
"%lg\t%lg\t%lg\t" /* pos */
"%lg\t%lg\t%lg\t" /* vel */
"%lg\t%lg\t%lg\t" /* acc */
t@@ -143,6 +146,8 @@ grain_read(char *line)
"%lg\t" /* youngs_modulus */
"%lg\t" /* poissons_ratio */
"%lg\t" /* friction_coeff */
+ "%lg\t" /* damping_n */
+ "%lg\t" /* damping_t */
"%lg\t" /* tensile_strength */
"%lg\t" /* shear_strength */
"%lg\t" /* fracture_toughness */
t@@ -151,7 +156,7 @@ grain_read(char *line)
"%lg\t%lg\t%lg\t" /* contact_stress */
"%lg\t" /* thermal_energy */
"%d\n", /* color */
- &diameter,
+ &g->diameter,
&g->pos[0], &g->pos[1], &g->pos[2],
&g->vel[0], &g->vel[1], &g->vel[2],
&g->acc[0], &g->acc[1], &g->acc[2],
t@@ -170,6 +175,8 @@ grain_read(char *line)
&g->youngs_modulus,
&g->poissons_ratio,
&g->friction_coeff,
+ &g->damping_n,
+ &g->damping_t,
&g->tensile_strength,
&g->shear_strength,
&g->fracture_toughness,
t@@ -177,11 +184,9 @@ grain_read(char *line)
&g->ncontacts,
&g->contact_stress[0], &g->contact_stress[1], &g->contact_stress[2],
&g->thermal_energy,
- &g->color) != 53)
+ &g->color) != 55)
errx(1, "%s: could not read line: %s", __func__, line);
- g->radius = 0.5 * diameter;
-
if (grain_check_values(g))
errx(1, "%s: invalid values in grain line: %s", __func__, line);
t@@ -194,7 +199,7 @@ grain_check_values(const struct grain *g)
int i;
int status = 0;
- check_float_positive("grain->radius", g->radius, &status);
+ check_float_positive("grain->diameter", g->diameter, &status);
for (i = 0; i < 3; i++) {
check_float("grain->pos", g->pos[i], &status);
t@@ -224,6 +229,12 @@ grain_check_values(const struct grain *g)
check_float_non_negative("grain->friction_coeff",
g->friction_coeff,
&status);
+ check_float_non_negative("grain->damping_n",
+ g->damping_n,
+ &status);
+ check_float_non_negative("grain->damping_t",
+ g->damping_t,
+ &status);
check_float_non_negative("grain->tensile_strength",
g->tensile_strength,
&status);
t@@ -247,7 +258,7 @@ grain_check_values(const struct grain *g)
double
grain_volume(const struct grain *g)
{
- return 4.0 / 3.0 * PI * pow(g->radius, 3.0);
+ return 4.0 / 3.0 * PI * pow(g->diameter / 2.0, 3.0);
}
double
t@@ -259,7 +270,7 @@ grain_mass(const struct grain *g)
double
grain_moment_of_inertia(const struct grain *g)
{
- return 2.0 / 5.0 * grain_mass(g) * pow(g->radius, 2.0);
+ return 2.0 / 5.0 * grain_mass(g) * pow(g->diameter / 2.0, 2.0);
}
void
t@@ -348,7 +359,7 @@ grain_register_contact(struct grain *g, size_t i, size_t j,
}
/* second pass: register as new contact if overlapping */
- if (overlap < 0.0)
+ if (overlap > 0.0)
for (ic = 0; ic <= MAXCONTACTS; ic++) {
if (ic == MAXCONTACTS)
err(1, "%s: contact %zu exceeds MAXCONTACTS (%d)",
(DIR) diff --git a/grain.h b/grain.h
t@@ -8,7 +8,7 @@
struct grain {
- double radius;
+ double diameter;
double pos[3];
double vel[3];
double acc[3];
t@@ -27,6 +27,8 @@ struct grain {
double youngs_modulus;
double poissons_ratio;
double friction_coeff;
+ double damping_n;
+ double damping_t;
double tensile_strength;
double shear_strength;
double fracture_toughness;
(DIR) diff --git a/grains.c b/grains.c
t@@ -1,6 +1,7 @@
#include <stdio.h>
#include <math.h>
#include "grain.h"
+#include "arrays.h"
#define VTK_FLOAT_FMT "%.17g "
t@@ -59,15 +60,7 @@ grains_print_vtk(FILE *stream, const struct grain *grains, size_t ng)
"\t\t\t<PointData Scalars=\"Diameter [m]\" "
"Vectors=\"Angular position [-]\">\n");
- fprintf(stream,
- "\t\t\t\t<DataArray type=\"Float64\" Name=\"Diameter [m]\" "
- "NumberOfComponents=\"1\" format=\"ascii\">\n");
- for (i = 0; i < ng; i++)
- fprintf(stream, VTK_FLOAT_FMT, grains[i].radius * 2.0);
- fprintf(stream,
- "\n"
- "\t\t\t\t</DataArray>\n");
-
+ VTK_XML_SCALAR(diameter, "Diameter [m]", "Float64", VTK_FLOAT_FMT);
VTK_XML_VECTOR(vel, "Velocity [m/s]", "Float64", VTK_FLOAT_FMT);
VTK_XML_VECTOR(acc, "Acceleration [m/s^2]", "Float64", VTK_FLOAT_FMT);
VTK_XML_VECTOR(acc_lock, "Acceleration lock [-]", "Int64", "%d ");
t@@ -85,6 +78,8 @@ grains_print_vtk(FILE *stream, const struct grain *grains, size_t ng)
VTK_XML_SCALAR(youngs_modulus, "Young's modulus [Pa]", "Float64", VTK_FLOAT_FMT);
VTK_XML_SCALAR(poissons_ratio, "Poisson's ratio [-]", "Float64", VTK_FLOAT_FMT);
VTK_XML_SCALAR(friction_coeff, "Friction coefficient [-]", "Float64", VTK_FLOAT_FMT);
+ VTK_XML_SCALAR(damping_n, "Contact-normal damping constant [1/s]", "Float64", VTK_FLOAT_FMT);
+ VTK_XML_SCALAR(damping_t, "Contact-tangential damping constant [1/s]", "Float64", VTK_FLOAT_FMT);
VTK_XML_SCALAR(tensile_strength, "Tensile strength [Pa]", "Float64", VTK_FLOAT_FMT);
VTK_XML_SCALAR(shear_strength, "Shear strength [Pa]", "Float64", VTK_FLOAT_FMT);
VTK_XML_SCALAR(fracture_toughness, "Fracture toughness [Pa]", "Float64", VTK_FLOAT_FMT);
t@@ -109,8 +104,8 @@ grains_minimum_grain_edge(const struct grain *grains, size_t ng, int d)
double res = INFINITY;
for (i = 0; i < ng; i++)
- if (res > grains[i].pos[d] - grains[i].radius)
- res = grains[i].pos[d] - grains[i].radius;
+ if (res > grains[i].pos[d] - grains[i].diameter / 2.0)
+ res = grains[i].pos[d] - grains[i].diameter / 2.0;
return res;
}
t@@ -122,8 +117,83 @@ grains_maximum_grain_edge(const struct grain *grains, size_t ng, int d)
double res = -INFINITY;
for (i = 0; i < ng; i++)
- if (res < grains[i].pos[d] + grains[i].radius)
- res = grains[i].pos[d] + grains[i].radius;
+ if (res < grains[i].pos[d] + grains[i].diameter / 2.0)
+ res = grains[i].pos[d] + grains[i].diameter / 2.0;
return res;
}
+
+/* Hertz-Mindlin contact model for compressible spheres */
+void
+grains_interact(struct grain *g_i, struct grain *g_j, int ic)
+{
+ int d;
+ double f_n_ij[3], f_t_ij[3],
+ v_ij[3], v_n_ij[3], v_t_ij[3], v_n_dot_n_ij,
+ delta_ij,
+ n_ij[3], u_t_ij[3],
+ d_i, d_j,
+ r_ij[3], r_ij_norm,
+ m_i, m_j, m_ij,
+ E_ij, nu_ij, mu_ij,
+ gamma_n_ij, gamma_t_ij,
+ k_n_ij, k_t_ij, A_ij,
+ angvel_ij[3], *angvel_cross_r_ij, *f_t_cross_r_ij;
+
+ delta_ij = g_i->contacts[ic].overlap;
+ d_i = g_i->diameter;
+ d_j = g_j->diameter;
+ m_i = grain_mass(g_i);
+ m_j = grain_mass(g_j);
+ m_ij = m_i * m_j / (m_i + m_j);
+ r_ij_norm = euclidean_norm(g_i->contacts[ic].centerdist, 3);
+
+ for (d = 0; d < 3; d++) {
+ r_ij[d] = g_i->contacts[ic].centerdist[d];
+ n_ij[d] = r_ij[d]/r_ij_norm;
+ v_ij[d] = g_i->vel[d] - g_j->vel[d];
+ angvel_ij[d] = g_i->angvel[d] + g_j->angvel[d];
+ }
+
+ v_n_dot_n_ij = dot(v_ij, n_ij, 3);
+ angvel_cross_r_ij = cross(angvel_ij, r_ij);
+
+ for (d = 0; d < 3; d++) {
+ v_n_ij[d] = v_n_dot_n_ij * n_ij[d];
+ v_t_ij[d] = v_ij[d] - v_n_ij[d] - 0.5 * angvel_cross_r_ij[d];
+ }
+
+
+ E_ij = fmin(g_i->youngs_modulus, g_j->youngs_modulus);
+ nu_ij = fmin(g_i->poissons_ratio, g_j->poissons_ratio);
+ mu_ij = fmin(g_i->friction_coeff, g_j->friction_coeff);
+ gamma_n_ij = fmin(g_i->damping_n, g_j->damping_n);
+ gamma_t_ij = fmin(g_i->damping_t, g_j->damping_t);
+
+ k_n_ij = 2.0 / 3.0 * E_ij / (1.0 - pow(nu_ij, 2.0));
+ k_t_ij = 2.0 * E_ij / (1.0 + nu_ij) * (2.0 - nu_ij);
+
+ A_ij = sqrt(delta_ij) * sqrt(d_i * d_j / (2.0 * (d_i + d_j)));
+
+ for (d = 0; d < 3; d++) {
+ f_n_ij[d] = A_ij * (k_n_ij * delta_ij * n_ij[d]
+ - m_ij * gamma_n_ij * v_n_ij[d]);
+ f_t_ij[d] = A_ij * (-k_t_ij * u_t_ij[d]
+ - m_ij * gamma_t_ij * v_t_ij[d]);
+ }
+
+ if (euclidean_norm(f_t_ij, 3) >= euclidean_norm(f_n_ij, 3) * mu_ij)
+ return; /* TODO: limit f_t_ij according to Coulomb criterion */
+
+ f_t_cross_r_ij = cross(f_t_ij, r_ij);
+
+ for (d = 0; d < 3; d++) {
+ g_i->force[d] += f_n_ij[d] + f_t_ij[d];
+ g_j->force[d] -= f_n_ij[d] + f_t_ij[d];
+ g_i->torque[d] += -0.5 * f_t_cross_r_ij[d];
+ g_j->torque[d] -= -0.5 * f_t_cross_r_ij[d];
+ }
+
+ free(angvel_cross_r_ij);
+ free(f_t_cross_r_ij);
+}
(DIR) diff --git a/grains.h b/grains.h
t@@ -8,5 +8,6 @@ void grains_print_vtk(FILE *stream, const struct grain *grains, size_t ng);
double sim_minimum_grain_edge(const struct grain *grains, size_t ng, int d);
double sim_maximum_grain_edge(const struct grain *grains, size_t ng, int d);
+void grains_interact(struct grain *g_i, struct grain *g_j, int nc);
#endif
(DIR) diff --git a/granular.5 b/granular.5
t@@ -98,6 +98,10 @@ Poisson's ratio [-]
.It
friction coefficient [-]
.It
+contact-normal damping constant [1/s]
+.It
+contact-tangential damping constant [1/s]
+.It
tensile strength [Pa]
.It
shear strength [Pa]
(DIR) diff --git a/packing.c b/packing.c
t@@ -10,18 +10,18 @@
size_t
rectangular_packing(struct grain **grains,
size_t n[3],
- double radius_min, double radius_max,
+ double diameter_min, double diameter_max,
double (*gsd)(double min, double max),
double padding_factor,
double origo[3])
{
size_t i, j, k, l, ng, N = 0;
- double dx_padding = radius_max * 2.0 * padding_factor;
- double dx = radius_max * 2.0 + dx_padding;
+ double dx_padding = diameter_max * padding_factor;
+ double dx = diameter_max + dx_padding;
- if (radius_max < radius_min)
- errx(1, "%s: radius_max (%g) is smaller than radius_min (%g)",
- __func__, radius_max, radius_min);
+ if (diameter_max < diameter_min)
+ errx(1, "%s: diameter_max (%g) is smaller than diameter_min (%g)",
+ __func__, diameter_max, diameter_min);
if (!(*grains = calloc(n[0] * n[1] * n[2], sizeof(struct grain))))
err(1, "%s: grains calloc", __func__);
t@@ -32,7 +32,7 @@ rectangular_packing(struct grain **grains,
N++;
ng = idx3(i, j, k, n[0], n[1]);
grain_defaults(&(*grains)[ng]);
- (*grains)[ng].radius = gsd(radius_min, radius_max);
+ (*grains)[ng].diameter = gsd(diameter_min, diameter_max);
(*grains)[ng].pos[0] = i * dx + 0.5 * dx + origo[0];
(*grains)[ng].pos[1] = j * dx + 0.5 * dx + origo[1];
(*grains)[ng].pos[2] = k * dx + 0.5 * dx + origo[2];
t@@ -51,19 +51,19 @@ rectangular_packing(struct grain **grains,
size_t
triangular_packing(struct grain **grains,
size_t n[3],
- double radius_min, double radius_max,
+ double diameter_min, double diameter_max,
double (*gsd)(double min, double max),
double padding_factor,
double origo[3])
{
size_t i, j, k, l, ng, N = 0;
- double dx_padding = radius_max * 2.0 * padding_factor;
- double dx = radius_max * 2.0 + dx_padding;
+ double dx_padding = diameter_max * padding_factor;
+ double dx = diameter_max + dx_padding;
double dy = dx * sin(PI/3.0);
- if (radius_max < radius_min)
- errx(1, "%s: radius_max (%g) is smaller than radius_min (%g)",
- __func__, radius_max, radius_min);
+ if (diameter_max < diameter_min)
+ errx(1, "%s: diameter_max (%g) is smaller than diameter_min (%g)",
+ __func__, diameter_max, diameter_min);
if (!(*grains = calloc(n[0] * n[1] * n[2], sizeof(struct grain))))
err(1, "%s: grains calloc", __func__);
t@@ -74,7 +74,7 @@ triangular_packing(struct grain **grains,
N++;
ng = idx3(i, j, k, n[0], n[1]);
grain_defaults(&(*grains)[ng]);
- (*grains)[ng].radius = gsd(radius_min, radius_max);
+ (*grains)[ng].diameter = gsd(diameter_min, diameter_max);
(*grains)[ng].pos[0] = i * dx + 0.5 * dx + origo[0];
if (j % 2 == 0)
(*grains)[ng].pos[0] += 0.5 * dx;
(DIR) diff --git a/simulation.c b/simulation.c
t@@ -97,8 +97,8 @@ sim_grain_pair_overlap(struct simulation *sim, size_t i, size_t j)
for (d = 0; d < 3; d++)
pos[d] = sim->grains[i].pos[d] - sim->grains[j].pos[d];
- overlap = euclidean_norm(pos, 3)
- - (sim->grains[i].radius + sim->grains[j].radius);
+ overlap = (sim->grains[i].diameter + sim->grains[j].diameter) / 2.0
+ - euclidean_norm(pos, 3);
return overlap;
}
t@@ -118,8 +118,8 @@ sim_check_add_contact(struct simulation *sim, size_t i, size_t j)
centerdistinv[d] = sim->grains[j].pos[d] - sim->grains[i].pos[d];
}
- overlap = euclidean_norm(centerdist, 3)
- - (sim->grains[i].radius + sim->grains[j].radius);
+ overlap = 0.5*(sim->grains[i].diameter + sim->grains[j].diameter)
+ - euclidean_norm(centerdist, 3);
if (sim->grains[i].enabled)
grain_register_contact(&sim->grains[i], j, i, centerdist, overlap);
t@@ -134,9 +134,10 @@ sim_detect_contacts(struct simulation *sim)
size_t i, j;
for (i = 0; i < sim->ng; i++)
- for (j = 0; j < sim->ng; j++)
- if (i < j)
- sim_check_add_contact(sim, i, j);
+ if (sim->grains[i].enabled)
+ for (j = 0; j < sim->ng; j++)
+ if (i < j)
+ sim_check_add_contact(sim, i, j);
}
void
t@@ -147,8 +148,11 @@ sim_resolve_interactions(struct simulation *sim)
for (i = 0; i < sim->ng; i++)
for (ic = 0; ic < MAXCONTACTS; ic++)
- if (sim->grains[i].contacts[ic].active)
- grain_interact(i);
+ if (sim->grains[i].contacts[ic].active &&
+ i < sim->grains[i].contacts[ic].j)
+ grains_interact(&sim->grains[i],
+ &sim->grains[sim->grains[i].contacts[ic].j],
+ ic);
}
void