Influence of Different Altitudes on the Solid Rocket Contrail Formation in the Near Field

. Detecting the infrared characteristics of the contrails is a reasonable approach to tracing the rocket, and the particle properties of the contrails are the basis of the infrared analysis. The conventional numerical approach to obtaining the particle properties is a Euler/Lagrange method or a simple Euler/Euler method, di ﬃ cultly obtaining more accurate results because it ignores the particle size distribution in parcels or cells. A modi ﬁ ed Euler/Euler method is applied to simulate the contrail formation in the near ﬁ eld of a solid rocket motor at di ﬀ erent altitudes, which considers the size distribution by adding the ﬁ rst-to second-order particle radius moments based on the simple Euler/Euler method. The simulation results show that the crystals are generated at altitudes from 10 km to 20 km and that the contrails are visible at altitudes from 10 km to 15 km, where the radii of the crystals are from 0.1 μ m to 0.3 μ m. The visible contrails indicate that aviation vehicles are cruising at altitudes from 10km to 15km, and the smaller crystals indicate that the contrails are generated by rockets, not aircraft. Our work can provide important insight for the follow-up infrared analysis of the contrails based on the obtained particle radii.


Introduction
The contrail formation is due to the condensation of water vapor emitted by aviation vehicles [1]. Recently, with the increase of aircraft traffic, contrails produced by aircraft have contributed greatly to the global greenhouse effect, which has received much attention [2]. However, the ability to predict the contrail formation at both low and high altitudes for rockets is of interest for a variety of logistical reasons; for example, detecting the infrared characteristics of the contrails is a reasonable approach to tracing the rocket. Or rather, the information of the early contrails, such as the distances between the aviation vehicles and the contrail formation points, is most likely to expose the locations of the aviation vehicles. The objects of the recent research on the contrail formation are mainly the aircraft [3].
Liquid and solid particles in the plumes of jet aircraft cruising in the upper troposphere and lower stratosphere are the cores of the contrail formation [4], so the formation involves the wake flow dynamics, the microphysical dynamics, and the flow-microphysics couple process. In terms of the microphysical dynamics, the heterogeneous nucleation model is generally used to describe the formation process of ice crystals in contrails, and Fick's diffusion law is commonly adopted to simulate the growth process of ice crystals in contrails [5]. Smolders et al. used the model to simulate the wave-induced condensation and evaporation in a shock tube, from which numerical results obtained were close to the experimental results [6]. Khou et al. applied the model to the simulation of the ice crystal formation on soot particle cores emitted by the commercial aircraft, analyzing the distribution of the ice crystal radii [7,8].
The methods of simulating the wake flow dynamics coupled with the microphysical dynamics are mainly Euler/ Lagrange method and Euler/Euler method [9]. Paoli et al. applied the Euler/Lagrange method to simulate the contrail formation considering the interaction between the jet flow and the wing-tip trailing vortex and analyzed the characteristics of the ice particles in contrails [10]. Khou et al. employed a simple Euler/Euler method to assess the shape of the early aircraft contrail, which only considered the transport equation of the particle number density (the zeroth particle radius moment) [7,8].
A few scholars focus on rocket contrails (shown in Figure 1), and fewer scholars focus on rocket contrail formation. Voigt et al. used the Schmidt-Appleman criterion to judge the rocket contrail formation altitude, indicating that the contrails possibly formed only in the upper troposphere/lower stratosphere region [11]. Platov et al. investigated the rocket contrail formation in the upper atmosphere by solving the equations of thermal balance and mass balance of condensing particles, who found that the particle radius increase was only more than 7 nm which was much smaller than the increase in the aircraft cases [12,13]. Chenoweth et al. used the heterogeneous nucleation model to simulate the rocket contrail formation based on the CRAFT CFD code but did not present a specific approach [14]. Dalin et al. used detailed photographic imaging taken from the ground to analyze the dynamics of Soyuz rocket exhaust products at the initial trajectory in the mesosphere and lower thermosphere, which only focused on the later stage development of the formed contrail [15].
The above investigations of the rocket contrails could not completely describe the characteristics of the contrails as the study of aircraft contrails due to the limited one-dimensional model. The jet plumes produced by rocket motors have higher temperatures than those produced by aircraft engines [17], which makes the contrails more difficult to form, and the different shock waves make the plume shapes different, which affects the contrail shapes. To overcome the defects of Khou et al.'s method which ignores the particle size distribution in cells, we apply a modified Euler/Euler method with the heterogeneous nucleation model to simulate the contrail formation in the near field of the solid rocket motor, which extends Khou et al.'s method by calculating the transport equations of 0th to 2nd moments of the particle size. The modified method has a higher accuracy for calculating the mass transfer rate due to the high order moments which can introduce the particle size distribution in cells. The simulation applies the real gas components in the solid rocket motor combustor and considers the influence of the fuselage shape on the jet plume.

