Iteration Methods for Coupled Geomechanics and Nonlinear Multiphase Flow¶

Alejandro Queiruga¶

Lawrence Berkeley Lab¶

SIAM Conference on Mathematical & Computational Issues in the Geosciences, 2019¶

Outline¶

Focus: Improving staggered iteration between multiphase reaction transport and geomechanics

Problems we were trying to solve
Challenges we encountered
Numerical algorithms
Looking at poroelasticity as inspiration
A new preconditioned staggered iteration
Numerical experiments
Problems we can now solve

Intro¶

• Motivating problems are hydrate dissociation problems
• Wildly varying fluid compositions and mechanical properties
• Large changes causes drastic changes to the fluid mixture properties
• Properties vary greatly in space and time: can't estimate a priori
• Very soft sediments with soft gases and stiff liquids and ices

General System of Equations¶

Coupled transport, kinetics, and momentum balances:

Solve for $p,T,X$ , and $u$ such that: $\begin{eqnarray} \frac{\mathrm{d}m}{\mathrm{d}t} = & \nabla \cdot q + s & \} \text{balance of mass}\\ \frac{\mathrm{d}X}{\mathrm{d}t} = & r & \}\text{reactions (possibly equilibrium)} \\ 0 = & \nabla \cdot \sigma + f & \} \text{balance of momentum (quasistatic)} \end{eqnarray}$ Plus statemachine logic for phase change: $\begin{eqnarray} state = \text{transition}(state; p,X,T) \end{eqnarray}$

Solving the Algebraic Equations¶

At each timestep, we're solving two nonlinear algebraic problems:

$\begin{align} F(p,T,X,...)= & 0\left.\right\} \text{flow code}\\ M(u)= & 0\left.\right\} \text{mechanical code} \end{align}$

To really do Newton's method we would need: $\begin{align} \left[\begin{array}{cc} \frac{\partial F}{\partial p} & \frac{\partial F}{\partial u}\\ \frac{\partial M}{\partial p} & \frac{\partial M}{\partial u} \end{array}\right]\left\{ \begin{array}{c} \Delta p^{k}\\ \Delta u^{k} \end{array}\right\} =-\left\{ \begin{array}{c} F^{k}\\ M^{k} \end{array}\right\} \end{align}$ E.g. in hydraulic fracturing the fluid pressure must be fully coupled like this or else it won't converge.

Hard to formulate those cross-terms:

Very different discretizations on different meshes
Two codes due to seperate development legacies
Bigger matrix systems, harder to partition and solve in parallel
Won't necessarily yield better performance

General Staggered Algorithm¶

For each Newton loop $k$ in the timestep:

Solve mechanical stiffness using $K_d$ $\begin{equation} \Delta u^{k+1} = [\mathtt{K}(K_d)] \,\,\backslash\,\, \{\mathtt{R(p^k,\sigma^k)} \} \end{equation}$
Increment Stress with $K_d$ $\begin{equation} \Delta \sigma^{k+1} = \left(\left(K_d+\frac{2}{3}G\right) \nabla \cdot + 2G\nabla^s \right) : \Delta u^{k+1} + \alpha \mathbf{I} \Delta p^{k} \end{equation}$
Solve flow Jacobian with new volume $\begin{equation} \Delta p^{k+1} = [\mathtt{J}(V^{k+1})]\,\,\backslash\,\,\{\mathtt{R}(V^{k+1})\} \end{equation}$
Do state transitions $\begin{equation} state^{k+1} = \text{transition}(state^{k},p^{k+1}) \end{equation}$
Repeat until $\left|\Delta p\right|<\mathtt{tol\_flow}$ and $\left|\Delta u\right|< \mathtt{tol\_mech}$

General Algorithm¶

