Numerical Analysis of Self-Excited Combustion Instabilities in a Small MMH/NTO Liquid Rocket Engine

Combustion instabilities in a small MMH/NTO liquid rocket engine used for satellite attitude and course control are numerically investigated. A three-dimensional Navier-Stokes code is developed to simulate two-phase spray combustion for cases with five different droplet Sauter Mean Diameters. As the droplet size increases from 30 microns to 80 microns, pressure oscillations are stronger with larger amplitudes. But an increase of the droplet size in the range of 80 microns to 140 microns indicates a reduction in the amplitudes of pressure oscillations. This trend is the same as the Hewitt criterion. The first tangential (1T) mode and the first longitudinal (1L) mode self-excited combustion instabilities are captured in the 60-micron and 80-micron cases. Abrupt spikes occur in the mass fraction of MMH and coincide with abrupt spikes in the mass fraction of NTO at the downstream regions just adjacent to the impinging points. Thus, local combustible high-dense mixtures are formed, which result in quasiconstant volume combustion and abrupt pressure spikes. The propagation and reflection of pressure waves in the chamber stimulate the combustion instability. When the droplet size is too small or too large, it is difficult to form local high-dense premixtures and combustion is stable in the chamber.


Introduction
A small MMH/NTO liquid rocket engine (LRE) is developed for the purpose of attitude control and course corrections for tactical missiles and artificial satellites. High-frequency combustion instability has plagued the development of LRE since the 1940s, in which organized oscillations with high amplitude greater than 10% of mean chamber pressure would be induced [1,2]. Such oscillations may result in poor performance, unacceptable vibrations, or even catastrophic events. Thus, combustion instability must be considered in a priori. Because of the complexity of this problem and limited computational power to capture all the physical processes in the numerical simulations, combustion instability is still a challenge to the development of a rocket engine.
Hot fire tests were conducted to study combustion instability in the chamber [3][4][5][6]. But the quantity of tests is needed due to the notoriously probability and uncertainty of combustion instability. Thus, the prediction of combus-tion instability by numerical simulation is a major issue. Various levels of models have been developed to deal with combustion instability in the chamber. Wave equations [7][8][9][10], Finite Element Method (FEM) [11], and Linear Euler Equations (LEE) [12,13] have been used to predict combustion instability. But transfer functions [14,15] need to be provided and calibrated with high fidelity data. Computational Fluid Dynamics (CFD) method is regarded as an effective way to investigate combustion instability, which can provide insight into more detailed physical processes. Large Eddy Simulation (LES) is a good candidate, and it has been used to predict combustion instability successfully [16][17][18][19][20][21]. Detached Eddy Simulation (DES) featuring little demands on computational resources has been developed to investigate combustion instability [22][23][24]. However, liquid droplet dynamic was not considered and the propellants were gasified when they left the injectors in these simulations. For two-phase turbulent reactive flow, Unsteady Reynolds-Averaged Navier-Stokes (URANS) is always employed [25,26] and has been applied to study combustion instability in LRE [27][28][29]. Tucker et al. [30,31] compared LES and URANS simulations of a combustor. The results showed the highest-fidelity LES and URANS best predicted the instantaneous temperature, vorticity fields, and wall heat flux, while the two lower-order LES provided a poor prediction. Therefore, URANS is employed in this paper due to the enormous demands of computer resources for DES or LES considering two-phase interactions.
As is well known, high-frequency combustion instability results from the coupling between unsteady heat release and acoustic pressure oscillations. Many researches focused on the simulation of unsteady heat release [17-20, 22, 32] and a heat release response function [33]. Raleigh index [34] was used to quantitatively evaluate the coupling between unsteady heat release and acoustic pressure oscillations, which represented a source term in the balance of acoustic energy [17]. When combustion instability occurs, the Raleigh index is greater than 0. But when the Raleigh index is greater than 0, there are cases shown that combustion instability does not occur. Moreover, it says nothing about how the initial heat release oscillation is produced and what is the source of pressure oscillation. High-frequency combustion instabilities can be classified as intrinsic and injectioncoupled instabilities [2]. Injection-coupled instabilities arise when injected flow fluctuations couple with chamber acoustic modes, which mainly occur in LRE equipped with coaxial injectors. Mass flow rate oscillations in the injectors are the sources of heat release oscillations [35]. It was found that the acoustic mode of the chamber was coupled with that of injectors when combustion instability occurred [36][37][38]. For the cases in this work, like-doublet injectors are used and the mass flow rate is constant without oscillation. Intrinsic instabilities are governed by the subprocesses that occur following propellant injection, atomization, vaporization, mixing, and chemical reaction, in which injectors are passively involved. Sirignano and Duvvur [39,40] summarized that the vaporization process was a key factor in combustion instability in the liquid rocket chamber. Lei and Turan [41] declared the evaporation process had a great influence on the nonlinear and intensity behavior of combustion instability. The vaporization rate depends on the flow field around the droplet and droplet size. Hewitt criterion [42,43] was used to predict the highest mode of combustion instability in the chamber equipped with impinging injectors. Combustion instability is related to a stability correlating parameter d 0 /U j , where d 0 is the injector's orifice diameter and U j is the injection velocity of the least volatile propellant. A reduction in d 0 /U j , indicating decreased stability margin, coincides with a reduction in the droplet size. Therefore, the droplet size may be a key factor in combustion instability. Few studies focused on the effects of droplet sizes on combustion instability in small thrust LREs. On the other hand, in previous studies, how the initial pressure oscillation is produced and what is the source of pressure oscillation are still not clearly explained.
In this paper, a robust three-dimensional two-phase reaction flow model based on URANS is developed. The spray combustion process for cases with different droplet sizes in a small MMH/NTO LRE are investigated. Combustion instabilities are captured, and how the initial disturbance is induced and evolved combustion instability is analyzed. The effects of droplet sizes on combustion instability are discussed. The structure of this paper is in the following manner. First, a numerical model is described briefly in the next section followed by numerical method and boundary conditions. Then, physical model and calculation cases are given. Finally, the results of these cases are presented, in which the instability phenomenon is performed. The instability mechanism is analyzed further. And the effects of droplet sizes on combustion instability are discussed later.

