A new particle-like method for high-speed flows with chemical non-equilibrium

The present work is concerned with the numerical simulation of hypersonic blunt body flows with chemical non-equilibrium. New theoretical and numerical formulations for coupling the chemical reaction to the fluid dynamics are presented and validated. The fluid dynamics is defined for a stationary unstructured mesh and the chemical reaction process is defined for “finite quantities” moving through the stationary mesh. The fluid dynamics is modeled by the Euler equations and the chemical reaction rates by the Arrhenius law. Ideal gases are considered. The thermodynamical data are based on JANNAF tables and Burcat’s database. The algorithm proposed by Liou, known as AUSM+, is implemented in a cell-centered based finite volume method and in an unstructured mesh context. Multidimensional limited MUSCL interpolation method is used to perform property reconstructions and to achieve second-order accuracy in space. The minmod limiter is used. The second order accuracy, five stage, Runge-Kutta time-stepping scheme is employed to perform the time march for the fluid dynamics. The numerical code VODE, which is part of the CHEMKIN-II package, is adopted to perform the time integration for the chemical reaction equations. The freestream reacting fluid is composed of H2 and air at the stoichiometric ratio. The emphasis of the present paper is on the description of the new methodology for handling the coupling of chemical and fluid mechanic processes, and its validation by comparison with the standard time-splitting procedure. The configurations considered are the hypersonic flow over a wedge, in which the oblique detonation wave is induced by an oblique shock wave, and the hypersonic flow over a blunt body. Differences between the solutions obtained with each formulation are presented and discussed, including the effects of grid refinement in each case. The primary objective of such comparisons is the validation of the proposed methodology. Moreover, for the hypersonic flow over a blunt body, solutions obtained for two meshes are shown, compared and analyzed. The numerical solutions are also compared with experimental data.


