Analysis of a fractional-step scheme for the P1 radiative diffusion model

TL;DR

This paper develops and analyzes a finite volume fractional-step scheme for the P1 radiative diffusion model, proving existence, uniqueness, and convergence of solutions with uniform bounds.

Contribution

It introduces a novel fractional-step scheme for the nonlinear radiative diffusion model and provides rigorous proofs of stability, existence, uniqueness, and convergence of the discrete solutions.

Findings

Discrete solutions satisfy a priori L-estimates

Existence and uniqueness of solutions are established

Convergence to the continuous problem is proven

Abstract

We address in this paper a nonlinear parabolic system, which is built to retain the main mathematical difficulties of the P1 radiative diffusion physical model. We propose a finite volume fractional-step scheme for this problem enjoying the following properties. First, we show that each discrete solution satisfies a priori L -estimates, through a discrete maxi- mum principle; by a topological degree argument, this yields the existence of a solution, which is proven to be unique. Second, we establish uniform (with respect to the size of the meshes and the time step) L2 -bounds for the space and time translates; this proves, by the Kolmogorov theorem, the relative compactness of any sequence of solutions obtained through a sequence of discretizations the time and space steps of which tend to zero; the limits of converging subsequences are then shown to be a solution to the continuous problem. Estimates of time translates of the discrete solutions are obtained through the formalization of a generic argument, interesting for its own sake.

Equations

$$\left|\begin{array}[]{ll}\displaystyle\frac{\partial u}{\partial t}-\Delta u+u^{4}-\varphi=0&\displaystyle\mbox{for }(x,t)\in\Omega\times(0,T),\\[4.30554pt] \displaystyle\varphi-\Delta\varphi-u^{4}=0&\displaystyle\mbox{for }(x,t)\in\Omega\times(0,T),\\[4.30554pt] \displaystyle u(x,0)=u_{0}(x)&\displaystyle\mbox{for }x\in\Omega,\\[4.30554pt] \displaystyle u(x,t)=0&\displaystyle\mbox{for }(x,t)\in\partial\Omega\times(0,T),\\[4.30554pt] \displaystyle\nabla\varphi\cdot n=0&\displaystyle\mbox{for }(x,t)\in\partial\Omega\times(0,T),\end{array}\right.$$

$${\theta }_{\mathcal{M}}:=\underset{K\in \mathcal{M}}{\min }\frac{{\rho }_K}{h_K}$$

$$[u,v{]}_{\mathcal{M}}:=\sum _{\sigma \in {\mathcal{E}}_{\mathrm{int}}(\sigma =K|L)}\frac{|\sigma |}{d_{\sigma }}(u_L-u_K)(v_L-v_K)+\sum _{\sigma \in {\mathcal{E}}_{\mathrm{ext}}(\sigma \in \mathcal{E}(K))}\frac{|\sigma |}{d_{\sigma }}{u}_K{v}_K.$$

$$\Vert u{\Vert }_{1,\mathcal{M}}:=[u,u{]}_{\mathcal{M}}^{1/2}.$$

$$\Vert u{\Vert }_{L^2(\Omega )}\le \mathrm{diam}(\Omega )\Vert u{\Vert }_{1,\mathcal{M}}\forall u\in {\mathrm{H}}_{\mathcal{M}}(\Omega ).$$

$$<u,v{>}_{\mathcal{M}}:=\sum _{\sigma \in {\mathcal{E}}_{\mathrm{int}}(\sigma =K|L)}\frac{|\sigma |}{d_{\sigma }}(u_L-u_K)(v_L-v_K),|u{|}_{1,\mathcal{M}}:=<u,u{>}_{\mathcal{M}}^{1/2}.$$

$$\Vert u{\Vert }_{-1,\mathcal{M}}:=\underset{\begin{array}{c}v\in {\mathrm{H}}_{\mathcal{M}}(\Omega )\\ v\ne 0\end{array}}{\sup }\frac{1}{\Vert v{\Vert }_{1,\mathcal{M}}}{\int }_{\Omega }uv\mathrm{d}x.$$

$$\begin{array}[]{l}\displaystyle\forall\psi\in{\rm H}_{\mathcal{M}}(\Omega),\\[8.61108pt] \qquad\begin{array}[]{l}\displaystyle(\Delta_{\mathcal{M},{\rm N}}(\psi))_{K}=\frac{1}{|K|}\sum_{\sigma=K|L}\frac{|\sigma|}{d_{\sigma}}\ (\psi_{L}-\psi_{K}),\\ \displaystyle(\Delta_{\mathcal{M},{\rm D}}(\psi))_{K}=(\Delta_{\mathcal{M},{\rm N}}(\psi))_{K}+\frac{1}{|K|}\sum_{\sigma\in\mathcal{E}_{{\rm ext}}\cap\mathcal{E}(K)}\frac{|\sigma|}{d_{\sigma}}\ (-\psi_{K}).\end{array}\end{array}$$

$$\begin{array}[]{ll}\displaystyle\forall\psi\in{\rm H}_{\mathcal{M}}(\Omega),&\displaystyle\sum_{K\in\mathcal{M}}-|K|\ \psi_{K}\ (\Delta_{\mathcal{M},{\rm N}}(\psi))_{K}=<\psi,\psi>_{\mathcal{M}}\\[17.22217pt] \displaystyle\mbox{and}&\displaystyle\sum_{K\in\mathcal{M}}-|K|\ \psi_{K}\ (\Delta_{\mathcal{M},{\rm D}}(\psi))_{K}=[\psi,\psi]_{\mathcal{M}}.\end{array}$$

$$\forall K\in \mathcal{M},$$

$$\frac{u_K^{n+1}-u_K^n}{\delta t}-({\Delta }_{\mathcal{M},\mathrm{D}}(u^{n+1}){)}_K+|u_K^{n+1}|(u_K^{n+1}{)}^3-{\phi }_K^n=0,\text{mbox}forn=0,\dots ,N-1,$$

$${\phi }_K^{n+1}-({\Delta }_{\mathcal{M},\mathrm{N}}({\phi }^{n+1}){)}_K-(u_K^{n+1}{)}^4=0\text{mbox}forn=0,\dots ,N-1,$$

