Subscriber access provided by Nottingham Trent University
Combustion
Numerical investigation of gaseous hydrogen and liquid oxygen combustion under subcritical condition Amir Mardani, Arash Ghasempour Farsani, and Mohammad Farshchi Energy Fuels, Just Accepted Manuscript • DOI: 10.1021/ acs.energyfuels.9b02050 • Publication Date (Web): 26 Aug 2019 Downloaded from pubs.acs.org on August 26, 2019
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
Numerical investigation of gaseous hydrogen and liquid oxygen combustion under subcritical condition Amir Mardania,Arash Ghasempour Farsania, Mohammad Farshchia a
Sharif university of Technology, Azadi Ave., Tehran, Iran, P.O. Box: 11365-11155 Coresponding Athour:
[email protected] Abstract
This study is on combustion modelling of gaseous hydrogen and cryogenic liquid oxygen at subcritical condition for well-known Mascotte laboratory combustor. The proposed strategy relies on hybrid EulerianLagrangian framework in which the continuous phase is evaluated by Reynolds Average Navier-Stokes (RANS) equations and quick discretization method. The dispersed phase of combustion field is evaluated by Discrete Phase Method (DPM). The Eddy Dissipation Concept (EDC) has been performed for combustion-turbulence interaction modelling. Effects of turbulence model, chemical kinetic mechanism, equation of state, and inlet momentum jet flux are investigated in terms of combustion field characteristics. The accuracy of predictions is identified through a detailed comparison with the experimental databases gathered on the ONERA Mascotte test bench. Results show that predictions using the turbulence model family, consists of Standard and RNG models, are closer to the measurements. In order to investigate effect of compressibility on combustion flow field, compressible and incompressible forms of ideal gas equation of state have been utilized. Mach number and compressibility factor values in vicinity of the injector post emphasize to apply compressible form of ideal gas. Four reduced chemical kinetic mechanisms are performed to compute chemical properties. It is found that the result of all mechanisms are highly similar, however Burke mechanism can be selected as the best. Moreover, C10 test case has been modeled to study the effect of mixture ratio change. WSR modelling shows no differences in H2 oxidation neither routs, nor reaction rates; and all differences in results refer to less H2 / O2 momentum flux and unburnt fuel for C10. Keywords: Turbulent Combustion; Multiphase Flow; Cryogenic temperature; Subcritical Condition; Liquid Oxygen
1.
Introduction
Over the last decades, liquid propellant engines have been used as a main part of space vehicles to explore the new horizons of space. In this way, hydrogen propellant with its desirable performance has been widely assigned in spacecraft engines. Consideration of some major shortcomings include complexity of phenomena and development costs conducts to progress analytical and Computational Fluid Dynamic (CFD) design tools beside of experimental approaches. Nowadays, the critical parameter is how to build such system which gains lower budget and efficient performance in a short time interval. From thermodynamic point of view, complexity of cryogenic combustion necessitates to review the phenomena. As it is known, the operation condition of combustor imposes particular effects especially on spray and droplet size distribution. Injection of cryogenic liquid oxygen in combustor with pressure higher than its critical value is the main process that happens at trans- and supercritical conditions [1]. In these situations, the fluid behaves gas-like in terms of density and liquid-like in terms of transport properties [2-4]. Consequently, such special behavior causes uncommon variations in transport properties, density, surface tension, and evaporation enthalpy [5, 6]. But under lower pressures, there is classical standard behavior associated with spray, especially breakup and atomization processes. It means that there is an interaction between inertia, viscous, and liquid surface tension forces which leads to dissociation of liquid jet for the cryogenic Lox / GH2 injection at subcritical condition in the vicinity of injector plate. Research community has attempted to study multiphase flows in different manners of numerical modelling. There are two major approaches to study multiphase flow field. One of them is Eulerian-Lagrangian approach which is based on calculation of tracking of each droplet trajectory in a particular reference frame [7-10] . The trajectory is obtained by setting the equation of force balance for that specific droplet in every arbitrary moment [11, 12]. The other approach evaluates both gaseous and dispersed phases in segregated continuum frameworks and obtains the specifications of dispersed phase in average quantities [13]. With respect to the type of
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
2
algorithm, the dispersed phase specifications are averaged in a volumetric or ensemble definition, and evaluated through the Probability Density Function (PDF) [14, 15]. Very limited facilities study cryogenic fluid flows under various conditions. In order to understand the processes involved in cryogenic combustion, Office National d'Etudes et de Recherches Aérospatiales (ONERA) has developed Mascotte facility [16, 17] at Palaiseau center. Since 1991, three versions of combustors have been built to achieve higher pressures and mass flow rates. Accessible chamber pressure up to 10Mpa, Lox mass flow rate range of 100gr/s to 400gr/s, sufficient test time duration, also optical and laser-based diagnosis, allow research groups to investigate processes of cryogenic combustion under sub/trans/supercritical conditions. P8 [18, 19] combustor is a high pressure test facility which was invested by French and German partners, located at Deutsches Zentrum für Luft- und Raumfahrt (DLR) research center of Lampolshuasen, since 1995. The interface pressure range of 2.5-36Mpa, GH 2 mass flow rate of 0.5-1.5kg/s, and Lox mass flow rate of 0.2-8.0kg/s provide diverse conditions to investigate in P8 combustor. Another test facility at DLR Lampoldshuasen is M3 combustor, [19, 20] used to investigate atomization and combustion of cryogenic hydrogen and oxygen mixture at pressure range of 0.1-2Mpa. Cryogenic Combustion Laboratory (CCL) [21-23] is located in the Propulsion Engineering Research Center (PERC) at the Pennsylvania State University which was funded by National Aeronautics and Space Administration (NASA) at 1990. Data is generated by laser-based and optical diagnostics in relatively wide pressure range of 200-800 psi, and oxidizer to fuel proportion range (O/F) of 0.5-170 for several types of injectors [24]. Facility of University of Florida [25] provides chamber pressure up to 60atm and O/F proportion range of 0.44-8.25 for different type of shear coaxial injectors with accessible optical diagnostics. Mascotte is the best choice from the above list due to significant experimental data, quantitatively and qualitatively, which it operates at subcritical condition with A10 operation condition. There are several numerical studies which have investigated phenomenology of cryogenic combustion in ONERA Mascotte under subcritical condition with respect to multiphase frameworks. In some, a single Eulerian approach has been performed to evaluate multiphase flow field. Chen et al. [26] used ONERA Mascotte combustor data to validate their computational method which it was performed to model the combustion process of a hybrid system based on propellant mixture of N2O-HTPB. In some other studies, two Eulerian descriptions have been utilized. Farmer and Cheng [27, 28] modeled injection of liquid oxygen and gaseous hydrogen by two Eulerian frameworks in terms of Volume of Fluid (VOF) model [29] which it simulates primary breakup and atomization processes to capture interface between gas and liquid phases. By applying a distinct transport equation for density of liquid, Jay et al. [30] investigated the two phase Lox / GH2 injection flow field and study atomization and breakup numerically. Bodèle et al. [31] used such an approach to simulate the effects of secondary breakup truncation. Some other study use Eulerian-Lagrangian approach for investigation. Farmer and Cheng [27, 28] modeled two phases reacting flow field of RCM-2 test case which it tracks cluster of droplets based on International Workbench Rocket Combustion Modeling (IWRCM) [32] recommendations. Pourouchottamane et al. [33] employed a Lagrangian solver to track droplets with one/two way coupling feature. Izard and Mura [34] characterized the dispersed phase by performing a PDF evaluated by Williams spray equation [35]. Comparison of steady and unsteady Reynold Average Navier-Stokes (RANS) approaches was studied by Lemke et al. [24] on PennState preburner combustor configuration alongside the ONERA Mascotte. Gomet et al. [36] considered additional source terms in transport equation of mixture fraction induced by evaporation of liquid phase . Moreover, there are some CFD simulations which have been performed Large Eddy Simulation (LES) [37-40] approach to evaluate turbulent closure, and none of them has taken account second phase to model. Furthermore, a few Direct Numerical Simulation (DNS) [41-43] approaches are utilized to model the zone near the injector. Table 1 shows previous modeling strategies which were used in some important numerical investigations of Mascotte A10. With a quick glance, it can be understood that there is limited list of numerical approaches which have been employed. It should be considered that most studies were mainly performed two-equation models, especially Standard, to evaluate turbulent closure. Investigation of different turbulent models performance could also be contemplated. To considering details of primary atomization process, experimental correlations recommended by IWRCM [32] has been used in entire simulation. Because of spray quality and the quantity of volume fraction (which has been explained in section 3.2), most of studies neglected the secondary breakup process. Regarding the chemical kinetic mechanisms, some more detailed and comprehensive ones could still be used to represent the reaction more precisely. Moreover, to predict turbulent combustion process precisely, turbulent-
ACS Paragon Plus Environment
Page 2 of 31
Page 3 of 31
Energy & Fuels
1 2 chemistry interaction is an open area of research which could be studied in more, because are less 3 addressed in modelling of Lox / GH2 combustion. Indeed, most of developing numerical investigations of 4 subcritical combustion concentrate on high precise approaches with massive resolution in computational 5 domain without considering to budget timely and costly. In this way, it is worthy to establish a modeling 6 strategy with logical budget, and is still a challenge in the process of parametric study. 7 8 Table 1 Sub-models of some numerical studies, used to model multiphase combustion of Mascotte A-10 test case 9 10 Sub-models 11 Study Primary Secondary Multiphase Multiphase Continuous Turbulent Dispersed Turbulent Equation of Chemical Reaction Atomization Atomization 12 Framework Model Model Model State Model Model Model Eulerian13 HBMS Finite rate Cheng and Farmer[27, 44] Eulerian VOF[45] Cloud Model[12, 46] ignored TAB[47] equation[48, 49] Chemistry model description 14 EulerianTwo equations HBMS Eulerian 15 Chen et al.[26] model equation description Eulerian16 Kelvin Jay et al.[30] Eulerian TAB Standard Helmhotz description 17 EulerianDeveloped IWRCM EBU combustion 18 Bodèle et al.[31] Eulerian experimental Standard Guideline model description correlation[31] 19 EulerianDroplet track Eddy Life Time IWRCM Pourouchottamane et al. [33] Lagrangian ignored Perfect gas CLE model model (DLS) dispersion model Guideline 20 description EulerianPDF based on PDF based on IWRCM 21Izard and Mura [29] Lagrangian Stochastic Particle ignored PDF model Williams equation Standard Guideline description method 22 EulerianParticle Tracking IWRCM Cubic equation Finite rate Lagrangian ignored ignored q model Guideline of state Chemistry model 23 Lemke et al. [24] discription EulerianPDF Based on PDF based on IWRCM Compressible 24 Gomet et al. [36] Lagrangian Stochastic Particle ignored PDF model Williams equation Standard Guideline Model description method[52] 25 EulerianParcels Tracking Wilcox IWRCM Cubic equation Finite rate Lagrangian ignored model Guideline of state Chemistry model 26 Lemke et al.[53] description 27 In this paper, several numerical algorithms and computational models have been used to improve diversity 28 of significant sub-models of previous studies [24, 26-28, 30, 31, 33, 34, 36, 44, 53] which could been 29 performed to model subcritical cryogenic combustion of Mascotte A-10. For this purpose, following road map has been traced: first, the geometry, initial, and boundary conditions of Mascotte test facility are 30 presented to create an optimum computational domain. At second, numerical methods to evaluate 31 properties of multiphase reacting flow are described in section 3. For the sake of comparison, several 32 substitutions have been considered in sub-models. To evaluation of turbulent closure, two-equation models 33 of Standard, RNG, SST, and four-equation model of Transition SST have been 34 employed. To study the effects of compressibility, ideal gas equation of state has been examined in two 35 compressible and incompressible forms. Moreover, four different chemical mechanisms have been utilized 36 to analyze effect of chemical kinetics on combustion flow field properties. Also, to comprehend the effect 37 of various momentum flux of injector, two different operation conditions of A10 and C10 have been 38 investigated. Finally, the results of performed numerical strategies are explained in last section. 39 40 2. Experimental Setup 41 42 Mascotte, as a pioneer test facility, provides a wide range from subcritical to supercritical operation 43 conditions for experimental investigation of cryogenic combustion [16, 54]. In this study, the focus is on the 44 A10 which it was proposed to investigate phenomenology of subcritical non-premixed Lox / GH2 45 combustion. The test case of A10 corresponds to the case of RCM-2 represented by IWRCM [32]. Also 46 for parametric study of mixture ratio, the case of C10 has been modelled and compared with A10. Detailed 47 conditions of A10 and C10 have been reviewed in Table 2. 48 49 Table 2 Mascotte A10 and C10 operation conditions 50 Operation conditions A-10 C-10 51 H2 O2 H2 O2 52 Mixture ratio 2.1 3.2 53 Momentum Flux 14.5 6.5 54 55 ACS Paragon Plus Environment 56 57 58
3
Chemical Mechanism 9 step reactions[50, 51] 9 steps reactions 1 reaction 19 steps reaction 19 steps reaction
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 4 of 31
4
Pressure [MPa] Mass flow rate [gr/s] Temperature [K] density [kg/m3]
1.05 23.7 287 0.84
1.05 50 85 1170
0.95 15.8 287 0.84
0.95 50 85 1170
As seen in Figure 1, Mascotte combustor consists of a single hollow post injector in a chamber with square cross section, and operates at pressure of 1 MPa for A10 and C10 operation conditions. For the reason that the critical pressure of oxygen is 5.04 MPa, accordingly the test cases is under subcritical condition. The internal orifice for injection of Lox with diameter of 3.2 mm is surrounded by GH 2 tube with inner diameter of 5.6 mm and outer diameter of 12 mm. As mentioned, the chamber has a square section with internal dimensions of 50 50 square millimeter and the distance from injector plate to outlet is equal to 480 mm. To gathering experimental data, two silicon windows have been provided close to the injector plate. These experimental data are obtained in significant variety by specific optical diagnostic methods [2, 16, 17, 55, 56]. Phase Doppler Particle Analyzer (PDPA) method has been utilized to visualize droplet breakup and atomization in terms of particle size and velocity. Coherent Anti-Stokes Raman Spectroscopy (CARS) has been used to capture local temperature of continuous phase. As another observation diagnosis, Laser Induced Fluorescence (LIF) has been employed to visualize hydroxyl spectrometric image in Mascotte facility. In this study, some of these experimental data are used to validate numerical modelling.
Figure 1 Geometrical details and boundary conditions of Mascotte combustor
3. Numerical Method Figure 2 illustrates all geometrical characteristics of A10 test case used to create the computational domain. The geometry of Mascotte combustor has a cross section of 50 50 mm . Nevertheless, the assumption of two dimensional axisymmetric domain is considered to simplify and gain lower calculation cost, such as the similar simulations [24, 26-28, 30, 31, 33, 34, 36, 44, 53]. In this manner, the RANS approach has been followed for modelling the continuous phase. Conservation equations of mass, momentum, energy, and species as well as turbulent equations have been solved. Two-equation models of Standard, RNG, SST, and four-equation model of Transition SST have been implemented separately to study on the effect of turbulence model. It seems that incompressible form of ideal gas, as the equation of state, is able to predict thermodynamic properties sufficiently, because of low reduced pressure and high temperature conditions in chamber atmosphere. But to investigate the effect of compressibility, compressible form of ideal gas equation has been examined too. Finite rate chemistry approach is utilized for the chemical reactions modelling of combustion flow field. Exclusively the net rate of production of each species are obtained by turbulence-chemistry interaction model called Eddy Dissipation Concept (EDC) [57] with assumption of reaction occurring in very small scale of turbulent structures [58]. Four multi-step chemical kinetic mechanisms of Burke et al. [59], O’conaire et al. [60], Konnov [61], and Alekseev et al. [62] have been used. Pressure-based algorithm with Simple solution scheme of pressure-velocity coupling is utilized with discretization of weighted average of upwind second order all over the modeling. Also to evaluate mass and heat diffusivity, kinetic theory has been used [63]. Student version of Ansys-Fluent 19.1 has been adopted to perform all mentioned models and schemes. As a suitable procedure in this manuscript, the primary atomization has not been considered to model. There are some recommendations from IWRCM about modelling of Mascotte test case such as hydraulic 2
ACS Paragon Plus Environment
Page 5 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
5
diameter and turbulent intensity which are used in this study [32]. Furthermore, the guideline about the size and velocity of droplets right after injection surface is also highly applicable. Presuming a conical surface for the core of Lox where a specific distribution of droplets spread from, avoid us to struggle with primary breakup and atomization modelling. Based on the momentum ratio and an empirical correlation, the length of conical surface along the injector axis equals to 7.8 mm [64]. Experimental data gathered from PDPA establishes a form of Rosin-Ramler distribution for particles scattered from the presumptive conical surface as follows: D m F ( D ) 1 exp (1) Where the two empirical constants of and m are 130 and 2.25, respectively. Also with respect to experimental data, the Sauter Mean Diameter (SMD) or DSMD is 82 micro meter [65]. The liquid oxygen particles are scattered with constant velocity component of 10 m/s from the presumptive conical surface in various angles depending on radial position (see for instance Figure 2) as follows: r tan ( r ) r R (2) Lc ( 1 ) i Ri tan 0 Where Lc is length of Lox core and equals to 7.8 mm, R i is the radius of Lox injector, and tan0 is u defined as 0.68 gas 1 gas , based on an empirical correlation of velocities and densities of Lox and u liq liq GH 2 . Full-Structured Computational Domain Pressure Farfield Condtion (Atmospheric) Wall Condition (Adiabatic, No Slip)
00
0.05 0.05
0.1
Inlet Condition
0.1
0.15
0.15
0.2
0.2
0.25
0.25
0.3
0.3
0.35
0.35
0.4
0.4 0.45
0.45 0.5
AxialSymmetry SymmetryCondition Condition Axial
Figure 2 geometry and boundary conditions of computational domain (top) and injector (bottom).
Additionally, the secondary atomization has been neglected. the parameter of volume fraction, , classifies the method to setup secondary breakup [12] modelling. A range of 103 indicates dilute spray which
ACS Paragon Plus Environment
Energy & Fuels
6
particles interaction and collisions can be neglected. For collision-dominated and contact-dominated sprays, the range of volume fraction are 103 101 and 101 , respectively. In section 4.1 and based on such criterion, the variation of volume fraction approves the assumption of ignoring secondary breakup (see for instance Figure 10). The Abramzon-Sirignano model [66] has been performed to evaluate droplet vaporization. Based on Lagrangian approach, the code predicts trajectory of an arbitrary droplet by integrating force balance which is written in an exclusive reference coordinate located on that droplet. Inertia of droplet and other forces participate in the equation of force balance. Also, three heat transfer and mass diffusion relations have been considered based on integration of energy balance to evaluate heating and mass exchange. In this study, prediction of droplet dispersing in the turbulent flow field is attained by performing assumption of equality of eddy life time and Lagrangian integral time scale. Eddy life time model was developed in concern with the turbulence models [67] such as all types of and , as follows: T 0.15 For mod els u p ( t )u p ( t ) T d (3) 2 0 T 0.15 1 For mod els u p
3.1 Modeling Verification Three different full structured quadratic grids of 55,000, 80,000 and 110,000 elements were generated and grid study has been done. Figure 3 illustrates profiles of temperature and mass fractions of OH, O2 , and H 2O for three different computational domains. As it can be seen, prediction of variables for grid with 80,000 elements has a reasonable agreement with high quality grid of 110,000 elements.
2
4
3000
Grid 1 (55000) Grid 2 (80000) Grid 3 (110000)
5
0.2
Axial temperature
2000
Radial temperature
0.1
Temperature [K]
3
0.4
r/DLOX
0.3
1
Mass Fraction
0
1000
Radial OH mass fraction Axial OH mass fraction
0
20
40
x/DLOX
60
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 6 of 31
80
Figure 3 Grid dependency study of three candidate computational domains in axial section of r / DLOX 0 and radial section of x / DLOX 1
Therefore, the optimum grid is made of approximately 80,000 elements and selected for all modelling due to the computational cost. This computational domain has been utilized to model different cases. In this study, eleven different test cases have been set to investigate effect of turbulence models, compressibility, chemical kinetic mechanisms, and mixture ratio which are described in Table 3. The results of modelling have been compared with available experimental data [55] to assess level of agreement. Figure 4 shows the variation of temperature in two different axial and radial sections. In addition to experimental data which are shown with black square, values evaluated by other numerical studies [24, 27, 28, 31, 33, 34, 36, 44, 53] are also shown with different symbols. It is obvious that the
ACS Paragon Plus Environment
Page 7 of 31
7
3000
3000
precision of prediction has been improved through current modelling. The other investigations have good agreement in radial section with experimental data, nevertheless inappropriate error can be seen in the axial section and vice versa. For example, the most agreement in axial section belongs to Cheng et al. [44] with Mean Absolute Percentage Error (MAPE) of 9.1%. However, in radial section, its MAPE is equal to 114.3%. Farmer et al. [27] has the MAPE of 11.0% as the best at the radial section, but in the axial section, MAPE is 20.7%. For present study, the MAPE at the axial and radial sections are 3.91% and 17.0%, respectively. Distribution of OH mass fraction is illustrated alongside Abel transformed of OH* emission in Figure 7. With respect to such resemblances, it can be concluded that the modelling procedure has significant performance.
40
x/DLOX
60
80
2000
# # #
#
#
*
*
*
3
4
1000
Temperature [K] 20
0
2000
Temperature [K]
1000 0
Present study EXP [55] Lempke et al. [24] Farmer et al. [27] Cheng et al. [28] Bodèle et al. [31] Pourouchottamane et al. [33] Izard et al. [34] Gomet et al. [36] Cheng et al. [44] Lempke et al. [53]
*
*
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
0
1
2
r/DLOX
5
Figure 4 Comparison of computed axial and radial profiles of temperature through current work and other reports [24, 27, 28, 31, 33, 34, 36, 44, 53] with CARS measurements from [55].
Figure 5 Comparison of OH mass fraction distribution obtained from numerical modelling of cryogenic combustion of Mascotte A10 (Case no. 1 Table 3) (top) and Abel transformed OH* emission image [68] (bottom) Table 3 Numerical test cases Turbulent model
no.
Compressibility
Std. *
RNG
SST
*
standard
†
Transition SST Incompressible ideal gas
‡
Trn. SST †
I.I.G. ‡
Ideal gas
Chemical kinetic mechanism Nonreactive
Global
Burke[59]
ACS Paragon Plus Environment
O’conaire[60]
Konnov[61]
Mixture ratio Alekseev[62]
A10
C10
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 8 of 31
8 1 2 3 4 5 6 7 8 9 10 11
*
* * * *
* * * * * * * * * *
* * * * *
* * * * * * *
* * * * *
* * * * * * * * * *
*
4. Numerical Results In this study, results of numerical modelling are described in five different sections. First, some major physics related to combustion field and flame structure are explained. Second, the effect of compressibility on the flow field are checked by performing suitable equation of state. Third, the structure of flame and its variations are investigated for different turbulence models. Fourth, the effects of using four different chemical kinetic mechanisms on modelling are studied. Finally, difference in combustion flow field of A10 and C10 test cases are investigated at a constant chamber pressure. 4.1 Structure of Flame and Combustion Field In this part, the focus is on explanation of combustion field properties from dynamics, thermodynamics, and chemical points of view. It should be mentioned that the all results used in this section are obtained from case no.1 in Table 3 which has the best agreement with the experimental data. For case no.1, the Standard; Burke et al. [59], and compressible format of the ideal gas equation of state are performed. Detail of H2 / O2 mechanism introduced by Burke et al. is presented in Table 4. To analyse the reacting flow structure, the flow field is classified to four major zones as titled in Figure 6. The figure shows that in vicinity of the injector post, anchor vortex stabilizes the flame near the lip and shapes the spreading angle in very early position. In this zone, the pressure is in the maximum value relative to the other parts of chamber due to circulation of high velocity hydrogen stream over the lip. By trapping hydrogen, the anchor vortex creates a situation to mix hydrogen and oxygen right after the lip. The very early ignition reactions occur in the region of the fuel-rich mixture, and the initial part of flame front forms right in this spot. Temperature increment in turbulent mixing layer increases the net rate of evaporation of low speed liquid oxygen excessively at any radial section, especially 0 r / DLOX 2 . Therefore, an expansion occurs, and causes to deviate hydrogen jet upward to the chamber wall. Impact of the deflected hydrogen jet to the chamber wall causes to create a high pressure stagnation zone. Figure 7 shows that the intersection location of number 1 Streamline and the wall is the stagnation point position. Part of impacted hydrogen stream continues to its previous path towards the outlet; but part of it circulates and turns back where creates a zone with negative velocity component. Existence of a region with zero axial and radial velocity causes to create a special vortex, and therefore recirculation zone. The center of the vortex is identified by intersection of zero axial and radial velocity isolines, as it is illustrated with number 2 streamline in Figure 7. However the most important effect of recirculation zone on combustion flow field is to configure flame shape in axial distance of 6 x / DLOX 12.5 . The effect of vortex of recirculation zone causes to generate a zone with minus axial velocity called backward zone. Therefore, the suction effect of recirculation zone shape a bulge in flame front and stabilizes middle zone of flame, as it could be seen in Figure 6 and Figure 7. As seen in Figure 7, drastic enhancement occurs in thickness and intensity of OH mass fraction due to reverse of gaseous oxygen in backward zone. Number 3 streamline shows how the gaseous oxygen turns back under effect of backward zone. Limited area of backward zone is controlled by two parameters. First is vorticity magnitude and suction effect of recirculation zone that specify the rear side of backward zone. Second is proximity of front zero radial velocity line to backward zone border leads to impose positive velocity gradient to backward stream. Two streamlines of number 4 &5 are tangent to both sides of backward zone and show zero axial velocity. In further axial locations, only the consumption and production rate of gaseous oxygen configures the shape of flame front. To understand level of agreement of
ACS Paragon Plus Environment
*
Page 9 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
9
numerical results with experimental data, Figure 9 illustrates superimposed streamlines on Abel transformed of OH* image. Due to distinct behaviour of fluid in range of subcritical to supercritical pressure, different numerical methods should be used to investigate each phenomenon. But different phenomenon and different numerical methods might be resembling in some physical concepts. Experimental [69] and numerical [70] investigations of supercritical combustion reported some similar concepts such as anchor vortex, recirculation zone, bulge shape in flame front and backward zone. All such similarities refer to specific volume of fluid which is related to density variations directly due to thermal expansion of flow field. Mardani et al. [70] evaluated decrease of density close to the injector which it causes to increase volume of fluid and diverts hydrogen jet to chamber wall. This physical phenomenon is the reason of subsequent resemblances that happen because of that deflection. But the differences between sub and supercritical can be clarified by considering the order of magnitude of properties value specially density. For example, change of density at subcritical condition is about 10 times, while change of density at supercritical condition is about 1000 times. In other word, the change of order of magnitude for density at sub and supercritical conditions are equal to 1 and 3, respectively. Table 4 H2 / O2 reduced chemical kinetic mechanism of Burke et al. [59]. Reaction 1. H O2 O OH 2. O H2 H OH O H2 H OH
3. H2 OH H2O H 4. OH OH O H2O 5. H2 M H H M 6. H2 Ar H H Ar 7. H2 He H H He 8. O O M O2 M 9. O O Ar O2 Ar 10. O O He O2 He 11. O H M OH M 12. O2H M H OH M 13. O2H O2H H OH O2H 14. H O2 (M) HO2 (M) 15. HO2 H H2 O2 16 HO2 H OH OH 17. HO2 O O2 OH 18. HO2 OH H2O O2 19. HO2 HO2 H2O2 O2
HO2 HO2 H2O2 O2 20. H2O2 (M) OH OH(M) 21. H2O2 H H2O OH 22. H2O2 H HO2 H2 23. H2O2 O OH HO2 24. H2O2 OH HO2 H2O H2O2 OH HO2 H2O
A H 2 / O 2 Chain Reactions 1.04 1014 3.82 1012
8.792 1014 2.16 108 3.34 104
n
E
0 0 0 1.51 2.42
1.53 104 7.95 103 1.917 104 3.43 103 1.93 103
H 2 / O 2 Dissociation/Recombination Reactions 4.58 1019 1.4 1.1 5.84 1018 18 1.1 5.84 10 15 6.16 10 0.5 1.89 1013 0 1.89 1013 0 18 4.71 10 1.0 27 3.32 6.06 10 1.01 10 26 2.44 Formation and Consumption of HO2 4.65 1012 0.44 2 .9 2.75 106 13 7.08 10 0 10 2.85 10 1 2.89 1013 0
Formation and Consumption of H2O2 4.20 1014
1.30 1011
2.0 1012
2.41 1013 4.82 1013
9.55 106
1.74 1012 7.59 1013
0 0 0 .9 0 0 2 .0 0 0
ACS Paragon Plus Environment
1.04 105 1.04 105 1.04 105 0 1.79 103 1.79 103 0 1.79 103 1.79 103 0 1.45 103 2.95 102 7.24 102 4.97 102
1.20 104 1.6293 103 4.88 104 3.97 103 7.95 103 3.97 103 3.18 102 7.27 103
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 10 of 31
10
Figure 6 Streamlines of combustion flow field, colored by static pressure distribution (Case no. 1 Table 3).
With respect to lack of sufficient experimental data for A10 test case, some geometrical parameters have been defined to compare numerical methods such as turbulent models, chemical mechanisms, and etc. Figure 8 shows OH* Abel transformed image in radial-axial coordinate system which is normalized by oxygen injector post diameter and some defined geometrical parameters. L is considered as a nondimensional axial location with constant spreading angle of maximum OH*. The measured magnitude of L is 5.7. In further locations, the spreading angle of OH* emission has been changed. is the spreading angle of maximum OH* emission intensity starts from oxygen injector lip to L and has the amount of 6.95°. Also the maximum spreading angle occurs at the point of M and two geometrical parameters of L x and L r are defined by this point. Their values are 7.2 and 2.1, respectively. Therefore, by using such geometrical parameters, there is a good approach to study and compare different models predictions. Additionally, some parameters have been defined in Figure 7 based on numerical solutions. RZ (recirculation zone) point is the location of center of recirculation. To compare the power of vortex in recirculation zone for different models, vorticity magnitude called VM (vorticity magnitude) is measured in mentioned location. SP (stagnation point) is the location of stagnation point where hydrogen jet strikes right into the chamber wall and velocity magnitude reaches to zero. Two points of SBZ (start of backward zone) and EBZ (end of backward zone) are used to locate entrance and exit positions of backward zone where axial velocity component has negative sign. Comparisons of predictions using different models in terms of such parameters help to understand the performances of different models on combustion flow field modeling.
Figure 7 Computed OH mass fraction, locations of zero axial (dashed-line), zero radial (Dashed-dotted line) velocities, and defined parameters (top), OH* Abel transformed image [68] and superimposed computed streamlines (bottom); (Case no. 1 Table 3).
ACS Paragon Plus Environment
Page 11 of 31
11
Figure 8 Experimental distribution of OH* emission intensity [68] and defined geometrical parameters.
To get a more details about understanding of the flow field, Figure 9a, b, c & d describe O2 mass fraction, static pressure, dynamic pressure, axial velocity, density, constant pressure specific heat, and temperature along the axis. As seen at about x / DLOX 30 , there are very small changes in axial velocity, C p , and temperature, but static and dynamic pressure smoothly increases and decreases, respectively, and density decrement and O2 mass fraction increment are done sharply. By moving to the downstream at x / DLOX 30 , due to drastic variation in temperature and in the location where the O2 starts to consume until O2 completely vanished, velocity and dynamic pressure increase monotonically while static pressure decreases. Indeed, all thermodynamic properties behave with respect to ideal gas equation of state. (a)
5000 60
80
0
40
x/DLOX
3000 2000 20
40
60
x/DLOX
80
Temperature [K]
CP Temperature
4000
(d)
6000
CP [J/kg.k]
2000 0
80
10000 15000 20000
1.03E+06
Static Pressure [Pa] 50 30 20 10
60
0
40
x/DLOX
Density [kg/m3]
40
200 100
Axial Velocity [m/s]
20
Axial velocity Density
20
1.02E+06
80
1000
60
0
40
x/DLOX
8000
20
4000
0.8 0.6 0.4 0
0.2
O2 Mass Fraction
Dynamic pressure Static Pressure
Dynamic Pressure [Pa]
1
(b) r/DLOX=0 r/DLOX=1 r/DLOX=2 r/DLOX=3 r/DLOX=4
(c)
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
Figure 9 Distributions of O 2 mass fraction in different axial directions (a), dynamic Pressure and static pressure (b), axial velocity and density (c), specific heat in constant pressure and temperature (d) in the axial direction; (Case no. 1 Table 3).
To focus on the discrete phase behaviour, distribution of spray volume fraction are shown at different axial sections in Figure 10a. As seen, the maximum values happen at sections of r / DLOX 0 & 1 while their magnitudes are less than the criterion of 103 . Thus, the dilute spray assumption to neglect the secondary breakup is approved. Characteristic diameters of SMD and cumulative undersize distribution at 10% ( Dvol .10% ) and 90% ( Dvol .90% ) are used to study variation of droplet diameters along the chamber which is shown in Figure 10b. Initial values of these three diameters can be obtained by particle distribution function which is represented by IWRCM in Rosin-Ramler form described in section 3. In view of the fact, SMD,
ACS Paragon Plus Environment
Energy & Fuels
12
Dvol .10% , and Dvol .90% are 82, 48, and 187 micrometers, respectively. As seen in this graph, three major events take place in terms of characteristic diameters. First, at x / DLOX 15 the characteristic diameters decrease slightly due to big size droplets relative to further axial locations which leads to absorb more energy to heat up and evaporation. Second, at 15 x / DLOX 55 the characteristic diameters experience severe gradients and decrease rapidly. Heating up to reach evaporation threshold for most of particles is done at x / DLOX 15 . Beyond this location and by passing through the flame, reduction of characteristic diameters grow excessively up to axial location of x / DLOX 35 . Afterwards, effects of flame reduce gradually and ratio of diameters reduction decreases to x / DLOX 55 . From the mathematical point of view, there is a turning point at 15 x / DLOX 55 where the sign of derivation of characteristic diameter curve line changes. The point might refer to position where the effect of flame on droplet evaporation starts to reduce at about x / DLOX 35 . Third, there is no footprint of significant exothermic reactions at 55 x / DLOX 80 , therefore intensity of diameter reduction decreases severely. Finally, all droplets disappear at x / DLOX 80 . Characteristic diameters variation of SMD, Dvol .10% and Dvol .90% were studied by Lempke et al. [24] who predicted the variation at x / DLOX 15 similarly. Also they reported that droplets are vanished at x / DLOX 80 . But major difference happens at 15 x / DLOX 55 where Lempke et al. didn’t report high rate of evaporation and the turning point. It might happen due to performing 19 steps kinetic mechanism and neglecting some exothermic reactions in contrast with 27 steps kinetic mechanism in this study. But these none-similar predictions should be verified by more detail experimental studies.
0
20
40
x/DLOX
60
80
150
200
Dvol.10% DSMD Dvol.90%
100
Characteristic Diameter [m]
250
(b)
50 0
4.0E-04
r/DLOX=0 r/DLOX=1 r/DLOX=2 r/DLOX=3
2.0E-04
Volume Fraction
(a)
0.0E+00
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 12 of 31
0
20
40
60
80
x/DLOX
Figure 10 Spray volume fraction distribution in different axial directions (a), and distribution of characteristic diameters in the axial direction of r / D LOX 1 (b); (Case no. 1 Table 3).
Following the combustion field through the chemical kinetic aspects could also be worthy. In this way, Figure 11 shows the oxidation pathways of H2 in a well-stirred reactor (WSR) constructed for the position with stoichiometric condition at temperature of 3386 K. Each arrow indicates the conversion of one species to another from tail to head. The related reaction and its rate are shown along the arrow with number and value in parenthesis. The number of reactions are according to Table 4. The reactions are filtered by rate of reaction limitation while the reactions with rates less than 1020 gmol / cm3 s are truncated. Furthermore, to better understanding the procedure of species conversion and pathways, arrows are visualized by different colors based on reaction rate values. This figure illustrates that H2 and every radical species have an exclusive direct pathway to produce H 2O . However it can be seen, two intermediate species of OH and H convert to H 2O with higher rates than the others. Li[71] and Burke[59] reported that the one of the most important reaction of H2 / O2 combustion is H O2 . This reaction has different behaviors for different pressures. When pressure is low and temperature is high, O2 turns to two radials of O. One of these two reacts with H and produces the determinant species of OH. Moreover, there are two primary
ACS Paragon Plus Environment
Page 13 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
13
routes specify the conversion of H2 to H 2O which are H HO2 H2O and OH H2O2 H2O . Indeed, conversions of radial species to each other construct some secondary routes as it well seen.
Figure 11 Oxidation pathways of H 2 with oxidizer of O 2 in a WSR constructed at stoichiometric condition for A10 test case while it is at T=3386 K and P=1 MPa. Reaction numbers refer to the presented reactions in Table 4 and the reaction rates are shown in parentheses in gmol / cm3 s .
According to Table 3, Burke chemical kinetic mechanism has been performed in reference case of no.1. Ignition reactions of Burke mechanism and their evaluated rate of reactions are described in Table 5. As seen, the reaction R1 has the maximum rate of reaction value in contrast to the others by a wide margin. Therefor the reaction R1 performs the primary role to initiate the chain of reactions. On the other side, high rate of reaction leads to the small chemical time scale which means reaction no.1 is important for the ignition time scale in such combustion field. The contour of rate of reaction R1 is plotted in Figure 12. It is apparent that in the chain of reactions, the ignition reaction R1 has contribution in most part of combustion field such as near the injector lip in location of anchor vortex. Table 5 Ignition reactions of Burke chemical kinetic mechanism and related rate of reaction margins (case no. 1 in Table 3). no. R1 R2 R5 R12
Reaction H O2 O OH O H2 H OH H2 M H H M O2 H M H OH M
Rate of reaction ( kg / m 3 .s ) 214 R .R . 60632 48 R .R . 3649 1.9E 26 R .R . 9.2E 28 1007 R .R . 8.0E 08
Figure 12 Rate of reaction R1 distribution as ignition reaction (case no. 1 in Table 3).
4.2 Effect of Compressibility ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 14 of 31
14
It is well known when pressure is low and temperature is too high, the assumption of incompressibility is approved, mainly because the intermolecular interactions can be ignored. Therefore, in the area of validity of ideal gas equation of state, two formats of incompressible and compressible may be used as done in some commercial CFD software packages such as Ansys-Fluent. If the density depends only on the operating pressure and not on the local relative pressure field, an incompressible ideal gas equation can be satisfied the physical phenomenon. Such a condition is relevant for many regular and industrial reacting flow fields. In this form, the flow field pressure is considered as a constant value in ideal gas EOS formulation to evaluate density field. In other words, the density depends only on the operating pressure and is not dependency of the local relative pressure field. In the subcritical condition and highly local Mach number, dependency of density on local pressure field might be significant. Therefore, it causes to involve local pressure values in ideal gas EOS formulation. However, in order to understand the effect of compressibility on flow field properties of current reacting test case, ideal gas equation is performed in two compressible and incompressible forms. All other conditions are the same in both modelling. These two exactly refer to cases no. 1&5 in Table 3. Contours of Mach number and compressibility factor are illustrated in Figure 13. With respect to Mach number distribution, Presence of flow near the threshold of M 0.3 necessitates to investigate the effect of compressibility, particularly in vicinity of the injector post. Figure 14a shows the compressibility factor variations at different axial sections while compressible form of ideal gas has been employed. Maximum deviation of compressibility factor from unity is equal to 4%, where happens at about r / DLOX 1 . Moreover, the maximum deviation at every axial section happens in the first eighth distance ( 0 x / DLOX 10 ). So it is expected that the effect of compressibility on flow properties are predicted in mentioned locations. The effect of such deviation can be seen in density and temperature as Figure 14b. For incompressible ideal gas, the density is deviated -9% from compressible form prediction and occurs in the first eighth distance at section of r / DLOX 1 . Afterwards, drastic differences in density profile disappear and both diagrams approach to similar values. Evaluated temperature profiles (Figure 14b) emphasize the resemblance of compressible and incompressible EOS results, even at the first eighth locations. Finally, it can be concluded that the assumption of incompressible ideal gas is suitable for the most part of the chamber. But in vicinity of the injector post, the story changes. Therefore, compressible form of ideal gas is used all over the modelling.
Figure 13 Distributions of computed Mach number (top) and compressibility factor (bottom); (Case no. 1 Table 3)
ACS Paragon Plus Environment
Page 15 of 31
15
40
60
80
20
x/DLOX
40
37
60
[k] Velocity [m/s] Temperature Axial
4000 3000 4000 3000 2000 2000 1000 1000
36
10
38
80
0
10 20
5
0
0
Temperature
Density
3300 3350
30
40
Incompressible Compressible
20
1 0.995
Density [Kg/m3]
1.02 1
r/DLOX=0 r/DLOX=1 r/DLOX=2 r/DLOX=3
0.98
Compressibility Factor
(b)
50
(a)
0.96
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
x/DLOX
Figure 14 Comparison of compressibility factor in different axial directions based on compressible EOS (a) (Case no. 1 Table 3), comparisons of density and temperature between compressible and incompressible EOS in the axial direction (b) (Case no. 1&5 Table 3).
4.3 Effect of turbulent model As it well seen in literature review according to the Table 1, the most popular turbulent models of are performed to study A10 combustion flow field till now in the other reports [27, 28, 30, 31, 33, 34, 36, 44]. In spite of the fact that some other models such as was also employed [24, 53], none of these studies concentrated on effects of different turbulent models on combustion and flame structure quantitatively and qualitatively. In order to study turbulence models effects, two-equation models of Standard, RNG, SST, and four-equation models of transition SST are used which refer to cases no. 1-4 in Table 3. Computed OH mass fraction based on different turbulence models are plotted alongside the Abel transformed OH* image in Figure 15 where the maximum values of prediction are shown with dashed-line. As it is seen, models belonged to family have more agreement to OH* image with respect to superimposed maximum values dashed-line. However, there are significant similarities between Standard and RNG. Such similarities happen for SST and Transition SST models, too. Comparison between predictions by using different turbulent models is done by implementing the parameters which are defined in Figure 8. In Table 6, three types of values are described. First is absolute measurements of defined parameters which have been obtained by performing different turbulence models. Second is absolute deviation percentage of mentioned parameters which are computed in contrast with experimental data. Last is the overall error percentage of each turbulence model which has been evaluated by mean average. Computed values describe that all turbulent models have specific error. Based on, it is obviously concluded that family, consist of Standard and RNG, have better performance, especially Standard model. In addition, with respect to defined parameters in Figure 7, locations of center of recirculation zone (RZ), stagnation point (SP), start of backward zone (SBZ), end of backward zone (EBZ), and also vorticity magnitude (VM) are shown for different turbulent models in Table 7. Due to lack of sufficient experimental data about such parameters, related deviations are computed from the best case of no.1 in Table 3. It is seen, deviation of RNG as a member of family is approximately half of family related values.
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 16 of 31
16
Figure 15 Computed OH mass fraction distributions based on different turbulence models (top) (Case no. 1-4 Table 3) in comparison with experimental Abel transformed OH*image [68] (bottom). Table 6 Comparison between defined parameters obtained by performing different turbulent models (Case no. 1-4 Table 3) and experimental transformed OH* image [68] in terms of values and deviations. Absolute parameter values Turbulence model Experiment Standard RNG SST Trans. SST
6.95 7.61 7.95 13.69 12.79
L
Lx
5.7 5.3 5.5 4.7 4.2
7.2 7.3 7.7 6.0 6.0
Lr
2.1 2.2 2.4 2.5 2.4
%
Absolute deviation
%L
9.5 14.3 96.9 84.4
%L x
7.0 3.5 17.5 26.3
Average deviation
%error
%Lr
1.3 6.9 16.6 16.6
4.7 14.2 19.0 14.2
5.6 9.72 37.5 35.4
Table 7 Comparison between turbulent models (Case no. 1-4 Table 3) in terms of measures, and related deviations from Case no.1 parameter values Turbulence model Standard RNG SST Trans. SST
VM ( S 1 ) 23691.9 26974.3 25967.8 27685.5
Point RZ (5.56,3.93) (6.11,4.02) (4.13,4.07) (4.32,4.04)
Point SP (7.70,5) (8.03,5) (6.04,5) (6.24,5)
Point SBZ (5.98,1.31) (5.66,1.24) (5.43,1.85) (5.44,1.83)
Absolute deviation Point EBZ (12.31,1.92) (13.85,1.79) (10.64,2.36) (10.69,2.28)
%VM 13.8 9.6 16.8
%RZ 7.3 14.8 13.2
%SP 3.0 14.5 12.9
%SBZ 5.3 6.3 6.2
%EBZ 12.0 12.5 12.3
Average deviation %dev. 8.3 11.5 12.3
To complete the discussion, temperature and OH mass fraction distributions are illustrated in Figure 16 based on different turbulent models. The figure depicts differences in temperature of gaseous oxygen zone which family evaluates higher values. Maximum temperature zone for family has an entirely solid border. Nevertheless, it is represented by two detached zone for family. Also the axial location of maximum temperature zone is minimum for Standard. The length and thickness of maximum temperature zone directly depends on chemical properties. Due to turbulence-chemistry interaction, concentration of species are estimated dissimilarly for different turbulence models (see for instance Figure 16b) and it goes to different temperature field. Also another reason is the effect of turbulent models on prediction of the flow field properties. Streamlines visualize flow field behavior for different turbulent models which has been explained in Table 6 and Table 7, quantitatively and qualitatively for each models.
ACS Paragon Plus Environment
Page 17 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
17
Figure 16 Effect of performing different turbulent models on temperature (a) and OH mass fraction (b) of combustion field; (Case no. 1-4 Table 3)
Figure 17 describes the axial distributions of Temperature, O2 mass fraction, density, OH mass fraction, axial velocity, and specific heat in constant pressure for different employed turbulent models. Temperature distribution is shown in range of 0 x / DLOX 80 to illustrate variations in overall axial distance of the chamber. As it can be seen, major difference between predictions using various turbulent models happens close to the injector and in fact, source of turbulent energy. Therefore, to investigate temperature variation in that location, other properties are plotted in 0 x / DLOX 25 . In this figure, temperature distribution approves recent conclusion which Standard has the most precise prediction to experimental data. From turbulence point of view, it seems that the energy of turbulent eddies and so the combustion flow behaviour is controlled by transition of jet flow stream which relates to the type of employed turbulent model. Consequently, different turbulent models have two major effects on the flow field. First direct effect is due to turbulent-chemical interaction which it causes to different chemical field and species concentrations. Second is the effect of turbulence model on flow field properties due to different source terms in conservation equations which refers to different turbulent viscosity ratio and mixing level. With respect to the Figure 17, temperature profiles could be described at three important axial ranges. The first difference between and SST turbulence model families happens at 3.5 x / DLOX 6.5 , where the minimum temperature value is located in earlier axial position for SST models. It is (may be) due to type of turbulent model and its source terms effects on flow field properties such as axial velocity and density at such range. In little further locations, a drastic positive gradient causes temperature gains maximum value, but in different axial positions. In the location of maximum temperature, because of turbulent-chemistry interactions, SST family predicts double peaks for OH mass fraction and therefore, its consequences can be seen in temperature variations. Nevertheless, distribution of C P doesn’t emphasize such treatment and it means that prediction of double peaks is completely relates to turbulent-chemistry interactions and chemical field properties. Direct effect of C P can be seen at 15 x / DLOX 30 , where SST family with lower C P conducts higher temperature in contrast with family. As the distance from the injector and source of turbulent energy increases, effect of different turbulent models decreases, and evaluation of flow
ACS Paragon Plus Environment
Energy & Fuels
18
15
20
25
5
10
15
20
25
5
10
15
20
x/DLOX
25
0
5
10
15
20
25
5
10
15
20
25
x/DLOX
CP [J/kg.K]
150
8000 4000
0 0
0
0
12000
x/DLOX
0.6
0.7 10
0.5
Density [kg/m3] 5
100
Axial Velocity [m/s]
0
0.4
0.3 0.2
80
0.1
x/DLOX
60
0
10
40
0.05
OH Mass Fraction
20
50
3150
0.1
3300
2000
5
0
O2 Mass Fraction
3000
EXP [55] Standard RNG SST Trans. SST
1000
Temperature [K]
field properties is roughly independent from turbulent model. So, the temperature distributions approach to similar values at the range of x / DLOX 30 for different turbulence models approximately.
x/DLOX
0
x/DLOX
Figure 17 Axial profile distributions of temperature and related experimental data [55], O 2 mass fraction, density, OH mass fraction, axial velocity, and specific heat in constant pressure for different turbulent models (Case no. 1-4 Table 3).
By considering the maximum intensity position of predicted OH mass fraction as the flame front, Figure 18 illustrates the comparison between the flame shapes of four cases employed different turbulent models and maximum intensity of Abel transformed OH* image for A10 test case. In general, comparison in axial range of 0 x / DLOX 15.6 , the standard seems to have the most agreement to experimental data. After that, the RNG has proper result. It is obviously seen that the experimental data is limited to axial range of 0 x / DLOX 15.6 . Therefore, due to such limitation, there is no way to compare length of flames with experimental data. However, the flame lengths obtained by Standard is 36D. By applying the models of RNG, SST, and transition SST, the flame length to be 8.8%, 14.4%, and 13%, respectively, more than Standard. 5
r/DLOX
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 18 of 31
EXP [68] Standard RNG SST Trans. SST
0
-5
10
20
x/DLOX
30
40
50
Figure 18 Comparison between position of maximum intensity of Abel transformed of OH* image [68] and computed OH mass fraction for different turbulent models (Case no. 1-4 Table 3).
ACS Paragon Plus Environment
Page 19 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
19
4.4 Effect of chemical kinetic mechanism In this section, the effect of different chemical kinetic mechanisms is investigated. Four different multi-step chemical mechanisms with various range of operating conditions and details have been employed. Also, the one-step global reaction mechanism is used to specify importance of using detailed mechanisms to predict combustion field properties accurately in the Case no. 7 of Table 3. In Figure 19, the computed temperature fields are tabulated for one-step and Burke [59] multi-step reaction mechanisms. It is well seen that the structure of combustion flow field for the two types of mechanisms are approximately the same from velocity field aspects such as spreading angle, point M location, length of flame, etc., but the temperature fields show great differences quantitatively. The one-step mechanism temperature has a maximum of 4400 K. However, such value for multistep mechanism is about 3400 K.
Figure 19 Distribution of temperature in combustion flow field for Burke [59] multi-step reaction mechanism (top) and one-step reaction mechanism (bottom); (Case no. 1&7 Table 3).
The Konnov [61] mechanism was presented in 2008 to specify different channels of reactions between H and HO2 radicals in various range of pressure. The H2 / O2 Konnov mechanism has 31 reactions with 10 chemical species to cover a range of 950 to 2700 K for temperature and 0.3 up to 87 atm for pressure. In 2014, the Alekseev [62] represented an update on Konnov mechanism to improve its performance for jet stirred, flow reactor, oxidation, decomposition, and shock tube studies in a wide range of equivalence ratios. In the Alekseev mechanism, 30 step reactions with 10 different species are used for predictions of the experiments in temperature range of 900 to 2700 K and pressure range of 0.3 up to 87 atm. The detailed kinetic mechanism of Óconaire [60] was developed in 2004 to simulate the H2 / O2 combustion in wide range of temperature and pressure of 298 to 2700 K and 0.05 to 87 atm, respectively. Óconaire mechanism consists of 21 step reactions and 10 chemical species. In 2012, Burke [59] introduced an update version of Li H2 / O2 combustion mechanism for a wider operation condition of 0.3 to 87 atm for pressure range, 298 up to 3000 K for temperature, and 0.25 to 5 for equivalence ratio to handle chemical problems such as shock tubes, flow reactors, and laminar premixed flames which showed a good agreement with experimental data. The Burke reaction model has 13 different species used in 27 step reactions which has been employed as the major mechanism in this study (Case no.1 Table 3). Figure 20 illustrates temperature, OH mass fraction, and superimposed streamlines which are computed by employing four mentioned chemical mechanisms. Comparisons of results show significant similarities, especially for three models of Burke, Konnov, and Óconaire. For Alekseev mechanism, the zone where liquid oxygen evaporates, is longer than the others. Also, the 3200 K isothermal cut-off line shows some differences for Alekseev mechanism in Figure 20a. In order to make a precise comparison, deviations of mentioned geometrical characteristics (in Figure 7) are tabulated in Table 8 and Table 9. It can be obviously realized that the Burke mechanism has slightly better performance than the Konnov and Óconaire for predictions of the measurements. However, the results based on Alekseev mechanism conduct the maximum error.
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 20 of 31
20
Figure 20 Effects of employing different multistep reaction mechanisms on temperature (a) and OH mass fraction (b) fields; (Case no. 1, 8-10 in Table 3). Table 8 Comparison of defined parameters obtained from numerical modelling using different chemical kinetic mechanisms (Case no. 1, 8-10 Table 3) and experimental transformed OH* image [68] in terms of values and deviations. Absolute parameter values Mechanism
L
Lx
Experiment Burke Alekseev Konnov Óconaire
6.95 7.61 7.96 7.68 7.71
5.7 5.3 5.0 5.2 5.3
7.2 7.3 8.0 7.3 7.3
Absolute deviation
Average deviation
Lr
%
% L
% Lx
% Lr
%error
2.1 2.2 2.0 2.2 2.2
9.5 14.5 10.5 10.9
7.0 12.2 8.8 7.0
1.3 11.1 1.3 1.3
4.7 4.7 4.7 4.7
5.6 10.6 6.3 6.0
Table 9 Comparison of defined parameters obtained from numerical modelling using different chemical kinetic mechanisms (Case no. 1, 8-10 Table 3) in terms of values, and related deviations from Case no.1. parameter values Turbulence model Burke Alekseev Konnov Óconaire
VM ( S 1 ) 23691.9 21984.0 23900.6 24120.5
Average deviation
Absolute deviation
Point RZ
Point SP
Point SBZ
Point EBZ
%VM
%RZ
%SP
%SBZ
%EBZ
%dev.
(5.56,3.93) (6.01,3.89) (5.65,3.92) (5.64,3.93)
(7.70,5) (8.32,5) (7.77,5) (7.77,5)
(5.98,1.31) (6.54,1.27) (6.03,1.28) (5.96,1.27)
(12.31,1.92) (12.53,1.82) (12.53,1.85) (12.46,1.89)
7.2 0.9 1.8
5.1 1.0 0
5.8 0 0
1.0 0 0
1.6 1.6 1.1
2.98 0.7 0.58
Variations of predicted radial temperatures by different multi-step and the one-step global mechanisms are shown in Figure 21a. First of all, maximum temperature for global reaction has been estimated about 1000 K more in comparison with multi-step mechanisms. Moreover, three mechanisms of Burke, Konnov, and Óconaire have predicted the maximum value of temperature similarly, but evaluated maximum
ACS Paragon Plus Environment
Page 21 of 31
21
temperature by Alekseev is 100 K less. The mass fractions of OH and H2 species are plotted in Figure 21b. OH mass fractions obtained by Burke, Konnov, and Óconaire have a maximum 2% deviation from each other. But for Alekseev, the deviation from the others is around to 13%. Also, the H2 mass fraction profiles emphasize that maximum deviation of multi-step mechanisms predictions from each other is around 2% which it belongs to Alekseev. Positions of maximum OH mass fraction are plotted for different employed mechanisms in Figure 22 and compared with the measurements. As it is well seen, Burke mechanism has the most agreement with available experimental data. In addition, the maximum deviation from Burke’s length of flame is about 2% which is for Alekseev. Finally, it can be concluded that Burke, Konnov, and Óconaire deliver a significant similar performance to predict A10 combustion flow field properties.
(a)
1
2
4
5
0
1
r/DLOX
2
3
0.4 0.2
0.8
H2 Mass Fraction
0.6
0.95 0.94 0.93 0.92
0.7
0.91
OH
0.108 0.11 0.112
0.8
0.1
OH Mass Fraction
0.02 0.04 0.06 0.08
0.8
3
H2
0
0
0.7
0
2000
3150 3225
3000
EXP [55] Global Burke Alekseev Konnov Oconaire
1000
Temperature [K]
4000
1
(b)
4
5
r/DLOX
Figure 21 Radial profiles of temperature and related experimental data [55] (a), OH and H 2 mass fractions (b), for one/multistep reaction mechanisms; (Case no. 1, 8-10 in Table 3).
5 EXP [68] Burke Alekseev Konnov Óconaire
r/DLOX
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
0
-5
10
20
x/DLOX
30
40
Figure 22 Comparison between position of maximum intensity of Abel transformed of OH* image [68] and computed OH mass fraction for different chemical kinetic mechanisms (Case no. 1, 8-10 Table 3).
4.5 Effect of mixture ratio One of the most important parameter which has been investigated in experimental studies is the effect of mixture ratio changes at a constant pressure on the flame structure. In this way, A10 and C10 are the two test cases which have been represented to comprehend the effect of mixture ratio change on flame structure at subcritical condition. Operation conditions described in Table 2Error! Reference source not found. show
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 22 of 31
22
the fuel-rich condition for A10 and C10 with the equivalence ratios of 3.8 and 2.5, respectively. In this section, the flame structures for different mixture ratios of A10 and C10 are numerically investigated. It should be mentioned that Jay et al. studied C10 by an Eulerian framework [30]. However, there is no numerical modeling of C10 using a Lagrangian point of view. With respect to best methods qualified in previous sections, C10 is modeled as it is described in Table 3. Similar to case A10, there are few experimental data to compare and validate numerical results for C10. Figure 23 shows Abel transformed of OH* image alongside the OH mass fraction computed from numerical modeling for C10. As it can be seen, there is a good agreement with the available experimental data.
Figure 23 Comparison of OH mass fraction distribution obtained from numerical modelling of cryogenic combustion of Mascotte C10 (test case no. 11) (top) and Abel transformed OH* emission image [68] (bottom)
However to get a more precise description, a detailed comparison is shown in Table 10 with respect to defined parameters in Figure 8. First of all, comparisons of , L , L x and L r for experimental A10 and C10 (Case no. 1&11 in Table 3) studies emphasize some physical differences. Spreading angle of C10 test case is about 2˚ more than the A10’s due to lower H2 to O2 momentum flux. H2 stream with lower injected velocity is deflected easier under expansion of evaporated O2 for C10. Also, due to such reason and measurements of L x and L r , Bulge of C10 is located upward and rearward of A10’s. Computed spreading angle has a deviation of 6.6% in contrast to Abel transformed OH* image from Table 10. Based on L , L x , and L r values and related deviations in that table, it can be seen that the modeling of C10 has good agreement with experimental data. Figure 24 shows a comparison between predicted radial distribution of temperature with available experimental data. Both mean values and standard deviation of temperature gathered from measurements are shown in this figure. Comparison of computed temperature from numerical modelling with experimental data emphasizes proper agreement. Table 10 Comparison of defined parameters obtained from numerical modelling of C10 combustion (Case no.11 Table 3) with related experimental transformed OH* image [68] in terms of values and deviations Absolute parameter values
C10 Experiment C10 Prediction
Absolute deviation
Average deviation
L
Lx
Lr
%
%L
%Lx
%L r
%error
9.0 10.4
6.11 5.96
6.73 6.88
2.21 2.35
15.5
2.4
2.2
6.3
6.6
ACS Paragon Plus Environment
Page 23 of 31
500 0
1
2
3
r/DLOX
4
400
450
2000 1000
Temperature [K]
3000
EXP [55] EXP Dev. [55] Present study
Standard Deviation [K]
23
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
5
Figure 24 Comparison of computed radial profile of temperature with mean and standard deviation of CARS experimental results [55].
Figure 25 shows the distributions of temperature and OH mass fraction and superimposed streamlines in A10 and C10 flow fields. Comparisons of results show that the flame of C10 has longer length due to effect of less fuel-rich in contrast to A10. In other words, existence of unburnt fuel causes to temperature reduction due to cooling effect. So, the temperature decreases drastically right after consumption of O2 in A10 test case. Furthermore, based on comparison done in Table 11, vortex strength of C10 is 29.7% weaker than the A10’s which it causes to increase length of flame. At x / DLOX 12 , the maximum temperature zone is thicker for C10 because of less H2 to O2 momentum flux ratio, less fuel-rich condition, and less cooling effect as it is discussed before. Moreover, the length and thickness of OH contour have similar behavior to temperature one in combustion flow field. Due to weaker strength of vortex in C10, it is driven more to the corner of chamber by H2 stream and also expansion effect of evaporated O2 . To get more precise information about differences between A10 and C10 streams, Table 11 describes the parameters defined in Figure 7. By considering the coordinates of RZ and SP points, it can be realized that the center of recirculation zone and stagnation point for C10 locate 5.1% and 4.9%, relatively closer to the injector plate than the A10’s. Also SBZ and EBZ coordinates for C10 emphasize that the backward zone happens 3.9% upper and 5.2% closer to the injector plate.
Figure 25 Effects of different mixture ratio of A10 and C10 test cases on temperature (a) and OH mass fraction (b) fields; (Case no. 1&11 in Table 3)
ACS Paragon Plus Environment
Energy & Fuels
24 Table 11 Comparison of defined parameters obtained from numerical modelling of A10 and C10 combustion (Case no. 1&11 Table 3) in terms of values and deviations parameter values Test Case A10 C10
VM ( s 1 ) 23691.9 16640.2
Average deviation
Absolute deviation
Point RZ
Point SP
Point SBZ
Point EBZ
%VM
%RZ
%SP
%SBZ
%EBZ
%dev.
(5.56,3.93) (6.01,3.89)
(7.70,5) (7.16,5)
(5.98,1.31) (5.71,1.41)
(12.31,1.92) (11.63,2.07)
29.7
5.1
4.9
3.9
5.2
9.76
40
60
80
30 20
3
80
10
Density [kg/m ]
60
0
40
x/DLOX
20
40
60
80
40
60
80
x/DLOX
0.1
6000
20
4000
CP [J/kg.K]
Temperature [K]
A10 C10
2000
0.05
OH Mass Fraction
0
1000 2000 3000 4000
In Figure 26 axial profiles of temperature, density, OH mass fraction, and specific heat in constant pressure are shown at the axis direction. Maximum temperature of A10 and C10 are 3333 K and 3361 K, respectively. Therefore, by O/F mixture ratio increase under fuel-rich condition, not only does the location of maximum temperature happen at early location, but also its value increases and vice versa. By mixture ratio increase, C P profiles moves backwards in axial direction. In the same attitude, it happens for density and OH mass fraction axial profiles. To understand such behavior, it is necessary to focus on flame front and its effective characteristics, especially vortex strength as it is discussed above. However, C P starts to rise at x / DLOX 22.7 and x / DLOX 30.6 for A10 and C10, respectively. Thus, in these locations, a rapid rise in temperature happens. From a chemical point of view, the location of maximum temperature can be anticipated from OH mass fraction profiles. The peak values of OH mass fraction for A10 and C10 are 0.114 and 0.117, respectively, which cause to a slight increment in C10 maximum temperature. Moreover, temperature variations are in relation with density profiles conversely, in other words, the higher the temperature, the lower the density.
0
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 24 of 31
20
x/DLOX
20
x/DLOX
Figure 26 Axial profile distributions of temperature ,density, OH mass fraction, and specific heat in constant pressure for A10 and C10 combustion test cases (Case no. 1&11 Table 3).
Figure 27 illustrates the comparisons of mass fraction of H, O2 , H2O2 , and HO2 for A10 and C10 test cases. The effect of change of mixture ratio have been investigated by well-stirred reactor modelling at stoichiometric condition for C10 such as the A10 test case. The results consist of pathways and reaction rates resemble for both test cases closely. Therefore, it can be said that the species concentration differences should be explored in flow dynamics characteristics of A10 and C10. With respect to Figure 11, H species is produced by H2 through the pathway of H2 H by R2 ( H2 O H OH ) and R3 ( H2 OH H2O H ) reactions with rates of 2.25 103 and 4.04 103 gmol / m3 s , respectively. It is expected that the H radical is produced in vicinity of H2 and evaporated O2 interaction layer. The axial spreading distance of H radical can be explained by getting focus on flow dynamics characteristics such as vortex strength.
ACS Paragon Plus Environment
Page 25 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
25
H2O2 is produced by the reaction of R20 ( H2O2 ( M) OH OH( M) ) through the pathway of
OH H2O2 with rate of 6.50 103 gmol / m3 s . Due to less H2 mass flow rate in C10 test case, fraction of unreacted evaporated O2 enters in recirculation zone. Based on oxidation pathway diagram in
Figure 11, the only source of H2O2 production is OH which is produced by three pathways of H2 OH , H OH , and HO2 OH . Between A10 and C10, H and H2 distributions have no major differences in the recirculation zone, but HO2 has sufficient mass fraction in that zone for C10. It can be concluded that the rout of HO2 OH H2O2 emphasizes the presence of H2O2 in recirculation zone. Moreover, in Figure 27, HO2 concentration is plotted for A10 and C10 test cases. There are two major differences with various sources in HO2 distribution. The first is presence of more HO2 in recirculation zone and the second is longer axial spreading distance, both for C10. According to the Figure 11, three different pathways lead to produce HO2 which are H HO2 , OH HO2 , and H2O2 HO2 . The first pathway of H HO2 is done through R14 ( H O2 ( M) HO2 ( M) ) and R22 ( H2O2 H HO2 H2 ) with rates of 2.95 102 and 7.71 104 gmol / m3 s , respectively. As it can be seen, presence of H radical and evaporated O2 causes to activate R14 with significant reaction rate which leads to produce HO2 in recirculation zone. The area related to simultaneous presence of R14 reactants locates the axial spreading of produced HO2 . In addition, HO2 is produced by the reaction R24 ( H2O2 OH HO2 H2O ) through the second pathway of OH HO2 with a rate of 2.97 103 gmol / m3 s . As it can be seen in Figure 27, existence of H2O2 in recirculation zone activates the R24 reaction to produce HO2 . The final sources of HO2 production are reactions of R22 ( H2O2 H HO2 H2 ), R23 ( H2O2 O OH HO2 ), and R24 ( H2O2 OH HO2 H2O ) through
the third pathway of H2O2 HO2 with rates of 7.71 104 , 1.47 103 , and 2.97 103 gmol / m3 s , respectively. The area related to H2O2 presence explains HO2 production by R22 and R23 reactions.
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 26 of 31
26
Figure 27 Comparison of mass fractions of H, O 2 , H 2 O2 , and HO2 species between A10 (top) and C10 (bottom) test cases (Case no. 1&11 Table 3).
5. Conclusion In this study, the combustion of gaseous hydrogen and liquid oxygen has been investigated numerically at subcritical condition by means of eleven different test cases. For this purpose, Mascotte test facility has been selected to study. Eulerian-Lagrangian framework has been used to specify equations of two-phase flow in a 2D axisymmetric domain. Continues phase has been modelled by employing RANS approach under steady state condition. To consider effect of compressibility, compressible ideal gas equation has been used to compute density field. Finite rate chemistry approach has been utilized for elementary reaction rates calculations. The net rate of production of each species has been obtained by turbulence-chemistry interaction model of Eddy Dissipation Concept (EDC). To evaluate turbulent closure, Standard are employed. The reduced kinetic mechanism of Burke et al. was selected which consists of 13 species and 27 step reactions. Also Discrete phase model (DPM) has been performed to evaluate dispersed phase properties. Droplets dispersion in turbulent flow field have been predicted by performing eddy life time model. It is cleared that two major parameters of vortex strength in recirculation zone and turbulence-chemistry interaction effects explain the flame configuration. In addition, anchor vortex stabilizes the flame close to the injector lip. It is found that two primary routes of H HO2 H2O and OH H2O2 H2O specify the conversion of H2 into H 2O at stoichiometric condition. Also the ignition reaction of H O2 O OH performs the main role to initiate the chain of reactions. The results of study could further categorize as below. Effect of compressibility. In order to understand the effect of compressibility, incompressible ideal gas EOS was substituted by compressible form of ideal gas EOS. The maximum deviation of compressibility factor from unity is about 4% happens near the gaseous hydrogen injection surface. Due to significant deviation in
ACS Paragon Plus Environment
Page 27 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
27
prediction of density, about 9%, the assumption of incompressible was rejected and compressible ideal gas equation was performed all over the modelling. Effects of turbulent models. Two-equation models of Standard, RNG, SST, and fourequation model of Transition SST were examined, and the results showed that the Standard turbulence model has the most agreement in comparison to the available experimental data (i.e., OH* image, related defined parameters, and temperature distribution). However, it can be said that the results of family, consists of Standard and RNG models, are close to each other. Such similarity happened for family, consists of SST and Transition models. Major differences for family results are shown in flame length where due to lack of sufficient experimental data, there is no chance to judge about. Effect of chemical kinetic mechanism. Four reduced H2 / O2 chemical kinetic mechanisms of Burke et al. ,O’conaire et al., Konnov, and Alekseev et al. have been used which consist of several elementary reactions and different pre-exponential factors depend on pressure and temperature. It seems that the performance of three mechanisms of Burke, Konnov, and O’conaire are the same in aspects of flame structure and flow field properties. Effect of mixture ratio. Two test cases of A10 and C10 are used to investigate numerically the effect of different mixture ratios at a constant pressure. There is good agreement between predictions and available experimental data (i.e., OH* image, related defined parameters, and temperature distribution). Spreading angle of C10 test case is more than the A10’s due to low H2 to O2 momentum flux ratio. Also based on measurements of L x and L r , Bulge of C10 is located upward and rearward of A10’s. The flame of C10 has longer length due to effect of less unburnt fuel in contrast to A10. In addition, the vortex of recirculation zone is driven more to the corner of chamber by H2 stream and expansion effect of evaporated O2 for C10. Furthermore, it is found that the oxidation pathways for C10 has no differences with A10 neither routes, nor rates.
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 28 of 31
28
Acknowledgements The authors would like to thank Mr. Ehsan Barani for his help regarding the current research. References 1. 2. 3. 4. 5.
6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16.
17. 18.
19. 20.
21.
Newman, J.A. and T.A. Brzustowski, Behavior of a Liquid Jet near the Thermodynamic Critical Region. AIAA Journal, 1971. 9(8): p. 1595-1602. Candel, S., et al., Structure and Dynamics of Cryogenic Flames at Supercritical Pressure. Combustion Science and Technology, 2006. 178(1-3): p. 161-192. Mayer, W. and H. Tamura, Propellant injection in a liquid oxygen/gaseous hydrogen rocket engine. Journal of Propulsion and Power, 1996. 12(6): p. 1137-1147. Oschwald, M., et al., Injection of Fluids into Supercritical Environment. Combustion Science and Technology, 2006. 178(1-3): p. 49-100. H. Mayer, W.O., et al., Atomization and Breakup of Cryogenic Propellants Under HighPressure Subcritical and Supercritical Conditions. Journal of Propulsion and Power, 1998. 14(5): p. 835-842. Sato, J.i., et al., Effects of natural convection on high-pressure droplet combustion. Combustion and Flame, 1990. 82(2): p. 142-150. Sankaran, V. and S. Menon †, LES of spray combustion in swirling flows. Journal of Turbulence, 2002. 3: p. N11. Ashgriz, N., Handbook of Atomization and Sprays: Theory and Applications. 2011: Springer US. Faghri, A. and Y. Zhang, Transport Phenomena in Multiphase Systems. 2006: Elsevier Academic Press. Lefebvre, A., Atomization and Sprays. 1988: Taylor & Francis. Maxey, M.R. and B.K. Patel, Localized force representations for particles sedimenting in Stokes flow. International Journal of Multiphase Flow, 2001. 27(9): p. 1603-1626. Crowe, C.T., et al., Multiphase Flows with Droplets and Particles. 1997: Taylor & Francis. Ferry, J. and S. Balachandar, A fast Eulerian method for disperse two-phase flow. International Journal of Multiphase Flow, 2001. 27(7): p. 1199-1226. Whitaker, S., The Method of Volume Averaging. Springer Netherlands, 1999: p. XVI, 210. Zhang, D.Z. and A. Prosperetti, Averaged equations for inviscid disperse two-phase flow. Journal of Fluid Mechanics, 1994. 267: p. 185-219. Grisch, F., et al., CARS Measurements at High Pressure in Cryogenic LOX/GH2 Jet Flames, in Liquid Rocket Thrust Chambers. 2004, American Institute of Aeronautics and Astronautics. p. 369-404. Vingert, L., et al., Atomization in Coaxial-Jet Injectors, in Liquid Rocket Thrust Chambers. 2004, American Institute of Aeronautics and Astronautics. p. 105-140. Froehlke, K., et al., First hot fire test campaign at the French/German research facility P8, in 33rd Joint Propulsion Conference and Exhibit. 1997, American Institute of Aeronautics and Astronautics. Haidn, O.J. and M. Habiballah, Research on high pressure cryogenic combustion. Aerospace Science and Technology, 2003. 7(6): p. 473-491. Haidn, O., et al., CFD Analysis of a Model Combustor Ignition and Comparison with Experimental Results, in 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. 2007, American Institute of Aeronautics and Astronautics. Santoro, R., et al., Testing at university facilities, in 39th Aerospace Sciences Meeting and Exhibit. 2001, American Institute of Aeronautics and Astronautics.
ACS Paragon Plus Environment
Page 29 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
29
22.
23.
24.
25. 26. 27.
28.
29.
30. 31.
32. 33. 34. 35. 36. 37. 38. 39. 40.
41. 42.
43.
Boyer, J.E., et al., Propulsion Education at the Pennsylvania State University, in 53rd AIAA/SAE/ASEE Joint Propulsion Conference. 2017, American Institute of Aeronautics and Astronautics. Santoro, R., University Propulsion Programs at Penn State, in 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit. 2004, American Institute of Aeronautics and Astronautics. Lempke, M., et al., Steady and Unsteady RANS Simulations of Cryogenic Rocket Combustors, in 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition. 2011, American Institute of Aeronautics and Astronautics. Conley, A., A. Vaidyanathan, and C. Segal, Heat Flux Measurements for a GO2/GH2 SingleElement, Shear Injector. Journal of Spacecraft and Rockets, 2007. 44(3): p. 633-639. Chen, Y.-S., et al., Multiphysics simulations of rocket engine combustion. Computers & Fluids, 2011. 45(1): p. 29-36. Farmer, R., G. Cheng, and Y.-S. Chen, CFD Simulation of Liquid Rocket Engine Injectors: Part 2. Simulations of the RCM-2 Experiment. Proceedings of the Second International Workshop on Rocket Combustion Modeling, 2001. Cheng, G. and R. Farmer, CFD spray combustion model for liquid rocket engine injector analyses, in 40th AIAA Aerospace Sciences Meeting & Exhibit. 2002, American Institute of Aeronautics and Astronautics. Liang, P. and R. Ungewitter, Multi-phase simulations of coaxial injector combustion, in 30th Aerospace Sciences Meeting and Exhibit. 1992, American Institute of Aeronautics and Astronautics. Jay, S., F. Lacas, and S. Candel, Combined surface density concepts for dense spray combustion. Combustion and Flame, 2006. 144(3): p. 558-577. Bodèle, E., et al., Modeling of MASCOTTE 10 Bar Case with THESEE with andWithout a Secondary Atomization Model. Proceedings of the Second International Workshop on Rocket Combustion Modeling, 2001. Haiden, O.J. 2ND International Workshop on Rocket Combustion Modeling: Atomization Combustion and Heat Transfer held in Lampoldshausen. 2001. Pourouchottamane, M., et al., Numerical Analysis of the 10 Bar MASCOTTE Flow Field. Proceedings of the Second International Workshop on Rocket Combustion Modeling, 2001. Izard, J.-F. and A. Mura, Lagrangian modeling of turbulent spray combustion: application to rocket engines cryogenic conditions. Progress in Propulsion Physics, 2011. 2: p. 207-224. Williams, F.A., Spray Combustion and Atomization. 1958. 1(6): p. 541-545. Gomet, L., V. Robin, and A. Mura, Lagrangian modelling of turbulent spray combustion under liquid rocket engine conditions. Acta Astronautica, 2014. 94(1): p. 184-197. Park, T.S., LES and RANS simulations of cryogenic liquid nitrogen jets. The Journal of Supercritical Fluids, 2012. 72: p. 232-247. Juniper, M., N. Darabiha, and S. Candel, The extinction limits of a hydrogen counterflow diffusion flame above liquid oxygen. Combustion and Flame, 2003. 135(1): p. 87-96. Masquelet, M., et al., Simulation of unsteady combustion in a LOX-GH2 fueled rocket engine. Aerospace Science and Technology, 2009. 13(8): p. 466-474. Oefelein, J., Large Eddy Simulation of a Shear-Coaxial LOX-H2 Jet at Supercritical Pressure, in 38th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit. 2002, American Institute of Aeronautics and Astronautics. Bellan, J., Theory, Modeling and Analysis of Turbulent Supercritical Mixing. Combustion Science and Technology, 2006. 178(1-3): p. 253-281. Okong'o, N., K. Harstad, and J. Bellan, Direct numerical simulations of O2/H2 temporal mixing layers under supercritical conditions, in 40th AIAA Aerospace Sciences Meeting & Exhibit. 2002, American Institute of Aeronautics and Astronautics. Foster, J. and R. Miller, A Priori Analysis of Subgrid Mass Flux Vectors from Massively Parallel Direct Numerical Simulations of High Pressure H2/O2 Reacting Shear Layers. American Physical Society, 64th Annual Meeting of the APS Division of Fluid Dynamics, 2011: p. D23.009.
ACS Paragon Plus Environment
Energy & Fuels 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Page 30 of 31
30
44. 45.
46. 47. 48. 49. 50. 51. 52. 53. 54. 55. 56. 57.
58.
59. 60. 61. 62.
63. 64. 65.
66. 67. 68.
Cheng, G.C. and R. Farmer, Real Fluid Modeling of Multiphase Flows in Liquid Rocket Engine Combustors. Journal of Propulsion and Power, 2006. 22(6): p. 1373-1381. Chen, Y., H. Shang, and P. Liaw, A fast algorithm for transient all-speed flows and finite-rate chemistry, in Space Programs and Technologies Conference. 1996, American Institute of Aeronautics and Astronautics. Shang, H.M., Numerical Studies of Spray Combustion in Liquid-Fueled Engines. 1992, University of Alabama in Huntsville. O'Rourke, P.J. and A.A. Amsden, The Tab Method for Numerical Calculation of Spray Droplet Breakup. 1987, SAE International. Hirschfelder, J.O., et al., Generalized Equation of State for Gases and Liquids. Industrial & Engineering Chemistry, 1958. 50(3): p. 375-385. Hirschfelder, J.O., et al., Generalized Thermodynamic Excess Functions for Gases and Liquids. Industrial & Engineering Chemistry, 1958. 50(3): p. 386-390. Edelman, R.B. and P.T. Harsha, Laminar and turbulent gas dynamics in combustors—current status. Progress in Energy and Combustion Science, 1978. 4(1): p. 1-62. W.C., J.G., Gas-Phase Combustion Chemistry. 2 ed. 2000: Springer-Verlag New York. XIII, 543. O'Rourke, P.J., Statistical properties and numerical implementation of a model for droplet dispersion in a turbulent gas. Journal of Computational Physics, 1989. 83(2): p. 345-360. Lempke, M., et al., Unsteady High-Order Simulation of a Liquid Oxygen/Gaseous Hydrogen Rocket Combustor. Journal of Propulsion and Power, 2015. 31(6): p. 1715-1726. Popp, M., et al., Liquid Rocket Thrust Chambers. Progress in Astronautics and Aeronautics. 2004: American Institute of Aeronautics and Astronautics. -1. Grisch, F., P. Bouchardy, and W. Clauss, CARS thermometry in high pressure rocket combustors. Aerospace Science and Technology, 2003. 7(4): p. 317-330. Singla, G., et al., Planar laser-induced fluorescence of OH in high-pressure cryogenic LOx/GH2 jet flames. Combustion and Flame, 2006. 144(1): p. 151-169. Magnussen, B., On the structure of turbulence and a generalized eddy dissipation concept for chemical reaction in turbulent flow, in 19th Aerospace Sciences Meeting. 1981, American Institute of Aeronautics and Astronautics. Gran, I.R. and B.F. Magnussen, A Numerical Study of a Bluff-body Stabilized Diffusion Flame. Part 1. Influence of Turbulence Modeling and Boundary Conditions. Combustion Science and Technology, 1996. 119(1-6): p. 171-190. Burke, M.P., et al., Comprehensive H2/O2 kinetic model for high-pressure combustion. 2012. 44(7): p. 444-474. Ó Conaire, M., et al., A comprehensive modeling study of hydrogen oxidation. 2004. 36(11): p. 603-622. Konnov, A.A., Remaining uncertainties in the kinetic mechanism of hydrogen combustion. Combustion and Flame, 2008. 152(4): p. 507-528. Alekseev, V.A., M. Christensen, and A.A. Konnov, The effect of temperature on the adiabatic burning velocities of diluted hydrogen flames: A kinetic study using an updated mechanism. Combustion and Flame, 2015. 162(5): p. 1884-1898. McGee, H.A., Molecular engineering. 1991: McGraw Hill. Rehab, H., E. Villermaux, and E.J. Hopfinger, Flow regimes of large-velocity-ratio coaxial jets. Journal of Fluid Mechanics, 1997. 345: p. 357-381. Vingert, L. and M. Habiballah, Test case RCM 2: Cryogenic spray combustion at 10 bar at Mascotte. Proceedings of the 2nd International Workshop on Rocket Combustion Modeling, 2001. Abramzon, B. and W.A. Sirignano, Droplet vaporization model for spray combustion calculations. International Journal of Heat and Mass Transfer, 1989. 32(9): p. 1605-1618. Daly, B.J. and F.H. Harlow, Transport Equations in Turbulence. 1970. 13(11): p. 2634-2649. Candel, S., et al., Experimental Investigation of Shear Coaxial Cryogenic Jet Flames. Journal of Propulsion and Power, 1998. 14(5): p. 826-834.
ACS Paragon Plus Environment
Page 31 of 31 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58
Energy & Fuels
31
69. 70. 71.
Habiballah*, M., et al., Experimental Studies of High-Pressure Cryogenic Flames on The Mascotte Facility. Combustion Science and Technology, 2006. 178(1-3): p. 101-128. Mardani, A. and E. Barani, Numerical Investigation of Supercritical Combustion of H2–O2. Energy & Fuels, 2018. 32(3): p. 3851-3868. Li, J., et al., An updated comprehensive kinetic model of hydrogen combustion. 2004. 36(10): p. 566-575.
ACS Paragon Plus Environment
