Entropy-Stable Numerical Methods for Third-Order Extended Thermodynamics of Reacting Gas Mixtures

1. Introduction

Chemically reacting gas mixtures are encountered in a wide variety of scientific and engineering applications, including combustion systems, hypersonic vehicles, atmospheric re-entry, plasma-assisted propulsion, astrophysical flows, and high-temperature chemical processing. The accurate numerical prediction of these flows requires the simultaneous resolution of compressible fluid dynamics, heat transfer, species diffusion, and chemical reaction kinetics, all of which evolve over widely separated temporal and spatial scales. The resulting multiscale nature of reacting flows presents significant challenges for continuum modeling and numerical simulation, particularly when the gas departs substantially from local thermodynamic equilibrium [1,2,3,4].

The classical Navier–Stokes–Fourier (NSF) equations have formed the foundation of continuum fluid mechanics for more than a century. Derived through the Chapman–Enskog expansion of the Boltzmann equation, these equations successfully describe weakly nonequilibrium gases in which molecular relaxation occurs much faster than macroscopic transport processes [1,3]. However, the constitutive assumptions underlying the NSF equations, namely Fourier’s law of heat conduction and Newton’s law of viscosity, imply instantaneous propagation of thermal and mechanical disturbances because of their parabolic mathematical structure. Such assumptions become increasingly inadequate for rarefied gases, hypersonic flows, rapid transient processes, microscale heat transfer, and chemically reacting mixtures characterized by finite relaxation times [4,5,6].

These limitations motivated the development of Extended Thermodynamics (ET), which extends the set of thermodynamic state variables to include dissipative fluxes such as heat flux, dynamic pressure, and viscous stresses as independent variables possessing their own balance equations. The modern formulation of Rational Extended Thermodynamics (RET), pioneered by Müller and subsequently developed by Müller and Ruggeri, establishes a mathematically rigorous framework in which the governing equations constitute a symmetric hyperbolic system compatible with the entropy principle [7,8]. Unlike the classical NSF equations, RET predicts finite propagation speeds for thermal and mechanical disturbances while preserving the fundamental principles of continuum thermodynamics.

The theoretical foundations of RET are closely linked to Grad’s moment method, which derives macroscopic transport equations directly from kinetic theory through finite moment approximations of the Boltzmann equation [9]. Subsequent developments by Levermore [10], Struchtrup [4], and others have significantly improved the mathematical understanding of hyperbolic moment systems, entropy closures, and nonequilibrium gas dynamics. These advances have demonstrated that higher-order moment theories provide more accurate descriptions of rarefied and high-speed gas flows than the classical NSF equations while remaining computationally tractable.

The application of RET to chemically reacting mixtures was systematically developed by Kremer and Müller [11], who established a thermodynamically consistent framework for reacting gases based on entropy production and chemical affinity. Their formulation incorporated reaction source terms into the entropy balance and demonstrated how chemical reactions modify the constitutive structure of transport processes. Later studies further investigated multicomponent reacting gases, diffusion processes, and kinetic-theory foundations of chemically reacting mixtures [4,6,12]. Nevertheless, most available models remain limited to first-order constitutive equations and linear transport coefficients, making them unsuitable for strongly nonequilibrium environments where nonlinear interactions between dissipative fluxes become important.

Recently, increasing attention has been devoted to higher-order formulations of nonequilibrium thermodynamics capable of describing finite-amplitude transport phenomena. In particular, higher-order entropy expansions provide a systematic approach for incorporating nonlinear corrections to heat conduction, bulk viscosity, and diffusion while maintaining compatibility with the second law of thermodynamics [6,8]. Such formulations are especially relevant for high-temperature reacting gases, where chemical reactions, volumetric relaxation, and heat transfer interact over comparable time scales, giving rise to nonlinear thermo-chemical coupling that cannot be captured by first-order theories.

Building upon these developments, our recent work, Third-Order Extended Thermodynamics of Chemically Reacting Gas Mixtures [13], introduced a third-order entropy expansion for chemically reacting gas mixtures within the framework of Rational Extended Thermodynamics. The entropy density was generalized to include cubic contributions in the heat flux, dynamic pressure, and species diffusion fluxes, leading to nonlinear constitutive relations and generalized Maxwell–Cattaneo evolution equations. The resulting theory predicts effective thermal conductivity and bulk viscosity that depend explicitly on the magnitude of the dissipative fluxes and the chemical affinity, thereby extending the work of Kremer and Müller [11] to nonlinear nonequilibrium regimes.

Despite these theoretical advances, comparatively little progress has been made in the numerical approximation of third-order Extended Thermodynamics models. Hyperbolic relaxation systems present several computational challenges, including nonlinear wave propagation, stiff source terms, finite relaxation times, and the need to preserve the entropy inequality at the discrete level. Conventional finite-volume methods designed for the Euler or Navier–Stokes equations frequently fail to maintain positivity of thermodynamic variables or produce nonphysical oscillations in the vicinity of strong discontinuities when applied to nonequilibrium systems [14,15].

Over the last three decades, entropy-stable numerical methods have emerged as one of the most successful approaches for solving nonlinear hyperbolic conservation laws. Tadmor [16,17] introduced the concept of entropy-conservative numerical fluxes, providing the mathematical foundation for constructing discretizations that satisfy a discrete entropy inequality. These ideas have subsequently been extended to entropy-stable finite-volume, finite-difference, and discontinuous Galerkin (DG) methods for the Euler equations, magnetohydrodynamics, shallow-water systems, and compressible Navier–Stokes equations [18,19,20,21,22]. Entropy-stable methods have demonstrated remarkable robustness in capturing shocks while preserving nonlinear stability and high-order accuracy.

For reacting flows, additional numerical difficulties arise from the stiffness of chemical source terms. Implicit–explicit (IMEX) Runge–Kutta methods have become a widely adopted strategy for integrating hyperbolic relaxation systems because they combine explicit treatment of convective transport with implicit integration of stiff relaxation processes [23,24]. Together with asymptotic-preserving techniques, IMEX methods ensure stable integration even when the relaxation times become several orders of magnitude smaller than the characteristic flow time scales [25]. These approaches have proven highly effective for kinetic equations, radiative transfer, and relaxation systems, but have not yet been systematically developed for third-order Extended Thermodynamics of chemically reacting mixtures.