$$u_K^0=\frac{1}{|K|}{\int }_K{u}_0(x)\mathrm{d}x,$$

$${\phi }_K^0-({\Delta }_{\mathcal{M},\mathrm{N}}({\phi }^0){)}_K-(u_K^0{)}^4=0$$

$$\forall K\in \mathcal{M},0\le {\phi }_K^n\le {\left[\underset{L\in \mathcal{M}}{\max }{u}_L^n\right]}^4.$$

$$\forall K\in \mathcal{M},0\le u_K^n\le \underset{L\in \mathcal{M}}{\max }{u}_L^{n-1}.$$

$${\phi }_K^n-({\Delta }_{\mathcal{M},\mathrm{N}}({\phi }^n){)}_K=(u_K^n{)}^4.$$

$$\left[\frac{1}{\delta t}+|u_K^{n+1}|(u_K^{n+1}{)}^2\right]u_K^{n+1}-({\Delta }_{\mathcal{M},\mathrm{D}}(u^{n+1}){)}_K=\frac{1}{\delta t}{u}_K^n+{\phi }_K^n.$$

$$({\phi }^{n+1}-\bar{\phi }{)}_K-({\Delta }_{\mathcal{M},\mathrm{N}}({\phi }^{n+1}-\bar{\phi }){)}_K=(u_K^{n+1}{)}^4-{\bar{\phi }}_K.$$

$${\phi }_K^{n+1}\le (\underset{L\in \mathcal{M}}{\max }{u}_L^{n+1}{)}^4,\forall K\in \mathcal{M}.$$

$$-({\Delta }_{\mathcal{M},\mathrm{D}}(u^{n+1}){)}_{K_0}\ge 0.$$

$$\frac{1}{\delta t}(u_{K_0}^{n+1}-u_{K_0}^n)-({\Delta }_{\mathcal{M},\mathrm{D}}(u^{n+1}){)}_{K_0}+\left[|u_{K_0}^{n+1}|(u_{K_0}^{n+1}{)}^3-{\phi }_{K_0}^n\right]=0.$$

$$\left|\begin{array}[]{l}\begin{array}[]{lcl}F:{\mathbb{R}}^{{\rm card}(\mathcal{M})}\times[0,1]&\longrightarrow&{\mathbb{R}}^{{\rm card}(\mathcal{M})},\\[4.30554pt] \displaystyle(U,\alpha)=\left((u_{K})_{K\in\mathcal{M}},\ \alpha\right)&\mapsto&F(U,\alpha)=\displaystyle(v_{K})_{K\in\mathcal{M}}\quad\mbox{ such that:}\end{array}\\ \displaystyle\hskip 30.1388pt\forall K\in\mathcal{M},\quad v_{K}=\frac{1}{\delta t}\,(u_{K}-u_{K}^{n})-(\Delta_{\mathcal{M},{\rm D}}(u))_{K}+\alpha\left[|u_{K}|\,(u_{K})^{3}-\varphi_{K}^{n}\right].\end{array}\right.$$

$$F((u_K^{n+1}{)}_{K\in \mathcal{M}},1)=0.$$

$$\forall K\in \mathcal{M},\left[\frac{1}{\delta t}+\left((u_K^{n+1}{)}^3+(u_K^{n+1}{)}^2{v}_K^{n+1}+u_K^{n+1}(v_K^{n+1}{)}^2+(v_K^{n+1}{)}^3\right)\right]\delta u_K-({\Delta }_{\mathcal{M},\mathrm{D}}(\delta u){)}_K=0.$$

$${\partial }_{t,\mathcal{D}}(u)(x,t)={\partial }_{t,\mathcal{D}}(u{)}^n(x)\text{mbox}for{t}^n\le t<t^{n+1},0\le n\le N-1$$

$$\partial_{t,\mathcal{D}}(u)^{n}(x):=\frac{u^{n+1}(x)-u^{n}(x)}{\delta t}\qquad\bigl{(}\mbox{\emph{i.e.\/}\ }\partial_{t,\mathcal{D}}(u)^{n}_{K}=\frac{u^{n+1}_{K}-u^{n}_{K}}{\delta t},\quad\forall K\in\mathcal{M}\bigr{)}.$$

$$\begin{array}[]{l}\displaystyle\|u\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}^{2}:=\delta t\ \sum_{n=0}^{N}\|u^{n}\|_{1,\mathcal{M}}^{2},\\[12.91663pt] \displaystyle\|u\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{-1})}^{2}:=\delta t\ \sum_{n=0}^{N-1}\|u^{n}\|_{-1,\mathcal{M}}^{2},\\[12.91663pt] \displaystyle|u|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}^{2}:=\delta t\ \sum_{n=0}^{N}|u^{n}|_{1,\mathcal{M}}^{2}.\end{array}$$

$$\begin{array}[]{l}\displaystyle\|u\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}+\|\partial_{t,\mathcal{D}}(u)\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{-1})}+\|\varphi\|_{L^{2}((0,T)\times\Omega)}\\[8.61108pt] \displaystyle\hfill+|\varphi|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}+\|\partial_{t,\mathcal{D}}(\varphi)\|_{L^{2}((0,T)\times\Omega)}\leq c_{\rm e},\end{array}$$

$$-\sum _{K\in \mathcal{M}}|K|u_K({\Delta }_{\mathcal{M},\mathrm{D}}(u){)}_K=\Vert u{\Vert }_{1,\mathcal{M}}^2,$$

$$-\sum _{K\in \mathcal{M}}|K|u_K({\Delta }_{\mathcal{M},\mathrm{N}}(u){)}_K=|u{|}_{1,\mathcal{M}}^2.$$

$$\begin{array}[]{l}\displaystyle\|u^{n+1}\|_{L^{2}(\Omega)}^{2}+\|u^{n+1}-u^{n}\|_{L^{2}(\Omega)}^{2}-\|u^{n}\|_{L^{2}(\Omega)}^{2}+2\,\delta t\,\|u^{n+1}\|_{1,\mathcal{M}}^{2}\\[8.61108pt] \displaystyle\hfill+2\,\delta t\,\int_{\Omega}(u^{n+1})^{5}\hskip 2.15277pt{\rm d}x=2\,\delta t\,\int_{\Omega}\varphi^{n}u^{n+1}\hskip 2.15277pt{\rm d}x.\end{array}$$

