Analysis of an alternative Navier-Stokes system: Weak entropy solutions and a convergent numerical scheme

Magnus Svärd Dept. of Mathematics, University of Bergen, P.O. Box 7803, 5020 Bergen, Norway,
[email protected]

Abstract

We consider an alternative Navier-Stokes model for compressible viscous ideal gases, originally proposed in [Svä18]. We derive a priori estimates that are sufficiently strong to support a weak entropy solution of the system. Guided by these estimates, we propose a finite volume scheme, derive the analogous estimates and demonstrate grid convergence towards a weak entropy solution. Furthermore, this existence proof is valid for “large” initial data and no a priori assumptions on the solution are needed.

1 Introduction

The importance of modelling viscous and heat conducting compressible fluids in physics and engineering can not be understated. Such fluids are typically modelled by the well-known Navier-Stokes-Fourier (NSF) equations. The complexity of these equations often necessitates numerical methods to generate approximate solutions. However, the design of numerical schemes is severely hampered by the lack of well-posedness results for the NSF system. This invariably leads to uncertainties pertaining to the validity of numerical results obtained with state-of-the-art codes. (See e.g. [GSH21].)

The NSF equations model the dynamics of (Lagrangian) mass elements by Newton’s 2nd Law, which by definition precludes mass diffusion. (See [Svä18, Sto45].) The effect of the random motions at the molecular level enters the NSF system via the stress tensor and a model of viscosity that exerts a force on each mass element. Furthermore, heat diffusion is modelled by Fourier’s law in the energy equation. A consequence of the Lagrangian perspective is that the continuity equation becomes hyperbolic, which poses a significant mathematical challenge when trying to develop a well-posedness theory.

In recent years, many authors have, for physical reasons, proposed alternative models. See e.g. [RDO⁺19] where a class of mass-diffusive continuum models is derived and demonstrated to be more accurate than the NSF system in many cases. Another well-known mass-diffusive modification of the NSF system was given in [Bre05a, Bre05b]. For a constant mass-diffusive coefficient, this system was proven to have weak solutions in [FV10].

Herein, we consider the alternative Navier-Stokes model proposed in [Svä18]. Let ${\bf x}=(x,y,z)$ be the Cartesian coordinates and $\Omega$ a bounded spatial domain. Then the model is given by

$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho{\bf v})$	$\displaystyle=\nabla\cdot(\nu\nabla\rho),$
$\displaystyle\partial_{t}(\rho{\bf v})+\nabla\cdot(\rho{\bf v}\otimes{\bf v})+\nabla p$	$\displaystyle=\nabla\cdot(\nu\nabla\rho{\bf v})),\qquad t\in[0,\mathcal{T}],\,{\bf x}\in\Omega\subset\mathbb{R}^{3}$	(1)
$\displaystyle\partial_{t}(E)+\nabla\cdot(E{\bf v}+p{\bf v})$	$\displaystyle=\nabla\cdot(\nu\nabla E),$
$\displaystyle p$	$\displaystyle=\rho RT,\quad\textrm{(ideal gas law).}$

$\rho,\rho{\bf v},E$ are the conserved variables density, momentum and total energy. ${\bf v}=(v_{1},v_{2},v_{3})^{T}=(u,v,w)$ is the velocity vector. $p,T$ are pressure and temperature. Furthermore, the total energy is given by $E=\frac{p}{\gamma-1}+\frac{\rho|{\bf v}|^{2}}{2}$ where $\gamma=c_{p}/c_{v}$ and $c_{p,v}$ are the heat capacities at constant pressure or volume. The model is (at least) valid for ideal gases with $1<\gamma\leq 5/3$ . The conservative variables are collected in the vector ${\bf u}=(\rho,{\bf m}^{T},E)^{T}$ where ${\bf m}=(m_{1},m_{2},m_{3})^{T}=\rho{\bf v}$ is the momentum vector. $R$ is the gas constant and $\nu=\nu(\rho,T)$ is the diffusion coefficient.

An ideal gas is one where the molecules bounce elastically against each other and have no rotational energy. Furthermore, the mean free path is large relatively to their sizes. In such a gas, the range of the diffusion is proportional to the mean free path which scales as the inverse of density. Hence, it was proposed in [Svä18] that,

\displaystyle\nu=\frac{\mu_{0}}{\rho},

(2)

where $\mu_{0}$ is the dynamic viscosity. This value results in viscous terms that resemble those of the Navier-Stokes equations, and Fourier’s law appears as a diffusive term in the energy equation of (1). In [DS21] and [SDP21], it was shown, in a suite of problems, that the new system (1)-(2) produces solutions that are next to indistinguishable from those of the NSF system. (The differences are much less than what can be measured in experiments.)

Remark 1.

With $\mu_{0}$ being the dynamic viscosity, the value of $\nu$ is identical to the kinematic viscosity. However, we emphasise that physically $\nu$ is a diffusion coefficient, and not a viscosity coefficient.

Unlike the NSF system, (1) is derived in an Eulerian¹¹1The system (1), was termed “Eulerian model” in [Svä18]. Here, we call it an “alternative Navier-Stokes model” since the original name has caused some confusion that it is a (numerical) regularisation of the Euler equations rather than an alternative to the standard Navier-Stokes model. frame by directly appealing to the principles of conservation of mass, momentum, and energy. In the Eulerian perspective, the random motions of molecules are naturally modelled as diffusion in all three conservative variables. (Indeed, (1) differs from the NSF system only in the right-hand side.)

Our goal is to demonstrate that the system (1) has weak entropy solutions by proving that solutions to a finite volume scheme converge as the grid is refined. In this effort, we have benefited greatly from [FV10] and have used many similar techniques and arguments.

1.1 Adding more physics to the model

Although (1) is accurate in the regime it is intended to model, considerations for reduced models (isothermal and isentropic) reveal that (1) appears to lack mechanisms that prevent concentration of mass. Hence, we propose some minor modifications to (1) below. These modifications do not extend validity range of the model and can thus alternatively be viewed as technical assumptions.

Nevertheless, we continue with a short physical discussion: The choice $\nu=\mu_{0}/\rho$ models diffusion whereby molecules of the ideal gas travel relatively freely in space. Their random crossings of (Eulerian) control-volume boundaries, is modelled as diffusion on the global scale and the more rarefied the gas is, the more significant is the diffusion relatively to the convection. This is consistent with the diffusion coefficient being proportional to the mean free path, i.e., $\rho^{-1}$ . However, this mechanism is not prevalent when the gas becomes denser since molecules will bounce into other molecules before they have travelled any significant distance. The latter process is conductive rather than diffusive and not accounted for with $\nu=\mu_{0}/\rho$ . Hence, the physics suggests that the model should be augmented with collisional/frictional diffusion.

There are several options to account for collisional diffusion: The dynamic viscosity coefficient appearing in the Navier-Stokes equations is often assumed to depend on temperature as $T^{r}$ , for some $0<r<1$ . Although $\nu$ has a very different physical interpretation it might be tempting to try $\nu\sim T^{r}/\rho$ . However, some preliminary (unpublished) numerical investigations suggest that a temperature dependent diffusion is less accurate for (1) than $\nu=\mu_{0}/\rho$ where $\mu_{0}=constant$ . Hence, we will not consider models with $\nu\sim T^{r}/\rho$ .

In a dense gas, the molecules will constantly experience the repelling forces of neighbouring molecules, which will randomly deviate their paths, which in turn is a diffusive process at the global scale. Moreover, the more densely packed the gas is, the more prominent will the internal friction be in relation to the convective process. Hence, it seems plausible that there is some $\rho^{s}$ , $s>0$ dependence of $\nu$ . We take $s=1$ . Since $\nu=\mu_{0}/\rho$ , gives accurate results under normal conditions, we assume that the diffusion coefficient is given by,

\displaystyle\nu=\frac{\mu_{0}}{\rho}+\mu_{1}\rho,\quad\mu_{0}>>\mu_{1}>0.

(3)

Finally, to control pressure in the full system, a sufficiently strong estimate of temperature is also needed. Such estimate is unattainable for the base model and we add a Fourier-type heat flux that models radiation. That is, $\kappa\nabla T=4\kappa_{r}T^{3}\nabla T$ , where $\kappa_{r}$ is a material dependent constant. Under normal conditions (modest temperatures), $4\kappa_{r}T^{3}$ should be orders of magnitude less than $\mu_{0}$ such that the radiation has a negligible effect.

Remark 2.

The extra radation term does not increase the physical validity range. For that, more aspects of radiation physics would have to be included in the model. In fact, if the radiation diffusion is significant, the solution is already outside the validity range. In this sense, one can view the extra temperature diffusion as a technical modification.

In the remainder of this text, we consider the following alternative Navier-Stokes system:


$\displaystyle\partial_{t}\rho+\nabla\cdot(\rho{\bf v})$	$\displaystyle=\nabla\cdot(\nu\nabla\rho),$	(4a)
$\displaystyle\partial_{t}(\rho{\bf v})+\nabla\cdot(\rho{\bf v}\otimes{\bf v})+\nabla p$	$\displaystyle=\nabla\cdot(\nu\nabla\rho{\bf v})),\qquad t\in[0,\mathcal{T}]$	(4b)
$\displaystyle\partial_{t}(E)+\nabla\cdot(E{\bf v}+p{\bf v})$	$\displaystyle=\nabla\cdot(\nu\nabla E)+\nabla\cdot(\kappa_{r}\nabla T^{4}),$	(4c)
$\displaystyle p$	$\displaystyle=\rho RT,\quad\textrm{ideal gas law,}$	(4d)

with $\nu$ given by (3). Finally, we take $\Omega\subset\mathbb{R}^{3}$ to be bounded, and its boundary, $\partial\Omega$ , models an adiabatic wall, which is given by the no-slip boundary condition

\displaystyle{\bf v}=0,

(5)

along with

\displaystyle\partial_{n}T=\partial_{n}\rho=0,

(6)

where $\partial_{n}={\bf n}\cdot\nabla$ and ${\bf n}$ is the wall normal. (See [Svä18] for a discussion on boundary conditions.)

1.2 Outline

•

In Section 2, auxiliary balance laws are derived.
•

Using the equations and the auxiliary balance laws, we derive a priori estimates in Section 3. As a consequence of the estimates, we are able to demonstrate positivity of $\rho$ and $T$ in Section 3.3.
•

In Section 3.8, we formally demonstrate that the a priori estimates supports a weak entropy solution.
•

In Section 3.9, we define weak entropy solutions.
•

Next, we propose a numerical scheme in Section 4.1 and give the main result of this paper: Existence of a weak entropy solution to (4). The remainder of the paper proves this result.
•

Some auxiliary (numerical) balance laws are derived in Section 5.
•

We derive a priori estimates for the numerical scheme in Section 6.
•

In Section 7, we use the global estimates to infer existence of weak solutions to the full system by proving weak convergence of the numerical scheme.
•

Section 8, contains some concluding remarks.

1.3 Notation

The following notation will be used throughout the article.

•

The space $L^{p}(0,\mathcal{T};L^{q}(\Omega))$ is denoted $L^{p}_{q}$ for short, when there is no risk of confusion.
•

Similarly, we write $S_{1}(S_{2})=S_{1}(0,\mathcal{T};S_{2}(\Omega))$ where $S_{1,2}$ are two generic function spaces.
•

$W^{r,p}(\Omega)$ is the Sobolev space where all derivatives of order $0,1,...,r$ are bounded in $L^{p}(\Omega)$ . The space $W^{1,2}$ is denoted $H^{1}$ . Furthermore, the norm of $W^{r,p}(\Omega)$ is denoted $\|\cdot\|_{r,p}$ .
•

Throughout, we let $\mathcal{C}$ denote a generic a priori bounded constant obtained from the initial data. Likewise, $\epsilon$ and $\delta$ are used as small positive constants. Note that all three constants may change their actual values throughout a calculation.
•

We denote the $i$ th component of a vector ${\bf a}$ as ${\bf a}_{i}$ .
•

In the remainder of the paper, we will use a number of standard results which have been collected in Appendix I for the reader’s convenience.

2 Additional balance laws

2.1 Entropy balance

The first relation we derive is an entropy equation, which is an expression of the Second Law of Thermodynamics. The specific entropy is $S=\log(p/\rho^{\gamma})$ and here we take the standard approach and derive an equation for the entropy function $U=-\rho S$ . Associated with the entropy function $U$ are the entropy fluxes ${\bf F}_{i}=-m_{i}S$ , $i=1,2,3$ , and the entropy variables,

U_{\bf u}={\bf w}^{T}=(-(S-\gamma)-\frac{|{\bf v}|^{2}}{2c_{v}T},\frac{{\bf v}}{c_{v}T},-\frac{1}{c_{v}T}).

Contracting (4) with the entropy variables gives,

\displaystyle U_{t}+\nabla\cdot{\bf F}=\sum_{i=1}^{5}{\bf w}_{i}\nabla\cdot(\nu\nabla{\bf u}_{i})+{\bf w}_{5}\nabla\cdot(\kappa_{r}\nabla T^{4}),

which can be recast as,

\displaystyle U_{t}+\nabla\cdot{\bf F}=\sum_{i=1}^{5}\left(\nabla\cdot({\bf w}_{i}\nu\nabla{\bf u}_{i})-(\nabla{\bf w}_{i}^{T})\nu(\nabla{\bf u}_{i})\right)+\nabla\cdot({\bf w}_{5}\kappa_{r}\nabla T^{4})-\nabla{\bf w}_{5}\cdot(\kappa_{r}\nabla T^{4}),

and finally,

$\displaystyle(-\rho S)_{t}+\nabla\cdot(-{{\bf m}}S)$	$\displaystyle=$
$\displaystyle\nabla\cdot\left(-(S-\gamma)\nu\nabla\rho+\frac{\nu}{c_{v}T}\nabla\left(\frac{\rho\|{\bf v}\|^{2}}{2}\right)-\frac{1}{c_{v}T}\nu\nabla e-\frac{1}{c_{v}T}\kappa_{r}\nabla T^{4}\right)$		(7)
$\displaystyle-\nu(\gamma-1)\rho\|\nabla\ln\rho\|^{2}-\nu\rho\|\nabla\ln T\|^{2}-\frac{\nu\rho}{c_{v}T}\sum_{i=1}^{3}\|\nabla v_{i}\|^{2}-4\kappa_{r}T^{3}\frac{\|\nabla T\|^{2}}{T^{2}},$

which is the entropy balance.

2.2 Kinetic Energy balance

To obtain the kinetic energy equation for the model (4), we multiply the momentum equation by ${\bf v}$ . For the velocity terms on the left-hand side, we obtain

	$\displaystyle T1={\bf v}\cdot(\partial_{t}(\rho{\bf v})+\nabla\cdot(\rho{\bf v}\otimes{\bf v}))$	$\displaystyle=$
	$\displaystyle(\partial_{t}\rho+\nabla\cdot(\rho{\bf v}))\|{\bf v}\|^{2}+\rho(\partial_{t}{\bf v}+{\bf v}\cdot\nabla{\bf v})\cdot{\bf v}$	$\displaystyle=$
	$\displaystyle(\partial_{t}\rho+\nabla\cdot(\rho{\bf v}))\|{\bf v}\|^{2}+\rho\partial_{t}\frac{\|{\bf v}\|^{2}}{2}+\rho{\bf v}\cdot\nabla\frac{\|{\bf v}\|^{2}}{2}$	$\displaystyle=$
	$\displaystyle\partial_{t}\left(\frac{\rho\|{\bf v}\|^{2}}{2}\right)+\nabla\cdot\left(\frac{\rho\|{\bf v}\|^{2}}{2}{\bf v}\right)+(\partial_{t}\rho+\nabla\cdot(\rho{\bf v}))\frac{\|{\bf v}\|^{2}}{2}.$

Using the continuity equation (4a) in the last line, we obtain

\displaystyle T1=\partial_{t}\left(\frac{\rho|{\bf v}|^{2}}{2}\right)+\nabla\cdot\left(\frac{\rho|{\bf v}|^{2}}{2}{\bf v}+\nu\frac{|{\bf v}|^{2}}{2}\nabla\rho\right)-\nu\nabla\rho\cdot\nabla\frac{|{\bf v}|^{2}}{2}.

(8)

Next, we turn to the remaining terms of the momentum equation:

	$\displaystyle T2={\bf v}\cdot(-\nabla p+\nabla\cdot(\nu\nabla\rho{\bf v})))$	$\displaystyle=$
	$\displaystyle-\nabla\cdot(p{\bf v})+p\nabla\cdot{\bf v}+{\bf v}\cdot(\nabla\cdot(\nu\nabla\rho{\bf v}))$	$\displaystyle=$
	$\displaystyle-\nabla\cdot(p{\bf v})+p\nabla\cdot({\bf v})+\nabla\cdot\left(\nu\nabla\rho\|{\bf v}\|^{2}+\nu\rho\nabla\frac{\|{\bf v}\|^{2}}{2}\right)-\nu\rho\|\nabla{\bf v}\|^{2}-\nu\nabla\rho\nabla\frac{\|{\bf v}\|^{2}}{2}$	.

Combining $T1$ and $T2$ , yields the kinetic energy balance,

	$\displaystyle\partial_{t}\left(\frac{\rho\|{\bf v}\|^{2}}{2}\right)+\nabla\cdot\left(\frac{\rho\|{\bf v}\|^{2}}{2}{\bf v}+\nu\frac{\|{\bf v}\|^{2}}{2}\nabla\rho\right)$	$\displaystyle=$
	$\displaystyle-\nabla\cdot(p{\bf v})+p\nabla\cdot({\bf v})+\nabla\cdot\left(\nu(\nabla\rho)\|{\bf v}\|^{2}+\nu\rho\nabla\frac{\|{\bf v}\|^{2}}{2}\right)-\nu\rho\|\nabla{\bf v}\|^{2},$

that simplifies to

	$\displaystyle\partial_{t}\left(\frac{\rho\|{\bf v}\|^{2}}{2}\right)+\nabla\cdot\left(\frac{\rho\|{\bf v}\|^{2}}{2}{\bf v}\right)$	$\displaystyle=$
	$\displaystyle-\nabla\cdot(p{\bf v})+p\nabla\cdot({\bf v})+\nabla\cdot\left(\nu\nabla\left(\rho\frac{\|{\bf v}\|^{2}}{2}\right)\right)-\nu\rho\|\nabla{\bf v}\|^{2}.$			(9)

2.3 Specific volume balance

Multiply the continuity equation (4a) by $-\rho^{-2}$ , to obtain

\displaystyle\frac{1}{2}\partial_{t}\rho^{-1}+\nabla\cdot(\rho^{-1}{\bf v})+\nabla(\rho^{-2})\cdot\rho{\bf v}=-\nabla(\rho^{-2}\nabla(\nu\rho))-\frac{2}{\rho^{3}}\nabla\rho\cdot(\nu\nabla\rho),

or,

\displaystyle\frac{1}{2}\partial_{t}\rho^{-1}+\nabla\cdot(\rho^{-1}{\bf v})+2\nabla(\rho^{-1})\cdot{\bf v}=-\nabla(\rho^{-2}\nabla(\nu\rho))-2\nu\rho|\nabla\rho^{-1}|^{2}.

(10)

2.4 Renormalised internal energy balance

Subtracting the kinetic energy balance (9), from the total energy equation (4c), gives the internal energy balance

\displaystyle\partial_{t}(c_{v}\rho T)+\nabla\cdot(c_{v}\rho T{\bf v})

\displaystyle=-p\nabla\cdot({\bf v})+\nabla\cdot(\nu\nabla(c_{v}\rho T))+\nu\rho|\nabla{\bf v}|^{2}+\nabla\cdot(\kappa_{r}\nabla T^{4}).

(11)

We follow [FV10] and introduce $H(T)=(1+T)^{1-\omega}$ , $\omega>0.$ Then we renormalise the equation by multiplying (11) by $H^{\prime}(T)$ . We carry out the procedure step by step and begin with the left-hand side. First,

	$\displaystyle H^{\prime}(T)(c_{v}\rho_{t}T+c_{v}\rho T_{t})$	$\displaystyle=$
	$\displaystyle c_{v}(H^{\prime}(T)\rho_{t}T+\rho H_{t})+\rho_{t}c_{v}H-\rho_{t}c_{v}H$	$\displaystyle=$
	$\displaystyle c_{v}(\rho H)_{t}+c_{v}(H^{\prime}(T)T-H)\rho_{t}.$

In the same way,

\displaystyle H^{\prime}(T)(c_{v}\nabla\cdot(\rho{\bf v}T))

\displaystyle=c_{v}\nabla\cdot(\rho{\bf v}H)+c_{v}(H^{\prime}(T)T-H)\nabla\cdot(\rho{\bf v}).

Using the continuity equation, the left-hand side of the renormalised equation becomes,

\displaystyle c_{v}(\rho H)_{t}+c_{v}\nabla\cdot(\rho{\bf v}T))+c_{v}(H^{\prime}(T)T-T)\nabla\cdot(\nu\nabla\rho).

Turning to the right-hand side, we begin with

	$\displaystyle H^{\prime}\nabla\cdot(\nu\nabla(c_{v}\rho T))=c_{v}H^{\prime}\nabla\cdot(\nu T\nabla\rho+\nu\rho\nabla T))$	$\displaystyle=$
	$\displaystyle c_{v}H^{\prime}T\nabla\cdot(\nu\nabla\rho)+c_{v}H^{\prime}\nu\nabla\rho\cdot\nabla T+c_{v}\nabla\cdot(H^{\prime}\nu\rho\nabla T))-c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}$	$\displaystyle=$
	$\displaystyle c_{v}H^{\prime}T\nabla\cdot(\nu\nabla\rho)+c_{v}\nu\nabla\rho\cdot\nabla H+c_{v}\nabla\cdot(\nu\rho\nabla H))-c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}$	$\displaystyle=$
	$\displaystyle c_{v}H^{\prime}T\nabla\cdot(\nu\nabla\rho)+c_{v}\nabla\cdot(\nu H\nabla\rho)-c_{v}H\nabla\cdot(\nu\nabla\rho)+c_{v}\nabla\cdot(\nu\rho\nabla H))-c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}$	$\displaystyle=$
	$\displaystyle c_{v}(H^{\prime}T-H)\nabla\cdot(\nu\nabla\rho)+c_{v}\nabla\cdot(\nu H\nabla\rho)+c_{v}\nabla\cdot(\nu\rho\nabla H))-c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}$	.

Collecting the terms, we obtain the renormalised internal energy balance,

$\displaystyle c_{v}(\rho H)_{t}+c_{v}\nabla\cdot(\rho{\bf v}H))$	$\displaystyle=$
$\displaystyle c_{v}\nabla\cdot(\nu H\nabla\rho)+c_{v}\nabla\cdot(\nu\rho\nabla H))-c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}$		(12)
$\displaystyle+H^{\prime}p\nabla\cdot{\bf v}+H^{\prime}\nu\rho\|\nabla{\bf v}\|^{2}+H^{\prime}\nabla(\kappa_{r}\nabla T^{3}).$

3 A priori estimates for the full system

We begin with the standard conservation properties. That is, we assume that $p,\rho>0$ and (separately) integrate (4a) and (4c). Using the boundary conditions (5) and (6) we obtain the bounds,

	$\displaystyle\rho\in L^{\infty}(0,\mathcal{T};L^{1}(\Omega)),$	$\displaystyle\quad\quad E\in L^{\infty}(0,\mathcal{T};L^{1}(\Omega)),$		(13)
	$\displaystyle\sqrt{\rho}\|{\bf v}\|\in L^{\infty}(0,\mathcal{T};L^{2}(\Omega)),$	$\displaystyle\quad\quad p\in L^{\infty}(0,\mathcal{T};L^{1}(\Omega)).$

3.1 Entropy estimate

Integrating (7) in space and time, using the boundary conditions (5) and (6), leads to

	$\displaystyle\int_{\Omega}(-\rho S)\|_{t=\mathcal{T}}-(-\rho S)\|_{t=0}\,d{\bf x}$	$\displaystyle=$		(14)
	$\displaystyle-\int_{0}^{\mathcal{T}}\int_{\Omega}\nu(\gamma-1)\rho\|\nabla\ln\rho\|^{2}+\nu\rho\|\nabla\ln T\|^{2}+\frac{\nu\rho}{c_{v}T}\sum_{i=1}^{3}\|\nabla v_{i}\|^{2}+4\kappa_{3}T^{3}\frac{\|\nabla T\|^{2}}{T^{2}}\,d{\bf x}\,dt.$

The second term is bounded by initial data. Provided that we can control the first, we obtain the bounds

\displaystyle\int_{0}^{\mathcal{T}}\int_{\Omega}\nu(\gamma-1)\rho|\nabla\ln\rho|^{2}+\nu\rho|\nabla\ln T|^{2}+\frac{\nu\rho}{c_{v}T}\sum_{i=1}^{3}|\nabla v_{i}|^{2}+4\kappa_{r}T|\nabla T|^{2}\,d{\bf x}\,dt\leq\mathcal{C}.

(15)

Hence, we proceed to analyse the first integral of (14). (This was shown in [FV10] and we repeat their arguments here.) We introduce $[z]^{+}=\max(z,0)$ and $[z]^{-}=\min(z,0)$ and carry out the following manipulations:

	$\displaystyle\int_{\Omega}(-\rho S)\,d{\bf x}=\int_{\Omega}-\rho\log(RT)+(\gamma-1)\rho\log(\rho)\,d{\bf x}=$
	$\displaystyle\int_{\Omega}-\rho([\log(RT)]^{+}+[\log(RT)]^{-})+(\gamma-1)\rho([\log(\rho)]^{+}+[\log(\rho)^{-}]\,d{\bf x}.$		(16)

We need a bound on the negative terms. $\int_{\Omega}-\rho[\log(RT)]^{+}\,d{\bf x}$ is bounded by $p\in L^{\infty}_{1}$ given by (13). Since $|\rho[\log(\rho)]^{-}|\leq 1$ and the domain is bounded, the last term is also bounded. Hence, we have deduced that

	$\displaystyle\rho\log(\rho)$	$\displaystyle\in L^{\infty}(0,\mathcal{T},L^{1}(\Omega)),$		(17)
	$\displaystyle\rho\log(T)$	$\displaystyle\in L^{\infty}(0,\mathcal{T},L^{1}(\Omega)),$

and, with $\nu\sim\rho^{-1}+\rho$ , we obtain from the diffusive terms (15),

$\displaystyle\nabla(\log(T))$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2})),$
$\displaystyle\rho\nabla(\log(T))$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2})),$
$\displaystyle\nabla(\log(\rho))$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2})),$	(18)
$\displaystyle\nabla\rho$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2})),$
$\displaystyle\frac{\rho}{\sqrt{T}}\|\nabla{\bf v}\|$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2})),$
$\displaystyle\frac{1}{\sqrt{T}}\|\nabla{\bf v}\|$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2})).$

From the radiation dissipation (the last term of (15)), we obtain

\displaystyle\nabla T^{3/2}

\displaystyle\in L^{2}(0,\mathcal{T};L^{2})).

(19)

From the $\nabla\rho$ estimate in (18), the $\rho\in L^{\infty}_{1}$ bound given in (13) and Poincare’s inequality, we also have,

\displaystyle\rho\in L^{2}_{2}.

From the current estimates, it is possible to obtain a better estimate on temperature. We begin by showing that $\rho>\delta>0$ on a subset of finite measure. (Again, we repeat the arguments laid out in [FV10].)

That is, we wish to show that

\displaystyle|\mathcal{B}|=|\{x\in\Omega:\,|\rho({\bf x},t)>\delta\}|>m>0.

(20)

First, let $M_{0}$ be the total mass, which is constant thanks to the boundary conditions. Then consider the bound (17). Since $|\rho\log(\rho)|>\rho$ for $0<\rho<<1$ and $\rho>>1$ , there must exist an $\alpha>0$ such that,

\displaystyle\int_{\{\rho(\cdot,t)>\alpha\}}\rho(\cdot,t)\,d{\bf x}\leq\frac{M_{0}}{3}.

We write,

	$\displaystyle M_{0}=\int_{\Omega}\rho(\cdot,t)\,d{\bf x}$	$\displaystyle=$
	$\displaystyle\int_{\{\rho(\cdot,t)\leq\delta\}}\rho(\cdot,t)\,d{\bf x}+\int_{\{\delta<\rho(\cdot,t)\leq\alpha\}}\rho(\cdot,t)\,d{\bf x}+\int_{\{\rho(\cdot,t)>\alpha\}}\rho(\cdot,t)\,d{\bf x}$	$\displaystyle\leq$
	$\displaystyle\delta\|\Omega\|+\alpha m+\frac{M_{0}}{3}.$

We choose $\delta=M_{0}/(3|\Omega|)$ and deduce that $m\geq M_{0}/(3\alpha)>0$ .

Next, we use $p\in L^{\infty}_{1}$ from (13). Thanks to $\mathcal{B}$ having non-zero measure, we conclude that $T\in L^{\infty}(0,\mathcal{T};L^{1}(\mathcal{B}))$ . Hence, by a generalised Poincare inequality (Theorem 2, Eqn. (165) in Appendix I) and (19), we have

	$\displaystyle T^{3/2}$	$\displaystyle\in L^{2}(H^{1})\quad\textrm{and by Sobolev embedding},$
	$\displaystyle T$	$\displaystyle\in L^{3}(L^{9}).$		(21)

3.2 Kinetic energy estimate

Integrating the balance (9), using the boundary conditions (5) and (6), yields,

\displaystyle\int_{\Omega}\partial_{t}\left(\frac{\rho|{\bf v}|^{2}}{2}\right)\,d{\bf x}

\displaystyle=\int_{\Omega}\left(p\nabla\cdot({\bf v})-\nu\rho|\nabla{\bf v}|^{2}\right)\,d{\bf x}.

(22)

We need to control $p\nabla\cdot({\bf v})$ , which is accomplished by

\displaystyle\int_{0}^{T}\int_{\Omega}p\nabla\cdot({\bf v})\,d{\bf x}dt\leq\int_{0}^{T}\int_{\Omega}\eta^{-1}T^{2}+\eta\rho^{2}|\nabla\cdot({\bf v})|^{2}\,d{\bf x}dt.

(23)

The first term on the right-hand side is controlled by (21). Since $\nu\sim\mu_{1}\rho$ , we can control the last term by choosing $\eta<\mu_{1}/2$ and from (22) we thus obtain the estimates,

	$\displaystyle\int_{0}^{T}\\|\rho\nabla{\bf v}\\|^{2}_{2}\,dt\leq\mathcal{C},$		(24)
	$\displaystyle\int_{0}^{T}\\|\nabla{\bf v}\\|^{2}_{2}\,dt\leq\mathcal{C}.$

Since $\sqrt{\rho}{\bf v}\in L^{2}_{2}$ by (13), and $\rho>\delta$ on $\mathcal{B}$ , we have that ${\bf v}\in L^{2}(0,\mathcal{T};L^{2}(\mathcal{B}))$ . Using the generalised Poincare inequality (165), we obtain

	$\displaystyle{\bf v}$	$\displaystyle\in L^{2}(H^{1}),\quad\textrm{and by Sobolev embedding}$		(25)
	$\displaystyle{\bf v}$	$\displaystyle\in L^{2}_{6}.$

3.3 Estimates of the specific volume

Integrating (10) in space and time, using the boundary conditions (5) and (6), gives,

\displaystyle\int_{0,\Omega}^{\mathcal{T}}\frac{1}{2}\partial_{t}\rho^{-1}\,d{\bf x}\,dt+2\int_{0,\Omega}^{\mathcal{T}}\nabla(\rho^{-1})\cdot{\bf v}\,d{\bf x}\,dt=-\int_{0,\Omega}^{\mathcal{T}}2\nu\rho|\nabla\rho^{-1}|^{2}\,d{\bf x}\,dt,

(26)

and using Young’s inequality, gives,

\displaystyle\int_{0,\Omega}^{\mathcal{T}}\frac{1}{2}\partial_{t}\rho^{-1}d{\bf x}\,dt\leq 2\int_{0,\Omega}^{\mathcal{T}}\eta|\nabla(\rho^{-1})|^{2}+\eta^{-1}|{\bf v}|^{2}d{\bf x}\,dt-2\int_{0,\Omega}^{\mathcal{T}}\nu\rho|\nabla\rho^{-1}|^{2}d{\bf x}\,dt.

(27)

Using (25), choosing $\eta\leq\mu_{0}/2$ and noting that $\rho>0$ , we obtain the bounds,

	$\displaystyle\sup_{t\in[0,\mathcal{T}]}\\|\rho^{-1}\\|_{1}\leq\mathcal{C},$		(28)
	$\displaystyle\int_{0}^{T}\\|\nabla(\rho^{-1})\\|^{2}_{2}\,dt\leq\mathcal{C}.$

By Poincare’s inequality (165) and Sobolev embedding, we also have

	$\displaystyle\rho^{-1}\in L^{2}(H^{1}),$		(29)
	$\displaystyle\rho^{-1}\in L^{2}(L^{6}).$

Consequently, we conclude that $\rho>0$ a.e. Furthermore, from (17) we have a bound on $\rho\log(T)$ . Hence, since $\rho>0$ a.e., we have that $T>0$ a.e.

Remark 3.

Being able to establish the absence of large vacuum regions is a remarkable feature of (4) and it is a key property for the existence proof below. Note that it is a consequence of the mass diffusion.

We can also obtain a useful estimate on temperature. By using the gas law, $p=\rho RT$ , we have $T^{1/2}\sim p^{1/2}\rho^{-1/2}$ . Since both $p^{1/2}$ and $\rho^{-1/2}$ are bounded in $L^{\infty}_{2}$ by (13) and (28), we have

\displaystyle\sqrt{T}\in L^{\infty}_{1}.

(30)

3.4 Estimates from the continuity equation

