# Positivity-preserving and entropy-bounded discontinuous Galerkin method for
the chemically reacting, compressible Euler equations. Part I: The one-
dimensional case

Eric J. Ching, Ryan F. Johnson, and Andrew D. Kercher  Laboratories for
Computational Physics and Fluid Dynamics, U.S. Naval Research Laboratory, 4555
Overlook Ave SW, Washington, DC 20375

###### Abstract

In this paper, we develop a fully conservative, positivity-preserving, and
entropy-bounded discontinuous Galerkin scheme for simulating the
multicomponent, chemically reacting, compressible Euler equations with complex
thermodynamics. The proposed formulation is an extension of the fully
conservative, high-order numerical method previously developed by Johnson and
Kercher [_J. Comput. Phys._ , 423 (2020), 109826] that maintains pressure
equilibrium between adjacent elements. In this first part of our two-part
paper, we focus on the one-dimensional case. Our methodology is rooted in the
minimum entropy principle satisfied by entropy solutions to the
multicomponent, compressible Euler equations, which was proved by Gouasmi et
al. [_ESAIM: Math. Model. Numer. Anal._ , 54 (2020), 373–389] for nonreacting
flows. We first show that the minimum entropy principle holds in the reacting
case as well. Next, we introduce the ingredients, including a simple linear-
scaling limiter, required for the discrete solution to have nonnegative
species concentrations, positive density, positive pressure, and bounded
entropy. We also discuss how to retain the aforementioned ability to preserve
pressure equilibrium between elements. Operator splitting is employed to
handle stiff chemical reactions. To guarantee discrete satisfaction of the
minimum entropy principle in the reaction step, we develop an entropy-stable
discontinuous Galerkin method based on diagonal-norm summation-by-parts
operators for solving ordinary differential equations. The developed
formulation is used to compute canonical one-dimensional test cases, namely
thermal-bubble advection, multicomponent shock-tube flow, and a moving
hydrogen-oxygen detonation wave with detailed chemistry. We demonstrate that
the developed formulation can achieve optimal high-order convergence in smooth
flows. Furthermore, we find that the enforcement of an entropy bound can
considerably reduce the large-scale nonlinear instabilities that emerge when
only the positivity property is enforced, to an even greater extent than in
the monocomponent, calorically perfect case. Finally, mass, total energy, and
atomic elements are shown to be discretely conserved.

###### keywords:

Discontinuous Galerkin method; Combustion; Detonation; Minimum entropy
principle; Positivity-preserving; Entropy stability; Summation-by-parts

††footnotetext:
DISTRIBUTION STATEMENT A. Approved for public release. Distribution is
unlimited.

##  1 Introduction

The discontinuous Galerkin (DG) method [1, 2, 3, 4, 5] has recently gained
considerable attention in the computational fluid dynamics community [6].
Several desirable properties, such as local conservation, arbitrarily high
order of accuracy on unstructured grids, and suitability for heterogeneous
computing systems, highlight its great potential to accurately and efficiently
simulate complex fluid flows. However, it is well-known that nonlinear
instabilities are easily introduced in underresolved regions and near flow-
field discontinuities. This issue is exacerbated when realistic thermodynamics
(e.g., the thermally perfect gas model) and multispecies chemical reactions
are incorporated. For example, fully conservative schemes (not just the DG
method) are known to generate spurious pressure oscillations in moving
interface problems [7, 8, 9]. To remedy this issue, quasi-conservative methods
are often employed, such as the double-flux scheme [10], in which the equation
of state is recast based on a calorically perfect gas model using frozen,
elementwise-constant auxiliary variables. The double-flux method has been
utilized in a number of studies involving DG simulations of multicomponent
flows [11, 12, 13]. While effective at eliminating the aforementioned pressure
oscillations, the double-flux approach violates energy conservation, which can
be crucial for the prediction of shock speeds and locations, as well as heat
release in combustion processes. As a compromise, Lv and Ihme [12] proposed a
hybrid double-flux strategy wherein a fully conservative method is employed at
shocks.

On the other hand, Johnson and Kercher recently proposed an explicit, fully
conservative, high-order method that does not generate spurious pressure
oscillations in smooth flow regions or across material interfaces when the
temperature is continuous [14]. This is done via (a) exact evaluation of the
thermodynamics and (b) calculation of the inviscid and viscous fluxes in a
consistent manner that maintains pressure equilibrium, which is drastically
simplified through a particular choice of nodal basis. It is worth noting that
their proposed strategy is not limited to DG schemes, but can be applied to
other numerical methods as well. Stiff chemical reactions were handled via
operator splitting. Efficient and accurate integration of the chemical source
terms was achieved via an $hp$-adaptive DG method for solving ordinary
differential equations, termed _DGODE_. Optimal high-order accuracy was
demonstrated for smooth flows, and a suite of complex multicomponent reacting
flows was computed. The high-order calculation of a three-dimensional reacting
shear flow in the presence of a splitter plate did not require any additional
stabilization. Also computed was a two-dimensional, moving detonation wave.
Artificial viscosity was used to stabilize the shock fronts present in the
solution. However, a very fine mesh was required to maintain robustness while
achieving correct prediction of the cellular structures behind the shock,
highlighting the difficulty of robustly and accurately simulating
multidimensional detonation waves on coarse meshes, especially for high-order
methods. Even if spurious pressure oscillations are sufficiently minimized,
the wide range of complex flow features characterizing detonations is
difficult to capture [15]. Such features include thin reaction zones,
traveling pressure waves, shock-shock interactions, Kelvin-Helmholtz
instabilities, vortical structures, and triple points. As previously
discussed, underresolution of flow features induces instabilities that can
cause solver divergence, and these instabilities may be amplified by the added
nonlinearity of the variable thermodynamics, multicomponent flow, and stiff
chemical reactions. Another difficulty associated with multicomponent flow is
the frequent occurrence of negative species concentrations, especially since
initial and boundary conditions often specify the mole fractions of certain
species to be zero. Many reacting-flow solvers simply “clip” negative
concentrations to zero, which violates conservation and introduces low-order
errors.

In light of the above, our primary objective in this study is to develop a
positivity-preserving and entropy-bounded DG scheme for simulating the
multicomponent, chemically reacting Euler equations with exact thermodynamics
for mixtures of thermally perfect gases. Specifically, we build upon the
aforementioned fully conservative high-order method that can maintain pressure
equilibrium [14]. In this first part of our two-part paper, we focus on the
one-dimensional case. The unique challenges posed by realistic thermodynamics
and stiff chemical source terms are discussed and addressed. _Entropy-bounded_
in this context means that the specific thermodynamic entropy of the discrete
solution is bounded from below, an idea rooted in the minimum entropy
principle satisfied by entropy solutions to the compressible, multicomponent
Euler equations. This principle was recently proved for the nonreacting case
by Gouasmi et al. [16], which is an important prerequisite of this work. The
developed formulation in general preserves order of accuracy for smooth
solutions. Our main contributions are as follows:

  * 1.

A minimum entropy principle for the compressible, multicomponent, chemically
reacting Euler equations is demonstrated, which follows naturally from the
proof in [16].

  * 2.

In a DG framework, we extend the fully conservative high-order discretization
in [14] to be positivity-preserving and entropy-bounded. We also discuss a
different local entropy bound that is less restrictive than that previously
introduced in [17] in the context of the monocomponent Euler equations.

  * 3.

To provably guarantee satisfaction of the minimum entropy principle in the
temporal integration of stiff chemical source terms, we extend DGODE by
developing an entropy-stable DGODE based on diagonal-norm summation-by-parts
(SBP) operators. This involves deriving a new entropy-conservative two-point
numerical state function (note that similar entropy-conservative numerical
state functions were derived for the monocomponent Euler, shallow-water, and
ideal magnetohydrodynamics equations by Friedrich et al. [18] in the context
of an entropy-stable space-time DG discretization).

  * 4.

We employ the proposed entropy-bounded DG method to robustly and accurately
compute a series of canonical one-dimensional test cases. Mass, energy, and
atom conservation are maintained. The proposed formulation more effectively
suppresses spurious oscillations than the positivity-preserving DG scheme. In
particular, we find that the relative benefit of enforcing an entropy bound is
significantly greater in the multicomponent, thermally perfect setting than in
monocomponent, calorically perfect setting.

In Part II [19], we extend the entropy-bounded DG scheme to multiple
dimensions on arbitrary elements. Our multidimensional extension is a further
generalization of the multidimensional positivity-preserving/entropy-bounded
schemes currently in the literature [20, 21, 17, 22]. Specifically,
restrictions on the numerical flux, physical modeling, element shape,
polynomial order of the geometric approximation, and/or quadrature rules are
relaxed. Complex detonation waves in two and three dimensions are computed. We
find that whereas the positivity-preserving DG scheme is often not
sufficiently stable (even with artificial viscosity to stabilize strong
discontinuities), enforcing an entropy bound enables robust calculations on
relatively coarse meshes.

###  1.1 Background

Various stabilization strategies for high-order DG schemes have been
introduced in the literature. Artificial viscosity is a popular approach that
very effectively suppresses oscillations and is perfectly compatible with high
polynomial orders, arbitrary elements, and general equation sets [23, 24, 25].
However, there are certain limitations that discourage an overreliance on
artificial viscosity for suppressing _all_ instabilities. First, it can
significantly pollute accuracy, especially in smooth regions of the flow. As
such, it should ideally be added only where necessary (e.g., strong shocks
that are otherwise difficult to robustly capture). Second, design of a shock
sensor that can reliably detect discontinuous flow features for general
configurations remains an open problem. Even in a single flow configuration
that involves discontinuities of varying strengths (as is the case for
detonations), it is difficult to detect all such discontinuities and add the
“right” amount of artificial viscosity. Third, there is typically a very
strong dependence on tunable parameters.

A common alternative to artificial viscosity is limiting. WENO-type [26, 27],
TVD/TVB [28, 29], and moment limiters [30] are well-known examples. However,
it can be difficult to extend these limiters to arbitrary polynomial orders
for both the solution and geometric approximations, as well as to general
equation sets. Furthermore, these limiters are not guaranteed to yield
physically admissible solutions (i.e., positive density and pressure), a
drawback of artificial viscosity as well. An alternative approach for
maintaining robustness is to align the grid with discontinuities. Recent
formulations that do so in an implicit manner by treating the grid as a
variable were developed by Corrigan et al. [31] and Zahr and Persson [32]. An
encouraging preliminary effort to apply implicit shock tracking to reacting
flow is discussed in [33], in which supersonic inviscid reacting flow over a
two-dimensional wedge with simple thermodynamics and chemistry was computed.

Recently, positivity-preserving DG schemes [20, 34, 21, 35] have seen
considerable success in solving the monocomponent, nonreacting Euler equations
with explicit time stepping. These schemes prevent the occurrence of negative
densities and pressures under a constraint on the time step size, a
positivity-preserving numerical flux, and a limiting procedure consisting of a
simple linear “squeezing” of the solution towards its cell average. The
limiting operator is conservative and maintains order of accuracy for smooth
solutions. However, the limiter is not very effective at dampening
oscillations. The positivity-preserving DG method was extended to the entropy-
bounded DG scheme by Zhang and Shu [36], in which an additional limiting step
based on a Newton search was introduced to enforce a global entropy bound..
Note that it is implicitly assumed that an entropy-bounded scheme is also
positivity-preserving. Under their invariant-region-preserving DG framework,
Jiang and Liu [36, 22] introduced a more straightforward limiter to enforce
the entropy principle in an algebraic manner, which is particularly desirable
in the multicomponent, thermally perfect case due to the cost of evaluating
complex thermodynamics. Lv and Ihme [17] further extended the entropy-bounded
DG scheme (for the monocomponent Euler equations) by relaxing restrictions on
the geometry and quadrature rules. They also introduced a local entropy bound.
Numerical tests demonstrated the superiority of the entropy-bounded DG scheme
for suppressing spurious oscillations, compared to the positivity-preserving
DG method. Lv and Ihme later applied entropy bounding to their reacting flow
DG solver [37], which utilizes the double-flux approach discussed above. The
frozen thermodynamics and relaxation towards a calorically perfect gas model
circumvent the difficulties of extending the entropy-bounded DG method to
reacting flow with exact thermodynamics and non-calorically-perfect gases.
Furthermore, the physico-mathematical validity of combining the double-flux
model with enforcement of a discrete minimum entropy principle is not
immediately clear. We also note that these positivity-preserving and entropy-
bounded DG methods are related to the recent geometric quasilinearization
framework by Wu and Shu [38], as well as the invariant-domain-preserving
schemes based on graph viscosity by, for example, Guermond et al. [39] and
Pazner [40], which employ a convex limiting procedure relying on an iterative
line search.

Before concluding this section, we note several other efforts to compute
chemically reacting flows with DG schemes in a stable manner. Gutierrez-
Jorquera and Kummer [41] computed steady-state diffusion flames using a low-
Mach pressure-based solver and a one-step reaction mechanism. May et al. [42]
simulated steady hypersonic flows with some high-enthalpy effects using a
hybridized DG solver that employs artificial viscosity for shock capturing.
Papoutsakis et al. [43] computed chemically reacting hypersonic flow over a
double cone with a TVB limiter; however, only a linear polynomial
approximation of the solution was employed, and there were discrepancies with
finite-volume and experimental results. A number of positivity-preserving DG
schemes for the reacting Euler equations have also been developed. For
example, Zhang and Shu [34] and Wang et al. [44] extended the positivity-
preserving DG scheme to handle source terms with explicit time stepping.
Complex thermodynamics and stiff source terms were not addressed. Du and Yang
[45] and Du et al. [46] presented a positivity-preserving DG method based on a
new explicit, exponential Runge-Kutta (RK)/multistep time integration scheme
that can handle stiff source terms [47]. However, the temporal order of
accuracy is strongly dependent on the initial conditions. Other positivity-
preserving time integrators compatible with stiff source terms include
implicit Patankar-type RK schemes [48, 49, 50], which have been applied to
finite difference [48, 49] and finite volume [50] discretizations of the
reacting Euler equations. These schemes are still undergoing development, and
there may be issues with extending them to DG discretizations [48]. An
additional obstacle of both exponential RK/multistep and Patankar-type RK
schemes is proper enforcement of the minimum entropy principle. For these
reasons and due to its proven past success, we elect to still employ operator
splitting to deal with stiff source terms. Nevertheless, exponential
RK/multistep and Pantankar-type schemes are indeed worthy of future
consideration.

That we use an entropy-stable DG discretization to temporally integrate the
chemical source terms may cause readers to question why an entropy-stable DG
discretization is not also employed for the transport step. Note that entropy-
stable schemes guarantee the global integral of mathematical entropy to be
nonincreasing in time (assuming periodic or entropy-stable boundary
conditions), whereas entropy-bounded schemes enforce, in a pointwise fashion,
the specific thermodynamic entropy to be greater than some lower bound.
Entropy-stable DG schemes for multicomponent (nonreacting) flows have emerged
only recently [51, 52]. For general hyperbolic conservation laws, there are
various ways to achieve entropy stability; two of the most well-known are as
follows: (a) the entropy variables, instead of the conservative variables, are
directly solved for, and (b) SBP operators are used to approximate discrete
derivatives. The first approach, in which the conservative variables depend
implicitly on the polynomial approximation of the entropy variables, is not
appropriate for explicit time stepping, which we employ in this work (for the
transport terms). Furthermore, with the typical choice of entropy, the entropy
variables are undefined when any of the partial densities vanishes [51].
Conversely, the second approach, which has surged in popularity in recent
years, is compatible with explicit time stepping and does not assume exact
integration to provably guarantee entropy stability. Belonging to this
category are the methods by Renac [52] and Peyvan et al. [53]. On unstructured
grids, particularly those with (possibly curved) simplicial elements, current
SBP-based entropy-stable methods can become suboptimal [54], require direct
use of the entropy variables (even if the conservative variables are the
unknowns) [55, 56], and/or significantly increase in complexity [55, 57, 56,
58]. Note that these issues are not present in the developed entropy-stable
DGODE (except a moderate increase in complexity), which entails a one-
dimensional discretization in time. Due to the above factors, as well as the
relative simplicity of the entropy-bounded DG scheme and its compatibility
with the pressure-equilibrium-maintaining discretization by Johnson and
Kercher [14], we choose to employ the proposed entropy-bounded method for the
transport step. Nevertheless, we emphasize that we are not ruling out entropy-
stable schemes; advancements to these formulations are rapid and their
potential is evident [59]. Finally, it should be noted that although the
construction of entropy-stable and entropy-bounded methods rely on different
techniques, entropy stability and entropy boundedness are not necessarily
mutually exclusive. Depending on the type of entropy-stable scheme, it should,
in principle, be feasible to achieve discrete satisfaction of both properties.

The remainder of this paper is organized as follows. Sections 2 and 3
summarize the governing equations and basic DG discretization, respectively.
Section 4 reviews the minimum entropy principle associated with the
compressible, multicomponent, nonreacting Euler equations and extends it to
the reacting Euler equations. The next section presents the positivity-
preserving and entropy-bounded DG formulation for the transport step. We then
discuss the entropy-stable DG discretization for the reaction step in Section
6. Results for fundamental nonreacting test cases and one-dimensional
detonation-wave simulations are given in the following section. We close the
paper with concluding remarks.

##  2 Governing equations

The compressible, multicomponent, chemically reacting Euler equations are
given as

$\frac{\partial y}{\partial
t}+\nabla\cdot\mathcal{F}\left(y\right)-\mathcal{S}\left(y\right)=0$ (2.1)

where $t\in\mathbb{R}^{+}$ is time,
$y(x,t):\mathbb{R}^{d}\times\mathbb{R}^{+}\rightarrow\mathbb{R}^{m\times d}$
is the conservative state vector (with $x=(x_{1},\ldots,x_{d})$ denoting the
physical coordinates),
$\mathcal{F}(y):\mathbb{R}^{m}\rightarrow\mathbb{R}^{m\times d}$ is the
convective flux, $\mathcal{S}(y):\mathbb{R}^{m}\rightarrow\mathbb{R}^{m}$ is
the chemical source term. The state vector is expanded as

$y=\left(\rho v_{1},\ldots,\rho v_{d},\rho
e_{t},C_{1},\ldots,C_{n_{s}}\right)^{T},$ (2.2)

