Nonlinear analysis of low-frequency combustion instabilities in liquid rocket engines

Low-frequency combustion instabilities are here studied taking advantage of the software EcosimPro. A specific module has been implemented based on the double time lag model and the coupling of combustion chamber and feed line oscillations were investigated by using a complete set of nonlinear equations. The characteristic time lags have been identified following two approaches: (i) a constant time lag approach; and (ii) a variable time lag approach based on correlations available in open literature. To prove the module capabilities, an experimental setup was reproduced and a stability map was generated, comparing the obtained results with literature data from both experiments and a linear double time lag model. The stability boundaries obtained with the chugging module are in good agreement with those obtained in open literature and the first characteristic frequency of the engine is well predicted. Furthermore, the model proves its capability in reconstructing the reversal in the slope of the stability boundary at low fuel injector pressure drops and in detecting the high-frequency content typically observed in presence of multimode oscillations. However, in the calculations, the higher frequency does not dominate the instabilities, that is, in the unstable regime, the model diverges with a frequency equal to the first characteristic frequency. In the last part of the paper, the variable time lag approach is used to investigate a portion of the aforementioned stability map. Thanks to the semiempirical correlations, the present authors managed to improve the prediction of the first characteristic frequency, whereas the stability boundary does not change significantly and remains comparable with the one predicted by the constant double time lag approach.

Low-frequency combustion instabilities are here studied taking advantage of the software EcosimPro. A speci¦c module has been implemented based on the double time lag model and the coupling of combustion chamber and feed line oscillations were investigated by using a complete set of nonlinear equations. The characteristic time lags have been iden-ti¦ed following two approaches: (i) a constant time lag approach; and (ii) a variable time lag approach based on correlations available in open literature. To prove the module capabilities, an experimental setup was reproduced and a stability map was generated, comparing the obtained results with literature data from both experiments and a linear double time lag model. The stability boundaries obtained with the chugging module are in good agreement with those obtained in open literature and the ¦rst characteristic frequency of the engine is well predicted. Furthermore, the model proves its capability in reconstructing the reversal in the slope of the stability boundary at low fuel injector pressure drops and in detecting the high-frequency content typically observed in presence of multimode oscillations. However, in the calculations, the higher frequency does not dominate the instabilities, that is, in the unstable regime, the model diverges with a frequency equal to the ¦rst characteristic frequency. In the last part of the paper, the variable time lag approach is used to investigate a portion of the aforementioned stability map. Thanks to the semiempirical correlations, the present authors managed to improve the prediction of the ¦rst characteristic frequency, whereas the stability boundary does not change signi¦cantly and remains comparable with the one predicted by the constant double time lag approach.

