When dealing with systems with variable mass, one needs to consider not only the "reactive forces", but also the shifting of the center of mass of a system due to its mass redistribution. Let us find an equation of motion for the center of mass of an open and not isolated system A.
In order that Newton's Laws be applicable, we must consider both the closure and isolation of this system A + B by identifying all the particles that flow out and into the system during a time interval \Delta t as well as all particles that exert forces on the particles.
The acceleration of the center of mass of the system is given by:
<br />
\mathbf{a}^{(A)}_{\mathrm{C M}}(t) = \frac{d \mathbf{v}^{(A)}_{\mathrm{C M}}(t)}{dt}<br />
where the velocity of the center of mass is by definition given by:
<br />
\mathbf{v}^{(A)}_{\mathrm{C M}}(t) \equiv \frac{\mathbf{P}^{(A)}(t)}{M^{(A)}(t)}<br />
where:
<br />
\mathbf{P}^{(A)}(t) = \sum_{a \in A(t)}{\mathbf{p}_{a}(t)}<br />
is the total momentum of the particles in the system A at the instant t and:
<br />
M^{(A)}(t) = \sum_{a \in A(t)}{m_{a}}<br />
is their total mass. As was pointed out by D-H's derivation in post #37, these quantities can be time dependent in two different ways: Either because the quantities that enter the sum change with time, or because the set over which we sum changes itself, or both, of course.
According to the rules of Calculus, we may write:
<br />
\mathbf{a}^{(a)}_{\mathrm{C M}})(t) = \frac{1}{M^{(A)}(t)} \, \frac{d \mathbf{P}^{(A)}(t)}{d t} - \frac{\mathbf{P}^{(A)}(t)}{[M^{(A)}(t)]^{2}} \, \frac{d M^{(A)}(t)}{d t}<br />
or:
<br />
M^{(A)}(t) \, \mathbf{a}^{(A)}_{\mathrm{C M}}(t) = \frac{d \mathbf{P}^{(A)}(t)}{d t} - \frac{\mathbf{P}^{(A)}(t)}{M^{(A)}(t)} \, \frac{d M^{(A)}(t)}{d t} = \frac{d \mathbf{P}^{(A)}(t)}{d t} - \mathbf{V}^{(A)}_{\mathrm{C M}}(t) \, \frac{d M^{(A)}(t)}{d t}<br />
Let us find the derivatives on the rhs of this equality. In doing so, we will adopt the following notation. Let the index a enumerate the particles that were inside the system A at time t; let the index i enumerate the particles that exited system A into system B and let j enumerate the particles that entered from B to A during the time interval \Delta t. Then we have:
<br />
\mathbf{P}^{(A)}(t + \Delta t) = \sum_{a}{\mathbf{p}_{a}(t + \Delta t)} - \sum_{i}{\mathbf{p}_{i}(t + \Delta t)} + \sum_{j}{\mathbf{p}_{j}(t + \Delta t)}<br />
We will keep quantities up to order O(\Delta t) in the above sum. It is important to realize that the number of particles that enter or exit (and, therefore, both the mass and momentum they carry with them) is a quantity of the order O(\Delta t). That is why we can exchange the argument from t + \Delta t to t in the last two sums. Similarly, we can write:
<br />
p_{a}(t + \Delta t) = p_{a}(t) + \Delta t \, \sum_{b \in A + B, b \neq a}{\mathbf{F}_{b a}(t)} + o(\Delta t) = p_{a}(t) + \Delta t \, \left[\sum_{a' \in A, a' \neq a}{\mathbf{F}_{a' a}(t)} + \sum_{b \in B}{\mathbf{F}_{b a}(t)}\right] + o(\Delta t)<br />
Collecting everything together, we may write:
<br />
\mathbf{P}^{(A)}(t + \Delta t) = \mathbf{P}^{(A)}(t) + \Delta t \, \left[\sum_{a, a' \in A, a' \neq a}{\mathbf{F}_{a' a}(t)} + \sum_{a \in A, b \in B}{\mathbf{F}_{b a}(t)}\right] - \sum_{i \in \Delta N_{\mathrm{out}}}{m_{i} \mathbf{v}_{i}(t)} + \sum_{j \in \Delta N_{\mathrm{in}}}{m_{j} \mathbf{v}_{j}(t)} + o(\Delta t)<br />
<br />
\frac{d\mathbf{P}^{(A)}(t)}{d t} = \sum_{a, a' \in A, a' \neq a}{\mathbf{F}_{a' a}(t)} + \sum_{a \in A, b \in B}{\mathbf{F}_{b a}(t)} - \lim_{\Delta t \rightarrow 0}{\left[\frac{1}{\Delta t} \sum_{i \in \Delta N_{\mathrm{out}}}{m_{i} \, \mathbf{v}_{i}(t)} \right]} + \lim_{\Delta t \rightarrow 0}{\left[\frac{1}{\Delta t} \sum_{j \in \Delta N_{\mathrm{in}}}{m_{j} \, \mathbf{v}_{j}(t)} \right]}<br />
The first term in this equation is exactly equal to zero due to Third Newton's Law (the summation is over the same system A and \mathbf{F}_{a a'} = -\mathbf{F}_{a' a}) The second term is the sum of all the instantaneous external forces that act on the system \sum{\mathbf{F}_{\mathrm{ext}}(t)}. The third and the fourth term are the outgoing and incoming momentum flux, respectively, these are frame dependent quantities (because velocity is a frame dependent quantity). Although expressed as limits, these are actually finite quantities (if the particle flux is finite).
The other derivative is:
<br />
M^{(A)}(t + \Delta t) = \sum_{a}{m_{a}} - \sum_{i}{m_{i}} + \sum_{j}{m_{j}} = M^{(A)}(t) - \sum_{i \in \Delta N_{\mathrm{out}}}{m_{i}} + \sum_{j \in \Delta N_{\mathrm{in}}}{m_{j}}<br />
<br />
\frac{d M^{(A)}(t)}{d t} = -\lim_{\Delta t \rightarrow 0}{\left \frac{1}{\Delta t} \, \sum_{i \in \Delta N_{\mathrm{out}}}{m_{i}}\right]} + \lim_{\Delta t \rightarrow 0}{\left \frac{1}{\Delta t} \, \sum_{j \in \Delta N_{\mathrm{in}}}{m_{j}}\right]}<br />
Again, both of these terms are finite.
Combining everything together, we can finally write:
<br />
M^{(A)}(t) \, \mathbf{a}^{(A)}(t) = \sum{\mathbf{F}_{\mathrm{ext}}}(t) - \mathbf{\Pi}_{\mathrm{out}}(t) + \mathbf{\Pi}_{\mathrm{in}}(t)<br />
where \mathbf{\Pi}_{\mathrm{in/out}}(t) is the incoming (outgoing) momentum flux in the center-of-mass reference frame (a frame independent quantity, because relative velocities are Galilean invariants):
<br />
\mathbf{\Pi}_{\mathrm{out}}(t) = \lim_{\Delta t \rightarrow 0}{\left[\frac{1}{\Delta t} \, \sum_{i \in \Delta N_{\mathrm{out}}}{m_{i} \, \mathbf{u}_{i}}\right]}, \; \mathbf{u}_{i} = \mathbf{v}_{i} - \mathbf{v}^{(A)}_{\mathrm{C M}}<br />
<br />
\mathbf{\Pi}_{\mathrm{in}}(t) = \lim_{\Delta t \rightarrow 0}{\left[\frac{1}{\Delta t} \, \sum_{j \in \Delta N_{\mathrm{in}}}{m_{j} \, \mathbf{u}_{j}}\right]}, \; \mathbf{u}_{j} = \mathbf{v}_{j} - \mathbf{v}^{(A)}_{\mathrm{C M}}<br />
By specifying the number of particles that leave (enter) the system in unit time \varphi_{\mathrm{out}/\mathrm{in}}(m, \mathbf{u}, t) with a relative velocity in a unit velocity-space volume around \mathbf{u} and in unit interval around the mass m, we can express the above quantities as:
<br />
\mathbf{\Pi}_{\mathrm{out}/\mathrm{in}}(t) = \int{m \, \mathbf{u} \, \varphi_{\mathrm{out}/\mathrm{in}}(m, \mathbf{u}, t) \, d^{3}u \, dm}<br />
Also, the mass rate of change is:
<br />
\frac{d M^{(A)}(t)}{d t} = -\mu_{\mathrm{out}}(t) + \mu_{\mathrm{in}(t)}<br />
where the incoming (outgoing) mass flux rate is given by:
<br />
\mu_{\mathrm{out}/\mathrm{in}}(t) = \int{m \varphi_{\mathrm{out}/\mathrm{in}}(m, \mathbf{u}, t) \, d^{3}u \, dm}<br />
