Prediction of aerodynamic coefficients for irregularly shaped body using numerical simulations

Article history: Received 20 February 2018 Received in revised form 26 April 2018 Accepted 2 May 2018 An evaluation of aerodynamic forces and moments that act on an irregularly shaped body (fragments of HE projectile) moving at high velocities was made, using numerical simulations, analytical models, and CAD methods. Using the results obtained for aerodynamic forces and moments, and known values of body exposed area, aerodynamic drag and lift coefficients were determined for different body orientations and different flow velocity. Analysis of the influence of the body front surface and body slenderness on the position of the maximum CD value (on the CD (Ma) curve) was performed.


Introduction
*Aerodynamics studies the laws of air movement and its interaction with solid bodies in it. Aerodynamics is closely related to fluid mechanics, thermodynamics, and gas dynamics.
External aerodynamics deals with the prediction of force, moment and heat transfer to bodies that move through the air. It studies the generation of force and moment on various aero profiles, wings, aircraft configurations and other bodies moving through the air (Anderson, 2016).
Irregularly shaped high-speed bodies (such as HE projectile fragments) are characterized by high velocities, viscous and compressible flow, dominant resistance due to pressure, shock waves in front, laterally and behind the body, predominantly turbulent flow due to high velocities, unsteady flow fields of velocities and pressures, and most often by separating the boundary layer from the surface of the body during its movement through the air (Buresti, 2000).
The flow in which the density of the fluid is constant is incompressible. All flows are more or less compressible, since completely incompressible flows do not exist in nature. Flows in which the Mach number is less than 0.3 are considered incompressible (Anderson, 2016) and there are several fluid flow modes using this criterion [1]: subsonic flow (M <1 throughout the flow), transonic flow (mixed regions where M <1 and M> 1), the supersonic flow (M> 1 in the entire flow) and the hypersonic flow (high supersonic velocities, M> 5).

Aerodynamic forces, moments and coefficients
Aerodynamic forces and moments arise as a result of pressure distribution and tangential stress on the surface of the body moving through the atmosphere. The pressure (p) acts normally on the surface of the body, while the tangential stress () acts perpendicular to the surface of the body. The tangential stress is caused by friction between the body and the air that flows around it. The main goal of aerodynamics is to determine p and  for given body shape and free flow conditions, and with the help of obtained values -determine aerodynamic forces and moments (Anderson, 2016).
The total effect of the pressure and tangential stress, integrated across the entire surface of the body, is the total aerodynamic force FR and total aerodynamic moment MR. The resulting force FR acts in the so-called center of pressure cp. In general case, the body center of mass cm is not located in the center of pressure cp. The resulting force MR acts in the body center of mass and causes its instability during movement.
The total aerodynamic force FR that acts on the body equals the sum of all elemental forces on the total surface of the body. Using the dimensional analysis, the aerodynamic force FR can be presented in the general form as a function of the following variables (Anderson, 2016): where ρ is the free flow density, v is the free flow velocity, l is the reference length, μ is the viscosity coefficient, and a is the free-flow sound velocity. After applying dimensional analysis to expression 1, it can be written as (Anderson, 2016): ( 0,5 2 , , ) = 0 ( 2) where S is the body reference area. The reference area of the body S and the reference length of the body l are chosen arbitrarily. For different body forms, S and I can be different things. For an axisymmetric body (e.g., spheres) S is the crosssectional area, and l-diameter (Anderson, 2016).
In expression 2 there are non-dimensional coefficients that are important in aerodynamics (Anderson, 2016): Mach number: = If the first part of the expression 2 is written in the form = /0,5 2 , where is the nondimensional coefficient of the total aerodynamic force FR, then it follows that (Anderson, 2016): This means that the coefficient of total aerodynamic force is only a function of Reynold and Mach numbers. Thus, expression 2 is reduced from five to two dependent variables, Re and M.
In addition to aerodynamic forces and moments, in aerodynamics dimensionless aerodynamic coefficients are very important. Aerodynamic coefficients are used to determine the aerodynamic characteristics of a body (i.e., aircraft, projectile, wing). They bring all bodies into the same arena by having a ratio of forces rather than simply using the forces. So a large body might have more lift than a small body but have a smaller lift coefficient. Using the aerodynamic coefficients, efficiencies can quickly be compared.
If ρ and v define the density and free-flow velocity, the expression for the dynamic pressure q in the free flow becomes (Anderson, 2016): Since the parameters S and l are defined as the reference area and the reference body length, the dimensional coefficients of the force (Fi) and the moment (Mi) can finally be expressed in the following way (Anderson, 2016): Mach number: = Although aerodynamic coefficients are significant for estimating the body trajectory, there is very little available data on their values for irregularly shaped bodies. A review of the literature in this field showed that very few researchers were engaged in more serious research in order to estimate these coefficients taking into account the stochastic shape of the body, its velocity and its real exposed area.
In general, there are very few publicly available data on the values of the drag coefficient CD. Investigation of the drag coefficient during the flight of irregularly shaped bodies through the atmosphere is a very complex, and experimental tests used are expensive.
In most studies (Dunn and Porter, 1955;Ramsey et al., 1978;AASTP, 2006;Moga and Kisielewski, 1979;McDonald, 1980;McCleskey, 1988;Miller, 1990;Haverdings, 1994) where experimental tests were performed to determine the value of the drag coefficient CD values for projectile fragments, the aim was to determine the value of CD in the subsonic regime because the value CD in that zone is critical to determining the total range of fragments to determine the real zone of danger after a possible explosion of military warehouses or accidents in open-air ammunition stockpiles.
A large number of researchers in their papers indicate that the drag coefficient for fragments can be taken as a constant (Crull, 1998;Zehrt Jr and Crull, 1998;Crull and Swisdak, 2003), arguing that the fragments in the initial phase of their trajectory move at speeds up to several times higher than the local sound velocity (CD(Ma) curves usually changes very little in supersonic motion regime) However, the drag coefficient CD variation may have a significant effect on body movement, and this fact is rarely emphasized in the literature.
Analysis of the literature dealing with drag coefficient did not find any research involving numerical simulations of airflow around real fragments (produced by the detonation of HE projectiles), where the phenomenology of the movement of the actual fragment in the atmosphere and its aerodynamic parameters would be studied more seriously.
No experimental data on aerodynamic lift coefficients CL were found for irregularly shaped fragments in the available literature, but it is possible to find in the literature the data on the CL coefficient for some other forms (for example, data for a flat plate, parallelepipedic form) and draw certain conclusions.

