Mixed and Hybrid Finite Element Methods for Convection-Di ff usion Problems and Their Relationships with Finite Volume : The Multi-Dimensional Case

We introduced in (Fortin & Serghini Mounim, 2005) a new method which allows us to extend the connection between the finite volume and dual mixed hybrid (DMH) methods to advection-diffusion problems in the one-dimensional case. In the present work we propose to extend the results of (Fortin & Serghini Mounim, 2005) to multidimensional hyperbolic and parabolic problems. The numerical approximation is achieved using the Raviart-Thomas (Raviart & Thomas, 1977) finite elements of lowest degree on triangular or rectangular partitions. We show the link with numerous finite volume schemes by use of appropriate numerical integrations. This will permit a better understanding of these finite volume schemes and the large number of DMH results available could carry out their analysis in a unified fashion. Furthermore, a stabilized method is proposed. We end with some discussion on possible extensions of our schemes.


Introduction
The aim of this paper is to extend to multidimensional hyperbolic and parabolic problems the results of (Fortin & Serghini Mounim, 2005).In (Fortin & Serghini Mounim, 2005) we presented in a mixed finite element framework, two nonstandard formulations in one space dimension allowing, among other things, to recover in a unique fashion, numerous finite volume schemes.We shall present in this paper a two-dimensional extension.To achieve the discretization we shall use mixed finite elements to approximate the convective and diffusive fluxes.We concentrate on the case of the lowest order Raviart-Thomas spaces (Raviart & Thomas, 1977), though many of our results will be more general.In particular, to obtain higher-order accuracy finite volume schemes less sensitive to partitions of domain, we can follow the same procedure but utilize others usual mixed finite element approximation subspaces of H(div, Ω) such as BDM spaces (Brezzi, Douglas & Marini, 1986).This paper is organized as follows.In the next section, we systematically build the extension of the formulations of (Fortin & Serghini Mounim, 2005) to the multi-dimensional case.Results of existence and uniqueness are presented.In Section 3, we build the numerical schemes on triangular or rectangular partitions of the domain Ω using the Raviart-Thomas finite elements of lowest degree (Rarviart & Thomas, 1977) to approximate the fluxes.In Section 4, we show the link between the finite volume methods and our formulations using the suitable quadrature formulas suggested in (Baranger, Maitre, & Oudin, 1996) (triangular case) or trapezoidal quadrature (rectangular case) to diagonalize the element mass matrix and numerical integration formulas to approximate the term ∫ K f (u).q dx (here f (u) is the convective flux).Moreover, we focus our attention to the multidimensional extension of some recent schemes.Finally, in Section 5 we discuss how we can also establish the extension of the second method to both convection-diffusion equations and system of equations.We end with some discussion on possible extensions of our schemes.The DMH formulation requires the solution of a linear system; the approach adopted here was suggested by Arnold and Brezzi (Arnold & Brezzi, 1985).In their method, one eliminates the velocity and the flux unknowns to leave a system for the Lagrange multipliers alone.The lowest order Raviart-Thomas spaces have one Lagrange multiplier unknown per edge.

