Numerical Assessment of the Convective Heat Transfer in Rotating Detonation Combustors Using a Reduced-Order Model

: The pressure gain across a rotating detonation combustor offers an efﬁciency rise and potential architecture simpliﬁcation of compact gas turbine engines. However, the combustor walls of the rotating detonation combustor are periodically swept by both detonation and oblique shock waves at several kilohertz, disrupting the boundary layer, resulting in a rather complex convective heat transfer between the ﬂuid and the solid walls. A computationally fast procedure is presented to calculate this extraordinary convective heat ﬂux along the detonation combustor. First, a numerical model combining a two-dimensional method of characteristics approach with a monodimensional reaction model is used to compute the combustor ﬂow ﬁeld. Then, an integral boundary layer routine is used to predict the main boundary layer parameters. Finally, an empirical correlation is adopted to predict the convective heat-transfer coefﬁcient to obtain the bulk and local heat ﬂux. The procedure has been applied to a combustor operating with premixed hydrogen–air fuel. The results of this approach compare well with high-ﬁdelity unsteady Reynolds-averaged Navier–Stokes three-dimensional simulations, which included wall reﬁnement in an unrolled combustor. The model demonstrates that total pressure has an important inﬂuence on heat ﬂux within the combustor and is less dependent on the inlet total temperature. For an inlet total pressure of 0.5 MPa and an inlet total temperature of 300 K, a peak time-averaged heat ﬂux of 6 MW/m 2 was identiﬁed at the location of the triple point, followed by a decrease downstream of the oblique shock, to about 4 MW/m 2 . Maximum discrepancy between the reduced-order model and the high-ﬁdelity solver was 16%, but the present reduced-order model required a computational time of only 200 s, that is, about 7000 times faster than the high-ﬁdelity three-dimensional unsteady solver. Therefore, the present tool can be used to optimize the combustor cooling system.


