Effects of LOX Particle Diameter on Combustion Characteristics of a Gas-Liquid Pintle Rocket Engine

LOX/GCH4 pintle injector is suitable for variable-thrust liquid rocket engines. In order to provide a reference for the later design and experiments, three-dimensional numerical simulations with the Euler-Lagrange method were performed to study the effect of the initial particle diameter on the combustion characteristics of a LOX/GCH4 pintle rocket engine. Numerical results show that, as the momentum ratio between the radial LOX jet and the axial gas jet is 0.033, the angle between the LOX particle trace and the combustor axial is very small. Due to the large recirculation zones, premixed combustion mainly occurs in the injector wake region. As the initial LOX particle diameter increases, the LOX evaporation rate and the combustion efficiency decrease until the combustion terminates with the initial LOX particle diameter greater than 110 μm. The oscillation amplitude of the combustor pressure increases significantly along with the increase of the initial LOX particle diameter, and the low-frequency unstable combustion occurs when the initial LOX particle diameter exceeds 60 μm. The combustor pressure oscillation at about 40Hz couples with the swinging process of spray and flame, while the unsteady LOX evaporation amplifies the combustor pressure oscillations at 80Hz and its harmonic frequency.


Introduction
Space activities have been going on for more than half a century, and the rocket engines using toxic propellants have achieved good performance and stability in a large number of launch missions. However, toxic propellants always have the disadvantages of toxicology, high cost, and environmental pollution. With the enhancement of social awareness of environmental protection and the continuous development of manned space activities, nontoxicity propellants will inexorably become the mainstream in the future [1]. Many research results indicate that methane has advantage in various aspects, such as specific impulse, cooling capacity, and maintainability. Methane overcomes the weakness of low density of hydrogen and small specific impulse of kerosene but combines the advantages of them [2,3]. So the overall performance of methane is between kerosene and hydrogen. Additionally, the possibility of refining methane on Mars and Titan is conducive to the on-site refueling and reuse of aircraft in interstellar exploration [4]. Thus, LOX/GCH 4 bipropellant has become an important research direction of chemical propulsion.
The pintle injector has a wide range of applications, such as variable thrust adjustment and continuous rotating detonation [5]. Compared to conventional coaxial [6] and swirl [7] injectors, the pintle injector has distinguishing features of simple structure, good variable-condition adaptability, and stable combustion [8]. The abovementioned features have made a series of successful applications of the pintle injector, including the Apollo lunar descent engine (thrust ratio 10 : 1) [9], the main engine of the Chang'e-3 lunar probe mission (thrust ratio 5 : 1) [10,11], and the Merlin 1D engine of the Falcon 9 rocket (thrust ratio 2 : 1) [12]. In the course of decades of development, pintle injectors have successively experienced three development directions of linear adjustment, low temperature, and nontoxicity [13]. Therefore, using pintle injectors in LOX/GCH 4 variable-thrust rocket engines has important research significance and broad application prospects.
The pintle rocket engine contains the pintle structure that extends into the combustor, and its combustor flow field differs greatly from traditional engines [8,14,15]. There are three types of main design parameters: (a) geometric parameters, including D p , h, d, and L s ; (b) dimensionless parameters, including TMR (the ratio of radial jet to axial jet), BF (the ratio of the sum of orifice diameters to the pintle diameter), the skip distance-to-pintle diameter ratio, and the combustor diameter-to-pintle diameter; and (c) injection parameters, including injection speed and angle. Although the pintle injector has many design parameters and mechanisms, there are few publications about the fundamental researches of the injector. Son et al. [16][17][18] mainly studied the effects of TMR and Weber number on the spray cone angle and D SMD of gas-liquid pintle injectors. A spray simulation method and a set of pintle injector design procedures were proposed on basis of test results. Lee et al. [19] put the primary focus on spray tests of a pintle injector for improving spray uniformity in a 400 N class engine. The empirical correlations of spray angle and D SMD were derived, and then, the optimum control range for the annular flow's orifice area was estimated from the perspective of atomization performance. Zhang et al. [20] analyzed the results of three cold tests and found that increasing the Weber number helps to reduce D SMD and make its distribution more uniform. Sakaki et al. [21][22][23][24] carried out a large number of firing tests for LOX/alcohol pintle rocket engines and found that the combustion efficiency decreases as TMR increases. Besides, the coupling between the heat release rate and the combustion chamber acoustics occurs with the chamber acoustic natural frequencies when TMR is smaller than unity. Cheng et al. [7,[25][26][27][28] researched the spray characteristics of pintle injector through cold tests and theoretically derived the prediction formula of its spray cone angle. Fang and Shen [29] adopted the D SMD obtained in cold tests as the initial condition for two-dimensional numerical simulations and studied the effects of various structural parameters on the combustion performance of the LOX/GCH 4 pintle rocket engine thoroughly.
The combustion flow fields are important to understand the combustion characteristics of liquid engines. Asakawa et al. [30] observed spray distribution and evolution of flame front in a LOX/LNG rocket chamber with single coaxial injector element. Lux and Haidn [31] provided visualization of combustion flow field in a LOX/GCH 4 subscale thruster and studied the effect of a recessed LOX tube in shear coaxial injection on the flame roughness and stability. As for the pintle injector, Sakaki et al. [21] achieved the CH * chemiluminescence and backlit spray measurement of a LOX/alcohol pintle rocket engine, but the images are not clear. The related numerical works of pintle engines mostly concerned twodimensional simulations.
Generally, fundamental research about pintle injectors is little, and the published literature mainly focuses on TMR and structural parameters. In the present study, the effect of initial particle diameter on the combustion characteristics of a LOX/GCH 4 pintle rocket engine was studied through three-dimensional numerical simulations. The simulation results provide a reference for the design optimization and firing tests of pintle injectors.

