1d_fd_simple_shear_transient

commit c0199d862f6b626602fe7b904b838fa6354d4672
parent d94d9afa17f8d1ad7bc2e903aa9e1f487bd7eb2e
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Fri,  8 Nov 2019 14:23:11 +0100

Begin implementing transient solution

Diffstat:
M
Makefile
 | 
2
+-
M
main.c
 | 
12
++++++++++--
M
parameter_defaults.h
 | 
4
++++
M
simulation.c
 | 
87
++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---
M
simulation.h
 | 
24
++++++++++++++++++++++++

5 files changed, 123 insertions(+), 6 deletions(-)
diff --git a/Makefile b/Makefile
@@ -3,7 +3,7 @@ LDFLAGS = -lm
 SRC = arrays.c fluid.c main.c simulation.c
 OBJ = $(SRC:.c=.o)
 HDR = arrays.h fluid.h parameter_defaults.h simulation.h
-BIN = ./1d_fd_simple_shear
+BIN = ./1d_fd_simple_shear_transient
 
 PREFIX ?= /usr/local
 INSTALL ?= install
diff --git a/main.c b/main.c
@@ -7,8 +7,8 @@
 #include "simulation.h"
 #include "fluid.h"
 
-#define VERSION "0.4.6"
-#define PROGNAME "1d_fd_simple_shear"
+#define VERSION "0.1.0"
+#define PROGNAME "1d_fd_simple_shear_transient"
 
 #include "parameter_defaults.h"
 
@@ -347,12 +347,19 @@ main(int argc, char* argv[])
 
 		stressiter = 0;
 		do {
+			compute_permeability(&sim);
+
 			if (sim.fluid) {
 				if (darcy_solver_1d(&sim, MAX_ITER_DARCY, RTOL_DARCY))
 					exit(1);
 			}
 
 			compute_effective_stress(&sim);
+			compute_inertia_number(&sim);
+			compute_critical_state_porosity(&sim);
+			compute_tan_dilatancy_angle(&sim);
+			compute_critical_state_shear_stress(&sim);
+			compute_shear_stress(&sim);
 			compute_friction(&sim);
 			compute_cooperativity_length(&sim);
 
@@ -362,6 +369,7 @@ main(int argc, char* argv[])
 
 			compute_shear_strain_rate_plastic(&sim);
 			compute_shear_velocity(&sim);
+			compute_porosity_change(&sim);
 
 			if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
 				if (!isnan(sim.v_x_limit)) {
diff --git a/parameter_defaults.h b/parameter_defaults.h
@@ -43,6 +43,10 @@ struct simulation init_sim(void)
 	sim.phi = initval(0.25, 1);
 	sim.d = 0.04;        /* Damsgaard et al 2013 */
 
+	sim.phi_min = 0.2;
+	sim.phi_max = 0.4;
+	sim.K = 4.09;        /* Pailha and Pouliquen 2009 */
+
 	/* Iverson et al 1997, 1998: Storglaciaren till */
 	/* sim.mu_s = tan(DEG2RAD(26.3)); */
 	/* sim.C = 5.0e3; */
diff --git a/simulation.c b/simulation.c
@@ -14,17 +14,23 @@ prepare_arrays(struct simulation* sim)
 	                  sim->nz);
 	sim->dz = sim->z[1] - sim->z[0];   /* cell spacing */
 	sim->mu = zeros(sim->nz);          /* stress ratio */
+	sim->mu_c = zeros(sim->nz);        /* critical-state stress ratio */
+	sim->tau = zeros(sim->nz);         /* shear stress */
+	sim->tau_c = zeros(sim->nz);       /* critical-state shear stress */
 	sim->sigma_n_eff = zeros(sim->nz); /* effective normal stress */
 	sim->sigma_n = zeros(sim->nz);     /* normal stess */
 	sim->p_f_ghost = zeros(sim->nz+2); /* fluid pressure with ghost nodes */
 	free(sim->phi);
 	sim->phi = zeros(sim->nz);         /* porosity */
+	sim->phi_c = zeros(sim->nz);       /* critical-state porosity */
 	free(sim->k);
 	sim->k = zeros(sim->nz);           /* permeability */
 	sim->xi = zeros(sim->nz);          /* cooperativity length */
 	sim->gamma_dot_p = zeros(sim->nz); /* shear velocity */
 	sim->v_x = zeros(sim->nz);         /* shear velocity */
 	sim->g_ghost = zeros(sim->nz+2);   /* fluidity with ghost nodes */
+	sim->I = zeros(sim->nz);           /* inertia number */
+	sim->tan_dilatancy_angle = zeros(sim->nz); /* tan(dilatancy_angle) */
 }
 
 void