$$\Vert u^{n+1}{\Vert }_{L^2(\Omega )}^2-\Vert u^n{\Vert }_{L^2(\Omega )}^2+2\delta t\Vert u^{n+1}{\Vert }_{1,\mathcal{M}}^2\le 2\delta t{\bar{u}}_0^4{\int }_{\Omega }{u}^{n+1}\mathrm{d}x.$$

Full text

1. 1 Introduction
2. 2 The finite volume scheme
3. 3 A priori ${\rm L}^{\infty}$ estimates, existence and uniqueness of the discrete solution
4. 4 Convergence to a solution of the continuous problem
5. 5 Conclusion
6. A Estimation of time translates
7. B Regularity of the limit

1. Introduction

We address in this paper the following nonlinear parabolic system:

[TABLE]

where $\Omega$ is a connected bounded subset of ${\mathbb{R}}^{d}$, $d=2$ or $d=3$, which is supposed to be polygonal ($d=2$) or polyhedral ($d=3$), $T<\infty$ is the final time, $u$ and $\varphi$ are two real-valued functions defined on $\Omega\times[0,T)$ and $\partial\Omega$ stands for the boundary of $\Omega$ of outward normal $n$. The initial value for $u$, denoted by $u_{0}$, is supposed to lie in ${\rm L}^{\infty}(\Omega)\cap{\rm H}^{1}_{0}(\Omega)$ and to satisfy $u_{0}(x)\geq 0$ for almost every $x\in\Omega$.

This system of partial differential equations is inspired from a simplified radiative transfer physical model, the so-called P1 model, sometimes used in computational fluid dynamics for the simulation of high temperature optically thick flows, as encountered for instance in fire modelling (see e.g. [13, 10] for expositions of the theory or the documentation of the CFX or FLUENT commercial codes for a synthetic description).

In this context, the unknown $u$ stands for the temperature, $\varphi$ for the radiative intensity and the first equation is the energy balance. System (1) has been derived with the aim of retaining the main mathematical difficulties of the initial physical model; in particular, adding a convection term in the first equation would only require minor changes in the theory developed hereafter.

In this short paper, we give a finite volume scheme for the discretization of (1) and prove the existence and uniqueness of the discrete solution and its convergence to a solution of (1), thus showing that this problem indeed admits a solution, in a weak sense which will be defined.

2. The finite volume scheme

Even though the arguments developed in this paper are valid for any general admissible discretization in the sense of [6, Definition 9.1, p. 762], we choose, for the sake of simplicity, to restrict the presentation to simplicial meshes. We thus suppose given a triangulation $\mathcal{M}$ of $\Omega$, that is a finite collection of d-simplicial control volumes $K$, pairwise disjoint, and such that $\displaystyle\bar{\Omega}=\cup_{K\in\mathcal{M}}\bar{K}$; the mesh is supposed to be conforming in the sense that two neighbouring simplices share a whole face (i.e. there is no hanging node). In addition, we assume that, for any $K\in\mathcal{M}$, the circumcenter $x_{K}$ of $K$ lies in $K$; note that, for each neighbouring control volumes $K$ and $L$, the segment $[x_{K},x_{L}]$ is orthogonal to the face $K|L$ separating $K$ from $L$.

For each simplex $K$, we denote by $\mathcal{E}(K)$ the set of the faces of $K$ and by $|K|$ the measure of $K$. The set of faces of the mesh $\mathcal{E}$ is split into the set $\mathcal{E}_{{\rm int}}$ of internal faces (i.e. separating two control volumes) and the set $\mathcal{E}_{{\rm ext}}$ of faces included in the domain boundary. For each internal face $\sigma=K|L$, we denote by $|\sigma|$ the (d-1)-dimensional measure of $\sigma$ and by $d_{\sigma}$ the distance $d(x_{K},x_{L})$; for an external face $\sigma$ of a control volume $K$, $d_{\sigma}$ stands for distance from $x_{K}$ to $\sigma$. The regularity of the mesh is characterized by the parameter $\theta_{\mathcal{M}}$ defined by:

[TABLE]

where $\rho_{K}$ and $h_{K}$ stands for the diameter of the largest ball included in $K$ and the diameter of $K$, respectively.

We denote by ${\rm H}_{\mathcal{M}}(\Omega)\subset{\rm L}^{2}(\Omega)$ the space of functions which are piecewise constant over each control volume $K\in\mathcal{M}$. For all $u\in{\rm H}_{\mathcal{M}}(\Omega)$ and for all $K\in\mathcal{M}$, we denote by $u_{K}$ the constant value of $u$ in $K$. The space ${\rm H}_{\mathcal{M}}(\Omega)$ is equipped with the following classical Euclidean structure: for $(u,v)\in({\rm H}_{\mathcal{M}}(\Omega))^{2}$, we define the following inner product:

[TABLE]

Thanks to the discrete Poincaré inequality (5) given below, this scalar product defines a norm on ${\rm H}_{\mathcal{M}}(\Omega)$:

[TABLE]

The following discrete Poincaré inequality holds (see Lemma 9.1, p. 765, in [6]):

[TABLE]

We also define the following semi-inner product and semi-norm:

[TABLE]

These inner products can be seen as discrete analogues to the standard ${\rm H}^{1}$-inner product, with an implicitly assumed zero boundary condition in the case of the inner product $[\cdot,\cdot]_{\mathcal{M}}$.

For any function $u\in{\rm H}_{\mathcal{M}}(\Omega)$, we also define the following discrete ${\rm H}^{-1}$-norm:

[TABLE]

By inequality (5), the $\|\cdot\|_{-1,\mathcal{M}}$-norm is controlled by the ${\rm L}^{2}(\Omega)$-norm.

The discrete Laplace operators associated with homogeneous Dirichlet and Neumann boundary conditions, denoted by $\Delta_{\mathcal{M},{\rm D}}$ and $\Delta_{\mathcal{M},{\rm N}}$ respectively, are defined as follows:

