On the development of an unstructured grid solver for inert and reactive high speed flow simulations

On the Development of an Unstructured Grid Solver for Inert and Reactive High Speed flow Simulations

Luis Fernando Figueira da Silva

Laboratoire de Combustion et de Détomique ^_ UPR 9028 du CNRS

86960 Futuroscope France

[e-mail]

João Luiz F. Azevedo

Heidi Korzenowsk

CTA/IAE/ASE-N

12228-900 São José dos Campos, SP Brazil

[e-mail]

Abstract

An unstructured grid Euler solver for reactive compressible flow applications is presented. The method is implemented in a cell centered, finite volume context for unstructured triangular grids. Three different schemes for spatial discretization are implemented and analyzed. Time march is implemented in a time-split fashion with independent integrators for the flow and chemistry equations. The capability implemented is tested for inert flows in a hypersonic inlet and for inert and reactive supersonic flows over a 2-D wedge. The results of the different schemes are compared with each other and with independent calculations using a structured grid code. The strengths and the possible weaknesses of the proposed methods are discussed.

Keywords: Unstructured Grids, Finite Volume Methods, High Speed Flow, Reactive Flow, Shock Waves.

Introduction

The use of unstructured grids has received considerable attention of the Computational Fluid Dynamics (CFD) community in the past several years due to the desire of treating complex flow topologies. As flowfields and configurations of interest get to be more complex, it is fairly well accepted that only triangular grids in 2-D, or tetrahedral grids in 3-D, possess the required degree of flexibility necessary to efficiently discretize the computational domain. In the case of chemically reacting gaseous mixtures, complex flow structures may develop even if the geometry of the domain is simple. As a result of the inherent stiffness of the chemical source terms, the time and length scales present in the flowfields may span over several orders of magnitude, immediately leading to an interest on adaptive grids. In such a context, a definite advantage of unstructured grids is the fact that they allow for a far more natural and efficient way of implementing solution adaptive refinement procedures.

The present work is concerned with developing the capability of simulating two-dimensional supersonic flows with chemical reactions using an unstructured grid methodology. From a fluid dynamics point of view, the flowfields of interest in the present case are simulated using the 2-D Euler equations, which are discretized in a cell centered based finite volume procedure on unstructured triangular meshes. Interface fluxes are calculated by three different algorithms, including a central difference type scheme (Jameson, Schmidt and Turkel, 1981) and a van Leer flux vector splitting scheme (van Leer, 1982).

The van Leer scheme is implemented both as a 1st-order and as a 2nd-order accurate scheme. The 2nd-order scheme uses MUSCL extrapolation (Anderson, Thomas and van Leer, 1986) of primitive variables as opposed to actual flux extrapolation. The central difference scheme is always 2nd order for sufficiently smooth meshes. In this case, the interface fluxes are obtained by straightforward arithmetic averages of the properties on both sides of the interface. Since this approach provides no numerical dissipation terms to control nonlinear instabilities, an appropriate blend of undivided Laplacian and biharmonic operators (Azevedo, 1992) are explicitly added as artificial dissipation terms.

The chemical kinetics scheme used involves 9 species, H₂, O₂, H₂O, OH, O, H, HO₂, H₂O₂ and N₂ and is due to Balakrishnan and Williams (1994). The chemical production rates are calculated using the CHEMKIN II package (Kee, Rupley and Miller, 1991). The use of such a detailed chemistry mechanism is needed in order to obtain a correct prediction of the different time and length scales present in the reactive flowfields of interest, as well as their variation with the flow conditions.

Time march uses two different integrators for the flow and chemical kinetics parts of the formulation in a time-split fashion. The flow integrator is a fully explicit, 2nd-order accurate, five-stage Runge-Kutta time stepping scheme. This scheme allows for the use of variable time stepping and implicit residual smoothing procedures to accelerate convergence to steady state cases. In the present context, these convergence acceleration techniques are only used for the perfect gas computations here presented. The integration of the chemical kinetics is achieved through the use of a specially tailored stiff ODE solver (Byrne and Dean, 1993).

The schemes here discussed are initially applied to the solution of supersonic and hypersonic inlet flows as a preliminary validation phase. The fluid is treated as a perfect gas for these calculations and a 2-D inlet configuration which is representative of some proposed inlet geometries for a typical transatmospheric vehicle is considered (Azevedo, 1992). The entrance conditions are varied from a freestream Mach number M_¥ = 4 up to M_¥ = 16 in order to test the two schemes implemented for a wide range of possible inlet operating conditions. The majority of the inlet simulations are not included in the present paper, since they have already been reported elsewhere (Azevedo, 1992, Azevedo and Mitchell, 1995, and Azevedo and Korzenowski, 1996). Only the inlet solutions with the 2nd-order van Leer scheme are new and are included here.

