# L2-stability independent of diffusion for a finite element
– finite volume discretization of a linear
convection-diffusion equation

Paul Deuring, Robert Eymard, Marcus Mildner

► To cite this version:

Paul Deuring, Robert Eymard, Marcus Mildner. L2-stability independent of diffusion for a finite element – finite volume discretization of a linear convection-diffusion equation. SIAM Journal on Numerical Analysis, 2015, 10.1137/140961146. hal-01060065v2

HAL Id: hal-01060065

https://hal.science/hal-01060065v2

Submitted on 13 Nov 2014

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
L²-stability independent of diffusion for a finite element – finite volume discretization of a linear convection-diffusion equation.

Paul Deuring,* Robert Eymard†and Marcus Mildner*

November 9, 2014

Abstract

We consider a time-dependent and a steady linear convection-diffusion equation. These equations are approximately solved by a combined finite element – finite volume method: the diffusion term is discretized by Crouzeix-Raviart piecewise linear finite elements on a triangular grid, and the convection term by upwind barycentric finite volumes. In the un- steady case, the implicit Euler method is used as time discretization. This scheme is shown to be unconditionally L²-stable, uniformly with respect to the diffusion coefficient.

**AMS subject classifications.** 65M30, 65M60, 76M10, 76M12.

**Key words.** convection-diffusion equation, combined finite element - finite volume method, Crouzeix-Raviart finite elements, barycentric finite volumes, upwind method, stability.

# 1. Introduction

We consider the convection-diffusion equation

$$ \partial_t u - \nu \Delta u + \beta \cdot \nabla u = g \quad \text{in } \Omega \times (0, T), \qquad (1.1) $$

supplemented by the initial and boundary conditions

$$ u(x, 0) = u^{(0)}(x) \quad \text{for } x \in \Omega, \qquad u|_{\partial \Omega} \times (0, T) = 0, \qquad (1.2) $$

respectively. Here $\Omega \subset \mathbb{R}^2$ is a bounded open polygon with Lipschitz boundary, $\nu$ and $T$ are positive reals, and $\beta: \Omega \to \mathbb{R}^2$, $u^{(0)}: \Omega \to \mathbb{R}$ as well as $g: \Omega \times (0, T) \to \mathbb{R}$ are given functions. Our key assumptions concern the advective velocity $\beta$: we require that $\beta \in H^1(\Omega)^2$,

$$ \nabla \cdot \beta = 0, \quad -\beta \cdot \nabla \varphi \ge \underline{\beta} \quad \text{in } \Omega \qquad (1.3) $$

for some function $\varphi \in C^1(\bar{\Omega})$ and some constant $\underline{\beta} > 0$. In the case where $\beta(x) = \beta_0$ for all $x \in \bar{\Omega}$ and for some $\beta_0 \in \mathbb{R}^2 \setminus \{0\}$, a suitable function $\varphi$ is given by $\varphi(x) := -\beta_0 \cdot x$.

<footnote>*Univ Lille Nord de France, 59000 Lille, France; ULCO, LMPA, 62228 Calais, France.</footnote>

<footnote>†Université Paris-Est Marne-la-Vallée, 5 boulevard Descartes, 77454 Marne-la-Vallée, France.</footnote>
We also consider the steady variant of problem (1.1), (1.2), that is,

$$-\nu \Delta U + \beta \cdot \nabla U = G \quad \text{in } \Omega, \qquad (1.4)$$

$$U|_{\partial \Omega} = 0, \qquad (1.5)$$

where $G: \Omega \to \mathbb{R}$ is another given function. These problems are of particular interest in the convection-dominated regime, that is, if $\nu \ll |\beta|$, an interest that seems to be motivated by the belief that the preceding problems in the convection-dominated case show some affinity (although distant) with the Navier-Stokes system in the same regime. In this spirit, numerical schemes working well for that latter system are sometimes reduced to problem (1.1), (1.2) or (1.4), (1.5) so that they may be accessible to theoretical studies regarding stability or accuracy.

In the work at hand, we consider a discretization of (1.1), (1.2) and (1.4), (1.5), respectively, that is motivated in this way. This scheme may be described as follows: the diffusion term in (1.1) and (1.4) is discretized by piecewise linear Crouzeix-Raviart finite elements, and the convective term by an upwind finite volume method based on barycentric finite volumes on a triangular grid. Choosing an explicit time discretization, Feistauer e.a. [15, Section 7], [23, Chapter 4.4] tested this FE-FV method in the case of high-speed compressible Navier-Stokes flows in complex geometries and obtained very satisfactory results.

In [12], we applied this FE-FV method to problem (1.1), (1.2), using the implicit Euler method as time discretization. Under the assumption that $\beta$ is constant, we showed that the approximate solution provided by this approach may be estimated in the $L^2$-norm against the data, with the constant in this estimate being independent of the diffusion parameter $\nu$. An analogous result was established with respect to problem (1.4), (1.5). However, we had to impose a condition on the grid ([12, (3.9)]) that goes beyond standard assumptions in a finite element or finite volume context, and which is not satisfied by certain structured regular meshes. Additional details in this respect may be found in [34].

It is the aim of the work at hand to establish the stability estimates from [12] without this inconvenient condition. In fact, we will show that even if the grid is only required to be shape-regular (minimum angle condition), an upper bound independent of $\nu$ and only involving the data may still be constructed for the $L^2$-norm of the approximate solution of (1.1), (1.2) and (1.4), (1.5), respectively, obtained by the FE-FV approach described above, with the implicit Euler method being used as time discretization in the unsteady case. As an additional generalization with respect to [12], the function $\beta$ need no longer be constant. Instead it will only be supposed to satisfy (1.3). Our stability constants depend polynomially on $\beta^{-1}$, $\|\beta\|_{1,2}$ and an upper bound of $\varphi$ and $\nabla\varphi$, respectively. For a detailed statement of our results, we refer to Theorem 2.1 and 2.2 below.

We remark that our approach also carries through if the functions $\beta$ and $\varphi$ in (1.1) depend on time, provided that the spatial $H^1$-norm of $\beta$ is bounded uniformly with respect to time, and (1.3) holds at any time with the same constant $\beta$. Moreover the diffusion term $\nu\Delta u$ may be replaced by an elliptic operator in divergence form $\nabla \cdot (A \cdot \nabla u)$, where $A = A(x)$ or $A = A(x, t)$ is a symmetric matrix in $\mathbb{R}^{2 \times 2}$ which is positive definite uniformly with respect to $x$ and (in the unsteady case) $t$. However, since our presentation would become more complicated while our proofs would remain without essential modification, we will not make these generalizations explicit.
Our theory should be expected to hold in the 3D case as well. Poincaré's inequality (2.5) is available in that case as well ([31], [5], [45]), so there should be no obstacle in this respect. There is only one point that might need some effort, that is, to prove an analogue or find a replacement of equation (2.6) pertaining to the discrete $L^2$-scalar product.

Condition $(1.3)_1$ on $\beta$ cannot be expected to be sufficient for the validity of our stability results, even if $\beta \in C^1(\bar{\Omega})^2$. This will be explained in Section 5. So inevitably an additional assumption has to be imposed on $\beta$. The one we have chosen – inequality $(1.3)_2$ – may look ad hoc. However, in Section 5 we will present a strong heuristic argument implying that $(1.3)_2$ is necessary for Theorem 2.1 and 2.2 to hold. This observation means in particular that our result is not valid in the purely diffusive case ($\beta = 0$), that is, for Poisson's equation $\nu \Delta U = G$ in the steady and the heat equation in the evolutionary case. This is borne out by the counterexamples we present in Section 5.

Assumption $(1.3)_2$ means that the advective velocity $\beta$ exhibits neither closed curves nor stationary points (points $x \in \Omega$ with $\beta(x) = 0$). In fact, it is shown in [13] that if $\beta$ is smooth and presents these geometrical properties, there is a function $\varphi$ with $(1.3)_2$. This result was generalized to the case $\beta \in W^{1,\infty}(\Omega)^2$ in [3], a reference we will come back to below.

The key argument allowing us to improve the theory from [12] consists of a weak BV estimate providing an upper bound for a certain weighted variation of our approximate solution. This upper bound involves the $L^2$-norm of the solution and of the data. For technical reasons, the estimate in question will not appear explicitly, but will be split into two parts (Lemma 3.3 and (4.1) in the unsteady case, Lemma 3.3 and (4.11) in the steady one). Usually this type of inequality is exploited to pass to the limit in numerical schemes for the approximation of scalar hyperbolic equations; see [42] or [9] for example. Here, however, it is a critical part of the sequence of inequalities leading to our stability estimate.

A very large body of work deals with convection-diffusion-reaction equations, that is, with equations of the form

$$-\nu \Delta U + (\beta \cdot \nabla)U + \mu U = G \quad \text{in } \Omega \qquad (1.6)$$

and their time-dependent counterparts, under the assumption that

$$\mu - \nabla \cdot \beta/2 \ge \gamma_0 \quad \text{for some } \gamma_0 > 0. \qquad (1.7)$$

As is well known, in this situation various discretizations allow error and stability estimates with constants independent of $\nu$. We refer to the monographs [37, Chapter 8 and 12], [16, Section 5.2.3, 5.4.4] and [39] for more details. In the case of (1.1) and (1.4), however, the condition $\mu - \nabla \cdot \beta/2 \ge \gamma_0$ reduces to the inequality $-\nabla \cdot \beta \ge \gamma_0$, which, in view of (1.3), is not satisfied.

