{'fraction_non_alphanumeric': 0.03559704863725342, 'fraction_numerical': 0.015810871856648097, 'mean_word_length': 4.638648327390049, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Conservative Discontinuous Galerkin Discretization for the Total Energy Formulation of the Reacting Navier Stokes Equations

This paper describes the total energy formulation of the compressible reacting Navier-Stokes equations which is solved numerically using a fully conservative discontinuous Galerkin finite element method (DG). Previous applications of DG to the compressible reacting Navier-Stokes equations required nonconservative fluxes or stabilization methods in order to suppress unphysical oscillations in pressure that led to the failure of simple test cases. In this paper, we demonstrate that material interfaces with a temperature discontinuity result in numerical unphysical pressure oscillations if the species internal energy is nonlinear with respect to temperature. We demonstrate that a temperature discontinuity is the only type of material interface that results in unphysical pressure oscillations for a conservative discretization of the total energy formulation. Furthermore, we demonstrate that unphysical pressure oscillations will be generated at any material interface, including material interfaces with at which the temperature is continuous, if the thermodynamics are frozen during the temporal integration of the conserved state. Additionally, we demonstrate that the oscillations are amplified if the specific heat at constant pressure is incorrectly evaluated directly from the NASA polynomial expressions. Instead, the mean value, which we derive in this manuscript, should be used to compute the specific heat at constant pressure. This can reduce the amplitude of, but not prevent, unphysical oscillations where the species concentrations numerically mix. We then present solutions to several test cases using the total energy formulation and demonstrate spurious pressure oscillations were not generated for material interfaces if the temperature is continuous and that it is better behaved than frozen thermodynamic formulations if the temperature is discontinuous.


Background
The discontinuous Galerkin finite element method (DG) [1,2,3,4,5,6,7,8,9,10,11,12], has been applied to the field of computational fluid dynamics with great success over the past two decades. The method is fully conservative, able to achieve high-order accuracy on unstructured grids, and has a welldeveloped theory of adjoint consistency [13,14,12,15], which makes it a powerful tool for adjoint-based optimization. Furthermore, an extension of DG, the Moving Discontinuous Galerkin Method with Interface Condition Enforcement (MDG-ICE) [16,17], maintains high order accuracy for flows with interfaces and sharp gradients. Recent work has shown success in modeling combustion systems using DG [18,19]. In these previous works, nonphysical pressure oscillations were generated for continuous profiles of species and temperature with constant velocity and pressure. These oscillations have been attributed to variations in the thermodynamic properties of multicomponent gases and it was previously concluded that any fully conservative Godunov-type scheme would be unable to maintain a pressure equilibrium across the material fronts [20]. This was shown conclusively for formulations that used variable ratio of specific heats, γ, where γ was a function of species concentrations [21,22].
A nonconservative flux, referred to as the double flux method, was implemented to avoid these oscillations in multicomponent reacting flows [23,24,18,19]. The method assumes consistent fluid thermal properties through a material interface, and thereby breaks energy conservation to achieve the desired stability. This method has been successfully applied to multidimensional reacting flow problems, including detonations, and has been shown to be stable [25,19]. Other methods, such as increasing the size of the state and solving for additional transport equations, have also been employed to avoid these unwanted oscillations in multicomponent flows [26,27,28,29,30].
In this work, we present the total energy formulation of the reacting Navier-Stokes equations to simulate multicomponent flows that is suitable for conservative DG discretizations. The formulation avoids unwanted pressure oscillations without the need for nonconservative fluxes or other stabilization methods by solving for temperature at each degree of freedom such that the internal energy of the conserved fluid is consistent with the mixture-averaged, temperature-based polynomial expression for internal energy. This formulation is used to approximate solutions to various test cases with material interfaces. The results confirm that pressure oscillations are only generated at material interfaces if the temperature is discontinuous. Additionally, we demonstrate that unphysical pressure oscillations will be generated at any material interface, including material interfaces where the temperature is continuous, if the thermodynamics are frozen during the temporal integration of the conserved state. Furthermore, we present verification for a steady flame, demonstrating the ability of the method to model the desired physics in combustion systems with high order accuracy without generating pressure oscillations.

The Total Energy Formulation
The nonlinear conservation law, given in strong form, defined for piecewise smooth, R m -valued functions y, and gradient ∇y, is ∂y ∂t + ∇ · F (y, ∇y) − S (y) = 0 in Ω, (2.1) y (·, t 0 ) − y 0 = 0 in Ω, (2.2) n · F (y, ∇y) − n · F Γ (y Γ , ∇y) = 0 in Γ, is defined in terms of the convective flux F c (y), which is only a function of the state y, and viscous flux F v (y, ∇y), which is a function of the state and the gradient, ∇y. The reacting Navier-Stokes flow state variable is given by y = (ρv 1 , . . . , ρv n , ρe t , C i , . . . , C ns ) ∈ R m , (2.5) where m = n + n s + 1, n is the number of spatial dimensions, n s is the number of thermally perfect species, ρ is density, (v 1 , . . . , v n ) ∈ R n is velocity, e t is the specific total energy, and C i is the concentration of species i. The density is calculated from the concentrations as where W i is the molecular weight of species i.
The k-th spatial convective flux component is given by The pressure is calculated from the equation of state, where T is the temperature and R o is the universal gas constant, 8314.4621 J/Kmol/K. The total energy in Eq. (2.5) is related to kinetic and internal energy by where ρu is a mass weighted sum of thermally perfect species specific internal energies that are n p -order polynomials with respect to temperature, The k-th spatial component of the viscous flux is given by where λ is the thermal conductivity, τ ∈ R n×n is the viscous stress tensor, h i is the specific enthalpy of species i, and (V i1 , . . . , V in ) ∈ R n is the diffusion velocity of species i. The k-th spatial component of the viscous stress tensor is given by where µ is the dynamic viscosity coefficient and (x 1 , . . . , x n ) ∈ R n are spatial coordinates. The transport properties are calculated using mixture averaged properties. The species diffusion velocity is calculated using the mixture averaged diffusion, where D i is the mixture averaged diffusion coefficient of species i from [31], p atm = 101325 Pa, X j is the mole fraction of species j, D ij is the diffusion coefficient of species i to species j, andW is the mixture molecular weight,W = ρ/ ns i=1 C i . The mole fractions can be calculated directly from concentrations, , (2.15) defined in terms of where µ i and µ j are the species specific viscosities for species i and j respectively. The Mathur model [33] is used to calculate conductivity, where λ i is the conductivity of species i. Finally, the source term is given by where ω i is the production rate of species i. The production rate comes from the sum of the progress reaction rates from any arbitrary number of reactions and reaction types.

The Ratio of Specific Heats Formulation
Here we present how the ratio of specific heats, γ, can be used in the total energy formulation, described in Section 2.1, and we present the necessary constraints to keep the two formulations consistent. We seek a formulation analogous to the calorically perfect gas formulation, ρu = p/ (γ − 1). To do so we use the definition of internal energy in terms of enthalpy and pressure, ρu = ρh − p. (2.18) Here the enthalpy is where Y i is the mass fraction of species i, Y i = W i C i /ρ, C p,i is the specific heat at constant pressure of species i, and h i is the species specific enthalpy polynomial of temperature that is degree n p , h i = np k=0 a ik T k + R o T . We reduce the definition of internal energy to achieve the equivalent formulation that contains a similar expression to ρu = p γ−1 by introducing the mean value from reference temperature, T 0 , to current temperature, T , ofC p,i andC p , where h 0 i is the species specific enthalpy at T 0 . Using the following definitions the inviscid total energy conservation without reactions becomes is eliminated from Eq. (2.25) by fixing T 0 to 0 K and multiplying the non-reacting inviscid form of the species conservation equations from Eq. (2.1) by W i h 0 i and summing over all species conservation equations. Eq. (2.25) is equivalent to the non-reacting inviscid form of the conservation of energy from Eq. (2.1). Eq. (2.25) has the same form of the compressible Euler equations, and therefore is convenient for routines that require the ratio of specific heats, e.g., characteristic boundary conditions [34]. However, those routines require the flow to be non-reacting with constant thermodynamic properties asγ is assumed to be constant.

Material Discontinuities in Multicomponent 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 section we exact solutions for problems involving material interfaces by considering the non-reacting inviscid formulation of Eqs. 3) if the jump in the flux is equal to the product of the jump in the state and the material interface velocity [35], 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, v =v, 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 Eq. (2.5) is therefore Substituting the fluid state from Eq. (3.2) in Eq. (3.1) we arrive at the following condition Therefore a material discontinuity where velocity and pressure are constant and the temperature is discontinuous satisfies Eqs. (2.1)-(2.3) with v s =v. A diagram of the space-time solution is shown in Fig. 3.1a.
We now present the effect of a linear discretization on the same two species discontinuity. Using the notation from Abgrall and Karni [23], the inviscid non-reacting conservation equations can be written as ∆x where ∆t is the chosen time step and ∆x is the spatial distance across the interface. The material interface is initially between two nodes, j and j − 1, as depicted at time t n in Fig. 3.1b. Specifically, the initial flow state at t n is 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 n j = p n j−1 =p, and the equation of state, Eq. (2.8), The species conservation, Eq. (3.9), gives the concentrations at t n+1 in terms of the initial species 1 concentration, Eqs. (3.13) and (3.14) show that there is numerical mixing of the species at time t n+1 and node j, as depicted in Fig. 3.1b. This is a departure from the exact solution that satisfies the interface condition, depicted in Fig. 3.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.