INTRODUCTION
This work is part of a continuing effort that is being carried out at Instituto de Aeronáutica e Espaço of Departamento de Ciência e Tecnologia Aeroespacial (DCTA/IAE) to develop numerical methods to simulate fluid and gasdynamic flows.The objective of this effort is to conceive comprehensive codes and to test and validate them over simple configurations in which experimental, analytical or numerical solutions are available in the literature.The present paper addresses specifically hypersonic flow problems in which chemical non-equilibrium is present.Hence, there are effects arising from the chemical nonequilibrium condition that cannot be neglected.
Hypersonic flow problems including chemical reactions have become important since the end of World War II with the advent of rockets, hypersonic aircrafts and atmospheric reentry vehicles.For reentry vehicles, the gasdynamic problems become all the more challenging since the flight trajectories traverse a wide range of Mach, Knudsen, Reynolds and Damköhler numbers (Sarma, 2000).The full range of molecular physical phenomena of reentry flights encompasses rarefaction, ionization, radiation, relaxation and thermal and compositional nonequilibrium (Sarma, 2000).For hypersonic aircrafts, the combustion in the engine may take place at supersonic speeds, as in Scramjets (Supersonic Combustion Ramjets) and Shramjets (Shock-Induced Combustion Ramjets) (Ahuja, Tiwari and Singh, 1992).DCTA/IAE is developing a recoverable orbital platform for experimentation in microgravity environment, called SARA (acronym in Portuguese for Satélite de Reentrada Atmosférica).
Development of a solver for reactive hypersonic flow applications started in DCTA/IAE at the end of 1990s.Successful validations of the numerical formulations and the physico-chemical data were carried out over simulations for reactive flows composed of a stoichiometric H 2 -Air mixture (Azevedo and Figueira da Silva, 1997;Korzenowski, 1998;Figueira da Silva, Azevedo and Korzenowski, 1999;Pimentel, 2000;and Figueira da Silva, Azevedo and Korzenowski, 2000).Hypersonic flows around wedges were considered.In these configurations, the oblique shock wave triggers the combustion that eventually takes the form of an oblique detonation wave.Comparisons between numerical and experimental results showed a good agreement concerning the overall structure of the flow.Nevertheless, the first attempt to simulate a reactive flow composed of the same mixture around a blunt projectile was not successful (Guzzo, 2003).The configuration simulated was Lehr's experiment with a steady flow around a sphere-cylinder (Lehr, 1972).The experiment shadowgraph showed that normal shock and the reaction front near the stagnation line seem coupled.These two fronts move away from each other downstream, and such detachment becomes more noticeable as the shock wave becomes more oblique.The first simulation obtained was considered unsuccessful since the detachment was not observed in the numerical solution (Guzzo, 2003).The reaction front remained visibly coupled to the shock wave for the overall solution.The differences between the numerical simulations and the experimental data can be attributed either to numerical errors or to improperly modeled chemical kinetics.The physics of these flows are predominantly driven by reaction kinetics and convection phenomena (Wilson and MacCormack, 1990).The uncertainties of diffusion and mixing can be neglected.
Several authors investigated and simulated both Lehr's steady and unsteady cases (Lee and Deiwert, 1989;Yungster, Eberhardt and Bruckner, 1989;Wilson and MacCormack, 1990;Matsuo and Fujiwara, 1991;Sussman and Wilson, 1991;Ahuja, Tiwari and Singh, 1992;Singh, Ahuja and Carpenter, 1992;Ess and Allen, 2005).Lee and Deiwert (1989) and Yungster, Eberhardt and Bruckner (1989) simulated Lehr's experimental data for freestream Mach numbers 4.18,5.11 and 6.46.They used the Euler equations, a 6-species model plus an inert gas, such as argon or nitrogen, and a 8-reaction chemistry model.The flow field for Mach 4.18 and 5.11 was found to be steady in contrast to the experimental evidence that the flow field is, indeed, unsteady.The grids used were not sufficient to resolve the flow field correctly (Ahuja, Tiwari and Singh, 1992).Wilson and MacCormack (1990) used the Euler equations, and a 13-species and 33-reaction chemistry model.They showed the validity of the reaction model and the importance of adequate grid resolution to properly model the flow physics.The calculations were not time accurate and, hence, the unsteady behavior was not captured (Ahuja, Tiwari and Singh, 1992).Sussman and Wilson (1991) also used the Euler equations, and a 13-species and 33-reaction chemistry model to simulate the flow field for Mach number 4.79.They proposed a new formulation based on logarithmic transformation.Such formulation greatly reduced the number of grid points needed to properly resolve the induction zone.They successfully simulated the unsteady case (Ahuja, Tiwari and Singh, 1992).However, frequency was slightly underpredicted.Ahuja, Tiwari and Singh (1992) used the Navier-Stokes equations, and a 7-species and 7-reaction chemistry model to simulate the flow field for freestream Mach numbers 5.11 and 6.46.They successfully simulated the unsteady case and the frequency oscillation was found to be in a good agreement with the experimentally observed frequency (Ahuja, Tiwari and Singh, 1992).Matsuo, Fujii and Fujiwara (1995) used the Euler equations and a 8-species and 19-reaction chemistry model to simulate the flow field for Mach numbers 4.18 and 4.79.A fine grid distribution, such as 161 x 321 points, was used, and the authors have obtained successful solutions with a good agreement concerning the frequency of the oscillations.
There is a major difference between simulating the steady and the unsteady cases.For the flow field with an unstable reaction front, the critical zone that requires accurate solution is the region near the stagnation line, between the normal shock wave and the projectile leading edge.The unsteady pattern of the overall flow field is due to the convection trigged by the coupled oscillatory reaction front and the convective phenomena in this region of the flow.On the other hand, for the steady case, the numerical solution must be accurate in the region of the induction zone just before the oblique reaction front.In this region, as the shock wave becomes more oblique and the temperature increment diminishes, the detonation wave gradually aligns to the local flow.As long as the chemical kinetics is strongly nonlinear with the species mass fraction, the simulation must be able to properly represent the induction zone in a region of flow that progressively tightens to reacted gases.Ahuja and Tiwari (1994) used the Navier-Stokes equations and a 9-species and 18-reaction chemistry model to simulate Lehr's experiment for Mach number 6.46, which is a steady flow.The detachment between the shock wave and the reaction front was not clearly noticeable.Ess and Allen (2005) simulated Lehr's experiments for freestream Mach numbers 4.48 and 6.46.They used the Navier-Stokes equations and a 9-species and 19-reaction chemistry model.They successfully simulated the unsteady flow field for Mach 4.48 and the oscillation frequency compares well with the frequency obtained experimentally (Ess and Allen, 2005).In the steady flow case, for Mach number 6.46, they used two meshes in order to demonstrate the sensitivity of the physical problem to spatial resolution.In the coarse mesh, with half the number of cells in both directions, the shock and the reaction front remained visually attached.In the finer mesh, the detachment is barely observable.Wilson and Sussman (1993) also simulated Lehr's experiment for Mach number 6.46.They obtained good agreement with the experiment concerning the position of the shock wave and of the reaction front.They used a coarse mesh, such as 52 x 52 grid points.Probably, the topology of the mesh may have contributed to reduce the numerical diffusion of the species mass fractions.In the oblique detonation wave region, the diffusion of the mass fraction is minimized if the faces of the cells, that separate unreacted gases to reacted gases, are aligned to the local flow.The mesh influence in the solution may have been significant in the success of the simulation.Wilson and MacCormark (1992) also successfully simulated Lehr's experiment for Mach number 6.46 using an adaptative mesh algorithm.The final mesh topology probably contributed in the reduction of the numerical diffusion and in the success of the simulation as well.
With the objective of better simulating the flowfield in the region near the oblique detonation wave that progressively aligns to the local flow, a new formulation is proposed in the present paper to couple the chemical kinetics to the fluid dynamics.It is a simple Eulerian/Lagrangian hybrid methodology.The fluid dynamics is defined for a stationary unstructured mesh and the chemical kinetics is defined for finite quantities, or particles, moving through the stationary mesh.The objective of adopting Lagrangian particles is to avoid that the chemical reaction process considers average species mass fractions.The particles carry the chemical composition information and the reaction is performed over each particle independently.The interface between Eulerian and Lagrangian models is conducted in a such way that the density, speed and the internal energy values of the particles are defined as being equal to the corresponding values of the control volume that contains the particle, whereas the mass fractions of the control volume are defined as being equal to the weighted averaged mass fractions of the chemical composition of the particles contained in the control volume.
The paper, initially, concentrates in the comparisons of the hybrid methodology and the conventional Eulerian methodology.Solutions obtained with the two formulations are presented, compared and investigated.For the Eulerian methodology, besides the continuity, momentum and energy equations, the species equations are also defined for stationary meshes.The coupling is achieved by Strang's time-step splitting procedure (Strang, 1968).Lehr's experiment for Mach number equal to 6.46 (Lehr, 1972) is simulated with the Eulerian methodology.Three meshes, different in topology and refinement, are employed.Although the formulation is for unstructured grids, only structured quadrilateral meshes are used.Guzzo and Azevedo (2006) have shown that the regularity and the smoothness of the computational meshes have a very strong positive effect on the quality of solutions.The present calculations are compared to experimental data.An investigation specifically concerned with the reactive mixture is carried out to clarify the challenges of simulating, with the Eulerian methodology, the oblique reaction front that progressively aligns to the local flow.Afterwards, a reactive hypersonic flow around a wedge is simulated using both formulations.Two meshes, different in topology and refinement, are employed.The solutions are presented, compared and analyzed.Another investigation, specifically targeting the reactive mixture, is carried out in order to clarify the characteristics of the solutions obtained with the Eulerian methodology corresponding to the position of the reaction front.Finally, for the hybrid methodology, the solutions are shown for Lehr's experiment for Mach number 6.46, performed over two different meshes.Analyses of the solutions obtained and comparisons between the solutions and the experimental data are discussed.

