## 2 Mathematical Modeling
### Computational Fluid Dynamics
#### 2.4.1 Continuous Phase
For the spherical particles suspended in the fluid, the drag force is given as follows:  
\[F_{D}=\frac{18\mu C_{d}}{\rho_{p}d_{p}^{2}}\cdot\frac{Re_{p}}{24} \tag{14}\]  
Here, particle diameter is denoted as \(d_{p}\), fluid molecular viscosity as \(\mu\), and \(Re_{p}\) is the Reynolds number of the particle, given as follows:  
\[Re_{p}=\frac{\rho d_{p}\big{(}\ u_{p_{i}}-\ u_{i}\big{)}}{\mu} \tag{15}\]  
The drag on smooth particles is due to [PERSON] and [PERSON] [37], given as follows:  
\[C_{D}=a_{1}+\frac{a_{2}}{Re_{p}}+\frac{a_{3}}{Re_{p}^{2}} \tag{16}\]  
where the coefficients \(a_{1}\), \(a_{2}\), and \(a_{3}\) are constants. The instantaneous velocity for the fluid flow is expressed as \(u_{i}=\overline{u}_{i}+u^{\prime}_{i}\); where \(u_{i}\) represents the instantaneous velocity, \(\overline{u}_{i}\) is the mean velocity component, and \(u^{\prime}_{i}\) is the fluctuating velocity component. Since the LES works on the filtering approach, the directly resolved LES field very much resembles \(u_{i}\) (same as used in Equation (13)), thereby requiring no approximation to account for the particle dispersion due to velocity fluctuations. However, for the unresolved part of the flow, stochastic modeling is required to be introduced using the discrete random walk (DRW) model--taken as discrete piecewise constant functions of time--and is given as follows:  
\[u^{\prime}_{i}=\zeta\sqrt{\frac{2}{3}k} \tag{17}\]  
where \(\zeta\) is a normal distributed random number and \(k\) is the sgs turbulent kinetic energy. Meshing  
The multi-block approach has been used, wherein the blocks are split to fill the entire domain with the industry-leading tool ICEM CFD. The grids obtained are non-orthogonal and fit well to the complex geometry. To minimize the effects of numerical diffusion, the cell faces are aligned in the direction of the fluid flow. The mesh closer to the boundaries and near the geometric axis is refined to capture the near-wall flow physics and impacts of the precessing vortex core near the geometrical axis, respectively. Figure 2 represents the details of the surface mesh of cyclones A, B, and C, consisting of 1.287 million, 1.092 million, and 1.114 million hexahedra, respectively. We have checked the results at two levels of the grid for all four models. The refinement (by more than 2.5 times that of the coarse mesh) was performed on a global basis, which also included further refinement of the grid at the boundaries at higher mesh counts. The maximum aspect ratio was around 68 at the conical walls of the refined mesh. When comparing the two levels of the grid, a maximum difference of 2.23% in pressure drop and 4.17% in the cut-off particle size was observed in cyclone model C. The results from the coarse mesh have only been presented. Numerical Settings  
The RNG \(k\)-\(\epsilon\) turbulence model was initially used for the steady-state simulations until the residuals became (quasi-) steady with the number of iterations. The standard wall function was selected, and the solver tolerance was set at \(10^{-5}\). This method significantly reduced the computational efforts, and the resulting flow field served as a starting point for LESs. The simulations were run till the monitored flow variables became stable with time. To couple the pressure and velocity, we use the NITA (Non-iterative Time Advancement) algorithm (in accordance with [38]), followed by the fractional step scheme. NITA utilizes more memory compared to the iterative time advancement technique, but  
Figure 2: Representation of the surface mesh of models A–C (from left to right), consisting of hexahedra cells. it saves significant computational time. The pressure was discretized using the PRESTO! Interpolation scheme (as in [39]), which is optimized for the swirling flows with steep pressure gradients. For momentum, a bounded central differencing scheme was used, and for time enhancement, a bounded second-order implicit scheme was used.