GopherProxy

	tallow contacts to break and keep track of contact age - granular - granular dy…
	git clone git://src.adamsgaard.dk/granular
	Log
	Files
	Refs
	README
	LICENSE
	---
	commit 0fe34fa5175dd35819839ff25ac866300fa82713
	parent 35334244861f0274779e5f5df8936931e8649c49
	Author: Anders Damsgaard <[email protected]>
	Date: Thu, 22 Apr 2021 09:52:03 +0200

	allow contacts to break and keep track of contact age

	Diffstat:
	M grains.c \| 23 ++++++++++++++++++-----
	M grains.h \| 2 +-
	M simulation.c \| 10 +++-------

	3 files changed, 22 insertions(+), 13 deletions(-)
	---
	diff --git a/grains.c b/grains.c
	t@@ -125,22 +125,30 @@ grains_maximum_grain_edge(const struct grain *grains, si…

	/* Hertz-Mindlin contact model for compressible spheres */
	void
	-grains_interact(struct grain g_i, struct grain g_j, int ic)
	+grains_interact(struct grain g_i, struct grain g_j, int ic, double dt)
	{
	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],
	+ n_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;
	+ angvel_ij[3], u_t_dot_v_ij,
	+ angvel_cross_r_ij, f_t_cross_r_ij;

	delta_ij = g_i->contacts[ic].overlap;
	+
	+ if (delta_ij < 0.0) { /* TODO: implement tensile strength */
	+ g_i->contacts[ic].active = 0;
	+ g_i->ncontacts--;
	+ return;
	+ }
	+
	d_i = g_i->diameter;
	d_j = g_j->diameter;
	m_i = grain_mass(g_i);
	t@@ -163,7 +171,6 @@ grains_interact(struct grain g_i, struct grain g_j, int …
	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);
	t@@ -178,7 +185,7 @@ grains_interact(struct grain g_i, struct grain g_j, int …
	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 * …
	- f_t_ij[d] = A_ij * (-k_t_ij * u_t_ij[d]
	+ f_t_ij[d] = A_ij * (-k_t_ij * g_i->contacts[ic].tandisp[d]
	- m_ij * gamma_t_ij * v_t_ij[d]);
	}

	t@@ -194,6 +201,12 @@ grains_interact(struct grain g_i, struct grain g_j, int…
	g_j->torque[d] -= -0.5 * f_t_cross_r_ij[d];
	}

	+ g_i->contacts[ic].age += dt;
	+ u_t_dot_v_ij = dot(g_i->contacts[ic].tandisp, v_ij, 3);
	+ for (d = 0; d < 3; d++)
	+ g_i->contacts[ic].tandisp[d] += v_t_ij[d] * dt
	+ - u_t_dot_v_ij / (r_ij_norm * r_…
	+
	free(angvel_cross_r_ij);
	free(f_t_cross_r_ij);
	}
	diff --git a/grains.h b/grains.h
	t@@ -8,6 +8,6 @@ void grains_print_vtk(FILE stream, const struct grain grains…

	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);
	+void grains_interact(struct grain g_i, struct grain g_j, int nc, double dt);

	#endif
	diff --git a/simulation.c b/simulation.c
	t@@ -107,24 +107,20 @@ static void
	sim_check_add_contact(struct simulation *sim, size_t i, size_t j)
	{

	- double overlap, centerdist[3], centerdistinv[3];
	+ double overlap, centerdist[3];
	int d;

	if (i >= j)
	return;

	- for (d = 0; d < 3; d++) {
	+ for (d = 0; d < 3; d++)
	centerdist[d] = sim->grains[i].pos[d] - sim->grains[j].pos[d];
	- centerdistinv[d] = sim->grains[j].pos[d] - sim->grains[i].pos[…
	- }

	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, over…
	- if (sim->grains[j].enabled)
	- grain_register_contact(&sim->grains[j], i, j, centerdistinv, o…
	}

	/* TODO: add grid sorting */
	t@@ -152,7 +148,7 @@ sim_resolve_interactions(struct simulation *sim)
	i < sim->grains[i].contacts[ic].j)
	grains_interact(&sim->grains[i],
	&sim->grains[sim->grains[i].co…
	- ic);
	+ ic, sim->dt);
	}

	void