# Outline¶

Focus: Improving staggered iteration between multiphase reaction transport and geomechanics

1. Problems we were trying to solve
2. Challenges we encountered
3. Numerical algorithms
4. Looking at poroelasticity as inspiration
5. A new preconditioned staggered iteration
6. Numerical experiments
7. 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

# Challenges:¶

• Difficult to marry two physical descriptions
• One code for flow simulation, one code for geomechanics
• Slow convergence of staggered iteration
• Iteration can diverge when pore pressure dominates skeleton in stiffness
• Iteration can fail due to small valid ranges of in multiphase relations

## What we have to work with:¶

• Ideally, one would develop a fully coupled simulator from the ground up
• That is easier said than done
• Disjoint development between multiphase chemistry and computational mechanics
• The problem is hard, not just writing the software
• Current boat is having two distinct solvers
• Working in TOUGH+ and Millstone
• Multiphase flow is very tricky: make mechanics package do the heavy-lifting
• Fortran calls C calls Python calls C calls Python-generated C

# 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}

# First: How We Solve It Now¶

Pick numerical methods that are best for each set of equations.

• Reaction-transport in reservoirs needs a very fine mesh
• Integral Finite Difference Method with TOUGH+ for reaction transport
• Finite Element Method with Millstone for mechanics
• One-to-one aspect ratio problem for mechanics solved by using two meshes
• Axisymmetry in mechanics reduces unknowns and improves numerical performance

# Problem fields are on two meshes¶

Have an additional interpolation step

Queiruga, Moridis, and Reagan, "Simulation of Gas Production from Multilayered Hydrate-Bearing Media with Fully Coupled Flow, Thermal, Chemical and Geomechanical Processes Using TOUGH+Millstone. Parts 1-3," Transport in Porous Media, 2019

# 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:

1. Solve mechanical stiffness using $K_d$ $$\Delta u^{k+1} = [\mathtt{K}(K_d)] \,\,\backslash\,\, \{\mathtt{R(p^k,\sigma^k)} \}$$
2. Increment Stress with $K_d$ $$\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}$$
3. Solve flow Jacobian with new volume $$\Delta p^{k+1} = [\mathtt{J}(V^{k+1})]\,\,\backslash\,\,\{\mathtt{R}(V^{k+1})\}$$
4. Do state transitions $$state^{k+1} = \text{transition}(state^{k},p^{k+1})$$
5. 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: $$\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\}$$ 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: $$F(p_{bad}) = \text{raise out of range error}$$ Physical descriptions have limits! May never be able to solve the system.

Can we find some matrix $C$ that helps? $$\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\}$$

# Poroelasticity¶

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

For the total stress, the constitutive response is, $$\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}.$$

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: $$\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}$$ where the undrained bulk modulus is: $$K_{u} = K_d + \alpha^2 M$$

Making the analogy to our matrix system, we had (note: more variables than $p$) $$\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\}$$ and are trying to do: $$\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\}$$ 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: $$\frac{\partial}{\partial t} \rho = \nabla \cdot \mathbf{q}$$ where we assumed: $$\rho = \rho_0 ( 1+ p / K_f )$$ 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

### Very complicated...¶

