QUASI-TWO-DIMENSIONAL APPROXIMATION FOR NUMERICAL SIMULATION OF FLOWS IN ENGINE DUCTS

New quasi-two-dimensional (2.5D) approach to description of three-dimensional (3D) §ows in ducts is proposed. It generalizes quasi-one-dimensional (quasi-1D, 1.5D) theories. Calculations are performed in the ( x ; y ) plane, but variable width of duct in the z direction is taken into account. Derivation of 2.5D approximation equations is given. Tests for veri¦cation of 2.5D calculations are proposed. Parametrical 2.5D calculations of §ow with hydrogen combustion in an elliptical combustor of a high-speed aircraft, investigated within HEèAFLY-INT international project, are described. Optimal scheme of fuel injection is found and explained. For one regime, 2.5D and 3D calculations are compared. The new approach is recommended for use during preliminary design of combustion chambers.


INTRODUCTION
Turbulent §ows of viscous gas in ducts (especially, in aircraft power plants) are, as a rule, essentially 3D, complex, multiscale phenomena that can contain boundary layers, recirculation zones, mixing layers and jets, compression and rarefaction waves, ¦nite-rate reactions and unsteady processes. Calculation of such §ows on the basis of full 3D Reynolds-averaged NavierStokes (RANS) equations requires huge computer resources even in the case of parallel computing. As a result, multiple parametrical calculations, which are necessary on the stage of combustor design, are impossible. Even now, simple engineering estimations, based on empirical relations, and calculations in quasi-1D approximation are often used for prediction of combustor characteristics for the design purposes. However, quasi-1D approximation cannot take into account some spe-ci¦c features resulting from non-1D character of §ow, e. g., nonuniformity of heat release depending on mixing of fuel with air §ow and on the heat turbulent transport across the duct. To improve accuracy and informativeness of approximate parametrical calculation on the stage of design, approximation of 2.5D §ow may be used.
The term ¤2.5D approximation¥ means that for the description of 3D processes, the 2D equations, taking into account §ow nonuniformity in the third spatial direction, are used. In this sense, classical equation system for axisymmetric §ows is also 2.5D approach for the description of 3D §ow. But it is applicable only if the problem geometry is axisymmetric and if assumption about the absence of §ow swirl is acceptable. Literature review has allowed to ¦nd other examples of 2.5D approach to description of 3D §ows. Probably, the best known is the ¤method of §at sections¥ for the description of §ow around in¦nite swept wing (see, for example, [14]). The 2.5D approach is used in aeroacoustics for the description of 3D sonic wave propagation, when the basic aerodynamic §ow has two-dimensional (2D) character [5,6]. There are techniques, analogous to method of §at sections, for description of strati¦ed §ows with given character of §ow variation from one layer to another [7,8]. But all these approaches cannot be used for the description of §ows in ducts. In this work, new 2.5D approach to the description of 3D §ows in ducts is proposed. This approach generalizes quasi-1D (¤1.5D¥) theories. In 2.5D approximation, the real 3D §ow is replaced by a §ow where all parameters are constant along z axis. It may be treated as a result of 3D §ow averaging along z. Calculations are performed in the (x; y) plane but variable width of duct in z direction is taken into account.

