A cut finite element method for a model of pressure in fractured media

1 Introduction

The numerical modelling of flow in fractured porous media is important both in environmental science and in industrial applications. It is therefore not surprising that it is currently receiving increasing attention from the scientific computing community. Here we are interested in models where the fractures are modelled as embedded surfaces of dimension $d - 1$ in a $d$ dimensional bulk domain. Models on this type of geometries of mixed dimension are typically obtained by averaging the flow equations across the width of the fracture and introducing suitable coupling conditions for the modelling of the interaction with the bulk flow. Such reduced modelled have been derived for instance in [20, 24, 1]. The coupling conditions in these models typically take the form of a Robin type condition. The physical properties of the coupling enters as parameters in this interface condition. The size of these parameters can vary with several orders of magnitude depending on the physical properties of the crack and of the material in the porous matrix. This makes it challening to derive methods that both are flexible with respect to mesh geometries and robust with respect to coupling conditions. A wide variety of different strategies for the discretisation of fractured porous media flow has been proposed in the literature. One approach is to use a method that allows for nonconforming coupling between the bulk mesh and the fracture mesh [3], or even arbitrary polyhedral elements in the bulk mesh in order to be able to mesh the fractures easily. This latter approach has been developed using discontinuous Galerkin methods [2], virtual element methods [15] and high order hybridised methods [13].

Herein we will consider an unfitted approach, drawing on previous work [4, 12, 10] where flow in fractured porous media was modelled in the situation where the pressure is a globally continuous function. When using unfitted finite element methods, the bulk mesh can be created completely independently of the fractures. Instead the finite element space is modified locally to allow for discontinuities across fractures and interface conditions are typically imposed weakly, or using methods similar to Nitsche’s method. For other recent work using unfitted methods we refer to [22], where a stabilized Lagrange multiplier method is considered for the interface coupling and [14] where a mixed method is considered for the Darcy’s equations both in the bulk and on the surface.

The upshot here, compared to [10] is that the pressure in the crack has its own pressure, allowing for the accurate approximation of problems where the pressure is discontinuous between the bulk and the fracture, and that the interface conditions are imposed in a way allowing for the full range of parameter values in the Robin condition, without loss of stability. We use the variant of the interface modelling considered in [24], that was also recently applied for the numerical modelling in [2]. In these models we may obtain a wide range of parameter values in the interface condition and we therefore develop a method which handle the full range of values and produces approximations with optimal order convergence. The approach is inspired by the work of Stenberg [21] and may be viewed as a version of the Nitsche method that can handle Robin type conditions and which converges to the standard Nitsche method when the Robin parameter tends to infinity. Previous applications of this approach in the context of fitted finite elements include [18] and [28].

The finite element spaces are constructed starting from a standard mesh equipped with a finite element space. For each geometric domain (subdomains and interface) we mark all the elements intersected by the domain and then we restrict the finite element space to that set to form a finite element space for each domain. This procedure leads to cut finite elements and we use stabilization to ensure that the resulting form associated with the method is coercive and that the stiffness matrix is well conditioned. The stabilization is of face or ghost penalty type [5, 6, 23] , and is added both to the bulk and interface spaces. Previous related work work on cut finite element methods include the interface problem [17]; overlapping meshes [19]; coupled bulk-surface problems [11, 7, 10] and [16]; mixed dimensional problems [8]

, and surface partial differential equations

[6, 25]. For a general introduction to cut finite element methods we refer to [4].

The outline of the paper is as follows: In Section 2 we introduce the model problem and discuss the relation between our formulation of the interface conditions and previous work; in Section 3 we formulate the cut finite element method; in Section 4 we prove the basic properties of the formulation and in particular an optimal order a priori error estimate which is uniform in the full range of interface parameters; and in Section 5 we present numerical results.

2 The Model Problem

2.1 Governing Equations

Let $Ω$ be a convex polygonal domain in $R^{d}$ , $d = 2$ or $3$ , with boundary $\partial Ω$ and exterior unit normal $n$ . Let $Γ$ be a smooth embedded interface in $Ω$ , which partitions $Ω$ into two subdomains $Ω_{1}$ and $Ω_{2}$ with exterior unit normals $n_{1}$ and $n_{2}$ . We assume that $Γ$ is a closed surface without boundary residing in the interior of $Ω$ , more precisely we assume that there is $δ_{0} > 0$ such that the distance between $Γ$ and $\partial Ω$ is larger than $δ$ . We consider for simplicity the case with homogeneous Dirichlet conditions on $\partial Ω$ .

The problem takes the form: find $u_{i} : Ω_{i} \to R$ and $u_{Γ} : Γ \to R$ such that

$- \nabla \cdot A_{i} \nabla u_{i}$	$= f_{i}$	in $Ω_{i}$	(2.1)
$- \nabla_{Γ} \cdot A_{Γ} \nabla_{Γ} u_{Γ}$	$= f_{Γ} + ⟦ n \cdot A \nabla u ⟧$	on $Γ$	(2.2)
$n \cdot A \nabla u + B (u - u_{Γ})$	$= 0$	on $Γ$	(2.3)
$u$	$= 0$	on $\partial Ω$	(2.4)

