commit 2cb92cd64835afb9dc1a294082002be1d68be92b
parent 6391e0cc4a3206ad3e68fc78031e6dae4d9b9d47
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Thu, 9 Jul 2020 17:11:03 +0200
Add option to override individual hydraulic parameters with bulk diffusivity
Diffstat:
7 files changed, 75 insertions(+), 32 deletions(-)
diff --git a/1d_fd_simple_shear.1 b/1d_fd_simple_shear.1
@@ -14,6 +14,7 @@
.Op Fl b Ar grain-rate-dependence
.Op Fl C Ar fluid-compressibility
.Op Fl c Ar grain-cohesion
+.Op Fl D Ar fluid-diffusivity
.Op Fl d Ar grain-size
.Op Fl e Ar end-time
.Op Fl F
@@ -73,6 +74,11 @@ Only relevant with fluid dynamics enabled
.Fl ( F ) .
.It Fl c Ar grain-cohesion
Granular material cohesion [Pa] (default 0).
+.It Fl D Ar fluid-diffusivity
+Fluid diffusion coefficient [m^2/s] (default -1). Overrides fluid
+permeability (-k), grain compressibility (-P), fluid compressibility
+(-C), and fluid viscosity (-i). Do not use for transient simulations
+(-T). Disabled when set to a negative value.
.It Fl d Ar grain-size
Granular material representative grain size [m] (default 0.04).
.It Fl e Ar end-time
diff --git a/1d_fd_simple_shear.c b/1d_fd_simple_shear.c
@@ -28,6 +28,7 @@ usage(void)
"[-C fluid-compressibility] "
"[-c grain-cohesion] "
"[-d grain-size] "
+ "[-D fluid-diffusivity] "
"[-e end-time] "
"[-F] "
"[-f applied-shear-friction] "
@@ -108,6 +109,9 @@ main(int argc, char *argv[])
case 'd':
sim.d = atof(EARGF(usage()));
break;
+ case 'D':
+ sim.D = atof(EARGF(usage()));
+ break;
case 'e':
sim.t_end = atof(EARGF(usage()));
break;
diff --git a/fluid.c b/fluid.c
@@ -16,6 +16,14 @@ hydrostatic_fluid_pressure_distribution(struct simulation *sim)
* (sim->L_z - sim->z[i]);
}
+static double
+diffusivity(struct simulation *sim, int i) {
+ if (sim->D > 0.0)
+ return sim->D;
+ else
+ return sim->k[i]/((sim->alpha + sim->phi[i] * sim->beta_f) * sim->mu_f);
+}
+
/* Determines the largest time step for the current simulation state. Note
* that the time step should be recalculated if cell sizes or
* diffusivities (i.e., permeabilities, porosities, viscosities, or
@@ -43,7 +51,7 @@ set_largest_fluid_timestep(struct simulation *sim, const double safety)
diff_max = -INFINITY;
for (i = 0; i < sim->nz; ++i) {
- diff = sim->k[i] / ((sim->alpha + sim->phi[i] * sim->beta_f) * sim->mu_f);
+ diff = diffusivity(sim, i);
if (diff > diff_max)
diff_max = diff;
}
@@ -113,41 +121,49 @@ darcy_pressure_change_1d(const int i,
const double dz,
const double beta_f,
const double alpha,
- const double mu_f)
+ const double mu_f,
+ const double D)
{
- double k_ = k[i], div_k_grad_p, k_zn, k_zp;
-
- if (i == 0)
- k_zn = k_;
- else
- k_zn = k[i - 1];
- if (i == nz - 1)
- k_zp = k_;
- else
- k_zp = k[i + 1];
+ double k_, div_k_grad_p, k_zn, k_zp;
+
+ if (D > 0.0)
+ return D*(p_f_ghost_in[i + 2] - 2.0 * p_f_ghost_in[i + 1] - p_f_ghost_in[i])/(dz*dz);
+
+ else {
+
+ k_ = k[i];
+ if (i == 0)
+ k_zn = k_;
+ else
+ k_zn = k[i - 1];
+ if (i == nz - 1)
+ k_zp = k_;
+ else
+ k_zp = k[i + 1];
#ifdef DEBUG
- printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n",
- i, i + 1,
- p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
- k_zn, k_, k_zp);
+ printf("%d->%d: p=[%g, %g, %g]\tk=[%g, %g, %g]\n",
+ i, i + 1,
+ p_f_ghost_in[i], p_f_ghost_in[i + 1], p_f_ghost_in[i + 2],
+ k_zn, k_, k_zp);
#endif
- div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
- * (p_f_ghost_in[i + 2]
- - p_f_ghost_in[i + 1]) / dz
- - 2.0 * k_zn * k_ / (k_zn + k_)
- * (p_f_ghost_in[i + 1]
- - p_f_ghost_in[i]) / dz
- ) / dz;
+ div_k_grad_p = (2.0 * k_zp * k_ / (k_zp + k_)
+ * (p_f_ghost_in[i + 2]
+ - p_f_ghost_in[i + 1]) / dz
+ - 2.0 * k_zn * k_ / (k_zn + k_)
+ * (p_f_ghost_in[i + 1]
+ - p_f_ghost_in[i]) / dz
+ ) / dz;
#ifdef DEBUG
- printf("darcy_pressure_change [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
- i, phi[i], div_k_grad_p, phi_dot[i]);
+ printf("darcy_pressure_change [%d]: phi=%g\tdiv_k_grad_p=%g\tphi_dot=%g\n",
+ i, phi[i], div_k_grad_p, phi_dot[i]);
#endif
- /* TODO: add advective term */
- return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
- - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
+ /* TODO: add advective term */
+ return 1.0 / ((alpha + beta_f * phi[i]) * mu_f) * div_k_grad_p
+ - 1.0 / ((alpha + beta_f * phi[i]) * (1.0 - phi[i])) * phi_dot[i];
+ }
}
int
@@ -207,7 +223,8 @@ darcy_solver_1d(struct simulation *sim,
sim->dz,
sim->beta_f,
sim->alpha,
- sim->mu_f);
+ sim->mu_f,
+ sim->D);
}
if (epsilon > 0.0) {
/* compute implicit solution with Jacobian iterations */
@@ -236,7 +253,8 @@ darcy_solver_1d(struct simulation *sim,
sim->dz,
sim->beta_f,
sim->alpha,
- sim->mu_f);
+ sim->mu_f,
+ sim->D);
for (i = 0; i < sim->nz - 1; ++i) {
#ifdef DEBUG
diff --git a/max_depth_simple_shear.1 b/max_depth_simple_shear.1
@@ -8,6 +8,7 @@
.Nm
.Op Fl a Ar fluid-pressure-ampl
.Op Fl C Ar fluid-compressibility
+.Op Fl D Ar fluid-diffusivity
.Op Fl g Ar gravity-accel
.Op Fl h
.Op Fl i Ar fluid-viscosity
@@ -35,6 +36,11 @@ and are as follows:
Amplitude of fluid-pressure perturbations [Pa] (default 0).
.It Fl C Ar fluid-compressibility
Fluid adiabatic compressibility [Pa^-1] (default 3.9e-10).
+.It Fl D Ar fluid-diffusivity
+Fluid diffusion coefficient [m^2/s] (default -1). Overrides fluid
+permeability (-k), grain compressibility (-P), fluid compressibility
+(-C), and fluid viscosity (-i). Disabled when set to a negative
+value.
.It Fl g Ar gravity-accel
Gravity magnitude [m/s^2] (default 9.81).
.It Fl h
diff --git a/max_depth_simple_shear.c b/max_depth_simple_shear.c
@@ -29,6 +29,7 @@ usage(void)
fprintf(stderr, "usage: %s "
"[-a fluid-pressure-ampl] "
"[-C fluid-compressibility] "
+ "[-D fluid-diffusivity] "
"[-g gravity-accel] "
"[-h] "
"[-i fluid-viscosity] "
@@ -47,8 +48,10 @@ usage(void)
static double
skin_depth(const struct simulation *sim)
{
- return sqrt(sim->k[0] /
- (sim->mu_f * (sim->alpha + sim->phi[0] * sim->beta_f) * M_PI * sim->p_f_mod_freq));
+ if (sim->D > 0.0)
+ return sqrt(sim->D * M_PI * sim->p_f_mod_freq);
+ else
+ return sqrt(sim->k[0] / (sim->mu_f * (sim->alpha + sim->phi[0] * sim->beta_f) * M_PI * sim->p_f_mod_freq));
}
/* using alternate form: sin(x) + cos(x) = sqrt(2)*sin(x + pi/4) */
@@ -186,6 +189,9 @@ main(int argc, char *argv[])
case 'C':
sim.beta_f = atof(EARGF(usage()));
break;
+ case 'D':
+ sim.D = atof(EARGF(usage()));
+ break;
case 'g':
sim.G = atof(EARGF(usage()));
break;
diff --git a/simulation.c b/simulation.c
@@ -107,6 +107,8 @@ init_sim(struct simulation *sim)
sim->alpha = 1e-8;
+ sim->D = -1.0; /* disabled when negative */
+
/* Damsgaard et al 2015 */
sim->k = initval(1.9e-15, 1);
diff --git a/simulation.h b/simulation.h
@@ -98,6 +98,7 @@ struct simulation {
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 */
/* arrays */
double *mu; /* static yield friction [-] */