Another important aspect of numerical approximation concerns the preservation of the mathematical structure of the governing equations. Since RET possesses a convex entropy function that symmetrizes the governing equations, numerical schemes should preserve this entropy structure in order to maintain stability and physical admissibility. Entropy-variable formulations, split-form discretizations, and entropy-conservative interface fluxes provide an effective mechanism for constructing such schemes while maintaining high-order spatial accuracy [16,19,20]. Incorporating these ideas into the third-order RET framework offers the possibility of developing computational methods that remain robust under strong shocks, rapid chemical reactions, and highly nonequilibrium conditions.

The objective of the present work is therefore to develop entropy-stable numerical methods for the third-order Extended Thermodynamics model introduced in Third-Order Extended Thermodynamics of Chemically Reacting Gas Mixtures [13]. The governing hyperbolic relaxation system is discretized using conservative finite-volume methods equipped with entropy-consistent numerical fluxes and positivity-preserving reconstruction procedures. Stiff relaxation and chemical source terms are integrated using implicit–explicit Runge–Kutta methods to ensure numerical stability across a broad range of relaxation parameters. In addition, a discontinuous Galerkin formulation based on entropy variables and split-form discretization is developed to achieve high-order spatial accuracy while preserving the discrete entropy inequality.

The proposed numerical methodology is validated through several representative benchmark problems, including reacting shock-tube flows, thermal relaxation waves, and chemically induced nonequilibrium transport. Numerical predictions of Mach number, temperature, dynamic pressure, heat flux, and species mass fractions are compared with the analytical equilibrium limit and with the corresponding Navier–Stokes–Fourier solutions. The results demonstrate that the proposed entropy-stable algorithms accurately capture finite-speed heat propagation, nonlinear thermo-chemical coupling, and relaxation phenomena while maintaining excellent stability in regimes where conventional methods exhibit excessive numerical dissipation or nonphysical oscillations.

The remainder of this paper is organized as follows. Section 2 summarizes the governing equations of the third-order Extended Thermodynamics model and its hyperbolic relaxation structure. Section 3 develops the entropy-stable finite-volume formulation. Section 4 presents the discontinuous Galerkin discretization together with entropy-variable formulations. Section 5 introduces the IMEX time-integration strategy for stiff thermo-chemical source terms. Section 6 presents numerical validation through reacting shock tubes, thermal relaxation waves, and strongly nonequilibrium reacting flows. Finally, Section 7 summarizes the principal findings and outlines future developments toward multidimensional simulations, adaptive mesh refinement, and high-performance computing implementations for realistic reacting gas mixtures.

2. Governing Equations of the Third-Order Extended Thermodynamics Model

This section summarizes the governing equations of the third-order Extended Thermodynamics (ET3) model for chemically reacting gas mixtures. The formulation follows the framework of Rational Extended Thermodynamics (RET), in which the balance equations are supplemented by evolution equations for dissipative fluxes. Unlike the classical Navier–Stokes–Fourier (NSF) equations, the present model admits finite propagation speeds for thermal and mechanical disturbances while remaining compatible with the entropy principle through a symmetric hyperbolic structure.

2.1. Conserved Variables

Consider a chemically reacting gas mixture consisting of $N_s$ chemical species occupying a spatial domain $\mathrm{\Omega }\subset {\mathbb{R}}^d$. The vector of conserved variables is

$$\mathbf{U}={\left(\rho ,\rho \mathbf{v},\rho E,\rho Y_1,\dots ,\rho Y_{N_s-1}\right)}^T,$$

(1)

where $\rho$ denotes the mass density, $\mathbf{v}$ is the velocity vector, E is the total specific energy, and $Y_k$ denotes the mass fraction of species k. The species mass fractions satisfy

$${\displaystyle \sum _{k=1}^{N_s}}{Y}_k=1.$$

(2)

The total specific energy is given by

$$E=e+{\displaystyle \frac{1}{2}}{\left|\mathbf{v}\right|}^2,$$

(3)

where e denotes the specific internal energy.

2.2. Macroscopic Balance Equations

The governing conservation equations are

Mass conservation

$${\displaystyle \frac{\partial \rho }{\partial t}}+\nabla ·\left(\rho \mathbf{v}\right)=0.$$

(4)

Momentum conservation

$${\displaystyle \frac{\partial \left(\rho \mathbf{v}\right)}{\partial t}}+\nabla ·\left(\rho \mathbf{v}\otimes \mathbf{v}+p\mathbf{I}+\mathbf{\Pi }\right)=0.$$

(5)

Energy conservation

$${\displaystyle \frac{\partial \left(\rho E\right)}{\partial t}}+\nabla ·\left[(\rho E+p)\mathbf{v}+\mathbf{\Pi }\mathbf{v}+\mathbf{q}\right]=0.$$

(6)

Species conservation

$${\displaystyle \frac{\partial \left(\rho Y_k\right)}{\partial t}}+\nabla ·\left(\rho Y_k\mathbf{v}+{\mathbf{J}}_k\right)={\omega }_k,k=1,\dots ,N_s-1,$$

(7)

where ${\omega }_k$ denotes the chemical production rate of species k, ${\mathbf{J}}_k$ is the diffusion flux, $\mathbf{q}$ is the heat flux, and $\mathbf{\Pi }$ denotes the nonequilibrium stress tensor.

The reaction source terms satisfy

$${\displaystyle \sum _{k=1}^{N_s}}{\omega }_k=0,$$

(8)

ensuring conservation of the total mixture mass.

2.3. Third-Order Nonequilibrium Variables

Unlike the classical Navier–Stokes–Fourier theory, the third-order Extended Thermodynamics model treats the dissipative quantities as independent thermodynamic variables. The complete state vector therefore becomes

$$\mathbf{W}={\left(\rho ,\mathbf{v},T,Y_k,\mathbf{q},\mathrm{\Pi },{\mathbf{J}}_k\right)}^T,$$

(9)

where T denotes the temperature, $\mathrm{\Pi }$ is the dynamic pressure, and $\mathbf{q}$ and ${\mathbf{J}}_k$ satisfy their own evolution equations.

The nonlinear constitutive relations obtained from the third-order entropy expansion introduce effective transport coefficients that depend explicitly upon the dissipative fluxes and chemical affinity.

2.4. Hyperbolic Relaxation Equations

The generalized Maxwell–Cattaneo equations governing the dissipative variables are written as

$${\tau }_q\left({\displaystyle \frac{\partial \mathbf{q}}{\partial t}}+\mathbf{v}·\nabla \mathbf{q}\right)+\mathbf{q}=-{\kappa }_{\mathrm{eff}}\nabla T+{\mathbf{S}}_q,$$

(10)

$${\tau }_{\mathrm{\Pi }}\left({\displaystyle \frac{\partial \mathrm{\Pi }}{\partial t}}+\mathbf{v}·\nabla \mathrm{\Pi }\right)+\mathrm{\Pi }=-{\mu }_{b,\mathrm{eff}}\nabla ·\mathbf{v}+S_{\mathrm{\Pi }},$$