INTRODUCTION
In this paper, a numerical analysis of low-frequency combustion instabilities, frequently referred to as chugging, in Liquid Rocket Engine (LRE) is proposed taking advantage of an unsteady nonlinear approach. Low-frequency combustion instabilities in LREs are characterized by oscillations in the mass §ow rate of propellants that result in pulsations of the chamber pressure. The response of the combustion chamber to the upstream oscillations is not immediate. The delay is due to the di¨erent characteristic times associated with the processes of injection, atomization, vaporization, mixing, and combustion, taking place in quick succession inside the combustion chamber. Moreover, the mass §ow rate is not immediately responsive to pressure oscillations because of the propellant inertia in the feed lines. Hence, both feed lines and combustion chamber physics must be taken into account. Low-frequency instabilities have been usually investigated with the help of linearized models [15], whereas only a few attempts have been made to detail the phenomena with nonlinear models [6,7]. In the present study, an unsteady nonlinear analysis of chugging instabilities is proposed with particular attention to the appropriate modeling of delays induced by atomization, vaporization, and mixing processes. The analysis of chug instabilities requires to model the subsystems of an LRE to investigate the mutual in §uence of each component. Such a comprehensive simulation inherently needs simpli¦ed submodels able to give a satisfactory level of reliability along with a reasonable computational time. In this context, an object oriented simulation tool is a valid option to handle the complexity of a propulsion system. EcosimPro [8] is a simulation platform focused on the analysis of complete systems that spans di¨erent engineering ¦elds. Thanks to the object oriented philosophy, single components can be connected to each other in order to analyze systems at di¨erent complexity levels. Propulsion systems can be investigated by the European Space Propulsion System Simulation (ESPSS) library [9,10], in which the typical components of an LRE (tanks, turbomachinery, feed lines, valves, gas generator/preburner, combustion chamber, etc.) are modeled to study both steady-state and transient phases.
To detect low-frequency instabilities, it is important to properly model the characteristic times associated with both liquid and gases injected into the combustion chamber. Chug instabilities arise when the total time lag is such that oscillations in the rate of propellant injection are coupled with the corresponding oscillations of the chamber pressure. Di¨erent studies have been performed during the past decades to investigate which time lag is the most important for the generation of low-frequency instabilities. Most of the proposed chug models focused on steady linear analyses are based on the linearization of the continuity equation in the combustion chamber. In the attempt of explaining damages caused by chugging in LRE, Gunder and Friant [1] introduced the time lag concept. They linearized the energy conservation equation and solved the relative delay di¨erential equation. The time lag approach was further extended by Summer¦eld [2] who retrieved a stability analysis from the linearization of the continuity equation introducing the in §uence of the feed lines. He investigated the e¨ects of feed lines, injector pressure drops, and combustion chamber length on chug instabilities.
A comprehensive theoretical approach based on the concept of sensitive and insensitive time lag was then proposed by Crocco and Cheng [3,11,12]. Their approach is based on the concentrated combustion model and applies to both monopropellant and bipropellant cases. Later, Wenzel and Szuch [4,1315] developed a model that incorporates a double time lag approach, hence extending the potentialities of the original Crocco£s model to the cases in which the two propellants have di¨erent time lags. They applied the model to analyze a system ignoring the in §uence of the feed lines and considering only the injector impedances. They used a useful engineering approach to express the solutions, highlighting the importance of injector pressure drops on instabilities. Recently, Casiano [5] extended the Wenzel and Szuch£s model by formulating the problem in a parametric form and adding the in §uence of the feed lines. Linear models have been also developed to study Pogo instabilities that are combustion instabilities characterized by the coupling of combustion and structural dynamics and that usually bring to structural as well as to thrust oscillations. Pogo instabilities arise at low frequencies well below 100 Hz, thus making Pogo models appealing for the study of chug instabilities. Such is the case of Ordonneau et al. work, where a Pogo model, originally meant to study Pogo instabilities in the VULCAIN engine, has been adapted to both analyze chug instabilities of VULCAIN and predict those of VULCAIN 2 engine [16,17]. Beside these linear approaches, a few attempts have been made in the ¦eld of nonlinear analysis, mainly because the computational cost of such approaches was too high for the hardware resources. Webber£s [6] work was the ¦rst attempt of performing an unsteady analysis based on nonlinear equation for both the feed lines and the combustion chamber. He focused his attention towards the modeling of the liquid phase in the combustion chamber, showing the e¨ects of droplet sizes on the amplitude of the low-frequency oscillations. Bartrand [7] proposed an extension of this model, including the e¨ects of ¦nite rate chemistry and showed that in the Space Shuttle Main Engine preburner case, the assumption of ¦nite combustion velocity a¨ects the magnitude of the damping rather that the frequency of oscillations.
The analyses conducted in this paper take advantage of the ESPSS feed lines modeling approach, whereas the combustion chamber component is suitably improved by means of a chug dedicated module. It is well known that the shift, or time lag, between disturbances coming from injectors and pressure response is the sum of ¦ve time lags [18], each related to a speci¦c process taking place in the combustion chamber. The ¦ve processes are, in order: injector response, atomization, vaporization, mixing, and combustion. It is commonly accepted that the most important delays are those related with atomization, vaporization, and mixing [7,19], whereas spatial droplet distribution produced by injection has been shown to have a great e¨ect on high-frequency instabilities but has a little or no e¨ect on chugging [15]. Moreover, in the case of liquid oxygen and gaseous hydrogen propellant combination, combustion can be considered faster than other processes and the related time lag can be ignored. On the contrary, the vaporization rate is of great importance for the evaluation of chugging and is strictly related with the atomization process (droplet size and surface area). Therefore, the present authors focus their attention on atomization, vaporization, and mixing phenomena. Two di¨erent approaches are foreseen: a constant time lag approach and a variable time lag approach. In the constant approach, both atomization and mixing time lag remains unchanged throughout the entire simulation. Instead, for the variable approach both atomization [20] and vaporization [20,21] will be modeled with the aim of semiempirical correlations while mixing times will be retrieved from an empirical curve [14]. In this way, the characteristic times depend on the main combustion chamber parameters, thus changing during the simulation. The chug module with the two approaches will be validated against experimental and numerical data available in open literature [22] and the di¨erences between the two methods will be discussed.