Numerical Model
In a small MMH/NTO LRE, the propellants undergo injection, atomization, vaporization, mixing, chemical reaction, and expansion process. The Eulerian-Lagrangian method is employed to solve two-phase reactive flow in the combustion chamber, in which the Navier-Stokes (N-S) equations and the Lagrangian equations are for the gas phase and the liquid particles, respectively. The interactions between the two phases are mathematically expressed in source terms in the gas phase equations.

Continuous Phase Model.
Three dimensional, multispecies gas-phase Navier-Stokes equations are used to describe the flow in the chamber. The standard k − ε two-equation turbulence model is used. The governing equations of mass, momentum, energy, species mass fraction, turbulence kinetic energy, and turbulent dissipation can be expressed in the unified form [44]: where φ, Γ φ , S φ , and u j represent conservation variables, diffusion coefficients, convective source terms, and the velocity in three directions. S pφ and S cφ are the source terms determined by droplet vaporization and chemical reaction. The detailed variables are presented in Reference [44].

Atomization
Model. The atomization of a liquid impinging jet is considered through a correction formula, by which spray properties are directly determined [45]. The droplets are introduced into the chamber with a N-T distribution [45] at given locations, where are the locations of impinging points of the like-doublet injectors. Initial velocities, mass flow rates, and Sauter Mean Diameters (SMD) of the droplets are given as boundary conditions. The droplets are injected in a fan or cone angle in spatial distribution. The initial velocity vector of droplets is determined by a random process. The number of parcel injected in each time step is obtained by the mass flow rate dividing the mass of droplet particle. International Journal of Aerospace Engineering droplets and the surrounding gas are included in the source terms in gas-phase equations. The trajectories of discrete droplets are predicted by integrating the force balance on the droplets in a Lagrangian frame. The droplet gravity and other forces are neglected. Only the drag force is considered. The momentum equations for a droplet are as follows: where the subscript p represents the droplet. The subscript j denotes droplet species. C D,j is the drag coefficient of the droplets. u ! is the mean gas velocity; u ! p,j is the velocity of the j droplet. u ! ′ is the turbulence velocity fluctuation.
The source terms for the two-phase interactions are described in detail in Reference [44].

Chemical Reaction Model.
To make this problem simple, only the gas-phase chemical reaction is considered, while the liquid-phase reaction is neglected. The kinetically chemical reaction rate is calculated by the Arrhenius expression of four-step chemical reactions. The MMH-NTO four-step chemical reactions are expressed:

Computational Methodology
The numerical methods employed to solve this problem include a temporal discretization and a spatial discretization. The finite volume method is employed to discrete Equation (1) in space, while the Euler scheme is used for the time terms. The diffusion terms and convection terms are approximated by a central difference scheme and a second-order upwind scheme, respectively. First, the source terms due to evaporation and two-phase interaction are calculated in a Lagrangian phase. Then, the field coupling pressure with velocity is solved by the Semi-Implicit Method for Pressure Linked Equations (SIMPLE) method. The convection fluxes are calculated due to the movement of mesh in the Euler phase finally. This is the so-called Arbitrary Lagrangian-Eulerian (ALE) technique [46]. The grid size is 2 mm. Further refined grid does not improve the results significantly.
Adiabatic wall conditions are applied for the chamber walls including the injector face. Wall function treatment [47] is employed to the zone near the wall. The characteristic boundary (outflow boundary) is used at the exit, where the back pressure is specified and the velocity components are set equal to those of the logical inside neighbor vertexes due to supersonic flow at the exit.

