commit 593419284792a4d195e59bd8fda8c22c45873dbb
parent 3cabd0bed45a1f2f191cd53c3cc8ae807c34b304
Author: Anders Damsgaard <anders@adamsgaard.dk>
Date: Mon, 24 Jun 2019 13:12:34 +0200
Add fluid methods, use tabs for indentation
Diffstat:
1 file changed, 33 insertions(+), 26 deletions(-)
diff --git a/continuum-granular-manuscript1.tex b/continuum-granular-manuscript1.tex
@@ -62,40 +62,39 @@ We expand a steady-state continuum model for granular flow by \citet{Henann2013}
In our model, the sediment deforms as a highly nonlinear Bingham material with yield beyond the Mohr-Coulomb failure limit for a cohesionless material.
We assume that the elasticity is negligible and set the total shear rate $\dot{\gamma}$ to consist of a plastic contribution $\dot{\gamma}^\text{p}$:
\begin{equation}
- \dot{\gamma} \approx \dot{\gamma}^\text{p} = g(\mu_\text{c}, \sigma_\text{n}') \mu,
- \label{eq:shear-strain-rate}
+ \dot{\gamma} \approx \dot{\gamma}^\text{p} = g(\mu_\text{c}, \sigma_\text{n}') \mu,
+ \label{eq:shear-strain-rate}
\end{equation}
where $\mu = \tau/\sigma_\text{n}'$ is the dimensionless ratio between shear stress ($\tau$ [Pa]) and effective normal stress ($\sigma_\text{n}' = \sigma_\text{n} - p_f$ [Pa]).
Water pressure is $p_\text{f}$ [Pa] and $g$ [s$^{-1}$] is the granular fluidity. The fluidity consists of local and non-local components.
The local fluidity is defined as:
\begin{equation}
- g_\text{local}(\mu, \sigma_\text{n}') =
- \begin{cases}
- \sqrt{d^2 \sigma_\text{n}' / \rho_\text{s}} (\mu - \mu_\text{s})/(b\mu) &\text{if } \mu > \mu_\text{s} \text{, and}\\
- 0 &\text{if } \mu \leq \mu_\text{s}.
- \end{cases}
- \label{eq:g_local}
+ g_\text{local}(\mu, \sigma_\text{n}') =
+ \begin{cases}
+ \sqrt{d^2 \sigma_\text{n}' / \rho_\text{s}} (\mu - \mu_\text{s})/(b\mu) &\text{if } \mu > \mu_\text{s} \text{, and}\\
+ 0 &\text{if } \mu \leq \mu_\text{s}.
+ \end{cases}
+ \label{eq:g_local}
\end{equation}
For steady flow the non-locality is determined by the cooperativity length $\xi$:
\begin{equation}
- \nabla^2 g = \frac{1}{\xi^2(\mu)} (g - g_\text{local}),
- \label{eq:g}
+ \nabla^2 g = \frac{1}{\xi^2(\mu)} (g - g_\text{local}),
+ \label{eq:g}
\end{equation}
where
\begin{equation}
- \xi(\mu) = \frac{Ad}{\sqrt{|\mu - \mu_\text{s}|}}.
- \label{eq:cooperativity}
+ \xi(\mu) = \frac{Ad}{\sqrt{|\mu - \mu_\text{s}|}}.
+ \label{eq:cooperativity}
\end{equation}
-
Unlike \citet{Pailha2009} we do not implicitly prescribe the viscous drag during dilation and equation, and instead solve for the fluid pressure.
\subsection{Fluid-pressure evolution}%
\label{sub:fluid_pressure_evolution}
The transient evolution of pore-fluid pressure ($p_\text{f}$) is governed by Darcian pressure diffusion \citep[e.g.]{Goren2010, Goren2011, Damsgaard2017}:
\begin{equation}
- \frac{\partial p_\text{f}}{\partial t} = \frac{1}{\phi\mu_\text{f}\beta_\text{f}} \nabla \cdot (k \nabla p_\text{f})
- \label{eq:p_f}
+ \frac{\partial p_\text{f}}{\partial t} = \frac{1}{\phi\mu_\text{f}\beta_\text{f}} \nabla \cdot (k \nabla p_\text{f})
+ \label{eq:p_f}
\end{equation}
where $\mu_\text{f}$ denotes dynamic fluid viscosity [Pa s], $\beta_\text{f}$ is adiabatic fluid compressibility [Pa$^{-1}$], and $k$ is intrinsic permeability [m$^2$].
The sediment is assumed to be in the critical state throughout the domain, as in the original formulation by \citet{Henann2013}.
@@ -116,24 +115,32 @@ The fluidity field $g$ is solved for a set of mechanical forcings ($\mu$, $\sigm
We rearrange Eq.~\ref{eq:g} and split the Laplace operator ($\nabla^2$) into a 1D central finite difference 3-point stencil.
An iterative scheme is applied to relax the following equation at each grid node $i$:
\begin{equation}
- g_i = {\left(1 + \alpha_i\right)}^{-1}
- \left(\alpha_i g_\text{local}(\sigma_{\text{n},i}', \mu_i)
- + \frac{g_{i+1} + g_{i-1}}{2}
- \right)
- \label{eq:g_i}
+ g_i = {\left(1 + \alpha_i\right)}^{-1}
+ \left(\alpha_i g_\text{local}(\sigma_{\text{n},i}', \mu_i)
+ + \frac{g_{i+1} + g_{i-1}}{2}
+ \right)
+ \label{eq:g_i}
\end{equation}
where
\begin{equation}
- \alpha_i = \frac{\Delta z^2}{2\xi^2(\mu_i)}
- \label{eq:alpha}
+ \alpha_i = \frac{\Delta z^2}{2\xi^2(\mu_i)}
+ \label{eq:alpha}
\end{equation}
-Similarly, we also use operator splitting and finite differences to solve the equation for pore-pressure diffusion (Eq.~\ref{eq:p_f}).
The pore-pressure solution is constrained by a zero pressure gradient at the bottom ($dp_\text{f}/dz (z=0) = 0$), and a sinusoidal pressure forcing at the top ($p_\text{f}(z = L_z) = A \sin(2\pi f t) + p_{\text{f},0}$).
Here, $A$ is the forcing amplitude [Pa], $f$ is the forcing frequency [1/s], and $p_{\text{f},0}$ is the mean pore pressure over time [Pa].
-
-For each time step $\Delta t$, a pore-pressure solution is found by explicit temporal integration.
-We then use Jacobian iterations to find an implicit solution through underrelaxation over the same time step.
+We also use operator splitting and finite differences to solve the equation for pore-pressure diffusion (Eq.~\ref{eq:p_f}):
+\begin{equation}
+ \Delta p_{\text{f},i} = \frac{\Delta t}{\beta_\text{f} \phi_i \mu_\text{f}}
+ \frac{1}{\Delta z}
+ \left(
+ \frac{2 k_{i+1} k_i}{k_{i+1} + k_i} \frac{p_{i+1} - p_i}{\Delta z} -
+ \frac{2 k_{i-1} k_i}{k_{i-1} + k_i} \frac{p_i - p_{i-1}}{\Delta z} -
+ \right)
+ \label{eq:p_f_solution}
+\end{equation}
+For each time step $\Delta t$, a solution to Eq.~\ref{eq:p_f_solution} is first found by explicit temporal integration.
+We then use Jacobian iterations to find an implicit solution to the same equation through underrelaxation.
For the final pressure field at $t + \Delta t$ we mix the explicit and implicit solutions with equal weight, which is known as the Crank-Nicholson method \citep[e.g.][]{Patankar1980, Ferziger2002, Press2007}.
The method is unconditionally stable and second-order accurate in time and space.