NUMERICAL FORMULATION
The major emphasis of the present paper concerns the discussion of the Eulerian/Lagrangian hybrid methodology.For such methodology, the fluid dynamics is defined for a stationary mesh and the chemical kinetics is calculated for particles that move through the stationary mesh.This section, therefore, initially presents the hybrid methodology, discussing both its Eulerian and Lagrangian components, as well as the coupling of such components.Afterwards, brief comments are presented for the fully Eulerian formulation.

Eulerian method
The fluid dynamics is modeled by the axisymmetric Euler equations applied over stationary unstructured meshes.Four equations are included: continuity, two-component momentum and energy equations.The species equations are not considered.The integral form of the Euler equations are written for the i-th volume as where r is the radial coordinate, z is the axial coordinate, and V i is the volume of the i-th cell.Q is the vector of the conserved variables, E and F are the inviscid flux vector, and H is the axisymmetric source term.The definition of these terms is where ρ is the density, u r and u z are the velocity components in the radial and axial coordinates, ε is the specific total energy and p is the pressure.
The reinterpretation of the formulation for unstructured meshes follows Azevedo and Korzenowski (1998).The 2nd-order spatial accuracy is obtained using MUSCL reconstruction (van Leer, 1979) and the minmod limiter (Hirsch, 1988), in order to extrapolate primitive variables from the centroid to the faces of the cells.The secondorder accuracy, five-stage Runge-Kutta (Mavriplis, 1988;Mavriplis, 1990) scheme is employed to perform the time march.The boundary conditions are implemented with the use of "ghost", or "slave", cells.For solid wall boundaries, the velocity component normal to the wall in the ghost volume has the same magnitude and opposite sign of its value in the adjacent interior volume, whereas the velocity component tangent to the wall is equal to its internal cell counterpart.For the other properties, zero normal gradients are assumed at the wall.For the entrance boundary, the flow variables receive the freestream values.For the subsonic exit boundary, three properties are extrapolated from interior information.For a supersonic exit boundary, all properties are extrapolated from interior information.For the initial condition, freestream conditions are specified throughout the flowfield.Further details on the boundary and initial conditions are described by Guzzo (2006).

