# Modelling of a spherical deflagration at constant speed

D Grapsas, R Herbin, J.-C Latché, Y Nasseri

► To cite this version:

D Grapsas, R Herbin, J.-C Latché, Y Nasseri. Modelling of a spherical deflagration at constant speed. Computational & Applied Mathematics, 2021, 40 (163), 10.1007/s40314-021-01543-1. hal-02967980v2

HAL Id: hal-02967980

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

Submitted on 13 Apr 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
MODELLING OF A SPHERICAL DEFLAGRATION AT CONSTANT SPEED

D. GRAPSAS<fnref n="1" />, R. HERBIN<fnref n="2" />, J.-C. LATCHÉ<fnref n="3" /> AND Y. NASSERI<fnref n="4" />

**Abstract.** In this paper, we build a numerical solution procedure to compute the flow induced by a spherical flame expanding from a point source at a constant expansion velocity, with an instantaneous chemical reaction. The solution is supposed to be self-similar and the flow is split in three zones: an inner zone composed of burnt gases at rest, an intermediate zone where the solution is regular and the initial atmosphere composed of fresh gases at rest. The intermediate zone is bounded by the reactive shock (inner side) and the so-called precursor shock (outer side), for which Rankine-Hugoniot conditions are written; the solution in this zone is governed by two ordinary differential equations which are solved numerically. We show that, for any admissible precursor shock speed (and, in some cases, for a reaction heat large enough), the construction combining this numerical resolution with the exploitation of jump conditions is unique, and yields decreasing pressure, density and velocity profiles in the intermediate zone. In addition, the reactive shock speed is larger than the velocity on the outer side of the shock, which is consistent with the fact that the difference of these two quantities is the so-called flame velocity, i.e. the (relative) velocity at which the chemical reaction progresses in the fresh gases. Finally, we also observe numerically that the function giving the flame velocity as a function of the precursor shock speed is increasing; this allows to embed the resolution in a Newton-like procedure to compute the flow for a given flame speed (instead of for a given precursor shock speed). The resulting numerical algorithm is applied to stoichiometric hydrogen-air mixtures.

**2010 AMS Subject Classification.** Primary 76N15.

April 13, 2021.

# 1. PROBLEM POSITION

Explosions, and, more specifically, deflagration hazards are one of the main threats in industry processes. For instance, in the context of electricity production by nuclear energy, the accident at Japan's Fukushima Daiichi power plant in March 2011 demonstrated the devastating nature of a hydrogen explosion, particularly for the containment of radioactive materials. Since then, IRSN has been developing its own simulation tool, the P²REMICS software, to better handle the complexity of the modeled phenomena and the relevance of the performed calculations. The computer code P²REMICS can simulate three-dimensional, incompressible (explosive atmosphere formation phase) or compressible (deflagration and shock wave propagation), turbulent and reactive flows. Turbulence can be handled by statistical models (RANS), by large scale simulation (LES), or by so-called hybrid models, achieving a continuous transition between these two approaches. For the treatment of deflagration in large rooms, P²REMICS proposes an explicit flame front tracking model, requiring a less fine discretization in space than the resolution of the primitive conservation equations of reactive fluid mechanics. The verification of P²REMICS is based on a large number of analytical tests. One of them is a mathematically constructed exact solution of a model constant-speed flame propagation in a reactive infinite atmosphere of constant initial composition, which is the aim of this paper.

The ignition is supposed to occur at a single point (chosen to be the origin of $\mathbb{R}^3$) and the flow is supposed to satisfy a spherical symmetry property: the density $\rho$, the pressure $p$, the internal energy $e$ and the entropy $s$ only depend on the distance $r$ to the origin and the velocity reads $u = ur/r$, where $r$ stands for the position

<footnote>Keywords and phrases: spherical flames, reactive Euler equations, Riemann problems</footnote>

<footnote n="1">¹ I2M UMR 7373, Aix-Marseille Université, CNRS, Ecole Centrale de Marseille. 39 rue Joliot Curie. 13453 Marseille, France.
(dionysis.grapsas@gmail.com)</footnote>

<footnote n="2">² I2M UMR 7373, Aix-Marseille Université, CNRS, Ecole Centrale de Marseille. 39 rue Joliot Curie. 13453 Marseille, France.
(raphaele.herbin@univ-amu.fr)</footnote>

<footnote n="3">³ IRSN, BP 13115, St-Paul-lez-Durance Cedex, France (jean-claude.latche@irsn.fr)</footnote>

<footnote n="4">⁴ I2M UMR 7373, Aix-Marseille Université, CNRS, Ecole Centrale de Marseille. 39 rue Joliot Curie. 13453 Marseille, France.
(youssouf.nasseri@univ-amu.fr)</footnote>
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Structure of the solution.</figcaption></figure>

vector. The flame is supposed to be infinitely thin and to move at a constant speed. The flow is governed by the Euler equations, and we seek a solution with the following structure:

- the solution is self-similar, i.e. the quantities $\rho, p, e, s$ and $u$ are functions of the variable $x = r/t$ only.

- the flow is split in three zones, referred to as the inner, intermediate and outer zones. The inner zone stands for the burnt zone while, in the other two zones, the gas is supposed to be in its initial (referred to as fresh or unburnt) composition. Burnt and fresh gases differ by the expression of the total energy:

$$E = \frac{1}{2}u^2 + e - \zeta_b Q, \quad \zeta_b = 1 \text{ in the burnt zone, } \zeta_b = 0 \text{ in the fresh zone,} \quad (1)$$

with $Q > 0$ the chemical heat reaction. Both burnt and unburnt gases are considered as ideal gases, possibly with different heat capacity ratios:

$$p = (\gamma - 1)\rho e, \quad \gamma = \gamma_b \text{ for burnt gases, } \gamma = \gamma_u \text{ for unburnt gases.}$$

- In the burnt zone, the solution is supposed to be constant; this constant state is denoted by $W_b = (\rho_b, u_b, p_b)$. For symmetry reasons, the fluid is at rest in this zone, i.e. $u_b = 0$.

- The burnt and intermediate zones are separated by a shock, which coincides with the flame front. This shock is called the reactive shock, and travels at a constant speed $\sigma_r$. The outer state of the shock is denoted by $W_2 = (\rho_2, u_2, p_2)$. Note that the usual Rankine-Hugoniot jump conditions apply at the reactive shock, up to the fact that the expression of the total energy in the inner and outer states differs (see Equation (1)).

- The intermediate and the outer zones are separated by a 3-shock, referred to as the precursor shock, and travelling at a velocity denoted by $\sigma_p$. We denote by $W_1 = (\rho_1, u_1, p_1)$ the inner state of the precursor shock, and, since the usual jump conditions for the Euler equations apply, we have $\sigma_p \ge u_1$. In the outer zone, conditions are constant and equal to the initial condition $W_0 = (\rho_0, u_0, p_0)$; the fluid is at rest, i.e. $u_0 = 0$.

- In the intermediate zone, the states $W_2$ and $W_1$ are supposed to be linked by a regular solution.

In addition, for physical reasons, we expect that

$$u_2 > 0 \text{ and } \sigma_r = u_2 + u_f \text{ with } u_f > 0. \quad (2)$$

Indeed, the velocity $u_f$ is the velocity at which the chemical reaction progresses in the fresh gases; these are pushed away from the origin by the expansion of the burnt gases, and therefore $u_2 > 0$.
The aim of this paper is to build a numerical procedure to compute a solution with the above described structure. More precisely speaking, we present the two following developments:

- First, for a given precursor shock speed $\sigma_p$, we derive a solution with the desired structure in a constructive way (and this construction yields a unique solution), and propose a simple numerical scheme to compute it. Moreover, the constructed solution is such that the inequalities (2) are satisfied (in fact, we obtain that $\sigma_r - u_2 > 0$ since we seek and find a solution such that $x - u(x) > 0$ in the whole intermediate zone), and thus yields a physically meaningful flame velocity $u_f$.