Sometimes, when equations (1.1), (1.4) or (1.6) are solved numerically, anisotropic meshes ("Shishkin meshes"; see [39]) are used. They take account of boundary layers, which, in turn, depend on $\nu$. Construction of such meshes involves the advective velocity $\beta$ and becomes difficult if $\Omega$ does not exhibit a simple geometry. Moreover, for this type of mesh, the constant $\sigma_0$ in condition (2.1) relative to the grid shape is influenced by the diffusion coefficient. In view of these features, we think Shishkin's approach works best in situations that are not the main interest of the work at hand. So, when we compare our theory with what is available in existing
literature, we look for stability estimates for discretizations based on conventional grids and pertaining to scalar convection-diffusion equations whose structure is similar to that of (1.1) or (1.4) and which do not satisfy (1.7).

In this perspective, our work is similar in spirit to the theory developed by Ayuso, Marini in [3]. These authors consider problem (1.6) with nonhomogeneous mixed boundary conditions, applying various discontinuous Galerkin methods in order to obtain numerical solutions. Starting from $(1.3)_2$ and the inequality $\mu - \nabla \cdot \beta/2 \ge 0$, they establish not only $\nu$-uniform stability estimates, but also error estimates with constants independent of $\nu$. Still we think the theory in [3] is not completely satisfying. Besides the fact that the DG methods in question involve a parameter that has to be chosen appropriately, a strange technical condition ([3, (H2)]) is imposed on $\beta$ in addition to $(1.3)_2$. Moreover the constants in stability and error estimates depend exponentially on $\max \varphi - \min \varphi$ ([3, Lemma 4.1]). Apart from [3], we do not know of any reference dealing with stability or error estimates of numerical solutions to (1.1), (1.4) or (1.6) on the basis of condition $(1.3)_2$.

Combined FE-FV methods similar to the one considered in the work at hand were already studied in a number of earlier articles besides reference [12] already mentioned. Ohmori and Ushijima [36] treated the steady problem (1.4), (1.5) under assumptions on $\beta$ that are weaker than ours. They discretized (1.4) by Crouzeix-Raviart finite elements, using an upwind version of this element in order to approximate the convective term. For constant $\beta$, this approach reduces to the FE-FV discretization of (1.4), (1.5) we study in the work at hand. The authors provided a $\nu$-uniform $L^\infty$-bound for their discrete solutions, under the assumption that the right-hand side $G$ in (1.4) is negative and the underlying triangulation is of acute type. Also an estimate of the $H^1$-error is derived. The constant in this estimate depends on $\nu$, which is inevitable for this type of error. Feistauer e.a. [1], [15] considered a scalar time-dependent nonlinear conservation law with a diffusion term. Discretizing this equation by the combined FE-FV scheme described above, with a rather general numerical flux adapted to the nonlinearity in question, and with a semi-implicit Euler method as time discretization, they established $L^\infty$-stability independent of $\nu$. However, they required their triangulations to be of weakly acute type, and they supposed that the ratio of the time step to the grid size is small with respect to a quantity depending on the lifetime $T$ of the solution, on the $L^\infty$-norm of the initial data and the $L^\infty$-norm of the right-hand side in the differential equation. In [1] the stability result in question is a key element in the proof of a convergence result, whereas in [15], it is applied for deriving $L^2(L^2)$- and $L^2(H^1)$-error estimates, which involve constants with exponential dependence on $1/\nu$. The articles [21], [22] and [25] establish results analogous to those in [1] and [15], but for a combined finite element - finite volume method involving piecewise linear conforming finite elements and dual finite volumes. Similar $L^2(L^2)$- and $L^2(H^1)$-error estimates as in [15], [22] and [25] are shown in [24], but pertaining to various discontinuous Galerkin discretizations.

Reference [18] deals with a degenerate time-dependent convection-diffusion equation which is much more general than (1.1) in some respects. A discretization is performed with a similar finite element - finite volume method as in the present article. The theory in [18] yields $\nu$-uniform $L^\infty$-bounds for the discrete solution, but covers (1.1) only in the case $g=0$. This means in particular the approach in [18] does not provide an access to the steady problem (1.4). In [19], a convergence proof for a mixed FE-FV method is carried through for a similar problem as in [18]. The method in question involves conforming P1 finite elements and cell-centered
finite volumes.

The stability estimates in [36], [1], [15], [21], [22], [25], [18] require rather severe assumptions because they are deduced from the maximum principle. In fact, even in the case of the continuous steady equation (1.4), this principle introduces a factor of the form $e^{\text{diam} \Omega \|\beta\|_{\infty}/\nu}$ if it is applied in order to yield $L^{\infty}$-bounds; see [28, Theorem 3.7]. The same factor arises for the parabolic problem (1.1) if $g \neq 0$ ([26, (3.3.10)]). So the assumptions in the references just mentioned are needed in order to circumvent this phenomenon.

In [11] we derived $L^2$- and $H^1$-error estimates for the same discretizations of (1.1), (1.2) and (1.4), (1.5) as considered here, under the assumption $\beta \in \mathbb{R}^2 \setminus \{0\}$. The constants in our estimates are analyzed with respect to their dependence on $\nu$. This analysis, which includes the norms of the exact solution appearing in our error bounds, reveals a polynomial dependence on $\nu^{-1/2}$.

$L^2$-stability and error bounds independent of $\nu$ are constructed in [38] for finite volume approximate solutions to (1.1). Article [38] is the most recent one in a series of papers dealing with $\nu$-independent stability and error estimates, with previous publications treating Euler-Lagrangian and finite element methods. The corresponding references may be found in [38]. However, although constants independent of $\nu$ constitute the best possible case one might hope for, the type of theory developed in [38] and its predecessor papers also has drawbacks. It requires periodic boundary conditions instead of Dirichlet ones, it assumes $\Omega$ to be a rectangle and the grid to be structured according to this geometry, and it relies on Gronwall's inequality in a way which leads to exponential dependence of stability and error bounds on the $W^{1,\infty}$-norm of the convective velocity. Moreover, using Gronwall's inequality leaves open the question of how to deal with the steady case.

In [2], a mixed finite element method is applied to an unsteady nonlinear degenerate convection-diffusion equation. An argument based on Gronwall's inequality leads to a stability constant with exponential dependence on the reciprocal value of the smallest eigenvalue of the diffusion matrix. Li [33] studied the steady problem (1.6), assuming that $\mu - \nabla \cdot \beta/2 \ge 0$. He further supposed the domain $\Omega$ to be a rectangle and admitted mixed boundary conditions. After transforming equation (1.6) into a system of first-order partial differential equations, he discretized this system by an upwind finite element method based on Raviart-Thomas finite elements on a structured grid. $H^1$-error estimates are derived, with constants depending on $\nu$ in an explicit way and via certain norms of the exact solution. Various least-squares finite element methods applied to (1.4) - again written as first-order system - were studied by Hsieh and Yang [30]. Requiring that $\beta$ is constant, these authors established error bounds whose dependence on $\nu$ is explicit and also indirect through some norms of the exact solution. Buffa, Hughes, Sangalli [6] proved stability and error estimates for a stabilized discontinuous Galerkin method applied to (1.4). Their stability bounds are $\nu$-uniform, and their error bounds depend on $\nu$ only insofar as certain norms of the exact solutions are concerned. However their estimates pertain to norms weighted by $\nu$ and $\beta$. Using norms weighted in a similar way, Verfürth [43] derived a posteriori error estimates for finite element discretizations of (1.6) in the case $\mu - \nabla \cdot \beta/2 \ge 0$. The unsteady case is treated in [44]. A lowest order discontinuous Petrov-Galerkin formulation with flux-upwind stabilization was applied to (1.4) by Causin, Sacco, Bottasso [8]. Since their discrete bilinear form associated with the convection term is positive, the authors obtained $H^1$-stability. The related constant inevitably depends on $\nu$. In addition, $H^1$-error estimates are shown. A
streamline-diffusion finite element method with a stability parameter incorporating the flow direction is proposed for (1.4) in [35]. An $H^1$-error estimate is proved. The unsteady variant of (1.6) with $\mu \in (0, \infty)$ and $\nabla \cdot \beta = 0$ is discretized by a least-squares finite element method in [32]. Stability holds with a constant which, due to Gronwall's lemma, depends exponentially on the reciprocal value of the smallest eigenvalue of the diffusion matrix.

There is still another way of deriving $L^2$-stability for certain numerical solutions to the unsteady problem (1.1). In fact, if a discrete analogue to the equation $\int_{\Omega}(\beta \cdot \nabla)u u dx = 0 (u \in H_0^1(\Omega))$ is valid, Gronwall's inequality yields a stability bound of the form $CT \|f\|_2$ on the time interval $[0, T]$, with $C$ independent of $\nu$. Burman and Fernández [7] applied this idea to discrete solutions of (1.1) constructed via edge-stabilization, local projection stabilization or quasi-static subscaling in standard finite element spaces of continuous piecewise polynomial functions. An extension to a nonlinear local projection stabilization is indicated in [4, Remark 4.5]. But of course, due to the factor $T$ in the stability constant, this approach does not yield estimates uniform in $t \in (0, \infty)$, nor does it provide an access to the steady problem (1.4).

We mention that numerical tests of various FE-FV discretizations of compressible flows are reported in [10], [23, Chapter 4.4] and [29]. Of course, combined FE-FV methods are not only applied to scalar convection-diffusion equations, but also to other problems. As examples, we indicate that in [17], [27], combined FE-FV methods based on the Crouzeix-Raviart finite element yield approximate solutions to a nonlinear version of the steady compressible Stokes system. These articles present convergence results. Reference [41] uses the scheme from [36] described above, applying it to the 2D steady incompressible Navier-Stokes system. Error estimates are presented for the case of Reynolds number 1.

## 2. Notation. FE-FV discretization of (1.1), (1.2) and (1.4), (1.5), respectively. Statement of main results.

