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

Nonlinear Dynamics, Jul 2018

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.

A PDF file should load here. If you do not see its contents the file may be temporarily unavailable at the journal website or you do not have a PDF plug-in installed and enabled in your browser.

Alternatively, you can download the file locally and open with any standalone PDF reader:

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

Nonlinear Dynamics pp 1–14 | Cite as Geometry of motion and nutation stability of free axisymmetric variable mass systems AuthorsAuthors and affiliations Angadh Nanjangud Open Access Original Paper First Online: 26 July 2018 Received: 16 April 2018 Accepted: 17 July 2018 141 Downloads Abstract 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. KeywordsMass variation Spacecraft attitude dynamics Classical mechanics  List of symbols B Rigid base of system \(\mathbf{b}_i (i=1,2,3)\) Dextral unit vectors attached to B C Rigid massless control volume F Fluid phase of variable mass system \(\mathbf{H}^*\) Central angular momentum of variable mass system h Initial half length of cylinder \(\mathbf{I}^*\) Central moment of inertia of the variable mass system I Transverse moment of inertia scalar J Spin moment of inertia scalar m Instantaneous mass of the system \(\mathbf{M}^*\) Resultant external moment about \(S^*\) \(\mathbf n\) Unit normal directed outward from the exit face of C \(\mathbf{n}_i (i= f,g,h)\) Dextral unit vectors attached to an inertial frame O Point on B \(\mathbf{p}\) Position vector from \({S}^*\) to P P A generic particle within C \(\mathbf{r}\) Position vector from O to P \(\mathbf{r}^*\) Position vector from O to \({S}^*\) R Radius of exit plane, and initial radius of cylindrical body \(\text {S}^*\) Mass center of system \(V_0\) Volume of C \(\mathbf{v}\) Inertial velocity of P \(\mathbf{v}_r\) Velocity of P relative to B z Instantaneous half length of cylinder \(z_e\) Shortest distance between \({S}^*\) and exit \({\varvec{\alpha }}\) Angular acceleration of the system \({\varvec{\omega }}\) Angular velocity of the system \(\rho \) Density of matter within C 1 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 nearly-symmetric 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. 2 Equations of motion and the angular velocity of an axisymmetric systems Figure 1 is of a general variable mass system. It comprises a consumable rigid base B and a fluid phase F. 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 $$\begin{aligned} \mathbf{M}^*= & {} \mathbf{I}^* \cdot {\varvec{\alpha }} + {\varvec{\omega }} \times (\mathbf{I}^* \cdot {\varvec{\omega }}) + \frac{^B\mathrm{d} \mathbf{I}^*}{\mathrm {d}t} \cdot {\varvec{\omega }}\nonumber \\&+ \int _{S_0} \rho [\mathbf{p} \times ({\varvec{\omega }} \times \mathbf{p})] (\mathbf{v}_\mathbf{r} \cdot {\mathbf{n}}) \, \mathrm{d}S \nonumber \\&+ \int _{V_0} \rho [{\varvec{\omega }} \times (\mathbf{p} \times \mathbf{v}_\mathbf{r})] \, \mathrm{d}V \nonumber \\&+ \frac{^B\mathrm{d}}{\mathrm {d}t} \int _{V_0}\rho (\mathbf{p} \times \mathbf{v}_\mathbf{r}) \, \mathrm{d}V\nonumber \\&+ \int _{S_0} \rho (\mathbf{p} \times \mathbf{v}_\mathbf{r})(\mathbf {v}_\mathbf{r} \cdot {\mathbf n}) \, \mathrm{d}S. \end{aligned}$$ (1) 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 of attitude motion becomes $$\begin{aligned} \mathbf{M}^*= & {} \mathbf{I}^* \cdot {\varvec{\alpha }} + {\varvec{\omega }} \times (\mathbf{I^*} \cdot {\varvec{\omega }}) + \frac{^B\mathrm{d} \mathbf{I}^*}{\mathrm {d}t} \cdot {\varvec{\omega }}\nonumber \\&+ \int _{S_0} \rho {\bigg [\mathbf{p} \times ({\varvec{\omega }} \times \mathbf{p})\bigg ] (\mathbf{v}_\mathbf{r} \cdot {\mathbf{n}})} \, \mathrm{d}S. \end{aligned}$$ (2) Open image in new window Fig. 1 General variable mass system The remainder of this paper concerns the torque-free motions of axisymmetric variable mass systems so \(\mathbf {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, \(\mathbf{b}_1, \mathbf{b}_2, \mathbf{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 $$\begin{aligned} \mathbf{I}^* \triangleq I \, \, \mathbf{b}_1 \mathbf{b}_1 + I \, \, \mathbf{b}_2 \mathbf{b}_2 + J \, \, \mathbf{b}_3 \mathbf{b}_3 \end{aligned}$$ (3) 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 \(\mathbf{b}_3\). Further, the inertial angular velocity of B at any instant is $$\begin{aligned} {\varvec{\omega }} \triangleq \omega _1 \mathbf{b}_1 + \omega _2 \mathbf{b}_2 + \omega _3 \mathbf{b}_3 \end{aligned}$$ (4) and, as a result, the inertial angular acceleration is $$\begin{aligned} {\varvec{\alpha }} = \dot{\omega }_1 \mathbf{b}_1 +\dot{\omega }_2 \mathbf{b}_2 +\dot{\omega }_3 \mathbf{b}_3. \end{aligned}$$ (5) The first three terms on the right-hand side of Eq. (2) evaluate to $$\begin{aligned}&\mathbf{I}^* \cdot {\varvec{\alpha }} = I \dot{\omega }_1\mathbf{b}_1 + I\dot{\omega }_2\mathbf{b}_2 + J \dot{\omega }_3 \mathbf{b}_3, \end{aligned}$$ (6) $$\begin{aligned}&{\varvec{\omega }} \times (\mathbf{I}^* \cdot {\varvec{\omega }}) = \bigg (J - I\bigg )\omega _2 \omega _3 \mathbf{b}_1 - \bigg (J - I\bigg ) \omega _1 \omega _3\mathbf{b}_2, \end{aligned}$$ (7) and $$\begin{aligned} \frac{^B\mathrm{d} \mathbf{I}^*}{\mathrm{d}t} \cdot {\varvec{\omega }} = \dot{I} \omega _1 \mathbf{b}_1 + \dot{I} \omega _2\mathbf{b}_2 + \dot{J} \omega _3 \mathbf{b}_3, \end{aligned}$$ (8) respectively. The assumed steady symmetry of the internal fluid flow is further constrained to be uniform along the exit plane so that \(\mathbf{v_r \cdot n}= u = {\text{ constant }}\), which then yields the following expression for the surface integral in Eq. (2) $$\begin{aligned}&\int _{S_o} \rho \mathbf{{p \times ({\varvec{\omega }} \times p)(v_r \cdot {n})}} \, \mathrm{d}S \nonumber \\&\quad = -\dot{m} \bigg [\left( z_e^2 + \frac{R^2}{4}\right) (\omega _1 \mathbf{b_1} + \omega _2 \mathbf{b_2}) + \frac{R^2}{2}\omega _3\mathbf{b_3}\bigg ]. \end{aligned}$$ (9) where \(\dot{m}\), the rate of mass loss, is given by $$\begin{aligned} \dot{m} = -\rho _\mathrm{exit} \pi R^2 u = \hbox {constant}. \end{aligned}$$ (10) \(\rho _\mathrm{exit}\) is the assumed density of exhaust gases on the exit plane. Substituting Eqs. (6), (7), (8), and (9) into Eq. (2) yields $$\begin{aligned} \dot{\omega }_1= & {} \frac{I - J}{I} \omega _2\omega _3 - \frac{1}{I}\bigg [\dot{I} - \dot{m}\left( z_e^2 + \frac{R^2}{4}\right) \bigg ] \omega _1 \end{aligned}$$ (11) $$\begin{aligned} \dot{\omega }_2= & {} -\frac{I - J}{I} \omega _3\omega _1 - \frac{1}{I}\bigg [\dot{I} - \dot{m}\left( z_e^2 + \frac{R^2}{4}\right) \bigg ] \omega _2 \end{aligned}$$ (12) and $$\begin{aligned} \dot{\omega }_3 = - \frac{1}{J}\left[ \dot{J} - \dot{m}\left( \frac{R^2}{2}\right) \right] \omega _3. \end{aligned}$$ (13) 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; \(\omega _1\) and \(\omega _2\) are transverse rates and \(\omega _3\) is the spin rate. The attitude motion of a variable system is characterized by solving these equations. 2.1 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 $$\begin{aligned} \omega _3 = \omega _{30} \text {exp}(\varPhi ) \end{aligned}$$ (14) where $$\begin{aligned} \varPhi = - \int _0^t \bigg [\bigg (\dot{J} - \dot{m} R^2/2\bigg )/J \bigg ] \mathrm {d}t. \end{aligned}$$ (15) If \(J \triangleq m k_3^2\), where \(k_3\) is the time-varying axial radius of gyration, then \(\varPhi \) evaluates to $$\begin{aligned} \varPhi = \int _{m(0)}^{m} \frac{(R^2/2)}{k_3^2} \frac{\mathrm {d}m}{m} + \mathrm {ln}\bigg (\frac{J(0)}{J}\bigg ). \end{aligned}$$ (16) Substituting for \(\varPhi \) in Eq. (14) gives the following expression for the spin rate $$\begin{aligned} \omega _3 = \omega _{30} \frac{J(0)}{J} \text {exp}\bigg (\int _{m(0)}^{m} \frac{(R^2/2)}{k_3^2} \frac{\mathrm {d}m}{m} \bigg ) \end{aligned}$$ (17) 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_3^2\) is assumed to be a constant in Eq.  (17), then $$\begin{aligned} \omega _3 = \omega _{30} \bigg ( \frac{m(0)}{m}\bigg ) ^{1 - \frac{R^2}{2k_3^2}}. \end{aligned}$$ (18) 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]. 2.2 Transverse rate The transverse angular speeds can be solved for in a variety of different ways. Defining a complex variable \(\omega _c\) as $$\begin{aligned} \omega _c \triangleq \omega _1 + j \omega _2 \end{aligned}$$ (19) allows Eqs. (11) and (12) to be combined to yield a differential equation linear in \(\omega _c\): $$\begin{aligned} \dot{\omega }_c = \Bigg \{-j \frac{I - J}{I} \omega _3 - \frac{1}{I}\bigg [\dot{I} - \dot{m}\left( z_e^2 + \frac{R^2}{4}\right) \bigg ]\Bigg \} \omega _c. \end{aligned}$$ (20) Integrating Eq. (20) yields $$\begin{aligned} {\omega }_c = \omega _{co} \cdot \varGamma \cdot \text {exp} \bigg (-j\chi \bigg ) \end{aligned}$$ (21) where $$\begin{aligned} \varGamma = \text {exp} \bigg [ -\int _0^t \frac{1}{I}\bigg [\dot{I} - \dot{m}\left( z_e^2 + \frac{R^2}{4}\right) \bigg ] \mathrm {d}t \bigg ] \end{aligned}$$ (22) and $$\begin{aligned} \chi = \int _0^t \left( 1 - \frac{J}{I}\right) \omega _3 \mathrm {d}t. \end{aligned}$$ (23) A consequence of Eq. (23) is $$\begin{aligned} \dot{\chi } = \left( 1 - \frac{J}{I}\right) \omega _3 \end{aligned}$$ (24) where \(\dot{\chi }\) has the units of angular speed. From Eqs. (21) and (19), the transverse speeds are given by $$\begin{aligned} \omega _1= & {} \varGamma (\omega _{10} \cos \chi + \omega _{20}\sin \chi ) \end{aligned}$$ (25) $$\begin{aligned} \omega _2= & {} \varGamma (-\,\omega _{10}\sin \chi + \omega _{20}\cos \chi ) \end{aligned}$$ (26) and are oscillatory in nature. Further, a transverse angular velocity vector \({\varvec{\omega }}_{t}\) which lies in the plane made by the principal directions \(\mathbf{b}_1\) and \(\mathbf{b}_2\) can also be defined as $$\begin{aligned} {\varvec{\omega }}_{t} \triangleq \omega _1 \mathbf{b}_1+ \omega _2 \mathbf{b}_2, \end{aligned}$$ (27) whose magnitude is $$\begin{aligned} \omega _{12} = ||{\varvec{\omega }}_{t}||= \sqrt{\omega _1^2 + \omega _2^2} = \sqrt{\varGamma ^2 (\omega _{10}^2 +\omega _{20}^2)} = \varGamma \omega _0 \end{aligned}$$ (28) where \(\omega _0 \triangleq \sqrt{\omega _{10}^2 + \omega _{20}^2} = \text {constant}\). The principal directions may be chosen such that \(\omega _{10} = 0\) and \(\omega _{20} = \omega _0\). This simplifies the appearance of the solutions in Eqs. (25) and (26) to $$\begin{aligned} \omega _1= & {} \omega _0 \varGamma \sin \chi \end{aligned}$$ (29) $$\begin{aligned} \omega _2= & {} \omega _0 \varGamma \cos \chi . \end{aligned}$$ (30) Equation (27) can then be written as $$\begin{aligned} {\varvec{\omega }}_{t} = \omega _0 \varGamma [\sin \chi \mathbf{b}_1+ \cos \chi \mathbf{b}_2] = \omega _0 \varGamma \mathbf{b}_{12} = \omega _{12} \mathbf{b}_{12}, \end{aligned}$$ (31) where \(\mathbf{b}_{12} = \sin \chi \mathbf{b}_1+ \cos \chi \mathbf{b}_2\); \(\mathbf{b}_{12}\) is clearly a unit vector that is orthogonal to \(\mathbf{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 \(\omega _0\) by \(\varGamma \), a function of time that can either grow, decay, or fluctuate. \(\mathbf{b}_{12}\) is a unit vector that rotates in the plane of \(\mathbf{b}_1\) and \(\mathbf{b}_2\) at the angular speed \(\dot{\chi }\). As revealed by Eq. (24), \(\dot{\chi }\) may be positive or negative (assuming that \(\omega _3\) is always positive), depending on the relative values of \(k_1\) and \(k_3\); if \(k_3 > k_1\), then \(\mathbf{b}_{12}\) rotates negatively in the body, advancing from \(\mathbf{b}_1\) to \(-\,\mathbf{b}_2\) and around; if \(k_1 > k_3\), then \(\mathbf{b}_{12}\) rotates positively in the body, coinciding with \(\mathbf{b}_1\) and then with \(\mathbf{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 \(\chi \). Its value can be obtained numerically from Eq. (23). An expression for \(\varGamma \) can be developed involving the system’s transverse radius of gyration, \(k_1\). If \(I \triangleq m k_1^2\), \(\varGamma \) may be reformulated in Eq. (22) as $$\begin{aligned} \varGamma = \frac{I(0)}{I}\text {exp}\bigg [ \int _{m(0)}^m \frac{\bigg [z_e^2 + \frac{R^2}{4}\bigg ]}{k_1^2} \frac{\mathrm {d}m}{m}\bigg ] \end{aligned}$$ (32) 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 \(\varGamma \) into Eq. (28) yields the following expression for the magnitude of the transverse angular velocity: $$\begin{aligned} \omega _{12} = \omega _{0} \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\bigg [ \int _{m(0)}^m \frac{\bigg [z_e^2 + \frac{R^2}{4}\bigg ]}{k_1^2} \frac{\mathrm {d}m}{m}\bigg ]. \end{aligned}$$ (33) 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 \(\omega _{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: $$\begin{aligned} {\varvec{\omega }} = \omega _0 \varGamma \mathbf{b}_{12} + \omega _3 \mathbf{b}_3, \end{aligned}$$ (34) the form of which mimics the constant mass rigid body case, for which \(\varGamma = 1\). 3 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 body-fixed 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 as $$\begin{aligned} \mathbf{H}^* = \int _{V_0} \rho [\mathbf{p} \times [\mathbf{v}_r + ({\varvec{\omega }} \times \mathbf{p}) ]\, \mathrm {d}V. \end{aligned}$$ (35) Retaining the assumption of steady and axisymmetric relative internal mass flow reduces Eq. (35) to $$\begin{aligned} \mathbf{H}^* = \int _{V_0} \rho [\mathbf{p} \times ({\varvec{\omega }} \times \mathbf{p}) ] \, \mathrm {d}V = \mathbf{I}^* \cdot {\varvec{\omega }}. \end{aligned}$$ (36) With the definitions of \(\mathbf{I}^*\) and \(\varvec{\omega }\) from Eqs. (3) and (4), respectively, the component form of the angular momentum vector is $$\begin{aligned} \mathbf{H}^* = I (\omega _1 \mathbf{b}_1 + \omega _2\mathbf{b}_2) + J \omega _3 \mathbf{b}_3 \end{aligned}$$ (37) or $$\begin{aligned} \mathbf{H}^* = I {\varvec{\omega }}_{t} + J \omega _3 \mathbf{b}_3. \end{aligned}$$ (38) Using the component form of \({\varvec{\omega }}_{t}\) from Eq. (31) in Eq. (38) gives $$\begin{aligned} \mathbf{H}^* = I \omega _0 \varGamma \mathbf{b}_{12} + J \omega _3 \mathbf{b}_3 \end{aligned}$$ (39) Open image in new window Fig. 2 Geometry of motion setup of vectors 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); \(\mathbf{H}^*\) and \({\varvec{\omega }}\) are in the same plane formed by \(\mathbf{b}_{12}\) and \(\mathbf{b}_3\) and this is shown in Fig. 2. In this figure, \(\theta \) is the angle between \(\mathbf{H}^*\) and \(\mathbf{b}_3\), and \(\beta \) is the angle between \(\mathbf{b}_3\) and \(\varvec{\omega }\). These angles can be computed from the figure as $$\begin{aligned} \theta = \tan ^{-1}\bigg (\frac{I \omega _{12}}{J \omega _3} \bigg ) = \tan ^{-1}\bigg (\frac{\varGamma k_1^2 \omega _{0}}{k_3^2 \omega _3} \bigg ) \end{aligned}$$ (40) and $$\begin{aligned} \beta = \tan ^{-1}\bigg (\frac{\omega _{12}}{\omega _3} \bigg ) = \tan ^{-1} \bigg (\frac{\varGamma \omega _{0}}{\omega _3} \bigg ). \end{aligned}$$ (41) Evaluation of these angles \(\theta \) and \(\beta \) 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, \(\varGamma = 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 \({\varvec{\omega }}\) results in the formation of two cones: as it rotates about \(\mathbf{b}_3\) and another as it rotates about \(\mathbf{n}_h\), a unit vector parallel to the inertially fixed \(\mathbf{H}^*\). The two cones are appropriately named the body cone and space cone, respectively. In the case of axisymmetric variable mass systems, angles \(\theta \) and \(\beta \) are typically functions of time. The angles may grow or decay as dictated by the instantaneous values of the radii of gyration, \(\varGamma \), and the spin rate. It is also possible, in the case of the same variable mass system, for \(\beta \) to exceed \(\theta \) at certain instants, while it may not at certain other instants. Figure 2 represents an instant where \(\beta >\theta \), which occurs when \(k_3>k_1\). Also, \({\varvec{\omega }}\) does not trace cones as it rotates about \(\mathbf{b}_3\) and \(\mathbf{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, \(\theta \) is referred to as the nutation angle [18]. In classical mechanics, \(\theta \) 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 \(\theta \) 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 closed-form 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. 4 The axisymmetric, variable mass cylinder Open image in new window Fig. 3 Variable mass cylinder 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. 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 \(\omega _0=0.2\), and \(\omega _{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. 4.1 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 $$\begin{aligned} I = m \bigg (\frac{R^2}{4} + \frac{h^2}{3}\bigg ) \end{aligned}$$ (42) and $$\begin{aligned} J = m\frac{R^2}{2}. \end{aligned}$$ (43) This is an example of a variable mass system where the radii of gyration are constant where \(k_1^2 = R^2/4 + h^2/3\) and \(k_3^2 = R^2/4\). The time derivatives of the moment of inertia are $$\begin{aligned} \begin{aligned} \dot{I}&= \dot{m}\left( \frac{R^2}{4} + \frac{h^2}{3}\right) \\ \dot{J}&= \dot{m} \frac{R^2}{2} \end{aligned} \end{aligned}$$ (44) The resulting spin and transverse rates are obtained from Eqs. (17) and (28) as $$\begin{aligned} \omega _3= & {} \omega _{30} \end{aligned}$$ (45) $$\begin{aligned} \omega _{12}= & {} \omega _0 \bigg (\frac{m}{m_0}\bigg )^{\frac{2h^2}{3k_1^2}} = \omega _0 \varGamma (t) \end{aligned}$$ (46) where \(\varGamma (t) \triangleq (m/m_0)^{2h^2/3k_1^2}\). Prior to the start of the burn \(m = m_0\) so \(\varGamma = 1\). As the burn proceeds, \(\varGamma \) 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 $$\begin{aligned} \theta (t) = \tan ^{-1}( K \varGamma (t) ), \end{aligned}$$ (47) where \(K = k_1^2 \omega _0/k_3^2 \omega _{30}\). Since K is constant and \(\varGamma (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 \(\varvec{\omega }\) about \(\mathbf{b}_3\) and \(\mathbf{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 \(\mathbf{b}_3\) in a counterclockwise direction and eventually converges to the symmetry axis. Open image in new window Fig. 4 Uniform burn: body surface Open image in new window Fig. 5 Uniform burn: space surface The rotation of \(\omega \) in the inertial space (i.e., about \(\mathbf{n}_h\) since it is inertially fixed) requires knowledge of the angle, \(\beta \), in Eq. (41) $$\begin{aligned} \beta = \tan ^{-1}\bigg (\frac{\omega _{0}}{\omega _{30}} \bigg (\frac{m}{m(0)}\bigg )^{\frac{2h^2}{3k_1^2}} \bigg ). \end{aligned}$$ (48) 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 \(\theta \), and \(\beta \) along with the angular velocity permit the construction of the space surface as shown in Fig. 5. Note that \(\mathbf{n}_f\), \(\mathbf{n}_g\), and \(\mathbf{n}_h\) represent a dextral set of unit vectors which are inertially fixed. 4.2 End burn Open image in new window Fig. 6 End-burning cylinder 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, given by \(k_3^2 = R^2/4\). The corresponding axial moment of inertia is $$\begin{aligned} J = m R^2/4 \quad \end{aligned}$$ (49) As in the case of uniform burn, the spin rate can thus be determined from Eq. (18) to be constant, $$\begin{aligned} \omega _3 = \omega _{30}. \end{aligned}$$ (50) The more interesting development occurs in the transverse directions. The transverse radius of gyration is \(k_1^2 = \frac{R^2}{4} + \frac{z^2}{3}\). Due to the changing length of the cylinder, \(k_1\) is not a constant. Then, the transverse moment of inertia is $$\begin{aligned} I = m\left( \frac{R^2}{4} + \frac{z^2}{3}\right) \end{aligned}$$ (51) If z is the instantaneous half length of the unburned cylinder, then its instantaneous mass is given by $$\begin{aligned} m = \pi R^2 (2z) \end{aligned}$$ (52) The rate of mass loss is then $$\begin{aligned} \dot{m} = 2 \pi R^2 \dot{z}. \end{aligned}$$ (53) which can also be written as $$\begin{aligned} \mathrm {d}m = 2 \pi R^2 \mathrm {d}z. \end{aligned}$$ (54) From Eqs. (52) and (54), it is possible to recast Eq. (33) as $$\begin{aligned} \omega _{12} = \omega _0 \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\left[ \int _{z_0}^z \frac{\bigg [z_e^2 + \frac{R^2}{4}\bigg ]}{\frac{R^2}{4} + \frac{z^2}{3}} \frac{\mathrm {d}z}{z}\right] . \end{aligned}$$ (55) Substituting \(z_e = 2h - z\) into Eq. (55) gives $$\begin{aligned} \omega _{12} = \omega _0 \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\left[ \int _{z_0}^z \frac{\left[ z^2 + 4h^2 -4hz + \frac{R^2}{4}\right] }{\frac{R^2}{4} + \frac{z^2}{3}} \frac{\mathrm {d}z}{z}\right] \end{aligned}$$ (56) which can also be reformulated as $$\begin{aligned} \omega _{12}= & {} \omega _0 \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\left[ \int _{z_0}^z \frac{3z - 12h}{\frac{3R^2}{4} + z^2} \mathrm {d}z \right. \nonumber \\&\left. + \bigg (12h^2 + \frac{3R^2}{4} \bigg )\int _{z_0}^z \frac{\mathrm {d}z}{z\bigg (\frac{3R^2}{4} + z^2\bigg )} \right] . \end{aligned}$$ (57) A partial fraction expansion of the second integral in Eq. (57) then gives $$\begin{aligned} \omega _{12}= & {} \omega _0 \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\bigg [ \int _{z_0}^z\frac{3z - 12h}{\frac{3R^2}{4} + z^2} \mathrm {d}z \nonumber \\&+ \bigg (\frac{16h^2}{R^2} + 1 \bigg ) \int _{z_0}^z \bigg (\frac{1}{z} - \frac{z}{\frac{3R^2}{4} + z^2}\bigg ) \mathrm {d}z \bigg ] \end{aligned}$$ (58) or, after some algebra, $$\begin{aligned} \omega _{12}= & {} \omega _0 \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\bigg [ -\int _{z_0}^z \frac{12h}{ \frac{3R^2}{4} + z^2 } \mathrm {d}z \nonumber \\&- \bigg (\frac{16h^2}{R^2} - 2 \bigg ) \int _{z_0}^z \frac{z}{\frac{3R^2}{4} + z^2} \mathrm {d}z\nonumber \\&+ \bigg (\frac{16h^2}{R^2} + 1 \bigg ) \int _{z_0}^z\frac{\mathrm {d}z}{z} \bigg ]. \end{aligned}$$ (59) The solution to the transverse rate is $$\begin{aligned} \omega _{12}= & {} \omega _0 \bigg (\frac{I(0)}{I}\bigg ) \text {exp}\bigg [ \frac{8 \sqrt{3}h}{R} \tan ^{-1}\bigg ( \frac{2\sqrt{3}R(z-h)}{3R^2 + 4zh} \bigg ) \nonumber \\&- \bigg (\frac{8h^2}{R^2} - 1 \bigg ) \log \bigg ( \frac{\frac{3R^2}{4} + z^2}{\frac{3R^2}{4} + h^2}\bigg ) \nonumber \\&+ \bigg (\frac{16h^2}{R^2} + 1 \bigg ) \log {\frac{z}{h}} \bigg ]. \end{aligned}$$ (60) which can also be expressed as $$\begin{aligned} \omega _{12} =&\omega _0 \bigg (\frac{I(0)}{I}\bigg ) \bigg ( \frac{\frac{3R^2}{4} + h^2}{\frac{3R^2}{4} + z^2}\bigg )^{ \bigg (\frac{8h^2}{R^2} - 1 \bigg )} \bigg (\frac{z}{h}\bigg )^ {\bigg (\frac{16h^2}{R^2} + 1 \bigg )} \nonumber \\&\text {exp}\bigg [\frac{8 \sqrt{3}h}{R} \tan ^{-1}\bigg ( \frac{2\sqrt{3}R(z-h)}{3R^2 + 4zh} \bigg ) \bigg ]. \end{aligned}$$ (61) The above solution can also be expressed in terms of the transverse radius of gyration as $$\begin{aligned} \omega _{12} = \omega _0&\bigg (\frac{k_1(0)}{k_1}\bigg )^{\frac{8h^2}{R^2}} \bigg (\frac{z}{h}\bigg )^{\frac{16h^2}{R^2}}\nonumber \\&\text {exp}\bigg [\frac{8 \sqrt{3}h}{R} \tan ^{-1}\bigg ( \frac{2\sqrt{3}R(z-h)}{3R^2 + 4zh} \bigg ) \bigg ]. \end{aligned}$$ (62) Equation (62) explains that, for an end-burning cylinder, the transverse rate is regulated by three time-varying 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 $$\begin{aligned} \theta (t) = \tan ^{-1}( K k_1^2 \varGamma (t) ), \end{aligned}$$ (63) where \(K = \omega _0/k_3^2 \omega _{30}=\) constant. As both \(\varGamma (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 \(\mathbf{b}_3\) and \(\mathbf{n}_h\). Initially, it rotates in a clockwise direction and then switches to a counterclockwise direction, eventually converging to a pure spin about \(\mathbf{b}_3\) in the body frame and \(\mathbf{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 \(\dot{\chi }\) 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. Open image in new window Fig. 7 End burn: body surface (\(J_0>I_0\)) Open image in new window Fig. 8 End burn: space surface (\(J_0>I_0\)) Open image in new window Fig. 9 End burn: body surface (\(J_0<I_0\)) Open image in new window Fig. 10 End burn: space surface (\(J_0<I_0\)) Open image in new window Fig. 11 Radially burning cylinder 4.3 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 radially burning cylinder are $$\begin{aligned} I = m\left( \frac{R^2 + r^2}{4} + \frac{h^2}{3}\right) \end{aligned}$$ (64) and $$\begin{aligned} J = m\frac{R^2 + r^2}{2}. \end{aligned}$$ (65) Their corresponding time derivatives are $$\begin{aligned} \begin{aligned} \dot{I}&= \dot{m}\left( \frac{R^2}{2} + \frac{h^2}{3}\right) \\ \dot{J}&= \dot{m} {r^2}. \end{aligned} \end{aligned}$$ (66) The angular speeds evaluate to $$\begin{aligned} \omega _3 = \omega _{30} \frac{R^4}{(R^2 + r^2)\sqrt{(R^2 + r^2)(R^2 - r^2)}} \end{aligned}$$ (67) $$\begin{aligned}&\omega _{12} = \omega _0 \bigg (\frac{R^2 + {4h^2}/3}{R^2 + {4h^2/3 + r^2}} \bigg )^\frac{(3R^2 + 16h^2/3)}{(2R^2 + 4h^2/3)}\nonumber \\&\quad \qquad \qquad \bigg (\frac{R^2 - r^2}{R^2} \bigg )^\frac{(-R^2 + 8h^2/3)}{(2R^2 + 4h^2/3)}. \end{aligned}$$ (68) 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. Open image in new window Fig. 12 Radial burn: body surface Open image in new window Fig. 13 Radial burn: space surface 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 \ge \sqrt{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 \(\theta \). 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 $$\begin{aligned} \omega _{12} = \omega _0 \frac{\omega _3}{\omega _{30}} \end{aligned}$$ (69) 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 $$\begin{aligned} \theta = \tan ^{-1}\bigg (\frac{\omega _0}{2\omega _{30}} \bigg ) = \text {constant}. \end{aligned}$$ (70) 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. 5 Conclusion A graphical approach to study the coning motion, also known as the geometry of motion, of a torque-free, 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 time-varying 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. Notes Acknowledgements The author would like to acknowledge the University of California Davis’ Drake Fellowship award and the Mechanical and Aerospace Engineering Department’s N. & M. SarigulKlijn Space Engineering/Flight Research award that supported this work. The author would also like to thank Prof. Fidelis Eke and A Elbakyan for supporting his research efforts. Compliance with ethical standards Conflict of interest The authors declare that they have no conflict of interest. References 1. Barkley, R.J., Newton, R.R., Gross, G.L.: Mathematical Theory of Rocket Flight. McGraw-Hill, New York (1947)Google Scholar 2. Reiter, G., Thomson, W.: Jet damping of a solid rocket-theory and flight results. AIAA J. 3, 413417 (1965)Google Scholar 3. Thomson, W.T.: Equations of motion for the variable mass system. AIAA J. 4, 766768 (1966)Google Scholar 4. Thomson, W.T.: Introduction to Space Dynamics. Dover, New York (1986)Google Scholar 5. Nanjangud, A., Eke, F.O.: Lagrange’s equations for rocket-type variable mass systems. Int. Rev. Aerosp. Eng. 5(5), 256–260 (2012)Google Scholar 6. Djerassi, S.: Algorithm for simulation of motions of variable-mass systems. J. Guid. Control Dyn. 21, 427434 (1998)Google Scholar 7. Banerjee, A.K.: Dynamics of a variable-mass, flexible-body system. J. Guid. Control Dyn. 23, 501–508 (2000)CrossRefGoogle Scholar 8. Hu, P., Ren, G.: Multibody dynamics of flexible liquid rockets with depleting propellant. J. Guid. Control Dyn. 36, 1840–1849 (2013)CrossRefGoogle Scholar 9. Cveticanin, L., Djukic, D.: Motion of body with discontinual mass variation. Nonlinear Dyn. 52(3), 249–261 (2008). CrossRefMATHGoogle Scholar 10. Cveticanin, L.: Dynamics of body separationanalytical procedure. Nonlinear Dyn. 55(3), 269–278 (2007)MathSciNetCrossRefGoogle Scholar 11. Eke, F.O., Wang, S.-M.: Attitude behavior of a variable mass cylinder. J. Appl. Mech. 62(4), 935–940 (1995)CrossRefMATHGoogle Scholar 12. Wang, S.-M., Eke, F.O.: Rotational dynamics of axisymmetric variable mass systems. J. Appl. Mech. 62(4), 970–974 (1995)CrossRefMATHGoogle Scholar 13. Eke, F.O., Mao, T.-C., Morris, M.J.: Free attitude motions of a spinning body with substantial mass loss. J. Appl. Mech. 71(2), 190–194 (2004)CrossRefMATHGoogle Scholar 14. Flandro, G.A., Van Moorhem, W.K., Shorthill, R., Chen, K., Woolsey, M.: Fluid Mechanics of Spinning Rockets. U.S. Air Force Rocket Propulsion Lab. Rept. TR-86-072, Edwards AFB, CA (1987)Google Scholar 15. Ha, J.C.V.D., Janssens, F.L.: Jet-damping and misalignment effects during solid rocket motor burn. J. Guid. Control Dyn. 28, 412420 (2005)Google Scholar 16. Nanjangud, A., Eke, F.O.: Approximate solution to the angular speeds of a nearly-symmetric mass-varying cylindrical body. J. Astronaut. Sci. 64(2), 99–117 (2017)CrossRefGoogle Scholar 17. Likins, P.W.: Elements of Engineering Mechanics. McGraw-Hill, New York (1973)MATHGoogle Scholar 18. Kaplan, M.H.: Modern Spacecraft Dynamics and Control. Wiley, New York (1976)Google Scholar 19. Ayoubi, M.A., Martin, K.M., Longuski, J.M.: Analytical solution for spinning thrusting spacecraft with transverse ramp-up torques. J. Guid. Control Dyn. 37(4), 1272–1282 (2014)CrossRefGoogle Scholar 20. Martin, K.M., Longuski, J.M.: Velocity pointing error reduction for spinning, thrusting spacecraft via heuristic thrust profiles. J. Spacecr. Rockets 52(4), 1268–1272 (2015)CrossRefGoogle Scholar 21. Javorsek, D., Longuski, J.M.: Velocity pointing errors associated with spinning thrusting spacecraft. J. Spacecr. Rockets 37(3), 359–365 (2000)CrossRefGoogle Scholar 22. Nanjangud, A., Eke, F.O.: Angular momentum of free variable mass systems is partially conserved. Aerosp. Sci. Technol. (2018). Google Scholar 23. Nanjangud, A.: On the rotational dynamics of variable mass systems. Ph.D. dissertation, University of California, Davis (2016)Google Scholar 24. Snyder, V.W., Warner, G.G.: A re-evaluation of jet damping. J. Spacecr. Rockets 5, 364–366 (1968)CrossRefGoogle Scholar 25. Meyer, R.X.: Convective instability in solid propellant rocket motors. Adv. Astronaut. Sci. 54, 173 (1983). (AAS Paper 83-368)Google Scholar 26. Cochran Jr., J.E., Kang, J.Y.: Nonlinear stability analysis of the attitude motion of a spin-stabilized upper stage. Adv. Astronaut. Sci. 75(1), 345–364 (1991)Google Scholar 27. Or, A.C.: Rotor-pendulum model for the perigee assist module nutation anomaly. J. Guid. Control Dyn. 15(2), 297–303 (1992)CrossRefGoogle Scholar 28. Mingori, D., Yam, Y.: Nutational Stability of a Spinning Spacecraft with Internal Mass Motion and Axial Thrust. AIAA Paper 86-2271 AIAA Astrodynamics Conference Proceedings (Williamsburg, VA), pp. 367-375. AIAA, Washington, DC (1986)Google Scholar 29. Halsmer, D.M., Mingori, D.L.: Nutational stability and passive control of spinning rockets with internal mass motion. J. Guid. Control Dyn. 18(5), 1197–1203 (1995)CrossRefMATHGoogle Scholar 30. Webster, E.: Active Nutation Control for Spinning Solid Motor Upper Stages. American Institute of Aeronautics and Astronautics, Monterey, CA (1985)CrossRefGoogle Scholar 31. Mao, T.C., Eke, F.O.: Attitude dynamics of a torque-free variable mass cylindrical body. J. Astronaut. Sci. 48, 435–448 (2000)Google Scholar Copyright information © The Author(s) 2018 Open AccessThis article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. Authors and Affiliations Angadh Nanjangud1Email authorView author's OrcID profile1.University of California DavisDavisUSA

This is a preview of a remote PDF:

Angadh Nanjangud. Geometry of motion and nutation stability of free axisymmetric variable mass systems, Nonlinear Dynamics, 2018, 1-14, DOI: 10.1007/s11071-018-4485-6