{'fraction_non_alphanumeric': 0.03145611364789447, 'fraction_numerical': 0.016344132782488946, 'mean_word_length': 4.809473684210526, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Investigation of temporal variation of combustion instability intensity in a back step combustor using LES

Lean combustion is susceptible to combustion instabilities, especially pressure oscillation, and this severely damages the combustor. To prevent this, many studies have attempted to elucidate this phenomenon. However, the physical understanding of the mechanism is still incomplete owing to its complexity. In this study, combustion instability in turbulent spray combustion in lean fuel condition is investigated using Large Eddy Simulation (LES). The combustion chamber has a backward facing step, and kerosene fuel droplets are injected vertically just upstream from the step. The dynamic thickened flame model is used as the turbulent combustion model, and a two-step global reaction model is used for kerosene-air flames. The influence of pressure oscillation on fuel atomization behavior is considered using a model, which predicts the Sauter Mean Diameter of injected the fuel droplets with time. The equivalence ratio is set as 0.6, 0.8, and 1.0. The results show that combustion instability is observed in all cases, and the intensity of combustion instability decreases with a decrease in the equivalence ratio; however a unique phenomenon is observed for the lowest equivalence ratio of 0.6. The unique phenomenon here is the time variation of amplitude of pressure oscillation, whose behavior has been also observed in some other studies. At this condition, while the frequency of pressure oscillation is temporally constant, the frequency of heat release rate varies with time, which is investigated with newly proposed index Time Gap. In addition, the spatial distribution of heat release rate temporally changes with varying frequency. They cause the time variations of correlation between pressure and heat release rate, and finally the amplitude of pressure oscillation temporally varies.


Introduction
Global warming and air pollution are serious global problems. However, owing to internationalization and development, the demand for air transportation is increasing. As resolved in ICAO (2010), the regulation, especially for NOx, is becoming stricter; therefore, the development of jet engines that emit less NOx is necessary. There are several approaches for reducing NOx emissions. An approach, that is attracting social attention is Lean Premixed Prevaporized combustion method (JAXA, 2007). By using this method, the combustion temperature is decreased with a decrease in fuel amount, and then NOx emission is reduced. However, this method can result in unstable combustion and combustion instability can occur (Lieuwen and McManus, 2003). Combustion instability is one of the phenomena of unstable combustion, and it causes big pressure oscillation in the combustion chamber, severely damaging the chamber or result in significant combustion noise. These problems are critical for an airplane; therefore precise prediction of combustion instability is indispensable considering the safety and the efficiency of the jet engines.
Combustion instability is caused by the interaction of pressure oscillation and heat release rate oscillation (Rayleigh, 1878). However, there are many components that affect or induce combustion instability, such as turbulent flow, fuel particle movement, evaporation, and interaction with wall; hence, it is not well understood. Therefore, to elucidate Nagao, Pillai and Kurose, Journal of Thermal Science and Technology, Vol.15, No.3 (2020) the mechanism of combustion instability precisely, many experiments and numerical simulations have been conducted (Gotoda et al, 2011, Kabiraj et al, 2012, Tachibana et al, 2015, Kitano et al, 2016. Experimental studies have suggested some solutions and characteristics of combustion instability; however they did not fully elucidate its mechanism yet because of the limitation of data sampling. In contrast, numerical simulations have investigated abundant data in detail owing to the development of the scheme to calculate the combustion instability precisely (Sato et al,2009, Wolf et al, 2012, Franzelli et al, 2012, Tachibana et al, 2015, Kitano et al, 2016, although spray combustion is complex due to a lot of components that affect the combustion characteristics. In spray combustion, the behavior of injected fuel is affected by the pressure oscillation, and then atomization process is also influenced. Moreover, the regions where fuel droplets evaporate are disturbed by complex turbulent flow under combustion instability. Further, in proceeding numerical studies, Kitano et al. (2016) and Pillai et al. (2020) investigated the influence of droplet size for spray combustion instability; however, in their study, the equivalence ratio was set as 1.0 for simplification, and the effect of equivalence ratio on combustion instability was not investigated.
Therefore, in this study, the equivalence ratio is varied 0.6, 0.8, and 1.0, and the effect of equivalence ratio on the combustion instability in turbulent spray combustion is investigated in a combustion chamber with back-step using LES. The dynamic thickened flame model is used as the turbulent combustion model, and a two-step global reaction model is used for kerosene-air flames. The influence of pressure oscillation on fuel atomization behavior is considered using a model, which can predict the fuel droplets size with time.

LES of spray combustion
The governing equations for gas phase are the conservation equations of mass, momentum, energy, mass fraction of chemical species and equation of state for ideal gas (Kitano et al, 2016), and these equations are solved with semi-implicit compressible solver (Moureau et al, 2007) as below.
Here, the overbar¯denotes the filtered mean value of the physical quantity in gird scale for LES, and the tilde denoets the Favre-averaged value. Each value ρ, P, h, Y k , R, T is density, pressure, enthalpy, mass fraction of species k and gas constant, temperature respectively. u is velocity, and t is stress tensor considering sub grid scale (SGS) stress. Each D h , D k and D t is thermal diffusivity, mass diffusion coefficient of species k which is given under the unity Lewis number assumption and turbulence diffusion coefficient respectively. The value D h and D k are expressed as ρD h = λ/c p and ρD k = λ/c p , and here, λ is heat conductivity and c p is specific heat capacity. These SGS terms are calculated using the Dynamic Smagorinsky model (Moin et al, 1991, Pierce andMoin, 1998). The effect of combustion product is taken into account as S comb,k . E and F are the efficiency function and flame thickening factor and Ω is flame sensor from Dynamically thickened flame model (Legier et al, 2000, Charlette et al, 2002, Strakey and Eggenspieler, 2010, Proch and Kempf, 2014, Rittler et al, 2015. This model ensure that there are at least 8 grid points across the flame thickness. The interaction of fuel droplet and gas flow is considered with S ρ , S ρu , S ρh and S ρY k , with Particle-Source-In-Cell (Crowe et al, 1977) method.
The behavior of fuel droplets is captured by a lagrangian specification, as below.
Nagao, Pillai and Kurose, Journal of Thermal Science and Technology, Vol.15, No.3 (2020) Here, T is the gas temperature, c p,d is the specific heat of droplet, L V is the latent heat of evaporation at T d and τ d is the particle response time. Pr is Prandtl number and Nu is Nusselt number. f s is the correction of the Stokes drag and f t is the correction of heat transfer for an evaporating fuel droplet.
As evaporation model, Langmuir-Knudsen model is used (Bellan and Harstad, 1987, Miller et al, 1998, Miller and Bellan, 1999, Kitano et al, 2014. As fuel, kerosene (hereafter referred to as KERO) that is represented by three components of C 10 H 22 (76.7 wt%), C 9 H 12 (13.2 wt%) and C 9 H 18 (10.1 wt%) is used, and for the calculation of reaction, a two-step global reaction model (Franzelli et al, 2010) is used. With this reaction model the important combustion characteristics, such as laminar flame speed and adiabatic flame temperature are captured accurately in a wide range of equivalence ratio.
Moreover, under the occurrence of combustion instability, the inflow gas velocity oscillates in tune with the pressure oscillation, which influences the atomization of fuel spray. In order to consider the effect on fuel spray atomization, an atomization model proposed by Lee et al (2017) that is theoretically derived from the conservation equations of momentum and energy is adopted. With this model, Sauter Mean Diameter (hereafter referred to as SMD) can be calculated based on the condition of combustor and the effect of fluctuations are able to be captured.  Figure 1 shows the computational domain of only the combustor and conditions for spray combustion, and the computational domain for whole region including inflow region and outflow region is almost same as previous study (refer to Fig.1 in previous study (Kitano et al, 2016)). In combustor, it has 230 grid points in x direction, 80 grid points in y direction and 80 grid points in z direction. Moreover, in whole computational domain including inflow and outflow region, 1120 grid points in x direction, 160 grid points in y direction and 120 grid points in z direction. Air is injected from the inlet, and fuel (KERO) is injected into the combustor from positions situated 5 mm upstream of the back-step, with a velocity of 2 m/s. The injected droplet size distribution is calculated with a modified Nukiyama Tanazawa function (González-Tello et al) that requires SMD calculated by the atomization model mentioned above, as a parameter. The initial air temperature is set as 760 K, the initial pressure is set to the atmospheric pressure of 0.1 MPa, and the initial temperature of injected fuel droplets is 300 K. The equivalence ratio of the 3 cases investigated in this study, varies from 1.0 to 0.6 (i.e., Φ=1.0, 0.8 and 0.6). The equivalence ratio for each cases is defined based on the fuel injection rate and the initial air mass flow rate (which is same among all 3 cases). The value of pressure is investigated 3 mm just above the step at the center plane, and the value of heat release rate is represented by the integral of whole domain in the combustor. The LESs are preformed using an in-house code FK 3 (Kitano et al, 2016, Kitano et al, 2014 with the finite difference formulation in Cartesian coordinate system. There are approximately 2.4 million grid points within the combustor, while the total number of grid points in the entire computational domain is about 21.5 million. The minimum grid size is about 0.1 mm and the reliability of the grid resolution for pressure and heat release oscillations was validated by performing LES with double coarser grid resolution. The physical time for each case is about 730 hours for Φ=1.0 and Φ=0.8 and 1460 hours for Φ=0.6 by parallel computation using 1024 cores on a CRAY-XC40 supercomputer at the Academic Centre for Computing and Media Studies (ACCMS), Kyoto University. This computational time differences derive from the unsteady oscillation at Φ=0.6. For equivalence ratio=0.8 and 1.0, after the amplitude of pressure grows, the amplitude is almost constant, and then the statistical data is collected. However for equivalence ratio=0.6, after the amplitude grows, the amplitude gets smaller. In order to confirm if this trend is going to be observed repeatedly or not, the simulation is continued to see at least two of the same trend. This is why only for equivalence ratio=0.6 take so long time.

Pressure oscillation under combustion instability
In this section, the effect of equivalence ratio on the behavior of combustion instability is investigated. Figure 2 shows the external view of combustor under combustion instability at Φ=1.0. The flame is represented by the red surface which is the isothermal surface of 1800 K and fuel droplets are represented by green dots. It is observed that the droplets flow downstream with air inflow and get burned.
For all equivalence ratios, combustion instability is observed; however, the behavior of pressure oscillation varies among them. The time variations of pressure, P, at each equivalence ratio are shown in Fig.3. The frequencies of pressure oscillations are 699Hz for Φ=1.0 and Φ=0.8, and 677Hz for Φ=0.6 "Time Gap" in this figure is discussed below. Here, only for Φ=0.6, time range is longer than the others, and this is because to capture the peculiar phenomena that is observed only at Φ=0.6. In all cases, the pressure oscillates periodically, and the amplitude of pressure oscillation is in the range of approximately 2-10%. At Φ=1.0, the pressure oscillates stably, and the amplitude is almost constant. At Φ=0.8, the amplitude varies slightly, but it is almost constant like Φ=1.0. However, at Φ=0.6, the amplitude is not constant and exhibits repeated small and large variations. This unique phenomenon is not observed for other equivalence ratios, and this fluctuating amplitude was also observed in some experimental studies (Gotoda et al, 2011, Kabiraj et al, 2012. However, the detailed mechanism for this phenomenon is not still fully elucidated.

Time Gap
As mentioned above, at Φ=0.6, the amplitude of pressure oscillation fluctuates with the lower frequency compared with the frequency of pressure oscillation. It is known that the characteristics of combustion instability are mainly determined by pressure, heat release rate, and their correlation. Therefore, as index to determine the degree of correlation, Rayleigh Index (RI) is often adopted. This index shows the spatial correlation between pressure and heat release rate, but it is difficult to investigate the time variation of correlation. At Φ=0.6, the amplitude of pressure oscillation varies with time, and there should be the time variation of correlation between pressure and heat release rate. To understand the time variation of the correlation, the quantity "Time Gap" is proposed as an index, and its time variation is analyzed. The Time Gap is calculated as follows and Fig.4 shows the conceptual diagram of Time Gap calculation. First, cross correlation for a certain time range is calculated as Here, P ′ and q ′ are variations of pressure and heat release rate, respectively, and t p is the one cycle time of pressure oscillation. τ is the imaginary time delay of pressure from original pressure, and τ varies from −t p /2 to t p /2 to evaluate at least one cycle time of pressure oscillation. The value of cross correlation represents to what extent the oscillation of imaginary time delayed pressure and one of heat release rate overlap, and this value is changed with varying τ because the phase difference between pressure and heat release rate changes. After the calculation of the cross correlation, among the τ that provides the maximum value of Cross(τ), the one closest to zero is defined as the Time Gap. Here, maximum Cross(τ) means that the oscillation of pressure and heat release rate most closely overlap, and τ that gives maximum Cross(τ) represents how close the phase of original pressure and heat release rate are. Therefore, the smaller τ means the closer phase difference between pressure and heat release rate. The reason to choose the τ closest to zero is to consider only one cycle time. For example, when pressure has maximum value at the phase when heat release rate has minimum value, Cross(τ) could have two maximum value with τ being about ±t p /2. In that case, the closer τ to zero, which means the phase of pressure is closer to one of heat release rate, is adopted as Time Gap. (c) Φ=0.6  Nagao, Pillai and Kurose, Journal of Thermal Science and Technology, Vol.15, No.3 (2020) By using this technique, the Time Gap are calculated, and Fig.3 shows the time variations of the Time Gaps and pressure oscillation for each equivalence ratio. At Φ=1.0 and Φ=0.8, Time Gap does not vary and fall between 0.0 and -0.2ms (which correspond to ϕ=0 • and ϕ=50 • and in this study the phase ϕ=0 • corresponds to the time instance when pressure has minimum value). However, focusing on Φ=0.6, Time Gap varies with time, and when the amplitude of pressure gets the maximum, Time Gap lies between 0.0 and -0.2ms, and when the amplitude gets minimum, Time Gap has a wide time range and at maximum it reaches 0.8ms which corresponds to approximately ϕ=180 • . Along with the time variation of Time Gap, the amplitude of pressure oscillation fluctuates. Moreover, Time Gap lies between 0.0 and -0.2ms before amplitude gets maximum value. This implies that the pressure oscillation is increased with the concentrated correlation between pressure and heat release rate. Further this temporal proceeding of Time Gap because of the fluctuating amplitude of pressure could be used to predict the occurrence of combustion instability. By using the measured data of pressure and heat release rate, the Time Gap is calculated, and the time when the Time Gap is in the narrower range of time than before, it could be the premonitory symptom of increasing combustion instability.
Moreover, the Time Gap is approximately between 0.0 and -0.2ms when it is in the narrower range of time, which implies that when the phase of pressure slightly delays the phase of heat release rate, the correlation gets stronger, and it is not the time when the phase of pressure completely matches the phase of heat release rate, which means the Time Gap is 0.0ms. This time variation of the Time Gap is strongly affected by the time variation of pressure and heat release rate. Therefore, these factors are investigated below.

Pressure and heat release rate
The amplitude of pressure oscillation fluctuates for Φ=0.6, and the characteristics of combustion could be different for each amplitude; therefore, based on the amplitude of pressure oscillation, the series of oscillation phenomenon is partitioned into four time intervals; a, b, c, d, and the characteristics of combustion instability are investigated in detail in each time interval. Here, each time interval is approximately 0.01s long as shown in Fig.5. First, the effect of pressure is investigated. Considering the time for one cycle of pressure oscillation, each cycle time is almost the same with a small variation less than 0.1ms. In contrast, the Time Gap varies from -0.6 to -0.1ms, which implies that the range of variation of Time Gap is 0.5ms. The range of variation of pressure is smaller than that of Time Gap; thus, the time variation of Time Gap is more strongly affected by the time variation of a phase of heat release rate than that of pressure. Hereafter, this time variation is expressed as temporal fluctuation. Moreover, Fig.6 shows the streamwise distributions of time-averaged pressure for each time interval.
The position of transition from low pressure to high pressure is the same for all time intervals; only the magnitudes differ. Around the step, there is an anti-node of pressure oscillation, and around the exit of the combustor, there is a node. As mentioned above, each cycle time of pressure does not fluctuate considerable, and as shown in Fig.6, streamwise distributions of pressure in all time intervals are qualitatively similar in shape. The only difference is the magnitude of pressure oscillation and this is the only factor affecting the correlation and hence the Time Gap variation. However, this factor alone is not sufficient for the investigation of time variation of Time Gap; therefore, the heat release rate is also investigated below. Figure 7 shows the streamwise distribution of the time interval's averaged and cross-sectional averaged heat release rate. From this figure, it can be seen that the peaks of heat release rate of time intervals b and d, i.e., the time intervals with larger amplitudes of pressure oscillation, are bigger than those of heat release rate of time intervals a and c, i.e., the time intervals with smaller amplitudes of pressure oscillation. Additionally, the peak heat release rate positions of time intervals b and d lie further upstream than those of time intervals a and c, and their second peak lies in the downstream region. Considering the pressure, the upstream region of the combustor is the antinode of pressure oscillation; therefore,  when the distributions of heat release rate shift upstream, regions of high pressure and high heat release rate overlap and enhance the correlation. Therefore, the combustion instability becomes stronger.
To clarify the reason as to why the heat release rate shifts upstream during time intervals b and d, the streamwise distributions of time averaged and cross-sectional averaged evaporation rate, and the x-y plane instantaneous distributions of heat release rate at the phase ϕ=180 • , which is defined as the time instance when pressure has maximum value, are shown in Fig.8 and Fig.9, respectively. Figure 8 shows that the distributions of the evaporation rate do not change appreciably among the different time intervals, which implies that the distributions of the evaporation rate does not affect the transition of combustion instability. The reason why the distributions of evaporation do not change can be explained by the droplet size distribution and the Nagao, Pillai and Kurose, Journal of Thermal Science and Technology, Vol.15, No.3 (2020)  position where droplets get evaporated. Figure 10 shows the time averaged Probability Density Function (PDF) at each time interval. As this figure shows, the droplet size distributions do not change among each time interval. Moreover most droplets evaporate at x<0.02m, and the shapes of vortex behind step change drastically x>0.02m among each time interval as shown if Fig.9. Moreover, the average inflow velocity is same for all time interval, therefore, there is not any factor that changes droplet evaporation. On the other hand, the x-y plane distributions of the heat release rate varies among the time intervals. During time intervals b and d, the spatial distribution of the heat release rate flaps up and down, and is concentrated around x=20-30mm. To investigate this, x-y plane instantaneous distributions of fuel mass fraction and heat release rate at some phases in a single period of pressure oscillation are shown in Fig.11. At the phase ϕ=0 • , when the pressure has minimum value, fuel mass fraction stretches along the downstream direction, and as a result, the heat release rate is distributed in a similar manner to fuel mass fraction. Subsequently, with increasing pressure, the distribution of fuel mass fraction separates into two regions upstream and downstream. Then, at ϕ=180 • , when pressure has maximum value, the upstream fuel is pushed back further upstream, and, the downstream fuel diffuses further downstream. Consequently, because the distribution of the fuel mass fraction is shorter than that at ϕ=0 • , the heat release rate has high value locally in the upstream region of the chamber. Above ϕ=180 • , with a decrease in pressure, fuel that was congested upstream, expands downstream with swirling and is burned; therefore, considerable heat is released in the entire combustion chamber. This phenomenon; congestion, diffusion, and combustion that is recognized as heat fluctuation, are hardly observed at time intervals a and c. Further, for each phase the distributions of fuel do not change considerably from 0mm to 20mm. Therefore, the stream-wise distributions of evaporation rate do not change among different time interval. This spatial fluctuation of heat release rate distribution is triggered by the non-periodical turbulence of flow and flame.
The Time Gap has large values in the time intervals a and c, and small values in the time intervals b and d. This difference is caused by the heat release rate's temporally and spatially fluctuation. In time intervals a and c, the distribution of heat release rate does not vary, which implies that when the pressure gets its minimum, the heat release rate cannot decrease and keep in high value. Consequently, the correlation of pressure and heat release rate gets weaker. In contrast, in the time intervals b and d, the distribution of heat release rate fluctuates, and its maximum value also fluctuate stronger than in the time intervals a and c. Therefore, the correlation gets stronger, and combustion instability increases.
As mentioned above, the heat release rate distribution temporally fluctuates and degree of correlation also fluctuates; therefore, the spatial correlation should be different. Hence, to investigate the spatial correlation among time intervals, RI is shown in Fig.12. RI is calculated as, Here, P is pressure, Q is heat release rate, t s is the sampling time, dash, ', represents the variation, and ave represents Nagao, Pillai and Kurose, Journal of Thermal Science and Technology, Vol.15, No.3 (2020) Fig. 11 Comparison of instantaneous distributions fuel mass fraction and heat release rate on the x-y plane among different phases of time interval d.  Nagao, Pillai and Kurose, Journal of Thermal Science and Technology, Vol.15, No.3 (2020) the averaged value. The region where RI is positive represents the region where combustion instability is encouraged, and vice versa. As the figure shows, the positive area is larger in the time intervals b and d than the time intervals a and c, and the maximum RI value is also higher. This difference of RI distribution is one of the reasons of the difference in the amplitude of pressure oscillation. At approximately x=30mm, strong positive regions are distributed. This is because of the spatial fluctuation of the heat release rate as shown in Fig.11. When the heat release rate stagnates upstream, RI has a high value around x=30mm, because this region is the anti-node of pressure oscillation and in this phase, the pressure gets its maximum.

Conclusions
In this study, combustion instability of spray combustion in a combustion chamber with a back-step was investigated using LES employing a dynamic thickened flame model. In this setup, three equivalence ratios (Φ=1.0, 0.8, and 0.6) were simulated. The combustion instability was observed in all of cases, and the intensity of combustion instability decreased with a decrease in the equivalence ratio. This is mainly because whole heat energy decreased with the decreasing amount of fuel injected into combustor and the unstable combustion phenomena such as blow off did not occur in this setup. However, for the lowest equivalence ratio of Φ=0.6, the amplitude of the pressure oscillation fluctuated temporally, and its tendency was quite different from that observed under combustion instability keeping a certain amplitude of pressure oscillation as for Φ=0.8 and 1.0. In order to investigate the time variations of correlation, a new index "Time Gap" was proposed and confirmed useful. With using Time Gap, it was elucidated that the frequency of heat release rate varied with time, and the phase difference between pressure and heat release rate got longer or shorter. In addition, it was also observed that, when the phase difference got longer, the flame structure was stable without fluctuating, and when the phase difference got shorter, the flame structure began to fluctuate and made the spatial distribution of heat release rate vary with time. Because of these temporal and spatial fluctuations of heat release rate which induced the time variation of the correlation between pressure and heat release rate, the amplitude of pressure oscillation temporally fluctuated.