Determination of the center of pressure and dynamic stability for irregularly shaped bodies

Article history: Received 23 June 2017 Received in revised form 3 August 2017 Accepted 5 August 2017 Irregularly shaped rigid bodies, which move with high velocity through a resistive medium, take specific stable orientations during the motion. The paper shows physical models for estimation of the center of pressure for irregularly shaped bodies and for assessing the dynamic stability of such bodies. Using the developed physical models and the results of numerical simulations for the body of an irregular shape (a geometry similar to the fragments of high explosive projectile), a center of pressure and stability analysis was performed for the given body. It has been shown, based on rotation about one axis, that there is a zone of stability for the body in certain orientations, but to determine it more precisely it is necessary to continuously change two angles about two mutually perpendicular axes that perpendicular to the velocity vector of the center of mass. Generally, for each body, there is an orientation in which the body is in a stable zone (oscillation around the center of mass), but it is difficult to find that zone since one has to perform an extremely large number of numerical simulations for each body of a different shape.


Introduction
*Studying stability of the body moving through the atmosphere is complex, especially for irregularly shaped bodies. In the analysis of the stability of the body, one of the tasks is to determine the position of the center of pressure relative to the center of mass cm. The total aerodynamic force acts at the point called the center of pressure cp. The effect of force at one point can be replaced by the action of the same force in the center of mass and moment of force with respect to the center of mass. The center of pressure is the point where the total aerodynamic moment is zero.
The magnitude, direction, and position of the center of pressure depend on the shape, body dimensions, orientation to the air flow and the characteristics of the air flow (density, velocity, degree of compressibility, etc.). When the attack angle of the body changes, the pressure field around the body changes. Because of this, the center of pressure changes also with the change of attack angle (change of body orientation).
Stability is the property of a body that causes it when disturbed from a condition of equilibrium or steady motion to develop forces or moments that restore the original condition. Generally, we distinguish between static and dynamic stability, with subgroups: stable, neutral and unstable. Static stability refers to the ability of the body to return to its original position or orientation, i.e. to maintain its original orientation without any other dynamic movement. Dynamic stability refers to the ability of the body to return to its original position, i.e. to maintain its original orientation through the interaction of disturbed movement with other movements of the body.
Static stability depends on body design and there are three types of static stability ( Fig. 1):  Stable -the body tries to maintain itself in the initial conditions,  Marginally stable -the body remains in a disordered state, and  Unstable -the body strives to move away from the initial conditions in the direction of the disorder.
Depending on the relative position of the center of pressure (cp) and center of mass (cm), there are the following cases ( Fig. 2):  cm in front of cp (stabilizing moment -statically stable),  cm behind cp (destabilizing moment -statically unstable), and  cm and cp coincide -a neutral case. The assessment of the center of pressure and the stability of the irregularly shaped body that moves through the atmosphere is complex, and in the literature no research has been found that deals with this problem.

