Research on particle collision forward search algorithm and CFD–DEM variable time step coupling calculation method

In CFD–DEM coupling calculation, the calculation accuracy will decrease if the particle calculation time step is too large, while the calculation efficiency will decrease if the particle calculation time step is too small. In this study, first, the displacement of the center particle and the fastest particle in the computational domain is added to construct a search sphere. Subsequently, the particles that may collide are screened to create a search list, and particle collisions are determined using a forward search method. Finally, an improved DEA method is proposed; the time step of particle calculation is automatically adjusted according to the collision time, which solves the contradiction between selecting particle calculation time step, accuracy, and efficiency. The relative error between the particle collision simulation results and the analytical solution obtained by this method is less than 3.00%. In this study, three selected calculation time steps can ensure excellent calculation accuracy and efficiency. In multi-particle and fluid coupling simulation analysis, when the calculation time step is 10–5s, the calculation efficiency is improved by 19.8%.


Introduction
Solid-liquid two-phase flow exists in numerous chemical, petroleum, agricultural, biochemical, and food production processes [1,2]. Numerical method is widely used in the improvement in process equipment [3,4]. At present, there are two commonly utilized numerical models to explain the solid-liquid two-phase flow [5]. The one is the two-fluid model (TFM), which regards fluids and solids as fluids or continuous media that permeate each other [6]. The other is computational fluid dynamics-discrete element method (CFD-DEM), which regards particles as discrete media and uses discrete element method to solve the problem. The particle trajectory [7] is calculated using Newton's second law of motion. Flow is described using the Navier-Stokes and continuity equations [8], with the energy equation sometimes applied. The coupling interface between CFD and DEM realizes the interaction between particles and fluid, in which CFD simulates fluid flow and DEM simulates particle motion and collision [9]. The discrete particle model (CFD-DEM) is different from the two-fluid model (TFM) in that it can track each moving particle [10]. Therefore, the discrete particle model (CFD-DEM) is usually used when considering the particle collision effect.
In traditional CFD-DEM coupling calculations, the time step of the fluid calculation is set as an integral multiple of the time step of particle calculation ( t f j · t p ) to reconcile the disparity between the fluid calculation efficiency and the particle calculation accuracy [11]. The fluid as a continuous medium, a relatively large calculation time step is selected in order to ensure the efficiency of fluid calculation. The particles are set with a small time step in order to avoid the collision of missing particles and ensure the accuracy of particle calculation. A fixed time step t p is set, and the collision between any two particles is examined after solving each calculation time step. If a collision occurs, it is analyzed using either the soft-sphere or hard-sphere model. The problem with this method is that two particles may collide, overlap, and then separate again in a time step, but leaving no evidence of collision at the end of the time step. Because the particle has different running time before the collision, it is difficult to ensure that the calculation time step concludes at the same instant as the collision. Reducing the time step for particle calculation cannot resolve the above problem, and it increases the number of particle calculation iterations and extends the particle calculation time. The particle collision algorithm is the key to enhancing particle calculation accuracy and efficiency. Numerous researchers have investigated particle collision algorithms, such as the grid element method suggested by Hockney and Eastwood [12]. In this method, the analysis space is divided into regular grids, and each particle is in a certain grid. For a particle, the particles or boundaries in its grid and its adjacent grid are the neighbor units of the particle. It is sufficient to detect only the particle's unit and nearby unit when detecting particle collisions. Allen and Tildesley [13] proposed a neighbor chain list method. A ball is drawn with a particle as the center and a specific length r (r > 1.5R max , where R max is the radius of the largest particle in this system) as the radius. Then, other particles in the ball are employed as the neighboring units of the particle, and collision detection is required for them. Cohen et al. [14], Baraff [15], and Lin [16] suggested the bounding box method. The parallel to the coordinate plane and related to the particles bounding boxes are projected vertically to the x and y coordinate axes. The projection of boundary box on both coordinate axes must overlap if two particles collide.
On the basis of grid element method, neighbor chain list method, and boundary box algorithm, researchers have investigated the accuracy [17], speed [18,19], and adaptability [20] of grid optimization [21], size, and shape. Nezami et al. [22] suggested a new fast common plane (FCP) collision search algorithm for resolving the common plane between polygon particles. Munjiza and Andrews [23] presented the no-binary search (NBS) contact collision detection algorithm, in which the entire detection time is proportional to the total number of particles. Chen et al. [24] suggested an approach for optimizing direct simulations of particle collisions in three-dimensional space. Pischke et al. [25] suggested a hybrid collision algorithm combining deterministic and random collisions. Wang et al. [26] presented a two-step accurate collision detection fast algorithm for 'bounding sphere-maximum detection area' preprocessing. The above determination of particle collision is calculated by fixed time step. This study suggests a particle collision forward search algorithm to resolve the discrepancy among particle time step selection, calculation accuracy, and efficiency in CFD-DEM coupling calculation. In this algorithm, the initial particle calculation time step does not affect the particle calculation accuracy, and the particle calculation time step is determined by the particle collision time step. On the premise of accurately describing particle collisions, the algorithm realizes the calculation of large time step, and the number of fluid-particle coupling is determined by the number of particle collisions. Therefore, it will be an efficient and accurate calculation method of fluid-solid coupling.
2 Particle-fluid coupling calculation method

