Analysis and approximation of a vorticity-velocity-pressure formulation for the Oseen equations

# Analysis and approximation of a vorticity-velocity-pressure formulation for the Oseen equations

Verónica Anaya,  Afaf Bouharguane,   David Mora,   Carlos Reales,
Ricardo Ruiz-Baier,  Nour Seloula,  and Hector Torres
GIMNAP, Departamento de Matemática, Universidad del Bío-Bío, Concepción, Chile. E-mail: vanaya@ubiobio.cl.Institut de Mathématiques de Bordeaux, CNRS UMR 5251, Université de Bordeaux, 33405 Talence, France. E-mail: afaf.bouharguane@math.u-bordeaux.fr.GIMNAP, Departamento de Matemática, Universidad del Bío-Bío, Concepción, Chile and Centro de Investigación en Ingeniería Matemática (CIMA), Universidad de Concepción, Concepción, Chile. E-mail: dmora@ubiobio.cl.Departamento de Matemáticas y Estadísticas, Universidad de Córdoba, Montería, Colombia. E-mail: creales@correo.unicordoba.edu.co.Mathematical Institute, University of Oxford, OX2 6GG Oxford, UK. E-mail: ruizbaier@maths.ox.ac.uk.LMNO, CNRS UMR 6139, Université de Caen, 5186 Caen, France. E-mail: nour-elhouda.seloula@unicaen.fr.Departamento de Matemáticas, Universidad de La Serena, La Serena, Chile. E-mail: htorres@userena.cl.
July 13, 2019
###### Abstract

We introduce a family of mixed methods and discontinuous Galerkin discretisations designed to numerically solve the Oseen equations written in terms of velocity, vorticity, and Bernoulli pressure. The unique solvability of the continuous problem is addressed by invoking a global inf-sup property in an adequate abstract setting for non-symmetric systems. The proposed finite element schemes, which produce exactly divergence-free discrete velocities, are shown to be well-defined and optimal convergence rates are derived in suitable norms. In addition, we establish optimal rates of convergence for a class of discontinuous Galerkin schemes, which employ stabilisation. A set of numerical examples serves to illustrate salient features of these methods.

Key words: Oseen equations; vorticity-based formulation; mixed finite elements; exactly divergence-free velocity; discontinuous Galerkin schemes; numerical fluxes; a priori error bounds.

Mathematics subject classifications (2000): 65N30, 65N12, 76D07, 65N15.

## 1 Introduction

The Oseen equations stem from linearisation of the steady (or alternatively from the backward Euler time-discretisation of the transient) Navier-Stokes equations. Of particular appeal to us is their formulation in terms of fluid velocity, vorticity vector, and pressure. A diversity of discretisation methods is available to solve incompressible flow problems using these three fields as principal unknowns. Some recent examples include spectral elements [3, 8] as well as stabilised and least-squares schemes [2, 9] for Navier-Stokes; also several mixed and augmented methods for Brinkman [4, 5, 7], and a number of other discretisations specifically designed for Stokes flows [6, 22, 23, 25, 30].

Both the implementation and the analysis of numerical schemes for Navier-Stokes are typically based on the Oseen linearisation. A few related contributions (not only restricted to the velocity-pressure formulation) include for instance [10], that presents a least-squares method for Navier-Stokes equations with vorticity-based first-order reformulation, and whose analysis exploits the elliptic theory of Agmon-Douglas-Nirenberg. Conforming finite element methods exhibit optimal order of accuracy for diverse boundary conditions. We also mention the non-conforming exponentially accurate least-squares spectral method for Oseen equations proposed in [27], where a suitable preconditioner is also proposed. In [32] the authors introduce a velocity-vorticity-pressure least-squares finite element method for Oseen and Navier-Stokes equations with velocity boundary conditions. They derive error estimates and reported a degeneracy of the convergence for large Reynolds numbers. A div least-squares minimisation problem based on the stress-velocity-pressure formulation was introduced in [13]. The study shows that the corresponding homogeneous least-squares functional is elliptic and continuous in suitable norms. Several first-order Oseen-type systems are analysed in [15], also including vorticity and total pressure in the formulation.