[TABLE]

The links between these operators and the above defined inner products is clarified by the following identities:

[TABLE]

Finally, we suppose given a partition of the interval $(0,T)$, which we assume regular for the sake and simplicity, with $t^{0}=0,\dots,t^{n}=n\,\delta t,\dots t^{N}=T$.

A each time $t^{n}$, an approximation of the solution $(u^{n},\varphi^{n})\in{\rm H}_{\mathcal{M}}(\Omega)\times{\rm H}_{\mathcal{M}}(\Omega)$ is given by the following finite volume scheme:

[TABLE]

Note that in (7a), the term $u^{4}$ is discretized as $|u_{K}^{n+1}|\ (u_{K}^{n+1})^{3}$ to ensure positivity (see proof of Proposition 3.1 and Remark 3.1 below).

3. A priori ${\rm L}^{\infty}$ estimates, existence and uniqueness of the discrete solution

We prove in this section that the numerical scheme (7) admits a unique solution. The proof is based on a topological degree argument, which requires some a priori estimates on possible solutions to (7). Uniqueness is based on the positivity property of the scheme.

Proposition 3.1** (Existence and uniqueness of the approximate solution).**

1. (1)

The scheme (7) has a unique solution. 2. (2)

For $0\leq n\leq N$, the unknown $\varphi^{n}$ satisfies the following estimate:

[TABLE] 3. (3)

For $1\leq n\leq N$, the unknown $u^{n}$ satisfies the following estimate:

[TABLE]

Proof.

Step one: positivity of the unknowns.

We first observe that, from its definition (7c) and thanks to the fact that $u_{0}$ is non-negative, $u^{0}$ is a non-negative function. Let us suppose that this property still holds at time step $n$. We write the second equation of the scheme (7b) as:

[TABLE]

This set of relations is a linear system for $\varphi^{n}$, and the matrix of this system is clearly an M-matrix: indeed, from the definition (6) of $\Delta_{\mathcal{M},{\rm N}}$, it can be easily checked that its diagonal is strictly dominant and has only positive entries, and all its off-diagonal entries are non-positive. Since the right-hand side of this equation is non-negative, $\varphi^{n}$ is also non-negative. The first equation (7a) of the scheme can now be recast as:

[TABLE]

This set of relations can be viewed as a (nonlinear) matrix system for the unknown $u^{n+1}$; the matrix of this system depends on $u^{n+1}$ but is also an M-matrix for any values of $u^{n+1}$. Since we know that $\varphi^{n}\geq 0$, the right-hand side of this equation is by assumption non-negative and so is $u^{n+1}$.

Step two: upper bounds.

Let $\bar{\varphi}$ be a constant function of ${\rm H}_{\mathcal{M}}(\Omega)$. From the definition (6) of $\Delta_{\mathcal{M},{\rm N}}$, we see that $\Delta_{\mathcal{M},{\rm N}}(\bar{\varphi})=0$; the second relation of the scheme (7b) thus implies:

[TABLE]

Choosing for the constant value of $\bar{\varphi}$ the quantity $(\max_{K\in\mathcal{M}}u_{K}^{n+1})^{4}$ yields a non-positive right-hand side, and so, from the above mentioned property of the matrix of this linear system, $(\varphi^{n+1}-\bar{\varphi})_{K}\leq 0,\ \forall K\in\mathcal{M}$, which equivalently reads:

[TABLE]

Let us now turn to the estimate of the first unknown $u^{n+1}$. Let $K_{0}$ be a control volume where $u^{n+1}$ reaches its maximum value. From the definition (6) of $\Delta_{\mathcal{M},{\rm D}}$, it appears that:

[TABLE]

The first equation (7) of the scheme for the unknown $(u^{n+1}))_{K_{0}}$ reads:

[TABLE]

By the inequality (9), we see that supposing that $u_{K_{0}}^{n+1}>\max_{K\in\mathcal{M}}u_{K}^{n}$ yields that the first and third term of the preceding relation are positive, while the second one is non-negative, which is in contradiction with the fact that their sum is zero.

Step three: existence of a solution.

Let us suppose that we have obtained a solution to the scheme up to time step $n$. Let the function $F$ be defined as follows:

[TABLE]

It is clear that $u^{n+1}$ is a solution to (7a) if and only if:

[TABLE]

First, we observe that the function $F_{0}$, which maps ${\mathbb{R}}^{{\rm card}(\mathcal{M})}$ onto ${\mathbb{R}}^{{\rm card}(\mathcal{M})}$ and is defined by $F_{0}(U)=F(U,0)$, is affine and one-to-one. Second, a straightforward adaptation of steps one and two yields that the estimates on any possible solution $u^{n+1}$ to (11) hold uniformly for $\alpha\geq 0$. The existence of a solution to (11) then follows by a topological degree argument (see e.g. [3]).

Finally, the existence (and uniqueness) of the solution to the second equation (7b) of the scheme, which is a linear system, follows from the above mentioned properties of the associated matrix.

Step four: uniqueness of the solution.

Let us suppose that the solution is unique up to step $n$ and that there exist two solutions $u^{n+1}$ and $v^{n+1}$ to (7a). By the identity $a^{4}-b^{4}=(a-b)\,(a^{3}+a^{2}b+ab^{2}++b^{3})$, the difference $\delta u=u^{n+1}-v^{n+1}$ satisfies the following system of equations:

[TABLE]

Since we know from the precedent analysis that both $u_{K}^{n+1}$ and $v_{K}^{n+1}$ are non-negative, this set of relations can be seen as a matrix system for $\delta u$; the matrix of this system is again an M-matrix; we thus get $\delta u=0$, which proves the uniqueness of the solution. ∎

Remark 3.1** (On the non-negativity).**

We see from Relation (8) that the discretization of $u^{4}$ as a product of a positive quantity (here $|u^{n+1}|^{p}$) and $(u^{n+1})^{q}$ (where $p+q=4$) with $q$ odd (i.e. $q=3$ or $q=1$) is essential to prove the non-negativity of $u^{n+1}$; note that we have indeed observed in practice some (non-physical) negative values when this term is discretized as $(u^{n+1})^{4}$.