Multiplying (4a) by $\rho$ and integrating in space, using conditions (5) and (6), yield

\displaystyle\frac{1}{2}d_{t}\|\rho\|_{2}^{2}-\int_{\Omega}\sqrt{\rho}\nabla\rho\cdot\sqrt{\rho}{\bf v}\,d{\bf x}=-\int_{\Omega}\nu(\nabla\rho)^{T}(\nabla\rho)\,d{\bf x}.

Using Young’s inequality, we obtain,

\displaystyle\frac{1}{2}d_{t}\|\rho\|_{2}^{2}\leq-\int_{\Omega}\left(\sqrt{\nu}|\nabla\rho|\right)^{2}\,d{\bf x}+\int_{\Omega}\sqrt{\eta}(\sqrt{\rho}|\nabla\rho|)^{2}+\eta^{-1/2}(\sqrt{\rho}|{\bf v}|)^{2}\,d{\bf x}.

Choose, $\eta\leq\mu_{1}/2$ and use $\sqrt{\rho}{\bf v}\in L^{\infty}_{2}$ from (13) to obtain the bounds,

	$\displaystyle\sup_{t\in[0,\mathcal{T}]}\\|\rho\\|_{2}^{2}\leq\mathcal{C}.$		(31)
	$\displaystyle\int_{0}^{T}\int_{\Omega}\left(\sqrt{\nu}\|\nabla\rho\|\right)^{2}\,d{\bf x}<\mathcal{C}.$

By (31), Poincare, Sobolev embedding, and a standard interpolation inequality (see Appendix I, eqn (161)), we have,

\displaystyle\rho\in L^{\infty}_{2},\quad\quad\rho\in L^{2}(H^{1}),\quad\quad\rho\in L^{2}_{6},\quad\quad\rho\in L^{10/3}_{10/3}.

(32)

With (32) at hand, even stronger estimates on density are obtained by testing the continuity equation against $\rho^{2}$ ,

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}+\int_{\Omega}\rho^{2}(\nabla\cdot(\rho{\bf v}))\,d{\bf x}=\int_{\Omega}\rho^{2}\,\nabla\cdot(\nu\,\nabla\rho)\,d{\bf x}.

Integrating by parts with the use of the boundary conditions (5) and (6), results in

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}-\int_{\Omega}\nabla\rho^{2}\cdot(\rho{\bf v})\,d{\bf x}=-\int_{\Omega}\nabla\rho^{2}\cdot(\nu\,\nabla\rho)\,d{\bf x},

and

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}-\int_{\Omega}(\frac{2}{3}(\nabla\rho^{3})\cdot{\bf v})\,d{\bf x}=-\int_{\Omega}(2\rho\nabla\rho)\cdot(\nu\,\nabla\rho)\,d{\bf x}.

Partial integration and no-slip, lead to

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}+\int_{\Omega}(\frac{2}{3}\rho^{3}\nabla\cdot{\bf v})\,d{\bf x}=-\int_{\Omega}(2\rho\nabla\rho)\cdot(\nu\,\nabla\rho)\,d{\bf x}.

Using Young’s inequality, we split the term, $\rho^{3}\nabla\cdot{\bf v}<\eta\rho^{4}+\eta^{-1}(\rho\nabla\cdot{\bf v})^{2}$ , $\eta>0$ and use that the last term is bounded by (24). We have

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}\leq\int_{\Omega}\eta\frac{2}{3}\rho^{4}d{\bf x}-\int_{\Omega}(2\rho\nabla\rho)\cdot(\nu\,\nabla\rho)\,d{\bf x}+\mathcal{C},

and using (3),

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}\leq\int_{\Omega}\eta\frac{2}{3}\rho^{4}d{\bf x}-\int_{\Omega}(2\rho\nabla\rho)\cdot((\mu_{0}\rho^{-1}+\mu_{1}\rho))\,\nabla\rho)\,d{\bf x}+\mathcal{C}.

Now we intend to use (31) and the generalised Poincare inequality (165). Hence, we add and subtract the bounded term $\|\rho\|_{2}^{2}$ to the right-hand side.

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}-\int_{\Omega}\eta\frac{2}{3}\rho^{4}d{\bf x}+\int_{\Omega}(\mu_{1}|\nabla\rho^{2}|^{2})+2\mu_{0}\,|\nabla\rho|^{2}\,d{\bf x}+\|\rho\|_{2}^{2}-\|\rho\|_{2}^{2}\leq\mathcal{C}.

Next, we use the generalised Poincare inequality to conclude that, $\epsilon_{1}\|\rho^{2}\|_{2}^{2}\leq\|\rho\|_{2}^{2}+\int_{\Omega}\frac{1}{2}(\mu_{1}|\nabla\rho^{2}|^{2})d{\bf x}$ for some $\epsilon_{1}>0$ , such that

\displaystyle\frac{1}{3}d_{t}\|\rho\|_{3}^{3}-\int_{\Omega}\eta\frac{2}{3}\rho^{4}d{\bf x}+\epsilon_{1}\|\rho^{2}\|_{2}^{2}+\int_{\Omega}\frac{1}{2}(\mu_{1}|\nabla\rho^{2}|^{2})+2\mu_{0}\,|\nabla\rho|^{2}\,d{\bf x}\leq\|\rho\|_{2}^{2}+\mathcal{C}.

By choosing $0<\eta<\epsilon_{1}$ , we obtain the bounds,

\displaystyle\rho\in L^{\infty}_{3},\quad\quad\rho^{2}\in L^{2}(H^{1}),\quad\quad\rho\in L^{4}_{12},\quad\quad\rho\in L^{6}_{6}.

(33)

(The last bound by (164).)

Given these bounds, we can get even better bounds. Multiply the continuity equation by $\rho^{3}$ to obtain,

\displaystyle\int_{0,\Omega}^{\mathcal{T}}\rho^{3}\rho_{t}\,d{\bf x}\,dt+\int_{0,\Omega}^{\mathcal{T}}\rho^{3}(\nabla\cdot\rho{\bf v})\,d{\bf x}\,dt=\int_{0,\Omega}^{\mathcal{T}}\rho^{3}\nabla\cdot(\nu\nabla\rho)\,d{\bf x}\,dt.

(34)

We focus on the second integral. By repeatedly integrating by parts and using the boundary conditions (5) and (6), we have

	$\displaystyle\int_{0,\Omega}^{\mathcal{T}}\rho^{3}(\nabla\cdot\rho{\bf v})\,d{\bf x}\,dt=-\int_{0,\Omega}^{\mathcal{T}}(\nabla\rho^{3})\rho\cdot{\bf v}\,d{\bf x}\,dt$	$\displaystyle=$
	$\displaystyle-\int_{0,\Omega}^{\mathcal{T}}\frac{3}{4}(\nabla\rho^{4}\cdot{\bf v})\,d{\bf x}\,dt=\int_{0,\Omega}^{\mathcal{T}}\frac{3}{4}(\rho^{3}\cdot\rho\nabla{\bf v})\,d{\bf x}\,dt$	$\displaystyle\leq$
	$\displaystyle\int_{0,\Omega}^{\mathcal{T}}(\frac{3}{4}\rho^{3})^{2}+(\rho\nabla\cdot{\bf v})^{2}\,d{\bf x}\,dt.$

These integrals are controlled by (24) and (33). From (34), we thus obtain the following estimates:

\displaystyle\rho\in L^{\infty}_{4},\quad\quad\rho^{3/2}\nabla\rho\in L^{2}_{2},\quad\quad\rho\in L^{5}_{15},\quad\quad\rho\in L^{23/3}_{23/3}.

(35)

Remark 4.

It is possible to bootstrap further and obtain a little better integrability of density but the estimates above will serve the current purposes.

3.5 Renormalised internal energy

By integrating (12) in space and employing the boundary conditions, we obtain,

	$\displaystyle\int_{\Omega}c_{v}(\rho H)_{t}\,d{\bf x}=-\int_{\Omega}c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}\,d{\bf x}+\int_{\Omega}\left(H^{\prime}p\nabla\cdot{\bf v}\ +H^{\prime}\nu\rho\|\nabla{\bf v}\|^{2}\right)\,d{\bf x}$
	$\displaystyle-\int_{\Omega}H^{\prime\prime}4\kappa_{r}T^{3}\|\nabla T\|^{2}\,d{\bf x},$

which we integrate in time to obtain

	$\displaystyle-\int_{\Omega,0}^{\mathcal{T}}H^{\prime\prime}4\kappa_{r}T^{3}\|\nabla T\|^{2}\,d{\bf x}\,dt-\int_{\Omega,0}^{\mathcal{T}}c_{v}\rho\nu H^{\prime\prime}\|\nabla T\|^{2}\,d{\bf x}\,dt$	$\displaystyle=$
	$\displaystyle-\int_{\Omega}c_{v}(\rho H)\|_{t=\mathcal{T}}\,d{\bf x}+\int_{\Omega}c_{v}(\rho H)\|_{t=0}\,d{\bf x}+\int_{\Omega,0}^{\mathcal{T}}\left(H^{\prime}p\nabla\cdot{\bf v}\ +H^{\prime}\nu\rho\|\nabla{\bf v}\|^{2}\right)\,d{\bf x}\,dt.$			(36)

The particular choice, $H(T)=(1+T)^{1-\omega}$ , $0<\omega<1$ , implies:

•

$H(T)<T+1$ .
•

$0<H^{\prime}(T)=(1-\omega)(1+T)^{-\omega}<(1-\omega)$ for all $T\geq 0$ .
•

$H^{\prime\prime}(T)=-\omega(1-\omega)(1+T)^{-\omega-1}<0$ for all $T>0$ .

Hence, the left-hand side is positive. The second integral on the right-hand side of (36) is bounded by initial data and the first by (13). (That is, use the $L^{\infty}_{1}$ estimates of $p$ and $\rho$ and note that $\rho H\lesssim\rho+p$ .) Finally, the third integral on the right-hand side can be bounded as in (23) and by (24).

We conclude that the left-hand side is bounded and the strongest bound is obtained from the first integral,

\displaystyle\int_{\Omega,0}^{\mathcal{T}}(1+T)^{-\omega-1}T^{3}|\nabla T|^{2}\,d{\bf x}\,dt\leq\mathcal{C}.

From this bound, we have, $T^{1-\epsilon}\nabla T\in L^{2}_{2}$ , or $\nabla(T^{2-\epsilon})\in L^{2}_{2}$ . By (21), the generalised Poincare inequality and Sobolev embedding, we have $T^{2-\epsilon}\in L^{2}_{6}$ for any $\epsilon>0$ . Here, we recast the bound as,

\displaystyle\sqrt{T}\in L^{8-\epsilon}_{24-\epsilon}\quad\textrm{ for any}\quad\epsilon>0.

(37)

Next, we need the following interpolation inequality (see Appendix I),

\displaystyle\|u\|^{r}_{r}\leq\|u\|^{r-r\theta}_{p}\|u\|^{r\theta}_{q},

where $1=(r-r\theta)/p+r\theta/q$ . We intend to use (37) and (30). Therefore, we choose $r\theta=8$ , $q=24$ and $p=1$ and obtain $r=8+2/3$ . We conclude that $\sqrt{T}\in L^{8+2/3-\epsilon}_{8+2/3-\epsilon}$ . By choosing $\omega$ sufficiently small, we get $\delta=2/3-\epsilon>0$ such that,

\displaystyle T^{4}\in L^{1+\delta}_{1+\delta}\quad\textrm{for}\quad 1/12>\delta>0.

(38)

3.6 Improved estimates

By (21), we have $T\in L^{3}_{3}$ , which we use in

\displaystyle\int_{0}^{\mathcal{T}}\|p\|_{2}^{2}\,dt=R^{2}\int_{0,\Omega}^{\mathcal{T}}(\rho T)^{2}\,d{\bf x}\,dt\leq\mathcal{C}\int_{0}^{\mathcal{T}}\|\rho^{2}\|_{3}\|T^{2}\|_{3/2}\,dt\leq\mathcal{C}\int_{0}^{\mathcal{T}}\|\rho^{2}\|_{3}^{3}+\|T^{2}\|_{3/2}^{3/2}\,dt.

Since $L^{6}_{6}$ is embedded in $L^{23/3}_{23/3}$ , the density term is bounded by (35). Hence,

\displaystyle p\in L^{2}_{2}.

(39)

However, using the full integrability of $\rho$ and $T(\in L^{4}_{4})$ we have that

\displaystyle p\in L^{2+\delta}_{2+\delta}.

(40)

for some $\delta>0$ .

Next, we turn to velocity,

\displaystyle\|{\bf v}\|_{1}\leq\|\rho^{-1/2}\|_{2}^{2}+\|\sqrt{\rho}|{\bf v}|\|_{2}^{2}.

The right-hand side is bounded by (28) and (13), and we have,

\displaystyle{\bf v}\in L^{\infty}_{1}.

(41)

Next, we use Nash’ inequality,

\displaystyle\|u\|_{2}^{1+2/n}\leq\mathcal{C}\|u\|_{1}^{2/n}\|Du\|_{2},

(42)

where $u$ is a function, $n$ is the number of spatial dimensions, and $D$ the differential operator. Applying (42) with $n=3$ to velocity gives,

\displaystyle\int_{0}^{\mathcal{T}}(\|{\bf v}\|_{2}^{5/3})^{2}\,dt\leq\mathcal{C}\int_{0}^{\mathcal{T}}\|{\bf v}\|_{1}^{4/3}\|\nabla{\bf v}\|^{2}_{2}\,dt\leq\mathcal{C}\sup_{t\in[0,\mathcal{T}]}\|{\bf v}\|_{1}^{4/3}\int_{0}^{\mathcal{T}}\|\nabla{\bf v}\|^{2}_{2}\,dt,

(43)

which is bounded by (41) and (25).

Next, we apply the interpolation inequality $\|u\|_{3}\leq\|u\|_{1}^{1/5}\|u\|_{6}^{4/5}$ to obtain,

\displaystyle\int_{0}^{\mathcal{T}}\|{\bf v}\|_{3}^{10/4}\,dt\leq\int_{0}^{\mathcal{T}}\|{\bf v}\|_{1}^{10/4\cdot 1/5}\|{\bf v}\|_{6}^{2}\,dt,

(44)

which is bounded by (41) and (25). We summarise (44) and (43),

\displaystyle{\bf v}\in L^{10/3}_{2},\quad\quad{\bf v}\in L^{10/4}_{3}.

(45)

Turning to momentum, we consider $\|\rho^{2}|{\bf v}|^{2}\|_{1}\leq\|\rho^{2}\|_{3/2}\|{\bf v}\|_{6}^{2}$ , and by (25) and (33), we conclude that

\displaystyle\rho{\bf v}\in L^{2}_{2}.

(46)

Next, we will improve (46) slightly:

\displaystyle(\rho^{2}u^{2})^{1+\delta}=\rho^{2(1+\delta)}u^{2\epsilon(1+\delta)}u^{2(1-\epsilon)(1+\delta)}.

This yields

\displaystyle\|(\rho^{2}u^{2})^{1+\delta}\|_{1}\leq\|\rho^{2(1+\delta)}u^{2\epsilon(1+\delta)}\|_{3/2}\|u^{2(1-\epsilon)(1+\delta)}\|_{3}.

To bound the last factor, we require that $2(1-\epsilon)(1+\delta)\leq 2$ . This is satisfied for $\delta\leq\epsilon$ . Set $r=2(1+\delta)$ and consider the remaining factor

	$\displaystyle\\|\rho^{r}u^{r\epsilon}\\|_{3/2}=\\|\rho^{r-\epsilon r/2}\rho^{\epsilon r/2}u^{r\epsilon}\\|_{3/2}=\\|\rho^{r-\epsilon r/2}(\sqrt{\rho}u)^{\epsilon r}\\|_{3/2}$	$\displaystyle\leq$
	$\displaystyle\\|\rho^{\frac{3}{2}(r-\epsilon r/2)}\\|_{q}\\|(\sqrt{\rho}u)^{\frac{3}{2}\epsilon r}\\|_{p}$	.

Set $\frac{3}{2}\epsilon rp=2$ , i.e. $p=\frac{2}{3\epsilon(1+\delta)}$ . It is clear that we can choose $\epsilon>0$ so small, while keeping $\delta<\epsilon$ , such that $\frac{3}{2}(r-\epsilon r/2)q\leq 4$ . (We use (13),(35) and (25) to control all factors.) Hence,

\displaystyle\rho u\leq L^{2+\delta}_{2+\delta},\quad\textrm{for some}\quad\delta>0.

(47)

Furthermore, by $\|(\rho|{\bf v}|)^{8/5}\|\leq\|\sqrt{\rho}^{8/5}\|_{5}\|(\sqrt{\rho}|{\bf v}|)^{8/5}\|_{5/4}=\|\sqrt{\rho}\|^{8/5}_{8}\|(\sqrt{\rho}|{\bf v}|)\|^{8/5}_{2}\leq\|\sqrt{\rho}\|^{16/5}_{8}+\|(\sqrt{\rho}|{\bf v}|)\|^{16/5}_{2}$ . Using (35) and (13), we have

\displaystyle\rho{\bf v}\in L^{\infty}_{8/5}.

(48)

Moreover, using (48) together with ((45)) and the generalised Hölder inequality imply,

\displaystyle\|\rho({\bf v}\otimes{\bf v})\|_{24/23}\leq\mathcal{C}\|\rho{\bf v}\|_{8/5}\|{\bf v}\|_{3},

for some bounded constant $\mathcal{C}$ , such that

\displaystyle\rho({\bf v}\otimes{\bf v})\in L^{10/4}_{24/23}.

(49)

Next, we turn to the terms appearing in the energy flux. We will need integrability a little better than $L^{1}_{1}$ and therefore we take $\delta,\epsilon>0$ and consider

	$\displaystyle\\|\rho\|{\bf v}\|\|{\bf v}\|^{2}\\|_{1+\delta}^{1+\delta}=\\|(\rho\|{\bf v}\|\|{\bf v}\|^{2})^{1+\delta}\\|$	$\displaystyle=$
	$\displaystyle\\|((\sqrt{\rho})^{1-\epsilon}(\sqrt{\rho})^{1+\epsilon}\|{\bf v}\|^{1+\epsilon}\|{\bf v}\|^{2-\epsilon})^{1+\delta}\\|$	$\displaystyle=$
	$\displaystyle\\|(\sqrt{\rho})^{(1+\delta)(1-\epsilon)}(\sqrt{\rho}\|{\bf v}\|)^{(1+\delta)(1+\epsilon)}\|{\bf v}\|^{(2-\epsilon)(1+\delta)}\\|$	$\displaystyle\leq$
	$\displaystyle\\|(\sqrt{\rho})^{(1+\delta)(1-\epsilon)}(\sqrt{\rho}\|{\bf v}\|)^{(1+\delta)(1+\epsilon)}\\|_{3/2}\\|\|{\bf v}\|^{(2-\epsilon)(1+\delta)}\\|_{3}$	.

To bound the last factor in $L^{2}_{6}$ using (25), we demand that $(2-\epsilon)(1+\delta)\leq 2$ . We proceed with the choice $\epsilon=1/100$ and $\delta=1/200$ .

\displaystyle\|\rho|{\bf v}||{\bf v}|^{2}\|_{1+\delta}^{1+\delta}

\displaystyle\leq\|(\sqrt{\rho})^{197/198}(\sqrt{\rho}|{\bf v}|)^{607/598}\|_{3/2}\||{\bf v}|\|_{r}^{s},

where $0<r<6$ and $0<s<2$ . Next, we estimate

	$\displaystyle\\|(\sqrt{\rho})^{(197/198)\frac{3}{2}}(\sqrt{\rho}\|{\bf v}\|)^{(607/598)\frac{3}{2}}\\|_{1}$	$\displaystyle\leq$
	$\displaystyle\\|(\sqrt{\rho})^{(197/198)\frac{3}{2}}\\|_{p}\\|(\sqrt{\rho}\|{\bf v}\|)^{(607/598)\frac{3}{2}}\\|_{q}$	,		(50)

and choose $\frac{607}{598}\frac{3}{2}q=2$ , i.e., $q=2392/1821$ and $p=2392/571$ . The second factor in (50) is controlled by (13). The first is,

\displaystyle\|(\sqrt{\rho})^{(197/198)\frac{3}{2}}\|_{2392/571}\leq\|(\sqrt{\rho})\|^{197/132}_{52973/8473}\leq\|\sqrt{\rho}\|_{8}^{\infty}=\|\rho\|^{\infty}_{4},

where the right-hand side is controlled by (35). Hence,

\displaystyle\rho|{\bf v}|^{2}{\bf v}\in L^{1+1/200}_{1+1/200}.

(51)

We also need the following standard result: $\int_{\Omega,0}^{\mathcal{T}}|\log\rho|^{p}\,d{\bf x}\,dt\leq C_{p}\left(\|\rho\|_{2}^{2}+\|\rho^{-1}\|_{2}^{2}\right),$ where $C_{p}$ is a bounded constant for any $p<\infty$ . Hence,

\displaystyle\log\rho\in L^{p}_{p},

(52)

for any arbitrary $p$ , since $\rho,\rho^{-1}\in L^{2}_{2}$ by ((29) and (32))

3.7 Compactness of $\rho,{\bf v},T$

Assuming that we have a sequence of approximate solutions satisfying the a priori estimates, we will demonstrate that the primary variables converge a.e. (Such a sequence will be generated by means of a numerical scheme below.)

Density: We intend to apply Aubin-Lions Lemma (see Lemma 2 in Appendix I) to the continuity equation. From (32), we know that $\rho\in L^{2}(H^{1})$ , and $H^{1}$ is compactly embedded i $L^{2}$ .

We need a that $\rho_{t}$ is bounded in a space $X$ , such that $L^{2}$ is continuously embedded in $X$ . We try the negative Sobolev space $W^{-3,2}$ which is the set of all distributions, $u$ , such that $\int_{\Omega}u\phi\,d{\bf x}$ is bounded for all $\phi\in W^{3,2}$ .

Remark 5.

$L^{2}$ is continuously embedded in $W^{-3,2}$ if $\|u\|_{W^{-3,2}}\leq\mathcal{C}\|u\|_{L^{2}}$ , $u\in L^{2}$ , Here, $\|u\|_{W^{-3,2}}=\sup_{\phi\in W^{3,2}(\Omega),\|\phi\|_{3,2}=1}\frac{|<\phi,u>|}{\|\phi\|_{3,2}}$ . Since $|<\phi,u>|$ is bounded for all $u\in L^{2}$ , we conclude that $L^{2}$ functions, are bounded in $W^{-3,2}$ .

We proceed with the test functions $\phi\in W^{3,2}$ , for which $\nabla\phi\in L^{\infty}$ , and use them to test if $\rho_{t}\in L^{2}(W^{-3,2})$ .

	$\displaystyle\|<\partial_{t}\phi,\rho>\|=\|-<\nabla\phi,\nu\nabla\rho>-<\rho{\bf v},\nabla\phi>\|\leq$
	$\displaystyle\\|\nu\nabla\rho\\|_{1}\\|\nabla\phi\\|_{\infty}+\\|\rho{\bf v}\\|_{1}\\|\nabla\phi\\|_{\infty}.$

Furthermore, we obtain

	$\displaystyle\int_{0}^{\mathcal{T}}\sup_{\phi\in W^{3,2},\\|\phi\\|=1}\|<\partial_{t}\phi,\rho>\,dt\|^{2}$	$\displaystyle\leq$
	$\displaystyle\int_{0}^{\mathcal{T}}\sup_{\phi\in W^{3,2},\\|\phi\\|=1}\left(\\|\nu\nabla\rho\\|_{1}\\|\nabla\phi\\|_{\infty}+\\|\rho{\bf v}\\|_{1}\\|\nabla\phi\\|_{\infty}\right)^{2}dt$	.

First, we note that $\|\nabla\phi\|_{\infty}\leq\mathcal{C}\|\phi\|_{3,2}(\leq\mathcal{C})$ . Furthermore, (46), (35)/(33) and (18) bounds the fluxes on the right-hand side in $L^{2}_{2}$ and we can use Aubin-Lions Lemma, to conclude (by Sobolev embedding) that $\rho$ is compact in $L^{10/3-\epsilon}_{10/3-\epsilon}$ . Hence, a subsequence converges a.e. and since $\rho\in L^{23/3}_{23/3}$ by (35), we have strong convergence of $\rho$ in $L^{r}_{r}$ for $r<23/3$ .

Temperature: We assume that we have two sequences of approximate solutions $\rho^{n},T^{n}$ satisfying the a priori estimates obtained from the equation: $\rho^{n}\in L^{23/2}_{23/2};\quad T^{n}\in L^{4+\epsilon}_{4+\epsilon}$ (see (35) and (38)).

We have strong convergence of $\rho_{n}$ in $L^{p}_{p}$ , $p<23/3$ and (subsequential) weak convergence of $T_{n}$ in $L^{4}_{4}$ .

From this, we have weak convergence of $\rho^{n}T^{n}$ in $L^{2}_{2}$ since for any $\phi\in L^{2}_{2}$

\displaystyle<\phi,\rho^{n}T^{n}-\rho T>=<\phi T^{n},\rho^{n}-\rho>+<\phi\rho,T^{n}-T>.

(53)

For the first term to vanish $\phi T^{n}$ must be bounded in $L^{23/20+\epsilon}_{23/20+\epsilon}$ , which it is. For the other term, we need $\rho\phi\in L^{4/3}_{4/3}$ , which also holds.

Next, consider convergence of the $L^{2}$ -norm:

\displaystyle\int_{0}^{T}\int_{\Omega}(\rho^{n}T^{n})^{2}-(\rho T)^{2}\,dx\,dt=\int_{0}^{T}\int_{\Omega}((\rho^{n})^{2}-\rho^{2})(T^{n})^{2}+\rho^{2}((T^{n})^{2}-T^{2})\,dx\,dt,

which approaches zero thanks to the strongly and weakly converging sequences.

Next, we consider convergence in norm:

	$\displaystyle\int_{0}^{\mathcal{T}}\\|\rho^{n}T^{n}-\rho T\\|_{2}^{2}\,dt$	$\displaystyle=$
	$\displaystyle\int_{0}^{\mathcal{T}}<\rho^{n}T^{n}-\rho T,\rho^{n}T^{n}-\rho T>\,dt$	$\displaystyle=$
	$\displaystyle\int_{0}^{\mathcal{T}}(\rho^{n})^{2}(T^{n})^{2}-\rho^{2}T^{2}-2(\rho^{n}T^{n}-\rho T)\rho T\,dt$	.

Convergence of the norm and the weak convergence in $L^{2}_{2}$ ensures that the last integral vanishes. Hence, we have strong convergence of $\rho T$ in $L^{2}_{2}$ and a.e. convergence of a subsequence. Since $\rho$ is a.e. convergent, there is a subsequence of $T_{n}$ that converges a.e.

Velocity: Here we use that the sequence $\sqrt{\rho^{n}}u^{n}$ is bounded in $L^{2}_{2}$ (13); we have $u^{n}$ bounded and weakly convergent in $L^{2.5}_{2.5}$ , (45); and $\sqrt{\rho^{n}}$ is strongly convergent up to $L^{p}_{p}$ , $p<46/3$ , thanks to (35).

Weak convergence of $\sqrt{\rho^{n}}u^{n}\in L^{2}_{2}$ requires that for $\phi\in L^{2}_{2}$ , we have $\sqrt{\rho^{n}}\phi\in L^{5/3}_{5/3}$ and $\phi u\in L^{46/43+\epsilon}_{46/43+\epsilon}$ . Both conditions can be confirmed with the available estimates.

Next, we check convergence of the $L^{2}_{2}$ norm:

\displaystyle\int_{0}^{\mathcal{T}}\int_{\Omega}\rho^{n}(u^{n})^{2}-\rho u^{2}\,dx\,dt=\int_{0}^{T}\int_{\Omega}(\rho^{n}-\rho)(u^{n})^{2}+\rho((u^{n})^{2}-u^{2})\,dx\,dt.

(54)

Since $(u^{n})^{2}\in L^{5/4}_{5/4}$ , we need $\rho$ to be strongly convergent in at least $L^{5}_{5}$ , which indeed we have.

Hence, one can in the same way as for temperature deduce strong convergence of $\sqrt{\rho^{n}}u^{n}$ in $L^{2}_{2}$ and draw an a.e. converging subsequence and get a.e. convergence of $u^{n}$ .

Remark 6.

It is also possible to prove that $\rho u$ is compact in $L^{3/2}_{3/2}$ via Aubin-Lions Lemma. The fluxes are bounded in $L^{1+\epsilon}_{1+\epsilon}$ which is enough to show that $\partial_{t}(\rho u)\in L^{1+\epsilon}(W^{-3/2})$ . Furthermore, $\rho{\bf v}$ and $\nabla\cdot(\rho{\bf v})$ are both in $L^{1+\epsilon}_{1+\epsilon}$ and $L^{3/2}$ is compactly embedded in $W^{1,1}$ . Hence, we have a.e. convergence of $\rho{\bf v}$ and therefore also of ${\bf v}$ .

Magnitude of velocity gradients:

We apply the analogous arguments as above to the sequence $\sqrt{\rho^{n}}\sqrt{|u^{n}_{x}|}$ . We have strong convergence of $\sqrt{\rho^{n}}\in L^{p<46/3}_{p<46/3}$ and weak convergence of $\sqrt{|u^{n}_{x}|}\in L^{4}_{4}$ . This is enough to get
$\lim_{n\rightarrow\infty}\int_{0}^{\mathcal{T}}\left(\|\sqrt{\rho^{n}}\sqrt{|u^{n}_{x}|}\|_{2}^{2}-\|\sqrt{\rho}\sqrt{|u_{x}|}\|_{2}^{2}\right)\,dt\rightarrow 0$ .

Weak convergence of $\sqrt{\rho^{n}}\sqrt{|u^{n}_{x}|}\in L^{2}_{2}$ can be verified with the help of the estimates. Hence, we can infer strong convergence of $\sqrt{\rho^{n}}\sqrt{|u^{n}_{x}|}\in L^{2}_{2}$ and draw an a.e. converging subsequence. We obtain that $|u^{n}_{x}|$ converges a.e. and strongly up to $L^{2-\epsilon}_{2-\epsilon}$ for $\epsilon>0$ .

Magnitude of temperature and density gradients: The argument can be mimicked for the sequences $\sqrt{\rho^{n}}\sqrt{|T^{n}_{x}|}$ and $\sqrt{\rho^{n}}\sqrt{|\rho^{n}_{x}|}$ that satisfy equally strong estimates.

Summary: $\rho,{\bf v},T$ and $|\nabla{\bf v}|,|\nabla T|,|\nabla\rho|$ are all converging a.e.

3.8 Weak sequential convergence

Assuming that we have a sequence of approximate solutions that satisfies the a priori estimates, we will demonstrate that the estimates are sufficient to infer convergence to an $L^{1}_{1}$ weak solution. To this end, we need that all inviscid fluxes are bounded in $L^{1+\epsilon}_{1+\epsilon}$ , $\epsilon>0$ and a.e. convergence of the principal variables. The convergence of the diffusive fluxes are discussed separately.

Throughout, we assume that initial data is appropriately bounded. For instance, all variables bounded in $L_{\infty}$ and $\rho,T>0$ and bounded away from zero, would do.

Continuity equation: We begin by testing the continuity equation against $\phi\in C_{0}^{1}(\bar{\Omega}\times[0,\mathcal{T}))$ ,

\displaystyle\int_{0}^{\mathcal{T}}\left(<\phi_{n},\rho_{t}>+<\phi,\nabla(\rho_{n}{\bf v}_{n})>\right)\,dt=\int_{0}^{\mathcal{T}}<\phi,\nabla\cdot(\nu_{n}\nabla\rho_{n})>\,dt.

Integration by parts results in

	$\displaystyle-<\phi,\rho_{n}>\|_{t=0}-\int_{0}^{\mathcal{T}}\left(<\partial_{t}\phi,\rho_{n}>-<\nabla\cdot\phi,\rho_{n}{\bf v}_{n}>\right)\,dt$	$\displaystyle=$
	$\displaystyle-\int_{0}^{\mathcal{T}}<\nabla\phi,(\nu_{n}\nabla\rho_{n})>\,dt$	,

which we require to hold for any $\phi\in C_{0}^{1}(\bar{\Omega}\times[0,\mathcal{T}))$ .

The first integral is well-defined thanks to the strong estimates of $\rho_{n}$ . In the second, the momentum is bounded in $L^{2}_{2}$ by (47) and both $\rho_{n}$ and ${\bf v}_{n}$ converge a.e. Hence, we have strong convergence of the inviscid flux term in $L^{2}_{2}$ . The last integral is controlled by (18) and (33). Since $\nu_{n}$ is a function of $\rho_{n}$ , it converges a.e.. The weak convergence in $L^{2}_{2}$ of the density gradient, allows us to conclude that the product converges weakly to the correct limit.

Momentum equations: The weak form is given by,

	$\displaystyle-<\phi,\rho_{n}{\bf v}_{n}>\|_{t=0}-\int_{0}^{\mathcal{T}}\left(<\partial_{t}\phi,\rho_{n}{\bf v}_{n}>-<\nabla\cdot\phi,(\rho_{n}({\bf v}_{n}\otimes{\bf v}_{n})+p_{n})\right)\,dt$	$\displaystyle=$
	$\displaystyle-\int_{0}^{\mathcal{T}}<\nabla\phi,(\nu_{n}\nabla(\rho_{n}{\bf v}_{n}))\,dt,$

where $\phi\in C_{0}^{1}(\Omega\times[0,\mathcal{T});\mathbb{R}^{3})$ .

The first integral is controlled by (46). The second integral is controlled by (49) which, along with a.e. convergence of $\rho_{n},{\bf v}_{n}$ , implies weak convergence in $L^{1+\epsilon}_{1+\epsilon}$ . Convergence of the pressure term follows from (40).