where $n_{s}$ is the number of species (which yields $m=d+n_{s}+1$), $\rho$ is
density, $v=\left(v_{1},\ldots,v_{d}\right)$ is the velocity, $e_{t}$ is the
mass-specific total energy, and $C=\left(C_{1},\ldots,C_{n_{s}}\right)$ are
the species concentrations. The density is computed from the species
concentrations as

$\rho=\sum_{i=1}^{n_{s}}\rho_{i}=\sum_{i=1}^{n_{s}}W_{i}C_{i},$

where $\rho_{i}$ is the partial density and $W_{i}$ is the molecular weight of
the $i$th species. The mass fraction and mole fraction of the $i$th species
are defined as

$Y_{i}=\frac{\rho_{i}}{\rho}$

and

$X_{i}=\frac{C_{i}}{\sum_{i=1}^{n_{s}}C_{i}},$

respectively.

The $k$th spatial convective flux component is written as

$\mathcal{F}_{k}^{c}\left(y\right)=\left(\rho
v_{k}v_{1}+P\delta_{k1},\ldots,\rho v_{k}v_{d}+P\delta_{kd},v_{k}\left(\rho
e_{t}+P\right),v_{k}C_{1},\ldots,v_{k}C_{n_{s}}\right)^{T},$ (2.3)

where $P$ is the pressure. The mass-specific total energy is the sum of the
specific internal and kinetic energies, given by

$e_{t}=u+\frac{1}{2}\sum_{k=1}^{d}v_{k}v_{k},$

where the (mixture-averaged) mass-specific internal energy, $u$, is the mass-
weighted sum of the mass-specific internal energies of each species:

$u=\sum_{i=1}^{n_{s}}Y_{i}u_{i}.$

This work assumes thermally perfect gases, with $u_{i}$ given by [60]

$u_{i}=h_{i}-R_{i}T=h_{\mathrm{ref},i}+\int_{T_{\mathrm{ref}}}^{T}c_{p,i}(\tau)d\tau-
R_{i}T,$

where $h_{i}$ is the mass-specific enthalpy of the $i$th species,
$R_{i}=R^{0}/W_{i}$ (with
$R^{0}=8314.4621\,\mathrm{JKmol}^{-1}\mathrm{K}^{-1}$ denoting the universal
gas constant), $T$ is the temperature, $T_{\mathrm{ref}}$ is the reference
temperature (298.15 K), $h_{\mathrm{ref},i}$ is the reference-state species
formation enthalpy, and $c_{p,i}$ is the mass-specific heat at constant
pressure of the $i$th species. $c_{p,i}$ is computed from an $n_{p}$-order
polynomial as

$c_{p,i}=\sum_{k=0}^{n_{p}}a_{ik}T^{k},$ (2.4)

based on the NASA coefficients [61, 62]. The mass-specific thermodynamic
entropy of the mixture is defined as

$s=\sum_{i=1}^{n_{s}}Y_{i}s_{i},$

with $s_{i}$ given by

$s_{i}=s_{i}^{o}-R_{i}\log\frac{P_{i}}{P_{\mathrm{ref}}},\quad
s_{i}^{o}=s_{\mathrm{ref},i}^{o}+\int_{T_{\mathrm{ref}}}^{T}\frac{c_{p,i}(\tau)}{\tau}d\tau,$

where $s_{\mathrm{ref},i}^{o}$ is the species formation entropy at the
reference temperature and reference pressure ($P_{\mathrm{ref}}=1\text{
atm}$), $s_{i}^{o}$ denotes the species entropy at atmospheric pressure, and
$P_{i}=C_{i}R^{0}T$ is the partial pressure. $s_{i}$ can also be expressed as
[60, 16, 51]

$s_{i}=s_{\mathrm{ref},i}^{o}+\int_{T_{\mathrm{ref}}}^{T}\frac{c_{v,i}(\tau)}{\tau}d\tau-
R_{i}\log\frac{C_{i}}{C_{\mathrm{ref}}},$

where $C_{\mathrm{ref}}=P_{\mathrm{ref}}/R^{0}T_{\mathrm{ref}}$ is the
reference concentration and $c_{v,i}=c_{p,i}-R_{i}$ is the mass-specific heat
at constant volume of the $i$th species. Summing up the partial pressures
yields the equation of state for the mixture:

$P=R^{0}T\sum_{i=1}^{n_{s}}C_{i}.$ (2.5)

$u_{i}$, $h_{i}$, and $s_{i}^{o}$ are computed by integrating Equation (2.4)
and incorporating the integration constants in [61] and [62]. For example,
$u_{i}$ is calculated as

$u_{i}=b_{i0}+\sum_{k=0}^{n_{p}}\frac{a_{ik}}{k+1}T^{k+1}-R_{i}T=\sum_{k=0}^{n_{p}+1}b_{ik}T^{k},$
(2.6)

where $b_{i0}$ is the integration constant and

$b_{ik}=\begin{cases}\frac{a_{i,k-1}}{k},&k>1\\\
a_{i0}-R_{i},&k=1.\end{cases}$

###  2.1 Chemical reaction rates

The source term in Equation (2.1) is a smooth function of the state variables,
written as [63]

$\mathcal{S}\left(y\right)=\left(0,\ldots,0,0,\omega_{1},\ldots,\omega_{n_{s}}\right)^{T},$
(2.7)

where $\omega_{i}$ is the production rate of the $i$th species, which
satisfies mass conservation:

$\sum_{i=1}^{n_{s}}W_{i}\omega_{i}=0.$ (2.8)

The production rate is computed as

$\omega_{i}=\sum_{j=1}^{n_{r}}\nu_{ij}q_{j}.$

$n_{r}$ is the number of reactions, $\nu_{ij}=\nu_{ij}^{r}-\nu_{ij}^{f}$ is
the difference between the reverse ($\nu_{ij}^{r}$) and the forward
stoichiometric coefficients ($\nu_{ij}^{f}$), and $q_{j}$ is the rate of
progress of the $j$th reaction, computed as

$q_{j}=k_{j}^{f}\prod_{i=1}^{n_{s}}C_{i}^{\nu_{ij}^{f}}-k_{j}^{r}\prod_{i=1}^{n_{s}}C_{i}^{\nu_{ij}^{r}},$
(2.9)

where $k_{j}^{f}$ and $k_{j}^{r}$ are the forward and reverse rate constants,
respectively, of the $j$th reaction. The forward and reverse rate constants
are related via the equilibrium constant,

$K_{j}^{e}=\exp\left(-\frac{\Delta
G_{j}^{\prime}}{R^{0}T}\right)\left(\frac{P_{\mathrm{ref}}}{R^{0}T}\right)^{\sum_{i}\nu_{ij}},$
(2.10)

where $\Delta G_{j}^{\prime}$ is the change in reference-state Gibbs free
energy for the $j$th reaction, given as

$\Delta
G_{j}^{\prime}=\sum_{i=1}^{n_{s}}\nu_{ij}W_{i}h_{i}-T\sum_{i=1}^{n_{s}}\nu_{ij}W_{i}s_{i}^{\prime}.$

Introducing the reduced chemical potentials of the $i$th species,

$\displaystyle\mu_{i}$ $\displaystyle=\frac{g_{i}}{R^{0}T},$
$\displaystyle\mu_{i}^{\mathrm{u}}$ $\displaystyle=\mu_{i}-\frac{1}{W_{i}}\log
C_{i},$

where $g_{i}=h_{i}-Ts_{i}$ is the Gibbs function of the $i$th species, the
equilibrium constant can also be written as [60, Chapter 6.4]

$K_{j}^{e}=\exp\left(-\sum_{i=1}^{n_{s}}\nu_{ij}W_{i}\mu_{i}^{\mathrm{u}}\right),$
(2.11)

There exist various models for approximating the forward rate constants in
Equation (2.9), several of which will be briefly discussed next.

####  2.1.1 Arrhenius reactions

The Arrhenius form is the most common model for approximating reaction rates.
The forward rate constants are computed as

$k_{j}^{f}=A_{j}T^{b_{j}}\exp\left(-\frac{E_{j}}{R^{0}T}\right),$

where $A_{j}>0$ and $b_{j}$ are parameters and $E_{j}\geq 0$ is the activation
energy [60, 63].

####  2.1.2 Three-body reactions

These reactions require a “third body” in order to proceed. Dissociation and
recombination reactions are often of this type. The rate of progress is scaled
by a prefactor as [63]

$q_{j}=\left(\sum_{i=1}^{n_{s}}\alpha_{ij}C_{i}\right)\left(k_{j}^{f}\prod_{i=1}^{n_{s}}C_{i}^{\nu_{ij}^{f}}-k_{j}^{r}\prod_{i=1}^{n_{s}}C_{i}^{\nu_{ij}^{r}}\right),$

where $\alpha_{ij}$ are the third-body efficiencies.

####  2.1.3 Unimolecular/recombination fall-off reactions

Unimolecular/recombination fall-off reactions incorporate a dependence on
pressure. In general, this model predicts an increase in the reaction rate
with increasing pressure. For brevity, we drop the $j$ subscript and $f$
superscript. Given Arrhenius-type low-pressure and high-pressure limits for
the rate constant ($k_{0}$ and $k_{\infty}$, respectively), $k$ is computed as

$k=k_{\infty}\left(\frac{P_{r}}{1+P_{r}}\right)F,$ (2.12)

where $P_{r}$ is the reduced pressure, defined as

$P_{r}=\frac{k_{0}}{k_{\infty}}\sum_{i=1}^{n_{s}}\alpha_{i}C_{i}.$

There are different ways to compute $F$ in Equation (2.12). With the Lindemann
[64] approach, $F$ is simply unity. In the Troe [65] form, $F$ is given by

$\log F=\frac{\log F_{\mathrm{cent}}}{1+\left[\frac{\log
P_{r}+c_{1}}{c_{2}-c_{3}\left(\log P_{r}+c_{1}\right)}\right]^{2}},$

where the definitions of $c_{1}$, $c_{2}$, $c_{3}$, and $F_{\mathrm{cent}}$
are given in [65].

####  2.1.4 Chemically activated bimolecular reactions

Reactions of this type are also pressure-dependent, but the reaction rates
typically decrease with increasing pressure. The rate constants are computed
as [63]

$k=k_{0}\left(\frac{1}{1+P_{r}}\right)F,$

where $k_{0}$, $P_{r}$, and $F$ are calculated as in Section 2.1.3.

##  3 Discontinuous Galerkin discretization

In this section, we briefly describe the DG discretization of the governing
equations and review the techniques proposed in [14] to prevent spurious
pressure oscillations in smooth regions of the flow.

Let $\Omega\subset\mathbb{R}^{d}$ be the $d$-dimensional computational domain
partitioned by $\mathcal{T}$, which consists of non-overlapping cells $\kappa$
with boundaries $\partial\kappa$. Let $\mathcal{E}$ denote the set of
interfaces $\epsilon$, with
$\cup_{\epsilon\in\mathcal{E}}\epsilon=\cup_{\kappa\in\mathcal{T}}\partial\kappa$.
$\mathcal{E}$ consists of the interior interfaces,

$\epsilon_{\mathcal{I}}\in\mathcal{E_{I}}=\left\\{\epsilon_{\mathcal{I}}\in\mathcal{E}\,\middle|\,\epsilon_{\mathcal{I}}\cap\partial\Omega=\emptyset\right\\},$

and boundary interfaces,

$\epsilon_{\partial}\in\mathcal{E}_{\partial}=\left\\{\epsilon_{\partial}\in\mathcal{E}\,\middle|\,\epsilon_{\partial}\subset\partial\Omega\right\\},$

such that $\mathcal{E}=\mathcal{E_{I}}\cup\mathcal{E}_{\partial}$. At interior
interfaces, there exists $\kappa^{+},\kappa^{-}\in\mathcal{T}$ such that
$\epsilon_{\mathcal{I}}=\partial\kappa^{+}\cap\partial\kappa^{-}$. $n^{+}$ and
$n^{-}$ denote the outward facing normal of $\kappa^{+}$ and $\kappa^{-}$,
respectively, with $n^{+}=-n^{-}$. The discrete (finite-dimensional) subspace
$V_{h}^{p}$ over $\mathcal{T}$ is defined as

$\displaystyle V_{h}^{p}$ $\displaystyle=$
$\displaystyle\left\\{\mathfrak{v}\in\left[L^{2}\left(\Omega\right)\right]^{m}\,\middle|\,\forall\kappa\in\mathcal{T},\left.\mathfrak{v}\right|_{\kappa}\in\left[\mathcal{P}_{p}(\kappa)\right]^{m}\right\\},$
(3.1)

where, for $d=1$, $\mathcal{P}_{p}(\kappa)$ is the space of polynomial
functions of degree no greater than $p$ in $\kappa$. For $d>1$, the choice of
polynomial space typically depends on the element type [66].

The semi-discrete form of the governing equations (Equation (2.1)) is given
as: find $y\in V_{h}^{p}$ such that

$\displaystyle\sum_{\kappa\in\mathcal{T}}\left(\frac{\partial y}{\partial
t},\mathfrak{v}\right)_{\kappa}-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{F}\left(y\right),\nabla\mathfrak{v}\right)_{\kappa}+\sum_{\epsilon\in\mathcal{E}}\left(\mathcal{F}^{\dagger}\left(y^{+},y^{-},n\right),\left\llbracket\mathfrak{v}\right\rrbracket\right)_{\mathcal{E}}-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{S}\left(y\right),\mathfrak{v}\right)_{\kappa}=0\qquad\forall\>\mathfrak{v}\in
V_{h}^{p},$ (3.2)

where $\left(\cdot,\cdot\right)$ denotes the inner product,
$\mathcal{F}^{\dagger}\left(y^{+},y^{-},n\right)$ is the numerical flux, and
$\left\llbracket\cdot\right\rrbracket$ is the jump operator, given by
$\left\llbracket\mathfrak{v}\right\rrbracket=\mathfrak{v}^{+}-\mathfrak{v}^{-}$
at interior interfaces and
$\left\llbracket\mathfrak{v}\right\rrbracket=\mathfrak{v}^{+}$ at boundary
interfaces. Applying a standard, fully explicit time stepping scheme to
Equation (3.2) would yield an exceedingly small time step due to the stiff
chemical source terms. As such, operator splitting is employed to decouple the
temporal integration of the convection operator from that of the source term.
Specifically, we apply Strang splitting [67] over a given interval
$(t_{0},t_{0}+\Delta t]$ as

$\displaystyle\frac{\partial y}{\partial
t}+\nabla\cdot\mathcal{F}\left(y\right)=0$ $\displaystyle\textup{ in
}\Omega\times\left(t_{0},t_{0}+\nicefrac{{\Delta t}}{{2}}\right],$ (3.3)
$\displaystyle\frac{\partial y}{\partial t}-\mathcal{S}\left(y\right)=0$
$\displaystyle\textup{ in }\left(t_{0},t_{0}+\Delta t\right],$ (3.4)
$\displaystyle\frac{\partial y}{\partial
t}+\nabla\cdot\mathcal{F}\left(y\right)=0$ $\displaystyle\textup{ in
}\Omega\times\left(t_{0}+\nicefrac{{\Delta t}}{{2}},t_{0}+\Delta t\right],$
(3.5)

where Equations (3.3) and (3.5) are integrated in time with an explicit RK-
type scheme, while Equation (3.4) is solved using a fully implicit, temporal
DG discretization for ODEs (DGODE). Details on DGODE and its extension to
entropy-stable DGODE are given in Section 6. More sophisticated operator-
splitting schemes can be employed as well [68].

Unless otherwise specified, the volume and surface terms in Equation (3.2) are
evaluated using a quadrature-free approach [69, 70]. Throughout this work, we
employ a nodal basis, such that the element-local polynomial approximation of
the solution is expanded as

$y_{\kappa}=\sum_{j=1}^{n_{b}}y_{\kappa}(x_{j})\phi_{j},$ (3.6)

where $n_{b}$ is the number of basis functions,
$\left\\{\phi_{1},\ldots,\phi_{n_{b}}\right\\}$ are the basis functions, and
$\left\\{x_{1},\ldots,x_{n_{b}}\right\\}$ are the node coordinates. Its
average over $\kappa$ is given by

$\overline{y}_{\kappa}=\frac{1}{\left|\kappa\right|}\int_{\kappa}ydx,$ (3.7)

where $\left|\kappa\right|$ is the volume of $\kappa$. In the evaluation of
the second and third integrals in Equation (3.2), the nonlinear convective
flux can be approximated as

$\mathcal{F_{\kappa}}\approx\sum_{k=1}^{n_{c}}\mathcal{F}\left(y_{\kappa}\left(x_{k}\right)\right)\varphi_{k},$
(3.8)

where $n_{c}\geq n_{b}$ and
$\left\\{\varphi_{1},\ldots,\varphi_{n_{c}}\right\\}$ is a set of (potentially
different) polynomial basis functions. If $n_{c}=n_{b}$ and the integration
points are included in the set of solution nodes (e.g., Gauss-Lobatto points
for tensor-product elements), then pressure equilibrium is trivially
maintained [14]. However, over-integration (i.e., $n_{c}>n_{b}$) is often
necessary to minimize aliasing errors and improve stability. Unfortunately,
standard over-integration, as defined in Equation (3.8), causes a loss of
pressure equilibrium and generation of spurious pressure oscillations at
material interfaces [14]. Instead, Johnson and Kercher [14] proposed the
following approximation of the convective flux:

$\mathcal{F_{\kappa}}\approx\sum_{k=1}^{n_{c}}\mathcal{F}\left(\widetilde{y}_{\kappa}\left(x_{k}\right)\right)\varphi_{k},$
(3.9)

where $\widetilde{y}:\mathbb{R}^{m}\times\mathbb{R}\rightarrow\mathbb{R}^{m}$
is a modified state defined as

$\widetilde{y}\left(y,\widetilde{P}\right)=\left(\rho v_{1},\ldots,\rho
v_{d},\widetilde{\rho
u}\left(C_{1},\ldots,C_{n_{s}},\widetilde{P}\right)+\frac{1}{2}\sum_{k=1}^{d}\rho
v_{k}v_{k},C_{1},\ldots,C_{n_{s}}\right)^{T}.$ (3.10)

$\widetilde{P}$ is a polynomial approximation of the pressure that
interpolates onto the span of $\left\\{\phi_{1},\ldots,\phi_{n_{b}}\right\\}$
as

$\widetilde{P}_{\kappa}=\sum_{j=1}^{n_{b}}P\left(y_{\kappa}\left(x_{j}\right)\right)\phi_{j},$