- As a by-product, we numerically obtain the velocity $u_f$ as a function of $\sigma_p$, i.e. we construct a function $\tilde{G}$ such that $u_f = \tilde{G}(\sigma_p)$, and observe that this function $\tilde{G}$ is strictly increasing, which was expected from physical reasons (the faster the combustion, the stronger the generated shock-wave). It is thus easy to build an iteration to compute the flow associated to a given $u_f$, which is generally the problem of physical interest.

Finally, this process is implemented in the free CALIF³S software [2] developed at the French Institut de Radio-protection et de Sûreté Nucléaire (IRSN) and applied to obtain solutions as a function of $u_f$ for a stoichiometric mixture of hydrogen and air.

The derivation of a solution for the same problem may be found in [7]; however, the techniques used in this latter paper are different (the solution is performed in the phase space), and the uniqueness of the construction together with the proof of the decreasing properties of the solution are not explicit. The developments in [7] are built upon techniques developed for non-reactive problems in [9, 11]. Approximate solutions in closed form are given in [3, 6], and extensions to accelerating flames may be found in [4, 10]. Finally, the complete solution of the plane case (i.e. the one-dimensional case in cartesian coordinates), in closed form, is given in [1].

## 2. SOLUTION FOR A GIVEN PRECURSOR SHOCK SPEED

We first propose a constructive derivation of the solution for a given precursor shock speed, and then state the numerical scheme to compute it.

### 2.1. Derivation of the solution

Since the fluid state in the outer zone $W_0$ is given, we may equivalently use hereafter either $\sigma_p$ or the precursor shock Mach number defined by $M_p = \sigma_p/c_0$, with $c_0 = (\gamma_u p_0/\rho_0)^{1/2}$ the speed of sound in the outer zone. We recall that $M_p > 1$ thanks to the Lax entropy condition. By a standard computation, the Rankine-Hugoniot conditions at the precursor shock yields the state $W_1$ (see Appendix B).

$$ \rho_1 = \frac{\gamma_u + 1}{\gamma_u - 1 + \frac{2}{M_p^2}} \rho_0, \quad u_1 = (1 - \frac{\rho_0}{\rho_1}) \sigma_p, \quad p_1 = p_0 + (1 - \frac{\rho_0}{\rho_1}) \rho_0 \sigma_p^2. \qquad (3) $$

In addition, using the fact that the fluid is at rest in the burnt zone, we show in Appendix C that State 2 must satisfy the following relation:

$$ \mathcal{F}_r(\sigma_r) := \frac{1}{2} u_2^2 + \frac{1}{\gamma_b - 1} u_2 \sigma_r + \left( \frac{\gamma_u}{\gamma_u - 1} - \frac{\gamma_b}{\gamma_b - 1} \frac{\sigma_r}{\sigma_r - u_2} \right) \frac{p_2}{\rho_2} + Q = 0, \qquad (4) $$

where $\sigma_r$ stands for the (unknown) reactive shock speed. In addition, the same jump conditions give the burnt state $W_b$ as a function of $W_2$:

$$ \rho_b = \rho_2(\frac{\sigma_r - u_2}{\sigma_r}), \quad p_b = p_2 - \rho_2 u_2(\sigma_r - u_2). \qquad (5) $$

To complete the derivation of the solution, we must now show that the following program makes sense: starting from $x = \sigma_p$, solve the Euler equations for $x \le \sigma_p$ until the point $x = \sigma_r$ where Equation (4) is verified. The solution at this point is equal to $W_2$ and Equations (5) yield the burnt state $W_b$. Let us now embark on this development.
Governing equations in the intermediate zone. Since we suppose that the solution is regular in this zone, we may replace the total energy balance in the Euler equations by the entropy equation, which, under the spherical symmetry assumption, yields the following system (see Appendix A):

$$ \partial_t(r^2\rho) + \partial_r(r^2\rho u) = 0, \qquad (6a) $$

$$ \partial_t(r^2\rho u) + \partial_r(r^2(\rho u^2 + p)) = 2rp, \qquad (6b) $$

$$ \partial_t(r^2\rho s) + \partial_r(r^2\rho su) = 0. \qquad (6c) $$

The mass balance equation (6a) may be developed to obtain:

$$ \partial_t \rho + u \partial_r \rho + \rho \partial_r u + \frac{2}{r} \rho u = 0. \qquad (7) $$

In addition, thanks to the mass balance equation and for a regular function $f = f(t, r)$, we have:

$$ \partial_t(r^2\rho f) + \partial_r(r^2\rho fu) = r^2\rho(\partial_t f + u \partial_r f). $$

Using this identity in the momentum and entropy balances, i.e. Equation (6b) and (6c) respectively, we get:

$$ \begin{aligned} \partial_t u + u \partial_r u + \frac{1}{\rho} \partial_r p &= 0, \\ \partial_t s + u \partial_r s &= 0. \end{aligned} \qquad (8) $$

We now use the fact that, if a regular function $\varphi(t, r)$ only depends on $x = r/t$, which means that there exists $\tilde{\varphi}: \mathbb{R} \to \mathbb{R}$ such that $\tilde{\varphi}(x) = \varphi(t, r)$, we have

$$ \partial_t \varphi(t, r) = -\frac{r}{t^2} \tilde{\varphi}'(x) \quad \text{and} \quad \partial_r \varphi(t, r) = \frac{1}{t} \tilde{\varphi}'(x). $$

Since we look for a self-similar solution, we may apply this identity to (7) and (8). Keeping the same notation for functions of the pair $(t, r)$ and $x$ for short, we obtain the following system:

$$ \frac{-x + u}{\rho} \rho'(x) + u'(x) + \frac{2u(x)}{x} = 0, \qquad (9a) $$

$$ (-x + u(x)) u'(x) + \frac{1}{\rho(x)} p'(x) = 0, \qquad (9b) $$

$$ (u(x) - x) s'(x) = 0. \qquad (9c) $$

Let us now suppose that $u < x$ in the intermediate zone. Note that the fact that the precursor shock is a 3-shock implies that $u_1 < \sigma_p$, so that the assumed inequality is true in the outer boundary of the intermediate zone, and the assumption amounts to state that the intermediate zone may be supposed to end, when constructing the solution, before $u = x$ occurs; this is checked a posteriori. The last relation thus implies that the entropy remains constant over the zone:

$$ s = \frac{p}{\rho^{\gamma_u}} = s_1 = \frac{p_1}{\rho_1^{\gamma_u}}, \qquad (10) $$

and this is a known value thanks to (3). We thus have $p' = \gamma_u s_1 \rho^{\gamma_u - 1} \rho'$; using $c^2 = \gamma_u p/\rho = \gamma_u s_1 \rho^{\gamma_u - 1}$, we thus get $p' = c^2 \rho'$. Substituting this expression in (9a)-(9b) and solving for $\rho'$ and $u'$, we get:

$$ \rho'(x) = - \frac{2u(x)(u(x) - x)}{x((u(x) - x)^2 - c(x)^2)} \rho, \qquad (11a) $$

$$ u' = \frac{2c(x)^2}{x((u(x) - x)^2 - c(x)^2)} u. \qquad (11b) $$
This system of coupled ODEs is complemented by initial conditions, which consist in the data of the velocity and the density at the precursor shock, i.e. at the outer boundary of the intermediate zone $x = \sigma_p$:

$$ \rho(\sigma_p) = \rho_1, u(\sigma_p) = u_1. \qquad (12) $$

Existence, uniqueness and properties of the solution. We begin by proving an a priori property of the solution, namely the fact that $\rho(x)$ and $u(x)$ are necessarily decreasing functions in the intermediate zone. To this end, we will invoke the following easy lemma, which is a consequence of the mean value theorem.

**Lemma 2.1.** Let $h$ be a continuously differentiable real function, let us suppose that there exists $a > 0$ such that $h(a) > 0$, and that $h$ satisfies the property $h'(x) \le 0$ if $h(x) > 0$. Then $h(x) \ge h(a)$, for all $x \le a$.

