Improved Estimation of Air-to-Air Missile Aerodynamic Characteristics with Using CFD Method

: The paper describes a solving procedure of the special Finite Volume Method (FVM) problem dealing with calculation of a flow field around air-to-air missile at subsonic, transonic and supersonic velocity range in order to obtain a coefficient of drag (c D ). Some emphasis is put on numerical solution and fundamentals of FVM and use of a k- ω turbulence model for simulation of subsonic and supersonic flow fieldsaround air-to-air missile AIM-9 Sidewinder version M. Obtained results of missile aerodynamics characteristics are verified using semi-empirical software and by analysis of flow fields.


Introduction
Computational Fluid Dynamics (CFD) is already a classical branch of a contemporary fluid dynamics. In the past few decades, the CFD used personal computers or work stations with appropriate software based on data structures and numerical analysis to solve problems involving interaction between fluid and solid object in defined integration domain.
CFD based on the Finite Volume Method (FVM) is also increasingly applied inside development cycle of new generations of air-to-air missiles. The growing importance of FVM reflects increased requirements for higher missile velocity, higher maneuver ability or other tactical demands. Another factor in rising role of FVM is improved engineering effectivity in scope of design and development processes of missiles systems, as well as their analysis. When properly applied, the discipline of CFD should play an alternative role in replacement of some empirical and experi-mental approaches in solving complex problems of flow and heat transfer domains. In any case, the FVM is a low-cost tool in comparison to classical flight or wind tunnel tests in aviation or aeroballistics. This paper will describe the procedure for solving the special FVM problem which deals with calculation of flow field around air-to-air missile at subsonic, transonic and supersonic speed range in order to obtain Coefficient of Drag (cD). cD is the most required aerodynamic coefficient, which is required in missile performance calculations. In a limited scope of this article, we will deal with a brief description of the fundamentals of FVM. Some emphasis will be given on the explanation of k-ω turbulence model used to solve this problem.
In order to solve any FVM problem, it will be necessary to create: • the 3-D body of examined object in suitable Computer Aided Drawing (CAD) tool, • the integration volume, in which the fluid medial will act on to surface of examined object.

Governing Equations of Fluid Flow
Missile aerodynamics, from subsonic to hypersonic velocity, is governed by fundamental equations of fluid dynamics; see below the list of governing equations for compressible and unsteady flow of Newtonian fluid. These equations are mathematical statements for conservation of mass, momentum and energy. The basic set together with the equations of state relating pressure, density and temperature of the fluid, and also with equations described the properties of the fluid and formulas for numbers of similarity [1][2][3][4][5][6], provide the governing equations of motion of volume of fluid continuum. Continuity equation (mass conservation law): Navier-Stokes equation (momentum conservation law) for x-component of velocity: Navier-Stokes equation for y-component of velocity: Navier-Stokes equation for z-component of velocity: where p -the pressure; µ -the dynamic viscosity of flowing fluid; fox and latter foy and foz represent external volumetric force effects, e.g. in form of centripetal acceleration or acceleration due to gravity. Fourier-Kirchhoff equation (energy conservation law): where cp -the heat capacity at constant pressure; T -the temperature in Kelvin scale; FT denotes so-called equation of thermal state -special function developed for every fluid. λ -the thermal conductivity of liquid, and qg -the specific intensity of heat source. State equation of ideal gas: Function of the specific heat capacity: Function of the dynamic viscosity: Function of the thermal conductivity: ( ) The introduced set of partial differential and algebraic equations gives us a good insight to the measure of its complexity. Direct solution of this set of equations for a general flow problem is not possible because of its nonlinear nature. In the past, engineers made further approximations and simplifications to the equation set until they had a set of equations that they could solve. Current high-speed computers are used to solve appropriate modified governing equations (Eqs (1)-(9)) using a variety of techniques grouped to FVM. This is current scientific and industrial standard which is realized in various commercial software packages, e.g. Ansys ® , Comsol-Multyphysics ® , STAR-CCM+ ® , etc.

