Geometry of motion and nutation stability of free axisymmetric variable mass systems

In classical mechanics, the ‘geometry of motion’ refers to a development to visualize the motion of freely spinning bodies. In this paper, such an approach of studying the rotational motion of axisymmetric variable mass systems is developed. An analytic solution to the second Euler angle characterizing nutation naturally falls out of this method, without explicitly solving the nonlinear differential equations of motion. This is used to examine the coning motion of a free axisymmetric cylinder subject to three idealized models of mass loss and new insight into their rotational stability is presented. It is seen that the angular speeds for some configurations of these cylinders grow without bounds. In spite of this phenomenon, all configurations explored here are seen to exhibit nutational stability, a desirable property in solid rocket motors.

Position vector from S * to P P A generic particle within C r Position vector from O to P r * Position vector from O to S * R Radius of exit plane, and initial radius of cylindrical body S * Mass center of system V 0 Volume of C v Inertial velocity of P v r Velocity of P relative to B z Instantaneous half length of cylinder z e Shortest distance between S * and exit

Introduction
Recent examinations into the dynamics of variable mass systems can be divided into two distinct periods; the first begins in the mid-1940s and spans the early 1970s, while the current generation of studies begins in the 1990s. The main interest during the first stage appears to be interest in the rocket modeling problem; here, the contributions of Rosser et al. [1] and Thomson [2][3][4] are essential reading on the topic. Whereas Rosser et al. were the first to present the concept of jet damping of rotating rockets, Thomson added to the body of work on more general variable mass systems [2] and also examined jet damping in the specific case of solid rocket motors [3]. His classic text [4] is notable for presenting one of the first broadly accessible discussions on attitude dynamics of rocket-type systems.
The second generation of studies was more expansive. Apart from revisiting the derivation of the equations of motion of continuous mass-varying systems [5], attention has also been paid to algorithm development for simulating rigid [6] and flexible systems with mass variation [7,8], as well as analyses of systems with discontinuous mass variation [9,10]. Yet another aspect that has been scrutinized in this period is the stability of the rotational (or attitude) motion of spinning rigid systems with continuous mass loss. These investigations into the stability of axisymmetric rigid systems with mass loss [11][12][13] were motivated by the anomalous coning motion seen in some solid rocket motors (SRM) [14]. The resulting growth in the cone angle or nutation angle leads to a velocity pointing error and is referred to as a coning or nutation instability. The set of papers by Eke and Wang [11,12] began a systematic re-examination of attitude dynamics of axisymmetric variable mass systems. In [11], axisymmetric cylindrical systems with mass loss are examined; the work presented in the present manuscript most closely relates to their work. In [12], a more general version of the axisymmetric system with mass loss is examined where they show that accounting for the whirl component of internal flow in these systems does not affect the magnitude of the transverse angular velocity. This finding informs the approach of the current paper where a similar non-whirling flow pattern is assumed. Following up, Eke and Mao [13] studied a composite variable mass system: a constant mass spacecraft with mass-varying cylindrical propellant. More recently, Ha and Janssens [15] investigated the CONTOUR mishap using an identical model for mass variation as the aforementioned papers. As the CONTOUR spacecraft exhibits slight variations in its transverse moments of inertia, an approach to analytically determine the angular speeds of such a nearlysymmetric system with mass loss has recently been developed and numerically validated [16].
A limitation of some of these studies on variable mass cylinders [11,13,16] is that stability information is heuristically interpreted by examining the evolution of the angular velocity. This paper develops an alternative graphical approach to assess the coning motion by generating a solution to the second Euler angle. This is used to evaluate the system's nutational stability. The approach is also referred to as the 'geometry of motion' in classical textbooks on mechanics [17] and spacecraft dynamics [18]. Apart from allowing an analytic solution to the Euler angle characterizing nutation, it also permits visualizing the motion of a system's angular velocity vector from body-fixed and inertial frames; the document also covers this aspect. As shown in this paper, this technique offers a more mathematically rigorous explanation of the coning motion than one obtained from analyzing the angular speeds alone. It has previously been shown [11,13] that some idealized burns of axisymmetric cylinders can lead to unbounded growths in the transverse angular speeds. In this paper, it is shown that these burns still damp out the coning motion in spite of unbounded transverse speeds. A limitation of this work and prior studies that will need to be explored in a future study is the effect of thrust misalignment effects on the coning motion; such effects have been studied in greater detail recently for spinning rigid systems of constant mass [19][20][21].
The structure of this paper is as follows. First, the scalar equations of attitude motion for a general axisymmetric system are formulated and analytic solutions to these equations are obtained. Then, based on a recent result concerning the inertial fixedness of the central angular momentum vector of such systems [22,23], the approach to visualize the rotations of these systems from both body-fixed and inertial frames is developed. Finally, this theory is used to evaluate the motion of axisymmetric rigid cylinders subject to three idealized models of mass loss. From analysis, it is deduced that all of these freely spinning cylindrical configurations are nutationally stable, i.e., the coning motion damps out or remains bounded. Figure 1 is of a general variable mass system. It comprises a consumable rigid base B and a fluid phase F.