Theoretical Modeling
Three models are applied in the investigation, including the fluid flow model, the ice formation microphysical model, and the gas-particle coupling model. The fluid flow model is used to simulate the mixing of the hot jet exhaust with the ambient air disturbed by the turbulence of the fuselage. The ice microphysical model is used to describe the birth and growth of ice particles in the contrails. The gasparticle coupling model is used to describe and solve the jet flow coupling with the particle motion.
2.1. Fluid Flow Model. The plume mixture is assumed to be composed of air, gas, and soot particles. The relaxation time of the particles is far less than the characteristic time of the jet flow because the particle sizes are about 1 μm, so the jet and particles are assumed to be in dynamic and thermal equilibrium [18]. Thus, the particle phase has the same velocity as the carrier gas phase. The mass fraction of the particles is small enough to assume that the particles are a species of the fluid without impacting the fluid flow characteristics. The axisymmetric, multispecies, and compressible governing equations used in the flow model are Figure 1: The rocket contrail [16].
2 International Journal of Aerospace Engineering where x is the axial coordinate, y is the radial coordinate, ðu, vÞ T is the velocity vector, p is the pressure of the mixture (including gas and particles), E is the total energy of the mixture, ρ is the mixture density, Y i is the mass fraction of the i -th species, S is the source term vector described by the gasparticle coupling model, and D is the diffusion term vector including the viscosity effects, the heat conduction, and the diffusion of species with the turbulence, detailed in Ref. [7]. The turbulence is simulated by the spatial RANS model to reduce the computation cost without reducing the accuracy. The Reynolds tensors in the diffusion term vector are given by the Boussinesq hypothesis and the two equations realizable k-ε model to adapt to the flow separation on the nozzle exit with the supersonic flow [19], detailed in Ref. [20].

Microphysical Model.
Before describing the microphysical model, some assumptions are introduced. Inert particles are assumed to be spherical for the sake of simplicity, as their fractal-like shapes are too complex to be counted. The condensation is supposed to be immersion freezing from a thin liquid layer, and the turbulence is assumed to not promote the ice growth in any particular direction.
The main assumption used in this work is that contrail formation is driven by ice heterogeneous nucleation on inert particles. Although all inert particles do not necessarily freeze, an upper limit case can be selected by first assuming that they are all activated. All emitted inert particles may therefore be considered as ice nuclei. A kinetic model has been used for the deposition of water vapor onto ice particles. The particle growth is evaluated using a modified Fick law dedicated to mass transfer on particles of which the radius r p is of the order of the mean free path λ. The ice generation rate _ ω ice is given by where N p is the ice particle number in the unit volume, r p is the radius of the ice particle, M w is the water molar mass, R is the ideal gas constant, D vap is the diffusion coefficient of water vapor to the air, p vap is the water vapor partial pressure in the plume, p sat/ice vap is the saturation vapor pressure above the ice, p sat/liq vap is the saturation vapor pressure above the liquid, Π is the function judging whether the vapor mass transfer onto the dry particles, and Gðr p Þ is the function taking into account the transition of the uptake from the gas kinetic to diffusion regimes near the nanometric particle, given by , Π = 0 when p vap ≤ p sat/liq vap and r p = r s , 1 when p vap ≥ p sat/ice vap or r p > r s , where α is the deposition coefficient of water molecules on ice, α = 0:1, r s is the initial particle radius, and Knðr p Þ is the particle Knudsen number, Knðr p Þ = λ/r p where λ is the mean free path of a gas molecule. The rest functions, pressures, and coefficients are detailed in Ref. [7]. The ice crystal growth rate is the function of r p , given by 2.3. Gas-Particle Coupling Model. The Eulerian bulk method, solving for the moments M k of the size distribution, is used as the gas-particle coupling model. The moments are defined as by assuming a certain shape for the size distribution f that depends on a limited number of parameters, including the local particle number density, the mean particle radius, and the standard deviation of the particle radii. According to the moment definition, the zeroth moment M 0 is the local particle number density, the 1st moment M 1 is the half total length of the particles per unit volume, the 2nd moment M 2 is 1/ð4πÞ times the particle total surface area per unit volume, and the 3rd moment M 3 is 3/ð4πÞ times the particle phase volume fraction. The 1st to 3rd moments are the particle radius statistics, which can better reflect the information of the particle size distribution in cells. Khou et al.'s method only considered the zeroth moment and ignored the radius distribution in cells throughout the simulation, so the accuracy of the simulation was relatively low. To overcome the defects of his method, the vectors in the governing equations extending the transport equations for the moments are defined by International Journal of Aerospace Engineering by ignoring the effect of drag force in atmospheric applications and only considering the nucleation source. The diffusion terms in the transport equations for the moments are 0. The detailed forms of the source terms are by assuming that the nucleation rate is 0 due to the initial particle radii greater than the nucleation critical radius r * obtained by the Kelvin equation, where _ r ≡ dr/dtð rÞ, ρ ice , and L ice are, respectively, the ice crystal condensation/sublimation rate, the ice density, and the latent heat of condensation, related with the local temperature, with the "surface-averaged" radius r ≡ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi

Numerical Setup
The investigation involves modeling and simulating the external flow field of a rocket, whose performance parameters are provided by the third party [22]. The axisymmetric profile of the rocket is shown in Figure 2. The rocket has a Laval nozzle with a combustion chamber and a fuselage without wings. The performance parameters of the nozzle are designed as follows: a trapezoid profile of the expansion, the area expansion ratio being 16.2, and the nozzle being conical with an exit half angle of 18°. The profile of the fuselage is a combination of a rectangle and a cone. According to the axisymmetric rocket, the computational domain is modeled as a 2D axisymmetric domain, as shown in Figure 3. The domain dimensions are given by referring to the length and diameter of the fuselage. The domain shape is a rectangle with the radial dimension being 45 times the diameter of the fuselage and the axial dimension being 15 times the length of the fuselage.
Based on the 2D axisymmetric domain, the domain mesh is partitioned by some blocks paving unstructured quadrilaterals to adapt to the irregular shape, as shown in Figure 3. The internal and adjacent mesh of the nozzle domain is refined to capture the shock wave generated by the jet, where the maximum grid spacing Δ is set to 10 -2 m. The mesh far from the rocket is gradually coarsened to improve calculation efficiency, where the maximum grid spacing Δ is set to 0.5 m. The total number of cells is 0.35 million. To ensure the grid independence, four sets of meshes (0.25 million, 0.30 million, 0.35 million, and 0.40 million) are calculated with a simple equivalent gas model without reactions. The axial vapor profiles are compared in Figure 4. It has been indicated that the derivation is less than 8% between 0.35 million and 0.40 million cell meshes, so the grid with 0.35 million cells is applied.
The initial and boundary conditions are set in the 2D axisymmetric solver. The boundary conditions, shown in Figure 3, are enforced as follows: the nozzle inlet set as the pressure inlet, the domain exit set as the pressure outlet, the domain inlet set as the pressure far field, the fuselage and the nozzle wall set as the adiabatic wall, and the International Journal of Aerospace Engineering symmetry axis. The pressure inlet is also the particle inlet. The pressure inlet conditions depend on the work state of the combustion chamber. The pressure outlet conditions and the initial conditions correspond to cruise conditions in the upper troposphere. The ambient temperature, pressure, and sound speed of the upper troposphere ranging from 5 km to 25 km are listed in Table 1 by using the NACA Standard Atmosphere [23]. The initial conditions of the domain other than the chamber zone are set as the ambient conditions, and the initial conditions of the chamber zone are set as the pressure inlet conditions of the beginning, listed in Table 2, where the particle density is detailed in Ref. [24]. The component equilibrium equations of double-base propellant combustion products are solved by the principle of minimum Gibbs free energy [25], and the mass fractions of the major combustion products in the combustion chamber are listed in Table 3. The particle radii in the      [1] indicated that the threshold temperature was less sensitive to the relative humidity when the relative humidity (RH) was under 60%, and Khou et al.'s work showed that the locations of supersaturated areas with respect to liquid water were similar at RH = 30% and RH = 60% [7]. We assume that RH keeps constant at 20% and that the cruise speed is set as Ma = 1:6 at five altitudes, including 5 km, 10 km, 15 km, 20 km, and 25 km.

