Simulation Study of the Swirl Spray Atomization of a Bipropellant Thruster under Low Temperature Conditions

Abstract

The spray atomization of an injector significantly influences the performance and working life span of a bipropellant thruster of a spacecraft. Deep space exploration requires the thruster to be able to operate reliably at a low temperature range from −40 °C to 0 °C, so the effect of low temperature conditions on the atomization characteristics of injector spray is motivated to be comprehensively investigated. To study the swirl atomization characteristics of MMH (methylhydrazine), which is more difficult to atomize than NTO (nitrogen tetroxide), numerical simulations were conducted, employing the methods of VOF (volume of fluid) and LES (large eddy simulation) under low temperature conditions. The physical model with a nozzle size of 0.5 mm and boundary conditions with a velocity inlet of 3.89 m/s both follow the actual operation of thrusters. The development of spray atomization at low temperatures was observed through parametric comparisons, such as spray velocity, liquid total surface area, droplet particle size distribution, spray cone angle and breakup distance. When the temperature decreased from 20 °C to −40 °C at the same condition of flowrate inlet, those atomization characteristics of MMH propellant vary following these rules: the spray ejection velocity of MMH is significantly reduced by 7.7%, and gas-liquid disturbance sequentially decreases; the liquid film development is more stable, with a negative influence on atomization quality, causing difficulties for primary and secondary breakup, so the total surface area of droplets also decreases by 6.4%; the spatial distribution characteristics, spray cone angle and breakup distance vary less than 5%. Therefore, the low temperature condition can directly lower the combustion efficiency of thrusters with obvious performance degradation, but there are no significant changes in the propellant mixing and liquid film cooling. It is concluded that the bipropellant thruster can reliably work at low temperatures around −40 °C for deep space probe operation.

1. Introduction

Swirl injector spray plays a very decisive role in a bipropellant thruster. After fine atomization and organization, the propellant ejected from the swirl injector can effectively evaporate, mix and then combustion in a downstream chamber to power a spacecraft. The process of propellant atomization and evaporation directly affects the combustion efficiency and the cooling of the combustion chamber wall. Therefore, the atomization quality of the injector spray field significantly influences the performance and working life span of bipropellant thrusters [1,2]. The simple structure of the swirl injector benefits the reliability and stability of thrusters [3]. As a typical and simple design of a mechanical pressure atomizer, the swirl injector has been widely used in spacecraft thrusters [4]. Although the structure of a swirl injector is simple, its spray process is rather complicated, and it is also affected by many factors such as structure parameters, working condition parameters, manufacturing quality and environmental conditions [5].

With the continuous development of deep space exploration technology, the detection of more distant planets, such as Mars, Jupiter and Saturn, has been carried out [6]. With the increase in the exploration distance between the target of detection and the sun, the solar radiation energy density decreases rapidly [7], which causes a huge challenge for the energy acquisition. This problem makes the temperature control of the propellant more and more difficult. The spacecrafts mainly use MMH as fuel and NTO as an oxidizer for propellant combination [8,9], which has the advantages of mature technology, stability and reliability. The freezing point of MMH is −52.2 °C, but the freezing point of NTO is only −11.2 °C. The freezing point of the oxidizer can be reduced to about −54 °C by adding 25% Nitric Oxide into NTO, thereby satisfying the temperature control requirements of the propulsion system under low temperature conditions. However, under these harsh environmental conditions, the viscosity and the surface tension of the liquid propellant rises dramatically. The quality of atomization deteriorates accordingly, and the primary and secondary breakup also become difficult, with a negative influence on the subsequent evaporation and combustion. Especially for the MMH propellant with a larger viscosity and surface tension, this phenomenon becomes more significant.