Here the jump (or sum) of the normal fluxes is defined by

⟦ n \cdot A \nabla v ⟧ = 2 \sum i = 1 n_{i} \cdot A_{i} \nabla v_{i},

(2.5)

In the interface condition (2.3), $B$ is a $2 \times 2$

symmetric matrix with eigenvalues

λ_{i}

such that

λ_{i} \in [0, \infty)

and we used the notation

n \cdot A \nabla v = [\begin{matrix} n_{1} \cdot A_{1} \nabla v_{1} n_{2} \cdot A_{2} \nabla v_{2} \end{matrix}], v - v_{Γ} = [\begin{matrix} v_{1} - v_{Γ} v_{2} - v_{Γ} \end{matrix}]

(2.6)

and thus in component form (2.3) reads

[\begin{matrix} n_{1} \cdot A_{1} \nabla u_{1} n_{2} \cdot A_{2} \nabla u_{2} \end{matrix}] + B [\begin{matrix} u_{1} - u_{Γ} u_{2} - u_{Γ} \end{matrix}] = [\begin{matrix} 00 \end{matrix}]

(2.7)

The coefficients $A_{1}$ , $A_{2}$ , are smooth uniformly positive definite symmetric $d \times d$ matrices, $A_{Γ}$ is smooth tangential to $Γ$ and uniformly positive definite on the tangent space of $Γ$ , so that

2 \sum i = 1 ∥ \nabla v_{i} ∥_{Ω_{i}}^{2} + ∥ \nabla_{Γ} v_{Γ} ∥_{Γ}^{2} ≲ 2 \sum i = 1 (A_{i} \nabla v_{i}, \nabla v_{i})_{Ω_{i}} + (A_{Γ} \nabla_{Γ} v_{Γ}, \nabla_{Γ} v_{Γ})

(2.8)

where $≲$ denotes less or equal up to a constant. Finally we assume $f_{i} \in L_{2} (Ω_{i})$ and $f_{Γ} \in L_{2} (Γ)$ .

Remark 2.1

Several generalizations are possible on the external boundary. For instance, we may let the interface intersect the boundary of $Ω$ . In this case we let $ν$ denote the unit exterior conormal to $Γ \cap \partial Ω$ , i.e. $ν$ is tangent to $Γ$ and normal to $\partial Ω \cap Γ$ , and we assume that $ν \cdot n \geq c > 0$ for some constant $c$ so that the interface is transversal to $\partial Ω$ . We may then enforce the Dirichlet condition $u_{Γ} = g_{Γ}$ on $\partial Ω \cap Γ$ (see [9]) or some other standard boundary condition.

Remark 2.2

In practical modeling we may want to take the thickness of the interface inte account. Assuming that the permeability matrix in an interface of thickness $t$ takes the form

A |_{U_{t / 2} (Γ)} = A_{Γ}^{e} + a_{Γ}^{e} n_{Γ} \otimes n_{Γ}

(2.9)

where $n_{Γ}$

is a unit normal vector field to

Γ

U_{t / 2} (Γ)

is the set of points with distance less than

t / 2

Γ

v^{e}

denotes the extension of a function

v

Γ

that is constant in the normal direction,

A_{Γ}

is the tangential tangential permeability tensor, and finally

a_{Γ, n}

is the permeability across the interface. Also assuming that

f = f_{Γ}^{e}

and

u = u^{e}

U_{t / 2} (Γ)

, the equation on the interface (2.2) may be modelled as follows

- \nabla_{Γ} \cdot t A_{Γ} \nabla_{Γ} u_{Γ}

= t f_{Γ} + ⟦ n \cdot A \nabla u ⟧

Γ

(2.10)

Note that the last term on the right hand side does not scale with $t$ since it accounts for flow into the crack from the bulk domains.

Remark 2.3

We comment on how our interface condition (2.3) relates to the condition in [24] and later reformulated, see [2], in terms of averages and jumps of the bulk fields across the interface. The interface conditions in [24], equations (3.18) and (3.19), takes the form

	$ξ n_{1} \cdot A_{1} \nabla v_{1} - (1 - ξ) n_{2} \cdot A_{2} \nabla v_{2}$	$= α (v_{Γ} - v_{1})$		(2.11)
	$ξ n_{2} \cdot A_{2} \nabla v_{2} - (1 - ξ) n_{1} \cdot A_{1} \nabla v_{1}$	$= α (v_{Γ} - v_{2})$		(2.12)

where $ξ$ and $α$ are parameters. The parameter $α$ is related to physical properties of the interface as follows

α = \frac{2 a_{Γ, n}}{d}

(2.13)

where $a_{Γ, n}$ is the permeability coefficient across the interface $Γ$ and $d$ is the thickness of the interface, see (3.8) in [24] In matrix form we obtain