From the expression (3), we know that $0 < u_1 < \sigma_p$ and $\rho_1 > 0$. Let us now introduce $\sigma_\ell$ as the smallest real number in $[0, \sigma_p)$ such that, for $x \in [\sigma_\ell, \sigma_p]$, $0 \le u \le x$ and $\rho \ge 0$. Note that such a closed interval exists by the continuity (assumed in this zone) of $\rho$ and $u$. We are now in position to state the following result.

**Lemma 2.2 (Variations of the solution).** Let us suppose that the pair $(\rho, u)$ satisfies (11). Then $\rho$ and $u$ are two decreasing functions over $[\sigma_\ell, \sigma_p]$. Consequently, $\rho \ge \rho_1 > 0$ and $u \ge u_1 > 0$ over $[\sigma_\ell, \sigma_p]$.

*Proof.* Let us consider the function $h: \mathbb{R}_+ \to \mathbb{R}$ defined by $h(x) = u(x) + c(x) - x$. First, we remark that, by the Lax entropy condition, we have $h(\sigma_p) = u_1 + c_1 - \sigma_p > 0$ (and this property may be checked using the expressions (3) of $W_1$). Second, if $h(x) > 0$, since by assumption $u(x) \le x$, $\rho \ge 0$ and $c \ge 0$, we have:

$$ (u(x) - x)^2 - c(x)^2 = (u(x) - x - c(x))(u(x) - x + c(x)) \le 0. $$

Equations (11a) and (11b) thus readily imply that $\rho' \le 0$ and $u' \le 0$ since, still by assumption, $u \ge 0$ and $\rho \ge 0$. The function $h$ is thus the sum of three non-increasing functions, and is hence non-increasing itself. Lemma 2.1 applies, and yields $h(x) = u(x) - x - c(x) \ge h(\sigma_p) > 0$ over the whole interval $[\sigma_\ell, \sigma_p]$, which in turn implies $\rho' \le 0$ and $u' \le 0$. We thus have $\rho \ge \rho_1 > 0$ and $u \ge u_1 > 0$, which finally yields $\rho' < 0$ and $u' < 0$ over $[\sigma_\ell, \sigma_p]$. $\square$

Note that the inequality $h(x) \ge h(\sigma_p) > 0$ derived in this proof implies that the denominator in Equations (11a) and (11b) does not vanish in the interval $[\sigma_\ell, \sigma_p]$. The right-hand side of System (11) is thus a $C^\infty$ function of $\rho$, $u$ and $x$, and the existence and uniqueness of a solution follows by the Cauchy-Lipschitz theorem. This result is stated in the following lemma.

**Lemma 2.3 (Existence and uniqueness of the solution).** There exists one and only one solution $(\rho, u)$ of System (11)-(12) over the interval $[\sigma_\ell, \sigma_p]$.

In addition, Lemma 2.2 allows to characterize $\sigma_\ell$. Indeed, since $u \ge u_1 > 0$ and $\rho \ge \rho_1 > 0$ over $[\sigma_\ell, \sigma_p]$, by definition of $\sigma_\ell$, either $\sigma_\ell = 0$ or $u(\sigma_\ell) = \sigma_\ell$. Since $u \ge u_1 > 0$ and $u(x) \le x$ over $[\sigma_\ell, \sigma_p]$, the first option cannot hold, and we get

$$ u(\sigma_\ell) = \sigma_\ell. \qquad (13) $$

To complete the construction of a solution, it now remains to show the existence of a real number $\sigma_r$ satisfying (4), i.e. the fact that there exists $x \in (\sigma_\ell, \sigma_p)$ such that $W(x) = (\rho(x), u(x), p(x))$ satisfies the condition $\mathcal{F}_r(x) = 0$, where $\mathcal{F}_r$ is given in (4):

$$ \mathcal{F}_r(x) = \frac{1}{2}u(x)^2 + \frac{1}{\gamma_b - 1}xu(x) + \left(\frac{\gamma_u}{\gamma_u - 1} - \frac{\gamma_b}{\gamma_b - 1}\frac{x}{x - u(x)}\right)\frac{p(x)}{\rho(x)} + Q. \qquad (14) $$

The existence of $\sigma_r$ is stated in the following lemma.

**Lemma 2.4 (Existence of $\sigma_r$).** The function $\mathcal{F}_r$ is defined and continuously differentiable on the interval $(\sigma_\ell, \sigma_p]$, and $\lim_{x \to \sigma_\ell^+} \mathcal{F}_r(x) = -\infty$. In addition, $\mathcal{F}_r(\sigma_p) > 0$ when $\gamma_u = \gamma_b$ or when the reaction heat $Q$ is large enough, according to Equation (16) below. Consequently, under one of these conditions, the set $S_r = \{x \in (\sigma_\ell, \sigma_p) \text{ such that } \mathcal{F}_r(x) = 0\}$ is a non-empty closed subset of $(\sigma_\ell, \sigma_p)$ which admits a maximal element $\sigma_r$.
Proof. When $x$ tends to $\sigma_\ell$, we have seen that $u(\sigma_\ell)$ tends to $\sigma_\ell$ and thus $\mathcal{F}_r$ tends to $-\infty$. When $x = \sigma_p$, if $\gamma_u = \gamma_b$, we observe that:

$$ \mathcal{F}_r(\sigma_p) = Q > 0. $$

This relation is a consequence of the fact that, with $Q = 0$, $\mathcal{F}_r(\sigma) = 0$ is the relation satisfied by one adjacent state of a shock travelling at the speed $\sigma$ when the velocity in the other adjacent state is zero, which is precisely the case of the state $W_1$ with $\sigma = \sigma_p$. It may also be checked by injecting the relations (3) in the definition of $\mathcal{F}_r$. When $\gamma_u \neq \gamma_b$, we thus get:

$$ \mathcal{F}_r(\sigma_p) = \left[ \frac{1}{\gamma_b - 1} - \frac{1}{\gamma_u - 1} \right] \sigma_p u_1 - \left[ \frac{\gamma_b}{\gamma_b - 1} - \frac{\gamma_u}{\gamma_u - 1} \right] \frac{\sigma_p}{\sigma_p - u_1} \frac{p_1}{\rho_1} + Q. \quad (15) $$

Let us recast Relations (3) as

$$ \rho_1 = \frac{1}{\alpha} \rho_0, \quad u_1 = (1-\alpha)\sigma_p, \quad p_1 = p_0 + (1-\alpha)\rho_0 \sigma_p^2, \quad \text{with } \alpha = \frac{\gamma_u - 1 + \frac{2}{M_p^2}}{\gamma_u + 1}. $$

Using these relations in (15), we get:

$$ \mathcal{F}_r(\sigma_p) = \left[ \frac{1 - (1 - \alpha) \gamma_b}{\gamma_b - 1} - \frac{1 - (1 - \alpha) \gamma_u}{\gamma_u - 1} \right] (1 - \alpha) c_0^2 M_p^2 - \left[ \frac{\gamma_b}{\gamma_b - 1} - \frac{\gamma_u}{\gamma_u - 1} \right] \frac{p_0}{\rho_0} + Q. \quad (16) $$

Thanks to the Lax entropy conditions, $M_p > 1$, and, for $M_p \in (1, +\infty)$, the function $\alpha(M_p)$ decreases from 1 to $(\gamma_u - 1)/(\gamma_u + 1)$. Depending on the values of $\gamma_b$ and $\gamma_u$, the quantity $1 - (1 - \alpha)\gamma_b$ may become negative when $M_p \to +\infty$ for admissible values of $\alpha$, and thus the first term may tends to $-\infty$; however, for a given $M_p$, this term is finite. The second term may be negative (still according to the values of $\gamma_b$ and $\gamma_u$) but does not depend on $M_p$. Hence, for any given $M_p$, for $Q$ large enough, the condition $F_r(\sigma_p) > 0$ is satisfied. $\square$ $\square$