Discontinuous Galerkin (DG) methods have also been used to solve the Oseen problem, as for example, in [19, 18] for the case of Dirichlet boundary conditions. Compared with conforming finite elements, discretisations based on DG methods have a number of attractive, and well-documented features. These include high order accuracy, being amenable for -adaptivity, relatively simple implementation on highly unstructured meshes, and superior robustness when handling rough coefficients. We also mention the a priori error analysis of hybridisable DG schemes introduced in [16] for the Oseen equations. The family of DG methods we propose here has resemblance with those schemes, but exploits a three-field formulation described below.

This paper is concerned with mixed non-symmetric variational problems which will be analysed using a global inf-sup argument. To do this, we conveniently restrict the set of equations to the space of divergence-free velocities, and apply results from [24] in order to prove that the equivalent resulting non-symmetric saddle-point problem is well-posed. For the numerical approximation, we first consider Raviart-Thomas elements of order for the velocity field, Nédélec elements or order for the vorticity, and piecewise polynomials of degree without continuity restrictions, for the Bernoulli pressure. We prove unique solvability of the discrete problem by adapting the same tools utilised in the analysis of the continuous problem. In addition, the proposed family of Galerkin finite element methods turns out to be optimally convergent, under the common assumptions of enough regularity of the exact solutions to the continuous problem. The method produces exactly divergence-free approximations of the velocity by construction; thus preserving, at the discrete level, an essential constraint of the governing equations. Next, inspired by the methods presented in [20, 19], we present another scheme involving the discontinuous Galerkin discretisation of the - and - operators. We prove the well-posedness of the DG scheme and derive error estimates under some solution regularity assumptions.

We have structured the contents of the paper in the following manner. Notation-related preliminaries are stated in the remainder of this Section. We then present the model problem as well as the three-field weak formulation and its solvability analysis in Section 2. The finite element discretisation is constructed in Section 3, where we also derive the stability and convergence bounds. In Section 4, we present the mixed DG formulation for the model problem. The well-posedness of the method and the error analysis are established in the same section. We close in Section 5 with a set of numerical tests that illustrate the properties of the proposed numerical schemes in a variety of scenarios.

Let be a bounded domain of with Lipschitz boundary . Moreover, we assume that admits a disjoint partition . For any , the symbol denotes the norm of the Hilbertian Sobolev spaces or , with the usual convention . For , we recall the definition of the space

 Hs(curl;Ω):={θ∈Hs(Ω)3: curlθ∈Hs(Ω)3},

endowed with the norm , and will denote . Finally, and , with subscripts, tildes, or hats, will represent a generic constant independent of the mesh parameter .

## 2 Statement and solvability of the continuous problem

### Oseen problem in terms of velocity-vorticity-pressure.

A standard backward Euler time-discretisation of the classical Navier-Stokes equations, or a linearisation of the steady version of the problem combined with standard curl-div identities, leads to the following set of equations, known as the Oseen equations (see [29, 26]):

 σu−νΔu+curlu×β+∇p=f in Ω,divu=0 in Ω, (2.1)

where is the kinematic fluid viscosity, is inversely proportional to the time-step, is an adequate approximation of velocity to be made precise below, and the vector of external forces also absorbs the contributions related to previous time steps, or to fixed states in the linearisation procedure of the steady Navier-Stokes equations. As usual in this context, in the momentum equation we have conveniently introduced the Bernoulli (also known as dynamic) pressure , where is the actual fluid pressure.

The structure of (2.1) suggests to introduce the rescaled vorticity vector as a new unknown. Furthermore, in this study we focus on the case of zero normal velocities and zero tangential vorticity trace imposed on a part of the boundary , whereas a non-homogeneous tangential velocity and a fixed Bernoulli pressure are set on the remainder of the boundary . Therefore, system (2.1) can be recast in the form

 σu+√νcurlω+ν−1/2ω×β+∇p=f,ω−√νcurlu=0,anddivu =0 in Ω, u⋅n=0 and ω×n =0 on Γ, (2.2) p=pΣ and u×n =uΣ on Σ,

where stands for the outward unit normal on . Should the boundary have zero measure, the additional condition is required to enforce uniqueness of the Bernoulli pressure.

### Defining a weak formulation.

Let us introduce the following functional spaces

 H:={v∈H(div;Ω):v⋅n=0onΓ},Z:={θ∈H(curl;Ω):γt(θ)=0onΓ},and Q:=L2(Ω),