Species 2
Species 1 The speed of the material interface is the slope of the trajectory vs =v = ∠. Substituting the concentrations from Eqs. (3.13) and (3.14) into Eq. (2.6) gives the density at t n+1 , Using Eqs. (3.13)-(3.15) and Eq. (2.6) for the initial density and substituting into Eq. (3.7) reveals that the velocity remains constant, Using Eqs. (3.13)-(3.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 Eq. (3.7) by 1 2v , and we derive a relationship for pressure by noting pv is constant across the interface at t n , ∆ (pv) = 0. We substitute Eq. (2.10), Eq. (3.12), and Eqs. (3.13)-(3.15) in Eq. (3.19) and arrive at n p expressions for the temperature at time t n+1 by collecting like terms, (3.20) Finally, the change in pressure is given as From analyzing Eqs. (3.20) and (3.21) we come to similar conclusions to those of Jenny et al [36], that pressure oscillations, p n+1 j − p n j = 0, do not exist if one of the following conditions is true 1. The temperature is continuous, η = 1. 2. The contact discontinuity remains grid aligned, νv = 1. 3. The contact discontinuity is stationary,v = 0. 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, Eqs. (3.13) and (3.14), inside the cell does not cause a pressure oscillation as both species are at the same temperature despite having different internal energies. When η = 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 [16,17,37,38] 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 Eq. (2.25) we can arrive at a similar relationship for the internal energy based onγ and p, whereγ 1 andγ 2 are the known specific heat ratios of the right and left hand side based on C 0 1 at temperaturē T and C 0 2 at temperature ηT , respectively. The equivalent process for the ratio of specific heats formulation would be to use the definition of pressure andγ in terms of known concentrations, C n+1 i , and C 0 2 , and temperatures, ηT andT , 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 the total energy formulation apply to the ratio of specific heats formulation.

Numerical Methods
In order to discretize the Eqs. (2.1)-(2.3), the DG method assumes that Ω can be subdivided into a mesh T h = {κ}, consisting of disjoint cells κ such thatΩ = ∪ κ∈T hκ . We also consider I h , consisting of interior interfaces defined between pairs of cells, so that for each ∈ I h , there exists a pair of κ + , κ − ∈ T h such that = κ + ∩ κ − . A discrete subspace, consisting of piecewise polynomials of degree p, is defined over T h , where h (y + , y − , n) is the numerical flux, choose to be the HLLC approximate Riemann [39], see also Appendix B of [19], is the average viscous flux, and δ (y + , y − , n) is a penalty term that is required for stability and is implemented via the BR2 formulation [40,41,42]. On the exterior interfaces, ∂κ / ∈ I h , the numerical flux is defined consistently with the imposed boundary condition Eq. (2.3) and The DG space semi-discretization is integrated temporally with a strong-stability-preserving Runge-Kutta method (SSP-RK2) [43]. This discretization of convection and diffusion is combined with a stiff ODE solver for the reaction term via Strang splitting [44].

Implementation of Thermodynamics
It is important to note thatC p,i , Eq. (2.20), andC p , Eq. (2.21), from Section 2.2 are not equivalent to the evaluations of C p,i and C p from polynomial expressions, e.g., the analytic form of NASA's polynomial representations [45]. Rather,C p,i andC p can be viewed as the mean value of C p,i and C p , whereC p approaches C p as T −T o → 0. This difference can be shown mathematically by using the following polynomial definitions of specific heat at constant pressure, and mixture averaged C p , We arrive at the total enthalpy in polynomial form by integrating of Eqs. , By substituting Eq. (4.5) in Eq. (2.21), we arrive atC p in terms of the C pi polynomial coefficients and temperature,C Therefore,C p and C p are only equivalent if C p is constant with respect to temperature, i.e., n p = 0.  In Section 5 three different methods are used to evaluate the thermodynamics in the convective flux, F c k (y), for various test cases. Table 1 shows the steps required to calculate the convective fluxes. The first column corresponds to the total energy formulation that requires the solution of a nonlinear equation for the temperature. The solution ensures equivalency between the computed internal energy and the internal energy given by the evaluation of the corresponding polynomial expression, i.e., find T such that where the internal energy, ρu = ρe t − 1 2 n k=1 ρv k v k , and concentrations, C i , are computed from the known state, y. In this work the temperature is computed such that the following is satisfied to machine precision for a given an initial temperature: where δT is the temperature decrement corresponding to Newton's method and is the partial derivative of internal energy with respect to temperature. In practice, we observe that the temperature converges within five nonlinear iterations. Once the temperature and state are known, pressure is computed by evaluating Eq. (2.8). The second column of Table 1 corresponds to a frozen thermodynamics evaluation of convective flux where the ratio of specific heats formulation is used andγ andC p are held constant throughout each time step. The third row for the frozen thermodynamics is split into two options.
The left corresponds to usingC p in the evaluation of the flux and the right corresponds to using C p , which is evaluated from NASA polynomials expressions. Calculate and freeze newγ fromC p and R Ifγ is frozen as presented in the second column in Table 1 then γ n+1 j in Eq. (3.22) isγ 1 . A pressure oscillation is now generated regardless of the type of material interface, The freezing of thermodynamics would therefore require other methods to identify and stabilize the nonphysical pressure oscillations even in smooth regions of the flow where η = 1. Previous work showed that pressure oscillations develop even in the presence of continuous temperature profiles which was due to the freezing of thermodynamics [24,18,25,19,30]. In the following section we explore the magnitude of these pressure oscillations due to frozen thermodynamics and compare to the behavior of the total energy formulation.

Test Cases and Verification
In this section we solve several test cases with the total energy formulation and verify the analysis from Section 3. Furthermore, we compare the stability of the total energy formulation to the ratio of specific heats formulation using both frozen C p and frozenC p . We then verify the total energy formulation by comparing to a one-dimensional premixed hydrogen flame simulation generated by Cantera's free-flame solver [46]. In the test cases below, the time step is restricted by the Courant-Friedrichs-Lewy number, CFL, defined as where p is the polynomial degree and c is the speed of sound, defined as c = √ γRT . Here γ is the ratio of specific heats which is defined as γ = Cp Cp−R for both the total energy formulation and the ratio of specific heats formulation with frozen C p and γ =γ for the ratio of specific heats formulation with frozenC p .

Discontinuities
We approximate exact solutions corresponding to the discontinuities described in Section 3 by solving the non-reacting, inviscid formulation of the Navier-Stokes Eqs.  Table 1.

Species Discontinuities at Constant Temperature
A species discontinuity at constant temperature, pressure, and velocity, is imposed with the following initial conditions v = 10 m/s, Figs. (5.1a), (5.1b), and (5.1c) show the species mass fractions for the DG(p = 1), DG(p = 2), and DG(p = 3) solutions using the total energy formulation after one cycle. All three solutions present numerical overshoots and mixing of the species mass fractions across the discontinuities. Some of these numerical instabilities cause the mass fractions to be greater than one or negative. As expected, the higher order solutions are more oscillatory for the species mass fractions. Furthermore, larger overshoots are present for the right hand side discontinuity as compared to the left hand side discontinuity.    2) shows the pressure after one cycle for the total energy formulation and ratio of specific heats formulations using frozenC P and frozen C p . The frozen C p solution has pressure oscillations that are on the order of three percent error. Small perturbations of pressure caused by freezingC p are shown in Fig. (5.2b). The error caused by the frozenC p formulation is reduced as the approximation order is increased from DG(p = 1) to DG(p = 3). In contrast, the pressure for each solution using the total energy formulation maintains a flat profile, which is the expected result from Section 3. The pressure in the total energy solutions never exceed an error of 1e−7 atm.
The corresponding temperature solutions are shown in Figs. (5.3a) and (5.3b). The frozen C p solution gives temperature fluctuations on the order of 1 K whereas the frozenC p solution fluctuates less than 0.025 K. The total energy solution remains flat and does not exceed 1e−4 K from the expected 300 K.        The larger time steps exacerbate the instability introduced by freezingC p to exceed 0.001 atm which is still an order of magnitude less than the error of the frozen C p solutions. The solutions for the total energy formulation is unaffected by the time step and therefore the corresponding results are not shown. Figure (5.4b) shows the pressure solution of the DG(p = 1) solution for the frozenC p formulation and total energy formulation after 100 cycles, i.e. t = 1 s. The pressure oscillations for the 100 cycle solution using the frozenC p formulation grow in time to be on the order of the one cycle frozen C p formulation, whereas the total energy formulation remains constant after 100 cycles.