The problem of primary interest in the present case is the flow of a reactive mixture of hydrogen and air over a 2-D wedge. Both inert and reacting gas flows over this configuration are considered. The onset of combustion by an oblique shock wave stabilized by a wedge has been studied previously by various authors (Li, Kailasanath and Oran, 1994, Figueira da Silva and Deshaies, 1995, and Viguier et al., 1996) using structured mesh solvers. The cited references show that, even though the geometry of the problem considered is extremely simple, the transition that occurs between the oblique shock wave stabilized by the wedge and an oblique detonation wave involves fairly complex flowfields. Comparisons are made between structured and unstructured mesh solver results.

It is important to stress that, although unstructured grid methods have now been used for several years in the aerodynamics community, the same is not true for flows with combustion. Typical exothermic reactive flow simulations found in the literature consider structured grid methods. Therefore, the most relevant contribution of the present paper is in the development of an unstructured grid solver which is able to compute compressible reactive flows. Moreover, since the applications of interest here consider supersonic/hypersonic combustion flows, it is important to have schemes which are able to handle strong shock waves. Hence, the ability to obtain 2nd-order spatial accuracy with an unstructured upwind scheme is also an important contribution. Furthermore, although the present effort does not specifically implement an adaptive grid capability, it sets up the framework within which such a capability can be introduced in the future.

Theoretical Formulation

The 2-D Euler equations for a chemically reacting gas mixture can be written in integral form for a 2-D Cartesian coordinate system (Azevedo, 1992) as

(1)

Here, V represents the area of the control volume and S is its boundary. For a stationary mesh, the vector of conserved quantities Q and the convective flux vectors are given by

(2)

The chemical source vector can be written as

(3)

The nomenclature used here is the standard one, such that r is the density, u and v are Cartesian velocity components, e is the total energy per unit of volume and Y_k, w_k and W_k are the mass fraction, the molar production rate and the molecular weight of species k, respectively. The mass fraction of species K is calculated b: . Equation (1) must be suplemented by the equations of state

(4)

(5)

(6)

In these equations, T is the static temperature, R is the universal gas constant and e_k, and c_rk are the internal energy, the standard-state enthalpy and the specific heat at constant pressure per unit mass of species k, respectively. The molar production rates w_k are given by Arrhenius law, and calculated using the CHEMKIN-II package (Kee, Rupley and Miller, 1991). If one considers a homogeneous mixture of gases with constant specific heat c_p, the pressure r can be computed as

(7)

where g is the ratio of specific heats.

The governing equations are discretized in a cell centered context in which the discrete vector of conserved variables for the i-th cell is defined as

(8)

The Euler equations can, then, be rewritten for each i-th control volume as

(9)

It should be observed that, for a cell centered approach, the control volume used for the integration of the governing equations is formed by each triangular cell itself (Batina, 1991). It should be also emphasized that the previous formulation can be used with essentially any chemical kinetics scheme which is more suitable for the problem at hand. In the present work, all reacting flow calculations have considered stoichiometric hydrogen-air mixtures. The chemical kinetics scheme of Balakrishnan and Williams (1994) was used, and this considers 9 species, H₂, O₂, H₂O, OH, O, H, HO₂, H₂O₂ and N₂. For the reactive flowfields of interest here, it is necessary to use such a detailed chemistry mechanism in order to obtain a correct prediction of the different time and length scales present, as well as their variation with the flow conditions.

Spatial Discretization Algorithms

In a cell centered finite volume context, the standard procedure for flux calculation consists of a loop over the control volumes which adds up the contribution of each edge, or side, to form the flux balance for that particular volume. This is usually called a volume-based data structure, which is the equivalent in the present case of an element-based data structure for the finite element community. Although "natural" and straightforward to implement, this procedure is not the most efficient because fluxes end up being computed twice for each edge of the control volume. For an explicit scheme, this means that the code could theoretically run twice as fast simply by implementing some procedure that would avoid computing the same flux twice.

One of the possibilities for solving this problem is to implement a so-called edge-based (or side-based, for the 2-D case) data structure. In this case, the idea is to index the code computations based on the control volume edges. All three schemes implemented in the present work used an edge-based data structure. The basic description of how the connectivity information is handled in this case was presented in Azevedo and Mitchell (1995) and Azevedo and Korzenowski (1996). The specific additional details necessary in the 2nd-order van Leer scheme will be discussed in the corresponding section.

Centered Scheme

In the centered scheme, the convective operator, C(Q_i), which approximates the surface integral in Eq. (9), is defined as

(10)

In this expression, Q_ik is the arithmetic average of the conserved properties in the cells which share the ik interface, i.e.,

(11)

Moreover, in the previous expression (x_k1, y_k1) and (x_k2, y_k2) are the vertices which define the interface between cells i and k. These points are always ordered in the counterclockwise direction for each i-th control volume. An alternate approach for the calculation of the convective operator is to define

(12)

i.e., the interface flux is computed as the arithmetic average of the cell fluxes themselves. Both implementations have been previously tested (Azevedo, 1992, and Azevedo and Mitchell, 1995) and the results were quite similar. In the present case, the definition implied by Eq. (10) is adopted.

