cngf-pf

commit 1d8e369516e373fbc0ed0bebd2dcbaad34d08e86
parent 7a0546b0c631a1241defb79e95c29cc3c7c91dce
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Mon,  6 Apr 2020 11:45:42 +0200

Fix calculation of friction, shear stress, and inertia number

Diffstat:
M
1d_fd_simple_shear.c
 | 
3
++-
M
simulation.c
 | 
12
+++++++-----
M
simulation.h
 | 
1
-

3 files changed, 9 insertions(+), 7 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -272,8 +272,8 @@ main(int argc, char* argv[])
 				compute_critical_state_porosity();	/* Eq. 2 */
 				/* step 3 */
 				compute_tan_dilatancy_angle();		/* Eq. 5 */
-				compute_critical_state_friction();	/* Eq. 7 */
 			}
+			compute_critical_state_friction();	/* Eq. 7 */
 
 			/* step 4 */
 			if (sim.transient) {
@@ -301,6 +301,7 @@ main(int argc, char* argv[])
 
 			/* step 9 */
 			compute_shear_strain_rate_plastic();	/* Eq. 8 */
+			compute_inertia_number();
 			compute_shear_velocity();
 
 			if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
diff --git a/simulation.c b/simulation.c
@@ -40,7 +40,6 @@ prepare_arrays()
 	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.tau = zeros(sim.nz);         /* shear stress */
 	sim.tan_psi = zeros(sim.nz);     /* tan(dilatancy_angle) */
 }
 
@@ -63,7 +62,6 @@ free_arrays()
 	free(sim.g_local);
 	free(sim.g_ghost);
 	free(sim.I);
-	free(sim.tau);
 	free(sim.tan_psi);
 }
 
@@ -304,8 +302,12 @@ void
 compute_friction()
 {
 	int i;
-	for (i=0; i<sim.nz; ++i)
-		sim.mu[i] = sim.tau[i]/sim.sigma_n_eff[i];
+	if (sim.transient)
+		for (i=0; i<sim.nz; ++i)
+			sim.mu[i] = sim.tan_psi[i] + sim.mu_c[i];
+	else
+		for (i=0; i<sim.nz; ++i)
+			sim.mu[i] = sim.mu_c[i];
 }
 
 void
@@ -612,7 +614,7 @@ print_output(FILE* fp, const int norm)
 		        sim.gamma_dot_p[i],
 				sim.phi[i],
 				sim.I[i],
-				sim.tau[i]);
+				sim.mu[i]*sim.sigma_n_eff[i]);
 
 	free(v_x_out);
 }
diff --git a/simulation.h b/simulation.h
@@ -112,7 +112,6 @@ struct simulation {
 	double* g_local;      /* local fluidity */
 	double* g_ghost;      /* fluidity with ghost nodes */
 	double* I;            /* inertia number [-] */
-	double* tau;          /* shear stress [Pa] */
 	double* tan_psi;      /* tan(dilatancy_angle) [-] */
 };