Results and Discussion
4.1. Physic Model. The schematic diagram of the small MMH/NTO thrust chamber is presented in Figure 1. The diameter of the chamber is 60 mm, while the diameter of the throat is 29.3 mm. The length of the cylindrical part is 59.5 mm, while the total length of the chamber is 146.8 mm. The nozzle expansion ratio is 1.55. There are 60 likedoublet injectors and 10 cooling injectors shown as Figure 1(a). The like-doublet injectors are in the pattern of fuel-oxidizer-fuel-oxidizer. The atomization process is neglected, and the droplets are injected in the impinging points shown as Figure 1(b). The like-doublet injectors (impinging points) are distributed in two rings. The first ring is located at r = 10:3 mm, while the second ring is located at r = 19:4 mm. 20 like-doublet injectors are uniformly distributed in the first ring every 18 degrees in the circumference direction, while 40 like-doublet injectors are uniformly distributed in the second ring every 9 degrees. The cooling injectors are located at r = 23 mm and distributed every 36 degrees. 30% of the fuel injected by the cooling injectors are used to cool the chamber wall. The mixture ratio of NTO to MMH (O/F ratio) is 1.65, and the total mass flow rate is 484.5 g/s. The initial temperature of MMH and NTO droplets is 300 K. The spray cone angle is 75 degrees. The injection velocities of MMH and NTO are 23.5 m/s and 28.5 m/s, respectively. The mean chamber pressure is 1.09 MPa. At the initial time, the chamber is filled with N 2 , while the chamber pressure and chamber temperature are 1.09 MPa and 3000 K, respectively.

International Journal of Aerospace Engineering
The theoretical characteristic frequencies of the chamber can be estimated by simplifying the chamber as a cylinder with the same diameter and equivalent length, which equals the sum of the length of the cylindrical part and two-thirds of the contraction part of the chamber [5]. The detailed values are presented in Table 1. The sound speed is taken as 1020 m/s according to the NASA-CEA.

The Hewitt
Criterion. The stability correlating parameter d o /U j had been successfully used to predict combustion instability in the combustor with impinging jet injectors. The Hewitt Stability Correlation shown as Figure 2 presented the stability margin with d o /U j . d o is the injector's orifice diameter, while U j is the velocity of the least volatile propellant. A reduction in d o /U j indicates a decreased stability margin. It can be concluded that when f d o /U j is less than 0.1, combustion instability may occur.
Hot fire tests for the small MMH/NTO thrust chamber shown as Figure 1 are conducted to study whether the Hewitt criterion is suitable for the combustor. The parameters of the combustor and injectors are the same as those in Section 4.1, while the mass flow rate and parameter of propellant droplets are different. The parameters d 0 /U oðf Þ are controlled by adjusting velocities of fuel (oxidizer). Thus, different d 0 / U oðf Þ and the mass flow rate of fuel (oxidizer) are obtained. Pressure oscillations are monitored at the head of the chamber. The stable and unstable cases are performed in Figure 3 with different d o /U j . It can be found that the Hewitt criterion is applicable to this combustor. When the parameter d o /U j is less than a determined value, combustion instability occurs.
Anderson et al. [42,43] declared that an increase in d o /U j implied an increase in the mean droplet size. In the atomization test for like-doublet injectors with water, it can be concluded that the SMD of the droplets can be expressed as follows: where the parameter σ is surface tension and θ is the spray cone angle. The parameter ρ p is the density of the droplet. The SMD is proportional to d o /U j . Thus, when the SMD is less than a determined value, combustion instability may occur.
According to the Hewitt criterion, the first longitudinal (1L) and the first tangential (1T) combustion instabilities would occur when d o /U j is less than 1e-5. It can be estimated that when d o /U j equals to 1e-5 according to Equation (4), the SMD is about 80 microns. Thus, the SMD of droplet sizes investigated for MMH and NTO are 30 μm, 60 μm, 80 μm, 100 μm, and 140 μm. The detailed cases are presented in Table 2.

