commit 6849a3075184bf86a78f553945099b8c3c250d58
parent 4d51b11991631dd7cb45ee6395c5288268bab329
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 8 Jul 2019 12:19:07 +0200
Add cohesion as argument and in granular flow computations
Diffstat:
4 files changed, 41 insertions(+), 11 deletions(-)
diff --git a/main.c b/main.c
@@ -30,6 +30,7 @@ usage(void)
" -A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default %g)\n"
" -b, --rate-dependence VAL rate dependence beyond yield [-] (default %g)\n"
" -f, --friction-coefficient VAL grain friction coefficient [-] (default %g)\n"
+ " -C, --cohesion VAL grain cohesion [Pa] (default %g)\n"
" -p, --porosity VAL porosity fraction [-] (default %g)\n"
" -d, --grain-size VAL representative grain size [m] (default %g)\n"
" -r, --density VAL grain material density [kg/m^3] (default %g)\n"
@@ -57,6 +58,7 @@ usage(void)
sim.A,
sim.b,
sim.mu_s,
+ sim.C,
sim.phi[0],
sim.d,
sim.rho_s,
@@ -102,7 +104,7 @@ main(int argc, char* argv[])
norm = 0;
- optstring = "hvNn:G:P:m:V:A:b:f:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T:D:I:";
+ optstring = "hvNn:G:P:m:V:A:b:f:C:Fp:d:r:o:L:c:i:R:k:O:a:q:H:t:T:D:I:";
const struct option longopts[] = {
{"help", no_argument, NULL, 'h'},
{"version", no_argument, NULL, 'v'},
@@ -114,6 +116,7 @@ main(int argc, char* argv[])
{"nonlocal-amplitude", required_argument, NULL, 'A'},
{"rate-dependence", required_argument, NULL, 'b'},
{"friction-coefficient", required_argument, NULL, 'f'},
+ {"cohesion", required_argument, NULL, 'C'},
{"porosity", required_argument, NULL, 'p'},
{"grain-size", required_argument, NULL, 'd'},
{"density", required_argument, NULL, 'r'},
@@ -180,6 +183,9 @@ main(int argc, char* argv[])
case 'f':
sim.mu_s = atof(optarg);
break;
+ case 'C':
+ sim.C = atof(optarg);
+ break;
case 'p':
new_phi = atof(optarg);
break;
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -32,6 +32,7 @@ struct simulation init_sim(void)
sim.b = 0.9377; /* Henann and Kamrin 2016 */
sim.mu_s = atan(DEG2RAD(22.0)); /* Damsgaard et al 2013 */
+ sim.C = 0.0; /* Damsgaard et al 2013 */
sim.phi = initval(0.25, 1); /* Damsgaard et al 2013 */
diff --git a/simulation.c b/simulation.c
@@ -99,6 +99,11 @@ check_simulation_parameters(const struct simulation* sim)
warn_parameter_value("sim.mu_s is negative", sim->mu_s,
&return_status);
+ check_float("sim.C", sim->C, &return_status);
+ if (sim->C < 0.0)
+ warn_parameter_value("sim.C is negative", sim->C,
+ &return_status);
+
check_float("sim.d", sim->d, &return_status);
if (sim->d <= 0.0)
warn_parameter_value("sim.d is not a positive number", sim->d,
@@ -163,6 +168,11 @@ check_simulation_parameters(const struct simulation* sim)
if (sim->p_f_mod_ampl < 0.0)
warn_parameter_value("sim.p_f_mod_ampl is not a zero or positive",
sim->p_f_mod_ampl, &return_status);
+
+ if (sim->P_wall - sim->p_f_mod_ampl < 0.0)
+ warn_parameter_value("sim.P_wall - sim.p_f_mod_ampl is negative",
+ sim->P_wall - sim->p_f_mod_ampl,
+ &return_status);
check_float("sim.p_f_mod_freq", sim->p_f_mod_freq, &return_status);
if (sim->p_f_mod_freq < 0.0)
@@ -261,13 +271,15 @@ compute_effective_stress(struct simulation* sim)
sim->sigma_n_eff[i] = sim->sigma_n[i];
}
-double
+static double
cooperativity_length(const double A,
const double d,
const double mu,
- const double mu_s)
+ const double p,
+ const double mu_s,
+ const double C)
{
- return A*d/sqrt(fabs(mu - mu_s));
+ return A*d/sqrt(fabs((mu - C/p) - mu_s));
}
void
@@ -275,22 +287,27 @@ compute_cooperativity_length(struct simulation* sim)
{
int i;
for (i=0; i<sim->nz; ++i)
- sim->xi[i] = cooperativity_length(sim->A, sim->d, sim->mu[i],
- sim->mu_s);
+ sim->xi[i] = cooperativity_length(sim->A,
+ sim->d,
+ sim->mu[i],
+ sim->sigma_n_eff[i],
+ sim->mu_s,
+ sim->C);
}
-double
+static double
local_fluidity(const double p,
const double mu,
const double mu_s,
+ const double C,
const double b,
const double rho_s,
const double d)
{
- if (mu <= mu_s)
+ if (mu - C/p <= mu_s)
return 0.0;
else
- return sqrt(p/rho_s*d*d) * (mu - mu_s)/(b*mu);
+ return sqrt(p/rho_s*d*d) * ((mu - C/p) - mu_s)/(b*mu);
}
void
@@ -301,6 +318,7 @@ compute_local_fluidity(struct simulation* sim)
sim->g_ghost[idx1g(i)] = local_fluidity(sim->sigma_n_eff[i],
sim->mu[i],
sim->mu_s,
+ sim->C,
sim->b,
sim->rho_s,
sim->d);
@@ -339,7 +357,7 @@ set_bc_dirichlet(double* g_ghost,
}
}
-void
+static void
poisson_solver_1d_cell_update(int i,
const double* g_in,
double* g_out,
@@ -349,6 +367,7 @@ poisson_solver_1d_cell_update(int i,
const double* p,
const double* xi,
const double mu_s,
+ const double C,
const double b,
const double rho_s,
const double d)
@@ -359,7 +378,7 @@ poisson_solver_1d_cell_update(int i,
coorp_term = dz*dz/(2.0*pow(xi[i], 2.0));
gi = idx1g(i);
g_out[gi] = 1.0/(1.0 + coorp_term)*(coorp_term*
- local_fluidity(p[i], mu[i], mu_s, b, rho_s, d)
+ local_fluidity(p[i], mu[i], mu_s, C, b, rho_s, d)
+ g_in[gi+1]/2.0
+ g_in[gi-1]/2.0);
@@ -409,6 +428,7 @@ implicit_1d_jacobian_poisson_solver(struct simulation* sim,
sim->sigma_n_eff,
sim->xi,
sim->mu_s,
+ sim->C,
sim->b,
sim->rho_s,
sim->d);
diff --git a/simulation.h b/simulation.h
@@ -33,6 +33,9 @@ struct simulation {
/* bulk and critical state static yield friction coefficient [-] */
double mu_s;
+ /* material cohesion [Pa] */
+ double C;
+
/* representative grain size [m] */
double d;
