Numerical Simulation of Combustion and Rotor-Stator Interaction in a Turbine Combustor

This article presents the development of a numerical algorithm for the computation of flow and combustion in a turbine combustor. The flow and combustion are modeled by the Reynolds-averaged Navier-Stokes equations coupled with the species-conservation equations. The chemistry model used herein is a two-step, global, finite-rate combustion model for methane and combustion gases. The governing equations are written in the strong conservation form and solved using a fully implicit, finite-difference approximation. The gas dynamics and chemistry equations are fully decoupled. A correction technique has been developed to enforce the conservation of mass fractions. The numerical algorithm developed herein has been used to investigate the flow and combustion in a one-stage turbine combustor.

signers are facing the fact that the combustor residence time can become shorter than the time required to complete combustion.Consequently, combustion would continue in the turbine, which until recently has been considered undesirable.The process of combustion in the turbine is called in situ reheat, and the turbine in which the combustion takes place is called the turbine burner.Herein, fuel is injected at the trailing edge of the stator.A thermodynamic cycle analysis demonstrates performance gains in a turbojet engine with a turbine burner (Sirignano and Liu, 1999).Even better performance gains in specific power and thermal efficiency are predicted for a powergeneration gas-turbine engine when the turbine is coupled with a heat regenerator.
Several challenges, however, are associated with combustion in the turbineburner: mixed subsonic and supersonic flows; flows with large unsteadiness caused by the rotating blades; hydrodynamic instabilities; and intense straining of the flow due to the very large three-dimensional acceleration and stratified mixtures (Sirignano and Liu, 1999).The obvious drawback associated with the strained flows in the turbine burner is that widely varying velocities can result in widely varying residence time for different flow paths, and as a result there are flammability difficulties in regions with shorter residence times.In addition, transverse variation in velocity and kinetic energy can cause variations in entropy and stagnation entropy that impact heat transfer.The heat transfer and mixing may be enhanced by increasing the interface area caused by the strained flows.
Turbine aerodynamics might be drastically modified by strong exothermic combustion processes in a turbine burner.Thermal expansion due to combustion could significantly change the pressure variation and the shock strength and location.As a result, the blade loading would be modified.There is evidence that in a low-pressure turbine without in situ reheating, the temperature nonuniformities can generate strong entropic and vortical waves.These waves produce excitations large enough to 364 D. D. ISVORANU AND P. G. A. CIZMAS generate unsteady loadings and stresses on the fifth stage of a low-pressure turbine, stresses sufficient to cause high-cycle fatigue failures (Manwaring and Kirkeng, 1997).The danger of high-cycle fatigue is even greater in a turbine burner because larger temperature nonuniformities are likely to produce stronger entropic and vortical waves.
Experimental data for conventional (i.e., without in situ reheating) gas turbines have shown the existence of large radial and circumferential temperature gradients downstream of the combustor (Dills and Follansbee, 1979;Elmore, Robinson, and Watkins, 1983).These temperature nonuniformities, called hot streaks, have a significant impact on the secondary flow and wall temperature of the entire turbine.Since the combustor exit flow may contain regions where the temperature exceeds the allowable metal temperature by 260 • C to 520 • C (Butler et al., 1989), understanding the effects of temperature nonuniformities on the flow and heat transfer in the turbine is essential for increasing vane and blade durability.It is estimated that an error of 55 • C in predicting the time-averaged temperature in a turbine rotor can result in an order-of-magnitude change in a blade's life (Graham, 1980;Kirtley et al., 1993).
Temperature nonuniformities are amplified in a turbineburner.Consequently, it is expected that not only will the secondary flow and wall temperature be affected but also the blade loading, as a result of the modified pressure distribution.Temperature nonuniformities in a turbine burner also affect the locations of hot spots on airfoils and as a result impact the internal and film-cooling schemes.
Extensive experimental studies (Butler et al., 1989;Schwab et al., 1983;Shang et al., 1995;Sharma et al., 1992;Stabe et al., 1984;Whitney et al., 1980) and numerical studies (Dorney et al., 1990(Dorney et al., , 1996(Dorney et al., , 1999(Dorney et al., , 2000;;Krouthen and Giles, 1988;Rai and Dring, 1990;Shang and Epstein, 1996;Takahashi and Ni, 1991) have been made of the influence of temperature nonuniformities on the flow and heat transfer in a conventional turbine.To the best knowledge of the authors, however, no data are available in the open literature concerning the effect of in situ reheating in turbine burners.
The objective of this article is to present the development of an integrated fluid-flow and combustion-numerical model.This model is used to simulate numerically the influence of in situ reheating on turbine aerodynamics.The development of this numerical model is crucial to the development of turbine burners which, in spite of their challenges, can provide significant performance gains for turbojet engines and power-generation gas turbine engines.
The next section presents the physical model used for the simulation of flow and combustion in the turbine combustor.The governing equations and the chemistry model are presented.The third section describes the numerical model.That section includes information about the grid generation, boundary conditions, numerical method, and parallel algorithm.The results are presented in the fourth section.

PHYSICAL MODEL
The flow through and combustion in a multirow turbine burner with arbitrary blade counts are modeled by the Reynoldsaveraged Navier-Stokes equations and the species conservation equations.To reduce the computation time, the flow and combustion are modeled as being quasi-three-dimensional.This section presents the details of the governing equations and the chemistry model.

Governing Equations
The unsteady, compressible flow through a turbine combustor is modeled by the Reynolds-averaged Navier-Stokes equations.The flow is assumed to be fully turbulent and the kinematic viscosity is computed using Sutherland's law.The eddy viscosity is computed using the Baldwin-Lomax algebraic turbulence model.The Reynolds-averaged Navier-Stokes equations and species conservation equations are simplified by using the thin-layer assumption.For turbulent flows, the ratio of the neglected to the retained terms is of the order of δ/L, where δ is the thickness of the shear layer at distance L from the origin (Bradshaw, 1996).Typically, the neglected streamwise molecular diffusion and turbulent diffusion terms are of the order of 5.3Re −1 and 0.37Re −1/5 for Re ≤ 10 7 , respectively (Peters et al., 1986).
The dimensional variables are nondimensionalized using free stream quantities: In the hypothesis of unity Lewis number, both the Reynoldsaveraged Navier-Stokes and the species equations can be written according to Balakrishnan (1987): Note that Equation (1) is written in the body-fitted curvilinear coordinate system (ξ, η, τ ).The state and flux vectors in the (ξ, η, τ ) coordinates are written as a function of the state and flux vectors in the Cartesian coordinates (x, y, t): where the J is the Jacobian of the coordinate transformation, J = 1/(ξ x η y − η x ξ y ).The state and flux vectors of the Reynolds-averaged Navier-Stokes equations in the Cartesian coordinates are The state and flux vectors of the species conservation equations in the Cartesian coordinates are The viscous terms are: and where The chemical source terms are: where ŵi is the species net production rate.The set of equations is completed by the equation of state, p = ρ RT /γ , the caloric equation, e = p/(γ − 1) + 0.5ρ(u 2 + v 2 ), the Sutherland's formula for computing molecular viscosity, µ = µ 0 (T 0 + S 1 )/(T + S 1 )(T /T 0 ) 3/2 , and the definition of the adiabatic exponent γ = The temperature values used herein are T 0 = 273 K and T ref = 298 K.The specific heats c p and c v are modeled by a fourth-degree temperature polynomial (Kee et al., 1989).

Chemistry Model
The chemistry model used herein to simulate the in situ reheating is a two-step, global, finite-rate combustion model for methane and combustion gases (Hautman et al., 1981;Westbrook and Dryer, 1981 The rate of progress (or Arrhenius-like reaction rate) for methane oxidation is given by: where The reaction rate for the CO/CO 2 equilibrium is: with A 2 = 2.249 • 10 12 (m 3 /kmol) 0.75 s −1 and E 2 /R M = 20130 K.The symbols in the square brackets represent local molar concentrations of various species.The net formation/destruction rate of each species due to all reactions is: where ν ik are the generalized stoichiometric coefficients.Note that the generalized stoichiometric coefficient is ν ik = ν ik − ν ik , where ν ik and ν ik are stoichiometric coefficients for species i in reaction k appearing as a reactant or as a product.Equation (3) models the flame speed variation with pressure and equivalence ratio φ, in the range of 1-25 atmospheres and 0.5-1.5, respectively (Westbrook and Dryer, 1981).The negative exponent of the methane molar concentration accounts for the ignition inhibition, as seen in shock tubes (Bowman, 1970).The negative exponent of the methane molar concentration may create numerical problems because the rate of methane consumption exhibits an unbounded increase as fuel concentration approaches zero.There are several possible solutions for this problem.One option is to use a reverse reaction that would provide an equilibrium fuel concentration at some small level.Another option is to artificially truncate the rate of progress at some predetermined fuel concentration (Westbrook and Dryer, 1981).The latter option has been used herein.
This chemistry model does not take into account the influence of turbulence on the reaction rates.The influence of the flow field on species transport equations is modeled only by eddy diffusivity.The nonunity exponents of fuel and oxidizer concentration, however, make this reaction mechanism useful even for turbulent flows (Westbrook and Dryer, 1981).Successful simulation of turbulent combustion using a single-step global reaction rate has been reported (Butler et al., 1980;Gupta et al., 1980).Consequently, it is expected that the two-step, global, finite-rate combustion model used herein is adequate to capture the salient features of the in situ reheating in the turbine.

NUMERICAL MODEL
The numerical model used herein is based on an existing algorithm developed for unsteady flows in turbomachinery (Cizmas and Subramanya, 1997).The Reynolds-averaged Navier-Stokes equations and the species equations are written in the strong conservation form.The fully implicit, finite-difference approximation is solved iteratively at each time level, using an approximate factorization method.Three Newton-Raphson subiterations are used to reduce the linearization and factorization errors at each time step.The convective terms are evaluated using a thirdorder accurate, upwind-biased Roe scheme.The viscous terms are evaluated using second-order accurate central differences.The scheme is second-order accurate in time.

Grid Generation
Two types of grids are used to discretize the flow field surrounding the rotating and stationary airfoils, as shown in Figure 1.An O-grid is used to resolve the governing equations near the airfoil, where the viscous effects are important.An H-grid is used to discretize the governing equations away from the airfoil.The O-grid is generated using an elliptical method.The H-grid is algebraically generated.The O-and H-grids are overlaid.The flow variables are communicated between the O-and H-grids through bilinear interpolation.The H-grids corresponding to consecutive rotor and stator airfoils are allowed to slip past each other to simulate the relative motion.

Discretization of Governing Equations
The gas-dynamics and chemistry governing equations can be solved coupled or decoupled.The decoupled chemistry and

FIGURE 1
Detail of the coarse grid; every other grid point in each direction is shown.
gas-dynamics governing equations are solved in a successive sequence.The decoupled governing equations can be either fully or partially decoupled (Yee, 1987).In the latter technique, the flow and chemistry equations are solved separately while retaining the local characteristics of the fully coupled methods.Under certain conditions, the fully decoupled methods, also called chemistry split methods (Balakrishnan, 1987), exhibit mass fractions nonconservative behavior (Eberhardt and Brown, 1986).To prevent this, either correction techniques (Eberhardt and Brown, 1986) or block techniques (Li, 1987) are employed.Herein, the fully decoupled implicit algorithm has been chosen.A correction technique has been developed to enforce the conservation of mass fractions.The governing equations are discretized using an implicit, approximate-factorization, finite-difference scheme in delta form (Warming and Beam, 1978).The discretized operational form of both the Reynolds-averaged Navier-Stokes and the species conservation equations, combined in a Newton-Raphson algorithm (Rai and Chakravarthy, 1986), is: where A and B are the flux Jacobian matrixes Note that the flux Jacobian matrixes are split into A = A + + A − , where A ± = P ± P −1 . is the spectral matrix of A and P is the modal matrix of A. The spectral matrix is split into = + + − , where the components of + and − are λ − i = 0.5(λ i −|λ i |) and λ + i = 0.5(λ i +|λ i |), respectively (Steger and Warming, 1981).The same flux vector splitting approach is applied to matrix B. In Equation ( 5), , ∇, and δ are forward, backward, and central differences operators, respectively.Q p is an approximation of Q n+1 .At any time step n, the value of Q p varies from Q n at the first internal iteration, when p = 0, to Q n+1 , when the integration of Equation ( 5) has converged.
The intercell numerical flux is expressed according to Toro (1999): where Roe's approximate Riemann solver is used to compute cellface numerical flux differences F (Roe, 1981).Averaged cell interface values are determined for density, velocity, enthalpy, and mass fractions: The numerical flux differences are expressed in terms of the primitive variables' jump across the intercell boundary and the right eigenvectors Pk : where m represents the number of eigenvalues and the ᾱk parameters are determined from the linear system Q = m k=1 ᾱk Pk .The appropriate coefficients of the linear system are calculated for either the Navier-Stokes or the species conservation equations.
Roe's numerical approach is not consistent with the entropy inequality and might converge to some nonphysical solution.To avoid this problem, the coefficient of the numerical viscosity term is adjusted by modifying the flux difference term (Yee and Harten, 1987;Yee and Shinn, 1989).Note that the elements of the Jacobian matrix in Equation ( 7) depend on the system's eigenvalues.To correct the problem of entropy violation, the Jacobian matrix Ā( λk ) is replaced by Ā( λk corr ), where Here, ε is a small number.A correction algorithm is proposed herein to enforce the mass fraction conservation.This correction algorithm is similar in principle to the methodology presented by Eberhardt and Brown (1986).The proposed algorithm corrects the mass fractions of the reacting species such that the sum of nonreacting species is constant.Here, y nr represents the sum of nonreacting species, and y c i is the corrected mass fractions.

Boundary Conditions
Because multiple grids are used to discretize the governing equations, two classes of boundary conditions must be enforced on the grid boundaries: natural boundary conditions and zonal boundary conditions.The natural boundaries include inlet, outlet, periodic, and the airfoil surfaces.The zonal boundaries include the patched and the overlaid boundaries.
The inlet boundary conditions include the specification of the flow angle, average total pressure, and downstream propagating Riemann invariant.The upstream propagating Riemann invariant is extrapolated from the interior of the domain.At the outlet, the average static pressure is specified, while the downstream propagating Riemann invariant, circumferential velocity, and entropy are extrapolated from the interior of the domain (Jameson and Yoon, 1985;Rai and Chaussee, 1984).Periodicity is enforced by matching flow conditions between the lower surface of the lowest H-grid of a row and the upper surface of the topmost H-grid of the same row.At the airfoil surface, the following boundary conditions are enforced: the "no slip" condition, the adiabatic wall condition, and the zero normal pressure gradient condition.
For the zonal boundary conditions of the overlaid boundaries, data are transferred from the H-grid to the O-grid along the Ogrid's outermost grid line.Data are then transferred back to the H-grid along its inner boundary.At the end of each iteration, an explicit, corrective, interpolation procedure is performed.The patch boundaries are treated similarly, using linear interpolation to update data between adjoining grids (Rai, 1985).
Additional boundary conditions are imposed at the vane's trailing edge to model the plane jet of pure fuel that is injected in the turbine combustor.The following variables are specified at the fuel injection hole: normal velocity, temperature, pressure, and species concentrations.The species concentrations are also specified at the inlet in the turbine burner.

Parallel Computation
The parallel code uses message-passing interface libraries and runs on different parallel platforms, from the Beowulftype PC cluster to the Cray T3E.One processor is allocated for each inlet and outlet H-grid.One processor is allocated for the O-and H-grids corresponding to each airfoil.Interprocessor  communication is used to match boundary conditions between grids.Periodic boundary conditions are imposed by cyclic communication patterns within rows.Inter-blade-row boundary conditions are imposed by gather-send-receive-broadcast communication routines between adjacent rows.Load imbalance issues need to be considered at grid generation time so as to reduce synchronization overhead (Cizmas and Subramanya, 2001).

FIGURE 3
Skin friction variation on the coarse, medium, and fine grids.

RESULTS
The computation algorithm presented in the previous section was implemented in the CoRSI code.This code was built on an existing Reynolds-averaged Navier-Stokes solver (Cizmas and Subramanya, 1997).Herein, the code was used to simulate the unsteady flow and combustion in a one-stage turbine burner.This section begins with the validation of the accuracy of the  numerical results.This is followed by a set of results that illustrate the salient features of the interaction between combustion and rotor-stator interaction in a turbine burner.

Geometry and Flow Conditions
The one-stage turbine burner has 32 vanes and 49 blades.A dimensionally accurate computation of this geometry would require full-annulus simulation.To reduce the computational effort, it was assumed that there were an equal number of airfoils (32) in each turbine row.As a result, the rotor airfoils were rescaled by a 49/32 factor.An investigation of the influence of the airfoil count on the turbine flow showed that the unsteady effects were amplified when a simplified airfoil count of 1:1 was used (Cizmas, 1999).Consequently, the results obtained using the simplified airfoil count represent an upper limit of the unsteady effects.
The inlet temperature in the turbine combustor exceeds 1800 K and the inlet Mach number is 0. The fuel is injected at the trailing edge of the stator airfoil.The width of the injection hole is 0.0098 of the axial chord length.The Reynolds number based on the injection conditions at the trailing edge is 5700.The inlet temperature of the jet is 563 K and the axial velocity is 0.098 of the stage inlet velocity.

Accuracy of Numerical Results
To validate the accuracy of the numerical results it was necessary to show that the results were independent of the grid, which discretizes the computational domain.Three grids were used to assess the grid's independence of the solution.The coarse grid had 35 grid points normal to the airfoil and 150 grid points along the airfoil in the O-grid; there were 50 grid points in the axial direction and 50 grid points in the circumferential direction in the H-grid.The stator airfoil and rotor airfoil had the same number of grid points.The number of grid points in the medium and fine grids is presented in Table 1.The coarse grid is presented in Figure 1, where for clarity every other grid point in each direction is shown.The distance between the grid points on the airfoil and the next layer of grid points around the airfoil was the same for the coarse, medium, and fine grids in order to have the same y + number.The grid was generated such that, for the given flow conditions, the y + number was less than 1.Approximately 20 grid points were used to discretize the boundary layers.

FIGURE 7
Variation of Mach number.
The results presented in this article were computed using three Newton subiterations per time-step and 3000 time steps per cycle.Here, a cycle is defined as the time required for a rotor to travel a distance equal to the pitch length at midspan.To ensure time periodicity, each simulation was run in excess of 80 cycles.
The flow in the rotor row included the influences of the upwind stator row.Consequently, any differences that exist between the results caused by different grid sizes would be the largest in the rotor row.For this reason, the rotor row of the turbine burner was used to assess the grid's independence of the numerical results.The nondimensional skin friction ) was used to validate the independence of the grid solution.
Before validating the grid's independence of the numerical results, we had to verify that the unsteady solution was periodic.Solution periodicity was assessed by comparing the results of consecutive cycles, as shown in Figure 2. Since the values of the skin friction were almost identical for the three consecutive cycles, it was concluded that the solution was periodic.
To validate the grid's independence, three values of the skin friction were compared: the averaged, the minimum, and the maximum over one period.The comparison of the skin friction coefficients computed using the three grids is shown in Figure 3.
Good agreement was obtained between the results corresponding to the medium and fine grids, the maximum difference being less than 8% of the maximum skin friction value.The difference between the coarse grid skin friction and the medium or fine grid skin friction was significantly larger.The largest difference was approximately 17% of the maximum skin friction value, shown on the pressure side of the maximum skin friction.As a result, it was concluded that the numerical results were independent of the grid only for the medium and fine grids.The medium grid was used for the computation of all subsequent results throughout the study.

Variation in Gas-Dynamics and Chemistry Variables
The ignition of the cold, slow fuel jet diffused into the hot, fast stream of oxidizer is delayed until the mixture achieves both the necessary local oxidizer/fuel ratio and the temperature to boost the chemical mechanism.The enhanced turbulence in the wake of the stator airfoil's trailing edge accelerates the mixing between the two streams.The flame front starts just before entering the rotor and moves downstream.The flame is broken into patches of burning mixture by the intermittent passage of the rotor's leading edges.The patches of burning mixture slide along the rotor airfoil, continuing to burn due to the lower velocities in the boundary layer.These patches of burning mixture expand their volume because of the rise in temperature and the drop in pressure.Coherent flame structures are formed downstream of rotor's trailing edge, as shown in Figures 4 and 5.

CONCLUSIONS
The complexity of the reacting flow in the turbine-combustor stage makes it one of the most challenging numerical simulation problems.The large unsteadiness and straining of the flow along with the wide range of variation in velocity are most likely to lead to a wide spread of local characteristic time scales that will strongly impact the ongoing reactions.As the first step in the numerical simulation of the in situ reheating, a simple chemistry mechanism was considered.The numerical simulation proved that the two-step, global, finite-rate combustion model used herein captures the salient features of the in situ reheating in the turbine combustor.This numerical simulation is currently being used to analyze the influence of different parameters on the flow and combustion in a turbine combustor.

FIGURE 2
FIGURE 2Skin friction at three consecutive cycles.
FIGURE 4Variation of CO.
155.The inlet flow angle is 0 degrees and the inlet Reynolds number is 65799 per inch, based on the axial chord of the first-stage stator.The values of the species concentrations at the inlet in the turbine burner are: y CO 2 = 0.0775, y H 2 O = 0.068, y CO = 5.98 • 10 −06 , y H 2 = 2.53 • 10 −07 , y O 2 = 0.1131, y N 2 = 0.7288, and y Ar = 0.0125.The rotational speed of the test turbine burner is 3600 rpm.

TABLE 1
Grid Points for the Single-Stage Turbine Burner