NUMERICAL METHOD AND APPROACH
The LRE system is here modeled with EcosimPro, taking advantage of the components included in the ESPSS library, for which a complete description can be found in di¨erent works [9,10,23]. In the following, the feed lines and injector components are brie §y introduced, whereas more details are provided for the combustion chamber component. The modeling approach for the gas phase in the combustion chamber is particularly important for the scope of the present paper, because it is the frame in which the vaporization, atomization, and mixing models are included.

Feed line and injector components
The tube component is modeled with a quasi-one-dimensional, unsteady approach for two phase §ow. The two phase §ow evolution is described by the Homogeneous Equilibrium Model (HEM) [24], which is based on two assumptions: (1) thermodynamic equilibrium between the two phases; and (2) same velocity for both liquid and vapor.
Thanks to these hypothesis of §ow in equilibrium, it brings to a set of equations whose structure is close to the Euler equations and, more importantly, that are unconditionally hyperbolic. The two-phase §uid can be either a simpli¦ed liquid or a real §uid. The governing equations consider one continuity equation for the noncondensable perfect gas, one for the noncondensable perfect gas diluted in the liquid phase, one for the overall mixture, one momentum equation, and one energy equation for the overall mixture. The set is completed by N equations, one for each of the chemical species convected in the pipe. These latter continuity equations are needed when the pipe is located downstream of a component in which reactions take place, e. g., preburner and gas generator, and the convection of the combustion gases must be taken into account. It is worth noticing that both the noncondensable gas and the two-phase §uid can be two di¨erent real §uids, whose properties are retrieved from a dedicated database. Di¨erent §ow mixture con¦gurations are possible, e. g., pure liquid, liquid and vapor, and vapor and noncondensable. Therefore, §ow properties must be calculated according to the mixture state and proper mixing rules have to be adopted. Injectors and injector cavities are located between the feed lines and the combustion chamber. Injector cavities are modeled by means of a two-phase set of equations for a nonadiabatic variable volume. Two mass conservation equations and one energy equation describe the evolution of the two-phase mixture. The velocity in the cavity is calculated as an average velocity retrieved from the mass §ow rate of the incoming §uid (calculated by the tube component) and the §uid going to the combustion chamber (calculated by the injector component); hence, no momentum equation is solved. Finally, the injector component is based on a momentum equation and calculates the mass §ow across the injector in both turbulent and laminar regimes. The momentum equation accounts for concentrated pressure losses and for the inertia associated with the connected components. Once the mass §ow rate is known, the energy §ux is simply calculated as the product of the injector mass §ow rate and the stagnation enthalpy of the §uid in the upstream cavities. In this way, both mass and energy §uxes are provided to the combustion chamber component.

