Abstract
The development of high-performance aeroengine combustion chambers strongly depends on the accuracy and reliability of efficient numerical models. In the present work, a reacting solver with a steady laminar flamelet model and spray model has been developed in OpenFOAM and the solver details are presented. The solver is firstly validated by Sandia/ETH-Zurich flames. Furthermore, it is used to simulate nonpremixed kerosene/air spray combustion in an aeroengine combustion chamber with the RANS method. A comparison with available experimental data shows good agreement and validates the capability of the new developed solver in OpenFOAM.
1. Introduction
Turbulent combustion modeling has been successfully used to facilitate the combustion chamber design of an aeroengine for many years now. However, during the development phase of a next-generation combustion chamber, the current CFD software is commonly used to understand the experimental phenomena and only provide qualitative guide. The prediction of the interaction between two-phase flow and combustion is one of the ultimate goals of the simulation of the combustion chamber. A lot of turbulent combustion models of varying complexity and accuracy have been proposed and applied to the aeroengine combustion simulation. An excellent overview of recent advances in a turbulent combustion model is given by Poinsot and Veynante [1]. And the advantages and potentials of different combustion models in aeroengine combustion chamber simulation are summarized by Mongia [2].
In this work, the open source code OpenFOAM [3, 4] is chosen as the computational platform due to its flexible framework and massive numerical libraries. Although OpenFOAM has many different solvers and physical models in the latest edition, it is still not widely used in a turbulent combustion area especially in engineering applications. This is because the solvers and models concerned on turbulent combustion in the OpenFOAM distribution are still far from enough compared to the commercial CFD software, such as ANSYS Fluent.
Various researchers have developed their individual combustion solvers based on OpenFOAM for their respective combustion cases recently. Marzouk and Huckaby [5] applied different kinetic models with a modified reacting solver reactingFoam to simulate the lab scale turbulent flame of syngas. The turbulence-chemistry interaction in this solver was modeled by the partially stirred reactor (PaSR) concept. Nogenmyr et al. [6] proposed a new version of reactingFoam based on an assumption of low Mach number. The new developed solver was addressed as reactingLMFoam and was used to simulate a jet flame together with large eddy simulations. Furthermore, Pang et al. [7] developed and validated a local time stepping-based PaSR solver called rareLTSFoam, which also includes the chemical reaction, the radiative heat transfer, and the soot formation. Apart from these, Lysenko et al. [8, 9] used the eddy dissipation concept (EDC) model for the simulations of the Sandia flame D and Sydney bluff-body flame HM1E. The radiation heat transfer was treated by the P1-approximation. Numerical results including temperature and composition distribution showed good accordance with the experimental measurements. Lilleberg et al. [10] incorporated a local extinction database into a new EDC solver and simulated the Sandia flames D and E. The local extinction approach improved the simulation results with a fairly small computational cost increase in comparison with the fast chemistry approach.
To reduce the chemistry efficiently, Ottino et al. [11] successfully applied the flamelet-generated manifold (FGM) technique with both RANS and LES methods resulting in high-quality solution for the turbulent combustion modeling. An improved accuracy is obtained in terms of the recirculation region and flame stabilization predictions. Mohd Yasin et al. [12] simulated a swirl-stabilized flame of biodiesel in OpenFOAM using real biodiesel physical properties. The spray structures with the reaction zone were obtained, and the distribution of droplets was analyzed. Bulat et al. [13] investigated a natural gas-air premixed swirl flame in an industrial gas turbine combustion chamber and well reproduced the experimental temperature and velocity distribution using the LES/PaSR model.
The flamelet approach has been the workhorse in the combustion simulation for large meshes and complex configurations of practical industrial applications due to its simplicity and low computational time. It is commonly used to fast predict the performance of the combustion chamber for gas turbines. However, it is not yet available among the combustion models in the OpenFOAM distribution package provided by OpenCFD Ltd. Muller et al. [14] have recently developed an extension edition which provided all steps necessary to use a steady laminar flamelet model (SLFM) with OpenFOAM. The solver of Muller was a transient solver combining LES with the SLFM and only limited to gas fuel combustion. Besides this, Cuoci et al. [15, 16] developed the libOpenSMOKE code with OpenFOAM. The flamelet concept was employed for the modeling of turbulent nonpremixed flame. Various soot formation models were also provided in the related library. In addition, Messig et al. [17] implemented the flamelet solver including several radiation models and simulated laminar diffusion flames.
In the current work, a new solver named SprayFlameletFoam is developed and validation is conducted through comparisons with experimental data. Our long-term goal is to develop LES methods with an advanced turbulent combustion model to simulate spray combustion encountered in high-temperature-rise aeroengine combustion chambers. Previously, this solver has been successfully used to simulate the afterburner flow field of a solid ducted rocket [18]. This paper mainly describes the implementation of the SLFM with a steady-state RANS method in OpenFOAM and its application to the simulation of spray combustion in an aeroengine. The need to perform spray combustion simulation motivated the implementation of a new SLFM OpenFOAM extension. In the next section, the formulation of the implemented SLFM is shown. The results of the CFD simulation of an aeroengine combustion chamber with a counter-rotating swirl cup are presented.
2. Model Formulation
In order to improve the reliability of the combustion simulation, a liquid fuel atomization should not be neglected in the solver. Therefore, the current solver SprayFlameletFoam is based on an existing solver simpleReactingParcelFoam. Similar to the most widely used solver reactingFoam, simpleReactingParcelFoam is also applicable for reacting flow but is a steady-state solver using a PaSR combustion model for compressible flow. The greatest advantage of choosing simpleReactingParcelFoam as the baseline solver is that a discrete phase Lagrangian model for liquid fuel and turbulent models are both taken into account by default. As a result, only the SLFM is introduced here.
The SLFM views turbulent flame as an ensemble of stretched laminar flames. Its main advantage is related to few variables which need to be solved in the CFD code, since it does not require solving a large number of species transport equations. The SLFM methods can be derived from the equations describing the counter flow flame as shown in Figure 1, by adopting a coordinate transformation from the physical space to the mixture fraction space. The mixture fraction is a conserved scalar, which significantly simplified the solution of species and temperature, particularly for nonpremixed flames. It is defined as the ratio of mass of material originating in fuel stream and mass of mixture.
A key assumption in the development of SLFM equations is that the Lewis number is unity. It means that the multicomponent diffusion coefficient equals the thermal diffusivity. The conservation equations of species and temperature in the mixture fraction space then can be written as where is the density, is the chemical species mass fraction, is the temperature, is the mixture specific heat at constant pressure, and and are the reaction rate and the specific enthalpy of specie , respectively. The influence of flow field on local flame structure is represented by the instantaneous scalar dissipation rate in the above equations. The value of is defined by the expression below.
For the counter flame configuration, the instantaneous stoichiometric scalar dissipation plays a role of a nonequilibrium parameter. The value of zero means the chemical equilibrium state, and as it increases due to aerodynamic straining, the nonequilibrium increases.
The flamelet library is built with the open-source chemistry software Cantera [19]. Starting with a low mass flow rate, flamelets can be calculated by gradually increasing the mass flow until the local quenching occurs.
Since the species mass fraction and temperature in the laminar flamelets are completely parameterized by the mixture fraction and the stoichiometric scalar dissipation , the mean species mass fractions and temperature in the turbulent flame brush can be determined from the joint PDF of and as where represents the species mass fractions or temperature.
Since and are considered statistically independent, the joint PDF can be simplified as the product of and . Assuming that fluctuations in are neglected, the PDF of is a delta function. For the mixture fraction, the shape of a presumed β-PDF is used.
Here, is the gamma function and and are the β-PDF parameters determining its shape and depending on the mean mixture fraction and its variance.
In order to obtain α and β, the transport equations of the mean mixture fraction and its variance are solved here.
The turbulent viscosity in the above equations is calculated with a k-ε turbulence model, and in the mixing fraction equation represents the source term due to the evaporation of the disperse phase. Both of them have already been contained in the simpleReactingParcelFoam solver and OpenFOAM libraries by default. The mean scalar dissipation rate in the SLFM can be calculated as follows [20]. where is only a constant of 2.0.
3. Validation of Sandia/ETH-Zurich Flames
The Sandia/ETH-Zurich CO/H2/N2 flames from the international workshop on measurement and computation of turbulent nonpremixed flames (TNF) are considered the validation cases. The species concentration and velocity measurements are accomplished at the Sandia National Laboratories and ETH Zurich, respectively [21]. It consists of two different, unconfined jet flames which are different in fuel nozzle geometry and fuel discharge velocity. They are named flames A and B, respectively. Both the two simple jet flames are employed to validate the SLFM used in this research. The velocity fields of the two flames are measured by LDV while temperature and species mass fractions are measured by CARS and Rayleigh scattering techniques.
The simulations are conducted on structured, nonuniform, and axisymmetric computational mesh, which is generated by the blockmesh utility in OpenFOAM. There are 160 cells along the axial direction and 64 cells along the radial direction. The geometric dimensions and boundary conditions are defined according to the experimental conditions and given in Table 1. The two flames have the same Reynolds number at the exit of the fuel nozzle. The synthesis fuel velocity is 76 m/s for flame A, compared to 45 m/s for the second case. The internal flow in the fuel nozzle is modeled to achieve the fully developed turbulent pipe flow.
According to the research of [22, 23], the realizable k-ε model shows better agreement with the experimental results than the standard model. Therefore, RANS coupled with the realizable k-ε model and SLFM is implemented. The combustion mechanism including 14 species and 34 elementary reactions proposed by Frassoldati et al. [24] is employed in the current research. Heat transfer by radiation is modelled using the P1 model, in which the absorption and emission coefficients of the mixture are obtained using a weighted sum of a grey gas model. Bounded Gauss limited linear schemes are used for the convective terms. The Gauss linear corrected difference schemes are used for diffusion terms and mass fluxes at face centers from cell values. This solver is based on a steady SIMPLE algorithm and has two pressure correctors and two momentum correctors per time step. The preconditioner of diagonal incomplete LU and the biconjugate gradient method (PBICG) by Jacobs [3] is used for solving momentum equations with a local accuracy of 10−8 at each step. The flow for the two cases reaches the steady-state conditions after approximately 4000 steps. All results are reported after 10000 steps to ensure that the reacting flow has totally reached the steady-state conditions.
A comparison between experimental and predicted temperature, mixing fraction, and axial velocity in an axial direction for two cases is shown in Figures 2 and 3. The three profiles show the same distribution tendency and similar peak value with the experimental results. Interestingly, the SLFM underpredicts all the profiles especially in the postflame zone for flame A. On the contrary, it overpredicts all the profiles for flame B. Possible explanations for the discrepancies include the turbulent model and inherent limitations of the SLFM.
(a) Temperature
(b) Mixing fraction
(c) Axial velocity
(a) Temperature
(b) Mixing fraction
(c) Axial velocity
To provide further validation, the radial profiles of axial velocity, H2O mass fraction, and temperature at 20 jet diameters for two cases are given in Figures 4 and 5. It can be seen that the predicted profiles of all the three parameters are nearly identical with the measurements for flame A. On the contrary, the predicted profiles shift slightly upward from the measurement for flame B. It indicates that the solver underpredicts the spreading rate of the fuel jet, causing insufficient mixing of fuel and coflow air in the downstream of the jet.
(a) Axial velocity
(b) H2O mass fraction
(c) Temperature
(a) Axial velocity
(b) H2O mass fraction
(c) Temperature
The two validation cases demonstrate that the new developed solver using the SLFM in OpenFOAM reasonably predicts the temperature and species mass fraction distribution.
4. Test Case Description
The test case of spray combustion in this paper is a typical single-head annular combustion chamber model, which is shown in Figure 6. The flow field through a diffuser and the annulus are simulated together with the flame tube, so we do not need to specify the boundary conditions for the swirlers and holes on the flame tube surface.
The configuration is symmetric except the swirl cup and has two quartz windows on the lateral sides. The flame tube contains primary and secondary dilution holes and eight film cooling strips on each side. The dome of the flame tube mainly consists of cooling holes, an inlet cowl, and a swirl cup, which is shown in Figure 7. The fuel is introduced from a simplex atomized injector which is located at the center of the swirl cup. The swirl cup provides two coaxial counter-swirling streams. Eight chamfered holes are used to form the inner swirling airstream, and an outer swirler including twenty curved vanes issues the counter-rotating airstream into the dome in the radial direction.
The swirl cup is widely used in aeroengines such as CFM56 [25]. In operation, the majority of droplets formed by the simplex fuel injector impinge on the Venturi, forming a thin fuel film on its surface. At the exit of the Venturi, the fuel film breaks into ligaments and further shatters the droplets due to a shear layer generated by the swirling streams.
In the SprayFlameletFoam solver, the primary breakup phenomena cannot be simulated with the discrete phase model, so the initial diameter and velocity of fuel droplets are directly specified using the pressure-swirl atomizer model according to the physical and working parameters of the fuel injector. All the droplets are released from the middle zone of the Venturi. In order to define the starting positions and velocities of the spray droplets, the properties of a pressure-swirl atomizer need to specify certain parameters including the temperature, mass flow rate, drag coefficient, spray angle, and type of droplet size distribution, as shown in Table 2.
In the current simulation, the following Rosin Rammler equation is applied to describe the distribution of the droplet size. where is the fraction of the total number of the droplets whose diameter is less than . The scale factor is set to 3.0, and parameter is set to 74 μm. For a smooth distribution of the droplets, 10° dispersion angle is defined to give a hollow spray cone.
Only the liquid droplet evaporation and the breakup are considered for the current simulation, and coalescence of fuel droplets is not considered here. The ReitzDiwakar breakup model is included in simpleReactingParcelFoam by default.
An experimental study of the combustion chamber model is also conducted in our research group. The temperature near the primary holes was measured by means of the coherent anti-Stokes Raman scattering (CARS) microscopic imaging technique. The predicted temperature distribution will be compared with measurement data.
The geometry of the simulation zone is shown in Figure 8. The grids are generated by commercial software ICEM. Tetrahedral grids are generated in/around the dome and the flame tube. The other parts are generated in hexahedral grids. The scale of the grid is about 1 mm, and the total number of the grids is about 2.7 million. The mesh is refined inside the swirl cup and the dome while it is gradually coarsened towards the chamber outlet. The computational grid on the surface of the swirl cup and the flame tube is shown in Figures 9 and 10. The origin of the coordinates is placed at the center of the fuel injector orifice. The axis of the combustion chamber is along the x-axis while the y- and z-axes are along the vertical and horizontal directions, respectively.
Boundary conditions in OpenFOAM are split into boundary faces in a set of patches [3]. The type of these patches is specified using certain keywords in OpenFOAM. For each independent property field, a file exists where initial numerical conditions are specified. In the present case, an airflow rate of 0.44 kg/s is introduced through the air inlet of the diffuser. The initial pressure field is 5 bars for the air inlet boundary. Air inlet temperature is set to 870 K. The value of the turbulent intensity is 0.05 at the air inlet. A uniform velocity profile is specified for the air inlet, and the velocity is parallel to the x-axis. The other two velocity components are set to be zero. The zero-gradient pressure is imposed at the outlet boundary. A nonslipping condition for velocity is applied to all the walls. The temperature at the walls is calculated according to the zero-gradient boundary condition. The positions and diameter of fuel droplets are defined according to the pressure-swirl atomizer model contained in the OpenFOAM platform and released at the middle zone of the Venturi. In the present study, the cold-flow field is first simulated and used as the initial condition of the combustion case. The gas phase is described using the Eulerian formulation, while the liquid fuel droplets are tracked by the Lagrangian formulation. The 16-species 23-step kerosene (C12H23) reaction developed by Kundu et al. [26] is used in the study.
5. Results and Discussion
5.1. Analysis of Velocity Fields
The streamlines and axial velocity contours at the combustion chamber center plane are given in Figure 11. The incoming air is uniform and gradually diverges in the diffuser. Then, it gradually enters the flame tube of the combustion chamber from the swirl cup, dilution holes, and cooling holes along the axial direction. Two small vortices form in the corner region of the annulus passage between the combustion casing and the flame tube.
(a) The streamlines in the whole combustion chamber
(b) The streamlines in the primary zone
(c) Axial velocity contours
The flow in the flame tube is complex and three dimensional with high degrees of swirl. The full flow field in the flame tube may be divided into 3 regions: downstream of the swirl cup, near field of the primary holes, and fully developed regions. The swirling jets expand abruptly downstream of the swirl cup and impinge on the flame tube. The flow departing from the swirl cup is attached to the wall due to the Coanda effect [27]. The gas phase streamlines in the flame tube show that two vortices form in the primary zone and extend downstream of the primary holes. The black line in Figure 11(c) is the boundary of the center recirculation zone, which means that the axial velocity equals zero on the line. The existence of the central recirculation zone is the main feature of the primary zone. It plays a vital role of anchoring the flame of the combustion chamber. The position and shape of the recirculation zone are the most significant aspects of combustion chamber design and may vary with both physical and operating conditions. The high-temperature combustion products recirculate upstream in the center recirculation zone. The local low velocity in the vicinity of the recirculation zone boundary matches the flame propagation speed. At the same time, the recirculation zone increases the fuel air mixture residence times and enhances the combustion stability, thus increasing combustion efficiency.
The flow field near the primary holes is the secondary region, where the average flow is not uniform in different vertical positions. The primary jets terminate the center recirculation zone and impinge each other at the centerline. The last region covers the fully developed flow and extends up to the combustion chamber outlet. The presence of the flame tube walls causes the air in the dilution zone of acceleration. The airflow progressively becomes uniform and parallel to the axial direction.
The phenomenon of the counter-swirling flow is observed at the downstream of the swirl cup. Figure 12 shows the streamlines and axial velocity contour for the plane cutting through the primary swirler and the secondary swirler, respectively. The primary swirler produces the clockwise flow while the secondary swirler generates the counter-swirling flow. However, at the downstream of the swirl cup, the secondary swirl dominates the flow field.
(a) Primary swirler exit
(b) Secondary swirler exit
5.2. Analysis of Temperature Contours
The predicted temperature distribution and CARS measurement positions in the experiment are shown in Figure 13. The new developed solver SprayFlameletFoam successfully predicts the position of a high-temperature zone in the primary zone. The temperature shows a steep rise from the swirl cup exit due to the intensive reaction in the combustion zone. At the downstream of the swirl cup, temperature increases gradually from 870 K to a maximum value of 2400 K in the central region due to efficient combustion. The temperature in the central recirculation zone is about 2000 K. Upstream of the primary air jets, the minimum temperature exists in the vicinity of the wall while the maximum value exists near the centerline.
(a) OpenFOAM results
(b) CARS measurement positions
It should be noted that the asymmetric behavior of the flame with high temperatures in the primary zone is also found during the experiments. With the introduction of primary air jets, the temperature falls abruptly and increases again behind the primary air jets. After the secondary dilution holes, the gas diffuses away from the centerline to achieve entrainment of more air along the centerline of the flame tube, thus reducing the temperature in the dilution zone.
A comparison of the predicted and experimental temperatures is displayed in Table 3. Due to the measuring uncertainty of the scanning CARS system, the margin of error is below 4%. The difference between the measurement and simulation is acceptable. The main features such as the location of the flame zone as well as the highest temperature value are in good agreement with the experimental findings. The high temperature value emerges in the edge of the central recirculation zone and the recirculation zones behind the primary jets such as positions D, E, and H. The distribution tendency of SprayFlameletFoam is in agreement with the CARS measurement.
Figure 14 shows the temperature contour at two selected planes. Temperature contour at the plane cutting through the primary holes depicts an annular high-temperature zone in the outer region and an ellipse high-temperature zone in the inner region. The annular high-temperature zone is formed due to the reaction in the shear layer between the swirling jets and the recirculating zone. Temperature contour at the plane cutting through the secondary holes reveals that the reaction rate also increases at the downstream of the primary jet, causing the secondary high-temperature zone due to better mixing of air and fuel. The temperature distribution becomes more uniform in the central region, which may result in better match with the work of nozzle guide vanes.
(a) The plane cutting through primary holes
(b) The plane cutting through secondary holes
5.3. Analysis of Species Concentration
The mass fraction of different species at the center plane is analyzed. Figures 15(a)–15(d) show the contour of mass fraction of C12H23, CO2, H2O, and OH, respectively. Most of the fuel distributes in the shear layer of the swirl streams and decreases drastically downstream due to efficient mixing and combustion with air. Then, the fuel is consumed in the initial region of the flame tube. The concentration of CO2 is less than 0.12 but increases along the axial direction. The mass fraction contour of CO2 is almost symmetric, and most are generated behind the primary jets due to sufficient air addition. The mass fraction of H2O distribution is just similar to that of CO2 as it is also the final product of hydrocarbon fuel combustion. Difference is only in magnitude while plenty of H2O is formed in the recirculation zone. Figure 15(d) is particularly helpful for determining the location of the flame zone since the OH-to-H2O conversion is the major heat release reaction in the combustion process [28]. The OH is largely formed in a hook-shaped region attached to the inner edge of the swirl cup exit.
(a) The mass fraction of C12H23
(b) The mass fraction of CO2
(c) The mass fraction of H2O
(d) The mass fraction of OH
5.4. Droplet/Turbulence Interaction
The distribution of the spray droplets influences the heat release rate and exhaust gas emission significantly [25]. Figure 16 shows the droplet diameter distribution throughout the spray field within the flame tube compared with the experimental results from the laser hologram measurement. The gas phase is almost unaffected by the presence of the droplets. The complex flow field generated by the swirl cup leads to a much complicated spray structure. The spray is of a hollow cone structure with very few droplets found at the centerline. The droplets are found in large numbers in the gas phase shear layer, where the droplets break up, evaporate, and mix with the air. Smaller droplets may be recirculated by the gas phase. Instead, the larger droplets are not recirculated and flow away from the centerline.
(a) OF results
(b) Experimental results
In order to provide the insight into the droplet spatial distribution, the relative number of droplet distribution from the OF results and experiment is shown in Figure 17. The maximum diameter of fuel droplets is less than 150 μm under the operating conditions in the current work. In the measurement, the droplets of a diameter less than 30 μm have not been distinguished. Nevertheless, the distribution of the relative number of droplets is in good agreement. The diameter corresponding to the most relative number of droplets from the OF results is 84 μm, which is slightly larger than that from the experimental results.
(a) OF results
(b) Experimental results
6. Conclusions
(1)A new reacting solver named SprayFlameletFoam using the steady laminar flamelet model (SLFM) together with the spray model is developed based on the existing reacting solver simpleReactingParcelFoam in OpenFOAM. The validation of Sandia/ETH-Zurich flames demonstrates that the new solver reasonably predicted the temperature and species mass fraction distribution.(2)The SprayFlameletFoam solver is then used for modeling spray combustion of an aeroengine combustion chamber. The comparison with experimental measurements confirms the feasibility and capability of the SprayFlameletFoam solver for the industrial applications. The current work shows promise for using OpenFOAM to predict combustion chamber performance and provide insight into the structure of the flow and concentration fields. More validation cases with engineering applications will be performed in the future.
Nomenclature
| : | Density (kg/m3) |
| : | Mass fraction of the chemical species |
| : | The temperature (K) |
| : | The mixture specific heat at constant pressure (kJ/kg·K) |
| : | The reaction rate of specie (kg/m3·s) |
| : | The specific enthalpy of specie (kJ/kg) |
| : | Scalar dissipation rate |
| : | Mixture fraction |
| : | Stoichiometric scalar dissipation |
| : | Source term due to the evaporation |
| : | Dynamic viscosity (Pa·s) |
| : | Turbulence kinetic energy (m2/s2) |
| ε: | Turbulence dissipation rate (m2/s3) |
| : | Molecular diffusion coefficient (m2/s) |
| PaSR: | Partially stirred reactor |
| CFD: | Computational fluid dynamics |
| EDC: | Eddy dissipation concept |
| FGM: | Flamelet-generated manifolds |
| SLFM: | Steady laminar flamelet model |
| TNF: | Turbulent nonpremixed flames |
| CARS: | Coherent anti-Stokes Raman scattering |
| OF: | OpenFOAM. |
Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.
Acknowledgments
This work has been funded by the National Natural Science Foundation of China under Grant 51576164.