Numerical simulation method
Methods of numerical simulations using CFD (Computational Fluid Dynamics) are an important aspect of modern research because they are complemented by experiments and analytical models, reducing total time and labor costs.
In order to determine the values of the components of the aerodynamic force, torque and exposed area of the body, on the basis of which the aerodynamic coefficients are determined, the numerical simulations (for different velocities of the three-dimensional, compressible, turbulent, steady flow around the irregularly shaped body) were performed in the Ansys Fluent CFD package for fragment in different orientations.

Basic equations
The basic equations of flow represent continuity equation, momentum equation and energy equation.
Physical principle behind continuity equation is that mass can be neither created nor destroyed. Continuity equation can be written in the form of a partial differential equation (Anderson, 2016): This equation relates the flow field variables at a point in the flow. Equation (9) hold in general for the three-dimensional, unsteady flow of any type of fluid, inviscid or viscous, compressible or incompressible.
Physical principle behind momentum equation is that force is equal to time rate of change of momentum and can be written in differential form as (Anderson, 2016): The momentum equations for a viscous flow [Equations 11 to 13] are called the Navier-Stokes equations. For an incompressible flow, where ρ is constant, the primary flow-field variables are p and V. The continuity and momentum equations obtained earlier are two equations in terms of the two unknowns p and V. Hence, for a study of incompressible flow, the continuity and momentum equations are sufficient tools to do the job. However, for a compressible flow, ρ is an additional variable, and therefore we need an additional fundamental equation to complete the system. This fundamental relation is the energy equation (Anderson, 2016).
Physical principle behind energy equation is that energy can be neither created nor destroyed; it can only change in form. This physical principle is embodied in the first law of thermodynamics. It can be written as (Anderson, 2016): where ̇′ and ̇′ represent the proper forms of the viscous terms.
With the energy equation, we have introduced another unknown flow-field variable e. We now have three equations, continuity, momentum, and energy, which involve four dependent variables, ρ, p, V, and e. A fourth equation can be obtained from a thermodynamic state relation for e. If the gas is calorically perfect, then (Anderson, 2016): where c ϑ is the specific heat at constant volume. Equation (15) introduces temperature as yet another dependent variable. However, the system can be completed by using the perfect gas equation of state (Anderson, 2016): where R is the specific gas constant. Therefore, the continuity, momentum, and energy equations, along with Equations (15) and (16) are five independent equations for the five unknowns, ρ, p, V, e, and T .