Combustion chamber
The thrust chamber component consists of two subcomponents. The ¦rst one is an unsteady component that models the combustion chamber up to the throat. In this part of the combustion chamber, the §ow is considered either in chemical equilibrium or in a ¤time delayed¥ equilibrium state that is a simple approach mimicking ¦nite rate models. In the present analyses, we will consider the component as based on the chemical equilibrium assumption. Instead, the divergent part is a quasi-steady component that mimics the supersonic expansion by means of analytical correlations for a §ow in either frozen or equilibrium conditions. The ¦rst part of the combustion chamber is based on a quasi-one-dimensional formulation of the governing equations. The component is nonadiabatic, hence allowing for thermal connections with the external components, i. e., cooling channels. The unsteady conservation equations are: where ρ must be intended as the density of the mixture of gases and the subscript ¢bu£ refers to the burned gases that are calculated only when the ¤time delayed¥ rate approach is used. No further details are here given on these terms since we do not use such a modeling approach. The ¦rst source term on the right-hand side of the momentum equation models the friction e¨ects, while the second is related to the e¨ects of a variable area. The ¦rst source term in the energy equation ' q w is the heat exchanged with the combustion chamber walls while the second term ' q unst is an excitation term introduced in the present analysis to mimic the roughness of the combustion processes. Finally, the vaporized mass source term in the continuity equation ' m vap is expressed according to the vaporization models described in the next section.
The nozzle component is a quasi-steady component able to describe the evolution of a §ow in a duct with variable cross section. The nozzle takes the upstream thermodynamic state from the last computational cell of the combustion chamber and characterises the throat and the diverging section conditions. It should be noted that the modeling of the supersonic evolution is not directly needed to properly model chug instabilities since the only boundary condition that a¨ects the detection of chug instabilities is the one giving the mass §ow evolution at the throat section. However, the component is left in its original formulation as the presence of the diverging nozzle does not limit the chug model capabilities.

Time lag approach
Once injected into the combustion chamber, the liquid subsequently undergoes atomization, vaporization, mixing, and combustion. The time lags of each process must, hence, be characterized in an appropriate way in order to catch the coupling between the feed lines and the combustion chamber. The general numerical procedure is the following: once injected in the combustion chamber, the mass associated with the liquid droplet is stored in a memory array and conserved for a time equal to the total time lag (vaporization, atomization, and mixing) associated with the droplet. After the ascribed time lag, the incoming mass is injected as a vaporized source term in the continuity equations (1). Two di¨erent cases are foreseen. In the ¦rst one, the time lags are user-de¦ned input parameters and remain constant throughout the simulation. The second case is, instead, based on a time delayed approach in which the vaporization and atomization processes are described by time lags based on chamber conditions and liquid phase properties. The mixing times are retrieved from an empirical curve [14], while the combustion is considered as su©ciently fast and its time lag is neglected. In the present study, we decided to focus our attention towards the widely used oxygenhydrogen combination, hence the correlations applied to the characteristic times hold only for this propellants. In particular, the hydrogen is injected in a gaseous phase; therefore, the delay characterizing this propellant is only the mixing time lag. The oxidizer is instead injected in a liquid phase; hence, it undergoes also vaporization and atomization. The time lags and the related parameters are calculated as follows. For the atomization processes, we apply the correlation [20]: We 0.03 D l is the liquid post inner diameter; u l is the liquid velocity at injection; σ l is the liquid surface tension; Re l is the Reynolds number for the injected droplet; and We g is the Weber£s number for the combustion chamber gases. Vaporization characteristic times for droplets when in presence of convection [20,2527] are usually retrieved as a correction of those obtained for droplets immersed in a quiescent gaseous environment [28]. We here make use of the correlation [20]: We 0.03 where p c is the chamber pressure; MR is the mixture ratio; T ∞ is the temperature of the combustion chamber gases; T cr is the liquid critical temperature; and D 0 is the initial droplet diameter. The initial droplet diameter D 0 is here calculated from a correlation speci¦cally developed for coaxial injectors and liquid oxygen [21]: with V r being the relative velocity between injected liquid and gases in the combustion chamber, D g the annular gap outer diameter, R post the recess length, C inj the parameter taking into account the injector design ranging from 0.4 to 1.2, and K the correction factor to be applied when the experimental simulants di¨er from the actual propellants [29].

