A block preconditioner for thermo-poromechanics with frictional deformation of fractures

TL;DR

This work presents a robust block preconditioner for fully implicit thermo-poromechanics with frictional, mixed-dimensional fractures. It combines a preprocessing linear transformation to address the saddle-point structure from contact mechanics with a fixed-stress decoupling for elasticity, and two options (CPR or System-AMG) to solve the coupled energy–mass subproblem. Nested Schur complement reductions (S^1, S^2, S^3) enable decoupling of mechanics, momentum, and energy transport, yielding scalable GMRES convergence across 2D and 3D fracture networks, various fracture states, and a range of Peclet numbers. The CPR variant is typically cheaper and robust in injection-dominated regimes, while System-AMG offers superior scaling in diffusion-dominated long-term simulations, with the combination validated on benchmark-like fracture geometries and provided through open tools.

Abstract

The numerical modeling of fracture contact thermo-poromechanics is crucial for advancing subsurface engineering applications, including CO2 sequestration, production of geo-energy resources, energy storage and wastewater disposal operations. Accurately modeling this problem presents substantial challenges due to the complex physics involved in strongly coupled thermo-poromechanical processes and the frictional contact mechanics of fractures. To resolve process couplings in the resulting mathematical model, it is common to apply fully implicit time stepping. This necessitates the use of an iterative linear solver to run the model. The solver's efficiency primarily depends on a robust preconditioner, which is particularly challenging to develop because it must handle the mutual couplings between linearized contact mechanics and energy, momentum, and mass balance. In this work, we introduce a preconditioner for the problem based on the nested approximations of Schur complements. To decouple the momentum balance, we utilize the fixed-stress approximation, extended to account for both the porous media and fracture subdomains. The singularity of the contact mechanics submatrix is resolved by a linear transformation. Two variations of the algorithm are proposed to address the coupled mass and energy balance submatrix: either the Constrained Pressure Residual or the System-AMG approach. The preconditioner is evaluated through numerical experiments of fluid injection into fractured porous media, which causes thermal contraction and subsequent sliding and opening of fractures. The experiments show that the preconditioner performs robustly for a wide range of simulation regimes governed by various fracture states, friction coefficients and Peclet number. The grid refinement experiments demonstrate that the preconditioner scales well in terms of GMRES iterations, in both two and three dimensions.

Figures (7)

• Figure 1: Rectangular 2D porous medium with three embedded 1D fractures and two 0D intersections. Symbols inside the circles represent the primary variables and their domains of definition. The magnifying glass on the right shows the variables defined within the fracture and on the interfaces between the fracture and the porous medium.
• Figure 2: The block structure of the matrices at different stages of the preconditioner algorithm. White cells indicate empty submatrices, light-gray cells represent non-empty submatrices from the original Jacobian, and dark-gray cells denote modified submatrices. The letter "S" marks a singular submatrix. From left to right, the stages include: the original Jacobian, the Jacobian after the linear transformation (\ref{['sec:linear_transformation']}), the first-level Schur complement (\ref{['sec:eliminating_contact_mechanics']}), the second-level Schur complement (\ref{['sec:eliminating_force_balance']}), and the third-level Schur complement (\ref{['sec:eliminating_intf_flow']}). Column 4 of the matrix refers to all the interface fluxes: $\xi_j \in \{ v_j, w_j, q_j\}$.
• Figure 3: State of the 2D model using medium coarse grid with 3e+4 degrees of freedom at $t \approx 375$ days. Left: The computational grid with bold 1D fractures and maximum fracture slip. The red dot marks the injection cell. Center: Temperature distribution. Right: Volumetric strain distribution.
• Figure 4: Number of GMRES iterations in the 2D simulation, plotted against the simulation time, for the coarsest (left) and finest (right) grids. Results with CPR and System-AMG are shown for both grids. The FGMRES solver, applied only on the coarsest grid, has subproblem solvers that are close to being exact. The observed increase in GMRES iteration counts is due to the adaptive time-stepping strategy, wherein the simulation starts with a time step of 1e-3 seconds and progressively grows to as large as 3 years.
• Figure 5: State of the 3D model using grid with 6e+5 degrees of freedom at $t \approx 30$ years. Left: Computational grid of the fracture subdomains and maximum fracture slip. The upper-left fracture is disconnected from the remainder of the fracture network with a small gap. Right: Temperature distribution within fractures and in the slice of the porous medium.