Contribution to Compressibility Modelling in the Estimation of Forces Acting on Projectile Fragments

The compressibility model was developed and applied to generalized model for prediction of aerodynamic forces acting on irregularly shaped body, such as HE projectile fragments. Model assumes adiabatic compression of air in front of high‐velocity fragment since the motion of the fragment is an extremely fast process relative to the heat transfer process. The equation of the state of the ideal gas is adopted. Analysis of results and comparison with results obtained by numerical simulation (CFD software) and experimental data show this correction model significantly reduces the relative error of aerodynamic force modelling in the relevant area of the high velocity of the fragment, when the air compression is significant.


Introduction
Complete aerodynamic characteristics generally do not exist for the potentially wide variety of body shapes. For the calculation of the body trajectory, apart from the data on projected (exposed) surface area of the body, velocity data and data on the density of the fluid through which the body is moving, one needs also the data on the values of aerodynamic force acting on the body at any given moment. Values of force are crucial in the calculation of fragment trajectory and its movement through the air.
It is not practical to determine the values of aerodynamic forces and moments using numerical simulations (CFD) for each potential shape of a body, so it is useful to define a generalized model to estimate the values of total aerodynamic force and aerodynamic moment acting on a body with an arbitrary shape.
Examples of irregularly shaped bodies are: • primary fragments made after the detonation of high explosive projectiles (natural or controlled fragmentation), missiles, bombs, • primary and secondary fragments originating from improvised explosive devices (IED), • splinters formed after the rupture of various structures (different materials) due to the effects of severe storms (tornadoes and hurricanes, for example, both contain strong rotating winds that can cause significant damage), • primary and secondary fragments resulting from the explosion of ammunition storage sites, • meteoroid bodies coming from outer space.
Generally, most objects in nature are irregularly shaped, so the usefulness of this model (in regard of estimation of a trajectory, velocity, energy, and stability of these bodies) could be significant for many areas of research. Usually, determination of force required for calculation of body centre of mass trajectory (there is no mentioning of moment coefficients in the literature), is done using aerodynamic force coefficients.
These coefficients (coefficients of total aerodynamic force components) are usually determined analytically or experimentally. Twisdale [1] used cross-flow theory in order to analytically estimate the aerodynamic coefficients for a uniformly random spatial orientation of the body, knowing the aerodynamic coefficient values for precisely defined directions of the body with continuous geometry. This approach has been used earlier to develop the wind axis aerodynamic forces as a function of an angle of attack for slender cylinders knowing only the drag force coefficients for the body in normal flow to the major body axes. Twisdale presents cross-flow aerodynamics for rectangular parallelepipeds (approximation of irregularly shaped body), where for different parallelepiped lengths (slenderness ratio), drag, lift and side force coefficients are estimated as a function of attack and roll angle, using analytical expressions, where different corrective factors (skin friction and aspect-ratio correction) are also used. It is not explained in the paper how the exposed area of a body is determined for an arbitrary orientation of a body. Based on the value of the drag, lift and side force coefficients, and using the values of dynamic pressure and the reference area of the body, a total aerodynamic force is estimated.
As far as experimental data is concerned, the first thing to notice in the literature is that there is no experimental data on aerodynamic lift coefficients for bodies with irregular shapes. Furthermore, one of the disadvantages noticed is that most of the experimental tests contain very little data on the drag coefficients values in the supersonic regime. Another disadvantage is that no (available) experimental research has given specific details of how the exposed surface of fragments is determined, as well as the data on dimensions and the mass of fragments. These data significantly affect the drag coefficient values obtained in experimental tests. Most experimental tests for determination of drag coefficients for irregularly shaped bodies were conducted in the period from 1955 to 1995, and no data are available on whether more recent tests have been conducted with modern data acquisition equipment. Some researchers [2] state that the drag coefficient in the subsonic and supersonic flow is assumed constant, which in reality is not the case, and it is probably a certain approximation of the drag coefficient, due to the lack of data. Other researchers suggest that the drag coefficient for fragments can be taken as a constant [3], arguing that the fragments in the initial phase of their trajectory move at speeds up to several times higher than the local sound velocity. However, the drag coefficient depends on the shape of a body, as well as its velocity, and coefficient variation may have a significant effect on body trajectory, a fact that is rarely emphasized in the literature. Generally speaking, the task of estimation of aerodynamic forces and moments acting on a body of an irregular shape is quite complex since these irregularly shaped bodies do not have a continuous surface (their shape is stochastic). For this reason, in our model, a tri-axial ellipsoid was used as a shape that can approximate an arbitrary body shape of different dimensions. This model is also a follow-up on our previous work where the exposed area of the irregularly shaped body was estimated.
In our model [4], irregularly shaped body geometry can be defined parametrically, and it is possible to estimate the values of total force and moment acting on these bodies in a very short amount of time (in a matter of seconds). This is a significant advantage of the model in comparison with the method of determination of aerodynamic force and moment using CFD methods.
In the available literature, no similar model for the estimation of aerodynamic forces and moments for the irregularly shaped body was found, so our model presents a step towards a better understanding of flight dynamics of irregularly shaped bodies.
In the previous paper [5] of the authors, a generalized model for the estimation of aerodynamic forces and moments for irregularly shaped bodies (IE fragments of HE projectiles, secondary debris) was described, where irregularly shaped body was approximated with a triaxial ellipsoid. This model is a follow up on the previous work [6] of the authors where the exposed area of the irregularly shaped body is estimated. These models can be successfully used for the estimation of a trajectory parameters for a body with an irregular shape (i.e. shrapnels), and the obtained data can be used for the prediction of total range of fragments and lethal zones of HE projectiles.
Our general model for obtaining aerodynamic force is based on a change in the fluid momentum in the surrounding of a fragment. The energy that is spent on the movement of fluid (air) is directly taken from the kinetic energy of the fragment. In addition to the energy that is spent on fluid motion, part of the kinetic energy of the fragment is spent on air compression in front of the fragment. If the adiabatic change of the gas state is assumed at this compression, then the work spent on this compression is directly transformed into internal energy by heating the gas (air), according to the first law of thermodynamics. After the passing of fragment, this energy, originally in the form of the kinetic energy of the fragment, is further dissipated into the surrounding air.
The idea of a compressibility model presented in this paper is that by the analogy of a compression spring, as shown in Fig. 1, the mechanical work of air compression per unit of air mass (specific work) is calculated by determining one parameter of state (i.e. pressure) and assuming the adiabatic change in the gas state (since it is an extremely fast process in relation to heat transfer). Pressure can be estimated for each point directly, close to the fragment surface, based on the previously developed model [5] and used to determine the specific volume (according to the equation of adiabatic change), and to determine the size of air compression at that location. In this way, a specific work of adiabatic compression can be determined. Considering this specific work in the unit of time and integrating it in the domain of the disturbed air in front of the fragment, it is possible to determine the equivalent correction of the aerodynamic force which dissipates the same amount of mechanical (kinetic) energy of a fragment.