Lagrangian method
The chemical kinetics is defined for particles that move with the fluid.The chemical composition of the p-th particle is defined as where Y is the vector of the chemical species mass fractions and Ω is the chemical source vector, which are written as Here, Y k is the mass fraction and w .k is the mass production term of the k-th chemical species.The chemical reaction kinetics is modeled by Arrhenius law.The VODE numerical code (Byrne and Dean, 1993) is adopted to perform the time integration for the chemical reactions.The reaction rate for each species is determined with the use of the CHEMKIN-II code (Kee, Rupley and Miller, 1991).The variation of the particle position is defined as A simple first-order accuracy, one-stage method is employed to update the position of the particles from stage n to stage n+1, For the initial condition, one particle is added in each cell and a weight, wp, is attributed to it, defined as where r i is the radial coordinate of the i-th cell center, r il is the radial coordinate of its l-th node and (N n )i is the number of nodes.For the entrance boundary, one particle is created periodically in the center of each cell that shares a face with the boundary.A fixed and unique time value is employed for the overall boundary cells.The weight of the particles created in the entrance boundary is defined as where r il and r i2 are the radial coordinates of the two nodes of the i-th cell located at the boundary.When a particle crosses the exit boundary, this particle is discarded.When a particle crosses a solid wall or an axisymmetric boundary, this particle is reallocated inside the flow field.
The particle is positioned at the boundary surface, at the nearest point of the original outside position.

Eulerian-Lagrangian interface
The chemical composition of the control volume is defined as being equal to the weighted averaged chemical composition of the particles contained in the volume, i.e., where (N p ) i is the number of particles contained in the i-th control volume, (Y k ) ji is the mass fraction of the k-th chemical species of the j-th particle contained in the i-th volume and wp ji is the weight attributed to that particle.If a control volume contains no particles, the particles contained in the neighbor cells are considered, where (N v ) i is the number of cells that share one face with the i-th volume and (N p ) l is the number of particles contained in the l-th neighbor.If there are no particles in the control volume and in its neighbors, the chemical composition of the control volume is not updated.The density, speed and the specific internal energy of the j-th particle contained in the i-th control volume are defined as being equal to the values of the control volume, i.e., A temperature limiter was employed.In the simulations performed, when the fluid crosses the oblique detonation wave, as the mixture reacts, chemical energy is released and the internal energy of the gas decreases due to the gas expansion.The internal energy is considered to be composed of chemical and thermal energies.If a particle crosses the oblique reaction front without reacting, its temperature drops and it remains unreacted.The temperature limiter is a simple solution used to avoid that particles very close to each other carry very different chemical composition.The interested reader is referred to the work by Guzzo (2006) for further details on the temperature limiter.The temperature of the particle, considering the limiter, is defined as where T ji is the temperature of the particle without the limiter, T i is the temperature of the control volume that contains the particle, and T il is temperature of its l-th neighbor volume.
The time advancement of the solution from the n-th stage to the (n+1)-th stage is achieved by three independent steps: one step for the fluid dynamics, one for the chemical reaction process, and another for the particle positioning.These steps can be written as where L is the operator of the fluid dynamics, Q is the chemistry operator, and f is the particle positioning operator.Here, r ji is the position vector of the particles, which is written as Δt n is the time step of the n-th stage.A time step for each cell is calculated, Δt i , and a single value, Δt n , is used, where Δt n is the minimum value of the Δt i .The time step for each cell is calculated such that the CFL number is kept approximately constant throughout the field (Figueira da Silva, Azevedo and Korzenowski, 1999).

Constitutive equations
The Balakrishnan and Williams (1993) chemical kinetics mechanism is selected to model the H 2 and air reactive mixture.The thermodynamic properties are defined using the polynomial tables of JANNAF (Anon., 1970) and Burcat's database (Burcat, 1984).

Fully Eulerian formulation
The standard approach for handling flows with chemical reactions for finite difference or finite volume methods consists in also solving the species mass fraction equations using the same Eulerian computational grid developed for the fluid dynamics equations.Hence, in the present case, the Euler equations, including continuity, momentum and energy equations are (loosely) coupled to the species equations.In this context, Strang's time-step splitting procedure (Strang, 1968) is used for time advancement.
As before, the VODE numerical code (Byrne and Dean, 1993) is adopted to perform the time integration for the chemical reactions.The reaction rate for each species is determined with the use of the CHEMKIN-II code (Kee, Rupley and Miller, 1991).The thermodynamic properties are also defined using polynomial tables of JANNAF (Anon., 1970) and Burcat's database (Burcat, 1984).The Balakrishnan and Williams (1993) chemical kinetics mechanism is selected to model the H 2 and air reactive mixture.The same numerical method previously described is used for the solution of the fluid dynamics portion of the problem.The interested reader is referred to the work by Pimentel (2000) and by Pimentel et al. (2002) for further details on the fully Eulerian formulation.