In addition, when $\gamma_u = \gamma_b$, we are able to prove that $\mathcal{F}_r(\sigma_p)$ is an increasing function over $(\sigma_\ell, \sigma_p)$, and therefore the set $S_r$ contains a single point; the proof of this result is given in Appendix D.

Finally, note that $\sigma_r > \sigma_\ell$; since $u$ is a decreasing function, this yields that $\sigma_r - u(\sigma_r) > 0$. As mentioned in (2), this was expected, from a physical point of view, since the quantity $\sigma_r$ is nothing else that the flame velocity $u_f$.

## 2.2. Numerical approximation of the solution in the intermediate zone

The problem tackled in this section is twofold: first, we need to solve the system of ODEs (11)-(12) numerically, and second, we need to determine the speed of the reactive shock $\sigma_r$. To this purpose, we solve (11)-(12) by an explicit Euler scheme, starting at $N \in \mathbb{N}$ and $x^N = \sigma_p$ and, for indices $n$ decreasing from $N$, performing steps of size $-\delta x$, with $\delta x = \sigma_p/N$; at each new step $n$ associated to $x^n = n\delta x$, we obtain $W^n$ and we evaluate the function $\mathcal{F}_r$, until we obtain $\mathcal{F}_r(x^n) \le 0$. Then, the algorithm stops and we know that the computed approximation $\sigma_r^{\text{app}}$ of $\sigma_r$ satisfies $x^n < \sigma_r^{\text{app}} < x^{n+1}$; for $\delta x$ small enough, $x^{n+1}$ may thus be considered as a reasonable approximation of $\sigma_r$; this is indeed the way it is computed in the numerical experiments described below.
The scheme thus reads:

$$ \text{for } n = N, u^N = u_1, \rho^N = \rho_1, $$

for $n = N - 1$ to 0 and while $\mathcal{F}_r(x^{n+1}) > 0$,

$$ (c^{n+1})^2 = \gamma_u s_1 (\rho^{n+1})^{\gamma_u - 1}, \qquad (17) $$

$$ \rho^n = \rho^{n+1} + \delta x \frac{2 u^{n+1} (u^{n+1} - x^{n+1})}{x^{n+1} ((u^{n+1} - x^{n+1})^2 - (c^{n+1})^2)} \rho^{n+1}, $$

$$ u^n = u^{n+1} - \delta x \frac{2 (c^{n+1})^2}{x^{n+1} ((u^{n+1} - x^{n+1})^2 - (c^{n+1})^2)} u^{n+1}. $$

Then, for any valid value of $n \le N$, the pressure is given by

$$ p^n = s_1 (\rho^n)^{\gamma_u}. $$

Since the algorithm stops as soon as $\mathcal{F}_r(x^{n+1})$ becomes negative, from the expression of this latter function, we have $u^n < x^n$ in all the performed steps $n$. We thus have $u^n > 0$, $\rho^n > 0$, $u^n \ge u^{n+1}$ and $\rho^n \ge \rho^{n+1}$ at all steps.

### 3. SOLUTION FOR A GIVEN FLAME SPEED

The construction performed in the previous section shows that, to any precursor shock velocity $\sigma_p$ greater than the speed of sound $c_0$ in the outer zone of the fresh atmosphere, we are able to associate a positive flame velocity $u_f$ given by $u_f = \sigma_r - u_2$. In addition, even if we have no proof, physical arguments suggest that $u_f$ is an increasing function of $\sigma_p$ (or equivalently of the Mach number $M = \sigma_p/c_0$, considering a family of problems with the same initial atmosphere and thus $c_0$ as a fixed parameter); this behaviour is confirmed by numerical experiments (see Section 4). Computing the flow for a given $u_f$, which is in fact usually the engineering problem to be tackled, amounts to invert the function $u_f = \tilde{G}(M)$, and this equation for $M$ thus should have one and only one solution, at least for reasonable values of $u_f$. To compute this solution, we define $\mathcal{G}$ by

$$ \mathcal{G}(M) := \tilde{G}(M) - u_f \qquad (18) $$

and search for $M$ such that $\mathcal{G}(M) = 0$ with the following iterative algorithm depending on the parameters $M_0$, $\delta$ and $\epsilon$:

initialization: let $M_0$ be given, and compute $\mathcal{G}(M_0)$,

let $M_1 = M_0 + \delta$, and compute $\mathcal{G}(M_1)$

current iteration: For $k \ge 2$, let $M_k = M_{k-1} - \frac{M_{k-1} - M_{k-2}}{\mathcal{G}(M_{k-1}) - \mathcal{G}(M_{k-2})} \mathcal{G}(M_{k-1})$, and compute $\mathcal{G}(M_k)$.

stopping criteria: stop when $\mathcal{G}(M_k) \le \epsilon$.

This algorithm is used in the following section with $M_0 = 1.0001$, $\delta = 0.001$ and $\epsilon = 10^{-5}$. Convergence is obtained for all cases, provided that the number of cells used in the numerical computation of the solution in the intermediate zone is large enough; otherwise, the error on $\sigma_r$ is too large and the prescribed tolerance threshold for the value of $\mathcal{G}$ cannot be reached.

### 4. APPLICATION TO HYDROGEN DEFLAGRATIONS

We now apply the above developed procedure for a flame propagating in a stoichiometric mixture of hydrogen and air. We consider a unique total and irreversible chemical reaction, which reads:

$$ 2 \text{H}_2 + \text{O}_2 \rightarrow 2 \text{H}_2\text{O}. $$
Supposing that air is composed of 1/5 of oxygen and 4/5 nitrogen (molar or volume proportions), the molar fractions of hydrogen, oxygen and nitrogen in the considered stoichiometric mixture are thus equal to 2/7, 1/7 and 4/7 respectively. The mass fractions of these constituents are thus easily deduced from these values:

$$y_{H_2} = \frac{2W_{H_2}}{W_t}, \quad y_{O_2} = \frac{W_{O_2}}{W_t}, \quad y_{N_2} = \frac{4W_{N_2}}{W_t}, \quad W_t = 2W_{H_2} + W_{O_2} + 4W_{N_2},$$

where $W_{H_2} = 0.002$ Kg, $W_{O_2} = 0.032$ Kg and $W_{N_2} = 0.028$ Kg stand for the molar mass of the hydrogen, oxygen and nitrogen molecules respectively. Since $H_2$ and $O_2$ are pure substances (and thus their formation enthalpy is equal to zero), the chemical reaction heat reads:

$$Q = y_{H_2O} \Delta H_0^f = (y_{H_2} + y_{O_2}) \Delta H_0^f,$$

where $\Delta H_0^f = 1.3255 \times 10^7$ J/Kg stands for the formation enthalpy of steam. The initial pressure is $p = 10^5$ Pa, the initial temperature is $T = 283$°K, the initial density is given by the Boyle-Mariotte law and the heat capacity ratio is $\gamma_u = \gamma_b = 1.4$.

We plot on Figures 2 and 3 the density, velocity and pressure profiles obtained for $u_f = 32$ m/s and $u_f = 4$ m/s respectively.

The solution in the intermediate zone is obtained with a regular mesh, splitting the interval between 5 m and the position of the precursor shock (which is unknown up to the last solution step of the algorithm) in 5000 equal subintervals.

Then we show on Figures 4-7 the evolution of the states $W_1$, $W_2$ and $W_b$ as a function of the flame velocity. The temperature in the burnt state, not shown here, is close to $T = 3050$°K for all the values of the flame velocity $u_f$.

We observe that the precursor shock is of very weak amplitude for low values of the flame velocity; in fact, it becomes visible only when $u_f$ reaches 20 m/s. For $u_f = 4$, the computed velocity at state $W_1$ is lower than $10^{-6}$ m/s, while it reaches values greater than 30 m/s at state $W_2$. Since the ordinary differential equation governing the velocity in the intermediate zone (11b) is of the form