and the modified internal energy, $\widetilde{\rho u}$, is evaluated from the
modified pressure and unmodified species concentrations. Interpolating the
convective flux (including the numerical flux function) in this manner
achieves pressure equilibrium both internally and between adjacent elements.
Of course, with finite resolution, slight deviations from pressure equilibrium
are inevitable; nevertheless, apart from severely underresolved computations,
these deviations remain small and do not generate large-scale pressure
oscillations that cause the solver to crash, which is not the case if standard
flux interpolation (3.8) is employed. Additional information on the basic DG
discretization, enforcement of boundary conditions, and nonlinear flux
interpolation, as well as a detailed discussion of the conditions under which
pressure oscillations are generated, can be found in [14].

##  4 Minimum entropy principle

It is well-known that in the presence of discontinuities, weak solutions to
general systems of hyperbolic conservation laws, including the multicomponent
Euler equations, are not unique [71]. As such, physical solutions are
typically identified as those that satisfy entropy conditions, written as (in
the absence of source terms)

$\frac{\partial U}{\partial t}+\nabla\cdot\mathcal{F}^{s}\leq 0,$ (4.1)

where $U(y):\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a given convex
(mathematical) entropy function and
$\mathcal{F}^{s}(y):\mathbb{R}^{m}\rightarrow\mathbb{R}^{d}$ is the
corresponding spatial entropy flux satisfying

$\mathsf{v}^{T}\frac{\partial\mathcal{F}}{\partial
y}=\frac{\partial\mathcal{F}^{s}}{\partial y},$

with $\mathsf{v}$, the entropy variables, defined as

$\mathsf{v}=\left(\frac{\partial U}{\partial y}\right)^{T}.$

The mapping from the conservative variables to the entropy variables is one-
to-one and symmetrizes the system [72]. _Entropy solutions_ are weak solutions
that satisfy (4.1) for all entropy/entropy-flux pairs. We also introduce the
entropy potential and the corresponding entropy flux potential:

$\left(\mathcal{U},\mathcal{F}^{p}\right)=\left(\mathsf{v}^{T}y-U,\mathsf{v}^{T}\mathcal{F}-\mathcal{F}^{s}\right),$
(4.2)

which will be used in Section 6. For the multicomponent Euler equations,
$U=-\rho s$ and $\mathcal{F}^{s}=-\rho sv$ form a common admissible
entropy/entropy-flux pair [16, 60], assuming $C_{i}>0$ and $T>0$. Note the
distinction between the mathematical entropy and the thermodynamic entropy
(per unit volume); they are typically of opposite sign. The conservation
equation (for smooth solutions) for $U=-\rho s$ can be obtained by first
combining the Gibbs relation,

$Tds=du-\sum_{i=1}^{n_{s}}g_{i}dY_{i}-\frac{P}{\rho^{2}}d\rho,$

with the transport equations for species mass fractions, density, and internal
energy [16, 60],

$\displaystyle\frac{DY_{i}}{Dt}$ $\displaystyle=0,\;i=1,\ldots,n_{s},$
$\displaystyle\frac{D\rho}{Dt}$ $\displaystyle+\rho\nabla\cdot v=0,$
$\displaystyle\frac{Du}{Dt}$ $\displaystyle+\frac{P}{\rho}\nabla\cdot v=0,$

to yield the specific entropy transport equation,

$\frac{Ds}{Dt}=0.$ (4.3)

Then, using conservation of mass,

$\frac{\partial\rho}{\partial t}+\nabla\cdot\left(\rho v\right)=0,$

the conservation equation for $U=-\rho s$ is obtained:

$\frac{\partial\rho s}{\partial t}+\nabla\cdot\left(\rho sv\right)=0.$

###  4.1 Review: Minimum entropy principle in the compressible,
multicomponent, nonreacting Euler equations

Gouasmi et al. [16] recently proved a minimum entropy principle satisfied by
entropy solutions to the multicomponent, nonreacting Euler equations, which
means that the spatial minimum of the specific thermodynamic entropy is a
nondecreasing function of time. In this subsection, we summarize the main
steps of the proof, which itself builds on the proof by Tadmor [73] of a
minimum entropy principle in the monocomponent Euler equations.

Integrating the inequality (4.1) over $\Omega$ and assuming nonnegative net
spatial entropy outflux across $\partial\Omega$ yields

$\frac{d}{dt}\int_{\Omega}U\left(y(x,t)\right)dx\leq 0.$

Integrating in time then gives

$\int_{\Omega}U\left(y(x,t)\right)dx\leq\int_{\Omega}U\left(y(x,0)\right)dx.$

Tadmor [74], however, showed that by integrating (4.1) over the truncated cone
$\mathsf{C}=\left\\{\left|x\right|\leq R+v_{\max}(t-\tau)|0\leq\tau\leq
t\right\\}$, a more local inequality can be written:

$\int_{\left|x\right|\leq R}U\left(y(x,t)\right)dx\leq\int_{\left|x\right|\leq
R+v_{\mathrm{max}}t}U\left(y(x,0)\right)dx,$ (4.4)

where $v_{\max}$ is the maximum speed in the domain at $t=0$. If we consider
entropy/entropy-flux pairs of the form

$\left(U,\mathcal{F}^{s}\right)=\left(-\rho f(s),-\rho vf(s)\right),$ (4.5)

where $f$ is a smooth function of $s$, (4.4) then becomes

$\int_{\left|x\right|\leq R}\rho(x,t)\cdot
f\left(s\left(y(x,t)\right)\right)dx\geq\int_{\left|x\right|\leq
R+v_{\mathrm{max}}t}\rho(x,0)\cdot f\left(s\left(y(x,0)\right)\right)dx.$
(4.6)

Consider the following choice for $f(s)$:

$f_{0}(s)=\min\left\\{s-s_{0},0\right\\},$

where $s_{0}$ is the essential infimum of the specific thermodynamic entropy
in the subdomain $\Omega_{R}=\left\\{\left|x\right|\leq R+v_{\max}t\right\\}$,

$s_{0}=\underset{\left|x\right|\leq R+v_{\max}t}{\text{Ess inf}}s(x,0).$

Although $f_{0}(s)$ is not a smooth function of $s$, it can be written as the
limit of a sequence of smooth functions,
$f_{0}(s)=\underset{\epsilon\rightarrow 0}{\lim}f_{\epsilon}(s)$, where
$f_{\epsilon}(s)$ is defined as [16]

$f_{\epsilon}(s)=\int_{-\infty}^{\infty}f_{0}(s-\mathfrak{s})g_{\epsilon}(\mathfrak{s})d\mathfrak{s},$

with

$g_{\epsilon}(\mathfrak{s})=\frac{1}{\epsilon}\frac{\exp\left(-\frac{\mathfrak{s}^{2}}{\epsilon^{2}}\right)}{\sqrt{\pi}},\;\epsilon>0.$

The first and second derivatives of $f_{\epsilon}(s)$ satisfy the following
conditions:

$\frac{df_{\epsilon}}{ds}>0,\frac{d^{2}f_{\epsilon}}{ds^{2}}<0.$

Gouasmi et al. [16] proved the key result that the entropy/entropy-flux pairs
$\left(U,\mathcal{F}^{s}\right)=\left(-\rho f_{\epsilon}(s),-\rho
vf_{\epsilon}(s)\right)$, with $\epsilon>0$, are admissible. In particular,
they showed that a conservation equation (for smooth solutions) for said pairs
can be obtained and that the entropy functions are convex with respect to the
conservative variables.

With $U=-\rho f_{0}(s)$, the inequality (4.6) can then be written as

$\int_{\left|x\right|\leq
R}\rho(x,t)\cdot\min\left\\{s(x,t)-s_{0},0\right\\}dx\geq\int_{\left|x\right|\leq
R+v_{\mathrm{max}}t}\rho(x,0)\cdot\min\left\\{s(x,0)-s_{0},0\right\\}dx,$

where the RHS is zero by definition of $s_{0}$ and the LHS is nonpositive,
yielding, for $\left|x\right|\leq R$,

$\min\left\\{s(x,t)-s_{0},0\right\\}=0.$ (4.7)

As a direct result, we obtain, for $\left|x\right|\leq R$,

$s(x,t)\geq s_{0}=\underset{\left|x\right|\leq R+v_{\max}t}{\text{Ess
inf}}s(x,0),$ (4.8)

which is the minimum entropy principle for the compressible, multicomponent
Euler equations. Note that only the entropy inequalities associated with
$U=-\rho f_{\epsilon}(s)$ need to be satisfied for a minimum entropy principle
to hold.

###  4.2 Minimum entropy principle in the compressible, multicomponent,
reacting Euler equations

We now extend the result in the previous subsection to the reacting Euler
equations, where the only difference is the inclusion of the chemical source
terms (Equation (2.7)). The presence of the chemical source terms, which are
smooth functions of only the state variables, modifies the entropy inequality
(4.1) satisfied by entropy solutions as [75, 76, 77]

$\frac{\partial U}{\partial
t}+\nabla\cdot\mathcal{F}^{s}\leq\mathsf{v}^{T}\mathcal{S},$ (4.9)

where the RHS represents the mathematical entropy production rate due to the
source term. If $\mathsf{v}^{T}\mathcal{S}\leq 0$ (i.e., the entropy
production is nonpositive), the entropy inequality in (4.1) for the
homogeneous system is recovered. The local inequality in (4.4) can then be
obtained by again integrating over $\mathsf{C}=\left\\{\left|x\right|\leq
R+v_{\max}(t-\tau)|0\leq\tau\leq t\right\\}$. If, in particular,
$\mathsf{v}^{T}\mathcal{S}\leq 0$ for entropy functions of the form $U=-\rho
f_{\epsilon}(s),\>\forall\epsilon>0$, the remaining arguments in Section 4.1
can be used to establish the same minimum entropy principle in Equation (4.8)
for the reacting Euler equations. As such, we focus on showing
$\mathsf{v}^{T}\mathcal{S}\leq 0$ with $U=-\rho
f_{\epsilon}(s),\>\forall\epsilon>0$.

Consider again entropy/entropy-flux pairs of the form
$\left(U,\mathcal{F}^{s}\right)=\left(-\rho f(s),-\rho vf(s)\right)$. The
corresponding entropy variables are given as

$\mathsf{v}=\begin{pmatrix}\frac{df}{ds}\frac{v_{1}}{T}\\\ \vdots\\\
\frac{df}{ds}\frac{v_{d}}{T}\\\ -\frac{df}{ds}\frac{1}{T}\\\
W_{1}\frac{df}{ds}\left(\frac{g_{1}-\frac{1}{2}\sum_{k=1}^{d}v_{k}v_{k}}{T}+s\right)-W_{1}f\\\
\vdots\\\
W_{n_{s}}\frac{df}{ds}\left(\frac{g_{n_{s}}-\frac{1}{2}\sum_{k=1}^{d}v_{k}v_{k}}{T}+s\right)-W_{n_{s}}f\end{pmatrix},$

which differ slightly from the entropy variables derived by Gouasmi et al.
[16] since they used partial densities instead of species concentrations in
the vector of state variables. The entropy production rate due to chemical
reactions, $\mathsf{v}^{T}\mathcal{S}$, is then written as

$\displaystyle\mathsf{v}^{T}\mathcal{S}$
$\displaystyle=\sum_{i=1}^{n_{s}}\left[W_{i}\frac{df}{ds}\left(\frac{g_{i}-\frac{1}{2}\sum_{k=1}^{d}v_{k}v_{k}}{T}+s\right)-W_{i}f\right]\omega_{i}$
$\displaystyle=\frac{df}{ds}\sum_{i=1}^{n_{s}}W_{i}\omega_{i}\frac{g_{i}}{T}+\left(-\frac{1}{2}\frac{\sum_{k=1}^{d}v_{k}v_{k}}{T}\frac{df}{ds}+s\frac{df}{ds}-f\right)\sum_{i=1}^{n_{s}}W_{i}\omega_{i}$
$\displaystyle=\frac{df}{ds}\sum_{i=1}^{n_{s}}W_{i}\omega_{i}\frac{g_{i}}{T},$

where the last equality is due to mass conservation, as given by Equation
(2.8). Since $\frac{df_{\epsilon}}{ds}>0,\>\forall\epsilon>0$, a minimum
entropy principle holds under the condition

$\sum_{i=1}^{n_{s}}W_{i}\omega_{i}\frac{g_{i}}{T}\leq 0.$ (4.10)

The term on the LHS, $\sum_{i=1}^{n_{s}}W_{i}\omega_{i}g_{i}/T$, is precisely
the entropy production rate for $U=-\rho s$, which Giovangigli [60, Chapter
6.4] already showed to be nonpositive. For completeness, we review the proof
here.

Let $\mu$ and $\mu^{\mathrm{u}}$ denote the following vectors of the reduced
chemical potentials:

$\displaystyle\mu$
$\displaystyle=\left(\mu_{1},\ldots,\mu_{n_{s}}\right)^{T},$
$\displaystyle\mu^{\mathrm{u}}$
$\displaystyle=\left(\mu_{1}^{\mathrm{u}},\ldots,\mu_{n_{s}}^{\mathrm{u}}\right)^{T}.$

We also define $\mathcal{W}$ as the matrix with the molecular weights along
the diagonal:

$\mathcal{W}=\mathrm{diag}\left(W_{1},\ldots,W_{n_{s}}\right).$

The equilibrium constant can then be rewritten as

$\log
K_{j}^{e}=-\left(\mathcal{W}\nu_{j}\right)^{T}\mu^{\mathrm{u}}=\left(\mathcal{W}\nu_{j}^{f}\right)^{T}\mu^{\mathrm{u}}-\left(\mathcal{W}\nu_{j}^{r}\right)^{T}\mu^{\mathrm{u}},$
(4.11)

where

$\displaystyle\nu_{j}^{f}$
$\displaystyle=\left(\nu_{1j}^{f},\ldots,\nu_{n_{s}j}^{f}\right)^{T},$
$\displaystyle\nu_{j}^{r}$
$\displaystyle=\left(\nu_{1j}^{r},\ldots,\nu_{n_{s}j}^{r}\right)^{T},$
$\displaystyle\nu_{j}$ $\displaystyle=\nu_{j}^{r}-\nu_{j}^{f}.$

Using Equation (4.11), Giovangigli [60, Chapter 6.4] introduced a new reaction
constant, given by

$\log K_{j}^{s}=\log
K_{j}^{f}-\left(\mathcal{W}\nu_{j}^{f}\right)^{T}\mu^{\mathrm{u}}=\log
K_{j}^{r}-\left(\mathcal{W}\nu_{j}^{r}\right)^{T}\mu^{\mathrm{u}}.$

The rate of progress of the $j$th reaction can then be expressed as

$\displaystyle q_{j}$
$\displaystyle=K_{j}^{s}\left[\exp\left(\left(\mathcal{W}\nu_{j}^{f}\right)^{T}\mu^{\mathrm{u}}\right)\prod_{i=1}^{n_{s}}C_{i}^{\nu_{ij}^{f}}-\exp\left(\left(\mathcal{W}\nu_{j}^{r}\right)^{T}\mu^{\mathrm{u}}\right)\prod_{i=1}^{n_{s}}C_{i}^{\nu_{ij}^{r}}\right],$
$\displaystyle=K_{j}^{s}\left[\exp\left(\left(\mathcal{W}\nu_{j}^{f}\right)^{T}\mu\right)-\exp\left(\left(\mathcal{W}\nu_{j}^{r}\right)^{T}\mu\right)\right].$

We rewrite the entropy production rate for $U=-\rho s$ as

$\displaystyle\sum_{i=1}^{n_{s}}W_{i}\omega_{i}\frac{g_{i}}{T}$
$\displaystyle=R^{0}\sum_{i=1}^{n_{s}}\sum_{j=1}^{n_{r}}W_{i}\mu_{i}\nu_{ij}q_{j}$
$\displaystyle=-R^{0}\sum_{j=1}^{n_{r}}\mu^{T}\mathcal{W}\left(\nu_{j}^{f}-\nu_{j}^{r}\right)q_{j}$
$\displaystyle=-R^{0}\sum_{j=1}^{n_{r}}K_{j}^{s}\left(\mu^{T}\mathcal{W}\nu_{j}^{f}-\mu^{T}\mathcal{W}\nu_{j}^{r}\right)\left[\exp\left(\mu^{T}\mathcal{W}\nu_{j}^{f}\right)-\exp\left(\mu^{T}\mathcal{W}\nu_{j}^{r}\right)\right],$
$\displaystyle\leq 0.$

Note that this is equivalent to a nonnegative production rate of
_thermodynamic_ entropy per unit volume. Since the condition in (4.10) is
satisfied, a minimum entropy principle in the compressible, multicomponent,
reacting Euler equations holds. Another consequence of (4.10) is that the
entropy production rate for any entropy function of the form $U=-\rho f(s)$ is
nonpositive provided that $\frac{df}{ds}\geq 0$.

##  5 Transport step: Entropy-bounded discontinuous Galerkin scheme in one
dimension

In this section, we detail the entropy-bounded DG methodology for solving
Equations (3.3) and (3.5) (i.e., the explicit time integrations without source
terms in the operator splitting strategy) while accounting for the modified
flux interpolation in Equation (3.9). We build on related entropy-bounded DG
schemes for the monocomponent Euler equations [36, 17, 22]. These schemes are
formulated as extensions of positivity-preserving DG methods [20, 34, 35]
since the thermodynamic entropy is well-defined only for positive densities
and pressures and enforcement of an entropy constraint can be
straightforwardly incorporated into the positivity-preserving framework. The
one-dimensional entropy-bounded DG scheme is presented in Section (5.2), then
extended to multiple dimensions in Part II [19].

###  5.1 Preliminaries

Let $\mathcal{G}_{\sigma}$ denote the following set:

$\mathcal{G_{\sigma}}=\left\\{y\mid C_{1}>0,\ldots,C_{n_{s}}>0,\rho
u^{*}>0,s\geq\sigma\right\\},$ (5.1)

where $\sigma\in\mathbb{R}$ and $u^{*}$ is the “shifted” internal energy [48],
computed as

$u^{*}=u-u_{0}=u-\sum_{i=1}^{n_{s}}Y_{i}b_{i0},$ (5.2)