Determination of the center of pressure
In order to estimate the position of the center of pressure for an irregularly shaped body, it is necessary to know the values of the resulting aerodynamic force and moment and their components. For an irregularly shaped body this is possible using numerical simulations because simulations take into account their real geometry. Fig. 3 shows the center of pressure cp and body center of mass cm for an irregularly shaped body. The same image shows the resultant aerodynamic force vector ⃗ and an elementary aerodynamic force vector ⃗ on the body element dA. These vectors are, in general, arbitrarily oriented in space. Fig. 3 also shows the direction of the force vector those changes during the movement, which means that the position of the center of pressure changes also during the movement. Velocity rel v is the relative fluid velocity (relative to the body). In general, vectors ⃗ and ⃗ are not parallel. Fig. 4 gives a schematic representation of the aerodynamic force and moment (mutually perpendicular vectors) on the irregularly shaped body, as well as the radius vector r of the elementary surface dA. The total moment of aerodynamic forces for the center of mass is: (1) To determine the center of pressure for the irregularly shaped body let us consider Fig. 3. The aerodynamic force moment ⃗⃗⃗ can be defined for the arbitrary selected coordinate system by a cross product: where vector ⃗ (Fig. 3) is Vector of moment ⃗⃗⃗ is perpendicular to the plane made by vectors ⃗ and ⃗ . In (3) ⃗ * is the position vector of the center of pressure cp*, ⃗ is the position vector of the center of the mass of cm relative to the arbitrarily selected coordinate system. If the position of the center of pressure is determined relative to the center of mass (where the coordinate system is most often defined), the position vector ⃗ in (3) can be neglected.
The cross product in (2) is reduced to a system of equations with three unknowns, vector ⃗ components dx, dx and dz. The magnitude of vector ⃗ represents the distance between the center of pressure and the center of mass of the body: The system of equations from (2) has a nontrivial solution (singular matrix) since it is linearly dependent because the center of pressure can be anywhere in the direction of aerodynamic force ⃗ (Fig. 3), i.e. any position of the center of pressure in this direction satisfies the condition (2).
In order to solve the system of equations (2), it is necessary to define an additional condition, i.e. additional equation. This equation is obtained from the condition that the vector ⃗ is perpendicular to the direction of force ⃗ (the scalar product of two perpendicular vectors is equal to zero), which can be written: Now, the final system of equations can be written as: The system of four equations (6a-6d) with three unknowns (dx, dy and dz) is linearly dependent and could be unambiguously solved by the substitution method if there were no numeric errors (in numerical simulations) manifested by a small deviation of vectors ⃗ and ⃗⃗⃗ from their mutual orthogonality. After each simulation has been completed (for each orientation of the body for given Mach number), the orthogonality condition of these two vectors should be satisfied, that is, their scalar product should be zero. Only when the condition of orthogonality for vectors ⃗ and ⃗⃗⃗ is satisfied then any combination of three of the four equations (6a-6d) will give a vector ⃗ , with equation (6d) being included.
In case when the condition of orthogonality for vectors ⃗ and ⃗⃗⃗ (obtained via simulations) is not satisfied, optimization of system of equations (6a-6d) should be done using optimization function S: The method used in Eq. 7 is similar to the least squares method but is not the same as this is not interpolation or approximation of the function. From (7), the parameters dx, dy and dz (vector components) in which the function S achieves its minimum could be determined. The conditions of the local extreme (using partial derivatives) for function S are: where we now get a new system with three equations and three unknowns: After removing repeating parts, the above expressions are reduced to the following system: where 2 = 2 + 2 + 2 (squared magnitude of total force vector ⃗ ), and ( ⃗ × ⃗⃗⃗ ) , projections of the cross product ⃗ × ⃗⃗⃗ onto axes. Finally, (10) can be given in a compact vector form with which it is possible (based on data from numerical simulations) to determine the vector ⃗ components and its magnitude (distance cp* from cm). The term (11) is simple and does not require solving a new equation system, as would be the case if the optimization was done by the Lagrange multipliers method. The actual center of pressure is in the same direction as cp*, but it cannot be obtained directly from the moment equations, but an additional condition should be set.
It is important to note that, if the orthogonality condition is precisely satisfied, the intensity magnitude of the vector ⃗ could be obtained by simply dividing the magnitude of the moment vector ⃗⃗⃗ by the force vector ⃗ magnitude, but in this case the position of the vector (components dx, dy and dz) in the space is unknown. The vector ⃗ is an important parameter because it represents the measure of the aerodynamic moment and shows the position of the total aerodynamic force.

Condition of stability for irregularly shaped bodies
The question is how the body will orient itself during the flight after the kinetic energy of the rotation is completely dissipated through the work of aerodynamic moment (Kljuno and Catovic, 2017). The underlying precondition for such a body orientation is that there are cases when: Condition (12) is necessary, but not sufficient. To understand this better, let us use an analogy with a pendulum, as in Fig. 5. In the case of a) and b) vector ⃗ magnitude is zero (d = 0), but in the first case, the pendulum is in an unstable orientation, whereas in the second case ( Fig. 5b) it is in a stable orientation. Therefore, it is necessary to introduce an additional (sufficient) condition of stability of orientation.
The basic idea is to use work of moment ⃗⃗⃗ to define this additional stability condition. Before we begin the analysis of the irregularly shaped body, let us consider again the case of the stability of the pendulum, now in Fig. 6. If we perform pivoting (rotation) of the pendulum from the upper vertical position for the small angle φ in the direction indicated in Fig. 6a, then the work of the moment due to the force ⃗ is positive (W> 0). Conversely, the work of the moment due to the force ⃗ over a rotation for a small angle φ around the lower (stable) vertical position is negative (W <0).
Based on this simplified analysis, the conditions of body stability can be set as follows: in every neighborhood of equilibrium position.
The condition (14), in the case of an irregularly shaped body flight, refers to the arbitrary rotation axis which passes through the center of mass.
In the case of a pendulum, the rotation was around an axis perpendicular to the direction of the force. However, in the case of an irregularly shaped body, there are an infinite number of rotation axes or the ways in which the body can be pushed away from equilibrium orientation, so it is necessary to present the condition (14) in a more appropriate form.  For easier understanding, we can again use an analogy.
Let's imagine that the pendulum joint is spherical as shown in Fig. 7. Then it is possible to pivot the pendulum away from its the equilibrium orientation by rotating it around any axis perpendicular to the direction of the force ⃗.
In Fig. 7a the pendulum is shown in the equilibrium orientation, while in Fig. 7b the pendulum is rotated for angle φ in relation to the equilibrium position around an arbitrary axis, whereby ⃗⃗ ⊥ ⃗ (the force ⃗ is moved to the joint and the moment is added). The work of a moment ⃗⃗⃗ ⃗⃗ for a small angle of rotation (for pendulum) can be determined as: where the work done by the gravity force, expressed via the work of the moment about the joint, has to be negative in the neighborhood of the angle = 0, for which the moment is zero, i.e.
( ) < 0, for a small angle in the neighborhood of = 0.
The condition (16) can be generalized in the case of an irregularly shaped body, as in Figs. 8 and 9. Fig.  8 is an illustration of a body stable orientation, wherein the vectors ⃗ and ⃗ are in general not parallel, while in Fig. 9 the body was taken away from stable orientation by rotating it with angle about arbitrary axis, where ⃗⃗ ⊥ ⃗ . The work of a moment ⃗⃗⃗ for rotation with an angle φ is Unit vector ⃗⃗ can be expressed as ⃗⃗ = cos ⃗ + cos ⃗ + cos ⃗⃗ .
The work of a moment ⃗⃗⃗ (17) for angle φ can be approximately determined with the equation: Axis Ω -Ω is set perpendicular to ⃗ as there is no work done over a rotation about an axis in the direction parallel to ⃗⃗⃗⃗ . If we go back to the pendulum analogy, this axis would be the analogous to axis of the pendulum, in the direction of force ⃗, and then the moment for point A remains zero (Fig.  7), so there is no work on that rotation.
The condition of orthogonality for the rotation axis Ω -Ω and ⃗ is Only one of three angles of rotation axis Ω -Ω is arbitrary, while the other two are determined from (19) and (21). Based on (16) and (17) the condition can be set for the general case The work function (23) is zero at = 0, hence the condition (24a) and the (de)stabilization moment ⃗⃗⃗ is zero at the equilibrium orientation due to ⃗ = 0. So, given that ⃗⃗⃗ for = 0, the sufficient condition of stability is finally