[\begin{matrix} ξ & ξ - 1 ξ - 1 & ξ \end{matrix}] [\begin{matrix} n_{1} \cdot A_{1} \nabla v_{1} n_{2} \cdot A_{2} \nabla v_{2} \end{matrix}] + [\begin{matrix} α & 0 0 & α \end{matrix}] [\begin{matrix} v_{1} - v_{Γ} v_{2} - v_{Γ} \end{matrix}] = 0

(2.14)

which leads to

B = \frac{1}{2 ξ - 1} [\begin{matrix} ξ & 1 - ξ 1 - ξ & ξ \end{matrix}] [\begin{matrix} α & 0 0 & α \end{matrix}] = \frac{α}{2 ξ - 1} [\begin{matrix} ξ & 1 - ξ 1 - ξ & ξ \end{matrix}]

(2.15)

We note that we have the eigen pairs

B e_{1} = \frac{α}{2 ξ - 1}      λ_{1} e_{1}, B e_{2} = α    λ_{2} e_{2}

(2.16)

with the corresponding eigen vectors defined by

and thus $B$ is positive definite for $ξ > 1 / 2$ , singular for $ξ = 1 / 2$ , and indefinite for $ξ < 1 / 2$ . It is therefore natural to consider the case when $α > 0$ and $ξ > 1 / 2$ . We remark that when $α$ tends to infinity both eigenvalues tend to infinity and when $ξ$ tends to $1 / 2$ from above one eigenvalue tends to infinity. It is therefore important to construct a method which is robust in the full range $λ_{i} \in (0, \infty)$ of possible values for the two eigenvalues

To see the relation to the formulation of the interface conditions in [2] we first note that we have the expansions

[\begin{matrix} n_{1} \cdot A_{1} \nabla v_{1} n_{2} \cdot A_{2} \nabla v_{2} \end{matrix}] = 2^{- 1 / 2} ⟦ n \cdot A \nabla ⟧ e_{1} + 2^{1 / 2} ⟨ ⟨ n \cdot A \nabla v ⟩ ⟩ e_{2}

(2.17)

[\begin{matrix} v_{1} - v_{Γ} v_{2} - v_{Γ} \end{matrix}] = 2^{1 / 2} (⟨ ⟨ v ⟩ ⟩ - v_{Γ}) e_{1} + 2^{- 1 / 2} ⟦ v ⟧ e_{2}

(2.18)

where the jumps and averages of the bulk fields across the the interface are defined by

⟦ n \cdot A \nabla v ⟧ = 2 \sum i = 1 n_{i} \cdot A_{i} \nabla v_{i}, ⟦ v ⟧ = v_{1} - v_{2}

(2.19)

⟨ ⟨ n \cdot A \nabla v ⟩ ⟩ = \frac{1}{} 2 (n_{1} \cdot A_{1} \nabla v_{1} - n_{2} \cdot A_{2} \nabla v_{2}), ⟨ ⟨ v ⟩ ⟩ = \frac{1}{2} (v_{1} + v_{2})

(2.20)

Using the expansions (2.17) and (2.18) together with (2.16

) and matching the coefficients associated with each eigenvector we obtain the interface conditions

	$⟦ n \cdot A \nabla v ⟧ + \frac{2 α}{2 ξ - 1} (⟨ ⟨ v ⟩ ⟩ - v_{Γ}) = 0$		(2.21)
	$⟨ ⟨ n \cdot A \nabla v ⟩ ⟩ + \frac{α}{2} ⟦ v ⟧ = 0$		(2.22)

which are precisely the conditions used in [2].

2.2 Weak Form

Define the function spaces

$V$	$= V_{1} \oplus V_{2} \oplus V_{Γ}$	(2.23)
$V_{i}$	$= {v_{i} \in H^{1} (Ω_{i}) : v = 0 on \partial Ω \cap \partial Ω_{i}} i = 1, 2$	(2.24)
$V_{Γ}$	$= {v_{Γ} \in H^{1} (Γ) : v = 0 on \partial Ω \cap Γ}$	(2.25)

and let $v \in V$ denote the vector $v = (v_{1}, v_{2}, v_{Γ})$ . We will also use the notation $~ V$ for functions $v \in V$ such that $v_{i} \in H^{\frac{3}{2} + ϵ} (Ω_{i})$ , $i = 1, 2$ , and $v_{Γ} \in H^{\frac{3}{2} + ϵ} (Γ)$ , with $ϵ > 0$ . Using partial integration on $Ω_{i}$ we obtain

	$2 \sum i = 1 (f_{i}, v_{i})_{Ω_{i}}$	$= 2 \sum i = 1 (- \nabla \cdot A_{i} \nabla u_{i}, v_{i})_{Ω_{i}}$
		$= 2 \sum i = 1 (A_{i} \nabla u_{i}, \nabla v_{i})_{Ω_{i}} - (n_{i} \cdot A_{i} \nabla u_{i}, v_{i})_{\partial Ω_{i}}$
		$= 2 \sum i = 1 (A_{i} \nabla u_{i}, \nabla v_{i})_{Ω_{i}} - (n_{i} \cdot A_{i} \nabla u_{i}, v_{i} - v_{Γ})_{\partial Ω_{i}} - (n_{i} \cdot A_{i} \nabla u_{i}, v_{Γ})_{\partial Ω_{i}}$
		$= 2 \sum i = 1 (A_{i} \nabla u_{i}, \nabla v_{i})_{Ω_{i}} - (n \cdot A \nabla u, v - v_{Γ})_{Γ} - (⟦ n \cdot A \nabla u ⟧, v_{Γ})_{Γ}$
		$= 2 \sum i = 1 (A_{i} \nabla u_{i}, \nabla v_{i})_{Ω_{i}} + (B (u - u_{Γ}), v - v_{Γ})_{Γ}$
		$+ (A_{Γ} \nabla_{Γ} u_{Γ}, \nabla_{Γ} v_{Γ})_{Γ} - (f_{Γ}, v_{Γ})_{Γ}$

