Aerodynamic Coefficients Prediction for 122 mm Rocket by Using Computational Fluid Dynamics

In this study, the Computational Fluid Dynamics code in SolidWorks was used to predict the aerodynamic coefficients of 122 mm rocket in order to get more accurate aerodynamic coefficients and derivatives for flight simulation. This paper presented the governing equations of fluid dynamics applied in SolidWorks and developed the computational model of the rocket. The results were compared against calculated data from an empirical method. The conducted analyses showed that the comparison results are nearly the same. However, the empirical method provided a faster way to obtain results.


Introduction
This study focuses on the aerodynamics calculation and analysis of small rockets. The parameters of airflow characteristics, such as the drag force, pressure, and airflow velocity, have significant influence on the exterior ballistics characteristics of the rocket, and thus affect the effect and accuracy of the rocket. The essential thrust characteristics of the rocket engine are directly related to the drag coefficient and the drag force, which are the foremost aerodynamic parameters in the study of exterior ballistics. Computational fluid dynamics (CFD) is a currently pervasive application in aviation for evaluating aerodynamic performance during the conceptual stage and preliminary design stage in order to reduce design cycle time and minimize the experimental validation associated costs [1]. The purposes of this paper are to compare the aerodynamic coefficients of the static stability coefficient and the axial force coefficient of 122 mm rocket obtained by SolidWorks and an empirical method and to get more accurate aerodynamic coefficients and derivatives for flight simulation.

Axial Force Coefficient Prediction
Drag has two fundamental sources of pressure distribution and surface friction around the rocket. Pressure distribution is further divided into base drag, parasitic pressure drag attributed to prominences such as emitter lug and body pressure drag including supersonic shock wave. The total axial force applied to the rocket can be described by: Where A C is the axial force coefficient, V is the absolute velocity, and ref A is the reference area. 1.5ln 5.6 (10 ) Where s R represents the surface roughness in m  , and ref L represents the reference length in meter.
The turbulent skin friction drag coefficient is corrected for compressibility with: [2, 3] Where M represents the Mach number. The skin friction drag coefficient based on the reference area can be found by [2, 3] Where w A represents the wetted area of the rocket.
The wetted area of fins is twice its planform area, the skin friction drag coefficient of the tail is: [2, 3] 2 The subscript ⊥ symbolizes that the leading edge is perpendicular to the flow.
The pressure drag at the leading edge of a slanted fin is derived from the principle of cross flow as [2] ( ) ( ) Where L  represents the leading edge angle. Note that both coefficients in the equation are relative to the face area of the cylinder, so their reference area ratio is still cos L  . The fin base drag coefficient of the square profile fin is equal to the body base drag coefficient: The total fin pressure drag represents the sum of the trailing edge drag and the leading edge drag: 2.1.3. Body Pressure Drag. The pressure drag is generated by the air around the rocket. This paper describes the estimation of the pressure drag of the nose cone, ignoring other body pressure resistance. The pressure drag of the streamlined nose cone is notably less than the skin friction drag at subsonic speeds. Figure 1 presents the measured pressure drag coefficients for different nose cone shapes. (13) In this paper, a simple semi-empirical method is used to estimate the transonic pressure drag of the nose cone, rather than the general segmentation calculation. This paper adopts a semi-empirical method under the circumstance of Mach number over 0.8. The pressure drag is interpolated between the discreet value of Equation (13) and the transonic method at high subsonic speeds. Assume that the pressure drag does not decrease in the subsonic zone, and the derivative is zero at 0 M = . An interpolation function similar to the Prandtl factor is proposed as: (14) Where a and b are calculated for the fitting analysis of the drag coefficient as well as its derivative at the lower limit of the transonic method. Equation (15) is to calculate the pressure drag of the conical nose cone at supersonic speeds (M > 1.3). The drag between Mach 1 and 1.3 is interpolated with polynomials using the boundary conditions: Where the drag coefficient is approximately  at 1 M = and =1.4  is the air specific heat ratio. At supersonic and transonic speeds, the pressure drag of ogival nose cones with identical length with shape factor correction and diameter is equal to that of conical nose cones as: The shape factor is 1 at 0.1  = and 0.82 at 0.5

 =
. For ogival nose cone, the parameter  is one.

Base Drag.
Base drag is generated by anywhere where the body diameter rapidly decreases or the low pressure area at the bottom of the rocket. An empirical formula can be used to estimate the magnitude of the base resistance as: A reasonable approximation is obtained by reducing the propulsion motor area from the reference region [2]. Therefore, the base drag will not be generated with the same base size. In addition, with a large base and a single small motor at the center, the basic drag is roughly the same as when sliding.

Total Drag.
The pressure drag is generated by the air around the rocket. This paper describes the estimation of the pressure drag of the nose cone, ignoring other body pressure resistance.