Preliminaries
In this section, we start by presenting a few properties relative to the decomposition of Ω .The proofs can be found in Thomas (Thomas, 1977) or Brezzi-Fortin (Brezzi & Fortin, 1991).Let us consider the following partial differential equation for u: with given initial data u(x, 0) = u 0 (x) in Ω, and corresponding suitable boundary conditions on Γ = ∂Ω.Here the velocity field and the diffusion coefficient (which is positive) are denoted by a(x, t) and ν respectively.Moreover, for simplicity and without loss of generality, we restrict our attention to the case where div(a) = 0; taking into account a.n = 0 on Γ×]0, T [, and assuming that a(., t) ∈ L ∞ (Ω); also, that u = 0 on Γ×]0, T [.Let then T h = {K i } NT i=1 be the usual non-overlapping finite element triangulation of the domain Ω = ∪ K∈T h K.We also consider for 1 ≤ t < ∞ the spaces and which are both Banach spaces endowed with the usual norms and The following proposition, showing that the space W t (div; Ω) can be characterized as a subspace of X t (Ω), is well known: for all v ∈ W 1,s 0 ( Ω ), where t is the conjugate of s.Finally, let us recall that for all µ ∈ W 1/t,s (Γ) and µ * ∈ W −1/t,t (Γ) we have and Here ⟨., .⟩denotes the duality pairing between W −1/t,t (Γ) and W 1/t,s (Γ).Now, we are able to formulate precisely the parabolic problem (1).

Dual Mixed and Hybrid Finite Element Method DMH1
In this section, we deal with the first variational formulation of the convection-diffusion problem previously considered.The principal objective of the following is to extend to the 2D case the first formulation proposed in (Fortin & Serghini Mounim, 2005).
As in (Fortin & Serghini Mounim, 2005), we shall begin from equation (1).Next, we introduce convective and diffusion fluxes as auxiliary variables to obtain the following first order system: Boundary conditions and initial data are added to these equations.
Here the boundary of Ω is supposed polygonal.Introducing λ, a Lagrange multiplier to enforce the continuity constraint of q.n at elements' interfaces, we obtain We shall evidently have to make a proper choice for the space of λ.This technique was first used by Fraejis DE Veubeke (Fraejis DE Veubeke, 1965;Fraejis DE Veubeke, 1975).In the same way, we shall impose the continuity of the convective flux p at the interfaces of elements K i , with the aid of a Lagrange multiplier noted λ.Hence, the last equation of the system (9) can now be rewritten as: Now, to make the variational formulation precise we have to define the functional spaces.Let us consider an exponent 4 3 < s < 2 and its conjugate t.In addition, we set and M 1 = L t (Ω).Hence, the system (9) could be written as follows: The last two equations of system (12) insure the continuity of p.n and p.n at the interfaces of T h .

Remark 2
The proof of the existence of Lagrange multipliers λ and λ, which insure the continuity of normal traces is similar to the one given in Brezzi-Fortin (Brezzi & Fortin, 1991).
Theorem 1 Let u be the solution of the convection-diffusion problem (1).If p = grad u and p = au verify then ((p, p, λ); (u, u| ∂K )) is the unique solution of (12).

Dual Mixed and Hybrid Finite Element DMH2
Let us consider the following partial differential equation for u: where a(x, t) is the velocity.
The dual mixed and hybrid formulation of the problem ( 14) can be achieved proceeding in a similar way as in (Fortin & Serghini Mounim, 2005).First, we introduce again the convective flux as an auxiliary variable to obtain: Boundary conditions and initial data are added to these equations.
Next, we assume that the solution has discontinuities at each boundary of every K i .Finally, we suggest that the "Rankine-Hugoniot" jump condition should be imposed on the numerical flux through the elements' interfaces.With this aim, we employ a Lagrange multiplier to relax these jump conditions at the interelements. and We denote by λ the Lagrange multiplier and by a i ≥ 0 some approximation of the local propagation speeds.As we will see, it plays a stabilization parameter role, and for an appropriate choice of a i , we can recover in particular the well-known upwind numerical fluxes (approximate Riemann solvers).Here, we have set [ p.n] = pi .n− pe .n,where "the interior normal trace " pi .n(resp."the exterior normal trace" pe .n) is defined as the normal trace of the restriction p| K (resp.p| (Ω\K) ) on ∂K, and with [u] = u i − u e stands for the jump of u, where "the interior trace" u i and "the exterior trace" u e of u are defined from the trace of u as above.In order to obtain the non-standard mixed hybrid finite element formulation, we keep the same notation for functional spaces as in the previous section.Thus, the variational formulation of ( 14) is to find (u(t), λ(t)) ∈ M 1 × M 2 and p(t) ∈ X such that: The last equation of ( 18) imposes jump conditions to the flux p.n at the interelement interfaces of T h .We are now able to obtain a result of existence and uniqueness of the solution for problem (18), using an argument of regularity as in Theorem 2.23.

Spatial Approximation
Let us now deal with the finite element approximation of the above dual variational problem.To this aim, we assume that Ω is a polygonal bounded domain of I R. Let T h be a standard finite element partition of Ω.

Spatial Approximation of DMH1
To define a discretization of problem ( 12), we first define the approximation spaces X h , M 1h and M 2h of X, M 1 and M 2 , respectively.We thus build an appropriate approximation space for q ∈ X in the standard way as follows: where RT 0 (K) is the local Raviart-Thomas space of lowest degree.
To obtain the approximation of the spaces M 1 and M 2 , we define and where R 0 (∂K) = { φ; φ| e ∈ I P 0 (e) ∀ e ∈ ∂K } .
Finally, the discrete version of ( 12) is given by: To these equations we have to add given initial data, and corresponding suitable boundary conditions.
Theorem 3 Problem ( 23) has a unique solution.
Proof.Solving the problem (23) leads to solve a square finite-dimensional system of linear algebraic equations.Hence, we only have to prove the uniqueness of the solution.
If we consider the homogeneous problem associated to (23), f=0, and if the boundary conditions and initial value vanish, the proof is similar to that of the Theorem 1, in regard to u h , p h and ph .For the uniqueness of λ h and λh , we consider and we then deduce that λ h = 0, and similarly we have λh = 0.

Spatial Approximation of DMH2
To define the associated discrete problem of (18), we may now use the same suitable subspaces as above, and we also proceed in the same manner.Hence, we obtain the following discrete problem: a i is again the approximation of the local velocity at ∂K i , (a i ≥ 0).In a similar way as we did for Theorem 3, we can prove the following result: Theorem 4 The problem ( 24) has a unique solution.

Analogy with the Finite Volume Method
The contents of this section are discussed in more details in (Serghini Mounim, 2000).Notice that the most attractive property of finite volume methods is that the conservation of quantities such as mass, momentum and energy is exactly satisfied on each finite volume and also on the whole domain.To make the link explicit with finite volume methods, we use the RT vector finite elements to approximate the diffusion and convective fluxes.We point out this relationship, using a suitable quadrature formula to approximate the terms .q dx, f (u) denotes here the convective flux.

Rectangular Case
We consider a quadrangulation T h of Ω, and we denote by K the generic rectangle of T h defined by [x K , x K + ∆x K ] × [y K , y K + ∆y K ], and the degrees of freedom introduced by Raviart-Thomas.The usual RT basis functions of K are defined by: where |K| is the area of K.
In order to diagonalize the mass matrix ∫ K p.q dx, still retaining the good numerical properties of the mixed approximation, we employ the trapezoidal quadrature rule to obtain the following diagonal matrix: For simplicity, we suppose f = 0, and in what follows we shall write u, p, p instead u h , p h , ph , ... and omit the dependance of index sets u, p, p, ... on h.Now, let us fix our attention on two adjacent rectangles K and K i , i ∈ {S , E, N, W}, and denote the values of u in K and K i by u K and u i respectively.Next, as in the 1D case (see (Fortin & Serghini Mounim, 2005)) a direct computation provides: from which we can deduce the approximation of the convective flux at the faces of where the outward normal to the faces e i of K is denoted by n i , with |e i | stands for the length of e i for i ∈ {S , E, N, W}.
Then, we observe that the expressions of αi are equivalent to a finite volume method.Indeed, we obtain a "central difference scheme" (CDS).This scheme is not spatially stable; it does not prevent oscillations near discontinuities and at shocks, (for the properties of CDS see Lazarov et al (Lazarov, & Vassilevski, 1996)).

Other Formulation and Upwind Schemes
If we impose the continuity of convective flux only at inflow boundary of each control volume K, we can obtain an upwind scheme.For this purpose, we utilize the following form: ∂K − denotes the inflow boundary of K, which is defined by Using the same techniques as in the previous subsection, we can derive easily the numerical convective flux at the interface between adjacent cells: where u represents some average state between u i and u K .
For the linear convection case ( f (u) = au where a denotes again the convective field), the above expression then reads: where i ∈ { S , E, N, W } , and This scheme coincides with the standard UDS (upwind difference scheme) see for example (LeVeque, 2002).Before considering the generalization of the previous schemes, we note that one advantage of the use of λ is that no (approximate) Riemann solvers are required.

Remark 3
Since u is discretized using piecewise polynomial functions and our formulation involves the boundary integrals, we can introduce an upwinding in the formulation.To achieve the generalization of the schemes above, we propose the following form to replace the form already used n ds in the system (23).The following notation is used: Taking the same approach as above, this results in the generalized upwind scheme: ) .
Here f = ( f 1 , f 2 ), and β is the dissipation parameter (0 ≤ β ≤ 1) which controls the combination of fully upwind and centered scheme.We have supposed without loss of generality that the inflow boundary of K is e S ∪ e W . Notice that for β = 1, we recover the scheme ( 29), (30).