such that $u^{*}>0$ if and only if $T>0$, provided
$c_{v,i}>0,\>i=1,\ldots,n_{s}$ [60]. Pressure is then also positive. Note that
$C_{i}>0,\;\forall i$, implies $\rho>0$. This set is similar to that in [45,
46], except with the additional entropy constraint. In A, we review two
important lemmas from [36], the first of which establishes the quasi-concavity
of $s(y)$ (which holds in the multicomponent case as well) and the second of
which states
$s\left(\overline{y}_{\kappa}\right)\geq\min_{x\in\kappa}s\left(y(x)\right)$.
In B, we show that $\rho u^{*}(y)$ is a concave function of the state. For a
given $\sigma$, $\mathcal{G}_{\sigma}$ is then a convex set [34, 35, 48, 78].
We assume that the exact solution to the classical Riemann problem with
initial data

$y\left(x,0\right)=\begin{cases}y_{1},&x<0\\\ y_{2},&x>0\end{cases}$

is an entropy solution that preserves positivity. Then, by Lemmas 10 and 11,
$\mathcal{G}_{\sigma}$ is an invariant set [39, 79], i.e.,
$y_{1}\in\mathcal{G}_{\sigma}$ and $y_{2}\in\mathcal{G}_{\sigma}$ imply that
the average of the exact Riemann solution over a domain that includes the
Riemann fan [79] is in $\mathcal{G}_{\sigma}$.

Consider the following three-point system arising from a $p=0$, element-local
DG discretization with forward Euler time stepping:

$\displaystyle y_{\kappa}^{j+1}=y_{\kappa}^{j}-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(y_{\kappa}^{j},y_{\kappa_{L}}^{j},-1\right)+\mathcal{F}^{\dagger}\left(y_{\kappa}^{j},y_{\kappa_{R}}^{j},1\right)\right],$
(5.3)

where $j$ indexes the time step, $h$ is the element size, and $\kappa_{L}$ and
$\kappa_{R}$ are the elements to the left and right of $\kappa$, respectively.
Let $\lambda$ be an upper bound on the maximum wave speed of the system. Under
the following condition,

$\frac{\Delta t\lambda}{h}\leq\frac{1}{2},$ (5.4)

$y_{\kappa}^{j},y_{\kappa_{L}}^{j},y_{\kappa_{R}}^{j}\in\mathcal{G}_{\sigma}$
implies $y_{\kappa}^{j+1}\in\mathcal{G}_{\sigma}$ if certain _invariant-
region-preserving_ numerical fluxes are employed [22]. In particular,
$y_{\kappa}^{j+1}$ satisfies [22, 78]

$s\left(y_{\kappa}^{j+1}\right)\geq\min\left\\{s\left(y_{\kappa_{L}}^{j}\right),s\left(y_{\kappa}^{j}\right),s\left(y_{\kappa_{R}}^{j}\right)\right\\}.$
(5.5)

The Godunov, Lax-Friedrichs, HLL, and HLLC fluxes qualify [22] (see also [73,
17, 16, 80] for additional proofs regarding the Godunov and/or Lax-Friedrichs
fluxes), some of which allow for less restrictive time-step-size constraints
than (5.4). The proofs often rely on the notion that $\mathcal{G}_{\sigma}$ is
an invariant set, which itself invokes the aforementioned assumption that the
exact Riemann solution is an entropy solution satisfying the positivity
property. Two exceptions are the proof by Zhang and Shu [20] of the positivity
property of the Lax-Friedrichs flux and the proof by Lax [81] that the Lax-
Friedrichs flux satisfies a discrete cell entropy inequality for all
entropy/entropy-flux pairs, from which a sharper version of the inequality
(5.5) follows [73, 16] (note that the corresponding time-step-size constraints
for these two proofs are not necessarily the same as (5.4)). Unless otherwise
specified, we employ the HLLC numerical flux [82].

The three-point system (5.3) will be crucial in the construction of a
positivity-preserving and entropy-bounded DG scheme for $p>0$. Specifically,
we will show in the following subsection that the element average of the
solution (for $p>0$) at the next time step, $\overline{y}_{\kappa}^{j+1}$, can
be expressed as a convex combination of both pointwise values of
$y_{\kappa}^{j}(x)$ and three-point systems involving pointwise values of
$y_{\kappa}^{j}(x)$. If all of said pointwise values of $y_{\kappa}^{j}(x)$
are in $\mathcal{G}_{\sigma}$, then $\overline{y}_{\kappa}^{j+1}$ will also be
in $\mathcal{G}_{\sigma}$ under a time-step-size constraint. A simple limiter,
described in Section 5.2.1, is applied to ensure that the pointwise values of
$y_{\kappa}^{j}(x)$ are in $\mathcal{G}_{\sigma}$.

###  5.2 Entropy-bounded, high-order discontinuous Galerkin method in one
dimension

Suppose $\kappa=\left[x_{L},x_{R}\right]$. Let $x_{q}$ and $w_{q}$ denote the
quadrature points and weights, respectively, of a quadrature rule with
$x_{q}\in\kappa$ , $w_{q}>0$, and $\sum_{q=1}^{n_{q}}w_{q}=1$, where $n_{q}$
is the number of quadrature points/weights. This set of quadrature points does
not need to include the endpoints, and the quadrature rule need not be
explicitly used to evaluate any integrals in Equation (3.2). For now, we
assume that standard flux interpolation, as in Equation (3.8), is employed;
Section 5.2.2 discusses how to account for the modified flux interpolation in
Equation (3.9). As in [17], provided that the quadrature rule is sufficiently
accurate, the element-averaged solution in Equation (3.7) can be expanded as

$\displaystyle\overline{y}_{\kappa}$
$\displaystyle=\sum_{q=1}^{n_{q}}w_{q}y_{\kappa}(x_{q})$
$\displaystyle=\sum_{q=1}^{n_{q}}\theta_{q}y_{\kappa}(x_{q})+\theta_{L}y_{\kappa}(x_{L})+\theta_{R}y_{\kappa}(x_{R}).$
(5.6)

If the set of quadrature points includes the endpoints, then we can simply
take

$\theta_{q}=\begin{cases}w_{q}&x_{q}\neq x_{L},x_{q}\neq x_{R}\\\
0&\mathrm{otherwise}\end{cases}$

and

$\theta_{L}=w_{L},\quad\theta_{R}=w_{R},$

where $w_{L}$ and $w_{R}$ are the quadrature weights at the left and right
endpoints, respectively. If the set of quadrature points does not include the
endpoints, then we can instead take

$\theta_{q}=w_{q}-\theta_{L}\psi_{q}\left(x_{L}\right)-\theta_{R}\psi_{q}\left(x_{R}\right),$

where $\left\\{\psi_{1},\ldots,\psi_{n_{d}}\right\\}$, with $n_{b}\leq
n_{d}\leq n_{q}$, is a set of Lagrange basis functions whose nodes are located
at a subset of the quadrature points, while $\psi_{v}=0$ for
$v=n_{d}+1,\ldots,n_{q}$, such that [17]

$\displaystyle\sum_{q=1}^{n_{q}}\theta_{q}y_{\kappa}(x_{q})$
$\displaystyle=\sum_{q=1}^{n_{q}}\left[w_{q}-\theta_{L}\psi_{q}\left(x_{L}\right)-\theta_{R}\psi_{q}\left(x_{R}\right)\right]y_{\kappa}(x_{q})$
$\displaystyle=\sum_{q=1}^{n_{q}}w_{q}y_{\kappa}(x_{q})-\theta_{L}\sum_{q=1}^{n_{q}}y_{\kappa}(x_{q})\psi_{q}\left(x_{L}\right)-\theta_{R}\sum_{q=1}^{n_{q}}y_{\kappa}(x_{q})\psi_{q}\left(x_{R}\right)$
$\displaystyle=\sum_{q=1}^{n_{q}}w_{q}y_{\kappa}(x_{q})-\theta_{L}y_{\kappa}(x_{L})+\theta_{R}y_{\kappa}(x_{R}).$

$\theta_{L}$ and $\theta_{R}$ will be related to a constraint on the time step
size later in this section. The positivity of the quadrature weights
guarantees the existence of positive $\theta_{L}$ and $\theta_{R}$ that yield
$\theta_{q}\geq 0,\;q=1,\ldots,n_{q}$ [17]. Furthermore,
$\sum_{q}\theta_{q}+\theta_{L}+\theta_{R}=1$ since $\sum_{q}\psi_{q}=1$.
Define $\partial\mathcal{D}_{\kappa}=\left\\{x_{L,}x_{R}\right\\}$, and let
$\mathcal{D}_{\kappa}$ denote the set of points at which the state is
evaluated in Equation (5.6):

$\mathcal{D_{\kappa}}=\partial\mathcal{D}_{\kappa}\bigcup\left\\{x_{q},q=1,\ldots,n_{q}\right\\}=\left\\{x_{L,}x_{R},x_{q},q=1,\ldots,n_{q}\right\\}.$

Applying forward Euler time stepping to Equation (3.2) and taking
$\mathfrak{v}$ to be a vector of ones (i.e., $\mathfrak{v}\in V_{h}^{0}$)
gives the fully discrete scheme satisfied by the element averages [36, 17]:

$\displaystyle\overline{y}_{\kappa}^{j+1}$ $\displaystyle=$
$\displaystyle\overline{y}_{\kappa}^{j}-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{L}),y_{\kappa_{L}}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{R}),y_{\kappa_{R}}^{j}(x_{R}),1\right)\right]$
(5.8) $\displaystyle=$
$\displaystyle\sum_{q=1}^{n_{q}}\theta_{q}y_{\kappa}^{j}(x_{q})+\theta_{L}y_{\kappa}^{j}(x_{L})-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{L}),y_{\kappa_{L}}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{L}),y_{\kappa}^{j}(x_{R}),1\right)\right]$
$\displaystyle+\theta_{R}y_{\kappa}^{j}(x_{R})-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{R}),y_{\kappa}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{R}),y_{\kappa_{R}}^{j}(x_{R}),1\right)\right],$

where the second equality is due to the conservation property of the numerical
flux:

$\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{L}),y_{\kappa}^{j}(x_{R}),1\right)=-\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{R}),y_{\kappa}^{j}(x_{L}),-1\right).$

Note that Equations (5.8) and (5.8) hold regardless of whether the integrals
in Equation (3.2) are evaluated using conventional quadrature or a quadrature-
free implementation [69, 70]. Though the forward Euler time integration scheme
is used here, strong-stability-preserving Runge-Kutta (SSPRK) methods [83,
84], which are convex combinations of forward Euler steps, are compatible as
well. Equation (5.8) then leads to the following theorem, where we use
$y_{\kappa}^{-}$ to denote the exterior state along $\partial\kappa$.

######  Theorem 1 ([20, 36, 17]).

If $y_{\kappa}^{j}(x)\in\mathcal{G}_{\sigma},\;\forall
x\in\mathcal{D_{\kappa}}$, and
$y_{\kappa}^{-,j}\in\mathcal{G}_{\sigma},\;\forall
x\in\partial\mathcal{D}_{\kappa}$, with

$\sigma\leq\min\left\\{\min\left\\{s\left(y_{\kappa}^{j}(x)\right)|x\in\mathcal{D_{\kappa}}\right\\},\min\left\\{s\left(y_{\kappa}^{-,j}(x)\right)|x\in\mathcal{\partial
D_{\kappa}}\right\\}\right\\},$ (5.9)

then $\overline{y}_{\kappa}^{j+1}$ in Equation (5.8) is also in
$\mathcal{G}_{\sigma}$ under the constraint

$\frac{\Delta
t\lambda}{h}\leq\frac{1}{2}\min\left\\{\theta_{L},\theta_{R}\right\\}$ (5.10)

and the conditions

$\theta_{L}>0,\theta_{R}>0,\theta_{q}\geq 0,q=1,\ldots,n_{q}.$ (5.11)

###### Proof.

The proof follows the same procedure as in [20], [36], [17], and related
papers, which we review here. We first rewrite Equation (5.8) as

$\overline{y}_{\kappa}^{j+1}=\sum_{q=1}^{n_{q}}\theta_{q}y_{\kappa}^{j}(x_{q})+\theta_{L}y_{\kappa,s1}^{j+1}+\theta_{R}y_{\kappa,s2}^{j+1},$

where

$\displaystyle y_{\kappa,s1}^{j+1}$
$\displaystyle=y_{\kappa}^{j}(x_{L})-\frac{\Delta
t}{\theta_{L}h}\left[\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{L}),y_{\kappa_{L}}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{L}),y_{\kappa}^{j}(x_{R}),1\right)\right],$
$\displaystyle y_{\kappa,s2}^{j+1}$
$\displaystyle=y_{\kappa}^{j}(x_{R})-\frac{\Delta
t}{\theta_{R}h}\left[\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{R}),y_{\kappa}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(y_{\kappa}^{j}(x_{R}),y_{\kappa_{R}}^{j}(x_{R}),1\right)\right].$

As such, $\overline{y}_{\kappa}^{j+1}$ is a convex combination of
$y_{\kappa}^{j}(x_{q})$ evaluated at $x_{q}$ and two three-point systems of
the type (5.3). Under the conditions (5.10) and (5.11), $y_{\kappa,s1}^{j+1}$
and $y_{\kappa,s2}^{j+1}$ are both in $\mathcal{G}_{\sigma}$. It then follows
that $\overline{y}_{\kappa}^{j+1}\in\mathcal{G}_{\sigma}$.

∎

######  Remark 2.

A direct result of Theorem 1 is that

$s\left(\overline{y}_{\kappa}^{j+1}\right)\geq\min\left\\{\min\left\\{s\left(y_{\kappa}^{j}(x)\right)|x\in\mathcal{D_{\kappa}}\right\\},\min\left\\{s\left(y_{\kappa}^{-,j}(x)\right)|x\in\mathcal{\partial
D_{\kappa}}\right\\}\right\\}.$

According to the inequality (5.10), the upper bound on the time step size is
proportional to $\min\left\\{\theta_{L},\theta_{R}\right\\}$. For Gauss-
Lobatto quadrature, $\theta_{L}$ and $\theta_{R}$ are both equal to the
quadrature weight corresponding to either endpoint. See [17] for information
about Gauss-Legendre quadrature, as well as a discussion on how to find the
maximum allowable value of $\min\left\\{\theta_{L},\theta_{R}\right\\}$ for
general quadrature rules.

To complete the construction of an entropy-bounded, high-order DG scheme, we
need to enforce not only the positivity of $y_{\kappa}^{j+1}(x),\>\forall
x\in\mathcal{D_{\kappa}}$, for all $\kappa\in\Omega$, but also
$s\left(y^{j+1}(x)\right)>s_{b},\>\forall x\in\mathcal{D_{\kappa}}$, for all
$\kappa\in\Omega$, where $s_{b}$ is a lower bound on the specific
thermodynamic entropy. $s_{b}$ can vary among elements and time steps, and the
prescription of $s_{b}$ should be motivated by the physical principles
examined in Section 4. We will discuss $s_{b}$ in more detail in Section
5.2.3. In other words, we enforce
$y_{\kappa}^{j+1}(x)\in\mathcal{G}_{s_{b}},\>\forall
x\in\mathcal{D_{\kappa}}$, which is done via a simple limiting procedure that
will be described in Section 5.2.1.

In practice, we relax some of the requirements introduced thus far. First,
$\lambda$ is computed in a local (instead of global) manner and is calculated
as the maximum value of $\left|v\right|+c$, where $c$ is the speed of sound,
over the points of interest. However, $\left|v\right|+c$ does not bound the
wave speeds arising from the interactions between states at interfaces. A
similar remark can be made for most invariant-region-preserving numerical flux
functions, which typically require wave-speed estimates. Simple algorithms for
bounding the wave speeds in the monocomponent case have been developed [85,
86]; extending these to the multicomponent Euler equations may indeed be
worthy of future investigation, with [87] as one example. Second, we introduce
$\chi_{\sigma}=\rho s-\rho\sigma$, which is concave with respect to the state
[22], and revise the definition of $\mathcal{G}_{\sigma}$ as

$\mathcal{G}_{\sigma}=\left\\{y\mid\rho>0,\rho u^{*}>0,C_{1}\geq
0,\ldots,C_{n_{s}}\geq 0,\chi_{\sigma}\geq 0\right\\},$ (5.12)

where $\chi_{\sigma}\geq 0$ is a reformulation of $s\geq\sigma$ and the
species concentrations are now allowed to be equal to zero. From a practical
standpoint, allowing $C_{i}=0$ is necessary since the concentrations are
frequently zero in many reacting flow problems of interest. Unfortunately,
entropy functions of the form $U=-\rho f_{\epsilon}(s)$ and $U=-\rho s$ are no
longer convex if any of the concentrations is zero [16, 51]. Furthermore, the
specific thermodynamic entropy becomes ill-defined. Nevertheless, by making
use of $0\log 0=0$ [60, Chapter 6], $\rho s$ and thus $\chi_{\sigma}$ remain
well-defined. The entire methodology developed here also remains well-defined,
unlike entropy-stable schemes that rely on the entropy variables associated
with $U=-\rho s$. Throughout this work, we did not encounter any major issues
associated with relaxing the two aforementioned requirements. One potential
reason is that $\overline{y}_{\kappa}^{j+1}$ can be in $\mathcal{G}_{\sigma}$
even if $\mathcal{G}_{\sigma}$ is not convex and/or some of the conditions in
Theorem 1 are not satisfied. Furthermore, in this work, we choose
$\mathrm{CFL}=0.1$ to maintain low temporal errors, which in general yields a
smaller time step size than necessary. Should issues emerge in future work,
they can likely be alleviated by adaptively decreasing the time step size
[35].

Finally, we remark that $\mathcal{D}_{\kappa}$ is simply the set of points at
which limiting should be applied. Specifically, the interior quadrature points
in $\mathcal{D}_{\kappa}$ need not be explicitly used in numerical
integrations in Equation (3.2); if they are indeed not, the actual integration
points are added to $\mathcal{D}_{\kappa}$ as well [35].

####  5.2.1 Limiting procedure

In this subsection, we describe the limiting procedure to ensure
$y_{\kappa}^{j+1}(x)\in\mathcal{G}_{s_{b}},\>\forall
x\in\mathcal{D_{\kappa}}$. It is assumed that
$\overline{y}_{\kappa}^{j+1}(x)\in\mathcal{G}_{s_{b}}$. For brevity, we drop
the $j+1$ superscript and $\kappa$ subscript in this discussion. The limiting
operator is of the same form as in [44], [35], [22], [78], and related papers.

  1. 1.