We also wish to emphasize that the non-negativity is obtained here thanks to the fact that the discretization is performed with a two-point flux finite volume scheme, which yields an $M$-matrix for the discrete Laplace operator. Hence we have only considered here the admissible meshes in the sense of [6, Definition 9.1, p. 762], since these meshes satisfy an orthogonality condition that yields the consistency of the two-point flux approximation. On general meshes, several consistent schemes have been designed, but these schemes do not ensure the positivity of the solution naturally; more on this subject may be found in e.g. [5].

4. Convergence to a solution of the continuous problem

Let ${\rm H}_{\mathcal{D}}$ be the space of piecewise constant functions over each $K\times I^{n}$, for $K\in\mathcal{M}$ and $I^{n}=[t^{n},\ t^{n+1}),\ 0\leq n\leq N-1$. To each sequence $(u^{n})_{n=0,N}$ of functions of ${\rm H}_{\mathcal{M}}(\Omega)$, we associate the function $u\in{\rm H}_{\mathcal{D}}$ defined by $u(x,t)=u^{n}(x)$ for $t^{n}\leq t<t^{n+1},\ 0\leq n\leq N-1$. In addition, for any $u\in{\rm H}_{\mathcal{D}}$, we define $\partial_{t,\mathcal{D}}(u)\in{\rm H}_{\mathcal{D}}$ by

[TABLE]

where the function $\partial_{t,\mathcal{D}}(u)^{n}\in{\rm H}_{\mathcal{M}}(\Omega)$ is defined by:

[TABLE]

For any function $u\in{\rm H}_{\mathcal{D}}$, we define the following norms and semi-norms:

[TABLE]

The norms $\|\cdot\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ and $|\cdot|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ can be seen as discrete equivalents of the ${\rm L}^{2}(0,T;{\rm H}^{1}(\Omega))$-norm, and $\|\cdot\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{-1})}$ may be considered as a discrete ${\rm L}^{2}(0,T;{\rm H}^{-1}(\Omega))$-norm.

Let us now derive some estimates on the solution to the scheme (7), which will be useful to get some compactness on sequences of approximate solutions.

Proposition 4.1** (Estimates in energy norms).**

Let $u$ and $\varphi$ be the functions of ${\rm H}_{\mathcal{D}}$ associated to the solution $(u^{n})_{0\leq n\leq N}\in{\rm H}_{\mathcal{M}}(\Omega)^{N+1}$ and $(\varphi^{n})_{0\leq n\leq N}\in{\rm H}_{\mathcal{M}}(\Omega)^{N+1}$ of the scheme (7). Then the following estimate holds:

[TABLE]

where the real number $c_{\rm e}$ only depends on $\Omega$, the initial data $u_{0}$ and (as a decreasing function) on the parameter $\theta_{\mathcal{M}}$ characterizing the regularity of the mesh, defined by (2).

Proof.

First, we recall that, by a standard reordering of the summations, we have, for any function $u\in{\rm H}_{\mathcal{M}}(\Omega)$:

[TABLE]

Multiplying Equation (7a) by $2\,\delta t\,|K|\,u_{K}^{n+1}$, using the equality $2\,a\,(a-b)=a^{2}+(a-b)^{2}-b^{2}$, summing over each control volume of the mesh and using the first of the preceding identities yields, for $0\leq n\leq N-1$:

[TABLE]

By Proposition 3.1, $u^{n+1}$ is non-negative and $\varphi^{n}_{K}\leq\bar{u}_{0}^{4},\ \forall K\in\mathcal{M}$, where $\bar{u}_{0}$ stands for $\max_{K\in\mathcal{M}}u_{K}^{0}$ and is thus bounded by the ${\rm L}^{\infty}$-norm of on the initial data $u_{0}$. Therefore, we get:

[TABLE]

By Cauchy-Schwarz’ inequality, Young’s inequality and the discrete Poincaré inequality (5), we thus obtain:

[TABLE]

Summing from $n=0$ to $n=N-1$, we get:

[TABLE]

Since, by assumption, $u_{0}\in{\rm H}^{1}_{0}(\Omega)$, the discrete ${\rm H}^{1}$ norm of $u^{0}$, $\|u^{0}\|_{1,\mathcal{M}}$ is bounded by $c\|u_{0}\|_{{\rm H}^{1}(\Omega)}$ where the real number $c$ only depends on $\Omega$ and, in a decreasing way, on the parameter $\theta_{\mathcal{M}}$ characterizing the regularity of the mesh (see e.g. Lemma 3.3 in [9]). Together with the preceding relation, this provides the control of the first term in (12).

We now turn to the estimate of $\|\partial_{t,\mathcal{D}}(u)\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{-1})}$. Let $v$ be a function of ${\rm H}_{\mathcal{M}}(\Omega)$; multiplying (7a) by $|K|\,v_{K}$ and summing over $K\in\mathcal{M}$, we get for $0\leq n\leq N-1$:

[TABLE]

Both $(u^{n+1})^{4}$ and $\varphi^{n}$ are non-negative functions which are bounded by $\bar{u}_{0}^{4}$, aand therefore the difference $(u^{n+1})^{4}-\varphi^{n}$ is itself bounded by $\bar{u}_{0}^{4}$; using the Cauchy-Schwarz inequality and the discrete Poincaré inequality (5), we thus get:

[TABLE]

and so:

[TABLE]

which, by the bound of $\|u\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$, yields the control of the second term in (12).

As far as $\varphi$ is concerned, the second equation (7b) of the scheme and the initialization (7c) yield for $0\leq n\leq N$:

[TABLE]

Thus, by Young’s inequality, we get:

[TABLE]

Multiplying by $\delta t$ and summing over the time steps, this provides the estimates of $\|\varphi\|_{L^{2}((0,T)\times\Omega)}$ and $|\varphi|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ we are seeking.

