CoM, Inertia and Mass

If a density function \(\rho(r, \theta, y)\) is applied throughout the cylinder boundary illustrated in the figure above the mass will be

\[ m =\int_{y=0}^{h}\int_{\theta=0}^{2\pi}\int_{r=0}^{a}\rho(r, \theta, y)rdrd\theta dy \]

(Michel van Biezen 2017a) and the centre of mass in the \(y\) axis will be

\[ C_y = \frac{\int_{y=0}^{h}\int_{\theta=0}^{2\pi}\int_{r=0}^{a}y\rho(r, \theta, y)rdrd\theta dy}{m} \]

(Michel van Biezen 2017b). The centre of mass in the \(x\) and \(z\) axes will be

\[ C_x = \frac{\int_{y=0}^{h}\int_{\theta=0}^{2\pi}\int_{r=0}^{a}rcos(\theta)\rho(r, \theta, y)rdrd\theta dy}{m} \]

\[ C_z = \frac{\int_{y=0}^{h}\int_{\theta=0}^{2\pi}\int_{r=0}^{a}rsin(\theta)\rho(r, \theta, y)rdrd\theta dy}{m} \]

((https://math.stackexchange.com/users/250220/saketh-malyala), n.d.). The inertia tensor components of the shape created within the cylinder boundary are going to be

\[ \begin{aligned} I_{xx} &= \int_m \left( y^2 + z^2 \right) \, dm , \quad I_{yy} = \int_m \left( x^2 + z^2 \right) \, dm , \\ I_{zz} &= \int_m \left( x^2 + y^2 \right) \, dm , \quad I_{xy} = \int_m xy \, dm , \\ \;\, I_{xz} &= \int_m xz \, dm, \quad \;\, I_{yz} = \int_m yz \, dm \end{aligned} \]

(formula sheet database 2019). When considering the point of rotation to be at the centre of mass the components become

\[ \begin{aligned} I_{xx} &= \int_m \left((y-C_y)^2 + (z-C_z)^2 \right) \, dm , \\\ I_{yy} &= \int_m \left( (x-C_x)^2 + (z-C_z)^2 \right) \, dm , \\ I_{zz} &= \int_m \left( (x-C_x)^2 + (y-C_y)^2 \right) \, dm , \\ I_{xy} &= \int_m (x-C_x)(y-C_y) \, dm , \\ I_{xz} &= \int_m (x-C_x)(z-C_z) \, dm, \\ I_{yz} &= \int_m (y-C_y)(z-C_z) \, dm \end{aligned} \]

instead. Where

\[\int_m = \int_{y=0}^{h}\int_{\theta=0}^{2\pi}\int_{r=0}^a\]

\[dm = \rho\left(r, \theta, y\right)r dr d\theta dy\]

\[x = rcos(\theta),\;z = rsin(\theta)\]

E.O.M

Derivation

• Transformation
  • \(\theta \hat{n}\) – axis-angle representation of orientation.
  • \(x, y, z\) – position, also centre of mass.
• Engine
  • \(||F||\) – force (thrust).
  • \(\alpha, \beta\) – gimbal angles.
• Geometry
  • \(h\) – height.
  • \(a\) – radius.
• \(||r||\) – torque lever arm magnitude.

The linear kinetic energy is

\[ E_{lin} = \frac{1}{2}m\dot{x}^2 + \frac{1}{2}m\dot{y}^2 + \frac{1}{2}m\dot{z}^2 \]

and according to Fitzpatrick the rotational kinetic energy is

\[ \begin{aligned} E_{rot} &= \frac{1}{2}\left(I_{xx}\omega_x^2+I_{yy}\omega_y^2+I_{zz}\omega_z^2+2I_{xy}\omega_x\omega_y\right. \\ &+ \left.2I_{yz}\omega_y\omega_z+2I_{xz}\omega_x\omega_z\right) \end{aligned} \]

(Richard Fitzpatrick 2011). Lastly the potential energy is simply \(E_p = mgy\). Therefore the lagrangian of the rocket is

\[\mathcal{L} =E_{lin} + E_{rot} - E_p\]

which expands into

\[\begin{aligned} \mathcal{L} &= \frac{1}{2}m\dot{x}^2 + \frac{1}{2}m\dot{y}^2 + \frac{1}{2}m\dot{z}^2 \\ &+ \frac{1}{2}I_{xx}\omega_x^2+\frac{1}{2}I_{yy}\omega_y^2+\frac{1}{2}I_{zz}\omega_z^2\\ &+ I_{xy}\omega_x\omega_y + I_{yz}\omega_y\omega_z+I_{xz}\omega_x\omega_z \\ &- mgy \end{aligned} \]

The partial derivatives of the linear terms are

\[ \begin{aligned} \frac{\partial \mathcal{L}}{\partial \dot{x}} &= m\dot{x},\; \frac{\partial \mathcal{L}}{\partial \dot{y}} = m\dot{y},\; \frac{\partial \mathcal{L}}{\partial \dot{z}} = m\dot{z}, \\ \frac{\partial \mathcal{L}}{\partial x} &= 0,\; \frac{\partial \mathcal{L}}{\partial y} = -mg,\; \frac{\partial \mathcal{L}}{\partial z} = 0 \end{aligned} \]

and the angular terms are

\[ \begin{aligned} \frac{\partial \mathcal{L}}{\partial \omega_x} &= I_{xx}\omega_x+I_{xy}\omega_y+I_{xz}\omega_z \\ \frac{\partial \mathcal{L}}{\partial \omega_y} &= I_{yy}\omega_y+I_{xy}\omega_x+I_{yz}\omega_z \\ \frac{\partial \mathcal{L}}{\partial \omega_z} &= I_{zz}\omega_z+I_{yz}\omega_y+I_{xz}\omega_x \\ \end{aligned} \]

Keep in mind that \(\omega\) is angular velocity, that is why there is no dot written out in the partial derivatives above. Having that said, the partial derivative of just the angles, which would in this case be the partial derivative with respect to each component in a axis-angle vector is zero. This is because the terms are not present in the lagrangian. Thus, for a axis-angle vector \(\theta \hat{n}\) \[ \frac{\partial \mathcal{L}}{\partial \theta n_x} = \frac{\partial \mathcal{L}}{\partial \theta n_y} = \frac{\partial \mathcal{L}}{\partial \theta n_z} = 0 \]

The equations of motion can now be written as

\[ \begin{aligned} \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{x}}\right) - \frac{\partial \mathcal{L}}{\partial x} &= \dot{m}\dot{x}+m\ddot{x} = F_x \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{y}}\right) - \frac{\partial \mathcal{L}}{\partial y} &= \dot{m}\dot{y}+m\ddot{y} + mg= F_y \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \dot{z}}\right) - \frac{\partial \mathcal{L}}{\partial z} &= \dot{m}\dot{z}+m\ddot{z} = F_z \\ \end{aligned} \]