First, positivity of density is enforced. Specifically, if
$\rho(x)>\epsilon,\>\forall x\in\mathcal{D}_{\kappa}$, where $\epsilon$ is a
small positive number (e.g., $\epsilon=10^{-10}$), then set
$C_{i}^{(1)}=C_{i}=\sum_{j=1}^{n_{b}}C_{i}(x_{j})\phi_{j},i=1,\ldots,n_{s}$;
otherwise, compute

$C_{i}^{(1)}=\overline{C}_{i}+\theta^{(1)}\left(C_{i}-\overline{C}_{i}\right),\;i=1,\ldots,n_{s},$

with

$\theta^{(1)}=\frac{\rho(\overline{y})-\epsilon}{\rho(\overline{y})-\underset{x\in\mathcal{D}}{\min}\rho(y(x))}.$

  2. 2.

Next, nonnegativity of the species concentrations is enforced. If
$C_{i}^{(1)}(x)\geq 0,\>\forall x\in\mathcal{D}_{\kappa}$, then set
$C_{i}^{(2)}=C_{i}^{(1)},i=1,\ldots,n_{s}$; otherwise, compute

$C_{i}^{(2)}=\overline{C}_{i}+\theta^{(2)}\left(C_{i}^{(1)}-\overline{C}_{i}\right),\;i=1,\ldots,n_{s},$

with

$\theta^{(2)}=\frac{\overline{C}_{i}}{\overline{C}_{i}-\underset{x\in\mathcal{D}}{\min}C_{i}^{(1)}(x)}.$

Let $y^{(2)}=\left(\rho v_{1},\ldots,\rho v_{d},\rho
e_{t},C_{1}^{(2)},\ldots,C_{n_{s}}^{(2)}\right)$.

  3. 3.

Positivity of $\rho u^{*}(y)$ is then enforced. If $\rho
u^{*}\left(y^{(2)}(x)\right)>\epsilon,\>\forall x\in\mathcal{D}_{\kappa}$,
then set $y^{(3)}=y^{(2)}$; otherwise, compute

$y^{(3)}=\overline{y}+\theta^{(3)}\left(y^{(2)}-\overline{y}\right),$

with

$\theta^{(3)}=\frac{\rho u^{*}(\overline{y})-\epsilon}{\rho
u^{*}(\overline{y})-\underset{x\in\mathcal{D}}{\min}\rho u^{*}(y^{(2)}(x))}.$

It can be shown that $\rho u^{*}(y^{(3)}(x))>0,\>\forall
x\in\mathcal{D}_{\kappa}$ by concavity [44, 35]. The “positivity-preserving
limiter” refers to the limiting procedure up to this point. The “entropy
limiter” corresponds to the following step (in addition to the above steps).

  4. 4.

Finally, the entropy constraint is enforced. If
$\chi\left(y^{(3)}(x)\right)\geq 0,\>\forall x\in\mathcal{D}_{\kappa}$, then
set $y^{(4)}=y^{(3)}$; otherwise, compute

$y^{(4)}=\overline{y}+\theta^{(4)}\left(y^{(3)}-\overline{y}\right),$

with

$\theta^{(4)}=\frac{\chi(\overline{y})}{\chi(\overline{y})-\underset{x\in\mathcal{D}}{\min}\chi(y^{(3)}(x))}.$

It can be shown that $s\left(y^{(4)}(x)\right)\geq s_{b},\>\forall
x\in\mathcal{D}_{\kappa}$, by concavity of $\chi$ [22, 78].

$y^{(4)}$ then replaces $y$ as the solution. The limiting operator is
conservative, maintains stability, and in general preserves the formal order
of accuracy for smooth solutions [20, 35, 36, 17, 22]. There is extensive
empirical evidence demonstrating preservation of accuracy (see previously
cited references, as well as [46], [78], and related papers). However, the
order of accuracy can potentially deteriorate when the element average is
close to the boundary of $\mathcal{G}_{s_{b}}$ [35, 22]. Furthermore, the
linear scaling is not expected to suppress all oscillations [20, 22, 17, 78].
The limiting procedure described here is applied at the end of every RK stage.

####  5.2.2 Modified flux interpolation

We now discuss how to account for over-integration with the modified flux
interpolation in Equation (3.9). The scheme satisfied by the element averages
becomes

$\displaystyle\overline{y}_{\kappa}^{j+1}$ $\displaystyle=$
$\displaystyle\overline{y}_{\kappa}^{j}-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{L}),\widetilde{y}_{\kappa_{L}}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{R}),\widetilde{y}_{\kappa_{R}}^{j}(x_{R}),1\right)\right]$
(5.13) $\displaystyle=$
$\displaystyle\sum_{q=1}^{n_{q}}\theta_{q}y_{\kappa}^{j}(x_{q})+\theta_{L}y_{\kappa}^{j}(x_{L})-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{L}),\widetilde{y}_{\kappa_{L}}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{L}),\widetilde{y}_{\kappa}^{j}(x_{R}),1\right)\right]$
$\displaystyle+\theta_{R}y_{\kappa}^{j}(x_{R})-\frac{\Delta
t}{h}\left[\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{R}),\widetilde{y}_{\kappa}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{R}),\widetilde{y}_{\kappa_{R}}^{j}(x_{R}),1\right)\right].$

Equation (5.13) can be rewritten as

$\overline{y}_{\kappa}^{j+1}=\sum_{q=1}^{n_{q}}\theta_{q}y_{\kappa}^{j}(x_{q})+\theta_{L}y_{\kappa,s3}^{j+1}+\theta_{R}y_{\kappa,s4}^{j+1},$

where

$\displaystyle y_{\kappa,s3}^{j+1}$
$\displaystyle=y_{\kappa}^{j}(x_{L})-\frac{\Delta
t}{\theta_{L}h}\left[\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{L}),\widetilde{y}_{\kappa_{L}}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{L}),\widetilde{y}_{\kappa}^{j}(x_{R}),1\right)\right],$
$\displaystyle y_{\kappa,s4}^{j+1}$
$\displaystyle=y_{\kappa}^{j}(x_{R})-\frac{\Delta
t}{\theta_{R}h}\left[\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{R}),\widetilde{y}_{\kappa}^{j}(x_{L}),-1\right)+\mathcal{F}^{\dagger}\left(\widetilde{y}_{\kappa}^{j}(x_{R}),\widetilde{y}_{\kappa_{R}}^{j}(x_{R}),1\right)\right],$

which are not necessarily of the type (5.3) since in general,
$y_{\kappa}^{j}(x_{L})\neq\widetilde{y}_{\kappa}^{j}(x_{L})$ and
$y_{\kappa}^{j}(x_{R})\neq\widetilde{y}_{\kappa}^{j}(x_{R})$. The
incompatibility is a result of expressing $\overline{y}_{\kappa}$ as a convex
combination of pointwise values of $y_{\kappa}(x)$ (as opposed to
$\widetilde{y}_{\kappa}(x)$). Unfortunately, the element average of
$\widetilde{y}_{\kappa}$, denoted $\overline{\widetilde{y}}_{\kappa}$, is not
necessarily equal to $\overline{y}_{\kappa}$; consequently,
$\overline{y}_{\kappa}$ cannot be directly written as a convex combination of
pointwise values of $\widetilde{y}_{\kappa}(x)$. However, if the set of
solution nodes includes the endpoints (e.g., equidistant or Gauss-Lobatto
points), then Equation (5.13) recovers Equation (5.8) since
$y_{\kappa}^{j}(x_{L})=\widetilde{y}_{\kappa}^{j}(x_{L})$ and
$y_{\kappa}^{j}(x_{R})=\widetilde{y}_{\kappa}^{j}(x_{R})$. The previous
analysis, including Theorem 1, then holds. As such, for a nodal set that
includes the endpoints, the modified flux interpolation in Equation (3.9) does
not introduce any additional difficulties to the discussed framework. Note
also that the second term in Equation (3.2) (i.e., the volumetric flux
integral) does not factor into the scheme satisfied by the element averages,
so the modified flux interpolation can be freely employed in said integral.

####  5.2.3 Lower bound on specific thermodynamic entropy

We consider two options for specifying $s_{b}$. The first option is a _global_
entropy bound [36, 22, 78]:

$s_{b}(y)=\min\left\\{s\left(y(x)\right)|x\in\Omega\right\\},$ (5.14)

which can be evaluated once based on the initial condition, $y_{0}(x)$, or
updated at each time step [16]. Instead of calculating the true minimum via,
for example, Newton’s method, we either rely on user-specified information or
compute

$s_{b}(y)=\min\left\\{s\left(y(x)\right)|x\in\bigcup_{\kappa\in\mathcal{T}}\mathcal{D_{\kappa}}\right\\}.$
(5.15)

The second option is a _local_ entropy bound, which should satisfy

$s_{b,\kappa}^{j+1}(y)\leq\min\left\\{\min\left\\{s\left(y_{\kappa}^{j}(x)\right)|x\in\mathcal{D_{\kappa}}\right\\},\min\left\\{s\left(y_{\kappa}^{-,j}(x)\right)|x\in\mathcal{\partial
D_{\kappa}}\right\\}\right\\},$

in order to ensure compatibility with Theorem 1 and the limiting procedure,
which enforces
$y_{\kappa}^{j+1}(x)\in\mathcal{G}_{s_{b,\kappa}^{j+1}},\>\forall
x\in\mathcal{D_{\kappa}}$. Lv and Ihme [17] introduced the following local
entropy bound:

$s_{b,\kappa}^{j+1}(y)=\min\left\\{\min\left\\{s\left(y_{\kappa}^{j}(x)\right)|x\in\kappa\right\\},\min\left\\{s\left(y_{\kappa}^{-,j}(x)\right)|x\in\partial\kappa\right\\}\right\\}.$
(5.16)

They demonstrated that the local entropy bound (5.16) can more effectively
dampen overshoots and undershoots when the entropy varies significantly
throughout the domain (e.g., when multiple discontinuities are present).
However, we find that Equation (5.16) is often too restrictive, as will be
illustrated in Section 7.1.

Here, we employ a different local entropy bound,

$s_{b,\kappa}^{j+1}(y)=\min\left\\{s\left(y^{j}(x)\right)|x\in\kappa\cup\kappa_{L}\cup\kappa_{R}\right\\},$
(5.17)

which is based on the local minimum entropy principle satisfied by exact
entropy solutions (or discrete entropy solutions in the limit of infinite
resolution). Specifically, the inequality

$s(y(x,t))\geq\min\left\\{s\left(y\left(x,t_{0}\right)\right)|x\in\kappa\cup\kappa_{L}\cup\kappa_{R}\right\\},\quad
t\in\left[t_{0},t_{0}+\Delta t\right],$ (5.18)

can be considered the semi-discrete analog of (4.8) with $R=h/2$ and
$\kappa=\left[-h/2,h/2\right]$, under the condition

$\frac{\Delta tv_{\max}}{h}\leq 1.$ (5.19)

The RHS of Equation (5.17) is the fully discrete analog of the RHS of the
inequality (5.18). This property is imposed on the discrete solution as a
means to suppress instabilities. A similar bound was employed by Dzanic and
Witherden [88] in their entropy-based filtering framework. Note that if
$\Delta t$ satisfies (5.10), then it also satisfies (5.19) since

$\Delta
t\leq\frac{h}{2\lambda}\min\left\\{\theta_{L},\theta_{R}\right\\}\leq\frac{h}{2\lambda}\leq\frac{h}{\lambda}\leq\frac{h}{v_{\max}}.$

An alternate viewpoint draws from the generalized Riemann problem (GRP) [82],
which differs from the _classical_ Riemann problem by allowing for source
terms and piecewise smooth (as opposed to piecewise constant) initial
conditions. Suppose that $y_{\kappa_{L}}$ and $y_{\kappa}$ form the initial
conditions of a GRP centered at $x_{L}$, while $y_{\kappa}$ and
$y_{\kappa_{R}}$ form the initial conditions of a GRP centered at $x_{R}$.
Under the (potentially strong) assumption that exact solutions
($y_{L}^{\mathrm{GRP}}$ and $y_{R}^{\mathrm{GRP}}$, respectively) to the GRPs
exist, as well as the assumption that those solutions satisfy all entropy
inequalities, the exact solution in $\kappa$, $y_{\kappa}^{\mathrm{ex}}$,
arising from the GRPs satisfies [71]

$\int_{\kappa}U\left(y_{\kappa}^{\mathrm{ex}}(x,t)\right)\leq\int_{\kappa}U\left(y_{\kappa}(x,t_{0})\right)dx-\int_{t_{0}}^{t}\mathcal{F}^{s}\left(y_{R}^{\mathrm{GRP}}\left(x_{R},\tau\right)\right)d\tau+\int_{t_{0}}^{t}\mathcal{F}^{s}\left(y_{L}^{\mathrm{GRP}}\left(x_{L},\tau\right)\right)d\tau,$

at least before the local Riemann problems interact [89]. With
$\left(U,\mathcal{F}^{s}\right)=\left(-\rho f_{0}(s),-\rho
v_{1}f_{0}(s)\right)$ and
$s_{0}=\min\left\\{s\left(y\left(x,t_{0}\right)\right)|x\in\kappa\cup\kappa_{L}\cup\kappa_{R}\right\\}$,
the above inequality becomes

$\displaystyle\int_{\kappa}\rho\left(y_{\kappa}^{\mathrm{ex}}(x,t)\right)f_{0}\left(s\left(y_{\kappa}^{\mathrm{ex}}(x,t)\right)\right)dx$
$\displaystyle\geq$
$\displaystyle\int_{\kappa}\rho\left(y_{\kappa}(x,t_{0})\right)\min\left\\{s\left(y_{\kappa}(x,t_{0})\right)-s_{0},0\right\\}dx$
$\displaystyle-\int_{t_{0}}^{t}\rho
v_{1}\left(y_{R}^{\mathrm{GRP}}\left(x_{R},\tau\right)\right)\min\left\\{s\left(y_{R}^{\mathrm{GRP}}\left(x_{R},\tau\right)\right)-s_{0},0\right\\}dx$
$\displaystyle+\int_{t_{0}}^{t}\rho
v_{1}\left(y_{L}^{\mathrm{GRP}}\left(x_{L},\tau\right)\right)\min\left\\{s\left(y_{L}^{\mathrm{GRP}}\left(x_{L},\tau\right)\right)-s_{0},0\right\\}dx,$

where the first term on the RHS clearly vanishes and the second and third
terms vanish due to Equation (4.8). As such, we have
$s\left(y_{\kappa}^{\mathrm{ex}}(x,t)\right)\geq s_{0}$. This property is then
imposed on the discrete solution via the local entropy bound in Equation
(5.17).

To our knowledge, it is unclear whether the assumption that exact solutions to
the GRP exist is valid (even for calorically perfect gases, let alone mixtures
of thermally perfect gases). Semi-analytical methods for solving GRPs have
been developed; see [82, Chapter 19] for a detailed description and [90] for a
comparison of such methods. Discussions on exact solutions to the classical
Riemann problem for mixtures of thermally perfect gases can be found in [91,
92].

To avoid finding the true minimum over a given element via an iterative
procedure, we relax Equation (5.17) and instead compute

$s_{b,\kappa}^{j+1}(y)=\mathsf{c}\min\left\\{s\left(y^{j}(x)\right)|x\in\mathcal{D}_{\kappa}\cup\mathcal{D}_{\kappa_{L}}\cup\mathcal{D}_{\kappa_{R}}\right\\},$
(5.20)

where $\mathsf{c}\in(0,1]$ is a relaxation parameter. The entropy limiter in
Section 5.2.1 remains valid since

$s\left(\overline{y}_{\kappa}^{j+1}\right)\geq\min\left\\{s\left(y^{j}(x)\right)|x\in\mathcal{D_{\kappa}}\cup\mathcal{D}_{\kappa}^{-}\right\\}\geq\min\left\\{s\left(y^{j}(x)\right)|x\in\mathcal{D}_{\kappa}\cup\mathcal{D}_{\kappa_{L}}\cup\mathcal{D}_{\kappa_{R}}\right\\}\geq
s_{b,\kappa}^{j+1}(y).$

For $\mathsf{c}<1$, it can be useful to simultaneously account for a well-
defined global entropy bound, $s_{b}$, as

$s_{b,\kappa}^{j+1}(y)=\max\left\\{s_{b},\mathsf{c}\min\left\\{s\left(y^{j}(x)\right)|x\in\mathcal{D}_{\kappa}\cup\mathcal{D}_{\kappa_{L}}\cup\mathcal{D}_{\kappa_{R}}\right\\}\right\\}.$
(5.21)

In this work, we simply choose $\mathsf{c}=1$. An alternative approach for
estimating the true minimum in an algebraic manner can be found in [17]. Note
that Equation (5.20) with $\mathsf{c}=1$ yields the true minimum for $p=1$ if
the element endpoints are the solution nodes and included in
$\mathcal{D}_{\kappa}$. For $\kappa=\left[x_{L},x_{R}\right]$, the basis
functions are given by

$\phi_{1}(x)=\frac{x_{R}-x}{x_{R}-x_{L}},\phi_{2}(x)=\frac{x-x_{L}}{x_{R}-x_{L}},$

such that $\phi_{1}+\phi_{2}=1$ and $\phi_{i}(x)\in(0,1),\>\forall
x\in\left(x_{L},x_{R}\right)$. Therefore,
$y_{\kappa}(x)=y_{\kappa}\left(x_{1}\right)\phi_{1}(x)+y_{\kappa}\left(x_{2}\right)\phi_{2}(x)$
is a convex combination of the endpoint values of the solution. By Lemma 9,
$\min_{x}s\left(y_{\kappa}(x)\right)=\min\left\\{s\left(y_{\kappa}\left(x_{1}\right)\right),s\left(y_{\kappa}\left(x_{2}\right)\right)\right\\}$.
This is similarly true in two and three dimensions.

####  5.2.4 Artificial viscosity

As will be demonstrated in Section 7.2, the proposed entropy-bounded DG method
does not completely suppress smaller-scale oscillations, especially in the
presence of flow-field discontinuities. Therefore, artificial viscosity is
employed in certain test cases in Section 7 to more effectively dampen such
oscillations. Specifically, the following dissipation term is added to the LHS
of Equation (3.2) [66]:

$-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{F}^{\mathrm{AV}}\left(y,\nabla
y\right),\nabla\mathfrak{v}\right)_{\kappa},$ (5.22)

where

$\mathcal{F}^{\mathrm{AV}}(y,\nabla y)=\nu_{\mathrm{AV}}\nabla y,$