Numerical Solution
In order to find a solution of the flow field variables, the pressure, temperature and the velocity of gas, we have to discretize the fluid flow volume. The most of CFD software uses numerical methods of approximate solution of the partial differential equations applied inside the discretized volume -so-called Finite Volume Method (FVM). The first mention of this method appeared in Professor A. Jameson's presentation in 1977 [7].
The appropriate modification of the governing equations mentioned above involves their discretization in space and time. The flow of the physical quantity is given by the integral sum over the four boundaries of the control volume in case of twodimensional tasks, or six areas in case of three-dimensional tasks [8].
In the solution, the whole considered volume is divided into a finite number of small control volumes using a defined computing network also known as a mesh. The computational mesh is composed from polyhedrons and therefore it is possible to define boundary nodes, edges and faces for each of these polyhedrons. Computational node or control volume is also called a cell. The individual control volumes do not overlap each other. Since the 1990s, it has been possible to use unstructured computing meshes that can be created without Cartesian indexing. Control volumes of these meshes have no shape constraints, they provide considerably lower computational memory requirements, as well as they increase computational speed when dealing with more demanding tasks [8].
Values of fluid flow scalar and vector quantities are determined in geometric centers of control volumes. Different interpolation types are then used for the boundary values determination. The first type of interpolation is countercurrent 1 st order interpolation. This method assumes that the value of the physical quantity flow is equal to the value in the center of the cell lying upstream. For countercurrent 2 nd order interpolation, the value is determined from two upstream cell centers. The last possibility of the values interpolation at the control volume limits at a given setting is the central difference use. In this case, the flow value on the wall is determined by linear interpolation between the values at the centers of adjacent cells [8].
In most engineering applications of turbulent flow, where the missile flight simulation belongs, turbulence models based on the Reynolds Averaged Navier-Stokes Equations (RANS) method are used. This method uses statistical averaging methods in order to solve and simplify basic equations. The flow is assumed stable, so we can create the average value of a given quantity in different time periods from different time records and always get the same value. This phenomenon was first noticed by Osborn Reynolds who suggested decomposition of turbulent stream into the mean value and fluctuation in time. The main problem in the simulation of turbulent flow is the presence of Reynolds stress in the centered Navier-Stokes equations, which, unlike the laminar flow, leads to an unclosed system of mathematical equations. Therefore, additional equations and empirical relationships that form a model of turbulence are added to the equations of motion [9].

K-Omega Turbulence Model
The basis of the whole group of turbulent models is the so-called Boussinesque hypothesis. This hypothesis enables to define the shear stress using 9 components (three equations in three spatial directions), which will replace only one variable named turbulent viscosity. The turbulent viscosity does not depend on the flowing substance However, it depends on the properties of the flow mode only. It is equal to zero in laminar flow [8].
Turbulence models, which are based on the Boussinesque hypothesis of turbulent viscosity, solve the turbulent viscosity value using additional equations. Depending on the number of differential equations used to define the turbulent viscosities, the model name is used [10].
CFD simulations in this work were carried out using a two-equation model. Twoequation models solve differential equations for turbulent kinetic energy and differential equations for length measure of turbulence, which depends on transport and process history, as well as on kinetic energy. The swirls are carried downstream and their size at each location depends on their original size. A particular k-ω model is used in this article which uses a quantity called turbulence to describe the length scale of turbulence. This quantity expresses the rate of liquid rotation at a given point of the flow field. The turbulent viscosity in this model is described by the relation: where vt -the turbulent viscosity, k -the turbulent kinetic energy, ω -the vorticity.