Numerical simulations procedure
The method of numerical simulations of air flow around an irregularly shaped body (fragment) consisted of: a) Digitalization of the fragment model b) Physical domain discretization c) Characterization of materials d) Initial and boundary conditions e) Solver and turbulence model selection f) Aerodynamic force, moment and exposed area od fragment estimation (UDF Script)

Digitalization of the fragment model
The irregularly shaped body with which the numerical simulations were performed for different orientations is presented in Fig. 1.
Since the 3D scanner was not available, a body was digitized using CAD software. The three dimensional model of the body (Fig. 2) was made using CAD methods, modeling a real body in three projections, and then manipulating it: extruding in the direction of three coordinate axes, connecting extruded projections and determining their cross section as the final 3D model. This 3D model was later imported in Ansys Fluent.

Physical domain discretization
Since these type of bodies have a stochastic and irregular shape, it would be very difficult to define their geometry using only point by point technique in one of the older preprocessors such as, for example, Gambit.
For this reason, Ansys Fluid Flow Design Modeler was used, which enables the introduction of already finished CAD models into Ansys Fluent software.
The mesh of the numeric model was unstructured (Fig. 3), consisting of around 626 000 polyhedral elements. The mesh was especially fine around the body where high gradients are expected (Fig. 5) to reduce the numerical error. The Ansys Fluent coordinate system was placed in the body center of mass, with the initial orientation of the body adopted to coincide with the principal axes of inertia.

Characterization of the resistive medium
Air is modeled as homogeneous, isotropic, ideal gas with pressure-temperature dependent density , specific heat Cp, thermal conductivity k and dynamic viscosity .
During the high-velocity motion of the fragment, the pressure nearby the fragment boundary can reach values multiple times higher than the atmospheric pressure which results in significant local change in density. Fluent used following form of ideal gas law for compressible flows (Fluent, 2010): where: patm -atmospheric pressure, p -relative pressure, R -universal gas constant, M -molar mass, T -temperature determined from the equation of energy.
For material (ideal gas) following parameters were used: specific heat Cp = 1006.43 J/kgK, thermal conductivity k = 0.0242 W/mK, and the molar mass M = 28.966 kg/kmol.
The influence of temperature T on the dynamic viscosity of air can be significant for large variations in temperature. It was determined by the Sutherland model (Fluent, 2010): where 0 -dynamic viscosity, T0 -reference temperature, S -Sutherland constant. For air following values were used: 0 = 1,71610 -5 kg/ms, T0 = 273,11 K, and S = 110,56 K (Fluent, 2010).

Initial and boundary conditions
In the simulations, fragment was considered stationary and airflow around it was analysed. Numerical simulations for 24 fragment orientations were performed for angles 0-360 with 15 angle increments. Fig. 4 shows the schematic position of the fragment in numerical simulations. The velocity vector was directed in the positive direction of axis X of the coordinate system set in the body center of mass (Fig. 4).
Since the flow velocity around the body is always supersonic at the beginning of its flight, significant effects of compressibility and viscosity, as well as shock waves occur.
For the compressible and isentropic (reversible adiabatic thermodynamic process) flow of the ideal gas, expressions from the compressible fluid mechanics (Fluent, 2010) are used: where: p0 -isentropic (stagnation) pressure, pstatic pressure, T0 -isentropic (stagnation) temperature, T -static temperature,  -adiabatic exponent (ratio of specific heats). The value of  for the air is 1.4.
At the end of the domain, Pressure Farfield condition was used, which is commonly used in Fluent in aerodynamic simulations, where the effect of compressibility is dominant.

Fig. 4: Schematic position of the body (fragment) in numerical simulations
The No Slip condition is defined on the surface of the body, which means that the relative flow velocity on the surface of the body is equal to zero.
Boundary condition -the wall is used in case when the viscous effects cannot be ignored and is relevant to most fluid flow situations (Fluent, 2010).

