unreacted medium, as seen in solution time $40$-$80$ $\mu\mathrm{s}$.
Eventually the reaction zone behind the front begins to influence the shock
and an overdriven detonation event occurs. The detonation velocity, $v_{det}$,
falls to the steady value of $1618.8$ m/s which was calculated using linear
regression with $r^{2}=0.999994$. The calculated front location and the
approximated linear propagation using the computed detonation velocity are
shown in Figure 5.4. The Chapman-Jouget detonation velocity,
$v_{det}=D_{CJ}=1617.5$ m/s, as well as the post shock profiles for
temperature, pressure, and species concentrations can be calculated using
Shock and Detonation Toolbox [67]. Figures 5.7a and 5.7b show the solution at
time $t=235$ $\mu\mathrm{s}$ in direct comparison to the Shock Detonation
Toolbox computed results (dashed black lines) for temperature, pressure, and
species mass fractions. Additionally, Table 2 summarizes the post shock values
and detonation velocity corresponding to the current simulation and compares
them with the values calculated using the Shock Detonation Toolbox and the
values reported by Houim and Kuo [27] as well as Lv and Ihme [16]. The results
presented in this work agree well with all past reported solutions and the
Shock Detonation Toolbox in particular.

(a) Velocity of the leading shock front as function of time, the front
location is defined as the right most location where the jump in temperature
from $300$ K was larger than $10$ K and a second order finite differencing
scheme was used to calculate the instantaneous velocity. (b) The calculated
front location and the approximated linear propagation using the calculated
velocity of $1618.8$ m/s using linear regression with $r^{2}=0.999994$.

Figure 5.6: Front location and velocity as a function of physical time.

(a) Comparison of the temperature and pressure profiles corresponding to
$\mathrm{DG}\left(p=1\right)$ (solid red lines) and the solution computed via
the Shock Detonation Toolbox (dashed black lines) for the one-dimensional
hydrogen detonation described in Section 5.3. (b) Comparison between Shock
Detonation Toolbox (dashed black lines) and current simulations (solid colored
lines) for species mass fractions.

Figure 5.7: Comparison of the $\mathrm{DG}\left(p=1\right)$ solution for $h=9\times 10^{-5}$ m and the solution computed via the Shock Detonation Toolbox at at time $t=235$ $\mu\mathrm{s}$ for the one-dimensional hydrogen detonation described in Section 5.3. The discontinuous solution at the detonation front is stabilized with the addition of artificial viscosity. Table 2: Table of post-shock values for one-dimensional detonation wave case. Reference |  $p_{vN}$ [kPa] |  $T_{vN}$ [K] |  $v_{det}$ [m/s]
---|---|---|---
Shock and Detonation Toolbox [67] | 173.6 | 1904.7 | 1617.5
Current Simulations | 173.5 | 1904.7 | 1618.8
Houim and Kuo [27] | 173.9 | 1915.2 | 1619.8
Lv and Ihme [16] | 179.4 | 1926.0 | 1634.6

Spurious pressure oscillations were not generated during the simulation of
this test case, indicating the formulation can stably compute detonations. For
this test case, the integrated total energy, $\int_{\Omega}\left(\rho
e_{t}\right)_{0}d\Omega$, and the density,
$\int_{\Omega}\left(\rho\right)_{0}d\Omega$, are conserved as the domain has
walls at both the left and right boundaries. Figure 5.5 shows the percent loss
of the integrated total energy, $\left(\int_{\Omega}\left(\rho
e\right)_{0}d\Omega-\int_{\Omega}\left(\rho
e\right)_{t}d\Omega\right)/\int_{\Omega}\left(\rho e\right)_{0}d\Omega\times
100$, as well as the percent loss of the integrated density,
$\left(\int_{\Omega}\left(\rho\right)d\Omega-\int_{\Omega}\left(\rho\right)_{t}d\Omega\right)/\int_{\Omega}\left(\rho\right)_{0}d\Omega\times
100$, as a function of time. The error remains on the order of
${1}\mathrm{e}{-13}$% for the entire simulation which is expected for a
conservative method.

Figure 5.8: Percent loss of conservation quantities,
$\left(\int_{\Omega}\left(\phi\right)_{0}-\int_{\Omega}\left(\phi\right)_{t}\right)/\int_{\Omega}\left(\phi\right)_{0}\times
100$, for the one-dimensional detonation computed using
$\mathrm{DG}\left(p=1\right)$ for $h=9\times 10^{-5}$ m as a function of
physical time computed using the initial integrated quantity,
$\int_{\Omega}\left(\phi\right)_{0}d\Omega$, and current integrated quantity,
$\int_{\Omega}\left(\phi\right)_{t}d\Omega$, where $\phi=\rho e_{t}$ for total
energy and $\phi=\rho$ for density.

###  5.4 Two-dimensional detonation wave

Here we present the two-dimensional extension of the test case presented in
Section 5.3 with a $0.45$ m long domain that has a channel height of $0.06$ m
,
$\text{$\Omega$}=\left(0,0.45\right)\mathrm{m}\times\left(0,0.06\right)\mathrm{m}$,
with a uniform quadrilateral mesh with grid spacing $h=9\times 10^{-5}$ m.
Each quadrilateral element was split into two triangular elements. The results
presented here used a $\text{DG}\left(p=1\right)$ approximation. The top and
bottom boundaries are simulated as walls. The initial conditions are the same
as Section 5.3 with the exception of two additional high pressure and high
temperature regions, located within the circles
$\sqrt{\left(x-0.019\right)^{2}+\left(y-0.015\right)^{2}}=0.0025$ m and
$\sqrt{\left(x-0.010\right)^{2}+\left(y-0.044\right)^{2}}=0.0025$ m, with
conditions

$\begin{array}[]{ccc}\left(v_{1},v_{2}\right)&=&\left(0,0\right)\text{
m/s},\\\ X_{Ar}:X_{H_{2}O}:X_{OH}&=&8:2:0.01,\\\ p&=&{5.5}\mathrm{e}{5}\text{
Pa},\\\ T&=&3500\text{ K}.\end{array}$ (5.4)

These circular regions were added to perturb the detonation and develop
cellular structures and they where aligned with the unstructured grid
interfaces at initialization as shown in Figure (5.9). This test case is
similar to one performed by Oran et al. [30] with the exception that they
initialized their solution from a one-dimensional detonation and their domain
was $0.10$ m longer,
$\text{$\Omega$}=\left(0,0.55\right)\mathrm{m}\times\left(0,0.06\right)\mathrm{m}$.
Their simulations revealed two detonation cells in the vertical direction with
each cell computed to be
$\left(0.055\right)\mathrm{m}\times\left(0.03\right)\mathrm{m}$. Their results
were consistent with experiments performed by Lefebvre et al. [75]. Houim and
Kuo [27], as well as Lv and Ihme [16], ran similar test cases but the domain
was half the height of the one considered in this work. Houim and Kuo [27]
shifted their simulation to be in the detonation frame of reference whereas Lv
and Ihme [16] ran the simulation in the laboratory frame of reference. In both
cases, the solutions were initialized from the corresponding one-dimensional
detonation solutions. Their simulations revealed one detonation cell, which is
is consistent with domain height considered in their respective works.

Figure (5.9) presents the temperature solutions as the detonation front
progresses through the simulation domain. As the initial shock progresses into
the domain, it collides with the two circular high pressure and high
temperature regions before $t=30$ $\mu\mathrm{s}$, the detonation is then
established and the perturbations lead to transverse waves traveling in the
vertical directions that reflect off the top and bottom walls. To the left of
the detonation, complex, two-dimensional structures develop, including
Kelvin–Helmholtz instabilities as well as the triple point that connects the
Mach stems and incident shock.

Figure 5.9: Temperature contours for two-dimensional detonation wave. The
discontinuous solution at the detonation front is stabilized with the addition
of artificial viscosity. The corresponding one-dimensional initialization
given by (5.2) is augmented with two additional circular two additional high
pressure and high temperature regions (5.3).

Figure (5.10) shows the adapted polynomial degree for the DGODE chemistry
solve at each degree of freedom as the detonation propagates through the
simulation domain. The adapted polynomial degree is largest in the vicinity of
the detonation, reaching $p=3$ in the regions behind the detonation front and
at the triple points. The majority of the adapted polynomial degree reduces to
$p=0$ downstream of the shock. Thin regions of $p=1$ exist downstream of the
shock along the vertically traveling transverse waves. This indicates that the
traveling waves are strong enough to push the chemistry out of equilibrium and
require high-order integration to solve the chemistry. These thin regions of
higher order are similar to the spike in the adapted polynomial degree as
shown in Section 5.3’s one-dimensional detonation results at time $t=94$
$\mu\mathrm{s}$, see Figure 5.5.

Figure 5.10: The adapted polynomial degree of the DGODE approximation for two-
dimensional detonation wave.

Figure 5.11 shows the maximum pressure experienced at each location as the
solution evolved, $max(p,p_{max})$. This metric is a good indicator of the
detonation cell structures and shows two distinct detonation cells in the
vertical direction. A large burst of pressure exists between $x=0.04$ m and
$x=0.1$ m from the over driven detonation as the shock transitions to a
detonation in the early stages of development. The last solution, at $t=240$
$\mu\mathrm{s}$, displays two dashed white lines that indicated the measured
horizontal and vertical size of a detonation cell. For this solution the
detonation cell was computed to be
$\left(0.055\right)\mathrm{m}\times\left(0.03\right)\mathrm{m}$ which is
consistent with experimental work [75] and previous numerical results [30].

Figure 5.11: Maximum pressure, $\max(p,p_{max})$, revealing two distinct
detonation cellular structure, cf. [30] and [75].

No spurious pressure oscillations were noticed for the entire simulation of
this test case, indicating that the formulation presented in this work is
stable for two-dimensional detonations. For this test case, the integrated
total energy, $\int_{\Omega}\left(\rho e_{t}\right)_{0}d\Omega$, and the
density, $\int_{\Omega}\left(\rho\right)_{0}d\Omega$, are conserved as the
domain has walls at all boundaries. Figure 5.12 shows the percent loss of the
integrated total energy, $\left(\int_{\Omega}\left(\rho
e\right)_{0}d\Omega-\int_{\Omega}\left(\rho
e\right)_{t}d\Omega\right)/\int_{\Omega}\left(\rho e\right)_{0}d\Omega\times
100$, as well as the percent loss of the integrated density,
$\left(\int_{\Omega}\left(\rho\right)d\Omega-\int_{\Omega}\left(\rho\right)_{t}d\Omega\right)/\int_{\Omega}\left(\rho\right)_{0}d\Omega\times
100$, as a function of time. The error remains on the order of $10^{-13}$ %
for the entire simulation which is on the order of the expected machine error
and assures that the method remains conservative.

Figure 5.12: Percent loss of conservation quantities,
$\left(\int_{\Omega}\left(\phi\right)_{0}-\int_{\Omega}\left(\phi\right)_{t}\right)/\int_{\Omega}\left(\phi\right)_{0}\times
100$, for the two-dimensional detonation computed computed using the initial
integrated quantity, $\int_{\Omega}\left(\phi\right)_{0}d\Omega$, and current
integrated quantity, $\int_{\Omega}\left(\phi\right)_{t}d\Omega$, where
$\phi=\rho e_{t}$ for total energy and $\phi=\rho$ for density.

