Prediction of the trajectory of an irregularly shaped body moving through a resistive medium with high velocities

Article history: Received 29 June 2017 Received in revised form 4 September 2017 Accepted 5 September 2017 A model was made for prediction of the trajectory of an irregularly shaped body moving through a resistive medium with high velocities, using data for aerodynamic forces obtained from numerical simulations. For a different orientation of the body with respect to the velocity vector of the center of mass, the aerodynamic resistance force is different for two reasons: the exposed surface area is different and the shape is different. In this regard, 216 numerical simulations of airflow around of the body of an irregularly shaped body in different orientations were carried out, for one full rotation (around one axis of rotation) of the body, with angular increments of 15 (0 360), for the following velocities: 0.6, 0.8, 1, 1.2, 1.3, 1.5, 2, 3 and 4 Mach. The outcome of these simulations is the resistance forces and aerodynamic moments as the result of motion of the body in various directions relatively to the body. After the simulations had been performed, the results of the resistance forces and aerodynamic moments were used to integrate the equations of motion with an assumption that the irregularly shaped body had a continuous rotation all the way along the trajectory with relatively high angular velocities. With this assumption, an effective aerodynamic force was calculated which takes into consideration that the aerodynamic force varies due to the rotation of the body. The results show that the trajectory of an irregularly shaped body is curved in space because the side component of the aerodynamic force cannot be ignored because of the irregular shape of the body, which leads to significant lateral movement of the body from the initial direction of flight.


Introduction
*An estimate of trajectory and velocity during flight for irregularly shaped bodies is often of great practical interest. One example of such bodies are the fragments (with their stochastic and irregular shape) created by the detonation of high explosive projectile that move at high initial velocities (3)(4)(5), and are unstable since the beginning of movement due to initial disturbances (extremely high-pressure detonation products, of an order of 20-40 GPa, are responsible for high initial angular and translational velocity).
In Catovic (2017a, 2017b), the initial instability of such fragments and the stability zone around the equilibrium orientation was analyzed. It was found that after the decrease of their kinetic energy of rotation, the fragments may come into the state of oscillation (around the center of mass), and also there are rotations around other axes due to the aerodynamic moment, so that the body can still be overturned during the flight. All this points to the complexity of estimating trajectories of such bodies.
Estimation of the trajectory, velocity and kinetic energy of the HE projectile fragments is significant from the point of view of the lethal zone estimation around the HE projectile, but also for estimating the effective range of these fragments, which may be of practical interest in the case of an unexpected explosion of military depots.
In addition, in recent years there have been more and more attacks on civilians using improvised explosive devices, so exploration of motion, stability and kinetic energy of bodies resulting from the explosion of such devices is justified.
For the body of an irregular shape it is necessary to predict the trajectory of the center of mass; however, for each different orientation of this body relative to the airflow velocity, there is a different aerodynamic force which in general deviates from the air-flow velocity relative to the body of an irregular shape.
Since there is generally a relatively large number of body revolutions along the trajectory of the center of mass (i.e., for the fragment of the HE projectile), it is necessary to define a method for the approximate calculation of the trajectory, which will account for the rotation of that body or variation of the aerodynamic force at the rotation of the body.

Physical model
The idea is to use an effective vector of the aerodynamic force in one rotation of the body for the prediction of the trajectory, and that the trajectory is obtained incrementally, where each increment will correspond to one rotation of the body.
When the body rotates with high angular velocities, a gyroscopic effect appears, or "resistance" to the changes of the axis of rotation direction, as in Fig. 1.
The following assumptions have been adopted in the physical model:  Assuming there are relatively high angular velocities, the axis Ω-Ω ( Fig. 1) retains its direction at a full rotation of the body.  Every orientation of the body at rotation around the axis Ω-Ω is evenly represented, ie it acts in equal time during the full rotation of the body. Under the assumptions made, the effective value of the aerodynamic force can be represented in following way: where n is the number of incremental body orientations in one full rotation. The effective value of force is that force whose impulse will be the same over a time period of one full rotation: where T is a time period of a full rotation. If the integral on the right side of Eq. 2 is calculated approximatively (numerical methods) by sum: we obtain the same expression as (1) given earlier that is used as the basis for estimating the trajectory.
Generally speaking, the trajectory of the center of mass is precisely obtained using the law of the center of mass movement and the law of the momentum change for the center of mass: where ⃗⃗⃗ − is the aerodynamic moment for the center of mass.
Similar to the search for an effective aerodynamic force (1), the effective aerodynamic moment value can be written as: or, similiar as for the aerodynamic force: Depending on the magnitude of the aerodynamic moments, there will be a certain change of rotation axis Ω-Ω ( Fig. 1) at the end of each full rotation of the body. However, within one increment of the trajectory, the orientation of axis Ω-Ω is considered unchanged. Using the aforementioned assumptions and Eqs. 1, 4, 5, and 7, a basic model for estimating the trajectory can be set. Projections of the expression (4) on the inertial coordinate system axes are: where ⃗ ∥ ⃗⃗ is assumed, and vector ⃗⃗ is a unit vector of z axis. By integrating (8), (9), and (10) for the period of a single rotation, we obtain: By approximate calculation (numerically), (11), (12) and (13) become: where p denotes the end of the p-th trajectory increment or the end of the p-th rotation. By second integration of the terms (14), (15) and (16) we obtain: The force ⃗ − should be calculated for each increment of the trajectory, as the velocity of the body changes, as well as the axis of rotation Ω-Ω ( Fig.  1) relative to the body and relative to the airflow. To determine the movement of an axis Ω-Ω, a differential vector Eq. 5 would need to be solved, which in general is reduced to three scalar differential equations.
The idea is that the direction of rotation is considered constant along one increment and that the change of direction is obtained at the end of each increment, using the effective value of the aerodynamic moment and the integrated expression (5): The new axis of the body rotation for the next increment p is determined from the direction obtained from angular momentum ⃗⃗ . If the angular velocity of the initial rotation is relatively high or if the time of the body flight to the obstacle is relatively short (such as in the case of HE projectile fragments) then the Ω-Ω axis will not significantly move during body flight and can be considered constant. Then the force ⃗ − will only depend on the change in the fluid relative airflow velocity along the trajectory of the body center of mass.