$$u' = f(\rho, u) u,$$

one may anticipate such a low value as initial (right) condition to lead to severe accuracy problems. In this respect, the computed value which seems to be the most affected is the velocity at state $W_2$: the convergence value seems to be close to 33.00 m/s, we obtain $u_2 \simeq 34.5$ m/s with $n = 5 \times 10^3$ cells and $u_2 \in (32.95$ m/s, 33 m/s) for $n = 8 \times 10^4$, $n = 16 \times 10^4$, $n = 32 \times 10^4$ and $n = 64 \times 10^4$. As expected, convergence is easier when the precursor shock has a significant amplitude: $u_2 \simeq 243.0$ m/s for $n = 5000$, for a convergence value in the range of 243.8 m/s.

## 5. CONCLUSION

In this paper, we were able to construct an iterative procedure in order to obtain a reference solution for the reactive Euler equations in spherical coordinates. This reference solution is now used as a test for the development and improvement of the P²REMICS software at IRSN which simulates three-dimensional, incompressible or compressible turbulent and reactive flows, including deflagration. Note that P²REMICS has also been tested on several experimental data. The source code, the documentation and validation tests are all freely available on the gforge platform at IRSN [8].
<figure><img src="image_3.png" /><figcaption>FIGURE 2. Density (kg m⁻³), velocity (m s⁻¹) and pressure (Pa) profiles obtained for a flame velocity of u_f = 32 m/s.</figcaption></figure>

# Appendix

## APPENDIX A. EULER EQUATIONS IN SPHERICAL COORDINATES

The Euler equations read in Cartesian coordinates:

$$ \partial_t \bar{\rho} + \text{div}(\bar{\rho} \bar{u}) = 0, \qquad (19a) $$

$$ \partial_t (\bar{\rho} \bar{u}) + \text{div}(\bar{\rho} \bar{u} \otimes \bar{u}) + \nabla \bar{p} = 0, \qquad (19b) $$

$$ \partial_t (\bar{\rho} \bar{E}) + \text{div}(\bar{\rho} \bar{E} \bar{u} + \bar{p} \bar{u}) = 0, \qquad (19c) $$

where $\bar{\rho} = \bar{\rho}(t, x) \in \mathbb{R}$ the density, $\bar{u} = \bar{u}(t, x) \in \mathbb{R}^3$ the velocity, $\bar{p} = \bar{p}(t, x) \in \mathbb{R}$ the pressure and $\bar{E} = \bar{E}(t, x) \in \mathbb{R}$ the total energy for all $t \in \mathbb{R}$ and $x = (x_1, x_2, x_3) \in \mathbb{R}^3$. This system is closed by the equation of state, which for a perfect gas, is given by

$$ \bar{E} = \frac{1}{2} |\bar{u}|^2 + \bar{e}, \quad \text{with} \quad \bar{p} = (\gamma - 1) \bar{\rho} \bar{e}, \qquad (20) $$
<figure><img src="image_4.png" /><figcaption>FIGURE 3. Density (kg m⁻³), velocity (m s⁻¹) and pressure (Pa) profiles obtained for a flame velocity of u_f = 4 m/s.</figcaption></figure>

where $\bar{e} = \bar{e}(t, \mathbf{x})$ the internal energy and $\gamma > 1$ the heat capacity ratio. We suppose that the flow satisfies a spherical symmetry assumption, so the solution of equations (19)-(20) may be recast as:

$$ \bar{\rho}(t, \mathbf{x}) = \rho(t, r), \quad \bar{p}(t, \mathbf{x}) = p(t, r), \quad \bar{E}(t, \mathbf{x}) = E(t, r) \quad \text{and} \quad \bar{u}(t, \mathbf{x}) = u(t, r) \frac{\mathbf{x}}{r}, \qquad (21) $$

with $r = |\mathbf{x}|$ and where $(\rho, u, p, E)(t, r) \in \mathbb{R}^4$ are scalar functions, i.e. $(\rho, u, p, E) \in \mathbb{R}^4$. The aim of this section is to derive the system of equations satisfied by $(\rho, u, p, E)$. We suppose first that these functions are regular, so we obtain the strong form of the so-called Euler equations in spherical coordinates; then we derive the Rankine-Hugoniot conditions that are satisfied by discontinuous spherical solutions.

Let us first note that, denoting by $\partial_i$ the derivative with respect to the i-th space coordinate, one has $\partial_i r = \frac{x_i}{r}$ for $i = 1, 2, 3$. As a consequence, if $\bar{f}: \mathbb{R}^2 \to \mathbb{R}$ is a function such that $\bar{f}(\mathbf{x}) = f(r)$, one as

$$ \nabla \bar{f}(\mathbf{x}) = \partial_r f(r) \frac{\mathbf{x}}{r}. \qquad (22) $$
<figure><img src="image_5.png" /><figcaption>FIGURE 4. Density at state 1 and state 2 as a function of the flame velocity.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>FIGURE 5. Density in the burnt zone as a function of the flame velocity.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIGURE 6. Velocity at state 1 and state 2 and speed of the precursor shock as a function of the flame velocity.</figcaption></figure>

Furthermore,

$$ 
\begin{aligned}
\mathrm{div}(\bar{f}\bar{u})(\boldsymbol{x}) &= \sum_{i=1}^{3} \partial_i (f(r)u(r)\frac{x_i}{r}) \\
&= f(r)u(r)\sum_{i=1}^{3}\left(\frac{1}{r} - \frac{x_i^2}{r^3}\right) + \sum_{i=1}^{3}\frac{x_i}{r}\partial_i r \partial_r(fu)(r) \\
&= \frac{2}{r}f(r)u(r) + \partial_r(fu)(r) = \frac{1}{r^2}\partial_r(r^2f(r)u(r)).
\end{aligned} 
\quad (23) $$
<figure><img src="image_8.png" /><figcaption>FIGURE 7. Pressure at state 1, at state 2 and in the burnt zone as a function of the flame velocity.</figcaption></figure>

**Lemma A.1** (The Euler equations in spherical coordinates). Suppose that $(\bar{\rho}, \bar{u}, \bar{p}, \bar{E})$ is solution of (19) satisfying (21); then $(\rho, u, p, E)$ satisfies:

$$ \partial_t(r^2\rho) + \partial_r(r^2\rho u) = 0, \qquad (24a) $$

$$ \partial_t(r^2\rho u) + \partial_r(r^2(\rho u^2 + p)) = 2rp, \qquad (24b) $$

$$ \partial_t(r^2\rho E) + \partial_r(r^2(\rho u E + p u)) = 0, \qquad (24c) $$

with $E = \frac{1}{2}u^2 + e$ and $p = (\gamma - 1)\rho e$.

Furthermore, for a regular solution, the following entropy equation holds

$$ \partial_t(r^2\rho s) + \partial_r(r^2\rho s u) = 0. \qquad (25) $$

*Proof.* The mass balance equation (24a) is a straightforward consequence of (23). Let us then turn to the momentum balance; by definition of $u_i$,

$$ \mathrm{div}(\rho u_i \boldsymbol{u}) = \mathrm{div}(\frac{\rho u}{r} x_i \boldsymbol{u}) = \mathrm{div}(\frac{\rho u}{r} \boldsymbol{u}) x_i + \frac{\rho u}{r} \boldsymbol{u} \cdot \nabla x_i $$

Hence, thanks to (23) with $f(r) = \rho \frac{u}{r}$ for the first term in the sum, and noting that $\boldsymbol{u} \cdot \nabla x_i = \frac{u}{r} x_i$ for the second one, we get

$$ \mathrm{div}(\rho u_i \boldsymbol{u}) = \frac{x_i}{r^2} (\partial_r(\rho u^2 r) + \rho u^2) = \frac{x_i}{r^2} (r \partial_r(\rho u^2) + 2\rho u^2) = \frac{x_i}{r^3} \partial_r(r^2 \rho u^2). $$

