CFD Study of Base Drag of the Grot Rocket

The rocket used for the study is a sounding rocket developed by the Students’ Space Association, known as Grot. Grot is a 2-m-long dart-type rocket (Fig. 1), which comprises two stages. The first stage accelerates the dart (second stage) to an altitude of 2 km, where aerodynamic separation takes place and the motor (first stage) plunges back to the ground, while the dart continues flying to an altitude of 18.5 km. Key parameters describing the Grot sounding rocket can be found in Table 1.

Figure 1.

Various types of nozzles have found use in rocketry; most configurations consist of converging and diverging sections, with the throat located at the minimal cross-sectional area. The variations between the types usually concern the diverging section, as it is determined to be the most crucial for achieving high performance.

The conical nozzle (Fig. 2) is one of the simplest and oldest geometries. One can define a theoretical correction factor λ [7] for the ideal conical rocket nozzle: 1 λ=12(1+cosα) $$\lambda \, = \,{1 \over 2}(1 + \cos \alpha )$$ where α is the divergence cone angle. The correction factor determines the ratio between the momentum of the gases in a nozzle with a nozzle angle 2α and the momentum of an ideal nozzle with all gases flowing in the axial direction. It should be noted that it applies only to the momentum thrust, not the pressure thrust.

Figure 2.

Large divergence angles enable light and short designs at the expense of performance, while small angles lead to the majority of the momentum being axial, which in turn improves the specific impulse, but the increased length influences the complexity of the design and the weight of the propulsion system and vehicle. The optimal setting is usually found to have a half-angle between 12° and 18°, depending on the application and mission profile. The Grot rocket uses a conical nozzle with a half-angle of 15°.

Thrust is a reaction force generated when an exhaust plume is ejected from a motor; this in turn applies a force pushing the motor in the opposite direction, as described by Newton’s third law. The inertial term of the thrust equation can be derived from Newton’s second law, as shown in Eq. (2): 2 F→=dp→dt→F→=d(mv→)dt→F=m˙v→. $$\vec F\, = \,{{d\vec p} \over {dt}} \to \,\vec F\, = \,{{d(m\vec v)} \over {dt}} \to \,F\, = \,\dot m\vec v.$$

An additional source of thrust is the pressure force acting on the nozzle exit. The full thrust force equation for a rocket motor can be seen in Eq. (3). 3 T→=m˙⋅v→e+(p−pa)Ae→T=p⋅v→e⋅|v→e|⋅Ae+(p−pa)Ae. $$\vec T\, = \,\dot m\, \cdot \,{\vec v_e}\, + \,(p - {p_a}){A_e} \to \,T\, = \,p \cdot \,{\vec v_e} \cdot |{\vec v_e}| \cdot \,{A_e}\, + \,(p - {p_a}){A_e}.$$

In the more accurate integral form, Eq. (3) can be presented as Eq. (4): 4 T→=−∫∂τρ⋅v→e⋅|v→e|dS−∫∂τ(p−pa)dS. $$\mathop T\limits^ \to \, = \, - \mathop \smallint \limits_{\partial \tau } \rho \, \cdot \,{\mathop v\limits^ \to _e}\, \cdot \,|{\mathop v\limits^ \to _e}|dS - \,\mathop \smallint \limits_{\partial \tau } (p - {p_a})\,dS.$$

For an ideally expanded nozzle, the second term tends to zero; nozzle expansions are illustrated in Fig. 3.

Figure 3.

The Grot sounding rocket with the directions of the relevant axial forces is shown in Fig. 4. In the following discussions, it shall always be assumed that any component of drag is defined as positive when acting in the direction of the flow velocity v→∞${\vec v_\infty }$ (i.e. slowing the rocket down).

Figure 4.

The axial forces present during the flight of the rocket vary depending on the Mach number and operation status of the propulsion system. Therefore, two cases were distinguished, the first being the non-operational propulsion system and the second being the operational propulsion system.

The National Aeronautics and Space Administration’s (NASA) computer programme Chemical Equilibrium with Applications (CEA) was designed to calculate the chemical equilibrium compositions and properties of complex mixtures [8]. For this study, the theoretical rocket performance mode was utilised, knowing the product concentrations from the set of reactants that made up the propellant in the form of fuel and oxidiser. The thermodynamic and transport properties can be obtained for the product mixture at the exhaust of the motor. The motor uses a typical AP+HTPB+AL composition, with the exact composition being confidential.