where the operator is the tangential trace operator on , defined by: . Let us endow and with their natural norms. For the space however, we consider the following viscosity-weighted norm:

 ∥θ∥Z:=(∥θ∥20,Ω+ν∥curlθ∥20,Ω)1/2.

From now on, we will assume that the data are regular enough: and . We proceed to test (2.2) against adequate functions and to impose the boundary conditions in such a manner that we end up with the following formulation: Find such that

 a(u,v) + b1(v,ω) +b2(v,p)+c(ω,v) =F(v) ∀v∈H, b1(u,θ) − d(ω,θ) =G(θ) ∀θ∈Z, (2.3) b2(u,q) =0 ∀q∈Q,

where the bilinear forms , , , , , and the linear functionals , and are specified as follows

 a(u,v) :=σ∫Ωu⋅vdx,b1(v,θ):=√ν∫Ωcurlθ⋅vdx,b2(v,q):=−∫Ωqdivvdx, d(ω,θ) :=∫Ωω⋅θdx,c(θ,v):=1√ν∫Ω(θ×β)⋅vdx, F(v) :=∫Ωf⋅vdx−⟨v⋅n,pΣ⟩Σ,G(θ):=−√ν⟨uΣ,θ⟩Σ,

for all , , and .

### Solvability analysis.

In order to analyse the variational formulation (2.3), let us introduce the Kernel of the bilinear form and its classical characterisation

 X:={v∈H:b2(v,q)=0,∀q∈Q}={v∈H:divv=0inΩ},

and let us recall that satisfies the inf-sup condition:

 sup\lx@stackrelv∈Hv≠0|b2(v,q)|∥v∥H≥β2∥q∥0,Ω∀q∈Q, (2.4)

with an inf-sup constant only depending on (see e.g. [24]).

We will now address the well-posedness of (2.3). To that end, it is enough to study its reduced counterpart, defined on : Find such that

 a(u,v) + b1(v,ω) +c(ω,v) =F(v) ∀v∈X, b1(u,θ) − d(ω,θ) =G(θ) ∀θ∈Z. (2.5)

The equivalence between (2.3) and (2.5) is established in the following result, whose proof follows [26] and it is basically a direct consequence of the inf-sup condition (2.4).

###### Lemma 2.1

If is a solution of (2.3), then and also solves (2.5). Conversely, if is a solution of (2.5), then there exists a unique such that solves (2.3).

The abstract setting that will permit the analysis of (2.3) is stated in the following general result [24, Theorem 1.2].

###### Theorem 2.1

Let be a bounded bilinear form and a bounded functional, both defined on the Hilbert space . If there exists such that

 supy∈X∖{0}A(x,y)∥y∥X≥α∥x∥X∀x∈X, (2.6)

and

 sup\lx@stackrelx∈Xy≠0A(x,y)>0∀y∈X, (2.7)

then there exists a unique solution to the problem

 A(x,y)=G(y)∀y∈X.

Furthermore, there exists (independent of ) such that

 ∥x∥X≤1α∥G∥X′.
###### Lemma 2.2

Let us assume that

 2∥β∥2∞,Ωνσ<1, (2.8)

and let us define the bilinear form specified as

 A((u,ω),(v,θ)):=a(u,v)+b1(v,ω)+b1(u,θ)−d(ω,θ)+c(ω,v).

Then, there exist such that

 |A((u,ω),(v,θ))|≤α1∥(u,ω)∥X∥(v,θ)∥X, (2.9)

and

 sup\lx@stackrel(v,θ)∈X(v,θ)≠0A((u,ω),(v,θ))∥(v,θ)∥X≥α2∥(u,ω)∥X∀(u,ω)∈X, (2.10)

where , endowed with the corresponding product norm, is a Hilbert space.

Proof. As a consequence of the boundedness of , , and , the bilinear form is bounded and so condition (2.9) readily follows.

Concerning the satisfaction of the inf-sup condition (2.10), for a given , we can define

 ~θ:=−ω∈Z,and~v:=(u+^c√νcurlω)∈X,