Staggering the two solvers is the same as this: $\begin{equation} \left[\begin{array}{cc} \frac{\partial F}{\partial p} & 0\\ 0 & \frac{\partial M}{\partial u} \end{array}\right]\left\{ \begin{array}{c} \Delta p^{k}\\ \Delta u^{k} \end{array}\right\} =-\left\{ \begin{array}{c} F^{k}\\ M^{k} \end{array}\right\} \end{equation}$ Problematic when $\left| \frac{\partial M}{\partial p} \right| > \left| \frac{\partial M}{\partial u}\right|$ : fluid stiffer than skeleton.

Easy to implement.
More iterations (no quadratic convergence!)
Potentially a lot of over-shoot during the staggering.
Exacerbates state transition logic.

This is catastrophic in large-range multiphase flow problems: $\begin{equation} F(p_{bad}) = \text{raise out of range error} \end{equation}$ Physical descriptions have limits! May never be able to solve the system.

Can we find some matrix $C$ that helps? $\begin{equation} \left[\begin{array}{cc} \frac{\partial F}{\partial p} & 0\\ 0 & \frac{\partial M}{\partial u}+C \end{array}\right]\left\{ \begin{array}{c} \Delta p^{k}\\ \Delta u^{k} \end{array}\right\} =-\left\{ \begin{array}{c} F^{k}\\ M^{k} \end{array}\right\} \end{equation}$

Poroelasticity¶

The poroelastic constitutive response for the fluid pore pressure is $\begin{equation} \mathrm{d}p=-M\alpha\nabla\cdot\mathrm{d}\mathbf{u}+M\nabla\cdot\left(\mathbf{q}\mathrm{d}t\right) \end{equation}$ where the poroelastic coefficient is $\begin{equation} M^{-1} = \phi K_f^{-1} + (\alpha - \phi) K_s^{-1} \end{equation}$ where $\mathbf{q}\mathrm{d}t$ is an differential fluid volume into the control volume.

For the total stress, the constitutive response is, $\begin{equation} \mathrm{d}\sigma=-\alpha\mathrm{d}p\mathbf{I}+\left(K_{d}-\frac{2}{3}G\right)\nabla\cdot\mathrm{d}\mathbf{u}\mathbf{I}+2G\nabla^{s}\mathrm{d}\mathbf{u}. \end{equation}$

Note: We use the law for $\sigma$ , but not $p$

Have nonlinearity in composition-dependent properties; e.g. in hydrate, $K$ and $G$ are function of saturation $S$ : $\begin{eqnarray*} K = K_{0}(1-S)+K_{1}S, \quad G = G_{0}(1-S)+G_{1}S \end{eqnarray*}$

Reworking Poroelasticity¶

We can substitute $\mathrm{d}p$ for a familiar form: $\begin{equation} \mathrm{d}\sigma=\left(K_{u}-\frac{2}{3}G\right)\nabla\cdot\mathrm{d}\mathbf{u}+2G\nabla^{s}\mathrm{d}\mathbf{u} + \alpha^2 M \nabla\cdot\left(\mathbf{q}\mathrm{d}t\right) \mathbf{I} \end{equation}$ where the undrained bulk modulus is: $\begin{equation} K_{u} = K_d + \alpha^2 M \end{equation}$

Making the analogy to our matrix system, we had (note: more variables than $p$ ) $\begin{equation} \left[\begin{array}{cc} 1 & M\alpha\mathrm{tr}\\ \alpha\mathbf{I} & K_{d} \end{array}\right]\left\{ \begin{array}{c} \mathrm{d}p\\ \mathrm{d}u \end{array}\right\} =-\left\{ \begin{array}{c} f\\ m \end{array}\right\} \end{equation}$ and are trying to do: $\begin{equation} \left[\begin{array}{cc} 1 & M\alpha\mathrm{tr}\\ 0 & K_{d}+\alpha^{2}M \end{array}\right]\left\{ \begin{array}{c} \mathrm{d}p\\ \mathrm{d}u \end{array}\right\} =-\left\{ \begin{array}{c} f\\ m - \alpha f \mathbf{I} \end{array}\right\} \end{equation}$ Triangularized the system so solving $\mathrm{d}u$ first would substitute into $\mathrm{d}p$ equation closer to the answer.

General Balance of Mass¶

