CFD-DDPM coupled with an agglomeration model for simulation of highly loaded large-scale cyclones: Sensitivity analysis of sub-models and model parameters

The Dense Discrete Phase Model coupled with an agglomeration model is developed and validated for the simulation of industrial cyclones with high solid loads. The performance of the model is influenced by sub-models, model parameters, and numerical parameters. To optimize the performance of the present CFD model, an extensive sensitivity analysis was performed, varying one sub-model or parameter at a time, and systematically assessing the effect on the results through comparisons with measured pressure drop and sepa- ration efficiency of a highly loaded pilot-scale cyclone. The investigation shows that the turbulence model and particle-particle restitution coefficient have the strongest influence. This study concludes with the recommen- dation of a set of sub-models, model parameters, and numerical parameters providing the best prediction of the hydrodynamics of large-scale highly loaded cyclones. In addition, the impact of some operating conditions on the performance of a large-scale highly loaded cyclone were examined.


H I G H L I G H T S G R A P H I C A L A B S T R A C T
• Simulation of highly loaded cyclones: effects of sub-models and model parameters. • Particle-particle restitution coefficient and turbulence model influence the most. • The impact of particle-particle interaction force modeling is minor. • Most suitable sub-models, model parameters, and numerical parameters are suggested. The Dense Discrete Phase Model coupled with an agglomeration model is developed and validated for the simulation of industrial cyclones with high solid loads. The performance of the model is influenced by submodels, model parameters, and numerical parameters. To optimize the performance of the present CFD model, an extensive sensitivity analysis was performed, varying one sub-model or parameter at a time, and systematically assessing the effect on the results through comparisons with measured pressure drop and separation efficiency of a highly loaded pilot-scale cyclone. The investigation shows that the turbulence model and particle-particle restitution coefficient have the strongest influence. This study concludes with the recommendation of a set of sub-models, model parameters, and numerical parameters providing the best prediction of the hydrodynamics of large-scale highly loaded cyclones. In addition, the impact of some operating conditions on the performance of a large-scale highly loaded cyclone were examined.