Species Discontinuities and a Continuous Temperature Variation
A sinusoidal variation in temperature is introduced to the initial profile described in Section (5.1.1). The initial conditions are given as follows v = 10 m/s, Figs. (5.5a), (5.5b), and (5.5c) show the species mass fractions for DG(p = 1), DG(p = 2), and DG(p = 3) solutions after one cycle. Similar to the previous test case, the three solutions present numerical overshoots and mixing of the species mass fractions across the discontinuities. Here the left side discontinuity has larger overshoots in the higher temperature region as compared to the right hand side discontinuity. Figure (5.6) shows the pressure for all three solutions using the total energy formulation and ratio of specific heats formulation with frozen C p and frozenC p after one cycle. Fig. (5.6a) shows the large oscillations like the previous test case for the frozen C p formulation and smaller oscillations using the frozen C p formulation. The oscillations caused by freezingC p are reduced as the accuracy of the approximation is increased from DG(p = 1) to DG(p = 3). Again, the pressure for each solution using the total energy formulation maintains a flat pressure profile, even with a spatially varying temperature profile. This is expected as the temperature is continuous through the domain. Figure (5.7) shows the computed temperature from solutions using the total energy formulation and the ratio of specific heats formulation with frozen C p and frozenC p after one cycle. For the total energy and ratio of specific heats formulation with frozenC p the temperature is within 1e−4 K of the analytical sinusoidal solution. The ratio of specific heats formulation with frozen C p is also shown and departs from the analytical result with the largest deviation of 5 K in the higher temperature region