Caveat: Poroelasticity only holds in the linear regime! (Small changes)

Actual balance of mass that the reservoir simulator solves is: $\begin{equation} \frac{\partial}{\partial t} \rho = \nabla \cdot \mathbf{q} \end{equation}$ where we assumed: $\begin{equation} \rho = \rho_0 ( 1+ p / K_f ) \end{equation}$ Flow simulators use very complicated nonanalytic formulations for $\rho$ .

$K_f$ and $M$ varies greatly especially at phase changes
Don't necessarily have relations for $K_f$ for every possible scenario
Varies spatially as well as temporally

Control Volume Simulations¶

Solving balance with fluxes: $\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}M_{comp}=\sum_{facets}F_{comp} \end{equation}$ Masses and fluxes are integrals over discrete domains and facets: $\begin{equation} M_{comp}=\int_{\Omega}\phi\sum_{phase}S^{phase}\rho_{comp}^{phase}X_{comp}^{phase}\mathrm{d}^{3}x \end{equation}$ and: $\begin{equation} F_{comp}=\int_{\Gamma}\sum_{phase}\mathbf{q}_{comp}^{phase}\cdot\mathbf{n}\mathrm{d}^{2}x \end{equation}$

Flow and mechanics are coupled through volume of $\left|\Omega\right| = V(t)$ $\begin{equation} \frac{\mathrm{d}V}{\mathrm{d}t} = V(t)\frac{\mathrm{d}}{\mathrm{d}t} \mathrm{tr} \varepsilon \end{equation}$ Poroelastic formulation uses the product rule to extract this term: $\begin{equation} \frac{\mathrm{d}\rho V}{\mathrm{d}t}= V \frac{\partial \rho}{\partial p}\frac{\mathrm{d}p}{\mathrm{d}t} + \rho \frac{\mathrm{d}V}{\mathrm{d}t} \end{equation}$ (Derive from transformation of $\nabla_x\cdot$ to $\nabla_X\cdot$ given $x(t)=X+u(t)$ )

Slightly Improved Staggered Algorithm¶

For each Newton loop $k$ in the timestep:

Solve mechanical stiffness using $K_u^*$ $\begin{equation} \Delta u^{k+1} = [\mathtt{K}(K_d+\alpha^2 M)] \,\,\backslash\,\, \{\mathtt{R(p^k,\sigma^k)} \} \end{equation}$
Increment Stress with $K_d$ $\begin{equation} \Delta \sigma^{k+1} = \left(\left(K_d+\frac{2}{3}G\right) \nabla \cdot + 2G\nabla^s \right) : \Delta u^{k+1} + \alpha \mathbf{I} \Delta p^{k} \end{equation}$
Solve flow Jacobian with new volume $\begin{equation} \Delta p^{k+1} = [\mathtt{J}(V^{k+1})]\,\,\backslash\,\,\{\mathtt{R}(V^{k+1})\} \end{equation}$
Do state transitions $\begin{equation} state^{k+1} = \text{transition}(state^{k},p^{k+1}) \end{equation}$
Repeat until $\left|\Delta p\right|<\mathtt{tol\_flow}$ and $\left|\Delta u\right|< \mathtt{tol\_mech}$

Use prior knowledge of how the flow solver will behave (even though the physical description is different) to precondition the iteration.

Adaptive Preconditioning¶

At each update, the flow simulator yields a $\Delta p$ based on the current composition of the fluid.

Assume we're in the linear regime within the iteration: $\begin{equation} \Delta p=-M\alpha\nabla\cdot\Delta\mathbf{u}+\nabla\cdot\left(\mathbf{q}\Delta t\right) \end{equation}$ We can estimate the effective $M$ : $\begin{equation} M^{*}=-\frac{\alpha\nabla\cdot\Delta\mathbf{u}+\nabla\cdot\left(\mathbf{q}\Delta t\right)}{\Delta p} \end{equation}$

We can back out an estimate of the bulk modulus of the fluid: $\begin{equation} K_{f}^{*}=\phi\left(\frac{1}{M^{*}}-\frac{\left(\alpha-\phi\right)}{K_{s}}\right)^{-1} \end{equation}$

