Large Eddy Simulation of Kerosene/LOx Supercritical Mixing
5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS) Large Eddy Simulation of Kerosene/LOx Supercritical Mixing Characteristics of a Swirl Injector using Various Cubic Equation of State Junyoung Heo*, Kukjin Kim*, and Hong-Gye Sung* *Korea Aerospace University Goyang Gyeonggi, 412-791, South Korea Abstract The turbulent mixing of a kerosene and liquid oxygen coaxial swirl injector under supercritical pressures has been numerically investigated in three dimensional system. Turbulent numerical model is based on LES (Large Eddy Simulation) with real-fluid transport and thermodynamics over the entire pressure range; Cubic Equations of State (SRK, PR, and RK-PR), Chung’s model for viscosity/conductivity, and Fuller’s theorem for diffusivity to take account Takahashi’s compressible effect. The effects of real-gas equation of state and swirl intensity on flow structures in the swirl injector are investigated and power spectral densities of pressure fluctuations and mixing efficiency of the injector are analyzed. Nomenclature D E f p q t u x Z = = = = = = = = = Diffusivity Specific total energy Mixture fraction Pressure Heat flux Time Velocity components Spatial coordinate Pseudo-time variable vector = = = = Kronecker delta Preconditioning matrix or circulation Density Viscous stress tensor Greek Superscripts sgs = Resolved-scale = Favre-averaged resolved-scale = Subgrid-scale 1. Introduction An advanced liquid rocket engine is operating at high pressure in order to increase combustion efficiency and thrust, thus the propellants are exposed to supercritical conditions in the combustion chamber. The supercritical conditions are the state that arbitrary substance is enclosed by the environment over its critical point. In such Copyright 2013 by Heo, J. Y., Kim, K. J., and Sung, H. G. Published by the EUCASS association with permission. Heo, J. Y., Kim, K. J., and Sung, H. G. conditions, the propellants experience the transition of thermophysical properties different from phase change process at subcritical conditions. After the liquid propellant is introduced into the injector and combustion chamber, its properties are varied continuously and the distinction between gas and liquid is ambiguous. At near the critical point, propellant properties exhibit gas-like diffusivities and liquid-like densities so that the engine performance can be achieved with high energy density. Hence, the study for the mixing and combustion process of the propellants at supercritical conditions is an essential subject. Because of the peculiarities of substances at supercritical conditions, the conventional methods for property analysis are inappropriate so that the applicable thermodynamic relationships for a wide range of pressure must be employed. For understanding of the atomization, mixing, and combustion phenomena in advanced liquid rocket engines, experimental and numerical studies are performed in the wide range of pressure and temperature including the supercritical conditions. In research of numerical analysis, the mixing and combustion processes in shear coaxial injectors using propellants as methane/liquid oxygen and hydrogen/liquid oxygen at supercritical conditions have been investigated. The thermophysical characteristics of hydrogen and liquid oxygen flame in a shear coaxial injector at supercritical pressure have been investigated by Oefelein using direct numerical simulation [1]. Zong and Yang applied the large eddy simulation for the identification of cryogenic fluid dynamics in a swirl oxidizer injector operating at supercritical pressure [2]. Giorgi and Leuzzi performed the numerical analysis of mixing and combustion in methane and liquid oxygen shear coaxial injector at supercritical conditions [3]. Recently, Giorgi et al. simulated combustion characteristics of methane and liquid oxygen and compared between the results in accordance with combustion models, reaction mechanisms [4], and equations of state at supercritical conditions. Sierra et al. numerically investigated reaction mechanisms and combustion characteristics for the recess number using hydrogen and liquid oxygen as propellants [5]. Hiroshi et al. analysed nitrogen injection characteristics and compared between the numerical results to pressure and temperature in combustor [6]. H. Huo et al. simulated combustion characteristics in a kerosene and liquid oxygen supercritical injector applying the flamelet library technique as a combustion model [7]. In this study, on the basis of previous study [8], the turbulent mixing of a kerosene and liquid oxygen coaxial swirl injector under supercritical pressures has been numerically investigated in three dimensional system. Turbulent numerical model is based on LES with real-fluid transport and thermodynamics; SRK (Soave modification of Redlich-Kwong) [9], PR (Peng-Robinson) [10], RK-PR (Redlich-Kwong-Peng-Robinson) [11] equations of state, Chung’s model for viscosity/conductivity, and Fuller’s theorem for diffusivity to take account Takahashi’s compressible effect. The effect of real-gas equation of state and swirl intensity are investigated and power spectral densities of pressure fluctuations and mixing efficiency are analyzed. 2. Numerical Method 2.1 Filtered Transport Equations The theoretical formulation is based on the filtered Favre averaged mass, momentum, energy, and mixture fraction conservation equations in Cartesian coordinates. Turbulent closure is achieved by means of a LES technique in which large scale turbulent structures are directly computed and the unresolved small scale structures are treated by using the analytic or empirical modeling. The governing equations can be written as: u j 0 t x j (1) sgs sgs ui ( ui u j p ij ) ( ij ij Dij ) t x j x j (2) E ( E p)ui qi u j ij Qisgs Hisgs ijsgs t xi xi (3) f ( ui f ) f isgs D t xi xi xi (4) The unclosed sub-grid scale terms are defined 2 5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS) ij sgs ui u j ui u j (5) Dij sgs ij ijj (6) Qi sgs qij qijj (7) Hi sgs Eui Eui pui pui (8) i sgs u j ij u j ijj (9) (10) i sgs ui f ui f ij sgs , Dij sgs , Qij sgs , H ij sgs , ij sgs and ij sgs are the sgs stress, nonlinearity of viscous stress term, heat flux, energy flux, viscous work, and conserved scalar flux, respectively. The dynamic Smagorinsky model is employed to close those sgs terms [12]. 2.2 General-Fluid Thermodynamics and Transport Intermolecular forces and the molecular body volume are considered as important factors in the real-fluid. In order to understand the real-fluid phenomenon, effects such as compressibility, van der Waals force, molecular dissociation, and non-equilibrium thermodynamic quantity must be taken into account. An equation of state capable of handling dynamics of real-fluid at high pressure and low temperature regime is required for thermophysical property evaluation of real-fluid. The most common equations of state utilized for calculating high pressure fluid properties are the cubic equations of state, which include SRK, PR, and RK-PR equations of state. SRK equation of state: p Ru T a 2 M w M w b (11) a 2 M w b M w b 2b M w (12) a 2 M w 1b M w 2b (13) M w b PR equation of state: p RuT M w b RK-PR equation of state: p Ru T M w b 2.3 Numerical Scheme The governing equations are numerically solved by means of a finite volume method. This method allows for the treatment of arbitrary geometry. The spatial discretization employs a fourth order central differencing scheme in generalized coordinates. The temporal discretization is obtained by second order backward differencing scheme using a fourth order Runge-Kutta scheme for integration of the real time term. In the regions of low Mach number flows, the energy and momentum equations are practically decoupled and the system of conservation equations becomes stiff. Pressure decomposition and preconditioning techniques with dual time stepping are applied to circumvent the round-off errors associated with the calculation of the pressure gradient and the convergence difficulties for the low Mach number flows in the momentum equation. First, a rescaled pressure term is used in the momentum equation to circumvent the singular behavior of pressure at low Mach numbers. Second, a 3 Heo, J. Y., Kim, K. J., and Sung, H. G. preconditioning technique with dual time-stepping integration procedure is established. A unified treatment of general fluid thermodynamics is incorporated into a preconditioning scheme. The pseudo-time derivative may be chosen to optimize the convergence of the inner iterations through the use of an appropriate preconditioning matrix that is tuned to rescale the eigenvalues to render the same order of magnitude to maximize convergence. To unify the conserved flux variables, a pseudo-time derivative of the form Z / can be added to the conservation equation. Since the pseudo-time derivative term disappears as the solution is converged, a certain amount of liberty can be taken in choosing the variable, Z. We take advantage of this by introducing a pressure p as the pseudo-time derivative term in the continuity equation. The code is paralleled using an MPI library for more effective calculation. 3. Model Description Figure 1 shows a coaxial swirl injector geometry and the swirl numbers of liquid oxygen and kerosene are 1.598 and 6.548, respectively. The injector includes three major parts: a tangential inlet, a vortex chamber, and a discharge nozzle. The kerosene (outer) and liquid oxygen (inner) injects into the swirl injector through the tangential passage. In case of quasi-3D simulation, the radial and azimuthal velocities are determined from the tangential inlet port angle and these velocity components are factor to set the swirl strength and mass flow rate, respectively. LOx Fuel dt,f or dt,o Lr dn,o So 4d s , o d n , o nn,o dt ,o 2 1.598 , doo Sf dn,f 4 d n, f do,o d s , f n f dt , f 2 ds,f or ds,o 6.548 Figure 1: Injector geometry and swirl number Table 1: Operation conditions Injector Chamber Pressure Fuel Oxidizer Mass Flow Rate Fuel Coaxial Swirl 100 bar Kerosene; 350 K Liquid Oxygen; 103 K 0.084 kg/s Oxidizer 0.232 kg/s Table 2: Computation conditions Case 1 2 3 4 5 4 Equation of State SRK PR RK-PR RK-PR RK-PR Dimension Quasi-3D Quasi-3D Quasi-3D Quasi-3D Full 3D Tangential Velocity x1 x1 x1 x2 Determined by Inlet Geometry 5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS) Figure 2: Computational domain (Quasi-3D) Figure 3: Computational domain (Full 3D) Operation and computation conditions are listed in Table 1 and 2, respectively. Quasi-3D flow fields with periodic boundaries is treated herein because of the enormous computational time for simulating of full 3D using the large eddy simulation. The tangential inlets which are circular ports connected to the injector are approximated with a thin slit. Figure 2 shows a computational domain divided into 45 blocks, with each calculated on a single processor of a distributed computing facility. The total grid points are about 765 000 nodes. The computational domain of full 3D calculation is shown as Fig. 3 and grid points are 6 561 000 nodes which is divided into 146 blocks. Fuel and Oxygen are supplied through the projected faces like Fig. 3. The disturbances at all of calculation conditions are generated by a Gaussian random-number generator with an intensity of 5% of the mean quantity. 4. Results and Discussion SRK and PR equations of state usually have been used to predict the thermophysical properties of real-gas. The SRK equation of state has high accuracy for low molecular weight compounds while the PR equation of state does for high weight. Considering kerosene and liquid oxygen used as propellants, the SRK equation of state can predict the supercritical oxygen properties properly and the PR equation of state is accurate for kerosene properties at supercritical conditions. However, the RK-PR equation of state can be properly applied to simulate the properties of both kerosene and liquid oxygen [13-14]. The fuel and oxidizer having the only tangential and azimuthal velocity components are supplied for the quasi-3D calculation of the swirl coaxial injector. In that calculation, the tangential and azimuthal velocity determine the swirl strength and mass flow rate, respectively. The swirl strength is controlled by changing the tangential velocity in this study. 5 Heo, J. Y., Kim, K. J., and Sung, H. G. 4.1 The Effect of Cubic Equation of State (a) Kerosene (b) Oxygen Figure 4: Density distributions for different equation of state at 100 bar (a) SRK (b) PR (c) RK-PR Figure 5: Instantaneous density contours according to real-gas equation of state Figure 4 represents density distributions to different equation of state as compared to that of the principle of ECS (Extended Corresponding States) at 100 bar. The principal of ECS is that the properties of any interest fluid are represented by scaling that of the reference fluid and may also be functions of density and temperature. The principle of ECS is a very accurate method for prediction of thermophysical properties, applicable to the entire range of fluid states, from dense liquid to dilute gas, as well as to the supercritical fluid conditions [15-17]. As shown in the results, the RK-PR equation of state provides the most suitable results for kerosene and oxygen in the wide range of temperature. Figure 5 shows comparison results of density fields to the defferent equations of state with the same inlet conditions of fuel and oxidizer. The SRK and PR equations of state predict lower kerosene density and higher oxygen density respectively. On the other hand, the RK-PR equation of state has very high accuracy in the prediction results of both kerosene and liquid oxygen. This study also shows that the SRK and PR equations of state have lower density for kerosene and higher density for liquid oxygen, respectively, than the results by the RK-PR equation of state. As a result, it is investigated that the liquid oxygen core sustains its length as long and the CTRZ (Central Toroidal Recirculation Zone) is formed in the rear region compared with the other results in the results applied the PR equation of state due to the high density of liquid oxygen. Also, the mixing layer of the fuel/oxidizer gets longer and the CTRZ possesses the wider area since the kerosene density results by the RK-PR equations of state are higher than those by the SRK equation of state. 6 5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS) 4.2 Swirl Intensity Effect Two tangential velocities are applied to investigate the swirl strength effect on mixing characteristics. Figure 6 represents the density fields for each swirl strength. The formed liquid oxygen film is closer to the injector wall since a centrifugal force of swirl flow becomes stronger as the tangential velocity increases. As the swirl strength is stronger, the gas phase core is formed in the central region of oxidizer injector which affects the film boundary area enlarged in the injector and facilitates vorticity breakup. Figure 6: Instantaneous density contour according to swirl intensity Figure 7: Mixing efficiency at different swirl intensity To quantify the level of mixing, mixing efficiency is considered using the following formula: mix Y ( x) Y Fuel u dA Fuel udA 1 1 / global local ( local global ) ( local global ) (14) 7 Heo, J. Y., Kim, K. J., and Sung, H. G. Where ρ and YFuel are density and mass fraction of kerosene, respectively. The u is the velocity component normal to the area, dA , local and global are local and global equivalence ratio, respectively. The parameter is decided by local equivalence ratio. This parameter means that if the local equivalence ratio is equal or less than global equivalence ratio, the mixing is regarded as completed. On the other hand, local efficiency is less than 1 when local equivalence ratio is large. Therefore, the mixing efficiency represents how much fuel spread into combustion chamber, in other words, the efficiency becomes 1 when local equivalence ratio is equal to global equivalence ratio through whole cross sectional area. Figure 7 shows mixing efficiencies on the time averaged fields. The efficiency increases as moves to downstream after the oxidizer injector post where the mixing starts. As shown in Fig. 7, the stronger swirl intensity is, the earlier the mixing efficiency reaches 1, meaning that most regions become lean burn environment because the mixing completely takes place not much far from the oxidizer injector post. Then, when increasing the swirl strength, the film is formed inside the oxidizer injector by the strong swirl flow, rapid mixing of fuel and the oxidizer is performed. Figure 8: Probe location Figures 8 shows the locations of pressure measurement: inner side of injector (probe 1, 2), near the oxidizer injector post (probe 3), phase-change region (probe 4), mixing layer (probe 5), and front of the CTRZ (probe 6). Figure 9 represents the power spectral densities of the pressure fluctuations at the locations. The Fast Fourier Transform (FFT) technique was implemented for the spectral analysis. (a) Case 3 (b) Case 4 Figure 9: Power spectral densities of pressure fluctuations according to swirl intensity As shown in case 3, a frequency of 1.05 kHz is dominantly estimated inside the oxidizer injector filled with high density fluid and its amplitude is higher than the other locations. Frequencies of 1.05 and 2.75 kHz are observed at the other locations too. The oscillation at 2.75 kHz is same as a cycle of the vortex shedding at the oxidizer injector exit. Also, high amplitude of the pressure oscillation at 2.75 kHz is observed in kerosene injector due to the Doppler 8 5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS) effects. In case 4, a frequency of 2.23 kHz which is same as a cycle of the vortex shedding at the tip of oxidizer injector post is dominantly estimated for all locations while the low frequency observed inside the oxidizer injector did not occurred as compared with those of case 3. The initial turbulence flow at the tangential inlet generates a succession of wavy structures and is convected downstream through the dense liquid oxygen. In a swirl injector, those waves are generated by the presence of oscillations in the chamber or injector post. When the swirl intensity is strong, the core boundary (film) is formed in propellant feed-line. Also, gasification of the liquid oxygen is promoted from the core boundary and the mixing is rapidly performed in order of the vortex formed at the interface of the liquid phase. The pressure wave of the low frequency is not generated inside oxidizer injector due to the film. 4.3 3D Effect The swirl injector flow structures computed in full 3D are presented in Fig. 10. Liquid oxygen and kerosene are introduced into the injector through slots. Since the inlet ports consist of the tangential slit, the strong radial pressure gradient is formed by centrifugal force. Because the chamber pressure is higher than the critical pressure, there is no obvious boundary between the injected fluid and the ambience. Near the oxidizer injector post, a slowly swirling gaseous region is formed. The adverse pressure gradient in this region gives rise to recirculation flow and finally produces a central toroidal recirculation zone (CTRZ), resulting from vortex breakdown, which plays an important role in determining the flame stabilization characteristics. The CTRZ accelerates the earlier burst of vortices near the fuel/oxygen mixing zone. In the CTRZ, these vortical structures move downstream and upstream, and are rapidly dissipated by the turbulent diffusion and viscous damping. In the full 3D results, unlike the quasi-3D calculation, the core of the gas phase is formed in the oxidizer injector. Pressure oscillations in the chamber are reflected in the gas phase region inside injector interacted acoustically with the chamber, which deforms the boundary of the liquid and gaseous phases. As shown in Fig. 10, gasification of the liquid oxygen is promoted by the pressure wave transmitted to the inside of the injector, the dynamic flow of the injection is amplified in the formation of vortices at the film boundary. It may cause a significant effect on the mixing structure in the chamber. Figure 10: Flow structure of coaxial swirl injector Figure 11 represents the temporal evolution of pressure and temperature field over one cycle of the dominant pressure oscillation in oxidizer injector. As the propellants are discharged into the chamber, sinuous KelvinHelmholtz waves develop along the outer wall of the injector. As these waves move downstream, a hairpin vortice is generated and grows up. The oscillation is same as a cycle of the vortex shedding caused by the hairpin vortices at the LOx post. The pressure fluctuations from rotating in the tangential direction by the full 3D effect, these pressure waves affect the change of film thickness in the tangential direction. Similar to the quasi-3D results, fuel - and oxidizer injectors 9 Heo, J. Y., Kim, K. J., and Sung, H. G. have a dominant frequency about 1.4 kHz and 3.4 kHz, respectively, i.e. the dominant frequency in the fuel injector is higher than a frequency in the oxidizer injector due to the thin flow passage. θ=0˚ θ = 120 ˚ θ = 240 ˚ (a) Pressure (b) Temperature Figure 10: Temporal evolution of pressure and temperature field over one cycle of the dominant pressure oscillation in the oxygen injector 5. Conclusions The mixing characteristics of co-axial swirl injector are investigated according to real-gas equation of state and swirl intensity under supercritical condition. The RK-PR equation of state provides more reasonable flow and mixing structures than the SRK, PR equation of state. In order to compare the mixing efficiency according to the swirl intensity, inlet conditions are controlled by different tangential velocities. The case increasing tangential velocity up to twice reaches mixing efficiency 1 about 12% ahead of a baseline case. As the swirl intensity increases, the gaseous core is obviously formed inside oxidizer injector. Also, gasification of the liquid oxygen is facilitated from the core boundary so that the low frequency, 1.05~1.4 kHz, inside oxidizer injector disappears due to the active oscillation on the film boundary. Full 3D simulation has been calculated in order to apply accurate tangential inlet ports. The gasphase core is formed in the oxidizer injector unlike the quasi-3D calculation. The vertical structures are amplified at the core interface, causing a significant effect on the mixing structure in the chamber. 10 5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS) Acknowledgement This research was supported by National Space Laboratory (NSL) Program (No. 2008-2006287) through the National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology. References [1] Oefelein, J. C. 2005. Thermophysical Characteristics of Shear Coaxial LOx-H2 Flames at Supercritical Pressure. Proceedings of the Combustion Institute. Vol. 30. Issue 2. 2929-2937. [2] Zong, N. and Yang, V. 2008. Cryogenic Fluid Dynamics of Pressure Swirl Injectors at Supercritical Conditions. Physics of Fluids. Vol. 20. Issue 5. 056103-1-14. [3] Giorgi, M. G. D. and Leuzzi, A. 2009. CFD Simulation of Mixing and Combustion in LOx/CH4 Spray Under Supercritical Conditions. 39th AIAA Fluid Dynamics Conference. [4] De Giorgi, M., G., Sciolti, A., Ficarella, A. 2011. Comparisons between Different Combustion Models for High Pressure LOx/CH4 Jet Flame. AIAA Fluid Dynamics Conference and Exhibit. [5] Sierra, P., Masquelet, M., Menon, S. 2013. Large-Eddy Simulation of a Reactive Shear Coaxial Injector Configuration. AIAA Aerospace Sciences Meeting. [6] Hiroshi, T., Mitsuo, K. 2013. Characterization of Cryogenic Nitrogen Jet Mixings under Supercritical Pressures. AIAA Aerospace Sciences Meeting. [7] Huo, H., Yang, V. 2013. Large-Eddy Simulation of Supercritical Combustion of Liquid Oxygen and Kerosene of a Bi-Swirl Injector. AIAA Aerospace Sciences Meeting. [8] Heo, J. Y., Kim, K. J., Sung, H. G. 2011. Numerical Study on Kerosene/LOx Supercritical Mixing Characteristics of a Swirl Injector. AIAA Aerospace Sciences Meeting. [9] Soave, G. 1972. Equilibrium Constants from a Modified Redlich-Kwong Equation of State. Chemical Engineering Science. Vol. 27. Issue 6. 1197-1203. [10] Peng, D. and Robinson D. B. 1976. A New Two-Constant Equation of State. Industry Engineering Chemical Fundamentals. Vol. 15. No. 1. 59-64. [11] Cismondi, M. and Mollerup, J. 2005. Development and Application of a Three-Parameter RK-PR Equation of State. Fluid Phase Equilibria. Vol. 232. Issue 1-2. 74-89. [12] Germano, M. 1991. A Dynamic Subgrid-scale Eddy Viscosity Model. Phys. Fluids A. Vol. 3. No. 7. 1760-1765. [13] Kim, S. K., Choi, H. S., Kim, Y. M. 2012. Thermodynamic Modeling Based on a Generalized Cubic Equation of State for Kerosene/LOx Rocket Combustion. Combustion and Flame. Vol. 159. Issue 3. 1351-1365. [14] Kim, K. J., Heo, J. Y., Sung, H. G. 2012. Study on Thermophysical Property Characteristics of Kerosene Surrogates in a Swirl Injector at Supercritical Pressure Condition. The Korean Society of Propulsion Engineers Spring Conference. [15] Ely, J. F. and Hanley, H. J. M. 1981. Prediction of Transport Properties. 1. Viscosity of Fluids and Mixtures. Industrial and Engineering Chemistry Fundamentals. Vol. 20. No. 4. 323-332. [16] Ely, J. F. and Hanley, H. J. M. 1983. Prediction of Transport Properties. 2. Thermal Conductivity of Pure Fluids and Mixtures. Industrial and Engineering Chemistry Fundamentals. Vol. 22. No. 1. 90-97. [17] Yang, J. C., Vazquez, I., Boyer, C. I., Huber, M. L., and Weber, L. 1997. Measured and Predicted Thermodynamic Properties of Selected Halon Alternative/Nitrogen Mixtures. International Journal of Refrigeration. Vol. 20. No. 2. 96-105. 11
© Copyright 2024