Species Discontinuities with Temperature Discontinuities
The following test case uses the initial conditions described in Section 5.1.1, except the initial temperature profile is piecewise constant instead of constant. The initial conditions are given as follows v = 10 m/s, Figs. 5.8a, 5.8b, and 5.8c show the species mass fractions for DG(p = 1), DG(p = 2), and DG(p = 3) solutions after one cycle. Similar to the previous test case, the three solutions present numerical overshoots and mixing of the species mass fractions across the discontinuities. Figs. 5.9a and 5.9b show the pressure and temperature, respectively, after one cycle for the total energy formulation and the ratio of specific heats formulation with frozenC P . The ratio of specific heats formulation with frozen C p simulation fails before one complete cycle. The solution corresponding to frozen C p is shown after 100 time steps. The pressure for each solution using the total energy formulation causes pressure oscillations that are an order of magnitude less than the frozenC p simulations. The oscillations in the total energy formulation are expected based on the discussion in Section 3. Both the total energy formulation and the ratio of specific heats formulation with frozenC p produce overshoots and undershoots at the temperature discontinuities. The temperature oscillations associated with the ratio of specific heats formulation with frozenC p are larger than the oscillations corresponding to the total energy formulation.  Figs. 5.10a and 5.10b show the temperature and pressure solutions, respectively, for the total energy formulation and the ratio of specific heats formulation with frozenC P after 100 cycles, i.e. t = 1 s. The temperature profiles become more diffuse with the larger number of cycles (see Fig. 5.9b for the comparison of 1 cycle). The pressure oscillations for the total energy formulation cause the ambient pressure to fall below 1 atm. This departure from ambient was improved by increasing the approximation order from DG(p = 1) to DG(p = 3). The pressure oscillations for the 100 cycle solution using the frozenC p formulation grow in time regardless of approximation order (see Fig. 5.9a for the one cycle solution for frozenC P ). Furthermore, the frozenC p formulation oscillations are an order of magnitude larger than the oscillations present in the total energy formulation.