Triangular Case
We consider now the case where T h is a standard triangulation of Ω (here, we suppose again that Ω is a polygonal bounded domain of I R).We begin with a brief description of the notation used here.For each triangle K, we have: • a r : r th vertex of K with coordinates (x r , y r ); • e r : face opposite to vertex a r , with length | e r |; • h r : the length of perpendicular dropped from the vertex a r ; • n r : unit exterior normal to face e r ; • θ r : angle at vertex a r ; for r ∈ {i, j, k}.
We also consider the local shape functions {p r } r∈{i, j,k} of RT 0 (K) for every triangle K.The degrees of freedom are the fluxes through the faces {e r } r∈{i, j,k} , and are defined by: from which, for every local vectorial function p r of RT 0 (K) the well known properties follow: where δ rs denotes the Kronecker symbol We introduce the bilinear form a(p, q) ≡ ∫ K p.q dx for all p, q in RT 0 (K).Then, a(p r , p s ) is a 3 × 3 symmetric, positive definite matrix, with c r ≡ cot θ r , for r, s ∈ {i, j, k}.We employ the quadrature formula proposed in (Barranger, Maitre, & Oudin, 1996), to diagonalize the local mass matrix (34).We obtain the matrix: From this, after temporal approximation, and considering the neighbor elements of K which appear in the formula, the equation expressing the discrete conservation law becomes: and thanks to the following relation, where we have set β r = 1 2 cotθ r , with r = i, j, k.From the above equation, in conjunction with the continuity of the normal trace of the diffusive flux, we deduce Using ( 37) and (38), the numerical diffusion flux across edge e r reads We note that β r = d r |e r | , r = i, j, k (see the annexe in (Bank & Rose, 1987;Kerkhoven-Jerome, 1990)), where d r denotes the distance between the circumscribed circle of K and the center m r on sides e r (see previous and below figures).We mention that β r is also given by , where we have set s = I(r), s ′ = I(s), with I defined as and I(k) = i, and t r stands for the unit tangent to e r for r = i, j, k.
With the aid of the new expression of β r the system (39) becomes where d r + d r ′ is the distance between the circumcenters C, C r of two adjoining triangles K and K r .The above distance does not vanish, in the case of a strictly acute triangle.
We now deal with the convective flux using the same notation and procedure as above.We thus start from the equation We approximate the mass matrix with the same quadrature formula as in the diffusive case, in addition, we use the following numerical integration which is exact for polynomials of degree 0: where H K is the orthocenter of the triangle K which is considered acute, and |K| denotes the area of K as illustrated in Figure 3.
Thanks to the above formula (41), the integral Finally, these approximate relations allow us to write the numerical convective flux as It is precisely the central finite volume scheme "CDS".