Blunt body problem
Lehr's experiment for Mach 6.46 (Lehr, 1972) is simulated with the fully Eulerian formulation.The flow is composed of a stoichiometric H 2 -air mixture.The geometry is a sphere-cylinder of 15 mm diameter.The free stream conditions are T ∞ = 292 K and p ∞ = 320 mmHg.Simulations are carried out over three structured quadrilateral meshes.
The boundaries of the meshes are presented in Fig. 1.
Grid 1 is composed of 20,000 nodes and 19,671 cells, in which the distribution is 80 × 250 grid points in the normal and longitudinal directions, respectively.The nodes are equally spaced in each direction.Grid 2 is composed of 30,000 nodes and 29,631 cells, in which the distribution is 120 × 250 points, respectively, in the wall-normal and longitudinal directions.The nodes are also equally spaced in each direction.Grid 3 is composed of 60,000 nodes and 59,451 cells, in which the distribution is 150 × 400 grid points.Specifically for Grid 3, the nodes are not equally spaced in each direction.The density of nodes is greater in the region of the shock wave and the reaction front.In the direction normal to the body, from the entrance boundary to the wall, the distance between the nodes is linearly reduced, from the first to the 75th node, to 1/5 of the initial distance.From the 75th to the 150th node, the mesh spacing is kept constant.In the longitudinal direction, from the exit to the symmetry boundary, the spacing between the nodes is linearly reduced from the first to the last node to 1/3 of the initial distance.Grids 1 and 2 are similar in refinement, but differ considerably in topology, as one can clearly see in Figure 3 shows the density and the H 2 O mass fraction contours for the solutions obtained with the three meshes.
For the solutions obtained with Grid 1 and Grid 2, the pattern observed in Lehr's experiment (Lehr, 1972) is not seen.It is not clear that the reaction front separates from the shock wave as the two fronts become more oblique.On the other hand, for Grid 3 results, it is possible to notice the beginning of the detachment of the two fronts.As discussed in the next paragraph, however, such detachment is still quite shy of what is observed in the actual experiments.In any event, with the objective of further characterizing the results for the fully Eulerian formulation, Fig. 4 presents the results obtained with Grid 3 in terms of colored contours of density, temperature and molecular oxygen mass fraction, and it also shows the incipient detachment of shock wave and reaction front.Figure 5 presents a summary of the positions of the shock waves and the reaction fronts for the numerical solutions and for the experimental data.For the numerical solutions, in order to indicate the position of the fronts, two temperature contours are included.The differences in position of the oblique fronts for Grid 1 and Grid 2 evidence the influence of the mesh topology in the results.For Grid 1, for instance, the direction of the oblique reaction front is the same as the one of the longitudinal faces of the cells.For Grid 3, despite  the fact that the shock wave and the reaction front are slightly dislocated upstream, and the distance between them is underpredicted, the positions of the oblique fronts compare better to the experimenal data.The effect of the mesh refinement is very similar to the one observed by Ess and Allen (2005).These authors have simulated Lehr's experiment for Mach 6.46 using two meshes.In the coarse mesh, the shock and the reaction front remained visually attached.In the finer grid, the detachment is barely observable, which is essentially the result here obtained.

Chemical reaction analysis 1
In Lehr's experiment for Mach number 6.46, as the fluid crosses the normal shock wave, the temperature increment is such that the induction time is reduced in a way that it is not possible to distinguish the induction zone.As the shock wave becomes more oblique and the temperature increment diminishes, the induction time increases and the detonation wave starts to move away from the shock wave.At a certain point, the temperature increment is not sufficient to ignite the mixture and the detonation wave progressively aligns to the local flow.Some control volumes in this region of the flow are inevitably composed of gases with very different chemical compositions.In order to clarify the effect of using average species mass fraction for the condition of these control volumes, simulations of chemical reactions under constant volume mode are carried out.Three initial chemical compositions are considered.For the first test, the initial mixture is solely composed of unreacted gases.
For the second one, the mixture is composed of 99% of unreacted gases and 1% of reacted gases.For the third test, the mixture is composed of 95% of unreacted gases and 5% of reacted gases.Figure 6 shows the temperature evolution for the three simulations.The initial conditions are a stoichiometric H 2 -air mixture, ρ = 0.002 g/cm 3 and internal energy such that T = 1200 K for the unreacted mixture.The results demonstrate that even a small variation of the initial chemical composition reduces the induction time considerably.These simulations indicate that the use of average species mass fraction may have anticipated the oblique reaction front for the numerical solutions obtained with the fully Eulerian formulation, as presented in Figs. 3, 4 and 5.