Thanks to this relation and using (22) for the pressure gradient, for $i = 1, 2, 3$, the $i$-th component of the momentum equation reads:

$$ \frac{x_i}{r} [\partial_t(\rho u) + \frac{1}{r^2}\partial_r(r^2\rho u^2) + \partial_r p] = 0. $$

The three components of this vector balance equation thus boil down to the single relation:

$$ \partial_t(\rho u) + \frac{1}{r^2}\partial_r(r^2\rho u^2) + \partial_r p = 0. $$

Multiplying by $r^2$, we obtain the conservative form (24b).

Let us now turn to the energy balance, noting that $\bar{E}(t, \boldsymbol{x}) = \frac{1}{2}|\boldsymbol{u}(\bar{t}, \boldsymbol{x})|^2 + e(\bar{t}, \boldsymbol{x}) = \frac{1}{2}u(t, r)^2 + e(t, r) = E(t, r)$ and using (23) once again, we get that

$$ \mathrm{div}(\bar{\rho} \bar{E} \bar{u}) = \frac{1}{r^2} \partial_r (\rho E u r^2) \quad \text{and} \quad \mathrm{div}(\bar{p} \bar{u}) = \frac{1}{r^2} \partial_r (p u r^2), $$

and (24c) follows by adding the time derivative of $(\bar{\rho} \bar{E})$.
Finally, we apply the same process for the entropy balance equation. In Cartesian coordinates, the entropy reads $\bar{s}(t, \mathbf{x}) = \frac{\bar{p}(t, \mathbf{x})}{\bar{\rho}(t, \mathbf{x})^{\gamma}}$, so that, under the spherical symmetry conditions (21),

$$ \bar{s}(t, \mathbf{x}) = s(t, r) \text{ with } s = \frac{p}{\rho^{\gamma}}, $$

The entropy relation in Cartesian coordinates reads, for regular solutions:

$$ \partial_t(\bar{\rho}\bar{s}) + \operatorname{div}(\bar{\rho}\bar{u}\bar{s}) = 0, \quad \text{with } \bar{s} = \frac{\bar{p}}{\bar{\rho}^{\gamma}}, $$

with $\bar{s} = \bar{s}(t, \mathbf{x})$. A straightforward application of (23) with $f = \rho s$ then yields (25). □

Let us now turn to discontinuous solutions and show how the Rankine Hugoniot conditions read in spherical coordinates. For the sake of simplicity, we consider here a weak solution $\bar{\rho}$ of the mass equation $\partial_t\bar{\rho} + \operatorname{div}(\bar{\rho}\bar{u}) = 0$, with initial condition $\bar{\rho}_0(\mathbf{x}) = \bar{\rho}(0, \mathbf{x})$ i.e.

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}} \bar{\rho} \partial_t \bar{\varphi} + \bar{\rho} \bar{u} \cdot \nabla \bar{\varphi} \, d\mathbf{x} \, dt + \int_{\mathbb{R}^d} \bar{\rho}_0(\mathbf{x}) \bar{\varphi}(0, \mathbf{x}) \, d\mathbf{x} = 0, \quad \forall \bar{\varphi} \in C_c^\infty(\mathbb{R}_+ \times \mathbb{R}^d, \mathbb{R}), \quad (26) $$

The aim is to show that the Rankine Hugoniot condition holds for a discontinuous function under spherical form.

**Lemma A.2** (Rankine Hugoniot condition in spherical coordinates). Let $\Omega$ be an open subset of $\mathbb{R}_+^*$ et let $\bar{\rho}$ satisfying (26) and such that $\bar{\rho}(t, \mathbf{x}) = \rho(t, r)$ and $\bar{u}(t, \mathbf{x}) = u(t, r)\frac{\mathbf{x}}{r}$, Let $\bar{\varphi} \in C_c^\infty(\mathbb{R}_+ \times \mathbb{R}^d, \mathbb{R})$ and $\varphi \in C_c^\infty(\mathbb{R}_+ \times \mathbb{R}_+^*, \mathbb{R})$ such that $\bar{\varphi}(t, \mathbf{x}) = \varphi(t, r)$. Then $\rho$ satisfies the following weak form of the mass equation in spherical coordinates.

$$ \int_{\mathbb{R}_+} \int_{\mathbb{R}_+^*} (r^2 \rho \partial_t \varphi + r^2 \rho u \partial_r \varphi) \, dr \, dt + \int_{\mathbb{R}_+^*} r^2 \rho_0(r) \varphi(0, r) \, dr = 0. \quad (27) $$

Furthermore, let $\rho$ and $u$ be piecewise regular functions of the form:

$$ \rho(t, r) = \begin{cases} \rho_L(t, r) & \text{if } r \in \Omega^{-} \\ \rho_R(t, r) & \text{if } r \in \Omega^{+} \end{cases} \quad \text{with } u(t, r) = \begin{cases} u_L(t, r) & \text{if } r \in \Omega^{-} \\ u_R(t, r) & \text{if } r \in \Omega^{+} \end{cases} \quad (28) $$

where $\rho_L, u_L$ (resp. $\rho_R, u_R$) are smooth functions satisfying (24a) on $\Omega^{-} = \{r \in \Omega, r \le t\sigma\}$ (resp. $\Omega^{+} = \{r \in \Omega, r > t\sigma\}$), with $\sigma \in \mathbb{R}$. Then the following Rankine-Hugoniot condition holds:

$$ \sigma[\rho] = [\rho u], \quad (29) $$

where $[\cdot]$ denotes the jump of a given possibly discontinuous function at the interface $r = \sigma t$ between $\Omega^{-}$ and $\Omega^{+}$.

**Proof.** Let $\bar{\varphi} \in C_c^\infty(\mathbb{R}_+ \times \mathbb{R}^d, \mathbb{R})$ and let $\varphi \in C_c^\infty(\mathbb{R}_+ \times \mathbb{R}_+^*, \mathbb{R})$ be defined by $\bar{\varphi}(t, \mathbf{x}) = \varphi(t, r)$. Then

$$ \begin{aligned} & \int_{\mathbb{R}_+} \int_{\mathbb{R}^d} \bar{\rho} (\partial_t \bar{\varphi} + \bar{\rho} \bar{u} \cdot \nabla \bar{\varphi}) \, d\mathbf{x} \, dt \\ &= \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} (\rho(t, r) \partial_t \varphi(t, r) + \rho(t, r) u(t, r) \frac{\mathbf{x}}{r} \cdot \nabla \varphi(t, r)) r^2 \, dr \, dt \\ &= \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} (r^2 \rho(t, r) \partial_t \varphi(t, r) + r^2 \rho(t, r) u(t, r) \frac{\mathbf{x}}{r} \cdot \frac{\mathbf{x}}{r} \partial_r \varphi(t, r)) \, dr \, dt \\ &= \int_{\mathbb{R}_+} \int_{\mathbb{R}_+} (r^2 \rho(t, r) \partial_t \varphi(t, r) + r^2 \rho(t, r) u(t, r) \partial_r \varphi(t, r)) \, dr \, dt. \end{aligned} $$
We also have that

$$\int_{\mathbb{R}^d} \bar{\rho}_0(\mathbf{x}) \bar{\varphi}(0, \mathbf{x}) d\mathbf{x} = \int_{\mathbb{R}_+} r^2 \rho_0(r) \varphi(0, r) dr$$

Thus (26) reads:

$$\int_{\mathbb{R}_+} \int_{\mathbb{R}_+} (r^2 \rho(t, r) \partial_t \varphi(t, r) + r^2 \rho(r, t) u(t, r) \partial_r \varphi(t, r)) dr dt + \int_{\mathbb{R}_+} r^2 \rho_0(r) \varphi(0, r) dr = 0, \\ \forall \varphi \in C_c^\infty(\mathbb{R}_+ \times \mathbb{R}_+, \mathbb{R}). \quad (30)$$