Model Validation
The investigation validates the theoretical models by simulating the wave-induced condensation and evaporation. Details of the wave-induced experimental setup, conditions, and procedures are to be found in Ref. [6]. The validating model, with two valves whose sizes are the same as the tube, is simplified by generating the rarefaction waves of the same intensity as the experiment, shown in Figure 5 The droplet radii of the simulation and the experiment, obtained at the observation window located at the lowpressure chamber of 6.63 m from valve 2, are shown in Figure 6. Figure 6 indicates that the virtually same droplet growth trends are obtained by the numerical method and the experimental method. The mixture is initially unsaturated, and the particles keep the original size. The saturation ratios at the observation window start to rise with the passage of the first rarefaction wave. When the saturation ratios rise to the critical saturation ratios, vapor starts to condense on the particles, and the particles grow up. Vapor is consumed and the heat is released with the condensation, resulting in the fact that the mixture tends to an equilibrium state where droplet radii are about 0.8 μm without increasing. The reflected wave, generated through the front of the rarefaction wave impinging valve 1, again disturbs the state of the mixture, and the droplets continue to grow up to the mean radius of approximately 1 μm. The state of the gaseous carrier changes discontinuously to an unsaturated state, and the droplets start to evaporate with the shock wave passing the observation window. The pressures and temperatures of the mixture change continuously until the droplets have completely disappeared.
The maximum relative error between the numerical and experimental droplet radius occurs in the balance state, with a value of under 10%. The maximum numerical mean droplet radius is 1.05 μm, slightly greater than 1 μm of the experimental mean radius. The virtually same droplet growth trend and the enough small maximum relative error indicate that the models are validated and adapted to the investigation.