To obtain a control on $\partial_{t,\mathcal{D}}(\varphi)$, we need a sharper estimate on $\partial_{t,\mathcal{D}}(u)$. Our starting point is once again Equation (7)-$(i)$, which we multiply this time by $|K|\ \partial_{t,\mathcal{D}}(u)^{n}$ before summing over $K\in\mathcal{M}$, to get for $0\leq n\leq N-1$, once again by the identity $a^{2}-b^{2}\leq a^{2}+(a-b)^{2}-b^{2}=2\,a\,(a-b)$ and invoking the Cauchy-Schwarz inequality:

[TABLE]

so, by Young’s inequality:

[TABLE]

Multiplying by $\delta t$ and summing over the time steps yields:

[TABLE]

which provides an estimate for $\|\partial_{t,\mathcal{D}}(u)\|_{L^{2}((0,T)\times\Omega)}$. Taking now the difference of the second equation (7b) of the scheme at two consecutive time steps and using (7c) for the first step, we obtain for $0\leq n\leq N-1$:

[TABLE]

Multiplying by $|K|\ \partial_{t,\mathcal{D}}(\varphi)_{K}^{n}$ over each control volume of the mesh and summing yields:

[TABLE]

which, using Young’s inequality, multiplying by $\delta t$ and summing over the time steps yields the desired estimate for $\|\partial_{t,\mathcal{D}}(\varphi)\|_{L^{2}((0,T)\times\Omega)}$, thanks to (13). ∎

Let us now prove that sequences of solutions of the scheme (7) converge, as the space and time steps tend to 0 (up to a subsequence), to a weak solution of Problem (1), thereby proving the existence of this weak solution.

Theorem 4.2** (Convergence of the finite volume scheme).**

Let $(u^{(m)})_{m\in{\mathbb{N}}}$ and $(\varphi^{(m)})_{m\in{\mathbb{N}}}$ be a sequence of solutions to (7) with a sequence of discretizations such that the space and time step, $h^{(m)}$ and $\delta t^{(m)}$ respectively, tends to zero. We suppose that the parameters $\theta_{\mathcal{M}^{(m)}}$ characterizing the regularity of the meshes of this sequence are bounded away from zero, i.e. $\theta_{\mathcal{M}^{(m)}}\geq\theta>0,\ \forall\,m\in{\mathbb{N}}$. Then there exists a subsequence, still denoted by $(u^{(m)})_{m\in{\mathbb{N}}}$ and $(\varphi^{(m)})_{m\in{\mathbb{N}}}$ and two functions $\tilde{u}$ and $\tilde{\varphi}$ such that:

1. (1)

$u^{(m)}$* and $\varphi^{(m)}$ tends to $\tilde{u}$ and $\tilde{\varphi}$ respectively in ${\rm L}^{2}((0,T)\times\Omega)$,* 2. (2)

$\tilde{u}$* and $\tilde{\varphi}$ are solution to the continuous problem (1) in the following weak sense:*

[TABLE]

Proof.

The following estimates of the space translates can be found in [6], Lemma 9.3, p. 770 and Lemma 18.3, p. 851:

[TABLE]

where $\hat{v}$ stands for the extension by zero of $v$ to ${\mathbb{R}}^{d}$ and the real number $c(\Omega)$ only depends on the domain.

For $m$ given, let $\hat{u}^{(m)}$ and $\hat{\varphi}^{(m)}$ be the functions of ${\rm L}^{\infty}({\mathbb{R}}^{d}\times{\mathbb{R}})$ obtained by extending $u^{(m)}$ and $\varphi^{(m)}$ by [math] from $\Omega\times[0,T)$ to ${\mathbb{R}}^{d}\times{\mathbb{R}}$. The estimates of Proposition 4.1 of $\|u^{(m)}\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ and $|\varphi^{(m)}|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$, together with the ${\rm L}^{\infty}$-bound for $\varphi^{(m)}$, thus allow to bound independently of $m$ the space translates of $\hat{u}^{(m)}$ and $\hat{\varphi}^{(m)}$ in the ${\rm L}^{2}((0,T)\times\Omega)$-norm. In addition, Theorem A.1 applied with $\|\cdot\|_{\ast}$ equal to $\|\cdot\|_{1,\mathcal{M}}$ for $u^{(m)}$ and with $\|\cdot\|_{\ast}$ equal to $\|\cdot\|_{{\rm L}^{2}(\Omega)}$ for $\varphi^{(m)}$, together with the estimates of $\|u^{(m)}\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$, $\|\partial_{t,\mathcal{D}}(u)^{(m)}\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{-1})}$, $\|\varphi\|_{L^{2}((0,T)\times\Omega)}$ and $\|\partial_{t,\mathcal{D}}(\varphi)\|_{L^{2}((0,T)\times\Omega)}$ of Proposition 4.1 allows to bound the time translates of $\hat{u}^{(m)}$ and $\hat{\varphi}^{(m)}$, still independently of $m$ and in the ${\rm L}^{2}((0,T)\times\Omega)$-norm. Furthermore, by Proposition 3.1, the sequences $(\hat{u}^{(m)})_{m\in{\mathbb{N}}}$ and $(\hat{\varphi}^{(m)})_{m\in{\mathbb{N}}}$ are uniformly bounded in ${\rm L}^{\infty}((0,T)\times\Omega)$, and thus in ${\rm L}^{2}((0,T)\times\Omega)$. By Kolmogorov theorem (see e.g. [6], Theorem 14.1, p. 833), these sequences are relatively compact and strongly converge in ${\rm L}^{2}((0,T)\times\Omega)$ to, respectively, $\tilde{u}$ and $\tilde{\varphi}$.

