Prediction and Analysis of the Aerodynamic Characteristics of a Spinning Projectile Based on Computational Fluid Dynamics

Numerical simulations of a spinning projectile with a diameter of 120mm were conducted to predict the aerodynamic coefficients, and the CFD results were compared with the semiempirical method, PRODAS. Six coefficients or coefficient derivatives, including zero and the quadratic drag coefficient, lift force coefficient derivative, Magnus force coefficient derivative, overturning moment coefficient, and spinning damping moment coefficient, which are important parameters for solving the equations of motion of the spinning projectile, were investigated. Additionally, the nonlinear behavior of these coefficients and coefficient derivatives were analyzed through the predicted flow fields. The considered Mach number ranges from 0.14 to 1.2, and the nondimensional spinning rate (PD/2V) is set to 0.186. To calculate the coefficient derivative of the corresponding force or moment, additional simulations were conducted at the angle of attack of 2.5 degrees. The simulation results were able to predict nonlinear behavior, the especially abrupt change of the predicted coefficients and derivatives at the transonic Mach number, 0.95. The simulation results, including the skin friction, pressure, and velocity field, allow the characterization of the nonlinear behavior of the aerodynamic coefficients, thus, enabling better predictions of projectile trajectories.


Introduction
Common methods to predict the aerodynamic coefficients of a spinning projectile in the design phase are to adopt the semiempirical method and conduct a wind tunnel test of the model. The representative software, based on the semiempirical method, is PRODAS [1] and DATCCOM [2], respectively, developed by Arrow Tech Associates Inc., the USAF. These methods very quickly obtain the aerodynamic coefficients with a given geometry and flow conditions, but when the geometry and flow conditions are out of the recommended data range, the prediction accuracy cannot be guaranteed. The wind tunnel test can be an alternative method to obtain the aerodynamic coefficients of the spinning projectile. However, this method still assumes that the tested model is similar to the potentially supersonic projectile and is expensive. Recently, as computing power has increased with a lower cost and the simulation algorithms have been more fully developed, computational fluid dynamics (CFD) has been increasingly adopted to predict the aerodynamic coefficients of the spinning projectile. The CFD methodology is able to accurately predict the aerodynamic coefficient, including static and dynamic load, and improve the performance through the analysis of the flow field around the projectile and its control surfaces.
DeSpirito and Heavey [3] conducted CFD simulations to predict the aerodynamic coefficients and flow field over a spinning 25 mm projectile. Various turbulence models such as the realizable k-ε, R, k-ε-R, and DES turbulence models were adopted, and the results were compared with experimental data and the PRODAS results. They found that a steady-state simulation with the traditional k-ε equations model is adequate for the prediction of aerodynamic force, pitching moment, and rolling damping, but an unsteady detached-eddy simulation is able to predict the Magnus moment more accurately. Additionally, DeSpirito [4] investigated the aerodynamic characteristics of a 155 mm projectile at a high angle of attack using the steady Reynolds-Averaged Navier-Stokes (RANS) and time-accurate RANS/Large Eddy Simulation (LES) methods. He found that RANS is adequate for predicting all aerodynamic coefficients below an angle of attack of 20°, except the Magnus coefficient and unsteady RANS/LES are required for the correct prediction of an angle of attack greater than 20 degrees.
The aerodynamic coefficient of a 0.50 caliber spinning projectile was studied numerically from the subsonic to supersonic region by Silton [5,6]. The CFD results of the smooth geometry without considering rifling showed good agreement with the static aerodynamic coefficients of the reference data [5]. However, the predicted derivative coefficients such as rolling damping and the Magnus moment are not well-matched, though CFD results were not worse than the aeroprediction code AP02. He emphasized that the spinning projectile from a rifled barrel should be considered for a better prediction of derivative coefficients [5]. Additionally, the nonlinearity in the Magnus force could not be correctly predicted in the steady or unsteady state CFD [6].
Simon et al. [7] simulated a 6 caliber projectile with a Secant-Ogive-Cylinder-BoatTail (SOCBT) geometry using the Spalart-Allmaras turbulence model. They showed that the simulated Magnus forces were in good agreement with the experimental data and investigated the main part between the ogive and boattail to reduce the Magnus force in the supersonic region.
In the present work, we conducted a numerical simulation on the geometry of a spinning projectile with a diameter of 120 mm and compared the predicted aerodynamic coefficients with the results of the semiempirical method, PRO-DAS. Specifically, six aerodynamic coefficients comprising the dominant parameters to solve the equation of motion of the spinning projectile are considered from the subsonic to supersonic region. Further, the mechanism of the nonlinear behaviors of the aerodynamic coefficients in the subsonic and transonic Mach number regions is analyzed through a detailed investigation of the flow fields.