QUASI-TWO-DIMENSIONAL FLOW EQUATION SYSTEM
The §ow in a duct of arbitrary geometry is considered (Fig. 1). The duct geometry is projected on some plane. In this plane, Cartesian coordinate system (x; y) is introduced. Spatial direction of the third coordinate axis, z, is determined using the right-hand screw rule. Computational domain is projection of the duct geometry on the (x; y) plane. It is covered by computational grid (see Fig. 1). Spatial curve on the duct surface, corresponding to the contour of the duct projection on the (x; y) plane, subdivides the duct surface into two halves. These halves can be described by functions z + (x, y) (¤front¥ surface) and z − (x, y) (¤back¥ surface). For each pair of values (x; y), z + ≥ z − . If the duct is symmetrical about vertical plane, z + ≥ 0 and z − = −z + . Reynolds-averaged NavierStokes equation system for multicomponent gas with ¦nite-rate reactions, closed by di¨erential model of turbulence and by model Here, a is the quantity of some physical parameter (mass, momentum, energy, turbulence parameters, or mass of the reacting mixture chemical component) per unit volume of gas; F i is the §ux of a value along x i axis (x 1 ≡ x, x 2 ≡ y, and x 3 ≡ z); and W are the local sources of a. Inner space of the duct is subdivided into following elementary volumes: is the duct width in z direction. One such element, corresponding to an arbitrary cell of computational grid [x − -x/2, x + -x/2] × [y − -y/2, y + -y/2], is shown in Fig. 1. The 2.5D approximation is obtained, if one assumes that §ow parameters in each element are constant along z axis. Characteristics of 2.5D §ow in this cell may be treated as the result of real 3D §ow averaging along z direction: a(x, y, z) dz .
After that, let us integrate Eq. (1) over the constant volume de¦ned by Eq. (2), assuming §ow parameters to be constant along z direction and applying GaussOstrogradsky theorem: Here, � S + = S + x ; S + y ; S + z is the vector of outward normal to the volume ele-ment£s front side (coordinates of this side center are (x, y, z + (x, y))), and is the vector of outward normal to the volume element£s back side (coordinates of this side center are (x, y, z − (x, y))). Magnitudes of these vectors are equal to the areas of corresponding sides; magnitudes of the components of these vectors are equal to the areas of projections of these sides on planes, perpendicular to coordinate axes. Therefore, Then, let us substitute these relations into Eq. (4), divide it by -x-y, and consider the limit -x → 0, -y → 0: This is 2.5D analogue of Eq. (1). By writing Eq. (5) for a = ρ, ρu, ρv, ρE, ρp t k ; and ρY m , one gets the equation system for 2.5D analogue of 3D §ow in the duct. Here and below, ρ is the density; u and v are the velocity components along x and y axes, respectively; E is the total energy per unit mass of gas: Y m e m (T ) ; p t k , k = 1, . . . , N turb , are the parameters of the considered model of turbulence; Y m , m = 1, . . . , N sp , are the mass fractions of the reactive mixture components; k is the averaged kinetic energy of turbulence; e m (T ) is the internal energy of the mth component of gas mixture; and T is the temperature. In this work, qω turbulence model [911] is considered, for which N turb = 2, p t 1 ≡ q ≡ √ k is the characteristic value of turbulent §uctuations of velocity, and p t 2 ≡ ω ≡ ε/k is the characteristic frequency of turbulent §uctuations (ε is the average rate of turbulent kinetic energy dissipation). Chemical kinetics model, consisting of N sp = 9 components (H, O, OH, H 2 O, O 2 , H 2 , CO, CO 2 , and N 2 ) is used.
Essential moment in construction of 2.5D analogue of 3D §ow in the duct is the way to determine the §uxes through lateral sides F − i and F + i that are placed in additional source terms in the left-hand side of Eq. (5). As well as in construction of quasi-1D (¤1.5D¥) §ow theories, these §uxes should be determined taking into account the real 3D §ow structure near the duct walls.
Vector of §uxes along x i axis (x 1 = x, x 2 = y, and x 3 = z) on solid impermeable surface without slip can be represented as follows: Here, n is the component of unit vector of local outward normal to the duct wall � n ± (x, y) = � n ± x ; n ± y ; n ± z � ; V τ ≡ � V · � τ is the tangential-to-wall component of velocity (direction of � τ vector coincides with limit direction of velocity vector, when the point approaches to the wall); index ¤W¥ marks the values of parameters on the duct wall. Components of the unit outward normal to wall � n ± can be expressed as follows: Components of unit vector � τ , which is parallel to tangential velocity near wall, can be chosen on the basis of assumption that projection of the tangen-tial velocity vector on (x; y) plane is codirectional to velocity vector of 2.5D (z-averaged) §ow that has components (u(x, y), v(x, y)): (τ x , τ y )�(u, v). From condition � τ · � n = 0, one may ¦nd τ z = −(τ x n x + τ y n y )/n z . Finally, let us normalize � τ vector and get: In determination of wall §uxes, let us use the model of 3D §ow consisting of inviscid core, where the pressure is constant along z direction, and of boundary layers, where the §ow is decelerated to zero velocity. Due to the fact that the pressure is practically constant across boundary layer, let us take wall pressure to be equal to pressure in the inviscid core and to z-averaged value of pressure (i. e., with pressure in 2.5D analogue of this 3D §ow). In this case, the wall pressure p W is equal to pressure of 2.5D §ow in current spatial element p(x, y).
Molecular di¨usive §uxes of momentum, heat, and turbulence parameters in the direction of the local normal to the wall (in Eq. (6), they have form µ W ∂f /∂n) are determined by local structure of boundary layer.
In 2.5D calculation, boundary layers appear only near upper and lower contours of the duct. Let us estimate the local values of µ W ∂f /∂n through linear interpolation in y between upper and lower contours of the duct. Therefore, let us take where y + (x) and y − (x) are the coordinates of the duct contours in the computational domain. It is necessary to underline that n-direction is determined by local outward normal to the duct wall; it is di¨erent for upper contour of duct, for its lower contour, and for current point of computational domain (x, y).
Other terms of Eq. (5) will be calculated using the same formulas as the corresponding terms in the initial equation system of 3D §ow. But we shall substitute parameters of 2.5D §ow (in fact, z-averaged parameters of real §ow) into these formulas. Naturally, this approach to description of 3D §ow is approximate, because we assume that z-averaged §uxes and sources in equation of gas motion can be calculated by usual formulas, where z-averaged gas parameters are substituted. Average of nonlinear function does nÏt coincide with the value of this function after substitution of average values of its arguments. Therefore, §uxes and source terms will be determined with some errors. These errors, together with errors of approximate formulas (7), constitute the inaccuracy of 2.5D approximation.
The most important target function in simulation of §ow in combustion chamber is the longitudinal force R. This is x-component of the integral force that is applied to the duct surface. By de¦nition, this force is equal to integral of momentum x-component over the duct surface. Accordingly, this force should be calculated as follows: In this formula, F i (ρu) is the §ux of momentum x-component along x i axis.
The ¦rst integral is calculated over the lower contour of the duct and the second integral ¡ over its upper contour. In these integrals, � dl = (dl x ; dl y ) is the vector with magnitude equal to the length of the duct contour element, this vector is codirectional with outward normal to the contour; h z dl x is the projection of the duct surface part, corresponding to this contour element, on the plane that is perpendicular to x axis; Ánd h z dl y is its projection on the plane that is perpendicular to y axis. In practice, sides of near-boundary computational cells are used as the contour elements.
It is worth to note that h z � = 0, only if the duct has upper or lower walls parallel to the z axis. In the example, which is shown in Fig. 1, such walls are absent, and ¦rst two integrals in Eq. (8) are equal to zero.
The third (surface) integral is taken over the inner part of the computational domain (of the duct projection on the (x; y) plane). In this integral, ds(x, y) is the area of grid cell; and (∂z − /∂x) ds, (∂z − /∂y) ds, and (−ds) are the components of the vector of the duct back surface element that corresponds to this element of computational domain. The length of this vector is equal to the area of the surface element and its direction coincides with direction of outward normal to surface; and (−(∂z + /∂x) ds), (−(∂z − /∂y) ds), and ds are the components of the vector of the duct front surface element that corresponds to this grid cell.