The integrand in the last term contains,

\displaystyle\rho(\nabla\rho)\cdot{\bf v},\quad\rho^{2}\nabla\cdot{\bf v},\quad\rho^{-1}(\nabla\rho)\cdot{\bf v},\quad\nabla\cdot{\bf v}.

In the first term, we have $\rho\nabla\rho\in L^{2}_{2}$ by (33). Hence, a subsequence converges weakly in $L^{2}_{2}$ . By a.e. convergence of ${\bf v}$ and ${\bf v}\in L^{10/4}_{3}$ by (45), we have strong convergence of ${\bf v}\in L^{p}_{p}$ , $1\leq p<2.5$ . Hence, the product converges weakly in $L^{1+\epsilon}_{1+\epsilon}$ for some $\epsilon>0$ .

The term $\rho^{2}\nabla\cdot{\bf v}$ is handled in the same way. A subsequence of $\nabla\cdot{\bf v}$ converges weakly by (25). By (35) and a.e. convergence, $\rho$ converges strongly in $\rho\in L^{23/3-\epsilon}_{23/3-\epsilon}$ . Hence, we have weak convergence of the second term in $L^{1+\epsilon}_{1+\epsilon}$ for some $\epsilon>0$ .

The third term is handled like the first noting that $\nabla\log\rho\in L^{2}_{2}$ by (18). Weak convergence of the fourth term is immediate (for a subsequence) by (25).

The total energy equation: Finally, we turn to the weak form of the energy equation:

	$\displaystyle-<\phi,E_{n}>\|_{t=0}-\int_{0}^{\mathcal{T}}\left(<\partial_{t}\phi,E_{n}>-<\nabla\phi,(E_{n}{\bf v}_{n}+p_{n}{\bf v}_{n})\right)\,dt$	$\displaystyle=$
	$\displaystyle-\int_{0}^{\mathcal{T}}<\nabla\phi,(\nu_{n}\nabla E_{n})\,dt-\int_{0}^{\mathcal{T}}<\nabla\phi,(\kappa_{r}\nabla T^{4}_{n})\,dt$	,

with test functions $\phi\in C_{0}^{2}(\Omega\times[0,\mathcal{T}))$ .

The first integral is handled in the same way as the inviscid flux in the momentum equation. Weak convergence of $p_{n}{\bf v}_{n}$ appearing in the second integral follows from (40) and (45). The $\rho|{\bf v}|^{2}{\bf v}$ term is bounded in $L^{1+\epsilon}_{1+\epsilon}$ according to (51). Furthermore, it converges a.e. Hence, it converges weakly in $L^{1}_{1}$ to the correct limit.

The radiation term is recast as $<\nabla\phi,\nabla T^{4}>=-<\Delta\phi,T^{4}>$ . Strong convergence (of a subsequence) is ensured by (38) and a.e. convergence.

The pressure term in the diffusive flux is handled as follows,

\displaystyle\nu\nabla p\sim(\rho+\rho^{-1})\nabla p\sim\rho T\nabla\rho+\rho^{2}\nabla T+T\nabla\log\rho+\nabla T.

These are easily handled by using (18), (32), (40) and (33), along with a.e. convergence of $\rho$ and $T$ , and noting that by combining (19) and (18), we get $\nabla T\in L^{2}_{2}$ .

Finally, we consider

	$\displaystyle<\nabla\phi,\nu_{n}\nabla(\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2})>=-<\nabla\cdot((\nabla\phi)\nu_{n}),(\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2})>$	$\displaystyle=$
	$\displaystyle-<\Delta\phi,\nu_{n}\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2}>-<\nabla\phi\cdot\nabla\nu_{n},(\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2})>$	.

Dropping the subscript we have terms of the form:

\displaystyle|{\bf v}|^{2},\quad\rho^{2}|{\bf v}|^{2},\quad\rho|{\bf v}|^{2}\nabla(\rho+\rho^{-1}).

(55)

Strong convergence of the first follows directly from ${\bf v}\in L^{2.5}_{2.5}$ (see (45)) and a.e convergence. In the second, we have that $\rho{\bf v}$ is converging strongly in $L^{2+\delta}_{2+\delta}$ (see (47)).

The fourth is,

\displaystyle-\nabla(\rho^{-1})\rho|{\bf v}|^{2}=(\nabla\log\rho)|{\bf v}|^{2}=\nabla(\log\rho|{\bf v}|^{2})-2\log\rho({\bf v}\cdot\nabla{\bf v}).

In the first term of the last expression, we can move the derivatives onto the test function. The bound $\log\rho|{\bf v}|^{2}$ follows since $\log\rho\in L^{p}_{p}$ for any arbitrary $p$ (see (52)) and $|{\bf v}|^{2}\in L^{5/4}_{5/4}$ .

We are left to bound $\log\rho({\bf v}\cdot\nabla{\bf v})$ . We have that $\nabla{\bf v}$ is converging weakly in $L^{2}_{2}$ . Hence, we need that $\log(\rho){\bf v}$ is converging strongly in $L^{2}_{2}$ . Since both $\rho$ and ${\bf v}$ are a.e. convergent, we only need a bound on the product in a little better space than $L^{2}_{2}$ . This follows since we have ${\bf v}\in L^{10/4}_{10/4}$ and $\log(\rho)\in L^{11}_{11}$ .

Finally, we must deduce convergence of the third term of (55):

\displaystyle(\nabla\rho)\rho|{\bf v}|^{2}=\frac{1}{2}(\nabla\rho^{2})|{\bf v}|^{2}=\frac{1}{2}\nabla(\rho^{2}|{\bf v}|^{2})-\rho^{2}{\bf v}\nabla{\bf v}.

(56)

In the first term, we can move the derivative to the test function and are left with $\rho^{2}|{\bf v}|^{2}$ which is bounded in $L^{1+\epsilon}_{1+\epsilon}$ and a.e. convergent. The last step is to prove convergence of

\displaystyle\rho^{2}{\bf v}\nabla{\bf v}=(\rho{\bf v})(\rho\nabla{\bf v}).

(57)

We have strong convergence of $(\rho{\bf v})\in L^{2+\delta}_{2+\delta}$ for some fixed $\delta>0$ . Hence, if we have weak convergence of $(\rho\nabla{\bf v})\in L^{2-\epsilon}_{2-\epsilon}$ for $\epsilon>0$ , (57) is converging weakly in $L^{1}_{1}$ .

By (24), we know that $(\rho^{n}\nabla{\bf v}^{n})\in L^{2}_{2}$ and hence weakly converging to some limit $(\rho\nabla{\bf v})^{*}$ . Using the strong and weak convergence of $\rho^{n}$ and $\nabla{\bf v}^{n}$ , we get that $\rho^{n}\nabla{\bf v}^{n}\rightharpoonup\rho\nabla{\bf v}\in L^{1}_{1}$ . That is when we test it against an $L^{\infty}_{\infty}$ function. Since the sequence converges to $(\rho{\bf v})^{*}$ for all test functions in $L^{2}_{2}$ including those that are also bounded in $L^{\infty}_{\infty}$ , we conclude that $(\rho\nabla{\bf v})^{*}=\rho\nabla{\bf v}$ by uniqueness of the weak limit.

The entropy inequality: $\rho S$ is convergent if $\rho\log(T)$ and $\rho\log(\rho)$ are both convergent in $L^{1+\epsilon}_{1+\epsilon}$ . This follows from $\log\rho$ and $\log T$ being bounded in $L^{2}(H^{1})$ , which together with the strong estimates of $\rho$ along with a.e. convergence, ensures weak convergence of $\rho S\in L^{1+\epsilon}_{1+\epsilon}$

Convergence of the entropy flux follows from the $\log\rho$ and $\log T$ bounds and (47), i.e., $\rho{\bf v}\in L^{2+\epsilon}_{2+\epsilon}$ and a.e. convergence of $\rho,T,{\bf v}$ . ( $\rho{\bf v}$ converges a.e. and hence strongly in $L^{2+\delta}_{2+\delta}$ for $\delta<\epsilon$ .)

Finally, the right-hand side of the entropy equation (7) is summable and non-positive ensuring that the entropy inequality, $(-\rho S)_{t}+\nabla\cdot(-\rho{\bf v}S)\leq 0$ , is satisfied in a distributional sense in $L^{1}_{1}$ .

However, it is challenging to verify these relations for the numerical scheme below, and in this treatise we will demonstrate that entropy is diffused globally. That is,

\displaystyle\int_{\Omega}-\rho S\,d{\bf x}\leq\int_{\Omega}-\rho S\,d{\bf x}|_{t=0}.

3.9 Weak entropy solution

The weak sequential compactness above, demonstrates that it should be possible to construct weak entropy solutions to (4). Herein, we will utilise a numerical scheme to produce a sequence of approximate solutions. We will prove that the sequence converges to a weak solution in the following sense:

Definition 1.

We say that the triple $\{\rho,{\bf v},T\}$ is a weak entropy solution to (4) with boundary conditions (5) and (6), if $\rho,T>0$ a.e. in $\Omega\times[0,\mathcal{T}]$ and the triple satisfies:

•

For all $\phi\in C^{1}_{0}([0,\mathcal{T})\times\bar{\Omega})$ , the continuity equation is satisfied weakly,

\displaystyle-<\phi,\rho>|_{t=0}-\int_{0}^{\mathcal{T}}\left(<\partial_{t}\phi,\rho>+<\nabla\phi,\rho{\bf v}>\right)\,dt

\displaystyle=-\int_{0}^{\mathcal{T}}<\nabla\phi,(\nu\nabla\rho)>\,dt.

•

For all $\phi\in C^{1}_{0}([0,\mathcal{T})\times\Omega)$ , the momentum equation is satisfied weakly,

	$\displaystyle-<\phi,\rho{\bf v}>\|_{t=0}-\int_{0}^{\mathcal{T}}\left(<\partial_{t}\phi,\rho{\bf v}>+<\nabla\phi,(\rho({\bf v}\otimes{\bf v})+p)\right)\,dt$	$\displaystyle=$
	$\displaystyle-\int_{0}^{\mathcal{T}}<\nabla\phi,\nu\nabla\rho{\bf v})\,dt.$

•

For all $\phi\in C^{2}_{0}([0,\mathcal{T})\times\Omega)$ , the internal energy equation is satisfied weakly,

	$\displaystyle-<\phi,E>\|_{t=0}-\int_{0}^{\mathcal{T}}\left(<\partial_{t}\phi,E>+<\nabla\phi,(E+p){\bf v}>\right)\,dt$	$\displaystyle=$
	$\displaystyle-\mathcal{D}+\int_{0}^{\mathcal{T}}<\Delta\phi,(\kappa_{r}T^{4})\,dt,$

where the diffusive term is interpreted as:

\displaystyle\mathcal{D}=<\nabla\phi,\nu\nabla E>=<\nabla\phi,\nu_{n}\nabla(\frac{1}{2}\rho_{n}|{\bf v}_{n}|^{2}+p)>,

and

	$\displaystyle<\nabla\phi,\nu_{n}\nabla(\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2})>$	$\displaystyle=$
	$\displaystyle-<\Delta\phi,\nu_{n}\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2}>-<\nabla\phi,(\nabla\nu_{n})\frac{1}{2}\rho_{n}\|{\bf v}_{n}\|^{2}>$	,

and

\displaystyle\mu_{1}<\nabla\phi,(\nabla\rho_{n})\frac{1}{2}\rho_{n}|{\bf v}_{n}|^{2}>=-\mu_{1}<\Delta\phi,\frac{1}{4}\rho_{n}^{2}|{\bf v}_{n}|^{2}>-\mu_{1}<\nabla\phi,\frac{1}{2}\rho_{n}^{2}{\bf v}_{n}\cdot\nabla{\bf v}>,

\displaystyle\mu_{0}<\nabla\phi,(\nabla\rho_{n}^{-1})\frac{1}{2}\rho_{n}|{\bf v}_{n}|^{2}>=\frac{\mu_{0}}{2}<\Delta\phi,\log\rho_{n}|{\bf v}_{n}|^{2}>+\mu_{0}<\nabla\phi,(\log\rho_{n}){\bf v}_{n}\cdot\nabla{\bf v}>.

•

The entropy inequality

$\displaystyle\int_{\Omega}-\rho S\,d{\bf x}\leq\int_{\Omega}-\rho S\,d{\bf x}|_{t=0},$

is satisfied weakly.

∎

4 The numerical approximation scheme

To generate a sequence of approximate solutions to (4), we use a finite volume scheme and consider the domain $\Omega=(0,1)^{3}$ . We discretise $\Omega$ with $N+1$ grid points in each direction: $x_{i}$ , $y_{j}$ , $z_{k}$ , $i,j,k=0...N$ . The discrete domain is denoted $\Omega_{h}$ . Furthermore, let $x_{i+1/2}=\frac{x_{i+1}+x_{i}}{2}$ and similarly for $y_{j+1/2},z_{k+1/2}$ .

Each grid point is associated with a control volume. In the interior of the domain, the control volumes are centred around the grid points: $\mathcal{V}_{ijk}=(x_{i-1/2},x_{i+1/2}]\times(y_{j-1/2},y_{j+1/2}]\times(z_{k-1/2},z_{k+1/2}]$ and, $V_{ijk}=|\mathcal{V}_{ijk}|$ denotes the measure of the control volume. At the boundaries, the control volumes are defined in the same way (towards other grid points) and are closed with the domain boundary. (We use the standard node-centred finite volume scheme on a Cartesian grid.)

The variables are constant within each control volume such that they form piecewise-constant functions in space. We adopt the same notation for the approximations that we used for the continuous equation. For instance, $p_{ijk}$ is the value of the piecewise constant pressure and ${\bf v}_{ijk}=(u_{ijk},v_{ijk},w_{ijk})$ is the velocity vector in $\mathcal{V}_{ijk}$ .

Furthermore, we need the control-volume-boundary areas

	$\displaystyle{\mathcal{S}}^{x}_{jk}$	$\displaystyle=(y_{j+1/2}-y_{j-1/2})(z_{k+1/2}-z_{k-1/2}),$
	$\displaystyle{\mathcal{S}}^{y}_{ik}$	$\displaystyle=(x_{i+1/2}-x_{i-1/2})(z_{k+1/2}-z_{k-1/2}),$
	$\displaystyle{\mathcal{S}}^{z}_{ij}$	$\displaystyle=(x_{i+1/2}-x_{i-1/2})(y_{j+1/2}-y_{j-1/2}),$

In Fig 1, a two-dimensional slice of control volumes is depicted (dashed lines). ${\mathcal{S}}^{x},{\mathcal{S}}^{y}$ that make up the boundary of the control volume are also indicated for the point $(x_{i},y_{j},z_{k})$ .

Refer to caption — Figure 1: A two-dimensional slice of a grid. $(i,j,k)$ denotes a point in the primal grid (solid lines) and the dashed lines indicate the dual grid (control volumes). $V_{ijk}$ , ${\mathcal{S}}^{x}_{jk},{\mathcal{S}}^{y}_{ik}$ are the control volume and control-volume-boundary areas associated with $(i,j,k)$ .

To simplify notation, we use the indices $-1/2$ and $N+1/2$ to denote the values at index $0$ and $N$ , respectively. For instance, $x_{-1/2}=x_{0}$ and $p_{N+1/2,j,k}=p_{N,j,k}$ .

We need the following notation for a scalar $a$ :

	$\displaystyle\Delta_{-}^{x}a_{I+1jk}$	$\displaystyle=\Delta_{+}^{x}a_{Ijk}=a_{I+1jk}-a_{Ijk},$
	$\displaystyle\Delta_{-}^{y}a_{iJ+1k}$	$\displaystyle=\Delta_{+}^{y}a_{iJk}=a_{iJ+1k}-a_{iJk},$
	$\displaystyle\Delta_{-}^{z}a_{ijK+1}$	$\displaystyle=\Delta_{+}^{z}a_{ijK}=a_{ijK+1}-a_{ijK},$

where $I,J,K$ signify either a “whole” or “half” index and the boundary-index convention $-1/2\rightarrow 0$ etc, is applied. The difference operators are given as

$\displaystyle D^{x}_{-}a_{I+1jk}$	$\displaystyle=\frac{\Delta^{x}_{-}a_{I+1jk}}{\Delta^{x}_{-}x_{I+1jk}},$
$\displaystyle D^{y}_{-}a_{iJ+1k}$	$\displaystyle=\frac{\Delta^{y}_{-}a_{iJ+1k}}{\Delta^{y}_{-}y_{iJ+1k}},$	(58)
$\displaystyle D^{z}_{-}a_{ijK+1}$	$\displaystyle=\frac{\Delta^{z}_{-}a_{ijK+1}}{\Delta^{z}_{-}z_{ijK+1}},$

and similarly for $D_{+}^{x,y,z}$ .

Furthermore, we need the two-point averages

$\displaystyle\bar{a}_{i+1/2jk}$	$\displaystyle=\frac{a_{ijk}+a_{i+1jk}}{2},$
$\displaystyle\hat{a}_{i+1/2jk}$	$\displaystyle=\left\{\begin{array}[]{cc}\frac{\Delta^{x}_{-}a_{i+1jk}}{\Delta_{-}^{x}\log a_{i+1jk}}&a_{ijk}\neq a_{i+1jk}\\ a_{ijk}&a_{ijk}=a_{i+1jk}\end{array}\right.\quad\textrm{(log mean)},$	(61)
$\displaystyle\check{a}_{i+1/2jk}$	$\displaystyle=\sqrt{a_{i+1jk}a_{ijk}}\quad\textrm{(geometric mean)},$

and for $a_{i}\geq 0$ , the following elementary relations hold:

\displaystyle\check{a}_{i+1/2}\leq\hat{a}_{i+1/2}\leq\bar{a}_{i+1/2}.

(62)

Remark 7.

In (61), the half index indicates in which direction the average is taken. See [IR09] for a stable numerical implementation of the log mean.

Remark 8.

The approximate derivative operators (58) are the same as the standard node-centred summation-by-parts (SBP) finite-volume scheme. (See [NFAE03].) In fact, since the grid is Cartesian, they coincide with the standard second-order SBP finite difference scheme. (See e.g. [SN14].)

The following relation between arithmetic averages is frequently needed:

\displaystyle\overline{ab}_{i+1/2}

\displaystyle=\bar{a}_{i+1/2}\bar{b}_{i+1/2}+\frac{\Delta_{+}a_{i}\Delta_{+}b_{i}}{4}.

(63)

The measure of the control volumes satisfies

\displaystyle V_{ijk}=\Delta_{-}^{x}x_{i+1/2}{\mathcal{S}}^{x}_{jk}=\Delta_{-}^{y}y_{j+1/2}{\mathcal{S}}^{y}_{ik}=\Delta_{-}^{z}z_{k+1/2}{\mathcal{S}}^{z}_{ik},

(64)

and we define the $L^{2}$ -equivalent norm,

\displaystyle\|a\|^{2}=\sum_{ijk=0}^{N}V_{ijk}a_{ijk}^{2}.

Throughout, we assume a uniform refinement of the grid such that $\Delta^{x}_{-}x_{i+1/2},\Delta^{y}_{-}y_{j+1/2},\Delta^{z}_{-}z_{k+1/2}$ are less than $h$ for all $i,j,k$ . In this notation, all ${\mathcal{S}}\sim\mathcal{O}(h^{2})$ and all $V\sim\mathcal{O}(h^{3})$ .

The discrete $L^{2}$ -norms of the derivatives are given by

	$\displaystyle\\|D^{x}_{+}a\\|^{2}=\sum_{ijk=0}^{N-1,N,N}V_{ijk}(D_{+}^{x}a_{ijk})^{2},$
	$\displaystyle\\|D^{y}_{+}a\\|^{2}=\sum_{ijk=0}^{N,N-1,N}V_{ijk}(D_{+}^{y}a_{ijk})^{2},$		(65)
	$\displaystyle\\|D^{z}_{+}a\\|^{2}=\sum_{ijk=0}^{N,N,N-1}V_{ijk}(D_{+}^{z}a_{ijk})^{2}.$

4.1 The numerical scheme

The scheme is given by,

	$\displaystyle V_{ijk}({\bf u}_{ijk})_{t}+{\mathcal{S}}^{x}_{jk}\Delta^{x}_{-}{\bf f}_{i+1/2jk}+{\mathcal{S}}^{y}_{ik}\Delta^{y}_{-}{\bf g}_{ij+1/2k}+{\mathcal{S}}^{z}_{ij}\Delta^{z}_{-}{\bf h}_{ijk+1/2}$	$\displaystyle=0$		(66)
	$\displaystyle 0\leq i,j,k,\leq N,$

where ${\bf u}_{ijk}=(\rho_{ijk},\rho_{ijk}{\bf v}_{ijk},E_{ijk})^{T}$ . The fluxes consist of a convective and a diffusive part, ${\bf f}_{i+1/2jk}={\bf f}^{c}_{i+1/2jk}-{\bf f}^{d}_{i+1/2jk}$ , and similarly for ${\bf g},{\bf h}$ . (These fluxes will be defined below.)

Remark 9.

With the operators defined above, we can equivalently state the scheme on finite difference form as,

	$\displaystyle({\bf u}_{ijk})_{t}+D_{-}^{z}{\bf f}_{i+1/2jk}+D_{-}^{y}{\bf g}_{ij+1/2k}+D_{-}^{z}{\bf h}_{ijk+1/2}$	$\displaystyle=0$		(67)
	$\displaystyle 0\leq i,j,k,\leq N$

To define the fluxes, we use the notation:

	$\displaystyle\overline{\|{\bf v}\|^{2}}_{i+1/2jk}$	$\displaystyle=(\overline{u^{2}})_{i+1/2jk}+(\overline{v^{2}})_{i+1/2jk}+(\overline{w^{2}})_{i+1/2jk},$
	$\displaystyle\overline{\|{\bf v}\|}^{2}_{i+1/2jk}$	$\displaystyle=(\bar{u}_{i+1/2jk})^{2}+(\bar{v}_{i+1/2jk})^{2}+(\bar{w}_{i+1/2jk})^{2}\quad\textrm{and analogously for y,z-means},$

The convective fluxes (defined below) are inspired by the entropy conservative and kinetic energy preserving flux proposed by Chandrashekar in [Cha13] and we adopt his notation $\beta=\frac{1}{2RT}=\frac{\rho}{2p}$ .

For $i=0...N-1$ and $j,k=0...N$ :

$\displaystyle p_{i+1/2jk}=$	$\displaystyle\frac{\bar{\rho}_{i+1/2jk}}{2\bar{\beta}_{i+1/2jk}}$
$\displaystyle{\bf f}^{c,1}_{i+1/2jk}=$	$\displaystyle\overline{\rho u}_{i+1/2jk}$
$\displaystyle{\bf f}^{c,2}_{i+1/2jk}=$	$\displaystyle\bar{u}_{i+1/2jk}\overline{\rho u}_{i+1/2jk}+p_{i+1/2jk}$
$\displaystyle{\bf f}^{c,3}_{i+1/2jk}=$	$\displaystyle\bar{v}_{i+1/2jk}\overline{\rho u}_{i+1/2jk}$	(68)
$\displaystyle{\bf f}^{c,4}_{i+1/2jk}=$	$\displaystyle\bar{w}_{i+1/2jk}\overline{\rho u}_{i+1/2jk}$
$\displaystyle{\bf f}^{c,5}_{i+1/2jk}=$	$\displaystyle\frac{1}{2(\gamma-1)\hat{\beta}_{i+1/2jk}}\overline{\rho u}_{i+1/2jk}-\frac{\overline{\|{\bf v}\|^{2}}_{i+1/2jk}}{2}\overline{\rho u}_{i+1/2jk}$
	$\displaystyle+\overline{\|{\bf v}\|}^{2}_{i+1/2jk}\overline{\rho u}_{i+1/2jk}+p_{i+1/2jk}\bar{u}_{i+1/2jk},$

For $j=0...N-1$ and $i,k=0...N$ :

$\displaystyle p_{ij+1/2k}=$	$\displaystyle\frac{\bar{\rho}_{ij+1/2k}}{2\bar{\beta}_{ij+1/2k}}$
$\displaystyle{\bf g}^{c,1}_{ij+1/2k}=$	$\displaystyle\overline{\rho v}_{ij+1/2k}$
$\displaystyle{\bf g}^{c,2}_{ij+1/2k}=$	$\displaystyle\bar{u}_{ij+1/2k}\overline{\rho v}_{ij+1/2k}$
$\displaystyle{\bf g}^{c,3}_{ij+1/2k}=$	$\displaystyle\bar{v}_{ij+1/2k}\overline{\rho v}_{ij+1/2k}+p_{ij+1/2k}$	(69)
$\displaystyle{\bf g}^{c,4}_{ij+1/2k}=$	$\displaystyle\bar{w}_{ij+1/2k}\overline{\rho v}_{ij+1/2k}$
$\displaystyle{\bf g}^{c,5}_{ij+1/2k}=$	$\displaystyle\frac{1}{2(\gamma-1)\hat{\beta}_{ij+1/2k}}\overline{\rho v}_{ij+1/2k}-\frac{\overline{\|{\bf v}\|^{2}}_{ij+1/2k}}{2}\overline{\rho v}_{ij+1/2k}$
	$\displaystyle+\overline{\|{\bf v}\|}^{2}_{ij+1/2k}\overline{\rho v}_{ij+1/2k}+p_{i+1/2jk}\bar{v}_{ij+1/2k}.$

For $k=0...N-1$ and $i,j=0...N$ :

$\displaystyle p_{ijk+1/2}=$	$\displaystyle\frac{\bar{\rho}_{ijk+1/2}}{2\bar{\beta}_{ijk+1/2}}$
$\displaystyle{\bf h}^{c,1}_{ijk+1/2}=$	$\displaystyle\overline{\rho w}_{ijk+1/2}$
$\displaystyle{\bf h}^{c,2}_{ijk+1/2}=$	$\displaystyle\bar{u}_{ijk+1/2}\overline{\rho w}_{ijk+1/2}$
$\displaystyle{\bf h}^{c,3}_{ijk+1/2}=$	$\displaystyle\bar{v}_{ijk+1/2}\overline{\rho w}_{ijk+1/2}$	(70)
$\displaystyle{\bf h}^{c,4}_{ijk+1/2}=$	$\displaystyle\bar{w}_{ijk+1/2}\overline{\rho w}_{ijk+1/2}+p_{ijk+1/2}$
$\displaystyle{\bf h}^{c,5}_{ijk+1/2}=$	$\displaystyle\frac{1}{2(\gamma-1)\hat{\beta}_{ijk+1/2}}\overline{\rho w}_{ijk+1/2}-\frac{\overline{\|{\bf v}\|^{2}}_{ijk+1/2}}{2}\overline{\rho w}_{ijk+1/2}$
	$\displaystyle+\overline{\|{\bf v}\|}^{2}_{ijk+1/2}\overline{\rho w}_{ijk+1/2}+p_{ijk+1/2}\bar{w}_{ijk+1/2}.$

Remark 10.

Above, we have used the convention that a numeral superscript signifies a component of a flux.

The diffusion coefficient is approximated as,

\displaystyle\nu_{i+1/2jk}

\displaystyle=\frac{\mu_{0}}{\hat{\rho}_{i+1/2jk}}+\mu_{1}\bar{\rho}_{i+1/2jk},

(71)

and similarly in the y- and z-directions.

Furthermore, we augment the diffusion with an with artificial (vanishing) component:

	$\displaystyle\lambda_{i+1/2jk}$	$\displaystyle=\|\bar{u}_{i+1/2jk}\|R^{*}_{i+1/2jk}+\frac{\|\Delta^{x}_{+}u_{ijk}\|}{4},$		(72)
	$\displaystyle R^{*}_{i+1/2jk}$	$\displaystyle=\max(\frac{1}{2},\|\Delta^{x}_{+}\log\rho_{ijk}\|).$

We need the following Lemma:

Lemma 1.

	$\displaystyle R^{*}_{i+1/2jk}\geq$
	$\displaystyle\max\left(\frac{1}{2},\left\|\frac{\Delta^{x}_{+}\rho_{i}}{12\bar{\rho}_{i+1/2}}\right\|,\left\|\frac{1}{2}\frac{(\sqrt{\rho_{i+1}}-\sqrt{\rho_{i}})}{\sqrt{\rho_{i+1}}+\sqrt{\rho_{i}}}\right\|,\frac{\bar{\rho}_{i+1/2}}{2}\left(\frac{\|\rho_{i+1}-\rho_{i}\|}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}\right),\frac{\|\rho_{i+1}-\rho_{i}\|}{\hat{\rho}_{i+1/2}}\right)_{jk}.$

(As usual the obvious generalisations to the other two dimensions holds.)

Proof.

We suppress the common $jk$ indices and check one term at the time. We use that $\rho_{i}>0$ (for all $i$ ) and the mean inequalities (62).

1. Trivial.

\displaystyle\frac{|\Delta_{+}^{x}\rho_{i}|}{\bar{\rho}_{i+1/2}}\leq\frac{|\Delta_{+}^{x}\rho_{i}|}{\hat{\rho}_{i+1/2}}=|\Delta_{+}^{x}\log\rho_{i}|.

	$\displaystyle\frac{(\sqrt{\rho_{i+1}}-\sqrt{\rho_{i}})}{\sqrt{\rho_{i+1}}+\sqrt{\rho_{i}}}=\frac{(\rho_{i+1}-\rho_{i})}{(\sqrt{\rho_{i+1}}+\sqrt{\rho_{i}})^{2}}=\frac{(\rho_{i+1}-\rho_{i})}{\rho_{i+1}+\rho_{i}+2\sqrt{\rho_{i+1}}\sqrt{\rho_{i}}}$	$\displaystyle\leq$
	$\displaystyle\frac{(\rho_{i+1}-\rho_{i})}{\rho_{i+1}+\rho_{i}}$	$\displaystyle\leq\frac{1}{2}\|\Delta_{+}^{x}\log\rho_{i}\|.$

	$\displaystyle\frac{\bar{\rho}_{i+1/2}}{2}\left(\frac{\|\rho_{i+1}-\rho_{i}\|}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}\right)\leq$
	$\displaystyle\frac{\bar{\rho}_{i+1/2}}{2}\left(\frac{\|\rho_{i+1}-\rho_{i}\|}{\frac{1}{2}(\rho_{i+1}+\rho_{i})^{2}}\right)=\frac{\bar{\rho}_{i+1/2}}{2}\left(\frac{\|\rho_{i+1}-\rho_{i}\|}{2\bar{\rho}_{i+1/2}^{2}}\right)$	$\displaystyle\leq\frac{1}{4}\|\Delta_{+}^{x}\log\rho_{i}\|.$

\displaystyle\frac{|\rho_{i+1}-\rho_{i}|}{\hat{\rho}_{i+1/2}}

\displaystyle=|\Delta_{+}^{x}\log\rho_{i}|.

∎

Next, we define

$\displaystyle\tilde{\nu}_{i+1/2jk}$	$\displaystyle=\nu_{i+1/2jk}+(\Delta_{+}^{x}x_{i})\lambda_{i+1/2jk},$
$\displaystyle\tilde{\nu}_{ij+1/2k}$	$\displaystyle=\nu_{ij+1/2k}+(\Delta_{+}^{y}y_{j})\lambda_{ij+1/2k},$	(73)
$\displaystyle\tilde{\nu}_{ijk+1/2}$	$\displaystyle=\nu_{ijk+1/2}+(\Delta_{+}^{z}z_{k})\lambda_{ijk+1/2}.$

The diffusive fluxes for $i,j,k=0...N-1$ are given by,

$\displaystyle{\bf f}^{d,1}_{i+1/2jk}$	$\displaystyle=\tilde{\nu}_{i+1/2jk}D^{x}_{+}\rho_{ijk},$
$\displaystyle{\bf g}^{d,1}_{ij+1/2k}$	$\displaystyle=\tilde{\nu}_{ij+1/2k}D^{y}_{+}\rho_{ijk},$	(74)
$\displaystyle{\bf h}^{d,1}_{ijk+1/2}$	$\displaystyle=\tilde{\nu}_{ijk+1/2}D^{z}_{+}\rho_{ijk},$

and

$\displaystyle{\bf f}^{d,\{2,3,4\}}_{i+1/2jk}$	$\displaystyle=\tilde{\nu}_{i+1/2jk}D^{x}_{+}(\rho_{ijk}\{u_{ijk},v_{ijk},w_{ijk}\}),$
$\displaystyle{\bf g}^{d,\{2,3,4\}}_{ij+1/2k}$	$\displaystyle=\tilde{\nu}_{ij+1/2k}D^{y}_{+}(\rho_{ijk}\{u_{ijk},v_{ijk},w_{ijk}\}),$	(75)
$\displaystyle{\bf h}^{d,\{2,3,4\}}_{ijk+1/2}$	$\displaystyle=\tilde{\nu}_{ijk+1/2}D^{z}_{+}(\rho_{ijk}\{u_{ijk},v_{ijk},w_{ijk}\}),$

and finally,

$\displaystyle{\bf f}^{d,5}_{i+1/2jk}=$	$\displaystyle\tilde{\nu}_{i+1/2jk}\left(\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1}+\frac{1}{2}D^{x}_{+}(\rho_{ijk}\|{\bf v}\|^{2}_{ijk})+(\overline{\|{\bf v}\|}^{2}_{i+1/2jk}-\overline{\|{\bf v}\|^{2}}_{i+1/2jk})D^{x}_{+}\rho_{ijk}\right)$
	$\displaystyle+\kappa_{r}D^{x}_{+}T^{4}_{ijk}$