Numerical Method and Simulation
Setup. The configuration in the present work is designed optimally to increase the maximum range of the projectile by the company. The computational model, which is shown in Figure 1, is about a 120 mm diameter (D) projectile, with a length of 4.9 D. The center of gravity is located at 2.9D from the nose tip. The rotating band is located at the position of 3.7 D from the nose tip with a width of 0.2D. Figure 2 shows the computational domain of circular conical type with a distance of 5 L upstream and 10 L downstream, where L is the length of the projectile. The computational domain is divided into two zones: one is a rotating zone around the projectile and the other is a stationary one outside the inner rotating zone. The unstructured grid is generated with the commercial software, ICEM CFD of ANSYS [8]. To ensure that the wall unit y + is less than 1.0 for the SST k-ω turbulence model, the first grid point off the wall is set to 2:0 × 10 −6 , as shown in Figure 3.
The ANSYS Fluent software V18.0 [9] is used to predict the flow around the spinning projectile from the 0.14 subsonic and 1.2 supersonic Mach numbers. The pressurebased solver and steady solution are set in the solver type, and the coupled scheme for pressure-velocity coupling is selected for the solution methods. The gradient term is computed by the least square cell-based method and the remaining terms in the momentum, energy, and turbulence equations are discretized in the second-order upwind scheme. To obtain faster and more robust solutions, the pseudo-transient solution method is used with implicit underrelaxation. The turbulent eddy viscosity is calculated from the SST k-ω turbulence model. The multiple reference frame (MRF) model, which considers the absolute velocity formulation and an extra source term related to the additional acceleration, is adopted to simulate the spinning object. Table 1 shows the Mach number and resulting roll rates with the nondimensional rotating velocity parameter PD/2 V = 0:186 in the present simulation. The analyzed Mach numbers range from 0.14 in the subsonic region to the supersonic value of 1.2, with a dense interval near M = 1:0. In general, most of the flight Mach numbers of the projectile are larger than 0.4, however, in the high angle fire, the projectile may cross into the condition of the low subsonic region near its maximum height. The considered angles of attack are 0.0 for the base drag and 2.5°to calculate the derivative of force and momentum. Other fluid properties, such as density and temperature, are set to the standard sea-level condition. The outer boundaries are set to the pressure far-field condition. The Reynolds number based on the velocity and the diameter of the projectile is from 3:9 × 10 5 to 3:3 × 10 6 according to the Mach number from 0.14 to 1.2 in the present simulation.

Force and Moment Coefficients.
In the present work, six dominant force and moment coefficients are calculated to solve the motion equations of the spinning projectile. The notation and definition of each coefficient are summarized in Table 2 and basic rules of direction and notation follow the reference of McCoy [10]. Forces are nondimensionalized using 1/2ρ ∞ V 2 ∞ S and S is the reference area as S = πD 2 /4. Similarly, moment coefficients are obtained by dividing 1/2 where P is the spinning rate in rps (revolutions per second) and PD/2V ∞ is the nondimensional spinning rate.  International Journal of Aerospace Engineering In the reference of McCoy [10], the drag coefficient is usually approximated by where C D 0 is the zero-yaw drag coefficient, C D δ 2 yaw drag coefficient, and δ = sin α t . The total yaw angle α t can be approximated with the angle of attack α and angle of sideslip . The difference between its approximation and exact definition may be insignificant when the total angle is less than 15° [10]. The lift coefficient can be expressed in the sum of linear and cubit coefficients with a term proportional to δ = sin α t and δ 3 . In the present work, the linear term will be considered through nondimensionalization by 1/2ρV 2 ∞ S sin α t . The Magnus force is produced by unequal pressures on opposite sides of a spinning object and is considered the important parameter for the prediction of the final draft range of a spinning projectile. The Magnus force coefficient derivative will be calculated based on the side force at two angles of attack (0°and 2.5°) through CFD. This parameter will be nondimensionalized by 1/2ρ ∞ V 2 ∞ SðPD/2V ∞ Þ sin α t . The overturning moment and spin damping moment are also dominant parameters in solving the motion equations of the projectile. The overturning moment is related to the lift force and can be referred to as the pitching moment. The spin damping moment, which tends to reduce the axial spin, should be negative, and the direction of the moment is x ! .
All coefficients which are considered in the present work are presented in Table 2 with the corresponding mathematical formulae. The last symbol in the nomenclature is α, which is described as the derivative of force or momentum coefficient and is calculated using the difference of the conditions of 2.5°and 0.0°and then divided by sin α t .

