A cut finite element method for elliptic bulk problems with embedded surfaces

We propose an unfitted finite element method for flow in fractured porous media. The coupling across the fracture uses a Nitsche type mortaring, allowing for an accurate representation of the jump in the normal component of the gradient of the discrete solution across the fracture. The flow field in the fracture is modelled simultaneously, using the average of traces of the bulk variables on the fractures. In particular the Laplace–Beltrami operator for the transport in the fracture is included using the average of the projection on the tangential plane of the fracture of the trace of the bulk gradient. Optimal order error estimates are proven under suitable regularity assumptions on the domain geometry. The extension to the case of bifurcating fractures is discussed. Finally the theory is illustrated by a series of numerical examples.


The model problem
In this section we introduce our model problem. First we present the strong form of the equations and then we derive the weak form that is used for the finite element modelling. We discuss the regularity properties of the solution and show that if the fracture is sufficiently smooth the problem solution, restricted to the subdomains partitioning the global domain, has a regularity that allows for optimal approximation estimates for piecewise affine finite element methods (Fig. 1).

Strong and weak formulations
Let Ω be a convex polygonal domain in R d , with d = 2 or 3. Let Γ be a smooth embedded interface in Ω, which partitions Ω into two subdomains Ω 1 and Ω 2 . We consider the problem: find the pressure u : Ω → R such that Here [v] = v 1 − v 2 , n · a∇v = n 1 · a 1 ∇v 1 + n 2 · a 2 ∇v 2 (2.5) where v i = v| H 1 (Ω i ) , n i is the exterior unit normal to Ω i , a i are positive bounded permeability coefficients, for simplicity taken as constant, and 0 ≤ a Γ < ∞ is a constant permeability coefficient on the interface. Note that it follows from (2.3) that v is continuous across Γ while from (2.2) we conclude that the normal flux is in general not continuous across Γ . Note also that taking a Γ = 0 and f Γ = 0 corresponds to a standard Poisson problem with possible jump in permeability coefficient across Γ .
To derive the weak formulation of the system we introduce the L 2 -scalar product over a domain X ⊂ R d , or X ⊂ R d−1 . For u, v ∈ L 2 (X ) let (u, v) with the associated norm u X = (u, u) integrating by parts over Ω i , and using (2.2), we obtain Γ (2.11) We thus arrive at the weak formulation: find u ∈ V such that Observing that V is a Hilbert space with scalar product a(v, w) = (a∇v, ∇w) Ω + (a Γ ∇ Γ v, ∇ Γ w) Γ (2.13) and associated norm v 2 a = a (v, v) it follows from the Lax-Milgram lemma that there is a unique solution to (2.12) in V for f ∈ H −1 (Ω) and f Γ ∈ H −1 (Γ ).

Regularity properties
To prove optimality of our finite element method we need that the exact solution is sufficiently regular. However since the normal fluxes jump over the interface the solution can not have square integrable weak second derivatives. If the interface is smooth however we will prove that the solution restricted to the different subdomains Ω 1 , Ω 2 and Γ is regular. The upshot of the unfitted finite element method is that this local regularity is sufficient for optimal order approximation. More precisely we have the elliptic regularity estimate Using (2.16) we conclude that Furthermore, using that u Γ ∈ H 1 (Γ ), which follows from the fact that u Γ ∈ V it follows that u Γ | Ω i ∈ H 3/2 (Ω i ), i = 1, 2, and thus n · a∇u Γ ∈ H 1/2 (Γ ) (2.21) Since the right hand side of (2.18) is in L 2 (Γ ) we may use elliptic regularity for the Laplace Beltrami operator to confirm that Collecting the bounds we obtain the refined regularity estimate where we note that we have stronger control of u Γ on the subdomains.

The mesh and finite element space
Let T h be a quasi-uniform conformal mesh, consisting of shape regular elements with mesh parameter h ∈ (0, h 0 ], on Ω and let be the active meshes associated with Ω i , i = 1, 2. Let V h be a finite element space consisting of continuous piecewise polynomials on T h and define In general, we simplify the notation and write v = v. Finally, we use T h (Γ ) to denote the set of elements intersected by Γ .