with $\nu_{\mathrm{AV}}\geq 0$ denoting the artificial viscosity. The
artificial viscosity is computed as [14]

$\nu_{\mathrm{AV}}=\left(C_{\mathrm{AV}}+S_{\mathrm{AV}}\right)\left(\frac{h^{2}}{p+1}\left|\frac{\partial
T}{\partial y}\cdot\frac{\mathcal{R}\left(y,\nabla
y\right)}{T}\right|\right).$

$C_{\mathrm{AV}}$ is a user-defined coefficient, $S_{\mathrm{AV}}$ is a shock
sensor based on intra-element variations [25], and $\mathcal{R}\left(y,\nabla
y\right)$ is the strong form of the residual (2.1). Note that Equation (5.8)
(the scheme satisfied by the element averages) remains the same since the
dissipation term (5.22) vanishes for $\mathfrak{v}\in V_{h}^{0}$. Therefore,
Theorem 1 still holds. This type of artificial viscosity was found to
effectively suppress spurious oscillations in the vicinity of flow-field
discontinuities in multicomponent reacting flows [14]. However, we remark that
the artificial-viscosity formulation presented here is not the focus of this
paper; other types of artificial viscosity or limiters can be employed to
dampen the small-scale instabilities that the linear-scaling limiter fails to
cure, provided that the element-local averages are unmodified.

##  6 Reaction step: Entropy-stable discontinuous Galerkin method for ODE
integration

In this section, we describe the entropy-stable DG discretization of Equation
(3.4) (i.e., the ordinary differential equation (ODE) with stiff chemical
source terms). We build on DGODE, the (non-entropy-stable) DG method for ODE
integration described in [14]. The local, semi-discrete integral form of
Equation (3.4) is given by

$\displaystyle\int_{\kappa}\mathfrak{v}^{T}\frac{\partial y}{\partial
t}dx-\int_{\kappa}\mathfrak{v}^{T}\mathcal{S}(y)dx=0.$ (6.1)

Approximating $\mathcal{S}\left(y\right)$ locally as a polynomial in
$V_{h}^{p}$,

$\mathcal{S}_{\kappa}\approx\sum_{j=1}^{n_{b}}\mathcal{S}\left(y\left(x_{j}\right)\right)\phi_{j},$

we can write

$\frac{d}{dt}y_{\kappa}\left(x_{j},t\right)-\mathcal{S}\left(y_{\kappa}\left(x_{j},t\right)\right)=0,\quad
j=1,\ldots,n_{b},$

which is a spatially decoupled system of ODEs advanced at the solution nodes
from $t=t_{0}$ to $t=t_{f}$. Our goal here is to ensure

$\displaystyle y_{\kappa}\left(x_{j},t_{f}\right)$
$\displaystyle\in\mathcal{G}_{s\left(y_{\kappa}\left(x_{j},t_{0}\right)\right)},\quad
j=1,\ldots,n_{b}.$

Assuming a Gauss-Lobatto nodal set, since $\overline{y}_{\kappa}(t_{f})$ is a
convex combination of the nodal values, we have
$\overline{y}_{\kappa}(t_{f})\in\mathcal{G}_{s_{b}}$, where $s_{b}$ is now
given by

$s_{b}=\min_{j=1,\ldots,n_{b}}s\left(y_{\kappa}\left(x_{j},t_{0}\right)\right).$

The limiting procedure described in Section 5.2.1 can then be applied to
enforce $y_{\kappa}\left(x,t_{f}\right)\in\mathcal{G}_{s_{b}},\>\forall
x\in\mathcal{D}_{\kappa}$ (unless
$\mathcal{D}_{\kappa}=\left\\{x_{j},j=1,\ldots,n_{b}\right\\}$, in which case
the limiting procedure is unnecessary).

In the following, we drop the “$\kappa$” and “$j$” subscripts, yielding

$\frac{dy}{dt}-S(y)=0,$ (6.2)

which is the system of ODEs solved at each node. Note that the formulation
described here is slightly different from that in [12].

###  6.1 Review: DGODE

We first briefly review DGODE, referred to as “standard DGODE,” which deals
with the following one-dimensional DG discretization in time of Equation
(6.2):

$\displaystyle
N_{h}\left(y,\mathfrak{v}\right)=\sum_{\epsilon\in\mathcal{E}}\left(y^{\dagger}\left(y^{+},y^{-},n\right),\left\llbracket\mathfrak{v}\right\rrbracket\right)_{\mathcal{E}}-\sum_{\kappa\in\mathcal{T}}\left(y,\frac{d\mathfrak{v}}{dt}\right)_{\kappa}-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{S}\left(y\right),\mathfrak{v}\right)_{\kappa}=0\qquad\forall\mathfrak{v}\in
V_{h}^{p},$ (6.3)

where $\epsilon$, $\mathcal{E}$, $\kappa$, $\mathcal{T}$, and $V_{h}^{p}$ are
temporal analogs of the spatial counterparts defined in Section 3.
Specifically, $\mathcal{E}$ is the set of temporal interfaces $\epsilon$,
$\mathcal{T}$ is the set of cells $\kappa$ that partitions the computational
domain, $\Omega=\left(t_{0},t_{f}\right)=\left(t_{0},t_{0}+\Delta t\right)$,
and $V_{h}^{p}$ is the discrete subspace defined similarly to Equation (3.1).
Note that “cell” in this context corresponds to a sub-time-step (referred to
in this section as simply “time step”) of size $h\in(0,\Delta t]$. Equation
(6.3) is obtained by integrating Equation (6.2) over each cell, performing
integration by parts on the time-derivative terms, and summing over the
domain. On interior faces, the numerical temporal flux (henceforth referred to
as the “numerical state”), $y^{\dagger}$, is defined as the upwind flux
function,

$\displaystyle
y^{\dagger}\left(y^{+},y^{-},n\right)=\begin{cases}y^{+}&\textup{ if }n\geq
0\\\ -y^{-}&\textup{ if }n<0\end{cases},$ $\displaystyle\textup{ on
}\epsilon\qquad\forall\epsilon\in\mathcal{E_{I}}.$ (6.4)

On exterior interfaces, the numerical state is defined as

$\displaystyle
y^{\dagger}\left(y^{+},y^{-},n\right)=\begin{cases}y^{+}&\textup{ if }n\geq
0\\\ -y_{\partial}\left(y^{+}\right)&\textup{ if }n<0\end{cases},$
$\displaystyle\textup{ on
}\epsilon\qquad\forall\epsilon\in\mathcal{E}_{\partial},$ (6.5)

where $y_{\partial}\left(y^{+}\right)$ is a prescribed boundary state. At the
inflow interface, located at $t=t_{0}$, $y_{\partial}=y_{0}$, which is the
initial condition. At the outflow interface, located at $t=t_{f}$,
$y_{\partial}=y^{+}$ (i.e., no boundary condition is imposed).

Let $m=\dim V_{h}^{p}$ and $\left(\phi_{1},\ldots,\phi_{m}\right)$ be a basis
for $V_{h}^{p}$. The discrete residual,
$\mathcal{R}=\left(\mathcal{R}_{1}\left(y\right),\ldots,\mathcal{R}_{m}\left(y\right)\right)$
is defined as

$\mathcal{R}_{i}\left(y\right)=\sum_{\epsilon\in\mathcal{E}}\left(y^{\dagger}\left(y^{+},y^{-},n\right),\left\llbracket\phi_{i}\right\rrbracket\right)_{\mathcal{E}}-\sum_{\kappa\in\mathcal{T}}\left(y,\frac{d\phi{}_{i}}{dt}\right)_{\kappa}-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{S}\left(y\right),\phi_{i}\right)_{\kappa}$
(6.6)

for $i=1,\ldots,m$. We can then recast (6.3) as

$\mathcal{R}\left(y\right)=0,$

which is solved for $y$ via Newton’s method. With initial guess $y^{0}$, the
$k$th update, $y^{k}$, is computed by solving the linear system

$\frac{d}{dy}\mathcal{R}\left(y^{k}\right)\left(y^{k+1}-y^{k}\right)=-\mathcal{R}\left(y^{k}\right)$
(6.7)

recursively until a convergence criterion is satisfied.

Johnson and Kercher [14] introduced a simple yet effective $hp$-adaptation
strategy to control the time step, $h$, and polynomial degree, $p$. They
applied the strategy to efficiently and accurately integrate chemical systems
with complex mechanisms. Here, we focus on $h$-adaptation with uniform $p$. W
define the norm of the local error estimate as

$\mathrm{err}_{h}=\left\|\frac{\mathcal{R}(y)}{\epsilon_{\mathrm{abs}}+\epsilon_{\mathrm{rel}}\left|y\right|}\right\|,$

where $\epsilon_{\mathrm{abs}}$ and $\epsilon_{\mathrm{rel}}$ are user-
specified absolute and relative tolerances, respectively. The convergence
criterion is

$\mathrm{err}_{h}<1.$ (6.8)

If (6.8) is satisfied within a user-specified number of Newton iterations, the
solution is updated and a new time step is determined using Gustafsson’s
method [93]. Otherwise, the time step is reduced by a factor of ten.

In this work, given that implicit time stepping schemes with order greater
than one are typically not unconditionally positivity-preserving [94], we
introduce an additional convergence criterion that improves stability and more
robustly maintains conservation of mass:

$C_{i}\geq 0,\quad i=1,\ldots,n_{s}.$ (6.9)

Specifically, we require the solution at the end of each time step to satisfy
(6.9). In our experience, this additional criterion does not significantly
restrict the time step size. Furthermore, the above criterion can be relaxed.
Instead of requiring pointwise __ nonnegativity of the species concentrations,
we can simply require that the spatial averages of the concentrations over the
element be nonnnegative. The positivity-preserving limiter in Section 5.2.1
can then be applied to guarantee nonnegativity of the species concentrations
in a pointwise manner. Though not pursued in this work, the conservative and
positivity-preserving projection method by Sandu [95] can also be employed
when (6.9) is violated. Another possible convergence criterion is $T>0$, but
this is almost never a concern.

###  6.2 Entropy stability: Preliminaries

Entropy stability in this ODE setting is defined as

$U\left(y\left(t_{f}\right)\right)\leq U\left(y_{0}\right).$ (6.10)

Here, we consider the entropy function $U=-\rho s$. Combined with discrete
mass conservation (i.e., $\rho$ is constant), we have

$s\left(y\left(t_{f}\right)\right)\geq s\left(y_{0}\right).$ (6.11)

In other words, entropy stability implies that the specific thermodynamic
entropy is nondecreasing in time, which is in line with the minimum entropy
principle.

Standard DGODE, as in Section 6.1, does not necessarily satisfy the inequality
(6.10). In the following, we introduce a modified DG discretization that is
guaranteed to be entropy-stable.

###  6.3 Entropy-stable DGODE with summation-by-parts property

We work with the following strong-form discretization, obtained by performing
integration by parts on the temporal-derivative term in Equation (6.3):

$\displaystyle
N_{h}\left(y,\mathfrak{v}\right)=\sum_{\kappa\in\mathcal{T}}\left[\left(y^{\dagger}\left(y^{+},y^{-},n\right)-n\cdot
y^{+},\mathfrak{v}^{+}\right)_{\partial\kappa}+\left(\frac{dy}{dt},\mathfrak{v}\right)_{\kappa}\right]-\sum_{\kappa\in\mathcal{T}}\left(\mathcal{S}\left(y\right),\mathfrak{v}\right)_{\kappa}=0\qquad\forall\mathfrak{v}\in
V_{h}^{p}.$ (6.12)

It is well-known that a collocated DG scheme (i.e., solution nodes and
integration points are the same) with Gauss-Lobatto points possesses the
diagonal-norm summation-by-parts (SBP) property [54, 96]. Consider the
discrete mass matrix and discrete derivative matrix derived from Gauss-Lobatto
collocation, given by

$\mathsf{M}_{ij}=w_{i}\delta_{ij},\quad\mathsf{D}_{ij}=\ell^{\prime}_{j}(\xi_{i}),$

where $w_{i}$ is the $i$th Gauss-Lobatto weight, $\ell_{j}$ is the $j$th
Lagrange basis polynomial, and $\xi\in[0,1]$ is the reference coodinate.
$\mathsf{M}$ and $\mathsf{D}$ satisfy the SBP property,

$\mathsf{Q}+\mathsf{Q}^{T}=\mathsf{B},$

where $\mathsf{Q}=\mathsf{MD}$ and
$\mathsf{B}=\mathrm{diag}(-1,0,\ldots,0,1)$. The following element-local
discrete form is then obtained [18]:

$\displaystyle\left[y^{\dagger}\left(y_{\kappa}^{+},y_{\kappa}^{-},n\right)-n\cdot
y_{\kappa}^{+}\right]^{T}\mathfrak{v}_{\kappa}^{+}\biggr{|}_{\xi=0}+\left[y^{\dagger}\left(y_{\kappa}^{+},y_{\kappa}^{-},n\right)-n\cdot
y_{\kappa}^{+}\right]^{T}\mathfrak{v}_{\kappa}^{+}\biggr{|}_{\xi=1}$ (6.13)
$\displaystyle+\sum_{i=1}^{n_{b}}w_{i}\left[\sum_{j=1}^{n_{b}}\mathsf{D}_{ij}y_{\kappa}(t_{j})-h\mathcal{S}\left(y_{\kappa}(t_{i})\right)\right]^{T}\mathfrak{v}_{\kappa}(t_{i})=0.$

Note that in [14], standard DGODE was solved using a quadrature-free approach
[69, 70]. The discrete form (6.13), though similar, specifically invokes
quadrature in order to exploit the SBP property.

We replace the temporal-derivative interpolation operator in Equation (6.13)
with a specific temporal-derivative projection operator as [18]

$\sum_{j=1}^{n_{b}}\mathsf{D}_{ij}y_{\kappa}(t_{j})\rightarrow
2\sum_{j=1}^{n_{b}}\mathsf{D}_{ij}y^{\ddagger}(y_{\kappa}(t_{i}),y_{\kappa}(t_{j})),$

where $y^{\ddagger}(y_{1,}y_{2})$ is a two-point numerical state function
(distinct from $y^{\dagger}$) that is consistent and symmetric. Equation
(6.13) thus becomes

$\displaystyle\left[y^{\dagger}\left(y_{\kappa}^{+},y_{\kappa}^{-},n\right)-n\cdot
y_{\kappa}^{+}\right]^{T}\mathfrak{v}_{\kappa}^{+}\biggr{|}_{\xi=0}+\left[y^{\dagger}\left(y_{\kappa}^{+},y_{\kappa}^{-},n\right)-n\cdot
y_{\kappa}^{+}\right]^{T}\mathfrak{v}_{\kappa}^{+}\biggr{|}_{\xi=1}$ (6.14)
$\displaystyle+\sum_{i=1}^{n_{b}}w_{i}\left[2\sum_{j=1}^{n_{b}}\mathsf{D}_{ij}y^{\ddagger}(y_{\kappa}(t_{i}),y_{\kappa}(t_{j}))-h\mathcal{S}\left(y_{\kappa}(t_{i})\right)\right]^{T}\mathfrak{v}_{\kappa}(t_{i})=0.$

With the mean-value numerical state,

$y^{\ddagger}(y_{L},y_{R})=\frac{y_{L}+y_{R}}{2},$ (6.15)

Equation (6.14) recovers Equation (6.13) [96, 54].

The discretization (6.14) combined with the mean-value numerical state is not
necessarily entropy-stable. Although, as will be shown in Section 6.3.2, the
source-term discretization results in destruction of mathematical entropy, the
temporal-derivative operator may cause production of mathematical entropy.
Entropy stability is achieved only if the entropy destruction due to the
source term outweighs the entropy production due to the temporal-derivative
operator. In the following, by replacing the mean-value numerical state with a
more appropriate numerical state, we will obtain a temporal DG scheme that is
guaranteed to be entropy-stable.

####  6.3.1 Entropy-conservative numerical state function

A key ingredient in the development of an entropy-stable DG scheme is a two-
point numerical state function that satisfies the following condition [18]:

$(\mathsf{v}_{R}-\mathsf{v}_{L})^{T}y^{\ddagger}(y_{L},y_{R})=\mathcal{U}_{R}-\mathcal{U}_{L}.$
(6.16)

A numerical state that satisfies (6.16) is _entropy-conservative_. Note that
an entropy-conservative numerical _flux_ satisfies an analogous condition:

$(\mathsf{v}_{R}-\mathsf{v}_{L})^{T}\mathcal{F}^{\ddagger}(y_{L},y_{R})=\mathcal{F}_{R}^{p}-\mathcal{F}_{L}^{p}.$
(6.17)

Gouasmi et al. [51] derived a simple, closed-form entropy-conservative
numerical flux for the multicomponent Euler equations with the entropy
function $U=-\rho s$. They built on the techniques by Roe [97] originally used
to construct an entropy-conservative numerical flux for the monocomponent
Euler equations. Said techniques rely on (6.17) as a starting point. An
analogous procedure, instead using (6.16) as a starting point, is employed
here to derive an entropy-conservative numerical state function.

With $\left\llbracket\cdot\right\rrbracket$ denoting the jump operator, the
entropy conservation condition (6.16) can be expressed as

$\left\llbracket\mathsf{v}\right\rrbracket{}^{T}y^{\ddagger}=\left\llbracket\mathcal{U}\right\rrbracket.$
(6.18)

For consistency with [51], we first work with a re-ordered state vector where
the species concentrations are replaced with partial densities:

$y=\left(\rho_{1},\ldots,\rho_{n_{s}},\rho v_{1},\ldots,\rho v_{d},\rho
e_{t}\right)^{T}.$ (6.19)

We employ the notation

$y^{\ddagger}=\left(y_{1,1}^{\ddagger},\ldots,y_{1,n_{s}}^{\ddagger},y_{2,1}^{\ddagger},\ldots,y_{2,d}^{\ddagger},y_{3}^{\ddagger}\right),$

where the first $n_{s}$ components correspond to the partial densities, the
next $d$ components correspond to the momentum, and the last component
corresponds to the total energy. Let $z$ denote the vector

$z=\left(\rho_{1},\ldots,\rho_{n_{s}},v_{1},\ldots,v_{d},1/T\right)=\left(z_{1,1},\ldots,z_{1,n_{s}},z_{2,1},\ldots,z_{2,d},z_{3}\right).$

With the entropy function $U=-\rho s$, the entropy variables and entropy
potential are given by