\[ \begin{aligned} \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \omega_x}\right) - \frac{\partial \mathcal{L}}{\partial \theta n_x} &= \dot{I_{xx}}\omega_x+I_{xx}\dot{\omega_x} \\ &+\dot{I_{xy}}\omega_y + I_{xy}\dot{\omega_y} \\ &+\dot{I_{xz}}\omega_z + I_{xz}\dot{\omega_z} = \tau_x \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \omega_y}\right) - \frac{\partial \mathcal{L}}{\partial \theta n_y} &= \dot{I_{yy}}\omega_y+I_{yy}\dot{\omega_y} \\ &+ \dot{I_{xy}}\omega_x + I_{xy}\dot{\omega_x} \\ &+ \dot{I_{yz}}\omega_z + I_{yz}\dot{\omega_z} = \tau_y \\ \frac{d}{dt}\left(\frac{\partial \mathcal{L}}{\partial \omega_z}\right) - \frac{\partial \mathcal{L}}{\partial \theta n_z} &= \dot{I_{zz}}\omega_z+I_{zz}\dot{\omega_z} \\ &+ \dot{I_{yz}}\omega_y + I_{yz}\dot{\omega_y} \\ &+ \dot{I_{xz}}\omega_x + I_{xz}\dot{\omega_x} = \tau_y \end{aligned} \]

Solution

The solutions for the linear equations of motion are as following

\[ \begin{aligned} \ddot{x} &= \frac{F_x - \dot{m}\dot{x}}{m} \\ \ddot{y} &= \frac{F_y - \dot{m}\dot{y} - mg}{m} \\ \ddot{z} &= \frac{F_z - \dot{m}\dot{z}}{m} \end{aligned} \]

Similarly the angular solutions are

\[ \begin{aligned} \dot{\omega_x} = \frac{\tau_x - \dot{I_{xx}}\omega_x-\dot{I_{xy}}\omega_y - I_{xy}\dot{\omega_y} - \dot{I_{xz}}\omega_z - I_{xz}\dot{\omega_z}}{I_{xx}} \\ \dot{\omega_y} = \frac{\tau_y - \dot{I_{yy}}\omega_y - \dot{I_{xy}}\omega_x - I_{xy}\dot{\omega_x} - \dot{I_{yz}}\omega_z + I_{yz}\dot{\omega_z}}{I_{yy}} \\ \dot{\omega_z} = \frac{\tau_z - \dot{I_{zz}}\omega_z - \dot{I_{yz}}\omega_y - I_{yz}\dot{\omega_y} - \dot{I_{xz}}\omega_x + I_{xz}\dot{\omega_x} }{I_{zz}} \end{aligned} \]

Forces and Torques

The local force vector, that is to say the force vector of the engine relative to the rocket is simply expressed as

\[ F_l = R_x(\beta)R_y(\alpha)||F||\hat{y} \]

\[R_x(\theta) = \begin{bmatrix} 1 & 0 & 0 \\ 0 & cos(\theta) & -sin(\theta) \\ 0 & sin(\theta) & cos(\theta)\end{bmatrix}\] \[R_z(\theta) = \begin{bmatrix} cos(\theta) & -sin(\theta) & 0 \\ sin(\theta) & cos(\theta) & 0 \\ 0 & 0 & 1 \end{bmatrix}\]

(Rotation matrix 2019).

From discretizing the terms \(\dot{\omega_x}\), \(\dot{\omega_y}\) and \(\dot{\omega_z}\) a axis-angle representation of the orientation of the rocket is obtained. Rotating the local force vector with this rotation gives the global force vector of the engine. That is to say the forces acting on the centre of mass aligned to the global axes. According to Jorge a vector can be rotated by an axis-angle representation in the following way

\[ v'(v) = vcos(\theta)+\left(v \cdot \hat{n}\right)\hat{n}\left(1-cos(\theta)\right)+\left(\hat{n} \times v\right)sin(\theta) \]

Thus the global force vector is simply

\[ F = v'(F_l) \]

(Jorge Rodriguez 2019). When it comes to torque, the lever arm vector is defined similarly

\[ r = v'(-||r||\hat{y}) = v'(-C_y\hat{y}) \]

from which the torque is defined as

\[ \begin{aligned} \tau_x &= r_yF_z + r_zF_y \\ \tau_y &= r_xF_z + r_zF_x \\ \tau_z &= r_xF_y + r_yF_x \\ \end{aligned} \]

References