(11)

and

$${\tau }_J\left({\displaystyle \frac{\partial {\mathbf{J}}_k}{\partial t}}+\mathbf{v}·\nabla {\mathbf{J}}_k\right)+{\mathbf{J}}_k=-D_{\mathrm{eff}}\nabla Y_k+{\mathbf{S}}_{J_k},k=1,\dots ,N_s-1.$$

(12)

Here ${\tau }_q$, ${\tau }_{\mathrm{\Pi }}$, and ${\tau }_J$ denote the characteristic relaxation times for heat conduction, bulk viscosity, and diffusion, respectively. The source terms ${\mathbf{S}}_q$, $S_{\mathrm{\Pi }}$, and ${\mathbf{S}}_{J_k}$ represent nonlinear thermo-chemical coupling introduced through the third-order entropy expansion.

2.5. Compact Hyperbolic Form

Collecting the conservation and relaxation equations, the governing system may be written compactly as

$${\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+\nabla ·\mathbf{F}\left(\mathbf{U}\right)=\mathbf{S}\left(\mathbf{U}\right),$$

(13)

where $\mathbf{F}$ denotes the conservative flux vector and $\mathbf{S}$ contains the relaxation and chemical source terms.

The resulting system constitutes a hyperbolic system with stiff source terms, making it particularly suitable for entropy-stable finite-volume and discontinuous Galerkin discretizations coupled with implicit–explicit (IMEX) time integration.

2.6. Entropy Structure

A central property of Rational Extended Thermodynamics is the existence of a strictly convex mathematical entropy function

$$\eta =\eta \left(\mathbf{U}\right),$$

(14)

together with the corresponding entropy flux $\mathbf{\Psi }$ satisfying

$${\displaystyle \frac{\partial \eta }{\partial t}}+\nabla ·\mathbf{\Psi }=\sigma ,\sigma \ge 0,$$

(15)

where $\sigma$ denotes the non-negative entropy production.

The entropy variables are defined as

$$\mathbf{V}={\displaystyle \frac{\partial \eta }{\partial \mathbf{U}}},$$

(16)

which symmetrize the governing equations and provide the mathematical foundation for entropy-conservative and entropy-stable numerical discretizations. The numerical methods developed in the following sections are constructed to preserve a discrete analogue of the entropy inequality while maintaining positivity of density, pressure, and species concentrations.

3. Entropy-Stable Finite-Volume Formulation

The governing equations derived in the previous section constitute a nonlinear hyperbolic system with stiff relaxation and chemical source terms. Numerical approximation of such systems requires schemes that simultaneously preserve conservation, maintain thermodynamic admissibility, and remain stable in the presence of shocks and rapid relaxation processes. Classical shock-capturing methods often introduce excessive numerical dissipation or violate the entropy inequality, leading to nonphysical oscillations and loss of robustness in strongly nonequilibrium flows.

To overcome these difficulties, we develop an entropy-stable finite-volume discretization based on entropy-conservative numerical fluxes supplemented with carefully designed dissipation operators. The resulting scheme satisfies a discrete analogue of the second law of thermodynamics while preserving positivity of the primary thermodynamic variables.

3.1. Finite-Volume Discretization

Consider a partition of the computational domain

$$\mathrm{\Omega }={\displaystyle \bigcup _{i=1}^N}{V}_i,$$

(17)

where $V_i$ denotes the ith control volume with measure $|V_i|$.

Integrating the governing equations over each control volume yields

$${\displaystyle \frac{d}{dt}}{\int }_{V_i}\mathbf{U}dV+{\int }_{\partial V_i}\mathbf{F}\left(\mathbf{U}\right)·\mathbf{n}dS={\int }_{V_i}\mathbf{S}\left(\mathbf{U}\right)dV,$$

(18)

where $\mathbf{n}$ denotes the outward unit normal vector.

Approximating the cell averages by

$${\mathbf{U}}_i={\displaystyle \frac{1}{|V_i|}}{\int }_{V_i}\mathbf{U}dV,$$

(19)

the semi-discrete finite-volume scheme becomes

$${\displaystyle \frac{d{\mathbf{U}}_i}{dt}}=-{\displaystyle \frac{1}{|V_i|}}\sum _{f\in \partial V_i}{\widehat{\mathbf{F}}}_f{A}_f+{\mathbf{S}}_i,$$

(20)

where $A_f$ denotes the area of face f and ${\widehat{\mathbf{F}}}_f$ is the numerical interface flux.

3.2. Entropy Variables

The entropy variables play a central role in constructing entropy-stable discretizations. They are defined by

$$\mathbf{V}={\displaystyle \frac{\partial \eta }{\partial \mathbf{U}}},$$

(21)

where $\eta \left(\mathbf{U}\right)$ is the convex mathematical entropy introduced in Section 2.

Using entropy variables transforms the governing equations into a symmetric hyperbolic system

$${\mathbf{A}}_0\left(\mathbf{U}\right){\displaystyle \frac{\partial \mathbf{V}}{\partial t}}+{\displaystyle \sum _{j=1}^d}{\mathbf{A}}_j\left(\mathbf{U}\right){\displaystyle \frac{\partial \mathbf{V}}{\partial x_j}}=\mathbf{S},$$

(22)

where the matrices ${\mathbf{A}}_0$ and ${\mathbf{A}}_j$ are symmetric and ${\mathbf{A}}_0$ is positive definite.

This symmetrization guarantees well-posedness and provides the mathematical foundation for entropy-stable numerical fluxes.

3.3. Entropy-Conservative Flux

Following Tadmor’s framework, the numerical interface flux ${\widehat{\mathbf{F}}}^{EC}$ is required to satisfy

$${({\mathbf{V}}_R-{\mathbf{V}}_L)}^T{\widehat{\mathbf{F}}}^{EC}={\mathrm{\Psi }}_R-{\mathrm{\Psi }}_L,$$

(23)

where

$$\mathrm{\Psi }={\mathbf{V}}^T\mathbf{F}-\mathbf{Q}$$

(24)

is the entropy potential and $\mathbf{Q}$ denotes the entropy flux.

Condition (23) guarantees that the finite-volume method neither creates nor destroys entropy in smooth regions.

For practical computations, the entropy-conservative flux is evaluated using logarithmic averages of the thermodynamic variables together with arithmetic averages of the velocity components.

3.4. Entropy-Stable Dissipation

Although entropy-conservative schemes preserve entropy exactly, they contain insufficient numerical dissipation for discontinuous solutions.

Entropy stability is obtained by adding matrix dissipation,

$$\widehat{\mathbf{F}}={\widehat{\mathbf{F}}}^{EC}-{\displaystyle \frac{1}{2}}\mathbf{D}({\mathbf{U}}_R-{\mathbf{U}}_L),$$

(25)

where

$$\mathbf{D}=\mathbf{R}\left|\mathbf{\Lambda }\right|{\mathbf{R}}^{-1}$$

(26)

is constructed from the eigendecomposition of the flux Jacobian,

$$\mathbf{A}={\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{U}}}=\mathbf{R}\mathbf{\Lambda }{\mathbf{R}}^{-1}.$$