Solver and turbulence model selection
Solver settings include solver type selection, discretization scheme, solution initialization, and convergence monitoring. There are two basic types of solver in Fluent (Fluent, 2010): pressure based and density-based solver. According to the recommendation (Fluent, 2010) for use with compressible flows, a density-based solver was selected in the simulations, where mass, flow and energy equations are determined as the Navier-Stokes equation system in integral form for an arbitrary control volume.
Using the Navier-Stokes equation [equations 11-13] in the density-based solver in certain cases (when there is a large difference between the velocity of the flow and the local sound velocity) results in lower convergence, and in this case the preconditioning technique (Fluent, 2010) is used.
By determining the Reynolds number for the air flow around the body (kinematic viscosity of air is 1,5110 -5 m 2 /s for density of 1,2 kg/m 3 ) for the velocity of the flow above 0,6 Mach, the Re number exceeds 10 5 , which means that the flows around the fragment with these velocities are turbulent, although the time of movement of this body through the atmosphere is relatively short.
According to the recommendations (Fluent, 2010), the Spalart-Allmaras turbulence model was used in the simulations. This is a relatively new physical model of turbulence. It has been developed specifically for aerodynamic applications (especially in the aerospace industry) and has proven to be effective for the boundary layers with high pressure gradients, and has been particularly effective for transonic flows around the aero profiles, including flows with significant separation of the boundary layer (Pope, 2000).

Aerodynamic moment estimation (UDF Script)
A program (in C programming language) was written that determines the aerodynamic forces, moments and exposed areas for fragment (perpendicular to airflow velocity vector) for all three coordinate axes. For each cell on the body, the forces are determined in three directions, as follows (Kljuno and Catovic, 2017a): where: Fxi, Fyi i Fzi-force components for each body cell, ps -pressure on the surface of the wall, and Sxi, Syi and Szi -the elementary exposed body surfaces vertical to the respective coordinate axes.
Then the components of the total pressure force Fx, Fy and Fz are determined as (Kljuno and Catovic, 2017b): The aerodynamic moment is determined for each cell of the discrete model as the vector product of the cells vector radiuses and the aerodynamic force acting on that cell (Kljuno and Catovic, 2017c): The components of moment M for each cell of the body using the function NV_CROSS (ri, Fi) are determined by the program, and then total moment for all three axes is obtained by summing all elementary moments.
The moment of the forces tangential to the surface ("friction forces") of each boundary cell is neglected with respect to the moment due to the dominant pressure forces.
Exposed areas S of fragment, perpendicular to axes x, y and z are determined as (Kljuno and Catovic, 2017a): where Sxi, Syi i Szi are exposed areas of each individual cell of fragment, perpendicular to x, y and z axes, respectively. Exposed area of fragment is determined by summing all individual area of each cell using the function F_AREA in C program. Solver of Ansys Fluent is set by dynamically loading the developed UDF (user defined function) program and executing the commands entered into it. Aerodynamic forces and moments are determined using UDF script, specifically using the DEFINE_EXECUTE_AT_END commands (general type macro, executed in the simulation) for each cell of the model.
The UDF program is written so that the results are printed as a table in a separate .txt document.
A numerical model was validated using available experimental data for drag coefficient CD of the cube, determined using expression (Anderson, 2016): where: FD is the drag force, S is the reference area and = 0,5 2 is dynamic pressure. The experimental data (Schamberger, 1971;Hoerner, 1965) of CD for cube flat-on flow (flow perpendicular to the cube side) were compared with values of CD obtained using numerical simulations.
In the process of validation of results, the discretization of space and time, solver and initial and boundary conditions in the case of simulation of the airflow around 3D cube models were the same as in the numerical model of flow around the 3D model of the irregularly shaped body. Fig. 5 gives a comparison of numerical simulation results with experimental data (Schamberger, 1971;Hoerner, 1965) for flat-on orientation of cubes. The difference between values of the CD from simulations and experiments were less than 10%.