Wedge flow
Hypersonic flow around a wedge is simulated with the Eulerian and with the hybrid Eulerian/Lagrangian formulations.The freestream conditions are a stoichiometric H 2 -air mixture, pressure equal to 0.266 atm, temperature equal to 300 K and Mach number equal to 7. The wedge semi-angle is equal to 35°.Simulations were carried out over two structured quadrilateral meshes.The meshes differ mainly in refinement.The boundaries of the meshes are presented in the Fig. 7. Grid 1 is composed of 34,200 nodes and 33,787 cells, in which the distribution is 114 × 300 points in the transversal and longitudinal directions, respectively.The nodes are equally spaced in each direction.Grid 2 is composed of 90,000 nodes and 89,251 cells, in which the distribution is 150 × 600 points.
The nodes are also equally spaced in each direction.
Position coordinates are in centimeters in Fig. 7.  Figure 8 shows the density and the H 2 O mass fraction contours of the solutions obtained.The overall structure of the flow is similar for all simulations.The distance that separates the shock wave from the detonation wave is progressively reduced up to the triple point and, then, the two fronts remain coupled downstream.Moreover, the same structure of the transition region, that exists between the initial oblique shock wave and the resulting oblique detonation wave, is observed for all solutions.The effect of the mesh is more significant for the Eulerian formulation, since the triple point is better defined for the finer mesh.Further comparison of the results can be seen in Fig. 9, which shows density, temperature and molecular oxygen mass fraction contours for the simulations  in the fine grid (Grid 2) and with both formulations.These flow visualization figures indicate that, in the fine grid, the calculations with the different formulations yield quite similar results, which then provides validation for the presently proposed formulation.It must be emphasized that the wedge flow calculations, using the fully Eulerian formulation implemented in the present computational code, have been extensively validated in the work of Pimentel (2000) and Pimentel et al. (2002), among other references.Hence, such results provide confidence in the proposed hybrid formulation.

Chemical reaction analysis 2
A second investigation of the reactive mixture is carried out in order to clarify the effects regarding the use of average species mass fraction in the induction zone.The previous analysis of the chemical reaction addressed the influence in the oblique reaction front that is approximately aligned to the local flow, in which some volumes may be composed of gases with very different chemical composition.The analysis of the present section addresses the reaction front approximately normal to the local flow.In this condition, the mixture of the control volumes in the induction zone is not composed of gases with very different chemical compositions.
Simple one-dimensional simulations of chemical reactions under constant volume conditions are carried out.The Eulerian mesh used is presented in Fig. 11.The density, energy and velocity are kept constant in space and time.The freestream conditions consider unreacted stoichiometric H 2 -air mixture, ρ = 0.002 g/cm 3 and T = 1200 K.The simulations are performed from the instant t 0 = 0 μs to the instant t N = 20 μs.For a number of volumes equal to N, the time that the fluid takes to traverse a volume, Δt vol , may written as For simplicity, two independent steps are considered, the first one for the fluid dynamics and the second one for the chemical reaction, in this order, i.e., where C is the chemistry operator and Y is the vector of the chemical species mass fractions.Y is defined as f Δt is the time step divided by the time that the fluid remains within each cell, i.e.,  The numerical code VODE (Byrne and Dean, 1993) is adopted to time integrate the chemical composition from the instant t to the instant t + Δt.For the initial and the entrance boundary conditions, the freestream properties are defined for all the volumes.Convergence is achieved when the maximum absolute variation of the temperature in a volume from the instant t to t + Δt becomes lower than a certain tolerance.As expected, before combustion starts, the variation of the temperature is insignificant, and convergence is only allowed when the maximum temperature variation starts to decrease.Figure 12 presents the converged solutions for three meshes, the first one containing 20 cells, as shown in Fig. 12(a), the second one containing 100 cells, as shown in Figs.12(b1) and (b2), and the last one containing 1000 cells, as shown in Figs.12(c1) and (c2).Figures 12(b2) and 12(c2) are simply enlargements of their corresponding Figs.12(b1) and 12(c1), respectively.Six different time step values are considered, f Δt = 0.1, 0.3, 0.5, 0.7, 0.9 and 1.0.The simulation performed with f Δt equal to 1.0 corresponds to a simulation without considering averaged species mass fraction, since the integral information contained in a cell is passed to the following cell without mixing.The solutions presented in Fig. 12 have evidenced that the induction time is reduced when averaged species mass fraction is considered.This effect is minimized by the mesh refinement and by increasing the time step.These results and observations explain why the reaction front and the triple point moved downstream with the grid refinement for the Eulerian formulation in the simulation of the flow over a wedge, as shown in the Fig. 10(a).Moreover, such results also suggest that, even for the finer mesh, the distance from the detonation wave to the shock wave may have be underpredicted, as shown in the Fig. 10(c).

