Study of a finite volume - finite element scheme for a nuclear transport model

TL;DR

This paper develops and analyzes a finite volume - finite element scheme for a complex, nonlinear, coupled model of nuclear waste contamination that includes thermal effects, proving its stability and convergence.

Contribution

It introduces a novel combined finite volume and finite element scheme for nonlinear coupled convection-diffusion and Darcy equations with proven stability and error estimates.

Findings

Scheme is stable under certain conditions.

Convergence and error estimates are established.

Applicable to nonlinear, coupled thermal and contaminant transport models.

Abstract

We consider a problem of nuclear waste contamination. It takes into account the thermal effects. The temperature and the contaminant's concentration fulfill convection-diffusion-reaction equations. The velocity and the pressure in the flow satisfy the Darcy equation, with a viscosity depending on both concentration and temperature. The equations are nonlinear and strongly coupled. Using both finite volume and nonconforming finite element methods, we introduce a scheme adapted to this problem. We prove the stability and convergence of this scheme and give some error estimates.

Equations

$$\varphi R_i{\partial }_t{a}_i+div(a_i\mathbf{u})-div(\varphi D_c\nabla a_i)=s_i-s{a}_i-{\lambda }_i{R}_i\varphi a_i+\sum _{j=1,j\ne i}^{N_r-1}{k}_{ij}{\lambda }_j{R}_j\varphi a_j,$$

$$\varphi {\mathcal{C}}_p{\partial }_t\theta +div(\theta \mathbf{u})-div(\varphi D_{\theta }\nabla \theta )=-s_{\theta }-s(\theta -{\theta }_{\ast }),$$

$$\displaystyle{\hbox{div}}\,{\bf u}=s,\quad{\bf u}+\frac{K}{\mu\big{(}(a_{j})_{j=1}^{N_{r}-1},\theta\big{)}}\,\nabla p={\bf f}.$$

$$\mu (a,\theta )={\mu }_R(a)\exp \left(\frac{1}{\theta }-\frac{1}{{\theta }_{\ast }}\right)$$

$${\mu }_R(a)=\mu (0)(1+(M^{1/4}-1)a_1{)}^{-4}$$

$$c_1=a_1,c_i=a_i+\sum _{j=1}^{i-1}\left(\prod _{l=j}^{i-1}\frac{y_{l+1}{\lambda }_l}{{\lambda }_l-{\lambda }_i}\right)a_j\mathrm{for}i=2,..,N_r.$$

$${\partial }_t{c}_i+div(c_i\mathbf{u})-D_c\Delta c_i=s_{c_i}-s{c}_i-{\lambda }_i{c}_i,$$

$${\partial }_t\theta +div(\theta \mathbf{u})-D_{\theta }\Delta \theta =-s_{\theta }-s(\theta -{\theta }_{\ast }),$$

$$div\mathbf{u}=s,\mathbf{u}+\kappa (c,\theta )\nabla p=\mathbf{f},$$

$$\nabla c_i\cdot \nu =0,\nabla \theta \cdot \nu =0,\mathbf{u}\cdot \mathbf{n}=0,$$

$$c_i{|}_{t=0}=c_0^i,\theta {|}_{t=0}={\theta }_0.$$

$$\kappa \in W^{1,\infty }((0,1)\times (0,\infty )),s,s_c,s_{\theta }\in L^2,\mathbf{f}\in \mathcal{C}(0,T;{\mathbf{L}}^2).$$

$$0\le c_0(\mathbf{x})\le 1,{\theta }_{-}\le {\theta }_0(\mathbf{x})\le {\theta }_{+}$$

$$2s(\mathbf{x})+\lambda \ge s_c(\mathbf{x})\ge 0,2({\theta }^{-}-{\theta }^{\ast })s(\mathbf{x})+s_{\theta }(\mathbf{x})\le 0,2({\theta }^{+}-{\theta }^{\ast })s(\mathbf{x})+s_{\theta }(\mathbf{x})\ge 0.$$

$$\frac{d}{dt}(\varphi R_j{\hat{a}}_j)=-\varphi ({\lambda }_j{R}_j+k^2{D}_c-ik\mathbf{u}){\hat{a}}_j+y_i{R}_{j-1}\varphi {\hat{a}}_{j-1}=-{\lambda }_j^{\prime }{R}_j\varphi {\hat{a}}_j+y_i{R}_{j-1}\varphi {\hat{a}}_{j-1}$$

$$\forall \sigma \in {\mathcal{E}}_h,d_{\sigma }\ge C|\sigma |\mathrm{and}|\sigma |\ge Ch.$$

$$\forall \sigma \in {\mathcal{E}}_h^{int},{\tau }_{\sigma }=|\sigma |/d_{\sigma }\ge C.$$

$$P_0=\{q\in L^2;\forall K\in {\mathcal{T}}_h,q{|}_{K}is\; a\; constant\}.$$

$$P_1^d$$

$$P_1^c$$

$$P_1^{nc}$$

$$H_d^1=\{q\in L^2;\forall K\in {\mathcal{T}}_h,q{|}_K\in H^1(K)\}.$$

$$\forall q_h\in H_d^1,\forall K\in {\mathcal{T}}_h,{\tilde{\nabla }}_h{q}_h{|}_K=\nabla q_h{|}_K$$

$$\forall q_h\in H_d^1,\Vert q_h{\Vert }_{1,h}=(|q_h{|}^2+|{\tilde{\nabla }}_h{q}_h{|}^2{)}^{1/2}.$$

$$\|q_{h}\|_{h}=\Bigl{(}\sum_{\sigma\in{\cal{E}}^{int}_{h}}\tau_{\sigma}\,(q_{L_{\sigma}}-q_{K_{\sigma}})^{2}\Bigr{)}^{1/2}\quad\mathrm{and}\quad\|q_{h}\|_{-1,h}=\sup_{\psi_{h}\in P_{0}}\frac{(q_{h},\psi_{h})}{\|\psi_{h}\|_{h}}$$

$$\mathbf{R}{\mathbf{T}}_0^{\mathbf{d}}$$

$$\mathbf{R}{\mathbf{T}}_0$$

$$\forall\,q\in H^{1}_{d}\,,\hskip 28.45274pt\forall\,\phi_{h}\in P^{c}_{1}\,,\hskip 28.45274pt\big{(}\nabla(\Pi_{P^{c}_{1}}q),\nabla\phi_{h}\big{)}=(\nabla q,\nabla\phi_{h}).$$

$$({\Pi }_{P_0}p{)}_K=\frac{1}{|K|}{\int }_{K}p(\mathbf{x})d\mathbf{x},({\tilde{\Pi }}_{P_0}q{)}_K=q({\mathbf{x}}_K),$$

$$\forall {\psi }_h\in P_1^{nc},({\tilde{\Pi }}_{P_1^{nc}}p,{\psi }_h)=(p,{\psi }_h);\forall \sigma \in {\mathcal{E}}_h,{\int }_{\sigma }({\Pi }_{P_1^{nc}}q)d\sigma ={\int }_{\sigma }qd\sigma .$$

$$\forall {\mathbf{w}}_h\in \mathbf{R}{\mathbf{T}}_0,({\tilde{\Pi }}_{\mathbf{R}{\mathbf{T}}_0}\mathbf{v},{\mathbf{w}}_h)=(\mathbf{v},{\mathbf{w}}_h);\forall \sigma \in {\mathcal{E}}_h^{int},{\int }_{\sigma }(\mathbf{w}-{\Pi }_{\mathbf{R}{\mathbf{T}}_0}\mathbf{w})\cdot {\mathbf{n}}_{K_{\sigma },\sigma }d\sigma =0.$$

$$|q-{\Pi }_{P_0}q|\le Ch\Vert q{\Vert }_1,|\mathbf{v}-{\Pi }_{\mathbf{R}{\mathbf{T}}_0}\mathbf{v}|\le Ch\Vert \mathbf{v}{\Vert }_1.$$

$$|p-{\Pi }_{P_1^{nc}}p|\le Ch\Vert p{\Vert }_1,|{\tilde{\nabla }}_h(q-{\Pi }_{P_1^{nc}}q)|\le Ch\Vert q{\Vert }_2.$$

$$|q-{\Pi }_{P_1^c}q|\le Ch\Vert q{\Vert }_{1,h}.$$

Full text

1 Introduction

A part of the high-level nuclear waste is now stored in environmentally safe locations. One has to consider the eventuality of a leakage through the engineering and geological barriers. It may cause the contamination of underground water sources far away from the original repository’s location. In the present paper, we consider such a problem of nuclear waste contamination in the basement. We take into account the thermal effects. The evolution in time of the temperature and of the contaminants concentration is then governed by convection-diffusion-reaction equations. The velocity and the pressure in the flow satisfy the Darcy equation, with a viscosity depending on both concentrations and temperature in a nonlinear way. The velocity satisfies an incompressibility constraint. We introduce a scheme adapted to this problem. We use both finite volume and nonconforming finite element methods. It ensures that a maximum principle holds and that the associated linear systems have good-conditioned matrices. We prove the stability and convergence of the scheme and give some error estimates.

Let us briefly point out some previous works. A complete model coupling concentrations and pressure equations is very seldom studied, since the system is strongly coupled. Instead, each equation is considered separately. In the general context of convection-diffusion-reaction equations, numerous schemes are available (see ? or ? and the references therein). Finite difference schemes are sometimes used for the convective term (in ? for instance). But they are not adapted for the complex geometry of a reservoir. More recently, finite volume methods were developed and analysed. Let us just cite the book ? or ? and the references therein. Finite elements (for the diffusive term) and finite volumes (for the convective term) are coupled for instance in ?? . In convection dominant problems, the equations are of degenerate parabolic type. This setting is considered in ?? . The reaction terms are specifically studied through operator splitting methods in ? . Now, in the specific context of porous media flow, we mention ? who consider only the evolution of the pressure. In ?? a more complete set of equations is used, and a mixed finite element approximation is developed. We stress that in all theses works, as in most, the thermic effects are not taken into account.

