IJCPE Vol.9 No.3 (2008) Iraqi Journal of Chemical and Petroleum Engineering Vol.9 No.3 (September 2007) 37-41 ISSN: 1997-4884 COMPUTATIONAL ANALYSIS OF THE MIXING ZONE IN THE COMBUSTION CHAMBER OF RAMJET Adil A. Al-Hemiri * and Sa’ad A. Fa’ek * Chemical Engineering Department - College of Engineering - University of Baghdad – Iraq Abstract A theoretical analysis of mixing in the secondary combustion chamber of ramjet is presented. Theoretical investigations were initiated to insight into the flow field of the mixing zone of the ramjet combustor and a computer program to calculate axisymmetric, reacting and inert flow was developed. The mathematical model of the mixing zone of ramjet comprises differential equations for: continuity, momentum, stagnation enthalpy, concentration, turbulence energy and its dissipation rate. The simultaneous solution of these equations by means of a finite-difference solution algorithm yields the values of the variable at all internal grid nodes. The results showed that increasing air mass flow (0.32 to 0.64 kg/s) increases the development of velocity profile due to the high turbulence generated resulting in very fast mixing and homogenous flow. And the occurrence of chemical reaction causes higher local temperature and composition resulting in faster development of the velocity profile Keywords: oil Non-premixed chemically reacting ducted flows, Mixing of two streams, Ramjet modeling. Introduction Increased interest in ramjet propulsion system with higher performance requirements and tighter constraints on system size and weight has led to the need for improved techniques for analyzing and designing such system. A critical requirement for achieving high system performance within specified geometric limits is an accurate description of the secondary combustor flow field, including the effect of the mixing zone and chemical reaction. In a conventional application (fig.1), the mixing zone is necessary for more complete combustion before the gases exit through the nozzle. It is important to be able to determine the combustion behavior within this region as a function of upstream conditions. A recent study (Faek, 2004) gave a comprehensive survey on the subject and considered the effect of various parameters on the velocity profile. These parameters were: air mass flow rate, air to fuel ratio, working pressure, air and fuel temperatures. He showed that: Increasing air mass flow increases the development of velocity profile due to the high turbulence generated resulting in very fast mixing and homogenous flow. Increasing air and fuel temperature reduces the cooling effect of air and gives higher mixing temperature because it increases the thermal energy exchange between the air and fuel streams. Increasing upstream combustor pressure increases the development of velocity profile due to the development of highly recirculating region. The fuel to air ratio has little effect on the mixing zone characteristics, since it only effect the inlet velocities of air and fuel thus it give the same profile and the mixing temperature. The best air mass flow was reported to be 0.64kg/s, upstream combustor pressure of 8 bar, air inlet temperature of 500K, fuel inlet temperature of 1100K and fuel to air ratio of 0.125 (Faek,2004). University of Baghdad College of Engineering Iraqi Journal of Chemical and Petroleum Engineering Removal of asasbn sdfs fs f sf sdfsdfsdf fsdfsdfsdfs fsdfsdfsd 2 IJCPE Vol.9 No.3 (2008) The purpose of this paper is to study the effect of air flow rate on the mixing zone velocity profile (for the two cases of reacting and non-reacting systems) by using a governing equations and solving them by a computer program developed for this goal. MODEL OVERVIEW The flow was assumed to be steady, two-dimensional, and subsonic. For simplicity, the value of specific heat (CP) was assumed to be constant although its dependence on temperature and/or composition could be easily included. A modified Jones-Launder,( Jones and Launder,1972; Launder and Spalding,1974) two- parameter turbulence model was incorporated to calculate the effective viscosity. It uses five empirical constants (Table 1) and requires that two additional variables [turbulence kinetic energy (k) and turbulence dissipation rate ()] be evaluated. Effective viscosity was calculated using the formulas μμμ teff  (1) Where, ερ kCμ 2 μt  (2) For reacting flows, four species were considered: oxygen, nitrogen, fuel, and products. Simple, one step, infinitely fast kinetics was assumed in which 1kg of fuel combines with i kg of oxidant to form (1+ i) kg of product without intermediaries (Moult and Srivatsa, 1977). Fuel and oxygen, therefore, could not exist simultaneously, and the combustion process was mixing limited. The turbulent Prandtl and Schmidt numbers were taken equal to unity, and therefore, the turbulent Lewis number was unity; the laminar Prandtl number was also taken to be unity. The conservation equations for axisymmetrical flows with no tangential variations can be put into general form (Khalil, E.E., Spalding and Whitelaw,1975).        S r r rrxx Vr rr U x                             11 (3) where ψ stands for the dependent variables (U,V,k,,h, etc….) being considered (ψ=1 for continuity equation), ψ is the appropriate effective exchange coefficient for turbulent flow, and Sψ as the "source term" (Table 2). The energy equation in terms of stagnation enthalpy has no source terms since the turbulent Prandtl and Schmidt numbers were chosen equal to unity and radiative transport was neglected. Thus, the stagnation enthalpy is given by   2VUhh 22  (4) where for non-reacting flows TCh p (5) and for reacting flows  refPox TTCiΔΗmh  (6) Table (1) Turbulence model constants (Launder and Spalding, 1974) Constant Value C1 1.44 C2 1.92 Cμ 0.09 K 0.42 E 8.8   22.1CCCKσ 9.0σσσ 2 1 μ12 2 ε fhk          IJCPE Vol.9 No.3 (2008) Table (2) Governing equations parameters (Launder and Spalding, 1974) Conse rvatio n of ψ Γψ Sψ Mass 1 0 0 Axial mome ntum U  eff                            x V r rrx P x U x effeff  1 Radia l mome ntum V  eff 2 2 1 r V r U xr P r V r rr effeffeff                            Kineti c energ y k   k eff G Dissi pation rate     eff    21 CGC k  Stagn ation enthal py h   h eff 0 Mass fracti on of fuel m i   fu eff 0 Mixtu re fracti on f f eff   0 Special Treatment near A Wall To avoid the need for detailed calculations in the near- wall grid regions, algebraic relations are employed for the near wall grid nodes, which have to be spaced at such a distance from the neighboring walls that they lie within the so-called logarithmic layer. Such relations are termed wall functions which are derived so as to reproduce identically the implications of the logarithmic profiles with uniform shear stress prevailing up to the near wall grid node, with generation and dissipation of energy are in balance at this locality (Patankar and Spalding,1970). The particular ones employed for the equation of the three velocity components, k and , are:      ν/ykECln κ 1 ρτ kC U P 21 P 41 μ w 21 P 41 μ P  (7)   21 μwP Cρτk   (8) P 23 P 43 μ P yκ kC ε  (9) The Solution Procedure The set of differential equations are first reduced to finite-difference equations exhibiting "hybrid formulation of the coefficients", i.e., coefficient that contain combinations of the convective flux per unit mass F and the diffusive conductance. And then solved iteratively by "SIMPLE" procedure (Patankar and Spalding, 1972; Patankar, 1980; Versteeg and Malalasekera, 1995). The grid arrangement for discretisation is as shown in figure (2). And the logic diagram for SIMPLE algorithm is given by figure (3), where the procedure steps are as given below: - Start - Initial guess of p*, u*, v*, Φ*. - STEP! Solving discretised momentum equations, ai,j U*I,j = Σ anb U*nb + ( P*I-1,J - P*I,J ) Ai,J + bi,J aiIj V*I,j = Σ anb V*nb + ( P*I,J-1 - P*I,J ) AI,j + bI,j - STEP 2: Solving pressure correction equation, aI,J P'I,J = aI+1,J P'I+1,J + aI-1,J P'I-1,J + aI,J+1 P'I,J+1 + aI,J-1 P'I,J-1 + bI,J - STEP 3: Correction of pressure and velocities, PI,J = P*I,J + P'I,J UI,j = U*I,j + di,J (P'I-1,J - P'I,J ) VI,j = V*I,j + dI,j (P'I-1,J - P'I,J ) - STEP 4: Solving the discretised equations for ρ and T, Removal of asasbn sdfs fs f sf sdfsdfsdf fsdfsdfsdfs fsdfsdfsd 4 IJCPE Vol.9 No.3 (2008) aI,J ρ I,J = aI+1,J ρ I+1,J + aI-1,J ρ I-1,J + aI,J+1 ρ I,J+1 + aI,J-1 ρ I,J-1 + bI,J aI,J TI,J = aI+1,J TI+1,J + aI-1,J TI-1,J + aI,J+1 TI,J+1 + aI,J-1 TI,J-1 + bI,J - Test for convergence:, If yes STOP program, If no SET: P* = P; U* = U; V* = V; ρ* = ρ; T* = T. - Return to STEP 1 Results and Discussion The effect of air mass flow rate on the combustor mixing zone was studied by using the computer program, under two conditions (inert and reacting flows), in which pressure values used was (8 bar) and inlet air mass flow rate are (0.32, 0.64 kg/s) inlet air temperature was about (500K), fuel inlet temperature was (1100K) and fuel to air ratio was (0.125). The results are shown in vector forms in figures below, where the radial (r) and axial (x) distances are in meter units. From Fig. (4) it can be seen that the deflection of velocity is from the lower velocity stream to the higher velocity stream, i.e., from air stream to fuel stream where the fuel stream behaves as the main stream. In Fig. (5), the fuel stream is deflected to air stream that give a recirculation region downstream of fuel and air inlets and the velocity profile has its nearly constant shape after the end of the recirculation region. This, also, leads to low local values of temperature profile due to the large local quantity of the air stream and because there is no chemical reaction. From this comparison it can see the effect of inlet air mass flow rate on the ramjet mixing zone flow properties at inert flow system and at constant other inlet flow conditions. Considering Fig. (6), it can be seen that the air stream is deflected to fuel stream and the fuel stream behaved as the main stream and because there is no recirculation region the velocity profile developed faster and began to have its constant shape. In Fig. (7) it can be seen that the fuel stream is deflected to air stream and air stream behaved as the main stream in the flow field and recirculation region occurred down stream of the fuel inlet and the velocity profile begin to have its nearly constant shape at the end of the recirculation region. The chemical reaction gives a higher local temperature (since combustion is exothermic reaction) and a higher local axial velocity (due to thermal expansion of the product gases) than the inert flow case. It can be seen from Fig. (4) and Fig. (6) that the velocity profile is the same but the flow conditions are different because in Fig.(4) the system is inert and there is no changing in the density of the mixing zone but for Fig. (6) the reaction occurring give a high change in the temperature and composition that a homogenous distribution of the product makes the mixing zone length shorter. Similar effect may be deduced by comparing Fig. (5) and (7). Fig.1: Flow boundaries with symmetry for a ramjet. Fig.2: Grid arrangement N,n: North; S,s: south; W,w:west;E,e:east Top-wall boundary air fue l S id e – w a ll b o u n d a ry O u tle t b o u n d a ry Symmetry line boundary IJCPE Vol.9 No.3 (2008) Step 1: START Step 2 Step 3 Step 4 Convergence ? STOP Yes No p=p*, u=u*, v=v* ρ,= ρ*, T*=T Figure (3); Logic Diagram For Program “SIMPLE” Initial guess: p*,u*,v*,ρ*, T* ρ,T p, u, v, ρ*,T* p* u*, v* Fig.4: velocity vector plot at =0.125, ma =0.64 kg/s, pc=8bar, Tf=1100k and Ta= 500k, and inert flow system. Fig.5:velocity vector plot at =0.125, ma =0.32 kg/s, pc=8bar, Tf=1100k and Ta= 500k, and inert flow system. Fig.6: velocity vector plot at =0.125, ma =0.64 kg/s, pc=8bar, Tf=1100k, Ta= 500k, and reacting flow system. (2D) ½ 23 Sep 2004 ½ 0 0.1 0.2 0.3 x 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 r (2D) ½ 23 Sep 2004 ½ (2D) ½ 22 Sep 2004 ½ 0 0.1 0.2 x 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 r (2D) ½ 22 Sep 2004 ½ (2D) ½ 22 Sep 2004 ½ 0 0.05 0.1 0.15 0.2 0.25 0.3 x 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 r (2D) ½ 22 Sep 2004 ½ Removal of asasbn sdfs fs f sf sdfsdfsdf fsdfsdfsdfs fsdfsdfsd 6 IJCPE Vol.9 No.3 (2008) Figure 7 velocity vector plot at =0.125, ma =0.32 kg/s, pc=8bar, Tf=1100k, Ta= 500k, and reacting flow system. CONCLUSIONS A computer program has been developed for the analysis of ramjet mixing zone combustor. And comparison between the results at different input has shown the effect of the parameters selected in this study on the mixing zone characteristics of ramjet combustor. Reducing air mass flow rate will reduce the cooling effect of air stream and give higher mixing temperature. And increasing it will cause the velocity profile to develop faster. The difference between the performance with inert or reacting systems is that for the reacting system, the reaction causes changes in temperature and composition leading to a homogenous distribution and hence a faster development of velocity profile, i.e., shorter length of mixing zone. NOMENCLATURE C1,C2, Cμ k-ε model empirical constants (Table 1) CP Specific heat at constant pressure. E Integration constant in wall function=9 G Generation term for kinetic energy of turbulence. h Stagnation enthalpy h  enthalpy k Turbulence kinetic energy kP Wall function turbulence kinetic energy ma Inlet air mass flow rate pc Up stream combustor pressure r Radial distance S Source term of the general governing equation Ta Inlet air temperature Tf Inlet fuel temperature U Axial velocity UP Wall function axial velocity V Radial velocity x Axial distance yp Wall function distance variable GREEK LETTERS Γ Effective diffusion coefficient  Dissipation rate of turbulence kinetic energy  p Wall function dissipation rate К von Karman constant  Air to fuel ratio μ Viscosity μ t Turbulent viscosity μ eff Effective viscosity ψ General dependent variable ρ Density σ Prandtl or Schmidt number τ Shear stress REFERENCE 1. Faek, S. A., "Analysis of the Mixing Zone in the Combustion Chamber of Ramjet" M. Sc thesis, Baghdad University, 2004. 2. Jones, W.P., and Launder, B.E., "The Prediction of Laminarization with a Two-Equation Model of Turbulence", International Journal of Heat and Mass Transfer, Vol.15, pp.301-314, Feb. 1972. 3. Khalil, E.E., Spalding, D.B., and Whitelaw, J.H., "The Calculation of Local Flow Properties in Two Dimensional Furnaces", International Journal of Heat and Mass Transfer, Vol. 18, pp. 775-791, 1975. 4. Launder, B.E., and Spalding, D.B., "The Numerical Calculation of Turbulent flows", Computer Methods in Applied Mechanics and Engineering, Vol.3, No.2, pp. 269-289, 1974. 5. Moult, A. and Srivatsa, S.K., "A Computer Code for axi-symmetric Combustion Chambers", CHAM, 1977. 6. Patankar, S.V., and Spalding, D.B., "Heat and Mass Transfer in Boundary Layers", 2nd edition, London, 1970. 7. Patankar, S.V., and Spalding, D.B.,"A Calculation Procedure for Heat, Mass, and Momentum Transfer in Three-Dimensional Parabolic Flows", International Journal of Heat and Mass Transfer, Vol.15, pp.1787- 1806, 1972. 8. Patankar, S,V., "Numerical Heat Transfer and Fluid flow", Hemisphere Publishing Corporation, Taylor & Francis Group, New York, 1980. 9. Versteeg, H.K., and Malalasekera, W., "An Introduction to Computational Fluid Dynamics: The Finite-Volume Method", Longman Group Ltd, 1995. (2D) ½ 22 Sep 2004 ½ 0 0.1 0.2 x 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0.11 r (2D) ½ 22 Sep 2004 ½