Role of FVM in Missile Design and Investigation
FVM has many applications during the process of design and development of missiles. We can use it as a significant help at: • the determination of aero-thermodynamic load acting on missile airframe (i.e. temperature and pressure fields or heat flow fields) in order to support structural analysis, • the determination of aerodynamic coefficients (i.e. coefficient of drag, lift, pitching moment, …) as irreplaceable inputs o into the missile flight performance, stability and maneuverability analysis, o into the advanced 6-degree of freedom flight simulations, • the determination of thrust characteristics by aerothermal analysis of hot gases flow field inside of the rocket engine, • the determination of airframe interaction between the missile body and launcher or other parts of airborne carrier, • the preparation and tuning of preliminary wind tunnel testing or final flight testing.
We need these aerodynamic data very extensively, because air loads vary with Mach and Reynolds numbers, altitude, angle of attack, angle of sideslip, aerodynamic roll angle, accelerations, angular rates and accelerations, control surface deflections and thrust operation. It is clear that there are many variables involving aerodynamic forces and moments generation. Classical approach involves the building of a large database based on the wind tunnel measurements, which provides the closest results for many flow problems. However, this process is very expensive and time consuming.
FVM gives us a new added value in the phase of preliminary estimation of expected aerodynamic coefficients, pressure and temperature distribution as a result of the case study, as well as the supporting information to planning wind tunnel experiment, i.e. sensor selection and placement, estimation of test chamber walls effect, flow asymmetry, etc. Suitable application of FVM can even assist us in the design of entire test matrix. Also flow field visualization can offer us a valuable insight to the various aerodynamic phenomena and it can support wind tunnel test planning.
It seems that FVM will not replace a well-planned and well-executed wind tunnel test series, but on the other hand, the impact of FVM into the wind tunnel test process can save development costs and it also eliminates an unexpected cost of repeating poorly planned wind tunnel tests as well as flight tests.

Preparation of Canard Missile 3D Geometry into the CFD Processing
To achieve valuable results demonstrating CFD ability we will focus on the frequent aeroballistics task -the determination of the dependency of drag coefficient on Mach number, cD(M). With the help of this isolated aerodynamic characteristics, we can perform a simple but very useful estimation of missile trajectory and its elements which is the most frequent task in exterior ballistics. In order to calculate the trajectory, we need to estimate cD (M) dependency in a large range of velocities. The highest value, approx. M = 2.5-3, represents the missile velocity at the end of the powered phase of the missile flight. The lowest value, approx. M = 0.75-0.9, denotes a high subsonic velocity regime where the missile usually will lose the ability to self-control and it is also going to be without internal source of thermal battery power. For the sake of simplicity, we can assume a stable flight without significant maneuver where the angle of attack and sideslip angle will be negligible − approximately equal to zero. Practically, there is only one significant variable, which is the Mach number. Thus, our task still covers three characteristics of flow phenomena -high subsonic, transonic and supersonic flow regimes. This causes our task being sensitive to correct tuning of the effect of turbulence which not simple in such a wide range of velocities.

Representative Missile Geometry Specification
We have chosen the AIM-9 Sidewinder as a typical and quite simple representative of the state of the art of heat-seeking, low-range, canard missile design. In order to calculate aerodynamic characteristics, a simplified model of the simulated body had to be created in the parametric CAD (Computer-Aided Design) modeling software at first. Geometry of the AIM-9M Sidewinder had to be obtained from technical drawings available on the internet [11], see the example in Fig. 1.

Fig. 1 Schematic diagram of AIM-9 Sidewinder with dimensions
Due to the incomplete documentation, some dimensions had to be estimated based on photos, see an example in Fig. 2. The resulting virtual missile assembly (Fig. 3) consisted of the body of the missile, four control surfaces and four wings.

Integration Volume
After the body import into the STAR-CCM +, the CFD simulation software, a one-half of cylindrical volume around the missile was created. This useful simplification was accepted because of the symmetry of the missile body; in addition, the axial flows are also symmetrical. The axis of this half-cylinder is identical with longitudinal axis of the missile. The resulting three-dimensional integration volume was created using the "Boolean" function by subtracting the volume of the guided missile from the cylinder volume. The integration volume was subsequently set in this volumes difference, see Fig. 4. The symmetry plane boundary condition was set for the longitudinal flat wall of the integration volume while the inlet and the outer cylindrical surface was set as the Free Stream and the remaining cylinder bottom as the Pressure Outlet.