Fluid control equation
The position of particles is projected onto the CFD calculation unit grid. The fluid volume fraction occupied by the residual fluid in the computational element mesh α f is calculated during the coupling calculation of particles and fluid. Next, the fluid force and torque are calculated utilizing the fluid velocity u f , particle velocity u p , fluid density, ρ f , and dynamic fluid viscosity μ. The iterative process of the particle-fluid coupling calculation is shown in Fig. 1. The fluid volume fraction in the unit [27] α f is calculated as follows: where α p is the volume fraction of particles in a fluid unit, V E the fluid unit volume, and the ith particle volume in the fluid unit. S f is the fluid and particle momentum exchange force source term [28] calculated as where F (i) f is the fluid force exerted on the i particle in the fluid unit and n p the number of particles in the fluid unit.
Fluid continuity equation [29]: The conservation equation of fluid momentum is: where ρ f is the fluid density, u f is there presents fluid velocity vector, τ is the fluid viscous stress tensor. The expression τ can be defined as follows [30]: where μ f is the fluid viscosity coefficient, λ f is the fluid volume viscosity, I is the fluid viscous stress tensor.
is the pressure gradient on the surface of the microelement, ω p is the particle's rotation angular velocity, o(Re) is the lowest order of Re is the series remainder that is not explicitly written as 1, and D p is the particle diameter

Particle dynamics equation and force analysis
(1) Particle dynamics equation: Particle motion comprises slip and rotation motions, based on Newton's second law [31,32]: where m i is the quality of the particle i, u i is the translational velocity of the particle, F cn, i j and F ct, i j are normal and tangential forces of particle collision, F f is the resultant fluid force on particles, F g represents particle gravity, I i is the inertia moment of particle i, ω i is the particle's angular velocity while and M c is the torque caused by particle collisions.
The particles are inside the fluid and follow the fluid motion. The unstable flow inside the fluid produces various forces on the particles, and the resultant fluid force on the particles is F f calculated as follows: where F d , F p , F vm , F saff , F ml and F b is the resistance of particles (drag), pressure gradient force, virtual mass force, Saffman shear lift, Magnus rotation lift, and Basset force, respectively ( Table 1). As stated in previous reports [33,34], resistance is the main force. The drag force computation formula suggested by Di Felice is as follows [39]: Re where C d is the resistance coefficient, ρ p is the particle density, and Re p is the Reynolds number of the particle.

