Rotational equation for variable mass in Matlab/Simulink ignores relative angular velocity of the mass flow?

Question

I am building a simulation using the Matlab/Simulink software. The software block I am using called 'Custom Variable Mass 6DOF (Euler Angles)'. This link provides the documentation of the simulation block which uses the 6 DoF.

I am trying to intuitively understand the two equations of this block, namely, the translational equations and rotational equations of motion. If you scroll down in the 'Algorithms' section on the linked webpage, you can see these two equations the translational and rotational equations given as follows.

In the translational equation, I see the term $V_{re}$. I understand that is the velocity of mass flow relative to body. But in the rotational equation, I don't see a similar corresponding term such as the 'relative angular velocity $\omega_{re}$'. I thought the mass flow has it own angular velocity and angular momentum when it is added or ejected from the body. Why does the rotational equation not include the relative angular velocity of the mass flow?

Answer 1

I think also that the equation are not correct or misspelling ?

first I write the "rocket" equations for translation and analog for rotation, then the equations of motion in body fixed system. the index r stay for rocket, the index f for fuel.

\begin{align*} &\textbf{Translation }\\ \boldsymbol{p}_t&=(m_r+dm_f)\,\boldsymbol{v}\\ \boldsymbol{p}_{t+dt}&=m_r(\boldsymbol{v}+d\boldsymbol{v})+dm_f\,(\boldsymbol{v}+d\boldsymbol{v}+\boldsymbol{v}_{\text{rel}})\\ \boldsymbol{F}&=\frac{\boldsymbol{p}_{t+dt}-\boldsymbol{p}_t}{dt}= \frac{m_r(\boldsymbol{v}+d\boldsymbol{v})+dm_f\,(\boldsymbol{v}+d\boldsymbol{v}+\boldsymbol{v}_{\text{rel}}) -(m_r+dm_f)\,\boldsymbol{v}}{dt}\\&=m_r\,\frac{d\boldsymbol{v}}{dt}+\underbrace{ {\frac{dm_f\,d\boldsymbol{v}}{dt}}}_{=0} +\frac{dm_f}{dt}\,\boldsymbol{v}_{\text{rel}}\\ &\text{with:}\\ dm_f&=-dm_r\\\\ &\Rightarrow\\ \boldsymbol{F}&=m_r\frac{d\boldsymbol{v}}{dt}-\frac{dm_r}{dt}\,\boldsymbol{v}_{\text{rel}}\\\\ &\textbf{Rotation }\\ \boldsymbol{L}_t&=(\boldsymbol{I}_r+d\boldsymbol{I}_f)\,\boldsymbol{\omega}\\ \boldsymbol{L}_{t+dt}&=\boldsymbol{I}_r(\boldsymbol{\omega}+d\boldsymbol{\omega})+d\boldsymbol{I}_f\,(\boldsymbol{\omega}+d\boldsymbol{\omega}+\boldsymbol{\omega}_{\text{rel}})\\ \boldsymbol{M}&=\frac{\boldsymbol{L}_{t+dt}-\boldsymbol{L}_t}{dt}= \frac{\boldsymbol{I}_r(\boldsymbol{\omega}+d\boldsymbol{\omega})+d\boldsymbol{I}_f\,(\boldsymbol{\omega}+d\boldsymbol{\omega}+\boldsymbol{\omega}_{\text{rel}}) -(\boldsymbol{I}_r+d\boldsymbol{I}_f)\,\boldsymbol{\omega}}{dt}\\&=\boldsymbol{I}_r\,\frac{d\boldsymbol{\omega}}{dt}+\underbrace{ {\frac{d\boldsymbol{I}_f\,d\boldsymbol{\omega}}{dt}}}_{=0} +\frac{d\boldsymbol{I}_f}{dt}\,\boldsymbol{\omega}_{\text{rel}}\\ &\text{with:}\\ d\boldsymbol{I}_f&=-d\boldsymbol{I}_r\\\\ &\Rightarrow\\ \boldsymbol{M}&=\boldsymbol{I}_r\frac{d\boldsymbol{\omega}}{dt}- \frac{d\boldsymbol{I}_f}{dt}\,\boldsymbol{\omega}_{\text{rel}}\\\\ \end{align*} \begin{align*} &\textbf{The equations of motion in body fixed coordinate system}\\\\ &\textbf{Translation}\\ \boldsymbol{F}&=m\left(\frac{d\boldsymbol{v}}{d\tau}+\boldsymbol{\omega}\times\,\boldsymbol{v}\right)- \frac{dm}{d\tau}\,\boldsymbol{v}_{\text{rel}}\\\\ &\textbf{Rotation}\\ \boldsymbol{M}&=\boldsymbol{I}\frac{d\boldsymbol{\omega}}{d\tau} +\boldsymbol{\omega}\times (\boldsymbol{I}\,\boldsymbol{\omega})- \frac{d\boldsymbol{I}}{d\tau}\,\boldsymbol{\omega}_{\text{rel}} \end{align*} \begin{align*} \\\\\\ &\textbf{Euler angle}\\ &\textbf{rotation matrix}\\ &\boldsymbol S=S_z(\psi)\,S_x(\theta)\,S_z(\phi)\\ &\Rightarrow\\ &\boldsymbol \omega= \left[ \begin {array}{ccc} 0&\cos \left( \phi \right) &\sin \left( \theta \right) \sin \left( \phi \right) \\ 0&-\sin \left( \phi \right) &\sin \left( \theta \right) \cos \left( \phi \right) \\ 1&0&\cos \left( \theta \right) \end {array} \right] \,\dot{\boldsymbol{\varphi}}\\ &\dot{\boldsymbol{\varphi}}=\left[ \begin {array}{ccc} -{\frac {\cos \left( \theta \right) \sin \left( \phi \right) }{\sin \left( \theta \right) }}&-{\frac {\cos \left( \theta \right) \cos \left( \phi \right) }{\sin \left( \theta \right) }}&1\\ \cos \left( \phi \right) &-\sin \left( \phi \right) &0\\ {\frac {\sin \left( \phi \right) }{\sin \left( \theta \right) }}&{\frac {\cos \left( \phi \right) }{\sin \left( \theta \right) }}&0\end {array} \right] \boldsymbol{\omega} \end{align*}