Mesh Creation
The start of the mesh creation in the integration volume was possible after the boundary conditions setting. Many options have been tried, but finally the Trimmer mesh was chosen for evaluations. In this type of mesh, the control cells are cubic. The increase in cell size is achieved by joining eight cubes of the previous size into one. The number of cells of the same size then gives the cell growth gradient. The higher the number of cells of the same size, the smaller the growth gradient. The Trimmer computing mesh provides a robust computing mesh solution around simple and complex bodies. It combines the advantages of hexagonal mesh and cells with optimum skewness. The edges of the cells can be arbitrarily oriented according to the selected coordinate system. The Trimmer computing mesh is a very efficient tool for copying small radius of a flowed body and at the same time, it very quickly reduces the number of cells where there is a flat surface and increasing distance from the flowing body. Trimmer computational mesh, see Fig. 5, had an integration volume with external borders in the shape of half a cylinder. The settings included adjustments to the size of the control volume on the control surfaces, wings and missile body. Furthermore, five structured layers of the integration mesh were set in the boundary layer. In addition, two areas of smaller size control cells were added at the predicted shock wave locations around the missile nose, control surfaces, and wings. Next, a cylindrical area was added with finer cells in the wake area behind the missile. This integration mesh con-

Calculation of Missile Drag Coefficient
The aerodynamic characteristics were obtained continuously in the range of Mach number (flight velocity / local speed of sound) 0.4-3.2, which is a bit wider than necessary. The first case, M = 0.4, was calculated and the result was used as an initial estimation for cases with higher M. This solution scheme (previous result = initial estimation for higher M) was applied to the entire matrix of Mach numbers. The Mach number was used as a boundary condition at the inlet of the integration field. From this calculation, we obtained longitudinal aerodynamic force (XA) for particular M. At the zero angle of the attack, we can take this force equal to drag force (D = XA). The primary objective of our effort, coefficient of drag (cD), was calculated by the following equation [12].
where S -the cross-section area based on fuselage diameter as a referential dimension. Other forces like a lift (L) and a side force (Y) were equal to zero due to the symmetry of the missile body and flow. As a primary result of our CFD calculation, the missile drag curve was obtained, see Fig. 6. In the subsonic area, we can see an increase of the drag coefficient. Latter, sharp increase is clearly visible as it started from approximately M = 0.85. The increase of value of cD is caused by the growth and formation of shock waves around missile body. Drag coefficient rises continuously in transonic region even when the speed of sound is exceeded. The maximum of the drag coefficient corresponds to a velocity given by approx. M = 1.05. The calculated relationship corresponds to a typical shape of the drag coefficient curve, with the peak in transonic region at M ≥ 1. In a very small scale at those velocity range, we can see a little oscillation in cD(M) distribution. This local minimum and regrowth can be caused by complex shock waves formations from various parts of the missile or rough mesh scheme at the nose part of missile body. To determine the unequivocal cause, it would be useful to perform further calculations with different settings of mesh or different choice of model of turbulence. Ideally, it would be appropriate to compare the calculated curve of the drag coefficient with experimentally measured data.

Validation of Results
The problem of validation leads to an appropriate form of some experiment. Today, the technologies of the wind tunnel testing allow quite accurate measurements of the missile drag curve in the desired velocity range. However, the obvious price is proportional to the required accuracy. In any case, we cannot expect anything simple in low price. It is out of the scope of our investigation.
The next way of experimental verification is the flight test. We insisted that Raytheon Company -the manufacturer of missile, perform series of flight tests with results in form of aerodynamic characteristics. Those results are secret and therefore unpublished because of expectable defense industry restrictions.
Web database does not contain any hint of experimental data available for our objectives. Some sources refer to the result of calculations, but unfortunately with none or insufficient description of conditions. In this case, we can use some result as an imperfect but reliable CFD counterpart -calculation with help of semi-empirical methodology. We used freeware (www.rasaero.com) [12] to generate approximate missile geometry and perform calculation of cD(M) distribution at the example of AIM-9 Sidewinder. For the geometry of the missile, see Fig. 7.
The final result calculated by RASAERO freeware is shown in Fig. 8. The results show us a shift to higher (approx. 30%) values of cD(M) curve in subsonic flow regime and slightly higher (approx. 10%) in transonic regime as obtained by semi-empirical calculation. The semi-empirical calculation also indicates a tougher drop of this characteristic in supersonic flow regime. Fig. 7 Geometry of Missile AIM-9 used in RASAERO freeware [12] Advances in Military Technology, 2020, vol. 15, no. 2, pp. 301-315 311 Fig. 8 Drag coefficient curve of guided missile AIM-9 Sidewinder by RASAERO freeware [12]