The Instability Phenomenon.
In the hot test, pressure oscillations are monitored at Point A (28 mm, 0 mm, and 2 mm), which is near the chamber sidewall and shown in Figure 1(b). This point is located at the antinode of 1T mode and 1L mode. The test condition is the same as the simulation condition in Section 4.1, while the mean diameter and distribution of droplets are unknown. The FFT analysis of pressure oscillations is shown in Figure 4. There are two peak frequencies 4800 Hz and 10780 Hz.
Pressure oscillations in Case 1-Case 5 monitored at the same position (Point A) with the hot test are presented in Figure 5, which give an indication of the trends of time histories of pressure oscillations in the chamber including any growth or decay. The pressure data over time slice 5 ms to 10 ms are taken as FFT analysis, which are performed in Figure 6. The resonant frequencies of pressure oscillations are obtained by FFT analysis. There are also two peak frequencies around 4800 Hz and 9800 Hz in Case 1-Case 3. Compared with theoretical acoustics frequencies of the chamber, these resonant frequencies in the hot test and simulation results can be identified. The peak frequency around 4800 Hz is identified as 1L mode, while the peak frequency around 9800 Hz is identified as 1T mode. It is shown that the acoustic modes excited in the numerical simulation are in agreement with those in the hot test, which indicates the numerical method can predict acoustic modes of excited pressure oscillations in this paper. There are obvious differences in the amplitudes of acoustic modes due to some differences in the numerical simulation and hot test. Firstly, the atomization process is neglected in the numerical simulation. Secondly, no every droplet can be traced due to the limit of computational resources. Droplet parcels are employed to represent the quantity of droplets with the same trajectory. Thirdly, a global one-step chemical reaction is used in the simulation, while it is multistep complicated chemical reactions in the hot test. Moreover, it is a challenge for accurate measurements under a high-pressure and hightemperature condition in the hot test. In this paper, it focuses on whether combustion instability is acoustic combustion instability (high-frequency combustion instability). Therefore, the acoustic modes of pressure oscillations and the frequencies of excited acoustic modes are compared in the validation. It is difficult to compare their amplitudes quantitatively.
There are several features to note in Figure 5. The value of pressure oscillations is less than 0 at the initial time, although the initial chamber pressures are set as the mean chamber pressure 1.09 MPa. Because there is a time lag between the time when the propellant enters into the chamber and the time when they are burned and release their chemical energy, the pressure will drop at the initial time, which is the first feature of pressure oscillations. Upon initial, pressure oscillations rise to 0.3 MPa due to autoignition in Case 1, which is higher than that in Case 2. However, pressure oscillations drop to lower values in Case 3-Case 5. The pressure firstly declines to a lower value and then rises to a higher value from In Case 1, pressure oscillations decay within the amplitude of 10% of mean chamber pressure. The combustion is considered stable. For Case 2, the peak-to-peak amplitude of pressure oscillations steadies out approximately 0.7 MPa (64.2% of mean chamber pressure). Pressure waves are steep-fronted and not symmetric due to the nonlinear effect, which is the second feature of pressure oscillations. For Case 3, the peak-to-peak amplitude of pressure oscillations is 0.8 MPa (73.4% of mean chamber pressure). There exist many abrupt pressure spikes. For Case 4, pressure oscillations decline to -0.3 MPa at 0.5 ms and rise to 2.74 MPa at 1.75 ms. At 5.0 ms, pressure oscillations damp out to the steady mean chamber pressure. Thus, combustion instability is not triggered. In Case 5, the pressure drops to 0.55 MPa      There not exist high-amplitude pressure peaks. In Case 2 and Case 3, pressure oscillations are organized and periodic, which is the third feature of pressure oscillations. And the amplitudes are greater than 10% of mean chamber pressure. The 1T mode and 1L mode are excited. The amplitudes of 1L mode and 1T mode in Case 3 are the largest. Thus, combustion instability occurs in these two cases. Although the 1L mode and 1T mode are also excited in Case 1, the amplitudes of pressure oscillations are less than 10% of mean chamber pressure. Therefore, combustion instability is not excited in this case.
As the droplet size varies from 30 μm to 140 μm, combustion varies from stable to unstable, and to stable again. The 30 μm and 100 μm cases initially display high amplitude of pressure oscillations which damp over time and then keep stable. The 140 μm case displays no instability over 10 ms run time. There are obvious periodic high-amplitude pressure oscillations in 60 μm and 80 μm cases, and those in 80 μm are most violent. The amplitudes of pressure oscillations increase firstly and then decrease in the calculated range of droplet size. And the same trends are for the amplitude of 1L mode and 1T mode. Similar results were also provided by Keenan et al. [48] and Kim et al. [49]. It indicates the results are reasonable in this paper. Keenan et al. [48] investigated the effects of 20 μm, 60 μm, 100 μm, and 140 μm droplet sizes of LOX and RP-1 on stability. The 20 μm and 60 μm displayed transverse instability, while the 100 microns initially displayed transverse instability but damped over time. The 140 microns showed no instability. Kim et al. [49] found that droplets on the order of 10 microns were stable to bombing, while 50-micron droplets were not stable. But the 100micron droplets were stable to a simulated triggering mechanism.