Thus we arrive at the weak problem: find $u = (u_{1}, u_{2}, u_{Γ}) \in V$ such that

(2.26)

where the forms are defined by

	$A (u, v)$	$= 2 \sum i = 1 (A_{i} \nabla u_{i}, \nabla v_{i})_{Ω_{i}} + (A_{Γ} \nabla_{Γ} u_{Γ}, \nabla_{Γ} v_{Γ})_{Γ} + (B (u - u_{Γ}), v - v_{Γ})_{Γ}$		(2.27)
	$L (v)$	$= 2 \sum i = 1 (f_{i}, v_{i})_{Ω_{i}} + (f_{Γ}, v_{Γ})_{Γ}$		(2.28)

2.3 Existence and Uniqueness

Introducing the energy norm

⫴ v ⫴^{2} = 2 \sum i = 1 ∥ v ∥_{H^{1} (Ω_{i})}^{2} + ∥ v_{Γ} ∥_{H^{1} (Γ)}^{2} + ∥ v - v_{Γ} ∥_{Γ}^{2}

(2.29)

on $V$ , we directly find using a Poincaré inequality and the Cauchy-Schwarz inequality that the form $A$ is coercive and continuous

⫴ v ⫴^{2} ≲ A (v, v), A (v, w) ≲ ⫴ v ⫴ ⫴ w ⫴

(2.30)

Furthermore, $L$ is a continuous functional on $V$ and it follows from the Lax-Milgram Lemma that there is a unique solution in $V$ to (2.26).

In the case considered here where $Γ$ is a smooth, closed interface the model problem (2.26) satisfies the elliptic regularity estimate

∥ u_{1} ∥_{H^{2} (Ω_{1})} + ∥ u_{2} ∥_{H^{2} (Ω_{2})} + ∥ u_{Γ} ∥_{H^{2} (Γ)} ≲ ∥ f_{1} ∥_{Ω_{1}} + ∥ f_{2} ∥_{Ω_{2}} + ∥ f_{Γ} ∥_{Γ}

(2.31)

This follows in a straightforward manner from standard regularity theory. First note that since $u_{i} \in H^{1} (Ω_{i})$ , $i = 1, 2,$ and $u_{Γ} \in H^{1} (Γ)$ we have $B (u - u_{Γ}) |_{Γ} \in [H^{\frac{1}{2}} (Γ)]^{2}$ and using (2.3), $n \cdot A \nabla u \in [H^{\frac{1}{2}} (Γ)]^{2}$ . This means that the right hand side of (2.2) is in $L^{2}$ and hence $u_{Γ} \in H^{2} (Γ)$ by elliptic regularity. Considering once again (2.3) we see that in each subdomain the solution coincides with a single domain solution with a Robin condition with data in $H^{\frac{1}{2}} (Γ)$ on $Γ$ . By the elliptic regularity of the Robin problem we can then conclude that (2.31) holds.

3 A Robust Finite Element Method

3.1 The Mesh and Finite Element Spaces

To formulate the finite element method we introduce the following notation:

Let $T_{h, 0}$ be a quasiuniform mesh on $Ω$ with mesh parameter $h \in (0, h_{0}]$ . Define the active meshes

$T_{h, i} = {T \in T_{h, 0} : T \cap Ω_{i} \neq \emptyset} i = 1, 2, T_{h, Γ} = {T \in T_{h, 0} : T \cap Γ \neq \emptyset}$ (3.1)

associated with the bulk domains $Ω_{i}$ , $i = 1, 2,$ and interface $Γ$ , and the domains covered by the meshes

$O_{h, i} = \cup_{T \in T_{h, i}} i = 1, 2, O_{h, Γ} = \cup_{T \in T_{h, Γ}}$ (3.2)
Let $T_{h, i} (Γ) = {T \in T_{h, i} : T \cap Γ \neq \emptyset}$ and define $F_{h, i}$ as the set of all interior faces associated with an element in $T_{h, i} (\partial Ω_{i})$ .
Let $F_{h, Γ}$ be the set of all interior faces in $T_{h, Γ}$ and $K_{h, Γ} = {K = T \cap Γ : T \in T_{h, Γ}}$ .
Let $V_{h, 0}$ be the space of continuous piecewise linear functions on $T_{h, 0}$ and define

