A divergence preserving cut finite element method for Darcy flow

TL;DR

This work addresses accurate and stable unfitted discretizations for Darcy flow across interfaces by employing $\mathbf{RT}_k\times Q_k$ pairs within a CutFEM framework. It identifies that standard ghost penalty stabilization can destroy the divergence-free property and proposes a new mixed stabilization together with macro-element localization to preserve divergence, retain a saddle-point structure, and ensure well-posed linear systems. Theoretical results include stability, commutativity-based interpolation, and a priori error estimates showing $O(h^{k+1/2})$-type convergence for velocity and pressure, with pointwise divergence-free velocity fields for the proposed method. Numerical experiments in 2D and 3D with $\text{RT}_k\times Q_k$ and $\text{BDM}_k\times Q_k$ demonstrate optimal convergence, controlled condition numbers comparable to fitted meshes, and preserved divergence accuracy regardless of interface position, highlighting practical impact for fractured porous media simulations.

Abstract

We study cut finite element discretizations of a Darcy interface problem based on the mixed finite element pairs $\textbf{RT}_k\times Q_k$, $k\geq 0$. Here $Q_k$ is the space of discontinuous polynomial functions of degree less or equal to $k$ and $\textbf{RT}$ is the Raviart-Thomas space. We show that the standard ghost penalty stabilization, often added in the weak forms of cut finite element methods for stability and control of the condition number of the linear system matrix, destroys the divergence-free property of the considered element pairs. Therefore, we propose new stabilization terms for the pressure and show that we recover the optimal approximation of the divergence without losing control of the condition number of the linear system matrix. We prove that the method with the new stabilization term has pointwise divergence-free approximations of solenoidal velocity fields. We derive a priori error estimates for the proposed unfitted finite element discretization based on $\textbf{RT}_k\times Q_k$, $k\geq 0$. In addition, by decomposing the mesh into macro-elements and applying ghost penalty terms only on interior edges of macro-elements, stabilization is applied very restrictively and only where needed. Numerical experiments with element pairs $\textbf{RT}_0\times Q_0$, $\textbf{RT}_1\times Q_1$, and $\textbf{BDM}_1\times Q_0$ (where $\textbf{BDM}$ is the Brezzi-Douglas-Marini space) indicate that we have 1) optimal rates of convergence of the approximate velocity and pressure; 2) well-posed linear systems where the condition number of the system matrix scales as it does for fitted finite element discretizations; 3) optimal rates of convergence of the approximate divergence with pointwise divergence-free approximations of solenoidal velocity fields. All three properties hold independently of how the interface is positioned relative to the computational mesh.

Figures (13)

• Figure 1: The grey domain is $\Omega_{\mathcal{T}_{h,i}}$ and the yellow faces are faces in $\mathcal{F}_{h,i}$. Left: $i=1$. Right: $i=2$.
• Figure 1: Example 1: Solution obtained using Method 2 with the element pair $\mathbf{RT}_1\times Q_1$ and a macro-element partitioning with $\delta = 0.25$. Left: Magnitude of the approximated velocity. Right: Approximated pressure.
• Figure 2: Decomposition of the triangulation $\mathcal{T}_{h,1}$ into macro-elements for different $\delta_1$ (see \ref{['eq:largeel']}). The thick black lines illustrate the boundary of macro-elements consisting of more than one element. The yellow lines are interior faces of the macro-elements. Left: $\delta_1 = 1$. Right: $\delta_1 = 0.2$.
• Figure 2: Example 1: Results using element pairs $\mathbf{RT}_0\times Q_0$ (first row), $\mathbf{BDM}_1\times Q_0$ (second row) and $\mathbf{RT}_1\times Q_1$ (third row). Left: The $L^2$-error of the pressure versus mesh size h. Middle: The $L^2$-error of the velocity field versus mesh size h. Right: The spectral condition number versus mesh size h.
• Figure 3: Example 1: Left: The spectral condition number for different mesh sizes. Right: Condition number as a function of the smallest area of a cut region, $\min_{T \in \mathcal{G}_{h}} |T\cap \Omega_i|$. Blue diamonds represent the unstabilized method and red stars represent Method 2 with macro stabilization. \newlabelfig:condition number step0

Theorems & Definitions (20)

• Lemma 3.1
• Theorem 3.2
• Proof 1
• Lemma 4.1
• Lemma 4.2
• Lemma 4.3
• Proposition 4.4
• Proof 2
• Lemma 4.5
• Proof 3