Abrupt Pressure Spikes and Instability Mechanism.
In order to reveal how combustion instability occurs, Case 2 is analyzed in detail. Figure 7 shows several instants of pressure contours of the chamber in Case 2, which perform the propagation of pressure waves in the chamber. The transverse sections are 2 mm, 4 mm, 6 mm, and 20 mm away from the injector face, respectively. The acoustic modes of oscillations can be revealed by the patterns of time evolution of pressure field. Small local high-pressure regions appear at the head of the chamber shown as Figure 7(a). Pressure waves propagate from a high-pressure region to a low-pressure region. In the transverse section, the center of the section becomes a lowpressure region, while the region near the sidewall becomes a high-pressure region shown in Figure 7(b). And the pressure at the right part is high than that at the left part, which is the typical pressure distribution of the 1T mode. Pressure waves reach the contraction section. The head region of the chamber becomes a low-pressure region, while the contraction section becomes a high-pressure region shown as Figure 7(c). Pressure waves are reflected and propagate to the head region of the chamber. And the head region of the chamber becomes a high-pressure region shown as Figure 7(d). It can be concluded that the 1L mode is excited. Therefore, the 1T mode and 1L mode are stimulated in Case 2, which indicates the identification of the peak frequencies of the FFT results is correct. Moreover, it is worth noting that small local high-pressure regions appear at the transverse sections 2 mm and 4 mm occasionally shown as Figures 7(a) and 7(d), where are at the downstream region just adjacent to impinging points, which may be the source of pressure oscillations.
In order to reveal how the small local high-pressure regions are formed, a monitored point named Point B is set there. Point B is at the location (10 mm, 0 mm, and 2 mm) shown as in Figure 1(b), where is the downstream region just adjacent to the first ring like-doublet injectors' impinging point. Point A is at the location (28 mm, 0 mm, and 2 mm), at the same transverse section with Point B, but near the sidewall. Figure 8 shows pressure oscillations at Point B are obviously stronger than those at Point A. The peak-to-peak amplitudes of pressure oscillations are approximately 1.7 MPa at Point B, but only 0.6 MPa at Point A. Furthermore, pressure peaks at Point B appear earlier than those at Point A shown as Figure 8    International Journal of Aerospace Engineering Then, the time histories of pressure oscillations of the surrounding points of Point B are observed. Figure 9 presents pressure oscillations from the time 8.320 ms to the time 8.360 ms, which give an indication of that where the pressure spike is induced. The back point is at the injector face, not presented. The amplitude of pressure oscillation at Point B is the highest. An abrupt pressure peak is observed in a period from t 1 = 8:33830 ms to t 2 = 8:33906 ms at Point B. The t 1 and t 2 are marked by squares in Figure 9. For the short interval Δt = 0:76 μs, the pressure rises from 3.14 MPa to 4.58 MPa in this control volume. The pressure increase in the observed control volume is 1.44 MPa, and it is nearly   7 International Journal of Aerospace Engineering 50% of the pressure at t 1 , which indicates a pressure spike occurs. The pressure, temperature, density, density of MMH, evaporation amount of MMH, and combustion amount of MMH at the two instants in the Point B control volume and surrounding control volumes are shown in Table 3. The combustion amount of MMH at t 2 is the largest in the observed control volume B, while the evaporation amount of MMH at t 1 is the largest. This means violent combustion occurs from t 1 to t 2 . In conclusion, the pressure peak at Point B in this period cannot result from the propagation of pressure waves of surrounding points. It may be caused by violent combustion.
The total heat release produced by a chemical reaction is 0.85 J, and the heat release rate is about 1.12e6 w. The flow velocity in the control volume is 93 m/s. During the short interval, the corresponding propagation distance of fluid flow is 0.078 mm. Compared with the characteristic size of the control volume 2 mm, convection fluxes with the surrounding vertexes can be neglected. The heat release due to combustion is consumed by the temperature rise of mixture, evaporation process, and two-phase interactions in the control volume. Base on the constant volume combustion theory, the temperature rise of the mixture in the control volume is 585 K, while that predicted by the numerical simulation is 538 K. In the theoretical calculation, the energy used for two-phase interactions is not considered. Thus, the temperature increase of the mixture is a little higher than that in the numerical simulation. The average molecular weight varies from 34.6 to 29.7. According to the constant volume combustion theory, the pressure at t 2 can be estimated. It is 4.47 MPa, which is comparable to the result of 4.58 MPa in the numerical simulation. The comparison of results obtained by constant-volume combustion theory and numerical simulation results is shown in Table 4. For the short interval Δt = 0:76 μs, the propagation distance of pressure waves is about 0.875 mm with sound speed 973 m/s. The pressure peaks cannot propagate to surrounding control volumes in this short interval. The density in the control volume is 4.93e-3 g/cm 3 at t 1 , while it is 5.14e-3 g/cm 3 at t 2 . Therefore, the density can be regarded unchanged, while the rise of pressure is proportional to the rise of temperature. The above analysis shows that quasiconstant volume combustion happens in a short period, which causes the abrupt rise of the pressure.
The consumed amount of MMH due to combustion is greater than the amount of MMH at t 1 . Considering the propagation distance in this short interval, the fuel gas and oxidizer gas cannot flow from the surrounding control volumes. Thus, the evaporation amount has a great influence on combustion. Figure 10 shows the evaporation amount and consumed amount of MMH in the period from 8.335 ms to 8.345 ms. The abrupt enhancement of evaporation is followed by the enhancement of combustion. The abrupt rise of MMH consumed amount in the short interval means that violent combustion occurs. The time histories of mass fraction of fuel and oxidizer are given in Figure 11.  gas. The mass fraction of fuel gas and oxidizer gas decreases sharply from 8.338 ms to 8.339 ms, which indicates that violent combustion occurs and the reactant gas is burned. Thus, violent combustion occurs in the period t 1 to t 2 . Pressure waves cannot propagate to the surrounding control volumes in this short period, which leads to a high-amplitude increase of pressure, namely, pressure spikes. These results are also supported by the work of Keenan et al. [48], Zhang et al. [44], Grenda et al. [50], and Daimon et al. [51], which indicates violent combustion (bombing) may occur in LRE. Keenan et al. [48] found that injecting the LOX and RP-1 droplets along the same vector in a given angel to mimic fan atomization may trigger self-excited instabilities, because the LOX and RP-1 droplets can interact more rapidly. Zhang et al. [44] found that quasiconstant volume combustion resulted in bombing, which induced self-excited combustion instability. Grenda et al. [50] found that the simulated bombing occurred when the droplets were    [51] investigated the explosion mechanism of hypergolic propellant and found that the instantaneous gasification of the propellant resulted in the explosion.
As Figure 7 shows, the local bombing phenomenon and abrupt pressure spikes appear in the chamber, mainly the downstream region just adjacent to the like-doublet injectors impinging points. Furthermore, the spatial and temporal distribution of pressure spikes is analyzed. Four lines A, B, C, and D parallel with the axis of the chamber are selected, which locate near the axis of the chamber, the inner likedoublet injector impinging points, the outer like-doublet injector impinging points, and the chamber sidewall, shown as red lines in Figure 1(b). The time histories of the pressure of six control volumes with space about 2 mm in the axial direction are presented in Figure 12. For spatial distribution, the pressure spikes mainly appear at line B and line C within 4.1 mm away from the injector face, where are the corresponding downstream regions just adjacent to the likedoublet injectors impinging point. The fuel and oxidizer droplets almost disappear at the head of the chamber in the 60-microns case. The combustion mainly occurs at the head of the chamber, where are the downstream region just adjacent to the impinging points. For temporal distribution, the pressure spikes come up frequently and stochastically over the 10 ms run time. Thus, pressure spikes occur at the downstream region just adjacent to the like-doublet injectors frequently and occasionally. The propagation and reflection of pressure waves would stimulate periodical pressure oscillations. And the characteristic frequencies of pressure oscillations are equal to the theoretical characteristic frequencies of the chamber. Thus, combustion instability would be stimulated in the chamber.