$V_{h, i} = V_{h, 0} |_{T_{h, i}} i = 1, 2, V_{h, Γ} = V_{h, 0} |_{T_{h, Γ}}$ (3.3)

and

$V_{h} = V_{h, 1} \oplus V_{h, 2} \oplus V_{h, Γ}$ (3.4)

3.2 Standard Formulation

The standard finite element method takes the form: find $u_{h} = (u_{h, 1}, u_{h, 2}, u_{h, Γ}) \in V_{h} = V_{h, 1} \oplus V_{h, 2} \oplus V_{h, Γ}$ such that

A_{h}^{S} (u_{h}, v) = L (v) \forall v \in V_{h}

(3.5)

Here the form $A_{h}^{S}$ is defined by

A_{h}^{S} = A + s_{h}

(3.6)

where $s_{h}$ is a stabilization term of the form

s_{h} = s_{h, 1} + s_{h, 2} + s_{h, Γ}

(3.7)

with

s_{h, i} (v, w) = \sum F \in F_{h, i} h_{F} ∥ ζ (A_{i}) ∥_{\infty, F} (⟦ n \cdot \nabla v ⟧, ⟦ n \nabla w ⟧)_{F}, i = 1, 2

where $ζ (X)$ denotes the maximum eigenvalue of the matrix $X$ ,

	$s_{h, Γ} (v, w)$	$= \sum F \in F_{h, Γ} h_{F} ∥ ζ (A_{Γ}) ∥_{\infty, F \cap Γ} (⟦ n \cdot \nabla v ⟧, ⟦ n \nabla w ⟧)_{F_{h, Γ}}$
		$+ \sum T \in T_{h, Γ} h_{K}^{2} ∥ ζ (A_{Γ}) ∥_{\infty, K \cap Γ} (n_{Γ} \cdot \nabla v, n_{Γ} \cdot \nabla w)_{T \cap Γ} .$

3.2.1 Properties of the Stabilization Terms

The rationale for the design of the stabilizing terms is that they improve the stability, while remaining consistent for sufficiently smooth solutions.

Accuracy relies on the following consistency property that is immediate from the definitions above. For any function $v \in H^{\frac{3}{2} + ϵ} (O_{h, i})$ there holds $s_{h, i} (v, w) = 0$ for all $w \in V_{h, i} + H^{\frac{3}{2} + ϵ} (O_{h, i})$ , $i = 1, 2$ . For any function $v \in H^{\frac{3}{2} + ϵ} (O_{h, Γ})$ , such that $n_{Γ} \cdot \nabla v = 0$ on $Γ$ there holds $s_{h, Γ} (v, w) = 0$ for all $w \in V_{h, Γ} + H^{\frac{3}{2} + ϵ} (O_{h, Γ})$ .

The stability properties are well known and we collect them in the following Lemma.

Lemma 3.1

There are constants such that

∥ \nabla v ∥_{A_{i}, O_{h, i}}^{2} ≲ ∥ \nabla v ∥_{A_{i}, Ω_{i}}^{2} + ∥ v ∥_{s_{h, i}}^{2} i = 1, 2

(3.8)

and

∥ \nabla_{Γ} v ∥_{A_{Γ}, O_{h, Γ}}^{2} ≲ ∥ \nabla_{Γ} v ∥_{A_{Γ}, Γ}^{2} + ∥ v ∥_{s_{h, Γ}}^{2}

(3.9)

where we introduced the (semi) norm $∥ v ∥_{s_{h}}^{2} = s_{h} (v, v)$ .

Proof. See [5], [6], and [23], with minor modifications to account for the varying coefficients.

Remark 3.1

Observe that the hidden constants in Lemma 3.1 depend on the variation of the $A_{i}$ and $A_{Γ}$ .

3.3 Robust Formulation

The stabilizing terms ensure robustness irrespective of the intersection of the fracture and the mesh. They do not counter instabilities due to degenerate $B$ . Our aim is to design a formulation which is robust in the case when the eigenvalues of $B$ degenerate. Indeed as we saw above as $ξ$ approaches $1 / 2$ , $λ_{1}$ blows up. For clarity we recall the abstract boundary condition

n \cdot A \nabla v + B (v - v_{Γ}) = 0

(3.10)

where we now assume that the matrix $B$ is a positive definite symmetric $2 \times 2$ matrix with eigenvalues $λ_{i}$ and eigenvectors $e_{i}$ , such that $λ_{i} \in (0, \infty)$ and thus one or both eigenvalues may become very large or small. To handle this situation we instead enforce

B^{- 1} n \cdot A \nabla v + (v - v_{Γ}) = 0

(3.11)

weakly using a modified Nitsche method. This approach was originally developed in [21] where fitted finite element approximation of Robin conditions were considered.

Derivation of an Alternative Weak Form.