If $\epsilon > 0$ and $x \in \mathbb{R}^2$, we write $B_{\epsilon}(x)$ for the open disk in $\mathbb{R}^2$ centered in $x$ and with radius $\epsilon$. If $U \subset \mathbb{R}^2$, we put $diam(U) := \sup\{|x - y| : x, y \in U\}$. We write $|U|$ for the measure of a measurable set $U \subset \mathbb{R}^2$. For $u, v \in L^2(\Omega)$, we use the abbreviation $(u, v) := \int_{\Omega} uv dx$. We write $\| \cdot \|_2$ for the standard norm of $L^2(\Omega)$. The Sobolev space $H^1(\Omega)$ and its subspace $H_0^1(\Omega)$ are defined in the usual way. The standard norm of $H^1(\Omega)$ is denoted by $\| \cdot \|_{1,2}$.

As already indicated in Section 1, we assume that $\Omega \subset \mathbb{R}^2$ is a bounded open polygon with Lipschitz boundary. The functions $\beta$, $\varphi$, $g$, $u_0$ and $G$ and the constants $\nu$, $T$, $\underline{\beta} \in (0, \infty)$ were also introduced in Section 1, with $\beta \in H^1(\Omega)^2$, $\varphi \in C^1(\bar{\Omega})$ being such that (1.3) is satisfied. We assume that $g \in C^0([0, T], H^1(\Omega))$ and $u^{(0)}, G \in H^1(\Omega)$. The functions $g(\cdot, t)$ and $G$ are required to belong to $H^1(\Omega)$ instead of only to $L^2(\Omega)$ so that they admit traces on edges, in view of an interpolation operator we will introduce below.

By adding a constant to $\varphi$, we may suppose without loss of generality that $\varphi(x) \ge \varphi_0 (x \in \bar{\Omega})$ for some $\varphi_0 > 0$. For example, in the case $\beta = \beta_0$ for some $\beta_0 \in \mathbb{R}^2 \setminus \{0\}$, we may put $\varphi(x) := 2|\beta_0|diam(\Omega) - \beta_0 \cdot (x - x_0)$, where $x_0$ is an arbitrary but fixed point in $\Omega$. Obviously there is a constant $\varphi_1 > 0$ with $\varphi(x) \le \varphi_1$ and $|\nabla\varphi(x)| \le \varphi_1$ for $x \in \bar{\Omega}$. We further introduce a parameter $\sigma_0 \in (0, 1)$, which will appear in condition (2.1) below. The set $\Omega$, the functions $\beta$, $\varphi$, $g$, $u_0$ and
G as well as the numbers $\nu$, $T$, $\varphi_0$, $\varphi_1$, $\underline{\beta}$ and $\sigma_0$ will be kept fixed throughout.

By the symbol $\mathfrak{C}$, we denote constants that may depend on $\sigma_0$, $\mathrm{diam}\,\Omega$, $\underline{\beta}$, $\varphi_0$ and $\varphi_1$, with polynomial dependence on $\underline{\beta}^{-1}$, $\varphi_0$ and $\varphi_1$. This last feature is important because we want to control how our estimates are influenced by $\beta$. This influence not only manifests itself by the factor $1 + \|\beta\|_{1,2}$ appearing in Theorem 2.1 and 2.2 and their proof, but also via the three quantities $\underline{\beta}$, $\varphi_0$ and $\varphi_1$.

We consider triangulations $\mathfrak{T}$ of $\Omega$ with the following three properties: Firstly, $\mathfrak{T}$ is a finite set of open triangles $K \subset \mathbb{R}^2$ with $\bar{\Omega} = \cup\{\bar{K} : K \in \mathfrak{T}\}$. Secondly, if $K_1, K_2 \in \mathfrak{T}$ with $\bar{K}_1 \cap \bar{K}_2 \neq \emptyset$ and $K_1 \neq K_2$, then $\bar{K}_1 \cap \bar{K}_2$ is a common vertex or a common side of $K_1$ and $K_2$. And thirdly, for any $K \in \mathfrak{T}$, the relation

$$B_{\sigma_0 \mathrm{diam} K}(x) \subset K \qquad (2.1)$$

is valid for some $x \in K$. All estimates appearing in the following involve constants that do not depend on the grid except via the parameter $\sigma_0$ in (2.1). Therefore we may simplify our presentation by fixing $\mathfrak{T}$ already at this point. So let $\mathfrak{T}$ be given such that the preceding assumptions, in particular (2.1), are verified. Below, when we introduce some spaces, norms and functions related to $\mathfrak{T}$, we will indicate this relation by an index $h$. This parameter $h$ may be considered as a quantity related to the size of the triangles $K \in \mathfrak{T}$, for example $h := \max\{\mathrm{diam} K : K \in \mathfrak{T}\}$, but it will not play any role in our theory. As a consequence of (2.1), we have

$$(\mathrm{diam} K)^2 \leq \mathfrak{C}|K| \quad \text{for } K \in \mathfrak{T}. \qquad (2.2)$$