Effects of Droplet
Sizes on the Distribution of Abrupt Pressure Spikes. As shown in Section 4.3, the droplet size has a great influence on the amplitudes of pressure oscillations, which increase firstly and then decrease from 30 μm to 140 μm. Self-excited combustion instabilities are stimulated in the cases with 60 μm and 80 μm, while other cases are stable. In the detailed analysis of 60 μm case, it is found that quasiconstant volume combustion causing abrupt pressure spikes is the source to drive and sustain combustion instability. Therefore, the effects of droplet sizes on the distribution of abrupt pressure spikes are discussed in this section. Figures 13-17 show the acoustic pressure contours at two instants during 8 ms-10 ms in Case 1-Case 5, which provide an insight into the distribution of abrupt pressure spikes in the cases with different droplet sizes. The transverse sections are successively 2 mm, 4 mm, and 6 mm away from the injector face from left to right in Figures 13-17. The distribution of MMH and NTO droplets is shown in Figures 13(a)-17(a) and Figures 13(b)-17(b), respectively.
There are several features in these figures. Firstly, the MMH and NTO droplets are evaporated and disappear at the head of the chamber in Case 1-Case 3. The droplets go through the cylindrical part of the chamber in Case 4, while the droplets reach the throat of the thrust chamber in Case 5. As the droplet size increases, the distances the droplets go through are longer, because more time is needed to heat up the droplets and the evaporation rate is smaller with the increase of the droplet size. Furthermore, the pressure is almost uniform in the chamber in Case 1, Case 4, and Case 5, while there are large pressure gradients in Case 2 and Case 3. The pressure at the head of the chamber is obviously higher than that at the downstream region of the chamber, which implies that 1L mode combustion instability may occur. For transverse sec-tions, there are small high-pressure regions (pressure spikes) in Figures 14 and 15, whose amplitudes are greater than 50% of mean chamber pressure. The pressure spikes mainly occur at the head of the chamber. They come up at the downstream region of the impinging points shown as the transverse section in Figures 14 and 15    11 International Journal of Aerospace Engineering in these regions. For Case 1, the amplitudes of the red regions (high-pressure regions) are less than 10% of mean chamber pressure. It can be considered that there are not pressure spikes after the pressure spikes due to autoignition decay. There no exist obvious pressure spikes in Case 4. For Case 5, the pressure keeps stable without abrupt pressure spikes. The propellant droplets may go through the throat, which makes that mass flow rate of the throat section is not constant. Thus, mean chamber pressure changes slightly without strong pressure oscillations. Abrupt pressure spikes appear at the head of the chamber and combustion instabilities are induced with SMD 60 μm and 80 μm. When the SMD of droplets is too small or too large, abrupt pressure spikes are not observed and combustion is stable. It indicates combustion instabilities are involved with abrupt pressure spikes. And only when pressure spikes occur, combustion instabilities are excited. This conclusion is consistent with the instability mechanism. Pressure spikes may be the sources to drive and sustain combustion instability.