Derivation of the method
To derive the finite element method we follow the same approach as when introducing the weak formulation, but taking care to handle the boundary integrals that appear due to the discontinuities in the approximation space.
Testing the exact problem with v ∈ W h and integrating by parts over Ω 1 and Ω 2 we find that GEM -International Journal on Geomathematics (2019) 10:10 where in the last identity we symmetrized using the fact that [u] = 0. We also used the identity where the averages are defined by w = κ 1 w 1 + κ 2 w 2 , w * = κ 2 w 1 + κ 1 w 2 (3.12) with κ 1 + κ 2 = 1 and 0 ≤ κ i ≤ 1.
Introducing the bilinear forms (3.14) the above formal derivation leads to the following consistent formulation for discontinuous test functions w.
Observe that we have modified a h,Γ by introducing the average v * also in the left factor. This changes nothing when applied to a smooth solution, but will allow also to apply the form to the discontinuous discrete approximation space. The subscript h in the form indicates that it is modified to be well defined for the discontinuous approximation space. The definition of W is motivated by the fact that the trace terms should be well defined, for instance,

The finite element method
The finite element method that we propose is based on the formulation (3.16). However, using this formulation as it stands does not lead to a robust approximation method. Indeed we need to ensure stability of the formulation through the addition of consistent penalty terms. First we need to enforce continuity of the discrete solution across Γ . To this end we introduce an augmented version of a Ω , u, w). Secondly, to obtain stability independently of how the interface cuts the computational mesh and for strongly varying permeabilities a 1 , a 2 and a Γ we also need some penalty terms in a neighbourhood of the interface. We define where γ is a positive parameters and F h,i is the set of interior faces in T h,i that belongs to an element T ∈ T h,i which intersects Γ , see Fig. 6. Observe that for Collecting the above bilinear forms in the finite element method reads:

Analysis of the method
In this section we derive the basic error estimates that the solution of the formulation (3.21) satisfies. The technical detail is kept to a minimum to improve readability. In particular, we assume that the bilinear forms can be computed exactly and that Γ fulfils the conditions of Hansbo and Hansbo (2002). For a more complete exposition in a similar context we refer to Burman et al. (2016).

Properties of the bilinear form
For the analysis it is convenient to define the following energy norm where |v| s i = s i (v, v) 1/2 and c a is a constant fulfilling c a ∼ a −1 min , cf. Lemma 4.1 below. (3.20), satisfies the following bounds:

Lemma 4.1 The form A h , defined in
where W was introduced in (3.16).
provided β is large enough.
Proof The first estimate (4.2) follows directly from the Cauchy-Schwarz inequality.
To show (4.3) we recall the following inequalities: Hansbo and Hansbo 2002) (4.5) In (4.5) we used the notation Using (4.4) and (4.5) we obtain the following bound on the fluxes Now assume that κ i a i ≤ a min := min i∈{1,2} a i , for instance one may take κ 1 = a 2 /(a 1 + a 2 ) and κ 2 = a 1 /(a 1 + a 2 ) then It follows that The bound (4.3) now follows taking ε = 1/(2C) and β > 2Ca min and by applying once again (4.7), taking c a ∼ a −1 min .
A consequence of the bound (4.3) is the existence of a unique solution to (3.21).

Lemma 4.2 The linear system defined by the formulation (3.21) is invertible.
Proof Follows from Lax-Milgram's lemma.