where is a constant to be chosen later. We can then immediately assert that

 A((u,ω),(~v,~θ)) =σ∫Ωu⋅~vdx+√ν∫Ωcurlω⋅~vdx+√ν∫Ωcurl~θ⋅udx −∫Ωω⋅~θdx+1√ν∫Ω(ω×β)⋅~vdx ≥σ∥u∥20,Ω+^c√νσ∫Ωu⋅curlωdx+√ν∫Ωu⋅curlωdx+^cν∥curlω∥20,Ω −√ν∫Ωu⋅curlωdx+∥ω∥20,Ω+1√ν∫Ω(ω×β)⋅udx+^c∫Ω(ω×β)⋅curlωdx ≥σ∥u∥20,Ω−σ4∥u∥20,Ω−^c2σν∥curlω∥20,Ω+^cν∥curlω∥20,Ω+∥ω∥20,Ω−2σ3∥u∥20,Ω −3∥β∥2∞,Ω2νσ∥ω∥20,Ω−∥β∥2∞,Ω2νσ∥ω∥20,Ω−2^c2σν∥curlω∥20,Ω =σ12∥u∥20,Ω+^c(1−3^cσ)ν∥curlω∥20,Ω+(1−2∥β∥2∞,Ωνσ)∥ω∥20,Ω,

where we have used the bound . Choosing and exploiting (2.8), we arrive at

 A((u,ω),(~v,~θ))≥C∥(u,ω)∥2X,

with independent of . On the other hand, by construction we realise that and , and consequently

 sup\lx@stackrel(v,θ)∈X(v,θ)≠0A((u,ω),(v,θ))∥(v,θ)∥X≥A((u,ω),(~v,~θ))∥(~v,~θ)∥X≥α2∥(u,ω)∥X∀(u,ω)∈X,

which finishes the proof.

###### Lemma 2.3

Suppose that the bound (2.8) is satisfied. Then, there exists such that

 sup\lx@stackrel(u,ω)∈X(u,ω)≠(0,0)A((u,ω),(v,θ))>0∀(v,θ)∈X.

Proof. For all , we have that:

 A((v,−θ),(v,θ)) =σ∥v∥20,Ω+∥θ∥20,Ω−1√ν∫Ω(θ×β)⋅vdx ≥σ∥v∥20,Ω+∥θ∥20,Ω−σ2∥v∥20,Ω−2∥β∥2∞,Ωνσ∥θ∥20,Ω ≥σ2∥v∥20,Ω+(1−2∥β∥2∞,Ωνσ)∥θ∥20,Ω.

As a consequence of the previous lemmas, we have the following result.

###### Theorem 2.2

Let us assume (2.8). Then, the variational problem (2.5) admits a unique solution . Moreover, there exists such that

 ∥u∥H+∥ω∥Z≤C(∥f∥0,Ω+∥pΣ∥1/2,Σ+∥uΣ∥−1/2,Σ). (2.11)

Proof. It suffices to verify the hypotheses of Theorem 2.1. First, we define the linear functional

 G(v,θ):=F(v)+G(θ),

which is bounded on . Thus, the proof follows from Lemmas 2.2 and 2.3.

The following result establishes the corresponding stability estimate for the Bernoulli pressure.

###### Corollary 1

Let , be the unique solution of (2.5), with and satisfying (2.11). In addition, let be the unique pressure provided by Lemma 2.1, so that is the unique solution of (2.3). Then, there exists such that

 ∥p∥0,Ω≤C(∥f∥0,Ω+∥pΣ∥1/2,Σ+∥uΣ∥−1/2,Σ).

Proof. Combining the inf-sup condition (2.4) with the first equation in (2.3) gives the bound

 ∥p∥0,Ω≤1β2sup\lx@stackrelv∈Hv≠0|b2(v,p)|∥v∥H=1β2sup\lx@stackrelv∈Hv≠0|F(v)−a(u,v)−b1(v,ω)−c(ω,v)|∥v∥H,

which together with (2.11), and the boundedness of , , and , complete the proof.

###### Remark 2.1

An alternative analysis for the nonsymmetric variational problem (2.3) can be carried out using a fixed-point argument that allows a symmetrisation of the mixed structure. The resulting weak form could then be analysed using classical tools for saddle-point problems, for instance, following the similar treatment in [14]. Establishing inf-sup conditions for the off-diagonal bilinear forms in the original nonsymmetric formulation is, however, much more involved (see e.g. [28]).

###### Remark 2.2

Assumption (2.8) holds provided one chooses appropriately. As this parameter represents the inverse of the timestep, the aforementioned relation constitutes then a CFL-type condition. We also note that this bound for coincides with the hypotheses that yield solvability of least-squares formulations for the Oseen problem analysed in [15, 13].

