Coupled Geomechanics and Single Phase Flow
Mathematical Formulation
The governing equation for mechanical deformation of the solid-fluid system is
\[\boxed{\mathrm{div} \sigma + \rho_b g = 0}\]
where $\mathrm{div}$ is the divergence operator, $\sigma$ is the Cauchy total-stress
\[\sigma = \begin{bmatrix} \sigma_{xx} & \sigma_{xy}\\ \sigma_{xy} & \sigma_{yy} \end{bmatrix}\]
and $g\in\mathbb{R}^2$ is the gravity vector, $\rho_b=\phi \rho_f + (1-\phi)\rho_s$ is the bulk density, $\rho_f$ is total fluid density, $\rho_s$ is the density of the solid phase and $\phi$ is the true porosity.
The stress-strain relation for linear poroelasticity takes the form
\[\sigma = \sigma' - bp\mathbf{I},\quad \sigma' = \begin{bmatrix} \sigma'_{xx} & \sigma'_{xy}\\ \sigma'_{xy} & \sigma'_{yy} \end{bmatrix}, \quad\begin{bmatrix} \delta\sigma'_{xx}\\\delta\sigma'_{yy}\\\delta\sigma'_{xy} \end{bmatrix} = H\begin{bmatrix} \delta\varepsilon_{xx}\\\delta\varepsilon_{yy}\\2\delta\varepsilon_{xy} \end{bmatrix}\]
where $\mathbf{I}$ is the identity matrix, $p$ is the pressure, $b$ is the Biot coefficient, $D$ is the elasticity matrix
\[H = \frac{E}{(1-\nu^2)}\begin{bmatrix} 1 & \nu& 0\\ \nu & 1 & 0\\ 0 & 0 & 1-\nu \end{bmatrix}\]
Here $E$ is the Young modulus, $\nu$ is the Poisson ratio and $\varepsilon$ is the strain
\[\varepsilon = \begin{bmatrix} \varepsilon_{xx} & \varepsilon_{xy}\\ \varepsilon_{xy} & \varepsilon_{yy} \end{bmatrix}\]
The relation between $\sigma'$ and $\varepsilon$ may be nonlinear; that's why we only write $\delta \sigma'$ in terms of $\delta \varepsilon$.
The fluid mass convervation in terms of pressure and volumetric strain is given by
\[\boxed{\frac{1}{M}\frac{\partial p}{\partial t} + b\frac{\partial \varepsilon_v}{\partial t} + \mathrm{div}\mathrm{v} = f}\]
where $\varepsilon_v = \mathrm{tr} \varepsilon$, $f$ is a volumetric source term and
\[\mathbf{v} = -\frac{1}{B_f}\frac{k}{\mu}(\nabla p - \rho_f g)\]
where $k$ is the absolute permeability tensor, $\mu$ is the fluid viscosity and $B_f$ is the formation volume factor of the fluid.
The mechanical equation and fluid equation are coupled through $p$ and $\varepsilon$. In the drained split scheme, in each step $p$ is kept fixed while solving the mechanics equation and then the fluid equation is solved keeping $\varepsilon$ fixed. The drained scheme can be viewed as a Jacobian iteration of the fully coupled system.
The linear poroelasticity equations with $g=0$ can be expressed as [linear]
\[\begin{aligned} \mathrm{div}\sigma(u) - b \nabla p &= 0\\ \frac{1}{M} \frac{\partial p}{\partial t} + b\frac{\partial \varepsilon_v(u)}{\partial t} - \nabla\cdot\left(\frac{k}{B_f\mu}\nabla p\right) &= f(x,t) \end{aligned}\]
with boundary conditions
\[\begin{aligned} \sigma n = 0,\quad x\in \Gamma_{N}^u, \qquad u=0, \quad x\in \Gamma_D^u\\ -\frac{k}{B_f\mu}\frac{\partial p}{\partial n} = 0,\quad x\in \Gamma_{N}^p, \qquad p=g, \quad x\in \Gamma_D^p \end{aligned}\]
and the initial condition
\[p(x,0) = p_0,\ u(x,0) =0,\ x\in \Omega\]
Numerical Discretization
Mechanics
We discretize the domain $[0,(n-1)h]\times [0, (m-1)h]$ uniformly with step size $h$.
The finite element method is usually used to solve the mechanics equation, whose discretization reads
\[\int_{\Omega} \delta \varepsilon :\sigma'\mathrm{d}x - \int_\Omega b p \delta \varepsilon_v\mathrm{d}x = \int_{\Gamma} t\cdot\delta u\mathrm{d}s + \int_\Omega \rho_b g\cdot\delta u dx\]
where $t = \sigma n = \sigma' n - bpn$, $\Gamma$ is the part of $\partial \Omega$ with external traction, and $n$ is the unit normal vector pointing outwards. One each element $A$, define $u_A$ as the nodal values of the basis functions whose supports overlap $A$, then the strain at $(x,y)$ can be expressed as (see the figure for illustration)
\[\varepsilon_A = Bu_A, \quad \varepsilon_A =\begin{bmatrix} \varepsilon_{xx}\\ \varepsilon_{yy}\\ 2\varepsilon_{xy} \end{bmatrix},\quad u_A = \begin{bmatrix} u_1\\ u_2\\ u_3\\ u_4\\ v_1\\ v_2\\ v_3\\ v_4 \end{bmatrix}\]
where
\[B = \begin{bmatrix} \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial x} & \frac{\partial N_4}{\partial x} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & \frac{\partial N_1}{\partial y} & \frac{\partial N_2}{\partial y} & \frac{\partial N_3}{\partial y} & \frac{\partial N_4}{\partial y}\\ \frac{\partial N_1}{\partial y} & \frac{\partial N_2}{\partial y} & \frac{\partial N_3}{\partial y} & \frac{\partial N_4}{\partial y} & \frac{\partial N_1}{\partial x} & \frac{\partial N_2}{\partial x} & \frac{\partial N_3}{\partial x} & \frac{\partial N_4}{\partial x} \end{bmatrix} = \begin{bmatrix} -\frac{1-\eta}{h}&\frac{1-\eta}{h} &-\frac{\eta}{h} & \frac{\eta}{h} & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & -\frac{1-\xi}{h} & -\frac{\xi}{h} & \frac{1-\xi}{h} & \frac{\xi}{h}\\ -\frac{1-\xi}{h} & -\frac{\xi}{h} & \frac{1-\xi}{h} & \frac{\xi}{h} & -\frac{1-\eta}{h}&\frac{1-\eta}{h} &-\frac{\eta}{h} & \frac{\eta}{h} \end{bmatrix}\]
and
\[\xi = \frac{x-x_0}{h}\quad \eta = \frac{y-y_0}{h}\]
The terms in the weak form can be expressed as
\[\int_{A}\delta \varepsilon :\sigma'\mathrm{d}x = \int_A u_AB^TDB\delta u_A\mathrm{d}x\]
\[\int_A b p \delta \varepsilon_v\mathrm{d}x = \int_A bp [1,1,0]B^T\delta u_A\mathrm{d}x\]
Typically, the integration is computed using Gauss quadrature; for example, we have
\[\int_A u_AB^TDB\delta u_A\mathrm{d}x = u_A \left[\sum_{i=1}^{n_g} B(\xi_i, \eta_i)^T DB(\xi_i, \eta_i)h^2w_i\right]\delta u_A\]
where $(x_i, \eta_i)$ are Gauss quadrature points and $w_i$ is the corresponding weight.
We have the following convention for bdnode and bdedge, which denote the Dirichlet boundary conditions and the Neumann boundary conditions:
bdnode$\in \mathbf{R}^{d}$ and each entry inbdnoderepresents the node index. The corresponding row and column indices can be retrieved withfemidx.bdedge$\in \mathbf{R}^{d\times 2}$ and each row inbdedgerepresents indices of two end points of the edge. The corresponding cell row and column can be retrieved withfvmidx.
Fluid
The fluid equation is discretized using finite volume method.
\[\int_{A_i} \frac{1}{M}\frac{p_i^{n+1} - p_{i}^{n}}{\Delta t} \mathrm{d}x + \int_{A_i} b \frac{\varepsilon_v^{n+1}-\varepsilon_v^n}{\Delta t} \mathrm{d} x + \int_{A_i} \mathrm{div}\mathbf{v}\mathrm{d}x = \int_{A_i} f\mathrm{d}x\]
For the divergence term, we use the two-point flux approximation and we have (assuming $k$ is a constant scalar)
\[\int_{A_i} \mathrm{div}\mathbf{v} \mathrm{d}x = -\frac{k}{B_f\mu}\sum_{j=1}^{n_{\mathrm{faces}}} (q_j-q_i) = -\frac{k}{B_f\mu}\sum_{j=1}^{n_{\mathrm{faces}}} (p_j^{n+1} - p_i^{n+1}) + \frac{k\rho_f|g|}{B_f\mu}\sum_{j=1}^{n_\mathrm{faces}} (y_j-y_i)\]
where
\[q = p^{n+1} - \rho_f|g|y\]
Initial and Boundary Conditions
For the mechanial problem we consider
- Prescribed displacement: $u = \bar u$; or
- Prescribed traction: $\sigma\cdot n=\bar t$ (also called overburden).
For the flow problem we consider
- Prescribed pressure: $p=\bar p$; or
- Prescribed volumetric flux: $\mathbf{v}\cdot n=\bar v$ (called no flow if $\bar v=0$).
The initial displacement and strains are zero. The initial pressure is prescribed.
Verification
To verify our numerical scheme, we consider manufactured solution $u(x,y) = \begin{bmatrix} x^2+y^2\\ x^2-y^2 \end{bmatrix}t,\quad p(x,y) = x^2y^2(1-x)^2(1-y)^2e^{-t}$
Then we have
\[\begin{aligned} f(x,y,t)&= (-x^2y^2(x - 1)^2(y - 1)^2 - 2x^2y^2(x - 1)^2 - 2x^2y^2(y - 1)^2 - 8x^2y(x - 1)^2(y - 1) - 2x^2(x - 1)^2(y - 1)^2 - 8xy^2(x - 1)(y - 1)^2 - 2y^2(x - 1)^2(y - 1)^2 + 2(x - y)\exp(t))\exp(-t)\\ g(x,y,t)&= \begin{bmatrix} 3t-2xe^{-t}\\ -t + 2ye^{-t} \end{bmatrix} \end{aligned}\]
| Description | u displacement | v displacement | Pressure |
|---|---|---|---|
| Numerical Result |  |  |  |
| Error |  |  |  |
Benchmarks
Flooding
| u displacement | v displacement | Pressure | Von Mises Stress |
|---|---|---|---|
 |  |  |  |
Injection-Production in Homogenious Media
| u displacement | v displacement | Pressure |
|---|---|---|
 |  |  |
Injection-Production in Homogenious Media and with J2 plasticity
| Description | u displacement | v displacement | Pressure |
|---|---|---|---|
| $K = 0.5,\ \sigma_Y = 0.3$ |  |  |  |
| $K = 0,\ \sigma_Y = 0.3$ |  |  |  |
| $K = 0.5,\ \sigma_Y = 1.0$ |  |  |  |
Injection-Production in Heterogenious Media
| u displacement | v displacement | Pressure |
|---|---|---|
 |  |  |
- linearKolesov, Alexandr E., Petr N. Vabishchevich, and Maria V. Vasilyeva. "Splitting schemes for poroelasticity and thermoelasticity problems." Computers & Mathematics with Applications 67.12 (2014): 2185-2198.