As before we have the identity

	$L (v)$	$= 2 \sum i = 1 (A_{i} \nabla u_{i}, \nabla v_{i})_{Ω_{i}} + (A_{Γ} \nabla_{Γ} u_{Γ}, \nabla_{Γ} v_{Γ})_{Γ}      =: A_{1} (u, v) - (n \cdot A \nabla u, v - v_{Γ})_{Γ}$		(3.12)
		$= A_{1} (u, v) - (n \cdot A \nabla u, v - v_{Γ})_{Γ}$		(3.13)

where we introduced the bilinear form $A_{1}$ for brevity. To enforce the interface conditions we proceed as follows

	$L (v)$	$= A_{1} (u, v) - (n \cdot A \nabla u, v - v_{Γ})_{Γ}$
		$= A_{1} (u, v) + (n \cdot A \nabla u, B^{- 1} (n \cdot A \nabla v))_{Γ}$
		$- (n \cdot A \nabla u, B^{- 1} (n \cdot A \nabla v) + (v - v_{Γ}))_{Γ}$
		$= A_{1} (u, v) + (n \cdot A \nabla u, B^{- 1} (n \cdot A \nabla v))_{Γ}$
		$- (n \cdot A \nabla u, B^{- 1} (n \cdot A \nabla v) + (v - v_{Γ}))_{Γ}$
		$- (B^{- 1} (n \cdot A \nabla u) + (u - u_{Γ}), n \cdot A \nabla v)_{Γ}$
		$+ (B^{- 1} (n \cdot A \nabla u) + (u - u_{Γ}), τ (B^{- 1} (n \cdot A \nabla v) + (v - v_{Γ})))_{Γ}$

where the last two terms are zero due to the interface condition and the resulting form on the right hand side is symmetric. Furthermore, $τ$ is a stabilization parameter (a $2 \times 2$ matrix) of the form

τ = 2 \sum i = 1 τ_{i} e_{i} \otimes e_{i}, τ_{i} = \frac{λ_{i} β}{λ_{i} h + β} i = 1, 2

(3.14)

where $β$ is a positive parameter and we recall that $λ_{i}$ and $e_{i}$ are the eigenvalues and eigenvectors of $B$ . The parameter $β$ is chosen so that

∥ n ∥_{A, \infty, Γ}^{2} := 2 \sum i = 1 ∥ n_{i} ∥_{A_{i}, \infty, Γ}^{2} ≲ β

(3.15)

where $∥ w ∥_{A_{i}, \infty, Γ} := ∥ A_{i}^{\frac{1}{2}} w ∥_{\infty, ω}$ is the $A_{i}$ weighted $L^{\infty}$ norm over $Γ$ .

Remark 3.2

The choice of $τ_{i}$ can be further refined as follows

τ_{i} = \frac{λ_{i} β_{i}}{λ_{i} h + β_{i}} i = 1, 2

(3.16)

with

2 \sum j = 1 ∥ n_{j} ∥_{A_{j}, \infty, Γ}^{2} | e_{i j} |^{2} ≲ β_{i}

(3.17)

where $e_{i} = [e_{i 1} e_{i 2}]^{T}$ . This approach is beneficial in situations where the components of $e_{i}$ are very different and there is a large difference between the $∥ n_{j} ∥_{A_{j}, \infty, Γ}^{2}$ with $j = 1$ and $j = 2$ .

The Robust Finite Element Method.

Find $u_{h} \in V_{h}$ such that

A_{h}^{R} (u_{h}, v) := A^{R} (u_{h}, v) + s_{h} (u_{h}, v) = L (v) \forall v \in V_{h}

(3.18)

where

$A^{R} (v, w)$	$= A_{1} (v, w) + (n \cdot A \nabla v, B^{- 1} (n \cdot A \nabla w))_{Γ}$	(3.19)
	$- (n \cdot A \nabla v, B^{- 1} (n \cdot A \nabla w) + (w - w_{Γ}))_{Γ}$	(3.20)
	$- (B^{- 1} (n \cdot A \nabla v) + (v - v_{Γ}), n \cdot A \nabla w)_{Γ}$	(3.21)
	$+ (B^{- 1} (n \cdot A \nabla v) + (v - v_{Γ}), τ (B^{- 1} (n \cdot A \nabla w) + (w - w_{Γ})))_{Γ} .$	(3.22)

It follows by the design of $A^{R}$ that for a sufficiently smooth exact solution $u \in ~ V$ of the problem (2.26) there holds