The regularity of the limits, i.e. the fact that $\tilde{u}$ (resp. $\tilde{\varphi}$) lies in ${\rm L}^{2}(0,T;\,{\rm H}^{1}_{0}(\Omega))$ (resp. ${\rm L}^{2}(0,T;\,{\rm H}^{1}(\Omega))$ can be proven by invoking [6, Theorem 14.2 and Theorem 14.3] thanks to the uniform bounds on $\|u^{(m)}\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ and $|\varphi^{(m)}|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ given in Proposition 4.1, similarly to the technique of [8]; it can also be proven (again thanks to the estimates of Proposition 4.1) by invoking Theorem B.1 given in the appendix below, which gives a more straightforward way to obtain the regularity in a more general setting. The spaces which are considered in the application of Theorem B.1 are $B=L^{2}(\Omega)$, $X_{m}=H_{\mathcal{M}^{(m)}}$ and $X=H^{1}_{0}(\Omega)$.

To prove that $\tilde{u}$ and $\tilde{\varphi}$ are solution to the continuous problem, it remains to prove that (14) holds. This proof is rather standard (see e.g. the proof of Theorem 18.1 pp. 858–862 in [6] for a similar, although more complicated, problem) and we only give here the main arguments. Let $\psi$ be a function of ${\rm C}^{\infty}_{c}([0,T)\times\Omega)$. We define $\psi_{K}^{n}$ by $\psi_{K}^{n}=\psi(x_{K},t^{n})$. Multiplying the (7a) by $\delta t\,|K|\ \psi_{K}^{n+1}$ and summing up over the control volumes and the time steps, we get for any element of the sequence of discrete solutions:

[TABLE]

where, for enhanced readability, the superscript (m) has been omitted and:

[TABLE]

Reordering the summations and using the fact that $\psi(\cdot,T)=0$, we get for $T_{1}$:

[TABLE]

The first term of the right hand side reads:

[TABLE]

On one hand, $u^{0}$ converges to $u_{0}$ in ${\rm L}^{1}{(\Omega)}$ and $\psi(\cdot,0)\in{\rm C}^{\infty}_{c}(\Omega)$, so the second term of $T_{1,1}$ tends to zero with $h$; on the other hand, since $u_{K}^{0}\leq\bar{u}_{0}$, $\forall K\in\mathcal{M}$ and, from the regularity of $\psi$, $|R_{\psi}|\leq c_{\psi}\ (\delta t+h)$, the third term of $T_{1,1}$ also tends to zero with $\delta t$ and $h$. Let us now turn to the second term in the expression of $T_{1}$:

[TABLE]

with:

[TABLE]

and thus $|(R_{\psi})_{K}^{n}|\leq c_{\psi}\ (\delta t+h)$. Since $u_{K}^{n}\leq\bar{u}_{0}$, $\forall K\in\mathcal{M}$ and $0\leq n\leq N$, we get:

[TABLE]

Reordering the summations in $T_{2}$ and using the fact that $\psi(\cdot,T)=0$, we obtain:

[TABLE]

where the residual term $(R_{\psi})_{\sigma}^{n+1}$ is the difference of the mean value of $\nabla\psi\cdot n$ over $\sigma\times(t^{n+1},t^{n+2})$ and its finite volume approximation. The fact that the second term in the right hand side of this relation tends to zero is thus a classical consequence of the control of $\|u\|_{{\rm L}^{2}(0,T;\,{\rm H}_{\mathcal{M}}^{1})}$ and the consistency of the diffusive fluxes $(R_{\psi})_{\sigma}^{n+1}$ (see Theorem 9.1, pp. 772–776, in [6]) and yields, as $\tilde{u}$ is known to belong to ${\rm H}^{1}_{0}(\Omega)$:

[TABLE]

Finally, $T_{3}$ reads:

[TABLE]

with:

[TABLE]

The second term tends to zero by the ${\rm L}^{\infty}$-estimates for $u$ and $\varphi$ and the regularity of $\psi$. Since $u^{(m)}$ tends to $\tilde{u}$ in ${\rm L}^{2}((0,T)\times\Omega)$ and is bounded in ${\rm L}^{\infty}((0,T)\times\Omega)$, $u^{(m)}$ converges to $\tilde{u}$ in ${\rm L}^{2}(0,T;\,{\rm L}^{p}(\Omega))$, for any $p\in[1,+\infty)$; in addition, from the time translates estimates, $\varphi^{(m)}(\cdot,\cdot-\delta t)$ also converge to $\tilde{\varphi}$. We thus get:

[TABLE]

Gathering the results for $T_{1}$, $T_{2}$ and $T_{3}$, we obtain (14a). The relation (14b) is obtained using the same arguments; the convergence of the diffusion term in case of Neumann boundary conditions poses an additional difficulty which is solved in Theorem 10.3, pp. 810–815, in [6]. ∎

The uniqueness of the solution to the problem under consideration is beyond the scope of the present paper. Note however that such a result would imply, by a standard argument, the convergence of the whole sequence to the solution.

5. Conclusion

We propose in this paper a finite volume scheme for a problem capturing the essential difficulties of a simple radiative transfer model (the so-called P1 model), which enjoys the following properties: the discrete solution exists, is unique, and satisfies a discrete maximum principle; in addition, it converges (possibly up to the extraction of a subsequence) to a solution of the continuous problem, which yields, as a by-product, that such a solution indeed exists. For the proof of this latter result, we state and prove an abstract estimate allowing to bound the time translates of a finite volume discrete function, as a function of (possibly discrete) norms of the function itself and of its discrete time derivative; although this estimate is underlying in some already available analysis on finite volumes [6, Chapter IV], [7, 8] and discontinous Galerkin approximations [4], the present formulation is new and should be useful to tackle new problems. Variants of the presented numerical scheme are now successfully running for the modelling of radiative transfer in the ISIS free software [12] developed at IRSN and devoted to the simulation of fires in confined buildings [1], and in particular nuclear power plants.

Appendix A Estimation of time translates

The objective of this appendix is to state and prove an abstract result allowing to bound the time translates of a discrete solution.

We now introduce some notations. Let ${\rm H}_{\mathcal{M}}(\Omega)$ and ${\rm H}_{\mathcal{D}}$ be the discrete functional spaces introduced in section 2 and 4 respectively. We suppose given a norm $\|\cdot\|_{\ast}$ on ${\rm H}_{\mathcal{M}}(\Omega)$, over which we also define the dual norm $\|\cdot\|^{\ast}$ with respect to the ${\rm L}^{2}$-inner product:

[TABLE]

These two spatial norms may be associated to a corresponding norm on ${\rm H}_{\mathcal{D}}$ as follows:

[TABLE]

Theorem A.1**.**

Let $u$ be a function of ${\rm H}_{\mathcal{D}}$ and $\tau$ a real number. We denote by $\hat{u}$ the extension by zero of $u$ to ${\mathbb{R}}^{d}\times{\mathbb{R}}$. Then we have:

[TABLE]

Proof.

Let $u$ be a function of ${\rm H}_{\mathcal{D}}$ and $t\in{\mathbb{R}}$. Let $\tau$ be a real number that we suppose positive. The following identity holds:

[TABLE]

For $s\in{\mathbb{R}}$ we define $n(s)$ by: $n(s)=-1$ if $s<0$, $n(s)$ is the index such that $t^{n(s)}\leq s<t^{n(s)+1}$ for $0\leq s<t^{N}$, $n(s)=N+1$ for $t\geq t^{N}$. Let $n_{0}(t)$ and $n_{1}(t)$ be given by $n_{0}(t)=n(t)$, $n_{1}(t)=n(t+\tau)$. We adopt the convention $u^{-1}=u^{N}=0$. With this notation, we have for $u(\cdot,t+\tau)-u(\cdot,t)$ the following equivalent expression:

[TABLE]

and thus:

[TABLE]

Developping, we get:

[TABLE]

with:

[TABLE]

We first estimate the integral of $T_{1}(t)$ over ${\mathbb{R}}$. Since $\chi^{0}_{\tau}(t)$ is equal to $1$ in the interval $[-\tau,0)$ and [math] elsewhere, and since $u^{n_{0}(t)}=0$ for any negative $t$, we get:

[TABLE]

By the same arguments, we get the same bound for the integral of $T_{3}(t)$:

[TABLE]

From the definition of the $\|\cdot\|^{\ast}$ norm, we get:

[TABLE]

and thus, by Young’s inequality:

[TABLE]

Integrating over the time, we get:

[TABLE]

with

[TABLE]

Thans to Relation $(i)$ of the technical lemma A.2 given below, we have

[TABLE]

Since $u^{n_{0}(t)}=u^{m}$ for $t^{m}\leq t<t^{m+1}$, and thanks to Relation (16) of Lemma A.2 we get:

[TABLE]

Finally, $u^{n_{1}(t)}=u^{m}$ for $t^{m}-\tau\leq t<t^{m+1}-\tau$, and thus, by the same argument:

[TABLE]

This concludes the proof for positive $\tau$. The case of negative $\tau$ follows by symmetry. ∎

Lemma A.2**.**

Let $(t^{n})_{0\leq n\leq N}$ be such that $t^{0}=0$, $t^{n}=n\delta t$, $t^{N}=T$, $\tau$ be a positive real number and $\chi^{n}_{\tau}:\ {\mathbb{R}}\rightarrow{\mathbb{R}}$ be the function defined by $\chi^{n}_{\tau}(t)=1$ if $t<t^{n}\leq t+\tau$ and $\chi^{n}_{\tau}(t)=0$ otherwise. Then, for any family of real numbers $(\alpha_{n})_{n=1,N}$ and, respectively, for any real number $t$, we have the following identities:

[TABLE]

Proof.

The function $\chi^{n}_{\tau}(t)$ is equal to one for $t\in[t^{n}-\tau,\ t^{n})$, so we have:

[TABLE]

To obtain the inequality (16) we remark that $t\in[t^{n}-\tau,\ t^{n})$ is equivalent to $t-t^{n}\in[-\tau,\ 0)$ and so $\chi^{n}_{\tau}(t)=1$ is equivalent to $\chi^{0}_{\tau}(t-t^{n})=1$ (under the assumption $t^{0}=0$). We thus have:

[TABLE]

∎

Let us remark that the above theorem may be generalized to the Banach setting as in [8]. In fact, this has led to some adaptations of the Aubin-Simon compactness result to the discrete setting see [11, 2]. Note however that the compactness result is somewhat stronger here than in the general Banach setting, in the sense that the estimate on the time translates is bounded by a power of the translation. This is possible thanks to the fact that the norms $\|\cdot\|_{\ast}$ and $\|\cdot\|^{\ast}$ are dual.

Appendix B Regularity of the limit

Definition B.1** ($B$-limit included sequence of spaces).**

Let $B$ be a Banach space, $(X_{n})_{n\in{\mathbb{N}}}$ be a sequence of Banach spaces included in $B$ and $X$ be a Banach space included in $B$. We say that the sequence $(X_{n})_{n\in{\mathbb{N}}}$ is $B$-limit-included in $X$ if there exist $C\in{\mathbb{R}}$ such that if $u$ is the limit in $B$ of a subsequence of a sequence $(u_{n})_{n\in{\mathbb{N}}}$ verifying $u_{n}\in X_{n}$ and $\left\|u_{n}\right\|_{X_{n}}\leq 1$, then $u\in X$ and $\left\|u\right\|_{X}\leq C$.

Theorem B.1** (Regularity of the limit).**

Let $1\leq p<+\infty$ and $T>0$. Let $B$ be a Banach space, $(X_{n})_{n\in{\mathbb{N}}}$ be a sequence of Banach spaces included in $B$ and $B$-limit-included in $X$ (where $X$ is a Banach space included in $B$) in the sense of Definition B.1.

Let $T>0$ and, for ${n\in{\mathbb{N}}}$, let $u_{n}\in L^{p}((0,T),X_{n})$. Let us assume that the sequence $(\left\|u_{n}\right\|_{L^{p}((0,T),X_{n})})_{n\in{\mathbb{N}}}$ is bounded and that $u_{n}\to u$ in $L^{p}((0,T),B)$ as ${n\to+\infty}$. Then $u\in L^{p}((0,T),X)$.

Proof.

Since $u_{n}\to u$ in $L^{p}((0,T),B)$ as ${n\to+\infty}$, we can assume, up to subsequence, that $u_{n}\to u$ in $B$ a.e..

Then, since the the sequence $(X_{n})_{n\in{\mathbb{N}}}$ is $B$-limit-included in $X$, we obtain, with $C$ given by Definition B.1,

[TABLE]

Using now Fatou’s lemma, we have

[TABLE]

Then, since $(\left\|u_{n}\right\|_{L^{p}((0,T),X_{n})})_{n\in{\mathbb{N}}}$ is bounded, we conclude that $u\in L^{p}((0,T),X)$. ∎