The Condition of the Formation of Abrupt Pressure
Spikes. According to Section 4.5 in the results and discussions part, it can be known that combustion instabilities are related to the occurrence of pressure spikes, which are sensitive to droplet sizes. In order to reveal how droplet sizes affect the occurrence of pressure spikes, the mass fraction of fuel and oxidizer gas at different radius with the distance of 2 mm away from the injection face is presented in Figures 18-22. The red lines represent the NTO mass fraction, while the black lines represent the MMH mass fraction. At r = 6 mm, the mass fraction of MMH is greater than that of NTO, while at r = 10 mm-18 mm, the mass fraction of NTO is greater. It is fuel rich for r = 22 mm, because the cooling injectors are located there and spray MMH. Moreover, the quantity of mass fraction peaks of MMH and NTO occurs at r = 10 mm and r = 18 mm, where are the downstream regions just adjacent to impinging points.
In Case 1, the oscillations of the mass fraction of MMH and NTO are not ordered with a small amplitude.

International Journal of Aerospace Engineering
Abrupt spikes arise at the time series of mass fraction of MMH accompanied with abrupt spikes at the time series of mass fraction of NTO. Thus, local premixtures are formed with a high density of reactant gas. Quasiconstant volume combustion happens there, which results in a high-amplitude rise of pressure. The abrupt pressure peaks arise frequently and stochastically. For Case 4, the oscillations of mass fraction of MMH and NTO are ordered over time slice 2 ms to 5 ms. Pressure oscillations over this time slice also display organized shown as in Figure 5. After 5 ms, the oscillations of mass fraction of MMH are out of phase with those of mass fraction of NTO. For Case 5, there no exist abrupt peaks in Figure 22. The mass fraction of MMH and NTO changes slowly with time.
Abrupt spikes occur in the mass fraction of MMH and NTO at the same time frequently and stochastically at the downstream regions just adjacent to impinging points. Local combustible mixtures are formed there, which lead to quasiconstant volume combustion and abrupt pressure peaks. Their propagation and reflection in the chamber wall may excite self-triggered combustion instability.

