Augmented Lagrangian Solvers for Poroelasticity with Fracture Contact Mechanics

Abstract

In the subsurface, fractures and the surrounding porous rock can deform in interaction with fluid flow. Advanced mathematical models governing these coupled processes typically combine fluid flow, poroelasticity, and fracture contact mechanics. The resulting system of equations is complex and highly nonlinear. As a result, convergence issues with nonlinear solvers are common, causing a bottleneck for the numerical solution of such models. One source of difficulty for the nonlinear solvers comes from the fracture contact mechanics, due to its inherently nonsmooth character. In addition, depending on the chosen constitutive model, the degree of nonlinearity is increased through coupling of flow and contact mechanics. In this paper, we investigate solvers based on the augmented Lagrangian formulation of the frictional contact problem. This includes two classical solvers, namely the generalized Newton method (using complementarity functions) and the return map method (equivalent to an Uzawa method). In addition, we propose a new solver that combines features of both approaches. Numerical experiments in two and three dimensions, designed to simulate hydraulic stimulation of geothermal reservoirs, are conducted to assess the performance of the solvers on problems of poromechanics with fracture contact mechanics. The return map method has more difficulty handling the nonlinear coupling between flow and contact mechanics than the other solvers, in many cases not converging or using an excessive number of iterations. Our new combined solver performs the most robustly across the experiments, its performance being less sensitive to the value of the augmentation parameter than the other solvers.

Figures (10)

• Figure 1: Illustration of the mixed-dimensional geometry. The higher-dimensional subdomain $\Omega_h$ is coupled to the lower-dimensional subdomain $\Omega_l$ through the interfaces $\Gamma_j$ and $\Gamma_k$. Note that $\partial_j \Omega_h, \partial_k \Omega_h, \Gamma_j, \Gamma_k$ and $\Omega_l$ all coincide geometrically. Figure adapted from Stefansson_2021
• Figure 2: Left: Fracture network used in the two-dimensional simulations, with the injection well marked by a red dot. Right: Injection schedule for the two-dimensional simulations. The wellhead pressure is momentarily increased by 1 MPa every hour, and otherwise kept constant. The same injection schedule will also be used for the three-dimensional simulations in section \ref{['subsec:three_dim_ex']}
• Figure 3: Visualization of the simulations of two-dimensional injection. Results are shown at the end of each simulation phase using model C, and at the end of the final simulation phase using models A and B. Left: Pressure profile and contact states. Here, slip is measured in a cumulative manner; a cell is regarded as being in a "slip" state if the fracture is closed and the norm of the tangential displacement jump is positive, using the state before the injection (at $t=0$) as the reference. There is no opening of fractures. Right: Changes in volumetric aperture and displacement magnitude, relative to the respective values at $t=0$
• Figure 4: Bar plots showing the cumulative number of nonlinear iterations, summed over all time steps, for the example of section \ref{['sec:two_dim_sim']} with the different degrees of nonlinear couplings between flow and mechanics outlined in table \ref{['tab:models_experiments']}. The three different colors per bar correspond to the three different phases of the simulation, after each pressure increase in the injection well. Hatched black lines indicate steps where the nonlinear solver failed to converge, and red crosses indicate a simulation that was terminated. The simulation is terminated if one of three cases occur: The number of nonlinear iterations exceed 800, the nonlinear solver failed to converge for the smallest allowed time step, or after 6 unsuccessful recomputation attempts. Finally, the total number of iterations of the linear solver is marked by a black dot
• Figure 5: Residual norms at every nonlinear iteration for the simulation using model A and $c=10^3$ (left), and the simulation using model C and $c=10^{-3}$ (right). Only the first attempted time step is shown for all solvers