Auto-Stabilized Weak Galerkin Finite Element Methods for Biot's consolidation model on Non-Convex Polytopal Meshes

1. Problem Statement

The paper addresses the numerical simulation of Biot's consolidation model, a coupled system of partial differential equations (PDEs) describing fluid flow within deformable porous media. The model involves two primary variables:

• Solid displacement ($u$): A vector field representing the deformation of the porous skeleton.
• Fluid pressure ($p$): A scalar field representing the pore pressure.

Key Challenges:

• Numerical Instability: Traditional finite element methods often suffer from non-physical oscillations in the pressure field, particularly in regimes with low permeability or small time steps. This is often attributed to the violation of the inf-sup (Ladyzhenskaya-Babuška-Brezzi) condition or lack of monotonicity.
• Mesh Flexibility: Standard methods often struggle with complex geometries. There is a need for robust methods that work on non-convex polytopal meshes (general polygonal/polyhedral elements), which are common in real-world applications but difficult to handle with classical finite elements.
• Stabilization Complexity: Existing Weak Galerkin (WG) methods typically require manually tuned stabilization parameters or specific stabilizer terms to ensure stability, increasing implementation complexity.

2. Methodology

The authors propose a Fully Implicit Auto-Stabilized Weak Galerkin (WG) Finite Element Method.

A. Mathematical Formulation

• Spatial Discretization: The method uses the Weak Galerkin framework. Instead of standard continuous derivatives, it employs weak differential operators (weak gradient, weak divergence, and weak strain tensor).
  • The displacement $u$ is approximated by a pair $\{u_0, u_b\}$, where $u_0$ is the interior value and $u_b$ is the boundary value.
  • The pressure $p$ is approximated similarly by $\{p_0, p_b\}$.
  • Crucially, the boundary values ($u_b, p_b$) are treated as independent variables from the traces of the interior values, allowing for discontinuities across element interfaces.
• Temporal Discretization: A Backward Euler scheme is used for time stepping, ensuring a fully implicit and unconditionally stable formulation.
• Auto-Stabilization Strategy:
  • Unlike traditional WG methods that add explicit stabilizer terms (e.g., penalty terms on jumps), this method achieves stability without traditional stabilizers.
  • Stability is achieved by selecting a higher-degree polynomial space for the computation of the discrete weak differential operators.
  • Specifically, for an element $T$ with $N$ faces/edges:
    • If $T$ is convex, the weak operator degree is $r = k - 1 + N$.
    • If $T$ is non-convex, the degree is $r = k - 1 + 2N$.
  • This higher-degree space allows the method to capture the necessary physics and satisfy stability conditions naturally.

B. Theoretical Analysis

• Well-posedness: The authors prove that the resulting linear system at each time step is well-posed using the inf-sup condition and standard Babuška theory.
• Error Analysis: A rigorous error analysis is conducted. By utilizing bubble functions as a primary analytical tool, the authors derive optimal-order convergence rates for both displacement and pressure under standard regularity assumptions.
• Locking-Free: The formulation is proven to be locking-free (avoiding volumetric locking in nearly incompressible materials) and stable for both convex and non-convex geometries.

3. Key Contributions

1. Elimination of Stabilizers: The most significant innovation is the removal of traditional stabilizer terms. The method relies solely on the enrichment of the polynomial degree for weak operators to ensure stability, simplifying the implementation and removing the need for parameter tuning.
2. Non-Convex Mesh Capability: The method is explicitly designed and proven to work on non-convex polytopal meshes. This extends the applicability of WG methods to highly complex geometries where convexity cannot be guaranteed.
3. Oscillation-Free Pressure: The scheme produces stable, oscillation-free pressure approximations without requiring mass lumping or specialized post-processing, even in low-permeability scenarios.
4. Unified Framework: It provides a unified approach for the Biot system that handles both the mechanical (displacement) and hydraulic (pressure) fields with high flexibility.

4. Numerical Results

The authors validate the theory through extensive numerical experiments on the unit square domain using various non-convex polygonal grids (G1–G4).

• Convergence Rates:
  • Tests were performed using $P_k$ elements for displacement and pressure ($k=1, 2, 3$).
  • The results consistently demonstrated optimal-order convergence in both $L^2$ and energy norms for displacement and pressure.
  • Convergence rates were observed to be independent of the Poisson ratio ($\nu$), even in the nearly incompressible limit ($\nu \to 0.5$).
• Robustness:
  • Example 1 (Transient): Validated the time-dependent convergence and stability.
  • Example 2 (Steady-state with Discontinuous Permeability): Simulated a scenario with a sharp internal layer in permeability ($K_0 = 10^{-6}$). The method successfully captured the sharp pressure interface without typical numerical oscillations or interface smearing, demonstrating superior performance over standard methods in handling discontinuities.
• Efficiency: The method maintained the global sparsity and size of the stiffness matrix comparable to standard methods, despite the higher-degree internal computations for weak operators.

5. Significance

This paper represents a significant advancement in computational poroelasticity and finite element analysis:

• Practical Applicability: By supporting non-convex polytopal meshes, the method allows for more accurate modeling of complex geological formations and fractured media without the need for excessive mesh refinement or element splitting.
• Simplified Implementation: The "auto-stabilized" nature removes the burden of selecting stabilization parameters, making the method more robust and easier to implement for practitioners.
• Theoretical Rigor: The proof of optimal convergence and well-posedness on non-convex domains fills a gap in the theoretical understanding of WG methods for coupled multiphysics problems.
• Stability: The ability to handle low permeability and nearly incompressible materials without spurious oscillations makes this method highly valuable for geotechnical engineering, petroleum reservoir simulation, and biomechanics.

In summary, Wang and Zhang present a robust, parameter-free, and geometrically flexible numerical scheme that overcomes long-standing stability issues in Biot's consolidation modeling, offering a powerful tool for simulating complex poroelastic phenomena.