$\displaystyle{\bf g}^{d,5}_{ij+1/2k}=$	$\displaystyle\tilde{\nu}_{ij+1/2k}\left(\frac{({\mathfrak{p}}_{y})_{ij+1/2k}}{\gamma-1}+\frac{1}{2}D^{y}_{+}(\rho_{ijk}\|{\bf v}\|^{2}_{ijk})+(\overline{\|{\bf v}\|}^{2}_{ij+1/2k}-\overline{\|{\bf v}\|^{2}}_{ij+1/2k})D^{y}_{+}\rho_{ijk}\right)$
	$\displaystyle+\kappa_{r}D^{y}_{+}T^{4}_{ijk}$	(76)
$\displaystyle{\bf h}^{d,5}_{ijk+1/2}=$	$\displaystyle\tilde{\nu}_{ijk+1/2}\left(\frac{({\mathfrak{p}}_{z})_{ijk+1/2}}{\gamma-1}+\frac{1}{2}D^{z}_{+}(\rho_{ijk}\|{\bf v}\|^{2}_{ijk})+(\overline{\|{\bf v}\|}_{ijk+1/2}^{2}-\overline{\|{\bf v}\|^{2}}_{ijk+1/2})D^{z}_{+}\rho_{ijk}\right)$
	$\displaystyle+\kappa_{r}D^{z}_{+}T^{4}_{ijk}$

where

	$\displaystyle({\mathfrak{p}}_{x})_{i+1/2jk}$	$\displaystyle=\frac{1}{2\hat{\beta}_{i+1/2jk}}D_{+}^{x}\rho_{ijk}+\frac{\bar{\rho}_{i+1/2jk}}{2}D_{+}^{x}\frac{1}{\beta_{ijk}},$
	$\displaystyle({\mathfrak{p}}_{y})_{ij+1/2k}$	$\displaystyle=\frac{1}{2\hat{\beta}_{ij+1/2k}}D_{+}^{y}\rho_{ijk}+\frac{\bar{\rho}_{ij+1/2k}}{2}D_{+}^{y}\frac{1}{\beta_{ijk}},$
	$\displaystyle({\mathfrak{p}}_{z})_{ijk+1/2}$	$\displaystyle=\frac{1}{2\hat{\beta}_{ijk+1/2}}D_{+}^{z}\rho_{ijk}+\frac{\bar{\rho}_{ijk+1/2}}{2}D_{+}^{z}\frac{1}{\beta_{ijk}}.$

Furthermore, we isolate the artificial diffusion by splitting the diffusive fluxes into

	$\displaystyle{\bf f}^{d}_{i+1/2jk}$	$\displaystyle={\bf f}^{\nu}_{i+1/2jk}+{\bf f}^{\lambda}_{i+1/2jk},$
	$\displaystyle{\bf g}^{d}_{ij+1/2k}$	$\displaystyle={\bf g}^{\nu}_{ij+1/2k}+{\bf g}^{\lambda}_{ij+1/2k},$
	$\displaystyle{\bf h}^{d}_{ijk+1/2}$	$\displaystyle={\bf h}^{\nu}_{ijk+1/2}+{\bf h}^{\lambda}_{ijk+1/2},$

where the $\lambda$ -fluxes are obtained by replacing $\tilde{\nu}$ with $\lambda$ , and $\kappa_{r}=0$ , in (74),(75) and (76). (The $\nu$ -fluxes are thus given by the ${\bf f}^{d}$ with $\nu$ replacing $\tilde{\nu}$ and include the $\kappa_{r}$ part.)

Boundary conditions: We invoke the no-slip condition strongly. That is, $u,v,w$ are all zero at the boundaries. Therefore,

\displaystyle{\bf u}_{0jk}=(\rho_{0jk},0,0,0,\frac{p_{0jk}}{\gamma-1})^{T},\quad{\bf u}_{Njk}=(\rho_{Njk},0,0,0,\frac{p_{Njk}}{\gamma-1})^{T},

(77)

which are used in the computation of ${\bf f}_{1/2jk},{\bf f}_{N-1/2jk}$ . (Similarly in the other two directions.)

Since the velocities are zero at the boundary, there should be no scheme updating momentum at the boundary points. To to keep the simple structure of (66), we enforce this by stipulating that every partial flux, ${\bf f},{\bf f}^{c},{\bf f}^{d},{\bf f}^{\lambda},...$ satisfy

\displaystyle{\bf f}^{2,3,4}_{-1/2jk}={\bf f}^{2,3,4}_{1/2jk},\quad{\bf g}^{2,3,4}_{0j+1/2k}=0,\quad{\bf h}^{2,3,4}_{0jk+1/2}=0,

(78)

and similarly for all other boundaries. This choice amounts to $(\rho{\bf v})_{t}=0$ along the boundaries. Hence, if the momentum/velocity is initially zero, it will remain zero for all time.

Furthermore, invoking the no-slip in the remaining convective boundary flux components (for continuity and total energy), implies that they are all identically zero,

	$\displaystyle{\bf f}^{c,1}_{-1/2jk}={\bf g}^{c,1}_{-1/2jk}={\bf h}^{c,1}_{-1/2jk}=0,$			(79)
	$\displaystyle{\bf f}^{c,5}_{-1/2jk}={\bf g}^{c,5}_{-1/2jk}={\bf h}^{c,5}_{-1/2jk}=0,$

and similarly at the other boundaries.

At the boundaries, the diffusive momentum fluxes are given by (78) which ensures that the momentum terms remain zero (in accordance with the no-slip condition). Using no-slip and homogeneous Neumann approximations for the wall normal gradient of density and temperature, the remaining two diffusive flux components are:

$\displaystyle{\bf f}^{d,1}_{-1/2jk}={\bf f}^{d,1}_{N+1/2jk}=0,$	$\displaystyle\quad{\bf f}^{d,5}_{-1/2jk}={\bf f}^{d,5}_{N+1/2jk}=0,$
$\displaystyle{\bf g}^{d,1}_{i,-1/2k}={\bf g}^{d,1}_{i,N+1/2k}=0,$	$\displaystyle\quad{\bf g}^{d,5}_{i,-1/2k}={\bf g}^{d,5}_{i,N+1/2k}=0,$	(80)
$\displaystyle{\bf h}^{d,1}_{ij,-1/2}={\bf h}^{d,1}_{ij,N+1/2}=0,$	$\displaystyle\quad{\bf h}^{d,5}_{ij,-1/2}={\bf h}^{d,5}_{ij,N+1/2}=0$

Note that it is only the normal boundary flux that is affected by the Neumann boundary condition. The diffusive boundary fluxes (for the 1st and 5th components) tangent to the boundaries are not subject to any boundary conditions.

Initial conditions: We make the following assumptions on the initial data (t=0):

•

All variables are uniformly bounded.
•

Density and temperature are positive and bounded away from zero.
•

At the boundary points, we set ${\bf v}=0$ at $t=0$ (no-slip).
•

The initial data is projected onto the grid by e.g. computing the averages in each control volume. This also determines $\rho$ and $T$ at the boundary points. (We remark that injecting, instead of averaging, initial values to the grid points, does not change the overall formal accuracy of the method.)

It is possible to relax these assumptions but for practical purposes they are sufficient. Note that there is no smallness assumption on the initial data.

The following is the main result of this paper:

Theorem 1.

Assume that at $t=0$ the initial data $\rho^{0},{\bf v}^{0},T^{0}$ are bounded in $L^{\infty}(\Omega)$ and $\rho^{0},T^{0}\geq constant>0$ . Let $\{\rho_{h},{\bf v}_{h},T_{h}\}$ be the family of numerical solutions generated from the initial data by the scheme (66) (see Section 4.1 for all approximations) as $h\rightarrow 0$ . Then a subsequence converges strongly in the spaces,

	$\displaystyle\rho_{h}\rightarrow\rho$	$\displaystyle\quad\in L^{23-\epsilon}_{23-\epsilon}$
	$\displaystyle{\bf v}_{h}\rightarrow{\bf v}$	$\displaystyle\quad\in L^{10/4-\epsilon}_{10/4-\epsilon},\quad\quad(\epsilon>0),$
	$\displaystyle T_{h}\rightarrow T$	$\displaystyle\quad\in L^{4-\epsilon}_{4-\epsilon},$

and $\{\rho,{\bf v},T\}$ is a weak entropy solution of the problem (4) in the sense of Def. 1.

We summarise some of the estimates that resulting solution satisfies:

	$\displaystyle\rho\in L^{\infty}(L^{4}(\Omega))\cap L^{23/3}_{23/3}\cap L^{5}_{15},$	$\displaystyle\quad\rho^{5/2}\in L^{2}(H^{1}),$
	$\displaystyle\rho^{-1}\in L^{\infty}_{1}\cap L^{2}(H^{1})\cap L^{2}(L^{6}),$	$\displaystyle\quad\nabla(\log(\rho))\in L^{2}_{2},$
	$\displaystyle{\bf v}\in L^{\infty}_{1}\cap L^{2}(H^{1})\cap L^{2}(L^{6})\cap L^{10/4}_{3},$	$\displaystyle\quad\rho\nabla{\bf v}\in L^{2}_{2},$
	$\displaystyle T\in L^{4+\delta}_{4+\delta},\quad\textrm{for some}\quad\delta>0,$
	$\displaystyle T^{2-\epsilon}\in L^{2}(H^{1})\cap L^{2}(L^{6}),\quad\textrm{for}\quad\epsilon>0,$
	$\displaystyle\sqrt{T}\in L^{\infty}_{1},\quad p\in L^{\infty}_{2},$	$\displaystyle\quad\rho\log(T)\in L^{\infty}_{1}.$

The remainder of the paper contains a proof of the theorem.

5 Additional balance laws

Before we derive the discrete versions of the kinetic and internal energy balances, we state some auxiliary relations. For positive entities the following relations (exemplified with the positive density) hold:

\displaystyle\frac{1}{\check{\rho}_{i+1/2jk}}\geq\frac{1}{\hat{\rho}_{i+1/2jk}}\geq\frac{1}{\bar{\rho}_{i+1/2jk}}.

(81)

They follow directly from (62).

We shall also frequently use (63) to obtain the following alternative form of the first component of the internal inviscid fluxes and artificial diffusion (here exemplified with the x-flux).

	$\displaystyle{\bf f}^{c,1}_{i+1/2jk}-{\bf f}^{\lambda,1}_{i+1/2jk}=$	$\displaystyle\overline{\rho u}_{i+1/2jk}-\lambda_{i+1/2jk}\Delta_{+}^{x}\rho_{ijk}$
	$\displaystyle=$	$\displaystyle\bar{\rho}_{i+1/2jk}\bar{u}_{i+1/2jk}-\lambda^{a}_{i+1/2jk}\Delta_{+}^{x}\rho_{ijk},$		(82)

where

\displaystyle\lambda^{a}_{i+1/2jk}

\displaystyle=\lambda_{i+1/2jk}-\frac{\Delta^{x}_{+}u_{ijk}}{4}>0,

and

	$\displaystyle{\bf f}^{c,1}_{i+1/2jk}-{\bf f}^{\lambda,1}_{i+1/2jk}=$	$\displaystyle\overline{\rho}_{i+1/2jk}\bar{u}_{i+1/2jk}-\lambda^{a}_{i+1/2jk}\Delta_{+}^{x}\rho_{ijk}$
	$\displaystyle=$	$\displaystyle\check{\rho}_{i+1/2jk}\bar{u}_{i+1/2jk}-\lambda^{c}_{i+1/2jk}\Delta_{+}^{x}\rho_{ijk},$		(83)

where

\displaystyle\lambda^{c}_{i+1/2jk}

\displaystyle=\lambda_{i+1/2jk}-\frac{\Delta^{x}_{+}u_{ijk}}{4}+\bar{u}_{i+1/2jk}\frac{\bar{\rho}_{i+1/2}-\check{\rho}_{i+1/2jk}}{\Delta_{+}^{x}\rho_{ijk}}>0.

The last inequality follows from

	$\displaystyle\frac{\bar{\rho}_{i+1/2}-\check{\rho}_{i+1/2jk}}{\Delta_{+}^{x}\rho_{ijk}}=\frac{\frac{\rho_{i+1}+\rho_{i}}{2}-\sqrt{\rho_{i+1}\rho_{i}}}{\Delta_{+}^{x}\rho_{ijk}}$	$\displaystyle=$
	$\displaystyle\frac{1}{2}\frac{(\sqrt{\rho_{i+1}}-\sqrt{\rho_{i}})^{2}}{(\sqrt{\rho_{i+1}-\sqrt{\rho_{i}}})(\sqrt{\rho_{i+1}+\sqrt{\rho_{i}}})}=\frac{1}{2}\frac{(\sqrt{\rho_{i+1}}-\sqrt{\rho_{i}})}{\sqrt{\rho_{i+1}+\sqrt{\rho_{i}}}}$			(84)

and Lemma 1.

5.1 Kinetic and internal energy balances

We derive the kinetic energy balance that will subsequently be subtracted from the total energy equation to produce the internal energy balance.

\displaystyle\frac{d}{dt}K_{ijk}=-\frac{1}{2}|{\bf v}|^{2}_{ijk}\frac{d\rho_{j}}{dt}+u_{ijk}\frac{d(\rho u)_{ijk}}{dt}+v_{ijk}\frac{d(\rho v)_{ijk}}{dt}+w_{ijk}\frac{d(\rho w)_{ijk}}{dt}.

(85)

We need the Leibniz rule

\displaystyle a_{i}D_{-}f_{i+1/2}=D_{-}(\bar{a}_{i+1/2}f_{i+1/2})-\frac{1}{2}f_{i+1/2}D_{-}a_{i+1}-\frac{1}{2}f_{i-1/2}D_{-}a_{i}.

(86)

At the boundaries, we have the corresponding rules

\displaystyle a_{0}D_{-}f_{1/2}=D_{-}(\bar{a}_{1/2}f_{1/2})-\frac{1}{2}f_{1/2}D_{-}a_{1}

where, in accordance with our convention, $\bar{a}_{-1/2}=a_{0}$ . By adopting the convention that $D_{-}a_{0}\equiv 0$ , (86) applies to the boundaries as well, and the derivations below need not be repeated for the boundary points.

First we rewrite,

$\displaystyle u_{ijk}\frac{d(\rho u)_{ijk}}{dt}$	$\displaystyle=$
$\displaystyle-u_{ijk}\left(D_{-}^{x}{\bf f}^{2}_{i+1/2jk}+D_{-}^{y}{\bf g}^{2}_{ij+1/2k}+D_{-}^{z}{\bf h}^{2}_{ijk+1/2}\right)$	$\displaystyle=$
$\displaystyle-\left(D_{-}^{x}(\bar{u}_{i+1/2jk}{\bf f}^{2}_{i+1/2jk})-\frac{1}{2}(D_{-}^{x}u_{i+1jk}){\bf f}^{2}_{i+1/2jk}-\frac{1}{2}(D_{-}^{x}u_{ijk}){\bf f}^{2}_{i-1/2jk}\right)$
$\displaystyle-\left(D_{-}^{y}(\bar{u}_{ij+1/2k}{\bf g}^{2}_{ij+1/2k})-\frac{1}{2}(D_{-}^{y}u_{ij+1k}){\bf g}^{2}_{ij+1/2k}-\frac{1}{2}(D_{-}^{y}u_{ijk}){\bf g}^{2}_{ij-1/2k}\right)$		(87)
$\displaystyle-\left(D_{-}^{z}(\bar{u}_{ijk+1/2}{\bf h}^{2}_{ijk+1/2})-\frac{1}{2}(D_{-}^{z}u_{ijk+1}){\bf h}^{2}_{ijk+1/2}-\frac{1}{2}(D_{-}^{z}u_{ijk}){\bf h}^{2}_{ijk-1/2}\right),$

and similarly for the momentum in the y- and z-directions.

Turning to the contribution from the continuity equation, and again using (86), we have

$\displaystyle-\frac{1}{2}\|{\bf v}\|^{2}_{ijk}\frac{d\rho_{j}}{dt}$	$\displaystyle=$
$\displaystyle\frac{1}{2}\|{\bf v}\|^{2}_{ijk}\left(D_{-}^{x}{\bf f}^{1}_{i+1/2jk}+D_{-}^{y}{\bf g}^{1}_{ij+1/2k}+D_{-}^{z}{\bf h}^{1}_{ijk+1/2}\right)$	$\displaystyle=$
$\displaystyle\frac{1}{2}\left(D_{-}^{x}(\overline{\|{\bf v}\|^{2}}_{i+1/2jk}{\bf f}^{1}_{i+1/2jk})-\frac{1}{2}(D_{-}^{x}\|{\bf v}\|^{2}_{i+1jk}){\bf f}^{1}_{i+1/2jk}-\frac{1}{2}(D_{-}^{x}\|{\bf v}\|^{2}_{ijk}){\bf f}^{1}_{i-1/2jk}\right)$
$\displaystyle+\frac{1}{2}\left(D_{-}^{y}(\overline{\|{\bf v}\|^{2}}_{ij+1/2k}{\bf g}^{1}_{ij+1/2k})-\frac{1}{2}(D_{-}^{y}\|{\bf v}\|^{2}_{ij+1/2k}){\bf g}^{1}_{ij+1/2k}-\frac{1}{2}(D_{-}^{y}\|{\bf v}\|^{2}_{ijk}){\bf g}^{1}_{ij-1/2k}\right)$		(88)
$\displaystyle+\frac{1}{2}\left(D_{-}^{z}(\overline{\|{\bf v}\|^{2}}_{ijk+1/2}{\bf h}^{1}_{ijk+1/2})-\frac{1}{2}(D_{-}^{z}\|{\bf v}\|^{2}_{ijk+1/2}){\bf h}^{1}_{ijk+1/2}-\frac{1}{2}(D_{-}^{z}\|{\bf v}\|^{2}_{ijk}){\bf h}^{1}_{ijk-1/2}\right).$

Using the identity,

	$\displaystyle D_{-}^{x}\|{\bf v}\|^{2}_{ijk}=$	$\displaystyle D_{-}^{x}(u_{ijk}^{2}+v^{2}_{ijk}+w_{ijk}^{2})=$
		$\displaystyle 2\bar{u}_{i-1/2}D_{-}^{x}u_{ijk}+2\bar{v}_{i-1/2}D_{-}^{x}v_{ijk}+2\bar{w}_{i-1/2}D_{-}^{x}w_{ijk},$

and (87)-(88) in (85), we obtain an expression for the kinetic energy balance:

	$\displaystyle\frac{dK_{ijk}}{dt}+D_{-}^{x}{\mathfrak{F}}_{i+1/2jk}+D_{-}^{y}{\mathfrak{G}}_{ij+1/2k}+D_{-}^{z}{\mathfrak{H}}_{ijk+1/2}-(pD{\bf v})_{ijk}$	$\displaystyle=-{\mathfrak{D}}_{ijk},$		(89)
	$\displaystyle 0\leq i,j,k\leq N,$

where

$\displaystyle{\mathfrak{F}}_{i+1/2jk}$	$\displaystyle=-\frac{1}{2}\overline{\|{\bf v}\|^{2}}_{i+1/2jk}{\bf f}^{1}_{i+1/2jk}+\bar{u}_{i+1/2jk}{\bf f}^{2}_{i+1/2jk}+\bar{v}_{i+1/2jk}{\bf f}^{3}_{i+1/2jk}+\bar{w}_{i+1/2jk}{\bf f}^{4}_{i+1/2jk}$
$\displaystyle{\mathfrak{G}}_{ij+1/2k}$	$\displaystyle=-\frac{1}{2}\overline{\|{\bf v}\|^{2}}_{ij+1/2k}{\bf g}^{1}_{ij+1/2k}+\bar{u}_{ij+1/2k}{\bf g}^{2}_{ij+1/2k}+\bar{v}_{ij+1/2k}{\bf g}^{3}_{ij+1/2k}+\bar{w}_{ij+1/2k}{\bf g}^{4}_{ij+1/2k}$	(90)
$\displaystyle{\mathfrak{H}}_{ijk+1/2}$	$\displaystyle=-\frac{1}{2}\overline{\|{\bf v}\|^{2}}_{ijk+1/2}{\bf h}^{1}_{ijk+1/2}+\bar{u}_{ijk+1/2}{\bf h}^{2}_{ijk+1/2}+\bar{v}_{ijk+1/2}{\bf h}^{3}_{ijk+1/2}+\bar{w}_{ijk+1/2}{\bf h}^{4}_{ijk+1/2},$

and

	$\displaystyle(pD{\bf v})_{ijk}=$	$\displaystyle\quad\frac{1}{2}\left((D_{-}^{x}u_{i+1jk})p_{i+1/2jk}+(D_{-}^{x}u_{ijk})p_{i-1/2jk}\right)$
		$\displaystyle+\frac{1}{2}\left((D_{-}^{y}v_{ij+1k})p_{ij+1/2k}+(D_{-}^{y}v_{ijk})p_{ij-1/2k}\right)$
		$\displaystyle+\frac{1}{2}\left((D_{-}^{z}w_{ijk+1})p_{ijk+1/2}+(D_{-}^{z}w_{ijk})p_{ijk-1/2}\right).$

Using the boundary convention introduced in (86) and the no-slip condition (implying that velocity differences along boundaries are zero), we have, ${\mathfrak{F}}_{-1/2jk}\equiv 0$ , etc. and

\displaystyle(pD{\bf v})_{0jk}=

\displaystyle\quad\frac{1}{2}\left((D_{-}^{x}u_{1jk})p_{i+1/2jk}\right)

and similarly for the other boundaries. The diffusive terms are,

	$\displaystyle{\mathfrak{D}}_{ijk}$	$\displaystyle=\frac{1}{2}\left(\tilde{\nu}_{i+1/2jk}\bar{\rho}_{i+1/2jk}(\|D_{-}^{x}{\bf v}_{i+1jk}\|)^{2}+\tilde{\nu}_{i-1/2jk}\bar{\rho}_{i-1/2jk}(\|D_{-}^{x}{\bf v}_{ijk}\|)^{2}\right)$
		$\displaystyle+\frac{1}{2}\left(\tilde{\nu}_{ij+1/2k}\bar{\rho}_{ij+1/2k}(\|D_{-}^{y}{\bf v}_{ij+1k}\|)^{2}+\tilde{\nu}_{ij-1/2k}\bar{\rho}_{ij-1/2k}(\|D_{-}^{y}{\bf v}_{ijk}\|)^{2}\right)$
		$\displaystyle+\frac{1}{2}\left(\tilde{\nu}_{ijk+1/2}\bar{\rho}_{ijk+1/2}(\|D_{-}^{z}{\bf v}_{ijk+1}\|)^{2}+\tilde{\nu}_{ijk-1/2}\bar{\rho}_{ijk-1/2}(\|D_{-}^{z}{\bf v}_{ijk\|})^{2}\right),$

where,

\displaystyle(|D_{-}^{x}{\bf v}_{ijk}|)^{2}=(|D_{-}^{x}u_{ijk}|)^{2}+(|D_{-}^{x}w_{ijk}|)^{2}+(|D_{-}^{x}w_{ijk}|)^{2},etc.

At the boundary, most of the terms of ${\mathfrak{D}}_{ijk}$ are zero due to the no-slip condition and the boundary convention in (86). For instance,

\displaystyle{\mathfrak{D}}_{0jk}=\frac{1}{2}\tilde{\nu}_{1/2jk}\bar{\rho}_{1/2jk}(|D_{-}^{x}{\bf v}_{1jk}|)^{2}.

Subtracting the kinetic energy equation (89) from the total energy equation of (67) gives the internal energy equation:

$\displaystyle D_{t}\frac{p_{ijk}}{\gamma-1}+(pD{\bf v})_{ijk}-{\mathfrak{D}}_{ijk}$
$\displaystyle+D_{-}^{x}(\frac{\overline{\rho u}_{i+1/2jk}}{2(\gamma-1)\hat{\beta}_{i+1/2jk}})+D_{-}^{y}(\frac{\overline{\rho v}_{ij+1/2k}}{2(\gamma-1)\hat{\beta}_{ij+1/2k}})+D_{-}^{y}(\frac{\overline{\rho w}_{ijk+1/2}}{2(\gamma-1)\hat{\beta}_{ijk+1/2}})$	$\displaystyle=$	(91)
$\displaystyle D_{-}^{x}(\tilde{\nu}_{i+1/2jk}\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1})+D_{-}^{y}(\tilde{\nu}_{ij+1/2k}\frac{({\mathfrak{p}}_{y})_{ij+1/2k}}{\gamma-1})+$
$\displaystyle D_{-}^{z}(\tilde{\nu}_{ijk+1/2}\frac{({\mathfrak{p}}_{z})_{ijk+1/2}}{\gamma-1})$
$\displaystyle+\kappa_{r}\left(D_{-}^{x}D^{x}_{+}T^{4}_{ijk}+D_{-}^{y}D^{y}_{+}T^{4}_{ijk}+D_{-}^{z}\kappa_{r}D^{z}_{+}T^{4}_{ijk}\right)$	.

6 Stability

Our aim is to derive the same a priori estimates for the semi-discrete system as we found for the PDE system. To carry out the calculations, we need the following summation-by-parts formula,

\displaystyle\sum_{i=0}^{N}a_{i}\Delta_{-}b_{i+1/2}

\displaystyle=-a_{0}b_{-1/2}+a_{N}b_{N+1/2}-\sum_{i=0}^{N-1}(\Delta_{+}a_{i})b_{i+1/2},

(92)

and product rule

\displaystyle\Delta_{+}a_{i}b_{i}

\displaystyle=\bar{a}_{i+1/2}\Delta_{+}b_{i}+\bar{b}_{i+1/2}\Delta_{+}a_{i}.

(93)

6.1 Conservation

As in the continuous case, we obtain

\displaystyle\rho_{ijk},\,E_{ijk},p_{i}\in L^{\infty}_{1},\quad\sqrt{\rho_{ijk}}|{\bf v}_{ijk}|\in L^{\infty}_{2},

(94)

from conservation, positivity and the boundary conditions.

Remark 11.

In the notation above, $\rho_{ijk}$ is not a point value but symbolises the piecewise constant (in space) and continuous in time grid function on the entire domain that is bounded in $L^{\infty}(0,\mathcal{T};L^{1}(\Omega))$

6.2 Entropy

Here, we use the same scaling of the entropy as in [Cha13]. Dropping the $ijk$ indices, the entropy variables are:

\displaystyle{\bf w}=\left(\frac{\gamma}{\gamma-1}-\frac{s}{\gamma-1}-\beta|{\bf v}|^{2},2\beta u,2\beta v,2\beta w,-2\beta\right),

where $|{\bf v}|^{2}=u^{2}+v^{2}+w^{2}$ and the specific entropy,

\displaystyle s=\log\frac{p}{\rho^{\gamma}}=-\log\beta-(\gamma-1)\log\rho-\log 2.

Furthermore, the entropy potentials are: $\Psi^{x,y,z}=m_{1,2,3}=\rho\{u,v,w\}$ .

An entropy estimate is obtained by contracting the scheme with the entropy variables, summing in space and integrating in time. We begin by multiplying the diffusive fluxes (without artificial viscosity) by ${\bf w}_{ijk}^{T}$ and summing in space.

That is,

$\displaystyle\sum_{ijk=0}^{N}{\bf w}_{ijk}^{T}\left((\Delta_{-}^{x}{\bf f}^{\nu}_{i+1/2jk}){\mathcal{S}}^{x}_{jk}+(\Delta^{y}_{-}{\bf g}_{ij+1/2k}^{\nu}){\mathcal{S}}^{y}_{ik}+(D_{-}^{z}{\bf h}^{\nu}_{ijk+1/2}){\mathcal{S}}^{z}_{ij}\right)$	$\displaystyle=$
$\displaystyle\sum_{jk=0}^{N}({\bf w}_{Njk}^{T}{\bf f}^{\nu}_{N+1/2jk}-{\bf w}_{0jk}^{T}{\bf f}^{\nu}_{-1/2jk}){\mathcal{S}}^{x}_{jk}+\sum_{ik=0}^{N}({\bf w}_{iNk}^{T}{\bf g}_{iN+1/2k}^{\nu}-{\bf w}_{i0k}^{T}{\bf g}^{\nu}_{i,-1/2k}){\mathcal{S}}^{y}_{ik}$
$\displaystyle+\sum_{jk=0}^{N}({\bf w}_{ijN}^{T}{\bf h}^{\nu}_{ijN+1/2}-{\bf w}_{ij0}^{T}{\bf h}^{\nu}_{ij,-1/2}){\mathcal{S}}^{z}_{ij}$
$\displaystyle-\sum_{ijk=0}^{N-1,N,N}(\Delta_{+}^{x}{\bf w}_{ijk})^{T}{\bf f}^{\nu}_{i+1/2jk}{\mathcal{S}}^{x}_{jk}-\sum_{ijk=0}^{N,N-1,N}(\Delta^{y}_{+}{\bf w}_{ijk})^{T}{\bf g}_{ij+1/2k}^{\nu}{\mathcal{S}}^{y}_{ik}$		(95)
$\displaystyle-\sum_{ijk=0}^{N,N,N-1}(\Delta^{z}_{+}{\bf w}_{ijk})^{T}{\bf h}^{\nu}_{ijk+1/2}{\mathcal{S}}^{z}_{ij}$	.

The boundary conditions (77) and numerical boundary fluxes (80) make all boundary terms vanish.

To further, manipulate (95), we need the entropy-variable differences,

$\displaystyle\Delta_{+}{\bf w}^{1}_{i}=$	$\displaystyle\frac{\Delta_{+}\rho_{i}}{\hat{\rho}_{i+1/2}}+\left(\frac{1}{(\gamma-1)\hat{\beta}_{i+1/2}}-\overline{\|{\bf u}\|^{2}}_{i+1/2}\right)\Delta_{+}\beta_{i},$
	$\displaystyle-2\bar{u}_{i+1/2}\bar{\beta}_{i+1/2}\Delta_{+}u_{i}-2\bar{v}_{i+1/2}\bar{\beta}_{i+1/2}\Delta_{+}v_{i}-2\bar{w}_{i+1/2}\bar{\beta}_{i+1/2}\Delta_{+}w_{i},$
$\displaystyle\Delta_{+}{\bf w}^{2}_{i}=$	$\displaystyle 2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}+2\bar{u}_{i+1/2}\Delta_{+}\beta_{i},$	(96)
$\displaystyle\Delta_{+}{\bf w}^{3}_{i}=$	$\displaystyle 2\bar{\beta}_{i+1/2}\Delta_{+}v_{i}+2\bar{v}_{i+1/2}\Delta_{+}\beta_{i},$
$\displaystyle\Delta_{+}{\bf w}^{4}_{i}=$	$\displaystyle 2\bar{\beta}_{i+1/2}\Delta_{+}w_{i}+2\bar{w}_{i+1/2}\Delta_{+}\beta_{i},$
$\displaystyle\Delta_{+}{\bf w}^{5}_{i}=$	$\displaystyle-2\Delta_{+}\beta_{i}.$

The relations are stated generically to any difference direction. That is, $i$ represents one of the three indices and the other two are held constant.

Remark 12.

A derivation of the identities (96) is found in [Cha13]. We have verified their correctness and made the trivial extension to 3-D in the current treatise.

Furthermore, we use (93) to obtain the two auxiliary relations,

	$\displaystyle D_{+}^{x}(\rho_{ijk}u_{ijk})$	$\displaystyle=\bar{\rho}_{i+1/2}D_{+}^{x}(u_{ijk})+\bar{u}_{i+1/2jk}D_{+}^{x}(\rho_{ijk})$
	$\displaystyle\frac{1}{2}D_{+}^{x}(\rho_{i}u^{2}_{i})$	$\displaystyle=\frac{1}{2}(D_{+}^{x}\rho_{i})\overline{u^{2}}_{i+1/2}+\bar{\rho}_{i+1/2}\bar{u}_{i+1/2}D_{+}^{x}u_{i},$

which along with (96) is used to calculate the remaining (volume) diffusive terms in (95). There are a number of cancellations and the three directions are handled independently of each other. Suppressing the common $jk$ -indices, the terms in the x-direction from (95) become,

	$\displaystyle(\Delta_{+}^{x}{\bf w}_{i}){\bf f}^{\nu}_{i+1/2}=\nu_{i+1/2jk}\left(\frac{(\Delta_{+}^{x}\rho_{i})(D_{+}^{x}\rho_{i})}{\hat{\rho}_{i+1/2}}\right.$
	$\displaystyle+2\bar{\beta}_{i+1/2}\bar{\rho}_{i+1/2}\left((\Delta^{x}_{+}u_{i})(D^{x}_{+}u_{i})+(\Delta^{x}_{+}v_{i})(D^{x}_{+}v_{i})+(\Delta^{x}_{+}w_{i})(D^{x}_{+}w_{i})\right)$
	$\displaystyle\left.-\frac{\bar{\rho}_{i+1/2}}{\gamma-1}\Delta_{+}^{x}\beta_{i}D_{+}^{x}\left(\frac{1}{\beta_{ijk}}\right)\right)-\frac{\kappa_{r}}{4}2(\Delta_{+}^{x}\beta_{ijk})D_{+}^{x}T^{4}_{i}.$		(97)

Furthermore, by (64) and (65), we have

\displaystyle\sum_{ijk=0}^{N-1,N,N}(\Delta_{+}^{x}{\bf w}_{ijk}){\bf f}^{\nu}_{i+1/2jk}{\mathcal{S}}_{jk}=\sum_{ijk=0}^{N-1,N,N}(D_{+}^{x}{\bf w}_{ijk}){\bf f}^{\nu}_{i+1/2jk}V_{ijk},

and similarly in the other directions.

Turning to the inviscid fluxes, we must verify that the convective numerical fluxes, augmented with artificial diffusion, are entropy dissipative in the following sense:

$\displaystyle\Delta_{+}^{x}{\bf w}^{T}_{ijk}({\bf f}^{c}_{i+1/2jk}-{\bf f}^{\lambda}_{i+1/2jk})$	$\displaystyle\leq\Delta\Psi^{x}_{ijk}=\Delta^{x}_{+}(\rho u)_{ijk},$
$\displaystyle\Delta_{+}^{y}{\bf w}^{T}_{ijk}({\bf g}^{c}_{ij+1/2k}-{\bf g}^{\lambda}_{ij+1/2k})$	$\displaystyle\leq\Delta\Psi^{y}_{ijk}=\Delta^{y}_{+}(\rho v)_{ijk},$	(98)
$\displaystyle\Delta_{+}^{z}{\bf w}^{T}_{ijk}({\bf h}^{c}_{ijk+1/2}-{\bf h}^{d}_{ijk+1/2})$	$\displaystyle\leq\Delta\Psi^{z}_{ijk}=\Delta^{z}_{+}(\rho w)_{ijk}.$

(The relations (98) are sometimes called Tadmor’s shuffle conditions.)

To reduce notation, we demonstrate that (98) holds in the one-dimensional case. (It is straightforward to extend the analysis to the three-dimensional case.) To carry out the derivation, we temporarily relabel ${\bf w}^{5}$ as ${\bf w}^{3}$ and similarly for the fluxes. In view of (98), we begin by calculating

\displaystyle LHS^{c}=(\Delta_{+}{\bf w}^{1}_{i}){\bf f}^{1,c}_{i+1/2}+(\Delta_{+}{\bf w}^{2}_{i}){\bf f}^{2,c}_{i+1/2}+(\Delta_{+}{\bf w}^{3}_{i}){\bf f}^{3,c}_{i+1/2}

Using (96) and (68) (reduced to 1-D), we have

	$\displaystyle LHS^{c}=$
	$\displaystyle(\frac{\Delta_{+}\rho_{i}}{\hat{\rho}_{i+1/2}}+\left(\frac{1}{(\gamma-1)\hat{\beta}_{i+1/2}}-\overline{\|u\|^{2}}_{i+1/2}\right)\Delta_{+}\beta_{i}-2\bar{u}_{i+1/2}\bar{\beta}_{i+1/2}\Delta_{+}u_{i})(\overline{\rho u}_{i+1/2})$
	$\displaystyle+(2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}+2\bar{u}_{i+1/2}\Delta_{+}\beta_{i})(\bar{u}_{i+1/2}\overline{\rho u}_{i+1/2}+p_{i+1/2})$
	$\displaystyle-2\Delta_{+}\beta_{i}(\frac{1}{2(\gamma-1)\hat{\beta}_{i+1/2}}\overline{\rho u}_{i+1/2}-\frac{\overline{\|u\|^{2}}_{i+1/2}}{2}\overline{\rho u}_{i+1/2}$
	$\displaystyle+\overline{\|u\|}^{2}_{i+1/2}\overline{\rho u}_{i+1/2}+p_{i+1/2}\bar{u}_{i+1/2})$	.

After some further manipulations, we obtain

	$\displaystyle LHS^{c}=(\frac{\Delta_{+}\rho_{i}}{\hat{\rho}_{i+1/2}}+\left(\frac{1}{(\gamma-1)\hat{\beta}_{i+1/2}}-\overline{\|u\|^{2}}_{i+1/2}\right)\Delta_{+}\beta_{i}-2\bar{u}_{i+1/2}\bar{\beta}_{i+1/2}\Delta_{+}u_{i})(\overline{\rho u}_{i+1/2})$
	$\displaystyle+(2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}+2\bar{u}_{i+1/2}\Delta_{+}\beta_{i})\bar{u}_{i+1/2}\overline{\rho u}_{i+1/2}+(2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}+2\bar{u}_{i+1/2}\Delta_{+}\beta_{i})p_{i+1/2}$
	$\displaystyle-2\Delta_{+}\beta_{i}(\frac{1}{2(\gamma-1)\hat{\beta}_{i+1/2}}\overline{\rho u}_{i+1/2}-\frac{\overline{\|u\|^{2}}_{i+1/2}}{2}\overline{\rho u}_{i+1/2}+\overline{\|u\|}^{2}_{i+1/2}\overline{\rho u}_{i+1/2}+p_{i+1/2}\bar{u}_{i+1/2})$	.

Cancelling velocity terms, leads to

	$\displaystyle LHS^{c}=(\frac{\Delta_{+}\rho_{i}}{\hat{\rho}_{i+1/2}}+\left(\frac{1}{(\gamma-1)\hat{\beta}_{i+1/2}}\right)\Delta_{+}\beta_{i})(\overline{\rho u}_{i+1/2})$
	$\displaystyle+(2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}+2\bar{u}_{i+1/2}\Delta_{+}\beta_{i})p_{i+1/2}$
	$\displaystyle-2\Delta_{+}\beta_{i}(\frac{1}{2(\gamma-1)\hat{\beta}_{i+1/2}}\overline{\rho u}_{i+1/2}+p_{i+1/2}\bar{u}_{i+1/2})$	.

Cancelling pressure terms, results in

\displaystyle LHS^{c}=(\frac{\Delta_{+}\rho_{i}}{\hat{\rho}_{i+1/2}})(\overline{\rho u}_{i+1/2})+2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}p_{i+1/2}.

Similarly, the artificial diffusion part of the left-hand side of (98) leads to

\displaystyle LHS^{AD}=\lambda_{i+1/2}\left(\frac{(\Delta_{+}\rho_{i})^{2}}{\hat{\rho}_{i+1/2}}+2\bar{\beta}_{i+1/2}\bar{\rho}_{i+1/2}(\Delta_{+}u_{i})^{2}+\frac{\bar{\rho}_{i+1/2}}{\gamma-1}\Delta_{+}\beta_{i}(D_{+}\frac{1}{\beta_{i}})\right),

in the same way as (97) was obtained. The complete left-hand side of (98), is given by $LHS=LHS^{c}+LHS^{AD}$ . By substituting

\displaystyle\overline{\rho u}_{i+1/2}=\bar{\rho}_{i+1/2}\bar{u}_{i+1/2}+\frac{(\Delta_{+}u_{i})(\Delta_{+}\rho_{i})}{4}=\hat{\rho}_{i+1/2}\bar{u}_{i+1/2}+(\bar{\rho}_{i+1/2}-\hat{\rho}_{i+1/2})\bar{u}_{i+1/2}+\frac{(\Delta_{+}u_{i})(\Delta_{+}\rho_{i})}{4},

we obtain

	$\displaystyle LHS=(\frac{\Delta_{+}\rho_{i}}{\hat{\rho}_{i+1/2}})\left(\hat{\rho}_{i+1/2}\bar{u}_{i+1/2}+(\bar{\rho}_{i+1/2}-\hat{\rho}_{i+1/2})\bar{u}_{i+1/2}+\frac{(\Delta_{+}u_{i})(\Delta_{+}\rho_{i})}{4}\right)+2\bar{\beta}_{i+1/2}\Delta_{+}u_{i}p_{i+1/2}$
	$\displaystyle+\lambda_{i+1/2}(\frac{(\Delta_{+}\rho_{i})^{2}}{\hat{\rho}_{i+1/2}}+2\bar{\beta}_{i+1/2}\bar{\rho}_{i+1/2}(\Delta_{+}u_{i})^{2}+\frac{\bar{\rho}_{i+1/2}}{\gamma-1}\Delta_{+}\beta_{i}(D_{+}\frac{1}{\beta_{i}}))$	.

Use $p_{i+1/2}=\frac{\bar{\rho}_{i+1/2}}{2\bar{\beta}_{i+1/2}}$ (defined in (68)), to obtain

	$\displaystyle LHS=\bar{u}_{i+1/2}(\Delta_{+}\rho_{i})+\bar{\rho}_{i+1/2}(\Delta_{+}u_{i})+\frac{(\Delta_{+}\rho_{i})^{2}}{\hat{\rho}_{i+1/2}}\left(\frac{\bar{\rho}_{i+1/2}-\hat{\rho}_{i+1/2}}{\Delta_{+}\rho_{i}}\bar{u}_{i+1/2}+\frac{\Delta_{+}u_{i}}{4}\right)$
	$\displaystyle+\lambda_{i+1/2}(\frac{(\Delta_{+}\rho_{i})^{2}}{\hat{\rho}_{i+1/2}}+2\bar{\beta}_{i+1/2}\bar{\rho}_{i+1/2}(\Delta_{+}u_{i})^{2}+\frac{\bar{\rho}_{i+1/2}}{\gamma-1}\Delta_{+}\beta_{i}(D_{+}\frac{1}{\beta_{i}}))$	.

The first two terms combine to $\Delta(\rho u)_{i}=\Delta\Psi_{i}$ , which is the right-hand side of (98). The last terms are positive and the $\lambda_{i+1/2}$ bounds $|\frac{\bar{\rho}-\hat{\rho}}{\Delta_{+}\rho_{i}}\bar{u}+\frac{\Delta u}{4}|$ . (In this expression and the next, we have suppressed the indices since all averages and differences are taken between $i,i+1$ .) The latter is realised by the following calculation using (62),

	$\displaystyle\left\|\frac{\bar{\rho}-\hat{\rho}}{\Delta\rho}\right\|=\left\|\frac{\bar{\rho}-\frac{\Delta\rho}{\Delta\log\rho}}{\Delta\rho}\right\|=\left\|\frac{\bar{\rho}\Delta\log\rho-\Delta\rho}{\Delta\log\rho\Delta\rho}\right\|=\left\|\frac{1}{\Delta\log\rho}\left(\frac{\bar{\rho}}{\hat{\rho}}-1\right)\right\|$	$\displaystyle=$
	$\displaystyle\left\|\frac{1}{\Delta\log\rho}\frac{\bar{\rho}-\hat{\rho}}{\hat{\rho}}\right\|\leq\left\|\frac{1}{\Delta\log\rho}\frac{\Delta\rho}{2\hat{\rho}}\right\|=\left\|\frac{\hat{\rho}}{2\hat{\rho}}\right\|$	$\displaystyle=\frac{1}{2},$		(99)

Hence, the fluxes (68),(69) and (70), augmented with the artificial diffusion fluxes, satisfy (98).

Next, we use the above results when contracting the convective terms of the scheme (66) with ${\bf w}_{ijk}^{T}$ , and use (98), to obtain

$\displaystyle\sum_{ijk=0}^{N}{\bf w}_{ijk}^{T}\left((\Delta_{-}^{x}({\bf f}^{c}_{i+1/2jk}-{\bf f}^{\lambda}_{i+1/2jk})){\mathcal{S}}^{x}_{jk}+(\Delta_{-}^{y}({\bf g}_{ij+1/2k}^{c}-{\bf g}_{ij+1/2k}^{\lambda})){\mathcal{S}}^{y}_{ik}\right.$
$\displaystyle\left.+(\Delta_{-}^{z}({\bf h}^{c}_{ijk+1/2}-{\bf h}^{\lambda}_{ijk+1/2})){\mathcal{S}}^{z}_{ij}\right)$	$\displaystyle\leq$
$\displaystyle\sum_{jk=0}^{N}({\bf w}_{Njk}^{T}{\bf f}^{c}_{N+1/2jk}-{\bf w}_{0jk}^{T}{\bf f}^{c}_{-1/2jk}+\Psi^{x}_{0jk}-\Psi^{x}_{Njk}){\mathcal{S}}^{x}_{jk}$
$\displaystyle+\sum_{ik=0}^{N}({\bf w}_{iNk}^{T}{\bf g}_{iN+1/2k}^{c}-{\bf w}_{i0k}^{T}{\bf g}^{c}_{i,-1/2k}+\Psi^{y}_{i0k}-\Psi^{y}_{iNk}){\mathcal{S}}^{y}_{ik}$
$\displaystyle+\sum_{jk=0}^{N}({\bf w}_{ijN}^{T}{\bf h}^{c}_{ijN+1/2}-{\bf w}_{ij0}^{T}{\bf h}^{c}_{ij,-1/2}+\Psi^{z}_{ij0}-\Psi^{z}_{ikN}){\mathcal{S}}^{z}_{ij}$	$\displaystyle=0.$	(100)

All the boundary terms vanish thanks to the no-slip condition implying that $w^{2,3,4}=0$ ; the first and fifth components are zero since the flux components are zero (see (79)); the entropy potentials also vanish due to the no-slip condition.

Combining (100) and (95), we compute the entropy estimate

	$\displaystyle\sum_{ijk=0}^{N}V_{ijk}(U_{ijk})_{t}+\sum_{ijk=0}^{N-1,N,N}(\Delta_{+}^{x}{\bf w}_{ijk})^{T}{\bf f}^{d}_{i+1/2jk}{\mathcal{S}}^{x}_{jk}$			(101)
	$\displaystyle+\sum_{ijk=0}^{N,N-1,N}(\Delta^{y}_{+}{\bf w}_{ijk})^{T}{\bf g}_{ij+1/2k}^{d}{\mathcal{S}}^{y}_{ik}+\sum_{ijk=0}^{N,N,N-1}(\Delta^{z}_{+}{\bf w}_{ijk})^{T}{\bf h}^{d}_{ijk+1/2}{\mathcal{S}}^{z}_{ij}$	$\displaystyle\leq 0,$

where we have used ${\bf w}_{ijk}^{T}({\bf u}_{ijk})_{t}=U({\bf u}_{ijk})_{t}=(U_{ijk})_{t}\quad$ .

Integrating (101) and using (97) and (71), we obtain (among others) the following bounds

$\displaystyle\mathcal{C}\geq$	$\displaystyle\int_{0}^{\mathcal{T}}\sum_{ijk=0}^{N-1,N,N}\nu_{i+1/2jk}\frac{(D_{+}^{x}\rho_{ijk})^{2}}{\hat{\rho}}V_{ijk}\,dt$
$\displaystyle=$	$\displaystyle\int_{0}^{\mathcal{T}}\sum_{ijk=0}^{N-1,N,N}(\mu_{1}\bar{\rho}_{i+1/2jk}+\frac{\mu_{0}}{\hat{\rho}_{i+1/2jk}})\frac{(D_{+}^{x}\rho_{ijk})^{2}}{\hat{\rho}_{i+1/2jk}}V_{ijk}\,dt$
$\displaystyle=$	$\displaystyle\int_{0}^{\mathcal{T}}\sum_{ijk=0}^{N-1,N,N}(\mu_{1}\frac{\bar{\rho}_{i+1/2jk}}{\hat{\rho}_{i+1/2jk}}(D_{+}^{x}\rho_{ijk})^{2}+\mu_{0}(D_{+}^{x}\log\rho_{ijk})^{2})V_{ijk}\,dt,$	(102)

where $\mathcal{C}$ is a constant obtained from the initial data. Noting that, $\bar{\rho}/\hat{\rho}\geq 1$ , we obtain bounds on $D_{+}^{x}\rho_{ijk}$ and $D_{+}^{x}\log\rho_{ijk}$ . Along with (94), we have,

	$\displaystyle\rho_{ijk}$	$\displaystyle\in L^{2}(0,\mathcal{T};H^{1}(\Omega)),$		(103)
	$\displaystyle\nabla_{+}\log\rho_{ijk}$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2}(\Omega)),$

where $\nabla_{+}$ is short for the finite differences $(D_{+}^{x},D_{+}^{y},D_{+}^{z})$

In the same way as in the derivation following (16), and again by using that $\frac{\bar{\rho}_{i+1/2jk}}{\hat{\rho}_{i+1/2jk}}\geq 1$ , we obtain the discrete counterpart of (17) and (18):

$\displaystyle\rho_{ijk}\log(\rho_{ijk})$	$\displaystyle\in L^{\infty}(0,\mathcal{T},L^{1}(\Omega)),$
$\displaystyle\rho_{ijk}\log(T_{ijk})$	$\displaystyle\in L^{\infty}(0,\mathcal{T},L^{1}(\Omega)),$	(104)
$\displaystyle\nabla_{+}\log(T_{ijk})$	$\displaystyle\in L^{2}(0,\mathcal{T},L^{2}(\Omega)).$

(We remark that we get all the counterparts of (18) but limit ourselves to list the estimates necessary for the convergence proof below.)

Next, we turn to the temperature estimate from radiation in (101). We have,

\displaystyle\mathcal{C}\geq\int_{0}^{\mathcal{T}}\sum_{ijk=0}^{N}-V_{ijk}\kappa_{r}2(D_{+}^{x}\beta_{ijk})D_{+}^{x}T^{4}_{ijk}\,dt,

where

\displaystyle D_{+}^{x}\beta_{ijk}=-\frac{D_{+}^{x}T_{ijk}}{2R\check{T}_{i+1/2jk}^{2}},

and

	$\displaystyle D_{+}^{x}(T_{i+1/2jk})^{4}$	$\displaystyle=4\widetilde{T^{3}}_{i+1/2jk}D_{+}^{x}T_{i}$
	$\displaystyle 4\widetilde{T^{3}}_{i+1/2jk}$	$\displaystyle=(T_{i+1}^{3}+T_{i}^{3}+T_{i}^{2}T_{i+1}+T_{i+1}^{2}T_{i})_{jk}.$

Hence, we have a bound on

\displaystyle\frac{\sqrt{\widetilde{T^{3}}_{i+1/2jk}}}{\sqrt{T_{i+1jk}T_{ijk}}}D_{+}^{x}T_{i}\in L^{2}_{2}.

We wish to show that

\displaystyle D_{+}^{x}T^{3/2}_{ijk}\in L^{2}_{2}.

(105)

Hence, we need to demonstrate that

\displaystyle\frac{\sqrt{\widetilde{T^{3}}_{i+1/2jk}}}{\sqrt{T_{i+1jk}T_{ijk}}}\gtrsim\frac{D_{+}^{x}T^{3/2}_{ijk}}{D^{x}_{+}T_{ijk}}.

(106)

We note that $\frac{{\widetilde{T^{3}}_{i+1/2jk}}}{T_{i+1jk}T_{ijk}}\gtrsim T_{i+1jk}+T_{ijk}$ . That is, the left-hand side majorises the arithmetic mean (with an appropriate scaling) and twice the arithmetic mean majorises both $T_{ijk}$ and $T_{i+1jk}$ . Furthermore, by the mean-value theorem

\displaystyle\frac{(D_{+}^{x}T^{3/2}_{ijk})}{(D^{x}_{+}T_{ijk})}=\frac{3}{2}\sqrt{T^{*}},\quad T^{*}\in(T_{i},T_{i+1}),

showing that it is a monotone average. Since it is a monotone average, it is majorised by either end points and therefore by the left-hand side of (106). We conclude that (105), holds. (The same arguments apply in the y- and z-directions.)

By repeating the arguments in the continuous case, we have a $|\mathcal{B}|$ (as in (20)) of non-zero measure and $T_{i}\in L^{\infty}(0,\mathcal{T};L^{1}(\mathcal{B}))$ . Hence,

	$\displaystyle T_{ijk}^{3/2}$	$\displaystyle\in L^{2}(H^{1}),$
	$\displaystyle T_{ijk}$	$\displaystyle\in L^{3}(L^{9}).$		(107)

Finally, the bound on $\frac{(D_{+}^{x}T_{ijk})^{2}}{\check{T}^{2}_{i+1/2jk}}\widetilde{T^{3}_{i+1/2jk}}$ implies

\displaystyle T_{ijk}^{3/2}\nabla_{+}\log T_{ijk}

\displaystyle\in L^{2}(L^{2}).

(108)

6.3 Kinetic energy

Summing (89) by parts over the spatial domain, results in,

	$\displaystyle\sum_{ijk=1}^{n}V_{ijk}\frac{d}{dt}K_{ijk}=\sum_{jk=1}^{N,N}\left({\mathfrak{F}}_{-1/2jk}-{\mathfrak{F}}_{N+1/2jk}\right){\mathcal{S}}^{x}_{jk}$
	$\displaystyle+\sum_{ik=1}^{N,N}\left({\mathfrak{G}}_{i,-1/2k}-{\mathfrak{G}}_{i,N+1/2k}\right){\mathcal{S}}^{y}_{ik}+\sum_{ij=1}^{N,N}\left({\mathfrak{H}}_{ij,-1/2}-{\mathfrak{H}}_{ij,N+1/2}\right){\mathcal{S}}^{z}_{ij}$		(109)
	$\displaystyle-\sum_{ijk=1}^{N}V_{ijk}\left({\mathfrak{D}}_{ijk}-(pDv)_{ijk}\right).$

The fluxes were defined in (90) and vanish at the boundary thanks to no-slip.

As in the continuous case, we split the pressure term, using

\displaystyle-\sum_{ijk=1}^{N-1,N,N}V_{ijk}(D_{+}^{x}u_{ijk})\frac{\bar{\rho}_{i+1/2jk}}{2\bar{\beta}_{i+1/2jk}}\leq\|\rho_{i+1/2jk}D_{+}^{x}u_{ijk}\|^{2}_{2}+\|\frac{1}{2\bar{\beta}_{i+1/2jk}}\|_{2}^{2},

and similarly for the other terms. Since $\nu_{i+1/2jk}$ contains a $\bar{\rho}_{i+1/2jk}$ term (see (71)) the first term on the right-hand side is controlled by $\sum_{ijk=1}^{N}V_{ijk}{\mathfrak{D}}_{ijk}$ . The second term, is controlled by (107) via

\displaystyle\frac{1}{\bar{\beta}_{i+1/2}}=\frac{1}{\frac{1}{2RT_{i+1}}+\frac{1}{2RT_{i}}}=\frac{1}{\frac{T_{i}+T_{i+1}}{2RT_{i+1}T_{i}}}=\frac{RT_{i+1}T_{i}}{\frac{T_{i}+T_{i+1}}{2}}\leq\frac{R\check{T}_{i+1/2}^{2}}{\bar{T}_{i+1/2}}\leq R\bar{T}_{i+1/2},

where the last inequality follows from (62) and positivity of temperature.

Using that $\bar{\rho}/\hat{\rho}\geq 1$ , we conclude that

$\displaystyle\bar{\rho}_{i+1/2jk}\left(\|D_{+}^{x}u_{ijk}\|+\|D_{+}^{x}v_{ijk}\|+\|D_{+}^{x}w_{ijk}\|\right)$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2}),$
$\displaystyle\|D_{+}^{x}u_{ijk}\|+\|D_{+}^{x}v_{ijk}\|+\|D_{+}^{x}w_{ijk}\|$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2}),\quad\textrm{and}$	(110)
$\displaystyle\sqrt{\frac{\bar{\rho}_{i+1/2jk}}{\hat{\rho}_{i+1/2jk}}}\left(\|D_{+}^{x}u_{ijk}\|+\|D_{+}^{x}v_{ijk}\|+\|D_{+}^{x}w_{ijk}\|\right)$	$\displaystyle\in L^{2}(0,\mathcal{T};L^{2}),$

and similarly in the other two directions. These bounds are the discrete equivalents of (24) and also imply that

\displaystyle{\bf v}_{ijk}\in L^{2}(H^{1})\cap L^{2}(L^{6}).

(111)

We will also need the complete estimate obtained from ${\mathfrak{D}}_{ijk}$ :

\displaystyle\tilde{\nu}_{i+1/2jk}\bar{\rho}_{i+1/2jk}\left((D_{+}^{x}u_{ijk})^{2}+(D_{+}^{x}v_{ijk})^{2}+(D_{+}^{x}w_{ijk})^{2}\right)

\displaystyle\in L^{1}(0,\mathcal{T};L^{1}),

(112)

and the corresponding bounds in the other two directions.

6.4 Specific volume

We need the following identities

	$\displaystyle D_{+}^{x}\frac{1}{\rho_{i}^{2}}=\frac{1}{\Delta_{+}^{x}x_{i}}\left(\frac{1}{\rho_{i+1}}+\frac{1}{\rho_{i}}\right)\left(\frac{1}{\rho_{i+1}}-\frac{1}{\rho_{i}}\right)$	$\displaystyle=$
	$\displaystyle-\frac{1}{\Delta_{+}^{x}x_{i}}\frac{\rho_{i+1}+\rho_{i}}{\rho_{i+1}\rho_{i}}\frac{\rho_{i+1}-\rho_{i}}{\rho_{i+1}\rho_{i}}=-\frac{2\bar{\rho}_{i+1/2}}{\check{\rho}_{i+1/2}^{4}}D_{+}^{x}\rho_{i},$

and

\displaystyle D_{+}^{x}\frac{1}{\rho_{i}}=-\frac{1}{\check{\rho}_{i+1/2}^{2}}D_{+}^{x}\rho_{i},

where $\check{\rho}_{i+1/2}$ is given by (61).

Remark 13.

To reduce notation, we take the grid to be uniform in each direction, such that $\Delta_{+}^{x}x_{i}=\Delta_{+}^{x}x_{i+1/2}=h$ . This is not a necessity and all estimates are derivable on a non-uniform grid. In fact, the scheme (and the on-going convergence proof) is readily generalised unstructured to Voronoi grids as well.

Next, multiply the continuity equation of (66) by $-1/\rho_{ijk}^{2}$ ,

	$\displaystyle V_{ijk}(\rho^{-1}_{ijk})_{t}$	$\displaystyle=\frac{1}{\rho_{ijk}^{2}}\left(\Delta_{-}^{x}{\bf f}^{c\lambda,1})_{i+1/2jk}{\mathcal{S}}^{x}_{jk}+\Delta_{-}^{y}({\bf g}^{c\lambda,1})_{ij+1/2k}{\mathcal{S}}_{ik}^{y}+\Delta_{-}^{z}({\bf h}^{c\lambda,1})_{ijk+1/2}{\mathcal{S}}_{ij}^{z}\right)$
		$\displaystyle-\frac{1}{\rho_{ijk}^{2}}\left(\Delta_{-}^{x}(\nu_{i+1/2jk}D_{+}^{x}\rho_{ijk}){\mathcal{S}}_{jk}^{x}+\Delta_{-}^{y}(\nu_{ij+1/2k}D_{+}^{y}\rho_{ijk}){\mathcal{S}}_{ik}^{y}+\Delta_{-}^{z}(\nu_{ijk+1/2}D_{+}^{z}\rho_{ijk}){\mathcal{S}}_{ij}^{z}\right)$		(113)

where ${\bf f}^{c\lambda,1}={\bf f}^{c,1}-{\bf f}^{\lambda,1}$ etc. We rewrite the flux as,

\displaystyle{\bf f}^{c\lambda,1}_{i+1/2jk}=\overline{\rho u}_{i+1/2jk}-\lambda_{i+1/2jk}\Delta_{+}^{x}\rho_{ijk}=\rho^{*}_{i+1/2jk}\bar{u}_{i+1/2jk}-\lambda^{*}_{i+1/2jk}\Delta_{+}^{x}\rho_{ijk}

where $\lambda^{*}=\lambda-\frac{\rho^{*}-\bar{\rho}}{\Delta\rho}\bar{u}+\frac{\Delta u\Delta\rho}{4}\geq 0$ . (Again, we momentarily suppress the indices $ijk$ and $i+1jk$ .) Furthermore, we choose

\displaystyle\rho^{*}=\frac{\check{\rho}^{2}}{\sqrt{\hat{\rho}}\sqrt{\bar{\rho}}},

and calculate

\displaystyle|\rho^{*}-\bar{\rho}|

\displaystyle=\left|\frac{\check{\rho}^{2}}{\sqrt{\hat{\rho}}\sqrt{\bar{\rho}}}-\bar{\rho}\right|=\frac{|\check{\rho}^{2}-\bar{\rho}^{3/2}\sqrt{\hat{\rho}}|}{\sqrt{\hat{\rho}}\sqrt{\bar{\rho}}}.

Since $\check{\rho}\leq\bar{\rho}$ and $\hat{\rho}\leq\bar{\rho}$ , we have

\displaystyle\frac{|\rho^{*}-\bar{\rho}|}{\Delta\rho}\leq\frac{|\check{\rho}^{2}-\bar{\rho}^{2}|}{\sqrt{\hat{\rho}}\sqrt{\bar{\rho}}\Delta\rho}=\frac{|(\Delta\rho)^{2}|}{4\sqrt{\hat{\rho}}\sqrt{\bar{\rho}}\Delta\rho}\leq\frac{|(\Delta\rho)^{2}|}{4\hat{\rho}\Delta\rho}=\frac{1}{4}|\Delta\log\rho|.

In view of (72), the resulting artificial diffusion coefficient is

\displaystyle\lambda^{*}=\lambda-(\frac{\rho^{*}-\bar{\rho}}{\Delta\rho}\bar{u}-\frac{\Delta u}{4})\geq 0.

We introduce the notation, $\nu^{*}_{i+1/2jk}=\nu_{i+1/2jk}+(\Delta_{+}^{x}x_{i})\lambda^{*}_{i+1/2jk}$ in analogy with (73).

Next, we sum the x-terms of (113) by parts and use the boundary conditions. Furthermore, we choose $\eta$ such that $\mu_{0}/2>\eta>0$ , and obtain

	$\displaystyle\sum_{i=0}^{N}{\mathcal{S}}_{jk}^{x}\left(\frac{1}{\rho_{ijk}^{2}}\Delta_{-}^{x}(\rho^{}\bar{u})_{i+1/2jk}-\frac{1}{\rho_{ijk}^{2}}\Delta_{-}^{x}(\nu^{}_{i+1/2jk}D_{+}^{x}\rho_{ijk})\right)$	$\displaystyle=$
	$\displaystyle\sum_{i=0}^{N-1}{\mathcal{S}}_{jk}^{x}\Delta^{x}_{+}x_{i}\left(\frac{2\bar{\rho}_{i+1/2jk}}{\check{\rho}_{i+1/2jk}^{4}}(D_{+}^{x}\rho_{ijk})\rho^{}_{i+1/2jk}\bar{u}_{i+1/2jk}-\nu^{}_{i+1/2jk}\frac{2\bar{\rho}_{i+1/2jk}}{\check{\rho}_{i+1/2jk}^{4}}(D_{+}^{x}\rho_{ijk})^{2}\right)$	$\displaystyle\leq$
	$\displaystyle\sum_{i=0}^{N-1}V_{ijk}\left(-\frac{2\bar{\rho}_{i+1/2jk}\rho^{}_{i+1/2jk}}{\check{\rho}_{i+1/2jk}^{2}}(D_{+}^{x}\frac{1}{\rho_{ijk}})\bar{u}_{i+1/2jk}-\nu^{}_{i+1/2jk}2\bar{\rho}_{i+1/2jk}(D_{+}^{x}\frac{1}{\rho_{ijk}})^{2}\right)$	$\displaystyle\leq$
	$\displaystyle\sum_{i=0}^{N-1}V_{ijk}\left(\eta\|\frac{2\bar{\rho}_{i+1/2}\rho^{*}_{i+1/2jk}}{\check{\rho}_{i+1/2jk}^{2}}D_{+}^{x}\frac{1}{\rho_{ijk}}\|^{2}+\eta^{-1}\|\bar{u}_{i+1/2jk}\|^{2}-\frac{2\mu_{0}\bar{\rho}_{i+1/2jk}}{\hat{\rho}_{i+1/2jk}}\|D_{+}^{x}\frac{1}{\rho_{ijk}}\|^{2}\right)$	,

where we have used (71) in the last row. Since $\frac{\bar{\rho}_{ijk}}{\check{\rho}_{ijk}}>1$ , and by summing over $jk$ , we obtain the bounds

	$\displaystyle\sup_{0,\mathcal{T}}\\|\rho^{-1}_{ijk}\\|_{1}\leq\mathcal{C},$	$\displaystyle\quad\int_{0}^{\mathcal{T}}\\|\nabla_{+}\rho^{-1}_{ijk}\\|_{2}^{2}\,dt,\leq\mathcal{C},$		(114)
	$\displaystyle\rho_{ijk}^{-1}\in L^{2}(H^{1}),$	$\displaystyle\quad\rho_{ijk}^{-1}\in L^{2}(L^{6}).$

(Recall that $\nabla_{+}$ is short for the finite differences $(D_{+}^{x},D_{+}^{y},D_{+}^{z})$ .) These estimates correspond to (28)-(29)). Furthermore, we get

\displaystyle\sqrt{T_{ijk}}\in L^{\infty}_{1},

(115)

as in (30).

6.5 Estimates from the continuity equation

Multiply the continuity equation by $\rho_{ijk}$ and sum in space.

	$\displaystyle\sum_{i,j,k=0}^{N}\frac{1}{2}V_{ijk}(\rho^{2}_{ijk})_{t}=$
	$\displaystyle\sum_{i,j,k=0}^{N}-\rho_{ijk}V_{ijk}\left(D_{-}^{x}({\bf f}^{c,1})_{i+1/2jk}+D_{-}^{y}({\bf g}^{c,1})_{ij+1/2k}+D_{-}^{z}({\bf h}^{c,1})_{ijk+1/2}\right)$
	$\displaystyle+\sum_{i,j,k=0}^{N}\rho_{ijk}V_{ijk}\left(D_{-}^{x}\tilde{\nu}_{i+1/2jk}D_{+}^{x}\rho_{ijk}+D_{-}^{y}\tilde{\nu}_{ij+1/2k}D_{+}^{y}\rho_{ijk}+D_{-}^{z}\tilde{\nu}_{ijk+1/2}D_{+}^{z}\rho_{ijk}\right).$

We sum by parts and thanks to the boundary conditions, the boundary terms vanish. (Again, all three directions are handled analogously and we focus on the x-direction and denote the remaining terms by $YZ$ )

\displaystyle\sum_{i,j,k=0}^{N}\frac{1}{2}V_{ijk}(\rho^{2}_{ijk})_{t}

\displaystyle=\sum_{i,j,k=0}^{N-1,N,N}\left((D_{+}^{x}\rho_{ijk})\overline{\rho u}_{i+1/2jk}-(D_{+}^{x}\rho_{ijk})\tilde{\nu}_{i+1/2jk}D_{+}^{x}\rho_{ijk}\right)+YZ.

