An efficient preconditioner for mixed-dimensional contact poromechanics based on the fixed stress splitting scheme

Technical Summary: An Efficient Preconditioner for Mixed-Dimensional Contact Poromechanics Based on the Fixed Stress Splitting Scheme

Problem Statement
The numerical simulation of fracture contact poromechanics is critical for applications such as CO2 sequestration, geothermal energy production, and underground gas storage. These problems involve strongly coupled physics: fluid flow, matrix deformation (Biot poromechanics), and frictional contact mechanics within fractures. Accurate modeling is challenging due to the complex geometry of fractures (often represented as lower-dimensional objects embedded in a higher-dimensional domain) and the non-smooth nature of contact conditions (stick, slide, and open states governed by Coulomb friction and non-penetration constraints).

When solved using fully implicit time-stepping and Newton's method, these problems generate large, fully coupled linear systems. The efficiency and robustness of the solver depend heavily on the preconditioner. Existing methods face specific hurdles:

1. Saddle-Point Structure: The Augmented Lagrange formulation of contact mechanics introduces a saddle-point structure with a singular submatrix, which complicates the application of standard preconditioning techniques.
2. Coupling Complexity: The fixed stress splitting scheme, a well-established method for decoupling poromechanics, is difficult to extend to fractured media with contact mechanics. Standard saddle-point solvers often assume specific matrix structures (e.g., a Laplacian-like operator) that do not hold when the operator encapsulates the coupled poromechanical problem.
3. Scalability: Preconditioners must remain efficient as problem sizes grow to include large fracture networks typical of real-world reservoir simulations.

Methodology
The authors propose a robust and scalable preconditioner that combines nested Schur complement approximations with a linear transformation. The approach extends the fixed stress scheme to both matrix and fracture subdomains while addressing the singular nature of contact mechanics.

The methodology proceeds through the following stages:

1. Linear Transformation (Preprocessing): To handle the singularity in the contact mechanics submatrix (specifically for sticking fracture states), a right linear transformation matrix ($Q_r$) is applied to the Jacobian. This transformation renders the contact mechanics submatrix invertible by coupling it with the interface force balance equations, inspired by techniques used in other saddle-point problems.
2. Nested Schur Complement Decomposition: The preconditioner approximates the upper triangular factor of the block LDU decomposition of the transformed Jacobian through three levels of Schur complement elimination:
   • Level 1 (Contact Mechanics): Eliminates contact mechanics degrees of freedom. The transformed contact submatrix is inverted exactly (or via block-diagonal approximation) to decouple contact conditions from interface forces.
   • Level 2 (Momentum Balance): Applies a fixed stress-based approximation to decouple the momentum balance equations. This involves constructing a stabilization matrix ($D_{55}$) with specific coefficients for the ambient matrix ($l_{mat}$) and fracture subdomains ($l_{frac}$). The coupled momentum and interface force balance equations are solved monolithically using Algebraic Multigrid (AMG).
   • Level 3 (Interface Flow): Eliminates interface fluid flow by approximating the interface flow submatrix with its diagonal (ILU(0) is used for the interface flow system). The resulting Schur complement, representing the fluid mass balance, is solved monolithically with AMG.
3. Theoretical Analysis: The authors provide a theoretical justification for the fixed stress splitting scheme in this context. By analyzing a simplified linearized model, they prove that the sequential iterative scheme is a contraction under specific assumptions regarding material properties, friction coefficients, and the stability of the contact state. This analysis validates the choice of the stabilization coefficients used in the preconditioner.

Key Contributions

• Novel Preconditioner Design: The paper introduces a preconditioner that successfully integrates the fixed stress splitting scheme with contact mechanics in mixed-dimensional poromechanics. It overcomes the singularity of the contact submatrix via a linear transformation and handles the complex coupling through nested Schur complements.
• Theoretical Convergence Proof: The authors derive convergence guarantees for the fixed stress splitting scheme applied to fractured contact poromechanics. They demonstrate that the scheme is a contraction under physically reasonable assumptions, linking the convergence rate to friction coefficients and contact elasticity parameters.
• Scalability and Robustness: The method is designed to be scalable, utilizing AMG for the dominant subproblems (elasticity and fluid mass balance) and ILU for interface flows, making it suitable for large-scale problems with complex fracture networks.

Results
The preconditioner was validated through a series of numerical experiments using the PorePy framework:

• Parameter Sensitivity: Tests on a single-fracture model showed that the iteration count depends primarily on the friction coefficient ($F$), with higher friction requiring more iterations, consistent with the theoretical analysis. The fracture normal stiffness ($K_n$) had little effect on convergence.
• Grid Refinement (2D): In a 2D model with a single fracture, the number of GMRES iterations remained stable when using direct subsolvers. When replacing direct solvers with AMG, the iteration count increased slightly (by a factor of ~2) as the degrees of freedom increased by three orders of magnitude, demonstrating good scalability.
• Geometric Complexity (2D): Experiments with multiple intersecting fractures showed that the preconditioner maintained efficiency, with only a small increase in GMRES iterations for AMG subsolvers as the grid was refined.
• 3D Scalability: A 3D simulation with multiple fractures and nearly $10^6$ degrees of freedom demonstrated that the GMRES iteration count remained almost constant with grid refinement, confirming the method's scalability for large-scale problems. No linear solver breakdowns were observed in the 3D case.

Significance
The paper claims that the proposed preconditioner offers a robust and scalable solution for the challenging problem of frictional contact poromechanics in fractured media. By extending the fixed stress scheme to handle the saddle-point structure of contact mechanics via linear transformation and nested Schur complements, the method enables efficient simulation of complex subsurface processes. The authors emphasize that the approach is applicable to any discretization of contact (poro)mechanics that utilizes Lagrange multiplier-based approaches, providing a pathway for accurate and computationally feasible modeling of induced seismicity, fluid transport, and deformation in fractured reservoirs.