The considered geometry was prepared on the basis of the computer-aided design (CAD) model of the Grot sounding rocket and technical drawings of the nozzle (Fig. 5). Given that this case study involved an axisymmetric model, only half of the geometry was modelled in Ansys® DesignModeler. The outline of the Grot sounding rocket without fins in the axisymmetric view can be seen in Fig. 6.

Figure 5.

Figure 6.

In addition, a domain study was conducted in order to minimise and preferably eliminate the influence of domain size on the results; the selected domain size can be seen in Fig. 7.

Figure 7.

The overall quality of the mesh and the solver settings determine the accuracy of the results. However, the mesh used in a study may not require a large number of elements in order to obtain a fairly accurate solution. For this reason, a mesh independency study was conducted in order to determine the appropriate mesh settings so as to have minimal effect on the accuracy of results; as a result, the final mesh contained 626,458 elements. In addition, a structured mesh is favourable for increasing the efficiency of the calculations as it requires less memory. The nozzle and the plume region experience high gradients of pressure, temperature and Mach number; therefore, it was advantageous to utilise a dense-structured mesh in these areas; for this purpose, it was necessary for the fluid domain to be split into regions (Fig. 8).

Figure 8.

Secondly, performing external aerodynamic CFD simulations requires the modelling of the boundary layer during meshing; this is done in the form of inflation layers. For this to be done correctly, one must estimate the boundary layer thickness formed in the range of the Reynold numbers experienced during the flight of the vehicle [9] and model the inflation layers to contain the approximated boundary layer thickness. It was decided to keep the 30 ≤ y⁺ ≤ 300 to implement the wall function approach for the shear stress transport (SST) k-ω turbulent model within Ansys® Fluent, as this provides a shorter calculation time and lower number of mesh elements compared to y⁺ ≤ 1, while keeping the results relatively accurate. A similar choice of y+ has been observed in a previous paper [10]. For the nozzle, a bias factor was applied for capturing the boundary layer, followed with a dense mesh in the plume region. The mesh for the aforementioned domain regions was divided in the same manner (Figs 9 and 10).

Figure 9.

Figure 10.

By using the NASA CEA software, the physical properties of the exhaust gases were obtained as a function of motor operation duration. For this case study, the input chamber pressure was known beforehand. This case study made use of the species transport model, since it models the interaction between the exhaust gases and the free-stream air. The motor operating point chosen for the aerodynamic analysis is the one providing the highest value of the thrust force.

Compressibility effects are prevalent when the fluid accelerates beyond Mach 0.3; thus, it was necessary to use the density-based solver in this study. As mentioned in Subsection 3.2, the SST k-ω turbulent model was utilised, which was found to be more favourable in previous works [10,11]; alongside the use of the energy equation, the ideal gas model was used, supplemented with the Sutherland model for viscosity. In addition, the second-order upwind scheme was chosen for all convective term flow variables as this provides higher accuracy. The study focuses on atmospheric conditions at sea level with a far-field pressure boundary condition at the domain inlet (p_a = 101,325 Pa and T_a = 288.15 K [15°C]). In addition, turbulent intensity was 2.16% and turbulent length scale was 0.012 m. Finally, solution steering and full multi-grid (FMG) initialisation [12] were both activated.

For qualitative analysis of the obtained flow fields, the Mach number and static pressure, the contour plots for the operational and non-operational propulsion systems were analysed for the Mach number 1.2. The shock waves developed on the rocket appeared dimmer in the Mach number contour plot for the active motor than for the inactive motor (see Figs 11 and 12). This was due to the exhaust plume travelling at a Mach number three times greater than the free-stream air. Figures 13 and 14 present the static pressure for the same Mach number. Figures 15 and 16 show the closed-up static pressure contour plots for the aforementioned contour plots. It is important to note that the static pressure near the aft end (as defined by the red lines in Fig. 17) is lower in the engine-off case.

Figure 11.

Figure 12.

Figure 13.

Figure 14.

Figure 15.

Figure 16.

Figure 17.

To evaluate the thrust force, it was necessary to define the nozzle exit surface, as the calculations were made for an axisymmetric model. The exit surface of the nozzle is represented by a black line in Fig. 17.