Contact force of particles
The following assumptions are made to construct a method for calculating particle collisions: (1) The particles are rigid spherical bodies. (2) Consider only two-particle collisions (3) The particle densities are identical, the particle diameters may vary, and the initial distribution of the particle group is uniform. If any two particles i and j are taken, the diameter of the particle i is set as D i p , the mass is set to m i , the diameter of the particle j is set as D i p , and the mass is set as m j . According to the momentum theorem, the following equation must be fulfilled when two particles collide: where u i and u i are translational velocities of particle i before and after collision, ω i and ω i are rotational velocities of particle i before and after collision, k is the normal unit vector of particle relative velocity, J is the impulse when particles collide, and I i and I j are inertia moments of particle i and particle j, respectively. According to the theory of elastic collision mechanics [40], the collision time t c of the two particles is where e is the coefficient of restitution, u * i j is the relative velocity of particles after correction, u * i j (1 − χ)u i j + χ u i j , u i j , u i j is the relative velocity after collision, χ is the weight coefficient, M represents the equivalent mass of particle i and particle j, is the correlation coefficient of the collision force, α 1 (1 − μ 2 i )/π E i , α 2 (1 − μ 2 j )/π E j , μ i and μ j is the Poisson ratio of particles i and j, respectively, and E i and E j is the elastic modulus of the particles i and j, respectively.
The impulse theorem provides the following formula for calculating the collision force between two particles: First, based on the relative velocity u i j of collision particles i and j, the particle slip is determined, and the velocity after the collision is calculated. Then, utilizing the interpolation algorithm, the relative velocities before and after particle collision u i j and u i j are weighted and interpolated to acquire the u * i j particle collision time t c and collision force F c,ij by substituting them into Eqs. (15) and (16). Next, they are substituted into Eqs. (6) and (7) to get the velocity after a particle collision. Eventually, the collision calculation is complete if the difference between the set and acquired values satisfies the speed convergence tolerance; otherwise, these stages are repeated until the tolerance is achieved [41].
2.4 Fluid-particle coupling calculation method The search for particle collisions is a forward search. The collision of particle before t + t p is determined according to its position and velocity at t. If this occurs, the collision time step t i c can be calculated. If there is no collision, t p t p then it can be considered for the calculation. The particle collision calculation is performed using the collision time step as the calculation time step when particles collide to prevent collisions between missing particles. When there was no collision between the particles, the particles were calculated across a large time step with increased calculation efficiency. A higher calculation time step increases the particle displacement and makes the fluid volume fraction α f and momentum exchange force source term S f change too fast. The convergence criteria can be utilized to limit this change, ensuring that the fluid-particle coupling's calculation error does not vary with time step.
In calculating fluid-particle coupling, the effect of fluid on particles is determined by substituting the fluid force acting on particles into the particle motion equation to solve the particle motion. In contrast, the effect of particles on the fluid is determined by replacing the fluid volume fraction α f and momentum exchange force source term S f into the fluid control equation. The movement of particles directly alters the fluid volume fraction α f and momentum exchange force source term S f . Thus, the convergence criterion of fluid-particle coupling is established utilizing the variation in the fluid volume fraction α f and momentum exchange force source term S f .
where ε tol f , α is the convergence tolerance of the fluid volume fraction, ε tol f , S is the convergence tolerance of the momentum exchange force source term, S 0 f , and S N f is the starting and end momentum exchange force source terms, respectively, α 0 f and α N f is the starting and ending fluid volume fractions, respectively. The fluid control equation is solved using the SIMPLE algorithm, as shown in Fig. 2. First, the fluid force of particles is calculated using the relative velocities of the particles and the fluid. The information is then transmitted to the particles, and the DEM solver is initiated. The particle motion equation is used to solve the particle velocity at the current time and evaluate the particle collision in the initial time step of the calculation. The collision time step is used as the time step to advance the particle calculation. Based on the change in the particle position and force in the computational domain, the volume fraction α f and momentum exchange force source term of the fluid S f are calculated. If the convergence criterion is unsatisfied, DEM solver is stopped, and CFD solution is initiated.
The CFD solution is shown in Fig. 2a. The CFD solution is initialized, and the fluid is calculated in one time step with the initial calculation time step ( t f ), and the fluid force (F f ) on the particles is obtained. Once the CFD solver was over, the DEM solution was started.
The DEM solution is shown in Fig. 2b. The particles set the initial calculation time step t p t f . (1) The particle collision inside the initial calculation time step was determined based on the position and velocity at the beginning of the particle calculation. If the particle does not collide, the particle's forward time step is considered equivalent to the initial calculation time step of the particle t p t p . Following the calculation, convergence was determined again. If the calculation results converge, the DEM solution is completed. If it does not converge, the time step of particle calculation is shortened to t * p t p · k n , and DEM is solved again until convergence is achieved.
(2) If the particle collides, calculating the particle collision time step t i c requires judging the particle collision time step and the residual time of the initial calculation time step (i indicates the number of particle collisions, m If it is not tunable, the residual time in the initial time step of the particle is utilized as the calculation to solve. Following the calculation, convergence was determined again. If the calculation results converge, the DEM solution is completed. If it does not converge, the time step of particle calculation is reduced to t * p t p · k n , and DEM is solved again until convergence is completed. , it was established, the particle moved forward t p t i c with the collision time step as the particle calculation time step. After the particle collision calculation is achieved, convergence judgment is carried out. (5) If it is convergent, the particle collision search is continued until it is not convergent or the initial calculation time of the particle is completed. (6) If it does not converge, the particle collision time step t * c t i c · k n is shortened, and the DEM is solved again until convergence is completed. After DEM solution was completed, the fluid volume fraction α f and momentum exchange force source item S f were calculated. The fluid calculation time step t f t * p was adaptively modified in accordance with the calculation time step t * p completed by the particles. The fluid domain was resolved to determine whether t f + t f < T it is tunable. If it is tunable α f and S f are transferred to CFD, the CFD solution is restarted. Otherwise, the coupling calculation is terminated. In particle collision search, the center particles can only collide with the particles surrounding them, and a certain range surrounding the center particles is separated into collision search areas. The method used is as follows: Assuming that the center particle i collides with the fastest moving particle j in the computational domain after the computational time step, the center particle collision search space surrounding the sphere is established as the distance R between the center particle i and the particle j as the radius. As shown in Fig. 3, it is assumed that all particles in the sphere can collide with the center particle i: where L i is the motion displacement of the central particle. L i u i · t p , L max u j · t p is the motion displacement of the fastest moving particle, D i p is the diameter of the central particle, D j p is the diameter of the fastest moving particle, u i is the motion velocity of the particle i, u j is the motion velocity of the particle j.

The establishment of particle collision list
The particles within the sphere surrounding the central particle collision search space were further screened to maintain the possible collision particles and enhance the efficacy of the collision search. The time center particle coordinates were set at (x i , y i , z i ) time t. The coordinates of particles g were set as (x g , y g , z g ), making the distance between both particles L ig During the calculation time step, if the sum of the movement displacement of the central particle i and the particle g was more significant than distance L ig between the two particles at a time t, the particle g can collide with the central particle i: where L g u g · t p is the movement displacement of particles g. Therefore, the remaining particles established a central particle collision search list after screening.

Forward search algorithm for collision particles
According to the speed, direction, acceleration, and other characteristics of particles, the particle forward collision search algorithm calculates the moment of particle collision and determines if a collision happened in the initial time step. If a collision happens inside a time step, change the particle calculation time step based on the collision time. If no collision happens, the propulsion calculation uses the original time step. As shown in Fig. 4, t is the starting time, the initial position of the particle i is (x i , y i , z i ), and the particle j is (x j , y j , z j ), making the particle relative distance force of particles i and j were F i(t) and F j(t) , the particle accelerations were a i(t) m j , the velocities were u i(t) and u j(t) , and the spatial positions were r i and r j , respectively. k shows the unit vector k r j − r i The direction considered was from the center of the particle i to the center of the particle j.
When the distance among the particles at a time t is L i j > 1 2 (D i p + D j p ), the particles do not collide. In the process of particles moving forward, the velocity and acceleration produced components in the unit vector k, and the two particles gradually moved closer. At the time t + t c , when the distance between the two particles was L i j ≤ 1 2 (D i p + D j p ), the particles collide. The equation for resolving the particle collision time is given by: The equation can be simplified to a quadratic one in one variable: Two solutions t c, 1 , t c, 2 were acquired by solving the equation, where the collision time ranged from 0 to t p (0 < t c < t p ). If both solutions were within this range, the smaller solution was considered the particle calculation time step. If neither solution was in this range, the particles do not collide in the calculation time step.

Particle collision solving process
The particle collision solving procedure is shown in Fig. 5, which begins with loading the model's parameters, including the number of input particles, grid parameters, and boundary conditions. The particle velocity was calculated using the particle's position and accept force, and the fastest particle velocity in the current calculation time step was determined. Using the sum of the center particle displacement and the fastest moving particle displacement, the radius of the space-enclosing sphere was calculated. With n particles, n − 1 collision lists were established and updated before calculating the time step. The particles moving out of the space surrounding the sphere were removed from the original list, and new particles were input. Assuming that the collision search list was established after the center particle i searched and judged particle j, which was established as the center particle to establish lists; consequently, it was unnecessary to re-judge particle i, which may be directly mapped from the entire collision list.
As far as the calculation time step is concerned, the earliest collision time for each list was added to the summary table. As the calculation time in the component calculation domain progressed, and collisions occurred, the list was updated again; if no collisions happened, the whole list was updated by subtracting the finished calculation time from the original collision time.

CFD-DEM coupling analysis and particle collision search algorithm verification
The forward search algorithm for particle collisions was written using Delphi programming language on DEM platform. Additionally, the viability of the suggested algorithm was determined by comparing it to the existing DEM software. The simulation utilized an Intel Core i7-9700 @ 3.00 GHz eight-core processor, a main graphics card AMD Radeon RX 550 Series, and 16 GB of computer memory.  In order to prove the efficiency and accuracy of the particle collision forward search algorithm (DEM-M algorithm) is proposed in this study. The DEM-M algorithm and the existing DEM algorithm are used to simulate the collision process of two particles, and the results are compared and analyzed.
The distance between two circle centers as 50 mm is shown in Fig. 6, the radius of two particles A and B was 5 mm, the density was 2500 kg/m 3 , Poisson's ratio was 0.25, and the elastic modulus was 210 GPa. Working condition 1 was defined as: particles A and B move uniformly. The initial velocity of the particle A was u A 3 m/s in the positive x axial direction, and the initial velocity of particle B was u B 1 m/s also in the positive x-direction. In contrast, working condition 2 was defined as particles A and B show variable motion. The initial velocity of particle A along the positive x axial direction was u A 3 m/s, and the acceleration along the same direction was 10 m/s 2 . The initial velocity of particle B along the positive x axial direction was u B 1 m/s, and the acceleration along the same direction was 10 m/s 2 . The velocities of the A and B particles after the collision were considered as u A and u B , respectively, and the collision force of particles A and B was F cn .
The DEM and DEM-M methods simulated the two particles' chase and collision problems, respectively. The calculation time steps were taken as 10 −4 , 10 −5 , and 10 −6 s, and the calculation time was frequently set to 0.1 s; results of the collision time, collision force, collision velocity, collision position were acquired are shown in Table 2 together with the analytical solution as the reference value.
As shown in Fig. 7, when the DEM method calculated the time step of 10 −4 s under the condition of uniform velocity and the calculated results were compared with the analytical solution, the errors in collision velocity and collision time were both greater than 30%, and the error in collision force reached 138.61%. The maximum error in velocity ranged from 1.47%, − 2.96% when the calculation time step was 10 −5 s, the maximum error in collision position was 2.19%, and the error in collision force was 2.24%. The error in collision time was 2.00%. The errors in the simulation results were less than 3.00% for each parameter. The maximum error in velocity ranged from 0.40%, 0.16% when the calculation time step was 10 −6 s, the maximum error in collision position was 0.05%, the error in collision force was 1.07%, and the error in collision time was 0.50%. The simulation results have errors of less than 2.00% for each parameter. The working condition's A2 and A3 calculation timings were 27.73 s and 44.07 s, respectively. The aforementioned analysis reveals that the calculation time step significantly impacts the calculation accuracy of the existing DEM method. The calculation efficiency reduces significantly as the time step lowers. The errors in different calculation time steps between the calculation results and the analytical solution obtained using the DEM-M method were less than 2.00%, with the maximum error in velocity being 1.87%, the maximum error in collision position being 0.21%, the maximum error in collision force being 1.36%, the error in collision time being 1.50%, and B1, B2, and B3 calculation times of the working condition being 23.01%. The preceding analysis reveals that the calculation time step has no effect on the calculation accuracy of the improved DEM-M method of the particle collision forward search algorithm suggested in this paper. For example, even if the B2 and A2 calculation time steps were equal, the calculation time is reduced by 2.59 s (9.34%) when the suggested algorithm is used. Additionally, when the calculation time in step B1 was one order of magnitude greater than the calculation time in step A2, the calculation time was reduced by 4.72 s. (17.02%). As shown in Fig. 8 that when the calculation time step was 10 −4 s for DEM method under variable speed conditions, the error between the calculated results and the analytical solution was significant, and the velocity, collision force, and collision time errors were all greater than 30%. The calculation results of the physical parameters were less than 3.00% when the calculation times for steps C2 and C3 were 10 −5 and 10 −6 s, respectively, and the calculation times were 32.53 s and 46.21 s. As the time step is decreased, the calculation accuracy increases while the calculation efficiency decreases. The error between the calculation results and the analytical results acquired in this study was less than 2.00% in different calculation time steps, utilizing the improved DEM-M method, with the maximum error in velocity being − 1.37%, the maximum , which was similar to calculation time step C2, was reduced by 3.33 s (10.23%). The calculation time for D1, which was an order of magnitude higher than the C2 calculation time step, was reduced by 6.03 s (18.54%). The above analysis demonstrates that the traditional DEM method is consistent with the forward search algorithm's conclusions and laws, strengthening the DEM in simulating variable and uniform motion.
In summary, DEM-M method suggested in this study was utilized to simulate particle motion that was both uniform and variable numerically. The size of the chosen time step has little effect on the calculated results. When the calculation time step was one order of magnitude larger than that of DEM, uniform motion calculation efficiency increased by 17.02%, and variable motion calculation efficiency increased by 18.54%, effectively resolving the contradiction between time step selection and calculation accuracy and efficiency.

Multiple particle collision examples
The initial positions of the particles are shown in Fig. 9, where 12 particle spacings of 20 mm were arranged in a square, numbered as a, b, c, e, f, g, h, i, j, k, and l, respectively. The particle material properties were the same as those in the calculation example in Sect  (20, − 15), (20,15), (15,20), and (− 20, 15), respectively. Continuing particle collisions are too complicated for analytical solutions to be acquired.  Table 3, 37 particle collisions happened when the calculation time step was 10 −6 s, and 41 particle collisions occurred when the calculation time step was 10 −7 s and 10 −8 s. The analysis reveals that the DEM missed four collisions when the calculation time step was 10 −6 s, rendering unreliable results. In comparison, when the calculation time steps were 10 −7 and 10 −8 s, the first collision time of particles happened, and the error between the collision position and the analytical solution was less than 1.00% as shown in Fig. 10. The respective calculation times were 239.18 and 1532.56 s. Consequently, as the calculation time step decreased, the calculation efficiency of the numerical simulation of multi-particle collisions decreased significantly.
When the calculation time steps were 10 −5 , 10 −6 , and 10 −7 s, the improved DEM-M numerical simulation resulted from 41 collisions, Table 4 provides the results. As shown in Fig. 11, the error between the collision position and the analytical solution was less than 1.00% for the first collision time calculated using different time steps. The relative calculation times were 193.75 s, 198.05 s, and 247.23 s. When the calculation time steps of DEM and DEM-M were 10 −7 s, the calculation time of DEM-M was 8.05 s longer than that of DEM due to the collision time of forward searching particles in each calculation time step. When the calculation time step of DEM-M was increased by one order of magnitude (10 −6 s), the calculation time was reduced by 41.13 s, implying that the calculation time was shortened by 17.19%. When the calculation time step of DEM-M was increased by two orders of magnitude (10 −5 s), the calculation time was reduced by 45.43 s, implying that the calculation time was shortened by 18.99%. The particles with the higher number of collisions were k, e, g, i, a, and c. Particles a, k, e, g had 4, 12, 12, and 6 collisions, respectively. As shown in Fig. 12, particles i and c had two collision trajectories (the dotted line shows DEM-M and the solid line shows DEM). In the simulation of multi-particle collisions, the DEM-M method suggested in this study is capable of adopting long-time steps, which improves calculation efficiency and ensures calculation accuracy.
In conclusion, this study chose 10 −5 , 10 −6 , and 10 −7 s time steps to conduct numerical simulations of multiple particle collisions, and accurate calculation results were acquired. The selection of the calculation time step was one order of magnitude greater than DEM method, decreasing the calculation time by 18.99%.

CFD-DEM coupling analysis example
As shown in Sect. 4.2, CFD-DEM and CFD-DEM-M methods were adopted to numerically simulate the coupling of multiple particles and fluid. The calculation fluid domain was a 400 mm × 400 mm × 400 mm   cube. The fluid density ρ 1.0 g/cm 3 , viscosity μ f 0.1 g/(cm × s), and initial state of the fluid were static. Regular rectangular grids were used to discretize the computational domain. The upper boundary was set as an open boundary, while the other walls adopted non-slip boundary conditions. The geometric parameters, material, initial position, and calculation parameters of the particles were identical to those described in Sect. 4.2. The fluid calculation time step in CFD-DEM was chosen 10 times the particle calculation time step. The CFD-DEM-M particle calculation's initial time step was equal to the fluid calculation's initial time step. Based on the coupling convergence condition, the real calculation time step of the particles and fluid in the coupling calculation process was adjusted and corrected in real time. Table 5 summarizes different time steps set in the simulation. Results of the CFD-DEM calculation demonstrate that the number of particle collisions acquired in the calculation time steps A1, A2, and A3 was 31, 36, and 36, respectively. The analysis illustrated that the A1 working condition was unreliable due to the deceleration of particles caused by fluid resistance, which resulted in five missed particle collisions and a reduction in the number of collisions by five times compared with the example in Sect. 4.2. The calculation times for each working condition were 27, 53, and 109 min, respectively. Only when the particle calculation time step was less than or equal to 10 -7 s, the CFD-DEM coupling calculation could achieve an accurate value.
The CFD-DEM-M calculation results described that the number of particle collisions achieved by the calculation time steps B1, B2, and B3 was 36, and the calculation times were 42.5, 47.0, and 67.2 min, respectively. The change in the initial time step of particle calculation did not affect CFD-DEM-M coupling calculation accuracy.
The time step of the fluid calculation A2 was one order of magnitude greater than the initial time step of the fluid calculation B3. The fluid solver initiated 10 5 times, the particle solver started 106 times, and the fluid and particle solvers both initiated 1,001,783 times. The particle solver A2 started 1783 times more than the particle solver B3 because the DEM-M algorithm corrected the initial time step in accordance with the coupling convergence condition, and the real calculation time step was smaller than the initial time step with B3, taking 14.2 min more than the A2 calculation. The initial time step of the fluid calculation B2 was equal to the fluid calculation time step A2, while the initial time step of the particle calculation B2 was one order of magnitude larger than the particle calculation A2. The average number of B2 startups of fluid and particle solvers was 105,329. The number of solver startups was 5329, more than that of B2 and A2 fluid solvers.
Additionally, the number of B3 multi-starts of the solver increased from 1783 to 5329 compared with A2 owing to the expansion of the initial time step of particle calculation and the non-convergence of the particle solution. B2 took 6 min, less than the calculation time of A2, and the calculation time was reduced by 11.30%. The initial time step of the fluid calculation B1 was one order of magnitude larger than that of the fluid calculation A2. Moreover, the initial time step of the particle calculation was two orders of magnitude larger than that of the particle calculation A2. Both B1 fluid and particle solvers started 13,254 times, less than the A2 calculation time of 10.5 min, shortening the calculation time by 19.80%. Figure 13 illustrates A2 and B1 motion trajectories of particles, a, c, e, g, i, and k are basically consistent.
In conclusion, CFD-DEM-M method was utilized to numerically simulate the coupling of multiple particles and fluid. The initial time steps of fluid and particle calculation were 10 −5 , 10 −6 , and 10 −7 s, which could achieve accurate results, and the calculation time was 19.80% lower than that of CFD-DEM method.

Conclusions
When the particle calculation time step was too small in the CFD-DEM coupling calculation employed in the study, the calculation efficiency was low. However, a significant number of particle collisions were missed when the selected time step size was too large. The following conclusions were obtained through algorithm research and example verification: The study suggest an improved particle collision search algorithm in which the target particle is used as the center particle, and the displacement between the center particle and the fastest moving particle in the calculation domain is considered the radius to construct the search ball. The particles that could collide were screened to create the search list, which reduced the particle collision search time. The forward search method was utilized to evaluate the particle collision, avoiding the error produced by the huge time step selection in the traditional DEM and ensuring an accurate description of the particle collision. Two-particle and multi-particle simulations' motion and collision results revealled that changing the particle calculation time step does not affect the particle collision calculation accuracy. The calculation time of two particles under uniform and variable speed conditions was lowered by 17.02% and 18.54%, respectively, compared to the traditional DEM method. The reduction in the calculation time was equal to 18.99% in multi-particle numerical simulation.
A coupling algorithm was developed to automatically adjust the time step of the particle and fluid calculations. The initial time step for particle and fluid simulations can be identical. The fluid calculation time can be adjusted in real-time based on the particle propulsion time. The automatic adjustment of the fluid calculation time step with the particle calculation time was realized, and the particle and fluid coupling calculation efficiency was improved. The numerical simulation of the multi-particle and fluid coupling example demonstrated that the algorithm's calculation accuracy was less affected by the initial calculation time step of particles and fluid and that the algorithm's calculation efficiency can be improved by using a large particle calculation time step. The algorithm's calculation time was decreased by 19.80% compared to the traditional CFD-DEM method.
Author contributions XW took part in methodology, software, visualization, and writing of the original draft. SW involved in methodology, funding acquisition, and supervision. MW took part in data curation and formal analysis. MW and XL involved in data curation and language check. DS and CS involved in software and results validation.

Conflict of interest
The authors declare no competing financial interests or personal relationships that could have influenced the work reported in this paper.