###  5.5 One-dimensional premixed flame

In this section we consider a fully compressible multi-component reacting
Navier-Stokes flow corresponding to a one-dimensional premixed flame using the
Hydrogen-Air chemistry from [76]. The chosen domain is $6$ cm in length and
discretized with uniform linear line elements with $h=5\times 10^{-5}\,$m. The
initialization is given by

$\begin{array}[]{cccc}\left(p,T,X_{O_{2}},X_{H_{2}},X_{N_{2}}\right)&=&\left(1\,\mathrm{\text{bar}},300\,\mathrm{K},0.170,0.188,0.642\right)&x<5.8\,\mathrm{\text{cm}}\\\
y&=&\hat{y}&x>5.8\,\mathrm{\text{cm}}\end{array},$ (5.5)

where $\hat{y}$ corresponds to the fully reacted state from Cantera’s
homogeneous constant pressure reactor simulation calculated from the unreacted
conditions specified at $x<5.8$ cm [66]. The solution at the initially
discontinuous interface immediately diffuses to form a continuous, smooth,
solution. No additional stabilization, e.g., artificial viscosity, limiting or
filtering was needed for the duration of the simulation.

The right hand side boundary is a fixed pressure outflow characteristic
boundary condition with $p_{\infty}=1$ bar. The left hand side boundary is a
characteristic farfield wave boundary condition with $T_{\infty}=300$ K and
$v_{\infty}=0$ m/s. This boundary allows any pressure waves caused by the
initialization to exit the domain. The initialization contains a temperature
and species discontinuity which gives rise to a pressure wave that leaves the
system as the reaction front diffuses into the unreacted region. Eventually,
the solution converges to a propagating flame traveling at a unique velocity.

Table 3: Computed flame speeds corresponding to the one-dimensional premixed flame described in Section 5.5.  Method | Flame Speed [m/s]
---|---
Cantera | 0.643
$\mathrm{DG}\left(p=1\right)$ | 0.641
$\mathrm{DG}\left(p=2\right)$ | 0.643

Figures 5.13 and 5.14 present the $\mathrm{DG}\left(p=1\right)$ and
$\mathrm{DG}\left(p=2\right)$ solutions, respectively. The temperature and
species mass fraction profiles are compared to the Cantera flame solution with
$h=10^{-5}\,$ m and the profiles are shifted so that $T=400$ K at $x=0$ m. The
$\mathrm{DG}\left(p=1\right)$ solutions reach the correct reacted state but
cannot fully resolve the flame structures in the $-0.0005<x<0$ region. The
$\mathrm{DG}\left(p=2\right)$ solution overcomes these errors and is in good
agreement with the Cantera solution. Despite the under resolved profiles in
the $\mathrm{DG}\left(p=1\right)$ solution, the flame speeds corresponding to
both solutions compare well to the flame speed calculated in Cantera as
presented in Table 3. The Cantera flame speed, given as the inflow velocity
for the constant mass flow-rate, is 0.643 m/s. We considered the flame front
in the unsteady $\mathrm{DG}\left(p=1\right)$ and
$\mathrm{DG}\left(p=2\right)$ solutions to be the location corresponding to
$T=1000$ K. We tracked this location and computed a steady velocity of 0.641
m/s and 0.643 m/s for the $\mathrm{DG}\left(p=1\right)$ and
$\mathrm{DG}\left(p=2\right)$ solutions, respectively.

Figures 5.13d and 5.14d show the pressure through the entire computational
domain for $\mathrm{DG}\left(p=1\right)$ and $\mathrm{DG}\left(p=2\right)$,
respectively. For the $\mathrm{DG}\left(p=1\right)$ solution, there are small
oscillations, on the order of $0.25$% of the ambient pressure. These
oscillations are not present for the $\mathrm{DG}\left(p=2\right)$ solution,
where only a slight variation is seen through the flame front but is constant
on both sides of the flame within $0.1$% of the desired ambient pressure. The
lack of pressure oscillations in the higher order solution indicate that the
$\mathrm{DG}\left(p=1\right)$ solution is under-resolved.

(a) The temperature profile correspond to the $\mathrm{DG}\left(p=1\right)$
solution (grey lines) and the Cantera solution (solid lines) on a uniform
$h=5\times 10^{-5}\,$m grid. Both the $\mathrm{DG}\left(p=1\right)$ and
Cantera solutions have been shifted so that $T=400$ K at $x=0$. (b) Minor
species profiles correspond to the $\mathrm{DG}\left(p=1\right)$ solution
(grey lines) and the Cantera solution (solid lines) on a uniform $h=5\times
10^{-5}\,$m grid. Both the $\mathrm{DG}\left(p=1\right)$ and Cantera solutions
have been shifted so that $T=400$ K at $x=0$.  (c) Major species profiles
correspond to the $\mathrm{DG}\left(p=1\right)$ solution (grey lines) and the
Cantera solution (solid lines) on a uniform $h=5\times 10^{-5}\,$m grid. Both
the $\mathrm{DG}\left(p=1\right)$ and Cantera solutions have been shifted so
that $T=400$ K at $x=0$. (d) The pressure profile corresponding to the
$\mathrm{DG}\left(p=1\right)$ solution (grey line) for the entire simulation
domain. For reference the the pressure profile corresponding to a constant 1
bar solution is included (solid red line). The initialization is given by
(5.5).

Figure 5.13: The profiles corresponding to the $\mathrm{DG}\left(p=1\right)$
solution for the one-dimensional premixed flame, described in Section 5.5,
computed on a uniform grid with on a uniform $h=5\times 10^{-5}\,$m grid. The
initialization is given by (5.3). The flame speed of 0.641 m/s corresponding
to the $\mathrm{DG}\left(p=1\right)$ solution compares well with the flame
speed corresponding to the Cantera solution of 0.643 m/s. This smooth solution
of a multi-component chemically reacting Navier-Stokes flow was computed
without the use of additional stabilization, e.g., artificial viscosity,
limiting or filtering. The pressure oscillations visible in Figure 5.13d are
not present in the corresponding $\mathrm{DG}\left(p=2\right)$ solution shown
in Figure 5.14d, indicating the $\mathrm{DG}\left(p=1\right)$ solution is
under-resolved.

(a) The temperature profile correspond to the $\mathrm{DG}\left(p=2\right)$
solution (grey lines) and the Cantera solution (solid lines) on a uniform
$h=5\times 10^{-5}\,$m grid. Both the $\mathrm{DG}\left(p=2\right)$ and
Cantera solutions have been shifted so that $T=400$ K at $x=0$. (b) Minor
species profiles correspond to the $\mathrm{DG}\left(p=2\right)$ solution
(grey lines) and the Cantera solution (solid lines) on a uniform $h=5\times
10^{-5}\,$m grid. Both the $\mathrm{DG}\left(p=2\right)$ and Cantera solutions
have been shifted so that $T=400$ K at $x=0$. (c) Major species profiles
correspond to the $\mathrm{DG}\left(p=2\right)$ solution (grey lines) and the
Cantera solution (solid lines) on a uniform $h=5\times 10^{-5}\,$m grid. Both
the $\mathrm{DG}\left(p=2\right)$ and Cantera solutions have been shifted so
that $T=400$ K at $x=0$. (d) The pressure profile corresponding to the
$\mathrm{DG}\left(p=2\right)$ solution (grey line) for the entire simulation
domain. For reference the the pressure profile corresponding to a constant 1
bar solution is included (solid red line).

Figure 5.14: The profiles corresponding to the $\mathrm{DG}\left(p=2\right)$
solution for the one-dimensional premixed flame, described in Section 5.5,
computed on a uniform grid with on a uniform $h=5\times 10^{-5}\,$m grid. The
initialization is given by (5.3). The flame speed of 0.643 m/s corresponding
to the $\mathrm{DG}\left(p=2\right)$ solution compares well with the flame
speed corresponding to the Cantera solution of 0.643 m/s. This smooth solution
of a multi-component chemically reacting Navier-Stokes flow was computed
without the use of additional stabilization, e.g., artificial viscosity,
limiting or filtering.

###  5.6 Hydrogen-Air chemically reacting shear layer

We use the viscous reacting formulation of Equations (3.6)- (3.8) with
detailed chemical kinetics and transport to approximate the solution of a
chemically reacting shear layer using the Hydrogen-Air chemistry extracted
from the GRI Mech II mechanism [76]. This test case was designed to test the
stability and robustness of the formulation presented in this work and verify
that smooth multidimensional multi-component reacting flows can be computed
without artificial viscosity, limiting, or filtering. The entire simulation
domain is depicted in Figure 5.15. The bounding domain is $0.10$ m in length
and $0.05$ m in height,
$\text{$\Omega$}=\left(0,0.10\right)\mathrm{m}\times\left(0,0.05\right)\mathrm{m}$.
A splitter plate is used to separate the incoming streams, anchor the flame,
and to test slip and adiabatic wall boundary conditions described in Section
3.1. The plate is $0.001$ m thick and $0.005$ m long and is centered at
$y=0.025$ m. The flow is initialized with the following conditions

$\begin{array}[]{ccc}\left(M,T,p,Y_{O_{2}},Y_{N_{2}},Y_{H_{2}}\right)&=&\begin{cases}\left(0\text{.8},\>2000\text{
K},\>101325\text{ Pa},\>0.21,\>0.79,\>0\right)&y<0.025\text{ m}\\\
\left(0\text{.8},\>2000\text{ K},\>101325\text{
Pa},\>\phantom{.21}0,\>\phantom{.79}0,\>1\right)&y>0.025\text{
m}\end{cases},\end{array}$ (5.6)

where $M$ is the freestream Mach number.

The upstream half of the plate, $0\,\mathrm{m}<x<0.0025\,\mathrm{m}$, is
specified as a slip wall, whereas the downstream half of the plate,
$0.0025\,\mathrm{m}<x<0.005\,\mathrm{m}$, is an adiabatic no-slip wall. The
remaining boundaries, at $x=0\,\mathrm{m}$, $x=0.1\,\mathrm{m}$,
$y=0\,\mathrm{m}$ , and $y=0.05\,\mathrm{m}$, are specified as non-reflecting
characteristic boundary conditions where the freestream state is given by
(5.6). The characteristic conditions downstream of the plate at
$x=0.1\,\mathrm{m}$ allow for pockets of supersonic flow to exit the domain if
the heat release from reactions is large enough to transition the flow to
supersonic.

Initial simulations were run to determine the regions of refinement required
to adequately capture the reacting mixing layer. Figure 5.15 presents the
final grid corresponding to the entire domain, which was constructed using
Gmsh [77]. The finest resolution in the mixing region corresponds to a mesh
size of $200\times 10^{-6}\,\mathrm{m}$ and transitions to $500\times
10^{-6}\,\mathrm{m}$ at $x=0.1$ m. This region begins at a thickness of
$0.0035\,\mathrm{m}$, $0.02325\,\mathrm{m}<y<0.02675\,\mathrm{m}$, at
$x=0\,\mathrm{m}$, and grows to a thickness of $0.015\,\mathrm{m}$,
$0.0155\,\mathrm{m}<y<0.0355\,\mathrm{m}$, at $x=0.1\,\mathrm{m}$.