The spatial discretization procedure presented in Eqs.(10) or (12) is equivalent to a central difference scheme. Therefore, artificial dissipation terms must be added in order to control nonlinear instabilities (Jameson, 1987). In the present case, the artificial dissipation operator, D(Qi), is formed as a blend of undivided Laplacian and biharmonic operators (see, for instance, Jameson and Mavriplis, 1986, Mavriplis, 1988, and Mavriplis, 1990). These mimic, in an unstructured mesh context, the concept of using 2nd and 4th difference terms as presented by Jameson, Schmidt and Turkel (1981) and Pulliam (1986). Therefore, the artificial dissipation operator is written as

(13)

where d⁽²⁾(Qi) represents the contribution of the undivided Laplacian operator, and d⁽⁴⁾(Qi) the contribution of the biharmonic operator. The process of constructing the biharmonic operator starts by forming the four-point, undivided Laplacian operator at the i-th cell as

(14)

where the summation in k is taken over all the control volumes which have a common interface with the i-th cell. The biharmonic operator is, then, formed as

(15)

The undivided Laplacian artificial dissipation operator is constructed as

(16)

Here, the coefficient is defined as

(17)

where n_i is an scaling given by a normalized undivided Laplacian of the pressure, i.e.},

(18)

In order to switch off the biharmonic operator in the vicinity of shock waves, the coefficient is defined as

(19)

Based on information available in the literature (Jameson, Schmidt and Turkel, 1981, and Mavriplis, 1990), suggested ranges for the K⁽²⁾ and K⁽⁴⁾ constant coefficients are 1/4 £ K⁽²⁾£ 1/2 and 1/256 £ K⁽⁴⁾£ 3/256. Numerical experimentation has indicated that a good compromise is to use K⁽²⁾ = 1/4 and K⁽⁴⁾ = 3/256 for the present applications. The artificial dissipation scheme used in the present work is basically equivalent to the ones described in Jameson and Mavriplis (1986), Mavriplis (1988) and Mavriplis (1990), but with a different scaling of the various terms.

First-Order Upwind Scheme

The upwind version of the algorithm implemented in the present work computes interface fluxes using the van Leer flux-vector splitting scheme (van Leer, 1982). The convective operator, C(Qi), can be written in this case as

(20)

where Dx_ik = x_k2 - x_k1 and Dy_ik = y_k2 - y_k1. In the 1st-order case, the interface fluxes, E_ik and F_ik, are simply defined as

(21)

(22)

Here, and are the split fluxes calculated using van Leer's formulas and the conserved properties of the i-th control volume. The interested reader is referred to the work of van Leer (1982), Anderson, Thomas and van Leer (1986), and Hirsch (1990) for the detailed expressions for the evaluation of the split fluxes in the van Leer context.

As previously discussed, this scheme only gives 1st-order accuracy in space. As a result of that, unless the computational meshes are excessively fine, flow gradients are captured in a rather smeared fashion by the scheme. For the problems at hand, this can be a very serious drawback because it could completely destroy the accuracy with which the chemical process is reproduced by the numerical solution. Therefore, accuracy higher than 1st order seems to be a necessity for the present case. It is important to emphasize that the code has been tested with the 1st-order scheme and it is very robust in this case. However, the desired resolution is simply not there.

Second-Order Upwind Scheme

A 2nd-order van Leer flux vector splitting scheme was implemented in the present case using MUSCL extrapolation (Anderson, Thomas and van Leer, 1986) of the primitive variables (r, u, u, T, Y_k). The convective operator, C(Qi), in this case, can still be written as indicated in Eq. (20), where

(23)

(24)

Here, Q^- and Q⁺ are the "left" and "right"states at the ik interface, respectively. These states are determined using linear extrapolation of primitive variables, as previously indicated, together with some appropriate limiter (Hirsch, 1990).

There are two aspects of the unstructured grid implementation of such a scheme which deserve further consideration. The first aspect concerns the definition of "left" and "right" of a given cell interface. Since there is no concept similar to curvilinear coordinates in this case, the cell interfaces can have virtually any orientation and one must decide which way to "look" in order to construct left and right states. The other aspect is associated with deciding which second control volume will be used for the reconstruction process besides the volume immediately adjacent to the interface considered. The authors again emphasize that an edge-based data structure is being used in this development.

The procedure adopted in the present case to handle the second aspect previously mentioned is inspired in the work of Lyra (1994). The major difference between the present implementation and the cited reference lies in the direction in which the one-dimensional stencil is constructed. In Lyra (1994), the stencil for extrapolation is constructed along the direction of the edge. Here, since a cell centeredfinite volume method is of interest, the extrapolation stencil is constructed in a direction normal to the edge. In an attempt to reinterpret the 1-D ideas in the present unstructured grid context, a line is drawn normal to the edge and passing through the center of the inscribed circle. A third point is located over this line at a distance from the center of the inscribed circle equal to the diameter of the circle. The code, then, identifies in which control volume this 3rd point lies, and it uses the properties of this triangle in the linear reconstruction of the primitive variables.