Numerical simulations
In order to estimate the components of the aerodynamic force, numerical simulations around an irregularly shaped body were performed in Ansys Fluent.
 The method of numerical simulations of air flow around an irregularly shaped body consisted of the: digitalization (3D model in CAD) of the body model using real fragment ( Fig. 2),  domain discretization with around 700 000 cells (Fig. 2),  characterization of the resistive medium,  initial and boundary conditions (Fig. 2),  solver and turbulence model selection and  aerodynamic force components estimation (UDF Script).
The body was considered stationary (for a single orientation, Fig. 3) and the flow around it was analyzed. The reason for this approach is extremely high requirements for large solution domain, and therefore, high requirements for computer resources if the spatial motion of the fragment (6DOF numerical simulations) are considered.
Numerical simulations for 24 body orientations were performed for angles 0-360 with 15 angle increments. Fig. 3 shows the schematic position of the body (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. 3). The coordinate system in the initial position of the body coincides with the principal axes of inertia (defined in CAD and exported as iges format, together with the 3D body model).
For all body orientations ( Fig. 3), simulations of flow over the the body for 9 different velocities (0.6, 0.8, 1, 1.2, 1.3, 1.5, 2, 3 and 4 Mach) were carried out. Using the whole range of velocities in simulations, it is possible, together with the model, shown in Section 2, to estimate trajectory for the the body with an irregular shape. In simulations, air is modeled as homogeneous, isotropic, ideal gas with pressure-temperature dependent density , specific heat Cp, thermal conductivity k and dynamic viscosity .
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. 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 the the case when the viscous effects cannot be ignored and is relevant to most fluid flow situations (Fluent ANSYS, 2011). According to the recommendation (Fluent ANSYS, 2011) 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 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 ANSYS, 2011) is used.
According to the recommendations (Fluent ANSYS, 2011), 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. The everincreasing popularity of the Spalart-Allmaras model contributed to the rapid implementation of these models on unstructured meshes, unlike the classic aerodynamic turbulence models such as Baldwin-Lomax or Johnson-King (Anderson, 1991;Pope, 2000;Ramsey, 2012).
A program (in C programming language) was written that determines the aerodynamic forces for all three coordinate axes. For each cell on the body, the forces are determined in three directions, as follows: 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: Solver of Ansys Fluent is set by dynamically loading the developed UDF (user defined function) program and executing the commands entered into it. 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. This coefficient is determined using expression (Anderson, 1991): 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. 4 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%.  (Schamberger, 1971;Hoerner 1965) for cube

Program for prediction of trajectory of irregularly shaped bodies
MatLab was used to make a program for prediction of the trajectory of irregularly shaped bodies, using the physical model developed in Section 2 and the results of numerical simulations in Section 3.
Using the MatLab program developed in the research, the trajectory of an irregularly shaped body is estimated based on the effective value of the aerodynamic force on the trajectory increment, which means that the exact orientation of the body during the flight is not taken into consideration, and that the influence of the different body orientations is taken into account through the effective aerodynamic force. Input data of the program are (Fig. 5):  database with values of aerodynamic force for different Ma numbers (results of numerical simulations),  mass of the body,  initial velocity,  initial angular velocity, and  initial coordinates of the center of mass (initial body position).
Output data of the program are (Fig. 5):  coordinates of the center of mass (spatial trajectory), and body velocity. After completion of the calculation, the program makes the diagrams of the body velocity as a function of time and distance, as well as the visual representation of the spatial trajectory of the body.

Results and discussion
Based on the developed physical model, database (Appendix A) on aerodynamic forces acting on the body at different velocities (obtained from numerical simulations), and a developed program in Matlab (Appendix B) for the evaluation of the basic elements of the trajectory for an irregularly shaped bodies, an estimate of the trajectory (Fig. 6) for the real body of an irregular shape (shown in Fig. 3) is made. Fig. 6 shows the results of MatLab program for estimating the trajectory and velocity of the irregularly shaped body.  body mass (52.5g, numerical simulations were performed with body of this mass),  initial center of mass coordinates (0, 2m, 0),  the initial translational speed (1000m/s, 0, 0) -the initial velocity was in the direction of axis x and it coincides with a coordinate system for which the values of the aerodynamic force are obtained in numerical simulations),  initial angular velocity (0, 0, 200rad/s). Fig. 6 the following conclusions can be made:

From the diagram in
 The trajectory of an irregularly shaped body is not located in a single plane but represents a spatial trajectory (body moves in all three coordinate axes). The effect of lateral wind is not taken into account in this analysis.  Significant lateral (side) movement (in the direction of z axis) of the body is noticed, and for these initial conditions, the lateral movement was approximately 30m. Lateral movement, in this case, was solely due to the irregular shape of the body. Ratio of lateral movement (30m in direction of z axis) to straight movement (193m in direction of x axis) was around 15%.  Time of flight was around 0,42s.  The velocity of the body is decreasing relatively fast and by 0,42s it is already reduced from 1000m/s to around 270m/s. During this time the body traveled around 193m. This means that the kinetic energy of the body with an initial value of around 26KJ decreased by 0,42s to a value of around 1,9 KJ. In the case of HE projectile fragments, this kinetic energy of 1,9kJ at 193m is more than enough to kill a human since it is considered that kinetic energy of 80J (energy of lethal fragment) is enough for incapacitation of the human target (AASTP, 2010). Fig. 7 shows trajectories of an irregularly shaped body (from Fig. 3) for the same initial conditions as in Fig. 6, but now with simulations in the MatLab program being performed for several different initial velocities (4, 3, 2, 1 and 0.8 Mach). Fig. 7 shows that by increasing the initial velocity the range of the body also increases, but at the same time there is a greater deflection of the body in the direction of Z axis (sideways). Thus for example, at the initial velocity of 4 Mach, for given initial conditions, the range of the body is around 227 meters and the lateral deflection is around 49m. On the other hand, at the initial velocity of 0.8 Mach the range of the body is around 88m and the lateral movement was around 10m. The vertical axis is taken to be the y axis, such that the − term in Eq. 9 instead of Eq. 10.
These results point to the fact that the irregularly shaped body may have a significant lateral force component, which lead to a sideways motion of such irregular bodies. Results also show that fragments from high explosive projectiles can be lethal for humans even at large distances.

Conclusion
Currently, to the best of our knowledge, there is no effective exact analytical model for determining the trajectory of an irregularly shaped body. Although the general equations of motion for any rigid body are represented by (4) and (5), those differential equations of motion are not applicable for bodies that move with high velocities over relatively long distances and make extremely high number of rotations during the motion. Even if the resistance force is determined for high number of directions relative to the body, the numerical procedure would last extremely long to determine exact orientation of the body for every time step during the motion. In this paper, a model was developed for estimating trajectory and velocity of an irregularly shaped body using the data of the aerodynamic force obtained from numerical simulations. This way, the real body geometry was taken into account.
The translational component of motion was solved for incrementally; using the expanded form of (4), but the total aerodynamic force was modeled and implemented as an effective resistance force that is generated during the rotation over a full rotation. Therefore, the rotational component is not solved for exactly, but modeled with an assumption that the body will continuously rotate during the motion along the trajectory due to the initial angular velocity and the gyroscopic effect that keeps the axis of the rotation due to relatively high initial angular velocities.
In this analysis, the focus was on obtaining the body's center of mass trajectory, but since the resistance force depends on how the body is oriented, we modeled the rotational component to obtain the effective resistance force that has the same impulse over a time period of a full rotation.
This time period was taken as a time step for the rotational component of motion.
To calculate the effective force, that has the same effect as a number of the resistance forces in different directions, due to the rotation of the body along the trajectory, in total 216 numerical simulations of airflow around irregularly shaped body were performed for different body orientations, and for one full rotation (about one axis of rotation) of the body, with a 15° angle increment (0-360), for the following flow velocities: 0.6, 0.8, 1, 1.2, 1.3, 1.5, 2, 3 and 4 Mach.
The results show that the trajectory of an irregularly shaped body is spatial curve because for the irregularly shaped bodies we cannot neglect the lateral component of the aerodynamic force -which leads to the significant sideways motion of the body.
By increasing initial velocity of the body, for the same other initial conditions, the range of the body increases, but also its lateral movement.
Results show that fragments from HE projectiles (perhaps a best representative of the irregularly shaped bodies) can be lethal even at large distances.

Appendix A. Aerodynamic force data
In this appendix tables for required data regarding Aerodynamic force data in different velocities are provided.    (2,2,4), plot(xcm(1:ind),vv(1:ind)) title('Velocity as a function of distance') % printing of results in command window is the rest of the program