TEST CASE
In order to prove the chugging module capabilities, let us refer to an experimental test case speci¦cally developed to validate a linear double time lag model [22]. The experimental apparatus is a gaseous-hydrogen liquid-oxygen engine. The injector cavities and injector plate were speci¦cally designed to ensure the desired pressure drops and a constant pressure in the injector cavities. The main engine parameters are described in Table 1, while the EcosimPro schematic and its component descriptions are proposed, respectively, in Fig. 1 and in Table 2.
Following the approach proposed in [22], let us de¦ne a stability map for the engine by performing simulations with variable injector pressure drops. Changing both fuel and oxidizer pressure drops, it is possible to create a map in which stable and unstable regions can be identi¦ed. The combustion inside the chamber is a rough process; hence, spontaneous oscillations in the pressure values are   Thermal insulation (zero heat §ux) SZC CHUG Combustion chamber usually experienced. In order to mimic these oscillations and excite the chamber to eventually initiate the closed unstable loop between the camber and the upstream §ow, we decided to introduce an unsteady term in the energy equation. The purpose here is to reproduce the combustion roughness, phenomena commonly identi¦ed as the exciting mechanism for combustion instabilities, by means of a broadband signal. The broadband reported in: f 0 , -f , N , and ǫ being, respectively, the ¦rst excited frequency, the sampling frequency, the number of samples, and the magnitude of the perturbation, precisely selects the frequencies of interest and, in the hypothesis of small perturbations, enables the analysis of the chamber response at the selected frequencies.
In this way, frequencies between f 0 and f 0 + N -f are excited every -f . It is worth noticing that the broadband has a speci¦c frequency content and periodicity characterizing its evolution in time. The steady value ' q 0 is selected to mimic the losses measured in the combustion facility in order to match the experimental combustion e©ciency η c * = 0.75. The ¦nal value is kept constant throughout all the simulations and equal to ' q o = 1.933 MW. This gives a low value of the combustion chamber temperature T c = 2038 K as expected from the low combustion e©ciency.
Pressure boundary conditions on both feed lines are changed throughout the simulations to ensure constant mass §ow rates at steady state for di¨erent injectors pressure drops. It is worth saying that EcosimPro has its own component to ensure a constant mass §ow rate in the feed line. This component could have been replaced in the schematic in Fig. 1 by both main oxidizer and main fuel valves, thus imposing the desired mass §ow rates. However, it has been veri¦ed that the use of this component adds an arti¦cial damping preventing the system from falling in any instability regime, hence failing in detecting chugging instabilities.

