A higher order polytopal method for contact mechanics with Tresca friction

1. Problem Statement

The paper addresses the numerical simulation of contact mechanics in fractured porous media, specifically focusing on fracture networks modeled as co-dimension 1 interfaces within a 3D domain. The physical model involves:

• Linear Isotropic Elasticity: Governed by the equilibrium equation with a linear constitutive law (Hooke's law).
• Tresca Friction: A frictional contact model where the tangential traction threshold ($g$) is independent of the normal traction. This introduces non-smoothness and inequality constraints (non-penetration and friction laws).
• Mixed-Dimensional Framework: The fractures are treated as lower-dimensional manifolds embedded in the 3D matrix, requiring a Discrete Fracture-Matrix (DFM) approach.
• Challenges: The problem involves variational inequalities, non-linear friction constraints, and the need for stability in the quasi-incompressible limit (Poisson ratio $\nu \to 0.5$), where standard low-order methods often suffer from "locking" (loss of accuracy).

2. Methodology

The authors propose and analyze a second-order Discrete de Rham (DDR) scheme. This method belongs to the family of polytopal discretization methods, offering geometric flexibility for complex fracture networks.

Key Components of the Scheme:

• Discretization Space:
  • Displacement ($u_h$): Approximated using a vector-valued space ($X_{grad,h}$) where degrees of freedom (DOFs) are attached to vertices, edges, faces, and elements. This allows for the reconstruction of piecewise quadratic displacement fields.
  • Lagrange Multiplier ($\lambda_h$): Represents the traction at the fracture interface. It is approximated by piecewise constant vectors on each fracture face.
  • Discontinuity: The displacement DOFs are allowed to be discontinuous across fracture faces (representing the jump $[[u]]$) but remain continuous across non-fracture mesh interfaces.
• Reconstruction Operators:
  • The scheme utilizes local reconstruction operators to define discrete gradients ($G^1_h$) and displacements ($\Upsilon^2_h$) of higher order (quadratic) based on the DOFs.
  • A discrete symmetric gradient $\epsilon_h$ and stress tensor $\sigma_h$ are constructed from these reconstructions.
• Formulation:
  • A mixed variational formulation is used, enforcing contact conditions via a Lagrange multiplier.
  • The discrete problem seeks $(u_h, \lambda_h)$ satisfying a saddle-point system involving a stabilization bilinear form ($S_h$) to ensure coercivity.
• Stabilization: A specific stabilization term ($S_h$) penalizes the difference between the reconstructed polynomial and the local DOFs on faces, edges, and vertices, ensuring the discrete norm is well-defined.

3. Key Contributions

The paper makes several significant theoretical and practical contributions:

1. Discrete Korn Inequality: The authors prove a discrete Korn inequality adapted to DDR spaces in fractured domains. This is crucial for establishing the coercivity of the elastic energy norm, accounting for the discontinuities across fractures.
2. Discrete Inf-Sup Condition: They establish a discrete inf-sup condition between the displacement space and the Lagrange multiplier space. Notably, this is proven using a non-standard $H^{-1/2}$-like norm tailored to the fracture network. The proof leverages a specialized Clément-like interpolator that handles the minimal regularity required near fracture tips.
3. Locking-Free Property: A major theoretical achievement is the proof that the method is robust in the quasi-incompressible limit ($\nu \to 0.5$). This is achieved by designing a bespoke Fortin interpolator that corrects cell values to ensure commutation with the divergence operator, preventing the error constants from blowing up as $\nu$ approaches 0.5.
4. Error Analysis:
   • An abstract error estimate is derived without assuming high regularity.
   • Under the assumption of $H^{3/2+r}$ regularity, the method achieves an optimal convergence rate of $O(h^{1/2+r})$ in the energy norm.
   • The analysis distinguishes between frictionless ($g=0$) and frictional ($g \neq 0$) cases, noting that while theoretical rates for friction are conservative ($r \le 1/2$), numerical results suggest higher accuracy is achievable.
5. Extension to Polyhedral Meshes: The method naturally handles general polyhedral meshes (tetrahedra, hexahedra, cut cells), overcoming limitations of traditional Finite Element Methods (FEM) that struggle with non-simplicial elements in geological applications.

4. Numerical Results

Extensive numerical experiments validate the theoretical findings:

• Convergence Rates:
  • Frictionless Case ($g=0$): The scheme demonstrates convergence rates close to order 2 for displacement and gradient in $L^2$ norms, and order 1 for the Lagrange multiplier, often exceeding theoretical lower bounds (super-convergence observed on Cartesian meshes).
  • Tresca Friction Case ($g \neq 0$): The method maintains high accuracy, with gradient convergence rates approaching order 2 and Lagrange multiplier convergence of order 1, outperforming the conservative theoretical prediction of 0.5.
• Robustness (Incompressibility):
  • Tests with varying Lamé coefficients (up to $L = 10^6$) and Poisson ratios approaching 0.5 show that the error magnitude remains stable. The method is confirmed to be locking-free, unlike standard low-order nodal schemes.
• Efficiency:
  • Comparisons with a lowest-order DDR scheme (from previous literature) show that the higher-order scheme achieves the same accuracy with significantly fewer degrees of freedom.
  • The higher-order method also yields smoother traction profiles and requires fewer Newton iterations to converge in the nonlinear solver, demonstrating superior computational efficiency.
• Complex Geometry: A 3D Discrete Fracture-Matrix model with intersecting fractures confirms the method's ability to handle complex fracture networks and contact state transitions (stick/slip) accurately.

5. Significance

This work represents a significant advancement in computational mechanics for subsurface applications (e.g., CO2 sequestration, geothermal energy, seismic risk assessment).

• Geometric Flexibility: By utilizing polyhedral meshes, the method can accurately represent complex fracture networks with corners and intersections, which are difficult to mesh with traditional FEM.
• High Accuracy with Low Cost: The higher-order nature of the scheme provides superior accuracy per degree of freedom compared to low-order methods, making large-scale 3D simulations more feasible.
• Physical Fidelity: The robustness in the incompressible limit is critical for modeling rock mechanics where materials are often nearly incompressible.
• Theoretical Rigor: The paper provides a rigorous mathematical foundation (Korn inequality, inf-sup condition, error estimates) for applying DDR methods to contact mechanics, a field previously dominated by lower-order or Nitsche-based approaches.

In summary, the paper successfully extends the Discrete de Rham framework to high-order contact mechanics with friction, offering a mathematically sound, robust, and computationally efficient tool for simulating fractured porous media.