In order to make the nomenclature clear, the two triangles which are adjacent to the edge under consideration are denoted i and k. The second triangle identified by the previously described process and associated with triangle i is denoted . The corresponding one associated with k is denoted triangle m. This is illustrated in Fig. 1. Therefore, in the calculation of the E^±fluxes, the left state, Q^-, is defined using the properties of the i and triangles and the right state, Q⁺, is defined using those of the k and m triangles, if Dy_ik ³ 0. The reverse is true if Dy_ik < 0. Similarly, the definition of the F^± fluxes uses data at the i and triangles to define the left state and data at the k and m triangles to define the right state if Dx_ik £ 0, and vice-versa if Dx_ik > 0.

Fig. 1

In order to avoid oscillations, the extrapolated states must be limited. In the present case, a simple minmod limiter (Hirsch, 1990) was initially adopted. Computations were also run with the van Leer, van Albada and superbee limiters. However, the convergence of the solution, measured both in terms of a L_¥ and a L₂ norms of the residue, seems to stop after two to four orders of magnitude decay with all but the minmod limiter, for which machine zero convergence is achieved.

Time Discretization Methods

The Euler equations, fully discretized in space and assuming a stationary mesh, can be written as

(25)

where the D(Qi) operator is identically zero if the upwind spatial discretization is used.

Time advancement of the solution from time step n to n+1 is achieved by the use of Strang's time-step splitting procedure (Strang, 1968)

(26)

that separately integrates the fluid dynamics operator L and the chemistry operator C at each cell. This procedure is 2nd-order accurate and gives the flexibility of choosing specialized integrators for the chemical kinetics and the fluid dynamics. A more detailed discussion on this subject can be found in the work of LeVeque and Yee (1990).

Integration of Flow Equations

The present work uses a well-tested, fully explicit, 2nd-order accurate, 5-stage Runge-Kutta time-stepping scheme (Mavriplis, 1988) to advance the "fluid dynamics" part of the governing equations in time. The time integration scheme can, therefore, be written as

(27)

where the superscripts n and n+1 indicate that these are property values at the beginning and at the end of the n-th time step. The values used for the a_l coefficients were

(28)

It should be observed that the convective operator, C(Q), is evaluated at every stage of the integration process, but the artificial dissipation operator, D(Q), is only evaluated at the two initial stages (and, obviously, only for the central difference scheme).

For steady state, inert gas problems, a local time stepping option has been implemented. The objective in this case is to keep an approximately constant CFL number throughout the whole field. In the present work, the time step used to update the conserved variables in the i-th cell, i.e., Dt_i, is calculated as

(29)

where (Ds)_i is the characteristic length associated to the i-th cell, a is the local speed of sound, and

(30)

is the magnitude of the local velocity vector. The idea here is to use the magnitude of the maximum possible eigenvalue of the Jacobian matrices associated with the Euler equations as the characteristic speed with which to compute the CFL number. The least well defined quantity in Eq. (29) is certainly the characteristic length for a triangular mesh. In the present approach, (Ds)_i has been defined as the diameter of the inscribed circle for that triangle. In the reactive flow case, one cannot use a space-varying time step. However, some convergence acceleration can be achieved by recalculating a global time step at each iteration as the minimum of the local Dt_i, defined by Eq. (29).

Further steady state convergence acceleration can be obtained by performing implicit residual smoothing. The procedure implemented follows very closely the work described in Jameson and Baker (1983), and Jameson and Baker (1987). This procedure is only correct for inert flow cases and it has been used in the present study only for the ideal gas simulations. Essentially, the residual at the i-th cell is defined as

(31)

and the averaged residuals, , are calculated by the solution of

(32)

Here, Ѳ ( ) is the undivided Laplacian operator. For a cell centered scheme applied on a triangular mesh, this operator is defined as

(33)

where the k index is looped over all the neighbors of the i-th cell. Equation (32) is solved by a few Jacobi iterations, as usually recommended in the literature (Jameson and Mavriplis, 1986, and Batina, 1989). Typically, only two iterations are used, since there is actually no need to solve this equation accurately. All one is seeking is to obtain some averaging of the residuals at the i-th control volume. Moreover, it is sufficient to apply the smoothing at alternate stages of the Runge-Kutta time-stepping scheme (Jameson and Baker, 1983). Since the total number of stages is odd in the present case, smoothing should be applied at stages 1, 3 and 5.

Chemistry Integrator

As previously stated, the use of a time-step splitting procedure allows the adoption of a specialized solver for the integration of the chemistry operator C in Eq. (26). This corresponds to a separate integration of the ODE

(34)

It can be noticed from Eqs. (34) and (1) - (3) that such a time-step splitting procedure results in a constant density thermal explosion problem at each computational cell.

