cngf-pf

commit d808b3f1328fbb79c15780d77112b676320c51d5
parent 0bc432ae20a857af1cc3b8dfd9468e0f8b9154fc
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date:   Thu,  2 Apr 2020 16:50:05 +0200

Introduce functions for transient dynamics

Diffstat:
M
1d_fd_simple_shear.c
 | 
9
+++++++++
M
simulation.c
 | 
83
+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

2 files changed, 92 insertions(+), 0 deletions(-)
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -268,6 +268,13 @@ main(int argc, char* argv[])
 					exit(1);
 
 			compute_effective_stress();
+			if (sim.transient) {
+				compute_inertia_number();               /* new */
+				compute_critical_state_porosity();      /* new */
+				compute_tan_dilatancy_angle();          /* new */
+				compute_critical_state_shear_stress();  /* new */
+				compute_shear_stress();                 /* new */
+			}
 			compute_friction();
 			compute_cooperativity_length();
 
@@ -277,6 +284,8 @@ main(int argc, char* argv[])
 
 			compute_shear_strain_rate_plastic();
 			compute_shear_velocity();
+			if (sim.transient)
+				compute_porosity_change();              /* new */
 
 			if (!isnan(sim.v_x_limit) || !isnan(sim.v_x_fix)) {
 				if (!isnan(sim.v_x_limit)) {
diff --git a/simulation.c b/simulation.c
@@ -252,6 +252,89 @@ lithostatic_pressure_distribution()
 		                  (sim.L_z - sim.z[i]);
 }
 
+/* NEW FUNCTIONS START */
+
+void
+compute_inetria_number()
+{
+	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()
+{
+	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()
+{
+	int i;
+	if (sim.fluid)
+		for (i=0; i<sim.nz; ++i)
+			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_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()
+{
+	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()
+{
+	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()
+{
+	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()
+{
+	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()
+{
+	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()
+{
+	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 */
+
 void
 compute_friction()
 {