Validation.
The employed grid generation method and numerical scheme are validated through the numerical simulation of a 155 mm M107 projectile without the spin and the results are compared with the experimental results [10,11]. The geometry of M107 is shown in Figure 4. The total number of grids is 6:2 × 10 6 , and the first grid point off the wall is set to be less than y +~1 :0. The methodology for grid generation and the adopted numerical scheme are the same as the simulation for the modeling of a 120 mm projectile in the present work. In particular, two turbulence models, the SA    [12] and Shear Stress Transport (SST) k − ω turbulence model [13], are tested to ascertain the optimal turbulence model for this study. Figure 5 shows a plot of the drag coefficient with respect to the Mach number from 0.6 to 1.5 at an angle of attack of zero. The line corresponds to the experimental result, and the two additional symbols represent the results of the SA and SST k − ω turbulence models. The drag coefficient is flat in the subsonic region and then increases abruptly in the transonic region until M = 1:1 and then decreases gradually with higher Mach numbers. Since we are focusing on the flight with a Mach number of 1.2, the simulation is carried out until M = 1:5. The results of the SST k − ω model agree well with experimental data [10], but the SA model overpredicts the drag coefficient across the entire Mach number region. To compare the results of two turbulence models, pressure drag and viscous drag are extracted in the total drag at the Mach number of 1.2. The respective pressure drag coefficients are 0.3659 and 0.3873 in the SST k-ω model and SA model, with a 5% difference. However, the viscous drag coefficients are 0.0096 in the SST k − ω model and 0.0545 in the SA model, which has a one-order difference. The overprediction of the viscous drag coefficient in the SA model is the cause of the overprediction of the total drag coefficient of 0.4418 in the SA model, whereas the total drag coefficient is predicted to be 0.3791 in the SST k − ω model. This phenomenon is well known and related to the fact that the viscous boundary layer is poorly revolved when y + is less than 1.0 in the SA model. Therefore, we will employ the SST k − ω turbulence model to simulate the 120 mm spinning projectile in the present work.
2.2. Drag Force Coefficient. Figures 6(a) and 6(b) show the coefficients of zero yaw drag and quadratic yaw drag, which are compared with the results of the semiempirical method, PRODAS. The predicted zero yaw drag coefficient agrees well with reference data until a Mach number of 1.0 and is slightly overpredicted at higher Mach numbers. This tendency that the CFD simulations overpredict the drag coefficient for supersonic Mach numbers is consistent with the spinning projectile simulation results by DeSpirito and Heavey [3]. On the other hand, there is a discrepancy in the prediction of the quadratic yaw drag coefficient. In such situations, this coefficient is calculated by multiplying the square of δ for the total drag, its order is approximated as 10 -2 and its effect may not be very significant.
The total drag coefficients, comprised of the dominant pressure drag in red and skin friction drag in blue, are plotted from the calculation of each component at five different Mach number conditions in Figure 7. The numbers in parentheses correspond to the percentage of the total drag coefficient for each drag term. As the Mach number increases, the total drag coefficient increases after the transonic Mach number and the percentage by pressure drag also increases, whereas the skin friction drag coefficient decreases. The abrupt increase after the transonic Mach number is due to the shock wave around the projectile. Figure 8 shows the contours of the Mach number in the xy plane at the conditions of six different projectile Mach numbers. As expected, the shock wave begins to generate Table 2: Notation and formula of six aerodynamic coefficients.