Flow Field Calculation
The ability to calculate and display flow field in the vicinity of the submersed body is an exceptional advantage of FVM in scope of CFD. It is very useful to know spatial distribution of velocity, pressure, temperature and many other values used in aerodynamics. Moreover, the value of the drag coefficient, the scalar fields of absolute pressures, temperatures and Mach numbers were also calculated. In this paper, we present three typical examples of these scalar fields, from the subsonic to transonic and to supersonic regions:

Subsonic flow
In the first case, the second computational example, where M = 0.6, the subsonic flow field around missile is described, see Figs 9 and 10. In the case of scalar field of absolute pressure (Fig. 9), a significant increase can be seen at the leading point of the missile. This increase logically occurs in all other cases. In the scalar temperature field, Fig. 10, only weak (in the range of 20 K) temperature difference between the environment and the body-guided missile can be seen. This is because of moderate subsonic velocity of missile flight. Field of Mach number distribution is presented in Fig. 11 in order to better illustrate the relations between velocity, pressure and temperature field.

Transonic flow
The case of point of the highest drag coefficient occurred at velocity M ≥ 1 was also solved in order to illustrate pressure and temperature scalar fields in the transonic region. In case of scalar field of absolute pressure (Fig. 12), the emerging shock waves can be recognized. Temperature field (Fig. 13) indicates an increase in surface temperature of the missile body up to 50 K in contrary to previous case. Also the field of Mach number distribution is added (Fig. 14).

Supersonic flow
The last illustrative case shows a higher speed M = 2.8. In the scalar field of absolute pressure, see Fig. 15, a decreasing angle of the Mach cone around the leading point, control surfaces, and missile wings can be seen. The temperature difference between the missile control surface and the environment caused by aerodynamic heating reaches 250 K (Fig. 16). Mach number distribution is presented in Fig. 17.  Obtaining such complex image of missile environment in the vicinity of the body is technically much more difficult to reach by means of experimental aerodynamics. Therefore, these CFD methods give us an estimative insight to possible physical processes.

Conclusion
As mentioned in the introduction, FVM is "modern" application of CFD which represents quite a robust tool for obtaining aerodynamic properties of solid object with complex shape. Aerodynamic properties of bodies are fundamental input to many other technical simulations in scope of flight mechanics or ballistics. Of course, we always have to take a special care of experimental validation. Unfortunately, there was not a practical opportunity to realize any serious wind tunnel test in order to check the calculated data. In addition, there are no loosely available aerodynamic characteristics dealing with air-to-air missile bodies, which we can compare because of worldwide restrictions in defense industry.
This paper describes the calculation of one the most important aerodynamic characteristic of AIM-9M Sidewinder missile using FVM -the coefficient of drag related to Mach number. Of course, we can investigate all components of resultant aerodynamic force and moment, but the relation cD(M) is a representative and sufficient example of our effort. A matrix of calculations was performed in the range of input Mach number from 0.4 up to 3.2 with zero angle of attack. Aerodynamic drag force was obtained from these simulations and subsequently, it was recalculated to the drag coefficient in sense of virtual imitation of real measurement. The obtained distribution cD(M) gives us a basic idea of the course of possible flight performance of the missile AIM-9. Thereafter, one can continue with calculations of cD(M) at cases with different angle of attack at given range of M and so enhance herein presented results. In addition, we can perform the same task for another remaining aerodynamic characteristics determination, e.g. lift force or pitching moment dependency on angle of attack and Mach number. As mentioned above, we were not able to compare these characteristics with any available experimental data. Therefore, we used another possible kind of verification based on comparison with a calculation based on a different principle.
All data obtained here by CFD are consistent with the physical nature of the fluid flow around a slender body shaped as a missile. The calculated aerodynamic characteristics comply with common expectations in an acceptable range of results. Thus, we can consider the setting of the simulation as correct, but no fine. In order to achieve a more accurate tuning of this class of tasks, it will by very useful to repeat the above described process on an example of previously measured missile whose aerodynamic database is not restricted.