Many researchers have carried out in-depth research on liquid atomization under low temperature conditions. Ding [10] pointed out that cryogenic fluids show better physical properties than normal temperature fluids for enhanced spray breakup at low temperatures, indicating that fluid properties are important for low-temperature atomization. Zhang [11] employed the vibrating string method to test the kinematic viscosity characteristics of aviation fuel RP-3 under the same pressure condition at a different temperature. He found that the kinematic viscosity of RP-3 gradually increased as the temperature decreased from 0 °C to −30 °C. Both Wei [12] and Rostami [13] conducted experiments to study the effect of temperature on the atomization of aviation fuels. Their results show that the viscosity, molecular tension, breakup length and liquid film thickness increase at low temperatures, the SMD (Sauter Mean Diameter) is higher than that at room temperature and the spray cone angle decreases. Li [14] found that the temperature condition changes the physical properties of fuel, including density, surface tension and viscosity, and then affects the atomization dynamic characteristics of the nozzle by using a phase Doppler particle analyzer. Shahnaz [15] carried out experiments and found that, at lower fuel temperatures and fixed injection pressures, an increase in the mass flowrate and breakup length, a decrease in the spray cone angle and a liquid film breakup delay were observed due to the larger viscosity. Liu [16] conducted an experiment and simulation to research the influence of fuel temperature on the atomization characteristics of aviation kerosene under different pressures. He indicated that as the fuel temperature decreases from 50 °C to −50 °C, the fuel film thickness at the nozzle constantly increases, and the spray angle decreases under the same injection pressure. Zhao [17] used Fluent to simulate the steady spray field of the single-head model combustion chamber. The results show that when pressure and temperature decrease, the atomization cone angle tends to be smaller with the increased atomization droplets size because of the increase in the viscosity at low temperatures. Chen [18] used RP-3 aviation kerosene as the working medium to study the spray cone angle of a flared swirl nozzle. The VOF method with RSM (Reynolds stress model) was used to calculate the numerical flow, showing that the spray cone angle increases first and then decreases. Zhao [19] conducted multiphase simulation with the RSM turbulence model to numerically simulate the fuel swirl spray of an auxiliary nozzle. The mesh refinement near the nozzle, such as the swirl groove, can improve the simulation of fluid symmetry distribution that showed a good agreement with the experimental results.

In this paper, we mainly studied the changes in the atomization quality of the swirl spray of MMH propellants at low temperatures in order to comprehensively evaluate the feasibility of the low-temperature operation of bipropellant thrusters. The open code OpenFOAM was used to simulate the swirl spray characteristics of an MMH propellant at low temperatures. The simulation was carried out by the VOF method and LES turbulent model, and the gas–liquid interface in the computational domain was captured by the AMR (adaptive mesh refinement) method. The spray developments of the swirling atomization of the MMH at different temperatures, including the spray cone angle, liquid total surface area, breakup distance and droplets size distribution, were compared to determine if it can satisfy the requirement of a low-temperature propulsion system in the deep-space-exploration task.

2. Methodology

2.1. Governing Equation

The VOF model is employed to capture the gas–liquid interface when the liquid is sprayed out of the injector nozzle. The two-phase fluids in this paper are Newtonian fluids, and there is no heat exchange considered in the simulation. The governing mass and momentum conservation equations are presented as [20,21]:

$$\nabla \cdot \mathit{V}=0$$

(1)

$$\rho \frac{\partial \mathit{V}}{\partial \mathit{t}}+\rho \nabla \cdot (\mathit{V}\mathit{V})=-\nabla \mathit{p}+\mu \nabla \mathit{V}+\rho \mathit{g}+\mathit{F}$$

(2)

where V is the fluid velocity vector, ρ is the fluid density, p is the fluid pressure, $\mathit{F}$ is the surface tension, g is the gravity and μ is the fluid dynamic viscosity.

In the VOF model, ε_l represents the volume fraction to describe the interface between liquid and gas with a value between 0 and 1. The mixture density ρ and mixture viscosity μ can be calculated as:

$$\rho ={\epsilon }_l{\rho }_l+(1-{\epsilon }_l){\rho }_g$$

(3)

$$\mu ={\epsilon }_l{\mu }_l+(1-{\epsilon }_l){\mu }_g$$

