Instability estimation of irregularly shaped bodies moving through a resistive medium with high velocity

Article history: Received 8 June 2017 Received in revised form 27 July 2017 Accepted 29 July 2017 The purpose of the paper was to show an idea how numerical simulations of flow around a stationary irregularly shaped body can be used to estimate instability of the body during a real-world motion of such a body (e.g. a metal fragment). To the best of our knowledge, there is no evidence that such an analysis is available in literature for irregularly shaped bodies. The novelty is in the introduced method for the stability analysis and the fact that a realworld fragment shape was digitized and used for the analysis. However, the disadvantage is in necessity that real fragments need to be scanned and digitized for the analysis, but the future work should give improvements in this direction. The focus was on the rotational part of the motion, particularly on obtaining the period of the motion when the body rotates, but the solving for angles of rotation was not the objective. We showed an idea on how to estimate the period of instability when continuous rotation occurs after the initial projection of the fragment. We assumed that relatively high angular velocity occurs at the initial condition (initial projection of the fragment), which provided an opportunity to assume further that the axis of rotation remains unchanged during the motion. By analyzing the kinetic energy of rotation, we estimated the period of body rotation until it reached a stable orientation during the high velocity motion. To employ this approach that uses the mechanical energy, it was necessary to obtain the work done by the (aerodynamic) moments of resistance forces about the center of mass. These resistance (aerodynamic) moments were obtained for various orientations of the body using simulations of fluid flow around the real geometry of the body, which was obtained by scanning a real-world fragment, digitizing it, and importing it in a CAD software, which provided the inertial properties through moments of inertia. At each rotation, the kinetic energy of rotation is dissipated through work of the aerodynamic moment which was the basis for calculation when the body takes a steady orientation for the rest of the motion.


Introduction
*Aerodynamic forces and moments are the results of the distribution of pressure and tangential stress on the surface of the body moving through the atmosphere. The main goal of aerodynamics is to determine the pressure and the tangential stress for given body shape and free flow conditions, and with the help of obtained values-determination of aerodynamic forces and moments (Anderson, 1991).
The total effect of pressure and tangential stress, integrated across the entire surface of the body, is the resultant aerodynamic force of ⃗ and the resultant moment ⃗⃗⃗ . The resultant force acts in the center of pressure, whereby in general case the center of mass of a body is not in the center of pressure. The resulting moment acts in the center of the mass of a body and causes its instability regarding the body orientation during the highvelocity motion.
Motion of irregular shape bodies at high velocities through the atmosphere is complex and generally difficult for analysis.
Examples of such bodies include fragments of high explosive devices (projectiles, bombs, improvised explosive devices), shrapnel generated by fractures of various structures due to the effects of strong storms, meteoroids (smaller stones or metal bodies from the space), comets (icy bodies from the outer part of the solar system), asteroids (larger bodies of different structures that also come from the solar system).
While fragments of high explosive devices move at initial velocities of 2 to 6 Mach (dimensional quantity representing the ratio of flow velocity past a boundary to the local speed of sound), meteoroids, comets and asteroids move at considerably higher velocities when entering the atmosphere (about 50 Mach and more).
Irregular shape bodies that are moving with high velocities are characterized by viscous and compressible fluid flow, dominant resistance due to pressure drag, shock waves, turbulent flow, and by a significant separation of the boundary layer from the surface of the body during movement (Buresti, 2000).
A review of the literature showed that this topic was not adequately researched. Regarding the problem of stability of the body, researchers mainly deal with axisymmetric bodies (projectiles, rockets, truncated cones, cones with cylindrical bodies, etc.). There is no evidence that the problem of the stability of the orientation of the body with an irregular shape is analyzed and adequately solved.
In this paper, an analysis of initial instability of irregularly shaped bodies with high velocities is performed. Analytical methods (physical model) and numerical simulation methods were used.

The physical model
Initial instability of irregularly shaped bodies (an example given in Fig. 1) with high velocity occurs due to the eccentric action (does not act in the center of mass) of the force that causes their movement.
An example is the fragments created by detonation of high explosive projectiles. Due to the detonation of explosive charge in the projectile, pressure increases within the body (up to 400 000 bar), in a very short period of time (order of s).
The effect of the detonation products on the body of the projectile can be considered as an impulse load. The body of projectile begins to expand and cracks begin to emerge. When the internal pressure, due to the expansion of the detonation products, exceeds structural resistance of the body, the fragmentation of projectile body leads into a large number of irregularly shaped bodies, each body having a different shape (stochasticity of the process). Due to the nonuniform effect of the product of detonation pressure force, the fragments have always an initial angular velocity which generally has its axis of rotation oriented arbitrarily in the space.
In order to determine the minimum value of the initial angular velocity of the irregularly shaped body, which ensures its continuous rotation, it is necessary to assess the values of the aerodynamic moment acting on the body in the initial phase of movement, and accordingly from the expression of the body kinetic energy of the rotation to determine the required angular velocity threshold.
During the motion, the kinetic energy of rotation dissipates through the work of ⃗⃗⃗ , until it meets the threshold equal to the total work of the aerodynamic moment on a single rotation. The idea is to determine the rotational kinetic energy threshold by numerical integration using aerodynamic moment values obtained using simulations for various fragment orientations and use the threshold to estimate the time frame of continuous rotation before the fragment takes a stable orientation during the motion.  Fig. 1 shows the aerodynamic forces and moment acting on the irregularly shaped body (with a center of mass cm), and the vector of angular velocity the direction of which is located arbitrarily in the space.
Work of moment as a consequence of aerodynamic forces can be written as: (1) is: while relation for unit vector for angular velocity axis is: In (3) α, β and γ are angles between the direction of angular velocity vector and axes x, y and z, respectively. For angles α, β and γ following equation holds: cos 2 + cos 2 + cos 2 = 1.
Angles α, β and γ are determined using expression: where angular velocity magnitude is: Based on (1)-(3), mechanical work of the aerodynamic moment can now be written as: or approximately using expression: where S = (φ1-φ0)/φ. Angles φ1 and φ0 are angles of two orientations when aerodynamic moments are negligible i.e., if a coordinate system is set so that φ0 = 0, then φ1 = . An analogy with the rotating pendulum ( Fig. 2) can be used here. The pendulum in Fig. 2 is schematically shown in two different positions where the moment is zero, for the following angles: φ0 = 0 and φ1 = . Thus, the integration angle is φ1-φ0 = . However, in the case of the fragment motion, neither ⃗ , nor ⃗⃗⃗ is constant in direction or magnitude. In order to estimate the work of the aerodynamic moment from the expression (8), it is necessary to know the moment components for various body orientations and the position of angular velocity axis (using angles α, β and γ). The aerodynamic moment components are determined using numerical simulations for different body orientations, for angles of 0° to 180°, with an increment of rotation of 15°. When turned into radians, the initial orientation of the body corresponds to the angle φ0 = 0, and the last orientation of the body corresponds to the angle φ1 = , similar to the schematic diagram of the pendulum in Fig. 2.
The position of angular velocity axis ( Fig. 1) depends on angles α, β and γ, and they can be assumed in calculations. More precisely, two of these three angles can be assumed, and the third is calculated from the expression (4).
The kinetic energy of the body generally consists of translational and rotational components. In this case we are interested in the kinetic energy of rotation of the body: where is the moment of inertia of the body for the instantaneous axis - in the direction of the angular velocity. Since the position of the instantaneous axis constantly changes, I is variable, so it is more practical to express the kinetic energy using the moments of inertia for axes that are bound and together with the body move. If the origin of such a moving system, Oxyz, is adopted in the stationary point O, then the expression (9) can be written in the form: where x, y and z are projections of the current angular velocity on Oxyz system, while Ix, Iy, Iz, Ixy, Ixz, and Iyz are moments of inertia of the body for the Oxyz system. Expression (10) can also be written in the following form: If axes Oxyz are principal axes of inertia, then moments Ixy = Ixz = Iyz = 0, so (11) can be reduced to: = 1 2 2 ( cos 2 + cos 2 + cos 2 ).
The moments of inertia of the irregularly shaped body for coordinate axes x, y and z (bound to the body) can be determined in the CAD software, based on the three-dimensional body model.
The body will rotate continuously during the flight if the following condition is met: In (13), and are given by (8) and (11). Since part of the kinetic energy is dissipated with every rotation, the condition (13) should be satisfied for each subsequent rotation. Therefore, the condition of continuous rotation during flight is: ½ 2 ( cos 2 + cos 2 + cos 2 − 2( cos cos + cos cos + cos cos ))> (14) where α, β and γ can be considered constant for a single rotation, and S = (φ1-φ0)/φ.
The minimum angular velocity required for continuous rotation can be finally written as: It is possible to estimate the minimum number of revolutions of the body during the flight due to the initial kinetic energy of rotation: where Ek rot is initial kinetic energy of the rotation, and the W rot work of a moment for rotation between the two stationary orientations, while the symbol [ ] indicates that the first lower whole number is taken.
Here it is assumed that the moment function changes slightly between the two adjacent rotations so that the work of moment remains approximately the same.
After the decrease of the kinetic energy of rotation below the threshold of energy for a full rotation, the body can reach the state of oscillation around the center of mass, with the oscillation amplitude decreasing due to the dissipation of kinetic energy. Parallel with this, rotational vibration around other axes can occur due to the aerodynamic moment, so the body can still tumble and rotate within a specific spatial angle during the flight.

Numerical simulations
Instead of simulating the fragment spatial motion, the fragment was considered stationary and the flow around it was analyzed. The reason for this approach is the extremely high requirements for large solution domain, and therefore, high requirements for computer resources if the spatial motion of the fragment are considered.
In order to estimate the components of the aerodynamic moment (required in expression 8) for an irregularly shaped body, numerical simulations of flow around an irregularly shaped body were performed in the Ansys Fluent CFD software package.
Methods of numerical simulations, using computer fluid dynamics, are an important aspect of modern research as they are complemented by experiments and analytical models, thus reducing total time and labor costs.
Generally speaking, basic flow equations are represented by the equation of continuity, the momentum equation and the energy equation (Anderson, 1991).
The method of numerical simulations of air flow around an irregularly shaped body consisted of the following: a. Digitalization of the body model using real fragments b. Physical domain discretization c. Initial and boundary conditions d. Characterization of materials e. Solver and turbulence model selection f. Aerodynamic moment estimation (UDF Script)

Digitalization of the body model
The irregularly shaped body with which the numerical simulations were performed for different orientations is presented in Fig. 3.
Since the 3D scanner was not available, a body was digitized using CAD software. A 3D scanner will digitize the body more precise, but CAD technique can also very useful when scanner is not available.
The three dimensional model of the body (Fig. 4) is 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.

Physical domain discretization
Since these bodies have a stochastic and irregular shape, it would be very difficult to define their geometry using point by point technique in one of the older preprocessors such as Gambit.
For this reason, Ansys Design Modeler was used, which enables the introduction of finished CAD models.
The mesh of the numeric model was unstructured (Fig. 5), consisting of polyhedral elements: 625689 cells. 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.

Initial and boundary conditions
In the simulations, the velocity of fluid flow around the body was 3 Ma. This speed corresponds to the initial velocity of the fragments generated by the detonation of high-explosive projectiles.
The velocity vector ( Fig. 6) was directed in the positive direction of axis X of the coordinate system set in the body center of mass. 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 document, together with the model). Numerical simulations for 14 body orientations were performed: 0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 165 and 180. Fig. 6 shows the schematic position of the body in numerical simulations.
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, 2011) are used: where are: 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, so called 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 case when the viscous effects cannot be ignored and is relevant to most fluid flow situations (Fluent, 2011).

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, 2011): where: pat -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, 2011): 3 2 ( 0 + ) 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, 2011).

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, 2011): pressure based and density-based solver. According to the recommendation (Fluent, 2011) for use with compressible flows, a densitybased 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 densitybased 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, 2011) 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 of 3 Mach, the Re number exceeds 10 5 , which means that the flow around the body with this velocity is prevalently turbulent, although the time of movement of this body through the atmosphere is relatively short.
According to the recommendations (Fluent, 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 (Pope, 2000).
The ever-increasing 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, Rumsey, 2012.

Aerodynamic moment estimation (UDF Script)
A program (in C programming language) was written that determines the aerodynamic forces and moments 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: 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: The components of moment ⃗⃗⃗ 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. 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 command (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.

Validation of numerical model
Validation is a process of estimating errors and uncertainties of a numerical model by comparing obtained results with available experimental data. Validation is final stage of model checking process, which determines degree of agreement of the adopted model with the real physical phenomenon investigated in the analysis (Tu et al., 2008). The adopted numerical model in operation was verified with experimental data for the cube. Verification was performed by comparing the coefficient of drag in the simulation and for available experimental data. The cube is chosen because experimental data on drag coefficients are available for it, and its geometry is known.
For comparison, the experimental data (Schamberger, 1971) were used for the different positions of the cubes (Fig. 7) during the airflow around them: when the flow is facing a cube side (flat-on orientation) and when the flow is facing the cube edges (edge-on orientation).
Based on the obtained results on forces, moments and exposed surfaces (reference area), for different values of Mach numbers, the values of the CD (cube drag coefficient) is determined using formula (Anderson, 1991): = where: FD -drag force in the direction of velocity vector, S -reference area (projection of total body surface perpendicular to the velocity vector, also determined in UDF program). Dynamic pressure is defined by expression = 0,5 2 .  (Schamberger, 1971) Numerical simulations were carried out in such a way that the same conditions apply to experimental tests by Schamberger. The simulations were performed in two different positions of fluid flow over cubes, with the same Ma numbers for which the Schamberger tests were performed. 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. 8 gives a comparison of numerical simulation results with experimental data (Schamberger, 1971) for the cube for flat-on and edge-on orientation of cubes. The difference between values of the CD from simulations and experiments are smaller than 10%.  (Schamberger, 1971) for a cube

Results and discussion
In Fig. 9 and Fig. 10, the pressure field and streamlines are represented for flow around irregularly shaped body.
The pressure field (Fig. 9) around the body is complex and nonsymmetrical because the body is irregularly shaped with a large number of edges and dents. An overpressure zone is present in front of the body, while an under pressure zone appears behind it. Of great importance to the movement of the body of this velocity through is the occurrence of shock waves. Namely, the velocity of 3 Ma implies pronounced shock waves around the body. In this case, the shock wave is curved and usually called bow shock wave since the body is blunt. Generally speaking, the shape of the bow shock wave and the complete field of flow between the bow waves and the body depend on the Mach number and size and shape of the body. Bow shock waves significantly increase the drag of supersonic velocity bodies, which is, for example, used in the design of the Apollo capsule that required a high drag to prevent the capsule from burning when entering the atmosphere. Solving problems with bow shock waves is complex and regularly requires the use of numerical simulation methods (Anderson, 1991).
Regarding streamlines (Fig. 10), the relative movement of fluid from the region of high pressure to the low pressure region in combination with the movement of the fluid in the direction leads to the formation of a vortex flow behind the body, dissipating energy which significantly increases the drag. Vortices are usually in mutual interaction, they are movable and can exchange energy.
It is important to notice that the simulation considers that the fragment is stationary, but the surrounding medium is in motion. It is not exactly the same situation as in the case when the fragment is in motion through the stationary medium. However, the relative motion of the surrounding medium with respect to the fragment is the same. This means that the pressure field and the streamline field will be the same as if the fragment is in motion, meaning that all values of the surface forces and resulting moments can be used to estimate motion of the fragment through the stationary medium.
This approach provided the possibility to use relatively dense mesh in the relatively small solution domain around the fragment. Otherwise, a large solution domain would be needed, which would result in a very coarse mesh that cannot give sufficiently accurate results.
In order to estimate minimal angular velocity required for continuous rotation, the results of numerical simulations and analytical expressions (1)

Flat on flow
Edge on flow -(16) were used. The following assumptions have been adopted:  The analysis refers to body with mass of 52.5g because simulations have been performed with this body.  Aerodynamic moment components are determined by numerical simulations for several orientations of a body (0 -180, with a 15 increments, Table 2).  Coordinate axes bound to the body are principal axes (these are determined in CAD).  Moments of inertia for axes bound to the body are determined using CAD software (Table 1), and based on them moments of inertia for body for current axis of rotation  are determined.  It is assumed that the angles of the current or rotation of the body are the same (α = β = γ). Based on the expression (4), that angle is 54.7°. These angles can be set arbitrary, with the two angles being assumed, and the third being defined by the expression (4). Since the orientation of the initial rotation axis can be any axis through the center of mass, due to the nature of the impulse forces that cause the initial rotation of the body, then this analysis is not limited to any particular axis of initial rotation. However, the numerical simulation has specific given initial conditions.  It is assumed that the initial velocity of the body is 3 Mach (equivalent to initial velocity of fragments from some HE projectiles), so the moment component values obtained for velocity 3 Mach are taken.  It is assumed that, due to the intense gyroscopic effect, there will be no significant rotation around the other axis that significantly deviates from the axis of initial (dominant) rotation. Beside the significant gyroscopic effect, this assumption is further justified by the fact that the time frame of the body movement, before it hits an obstacle, is relatively short, so that the body does not have enough time to develop any significant rotation about any axis different than the axis of the initial rotation.
In Table 1, moments of inertia for principal axes (coordinate system bound to the body) of the body with a mass of 52.5g and moments of inertia for the current axis of rotation  are given.  Work of an aerodynamic moment W rot for a rotation between the two stationary orientations is determined using the expression (8). Table 2 shows data for the aerodynamic moment components and the calculated work performed by these moments for the angle of integration from 0 to .
As can be seen from Table 2, single work of moment (for rotation angle 15) can be negative or positive. It is negative when the aerodynamic moment contradicts the initial rotation of the body, and positive when the moment further enhances the initial rotation of the body. Based on the data of total work of moment between the two stationary orientations (Table 2), it is now possible using the expression (15) to estimate the minimum initial angular velocity required for continuous rotation of the body.
For this particular case (mass of the body was 52.5g), obtained minimum angular velocity required for continuous rotation of the body is 112,4 rad/s Pressure around the body (17,9 rev/s). Thus, if the initial angular velocity of the body is greater than 18 rev/s, the body rotation (instability) will occur immediately, since the body will rotate through the stabilization zone (similar to the pendulum in Fig. 2 which passes through the stabilization zone if the initial angular velocity magnitude is sufficient such that rotational kinetic energy is higher than the absolute amount of work done by the moment about the joint). It is important to note that research on estimate of initial angular velocity of fragments from detonating warheads is scarce. Recently Norwegian researchers (Moxnes et al., 2017) made numerical simulations of HE warhead detonation using IMPETUS Afea program, where they concluded that in their case (fragment with mass of 0,51 g) initial angular velocity of fragment was in the order of 23 000 rad/s (3650 rev/s). So significant angular velocity can be imparted to fragments in the initial stage of their flight.
In order to estimate the minimum number of body revolutions due to the initial kinetic energy of rotation, the initial angular velocity of the body of 200 rev/s was assumed. Initial angular velocity estimate given here is rather conservative since, as we mentioned, there are reports that initial angular velocity of HE fragments can be in order of around 23 000 rad/s (Moxnes et al., 2017).
Using Eq. 16, where the initial kinetic energy of rotation (Eq. 12) and work of the moment during the rotation between two stationary orientations (Eq. 8) are determined, we obtained minimum 125 rotations of the body around the center of mass due to the initial angular velocity. Since initial angular velocities of real fragments are much higher (Moxnes et al., 2017), it is reasonable to conclude that fragments could rotate due to the initial disturbance (initial angular velocity) during its whole trajectory, since travel time of these bodies is usually very short due to their high velocities.
It is interesting also to analyze the time at which the body will rotate (rotation period) and at what distance from the center of the explosion. Rotation period can be determined using known formula: In the case of an assumed initial angular velocity of 200 rad/s, the body rotation period is 31.4 ms. If the assumed initial angular velocity is 2000 rad/s, the body rotation period would be 3.14 ms. The distance of a body for a full rotation can be determined using the expression: assuming that the horizontal component of the initial velocity is vx0. For an initial velocity of 1000 m/s and an initial angular velocity of 200 rad/s, the (approximate) distance that the fragment travels for one full rotation is around 31,4 m. If the initial angular velocity was, for example, 2000 rad/s, then the distance traveled is 3,14m during a single rotation. One of the possible criteria for comparing the translational and rotational component of the motion of the fragment is the ratio of the kinetic energy of the translation and the kinetic energy of the rotation: where I is moment of inertia around arbitrary axis.
If the values of the kinetic energy of the translation and the kinetic energy of rotation are compared for a given fragment (m=52,5g, v0 = 1000 m/s, I = 11,3472 kgmm 2 , 0 = 200 rad/s) than = 26250J, and = 0,223J, which means that the kinetic energy of translation for this case is 115668 times greater than the kinetic energy of the rotation of the fragment. If the initial angular velocity of the fragment is, for example, 2000 rad/s, then the ratio h = 1157.
If the body still does not hit the obstacle (the flight to the obstacle is usually so short that the body does not at all reach a stage where its angular velocity will fall to a negligible value), after a body tumbling due to a dominant initial disturbance (initial angular velocity), temporary body stabilization can occur with oscillation about a stable orientation (oscillation around the center of mass), or (most likely) further rotation and tumbling of the body in the case of extended time intervals due to the translational velocity drop.

Conclusion
The study evaluated initial instability of the irregularly shaped body by analyzing the minimum angular velocity required for continuous rotation during motion. This angular velocity is estimated based on moment components (obtained by numerical simulations) that acts on the body with an irregular shape and the suggested method that takes into account the work of the aerodynamic moment and the initial kinetic energy of the rotation.
After the decrease of the kinetic energy of rotation, the body can get into the state of oscillation (around the center of mass) within the stability zone, within which the amplitude of the rotational oscillation decreases due to further dissipation of the kinetic energy, resulting in the final stable orientation of the body during the flight.

Future work
Further steps in the path of this work should contain the following:  Stability zone determination regarding fragment orientation during motion through resistive environment,  Determination of the center of pressure and direction of the resultant force vector for irregularly shaped bodies,  Determination of stability conditions on specific fragment orientation where the aerodynamic moment is zero,  Estimation of the fragment trajectory using aerodynamic forces for real fragments,  Numerical simulation of irregularly shaped bodies using 6DOF solver.