Results and Discussion
5.1. Wake Dynamics. The axisymmetric realizable k-ε model is used to resolve the turbulence induced by the supersonic flow passing through the fuselage, which disturbs the jets mixed with the ambient air. The simulations are set to transient due to the particle growth with time. The first-order implicit scheme is used to discretize the time, where the Courant number is set to one. The spatial scheme is set to the second-order upwind scheme. To compare the steady contrails at the different altitudes, we regard that a quasisteady state is reached when the fluid parameters are near the downstream edge of the computational zone level off over time. The following results are based on the quasisteady state, and the convergence curve is shown in Figure 7.
The vapor partial pressure is related to the ice generation rate, recalling from Equation (2). The mixture pressure fields at the 5 altitudes are firstly analyzed due to the vapor partial pressure depending on the mixture pressure, as shown in Figure 8. The nozzle pressure ratios (NPRs) only depend on the ambient pressures at the different altitudes, due to the nozzle throat flow keeping supersonic. The NPR ranges from 1.29 up to 27.36 due to the ambient pressures dropping rapidly with the altitude from 5 km to 25 km. The NPR criterion indicates that the jet regime develops from the lowly underexpanded jet to the super highly underexpanded jet with the rocket flying from 5 km to 25 km. The obvious deviations between the mixture pressures and the ambient pressures focus on the jet core zone, where the centerline jets expand and compress repeatedly. The rarefaction waves  Figure 5: The geometry of the validation case. 6 International Journal of Aerospace Engineering generated by the nozzle exit and the fuselage end make the jet and the inflow excessive expansion, so the pressures drop below the ambient pressures, and the dimensionless mixture pressures are negative in the vicinity of the nozzle exit. The areas of the pressure deviation zone (the zone where the local pressures are more than 110% of the ambient pressures) behind the nozzle exit are enlarged with the jet underexpansion aggravation. The pressures at the mixed layers of the jets recover to the ambient pressures through a step process due to the intercepting shock, as shown in the bottom curves of Figure 8. The pressures at the zones except the pressure deviation areas are approximately the ambient pressures, where the ice particles are likely to form. Thus, the mixture pressure of the ice generation is the ambient pressure. Given that the saturation vapor pressures above ice and liquid are the function of the temperature, the temperature fields at different altitudes are extracted as shown in the top half panels of Figure 9. The temperatures rapidly drop to the ambient temperatures at the edges of the jet, which forms the outer shear layer. The temperature drop along the centerline is similar to that in the aircraft cases, but the centerline temperatures are more difficult than the aircraft cases to recover to the ambient temperatures due to the much higher total temperature [26], as shown in the bottom curves of Figure 9. The above zero temperature zones (red zones in the top half panels of Figure 9) shrink with the downstream distance and broaden downstream with the altitude increase due to the jet underexpansion aggravation. Therefore, the contrail is not likely to form in the centerline zone, because the threshold temperature of the contrail formation is below zero.
To find the zones favorable to condensation, we overlap the isocontour lines of the vapor local density and the temperature contours together in the top half panels of Figure 9. The vapor local density at the lines is much greater than the ambient value, as shown in the vapor local density contours (shown in the bottom half panels of Figure 9), and the temperatures in the green zones of the top half panels are 7 International Journal of Aerospace Engineering much lower than the centerline values, lower than 240 K. The isocontour lines are overlapped with the green zones at the altitudes of 10 km, 15 km, and 20 km, while these are not overlapped at the altitudes of 5 km and 25 km. The contrail is not likely to form at the altitude of 5 km, in which the ambient temperatures are too high and the vapor local density at the ambient zone is too low, although the lines interact with the ambient temperature zone. The contrail is also not likely to form at the altitude of 25 km, in which all lines intersect with the red zone of the temperature field. Therefore, the contrail is only possible to form at the altitudes of 10 km, 15 km, and 20 km, in which the temperatures and the vapor local densities at the jet edges are appropriate.
The ice generation rates and areas of the ice generation zones (the white zones in Figure 10) decrease with the alti-tude in the range of [10,25] km, which is in consistency with the analysis of the temperature and the vapor local density. The vapor partial pressures p vap does not exceed p sat/liq vap everywhere at the altitudes 5 km and 25 km, while p vap only exceeds p sat/liq vap on the white zones at the altitudes from 10 km to 20 km. The particles grow when they are passing the white zones with the flow, and their radii exceed the initial radius r s when they leave the white zones, causing that the judging function Π equals to 1 on the zones except for the white zones. The vapor partial pressures p vap on the zones except the white zones is less than p sat/liq vap , so the ice generation rates on the zones are mostly negative where the crystal radii exceed r s . The vapor local densities at the isocontour lines intersected with the green zones at the    International Journal of Aerospace Engineering altitude of 15 km are lower than the values at the altitude of 10 km, causing that the ice generation rates in the white zones at the altitude of 15 km are less than the values at the altitude of 10 km. According to the Schmidt-Appleman criterion, the curve of threshold temperature T LC versus altitude predicts that threshold temperatures are higher than the local temperatures at the altitudes from 8 km to 22 km (shown in Figure 11), indicating that the ice particles can grow at these altitudes, which is in consistency with the results of the ice generation rate contours.