Constant Time Lag
In the case of constant time lag, the characteristic times are imposed as user input data. The values used here are inferred from the reference work [22]. In particular, the atomization time is, in this case, neglected and the vaporization lag is retrieved from the length needed to vaporize the 50 percent of the injected liquid mass l 50 [30]. Once the length is known, the vaporization time is calculated as the ratio of the vaporization length and the injection velocity of the liquid, τ vap = l 50 /v inj . Calculations for the examined engine con¦guration give a value of τ vap = 4.4 ms. The mixing and reaction times are then retrieved from the neutral stability condition. Starting from the experimentally observed chugging frequency at high injector pressure drops of 68 Hz, the phase requirement for stability gives an overall time delay of τ tot = 6.7 ms. Hence, the mixing time lag is obtained as a di¨erence of the previously calculated delays τ mix = 2.23 ms.
The results obtained with the chugging module are here compared to those coming from both experimental evidences and Szuch£s double time lag model. The tests are performed by ¦xing one of the two injector pressure drops, varying the other and measuring the response in terms of chamber pressure. The Root Mean Square (RMS) of the oscillations of the measured signal is then extracted and presented in terms of percentage values with respect to the mean chamber pressure, while the frequency content is analyzed by means of the signal Power Spectral Density (PSD). In this way, one is able to retrieve both amplitude of the response and characteristic frequencies of the system. It is worth noticing that the amplitude of the chamber response is strictly related to the amplitude of the input perturbation that, in turn, a¨ects the de¦nition of the stability region. In the experimental campaign, here referenced as [22], a 10 percent RMS criterion was used to identify stable and unstable operating regimes. In particular, all the con¦gurations that give a pressure signal with a RMS exceeding the 10 percent of the mean chamber pressure value are considered unstable, while all the others can be considered stable.
The same consideration applies to the results coming from the simulation performed by Szuch et al. with the double time lag model. Furthermore, Szuch et al. imposed a white noise signal on the mean combustion products values to mimic the roughness of the combustion processes. The white noise input perturbation a¨ects the entire spectrum but we are just interested in the lower part of the spectrum, since chug instabilities are associated with low frequencies. This consideration justi¦es the choice of a broadband input. In the subsequent analyses, the broadband signal in Eq.   the reduction of pressure drops. This is a well-known behavior since increasing the injector pressure drops has a stabilizing e¨ect on the overall system and decouples the oscillations in the feed lines of liquid oxygen from those in the combustion chamber, thus preventing the system from falling in the unstable regime. The predominant frequency resulting from the PSD in the single time lag case (55 Hz) di¨ers from the one of the double time lag case (60 Hz). The result of this prediction has to be discussed considering that the broadband sampling frequency -f introduces an uncertainty of ±5 Hz.
The di¨erence between the two approaches is related to the fact that in the double time lag, also, the fuel feed line couples dynamically with the oscillations in the combustion chamber, thus changing the phase relation between lines and combustion chamber. However, the small registered shift indicates that the chug frequency is mostly determined by the liquid oxygen, since the greater time lag is associated to this propellant and it dominates the phase relation between oscillations. Furthermore, in the single time lag case, the RMS values are always smaller than those of the double time lag approach and a plateau is observed between -p fu /p c = 0.6 and 0.7. Furthermore, both time lag approaches give RMS values of the chamber pressure that are lower than those registered in the reference test cases. This is related with the amplitude ǫ of the broadband in Eq. (5).
In Fig. 3, the in §uence of the amplitude of the broadband signal on the chamber pressure is shown. Increasing the amplitude of the input perturbation increases the chamber pressure response and the values reported in the reference paper are approached. It must be noted that in the double time lag model, the inherent stability limit is placed before the one reported in the reference work and this is true for all the selected level of amplitudes.
In both single and double time lag approaches, the system response to the broadband can be divided in three zones characterized by a di¨erent behavior in terms of both frequency content and signal shape. In Figs. 4 and 5, the time evolution of chamber pressure and the PSD of the same variable in the case of the double time lag model are, respectively, shown. Similar trends have been noted in the single time lag case. At high oxidizer pressure drops, the chamber pressure signal is proportional to the input signal, with a frequency content characterized by a small peak at 60 Hz. The frequency content is simply the one of the broadband perturbation and the dumping shown at this pressure drops must be considered as characteristic of the input signal instead of as system damping. When the oxidizer injector pressure drop is reduced, the 60-hertz frequency becomes predominant and the pressure signal assumes a characteristic damped behavior, now associated to the system. At -p fu /p c = 0.35, the system is still able to damp the perturbation and to prevent the oscillations from growing.  No limit cycle is observed and the chamber pressure oscillations are damped and tend to disappear if not sustained with the input broadband signal. Below the inherent stability limit, -p ox /p c < 0.35, the system is not able to damp the oscillations, and the pressure response grows to high values causing the simulation to fail as shown in Fig. 4d. In this case, the system oscillates at the system characteristic frequency of 60 Hz and the oscillations are self-sustained even when the input is deactivated. Furthermore, the system oscillates at the characteristic frequency regardless the frequency content of the input, that is: the system diverges with a 60-hertz frequency even if a monochromatic input of 20 Hz is imposed for a short time and then deactivated.
In Fig. 6, the stability map obtained for the selected test case with the double time lag model is shown. The results are compared with the boundaries identi¦ed by Szuch et al. It is worth noticing that the 10 percent RMS criterion is a practical criterion [22] to de¦ne a stability boundary. Since this criterion is directly related with the amplitude of the input perturbation, we choose to use a value of ǫ = 0.01 throughout the simulations and we make use of also a 50 percent RMS criterion to identify the stability boundary. The boundaries identi¦ed with both the 10 and the 50 percent RMS criterion are shown with white dashed lines, whereas the reference results are shown in black. When the boundary is iden-ti¦ed with the 50 percent criterion, the reversal in slope at lower -p fu /p c is well matched, and, di¨erently from the analog referenced computer program, no unstable regimes are identi¦ed at high oxidizer injector pressure drops when the fuel injector pressure drop is Figure 6 Stability map. Chamber pressure RMS dependency on oxidizer and fuel injector pressure drops: 1 ¡ 10% RMS; 2 ¡ 50% RMS; 3 ¡ analog Szuch; and 4 ¡ experimental lower than 0.15. The authors of the reference paper [22] jus-ti¦ed this discrepancy assuming a back §ow in the fuel injector element at lower pressure drops; however, this back §ow has never been registered in our simulations. Furthermore, the reversal in slope is due to the intersection of two boundaries: the low-frequency boundary and the high-frequency boundary.
At lower fuel injector pressure drops, characteristic system frequencies ranging from 147 to 210 Hz were observed both experimentally and numerically [22]. In order to investigate the presence of this second frequency, we extend the broadband input up to 200 Hz, measuring the resulting chamber pressure signal. With this input, it is possible to register a higher frequency content of 175 Hz that is well within the limits of the aforementioned highfrequency range, giving a ratio of 2.9 between higher and lower frequency that well approximates the 2.7 ratio noted during multimode oscillations [31]. However, when the system enters its unstable regime at fuel injector pressure drops lower than 0.15, the pressure diverges with a frequency of 60 Hz, thus meaning that in our analysis, the 60-hertz frequency is the characteristic frequency of the system even at lower fuel injector pressure drops. Furthermore, PSD shows that the power associated to the 60-hertz component is always greater than the one at 175 Hz, explaining the predominance of the ¦rst frequency throughout the entire stability map.
For the constant time lag model, we ¦nally investigate the in §uence of different input amplitudes for various fuel pressure drops. The scope is here to investigate how the stability boundaries are in §uenced by the level of the input perturbations. In particular, the results obtained when the amplitude ǫ in Eq. (5) is increased of 20% and 100% are presented. The results in terms of the RMS of chamber pressure oscillations are then normalized with those obtained with the reference case ǫ = 0.01 and divided by the respective percentage increment. In this way, we are able to understand whether the response of the system is linear. In fact, in the case of a lin- ear behavior, the plots should be equal to 1 for all the selected cases, whereas the di¨erences show the presence of nonlinearities related to the e¨ect of the injector pressure drops. In particular, for each of the fuel pressure drops, the response is about 20% higher at higher oxidizer pressure drops and is nearly linear in the lower oxidizer pressure drops domain. Furthermore, this trend indicates that the relative increase with respect to the reference test case with ǫ = 0.01 is not a¨ected by the amplitude of the input noise (both for ǫ = 0.012 and 0.02). Considering as example the points with -p ox /p c = 0.8 and -p fu /p c = 0.4 in Fig. 7c, one may ¦nd that the absolute root mean square values are 1.34 and 2.47 bar, respectively, for ǫ = 0.012 and 0.02, both showing a linear trend with respect to the 1.08-bar value obtained in the ǫ = 0.01 case. Therefore, the system shows that the pressure response to the input noise is nearly linear. Nevertheless, there are no experimental measures in the referenced paper for RMS values below 6 bar (≈ 13% of the operative pressure) (see Fig. 3), thus not providing any validation point for the model.
However, it must be pointed out that the amplitude of pressure oscillations is such that nonlinear e¨ects could be activated for these levels of input noise.
Furthermore, the results shown in Fig. 7 highlight that for di¨erent fuel injector pressure drops, the system response has a maximum at higher oxidizer injector pressure drops. Therefore, when the amplitude of the input perturbation is increased, the boundary of the stability regions experiences a greater shift if placed in the high-pressure drop region and a smaller shift if placed in the low-pressure drop region. It should be noted that increasing the input amplitude does not a¨ect the characteristic frequency of the system that remains equal to 60 Hz.

