2. Flight in a Vertical Plane
The tumbling motion of the body is considered in a vertical plane, and three sets of orthogonal axes are shown in
Figure 1: (i) Earth fixed axis, with altitude
Z and opposite gravity
g, and coordinate
X along flat ground; (ii) wind axis, with
along airspeed
U and opposite to drag, and
orthogonal along the lift
L; (iii) body axis with
along the aircraft datum and
orthogonal. The angle between
and
is the angle of attack
, and the angle between
and
X is the flight path angle
. Equations of motion in a vertical plane are written in wind axis, so that the wind velocity is included in the airspeed
V, and they state the balance of forces tangent—Equation (1a)—and normal—Equation (1b)—to the flight path and the pitching moment—Equation (1c):
where
U is the modulus of the total velocity,
is the total flight path angle,
L is the lift,
D is the drag,
m is the mass,
I is the pitch moment of inertia,
is the angle of attack,
is the pitch angle,
M is the pitching moment and
g is the acceleration of gravity. Since the motion of the tumbling body is unpowered, there is no thrust in the previous equations.
The Cartesian coordinates of the trajectory,
and
, in a vertical plane are decomposed into the sum of the following: (i) a parabolic free fall under uniform gravity; (ii) a perturbation,
and
, due to the aerodynamic forces:
In Equations (2a) and (2b), the pair corresponds to the initial position of the body and the pair is associated to the initial velocity components due to the gravitational force; however, it will be assumed that the last pair is and consequently the components of the total initial velocity are given by the pair .
The horizontal and vertical velocities are obtained by differentiating, respectively, Equations (2a) and (2b):
that imply Equation (4b) for the modulus of the total velocity (4a):
and is related to the modulus of the relative velocity
V in free fall coordinates by
The components of the total velocity in Earth fixed coordinates
agree with the total flight path angle:
and the components of the relative velocity (3a) and (3b) in free fall coordinates
lead to the relative flight path angle:
Substituting (5a) and (5b) in (4d) leads to the following:
whose roots specify the total velocity in Earth fixed reference frame
U in terms of the relative velocity
V in free fall coordinates
Here, the flight path angle
and initial velocities
, along with acceleration due to gravity
g and time
t, explicitly and implicitly appear in
.
The three coupled equations of motion—Equations (1a) to (1c)—form a fourth-order system. Since free fall coordinates (2a) and (2b) already include gravity, the approximation is made that, in the free fall reference frame, gravity can be omitted from Equations of motion. This free fall coordinate approximation will be checked in subsequent work, considering the same problem without this approximation. As a simpler preliminary approach to the problem, the free fall coordinate approximation is made, leading to a third-order system, as shown below. Thus, the free fall coordinate approximation includes the effect of gravity in the moving non-Galilean reference frame (2a) and (2b), but omits it in Equations of motion:
where
,
and
are respectively the lift, drag and pitching moment coefficients that are assumed to depend only on the angle of attack,
l is a reference length and
V is the velocity in free fall coordinates. The free fall coordinate approximation retains the angles of attack
and flight path angle
, with the supposition that the absolute and relative flight path angles are of the same order
, and assumes that the acceleration of the latter is much smaller than that of the angle of attack
, so that Equations (8a) to (8c) after simplification of (8c) become
Consequently the relative velocity is expressed in terms of the angle of attack by solving Equation (9c) as follows:
Substitution of Equation (
10) into Equation (
9a) leads to a differential equation for the angle of attack alone:
that is of the third-order and can be rewritten with an explicit third-order derivative:
In the last equation appear the following: (i) the radius of gyration
r that is related to the moment of inertia
I and mass
m by the relation:
(ii) the dimensionless mass factor defined by the fraction
as the ratio of the mass of a cylinder of air with density
, cross-section
S and length
l to the mass
m of the tumbling body. The mass factor is less than unity if the tumbling body is heavier than its volume of air that it displaces; this should be nearly always the case by some margin. Using Equations (
12a) and (
12b) in Equation (
11b), the third order differential equation becomes
This equation is nonlinear and remains valid for any form of drag coefficient
and pitching moment coefficient
, both of which depend solely on the angle of attack. The solution of Equation (
13) specifies the angle of attack as a function of time
and substitution in Equation (
10) specifies the modulus of the relative velocity
. The substitution of Equation (
10) in Equation (9b) leads to the following:
This equation specifies by integration the relative flight path angle
. The initial conditions for the third-order differential equation are the angle of attack and its first- and second-order time derivatives; the formula for the modulus of the relative velocity
V (
10), arranged in the way:
is used to evaluate the second-order time derivative of the angle of attack
, particularly at initial time.
So far, no assumptions have been made on the dependence of the drag, lift and pitching moment coefficients on the angle of attack.
3. Aerodynamic Coefficients in a Full Circle
In order to consider tumbling motion, the aerodynamic coefficients must be specified for the full range of angles of attack. Although the behaviours of lift, drag and pitching moment coefficients have been thoroughly studied for small angles of attack for aerofoils or streamlined shapes, this is not the case for the tumbling motion of irregular shaped bodies, considering the full revolution of angles of attack. In order to extend the “usual” aerodynamic theory to tumbling motion, some simple plausible assumptions are made.
For the lift coefficient, the starting point is the expression for the lift coefficient for aerofoils at small angles of attack [
5,
6], where
corresponds to the slope for small angles of attack and the angle is measured relative to zero lift:
For a tumbling body, the lift coefficient should obey three conditions: (i) for small angles of attack, the resulting expression should be the same (16) as the one obtained using small angle theory; (ii) when the body is perpendicular to the flow (
) there should be no lift; (iii) the lift coefficient variation must be periodic with period
. Equation (16) does not meet the third requirement. Multiplying Equation (16) by
leads to the following:
This equation meets all three requirements: (i) for small angle of attack, then and the expression (17) reduces to (16); (ii) if the body is perpendicular to the flow, then and the lift coefficient (17) is zero; (iii) the lift coefficient (17) is a periodic function with period , since .
The next aerodynamic coefficient to be considered is the pitching moment coefficient. Since the pitching moment is taken relative to the centre of mass of the rocket whereas the lift is applied in the centre of pressure, the pitching moment can be assumed to be proportional to the lift:
through the distance
a between the centre of pressure and the centre of mass of the rocket. The dependence of the pitching moment on the angle of attack is thus similar to Equation (17), with a distinct pitching moment slope in the next equation:
in which the derivative of the coefficient moment is proportional to the derivative of the lift coefficient with respect to the angle of attack through the following relation:
that follows from
The relation (17) also implies
The only aerodynamic coefficient still to be considered is the drag coefficient. The drag coefficient will be given by the sum of three terms, of which the first two are the usual in aerodynamic theory [
7]:
namely: (i) the first term, specifically
, is independent from the angle of attack and represents the friction drag; (ii) the second term is proportionate to the square of the lift coefficient, which ensures a parabolic lift–drag relationship and is commonly referred to as induced drag. To accommodate the wide angle of attack experienced by a tumbling body, a third term is introduced, representing the "wake drag" of a lamina at an angle
relative to a stream that features a "vacuous" region behind it. This region lies between the free streamlines originating from the lamina’s tips [
5,
6].
The expressions for the pitching moment coefficient (19) and its time derivative (22) appear in the non-linear third-order differential Equation (11b) for the angle of attack, whose square takes the following form:
whose right-hand side involves: (i) the drag and pitching moment coefficients that depend on angle of attack; (ii) two dimensionless parameters that do not depend on any of aerodynamic coefficients. The first dimensionless parameter
is the moment of inertia of the rocket divided by the product of the mass by the square of the lengthscale:
The lengthscale can be taken to be the rocket length, which is usually greater than the radius of gyration relative to the pitch axis, so that
is less than unity. The second dimensionless parameter
equals the first parameter (
25) multiplied by the mass ratio (
12b), which is also smaller than unity:
so that their product is also less than unity. Using the parameter
obtained in Equation (
26), the non-linear third order differential equation for the angle of attack (
24) becomes
in which all terms have derivatives with regard to time whose orders always add to six.
After determining the dependence of the angle of attack on time
as the solution of Equation (
27), the relative flight path angle
is obtained by integrating the Equation (
14). In this integration, the dimensionless parameter
defined in Equation (
26) can be substituted, resulting in the following:
The relative velocity (
4c) in free fall coordinates and flight path angle (
28) determine the total velocity (
7b) in Earth-fixed coordinates and integration of (5a) and (5b) specifies the trajectory. The third-order non-linear differential equation (
27) for the angle of attack has a unique solution satisfying three independent initial conditions at time zero. These are as follows: (i) the initial angle of attack:
(ii) its first-order time derivative at time zero:
(iii) its second-order time derivative at time zero:
which is specified by Equation (
15), knowing the initial velocity
at time
, repeated here for convenience.
The analytical solution is obtained for small angles of attack (
Section 4) followed by numerical solutions for unrestricted angles of attack (
Section 5).
4. Particular Analytic Solutions
The main innovation in the present paper is to consider the full circle of angles of attack for the flight trajectory of a tumbling body. It is useful to start by checking compliance with the usual aerodynamic theory for small angles of attack , in which case some simple analytic solutions are possible. This also applies to the initial motion in the particular case of release at small angles of attack. The analytical results provide initial insights into the problem, serving as a preliminary step before delving into the general case of release and motion at large angles of attack. The subsequent numerical treatment of this scenario constitutes the main focus of the paper.
By restricting the angle of attack domain to only small angles, all trigonometric functions disappear from the differential equations, and both the lift coefficient and pitching moment coefficient will vary linearly with the angle of attack relative to the angle of zero lift. The drag coefficient, which was originally a sum of three terms, will become a constant, due to the friction drag, and a linear term, due to the wake drag. The induced drag is neglected because it is a non-linear term and hence too small to be taken into account again in the particular case of small angles of attack. The lift (17), pitching moment (19), and drag (
23) coefficients simplify for small angles of attack,
, respectively to the following relations:
If
holds using the previous relations, the differential Equation (
27) for the angle of attack becomes
In all terms of the non-linear third-order differential Equation (
31), the sum of the orders of the time derivatives is 6 and the powers of the angle of attack is 4. Thus, a possible solution is a power of time with amplitude
A and exponent
:
This is considered as the simplest possible analytic solution in the particular case of small angle of attack. It will be compared subsequently with more general numerical solutions for unrestricted angles of attack. Substitution of the last Equation (
32) in (
31) leads to the following relation:
where
is a dimensionless parameter and is related to Equation (
26) by
For a stable tumbling body
; thus, the parameter from Equation (
34) is positive,
. This would apply if the whole rocket is the tumbling body, which is the usual worst case assumption. However, if the rocket breaks up into smaller pieces, each of these would have to be considered separately with the initial conditions at break-up time.
From Equation (
33), it follows that the amplitude
A is arbitrary whereas the exponent
can take four values. Since the differential Equation (
31) for the angle of attack is of the third order, the general solution involves three arbitrary constants. In Equation (
32), only one arbitrary constant appears,
A; thus, it is a particular solution. In total, Equation (
33) has four solutions. The first two solutions to Equation (
33) are
and
. The first solution,
, applies to a body falling without any kind of rotation, because the angle of attack is constant and the time derivatives are all zero for all time. The second solution,
, applies to a body falling at a constant rate of rotation, since the first derivative of the angle of attack is constant and both the second and third time derivatives are zero. In both cases,
or
, the angular acceleration is zero, the relative velocity is zero by Equation (
10) and the flight path angle is constant by Equation (
14), so the free fall velocity components
are nearly zero, and the total velocity components
are linear functions of time in (3a) and (3b), leading to a parabolic trajectory. The deviations from a ballistic parabola occur for the remaining two values of
, which are the roots of the term in square brackets in Equation (
33), leading to the quadratic equation
The roots of the last equation are
and depend only on the dimensionless parameter defined in Equation (
34), leading to four cases:
Figure 2 shows the analytic solutions for the four previous cases of the roots
and
of the term in square brackets in Equation (
33). The four cases in
Figure 2 correspond to the next four values of
:
,
,
and 100. In the first case, when
, the two roots
and
tend to be equal to the first two roots calculated from (
33), specifically
and
. In the third case, when
, the two roots
and
are complex numbers; therefore, each one leads to two distinct solutions of
, evaluated from (
32), one associated with the real part and the other with the imaginary part of
. Moreover, in the third case, when choosing the complex root
or
, the real part of
is the same, whereas the imaginary part of
is symmetric. Both solutions are plotted in
Figure 2. The existence of four possible cases is due to the competition of two different effects. The first effect is the angle of attack oscillation caused by the pitching moment (9c); in aeronautics, this effect is known as the short period mode. The second competing effect is due to the balance of forces, given by Equations (9a) and (9b), and corresponds to the phugoid mode in aeronautics. The relative importance of the two effects is specified only by the dimensionless parameter defined in Equation (
34). An example for the four sub-cases is provided in
Figure 2. In most cases, the positive root solution
from Equation (
36) leads to a faster mode (the angle of attack increases faster) than with the negative root
.
The power type solutions for the angle of attack (
32) with exponents (
36) are given by
lead for the relative velocity (
10) to:
where Equations (30b), (
38), (
35), (
34), (
26) and (
12b) were used successively in the last relations. The rate of change with time of the relative flight path angle (
28) is given by:
where Equations (30a), (30b), (
38), (
35) and (
34) are successively applied in the last relations. The integration of the last equation leads to the flight path angle as a function of time:
where
is the initial value.
The four cases (Equations (37a) to (37c)) are considered next in turn. The case from Equation (
37a) leads to
and
, which are the parabolic ballistic trajectories considered previously. In all remaining cases, the velocity (
39) diverges for zero time:
and decays to zero for a long time:
The angle of attack (
38) and variation of the flight path angle (
41) scale similarly with time: (i) in the case (37b) leading to
(ii) in the case (37d) leading to
In both cases, the angle of attack (and also the flight path angle) vary monotonically with time; that is, the body rotates as it falls along a modified parabolic trajectory.
In the case from Equation (37c), the roots
lead to oscillations of the angle of attack:
whose amplitude is zero initially:
and diverges for a long time:
The oscillation frequency (45b) vanishes in the transition
to the monotonic case (37d). The real solutions corresponding to (
46) are the real and imaginary parts:
In this case, the motion involves an oscillation in the angle of attack as the body falls.
Thus, the analytic solutions, although valid only for small angles of attack, show the existence of two regimes of fall depending on the parameter (
34) and (
26):
If the parameter (
49) is large (37d), the tumbling body rotates as it falls (
Figure 3a); this corresponds to a predominance of the phugoid mode. If the parameter (
49) is small (37c), the tumbling body oscillates as it falls (
Figure 3b); this corresponds to the predominance of the short-period mode. A simple analogy [
8] is the motion of a circular pendulum of mass
m, length
s in the gravity field
g, for which the difference in potential energy between the lowest
and highest
points is:
The total energy of the pendulum, comprising both kinetic and potential components, is given by
and remains conserved. This leads to three cases: (i) if the total energy (
50b) is less than the difference in potential energy between top and bottom (
50a),
, the pendulum cannot reach the top, the pendulum motion is oscillatory (
Figure 4b), similar to the short-period mode of aircraft, and the oscillating fall of tumbling body (
Figure 3b); (ii) if the total energy (
50b) is larger than the difference in potential energy between top and bottom (
50a),
, the pendulum goes through the top with non-zero velocity, the pendulum motion is circulatory (
Figure 4a), corresponds to the phugoid mode of aircraft and to the rotating fall of a tumbling body (
Figure 3a); (iii) in the intermediate border line case of total energy (
50b) equal to the difference (
50a) of potential energy at the top and bottom,
, the pendulum would stop at the top, corresponding to
in (
36), leading to
, while the angle of attack (
32) and flight path angle (
41) both increase with time as local velocity (
39) decreases.
Since the differential Equation (
31) for small angles of attack is non-linear, the principle of superposition does not hold and a linear combination of the solutions from Equations (
38) or (
46) is generally not a solution. Since the analytic solutions are restricted to small angles of attack, the limits of Equations (
43a), (44b) and (
47b) hold only for small values. However, the existence of monotonic rotation and oscillatory solutions is confirmed by the numerical solutions (
Section 5) that hold for unrestricted angles of attack.
7. Effect of the Rocket Parameters
Proceeding to the sensitivity analysis of the effects of the rocket parameters, the rocket mass was found to have a direct influence in both trajectory and fall mode. In terms of fall mode, smaller masses lead to smaller periods and amplitudes of oscillation, as is shown in
Table 6, which was obtained using the second configuration. For an extremely heavy body, the period of oscillations is so large that the fall mode becomes RM. Using the first configuration leads to the RM conclusions, independently of the mass of the tumbling body.
In terms of trajectory (
Figure 10), larger masses have the longest range and the smaller masses the shortest. Larger masses are less affected by aerodynamic forces and come closer to the ballistic trajectory while lighter masses are more affected by drag and travel less far. The usual assumption that the heaviest tumbling body is usually the worst case (that travels farther) is supported by the results obtained in this analysis. The effect of mass is less significant for heavier bodies, by observing that changing the value of mass from 500,000 kg to 1,500,000 kg does not change the range as much as changing the mass from 5000 kg to 500,000 kg.
Smaller masses also have a second effect other than reducing the trajectory range, namely introducing some wobbles to the usually ballistic shaped trajectory. These wobbles are a consequence of the following: (i) the flight path angle amplitude of oscillation increasing as the mass of the rocket is reduced, which supports the explanation that the flight path angle was seemingly constant because of the large mass value; (ii) large velocity oscillations as the mass is reduced.
The radius of gyration of the rocket also has an effect on both the trajectory and fall mode. Smaller values of the radius of gyration usually lead to larger ranges not only when the RM is dominant (
Figure 11), but also when the OM configuration is used (
Figure 12). An increasing radius of gyration or moment of inertia implies greater kinetic energy of oscillation in both modes and reduces translational energy and range.
In terms of fall mode, as can be seen in
Table 7, larger radius of gyration values work in favour of the OM fall by reducing the amplitude of the oscillations.
The height of the rocket was found to have effects on the fall mode and trajectory, as can be seen in
Table 8 and in
Figure 13 and
Figure 14. In terms of the fall mode, having a larger height will favour the RM. Increasing the height of the rocket maintaining the mass constant is the same as increasing only the moment of inertia of the rocket; therefore, a larger moment of inertia leads to a stronger tendency to rotate.
In terms of trajectory, the variation to the height of the rocket will lead to the same range changes as the increase of the gyration radius since both parameters contribute to calculating the inertia of the rocket.
The distance between the centre of pressure and centre of mass is very important toward determining the fall mode. The results of
Table 9 are inconclusive since it is difficult to predict what the fall mode is for a certain distance
a.
Table 9 also shows that the fall mode is dependent on the initial angle of attack, thus predicting that there are more parameters, not depicted on
Table 9, that significantly influence the fall mode.
Assuming that the initial angle of attack is close to
, when the distance
a is negative, the fall mode is RM, but when the sign of that distance changes, the movement becomes OM, generally with the period of the oscillation decreasing when the distance
a increases. When the initial angle of attack is close to
, the sign of the distance
a does not predict the fall mode; nevertheless,
Table 9 shows that in almost all cases the fall mode is OM.
Lastly, regarding the aerodynamic coefficients, while
has its most noticeable effect in the fall mode,
and
K can exert their influence on the trajectory and fall mode. Larger values of
decrease both the period and amplitude of oscillation as can be seen in
Table 10.
The results for
and
K are the same since an increase in either parameter increases the value of the drag coefficient. Larger drag values lead to trajectories with shorter ranges because if the drag coefficient is higher, then the resistance to the motion is also higher. This is illustrated in
Figure 15.
The increase in drag can also change the fall mode by itself. Larger values of friction drag contribute to an increase of the period of oscillations as it is shown in
Table 11. Using the first configuration leads to the RM conclusions, independently of the mass of drag coefficient.
Again, since a larger drag coefficient means more resistance to the movement of the rocket, it is only natural that, as the drag increases, it takes more time for the rocket to make one complete oscillation.