(4)

where the subscript l represents liquid, and g represents gas.

The CSF (continuous surface tension force) model is used in the VOF model. The surface tension force term is shown as:

$$F={\displaystyle {\int }_{S(t)}\sigma kn\delta (x-x)}dS\approx \sigma k\nabla {\epsilon }_l$$

(5)

where k is the surface curvature, and σ means the surface tension coefficient.

2.2. Turbulence Model

The Naver–Stokes equation is a second-order nonlinear partial differential equation, so it needs to be solved by a numerical method with a turbulence model. The typical solution methods include LES, RANS (Reynolds-averaged Navier–Stokes equations), DNS (direct numerical simulation), etc. In this simulation, the LES method is employed to solve governing equations.

The LES method divides the turbulence into large-scale turbulent eddies and small-scale turbulent eddies by filtering. The direct solution method is used to calculate large-scale turbulent eddies; the SGS (sub-grid scale) turbulence model is introduced to simulate and solve small-scale turbulent eddies. The sublattice stress term is added for calculation using the k-equation model. The k-equation model is a single-equation model based on the eddy viscosity assumption, and the turbulent kinematic viscosity ν_t is calculated from the sub-grid turbulent kinetic energy k_sgs.

$$\frac{\partial {\tilde{u}}_i}{\partial x_i}=0$$

(6)

$$\frac{\partial (p{\tilde{u}}_i)}{\partial t}+\frac{\partial (p{\tilde{u}}_i{\tilde{u}}_j)}{\partial x_j}=\frac{\partial }{\partial x_j}(\mu \frac{\partial {\tilde{u}}_i}{\partial x_i})-\frac{\partial \tilde{p}}{\partial x_i}-\frac{\partial {\tau }_{ij}}{\partial x_j}$$

(7)

The two-phase fluid was regarded to be incompressible, so the incompressible Naver–Stokes equation is filtered to be described as [22,23]: where the sub-grid stress tensor ${\tau }_{ij}$ is obtained from the k equation model as [24,25,26]:

$${\tau }_{ij}=-2{\nu }_t{\tilde{S}}_{ij}+\frac{2}{3}{\delta }_{ij}{k}_{sgs}$$

(8)

$${\nu }_t=C_v\Delta {(k_{sgs})}^{1/2}$$

(9)

where ${\delta }_{ij}$ is the Kronecker symbol, ${\tilde{S}}_{ij}$ is the known strain tensor, ν_t is the turbulent eddy viscosity and ∆ is the filter scale.

The sub-grid turbulent kinetic energy k_sgs is calculated from its transport equation as:

$$\frac{\partial k_{sgs}}{\partial t}+\frac{\partial {\tilde{u}}_j{k}_{sgs}}{\partial x}=\frac{\partial }{\partial x_j}({\nu }_t\frac{\partial k_{sgs}}{\partial x_j})-{\tau }_{ij}{\tilde{S}}_{ij}-\epsilon$$

(10)

where the right side respectively represents the diffusion, production and dissipation terms of k_sgs, and $\epsilon$ means the dissipation rate of k_sgs, described as,

$$\epsilon =C_{\epsilon }\frac{(k_{sgs}{)}^{3/2}}{\Delta }$$

(11)

where $C_{\epsilon }$ and $C_{\nu }$ are constant coefficients, which need to be assigned or dynamically obtained.

2.3. Physical Model and Fluid Properties Assumptions

The 3D geometric model is mainly composed of an inlet, swirl chamber and atomization field. The spray domain is a combination of a cylinder and a frustum with a length of 15 mm and a diameter of 12 mm, as shown in Figure 1. The boundary conditions were set with a fixed inlet velocity of 3.89 m/s, a total pressure of 100 Pa at the spray domain and a no-slip wall boundary. All boundary conditions are set to follow the actual working operation of thrusters.

