Research Article  Open Access
Jianxiu Qin, Huiqiang Zhang, "Numerical Analysis of SelfExcited Combustion Instabilities in a Small MMH/NTO Liquid Rocket Engine", International Journal of Aerospace Engineering, vol. 2020, Article ID 3493214, 17 pages, 2020. https://doi.org/10.1155/2020/3493214
Numerical Analysis of SelfExcited Combustion Instabilities in a Small MMH/NTO Liquid Rocket Engine
Abstract
Combustion instabilities in a small MMH/NTO liquid rocket engine used for satellite attitude and course control are numerically investigated. A threedimensional NavierStokes code is developed to simulate twophase 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 selfexcited combustion instabilities are captured in the 60micron and 80micron 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 highdense 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 highdense premixtures and combustion is stable in the chamber.
1. 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. Highfrequency 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–6]. But the quantity of tests is needed due to the notoriously probability and uncertainty of combustion instability. Thus, the prediction of combustion 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–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–21]. Detached Eddy Simulation (DES) featuring little demands on computational resources has been developed to investigate combustion instability [22–24]. However, liquid droplet dynamic was not considered and the propellants were gasified when they left the injectors in these simulations. For twophase turbulent reactive flow, Unsteady ReynoldsAveraged NavierStokes (URANS) is always employed [25, 26] and has been applied to study combustion instability in LRE [27–29]. Tucker et al. [30, 31] compared LES and URANS simulations of a combustor. The results showed the highestfidelity LES and URANS best predicted the instantaneous temperature, vorticity fields, and wall heat flux, while the two lowerorder 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 twophase interactions.
As is well known, highfrequency 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. Highfrequency combustion instabilities can be classified as intrinsic and injectioncoupled instabilities [2]. Injectioncoupled 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–38]. For the cases in this work, likedoublet 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 , where is the injector’s orifice diameter and is the injection velocity of the least volatile propellant. A reduction in , 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 threedimensional twophase 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.
2. Numerical Model
In a small MMH/NTO LRE, the propellants undergo injection, atomization, vaporization, mixing, chemical reaction, and expansion process. The EulerianLagrangian method is employed to solve twophase reactive flow in the combustion chamber, in which the NavierStokes (NS) 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.
2.1. Continuous Phase Model
Three dimensional, multispecies gasphase NavierStokes equations are used to describe the flow in the chamber. The standard twoequation 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 , , , and represent conservation variables, diffusion coefficients, convective source terms, and the velocity in three directions. and are the source terms determined by droplet vaporization and chemical reaction. The detailed variables are presented in Reference [44].
2.2. 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 NT distribution [45] at given locations, where are the locations of impinging points of the likedoublet 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.
2.3. Droplet Phase Model
The Discrete Droplet Model (DMD) is used to describe the trajectories of droplets with the assumption that the droplets are enough diluted so that the interactions of droplets are ignored. The interaction of droplets and the surrounding gas are included in the source terms in gasphase 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 represents the droplet. The subscript denotes droplet species. is the drag coefficient of the droplets. is the mean gas velocity; is the velocity of the droplet. is the turbulence velocity fluctuation.
The source terms for the twophase interactions are described in detail in Reference [44].
2.4. Chemical Reaction Model
To make this problem simple, only the gasphase chemical reaction is considered, while the liquidphase reaction is neglected. The kinetically chemical reaction rate is calculated by the Arrhenius expression of fourstep chemical reactions. The MMHNTO fourstep chemical reactions are expressed:
3. 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 secondorder upwind scheme, respectively. First, the source terms due to evaporation and twophase interaction are calculated in a Lagrangian phase. Then, the field coupling pressure with velocity is solved by the SemiImplicit 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 socalled Arbitrary LagrangianEulerian (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.
4. 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 likedoublet injectors are in the pattern of fueloxidizerfueloxidizer. The atomization process is neglected, and the droplets are injected in the impinging points shown as Figure 1(b). The likedoublet injectors (impinging points) are distributed in two rings. The first ring is located at mm, while the second ring is located at mm. 20 likedoublet injectors are uniformly distributed in the first ring every 18 degrees in the circumference direction, while 40 likedoublet injectors are uniformly distributed in the second ring every 9 degrees. The cooling injectors are located at 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.
(a)
(b)
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 twothirds 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 NASACEA.

4.2. The Hewitt Criterion
The stability correlating parameter 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 . is the injector’s orifice diameter, while is the velocity of the least volatile propellant. A reduction in indicates a decreased stability margin. It can be concluded that when is less than 0.1, combustion instability may occur.
(a)
(b)
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 are controlled by adjusting velocities of fuel (oxidizer). Thus, different 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 . It can be found that the Hewitt criterion is applicable to this combustor. When the parameter is less than a determined value, combustion instability occurs.
(a)
(b)
Anderson et al. [42, 43] declared that an increase in implied an increase in the mean droplet size. In the atomization test for likedoublet 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 is the density of the droplet. The SMD is proportional to . 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 is less than 15. It can be estimated that when equals to 15 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.

4.3. 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 1Case 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 1Case 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 onestep 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 highpressure and hightemperature condition in the hot test. In this paper, it focuses on whether combustion instability is acoustic combustion instability (highfrequency 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.
(a)
(b)
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 3Case 5. The pressure firstly declines to a lower value and then rises to a higher value from Case 3 to Case 5. Because the time lag becomes larger with the increase of the droplet size, it is important to note that pressure spikes due to the autoignition behavior may produce enough disturbances to trigger combustion instability. However, it is not a reason to induce combustion instability in this paper, because pressure oscillations decay and keep stable in Case 1 and Case 4, while highamplitude pressure oscillations are excited in Case 2 and Case 3.
In Case 1, pressure oscillations decay within the amplitude of 10% of mean chamber pressure. The combustion is considered stable. For Case 2, the peaktopeak amplitude of pressure oscillations steadies out approximately 0.7 MPa (64.2% of mean chamber pressure). Pressure waves are steepfronted and not symmetric due to the nonlinear effect, which is the second feature of pressure oscillations. For Case 3, the peaktopeak 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 and then rises to 1.1 MPa at 3.0 ms. Mean chamber pressure varies slightly with time. There not exist highamplitude 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 highamplitude 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 RP1 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 50micron droplets were not stable. But the 100micron droplets were stable to a simulated triggering mechanism.
4.4. 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 highpressure regions appear at the head of the chamber shown as Figure 7(a). Pressure waves propagate from a highpressure region to a lowpressure region. In the transverse section, the center of the section becomes a lowpressure region, while the region near the sidewall becomes a highpressure 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 lowpressure region, while the contraction section becomes a highpressure 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 highpressure 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 highpressure 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.
(a)
(b)
(c)
(d)
In order to reveal how the small local highpressure 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 likedoublet 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 peaktopeak 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(b). Pressure oscillations at Point A may result from the propagation of pressure waves produced at Point B. Pressure peaks at Point B may be the sources of pressure oscillations in the chamber.
(a)
(b)
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 ms to ms at Point B. The and are marked by squares in Figure 9. For the short interval μ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 50% of the pressure at , 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 is the largest in the observed control volume B, while the evaporation amount of MMH at is the largest. This means violent combustion occurs from to . 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.126 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 twophase 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 twophase 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 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 constantvolume combustion theory and numerical simulation results is shown in Table 4. For the short interval μ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.933 g/cm^{3} at , while it is 5.143 g/cm^{3} at . 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 . 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. The oscillations of mass fraction of the oxidizer are in phase with those of mass fraction of the fuel. The mass fraction of fuel gas and oxidizer gas increases first from 8.335 ms to 8.338 ms. The enhancement of evaporation of fuel and oxidizer droplets leads to the accumulation of fuel and oxidizer 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 to . Pressure waves cannot propagate to the surrounding control volumes in this short period, which leads to a highamplitude 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 RP1 droplets along the same vector in a given angel to mimic fan atomization may trigger selfexcited instabilities, because the LOX and RP1 droplets can interact more rapidly. Zhang et al. [44] found that quasiconstant volume combustion resulted in bombing, which induced selfexcited combustion instability. Grenda et al. [50] found that the simulated bombing occurred when the droplets were on the order of 100 microns and 200 microns. Daimon et al. [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 likedoublet 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 likedoublet 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 60microns 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 likedoublet 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.
(a)
(b)
(c)
(d)
4.5. 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. Selfexcited 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 ms10 ms in Case 1Case 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.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
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 1Case 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 sections, there are small highpressure 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. Violent chemical reactions happen in these regions. For Case 1, the amplitudes of the red regions (highpressure 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.
4.6. 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 mm, the mass fraction of MMH is greater than that of NTO, while at mm18 mm, the mass fraction of NTO is greater. It is fuel rich for mm, because the cooling injectors are located there and spray MMH. Moreover, the quantity of mass fraction peaks of MMH and NTO occurs at mm and 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. For Case 2 and Case 3, the oscillations of mass fraction of MMH and NTO are organized and periodic especially at mm, 18 mm, and 22 mm. From mm, it can note that the frequencies of the mass fraction of MMH and NTO are both estimated as 5000 Hz by counting the number of peaks, about five, over 1 ms time slice. The oscillations of mass fraction of MMH are in phase with those of mass fraction of NTO. 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 highamplitude 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 selftriggered combustion instability.
4.7. 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 highdense 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 80 microns, indicating an increased tendency towards instability, implies a reduction in the stability parameter in the Hewitt Stability Correlation. In the actual test, the mean droplet size is greater than 30 microns. The 30microns 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 highdense 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 highdense 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.
5. Conclusions
Combustion instabilities in a small MMH/NTO LRE are investigated numerically in this paper. A threedimensional twophase reaction turbulent flow is predicted by the EulerianLagrangian method, in which the URANS equations are for the gasphase flow and the DDM are for the trajectories of droplets. The twophase interactions are modeled by the mass, momentum, and energy source terms in the gasphase 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, peaktopeak pressure oscillations increase firstly and then decrease. Pressure oscillations are the strongest in 80micron case. The amplitude and FFT analysis of pressure oscillations show that the 1L and 1T mode selfexcited combustion instabilities occur in 60micron and 80micron 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 oscillations 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 highdense 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.
Acknowledgments
This work is supported by the National Natural Science Foundation of China under Grant No. 91841303.
References
 D. T. Harrje and F. H. Reardon, Liquid propellant rocket combustion instability, NASA, Washington, DC, 1972.
 V. Yang and W. E. Anderson, Liquid Rocket Engine Combustion Instability, AIAA, Washington, DC, 1995.
 F. Richecoeur, P. Scouflaire, S. Ducruix, and S. Candel, “HighFrequency transverse acoustic coupling in a multipleinjector cryogenic combustor,” Journal of Propulsion and Power, vol. 22, no. 4, pp. 790–799, 2006. View at: Publisher Site  Google Scholar
 K. Miller, J. Sisco, N. Nugent, and W. Anderson, “Combustion instability with a singleelement swirl injector,” Journal of Propulsion and Power, vol. 23, no. 5, pp. 1102–1112, 2007. View at: Publisher Site  Google Scholar
 B. Pomeroy and W. Anderson, “Transverse instability studies in a subscale chamber,” Journal of Propulsion and Power, vol. 32, no. 4, pp. 939–947, 2016. View at: Publisher Site  Google Scholar
 C. Ruan, F. Chen, W. Cai, Y. Qian, L. Yu, and X. Lu, “Principles of nonintrusive diagnostic techniques and their applications for fundamental studies of combustion instabilities in gas turbine combustors: A brief review,” Aerospace Science and Technology, vol. 84, pp. 585–603, 2019. View at: Publisher Site  Google Scholar
 A. P. Dowling and S. R. Stow, “Acoustic analysis of gas turbine combustors,” Journal of Propulsion and Power, vol. 19, no. 5, pp. 751–764, 2003. View at: Publisher Site  Google Scholar
 J. Portillo, J. Sisco, Y. Yu, V. Sankaran, and W. Anderson, “Application of a generalized instability model to a longitudinal mode combustion instability,” in 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, Cincinnati, OH, 2007, July. View at: Publisher Site  Google Scholar
 W. A. Sirignano and P. P. Popov, “Twodimensional model for liquidrocket transverse combustion instability,” AIAA Journal, vol. 51, no. 12, pp. 2919–2934, 2013. View at: Publisher Site  Google Scholar
 M. Sliphorst, S. Gröning, and M. Oschwald, “Theoretical and experimental identification of acoustic spinning mode in a cylindrical combustor,” Journal of Propulsion and Power, vol. 27, no. 1, pp. 182–189, 2011. View at: Publisher Site  Google Scholar
 S. K. Kim, D. Kim, and D. J. Cha, “Finite element analysis of selfexcited instabilities in a lean premixed gas turbine combustor,” International Journal of Heat and Mass Transfer, vol. 120, pp. 350–360, 2018. View at: Publisher Site  Google Scholar
 Y. C. Yu, J. C. Sisco, V. Sankaran, and W. E. Anderson, “Effects of mean flow, entropy waves, and boundary conditions on longitudinal combustion instability,” Combustion Science and Technology, vol. 182, no. 7, pp. 739–776, 2010. View at: Publisher Site  Google Scholar
 R. Smith, M. Ellis, G. Xia, V. Sankaran, W. Anderson, and C. L. Merkle, “Computational investigation of acoustics and instabilities in a longitudinalmode rocket combustor,” AIAA Journal, vol. 46, no. 11, pp. 2659–2673, 2008. View at: Publisher Site  Google Scholar
 X. Han, J. Li, and A. S. Morgans, “Prediction of combustion instability limit cycle oscillations by combining flame describing function simulations with a thermoacoustic network model,” Combustion and Flame, vol. 162, no. 10, pp. 3632–3647, 2015. View at: Publisher Site  Google Scholar
 J. Li, Y. Xia, A. S. Morgans, and X. Han, “Numerical prediction of combustion instability limit cycle oscillations for a combustor with a long flame,” Combustion and Flame, vol. 185, pp. 28–43, 2017. View at: Publisher Site  Google Scholar
 R. Garby, L. Selle, and T. Poinsot, “Simulation aux grandes echelles des instabilites de combustion dans un bruleur de longueur variable,” Comptes Rendus Mécanique, vol. 341, no. 12, pp. 220–229, 2013. View at: Publisher Site  Google Scholar
 S. Matsuyama, J. Shinjo, S. Ogawa, and Y. Mizobuchi, “LES of highfrequency combustion instability in a single element rocket combustor,” in 50th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, p. 1271, Nashville, Tennessee, January 2012. View at: Publisher Site  Google Scholar
 S. Matsuyama, J. Shinjo, and Y. Mizobuchi, “LES of highfrequency combustion instability in a rocket combustor,” in 51st AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, p. 0564, Grapevine (Dallas/Ft. Worth Region), Texas, January 2013. View at: Publisher Site  Google Scholar
 S. Srinivasan, R. Ranjan, and S. Menon, “Flame dynamics during combustion instability in a highpressure, shearcoaxial injector combustor,” Flow, Turbulence and Combustion, vol. 94, no. 1, pp. 237–262, 2015. View at: Publisher Site  Google Scholar
 C. Umphrey, M. E. Harvazinski, S. A. Schumaker, and V. Sankaran, “Largeeddy simulation of singleelement gascentered swirlcoaxial injectors for combustion stability prediction,” in 53rd AIAA/SAE/ASEE Joint Propulsion Conference, p. 4698, Jul 2017. View at: Publisher Site  Google Scholar
 C. Kraus, L. Selle, and T. Poinsot, “Coupling heat transfer and large eddy simulation for combustion instability prediction in a swirl burner,” Combustion and Flame, vol. 191, pp. 239–251, 2018. View at: Publisher Site  Google Scholar
 M. E. Harvazinski, W. E. Anderson, and C. L. Merkle, “Analysis of selfexcited combustion instabilities using two and threedimensional simulations,” Journal of Propulsion and Power, vol. 29, no. 2, pp. 396–409, 2013. View at: Publisher Site  Google Scholar
 M. E. Harvazinski, C. Huang, V. Sankaran et al., “Coupling between hydrodynamics, acoustics, and heat release in a selfexcited unstable combustor,” Physics of Fluids, vol. 27, no. 4, pp. 045102–045127, 2015. View at: Publisher Site  Google Scholar
 T. M. Nguyen, P. P. Popov, and W. A. Sirignano, “Longitudinal combustion instability in a rocket engine with a single coaxial injector,” Journal of Propulsion and Power, vol. 34, no. 2, pp. 354–373, 2018. View at: Publisher Site  Google Scholar
 M. J. Nusca, “Utility of computational modeling for the study of combustion instability in small MMHNTO liquid rocket engines,” in 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, p. 5562, Cincinnati, OH, July 2007. View at: Publisher Site  Google Scholar
 L. B. Zhang, M. Chu, and X. Xu, “Performance prediction of apogee attitude and orbit control thruster for MMH/NTO hypergolic bipropellant,” in 50th AIAA/ASME/SAE/ASEE Joint Propulsion Conference, p. 3572, Cleveland, OH, July 2014. View at: Publisher Site  Google Scholar
 M. Habiballah, D. Lourme, and F. Pit, “PHEDRE  Numerical model for combustion stability studies applied tothe Ariane Viking engine,” Journal of Propulsion and Power, vol. 7, no. 3, pp. 322–329, 1991. View at: Publisher Site  Google Scholar
 O. Knab, D. Preclik, and D. Estublier, “Flow field prediction within liquid film cooled combustion chambers of storable bipropellant rocket engines,” in 34th AIAA/ASME/SAE/ASEE Joint Propulsion Conference and Exhibit, p. 3370, Cleveland,OH,U.S.A., July 1998. View at: Publisher Site  Google Scholar
 F. Zhuang, W. Nie, and Q. Zou, “Numerical simulation of MMH/NTO rocket engine combustion instability.,” in 35th Joint Propulsion Conference and Exhibit, p. 2779, Los Angeles,CA,U.S.A., June 1999. View at: Publisher Site  Google Scholar
 P. Tucker, S. Menon, C. Merkle, J. Oefelein, and V. Yang, “An approach to improved credibility of CFD simulations for rocket injector design,” in 43rd AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, p. 5572, Cincinnati, OH, July 2007. View at: Google Scholar
 P. Tucker, S. Menon, C. Merkle, J. Oefelein, and V. Yang, “Validation of highfidelity CFD simulations for rocket injector design,” in 44th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, p. 5226, Hartford,CT, July 2008. View at: Publisher Site  Google Scholar
 R. Smith, G. Xia, W. A. Anderson, and C. L. Merkle, “Computational simulations of the Effect of backstep height on nonpremixed combustion instability,” AIAA Journal, vol. 48, no. 9, pp. 1857–1868, 2010. View at: Publisher Site  Google Scholar
 M. L. Frezzotti, F. Nasuti, C. Huang, C. L. Merkle, and W. E. Anderson, “Quasi1D modeling of heat release for the study of longitudinal combustion instability,” Aerospace Science and Technology, vol. 75, pp. 261–270, 2018. View at: Publisher Site  Google Scholar
 J. Rayleigh, “The explanation of certain acoustical phenomena^{1},” Nature, vol. 18, no. 455, pp. 319–321, 1878. View at: Publisher Site  Google Scholar
 V. G. Bazarov and V. Yang, “Liquidpropellant rocket engine injector dynamics,” Journal of Propulsion and Power, vol. 14, no. 5, pp. 797–806, 1998. View at: Publisher Site  Google Scholar
 S. Gröning, J. S. Hardi, D. Suslov, and M. Oschwald, “Injectordriven combustion instabilities in a hydrogen/oxygen rocket combustor,” Journal of Propulsion and Power, vol. 32, no. 3, pp. 560–573, 2016. View at: Publisher Site  Google Scholar
 A. Urbano, L. Selle, G. Staffelbach et al., “Exploration of combustion instability triggering using large eddy simulation of a multiple injector liquid rocket engine,” Combustion and Flame, vol. 169, pp. 129–140, 2016. View at: Publisher Site  Google Scholar
 W. Armbruster, J. S. Hardi, D. Suslov, and M. Oschwald, “Experimental investigation of selfexcited combustion instabilities with injection coupling in a cryogenic rocket combustor,” Acta Astronautica, vol. 151, pp. 655–667, 2018. View at: Publisher Site  Google Scholar
 A. Duvvur, C. H. Chiang, and W. A. Sirignano, “Oscillatory fuel droplet vaporizationdriving mechanism for combustion instability,” Journal of Propulsion and Power, vol. 12, no. 2, pp. 358–365, 1996. View at: Publisher Site  Google Scholar
 W. A. Sirignano and G. Wu, “Multicomponentliquidfuel vaporization with complex configuration,” International Journal of Heat and Mass Transfer, vol. 51, no. 1920, pp. 4759–4774, 2008. View at: Publisher Site  Google Scholar
 S. Lei and A. Turan, “Chaotic modelling and control of combustion instabilities due to vaporization,” International Journal of Heat and Mass Transfer, vol. 53, no. 21–22, pp. 4482–4494, 2010. View at: Publisher Site  Google Scholar
 W. E. Anderson, H. M. Ryan, and R. J. Santoro, “Combustion instability mechanisms in liquid rocket engines using impinging jet injectors,” in 1st Joint Propulsion Conference and Exhibit, p. 2357, San Diego,CA,U.S.A., July 1995. View at: Publisher Site  Google Scholar
 W. E. Anderson, K. L. Miller, H. M. Ryan, S. Pal, R. J. Santoro, and J. L. Dressler, “Effects of periodic atomization on combustion instability in liquidfueled propulsion systems,” Journal of Propulsion and Power, vol. 14, no. 5, pp. 818–825, 1998. View at: Publisher Site  Google Scholar
 H. Q. Zhang, Y. J. Ga, B. Wang, and X. L. Wang, “Analysis of combustion instability via constant volume combustion in a LOX/RP1 bipropellant liquid rocket engine,” Science China Technological Sciences, vol. 55, no. 4, pp. 1066–1077, 2012. View at: Publisher Site  Google Scholar
 W. S. Nie and S. J. Feng, Chemical Kinetics Model and Numerical Simulation of Liquid Rocket Engine, National Defense Industry Press, Beijing, 2011.
 C. W. Hirt, A. A. Amsden, and J. L. Cook, “An arbitrary LagrangianEulerian computing method for all flow speeds,” Journal of Computational Physics, vol. 135, no. 2, pp. 203–216, 1997. View at: Publisher Site  Google Scholar
 T. S. Wang and Y. S. Chen, “Unified NavierStokes flowfield and performance analysis of liquid rocket engines,” Journal of Propulsion and Power, vol. 9, no. 5, pp. 678–685, 1993. View at: Publisher Site  Google Scholar
 J. A. Keenan, M. Z. Pindera, and M. G. Giridharan, Simulation of spray combustion instability in liquid rocket engines, AIAA Paper, 1997.
 Y. Kim, C. Chen, J. Ziebarth, and Y. CHEN, “Prediction of high frequency combustion instability in liquid propellant rocket engines,” in 28th Joint Propulsion Conference and Exhibit, p. 3763, Nashville,TN,U.S.A., July 1992. View at: Publisher Site  Google Scholar
 J. M. Grenda, S. Venkateswaran, and C. L. Merkle, “Mutidimensional analysis of combustion instabilities in liquid rocket motors,” in 8th Joint Propulsion Conference and Exhibit, p. 3764, Nashville,TN,U.S.A., July 1992. View at: Publisher Site  Google Scholar
 W. Daimon, Y. Gotoh, and I. Kimura, “Mechanism of explosion induced by contact of hypergolic Liquids,” Journal of Propulsion and Power, vol. 7, no. 6, pp. 946–952, 1991. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2020 Jianxiu Qin and Huiqiang Zhang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.