Physical Model and Numerical Methodology
2.1. Configure. The three-dimensional structure of the LOX/GCH 4 pintle rocket engine prototype and its computational domain are shown in Figure 1. The engine comprises a gas-liquid pintle injector, a rectangular combustion chamber with an observation window, and a short contractionexpansion nozzle to choke the flow. The origin of the computational domain is the cross-section center of the combustion chamber inlet. The length, height, and width of the rectangular combustion chamber are 348 mm in the x-direction, 141 mm in the y-direction, and 45 mm in the z-direction, respectively. The heights at the nozzle throat and outlet are 8.2 mm and 11.9 mm, respectively. The corresponding subsonic angle of the combustion chamber is 31°. Multiple monitoring points are set in the computational domain to obtain the local pressure and temperature during unsteady numerical simulations, where P01 and P02 are located in the gas manifold and P1 to P8 in the combustion chamber. The spatial coordinates of these monitoring points are shown in Table 1. The designed mass flow rate of the engine is 598 g/s (112 g/s for LOX, 43 g/s for GCH 4 , and 443 g/s for N 2 , respectively), and the fuel is the premixed gas, methane and nitrogen. This helps to maintain the targeted combustor pressure (1.8 MPa) to facilitate ignition and reduce the combustor temperature with no additional thermal protection required.
As shown in Figure 2, LOX is radially injected from 12 injection orifices at the pintle tip after passing through the pintle central flow channel. The premixed fuel gas in the gas manifold flows through the axial injection channel and impinges on the radial LOX jets with an impact angle of 90°. Then, the atomization and combustion occur in the combustor. The structural parameters shown in Figure 2 are listed in Table 2. For the pintle injector, the designed pressure drops for LOX and premixed fuel gas are both 1 MPa.

