Is there a variable mass Verlete like integration algorithm

Question

I'm currently modeling the explosion of a star. For my simulation, I use a Verlete like integration algorithm. This is quite common in celestial mechanics modeling. The thing is that now that I have to account for mass change I don't know if the algorithm is still valid. The Newtons second law says
\begin{equation}
F=m\frac{d^2X}{dt^2}+V\frac{dm}{dt}
\end{equation}
The typical Verlete like algorithm is
\begin{equation}
V+=\frac{dV}{dt}\cdot dt
\end{equation}
\begin{equation}
X+=\frac{dX}{dt}\cdot dt
\end{equation}
Some of them multiply by some other constant and there are many variations. So the question is if I take into account mass variation the algorithm becomes something like this?
\begin{equation}
V+=\frac{dV}{dt}\cdot dt+\frac{V}{m}\frac{dm}{dt}\cdot dt
\end{equation}
\begin{equation}
X+=\frac{dX}{dt}\cdot dt
\end{equation}
Is this a valid approach? Or the algorithm stays the same?

Answer 1

Starting your first equation, I rewrite by shifting $dV/dt$ to the left hand side: \begin{align} F&=m\frac{d V}{dt}+V\frac{dm}{dt}\\ \frac{d V}{dt} &= \frac{F}{m} - \frac{V}{m}\frac{dm}{dt} \end{align}

Therefore,

\begin{align} V+ = \frac{d V}{dt} dt &= \frac{F}{m} dt - \frac{V}{m}\frac{dm}{dt} dt\\ &= \frac{F}{m} dt - V\frac{d (\ln m)}{dt} dt \end{align}

In the rocket motion, the second term is a positive acceleration, since the mass changing rate $\frac{dm}{dt}$ is negative.

But this is only valid for case of no explosion. That means the separated part $dm$ is of a same velocity with the main body $m(t)$. For case of internal explosion, there is a extra expelling force which is not taking into account in the $ \frac{d(m(t)v(t))}{dt}$ term. You have to refer to answer of Elis for internal explosion.

For extra expelling force or even not in the same direction, we have to start from the conservation of momentum:

Momentum before separation: $$ \vec P(t) = m(t) \vec V(t). $$

Momentum after separation: \begin{align} \vec P(t+\Delta t) &= m(t+\Delta t) \vec V(t+\Delta t) + \Delta m \left( \vec U(t)+ \vec V(t) \right)\\ &\approx m(t) \vec V(t) - \Delta m \vec V(t) + m(t) \Delta \vec{V} +\Delta m \left( \vec U(t)+ \vec V(t) \right) \\ &= m(t) \vec V(t) + m(t) \Delta \vec{V} +\Delta m \, \vec U(t) \end{align} Where $\Delta m = - \frac{dm}{dt} \Delta t$, and $\vec U$ is the velocity of the separated part w.r.t to the reference frame of the mass moving object, therefore the velocity of the separated w.r.t the ground is $\vec U(t) + \vec V(t)$. For linear rocket motion $\vec U(t)$ is in the opposite direction of $\vec V$.

For the conservation of momentum, $ P(t+\Delta t) - \vec P(t) = \vec F \Delta t$, we then have: $$ m(t)\Delta \vec V(t) + \Delta m\, \vec U(t) = \vec F \Delta t. $$ recall that $\Delta m = -\frac{dm}{dt} \Delta t$. We than have the differential form: \begin{align} & m(t) \Delta \vec V(t) -\frac{dm}{dt} \Delta t\, \vec U(t) = \vec F \Delta t.\\ & \Delta \vec V(t) =\left( \frac{\vec F}{m} + \frac{d (\ln m)}{dt} \vec U(t) \right)\Delta t \end{align}

This renders: $$ \vec V += \left( \frac{\vec F}{m} + \frac{d (\ln m)}{dt} \vec U(t) \right)\Delta t. $$

This may include the cases of explosion and side impact motion. For rocket motion. we consider $\vec U(t) = \vec U_0$ is a constant. But for mass attaching problem, it may be more appropriate to consider $\vec U(t) + \vec V(t)$ is a constant.