The present paper is organized as follows. Section 2 is devoted to the derivation of the model. In section 3, we introduce the discrete tools used in this paper. It allows us to define the numerical scheme of section 5. The analysis of the scheme uses the properties of section 4. We then prove the stability and convergence of the scheme, in sections 6 and 7 respectively. We conclude with some error estimates in section 8.

2 Model of contamination

The thickness of the medium is significantly smaller than its length and width. Hence it is reasonable to average the medium properties vertically and to describe the far-field repository by a polyhedral domain $\Omega$ of $\mathbb{R}^{2}$ with a smooth boundary $\partial\Omega$. It is characterized by a porosity $\phi$ and a permeability tensor $K$. The time interval of interest is $[0,T]$. We denote by $p$ the pressure, by $(a_{i})_{i=1}^{N_{r}}$ the concentrations of the $N_{r}$ radionuclides involved in the flow and by $\theta$ the temperature. The Darcy velocity is represented by ${\bf u}$. We assume a miscible and incompressible displacement. Due to the mass and energy conservation principles, the flow is governed by the following system satisfied in $\Omega\times[0,T]$, with $i=1,..,N_{r}-1$ (see ?).

[TABLE]

In (2.1) the retardation factors $R_{i}>0$ are due to the sorption mechanism. The real $\lambda_{i}^{-1}>0$ denotes the half life time of radionuclide $i$. The term $-\lambda_{i}\,R_{i}\,\phi\,a_{i}$ describes the radioactive decay of the i–th specy. Meanwhile, the quantity $\sum_{j\neq i}k_{ij}\,\lambda_{j}\,R_{j}\,\phi\,a_{j}$ is created by radioactive filiation. The molecular diffusion effects are given by the coefficient $D_{c}>0$. The contamination is represented by the source term $s_{i}$ and $s=\sum_{i=1}^{N_{r}}s_{i}$. In (2.2) the coefficient $\mathcal{C}_{p}>0$ is the relative specific heat of the porous medium. The thermic diffusion coefficient is denoted by $D_{\theta}>0$. The real $\theta_{*}>0$ is a reference temperature. The constitutive relation (2.3) is the Darcy law and ${\bf f}$ is a density of body forces. For a large range of temperatures $\mu$ has the form

[TABLE]

where $\mu_{R}$ is a nonlinear function. For instance, in the Koval model for a two-species mixture ?, we have

[TABLE]

where $M=\mu(0)/\mu(1)$ is the mobility ratio.

We notice that the equations (2.1)-(2.3) are strongly coupled. Moreover, every concentration equation (2.1) involves a different time scale. Therefore, it is difficult to build a numerical scheme that captures all the physical phenomena. We have to transform these equations. We first assume that only serial or parallel first-order reactions occur, so that $k_{ij}=y_{i}$ with $y_{1}=0$. Next, following ?, we assume that no two isotopes have identical decay rates and we set

[TABLE]

Lastly, without losing any mathematical difficulty (see remark 2.2 below), we set $R_{i}=1$ for $i=1,..,N_{r}-1$ and $\phi=1$, ${\mathcal{C}}_{p}=1$. We also set s_{c_{i}}=s_{i}+\sum_{j=1}^{i-1}\bigl{(}\prod_{l=j}^{i-1}\frac{y_{l+1}\lambda_{l}}{\lambda_{l}-\lambda_{i}}\bigr{)}\,s_{j} and \kappa(c,\theta)=K/\mu\big{(}(a_{i})_{i=1}^{N_{r}-1},\theta\big{)}. The contamination problem is now modelized by the following parabolic-elliptic system

[TABLE]

with $i=1,..,N_{r}-1$. These equations are completed with the boundary and initial conditions

[TABLE]

The pressure $p$ is normalized by $\int_{\Omega}p\,d{\bf x}=0$. Equations (2.5) are all similar. Thus, for the sake of simplicity, we will assume that there is only one. We set $N_{r}=2$ and $c=c_{1}$, $c_{0}=c^{1}_{0}$, $s_{c}=s_{c_{1}}$, $\lambda=\lambda_{1}$. The results of this paper readily extend to the general case.