Figure 5.15: Splitter plate mesh comprised of linear triangular elements with
increased resolution mixing region. This mesh was generated using Gmsh [77].

Figures 5.16a and 5.16b show the temporal evolution for $\mathrm{DG}(p=1)$ and
$\mathrm{DG}(p=2)$ temperature solutions, respectively. The initial interface
between H2 and air diffuses and reactions occur along the mid-line. As the
reactions progress, flow features develop that create mixing structures that
grow in the downstream direction. As the $\mathrm{DG}(p=1)$ and
$\mathrm{DG}(p=2)$ solutions evolve in time, they begin to diverge from each
other. Figure 5.17 shows a direct comparison of the solutions at $t=51$
$\mu$s. The profiles for the $\mathrm{DG}(p=2)$ solution are thinner and less
diffuse than the $\mathrm{DG}(p=1)$ solution which contributes to the
difference in the downstream features. Figures 5.18a and 5.18b show the
upstream and down stream $Y_{H_{2}O}$ profiles, respectively, for both the
$\mathrm{DG}(p=1)$ and $\mathrm{DG}(p=2)$ at $t=51$ $\mu$s. The
$\mathrm{DG}(p=2)$ solution resolves a sharper flame attachment to the
splitter wall than the $\mathrm{DG}(p=1)$ solution. In the downstream region a
large structure is apparent in both solutions, with the $\mathrm{DG}(p=2)$
solution resolving more fine scale features than the $\mathrm{DG}(p=1)$
solution.

(a)  Evolution of temperature profile using $\mathrm{DG}\left(p=1\right)$ for
H2-air splitter plate.  (b) Evolution of temperature profile using
$\mathrm{DG}\left(p=2\right)$ for H2-air splitter plate.

Figure 5.16:  Evolution of temperature profile corresponding to the
$\mathrm{DG}\left(p=1\right)$ and $\mathrm{DG}\left(p=2\right)$ solutions at
different times. The initialization is given by (5.6). This smooth solution of
a multi-component chemically reacting Navier-Stokes flow was computed without
the use of additional stabilization, e.g., artificial viscosity, limiting or
filtering. Figure 5.17: Comparison of splitter plate temperature solution for
$\mathrm{DG}\left(p=1\right)$ and $\mathrm{DG}\left(p=2\right)$ at $t=51$
$\mu$s.

(a)  Upstream structures of $H_{2}O$ for $\mathrm{DG}\left(p=1\right)$ and
$\mathrm{DG}\left(p=2\right)$ at $t=51$ $\mu$s. (b) Downstream structures of
$H_{2}O$ for $\mathrm{DG}\left(p=1\right)$ and $\mathrm{DG}\left(p=2\right)$
at $t=51$ $\mu$s..

Figure 5.18:  Comparison of splitter plate $H_{2}O$ mass fraction solution for
$\mathrm{DG}\left(p=1\right)$ and $\mathrm{DG}\left(p=2\right)$ at $t=51$
$\mu$s. The initialization is given by (5.6). This smooth solution of a multi-
component chemically reacting Navier-Stokes flow was computed without the use
of additional stabilization, e.g., artificial viscosity, limiting or
filtering.

We extended this two-dimensional shear layer problem to three dimensions by
extruding the mesh in the $z$-direction. To save on computational cost, the
domain was changed to $0.03$ m in height but maintained a length of $0.10$ m,
$\text{$\Omega$}=\left(0,0.10\right)\mathrm{m}\times\left(0,0.03\right)\mathrm{m}$.
The splitter plate is the same size but is centered at $y=0.015$ m. The mixing
region resolution and size is the same as the two dimensional problem but
shifted to center the new plate location. The mesh was constructed using Gmsh
and was extruded in the $z$-direction at a cellular width of $200\times
10^{-6}$ m for a total width of $0.0028$ m. The boundary conditions in the
resulting $xz$-planes and $yz$-planes were the same as the two dimensional
splitter plate shifted to be centered at $y=0.015$ m for the left hand side
boundary with the exception of a spatiotemporal variation introduced to the
momentum at the $x=0$ boundaries,

$\begin{array}[]{cccc}M&=&0.8,\\\
v_{1}&=&Mc+\frac{M}{10}c\sin\left({\frac{\pi}{w}z}\right)\sin\left({\frac{2\pi}{\tau}t}\right)&\text{
m/s, }\\\ v_{2}&=&0&\text{ m/s, }\\\
v_{3}&=&\frac{M}{5}c+\frac{M}{50}c\sin\left({\frac{\pi}{w}z}\right)\sin\left({\frac{2\pi}{\tau}t}\right)&\text{
m/s, }\end{array}$ (5.7)

where $c=\sqrt{\gamma RT}$ with $\gamma$ and $R$ calculated from the mixture
properties at the boundary. This inflow boundary condition was used to
encourage three-dimensional breakup and unsteadiness in the $z$-direction.
Periodicity is enforced between the resulting $xy$-planes at $z=0$ m and
$z=0.0028$ m. The initialization is the same as in two dimensional case, which
is given by (5.6) at time $t=0\,\mathrm{s}$. The $\mathrm{DG}(p=1)$
approximation was used for the first $100$ $\mu$s of the simulation, at which
point the approximate solution was interpolated to the space corresponding to
the $\mathrm{DG}(p=2)$ approximation. No artificial viscosity, limiting, or
filtering was required throughout the entire simulation as the physical
diffusion present in the solution provided adequate stabilization.

Figure 5.19 shows the $Y_{H_{2}O}$ solution for $xy$ -planes located at $z=0$
m and $z=0.0014$ m. The variation of $Y_{H_{2}O}$ in the z-direction is
greater downstream of $x=0.04$ m. This is seen particularly in the solution
for $y>0$, where the flow-field inertia is not as great since the species has
less mass and is therefore more susceptible to mixing due to the inflow
perturbations. Figure 5.20 shows the three dimensional isocontours for
$T=2200$ K between times $t=350$ and $t=370$ $\mu$s. The splitter plate is
shown as a black surfaces, where as the isocontours are colored by the
absolute value of the $z$-component of the isocontour surface normal, labeled
as $|n_{z}|$. Red and white surfaces represent areas where the flame is
wrinkled in the $z$-direction, whereas blue surfaces represent areas of the
flame that are planar in the $z$-direction, $|n_{z}|\rightarrow 0$. These
isocontours confirm that the flame becomes less uniform in the $z$-direction
in the region downstream of $x=0.04$ m, particularly in the top half of the
domain, $y>0.015$ m.

Figure 5.19: Three dimensional $\mathrm{DG}\left(p=2\right)$ splitter plate
$Y_{H_{2}O}$ xy-plane solutions at two different z-locations at $t=350$
$\mu$s. This smooth solution of a multi-component chemically reacting Navier-
Stokes flow was computed without the use of additional stabilization, e.g.,
artificial viscosity, limiting or filtering. Figure 5.20: Three dimensional
$\mathrm{DG}\left(p=2\right)$ splitter plate $T=2200$ K isocontour solution
from $t=350$ to $t=370$ $\mu$s. This smooth solution of a multi-component
chemically reacting Navier-Stokes flow was computed without the use of
additional stabilization, e.g., artificial viscosity, limiting or filtering.

These results indicate that the conservative DG formulation presented in this
work is capable of computing smooth multidimensional, multi-component,
chemically reacting Navier-Stokes flows without the need for artificial
viscosity, limiting, filtering. Furthermore, since the formulation maintains
discrete temperature and pressure equilibrium it is capable of incorporating
the viscous and diffusive operators without caveats or creating additional
concerns. As expected, the formulation behaves exactly like its single-
component counterpart for the case of viscous multi-component reacting flows.

##  6 Conclusion

We have presented a a fully conservative discontinuous Galerkin finite element
method for the multi-component chemically reacting Navier-Stokes equations
that retains the desirable properties of DG, namely discrete conservation and
high-order accuracy in smooth regions of the flow. In contrast to previous DG
implementations for multi-component chemically reacting flows, the formulation
presented in this work does not freeze the thermodynamics; instead a nonlinear
relationship for temperature is solved which ensures consistency between the
internal energy of the conserved state and the internal energy defined by a
mixture averaged polynomial expression based on temperature. Furthermore, the
discrete solution is represented in terms of a nodal basis with coefficients
defined on the element interfaces so that temperature and pressure equilibrium
between adjacent elements is maintained and the resulting weak form is
evaluated in terms of basis coefficients that are designed to maintain
pressure equilibrium. As such, unphysical pressure oscillations are not
generated at material interfaces if the temperature is continuous.
Furthermore, if the temperature is discontinuous, the unphysical pressure
oscillations remain small and do not lead to numerical failure. Other
numerical instabilities do arise at discontinuous interfaces, e.g., shocks and
detonations, but these are not specific to multi-component flows and they are
suppressed in this work by the addition of a residual based artificial
viscosity.

A coupled $hp$-adaptive ODE solver, DGODE, has been developed for the
integration of the stiff chemical source term. The solver adapts the local
time-step, $h$, as well as the polynomial degree, $p$, corresponding to the
local approximation in order to efficiently integrate the system of ODEs. The
solver was shown to successfully integrate the stiff source term corresponding
to the one- and two-dimensional $H_{2}$-$O_{2}$-$Ar$ detonation waves, where
the local polynomial degree adapted to the local stiffness of the chemical
system.

The generation of unphysical pressure oscillations was studied in the context
of the formulation described in this work. We found that pressure oscillations
generated for formulations based on frozen thermodynamics did not reach the
magnitudes previously reported in literature [25, 26, 15, 27, 16]. However,
when it was assumed that the mean specific heat at constant pressure was the
same as the NASA polynomial specific heat at constant pressure,
$\bar{c_{p}}=c_{p}$, the oscillations did reach the same levels as previously
reported, indicating that the severity of the pressure oscillations is
directly related to the method in which the thermodynamic quantities are
evaluated. Regardless, frozen thermodynamic formulations generate pressure
oscillations that grow in time as the solution evolves. In contrast,
evaluating the thermodynamics quantities exactly does not lead to pressure
oscillations at a material interface if the temperature is continuous. We
therefore concluded that the conservative DG discretization presented in this
work can stably compute solutions to smooth multi-component chemically
reacting Navier-Stokes flows without additional stabilization.

The discretization was applied to several multi-component non-reacting and
reacting test cases including: a one-dimensional $He$ and $N_{2}$ shock-tube,
a two dimension shock interaction with a helium bubble suspended in air, one-
dimensional $H_{2}$-$O_{2}$-$Ar$ detonation wave, and two-dimensional
$H_{2}$-$O_{2}$-$Ar$ detonation. We confirmed that additional stabilization is
not required for well-resolved smooth regions of the flow. Where applicable,
i.e., for wall-bounded problems, we computed the integrated percent loss of
the conserved quantities at each time-step. The value was shown to be on the
order of machine error while remaining constant for the duration of the
simulation. The solutions for the one- and two-dimensional
$H_{2}$-$O_{2}$-$Ar$ detonation simulations were compared to results
previously reported in the literature. The profile of the one-dimensional
solution in the region of the detonation was shown to be in close agreement
with the profile predicted by the detonation toolbox and two-dimensional
solution was shown to reproduce the correct cellular detonation structure.