Use $\overline{\rho u}_{i+1/2}=\sqrt{\rho_{i}}\sqrt{\rho_{i}}u_{i}+\sqrt{\rho_{i+1}}\sqrt{\rho_{i+1}}u_{i+1}$ which allows us to proceed as in the continuous case. We have $L^{\infty}_{2}$ bounds of $\sqrt{\rho_{i}}u_{i}$ . We need the diffusive term to provide a bound $\sqrt{\rho_{i}}D_{+}^{x}\rho_{i}\in L^{2}_{2}$ and $\sqrt{\rho_{i+1}}D_{+}^{x}\rho_{i}$ , which is what we get since $\tilde{\nu}_{i+1/2}\sim\bar{\rho}_{i+1/2}$ and $\sqrt{\bar{\rho}_{i+1/2}}\geq\frac{1}{2}(\sqrt{\rho_{i}}+\sqrt{\rho_{i+1}})$ . We collect the estimates,

	$\displaystyle\rho_{ijk}\in L^{\infty}_{2},$		(116)
	$\displaystyle\sqrt{\bar{\rho}_{i+1/2jk}}D_{+}^{x}\rho_{ijk}\in L^{\infty}_{2}.$

Next, we test the equation with $\rho_{i}^{2}$ . The procedure follows the same route as above and we turn directly to the convective terms. In the x-direction, we have

\displaystyle\sum_{i=0}^{N}\rho_{i}^{2}D_{-}^{x}{\bf f}^{c\lambda,1}_{i+1/2}=-\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{2}){\bf f}^{c\lambda,1}_{i+1/2},

where we have dropped the sum over $jk$ and suppressed those indices. As above, we denote ${\bf f}^{c\lambda,1}={\bf f}^{c,1}-{\bf f}^{\lambda,1}$ . The boundary terms vanish thanks to the boundary conditions. Some further manipulations yield,

	$\displaystyle\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{2}){\bf f}^{c\lambda,1}_{i+1/2}$	$\displaystyle=$
	$\displaystyle\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{2})\left(\rho^{}_{i+1/2}\bar{u}_{i+1/2}+\underbrace{\frac{\bar{\rho}_{i+1/2}-\rho^{}_{i+1/2}}{\Delta_{+}^{x}\rho_{i}}\bar{u}_{i+1/2}\Delta_{+}^{x}\rho_{i}+\frac{\Delta_{+}^{x}u_{i}}{4}\Delta_{+}^{x}\rho_{i}-\lambda_{i+1/2}\Delta_{+}^{x}\rho_{i}}_{RT_{i+1/2}}\right).$

Here, we choose the average:

\displaystyle\rho^{*}_{i+1/2}=\frac{2}{3}\frac{D_{+}^{x}\rho_{i}^{3}}{D_{+}^{x}\rho_{i}^{2}}.

Noting that $(D_{+}^{x}\rho_{i}^{2})=2\bar{\rho}_{i+1/2}D_{+}^{x}\rho_{i}$ , the rest terms ( $RT_{i+1/2}$ ) are bounded since

\displaystyle\lambda_{i+1/2}\geq|\frac{\bar{\rho}_{i+1/2}-\rho^{*}_{i+1/2}}{\Delta_{+}^{x}\rho_{i}}\bar{u}_{i+1/2}+\frac{\Delta_{+}^{x}u_{i}}{4}|.

The last inequality follows from the following calculation,

	$\displaystyle\frac{\bar{\rho}_{i+1/2}-\rho^{*}_{i+1/2}}{\Delta_{+}^{x}\rho_{i}}=\frac{\bar{\rho}_{i+1/2}-\frac{2}{3}\frac{\rho_{i+1}^{3}-\rho_{i}^{3}}{\rho_{i+1}^{2}-\rho_{i}^{2}}}{\Delta_{+}^{x}\rho_{i}}$	$\displaystyle=$
	$\displaystyle\frac{\bar{\rho}_{i+1/2}-\frac{2(\rho_{i+1}^{2}+\rho_{i}^{2}+\rho_{i+1}\rho_{i})}{3(\rho_{i+1}+\rho_{i})}}{\Delta_{+}^{x}\rho_{i}}=\frac{\frac{3(\rho_{i+1}+\rho_{i})^{2}-4(\rho_{i+1}^{2}+\rho_{i}^{2}+\rho_{i+1}\rho_{i})}{6(\rho_{i+1}+\rho_{i})}}{\Delta_{+}^{x}\rho_{i}}$	$\displaystyle=$
	$\displaystyle\frac{-\rho_{i+1}^{2}-\rho_{i}^{2}+2\rho_{i+1}\rho_{i}}{6(\rho_{i+1}+\rho_{i})\Delta_{+}^{x}\rho_{i}}=\frac{-(\rho_{i+1}-\rho_{i})^{2}}{6(\rho_{i+1}+\rho_{i})\Delta_{+}^{x}\rho_{i}}$	$\displaystyle=$
	$\displaystyle-\frac{(\rho_{i+1}-\rho_{i})}{6(\rho_{i+1}+\rho_{i})}$	,

and (72).

Using (92),

\displaystyle\sum_{i=0}^{N}\rho_{i}^{2}D_{-}^{x}{\bf f}^{c\lambda,1}_{i+1/2}=-\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{2}){\bf f}^{c\lambda,1}_{i+1/2}=-\frac{2}{3}\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3})\bar{u}_{i+1/2}+(D_{+}^{x}\rho_{i}^{2})RT_{i+1/2},

where $RT$ denotes a damping term from the rest terms and artificial diffusion and we have used $u_{0jk}=u_{Njk}=0$ in the summation-by-parts step. We drop the rest terms, as they have the correct sign and recast the sums using the no-slip condition ( $u_{0}=u_{N}=0$ ),

$\displaystyle\frac{2}{3}\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3})\bar{u}_{i+1/2}=\frac{2}{3}\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3})\frac{u_{i+1}+u_{i}}{2}=\frac{1}{3}\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3})u_{i+1}+\frac{1}{3}\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3})u_{i}$
$\displaystyle\frac{1}{3}\sum_{i=1}^{N-1}(D_{+}^{x}\rho_{i-1}^{3})u_{i}+\frac{1}{3}\sum_{i=1}^{N-1}(D_{+}^{x}\rho_{i}^{3})u_{i}=-\frac{1}{3}\sum_{i=1}^{N}\rho_{i-1}^{3}D_{-}^{x}u_{i}-\frac{1}{3}\sum_{i=1}^{N}\rho_{i}^{3}D_{-}^{x}u_{i}$	$\displaystyle=$
$\displaystyle-\frac{2}{3}\sum_{i=1}^{N}\overline{\rho^{3}}_{i-1/2}D_{-}^{x}u_{i}=-\frac{2}{3}\sum_{i=0}^{N-1}\overline{\rho^{3}}_{i+1/2}D_{+}^{x}u_{i}$	.	(117)

Since $\bar{\rho}_{i+1/2jk}D_{+}^{x}u_{ijk}\in L^{2}_{2}$ by (110), we can proceed as in the continuous derivation. First,

\displaystyle\rho_{i+1}^{3}+\rho_{i}^{3}=(\rho_{i+1}+\rho_{i})(\rho_{i+1}^{2}-\rho_{i}\rho_{i+1}+\rho_{i}^{2}).

We can now split the term $|\overline{\rho^{3}}_{i+1/2}D_{+}^{x}u_{i}|\leq\eta(\rho_{i+1}^{2}-\rho_{i}\rho_{i+1}+\rho_{i}^{2})^{2}+\eta^{-1}((\rho_{i+1}+\rho_{i})D_{+}^{x}u_{i})^{2}$ Hence, we need $\rho_{i+1}^{2}-\rho_{i}\rho_{i+1}+\rho_{i}^{2}\in L^{2}_{2}$ . This follows if $\rho_{i}\in L^{4}_{4}$ , which we will get from the diffusive term as in the continuous case.

To this end, we need to recast the diffusive term to bound the gradient of $\rho_{ijk}^{2}$ . In the x-direction, we have

	$\displaystyle\sum_{i=0}^{N}\rho_{i}^{2}D_{-}^{x}(\nu_{i+1/2}D_{+}^{x}\rho_{i})=-\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{2})\nu_{i+1/2}D_{+}^{x}\rho_{i}$	$\displaystyle=$
	$\displaystyle-\sum_{i=0}^{N-1}2\bar{\rho}_{i+1/2}(D_{+}^{x}\rho_{i})\nu_{i+1/2}D_{+}\rho_{i}$	.

Since $\nu_{i+1/2}\sim\bar{\rho}_{i+1/2}$ and seeing that $\bar{\rho}_{i+1/2}D_{+}\rho_{i}=\frac{1}{2}D_{+}\rho_{i}^{2}$ , we obtain the necessary estimate to bound $\rho_{i}\in L^{4}_{4}$ . We obtain all the estimates of (33):

\displaystyle\rho_{ijk}\in L^{\infty}_{3},\quad\quad\rho_{ijk}^{2}\in L^{2}(H^{1}),\quad\quad\rho_{ijk}\in L^{4}_{12},\quad\quad\rho_{ijk}\in L^{6}_{6}.

(118)

Next, we test against $\rho_{i}^{3}$ . As above, we begin with the convective term (and drop the common $jk$ indices),

\displaystyle\sum_{i=0}^{N}\rho_{i}^{3}D_{-}^{x}{\bf f}^{c\lambda,1}_{i+1/2}=-\sum_{i=0}^{N-1}D_{+}^{x}\rho_{i}^{3}{\bf f}^{c\lambda,1}_{i+1/2},

and rewrite the flux as,

\displaystyle{\bf f}^{c\lambda,1}_{i+1/2}=\rho^{*}_{i+1/2}\bar{u}_{i+1/2}+(\bar{\rho}_{i+1/2}-\rho^{*}_{i+1/2})\bar{u}_{i+1/2}-\lambda^{a}_{i+1/2}\Delta^{x}_{+}\rho_{i}.

(119)

This time we choose,

\displaystyle\rho^{*}_{i+1/2}=\frac{3}{4}\frac{\Delta^{x}_{+}\rho_{i}^{4}}{\Delta^{x}_{+}\rho_{i}^{3}}=\frac{3}{4}\frac{(\rho_{i+1}^{2}+\rho_{i}^{2})(\rho_{i+1}+\rho_{i})}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}.

(120)

Since, by the mean value theorem

\displaystyle\frac{\Delta^{x}_{+}\rho_{i}^{4}}{\Delta^{x}_{+}\rho_{i}^{3}}=\frac{(\rho^{3}_{i+1})^{4/3}-(\rho_{i})^{4/3}}{\rho_{i+1}^{3}-\rho_{i}^{3}}=\frac{4}{3}(\rho_{\star}^{3})^{1/3},

for some $\rho_{\star}\in(\rho_{i},\rho_{i+1})$ , we conclude that $\rho^{*}_{i+1/2}\in(\rho_{i},\rho_{i+1})$ . (That is, it is a monotone average.)

As before, we must verify that the artificial diffusion term dominates the error term in (119). That is, the following expression must be positive:

	$\displaystyle\lambda^{a}_{i+1/2}-\frac{\bar{\rho}_{i+1/2}-\rho^{*}_{i+1/2}}{\Delta^{+}_{x}\rho_{i}}\bar{u}_{i+1/2}$	$\displaystyle=$
	$\displaystyle\lambda_{i+1/2}-\frac{\Delta_{+}^{x}u_{i}}{4}-\frac{\bar{\rho}_{i+1/2}-\rho^{*}_{i+1/2}}{\Delta^{x}_{+}\rho_{i}}\bar{u}_{i+1/2}$	.		(121)

Using (120), we calculate

	$\displaystyle\bar{\rho}_{i+1/2}-\rho^{*}_{i+1/2}=\bar{\rho}_{i+1/2}-\frac{3}{4}\frac{(\rho_{i+1}^{2}+\rho_{i}^{2})(\rho_{i+1}+\rho_{i})}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}$	$\displaystyle=$
	$\displaystyle\bar{\rho}_{i+1/2}\left(1-\frac{3}{2}\frac{(\rho_{i+1}^{2}+\rho_{i}^{2})}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}\right)=\bar{\rho}_{i+1/2}\left(\frac{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}-\frac{3}{2}(\rho_{i+1}^{2}+\rho_{i}^{2})}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}\right)$	$\displaystyle=$
	$\displaystyle-\frac{\bar{\rho}_{i+1/2}}{2}\left(\frac{(\rho_{i+1}-\rho_{i})^{2}}{\rho_{i+1}^{2}+\rho_{i+1}\rho_{i}+\rho_{i}^{2}}\right).$

Using the last expression and the definition of the artificial diffusion coefficient (72), it is clear that (121) is positive.

Next, we recast convective term. Using (120) and the analogous derivation as (117), we get

\displaystyle-\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3})\rho^{*}_{i+1/2}\bar{u}_{i+1/2}=-\frac{3}{4}\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{4})\bar{u}_{i+1/2}=\frac{3}{4}\sum_{i=0}^{N-1}\overline{\rho^{4}}_{i+1/2}D_{+}^{x}u_{i}.

We proceed as in the previous estimate and factorise $\overline{\rho^{4}}_{i+1/2}=\bar{\rho}_{i+1/2}\frac{\overline{\rho^{4}}_{i+1/2}}{\bar{\rho}_{i+1/2}}$ and use (110) to bound one part. The other is bounded since $\frac{\overline{\rho^{4}}_{i+1/2}}{\bar{\rho}_{i+1/2}}$ can be majorised by $\max_{j=i,i+1}|\rho_{j}^{3}|$ , which in turn is bounded in $L^{2}_{2}$ by (118).

Turning to the diffusive terms, we have

	$\displaystyle\sum_{i=0}^{N}\rho_{i}^{3}D_{+}^{x}{\bf f}^{d,1}_{i+1/2}=-\sum_{i=0}^{N-1}(D_{+}^{x}\rho_{i}^{3}){\bf f}^{d,1}_{i+1/2}=-\sum_{i=0}^{N-1}D_{+}^{x}\rho_{i}^{3}\nu_{i+1/2}D_{+}\rho_{i}$	$\displaystyle=$
	$\displaystyle-\sum_{i=0}^{N-1}(\rho^{2}_{i+1}+\rho_{i}\rho_{i+1}+\rho_{i}^{2})\nu_{i+1/2}(D_{+}\rho_{i})^{2}$	.		(122)

We use that $\nu_{i+1/2}\sim\rho_{i+1/2}$ , which implies that $(\rho^{2}_{i+1}+\rho_{i}\rho_{i+1}+\rho_{i}^{2})\nu_{i+1/2}$ is dominating $c\cdot(\max(\rho_{i},\rho_{i+1}))$ for some $c>0$ . Hence, (122) dominates $D_{+}^{x}\rho_{i}^{5/2}$ since as, by the same argument as above

\displaystyle\frac{D_{+}^{x}\rho_{i}^{5/2}}{D_{+}^{x}\rho_{i}}=\frac{5}{2}\rho_{\star}^{3/2},

is an average and by the mean-value theorem, $\rho_{\star}\in(\rho_{i},\rho_{i+1})$ . Hence, we obtain

\displaystyle\rho\in L^{\infty}_{4},\quad\nabla_{+}\rho^{5/2}\in L^{2}_{2}\quad\rho\in L^{5}_{15},\quad\rho\in L^{23/3}_{23/3}.

(123)

6.6 Renormalised internal energy

The internal energy balance was derived in (91) and is repeated here for convenience,

$\displaystyle D_{t}\frac{p_{ijk}}{\gamma-1}+(pD{\bf v})_{ijk}-{\mathfrak{D}}_{ijk}$
$\displaystyle+D_{-}^{x}(\frac{\overline{\rho u}_{i+1/2jk}}{2(\gamma-1)\hat{\beta}_{i+1/2jk}})+D_{-}^{y}(\frac{\overline{\rho v}_{ij+1/2k}}{2(\gamma-1)\hat{\beta}_{ij+1/2k}})+D_{-}^{y}(\frac{\overline{\rho w}_{ijk+1/2}}{2(\gamma-1)\hat{\beta}_{ijk+1/2}})$	$\displaystyle=$	(124)
$\displaystyle D_{-}^{x}(\tilde{\nu}_{i+1/2jk}\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1})+D_{-}^{y}(\tilde{\nu}_{ij+1/2k}\frac{({\mathfrak{p}}_{y})_{ij+1/2k}}{\gamma-1})+$
$\displaystyle D_{-}^{z}(\tilde{\nu}_{ijk+1/2}\frac{({\mathfrak{p}}_{z})_{ijk+1/2}}{\gamma-1})$
$\displaystyle+\kappa_{r}\left(D_{-}^{x}D^{x}_{+}T^{4}_{ijk}+D_{-}^{y}D^{y}_{+}T^{4}_{ijk}+D_{-}^{z}\kappa_{r}D^{z}_{+}T^{4}_{ijk}\right)$	.

As in the continuous case, we multiply by $H^{\prime}(T_{i})=H^{\prime}_{i}$ and sum in space. Since $T\geq 0$ , we have

	$\displaystyle H(T)=(1+T)^{1-\omega}<1+T,\,\,\omega>0,$
	$\displaystyle 0<H^{\prime}(T)<1,\quad H^{\prime\prime}(T)<0,\quad\|H^{\prime\prime}\|<c<1.$

We carry out the calculations for a few of the terms separately and begin with the temporal term of (124):

\displaystyle H^{\prime}(T_{i})D_{t}\frac{p_{ijk}}{\gamma-1}=c_{v}(\rho_{ijk}H_{ijk})_{t}+c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})\partial_{t}\rho_{ijk}.

Next, we multiply by $V_{i}$ and sum. Using the continuity equation and the boundary conditions/fluxes, we obtain

$\displaystyle\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})D_{t}\frac{p_{ijk}}{\gamma-1}$	$\displaystyle=$
$\displaystyle\sum_{ijk=0}^{N}V_{ijk}c_{v}(\rho_{ijk}H_{ijk})_{t}+c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})(-\nabla_{-}\cdot({\bf f}^{c,1}_{i+1/2jk},{\bf g}^{c,1}_{ij+1/2jk},{\bf h}^{c,1}_{ijk+1/2})$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})(\nabla_{-}\cdot({\bf f}^{d,1}_{i+1/2jk},{\bf g}^{d,1}_{ij+1/2jk},{\bf h}^{d,1}_{ijk+1/2})$	,	(125)

where $\nabla_{-/+}=(D_{-/+}^{x},D_{-/+}^{y},D_{-/+}^{z})$ .

We adopt the short-hand notation ${\bf f}^{T}_{i+1/2jk}=\frac{\overline{\rho u}_{i+1/2jk}}{2(\gamma-1)\hat{\beta}_{i+1/2jk}}$ etc. for the discrete thermal fluxes and define the corresponding discrete renormalised thermal flux as, ${\bf f}^{RT}_{i+1/2jk}=c_{v}\frac{H(T_{i+1jk})+H(T_{ijk})}{2}\overline{\rho u}_{i+1/2jk}$ .

Remark 14.

As above, we adopt the convention that values “outside” the domain are identical with the boundary value. That is, ${\bf f}^{RT}_{-1/2jk}=c_{v}\frac{H(T_{0jk})+H(T_{-1jk})}{2}\overline{\rho u}_{-1/2jk}=c_{v}H(T_{0})(\rho u)_{0jk}=0$ .

Multiplying the first convective thermal flux of (124) by $V_{ijk}H^{\prime}(T_{ijk})$ , yields

$\displaystyle V_{ijk}H^{\prime}(T_{ijk})D_{-}^{x}(\frac{\overline{\rho u}_{i+1/2jk}}{2(\gamma-1)\hat{\beta}_{i+1/2jk}})=V_{ijk}H^{\prime}(T_{ijk})D_{-}^{x}({\bf f}^{T}_{i+1/2jk})$	$\displaystyle=$
$\displaystyle V_{ijk}D_{-}^{x}{\bf f}^{RT}_{i+1/2jk}$
$\displaystyle-V_{ijk}c_{v}D_{-}^{x}\frac{H(T_{i+1jk})+H(T_{ijk})}{2}\overline{\rho u}_{i+1/2jk}+V_{ijk}H^{\prime}(T_{ijk})\frac{({\bf f}^{T}_{i+1/2jk}-{\bf f}^{T}_{i-1/2jk})}{\Delta_{-}^{x}x_{i+1/2}}$	$\displaystyle=$
$\displaystyle V_{ijk}D_{-}^{x}{\bf f}^{RT}_{i+1/2jk}+E^{x}_{ijk},$		(126)

where

	$\displaystyle E^{x}_{ijk}=-\frac{V_{ijk}}{\Delta_{-}^{x}x_{i+1/2}}\left(c_{v}\frac{H(T_{i+1jk})+H(T_{ijk})}{2}-\frac{H^{\prime}(T_{ijk})}{2(\gamma-1)\hat{\beta}_{i+1/2jk}}\right)\overline{\rho u}_{i+1/2jk}$		(127)
	$\displaystyle+\frac{V_{ijk}}{\Delta_{-}^{x}x_{i+1/2}}\left(c_{v}\frac{H(T_{ijk})+H(T_{i-1jk})}{2}-\frac{H^{\prime}(T_{ijk})}{2(\gamma-1)\hat{\beta}_{i-1/2jk}}\right)\overline{\rho u}_{i-1/2jk}.$

We carry out the same operations for the thermal fluxes in the $y,z$ -directions resulting in the analogous definition of $E^{y,z}_{ijk}$ .

Next, we multiply (124) by $V_{ijk}H^{\prime}_{ijk}$ , sum in space and use (125) and (126) (along with the similar expressions obtained in the other two directions) to obtain

$\displaystyle\sum_{ijk=0}^{N}V_{ijk}\left(c_{v}(\rho_{ijk}H_{ijk})_{t}+D_{-}^{x}{\bf f}^{RT}_{i+1/2jk}+D_{-}^{x}{\bf g}^{RT}_{ij+1/2k}+D_{-}^{z}{\bf h}^{RT}_{ijk+1/2}\right)+\sum_{ijk=0}^{N}(E^{x}_{ijk}+E^{y}_{ijk}+E^{z}_{ijk})$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})\nabla_{-}\cdot({\bf f}^{d,1}_{i+1/2jk},{\bf g}^{d,1}_{ij+1/2jk},{\bf h}^{d,1}_{ijk+1/2})$
$\displaystyle-\sum_{ijk=0}^{N}V_{ijk}c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})\nabla_{-}\cdot({\bf f}^{c,1}_{i+1/2jk},{\bf g}^{c,1}_{ij+1/2jk},{\bf h}^{c,1}_{ijk+1/2})$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})((pD{\bf v})_{ijk}-{\mathfrak{D}}_{ijk})$	$\displaystyle=$	(128)
$\displaystyle\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})\nabla_{-}\cdot\left(\tilde{\nu}_{i+1/2jk}\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1},\tilde{\nu}_{ij+1/2k}\frac{({\mathfrak{p}}_{y})_{ij+1/2k}}{\gamma-1},\tilde{\nu}_{ijk+1/2}\frac{({\mathfrak{p}}_{z})_{ijk+1/2}}{\gamma-1}\right)$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})\kappa_{r}\left(D_{-}^{x}D^{x}_{+}T^{4}_{ijk}+D_{-}^{y}D^{y}_{+}T^{4}_{ijk}+D_{-}^{z}D^{z}_{+}T^{4}_{ijk}\right).$

Combining the error terms $(E^{x,y,z}_{ijk})$ and the convective “c,1” terms lead to,

$\displaystyle\sum_{ijk=0}^{N}V_{ijk}\left(c_{v}(\rho_{ijk}H_{ijk})_{t}+D_{-}^{x}{\bf f}^{RT}_{i+1/2jk}+D_{-}^{x}{\bf g}^{RT}_{ij+1/2k}+D_{-}^{z}{\bf h}^{RT}_{ijk+1/2}\right)+\sum_{ijk=0}^{N}({\mathcal{E}}^{x}_{ijk}+{\mathcal{E}}^{y}_{ijk}+{\mathcal{E}}^{z}_{ijk})$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})\nabla_{-}\cdot({\bf f}^{d,1}_{i+1/2jk},{\bf g}^{d,1}_{ij+1/2jk},{\bf h}^{d,1}_{ijk+1/2})$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})((pD{\bf v})_{ijk}-{\mathfrak{D}}_{ijk})$	$\displaystyle=$	(129)
$\displaystyle\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})\nabla_{-}\cdot\left(\tilde{\nu}_{i+1/2jk}\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1},\tilde{\nu}_{ij+1/2k}\frac{({\mathfrak{p}}_{y})_{ij+1/2k}}{\gamma-1}),\tilde{\nu}_{ijk+1/2}\frac{({\mathfrak{p}}_{z})_{ijk+1/2}}{\gamma-1}\right)$
$\displaystyle+\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})\kappa_{r}\left(D_{-}^{x}D^{x}_{+}T^{4}_{ijk}+D_{-}^{y}D^{y}_{+}T^{4}_{ijk}+D_{-}^{z}D^{z}_{+}T^{4}_{ijk}\right),$

where the error terms are given by

	$\displaystyle{\mathcal{E}}^{x}_{ijk}=-\frac{c_{v}V_{ijk}}{\Delta_{-}^{x}x_{i+1/2}}\left(\frac{H(T_{i+1jk})+H(T_{ijk})}{2}-\frac{H^{\prime}(T_{ijk})}{\widehat{(T^{-1})}_{i+1/2jk}}-(H(T_{ijk})-H^{\prime}(T_{ijk})T_{ijk})\right)\overline{\rho u}_{i+1/2jk}$		(130)
	$\displaystyle+\frac{c_{v}V_{ijk}}{\Delta_{-}^{x}x_{i+1/2}}\left(\frac{H(T_{ijk})+H(T_{i-1jk})}{2}-\frac{H^{\prime}(T_{ijk})}{\widehat{T^{-1}}_{i-1/2jk}}-(H(T_{ijk})-H^{\prime}(T_{ijk})T_{ijk})\right)\overline{\rho u}_{i-1/2jk},$

and analogously for ${\mathcal{E}}^{y,z}_{ijk}$ . Since $H^{\prime}$ is bounded, we have that

\displaystyle\frac{1}{\Delta_{-}^{x}x_{i+1/2}}|\frac{H(T_{i+1jk})+H(T_{ijk})}{2}-H(T_{ijk})|\leq C|D^{+}_{x}T_{ijk}|.

By (62), $\frac{1}{\Delta_{-}^{x}x_{i+1/2}}|\frac{1}{\widehat{T^{-1}}_{i+1/2jk}}-T_{ijk}|\leq|D^{+}_{x}T_{ijk}|$ . We conclude that,

\displaystyle|\sum_{ijk=0}^{N}{\mathcal{E}}^{x}_{ijk}|\leq\mathcal{C}\sum_{ijk=0}^{N}|(D_{x}^{+}T_{ijk})\overline{\rho u}_{i+1/jk}|\lesssim\|D_{+}T\|_{2}\|\overline{\rho u}\|_{2}.

Since $D^{+}_{x}T_{ijk}\in L^{2}(H^{1})$ by (104) and (107), and $\overline{\rho u}_{i+1/2jk}\in L^{2}_{2}$ by (123) and (111), we conclude that the error term $|\sum{\mathcal{E}}^{x}_{ijk}|$ is bounded. By repeating the arguments in the $y,z$ -directions, the remaining two error terms are also bounded.

Next, we consider the terms

\displaystyle\sum_{ijk=0}^{N}V_{ijk}H^{\prime}(T_{ijk})((pD{\bf v})_{ijk}-{\mathfrak{D}}_{ijk}),

in (129). Noting that $H^{\prime}$ is bounded, the pressure term is controlled just like in (109) and the viscous terms are bounded by the estimate obtained by (109).

We must also control the diffusive terms in (129). As usual, we present the arguments for the $x$ -term. Namely, we use

\displaystyle({\mathfrak{p}}_{x})_{i+1/2jk}

\displaystyle=\frac{1}{2\hat{\beta}_{i+1/2jk}}D_{+}^{x}\rho_{ijk}+\frac{\bar{\rho}_{i+1/2jk}}{2}D_{+}^{x}\frac{1}{\beta_{ijk}},

in the diffusive terms,

$\displaystyle\sum_{ijk=0}^{N}V_{ijk}\left(H^{\prime}(T_{ijk})D_{-}^{x}\left(\tilde{\nu}_{i+1/2jk}\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1}\right)-c_{v}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})D_{-}^{x}{\bf f}^{d,1}_{i+1/2jk}\right)$	$\displaystyle=$
$\displaystyle-\sum_{ijk=0}^{N-1,N,N}V_{ijk}(D_{+}^{x}H^{\prime}(T_{ijk}))\left(\tilde{\nu}_{i+1/2jk}\frac{1}{\gamma-1}\left(\frac{1}{2\hat{\beta}_{i+1/2jk}}D_{+}^{x}\rho_{ijk}+\frac{\bar{\rho}_{i+1/2jk}}{2}D_{+}^{x}\frac{1}{\beta_{ijk}}\right)\right)$
$\displaystyle+\sum_{ijk=0}^{N-1,N,N}V_{ijk}c_{v}(D_{+}^{x}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})){\bf f}^{d,1}_{i+1/2jk}$	$\displaystyle=$
$\displaystyle-\sum_{ijk=0}^{N-1,N,N}V_{ijk}(D_{+}^{x}H^{\prime}(T_{ijk}))\left(\tilde{\nu}_{i+1/2jk}\frac{1}{\gamma-1}\left(\frac{\bar{\rho}_{i+1/2jk}}{2}D_{+}^{x}\frac{1}{\beta_{ijk}}\right)\right)$		(131)
$\displaystyle+\sum_{ijk=0}^{N-1,N,N}V_{ijk}c_{v}\left((D_{+}^{x}(H^{\prime}_{ijk}T_{ijk}-H_{ijk}))-\frac{D_{+}^{x}H^{\prime}(T_{ijk})}{c_{v}(\gamma-1)2\hat{\beta}_{i+1/2jk}}\right){\bf f}^{d,1}_{i+1/2jk}$	,

where we have summed by parts and used the boundary conditions/fluxes.

Using,

\displaystyle\frac{1}{\hat{\beta}_{i+1/2}}=\frac{\Delta_{+}\log\beta_{i}}{\Delta\beta_{i}}=2R\frac{-\Delta_{+}\log T_{i}}{\Delta_{+}T_{i}^{-1}}=2RT_{i}T_{i+1}\frac{\Delta_{+}\log T_{i}}{\Delta_{+}T_{i}}=2R\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}

and $c_{v}(\gamma-1)=c_{v}(\frac{c_{p}}{c_{v}}-1)=R$ the last term of (131) can be rewritten as,

	$\displaystyle\frac{D_{+}^{x}H^{\prime}(T_{i})}{c_{v}(\gamma-1)2\hat{\beta}_{i+1/2jk}}{\bf f}^{d,1}_{i+1/2jk}=\frac{D_{+}^{x}H^{\prime}(T_{i})}{2R\hat{\beta}_{i+1/2jk}}{\bf f}^{d,1}_{i+1/2jk}$	$\displaystyle=$
	$\displaystyle(D_{+}^{x}H^{\prime}(T_{i}))\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}{\bf f}^{d,1}_{i+1/2jk}$	.

This allows us to rewrite the error term in (131), i.e., the last sum (denoted $E$ ). We make the following manipulations.

	$\displaystyle E\leq c_{v}\left\|\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}D_{+}^{x}H^{\prime}(T_{ijk})-D_{+}^{x}(H^{\prime}_{ijk}T_{ijk}-H_{ijk})\right\|\|{\bf f}^{d,1}_{i+1/2jk}\|$	$\displaystyle=$
	$\displaystyle c_{v}\left\|\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}D_{+}^{x}H^{\prime}_{ijk}-\overline{H^{\prime}}_{i+1/2jk}D_{+}^{x}T_{ijk}-\bar{T}_{i+1/2jk}D_{+}^{x}H^{\prime}_{ijk}+D_{+}^{x}H_{ijk}\right\|\|{\bf f}^{d,1}_{i+1/2jk}\|$	$\displaystyle=$
	$\displaystyle c_{v}\left\|\left(\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}-\bar{T}_{i+1/2jk}\right)D_{+}^{x}H^{\prime}_{ijk}-(\overline{H^{\prime}}_{i+1/2jk}-D^{T}H_{i+1/2jk})D_{+}^{x}T_{ijk}\right\|\|{\bf f}^{d,1}_{i+1/2jk}\|$	,

where $D^{T}H_{i+1/2jk}=\frac{H_{i+1jk}-H_{ijk}}{T_{i+1jk}-T_{ijk}}$ denotes the finite difference w.r.t. temperature. By the mean value theorem, we have $D^{T}H_{i+1/2jk}=H^{\prime}(\theta)$ for some $\theta\in[T_{i},T_{i+1}]$ and $|H^{\prime}(T)|\leq 1$ for all $T\geq 0$ . Furthermore,

\displaystyle(\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}-\bar{T}_{i+1/2jk})D_{+}^{x}H^{\prime}_{ijk}

\displaystyle=\frac{1}{\Delta_{+}x_{i-1/2}}(\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}-\bar{T}_{i+1/2jk})\Delta_{+}^{x}H^{\prime}_{ijk},

and $|\frac{1}{\Delta_{+}x_{i-1/2}}(\frac{\check{T}_{i+1/2}^{2}}{\hat{T}_{i+1/2}}-\bar{T}_{i+1/2jk})|\lesssim|D_{+}^{x}T_{ijk}|$ . Hence, we conclude that all temperature dependent factors appearing in $E$ scale as $H^{\prime}|D_{+}^{x}T_{ijk}|$ where $H^{\prime}$ is bounded. Hence,

\displaystyle E\lesssim|D_{+}^{x}T_{ijk}||{\bf f}^{d,1}_{i+1/2jk}|,