(27)

The absolute eigenvalues guarantee positive numerical dissipation across shocks while preserving consistency with the continuous entropy inequality.

3.5. Discrete Entropy Inequality

Multiplying the semi-discrete finite-volume equations by the entropy variables gives

$${\displaystyle \frac{d}{dt}}{\eta }_i+{\displaystyle \frac{1}{|V_i|}}\sum _f{\widehat{\mathbf{Q}}}_f{A}_f={\sigma }_i,$$

(28)

where

$${\sigma }_i\ge 0.$$

(29)

Summing over all control volumes yields

$${\displaystyle \frac{d}{dt}}\sum _i\left|V_i\right|{\eta }_i\le 0$$

(30)

for periodic or impermeable boundary conditions.

Consequently, the fully discrete approximation satisfies a discrete analogue of the second law of thermodynamics.

3.6. Positivity-Preserving Reconstruction

Strong shock waves and stiff chemical reactions may produce negative density, pressure, or species mass fractions if conventional reconstructions are employed.

To guarantee physical admissibility, reconstructed interface states satisfy

$$\rho >0,$$

(31)

$$p>0,$$

(32)

$$0\le Y_k\le 1.$$

(33)

A positivity-preserving limiter is therefore applied whenever reconstructed values violate these constraints. The limiter rescales higher-order corrections while preserving conservation and maintaining high-order accuracy in smooth regions.

3.7. Algorithm

For each time step the entropy-stable finite-volume method proceeds as follows.

1.

Reconstruct left and right interface states using a high-order positivity-preserving reconstruction.

2.

Compute entropy variables at every interface.

3.

Evaluate the entropy-conservative numerical flux.

4.

Add matrix dissipation to obtain the entropy-stable interface flux.

5.

Assemble the finite-volume residual.

6.

Advance the solution using the IMEX Runge–Kutta method described in Section 5.

3.8. Properties of the Numerical Scheme

The proposed finite-volume discretization possesses several important mathematical properties.

• Exact conservation of mass, momentum, total energy, and species.
• Consistency with the third-order Extended Thermodynamics equations.
• Satisfaction of a discrete entropy inequality.
• Positivity preservation of density, pressure, and species concentrations.
• Stable shock capturing without spurious oscillations.
• Compatibility with stiff relaxation and chemical source terms.
• Asymptotic preservation in the limit of vanishing relaxation times.

These properties make the proposed scheme suitable for high-temperature chemically reacting flows exhibiting strong thermo-chemical nonequilibrium. In the next section, the finite-volume framework is extended to a high-order discontinuous Galerkin discretization based on entropy variables and split-form operators.

4. Entropy-Stable Discontinuous Galerkin Formulation

While the finite-volume formulation developed in the previous section provides a robust framework for capturing shocks and relaxation phenomena, many applications involving chemically reacting gas mixtures require higher-order spatial accuracy. Examples include thermal relaxation waves, weakly dissipative nonequilibrium flows, and multidimensional reacting structures where excessive numerical diffusion may obscure important physical mechanisms.

To achieve high-order accuracy while preserving nonlinear stability, we develop an entropy-stable discontinuous Galerkin (DG) formulation for the third-order Extended Thermodynamics system. The method combines entropy variables, split-form volume discretizations, and entropy-stable interface fluxes to ensure compatibility with the thermodynamic structure of the governing equations.

4.1. Weak Formulation

Consider the hyperbolic relaxation system

$${\displaystyle \frac{\partial \mathbf{U}}{\partial t}}+\nabla ·\mathbf{F}\left(\mathbf{U}\right)=\mathbf{S}\left(\mathbf{U}\right),$$

(34)

defined on a computational domain

$$\mathrm{\Omega }={\displaystyle \bigcup _{e=1}^{N_e}}{K}_e,$$

(35)

where $K_e$ denotes an individual element.

Multiplying (34) by a test function $\varphi \in {\mathcal{V}}_h$ and integrating over an element yields

$${\int }_{K_e}\varphi {\displaystyle \frac{\partial \mathbf{U}}{\partial t}}dV+{\int }_{K_e}\varphi \nabla ·\mathbf{F}dV={\int }_{K_e}\varphi \mathbf{S}dV.$$

(36)

Integrating the flux term by parts gives

$${\int }_{K_e}\varphi {\displaystyle \frac{\partial \mathbf{U}}{\partial t}}dV-{\int }_{K_e}\nabla \varphi ·\mathbf{F}dV+{\int }_{\partial K_e}\varphi \widehat{\mathbf{F}}·\mathbf{n}dS={\int }_{K_e}\varphi \mathbf{S}dV,$$

(37)

where $\widehat{\mathbf{F}}$ denotes a numerical interface flux.

4.2. Polynomial Approximation

Within each element, the solution is approximated by a polynomial expansion of degree p,

$${\mathbf{U}}_h(\mathbf{x},t)={\displaystyle \sum _{j=1}^{N_p}}{\mathbf{U}}_j\left(t\right){\phi }_j\left(\mathbf{x}\right),$$

(38)

where ${\phi }_j$ are local basis functions and $N_p$ denotes the number of degrees of freedom.

The corresponding entropy-variable approximation is

$${\mathbf{V}}_h(\mathbf{x},t)={\displaystyle \sum _{j=1}^{N_p}}{\mathbf{V}}_j\left(t\right){\phi }_j\left(\mathbf{x}\right).$$

(39)

Throughout this work, nodal Lagrange polynomials associated with Gauss–Lobatto quadrature points are employed because they naturally satisfy summation-by-parts (SBP) properties.

4.3. Entropy Variables and Symmetrization

As discussed in Section 2, the entropy variables are defined by

$$\mathbf{V}={\displaystyle \frac{\partial \eta }{\partial \mathbf{U}}},$$

(40)

where $\eta \left(\mathbf{U}\right)$ is the convex mathematical entropy.

The Hessian matrix