@@ -32,15 +38,21 @@ free_arrays(struct simulation* sim)
 {
 	free(sim->z);
 	free(sim->mu);
+	free(sim->mu_c);
+	free(sim->tau);
+	free(sim->tau_c);
 	free(sim->sigma_n_eff);
 	free(sim->sigma_n);
 	free(sim->p_f_ghost);
 	free(sim->k);
 	free(sim->phi);
+	free(sim->phi_c);
 	free(sim->xi);
 	free(sim->gamma_dot_p);
 	free(sim->v_x);
 	free(sim->g_ghost);
+	free(sim->I);
+	free(sim->tan_dilatancy_angle);
 }
 
 static void
@@ -227,21 +239,90 @@ lithostatic_pressure_distribution(struct simulation* sim)
 		                  (sim->L_z - sim->z[i]);
 }
 
+/* NEW FUNCTIONS START */
+
 void
-compute_friction(struct simulation* sim)
+compute_inertia_number(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->I[i] = sim->gamma_dot_p[i]*sim->d
+				    /sqrt(sim->sigma_n_eff[i]/sim->rho_s);
+}
+
+void
+compute_critical_state_porosity(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->phi_c[i] = sim->phi_min + (sim->phi_max - sim->phi_min)*sim->I[i];
+}
+
+void
+compute_critical_state_friction(struct simulation* sim)
 {
 	int i;
 	if (sim->fluid)
 		for (i=0; i<sim->nz; ++i)
-			sim->mu[i] = sim->mu_wall/
+			sim->mu_c[i] = sim->mu_wall/
 			             (sim->sigma_n_eff[i]/(sim->P_wall - sim->p_f_top));
 	else
 		for (i=0; i<sim->nz; ++i)
-			sim->mu[i] = sim->mu_wall/
+			sim->mu_c[i] = sim->mu_wall/
 			             (1.0 + (1.0 - sim->phi[i])*sim->rho_s*sim->G*
 			             (sim->L_z - sim->z[i])/sim->P_wall);
 }
 
+void
+compute_friction(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->mu[i] = sim->tau[i]/sim->sigma_n_eff[i];
+}
+
+void
+compute_critical_state_shear_stress(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->tau_c[i] = sim->mu_c[i]*sim->sigma_n_eff[i];
+}
+
+void
+compute_shear_stress(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->tau[i] = sim->tan_dilatancy_angle[i]*sim->sigma_n_eff[i] + sim->tau_c[i];
+}
+
+void
+compute_tan_dilatancy_angle(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->tan_dilatancy_angle[i] = sim->K*(sim->phi_c[i] - sim->phi[i]);
+}
+
+void
+compute_porosity_change(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->phi[i] += sim->tan_dilatancy_angle[i]*sim->gamma_dot_p[i]*sim->phi[i]*sim->dt;
+}
+
+void
+compute_permeability(struct simulation* sim)
+{
+	int i;
+	for (i=0; i<sim->nz; ++i)
+		sim->k[i] = pow(sim->d, 2.0)/180.0 * pow(sim->phi[i], 3.0)/pow(1.0 - sim->phi[i], 2.0);
+}
+
+/* NEW FUNCTIONS END */
+
 double
 shear_strain_rate_plastic(const double fluidity, const double friction)
 {
diff --git a/simulation.h b/simulation.h
@@ -48,6 +48,15 @@ struct simulation {
 	/* grain material density [kg/m^3] */
 	double rho_s;
 
+	/* smallest observable porosity [-] (I -> 0) */
+	double phi_min;
+
+	/* porosity when I=0 [-] */
+	double phi_max;
+
+	/* dilatancy constant [-] */
+	double K;
+
 	/* nodes along z */
 	int nz;
 
@@ -92,15 +101,21 @@ struct simulation {
 
 	/* arrays */
 	double* mu;           /* static yield friction [-] */
+	double* mu_c;         /* critical-state static yield friction [-] */
+	double* tau;          /* shear stess [Pa] */
+	double* tau_c;        /* critical-state shear stress [Pa] */
 	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* phi_c;        /* critical-state 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 */
+	double* I;            /* inertia number [-] */
+	double* tan_dilatancy_angle; /* tan(dilatancy_angle) [-] */
 };
 
 
@@ -111,6 +126,15 @@ void check_simulation_parameters(const struct simulation* sim);
 
 void lithostatic_pressure_distribution(struct simulation* sim);
 
+void compute_inertia_number(struct simulation* sim);
+void compute_critical_state_porosity(struct simulation* sim);
+void compute_critical_state_friction(struct simulation* sim);
+void compute_critical_state_shear_stress(struct simulation* sim);
+void compute_shear_stress(struct simulation* sim);
+void compute_tan_dilatancy_angle(struct simulation* sim);
+void compute_porosity_change(struct simulation* sim);
+void compute_permeability(struct simulation* sim);
+
 void set_bc_neumann(double* g_ghost,
                     const int nz,
                     const int boundary,