The authors have chosen to perform the integration of Eq. (34) using VODE (Byrne and Dean, 1993), which is an ODE solver tailored for the solution of problems which include stiff source terms. VODE uses a variable time-step, variable order, backward differentiation formula, together with a modified Newton method whose Jacobian matrix is evaluated numerically. This last feature is of particular value when using different chemical kinetics schemes. A further advantage of using this stiff ODE solver is the fact that it has no stability limits on the choice of the time step Dt. As a consequence, only the fluid dynamics requirements constrain the choice of Dt.

Boundary Condition Implementation

All boundary conditions were implemented with the use of "ghost", or "slave", cells. These are fictitious control volumes which are defined outside the computational domain of interest and which serve the solely purpose of implementing boundary conditions. Solid wall boundary conditions consider the flow tangent to the wall in the inviscid case. This is enforced imposing that the velocity component normal to the wall in the ghost volume has the same magnitude and opposite sign of the normal velocity component in its adjacent interior volume, whereas the ghost volume velocity component tangent to the wall is exactly equal to its internal cell counterpart. Moreover, zero normal pressure, temperature and species mass fraction gradients are assumed at the wall.

Entrance and exit conditions use the concept of one-dimensional characteristic relations in order to determine the conditions that should be specified or extrapolated at those boundaries. For a supersonic entrance, all flow variables are considered known, and they are typically fixed to their freestram values. The exit static pressure is assumed known for a subsonic exit. Three properties are extrapolated from interior information in this case. For a supersonic exit, all properties are extrapolated from interior information.

Results and Discussion

This section will initially discuss the results of hypersonic flow calculations within a 2-D inlet obtained with the use of different spatial discretization strategies.These simulations were performed in the framework of an inert gas mixture with constant specific heats. Then, the effects of the inclusion of a thermal equation of state which takes into account a multi-species mixture with variable specific heats are analyzed for the case of the supersonic flow over a 2-D wedge. Lastly, the results of the transition from an oblique shock wave to an oblique detonation wave obtained with the present code are compared with those calculated using a structured grid program.

Ideal Gas Results

A 2-D inlet configuration which is representative of some proposed inlet geometries for a typical transatmospheric vehicle was used as one of the test cases for the present work. A typical mesh used for the inlet simulations is shown in Fig. 2. This particular mesh has 4204 nodes and 8006 triangles. Results for coarser meshes are also available, but will not be shown here. Numerical tests have indicated that the grid shown in Fig. 2 yields a good compromise between accuracy and convergence rate. It must be emphasized that, for the present inlet simulations, the fluid was treated as a perfect gas with constant specific heat and no chemistry was taken into account. From a physical standpoint, the present simulations are typical of cold gas flows which are usually achieved in experimental facilities such as gun tunnels. This is certainly not representative of actual flight conditions in which real gas phenomena are important for such devices, especially for the higher Mach number cases. However, the purpose of these simulations is to aid in the construction of a robust code for the combustion applications of major interest here.

Fig. 2

For comparison purposes, Figs. 3 and 4 present Mach number contours for the freestream Mach number case M_¥= 12 for the solutions with the centered scheme and with the 1st-order van Leer flux vector splitting scheme. It is clear from Fig. 3 that the flow features are well captured by the numerical solution, except that some oscillations can be seen upstream of the strong upper wall entrance shock. The oscillations are evidenced by the intermediate grey line just upstream and somewhat parallel to the upper wall entrance shock. For the cold gas simulations here performed, these oscillations do not cause too much inconvenience because one can always filter out the obvious numerical oscillations in a post-processing stage. On the other hand, for simulations with chemistry, these oscillations could completely destroy the correct capture of the physical phenomena that are found to occur in practice.

Fig. 3

Fig. 4

The results in Fig. 4 do not exhibit any oscillations. However, they present another problem which is typical of 1st-order schemes. Spatial gradients are smoothed out by the added dissipation of the 1st-order upwind scheme. This means that, in order to recover spatial resolution, one would have to resort to a much finer computational grid. In the present case, the added smearing of the various discontinuities was deemed not acceptable for the applications at hand, which has motivated the implementation of the 2nd-order upwind scheme.

Mach number contours for the solution with the 2nd-order van Leer flux vector splitting scheme are shown in Fig. 5. The minmod limiter was used for this calculation. The corresponding density contours are shown in Fig. 6. It is clear from these figures the much better resolution of discontinuities provided by the 2nd-order scheme, as compared to the 1st-order one. The shock waves present in the flow are very sharply defined in the solution shown in Figs. 5 and 6. Moreover, these results also indicate a solution free of oscillations, with a marked improvement over the calculations with the centered scheme.

Fig. 5

Fig. 6

Inert Gas Mixtures

Before presenting the results in which chemistry was taken into account, it is useful to examine the calculations performed with the different spatial discretization algorithms for an inert gas mixture. In the present case, the flow of a M_¥= 8.0, r_¥ = 0.1 MPa and T_¥ = 300 K, inert stoichiometric mixture of hydrogen and air around a wedge with a half-angle of 25 deg. is considered. The wall boundary condition assumed a wedge with an adiabatic, non-catalytic slip surface. The mesh used contains 1504 triangles and 806 nodes.

