cngf-pf

commit e31d3345807d0cb55f417a1955a601fa4ee8b14f
parent f36482758fdcf502bc8934b5a14aa13604d01d7b
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Wed, 18 Nov 2020 10:28:15 +0100

fix transient porosity computations

Diffstat:
M
cngf-pf.1
 | 
4
++--
M
cngf-pf.c
 | 
3
+--
M
simulation.c
 | 
105
+++++++++++++++++++++++++++++++++++++++++++++++--------------------------------
M
simulation.h
 | 
22
++++------------------

4 files changed, 70 insertions(+), 64 deletions(-)
diff --git a/cngf-pf.1 b/cngf-pf.1
@@ -108,8 +108,8 @@ Fluid dynamic viscosity [Pa*s] (1.787e-3).
 Only relevant with fluid dynamics enabled
 .Fl ( F ) .
 .It Fl K Ar dilatancy-constant
-Factor relating dilatancy and shear stress [TODO] (default 1.0).
-Only relevant with transient granular dynamics enabled
+Factor relating dilatancy angle to the volume fraction [-] (default
+4.09).  Only relevant with transient granular dynamics enabled
 .Fl ( T ) .
 .It Fl k Ar fluid-permeability
 Darcian intrinsic permeability of granular material [m^2] (default
diff --git a/cngf-pf.c b/cngf-pf.c
@@ -137,7 +137,7 @@ main(int argc, char *argv[])
 		sim.mu_f = atof(EARGF(usage()));
 		break;
 	case 'K':
-		sim.dilatancy_angle = atof(EARGF(usage()));
+		sim.dilatancy_constant = atof(EARGF(usage()));
 		break;
 	case 'k':
 		new_k = atof(EARGF(usage()));
@@ -280,7 +280,6 @@ main(int argc, char *argv[])
 			write_output_file(&sim, normalize);
 			filetimeclock = 0.0;
 		}
-		sim.t += sim.dt;
 		filetimeclock += sim.dt;
 		iter++;
 
diff --git a/simulation.c b/simulation.c
@@ -58,7 +58,7 @@ init_sim(struct simulation *sim)
 	sim->transient = 0;
 	sim->phi_min = 0.20;
 	sim->phi_max = 0.55;
-	sim->dilatancy_angle = 1.0;
+	sim->dilatancy_constant = 4.09;	/* Pailha & Pouliquen 2009 */
 
 	/* Iverson et al 1997, 1998: Storglaciaren till */
 	/* sim->mu_s = tan(DEG2RAD(26.3)); */
@@ -141,26 +141,27 @@ prepare_arrays(struct simulation *sim)
 	free(sim->phi);
 	free(sim->k);
 
-	sim->z = linspace(sim->origo_z,	/* spatial coordinates */
+	sim->z = linspace(sim->origo_z,	     /* spatial coordinates */
 			  sim->origo_z + sim->L_z,
 			  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->sigma_n_eff = zeros(sim->nz);	/* effective normal stress */
-	sim->sigma_n = zeros(sim->nz);	/* normal stess */
+	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->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 */
-	sim->phi = zeros(sim->nz);	/* porosity */
-	sim->phi_c = zeros(sim->nz);	/* critical-state porosity */
-	sim->phi_dot = zeros(sim->nz);	/* rate of porosity change */
-	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_local = zeros(sim->nz);	/* local fluidity */
-	sim->g_ghost = zeros(sim->nz + 2);	/* fluidity with ghost nodes */
-	sim->I = zeros(sim->nz);/* inertia number */
-	sim->tan_psi = zeros(sim->nz);	/* tan(dilatancy_angle) */
+	sim->p_f_dot = zeros(sim->nz);       /* fluid pressure change */
+	sim->phi = zeros(sim->nz);           /* porosity */
+	sim->phi_c = zeros(sim->nz);         /* critical-state porosity */
+	sim->phi_dot = zeros(sim->nz);       /* rate of porosity change */
+	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_local = zeros(sim->nz);       /* local fluidity */
+	sim->g_ghost = zeros(sim->nz + 2);   /* fluidity with ghost nodes */
+	sim->I = zeros(sim->nz);             /* inertia number */
+	sim->tan_psi = zeros(sim->nz);       /* tan(dilatancy_angle) */
 }
 
 void
@@ -172,6 +173,7 @@ free_arrays(struct simulation *sim)
 	free(sim->sigma_n_eff);
 	free(sim->sigma_n);
 	free(sim->p_f_ghost);
+	free(sim->p_f_dot);
 	free(sim->k);
 	free(sim->phi);
 	free(sim->phi_c);
@@ -317,10 +319,10 @@ check_simulation_parameters(struct simulation *sim)
 		warn_parameter_value("sim->phi_max is not within [0;1]",
 				     sim->phi_max, &return_status);
 
-	check_float("sim->dilatancy_angle", sim->dilatancy_angle, &return_status);
-	if (sim->dilatancy_angle < 0.0 || sim->dilatancy_angle > 100.0)
-		warn_parameter_value("sim->dilatancy_angle is not within [0;100]",
-				     sim->dilatancy_angle, &return_status);
+	check_float("sim->dilatancy_constant", sim->dilatancy_constant, &return_status);
+	if (sim->dilatancy_constant < 0.0 || sim->dilatancy_constant > 100.0)
+		warn_parameter_value("sim->dilatancy_constant is not within [0;100]",
+				     sim->dilatancy_constant, &return_status);
 
 	if (sim->fluid != 0 && sim->fluid != 1)
 		warn_parameter_value("sim->fluid has an invalid value",
@@ -407,7 +409,8 @@ 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];
+		sim->phi_c[i] = sim->phi_min
+		                + (sim->phi_max - sim->phi_min) * sim->I[i];
 }
 
 void