We also presented a one-dimensional premixed hydrogen-air flame and two- and
three-dimensional reacting shear layer with a splitter plate. The $H_{2}$-air
deflagration flame speed compared well to the speed computed via Cantera for
both the $\mathrm{DG}(p=1)$ and $\mathrm{DG}(p=2)$ solutions. Discrepancies
between the Cantera solution and the $\mathrm{DG}(p=1)$ solution were not
observed in the $\mathrm{DG}(p=2)$ solution, indicating the $\mathrm{DG}(p=2)$
solution better resolved the flame. The two- and three-dimensional shear layer
simulations confirmed, in a practical setting, that the conservative DG
formulation does not require additional stabilization for smooth multi-
component chemically reacting Navier-Stokes flows.

## Acknowledgements

This work is sponsored by the Office of Naval Research through the Naval
Research Laboratory 6.1 Computational Physics Task Area.

##

## References

  * [1] R. Johnson, A. Kercher, A conservative discontinuous Galerkin discretization for the total energy formulation of the reacting Navier Stokes equations, arXiv preprint arXiv:1910.10544.
  * [2] W. H. Reed, T. Hill, Triangular mesh methods for the neutron transport equation, Tech. rep., Los Alamos Scientific Lab., N. Mex.(USA) (1973).
  * [3] F. Bassi, S. Rebay, A high-order accurate discontinuous finite element method for the numerical solution of the compressible Navier–Stokes equations, Journal of Computational Physics 131 (2) (1997) 267–279.
  * [4] F. Bassi, S. Rebay, High-order accurate discontinuous finite element solution of the 2D Euler equations, Journal of Computational Physics 138 (2) (1997) 251–285.
  * [5] B. Cockburn, C.-W. Shu, The Runge–Kutta discontinuous Galerkin method for conservation laws V: multidimensional systems, Journal of Computational Physics 141 (2) (1998) 199–224.
  * [6] B. Cockburn, G. Karniadakis, C.-W. Shu, The development of discontinuous Galerkin methods, in: Discontinuous Galerkin Methods, Springer, 2000, pp. 3–50.
  * [7] D. Arnold, F. Brezzi, B. Cockburn, L. Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis 39 (5) (2002) 1749–1779.
  * [8] R. Hartmann, P. Houston, Adaptive discontinuous Galerkin finite element methods for the compressible Euler equations, Journal of Computational Physics 183 (2) (2002) 508–532.
  * [9] K. Fidkowski, T. Oliver, J. Lu, D. Darmofal, p-Multigrid solution of high-order discontinuous Galerkin discretizations of the compressible Navier–Stokes equations, Journal of Computational Physics 207 (1) (2005) 92–113.
  * [10] H. Luo, J. Baum, R. Löhner, A p-multigrid discontinuous Galerkin method for the euler equations on unstructured grids, Journal of Computational Physics 211 (2) (2006) 767–783.
  * [11] H. Luo, J. Baum, R. Löhner, A hermite weno-based limiter for discontinuous Galerkin method on unstructured grids, Journal of Computational Physics 225 (1) (2007) 686–713.
  * [12] H. Luo, J. Baum, R. Löhner, A discontinuous Galerkin method based on a taylor basis for the compressible flows on arbitrary grids, Journal of Computational Physics 227 (20) (2008) 8875–8893.
  * [13] P.-O. Persson, J. Peraire, Newton-GMRES preconditioning for discontinuous Galerkin discretizations of the Navier–Stokes equations, SIAM Journal on Scientific Computing 30 (6) (2008) 2709–2733.
  * [14] R. Hartmann, T. Leicht, Higher order and adaptive DG methods for compressible flows, in: H. Deconinck (Ed.), VKI LS 2014-03: 37th Advanced VKI CFD Lecture Series: Recent developments in higher order methods and industrial application in aeronautics, Dec. 9-12, 2013, Von Karman Institute for Fluid Dynamics, Rhode Saint Genèse, Belgium, 2014.
  * [15] G. Billet, J. Ryan, A Runge–Kutta discontinuous Galerkin approach to solve reactive flows: The hyperbolic operator, Journal of Computational Physics 230 (4) (2011) 1064 – 1083.  doi:https://doi.org/10.1016/j.jcp.2010.10.025.
  * [16] Y. Lv, M. Ihme, Discontinuous Galerkin method for multicomponent chemically reacting flows and combustion, Journal of Computational Physics 270 (2014) 105 – 137.  doi:https://doi.org/10.1016/j.jcp.2014.03.029.
  * [17] G. P. Smith, D. M. Golden, M. Frenklach, N. W. Moriarty, B. Eiteneer, M. Goldenberg, C. T. Bowman, Gri-mech 3.0, version 3.0 (2000).
  * [18] R. Abgrall, Generalisation of the Roe scheme for the computation of mixture of perfect gases, La Recherche Aérospatiale 6 (1988) 31–43.
  * [19] S. Karni, Multicomponent flow calculations by a consistent primitive algorithm, Journal of Computational Physics 112 (1) (1994) 31 – 43.  doi:https://doi.org/10.1006/jcph.1994.1080.
  * [20] R. Abgrall, How to prevent pressure oscillations in multicomponent flow calculations: A quasi conservative approach, Journal of Computational Physics 125 (1) (1996) 150 – 160.  doi:https://doi.org/10.1006/jcph.1996.0085.
  * [21] P. Jenny, B. Mueller, H. Thomann, Correction of conservative euler solvers for gas mixtures, Journal of Computational Physics 132 (1) (1997) 91 – 107.  doi:https://doi.org/10.1006/jcph.1996.5625.
  * [22] M. Zahr, P.-O. Persson, An optimization-based approach for high-order accurate discretization of conservation laws with discontinuous solutions, Journal of Computational Physics.
  * [23] A. Corrigan, A. Kercher, D. Kessler, A moving discontinuous Galerkin finite element method for flows with interfaces, International Journal for Numerical Methods in Fluids 89 (9) (2019) 362–406.  doi:10.1002/fld.4697.
  * [24] Z. Sun, M. Braack, J. Lang, An adaptive moving finite element method for steady low mach number compressible combustion problems, International Journal for Numerical Methods in Fluids n/a (n/a).  doi:10.1002/fld.4818.
  * [25] R. Abgrall, S. Karni, Computations of compressible multifluids, Journal of Computational Physics 169 (2) (2001) 594 – 623.  doi:https://doi.org/10.1006/jcph.2000.6685.
  * [26] G. Billet, R. Abgrall, An adaptive shock-capturing algorithm for solving unsteady reactive flows, Computers and Fluids 32 (10) (2003) 1473 – 1495.  doi:https://doi.org/10.1016/S0045-7930(03)00004-5.
  * [27] R. Houim, K. Kuo, A low-dissipation and time-accurate method for compressible multi-component flow with variable specific heat ratios, Journal of Computational Physics 230 (23) (2011) 8527 – 8553.  doi:https://doi.org/10.1016/j.jcp.2011.07.031.
  * [28] E. Johnsen, F. Ham, Preventing numerical errors generated by interface-capturing schemes in compressible multi-material flows, Journal of Computational Physics 231 (17) (2012) 5705 – 5717.  doi:https://doi.org/10.1016/j.jcp.2012.04.048.
  * [29] B. Thornber, M. Groom, D. Youngs, A five-equation model for the simulation of miscible and viscous compressible fluids, Journal of Computational Physics 372 (2018) 256 – 280.  doi:https://doi.org/10.1016/j.jcp.2018.06.028.
  * [30] E. S. Oran, E. I. Weber J. W.and Stefaniw, M. H. Lefebvre, J. D. Anderson, A numerical study of a two-dimensional h2-o2-ar detonation using a detailed chemical reaction model, Combustion and Flame 113 (1) (1998) 147 – 163.  doi:https://doi.org/10.1016/S0010-2180(97)00218-6.
  * [31] E. S. Oran, J. P. Boris, Numerical Simulation of Reactive Flow, 2nd Edition, Cambridge University Press, 2000.  doi:10.1017/CBO9780511574474.
  * [32] J. B. Chen, H. G. Im, Stretch effects on the burning velocity of turbulent premixed hydrogen/air flames, Proceedings of the Combustion Institute 28 (1) (2000) 211 – 218.  doi:https://doi.org/10.1016/S0082-0784(00)80213-1.
  * [33] E. R. Hawkes, R. Sankaran, S. J. C., J. H. Chen, Direct numerical simulation of turbulent combustion: fundamental insights towards predictive models, Journal of Physics: Conference Series 16 (2005) 65–79.  doi:10.1088/1742-6596/16/1/009.
  * [34] B. D. Taylor, D. A. Kessler, V. N. Gamezo, E. S. Oran, Numerical simulations of hydrogen detonations with detailed chemical kinetics, Proceedings of the Combustion Institute 34 (2) (2013) 2009 – 2016.  doi:https://doi.org/10.1016/j.proci.2012.05.045.
  * [35] B. Taylor, D. Schwer, A. Corrigan, Implementation of thermochemistry and chemical kinetics in a GPU-based CFD code, in: 53rd AIAA Aerospace Sciences Meeting Including The New Horizons Forum and Aerospace Exposition, 2015, 2015-0842.
  * [36] C. A. Kennedy, M. H. Carpenter, Several new numerical methods for compressible shear-layer simulations, Applied Numerical Mathematics 14 (4) (1994) 397 – 433\.  doi:https://doi.org/10.1016/0168-9274(94)00004-2.
  * [37] R. Deiterding, Parallel adaptive simulation of multi-dimensional detonation structures, Dissertation. de, 2003.
  * [38] A. Y. Poludnenko, B. D. Taylor, Turbulent chemical and thermonuclear flames: Intrinsic instability and anisotropic turbulence amplification, in: TSFP DIGITAL LIBRARY ONLINE, Begel House Inc., 2015.
  * [39] P. Colella, P. Woodward, The piecewise parabolic method (PPM) for gas-dynamical simulations, Journal of Computational Physics 54 (1) (1984) 174–201.
  * [40] B. J. McBride, M. J. Zehe, S. Gordon, NASA Glenn coefficients for calculating thermodynamic properties of individual species.
  * [41] R. J. Kee, J. A. Miller, G. H. Evans, G. Dixon-Lewis, A computational model of the structure and extinction of strained, opposed flow, premixed methane-air flames, Symposium (International) on Combustion 22 (1) (1989) 1479 – 1494.  doi:https://doi.org/10.1016/S0082-0784(89)80158-4.
  * [42] C. R. Wilke, A viscosity equation for gas mixtures, J. Chem. Phys 18 (1950) 517–519.  doi:10.1063/1.1747673.
  * [43] S. Mathur, P. K. Tondon, S. C. Saxena, Thermal conductivity of binary, ternary and quaternary mixtures of rare gases, Molecular Physics 12 (1967) 569–579.  doi:10.1080/00268976700100731.
  * [44] R. Kee, F. Rupley, J. Miller, Chemkin-ii: A Fortran chemical kinetics package for the analysis of gas-phase chemical kinetics.
  * [45] T. J. Poinsot, S. K. Lele, Boundary conditions for direct simulations of compressible viscous flows, Journal of Computational Physics 101 (1) (1992) 104 – 129.  doi:https://doi.org/10.1016/0021-9991(92)90046-2.
  * [46] J. D. Anderson, Hypersonic and High Temperature Gas Dynamics, McGraw-Hill series in aeronautical and aerospace engineering, American Institute of Aeronautics and Astronautics, 2000.
  * [47] A. Corrigan, D. Williams, A. Kercher, Weak formulation of a conservation law in reference space, Tech. rep., U.S. Naval Research Laboratory (March 2020).
  * [48] A. Kercher, A. Corrigan, D. Kessler, The moving discontinuous Galerkin finite element method with interface condition enforcement for compressible viscous flows, arXiv preprint arXiv:2002.12740 https://arxiv.org/pdf/2002.12740.pdf.
  * [49] E. Toro, Riemann solvers and numerical methods for fluid dynamics: a practical introduction, Springer Science & Business Media, 2013.
  * [50] F. Bassi, S. Rebay, An implicit high-order discontinuous Galerkin method for the steady state compressible Navier-Stokes equations, Computational fluid dynamics’98 (1998) 1226–1233.
  * [51] F. Bassi, S. Rebay, GMRES discontinuous Galerkin solution of the compressible Navier-Stokes equations, in: Discontinuous Galerkin Methods, Springer, 2000, pp. 197–208.
  * [52] F. Bassi, S. Rebay, Numerical evaluation of two discontinuous galerkin methods for the compressible navier–stokes equations, International journal for numerical methods in fluids 40 (1-2) (2002) 197–207.
  * [53] S. Gottlieb, C. Shu, E. Tadmor, Strong stability-preserving high-order time discretization methods, SIAM review 43 (1) (2001) 89–112.
  * [54] R. Spiteri, S. Ruuth, A new class of optimal high-order strong-stability-preserving time discretization methods, SIAM Journal on Numerical Analysis 40 (2) (2002) 469–491.
  * [55] Q. Zhang, C.-W. Shu, Error estimates to smooth solutions of Runge–Kutta discontinuous Galerkin methods for scalar conservation laws, SIAM Journal on Numerical Analysis 42 (2) (2004) 641–666.
  * [56] G. Strang, On the construction and comparison of difference schemes, SIAM Journal on Numerical Analysis 5 (3) (1968) 506–517.
  * [57] H. Atkins, C. Shu, Quadrature-free implementation of discontinuous galerkin methods for hyperbolic equations, ICASE Report 96-51, 1996, Tech. rep., NASA Langley Research Center, nASA-CR-201594 (August 1996).
  * [58] H. L. Atkins, C.-W. Shu, Quadrature-free implementation of discontinuous Galerkin method for hyperbolic equations, AIAA journal 36 (5) (1998) 775–782.
  * [59] T. Johnsen, T. Colonius, Implementation of WENO schemes in compressible multicomponent flow problems, Journal of Computational Physics 219 (2) (2006) 715 – 732.  doi:https://doi.org/10.1016/j.jcp.2006.04.018.
  * [60] E. Johnsen, On the treatment of contact discontinuities using WENO schemes, Journal of Computational Physics 230 (24) (2011) 8665 – 8668.  doi:https://doi.org/10.1016/j.jcp.2011.08.017.
  * [61] E. Ching, Y. Lv, P. Gnoffo, M. Barnhardt, M. Ihme, Shock capturing for discontinuous Galerkin methods with application to predicting heat transfer in hypersonic flows, Journal of Computational Physics 376 (2019) 54–75.
  * [62] P.-O. Persson, J. Peraire, Sub-cell shock capturing for discontinuous Galerkin methods, in: AIAA (Ed.), 44th AIAA Aerospace Sciences Meeting and Exhibit, 2006, AIAA-2006-112.  doi:10.2514/6.2006-112.
  * [63] R. Bauer, Discontinuous galerkin methods for ordinary differential equations, Ph.D. thesis, University of Colorado at Denver (June 1995).
  * [64] M. Delfour, W. Hager, F. Trochu, Discontinuous galerkin methods for ordinary differential equations, Math. Comp. 36 (1) (1981) 455–473.
  * [65] E. Hairer, G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, Vol. 14, 1996.  doi:10.1007/978-3-662-09947-6.
  * [66] D. G. Goodwin, H. K. Moffat, r. l. speth, cantera: an object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes, version 2.4.0 (2018).  doi:10.5281/zenodo.1174508.