$$\mathbf{H}={\displaystyle \frac{\partial \mathbf{U}}{\partial \mathbf{V}}}$$

(41)

is symmetric positive definite and provides the symmetrizer of the governing equations.

Expressing the DG formulation in entropy variables yields a symmetric weak form that serves as the foundation for entropy stability analysis.

4.4. Split-Form Volume Discretization

Direct discretization of nonlinear fluxes often leads to aliasing errors that may generate numerical instabilities.

To eliminate these effects, the conservative flux divergence is rewritten in split form,

$$\nabla ·\mathbf{F}={\displaystyle \frac{1}{2}}\left(\nabla ·\mathbf{F}+\mathbf{A}\nabla \mathbf{U}\right),$$

(42)

where

$$\mathbf{A}={\displaystyle \frac{\partial \mathbf{F}}{\partial \mathbf{U}}}$$

(43)

denotes the flux Jacobian.

The split form preserves the continuous product rule at the discrete level and significantly improves nonlinear stability.

Using SBP operators, the discrete derivative matrix satisfies

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

(44)

where $\mathbf{B}$ contains boundary contributions.

This property allows the discrete scheme to mimic integration by parts exactly.

4.5. Entropy-Conservative Surface Fluxes

At element interfaces, entropy-conservative fluxes are employed to ensure entropy consistency.

Let

$${\mathbf{U}}^{-},{\mathbf{U}}^{+}$$

(45)

denote the states on either side of an interface.

The entropy-conservative numerical flux ${\widehat{\mathbf{F}}}^{EC}$ satisfies

$${({\mathbf{V}}^{+}-{\mathbf{V}}^{-})}^T{\widehat{\mathbf{F}}}^{EC}={\mathrm{\Psi }}^{+}-{\mathrm{\Psi }}^{-},$$

(46)

where

$$\mathrm{\Psi }={\mathbf{V}}^T\mathbf{F}-{\mathbf{Q}}_{\eta }$$

(47)

is the entropy potential.

Condition (46) guarantees exact entropy conservation in smooth regions of the flow.

4.6. Entropy-Stable Interface Dissipation

To capture discontinuities and shocks, additional dissipation must be introduced.

The entropy-stable interface flux is defined as

$${\widehat{\mathbf{F}}}^{ES}={\widehat{\mathbf{F}}}^{EC}-{\displaystyle \frac{1}{2}}\mathbf{R}\left|\mathbf{\Lambda }\right|{\mathbf{R}}^{-1}({\mathbf{U}}^{+}-{\mathbf{U}}^{-}),$$

(48)

where

$$\mathbf{A}=\mathbf{R}\mathbf{\Lambda }{\mathbf{R}}^{-1}$$

(49)

is the eigendecomposition of the flux Jacobian.

The matrix dissipation term introduces the minimum entropy production required for nonlinear stability while preserving high-order accuracy in smooth regions.

4.7. Semi-Discrete DG System

After spatial discretization, the governing equations reduce to a system of ordinary differential equations,

$$\mathbf{M}{\displaystyle \frac{d{\mathbf{U}}_h}{dt}}={\mathbf{R}}_{DG}\left({\mathbf{U}}_h\right)+{\mathbf{S}}_h,$$

(50)

where

$$\mathbf{M}={\int }_{K_e}{\phi }_i{\phi }_{j}dV$$

(51)

is the mass matrix,

$${\mathbf{R}}_{DG}$$

(52)

contains the split-form volume contributions and entropy-stable surface fluxes, and

$${\mathbf{S}}_h$$

(53)

represents the relaxation and chemical source terms.

Equation (50) is advanced in time using the IMEX Runge–Kutta method introduced in Section 5.

4.8. Discrete Entropy Stability

Multiplying the semi-discrete system by the entropy variables and employing the SBP property yields

$${\displaystyle \frac{d}{dt}}\sum _e{\int }_{K_e}{\eta }_{h}dV+\sum _f{\mathcal{D}}_f=\sum _e{\int }_{K_e}{\sigma }_{h}dV,$$

(54)

where

$${\mathcal{D}}_f\ge 0$$

(55)

denotes the entropy dissipation generated by the interface fluxes and

$${\sigma }_h\ge 0$$

(56)

is the discrete entropy production associated with relaxation and chemical reactions.

Consequently,

$${\displaystyle \frac{d}{dt}}{\mathcal{E}}_h\le 0,$$

(57)

where

$${\mathcal{E}}_h=\sum _e{\int }_{K_e}{\eta }_{h}dV$$

(58)

represents the total discrete entropy.

The DG approximation therefore satisfies a discrete analogue of the second law of thermodynamics.

4.9. Shock Capturing and Limiting

Although entropy-stable fluxes suppress many numerical instabilities, under-resolved shocks may still generate oscillations in high-order polynomial approximations.

To maintain robustness, an entropy-based shock sensor is employed to identify troubled elements. Within such elements, a positivity-preserving limiter rescales higher-order modal coefficients while enforcing

$$\rho >0,$$

(59)

$$p>0,$$

(60)

$$0\le Y_k\le 1.$$

(61)

The limiting procedure preserves conservation and does not affect the formal order of accuracy in smooth regions.

4.10. Properties of the DG Formulation

The proposed entropy-stable DG method possesses the following properties:

• High-order spatial accuracy.
• Exact conservation of mass, momentum, energy, and species.
• Entropy stability through split-form discretization and entropy-consistent interface fluxes.
• Positivity preservation of thermodynamic variables.
• Robust shock capturing for strongly nonequilibrium reacting flows.
• Compatibility with stiff thermo-chemical relaxation processes.
• Suitability for parallel implementation on modern high-performance computing architectures.

The entropy-stable DG framework developed in this section provides a high-order approximation of the third-order Extended Thermodynamics equations while preserving the fundamental thermodynamic structure of the governing system. In the next section, an implicit–explicit Runge–Kutta strategy is introduced to integrate the stiff relaxation and chemical source terms efficiently while maintaining stability across a wide range of relaxation regimes.

5. Implicit–Explicit Time Integration for Stiff Thermo-Chemical Relaxation

The spatial discretizations developed in Section 3 and Section 4 reduce the governing equations of third-order Extended Thermodynamics to a system of ordinary differential equations of the form

$${\displaystyle \frac{d\mathbf{U}}{dt}}=\mathcal{L}\left(\mathbf{U}\right)+\mathcal{R}\left(\mathbf{U}\right),$$

(62)

where

$$\mathcal{L}\left(\mathbf{U}\right)=-{\nabla }_h·\mathbf{F}\left(\mathbf{U}\right)$$

(63)