Equations of motion and the angular velocity of an axisymmetric systems
A massless shell C of constant volume V 0 and constant surface area S 0 is rigidly attached to B. It is assumed that mass can enter or exit C through the region represented as a dashed circle of radius R. At every instant, the shell and everything within it is considered to be the system of interest. The vector equation of attitude motion of such a system can be derived from the laws of mechanics (such as Newton-Euler equations, Lagrange's method, etc.) and appropriately invoking Reynolds Transport Theorem. This approach has been discussed in Ref. [5]. The following is the equation of attitude motion of a general variable mass system Several other forms of the equation of attitude motion have been derived for a variable mass system in the literature. However, Eq. (1) is most tractable for examining the attitude evolution as it is decoupled from the translational dynamics. The volume integrals in Eq. (1) are typically hard to evaluate in closed form since the internal flow field is not generally known within a system. However, in the case of rockets, reasonable assumptions can be made about this internal flow. It is typically assumed that the internal flow of gases within a rocket relative to its casing is both steady and axisymmetric. Further, assuming that the flow lacks a whirl component force the last three terms of Eq. (1) to disappear, making it less cumbersome. Thus, the equation The remainder of this paper concerns the torquefree motions of axisymmetric variable mass systems so M * = 0. Note that though this non-whirling flow assumption is not realistic, it has been shown that the magnitude of the transverse angular rates are unaffected by such a flow assumption [12]. As stability information in this paper is derived from the magnitude of the angular rates, the non-whirling flow assumption is sufficient.
Referring back to Fig. 1, b 1 , b 2 , b 3 are a dextral set of unit vectors affixed in B and parallel to its principal directions. The instantaneous central inertia dyadic of this axisymmetric system is where I and J are the instantaneous central principal moments of inertia of the variable mass system. Note that, the mass center, S * is not necessarily fixed relative to B and is always located on the symmetry axis, parallel to b 3 . Further, the inertial angular velocity of B at any instant is and, as a result, the inertial angular acceleration is The first three terms on the right-hand side of Eq. (2) evaluate to and respectively. The assumed steady symmetry of the internal fluid flow is further constrained to be uniform along the exit plane so that v r · n = u = constant, which then yields the following expression for the sur- whereṁ, the rate of mass loss, is given bẏ ρ exit is the assumed density of exhaust gases on the exit plane. Substituting Eqs. (6), (7), (8), and (9) into Eq. (2) yieldṡ anḋ Equations (11)-(13) are the scalar equations of attitude motion for an axisymmetric, torque-free variable mass system with axisymmetric non-whirling internal mass flow. Equations (11) and (12) are commonly referred to as the equations of transverse motion, while Eq. (13) is the spin equation of motion; ω 1 and ω 2 are transverse rates and ω 3 is the spin rate. The attitude motion of a variable system is characterized by solving these equations.

Spin rate
As the spin equation of motion is decoupled from the transverse equations of motion, an expression for the spin rate can be easily obtained and is given by where If J mk 2 3 , where k 3 is the time-varying axial radius of gyration, then Φ evaluates to Substituting for Φ in Eq. (14) gives the following expression for the spin rate Examining Eq. (17), the spin rate may decay, stay constant, grow, or fluctuate. Another observation that can be added to this is that the spin rate retains the polarity of its initial condition; if the initial condition is positive, then the spin rate is always non-negative, and vice versa. Further analysis in this chapter assumes a positive value for the initial spin rate. Additionally, if k 2 3 is assumed to be a constant in Eq. (17), then Equation (18) asserts that the spin rate can not fluctuate for a system with constant axial radius of gyration; it can only grow or decay monotonically, or stay constant.
It can then also be inferred that a time-varying axial radius of gyration can lead to fluctuations in the spin rate. Thus, the radius of gyration has a crucial effect on the spin rate. These comments on the spin rate are in agreement with the work of Snyder and Warner [24], and Wang and Eke [12].

Transverse rate
The transverse angular speeds can be solved for in a variety of different ways. Defining a complex variable ω c as allows Eqs. (11) and (12) to be combined to yield a differential equation linear in ω c : Integrating Eq. (20) yields where and A consequence of Eq. (23) iṡ whereχ has the units of angular speed. From Eqs. (21) and (19), the transverse speeds are given by 10 sin χ + ω 20 cos χ) (26) and are oscillatory in nature. Further, a transverse angular velocity vector ω t which lies in the plane made by the principal directions b 1 and b 2 can also be defined as whose magnitude is where ω 0 ω 2 10 + ω 2 20 = constant. The principal directions may be chosen such that ω 10 = 0 and ω 20 = ω 0 . This simplifies the appearance of the solutions in Eqs. (25) and (26) to Equation (27) can then be written as where b 12 = sin χ b 1 + cos χ b 2 ; b 12 is clearly a unit vector that is orthogonal to b 3 . These developments concerning the transverse angular velocity vector will be useful in developments in Sect. 3 on the geometry of motion. From examining Eqs. (31) and (22), it is evident that the spin rate has no influence on the magnitude of the transverse angular velocity. Further, it shows that the magnitude of the transverse angular velocity vector is a rescaling of ω 0 by Γ , a function of time that can either grow, decay, or fluctuate. b 12 is a unit vector that rotates in the plane of b 1 and b 2 at the angular speedχ. As revealed by Eq. (24),χ may be positive or negative (assuming that ω 3 is always positive), depending on the relative values of k 1 and k 3 ; if k 3 > k 1 , then b 12 rotates negatively in the body, advancing from b 1 to − b 2 and around; if k 1 > k 3 , then b 12 rotates positively in the body, coinciding with b 1 and then with b 2 and so on. As the interest of this paper is in understanding the nutational stability, no attempt is made to obtain an analytically expression for χ . Its value can be obtained numerically from Eq. (23).
An expression for Γ can be developed involving the system's transverse radius of gyration, k 1 . If I mk 2 1 , Γ may be reformulated in Eq. (22) as Here, R is the constant radius of the exit plane while k 1 and z e can vary; these are parameters that are modeling constraints. Substituting for Γ into Eq. (28) yields the following expression for the magnitude of the transverse angular velocity: Some of the effects of the variations in z e and k 1 on the transverse rate will be studied in the subsequent passage on variable mass cylinders. However, it should be evident that ω 12 can grow, decay, or fluctuate but under no circumstance is it constant for an axisymmetric variable mass system. Finally, the inertial angular velocity of an axisymmetric variable mass system can now be written as: the form of which mimics the constant mass rigid body case, for which Γ = 1.