## 3 Finite element discretisation

In this section we introduce a Galerkin scheme for (2.3) and analyse its well-posedness by establishing suitable assumptions on the finite element subspaces involved. Error estimates are also derived.

### Defining the discrete problem.

Let be a shape-regular family of partitions of the polyhedral region , by tetrahedrons of diameter , with mesh size . In what follows, given an integer and a subset of , will denote the space of polynomial functions defined locally in and being of total degree .

Moreover, for any , we introduce the local Nédélec space

 Nk(T):=Pk(T)3⊕Rk+1(T),

where is a subspace of composed by homogeneous polynomials of degree , and being orthogonal to . With these tools, let us define the following finite element subspaces:

 Zh :={θh∈Z:θh|T∈Nk(T)∀T∈Th(Ω)}, (3.1) Qh :={qh∈Q:qh|T∈Pk(T)∀T∈Th(Ω)}, (3.2) Hh :={vh∈H:vh|T∈RTk(T)∀T∈Th(Ω)}, (3.3)

where is the Raviart-Thomas space defined locally in .

The proposed Galerkin scheme approximating (2.3) reads as follows: Find such that

 a(uh,vh) + b1(vh,ωh) +b2(vh,ph)+c(ωh,vh) =F(vh) ∀vh∈Hh, b1(uh,θh) − d(ωh,θh) =G(θh) ∀θh∈Zh, (3.4) b2(uh,qh) =0 ∀qh∈Qh.

### Solvability and stability of the discrete problem.

The analysis of the Galerkin formulation will follow the same arguments exploited in the continuous setting. Let us then consider the discrete kernel of :

 Xh:={vh∈Hh:b2(vh,qh)=0,∀q∈Qh}={vh∈Hh:divvh≡0inΩ}, (3.5)

where the characterisation is indeed possible thanks to the inclusion . Moreover, it is well-known that the following discrete inf-sup condition holds (see [24]):

 sup\lx@stackrelvh∈Hhvh≠0b2(vh,qh)∥vh∥H≥~β2∥qh∥0,Ω∀qh∈Qh. (3.6)

We again resort to a reduced version of the problem, now defined on the product space . Find such that

 a(uh,vh) + b1(vh,ωh)+c(ωh,vh) = F(vh) ∀vh∈Xh, b1(uh,θh) − d(ωh,θh) = G(θh) ∀θh∈Zh, (3.7)

and its equivalence with (3.4) is once more a direct consequence of the inf-sup condition (3.6).

###### Lemma 3.1

If is a solution of (3.4), then , and is also a solution of (3.7). Conversely, if is a solution of (3.7), then there exists a unique such that is a solution of (3.4).

In order to establish the well-posedness of (3.7), we will employ the following discrete version of Theorem 2.1.

###### Theorem 3.1

Assume (2.8). Let be an integer and let and be given by (3.5) and (3.1), respectively. Then, there exists a unique solution of the discrete scheme (3.7). Moreover, there exist positive constants independent of such that

 ∥uh∥H+∥ωh∥Z≤^C1(∥f∥0,Ω+∥pΣ∥1/2,Σ+∥uΣ∥−1/2,Σ), (3.8)

and

 ∥u−uh∥H+∥ω−ωh∥Z≤^C2inf(vh,θh)∈Xh×Zh(∥u−vh∥H+∥ω−θh∥Z), (3.9)

where is the unique solution to (2.5).

Proof. Let us define and reuse the forms and as in the proof of Lemma 2.2. The next step consists in proving that satisfies the discrete version of the inf-sup conditions (2.6)-(2.7), as in Lemmas 2.2 and  2.3. In order to assert (2.6), we consider , and define

 ~θh:=−ωh∈Zh,and ~vh:=(uh+√ν4σcurlωh)∈Xh.

Then, repeating exactly the same steps used in the proof of Lemma 2.2 the discrete version of (2.6) follows. Regarding the discrete version of (2.7), we once again repeat the same arguments given in the proof of Lemma 2.3. Finally, the Céa estimate follows from classical arguments.

We now state the stability and an adequate approximation property of the discrete pressure.

###### Corollary 2

Let be the unique solution of (3.7), with and satisfying (3.8). In addition, let be the unique discrete Bernoulli pressure provided by Lemma 3.1, so that is the unique solution of (3.4). Then, there exist positive constants , independent of and , such that

 ∥ph∥0,Ω≤¯C1(∥f∥0,Ω+∥pΣ∥1/2,Σ+∥uΣ∥−1/2,Σ),