The results in Fig. 7 show the converged temperature fields obtained for the different spatial discretization algorithms. It must be noticed that the good prediction of the temperature field is crucial to the onset of the combustion behind an oblique shock wave (Figueira da Silva and Deshaies, 1995, and Viguier et al., 1996). In this sense, the most relevant feature that can be observed in Fig. 7 is the slight discrepancy that exists for the value of the temperature behind the oblique shock wave, which obviously should not depend on the choice of the numerical method. It can be seen that the 2nd-order schemes exhibit post-shock temperatures that are 50-70 K higher than those of the 1st-order upwind scheme. Although this represents only 2% of the post-shock temperature, it will be shown in the following section that such a difference leads to modifications of the combustion induction length value behind the oblique shock wave. A possible explanation for this temperature difference can be the larger entropy increase across the shock wave for the more dissipative 1st-order scheme. It can also be seen in Fig. 7 that the highest wall temperatures are obtained in the case of the more compressive limiters, especially the superbee limiter. This again points in the direction of more dissipative schemes yielding smaller post-shock temperatures. Moreover, one can observe that the numerical shock wave thickness is slightly larger for the 1st-order scheme, as could be expected due to the more dissipative nature of this scheme.

Fig. 7

The L_¥ and the L₂ norms of the residue for the cases presented in Fig. 7 are plotted in Fig. 8. The fastest convergence is obtained for the 2nd-order upwind scheme with the minmod limiter. However, the 1st-order scheme seems to reduce the residue even further than the 2nd-order upwind scheme with the minmod limiter. The other 2nd-order upwind discretizations considered, which use the superbee, the van Leer and the van Albada limiters, do not converge to machine zero at all. The 2nd-order centered scheme has a slower convergence rate, but it does reach machine zero. These results are consistent with those of Venkatakrishnan (1995), who analyzed the role of different limiters on the convergence properties of 2nd-order upwind methods on unstructured meshes. Furthermore, it must be stressed that the temperature fields presented in Fig. 7 were fully invariant with further iteration.

Fig. 8

Based on the previous comparisons, the authors have chosen to retain both the 1st-order upwind scheme and the 2nd-order upwind scheme with a minmod limiter in order to illustrate the application of unstructured meshes to compressible reactive flows.

Results with Combustion

Results for a reactive flow case, considering a stoichiometric hydrogen-air mixture, are shown in Fig. 9. This figure presents the temperature fields obtained with the 1st-order upwind scheme, the 2nd-order upwind scheme using a minmod limiter and the baseline structured mesh code (Figueira da Silva and Deshaies, 1995). As previously mentioned, the present simulations consider a chemical kinetics scheme with 9 species, H₂, O₂, H₂O, OH, O, H, HO₂, H₂O₂ and N₂. Moreover, the particular kinetics scheme used is due to Balakrishnan and Williams (1994), and the chemical production rates were calculated using the CHEMKIN II package (Kee, Rupley and Miller, 1991).

Fig. 9

The unstructured mesh contains 6401 triangles, while the structured mesh contains 100 x 35 points uniformly distributed in both directions. This amounts to about the same overall spatial resolution on both meshes. It is clear that the resolution obtained with these meshes is not sufficient to describe the internal structure of the detonation wave, but the induction length and wave angles obtained with the structured mesh solver are not changed with further mesh refinement.

The structured mesh code is based on the AUSM+ flux vector splitting method (Liou, 1996). It is 2nd-order accurate both in space and time, and it has been successfully validated against experimental results for the oblique shock wave-oblique detonation wave transition and a simple shock/detonation polar analysis (Viguier et al., 1996). Therefore, the temperature field calculated with this code is presumed to be correct, even though the well known AUSM+ wall boundary condition problem is clearly evident in Fig. 9.

It is possible to observe from Fig. 9 that the 1st-order upwind (unstructured grid) solution agrees fairly well with the structured mesh solution, even if the former is highly smeared by the larger numerical diffusion. This smearing is responsible for a very bad resolution of the triple point - intersection between the shock and detonation waves (Figueira da Silva and Deshaies, 1995). The temperature field in the burned gases region also differs somewhat from the 2nd-order structured solution. The induction lengths of the chemical kinetics behind the oblique shock wave are virtually identical, as are the shock and detonation wave angles. It should be noticed that the AUSM+ method is supposed to have better shock capturing capabilities than the van Leer scheme (Liou, 1996). Current efforts, as a continuation of the present work, are directed towards implementation and testing of an unstructured version of the AUSM+ scheme.