URL http://www.cantera.org

  * [67] J. E. Shephard, Explosion dynamics laboratory: Shock and detonation toolbox - 2018 version, jES 9-19-2018 (2018).
URL http://shepherd.caltech.edu/EDL/PublicResources/sdt/

  * [68] A. Cuoci, A. Frassoldati, T. Faravelli, E. Ranzi, Numerical modeling of laminar flames with detailed kinetics based on the operator-splitting method, Energy & Fuels 27 (12) (2013) 7730–7753.  doi:10.1021/ef4016334.
URL https://doi.org/10.1021/ef4016334

  * [69] H. Wang, X. You, A. V. Joshi, S. G. Davis, A. Laskin, F. L. Egolfopoulos, C. K., Usc mech version ii. high-temperature combustion reaction model of h2/co/c1-c4 compounds (2007).
URL http://ignis.usc.edu/USC_mech_II.htm

  * [70] J. F. Haas, B. Sturtevant, Interaction of weak-shock waves, Journal of Fluid Mechanics 181 (1) (1987) 41 – 76.
  * [71] J. Quirk, S. Karni, On the dynamics of a shock-bubble interaction, Journal of Fluid Mechanics 318 (1996) 129–163.  doi:10.1017/S0022112096007069.
  * [72] A. Marquina, P. Mulet, A flux-split algorithm applied to conservative models for multicomponent compressible flows, Journal of Computational Physics 185 (1) (2003) 120 – 138.  doi:https://doi.org/10.1016/S0021-9991(02)00050-5.
  * [73] H. Terashima, G. Tryggvason, A front-tracking/ghost-fluid method for fluid interfaces in compressible flows, Journal of Computational Physics 228 (11) (2009) 4012 – 4037.  doi:https://doi.org/10.1016/j.jcp.2009.02.023.
  * [74] C. K. Westbrook, Chemical kinetics of hydrocarbon oxidation in gaseous detonations, Combustion and Flame 46 (1982) 191 – 210.  doi:https://doi.org/10.1016/0010-2180(82)90015-3.
  * [75] M. H. Lefebvre, J. Weber, J. W., E. S. Oran, Proceedings of the IUTAM Symposium (B. Deshaies and L. F. da Silva, eds.) (30) (1998) 347 – 358.
  * [76] M. Ó Conaire, H. J. Curran, J. M. Simmie, W. J. Pitz, C. K. Westbrook, A comprehensive modeling study of hydrogen oxidation, International Journal of Chemical Kinetics 36 (11) (2004) 603–622.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/kin.20036, doi:10.1002/kin.20036.
  * [77] C. Geuzaine, J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering (79(11)) (2009) 1310–1331.
  * [78] K. E. Niemeyer, N. J. Curtis, C.-J. Sung, pyjac: Analytical jacobian generator for chemical kinetics, Computer Physics Communications 215 (2017) 188 – 203.  doi:https://doi.org/10.1016/j.cpc.2017.02.004.
  * [79] F. Perini, E. Galligani, R. D. Reitz, An analytical jacobian approach to sparse reaction kinetics for computationally efficient combustion modeling with large reaction mechanisms, Energy & Fuels 26 (8) (2012) 4804–4822.  doi:10.1021/ef300747n.
URL https://doi.org/10.1021/ef300747n

  * [80] D. Kuzmin, R. Löhner, Flux-Corrected Transport: Principles, Algorithms, and Applications, 2nd Edition, Applied mathematical sciences, Springer-Netherlands, Springer-Verlag New York Inc, 2012.  doi:https://doi.org/10.1007/978-94-007-4038-9.
  * [81] A. A. Madja, Compressible fluid flow and systems of conservation laws in several space variables / A. Majda, Vol. 4 of Applied mathematical sciences, Springer-Verlag, Springer-Verlag New York Inc, 1984.

##  Appendix A Jacobian of the chemically reacting source term

For brevity we do not present details on the differentiation of the complex
reactions found in large chemical mechanisms and instead refer the reader to
[78] and [79]. Therefore, in the following we assume the derivatives of
chemical source term, $\omega_{i}$, with respect to concentrations and
temperature, $\frac{\partial\omega_{i}}{\partial C_{j}}$ and
$\frac{\partial\omega_{i}}{\partial T}$ are known for all $i$ and $j$.

The perturbation of the source term is found via the chain rule

$\mathcal{\delta S}=\mathcal{S}^{\prime}\delta
y=\frac{\partial\mathcal{S}}{\partial T}T^{\prime}\delta
y+\sum_{j=1}^{N_{s}}\frac{\partial\mathcal{S}}{\partial
C_{j}}C^{\prime}_{j}\delta y.$ (A.1)

The first term, $\frac{\partial\mathcal{S}}{\partial T}T^{\prime}\delta y$,
relies on temperature and its implicit relationship to the conserved state,
whereas the second term, $\frac{\partial\mathcal{S}}{\partial
C_{j}}C_{j}^{\prime}\delta y$, can be extracted from the species elements
corresponding to the perturbation,

$C_{j}^{\prime}\delta
y=\left(\begin{array}[]{c}0,\dots,0,\delta_{j,1},\dots\delta_{j,N_{s}}\end{array}\right)\left(\begin{matrix}\delta\rho
v_{1}\\\ \vdots\\\ \delta\rho v_{d}\\\ \delta\rho e_{t}\\\ \delta C_{1}\\\
\vdots\\\ \delta
C_{N_{s}}\end{matrix}\right)\begin{cases}\delta_{j,i}=1&i=j\\\
\delta_{j,i}=0&\text{otherwise}\end{cases}=\delta C_{j}.$ (A.2)

The temperature jacobian, $T^{\prime}$, in Equation A.1 is extracted by first
applying the chain rule to the definition total energy,

$\rho e_{t}=\rho u+\frac{1}{2}\sum_{k=1}^{d}\rho v_{k}v_{k}$ (A.3)

and grouping elements of the state perturbation together. Since temperature is
implicitly related to total energy through the internal energy, $\rho
u=\sum_{j=1}^{n_{s}}W_{j}C_{j}\sum_{k=1}^{n_{p}}a_{jk}T^{k}$, we extract the
temperature perturbation by using the derivatives of internal energy,
$\frac{\partial\rho u}{\partial
T}=\sum_{j=1}^{n_{s}}W_{j}C_{j}\sum_{k=1}^{n_{p}}ka_{jk}T^{k-1}$ and
$\frac{\partial\rho u}{\partial C_{j}}=W_{j}\sum_{k=1}^{n_{p}}a_{jk}T^{k}$,
combined with the perturbations of the kinetic energy,