Aerodynamic drag coefficient
After the completion of numerical simulations of airflow around 3D model of the fragment, the data on the components of the aerodynamic forces (the force components in the x, y and z axes), the components of the aerodynamic moments (components of the moment in the x, y and z axes), and the exposed areas of the fragment were determined, for each Ma number considered, using the UDF program [equations 19-22].  (Schamberger, 1971, Hoerner, 1965 for cube Based on the obtained results on forces, moments and exposed areas, the values of aerodynamic drag coefficients for the fragment are determined for different values of Mach numbers using expressions: The drag force FD from the expression (24) is in the direction of axis x; this axis is stationary in the simulations and is always in the direction of the velocity vector, Fig. 4). The lift force FL (expression 25) is in the direction of axis y in simulations (Fig. 4), always perpendicular to the velocity vector. The exposed surface of the fragment S is the projection of the total surface of the fragment on a plane perpendicular to the velocity vector, i.e. the plane YZ in Fig. 4. Dynamic pressure of free flow q is determined by the expression mentioned in equation 23.
In Table 1, values of drag coefficient CD for a 3D model of an irregularly shaped body (fragment) were presented, obtained using numerical simulations and expression (24), for different body orientations (0-360) and different velocities (0.6 M to 4 M). Angular increments of fragment orientation of 15 were used in plane XY (Fig. 4). Fig. 6 shows the values of drag coefficient CD, obtained using numerical simulations (Table 1) for different orientations of fragment (0 °-360 °).
Diagram in Fig. 6 shows that the maximum drag coefficient CD is obtained in the case when fragment is in the initial position (orientation of 0° or 360°, Fig. 4), and minimum when the fragment is rotated 90° relative to the initial position, i.e., when the fragment is in a position in which it represents relatively slender body, Fig. 4).
At velocites greater than 1.5 Ma, the fragment has minimum value of the coefficient CD in orientation when it is rotated 105 relative to the initial configuration (the lower red curve in Fig. 6). Other CD values (CD (Ma) curves for other fragment orientations) are in specified interval between CDmax and CDmin.
In Fig. 7 the results of the maximum and minimum value of drag coefficient CD are given, obtained using the simulations as the function of Mach number (CDmax -violet set of dotted points and CDmin -a green set of dotted points) together with the publicly available experimental data for the drag coefficient of fragments (Dunn and Porter, 1955;Ramsey et al., 1978;AASTP, 2006;Moga and Kisielewski, 1979;McDonald, 1980;McCleskey, 1988;Miller, 1990;Haverdings, 1994).
As far as experimental data is concerned, Fig. 7 indicates certain deficiencies that have been observed. One of the disadvantages is that most of these tests contain very little data on the CD values in the supersonic regime. Another disadvantage is that no research has given details of how the exposed surface of fragments is determined which significantly affects the final CD values, as well as the data on dimensions and the mass of fragments. These experimental tests were conducted in the period from 1955 to 1995 and no data are available on whether more recent tests have been conducted with modern measuring equipment. All these parameters play an essential role in determining the value of the coefficient CD, in particular the data on the assumed reference surface of the fragment. Some researchers (McDonald, 1980) state that the drag coefficient of in the subsonic and supersonic flow is assumed constant (i.e., red and blue curves representing the extreme values of CD in the experimental data sets in Fig. 7), which in reality is not the case, and it is probably a certain approximation of the CD, due to the lack of data.  Fig. 7: Comparation of max. and min. drag coefficients obtained with simulations and experimental drag coefficients (Dunn and Porter,1955;Ramsey et al., 1978;AASTP, 2006;Moga and Kisielewski, 1979;McDonald, 1980;McCleskey, 1988;Miller, 1990;Haverdings, 1994) Also, it has already been mentioned that no research has been conducted in which numerical simulations of airflow around the real 3D model of fragments have been found (mostly numerical simulations of flow around the symmetrical bodies are found in the literature) so the results of the simulations performed in this paper could not be compared to results of numerical simulations of other authors.
If the simulation results are compared with the experimental data, the first thing that is noticed in Fig. 7 is that the minimum values of a CD obtained by numerical simulations are considerably smaller than the minimum values of a CD obtained experimentally. This difference may be the result of a different way of defining the exposed surface of the fragment (some authors state that the entire outer surface of the fragment can be taken as its reference surface), so it is important for each researcher to emphasize which reference surface of fragment was used in his research, so the results could be compared.
The second concluson that can be deduced from the diagram in Fig. 7