@@ -426,40 +429,38 @@ compute_critical_state_friction(struct simulation *sim)
 				 (sim->L_z - sim->z[i]) / sim->P_wall);
 }
 
-void
+static void
 compute_friction(struct simulation *sim)
 {
 	int i;
 
 	if (sim->transient)
 		for (i = 0; i < sim->nz; ++i)
-			sim->mu[i] = sim->tan_psi[i] + sim->mu_c[i];
+			sim->mu[i] = sim->mu_c[i] - sim->tan_psi[i];
 	else
 		for (i = 0; i < sim->nz; ++i)
 			sim->mu[i] = sim->mu_c[i];
 }
 
-void
+static void
 compute_tan_dilatancy_angle(struct simulation *sim)
 {
 	int i;
 
 	for (i = 0; i < sim->nz; ++i)
-		sim->tan_psi[i] = sim->dilatancy_angle * (sim->phi_c[i] - sim->phi[i]);
+		sim->tan_psi[i] = sim->dilatancy_constant * (sim->phi_c[i] - sim->phi[i]);
 }
 
-void
+static void
 compute_porosity_change(struct simulation *sim)
 {
 	int i;
 
-	for (i = 0; i < sim->nz; ++i) {
+	for (i = 0; i < sim->nz; ++i)
 		sim->phi_dot[i] = sim->tan_psi[i] * sim->gamma_dot_p[i] * sim->phi[i];
-		sim->phi[i] += sim->phi_dot[i] * sim->dt;
-	}
 }
 
-void
+static void
 compute_permeability(struct simulation *sim)
 {
 	int i;
@@ -478,7 +479,7 @@ shear_strain_rate_plastic(const double fluidity, const double friction)
 	return fluidity * friction;
 }
 
-void
+static void
 compute_shear_strain_rate_plastic(struct simulation *sim)
 {
 	int i;
@@ -488,7 +489,7 @@ compute_shear_strain_rate_plastic(struct simulation *sim)
 								sim->mu[i]);
 }
 