$\displaystyle\delta\rho e_{t}$ $\displaystyle=$
$\displaystyle\frac{\partial\rho u}{\partial T}\delta
T+\sum_{j=1}^{n_{s}}\frac{\partial\rho u}{\partial C_{j}}\delta
C_{j}+\frac{1}{2}\left(\sum_{k=1}^{d}\left(\delta\rho
v_{k}\right)v_{k}+\sum_{k=1}^{d}\rho v_{k}\delta v_{k}\right).$ (A.4)

The internal energy derivatives are can be evaluated if the polynomial fits
for internal energy have been defined such that they are differentiable. The
perturbations of total energy, $\rho e_{t}$, species, $C_{k}$, and momentum,
$\delta\rho v_{i}$, all come from the perturbed state. In order to move
forward, we need the velocity perturbation which comes from the following
definition,

$\displaystyle v_{k}$ $\displaystyle=$ $\displaystyle\frac{\rho v_{k}}{\rho}$
$\displaystyle\delta v_{k}$ $\displaystyle=$ $\displaystyle\frac{\delta\rho
v_{k}}{\rho}-\frac{\rho v_{k}}{\rho^{2}}\delta\rho,$ (A.5)

and the perturbation of density,

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{n_{s}}W_{i}C_{i}$ $\displaystyle\delta\rho$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{n_{s}}W_{i}\delta C_{i}.$ (A.6)

By substituting A.5 and A.6 into A.4 we arrive at a total energy perturbation
dependent only on the state and temperature perturbations,

$\displaystyle\delta\rho e_{t}$ $\displaystyle=$
$\displaystyle\frac{\partial\rho u}{\partial T}\delta
T+\sum_{j=1}^{n_{s}}\frac{\partial\rho u}{\partial C_{j}}\delta
C_{j}+\sum_{k=1}^{d}\left(\delta\rho
v_{k}\right)v_{k}-\frac{1}{2\rho}\sum_{k=1}^{d}v_{k}v_{k}\sum_{j=1}^{n_{s}}W_{j}\delta
C_{j}.$

The temperature perturbation is then solved for,

$\displaystyle\delta T$ $\displaystyle=T^{\prime}\delta y=-\frac{\partial\rho
u}{\partial T}^{-1}$ $\displaystyle\left(\sum_{k=i}^{n_{s}}\frac{\partial\rho
u\left(C_{j,}T\right)}{\partial C_{i}}\delta
C_{i}+\sum_{k=1}^{d}\left(\delta\rho v_{k}\right)v_{k}\right.$ (A.7)
$\displaystyle\left.-\frac{1}{2\rho}\sum_{k=1}^{d}v_{k}v_{k}\sum_{j=1}^{n_{s}}W_{j}\delta
C_{j}-\delta\rho e_{t})\right)$

and the temperature jacobian is extracted by grouping like terms,

$\displaystyle\delta T$ $\displaystyle=T^{\prime}\delta y=-\frac{\partial\rho
u}{\partial T}^{-1}$
$\displaystyle\left(\begin{array}[]{c}v_{1},\dots,v_{d},-1,\frac{\partial\rho
u}{\partial
C_{1}}-\frac{1}{2\rho}\sum_{k=1}^{d}v_{k}v_{k}W_{1},\dots,\frac{\partial\rho
u}{\partial
C_{n_{s}}}-\frac{1}{2\rho}\sum_{k=1}^{d}v_{k}v_{k}W_{n_{s}}\end{array}\right)\left(\begin{matrix}\delta\rho
v_{1}\\\ \vdots\\\ \delta\rho v_{d}\\\ \delta\rho e_{t}\\\ \delta C_{1}\\\
\vdots\\\ \delta C_{n_{s}}\end{matrix}\right)$ (A.9) $\displaystyle
T^{\prime}$ $\displaystyle=-\frac{\partial\rho u}{\partial
T}^{-1}\left(\begin{array}[]{c}v_{1},\dots,v_{d},-1,\alpha_{1},\dots,\alpha_{n_{s}}\end{array}\right),$
(A.11)

where

$\alpha_{i}=\frac{\partial\rho u}{\partial
C_{i}}-\frac{1}{2\rho}\sum_{k=1}^{d}v_{k}v_{k}W_{i}.$ (A.12)

We now use A.11 and A.2 together to form the source term jacobian,

$\displaystyle\mathcal{S}^{\prime}=$
$\displaystyle\frac{d\mathcal{S}}{dT}T^{\prime}+\sum_{j=1}^{n_{s}}\frac{d\mathcal{S}}{dC_{j}}C_{j}^{\prime}=$
$\displaystyle-\frac{\partial\rho u}{\partial T}^{-1}\left(\begin{matrix}0\\\
\vdots\\\ 0\\\ 0\\\ \frac{d\omega_{1}}{dT}\\\ \vdots\\\
\frac{d\omega_{n_{s}}}{dT}\end{matrix}\right)\left(\begin{array}[]{c}v_{1},\dots,v_{d},-1,\alpha_{1},\dots,\alpha_{n_{s}}\end{array}\right)$
(A.16) $\displaystyle+\sum_{j=1}^{n_{s}}\left(\begin{matrix}0\\\ \vdots\\\
0\\\ \frac{\partial\omega_{1}}{\partial C_{j}}\\\ \vdots\\\
\frac{\partial\omega_{n_{s}}}{\partial
C_{j}}\end{matrix}\right)\left(\begin{array}[]{c}0,\dots,0,0,\delta_{j,1},\dots\delta_{j,n_{s}}\end{array}\right)\begin{cases}\delta_{j,i}=1&i=j\\\
\delta_{j,i}=0&\text{otherwise}\end{cases}$

which gives

$\mathcal{S}^{\prime}=\left(\begin{array}[]{ccccccc}0&\dots&0&0&0&\dots&0\\\
\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&\dots&0&0&0&\dots&0\\\
0&\dots&0&0&0&\dots&0\\\
-\beta_{1}v_{1}&\dots&-\beta_{1}v_{d}&\beta_{1}&\chi_{11}&\dots&\chi_{n_{s}1}\\\
\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
-\beta_{n_{s}}v_{1}&\dots&-\beta_{n_{s}}v_{d}&\beta_{N_{s}}&\chi_{1n_{s}}&\dots&\chi_{n_{s}n_{s}}\end{array}\right)$
(A.17)

where $\beta_{i}=\frac{\partial\rho u}{\partial T}^{-1}\frac{d\omega_{i}}{dT}$
, $\chi_{ij}=-\alpha_{i}\beta_{j}+\frac{\partial\omega_{j}}{\partial C_{i}}$,
and finally

$\delta\mathcal{S}=\left(\begin{array}[]{ccccccc}0&\dots&0&0&0&\dots&0\\\
\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\ 0&\dots&0&0&0&\dots&0\\\
0&\dots&0&0&0&\dots&0\\\
-\beta_{1}v_{1}&\dots&-\beta_{1}v_{d}&\beta_{1}&\chi_{11}&\dots&\chi_{n_{s}1}\\\
\vdots&\ddots&\vdots&\vdots&\vdots&\ddots&\vdots\\\
-\beta_{n_{s}}v_{1}&\dots&-\beta_{n_{s}}v_{d}&\beta_{N_{s}}&\chi_{1n_{s}}&\dots&\chi_{n_{s}n_{s}}\end{array}\right)\left(\begin{matrix}\delta\rho
v_{1}\\\ \vdots\\\ \delta\rho v_{n}\\\ \delta\rho e_{t}\\\ \delta C_{1}\\\
\vdots\\\ \delta C_{N_{s}}\end{matrix}\right)=\left(\begin{matrix}0\\\
\vdots\\\ 0\\\ 0\\\ \frac{\partial\omega_{1}}{\partial T}\delta
T+\sum_{i=1}^{n_{s}}\frac{\partial\omega_{1}}{\partial C_{k}}\delta C_{k}\\\
\vdots\\\ \frac{\partial\omega_{N_{s}}}{\partial T}\delta
T+\sum_{i=1}^{n_{s}}\frac{\partial\omega_{N_{s}}}{\partial C_{i}}\delta
C_{i}\end{matrix}\right)$ (A.18)

##  Appendix B Non-reflecting boundary conditions

The Riemann invariants for the calorically perfect Euler equations with
boundary normal, $n=\left(n_{1},\dots,n_{d}\right)$, and normal velocity
$v_{n}=\sum_{k=1}^{d}v_{k}\cdot n_{k}$ are specified in Kuzmin et al. [80]

$w=\left(v_{n}-\frac{2c}{\bar{\gamma}-1},\bar{c_{v}}\textup{log}\left(\frac{p}{\rho^{\bar{\gamma}}}\right),v_{\xi},v_{\eta},v_{n}+\frac{2c}{\bar{\gamma}-1}\right)^{T},$
(B.1)

where $v_{\xi}$ and $v_{\eta}$ are the velocity components tangential to the
boundary surface and $c=\sqrt{\gamma RT}$ is the speed of sound where $R$ is
calculated from Equation (2.23). The specific heat at constant volume is
$\bar{c}_{v}=\bar{c}_{p}-R$. These invariants have the corresponding
eigenvalues

$\Lambda=\textup{diag}\left\\{v_{n}-c,v_{n},v_{n},v_{n},v_{n}+c\right\\}$
(B.2)

which determine what is specified at a subsonic inflow and outflow boundary
conditions.

For a subsonic outlet $v_{n}>0$ and $M<1$ so $\Lambda_{1}=v_{n}-c$ is negative
and requires $w_{1}$ to be specified from the prescribed state, $y_{\infty}$,
and $w_{2}$-$w_{5}$ come from the interior state, $y^{+}$,

$\displaystyle w_{1}$ $\displaystyle=$ $\displaystyle
v_{n,\infty}-\frac{2c_{\infty}}{\bar{\gamma}_{\infty}-1}$ (B.3) $\displaystyle
w_{2}$ $\displaystyle=$
$\displaystyle\bar{c_{v}}^{+}\textup{log}\left(\frac{p^{+}}{\left(\rho^{\bar{\gamma}}\right)^{+}}\right)$
(B.4) $\displaystyle w_{3}$ $\displaystyle=$ $\displaystyle v_{\xi}^{+}$ (B.5)
$\displaystyle w_{4}$ $\displaystyle=$ $\displaystyle v_{\eta}^{+}$ (B.6)
$\displaystyle w_{5}$ $\displaystyle=$ $\displaystyle
v_{n}^{+}+\frac{2c^{+}}{\bar{\gamma}^{+}-1}.$ (B.7)

The species mass fractions and mole fractions are assumed to be from the
interior state which gives $c_{\infty}=\sqrt{\gamma^{+}R^{+}T_{\infty}}$ where
$T_{\infty}$ is specified but the gas properties come from the interior. It is
often inconvenient to specify temperature and velocity, which is required for
$w_{1}$ at an outflow, unless the boundary is specified far enough away to not
interfere with the flow. For a downstream outflow condition, the pressure,
$p_{\infty}$, can be specified and the temperature, $T_{\infty}$, can be
calculated from a constant entropy process assuming the mole fractions are
constant. Here entropy is

