commit 5e6b491c007eaf73cc937101f242f9c19227c277
parent c49326e4c9c1b8d919131befc31cb0cf2a48ef86
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 19 Aug 2019 17:35:54 +0200
Implement rate-controlled and rate-limited shear, bump version to 0.4.0
Diffstat:
| M | README.md | | | 14 | ++++++++++---- |
| M | main.c | | | 75 | +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------- |
| M | parameter_defaults.h | | | 2 | ++ |
| M | simulation.h | | | 146 | +++++++++++++++++++++++++++++++++++++++++-------------------------------------- |
4 files changed, 148 insertions(+), 89 deletions(-)
diff --git a/README.md b/README.md
@@ -23,11 +23,16 @@ Optional arguments:
-N, --normalize normalize output velocity
-G, --gravity VAL gravity magnitude [m/s^2] (default 9.81)
-P, --normal-stress VAL normal stress on top [Pa] (default 120000)
- -m, --stress-ratio VAL applied stress ratio [-] (default 0.4)
+ -m, --stress-ratio VAL applied stress ratio [-] (default 0.45)
+ -s, --set-shear-velocity VAL set top shear velocity to this value [m/s]
+ (default nan), overrides --stress-ratio
+ -l, --limit-shear-velocity VAL limit top shear velocity to this value [m/s]
+ (default nan), overrides --stress-ratio and
+ --set-shear-velocity
-V, --velocity-bottom VAL base velocity at bottom [m/s] (default 0)
-A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default 0.4)
-b, --rate-dependence VAL rate dependence beyond yield [-] (default 0.9377)
- -f, --friction-coefficient VAL grain friction coefficient [-] (default 0.366614)
+ -f, --friction-coefficient VAL grain friction coefficient [-] (default 0.404026)
-C, --cohesion VAL grain cohesion [Pa] (default 0)
-p, --porosity VAL porosity fraction [-] (default 0.25)
-d, --grain-size VAL representative grain size [m] (default 0.04)
@@ -38,8 +43,8 @@ Optional arguments:
Optional arguments only relevant with transient (fluid) simulation:
-F, --fluid enable pore fluid computations
- -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default 4.5e-10)
- -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default 0.001)
+ -c, --fluid-compressibility VAL fluid compressibility [Pa^-1] (default 3.9e-10)
+ -i, --fluid-viscosity VAL fluid viscosity [Pa*s] (default 0.001787)
-R, --fluid-density VAL fluid density [kg/m^3] (default 1000)
-k, --fluid-permeability VAL fluid intrinsic permeability [m^2] (default 1.9e-15)
-O, --fluid-pressure-top VAL fluid pressure at +z edge [Pa] (default 0)
@@ -49,6 +54,7 @@ Optional arguments only relevant with transient (fluid) simulation:
-t, --time VAL simulation start time [s] (default 0)
-T, --time-end VAL simulation end time [s] (default 1)
-I, --file-interval VAL interval between output files [s] (default 0.1)
+
```
## Results
diff --git a/main.c b/main.c
@@ -7,11 +7,14 @@
#include "simulation.h"
#include "fluid.h"
-#define VERSION "0.3.1"
+#define VERSION "0.4.0"
#define PROGNAME "1d_fd_simple_shear"
#include "parameter_defaults.h"
+/* relative tolerance criteria when shear velocity is restricted */
+#define RTOL 1e-3
+
static void
usage(void)
{
@@ -26,6 +29,11 @@ usage(void)
" -G, --gravity VAL gravity magnitude [m/s^2] (default %g)\n"
" -P, --normal-stress VAL normal stress on top [Pa] (default %g)\n"
" -m, --stress-ratio VAL applied stress ratio [-] (default %g)\n"
+ " -s, --set-shear-velocity VAL set top shear velocity to this value [m/s]\n"
+ " (default %g), overrides --stress-ratio\n"
+ " -l, --limit-shear-velocity VAL limit top shear velocity to this value [m/s]\n"
+ " (default %g), overrides --stress-ratio and\n"
+ " --set-shear-velocity\n"
" -V, --velocity-bottom VAL base velocity at bottom [m/s] (default %g)\n"
" -A, --nonlocal-amplitude VAL amplitude of nonlocality [-] (default %g)\n"
" -b, --rate-dependence VAL rate dependence beyond yield [-] (default %g)\n"
@@ -54,6 +62,8 @@ usage(void)
sim.G,
sim.P_wall,
sim.mu_wall,
+ sim.v_x_fix,
+ sim.v_x_limit,
sim.v_x_bot,
sim.A,
sim.b,
@@ -96,15 +106,15 @@ main(int argc, char* argv[])
int norm, opt, i;
struct simulation sim;
const char* optstring;
- unsigned long iter;
- double new_phi, new_k, filetimeclock;
+ unsigned long iter, stressiter;
+ double new_phi, new_k, filetimeclock, res_norm;
/* load with default values */
sim = init_sim();
norm = 0;
- 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:";
+ optstring = "hvNn:G:P:m:s:l: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'},
@@ -112,6 +122,8 @@ main(int argc, char* argv[])
{"gravity", required_argument, NULL, 'G'},
{"normal-stress", required_argument, NULL, 'P'},
{"stress-ratio", required_argument, NULL, 'm'},
+ {"set-shear-velocity", required_argument, NULL, 's'},
+ {"limit-shear-velocity", required_argument, NULL, 'l'},
{"velocity-bottom", required_argument, NULL, 'V'},
{"nonlocal-amplitude", required_argument, NULL, 'A'},
{"rate-dependence", required_argument, NULL, 'b'},
@@ -171,6 +183,12 @@ main(int argc, char* argv[])
case 'm':
sim.mu_wall = atof(optarg);
break;
+ case 's':
+ sim.v_x_fix = atof(optarg);
+ break;
+ case 'l':
+ sim.v_x_limit = atof(optarg);
+ break;
case 'V':
sim.v_x_bot = atof(optarg);
break;
@@ -273,22 +291,51 @@ main(int argc, char* argv[])
filetimeclock = 0.0;
iter = 0;
+ stressiter = 0;
while (sim.t <= sim.t_end) {
- if (sim.fluid) {
- if (darcy_solver_1d(&sim, 10000, 1e-5))
+ do {
+ if (sim.fluid) {
+ if (darcy_solver_1d(&sim, 10000, 1e-5))
+ exit(1);
+ }
+
+ compute_effective_stress(&sim);
+ compute_friction(&sim);
+ compute_cooperativity_length(&sim);
+
+ if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5))
exit(1);
- }
- compute_effective_stress(&sim);
- compute_friction(&sim);
- compute_cooperativity_length(&sim);
+ compute_shear_strain_rate_plastic(&sim);
+ compute_shear_velocity(&sim);
+
+ if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
+ if (!isnan(sim.v_x_limit)) {
+ res_norm = (sim.v_x_limit - sim.v_x[sim.nz-1])
+ /sim.v_x[sim.nz-1];
+ if (res_norm > 0.0)
+ res_norm = 0.0;
+ } else {
+ res_norm = (sim.v_x_fix - sim.v_x[sim.nz-1])
+ /sim.v_x[sim.nz-1];
+ }
+ sim.mu_wall *= 1.0 + (res_norm*1e-2);
+ }
- if (implicit_1d_jacobian_poisson_solver(&sim, 10000, 1e-5))
- exit(1);
+ if (++stressiter > 10000) {
+ fprintf(stderr, "error: stress solution did not converge:\n");
+ fprintf(stderr,
+ "v_x=%g, v_x_fix=%g, v_x_limit=%g, "
+ "res_norm=%g, mu_wall=%g\n",
+ sim.v_x[sim.nz-1], sim.v_x_fix, sim.v_x_limit,
+ res_norm, sim.mu_wall);
+ free_arrays(&sim);
+ return 10;
+ }
- compute_shear_strain_rate_plastic(&sim);
- compute_shear_velocity(&sim);
+ } while ((!isnan(sim.v_x_fix) || !isnan(sim.v_x_limit))
+ && fabs(res_norm) > RTOL);
sim.t += sim.dt;
filetimeclock += sim.dt;
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -20,6 +20,8 @@ struct simulation init_sim(void)
sim.P_wall = 120e3; /* larger normal stress deepens the shear depth */
sim.mu_wall = 0.45;
sim.v_x_bot = 0.0;
+ sim.v_x_fix = NAN;
+ sim.v_x_limit = NAN;
sim.nz = 100;
diff --git a/simulation.h b/simulation.h
@@ -9,90 +9,96 @@
/* Simulation settings */
struct simulation {
- /* simulation name to use for output files */
- char name[100];
+ /* simulation name to use for output files */
+ char name[100];
- /* gravitational acceleration magnitude [m/s^2] */
- double G;
+ /* gravitational acceleration magnitude [m/s^2] */
+ double G;
- /* wall parameters */
- double P_wall; // normal stress from top wall [Pa]
+ /* wall parameters */
+ double P_wall; // normal stress from top wall [Pa]
- /* bottom velocity along x [m/s] */
- double v_x_bot;
+ /* optionally fix top shear velocity to this value [m/s] */
+ double v_x_fix;
- /* stress ratio at top wall */
- double mu_wall;
+ /* optionally fix top shear velocity to this value [m/s] */
+ double v_x_limit;
- /* nonlocal amplitude [-] */
- double A;
+ /* bottom velocity along x [m/s] */
+ double v_x_bot;
- /* rate dependence beyond yield [-] */
- double b;
+ /* stress ratio at top wall */
+ double mu_wall;
- /* bulk and critical state static yield friction coefficient [-] */
- double mu_s;
+ /* nonlocal amplitude [-] */
+ double A;
- /* material cohesion [Pa] */
- double C;
+ /* rate dependence beyond yield [-] */
+ double b;
- /* representative grain size [m] */
- double d;
+ /* bulk and critical state static yield friction coefficient [-] */
+ double mu_s;
- /* grain material density [kg/m^3] */
- double rho_s;
+ /* material cohesion [Pa] */
+ double C;
- /* nodes along z */
- int nz;
+ /* representative grain size [m] */
+ double d;
- /* origo of axis [m] */
- double origo_z;
+ /* grain material density [kg/m^3] */
+ double rho_s;
- /* length of domain [m] */
- double L_z;
+ /* nodes along z */
+ int nz;
- /* array of cell coordinates */
- double* z;
+ /* origo of axis [m] */
+ double origo_z;
- /* cell spacing [m] */
- double dz;
+ /* length of domain [m] */
+ double L_z;
- /* current time [s] */
- double t;
+ /* array of cell coordinates */
+ double* z;
- /* end time [s] */
- double t_end;
+ /* cell spacing [m] */
+ double dz;
- /* time step length [s] */
- double dt;
+ /* current time [s] */
+ double t;
- /* interval between output files [s] */
- double file_dt;
+ /* end time [s] */
+ double t_end;
- /* output file number */
- int n_file;
+ /* time step length [s] */
+ double dt;
- /* Fluid parameters */
- int fluid; /* flag to switch fluid on (1) or off (0) */
- double p_f_top; /* fluid pressure at the top [Pa] */
- double p_f_mod_ampl; /* amplitude of fluid pressure variations [Pa] */
- double p_f_mod_freq; /* frequency of fluid pressure variations [s^-1] */
- double p_f_mod_phase; /* phase of fluid pressure variations [s^-1] */
- double beta_f; /* adiabatic fluid compressibility [Pa^-1] */
- double mu_f; /* fluid dynamic viscosity [Pa*s] */
- double rho_f; /* fluid density [kg/m^3] */
+ /* interval between output files [s] */
+ double file_dt;
- /* arrays */
- double* mu; /* static yield friction [-] */
- double* sigma_n_eff; /* effective normal pressure [Pa] */
- double* sigma_n; /* normal stress [Pa] */
- double* p_f_ghost; /* fluid pressure [Pa] */
- double* k; /* hydraulic permeability [m^2] */
- double* phi; /* porosity [-] */
- double* xi; /* cooperativity length */
- double* gamma_dot_p; /* plastic shear strain rate [1/s] */
- double* v_x; /* shear velocity [m/s] */
- double* g_ghost; /* fluidity with ghost nodes */
+ /* output file number */
+ int n_file;
+
+ /* Fluid parameters */
+ int fluid; /* flag to switch fluid on (1) or off (0) */
+ double p_f_top; /* fluid pressure at the top [Pa] */
+ double p_f_mod_ampl; /* amplitude of fluid pressure variations [Pa] */
+ double p_f_mod_freq; /* frequency of fluid pressure variations [s^-1] */
+ double p_f_mod_phase; /* phase of fluid pressure variations [s^-1] */
+ double beta_f; /* adiabatic fluid compressibility [Pa^-1] */
+ double mu_f; /* fluid dynamic viscosity [Pa*s] */
+ double rho_f; /* fluid density [kg/m^3] */
+
+ /* arrays */
+ double* mu; /* static yield friction [-] */
+ double* sigma_n_eff; /* effective normal pressure [Pa] */
+ double* sigma_n; /* normal stress [Pa] */
+ double* p_f_ghost; /* fluid pressure [Pa] */
+ double* k; /* hydraulic permeability [m^2] */
+ double* phi; /* porosity [-] */
+ double* xi; /* cooperativity length */
+ double* gamma_dot_p; /* plastic shear strain rate [1/s] */
+ double* v_x; /* shear velocity [m/s] */
+ double* g_ghost; /* fluidity with ghost nodes */
};
@@ -109,11 +115,10 @@ void set_bc_neumann(double* g_ghost,
const double df,
const double dx);
-void set_bc_dirichlet(
- double* g_ghost,
- const int nz,
- const int boundary,
- const double value);
+void set_bc_dirichlet(double* g_ghost,
+ const int nz,
+ const int boundary,
+ const double value);
void compute_cooperativity_length(struct simulation* sim);
void compute_shear_strain_rate_plastic(struct simulation* sim);
@@ -121,10 +126,9 @@ void compute_shear_velocity(struct simulation* sim);
void compute_effective_stress(struct simulation* sim);
void compute_friction(struct simulation* sim);
-int implicit_1d_jacobian_poisson_solver(
- struct simulation* sim,
- const int max_iter,
- const double rel_tol);
+int implicit_1d_jacobian_poisson_solver(struct simulation* sim,
+ const int max_iter,
+ const double rel_tol);
void write_output_file(struct simulation* sim, const int normalize);
void print_dry_output(FILE* fp, struct simulation* sim, const int normalize);