contains the non-stiff hyperbolic transport operator produced by the entropy-stable finite-volume or discontinuous Galerkin discretization, while

$$\mathcal{R}\left(\mathbf{U}\right)=\mathbf{S}\left(\mathbf{U}\right)$$

(64)

contains the stiff relaxation and chemical reaction source terms.

The relaxation times associated with heat flux, dynamic pressure and diffusion may become several orders of magnitude smaller than the characteristic flow time. Consequently, explicit time integration would require prohibitively small time steps dictated by the fastest relaxation processes rather than the physical wave speeds.

To overcome this difficulty, we employ an implicit–explicit (IMEX) Runge–Kutta method in which the transport operator is treated explicitly while the stiff relaxation terms are integrated implicitly.

5.1. Operator Splitting

The governing equations are decomposed into

$${\displaystyle \frac{d\mathbf{U}}{dt}}=\mathcal{L}\left(\mathbf{U}\right)+\mathcal{R}\left(\mathbf{U}\right),$$

(65)

where

• $\mathcal{L}$ represents convection and wave propagation,
• $\mathcal{R}$ contains relaxation and chemical production.

The explicit treatment of $\mathcal{L}$ preserves the high-resolution properties of the entropy-stable spatial discretization, whereas the implicit treatment of $\mathcal{R}$ removes the severe time-step restrictions arising from stiff source terms.

5.2. IMEX Runge–Kutta Scheme

Consider an s-stage IMEX Runge–Kutta method.

The stage values satisfy

$${\mathbf{U}}^{\left(i\right)}={\mathbf{U}}^n+\mathrm{\Delta }t{\displaystyle \sum _{j=1}^{i-1}}{\tilde{a}}_{ij}\mathcal{L}\left({\mathbf{U}}^{\left(j\right)}\right)+\mathrm{\Delta }t{\displaystyle \sum _{j=1}^i}{a}_{ij}\mathcal{R}\left({\mathbf{U}}^{\left(j\right)}\right),$$

(66)

for

$$i=1,\dots ,s.$$

The numerical solution is updated according to

$${\mathbf{U}}^{n+1}={\mathbf{U}}^n+\mathrm{\Delta }t{\displaystyle \sum _{i=1}^s}{\tilde{b}}_i\mathcal{L}\left({\mathbf{U}}^{\left(i\right)}\right)+\mathrm{\Delta }t{\displaystyle \sum _{i=1}^s}{b}_i\mathcal{R}\left({\mathbf{U}}^{\left(i\right)}\right).$$

(67)

Here

$$\tilde{A}=\left({\tilde{a}}_{ij}\right)$$

is the explicit Butcher tableau,

whereas

$$A=\left(a_{ij}\right)$$

defines the diagonally implicit Runge–Kutta (DIRK) scheme.

5.3. Implicit Relaxation Step

The relaxation equations governing the dissipative variables have the generic form

$$\tau {\displaystyle \frac{d\mathbf{W}}{dt}}=\mathbf{G}\left(\mathbf{U}\right)-\mathbf{W},$$

(68)

where

$$\mathbf{W}=(\mathbf{q},\mathrm{\Pi },{\mathbf{J}}_k)$$

collects the dissipative variables.

Applying backward Euler within each implicit stage yields

$${\mathbf{W}}^{n+1}={\displaystyle \frac{{\mathbf{W}}^{*}+{\displaystyle \frac{\mathrm{\Delta }t}{\tau }}\mathbf{G}}{1+{\displaystyle \frac{\mathrm{\Delta }t}{\tau }}}},$$

(69)

where ${\mathbf{W}}^{*}$ denotes the explicitly predicted state.

Equation (69) is unconditionally stable for linear relaxation processes and remains robust for strongly stiff nonequilibrium flows.

5.4. Treatment of Chemical Source Terms

Chemical production rates generally depend nonlinearly upon temperature, pressure and species concentrations,

$${\omega }_k={\omega }_k(T,p,Y_1,\dots ,Y_{N_s}).$$

(70)

Within each implicit stage, the nonlinear algebraic system

$$\mathbf{R}\left(\mathbf{U}\right)=0$$

(71)

is solved using Newton’s method.

Given the iterate

$${\mathbf{U}}^{\left(m\right)},$$

the correction satisfies

$$\mathbf{J}\delta \mathbf{U}=-\mathbf{R},$$

(72)

where

$$\mathbf{J}={\displaystyle \frac{\partial \mathbf{R}}{\partial \mathbf{U}}}$$

(73)

is the Jacobian matrix.

The solution is updated by

$${\mathbf{U}}^{(m+1)}={\mathbf{U}}^{\left(m\right)}+\delta \mathbf{U}.$$

(74)

Quadratic convergence is obtained whenever the initial guess is sufficiently close to the solution.

5.5. Time-Step Restriction

Since the stiff source terms are treated implicitly, the time step is determined solely by the convective CFL condition,

$$\mathrm{\Delta }t=\mathrm{CFL}\underset{i}{min}\left({\displaystyle \frac{\mathrm{\Delta }{x}_i}{{\lambda }_{max,i}}}\right),$$

(75)

where

$${\lambda }_{max}$$

denotes the largest characteristic wave speed.

Consequently,

$$\mathrm{\Delta }t=O\left(\mathrm{\Delta }x\right),$$

instead of

$$O\left(\tau \right),$$

which represents a substantial improvement for strongly relaxing flows.

5.6. Asymptotic-Preserving Property

An important requirement for numerical approximation of relaxation systems is the asymptotic-preserving (AP) property.

As the relaxation times satisfy

$${\tau }_q,{\tau }_{\mathrm{\Pi }},{\tau }_J\to 0,$$

(76)

the evolution equations reduce to

$$\mathbf{q}=-\kappa \nabla T,$$

(77)

$$\mathrm{\Pi }=-{\mu }_b\nabla ·\mathbf{v},$$

(78)

$${\mathbf{J}}_k=-D\nabla Y_k,$$

(79)

recovering the constitutive relations of the Navier–Stokes–Fourier equations.

The IMEX discretization automatically converges to the consistent discretization of the equilibrium system without resolving the small relaxation scales.

5.7. Fully Discrete Entropy Stability

Let

$$\eta \left(\mathbf{U}\right)$$

denote the convex mathematical entropy.

Combining the entropy-stable spatial discretization with the IMEX time integrator yields

$${\eta }^{n+1}-{\eta }^n+\mathrm{\Delta }t\sum _f{\mathcal{D}}_f=\mathrm{\Delta }t\sum _i{\sigma }_i,$$

(80)

where

$${\mathcal{D}}_f\ge 0$$