Introduction
Cyclones operating with high solid loads are a vital part of building material industries such as cement production [1][2][3][4], stone wool production [5,6] and also processes employing circulating fluidized beds (CFB) [7,8], and Fluid Catalytic Cracking applications [9]. The geometrical simplicity and the low costs in construction, operation, and maintenance have made cyclones the first option when gas-particle separation and/or heat exchange is required.
CFD simulation of such a device has become an invaluable tool for optimizing cyclone design and operation [10][11][12]. Despite the simplicity of highly loaded cyclones, the complex particle-particle and gas-particle interactions, highly turbulent flow, and particle agglomeration, make the CFD simulation quite complicated [13]. It becomes even more complicated in the simulation of large-scale cyclones where the use of coarse grids (that might be unable to resolve mesoscale structures) is inevitable due to the limitation of computational resources [13].
Conventionally there are two common strategies to model multiphase flow in CFD simulation of dense dispersed gas-solid flow: (1) Two Fluid Method (TFM) [14], and (2) coupling of CFD with the Discrete Element Method (DEM) [13]. Both strategies treat the gas phase as a continuum medium in an Eulerian framework while they are different in how the particles are treated: as a continuum in the Eulerian framework in the TFM method, and as discrete particles in the Lagrangian framework in CFD-DEM. However, they have limitations in the simulation of highly loaded large-scale cyclones. CFD-DEM is limited to lab-scale cases due to its limitation in coarse-graining, the number of grains to contain, and its necessity to simulate with very small time-steps to capture particle-particle collisions that make it unfeasible for simulation of highly loaded large-scale cyclones [15]. On the other hand, TFM cannot easily handle the polydispersity of particles which is crucial for the simulation of the hydrodynamics of cyclones and the prediction of their performance [14].
A more recent strategy called Dense Discrete Phase Model (DDPM) was developed to avoid these issues [16,17]. DDPM is a cost-effective solution for industrial dense particle-laden flows. Unlike DEM which directly resolves particle-particle interactions, DDPM models these interactions based on the Kinetic Theory of Granular Flow (KTGF) so that larger time steps and larger grain-to-particle ratio can be applied to accelerate the calculations tremendously [18][19][20][21]. In addition, DDPM can easily handle polydispersed particles [22]. However, there is a lack of criteria in the open literature to guide a proper choice for the particleparticle interaction forces modeling in DDPM, time-step size, and grainto-particle ratio in particular for the simulation of a highly loaded largescale cyclone.
The use of widely accepted turbulence models for cyclones, i.e., Reynolds Stress Model [23][24][25] and Large Eddy Simulation [26][27][28], is not practical for CFD simulation of industrial-scale cyclones under high loading of particles due to high computational cost and instability in case of unstructured mesh that is commonly used for large-scale simulations [22]. In addition to RSM and LES as widely accepted turbulence models for cyclones, novel Eddy Viscosity Models (EVMs) sensitized to  ω Specific dissipation rate p particle i ith direction pi ith particle 1 Collector particle 2 Contributor particle rotation and curvature have attracted much attention recently due to their superiority over LES and RSM in terms of computational cost and robustness. It was proven in several research works [22,[29][30][31][32][33] that this sensitization can overcome weaknesses of EVMs in the simulation of cyclones so that they can provide quite reliable results which is competitive with more complicated turbulence models, while these sensitized models are more robust and require much less computational resources [18][19][20][21][22]. Considering the industrial scale of the studied case, and the fact that in LES, relatively fine girds are required to resolve large eddies, LES is not favored. Also, LES cannot be coupled with DDPM in Ansys Fluent. On the other hand, regarding the robustness of EVMs and their superiority over RSM in terms of computational cost which are of great importance in our case (industrial-scale geometry), and its proven ability to capture mean flow field in cyclones when sensitized for flow rotation and curvature [1,[18][19][20][21][22], EVMs are preferred. All available EVMs (standard, RNG, RNG, and realizable versions of k-ε and standard and shear stress transport (sst) versions of k-ω) can be sensitized to flow rotation and curvature, but which perform better under highly loaded cyclones and how significant is the effect of this sensitization have yet to be studied. With multiphase models based on the Kinetic Theory of Granular Flow (KTGF), such as the DDPM, using conventional homogeneous drag models in coarse grids, as is required in large-scale simulations, fails to adequately capture the hydrodynamics of multiphase flow [34,35]. For this reason, heterogeneous drag models (also known as sub-grid drag models) that take into account mesoscale and heterogeneous structures are recommended [13,[34][35][36][37][38]. However, none of the sub-grid drag models presently available in the literature has been developed for cyclones, and there is a lack of guidance in choosing a suitable drag model for the simulation of highly loaded cyclones among the available drag models.
Agglomeration is a vital phenomenon in particle-laden systems with very fine particles e.g., in cyclones. Thus, a model accounting for agglomeration is essential if a valid prediction of cyclones is desired [39,40]. Furthermore, it is crucial to investigate the parameters affecting agglomeration extent, such as the particle-particle restitution coefficient, and how they impact simulation predictions.
The impact of particle-wall interaction modeling parameters on the performance of low-loading cyclones has been reported to be minor [26,28,41], but the extent to which it affects highly loaded cyclones has not been investigated. In addition, no criteria are provided in the literature to help determine whether including particle rotation, Saffman force, Magnus force, and particle turbulent dispersion is crucial for the simulation of highly loaded cyclones or not, while their relevance has been noted in the case of low loading cyclones [28,[42][43][44][45][46][47].
In the present work, sensitivity analyses were performed on submodels (models for particle-particle interaction force calculation, turbulence models, drag models and other gas-particle interactions), model parameters (particle-particle restitution coefficient and particle-wall rebound coefficient), and numerical parameters (parcel/grain size, time-step size and order of discretization) of a hybrid multiphase model developed and validated for simulation of large-scale highly loaded cyclones. The purpose of the sensitivity analysis is to assess the significance of each of the sub-models, model parameters, and numerical parameters, to clarify their impact, and to optimize the performance of the numerical model. In addition, the effect of some operating conditions, such as solid loading, particle size distribution and gas underflow on the performance of a highly loaded large-scale cyclone were studied.

Numerical model description
In the present work, a hybrid Eulerian-Lagrangian multiphase model, Dense Discrete Phase Model (DDPM) is used. The model implicitly approximates particle-particle interaction using the Kinetic Theory of Granular Flow (KTGF) and is coupled with an agglomeration model [39,48]. The gas phase is treated in an Eulerian framework where the Navier-Stokes equations including phase volume fraction are solved. The particle phase is treated in a Lagrangian framework while the particle-particle interaction is modeled based on the Kinetic Theory of Granular Flow (KTGF) [49]. The model is described as follows:

Gas phase
The continuity and momentum conservation equation for the continuous phase (gas phase) in DDPM multiphase modeling coupled with dispersed turbulence models are as follows: F Dc accounts for the force particles apply on fluid due to the drag and is modeled as follows:

Turbulence modeling
In the present work, we tested different versions of k-ε and k-ω turbulence models in their dispersed form: standard, RNG, RNG for swirl dominate flow, and realizable versions of k-ε and standard and shear stress transport (sst) versions of k-ω. In most cases, the used turbulence models included a term to sensitize the model for curvature correction and rotation. Here the turbulence model used as the base case (sst version of k-ω turbulence model sensitized for curvature correction, k-ω sst cc) is presented while details of other tested turbulence models can be found in the Ansys Fluent theory guide [16].
to calculate k c and μ t, c and consequently, turbulent stress tensor, τ c t , two transport equations (one for k and one for ω) must be solved: Definitions of terms and parameters in these equations and how curvature correction function (f r ) is modeled can be found Ansys Fluent theory guide [16]; in short, f r (curvature correction function) is modeled by a modified version of empirical function suggested by Spalart and Shur [50] to account for streamline curvature and system rotation effects.

Solid/particle phase
In DDPM multiphase model, the equation of motion of each parcel is expressed as follows: where F p− p is modeled based on the Kinetic Theory of Granular Flow (KTGF). In the KTGF model, instead of resolving all particle collisions, the interactions are averaged for the solid-gas flow. Note that α c in gas phase conservation equations are updated each time step by mapping particle location obtained from Eq. (7) (lagrangian framework) to the Eulerian framework. In addition to the translational particle motion equation, a particle rotation equation can be considered as follows: 2.3. Particle-particle interactions 2.3.1. Particle-particle interaction force based on KTGF The particle-particle interaction forces are obtained by the solids stress tensor (τ s ̿ ) as defined in Eqs. (9) and (10).
Estimation of these viscosities (bulk, collisional, kinetic and frictional) can be done using different closure models details of which can be found in the Ansys Fluent theory guide [16]. In short, granular viscosity, granular bulk viscosity, solid pressure, granular temperature, frictional viscosity, frictional pressure, and radial distribution terms must be modeled using a closure model. Closure models tested in the present work for each of these terms can be found in Table 2.

Agglomeration
In addition to particle-particle interaction force modeling based on KTGF, particle agglomeration that was found to be crucial in the case of cyclones handing very fine particles (d < 15 μm) [39,40] was modeled based on a stochastic Lagrangian inter-particle collision model proposed by Sommerfeld [39,48]. Details of the agglomeration model can be found in our published work [22]. In short, whether two particles in a cell collide or not is judged based on calculated collision probability and efficiency, and whether two collided particles agglomerate or rebound is judged by comparing the relative velocity of the particles with a critical velocity (agglomeration occurs if . This critical velocity is defined as follows:

Drag force
To close the momentum conservation equation (Eq. (2)) and for particle tracking (Eq. (7)) the drag force between fluid and particles must be modeled. In dense systems such as highly loaded cyclones, the drag force is influenced by the presence of other nearby particles which makes its modeling quite complex. In such cases, drag modeling is mainly based on empirical observations, and has limited ranges of validation.
Therefore, since in the literature there is not a model designed specifically for highly loaded cyclones, to find the best-suited drag model, in the present work, a set of models have been tested. This set includes five homogeneous drag models (Wen-Yu [51], Gidaspow [52], TGS [53], BVK [54], DiFelice [55]) and three heterogeneous drag models (Sarkar [34], Pirker [56], Igci [36]). Except for Wen-Yu and Gidaspow, all other models were implemented to Fluent using UDFs.

Other gas-solid interaction forces 2.4.2.1. Magnus force.
Rotating particles in a fluid generate a Magnus lift force due to pressure differential along their surfaces as follows [60]: where C RL is the rotational lift coefficient and is modeled based on Oesterle-Bui-Dinh [61] expression in the present work.

Saffman force.
Particles with low Reynolds numbers, especially particles finer than 1 μm, experience a Saffman lift force as a result of shear that is modeled by the following expression [62]:

Rotational drag force.
An additional drag force is generated for rotating particles represented as C ω in the particle rotation equation which is modeled by a function of the rotational Reynolds number [63]. The effect of these three additional gas-solid interaction forces on the gas phase is summarized as a source term (S) in the momentum conservation equation (Eq. (2)) and their effect on the particle phase is summarized as F other in particle motion equation (Eq. (7)).

Experimental setup and simulation details
The experimental setup described in Mirzaei et al. [22] was simulated in the present study; a pilot-scale cyclone (1.6 m in diameter) shown in Fig. 1(a) was tested under industrial operating conditions (particle load, material, and size shown in Fig. 1(b)) and room temperature. The experiments were conducted at various operating conditions with different gas flow rates and solid flow rates. Three of these operating conditions are used in the present study. In Table 1, the operating conditions of these three cases in addition to measured separation efficiency (defined in cyclones as separation efficiency (%) = × 100)) and pressure drop are summarized. Unless it is specifically mentioned, case-1 (base case) operating condition is used in sensitivity analyses. It is worth mentioning that the term "highly loaded cyclone" refers to a cyclone operating under a dense gas-solid flow [64]. The approximate borderline between dilute and dense gas-solid systems is a solid volume fraction of 0.001. Since in the cyclone studied here, the solid volume fraction, especially close to the walls, is >0.001 [65,66], the cyclone is considered a highly loaded cyclone.
In the present study, 4.2 million polyhedral cells with prism layers were used to discretize the geometry. This grid was found to provide the mesh-independent predictions in the present model and geometry. The solver employed in the present work was ANSYS Fluent 2019R3; details of the solver, numerical methods, and mesh independency study is discussed in detail in our previous work [22].
To find out the best-suited sub-models, model parameters, and numerical parameters of the developed model (DDPM coupled with agglomeration modeling) for the simulation of highly loaded large-scale cyclones, several sensitivity analyses were performed. Sensitivity Table 1 Operating conditions of conducted tests used in the present study. Unless it is mentioned case-1 (base case) operating condition is used in sensitivity analyses.  First order analyses include particle-particle interaction modeling (based on KTGF), turbulence models, drag models, particle-particle restitution, particlewall rebound coefficient, discretization method, grain-to-particle size ratio, grid-to-parcel size ratio, and time-step size. Table 2 summarizes the set of base case models and parameters chosen based on preliminary simulations. Starting from this set, one of them was varied at a time and its effect on the results was assessed. In addition to models and parameters in the base case, all other tested models and parameters are summarized in Table 2. To analyze the results, mainly pressure drop, and separation efficiency predictions were compared with experimental measurement and the predictions in the using base case model.

Results and discussion
The results obtained from sensitivity analyses on sub-models, model parameters, numerical parameters, and operating conditions are presented in this section. Note that throughout the study if the deviation of a predicted parameter from the base case is below 2.5%, the effect is considered minor.

Turbulence modeling
The sst (shear stress transport) version of the k-ω turbulence model sensitized for curvature correction (k-ω sst cc) has been validated for the simulation of highly loaded large-scale cyclones [22] where the use of more complicated turbulence models such as RSM or LES is not feasible due to their instability and computational cost. In addition, in a study of a highly loaded lab-scale case (discussed briefly in the supplementary material), the sst version of the k-ω turbulence model was proven to perform as well as RSM, capable of predicting performance parameters of a lab-scale highly loaded cyclone in agreement with experimental measurements (at least for pressure drop predictions that measurements were available) [67]. Therefore, this turbulence model was chosen as the base case in the present study.
In the sensitivity analysis of turbulence modeling in the present study, different versions of the k-ε turbulence model (standard, RNG, RNG for swirl dominate flow, and realizable) sensitized for curvature correction, the standard version of k-ω turbulence model sensitized for curvature correction, and sst version of k-ω turbulence model without curvature correction term were tested. The results were then compared with the base case (k-ω sst cc) and experimental measurements of pressure drop and separation efficiency (see Fig. 2). The following can be drawn from the comparison in Fig. 2: • The use of the curvature correction term to sensitize the turbulence flow for rotation is crucial as excluding it leads to highly under prediction of separation efficiency. • The sst version of the k-ω turbulence model is superior to the standard version as the standard version highly underpredicts the separation efficiency. • Among k-ε turbulence models, only two RNG versions provide reasonable predictions while the other two highly underpredict the separation efficiency.
The comparison in Fig. 2 shows that only k-ε RNG cc and k-ω sst cc provide reasonable predictions in particular in the prediction of separation efficiency. The sensitization of EVMs for flow rotation and curvature was reported to overcome weaknesses of EVMs in the simulation of cyclones [18][19][20][21][22], and the presented results confirm it in the case of the k-ω sst turbulence model as the one with the correction provides quite reliable predictions unlike the version without curvature correction. Also, The RNG version of k-ε was originally designed to include the effect of swirl on turbulence in k-ε and to improve accuracy for swirling flows and the presented results show that its impact is significant; unlike To further compare these models they were tested for two other operating conditions (operating conditions 2 and 3 with the same gas flow rates and different solid flow rates). As shown in Fig. 3, unlike k-ω sst cc, k-ε RNG cc is unable to capture the trend of improvement in separation efficiency when the solid load is increased (reported in the literature and observed in the measurements) that makes k-ω sst cc the superior turbulence model among all tested models. Even the other version of k-ε RNG cc that is modified for swirl dominate flows is not able to capture such a trend. In addition, it is worth mentioning that k-ω sst cc is almost 3 times faster than k-ε RNG as it requires fewer iterations to reach convergence in each time step.
The flow pattern in the cyclone obtained from each of the tested turbulence models are presented and compared in Fig. 4 and Fig. 5. Although all models capture the overall flow in cyclones (swirling down outer vortex and swirling up inner vortex), different models provide considerably different flow patterns, especially in the conical part of the cyclone. Flow pattern (tangential and axial velocity profiles) measurements of highly loaded cyclones that are not available for now can facilitate further comparison of turbulence models. It is worth mentioning that due to both instability and limited computational resources, RSM and LES have not been included in the present study but can be investigated in the future to further study the effect of turbulence modeling in highly loaded industrial-scale cyclones.

Particle-particle interaction modeling
As mentioned in the model description, several terms (granular viscosity, granular bulk viscosity, frictional viscosity, solid pressure and radial distribution) are required to be modeled to predict particleparticle interaction force and there are a few available options for each of these terms. Although the effect of using different models for each term was studied for different fluidization systems in the literature [68][69][70], such a study was missing for highly loaded cyclones. Therefore, a sensitivity analysis on each of these terms was done to figure out which model is more suitable and how they would affect the predictions in the case of highly loaded cyclones.
Simulation settings of cases used for particle-particle interaction force sensitivity analysis are summarized in Table 3, and the predictions of pressure drop and separation efficiency with different settings are compared in Fig. 6. The results show that the effect of using different closure models for each of the terms is minor both in the prediction of pressure drop and separation efficiency (<1% effect on pressure drop and <2% effect on separation efficiency). Even excluding the particleparticle interaction force (by deactivating all terms in the particleparticle interaction force modeling) did not affect the predictions  considerably, therefore, choosing between different available options to model particle-particle interaction force is not critical for the simulation of highly loaded cyclones (at least in the range studied in the present work <1 kg solid/kg air).
The solid volume fraction in near-wall regions exceeds 0.001 (The approximate borderline between dilute and dense flow) meaning that we have dense gas-solid flow there, so it was expected to have considerable particle-particle interactions near the walls affecting hydrodynamics of the multiphase flow (solid volume fraction distribution can be found in the supplementary material). However, the presented results suggest that particle-particle interaction forces (even if they are considerable) do not affect the performance parameter of highly loaded cyclones and their predictions as excluding them from the model did not affect the predictions notably.
In the present case, particle-particle interaction force calculations were not a computational load compared to other calculations, however, in cases where several solid phases are required to be included (e.g. using DDPM for polydispersed particles including energy conservation that requires additional solid phases to achieve energy balance), excluding particle-particle interaction force can speed up the computation significantly without affecting the predictions considerably.

Drag model
Several drag models including five homogeneous drag models (Wen-Yu, Gidaspow, TGS, BVK, and DiFelice) and three heterogeneous drag models (Sarkar, Pirker, and Igci) were examined. Comparisons of predictions of pressure drop and separation efficiency are presented in Fig. 7. The predicted pressure drop seems to be not affected by the drag model significantly (<2%), but the separation efficiency predictions are more sensitive (up to 4%). Except for the Sarkar drag model which is the base case model in the present work, all other drag models over-predict the separation efficiency.
The flow patterns in the cyclone obtained from each of the tested drag models are presented and compared in Fig. 8 and Fig. 9. Except for the Sarkar drag model, all other drag models provide similar flow patterns in the cyclone. Although all drag models provide reasonable flow patterns matched with expectations from cyclones, without flow pattern measurements (that are not available at the moment in case of highly loaded cyclones) it is impossible to comment on which model does a better job.
To further compare these models, they were tested for two other operating conditions (operating conditions 2 and 3 with the same gas flow rates and different solid flow rates). As shown in Fig. 10, although all models except the Sarkar drag model overpredict the separation efficiency, they all are able to capture the trend of improvement in separation efficiency when the solid load is increased. The extent of improvement is not the same in all models (see Fig. 10 (b)). The Igci drag model predicts the greatest improvement in the separation efficiency (closest to experimental measurement).
It is widely reported that in coarse-grid multiphase modeling based on the Kinetic Theory of Granular Flow (KTGF) (as it is done in the present study) using common homogeneous drag models, in some cases, leads to failure in the prediction of the multiphase flow hydrodynamics [6,[23][24][25][26][27][28][29][30][31]. However, the result presented here shows that the effect of

Table 3
Simulations settings for particle-particle interaction force analysis. unrealistic prediction of homogeneous drag force (caused by the inability of coarse grids to resolve mesoscale structures such as particle strands in cyclones) is not strong enough to influence predictions of performance parameters of highly loaded cyclones notably (compared to the turbulence modeling and particle-particle restitution coefficient that is discussed in a later section).

Fig. 6.
Effect of using various options for particle-particle interaction modeling on the overall pressure drop and separation efficiency. Arrows and corresponding percentages show the difference between CFD simulation and experimental measurement.  Fig. 11. Note that in heterogeneous drag models β β Wen− Yu equals the H d , i.e., heterogeneity index, value. It is widely reported that when multiphase models based on the Kinetic Theory of Granular Flow (KTGF) are used with coarse grids, then using common homogeneous drag models, such as Wen-Yu, over predict the drag force and that is why heterogeneous drag models should be used to overcome this over-   Fig. 11) in BVK and TGS drag models show that they are even worse than the Wen-Yu drag model. DiFelice drag model has values close to 1 which makes it similar to Wen-Yu which overpredicts the drag. However, heterogeneous drag models (Pirker, Sarkar, and Igci) have various amounts of particles with H d = β β Wen− Yu values lower than 1 (blue particles in Fig. 11) compensating for the overestimation of drag to different extents (among heterogeneous drag models, Pirker has the lowest number of particles with H d value lower than 1 and Sarkar has the highest).
When using sub-grid drag models that are a function of the filter size, the predictions might be affected by the filter-to-grid ratio so that in those cases the ratio should be optimized [71][72][73] or instead more novel models that make H d independent of filter size should be used [74]. For  the base case drag model (Sarkar), a sensitivity analysis of the filter-togrid ratio (the common values of 1, 2, and 3 were tested) was conducted showing that predictions are not affected by the filter-to-grid ratio considerably (see Fig. 12; the variation of pressure drop and separation efficiency is <1%), so the suggested value of 2 can be used with confidence.
To find out the extent of the effect of these terms in the present case, the same simulations but with the exclusion of one of these terms (models) at a time were performed and the predictions are compared in Fig. 13. The comparison shows that none of these terms has a major impact on the prediction of performance parameters of highly loaded cyclones Unlike Particle rotation, Saffman force, and Magnus force which have been neglected in several CFD simulation studies in the literature for their minor impact, particle turbulent dispersion has been included almost in all CFD simulation studies of cyclones with low solid loading supposedly for its considerable contribution. However, as observed in the present study, in the case of highly loaded cyclones, the contribution of particle turbulent dispersion is also minor (at least compared to other factors such as drag and agglomeration) and can be neglected.
Therefore, these sub-models can be dropped from the model to accelerate the calculations. Particularly, excluding particle rotation which necessitates one ODE to be solved at each time-step for each  The results indicate that the drag force is the dominant force for highloading cyclones, compared to Magnus, Saffman, and particle-particle interaction force.

Particle-particle restitution coefficient
Particle-particle restitution coefficient appears in some terms for particle-particle interaction force equations of KTGF (for estimation of some of the parameters in Eqs. (10) and (11)) and in the agglomeration model for calculation of critical velocity. Although particle-particle interaction force was found to have a minor effect and therefore   15. Comparison of different values of particle-particle restitution coefficient in capturing the trend of improvement in separation efficiency when the solid load is increased in the case of Sarkar drag model simulations.
particle-particle restitution coefficient in its terms would not be of importance, this coefficient plays an important role in agglomeration modeling as the critical velocity and therefore the possibility of colliding particles making agglomerates is a direct function of its value. In other words, the extent of agglomeration is governed and can be adjusted by the particle-particle restitution coefficient.
Higher values of the coefficient lead to lower critical velocities so only particles colliding with smaller relative velocities would make agglomerates ( ̅→ ⃒ ⃒ cos(ϕ) ≤ U crit ), also higher values of the coefficient mean higher particle rebound velocity and energy, reducing the chance of making agglomerates. On the other hand, lower values lead to higher critical velocity so more colliding particles would satisfy the condition for making agglomerates. It is worth mentioning that in the present case, the agglomeration mostly occurs in the injection pipe which can be observed by investigating how PSD and solid volume fraction distribution are changing throughout the domain (it is discussed in more detail in the supplementary material). To see how the particle-particle restitution coefficient affects the predictions, a range of values (0.5-0.95) were tested and compared in Fig. 14. As expected, lower values lead to more agglomeration and consequently larger agglomerates in the cyclone leading to higher separation efficiency. In addition, since larger particles require less energy to be carried by gas, lower values of particle-particle restitution coefficient and consequently larger agglomerates lead to lower pressure drop as can be observed in Fig. 14. Although such a trend was already expected, the extent of these changes is now confirmed for the PSD studied in the present work.
It should be noted that the extent of the effect of the particle-particle restitution coefficient on the predictions strongly depends on the operation conditions of the cyclone. For example, in a cyclone with larger particles where agglomeration is less intense, weaker dependency is expected, while in a cyclone with finer particles where agglomeration is more intense, stronger dependency is expected. Therefore, it should be noted that the level of dependency reported here is only valid for the present operating conditions.
As can be observed, separation efficiency is significantly influenced by the particle-particle restitution coefficient (or any other parameter defining the extent of agglomeration). Therefore, to see which value would better capture the trend of improvement in separation efficiency when the solid load is increased, different values of particle-particle restitution coefficient were tested and compared for two other operating conditions (operating conditions 2 and 3 with the same gas flow rates and different solid flow rates).
The comparison shown in Fig. 15 shows that in the case of using the Sarkar drag model, the value of 0.7 for particle-particle restitution   17. Comparison of using Sarkar drag model combined with particle restitution coefficient of 0.7 and using Igci drag model combined with particle restitution coefficient of 0.8 in the prediction of separation efficiency and in capturing the trend of improvement in separation efficiency when the solid load is increased. coefficient provides the best predictions; a better match with experimental measurement of separation efficiency and the highest (the closest to the experimental measurement) improvement in separation efficiency due to an increase in the solid load. Lower values of particleparticle restitution coefficient lead to over-prediction of separation efficiency and diminishing of the extent of improvement in separation efficiency due to an increase in solid load. On the other hand, higher values are unable to capture such an improvement.
Since the Igci drag model showed promising performance in capturing the trend of improvement in separation efficiency (better than all other drag models as can be observed in Fig. 10), a similar comparison (testing different values of particle-particle restitution coefficient for operating conditions 2 and 3 with the same gas flow rates and different solid flow rates) is done for this drag model as well (see Fig. 16). By increasing the particle-particle restitution coefficient to 0.8 the over-prediction of separation efficiency by the Igci drag model observed in Fig. 7 is reduced, also the extent of improvement in separation efficiency is slightly improved.
Overall, as presented in Fig. 17, using the Sarkar drag model combined with a particle restitution coefficient of 0.7 and using the Igci drag model combined with a particle restitution coefficient of 0.8 provide the best predictions of separation efficiency and the best trend of improvement in separation efficiency when the solid load is increased. As both Sarkar and Igci drag models provide accurate predictions, further comparison of these drag models requires flow pattern measurements since they provide considerably different axial and tangential velocity profiles (see Fig. 8 and Fig. 9).
Although the effect of drag modeling is minor, Igci and Sarkar drag models provided slightly better predictions as discussed. Their superiority over the other homogeneous drag models (Wen-Yu, Gidaspow, TGS, BVK, and DiFelice) may be explained by the fact that homogeneous drag models are not able to capture the effect of sub-grid scale heterogeneous solid structures. On the other hand, their superiority over the Pirker heterogeneous drag model might be explained by the fact that Igci and Sarkar drag models are derived from simulations in a periodic domain of a generic gas-solid system, while Pirker drag model is derived from simulations of a specific two-phase system (bubbling fluidized bed). Therefore, Igci and Sarkar drag models are able to be more general and perform properly in a wider range of gas-solid systems, including cyclones that contain dense and dilute solid regions at the same time.
Considering that properties such as particle restitution coefficient, are not available or easy to measure for a specific solid material, its value can be adjusted in a valid range (dependent on the physical properties of the solid material) to improve the predictions. The effect of the particle-particle restitution coefficient on the grade efficiency curve of the cyclone is presented in Fig. 18(a). The wellknown "fish-hook" behavior in highly loaded cyclones (the observation of minima in the grade efficiency curve) [80] is observed in all particle-particle restitution coefficients (except for 0.5) as there is at least one minimum in each grade efficiency curve in the range of 1-10 μm. Lower values of restitution coefficient lead to more agglomeration and consequently higher grade efficiencies. The number of these minima and their locations (particle size and fractional efficiency where a minimum is observed) depend on the solid properties (particle-particle restitution coefficient). Note that in all cases, particles larger than 45 μm are fully separated while the grade efficiencies of particles finer than 45 μm are dependent on the restitution coefficient meaning that up to this particle size agglomeration is playing a role in the performance of the studied cyclone.
Taking samples of the particles that escaped and separated from the cyclone to obtain the experimental grade efficiency can be valuable to have an understanding of the level of agglomeration happening in the cyclone and can help to choose an appropriate particle-particle restitution coefficient value as different values provide different grade efficiency curves.
In addition to grade efficiency curves, to better see the effect of particle-particle restitution coefficient on the formation of agglomerates, the PSD of particles leaving the domain (after agglomeration) are presented for all tested restitution coefficients in these three cases in Fig. 18(b) and are compared with PSD of injected particles. As the restitution coefficient decreases, overall larger particles are formed as the curve is shifted to the right. The lower value of the restitution coefficient, the more significant difference between the PSD of particles leaving the domain and the PSD of particles injected into the domain due to more agglomeration.

Particle-wall rebound coefficient
The particle-wall collision was modeled by normal and tangential rebound coefficients (i.e., the ratio between particle rebound velocity to particle impact velocity) that is a property dependent on particle and wall material (and collision angle). Like the particle-particle restitution coefficient, the particle-wall rebound coefficient is not available in the literature for various wall and particle materials and is not easy to measure. Therefore, a sensitivity analysis was done for this parameter in the range of 0.5-1 to see how it affects predictions of important performance parameters in the case of highly loaded cyclones.
As shown in Fig. 19, its impact on the pressure drop prediction is minor. However, the impact on the separation efficiency can be divided into two ranges; (1) 0.5-0.8 where the impact is minor, and (2) 0.8-1 where an increase leads to a reduction in separation efficiency up to 10%. Higher values of particle-wall rebound coefficient mean higher particle rebound velocity and a higher chance of passing the descending outer vortex and reaching ascending inner vortex (especially in the conical part of the cyclone where the outer vortex is thin) and eventually Fig. 19. Effect of particle-wall rebound coefficient on the overall pressure drop and separation efficiency. Arrows and corresponding percentages show the difference between CFD simulation and experimental measurement. Note that both tangential and normal rebound coefficients are changed at the same time.  a higher chance of escaping with gas that leads to a reduction in separation efficiency. Although values in the range of 0.9-1 might be unphysical, the rebound coefficient can be adjusted in the range of 0.8-0.9 to improve separation efficiency predictions slightly, in particular, in cases where pressure drop predictions are fine (in cases where both separation efficiency and pressure drop are off from experiments, adjusting particle-particle restitution coefficient seems more reasonable as it affects both separation efficiency and pressure drop).
In the present study, more complex (and more accurate) particle-wall models that relate the rebound coefficient to the angle of collision have not been investigated because the constants in such models are dependent on the physical properties of the wall and the particles and are not available in the literature for the present case.

Time step
One important superiority of DDPM over other common multiphase models such as DEM and TFM is that it can run simulations with relatively large time steps [18,19] which makes this model a feasible option for the simulation of large-scale systems. However, unlike DEM and TFM, in DDPM due to its relative novelty, there is no recommendation to decide the time step size. Also, to speed up calculations, it is important to find out how large the time-step size can be in DDPM without affecting the predictions at least in the case of highly loaded cyclones. Therefore, a sensitivity analysis is done covering 1-100 ms of time steps.
In this study, the time step is controlled by the courant number to make the analysis less dependent on the case and meshing. As shown in Fig. 20, in the range of 1-12 ms (courant number 25-300) the predictions are not notably sensitive to the time-step size (deviations from the base case are below 2%), while for larger time steps, the predictions especially the separation efficiency start to deviate. Therefore, it can be concluded that 12 ms (courant 300) is the largest time-steps providing time-step size independent results so it is used for all other simulations in the present work as it is the least computationally expensive.
The present study confirms that unlike other conventional multiphase models such as TFM which requires 0.1-1 ms time-steps to have a converging solution [32], DDPM can be run with greater time-steps without affecting the predictions which makes this model an interesting option for the simulation of large-scale cases [33,34]. Although the size of the time-step is case-dependent, the present study suggests a range (around a courant number of 300) as a starting point for future studies in other cases.

Parcel (grain) size
In the present study, the strategy of coarse-graining was to have the same parcel size (same parcel mass) for different particle sizes. In DDPM, to ensure convergence of gas-phase calculations, it is required to have parcel (grain) size finer than cell size to avoid zero gas volume fraction in conservation equations that leads to solution divergence. On the other hand, the largest size of the grain is desired, to have fewer grains in the domain to track and have faster particle tracking. It was found that the solver can handle having a minor fraction of cells smaller than the parcel size without causing a convergence issue for the solver, but the limit of this fraction is case-dependent and is not clear. Therefore, a sensitivity analysis on parcel size is done to find this limit for the present case and the parcel size that provides the fastest solution predicting parameters of interest similar to when finer parcels are used.
Seven parcel masses were tested and compared in Fig. 21 showing that up to a parcel mass of 80 mg the predictions are not sensitive to the parcel size while larger parcels provide notably diffident predictions. It can be observed that in the case of 320 mg parcels, the convergence is lost in a few time-steps (even with 100 iterations in each time-step), and there are more of these time steps (with no convergence) in cases with larger parcels. Therefore, it can be concluded that as far as convergence in each time-steps is achieved, the predictions are independent of parcel size.
In the case of 80 mg parcel mass (the largest tested parcel with no convergence issue and providing parcel size independent predictions) only 10 k cells out of 4.2 M cells (0.2%) are smaller than parcel size (80 mg), while in the case of 320 mg parcel mass (finest tested parcel with convergence issue) 36 k cells (0.9%) are smaller than parcel size (320 mg). Therefore, it seems that having up to 0.2% of cells smaller than the parcel size is the limit and can be handled by the solver.
It should be noted that although larger parcel size leads to faster particle tracking, it does not always speed up the whole calculations as  Fig. 22. Effect of the discretization method on the overall pressure drop and separation efficiency. Arrows and corresponding percentages show the difference between CFD simulation and experimental measurement.
when larger parcels are used, in some cases more iterations in Eulerian calculations will be needed to converge the solution. In the present study, among cases that provide parcel size independent predictions (10, 20, 40, 80 mg parcel) the one with 40 mg provides the fastest solution while it has more parcels in the domain to track compared to the case of 80 mg parcel. Therefore, a parcel size of 40 mg is chosen as the base case to be used for all other simulations. The present study confirms that DDPM (as it implicitly approximates particle-particle interaction forces) can be run with a larger parcel size compared to DEM (which directly resolves these forces). In the specific case of the present work which turned out that particle-particle interaction forces have a minor effect, the parcel size can be enlarged without affecting the predictions to the limit that solution convergence is achieved. However, this conclusion is limited to the present case, as in other DDPM CFD simulations where particle-particle interaction forces might be important, smaller parcels might be needed.

Discretization method
Another way to speed up calculations is to use lower-order discretization methods to have the convergence of the solution with fewer iterations. However, as the lower-order discretization means lower accuracy, the significance of its impact should be investigated. In the CFD solvers using finite volume formulation such as Ansys Fluent, there are different discretization schemes available for interpolation of each variable from the cell center to control volume surfaces. In the present work, two different sets of discretization schemes sumerized in Table 4 (one referred to as 1st order as most of the variables are discretized by First Order schemes and the other referred to as 2nd order as most of the variables are discretized by Second Order scheme) are tested and compared in Fig. 22.
As it can be observed, the impact is significant; unlike quite accurate predictions of the 2nd order scheme (base case), the 1st order scheme highly underpredicts both pressure drop and separation efficiency by 27% and 10% respectively. Therefore, although 2nd order scheme is more computationally expensive (it requires almost 1.5 times the number of iterations to converge), it is necessary to have accurate enough predictions.

Particle size distribution
It is known that the performance (pressure drop and separation efficiency) of a cyclone is directly dependent on the PSD of solid material. Although the measurements were conducted with one specific PSD (D 32 = 6 μm, base case), to see the effect of PSD on the performance parameters of a highly loaded cyclone, five other PSDs (D 32 = 1, 2, 6, 16, 39, and 98) presented in Fig. 23 were tested in CFD simulations and the predictions were compared in Fig. 24.
As expected, the higher the mean particle size, the higher the separation efficiency and the lower the pressure drop (as larger particles require less energy to be carried by gas). It can be observed that the predictions are highly dependent on the PSD, so it is important to measure the PSD of the injected particle as accurate as possible.
In Fig. 25(a), for each of these tested mean particle sizes, the PSD of injected particles are compared with the PSD of particles leaving the domain (after agglomeration). For finer particles (e.g., D 32 = 1 μm) greater differences between PSD of injected and leaving particles are  observed showing that the effect of agglomeration is more significant when finer particles are operated. While for larger particles (e.g., D 32 = 98 μm) the difference between the PSD of injected and leaving particles becomes minor compared to cases with finer particles as the agglomeration becomes less crucial.
In addition, grade efficiency curves obtained for each of these tested mean particle sizes are compared in Fig. 25(b). The larger the mean particle size, the harder to observe minima (fish-hook behavior) in the grade efficiency curve (no minimum observed in case of D 32 = 39 and 98 μm) as such behavior is due to agglomeration and with larger particles the impact of agglomeration becomes minor. Note that in all cases that the minima were observed it was located in the range of 1-10 μm.
Another point to mention is that even in the case of D 32 = 1, the grade efficiency curve does not go below 50% showing the high performance of the studied cyclone.

Solid load
In addition to PSD, another operating condition affecting the performance of cyclones is solid loading. The experiments in the present study were done in the range of 0.3-0.6 kg solid/kg air, but the range was extended in the CFD simulations to lower solid loadings (down to 0.001 kg solid/kg air) and the predictions were compared with wellknown analytical models valid in the range. Also, to highlight the significance of agglomeration and the importance of including it in the calculations, the predictions of CFD simulation and analytical models are compared with the predictions of CFD simulation excluding the agglomeration modeling.
The analytical model used for the prediction of separation efficiency in the cyclone is developed by Muschelknautz and co-workers [81][82][83] that uses an approach known as an equilibrium balance. The geometrical and operational (including solid load) features of the cyclone are both considered in this model. On the other hand, for the pressure drop a universal model accounting for four types of major losses in the cyclone (expansion loss at the cyclone inlet, contraction loss at the entrance of the outlet tube, swirling loss and dissipation loss of the gas dynamic energy in the outlet) proposed by Chen [84] that considers the effect of solid loading. It should be noted that, although the Chen model is reported to highly over-predict the pressure drop [7](the over-prediction is also observed in the present work), the trends are reliable.
The effect of solid loading in the range of 0.001-0.65 kg solid/kg air on the prediction of separation efficiency, overall (cyclone+riser) pressure drop, and cyclone pressure drop are presented in Fig. 26. As expected, in the whole range, with an increase in the solid load the separation efficiency is improved and the pressure drop is reduced which are observed both in the analytical models and the CFD simulation only when the agglomeration model is included. The analytical models suggest that the separation efficiency and pressure drop are more sensitive to the solid load at lower loadings, and the CFD model (with agglomeration) is capturing the same.
Another point to mention is that at low solid loads the deviation of CFD predictions with and without agglomeration is minor compared to higher loads (deviation increases with an increase in the solid load) suggesting that agglomeration plays a more important role at high solid loads as expected.

Gas underflow
In the present experimental setup, the dipleg was connected to a dust bin, which is supposed to have no leakage. However, in the case of industrial cyclones, in the usual cascade arrangement or cyclones in CFB units a minor amount of gas might leave the cyclone with separated particles (a few percentages, below 5%, of inlet gas flow rate flowing downward in the dipleg) [7]. It is worth mentioning that in some industrial applications the opposite might also happen (some gas might come to the cyclone from the dipleg due to the pressure difference) that is not studied here. To see how this underflow would affect the performance parameters of a highly loaded cyclone, a simulation case assuming 5% gas underflow was done and compared with a case with no gas underflow as the base case in the present study (see Fig. 27). The comparison shows that the impact of gas underflow is relatively notable as it affects the pressure drop and separation efficiency predictions by 3.9% and 3% respectively; gas underflow leads to higher separation efficiency and lower pressure drop.

Conclusions
In the present study, sensitivity analyses of sub-models, model parameters, and numerical parameters on a hybrid multiphase model coupled with an agglomeration model developed for the simulation of industrial highly loaded cyclones were performed. Different options for particle-particle interaction modeling (based on KTGF), several turbulence models, and several drag models (five homogeneous drag models and three heterogeneous drag models) were examined. Additionally, the particle-particle restitution coefficient (that affects particle agglomeration), and particle-wall rebound coefficient were studied. Several numerical parameters/schemes were investigated in this work, including the discretization method, grain-to-particle size ratio, grid-to-parcel size ratio, and time-step size. In addition, the effect of operating conditions including solid loading, PSD, and gas underflow were examined. We believe that the results of this study provide some guidelines that are valuable for future CFD simulations on large-scale cyclones with high solid loadings.
The following points are drawn from the present study as the main conclusions: • Regarding sub-models, turbulence modeling is the most crucial one.
Including a curvature correction term to sensitize the turbulence model for rotation is vital. The sst version of the k-ω turbulence model provides the most reliable predictions that captures the important trend of improvement in separation efficiency when the solid load is increased. The realizable and standard version of the k-ε model and the standard version of the k-ω model are not able to predict the separation efficiency accurately. Although the RNG version of the k-ε model provides more or less reasonable predictions, it is unable to capture the trend of improvement in separation efficiency when the solid load is increased. • Particle-particle interaction force modeling is not crucial and can be excluded to speed up calculations • Particle rotation, Saffman force, Magnus force, and particle turbulent dispersion can be excluded from the calculation as their impact is minor • Although the effect of drag modeling is minor compared to turbulence modeling or particle-particle restitution coefficient, Igci and Sarkar drag models provided the best predictions. • Regarding model parameters, the particle-particle restitution coefficient (directly affecting the extent of agglomeration) is the most crucial one. • Although the particle-particle restitution coefficient affects both separation efficiency and pressure drop in the whole range, only the separation efficiency is sensitive to the particle-wall rebound coefficient just in the range of 0.8-1. • The time-step corresponding to a courant number of 300 is the largest (the fasted) time-step providing independent results. • Parcel size can be enlarged to the limit that the solution has no convergence issue without affecting the predictions. It seems that (at least in the present case and mesh) having 0.2% of cells being smaller than the parcel size is fine. • 2nd order discretization scheme is essential to have sufficiently accurate predictions. • The finer the mean particle size and the higher the solid load, the more significant the effect of particle agglomeration. • Considering that different combinations of drag models and particleparticle restitution coefficients provide reasonable predictions of pressure drop and separation efficiency, LDA measurement providing flow patterns in a highly loaded cyclone is needed for further evaluation of drag models.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.

Acknowledgement
This work was supported by the ProBu project: Process Technology for Sustainable Building Materials Production (grant number: 8055- Fig. 27. Effect of gas underflow on the overall pressure drop and separation efficiency. Arrows and corresponding percentages show the difference between CFD simulation and experimental measurement.