Physical Model
The model for obtaining aerodynamic force, as reported in [5], is based on the law of momentum change (similar to the Newton model that is now used for hypersonic motion) in order to estimate the value of the force acting on the body. Fig. 1 is a schematic representation of the surface of the ellipsoidal body (approximating the fragment) to which the fluid flow with velocity V is directed, where the normal vector for the surface dA is denoted by n and the angle between the velocity vector and the normal vector is denoted by δ. Fig. 2 shows a schematic way of solving the problem (the output velocity vector is assumed to be perpendicular to the normal vector): fluid elements interact with each other immediately after impact on the fragment surface and their trajectory is such that the fluid flows around the fragment. Fluid flow simulations around the fragment, using CFD software packages [7], show that the fluid output velocity on the exposed side of the fragment surface is tangential to the surface, indicating that the approach shown in Fig. 2 is a good representation in terms of modelling.
According to model [1], the final expression for the total aerodynamic force is: Here ρ is the density of the fragment material, Vul is the input velocity of the flow, Vizl is the output velocity of the flow, γvz is the angle between the velocity vector and the plane of the projected surface of the fragment (whereby these unit vectors are determined by the method explained in reference [6]), and a1, b1 are half-axes of the ellipse defining a curve in the plane separating the exposed part (to the flow) of the fragment from the rest of the fragment [6]. The integration of expression (1) is performed using the x and y coordinates in the plane of the projected fragment surface, where the integration boundary is defined by the ellipse half-axes a1 and b1. Inlet velocity from expression (1) is equal to flow velocity: where Vul is magnitude of the inlet velocity and B is the vector defined as: Here a, b, and c are the largest dimensions of fragment (approximated with triaxial ellipsoid) in three perpendicular directions.

Fig. 2 Elemental control volume analysis
Regarding the introduction of compressibility in this model, as a first step specific work (per unit of mass) can be, generally, defined as: Here p is the pressure, and vs is the specific gas volume. If in Eq.
(3) the following thermodynamic equation is used: s Here the index 1 refers to the condition of the undisturbed flow and 2 to the state close to the surface of the fragment at relatively high velocities.
Let us consider an elemental fluid mass in front of the fragment (in the control volume in the form of a cylinder, whose axis is in the direction of the velocity vector, Fig. 3), which is still in state 1: ∆ is the element of fragment surface projected on a plane perpendicular to the direction of the velocity vector of fragment centre of mass, ∆L is the elementary length of the fragment motion for the time increment ∆t, and ρ1 is the specific mass of undisturbed air ( ) The work of the piston (see Fig. 3) in this case is:

Fig. 3 Analogy of an air compression
Similarly the compression work, carried out on the elementary volume (in which the elementary mass of air ∆m is located) in front of the fragment is: Now the total compression work on interval ∆t is: The average power that corresponds to this compression at interval ∆t is: On the other hand, the equivalent force that would produce this compression power is: From Eqs (9) and (10), it follows for corrective force: In Eq. (11) state 1 is a known state of undisturbed air (before the fragment is encountered), and Ap is the projection of the fragment surface on a plane perpendicular to the direction of the velocity vector of fragment centre of mass Vf. It should be noted here that vs1 and p1 are constants in relation to given integration over the surface, whereas vs2 varies as a function of spatial coordinates. It should also be noted that the pressures in this analysis are absolute pressures.
If the equation of the adiabatic change of the gas state is used 1 Now Eq. (11) becomes: This can be written as: or in the following form: In Eq. (15) p1 is the pressure of undisturbed air (p1 = p0 = 101 325 Pa, at temperature 15 °C), and p2 is the local variable pressure (in absolute amount) along the fragment surface.
The total aerodynamic force (necessary for calculation of fragment trajectory) acting on the fragment is equal to the vector sum of force obtained on the basis of the basic aerodynamic force model [5], presented in expression Eq. (1), and corrective force obtained on the basis of the concept of loss of kinetic energy due to compression presented in expression Eq. (15).

Validation of the Model
The comparison of the results obtained by numerical simulation (using CFD software) and experimental data shows that this model of compressibility significantly reduces the relative error of aerodynamic force modelling in the relevant area of the relatively high velocity of the fragment, i.e. in the supersonic velocity range (when the compression of the air is significant).
The model is first validated using aerodynamic force data obtained by numerical simulations of airflow around a triaxial ellipsoid, shown schematically in Fig. 4. Semi-axes of an ellipsoid (Fig. 4) with which the validation of the model was performed were as follows: a = 0.034 m, b = 0.00865 m and c = 0.006 m. These dimensions were selected in order for the results to be compared also with the results of aerodynamic force acting on irregularly shaped fragment.
The method of numerical simulations of airflow around an ellipsoid consisted of the: digitalization of the body model, domain discretization (around 1,5 million cells), characterization of the resistive medium, initial and boundary conditions, solver and turbulence model selection, and aerodynamic force and moment components determination (postprocessor).
The body was considered stationary and the flow around it was analysed. The velocity vector was in the direction of axes z and y (separately, Fig. 4). The coordinate system was set in the body centre of mass.

Fig. 4 Schematic representation of the ellipsoid with which validation of the aerodynamic force model was performed
Simulations of flow over the body for seven different velocities (1, 1.2, 1.3, 1.5, 2, 3 and 4 Mach) were carried out. In numerical simulations, air is modelled as a homogeneous, isotropic, ideal gas. At the end of the domain, "Pressure Far-field" condition was used, which is often used where the compressibility is significant. The "No-Slip" condition is defined on the surface of the ellipsoid, which means that the relative flow velocity on the surface of the body is equal to zero. Boundary condition ("The Wall") is generally used in the case when the viscous effects cannot be ignored and it is relevant to most fluid flow situations. A density-based solver was selected for the simulations, where mass, flow, and energy equations are determined as the Navier-Stokes equation system in integral form for an arbitrary control volume. The Spalart-Allmaras turbulence model was used in the simulations. This physical model of turbulence has been developed specifically for aerodynamic applications and has proven to be effective for the boundary layers with high-pressure gradients, as well as it has been very effective for transonic flows around the aero profiles, including flows with significant separation of the boundary layer.
Tab. 1 shows the comparisons of the results for the aerodynamic force acting on the ellipsoid (Fig. 4) obtained by numeric simulations in Ansys Fluent with the results for force obtained using model described here. The results show that there is no significant deviation (rel. diff. of 10.9% to 13.8% for the flow in the direction of y axis, and from 0.4% to 10.6% for the flow in the direction of z axis). The model is then further validated using also the experimental data [8][9][10] for the drag coefficient CD of a sphere (Fig. 5). Sphere was used because it was easy to determine its projected surface area, and some experimental data were also available for sphere. The drag coefficients CD of a sphere as function of Mach number were determined using Eq. CD=F/qAp, where F is the aerodynamic force component in the direction of the velocity vector (drag force), q is the dynamic pressure (0.5ρV 2 ) and Ap the projected surface of the sphere. The results show an excellent agreement between the results obtained using the model presented here and the experimental data. Relative differences between the experimental data and this force model which takes into account air compressibility were below 11%. The peaks of CD values in transonic zone show that shock waves are indirectly taken into account in the model presented here. Figs 6 and 7 show the pressure distribution on the ellipsoid (determined from the model described here) and from the results of numerical simulation (the flow over an ellipsoid in Ansys Fluent, for Mach number of 1.5). The static pressure shown in these figures is determined indirectly, based on the determination of the elemental aerodynamic force over the surface of the model. The elemental aerodynamic force is obtained on the basis of the change in the amount of fluid momentum, and this value of force is corrected due to the compressibility (and consequently the shock wave) effects, on the basis of the energy equation or energy consumption on the compression of the local fluid near the surface of the body. This process is considered to be relatively fast, so there is not enough time for heat transfer, and the adiabatic (isentropic) equation is used. [3][4][5] of drag coefficient CD for sphere with results obtained using aerodynamic force model presented here Fig. 6 shows that the pressure (obtained in the model presented here) is not constant over the surface of the body, which is similar to the numerical simulation results (Fig. 7). Regarding the values of overpressure, results in Fig. 6 show that, for the results obtained from the model, the overpressure varies along the body from 0.5 bar up to values of around 3 bar. Numerical simulation results (Fig. 7) show that the overpressure on the body varies from around 0.2 bar up to maximal 2.9 bar. The most important values in these figures are the maximal overpressure values, and the relative difference between the maximal values from the model and simulation is around 3.5%.

Fig. 5 Comparison of experimental values
Similar results for maximal overpressure could be obtained by applying a 1D analytical normal shockwave model [11], where for 1.5 Mach flow velocity, the maximal overpressure behind the shockwave is around 2.5 bar. However, we should bear in mind that the analytical shockwave 1D model [11] cannot be adequately applied to the body, because the pressure would be identical all around the body (since the flow velocity is not changed along the body in the 1D model). As noted in reference [11], for a full shockwave solution of real problems, one has to use numerical simulationcomplete Navier-Stokes equations, together with continuity equation, and energy equation (where compressibility is dominant). The underpressure zone behind the body was not taken into account in the model because it is assumed that the variations in pressure behind the fragment do not significantly affect the overall aerodynamic force acting on the body during the flight at relatively short distances for the following reasons: • this underpressure accounts for a small proportion of the total pressure distribution (underpressure/relative pressure behind the body may have values up to 1 bar maximum, which is negligible comparing to the overpressure values at the frontal surface of the body -surface exposed to incoming flow (usually by an order of magnitude greater; i.e. at 3 Ma overpressure in front of the body is about 11 bar -so the maximum relative difference between these pressures is < 10%), • given that the total aerodynamic force is calculated, the consideration of this underpressure behind the body cannot significantly influence the total amount of aerodynamic force, as an integral of the elemental forces due to pressure, • also, when one is mostly interested in the movement of the body up to 50 m from the centre of the explosion (for lethal zone calculation), high velocities and aerodynamic forces occur, and the dominance of the front overpressures relative to the pressures behind the fragment remain.

Application of the Model to an Irregularly Shaped Body
After the verification of the model, the results of the model for an estimation of aerodynamic force were compared with the results obtained by CFD numerical simulations for aerodynamic forces acting on an irregularly shaped fragment with jagged surface (Fig. 10). The methodology of numerical simulations for the fragment was similar as the ones used for an ellipsoid. Detailed explanations on the simulations can be found in the authors' earlier papers [7,12]. In this case the (numerical simulations) flow was directed towards positive and negative y and z axes (Fig. 10).

Fig. 10 CAD model of the body with an irregular shape (i.e. HE projectile fragment)
Tab. 2 shows the comparisons of the results for the aerodynamic force acting on the fragment (Fig. 10) obtained by numeric simulations with the results for force obtained using the aerodynamic force model presented here (correction due to compressibility of air).
The results in Tab. 2 show that there are no large deviations (relative difference of 5.3% to 23.6%) for the flow in the direction of z axis for this particular shape of the fragment (the flow towards the largest exposed area of fragment). In case of the flow towards the smaller exposed surface of the fragment (i.e. in the direction of the positive and negative y axis; Tab. 2, Fig. 10), the relative differences are somewhat larger (23.3%-59.6%). In addition, for different directions of the fluid flow, even in the same direction (same axis), the results may be significantly different (Tab. 2). It should be noted that every possible fragment (even with the same general dimensions) will have slightly different values of these forces because of the stochastic nature of the natural fragmentation process. The results from the numerical simulations for real fragments were also compared with the results of the aerodynamic force acting on fragment -determined using the "classical" model (F = 0.5ρACDV 2 ) where force is determined via drag force coefficient (to establish the order of errors), used by most authors. In this analysis, the fragment was approximated with rectangular parallelepiped (Fig. 11), as this approximation is frequently suggested by other authors. The largest projection of the exposed surface of the fragment was represented with Amax = ab (Fig. 11), the average projection of the exposed surface of the fragment is assumed Aaver = ac; smallest projection of the exposed surface of the fragment, in this case, is Amin = bc (Fig. 11). These projected surfaces were considered so the results can also be compared with the results obtained with the model presented here (shown in Tab. 2). Dimensions of the fragment were the same as in Fig. 10 (a = 34 mm, b = 8.65 mm, c = 6 mm). Air density was taken as standard value at sea level ρ = 1.225 kg/m 3 .

Fig. 11 Fragment approximated with a rectangular parallelepiped
Due to the lack of more recent data, experimental values (CD_max, CD_aver) from the reference [8] were taken as approximate CD values for this fragment (see Fig. 12). It is assumed that the values of CD_max correspond to the case when the fragment was exposed to the flow with the largest surface projection (Amax, Fig. 11), and CD_aver for the case when the fragment was exposed to the flow with an average surface projection (Aaver, Fig. 11). Fig. 12 Experimental data for drag force coefficient CD for fragments [8] Tab. 3 presents the comparison of the results for the aerodynamic force acting on the fragment (from Fig. 11), obtained by numerical simulations and using the "classical" model (force determined via drag force coefficient), for two directions of the flow (towards positive/negative y and z axis, Fig. 11).
Relative differences between the results for aerodynamic force obtained by the "classical" approach and numerical simulation (for a fragment) for the flows towards the y axis were from 192.8% to 386.3% and for the flow towards the z axis from 128.5% to 206.5% (Tab. 3). These errors can be larger or smaller depending on the adopted CD values, as well as the adopted values of the reference area (i.e. projected area of a body -projection of an exposed area of a body onto the plane perpendicular to velocity vector [2]).
As it can be seen, using the "classical" approach for the determination of aerodynamic force acting on a fragment can give significantly larger errors (192%-386%) than using physical model presented here (where rel. differences were < 60%), which can lead to erroneous (less accurate) results for fragment trajectory.
Tab. 3  Overall, given that model presented here significantly optimizes the processing and duration of the calculation, and given the simplicity of the model, the results are satisfactory. A significant advantage of the model presented here, compared to the classic approach, is the possibility of generalization by identifying certain parameters of the fragment (a, b, c semiaxes of a triaxial ellipsoid, or geometric relations a/b, a/c), and the generalization of on arbitrary fragment, which cannot be done through a classic approach using numerical simulations (Navier-Stokes equations in full form), regardless of their greater accuracy.
The developed model for assessing aerodynamic force and moment can be used, together with a developed model for estimating the projected surface of the fragment, to predict the elements of the trajectory of the fragments.

Conclusion
The model of compressibility was developed and applied to a generalized model for prediction of aerodynamic forces acting on the irregularly shaped body (such as primary fragments of naturally fragmenting high-explosive warheads, secondary projectile debris, fragments made from structures when strong tornadoes are formed, or any body with an irregular shape).
Changing the gas state at the compression of air in front of the fragment, formed immediately after the disintegration of a HE projectile, is assumed to be adiabatic, since the air compression process during the motion of the fragment is an extremely fast process (relative to the heat transfer process).
The equation of the state of the ideal gas is adopted, although it is possible to adopt some of the models of the real gas and carry out the same analysis.
Analysis of results and comparison with results obtained by simulation (using CFD software) and experiments show that this model significantly reduces the relative error of aerodynamic force modelling in the relevant area of the high velocity of the fragment, i.e. in the supersonic velocity range when the compression is particularly significant in the air.
This compression model, together with the authors' previous work can be used in a generalized 6DOF model for the estimation of all relevant kinematic parameters (trajectory, velocities, orientation) of an irregularly shaped body.