Remark 4
In view of the linear case f (u) = au, we can define a as a linear combination of the local shape functions p r of RT 0 (K), after some manipulation, we can exhibit the numerical flux through the face e r where For other approximations the reader can refer to (Serghini Mounim, 2000).
Remark 5 If we utilize the bilinear form (31) and the quadrature formula (41), the numerical convective flux (43) becomes for r ∈ {i, j, k}, with Ξ = min ( sign(a.nr ), 0 ) and β being again the upwinding parameter (0 ≤ β ≤ 1).Note that for β = 1 (or the continuity constraint imposed to p only at the inflow boundary of K ), we get the following upwind finite volume scheme: .n r , r = i, j, k, for the stability and error estimates of CDS and UDS on Voronoi meshes we refer to (Mishev, 1998).

DMH2 and Scheme2
The results of this section are derived in a similar fashion to DMH1 case.

Rectangular Case
Consider again the case of a rectangular mesh T h of Ω, and the local degrees of freedom associated to RT 0 (K).Here we denote by n e r (resp.n i r ) the exterior unit normal vector, i.e., oriented from K to K r (resp.interior unit normal vector, i.e., oriented from K r to K) along each edge e r = ∂K ∩ ∂K r ; clearly, the property n e r = −n i r holds.Now, following the same lines as in previous sections, we can state the result: Proposition 2 If we use the trapezoidal quadrature formula to diagonalize the element mass matrix, and exact computation of the spatial integral ∫ K f (u).q dx; then, the exterior αe r and the interior αi r fluxes approximating ∫ e r pe .ne r ds and ∫ e r pi .ne r ds respectively across edge e r are given by: .n e r + sign(a.ne r ) a r | e r (u r − u K ) ) , .n e r + sign(a.ne r ) a r | e r (u K − u r ) ) where λr is given by: where λr is given by: for r = i, j, k.
If we restrain to equilateral triangles, one then gets obviously ∥ and, in particular, the approximation of the exterior flux becomes: for r = i, j, k.