and

 ∥p−ph∥0,Ω≤¯C2inf(vh,θh,qh)∈Hh×Zh×Qh(∥u−vh∥H+∥ω−θh∥Z+∥p−qh∥0,Ω). (3.10)

Proof. The result follows using the same arguments considered in the proof of Corollary 1, but using the discrete inf-sup condition (3.6). We omit further details.

### A priori error estimates.

Let us introduce for a given , the Nédeléc global interpolation operator . From [1] we known that for all with , there exists independent of , such that

 ∥θ−Rhθ∥Z≤Chmin{s,k+1}∥θ∥Hs(curl;Ω). (3.11)

On the other hand, for the Raviart-Thomas interpolation , with , we recall (see e.g. [24]) that there exists , independent of , such that for all :

 ∥v−Πhv∥H≤Chmin{s,k+1}∥v∥Hs(div;Ω)∀v∈Hs(div;Ω)∩H. (3.12)

Finally we recall that the orthogonal projection from onto the finite element subspace , here denoted , satisfies the following error estimate for all :

 ∥q−Phq∥0,Ω≤Chmin{s,k+1}∥q∥s,Ω∀q∈Hs(Ω). (3.13)

These operators fulfil the following commuting diagram

 divRhv=Ph(divv)∀v∈Hs(Ω)3∩H(div;Ω).

The following result summarises the error analysis for our mixed finite element scheme (3.4).

###### Theorem 3.2

Let be an integer and let and be given by (3.1), (3.2), and (3.3). Let and be the unique solutions to the continuous and discrete problems (2.3) and (3.4), respectively. Assume that , , and , for some . Then, there exists independent of such that

 ∥u−uh∥H+∥ω−ωh∥Z+∥p−ph∥0,Ω≤^Chmin{s,k+1}(∥u∥Hs(div;Ω)+∥ω∥Hs(curl;Ω)+∥p∥s,Ω).

Proof. The proof follows from (3.9), (3.10), and standard interpolation estimates satisfied by the operators , and (see (3.11), (3.12) and (3.13), respectively).

## 4 Discontinuous Galerkin method

In this section, we propose and analyse a DG method for (2.2). We provide solvability and stability of the discrete scheme by introducing suitable numerical fluxes. A priori error estimates are also derived.

### Preliminaries.

Apart from the definitions laid out at the beginning of Section 3, let us denote by the set of internal faces, by the set of external faces on and by the set of external faces on . We set . We denote by the diameter of each face . Let and be two adjacent elements of and let (respectively ) be the outward unit normal vector on (respectively ). For a vector field , we denote by the trace of from the interior of . We define jumps

 [[v]]T:=v+×n++v−×n−,[[v]]N:=v+⋅n++v−⋅n−,[[q]]:=q+n++q−n−,

and averages

 {{v}}:=12(v++v−),{{q}}:=12(q++q−),

and adopt the convention that for boundary faces , we set , , , and .

Suitable finite dimensional spaces for vorticity and velocity that remove the restriction of continuity are defined by:

 ~Zh ~Hh :={vh∈L2(Ω)3:vh|T∈Pk+1(T)3∀T∈Th},

and we remark that the space for pressure approximation will coincide with the one used in Section 3, that is .

### Discrete formulation and solvability analysis.

Multiplying each equation in (2.2) by suitable functions, the resulting DG scheme consists in finding , such that for any test functions and for all elements in the partition

 σ∫Tuh⋅vhdx+√ν∫Tωh⋅curlvhdx+√ν∫∂Tˆωh⋅(vh×n)ds+1√ν∫T(ωh×β)⋅vhdx (4.3) − ∫Tphdivvhdx+∫∂Tˆphvh⋅nds=∫Tf⋅vhdx, ∫Tωh⋅θhdx=√ν∫Tuh⋅curlθhdx+√ν∫∂Tˆuωh⋅(θh×n)ds, −∫Tuh⋅∇qhdx+∫∂Tˆuph⋅nqds=0,

where and are numerical fluxes, which approximate the traces of , and on the boundary. The fluxes and are related to the - operator and are defined by

 ˆωh:=⎧⎪⎨⎪⎩{{ωh}}+C11[