A (u, v) = A^{R} (u, v) = L (v), \forall v \in (V \cap H^{2} (Ω_{1} \cup Ω_{2} \cup Γ) + V_{h} .

(3.23)

As a consequence we immediately get the Galerkin orthogonality

Lemma 3.2

Let $u \in ~ V$ be the solution of (2.26) and $u_{h} \in V_{h}$ the solution of (3.18) then there holds

A^{R} (u - u_{h}, v) = s_{h} (u_{h}, v) \forall v \in V_{h} .

(3.24)

Proof. The proof follows by combining (3.23) and (3.18).

4 Error Estimates

4.1 The Energy Norm

We introduce the energy norm

⫴ v ⫴_{h}^{2} = 2 \sum i = 1 ∥ \nabla v_{i} ∥_{A_{i}, Ω_{i}}^{2} + h ∥ \nabla v_{i} ∥_{A_{i}, Γ}^{2} + ∥ v ∥_{s_{h}}^{2} + ∥ \nabla_{Γ} v_{Γ} ∥_{A_{Γ}, Γ}^{2} + ∥ v - v_{Γ} ∥_{τ, Γ}^{2}

(4.1)

where $∥ w ∥_{ψ, ω}^{2} = \int_{ω} ψ w^{2}$ is the $ψ$ weighted $L^{2}$ norm over the set $ω$ .

4.2 Interpolation Error Estimates

We begin by introducing the interpolation operators and derive the basic approximation error estimates. Then collecting the estimates we show an interpolation error estimate in the energy norm (

4.1). Since the stabilization operator acts on the finite element solution outside its physical domain of definition, we must make sense of the solution it approximates also outside its physical domain of definition. We will show below show how this can be done using extensions from the physical geometry.

The Scott-Zhang Interpolant.

Given a mesh $T_{h}$ covering a domain $O_{h}$ and the space of piecewise linear continuous finite elements $W_{h}$ , the standard Scott-Zhang interpolation operator $π_{h, S Z} : H^{1} (Ω_{h}) \to W_{h}$ satisfies the element wise estimate

∥ v - π_{h, i, S Z} v ∥_{H^{m} (T)} ≲ h^{2 - m} ∥ v ∥_{H^{2} (N (T))}, m = 0, 1

(4.2)

where $N (T)$ is the set of all elements in $T_{h, i}$ that share a node with $T$ . Note also that the Scott-Zhang interpolant preserves homogeneous boundary conditions exactly. See [26] for further details.

Bulk Domain Fields.

It is shown in [27, Section 2.3, Theorem 5] that there is an extension operator $E_{i} : H^{s} (Ω_{i}) \to H^{s} (R^{d})$ , not dependent on $s \geq 0$ , which is stable in the sense that

∥ E_{i} v_{i} ∥_{H^{s} (R^{d})} ≲ ∥ v_{i} ∥_{H^{s} (Ω_{i})}

(4.3)

We define the interpolation operator $π_{h, i} : H^{1} (Ω_{i}) \to V_{h, i}$ by

π_{h, i} v_{i} = π_{h, i, S Z} E_{i} v

(4.4)

where $π_{h, i, S Z} : H^{1} (O_{h, i}) \to V_{h, i}$ is the Scott-Zhang interpolant and we recall that $O_{h, i} = \cup_{T \in T_{h, i}} T$ is the domain covered by $T_{h, i}$ . We then have the error estimate

∥ v_{i} - π_{h, i} v ∥_{H^{m} (Ω_{i})} ≲ h^{2 - m} ∥ v_{i} ∥_{H^{2} (Ω_{i})} m = 0, 1

(4.5)

Proof. Using the notation $ρ_{i} = v_{i} - π_{h, i} v_{i}$ we obtain

∥ ρ_{i} ∥_{H^{m} (Ω_{i})} ≲ ∥ ρ_{i} ∥_{H^{m} (O_{h, i})} ≲ h^{2 - m} ∥ E_{i} u_{i} ∥_{H^{2} (O_{h, i})} ≲ h^{2 - m} ∥ u_{i} ∥_{H^{2} (Ω_{i})}^{2}

where we used the fact that $Ω_{i} \subset O_{h, i}$ , the interpolation error estimate (4.2), and finally the stability (4.3) of the extension operator $E_{i}$ .

Interface Field.

Let $p_{Γ} : U_{δ} (Γ) \to Γ$ be the closest point mapping from the tubular neighborhood $U_{δ} (Γ) := {x : dist (x, Γ) < δ}$ to $Γ$ , which is well defined for all $δ \in (0, δ_{0}]$ for some $δ_{0} > 0$ . Define the extension operator $E_{Γ} : L^{2} (Γ) \to L^{2} (U_{δ} (Γ))$ by $E_{Γ} v = v \circ p_{Γ}$ . Since $Γ$ is smooth we have the stability estimate

∥ E_{Γ} v_{Γ} ∥_{H^{s} (U_{δ} (Γ))} ≲ δ^{1 / 2} ∥ v_{Γ} ∥_{H^{s} (Γ)} .

(4.6)

Observe also that since $n_{Γ} \cdot \nabla E_{Γ} v_{Γ} = 0$ by construction, and then assuming $s > 3 / 2$ in (4.6) we see that

s_{h, Γ} (E_{Γ} v_{Γ}, w) = 0, \forall w \in V_{h, Γ} + H^{\frac{3}{2} + ϵ} (O_{h, Γ}) .

(4.7)

To define the interpolant we first let

T_{h, δ, Γ} = {T \in T_{h, 0} : T \cap U_{δ} (Γ) \neq \emptyset}, O_{h, δ, Γ} = \cup_{T \in T_{h, δ, Γ}} T

(4.8)

for $δ \in (0, δ_{0} / 2]$ . Then $O_{h, Γ} \subset O_{h, δ, Γ}$ and there are $δ, δ^{'} \in (0, δ_{0}]$ such that $δ \sim δ^{'} \sim h$ and

U_{δ} (Γ) \subset O_{h, δ, Γ} \subset U_{δ^{'}} (Γ) \subset U_{δ_{0}} (Γ)

(4.9)

for all $h \in (0, h_{0}]$ , with $h_{0}$ small enough. We let $V_{h, δ, Γ} = V_{h, 0} |_{O_{h, δ, Γ}}$ and define $π_{h, Γ} : H^{1} (Γ) \to V_{h, Γ}$ by

π_{h, Γ} v_{Γ} = (π_{h, δ, Γ, S Z} E_{Γ} v_{Γ}) |_{O_{h, Γ}}

(4.10)

where $δ \sim h$ and $π_{h, δ, Γ, S Z} : H^{1} (O_{h, δ, Γ}) \to V_{h, δ, Γ}$ is the Scott-Zhang interpolant. We have the error estimate

∥ v - π_{h, Γ} v ∥_{H^{m} (Γ)} ≲ h^{2 - m} ∥ v ∥_{H^{2} (Γ)} m = 0, 1

(4.11)

Proof. Using the trace inequality

∥ v ∥_{Γ}^{2} ≲ δ^{- 1} ∥ v ∥_{U_{δ} (Γ)}^{2} + δ ∥ \nabla v ∥_{U_{δ} (Γ)}^{2} v \in H^{1} (U_{δ} (Γ))

(4.12)

where the hidden constant is independent of $δ$ , we obtain

	$∥ \nabla_{Γ}^{m} ρ ∥_{Γ}^{2}$	$≲ δ^{- 1} ∥ \nabla^{m} ρ ∥_{U_{δ} (Γ)}^{2} + δ ∥ \nabla^{m + 1} ρ ∥_{U_{δ} (Γ)}^{2}$
		$≲ δ^{- 1} ∥ \nabla^{m} ρ ∥_{O_{h, δ, Γ}}^{2} + δ ∥ \nabla^{m + 1} ρ ∥_{O_{h, δ, Γ}}^{2}$
		$≲ δ^{- 1} h^{2 (2 - m)} ∥ \nabla^{2} E_{Γ} v ∥_{O_{h, δ, Γ}}^{2} + δ h^{2 (1 - m)} ∥ \nabla^{2} E$

A cut finite element method for a model of pressure in fractured media

Authors

An adaptive high-order unfitted finite element method for elliptic interface problems

A novel hierarchical multiresolution framework using CutFEM

A cut finite element method for the Darcy problem

A stable and optimally convergent LaTIn-Cut Finite Element Method for multiple unilateral contact problems

Generate plane quad mesh with neural networks and tree search

A Nitsche-based cut finite element method for the coupling of incompressible fluid flow with poroelasticity

Monolithic simulation of convection-coupled phase-change - verification and reproducibility

1 Introduction

2 The Model Problem

2.1 Governing Equations

Remark 2.1

Remark 2.2

Remark 2.3

2.2 Weak Form

2.3 Existence and Uniqueness

3 A Robust Finite Element Method

3.1 The Mesh and Finite Element Spaces

3.2 Standard Formulation

3.2.1 Properties of the Stabilization Terms

Lemma 3.1

Remark 3.1

3.3 Robust Formulation

Derivation of an Alternative Weak Form.

Remark 3.2

The Robust Finite Element Method.

Lemma 3.2

4 Error Estimates

4.1 The Energy Norm

4.2 Interpolation Error Estimates

The Scott-Zhang Interpolant.

Bulk Domain Fields.

Interface Field.

A cut finite element method for a model of pressure in fractured media

Authors

Related Research

An adaptive high-order unfitted finite element method for elliptic interface problems

A novel hierarchical multiresolution framework using CutFEM

A cut finite element method for the Darcy problem

A stable and optimally convergent LaTIn-Cut Finite Element Method for multiple unilateral contact problems

Generate plane quad mesh with neural networks and tree search

A Nitsche-based cut finite element method for the coupling of incompressible fluid flow with poroelasticity

Monolithic simulation of convection-coupled phase-change - verification and reproducibility

This week in AI

1 Introduction

2 The Model Problem

2.1 Governing Equations

Remark 2.1

Remark 2.2

Remark 2.3

2.2 Weak Form

2.3 Existence and Uniqueness

3 A Robust Finite Element Method

3.1 The Mesh and Finite Element Spaces

3.2 Standard Formulation

3.2.1 Properties of the Stabilization Terms

Lemma 3.1

Remark 3.1

3.3 Robust Formulation

Derivation of an Alternative Weak Form.

Remark 3.2

The Robust Finite Element Method.

Lemma 3.2

4 Error Estimates

4.1 The Energy Norm

4.2 Interpolation Error Estimates

The Scott-Zhang Interpolant.

Bulk Domain Fields.

Interface Field.