Stable Mixed Finite Elements for Linear Elasticity with Thin Inclusions

TL;DR

This paper develops and analyzes stable mixed finite element methods for linear elasticity problems involving complex geometries with thin inclusions modeled as mixed-dimensional manifolds, ensuring stability and convergence.

Contribution

It introduces a novel mixed-dimensional formulation for elasticity with thin inclusions and proposes stable finite element discretizations with proven stability and convergence.

Findings

The mixed-dimensional system is well-posed with appropriate norms.

Finite element schemes conserve linear momentum locally.

Error estimates demonstrate stability and convergence.

Abstract

We consider mechanics of composite materials in which thin inclusions are modeled by lower-dimensional manifolds. By successively applying the dimensional reduction to junctions and intersections within the material, a geometry of hierarchically connected manifolds is formed which we refer to as mixed-dimensional. The governing equations with respect to linear elasticity are then defined on this mixed-dimensional geometry. The resulting system of partial differential equations is also referred to as mixed-dimensional, since functions defined on domains of multiple dimensionalities are considered in a fully coupled manner. With the use of a semi-discrete differential operator, we obtain the variational formulation of this system in terms of both displacements and stresses. The system is then analyzed and shown to be well-posed with respect to appropriately weighted norms. Numerical discretization schemes are proposed using well-known mixed finite elements in all dimensions. The schemes conserve linear momentum locally while relaxing the symmetry condition on the stress tensor. Stability and convergence are shown using a priori error estimates.

Tables

Table 1. Table 1: The finite element spaces chosen for n = 2 𝑛 2 n=2 and n = 3 𝑛 3 n=3 of lowest order within the family ( 6.1 ). The negative orders in the subscript denote discontinuous Lagrange elements. On the zero-dimensional manifolds, the polynomial order is redundant since any finite element space corresponds to point evaluation there.

$d$	$W_h$	${\Sigma }_h$	${𝑼}_h$	$R_h$
2	${\left(P_2\right)}^2$	${\left(BD{M}_1\right)}^2$	${\left(P_0\right)}^2$	$P_0$
1		${\left(P_2\right)}^2$	${\left(P_{-1}\right)}^2$
0			${\left(P_{-2}\right)}^2$

Table 2. Table 2: The lowest-order finite element spaces of the reduced family for n = 2 𝑛 2 n=2 and n = 3 𝑛 3 n=3 . The superscript minus indicates that the trace of the finite element space onto Γ h subscript Γ ℎ \Gamma_{h} is reduced by one order.

$d$	$W_h$	${\Sigma }_h$	${𝑼}_h$	$R_h$
2	${\left(P_2^{-}\right)}^2$	${\left(BD{M}_1^{-}\right)}^2$	${\left(P_0\right)}^2$	$P_0$
1		${\left(P_1\right)}^2$	${\left(P_0\right)}^2$
0			${\left(P_0\right)}^2$

Equations

$$I^d$$

$${\hat{J}}_i$$

$${\Gamma }^d$$

$${\partial }_u{\Omega }_i$$

$$\check{f}$$

$$A{\sigma }_i$$

$$A{\sigma }_i=\frac{1}{2\mu }\left({\sigma }_i-\frac{\lambda }{2\mu +n\lambda }\mathrm{Tr}{\sigma }_{i}I\right),$$

$$A{\sigma }_i$$

$$A{\sigma }_i-\nabla {\mathbf{u}}_i+{\chi }_i$$

$$\nabla \cdot {\sigma }_i$$

$$\mathrm{asym}{\sigma }_i$$

$$\mathbf{n}\cdot {\sigma }_i$$

$$V_i$$

$$({\widehat{ϵ}}_{\max }){|}_{{\Omega }_i}$$

$$ϵ$$

$$ϵ{\sigma }_i^{avg}={\sigma }_i={ϵ}^{-1}{\sigma }_i^{int}.$$

$$\sigma$$

$$\nabla \cdot {ϵ}_i{\sigma }_i-\sum _{j\in {\hat{J}}_i}(\mathbf{n}\cdot {ϵ}_{\hat{ȷ}}{\sigma }_{\hat{ȷ}}){|}_{{\Gamma }_j}$$

$$\sum _{j\in {\hat{J}}_i}(\mathbf{n}\cdot {ϵ}_{\hat{ȷ}}{\sigma }_{\hat{ȷ}}){|}_{{\Gamma }_j}$$

$$[[\varphi {]]}_i$$

$$\mathfrak{D}\cdot ϵ\sigma$$

$$\mathfrak{D}\mathbf{u}$$

$$\mathcal{A}({ϵ}^{-1}\sigma )$$

$$\mathfrak{A}\sigma$$

$${\mathfrak{A}}_{\Vert }\sigma$$

$${\chi }_i$$

$${ϵ}^{-1}\sigma$$

$$\mathbf{n}\cdot ({\widehat{ϵ}}^{-1}\hat{\sigma })$$

$$({\mathfrak{A}}_{\Vert }\sigma ){|}_{{\Omega }_i}$$

$$({\mathfrak{A}}_{\perp }\sigma ){|}_{{\Gamma }_j}$$

$$\mathrm{asym}ϵ\sigma$$

$$\mathfrak{A}\sigma -ϵ\mathfrak{D}\mathbf{u}+ϵ\chi$$

$$\mathfrak{D}\cdot ϵ\sigma$$

Full text

The governing equations with respect to linear elasticity are then defined on this mixed-dimensional geometry. The resulting system of partial differential equations is also referred to as mixed-dimensional, since functions defined on domains of multiple dimensionalities are considered in a fully coupled manner. With the use of a semi-discrete differential operator, we obtain the variational formulation of this system in terms of both displacements and stresses. The system is then analyzed and shown to be well-posed with respect to appropriately weighted norms.

Numerical discretization schemes are proposed using well-known mixed finite elements in all dimensions. The schemes conserve linear momentum locally while relaxing the symmetry condition on the stress tensor. Stability and convergence are shown using a priori error estimates.

1 Introduction

Thin inclusions in elastic materials arise in a variety of scientific fields, including geo-physics, bio-mechanics, and the study of composite materials. The subsurface, for example, typically includes rock layers with significantly larger horizontal extent compared to their height. Since it is often infeasible to resolve such small heights for large-scale simulations, we consider the setting where the layer, or aquifer, is represented by a lower-dimensional manifold. The governing equations on this manifold can be derived using vertical integration nordbottenbook , see e.g. bjornaraa2016vertically for an application with respect to CO2 storage.

Secondly, membranes occur frequently in the study of bio-mechanical systems. Examples range from cell walls in plants to the heart sac and dermal layer in human physiology. As a modeling assumption, each of these membranes can be represented by lower-dimensional manifolds. Their influence on the coupled mechanical system can then be incorporated by assigning significantly different material properties compared to the surroundings.

A third application concerns the study of composite or reinforced materials. In this context, the lower-dimensional manifolds correspond to the stiffer plates embedded in the material for strengthening purposes. This can be expanded to connected, two-dimensional objects such as H-beams and T-beams. The junctions are then considered one-dimensional manifolds, with inherited or separately defined material properties. We note that this work is limited to manifolds of codimension one and thus does not treat the case of embedded, one-dimensional rods in three dimensions.

The thin features are considered lower-dimensional and have elastic properties, yet a slightly different setting is presented than the conventional theory of thin shells ciarlet2000theory . The main difference is that we focus on a strong coupling of a thin inclusion with a surrounding, elastic medium. The interest of this work is therefore more closely aligned to caillerie1980effect , in which rigid, thin inclusions are considered.

The structure of the derived equations fits well with the mixed-dimensional framework derived in Nordbotten2017DD ; boon2017excalc . We aim to preserve this structure and retain a local conservation of linear momentum after discretization with the use of conforming, mixed finite elements. The construction of stable finite element pairs representing displacements and symmetric stresses is involved, and typically leads to higher-order elements for the stress space arnold2002mixed . By relaxing the symmetry condition on the stress tensor as in arnold2006differential ; awanou2013rectangular , these difficulties can be mitigated.

This article is structured as follows. Section 2 introduces the notational conventions and the decomposition of the geometry according to dimension. On this geometry, Section 3 introduces the governing equations of the model in all dimensions. We introduce the relevant function spaces and present the derivation of the variational formulation in Section 4. The resulting system of equations is proven to be well-posed in Section 5. Finally, we propose conforming discretization schemes in Section 6 for which stability and convergence are shown.

2 Geometry and Notation

In this section, we introduce the mixed-dimensional geometry and establish notation. Here, we follow the conventions introduced in Boon2018Robust ; boon2017excalc .

Let us consider an $n$-dimensional domain $Y$ that contains thin, embedded structures, represented by lower-dimensional manifolds. In general, we consider $n=3$, the two-dimensional case being simpler. Let $\Omega_{i}^{d_{i}}$ be such a manifold, with $i$ the unique index from a global set $I$ and $d_{i}$ its dimension. The superscript is frequently omitted for brevity. For $1\leq d\leq n-1$, we successively identify intersections between $d$-manifolds as $(d-1)$-manifolds. All $\Omega_{i}$ are open sets and mutually disjoint. For simplicity, we restrict this work to the case in which all $\Omega_{i}$ have zero curvature, i.e. are flat.

As an example, let us consider the two-dimensional set-up in Figure 1 (left). Here, an embedded H-beam is described using two zero-dimensional intersection points and five one-dimensional line segments. The open set corresponding to the surrounding medium is given by $\Omega_{8}^{2}=\operatorname{int}(Y\setminus\cup_{i=1}^{7}\Omega_{i}^{d_{i}})$.

We refer to the union $\bigcup_{i\in I}\Omega_{i}$ as the mixed-dimensional geometry $\Omega$. Let $I^{d}$ be the set of indices corresponding to $d$-manifolds and let $\Omega^{d}$ be the collection of such manifolds. In short, we denote

[TABLE]

The interface between manifolds of codimension one will play an important role, and we adopt a separate notation for these. Let $J$ be the set of indices such that each $j\in J$ corresponds to an interface $\Gamma_{j}$ between $\Omega_{i}$ for some $i\in I$ and an adjacent domain of dimension $(d_{i}+1)$. $\Gamma_{j}$ physically coincides with $\Omega_{i}$ and we assume that a unique ${\hat{\jmath}}\in I$ exists such that $\Gamma_{j}\subseteq\partial\Omega_{\hat{\jmath}}$.

To distinguish these different interfaces, we define the following index sets for $i\in I$:

[TABLE]

An example is shown in Figure 1 (right), which emphasizes that for $j_{1},j_{2}\in\hat{J}_{4}$ with $j_{1}\neq j_{2}$, we have ${\hat{\jmath}}_{1}={\hat{\jmath}}_{2}=8$. In other words, we allow for a manifold $\Omega_{\hat{\jmath}}$ to border on multiple sides of $\Omega_{i}$ and assign a unique index $j\in\hat{J}_{i}$ to each side. Finally, we remark that $\hat{J}_{i}$ is void for all $i\in I^{n}$, by definition.

Using the same summation convention per dimension as above, we denote

[TABLE]

Each $\Gamma_{j}$ is equipped with a unit normal vector $\bm{n}_{j}$, from the tangent space of $\Omega_{\hat{\jmath}}$, oriented outward with respect to $\Omega_{\hat{\jmath}}$. The subscript on $\bm{n}$ is omitted for brevity. In reference to the vector(s) normal to a manifold $\Omega_{i}$, we will use the check notation $\check{\bm{n}}_{i}$.

The boundary of the domain is given by the disjoint union $\partial_{\sigma}Y\bigcup\partial_{u}Y$ on which different boundary conditions will be imposed. In particular, we assume that the displacement is given on $\partial_{u}Y$ and the normal stress on $\partial_{\sigma}Y$. We denote for $i\in I$,

[TABLE]

For analysis purposes, we assume that $|\partial_{u}\Omega_{i}|>0$ for all $i\in I^{n}$, i.e. each subdomain of dimension $n$ is connected to a part of the boundary on which the displacement is prescribed. By omission of the subscript, we use $\partial_{u}\Omega$ and $\partial_{\sigma}\Omega$ to refer to the corresponding boundaries of the mixed-dimensional geometry.

Given a function $f$ defined on the mixed-dimensional geometry $\Omega$, let $f_{i}$ denote its restriction to $\Omega_{i}$, i.e. $f_{i}=(f)|_{\Omega_{i}}$. Furthermore, we employ the hat and check notation to distinguish instances of $f$ inherited from different domains onto the interface $\Gamma$. This means that on $\Gamma_{j}$ with $j\in\hat{J}_{i}$, we denote

[TABLE]

Note that the definition of $\hat{f}$ involves a trace of $f_{\hat{\jmath}}$ onto $\Gamma_{j}$.

3 Model Formulation

In this section, we consider the governing equations and introduce the model problem. Starting with the mathematical formulation of linear elasticity in the surrounding medium, we continue with the generalized equations on lower-dimensional manifolds to derive the strong form of the mechanics problem. The variational formulation is considered afterward in Section 4.

3.1 Governing Equations in the Surrounding Medium

Let us start by presenting the governing equations for linear elasticity in the surrounding medium $\Omega^{n}$. For $i\in I^{n}$, let $\sigma_{i}$ denote the elastic stress and $\bm{u}_{i}$ the displacement. Assuming infinitesimal strain, the stress-strain relationship has the general form:

[TABLE]

In case of homogeneous and isotropic media, the operator $A$ describes Hooke’s law and is given by

[TABLE]

in which $\lambda$ and $\mu$ are the Lamé parameters and Tr is the matrix trace operator.

In the variational formulation presented Section 4, the symmetry of the stress tensor $\sigma_{i}$ will be enforced in a weak sense. In preparation, we introduce the antisymmetric tensor $\chi_{i}={\operatorname{asym}\,}{\nabla\bm{u}_{i}}$ such that:

[TABLE]

With the addition of linear and angular momentum conservation, the following system of equations is formed in each $\Omega_{i}$ with $i\in I^{n}$.

[TABLE]

with $\bm{g}_{u}$ a given function. We limit the exposition to homogeneous stress boundary conditions, noting that this can readily be extended to the general case.

3.2 Geometrical Scaling and Constraints

The governing equations on the lower-dimensional manifolds will be significantly influenced by the small width of the thin inclusions. We therefore devote this section to defining two key parameters, $\gamma$ and $\epsilon$ and the constraints these parameters adhere to.

Let $\gamma$ be a virtual parameter representing the relatively small length from the interface of the higher-dimensional domain to the central plane (2D), line (1D), or point (0D) of the physical inclusion. We assume that on each interface $\Gamma_{j}$, $\gamma$ is constant and positive.

The small value of $\gamma$ will introduce a scaling in the equations. In the formulation of the problem, it is advantageous if this scaling appears in the coupling terms between variables, rather than on the main diagonal of the system. For this purpose, we introduce a key concept from the context of subsurface flow models arbogast2016linear ; Boon2018Robust . As shown in those works, the desired scaling of the system is achieved by employing an appropriately scaled flux variable. In analogy, we employ a scaled stress by introducing the scaling parameter $\epsilon$ on $\Omega$. This parameter is defined as the square root of the (small) cross-sectional length (2D), area (1D), or volume (0D) of the dimensionally reducible feature. $\epsilon$ is assumed to be constant and positive for each manifold $\Omega_{i}$. The following relationship is imposed

[TABLE]

with $V_{i}$ the cross-sectional measure of the feature corresponding to $\Omega_{i}$ and $\epsilon_{i}=1$ for $i\in I^{n}$. The relation $a\eqsim b$ (respectively $\lesssim$ and $\gtrsim$) implies that a constant $C>0$ exists, independent of $\epsilon$, $\gamma$, and the mesh size $h$, such that $a=Cb$ (respectively $\leq$ and $\geq$). Since $\epsilon$ is assumed to be small, we often use the relationship $\epsilon\lesssim 1$.

Next, we relate the different values of $\epsilon$ between the dimensions. For each $\Omega_{i}$, let $\hat{\epsilon}_{\max}$ be the maximal value of $\epsilon$ in the adjacent higher-dimensional manifolds:

[TABLE]

In $\Omega^{n}$, we define $\hat{\epsilon}_{\max}$ as unity. Using this definition, we add a constraint to the geometry by assuming that $\hat{\epsilon}_{\max}$ bounds $\epsilon$ from above, i.e.

[TABLE]

With the parameter $\epsilon$ defined, we continue with the scaling of the stress variable on the lower-dimensional manifolds. In $\Omega_{i}$ with $d_{i}<n$, let $\sigma^{avg}_{i}$ contain the columns of the Cauchy stress tensor associated with the tangent bundle of $\Omega_{i}$, averaged over the cross-section of the physical feature. The integrated stress tensor $\sigma^{int}_{i}$, on the other hand, is obtained after multiplication with the cross-sectional measure $V_{i}$.

Using the factor $\epsilon_{i}:=\sqrt{V_{i}}$, the scaled stress $\sigma_{i}$ on $\Omega_{i}$ is defined as

[TABLE]

The columns and first $d_{i}$ rows of $\sigma_{i}$ correspond to the basis vectors from the tangent bundle of $\Omega_{i}$. The final $(n-d_{i})$ rows relate to the directions normal to the manifold. Thus, $\sigma_{i}\in\mathbb{R}^{n\times d_{i}}$ by definition, making it undefined in the intersection points $\Omega^{0}$. The displacement $\bm{u}_{i}$ in $\Omega_{i}$ remains unscaled and is naturally in $\mathbb{R}^{n}$. Again, we omit the subscript $i$ to refer to the mixed-dimensional entities, i.e.

[TABLE]

We emphasize that the integrated and average stress quantities can readily be recovered from the scaled stress $\sigma$ after appropriate post-processing with the known quantity $\epsilon$.

3.3 Mixed-Dimensional Equations

With the given scaling from the previous subsection, let us consider the governing equations in the lower-dimensional manifolds. In this generalization to the mixed-dimensional geometry, a structure similar to the system (3.2) is uncovered. We start by introducing the linear momentum balance equation, followed by the stress-strain relationships and finish with the conservation of angular momentum.

The balance of linear momentum (3.2b) is generalized first. After integrating the conservation law in the direction(s) normal to the inclusion (see e.g. Boon2018Robust ; Roberts2 for the analogue in fracture flow models), we obtain

[TABLE]

Here, $\bm{f}_{i}$ is the body force acting on $\Omega_{i}$, averaged over the cross-section with measure $\epsilon^{2}$. The in-plane divergence $(\nabla\cdot)$ on $\Omega_{i}$ and the normal trace operator $(\bm{n}\cdot)$ onto $\Gamma_{j}$ are applied row-wise. Hence, the divergence in $\Omega_{i}$ maps from $\mathbb{R}^{n\times d_{i}}$ to $\mathbb{R}^{n}$ and the normal trace on $\Gamma_{j}$ maps $\mathbb{R}^{n\times d_{\hat{\jmath}}}$ to $\mathbb{R}^{n}$.

For the zero-dimensional manifolds, there are no divergence operator or $\sigma_{i}$ available and the balance law is completely given by the sum of forces from $\Omega^{1}$:

[TABLE]

To shorten notation, we introduce the jump operator $\llbracket\cdot\rrbracket_{i}$ which maps functions defined on the interface $\Gamma_{j}$ with $j\in\hat{J}_{i}$ to the central manifold $\Omega_{i}$ such that

[TABLE]

Following Boon2018Robust , we introduce the mixed-dimensional divergence operator ($\mathfrak{D}\cdot$) as $\mathfrak{D}\cdot\sigma:=\nabla\cdot\sigma-\llbracket\bm{n}\cdot\hat{\sigma}\rrbracket$ and rewrite the conservation equation to the concise form

[TABLE]

We now continue by defining the stress-strain relationships in the lower-dimensional manifolds in analogy with (3.2a). For that, we first introduce the gradient operator $\mathfrak{D}$ as

[TABLE]

We emphasize that the gradient $\nabla$ relates to the tangential direction(s) and is applied row-wise. Since we have $\mathfrak{D}\bm{u}$ defined on both $\Omega$ and $\Gamma$, we need to provide stress-strain relationships inside and on the boundaries of the domains.

The stress-strain relationship is then described by an operator $\mathfrak{A}$ acting on the averaged stress. Here, we pay special attention to the scaling with $\epsilon$. Thus, recalling that the averaged stress is denoted by $\sigma^{avg}=\epsilon^{-1}\sigma$, the stress-strain relationships are given by

[TABLE]

with $\chi$ to be defined. To obtain a symmetric system, we scale this equation with $\epsilon$. Noting that $\epsilon$ and $\mathcal{A}$ do not necessarily commute, we introduce $\mathfrak{A}:=\epsilon\mathcal{A}\epsilon^{-1}$ to obtain the generalized version of the stress-strain relationship:

[TABLE]

The restrictions of $\mathfrak{A}$ to the manifolds $\Omega$ and interfaces $\Gamma$ are respectively denoted by

[TABLE]

The variable $\chi$ is the generalization of the asymmetric $\chi_{i}$ from section 3.1, given by

[TABLE]

We interpret $(\chi)|_{\Gamma}=0$ and $\chi_{i}=0$ for $i\in I^{0}\cup I^{1}$.

Example 1

We provide an explicit example of $\mathfrak{A}$ using a fictitious material. In this material, we assume that the stress-strain relationships in tangential and normal directions are independent. This assumption leads to a model which captures in-plane shearing whereas out-of-plane stress components follow a one-dimensional Hooke’s law. This behavior is described by the following constitutive laws

[TABLE]

Here, $\mu$ and $\lambda$ (respectively $\mu_{\perp}$ and $\lambda_{\perp}$) are the Lamé parameters describing the stress-strain relationship tangential (and normal) to the manifold.

The inverse relations, mapping stresses to strains, are then given by

[TABLE]

Here, $I_{d}$ is the identity tensor in $\mathbb{R}^{d\times d}$. We remark that in this example, we have $\epsilon\mathfrak{A}\epsilon^{-1}=\mathfrak{A}$. ∎

Finally, we consider the symmetry of the stress tensor. Since the lower-dimensional manifolds model objects with finite width, the limit argument used to prove symmetry of the stress tensor is only valid within manifolds (and not transversely). Consequently, symmetry of the stress tensor is imposed within each manifold, expressed as:

[TABLE]

We remark that for $i\in I^{0}\cup I^{1}$, this equation is trivial since either $\sigma_{i}$ does not exist or is a vector. For $i\in I^{2}$, this equation evaluates the asymmetry of the in-plane components.

Gathering (3.5), (3.6), and (3.8), we arrive at the strong form of the generalized system of equations:

[TABLE]

System (3.9) has a structure similar to (3.2) in that it is composed of constitutive law(s) complemented with a differential and algebraic constraint. This structure is common in models concerning linear elasticity with relaxed symmetry arnold2006differential ; awanou2013rectangular . We will show in the next section that the system indeed corresponds to a symmetric saddle-point problem.

4 Variational Formulation

With the goal of obtaining a mixed finite element discretization, this section presents the weak formulation of the continuous problem. In order to do this, we introduce several analytical tools. First, the relevant function spaces are defined as well as the notational conventions concerning inner products. Next, we derive the variational formulation of (3.9) and show that it corresponds to a symmetric saddle point problem.

4.1 Function Spaces

The function spaces relevant for this problem are constructed as products of familiar function spaces on the $d$-dimensional manifolds. In particular, we define

[TABLE]

where $\Sigma$ denotes the function space for the stress, $\bm{U}$ contains the displacement, and $R$ is the function space for the Lagrange multiplier enforcing symmetry of the stress tensor. The exponent $k_{d}$ is given by $k_{d}=\left(\begin{smallmatrix}d\\ 2\end{smallmatrix}\right)=d(d-1)/2$ , see e.g. arnold2006differential ; awanou2013rectangular .

The mixed-dimensional $L^{2}$-inner products on $\Omega$ and $\Gamma$ are defined as the sum of inner products over all corresponding manifolds:

[TABLE]

Here, the implicit assumption is made that the contribution is zero for all manifolds on which $f$ is undefined. For example, for $\sigma,\tau\in\Sigma$, the inner product $(\sigma,\tau)_{\Omega}$ has no contribution on $\Omega^{0}$. Likewise for functions in $R$, the inner product is zero on manifolds $\Omega_{i}$ with $i\in I^{0}\cup I^{1}$.

For functions $\sigma,\tau\in\Sigma$, we note that they are defined on both $\Omega$ and $\Gamma$. For convenience, we introduce the combined inner product

[TABLE]

which, in the case of the operator $\mathfrak{A}$, is understood as

[TABLE]

The inner products naturally induce the $L^{2}$-type norms $\|\cdot\|_{\Omega}$, $\|\cdot\|_{\Gamma}$, and $\|\cdot\|_{\Omega\times\Gamma}$. With these norms, we assume that $\mathfrak{A}$ is continuous and coercive with respect to the norm $\|\cdot\|_{\Omega\times\Gamma}$. Thus, for all $\sigma,\tau\in\Sigma$, we have:

[TABLE]

We emphasize that the constants within these bounds are independent of $\epsilon$.

4.2 Identifying the Symmetric Saddle Point Problem

In this section, we make two key observations which allow us to derive a variational formulation of (3.9) which is symmetric. First let us consider the terms containing $u$ in the stress-strain relationships (3.9a) and (3.9a). We multiply these terms with $\tau\in\Sigma$ and $\bm{n}\cdot\tau\in L^{2}(\Gamma)$, respectively, and integrate to obtain the following integration by parts formula:

[TABLE]

Here $\mathfrak{D}\cdot$ is the mixed-dimensional divergence operator from (3.5).

Secondly, we introduce the operator ${\operatorname{skw}}$ which evaluates the asymmetric part of a matrix. More specifically, for a matrix $B\in\mathbb{R}^{n\times d}$ with components $b_{ij}$, let

[TABLE]

This operator is naturally lifted to ${\operatorname{skw}}:\Sigma\to R$. Next, we turn our attention to the term in (3.9a) containing the asymmetric variable $\chi$. Let us multiply this term with $\tau\in\Sigma$ and integrate over $\Omega^{2}\cup\Omega^{3}$. With the introduction of $r=\frac{1}{2}{\operatorname{skw}}\chi\in R$, we obtain

[TABLE]

By employing test functions $(\tau,\bm{v},s)\in\Sigma\times\bm{U}\times R$, the integration by parts formula (4.3), and the operator ${\operatorname{skw}}$, we obtain the following variational formulation of the problem (3.9):

Find $(\sigma,\bm{u},r)\in\Sigma\times\bm{U}\times R$ such that

[TABLE]

We identify system (4.4) as a saddle point problem by introducing the bilinear forms $a:\Sigma\times\Sigma\to\mathbb{R}$ and $b:\Sigma\times(\bm{U}\times R)\to\mathbb{R}$:

[TABLE]

The problem (4.4) can then be rewritten to the following, equivalent formulation:

Find $(\sigma,\bm{u},r)\in\Sigma\times\bm{U}\times R$ such that

[TABLE]

for all $(\tau,\bm{v},s)\in\Sigma\times\bm{U}\times R$.

5 Well-Posedness

In this section, we show well-posedness of the continuous formulation (4.4). The key is to associate appropriately weighted norms to the function spaces introduced in the previous section. In the mixed-dimensional setting considered here, let us endow $\Sigma$, $\bm{U}$, and $R$ with the following norms

[TABLE]

The proof of well-posedness consists of proving sufficient conditions on the bilinear forms $a$ and $b$ from (4.5) to invoke standard saddle-point theory. First, we show continuity of the operators, followed by ellipticity of $a$ and inf-sup on $b$.

Lemma 1 (Continuity)

The bilinear forms $a$ and $b$ from (4.5) are continuous with respect to the norms given by (5.1).

Proof

The continuity of $a$ follows from (4.2). For the blinear form $b$, we derive

[TABLE]

∎

Next, we focus on the bilinear form $a$. For the purposes of our analysis, it suffices to show that $a$ is elliptic on a specific subspace of $\Sigma$ generated by $b$. This is formally considered in the following lemma.

Theorem 5.1 (Ellipticity)

Given the bilinear forms $a$ and $b$ from (4.5). If $\sigma\in\Sigma$ satisfies

[TABLE]

then the following ellipticity bound holds

[TABLE]

Proof

We set $s=0$ in condition (5.2). The assumption holds for all $\bm{v}\in\bm{U}$, thus noting that $\mathfrak{D}\cdot\epsilon\sigma\in\prod_{i\in I}L^{2}(\Omega_{i})=\bm{U}$ and $\hat{\epsilon}_{\max}>0$, we obtain

[TABLE]

The proof is concluded by combining (5.3) with the coercivity of $\mathfrak{A}$ from (4.2). ∎

With the properties of $a$ proven, we continue by considering an inf-sup condition on the bilinear form $b$. This is shown in the following theorem, which relies on the constructions from Lemmas 2 and 4, presented afterwards.

Theorem 5.2 (Inf-Sup)

The bilinear form $b$ satisfies for all, $(\bm{u},r)\in\bm{U}\times R$,

[TABLE]

Proof

The proof consists of constructing a suitable $\tau\in\Sigma$ for a given pair $(\bm{u},r)\in\bm{U}\times R$. Its construction is based on constructing two auxiliary functions $\eta,\xi\in\Sigma$ using the techniques from Lemmas 2 and 4. Setting $\tau$ as the sum of these two functions then yields the result.

First, Lemma 2 allows us to construct $\eta\in\Sigma$ such that

[TABLE]

Secondly, we choose $\xi\in\Sigma$ using Lemma 4 with given $(r-\epsilon^{-1}{\operatorname{skw}}\eta)\in R$ such that

[TABLE]

By setting $\tau=\eta+\xi$, it follows that

[TABLE]

Furthermore, the bound on $\tau$ is derived using (5.4) and (5.5c)

[TABLE]

The proof is concluded by combining (5.6) and (5.7). ∎

Lemma 2

For each $\bm{u}\in\bm{U}$, a function $\eta\in\Sigma$ exists such that

[TABLE]

Proof

Considering $\bm{u}\in\bm{U}$ given, the function $\eta$ is constructed hierarchically. For each dimension $d$, we first set an interface value $\bm{\phi}$ on $\Gamma^{d}$, followed by a suitable extension into $\Omega^{d}$.

Given $i\in I^{0}$, we construct the adjacent interface functions $\bm{\phi}_{j}\in L^{2}(\Gamma_{j})$ such $\bm{\phi}_{j}=-\hat{\epsilon}_{\max}\bm{u}_{i}$ for a chosen $j\in\hat{J}_{i}$ with $(\hat{\epsilon})|_{\Gamma_{j}}=(\hat{\epsilon}_{\max})|_{\Omega_{i}}$ and zero for all other $j\in\hat{J}_{i}$. Repeating this construction for all $i\in I^{0}$, it follows that

[TABLE] 2. 1.

We continue with $i\in I^{1}$ and perform the following two steps. First, the function $\eta_{i}$ is constructed as the bounded $H(\operatorname{div},\Omega_{i})$-extension of the given $\bm{\phi}_{j}$ with $j\in\check{J}_{i}$. We use the extension operator as described in quarteroni1999domain (Section 4.1.2), giving us the properties

[TABLE]

Secondly, we further define $\bm{\phi}$ onto $\Gamma_{j}$ with $j\in\hat{J}_{i}$. We choose a single $j\in\hat{J}_{i}$ where $(\hat{\epsilon})|_{\Gamma_{j}}=(\hat{\epsilon}_{\max})|_{\Omega_{i}}$ and set $\bm{\phi}_{j}=-\hat{\epsilon}_{\max}\bm{u}_{i}+\hat{\epsilon}_{\max}^{-1}\nabla\cdot\epsilon\eta_{i}$. For all other $j\in\hat{J}_{i}$, we set $\bm{\phi}_{j}=0$. It then immediately follows that

[TABLE]

Repeating these two steps for all $i\in I^{1}$ gives us the bound

[TABLE]

in which the second and third inequalities follow from (3.4) and (5.10c). 3. 2.

For $n=3$, repeat the previous step for all $i\in I^{2}$ to obtain $\eta_{i}$ and $\bm{\phi}_{j}$ with $j\in\hat{J}_{i}$. 4. 3.

The construction of $\eta$ is finalized with its top-dimensional components $\eta_{i}$ with $i\in I^{n}$. Let the pair $(\eta_{i},\tilde{\bm{v}}_{i})\in(H({\operatorname{div}},\Omega_{i}))^{n}\times(L^{2}(\Omega_{i}))^{n}$ be the weak solution to the Poisson problem:

[TABLE]

This problem is solved for all $i\in I^{n}$. We then recall that $\epsilon=\hat{\epsilon}_{\max}=1$ in $\Omega^{n}$ and exploit the elliptic regularity of (5.13) (see e.g. evans1998partial ) to obtain

[TABLE]

With $\eta$ constructed, we consider its two main properties. First, by (5.9a), (5.11), and (5.13b), we have

[TABLE]

Secondly, we find the following bound from (5.9b), (5.10c), (5.12), and (5.14):

[TABLE]

thereby concluding the proof. ∎

Before introducing the second ingredient used in the proof of Theorem 5.2, we require several key analytical tools, organized in the following diagram:

[TABLE]

The function spaces ($\Theta$ and $W$) and mappings ($\Xi$, $\widetilde{\mathfrak{D}}\cdot$, and $\mathfrak{D}\times$) are defined next. Let the auxiliary space $\Theta$ be given by

[TABLE]

We emphasize that for an element $\theta\in\Theta$, this definition implies that $\theta_{i}$ is a 2-vector for $i\in I^{2}$ and a $3\times 3$ tensor for $i\in I^{3}$. Next, we follow awanou2013rectangular by introducing the mapping $\Xi$ and its right-inverse $\Xi^{-1}$ as:

[TABLE]

with $w_{i,\|}$ the tangential components of $w_{i}$ with respect to $\Omega_{i}$. $W$ is defined as the space of functions that lie in the image of the inverse operator and have a mixed-dimensional curl in $\Sigma$, i.e.

[TABLE]

We remark that for $w\in W$, $w_{i}$ is a 3-vector for $i\in I^{3}$ and a $3\times 3$ tensor for $i\in I^{3}$.

Next, we introduce a divergence-like operator $\widetilde{\mathfrak{D}}\cdot:\Theta\to R$ given by

[TABLE]

Here, $\check{\bm{n}}_{i}$ is the unique unit vector normal to $\Omega_{i}$ that forms a positive orientation with the chosen basis of the tangential bundle. By definition, this divergence operator maps from $\Theta$ to $R$, and we emphasize that $\widetilde{\mathfrak{D}}\cdot\theta$ is a vector for $d=3$ and a scalar for $d=2$.

Finally, the mixed-dimensional curl $(\mathfrak{D}\times)$ of $w\in W$ (see e.g. Licht ; boon2017excalc ) is given by

[TABLE]

Here, the superscript $\perp$ implies $[v_{1},v_{2}]^{\perp}=[-v_{2},v_{1}]$, e.g. $\nabla^{\perp}$ is a rotated gradient operator. We note that all differential operations are performed row-wise. Hence, for $n=3$, the mixed-dimensional curl maps to a $3\times 3$ tensor in $\Omega^{3}$, a $3\times 2$ tensor in $\Omega^{2}$ (in local coordinates) and a 3-vector in $\Omega^{1}$ (in local coordinates). Thus, an exact correspondence with the function space $\Sigma$ is obtained, as reflected in the diagram.

Lemma 3

The operators in diagram (5.19) enjoy the following two properties for all $w\in W$:

[TABLE]

Proof

The top row of (5.19) uses the differential operators from the mixed-dimensional De Rham complex boon2017excalc . The first equality then follows from the fact that exact forms are closed. It remains to show commutativity. By the definition of $\Xi$, see e.g. awanou2013rectangular ; Boffi , we have

[TABLE]

in which the subscript $d$ denotes a restriction of the operator to $\Omega^{d}$. Furthermore, we note that on $\Omega_{i}$ with $i\in I^{2}$, the skw operator evaluates the asymmetry with respect to the tangent bundle of $\Omega_{i}$. In turn, we have for $M\in\mathbb{R}^{3\times 3}$ that $({\operatorname{skw}}_{2}M)|_{\Omega_{i}}=\check{\bm{n}}_{i}\cdot{\operatorname{skw}}_{3}M$. This gives us

[TABLE]

∎

Lemma 4

Given $r\in R$, a function $\xi\in\Sigma$ exists such that

[TABLE]

Proof

We give the proof for $n=3$, the case $n=2$ being simpler. The strategy is to exploit the properties shown in Lemma 3 and first construct a bounded $\theta\in\Theta$ such that $\widetilde{\mathfrak{D}}\cdot\theta=\epsilon^{2}r$. Then, by setting $w=\Xi^{-1}\theta$ and $\xi=\epsilon^{-1}\mathfrak{D}\times w$, we obtain two of the desired properties

[TABLE]

The estimate will then follow from the boundedness of $\theta$.

The construction of $\theta$ proceeds according to the following three steps, consisting of an interface function $\phi\in H^{1}(\Gamma^{2})$ which serves as a source function for $\theta_{i}$ with $i\in I^{2}$ and a boundary condition for $\theta_{i}$ with $i\in I^{3}$.

We start by defining a scalar function $\phi$ in the trace space $H^{1}(\Gamma^{2})$. We let $\phi$ vanish at all intersections and extremities, i.e. $\phi$ is in the function space $\Phi$ given by

[TABLE]

Now, let $\phi$ be the solution to the following minimization problem:

[TABLE]

with $\Pi_{\mathbb{R}_{i}}$ the projection onto constants on $\Omega_{i}$. Due to the regularity of this problem and the imposed constraint, we have

[TABLE] 2. 2.

For each $i\in I^{2}$, we construct a function $\theta_{i}$ using $\phi$ as a source function. Specifically, let $(\theta_{i},p_{i})\in(H_{0}^{1}(\Omega_{i}))^{2}\times L^{2}(\Omega_{i})$ be the weak solution to the Stokes problem:

[TABLE]

The following bound is then satisfied from the regularity of (5.30a) (see e.g. evans1998partial ) combined with (5.29b)

[TABLE] 3. 3.

To finalize $\theta\in\Theta$, we create $\theta_{i}$ for $i\in I^{3}$ using $\phi$ from the first step as a boundary condition. Let $\theta_{i}\in(H^{1}(\Omega_{i}))^{3\times 3}$ and an auxiliary pressure variable $p_{i}\in(L^{2}(\Omega_{i}))^{3}$ be the weak solution to the following Stokes problem:

[TABLE]

Recall that $\check{\bm{n}}_{i}$, the unique normal vector of $\Omega_{i}$, and $\bm{n}_{j}$, the normal vector defined on $\Gamma_{j}$, are equal up to sign. This problem is well-posed since $\partial\Omega_{i}\setminus\Gamma^{2}$ has positive measure, for each $i\in I^{3}$, by assumption. We have the following bound due to the regularity of the Stokes problems and the fact that $\epsilon=1$ in $\Omega^{3}$

[TABLE]

Combining all $\theta_{i}$ from the final two steps gives us $\theta\in\Theta$. We first note that the properties (5.29a), (5.30b), and (5.32b) result in

[TABLE]

Hence, we have $\widetilde{\mathfrak{D}}\cdot\theta=\epsilon^{2}r$ and it follows from (5.26) that setting $w=\Xi^{-1}\theta$ and $\xi=\epsilon^{-1}\mathfrak{D}\times w$ provides the first two properties. The bound follows due to (5.29b), (5.31) and (5.33)

[TABLE]

as desired. ∎

With the proven properties of the bilinear forms $a$ and $b$, the main result of this section is summarized by the following theorem:

Theorem 5.3

Problem (4.4) is well-posed with respect to the norms (5.1). That is, a unique solution exists satisfying the bound

[TABLE]

Proof

It suffices to show continuity of the right-hand side of (4.4) with respect to the norms above. For that purpose, we derive the following bound on the first term using Cauchy-Schwarz and a trace inequality:

[TABLE]

Moreover, from (3.4), the second term is bounded as follows

[TABLE]

The result now follows readily from these two estimates, Theorems 5.1 and 5.2, and standard saddle point theory Boffi . ∎

6 Discretization

In this section, we discretize the continuous problem (4.4) using conforming finite elements. The notion of conformity and the choice of mixed finite element spaces is discussed in Section 6.1 and the resulting discretized problem is presented and analyzed in Section 6.2.

6.1 Discrete Spaces

For each $i\in I$, we introduce a shape-regular, simplicial grid $\Omega_{i,h}$ which tessellates $\Omega_{i}$. The union of meshes of a given dimension $d$ (with $0\leq d\leq n$) is denoted by $\Omega_{h}^{d}=\cup_{i\in I^{d}}\Omega_{i,h}$ and we define $\Omega_{h}=\cup_{i\in I}\Omega_{i,h}$. We let the grid respect all lower-dimensional features and be matching across all interfaces. The tesselation of $\Gamma$ is thus given by $\Gamma_{h}=\Gamma\cap\partial\Omega_{h}$. Moreover, there is an equivalence between $\Gamma_{j,h}$ and $\Omega_{i,h}$ for all $j\in\hat{J}_{i}$. The typical mesh size is denoted by $h$ and we use $h$ as a subscript to indicate the discretized counterpart of functions and function spaces.

With the aim of obtaining a stable and conforming method, we search for a discrete solution in subspaces of the function spaces defined in Section 4. We choose discrete function spaces $(\Sigma_{h},\bm{U}_{h},R_{h})$ on the grid $\Omega_{h}$ according to the following three conditions:

1. (S1)

The finite element spaces are conforming, i.e.

[TABLE] 2. (S2)

$\Sigma_{h}$ and $\bm{U}_{h}$ are such that for each $i\in I$:

[TABLE] 3. (S3)

A mixed-dimensional, finite element space $W_{h}$ exists such that

1. (a)

$\mathfrak{D}\times W_{h}\subseteq\Sigma_{h}$. 2. (b)

$(\Xi W_{i,h})\times R_{i,h}$ forms a stable pair for the two-dimensional Stokes problem for each $i\in I^{2}$. 4. (S4)

$\Sigma_{i,h}\times\bm{U}_{i,h}\times R_{i,h}$ forms a stable triplet for three-dimensional, mixed elasticity for each $i\in I^{3}$.

We provide an exemplary family of finite elements satisfying all four conditions. This choice is most concisely described using the notation of finite element exterior calculus AFW_FEEC . A translation to more conventional nomenclature is provided afterwards, for convenience. Given a polynomial degree $k\geq 0$, let

[TABLE]

In other words, for $n=3$: $\Sigma_{i,h}$ with $i\in I^{3}$ corresponds to three rows of Nedelec elements of the second kind ($N2^{f}_{k+1}$ Nedelec ) with degrees of freedom on the faces. For $i\in I^{2}$, it is three rows of Brezzi-Douglas-Marini elements ($BDM_{k+n-1}$ brezzi1985two ). Finally $\Sigma_{i,h}$ with $i\in I^{1}$ is given by a triplet of continuous Lagrange elements ($P_{k+2}$). The spaces $\bm{U}_{i,h}$ and $R_{i,h}$ are defined for $i\in I$ as three and $k_{d_{i}}$ rows, respectively, of discontinuous Lagrange elements ($P_{-(k+n-d_{i})}$).

In this case, the auxiliary space $W_{h}$ of (S3) is explicitly given by

[TABLE]

For $n=3$, $W_{i,h}$ is thus given by three rows of (first kind) edge-based Nedelec element ($N1_{k+2}^{e}$) for $i\in I^{3}$ and by three instances of Lagrange elements ($P_{k+n}$) for $i\in I^{2}$. The lowest order choice in this family, i.e. with $k=0$, is presented in Table 6.1.

A reduced family of finite elements arises by noting that all stability conditions remain valid after the polynomial order of the trace onto $\Gamma$ is reduced by one. Table 6.1 presents the lowest order member of this family.

[TABLE]

[TABLE]

Properties (S1)–(S3) can be verified with the use of the presented tables. Finally, these families correspond to the stable triplets analyzed in AFW_FEEC (Sections 11.6–11.7), hence (S4) holds as well.

6.2 Discrete Problem

Since the finite elements described above are contained in the continuous spaces from Section 5 by (S1), the discrete formulation of the model problem is a direct translation of (4.4):

Find $(\sigma_{h},\bm{u}_{h},r_{h})\in\Sigma_{h}\times\bm{U}_{h}\times R_{h}$ such that

[TABLE]

for all $(\tau_{h},\bm{v}_{h},s_{h})\in\Sigma_{h}\times\bm{U}_{h}\times R_{h}$. We note that the saddle-point structure of this problem has not changed, and can readily be uncovered using the bilinear forms from (4.5).

6.3 Stability

We continue with the analysis concerning the well-posedness of (6.20). Let us recall the norms from (5.1) for convenience

[TABLE]

Theorem 6.1 (Ellipticity)

If the discrete spaces satisfy (S2) and $\sigma_{h}\in\Sigma_{h}$ satisfies

[TABLE]

then the following ellipticity bound holds

[TABLE]

Proof

By (S2), we have $\mathfrak{D}\cdot\Sigma_{h}\subseteq\bm{U}_{h}$ and the same arguments are used as in Theorem 5.1. ∎

The next step is to consider the inf-sup condition for the bilinear form $b$ in the discrete case.

Theorem 6.2 (Inf-Sup)

If the discrete spaces satisfy conditions (S1)-(S4), then for all $(\bm{u}_{h},r_{h})\in\bm{U}_{h}\times R_{h}$,

[TABLE]

Proof

With $\bm{u}_{h}\in\bm{U}_{h}$ and $r_{h}\in R_{h}$ given, we follow a similar strategy as in the proof of Theorem 5.2. Here, we rely on Lemmas 5 and 6 to provide $\eta_{h}\in\Sigma_{h}$ to gain control of $\bm{u}_{h}$ and a divergence-free function $\xi_{h}\in\Sigma_{h}$ controlling $r_{h}\in R_{h}$.

In short, we choose $\eta_{h},\xi_{h}\in\Sigma_{h}$ such that

[TABLE]

and satisfy the bound

[TABLE]

Following the same steps as in Theorem 5.2, we define $\tau_{h}=\eta_{h}+\xi_{h}$ so that

[TABLE]

The proof is concluded by combining the above. ∎

Lemma 5

For each $\bm{u}_{h}\in\bm{U}_{h}$, a function $\eta_{h}\in\Sigma_{h}$ exists such that

[TABLE]

Proof

We use the same steps as in Lemma 2 to hierarchically construct $\eta_{h}\in\Sigma_{h}$. A concise exposition follows, starting with $d=0$. For each $i\in I^{d}$, we use (S2) to first construct a discrete $\bm{\phi}_{j,h}$ in the trace space $(\bm{n}\cdot\Sigma_{h})|_{\Gamma_{j}}$ for all $j\in\hat{J}_{i}$ such that

[TABLE]

in which $\nabla\cdot\epsilon\eta_{i,h}$ and $\bm{\phi}_{j,h}$ are understood as zero for $i\in I^{0}$, $j\in\check{J}_{i}$. For $i\in I^{d+1}$, the function $\eta_{i,h}$ is then defined as the bounded, discrete $H({\operatorname{div}},\Omega_{i,h})$-extension from quarteroni1999domain (Section 4.1.2). These steps are repeated by incrementing $d$ until $\bm{\phi}_{h}$ is completely defined on $\Gamma_{h}$. Finally, we solve a discrete Poisson problem for each $i\in I^{3}$, in analogy with (5.13), to complete $\eta_{h}\in\Sigma_{h}$. It follows by the same arguments as in Theorem 5.2 that the constructed $\eta_{h}$ satisfies (6.22). ∎

Lemma 6

Given $r_{h}\in R_{h}$, a function $\xi_{h}\in\Sigma_{h}$ exists such that

[TABLE]

with $\Pi_{R_{h}}$ the $L^{2}$ projection onto $R_{h}$.

Proof

In this proof, we make extensive use of the discrete space $W_{h}$ from stability requirement (S3). In particular, we will first introduce $w_{h}\in W_{h}$ such that $\mathfrak{D}\times w_{h}$ controls $r_{i,h}$ for $i\in I^{2}$. Then, a correction is introduced using (S4) in order to control $r_{i,h}$ for $i\in I^{3}$ as well. For brevity, we omit the subscript $h$ on all variables within this proof.

As in Lemma 4, we start with an interface function $\phi$ defined on the trace mesh $\Gamma_{h}^{2}$ which serves first as a source function and second as a boundary condition. We proceed according to the following four steps.

We consider functions on $\Gamma_{h}^{2}$ in the trace space of $(W_{h})|_{\Omega^{3}}$ that vanish at all intersections and extremities. Let us therefore introduce the function space $\Phi_{h}$ as

[TABLE]

It is important to note that the two instances of $\bm{n}$ in this definition are different. In particular, the former is defined as normal to $\partial\Omega_{i}$, of which $\Gamma_{j}$ is a subset, whereas the latter is normal with respect to the boundary $\partial\Gamma_{j}$.

We note that, since $\mathfrak{D}\times W_{h}\subseteq\Sigma_{h}$ by (S3), we have $W_{i,h}\subseteq(H({\operatorname{curl}},\Omega_{i}))^{3}$ for $i\in I^{3}$. In turn, it follows in the discrete setting that $\Phi_{j,h}\in(H({\operatorname{div}},\Gamma_{j}))^{3}$ for all $j\in\check{J}_{i}$, with the divergence tangential to $\Gamma_{j}$. Using this observation, we let $\phi$ solve the following minimization problem

[TABLE]

In other words, a finite-dimensional problem is solved for each $i\in I^{2}$ to obtain a bounded distribution that has an average asymmetry corresponding to $\epsilon^{2}r_{i}$. In particular, we obtain the following two properties

[TABLE] 2. 2.

We generate $w_{i}$ for $i\in I^{2}$ using the distribution $\phi$ from the first step as a source term. Introducing $\Theta_{i,h}=\Xi W_{i,h}$, it follows from (S3) that $\Theta_{i,h}\times R_{i,h}$ is a stable pair for the discretization of the Stokes problem (5.30a). Hence, we can find $\theta_{i}\in\Theta_{i,h}$ such that

[TABLE]

with $\Pi_{R_{i,h}}$ the $L^{2}$ projection onto $R_{i,h}$. We set $w_{i}=\Xi^{-1}\theta_{i}\in W_{i,h}$ for all $i\in I^{2}$. 3. 3.

For $i\in I^{3}$, let $w_{i}$ be given by any bounded extension of $\phi$ into $H({\operatorname{curl}},\Omega_{i})$, i.e. $w_{i}$ is chosen for all $i\in I^{3}$ such that

[TABLE] 4. 4.

Finally, we gain control of $r_{i}$ for $i\in I^{3}$. We recall that by stability condition (S4), the spaces $\Sigma_{i,h}\times\bm{U}_{i,h}\times R_{i,h}$ form a stable triplet for the mixed formulation of elasticity with relaxed symmetry. Considering $\Gamma^{2}$ as a zero traction boundary condition, we use the inf-sup condition associated to this stability to form a function $\psi_{i}\in\Sigma_{i,h}$ such that

[TABLE]

To complete the mixed-dimensional function $\psi\in\Sigma_{h}$, we set $\psi_{i}=0$ for $i\not\in I^{3}$.

Using the above ingredients, we set $\xi=\psi+\epsilon^{-1}\mathfrak{D}\times w$ and obtain the first two desired properties:

[TABLE]

The bound now follows by the estimates given in each step:

[TABLE]

∎

The two previous theorems provide the sufficient ingredients to show stability of the discretization, formally presented in the following theorem.

Theorem 6.3 (Stability)

If the discrete spaces satisfy the conditions (S1)-(S4), then the resulting mixed finite element method is stable, i.e. a unique solution exists satisfying the bound

[TABLE]

with a constant independent of the grid size $h$.

Proof

Using Theorems 6.1 and 6.2, this result follows from standard saddle point theory Boffi . ∎

6.4 Convergence

By consistency of the discretized problem (6.20) with respect to the continuous formulation (4.4) and stability from Theorem 6.3, we have shown that the proposed mixed finite element discretization is convergent. In turn, this section is devoted to obtaining the rates of convergence through a priori error estimation.

Let $\Pi_{R_{h}}$ and $\Pi_{U_{h}}$ be the $L^{2}$-projection operators onto the finite element spaces $R_{h}$ and $\bm{U}_{h}$. Under the assumption of sufficient regularity, we introduce the canonical projection operator $\Pi_{\Sigma_{h}}$ such that the commutativity properties hold for $\sigma\in\Sigma$:

[TABLE]

A direct consequence of these two properties is that for sufficiently regular $\sigma\in\Sigma$, we have the commuting property

[TABLE]

Using $\|\cdot\|_{\rho,\Omega}$ as short-hand notation for the $H^{\rho}(\Omega)$-norm, the projection operator $\Sigma_{h}$ has the following approximation properties for $0\leq\rho\leq k_{i}+1$:

[TABLE]

Note that for $\Omega^{0}$ and $\Gamma^{0}$, the projection operators correspond to the identity operator which makes (6.35c) and (6.35d) trivial there.

Theorem 6.4 (Convergence)

Given $(\sigma,\bm{u},r)$ as the solution to (4.4). Under the assumption of sufficient regularity, then the scheme converges optimally, i.e. the finite element solution $(\sigma_{h},\bm{u}_{h},r_{h})$ to (6.20) satisfies the following estimate

[TABLE]

Proof

We restrict our choice of test functions $(\tau_{h},\bm{v}_{h},s_{h})$ to the finite element spaces and subtract the discrete equations (6.20) from the continuous equations (4.4). We then obtain

[TABLE]

The projection operators from (6.35) are then used to project the true solution onto the finite element spaces. Introducing $\sigma_{\Pi}:=\Pi_{\Sigma_{h}}\sigma$, $\bm{u}_{\Pi}:=\Pi_{U_{h}}\bm{u}$, and $r_{\Pi}:=\Pi_{R_{h}}r$, we rewrite the above equation to

[TABLE]

Let $\tau_{b,h}$ be the discrete stress from the construction in Theorem 6.2 with the following properties

[TABLE]

The discrete test functions are then chosen to be

[TABLE]

with $\delta>0$ a constant to be determined later. Substituting this choice of functions into the left-hand side of (6.4) gives us

[TABLE]

Next, due to (6.34) and (S2), we have $\mathfrak{D}\cdot\epsilon(\sigma_{\Pi}-\sigma_{h})=\Pi_{U_{h}}\mathfrak{D}\cdot\epsilon(\sigma-\sigma_{h})=0$. Hence, the coercivity of $\mathfrak{A}$ from (4.2) allows us to bound $a(\sigma_{\Pi}-\sigma_{h};\sigma_{\Pi}-\sigma_{h})$ from below by $\|\sigma_{\Pi}-\sigma_{h}\|_{\Sigma}^{2}$. We then obtain the following bound with respect to the right-hand side of (6.4):

[TABLE]

Next, we use the continuity of the forms $a$ and $b$ from Lemma 1 to bound the right-hand side further

[TABLE]

An application of Young’s inequality and rearranging terms then gives

[TABLE]

Using (6.37b) and setting $\delta$ sufficiently small leads us to

[TABLE]

With the triangle inequality, we thus obtain the estimate

[TABLE]

An application of the approximation properties (6.35) then finishes the proof. ∎