Blunt body results using the hybrid formulation
As before, a supersonic flow composed of a stoichiometric H 2 -air mixture over a sphere-cylinder of 15 mm diameter is simulated.The configuration is the same of Lehr's experiment for Mach number equal to 6.46 (Lehr, 1972).The freestream conditions are T ∞ = 292 K and p ∞ = 320 mmHg.Two simulations are carried out over two different structured quadrilateral meshes.The two meshes differ from each other in refinement and topology.The boundaries of the meshes are presented in Fig. 13.The coarse mesh, Grid 1, is composed of 6,000 nodes and 5,831 cells, in which the distribution is 50 × 120 points in the wall-normal and longitudinal directions, respectively.The fine mesh, Grid 2, is composed of 25,000 nodes and 24,651 cells, in which the distribution is 100 × 250 points in the wall-normal and longitudinal directions, respectively.For both meshes, the nodes are equally spaced in each direction.Position coordinates are in centimeters in Fig. 13.It should be pointed out that these grids are considerably coarser than the ones used for the analysis of the blunt body problem using the fully Eulerian formulation.Moreover, Fig. 14 presents a visualization of the fine grid (Grid 2) for the present simulations.Both an overall view of the grid and an enlargement of the mesh near the front symmetry axis are shown in the figure.
observed in the numerical solutions.Near the stagnation line, it is not possible to distinguish the reaction front from the normal shock wave.As the shock wave becomes more oblique downstream, the two fronts separate from each other.One can observe that the contours for the coarse grid, both for density and for water mass fraction, are a bit wavy.This behavior occurs because Grid 1 is really too coarse.The property contours are much better defined in the solution for Grid 2, as one can clearly see in Fig. 15.However, as pointed out, even for the coarse grid results, the shock wave and the detonation front separate from each other as they move downstream.be pointed out that, in the final solution shown in Fig. 16, a total of approximately 218,000 particles are distributed throughout the computational domain.Moreover, the shock wave position is considered as the M = 6 contour and, in the same fashion, the detonation wave position was taken as the contour for which the O 2 mass fraction is equal to 0.2.
Figure 17 shows the position of the shock wave and the reaction front for the two solutions, i.e., for the two grids and using the hybrid formulation, and for the experimental data.
For the numerical solutions, in order to indicate the position of the fronts in the figure, four temperature contours were included.The effects of the mesh refinement and topology were not significant, since the numerical solutions compare well to each other.Comparisons between the numerical solution and the experiment shadowgraph showed good agreement concerning the shock wave and the reaction front positions, although the last was slightly dislocated.Although not completely coincident with the experiment, the solutions obtained indicated that the numerical errors may have been reduced with the use of the current formulation (for comparison purposes, see, for instance, the results in Guzzo, 2003, or the results in Fig. 5).Differences still remaining between numerical simulations and experimental data can be attributed either to the need of further grid refinement or to improperly modeled chemical kinetics.Both issues will be further addressed in the forthcoming discussion.In any event, and again comparing to the results shown, for instance, in Fig. 5, it is evident that a far better agreement with experiment is achieved with the present hybrid formulation.
Figure 18 shows zooms of two regions of the flow solution using Grid 2, that is, near the axis of symmetry, or normal reaction front, and near the oblique reaction front.This visualization evidences the challenges of the present simulations.Within the same flow field, the regime goes from frozen (Damköhler number « 1) to equilibrium (Damköhler number » 1).The Damköhler number is defined as the flow timescale over the chemical reaction timescale, and it is used to relate chemical reaction timescale to other phenomena occurring in a system.Figure 18(a) shows that, when the chemical reaction timescale is much lower than the convective time scale, the gas reacts nearly instantaneously, whereas Fig. 18(b) shows the opposite behavior away from the stagnation region.The streamtraces show how the reaction front aligns to the local flow as the shock and reaction fronts become more oblique.As long as the chemical kinetics is strongly nonlinear with the chemical composition, the use of average species mass fraction over an arbitrary mesh topology to simulate this region of the flow may require a prohibitively extensive grid refinement.The results presented in the paper have attempted to convey the idea that the use of average species mass fraction in the calculation of the chemical reaction processes introduces additional error, which the hybrid formulation here proposed has helped reduce considerably.However, there are still several other sources of errors, which are present and which must be adequately treated for a successful simulation.In order to avoid any confusion, it is important to stress, first of all, that the paper has not addressed any modeling errors.Modeling errors are associated with the ability of the governing equations to adequately represent the phenomena present in the flowfield or, for instance, the ability of the chemical kinetics mechanism (Balakrishnan and Williams, 1993) adopted to correctly model the H 2 and air reactive mixture behavior.The paper has focused on numerical errors.Nevertheless, it is clear that numerical errors can also have very diverse sources.
The paper has addressed effects that come about due to mesh refinement and mesh topology, and the influence of the different treatment of the chemical reaction equations, i.e., Lagrangian versus Eulerian approach.Other aspects of numerical errors have not been addressed.For instance, the effect of the spatial discretization scheme, for the convective terms, is a very important aspect which has not been discussed at all in the present paper.A detailed discussion of numerical errors associated with the spatial discretization schemes can be seen in Azevedo, Figueira da Silva and Strauss (2010).

