The integration of the Equations (27-28) is performed using backward difference formulation (BDF). BDFs as implicit schemes have less stability issues [74, for adaptive time stepping schemes as described in [71].

2.3. Tagging System

\label{tagging-system}
Due to the very high number of particles and necessity to use the binary collision, it is very important to manage the data efficiently for the sake of computational cost. In the following, a system of data communication is proposed to study collision of flows laden with numerous particles. The benefit of this system is that for each single particle, it does not require to check the entire list of particles to find the potential collision parameters. The following algorithm is based on the assumption that no more than two particles’ overlapping functions occupy each grid point of the domain:
  1. Each particle is tagged with an ID number. The walls also are tagged with the value of zero.
  2. The reduced mass variable \(m_{\text{ij}}\) and a new field variable \(P_{\text{id}}\) is declared for each grid point of the domain. \(P_{\text{id}}\) or particle ID, is a two-dimensional array whose the first argument represents the related cell and the second one with a size of two stands for ID number of the particle. For the particle \(i\), \(m_{\text{ij}}\) and \(P_{\text{id}}\) are stored to the grid points occupied by the particle when \(\psi_{i}>\varepsilon\). Here, for simplicity we do not include the first dimension related to the cells.
  3. At first \(m_{\text{ij}}=\infty\) and \(P_{\text{id}}=-1\) are set in the entire domain. Then, for the grid points nearby the walls \(P_{\text{id}}\left(1\right)=0\).
  4. For each particle, the values of \(P_{\text{id}}\) in the cells with \(\psi_{i}>\varepsilon\) is set to the particle ID number and also \(m_{\text{ij}}=m_{i}\). If \(P_{\text{id}}\neq-1\) for the first index, \(P_{\text{id}}\) of the second index is changed from \(-1\) to the particle ID number. Moreover, for those cells, \(m_{\text{ij}}=\left(\frac{1}{m_{i}}+\frac{1}{m_{j}}\right)^{-1}\)is used to update the reduced mass.
  5. Particles with \(P_{\text{id}}\left(2\right)\neq-1\) are in the process of collision and they are tagged to find the collision forces and moments. For these particles \(m_{\text{ij}}\) is obtained by using the data from grid points.