The chemical kinetics induction length which was obtained with the use of the 2nd-order unstructured method is smaller than those calculated with either the structured or the 1st-order upwind unstructured codes. This behavior was, actually, to be expected, if one considers the results of the previous section. On the other hand, the 2nd-order solution provides a better representation of the flow in the combustion products region. This can be confirmed by observing the lower right-hand corner of the corresponding temperature fields shown in Fig. 9. It is expected that further resolution of the physical aspects present in the oblique shock wave-oblique detonation wave transition problem will be possible with the use of adaptive mesh refinement techniques within the current unstructured grid context.

Concluding

The present paper has described in detail an unstructured mesh solver which was developed for the calculation of 2-D inviscid compressible reactive flows. Time advancement of the solution is achieved via a time step splitting procedure which consecutively integrates the fluid dynamics and the chemistry terms of the conservation equations. The calculation of the thermodynamic properties and chemical reaction rates is performed by the CHEMKIN-II package. A few space discretization schemes are included in the code. These consider both 1st- and 2nd-order spatial discretization algorithms. Furthermore, both centered and upwind discretizations are contemplated.

This unstructured grid code was tested both for inert and chemically reacting flow situations. The former included a representative air intake at hypersonic flow conditions and a supersonic wedge. For the hypersonic inlet case, all schemes were able to produce converged results, despite the very high Mach number conditions considered. However, the centered scheme presents oscillations in the converged solutions, especially for the higher Mach number cases. The 1st-order upwind scheme was able to produce solutions free of oscillations, at the expense of smearing the flow discontinuities. The 2nd-order upwind scheme has shown oscillation free solutions, with a marked improvement over the calculations with the other schemes for this case.

For the supersonic wedge case, both the 1st-order upwind and the 2nd-order centered schemes give a good prediction of the shock wave intensity. On the other hand, results obtained for this case with the 2nd-order upwind scheme exhibit a larger temperature increase, of the order of 2%, behind the oblique shock wave. Moreover, the 2nd-order upwind solutions converge to machine zero only when the minmod limiter is used. In this case, both the 1st-order upwind and the centered schemes also converge to machine zero.

Computational results for a supersonic reactive flow case were also presented and discussed. In this particular case, the simulation considered the flow over a wedge in which a detailed hydrogen-air chemistry is included. The solutions obtained with both the 1st- and the 2nd-order upwind schemes were compared with the one calculated using a well-known 2nd-order structured grid solver. The results for the 1st-order upwind method are in good agreement with those calculated using the later code, even if the solution is smeared in the first case due to the 1st-order spatial accuracy. The 2nd-order, unstructured upwind method exhibits, as it could be expected from the inert gas results, some discrepancies. In particular, the chemical kinetics induction length is smaller than observed in the structured grid solution or in the 1st-order unstructured grid results. This discrepancy is consistent with the observed higher temperatures for the 2nd-order unstructured solution in the inert gas case.

Acknowledgements

This work was accomplished in the framework of an international cooperation agreement between Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) and Centre National de la Recherche Scientifique (CNRS). The authors are indebted to the support of both institutions which has enabled a very fruitful interaction between the research groups represented in the present paper. The authors also gratefully acknowledge the additional partial support of CNPq under the Integrated Project Research Grants No. 530109/93-0 and 522413/96-0

Manuscript received: January 1998, Technical Editor: Angela Ourívio Nieckele.