Remark 7
In view of the linear case f (u) = au, we can express a as linear combination of the local shape functions p r of RT 0 (K).After some manipulations, we can exhibit the exterior numerical flux through the face e r : for r = i, j, k, and the approximation a h of a is defined by

Remark 8
To simplify, let K be an equilateral triangle with centroid G, and let K r , r = i, j, k be the neighbouring equilateral triangles, with centroids G r (r = i, j, k).If we assume u| K linear, and the least-squares gradient gradu(t)| K = (a K (t), b K (t)) for the triangle K is then chosen in order to minimize the functional Of course, the minimum is obtained when for details and the expressions of a K (t) and b K (t) see (Champier, 1992).Next, we define over all K and the exterior and interior traces of u at the face e r are now given by: , where x m r denote the coordinates of the center on sides e r (see Fig. 2).The resulting exterior flux is for r = i, j, k.The cell average flux f (u K ), (also f (u K r ), r = i, j, k) can be approximated by a second-order quadrature in time (such as the mid-point rule).The resulting scheme can be viewed as a natural new version of multidimensional extension of the second order central NT scheme (Nessyahu & Tadmor, 1990).For other extensions of the NT scheme, see (Arminjon & Viallon, 1995;Jiang & Tadmor, 1998).
We close by noting that we can get a 2D extension of the family of semi-discrete central schemes of Kurganov and Tadmor (Kurganov & Tadmor, 2000).Our recipe is similar to that above, taking into account the quadrature formula given below: From this, we can express the exterior numerical flux in the form: for r = i, j, k, we then obtain the extension to the triangular grid case of the semi-discrete KT scheme (Kurganov & Tadmor, 2000), and where a r (t) are the maximal local speeds.For more see (Kurganov & Tadmor, 2000).

Discussion
1-The extension to a system of equations can be carried out in the same fashion as (Fortin & Serghini Mounim, 2005), and again with the aid of Roe's construction (ROE, 1981).
2-The generalization to a convection-diffusion equation can also be performed using the extension to 2D of the variational formulation used in (Fortin & Serghini Mounim, 2005).We use, for this, the fact that the solutions to the general initial value problem of the viscous conservation laws converge, in the zero dissipation limit as ν goes zero, to the solutions of the hyperbolic conservation laws.Next, the "RH" jump conditions at the interelements are imposed again to the convective flux.Finally, we can establish the existence and uniqueness of the solution in this case taking the same approach to previous formulations, and using in particular an argument of regularity.
3-Note that the simplicity of the method makes it easy to extend to the three-dimensional case, in addition, if the triangulation of Ω is made of tetrahedrons, and the mesh is restricted to regular subdivisions of the domain, the quadrature formula given in Baranger et al (Baranger, Maitre, & Oudin, 1996) can be used to diagonalize the mass matrix.Moreover, the result of (Baranger, Maitre, & Oudin, 1996) cannot be extended to the general 3D case, for more details see (Thomas & Trujillo, 1997).
4-Notice that the methods presented here can also be viewed as an extension of the mortar finite elements method to transport problems.On the other hand, we mention that there is a possible link between our methods and the mimetic finite difference methods (see for example (Hyman & Shashkov, 1997)) that could be made by choosing an appropriate inner product.
Finally, the stability and the feasibility of the method had been successfully tested on numerical examples.For this, the dual mixed finite element method associated to the first-order accurate implicit Euler (for the time-discretization) has been experimented on structured and unstructured meshes see (Serghini Mounim, 2000).