Assemble the mechanical stiffness matrix using $K_{u}$ in place of $K_d$ with estimated $M^*$ : $\begin{equation} K_{u}^{*}=K_{d}+\alpha^{2}M^{*} \end{equation}$

Watch out: Not-a-number when $\Delta p=0$ , e.g. at boundary conditions.

Backing out flux divergence:¶

Estimate PDE $\nabla\cdot q$ from fluxes: $\begin{eqnarray} \nabla\cdot\mathbf{q}=&\frac{1}\\{V}\int_{\Omega}\nabla\cdot\mathbf{q}\mathrm{d}^{3}x\\=&\frac{1}\\{V}\oint_{\partial\Omega}\mathbf{q}\cdot\mathbf{n}\mathrm{d}^{2}x\\=&\frac{1}{V}\sum_{facets}\int_{\Gamma}\mathbf{q}\cdot\mathbf{n}\mathrm{d}^{2}x\\=&\frac{1}{V}\sum_{facets}F \end{eqnarray}$

Only modification needed to the flow solver; Millstone takes care of the rest.

Numerical Experiments: Setup¶

1D quasi-static consolidation with a ramp-up load filled with a mixture of $H_2O$ and $CH_4$

Verified against analytical solution for static poroelasticity
- Should get it exactly correct (if it's slow, stays linear, and we know $K_f$ )

Family of problems varied by:
- $K_ d$ skeletal stiffness
- $S_ {aqu}$ Water to gas ratio
Algorithm is permuted by:
- Choice: None, Fixed, Adaptive, Continuous
- Tolerance on convergence
Algorithm is evaluated by:
- Does it succeed
- Mean number of iterations for each set up
- Total runtime (too short)

Numerical Experiments: Results¶

Original program with no preconditioning fails when $K_d < K_f$
Fixed precondition solves that problem when $K_f$ was known
Adaptation of $K_f^*$ upon convergence improves it
Adapting during the time step iteration is robust for all compositions with reasonable tolerance
Adapting can perform worse with gaseous pore fluids
Algorithm estimates correct $K_f$ on this problem (not smooth in practice)

One more trick: under-relaxation¶

Algorithms work on 1D test problem, but still not 100\% robust in some reservoir problems.

Need to add under-relaxation: $\begin{equation} \Delta u = [\mathtt{K}(K_d+\alpha^2 M^*)] \backslash \{\mathtt{R}(\sigma)\} \end{equation}$ And update by $\begin{equation} u^{k+1} = u^k + \beta \Delta u \end{equation}$ With $0<\beta\leq 1$ . Smaller $\beta$ slows down convergence.

Also, watch out for pore ``collapse'' when gases appear ( $K_f \approx 50Pa$ vs $K_s \approx 10^{10} Pa$ )

Assumptions out of range for poroelastic porosity evolution?

Iteration Methods for Coupled Geomechanics and Nonlinear Multiphase Flow¶

Alejandro Queiruga¶

Lawrence Berkeley Lab¶

SIAM Conference on Mathematical & Computational Issues in the Geosciences, 2019¶

Outline¶

Intro¶

Challenges:¶

What we have to work with:¶

General System of Equations¶

First: How We Solve It Now¶

Problem fields are on two meshes¶

Solving the Algebraic Equations¶

General Staggered Algorithm¶

General Algorithm¶

Limited Range¶

Poroelasticity¶

Reworking Poroelasticity¶

General Balance of Mass¶

Very complicated...¶

Control Volume Simulations¶

Slightly Improved Staggered Algorithm¶

Adaptive Preconditioning¶

Backing out flux divergence:¶

Four Variations of the Algorithm:¶

Numerical Experiments: Setup¶

Four Algorithms:¶

Numerical Experiments: Results¶

Numerical Experiments: Results¶

One more trick: under-relaxation¶

Applications¶

Outro¶

References:¶

Practices:¶

Thank You¶

Acknowledgments¶