• Anderson, W. K., Thomas, J. L., and van Leer, B., 1986, "A Comparison of Finite Volume Flux Vector Splittings for the Euler Equations,"AIAA Journal, Vol. 24, No. 9, pp. 1453-1460.
• Azevedo, J. L. F., 1992, "On the Development of Unstructured Grid Finite Volume Solvers for High Speed Flows,"Report NT-075-ASE-N/92, Instituto de Aeronáutica e Espaço, São José dos Campos, SP, Brazil.
• Azevedo, J. L. F., and Korzenowski, H., 1996, "Analysis of Hypersonic Flows Using Unstructured Meshes,"(in Portuguese), Proceedings of the 6th Brazilian Congress of Engineering and Thermal Sciences and 6th Latin American Congress of Heat and Mass Transfer - ENCIT/LATCYM 96, Vol. I, Florianópolis, SC, Brazil, pp. 547-552.
• Azevedo, J. L. F, and Mitchell, V. C. G., 1995, "A Comparison of Spatial Discretization Schemes for Unstructured Grid Finite Volume Calculations,"Proceedings of the 16th Iberian Latin American Congress on Computational Methods for Engineering - CILAMCE 95, Vol. 1, Curitiba, PR, Brazil, pp. 518-526.
• Balakrishnan, G., and Williams, F. A., 1994, "Turbulent Combustion Regimes for Hypersonic Propulsion Employing Hydrogen-Air Diffusion Flames,"Journal of Propulsion and Power, Vol. 10, No. 3, pp. 434-436.
• Batina, J.T., 1989, "Unsteady Euler Airfoil Solutions Using Unstructured Dynamic Meshes,"AIAA Paper 89-0115, AIAA 27th Aerospace Sciences Meeting, Reno, NV, USA.
• Batina, J. T., 1991, "Unstructured-Grid Methods Development - Lessons Learned,"Proceedings of the 4th International Symposium on Computational Fluid Dynamics, Vol. I, Davis, CA, USA, pp. 91-96.
• Byrne, G. D., and Dean, A.M., 1993, "The Numerical Solution of Some Kinetics Models with VODE and CHEMKIN II,"Computers Chem., Vol. 17, No. 3, pp. 297-302.
• Figueira da Silva, L. F., and Deshaies, B., 1995, "Numerical Analysis of the Self-Ignition of Hydrogen-Air Supersonic Flows over a Wedge,"Proceedings of the Third Asian-Pacific International Symposium on Combustion and Energy Utilisation, Vol. II, Hong-Kong, Hong-Kong, pp. 673-678.
• Hirsch, C., 1990, Numerical Computation of Internal and External Flows. Vol. 2: Computational Methods for Inviscid and Viscous Flows, Wiley, New York, USA (Chapt. 20, pp. 408-443).
• Jameson, A., 1987, "Successes and Challenges in Computational Aerodynamics,"AIAA Paper 87-1184, American Institute of Aeronautics and Astronautics, USA.
• Jameson, A., and Baker, T. J., 1983, "Solution of the Euler Equations for Complex Configurations,"AIAA Paper 83-1929, American Institute of Aeronautics and Astronautics, USA.
• Jameson, A., and Baker, T. J., 1987, "Improvements to the Aircraft Euler Method,"AIAA Paper 87-0452, AIAA 25th Aerospace Sciences Meeting, Reno, NV, USA.
• Jameson, A., and Mavriplis, D., 1986, "Finite Volume Solution of the Two-Dimensional Euler Equations on a Regular Triangular Mesh,"AIAA Journal, Vol. 24, No. 4, pp. 611-618.
• Jameson, A., Schmidt, W., and Turkel, E., 1981, "Numerical Solution of the Euler Equations by Finite Volume Methods Using Runge-Kutta Time-Stepping Schemes,"AIAA Paper 81-1259, AIAA 14th Fluid and Plasma Dynamics Conference, Palo Alto, CA, USA.
• Kee, R. J., Rupley, F. M., and Miller, J.A., 1991, "CHEMKIN-II: A Fortran Chemical Kinetics Package for the Analysis of Gas Phase Chemical Kinetics,"Sandia National Laboratories, SAND89-8009B/UC-706, USA.
• LeVeque, R. J., and Yee, H. C., 1990, "A Study of Numerical Methods for Hyperbolic Conservation Laws with Stiff Source Terms,"Journal of Computational Physics, Vol. 86, No. 1, pp. 187-210.
• Li, C., Kailasanath, K., and Oran, E., 1994, "Detonation Structures Behind Oblique Shocks,"Physics of Fluids, Vol. 6, No. 4, pp. 1600-1611.
• Liou, M.-S., 1996, "A Sequel to AUSM: AUSM⁺,"Journal of Computational Physics, Vol. 129, No. 2, pp. 364-382.
• Lyra, P. R. M., "Unstructured Grid Adaptive Algorithms for Fluid Dynamics and Heat Conduction,"Ph.D. Thesis, Department of Civil Engineering, University of Wales Swansea, Swansea, Wales, U.K..
• Mavriplis, D. J., 1988, "Multigrid Solution of the Two-Dimensional Euler Equations on Unstructured Triangular Meshes,"AIAA Journal, Vol. 26, No. 7, pp. 824-831.
• Mavriplis, D. J., 1990, "Accurate Multigrid Solution of the Euler Equations on Unstructured and Adaptive Meshes,"AIAA Journal, Vol. 28, No. 2, pp. 213-221.
• Pulliam, T. H., 1986, "Artificial Dissipation Models for the Euler Equations,"AIAA Journal, Vol. 24, No. 12, pp. 1931-1940.
• Strang, G., 1968, "On the Construction and Comparison of Difference Schemes,"SIAM J. Num. Anal., Vol. 5, pp. 506-517.
• van Leer, B., 1982, "Flux-Vector Splitting for the Euler Equations,"Proceedings of the 8th International Conference on Numerical Methods in Fluid Dynamics, E. Krause, editor, Lecture Notes in Physics, Vol. 170, Springer-Verlag, Berlin, Germany, pp. 507-512.
• Venkatakrishnan, V., 1995, "Convergence to Steady State Solutions of the Euler Equations on Unstructured Grids with Limiters,"Journal of Computational Physics, Vol. 118, pp. 120-130.
• Viguier, C., Figueira da Silva, L. F., Desbordes, D., and Deshaies, B., 1996, "Onset of Oblique Detonation Waves: Comparison Between Experimental and Numerical Results for Hydrogen-Air Mixtures,"Proceedings of the Twenty-Sixth Symposium (International) on Combustion, The Combustion Institute, Vol. 2, Napoli, Italy, pp. 3023-3031.

Publication Dates

• Publication in this collection
  15 Oct 2001
• Date of issue
  Dec 1999

History

• Received
  Jan 1998