$\mathsf{v}=\left(\frac{g_{1}-\frac{1}{2}\sum_{k=1}^{d}v_{k}v_{k}}{T},\ldots,\frac{g_{n_{s}}-\frac{1}{2}\sum_{k=1}^{d}v_{k}v_{k}}{T},\frac{v_{1}}{T},\ldots,\frac{v_{d}}{T},-\frac{1}{T},\right)^{T},\quad\mathcal{U}=\sum_{i=1}^{n_{s}}W_{i}\rho_{i}.$

Next, we introduce the arithmetic mean, logarithmic mean, and product
operators,

$\displaystyle\left\\{\\!\\!\left\\{\alpha\right\\}\\!\\!\right\\}$
$\displaystyle=\frac{\alpha_{L}+\alpha_{R}}{2},$ $\displaystyle\alpha^{\ln}$
$\displaystyle=\begin{cases}\alpha_{L},&\text{if }\alpha_{L}=\alpha_{R}\\\
0,&\text{if }\alpha_{L}=0\text{ or }\alpha_{R}=0\\\
\frac{\alpha_{R}-\alpha_{L}}{\ln\alpha_{R}-\ln\alpha_{L}},&\text{otherwise},\end{cases}$
(6.20) $\displaystyle\alpha^{\times}$ $\displaystyle=\alpha_{L}\alpha_{R}$

which are equipped with the identities

$\displaystyle\left\\{\\!\\!\left\\{\alpha\beta\right\\}\\!\\!\right\\}$
$\displaystyle=\alpha\left\llbracket\beta\right\rrbracket+\beta\left\llbracket\alpha\right\rrbracket$
$\displaystyle\left\llbracket\ln\alpha\right\rrbracket$
$\displaystyle=\frac{\left\llbracket\alpha\right\rrbracket}{\alpha^{\ln}}$
$\displaystyle\left\llbracket\alpha\right\rrbracket$
$\displaystyle=-\alpha^{\times}\left\llbracket\frac{1}{\alpha}\right\rrbracket,\alpha\neq
0.$

Note that the three operators in (6.20) are symmetric; the first two are also
consistent (the product operator is consistent with $\alpha^{2}$). To compute
the logarithmic mean in a numerically stable manner, we employ the procedure
by Ismail and Roe [98]. The jump in the entropy potential can then be written
as

$\left\llbracket\mathcal{U}\right\rrbracket=\sum_{i=1}^{n_{s}}W_{i}\left\llbracket\rho_{i}\right\rrbracket=\sum_{i=1}^{n_{s}}W_{i}\left\llbracket
z_{1,i}\right\rrbracket,$ (6.21)

while the jump in the entropy variables, after some algebraic manipulation,
can be expressed as [51]