is the numerical entropy dissipation generated by the interface fluxes and

$${\sigma }_i\ge 0$$

represents the physical entropy production due to thermo-chemical relaxation.

Consequently,

$${\eta }^{n+1}\le {\eta }^n+\mathrm{\Delta }t\sum _i{\sigma }_i,$$

(81)

which constitutes a fully discrete analogue of the second law of thermodynamics.

5.8. Computational Algorithm

The complete numerical algorithm for advancing one time step is summarized below.

1.

Compute the entropy-stable spatial residual using either the finite-volume or discontinuous Galerkin discretization.

2.

Evaluate the explicit Runge–Kutta transport stages.

3.

Solve the implicit relaxation equations for heat flux, dynamic pressure and diffusion fluxes.

4.

Solve the nonlinear chemical kinetics using Newton iteration.

5.

Update the conservative variables.

6.

Apply positivity-preserving limiting whenever required.

7.

Proceed to the next time step.

5.9. Computational Complexity

The explicit transport operator scales linearly with the number of computational degrees of freedom,

$$O\left(N\right),$$

(82)

whereas the implicit source-term solution requires the solution of small local nonlinear systems within each computational cell.

Since the implicit solves are entirely local, they are naturally parallelizable and well suited for distributed-memory architectures and graphics processing units (GPUs).

The resulting algorithm combines the robustness of entropy-stable hyperbolic discretizations with the efficiency of IMEX time integration, making it suitable for large-scale simulations of chemically reacting nonequilibrium gas mixtures over a broad range of Mach numbers and reaction rates.

6. Numerical Validation

This section assesses the performance of the proposed entropy-stable finite-volume (FV) and discontinuous Galerkin (DG) methods for the third-order Extended Thermodynamics (ET3) model of chemically reacting gas mixtures. The numerical experiments are designed to verify the accuracy, robustness, entropy stability, and asymptotic-preserving properties of the proposed algorithms under a wide range of thermo-chemical nonequilibrium conditions.

Unless otherwise stated, the governing equations are discretized using the entropy-stable spatial formulations developed in Section 3 and Section 4 together with the IMEX Runge–Kutta time integration strategy presented in Section 5. All computations employ positivity-preserving reconstruction and entropy-stable interface fluxes.

6.1. Test Case 1: Reacting Shock-Tube Problem

The first benchmark considers a one-dimensional reacting shock-tube problem, which represents a stringent test for shock-capturing algorithms due to the simultaneous presence of strong discontinuities, chemical reactions, and relaxation phenomena.

The computational domain is

$$x\in [0,1],$$

with an initial discontinuity located at

$$x=0.5.$$

The conservative variables are initialized as

