15.2.3 Integration of Particle Equation of Motion

The trajectory equations, and any auxiliary equations describing heat or mass transfer to/from the particle, are solved by stepwise integration over discrete time steps. Integration of time in Equation  15.2-1 yields the velocity of the particle at each point along the trajectory, with the trajectory itself predicted by

Note that Equation  15.2-1 and Equation  15.2-50 are a set of coupled ordinary differential equations, and Equation  15.2-1 can be cast into the following general form

where the term $a$ includes accelerations due to all other forces except drag force.

This set can be solved for constant $u$, $a$ and $\tau_p$ by analytical integration. For the particle velocity at the new location $u_p^{n+1}$ we get

The new location $x_p^{n+1}$ can be computed from a similar relationship.

In these equations $u_p^{n}$ and $u^{n}$ represent particle velocities and fluid velocities at the old location. Equations  15.2-52 and 15.2-53 are applied when using the analytic discretization scheme.

The set of Equation  15.2-1 and Equation  15.2-50 can also be solved using numerical discretization schemes. When applying the Euler implicit discretization to Equation  15.2-51 we get

When applying a trapezoidal discretization to Equation  15.2-51 the variables $u_p$ and $u_n$ on the right hand side are taken as averages, while accelerations, $a$, due to other forces are held constant. We get

The averages $u_p^*$ and $u^*$ are computed from

$$u_p^{*} = \frac{1}{2} (u_p^n + u_p^{n+1})$$ (15.2-56)

$$u^{*} = \frac{1}{2} (u^n + u^{n+1})$$ (15.2-57)

$$u^{n+1} = u^n + \Delta t u_p^n \cdot \nabla u^n$$ (15.2-58)

The particle velocity at the new location $n+1$ is computed by

For the implicit and the trapezoidal schemes the new particle location is always computed by a trapezoidal discretization of Equation  15.2-50.

Equations  15.2-51 and 15.2-50 can also be computed using a Runge-Kutta scheme which was published by Cash and Karp [ 48]. The ordinary differential equations can be considered as vectors, where the left hand side is the derivative $\vec{y}'$ and the right hand side is an arbitrary function $\vec{f}(t, \vec{y})$.

We get

with

$$\vec{k}_1 = \Delta t \vec{f}(t, \vec{y}^n)$$

$$\vec{k}_2 = \Delta t \vec{f}(t+a_2 \Delta t, \vec{y}^n + b_{21}\vec{k}_1)$$

$$\vec{k}_3 = \Delta t \vec{f}(t+a_3 \Delta t, \vec{y}^n + b_{31}\vec{k}_1 + b_{32}\vec{k}_2 )$$

$$\vec{k}_4 = \Delta t \vec{f}(t+a_4 \Delta t, \vec{y}^n + b_{41}\vec{k}_1 + b_{42}\vec{k}_2 + b_{43}\vec{k}_3 )$$

$$\vec{k}_5 = \Delta t \vec{f}(t+a_5 \Delta t, \vec{y}^n + b_{51}\vec{k}_1 + b_{52}\vec{k}_2 + b_{53}\vec{k}_3 + b_{54}\vec{k}_4 )$$

$$\vec{k}_6 = \Delta t \vec{f}(t+a_6 \Delta t, \vec{y}^n + b_{61}\vec{k}_1 + b_{62}\vec{k}_2 + b_{63}\vec{k}_3 + b_{64}\vec{k}_4 + b_{65}\vec{k}_5 )$$

The coefficients $a_2 \ldots a_6$, $b_{21} \ldots b_{65}$, and $c_1 \ldots c_6$ are taken from Cash and Karp [ 48]

This scheme provides an embedded error control, which is switched off, when no Accuracy Control is enabled.

For rotating reference frames, the integration is carried out in the rotating frame with the extra terms described in Equation  15.2-6 and Equation  15.2-7, thus accounting for system rotation. Using the mechanisms available for accuracy control, the trajectory integration will be done accurately in time.

The analytic scheme is very efficient. It can become inaccurate for large steps and in situations where the particles are not in hydrodynamic equilibrium with the continuous flow. The numerical schemes implicit and trapezoidal, in combination with Automated Tracking Scheme Selection, consider most of the changes in the forces acting on the particles and are chosen as default schemes. The runge-kutta scheme is recommended of nondrag force changes along a particle integration step.

The integration step size of the higher-order schemes, trapezoidal and runge-kutta, is limited to a stable range. Therefore it is recommended to use them in combination with Automated Tracking Scheme Selection.

For the massless particle type, the particle velocity is equal to the velocity of the continuous phase, hence the solution of only the trajectory Equation  15.2-50 is required where the particle velocity $u_p = u$. The new particle location along the trajectory is always computed by Equations  15.2-58 and 15.2-60, with $u_p = u$.