Effects of Droplet
Sizes on Combustion Instability. In Case 1 to Case 3, abrupt pressure spikes come up more frequently with a higher amplitude, as the propellant droplet size increases. The evaporation of the droplets is faster with smaller diameters. The diffusion combustion is performed. The local mixed fuel and oxidizer gas are formed less frequently with smaller density. There are more abrupt peaks of the mass fraction of fuel and oxidizer with an increase in the droplet size in Case 1 to Case 3 shown in Figures 18-20. However, pressure spikes occur less frequently with a lower amplitude in Case 4. The vaporization of propellant droplets is slow because of large diameters. The peaks of the mass fraction of fuel and oxidizer occur less frequently. The mixture ratio O/F is large and it is hard to form local highdense premixtures. No pressure spike exists after 5 ms. For Case 5, the fuel droplets appear at the throat part of thrust chamber. The densities of reactant gas are low at the head of the chamber, and local high-dense premixtures cannot be formed.
The trend is the same as the Hewitt Stability Correlation. A reduction in the droplet size of the range of 140 microns to 14 International Journal of Aerospace Engineering 80 microns, indicating an increased tendency towards instability, implies a reduction in the stability parameter d o /U j in the Hewitt Stability Correlation. In the actual test, the mean droplet size is greater than 30 microns. The 30-microns droplet size case does not conclud in hot test. The fuel and oxidizer droplets are injected along the same vector in a given angle, which result in that they are distributed in the chamber similarly. When the evaporation of fuel and oxidizer droplets is enhanced and mass fraction of fuel and oxidizer rise to a peak at the same time, local premixtures are formed easily with high density. This results in pressure spikes. The spikes in fuel and oxidizer mass fractions are formed by a combination of three factors of rapid evaporation, high density of droplets, and well mixing of fuel and oxidant. The droplet sizes determine the evaporation rate and number density of the propellant droplets, which further affect the formation of the local premixtures. When the droplet size is too small, the evaporation rate and initial number density are large. The fuel and oxidizer gas are burned out as soon as they are gasified from the fluid phase. And fuel and oxidant do not mix very well for a short existing time of fuel droplet and oxidant droplet. Therefore, the abrupt peaks of mass fraction of reactant gas would not exist. It is difficult to form high-dense premixture. When the droplet size is too large, the evaporation rate and initial number density are small and mass fraction of reactant gas will not vary greatly. It is also difficult to form high-dense premixture. Only when the droplet size is in a definite range, the abrupt peaks of mass fraction of reactant gas will occur stochastically and frequently at the same time. And it is easy to form highdense premixture, which results in quasiconstant volume combustion and pressure spikes.

Conclusions
Combustion instabilities in a small MMH/NTO LRE are investigated numerically in this paper. A three-dimensional two-phase reaction turbulent flow is predicted by the Eulerian-Lagrangian method, in which the URANS equations are for the gas-phase flow and the DDM are for the trajectories of droplets. The two-phase interactions are modeled by the mass, momentum, and energy source terms in the gas-phase equations. The droplets are injected with given velocities and SMD at the determined positions. The 30 μm, 60 μm, 80 μm, 100 μm, and 140 μm cases are calculated.
As the droplet size increases from 30 microns to 140 microns, peak-to-peak pressure oscillations increase firstly and then decrease. Pressure oscillations are the strongest in 80-micron case. The amplitude and FFT analysis of pressure oscillations show that the 1L and 1T mode self-excited combustion instabilities occur in 60-micron and 80-micron cases. This shows when the droplet size is less than a value in a determined range, combustion instability would occur. This trend is the same as the Hewitt Stability Correction. But if the droplet size is too small such as 30 microns, pressure oscillations would decay and combustion keeps stable.
Further analysis shows that abrupt peaks occur in the mass fraction of MMH and NTO frequently at the downstream region just adjacent to impinging points. The oscilla-tions of mass fraction of MMH are in phase with those of NTO, which are organized and periodic. The local combustible mixtures are formed, and violent combustion arises. The chemical reaction is faster than pressure expansion. Thus, quasiconstant volume combustion and abrupt pressure spikes are induced, which are the sources of pressure oscillations. The propagation and reflection of pressure waves in the chamber stimulate combustion instability. The droplet size affects the mass fraction of MMH and NTO by evaporation rate and initial number density. When the droplet size is too small, the fuel and oxidizer gas are consumed as soon as they are transformed from the fluid phase because of fast evaporation. Local high-dense combustible premixture is hard to be formed. When the droplet size is too large, the evaporation rate is small and there are less abrupt peaks in the mass density of fuel and oxidizer. Local bombing phenomenon appears frequently within a definite range of droplet size.
This study shows where the initial disturbances that induced combustion instability are from and how combustion instabilities are triggered. Furthermore, it reveals that propellant droplet sizes have a great influence on combustion instability and provide guidelines on the atomization of propellant.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest. International Journal of Aerospace Engineering