and to control the last term in (131), we need $\nabla_{+}T_{ijk}\in L^{2}_{2}$ and the fluxes ${\bf f}^{d,1},{\bf g}^{d,1},{\bf h}^{d,1}$ all in $L^{2}_{2}$ . The temperature bound is obtained from (107). The diffusive flux (in the x-direction) is given by

\displaystyle{\bf f}^{d,1}_{i+1/2jk}=(\frac{\mu_{0}}{\hat{\rho}_{i+1/2jk}}+\mu_{1}\bar{\rho}_{i+1/2jk}+(\Delta_{+}^{x}x_{ijk})\lambda_{i+1/2jk})D_{+}^{x}\rho_{ijk}.

The first term is, $\mu_{0}D_{+}^{x}\log\rho_{ijk}$ which bounded in $L^{2}_{2}$ by (103). The second is in $L^{2}_{2}$ by (118). The third is proportional to $(R^{*}|\bar{u}|+\frac{\Delta u}{4})\Delta_{+}^{x}\rho_{ijk}$ , where $R^{*}=\max(\frac{1}{2},\Delta\log\rho)$ . Hence, we must control terms that are majorised by: $|u|\rho$ , $\log\rho|u|\rho$ , and $\rho\Delta u$ . By noting that $\log\rho$ is controlled in any $L^{p}_{p}$ space by $L^{2}_{2}$ estimates of $\rho$ and $\rho^{-1}$ (c.f. (52)), the estimates (111), (123) and (114) provide the necessary control. (As usual, the other two directions are handled analogously.)

We have now handled the last term of (131) and we turn to the second last,

\displaystyle-\sum_{ijk=0}^{N-1,N,N}V_{ijk}(D_{+}^{x}H^{\prime}(T_{ijk}))\left(\tilde{\nu}_{i+1/2jk}\frac{1}{\gamma-1}\left(\frac{\bar{\rho}_{i+1/2jk}}{2}D_{+}^{x}\frac{1}{\beta_{ijk}}\right)\right).

(132)

Using the mean-value theorem again, we have

\displaystyle D_{+}^{x}H^{\prime}(T_{ijk})=H^{\prime\prime}(\theta_{i+1/2jk})D_{+}^{x}T_{ijk},\quad\theta_{i+1/2jk}\in[T_{ijk},T_{i+1jk}].

Furthermore, $D_{+}^{x}\frac{1}{\beta_{ijk}}=2RD_{+}^{x}T_{ijk}$ such that (132) is a negative definite term. (Hence, it makes a negative contribution to the right-hand side of (129), which is what we wish.)

We now have all estimates, to repeat the argument made for the continuous a priori estimate. Hence, we integrate (129) in time and arrive at the desired bound

\displaystyle\mathcal{C}\geq-\int_{0}^{\mathcal{T}}\sum_{ijk=0}^{N-1,N,N}V_{ijk}\kappa_{r}D_{+}^{x}H^{\prime}_{i+1/2jk}D_{+}^{x}T^{4}_{ijk}\,dt.

(133)

Our goal is to mimic the continuous estimate. That is, we wish to demonstrate that

\displaystyle D_{+}T_{ijk}^{2-\omega}\in L^{2}_{2}.

follows from (133).

Dropping the recurring $jk$ indices, we use

\displaystyle D_{+}T_{i}^{4}=(T_{i+1}^{3}+T_{i}^{3}+T_{i}^{2}T_{i+1}+T_{i+1}^{2}T_{i})D_{+}T_{i},

to rewrite

\displaystyle-D_{+}H^{\prime}_{i}D_{i}T_{i}^{4}=-\frac{\Delta_{+}H^{\prime}_{i}}{\Delta T_{i}}(T_{i+1}^{3}+T_{i}^{3}+T_{i}^{2}T_{i+1}+T_{i+1}^{2}T_{i})(D_{+}T_{i})^{2}.

(134)

By the mean-value theorem, we have

\displaystyle-\frac{\Delta_{+}H^{\prime}_{i}}{\Delta T_{i}}=\omega(1-\omega)(1+\theta_{i+1/2})^{-\omega-1}>0,\quad\theta_{i+1/2}\in[T_{i},T_{i+1}].

We recast (134) as

\displaystyle-D_{+}H^{\prime}_{i}D_{i}T_{i}^{4}=\omega(1-\omega)\frac{(T_{i+1}^{3}+T_{i}^{3}+T_{i}^{2}T_{i+1}+T_{i+1}^{2}T_{i})}{(1+\theta_{i+1/2})^{\omega+1}}(D_{+}T_{i})^{2}.

(135)

Without restriction, assume that $T_{i}\geq T_{i+1}\geq 0$ . Then

\displaystyle(1+\theta_{i+1/2})^{\omega+1}\leq(1+T_{i})^{\omega+1}.

Hence,

\displaystyle\frac{T_{i}^{3}}{(1+\theta_{i+1/2})^{\omega+1}}\geq\frac{T_{i}^{3}}{(1+T_{i})^{\omega+1}}.

Noting that all terms in (135) are positive, we can extract the bound

\displaystyle T_{i}^{2-\omega}(D_{+}T_{i})^{2}=(T_{i}^{1-\omega/2}D_{+}T_{i})^{2}.

By the mean-value theorem,

\displaystyle D_{+}T_{i}^{2-\omega/2}=(2-\omega/2)\theta^{1-\omega/2}D_{+}T_{i},

(136)

for some $\theta\in[T_{i},T_{i+1}]$ . Since $T_{i}\geq\theta$ , $T_{i}^{1-\omega/2}D_{+}T_{i}$ bounds $D_{+}T_{i}^{2-\omega/2}$ . Finally, we note that the choice $T_{i+1}<T_{i}$ was arbitrary and we can reverse the argument. Hence, we have

\displaystyle D_{+}T_{i}^{2-\omega/2}\in L^{2}_{2}

and, as in the continuous case, Sobolev embeddings and an interpolation inequality ensures that

\displaystyle T_{i}\in L^{4+\delta}_{4+\delta},

(137)

for some $\delta>0$ .

6.7 Improved estimates

The derivations of the improved estimates are the same as in Section 6.7. We summarise the ones that are necessary to prove convergence:

	$\displaystyle p_{ijk}\in L^{2+\delta}_{2+\delta},\quad\delta>0,$		(138)
	$\displaystyle{\bf v}_{ijk}\in L^{\infty}_{1},\quad{\bf v}\in L^{10/3}_{2}\cap L^{10/4}_{3},$		(139)
	$\displaystyle\quad\rho_{ijk}\|{\bf v}_{ijk}\|\in L^{2+\delta}_{2+\delta},\quad\delta>0,$		(140)
	$\displaystyle(\rho({\bf v}_{l}){\bf v}_{m})_{ijk}\in L^{1+\delta}_{1+\delta},\quad\delta>0,\quad l,m=\{1,2,3\},$		(141)
	$\displaystyle(\rho\|{\bf v}\|^{2}{\bf v})_{ijk}\in L^{1+\delta}_{1+\delta},\quad\delta>0.$		(142)

We also add the following estimate that is a consequence of the bounds (114) and (123):

\displaystyle\log\rho\in L^{p}_{p},\quad\textrm{for any}\quad 0<p<\infty.

(143)

7 Convergence to a weak solution

We begin by introducing some notation. First, we denote discrete grid functions with a subscript $h$ , ie., $\rho_{h}$ is a field whose $ijk$ th component is $(\rho_{ijk})$ . A solution obtained from the numerical scheme is hence denoted ${\bf u}(t)_{h}=(\rho(t)_{h},{\bf m}(t)_{h},E(t)_{h})$ , where $({\bf m}_{h})_{ijk}=\rho_{ijk}(u_{ijk},v_{ijk},w_{ijk})$ . We will use the convention that the fields are multiplied component-wise such that e.g. $(m_{1})_{h}=\rho_{h}u_{h}$ and ${\bf m}_{h}=\rho_{h}{\bf v}_{h}$ . In this notation, $E_{h}=\frac{1}{2}\rho_{h}|{\bf v}_{h}|^{2}+\frac{p_{h}}{\gamma-1}$ with components $(E_{h})_{ijk}=\frac{1}{2}\rho_{ijk}|{\bf v}_{ijk}|^{2}+\frac{p_{ijk}}{\gamma-1}$ . Furthermore, differences will be collected into fields using the same the short-hand notation, e.g. $D_{+}^{x}\rho_{h}$ and $\lambda_{h}D_{+}^{x}\rho_{h}$ . Finally, it is often necessary to look at the precise structure of terms and then we may say that e.g. $a_{i+1/2jk}b_{ijk}$ is bounded in some space, by that we mean that the field $a_{h}b_{h}$ is bounded in said space.

7.1 Solvability of the semi-discrete scheme

The first step towards demonstrating the existence of a weak solution, is to establish that the semi-discrete system (66) can be solved in time. Given the strong a priori estimates and that $\rho,T>0$ (as have been demonstrated above), we can apply Carathéodory’s existence theorem, to conclude that there exists a unique and continuous solution to the semi-discrete system (66), for $0\leq t\leq\mathcal{T}$ and for any $h>0$ . Hence, we can generate a sequence of triplets $\{\rho_{h},{\bf v}_{h},T_{h}\}$ .

The next step is to show that this triplet converges (up to a subsequence) to a weak solution in the sense of Def. 1 as $h\rightarrow 0$ . To that end, we proceed as in the preliminary discussion on compactness in Section 3.7. A difference in the semi-discrete case, is that we must also show that the approximate numerical fluxes are consistent with the mathematical fluxes, i.e., that the error terms vanish. (In the subsequent analysis, we will tacitly draw subsequences where necessary.)

7.2 Convergence a.e.

First, we apply Aubin-Lions lemma (Lemma 2) to the continuity equation. The arguments are initially the same as in Section 3.7, partial sums replacing partial integration. Hence, we begin to verify that the numerical fluxes of the continuity equation reside in $L^{1}_{1}$ .

The total flux in the x-direction is given by (68) and (74):

$\displaystyle{\bf f}^{1}_{i+1/2jk}$	$\displaystyle=\frac{(\rho u)_{i+1jk}+(\rho u)_{ijk}}{2}-\tilde{\nu}_{i+1/2jk}D_{+}^{x}\rho_{ijk}$	(144)
$\displaystyle\tilde{\nu}_{i+1/2jk}$	$\displaystyle=\frac{\mu_{0}}{\hat{\rho}_{i+1/2jk}}+\mu_{1}\frac{\rho_{i+1jk}+\rho_{ijk}}{2}+(\Delta_{+}^{x}x)_{i+1/2}\lambda_{i+1/2jk}$
$\displaystyle\lambda_{i+1/2jk}$	$\displaystyle\sim\|\Delta_{+}^{x}u_{ijk}\|+\|\bar{u}_{i+1/2jk}\|\|\Delta_{+}^{x}\log\rho_{ijk}\|.$

The boundedness of the first two flux terms follows (140).

Next, we consider the artificial diffusion term. $\Delta^{x}_{+}u_{h}$ is bounded (and vanishes) by (111) in $L^{2}_{2}$ . $|\Delta_{+}^{x}\log\rho_{h}||\bar{u}_{h}|$ is bounded in (at least) $L^{2}_{2}$ by (139) and by $\log\rho_{h}$ being bounded in an arbitrary $L^{p}_{p}$ space (see (143)). Using the bound (103) on the density gradient, we have

	$\displaystyle\lim_{h\rightarrow 0}\int_{0}^{\mathcal{T}}\\|(\Delta_{+}^{x}x_{i+1/2})\lambda_{i+1/2jk}D_{+}^{x}\rho_{ijk}\\|_{1}\,dt$	$\displaystyle\leq$
	$\displaystyle\lim_{h\rightarrow 0}h\int_{0}^{\mathcal{T}}(\\|\lambda_{i+1/2jk}\\|_{2}^{2}+\\|D_{+}^{x}\rho_{ijk}\\|^{2}_{2})\,dt$	$\displaystyle\rightarrow 0.$

Likewise, $\frac{1}{\hat{\rho}_{i+1/2jk}}D_{+}^{x}\rho_{ijk}=D_{+}^{x}\log\rho_{ijk}$ , which is bounded in $L^{2}_{2}$ by e.g. (103), and $\left(\mu_{1}\frac{\rho_{i+1}+\rho_{i}}{2}\right)D_{+}^{x}\rho_{ijk}$ by (118). Hence, we may conclude, as in Section 3.7 using Aubin-Lions Lemma, that $\rho_{ijk}$ converges a.e. in the limit $h\rightarrow 0$ .

Using the analogous arguments as in Section 3.7, we obtain that $\rho_{i},{\bf v}_{i},p_{i},T_{i}$ all converges a.e. Likewise, we can use the arguments of Section 3.7 to conclude that the magnitude of the gradients converge a.e. and in particular that,

\displaystyle\{\,|\Delta_{+}^{x,y,z}\rho_{ijk}|,\quad|\Delta_{+}^{x,y,z}u_{ijk}|,\quad|\Delta_{+}^{x,y,z}T_{ijk}|\,\}\rightarrow 0\quad\textrm{a.e.}.

(145)

7.3 Weak convergence of fluxes

Before considering the numerical fluxes, we briefly explain the techniques we use and why convergence is a little more delicate than in Section 3.8. First $(m_{1})_{h}=\rho_{h}u_{h}\rightarrow\rho u$ strongly in $L^{2}_{2}$ by the same arguments as laid out in Section 3.8. (We remind that $\rho_{h}u_{h}$ is the pointwise multiplication of the two fields, i.e. $\rho_{ijk}u_{ijk}$ .)

This also ensures that e.g. $\overline{\rho u}_{i+1/2jk}=\frac{\rho_{i+1jk}u_{i+1jk}+\rho_{ijk}u_{ijk}}{2}$ converges to the correct limit in $L^{2}_{2}$ . However, it does not ensure that e.g. $\bar{\rho}_{i+1/2jk}\bar{u}_{i+1/2jk}$ converges to the same limit since,

\displaystyle\bar{\rho}_{i+1/2jk}\bar{u}_{i+1/2jk}=\overline{\rho u}_{i+1/2jk}-\frac{\Delta_{+}^{x}\rho_{ijk}\Delta_{+}^{x}u_{ijk}}{4}.

Here, $\overline{\rho u}_{i+1/2jk}$ yields the desired limit and the last term is an error term that must vanish in $L^{1}_{1}$ . We can see that it does, since $D_{+}^{x}\rho_{ijk}$ reside in $L^{2}_{2}$ it follows that $\Delta_{+}^{x}\rho_{ijk}\sim hD_{+}^{x}\rho_{ijk}$ vanishes in $L^{2}_{2}$ and the same is true for $D_{+}^{x}u_{ijk}$ .

