commit 30733e55ce95327404e4d55f8cd679f194c8447a
parent 32c8ce3cfe0cb095da8786a5eb2150b0777fc50a
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 30 Aug 2021 09:29:10 +0200
allow negative values of sigma_n_eff, but trunkate for granular solver
Diffstat:
1 file changed, 12 insertions(+), 7 deletions(-)
diff --git a/simulation.c b/simulation.c
@@ -15,6 +15,9 @@
/* tolerance criteria when in velocity driven or velocity limited mode */
#define RTOL_VELOCITY 1e-3
+/* lower limit for effective normal stress sigma_n_eff for granular solver */
+#define SIGMA_N_EFF_MIN 1e-1
+
/* Simulation settings */
void
@@ -430,7 +433,7 @@ compute_inertia_number(struct simulation *sim)
for (i = 0; i < sim->nz; ++i)
sim->I[i] = inertia_number(sim->gamma_dot_p[i],
sim->d,
- sim->sigma_n_eff[i],
+ fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
sim->rho_s);
}
@@ -451,8 +454,9 @@ compute_critical_state_friction(struct simulation *sim)
if (sim->fluid)
for (i = 0; i < sim->nz; ++i)
- sim->mu_c[i] = sim->mu_wall / (sim->sigma_n_eff[i]
- / (sim->P_wall - sim->p_f_top));
+ sim->mu_c[i] = sim->mu_wall
+ / (fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN)
+ / (sim->P_wall - sim->p_f_top));
else
for (i = 0; i < sim->nz; ++i)
sim->mu_c[i] = sim->mu_wall
@@ -551,7 +555,7 @@ compute_effective_stress(struct simulation *sim)
/ 2.0);
/* sim->sigma_n_eff[i] = sim->sigma_n[i] - sim->p_f_ghost[i + 1]; */
if (sim->sigma_n_eff[i] < 0)
- errx(1, "%s: sigma_n_eff[%d] is negative with value %g\n",
+ warnx("%s: sigma_n_eff[%d] is negative with value %g\n",
__func__, i, sim->sigma_n_eff[i]);
}
else
@@ -579,7 +583,7 @@ compute_cooperativity_length(struct simulation *sim)
sim->xi[i] = cooperativity_length(sim->A,
sim->d,
sim->mu[i],
- sim->sigma_n_eff[i],
+ fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
sim->mu_s,
sim->C);
}
@@ -605,7 +609,8 @@ compute_local_fluidity(struct simulation *sim)
int i;
for (i = 0; i < sim->nz; ++i)
- sim->g_local[i] = local_fluidity(sim->sigma_n_eff[i],
+ sim->g_local[i] = local_fluidity(fmax(sim->sigma_n_eff[i],
+ SIGMA_N_EFF_MIN),
sim->mu[i],
sim->mu_s,
sim->C,
@@ -770,7 +775,7 @@ print_output(struct simulation *sim, FILE *fp, const int norm)
sim->gamma_dot_p[i],
sim->phi[i],
sim->I[i],
- sim->mu[i] * sim->sigma_n_eff[i],
+ sim->mu[i] * fmax(sim->sigma_n_eff[i], SIGMA_N_EFF_MIN),
sim->d_x[i]);
free(v_x_out);