We conclude with some notations and hypothesis. Let $D$ be a bounded open set of $\mathbb{R}^{k}$ with $k\geq 1$. We denote by ${\mathcal{C}}^{\infty}_{0}(D)$ the set of functions that are continuous on $D$ together with all their derivatives, and have a compact support in $D$. For $p\in\{2,\infty\}$, we use the Lebesgue spaces \big{(}L^{p}(D),\|.\|_{L^{p}(D)}\big{)} and \big{(}{\bf L}^{p}(D),\|.\|_{{\bf L}^{p}(D)}\big{)} with ${\bf L}^{p}=(L^{p})^{2}$. We also use the Sobolev spaces $W^{p,q}(D)$ for $p\in[1,\infty[$ and $q\in[1,\infty[$. In the case $D=\Omega$ we use the following conventions. We drop the domain dependancy. We denote by $|.|$ (resp. $\|.\|_{\infty}$) the norms associated to $L^{2}=L^{2}(\Omega)$ and ${\bf L}^{2}={\bf L}^{2}(\Omega)$ (resp. $L^{\infty}$ and ${\bf L}^{\infty}$). We set $L^{2}_{0}=\{v\in L^{2};\ \int_{\Omega}v({\bf x})\,d{\bf x}=0\}$. For $p\in[1,\infty[$ we define $(H^{p},\|.\|_{p})$ and $({\bf H}^{p},\|.\|_{p})$ with $H^{p}=W^{p,2}$ and ${\bf H}^{p}=(H^{p})^{2}$. Now let $(X,|.|)$ be a Banach space. The functions $g:[0,T]\to X$ such that $t\to\|g(t)\|_{X}$ is continuous (resp. bounded and square integrable) form the set ${\mathcal{C}}(0,T;X)$ (resp. $L^{\infty}(0,T;X)$ and $L^{2}(0,T;X)$). The associated norm for the space $L^{\infty}(0,T;X)$ (resp. $L^{2}(0,T;X)$) is defined by $\|g\|_{L^{\infty}(0,T;X)}=\sup_{t\in[0,T]}\|g(t)\|_{X}$ (resp. $\|g\|_{L^{2}(0,T;X)}=\left(\int_{0}^{T}\|f(t)\|^{2}_{X}\,dt\right)^{1/2}$). Finally in all computations we use $C>0$ as a generic constant. It depends only on the data of the problem.

We assume the following regularities for the data in (2.5)–(2.7)

[TABLE]

We also assume that $\kappa\geq\kappa_{inf}$ with $\kappa_{inf}>0$. For the initial data, we assume that $c_{0}\in H^{1}$, $\theta_{0}\in H^{1}$, and that we have a.e. in $\Omega$

[TABLE]

with $\theta_{-}>0$. Finally we assume that we have a.e. in $\Omega$

[TABLE]

These conditions ensure a maximum principle (proposition 6.1 below).

Remark 2.1

We have assumed that first-order reactions occur, so that the coefficients $k_{ij}$ in (2.1) depend only on $i$. If $k_{ij}$ depends on $i$ and $j$, one can still uncouple the equations by iterating the transformation (2.4), provided that $k_{1j}=0$ for $j=2,..,N_{r}-1$. This assumption means that the first long-lasting isotope disappears and is not created anymore. It is satisfied by many radionuclides.

Remark 2.2

We have assumed that the retardation factors $R_{j}$ are identical. If it is not the case, the difficulty and the approach remain the same. Indeed, let us consider the Fourier transform of (2.1). For a Fourier mode $\hat{a}_{j}(k,t)$ we obtain

[TABLE]

with $\lambda^{\prime}_{j}=\lambda_{j}+(k^{2}\,D_{c}-ik\,{\bf u})/R_{j}$ for $j=1,..,N_{r}-1$. A transform analogous to (2.4) uncouples the problem. By taking the partial differential equation counterpart, we obtain an equation similar to (2.5).

3 Discrete tools

3.1 Mesh and discrete spaces

Let ${\cal{T}}_{h}$ be a triangular mesh of $\Omega$. The circumscribed circle of a triangle $K\in{\mathcal{T}}_{h}$ is centered at ${\bf x}_{K}$ and has the diameter $h_{K}$. We set $h=\max_{K\in{\mathcal{T}}_{h}}h_{K}$. We assume that all the interior angles of the triangles of the mesh are less than $\frac{\pi}{2}$, so that ${\bf x}_{K}\in K$. The set of the edges of the triangle $K\in{\mathcal{T}}_{h}$ is ${\cal{E}}_{K}$. The symbol ${\mathbf{n}}_{K,\sigma}$ denotes the unit normal vector to an edge $\sigma\in{\cal{E}}_{K}$ and pointing outward $K$. We denote by ${\cal{E}}_{h}$ the set of the edges of the mesh. We distinguish the subset ${\mathcal{E}}^{int}_{h}\subset{\mathcal{E}}_{h}$ (resp. ${\mathcal{E}}^{ext}_{h}$) of the edges located inside $\Omega$ (resp. on $\partial\Omega$). The middle of an edge $\sigma\in{\cal{E}}_{h}$ is ${\bf x}_{\sigma}$ and its length is $|\sigma|$. For each edge $\sigma\in{\cal{E}}^{int}_{h}$ let $K_{\sigma}$ and $L_{\sigma}$ be the two triangles having $\sigma$ in common; we set $d_{\sigma}=d({\bf x}_{K_{\sigma}},{\bf x}_{L_{\sigma}})$. For all $\sigma\in{\mathcal{E}}^{ext}_{h}$ only the triangle $K_{\sigma}$ located inside $\Omega$ is defined and we set $d_{\sigma}=d({\bf x}_{K_{\sigma}},{\bf x}_{\sigma})$. Then for all $\sigma\in{\mathcal{E}}_{h}$ we set $\tau_{\sigma}=\frac{|\sigma|}{d_{\sigma}}$. We assume the following on the mesh (see ?): there exists $C>0$ such that

[TABLE]

It implies that there exists $C>0$ such that

[TABLE]

We define on the mesh the following spaces. The usual space for finite volume schemes is

[TABLE]

For any function $q_{h}\in P_{0}$ and any $K\in{\cal{T}}_{h}$ we set $q_{K}=q_{h}|_{K}$. We also consider

[TABLE]

We have $P_{1}^{c}\subset H^{1}$. On the other hand $P^{nc}_{1}\not\subset H^{1}$, but $P^{nc}_{1}\subset H_{d}^{1}$ with

[TABLE]

Thus we define $\widetilde{\nabla}_{h}:H^{1}_{d}\to{\bf L}^{2}$ by setting

[TABLE]

and the associated norm $\|.\|_{1,h}$ is given by

[TABLE]

We then have the following Poincaré-like inequality for the space $P^{nc}_{1}\cap L^{2}_{0}$ (see ?).

Proposition 3.1

There exists $C>0$ such that $|q_{h}|\leq C\,|\widetilde{\nabla}_{h}q_{h}|$ for all $q_{h}\in P^{nc}_{1}\cap L^{2}_{0}$.

We also define discrete analogues of the norms $H^{1}$ and $H^{-1}$ for the space $P_{0}$ by setting

[TABLE]

for any function $q_{h}\in P_{0}$. Note that for any $p_{h}\in P_{0}$ and $q_{h}\in P_{0}$, $(p_{h},q_{h})\leq\|q_{h}\|_{-1,h}\,\|q_{h}\|_{h}$. The following Poincaré-like inequality holds for the space $P_{0}\cap L^{2}_{0}$ (see ?).

Proposition 3.2

There exists $C>0$ such that $|q_{h}|\leq C\,\|q_{h}\|_{h}$ for all $q_{h}\in P_{0}\cap L^{2}_{0}$.

Finally we set ${\bf P}_{0}=(P_{0})^{2}$, ${\bf P}_{1}^{d}=(P_{1}^{d})^{2}$ and use the Raviart-Thomas spaces ?

[TABLE]

For all ${\bf v}_{h}\in\mathbf{RT_{0}}$, $K\in{\cal{T}}_{h}$ and $\sigma\in{\cal{E}}_{K}$, we set $({\bf v}_{h}\cdot{\mathbf{n}}_{K,\sigma})_{\sigma}={\bf v}_{h}|_{K}\cdot{\mathbf{n}}_{K,\sigma}$.

3.2 Projection operators

We associate with the spaces of section 3.1 some projection operators. First, we define $\Pi_{P_{1}^{c}}:H^{1}_{d}\to P_{1}^{c}$ by setting

[TABLE]

Next, we consider the space $P_{0}$. Let ${\mathcal{C}}_{d}=\{q_{h}\in L^{2}\,;\hskip 8.5359ptq_{h}\hbox{ is equal a.e. to a continuous function}\}$. We define $\Pi_{P_{0}}:L^{2}\to P_{0}$ and $\widetilde{\Pi}_{P_{0}}:{\mathcal{C}}_{d}\to P_{0}$ by setting

[TABLE]

for all $p\in L^{2}$, $q\in{\mathcal{C}}_{d}$ and $K\in{\mathcal{T}}_{h}$. We also set $\Pi_{{\bf P}_{0}}=(\Pi_{P_{0}})^{2}$. For the space $P^{nc}_{1}$, we define $\widetilde{\Pi}_{P^{nc}_{1}}:L^{2}\to P^{nc}_{1}$ and $\Pi_{P^{nc}_{1}}:H^{1}\to P^{nc}_{1}$. For all $p\in L^{2}$ and $q\in H^{1}$, $\widetilde{\Pi}_{P^{nc}_{1}}p$ and $\Pi_{P^{nc}_{1}}q$ satisfy

[TABLE]

For the space $\mathbf{RT_{0}}$, we define $\widetilde{\Pi}_{\mathbf{RT_{0}}}:{\bf L}^{2}\to\mathbf{RT_{0}}$ and $\Pi_{\mathbf{RT_{0}}}:{\bf H}^{1}\to\mathbf{RT_{0}}$. For all ${\bf v}\in{\bf L}^{2}$ and ${\bf w}\in{\bf H}^{1}$, $\widetilde{\Pi}_{\mathbf{RT_{0}}}{\bf v}$ and $\Pi_{\mathbf{RT_{0}}}{\bf w}$ satisfy

[TABLE]

The operators $\Pi_{P_{0}}$, $\widetilde{\Pi}_{P^{nc}_{1}}$ (resp. $\Pi_{{\bf P}_{0}}$, $\widetilde{\Pi}_{\mathbf{RT_{0}}}$) are $L^{2}$ (resp. ${\bf L}^{2}$) projection operators. They are stable for the $L^{2}$ (resp. ${\bf L}^{2}$) norms. The operators $\widetilde{\Pi}_{P_{0}}$, $\Pi_{P^{nc}_{1}}$ and $\Pi_{\mathbf{RT_{0}}}$ are interpolation operators. The following estimates are classical (? p.109 and ?).

Proposition 3.3

There exists $C>0$ such that for all $q\in H^{1}$ and ${\bf v}\in{\bf H}^{1}$

[TABLE]

For all $p\in H^{1}$ and $q\in H^{2}$ we have

[TABLE]

For all $q\in H^{d}_{1}$ we have

[TABLE]

Finally, using the Sobolev embedding theorem, one checks that

[TABLE]

for all $q\in W^{1,r}$ with $r>1$ (see ?).

3.3 Discrete operators

Equations (2.5)–(2.7) use the differential operators gradient, divergence and laplacian. We have to define analogous operators in the discrete setting. The discrete gradient operator $\nabla_{h}:P^{nc}_{1}\to P_{0}$ is the restriction to $P^{nc}_{1}$ of the operator $\widetilde{\nabla}_{h}$ given by (3.2). The discrete divergence operator $\hbox{div}_{h}:{\bf P}_{0}\to P^{nc}_{1}$ is defined by

[TABLE]

for all ${\bf v}_{h}\in{\bf P}_{0}$. It is adjoint to $\nabla_{h}$ (proposition 4.1 below). The discrete laplacian operator $\Delta_{h}:P_{0}\to P_{0}$ is the usual one for finite volume schemes (see ?). For all $q_{h}\in P_{0}$ and $K\in{\mathcal{T}}_{h}$ we have

[TABLE]

Let us now consider the convection terms in (2.5) and (2.6). We define ${\widetilde{\hbox{\rm b}}}:{\bf H}^{1}\times H^{1}\to L^{2}$ by

[TABLE]

for all $q\in H^{1}$ and ${\bf v}\in{\bf H}^{1}$. In order to define a discrete counterpart to ${\widetilde{\hbox{\rm b}}}$ we use the classical upwind scheme (see ?). The discrete operator ${\widetilde{\hbox{\rm b}}}_{h}:\mathbf{RT_{0}}\times P_{0}\to P_{0}$ is such that

[TABLE]

for all ${\bf v}_{h}\in\mathbf{RT_{0}}$, $q_{h}\in P_{0}$ and $K\in{\mathcal{T}}_{h}$. We have set $a^{+}=\max(a,0)$ and $a^{-}=\min(a,0)$ for all $a\in\mathbb{R}$. Integrating by parts the convection terms also leads to consider ${\hbox{\rm b}}:{\bf L}^{2}\times L^{2}\times L^{\infty}\to\mathbb{R}$ defined by

[TABLE]

for all ${\bf v}\in L^{2}$, $p\in L^{2}$ and $q\in L^{\infty}$. The discrete counterpart is ${\hbox{\rm b}}_{h}:\mathbf{RT_{0}}\times P_{0}\times P_{0}\to\mathbb{R}$ with

[TABLE]

for all ${\bf v}_{h}\in\mathbf{RT_{0}}$, $p_{h}\in P_{0}$ and $q_{h}\in P_{0}$.

4 Properties of the discrete operators

The properties of the discrete operators are analogous to the ones satisfied by their continuous counterpart. The gradient and divergence operators are adjoint. For the operators $\nabla_{h}$ and $\hbox{div}_{h}$ we state in ? the following.

Proposition 4.1

For all ${\bf v}_{h}\in{\bf P}_{0}$ and $q_{h}\in P^{nc}_{1}$ we have $({\bf v}_{h},\nabla_{h}q_{h})=-(q_{h},\hbox{\rm div}_{h}{\bf v}_{h})$.

Let us now consider the convection terms. Let $q\in L^{\infty}\cap H^{1}$, ${\bf v}\in{\bf L}^{2}$ with $\hbox{div}\,{\bf v}\in L^{2}$ and $\hbox{div}\,{\bf v}({\bf x})\geq 0$ a.e. in $\Omega$. We obtain ${\hbox{\rm b}}({\bf v},q,q)=\int_{\Omega}(q^{2}/2)\,\hbox{div}\,{\bf v}\,d{\bf x}\geq 0$ by integration by parts. For ${\hbox{\rm b}}_{h}$ we state in ? a similar result.

Proposition 4.2

Let ${\bf v}_{h}\in\mathbf{RT_{0}}$ with $\hbox{\rm div}\,{\bf v}_{h}\geq 0$. We have ${\hbox{\rm b}}_{h}({\bf v}_{h},q_{h},q_{h})\geq 0$ for all $q_{h}\in P_{0}$ .

The following stability properties are used to prove the error estimates in section 8.

Proposition 4.3

There exists $C>0$ such that for all $p_{h}\in P_{0}$, $q_{h}\in P_{0}$ and ${\bf v}_{h}\in\mathbf{RT_{0}}$ with $\hbox{ \rm div}\,{\bf v}_{h}=0$

[TABLE]

There exists $C>0$ such that for all $p_{h}\in P_{0}$, $q_{h}\in P_{0}\cap L^{2}_{0}$, ${\bf v}_{h}\in\mathbf{RT_{0}}$

[TABLE]

Proof. For all $K\in{\cal{T}}_{h}$ and $\sigma\in{\cal{E}}_{K}\cap{\cal{E}}^{int}_{h}$ we write

[TABLE]

Thus (3.12) reads ${\hbox{\rm b}}_{h}({\bf v}_{h},p_{h},q_{h})=S_{1}+S_{2}$ with

[TABLE]

Using the Cauchy-Schwarz inequality we write

[TABLE]

Since ${\bf v}_{h}\in\mathbf{RT_{0}}\subset(P^{d}_{1})^{2}$ we have $h\,\|{\bf v}_{h}\|_{\infty}\leq C\,|{\bf v}_{h}|$ (? p. 112). Moreover (3.1) implies $\sum_{\sigma\in{\cal{E}}^{int}_{h}}(p_{L_{\sigma}}-p_{K_{\sigma}})^{2}\leq C\,\sum_{\sigma\in{\cal{E}}^{int}_{h}}\tau_{\sigma}\,(p_{L_{\sigma}}-p_{K_{\sigma}})^{2}=C\,\|p_{h}\|^{2}_{h}$ and $\sum_{\sigma\in{\cal{E}}^{int}_{h}}(q_{L_{\sigma}}-q_{K_{\sigma}})^{2}\leq C\,\|q_{h}\|^{2}_{h}$. Thus

[TABLE]

We now consider $S_{2}$. We have $({\bf v}_{h}\cdot{\mathbf{n}}_{K,\sigma})_{\sigma}=0$ for all $K\in{\cal{T}}_{h}$ and $\sigma\in{\cal{E}}_{K}\cap{\cal{E}}^{ext}_{h}$. Thus we write

[TABLE]

It gives the following relation.

[TABLE]

Thus if $\hbox{div}\,{\bf v}_{h}=0$ then $S_{2}=0$ and estimate (4.3) gives (4.1). Let us prove (4.2). Since $q_{h}\in P_{0}\cap L^{2}_{0}$ we can apply proposition 3.2. Using the Cauchy-Schwarz inequality we get

[TABLE]

This latter estimate together with (4.3) gives (4.2).

Lastly, we claim that ${\widetilde{\hbox{\rm b}}}_{h}$ is a consistent approximation of ${\widetilde{\hbox{\rm b}}}$ ?.

Proposition 4.4

Let $r>0$. There exists $C>0$ such that for all functions $q\in H^{2}$ and ${\bf v}\in{\bf H}^{1+r}$ with ${\bf v}\cdot{\mathbf{n}}|_{\partial\Omega}=0$

[TABLE]

Let us now consider the discrete laplacian. We have a coercivity and stability result.

Proposition 4.5

For all $p_{h}\in P_{0}$ and $q_{h}\in P_{0}$, we have

[TABLE]

Proof. Definition (3.8) implies

[TABLE]

Setting $q_{h}=p_{h}$ gives the first part of the result. Using the Cauchy-Schwarz inequality, we get the second one.

We also deduce from (4.4) the following property.

Proposition 4.6

For all $p_{h}\in P_{0}$ and $q_{h}\in P_{0}$ we have $(\Delta_{h}p_{h},q_{h})=(p_{h},\Delta_{h}q_{h})$.

Lastly, we state that $\Delta_{h}$ is a consistent approximation of the laplacian. The proof follows the lines of the one of proposition 1.14 in ?.

Proposition 4.7

There exists $C>0$ such that for all $q\in H^{2}$ with $\nabla q\cdot{\mathbf{n}}|_{\partial\Omega}=0$ we have

[TABLE]

5 The finite element-finite volume scheme

We now introduce the scheme for (2.5)-(2.9). The interval $[0,T]$ is split with a constant time step $k=T/N$. We set $[0,T]=\bigcup^{N-1}_{m=0}[t_{m},t_{m+1}]$ with $t_{m}=m\,k$. The time derivatives are approximated using a first order Euler scheme. The convection terms are discretized semi-implicitly in time and the other ones in an implicit way. We set $s_{h}=\Pi_{P_{0}}s$, $s_{h}^{c}=\Pi_{P_{0}}s_{c}$, $s_{h}^{\theta}=\Pi_{P_{0}}s_{\theta}$ and ${\bf f}^{m}_{h}=\Pi_{{\bf P}_{0}}{\bf f}(t_{m})$ for all $m\in\{0,\dots,N\}$. Since $\Pi_{P_{0}}$ (resp. $\Pi_{{\bf P}_{0}}$) is stable for the $L^{2}$ (resp. ${\bf L}^{2}$) norm we have

[TABLE]

The initial values are $c^{0}_{h}=\Pi_{P_{0}}c_{0}$ and $\theta^{0}_{h}=\Pi_{P_{0}}\theta_{0}$. Then for all $n\in\{0,\dots,N-1\}$, the quantities $c^{n+1}_{h}\in P_{0}$, $\theta^{n+1}_{h}\in P_{0}$, $p^{n+1}_{h}\in P^{nc}_{1}\cap L^{2}_{0}$, ${\bf u}^{n+1}_{h}\in\mathbf{RT_{0}}$ are the solutions of the following problem.

[TABLE]

with $\kappa^{n+1}_{h}=\kappa(c^{n+1}_{h},\theta^{n+1}_{h})\in P_{0}$. This term is defined thanks to proposition 6.1 below. Note also that the boundary conditions are implicitly included in the definition of the discrete operators (section 3.3). The existence of a unique solution to (5.2) and (5.3) is classical (see ?). Since $\kappa^{n+1}_{h}\geq\kappa_{min}>0$ and $p^{n+1}_{h}\in L^{2}_{0}$ equation (5.4) also has a unique solution (see ?). We have a discrete equivalent for the divergence condition (2.7).

Proposition 5.1

For all $m\in\{1,\dots,N\}$ we have $\hbox{\rm div}\,{\bf u}^{m}_{h}=s_{h}$ .

Proof. Let $m\in\{1,\dots,N\}$ and $n=m-1$. We compare the solution of (5.4)–(5.5) with the solution of the following mixed hybrid problem. Let ${\cal{E}}_{0}=\left\{\mu_{h}:\cup_{\sigma\in{\cal{E}}_{h}}\to\mathbb{R}\,;\;\forall\,\sigma\in{\cal{E}}_{h}\,,\;\mu_{h}|_{\sigma}\hbox{ is constant}\right\}$. Then ${\tilde{\mathbf{u}}}^{m}_{h}\in\mathbf{RT^{d}_{0}}$, $\overline{p}^{m}_{h}\in P_{0}$ and $\lambda^{m}_{h}\in{\cal{E}}_{0}$ are the solution of (see ?)

[TABLE]

and $\widetilde{p}^{m}_{h}\in P^{nc}_{1}$ is defined by $\int_{\sigma}\widetilde{p}^{m}_{h}\,d\sigma=\lambda^{m}_{\sigma}$ for all $\sigma\in{\cal{E}}_{h}$. Let $\sigma\in{\cal{E}}_{h}$. We define $\phi_{\sigma}\in P^{nc}_{1}$ by setting $\phi_{\sigma}({\bf x}_{\sigma})=1$ and $\phi_{\sigma}({\bf x}_{\sigma^{\prime}})=0$ for all $\sigma^{\prime}\in{\cal{E}}_{h}\backslash\{\sigma\}$. We set ${\bf v}_{h}=\nabla_{h}\phi_{\sigma}\in{\bf P}_{0}\subset\mathbf{RT^{d}_{0}}$ in (5.6). We have

[TABLE]

and according to the gradient formula

[TABLE]

Thus we get from (5.6)

[TABLE]

The first term in (5.8) is treated as follows. Integrating by parts we get

[TABLE]

Since (5.7) implies that ${\tilde{\mathbf{u}}}^{m}_{h}\in\mathbf{RT_{0}}$, we have

[TABLE]

Thus $({\tilde{\mathbf{u}}}^{m}_{h},\nabla_{h}\phi_{\sigma})=-(\phi_{\sigma},\hbox{div}\,{\tilde{\mathbf{u}}}^{m}_{h})$. Then, using (5.7), we get

[TABLE]

Furthermore, according to proposition 4.1, we have (\kappa^{m}_{h}\nabla_{h}\widetilde{p}^{m}_{h},\nabla_{h}\phi_{\sigma})=-\big{(}\phi_{\sigma},\hbox{div}_{h}(\kappa^{m}_{h}\,\nabla_{h}\widetilde{p}^{m}_{h})\big{)} and $({\bf f}^{m}_{h},\nabla_{h}\phi_{\sigma})=-(\phi_{\sigma},\hbox{div}_{h}{\bf f}^{m}_{h})$. Hence we deduce from (5.8) that

[TABLE]

Since $(\phi_{\sigma})_{\sigma\in{\cal{E}}_{h}}$ is a basis of $P^{nc}_{1}$, we get $\hbox{div}_{h}(\kappa^{m}_{h}\nabla_{h}\widetilde{p}^{m}_{h})=\hbox{div}_{h}{\bf f}^{m}_{h}-\widetilde{\Pi}_{P^{nc}_{1}}s_{h}^{m}$. Thus, by (5.4), there exists a real $C$ such that $\widetilde{p}^{m}_{h}=p^{m}_{h}+C$. We now compare ${\tilde{\mathbf{u}}}^{m}_{h}$ with ${\bf u}^{m}_{h}$. Since for all ${\bf v}_{h}\in\mathbf{RT_{0}}$ we have

[TABLE]

it follows from (5.6) that $({\tilde{\mathbf{u}}}^{m}_{h},{\bf v})=({\bf f}^{m}_{h}-\kappa^{m}_{h}\nabla_{h}\widetilde{p}^{m}_{h},{\bf v})$ for any ${\bf v}\in\mathbf{RT_{0}}$. It means that

[TABLE]

Thus ${\bf u}^{m}_{h}={\tilde{\mathbf{u}}}^{m}_{h}$ satisfies (5.7) and $\hbox{div}\,{\bf u}^{m}_{h}=s_{h}$.

6 Stability analysis

We first check that a maximum principle holds. It ensures that the computed concentration and temperature are physically relevant.

Proposition 6.1

For any $m\in\{0,\dots,N\}$ we have $0\leq c^{m}_{h}\leq 1$ and $\theta_{-}\leq\theta^{m}_{h}\leq\theta_{+}$.

Proof. We prove the result by induction. Since $c^{0}_{h}=\Pi_{P_{0}}c_{0}$ and $\theta^{0}_{h}=\Pi_{P_{0}}\theta_{0}$ the result holds for $m=0$ thanks to (2.11) and (3.4). Let us assume that it is true for $m=n\in\{0,\dots,N-1\}$. Let $K\in{\cal{T}}_{h}$. Equation (5.2) implies

[TABLE]

We consider the last term of this relation. Since for any $\sigma\in{\cal{E}}_{K}\cap{\cal{E}}^{int}_{h}$ we have

[TABLE]

We deduce from (3.10)

[TABLE]

Since ${\bf u}^{n}_{h}\in\mathbf{RT_{0}}$, $({\bf u}^{n}_{h}\cdot{\mathbf{n}}_{K,\sigma})_{\sigma}=0$ for any $\sigma\in{\cal{E}}_{K}\cap{\cal{E}}^{ext}_{h}$. It implies that $\sum_{\sigma\in{\cal{E}}_{K}\cap{\cal{E}}^{int}_{h}}|\sigma|\,({\bf u}^{n}_{h}\cdot{\mathbf{n}}_{K,\sigma})_{\sigma}=\sum_{\sigma\in{\cal{E}}_{K}}|\sigma|\,({\bf u}^{n}_{h}\cdot{\mathbf{n}}_{K,\sigma})_{\sigma}$. Thus using the divergence formula and proposition 5.1 we obtain

[TABLE]

Therefore we get

[TABLE]

We consider $K_{i}\in{\cal{T}}_{h}$ such that $c^{n+1}_{K_{i}}=\min_{K\in{\cal{T}}_{h}}c^{n+1}_{K}$. According to hypothesis (2.12) and definition (3.4) we have $2\,s_{K_{i}}+\lambda\geq 0$ and $s^{c}_{K_{i}}\geq 0$. Thus, using the induction hypothesis, we deduce from (6.1)

[TABLE]

We now consider $K_{s}\in{\cal{T}}_{h}$ such that $c^{n+1}_{K_{s}}=\max_{K\in{\cal{T}}_{h}}c^{n+1}_{K}$. Using again hypothesis (2.12) we have $2\,s_{K_{s}}+\lambda\geq s^{c}_{K_{s}}\geq 0$. Thus, using the induction hypothesis, we deduce from (6.1)

[TABLE]

A similar work for equation (5.3) proves that $\theta_{-}\leq\min_{K\in{\cal{T}}_{h}}\theta^{n+1}_{K}$ and $\max_{K\in{\cal{T}}_{h}}\theta^{n+1}_{K}\leq\theta_{+}$. Thus the induction hypothesis still holds for $m=n+1$.

We now state the stability of the scheme (5.2)-(5.5).

Proposition 6.2

For any $1\leq m\leq N$ we have

[TABLE]

Proof. Let $0\leq n\leq N-1$. Multiplying (5.2) by $2\,k\,c^{n+1}_{h}$ we get

[TABLE]

We have $(c^{n+1}_{h}-c^{n}_{h},2\,c^{n+1}_{h})=|c^{n+1}_{h}|^{2}-|c^{n}_{h}|^{2}+|c^{n+1}_{h}-c^{n}_{h}|^{2}$. Thanks to propositions 4.2 and 4.5

[TABLE]

Using the Cauchy-Schwarz and Young inequalities we write

[TABLE]

Finally thanks to (2.12) and (3.4) we have $s_{h}\geq 0$. Thus we obtain

[TABLE]

Let $m\in\{1,\dots,N\}$. Summing up the latter relation from $n=0$ to $m-1$ we get

[TABLE]

thanks to proposition 6.1. With a similar work on equation (5.3), we get (6.2). We now prove (LABEL:eq:estud). Let $n=m-1\in\{0,\dots,N-1\}$. Multiplying equation (5.4) by $-p^{n+1}_{h}$ and using proposition 4.1, we get

[TABLE]

The left-hand side term satisfies $(\kappa^{n+1}_{h}\nabla_{h}p^{n+1}_{h},\nabla_{h}p^{n+1}_{h})\geq\kappa_{inf}\,|\nabla_{h}p^{n+1}_{h}|^{2}$. We now consider the right-hand side. Using (5.1), the Cauchy-Schwarz and Young inequalities we write

[TABLE]

Also, the stability of $\widetilde{\Pi}_{P^{nc}_{1}}$ for the $L^{2}$-norm, proposition 3.1 and the Young inequality lead to

[TABLE]

Thus we deduce from (6.4) that $|\nabla_{h}p^{n+1}_{h}|^{2}=|\nabla_{h}p^{m}_{h}|^{2}\leq C$. Then (5.5) imply

[TABLE]

Estimate (LABEL:eq:estud) is proven.

7 Convergence analysis

Let $\varepsilon=\max(h,k)$. In this section we study the behavior of the scheme (5.2)-(5.5) as $\varepsilon\to 0$. We first define the applications $c_{\varepsilon}:\mathbb{R}\to P_{0}$, $\tilde{c}_{\varepsilon}:\mathbb{R}\to P_{0}$, $\theta_{\varepsilon}:\mathbb{R}\to P_{0}$, $p_{\varepsilon}:\mathbb{R}\to P^{nc}_{1}$, $s_{\varepsilon}:\mathbb{R}\to P_{0}$, $s^{c}_{\varepsilon}:\mathbb{R}\to P_{0}$ and ${\bf u}_{\varepsilon}:\mathbb{R}\to\mathbf{RT_{0}}$, ${\bf f}_{\varepsilon}:\mathbb{R}\to{\bf P}_{0}$ by setting for all $n\in\{0,\dots,N-1\}$ and $t\in[t_{n},t_{n+1}]$

[TABLE]

and for all $t\not\in(0,T)$

[TABLE]

We recall that the Fourier transform $\widehat{f}$ of a function $f\in L^{1}(\mathbb{R})$ is defined for any $\tau\in\mathbb{R}$ by

[TABLE]

We begin with the following estimate.

Proposition 7.1

Let $0<\gamma<\frac{1}{4}$. There exists $C>0$ such that for all $\varepsilon>0$

[TABLE]

Proof. Since equations (5.2) and (5.3) are similar we only prove the estimate on $\widehat{c}_{\varepsilon}$. We first define $g_{\varepsilon}:\mathbb{R}\to P_{0}\cap L^{2}_{0}$ as the solution of

[TABLE]

Multiplying this equation by $-g_{\varepsilon}$ we obtain

[TABLE]

Proposition 4.5 allows us to write

[TABLE]

Thanks to the Cauchy-Schwarz inequality, (5.1) and proposition 3.2 we have

[TABLE]

According to proposition 4.3, then proposition 6.1 and (LABEL:eq:estud), we have

[TABLE]

Let us plug these estimates into (7.2) and integrate from 0 to $T$. We get

[TABLE]

because of proposition 5.1 and (6.2). Definition (7.1) then leads to

[TABLE]

We now use this estimate to prove (7.1). Equation (5.2) reads

[TABLE]

where $\delta_{0}$ and $\delta_{T}$ are Dirac distributions respectively localized in [math] and $T$. Let $\tau\in\mathbb{R}$. Applying the Fourier transform to the latter equation we obtain

[TABLE]

Let us take the scalar product of this relation with $i\hbox{sign}(\tau)\,\widehat{\tilde{c}}_{\varepsilon}(\tau)$. Applying propositions 3.2 and 4.5 leads to

[TABLE]

We assume that $\tau\neq 0$ and multiply this estimate by $|\tau|^{2\,\gamma-1}$. Using proposition 6.1 and (7.3) we get

[TABLE]

Using the Young inequality and integrating over $\{\tau\in\mathbb{R}\,;|\tau|>1\}$, we obtain

[TABLE]

For $|\tau|\leq 1$, we have $|\tau|^{2\,\gamma}\,|\widehat{\tilde{c}}_{\varepsilon}(\tau)|^{2}\leq|\widehat{\tilde{c}}_{\varepsilon}(\tau)|^{2}\leq C\,\|\widehat{\tilde{c}}_{\varepsilon}(\tau)\|^{2}_{h}$ according to proposition 3.2. Thus

[TABLE]

By combining the bounds for $|\tau|>1$ and $|\tau|\leq 1$ we get

[TABLE]

Since $4\,\gamma-2<-1$, we have $\int_{|\tau|>1}|\tau|^{4\,\gamma-2}\,d\tau\leq C$. Thanks to the Parseval theorem and (6.2)

[TABLE]

because $\|c^{0}_{h}\|_{h}=\|\Pi_{P_{0}}c_{0}\|_{h}\leq C\,\|c_{0}\|_{1}$ (see ? p. 776). Hence the result.

We can now prove the following convergence result.

Proposition 7.2

There exists a subsequence of $(c_{\varepsilon},\theta_{\varepsilon},p_{\varepsilon},{\bf u}_{\varepsilon})_{\varepsilon>0}$, not relabeled for convenience, such that the following convergences hold for $\varepsilon\to 0$

[TABLE]

The limits $(c,\theta,p,{\bf u})$ satisfy the following properties. We have $c\in L^{2}(0,T;H^{1})$, $\theta\in L^{2}(0,T;H^{1})$, $p\in L^{\infty}(0,T;H^{1})$ and ${\bf u}\in L^{\infty}(0,T;{\bf L}^{2})$. We also have $0\leq c({\bf x},t)\leq 1$ and $\theta^{-}\leq\theta({\bf x},t)\leq\theta^{+}$ a.e. in $\Omega\times[0,T]$. For all $\phi\in{\mathcal{C}}^{\infty}_{0}(\Omega\times(-1,T))$, $c$ and $\theta$ satisfy

[TABLE]

Lastly we have

[TABLE]

Proof. In what follows, the convergence results hold for extracted subsequences. They are not relabeled for convenience. We begin by proving (7.4). According to proposition 6.1, the sequence $(c_{\varepsilon})_{\varepsilon>0}$ is uniformly bounded in $L^{\infty}(0,T;L^{2})$. Thus there exists $c\in L^{\infty}(0,T;L^{2})$ such that

[TABLE]

Using the Fourier transform, we prove that this convergence is strong. Let $d_{\varepsilon}=c_{\varepsilon}-c$ and $M>0$. We use the following splitting

[TABLE]

Since $|\widehat{d}_{\varepsilon}(\tau)|^{2}\leq 2|\widehat{c}_{\varepsilon}(\tau)|^{2}+2|\widehat{c}(\tau)|^{2}$ we have

[TABLE]

Using proposition 7.1 we write

[TABLE]

Hence

[TABLE]

This implies that for all $\varepsilon>0$, $I^{M}_{\varepsilon}\to 0$ when $M\to\infty$. We now consider $J^{M}_{\varepsilon}$. Let $\tau\in[-M,M]$. Since $c_{\varepsilon}(t)\in P_{0}$ for all $t\in\mathbb{R}$, and $c_{\varepsilon}\rightharpoonup c$ weakly in $L^{2}(0,T;L^{2})$, we deduce from (7.1) that $\widehat{c}_{\varepsilon}(\tau)\in P_{0}$ and $\widehat{c}_{\varepsilon}(\tau)\rightharpoonup\widehat{c}(\tau)$ weakly in $L^{2}$. Extanding $\widehat{c}_{\varepsilon}(\tau)$ by [math] outside $\Omega$, one checks (? p.811) that

[TABLE]

Then, using estimate (6.2), we deduce from ? (p.834) that $\widehat{c}_{\varepsilon}(\tau)\to\widehat{c}(\tau)$ strongly in $L^{2}$. Thus $\widehat{d}_{\varepsilon}(\tau)=\widehat{c}_{\varepsilon}(\tau)-\widehat{c}(\tau)\to 0$ in $L^{2}$, so that $J^{\varepsilon}_{M}\to 0$ when $\varepsilon\to 0$. Now, let us report the limits for $I^{\varepsilon}_{M}$ and $J^{\varepsilon}_{M}$ into (7.9). Using the Parseval identity we get

[TABLE]

Thus we have proven that $c_{\varepsilon}\to c$ in $L^{2}(0,T;L^{2})$. A similar work proves that $\theta_{\varepsilon}\to\theta$ in $L^{2}(0,T;L^{2})$ with $\theta\in L^{\infty}(0,T;L^{2})$. Hence (7.4) is proven. Moreover using proposition 6.1 we obtain $0\leq c({\bf x},t)\leq 1$ and $\theta^{-}\leq\theta({\bf x},t)\leq\theta^{+}$ a.e. in $\Omega\times[0,T]$. Lastly, using (6.2) and (7.9), we get as in ? (p.811) that $c\in L^{2}(0,T;H^{1})$ and $\theta\in L^{2}(0,T;H^{1})$.

Let us now consider the sequences $(p_{\varepsilon})_{\varepsilon>0}$ and $({\bf u}_{\varepsilon})_{\varepsilon>0}$. According to (3.3) and (LABEL:eq:estud) the sequence $(\Pi_{P^{c}_{1}}p_{\varepsilon})_{\varepsilon>0}$ is bounded in $L^{\infty}(0,T;H^{1})$. It implies that there exists $p\in L^{\infty}(0,T;H^{1})$ such that $\Pi_{P^{c}_{1}}p_{\varepsilon}\rightharpoonup p$ weakly in $L^{2}(0,T;H^{1})$. Using proposition 3.3 we get $p_{\varepsilon}\rightharpoonup p$ weakly in $L^{2}(0,T;H^{1})$. Moreover, according to (LABEL:eq:estud), the sequence $({\bf u}_{\varepsilon})_{\varepsilon>0}$ is bounded in $L^{\infty}(0,T;{\bf L}^{2})$ . Thus we have ${\bf u}_{\varepsilon}\rightharpoonup{\bf u}$ weakly in $L^{2}(0,T;{\bf L}^{2})$ with ${\bf u}\in L^{\infty}(0,T;{\bf L}^{2})$. We check the properties of ${\bf u}$. Using a Taylor expansion, the Cauchy-Schwarz inequality, and a density argument, we have

[TABLE]

Thus, using the strong convergence of the sequences $(c_{\varepsilon})_{\varepsilon>0}$ and $(\theta_{\varepsilon})_{\varepsilon>0}$, we have $\kappa(c_{\varepsilon},\theta_{\varepsilon})\to\kappa(c,\theta)$ in $L^{2}(0,T;L^{2})$. Since $\nabla_{h}p_{\varepsilon}\rightharpoonup\nabla p$ weakly in $L^{2}(0,T;{\bf L}^{2})$, we deduce from this

[TABLE]

Now let ${\bf v}\in L^{2}(0,T;({\mathcal{C}}^{\infty}_{0})^{2})$. According to (5.5) we have

[TABLE]

Using proposition 3.3 one checks easily that $({\bf f}_{\varepsilon},\Pi_{\mathbf{RT_{0}}}{\bf v})\to({\bf f},{\bf v})$ and $({\bf u}_{\varepsilon},\Pi_{\mathbf{RT_{0}}}{\bf v})\to({\bf u},{\bf v})$ in $L^{1}(0,T)$. Using moreover convergence (7.10) and a density argument, we deduce from this that ${\bf u}={\bf f}-\kappa(c,\theta)\,\nabla p$. And since $\hbox{div}\,{\bf u}_{\varepsilon}=s_{\varepsilon}$ by proposition 5.1, we also have $\hbox{div}\,{\bf u}=s$.

We finally prove that $c$ satisfies (7.6). For all $t\in(0,T)$ equation (5.2) reads

[TABLE]

Let $\psi\in{\cal{C}}^{\infty}_{0}(\Omega\times(-1,T))$ and $\psi_{h}=\widetilde{\Pi}_{P_{0}}\psi$. Multiplying the latter equation by $\psi_{h}$ and integrating over $[0,T]$, we obtain

[TABLE]

We now pass to the limit $\varepsilon\to 0$ in this equation. We begin with the term $\int^{T}_{0}{\hbox{\rm b}}_{h}({\bf u}_{\varepsilon},c_{\varepsilon},\psi_{h})\,dt$. We use the splitting ${\hbox{\rm b}}({\bf u},c,\psi)-{\hbox{\rm b}}_{h}({\bf u}_{\varepsilon},c_{\varepsilon},\psi_{h})=A^{\varepsilon}_{1}+A^{\varepsilon}_{2}+A^{\varepsilon}_{3}$ with

[TABLE]

According to definition (3.11)

[TABLE]

We know that $c\,\nabla\psi\in L^{2}(0,T;{\bf L}^{2})$. Since ${\bf u}_{\varepsilon}\rightharpoonup{\bf u}$ in $L^{2}(0,T;{\bf L}^{2})$ we get $\int^{T}_{0}A^{\varepsilon}_{1}\,dt\to 0$. We now consider $A^{\varepsilon}_{2}$. We have

[TABLE]

Using the Cauchy-Schwarz inequality we get

[TABLE]

Using a Taylor expansion, one checks that $\|\psi-\psi_{h}\|_{L^{\infty}(\Omega\times(0,T))}\leq h\,\|\nabla\psi\|_{L^{\infty}(\Omega\times(0,T))}$. Thus $\int^{T}_{0}A^{\varepsilon}_{2}\,dt\to 0$. Finally we estimate $A^{\varepsilon}_{3}$. For all triangles $K\in{\mathcal{T}}_{h}$ and $L\in{\mathcal{T}}_{h}$ sharing an edge $\sigma$, we set $c_{K,L}=c_{K}$ if ${\bf u}_{\varepsilon}\cdot{\mathbf{n}}_{K,\sigma}\geq 0$ and $c_{K,L}=c_{L}$ otherwise. Using the divergence formula, we deduce from definition (3.12)

[TABLE]

Using definition (3.5) this reads

[TABLE]

Using a Taylor expansion, one checks that $|\psi_{K_{\sigma}}-\psi_{L_{\sigma}}|\leq h\,\|\nabla\psi\|_{L^{\infty}(\Omega\times(0,T))}$. Moreover $|\sigma|\leq h$. Thus, using the Cauchy-Schwarz inequality, we have

[TABLE]

Using the assumption on the mesh, one checks that $|K|\geq C\,h^{2}$ for all $K\in{\mathcal{T}}_{h}$. Thus, thanks to a quadrature formula, we have

[TABLE]

We write $\Pi_{P^{nc}_{1}}c-c_{\varepsilon}=(\Pi_{P^{nc}_{1}}c-c)+(c-c_{\varepsilon})$ and we use proposition 3.3 . We obtain with (LABEL:eq:estud)

[TABLE]

Since $c_{\varepsilon}\to c$ in $L^{2}(0,T;L^{2})$ when $\varepsilon=\max(h,k)\to 0$, we conclude that $\int^{T}_{0}A^{\varepsilon}_{3}\,dt\to 0$. Gathering the limits for $A^{\varepsilon}_{1}$, $A^{\varepsilon}_{2}$, $A^{\varepsilon}_{3}$, we obtain

[TABLE]

We now consider the other terms in (7.11). Proposition 4.6 leads to

[TABLE]

According to proposition 4.7

[TABLE]

We then apply the Cauchy-Schwarz inequality and use (6.2). We obtain

[TABLE]

Moreover, since $c_{\varepsilon}\to c$ in $L^{2}(0,T;L^{2})$, we have $\int^{T}_{0}(c_{\varepsilon},\Delta\psi)\,dt\to\int^{T}_{0}(c,\Delta\psi)\,dt$. Thus we deduce from (7.12)

[TABLE]

We are left with two terms. First, using Taylor expansions, one checks that

[TABLE]

We know that $c_{\varepsilon}\to c$ in $L^{2}(0,T;L^{2})$. Thus

[TABLE]

Finally, integrating by parts the first term of (7.11), we get

[TABLE]

Since $\psi\in{\cal{C}}^{\infty}_{0}(\Omega\times(-1,,T))$ we have $(\tilde{c}_{\varepsilon},\psi_{h})_{t=T}=0$. Using proposition 3.3 one checks that $c^{0}_{h}=\Pi_{P_{0}}c_{0}\to c_{0}$ in $L^{2}$; using moreover (7.13) we get

[TABLE]

For the last term, one easily checks that $\|\tilde{c}_{\varepsilon}-c_{\varepsilon}\|_{L^{2}(0,T;L^{2})}\to 0$. Thus, since $c_{\varepsilon}\to c$ in $L^{2}(0,T;L^{2})$, we also have $\tilde{c}_{\varepsilon}\to c$ in $L^{2}(0,T;L^{2})$. Using moreover (7.13) we get $\int^{T}_{0}(\tilde{c}_{\varepsilon},\partial_{t}\psi_{h})\,dt\to\int^{T}_{0}(c,\partial_{t}\psi)\,dt$. Therefore

[TABLE]

By gathering the limits we have obtained in (7.11) we get (7.6). The relation (LABEL:eq:propconv4) for $\theta$ is proven in a similar way.

8 Error estimates

We have proven in section 7 that the problem (2.5)–(2.9) has a weak solution $(c,\theta,p,{\bf u})$. From now on, we assume the following regularity for this solution:

[TABLE]

with $r>0$ and $s>0$. We also assume that ${\bf f}\in{\mathcal{C}}(0,T;{\bf H}^{1})$.

8.1 Definitions

For all $m\in\{0,\dots,N\}$, we define the following errors

[TABLE]

We have the following splittings

[TABLE]

with the discrete errors

[TABLE]

and the interpolation errors

[TABLE]

The interpolation errors are estimated as follows. We write $|\eta^{m}_{h,c}|\leq|c(t_{m})-\Pi_{P_{0}}c(t_{m})|+|\Pi_{P_{0}}c(t_{m})-\widetilde{\Pi}_{P_{0}}c(t_{m})|$ and the same for $\eta^{m}_{h,\theta}$. Using proposition 3.3 and (3.7) we obtain

[TABLE]

According to proposition 3.3 we also have

[TABLE]

We now have to estimate the discrete errors.

Proposition 8.1

For all $n\in\{0,\dots,N-1\}$ and $\psi_{h}\in P^{nc}_{1}$ we have

[TABLE]

For all $m\in\{0,\dots,N\}$, the consistency errors $C^{m}_{h,1}$, $C^{m}_{h,2}$, $\Theta^{m}_{h,1}$, $\Theta^{m}_{h,2}$, ${\bf P}^{m}_{h}$ and $\mathbf{U}^{m}_{h}$ are defined in (8.9), (8.10), (8.14), (8.15) and the terms $\kappa^{m}_{h,1}$ and $\kappa^{m}_{h,2}$ are given by (8.13) below.

Proof. Let $n\in\{0,\dots,N-1\}$. Equation (2.5) for $t=t_{n+1}$ reads

[TABLE]

We introduce the time discretization by setting

[TABLE]

We get

[TABLE]

We apply $\Pi_{P_{0}}$ to this equation. By subtracting the result from (5.2) we get

[TABLE]

We now introduce the discrete errors as follows. Since $c(t_{n+1})-c(t_{n})=\int^{t_{n+1}}_{t_{n}}c_{t}(s)\,ds$ one checks that

[TABLE]

We also have

[TABLE]

Using the linearity of ${\widetilde{\hbox{\rm b}}}_{h}$, one easily checks that

[TABLE]

Lastly

[TABLE]

Using these relations in (8.8) we get (8.4). For any $m\in\{1,\dots,N\}$, the consistency errors $C^{m}_{h,1}\in P_{0}$ and $C^{m}_{h,2}\in P_{0}$ are given by

[TABLE]

A similar proof leads to (8.5) where the consistence errors $\Theta^{m}_{h,1}\in P_{0}$ and $\Theta^{m}_{h,2}\in P_{0}$ are defined for any $m\in\{1,\dots,N\}$ by

[TABLE]

We now consider the problem associated with the pressure. Let $n\in\{0,\dots,N-1\}$ and $\psi_{h}\in P^{nc}_{1}$. Multiplying equation (2.7) written for $t=t_{n+1}$ by $\psi_{h}$ and integrating by parts, we get

[TABLE]

On the other hand, using (5.4) and proposition 4.1, we have

[TABLE]

Since $\nabla_{h}\psi_{h}\in P_{0}$, one checks that $({\bf f}^{n+1}_{h},\nabla_{h}\psi_{h})=(\Pi_{{\bf P}_{0}}{\bf f}(t_{n+1}),\nabla_{h}\psi_{h})=({\bf f}(t_{n+1}),\nabla_{h}\psi_{h})$. According to (3.5) we also have $(\widetilde{\Pi}_{P^{nc}_{1}}s_{h},\psi_{h})=(s_{h},\psi_{h})$. Thus

[TABLE]

Substracting (8.11) from the latter relation, we obtain

[TABLE]

We split the left-hand side as

[TABLE]

Using a Taylor expansion, one can check that

[TABLE]

We have set for any $m\in\{0,\dots,N\}$ and $s\in[0,1]$

[TABLE]

We also have

[TABLE]

Plugging these relations into (8.12) we get (8.6). For all $m\in\{0,\dots,N\}$ we have

[TABLE]

We end with the equation associated with ${\bf u}$. Let $n\in\{0,\dots,N-1\}$. Applying the operator $\widetilde{\Pi}_{\mathbf{RT_{0}}}$ to (2.7) for $t=t_{n+1}$ we obtain

[TABLE]

Let us substract this equation from (5.5). Since ${\bf f}_{h}^{n+1}=\Pi_{{\bf P}_{0}}{\bf f}(t_{n+1})$ we get

[TABLE]

One easily checks that

[TABLE]

Thus we get (8.7). For all $m\in\{0,\dots,N\}$, we have

[TABLE]

This ends the proof of proposition 8.1.

8.2 Error estimates

We first estimate the consistency errors.

Proposition 8.2

For all $m\in\{1,\dots,N\}$ the consistency errors satisfy

[TABLE]

Proof. Let $n\in\{1,\dots,N\}$. Since the operator $\Pi_{P_{0}}$ is stable for the $L^{2}$-norm we have

[TABLE]

Using a Taylor expansion and the Cauchy-Schwarz inequality, we get

[TABLE]

On the other hand, since $\nabla c(t_{n})|_{\partial\Omega}=0$, we deduce from (3.10) by integrating by parts

[TABLE]

Using a Taylor expansion and the Cauchy-Schwarz inequality, we get

[TABLE]

Thus

[TABLE]

Thanks to the stability of $\Pi_{P_{0}}$ for the $L^{2}$-norm and to (8.1) we have

[TABLE]

The Cauchy-Schwarz inequality and (3.7) allow us to write

[TABLE]

By plugging these estimates into definition (8.9) we get

[TABLE]

Summing up the latter relation for $n=1$ to $m\in\{1,\dots,N\}$ and using a similar work on $\Theta^{m}_{h,1}$ we get (8.16). Now let $n\in\{1,\dots,N\}$. Using propositions 4.7 and 4.4 we have

[TABLE]

and

[TABLE]

Plugging these estimates into definition (8.1) and summing up from $n=1$ to $m$, we obtain

[TABLE]

A similar work on $\Theta^{m}_{h,2}$ then leads to (8.17). We finally prove (8.18). Let $m\in\{1,\dots,N\}$. On the one hand, we have by (8.14)

[TABLE]

Using estimates (8.1)–(LABEL:eq:esterrintu) we get

[TABLE]

On the other hand definition (8.15) leads to

[TABLE]

Using the stability of $\widetilde{\Pi}_{\mathbf{RT_{0}}}$ for the ${\bf L}^{2}$-norm and proposition 3.3 we have

[TABLE]

Using moreover (LABEL:eq:esterrintu) we obtain

[TABLE]

We have proven (8.18).

Using the former proposition we are now able to estimate the discrete errors.

Proposition 8.3

There exists some real $k_{0}>0$ such that for any $k<k_{0}$ and $m\in\{1,\dots,N\}$

[TABLE]

Proof. Multiplying (8.4) by $2\,k\,\varepsilon^{n+1}_{h,c}$, we obtain

[TABLE]

Using an algebraic identity we have

[TABLE]

We know by propositions 4.2 and 4.5 that

[TABLE]

We have

[TABLE]

Using the Young inequality, we also write

[TABLE]

We are left with the term {\hbox{\rm b}}_{h}\big{(}\hbox{\boldmath\epsilon^{\hbox{\unboldmath$\scriptstyle n$}}{\hbox{\unboldmath$\scriptstyle h,{\bf u}$}} \unboldmath},\widetilde{\Pi}_{P_{0}}c(t_{n+1}),\varepsilon^{n+1}_{h,c}\big{)}. We have \hbox{\boldmath\epsilon^{\hbox{\unboldmath$\scriptstyle n$}}{\hbox{\unboldmath$\scriptstyle h,{\bf u}$}} \unboldmath}\in\mathbf{RT_{0}}. Using the divergence formula, proposition 5.1 and (3.6), one easily checks that \hbox{div}\,\hbox{\boldmath\epsilon^{\hbox{\unboldmath$\scriptstyle n$}}_{\hbox{\unboldmath$\scriptstyle h,{\bf u}$}} \unboldmath}=0. Thus we can apply proposition 4.3 to get

[TABLE]

Let us first bound $\|\widetilde{\Pi}_{P_{0}}c(t_{n+1})\|_{h}$. We have

[TABLE]

Using an inverse inequality (see proposition 1.2 in ?) and (3.7)

[TABLE]

Moreover, according to ? (p. 776), we have $\|\Pi_{P_{0}}c(t_{n+1})\|_{h}\leq C\,\|c\|_{L^{\infty}(0,T;H^{1})}$. Thus $\|\widetilde{\Pi}_{P_{0}}c(t_{n+1})\|_{h}\leq C$. We now estimate |\hbox{\boldmath\epsilon^{\hbox{\unboldmath$\scriptstyle n$}}_{\hbox{\unboldmath$\scriptstyle h,{\bf u}$}} \unboldmath}|. Using the stability of $\widetilde{\Pi}_{\mathbf{RT_{0}}}$ for the ${\bf L}^{2}$-norm and the Cauchy-Schwarz inequality in (8.7) we get

[TABLE]

We bound $\varepsilon^{n+1}_{h,p}$ as follows. Setting $\psi_{h}=\varepsilon^{n+1}_{h,p}$ in (8.6) and using the Cauchy-Schwarz inequality, we get

[TABLE]

The left-hand side is such that $(\kappa(c^{n+1}_{h},\theta^{n+1}_{h})\nabla_{h}\varepsilon^{n+1}_{h,p},\nabla_{h}\varepsilon^{n+1}_{h,p})\geq\kappa_{inf}\,|\nabla_{h}\varepsilon^{n+1}_{h,p}|^{2}$. As for the right-hand side, we have $s_{h}=\Pi_{P_{0}}s$ and $\varepsilon^{n+1}_{h,p}\in P^{nc}_{1}\cap L^{2}_{0}$. Thus, according to propositions 3.1 and 3.3

[TABLE]

Finally $|{\bf P}^{n+1}_{h}|\leq C\,h$ thanks to (8.18). Therefore we obtain

[TABLE]

Let us plug this estimate into (8.23). Since $|\mathbf{U}^{n+1}_{h}|\leq C\,h$ thanks to (8.18), we get

[TABLE]

Now, plugging this bound into (8.23) and using the Young inequality, we obtain

[TABLE]

Now we have treated all the terms in (8.2). This equation implies

[TABLE]

Let $m\in\{1,\dots,N\}$. Let us sum up the latter estimate from $n=0$ to $m-1$. Thanks to (3.3)

[TABLE]

Using moreover estimates (8.16) and (8.17) we get

[TABLE]

Summing this relation with the one obtained by a similar work on (8.5) we obtain

[TABLE]

Using a discrete Gronwall lemma (see lemma 5.2 in ?) we get (8.19). Then (8.24) and (8.25) imply (8.20).

By combining proposition 8.3 with estimates (8.1)-(LABEL:eq:esterrintu), we obtain finally the following result.

Theorem 8.1

There exists a real $k_{0}>0$ such that for all $k<k_{0}$ and $m\in\{1,\dots,N\}$

[TABLE]

\bibname

• Achdou, Y., Bernardi, C. & Coquel, F. (2003) A priori and a posteriori analysis of finite volume discretizations of Darcy’s equations, Numer. Math., 96(1) 17–42.
• Afif, M.,Amaziane, B. (2002) Convergence of finite volume schemes for a degenerate convection-diffusion equation arising in flow in porous media, Comput. Methods Appl., 191 5265–5286.
• Angot, P., Doljši, V., Feistauer, M. & Felcman, J. (1998) Analysis of a combined barycentric finite volume-nonconforming finite element method for nonlinear convection-diffusion problems, Appl. Math., 43(4) 263–310.
• Bateman, H. (1910) The solution of a system of differential equations occuring in the theory of radioactive transformations, Proc. Camb. Philos. Soc., 15 423–427.
• Bear, J. (1972) Dynamics of Fluids in Porous Media, New York: American Elsevier.
• Brenner, S. C., Scott, L. R. & Ridgway, L. (2002) The mathematical theory of finite element methods, Texts in applied mathematics, 16, New-York: Springer-Verlag.
• Brezzi, F. & Fortin, M. (1991) Mixed and hybrid finite element methods, Springer Series in Computational Mathematics, 15, New-York: Springer-Verlag.
• Chavent, G., Jaffré, J. & Roberts, J. E. (1995) Mixed hybrid finite elements and cell-centered finite volumes for two-phase flow in porous media. In: A. P. Bourgeat and al, ed. Mathematical modeling of flow through porous media, London: World Scientific, 100–114.
• Chen, Z. (2000) Formulation and numerical methods of the black oil model in porous media, SIAM J. Numer. Anal., 38(2) 489–514.
• Choquet C. (2004) Existence result for a radionuclide transport model with unbounded viscosity , J. Math. Fluid Mech., 6(4) 365–388.
• Douglas J., Jr. & Spagnuolo, A. M. (2001) The approximation of nuclear contaminant transport in porous media, J. Korean Math. Soc., 38 723–761.
• Ewing, R. E., Yuan, Y. & Li, G. (1989) A time-discretization procedure for a mixed finite element approximation of contamination by incompressible nuclear waste in porous media, In: Mathematics for large scale computing, Lecture Notes in Pure and Appl. Math., 120, 127–145.
• Eymard, R., Gallouët, T. & Herbin, R. (2000) The finite volume method (chapter III), In P.G. Ciarlet and J. L. Lions, eds Handbook of numerical analysis, 7, 762–835.
• Geiser, J. (2002) Numerical simulation of a model for transport and reaction of radionuclides, In: Algoritmy 2002, Vysoké Tatry - Podbanské, Slovakia, Bratislava: Slovak university of technology.
• John, V., Maubach, J. M. & Tobiska, L. (1997) Nonconforming streamline-diffusion finite-element methods for convection-diffusion problems, Numer. Math., 78(2) 165–188.
• Koval, E. J. (1963) A method for predicting the performance of unstable miscible displacements in heterogeneous media, SPEJ trans. AIME, 228 145–154.
• Lin P. & Yang, D. (1998) An iterative perturbation method for the pressure equation in the simulation of miscible displacement in porous media, SIAM J. Sci. Comput., 19(3) 893–911.
• Michel, A. (1999) Convergence of a finite volume scheme for a nonlinear convection-diffusion problem, In: Proceedings of the Second international symposium on finite volumes for complex applications, Duisburg.
• Morton, K. W. (1996) Numerical solution of convection-diffusion problems, In: Applied Mathematics and mathematical Computation, 12, London: Chapman and Hall.
• Verwer, J. G., Blom, J. G. & Hundsdorfer, W. (1996) An implicit-explicit approach for atmospheric transport chemistry problems, Appl. Numer. Math., 20(1-2) 191–209.
• Zimmermann, S. (2006) Stability of a finite volume scheme for incompressible fluids, preprint.
• Zimmermann, S. (2006) Étude et implémentation de méthodes de volumes finis pour les fluides incompressibles. Thesis (PhD). Blaise Pascal university (France).