Equations of Motion v0
Introduction
This document describes the equations of motion which govern the flight
of a rocket. The equations are implemented in Python to simulate the
dynamics of a rocket flight, which can be used to determine its
stability, its dispersion as well as the altitude which the rocket
theoretically reaches.
We start by defining the rocket and the reference systems which will be
used throughout the derivation. Consider the following diagram:
The rocket is considered a composition of a rigid body (\(R\)), a
temporary solid phase (\(T\)) and a gas phase (\(G\)). Fixed to
the permanent rigid body \(R\) is a coordinate system \(B\)
whose axes are oriented considering a right-handed orthogonal base
\(b_{1}b_{2}b_{3}\). An inertial coordinate system \(A\) is also
defined. \(A\) is fixed to Earth, which is assumed inertial, and is
also oriented by a right-handed orthogonal base \(a_{1}a_{2}a_{3}\).
Velocity, Angular Velocity, and Reference Frame Transformation
Consider the center of mass \(*\) of \(R + T + G\). Its position
and velocity relative to the inertial reference frame \(A\) is
written as:
\[_{\ }^{A}{\overrightarrow{r}}_{*} = x{\overrightarrow{a}}_{1} + y{\overrightarrow{a}}_{2} + z{\overrightarrow{a}}_{3}\]
\[_{\ }^{A}{\overrightarrow{v}}_{*} = \frac{\text{dx}}{\text{dt}}{\overrightarrow{a}}_{1} + \frac{\text{dy}}{\text{dt}}{\overrightarrow{a}}_{2} + \frac{\text{dz}}{\text{dt}}{\overrightarrow{a}}_{3}\]
One may also choose to represent this in terms of
\(b_{1}b_{2}b_{3}\). To convert between both reference frames, a
transformation matrix shall be used. This transformation can be
represented in terms of Euler parameters which will help reach a concise
form of the equations of motion later. More about Euler parameters can
be seen in Analytical Mechanics by Haim Baruh in section 7.7. The
transformation is given by:
\[\begin{split}\begin{bmatrix}
{\overrightarrow{b}}_{1} \\
{\overrightarrow{b}}_{2} \\
{\overrightarrow{b}}_{3} \\
\end{bmatrix} = \begin{bmatrix}
e_{0}^{2} + e_{1}^{2} - e_{2}^{2} - e_{3}^{2} & 2\Bigl( e_{1}e_{2} + e_{0}e_{3} \Bigr) & 2(e_{1}e_{3} - e_{0}e_{2}) \\
2\Bigl( e_{1}e_{2} - e_{0}e_{3} \Bigr) & e_{0}^{2} - e_{1}^{2} + e_{2}^{2} - e_{3}^{2} & 2\Bigl( e_{2}e_{3} + e_{0}e_{1} \Bigr) \\
2(e_{1}e_{3} + e_{0}e_{2}) & 2\Bigl( e_{2}e_{3} - e_{0}e_{1} \Bigr) & e_{0}^{2} - e_{1}^{2} - e_{2}^{2} + e_{3}^{2} \\
\end{bmatrix}\begin{bmatrix}
{\overrightarrow{a}}_{1} \\
{\overrightarrow{a}}_{2} \\
{\overrightarrow{a}}_{3} \\
\end{bmatrix}\end{split}\]
The reverse transformation is found by inverting the rotation matrix,
which is the same as transposing it:
\[\begin{split}\begin{bmatrix}
{\overrightarrow{a}}_{1} \\
{\overrightarrow{a}}_{2} \\
{\overrightarrow{a}}_{3} \\
\end{bmatrix} = \begin{bmatrix}
e_{0}^{2} + e_{1}^{2} - e_{2}^{2} - e_{3}^{2} & 2\Bigl( e_{1}e_{2} - e_{0}e_{3} \Bigr) & 2(e_{1}e_{3} + e_{0}e_{2}) \\
2\Bigl( e_{1}e_{2} + e_{0}e_{3} \Bigr) & e_{0}^{2} - e_{1}^{2} + e_{2}^{2} - e_{3}^{2} & 2\Bigl( e_{2}e_{3} - e_{0}e_{1} \Bigr) \\
2(e_{1}e_{3} - e_{0}e_{2}) & 2\Bigl( e_{2}e_{3} + e_{0}e_{1} \Bigr) & e_{0}^{2} - e_{1}^{2} - e_{2}^{2} + e_{3}^{2} \\
\end{bmatrix}\begin{bmatrix}
{\overrightarrow{b}}_{1} \\
{\overrightarrow{b}}_{2} \\
{\overrightarrow{b}}_{3} \\
\end{bmatrix}\end{split}\]
Therefore, point O’s velocity relative to \(A\) expressed in
\(B\) is given by:
\[_{\ }^{A}{\overrightarrow{v}}_{*} = \frac{\text{dx}}{\text{dt}}\left\lbrack \Bigl( e_{0}^{2} + e_{1}^{2} - e_{2}^{2} - e_{3}^{2} \Bigr){\overrightarrow{b}}_{1} + \Bigl( 2\Bigl( e_{1}e_{2} - e_{0}e_{3} \Bigr) \Bigr){\overrightarrow{b}}_{2} + \Bigl( 2(e_{1}e_{3} + e_{0}e_{2}) \Bigr){\overrightarrow{b}}_{3} \right\rbrack + \frac{\text{dy}}{\text{dt}}\left\lbrack \Bigl( 2\Bigl( e_{1}e_{2} + e_{0}e_{3} \Bigr) \Bigr){\overrightarrow{b}}_{1} + \Bigl( e_{0}^{2} - e_{1}^{2} + e_{2}^{2} - e_{3}^{2} \Bigr){\overrightarrow{b}}_{2} + \Bigl( 2\Bigl( e_{2}e_{3} - e_{0}e_{1} \Bigr) \Bigr){\overrightarrow{b}}_{3} \right\rbrack + \frac{\text{dz}}{\text{dt}}\left\lbrack \Bigl( 2(e_{1}e_{3} - e_{0}e_{2}) \Bigr){\overrightarrow{b}}_{1} + \Bigl( 2\Bigl( e_{2}e_{3} + e_{0}e_{1} \Bigr) \Bigr){\overrightarrow{b}}_{2} + \Bigl( e_{0}^{2} - e_{1}^{2} - e_{2}^{2} + e_{3}^{2} \Bigr){\overrightarrow{b}}_{3} \right\rbrack\]
\[_{\ }^{A}{\overrightarrow{v}}_{*} = \Bigl( \frac{\text{dx}}{\text{dt}}\Bigl( e_{0}^{2} + e_{1}^{2} - e_{2}^{2} - e_{3}^{2} \Bigr) + \frac{\text{dy}}{\text{dt}}\Bigl( 2\Bigl( e_{1}e_{2} + e_{0}e_{3} \Bigr) \Bigr) + \frac{\text{dz}}{\text{dt}}\Bigl( 2(e_{1}e_{3} - e_{0}e_{2}) \Bigr) \Bigr){\overrightarrow{b}}_{1} + \Bigl( \frac{\text{dx}}{\text{dt}}\Bigl( 2\Bigl( e_{1}e_{2} - e_{0}e_{3} \Bigr) \Bigr) + \frac{\text{dy}}{\text{dt}}\Bigl( e_{0}^{2} - e_{1}^{2} + e_{2}^{2} - e_{3}^{2} \Bigr) + \frac{\text{dz}}{\text{dt}}\Bigl( 2\Bigl( e_{2}e_{3} + e_{0}e_{1} \Bigr) \Bigr) \Bigr){\overrightarrow{b}}_{2} + \Bigl( \frac{\text{dx}}{\text{dt}}\Bigl( 2(e_{1}e_{3} + e_{0}e_{2}) \Bigr) + \frac{\text{dy}}{\text{dt}}\Bigl( 2\Bigl( e_{2}e_{3} - e_{0}e_{1} \Bigr) \Bigr) + \frac{\text{dz}}{\text{dt}}\Bigl( e_{0}^{2} - e_{1}^{2} - e_{2}^{2} + e_{3}^{2} \Bigr) \Bigr){\overrightarrow{b}}_{3}\]
By contrast, the angular velocity of \(R\) with respect to the
inertial reference frame \(A\) written in terms of \(B\) is
given by:
\[_{\ }^{A}\overrightarrow{\omega} = \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3}\]
Where \(\omega_{1}\), \(\omega_{2}\) and \(\omega_{3}\) are
given in terms of Euler parameters and their derivatives:
\[\begin{split}\begin{bmatrix}
\omega_{1} \\
\omega_{2} \\
\omega_{3} \\
\end{bmatrix} = 2\begin{bmatrix}
- e_{1} & e_{0} &
e_{3} & {- e}_{2} \\
{- e}_{2} & {- e}_{3} &
e_{0} & e_{1} \\
{- e}_{3} & e_{2} &
- e_{1} & e_{0} \\
\end{bmatrix}\begin{bmatrix}
{\dot{e}}_{0} \\
{\dot{e}}_{1} \\
{\dot{e}}_{2} \\
{\dot{e}}_{3} \\
\end{bmatrix}\end{split}\]
Procedure to Derive the Equations of Motion
The equations of motion will be derived using Kane`s formalism, an
analytical mechanics method similar to Gibbs-Appell formulation. More
about Kane`s formalism can be seen in Dynamics, Theory and Applications
by Kane himself.
Kane`s equations of motion states that:
\[\boxed{\sum_{}^{}\Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} \cdot \overrightarrow{F_{i}} \Bigr) - \sum_{}^{}\Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} \cdot m_{\ }^{A}{\overrightarrow{a}}_{i} \Bigr) = 0}\]
Where the sub-index \(i\) indicates a particle in the system and
\(u_{r}\) is a generalized speed of choice. The summation is over
all particles and an equation is obtained for every \(u_{r}\) in the
system. It is simply a statement of d’Alembert Principle in terms of
generalized forces.
Kane’s equations are only valid in an inertial frame of reference and
for constant mass systems. Therefore, our strategy is to express
velocities and accelerations of every particle contained in
\(R + T + G\) with respect to reference frame \(A\). Since all
particles will be considered, we will treat the system
\(S = R + T + G\) as a constant mass one. Later, we shall use
Reynolds Transport Theorem to change from \(S\) to a control system
\(C\) which only includes the particles inside the rocket, ignoring
the ones which are leaving it.
Now, considering both the velocity of the center of mass of the rocket
and its angular velocity, we can choose generalized speeds such that:
\[_{\ }^{A}\overrightarrow{\omega} = u_{1}{\overrightarrow{b}}_{1} + u_{2}{\overrightarrow{b}}_{2} + u_{3}{\overrightarrow{b}}_{3}\]
\[_{\ }^{A}{{\overrightarrow{v}}_{*} = u_{4}{\overrightarrow{b}}_{1} + u_{5}{\overrightarrow{b}}_{2} + u_{6}{\overrightarrow{b}}_{3}}\]
This choice greatly facilitates the derivation of the equations of
motions using Kane’s method. Notice that the two equations above
represent the translational and angular velocity of \(S\) with
respect to \(A\) in terms of \(B\).
Kinematics
Notice that Kane’s equations involve two kinematic quantities, the
partial velocity and acceleration of every particle in \(S\).
Partial velocity, also known as quasi-velocity is simply
\(\partial\overrightarrow{v_{i}}/\partial u_{r}\).
The velocity of any particle in \(S\) can be written as:
\[_{\ }^{A}{\overrightarrow{v}}_{i} =_{\ }^{A}{\overrightarrow{v}}_{O} +_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - O} +_{\ }^{R}{\overrightarrow{v}}_{i}\]
Where \(_{\ }^{A}{\overrightarrow{v}}_{O}\) is the velocity of point
O which is part of rigid body R,
\(_{\ }^{A}\overrightarrow{\omega}\) is the angular velocity of
\(R\) with respect to \(A\),
\({\overrightarrow{r}}_{i - O}\) is the vector which goes from
\(O\) to \(i\). Therefore,
\(_{\ }^{A}{\overrightarrow{v}}_{O} +_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - O}\)
is the velocity of a point in the rigid body \(R\) located in the
same place as the particle \(i\) and
\(_{\ }^{R}{\overrightarrow{v}}_{i}\) is the relative velocity
between particle \(i\) and that point in \(R\). Notice how
\(_{\ }^{R}{\overrightarrow{v}}_{i}\) is always null, except for
particles in gas phase which can have arbitrary velocities.
Calculating the partial velocity of particle \(i\) with respect to a
generalized speed \(u_{r}\), we get:
\[\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{O}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - O} +_{\ }^{A}\overrightarrow{\omega} \times \frac{\partial{\overrightarrow{r}}_{i - O}}{\partial u_{r}} + \frac{\partial_{\ }^{R}{\overrightarrow{v}}_{i}}{\partial u_{r}}\]
However, the last two terms of the equation above are null. The reason
for this is because \({\overrightarrow{r}}_{i - O}\) is a vector
which is only a function of time, and its partial derivative with
respect to \(u_{r}\) is zero. The same goes for
\(_{\ }^{R}{\overrightarrow{v}}_{i}\), even though explaining why
\(_{\ }^{R}{\overrightarrow{v}}_{i}\) is only a function of time
requires a stronger assumption. We will be considering that
\(_{\ }^{R}{\overrightarrow{v}}_{i}\), which is the velocity of the
particles in the gas phase relative to the rigid body, is a known
quantity which does not depend on the partial velocities and only
changes with time.
Simplifying, we reach:
\[\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{O}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - O}\ \text{or}\ \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{O}}{\partial u_{r}} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} - \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - O}\]
Similarly, one may express the partial velocity of the center of mass
\(*\) of \(S\) as:
\[\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{O}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{* - O}\]
The derivation of this formula is the same as the derivation of an
arbitrary particle \(i\). However, one must be careful since the
center of mass is not a particle, but an abstract point. Even though it
may fall on \(R\), \(_{\ }^{R}{\overrightarrow{v}}_{*}\) is not
null, because the velocity of the center of mass depends on the
velocities of the particles in the gas phase.
Combining the last two equations, we reach:
\[\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} = \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} - \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - O} \Bigr) + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{* - O}\]
\[\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times \Bigl( {\overrightarrow{r}}_{i - O} - {\overrightarrow{r}}_{* - O} \Bigr)\]
\[\boxed{\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - *}}\]
Which gives the partial velocity of an arbitrary particle \(i\) in
\(S\) relative to the generalized speed
\(u_{r}(r = 1,\ 2,\ 3,\ 4,\ 5,\ or\ 6)\) in terms of the partial
velocity of the center of mass of \(S\), the angular velocity of
\(R\) and the vector \({\overrightarrow{r}}_{i - *}\) which goes
from the center of mass \(*\) to the particle \(i\).
The acceleration for an arbitrary particle \(i\) is given by the
time derivative of its velocity:
\[\frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial t} = \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{O}}{\partial t} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial t} \times {\overrightarrow{r}}_{i - O} +_{\ }^{A}\overrightarrow{\omega} \times \frac{\partial{\overrightarrow{r}}_{i - O}}{\partial t} + \frac{\partial_{\ }^{R}{\overrightarrow{v}}_{i}}{\partial t}\]
\[_{\ }^{A}{\overrightarrow{a}}_{i} =_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{i - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - O} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{i} +_{\ }^{R}{\overrightarrow{a}}_{i}\]
Likewise, the acceleration of the center of mass of \(S\) is given
by:
\[_{\ }^{A}{\overrightarrow{a}}_{*} =_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{*} +_{\ }^{R}{\overrightarrow{a}}_{*}\]
Combining both equations, we reach:
\[\boxed{_{\ }^{A}{\overrightarrow{a}}_{i} =_{\ }^{A}{\overrightarrow{a}}_{*} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{i - *} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times (_{\ }^{R}{\overrightarrow{v}}_{i} -_{\ }^{R}{\overrightarrow{v}}_{*}) + (_{\ }^{R}{\overrightarrow{a}}_{i} -_{\ }^{R}{\overrightarrow{a}}_{*})}\]
Application of Kane’s Equation of Motion
Using the results from the previous section on the partial velocity and
acceleration of an arbitrary particle, we can use Kane’s equation of
motion as follows:
\[\sum_{}^{}\Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} \cdot \overrightarrow{F_{i}} \Bigr) - \sum_{}^{}\Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{i}}{\partial u_{r}} \cdot m_{\ }^{A}{\overrightarrow{a}}_{i} \Bigr) = 0\]
\[\sum_{}^{}\Bigl( \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot \overrightarrow{F_{i}} \Bigr) = \sum_{}^{}\Bigl( \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot m{\overrightarrow{a}}_{i} \Bigr)\]
\[\sum_{}^{}\Bigl( \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} \Bigr) \cdot \overrightarrow{F_{i}} + \Bigl( \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot \overrightarrow{F_{i}} \Bigr) = \sum_{}^{}\Bigl( \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{r}} + \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{r}} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot m{\overrightarrow{a}}_{i} \Bigr)\]
Noting that:
\[_{\ }^{A}\overrightarrow{\omega} = u_{1}{\overrightarrow{b}}_{1} + u_{2}{\overrightarrow{b}}_{2} + u_{3}{\overrightarrow{b}}_{3}\]
\[_{\ }^{A}{{\overrightarrow{v}}_{*} = u_{4}{\overrightarrow{b}}_{1} + u_{5}{\overrightarrow{b}}_{2} + u_{6}{\overrightarrow{b}}_{3}}\]
\[\Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{1}} \Bigr) = \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{2}} \Bigr) = \Bigl( \frac{\partial_{\ }^{A}{\overrightarrow{v}}_{*}}{\partial u_{3}} \Bigr) = \Bigl( \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{1}} \Bigr) = \Bigl( \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{2}} \Bigr) = \Bigl( \frac{\partial_{\ }^{A}\overrightarrow{\omega}}{\partial u_{3}} \Bigr) = 0\]
We conclude that for \(r = 1,\ 2,\ 3,\ 4,\ 5,\ \text{or}\ 6\) we
have:
\[\sum_{}^{}{\Bigl( {\overrightarrow{b}}_{1} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot \overrightarrow{F_{i}}} = \sum_{}^{}{\Bigl( {\overrightarrow{b}}_{1} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot m{\overrightarrow{a}}_{i}}\]
\[\sum_{}^{}{\Bigl( {\overrightarrow{b}}_{2} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot \overrightarrow{F_{i}}} = \sum_{}^{}{\Bigl( {\overrightarrow{b}}_{2} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot m{\overrightarrow{a}}_{i}}\]
\[\sum_{}^{}{\Bigl( {\overrightarrow{b}}_{3} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot \overrightarrow{F_{i}}} = \sum_{}^{}{\Bigl( {\overrightarrow{b}}_{3} \times {\overrightarrow{r}}_{i - *} \Bigr) \cdot m{\overrightarrow{a}}_{i}}\]
\[\sum_{}^{}{{\overrightarrow{b}}_{1} \cdot \overrightarrow{F_{i}}} = \sum_{}^{}{{\overrightarrow{b}}_{1} \cdot m{\overrightarrow{a}}_{i}}\]
\[\sum_{}^{}{{\overrightarrow{b}}_{2} \cdot \overrightarrow{F_{i}}} = \sum_{}^{}{{\overrightarrow{b}}_{2} \cdot m{\overrightarrow{a}}_{i}}\]
\[\sum_{}^{}{{\overrightarrow{b}}_{3} \cdot \overrightarrow{F_{i}}} = \sum_{}^{}{{\overrightarrow{b}}_{3} \cdot m{\overrightarrow{a}}_{i}}\]
Rearranging and writing sum over mass in integral from assuming a
continuous system \(S\):
\[\sum_{}^{}{\Bigl( {\overrightarrow{r}}_{i - *} \times \overrightarrow{F_{i}} \Bigr) \cdot {\overrightarrow{b}}_{1}} = \iiint_{S}^{\ }{\Bigl( {\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i} \Bigr) \cdot {\overrightarrow{b}}_{1}\text{dm}}\]
\[\sum_{}^{}{\Bigl( {\overrightarrow{r}}_{i - *} \times \overrightarrow{F_{i}} \Bigr) \cdot {\overrightarrow{b}}_{2}} = \iiint_{S}^{\ }{\Bigl( {\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i} \Bigr) \cdot {\overrightarrow{b}}_{2}\text{dm}}\]
\[\sum_{}^{}{\Bigl( {\overrightarrow{r}}_{i - *} \times \overrightarrow{F_{i}} \Bigr) \cdot {\overrightarrow{b}}_{3}} = \iiint_{S}^{\ }{\Bigl( {\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i} \Bigr) \cdot {\overrightarrow{b}}_{3}\text{dm}}\]
\[\sum_{}^{}{{\overrightarrow{b}}_{1} \cdot \overrightarrow{F_{i}}} = \iiint_{S}^{\ }{{\overrightarrow{a}}_{i} \cdot {\overrightarrow{b}}_{1}\text{dm}}\]
\[\sum_{}^{}{{\overrightarrow{b}}_{2} \cdot \overrightarrow{F_{i}}} = \iiint_{S}^{\ }{{\overrightarrow{a}}_{i} \cdot {\overrightarrow{b}}_{2}\text{dm}}\]
\[\sum_{}^{}{{\overrightarrow{b}}_{3} \cdot \overrightarrow{F_{i}}} = \iiint_{S}^{\ }{{\overrightarrow{a}}_{i} \cdot {\overrightarrow{b}}_{3}\text{dm}}\]
Writing \(\sum_{}^{}\overrightarrow{F_{i}} = \overrightarrow{R}\)
and
\(\sum_{}^{}\Bigl( {\overrightarrow{r}}_{i - *} \times \overrightarrow{F_{i}} \Bigr) = {\overrightarrow{M}}_{*}\),
it simplifies to:
\[{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{1} = \iiint_{S}^{\ }{\Bigl( {\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i} \Bigr) \cdot {\overrightarrow{b}}_{1}\text{dm}}\]
\[{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{2} = \iiint_{S}^{\ }{\Bigl( {\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i} \Bigr) \cdot {\overrightarrow{b}}_{2}\text{dm}}\]
\[{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{3} = \iiint_{S}^{\ }{\Bigl( {\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i} \Bigr) \cdot {\overrightarrow{b}}_{3}\text{dm}}\]
\[\overrightarrow{R} \cdot {\overrightarrow{b}}_{1} = \iiint_{S}^{\ }{{\overrightarrow{a}}_{i} \cdot {\overrightarrow{b}}_{1}\text{dm}}\]
\[\overrightarrow{R} \cdot {\overrightarrow{b}}_{2} = \iiint_{S}^{\ }{{\overrightarrow{a}}_{i} \cdot {\overrightarrow{b}}_{2}\text{dm}}\]
\[\overrightarrow{R} \cdot {\overrightarrow{b}}_{3} = \iiint_{S}^{\ }{{\overrightarrow{a}}_{i} \cdot {\overrightarrow{b}}_{3}\text{dm}}\]
This can be written in a considerably reduced form as:
\[\overrightarrow{R} = \int_{S}^{\ }{{\overrightarrow{a}}_{i}\text{dm}}\]
\[{\overrightarrow{M}}_{*} = \int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times {\overrightarrow{a}}_{i}\text{dm}}\]
Now, expanding the integral we have:
\[\overrightarrow{R} = \int_{S}^{\ }{_{\ }^{A}{\overrightarrow{a}}_{*} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{i - *} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times (_{\ }^{R}{\overrightarrow{v}}_{i} -_{\ }^{R}{\overrightarrow{v}}_{*}) + (_{\ }^{R}{\overrightarrow{a}}_{i} -_{\ }^{R}{\overrightarrow{a}}_{*})\text{dm}}\]
\[{\overrightarrow{M}}_{*} = \int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}{\overrightarrow{a}}_{*} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{i - *} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times (_{\ }^{R}{\overrightarrow{v}}_{i} -_{\ }^{R}{\overrightarrow{v}}_{*}) + (_{\ }^{R}{\overrightarrow{a}}_{i} -_{\ }^{R}{\overrightarrow{a}}_{*}) \Bigr)\text{dm}}\]
Which after straight forward simplifications and remembering that
\(\iiint_{S}^{\ }{{\overrightarrow{r}}_{i - *}\text{dm}} = 0\)
becomes:
\[\overrightarrow{R} = m_{\ }^{A}{\overrightarrow{a}}_{*} \Rightarrow \overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{*} +_{\ }^{R}{\overrightarrow{a}}_{*} \Bigr)\]
\[{\overrightarrow{M}}_{*} = \int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times (_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{i - *}) + {\overrightarrow{r}}_{i - *} \times (_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr)) + {\overrightarrow{r}}_{i - *} \times (2_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} -_{\ }^{R}{\overrightarrow{v}}_{*} \Bigr)) + {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{R}{\overrightarrow{a}}_{i} -_{\ }^{R}{\overrightarrow{a}}_{*} \Bigr)\text{dm}}\]
Knowing that:
\[\int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{i - *} \Bigr)dm =}\int_{S}^{\ }{_{\ }^{A}\overrightarrow{\alpha} \cdot ({\overrightarrow{r}}_{i - *}^{2}\widetilde{U} - {\overrightarrow{r}}_{i - *}{\overrightarrow{r}}_{i - *})dm} = {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\alpha}\]
\[\int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times (_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr))dm =}\int_{S}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times ({\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr))dm =}_{\ }^{A}\overrightarrow{\omega}\ \times \ {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\omega}\]
Where \(\widetilde{U}\) is the identity dyadic and
\({\widetilde{I}}_{*}\) is the central inertia dyadic, we can
simplify the equations to:
\[\boxed{\overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{*} +_{\ }^{R}{\overrightarrow{a}}_{*} \Bigr)}\]
\[\boxed{{\overrightarrow{M}}_{*} = {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\alpha} +_{\ }^{A}\overrightarrow{\omega}\ \times \ {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\omega} + 2\int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times \Bigl( 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{i} \Bigr)\text{dm}} + \int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times_{\ }^{R}{\overrightarrow{a}}_{i}\text{dm}}}\]
Notice that when \(_{\ }^{R}{\overrightarrow{v}}_{i}\) and
\(_{\ }^{R}{\overrightarrow{a}}_{i}\) equals zero for every particle
in \(S\), the equations reduce to the classical equations in rigid
body dynamics.
The equations in this form are, however, not very useful since we would
have to deal with all particles in \(S\) and we would like to deal
with only particles which stay inside the rocket, that is, in a control
volume \(C\) that follows the rocket and includes all particles
inside it. Since \(S\) and \(C\) are the same in an instant, and
then become different due to particles leaving the rocket, we shall
apply Reynolds Transport Theorem to make the equations useful.
Application of Reynolds Transport Theorem
Let’s first work with the following equation:
\[\overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{*} +_{\ }^{R}{\overrightarrow{a}}_{*} \Bigr)\]
Knowing that:
\[m_{\ }^{R}{\overrightarrow{v}}_{*} = \int_{S}^{\ }{_{\ }^{R}{\overrightarrow{v}}_{i}\text{dm\ }} = \frac{d}{\text{dt}}\int_{S}^{\ }{{\overrightarrow{r}}_{i - O}\text{dm\ }}\text{and\ }\text{\ m}_{\ }^{R}{\overrightarrow{a}}_{*} = \int_{S}^{\ }{_{\ }^{R}{\overrightarrow{a}}_{i}\text{dm}} = \frac{d}{\text{dt}}\int_{S}^{\ }{_{\ }^{R}{\overrightarrow{v}}_{i}\text{dm}}\]
We can rewrite the equation as:
\[\overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times \frac{d}{\text{dt}}\int_{S}^{\ }{{\overrightarrow{r}}_{i - O}\text{dm\ }}\ \ + \frac{d}{\text{dt}}\int_{S}^{\ }{_{\ }^{R}{\overrightarrow{v}}_{i}\text{dm}}\]
Studying the integrals in \(S\) and using Reynolds Transport Theorem
when convenient, we have:
\[\frac{d}{\text{dt}}\int_{S}^{\ }{{\overrightarrow{r}}_{i - O}\text{dm\ }} = \frac{d}{\text{dt}}\int_{S}^{\ }{\rho{\overrightarrow{r}}_{i - O}\text{\ dV\ }} = \frac{d}{\text{dt}}\int_{C}^{\ }{\rho{\overrightarrow{r}}_{i - O}\text{\ dV\ }} + \int_{\partial C}^{\ }{\rho{\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
\[\frac{d}{\text{dt}}\int_{S}^{\ }{_{\ }^{R}{\overrightarrow{v}}_{i}\text{dm}} = \frac{d}{\text{dt}}\int_{S}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\text{\ dV}} = \frac{d}{\text{dt}}\int_{C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\text{\ dV}} + \int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
Where \(C\) is a control volume which includes \(R\), \(T\)
and all the particles of the gas phase \(G\) inside the rocket,
while \(S\) includes the particles which also leave the rocket. In
the instant of time \(t\) when \(S\) and \(C\) become
identical, Kane`s equations of motion are applied. Even though
calculating the integral in \(C\) or in \(S\) should give the
same result, the derivative is different. We want the derivative of the
integrals calculated in \(S\). Reynolds Transport Theorem gives us
just that, in terms of integrals done in the control volume \(C\)
and in its surface \(\partial C\).
We can rewrite the equation of motion as:
\[\overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) \Bigr) + 2_{\ }^{A}\overrightarrow{\omega} \times \Bigl( \frac{d}{\text{dt}}\int_{C}^{\ }{\rho{\overrightarrow{r}}_{i - O}\text{\ dV\ }} + \int_{\partial C}^{\ }{\rho{\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} \Bigr) + \Bigl( \frac{d}{\text{dt}}\int_{C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\text{\ dV}} + \int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} \Bigr)\]
This equation defines has no simplifications yet. However, calculating
derivatives of volume integrals concerning position and velocity of gas
phases is not exactly easy to do. However, in most rockets, the flow
which is established in the gas phase \(G\) becomes stationary rather
quickly. Therefore, this volume integrals are constant for the most part
and their derivatives become null.
With this assumption, the final form of the equation becomes:
\[\boxed{\overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) \Bigr) + 2\int_{\partial C}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times \rho{\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} + \int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}}\]
The first integral is generally referred to as the Coriolis term while
the second integral represents thrust. Notice how if
\(_{\ }^{R}{\overrightarrow{v}}_{i}\) is null in the nozzle of the
rocket, that is, if the gas does not leave the rocket, the equation
simplifies to the classical equations of rigid body motion.
Now we must deal with the following equation which describes rotational
motion:
\[{\overrightarrow{M}}_{*} = {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\alpha} +_{\ }^{A}\overrightarrow{\omega}\ \times \ {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\omega} + 2\int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times \Bigl( 2_{\ }^{A}\overrightarrow{\omega} \times_{\ }^{R}{\overrightarrow{v}}_{i} \Bigr)\text{dm}} + \int_{S}^{\ }{{\overrightarrow{r}}_{i - *} \times_{\ }^{R}{\overrightarrow{a}}_{i}\text{dm}}\]
Similarly, both integrals in this equation can be expanded using
Reynolds Transport Theorem. By using the simplification that the flow in
the gas phase is axis-symmetric, one finds that:
\[\boxed{{\overrightarrow{M}}_{*} = {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\alpha} +_{\ }^{A}\overrightarrow{\omega}\ \times \ {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\omega} + \Bigl( \frac{_{\ }^{B}d{\widetilde{I}}_{*}}{\text{dt}} \Bigr) \cdot_{\ }^{A}\overrightarrow{\omega} + \int_{\partial C}^{\ }{\rho\left\lbrack {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) \right\rbrack\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}}\]
Preparing Equations for Simulation
The following two equations have been derived to describe the complete
motion of a variable mass system:
\[\overrightarrow{R} = m\Bigl(_{\ }^{A}{\overrightarrow{a}}_{O} +_{\ }^{A}\overrightarrow{\alpha} \times {\overrightarrow{r}}_{* - O} +_{\ }^{A}\overrightarrow{\omega} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{* - O} \Bigr) \Bigr) + 2\int_{\partial C}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times \rho{\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} + \int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
\[{\overrightarrow{M}}_{*} = {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\alpha} +_{\ }^{A}\overrightarrow{\omega}\ \times \ {\widetilde{I}}_{*} \cdot_{\ }^{A}\overrightarrow{\omega} + \Bigl( \frac{_{\ }^{B}d{\widetilde{I}}_{*}}{\text{dt}} \Bigr) \cdot_{\ }^{A}\overrightarrow{\omega} + \int_{\partial C}^{\ }{\rho\left\lbrack {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) \right\rbrack\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
The only assumptions they make are that the internal gas flow is steady
and irrotational relative to the axis of symmetry.
We shall now consider how the terms in these two equations are written
in terms of the 6 degrees of freedom that we need to model in a rocket.
To do that, consider the following diagram:
A rocket without fins can be considered cylindrical in terms of
symmetry. Adding four fins to it separated by \(90{^\circ}\) does
break its cylindrical symmetry, however, important symmetry properties
remain if we choose our frame of reference \(B\) fixed to the rocket
with \({\overrightarrow{b}}_{1}\) and
\({\overrightarrow{b}}_{2}\) aligned to the fins. With this choice,
the central inertia dyadic of \(R\) is given by:
\[{\widetilde{I}}_{R} = R_{x}{\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + R_{y}{\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} + R_{z}{\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3}\]
In tensor form, considering
\({\overrightarrow{b}}_{1}{\overrightarrow{b}}_{2}{\overrightarrow{b}}_{3}\),
we have:
\[\begin{split}{\widetilde{I}}_{R} = \begin{bmatrix}
R_{x} & 0 & 0 \\
0 & R_{y} & 0 \\
0 & 0 & R_{z} \\
\end{bmatrix}\end{split}\]
Notice how \({\overrightarrow{b}}_{1}\),
\({\overrightarrow{b}}_{2}\) and \({\overrightarrow{b}}_{3}\)
are principal axis of inertia of \(R\). This is due to two facts:
first, \({\overrightarrow{b}}_{1}\) and
\({\overrightarrow{b}}_{2}\) are perpendicular to planes of
symmetries of \(R\), second, if \(R\) is rotated
\(180{^\circ}\) around \({\overrightarrow{b}}_{3}\) its mass
distribution with respect to \(B\) does not change. Consider
\(B\) position in such a way that \(R_{*}\), or the center of
mass of \(R\), is its origin to help visualize all of this.
For most rockets, we have that \(R_{x} = R_{y} = R_{I}\), therefore:
\[\begin{split}{\widetilde{I}}_{R} = R_{I}\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + R_{z}\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) = \begin{bmatrix}
R_{I} & 0 & 0 \\
0 & R_{I} & 0 \\
0 & 0 & R_{z} \\
\end{bmatrix}\end{split}\]
Similarly, the central inertial dyadic for the grains in the motor, that
is, the temporary solid phase \(T\), is given by:
\[\begin{split}{\widetilde{I}}_{T} = T_{I}\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + T_{z}\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) = \begin{bmatrix}
T_{I} & 0 & 0 \\
0 & T_{I} & 0 \\
0 & 0 & T_{z} \\
\end{bmatrix}\end{split}\]
By using the parallel axis theorem, we can write the central inertia
dyadic for \(S = R + T\):
\[{\widetilde{I}}_{*} = \Bigl( R_{I} + m_{R}a^{2} + T_{I} + m_{T}(b - a)^{2} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( R_{z} + T_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr)\]
\[\begin{split}{\widetilde{I}}_{*} = \begin{bmatrix}
R_{I} + m_{R}a^{2} + T_{I} + m_{T}(b - a)^{2} & 0 & 0 \\
0 & R_{I} + m_{R}a^{2} + T_{I} + m_{T}(b - a)^{2} & 0 \\
0 & 0 & R_{z} + T_{z} \\
\end{bmatrix}\end{split}\]
Where \(m_{R}\) is the mass of the rocket’s permanent rigid body
\(R\), \(m_{T}\) is the mass of the rocket’s propellant grains,
\(a\) is the distance between the center of mass of \(R\) and
the center of mass of \(S = T + R\), and \(b\) is the distance
between the center of mass of \(T\) and the center of mass of
\(R\).
By the center of mass relationship, we have that:
\[\Bigl( m_{R} + m_{T} \Bigr)a = m_{T}b\]
\[a = b\frac{m_{T}}{m_{R} + m_{T}}\]
Which lets us rewrite the central inertia dyadic of \(S\) as:
\[{\widetilde{I}}_{*} = \Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( R_{z} + T_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr)\]
Where \(\mu\) is the reduced mass of the system, given by:
\[\mu = \frac{m_{T}m_{R}}{m_{T} + m_{R}}\]
Like this, the central inertial dyadic of \(S = R + T\) is written
in terms of 6 parameters. Three of them are constant, \(R_{I}\),
\(R_{z}\) and \(b\), while the remaining three, \(T_{I}\),
\(T_{z}\) and \(\mu\) are functions of time.
Notice how \(S\) used to be defined as \(R + T + G\) when we
were deriving the equations of motion for the rocket and here, when
calculating the central inertia dyadic, it is defined as \(R + T\).
The reason is simple, we consider that the mass \(G\) has is
negligible when compared to \(R\) and \(T\).
Just as important as \({\widetilde{I}}_{*}\) is its derivative.
Since the derivative is taken with respect to frame of reference
\(B\), we have:
\[\frac{_{\ }^{B}d{\widetilde{I}}_{*}}{\text{dt}} = \frac{_{\ }^{B}d}{\text{dt}}\Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( R_{z} + T_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr)\]
\[\frac{_{\ }^{B}d{\widetilde{I}}_{*}}{\text{dt}} = \Bigl( \frac{dT_{I}}{\text{dt}} + b^{2}\frac{\text{dμ}}{\text{dt}} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( \frac{dT_{z}}{\text{dt}} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr)\]
The derivative of the reduced mass is given by:
\[\frac{\text{dμ}}{\text{dt}} = \frac{d}{\text{dt}}\Bigl( \frac{m_{T}m_{R}}{m_{T} + m_{R}} \Bigr) = {\dot{m}}_{T}\Bigl( \frac{m_{R}}{m_{R} + m_{T}} \Bigr)^{2} = \frac{{\dot{m}}_{T}}{m_{T}^{2}}\mu^{2}\]
Which lets us rewrite the previous equation as:
\[\frac{_{\ }^{B}d{\widetilde{I}}_{*}}{\text{dt}} = \Bigl( {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( {\dot{T}}_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr)\]
The temporal derivative of \(T_{I}\) and \(T_{Z}\) are function
of the grain geometry, and will be left as it is.
Next, we are left with three integrals which must be calculated. Since
we are dealing with the second equation, which describes only rotational
motion, lets first calculate the following:
\[\int_{\partial C}^{\ }{\rho\left\lbrack {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) \right\rbrack\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
First, remember that \(\partial C\), or the boundary of the control
volume in a rocket in which \(_{\ }^{R}{\overrightarrow{v}}_{i}\) is
not null is just its nozzle exit area and that, by flow and mass rate
conversion:
\[\rho\int_{\partial C}^{\ }{\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = \rho\pi r_{n}^{2}v_{e} = - {\dot{m}}_{T}\]
Where \(v_{e}\) is the propellant gas velocity relative to
\(R\), which is assumed uniform in the nozzle exit area, and
\(r_{n}\), as shown in the diagram is the nozzle exit radius.
With this information, one should be able to expand the integral
\(\int_{\partial C}^{\ }{\rho\left\lbrack {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) \right\rbrack\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\)
and solve it to get the following result:
\[\int_{\partial C}^{\ }{\rho\left\lbrack {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) \right\rbrack\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = \rho\pi r_{n}^{2}v_{e}\left\lbrack \Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( \frac{r_{n}^{2}}{2} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) \right\rbrack \cdot_{\ }^{A}\overrightarrow{\omega}\]
\[\int_{\partial C}^{\ }{\rho\left\lbrack {\overrightarrow{r}}_{i - *} \times \Bigl(_{\ }^{A}\overrightarrow{\omega} \times {\overrightarrow{r}}_{i - *} \Bigr) \right\rbrack\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = - {\dot{m}}_{T}\left\lbrack \Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( \frac{r_{n}^{2}}{2} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) \right\rbrack \cdot_{\ }^{A}\overrightarrow{\omega}\]
Where \(c\) is the distance between the center of mass of \(R\)
and the nozzle exit plane, and \(c - b\frac{\mu}{m_{r}} = c - a\) is
the distance between the center of mass of \(S\) and the nozzle exit
plane.
The other two integrals present in the other equation are easier to
calculate. Let’s start with the following:
\[\int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
This term can be easily recognized as the thrust produced by the motor
as exhaust is released at a high velocity \(v_{e}\) relative to the
rocket. In most cases, the thrust curve that a motor produces can be
determined by static firings in ground. With the results from the test,
other information can be approximated, such as \(v_{e}\) and
\({\dot{m}}_{T}\). We shall illustrate that next. The thrust
\(T_{\text{thrust}}\) a rocket motor produces is given by the thrust
equation from fluid dynamics (comes from the principle of momentum
conservation):
\[{\overrightarrow{T}}_{\text{thrust}} = \int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} + \pi r_{n}^{2}\Bigl( p_{e} - p_{\text{amb}} \Bigr)\overrightarrow{n}\]
Where \(p_{e}\) is the static pressure at nozzle exit and
\(p_{\text{amb}}\) is the outside ambient pressure. It is a general
assumption that, since \(p_{e} \approx p_{\text{amb}}\):
\[\int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} \gg \pi r_{n}^{2}\Bigl( p_{e} - p_{\text{amb}} \Bigr)\overrightarrow{n}\]
Which simplifies the thrust equation to:
\[{\overrightarrow{T}}_{\text{thrust}} = \int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}}\]
Once again, assuming a constant exhaust exit velocity relative to the
rocket \(v_{e}\) and that this velocity is in the direction normal
to the exit plane of the nozzle, we have:
\[{\overrightarrow{T}}_{\text{thrust}} = \text{ρπ}r_{n}^{2}v_{e} \cdot v_{e}\overrightarrow{n}\]
Remembering that we have the relation:
\[\text{ρπ}r_{n}^{2}v_{e} = - {\dot{m}}_{T}\]
We conclude that:
\[{\overrightarrow{T}}_{\text{thrust}} = - {\dot{m}}_{T}v_{e}\overrightarrow{n}\]
However, using the continuity relation we can express the thrust only as
a function of \(v_{e}\) and only as a function of
\({\dot{m}}_{T}\):
\[\left| {\overrightarrow{T}}_{\text{thrust}} \right| = \text{ρπ}r_{n}^{2}v_{e}^{2} \Rightarrow v_{e} = \sqrt{\frac{\left| {\overrightarrow{T}}_{\text{thrust}} \right|}{\text{ρπ}r_{n}^{2}}}\]
\[\left| {\overrightarrow{T}}_{\text{thrust}} \right| = \frac{{\dot{m}}_{T}^{2}}{\text{ρπ}r_{n}^{2}} \Rightarrow {\dot{m}}_{T} = - \sqrt{\left| {\overrightarrow{T}}_{\text{thrust}} \right|\text{ρπ}r_{n}^{2}}\]
Therefore, by having the thrust curve from a static firing, which gives
values of \(\left| {\overrightarrow{T}}_{\text{thrust}} \right|\) as
a function of time, we can have \({\dot{m}}_{T}\) and \(v_{e}\)
as a function of time. However, keep in mind that one must have
\(\rho\) as a function of time as well.
Not only did we derive relations to extract information about
\({\dot{m}}_{T}\) and \(v_{e}\) from
\({\overrightarrow{T}}_{\text{thrust}}\), we also derived the
relationship to express our integral:
\[\int_{\partial C}^{\ }{\rho_{\ }^{R}{\overrightarrow{v}}_{i}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = - \left| {\overrightarrow{T}}_{\text{thrust}} \right|\overrightarrow{n} = \left| {\overrightarrow{T}}_{\text{thrust}} \right|{\overrightarrow{b}}_{3}\]
The last integral which must be solved has to do the with the Coriolis
effect. If we chose \(O\) as the center of mass of the permanent
rigid body \(R\), we have:
\[2 \cdot \int_{\partial C}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times \rho \cdot {\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = 2\cdot \rho \cdot v_{e_{\ }^{A}}\cdot \overrightarrow{\omega} \times \int_{\partial C}^{\ }{{\overrightarrow{r}}_{i - R_{*}}\text{dS}}\]
\[2 \cdot \int_{\partial C}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times \rho \cdot {\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = 2\cdot \rho \cdot v_{e_{\ }^{A}}\cdot \overrightarrow{\omega} \times \Bigl( - c\pi r_{n}^{2}{\overrightarrow{b}}_{3} \Bigr)\]
\[2 \cdot \int_{\partial C}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times \rho \cdot {\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = 2\cdot c\cdot {\dot{m}}_{{T_{\ }}^{A}} \cdot \overrightarrow{\omega} \times \Bigl( {\overrightarrow{b}}_{3} \Bigr)\]
\[2 \cdot \int_{\partial C}^{\ }{_{\ }^{A}\overrightarrow{\omega} \times \rho \cdot {\overrightarrow{r}}_{i - O}\ \Bigl(_{\ }^{R}{\overrightarrow{v}}_{i} \cdot \overrightarrow{n} \Bigr)\text{dS}} = 2\cdot c\cdot {\dot{m}}_{T} \cdot(\omega_{2}{\overrightarrow{b}}_{1} - \omega_{1}{\overrightarrow{b}}_{2})\]
Having solved these integrals, we are ready to write our equations in
final form. The equation governing rotational motion is written as:
\[{\overrightarrow{M}}_{*} = \left\lbrack \Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( R_{z} + T_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) \right\rbrack \cdot \Bigl( \alpha_{1}{\overrightarrow{b}}_{1} + \alpha_{2}{\overrightarrow{b}}_{2} + \alpha_{3}{\overrightarrow{b}}_{3} \Bigr) + \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr)\ \times \ \left\lbrack \Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( R_{z} + T_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) \right\rbrack \cdot \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr)\ + \left\lbrack \Bigl( {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( {\dot{T}}_{z} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) \right\rbrack \cdot \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr) - {\dot{m}}_{T}\left\lbrack \Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr)\Bigl( {\overrightarrow{b}}_{1}{\overrightarrow{b}}_{1} + {\overrightarrow{b}}_{2}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( \frac{r_{n}^{2}}{2} \Bigr)\Bigl( {\overrightarrow{b}}_{3}{\overrightarrow{b}}_{3} \Bigr) \right\rbrack \cdot \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr)\]
We can also try a matrix representation:
\[\begin{split}\begin{bmatrix}
{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{1} \\
{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{2} \\
{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{3} \\
\end{bmatrix} = \begin{bmatrix}
R_{I} + T_{I} + b^{2}\mu & 0 & 0 \\
0 & R_{I} + T_{I} + b^{2}\mu & 0 \\
0 & 0 & R_{z} + T_{z} \\
\end{bmatrix}\begin{bmatrix}
\alpha_{1} \\
\alpha_{2} \\
\alpha_{3} \\
\end{bmatrix} + \begin{bmatrix}
\omega_{1} \\
\omega_{2} \\
\omega_{3} \\
\end{bmatrix} \times \begin{bmatrix}
R_{I} + T_{I} + b^{2}\mu & 0 & 0 \\
0 & R_{I} + T_{I} + b^{2}\mu & 0 \\
0 & 0 & R_{z} + T_{z} \\
\end{bmatrix}\begin{bmatrix}
\omega_{1} \\
\omega_{2} \\
\omega_{3} \\
\end{bmatrix} + \begin{bmatrix}
{\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} & 0 & 0 \\
0 & {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} & 0 \\
0 & 0 & {\dot{T}}_{z} \\
\end{bmatrix}\begin{bmatrix}
\omega_{1} \\
\omega_{2} \\
\omega_{3} \\
\end{bmatrix} - {\dot{m}}_{T}\begin{bmatrix}
\frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} & 0 & 0 \\
0 & \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} & 0 \\
0 & 0 & \frac{r_{n}^{2}}{2} \\
\end{bmatrix}\begin{bmatrix}
\omega_{1} \\
\omega_{2} \\
\omega_{3} \\
\end{bmatrix}\end{split}\]
However, the best representation for our purposes is in terms of its
components:
\[\begin{split}\begin{bmatrix}
M_{1} = \alpha_{1}\Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr) + \omega_{2}\omega_{3}\Bigl( {R_{Z} + T_{Z} - R}_{I} - T_{I} - b^{2}\mu \Bigr) + \omega_{1}\left\lbrack \Bigl( {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} \Bigr) - {\dot{m}}_{T}\Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr) \right\rbrack \\
M_{2} = \alpha_{2}\Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr) + \omega_{1}\omega_{3}\Bigl( R_{I} + T_{I} + b^{2}\mu - R_{z} - T_{z} \Bigr) + \omega_{2}\left\lbrack \Bigl( {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} \Bigr) - {\dot{m}}_{T}\Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr) \right\rbrack \\
M_{3} = \alpha_{3}\Bigl( R_{z} + T_{z} \Bigr) + \omega_{3}\left\lbrack \Bigl( {\dot{T}}_{z} \Bigr) - {\dot{m}}_{T}\Bigl( \frac{r_{n}^{2}}{2} \Bigr) \right\rbrack \\
\end{bmatrix}\end{split}\]
Where \(M_{1}\), \(M_{2}\) and \(M_{3}\) are given by:
\[\begin{split}\begin{bmatrix}
M_{1} \\
M_{2} \\
M_{3} \\
\end{bmatrix} = \begin{bmatrix}
{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{1} \\
{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{2} \\
{\overrightarrow{M}}_{*} \cdot {\overrightarrow{b}}_{3} \\
\end{bmatrix}\end{split}\]
The equation governing translational motion, remembering \(O\) was
chosen as the center of mass of the permanent rigid body \(R\), is
written as:
\[\overrightarrow{R} = m\Bigl( \Bigl( a_{1}{\overrightarrow{b}}_{1} + a_{2}{\overrightarrow{b}}_{2} + a_{3}{\overrightarrow{b}}_{3} \Bigr) + \Bigl( \alpha_{1}{\overrightarrow{b}}_{1} + \alpha_{2}{\overrightarrow{b}}_{2} + \alpha_{3}{\overrightarrow{b}}_{3} \Bigr) \times \Bigl( - b\frac{m_{T}}{m_{R} + m_{T}}{\overrightarrow{b}}_{1} \Bigr) + \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr) \times \Bigl( \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr) \times \Bigl( - b\frac{m_{T}}{m_{R} + m_{T}}{\overrightarrow{b}}_{1} \Bigr) \Bigr) \Bigr) + 2c{\dot{m}}_{T}\Bigl( \omega_{2}{\overrightarrow{b}}_{1} - \omega_{1}{\overrightarrow{b}}_{2} \Bigr) + \left| {\overrightarrow{T}}_{\text{thrust}} \right|{\overrightarrow{b}}_{3}\]
\[\overrightarrow{R} = m\Bigl( \Bigl( a_{1}{\overrightarrow{b}}_{1} + a_{2}{\overrightarrow{b}}_{2} + a_{3}{\overrightarrow{b}}_{3} \Bigr) + \Bigl( - b\frac{m_{T}}{m_{R} + m_{T}} \Bigr)\Bigl( - \alpha_{2}{\overrightarrow{b}}_{3} + \alpha_{3}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( \omega_{1}{\overrightarrow{b}}_{1} + \omega_{2}{\overrightarrow{b}}_{2} + \omega_{3}{\overrightarrow{b}}_{3} \Bigr) \times \Bigl( \Bigl( - b\frac{m_{T}}{m_{R} + m_{T}} \Bigr)\Bigl( \omega_{2}{\overrightarrow{b}}_{3} - \omega_{3}{\overrightarrow{b}}_{2} \Bigr) \Bigr) \Bigr) + 2c{\dot{m}}_{T}\Bigl( \omega_{2}{\overrightarrow{b}}_{1} - \omega_{1}{\overrightarrow{b}}_{2} \Bigr) + \left| {\overrightarrow{T}}_{\text{thrust}} \right|{\overrightarrow{b}}_{3}\]
\[\overrightarrow{R} = m\Bigl( \Bigl( a_{1}{\overrightarrow{b}}_{1} + a_{2}{\overrightarrow{b}}_{2} + a_{3}{\overrightarrow{b}}_{3} \Bigr) + \Bigl( - b\frac{m_{T}}{m_{R} + m_{T}} \Bigr)\Bigl( - \alpha_{2}{\overrightarrow{b}}_{3} + \alpha_{3}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( - b\frac{m_{T}}{m_{R} + m_{T}} \Bigr)\Bigl( \Bigl( - \omega_{2}^{2} - \omega_{3}^{2} \Bigr){\overrightarrow{b}}_{1} + \omega_{1}\omega_{2}{\overrightarrow{b}}_{2} + \omega_{1}\omega_{3}{\overrightarrow{b}}_{3} \Bigr) \Bigr) + 2c{\dot{m}}_{T}\Bigl( \omega_{2}{\overrightarrow{b}}_{1} - \omega_{1}{\overrightarrow{b}}_{2} \Bigr) + \left| {\overrightarrow{T}}_{\text{thrust}} \right|{\overrightarrow{b}}_{3}\]
\[\overrightarrow{R} = m\Bigl( a_{1}{\overrightarrow{b}}_{1} + a_{2}{\overrightarrow{b}}_{2} + a_{3}{\overrightarrow{b}}_{3} \Bigr) + \Bigl( - bm_{T} \Bigr)\Bigl( - \alpha_{2}{\overrightarrow{b}}_{3} + \alpha_{3}{\overrightarrow{b}}_{2} \Bigr) + \Bigl( - bm_{T} \Bigr)\Bigl( \Bigl( - \omega_{2}^{2} - \omega_{3}^{2} \Bigr){\overrightarrow{b}}_{1} + \omega_{1}\omega_{2}{\overrightarrow{b}}_{2} + \omega_{1}\omega_{3}{\overrightarrow{b}}_{3} \Bigr) + 2c{\dot{m}}_{T}\Bigl( \omega_{2}{\overrightarrow{b}}_{1} - \omega_{1}{\overrightarrow{b}}_{2} \Bigr) + \left| {\overrightarrow{T}}_{\text{thrust}} \right|{\overrightarrow{b}}_{3}\]
Or, in terms of its components:
\[\begin{split}\begin{bmatrix}
R_{1} \\
R_{2} \\
R_{3} \\
\end{bmatrix} = \begin{bmatrix}
ma_{1} + \Bigl( \omega_{2}^{2} + \omega_{3}^{2} \Bigr)bm_{T} + 2c{\dot{m}}_{T}\omega_{2} \\
ma_{2} - bm_{T}\Bigl( \alpha_{3} + \omega_{1}\omega_{2} \Bigr) - 2c{\dot{m}}_{T}\omega_{1} \\
ma_{3} + bm_{T}\Bigl( \alpha_{2} - \omega_{1}\omega_{3} \Bigr) + \left| {\overrightarrow{T}}_{\text{thrust}} \right| \\
\end{bmatrix}\end{split}\]
State-Space Representation
In order to apply numerical integrators to our equations, it is
convenient to write them in terms of state-space.
Consider the following variables: \(x\), \(y\), \(z\),
\(v_{x}\), \(v_{y}\), \(v_{z}\), \(e_{0}\),
\(e_{1}\), \(e_{2}\), \(e_{3}\), \(\omega_{1}\),
\(\omega_{2}\) and \(\omega_{3}\). Here, \((x,\ y,z)\)
represents the location of point \(O\), chosen as \(R_{*}\), in
the inertial coordinate system \(A\). On the other hand,
\(\Bigl( v_{x},\ v_{y},\ v_{z} \Bigr)\) represents the velocity of
point \(O\) relative to \(A\) written in terms of \(A\). As
expected, \(\Bigl( e_{0},\ e_{1},\ e_{2},\ e_{3} \Bigr)\) are the
quaternions used to represent the rocket’s attitude. Finally,
\(\Bigl( \omega_{1},\omega_{2},\omega_{3} \Bigr)\) represents the
rockets angular velocity relative to \(A\) but written in terms of
\(B\).
Consider the following vector:
\[\begin{split}\overrightarrow{u} = \begin{bmatrix}
x \\
y \\
z \\
v_{x} \\
v_{y} \\
v_{z} \\
e_{0} \\
e_{1} \\
e_{2} \\
e_{3} \\
\omega_{1} \\
\omega_{2} \\
\omega_{3} \\
\end{bmatrix}\end{split}\]
Defining that:
\[\begin{split}\begin{bmatrix}
a_{x} \\
a_{y} \\
a_{z} \\
\end{bmatrix} = \begin{bmatrix}
e_{0}^{2} + e_{1}^{2} - e_{2}^{2} - e_{3}^{2} & 2\Bigl( e_{1}e_{2} - e_{0}e_{3} \Bigr) & 2(e_{1}e_{3} + e_{0}e_{2}) \\
2\Bigl( e_{1}e_{2} + e_{0}e_{3} \Bigr) & e_{0}^{2} - e_{1}^{2} + e_{2}^{2} - e_{3}^{2} & 2\Bigl( e_{2}e_{3} - e_{0}e_{1} \Bigr) \\
2(e_{1}e_{3} - e_{0}e_{2}) & 2\Bigl( e_{2}e_{3} + e_{0}e_{1} \Bigr) & e_{0}^{2} - e_{1}^{2} - e_{2}^{2} + e_{3}^{2} \\
\end{bmatrix}\begin{bmatrix}
\frac{R_{1} - bm_{T}\Bigl( \omega_{2}^{2} + \omega_{3}^{2} \Bigr) - 2c{\dot{m}}_{T}\omega_{2}}{m} \\
\frac{R_{2} + bm_{T}\Bigl( \alpha_{3} + \omega_{1}\omega_{2} \Bigr) + 2c{\dot{m}}_{T}\omega_{1}}{m} \\
\frac{R_{3} - bm_{T}\Bigl( \alpha_{2} - \omega_{1}\omega_{3} \Bigr) + \left| {\overrightarrow{T}}_{\text{thrust}} \right|}{m} \\
\end{bmatrix}\end{split}\]
We have that the derivative of \(\overrightarrow{u}\) in time is
given by:
\[\begin{split}\frac{d}{\text{dt}}\overrightarrow{u} = \begin{bmatrix}
v_{x} \\
v_{y} \\
v_{z} \\
a_{x} \\
a_{y} \\
a_{z} \\
{\dot{e}}_{0} \\
{\dot{e}}_{1} \\
{\dot{e}}_{2} \\
{\dot{e}}_{3} \\
\alpha_{1} \\
\alpha_{2} \\
\alpha_{3} \\
\end{bmatrix} = \begin{bmatrix}
v_{x} \\
v_{y} \\
v_{z} \\
a_{x} \\
a_{y} \\
a_{z} \\
\frac{1}{2}\Bigl( - \omega_{1}e_{1} - \omega_{2}e_{2} - \omega_{3}e_{3} \Bigr) \\
\frac{1}{2}\Bigl( \omega_{1}e_{0} + \omega_{3}e_{2} - \omega_{2}e_{3} \Bigr) \\
\frac{1}{2}\Bigl( \omega_{2}e_{0} - \omega_{3}e_{1} + \omega_{1}e_{3} \Bigr) \\
\frac{1}{2}\Bigl( \omega_{3}e_{0} + \omega_{2}e_{1} - \omega_{1}e_{2} \Bigr) \\
\frac{\Bigl( M_{1} - \Bigl( \omega_{2}\omega_{3}\Bigl( {R_{Z} + T_{Z} - R}_{I} - T_{I} - b^{2}\mu \Bigr) + \omega_{1}\left\lbrack \Bigl( {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} \Bigr) - {\dot{m}}_{T}\Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr) \right\rbrack \Bigr) \Bigr)}{\Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr)} \\
\frac{\Bigl( M_{2} - \Bigl( \omega_{1}\omega_{3}\Bigl( R_{I} + T_{I} + b^{2}\mu - R_{z} - T_{z} \Bigr) + \omega_{2}\left\lbrack \Bigl( {\dot{T}}_{I} + {\dot{m}}_{T}b^{2}\frac{\mu^{2}}{m_{T}^{2}} \Bigr) - {\dot{m}}_{T}\Bigl( \frac{r_{n}^{2}}{4} + \Bigl( c - b\frac{\mu}{m_{R}} \Bigr)^{2} \Bigr) \right\rbrack \Bigr) \Bigr)}{\Bigl( R_{I} + T_{I} + b^{2}\mu \Bigr)} \\
\frac{\Bigl( M_{3} - \omega_{3}\left\lbrack \Bigl( {\dot{T}}_{z} \Bigr) - {\dot{m}}_{T}\Bigl( \frac{r_{n}^{2}}{2} \Bigr) \right\rbrack \Bigr)}{\Bigl( R_{z} + T_{z} \Bigr)} \\
\end{bmatrix}\end{split}\]