Variable Time Lag
The time lags are typically not known a priori; hence, when a new engine con-¦guration is analyzed, it is useful to retrieve the needed time delays by means of dedicated correlations. The bene¦t of introducing correlations for the time lags is dual: on one hand, the user is free from any a priori calculation of the time lags and, on the other hand, the time lags become time-varying functions of the combustion chamber variables. In this frame, we introduced semiempirical correlations to calculate atomization and vaporization time lags, considering the combustion as su©ciently fast. As mentioned in paragraph 2.1.3, the time lags for atomization and vaporization processes are retrieved, respectively, from Eqs. (3) and (2), while mixing time lag is evaluated from an empirical curve [14]. A double time lag approach is still used; hence, the total time delay for the liquid oxygen is retrieved as the sum of vaporization, atomization, and mixing time lags while only the mixing time lag characterizes the gaseous hydrogen delay. Some considerations must be done on the correlation for the initial droplet diameter Eq. (4). Due to the lack of informations on the experimental setup, the recess of the oxidizer post R post is an unknown data that greatly in §uences the value of the initial droplet diameter and, in turn, the vaporization time. Furthermore, the injector coe©cient C inj can range from 0.4 to 1.2, thus changing signi¦cantly the vaporization delay. We know, from experimental evidences [28], that the initial droplet diameters of liquid oxygen at high pressure range from 60 to 160 µm. Hence, stated that the recess length is a construction parameter of a few millimeters, we decided to use a combination of values that, while respecting all of the aforementioned constraints, ensures a vaporization time lag at steady state equal to τ vap = 4.4 ms. In this way, we retrieved a time varying time lag matching with the steady-state values suggested from the experimental evidences [22]. The chosen reference working point is characterized by the values R post = 0.003, C inj = 0.4, and K = 0.121.
The time-dependent approach is then tested on the engine con¦guration presented in section 3. A single fuel injector pressure drop is here investigated at di¨erent oxidizer injector pressure drops and results are compared in Fig. 8 against those obtained for the constant time lag approaches. The variable time lag approach gives a pressure response in line with the ¦xed time lag approaches. The amplitudes of the oscillations are less in §uenced by the reduction of the oxidizer injector pressure drop and a smoother growth is observed while passing from high to low injector pressure drops. Finally, the time-dependent time lags give a characteristic frequency of 65 Hz approaching the experimental measured value of 68 Hz, with an uncertainty of 5 Hz related to the broadband sampling frequency. Furthermore, the tendency of the double time lag approach in anticipating the insurgence of the inherent stability limit is here con¦rmed. In fact, in both double time lag approaches, the limit moves towards higher oxidizer pressure drops and for pressure drops lower than -p ox /p c = 0.23 the system is not able to damp the imposed perturbation causing the chamber pressure to diverge with a frequency equal to the characteristic frequency of 65 Hz.