Ma
Vertical tunnel (Moga, 1979) Air tunnel (Miller, 1990) Gas gun (Miller, 1990) Cdaver (Dunn, 1955) Cdaver (Ramsey, 1978) AASTP-1 (Heiser, 1979) Cdmin (McDonald, 1980) Cdmax (McDonald, 1980) Cdaver (McDonald, 1980) MISDAC (1994 in experimental data the maximum CD value is in the range of 1.2 (Haverdings, 1994) to 1.6 (McDonald, 1980. The curve of the maximum value of CD, obtained with numerical simulations (dotted line of purple color, Fig. 7), is closest to the experimental data (data in central part of the diagram in Fig. 7), indicating that in the literature the value of CD for the fragments is probably overestimated. In reality, for the same shape of fragment, when moving through the air, the value of the coefficient CD will vary between CDmax and CDmin (conclusion supported by recent papers from NASA (Murman et al., 2005;Murman, 2010), depending on the fragment orientation and its velocity, but these CD values are likely to be considerably less than the values of CD (for fragments) given in the literature.
The interesting thing is noted regarding CD(Ma) curves. Namely, the curves in Fig. 6 show that at certain orientations of the fragment (when it represents the relatively slender body, Fig. 4), the maximum value of the CD is shifted to larger Mach numbers. The maximum value of the drag coefficient CD (Table 1 and Fig. 6) is predominantly in the range of 1.2 M (for orientations: 0°, 15°, 30°, 45°, 105°, 135°, 150°, 180°, 195°, 210°, 225°, 315° and 345°) to 1.5 M (for 60°, 75°, 90°, 245° and 255°). At the fragment orientations of 120°, 165° and 300°, the maximum CD value was at 1.3 M. Interestingly, for the orientation of the fragment of 270°, the maximum CD value was at 3 Ma, and for 280°, this maximum was at 2 Ma.
The reason for this shift of the maximum CD value for a fragment (irregularly shaped body) is difficult to estimate because in different orientations the fragment has a different shape of the front surface and different values of the exposed surface. By examining the available literature, this phenomenology (in the case of fragments of irregular shape) was not considered.
Generally speaking, in order to obtain functional dependence between the position of the peak value od CD and some of the parameters (ie body shape, the shape of the reference area, etc.), only one variable should be considered and the others kept constant.
In this regard, an analysis of the influence of the body's frontal surface (through the variation of prism apex half-angle -idealized geometric shape) on the position of peak values of curve CD(M) was performed, where the reference area of the body was the same (constant) for all cases.
For this analysis, numerical simulations of airflow around the body have been made with different values of prism apex half-angle, while having the same reference area. Numerical simulations of aerodynamic flow were performed in Ansys Fluent. Fig. 8 is a schematic representation of a body with different frontal shape, with which numerical simulations have been performed. In the analysis prisms with the following apex half-angles were taken into consideration: 30°, 45°, 60° and 90°, as in Fig. 8). From Fig. 8 it can be seen that the shape and value of the reference area (rectangle shape) did not change in these bodies (dimensions of height b (50mm) and body width c (50mm) were constant for all shapes). Optimized (reduced cell size around the body) unstructured mesh (1000 000 cells) and boundary conditions wall and pressure far-field ( Fig. 9) were used in the numerical model. The flow velocity vector was directed in the positive direction of axis X of the coordinate system set in the body center of mass (Fig. 10). Enclosure mesh around body was a cube ( Fig. 9) with 1m long side. Air is modeled as homogeneous, isotropic, ideal gas with pressure-temperature dependent density, specific heat, thermal conductivity and dynamic viscosity. Simulations were carried out for a flow velocities of 0.9, 1, 1.1, 1.2, 1.3 and 1.5 Mach since peak values of CD are usually found in this interval.
A density-based solver was used in the simulations, where mass, flow and energy equations were determined as the Navier-Stokes equation system. Spalart-Allmaras turbulence model was used in the simulations since it has proven to be effective for the boundary layers with high pressure gradients, and is effective for transonic flows (Pope, 2000). Residual tolerance in numerial simulations was set to 10 -5 (Fluent, 2010). Table 2 gives numerical simulation results for examined prism shapes (prisms with four different apex half-angles) showing Mach number at which the peak value of CD is obtained for individual cases (for each apex half-angles of prisms).
The results from Table 2 show that the position of an extreme values of drag coefficient CD is at different Mach numbers for different body's frontal surface (at the same time with the equal value of the reference area -area perpendicular to flow velocity vector), which means that the shape of the front surface of the body affects the position of the maximum value of CD (M).  Fig. 11 shows the drag coefficient CD(M) curve for the prism with different apex half-angles, for tested velocities.
The diagram in Fig. 11, as expected, shows that the minimum values of the drag coefficient CD were obtained in the case of prism apex half-angle of 30, while for prism apex half-angle of 90° (which practically represents a parallelepiped body), the highest CD values for the full range of Mach numbers were obtained. These results can be presented differently, using the slenderness of a body. The body's slenderness can be defined (i.e., as for the axisymmetric projectiles) as the L/D ratio, where L is the dimension of the body in the direction of the velocity vector, and D is the equivalent diameter of a body (referring to axisymmetric bodies). The reference area was constant and of the same shape for tested bodies in the analysis.
Dimension L is easily determined for the prism (Fig. 8) from the trigonometry, but since the crosssection in this case is not a circle but rectangular, it is necessary to define a model to determine the value D for rectangular reference area. This can be done by first determining the polar moment of the inertia of the reference area S of the body with rectangular reference area about an axis of the velocity vector (x in this case). This axis passes through the center of mass of a body and is perpendicular to area S in Fig. 12.
Once the polar moment of inertia values are known for this rectangular area, equating these values with the values of polar moment of inertia for the circular reference area, the value of the equivalent diameter Dekv can be found. Fig. 12 shows a schematic representation of the dimensions (c and b) of rectangular reference area S. The polar moment of inertia of the rectangular area (i.e., area S, Fig. 12) is obtained using the expressions: where c is body dimension perpendicular to the velocity vector direction (dimension in the direction of axis y, Fig. 12), b -is body dimension perpendicular to the velocity vector direction (dimension in the direction of axis z, Fig. 12). The dimensions c and b did not change in this analysis, only the dimension L (length of a body) was varied by using different apex half-angle of the prism. The polar moment of inertia of the circular area is obtained using the expression: where D is the diameter of the cross-section (reference area) of the axisymmetric body. It is this parameter D that will be used as the equivalent diameter of the body with rectangular reference area -Dekv, by equating the expressions (10) and (11) where L is a dimension of a body in direction of velocity vector. The data from Table 2 can now be presented as a function of body (prism with different apex halfangles) slenderness. Table 3 gives the results of numerical simulations (Mach numbers for extreme values of CD) for prisms of different slenderness.  Table 3 shows that for the body slenderness 0.3405 (apex half-angle 30) the extreme value of CD is at 1.1 Mach. For body slenderness 0.5256 (apex half-angle 45) the peak value of CD is at a speed of 1.2 M. On the other hand, for a body slenderness 0.8463 (apex half-angle 60), the extreme value of CD is at a speed of 1.3 Mach and for a slenderness 0.0876 (apex half-angle 90°) the extreme value of CD is at a speed of 1 Mach.
Diagram in Fig. 13 shows the results obtained by numerical simulations, as shown in Table 3. The diagram in Fig. 13 shows that in the case of a rectangular reference area such as in a prism, by increasing the body slenderness (at the same time altering the frontal shape of the body), the values of the Mach numbers at which the extreme values CD are obtained, increased to the larger Mach numbers. The data obtained for the prisms of different slenderness were approximated by the linear function y = 0.402x + 0.9691 (Fig. 13), where the correlation coefficient was satisfactory (R 2 = 0,989). These results point to the fact that frontal shape of body and its slenderness can impact the position od the extreme value of CD on curve CD(Ma).

Aerodynamic lift coefficient
The lift coefficient CL is a dimensionless coefficient that relates the lift generated by a lifting body to the fluid density around the body, the fluid velocity and an associated reference area. The choice of the reference area should be specified since it is arbitrary.
The lift coefficient CL was also determined for fragment (irregularly shaped body) in different orientations (Fig. 4), depending on the angle of orientation(0-360).
Based on the results of the numerical simulations and the expression 25, in Table 4 are the values of drag coefficient CL for the 3D model of real fragment for different orientation of the fragment.  Fig. 14 shows the curves of drag coefficient CL for fragment as a function of the Mach number for the different orientation of the fragment (0-360), obtained from the results of numerical simulations.
The curves in Fig. 14 show that the same fragment may, in different orientations, have significantly different values of the lift coefficient. In order to interpret the results, CL data that were available in literature (for certain body shapes) are compared with the results obtained in numerical simulations. Baker et al. (2012) in the Manual of Explosion Hazards and Evaluation (Baker et al., 2012) lists data for the thin plate lift coefficient at different attack angles (different orientation of plate). The lift coefficient CL as a function of plate orientation (attack angle), as indicated by Baker et al. (2012) is given in Fig. 15. The author does not state which flow velocity refers to this curve, but it is assumed that the data relate to supersonic velocities as the author describes the bodies formed by explosion. Fig. 15 shows that for a zero attack angle, the lift force in symmetrical bodies is practically equal to zero. The lift force is approximately linearly proportional to the attack angle, and at higher angles of attack, the lift coefficient reaches its maximum at a certain angle (the so-called critical angle of attack), after which there is a significant separation of the flow from the upper surface of the body so the lift force decreases (and the drag force increases). For aero profiles (in aviation), for example, in this case, it says that there was a stall. After the critical attack angle, the lift coefficient CL decreases (Fig. 15) and reaches the zero value for the 90° attack angle (in this position, the body is also symmetrical, relative to the flow velocity vector). Twisdale and Vickery (1992) investigated models for predicting trajectories of primary and secondary fragments (fragments of walls, concrete, steel elements, pieces of soil, as well as fragments of missiles) during the explosion. In their paper they present the basics of several previously developed models (2D model without drag force; 2D model with drag force; 3D model with drag force; model with drag force, lift force and lateral force) for estimating the drag, lift and lateral force coefficients (for the body with the parallelepiped shape).
They used analytical models (the so-called Cross Flow theory -attempt to estimate the aerodynamic lift and lateral force coefficient as a function of the attack angle only on the basis of the known values of the body drag coefficient in three different directions. In the research, diagrams were presented (Fig. 16) on which they compared the results for the lift coefficient obtained analyticaly (Cross Flow theory), with experimental data for a paralelepiped body.
The obtained numerical simulations results for lift coefficient CL (Fig. 14) can be presented differently in order to qualitatively compare them with the diagrams in Figs. 15-16.   (Twisdale and Vickery, 1992) So in Fig. 17, the curves of the lift coefficient CL are given as a function of the fragment atack angle (orientation of a body in our case) -namely, two sets of data: for body orientations 0-90 and 90-180 (Fig. 4), for different Ma numbers (0.6-4 M).
By comparing the diagrams in Fig. 17 with the diagrams in Figs. 15-16, we can see that the curves of lift coefficients, obtained by numerical simulations (Fig. 17) have a similar trend as the lift coefficients curves in Figs. 15-16 (Baker et al., 2012;Twisdale and Vickery, 1992).
Based on the results of numerical simulations, for the orientations 0, 90, and 180, the obtained coefficient CL value is close to zero (not zero because the fragment is not a symmetrical body), similarly as in Figs. 15 and 16. Also, the maximum value of the coefficient CL is similar because in numerical simulations, the maximum CL is obtained for orientation between 40 and 55, while for example in Baker's case (Baker et al., 2012) this maximum CL is at atack angle of 40 for a flat flat plate, and in Twisdale's case (Twisdale and Vickery, 1992) it is between 40 and 60. Anderson (Anderson, 2016) states that this angle may be 55° in the case of hypersonic flow.
What is interesting in Fig. 17 is that the obtained diagrams have a certain symmetry (of course, not complete symmetry because the body is not symmetrical), i.e. the curve CL (Ma) from the upper set of results (orientation of the fragment: 0° -90°) are quite similar the curve CL (Ma) from the lower set of results (the orientation of the fragment: 90°-180°), only in the second case we have negative values of CL.
If we go back to Fig. 4 we can see the reason for this. Namely, if, for example, we compare the orientation of the fragment of 60 ° and 120 °, it can be seen that these orientations of the fragment are symmetrical in relation to the velocity vector direction (Fig. 4). This also applies to the orientations of the 15° and 165°, 30° and 150°, 45° and 135°, 75° and 105°, respectively. Thus, for these orientations, the CL values for fragment are similar but with different sign, so and CL(Ma) curves represent a kind of symmetry relative to the direction of airflow vector (Fig. 17).
Presented values of coefficients CD and CL for irregularly shaped bodies (such as fragments of HE projectiles), obtained in this paper, can be used when predicting trajectories of such bodies.

Conclusion
Using the analytical models, CAD methods and numerical simulation methods presented in the paper, the values of aerodynamic forces and moments that act on a high-velocity (HE projectile) fragment were determined. Based on the results, and  0  degrees  15  degrees  30  degrees  45  degrees  60  degrees  75  degrees  90  degrees  105  degrees  120  degrees  135  degrees  150  degrees  165  degrees  180  degrees  195  degrees  210  degrees  225  degrees  240  degrees by knowing the value of the exposed fragment area for each fragment orientation, an aerodynamic drag coefficient CD and aerodynamic lift coefficient CL for the given body were determined for different orientations of body and at different velocities. It was found that the same body (in this case a fragment of HE projectile) can have a whole range of CD and CL values. This fact is often neglected in the literature, and instead, some average value of these coefficients is taken into account usually. In reality, for the same shape of the fragment, during the motion through the atmosphere, the value of the coefficient CD will vary between CDmax and CDmin, depending on the fragment orientation and its velocity. Also, lift coefficient CD will also vary, which can lead to an unpredictable vertical motion of such a body.
In this paper, the minimum values of CD, obtained by numerical simulations, are found to be less than the minimal value of the CD obtained experimentally. The difference may be the result of a different way of defining the exposed area of the fragment.
Evaluation of the effect of the front surface shape of a body and slenderness of a body on the position of a maximum CD value (on CD (Ma) curve) was performed. It has been found (for an idealized geometric shape -prism) that the shape of the front surface of the body can affect the position of the maximum CD value on the curve CD (Ma). Similar results were obtained for fragment with an irregular shape.