The error term can also be handled in another way: From the estimates of $\rho_{h}$ and $u_{h}$ , we can immediately conclude that $\bar{\rho}_{h}\bar{u}_{h}$ is bounded in $L^{1+\epsilon}_{1+\epsilon}$ , i.e., the sequence is equi-integrable. Since that is also true for $(\rho u)_{h}$ , we conclude that the error term has to reside in $L^{1+\epsilon}_{1+\epsilon}$ . Then by a.e. convergence of the magnitude of the gradients, using (145) it follows that that $(|\Delta_{+}^{x}\rho|)_{h}$ as well as $(|\Delta_{+}^{x}u|_{h}$ vanish a.e. (Note that the last approach does not rely on an explicit bound of the error term via gradient estimates.)

We will employ both methods to demonstrate consistency of the numerical fluxes.

Continuity equation: We need to show that ${\bf f}^{1}_{h}$ in (144) converges weakly in $L^{1+\epsilon}_{1+\epsilon}$ , $\epsilon>0$ . (As usual, the $yz$ -fluxes are handled in the same way and omitted.)

As discussed above, $\overline{\rho u}_{i+1/2jk}$ converges strongly in (at least) $L^{2}_{2}$ , by (140) and a.e. convergence of $\rho_{i}$ and $u_{i}$ . The remaining terms are handled as in Section 3.8: First, we use a.e. convergence of $\rho$ and the estimates in (114), to deduce strong convergence of $\rho^{-1}_{h}$ in $L^{2+\delta}_{2+\delta}$ for some $\delta>0$ . Then we consider the diffusive $x$ -flux:

	$\displaystyle\tilde{\nu}_{i+1/2jk}D_{+}^{x}\rho_{ijk}$	$\displaystyle=$
	$\displaystyle\left((\Delta_{+}^{x}x_{ijk})(R^{*}\|\bar{u}_{i+1/2jk}\|+\frac{\|\Delta u\|}{4})+\frac{\mu_{0}}{\hat{\rho}_{i+1/2jk}}+\mu_{1}\bar{\rho}_{i+1/2jk}\right)D_{+}^{x}\rho_{ijk}.$

These terms are proportional to: $|\Delta\log\rho||\bar{u}|\Delta\rho,|\Delta u|\Delta\rho,D_{+}\rho^{2}$ and $D_{+}\log\rho$ , which are controlled in $L^{1+\epsilon}_{1+\epsilon}$ by (123), (103), (111), (143) and (139). Weak convergence of the two diffusion terms in $L^{1+\epsilon}_{1+\epsilon}$ is immediate and the artificial diffusion terms vanish since $|\Delta\rho|\rightarrow 0$ a.e. by (145) (and the remaining factors are functions of the primary variables that converge a.e.).

Momentum equations: We consider a few typical terms and the remaining can be handled in the same way. We show that all terms are bounded in $L^{1+\epsilon}_{1+\epsilon}$ , $\epsilon>0$ , and convergent a.e. to the correct limit.

Convergence of the temporal term follows immediately from the strong convergence of $\rho_{h}{\bf v}_{h}$ in $L^{2}_{2}$ . (See (140).)

The inviscid velocity terms are all of the form:

\displaystyle\bar{u}_{i+1/2jk}\overline{(\rho v)}_{i+1/2jk}=\overline{u\rho v}_{i+1/2jk}-\frac{1}{4}\Delta_{+}^{x}u_{ijk}\Delta_{+}^{x}(\rho v)_{ijk}.

$\overline{u\rho v}_{i+1/2jk}$ is equi-integrable thanks to $\rho_{h}{\bf v}_{h}\in L^{2}_{2}$ (c.f. (140)) and ${\bf v}_{h}\in L^{5/2}_{5/2}$ . Pointwise convergence of $\rho$ and ${\bf v}$ ensure the correct limit. Noting that $\bar{u}_{i+1/2jk}\overline{(\rho v)}_{i+1/2jk}$ is bounded in the same space by the same estimates, we conclude immediately that the error term $\frac{1}{4}\Delta_{+}^{x}u_{ijk}\Delta_{+}^{x}(\rho v)_{ijk}$ is bounded in $L^{1+\delta}_{1+\delta}$ . Furthermore, by (145) both differences approach zero a.e. Hence, we conclude that these terms converge weakly in $L^{1+\epsilon}_{1+\epsilon}$ for $0<\epsilon<\delta$ .

Next, we consider the pressure term (suppressing the $jk$ indices):

	$\displaystyle\frac{1}{R}\frac{\bar{\rho}_{i+1/2}}{2\bar{\beta}_{i+1/2}}=\frac{(\rho_{i+1}+\rho_{i})}{2}\frac{2}{\frac{1}{T_{i+1}}+\frac{1}{T_{i}}}$	$\displaystyle=$
	$\displaystyle\frac{\rho_{i+1}T_{i+1}+\rho_{i}T_{i}}{2}-\rho_{i}(\frac{T_{i}}{2}-\frac{1}{\frac{1}{T_{i}}+\frac{1}{T_{i+1}}})-\rho_{i+1}(\frac{T_{i+1}}{2}-\frac{1}{\frac{1}{T_{i}}+\frac{1}{T_{i+1}}}).$

The $\rho_{i}T_{i}$ terms converge to $p_{i}$ in $L^{2}_{2}$ thanks to a.e. convergence and (138). Moreover, $\bar{\rho}_{i+1/2}$ is controlled by (123) and $\frac{1}{2\bar{\beta}_{i+1/2}}\sim\frac{T_{i}T_{i+1}}{T_{i}+T_{i+1}}\lesssim T_{i+1}+T_{i}$ is controlled by (137). Hence, $\frac{\bar{\rho}_{i+1/2}}{2\bar{\beta}_{i+1/2}}$ is bounded in (a better space than) $L^{2}_{2}$ . Therefore, we conclude that the error terms are bounded in $L^{2}_{2}$ .

Consider the first error term:

	$\displaystyle\rho_{i}(\frac{T_{i}}{2}-(\frac{T_{i+1}T_{i}}{T_{i}+T_{i+1}})=\rho_{i}T_{i}(\frac{1}{2}-(\frac{T_{i+1}}{T_{i}+T_{i+1}})$	$\displaystyle=$
	$\displaystyle\rho_{i}T_{i}(\frac{T_{i}+T_{i+1}-2T_{i+1}}{2(T_{i}+T_{i+1})})=\rho_{i}T_{i}(\frac{T_{i}-T_{i+1}}{2(T_{i}+T_{i+1})})$

In the last expression, we note that $\rho_{i},T_{i},T_{i+1}$ are a.e convergent and $\Delta T_{i}\rightarrow 0$ a.e. Hence, the error terms vanish.

Turning to the diffusive terms, we must show that

\displaystyle\nu D_{+}^{x}\rho{\bf v}\sim(\bar{\rho}+\frac{1}{\hat{\rho}})(\bar{\rho}D_{+}^{x}{\bf v}+\bar{\bf v}D_{+}^{x}\rho),

converge weakly to the correct limits. (Here, and frequently below, we have dropped the indices altogether, since the calculations only involve to neighbouring points.) That is to say that the following error terms converge to zero:

	$\displaystyle\left\|\overline{\rho^{2}}D_{+}^{x}{\bf v}-(\bar{\rho})^{2}D_{+}^{x}{\bf v}\right\|,$
	$\displaystyle\left\|\frac{1}{2}\bar{\bf v}D_{+}^{x}\rho^{2}-\bar{\bf v}\bar{\rho}D_{+}^{x}\rho\right\|,$
	$\displaystyle\left\|\frac{\bar{\rho}}{\hat{\rho}}D_{+}^{x}{\bf v}-D_{+}^{x}{\bf v}\right\|,$
	$\displaystyle\left\|\frac{1}{\hat{\rho}}\bar{\bf v}D_{+}^{x}\rho-\bar{\bf v}D_{+}^{x}\log\rho\right\|.$

In the first, both terms are bounded in $L^{1+\delta}_{1+\delta}$ by the a priori estimates. (We have $D{\bf v}\in L^{2}_{2}$ and $\rho_{i}$ in a better space than $L^{4}_{4}$ .) Hence, also the difference is bounded in the same space. Furthermore, $|\overline{\rho^{2}}-(\bar{\rho})^{2}|=\frac{1}{2}|(\Delta_{+}^{x}\rho)^{2}|$ which vanishes a.e. such that it converges strongly in a space better than $L^{4}_{4}$ . Hence, the error is a product of a weakly and a strongly convergent sequence, such that the product converges weakly to zero.

The second error term is identically zero, since $D_{+}^{x}\rho_{i}^{2}=2\bar{\rho}_{i+1/2}D_{+}^{x}\rho_{i}$ and both terms are bounded in $L^{1+\delta}_{1+\delta}$ by (139) and (118).

In the third, we have a direct bound on $D_{+}^{x}{\bf v}\in L^{2}_{2}$ . Furthermore, $1/\hat{\rho}\in L^{2+\delta}_{2+\delta}$ thanks to the $L^{\infty}_{1}$ and $L^{2}(H^{1})$ bounds of $\rho_{i}^{-1}$ , and since $\hat{\rho}_{i+1/2}$ is a monotone average bounded below and above by $\rho_{i}$ and $\rho_{i+1}$ . As before, we have $\bar{\rho}D_{+}^{x}{\bf v}\in L^{2}_{2}$ and therefore $\frac{\bar{\rho}}{\hat{\rho}}D_{+}^{x}{\bf v}\in L^{1+\delta}_{1+\delta}$ . Hence, the difference is also bounded in $L^{1+\delta}_{1+\delta}$ . Seeing that $|\frac{\bar{\rho}}{\hat{\rho}}-1|=|\frac{\bar{\rho}-\hat{\rho}}{\hat{\rho}}|\leq\frac{|\Delta_{+}^{x}\rho|}{\hat{\rho}}$ approaches zero a.e., we have the desired convergence.

The last term is identically zero, since $\frac{D_{+}^{x}\rho}{\hat{\rho}}=D_{+}^{x}\log\rho$ .

Turning to the artificial diffusion, we consider,

\displaystyle\lambda_{i+1/2jk}\Delta_{+}^{x}(\rho u)_{ijk}.

(146)

This term is bounded in $L^{1+\epsilon}_{1+\epsilon}$ since $\lambda_{h}\sim|\log\rho_{h}||u|_{h}\in L^{2+\delta}_{2+\delta}$ (by (143) and (139)) and by using (140) to bound $\rho_{h}u_{h}$ . To see that (146) vanishes a.e., we recast it as,

\displaystyle|\lambda_{i+1/2jk}\Delta_{+}^{x}(\rho u)_{ijk}|\leq\lambda_{i+1/2jk}(|\bar{\rho}_{i+1/2jk}|\Delta_{+}^{x}u_{ijk}|+|\bar{u}_{i+1/2jk}||\Delta_{+}^{x}\rho_{ijk}|).

Since all factors converge a.e., and the differences to zero a.e., the artificial diffusion term is a.e. converging to zero. Hence, it vanishes in $L^{1+\epsilon}_{1+\epsilon}$ , for some $\epsilon>0$ .

Total energy:

The equation for total energy is the most difficult one. We consider some representative terms and demonstrate convergence one by one. The remaining terms follow the same pattern.

Time derivative [ $p_{ijk}/(\gamma-1)+\frac{1}{2}\rho_{ijk}|{\bf v}_{ijk}|^{2}$ ]: The pressure part converges strongly (at least in $L^{2}_{2}$ ) thanks to (138) and a.e. convergence of $\rho_{h},T_{h}$ . The kinetic energy is bounded by (141) and converges a.e., and therefore also strongly in $L^{1+\epsilon}_{1+\epsilon}$ ( $0<\epsilon<\delta$ ).

inviscid flux: The x-flux is given by:

	$\displaystyle{\bf f}^{c,5}_{i+1/2jk}=$	$\displaystyle\frac{1}{2(\gamma-1)\hat{\beta}_{i+1/2jk}}\overline{\rho u}_{i+1/2jk}-\frac{\overline{\|{\bf v}\|^{2}}_{i+1/2jk}}{2}\overline{\rho u}_{i+1/2jk}$
		$\displaystyle+\overline{\|{\bf v}\|}^{2}_{i+1/2jk}\overline{\rho u}_{i+1/2jk}+p_{i+1/2jk}\bar{u}_{i+1/2jk}$		(147)

Consider the internal-energy flux, $IE_{1}=\frac{1}{2(\gamma-1)\hat{\beta}_{i+1/2}}\overline{\rho u}_{i+1/2}$ (dropping the recurring $jk$ indices.). Our aim is to prove consistency with $\frac{\overline{pu}_{i+1/2}}{\gamma-1}$ . That is, we need to show that

\displaystyle E_{1}=|IE_{1}-\frac{\overline{pu}_{i+1/2}}{\gamma-1}|\rightarrow 0.

Or,

\displaystyle E_{1}\leq|IE_{1}-\frac{R\bar{T}_{i+1/2}\overline{\rho u}_{i+1/2}}{\gamma-1}|+|\frac{R\bar{T}_{i+1/2}\overline{\rho u}_{i+1/2}}{\gamma-1}-\frac{\overline{pu}_{i+1/2}}{\gamma-1}|\rightarrow 0.

The available a priori estimates provide bounds on all terms (in at least $L^{1+\epsilon}_{1+\epsilon}$ ). Hence, both differences are bounded in the same space. The second one is easily handled by using (63),

\displaystyle|\frac{R\bar{T}_{i+1/2}\overline{\rho u}_{i+1/2}}{\gamma-1}-\frac{\overline{pu}_{i+1/2}}{\gamma-1}|\lesssim|\Delta_{+}^{x}T_{i}||\Delta_{+}^{x}(\rho u)_{i}|.

This approaches zero a.e since $|\Delta_{+}^{x}T_{i}|$ does and $\rho u$ converges a.e.

Turning to the first part of the error, we need an estimate of

$\displaystyle A=\bar{T}-\frac{1}{\widehat{T^{-1}}}=\bar{T}-\frac{1}{\frac{\Delta 1/T}{\Delta\log T^{-1}}}$	$\displaystyle=$
$\displaystyle\bar{T}-\frac{1}{\frac{\Delta T^{-1}}{-\Delta\log T}}=\bar{T}+\frac{\Delta\log T}{\Delta T^{-1}}$	$\displaystyle=$
$\displaystyle\bar{T}+\frac{\Delta\log T}{-\frac{\Delta T}{\check{T}^{2}}}=\bar{T}-\frac{\check{T}^{2}}{\hat{T}}$	$\displaystyle=\frac{\bar{T}\hat{T}-\check{T}^{2}}{\hat{T}},$	(148)

where again we have suppressed the indices for brevity. We introduce $T_{*}^{2}=\bar{T}\hat{T}\in[\hat{T}^{2},\bar{T}^{2}]$ , to get,

\displaystyle A=\frac{(T_{*}+\check{T})(T_{*}-\check{T})}{\hat{T}}\leq\frac{(T_{*}+\check{T})|\Delta T|}{\frac{\Delta T}{\Delta\log T}}\leq(T_{*}+\check{T})|\Delta\log T|.

(149)

Hence,

\displaystyle|IE_{1}-\frac{R\bar{T}_{i+1/2}\overline{\rho u}_{i+1/2}}{\gamma-1}|\sim(T_{*}+\check{T})_{i+1/2}|\Delta_{+}\log T_{i}|(\overline{\rho u})_{i+1/2}.

We already know that this difference is equi-integrable and we conclude that it vanishes since all factors converge a.e, and in particular $|\Delta T|\rightarrow 0$ . (The term $p_{i+1/2jk}\bar{u}_{i+1/2jk}$ in (147) can be handled in the same way.)

Next, we turn to the convective part of the inviscid flux. The kinetic-energy fluxes appearing in (147) are bounded in $L^{1+\epsilon}_{1+\epsilon}$ by the same arguments leading to (51) (seeing that both terms are bounded by products of $\rho{\bf v}$ and $|{\bf v}|^{2}$ .).

The numerical kinetic-energy flux is compared with the simple average, $\frac{1}{2}\overline{\rho|{\bf v}|^{2}u}_{i+1/2jk}$ , which is bounded (as in (51)) and converges to the correct limit in the weak formulation thanks to the a.e. convergence of $\rho$ and ${\bf v}$ .

The error term that should vanish is thus:

\displaystyle E_{2}=\frac{1}{2}\overline{\rho|{\bf v}|^{2}u}_{i+1/2jk}-\left(-\frac{\overline{|{\bf v}|^{2}}_{i+1/2jk}}{2}\overline{\rho u}_{i+1/2jk}+\overline{|{\bf v}|}^{2}_{i+1/2jk}\overline{\rho u}_{i+1/2jk}\right).

Since both the numerical flux and the target flux are bounded in $L^{1+\epsilon}_{1+\epsilon}$ , the error $E_{2}$ is bounded in the same space.

Using (63), we recast the error,

\displaystyle E_{2}=\frac{1}{2}\overline{\rho u}\overline{|{\bf v}|^{2}}+\frac{\Delta\rho u\Delta|{\bf v}|^{2}}{8}-\left(-\frac{\overline{|{\bf v}|^{2}}}{2}\overline{\rho u}+\overline{|{\bf v}|}^{2}\overline{\rho u}\right).

Noting that $\Delta|{\bf v}|^{2}=2\bar{\bf v}D{\bf v}$ and $\overline{|{\bf v}|^{2}}-\overline{|{\bf v}|}^{2}=\frac{1}{4}(|\Delta{\bf v}|)^{2}$ , the error terms take the form

	$\displaystyle\frac{\Delta\rho u\Delta\|{\bf v}\|^{2}}{8}$	$\displaystyle=\frac{\Delta(\rho u)\bar{\bf v}\Delta{\bf v}}{4},$
	$\displaystyle\overline{\rho u}(\overline{\|{\bf v}\|^{2}}-\overline{\|{\bf v}\|}^{2})$	$\displaystyle=\frac{1}{4}\overline{\rho u}(\|\Delta{\bf v}\|)^{2}.$

Both of these terms can be bounded as products of $\rho{\bf v}$ and $|{\bf v}|^{2}$ and are thus equi-integrable. We see that they vanish a.e. since they are products of primitive variables and differences.

Diffusive flux:

We begin by considering the kinetic-energy part of the flux in the x-direction:

\displaystyle\tilde{\nu}_{i+1/2jk}\left(\frac{1}{2}(D_{+}^{x}\rho_{ijk}|{\bf v}_{ijk}|^{2})+(\overline{|{\bf v}|}^{2}_{i+1/2jk}-\overline{|{\bf v}|^{2}}_{i+1/2jk})D_{+}^{x}\rho_{ijk}\right)

(150)

The last two terms two terms approximate the same quantity and we begin by showing that their sum converges to zero (in $L^{1}_{1}$ ). Since $\overline{|{\bf v}|^{2}}-\overline{|{\bf v}|}^{2}=\frac{1}{4}(|\Delta{\bf v}|)^{2}$ , we get,

\displaystyle\tilde{\nu}(\overline{|{\bf v}|}^{2}-\overline{|{\bf v}|^{2}})D_{+}^{x}\rho=\frac{\tilde{\nu}}{4}(|\Delta_{+}^{x}{\bf v}|)^{2}D_{+}^{x}\rho=h\frac{\tilde{\nu}}{4}(|D_{+}^{x}{\bf v}|)^{2}\Delta_{+}^{x}\rho

Using (112), we conclude that the term vanishes in $L^{1}_{1}$ . (Note that the $h$ -factor also makes it equi-integrable.)

Next, we turn to the part of (150) that should converge to the correct average as given in Def. 1. That is, the term

\displaystyle\tilde{\nu}_{i+1/2jk}\frac{1}{2}D_{+}^{x}\rho_{ijk}|{\bf v}_{ijk}|^{2}.

(151)

First, we handle the terms arising from the $\bar{\rho}$ dependence of $\tilde{\nu}$ . We also limit the analysis to one velocity component as the others can be handled analogously. We have (when suppressing the indices),

\displaystyle\bar{\rho}D(\rho u^{2})=\bar{\rho}(\bar{\rho}Du^{2}+\overline{u^{2}}D\rho).

(152)

The first term of (152):

\displaystyle\bar{\rho}\bar{\rho}Du^{2}=\bar{\rho}\bar{\rho}\bar{u}Du=(\bar{\rho}\bar{u})(\bar{\rho}Du)=(\overline{\rho u}-\frac{\Delta\rho\Delta u}{4})\bar{\rho}Du.

(153)

We use that $\bar{\rho}Du\in L^{2}_{2}$ and $\overline{\rho u}\in L^{2+\delta}_{2+\delta}$ . We need the error term to vanish. To this end we note that $\Delta\rho\Delta u=h(\rho_{1}Du-\rho_{2}Du)$ and $\rho_{1,2}Du$ is bounded in $L^{2}_{2}$ . Furthermore, the $h$ -factor makes (153) both equi-integrable and ensures that it vanishes in $L^{1}_{1}$ .

Furthermore, $\overline{\rho u}\bar{\rho}Du$ approximates $\overline{\rho^{2}u}Du$ up to an error term,

\displaystyle\Delta(\rho u)\Delta\rho Du=\bar{u}(\Delta\rho)^{2}Du+\bar{\rho}(\Delta\rho)(\Delta u)Du.

The last term is the same as in (153) and handled in the same way. The second last is rewritten as:

\displaystyle\bar{u}(\Delta\rho)^{2}Du=\bar{u}\Delta\rho(\Delta\rho Du)=(\Delta(\rho u)-(\Delta u)\bar{\rho})(\Delta\rho Du).

The $\Delta(\rho u)(\Delta\rho Du)$ is equi-integrable by $\bar{\rho}Du\in L^{2}_{2}$ and $\rho u\in L^{2+\delta}_{2+\delta}$ (see 140) and (110)). It vanishes a.e. since $|(\Delta\rho u)(\Delta\rho Du)|\lesssim|\rho u||\Delta\rho||Du|$ where $|Du|$ and $\rho u$ converges a.e and $|\Delta\rho|\rightarrow 0$ a.e. Finally, the $(\bar{\rho}\Delta u)(\Delta\rho Du)$ is the same as a previously handled term.

The second term of (152) is $\overline{u^{2}}\bar{\rho}D\rho=\frac{1}{2}\overline{u^{2}}D\rho^{2}=\frac{1}{2}D(\rho^{2}u^{2})-\overline{\rho^{2}}\bar{u}Du$ . (This is the numerical counterpart of (56).) In view of (56) and by using (63), the $\overline{\rho^{2}}\bar{u}Du$ -term results in an error term,

\displaystyle-(\overline{\rho^{2}}\bar{u}-\overline{\rho^{2}u})Du=\frac{\Delta\rho^{2}\Delta u}{4}Du=\frac{1}{2}(\Delta\rho\Delta u)(\bar{\rho}Du).

We have already handled a term of this form and know that it vanishes.

The above considerations show that the $\bar{\rho}$ -dependent part of $\tilde{\nu}$ in the kinetic-energy flux can be recast as

\displaystyle\frac{1}{4}D(\rho^{2}u^{2})-\frac{1}{2}\overline{\rho^{2}u}Du+\textrm{error terms},

where the error terms all vanish. The remaining part converges to the correct limit by the same arguments as in Section 3.8. In short, in the term $\frac{1}{4}D(\rho^{2}u^{2})$ we move the difference to the test function which results in a momentum term that is equi-integrable and a.e. convergent. Moreover, recasting the second term as

	$\displaystyle\overline{\rho^{2}u}Du=\frac{1}{2}\left(\rho^{2}_{1}u_{1}+\rho^{2}_{2}u_{2}\right)Du$	$\displaystyle=$
	$\displaystyle\frac{1}{2}\left((\rho_{1}u_{1})(\rho_{1}Du)+(\rho_{2}u_{2})(\rho_{2})\right)Du$	,

( $1,2$ denotes the two points appearing in the average and differences) and using the same estimates as above, shows that it is equi-integrable. To demonstrate convergence to the correct limit, we argue that $\bar{\rho}Du$ converges weakly to the correct limit as in Section 3.8. Strong convergence of momentum, completes the argument.

Second, we handle the terms arising from the $(\hat{\rho})^{-1}$ dependence of $\tilde{\nu}$ in (151). Namely,

	$\displaystyle\frac{1}{\hat{\rho}}D(\rho u^{2})=\frac{1}{\hat{\rho}}\left(\overline{u^{2}}D\rho+2\bar{\rho}\bar{u}Du\right)=\overline{u^{2}}D\log\rho+2\frac{\bar{\rho}}{\hat{\rho}}\bar{u}Du$	$\displaystyle=$
	$\displaystyle D(u^{2}\log\rho)-2\overline{\log\rho}\bar{u}Du+2\frac{\bar{\rho}}{\hat{\rho}}\bar{u}Du.$

The difference in the $D(u^{2}\log\rho)$ is moved to the test function. It is bounded in $L^{1+\epsilon}_{1+\epsilon}$ and consistent with the corresponding term in Def. 1. In the second, $\bar{u}Du$ is bounded in $L^{1+\epsilon}_{1+\epsilon}$ and weakly convergent (by the analogous argument as $\rho\nabla{\bf v}$ is weakly convergent. See end of Section 3.7.) Since $\log\rho$ can be bounded in any $L^{p}_{p}$ space, and is a.e. convergent, it is also strongly convergent in a sufficiently good space. This ensures convergence to the right limit of $\overline{u\log\rho}Du$ . The error term $|\overline{\log\rho}\bar{u}Du-\overline{u\log\rho}Du|$ vanishes, since it is majorised by $|\Delta u||\Delta\log\rho||Du|$ . This is equi-integrable thanks to $\bar{u}Du\in L^{1+\epsilon}_{1+\epsilon}$ and the $L^{p}_{p}$ estimates of $\log\rho$ . It approaches zero a.e, since all factors are a.e. convergent and $|\Delta\log\rho|\rightarrow 0$ .

Finally, the last term should converge to $\bar{u}Du$ , which is weakly convergent in $L^{1}_{1}$ . Hence, we must show that the following error term vanishes:

	$\displaystyle\|(1-\frac{\bar{\rho}}{\hat{\rho}})\bar{u}Du\|=\|\frac{\hat{\rho}-\bar{\rho}}{\hat{\rho}}\bar{u}Du\|$	$\displaystyle\leq$
	$\displaystyle\|\frac{\Delta\rho}{\hat{\rho}}\bar{u}Du\|$	$\displaystyle=\|(\Delta\log\rho)\bar{u}Du\|.$

We have already shown that such a term is equi-integrable and vanishes.

Turning to the artificial diffusion part of (151), we need that

\displaystyle h\lambda D(\rho u^{2})\sim h(|u|+|\Delta\log\rho||u|+|\Delta u|)D(\rho u^{2}),

vanishes in $L^{1}_{1}$ . We begin with

\displaystyle h|\Delta u|D(\rho u^{2})=|\Delta u|\Delta(\rho u^{2}).

(154)

Let the subscripts $1,2$ denote the two points on which the flux term is based. Then,

	$\displaystyle(u_{1}-u_{2})(\rho_{1}u_{1}^{2}-\rho_{2}u_{2}^{2})=\rho_{1}u_{1}^{3}-\rho_{2}u_{2}^{3}-\rho_{1}u_{1}^{2}u_{2}-\rho_{2}u_{2}^{2}u_{1}$	$\displaystyle=$
	$\displaystyle\rho_{1}u_{1}^{3}-\rho_{2}u_{2}^{3}-\rho_{1}u_{1}(u_{1}u_{2})-\rho_{2}u_{2}(u_{2}u_{1})$	$\displaystyle\leq$
	$\displaystyle\rho_{1}u_{1}^{3}+\rho_{2}u_{2}^{3}-\rho_{1}u_{1}(u_{1}^{2}+u_{2}^{2})-\rho_{2}u_{2}(u_{2}^{2}+u_{1}^{2}),$

shows that (154) can be majorised by $\overline{\rho u}\overline{u^{2}}$ , which can be split in the same way leading to the bound (142). Hence, it is equi-integrable. Moreover, $|\Delta u|\rightarrow 0$ a.e. ensures that (154) vanishes a.e.

The same argument applies directly to $h|u|D(\rho u^{2})$ .

This leaves us with the artificial diffusion term,

\displaystyle|\Delta\log\rho||\bar{u}|\Delta(\rho u^{2}).

(155)

As previously shown, $|\bar{u}|\Delta(\rho u^{2})\in L^{1+\epsilon}_{1+\epsilon}$ and $\Delta\log\rho$ can be bounded in a sufficiently high $L^{p}_{p}$ space to obtain equi-integrability of the error term (155). Moreover, (155) vanishes a.e. since $\Delta(\rho u^{2})=(\Delta\rho)\overline{u^{2}}+2\bar{u}\bar{\rho}\Delta u$ and the differences vanish.

Diffusive flux, internal energy [ $\tilde{\nu}_{i+1/2jk}\frac{({\mathfrak{p}}_{x})_{i+1/2jk}}{\gamma-1}$ ]:

\displaystyle\tilde{\nu}_{i+1/2jk}({\mathfrak{p}}_{x})_{i+1/2jk}

\displaystyle=\tilde{\nu}_{i+1/2jk}\left(\frac{1}{2\hat{\beta}_{i+1/2jk}}D_{+}^{x}\rho_{ijk}+\frac{\bar{\rho}_{i+1/2jk}}{2}D_{+}^{x}\frac{1}{\beta_{ijk}}\right)

(156)

The “target terms” that, up to some constants, should be approximated are,

\displaystyle\rho TD\rho,\quad\rho^{2}DT,\quad TD\log\rho,\quad DT.

We begin with the terms associated with $\tilde{\nu}\sim\bar{\rho}$ . They generate error terms of the form:

	$\displaystyle\left(\frac{\bar{\rho}}{2\hat{\beta}}-\overline{R\rho T}\right)D\rho$	,		(157)
	$\displaystyle\left(\bar{\rho}^{2}-\overline{\rho^{2}}\right)DT$	$\displaystyle=-\frac{1}{4}\left(\rho_{1}-\rho_{2}\right)^{2}DT.$		(158)

The second term, (158), is bounded by $\rho\in L^{4}_{4}$ and $DT\in L^{2}_{2}$ . Furthermore, $DT$ converges weakly and $(\Delta\rho)^{2}$ strongly since $|\Delta\rho|\rightarrow 0$ a.e. Hence, the error term (158) vanishes in $L^{1}_{1}$ . The term $\overline{\rho^{2}}DT$ converges to the correct limit thanks to the strong convergence of $\rho$ and weak convergence of $DT$ .

Turning to (157), we make the following manipulations,

\displaystyle\left(\frac{\bar{\rho}}{2\hat{\beta}}-\overline{R\rho T}\right)=\left(\frac{\bar{\rho}}{2\hat{\beta}}-R\bar{\rho}\bar{T}+R\frac{\Delta\rho\Delta T}{4}\right)=\left(\bar{\rho}\left(\frac{1}{2\hat{\beta}}-R\bar{T}\right)+R\frac{\Delta\rho\Delta T}{4}\right)

(159)

Inserting (159) in (157), we see that we need to control $\Delta\rho\Delta TD\rho$ , or equivalently, $(\Delta\rho)^{2}DT$ . This is the same as term as (158), which we have already handled. This leaves us with $\bar{\rho}\left(\frac{1}{2\hat{\beta}}-R\bar{T}\right)D\rho$ (from (159) and (157)). To handle this term, we use (148) and (149), to obtain

\displaystyle|\frac{1}{2R\hat{\beta}}-\bar{T}|\lesssim\bar{T}|\Delta\log T|.

Hence, we must control $\bar{\rho}(D\rho)\bar{T}|\Delta\log T|$ . $\bar{\rho}D\rho=\frac{1}{2}D\rho^{2}\in L^{2}_{2}$ which converges weakly. We need strong convergence (to zero) of $|\bar{T}\Delta\log T|$ in a better space than $L^{2}_{2}$ . We note that it converges a.e. to zero and only need the bound. This follows by interpolating the bounds $T^{3/2}D\log T\in L^{2}_{2}$ (see (108)) and $D\log T\in L^{2}_{2}$ (see (104)) such that $TD\log T\in L^{2}_{2}$ . This implies that $|\bar{T}\Delta\log T|\in L^{2+\delta}_{2+\delta}$ , $\delta>0$ . Finally, the target entity $\overline{R\rho T}D\rho$ converges to the correct limit thanks to the strong convergence of $p\in L^{2+\delta}_{2+\delta}$ (following from (138) and a.e. convergence) and the weak convergence of $D\rho\in L^{2}_{2}$ .

Next, we consider error terms of (156) for $\tilde{\nu}\sim 1/\hat{\rho}$ :

	$\displaystyle\left\|\left(\frac{1}{2\hat{\beta}}-\overline{T}\right)\frac{D\rho}{\hat{\rho}}\right\|=\left\|\left(\frac{1}{2\hat{\beta}}-\overline{T}\right)D\log\rho\right\|$	$\displaystyle\lesssim\bar{T}\|\Delta\log T\|\|D\log\rho\|$
	$\displaystyle\left\|\left(1-\frac{\bar{\rho}}{\hat{\rho}}\right)DT\right\|\leq\left\|\left(\frac{\|\Delta\rho\|}{\hat{\rho}}\right)DT\right\|$	$\displaystyle\leq\|\Delta\log\rho\|\|DT\|\sim\|D\log\rho\|\|\Delta T\|.$

The first is handled in the same way as above, only now we are using that $D\log\rho\in L^{2}_{2}$ . The second is handled by the strong bound on $T$ and a.e. convergence, implying strong convergence of $|\Delta T|$ (to zero), and $D\log\rho\in L^{2}_{2}$ . Here, the target is $\bar{T}D\log\rho$ and it is easy to verify convergence to the correct limit.

The error terms of (156) for $\tilde{\nu}\sim h((\frac{1}{2}+|\Delta\log\rho|)|\bar{u}|+\Delta u)$ (artificial diffusion):

	$\displaystyle h((\frac{1}{2}+\|\Delta\log\rho\|)\|\bar{u}\|+\|\Delta u\|)\left(\frac{D\rho}{\hat{\beta}}+\bar{\rho}DT\right)$	$\displaystyle=$
	$\displaystyle((\frac{1}{2}+\|\Delta\log\rho\|)\|\bar{u}\|+\|\Delta u\|)\left(\frac{\Delta\rho}{\hat{\beta}}+\bar{\rho}\Delta T\right)$

The strong bounds imply that

	$\displaystyle\left(\frac{\Delta\rho}{\hat{\beta}}+\bar{\rho}\Delta T\right)$	$\displaystyle\in L^{2+\epsilon}_{2+\epsilon}\quad\textrm{and}$
	$\displaystyle(\frac{1}{2}+\|\Delta\log\rho\|)\|\bar{u}\|+\|\Delta u\|$	$\displaystyle\in L^{2}_{2}.$

(Note that $\Delta T\in L^{4}_{4}$ and $1/\hat{\beta}\in L^{4}_{4}$ .) This gives equi-integrability. The a.e. vanishing differences ensures that the error terms disappear.

Radiation diffusion [ $\kappa_{r}D_{-}^{x}D_{+}^{x}T^{4}_{ijk}$ ]: This term is straightforwardly handled by moving the differences to the test function. This results in a $T^{4}_{h}$ term, which is strongly convergent thanks to (137) and the a.e convergence of temperature.

Entropy inequality: The global entropy diffusivity condition in Def. 1 follows directly from (101), the strong bounds on the solution and the positive definiteness of the diffusive terms appearing in the left-hand side of (101).

8 Concluding remarks

The main result of this paper (Theorem 1) is the existence of weak entropy solutions to the alternative Navier-Stokes system (4). These weak solutions ensure positivity of temperature and density, except possibly on a set of Lebesgue measure zero. To the best of our knowledge, such results are not available for the standard Navier-Stokes-Fourier system.

The weak solutions are obtained as the limits of solutions to a finite volume scheme. We make some further remarks on the scheme:

•

As presented, the scheme is formally first-order accurate for smooth solutions (on regular grids). This is due to the first-order artificial-diffusion coefficent, namely the “ $1/2$ ” appearing in (72), which results in an “upwind”-type diffusion. As this is only needed in the entropy estimate (see (99)), we can replace $\lambda$ in (72) by

	$\displaystyle\lambda_{i+1/2jk}$	$\displaystyle=\|\bar{u}_{i+1/2jk}\|R^{\#}_{i+1/2jk}+\frac{\|\Delta_{+}^{x}u_{ijk}\|}{4}$
	$\displaystyle R^{\#}_{i+1/2jk}$	$\displaystyle=\max\left(\left\|\frac{\bar{\rho}_{i+1/2jk}-\hat{\rho}_{i+1/2jk}}{\Delta_{+}^{x}\rho_{ijk}}\right\|,\|\Delta_{+}^{x}\log\rho_{ijk}\|\right)$

The first entry in the maximum is $\mathcal{O}(|\Delta_{+}^{x}\rho_{ijk}|)$ . (Formally, this implies second-order accuracy.) Since $R^{\#}$ is dominated by $R^{*}$ (c.f. (99)), and converges to zero a.e., it does not affect the convergence proof.

However, in order to prevent that the finite arithmetic pollutes the numerical solutions, the first entry of $R^{\#}$ requires a well-conditioned numerical approximation. This should be obtainable in the same way that the approximation of the log mean was derived (see Appendix B of [IR09]), but we leave this as future work.

•

Unfortunately, the generalisation to high-order accuracy (three or higher for smooth solutions), discovered in [FC13] is effectively ruled out, since it requires the use of entropy-conservative fluxes, which in turn appears incompatible with the current approach.
•

The scheme is trivially generelisable to the unstructured finite-volume framework developed by Eymard et al. ([EGH97]) for Voronoi control volumes and two-point fluxes. The convergence proof holds in this case as well, thanks to the summation-by-parts property of such schemes.
•

The extension of this theory to allow some form of far-field boundary condition to close the domain for external flows, appears to be non-trivial. The ideas used in [Svä21] may be a possible route to this end.

Future work, also includes investigating if (4), like many other fluid systems, has some weak-strong uniqueness property.

9 Acknowledgements

I am very grateful to Prof. Josef Málek whose help has enabled me to write this article.

References

[Bre05a] H. Brenner. Kinematics of volume transport. Physica A, 349:11, 2005.
[Bre05b] H. Brenner. Navier-Stokes revisited. Physica A, 349:60, 2005.
[Cha13] P. Chandrashekar. Kinetic energy preserving and entropy stable finite volume schemes for compressible Euler and Navier-Stokes equations. Comm. Comp. Phys., 14(5):1252–1256, 2013.
[DS21] Vít Dolejší and Magnus Svärd. Numerical study of two models for viscous compressible fluid flows. J. Comp. Phys, 427:110068, 2021.
[EGH97] R. Eymard, T. Gallouët, and R. Herbin. Finite volume methods. In P.G. Ciarlet and J.L. Lions, editors, Handbook of Numerical Analysis, volume 7, pages 713–1020. Elsevier, 1997.
[Eva10] L. C. Evans. Partial Differential Equations, 2nd ed. American Mathematical Society, 2010.
[FC13] T. C. Fisher and M. H. Carpenter. High-order entropy stable finite difference schemes for nonlinear conservation laws: Finite domains. J. Comput. Phys., 252:518–557, 2013.
[FN17] E. Feireisl and A. Novotný. Singular limits in thermodynamics of viscous fluids, 2nd Ed. Birkhäuser, 2017.
[FV10] Eduard Feireisl and Alexis Vasseur. New Perspectives in Fluid Dynamics: Mathematical Analysis of a Model Proposed by Howard Brenner, pages 153–179. Birkhäuser Basel, Basel, 2010.
[GSH21] G.J. Gassner, M. Svärd, and F.J. Hindenlang. Stability issues of entropy-stable and/or split-form high-order schemes. to appear in Journal of Scientific Computing, 2021.
[IR09] F. Ismail and P.L. Roe. Affordable, entropy-consistent euler flux functions II: Entropy production at shocks. J. Comp. Phys., pages 5410–5436, 2009.
[NFAE03] J. Nordström, K. Forsberg, C. Adamsson, and P. Eliasson. Finite volume methods, unstructured meshes and strict stability for hyperbolic problems. Applied Numerical Mathematics, 45(4):453–473, June 2003.
[RDO⁺19] M. H. L. Reddy, S. K. Dadzie, R. Ocone, M. K. Borg, and J. M. Reese. Recasting Navier-Stokes equations. Journ. Phys. Comm., 3(10):1–25, 2019.
[SDP21] M. Sayyari, L. Dalcin, and M. Parsani. Development and analysis of entropy stable no-slip wall boundary conditions for the Eulerian model for viscous and heat conducting compressible flows. SN PDEA, 2, 2021.
[Sim86] J. Simon. Compact sets in the space $L^{p}(O,T;B)$ . Annali di Matematica Pura ed Applicata, 146:65–96, 1986.
[SN14] M. Svärd and J. Nordström. Review of summation-by-parts schemes for initial-boundary-value problems. J. Comput. Phys., 268:17–38, 2014.
[Sto45] G. G. Stokes. On the theories of internal friction of fluids in motion, and of the equilibrium and motion of elastic solids. Transactions of the Cambridge Philosophical Society, 8:287–305, 1845.
[Svä18] M. Svärd. A new Eulerian model for viscous and heat conducting compressible flows. Physica A, 506:350–375, 2018.
[Svä21] M. Svärd. Entropy stable boundary conditions for the Euler equations. J. Comput. Phys., 426:1–17, 2021.

APPENDIX

I Prerequisites

An elementary identity: $\|u^{r}\|_{p}^{q}=\|u\|_{rp}^{rq}$ .

A list of standard inequalities:

•

Cauchy-Schwarz inequality: If $u,v\in L^{2}$ then $<u,v>\leq\|u\|_{2}\|v\|_{2}$ .
•

Hölder’s inequality: $\|uv\|_{1}\leq\|u\|_{p}\|v\|_{q}$ where $(p,q)\in[1,\infty)$ are Hölder conjugates, $1/p+1/q=1$ .
•

A generalisation of Hölder’s inequality (see [FN17]): $\|uv\|_{r}\leq\|u\|_{p}\|v\|_{q}$ where $(p,q,r)\in[1,\infty)$ are Hölder conjugates, $1/p+1/q=1/r$ .
•

Minkowski (triangle) inequality: $\|u+v\|_{p}\leq\|u\|_{p}+\|v\|_{p}$ for $p\in[1,\infty)$ .
•

Young’s inequality: For $a,b\in\mathcal{R}$ , $a\cdot b\leq\frac{a^{p}}{p}+\frac{b^{q}}{q}$ , where $(p,q)$ is the Hölder conjugate.

•

Young’s inequality is often combined with Hölder in the following way:

\displaystyle\|uv\|_{1}\leq\epsilon\|u\|_{p}\epsilon^{-1}\|v\|_{q}\leq\frac{(\epsilon\|u\|)^{p}}{p}+\frac{(\epsilon^{-1}\|v\|)^{q}}{q}.

(160)

•

Ladyzhenskaya’s inequality in 3-D: $\|u\|_{4}\leq\mathcal{C}\|u\|^{1/4}_{2}\|\nabla u\|_{2}^{3/4}$ .
•

Nash’ inequality in $\mathcal{R}^{n}$ : $\|u\|_{2}^{1+2/n}\leq\mathcal{C}\|u\|^{2/n}_{1}\|Du\|_{2}$ .
•

Sobolev embedding (in 3-D): $\|u\|_{6}\leq\mathcal{C}(\|u\|_{2}+\|\nabla u\|_{2})$ .

•

Riesz-Thorin (interpolation) inequality:

	$\displaystyle\\|u\\|_{r}\leq\\|u\\|_{p}^{1-\theta}\\|u\\|^{\theta}_{q}$
	$\displaystyle\frac{1}{r}=\frac{1-\theta}{p}+\frac{\theta}{q},\quad 0<\theta<.1$

An important special case of this inequality is the “10/3”-rule,

\displaystyle\|u\|^{10/3}_{10/3}\leq\|u\|^{4/3}_{2}\|u\|^{2}_{6}

(161)

which is often used in combination with Sobolev embedding. Some other special cases:

$\displaystyle\\|u\\|^{20/3}_{20/3}$	$\displaystyle\leq\\|u\\|^{2/3}_{1}\\|u\\|^{6}_{18}$	(162)
$\displaystyle\\|u\\|^{3}_{3}$	$\displaystyle\leq\\|u\\|^{3/5}_{1}\\|u\\|^{12/5}_{6}$	(163)
$\displaystyle\\|u\\|^{6}_{6}$	$\displaystyle\leq\\|u\\|^{2/3}_{3}\\|u\\|^{4}_{12}$	(164)

•

A special case of the Gagliardo-Nirenberg inequality:

$\displaystyle\|u\|_{3}\leq\mathcal{C}\|u\|_{2}^{1/2}\|\nabla u\|_{2}^{1/2}$

A generalised Poincare inequality:

Theorem 2.

Let $1\leq p\leq\infty$ , $0<\Gamma<\infty$ , $V_{0}>0$ and let $\Omega\subset\mathcal{R}^{N}$ be a bounded Lipschitz domain. Then there exists a positive constant $c=(p,\Gamma,V_{0})$ such that,

\displaystyle\|v\|_{W^{1,p}}(\Omega)\leq c\left[\|\nabla v\|_{L^{p}(\Omega;R^{N})}+(\int_{V}|v|^{\Gamma}d{\bf x})^{\frac{1}{\Gamma}}\right]

(165)

for any measurable $V\subset\Omega$ , $|V|>V_{0}$ and any $v\in W^{1,p}(\Omega)$ .

Proof.

This is Theorem 11.20 in [FN17], where it is also proved. ∎

Aubin-Lions Lemma (see e.g. [Sim86]):

Lemma 2.

For Banach spaces $X\subset B\subset Y$ where the embedding of $X$ in $B$ is compact and $B$ in $Y$ is continuous. Let $U=\{u|u\in L^{p}(0,\mathcal{T};X)\,\textrm{and}\,u_{t}\in L^{q}(0,\mathcal{T};Y)\}$ where $1\leq p,q<\infty$ . Then $U$ is compactly embedded in $L^{p}(0,\mathcal{T};B)$ .

The following Lemma for weakly and strongly converging sequences is a standard result.

Lemma 3.

Let $1<p<\infty$ and $1/p+1/q=1$ and $\Omega\subset\mathcal{R}^{N}$ . If

	$\displaystyle u_{n}\rightharpoonup u\quad\textrm{in}\quad L^{p}(\Omega),$
	$\displaystyle v_{n}\rightarrow v\quad\textrm{in}\quad L^{q}(\Omega),$

then

\displaystyle u_{n}v_{n}\rightharpoonup uv\quad\textrm{in}\quad L^{1}(\Omega),

The following Lemma, that we state without proof, is a consequence of Vitali’s Convergence Theorem:

Lemma 4.

Let $\Omega\subset\mathcal{R}^{N}$ be bounded. If $u_{n}\rightarrow u$ a.e. in $\Omega$ and $u_{n}$ is bounded in $L^{p}(\Omega)$ , $p>1$ . Then, $u_{n}\rightarrow u$ in $L^{r}$ for all $1\leq r<p$ .

The next theorem is found in [Eva10] (Appendix D, Theorem 3).

Theorem 3.

Let $u_{n}\in L^{p}(\Omega)$ , $1<p<\infty$ be a uniformly bounded sequence. Then, there is a subsequence (still denoted $\{u_{n}\}$ ) and a function $u\in L^{p}(\Omega)$ , such that $u_{n}\rightharpoonup u$ in $L^{p}(\Omega)$ .

	$\displaystyle\\|\rho\|{\bf v}\|\|{\bf v}\|^{2}\\|_{1+\delta}^{1+\delta}=\\|(\rho\|{\bf v}\|\|{\bf v}\|^{2})^{1+\delta}\\|$	$\displaystyle=$
	$\displaystyle\\|((\sqrt{\rho})^{1-\epsilon}(\sqrt{\rho})^{1+\epsilon}\|{\bf v}\|^{1+\epsilon}\|{\bf v}\|^{2-\epsilon})^{1+\delta}\\|$	$\displaystyle=$
	$\displaystyle\\|(\sqrt{\rho})^{(1+\delta)(1-\epsilon)}(\sqrt{\rho}\|{\bf v}\|)^{(1+\delta)(1+\epsilon)}\|{\bf v}\|^{(2-\epsilon)(1+\delta)}\\|$	$\displaystyle\leq$
	$\displaystyle\\|(\sqrt{\rho})^{(1+\delta)(1-\epsilon)}(\sqrt{\rho}\|{\bf v}\|)^{(1+\delta)(1+\epsilon)}\\|_{3/2}\\|\|{\bf v}\|^{(2-\epsilon)(1+\delta)}\\|_{3}$	.

	$\displaystyle\\|D^{x}_{+}a\\|^{2}=\sum_{ijk=0}^{N-1,N,N}V_{ijk}(D_{+}^{x}a_{ijk})^{2},$
	$\displaystyle\\|D^{y}_{+}a\\|^{2}=\sum_{ijk=0}^{N,N-1,N}V_{ijk}(D_{+}^{y}a_{ijk})^{2},$		(65)
	$\displaystyle\\|D^{z}_{+}a\\|^{2}=\sum_{ijk=0}^{N,N,N-1}V_{ijk}(D_{+}^{z}a_{ijk})^{2}.$

$\displaystyle-\frac{1}{2}\|{\bf v}\|^{2}_{ijk}\frac{d\rho_{j}}{dt}$	$\displaystyle=$
$\displaystyle\frac{1}{2}\|{\bf v}\|^{2}_{ijk}\left(D_{-}^{x}{\bf f}^{1}_{i+1/2jk}+D_{-}^{y}{\bf g}^{1}_{ij+1/2k}+D_{-}^{z}{\bf h}^{1}_{ijk+1/2}\right)$	$\displaystyle=$
$\displaystyle\frac{1}{2}\left(D_{-}^{x}(\overline{\|{\bf v}\|^{2}}_{i+1/2jk}{\bf f}^{1}_{i+1/2jk})-\frac{1}{2}(D_{-}^{x}\|{\bf v}\|^{2}_{i+1jk}){\bf f}^{1}_{i+1/2jk}-\frac{1}{2}(D_{-}^{x}\|{\bf v}\|^{2}_{ijk}){\bf f}^{1}_{i-1/2jk}\right)$
$\displaystyle+\frac{1}{2}\left(D_{-}^{y}(\overline{\|{\bf v}\|^{2}}_{ij+1/2k}{\bf g}^{1}_{ij+1/2k})-\frac{1}{2}(D_{-}^{y}\|{\bf v}\|^{2}_{ij+1/2k}){\bf g}^{1}_{ij+1/2k}-\frac{1}{2}(D_{-}^{y}\|{\bf v}\|^{2}_{ijk}){\bf g}^{1}_{ij-1/2k}\right)$		(88)
$\displaystyle+\frac{1}{2}\left(D_{-}^{z}(\overline{\|{\bf v}\|^{2}}_{ijk+1/2}{\bf h}^{1}_{ijk+1/2})-\frac{1}{2}(D_{-}^{z}\|{\bf v}\|^{2}_{ijk+1/2}){\bf h}^{1}_{ijk+1/2}-\frac{1}{2}(D_{-}^{z}\|{\bf v}\|^{2}_{ijk}){\bf h}^{1}_{ijk-1/2}\right).$

Analysis of an alternative Navier-Stokes system: Weak entropy solutions and a convergent numerical scheme

Abstract

1 Introduction

Remark 1.

1.1 Adding more physics to the model

Remark 2.

1.2 Outline

1.3 Notation

2 Additional balance laws

2.1 Entropy balance

2.2 Kinetic Energy balance

2.3 Specific volume balance

2.4 Renormalised internal energy balance

3 A priori estimates for the full system

3.1 Entropy estimate

3.2 Kinetic energy estimate

3.3 Estimates of the specific volume

Remark 3.

3.4 Estimates from the continuity equation

Remark 4.

3.5 Renormalised internal energy

3.6 Improved estimates

3.7 Compactness of ρ,𝐯,T\rho,{\bf v},T

Remark 5.

Remark 6.

3.8 Weak sequential convergence

3.9 Weak entropy solution

Definition 1.

4 The numerical approximation scheme

Remark 7.

Remark 8.

4.1 The numerical scheme

Remark 9.

Remark 10.

Lemma 1.

Proof.

Theorem 1.

5 Additional balance laws

5.1 Kinetic and internal energy balances

6 Stability

6.1 Conservation

Remark 11.

6.2 Entropy

Remark 12.

6.3 Kinetic energy

6.4 Specific volume

Remark 13.

6.5 Estimates from the continuity equation

6.6 Renormalised internal energy

Remark 14.

6.7 Improved estimates

7 Convergence to a weak solution

7.1 Solvability of the semi-discrete scheme

7.2 Convergence a.e.

7.3 Weak convergence of fluxes

8 Concluding remarks

9 Acknowledgements

References

I Prerequisites

Theorem 2.

Proof.

Lemma 2.

Lemma 3.

Lemma 4.

Theorem 3.

3.7 Compactness of $\rho,{\bf v},T$