MCS formulation with weakly imposed symmetry

# 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

 (1)

where . By introducing a new variable where , equation (1) can be reformulated to

 (2a) 1νdev(σ)−ε(u) =0 in Ω, (2b) div(σ)−∇p =−f in Ω, (2c) div(u) =0 in Ω, (2d) u =0 on Γ.

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 and , respectively. Let denote the vector space of skew symmetric matrices, i.e., , and let .

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.

 curl(ϕ) =(−∂2ϕ,∂1ϕ)T, for ϕ∈D∗(Ω,R),d=2, curl(ϕ) =(∂2ϕ3−∂3ϕ2,∂3ϕ1−∂1ϕ3,∂1ϕ2−∂2ϕ1)T for ϕ∈D∗(Ω,R3),d=3,

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 ,

 (3) curl(ϕ) =−∂2ϕ1+∂1ϕ2, or (4) curl(ϕ) =(∂2ϕ1−∂1ϕ1∂2ϕ2−∂1ϕ2),

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

 (5) [H0(div,Ω)]∗=H−1(curl,Ω)={ϕ∈H−1(Ω,Rd):curl(ϕ)∈H−1(Ω,R~d)}.

In this work, the following space is important:

 H(curldiv,Ω) :={σ∈L2(Ω,M):div(σ)∈[H0(div,Ω)]∗},

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

 V :=H0(div,Ω), Q :=L20(Ω):={q∈L2(Ω):∫Ωq dx=0}.

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

 (6) dev(ν−1σ)=dev(ε(u))=ε(u)−νdtr(ε(u))Id=ε(u)−1ddiv(u)Id=ε(u).

Since and , we define the stress space as the following closed subspace of :

 Σsym:={τ∈H(curldiv,Ω):tr(τ)=0,τ=τT}.

Testing equations (2a) with a test functions and integrating over the domain, we have for the term including the identity

 ∫Ωε(u):τ dx =12∫Ω∇u:τ dx+12∫Ω(∇u)T:τ dx =12∫Ω∇u:τ dx+12∫Ω∇u:τ dx=∫Ω∇u:τ dx.

Using the knowledge that the velocity should be in we obtain

 (ν−1dev(σ),dev(τ))+⟨div(τ),u⟩H0(div,Ω) =0,

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

 (7) ⎧⎪ ⎪⎨⎪ ⎪⎩(ν−1dev(σ),dev(τ))+⟨div(τ),u⟩H0(div,Ω)=0 for all τ∈Σsym,⟨div(σ),v⟩H0(div,Ω)+(div(v),p)=−(f,v) for all v∈V,(div(u),q)=0 for all p∈Q.

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

 (8)

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

 (9a) 1νdev(σ)−∇u+ω =0 in Ω, (9b) div(σ)−∇p =−f in Ω, (9c) σ−σT =0 in Ω, (9d) div(u) =0 in Ω, (9e) u =0 on Γ.

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:

 Hm(Th):=∏T∈ThHm(T),Pk(Th):=∏T∈ThPk(T).

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

 Vh:=V∩RTk,Qh:=Q∩Pk(Th),Wh:=Pk(Th,K),

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

 (10) Q(qh)=^qh∘ϕ−1,P(^vh):=det(F)−1F(^vh∘ϕ−1),W(^ηh):=F−T(^ηh∘ϕ−1)F−1,

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

 (11) ⟨div(σ),v⟩H0(div,Ω)=∑T∈Th[(div(σ),v)T−⟨vn,σnn⟩H1/2(∂T)]

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

 (12) Σh :={τh∈Pk(Th,M):tr(τh)=0,[[(τh)nt]]=0 for all % F∈Finth}.

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

 (13a) B =3∑i=0λi−3λi−2λi−1∇λi⊗∇λi if d=3, (13b) B =λ0λ1λ2 if d=2,

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

 (14) δΣh:={dev(curl(curl(rh)B)):rh∈Pk⊥(Th,K)},

for and , with the understanding that in case, the outer curl is defined by (4), not (3). The total stress space is given by

 Σ+h:=Σh⊕δΣh,k≥1.

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]:

 (15) M(^σh):=1det(F)F−T(^σh∘ϕ−1)FT.

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

 (16) B:=F−T(^B∘ϕ−1)F−1.

### 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:

 a:L2(Ω,M)×L2(Ω,M)→R, b1:V×Q→R, a(σ,τ):=(ν−1dev(σ),dev(τ)), b1(u,p):=(div(u),p).

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

 (17) b2(τ,(v,η)) :=∑T∈Th∫Tdiv(τ)⋅v dx+∑T∈Th∫Tτ:η dx−∑F∈Fh∫F[[τnn]]vn ds.

Integrating the first integral by parts, we find the equivalent representation

 (18) b2(τ,(v,η)) =−∑T∈Th∫Tτ:(∇v−η) dx+∑F∈Fh∫Fτnt⋅[[vt]] ds.

Using these forms, we state the method. For any , the discrete MCS method with weakly imposed symmetry finds such that

 (19) ⎧⎪⎨⎪⎩a(σh,τh)+b2(τh,(uh,ωh))=0 % for all τh∈Σ+h,b2(σh,(vh,ηh))+b1(vh,ph)=(−f,vh) for all (vh,ηh)∈Uh:=Vh×Wh,b1(uh,qh)=0 % for all qh∈Qh.

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

 (20a) Si:=dev(∇λi+1⊗curl(λi+2)), (20b) Si0:=dev(∇λi+1⊗(∇λi+2×∇λi+3)),Si1:=dev(∇λi+2⊗(∇λi+3×∇λi+1)),

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

 ΦF(τ):=∫Fτnt⋅r ds

for each in any fixed basis for . The next group is the set of interior d.o.f.s, defined by

 Φ0(τ):=∫Tτ:ς dx

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

 (21) τ=2∑i=0μiλiSiorτ=1∑q=03∑i=0μqiλiSiq,

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

 ∫T2∑i=0μiλiSi:2∑i=0μiSi dx=∫Tλi∣∣∣2∑i=0μiSi∣∣∣2 dx=0,

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 , equalling . ∎

Using these d.o.f.s, a canonical local interpolant in can be defined as usual, by requiring that for all

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

 ∥σ−IΣhσ∥2+∑F∈Fhh∥(σ−IΣhσ)nt∥2F≲h2s∥σ∥2Hs(Th),

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.,

 ∥ε(v)∥2h:=∑T∈Th∥ε(v)∥2T,∥curl(γ)∥2h:=∑T∈Th∥curl(γ)∥2T, ∥v∥21,h,ε:=∥ε(v)∥2h+∑F∈Fh1h∥∥[[vt]]∥∥2F,

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 :

 ∥vh∥2Vh=∥vh∥21,h,ε,∥(vh,ηh)∥2Uh:=∥vh∥21,h,ε+∥κ(curlvh)−ηh∥2h.

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

 (22) ∥(vh,ηh)∥2Uh,∗:=∑T∈Th∥Πk−1Tdev(∇vh−ηh)∥2T+∑F∈Fh1h∥Π1F[[(vh)t]]∥2F.

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

 (23) ∥(vh,ηh,τh,qh)∥∗:=√ν||(vh,ηh)||Uh+1√ν(∥τ