Numerical
Method. Three-dimensional numerical simulations on the LOX/GCH 4 pintle rocket engine were conducted utilizing the Euler-Lagrange method in the commercial code ANSYS. The fluid phase is treated as a continuum by solving the Navier-Stokes equations, while the dispersed phase is solved by tracking a large number of particles through the calculated flow field. Flow variables in the computational domain were calculated through solving the discretized governing equations with the density-based implicit algorithm. A second-order, double-precision solver was adopted for spatial and temporal discretization with the finite volume method. The standard k-ε turbulence model was applied for modeling the turbulence, and the standard wall functions are used for the near-wall treatment. The published literature show that the numerical method and its applicability for combustion flow field simulation of liquid rocket engines have been demonstrated by other researchers [6,29].
In the Euler frame, the density-weighted filtered governing equations for continuous phase can be written as [32][33][34] ∂ ρ ∂t where the over bars (-) represent the time-filtered values and the over tildes (~) and quotes (") represent the densityweighted filtered values and their pulsation values.
In the Lagrange frame, the primary and secondary breakup processes of liquid jets are not taken into consider-ation. For simplicity, the atomization process of liquid jets is modeled as droplet particles by means of the DPM. LOX particles with a given diameter and velocity vector are injected into the computational domain from injection orifice surfaces, and no interaction between particles is yet considered. The governing equations for discrete phase can be expressed as particle trajectory equations, and unsteady Lagrange particle tracking method was used to describe the motion and evaporation of each particle. The discrete random walk model could predict the turbulent dispersion of particles by integrating the effect of instantaneous turbulent velocity fluctuations on the particle trajectories. The spherical drag law was adopted without considering the nonsphericity of the droplet shape. The interaction of the droplet with combustor walls utilized the reflect-type boundary condition with all normal and tangential momentum retained. During the calculation process, required information of the surrounding gas is provided by the governing equations for continuous phase, while the gas source terms are calculated and delivered to the latter. For the details of DPM, the reader should consult Ref. [35].
In the present simulation, the EDC was applied for turbulent combustion simulation, and the reaction rate of every reaction is calculated individually by formulating smallscale eddies. The eddy scale is governed by the turbulent kinetic energy and the turbulent dissipation rate in turbulence flow. Giorgi et al. [36] compared multiple turbulent combustion models for the coaxial rocket engine and revealed that the EDC was numerically cheaper to predict the combustor temperature and flame length. Besides, the  International Journal of Aerospace Engineering JL mechanism is used to model real chemical reactions, as listed in Table 3. Andersen et al. [37] combined the JL with the EDC and provided an accurate temperature field and CO distribution under the oxygen-methane condition.

Validation of the Euler-Lagrange Method.
Yuan [6] adopted the Euler-Lagrange method to simulate the combustion flow field of an air heater. The LOX spray length and cone angle obtained with numerical simulations are identical to the high-speed photography results. In this section, the optical measurement data from Ref. [38] was employed to further examine the validity of the Euler-Lagrange method in the study of spray, evaporation, and mixing. The corresponding computational domain is the rectangular duct, as shown in Figure 3. The red zone represents the inlet of heated air with a velocity of 120 m/s and a temperature of 750 K. The kerosene injection surface is located in the center of the red zone, from which liquid-kerosene jet (1 g/s) is injected into the rectangular duct at 453 K. The injected kerosene has a streamwise velocity component of 10 m/s, while its transverse velocity component is neglected. The outlet pressure and the operation pressure are set to 0.9 MPa and 0, respectively. The dimensions of the rectangular duct are 25 mm × 40 mm × 150 mm, and those of the kerosene injection surface are 5 mm × 1 mm. The pressure within the rectangular duct is constant at about 0.9 MPa. According to the precise spray boundary conditions provided by Ref. [38], the Rosin-Rammler distribution with the D SMD of 14.6 μm is applied to the initial kerosene particle diameter.
The spray trajectories of kerosene particles in the rectangular duct are shown in Figure 4, and these small spheres are colored by particle mass. Kerosene particles gradually evaporate as they flow downstream, and the evaporation products diffuse into heated air. This results in a decrease in the streamwise velocity component of the surrounding gas around the center of the rectangular duct, as shown in the velocity distribution of Figure 4. Because the transverse velocity component of kerosene particles at the injection surface is neglected in the numerical simulation, the dispersion ranges of kerosene particles in the x-and y-direction are slightly narrower than the experimental results shown in Ref. [38]. Figure 5 presents the quantitative comparisons of the numerical and experimental results [38], which include the D SMD , D 0:1 , and D 0:9 , as well as a cross-sectional averaged degree of evaporation. The degree of evaporation is defined as the ratio of the evaporated to the initial kerosene mass flow rate at one certain plane. Generally, the numerical results agree well with Ref. [38], while the difference is mainly due to the uncertainties of boundary conditions and kerosene physical properties. Thus, the Euler-Lagrange method can be applied to simulate the evaporation characteristics of liquid jets.

Boundary Conditions and Computational Mesh.
For the computational domain in Figure 1(b), the characteristic inflow condition with constant temperature (300 K), mass flow rate, and species mass fraction was utilized at the inlet of premixed fuel gas, methane and nitrogen. The LOX jets were assumed as fully atomized spray composed by spherical droplet particles at 80 K. The injection velocity of 33.48 m/s is perpendicular to the outlet surfaces of 12 LOX injection orifices. The initial particle diameter uses uniform distribution with D SMD . The supersonic outflow condition was employed at the nozzle outlet, and the outlet pressure was extrapolated from the solution. The adiabatic, nonslip wall condition was applied for the numerical studies, and the walls will reflect the LOX particles after collision. Figure 6 displays the continuous joint mesh of the LOX/GCH 4 pintle rocket engine. Unstructured mesh with tetrahedral cells was used for meshing the combustor and nozzle, while structured mesh with cuboidal cells for the International Journal of Aerospace Engineering gas manifold and axial injection channel. In order to better resolve the flames and shear layers, the corresponding grids were clustered in the surrounding region of the injector wake and combustor walls.

Results and Discussion
To study the combustion characteristics of the LOX/GCH 4 pintle rocket engine, cases with different initial particle diameters were considered, while the computational domain (as shown in Figure 1) and other numerical simulation setups remain in the same. The initial particle diameter of 10 μm was chosen as the basic case.
3.1. Basic Case. In this section, four sets of mesh were used for the steady simulation of the basic case. Mesh A, B, C, and D have a total of 503665, 740924, 988070, and 1260312 cells, respectively. As shown in Figure 7, the recirculation zones and spray trajectories obtained with 4 sets of mesh are similar    International Journal of Aerospace Engineering to each other. However, it can be seen from Figure 8 that the difference of local x-velocity distribution for mesh B, C, and D is small, but not for mesh A. This illustrates that mesh B, C, and D can properly capture the expansion feature of the premixed fuel gas in the downstream of the pintle tip.
Besides, Figure 8 shows that the cross-sectional area of gas-flow channel is reduced before the premixed fuel gas enters the combustor, while the axial velocity of subsonic gas flow increases gradually. Under the influence of the large recirculation zones between combustor walls and combustor axis in Figure 7, the premixed fuel gas continues to flow along the pintle wall with the same axial velocity after leaving the axial injection channel. The hemispherical shape of the pintle tip results into a local high-pressure zone with small axial velocity near the pintle tip, which is significantly different from that of the previous two-dimensional simulations [29,39]. The large recirculation zones and local high-pressure zone cause throttling area reduction when the premixed fuel gas approaches the pintle tip. Correspondingly, the axial velocity of the premixed fuel gas increases to the sound speed first and then decreases. As the TMR is small (0.033), LOX particles are convected downstream by the premixed fuel gas after being injected radially from the injection orifices. The angle between the LOX particle trace and the combustor axial is very small, and the particles are dispersed in a conical shape. The LOX particles move streamwise and are completely evaporated at x = 65 mm.
The temperature fields for different grid resolutions in Figure 7 are also similar. The high-temperature reaction products are transferred to the head of the combustor due   International Journal of Aerospace Engineering to the large recirculation zones, which contributes to the uniform temperature distribution in the combustor and the lowtemperature zone in the vicinity of pintle tip and combustor axis. Premixed combustion occurs at the border of the highand low-temperature zones. In addition, the asymmetry of the combustion flow field with respect to the y = 0 plane is mainly caused by the Reynolds number and the expansion ratio of the axial injection channel [40,41]. Under the numerical simulation setups, the Reynolds number of the premixed fuel gas in the axial injection channel is 7:93 × 10 6 , and the expansion ratio (the ratio of the height of the combustion chamber to the outer diameter of the axial injection channel) is 8.56. As the Reynolds number and expansion ratio are both large, the asymmetry of the combustion flow field appears with a certain randomness, but the main structure of the combustion flow field is similar. The quantitative comparison results for mesh independence analysis are shown in Figure 9. It can be seen that, compared with the results of mesh D, the relative deviations of mesh B and C are negligible, while those of mesh A is large.
In this study, the combustion efficiency is estimated with η c * , which is evaluated with Equations (2) and (3). c th * is calculated with NASA Chemical Equilibrium with Applications (NASA CEA) [42]. q m is the designed mass flow rate (0.598 kg/s) of the engine, and the A t is the geometric area (45 mm * 8:22 mm) of the rectangular-shaped throat.
As listed in Table 4, the simulation results for different mesh resolutions were compared with the designed values obtained with NASA CEA. Compared with the designed values, all the relative deviations of T c , p c , and c * for different   7 International Journal of Aerospace Engineering grid resolutions are less than 2%. It should be noted that p c and T c refer to the area-weighted average value of the entire combustor (0 ≤ x ≤ 348 mm). And all the combustion efficiencies are higher than 98%, which demonstrates the validity of the engine configuration and numerical simulation setup. However, the η c * for mesh A is close to 100%, and the streamwise length of the lowtemperature region is shorter than that of the others (as seen in Figure 7), indicating that the combustion reaction for mesh A is more violent. To the viewpoint of computational time, mesh B is suitable and sufficient for the subsequent calculation simulations.

Evaporation
Characteristics. The evaporation characteristics of the pintle injector are mainly evaluated by the spray trajectories and particle size distribution. Keeping other setups unchanged, 12 steady simulation cases were conducted with different initial LOX particle diameters (10~120 μm, changed per 10 μm). Spray trajectories for different initial LOX particle diameters are displayed in Figure 10, and the sphere color represents the particle diameter. It can be seen that the dispersion characteristics of LOX particles are sensitive to the initial particle diameter. As the initial particle size increases, the streamwise length of spray trajectories gradually increases until some particles leave the engine without evaporating completely. The transverse diffusion degree of LOX particles along the flow direction also alters significantly for different initial particle diameters. The spray trajectories contract to the combustor axis with the initial LOX particle diameter of 10 μm. The larger the initial particle diameter is, the more obvious the radial expansion trend of the spray trajectories becomes. When LOX particles enter the combustor with the initial particle diameter of 120 μm, more LOX particles have not completely evaporated in the combustor, and the combustion terminates as there is not enough oxygen to participate in chemical reaction.
As illustrated in Figure 11, under the same mass flow rate condition, the evaporation rate of the LOX particles with a smaller initial particle diameter is faster than those with larger initial particle diameter, and so does the decrease rate   Figure 12: Combustor pressure for different initial LOX particle diameters. 8 International Journal of Aerospace Engineering of D SMD . This is because the total surface area of the droplets is inversely proportional to the initial particle diameter, and the heat exchange from LOX particle surfaces with the surrounding gas is stronger with smaller initial particle diameter. The temperature of LOX particles gradually increases as the evaporation proceeds. And the decrease of LOX particle diameter leads to a reduced heat transfer area. Thus, the evaporation rate of LOX particles decreases along with the combustor axis. The initial particle diameter is an important factor that affects the performance of the LOX/CH 4 pintle engine. Figure 12 shows the change curve of the combustor pressure with respect to the initial LOX particle diameter. The increase of the initial LOX particle diameter leads to the decrease of the evaporation rate, and less time is left for oxygen to react with methane. The inadequate chemical reactions will lead to reduced combustion efficiency, which is proportional to the combustor pressure. When the initial LOX particle diameter is larger than 90 μm, the incomplete evaporation of some particles causes the further decline of the amount of oxygen required for the combustion reaction, and the combustion efficiency is accelerated to decrease until the combustion terminates. The initial LOX particle diameter of 110 μm is the limit for maintaining combustion under the fuel-rich condition (the residual oxygen coefficient is 0.65). So, it is necessary to improve the atomizing performance of the pintle injector.

Combustion Stability.
In this section, unsteady simulations with different initial LOX particle diameters, 10 μm, 60 μm, and 110 μm, were conducted to further investigate the effects of the latter on the combustion stability of the LOX/CH 4 pintle engine. The starting moments of unsteady simulations were set as 0 s of combustion process. Combus-tion stability was mainly evaluated in terms of the P ′ in the combustor during the steady stage, and the threshold of 3% is used to classify the stable and unstable combustion [43]. The locations of the monitoring points associated with the oscillation traces mentioned below are listed in Figure 1 and Table 1.
The time step Δt has a great effect on the simulation results and should meet the following condition generally: where Δt is the time step, u is the characteristic flow velocity, and Δx is the local grid size. The injection velocity of the premixed fuel gas is 236.75 m/s, and the minimum grid size near the pintle tip for mesh B is 0.47 mm. Thus, Δt should be smaller than 1:99 × 10 −6 s. In order to reduce the numerical costs of unsteady simulation, unsteady simulations with two different time steps, 1 × 10 −6 s and 5 × 10 −6 s, were used to demonstrate the time step convergence. It can be seen from Figure 13 that, under the initial LOX particle diameter of 60 μm, the time histories of the temperature oscillations at the monitoring point, P4, for the two different time steps are basically the same. So, 5 × 10 −6 s was selected as the time step for the following unsteady simulations. For different initial LOX particle diameters, Figure 14 shows the time histories of the pressure oscillations at P3 and P01 and the corresponding P ′ during the steady stage. As illustrated in Figure 14, the combustion stability is closely related to the initial LOX particle diameter. As the initial LOX particle diameter increases, the time-averaged value of the combustor pressure decreases, but the corresponding P′ increases. The P ′ of the combustor pressure is always lower than 5% due to the intrinsic combustion stability of the pintle injector. The unstable combustion (the P ′ of the combustor pressure > 3%) only occurs when the initial LOX particle diameter is greater than 60 μm. The pressure oscillation in the gas manifold is almost synchronous to that in the combustor, and the former is only slightly lagging behind the latter. However, the P′ of the injection pressure in the gas manifold is smaller than 0.5% and has a nonlinear relationship with the initial LOX particle diameter. Thus, no coupling occurs between the combustor pressure and gas injection pressure.
As displayed in Figure 15, the spectrogram plots of the combustor pressure fluctuations for different LOX particle diameters are similar, and strong peaks are observed in the vicinity of several main frequencies, 40 Hz, 80 Hz, and 160 Hz. As the initial LOX particle diameter increases, the main frequencies decrease slightly, while the corresponding amplitudes increase significantly.
From the evolution of the combustion flow field for different initial LOX particle diameters in one period of pressure oscillation in Figure 16, it can be seen that the spray and flame swing in the y-direction for different initial LOX particle diameters. The evaporation rate decreases along with the increase of the LOX particle diameter, and the reaction 9 International Journal of Aerospace Engineering region spreads downstream of the combustor with LOX spray. This exacerbates the degree of temperature field inhomogeneity in the combustor and the swinging scope of the spray and flame. Figure 17 shows the time histories of the LOX evaporation rate oscillations and their spectrogram plots during the steady stage for different initial particle diameters. The streamwise length of the spray trajectory is short when the initial LOX particle diameter is 10 μm, and the evaporation rate does not show regular oscillations. When the initial LOX particle diameter is 60 μm, the streamwise length of spray trajectory extends significantly, and the evaporation rate has peaked at 80 Hz and its harmonic frequency (for a detailed explanation, see Ref. [44]). The LOX evaporation rate oscillates almost in phase with the combustor pressure. According to the Rayleigh criterion, the unsteady evaporation process of LOX droplets maintains combustor pressure oscillation (80 Hz and its harmonic frequency).
The oscillation amplitude of the evaporation rate increases dramatically with the initial LOX particle diameter increasing to 110 μm, while the oscillation frequency decreases slightly compared to the initial LOX particle diameter of 60 μm. The phase difference between the oscillation waveforms of the LOX evaporation rate and the combustor pressure significantly increases, and the occurrence moments of the local high-temperature region ( Figure 16) and the combustor pressure peak (Figure 14(a)) indicate that the heat release rate is almost antiphase to the pressure oscillation. This means that the initial LOX particle diameter of 110 μm is the upper limit for maintaining low-frequency unstable combustion under such condition, as mentioned in Section 3.2, and a further increase will make the combustion terminate. In addition, the time-averaged evaporation rate becomes smaller as some particles leave the engine without complete evaporation.
Tamura et al. [45] conducted experimental research of the low-frequency combustion instability of the LOX/hydrogen rocket engine and found that there is coupling between the streamwise length of LOX spray trajectories and the combustor pressure. But they did not give the relevant flame behaviors. The present numerical study found that, for the LOX/GCH 4 pintle rocket engine with a small TMR of 0.033, the low-frequency unstable combustion occurs when the initial LOX particle diameter is larger than 60 μm. And the combustor flow field exhibits other typical characteristics: the streamwise length of spray trajectories is almost unchanged, while the swinging process of spray trajectories in the transverse direction is synchronous with the  Figure 15: Spectrogram plots of the pressure trace at P3 for different initial LOX particle diameters shown in Figure 14(a). 10 International Journal of Aerospace Engineering combustor pressure oscillation. Besides, with the appearance and disappearance of the local high-temperature region in the combustion chamber, the combustion reaction also enhances or weakens synchronously.
In order to perform a quantitative analysis of the combustor temperature field, the orthogonal decomposition method summarized briefly in the appendix was employed to separate the time and space information for the evolution  Figure 18, mode-1 contributes far more energy than other modes for different initial LOX particle diameters, which indicates that mode-1 dominates the evolution of the temperature distribution in the z = 0 plane. As shown in Figure 19(a), the spatial structures of mode-1 for different initial LOX particle diameters are symmetrical about the combustor axis, and the signs on both sides are opposite. The scope with larger absolute values is basically consistent with the reaction region in Figure 16, which expands along with the increase of the initial LOX particle diameter. Thus, mode-1 corresponds to the swinging process of the spray and flame in the y-direction. The Fourier transform was carried out for the time coefficient of mode-1, and strong peaks at 43 Hz are observed for different initial LOX particle diameters in Figure 19(b). With reference to Figure 15, it can be seen that the combustor pressure oscillation at 40 Hz is coupled with the combustion process. Under the effect of pressure difference between the upper and lower walls, the recirculation zone in the relatively high-pressure side continues to increase in size until spanning the entire combustor, while the recirculation zone in the other side is gradually squeezed to the head of the combustor. The spray and flame are also pushed to the relatively low-pressure side by the pressure difference. The enhanced combustion reaction in the relatively low-pressure side causes the increase of the local temperature and pressure, and the reverse change occurs on the other side. Accordingly, the spray and the flame move in the reverse direction with the size change of the recirculation zones near the upper and lower walls. It can be drawn that the reciprocating occurrence of the above process maintains the combustor pressure oscillation at 40 Hz.

Conclusions
Numerical simulation was performed to study the effect of the initial LOX particle diameter on the combustion   Figure 18: Energy contribution of the top five modes for different initial LOX particle diameters. 12 International Journal of Aerospace Engineering characteristics of a LOX/GCH 4 pintle rocket engine with a rectangular combustion chamber. The fuel is the premixed gas, methane and nitrogen. Numerical results agree well with the designed values, which indicates that the Euler-Lagrange method used in this paper is reasonable. Several conclusions drawn from the numerical analysis are gained: (1) As the momentum ratio between the radial LOX jet and the axial gas jet is 0.033, LOX particles with the initial particle diameter of 10 μm cannot penetrate the premixed fuel gas after being injected radially into the combustor and are completely evaporated at x = 65 mm. Large recirculation zones are formed in the transverse direction and contribute to the uniform temperature distribution in the combustor. And premixed combustion occurs in the surrounding region of the injector wake. Thus, the combustor flow field differs greatly from the previous research results (2) As the initial LOX particle diameter increases, the evaporation rate of the LOX particles decreases significantly, while the streamwise length of spray trajectories increases until some droplets leave the engine without evaporating completely. Accordingly, the combustor pressure stays away from the designed value, and the combustion terminates when the initial LOX particle diameter is greater than 110 μm (3) The oscillation amplitude of the combustor pressure increases significantly along with the increase of the initial LOX particle diameter, and the lowfrequency unstable combustion occurs when the initial LOX particle diameter exceeds 60 μm. And the main frequencies of the combustor pressure oscilla-tions are mainly concentrated in the vicinity of 40 Hz, 80 Hz, and 160 Hz. The combustor pressure oscillation at about 40 Hz couples with the swinging process of spray and flame in the transverse direction. The combustion reaction is synchronously enhanced or weakened in the swinging process with a large initial LOX particle diameter, while the streamwise length of spray trajectories remains nearly constant. With the same peaks at about 80 Hz and its harmonic frequency during the steady stage, it can be concluded that the unsteady LOX evaporation amplifies the combustor pressure oscillations The combustion efficiency and combustion stability of the LOX/GCH 4 pintle rocket engine can be dramatically improved by decreasing the initial LOX particle diameter. However, the initial LOX particle diameter is closely related to many factors. Subsequent experiments are required to verify the simulation results and further investigate the effects of related parameters, such as the injection orifice dimensions and injection pressure drop, on the combustion characteristics of the LOX/GCH 4 pintle rocket engine.
where Φ = ½Φ 1 ⋯ Φ n ⋯ Φ N is a set of maximum independent vectors in linear space R, i.e., a set of basis in R space;  Figure 19: Space structures (a) of mode-1 and the spectrograms (b) of the corresponding time coefficient for different initial LOX particle diameters.
13 International Journal of Aerospace Engineering Φ n is a mode of UðR, tÞ which includes spatial information; F = ½ f 1 ⋯ f n ⋯ f N is the projection of UðR, tÞ on the vector Φ, i.e., a linear expression relative to the set of basis in R space; f n ðtÞ is the time coefficient of Φ n and represents the projection of UðR, tÞ on the vector Φ n at time t.
It can be assumed that a set of time-discrete data U obtained in the unsteady simulations contains the data, U 1 ðRÞ, U 2 ðRÞ ⋯ U K ðRÞ, which corresponds to discrete time points, t 1 , t 2 ⋯ t K . The time-discrete data U can be decomposed as The error E between the original data U and the reconstructed information U * is Based on the idea of orthogonal decomposition, no two vectors in a set of basis in R space are linearly independent. In order to obtain the set of basis that minimizes the error E, the solution is as follows: (1) Solve the eigenvalue problem of the self-covariance matrix C = U T U: (2) Export the mode Φ k corresponding to the eigenvalue λ k through Equation (A.5): (3) Calculate the time coefficients of all modes using the following equation: (4) The magnitude of the eigenvalue λ k reflects the energy contribution of the corresponding mode Φ k . All the modes ½Φ 1 ⋯ Φ n ⋯ Φ K and their time coefficient ½ f 1 ⋯ f n ⋯ f K are rearranged in the order of eigenvalues from large to small λ 1 > λ 2 > ⋯>λ k Nomenclature A t : Throat area BF: Blocking factor c th * : Theoretical characteristic velocity c * : Characteristic velocity D 0:1 : 10% volume undersize diameters D 0:9 : 90% volume undersize diameters D m : Diffusion coefficient of species m DPM: Discrete phase model D p : Pintle diameter D SMD : Sauter average diameter d: Injection orifice diameter E: Total energy EDC: Eddy dissipation concept model F s,i : Interphase exchange terms h: Axial injection channel width JL: Jones-Lindstedt 4-step mechanism L s : Skip distance N: Injection orifice number P ′ : Pressure oscillation amplitude p: Pressure p c : Combustor pressure Q s : Interphase exchange terms q i : Heat flux in i-direction q m : Mass flow rate Sc t : Schmitt number S h : Source terms due to chemical reaction S s,m : Interphase exchange terms T c : Combustor temperature TMR: Total momentum ratio t: Physical time Δt: Time step u: Characteristic flow velocity u: Velocity in i-direction Δx: Local grid size Y m : Mass fraction of species m η c * : characteristic exhaust velocity efficiency μ t : Turbulent viscosity coefficient ρ: Density ρ s : Interphase exchange terms τ ij : Viscous stress ω m : Source terms due to chemical reaction.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
The authors declare no conflict of interest.