Static Stability Coefficient Prediction
A rocket is regarded to be: [4] • Statically unstable if the disturbance will keep the rocket away from equilibrium.
• Statically stable if small perturbations tend to bring the rocket back to equilibrium. Static stability is demonstrated in figure 2. The center of gravity (CG) is the point where the rocket rotates about when disturbed. The center of pressure (CP) is the point where synthetic aerodynamic forces are assumed to act. Note that CG and CP may not be on the longitudinal centerline of the rocket. In this case, the effect of axial force on stability should be considered. In short, a rocket can be considered statically stable when CP is behind CG, both measured from the rocket tip.
• The lift on the rocket body tube is negligible.
• The viscous forces can be neglected.
• The nose of the rocket is poised at a point.
• The rocket is relatively narrow in diameter.
• Steady air flow over the rocket without rapid changes. The normal force on the rocket can be calculated by: Where, N C is the normal force coefficient.
At very small attack angles, N C can be approximated as linear as  . Hence, it can be expressed as: Where  represents the attack angle in radian and N C  represents the slope of the normal force coefficient, also known as derivative of stability. The total derivative of stability can be derived from the sum of the derivatives of stability for each part.
Where R represents the whole rocket and P represents a rocket part.
The CP position of the rocket is calculated according to CP of every independent rocket component. Similar to the center of mass (CM), CP locates on the rocket's rolling axis, the CP location is defined as the distance ( CP X ) from CP at the tip of the nose cone. CP X can be calculated using Equation (25).
Where ( ) CP P X is the distance from CP of component P to the tip of the nose cone of the rocket.
Equations of N C  and CP X for conical variations in body diameter, fins and nose cones are as follows.
As long as the nose cone is parabolic, conical or ogive, the stability derivative of the nose cone is two.
The position of the center of pressure for the nose cone depends on the shape: The term m K is a Mach correction factor that is empirically derived from the body alone data: The margin of stability can be defined as: Where CG X is the center of gravity from the nose cone tip.

Background Theory of SolidWorks
SolidWorks Flow Simulation solves the Navier-Stokes equations.
[7] Continuity equation: Where u represents the vector of velocity and  represents the density.
Momentum equation: Where P is the pressure, S is the momentum source,  is the stress tensor.
Energy equation: Where H is the total enthalpy, related to the static enthalpy h by: The following energy equation is adopted for high-speed compressible flows and shock wave flows as:

Modeling and Mesh Generation of 122 mm Rocket
Note that the optimal size of the control volume is 5 to 10 times the length of model at the side and back and at least the length of model in the front [8]. In the paper, the control volume is 5 times the model length at the side and back, and one length in the front of the rocket. Geometric characteristics of the 122 mm rocket is shown in figure 4.  In the flow analysis, numbers of divisions are chosen 85 along X-axis and 75 along Y, Z-axis for initial mesh generation in Global Mesh setting. Local Mesh is added for the entire rocket to increase grid density. Grid generation on the 122 mm rocket is shown in figure 6.

Result and Discussions
Simulation data were obtained according to various Mach numbers ranging from M= 0.2 to M=2.2 with a step of 0.2 at 0° and 2° angles of attack. The computer to run SolidWorks Flow Simulation adopted a Quad-core Intel® Core™ i7-6500U processor running at 2.5 GHz with 12GB of RAM. The computer ran the 64-bit version of SolidWorks 2017. The average running time of one data point was approximately 1.5 hours. The airflow velocity and pressure distribution contours on the cross-section plane, when the Mach number is one at 0° angles of attack are shown in figure 7 and figure 8 respectively. High pressure regions occur on the rocket nose and fin leading edges due to air wave. Low pressure region occurs at the rocket base due to airflow separation from the rearward-facing step. These can be seen in figure 7. Due to low pressure condition, the air flow velocity is decreased behind the rocket base of the low pressure region. It can be seen in figure 8.   (47) and (48). In this paper, 0° and 2° are used for the computation for The comparison between the empirical and CFD computed axial force coefficient shows a difference of about 31 percent in the subsonic range, good agreement in the transonic region, and a difference increases in start point of the supersonic region, but a difference is gradually decreasing for increasing Mach number, as was shown in table 1 and figure 10.  The comparison of the normal force coefficient and normal force coefficient derivative calculated with CFD and Empirical method shows the good agreement in subsonic and transonic region, but a different of about 13% and 15% at M=0.8 and M=1.1, respectively. The results are shown in table 2  and table 3. Moreover, a difference in the supersonic region is a little more than in other regions. These are illustrated in figure 11 and figure 12.

Conclusions
The conducted analyses have shown that the comparison results are nearly the same. For the axial force coefficient, the differences are apparent in the subsonic region. However, for static stability AAME 2020 IOP Conf. Series: Materials Science and Engineering 816 (2020) 012010 IOP Publishing doi:10.1088/1757-899X/816/1/012010 10 coefficients and derivatives, the differences are noticeable in the supersonic region. For all coefficients, the comparison results show a good agreement in the transonic region. Generally, although more accurate results can be achieved by CFD method, the empirical method has provided a very fast way to obtain results. Based on the numerical solution of the Navier-Stokes equations, the CFD code took about 1.5 hours for one data point.