CONCLUDING REMARKS
A double time lag model to investigate low-frequency combustion instabilities has been implemented in the system analysis tool EcosimPro and tested against experimental results. Two di¨erent approaches have been followed: a constant time lag and a variable time lag approach based on semiempirical correlations. The constant time lag approach demonstrates its e¨ectiveness in reproducing the unstable behavior associated with the reduction of the injector pressure drops, in particular, no stable region has been identi¦ed for an oxidizer injector pressure drop lower than 0.35. The chug module is able to reproduce the reversal in slope of the stability boundary at low injector pressure drops, thus giving a result close to the one found in the reference experiment. However, it must be pointed out that the location of the stability boundary depends on both the criterion used to identify a regime as unstable and on the amplitude of the input perturbation. It is demonstrated that the boundaries can be easily shifted by changing the amplitude of the input perturbation and that this shift slightly depends on the selected injector pressure drop. Hence, in the computer program, the stability boundaries can be shifted in accordance with the level of the input perturbation. Despite its simplicity, the constant time lag approach is able to predict a characteristic frequency close to the one measured experimentally.
The registered high-frequency content is characterized by a ratio of higher to lower frequency close to 2.7, as reported in open literature. It must be noted that in the described case, high frequency does not dominate the dynamic response at low fuel injector pressure drops and the instability is still related to the ¦rst characteristic frequency. The inability in catching the suppression of the ¦rst mode could be related with a lack of data for parameters that usually in §uences the damping e¨ect and the associated system impedances, i. e., injector cavity volumes and length of the feed lines. These values have been selected to be consistent with the downstream injector geometry and mass §ow rates, but a ¦ne tuning with the real data is unavoidable to improve the results of the model. The chugging module has been ¦nally tested with semiempirical correlations to retrieve characteristic delays that, depending on the main §ow variables, are time-dependent functions of the main combustion chamber parameters. With the implemented correlations, the model retrieves a lower characteristic frequency that better approximates the experimental results, while the predicted RMS values of the chamber pressure remain comparable with those obtained with the constant double time lag model.