CONCLUDING REMARKS
A simple Eulerian/Lagrangian hybrid methodology is proposed and its formulation is presented.This methodology couples an Eulerian description of the gasdynamic equations with a particle formulation to represent the chemical kinetics.The objective of adopting a Lagrangian particle formulation for the chemical kinetics is to avoid the use of averaged species mass fractions in the chemical reaction process, which would necessarily occur if the standard, fully Eulerian formulation were used.In this latter approach, the coupling between fluid dynamics and chemical kinetics is achieved by Strang's time-step splitting procedure, in which fractional time steps are taken, alternatively, in the fluid dynamic equations and in the chemical equations, but using a single Eulerian grid for both sets of equations.The authors further emphasize that the present paper has not specifically addressed aspects related to computational costs.The reason for such approach rests on the fact that, as implemented, the two formulations here discussed have comparable computational costs and, hence, this is not an issue in the present case.It is clear, though, that a particle method may have additional advantages for numerical performance enhancements at almost no extra effort from the user, particularly in the multi-core machines that are currently becoming the rule in most computational environments.Such a discussion, however, is beyond the scope of the present work.
Simulations of hypersonic reactive flows over a blunt body and over a wedge are performed using the two different formulations to couple the chemical reaction to the fluid dynamics.For the simulations of hypersonic flows over a blunt body with the Eulerian formulation, the mesh format and refinement significantly influence the final solution.
For the simulations of hypersonic flows over a wedge with the Eulerian formulation, the distance from the reaction front to the shock wave is underpredicted.This effect is reduced with the refinement of the mesh.The effects of the mesh topology and refinement are considerably less significant for the hybrid formulation.The paper also discussed two independent investigations of the behavior of the reactive mixture, which have shown that the use of averaged species mass fractions may reduce the induction time and shorten the induction zone.The second of such investigations has also demonstrated that these effects may be minimized by mesh refinement and by increasing the time step.Finally, Lehr's experiment for freestream Mach number equal to 6.46 is successfully simulated using the presently proposed hybrid formulation using two different meshes.The grids differ in refinement and topology.The solutions compare well with each other, indicating that the numerical errors may have been reduced with the use of the current formulation.The overall structure of the flow shows good agreement with the experimental data.The shock wave position is also in good agreement with the experiment, and the reaction front presented the same inclination but a slight dislocation from the experimental position.The present results seem to indicate that the hybrid formulation is less influenced by grid parameters and, hence, is a more robust approach for the hypersonic problems of interest.

Fig. 1 .
Grid 3 is much finer and comprises only the flow field around the semi-spherical leading edge.Position coordinates are in centimeters in Fig. 1.An actual view of Grid 3 is shown in Fig. 2. In this figure, only one out of every five points is shown.

Figure 2 :
Figure 2: Visualization of Grid 3, with 150 × 400 points in the normal and longitudinal directions, respectively.Only one out of every five points is shown.

Figure 5 :
Figure 5: Position of the shock wave and reaction front.Comparison of the present computations and the experimental data (Lehr 1972).

Figure 6 :
Figure 6: Temperature versus time for a chemical reaction under constant volume mode.Effects due to a small quantity of reacted gas in the initial mixture.

Figure 9 :
Figure 9: Density [kg/m 3 ], temperature [K] and O 2 mass fraction contours for M ∞ = 7, p ∞ = 0.266 atm and T ∞ = 300 K flow of a stoichiometric H 2 -air mixture over a wedge.Solutions in (a), (b) and (c) use the standard Eulerian formulation, whereas results in (d), (e) and (f) use the hybrid formulation.

Figure 10
Figure 10 shows the positions of the shock wave and of the reaction front for the various solutions obtained.Two temperature contours are used to indicate the position of the fronts.Figure 10(a) presents the influence of mesh refinement in the position of the two fronts for the Eulerian formulation, whereas Fig. 10(b) does the same for the hybrid formulation.In a similar fashion, Fig. 10(c) shows the influence of the formulation in the simulations performed for Grid 2. The effect of the mesh is more significant for the Eulerian formulation, since the reaction front as well as the triple point position move downstream with the

Figure 10 :
Figure 10: Position of shock wave and reaction front: a) effect of grid refinement for the Eulerian formulation; b) effect of grid refinement for the hybrid formulation; c) effect of the formulation for Grid 2.

Figure 11 :
Figure 11: One-dimensional mesh used to simulate the time evolution of a reactive mixture under constant volume conditions.

Figure 14 :
Figure 14: Visualization of the fine mesh (Grid 2), with 100 × 250 points in the normal and longitudinal directions, respectively: (a) Complete grid; (b) Grid near the frontal stagnation point.

Figure 15
Figure15shows the density and the H 2 O mass fraction contours of the solutions obtained for the two meshes.The same pattern of Lehr's experiment(Lehr, 1972) is

Figure 16 :
Figure 16: Flow solution obtained with the hybrid formulation for Grid 2: a) density contours, ρ [kg/m 3 ]; b) temperature contours, T [K]; c) contours for O 2 mass fraction; d) position of shock wave and detonation wave.

Figure 17 :
Figure 17: Position of the shock wave and reaction front.Comparison to experimental data(Lehr, 1972).

Figure 18 :
Figure 18: Zoom over Grid 2 solution with streamtraces: a) normal reaction front; b) oblique reaction front.