Method of numerical simulation
Numerical simulations of high-speed external aerodynamic flow over irregularly shaped body were performed in Ansys Fluent. Body with irregular shape was scanned and digitized in AutoCAD software, based on its real geometry and imported in Ansys Workbench. Unstructured mesh around the body (700 000 cells) and boundary conditions (wall and pressure farfield) were used in the numerical model ( Fig. 10 and Fig. 11). 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. 11).
Air is modeled as homogeneous, isotropic, ideal gas with pressure-temperature dependent density , specific heat Cp, thermal conductivity k and dynamic viscosity . Simulations were carried out for a flow velocity of 3 Mach (assumed initial velocity of fragments from detonating high explosive projectiles), and body orientations for one full rotation (from 0 to 360), with angle increment of 15 (Fig. 12).
These orientations were used in simulation with the static body to estimate aerodynamic forces and moments that are generated during the flight. 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 particularly effective for transonic flows around the aero profiles, including flows with significant separation of the boundary layer (Pope, 2000). Residual tolerance was set to 10 -5 . Aerodynamic forces and moments, required for determination of the center of pressure and body stability, were determined using user-defined function (UDF) program, written in C programming language. Based on total pressure on the body (tangential stresses were ignored due to the dominant pressure forces) and area of every cell, forces were calculated for every cell of the body and then summed up to obtain total aerodynamic forces and their components. The aerodynamic moment is determined for each cell of the model as the cross product of the cells vector radiuses and the aerodynamic force acting on that cell, and these moments were also summed to obtain the total aerodynamic moment, as well as its components by axes. 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 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. 13 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 a cube

Analysis of results
In this paper center of pressure was determined for different orientation of an irregularly shaped body. Also, the stability for irregularly shaped body was assessed using numerical simulation results and physical model given.

Determination of center of pressure
Using a physical model developed in the research, using the eq. (11) and (12) for different body orientations (Fig. 12) of an irregularly shaped body, vector ⃗ component values (dx, dy i dz) and its magnitude (practically distance of the center of mass from the center of pressure) were determined. In Table 1 the values of moment components, obtained by numerical simulations are also given.  Table 1 shows that the magnitude of the vector ⃗ is not equal to zero for any orientation considered in simulations. Also, based on the vector ⃗⃗⃗⃗ component values, it is shown that the vector ⃗ is positioned arbitrarily in space, i.e., it is not bound to a single axis, such as for axi-symmetric bodies.
Generally speaking, the largest magnitude of the vector ⃗ were recorded in the case of a body orientation of 75° -105° and the minimum for orientations in the range of 135° to 225°.