VERIFICATION OF CODE FOR QUASI-TWO-DIMENSIONAL CALCULATIONS
New possibility to simulate §ows in ducts in 2.5D approximation has been realized in scienti¦c code SOLVER3 [12] that is intended for calculation of 2D §ows of multicomponent gas on the basis of RANS equations (previously, this code allowed to calculate only planar and axisymmetric §ows). This section describes the test cases that were used to verify the module of SOLVER3 code where the new technique was realized.
The ¦rst test is based on the fact that the new approach to simulation of 3D §ows generalizes the classical theory of quasi-1D §ows, which is described, e. g., in [13]. In particular, theory of Laval nozzle is based on the equations of quasi-1D §ow.
Theory of Laval nozzle considers a duct with given law of area variation F (x). Flow in the duct is assumed to be stationary and inviscid. Quasi-1D (¤1.5D¥) analogue of this §ow can be obtained, if one assumes that §ow parameters are constant in each cross section of the duct, but vary from one section to another. In fact, it means that the parameters of real §ow are averaged over the duct section.
Equations of the Laval nozzle theory can be considered as particular case of 2.5D §ow equations (5). To deduce these equations from (5), the 2.5D theory assumption about constant §ow parameters along z axis should be supplied with assumptions that stationary §ow of single-component perfect gas is considered in a duct that has constant width -y = const in (x, y) plane, while its lateral width is changed according to law With the use of these additional assumptions, Eqs. (5) can be rewritten as follows: Taking p W = p and h z (x) = F (x)/-y, after transformations, one can obtain the classical equations of Laval nozzle theory: Ordinary di¨erential equation with two closing relations (9) can be solved numerically with very high accuracy and can be used as reference. Calculation of §ow in duct with -y = const, h z (x) = F (x)/-y on the basis of Eqs. (5) should give the solution that coincides with the reference solution within the approximation errors.
As the second test for the veri¦cation of 2.5D code, calculations of viscous turbulent §ow in axisymmetric duct with supersonic §ow at the entrance have been used. In the (x, y) plane, uniform grid, containing 398 cells along the duct and 70 cells across the duct, have been constructed. Uniform §ow of air with speci¦c heat ratio γ = 1.4, with Mach number M ≈ 2.9, with stagnation parameters p 0 ≈ 15 atm and T 0 = 2150 K and with turbulence parameters q = 16 m/s and ω = 2000 Hz. On the duct walls, boundary condition of the class ¤wall functions¥ [12,14] was used (it allows to avoid extreme compression of grid to no-slip walls). Parameters in the exit section of the duct were taken from the near-boundary cells.  (1) reference ¦eld obtained by solution of axisymmetric §ow equations; (2) ¦eld obtained by averaging of the reference ¦eld along z coordinate (see formula (3)); and (3) result of the same duct calculation in 2.5D approximation.
One can see that results of 2.5D calculation are close to the parameters, obtained by averaging of the reference ¦eld along the lateral coordinate. Minor di¨erences of these ¦elds are resulted from the presence of 3D stationary wave structures in supersonic §ow. Such structures cannot be reproduced correctly, if §ow parameters are considered to be constant along z axis. Thickness of boundary layer in 2.5D calculation is higher than in the axisymmetric calculation. But averaging of axisymmetric ¦eld along z gives boundary layer of practically the same thickness as in the 2.5D calculation. Figure 3Á compares the pressure distributions along the duct walls, obtained in these calculations. The 2.5D calculation predicts the longitudinal distribution of loads on the duct walls with good accuracy.