According to the definition (28) and decomposing the space integral on $\Omega^+$ and $\Omega^-$, we obtain that

$$\int_{\mathbb{R}_+} \int_{\Omega^+} (r^2 \rho_R \partial_t \varphi + r^2 \rho_R u_R \partial_r \varphi) \, dr \, dt + \int_{\mathbb{R}_+} \int_{\Omega^-} (r^2 \rho_L \partial_t \varphi + r^2 \rho_L u_L \partial_r \varphi) \, dr \, dt = 0,$$

which classically implies (29).

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

## APPENDIX B. LEFT STATE OF A SHOCK AS A FUNCTION OF THE RIGHT STATE AND THE SHOCK VELOCITY

In this section, we recall a classical computation (see e.g. [5]) which consists in determining the left state of a shock as a function of the right state and the shock velocity.

**Lemma B.1.** Let $W = (\rho, u, p)$ be the left state of a shock travelling at the given speed $\sigma$. Let $W_R = (\rho_R, u_R, p_R)$ be the given right state, which is supposed to satisfy $u_R = 0$. Let $c_R$ be the speed of sound in the right state, i.e. $c_R^2 = \gamma p_R / \rho_R$, and let $M$ be the Mach number associated to the incident shock, defined by $M = \sigma/c_R$. Then $W$ is given by:

$$\rho = \frac{\gamma + 1}{\gamma - 1 + \frac{2}{M^2}} \rho_R, \quad (31a)$$

$$u = \left(1 - \frac{\rho_R}{\rho}\right) \sigma, \quad (31b)$$

$$p = p_R + \left(1 - \frac{\rho_R}{\rho}\right) \rho_R \sigma^2. \quad (31c)$$

*Proof.* We first change the coordinate system, in such a way that the shock is steady in the new coordinate system. The density and the pressure are left unchanged, while the velocity is now $u - \sigma$ and $-\sigma$ in the left and right state respectively. In this coordinate system, the Rankine-Hugoniot conditions imply that the jump of the fluxes vanishes, which reads for the Euler equations:

$$\rho(u - \sigma) = \rho_R(-\sigma), \quad (32a)$$

$$\rho(u - \sigma)^2 + p = \rho_R(-\sigma)^2 + p_R, \quad (32b)$$

$$\frac{1}{2} \rho (u - \sigma)^3 + \rho e (u - \sigma) + p (u - \sigma) = \frac{1}{2} \rho_R (-\sigma)^3 + \rho_R e_R (-\sigma) + p_R (-\sigma). \quad (32c)$$

This system must be complemented by the equation of state $p = (\gamma - 1)\rho e$ and $p_R = (\gamma - 1)\rho_R e_R$. Thanks to this relation, we may recast (32c) as:

$$\rho (\sigma - u)^3 + \xi p (\sigma - u) = \rho_R \sigma^3 + \xi \sigma p_R, \quad \xi = \frac{2\gamma}{\gamma - 1}. \quad (33)$$

Relation (32b) reads

$$\frac{1}{\rho} [\rho (u - \sigma)]^2 + p = \rho_R \sigma^2 + p_R,$$
so, using (32a):

$$ \frac{1}{\rho}(\rho_R \sigma)^2 + p = \rho_R \sigma^2 + p_R. $$

We thus obtain $p$ as a function of known quantities (i.e. $\sigma$ and the right state) and $\rho$ only:

$$ p = p_R + \left(1 - \frac{\rho_R}{\rho}\right) \rho_R \sigma^2. \qquad (34) $$

We now notice that substituting, in the jump condition associated to the energy balance (33), this expression for $p$ and $-(\rho_R/\rho)\sigma$ for $(u-\sigma)$, thanks once again to Equation (32a), we get an equation for $\rho$ only:

$$ \rho \left(\frac{\rho_R}{\rho} \sigma\right)^3 + \xi \left(p_R + \left(1 - \frac{\rho_R}{\rho}\right) \rho_R \sigma^2\right) \left(\frac{\rho_R}{\rho} \sigma\right) = \rho_R \sigma^3 + \xi \sigma p_R. $$

If $\sigma = 0$, the first jump condition implies $u = 0$ (excluding $\rho = 0$), then the second one yields $p = p_R$ and the third one is automatically satisfied: the considered discontinuity is a (stationary) contact. In such a case, the left state remains partially undetermined by the jump conditions: $\rho$ and $e$ may take any value satisfying $(\gamma - 1)\rho e = p_R$. Since we only consider the case of a shock, $\sigma \neq 0$ and the last relation may be simplified by $\sigma$. Reordering, we get:

$$ \rho_R \sigma^2 \left(\left(\frac{\rho_R}{\rho}\right)^2 - 1\right) + \xi p_R \left(\frac{\rho_R}{\rho} - 1\right) - \xi \rho_R \sigma^2 \frac{\rho_R}{\rho} \left(\frac{\rho_R}{\rho} - 1\right) = 0. $$

The case $\rho = \rho_R$ has no interest: it yields, by the first jump condition, $u = 0$, and the second one implies that $p = p_R$, which means in fine that $W = W_R$, i.e. that there is no discontinuity at all. We may thus simplify by $1 - \rho_R/\rho$, to obtain a linear equation for the ratio $\rho_R/\rho$. Solving this latter equation, we obtain:

$$ \rho = \frac{\gamma + 1}{\gamma - 1 + \frac{2}{M^2}} \rho_R, \quad \text{with } M = \frac{\sigma}{c_R}, \quad c_R^2 = \gamma \frac{p_R}{\rho_R}. $$

We thus obtain (31a). Relation (31b) is a straightforward consequence of (32a) and (31c) was already proven (see Relation (34)). The proof is thus complete. $\square$

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

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

For a 3-shock, the Lax entropy condition requires that $\sigma > c_R$, which is equivalent to $M > 1$. We thus have $\rho > \rho_R$, $p > p_R$ and $0 < u < \sigma$.

It is worth noting that it is now easy to relax the assumption $u_R = 0$ of Lemma B.1. Indeed, for the general case, we may work in the system of coordinates in translation at the velocity $u_R$ with respect to the initial one: in the new coordinate system, the right state is now at rest and Lemma B.1 applies, replacing $\sigma$ by $\sigma - u_R$ and $u$ by $u - u_R$ in the definition of the Mach number $M$ and in the system of equations (31).

## APPENDIX C. A RELATION SATISFIED BY THE RIGHT STATE OF A SHOCK WHEN THE LEFT STATE IS AT REST

In this section, we perform a technical computation motivated by the following arguments. Let us suppose that a shock travelling at the speed $\sigma_r$ separates a left state denoted by $W_b = (\rho_b, u_b, p_b)$ and a right state denoted by $W_2 = (\rho_2, u_2, p_2)$, and that $u_b = 0$. The Rankine-Hugoniot conditions yield three independent equations, and thus constitute a system in which $\rho_b$ and $p_b$ may be eliminated, to obtain a relation linking $W_2$ and $\sigma_r$ only. We now derive this relation, supposing that the following specific constitutive relations hold for both states:

$$ E_b = \frac{1}{2} u_b^2 + e_b - Q, \quad p_b = (\gamma_b - 1) \rho_b e_b, \quad E_2 = \frac{1}{2} u_2^2 + e_2, \quad p_2 = (\gamma_u - 1) \rho_2 e_2. \qquad (35) $$

Note that eliminating $\rho_b$ and $p_b$ consists in establishing an expression of these quantities (and thus of $W_b$, since $u_b = 0$) as a function of $W_2$ and $\sigma_r$. All of these relations, i.e. the equation linking $W_2$ and $\sigma_r$ and the expression of $W_b$ as a function of these variables, are gathered in the following lemma.
**Lemma C.1** (Some conditions at the reactive shock). The state $W_2$ and the shock speed $\sigma_r$ satisfy the following relation