Interpolation
For δ > 0 let E i : H s (Ω i ) → H s (Ω) be a continuous extension operator s > 0. We define the interpolation operator , 2, and π S Z h is the Scott-Zhang interpolation operator. We then have the interpolation error estimate where, with ρ Γ the signed distance function associated with Γ , and Proof To prove the estimate (4.13) we use a trace-inequality on functions in H 1 (T h (Γ )) (i.e., with · T h (Γ ) the broken H 1 -norm over the elements intersected by GEM -International Journal on Geomathematics (2019) 10:10 see Hansbo and Hansbo (2002), then interpolation on T h (Γ ) and finally we use the stability of the extension operator E i . First observe that by using the trace inequality (4.16) we obtain, with v = u − π h u Using standard interpolation for the Scott-Zhang interpolation operator we get the bound where we used the stability of the extension operator in the last inequality. The bound follows similarly using element-wise trace inequalities follows by interpolation (c.f. Burman and Hansbo 2012). The interpolation error estimate for the terms due to the Laplace-Beltrami operator on Γ is a bit more delicate. We use a trace inequality to conclude that (4.23) Observing that we may take δ ∼ h the estimate follows.
Comparing (4.13) with (2.23) we see that we have a small mismatch between the regularity that we can prove and that required to achieve optimal convergence. In view of this we need to assume a slightly more regular solution for the H 1 -error estimates below. The sub optimal regularity also interferes in the L 2 -error estimates. Here we need to use (2.23) on the dual solution and in this case the additional regularity of the estimate (4.13) is not available. Instead we need to find the largest ζ ∈ [0, 1] such that |||u − π h u||| h h ζ 2 i=1 u H 2 (Ω i ) , which will result in a suboptimality by a power of 1 − ζ in the convergence order in the L 2 -norm. Revisiting the analysis above up to (4.22) we see that

Theorem 4.1 If u is the solution to
, and u h is the finite element approximation defined by (3.21), then (4.26) (4.27) Proof (4.26). Splitting the error and using the interpolation error estimate we have Using coercivity (4.3), Galerkin orthogonality and continuity (4.2) the second term can be estimated as follows and thus applying the approximation result (4.13) we conclude that (4.33) (4.27). Consider the dual problem Setting v = e = u − u h and using Galerkin orthogonality, followed by the continuity (4.2) and the suboptimal approximation estimate (4.24) on |||φ − π h φ||| h we get In the last step we used the elliptic regularity estimate (4.35) for the dual problem. Setting ψ i = e i / e i Ω i and ψ Γ = e Γ / e Γ Γ estimate (4.27) follows.

Remark 4.1
As noted before the error estimate in the L 2 -norm is suboptimal with a power 1/2. To improve on this estimate we would need to sharpen the regularities required for the approximation estimate (4.13). This appears to be highly non-trivial since the interpolation of u and u Γ can not be separated when both are interpolated using the bulk unknowns. Therefore we did not manage to exploit the stronger control that we have on the harmonic extension of u Γ in (2.23). Note however that if separate fields are used on the fracture and in the bulk domains we would recover optimal order convergence in L 2 .

Remark 4.2
Using the stronger control of the regularity of the harmonic extension provided by (2.23) we may however establish an optimal order L 2 error estimate for the solution on Γ , (2.14)

Extension to bifurcating fractures
In the case most common in applications, fractures bifurcate, leading to networks of interfaces in the bulk. It is straightforward to include this case in the method above and we will discuss the method with bifurcating fractures below. The analysis can also be extended under suitable regularity assumptions, but becomes increasingly technical. We leave the analysis of the methods modelling flow in fractured media with bifurcating interfaces to future work.

The model problem
Description of the domain Let us for simplicity consider a two dimensional problem with a one dimensional interface. We define the following: -Let the interface Γ be described as a planar graph with nodes N = {x i } i∈I N and edges E = {Γ j } j∈I E , where I N , I E are finite index sets, and each Γ j is a smooth curve between two nodes with indexes I N ( j). Note that edges only meet in nodes. -For each i ∈ I N we let I E (i) be the set of indexes corresponding to edges for which x i is a node. For each i ∈ I N we let I E (i) be the set of indexes j such that The Kirchhoff condition The governing equations are given by (2.1)-(2.4) together with two conditions at each of the nodes x i ∈ N , the continuity condition ∀k, l ∈ I E (i) (5.1) and the Kirchhoff condition where t Γ j (x i ) is the exterior tangent unit vector to Γ j at x i . Note that in the special case when a node x i is an end point of only one curve the Kirchhoff condition becomes a homogeneous Neumann condition.

Forms associated with the bifurcating interface
We proceed as in the derivation (2.8)-(2.11) of the weak problem (2.12) in the standard case. However, when we use Green's formula on Γ we proceed segment by segment as follows where we changed the order of summation and used the Kirchhoff condition (5.2) to subtract the nodal average where 0 < κ Γ i , and j∈I E (i) κ Γ j = 1. Note that when a node x i is an end point of only one curve the contribution from x i is zero, because in that case we have v * i | x i − v * = 0 since there is only one element in I E (i), and thus we get the standard weak enforcement of the homogeneous Neumann condition.
Symmetrizing and adding a penalty term we obtain the form where β Γ is a stabilisation parameter with the same function as β. A similar derivation can be performed for a two dimensional bifurcating fracture embedded into R 3 , see Hansbo et al. (2017) for further details.
To ensure coercivity we add a stabilization term of the form where F h (x i ) is the set of interior faces in the patch of elements N h (T (x i )) and T (x i ) is an element such that x i ∈ T .
We finally define the form A h,Γ associated with the bifurcating crack by The method Define 6 Numerical examples

Implementation details
We will employ piecewise linear triangles and extend the implementation approach proposed in Hansbo and Hansbo (2002) to include also bifurcating fractures. Recall that T h (Γ ) denotes the set of elements intersected by Γ , where each side of the intersection belongs to Ω 1 and Ω 2 , respectively. For each element in T i ∈ T h (Γ ), we assign elements T i,1 ∈ T h,1 and T i,2 ∈ T h,2 by overlapping the existing element T i ∈ T h (Γ ) using the same nodes from the original triangulation. Elements T i,1 and T i,2 coincide geometrically, see Fig. 3. To ensure continuity, we used the same process on the neighboring elements and checked if new nodes had already been assigned. For each bifurcation point, two approaches can be adapted. Either by letting the bifurcation point coincide with a node or by the less straight-forward approach to overlap the existing element T i ∈ T h (Γ ) into T i,1 , T i,2 and T i,3 , see Fig. 4. For simplicity of implementation, we have here chosen to let the bifurcating point coincide with a node. The triangles T i / ∈ T h (Γ ) were handled in the usual way. The stabilization (3.19) was only applied to the cut sides of the elements which in all examples was sufficient for stability.

Example 1: No flow in fracture
We consider an example on Ω = (0, 1) × (0, 1), from Hansbo and Hansbo (2002). We solved the example with an added bifurcation point. For the added fracture, we  Fig. 4 The split of a triangle with bifurcation point denote the diffusion coefficient by a Γ 1 . The exact solution is given by where r = x 2 + y 2 . We chose r 0 = 3/4, a 1 = 1, a 2 = 1000 and a Γ = a Γ 1 = 0, with a right-hand side f = − 4 and f Γ = 0. The boundary conditions were symmetry boundaries at x = 0 and y = 0 and Dirichlet boundary conditions corresponding to the exact solution at x = 1 and y = 1. This example is outlined in Figs. 5 and 6. We give the elevation of the approximate solution in Fig. 7, on the last mesh in a sequence. The corresponding convergence of the L 2 -norm and the energy-norm is given in Fig. 8.

Example 2: Flow in the fracture
We considered a two-dimensional example on the domain Ω = 1, e 5/4 × 1, e 5/4 , from Burman et al. (2017). We solved the example with an additional fracture added, see Fig. 9. The exact solution is given by  Fig. 9 Active meshes with two embedded fractures, Example 2 u 1 = log (r ) 5 (4 + e) for 1 < r < e, u 2 = 4 − 4e 5 log (r ) − 5 4 + 1 for e < r < e 5/4 , where x 2 + y 2 := r = e. We chose a 1 = a 2 = a Γ = 1 and the right hand side to f = f Γ = 0. For the added crack we chose a Γ 1 = 0. The Dirichlet boundary conditions corresponding to the exact solution at x, y = 1 and x, y = e 5/4 . In Fig. 10, we give the elevation of the approximate solution. The corresponding L 2 -norm convergence and the energy-norm is given in Fig. 11.

Example 3: Flow in bifurcating fractures
We consider an example with two bifurcating points. The fractures are modeled using higher order curves. In Fig. 12 we show the fractures and construction of individual elements. On the domain Ω = (0, 1) × (0, 1), we chose a 1 = a 2 = 1, f Ω = 1 and f Γ = 0. We impose the Dirichlet boundary conditions u = 0 at x, y = 0 and u = 0 at x, y = 1. For the diffusion coefficient, we denote a Γ i for each fracture where a Γ i ∈ {0, 100} and assign an individual value for each Γ i , see Fig. 13. We report six different solutions by allowing different fractures to be active. The first result have been obtained with a Γ i = 0, thus no flow in the fractures is allowed, see Fig. 14.  figure: a Γ 1 = a Γ 2 = a Γ 3 = a Γ 4 = a Γ 5 = 0, and assigned values to the right figure: a Γ 1 = 100 and a Γ 2 = a Γ 3 = a Γ 4 = a Γ 5 = 0 Fig. 15 Elevation of the approximate solution using two bifurcating points, Example 3. Assigned value to the left figure: a Γ 1 = a Γ 2 = 100 and a Γ 3 = a Γ 4 = a Γ 5 = 0, and assigned values to the right figure: a Γ 1 = a Γ 2 = a Γ 3 = 100 and a Γ 4 = a Γ 5 = 0 Further, each solution has one additional fracture activated by changing the diffusion coefficent a Γ i = 100, see Figs. 14, 15 and 16. The last solution is presented with flow in all fractures.

Concluding remarks
We proposed a discontinuous finite element method using a one-field approach to modelling Darcy flow in a cracked medium. The pressure in the crack was modelled as an average of pressure on either side of the crack which, unlike our previous work (Burman et al. 2017), allows for pressure jumps across the crack. In particular, the case of bifurcating fractures was considered. Optimal order error estimates were proven and backed up by numerical experiments. Extension to other flow models in the crack have been considered in Burman et al. (2018).