$S=\sum_{i=1}^{n_{s}}X_{i}\sum_{k=0}^{n_{p}}s_{ik}T^{k}-X_{i}R\log
X_{i}-X_{i}R\ln{\frac{p}{p_{atm}}}$ (B.8)

where $s_{ik}$ are the polynomial coefficients for the species specific
entropies. The constant entropy process temperature, $T_{\partial}$, is
computed such that the following is satisfied to machine precision for a given
initial temperature guess:

$\delta
T=\frac{\sum_{i=1}^{n_{s}}X_{i}\sum_{k=0}^{n_{p}}s_{ik}\left(T^{+}\right)^{k}+R\ln{\frac{p_{\partial}}{p^{+}}}-\sum_{i=1}^{n_{s}}X_{i}\sum_{k=0}^{n_{p}}s_{ik}T_{\partial}^{k}}{\frac{\partial
S}{\partial T}},$ (B.9)

where $\delta T$ is the temperature decrement corresponding to Newton’s method
and

$\frac{\partial S}{\partial
T}=\sum_{i=1}^{n_{s}}X_{i}\sum_{k=1}^{K}ks_{ik}T_{\partial}^{k-1}.$ (B.10)

With $T_{\partial}$ known, $\bar{\gamma}_{\partial}$ and $\rho_{\partial}$ can
be calculated and the incoming characteristic can be specified,

$w_{1}=w_{5}-\frac{4}{\bar{\gamma_{\infty}}-1}\sqrt{\frac{\bar{\gamma_{\infty}}p_{\infty}}{\rho_{\infty}}}.$
(B.11)

The thermodynamic state for the boundary is then given by

$\displaystyle\rho^{*}$ $\displaystyle=$
$\displaystyle\left(\frac{c_{\infty}^{2}}{\bar{\gamma}^{+}}e^{-\frac{w_{2}}{\bar{c_{v}}^{+}}}\right)^{\frac{1}{\bar{\gamma}^{+}-1}}$
(B.12) $\displaystyle c^{*}$ $\displaystyle=$
$\displaystyle\left(\gamma^{+}-1\right)\left(w_{5}-w_{1}\right)$ (B.13)
$\displaystyle p^{*}$ $\displaystyle=$
$\displaystyle\frac{\rho^{*}\left(c^{*}\right)^{2}}{\bar{\gamma}^{+}}$ (B.14)
$\displaystyle T^{*}$ $\displaystyle=$
$\displaystyle\frac{p^{*}}{\rho^{*}\bar{R}^{+}}$ (B.15) $\displaystyle v^{*}$
$\displaystyle=$
$\displaystyle\frac{w_{1}+w_{5}}{2}\left(n_{1},\dots,n_{d}\right)+v^{+}-v_{n}\left(n_{1},\dots,n_{d}\right)$
(B.16) $\displaystyle Y_{i}^{*}$ $\displaystyle=$ $\displaystyle Y_{i}^{+}.$
(B.17)

Likewise, the subsonic inflow requires $w_{1}-w_{4}$ to be specified,

$\displaystyle w_{1}$ $\displaystyle=$ $\displaystyle
v_{n}-\frac{2c_{\infty}}{\bar{\gamma}_{\infty}-1},$ (B.18) $\displaystyle
w_{2}$ $\displaystyle=$ $\displaystyle
c_{v,\partial}\textup{log}\left(\frac{p_{\infty}}{\left(\rho^{\bar{\gamma}}\right)_{\infty}}\right),$
(B.19) $\displaystyle w_{3}$ $\displaystyle=$ $\displaystyle v_{\xi,\infty,}$
(B.20) $\displaystyle w_{4}$ $\displaystyle=$ $\displaystyle v_{n,\infty},$
(B.21)

where all the gas properties come from the specified exterior state. The last
invariant comes from the interior,
$w_{5}=v_{n}^{+}+\frac{2c^{+}}{\bar{\gamma}^{+}-1}$. This gives the
thermodynamic state at the boundary as

$\displaystyle\rho^{*}$ $\displaystyle=$ $\displaystyle\rho_{\infty}$ (B.22)
$\displaystyle c^{*}$ $\displaystyle=$
$\displaystyle\left(\gamma_{\infty}-1\right)\frac{\left(w_{5}-w_{1}\right)}{4}$
(B.23) $\displaystyle p^{*}$ $\displaystyle=$
$\displaystyle\frac{\rho^{*}\left(c^{*}\right)^{2}}{\bar{\gamma_{\infty}}}$
(B.24) $\displaystyle T^{*}$ $\displaystyle=$
$\displaystyle\frac{p^{*}}{\rho^{*}\bar{R}_{\infty}}$ (B.25) $\displaystyle
v^{*}$ $\displaystyle=$
$\displaystyle\frac{w_{1}+w_{5}}{2}\left(n_{1},\dots,n_{d}\right)+v_{\infty}-v_{n}\left(n_{1},\dots,n_{d}\right)$
(B.26) $\displaystyle Y_{i}^{*}$ $\displaystyle=$ $\displaystyle Y_{i,\infty}$
(B.27)

The estimates for the thermodynamic state are used to construct the boundary
value,

$y^{*}\left(y^{+},y_{\infty},n^{+}\right)=\left(\rho^{*}v_{1}^{*},\ldots,\rho^{*}v_{n}^{*},\rho
e_{t},\frac{\rho^{*}Y_{1}^{*}}{W_{i}},\ldots,\rho Y_{n_{s}}^{*}\right),$
(B.28)

where

$\rho e_{t}=\rho u+\frac{1}{2}\sum_{k=1}^{d}\rho^{*}v_{k}^{*}v_{k}^{*}$ (B.29)

with

$\rho
u=\sum_{i=1}^{n_{s}}W_{i}C_{i}\sum_{k=0}^{n_{p}}a_{ik}\left(T^{*}\right)^{k}.$
(B.30)

##  Appendix C Material discontinuities in multi-component flows

A material discontinuity is defined as a discontinuity across which there is
no mass flow. The velocity and pressure are constant across the discontinuity
but other material quantities are not. In this appendix we analyze problems
involving material interfaces for the formulation presented in this work by
considering the non-reacting inviscid formulation of Equations (3.6)-(3.8)
where $\mathcal{F}^{v}\left(y,\nabla y\right)=\left(0,\dots,0\right)$ in
Equation (2.2) and $\mathcal{S}\left(y\right)=\left(0,\ldots,0\right)$ in
Equation (3.6). A discontinuous solution, in one dimension satisfies, the
inviscid form of Equations (3.6)-(3.8) if the jump in the flux is equal to the
product of the jump in the state and the material interface velocity [81],

$\mathcal{F}\left(y_{r}\right)-\mathcal{F}\left(y_{l}\right)=v_{s}\left(y_{r}-y_{l}\right),$
(C.1)

where $y_{r}$ is the state on the right of the discontinuity, $y_{l}$ is the
state on the left of the discontinuity, and $v_{s}$ is the material velocity
normal to the interface.

Below we introduce a material discontinuity by considering a one-dimensional
two species discontinuity at $x_{j}$ where the velocity and pressure are
constant and the temperature is discontinuous,

$\displaystyle v$ $\displaystyle=$ $\displaystyle\bar{v},$ $\displaystyle
C_{1}$ $\displaystyle=$ $\displaystyle\begin{cases}C_{1}^{0}&\text{if
}x>x_{j}\\\ 0&\text{otherwise}\end{cases},$ $\displaystyle C_{2}$
$\displaystyle=$ $\displaystyle\begin{cases}0&\text{if }x<x_{j}\\\
C_{2}^{0}&\text{otherwise}\end{cases},$ $\displaystyle T$ $\displaystyle=$
$\displaystyle\begin{cases}\eta\bar{T}&\text{if }x<x_{j}\\\
\bar{T}&\text{otherwise}\end{cases},$ $\displaystyle p$ $\displaystyle=$
$\displaystyle\bar{p.}$

The species with index $i=1$, species 1, has molecular weight $W_{1}$ and the
species with index $i=2$, species 2, has molecular weight $W_{2}$. The initial
fluid state from Equation (2.3) is therefore

$y\left(x,t=0\right)=\begin{cases}\left(W_{1}C_{1}^{0}\bar{v},\frac{1}{2}W_{1}C_{1}^{0}\bar{v}^{2}+W_{1}C_{1}^{0}\sum_{k=0}^{n_{p}}a_{1k}\bar{T}^{k},C_{1}^{0},0\right)&x>x_{j}\\\
\left(W_{2}C_{2}^{0}\bar{v},\frac{1}{2}W_{2}C_{2}^{0}\bar{v}^{2}+W_{2}C_{2}^{0}\sum_{k=0}^{n_{p}}a_{2k}\left(\eta\bar{T}\right)^{k},0,C_{2}^{0}\right)&\text{otherwise}\end{cases}.$
(C.2)

Substituting the fluid state from Equation (C.2) in Equation (C.1) we arrive
at the following condition

$\displaystyle W_{1}C_{1}^{0}\bar{v}^{2}-W_{2}C_{2}^{0}\bar{v}^{2}$
$\displaystyle=$ $\displaystyle
v_{s}\left(W_{1}C_{1}^{0}v-W_{2}C_{2}^{0}v\right),$ (C.3)
$\displaystyle\left(\frac{1}{2}W_{1}C_{1}^{0}v^{2}+C_{1}^{0}W_{1}\sum_{k=0}^{n_{p}}a_{1k}T^{k}+\bar{p}\right)\bar{v}$
$\displaystyle-\left(\frac{1}{2}W_{2}C_{2}^{0}v^{2}+C_{2}^{0}W_{2}\sum_{k=0}^{n_{p}}a_{2k}\left(\eta
T\right)^{k}+\bar{p}\right)\bar{v}$ $\displaystyle=$ $\displaystyle
v_{s}\left(\frac{1}{2}W_{1}C_{1}^{0}v^{2}+W_{1}C_{1}^{0}\sum_{k=0}^{n_{p}}a_{1k}T^{k}\right.$
(C.4)
$\displaystyle\left.-\frac{1}{2}W_{2}C_{2}^{0}v^{2}-W_{2}C_{2}^{0}\sum_{k=0}^{n_{p}}a_{2k}\left(\eta
T\right)^{k}\right),$ $\displaystyle C_{1}^{0}\bar{v}$ $\displaystyle=$
$\displaystyle v_{s}\left(C_{1}^{0}\right),$ (C.5) $\displaystyle-
C_{2}^{0}\bar{v}$ $\displaystyle=$ $\displaystyle
v_{s}\left(-C_{2}^{0}\right).$ (C.6)

Therefore a material discontinuity where velocity and pressure are constant
and the temperature is discontinuous satisfies Equations (3.6)-(3.8) with
$v_{s}=\bar{v}$. A diagram of the space-time solution is shown in Figure C.1a.

We now present the effect of a linear discretization on the same two species
discontinuity. Using the notation from Abgrall and Karni [25], the inviscid
non-reacting conservation equations can be written as