Contrail Microphysics.
Heterogeneous condensation, the process that vapor condenses on inert particles when the local vapor saturation ratios reach the liquid water saturation ratios, is the only condensation type considered in the investigation, since homogeneous condensation only occurs under the conditions with much higher saturation ratios. According to the flow analysis, the ice particles hardly appear in the computational domain at the altitudes of 5 km and 25 km, due to the higher ambient temperatures for 5 km and the lower ambient pressures for 25 km, and appear at the edges of the jets at the altitudes from 10 km to 20 km, due to the intermediate condition.
The ice concentrations at altitudes of 10 km, 15 km, and 20 km are above zero, depicted in the top half panels of Figure 12. The ice concentrations at the altitude of 20 km are approximately zero, since the ice generation zone is too narrow and small and that the ice generation rate is too low. The ice concentrations at altitudes of 10 km and 15 km are approximately uniform in the nonzero ice concentration profiles, approximately 1.5e-5 kg/m 3 . The nonzero ice concentration profile at the altitude of 10 km starts from approximately 15De behind the nozzle exit and extends to approximately 90De downstream of the nozzle exit. At the altitude of 10 km, the maximum thickness of the profile is approximately 2De at x = −32De, and the maximum outer radius of the zone is close to 7De at the end of the zone.
Although the ice generation rates and areas of the ice generation zones at the altitude of 10 km are greater than those at the altitude of 15 km, the ice sublimation rates on the sublimation zones at the altitude of 10 km are much greater than those at the altitude of 10 km. When the crystals pass the sublimation zones, the crystal radii decrease rapidly (shown in the top half panels of Figure 13). So compared to the ice concentration profile at the altitude of 10 km, the profile at 15 km is wider and longer. At the altitude of 15 km, the nonzero concentration profile is approximately a triangle with a sharp angle from x = 20De to the end of the computational domain, whose maximum thickness and maximum outer radius are 11De and 14De, respectively, at the end of the domain.
According to Ref. [27], we use the least ice concentration as the contrail visible criterion. That is, the contrail is visible when the ice concentration is more than 4e-3 g/m 3 . Based on the ice concentration contours (shown in the top half panels of Figure 12), the visible contrail shapes are the white shapes shown in the bottom half panels of Figure 12. The contrail shapes are consistent with what can be observed in photographs of contrails at cruising altitudes [16]. The visible contrail is hardly formed at the altitude of 20 km, although the ice particles are formed. Visible contrails expand backward and sideways with the altitude from 10 km to 15 km, which is consistent with the results of Ref. [11]. The visible contrails are nonpersistent in the domain at the altitude of 10 km and persistent in the domain at the altitude of 15 km, indicating that the lengths of the visible contrails increase with the altitude from 5 km to 15 km and decrease with the altitude from 15 km to 25 km. Visible contrails begin to appear at the altitude of 10 km and grow at approximately 16 m behind the nozzle exit at the altitude of 10 km and grow at approximately 16De behind the nozzle exit at the altitude of 15 km, no more than the distance of 40De in the aircraft case [7]. The observed contrail diameters, the outer ring diameter of the visible contrails, are approximately The excess water vapor is immediately condensed into the ice particles after the formation of the ice nuclei. The particles grow until the local relative humidity with respect to the ice drops to 100%. The mean radii of the ice particles reach the maximum at the edges of the contrails, where the local relative humidity with respect to ice is slightly less than 100 percent. The mean ice particle radii decrease with the altitude increasing from 10 km to 25 km. The mean ice particle radii in the visible contrail are about 0.25 μm at the altitude of 10 km and about 0.15 μm at the altitude of 15 km, as shown in Figure 13. The mean radii of the ice particles exhibit nearly 1.4 μm in the aircraft cases [26], almost 10 times the ice particle radii in the investigation. The ice particle radius results indicate that detecting the maximum ice particle radius in the contrails can reveal the type of air vehicle. The maximum radius greater than 1 μm reflects that the air vehicle is an airplane, otherwise a rocket.

11
International Journal of Aerospace Engineering

Conclusion
In the present work, the contrail formation of the rocket motor is simulated by the modified Euler/Euler method at five altitudes from 5 km to 25 km, which considers the turbulence of the complete rocket geometry. The modified Euler/ Euler method is validated by the wave-induced condensation and evaporation experiment, and the simulation results are indirectly validated by the Schmidt-Appleman criterion. The wake results show that the flow parameters (temperatures and vapor partial pressures) are unfavorable to the vapor condensation near the axis zone. The contrail results show that the contrails can be visible at altitudes ranging from 10 km to 15 km and that the mean crystal radii range from 0.15 0.15 μm to 0.25 μm in the visible contrails. Compared to aircraft contrails, the rocket contrail results suggest that aviation vehicles are cruising at altitudes from 10 km to 15 km when the visible contrails are observed and that the contrails are generated by rockets when the smaller mean crystal radii are detected in the contrails.
Most notably, the modified Euler/Euler method extends Khou et al.'s method by calculating the transport equations of zeroth to 2nd moments, which improves the precision of the mass transfer rates. The present work demonstrates the feasibility of performing the simulations of the contrail formation for the solid rocket motor. The obtained results may preliminarily reveal the locations of the rocket and further serve as the input conditions for the infrared analysis of the contrails. The simulation is based on the known initial particle density and the initial particle radii with the uniform distribution, but the initial particle parameters vary with the motor work state. Therefore, the future study will consider the effect of the motor work state on the initial particle characteristics.