Continuous Simulations 5.2.1. Thermal Bubble
Here we present the one dimensional thermal bubble test case previously presented by [19].

One-Dimensional Premixed Flame
We use the total energy formulation with detailed chemical kinetics and transport in the viscous, reacting formulation of Eqs. (2.1)-(2.3) to approximate the solution to a one-dimensional premixed flame using the Hydrogen-Air chemistry from [47]. The chosen domain is 6 cm in length with h = 5e−5 m. We initialize the domain to an ambient pressure of 1 bar, temperature of 300 K, X H2 = 0.188, X O2 = 0.170, and X N2 = 0.642. A small section on the right hand side of the domain is initialized to the fully reacted state from Cantera's homogeneous constant pressure reactor simulation [46].
The right hand side boundary is fixed to constant pressure with remaining variables interpolated from the interior state. The left hand side boundary is a characteristic wave boundary condition that 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 oscillation, which is supported by the analysis in Section 3. The initial pressure oscillations leave the system as the reaction front diffuses into the unreacted region, which eventually creates a stable propagating flame.

Conclusion
We have presented the total energy formulation and the ratio of specific heats formulation of the reacting Navier-Stokes equations and examined pressure oscillations for material interfaces with discontinuities in species and temperature. We derived that if the species internal energy is nonlinear with respect to temperature then pressure oscillations are only generated at material interfaces with discontinuities in both species and temperature. Pressure oscillations are not generated at a material interface with a species discontinuity if the temperature is continuous through the interface unless the thermodynamics are frozen during the temporal integration of the conserved state. The total energy formulation 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. The converged temperature can then be use to evaluate other thermodynamic quantities so that they are consistent with the conserved state. As such, the total energy formulation can be integrated with a fully conservative method without generating pressure oscillations in regions of continuous pressure and temperature.
We presented several test cases to demonstrate the generation of unphysical pressure oscillations. We demonstrated that the pressure oscillations from the ratio of specific heats formulation with frozenC p did not reach the magnitudes previously reported in literature [23,24,18,25,19]. 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,C p = C p , the oscillations did reach the same levels as previously reported, indicating that the severity of the pressure oscillations is dependent on the correct evaluation of thermodynamic quantities. Regardless, the pressure should be constant across the material interfaces and whenC p is frozen the magnitude of pressure oscillations will grow in time as the solution evolves. Additionally, the expected numerical pressure oscillations due to temperature discontinuities at material interfaces were less for the total energy formulation than for the ratio of specific heats formulation with frozenC p .
Previous work showed that advecting continuous species and temperature profiles caused unstable oscillations [18,19]. We analyzed continuous profiles by presenting the thermal bubble test case for a multicomponent mixture of hydrogen and oxygen from [19]. The total energy formulation solution for this test case did not generate pressure oscillations. The solutions reached the magnitudes of pressure oscillations found in the previous work only when the ratio of specific heats with frozen C p was used. When the ratio of specific heats with frozenC p was used no pressure oscillations were generated due to the lack of numerical mixing for an advected continuous profile of species and temperature.
Finally, we presented the results for a fully reacting Navier-Stokes simulation of a continuous onedimensional flame. The solution, computed using a fully conservative DG method without additional stabilization, compared well with the Cantera solution. Flame speeds for both the DG(p = 1) and DG(p = 2) solutions were consistent with the Cantera flame speed. Discrepancies between the Cantera solution and the DG(p = 1) solution are not observed in the DG(p = 2) solution, indicating the DG(p = 2) solution better resolved the flame.
The total energy formulation does not require non-conservative methods or additional stabilization in smooth regions of the flow. This is an attractive feature since both non-conservative methods and artificial stabilization have associated inherent costs. Non-conservative methods often involve additional cell and face loops that increase the computational complexity. Additionally, artificial stabilization can prevent high order methods from achieving the formal order accuracy associated with the polynomial space of the approximate solution. We also demonstrated that the total energy formulation is equivalent to the ratio of specific heats formulation in non-reacting regions of the flow where thermodynamic quantities are not changing, which is a suitable alternative and may be convenient for methods that require the internal energy to be defined according to the ratio of specific heats, e.g., characteristic boundary conditions.