Geometry of torque-free motion and nutation angle solution
Explicit description of the inertia parameters coupled with the developments regarding the angular speeds from the previous section is sufficient to visualize the motion of the angular velocity vector from the bodyfixed frame. However, it is also possible to visualize its motion in an inertially fixed reference frame which is the focus of this section. These developments are then used to study the variable mass cylinder problem in the following section. Exploiting the fact that the variable mass system's angular momentum vector is fixed in inertial space [22,23] enables us to visually study the angular velocity evolution from an inertial frame. The angular momentum of a general variable mass system can be written Retaining the assumption of steady and axisymmetric relative internal mass flow reduces Eq. (35) to With the definitions of I * and ω from Eqs. (3) and (4), respectively, the component form of the angular momentum vector is or Using the component form of ω t from Eq. (31) in Eq. (38) gives Equation (39) gives an alternate form of the angular momentum of an axisymmetric variable mass system. Axisymmetric rigid systems (of constant or variable mass) are unique in that their angular momentum and angular velocity vectors are always in the same plane, evident by comparing Eqs. (34) to (39); H * and ω are in the same plane formed by b 12 and b 3 and this is shown in Fig. 2. In this figure, θ is the angle between H * and b 3 , and β is the angle between b 3 and ω. These angles can be computed from the figure as and Evaluation of these angles θ and β along with the solutions to the angular velocities permit a visual understanding of the motion of the body in both inertial space and the body-fixed coordinate system. In the case of axisymmetric rigid bodies of constant mass, Γ = 1 and the spin rate and inertia scalars are constant. Consequently, in Eq. (40) and (41), both angles evaluate to constants and visualizing the motion of ω results in the formation of two cones: as it rotates about b 3 and another as it rotates about n h , a unit vector parallel to the inertially fixed H * . The two cones are appropriately named the body cone and space cone, respectively.
In the case of axisymmetric variable mass systems, angles θ and β are typically functions of time. The angles may grow or decay as dictated by the instantaneous values of the radii of gyration, Γ , and the spin rate. It is also possible, in the case of the same variable mass system, for β to exceed θ at certain instants, while it may not at certain other instants. Figure 2 represents an instant where β > θ, which occurs when k 3 > k 1 . Also, ω does not trace cones as it rotates about b 3 and n h . Thus, generic names are used for the surfaces it traces; in this document, they are referred to as the body surface and space surface (which reduce to the 'body cone' and 'space cone' for the limiting case of a constant mass rigid body, respectively).
In the spacecraft dynamics literature, θ is referred to as the nutation angle [18]. In classical mechanics, θ may also be identified as the second Euler angle  [17] and, thus, Eq. (40) also gives the solution to said Euler angle. Growth in nutation angles results in a coning motion which is undesirable in spacecraft. Potential hazards of excessive nutation or coning motion might be poor initial conditions for subsequent mission operations, or possible loss of communications with a spacecraft due to incorrect pointing between its antenna and the base station. One such case of excessive nutation growth was observed in missions powered by the STAR-48 SRM [14,[25][26][27][28][29]. To mitigate this growth, an active nutation damper is used in these missions [30].
Closed-form expressions to θ are available when the angular speeds are themselves available in closed form. These closed-form expressions are invaluable to not only evaluate stability but also for appropriate control strategy formulation to temper any instabilities. Solutions to the angular speeds are derived in the next section for a special class of variable mass systems: the axisymmetric variable mass cylinder. These closedform solutions are then utilized in constructing the body and space surfaces. The choice of cylindrical systems as an idealization is based on their use in several prior studies on rocket flight; Snyder and Warner [24] and Eke et al. [11,13] are some examples of investigators who have incorporated cylindrical propellant grains in their studies on the attitude motions of variable mass space vehicles. Figure 3 is that of a system that is initially a cylinder as represented by the dashed lines. The mass and inertia of this system is allowed to vary with time. As in the case of rockets, combustion is one plausible cause for these time-varying properties. Thus, at a general instant, some parts of the system may be solid while other parts are fluid. The solid part of the system is denoted B. The dotted line represents a control region C of constant volume and surface area that has the same shape as the original cylinder and is rigidly attached to B. Only matter within C at any instant is considered to make up the system. Interest here is in visualizing the motion of variable mass cylinders and understanding their nutational stability.

The axisymmetric, variable mass cylinder
Equations (11)-(13) govern the attitude motion of the variable mass cylinder, while Eqs. (17), (29), and (30) are the solutions to the angular speeds. Eke and Wang [11] identified idealized burn scenarios of which three are examined here. They are named the uniform burn, end burn, and radial burn. Each of these burns is explained in their respective sections. Inertia change due to fluid mass is neglected in all burn cases. The assumed properties of the system in creating the results presented are as follows. For the presented results, the density of the rigid propellant is assumed to be 1000 kg/m 3 . The initial angular speeds are assumed to be ω 0 = 0.2, and ω 30 = 0.3 rad/s; note that in a realistic rocket scenario, the spin rate would be an order of magnitude greater than the transverse rates. The duration of the burn is assumed to be 100 s unless mentioned otherwise. Other assumptions concerning the geometry are discussed under each burn scenario.

Uniform burn
A uniformly burning cylinder is one in which the mass of the system varies as a function of time, while the dimensions of the unburned portion of the cylinder are the same as that of the original cylinder; radius and length of the cylinder are both assumed to be 1m. Thus, its instantaneous central moment of inertia scalars are and This is an example of a variable mass system where the radii of gyration are constant where k 2 1 = R 2 /4 + h 2 /3 and k 2 3 = R 2 /4. The time derivatives of the moment of inertia arė The resulting spin and transverse rates are obtained from Eqs. (17) and (28) as where Γ (t) (m/m 0 ) 2h 2 /3k 2 1 . Prior to the start of the burn m = m 0 so Γ = 1. As the burn proceeds, Γ decreases with time but is always non-negative. Thus, the transverse angular rate also decreases with time for a uniformly burning cylinder. Then, Eq. (40) gives the nutation angle for this system where K = k 2 1 ω 0 /k 2 3 ω 30 . Since K is constant and Γ (t) is decreasing with time, the nutation angle is also a decreasing parameter with time. This reduction in the transverse oscillations of a variable mass system by the exhaust gas is referred to as jet damping of a rocket.
The rotation of ω about b 3 and n h was discussed in the earlier section on generic axisymmetric variable mass systems. Figure 4 visualizes the first of these rotations for the uniformly burning cylinder, assuming k 3 > k 1 . The angular velocity vector traces a spiral in the transverse plane as it rotates about the symmetry axis b 3 in a counterclockwise direction and eventually converges to the symmetry axis.
The rotation of ω in the inertial space (i.e., about n h since it is inertially fixed) requires knowledge of the angle, β, in Eq. (41) This angle also decays with time for the uniformly burning cylinder, effectively reducing the system's motion to that of simple or pure spin, i.e., a motion in which the angular velocity and angular momentum vectors are parallel. The above expressions of θ , and β along with the angular velocity permit the construction of the space surface as shown in Fig. 5. Note that n f , n g , and n h represent a dextral set of unit vectors which are inertially fixed.

End burn
In this mechanism of mass variation, the cylinder burns from the circular face closest to the exit plane toward its other face. At any instant, the unburned part retains the shape of a right circular cylinder. Figure 6 represents a cylinder that burns from right to left. The axial radius of gyration for such a cylinder is constant, As in the case of uniform burn, the spin rate can thus be determined from Eq. (18) to be constant, The more interesting development occurs in the transverse directions. The transverse radius of gyration is k 2 1 = R 2 4 + z 2 3 . Due to the changing length of the cylinder, k 1 is not a constant. Then, the transverse moment of inertia is If z is the instantaneous half length of the unburned cylinder, then its instantaneous mass is given by The rate of mass loss is theṅ which can also be written as From Eqs. (52) and (54), it is possible to recast Eq. (33) The above solution can also be expressed in terms of the transverse radius of gyration as Equation (62) explains that, for an end-burning cylinder, the transverse rate is regulated by three timevarying functions: the ratio of the initial transverse radius of gyration to its instantaneous value, the ratio of the instantaneous half length to the initial half length, and an exponential function. The first of these functions grows as the burn progresses but is always bounded as its denominator is never zero. The latter two functions, however, always decay with time. Thus, irrespective of the nature of the cylinder, the transverse rates for the end burn are always decaying; the rate of decay is seen to be faster for prolate cylinders when compared to oblate cylinders, which makes sense in view of Eq. (62). This has also been observed in Ref. [11]. The nutation angle for an end-burning cylinder is given by where K = ω 0 /k 2 3 ω 30 = constant. As both Γ (t) and k 1 are generally decreasing with time, the nutation angle decays with time for any end-burning cylinder. In other words, jet damping is observed for such a mass-varying cylinder. Figures 7 and 8 show the body and space surfaces for end-burning cylinders where the J 0 > I 0 (initial L = 1 m and R = 0.8 m) whereas Figs. 9 and 10 are for end-burning cylinders where J 0 < I 0 (initial L = 1 m and R = 0.5 m). It is interesting to note in the second case (where J 0 < I 0 ) that the angular velocity vector changes its direction of rotation, about both b 3 and n h . Initially, it rotates in a clockwise direction and then switches to a counterclockwise direction, eventually converging to a pure spin about b 3 in the body frame and n h in the inertial frame. This transition point occurs when the ratio of the axial radius of gyration to the transverse radius of gyration falls below 1 and changes the sign ofχ as given by Eq. (24). In either case, the end burn is seen to exhibit a decay in the cone angle and is thus nutationally stable.

Radial burn
In radial burn, combustion starts along the symmetry axis and proceeds radially outwards. Thus, at any instant, the unburned portion resembles a hollow pipe and both the axial and transverse radius of gyration of the system are variable as shown in Fig. 11. For the results presented here, initial radius and length are 1 m. The instantaneous moment of inertia scalars for the and Their corresponding time derivatives arė The angular speeds evaluate to Evidently, the spin rate is not constant for a cylinder subject to radial burn. It varies such that it appears to damp out, slowly, for about two-thirds of the burn. Toward the end of the burn, i.e., as r approaches R, the spin rate grows without bounds. Analytically, this growth in spin rate has been found to begin at the instant r/R = 0.707 [31]. The transverse angular speed, on the other hand, can either grow, decay, or fluctuate depending on the shape of the cylinder before the start of the burn. The corresponding body and space surfaces for the radially burning cylinder are shown in Figs. 12 and 13, respectively, for a burn duration of 90 s; this duration is chosen purely to study a case where the end mass of the system approaches the mass of a payload. The case presented here is for R/ h = 2 and is a nutationally stable configuration; this is in agreement with Mao and Eke's work [31] where it was shown that transverse rates are unbounded for oblate radially burning cylinders with a R/ h ≥ √ 8/3. The cone angle does in fact monotonically decrease though this is not as evident from Fig. 12 when compared to the uniform and end burns; this is primarily due to the fact that the system's end mass is approximately 315 kg at the  end of this radial burn, whereas in the other burns it is nearly zero. One approach to verifying the nutational stability of the system is by plotting the time evolution of θ . However, the stability can also be verified from the angular speeds. For the extreme case of oblateness represented by a radially burning flat disk (i.e., R >> h), Equation (68) can be reformulated as by dropping the terms involving h from the solution to the transverse speed. Extending this rationale to the inertia scalars given by Eqs. (64) and (65), we also get I /J = 1/2. Then, the nutation angle from Eq. (40) can be expressed as which is clearly bounded. In other words, even the most oblate configuration of a radially burning cylinder results in a constant nutation angle and thus does not show growth in the coning motion. However, unlike the preceding burns, the radial burn would require a nutation damper system to attenuate the cone angle and bring the system into a state of pure spin.

Conclusion
A graphical approach to study the coning motion, also known as the geometry of motion, of a torquefree, axisymmetric variable mass system is developed. Based on this, a general solution to the second Euler angle is also presented. The theoretical development is seen to augment well with the work on constant mass rigid bodies. Following these developments, the coning motion is specifically explored for an axisymmetric cylinder with mass loss subject to three idealized burns: the uniform, end, and radial burns. The uniform and end-burning cylindrical configurations exhibit constant spin rates and bounded decaying transverse rates; they are seen to be nutationally stable. The radial burn generally exhibits unbounded spin rates, whereas the transverse rate is unbounded only for more oblate configurations. More interestingly, even these oblate configurations are seen to show stable coning motion. In relation to rocket dynamics, the results presented here show that oblate configurations that exhibit a "misbehavior" in their transverse rates could still be nutationally stable. The difference between the radial burn from the other burns examined is that the system does not always approach a state of pure spin. This work recommends that examining the transverse speed alone may be insufficient in determining the stability of a freely rotating variable mass system. The rotational behavior of a free axisymmetric variable mass system is remarkably different from its constant mass counterpart. For a free axisymmetric system with mass variation, spin rates are functions of time that either grow, decay, or remain constant whereas their transverse speeds can either grow or decay; on the other hand, the spin and transverse speeds of constant mass systems are of constant magnitude. As the inertia properties in the latter case are also constant, the corresponding nutation angle (or second Euler angle) is also unchanging which leads to the well known 'cones of motion' of the angular velocity. In the case of variable mass systems, the nutation angle is not as straightforward to generally characterize on account of the timevarying nature of both the angular velocity and its central inertia scalars. Thus, the results of the analyses presented clearly indicate that mass variation majorly impacts the dynamics of freely spinning systems.