The approach utilises the Ansys® Fluent post-processing surface integral reports and Eq. (3) presented in Subsection 2.3. However, to be able to use such an approach, one needs to define the inertial term from Eq. (3). For this purpose, a variable, ρ⋅v→e⋅| v→e |$\rho \cdot {\vec v_e} \cdot |{\vec v_e}|$, representing the custom field variable is adopted. Afterwards, surface integrals of the custom field variable and the static pressure, which is the second term in Eq. (3), were calculated for the nozzle exit area. By summing up these two results, the thrust was obtained for each engine operating time step. As stated in Subsection 3.3, the maximum thrust operation point was utilised as the reference for the drag calculations, and this point gives a total effective thrust force of 7,750 N. As expected, the thrust value does not vary significantly with the Mach number.

The raw data for the total drag coefficients for active and inactive propulsion systems, along with the difference between them, are presented in Table 2 and Figs 18 and 20.

M, free-stream Mach number; CDon$C_D^{on}$, drag coefficient computed for the operating propulsion system; CDoff$C_D^{off}$, drag coefficient computed for the non-operating propulsion system; ΔC_D, difference between drag coefficients of operational and non-operational propulsion systems.

Figure 18.

Figure 19.

Figure 20.

The total drag coefficients for both computational cases (active and inactive propulsion systems) are shown in Fig. 18. As expected, they are not equal; the difference is presented in Fig. 20. Further discussion of the results necessitates the introduction of some additional coefficients. The coefficient of drag acting on the nozzle exit (already defined in Eq. (6)) shall be denoted as CD,Noff$C_{D,N}^{off}$]. In the case of the aft end, the drag coefficients for the engine-on and engine-off cases will be denoted as CD,Aon$C_{D,A}^{on}$ and CD,Aoff$C_{D,A}^{off}$, respectively. In literature, both components (the drag acting on the nozzle exit and the aft end of the rocket) are being considered together and referred to as base drag. For the sake of comparison, the base drag coefficient of the Grot rocket can be defined as follows: 11 CD,Bon=CD,Aon $$C_{D,B}^{on} = C_{_{D,A}}^{on}$$ 12 CD,Boff=CD,Aoff+CD,Noff $$C_{D,B}^{off} = C_{_{D,A}}^{off} + C_{D,N}^{off}$$

In Eq. (11), the pressure acting on the nozzle exit is not taken into account, since it is already a component of thrust, as defined by Eq. (4). The C_D,B base drag coefficients computed during post-processing of the CFD study are directly compared with the analytical results obtained previously [4] (Fig. 19). It can be noted that the curves are in relatively good agreement in the supersonic region, especially given that they are not describing exactly the same rocket. In the transonic region, the curves diverge, although this is mainly due to the analytical model being based on second-order supersonic flow solutions and thus being valid only in the supersonic domain.

The established ΔC_D(M) dependency does not seem to be easily modelled by simple physical models. However, it can be demonstrated that it can be fully described by two factors: the difference in taking into account the nozzle exit area (as explained in Subsubsection 2.4.2) and the change of pressure in the aft region. This hypothesis can be represented as follows: 13 ΔCD≈CD,Noff+CD,Aoff+CD,Aon. $$\Delta {C_D} \approx C_{D,N}^{off} + C_{_{D,A}}^{off} + C_{D,A}^{on}.$$

Indeed, the plot of both sides of the approximate equality in Fig. 21 proves that they are indeed almost identical. Consequently, it is confirmed that the influence of propulsion system operation on drag can be fully explained by the pressure acting on the aft end of the rocket.

Figure 21.

Since the atmospheric parameters change with flight altitude, it is worth considering whether they affect the previously shown results. Additional numerical simulation for the flight altitude of 5 km was performed. For this altitude, the static pressure of 54,050 Pa and static temperature of 255.7 K (−17.45°C) were assumed. Figure 22 presents the comparison of ΔC_D between the results for this altitude and sea level. ΔC_D decreases with altitude in a proportional manner. The main reason for this difference is that the change in pressure and density influences the compressible behaviour of the flow. In addition, viscous effects are more pronounced due to the decrease of Reynolds number at a given Mach number for the 5 km case. However, upon analysing the CFD results, it became clear that the difference in viscous drag is negligible in comparison to the effects of pressure. Further analysis is required in order to model the dependence of ΔC_D on flight altitude.

Figure 22.