def gibbs_region1(T,p):
p1_star = 1.653e7
T1_star  = 1.386e3
n1 = [ 0.14632971213167e00, -0.84548187169114e00,
-3.7563603672040e+00,  3.3855169168385e+00,
-0.95791963387872e00,  0.15772038513228e00,
-1.6616417199501e-02,  8.1214629983568e-04,
2.8319080123804e-04, -6.0706301565874e-04,
-1.8990068218419e-02, -3.2529748770505e-02,
-2.1841717175414e-02, -5.2838357969930e-05,
-4.7184321073267e-04, -3.0001780793026e-04,
4.7661393906987e-05, -4.4141845330846e-06,
-7.2694996297594e-16, -3.1679644845054e-05,
-2.8270797985312e-06, -8.5205128120103e-10,
-2.2425281908000e-06, -6.5171222895601e-07,
-1.4341729937924e-13, -4.0516996860117e-07,
-1.2734301741641e-09, -1.7424871230634e-10,
-6.8762131295531e-19,  1.4478307828521e-20,
2.6335781662795e-23, -1.1947622640071e-23,
1.8228094581404e-24, -9.3537087292458e-26]
i1 = [ 0, 0,0, 0, 0, 0, 0, 0, 1,1, 1, 1,1, 1, 2, 2, 2, 2, 2, 3, 3, 3,
4, 4,4, 5, 8,  8, 21, 23, 29, 30, 31, 32  ]
j1 = [ -2, -1, 0, 1, 2, 3, 4, 5, -9, -7, -1, 0,1, 3, -3, 0, 1,3, 17, -4,
0, 6, -5, -2, 10, -8, -11, -6, -29, -31, -38, -39, -40, -41  ]
p_i = p/p1_star
t_i = T1_star/T
return R*T*sum([ n*(7.1-p_i)**I*(t_i-1.222)**J
for n,I,J in zip(n1,i1,j1)])
density_region1,enthalpy_region1 = density_enthalpy(gibbs_region1)


# Control Volume Simulations¶

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

Flow and mechanics are coupled through volume of $\left|\Omega\right| = V(t)$ $$\frac{\mathrm{d}V}{\mathrm{d}t} = V(t)\frac{\mathrm{d}}{\mathrm{d}t} \mathrm{tr} \varepsilon$$ Poroelastic formulation uses the product rule to extract this term: $$\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}$$ (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:

1. Solve mechanical stiffness using $K_u^*$ $$\Delta u^{k+1} = [\mathtt{K}(K_d+\alpha^2 M)] \,\,\backslash\,\, \{\mathtt{R(p^k,\sigma^k)} \}$$
2. Increment Stress with $K_d$ $$\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}$$
3. Solve flow Jacobian with new volume $$\Delta p^{k+1} = [\mathtt{J}(V^{k+1})]\,\,\backslash\,\,\{\mathtt{R}(V^{k+1})\}$$
4. Do state transitions $$state^{k+1} = \text{transition}(state^{k},p^{k+1})$$
5. 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.

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: $$\Delta p=-M\alpha\nabla\cdot\Delta\mathbf{u}+\nabla\cdot\left(\mathbf{q}\Delta t\right)$$ We can estimate the effective $M$: $$M^{*}=-\frac{\alpha\nabla\cdot\Delta\mathbf{u}+\nabla\cdot\left(\mathbf{q}\Delta t\right)}{\Delta p}$$

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

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

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¶

Each data point is one simulation; missing point reflects divergence.

# 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: $$\Delta u = [\mathtt{K}(K_d+\alpha^2 M^*)] \backslash \{\mathtt{R}(\sigma)\}$$ And update by $$u^{k+1} = u^k + \beta \Delta u$$ 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?

# Applications¶

All this work to get this type of problem to work:

# Outro¶

• Hit divergence on coupled problems with drastic fluid compositional changes
• Appearance of gases makes staggered iteration raise
• Used poroelastic theory to design a way to adaptively precondition
• Permuted idea to make multiple staggering algorithms
• Tried all permutations on a simple constructed test problem

• Look for places where our computer program is arbitrary and permute them.
• Design families of tiny test problems to program their limits

These sldies are hosted at: ocf.io/afq/slides/Siam2019

### References:¶

Queiruga, Moridis, and Reagan, "Simulation of Gas Production from Multilayered Hydrate-Bearing Media with Fully Coupled Flow, Thermal, Chemical and Geomechanical Processes Using TOUGH+Millstone. Parts 1-3," Transport in Porous Media, 2019

# Practices:¶

• Convergence difficult when $K_d \lessapprox K_f$ with naive approach
• Use unstructured quad meshes for mechanics
• True 2D axisymmetry for mechanics
• Continuous $M$ estimation for preconditioner
• Under relax $\Delta u$ by $\beta \approx 0.2-0.5$
• Timestep control: double below 25 and cutback at 100
• A lot of time per iteration, but stable on large timesteps.

# Thank You¶

### Acknowledgments¶

Colleagues: Matt Reagan, George Moridis

This work is supported by various projects under the Department of Energy, and a collaboration with Chevron.

In [ ]: