A mass conserving mixed stress formulation for Stokes flow with weakly imposed stress symmetry
Jay Gopalakrishnan, Philip L. Lederer, Joachim Sch\"oberl
TL;DR
This paper presents a novel finite element method for Stokes flow that weakly enforces stress symmetry, ensuring mass conservation, stability, and optimal convergence rates for velocity, pressure, and stress variables.
Contribution
It introduces a new discretization that directly approximates symmetric viscous stresses with weak enforcement, improving stability and accuracy over previous methods.
Findings
Achieves optimal convergence rates for pressure and stress.
Ensures exact mass conservation with $H(div)$-conforming velocity.
Method is pressure robust and stable.
Abstract
We introduce a new discretization of a mixed formulation of the incompressible Stokes equations that includes symmetric viscous stresses. The method is built upon a mass conserving mixed formulation that we recently studied. The improvement in this work is a new method that directly approximates the viscous fluid stress , enforcing its symmetry weakly. The finite element space in which the stress is approximated consists of matrix-valued functions having continuous "normal-tangential" components across element interfaces. Stability is achieved by adding certain matrix bubbles that were introduced earlier in the literature on finite elements for linear elasticity. Like the earlier work, the new method here approximates the fluid velocity using -conforming finite elements, thus providing exact mass conservation. Our error analysis shows optimal…
| ( | eoc ) | ( | eoc ) | ( | eoc ) | ( | eoc ) | ( | eoc ) | |
| 20 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 80 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 320 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1280 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 5120 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 20 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 80 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 320 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1280 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 5120 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 20 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 80 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 320 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1280 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 5120 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| ( | eoc ) | ( | eoc ) | ( | eoc ) | ( | eoc ) | ( | eoc ) | |
| 20 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 80 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 320 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1280 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 5120 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 20 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 80 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 320 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1280 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 5120 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 20 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 80 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 320 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1280 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 5120 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| ( | eoc ) | ( | eoc ) | ( | eoc ) | ( | eoc ) | ( | eoc ) | |
| 28 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 224 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1792 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 14336 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 114688 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 28 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 224 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1792 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 14336 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 114688 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 28 | ( | – ) | ( | – ) | ( | – ) | ( | – ) | ( | – ) |
| 224 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 1792 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 14336 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
| 114688 | ( | ) | ( | ) | ( | ) | ( | ) | ( | ) |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
A mass conserving mixed stress formulation for Stokes flow with weakly imposed stress symmetry
Jay Gopalakrishnan
Portland State University, PO Box 751, Portland OR 97207,USA
,
Philip L. Lederer
Institute for Analysis and Scientific Computing, TU Wien, Wiedner Hauptstraße 8-10, 1040 Wien, Austria
and
Joachim Schöberl
Institute for Analysis and Scientific Computing, TU Wien, Wiedner Hauptstraße 8-10, 1040 Wien, Austria
Abstract.
We introduce a new discretization of a mixed formulation of the incompressible Stokes equations that includes symmetric viscous stresses. The method is built upon a mass conserving mixed formulation that we recently studied. The improvement in this work is a new method that directly approximates the viscous fluid stress , enforcing its symmetry weakly. The finite element space in which the stress is approximated consists of matrix-valued functions having continuous “normal-tangential” components across element interfaces. Stability is achieved by adding certain matrix bubbles that were introduced earlier in the literature on finite elements for linear elasticity. Like the earlier work, the new method here approximates the fluid velocity using -conforming finite elements, thus providing exact mass conservation. Our error analysis shows optimal convergence rates for the pressure and the stress variables. An additional post processing yields an optimally convergent velocity satisfying exact mass conservation. The method is also pressure robust.
Key words and phrases:
mixed finite element methods; incompressible flows; Stokes equations; weak symmetry
Philip L. Lederer has been funded by the Austrian Science Fund (FWF) through the research program “Taming complexity in partial differential systems” (F65) - project “Automated discretization in multiphysics” (P10).
1. Introduction
In this work we introduce a new method for the discretization of steady incompressible Stokes system that includes symmetric viscous stresses. Let be a bounded domain with or having a Lipschitz boundary . Let and be the velocity and the pressure, respectively. Given an external body force and kinematic viscosity , the velocity-pressure formulation of the Stokes system is given by
[TABLE]
where . By introducing a new variable where , equation (1) can be reformulated to
[TABLE]
We shall call formulation (2) the mass conserving mixed formulation with symmetric stresses, or simply the MCS formulation. Although formulations (1) and (2) are formally equivalent, the MCS formulation (2) demands less regularity of the velocity field . Many authors have studied this formulation previously [15, 14, 13, 12], including us [18]. In [18], following the others, we introduced a new variable , which is in general nonsymmetric, and considered an analogous formulation (which was also called an MCS formulation). The main novelty in [18] was that was set in a new function space of matrix-valued functions whose divergence can continuously act on elements of . Accordingly, the appropriate velocity space there was not as in the classical velocity-pressure formulation.
In contrast to [18], in this work we set , not . Our goal is to apply what we learnt in [18] to produce a new method that provides a direct approximation to the symmetric matrix function . Being the viscous stress, this is of more direct practical importance (than ). We shall seek in the same function space that we considered in [18]. We have shown in [18] that matrix-valued finite element functions with “normal-tangential” continuity across element interfaces are natural for approximationg solutions in We shall continue to use such finite elements here. It is interesting to note that in the HDG (hybrid discontinuous Galerkin) literature [11, 16] the potential importance of such normal-tangential continuity was noted and arrived at through a completely different approach.
The main point of departure in this work, stemming from that fact that contains non-symmetric matrix-valued functions, is that we impose the symmetry of stress approximations weakly using Lagrange multipliers. This technique of imposing symmetry weakly is widely used in finite elements for linear elasticity [1, 2, 3, 14]. In particular, our analysis is inspired by the early work of Stenberg [30], who enriched the stress space by curls of local element bubbles. (In fact, this idea was even used in a Stokes mixed method [15], but their resulting method is not pressure robust.) These enrichment curls lie in the kernel of the divergence operator and are only “seen” by the weak-symmetry constraint allowing them to be used to prove discrete inf-sup stability. While in two dimensions – assuming a triangulation into simplices – this technique only increases the local polynomial order by , this is not the case in three dimensions. Years later [8, 17], it was realized that it is possible to retain the good convergence properties of Stenberg’s construction and yet reduce the enrichment space. Introducing a “matrix bubble,” these works added just enough extra curls needed to prove stability.
We shall see in later sections that the matrix bubble can also be used to enrich our discrete fluid stress space. This might seem astonishing at first. Indeed, an enrichment space for fluid stresses must map well when using a specific map that is natural to ensure normal-tangential continuity of the discrete stress space. Moreover, the enrichment functions must lie in the kernel of a realization of the distributional row-wise divergence used in MCS formulations (displayed in (11) below). It turns out that these properties are all fulfilled by an enrichment using a double curl involving matrix bubbles. Hence we are able to prove the discrete inf-sup condition. Stability then follows in the same type of norms used in [30] and is a key result of this work.
Some comments on the choice of the discrete velocity space and its implications are also in order here. As mentioned above, the velocity space within the MCS formulation is . One of the main features of the first MCS method [18], as well the new version with weakly imposed symmetry of this paper, is that we can choose a discrete velocity space using -conforming finite elements. Therefore, our method is tailored to approximate the incompressibility constraint exactly, leading to pointwise and exactly divergence-free discrete velocity fields. The use of such -conforming velocities in Stokes flow is by no means new: for the standard velocity-pressure formulation, once can find it in [9, 10], and for the Brinkman Problem in [20]. Therein, and also in the more recent works of [25, 24], the -conformity is treated in a weak sense and a (hybrid) discontinuous Galerkin method is constructed. When employing -conforming finite elements, one has the luxury of choice. In [18], we used the space [6] and added several local stress bubbles in order to guarantee stability. In contrast, in this paper, we have chosen to take the smaller Raviart-Thomas space [26] of order , denoted by . A similar choice was made also in the work of [16], where they presented a hybrid method for solving the Brinkman problem based off the work of [11]. Our current choice of the smaller space leads to a less accurate velocity approximation (compared to ), so in order to retain the optimal convergence order of the velocity (measured in a discrete -norm), we introduce a local element-wise post processing. Using the reconstruction operator of [21, 22] this post processing can be done retaining the exact divergence-free property.
The remainder of this paper is organized as follows. In Section 2, we define notation for common spaces used throughout this work and introduce an undiscretized formulation. Section 3 presents the MCS method for Stokes flow including symmetric viscous stresses. In Section 4, we present the new discrete method including the introduction of the matrix bubble. Section 5 proves a discrete inf-sup condition and develops a complete a priori error analysis of the discrete MCS system. In Section 6, we introduce a postprocessing for the discrete velocity. The concluding section (Section 7) reports various numerical experiments we performed to illustrate the theory.
2. Preliminaries
In this section, we introduce notation and present a weak formulation for Stokes flow that includes symmetric viscous stresses.
Let or denote the set of infinitely differentiable compactly supported real-valued functions on and let denote the space of distributions. To differentiate between scalar, vector and matrix-valued functions on , we include the co-domain in the notation, e.g., . Let denote the vector space of real matrices. This notation scheme is similarly extended to other function spaces as needed. Thus, denotes the space of square integrable -valued functions on , while analogous vector and matrix-valued function spaces are defined by L^{2}(\Omega,\mathbb{R}^{d}):=\left\{u:\Omega\to\mathbb{R}^{d}\big{|}\;u_{i}\in L^{2}(\Omega)\right\} and L^{2}(\Omega,{\mathbb{M}}):=\left\{\sigma:\Omega\to{\mathbb{M}}\big{|}\;\sigma_{ij}\in L^{2}(\Omega)\right\}, respectively. Let denote the vector space of skew symmetric matrices, i.e., , and let L^{2}(\Omega,\mathbb{K}):=\left\{\sigma:\Omega\to\mathbb{K}\big{|}\;\sigma_{ij}\in L^{2}(\Omega)\right\}.
Recall that the dimension in this work is either 2 or 3. Accordingly, depending on the context, certain differential operators have different meanings. The “curl” operator, depending on the context, denotes one of the differential operators below.
[TABLE]
where denotes the transpose and abbreviates For matrix-valued functions in both and cases, i.e., by we mean the matrix obtained by taking row wise. Unfortunately, this still does not exhaust all the curl cases. In the case, there are two possible definitions of for ,
[TABLE]
and we shall have occasion to use both. The latter will not be used until (14) below, so until then, the reader may continue assuming we mean (3) whenever we consider curl of vector functions in . The operator is to be understood from context as an operator that results in either a vector whose components are for or a matrix whose entries are for or a third-order tensor whose entries are for Finally, in a similar manner, we understand as either for vector-valued or the row-wise divergence for matrix-valued .
Let (so that and for and 3, respectively). In addition to the standard Sobolev space for any we shall use the well-known space By its trace theorem, is a well-defined closed subspace, where denotes the outward unit normal on . Its dual space , as proved in [18, Theorem 2.1], satisfies
[TABLE]
In this work, the following space is important:
[TABLE]
where the name results from (5): indeed a function fulfills .
Next, let us derive a variational formulation of the system (2), which is based on the mixed stress formulation (MCS) introduced in chapter 3 in the work [18]. The method is based on a weaker regularity assumption of the velocity as compared to the standard velocity-pressure formulation (1). The velocity and the pressure now belong, respectively, to the spaces
[TABLE]
Multiplying (2c) with a pressure test function and integrating over the domain ends up in the familiar equation which we write as the last equation of the final Stokes system (7) written below. Here and throughout, the inner product of a space is denoted by . When is the space of functions whose components are square integrable functions on , we abbreviate to simply , as done in (7) below. Similarly, while we generally denote the norm and seminorm on a Sobolev space by and , respectively, to simplify notation, we set , where denotes inner product for any and any subset . Moreover, when , we omit the subscript and simply write for .
To motivate the remaining equations of (7), let the deviatoric part of a matrix be defined by where denotes the identity matrix and denotes the matrix trace. Since , due to the incompressibility constraint we have the identity
[TABLE]
Since and , we define the stress space as the following closed subspace of :
[TABLE]
Testing equations (2a) with a test functions and integrating over the domain, we have for the term including the identity
[TABLE]
Using the knowledge that the velocity should be in we obtain
[TABLE]
which is the first equation in the system (7) below. Here and throughout, when working with elements of the dual space of a topological space , we denote the action of on an element by , where we may omit the subscript when its obvious from context. Finally we also test (2b) with and integrate the pressure term by parts. This results in the remaining equation of (7).
Summarizing, the weak problem is to find such that
[TABLE]
In the ensuing section, we shall focus on a discrete analysis of a nonconforming scheme based on (7). Although wellposedness of (7) is an interesting question, we shall not comment further on it here since it is of no direct use in a nonconforming analysis.
3. The new method
In [18], we introduced an MCS method where was an approximation to (the generally non-symmetric) instead of (the symmetric) considered above. Since there was no symmetry requirement in [18], there we worked with the space instead of . The finite element space for designed there can be reutilized in the current symmetric case (with some modifications), once we reformulate the symmetry requirement as a constraint in a weak form.
To do so, we need further notation. Let be defined by
[TABLE]
When represents the Stokes velocity, represents the vorticity. Since , introducing as a new variable, and the symmetry condition as a new constraint, we obtain the boundary value problem
[TABLE]
In the remainder of this section, we introduce a discrete formulation approximating (9).
The method will be described on a subdivision (triangulation) of consisting of triangles in two dimensions and tetrahedra in three dimensions. For the analysis later, we shall assume that the is quasiuniform. By we denote the maximum of the diameters of all elements . Quasiuniformity implies that for all mesh elements . Here and throughout, by we indicate that there exist two constants independent of the mesh size as well as the viscosity such . Similarly, we use the notation if there exists a constant such that . All element interfaces and element boundaries on are called facets and are collected into a set . This set is partitioned into facets on the boundary and interior facets . On each facet we denote by the standard jump operator. On a boundary facet the jump operator is just the identity. On all facets we denote by a unit normal vector. When integrating over boundaries of -dimensional domains, the orientation of is assumed to be outward. On a facet with normal adjacent to an mesh element , the normal and tangential traces of a smooth function are defined by and respectively. Similarly, for a smooth , the (scalar-valued) “normal-normal” and the (vector-valued) “normal-tangential” components are defined by and respectively.
For any integers , the following “broken spaces” are viewed as consisting of functions on without any continuity constraints across element interfaces:
[TABLE]
For we use the notation for the inner product of or its vector and tensor analogues such as Also let . Next for each element let denote the set of polynomials of degree at most on . The vector and tensor analogues such as have their components in . The broken spaces and are defined similarly. We shall also use the conforming Raviart-Thomas space (see [4, 27]),
3.1. Velocity, pressure, and vorticity spaces
For any , our method uses
[TABLE]
for approximating the velocity, pressure, and vorticity, respectively.
Standard finite element mappings apply for these spaces. Let be the unit simplex (for and ), which we shall refer to as the reference element, and let . Let be an affine homeomorphism and set . By quasiuniformity, and estimates that we shall use tacitly in our scaling arguments later. Such arguments proceed by mapping functions on to and from . Given a scalar-valued , a vector-valued , and a skew-symmetric matrix-valued on the reference element , we map them to using
[TABLE]
respectively, i.e., these are our mappings for functions in the pressure, velocity, and vorticity spaces, respectively. The first is the inverse of the standard pullback, the second is the standard Piola map, and the third is designed to preserve skew symmetry.
3.2. Stress space
The definition of our stress space is motivated by the following result, proved in [18, Section 4].
Theorem 1**.**
Suppose is in and for all elements . Assume that the normal-tangential trace is continuous across element interfaces. Then is in and moreover
[TABLE]
for all
Clearly, matrix finite element subspaces having normal-tangential continuity are suggested by Theorem 1. Technically, the theorem’s sufficient conditions for full conformity also include the condition This condition is very restrictive as it would enforce continuity at vertices and edges in two and three dimensions respectively. If this constraint is relaxed, much simpler, albeit nonconforming, elements can be constructed. This was the approach we adopted in [18]. We continue in the same vein here and define the nonconforming stress space
[TABLE]
As mentioned in the introduction, we must enrich the above stress space to guarantee solvability of the resulting discrete system due to the additional weak symmetry constraints. We follow the approach of [30] and its later improvements [8, 17] to construct the needed enrichment space.
Define a cubic matrix-valued “bubble” function as follows. On a -simplex with vertices , let denote the face opposite to , and let denote the unique linear function that vanishes on and equals one on , i.e., the th barycentric coordinate of . Following [8, 17], we define by
[TABLE]
where the indices on the barycentric coordinates are calculated mod in (13a). Let denote the -orthogonal complement of in for , and let For any , define
[TABLE]
for and , with the understanding that in case, the outer curl is defined by (4), not (3). The total stress space is given by
[TABLE]
That functions in this space have normal-tangential continuity is a consequence of the following property proved in [8, Lemma 2.3].
Lemma 2**.**
Let and . The products and have vanishing tangential trace on , so the function has vanishing normal trace on .
Lemma 3**.**
Any has vanishing and on all facets
Proof.
Since , this is a direct consequence of Lemma 2. ∎
We also need a proper mapping for functions in that preserves normal-tangential continuity. We shall continue to use the following map, first introduced in [18]:
[TABLE]
As shown in [18, Lemma 5.3], on each facet, is a scalar multiple of and if and only if Degrees of freedom are discussed in §3.4.
Remark 4*.*
Note that in (13), was given using barycentric coordinates as an expression that holds on any simplex. Let denote the function on the reference element obtained by replacing by reference element barycentric coordinates . Considering the obvious map that transforms to , we find that the matrix bubble on any simplex is given by
[TABLE]
3.3. Equations of the method
For the derivation of the discrete variational formulation we turn our attention back to the weak formulation (7) and identify these forms:
[TABLE]
The definition of the remaining bilinear form is motivated by the definition of the “distributional divergence” given by (11). To this end we define by
[TABLE]
Integrating the first integral by parts, we find the equivalent representation
[TABLE]
Using these forms, we state the method. For any , the discrete MCS method with weakly imposed symmetry finds such that
[TABLE]
Since and fulfills , the discrete velocity solution component satisfies point wise, providing exact mass conservation.
3.4. Degrees of freedom of the new stress space
We need degrees of freedom (d.o.f.s) for the stress space that are well-suited for imposing normal-tangential continuity across element interfaces. Since the bubbles in have zero normal-tangential continuity, we ignore them for this discussion and focus on d.o.f.s that control .
Consider on any mesh element . Letting denote the subspace of matrices satisfying we may identify with . Let us recall a basis for that was given in [18]. Define the following two sets of constant matrix functions, for and cases, respectively, by
[TABLE]
taking the indices mod 3 and mod 4, respectively. We proved in [18, Lemma 5.1] that the sets and form a basis of when and , respectively.
Our d.o.fs for are grouped into two. The first group is associated to the set of element facets ( subsimplices of ), namely, for each facet , we define the set of d.o.f.s
[TABLE]
for each in any fixed basis for . The next group is the set of interior d.o.f.s, defined by
[TABLE]
for all in any basis of . We proceed to prove that the set of these d.o.f.s, , is unisolvent.
Theorem 5**.**
The set is a set of unisolvent d.o.f.s for .
Proof.
Suppose satisfies for all d.o.f.s . We need to show that . From the facet d.o.f.s we conclude that vanishes on . By [18, Lemma 5.2], may be expressed as
[TABLE]
when or , respectively, where . The interior d.o.f.s imply that for any . Choosing for the expression on the right hand side in (21) omitting the , say for the case, we obtain
[TABLE]
yielding , and thus . A similar argument in case yields the same conclusion that .
To complete the proof, it now suffices to prove that equals the number of d.o.f.s, i.e., . Obviously, . The cardinality of equals the sum of the number of facet d.o.f.s and the number of interior d.o.f.s , which simplifies to (d^{2}-1)\big{(}\dim{\mathbb{P}}^{k-1}(T)+\dim{\mathbb{P}}^{k}(F)\big{)}, equalling . ∎
Using these d.o.f.s, a canonical local interpolant in can be defined as usual, by requiring that for all
Lemma 6**.**
For any we have
Proof.
This proceeds along the same lines as the proof of [18, Lemma 5.4]. ∎
The global interpolant is also defined as usual. On each element the global interpolant coincides with the local interpolant .
Theorem 7**.**
For any and any , the global interpolation operator satisfies
[TABLE]
for all .
Proof.
This follows from a standard Bramble-Hilbert argument using Lemma 6. ∎
4. A priori error analysis
In this section we first show the stability of the MCS method with weakly imposed symmetry by proving a discrete inf-sup condition (Theorem 21). We then prove consistency (Theorem 25), optimal error estimates (Theorem 26), and pressure robustness (Theorem 28). For simplicity, the analysis from now on assumes that is a constant.
4.1. Norms
In addition to the previous notation for norms (established in Section 2), hereon we also use to abbreviate , a notation that also serves to indicate that certain seminorms are defined using differential operators applied element by element, not globally, e.g.,
[TABLE]
for and . Recall that . Our analysis is based on norms of the type used in [30]. Accordingly, we will need to use the following norms for and :
[TABLE]
Lemma 15 below will show that the latter is indeed a norm.
On the discrete space , we will also need another norm defined using the following projections. On any mesh element , let denote the orthogonal projection onto where is determined from context to be an appropriate vector space such as or . When the element is clear from context, we shall drop the subscript in and simply write . Also, on each facet we introduce a projection onto the tangent plane : for any , the projection is defined by for all . Using these, define
[TABLE]
Lemma 14 below will help us go between this norm and .
The remaining spaces and are simply normed by the norm . The full discrete space is normed by
[TABLE]
for any .
4.2. Norm equivalences
Next, we use the finite element mappings introduced earlier –see (10) and (15)– to show several norm equivalences.
Lemma 8**.**
Let . Then
[TABLE]
Proof.
The first two follow by a simple scaling argument. For the third, see the proof of [18, Lemma 6.1]. ∎
In the proof of the next lemma, we use the space of rigid displacements For each element , let denote the projector defined in [5]. Then, for any the projection fulfills the properties (see [5, eq. (3.3), (3.11)])
[TABLE]
We shall also use a global discrete Korn inequality, implied by [5, Theorem 3.1]. Namely, there is an -independent constant such that
[TABLE]
Lemma 9**.**
For all ,
[TABLE]
Proof.
One side of the equivalence is obvious by the continuity of the . For the other direction first note that As we have again by the continuity of
[TABLE]
We conclude the proof using (28). ∎
The following well-known property of Raviart-Thomas spaces (see, e.g., [7, Lemma 3.1]) is needed at several points.
Lemma 10**.**
*Let and . Then is in . *
Lemma 11**.**
For all and
[TABLE]
Proof.
One side of the equivalence of (30) is obvious by the continuity of the . For the other direction, we use the following equivalence on the reference element :
[TABLE]
This follows by finite dimensionality, because by the Euler identity if any one of the above two terms is zero, then (see e.g., [23]). Consequently, given any , setting , the following problem is uniquely solvable: find such that
[TABLE]
Since , (34) implies that . Put . Then, due to the properties of the Piola map , is a function in satisfying in , and a scaling argument using (33) implies
[TABLE]
Let . Then and in . In particular, the former implies, by Lemma 10, that . Then we have
[TABLE]
This proves (30).
To prove (31), first note that due to the definition of we have . Thus, using the same decomposition as above, namely, ,
[TABLE]
As , the first term on the right vanishes. The last term satisfies
[TABLE]
due to (35). Hence (31) is proved.
The proof of (32) uses the same technique:
[TABLE]
where we have used that and (35). ∎
Remark 12*.*
The same technique shows that for all Raviart-Thomas functions . The technique allows one to control the gradient of the highest order terms of a velocity in the Raviart-Thomas space by . A similar estimate does not hold for in
Lemma 13**.**
For all and ,
[TABLE]
Proof.
The proof is based on a scaling argument and equivalence of norms on finite dimensional spaces on the reference element. Recall the map and . Calculations using the chain rule yield
[TABLE]
We continue with the case only (since case proceeds using (36b) analogously). With , standard estimates for yield
[TABLE]
Let and be such that and , where is as defined in (8). Then,
[TABLE]
In view of (37) and (38), to complete the proof, it suffices to establish the reference element estimate
[TABLE]
by proving that one side is zero if and only if the other side is zero. Note these two identities: and . If , then the latter identity implies , which when used in the former identity, yields . Combined with the obvious converse, we have established (39). ∎
Lemma 14**.**
For all and ,
[TABLE]
Proof.
Since the decomposition is orthogonal in the Frobenius inner product, so is Application of the deviatoric and preserves this orthogonality. Hence, by Pythagoras theorem,
[TABLE]
We shall now prove the result using (40) and Lemma 11.
Proof of “”: Since
[TABLE]
it suffices to prove that
[TABLE]
which we do next. Since the projection can be bounded using (40), we focus on the remainder .
[TABLE]
When this estimate for is used in and is bounded using (40), we obtain (41).
Proof of “”: The last term of the lemma obviously satisfies , while the first term satisfies by (40). It remains to bound . As , we obtain using an inverse inequality for polynomials
[TABLE]
where we used (27) in the last step. ∎
Lemma 15**.**
For any and
[TABLE]
While the first estimate in (42) involves only the local constants from Lemmas 13 and 14, using the global constant we also have
[TABLE]
Proof.
To prove the first estimate of (42),
[TABLE]
Taking infimum over , we obtain the lower estimate of (42). The upper bound of the first infimum obviously follows by choosing .
To prove the equality in (42), observe that the infimum over cannot be larger than because we may choose . The reverse inequality also holds since for any , so the equality must hold.
Finally, to prove (43), we use triangle inequality to get
[TABLE]
Applying the Korn inequality (29) and noting that the jump of the normal components are zero for functions in , the proof is complete. ∎
4.3. Stability analysis
The next three lemmas lead us to a discrete inf-sup condition.
Lemma 16**.**
Let for some and . Then for ,
[TABLE]
Proof.
If , then obviously . We claim that the converse is also true. Indeed, if , then putting , we have
[TABLE]
Taking divergence on both sides, we find that , so must be a constant on . Then, taking normal components of both sides of (44) on each facet, we find that , so . Hence , which in turn implies that . Therefore, by [8, Lemma 2.2], .
Applying this on the reference element for and where is in Remark 4, by finite dimensionality, we have
[TABLE]
We will now show that is related to by
[TABLE]
By the definition of ,
[TABLE]
as trace is preserved under similarity transformations. Focusing on the part of the last term inside the deviatoric, in the case,
[TABLE]
This proves that
[TABLE]
when . The same identity holds in the case: the argument is similar after changing the definitions of the curls and the mapping of appropriately. Thus, and (46) is proved.
Finally, the result follows from (46) by scaling arguments: indeed (45) implies, by (24) and (36) that
[TABLE]
from which the result follows. ∎
Lemma 17**.**
For any , there is a such that
[TABLE]
Furthermore, for any the same pair satisfies
[TABLE]
Proof.
Given a , set element by element by
[TABLE]
Clearly, is in Since is in we conclude that . Since is trace-free, which in turn implies, after integrating by parts and applying Lemma 2, .
In the case, this yields
[TABLE]
Noting that , where is the distance from the th vertex to the facet of the simplex opposite to it, and that the -norm of any matrix is equivalent to the sum of -norms of , a local scaling argument with and (49) imply
[TABLE]
Therefore, , by Lemma 16. This proves (47) in the case. In the case, the analogue of (49) gives where we have used Lemma 16 again. This completes the proof of (47).
To prove (48), we use (18). The last sum in
[TABLE]
vanishes due to Lemma 3. Hence by (47),
[TABLE]
To handle the last term, note that
[TABLE]
because . This follows by integrating one of the curls by parts, observing that the resulting volume term is zero (since ) and so is the resulting boundary term (due to Lemma 2). Continuing, we apply Cauchy-Schwarz inequality and an inverse inequality to get
[TABLE]
by Lemma 16. Returning to (50) and using this estimate, the proof is complete. ∎
Remark 18*.*
The message of Lemmas 16 and 17 is that it is possible to choose a in the form of a deviatoric of a curl of a bubble to bound (from below) the term arising from the weak symmetry constraint. If was just a curl, it would not be seen by the equilibrium equation and the bound in (48) would not have the -term, but our is a deviatoric (of a curl), thus necessitating this term.
Lemma 19**.**
For any there is a such that
[TABLE]
Proof.
We only present the proof in two dimensions, as the three dimensional case is similar. From the local element basis exhibited in (20) (see also [18, §5.5] for a more detailed discussion), its clear that on any facet , there exists a constant trace-free function with the property that , on the facet and equals on all other facets in . Given any , define
[TABLE]
where is the unique barycentric coordinate function on the element opposite to the facet (so that is an -bubble). Clearly, and are in . Using the norm equivalences stated in (26) and the mappings for and given in (10), a scaling argument yields
[TABLE]
Setting and selecting the constants appropriately, the rest of the proof proceeds along the same lines as the proof of [18, Lemma 6.5]. ∎
Remark 20*.*
It is interesting to contrast Lemma 19 with [18, Lemma 6.5]. The latter gives a similar LBB-condition. The differences are (i) the velocity space in [18] is (defined in Remark 12), (ii) the velocity norm is a discrete -norm defined using in place of , (iii) there is no weak symmetry constraint and no associated space , and (iv) the stress space in [18] equals the in (12) plus certain -bubbles of degree (different from our here). Lemma 19 shows that the inf-sup condition in [18, Lemma 6.5] continues to hold even if the -bubbles there are removed and is replaced by our Raviart-Thomas velocity space . This observation can be extended to prove the convergence of the MCS formulation in [18] with so modified spaces.
Theorem 21** (Discrete LBB-condition).**
Let and . Then,
[TABLE]
If is in the divergence-free subspace then
[TABLE]
Proof.
By Lemmas 17 and 19, for any given , there are satisfying
[TABLE]
Clearly, the same inequalities hold when and are scaled by any nonzero factor, so we may assume without loss of generality, that they have been scaled so that and Set where is to be chosen shortly. It follows from (53) and (54) that
[TABLE]
Next, we choose so that , where is another constant to be chosen shortly. Then (55) implies
[TABLE]
Choose any and . Then, using Young’s inequality for the last term,
[TABLE]
we can now conclude the proof of (51) using the norm equivalence of Lemma 14. The proof of (52) is similar (and in fact simpler since all terms involving vanish). ∎
4.4. Error estimates
In this subsection we show that the error in the discrete MCS solution converges at optimal order. As we have chosen polynomials of degree for the stress space , the optimal rate of convergence for is . However, the optimal rate for the velocity error in our discrete -like norm, namely, is only (since the Raviart-Thomas velocity space only contains within each mesh element ). Nevertheless, we are still able to prove optimal convergence rate of the stress error by using an appropriate interpolation operator and deducing that the stress error is independent of the velocity error. Another important property we shall conclude in this subsection is the pressure-robustness of the method.
Lemma 22** (Continuity).**
The bilinear forms and are continuous:
[TABLE]
Proof.
The continuity of and follow by the Cauchy Schwarz inequality. For we use (18) and to get
[TABLE]
Now, Cauchy-Schwarz inequality and (26) of Lemma 8 finishes the proof. ∎
Lemma 23** (Coercivity in the kernel).**
For all in the kernel
[TABLE]
we have \nu^{-1}\big{(}\,\left\|\tau_{h}\right\|+\|q_{h}\|\big{)}^{2}\lesssim\;a(\tau_{h},\tau_{h}).
Proof.
By [24, Theorem 2.2], for any there is a such that and a discrete -norm of is bounded by . The latter bound implies, in particular, that and also that satisfies This together with Lemma 22 implies
[TABLE]
yielding the needed bound for . ∎
We are now ready to conclude an inf-sup condition for
Corollary 24**.**
Let , , , and . There holds
[TABLE]
so, in particular, there is a unique solution for the discrete MCS system (19). Moreover, if is restricted to , we also have
[TABLE]
Proof.
The first inf-sup condition follows from the standard theory of mixed methods [4], using Theorem 21 (the inf-sup condition for and given by (51)), Lemma 22 (continuity of forms), and Lemma 23 (coercivity in the kernel).
The second inf-sup condition also follows in a similar fashion, but now using the other inequality (52) of Theorem 21. ∎
Theorem 25** (Consistency).**
The MCS method with weakly imposed symmetry (19) is consistent in the following sense. If the exact solution of the Stokes problem (9) is such that , , and , then
[TABLE]
*for all and *
The proof of Theorem 25 is easy (see, e.g., the similar proof of [18, Theorem 6.2]), so we omit it. We now have all the ingredients to prove the following convergence result. Let denote the standard Raviart-Thomas interpolator (see, e.g., [4]) and let
Theorem 26** (Optimal convergence).**
Let , , and be the exact solution of the mixed Stokes problem (9), let ,, and solve (19) and let . Then,
[TABLE]
Proof.
Let , , , (where the two occurrences of represent projections onto two different discrete spaces per our prior notation). Denoting the analogous approximation errors by , , , and , observe that Theorem 25 implies
[TABLE]
for any and . The right hand side above is a sum of five terms The second term vanishes: as . The third term also vanishes: since . The fourth term, due to (17), is
[TABLE]
where the last two terms vanish by the properties of the Raviart-Thomas d.o.f.s that define , i.e., The fifth term, due to (18), is
[TABLE]
Writing note that by the d.o.f.s of Theorem 5, the last term is zero, and moreover, . Incorporating these observations on each term into (59), we obtain
[TABLE]
We now proceed to estimate the right hand side of (60). By (42) and Lemma 11,
[TABLE]
Using these after an application of the Cauchy-Schwarz inequality, (60) yields
[TABLE]
where we have used Theorem 7 and the approximation property of in the last step.
To complete the proof, we apply triangle inequality starting from the left hand side of (58), to get
[TABLE]
again using Theorem 7. Bounding the last term above using (56) and (61),
[TABLE]
the proof is complete. ∎
Remark 27* (Convergence in standard norms).*
Using also Lemma 15’s estimate (43), a consequence of the global discrete Korn inequality, (58) implies
[TABLE]
under the assumptions of Theorem 26 for a sufficiently smooth solution. Note that even though the optimal rate for is only , (63) gives a superconvergent rate of for .
Theorem 28** (Pressure robustness).**
Under the same assumptions as Theorem 26,
[TABLE]
Proof.
Proceeding along the lines of the proof of Theorem 26, omitting the pressure error, we obtain, instead of (62),
[TABLE]
We may now complete the proof as before by using (57) instead of (56). ∎
5. Postprocessing
In this section we describe and analyze a postprocessing for the discrete velocity. While for the raw solution , we may only expect to go to zero at the rate we will show that a locally postprocessed velocity has error that converges to zero at the higher rate for sufficiently regular solutions. The key to obtain this enhanced accuracy, as in [30], is the -superconvergence of – see Remark 27. Finally, we shall also show that retains the prized structure preservation properties of exact mass conservation and pressure robustness.
The crucial ingredient is a reconstruction operator (see [21, 22]) whose properties are summarized in the next lemma. Let
[TABLE]
denote the BDM space (one order higher) and its “relaxed” analogue, respectively. The next result is a consequence of [21, Lemmas 3.3 and 4.8] and the Korn inequality (29).
Lemma 29**.**
There exists an operator whose application is computable element-by-element, satisfying
- (1)
* for al ,* 2. (2)
* for all , and* 3. (3)
whenever the local (element-wise) property holds for all and all , the global property holds.
A simple choice of is given by the classical BDM intepolant. This was used in [19]. Another choice of , given in [21], based on a simple averaging of coefficients, is significantly less expensive for high orders.
The postprocessed solution is given in two steps as follows. First, using the computed and , solve the local (see Remark 31) minimization problem
[TABLE]
Second, apply the reconstruction and set .
Theorem 30**.**
Suppose the assumptions of Theorem 26 hold. Then and for we have the pressure-robust error estimate
[TABLE]
Proof.
On any , the condition implies that the Raviart-Thomas d.o.f.s applied to and coincide. Hence, for all ,
[TABLE]
as . Thus, Lemma 29 implies that and .
It only remains to prove the error estimate. Let be the standard interpolator. Then, satisfies
[TABLE]
Since standard approximation estimates yield , we focus on the last term. A triangle inequality (where we add and subtract different functions in the element and facet terms) yields
[TABLE]
Naming the four sums on the right as and , respectively, we proceed to estimate each. Obviously by Theorem 28.
To bound , note that for any in the admissible set of the minimization problem (64), we have . We choose . Since implies , the chosen is in the admissible set. Hence,
[TABLE]
so a standard approximation estimate and Theorem 28 yield
The same standard approximation estimate for also gives . Hence it only remains to bound . Observe that because and . This implies, letting , the identity holds because (as ). Hence
[TABLE]
Since the first term can be bounded by Theorem 28, let us consider the last term. On any facet adjacent to a mesh element , a trace inequality yields h^{-1}\big{\|}{[\![a_{t}]\!]}\big{\|}_{F}^{2}\leq h^{-1}\|a_{t}\|_{\partial T}^{2}\lesssim\|\nabla a\|_{T}^{2}+h^{-2}\|a\|_{T}^{2}. Hence,
[TABLE]
where we have used the continuity properties of , scaling arguments, (27), and an estimate analogous to (28). Using triangle inequality and returning to (66),
[TABLE]
The last two terms are and , respectively. Hence the prior estimates, the standard approximation estimate for , and Theorem 28 shows ∎
Remark 31*.*
The restriction of the minimizer of (64) to an element , namely can be computed using the following Euler-Lagrange equations. Letting on all facets , the function is the unique function in , which together with and , satisfies
[TABLE]
for all and . The last two equations are another way to express the constraint in (64).
6. Numerical exampels
In this last section we present two numerical examples to verify our method. All examples were implemented within the finite element library NGSolve/Netgen, see [28, 29] and on www.ngsolve.org. The computational domain is given by and the velocity field is driven by the volume force determined by with the exact solution given by
[TABLE]
Here and defines a given potential in two and three dimensions respectively and we choose the viscosity .
In Tables 1(a) and 1(b) we report the errors in all the computed solution components for varying polynomial orders in the two and the three dimensional cases, respectively. As predicted by Theorem 26 and Theorem 30 the corresponding errors converge at optimal order. Furthermore, the -norm of error of the (postprocessed) velocity error converges at one order higher. Note that in three dimensions the errors are already quite small already on the coarsest mesh. It appears that to get out of the preasymptotic regime and see the proper convergence rate, it takes several steps.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] D. N. Arnold, F. Brezzi, and J. Douglas, Jr. , PEERS: a new mixed finite element for plane elasticity , Japan J. Appl. Math., 1 (1984), pp. 347–367.
- 2[2] D. N. Arnold, R. S. Falk, and R. Winther , Mixed finite element methods for linear elasticity with weakly imposed symmetry , Math. Comp., 76 (2007), pp. 1699–1723.
- 3[3] D. Boffi, F. Brezzi, and M. Fortin , Reduced symmetry elements in linear elasticity , Commun. Pure Appl. Anal., 8 (2009), pp. 95–121.
- 4[4] , Mixed Finite Element Methods and Applications , Springer Science & Business Media, 2013.
- 5[5] S. C. Brenner , Korn’s inequalities for piecewise H 1 superscript 𝐻 1 H^{1} vector fields , Math. Comp., 73 (2004), pp. 1067–1087.
- 6[6] F. Brezzi, J. Douglas Jr., and L. D. Marini , Two families of mixed finite elements for second order elliptic problems , Numerische Mathematik, 47 (1985), pp. 217–235.
- 7[7] B. Cockburn and J. Gopalakrishnan , A characterization of hybridized mixed methods for the Dirichlet problem , SIAM J. Numer. Anal., 42 (2004), pp. 283–301.
- 8[8] B. Cockburn, J. Gopalakrishnan, and J. Guzmán , A new elasticity element made for enforcing weak stress symmetry , Math. Comp., 79 (2010), pp. 1331–1349.