-void
+static void
 compute_shear_velocity(struct simulation *sim)
 {
 	int i;
@@ -525,7 +526,7 @@ cooperativity_length(const double A,
 	return A * d / sqrt(fabs((mu - C / p) - mu_s));
 }
 
-void
+static void
 compute_cooperativity_length(struct simulation *sim)
 {
 	int i;
@@ -554,7 +555,7 @@ local_fluidity(const double p,
 		return sqrt(p / rho_s * d * d) * ((mu - C / p) - mu_s) / (b * mu);
 }
 
-void
+static void
 compute_local_fluidity(struct simulation *sim)
 {
 	int i;
@@ -734,6 +735,22 @@ print_output(struct simulation *sim, FILE *fp, const int norm)
 	free(v_x_out);
 }
 
+static void
+temporal_increment(struct simulation *sim)
+{
+	int i;
+
+	if (sim->transient)
+		for (i = 0; i < sim->nz; ++i)
+			sim->phi[i] += sim->phi_dot[i] * sim->dt;
+
+	/*if (sim->fluid)
+		for (i = 1; i <= sim->nz; ++i)
+			sim->p_f_ghost[i] += sim->p_f_dot[i] * sim->dt;*/
+
+	sim->t += sim->dt;
+}
+
 int
 coupled_shear_solver(struct simulation *sim,
 		     const int max_iter,
@@ -753,7 +770,7 @@ coupled_shear_solver(struct simulation *sim,
 		do {		/* coupled iterations */
 
 			if (sim->transient) {
-				copy_values(sim->gamma_dot_p, oldval, sim->nz);
+				copy_values(sim->phi_dot, oldval, sim->nz);
 
 				/* step 1 */
 				compute_inertia_number(sim);	/* Eq. 1 */
@@ -791,26 +808,27 @@ coupled_shear_solver(struct simulation *sim,
 
 			/* step 9 */
 			compute_shear_strain_rate_plastic(sim);	/* Eq. 8 */
-			compute_inertia_number(sim);
+			compute_inertia_number(sim);	/* Eq. 1 */
 			compute_shear_velocity(sim);
 
 #ifdef DEBUG
-			for (i = 0; i < sim->nz; ++i) {
+			/* for (i = 0; i < sim->nz; ++i) { */
+			for (i = sim->nz-1; i < sim->nz; ++i) {
 				printf("\nsim->t = %g s\n", sim->t);
 				printf("sim->I[%d] = %g\n", i, sim->I[i]);
 				printf("sim->phi_c[%d] = %g\n", i, sim->phi_c[i]);
 				printf("sim->tan_psi[%d] = %g\n", i, sim->tan_psi[i]);
 				printf("sim->mu_c[%d] = %g\n", i, sim->mu_c[i]);
+				printf("sim->phi[%d] = %g\n", i, sim->phi[i]);
 				printf("sim->phi_dot[%d] = %g\n", i, sim->phi_dot[i]);
 				printf("sim->k[%d] = %g\n", i, sim->k[i]);
 				printf("sim->mu[%d] = %g\n", i, sim->mu[i]);
 			}
 #endif
-			/* stable velocity field == coupled solution
-			 * converged */
+			/* stable porosity change field == coupled solution converged */
 			if (sim->transient) {
 				for (i = 0; i < sim->nz; ++i)
-					r_norm[i] = fabs(residual(sim->gamma_dot_p[i], oldval[i]));
+					r_norm[i] = fabs(residual(sim->phi_dot[i], oldval[i]));
 				r_norm_max = max(r_norm, sim->nz);
 				if (r_norm_max <= rel_tol)
 					break;
@@ -858,6 +876,9 @@ coupled_shear_solver(struct simulation *sim,
 		free(r_norm);
 		free(oldval);
 	}
+
+	temporal_increment(sim);
+
 	return 0;
 }
 
diff --git a/simulation.h b/simulation.h
@@ -84,7 +84,7 @@ struct simulation {
 	double transient;
 	double phi_min;
 	double phi_max;
-	double dilatancy_angle;
+	double dilatancy_constant;
 
 	/* Fluid parameters */
 	int fluid;            /* flag to switch fluid on (1) or off (0) */
@@ -93,9 +93,9 @@ struct simulation {
 	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 p_f_mod_pulse_time; /* single pressure pulse at this time [s] */
-	int p_f_mod_pulse_shape; /* waveform for fluid-pressure pulse */
+	int p_f_mod_pulse_shape;   /* waveform for fluid-pressure pulse */
 	double beta_f;        /* adiabatic fluid compressibility [Pa^-1] */
-	double alpha;        /* adiabatic grain compressibility [Pa^-1] */
+	double alpha;         /* adiabatic grain compressibility [Pa^-1] */
 	double mu_f;          /* fluid dynamic viscosity [Pa*s] */
 	double rho_f;         /* fluid density [kg/m^3] */
 	double D;             /* diffusivity [m^2/s], overrides k, beta_f, alpha, mu_f */
@@ -106,6 +106,7 @@ struct simulation {
 	double *sigma_n_eff;  /* effective normal pressure [Pa] */
 	double *sigma_n;      /* normal stress [Pa] */
 	double *p_f_ghost;    /* fluid pressure [Pa] */
+	double *p_f_dot;      /* fluid pressure change [Pa/s] */
 	double *k;            /* hydraulic permeability [m^2] */
 	double *phi;          /* porosity [-] */
 	double *phi_c;        /* critical-state porosity [-] */
@@ -122,24 +123,9 @@ struct simulation {
 void init_sim(struct simulation *sim);
 void prepare_arrays(struct simulation *sim);
 void free_arrays(struct simulation *sim);
-
 void check_simulation_parameters(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_porosity_change(struct simulation *sim);
-void compute_permeability(struct simulation *sim);
-void compute_tan_dilatancy_angle(struct simulation *sim);
-
-void compute_cooperativity_length(struct simulation *sim);
-void compute_local_fluidity(struct simulation *sim);
-void compute_shear_strain_rate_plastic(struct simulation *sim);
-void compute_shear_velocity(struct simulation *sim);
 void compute_effective_stress(struct simulation *sim);
-void compute_friction(struct simulation *sim);
 
 void set_bc_neumann(double *a,
                     const int nz,