Let $\mathfrak{S}$ be the set of the sides of the triangles $K \in \mathfrak{T}$. Put $J := \{1, \dots, \#\mathfrak{S}\}$, where $\#\mathfrak{S}$ denotes the number of elements of $\mathfrak{S}$. Let $(S_i)_{i \in J}$ be a numbering of $\mathfrak{S}$, and denote the midpoint of $S_i$ by $Q_i$ ($i \in J$). Set $J^o := \{i \in J : Q_i \in \Omega\}$, so that $J \setminus J^o = \{i \in J : Q_i \in \partial\Omega\}$. Note that for $i \in J \setminus J^o$, we have $S_i \subset \partial\Omega$.

We further introduce a barycentric mesh $(D_i)_{i \in J}$ on the triangular grid $\mathfrak{T}$: If $i \in J^o$, there are two triangles in $\mathfrak{T}$, denoted by $K_i^1, K_i^2$, such that $\overline{K_i^1} \cap \overline{K_i^2} = S_i$. We join the barycenter of each of these triangles with the endpoints of $S_i$. In this way we obtain a closed quadrilateral containing $S_i$ (Fig. 1). This quadrilateral is denoted by $D_i$. If $i \in J \setminus J^o$ (hence $Q_i \in \partial\Omega$), let $D_i$ be the closed triangle whose sides are the segment $S_i$ and the segments joining the endpoints of $S_i$ with the barycenter of the (unique) triangle $K \in \mathfrak{T}$ with $S_i \subset \bar{K}$. If $i, j \in J$ with $i \neq j$ are such that the set $D_i \cap D_j$ contains more than one point, then this set is a common side of $D_i$ and $D_j$. In this case, the quadrilaterals $D_i$ and $D_j$ are called “adjacent”, and their common side is denoted $\Gamma_{ij}$. For $i \in J$, we set

$$s(i) := \{j \in J \setminus \{i\} : D_i \text{ and } D_j \text{ are adjacent}\}.$$

If $i \in J$ and $j \in s(i)$, let $n_{ij}$ denote the outward unit normal to $D_i$ on $\Gamma_{ij}$. This means that $n_{ij}$ points from $D_i$ into $D_j$. We will use the abbreviation

$$\Theta_{ij}^{+} := \int_{\Gamma_{ij}} \max\{\beta(x) \cdot n_{ij}, 0\} \, do_x \quad (i \in J, j \in s(i)).$$
<figure><img src="image_2.png" /><figcaption>Figure 1: A quadrangle $\overline{K_i^1} \cup \overline{K_i^2}$ (left) and the quadrangle $D_i$ inside $\overline{K_i^1} \cup \overline{K_i^2}$ (right).</figcaption></figure>

Since $\nabla \cdot \beta = 0$ (see (1.3)) and $n_{ij} = -n_{ji}$ ($i \in J, j \in s(i)$), we get

$$ \sum_{j \in s(i)} (\Theta_{ij}^{+} - \Theta_{ji}^{+}) = \sum_{j \in s(i)} \int_{\Gamma_{ij}} \beta(x) \cdot n_{ij} \, do_x = 0 \quad (i \in J^o). \quad (2.3) $$

Due to relation (2.1) and because $K_i^1$ and $K_i^2$ have a common side, we obtain

$$ \operatorname{diam} K_i^{\nu} \le \mathcal{C} \operatorname{diam} K_i^{\mu} \quad \text{for } i \in J^o, \nu, \mu \in \{1, 2\}. \quad (2.4) $$

Next we introduce two finite element spaces by setting

$$ X_h := \{v \in L^2(\Omega) : v|_{K \in P_1(K)} \text{ for } K \in \mathcal{T}_h, v \text{ continuous at } Q_i \text{ for } i \in J\}, $$

$$ V_h := \{v_h \in X_h : v_h(Q_i) = 0 \text{ for } i \in J \setminus J^o\}, $$

where $P_1(A)$, for $A \subset \mathbb{R}^2$, denotes the set of all polynomials of degree at most 1 over A. The spaces $X_h$ and $V_h$ are nonconforming finite element spaces based on the piecewise linear Crouzeix-Raviart finite element. For $i \in J$, let $w_i$ be the function from $X_h$ that is uniquely determined by the requirement that $w_i(Q_j) = \delta_{ij}$ for $j \in J$. The family $(w_i)_{i \in J}$ is a basis of $X_h$. For $v_h, w_h \in X_h$, we put

$$ ((v_h, w_h))_{X_h} := \sum_{K \in \mathcal{T}} \int_K \nabla(v_h|K) \cdot \nabla(w_h|K) \, dx, \quad \|v_h\|_{X_h} := ((v_h, v_h))_{X_h}^{1/2}. $$

Obviously

$$ ((v_h, v_h))_{X_h} = \|v_h\|_{X_h}^2, \quad |((v_h, w_h))_{X_h}| \le \|v_h\|_{X_h} \|w_h\|_{X_h} \quad \text{for } v_h, w_h \in X_h. $$

The discrete Poincaré inequality

$$ \|v_h\|_2 \le \mathcal{C} \|v_h\|_{X_h} \quad \text{for } v_h \in V_h \quad (2.5) $$

was proved in [14], [31], [5] and [45], with a detailed analysis of the constant in [45].
From [1, (3.29), (3.31), (3.33)], we take two formulas for the L²-scalar product of $v_h, w_h \in X_h$, that is,

$$ (v_h, w_h) = (1/3) \sum_{K \in \mathcal{T}} |K| \sum_{r=1}^{3} v_h(Q_K^r) w_h(Q_K^r), \quad (v_h, w_h) = \sum_{i \in J} v_h(Q_i) w_h(Q_i) |D_i|, \quad (2.6) $$

where $Q_K^r \in \mathbb{R}^2$ for $r \in \{1, 2, 3\}$ are the three vertices of $K \in \mathcal{T}$. Put $H^1(\Omega) \oplus X_h := \{v + w_h : v \in H^1(\Omega), w_h \in X_h\}$, and let $I_h : H^1(\Omega) \oplus X_h \to X_h$ be the interpolation operator introduced in [20, 8.9.79]; it is defined by

$$ I_h(v) := \sum_{i \in J} l_i^{-1} \int_{S_i} v(x) \, do_x w_i \quad \text{for } v \in H^1(\Omega) \oplus X_h, $$

where $l_i$ denotes the length of $S_i$ ($i \in J$). Note that a function $v \in H^1(\Omega)$ admits a trace on $S_i$ for $i \in J$, and a function $v_h \in X_h$ verifies the equation

$$ \int_{S_i} E(v_h|K_i^1) \, do_x = l_i v_h(Q_i) = \int_{S_i} E(v_h|K_i^2) \, do_x \quad (i \in J^o), $$

where $E(v_h|K_i^s)$ denotes the continuous extension of $v_h|K_i^s$ to $S_i$ ($s \in \{1, 2\}$). Thus the operator $I_h$ is well defined. By [20, Lemma 8.9.81], it satisfies the estimate

$$ \|I_h(v)\|_2 \le \mathcal{C} \|v\|_{1,2} \quad \text{for } v \in H^1(\Omega). $$

It will be useful to introduce another interpolation operator besides $I_h$. In fact, for $v \in L^2(\Omega)$ with $v|K \in C^0(K)$ for $K \in \mathcal{T}_h$, $v$ continuous at $Q_i$ for $i \in J$, we set $\varrho_h(v) := \sum_{i \in J} v(Q_i)w_i$. Next we define a discrete convection term $b_h$, which is to approximate the variational form

$$ b(v, w) := \int_{\Omega} (\beta \cdot \nabla)vw \, dx \quad (v, w \in H^1(\Omega)) $$

associated with the convection term $\beta \cdot \nabla u$ in (1.1) and (1.4). We put

$$ b_h(v_h, w_h) := \sum_{i \in J} w_h(Q_i) \sum_{j \in s(i)} \left( \Theta_{ij}^{+} v_h(Q_i) - \Theta_{ji}^{+} v_h(Q_j) \right) \quad \text{for } v_h, w_h \in X_h. $$

This definition means that we discretize $b$ by an upwind finite volume method on the barycentric grid $(D_i)_{i \in J}$. A motivation for this discretization may be found in [15, p. 308], where a general numerical flux $H$ is considered instead of the upwind method we use here. If $\beta$ is constant, our upwind method corresponds to a choice of $H$ given by

$$ H(a, b, n) = a \text{ if } n \cdot \beta > 0, \quad H(a, b, n) = b \text{ if } n \cdot \beta \le 0, \quad \text{for } a, b \in \mathbb{R}, n \in \mathbb{R}^2 \text{ with } |n| = 1. $$

The approach here and the one in [15], as any upwind method, is based on the idea of taking into account the direction from which the information arrives. In view of discretizing the time variable, we fix $N \in \mathbb{N}$ and choose $t_1, \dots, t_N \in (0, T)$ with $t_1 < \dots < t_N$. Put $t_0 := 0, t_{N+1} :=$
T, $\tau_k := t_k - t_{k-1}$ for $1 \le k \le N + 1$. For brevity, we put $G_h := I_h(G)$ and introduce functions $g_h^{(k)}: \bar{\Omega} \to \mathbb{R}$ by setting

$$g_h^{(k)}(x) := I_h(g(\cdot, t_k))(x) \quad \text{for } k \in \{0, \dots, N+1\}, x \in \bar{\Omega}.$$

Now we are in a position to introduce the finite element – finite volume discretization of problem (1.1), (1.2) and (1.4), (1.5), respectively, that we want to study in the work at hand. Concerning (1.1), (1.2), we consider functions $u_h^{(0)}$, ..., $u_h^{(N+1)} \in V_h$ with

$$\tau_k^{-1}(u_h^{(k+1)} - u_h^{(k)}, v_h) + \nu((u_h^{(k+1)}, v_h))_{X_h} + b_h(u_h^{(k+1)}, v_h) = (g_h^{(k+1)}, v_h) \quad (2.7)$$
for $v_h \in V_h$, $k \in \{0, \dots, N\}$, $u_h^{(0)} = I_h(u^{(0)})$.

This scheme is implicit because both the diffusion and the convection term are discretized implicitly. For the steady problem (1.4), (1.5), we consider an approximate solution $U_h \in V_h$ satisfying

$$\nu((U_h, v_h))_{X_h} + b_h(U_h, v_h) = (G_h, v_h) \quad \text{for } v_h \in V_h. \quad (2.8)$$

In view of (2.5), and because $b_h(v_h, v_h) \ge 0$ for $v_h \in V_h$ (Lemma 3.3) and $\dim V_h < \infty$, both problems admit a unique solution. Our main results may now be stated as follows.

**Theorem 2.1** Let $u_h^{(0)}$, ..., $u_h^{(N+1)} \in V_h$ be a system of functions satisfying (2.7). Then

$$\left(\sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_2^2\right)^{1/2} + \max_{1 \le l \le N+1} \|u_h^{(l)}\|_2 + \nu^{1/2} \left(\sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_{X_h}^2\right)^{1/2} \le \mathcal{C}(1 + \|\beta\|_{1,2}) \left[\left(\sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2^2\right)^{1/2} + \|u_h^{(0)}\|_2\right]. \quad (2.9)$$

**Theorem 2.2** Let $U_h \in V_h$ be the solution of (2.8). Then

$$\|U_h\|_2 + \nu^{1/2} \|U_h\|_{X_h} \le \mathcal{C}(1 + \|\beta\|_{1,2}) \|G_h\|_2,$$

### 3. Auxiliary results.

We begin with a simple observation.

**Lemma 3.1** Let $a_{ij} \in \mathbb{R}$ for $i, j \in J$. Then $\sum_{i \in J} \sum_{j \in s(i)} a_{ij} = \sum_{i \in J} \sum_{j \in s(i)} a_{ji}$.

**Proof:** For $i, j \in J$, we have $j \in s(i)$ if and only if $i \in s(j)$. As a consequence, we obtain $\{(i, j) \in J^2 : j \in s(i)\} = \{(i, j) \in J^2 : i \in s(j)\}$. Thus the left-hand side of the equation in Lemma 3.1 equals $\sum_{j \in J} \sum_{i \in s(j)} a_{ij}$. Lemma 3.1 now follows by renaming indices. $\square$

The ensuing lemma deals with the interpolation operator $\varrho_h$.
**Lemma 3.2** Let $v_h, w_h \in X_h$. Then

$$ \begin{aligned} \|v_h\|_2^2 &\le \mathcal{C}(v_h, \varrho_h(v_h\varphi)), \quad (v_h, \varrho_h(v_h\varphi)) \le \mathcal{C}\|v_h\|_2^2, \\ (v_h, \varrho_h(w_h\varphi)) &\le (v_h, \varrho_h(v_h\varphi))^{1/2} (w_h, \varrho_h(w_h\varphi))^{1/2}, \\ \|\varrho_h(v_h\varphi)\|_{X_h} &\le \mathcal{C}\|v_h\|_{X_h}. \end{aligned} $$

**Proof:** The first three inequalities stated in Lemma 3.2 are an immediate consequence of (2.6), the relation $0 < \varphi_0 \le \varphi(x) \le \varphi_1$ ($x \in \Omega$) and – as concerns the third estimate – of the Cauchy-Schwarz inequality.

In order to prove the fourth, take $K \in \mathcal{T}$. There are three elements $i_1, i_2, i_3 \in J$ such that the points $Q_{i_r}$ with $r \in \{1, 2, 3\}$ are the three vertices of $K$. Then $\varrho_h(v_h\varphi)|K = \sum_{r=1}^3 (v_h\varphi)(Q_{i_r}) w_{i_r}|K$ by the definition of $\varrho_h(v_h\varphi)$. As a consequence

$$ \|\nabla(\varrho_h(v_h\varphi)|K)\|_2^2 = \int_K \left| \sum_{r=1}^3 (v_h\varphi)(Q_{i_r}) \nabla(w_{i_r}|K) \right|^2 dx. $$

An estimate as in the proof of [12, Lemma 3.1] now yields that $\|\nabla(\varrho_h(v_h\varphi)|K)\|_2^2$ is bounded by

$$ \begin{aligned} &2 \int_K \left| \sum_{r=1}^3 v_h(Q_{i_r}) \int_0^1 \left( \nabla\varphi(Q_{i_1} + \vartheta(Q_{i_r} - Q_{i_1})) (Q_{i_r} - Q_{i_1}) \right) d\vartheta \nabla(w_{i_r}|K) \right|^2 dx \\ &+ 2\varphi(Q_{i_1})^2 \|\nabla(v_h|K)\|_2^2. \end{aligned} $$

Therefore we obtain the upper bound

$$ \mathcal{C} \|\nabla \varphi\|_{\infty}^2 \sum_{r=1}^{3} v_h(Q_{i_r})^2 (\text{diam } K)^2 \int_K |\nabla(w_{i_r}|K)|^2 dx + 2 \|\varphi\|_{\infty}^2 \|\nabla(v_h|K)\|_2^2 $$

for $\|\nabla(\varrho_h(v_h\varphi)|K)\|_2^2$. But in view of (2.2) $|\nabla(w_{i_r}|K)| \le \mathcal{C}(\text{diam } K)^{-1}$, so by (2.6) and because $\|\nabla\varphi\|_{\infty} \le \varphi_1$, $\|\varphi\|_{\infty} \le \varphi_1$,

$$ \|\nabla(\varrho_h(v_h\varphi)|K)\|_2^2 \le \mathcal{C}(\|v_h|K\|_2^2 + \|\nabla(v_h|K)\|_2^2). $$

The fourth estimate in Lemma 3.2 now follows with (2.5). $\square$

Next we prove a basic fact regarding the discrete convection term $b_h$. The proof in question, as well as the proof of Lemma 3.4 below, applies techniques similar to those used by Sacco, Saleri [40].

**Lemma 3.3** For $v_h \in V_h$, the equation $b_h(v_h, v_h) = \mathfrak{R}_h/2$ holds, where

$$ \mathfrak{R}_h := \mathfrak{R}_h(v_h) := \sum_{i \in J} \sum_{j \in S(i)} \Theta_{ji}^{+} (v_h(Q_i) - v_h(Q_j))^2. $$

In particular $b_h(v_h, v_h) \ge 0$.
**Proof:** By (2.3) and because $v_h(Q_i) = 0$ for $i \in J \setminus J^o$, we can write

$$b_h(v_h, v_h) = \sum_{i \in J} v_h(Q_i) \sum_{j \in s(i)} (v_h(Q_i) - v_h(Q_j)) \Theta_{ji}^{+}.$$

Using the equation $v_h(Q_i)(v_h(Q_i) - v_h(Q_j)) = \frac{1}{2}v_h(Q_i)^2 + \frac{1}{2}(v_h(Q_i) - v_h(Q_j))^2 - \frac{1}{2}v_h(Q_j)^2$, we get

$$b_h(v_h, v_h) = \frac{1}{2} \sum_{i \in J} \sum_{j \in s(i)} (v_h(Q_i) - v_h(Q_j))^2 \Theta_{ji}^{+} + \frac{1}{2} \sum_{i \in J} \sum_{j \in s(i)} (v_h(Q_i)^2 - v_h(Q_j)^2) \Theta_{ji}^{+}.$$

On the other hand, by (2.3), the relation $v_h(Q_i) = 0$ for $i \in J \setminus J^o$ and Lemma 3.1, we obtain

$$\sum_{i \in J} \sum_{j \in s(i)} (v_h(Q_i)^2 - v_h(Q_j)^2) \Theta_{ji}^{+} = \sum_{i \in J} \sum_{j \in s(i)} (\Theta_{ij}^{+} v_h(Q_i)^2 - \Theta_{ji}^{+} v_h(Q_j)^2) = 0.$$

This concludes the proof.

<figure><img src="image_3.png" /></figure>

In the remainder of this section, we prove some sort of coercivity inequality for $b_h$ (Corollary 3.1). We begin by establishing a lower bound for $b_h(v_h, \varrho_h(v_h \varphi))$.

**Lemma 3.4** Let $v_h \in V_h$. Then $b_h(v_h, \varrho_h(v_h \varphi)) \ge A_h/2$, with

$$A_h := A_h(v_h) := \sum_{j \in J} v_h(Q_i)^2 \sum_{j \in s(i)} (\Theta_{ji}^{+} \varphi(Q_i) - \Theta_{ij}^{+} \varphi(Q_j)).$$

**Proof:** By (2.3) and the equation $v_h(Q_i) = 0$ for $i \in J \setminus J^o$,

$$b_h(v_h, \varrho_h(v_h \varphi)) = \sum_{i \in J} v_h(Q_i) \varphi(Q_i) \sum_{j \in s(i)} (v_h(Q_i) - v_h(Q_j)) \Theta_{ji}^{+}.$$

Thus, with an argument as in the proof of Lemma 3.3,

$$b_h(v_h, \varrho_h(v_h \varphi)) = \frac{1}{2} \sum_{i \in J} \varphi(Q_i) \sum_{j \in s(i)} (v_h(Q_i) - v_h(Q_j))^2 \Theta_{ji}^{+} \\ + \frac{1}{2} \sum_{i \in J} \varphi(Q_i) \sum_{j \in s(i)} (v_h(Q_i)^2 - v_h(Q_j)^2) \Theta_{ji}^{+}.$$

The first term in the above right hand side is non-negative, and the second equals $A_h/2$. In fact, by Lemma 3.1,

$$\sum_{i \in J} \varphi(Q_i) \sum_{j \in s(i)} v_h(Q_j)^2 \Theta_{ji}^{+} = \sum_{i \in J} v_h(Q_i)^2 \sum_{j \in s(i)} \Theta_{ij}^{+} (Q_j).$$

This concludes the proof.

<figure><img src="image_4.png" /></figure>

**Lemma 3.5** Let $v_h \in V_h$, and put

$$B_h := B_h(v_h) := - \sum_{i \in J} v_h(Q_i)^2 \sum_{j \in s(i)} \int_{\Gamma_{ij}} \varphi(x) \beta(x) \cdot n_{ij} \, dx.$$

Then $B_h \ge \underline{\beta} \|v_h\|_2^2$.
**Proof:** By the Divergence theorem and (1.3),

$$ - \sum_{j \in s(i)} \int_{\Gamma_{ij}} \varphi(x) \beta(x) \cdot n_{ij} \, dx_x = - \int_{D_i} \beta(x) \cdot \nabla \varphi(x) \, dx \geq |D_i| \underline{\beta} \quad \text{for } i \in J^o. $$

The lemma follows with (2.6) and the equation $v_h(Q_i) = 0$ for $i \in J \setminus J^o$. $\square$

**Lemma 3.6** Let $v_h \in V_h$. Then

$$ B_h \leq 2b_h(v_h, \varrho_h(v_h\varphi)) + \mathfrak{C} \|\beta\|_{1,2}^{1/2} \|v_h\|_2 \mathfrak{R}_h^{1/2}, $$

where $B_h = B_h(v_h)$ and $\mathfrak{R}_h = \mathfrak{R}_h(v_h)$ were introduced in Lemma 3.5 and 3.3, respectively.

**Proof:** By Lemma 3.4, we have

$$ B_h = B_h - A_h + A_h \leq |A_h - B_h| + 2b_h(v_h, \varrho_h(v_h\varphi)), \quad (3.1) $$

with $A_h$ from that lemma. On the other hand, since $n_{ij} = -n_{ji}$ ($i \in J, j \in s(i)$),

$$ \begin{aligned} |A_h - B_h| &= \left| \sum_{i \in J} v_h(Q_i)^2 \sum_{j \in s(i)} \left[ \int_{\Gamma_{ij}} (\varphi(Q_i) - \varphi(x)) \max\{\beta(x) \cdot n_{ji}, 0\} \, dx_x \right. \right. \\ & \qquad \left. \left. - \int_{\Gamma_{ij}} (\varphi(Q_j) - \varphi(x)) \max\{\beta(x) \cdot n_{ij}, 0\} \, dx_x \right] \right| \\ &= \left| \sum_{i \in J} \sum_{j \in s(i)} (v_h(Q_i)^2 - v_h(Q_j)^2) \int_{\Gamma_{ij}} (\varphi(Q_i) - \varphi(x)) \max\{\beta(x) \cdot n_{ji}, 0\} \, dx_x \right|, \end{aligned} \quad (3.2) $$

where the last equation follows from Lemma 3.1. For $i \in J, j \in s(i)$ and $x \in \Gamma_{ij}$, we find with the relation $|\nabla\varphi(x)| \leq \varphi_1(x \in \bar{\Omega})$ that

$$ |\varphi(Q_i) - \varphi(x)| \leq \varphi_1 \max_{\nu \in \{1, 2\}} \mathrm{diam} K_i^{\nu}. $$

Here and in the following, if $i \in J \setminus J^o$, we use both the notation $K_i^1$ and $K_i^2$ for the unique triangle $K \in \mathfrak{T}$ with $Q_i \in K$. Now we may deduce from (3.2) that

$$ |A_h - B_h| \leq \varphi_1 \sum_{i \in J} \sum_{j \in s(i)} \Theta_{ji}^{+} \max_{\nu \in \{1, 2\}} \mathrm{diam} K_i^{\nu} |v_h(Q_j) + v_h(Q_i)| |v_h(Q_j) - v_h(Q_i)|, $$

hence by the Cauchy-Schwarz inequality,

$$ |A_h - B_h| \leq \varphi_1 \left( \sum_{i \in J} \sum_{j \in s(i)} \Theta_{ji}^{+} \left( \max_{\nu \in \{1, 2\}} \mathrm{diam} K_i^{\nu} \right)^2 (v_h(Q_j) + v_h(Q_i))^2 \right)^{1/2} \mathfrak{R}_h^{1/2}. $$

For $i \in J, j \in s(i)$, there is $\nu \in \{1, 2\}$ with $\Gamma_{ji} \subset K_j^{\nu}$, so

$$ \Theta_{ji}^{+} \leq \int_{\Gamma_{ji}} |\beta(x)| \, dx \leq (\mathrm{diam} K_j^{\nu})^{1/2} \|\beta|\Gamma_{ji}\|_2 \leq \mathfrak{C} \|\beta|D_j \cap K_j^{\nu}\|_{1,2} \leq \mathfrak{C} \cdot \|\beta\|_{1,2}, $$
where the second last inequality follows by a change of the domain of integration from $D_j \cap K_j^\nu$ to the unit triangle $\{x \in [0, 1]^2 : x_1 + x_2 \le 1\}$, an application of a standard trace inequality on that triangle, and then another change of the domain of integration back to $D_j \cap K_j^\nu$. Concerning that last transformation, note that $(\mathrm{diam} K_j^\nu)^2 \le \mathcal{C}|K_j^\nu|$ by (2.2) and $|K_j^\nu|/3 = |D_j \cap K_j^\nu|$. As a consequence, we arrive at the estimate

$$|A_h - B_h| \le \mathcal{C} \|\beta\|_{1,2}^{1/2} \left( \sum_{i \in J} \sum_{j \in s(i)} \left( \max_{\nu \in \{1,2\}} \mathrm{diam} K_i^\nu \right)^2 \left( v_h(Q_j)^2 + v_h(Q_i)^2 \right) \right)^{1/2} \mathfrak{K}_h^{1/2}. \quad (3.3)$$

On the other hand, with (2.4),

$$\sum_{i \in J} \sum_{j \in s(i)} \left( \max_{\nu \in \{1,2\}} \mathrm{diam} K_i^\nu \right)^2 \left( v_h(Q_j)^2 + v_h(Q_i)^2 \right) \\ \le \sum_{i \in J} \sum_{j \in s(i)} \left( (\mathrm{diam} K_j^1)^2 v_h(Q_j)^2 + (\mathrm{diam} K_i^1)^2 v_h(Q_i)^2 \right).$$

Therefore by Lemma 3.1,

$$\sum_{i \in J} \sum_{j \in s(i)} \left( \max_{\nu \in \{1,2\}} \mathrm{diam} K_i^\nu \right)^2 \left( v_h(Q_j)^2 + v_h(Q_i)^2 \right) \le 2 \sum_{i \in J} \sum_{j \in s(i)} (\mathrm{diam} K_i^1)^2 v_h(Q_i)^2.$$

Again we use that that $(\mathrm{diam} K_i^\nu)^2 \le \mathcal{C}|K_i^\nu|$ by (2.2) and $|K_i^\nu|/3 = |D_i \cap K_i^\nu|$ for $i \in J^o$, $\nu \in \{1, 2\}$. It follows with (2.6) that

$$\sum_{i \in J} \sum_{j \in s(i)} \left( \max_{\nu \in \{1,2\}} \mathrm{diam} K_i^\nu \right)^2 \left( v_h(Q_j)^2 + v_h(Q_i)^2 \right) \le \mathcal{C} \sum_{i \in J} v_h(Q_i)^2 |D_i| \le \mathcal{C} \|v_h\|_2^2.$$

Thus we obtain from (3.3) that $|A_h - B_h| \le \mathcal{C}\|\beta\|_{1,2}^{1/2}\|v_h\|_2 \mathfrak{K}_h^{1/2}$. Lemma 3.6 now follows from (3.1). $\square$

**Corollary 3.1** Let $v_h \in V_h$. Then

$$\|v_h\|_2^2 \le \mathcal{C}(1 + \|\beta\|_{1,2}) (b_h(v_h, \varrho(v_h, \varphi)) + b_h(v_h, v_h)). \quad (3.4)$$

**Proof:** Lemma 3.5, 3.6 and 3.3 imply

$$\|v_h\|_2^2 \le \mathcal{C} \left( b_h(v_h, \varrho(v_h, \varphi)) + \|\beta\|_{1,2}^{1/2} b_h(v_h, v_h)^{1/2} \|v_h\|_2 \right),$$

so that

$$\|v_h\|_2^2 \le \mathcal{C} \left( b_h(v_h, \varrho(v_h, \varphi)) + \|\beta\|_{1,2} b_h(v_h, v_h) \right) + \|v_h\|_2^2/2.$$

The corollary follows from this estimate. $\square$
## 4. Proof of Theorem 2.1 and 2.2.

The largest part of this section is taken by the

**Proof of Theorem 2.1:** Putting $v_h = u_h^{(k+1)}$ in (2.7), we get

$$ \begin{aligned} & \nu \|u_h^{(k+1)}\|_{X_h}^2 + \tau_{k+1}^{-1} \|u_h^{(k+1)}\|_{2}^2 + b_h(u_h^{(k+1)}, u_h^{(k+1)}) \\ & \le (g_h^{(k+1)}, u_h^{(k+1)}) + (2\tau_{k+1})^{-1} \|u_h^{(k+1)}\|_{2}^2 + (2\tau_{k+1})^{-1} \|u_h^{(k)}\|_{2}^2 \quad \text{for } 0 \le k \le N. \end{aligned} $$

Taking the sum with respect to $k \in \{0, ..., s\}$, for $s \in \{0, ..., N\}$, and using a simple shoestring argument, we find with $s = N$ that

$$ \sum_{l=1}^{N+1} \tau_l b_h(u_h^{(l)}, u_h^{(l)}) + \nu \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_{X_h}^2 \le \sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 + \|u_h^{(0)}\|_2^2/2, \quad (4.1) $$

and for $s \in \{0, ..., N\}$,

$$ \begin{aligned} & \sum_{l=1}^{s+1} \tau_l (b_h(u_h^{(l)}, u_h^{(l)}) + \nu \|u_h^{(l)}\|_{X_h}^2) + \|u_h^{(s+1)}\|_2^2/2 \\ & \le \sum_{l=1}^{s+1} \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 + \|u_h^{(0)}\|_2^2/2. \end{aligned} $$

The term $\|u_h^{(N+1)}\|_2^2/2$ was dropped on the left-hand side of (4.1) because it will not be needed. Since $b_h(v_h, v_h) \ge 0$ for $v_h \in V_h$ (Lemma 3.3), we may conclude that

$$ \max_{1 \le l \le N+1} \|u_h^{(l)}\|_2^2/2 \le \sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 + \|u_h^{(0)}\|_2^2/2. \quad (4.2) $$

Next we use (2.7) once more, this time with $v_h = \varrho_h(u_h^{(k+1)}\varphi)$, to obtain with Lemma 3.2 that

$$ \begin{aligned} & \nu \tau_{k+1} ((u_h^{(k+1)}, \tilde{u}_h^{(k+1)}))_{X_h} + \tau_{k+1} b_h(u_h^{(k+1)}, \tilde{u}_h^{(k+1)}) + (u_h^{(k+1)}, \tilde{u}_h^{(k+1)}) \\ & \le \tau_{k+1} (g_h^{(k+1)}, \tilde{u}_h^{(k+1)}) + (u_h^{(k)}, \tilde{u}_h^{(k)})^{1/2} (u_h^{(k+1)}, \tilde{u}_h^{(k+1)})^{1/2} \\ & \le \tau_{k+1} (g_h^{(k+1)}, \tilde{u}_h^{(k+1)}) + (u_h^{(k+1)}, \tilde{u}_h^{(k+1)})/2 + (u_h^{(k)}, \tilde{u}_h^{(k)})/2 \quad \text{for } k \in \{0, ..., N\}. \end{aligned} $$

Here and in the following, we use the abbreviation $\tilde{u}_h^{(k)} := \varrho_h(u_h^{(k)}\varphi)$ for $0 \le k \le N+1$. Now we again take the sum with respect to $k \in \{0, ..., N\}$. By a simple shoestring argument and because $(u_h^{(N+1)}, \tilde{u}_h^{(N+1)}) \ge 0$ according to Lemma 3.2, we then get

$$ \sum_{l=1}^{N+1} \tau_l b_h(u_h^{(l)}, \tilde{u}_h^{(l)}) \le \sum_{l=1}^{N+1} \tau_l |(g_h^{(l)}, \tilde{u}_h^{(l)})| + \nu \sum_{l=1}^{N+1} \tau_l |((u_h^{(l)}, \tilde{u}_h^{(l)}))_{X_h}| + (u_h^{(0)}, \tilde{u}_h^{(0)})/2. \quad (4.3) $$

On the other hand, from Corollary 3.1 with $v_h = u_h^{(l)}$,

$$ \|u_h^{(l)}\|_2^2 \le \mathcal{C}(1 + \|\beta\|_{1,2}) (b_h(u_h^{(l)}, \tilde{u}_h^{(l)}) + b_h(u_h^{(l)}, u_h^{(l)})) \quad \text{for } 1 \le l \le N+1. $$
Therefore, with (4.1) and (4.3),

$$ \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_2^2 \le \mathcal{C}(1 + \|\beta\|_{1,2}) \mathfrak{R}, \quad (4.4) $$

where

$$ \mathfrak{R} := \sum_{l=1}^{N+1} \tau_l |(g_h^{(l)}, \tilde{u}_h^{(l)})| + \nu \sum_{l=1}^{N+1} \tau_l |((u_h^{(l)}, \tilde{u}_h^{(l)}))_{X_h}| + (u_h^{(0)}, \tilde{u}_h^{(0)}) \\ + \sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 + \|u_h^{(0)}\|_2^2. $$

By adding the left- and right-hand side of (4.1) and (4.2) to respectively the left- and right-hand side of (4.4), and taking account of the fact that $b_h(v_h, v_h) \ge 0$ for $v_h \in V_h$ (Lemma 3.3), we get

$$ \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_2^2 + \max_{1 \le l \le N+1} \|u_h^{(l)}\|_2^2 + \nu \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_{X_h}^2 \le \mathcal{C}(1 + \|\beta\|_{1,2}) \mathfrak{R}. \quad (4.5) $$

But with Lemma 3.2,

$$ |(g_h^{(l)}, \tilde{u}_h^{(l)})| \le (g_h^{(l)}, \varrho_h(g_h^{(l)} \varphi))^{1/2} (u_h^{(l)}, \tilde{u}_h^{(l)})^{1/2} \le \mathcal{C} \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 \quad \text{for } 1 \le l \le N+1, \quad (4.6) $$

$$ (u_h^{(0)}, \tilde{u}_h^{(0)}) \le \mathcal{C} \|u_h^{(0)}\|_2^2, \quad (4.7) $$

$$ \nu \sum_{l=1}^{N+1} \tau_l |((u_h^{(l)}, \tilde{u}_h^{(l)}))_{X_h}| \le \nu \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_{X_h} \|\tilde{u}_h^{(l)}\|_{X_h} \le \mathcal{C} \nu \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_{X_h}^2 \quad (4.8) $$

$$ \le \mathcal{C} \left( \sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 + \|u_h^{(0)}\|_2^2 \right), $$

where the last inequality holds according to (4.1). Thus the left-hand side of (4.5) is bounded by $\mathcal{C}(1 + \|\beta\|_{1,2}) (\sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 + \|u_h^{(0)}\|_2^2)$. Noting that

$$ \mathcal{C}(1 + \|\beta\|_{1,2}) \tau_l \|g_h^{(l)}\|_2 \|u_h^{(l)}\|_2 \le \mathcal{C}(1 + \|\beta\|_{1,2})^2 \tau_l \|g_h^{(l)}\|_2^2 + \tau_l \|u_h^{(l)}\|_2^2/2 \quad (4.9) $$

for $1 \le l \le N+1$, we may deduce from (4.5) that the left-hand side of (4.5) is bounded by

$$ \mathcal{C}(1 + \|\beta\|_{1,2})^2 \left( \sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2^2 + \|u_h^{(0)}\|_2^2 \right) + \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_2^2/2. \quad (4.10) $$

Therefore again from (4.5), by a simple shoestring argument,

$$ \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_2^2/2 + \max_{1 \le l \le N+1} \|u_h^{(l)}\|_2^2 + \nu \sum_{l=1}^{N+1} \tau_l \|u_h^{(l)}\|_{X_h}^2 \\ \le \mathcal{C}(1 + \|\beta\|_{1,2})^2 \left( \sum_{l=1}^{N+1} \tau_l \|g_h^{(l)}\|_2^2 + \|u_h^{(0)}\|_2^2 \right). $$
This proves Theorem 2.1.

□

We finally turn to the

**Proof of Theorem 2.2:** From (2.8) with $v_h$ replaced by $U_h$, we get

$$ \nu \|U_h\|_{X_h}^2 + b_h(U_h, U_h) = (G_h, U_h) \le \|G_h\|_2 \|U_h\|_2. \quad (4.11) $$

Thus Corollary 3.1 with $U_h$ in the place of $v_h$ yields

$$ \|U_h\|_2^2 \le \mathcal{C}(1 + \|\beta\|_{1,2}) (b_h(U_h, W_h) + \|G_h\|_2 \|U_h\|_2), \quad (4.12) $$

where we used the abbreviation $W_h := \rho_h(U_h, \varphi)$. Moreover, from (2.8) with $v_h = W_h$,

$$ b_h(U_h, W_h) \le \nu |((U_h, W_h))_{X_h}| + |(G_h, W_h)|. \quad (4.13) $$

On the other hand, with Lemma 3.2,

$$ |(G_h, W_h)| \le (G_h, \rho_h(G_h, \varphi))^{1/2} (U_h, W_h)^{1/2} \le \mathcal{C} \|G_h\|_2 \|U_h\|_2, \quad (4.14) $$

and by Lemma 3.2, 3.3 and (4.11),

$$ \nu |((U_h, W_h))_{X_h}| \le \mathcal{C} \nu \|U_h\|_{X_h}^2 \le \mathcal{C} (\nu \|U_h\|_{X_h}^2 + b_h(U_h, U_h)) \le \mathcal{C} \|G_h\|_2 \|U_h\|_2. $$

This estimate together with (4.12) – (4.14) yield

$$ \|U_h\|_2^2 \le \mathcal{C}(1 + \|\beta\|_{1,2}) \|G_h\|_2 \|U_h\|_2, $$

hence with (4.11) and the relation $b_h(U_h, U_h) \ge 0$ (Lemma 3.3),

$$ \|U_h\|_2^2 + \nu |((U_h, W_h))_{X_h}| \le \mathcal{C}(1 + \|\beta\|_{1,2}) \|G_h\|_2 \|U_h\|_2 \le \mathcal{C}(1 + \|\beta\|_{1,2})^2 \|G_h\|_2^2 + \|U_h\|_2/2. $$

Theorem 2.2 follows from this inequality.

□

## 5. Final comments.

It is natural to ask what are the minimal hypotheses on $\beta$ such that Theorem 2.1 and 2.2 remain valid. In this respect we first remark that the assumptions $\beta \in C^1(\bar{\Omega})^2$, $\nabla \cdot \beta = 0$ are not sufficient for Theorem 2.2 to hold.

To see this, consider the case $\Omega = \{x \in \mathbb{R}^2 : |x| < 1\}$, $\beta(x) = (-x_2, x_1)$ for $x \in \bar{\Omega}$. Obviously $\beta \in C^1(\bar{\Omega})^2$ and $\nabla \cdot \beta = 0$. Note that since $\beta(x) \to 0$ for $|x| \to 0$, there is no function $\varphi \in C^1(\bar{\Omega})$ such that the second relation in (1.3) is valid. We further put $u_\nu(x) := |x|^2/\nu - 1/\nu$ for $x \in \bar{\Omega}$, $\nu \in (0, \infty)$. Then $-\nu \Delta u_\nu + \beta \cdot \nabla u_\nu = -4$ and $u_\nu|_{\partial\Omega} = 0$ for any $\nu \in (0, \infty)$. Since $\|u_\nu\|_2 = \sqrt{\pi/(3\nu)}$ for $\nu \in (0, \infty)$, it follows there cannot be a constant $C > 0$ such that $\|u\|_2 \le C \|f\|_2$ for $\nu \in (0, \infty)$, $f \in L^2(\Omega)$ and $u \in H^2(\Omega)$ the solution of the boundary value problem

$$ -\nu \Delta u + \beta \cdot \nabla u = f \quad \text{in } \Omega, \quad u|_{\partial\Omega} = 0. \quad (5.1) $$
With some additional effort, the preceding example may be extended to the case that $\Omega$ is an annular domain, for example $\Omega = \{y \in \mathbb{R}^2 : 1/2 < |y| < 1\}$, and to rotationally symmetric functions $f \in L^2(\Omega)$. In fact, for such a domain $\Omega$ and such a function $f$, the solution $u \in H^2(\Omega)$ of the boundary value problem

$$-\Delta u = f/\nu \quad \text{in } \Omega, \quad u|_{\partial\Omega} = 0, \qquad (5.2)$$

is rotationally symmetric as well. Due to our choice of $\beta$, this means that $u$ solves problem (5.1). On the other hand, $u$ as a solution of (5.2) may be represented as a series involving eigenfunctions of the Laplacian on $\Omega$. This representation shows that $\|u\|_2$ cannot be bounded uniformly in $\nu > 0$ when $\nu$ tends to zero.

It is true that in the above cases $\Omega$ is not a polygon, and that the exact solution of (5.1) and an approximate solution provided by (2.8) are not the same thing. Still the preceding observations on $L^2$-estimates of solutions to (5.1) indicate it cannot be expected that an inequality as in Theorem 2.2 holds with a constant $C$ independent of $h$ and $\nu$ if $\beta \in H^1(\Omega)^2$ is only supposed to satisfy the relation $\nabla \cdot \beta = 0$.

Going a step further, one might conjecture that if $\nabla \cdot \beta = 0$, the second condition in (1.3) is not only sufficient (as we have shown), but even necessary for existence of a constant $C$ as in Theorem 2.2. In fact, if such a constant exists, it should be expected there is $C > 0$ with $\|u\|_2 \le C\|f\|_2$ for $\nu \in (0, \infty)$, $f \in L^2(\Omega)$ and $u \in H^2(\Omega)$ the solution of (5.1). This means in particular that the preceding constant $C$ is independent of $\nu$. In this situation, let $\beta \in (0, \infty)$ and let $(\nu_n)$ be a sequence in $(0, \infty)$ with $\nu_n \to 0$. Let the solution of (5.1) with $f$, $\nu$ replaced by $\beta$, $\nu_n$, respectively, be denoted by $u_n$, for $n \in \mathbb{N}$. Then, due to the estimate $\|u\|_2 \le C\|f\|_2$ for $\nu$, $f$, $u$ as above. the sequence $(\|u_n\|_2)$ is bounded and thus possesses a subsequence which is weakly convergent in $L^2(\Omega)$. Let this subsequence be also denoted by $(u_n)$. Let $u \in L^2(\Omega)$ be the limit function of this sequence, and take $\phi \in C_0^\infty(\Omega)$. Choosing $f = -\beta$, $u = u_n$ in (5.1), multiplying this equation by $\phi$, integrating over $\Omega$ and performing some integrations by parts, we arrive at the equation

$$\int_{\Omega} u_n (-\nu_n \Delta \phi - \beta \cdot \nabla \phi + \underline{\beta} \phi) = 0 \quad \text{for } n \in \mathbb{N}.$$

Letting $n$ tend to infinity, we see that $\int_{\Omega} u(-\beta \cdot \nabla \phi + \underline{\beta} \phi) dx = 0$. This means that $u$ satisfies the equation $-\beta \cdot \nabla u = \underline{\beta}$ in a weak sense. Thus, if $u \in C^1(\bar{\Omega})$, the second condition in (1.3) would hold with $u$ in the role of $\varphi$. However, the question of how to prove $C^1$-regularity of $u$ is open.

## References

[1] P. Angot, V. Dolejší, M. Feistauer and J. Felcman, Analysis of a combined barycentric finite volume – nonconforming finite element method for nonlinear convection-diffusion problems. Appl. Math. 43, 263-310 (1998).

[2] T. Arbogast, M. F. Wheeler and N. Y. Zhang, A nonlinear mixed finite element method for a degenerate parabolic equation arising in flow in porous media. SIAM J. Numer. Anal. 33, 1669-1687 (1996).
[3] B. Ayuso and L. D. Marini, Discontinuous Galerkin methods for advection-diffusion-reaction problems. SIAM J. Numer. Anal. 47, 1391-1420 (2009).

[4] G. R. Barrenechea, V. John and P. Knobloch, A local projection stabilization finite element method with nonlinear crosswind diffusion for convection-diffusion-reaction equation. Math. Model. Numer. Anal. 47 (2013), 1335-1366.

[5] S. C. Brenner, Poincaré-Friedrichs inequalities for piecewise $H^1$-functions. SIAM J. Numer. Anal. 41, 306-324 (2003).

[6] A. Buffa, T. J. R. Hughes and G. Sangalli, Analysis of multiscale discontinuous Galerkin method for convection-diffusion problems. SIAM J. Numer. Anal. 44, 1420-1440 (2006).

[7] E. Burman and M. A. Fernández, Finite element methods with symmetric stabilization for the transient convection-diffusion-reaction equation. Comput. Methods Appl. Mech. Engrg. 198 (2009), 2508-2519.

[8] P. Causin, R. Sacco and C. L. Bottasso, Flux-upwind stabilization of the discontinuous Petrov-Galerkin formulation with Lagrangian multipliers for advection-diffusion problems. ESAIM: Math. Model. Numer. Anal. 39, 1087-1114 (2005).

[9] S. Champier, T. Gallouët and R. Herbin, Convergence of an upstream finite volume scheme for a nonlinear hyperbolic equation on a triangular mesh. Numer. Math. 66, 139-157 (1993).

[10] C. Debiez, A. Dervieux, K. Mer and B. Nkonga, Computation of unsteady flows with mixed finite volume/finite element upwind methods. Int. J. Numer. Methods Fluids 27, 193-206 (1998).

[11] P. Deuring and M. Mildner, Error estimates for a finite element – finite volume discretization of convection-diffusion equations. Applied Numer. Math. 61, 785-801 (2011).

[12] P. Deuring and M. Mildner, Stability of a combined finite element – finite volume discretization of convection-diffusion equations. Numer. Meth. Partial Diff. Equ. 28, 402-424 (2012).

[13] A. Devinatz, R. Ellis and A. Friedman, The asymptotic behavior of the first real eigenvalue of second order elliptic operators with a small parameter in the highest derivative. II, Indiana Univ. Math. J. 23, 991-1011 (1973).

[14] V. Dolejší, M. Feistauer and J. Felcman, On the discrete Friedrichs inequality for non-conforming finite elements. Numer. Funct. Anal. Optim. 20, 437-447 (1999).

[15] V. Dolejší, M. Feistauer, J. Felcman and A. Kliková, Error estimates for barycentric finite volumes combined with nonconforming finite elements applied to nonlinear convection-diffusion problems. Appl. Math. 47, 301-340 (2002).

[16] A. Ern and J.-L. Guermond, Theory and practice of finite elements. Applied Mathematical Sciences 159, Springer, New York e.a., 2004.
[17] R. Eymard, T. Gallouët, R. Herbin and J.-C. Latché, A convergent finite element - finite volume scheme for the compressible Stokes problem. Part II: the isentropic case. Math. Comp. 270, 649-675 (2010).

[18] R. Eymard, D. Hilhorst and M. Vohralík, A combined finite volume-nonconforming/mixed-hybrid finite element scheme for degenerate parabolic problems. Numer. Math. 105, 73-131 (2006).

[19] R. Eymard, D. Hilhorst and M. Vohralík, A combined finite volume - finite element scheme for the discretization of strongly nonlinear convection-diffusion-reaction problems on non-matching grids. Numer. Meth. Partial Diff. Equ. 26, 612-646 (2010).

[20] M. Feistauer, Mathematical methods in fluid dynamics. Pitman Monographs and Surveys in Pure and Applied Mathematics 67, Longman, Harlow, 1993.

[21] M. Feistauer, J. Felcman and M. Lukáčová-Medvidová, On the convergence of a combined finite volume – finite element method for nonlinear convection-diffusion problems. Numer. Meth. Partial Diff. Equ. 13, 163-190 (1997).

[22] M. Feistauer, J. Felcman, M. Lukáčová-Medvidová and G. Warnecke, Error estimates for a combined finite volume – finite element method for nonlinear convection-diffusion problems. SIAM J. Math. Anal. 36, 1528-1548 (1999).

[23] M. Feistauer, J. Felcman and I. Straškraba, Mathematical and computational methods for compressible flow. Clarendon Press, Oxford, 2003.

[24] M. Feistauer, V. Kučera, K. Najzan and J. Prokopová, Analysis of a space-time discontinuous Galerkin method for nonlinear convection-diffusion problems. Numer. Math. 117, 251-288 (2011).

[25] M. Feistauer, J. Slavík and P. Stupka, On the convergence of a combined finite volume – finite element method for nonlinear convection-diffusion problems. Explicit schemes. Numer. Meth. Partial Diff. Equ. 15, 215-235 (1999).

[26] A. Friedman, Partial Differential Equations of Parabolic Type. Dover, Mineola N.Y., 2008.

[27] T. Gallouët, R. Herbin and J.-C. Latché, A convergent finite element - finite volume scheme for the compressible Stokes problem. Part I: the isothermal case. Math. Comp. 267, 1333-1352 (2009).

[28] D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order. Springer, Berlin e.a., 1983.

[29] L. Hallo, C. Le Ribault and M. Buffat, An implicit mixed finite-volume-finite-element method for solving 3D turbulent compressible flows. Int. J. Numer. Meth. Fluids 25, 1241-1261 (1997).

[30] Po-Wen Hsieh and Suh-Yuh Yang, On efficient least-squares finite element methods for convection-dominated problems. Comp. Meth. Appl. Mech. Eng. 199, 183-196 (2009).
[31] P. Knobloch, Uniform validity of discrete Friedrich's inequality for general nonconforming finite element spaces. Numer. Funct. Anal. Optimiz. 22, 107-126 (2001).

[32] R. C. Leal Toledo and V. Ruas, Numerical analysis of a least-squares finite element method for the time-dependent advection-diffusion equation. J. Comp. Appl. Math. 235, 3615-3631 (2011).

[33] Zhitao Li, Convergence analysis of an upwind mixed element method for advection diffusion problems. Appl. Math. Comp. 212, 318-326 (2009).

[34] M. Mildner, Stabilité de l'équation d'advection-diffusion et stabilité de l'équation d'advection pour la solution du problème approché obtenue par la méthode upwind d'éléments-finis et volumes-finis avec des éléments de Crouzeix-Raviart. Thesis. Université du Littoral, Calais, 2013.

[35] A. S. Mounim, A stabilized finite element method for convection-diffusion problems. Numer. Meth. Part. Diff. Equ. 28, 1916-1943 (2012).

[36] K. Ohmori and T. Ushijima, A technique of upstream type applied to a linear nonconforming finite element approximation of convective diffusion equations. RAIRO Modél. Math. Anal. Numér. 18, 309-332 (1984).

[37] A. Quarteroni and A. Valli, Numerical Approximation of Partial Differential Equations. Springer Series in Computational Mathematics 23, Springer, Berlin e.a., 1994.

[38] Yongqiang Ren, Aijie Cheng and Hong Wang, A uniformly optimal-order estimate for finite volume method for advection-diffusion equation. Numer. Meth. Partial Diff. Equ. 30, 17-43 (2014).

[39] H.-G. Roos, M. Stynes and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations (2nd ed.). Springer Series in Computational Mathematics 24, Springer, New York e.a., 2008.

[40] R. Sacco and F. Saleri, Stabilized mixed finite volume methods for convection-diffusion problems. East West J. Numer. Math. 5, 291-311 (1997).

[41] F. Schieweck and L. Tobiska, A nonconforming finite element method of upstream type applied to the steady Navier-Stokes equation. RAIRO Modél. Math. Anal. Numér. 23, 627-647 (1989).

[42] A. Szepessy, Convergence of a streamline diffusion finite element method for scalar conservation laws with boundary conditions. RAIRO Modél. Math. Anal. Numér. 25, 749-782 (1991).

[43] R. Verfürth, Robust a posteriori error estimates for steady convection-diffusion equations. SIAM J. Numer. Anal. 43, 1766-1782 (2005).

[44] R. Verfürth, Robust a posteriori error estimates for unsteady convection-diffusion equations. SIAM J. Numer. Anal. 43, 1783-1802 (2005).
[45] M. Vohralík, On the discrete Poincaré-Friedrichs inequalities for nonconforming approximations of the Sobolev space $H^1$. Numer. Funct. Anal. Optimiz. 26, 925-952 (2005).