According to the research of Stephan [27] and NOAA (National Oceanic and Atmospheric Administration), the changes in the kinematic viscosity and surface tension of MMH over the range of temperature were obtained as shown in Figure 2. As the temperature decreased from 20 °C to −40 °C, the surface tension and kinematic viscosity of MMH gradually increased significantly. The surface tension increased from 54 mN/m to 64 mN/m, and the kinematic viscosity increased from 1.3 × 10⁻⁶ m²/s to 2.8 × 10⁻⁶ m²/s.

Due to the large-scale change in the spray domain, the extremely high requirement for the computing power makes it challenging to use structured grids with a very tiny size for capturing the liquid–gas interface in detail. Kevin [28] pointed out that an appropriate adaptive grid method can improve the computational efficiency in the spray atomization process. Therefore, the adaptive mesh refinement in OpenFOAM [29,30] was employed with a liquid–gas interface gradient, as shown in Figure 3. The initial grid side length is 0.1 mm after three refinements at the maximum, so the mesh size can be reduced to 0.0125 mm to adequately capture the droplets’ size of 50–100 μm, measured by the PDPA (phase doppler particle analyzer). For mesh independence validation, two initial mesh sizes of 0.12 mm and 0.08 mm were compared in MMH spray at 20 °C. The spray angles were the same, the maximum difference of the breakup length was 0.98% and the maximum difference of the droplets’ total surface area was only 0.42%. Figure 4 shows the comparison of several axial velocity distributions 3 mm away from the nozzle outlet, representing the initial mesh sizes of 0.08 mm, 0.1 mm and 0.12 mm at 20 °C. It can be seen that the axial velocity is symmetrically distributed along the central axis. The maximum deviation is less than 10% for the case with a 0.12 mm initial mesh. For cases with a 0.1 mm and 0.08 mm initial mesh, the different gap is reduced to less than 3%. Therefore, the initial mesh of 0.1 mm can be considered to meet the requirements of mesh independence.

In order to validate the accuracy of the simulation method and physical model, the simulation result was compared with the spray image by using a high-speed camera. Due to the toxicity of MMH, deionized water is used as a spray liquid to carry out the experiment instead. As shown in Figure 5, under the same inlet condition, the comparison shows an acceptable agreement in the conical liquid film of the umbrella-shape with a similar spray angle, breakup length and droplets distribution. For the spray angle, the tested angle is 69°, and the simulated one is 62°, with a 10.1% error. The error of the breakup length is only 5.31%. Therefore, the comparison indicates that the simulation method is suitable for atomization analysis [31].

2.4. Postprocessing Method

The total surface area of the liquid, denoted as S_spray, is obtained by integrating all areas of the spray interface. That is a very important parameter for evaluating the spray quality. For a fair comparison, the total area is normalized by a reference area. This reference area is chosen as the surface area of the spray cylinder jet with a closed end of S_ref with the diameter of the injector (D) and the penetration length of the spray (L) [32]. The ratio of the total surface area to the reference area is called the total surface area ratio, expressed as S.

$$S_{ref}=\pi DL+\pi D^2/4$$

(12)

$$S=S_{spray}/S_{ref}$$

(13)

Moreover, 24 slices around the central axis of the spray field are measured to capture the droplet size distribution and are calculated on average for testing the spray cone angle and the breakup distance, respectively. Figure 6 shows the distribution of 24 slices with an interval of 7.5 degrees.

3. Results and Discussion