Introduction
Rotating detonation combustors represent a potential efficiency leap for gas turbine power plants [1]. This boost in efficiency is due to the pressure gain obtained through the detonation process, which leads to a more efficient thermal cycle. A rotating detonation combustor is characterized by an annular channel where combustible mixture is injected and ignited by a self-sustained spinning detonation front at supersonic speeds. However, practical implementation is currently constrained by important thermal management issues due to the high-speed flow and high temperature levels reached inside the combustor. The combustor walls are periodically swept by shock waves at several kilohertz, which results in extraordinary heat-flux levels that limit the operability of the combustor.

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: ∂Q ∂t where Q is representing the dependent variables, F the inviscid flux vectors, G the viscous flux vectors and .
S the source terms. Below, Equation (2) represents the inviscid, source and viscous flux vectors of Equation (1).
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 H 2 , O 2 and 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 Fick's law of binary diffusion, which states that all species diffuse into one another in an equal way. The mixture consists of four species (H 2 , O 2 , N 2 and H 2 O). Since all species diffuse with the same diffusivity, mass is conserved and all mass fractions add up to 1 (σH 2 + σO 2 + σH 2 O + σN 2 = 1).
The thermodynamic properties were based on a real gas model with a temperature-dependent coefficient polynomial for each individual species with the coefficients found in Bride et al. [26]. The Sutherland law was used to model the viscosity for each individual species, and the constants are depicted in Table 1.  [11] was used. To solve the turbulent flow, the k-epsilon model was used for the 2D simulations, while for the three-dimensional simulations, the k-ω SST model was chosen for its better accuracy in the near-wall region and in free-shear layer flows. Before examining the combustor's flow field, the solver was verified with the Sod shock tube problem [27] to characterize the dynamic capabilities of the solver. The propagation of a detonation front was analyzed in a one-dimensional domain for a stoichiometric mixture of H 2 -air at 0.1 MPa and 273 K. The detonation was initiated through a narrow rectangle at high pressure and temperature located at the inlet of the tube (x = 0). The length of the tube was 0.1 m and the numerical cell width was 0.5 mm (Figure 1). The time step was 1 × 10 −8 s, similar to the simulations of Towery et al. [28]. The calculated detonation speed was around 2000 ms −1 and with a relative error below 1% compared to the theoretical Chapman-Jouguet velocity (1986 ms −1 [11]). Figure 1a depicts the pressure along the channel at two different time intervals (t = 2.0 × 10 −5 s and t = 3.0 × 10 −5 s, respectively). The pressure jump behind the detonation front and following pressure decrease related to the expansion fan are observed at positions x = 0.05 m and x = 0.07 m. Figure 1b represents iso-Mach contours where the detonation front is identified by a narrow and large Mach gradient, behind which the flow accelerates to M = 1.05. The temperature behind the shock attained 3000 K, comparable to the equilibrium value behind the detonation front. similar to the simulations of Towery et al. [28]. The calculated detonation speed was around 2000 ms −1 and with a relative error below 1% compared to the theoretical Chapman-Jouguet velocity (1986 ms −1 [11]). Figure 1a depicts the pressure along the channel at two different time intervals (t = 2.0 × 10 −5 s and t = 3.0 × 10 −5 s, respectively). The pressure jump behind the detonation front and following pressure decrease related to the expansion fan are observed at positions x = 0.05 m and x = 0.07 m. Figure 1b represents iso-Mach contours where the detonation front is identified by a narrow and large Mach gradient, behind which the flow accelerates to M = 1.05. The temperature behind the shock attained 3000 K, comparable to the equilibrium value behind the detonation front.  Figure 2a shows the static temperature flow-field of an unfolded two-dimensional RDC with a perimeter of 0.3 m and an axial length of 0.15 m. The numerical grid size was 0.5 mm in both x and y directions and the timestep was 1.0 × 10 −8 s. The outlet boundary condition was set to supersonic flow conditions, and at the inlet, a total pressure of 0.5 MPa and a total temperature of 330 K were applied. Periodicity was imposed at both sides of the RDC.
To evaluate the convergence of the numerical results, the total pressure gain (Equation (3)) was computed at the outlet of the combustor for every time step: The mass-flow-averaged total pressure was retrieved by integrating along the perimeter: The converged total pressure gain was 0.72 with fluctuations below 1% after 3.5 ms (Figure 2b). In the unwrapped RDC, the right triangle spanning across the bottom left is the refilling zone, shown in blue because this is the lowest-temperature region in the contour plot. The smallest leg in the right triangle is the detonation wave, which travels from right to left at a speed of 1840 ms −1 (Udetonation, Figure 2a, which was 6% lower than the theoretical Chapman-Jouguet speed predicted by the Zel'dovich-von Neumann Döring (ZND) model). At the upper side of the detonation wave, a triple point occurs, downstream of which an oblique shock and the slip region are formed. The slip region is the zone where new combustion products mix with the old detonation products from the previous detonation cycle.  To evaluate the convergence of the numerical results, the total pressure gain (Equation (3)) was computed at the outlet of the combustor for every time step: The mass-flow-averaged total pressure was retrieved by integrating along the perimeter: P total out, average = θ=2π θ=0 ρ out V axial P total, out dθ / θ=2π θ=0 ρ out V axial dθ.
The converged total pressure gain was 0.72 with fluctuations below 1% after 3.5 ms (Figure 2b). In the unwrapped RDC, the right triangle spanning across the bottom left is the refilling zone, shown in blue because this is the lowest-temperature region in the contour plot. The smallest leg in the right triangle is the detonation wave, which travels from right to left at a speed of 1840 ms −1 (U detonation , Figure 2a, which was 6% lower than the theoretical Chapman-Jouguet speed predicted by the Zel'dovich-von Neumann Döring (ZND) model). At the upper side of the detonation wave, a triple point occurs, downstream of which an oblique shock and the slip region are formed. The slip region is the zone where new combustion products mix with the old detonation products from the previous detonation cycle. Appl Pathlines characterize the properties of a physical particle within the RDC and provide insight for the boundary layer growth and heat flux. The pathlines were retrieved via integration of the properties along the streamlines for two consecutive passages of the shock with the following equation: ⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗⃗ .
The instantaneous streamlines at a certain timestep in the absolute frame of reference are depicted in white, and two pathlines (dashed lines, colored in black and brown) are also depicted in Figure 3a. The pressure trend for those two pathlines along the axial direction is shown in Figure 3b. The detonation front is marked by the steepest pressure jump close to the inlet, while the moving oblique shock crossing can be described by a smaller pressure jump further downstream.

Reduced-Order Model Based on the Method of Characteristics
This section details the reduced model implementation, followed by a validation with a commercial two-dimensional URANS solver and finally a grid sensitivity of the proposed model.
The reduced-order model that we used to retrieve the flow field consists of an iterative procedure to estimate injection velocity that certifies a stable operation of the RDC [20], in a similar way to the work of Lu and Braun [29]. The flow field is solved in the relative frame of reference of the rotating shock, and the numerical procedure relies on a one-dimensional chemical solver [30] to model the detonation front. In our approach, the complete flow field is computed with a twodimensional method of characteristics solver [31].
The initial inputs to the solver are inlet total pressure and inlet total temperature. This input (P0,1 and T0,1) is represented in Figure 4, which shows the triangular injection zone, where the supersonic shock moves at the velocity U, towards the swirl-free fresh reactants with velocity Vinj, a pressure P1 and a temperature T1. According to the velocity triangles, the relative incoming velocity (WCJ), traveling at the Chapman-Jouguet velocity, defines the inclination of the refilling line (β) and of the detonation front. Immediately behind the detonation wave, the relative Mach number is sonic (Mw = Pathlines characterize the properties of a physical particle within the RDC and provide insight for the boundary layer growth and heat flux. The pathlines were retrieved via integration of the properties along the streamlines for two consecutive passages of the shock with the following equation: The instantaneous streamlines at a certain timestep in the absolute frame of reference are depicted in white, and two pathlines (dashed lines, colored in black and brown) are also depicted in Figure 3a. The pressure trend for those two pathlines along the axial direction is shown in Figure 3b. The detonation front is marked by the steepest pressure jump close to the inlet, while the moving oblique shock crossing can be described by a smaller pressure jump further downstream. Pathlines characterize the properties of a physical particle within the RDC and provide insight for the boundary layer growth and heat flux. The pathlines were retrieved via integration of the properties along the streamlines for two consecutive passages of the shock with the following equation: The instantaneous streamlines at a certain timestep in the absolute frame of reference are depicted in white, and two pathlines (dashed lines, colored in black and brown) are also depicted in Figure 3a. The pressure trend for those two pathlines along the axial direction is shown in Figure 3b. The detonation front is marked by the steepest pressure jump close to the inlet, while the moving oblique shock crossing can be described by a smaller pressure jump further downstream.

Reduced-Order Model Based on the Method of Characteristics
This section details the reduced model implementation, followed by a validation with a commercial two-dimensional URANS solver and finally a grid sensitivity of the proposed model.
The reduced-order model that we used to retrieve the flow field consists of an iterative procedure to estimate injection velocity that certifies a stable operation of the RDC [20], in a similar way to the work of Lu and Braun [29]. The flow field is solved in the relative frame of reference of the rotating shock, and the numerical procedure relies on a one-dimensional chemical solver [30] to model the detonation front. In our approach, the complete flow field is computed with a twodimensional method of characteristics solver [31].
The initial inputs to the solver are inlet total pressure and inlet total temperature. This input (P0,1 and T0,1) is represented in Figure 4, which shows the triangular injection zone, where the supersonic shock moves at the velocity U, towards the swirl-free fresh reactants with velocity Vinj, a pressure P1 and a temperature T1. According to the velocity triangles, the relative incoming velocity (WCJ), traveling at the Chapman-Jouguet velocity, defines the inclination of the refilling line (β) and of the

Reduced-Order Model Based on the Method of Characteristics
This section details the reduced model implementation, followed by a validation with a commercial two-dimensional URANS solver and finally a grid sensitivity of the proposed model.
The reduced-order model that we used to retrieve the flow field consists of an iterative procedure to estimate injection velocity that certifies a stable operation of the RDC [20], in a similar way to the work of Lu and Braun [29]. The flow field is solved in the relative frame of reference of the rotating shock, and the numerical procedure relies on a one-dimensional chemical solver [30] to model the detonation front. In our approach, the complete flow field is computed with a two-dimensional method of characteristics solver [31].
The initial inputs to the solver are inlet total pressure and inlet total temperature. This input (P 0,1 and T 0,1 ) is represented in Figure 4, which shows the triangular injection zone, where the supersonic shock moves at the velocity U, towards the swirl-free fresh reactants with velocity V inj , a pressure P 1 and a temperature T 1 . According to the velocity triangles, the relative incoming velocity (W CJ ), traveling at the Chapman-Jouguet velocity, defines the inclination of the refilling line (β) and of the detonation front. Immediately behind the detonation wave, the relative Mach number is sonic (M w = 1) and the detonation jump conditions are computed (P 2 /P 1 , T 2 /T 1 ) based on the ZND theory with the shock and detonation tool box [32]. In the present reduced-order model, both the edge of the refilling zone and the detonation front are considered to be straight lines. The edge of the refilling zone follows the same direction as the streamlines (Wcj), and the detonation wave is perpendicular to the incoming fresh mixture [20]. Therefore, the triangular zone is characterized by a 90 degree angle between the detonation wave and the edge of the refilling zone. 1) and the detonation jump conditions are computed (P2/P1, T2/T1) based on the ZND theory with the shock and detonation tool box [32]. In the present reduced-order model, both the edge of the refilling zone and the detonation front are considered to be straight lines. The edge of the refilling zone follows the same direction as the streamlines (Wcj), and the detonation wave is perpendicular to the incoming fresh mixture [20]. Therefore, the triangular zone is characterized by a 90 degree angle between the detonation wave and the edge of the refilling zone. The model solves an unfolded rotating detonation combustor, illustrated in Figure 5a, in which a single period is represented. To analytically solve the triple point, which intersects the oblique shock wave, slip line and detonation front, continuity of pressure and flow direction (θ ' 3 = θ3) is guaranteed across the slip line, as introduced by Sichel and Foster [33] for planar detonations. However, in rotating detonation combustors, the flow conditions just in front of the oblique shock wave are not known beforehand and need to be updated each time the solver iterates over the injection velocity. The oblique shock equations are used to retrieve the correct oblique shock-wave angle and the slip line, which ensures pressure equality. This then provides the needed boundary conditions to initiate the method of characteristics [20]. At each iteration, the pressure at the inlet wall (Pwall) is evaluated, and convergence is obtained if the pressure at the start of injection of the fresh mixture matches the injection pressure. At the end of the iterative procedure, absolute flow properties are computed for all the points. Figure 5b shows how the reduced model solves the flow field. The blue lines depict the characteristic lines from the method of characteristics, while the green lines show the oblique shock streamlines downstream of the oblique shock (in red). The reduced model needs a CPU time of 200 s on a laptop, while the 2D URANS required 40 h on 4 cores of an Intel Xeon -E5-2660 v3 processor @ 2.6 Ghz (64 GB RAM) on the RICE purdue cluster (built through a partnership with HP and Intel in April 2015 [34]).
To verify that the characteristics net is sufficiently refined, Figure 6a shows a sensitivity analysis of the characteristics net that was solved via the method of characteristics. Three simulations were selected with three different characteristic net sizes (14, 20 and 29 characteristic lines) emerging from the detonation front. Figure 6b-left shows the inlet Mach number for the three cases, and the relative error between fine and medium mesh is within 2%. Figure 6b-right shows the Mach number at the outlet of the domain, and the medium mesh shows good agreement with the fine mesh. The model solves an unfolded rotating detonation combustor, illustrated in Figure 5a, in which a single period is represented. To analytically solve the triple point, which intersects the oblique shock wave, slip line and detonation front, continuity of pressure and flow direction (θ ' 3 = θ 3 ) is guaranteed across the slip line, as introduced by Sichel and Foster [33] for planar detonations. However, in rotating detonation combustors, the flow conditions just in front of the oblique shock wave are not known beforehand and need to be updated each time the solver iterates over the injection velocity. The oblique shock equations are used to retrieve the correct oblique shock-wave angle and the slip line, which ensures pressure equality. This then provides the needed boundary conditions to initiate the method of characteristics [20]. At each iteration, the pressure at the inlet wall (P wall ) is evaluated, and convergence is obtained if the pressure at the start of injection of the fresh mixture matches the injection pressure. At the end of the iterative procedure, absolute flow properties are computed for all the points. Figure 5b shows how the reduced model solves the flow field. The blue lines depict the characteristic lines from the method of characteristics, while the green lines show the oblique shock streamlines downstream of the oblique shock (in red). The reduced model needs a CPU time of 200 s on a laptop, while the 2D URANS required 40 h on 4 cores of an Intel Xeon -E5-2660 v3 processor @ 2.6 Ghz (64 GB RAM) on the RICE purdue cluster (built through a partnership with HP and Intel in April 2015 [34]).
To verify that the characteristics net is sufficiently refined, Figure 6a shows a sensitivity analysis of the characteristics net that was solved via the method of characteristics. Three simulations were selected with three different characteristic net sizes (14, 20 and 29 characteristic lines) emerging from the detonation front. Figure 6b-left shows the inlet Mach number for the three cases, and the relative error between fine and medium mesh is within 2%. Figure 6b-right shows the Mach number at the outlet of the domain, and the medium mesh shows good agreement with the fine mesh.  The accuracy of the reduced-order model was finally assessed through comparisons with the URANS results. Figure 7a depicts the total temperature evaluated with CFD++, while Figure 7b shows the results from the reduced-order model. Figure 7c shows the inlet pressure profile for four successive oblique shock traverses. The reduced-order model is able to accurately predict the onset of the refilling location as well as the speed of injection. In Figure 7d, the outlet total pressure is superimposed for both computations. There is a good agreement in the value of peak pressure, but the profile is distinct due to different expansions behind the slip zone.  The accuracy of the reduced-order model was finally assessed through comparisons with the URANS results. Figure 7a depicts the total temperature evaluated with CFD++, while Figure 7b shows the results from the reduced-order model. Figure 7c shows the inlet pressure profile for four successive oblique shock traverses. The reduced-order model is able to accurately predict the onset of the refilling location as well as the speed of injection. In Figure 7d, the outlet total pressure is superimposed for both computations. There is a good agreement in the value of peak pressure, but the profile is distinct due to different expansions behind the slip zone. The accuracy of the reduced-order model was finally assessed through comparisons with the URANS results. Figure 7a depicts the total temperature evaluated with CFD++, while Figure 7b shows the results from the reduced-order model. Figure 7c shows the inlet pressure profile for four successive oblique shock traverses. The reduced-order model is able to accurately predict the onset of the refilling location as well as the speed of injection. In Figure 7d, the outlet total pressure is superimposed for both computations. There is a good agreement in the value of peak pressure, but the profile is distinct due to different expansions behind the slip zone.

Heat-Transfer Analysis: URANS vs. Reduced-Order Model with Heat-Flux Correlation
The integral method established by Reshotko and Tucker [35] was used to compute the heat transfer between the flow boundary layer and wall, based on the estimation of the boundary layer pressure gradient from the MOC code [20]. The calculation of the aerodynamic properties of the boundary layer involves the momentum integral Equation (6) where is the transformed momentum thickness defined as follows: If a fully developed turbulence flow is considered, Equation (8) can be integrated to predict the value of the momentum thickness at a distance x from the leading edge.
Values with subscript e are defined at the edge of the boundary layer, while a0 and v0 are the speed of sound and kinematic viscosity based on the stagnation flow conditions. , is the incompressible form factor for a flat plate, while Tref is evaluated at the reference enthalpy: The adiabatic enthalpy (ℎ ) is calculated as:

Heat-Transfer Analysis: URANS vs. Reduced-Order Model with Heat-Flux Correlation
The integral method established by Reshotko and Tucker [35] was used to compute the heat transfer between the flow boundary layer and wall, based on the estimation of the boundary layer pressure gradient from the MOC code [20]. The calculation of the aerodynamic properties of the boundary layer involves the momentum integral Equation (6) where θ tr is the transformed momentum thickness defined as follows: If a fully developed turbulence flow is considered, Equation (8) can be integrated to predict the value of the momentum thickness at a distance x from the leading edge. where Values with subscript e are defined at the edge of the boundary layer, while a 0 and v 0 are the speed of sound and kinematic viscosity based on the stagnation flow conditions. H i, f p is the incompressible form factor for a flat plate, while T ref is evaluated at the reference enthalpy: The adiabatic enthalpy (h aw ) is calculated as: Finally, Reshotko and Tucker defined the Stanton number: The heat flux (q wall ) between the fluid and the wall is computed in the following way: A simplified numerical test case was selected to verify the accuracy of the present correlation. Figure 8a depicts the test case that comprises a two-dimensional tube with a channel height of 5 mm with a premixed H 2 -air at stoichiometric conditions, with a high-pressure and high-temperature zone that ignited the combustible at the start of the simulation. The lower and upper walls were isothermal with a grid that was refined at the wall to guarantee that the viscous sublayer was well resolved. The Mach number and total conditions were measured 0.4 mm from the wall and represent the properties at the edge of the boundary layer. The wall and stagnation properties were then loaded into the model to predict the boundary layer properties to solve Equation (14) and retrieve the heat flux. Figure 8b details the heat flux on the lower wall, predicted by the URANS solver and by the correlation of Reshotko and Tucker [35]. Through this test case, the value of H i,fp = 1.7 was set as the optimal value to correctly reproduce the heat flux peak in flows that are facing shock and detonation transients.
The heat flux (qwall) between the fluid and the wall is computed in the following way: A simplified numerical test case was selected to verify the accuracy of the present correlation. Figure 8a depicts the test case that comprises a two-dimensional tube with a channel height of 5 mm with a premixed H2-air at stoichiometric conditions, with a high-pressure and high-temperature zone that ignited the combustible at the start of the simulation. The lower and upper walls were isothermal with a grid that was refined at the wall to guarantee that the viscous sublayer was well resolved. The Mach number and total conditions were measured 0.4 mm from the wall and represent the properties at the edge of the boundary layer. The wall and stagnation properties were then loaded into the model to predict the boundary layer properties to solve Equation (14) and retrieve the heat flux. Figure 8b details the heat flux on the lower wall, predicted by the URANS solver and by the correlation of Reshotko and Tucker [35]. Through this test case, the value of Hi,fp = 1.7 was set as the optimal value to correctly reproduce the heat flux peak in flows that are facing shock and detonation transients. The dependence of the numerical accuracy with the first cell size (in terms of y+), time discretization and total number of cells was retrieved based on the Richardson extrapolation method [36]. The mesh independence study for the current manuscript was performed on a 2D mesh to determine the discretization requirements for the 3D test case. Figure 9a displays the log-log plot of the extrapolated numerical error of the maximum heat-flux level (at the shock front) as a function of the nondimensional wall distance, y+. For values of y+ smaller than 0.5, which relate to a dimensioned wall distance of 2.0 × 10 −8 m, the relative error is less than 10%. On Figure 9b, the extrapolated The dependence of the numerical accuracy with the first cell size (in terms of y+), time discretization and total number of cells was retrieved based on the Richardson extrapolation method [36]. The mesh independence study for the current manuscript was performed on a 2D mesh to determine the discretization requirements for the 3D test case. Figure 9a displays the log-log plot of the extrapolated numerical error of the maximum heat-flux level (at the shock front) as a function of the nondimensional wall distance, y+. For values of y+ smaller than 0.5, which relate to a dimensioned wall distance of 2.0 × 10 −8 m, the relative error is less than 10%. On Figure 9b, the extrapolated numerical error on the maximum heat flux is shown as a function of the timestep per iteration. For values smaller than 1.0 × 10 −8 s, the relative difference with regard to the 0.5 × 10 −9 s case is less than 4%. For these reasons, we opted for a timestep of 1.0 × 10 −8 s and a first wall distance of 2.0 × 10 −8 m.
Furthermore, a grid independence study was performed to determine the grid density based on a grid convergence study [36,37]. Figure 9c  values smaller than 1.0 × 10 −8 s, the relative difference with regard to the 0.5 × 10 −9 s case is less than 4%. For these reasons, we opted for a timestep of 1.0 × 10 −8 s and a first wall distance of 2.0 × 10 −8 m. Furthermore, a grid independence study was performed to determine the grid density based on a grid convergence study [36,37]. Figure 9c plots the extrapolated convergence error for three different grid sizes. The average heat flux on the bottom plate was chosen as the key variable. To obtain a feasible calculation time with an acceptable error, the middle mesh with 60 radial grid points was chosen with an extrapolation error of 10.95%. Finally, we applied the heat-flux methodology to our reduced-order model. First, the streamlines and the corresponding path lines were computed to track the static and the absolute stagnation conditions (Me, Te, a0 and T0) to retrieve the instantaneous heat flux at each spatial location via Equation (14).

Limitations of the Reduced-Order Model
The reduced-detonation-order model attempts to predict the macroscopic flow features inside the combustor, and consequently, several detailed features cannot be captured. In particular, due to the 2D assumption, local curvature is not taken into account, nor are unsteady shock-wave boundary layer interactions. Additionally, for simplicity, both the oblique shock wave and the slip line are assumed to behave as a straight line, which can slightly alter the local flow expansion after the detonation front. From the thermal analysis perspective, we assumed an isothermal wall, and therefore we neglect nonuniformities of the wall temperature. We selected the isothermal boundary conditions following the assumption that the combustor is operating under constant periodic conditions, and that the bulk wall temperature is maintained at constant temperature by the cooling system.

Navier-Stokes Solution vs. Reduced-Order Model
To resolve the boundary layer within the three-dimensional rotating detonation combustor, we modelled the combustor as an unrolled annulus, as shown in Figure 10a. The annulus has a thickness of 5 mm, and inlet stagnation conditions equal to 330 K and 0.5 MPa with a turbulent intensity of 5%. For this 3D simulation, turbulent equations were solved via the URANS k-omega SST model. The lateral walls were kept at a constant temperature of 450 K. The first grid point is 2.0 × 10 −8 m away from the wall, to ensure a y+ below one and avoid the need for wall functions. This resulted in a three-dimensional mesh of 6.9 million cells. In order to ensure a converged result, we plot the mean Finally, we applied the heat-flux methodology to our reduced-order model. First, the streamlines and the corresponding path lines were computed to track the static and the absolute stagnation conditions (Me, Te, a 0 and T 0 ) to retrieve the instantaneous heat flux at each spatial location via Equation (14).

Limitations of the Reduced-Order Model
The reduced-detonation-order model attempts to predict the macroscopic flow features inside the combustor, and consequently, several detailed features cannot be captured. In particular, due to the 2D assumption, local curvature is not taken into account, nor are unsteady shock-wave boundary layer interactions. Additionally, for simplicity, both the oblique shock wave and the slip line are assumed to behave as a straight line, which can slightly alter the local flow expansion after the detonation front. From the thermal analysis perspective, we assumed an isothermal wall, and therefore we neglect nonuniformities of the wall temperature. We selected the isothermal boundary conditions following the assumption that the combustor is operating under constant periodic conditions, and that the bulk wall temperature is maintained at constant temperature by the cooling system.

Navier-Stokes Solution vs. Reduced-Order Model
To resolve the boundary layer within the three-dimensional rotating detonation combustor, we modelled the combustor as an unrolled annulus, as shown in Figure 10a. The annulus has a thickness of 5 mm, and inlet stagnation conditions equal to 330 K and 0.5 MPa with a turbulent intensity of 5%. For this 3D simulation, turbulent equations were solved via the URANS k-omega SST model. The lateral walls were kept at a constant temperature of 450 K. The first grid point is 2.0 × 10 −8 m away from the wall, to ensure a y+ below one and avoid the need for wall functions. This resulted in a three-dimensional mesh of 6.9 million cells. In order to ensure a converged result, we plot the mean heat-flux level at consecutive timesteps in Figure 10b and observe that the heat-flux fluctuation is below 0.4%. Total computational time was around 400 h on 16 cores of an Intel Xeon -E5-2660 v3 processor @ 2.6 Ghz (64 GB RAM). To speed up the computational time, the simulations were initially run on a two-dimensional mesh (described in Section 2.1 until full convergence (3.5 ms, as shown in Figure 2). Afterwards, the solution was interpolated into a 3D mesh and run for another 0.5 ms, which represented three passages of the detonation wave.  Figure 11a illustrates the iso-heat-flux contour (in logarithmic scale) at a certain timestep. At the detonation front, the peak heat flux is around 200 MW/m 2 , followed by a smooth decrease set by the subsequent flow expansion. Since the temperature of the fresh mixture is below the wall temperature, a negative heat flux is measured during refill time of the mixture. The unsteady sweeping of the oblique shock wave augments the total temperature and total pressure, and hence increases the heat in this area. Figure 11b

Parametric Study
The impacts of different operating conditions on heat flux and heat-transfer coefficients were examined with the reduced model. Figure 12a depicts the time-averaged heat flux along the axial direction for several inlet total pressures and inlet temperatures. When only modifying the total pressure, the change in heat flux (Figure 12b) was directly correlated with variations in the convective heat-transfer coefficient. In Figure 12c, we show the impact of changing the inlet total temperature  Figure 11a illustrates the iso-heat-flux contour (in logarithmic scale) at a certain timestep. At the detonation front, the peak heat flux is around 200 MW/m 2 , followed by a smooth decrease set by the subsequent flow expansion. Since the temperature of the fresh mixture is below the wall temperature, a negative heat flux is measured during refill time of the mixture. The unsteady sweeping of the oblique shock wave augments the total temperature and total pressure, and hence increases the heat in this area. Figure 11b Figure 11a illustrates the iso-heat-flux contour (in logarithmic scale) at a certain timestep. At the detonation front, the peak heat flux is around 200 MW/m 2 , followed by a smooth decrease set by the subsequent flow expansion. Since the temperature of the fresh mixture is below the wall temperature, a negative heat flux is measured during refill time of the mixture. The unsteady sweeping of the oblique shock wave augments the total temperature and total pressure, and hence increases the heat in this area. Figure 11b compares the time-averaged heat flux from both models, at different streamwise positions. Our model predicts the same level of averaged heat flux in the region of the detonation front with a maximum value up to 6 MW/m 2 and a decrease towards the outlet, at 0.1 m, and the model under-predicts the 3D URANS by 16%. The bulk heat flux estimated by the URANS solution is 4.8 MW/m 2 and 5.1 MW/m 2 for the reduced-order model. Interestingly, the experimental data reported in [23] shows a similar magnitude of heat flux along the axis of the combustor. The maximum bulk heat-flux variation with the 2D MOC code was about 40% (0.4 MW/m 2 ) between the injection location and the triple point, while our 3D simulations predicted a variation of 33%.

Parametric Study
The impacts of different operating conditions on heat flux and heat-transfer coefficients were examined with the reduced model. Figure 12a depicts the time-averaged heat flux along the axial direction for several inlet total pressures and inlet temperatures. When only modifying the total

Parametric Study
The impacts of different operating conditions on heat flux and heat-transfer coefficients were examined with the reduced model. Figure 12a depicts the time-averaged heat flux along the axial direction for several inlet total pressures and inlet temperatures. When only modifying the total pressure, the change in heat flux (Figure 12b) was directly correlated with variations in the convective heat-transfer coefficient. In Figure 12c, we show the impact of changing the inlet total temperature while keeping the inlet total pressure at 0.5 MPa. The temperature variations appear to play a minor role in the overall heat-flux value (unlike the effects of pressure variations). Furthermore, when the total temperature is increased, the pressure jump across the detonation is smaller. Since this resulted in a lower total pressure with a consequent decrease in heat-transfer coefficient (Figure 12d), the heat-flux level is the lowest for the 500 K case.
Appl. Sci. 2018, 8, x 12 of 15 in a lower total pressure with a consequent decrease in heat-transfer coefficient (Figure 12d), the heatflux level is the lowest for the 500 K case.

Conclusions
This manuscript presents a fast methodology to quantify the convective heat flux in detonation combustors. The numerical approach considers the method of characteristics to compute the flow field and an integral boundary layer method to evaluate the local gas-wall heat exchange. With this tool we can estimate the average and local heat flux within the combustor at any desired operating condition. The parametric study exposed a high correlation between the convective heat-transfer coefficient and the inlet total pressure, and a lower correlation with the total inlet temperature. The model also identified the highest time-averaged heat fluxes at the location of the triple point, as equally observed in experimental data. Although the highest localized heat flux occurs all along the detonation front, this effect is then balanced by a negative heat flux from the cold fresh mixture. As a consequence, the average heat flux at the detonation-injection region is similar to the region downstream of the detonation zone. Further downstream, a slight decrease from 6 MW/m 2 to 4 MW/m 2 was observed, which agreed with existing experimental data. To evaluate the precision of the method, we performed 3D URANS on an unrolled combustor with wall refinement. The overall averaged combustor heat-flux value (5.1 MW/m 2 ) for both methods was within 4%, but the reduced model delivered the results in seconds, several orders of magnitude faster than the three-dimensional unsteady Reynolds-averaged Navier-Stokes solver. Thanks to the limited computational time of our reduced-order model, heat-flux distribution can be estimated in a timely manner. This is essential to optimize the design of cooling schemes that can ensure the integrity of future rotating detonation

Conclusions
This manuscript presents a fast methodology to quantify the convective heat flux in detonation combustors. The numerical approach considers the method of characteristics to compute the flow field and an integral boundary layer method to evaluate the local gas-wall heat exchange. With this tool we can estimate the average and local heat flux within the combustor at any desired operating condition. The parametric study exposed a high correlation between the convective heat-transfer coefficient and the inlet total pressure, and a lower correlation with the total inlet temperature. The model also identified the highest time-averaged heat fluxes at the location of the triple point, as equally observed in experimental data. Although the highest localized heat flux occurs all along the detonation front, this effect is then balanced by a negative heat flux from the cold fresh mixture. As a consequence, the average heat flux at the detonation-injection region is similar to the region downstream of the detonation zone. Further downstream, a slight decrease from 6 MW/m 2 to 4 MW/m 2 was observed, which agreed with existing experimental data. To evaluate the precision of the method, we performed 3D URANS on an unrolled combustor with wall refinement. The overall averaged combustor heat-flux value (5.1 MW/m 2 ) for both methods was within 4%, but the reduced model delivered the results in seconds, several orders of magnitude faster than the three-dimensional unsteady Reynolds-averaged Navier-Stokes solver. Thanks to the limited computational time of our reduced-order model, heat-flux distribution can be estimated in a timely manner. This is essential to optimize the design of cooling schemes that can ensure the integrity of future rotating detonation combustors.