Coefficient
Name Formula Magnus force coefficient derivative     International Journal of Aerospace Engineering from Mach number 0.9 at the position of a different radius of curvature at the ogive and near the rotating band. As the Mach number increases, an additional shock wave near the front of the nose is shown. The wave drag by the shock wave causes the drag coefficient to increase abruptly. A large difference of the results between the two methods is observed in two regions; the first is the low Mach number (0.1-0.6) region and is observed above the second transonic 0.95 Mach number region. Henceforth, other coefficients such as lift and overturning moment coefficients show a similar behavior, which means that there is a difference in the results between two methods in these two Mach number regions. The semiempirical method, PRODAS, predicts the nearly constant values of the lift and overturning moment coefficients before the subsonic Mach number, 0.4, whereas the CFD results show an abrupt change in this region. Another interesting point is that the CFD results show a sudden decrease and then increase of the quadratic yaw drag coefficient around the transonic Mach number region whereas PRODAS predicts a gradual increase of this coefficient after Mach 0.8. Figure 9 investigates the difference between the quadratic yaw drag coefficients at Mach numbers 0.14 and 0.6. In this figure, the skin friction contours at an angle of attack of zero and the skin friction differences (C f α=2:5 − C f α=0:0 ) at two angles of attack (0°and 2.5°), which is calculated from the subtraction of the skin friction at an angle of zero from that at 2.5 degrees, are plotted at Mach numbers 0.14, 0.6, 0.95, and 1.2. At the constant condition of a nondimensional spin rate (PD/2V = 0:186), the spin rate increases with the 7 International Journal of Aerospace Engineering projectile speed, which means the rotating effect is less dominant at a lower Mach number. Therefore, the skin friction at the zero angle of attack can be accurately predicted by its lower Mach number before, and the difference, shown in the right-hand column, is also estimated to have a higher value, in particular near the ogive. This pattern seems to increase the quadratic yaw drag coefficients at increasingly low Mach numbers.
At the transonic Mach number of 0.95, the skin friction difference shows a negative value region after the rotating band, which causes an abrupt decrease of the quadratic yaw drag coefficient. Figure 6(c) shows the normal force coefficient derivative, CLα of the two methods. PRODAS predicts that this coefficient is constant across low subsonic Mach numbers and then gradually increases after Mach number 0.6. However, the CFD results show an abrupt decrease at Mach number 0.95 and then another increase, which is a behavior similar to that of the quadratic yaw drag coefficient. The simulation results by Silton [5,6] and experimental data [10] show the same tendency with the lowest point reached shortly before Mach 1.0. Also, it was mentioned in [1] that the semiempirical Aeroprediction code AP02 [14] was not able to capture this peak in the transonic regime. Figure 10 shows the contours of the difference of pressure (P α=2:5 − P α=0:0 ) between two angles of attack for the calculation of the normal force coefficient derivative to investigate why the corresponding value (CLα) changes abruptly. The left figures correspond to the upper surface of the projectile and right ones provide a back view of the projectile base. The contours at four Mach numbers (0. 14  This phenomenon is also related to the abrupt increase of the overturning moment coefficient in the positive direction. A detailed explanation will be mentioned in the next section. The pattern of pressure coefficients at the base at Mach number 0.95 has different characteristics with a lower overall value than that at the other three Mach number cases. The pressure difference coefficient at the base region is related with the quadratic yaw drag in the axial direction and the lower level of pressure difference will be another reason for its abrupt decrease at the specified Mach number of 0.95.  Figure 11. The blue line corresponds to the pressure coefficient difference, and the red line is the zero-pressure coefficient difference for a better understanding of the asymmetric pattern. At the x/D = 2:25 position, the pressure difference at the upper part of the projectile is negative and the value at the lower part is positive. However, a symmetric pattern with respect to the vertical y-axis, which corresponds to the connected line at 90°and 270°in the figure, is shown. This is not able to produce the nonzero Magnus force. At the x/D = 4:16 position, even though the pressure coefficient differences appear close to zero at Mach numbers 0.65 and 1.2, they are not zero. This is due to the large scale interval of the pressure difference coefficient from -0.2 to 0.2. However, the large asymmetric behavior with respect to the vertical y-axis can be shown at x/D = 4:16 and at the 0.95 Mach number condition. The pressure difference at the right-upper part is positive and that at the left-lower part is negative. The integration of two differences results in a value with a negative z-direction. This big difference seems to produce the large Magnus force derivative at Mach 0.95.

Overturning Moment Coefficient
Derivative. The values of the overturning moment coefficient derivative, CMα, which is referred to as the pitching moment coefficient derivative in other documents, are compared in Figure 6(e). There is a peak near the transonic regime which then decreases outside this regime in both methods. This phenomenon is called the critical behavior in the transonic regime [5]. Overall, the CFD results underpredict the coefficient derivative relative to the prediction of PRODAS, except the peak at the Mach  International Journal of Aerospace Engineering number of 0.95. In the subsonic regime, PRODAS predicts the overturning moment coefficient as a constant value, which means that the pitching moment varies linearly with respect to Mach number. However, the CFD results show the increase of CMα with Mach number, and the pitching moment is no longer linear. Other simulations and experimental results [3,5] show the decrease of the overturning moment coefficient derivative as the Mach number decreases in the subsonic regime. The abrupt increase of the overturning moment coefficient can be explained with the different pressure distribution near the boattail at the corresponding Mach number condition, which is shown in Figure 10. The positive pressure difference at the upper part of the boattail and the reverse behavior at the lower part acts to abnormally increase the overturning moment.
2.6. Spin Damping Moment Coefficient. The rolling damping moment coefficients are predicted with big differences between the two methods even though these are predicted as negative values to reduce the axial spin. The semiempirical PRODAS method predicts this coefficient to range from -0.027 to -0.025 with approximately a 10% difference, whereas the CFD results show the value ranging from -0.016 to -0.010 with about a 50% difference. When it is considered that other research by Silton [5,6] and DeSpirito and Heavy [3] predict a 30-50% difference within the Mach number regime of the present work and increasing rate as the Mach number increases, the CFD results seem to be more reasonable than those of PRODAS. As the Mach number increases, the effect resulting in a reduction of the rolling moment decreases due to the increasing inertia by the high freestream velocity even though the spin rate increases based on the constant nondimensional spin rate. At a low Mach number of 0.14, this coefficient decreases abruptly. This phenomenon is consistent with the experiments and simulation results of other researchers [3,5]. Figure 12 shows the contours of each direction of velocity (u, v, and w) at three different positions (x/D = 0:75, 2.25, and 4.42) from the nose tip. For comparison with the same scale factor, the u-velocity in the x-direction is nondimensionalized using the freestream velocity and the other two velocities are nondimensionalized using PD/2. The low-velocity region in u-velocity contour can be shown near the boattail (Section 3) at a Mach number of 0.95, which is a different flow pattern when compared with the results at other Mach number conditions. This different flow field can be shown in vand w-velocity contours. Additionally, in every Mach number condition, the axisymmetric behavior of the predicted velocities is especially large near the boattail after the band and the flow field in this region has a very significant effect on the motion of the spinning projectile.

Conclusions
In the present work, numerical simulations on the geometry of a spinning projectile with a diameter of 120 mm are conducted, and the predicted aerodynamic coefficients are compared with the results of the semiempirical method, PRODAS. The six dominant aerodynamic coefficients needed to solve the equations of motion of the spinning projectile are investigated, and the mechanism of the nonlinear behavior of aerodynamic coefficients is analyzed through a detailed investigation of the flow fields.
The simulation results were able to predict nonlinear behavior, especially the abrupt decrease or increase of quadratic yaw drag, lift coefficient derivative, Magnum force coefficient derivative, and overturning moment coefficient derivative at a Mach number of 0.95. Additionally, different predictions between CFD simulations and the semiempirical method, PRODAS, are shown in the low subsonic region.
The increase of the quadratic yaw drag at a low Mach number of 0.14 and the abrupt decrease of this parameter at the transonic Mach number of 0.95 can be explained through the distribution of the skin friction coefficient around the body of the projectile.
The difference of the pressure coefficient between two angles of attack (0.0°and 2.5°) is presented to investigate the nonlinear behavior of the lift force coefficient derivative. The downward pressure difference distribution near the boattail at the transonic Mach number of 0.95 decreased the lift force coefficient derivative. Further, the axisymmetric distribution of the pressure difference near the boattail caused an abrupt change of the Magnus force coefficient derivatives.
The CFD results show the increase of the overturning moment coefficient derivative with Mach number and predicted a nonlinear pitching moment before the transonic Mach number region, whereas the semiempirical method, PRODAS, predicted a constant value of the derivative.
The contours of each direction of velocity (u, v, and w), which are nondimensionalized by the freestream velocity or PD/2, showed an obvious difference at a transonic Mach number of 0.95 when compared with the results at other Mach number conditions.
The predicted aerodynamic coefficients based on the CFD simulation results will be adapted to the coefficient matrix in equations of motion and increase the accuracy of the predicted trajectory of the spinning projectile.

Data Availability
"Data Availability" statement for the data in the present paper is corresponding to \"Data available on request\".

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this paper.