CALCULATIONS OF FLOW IN COMBUSTOR OF HIGH-SPEED AIRCRAFT
A 2.5D approach has been applied to parametrical calculations of §ow in the combustion chamber of a scramjet with hydrogen fuel for a hypothetic supersonic It is a duct with elliptical sections with two zones of fuel injection (Fig. 4). In the 1st zone, hydrogen is injected upwards from two short pylons, while in the 2nd zone, in the §ow direction from several holes in a high pylon placed in the symmetry plane of the duct. At the entrance of combustor, there is in §ow of air heated by a §ame heater and enriched with oxygen (to get the same mass fraction of oxygen as in air). Mach number in the §ow inviscid core M ∼ 2.6, pressure is about 0.5 atm. Temperature in the inviscid core at the entrance is close to 1200 K but temperature of the injected jets of hydrogen is equal to 163 K only. Initially, preliminary 3D calculation has been performed with hydrogen injection, but without chemical reactions. The 3D calculations were performed using scienti¦c code ZEUS-S3pp that is one component of EWT-TsAGI software package [19]. Calculation has been performed for the variant of fuel injection where 10% of the total hydrogen mass- §ow rate were injected through each pylon of the 1st injection zone and the rest 80% were injected through the holes of the 2nd injection zone pylon (scheme of fuel injection: 10%10%80%).
In Fig. 5, the ¦elds of oxidizer excess ratio, obtained in this calculation in several cross sections of the chamber, are shown. One can see that cross sections of fuel jets, injected from the central pylon, are spread in vertical direction. It is resulted from the fact that injected hydrogen enters into the wake past the central pylon. In the case of greater part of mass- §ow rate through the 1st injection zone pylons, cross sections of jets from the 1st injection zone pylons should also be spread in vertical direction due to the upward injection of fuel. Therefore, §ow parameters vary along y axis weaker than along the z axis. Consequently, there is sense to perform 2.5D calculations in the (x; z) plane instead of (x; y) plane. The 2.5D calculations in (x; y) plane give a picture of fuel mixing with air that di¨ers principally from the picture obtained in 3D calculations. So, below, only 2.5D calculations in (x; z) plane will be considered. These calculations have been performed with ¦nite-rate chemical reactions.
In Fig. 6a, the typical Mach number ¦eld, obtained in calculation for the fuel injection scheme 30%30%40%, is shown. To demonstrate details, only a part of combustor is shown in this ¦gure, and scale along the duct width is increased. This ¦eld shows that oblique shock waves arise ahead of the fuel jets, injected from the 1st injection zone pylons. At the place of their interaction with the wall, small separation of boundary layer arises. More intense shock wave is formed ahead the blunt leading edge of the central pylon. Its interaction with boundary layers on the duct walls leads to formation of separation zones of larger size. Oblique shocks, produced by these separations, intersect with leading shock wave from the central pylon in the region of jets from the 1st fuel injection zone, with lower Mach number. As a result, these shocks intersect irregularly, with formation of Mach disks, curved due to the §ow inhomogeneity in the hydrogen jets. Behind the Mach disks, there are regions of subsonic §ow. Further downstream, there are additional subsonic zones formed due to following intersections of shock waves. Figure 6b demonstrates for the same calculation the ¦eld of decimal logarithm of dimensional rate φ of heat release per unit length along a streamline: where N sp is the number of mixture components; N r is the number of reactions; W l is the molar rate of the lth reaction; v kl is the stoichiometric coe©cient of the kth component in the lth reaction; m k is the molecular weight of the kth component; h k (T ) is its static enthalpy (per unit mass); ρ is the mixture density; and � V is the magnitude of velocity vector. Derivation of formula (10) is given in [20]. Field of this parameter shows that in the jets of the 1st fuel injection zone, a weak heat release proceeds initially only in the outer regions of the jets (because of low temperature of hydrogen). Essential heat release starts only downstream from the 2nd hydrogen injection zone: growth of pressure and of temperature in shock-wave structures accelerates the reaction and decrease of velocity in this region leads to longer residence of fuel in the reaction zone. Downstream from the 2nd fuel injection zone, combustion proceeds in the regions of subsonic or transonic §ow. Curved leading shock wave ahead of the central pylon and curved Mach disks and, also, boundary layer separations generate vorticity and become strong generators of turbulence. Growth of turbulence enhances combustion downstream from the 2nd hydrogen injection zone through the transport of heat from combustion zones to cold §ow regions and through the transport of reagents to combustion zones. In the separations on the walls, there is no combustion because of the absence of fuel. Combustion practically stops at considerable distance upstream from the section, where the duct width begins to grow. When the heat release stops, the Mach number starts to grow. This e¨ect is related with growth of the duct area (height of the duct increases along x axis as it may be seen from the side view of the duct in Fig. 4); in addition, the turbulent di¨usion carries heat across the duct and diminishes the average temperature of §ow, enhancing increase of Mach number.
The best thrust characteristics have been obtained for the fuel injection scheme 30%30%40% and 33%33%34%. Analysis of §ow ¦elds obtained in 2.5D calculations has allowed to explain this result.
In Fig. 7, the static temperature ¦elds obtained in all calculations are shown. Because of low temperature of the injected hydrogen, in all schemes of fuel injection, the combustion upstream from the central pylon is possible only in outer regions of hydrogen jets, not inside these jets. Core of lateral jets remains to be too cold along the whole combustor in the scheme 50%50%0%, while the core of central jet is cold throughout in the scheme 10%10%80%. From the scheme 20%20%60%, core of central jet has enough time to start burning. But full development of combustion in the whole core of the central jet becomes possible only in the schemes 30%30%40% and 33%33%34%. In the case of scheme 30%30%40%, the most uniform distribution of heat across the duct is achieved and inner thrust of combustor is close to maximum. In the scheme 33%33%34%, heat release in the central jet is less due to lower quantity of fuel and in the scheme 40%40%20%, the heat release in lateral jets also begins to diminish because of insu©cient heating (thickness of jets increases with the growth of mass- §ow rate).  When in 2.5D calculations the most attractive regimes of §ow in combustor had been found, 3D calculations of these regimes have been performed. In Fig. 8, the block structure of computational grid for 3D calculation and ¦elds of temperature in several cross sections of combustor are shown. Figure 9 shows the ¦elds of temperature in the combustor sections by horizontal (Fig. 9a) and vertical (Figs. 9b and 9c) planes. The presented data correspond to the fuel injection scheme 30%30%40%. Figure 9Á demonstrates that it is possible to ¦nd horizontal section of the combustor, where the §ow ¦eld is qualitatively analogous to ¦elds, obtained in 2.5D calculations. But §ow structure in 3D calculation is naturally more complex. For example, it is impossible to reproduce the region of slow §ow in the wake after the pylons of the 1st injection zone. Moreover, the injected jets of hydrogen are placed at di¨erent heights. As a result, maximum of temperature that can be seen in 2.5D calculation at the duct symmetry plane is absent in the horizontal section in Fig. 9a. This maximum is reached on the outer boundary of the upper central jet and it is placed below the plane of Fig. 9Á. One can see this maximum in the combustor section by the vertical symmetry plane (see Fig. 9c).
In Fig. 3b, the distribution of static pressure along centerline of the duct upper wall, obtained in 3D calculation for the fuel injection scheme 30%30% 40%, is compared with the static pressure distribution along the centerline of computational domain in the corresponding 2.5D calculation. Agreement is quite satisfactory, especially if one takes into account that the plot from 3D calculation is obtained with local values of pressure, whereas in 2.5D calculation, the pressure is averaged along y axis.

CONCLUDING REMARKS
Parametrical study, described in section 4, would be impossible in the framework of quasi-1D theories, and 3D calculations would require too big resources of computer time and memory. In addition, physical analysis of §ow structure obtained in 2.5D calculations is much more simple than in the 3D case.
The 2.5D approach is approximate way of §ow analysis. It cannot take into account all features of 3D §ows but it provides much more information about §ow physics than the quasi-1D calculations. That is why, it may be recommended for use at the stage of preliminary design of combustion chambers. The 2.5D calculations can allow to make preliminary choice of most valuable variants of geometry and of §ow regime and to diminish essentially the quantity of expensive and prolonged 3D calculations (which are, of course, necessary at the ¦nal stage of the combustor design).