In order to study the atomization characteristics of the swirl spray under different temperature conditions, the temperature of MMH was set to 20 °C, 0 °C, −20 °C and −40 °C in this study. Figure 7 shows the atomization comparison of the spray structure shape and spray field velocity distribution at 4.5 ms. In this Figure, the section in red indicates the current liquid phase in the spray, and the background in the yellow-blue color represents the spray velocity distribution contour. Inside the nozzle of the injector, MMH forms a circular liquid film with a rotational velocity under the effect of centrifugal force. As liquid is ejected without the limit of a nozzle wall, it turns into a hollow conical liquid film [3,33]. The liquid film thickness decreases gradually, and the surface wave develops rapidly with the expansion of the film. The conical liquid film breaks up into ligaments, including longitudinal liquid ligaments and transverse liquid ligaments, and then it develops with the primary breakup of larger droplets and the secondary breakup of finer droplets under a strong aerodynamic interaction [20]. From the perspective of energy, the spray develops with liquid pressure potential energy transformation into kinetic energy and then into surface energy [34]. In the swirl chamber, the energy loss mainly comes from the internal friction caused by the liquid viscosity and the direct friction between the liquid and swirling wall. After being ejected from the nozzle, with the growth of the liquid film surface area, the liquid surface energy also increases. At this point, the energy loss is mainly determined by gas–liquid interaction and surface tension force.

In Figure 7, as the temperature decreases from 0 °C to −40 °C, it can be found that the penetration distance of the rotating hollow conical film and the ejection velocity drop obviously due to the greater energy loss, as mentioned above. Moreover, in the process of liquid film formation, the liquid velocity gradient in the nozzle changes dramatically, and the viscous stress plays an important role. Nie [33] found that the liquid surface wave significantly decreases as the spray velocity decreases from 42.8 m/s to 22.8 m/s when studying the influence of different parameters on the spray of the conical liquid film. Therefore, a lower liquid-gas velocity difference causes more of an energy loss with a slight aerodynamic interaction, and the combined effect of a higher surface tension and viscosity further weakens the surface wave fluctuation at the same time. Under this condition, the spray tends to maintain the original liquid shape at lower temperatures, resulting in the breakup difficulties of liquid film, ligaments and droplets.

The basic characteristic flow numbers influencing the spray development process are the Reynolds number, Weber number and Ohnesorge numbers, which relate the critical nature of inertia, the viscosity forces, the surface forces and the ratio of inertial forces influencing jet breakup. These dimensionless numbers are calculated as follows. As shown in

Table 1. Spray velocity, Reynolds number, Weber number and Ohnesorge numbers at the nozzle outlet at different temperatures.

Temperature (°C)	Nozzle Outlet Velocity (m/s)	Reynolds Number	Weber Number	Ohnesorge Number
20	8.01	4209.70	523.09	0.005433
0	7.70	2519.58	456.96	0.008484
−20	7.57	1590.81	421.59	0.012907
−40	7.39	1054.16	388.69	0.018702

, these differences in velocity and liquid properties at different temperatures cause a significant change in dimensionless numbers at the nozzle outlet. That is bound to have a significant impact on subsequent spray flows.

$$\mathit{Re}=\frac{{\overline{\rho }}_l\overline{U}D}{{\mu }_l}$$

(14)

$$We=\frac{{\overline{\rho }}_l{\overline{U}}^{2}D}{\overline{\sigma }}$$

(15)

$$Oh=\frac{{\mu }_l}{\sqrt{{\rho }_{l}D\sigma }}$$

(16)

The effects on the total surface area ratio at different liquid temperatures are shown in Figure 8. It shows that the total surface area ratio at 20 °C is obviously greater over time compared with that at 0 °C and −40 °C. The reason is that a higher spray velocity provides more initial energy to enhance the liquid film surface fluctuations and gas–liquid interaction. Moreover, the larger surface area benefits the better evaporation, mixing and combustion of the propellant. When the spray progresses to 4.5 ms, the surface area ratio at 20 °C rises to 3.2%, compared with that at 0 °C, and the increase can reach up to 6.9%, compared with the case at −40 °C.