Assessment of stability for irregularly shaped body
The physical model, developed in the research, has shown that stability conditions can be given in the form | ⃗ | = 0 and ( ) < 0.
It was necessary to find the orientation of the irregularly shaped body whereby these conditions were most satisfying for given velocity of the body (or relative air flow). The body used in the simulations represents a realistic fragment generated by the detonation of high explosive projectiles (initial velocity of these fragments is in the order of 3-6 Mach). Therefore, the air flow velocity of 3 Ma was of the greatest interest in numerical simulations.
Based on the results (shown in Table 1), the estimation of the stability zone for two orientations of the irregularly shaped body was conducted, where the sign of moment Maero-z (rotation is conducted around axis Z, Fig. 12) changes, and also where magnitude of vector ⃗ was lowest and moment intensity relatively low. For this body these orientations were: 150 i 195.
It was necessary to find the values and the sign of the work of a moment W (φ) around these equilibrium orientations ( ⃗⃗⃗ ≈ 0), using the eq. (20), for the angle of increment of 15, where for one side (relative to the equilibrium orientation) this angle increment is -15° (i.e. angle 135° relative to the equilibrium orientation of 150°), and on the other side this increment is 15 ° (i.e. angle 165° relative to the equilibrium orientation of 150°). Numerical simulations provided the values of the aerodynamic moment components (for given Mach number), which are needed to determine the work of a moment from the Eq. 20.
In the simulations, air flow velocity vector was oriented in positive X axis (Fig. 12), meaning that the angles of the rotation axis from the expression (20) were α = 90, β = 90 and γ = 0. In this way, rotation axis coincides with the rotation of the body in numerical simulations (rotation around Z axis, Fig.  12). Table 2 shows the results of the work of aerodynamic moment W (φ) in the equilibrium orientations of 150° and 195°, at the velocity of 3 Mach, obtained from the expression (20) and the numerical simulation data. As can be seen from Table 2, for the first equilibrium orientation (150) the body is in the unstable zone because the work of moment around these orientations (for angle increments of 15) is positive. On the other hand, in the orientation of 195, the body is in a stable zone because the work on this incremental rotation (15) with respect to the equilibrium orientation (195 0 ) is negative (inequality (25)).
The fact that the work of moment for the 195° orientation is negative indicates that the oscillation of the body around the center of its mass will occur with respect to the selected axis of the body rotation, but that does not mean that the body would not rotates around some another axis. However, the oscillation will occur within the stability zone in the neighborhood of the stable equilibrium orientation which the body takes during the motion through the resistive medium.
This stability zone refers to only one axis, the one around which the body was rotated in numerical simulations. In order to test this zone of stability more closely, it would be necessary to make an estimate of the work of moments for an axis perpendicular to the first axis (used in simulations). For each step around Z axis one should do 24 more steps around Y axis rigidly fixed to the body, which would require 24 2 (or 576) numerical simulations. But even then, this rotation does not have to give an exact solution of the stability zone. Instead, 24 numerical simulations based on rotation only around Z axis were made, and based on this; an estimate of where a stable orientation of the body is located was assessed. But, in general, this stability zone for body with an irregular shape should be iteratively examined in relation to rotations of mutually perpendicular axes, which is a long-term work.
It should be noted that these stability estimates apply only to the body of this particular shape. For the other body, it is necessary to make new simulations, but the developed method for prediction of center of pressure and stability is, of course, general, and based on it as a tool, it is possible to evaluate the stability of any body.
Finally, the stability zone is actually a spatial angle. This can be a conical surface, and it is necessary to estimate the boundary of this conical surface at rotation about one axis and then do the same for the second axis. Within this conical surface, the body will oscillate around the center of mas and will not overturn or tumble. If there is no initial angular velocity, the body will be in a stable orientation during the motion and will oscillate within this cone surface.

Conclusion
Physical models for estimating the center of pressure for irregularly shaped bodies and for assessing the stability of such bodies are given. Using the developed physical models and the results of numerical simulations (3D steady simulations, density based solver and Spalart-Allmaras turbulence model) for an irregularly shaped body, a center of pressure and stability analysis was performed for one specific body shape.
A method of determining equilibrium stable orientation based on the work of aerodynamic moment in the neighborhood of the equilibrium orientation is demonstrated in the general case of irregularly shaped body movement. An example of a practical estimation of stable orientation is given using simulations for a series of orientations, with the body rotating in increments of 15 around one of the central axis.
It has been shown that the body in specific orientations may be in the zone of stability, but this stability should be further examined iteratively with respect to rotations of mutually perpendicular axes.

Future work
Planed future work will be focused onto an estimation of the fragment trajectory using aerodynamic coefficients for real fragments. It is a challenging task since the general fragment motion has 6 DOF and relatively long trajectories due to high velocities. Particularly challenging task is solving for the body rotation due to the aerodynamic moment which is difficult to predict for an arbitrary orientation of the fragment during the flight. However, the plan is to introduce appropriate approximations and model the rotational degrees of freedom.