$\left\llbracket\mathsf{v}\right\rrbracket=\begin{pmatrix}\left\llbracket
z_{3}\right\rrbracket\mathsf{X}_{1}+\left\llbracket
z_{1,1}\right\rrbracket\frac{W_{1}}{\rho_{1}^{\ln}}-\sum_{k=1}^{d}\left\llbracket
z_{2,k}\right\rrbracket\left\\{\\!\\!\left\\{\frac{1}{T}\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}\\\
\vdots\\\ \left\llbracket
z_{3}\right\rrbracket\mathsf{X}_{n_{s}}+\left\llbracket
z_{1,n_{s}}\right\rrbracket\frac{W_{n_{s}}}{\rho_{n_{s}}^{\ln}}-\sum_{k=1}^{d}\left\llbracket
z_{2,k}\right\rrbracket\left\\{\\!\\!\left\\{\frac{1}{T}\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}\\\
\left\\{\\!\\!\left\\{z_{3}\right\\}\\!\\!\right\\}\left\llbracket
z_{2,1}\right\rrbracket+\left\\{\\!\\!\left\\{z_{2,1}\right\\}\\!\\!\right\\}\left\llbracket
z_{3}\right\rrbracket\\\ \vdots\\\
\left\\{\\!\\!\left\\{z_{3}\right\\}\\!\\!\right\\}\left\llbracket
z_{2,d}\right\rrbracket+\left\\{\\!\\!\left\\{z_{2,d}\right\\}\\!\\!\right\\}\left\llbracket
z_{3}\right\rrbracket\\\ -\left\llbracket z_{3}\right\rrbracket\end{pmatrix},$
(6.22)

where

$\mathsf{X}_{i}=b_{i0}+\frac{b_{i1}}{\left(1/T\right)^{\ln}}+\sum_{r=2}^{n_{p}+1}b_{ir}\left(\mathsf{f}_{r-1}(T)T^{\times}\right)-\frac{1}{2}\sum_{k=1}^{d}\left\\{\left\\{v_{d}^{2}\right\\}\right\\},$

with $\mathsf{f}_{r}(\alpha)$ denoting a special averaging operator,
consistent with $\alpha^{r-1}$, that satisfies
$\left\llbracket\alpha^{r}\right\rrbracket=r\mathsf{f}_{r}(\alpha)\left\llbracket\alpha\right\rrbracket$.
For $r=1,2,3,4$, $\mathsf{f}_{r}(\alpha)$ is defined as

$\displaystyle\mathsf{f}_{1}(\alpha)$ $\displaystyle=1,$
$\displaystyle\mathsf{f}_{2}(\alpha)$
$\displaystyle=\left\\{\\!\\!\left\\{\alpha\right\\}\\!\\!\right\\},$
$\displaystyle\mathsf{f}_{3}(\alpha)$
$\displaystyle=\frac{2}{3}\left\\{\\!\\!\left\\{\alpha\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{\alpha\right\\}\\!\\!\right\\}+\frac{1}{3}\left\\{\\!\\!\left\\{\alpha^{2}\right\\}\\!\\!\right\\},$
$\displaystyle\mathsf{f}_{4}(\alpha)$
$\displaystyle=\left\\{\\!\\!\left\\{\alpha\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{\alpha^{2}\right\\}\\!\\!\right\\}.$

$\mathsf{f}_{r}(\alpha)$ can be derived for $r>4$ as well [51]. Using
Equations (6.21) and (6.22), we rewrite the entropy conservation condition in
time (6.18) as a requirement that a linear combination of the jumps in the
components of $z$ equals zero:

$\displaystyle\sum_{i=1}^{n_{s}}\left\llbracket
z_{1,i}\right\rrbracket\left(y_{1,i}^{\ddagger}\frac{W_{i}}{\rho_{i}^{\ln}}-W_{i}\right)+\sum_{k=1}^{d}\left\llbracket
z_{2,k}\right\rrbracket\left(y_{2,k}^{\ddagger}\left\\{\\!\\!\left\\{z_{3}\right\\}\\!\\!\right\\}-\sum_{i=1}^{n_{s}}y_{1,i}^{\ddagger}\left\\{\\!\\!\left\\{\frac{1}{T}\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}\right)$
$\displaystyle+\left\llbracket
z_{3}\right\rrbracket\left(\sum_{i=1}^{n_{s}}y_{1,i}^{\ddagger}\mathsf{X}_{i}+\sum_{k=1}^{d}y_{2,k}^{\ddagger}\left\\{\\!\\!\left\\{z_{2,k}\right\\}\\!\\!\right\\}-y_{3}^{\ddagger}\right)=0.$

Invoking the independence of these jumps yields a system of $m$ equations:

$\displaystyle y_{1,i}^{\ddagger}\frac{W_{i}}{\rho_{i}^{\ln}}-W_{i}$
$\displaystyle=0,\quad i=1,\ldots,n_{s},$ $\displaystyle
y_{2,k}^{\ddagger}\left\\{\\!\\!\left\\{z_{3}\right\\}\\!\\!\right\\}-\sum_{i=1}^{n_{s}}y_{1,i}^{\ddagger}\left\\{\\!\\!\left\\{\frac{1}{T}\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}$
$\displaystyle=0,\quad k=1,\ldots,d,$
$\displaystyle\sum_{i=1}^{n_{s}}y_{1,i}^{\ddagger}\mathsf{X}_{i}+\sum_{k=1}^{d}y_{2,k}^{\ddagger}\left\\{\\!\\!\left\\{z_{2,k}\right\\}\\!\\!\right\\}-y_{3}^{\ddagger}$
$\displaystyle=0.$

Solving for the components of $y^{\ddagger}$ then gives

$y^{\ddagger}=\begin{pmatrix}\rho_{1}^{\ln}\\\ \vdots\\\ \rho_{n_{s}}^{\ln}\\\
\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left\\{\\!\\!\left\\{v_{1}\right\\}\\!\\!\right\\}\\\
\vdots\\\
\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left\\{\\!\\!\left\\{v_{d}\right\\}\\!\\!\right\\}\\\
\sum_{k=1}^{d}\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}+\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left[b_{i0}+\frac{b_{i1}}{\left(1/T\right)^{\ln}}+\sum_{r=2}^{n_{p}+1}b_{ir}\left(\mathsf{f}_{r-1}(T)T^{\times}\right)-\frac{1}{2}\sum_{k=1}^{d}\left\\{\\!\\!\left\\{v_{k}^{2}\right\\}\\!\\!\right\\}\right]\end{pmatrix}.$
(6.23)

Note that Equation (6.23) corresponds to the state vector (6.19). A simple re-
ordering and linear mapping [51] yields

$y^{\ddagger}=\begin{pmatrix}\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left\\{\\!\\!\left\\{v_{1}\right\\}\\!\\!\right\\}\\\
\vdots\\\
\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left\\{\\!\\!\left\\{v_{d}\right\\}\\!\\!\right\\}\\\
\sum_{k=1}^{d}\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}\left\\{\\!\\!\left\\{v_{k}\right\\}\\!\\!\right\\}+\sum_{i=1}^{n_{s}}\rho_{i}^{\ln}\left[b_{i0}+\frac{b_{i1}}{\left(1/T\right)^{\ln}}+\sum_{r=2}^{n_{p}+1}b_{ir}\left(\mathsf{f}_{r-1}(T)T^{\times}\right)-\frac{1}{2}\sum_{k=1}^{d}\left\\{\\!\\!\left\\{v_{k}^{2}\right\\}\\!\\!\right\\}\right]\\\
C_{1}^{\ln}\\\ \vdots\\\ C_{n_{s}}^{\ln}\end{pmatrix},$ (6.24)

which corresponds to the original state vector (2.2). We then have the
following theorem.

######  Theorem 3.

The two-point numerical state function in Equation (6.24) is entropy-
conservative, consistent, and symmetric.

###### Proof.

The numerical state (6.24) is entropy-conservative by construction. Since all
of the introduced operators are symmetric, the numerical state is symmetric.
Finally, recognizing that $\mathsf{f}_{r-1}(\alpha)\alpha^{\times}$ is
consistent with $\alpha^{r}$, taking the left and right states to be the same
yields

$y^{\ddagger}=\begin{pmatrix}\rho v_{1}\\\ \vdots\\\ \rho v_{d}\\\
\sum_{k=1}^{d}\rho
v_{k}v_{k}+\sum_{i=1}^{n_{s}}\rho_{i}\left[b_{i0}+b_{i1}T+\sum_{r=2}^{n_{p}+1}b_{ir}T^{r}-\frac{1}{2}\sum_{k=1}^{d}v_{k}^{2}\right]\\\
C_{1}\\\ \vdots\\\ C_{n_{s}}\end{pmatrix}=\begin{pmatrix}\rho v_{1}\\\
\vdots\\\ \rho v_{d}\\\ \rho u+\frac{1}{2}\rho\sum_{k=1}^{d}v_{k}v_{k}\\\
C_{1}\\\ \vdots\\\ C_{n_{s}}\end{pmatrix}=y.$

Therefore, it is also consistent. ∎

######  Remark 4.

The derived entropy-conservative numerical state function (6.24) can be
directly used in entropy-stable space-time DG schemes based on SBP operators
[18].

####  6.3.2 Discrete temporal entropy analysis

We analyze the entropy stability of the proposed temporal DG discretization in
the following theorem.

######  Theorem 5.

Consider the DG discretization (6.14). Assume that $y^{\dagger}$ is the upwind
numerical state and $y^{\ddagger}$ is the entropy-conservative numerical state
in (6.24). Then, for the entropy function $U=-\rho s$, the resulting DG
discretization is entropy-stable.

###### Proof.

Take $\mathfrak{v}_{\kappa}$ in (6.14) to be the polynomial interpolant of the
entropy variables, such that

$\mathfrak{v}_{\kappa}\left(t_{j}\right)=\mathsf{v}\left(y_{\kappa}\left(t_{j}\right)\right),\quad
j=1,\ldots,n_{b},$

which results in

$\begin{split}&\sum_{\kappa\in\mathcal{T}}\left[\left[y^{\dagger}\left(y_{\kappa}^{+},y_{\kappa}^{-},n\right)-n\cdot
y_{\kappa}^{+}\right]^{T}\mathsf{v}\left(y_{\kappa}^{+}\right)\biggr{|}_{\xi=0}+\left[y^{\dagger}\left(y_{\kappa}^{+},y_{\kappa}^{-},n\right)-n\cdot
y_{\kappa}^{+}\right]^{T}\mathsf{v}\left(y_{\kappa}^{+}\right)\biggr{|}_{\xi=1}\right]\\\
&+\sum_{\kappa\in\mathcal{T}}\left[\sum_{i=1}^{n_{b}}w_{i}\left(2\sum_{j=1}^{n_{b}}\mathsf{D}_{ij}y^{\ddagger}(y_{\kappa}(t_{i}),y_{\kappa}(t_{j}))-h\mathcal{S}\left(y_{\kappa}(t_{i})\right)\right)^{T}\mathsf{v}\left(y_{\kappa}(t_{i})\right)\right]=0.\end{split}$
(6.25)

We introduce $\mathcal{A}_{\kappa}$ and $\mathcal{B}_{\kappa}$, defined as

$\displaystyle\mathcal{A}_{\kappa}=$
$\displaystyle\left[y^{\dagger}\left(y^{+},y^{-},n\right)-n\cdot
y^{+}\right]^{T}\mathsf{v}\left(y_{\kappa}^{+}\right)\biggr{|}_{\xi=0}+\left[y^{\dagger}\left(y^{+},y^{-},n\right)-n\cdot
y^{+}\right]^{T}\mathsf{v}\left(y_{\kappa}^{+}\right)\biggr{|}_{\xi=1}$
$\displaystyle+\sum_{i=1}^{n_{b}}w_{i}\left[2\sum_{j=1}^{n_{b}}\mathsf{D}_{ij}y^{\ddagger}(y(t_{i}),y(t_{j}))\right]^{T}\mathsf{v}\left(y_{\kappa}(t_{i})\right)$
$\displaystyle\mathcal{B}_{\kappa}=$
$\displaystyle-\sum_{i=1}^{n_{b}}w_{i}hS\left(y_{\kappa}(t_{i})\right)^{T}\mathsf{v}\left(y_{\kappa}(t_{i})\right),$

such that Equation (6.25) can be rewritten as

$\sum_{\kappa\in\mathcal{T}}\mathcal{A}_{\kappa}+\sum_{\kappa\in\mathcal{T}}\mathcal{B}_{\kappa}=0.$
(6.26)

By invoking the SBP property and the fact that the upwind numerical state
function is an entropy-stable numerical state (i.e., it satisfies
$\left\llbracket\mathsf{v}\right\rrbracket{}^{T}y^{\dagger}\leq\left\llbracket\mathcal{U}\right\rrbracket$)
[18], we obtain the inequality

$U\left(y\left(t_{f}\right)\right)-U\left(y_{0}\right)\leq\sum_{\kappa\in\mathcal{T}}\mathcal{A}_{\kappa},$
(6.27)

the proof of which is very similar to that in [18, Theorem 1] (just without
any spatial component) and is therefore not included here. It remains to
analyze $\sum_{\kappa\in\mathcal{T}}\mathcal{B}_{\kappa}$. The quantity
$S\left(y_{\kappa}(t_{i})\right)^{T}\mathsf{v}\left(y_{\kappa}(t_{i})\right)$
is simply the pointwise entropy production rate due to the chemical source
terms, which, as demonstrated in Section 4.2, is nonpositive. Since $w_{i}>0$
and $h>0$, we have

$\sum_{\kappa\in\mathcal{T}}\mathcal{B}_{\kappa}\geq 0.$ (6.28)

Combining Equations (6.26), (6.27), and (6.28) gives the inequality (6.10),
which completes the proof. ∎

In the following, we refer to the proposed entropy-stable DG discretization as
simply “entropy-stable DGODE.”

######  Remark 6.

Entropy-stable DGODE is well-defined for zero concentrations. However, during
early iterations of Newton’s method, negative concentrations can occur, making
the formulation ill-defined. As a simple remedy, we use the solution obtained
with the mean-value numerical state (6.15) as an initial guess, which we find
to be sufficiently robust since this initial guess is generally close to the
entropy-stable DGODE solution (recall that taking $y^{\ddagger}$ to be the
mean-value numerical state essentially recovers standard DGODE). Other, more
sophisticated approaches can be employed as well.

######  Remark 7.

In general, entropy-stable DGODE is more expensive than standard DGODE.
Furthermore, when using the latter, we find that the thermodynamic entropy
produced by the source terms typically outweighs any thermodynamic entropy
destruction caused by the non-entropy-stable mean-value numerical state, at
least for the reaction mechanisms used in this study (though it is important
to note that this observation may not hold true for significantly different
mechanisms). As such, in practice, the following approach can be employed to
maximize efficiency:

  * 1.

Compute a solution with standard DGODE.

  * 2.

If $s\left(y\left(t_{f}\right)\right)\geq s\left(y_{0}\right)$, then proceed.
Otherwise, compute a solution with entropy-stable DGODE.

In Section 7.3, however, in order to test the formulation, we calculate
spatially one-dimensional detonation waves with entropy-stable DGODE alone.

######  Remark 8.

To reduce computational cost, it is common practice, especially for large-
scale simulations, to approximate a given chemical reaction as two
irreversible forward reactions (as opposed to having forward and reverse
reaction rates). However, this strategy can cause appreciable entropy
violations [60, 99]. A detailed investigation of whether the gains in speed
outweigh the loss in accuracy is outside the scope of this study. If such a
strategy is employed, then reverting to standard DGODE is the natural course
of action. Though not considered in this work, an alternative approach is to
use the least-squares-based method by Ream et al. [100] to generate chemical
mechanisms involving irreversible reactions that do not cause entropy
violations.

##  7 One-dimensional results

We consider three one-dimensional test cases. The first one involves the
advection of a thermal bubble. The second case is a shock-tube problem with
multiple flow discontinuities. These first two tests comprise nonreacting
multicomponent flows. The final case explores sustained detonations formed via
an overdriven initialization. Stiff chemical reactions are present in this
test. The SSPRK2 time integration scheme with $\mathrm{CFL}=0.1$ (based on the
linear-stability constraint) is employed throughout. All simulations are
performed using a modified version of the JENRE® Multiphysics Framework [101,
14] that incorporates the developments and extensions described in this work.

###  7.1 Thermal bubble

We use this smooth flow problem to assess the grid convergence of the entropy-
bounded DG method (without artificial viscosity). The order of accuracy of the
limiting procedure in Section 5.2.1 is of particular interest. We also compare
the local entropy bound in [17] to that proposed here (Equation (5.20)). The
initial conditions are as follows:

$\displaystyle v_{1}$ $\displaystyle=$ $\displaystyle 1\textrm{ m/s},$
$\displaystyle Y_{H_{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[1-\tanh\left(|x|-10\right)\right],$
$\displaystyle Y_{O_{2}}$ $\displaystyle=$ $\displaystyle 1-Y_{H_{2}},$ (7.1)
$\displaystyle T$ $\displaystyle=$ $\displaystyle
1200-900\tanh\left(|x|-10\right)\textrm{ K},$ $\displaystyle P$
$\displaystyle=$ $\displaystyle 1\textrm{ bar}.$

The computational domain is $\Omega=[-25,25]\>\mathrm{m}$, with periodic
boundaries. Four element sizes are considered: $h$, $h/2$, $h/4$, and $h/8$,
where $h=2\>\mathrm{m}$. Smaller element sizes are not investigated since the
limiters are not activated for such fine resolutions and optimal order of
accuracy without limiting was already demonstrated in [14]. The $L^{2}$ error
at $t=5\>\mathrm{s}$ is calculated in terms of the following normalized state
variables:

$\widehat{\rho v}_{k}=\frac{1}{\sqrt{\rho_{r}P_{r}}}\rho
v_{k},\quad\widehat{\rho e}_{t}=\frac{1}{P_{r}}\rho
e_{t},\quad\widehat{C}_{i}=\frac{R^{0}T_{r}}{P_{r}}C_{i},$

where $T_{r}=1000\,\mathrm{K}$, $\rho_{r}=1\,\mathrm{kg\cdot}\mathrm{m}^{-3}$,
and $P_{r}=101325\,\mathrm{Pa}$ are reference values. The results of the
convergence study for $p=1$ to $p=3$ are displayed in Figure 7.1. The dashed
lines represent the theoretical convergence rates. The “$\times$” marker
indicates that the positivity-preserving limiter is activated, the “$\Circle$”
marker indicates that the entropy limiter is activated, and the “$\triangle$”
marker indicates that neither limiter is activated. If both limiters are
activated, then the corresponding markers are superimposed as “$\otimes$”. For
$h$ and $h/2$, both the positivity-preserving and entropy limiters are engaged
regardless of $p$; for $h/4$ and $p=1$, only the positivity-preserving limiter
is engaged. That the limiters are no longer activated for well-resolved
solutions is a desirable property. In general, optimal order of accuracy is
recovered. Suboptimal accuracy is observed for the coarser resolutions with
$p=1$, likely because the asymptotic regime is not yet reached.

Figure 7.1: Convergence under grid refinement, with $h=2\>\mathrm{m}$, for the
one-dimensional thermal bubble test case. The $L^{2}$ error of the normalized
state with respect to the exact solution at $t=5\>\mathrm{s}$ is computed. The
dashed lines represent the theroretical convergence rates. The “$\times$”
marker indicates that the positivity-preserving limiter is activated, the
“$\Circle$” marker indicates that the entropy limiter is activated, and the
“$\triangle$” marker indicates that neither limiter is activated. If both
limiters are activated, then the corresponding markers are superimposed as
“$\otimes$”.

Figure 7.2 compares the temperature profile at $t=50\>s$ computed with the
local entropy bound in [17] (“Old”) to that computed with the local entropy
bound in Equation (5.20) (“New”). The exact solution is the same as the
initial condition. The element size is $h/4$ and the polynomial order is
$p=1$. The reason for choosing $p=1$ is to eliminate any ambiguity associated
with algebraically estimating the spatial minimum of the specific
thermodynamic entropy, as discussed in Section 5.2.3. Though oscillations are
present in both cases, they are significantly larger with the local entropy
bound in [17], suggesting that it may be overly restrictive. The local Lax-
Friedrichs flux is employed in this comparison. With the selected parameters,
the differences between the two local entropy bounds are more pronounced for
this numerical flux than the HLLC flux. Nevertheless, such discrepancies are
observed for the latter flux function in other configurations as well,
especially with higher polynomial orders.

Figure 7.2: Temperature profiles at $t=50\>s$ computed with the local entropy
bound in [17] (“Old”) and the local entropy bound in Equation (5.20) (“New”).
The exact solution is the same as the initial condition. The element size is
$h/4$ and the polynomial order is $p=1$.

###  7.2 Shock tube

This test case was first presented by Houim and Kuo [102] and computed as well
by Johnson and Kercher [14]. The goals here are to compare the stabilization
properties of the limiters and to illustrate the benefits of the local entropy
bound. The initial conditions are given by

$\left(v_{1},T,P,Y_{N_{2}},Y_{He}\right)=\begin{cases}\left(0\text{
m/s},300\text{ K},1\text{ atm},1,0\right),&x\geq 0.4\\\ \left(0\text{
m/s},300\text{ K},10\text{ atm},0,1\right),&x<0.4\end{cases}.$ (7.2)

The computational domain is $\Omega=[0,1]\text{ m}$, with walls at both ends.
Figure 7.3 shows the mass fraction, pressure, temperature, and entropy
profiles at $t=300\;\mu\mathrm{s}$ for $p=3$ and 200 elements. We present
results for only the positivity-preserving limiter (referred to as “PPL”) and
for both the positivity-preserving and entropy limiters with the local entropy
bound in Equation (5.20) (referred to as “local EB”). Also included is a
reference solution computed with $p=2$, 2000 elements, and artificial
viscosity, which corresponds to the configuration in [14]. Artificial
viscosity is not employed in the coarser cases in order to isolate the effects
of the limiters. The mass fraction profiles are well-captured for both types
of limiting. However, although the solver does not crash, the positivity-
preserving limiter by itself fails to suppress large-scale oscillations and
significant overshoots/undershoots in the pressure and temperature
distributions. While instabilities are still present with the local entropy
limiter, they are of substantially smaller magnitude. The entropy distribution
obtained with the local entropy limiter is very similar to that of the
reference solution, whereas the positivity-preserving limiter generates a
notable overshoot and undershoot at the shock. For the given flow conditions,
the global entropy bound yields very similar results (not shown for brevity)
to the local entropy bound since the specific thermodynamic entropy in the
vicinity of the shock, which is where much of the limiting occurs, is close to
the global minimum.

The instabilities observed in the “PPL” case in Figures 7.3b and 7.3c are
considerably larger than those typically observed in shock-tube solutions
obtained with the positivity-preserving limiter in the monocomponent,
calorically perfect case [20, 36, 35]. This difference reflects the numerical
challenges associated with not only multicomponent mixtures, but also variable
thermodynamics. In a similar vein, the relative benefit of the entropy limiter
(compared to the positivity-preserving limiter) seems significantly greater in
the multicomponent, thermally perfect case than in the monocomponent,
calorically perfect case.

(a) Mass fractions. (b) Pressure. (c) Temperature. (d) Specific thermodynamic
entropy.

Figure 7.3:  Results for $p=3$ solutions on 200 elements without artificial
viscosity for the one-dimensional, multicomponent shock-tube problem with
initialization in Equation (7.2). “PPL” corresponds to the positivity-
preserving limiter by itself, and “local EB” refers to both the positivity-
preserving and entropy limiters with the local entropy bound in Equation
(5.20). The reference solution [14] is computed with $p=2$, 2000 elements, and
artificial viscosity.

Figure 7.4 presents the percent error in conservation of mass, energy, and
atomic elements for the “local EB” case as a representative example,
calculated every $0.3\;\mu\mathrm{s}$ (for a total of 1000 samples).
$\mathsf{N}_{N}$ and $\mathsf{N}_{He}$ denote the total numbers of nitrogen
and helium atoms in the mixture. The error remains close to machine precision,
confirming that the proposed formulation is conservative. Also included is the
error in mass conservation (calculated every time step) for a solution
computed without the positivity-preserving and entropy limiters, but instead
with a simple clipping procedure in which negative species concentrations are
set to zero, a strategy employed by many reacting-flow solvers. The error
increases rapidly to non-negligible values until the solver diverges.

Figure 7.4: Percent error in conservation of mass, energy, and atomic elements
for the “local EB” case in Figure 7.3, computed with $p=3$ on 200 elements.
The initial conditions for this one-dimensional, multicomponent shock-tube
problem are given in Equation (7.2). Also included is the error in mass
conservation for a solution computed without the positivity-preserving and
entropy limiters, but instead with a simple clipping procedure in which
negative species concentrations are set to zero.

Next, we recompute this problem with artificial viscosity to confirm adequate
suppression of small-scale oscillations. Figure 7.5 presents the results for
$C_{\mathrm{AV}}=1$. The instabilities observed in Figure 7.3 are largely
eliminated by the artificial viscosity. As shown in Figure 7.5c, a temperature
undershoot at the shock emerges when only the positivity-preserving limiter is
used, but is suppressed by the entropy limiter. The artificial viscosity
causes some smearing of the solution at the contact. Note that without any
limiting, negative species concentrations occur for the $p=3$, 200-element
cases, even with artificial viscosity. As will be further discussed in Part II
[19], artificial viscosity alone, or even when combined with solely the
positivity-preserving limiter, does not provide sufficient stabilization in
simulations of complex detonation waves on relatively coarse meshes. In these
simulations, enforcement of the entropy principle is critical for robustness.

(a) Mass fractions. (b) Pressure. (c) Temperature. (d) Specfic thermodynamic
entropy.

Figure 7.5:  Results for $p=3$ solutions on 200 elements with artificial
viscosity for the one-dimensional, multicomponent shock-tube problem with
initialization in Equation (7.2). “PPL” corresponds to the positivity-
preserving limiter by itself, and “local EB” refers to both the positivity-
preserving and entropy limiters with the local entropy bound in Equation
(5.20). The reference solution [14] is computed with $p=2$, 2000 elements, and
artificial viscosity. The difference between these results and those in Figure
7.3 is the use of artificial viscosity in the (non-reference) solutions here.

To highlight discrepancies between the local and global entropy bounds, we
consider different initial conditions, given by

$\left(v_{1},T,P,Y_{N_{2}},Y_{He}\right)=\begin{cases}\left(0\text{
m/s},300\text{ K},1\text{ atm},0,1\right),&x\geq 0.4\\\ \left(0\text{
m/s},300\text{ K},10\text{ atm},1,0\right),&x<0.4\end{cases}.$ (7.3)

The only difference with the previous initial conditions is in the mass
fractions. Displayed in Figure 7.6 are the results obtained with the global
entropy bound (referred to as “global EB”) and the local entropy bound (again
referred to as “local EB”). These $p=2$ solutions are computed on 200 elements
without any artificial viscosity. The positivity-preserving limiter by itself
yields very similar results (not shown for brevity) to the entropy limiter
with the global entropy bound. The reference solution is again computed with
$p=2$, 2000 elements, and artificial viscosity. The mass fractions are well-
captured in all cases. The differences between the “global EB” and “local EB”
solutions here are not as large as those between the “PPL” and “local EB”
solutions in Figure 7.3. Nevertheless, the benefit of the local entropy bound
is evident, specifically at the shock. Spurious artifacts in the pressure
profile and especially the temperature profile near the shock are noticeably
larger for the global entropy bound. The discrepancies between the two
solutions are attributed to the considerable difference between the specific
thermodynamic entropy in the vicinity of the shock and the global minimum, as
illustrated in Figure 7.6d. As such, the local entropy bound is particularly
beneficial in flow problems with large variations in the specific
thermodynamic entropy throughout the domain.

(a) Mass fractions. (b) Pressure. (c) Temperature. (d) Specfic thermodynamic
entropy.

Figure 7.6:  Results for $p=2$ solutions on 200 elements without artificial
viscosity for the one-dimensional, multicomponent shock-tube problem with
initialization in Equation (7.3). “Global EB” refers to the entropy limiter
with the global entropy bound in Equation (5.15), and “local EB” refers to the
entropy limiter with the local entropy bound in Equation (5.20). The reference
solution is computed with $p=2$, 2000 elements, and artificial viscosity.

###  7.3 Detonation wave

The final one-dimensional test case is a hydrogen-oxygen detonation wave
diluted in Argon with initial conditions

$\begin{array}[]{cccc}\qquad\qquad\qquad\qquad v_{1}&=&0\text{ m/s},\\\
\quad\quad\quad\;X_{Ar}:X_{O_{2}}:X_{H_{2}}&=&7:1:2&\text{ }x>0.025\text{
m},\\\ X_{Ar}:X_{O_{2}}:X_{H_{2}}:X_{OH}&=&7:1:2:0.01&\text{ }0.015\text{
m}<x<0.025\text{ m,}\\\
\quad\quad\;X_{Ar}:X_{H_{2}O}:X_{OH}&=&8:2:0.0001&x<0.015\text{ m},\\\
\qquad\qquad\qquad\qquad P&=&\begin{cases}{5.50}\mathrm{e}{5}&\text{ Pa}\\\
{6.67}\mathrm{e}{3}&\text{ Pa}\end{cases}&\begin{array}[]{c}x<0.015\text{
m}\\\ x>0.015\text{ m}\end{array},\\\ \qquad\qquad\qquad\qquad
T&=&\begin{cases}300&\text{ K}\\\ 350&\text{ K}\\\ 3500\text{\hskip
10.00002pt\hskip 10.00002pt}&\text{ K}\end{cases}&\begin{array}[]{c}\text{
}x>0.025\text{ m}\\\ \text{ }0.015\text{ m}<x<0.025\text{ m}\\\ x<0.015\text{
m}\end{array}.\end{array}$ (7.4)

The domain is $\Omega=(0,0.45)$ m, with walls at the left and right
boundaries. In previous work [14], this case was computed with $p=1$ and mesh
spacing $h=9\times 10^{-5}$ m. At this resolution, artificial viscosity was
the only stabilization necessary to obtain an accurate solution while
maintaining conservation of mass and energy. Good agreement with the Shock and
Detonation Toolbox [103] was observed. Additional details can be found in
[14]. Here, our objective is to demonstrate that the proposed formulation,
specifically the entropy-bounded DG discretization of the convective operator
combined with entropy-stable DGODE for stiff chemical reactions, can compute
stable and accurate solutions with appreciably lower resolution. In light of
[14], we use a $p=1$, $h=9\times 10^{-5}$ calculation as a reference solution.
The only difference with [14] is that the reaction mechanism [104] is slightly
modified to solely contain reversible reactions.

Figure 7.7 presents $p=1$, $p=2$, and $p=3$ results at $t=235$
$\mu\mathrm{s}$. The mesh spacing in these simulations is $h=4.5\times
10^{-4}$ m, which is fives time larger than for the reference solution.
Artificial viscosity and the local entropy limiter are employed. At this mesh
spacing, artificial viscosity alone is not sufficient to stabilize the
solution. Good agreement in temperature and pressure with the reference
solution is observed. As shown in Figure 7.7b, which zooms in on the shock,
there are slight discrepancies in the predictions of the leading-shock-front
location and the induction region. However, these predictions improve with
increasing $p$. Also included in Figure 7.7b is a $p=1$ solution with standard
DGODE in order to ensure that entropy-stable DGODE is not the primary cause of
the disagreement. Figure 7.7d shows small-scale entropy oscillations behind
the leading shock front in the $p=3$ calculation. The presence of these
instabilities is not too surprising given the linear nature of the limiting
operator, coupled with the increased susceptibility to oscillations of high
polynomial orders. Nevertheless, these instabilities do not significantly
pollute the rest of the solution.

(a) Temperature. (b) Temperature, zoomed in on shock. (c) Pressure. (d)
Specific thermodynamic entropy.

Figure 7.7:  $p=1$, $p=2$, and $p=3$ results at $t=235$ $\mu\mathrm{s}$ for
the one-dimensional hydrogen detonation test case. The mesh spacing in these
simulations is $h=4.5\times 10^{-4}$ m, which is fives time larger than for
the reference solution. Artificial viscosity and the local entropy limiter are
employed.

Figure 7.8 compares mass-fraction profiles of selected species obtained with
$p=3$ to those of the reference solution. Marginal oscillations can be
observed behind the shock, particularly in the mass-fraction profiles of
$\mathrm{H_{2}O}$ and $\mathrm{HO_{2}}$. In Figure 7.8c, the $p=3$ solution
predicts slightly lower peaks in the mass fractions of $\mathrm{HO_{2}}$ and
$\mathrm{H_{2}O_{2}}$. Furthermore, the aforementioned discrepancy in the
leading-shock location is reflected in these results. Nevertheless, there is
very good agreement between the $p=3$ solution and the reference solution,
illustrating the ability of the developed formulation to obtain stable and
accurate detonation results on coarser meshes.

(a)  (b)  (c)

Figure 7.8:  Comparison of $p=3$ predictions of species mass fractions with
those of the reference solution for the one-dimensional hydrogen detonation
test case. The mesh spacing is $h=4.5\times 10^{-4}$ m, which is fives time
larger than for the reference solution.

Finally, Figure 7.9 gives the percent error in conservation of mass, energy,
and atomic elements for $p=1$ as a representative example, calculated every
$0.235\;\mu\mathrm{s}$ (for a total of 1000 samples). $\mathsf{N}_{O}$,
$\mathsf{N}_{H}$, and $\mathsf{N}_{Ar}$ denote the total numbers of oxygen,
hydrogen, and argon atoms in the mixture. The error remains negligible
throughout the simulation, confirming that the methodology is conservative.

Figure 7.9: Percent error in conservation of mass, energy, and atomic elements
for the $p=1$ calculation with $h=4.5\times 10^{-4}$ m. The initial conditions
for this one-dimensional hydrogen detonation problem are given in Equation
(7.4).

##  8 Concluding remarks