The atomization quality is further quantitatively described by calculating the droplets size distribution. At the moment when the atomized droplets fill the spray domain, the droplets number is statistically counted at different temperatures to generate a histogram. Then, the Boltzmann function is selected for fitting this histogram, with the correlation coefficient R² greater than 0.95 at different temperatures. The fitted curve is shown in Figure 9. Figure 9 shows that a lower temperature tends to produce a finer droplets distribution. In the range of 0–150 μm, the number of atomized droplets at 20 °C increased by 8.5% compared with that at 0 °C, and this number goes up to 25.6% compared with that at −40 °C. In the range of 150–250 μm, the number of atomized droplets at −40 °C was obviously greater than other cases. This result indicates that, at lower temperatures, the decrease in spray velocity and the increase in liquid surface tension cause insufficient kinetic energy and inefficiency in surface energy transformation. Therefore, a low temperature condition has a certain negative influence on the fineness of atomization [35].

A total of 24 slices around the central axis of the spray field are measured and calculated on average for testing the spray cone angle and the breakup distance, respectively. The results of the spray cone angle and the breakup distance over time are shown in Figure 10. The spray angle, which is mainly determined by liquid viscosity, surface tension, flow rate and spray pressure, can reflect the spray distribution characteristics in space, with a significant influence on the propellant mixing. As shown in Figure 10, the low temperature has a weak effect on the spray angle of MMH, with the smallest spray angle at −40 °C, 2.17% smaller than that at 20 °C. Therefore, it is concluded that the spray angle mainly depends on the structure of the injector design, with a limited difference over the low temperature range used in the deep space probe [36,37].

The breakup distance is one of the most important parameters, which is mainly related to the spray pressure, as shown in many experimental studies [38,39,40]. As shown in Figure 11, under the constant upstream pressure, the difference in the spray breakup distance at different temperatures is also not significant. This slight difference is mainly caused by the different spray velocity ejected from the injector nozzle. In the case at −40 °C, the spray velocity is also slower. So, there is a significant delay in the primary breakup compared with other cases, due to its insufficient kinetic energy with the weakest gas–liquid perturbation. As the liquid temperature rises from −40 °C to 20 °C, the breakup distance gradually decreases. When the temperature reaches 0 °C, the breakup distance is 1.3% shorter than that at −40 °C, and the difference in the breakup distance between 20 °C and −40 °C continually goes up to 3.1%.

4. Conclusions

The spray of MMH at low temperatures was investigated numerically based on the VOF method and the LES model, employing the open code OpenFOAM with the ARM method. A parametric comparison was conducted to investigate the spray development at different temperatures, including the velocity distribution, liquid total surface area, droplet size distribution, spray cone angle and breakup distance. The following conclusions are drawn:

(1)

The swirl spray of the MMH injector is distributed in the stable shape of a hollow cone spray, with a slight sensitivity to different temperatures. As the temperature decreases from 20 °C to −40 °C, the penetration distance and ejection velocity of the spray significantly reduce. This is because the increased kinematic viscosity and surface tension of MMH cause an extra loss in energy transformation from the pressure potential energy to the kinetic energy and then to the surface energy. This causes difficulties in the breakup of liquid films, ligaments and droplets.

(2)

The total surface area ratio of liquid increases obviously with rising temperatures. When the spray progresses to 4.5 ms, the surface area ratio at 20 °C rises to 3.2% compared with that at 0 °C, and the increase can reach up to 6.9% compared with the case at −40 °C. When the liquid temperature rises, a higher ejection velocity provides more energy for the liquid film surface wave fluctuation and gas–liquid dynamic interaction.

(3)

When atomized droplets fill the spray domain, the difference in droplets distribution is compared at different temperature. Under low temperature conditions, the number of small-sized droplets goes down, while the number of large-sized droplets grows significantly, and the spray quality is weakened.

(4)

The differences in the spray cone angle and breakup distance at different liquid temperatures are not significant, indicating that these characteristics are not sensitive to the change in physical properties and mainly depend on the injector structural parameters.

In summary, the low temperature condition has a certain impact on the spray breakup of the injector, which can lead to a degradation of the thruster performance. However, the spray shape, breakup distance and spray cone angle are slightly affected, with a limited influence on the subsequent mixing, combustion and liquid film cooling. Therefore, it shows that the bipropellant thruster can still work reliably to fulfill the deep-space-exploration task at low temperatures around −40 °C.