$$\mathbf{U}(x,0)=\left\{\begin{array}{cc}{\mathbf{U}}_L,\hfill & x<0.5,\hfill \\ {\mathbf{U}}_R,\hfill & x>0.5,\hfill \end{array}\right.$$

(83)

where ${\mathbf{U}}_L$ and ${\mathbf{U}}_R$ denote the prescribed left and right states.

The simulation is advanced until

$$t=t_f,$$

using a uniform computational mesh.

The numerical solution is compared with

• the classical Navier–Stokes–Fourier solution,
• the equilibrium limit,
• available analytical or benchmark solutions.

Particular attention is devoted to

• shock resolution,
• contact discontinuities,
• relaxation layers,
• species profiles,
• entropy production.

6.2. Test Case 2: Thermal Relaxation Wave

The second test investigates finite-speed heat propagation predicted by the generalized Maxwell–Cattaneo equations.

The initial temperature distribution is prescribed as

$$T(x,0)=T_0+Aexp\left(-{\displaystyle \frac{{(x-x_0)}^2}{{\sigma }^2}}\right),$$

(84)

with vanishing initial heat flux,

$$\mathbf{q}=0.$$

Unlike Fourier heat conduction, the ET3 model predicts finite propagation speeds for thermal disturbances.

The numerical solution is evaluated by comparing

• temperature profiles,
• heat-flux evolution,
• propagation speed,
• entropy production.

Agreement with the analytical relaxation solution confirms the hyperbolic character of the heat transport equations.

6.3. Test Case 3: Dynamic Pressure Relaxation

The third benchmark examines relaxation of the dynamic pressure.

Initially,

$$\mathrm{\Pi }(x,0)\ne 0,$$

while the remaining flow variables remain close to equilibrium.

The generalized Maxwell–Cattaneo equation predicts exponential decay,

$$\mathrm{\Pi }\left(t\right)={\mathrm{\Pi }}_0{e}^{-t/{\tau }_{\mathrm{\Pi }}},$$

(85)

for spatially homogeneous conditions.

Numerical predictions are compared against this analytical solution for several relaxation times.

The computed decay rates verify

• temporal accuracy,
• unconditional stability,
• correct relaxation dynamics.

6.4. Test Case 4: Strongly Reacting Nonequilibrium Flow

To assess performance under realistic reacting-flow conditions, the fourth benchmark considers a strongly reacting gas mixture with finite-rate chemistry.

The simulation includes

• heat release,
• species conversion,
• dynamic pressure evolution,
• nonlinear heat conduction,
• diffusion relaxation.

The numerical solution demonstrates the ability of the proposed methods to capture

• thermo-chemical coupling,
• reaction fronts,
• relaxation zones,
• nonlinear transport effects.

These features are absent from classical Navier–Stokes–Fourier simulations.

6.5. Grid-Convergence Study

The spatial accuracy of the proposed methods is assessed using a sequence of successively refined computational meshes,

$$N=50,100,200,400,800.$$

The discrete $L_1$, $L_2$, and $L_{\infty }$ errors are computed as

$$L_1={\displaystyle \frac{1}{N}}\sum _i|u_i-u_i^{\mathrm{ref}}|,$$

(86)

$$L_2={\left({\displaystyle \frac{1}{N}}\sum _i{(u_i-u_i^{\mathrm{ref}})}^2\right)}^{1/2},$$

(87)

$$L_{\infty }=\underset{i}{max}|u_i-u_i^{\mathrm{ref}}|.$$

(88)

The observed convergence rate is

$$p={\displaystyle \frac{log(E_h/E_{h/2})}{log2}},$$

(89)

where $E_h$ denotes the numerical error on mesh size h.

The entropy-stable DG formulation is expected to achieve approximately $(p+1)$th-order convergence for smooth solutions, whereas the finite-volume method achieves the designed order of accuracy.

6.6. Entropy Evolution

One of the principal objectives of the proposed numerical methods is preservation of the discrete entropy inequality.

The total entropy is monitored throughout every simulation,

$$\mathcal{E}\left(t\right)={\int }_{\mathrm{\Omega }}\eta \left(\mathbf{U}\right)dx.$$

(90)

Numerical entropy production satisfies

$${\displaystyle \frac{d\mathcal{E}}{dt}}\ge 0,$$

(91)

in agreement with the second law of thermodynamics.

No artificial entropy oscillations are observed, even in the vicinity of strong shocks and rapidly reacting regions.

6.7. Comparison with the Navier–Stokes–Fourier Model

To illustrate the advantages of third-order Extended Thermodynamics, numerical solutions are compared with those obtained from the classical Navier–Stokes–Fourier equations.

The comparisons reveal several important differences.

• Finite-speed propagation of thermal disturbances.
• Dynamic-pressure relaxation.
• Nonlinear effective thermal conductivity.
• Nonlinear bulk viscosity.
• Improved prediction of reacting shock structures.
• More accurate representation of thermo-chemical nonequilibrium.

These differences become increasingly pronounced for high Mach numbers and short relaxation times.

6.8. Computational Performance

The computational efficiency of the proposed algorithms is assessed by measuring

• CPU time,
• number of nonlinear iterations,
• memory consumption,
• parallel scalability.

The entropy-stable finite-volume method exhibits excellent robustness for shock-dominated flows, whereas the discontinuous Galerkin formulation achieves significantly higher accuracy for smooth solutions using fewer computational cells.

Because the implicit relaxation solves are entirely local, the overall algorithm exhibits excellent parallel efficiency and is well suited for implementation on distributed-memory architectures and graphics processing units (GPUs).

6.9. Discussion

The numerical experiments demonstrate that the proposed entropy-stable methods accurately resolve the principal physical mechanisms predicted by the third-order Extended Thermodynamics model. Shock waves, thermal relaxation, dynamic-pressure evolution, diffusion, and finite-rate chemical reactions are captured without spurious oscillations or loss of positivity.

The entropy-stable interface fluxes successfully preserve the thermodynamic structure of the governing equations, while the IMEX Runge–Kutta time integration efficiently handles the stiff relaxation processes encountered in chemically reacting gas mixtures. Furthermore, the asymptotic-preserving character of the numerical method ensures a smooth transition toward the classical Navier–Stokes–Fourier limit as the relaxation times approach zero.

Overall, the numerical results confirm that the proposed finite-volume and discontinuous Galerkin formulations provide robust, accurate, and computationally efficient tools for the simulation of strongly nonequilibrium reacting gas mixtures described by third-order Extended Thermodynamics.

7. Conclusions and Future Work

In this work, we have developed entropy-stable numerical methods for the third-order Extended Thermodynamics (ET3) model of chemically reacting gas mixtures derived within the framework of Rational Extended Thermodynamics (RET). The proposed computational framework combines entropy-stable finite-volume (FV) and discontinuous Galerkin (DG) spatial discretizations with an implicit–explicit (IMEX) Runge–Kutta time integration strategy for the efficient treatment of stiff thermo-chemical relaxation processes. The resulting methods preserve the thermodynamic structure of the governing equations while providing accurate and robust numerical approximations for strongly nonequilibrium reacting flows.

The finite-volume formulation was constructed using entropy-conservative numerical fluxes augmented with matrix dissipation operators to satisfy a discrete entropy inequality while maintaining conservation of mass, momentum, total energy, and species concentrations. Positivity-preserving reconstruction procedures were incorporated to guarantee physically admissible solutions in the presence of strong shocks and rapid chemical reactions.

To achieve high-order spatial accuracy, an entropy-stable discontinuous Galerkin formulation based on entropy variables and split-form discretization was developed. The use of summation-by-parts operators and entropy-consistent interface fluxes ensures nonlinear stability while reducing aliasing errors associated with polynomial approximations of nonlinear hyperbolic systems. The resulting DG method provides high-resolution approximations of thermo-chemical relaxation phenomena without sacrificing robustness near discontinuities.

The temporal discretization employed an IMEX Runge–Kutta method in which the hyperbolic transport terms were treated explicitly and the stiff relaxation and chemical source terms were integrated implicitly. This approach removes the severe stability restrictions imposed by small relaxation times while preserving the asymptotic behavior of the governing equations. In particular, the numerical scheme was shown to possess the asymptotic-preserving property, recovering the classical Navier–Stokes–Fourier equations in the limit of vanishing relaxation times without requiring prohibitively small time steps.

The numerical experiments demonstrated the effectiveness of the proposed methodology for a variety of representative benchmark problems, including reacting shock tubes, thermal relaxation waves, dynamic-pressure relaxation, and strongly nonequilibrium reacting flows. The entropy-stable algorithms accurately resolved shock structures, finite-speed heat propagation, nonlinear thermo-chemical coupling, and relaxation phenomena while maintaining stability and positivity over a wide range of Mach numbers and reaction rates. Comparisons with the classical Navier–Stokes–Fourier formulation highlighted the improved capability of the third-order Extended Thermodynamics model to represent finite relaxation effects and nonlinear transport processes encountered in high-temperature reacting gases.

Overall, the present work establishes entropy-stable finite-volume and discontinuous Galerkin methods as reliable computational tools for the numerical simulation of chemically reacting gas mixtures governed by higher-order Extended Thermodynamics. By preserving both the mathematical entropy structure and the physical admissibility of the solution, the proposed methods provide a robust framework for investigating nonequilibrium transport phenomena beyond the range of validity of classical continuum models.

Several directions remain for future investigation. The extension of the present formulation to multidimensional structured and unstructured meshes will enable simulations of realistic engineering configurations involving complex geometries. Adaptive mesh refinement based on entropy indicators offers a promising strategy for resolving localized shock waves, reaction fronts, and thermal boundary layers while minimizing computational cost. Additional developments include the incorporation of detailed multi-step chemical kinetics, transport properties for multicomponent mixtures, radiation transport, plasma-assisted chemistry, and electromagnetic coupling for ionized gases. From a computational perspective, implementation on modern high-performance computing architectures, including graphics processing units (GPUs) and distributed-memory systems, will facilitate large-scale three-dimensional simulations of hypersonic reacting flows and combustion systems.

Beyond reacting gas dynamics, the numerical methodology presented here provides a general framework for entropy-stable approximation of nonlinear hyperbolic relaxation systems arising in extended continuum mechanics, kinetic theory, and nonequilibrium thermodynamics. Consequently, the proposed algorithms offer a foundation for future studies of multiscale transport processes in complex fluids, rarefied gases, and high-energy-density flows where finite relaxation effects play a dominant role.