## 2 Aero-Thermal Modeling of Rotating Detonation Combustors  
### Unsteady Reynolds-Averaged Navier-Stokes Evaluation  
The commercial software CFD++ developed by Metacomp Technologies [25] was used for the numerical simulations and was based on a finite-volume density-based solver. The convective fluxes were calculated with the Harten-Lax-Van Leer-contact approximate Riemann scheme. Interpolation of the spatial variables was performed through a second-order total variation diminishing (TVD) polynomial interpolation with a continuous limiter for stability reasons. We used an implicit time-integration method based on dual-time stepping with fixed global timestep and local Courant-Friedrichs-Lewy (CFL) number variation. Local iteration was accelerated through a multigrid acceleration scheme based on 4 multigrid cycles. The residual convergence criterion for each inner timestep was 0.1. The governing Navier-Stokes equations are expressed as:  
\[\frac{\partial\mathbf{Q}}{\partial\mathbf{t}}+\frac{\partial(\mathrm{F_{1}}\ +\ G_{1})}{\partial\mathbf{x}}+\frac{\partial(\mathrm{F_{2}}\ +\ G_{2})}{\partial\mathbf{y}}+\frac{\partial(\mathrm{F_{3}}\ +\ G_{3})}{\partial\mathbf{z}}=\ \dot{S}, \tag{1}\]  
where \(\mathbf{Q}\) is representing the dependent variables, F the inviscid flux vectors, \(G\) the viscous flux vectors and \(\dot{S}\) the source terms. Below, Equation (2) represents the inviscid, source and viscous flux vectors of Equation (1). \[\mathbf{Q}=\left( \begin{array}{c}e\\ \rho\\ \rho u\\ \rho v\\ \rho v\\ \rho v\\ \rho v\\ \rho v_{H2}\\ \rho\sigma_{H2O}\\ \rho\sigma_{O2}\end{array} \right);\ \ F_{1}=\left( \begin{array}{c}u(e+p)\\ \rho u\\ \rho u^{2}+p\\ \rho vu\\ \rho vu\\ \rho vu\\ \rho vu_{H2}\\ \rho\sigma_{H2O}\\ \rho\sigma_{O2}\end{array} \right);\ \ F_{2}=\left( \begin{array}{c}v(e+p)\\ \rho v\\ \rho vu\\ \rho v^{2}+p\\ \rho vv\\ \rho v\sigma_{H2}\\ \rho v\sigma_{H2O}\\ \rho v\sigma_{O2}\end{array} \right);\ \ F_{3}=\left( \begin{array}{c}w(e+p)\\ \rho w\\ \rho w\\ \rho w^{2}+p\\ \rho w^{2}+p\\ \rho w\sigma_{H2}\\ \rho w\sigma_{H2O}\\ \rho w\sigma_{O2}\end{array} \right);\ \ \dot{S}=\left( \begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\\ \omega_{H2O}\\ \omega_{O2}\end{array} \right),\]  
\[G_{1}=\left( \begin{array}{c}\dot{q}_{x}-u\tau_{xx}-v\tau_{xy}-w\tau_{xz}\\ 0\\ -\tau_{xx}\\ -\tau_{xy}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{O2O}}{\partial\mathbf{x}}\end{array} \right);\ \ G_{2}=\left( \begin{array}{c}\dot{q}_{y}-u\tau_{yz}-v\tau_{xy}-w\tau_{yz} \\ 0\\ -\tau_{xy}\\ -\tau_{yz}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{y}}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{O2O}}{\partial\mathbf{x}}\end{array} \right);\ \ G_{3}=\left( \begin{array}{c}\dot{q}_{z}-u\tau_{xx}-v\tau_{xy}-w\tau_{zz} \\ 0\\ -\tau_{xz}\\ -\tau_{xz}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{H2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{O2O}}{\partial\mathbf{x}}\\ \rho\frac{\partial\sigma_{O2O}}{\partial\mathbf{x}}\end{array} \right) \tag{2}\]  
The first row represents the energy balance, the second row represents the continuity equation, rows 3-5 represent the Navier-Stokes equations, and the final three rows represent the species conservation equations for \(\mathrm{H_{2}}\), \(\mathrm{O_{2}}\) and \(\mathrm{H_{2}O}\). The diffusion coefficient is a scalar and depends on the local temperature, pressure and mixture viscosity and density via the Schmidt number. The diffusion term is written by assuming [PERSON]'s law of binary diffusion, which states that all species diffuse into one another in an equal way. The mixture consists of four species (\(\mathrm{H_{2}}\), \(\mathrm{O_{2}}\), \(\mathrm{N_{2}}\) and \(\mathrm{H_{2}O}\)).