$$ \frac{1}{2}u_2^2 + \frac{1}{\gamma_b - 1}u_2\sigma_r + \left(\frac{\gamma_u}{\gamma_u - 1} - \frac{\gamma_b}{\gamma_b - 1}\frac{\sigma_r}{\sigma_r - u_2}\right)\frac{p_2}{\rho_2} + Q = 0. \quad (36) $$

In addition, $W_b$ is given as a function of $W_2$ and $\sigma_r$ by:

$$ \rho_b = \rho_2\left(\frac{\sigma_r - u_2}{\sigma_r}\right), \quad p_b = p_2 - \rho_2 u_2(\sigma_r - u_2). \quad (37) $$

*Proof.* Using, as in the previous section, the standard change of coordinates to work in the coordinate system in which the shock is at rest, the Rankine-Hugoniot relationships (which boil down to the fact that the jump of the fluxes vanishes) through the shock for the mass and momentum balance equations read respectively:

$$ \rho_b \sigma_r = \rho_2 (\sigma_r - u_2), \quad \rho_b \sigma_r^2 + p_b = \rho_2 (\sigma_r - u_2)^2 + p_2. $$

These two relations readily yield the expression (37) of $W_b$ as a function of $W_2$ and $\sigma_r$ that we are looking for:

$$ \rho_b = \rho_2 \left( \frac{\sigma_r - u_2}{\sigma_r} \right) \quad \text{and} \quad p_b = p_2 - \rho_2 u_2 (\sigma_r - u_2). $$

Let us now write the Rankine-Hugoniot condition for the conservation equation of the total energy:

$$ \rho_b \sigma_r E_b + \sigma_r p_b = \rho_2 (\sigma_r - u_2) E_2 + (\sigma_r - u_2) p_2. \quad (38) $$

We may divide the left-hand side of this relation by $\rho_b\sigma_r$ and the right-hand side by $\rho_2(\sigma_r - u_2)$ (since these two expressions are equal by the jump condition associated to the mass balance equation), to obtain:

$$ E_b + \frac{p_b}{\rho_b} = E_2 + \frac{p_2}{\rho_2}, $$

that is, using the constitutive relations (35):

$$ \frac{1}{2}\sigma_r^2 + \frac{\gamma_b}{\gamma_b - 1}\frac{p_b}{\rho_b} - Q = \frac{1}{2}(\sigma_r - u_2)^2 + \frac{\gamma_u}{\gamma_u - 1}\frac{p_2}{\rho_2}. $$

Using expression (37) of $\rho_b$ and $p_b$ in this equation yields

$$ \frac{1}{2}\sigma_r^2 + \frac{\gamma_b}{\gamma_b - 1}\frac{p_2 - \rho_2 u_2(\sigma_r - u_2)}{\rho_2\frac{\sigma_r - u_2}{\sigma_r}} - Q = \frac{1}{2}(\sigma_r - u_2)^2 + \frac{\gamma_u}{\gamma_u - 1}\frac{p_2}{\rho_2}, $$

and thus, splitting the fraction at the left-hand side and reordering the terms:

$$ \frac{1}{2}(\sigma_r - u_2)^2 - \frac{1}{2}\sigma_r^2 + \frac{\gamma_b}{\gamma_b - 1}\sigma_r u_2 + \frac{\gamma_u}{\gamma_u - 1}\frac{p_2}{\rho_2} - \frac{\gamma_b}{\gamma_b - 1}\frac{\sigma_r}{\sigma_r - u_2}\frac{p_2}{\rho_2} + Q = 0, $$

which yields (36) and thus concludes the proof.

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

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

## APPENDIX D. VARIATIONS OF $\mathcal{F}_r$ WHEN $\gamma_b = \gamma_u$

In this appendix, we prove that the function $\mathcal{F}_r$ defined by (39) is increasing over $(\sigma_\ell, \sigma_p)$ if the heat capacity ratios $\gamma_b$ and $\gamma_u$ are equal. Let us denote by $\gamma$ the common heat capacity ratio, and recall the expression of $\mathcal{F}_r$:

$$ \mathcal{F}_r(x) = \frac{1}{2}u(x)^2 + \frac{1}{\gamma-1}xu(x) - \frac{\gamma}{\gamma-1}\frac{u(x)}{x-u(x)}\frac{p(x)}{\rho(x)} + Q. \quad (39) $$
Using the fact that $γp/ρ = c^2$ and the relation (10), we have

$$ \mathcal{F}'_r(x) = u(x)u'(x) + \frac{x}{\gamma-1}u'(x) + \frac{1}{\gamma-1}u(x) - \frac{1}{\gamma-1} \frac{xu'(x)-u(x)}{(x-u(x))^2}c^2(x) - \frac{1}{\gamma-1} \frac{u(x)}{x-u(x)}(c^2)'(x), $$

with

$$ (c^2)' = s_1\gamma(\rho^{\gamma-1})' = s_1\gamma(\gamma-1)\rho^{\gamma-2}\rho' = (\gamma-1)c^2\frac{\rho'}{\rho}, $$

where $s_1$ is the value of the (constant) entropy in the intermediate zone. Hence we get that $\mathcal{F}'_r(x) = T_1(x) + T_2(x) + T_3(x)$ with

$$ \begin{aligned} T_1(x) &= \frac{1}{\gamma-1} \left( 1 + \frac{c^2(x)}{(x-u(x))^2} \right) u(x), \\ T_2(x) &= \left( u(x) + \frac{x}{\gamma-1} - \frac{1}{\gamma-1} \frac{x c^2(x)}{(x-u(x))^2} \right) u'(x), \\ T_3(x) &= - \frac{u(x) c^2(x)}{(x-u(x))\rho} \rho'. \end{aligned} $$

By Lemma 2.2, $u > 0$ so that $T_1 > 0$. In addition, in the proof of the same Lemma 2.2, we showed that $u(x) + c(x) - x > 0$, so that $c(x) > x - u(x)$ and

$$ u(x) + \frac{x}{\gamma - 1} \left( 1 - \frac{c^2(x)}{(x - u(x))^2} \right) \le u(x). $$

Since $u' \le 0$, this implies that $T_2 \ge uu'$. Replacing $\rho'$ and $u'$ by their expressions given in (11), we obtain that $uu' + T_3 = 0$, which concludes the proof.

## REFERENCES

[1] A. Beccantini and E. Studer. The reactive Riemann problem for thermally perfect gases at all combustion regimes. *International Journal for Numerical Methods in Fluids*, 64:269–313, 2010.

[2] CALIF³S. A software components library for the computation of reactive turbulent flows. https://gforge.irsn.fr/gf/project/calif3s.

[3] P. Cambray and B. Deshaies. Ecoulement engendré par un piston sphérique : solution analytique approchée. *Acta Astronautica*, 5:611–617, 1977.

[4] B. Deshaies and J.C. Leyer. Flow induced by unconfined spherical accelerating flames. *Combustion and Flame*, 40:141–153, 1981.

[5] E. Godlewski and P.-A. Raviart. *Numerical approximation of hyperbolic systems of conservation laws*. Springer, 1996.

[6] C.M. Guirano, G.G. Bash, and John H. Lee. Pressure waves generated by spherical flames. *Combustion and Flame*, 27:341–351, 1976.

[7] A.L. Kuhl, M.M. Kamel, and A.K. Oppenheim. Pressure waves generated by steady flames. *Symposium (International) on Combustion*, 14:1201–1215, 1973.

[8] P²REMICS. A simulation tool for explosion hazards. https://gforge.irsn.fr/gf/project/calif3s.

[9] L.I. Sedov. On certain unsteady motions of compressible fluid. *Prikladnaya Matematika i Mekhanika*, 9:293–311, 1945.

[10] R.A. Strehlow, R.T. Luckritz, A.A. Adamczyk, and S.A. Shimpi. The blast wave generated by spherical flmaes. *Combustion and Flame*, 35:297–310, 1979.

[11] G.I. Taylor. The air wave surrounding an expanding sphere. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences*, 186:273–292, 1946.