Answer 2

Thanks to the OP for posting a great question. The short answer is that there is no $\vec{\omega}_{rel}$ term and the approximated dynamics used in the software is a good basic model. Unlike the translational dynamics, the simple variable mass rigid body dynamics does not involve a relative angular velocity term in the case of the attitude dynamics. In the translational case, $\vec{v}_\text{rel}:=\vec{v}_E^B$, that is the translational relative velocity of the exhaust gas $E$ measured using the reference frame attached to the variable mass rigid body system $B$ (say, say a flight vehicle). However, in the attitude case a relative angular velocity $\omega_\text{rel}$ (it is unclear what it would represent) would be an ad-hoc and unnecessary extra variable which does not arise in the dynamical analysis.

---

Prior to commencing the analysis we exemplify one of the many simplifying assumptions made in synthesizing this model. The center of mass of the variable mass system is assumed to be stationary with respect to the rigid body part of it, although in a real-world system which consumes fuel to generate exhaust, this assumption may be a good approximation but does not represent the system with complete accuracy. The notation used in the following is a combination of the tools available from the books by Zipfel and Etkin.

1. The translational dynamics are easy to obtain, but as can be seen from the some of the references in the paper cited in this answer, this simple model is missing some Coriolis terms, even in the translational case. Let us delve into some of the intricacies in the translational dynamics. Let $\frac{d^I}{dt}$ and $\frac{d^B}{dt}$ denote the time derivative of vectors calculated using the reference frames attached to an inertial $I$ and the body $B$ respectively. Let $0 < m_B(t)$ denote the mass of the system at time $0 < t$ with $\dot{m}_B < 0$. Let $\vec{F}_{B}^\text{ext}(t)$ denote the external forces acting on the system at time $t$. Applying Newton's second law of motion we have $$\vec{F}_B^\text{ext}(t) = \frac{1}{dt} [m_B(t+dt)\vec{v}_B^I(t+dt) - \dot{m}_B(t) dt \cdot \frac{d^I}{dt}\vec{s}_{EI} - m_B(t)\vec{v}_B^I(t)]$$ where $\vec{v}_P^R := \frac{d^R}{dt}s_{PR}$ denotes the time derivative of the displacement vector of point $P$ calculated using the reference frame $R$, wherein $\vec{s}_{PR}$ denotes the displacement vector of the point $P$ with respect to any fixed reference point of the reference frame $R$ (or a material particle of the rigid body of $R$) and the points $B$ and $E$ denotes the center of mass of the system at an instant of time and the location of the exhaust of materials from the body of system $B$. Further, in the equation above, the incremental decrease in the mass of the body $m_B$ given by $-\dot{m}_B \; dt$ represents the exhaust mass which has separated from the body at the time instant $t+dt$, which comes into existence after the time instant $t$. The Coriolis theorem implies that $$\frac{d^I}{dt}\vec{s}_{EI} = \frac{d^I}{dt}(\vec{s}_{BI} + \vec{s}_{EB}) = \vec{v}_B^I + \frac{d^I}{dt}\vec{s}_{EB} = \vec{v}_B^I + \vec{v}_E^B + \vec{\omega}\times\vec{s}_{EB}.$$ Putting the two equations obtained together, we have $$\vec{F}_B^\text{ext} = \frac{1}{dt} {\large [} (m_B + \dot{m}_B dt)(\vec{v}_B^I + \frac{d^I}{dt}\vec{v}_B^I dt) - \dot{m}_B dt \cdot (\vec{v}_B^I + \vec{v}_E^B + \vec{\omega}\times\vec{s}_{EB}) - m_B\vec{v}_B^I {\large ]}$$ which simplifies as $$\vec{F}_B^\text{ext} = m_B (\frac{d^B}{dt} \vec{v}_B^I + \vec{\omega}^{BI}\times\vec{v}_B^I) - \dot{m}_B \vec{v}_E^B - \dot{m}_B \vec{\omega}^{BI} \times \vec{s}_{EB}.$$ Let us denote the representation of a vector $\vec{a}$ in a coordinate system attached to the reference frame $R$ by $[a]^R$. Let us denote $[\vec{v}_B^I]^B := \vec{v}_B := [u \; v \; w]^T$, $[\vec{\omega}^{BI}]^B := \vec{\omega} := [p \; q\; r]^T$, $[\vec{v}_E^B]^B := \vec{v}_\text{rel}$ and $[\vec{F}^\text{ext}_B]^B := [F_x \; F_y \; F_z]^T$. In the body coordinate system we can then represent the translational equations of motion as $$[\vec{F}^\text{ext}_B]^B = m_B (\dot{\vec{v}}_B + \vec{\omega}\times\vec{v}_B) - \dot{m}_B \vec{v}_\text{rel} - \dot{m}_B \vec{\omega} \times [\vec{s}_{EB}],$$ where we recall that $\dot{m}_B < 0$ denotes the negative of the mass flow rate (the usual convention is to assume a positive value of mass flow rate, so that, for instance in a rocket, $\dot{m}_B < 0$ means a positive mass flow rate). On distinguishing the force due to the atmospheric and exhaust surface pressures ($0 < P_a$ and $0 < P_e$ respectively) on the exhaust mass from the other (usual) external forces, the expression $[\vec{F}_B^\text{ext}]^B$ in the equation above is replaced with $[\vec{F}_B^\text{ext}]^B - (P_e - P_a)\frac{\vec{v}_\text{rel}}{|\vec{v}_\text{rel}|}$ where these forces are expressed explicitly. In most of the literature on the rocket equation, the force due to pressure is usually not stated explicitly in the prior expression and a good reference for this modeling detail is the book by Wiesel. This effect of pressure has been neglected in the software model. Further, notice that the approximation in the software implementation is obtained by neglecting the Coriolis term as $$\dot{m}_B \vec{\omega} \times \vec{s}_{EB} \approx 0,$$ where $\vec{s}_{EB}$ denotes the displacement of the location of the exhaust on the body with respect to it's center of mass. Observe that this approximation is valid if $\|\vec{\omega}^{BI} \times \vec{s}_{EB}\| = \|\vec{\omega} \times [\vec{s}_{EB}]^B\| << \|\vec{v}_\text{rel}\| = \|\vec{v}_E^B\|,$ which includes the case $\|\vec{\omega}^{BI}\| \approx 0$ (variable mass body which does not rotate appreciably with the associated the rocket equation $[\vec{F}^\text{ext}_B]^B = m_B \vec{a}_B^I - \dot{m}_B \vec{v}_\text{rel}$ following immediately from our analysis).
2. The attitude dynamics can be obtained similarly (by avoiding more a detailed modeling approach as taken in the translational case earlier) as follows. Denote the moment of inertia tensor and the angular velocity of the variable mass system $B$ represented using the body coordinate system (coordinate system attached to body $B$) as $[\hat{I_B^B}]^B := \hat{I}$ and $[\vec{\omega}^{BI}]^B := \vec{\omega} := [p \; q \; r]^T$ respectively. The Coriolis theorem then implies that $$[\frac{d^I}{dt}H^{BI}_B]^B = [\vec{M}_B]^B := [L \; M \; N]^T = [\frac{d^I}{dt}(\hat{I}\vec{\omega})]^B = [\frac{d^B}{dt}(\hat{I}\vec{\omega})]^B + \vec{\omega} \times \hat{I} \vec{\omega},$$ where $\vec{H}_B^{BI}$ and $\vec{M}_B$ denote the angular momentum of the variable mass system calculated using the reference frame $I$ with respect to the center of mass $B$ and the moment of the external forces acting on it calculated with respect to the center of mass $B$ reference point respectively. In the variable mass rigid body case the above expression yields $\dot{\hat{I}}\vec{\omega} + \hat{I} \dot{\vec{\omega}}$ from the first term on the right hand side. The obtained result is thus $$[\vec{M}_B]^B = \dot{\hat{I}}\vec{\omega} + \hat{I} \dot{\vec{\omega}} + \vec{\omega} \times \hat{I} \vec{\omega}.$$ Indeed, the presented equation of motion for the attitude dynamics can be modified by including more terms which approximate the dynamics more precisely. See section 7.11 of the book by Thomson (further section 7.8 analyzes jet damping of transverse pitch-yaw damping due to exhaust) for one for one of the most complete analyses of the rotation dynamics of a variable mass vehicle body. However, a relative angular velocity $\vec{\omega}_\text{rel}$, unless imposed ad-hoc and unnecessarily on the dynamical formulation does not appear in the mathematical representation. It is tempting to attempt to analogously write the attitude dynamics by observing the translational dynamics via ad-hoc addition of such a term. However, this would result in an inaccurate model.

The equations which are used in the software constitute an approximate model for a simplistic simulation of the six degree of freedom motion of a variable mass system. Indeed, the equations can be modified by including additional terms for higher fidelity modeling. See this paper from JPL and it's references for further reading. An example of additional detail in variable mass system dynamical modeling (in comparison to the model used in the software) which is included in the referenced work is the motion of the center of mass relative to the rigid body component of the system (caused by, say, fuel consumption which is required to generate the exhaust).