$\displaystyle\delta\left(\rho v\right)+\nu\Delta\left(\rho v^{2}+p\right)=$
$\displaystyle 0,$ (C.7) $\displaystyle\delta\left(\rho
e_{t}\right)+\nu\Delta\left(\left(\rho e_{t}+p\right)v\right)=$ $\displaystyle
0,$ (C.8)
$\displaystyle\delta\left(C_{i}\right)+\nu\Delta\left(C_{i}v\right)=$
$\displaystyle 0\text{ for }i=1\dots n_{s},$ (C.9)

where the inviscid form of Equation (2.1) has been linearized with respect to
space and time. Here,
$\delta\left(\right)=\left(\right)_{j}^{n+1}-\left(\right)_{j}^{n}$ denotes
the temporal change of the state, $\Delta\left(\right)$ denotes spatial
variation across the interface
$\Delta\left(\right)=\left(\right)_{j}^{n}-\left(\right)_{j-1}^{n}$, and
$\nu=\frac{\Delta t}{\Delta x}$ where $\Delta t$ is the chosen time step and
$\Delta x$ is the spatial distance across the interface. As in the [25] and
[28], $\Delta\left(\right)=\left(\right)_{j}^{n}-\left(\right)_{j-1}^{n}$ is
assumed to have the following properties:
$\Delta\left(ab+c\right)=a\left(b\right)_{j}^{n}-a\left(b\right)_{j-1}^{n}+\left(c\right)_{j}^{n}-\left(c\right)_{j-1}^{n}$
for $a$ constant, and $b$ and $c$ not constant. We are deriving in the context
of first order approximations without reconstruction, however, it should be
noted that these properties do not always apply to all schemes. The material
interface is initially between two nodes, $j$ and $j-1$, as depicted at time
$t^{n}$ in Figure C.1b. Specifically, the initial flow state at $t^{n}$ is

$\displaystyle y_{j}^{n}$ $\displaystyle=$
$\displaystyle\left(W_{1}C_{1}^{0}\bar{v},\frac{1}{2}W_{1}C_{1}^{0}\bar{v}^{2}+W_{1}C_{1}^{0}\sum_{k=0}^{n_{p}}a_{1k}\bar{T}^{k},C_{1}^{0},0\right),$
(C.10) $\displaystyle y_{j-1}^{n}$ $\displaystyle=$
$\displaystyle\left(W_{2}C_{2}^{0}\bar{v},\frac{1}{2}W_{2}C_{2}^{0}\bar{v}^{2}+W_{2}C_{2}^{0}\sum_{k=0}^{n_{p}}a_{2k}\left(\eta\bar{T}\right)^{k},0,C_{2}^{0}\right).$
(C.11)

For simplification purposes, we define the initial concentration of species 2
in terms of the initial concentration of species 1 through the constant
initial pressure conditions, $p_{j}^{n}=p_{j-1}^{n}=\bar{p}$, and the equation
of state, Equation (2.6),

$R^{o}\bar{T}C_{1}^{0}=R^{o}\eta\bar{T}C_{2}^{0}\rightarrow
C_{2}^{0}=\frac{C_{1}^{0}}{\eta}.$ (C.12)

The species conservation, Equation (C.9), gives the concentrations at
$t^{n+1}$ in terms of the initial species 1 concentration,

$\displaystyle C_{1,j}^{n+1}$ $\displaystyle=$ $\displaystyle
C_{1}^{0}-\nu\bar{v}C_{1}^{0},$ (C.13) $\displaystyle C_{2,j}^{n+1}$
$\displaystyle=$
$\displaystyle\nu\bar{v}C_{2}^{0}=\nu\bar{v}\frac{C_{1}^{0}}{\eta}.$ (C.14)

Equations (C.13) and (C.14) show that there is numerical mixing of the species
at time $t^{n+1}$ and node $j$, as depicted in Figure C.1b. This is a
departure from the exact solution that satisfies the interface condition,
depicted in Figure C.1a, and we continue in this section by examining the
effect that the numerical mixing of the species concentrations has on the
stability of the material interface.

(a) The speed of the material interface is the slope of the trajectory
$v_{s}=\bar{v}=\angle$. (b) Numerical diffusion of material interface at
discrete times.

Figure C.1: Diagrams of moving material interfaces.

Substituting the concentrations from Equations (C.13) and (C.14) into Equation
(2.4) gives the density at $t^{n+1}$,

$\rho_{j}^{n+1}=\left(W_{1}+\nu\bar{v}\left(\frac{W_{2}}{\eta}-W_{1}\right)\right)C_{1}^{0}.$
(C.15)

Using Equations (C.13)-(C.15) and Equation (2.4) for the initial density and
substituting into Equation (C.7) reveals that the velocity remains constant,

$\displaystyle\left(W_{1}+\nu\bar{v}\left(\frac{W_{2}}{\eta}-W_{1}\right)\right)C_{1}^{0}v^{n+1}-W_{1}C_{1}^{0}\bar{v}+\nu\bar{v}^{2}\left(W_{1}C_{1}^{0}-W_{2}\frac{C_{1}^{0}}{\eta}\right)$
$\displaystyle=$ $\displaystyle 0\rightarrow v_{j}^{n+1}=\bar{v}.$ (C.16)

Using Equations (C.13)-(C.16) we consider the change in total energy to
analyze the stability of the material interface. We derive a relationship for
kinetic energy by multiplying Equation (C.7) by $\frac{1}{2}\bar{v}$,

$\displaystyle\left(\frac{1}{2}\rho_{j}^{n+1}\bar{v}^{2}-\frac{1}{2}\rho_{j}^{n}\bar{v}^{2}+\nu\bar{v}\left(\frac{1}{2}\rho_{j}^{n}\bar{v}^{2}-\frac{1}{2}\rho_{j-1}^{n}\bar{v}^{2}\right)\right)$
$\displaystyle=$ $\displaystyle 0,$ (C.17)

and we derive a relationship for pressure by noting $pv$ is constant across
the interface at $t^{n}$,

$\displaystyle\Delta\left(pv\right)$ $\displaystyle=$ $\displaystyle 0.$
(C.18)

Combining Equation (C.17) and Equation (C.18) with Equation (C.8) we remove
the kinetic energy, contained in $\rho e_{t}=\rho v^{2}/2+\rho u$, and the
pressure term to yield a linear relationship for the internal energy across
the interface,

$\displaystyle\delta\rho u+\nu v\Delta\left(\rho u\right)$ $\displaystyle=$
$\displaystyle 0.$ (C.19)

We substitute Equation (2.8), Equation (C.12), and Equations (C.13)-(C.15) in
Equation (C.19) and arrive at $n_{p}$ expressions for the temperature at time
$t^{n+1}$ by collecting like terms,

$\displaystyle T_{j}^{n+1}$ $\displaystyle=$
$\displaystyle\bar{T}\frac{\left(W_{1}a_{11}-\nu\bar{v}\left(W_{1}a_{11}-W_{2}a_{21}\right)\right)}{W_{1}\left(1-\nu\bar{v}\right)a_{11}+\frac{W_{2}}{\eta}\left(\nu\bar{v}\right)a_{21}}$
$\displaystyle\vdots$ $\displaystyle\vdots$ $\displaystyle\vdots$
$\displaystyle\left(T_{j}^{n+1}\right)^{n_{p}}$ $\displaystyle=$
$\displaystyle\left(\bar{T}\right)^{n_{p}}\frac{\left(W_{1}a_{1n_{p}}-\nu\bar{v}\left(W_{1}a_{1N_{p}}-W_{2}a_{2N_{p}}\eta^{n_{p}-1}\right)\right)}{W_{1}\left(1-\nu\bar{v}\right)a_{1n_{p}}+\frac{W_{2}}{\eta}\left(\nu\bar{v}\right)a_{2n_{p}}}.$
(C.20)

Finally, the change in pressure is given as

$\displaystyle p_{j}^{n+1}-p_{j}^{n}$ $\displaystyle=$ $\displaystyle
R^{o}T_{j}^{n+1}\left(C_{1}^{0}-\nu\bar{v}C_{1}^{0}+\nu\bar{v}\frac{C_{1}^{0}}{\eta}\right)-R^{o}\bar{T}C_{1}^{0}.$
(C.21)

From analyzing Equations (C.20) and (C.21) we come to the same conclusions to
those of Jenny et al [21], that pressure oscillations,
$p_{j}^{n+1}-p_{j}^{n}\neq 0$, do not exist if one of the following conditions
is true

  1. 1.

The temperature is continuous, $\eta=1$.

  2. 2.

The contact discontinuity remains grid aligned, $\nu\bar{v}=1$.

  3. 3.

The contact discontinuity is stationary, $\bar{v}=0$.

  4. 4.

The internal energies are linear, $n_{p}=1$, with respect to temperature _and_
the species are the same across the interface, i.e., molecular weights are
constant across the interface, $W_{1}=W_{2}$, and the internal energies are
the same across the interface, $a_{1k}=a_{2k}$.

For condition (1), the numerical mixing of species concentrations, Equations
(C.13) and (C.14), inside the cell does not cause a pressure oscillation as
both species are at the same temperature despite having different internal
energies.

When $\eta\neq 1$ the temperature is discontinuous and stabilization, e.g.,
artificial viscosity, would be required if (2)-(4) were not satisfied.
Satisfaction of condition (2) would requires an interface fitting method [22,
23] that dynamically fits a priori unknown discontinuities and is therefore
beyond the scope of this manuscript. Condition (3) is a trivial case.
Condition (4) applies to ideal gases that have a linear relationship between
temperature and internal energy and are assumed to be the same species in all
regions of the flow.

Applying the same linearization to Equation (2.31) we can arrive at a similar
relationship for the internal energy based on $\bar{\gamma}$ and $p$,

$\displaystyle\frac{p_{j}^{n+1}}{\bar{\gamma}_{j}^{n+1}-1}$ $\displaystyle=$
$\displaystyle\frac{\bar{p}}{\bar{\gamma_{1}}-1}-\nu
v\left(\frac{\bar{p}}{\bar{\gamma_{1}}-1}-\frac{\bar{p}}{\bar{\gamma}_{2}-1}\right),$
(C.22)

where $\bar{\gamma}_{1}$ and $\bar{\gamma}_{2}$ are the known specific heat
ratios of the right and left hand side based on $C_{1}^{0}$ at temperature
$\bar{T}$ and $C_{2}^{0}$ at temperature $\eta\bar{T}$, respectively. The
equivalent process for this formulation would be to use the definition of
pressure and $\bar{\gamma}$ in terms of known concentrations, $C_{1}^{n+1}$,
$C_{2}^{n+1}$, $C_{i}^{0}$, and $C_{2}^{0}$, and temperatures, $\eta\bar{T}$
and $\bar{T}$, to solve for $T^{n+1}$. This results in similar nonlinear
relationships for temperature but instead from the $h_{i}$ polynomials. It
follows that the same stability properties found for formulation presented in
this work also apply when the formulation is written in terms of specific
heats that are computed exactly, i.e., not frozen. However, when
$\bar{\gamma}$ is frozen and the temperature is not solved for exactly,
pressure oscillations occur according to the analysis of [25].