A hybrid multiphase model accounting for particle agglomeration for coarse-grid simulation of dense solid flow inside large-scale cyclones

A hybrid multiphase model (Dense Discrete Phase Model, DDPM) coupled with an agglomeration model and a sub-grid drag model was developed for simulation of industrial-scale cyclones with high solid loading. The model is validated by experimental results of pressure drop and separation ef ﬁ ciency from a pilot-scale cyclone withadiameterof1.6m.Keytrendssuchasimprovementinseparationef ﬁ ciencyandreductioninpressuredrop of cyclone due to an increase in particle load are well captured by the model. It is concluded that including the agglomeration model is crucial in particular in cases involving very ﬁ ne particles (d < 15 μ m) for accurate predictions of pressure drop and separation ef ﬁ ciency, while using the sub-grid drag modi ﬁ cation improves the prediction of separation ef ﬁ ciency.


Introduction
Industrial cyclones operating under high solid loads (greater than 0.1 kg solid/kg air) are widely used as a key component in building material industries such as cement production and also in processes employing Circulating Fluidized Beds (CFB). However, there is still a lack of consensus on the key mechanisms involved in operation of highly-loaded cyclone separators, considering the limited number of experimental and modeling studies focusing on them. While conducting informative experiments might be unfeasible on industrial cyclones, a reliable CFD model may provide better knowledge on the mechanisms occurring in such industrial cyclones and thereby helping for better scaling-up, design, optimization and troubleshooting.
Based on the framework by which the dispersed phase is treated, four-way coupling multiphase modeling of dense dispersed particleladen flows is conventionally divided into two types: Discrete Element Method (DEM) [38], and Eulerian-Eulerian or granular Two-Fluid Model (TFM) [39]. In addition to these two, a more novel type, known as the hybrid multiphase model [24], was developed to benefit from the advantages of two conventional types and avoid their disadvantages. All these three types have been used for the simulation of dense particle-laden flow in cyclones [18][19][20][21][22][23][24][25][26][27][28][30][31][32][33][34]36,40,41]. A version of the hybrid multiphase model introduced by Ansys Fluent, known as the Dense Discrete Phase Model (DDPM), seems to be a good option for simulation of industrial-scale highly loaded cyclones. DDPM is computationally affordable for industrial dense particle-laden flows as it models particle-particle interactions unlike DEM that directly resolves them. In addition, unlike TFM, DDPM easily handle polydispersed particles which is crucial for the hydrodynamics and the performance (separation efficiency) of cyclones. However, DDPM requires some modifications for simulation of highly loaded industrial cyclones working with fine particles.
Accurate prediction of drag force as the dominant force that governs particle motion in cyclones is vital if reasonable predictions of performance parameters of cyclones are desired [24]. It is widely reported that when multiphase models based on Kinetic Theory of Granular Flow (KTGF) (and the DDPM is one of them) are used with coarse grids, then using common homogeneous drag models leads to failure in the prediction of the multiphase flow hydrodynamics [42,43]. This failure is due to the unrealistic prediction of drag force caused by the inability of coarse grids to resolve mesoscale structures (such as particle strands in cyclones) [24,[42][43][44][45][46][47][48][49][50]. To overcome the deficiencies of homogeneous drag models for industrial-scale simulations heterogeneous drag models (also called sub-grid drag models) should be used to account for the effect of small mesoscale and heterogeneous structures on the drag correlations [24,[42][43][44][45][46][47][48][49][50].
The size distribution of particles in industrial cyclones often includes very fine (even sub-micron) particles. When a particle-laden system includes very fine particles (usually finer than 10 μm) particle agglomeration is an important mechanism [51,52] especially in cyclones that the separation efficiency is highly sensitive towards the particle size. In fact, through agglomeration, larger particles (agglomerates) are formed from fines, which easier separates from the gas and thereby improves the separation efficiency of the cyclone [24,[53][54][55]. In addition, the improvement in separation efficiency with increased solid load may be adequately explained by the agglomeration phenomenon in cyclones. This explanation proposed by Mothes et al. [56] and supported by others [53,54] elaborates that higher particle load leads to more particle collisions and consequently more agglomeration and larger agglomerates. This phenomenon is believed to be significant enough to cause the observed increase in the overall separation efficiency in the cyclone [57].
Another well-known observation called "fish-hook" behavior; the observation of minima in the grade efficiency curve of the cyclone at an intermediate particle size (meaning that finer particles are more easily separated compared to intermediate size particles) can be explained by agglomeration phenomenon. Finer particles are more prone to agglomeration compared to intermediate size particles, so they make more and larger agglomerates and have high separation efficiency. The observation of such minima in many experimental studies [54,55,[57][58][59][60] indicates the significance of the agglomeration phenomenon.
It is worth mentioning that the observation of such minima in the grade efficiency curve is usually interpreted as evidence of agglomeration, however occurring of agglomeration in a cyclone does not guarantee having a minima in the curve [54,55]. The presence or absence of a minima and its location (particle size and fractional efficiency where minima is observed) depends on the operating conditions (most importantly gas and solid flow rate), the geometry of the cyclone, and the properties and size distribution of particles; therefore, such a minima might not be observed in all cyclone cases even though agglomeration is occurring in them [54,55].
Considering the aforementioned reasons, if a reliable prediction of separation efficiency in the simulation of cyclones is required, including a model to account for particle agglomeration is necessary. This becomes more important in cyclones operating with very fine particles as such particles are more likely to agglomerate. Also, in cyclones operating under high solid loads where particle collision is more frequent, particle agglomeration plays a more important role.
Most CFD models of cyclones presented in the literature are adequate only for low and moderate loaded and lab-scale cyclones as they neglect particle-particle interaction and particle agglomeration and use quite fine mesh which is inapplicable in industrial cases. A few models focused on highly loaded cyclones mostly considering large particles [28,33,41] where agglomeration does not play a significant role, and they were mostly implemented on lab-scale cases [18,26,28,33] with fine grids where homogeneous drag models are adequate. In addition there is lack of validation of models in prediction of separation efficiency and the effect of solid load on the performance of cyclones.
In order to provide a better understanding of phenomena occurring inside highly loaded cyclones operating under industrial conditions, a pilot-scale cyclone running with limestone particles is simulated in this work. A hybrid Eulerian-Lagrangian multiphase model (DDPM) combined with a stochastic model to account for particle agglomeration is used with a sub-grid drag model to account for the effect of small structures not resolved by a coarse grid. The simulations are carried out using ANSYS FLUENT 2019R3 and validated against measurements of separation efficiency and pressure drop conducted under several operating conditions of the pilot-scale cyclone and thereby evaluating the reliability, and predictive capability of the proposed model. Fig. 1 shows the configuration of a pilot-scale cyclone (1.6 m in diameter) located in FLSmidth R&D Center Dania tested under industrial operating conditions (particle load, material, and size, and hydrodynamic conditions are as the industrial cyclones but the setup works at room temperature) to provide validation data for the present CFD study. The configuration is made up of a cyclone-like swirl maker (2), the cyclone (5) with a dust bin (7) and the outlet pipe, a riser (4) connecting the swirl maker to the cyclone, and a particle injection pipe (3). The swirl maker is designed to mimic the actual flow pattern (swirling motion in the riser) in industrial configurations of cascades of cyclones. The gas stream is introduced to the swirl maker to have a swirling motion, then the particles are injected in the riser where the upward swirling motion of gas carries all the particles to the cyclone.

Experimental study for model validation
The particles (mineral raw meal used for cement production mostly comprised of limestone and clay) used in the experiments had a wide PSD (Particles Size Distribution) similar to the PSD used in cement plants. It is seen in Fig. 2 that most particles are between 1 and 100 μm). The particles had a sauter mean diameter of 5.0 μm with a density of 3100 kg/m 3 , which is classified as Geldart C type of particle (based on Geldart's particle classification [61]); and known as cohesive particles which are very prone to agglomeration.
Static pressure was measured at locations (red dots) illustrated in Fig. 1 to measure the pressure drop over the riser and cyclone. The separation efficiency was also measured based on the variation of accumulated mass in the dust bin, and the variation of mass in the silo connected to the particle injection pipe. Both parameters were measured and recorded for 5-7.5 min to ensure the stability of operation.
The experiments were conducted at various operating conditions with different gas flow rates and solid flow rates (covering overall volume fraction of 0.0001 to 0.0003 as the common load in industrial cement cyclone preheaters, corresponding to solid load ratio from 0.34 to 0.66 kg solid/ kg air). Seven of these operating conditions are used in the present study as validation data. In Table 1, the operating conditions of these seven cases in addition to measured separation efficiency and pressure drop are summarized.

Model description
The hybrid Eulerian-Lagrangian multiphase model used in the present study (DDPM) is a Lagrangian based hybrid model that implicitly approximates particle-particle interaction using the Kinetic Theory of Granular Flow (KTGF). In DDPM, 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 Kinetic Theory of Granular Flow (KTGF) [62]. This approach was recently used for the study of cyclones working under high solid loads [25][26][27][28]41]. The equations describing this model including continuous phase and dispersed phase equations can be found in the literature [28]. To evaluate solids stress tensor required in DDPM for modeling of particleparticle interactions, the model uses several closure law equations based on the KTGF to estimate some required parameters such as bulk, collisional, kinetic and frictional viscosities. Those used in the present study (summarized in Table S-1 in the Supplementary Materials) were the same as other studies [26,28,41] in the literature using DDPM to simulate cyclones.
Regarding turbulence modeling, in addition to Reynolds Stress Model (RSM) [1,8,33,35,[63][64][65][66] and Large Eddy Simulation (LES) [6,8,11,13,20,63,[67][68][69][70] as widely accepted turbulence models in the literature for simulation of cyclones, a novel Eddy Viscosity Models (EVMs) sensitized to flow rotation and curvature such as the one available in Ansys Fluent (k − ω SST − CC) have attracted attention recently. Such sensitization can overcome weaknesses of EVMs in the simulation of cyclones providing quite reliable results, competitive with more complicated turbulence models [71][72][73][74][75]. LES requires very fine girds that is unfeasible in industrial-scale cases. On the other hand, EVMs are way more computationally efficient and robust compared to RSM that is crucial in simulation of industrial cases. These make k − ω SST − CC an optimal option for simulation of turbulence in cyclones. The equations describing the turbulence model used in the present work (a dispersed version of k − ω sst turbulence model sensitized to rotation and curvature) can be found in Ansys Fluent theory guide [76].
In the present study, for particle tracking, in addition to translational particle motion, rotational particle motion including a rotational drag force [77] was considered to provide more realistic particle motion in a cyclone as suggested by several authors [19,20,33,70,78]. Also, the Magnus force was added to the translational particle motion by using Oesterle-Bui-Dinh correlation [79].
Particle-wall collision was modeled by normal and tangential reflection coefficients (i.e. the ratio between particle rebound velocity to particle impact velocity) of 0.8 as used by many others in cyclones. Moreover, considering the proven impact of wall roughness in particle motion in confined cases such as cyclone [57,80] a rough wall model [81][82][83] available in Ansys Fluent with roughness parameters of a wall with ISO no. 12 (ISO 4287) was used to obtain more realistic particle trajectories.

Drag model
Various drag models were designed and validated for different flow regimes and unit operations, but none of them was designed or validated in cyclones. However, a model known as the revised Sarkar et al. drag model [42] successfully verified and validated for a wide range of applications such as bubbling fluidized bed, turbulent fluidized bed, dense circulating fluidized bed, dilute circulating fluidized bed [42], dual-circulating fluidized bed [50], CFB risers [84,85], and FCC riser [86] seems to be a suitable option to be tested in cyclones as well.
None of the sub-grid drag models available in the literature is intended for the cyclones where the flow rotation is dominant. However, Sarkar et al. drag model [42] was chosen over others due to its wide range of verified applications, and also because it provided reliable predictions in pre-tests simulations. The revised Sarkar et al. drag model [42] provided predictions in slightly better agreement with experimental results when compared with other tested sub-grid drag models such as Pirker et al. [87], Radl et al. [47], and Igci et al. [45]. In particular, it better captured the trend of improvement in the separation efficiency when the solid load is increased.
Sub-grid drag models basically formulate a drag reduction resulting from sub-grid level structures such as clustering effects that leads to a more accurate prediction of drag force in coarse-gird simulation. Normally a heterogeneity index H d , accounting for the hydrodynamic disparity between homogeneous and heterogeneous drag is defined based on a specific homogeneous drag model. In the present model heterogeneity index is defined based on the Wen-Yu drag model in Eq. (1): In the revised Sarkar et al. drag model [42], the heterogeneity index H d is formulated as follows by Eqs. (2) to (11):   2)-(10)) [42]. Note that using sub-grid drag models constituted based on Eulerian-Eulerian simulations (such as the revised Sarkar et al. drag model) in an Eulerian-Lagrangian model (such as DDPM in present work) is an approximation as the particle information needs be mapped to the Eulerian grid before using the correlations [48,88]. An approximation of this kind is reasonable and the drag models deduced for Eulerian-Eulerian simulations can be used in Eulerian-Lagrangian simulations (such as DDPM and MP-PIC) and the other way around [48,88] (Particularly, the revised Sarkar et al. drag model has been used successfully in Eulerian-Lagrangian models in the literature [50,84]).
The filter-to-grid ratio of 2 (Δ f = 2(V cell ) 1/3 ) was suggested in the revised Sarkar model [42]. In general, when using sub-grid drag models that are a function of the filter size, the predictions might be affected by filter-to-grid ratio so that in those cases the ratio should to be optimized [89][90][91] or instead more novel models that make H d independent of filter size should be used [92]. However, the sensitivity analysis of the filter-to-grid ratio (the common values of 1, 2, and 3 were tested) conducted for case-2 of operating conditions of the present study showed that the predictions are not considerably affected by the filter-to-grid ratio (the variation of pressure drop and separation efficiency is less than 2% and 1% respectively), so the suggested value of 2 could be used.
The revised Sarkar et al. drag model was implemented to the Ansys Fluent solver using "DEFINE DPM DRAG" user-defined function. It should be noted that Wen-Yu drag model is also tested in the present work for comparison.  Table 2 List of equations for agglomeration modeling.
XX and ZZ are randomly generated numbers in the range [0, 1] so that L a < 1 (17) Constants are summarized in Table 3; values proposed by Sommerfeld [51] (20)

Agglomeration model
The agglomeration was described by a stochastic Lagrangian interparticle collision model proposed by Sommerfeld [51,93], assuming that two particles may collide only if they are located in the same Eulerian cell. The model is based on the algorithm of O'Rourke [94] to efficiently reduce the computational cost.
The calculation algorithm of the agglomeration model is summarized in Fig. 3 and details of the model equations are listed in Table 2. The probability of two parcels in a cell colliding N times in a time step Δt, P(N), can be estimated using eqs. 12 and 13. Considering that the probability of not having any collision at all is P 0 ð Þ ¼ exp −n ð Þ, to at least have one probable collision a randomly generated number RN in the range [0, 1] must be greater than P(0) (eq. 14). On the other hand, the collision efficiency, η p , defined as the ratio of the number of particles hitting the collector to the number which would strike it if the streamlines were not diverted by the collector can be estimated using eqs. 15 and 16. A collision is only occurring if the lateral displacement L (defined using a random process in eq. 17) is smaller than the value of Y c (the radial distance defined in eq. 18) plus the small particle radius (eq. 19). To have a collision both criteria in eqs. 14 and 19 must be met.
When a collision takes place, two collided particles either stick (and form an agglomerate) or rebound. Only van der Waals force was considered as the interparticle force. The possibility of particles sticking or rebound can be evaluated using a critical velocity based on an energy balance, which is given in eq. 20. When the normal relative velocity of two colliding particles is lower than this critical velocity, the restitution energy after the impact does not exceed the van der Waals energy between collided particles, and agglomeration will occur (eq. 21).
When two collided particles form an agglomerate, the size of the agglomerate and its velocity, was calculated based on mass and momentum conservation. To simplify, spherical agglomerates (with the mass of two collided particles) without any void (which is referred to as the sphere model [20]) was considered as the outcome of the agglomeration process. The assumption and derivation of the equations are described in detail in the literature [51,93].
This algorithm was implemented to the Ansys Fluent solver using "DEFINE DPM SPRAY COLLIDE" user-defined function and was called whenever particle positions were updated.

Rresults and discussion
The suggested model performs numerically stable and provides reasonable flow pattern, particle motion, and pressure drop and separation efficiency predictions that will be discussed in detail in the following sections.

Mesh independence
The mesh independence study was performed for two of the operating conditions (case-4 and case-6) using seven levels of mesh resolution (polyhedral mesh with prism layers described with more detail in Supplemental materials). Some information about these seven levels of mesh including average of cell volume, grid resolution, and grid to parcel size ratio are presented in Table 4. It is worth mentioning that considering the values of =d p being in order of 10 3 , the present study is considered as a coarse-grid simulation (normally when the grid size is significantly greater than 10 particle diameters, the simulation is referred as coarse-grid simulation [95]). The two cases were selected as two cases with the same gas flow rate (9.4 kg/s) but with different solid loads (case-4 with loading ratio of 0.3 kg solid/kg and case-6 with loading ratio of 0.6 kg solid/kg) to investigate whether grid size could influence trends of cyclone performance parameters when the solid load is increased, in particular the improvement in separation efficiency.
Predicted pressure drop and separation efficiency obtained from these levels are compared in Fig. 4 and Fig. 5. As it can be observed from these figures the predictions did not deviate considerably from the prediction of finest mesh, when 4.2 million cells or finer was used; The deviations were less than 1% in case of separation efficiency and less than 2% in case of pressure drop, suggesting that predicted results were independent of the characteristics of the mesh size if at least 4.2 million computational cells were used. In addition, the improvement in separation efficiency when the solid load was increased from case-4 to case-6 as a well-known trend also observed in our measurements could just be captured if at least 4.2 million cells were used while coarser girds were unable to capture such an improvement. Therefore, the rest of the simulations and studies were done using this mesh resolution.
Also, to quantitatively evaluate the grid convergence and the error due to grid resolution [96], convergence calculations for the three finest grid levels using the approach suggested by Roache [97,98] are presented in Table 5. As the values for α are close to unity, it can be concluded that the solution is within the asymptotic range of convergence. Using these calculations, the values when the grid spacing vanishes can be estimated based on the Richardson extrapolation [97] with an error band (for example the pressure drop in case 4 is estimated to be 1211 Pa with an error band of 2.979%).
Furthermore, the comparison of flow pattern distributions in seven levels of mesh presented in Fig. 6 shows minor differences in flow patterns obtained by different grids. These differences were mainly in the core region and in quantifying the peak of velocities. However, in the near-wall region, the region with the greatest impact on the performance of a cyclone, very similar patterns were obtained by different grids. Variations in flow patterns explain the minor variation in predicted pressure drop and separation efficiency. Table 3 Particle properties for calculation of the critical impact velocity.

Parameter Value
Hamaker constant H (J) 5 × 10 −19 Minimum contact distance z 0 (m) 4 × 10 −10 Yield pressure p c (Pa) 5 × 10 9 Energetic restitution coefficient e 0.7 Table 4 Levels of grids in mesh independence study. dp is the sauter mean diameter of particles and equals 5 × 10 −6 m. Level 3 is used for the rest of simulations. In the present study, the strategy of coarse-graining was to have same parcel size (same parcel mass) for different particle sizes. Three parcel sizes (1, 2 and 4 micrograms) were tested to perform parcel size independence. The deviations of results were negligible (below 0.5% for both pressure drop and separation efficiency).
In DDPM, the use of relatively large time steps (up to 0.1 s) is reported [28,99]   Simulations were carried out on a Linux cluster using 32 processors (Intel(R) Xeon(R) Gold 6226R @ 2.90 GHz and 64 GB RAM with highspeed Infiniband network connection). For chosen mesh, parcel size, and time step, it took less than 2.5 h to have 1 s progress in flow-time (meaning that it will not take longer than 3 days to ensure steady results; tracking parameters show that 20-30 s of flow-time guarantee the steady condition, however, simulations continued up to 50 s). Therefore, considering available computational resources, the presented model is easily applicable and affordable for industrial-scale cyclone simulation.

Importance of agglomeration model
To show the significance of including an agglomeration model in the simulation of a highly-loaded cyclone, three of the cases were also simulated without implementing the proposed agglomeration model. The results are presented in Fig. 7 and Fig. 8 and are compared with experimental measurements and also model predictions where the agglomeration model was included. As shown in Fig. 8 if an agglomeration model is not included, the separation efficiency will be highly under-predicted as very fine particles with low inertia remain in the domain which mostly escape from the cyclone vortex finder. Without agglomeration, particles enter the cyclone with a lower mean size ending up with a lower separation efficiency compared to the case with agglomeration.
Without an agglomeration model the improvement in separation efficiency due to an increase in solid load cannot be captured (comparing case-2 and case-4). An increase in solid load dampens swirling motion and consequently lower the centrifugal force acting on particles while PSD is not altered as agglomeration is not included, so separation efficiency predictions will be lower when the solid load is increased.
In addition to separation efficiency predictions, the prediction of pressure drop is also erroneous when agglomeration is neglected. Without agglomeration, the gas phase must carry particles finer than the actual size of particles (that experience agglomeration) which leads to more energy loss as finer particles need more energy to be carried as they experience a greater drag force. While measurements and simulation with the agglomeration model did not report a significant change in Table 5 Grid convergence calculations using GCI method [97,98] for the three finest grid levels. i = 0 is the value at zero grid space.  It is also worth mentioning that excluding the agglomeration model, led to unphysical particle motion so that accumulation of particles in the cyclone was be observed.
Overall, as it can be observed in Fig. 7 and Fig. 8 to accurately predict separation efficiency and pressure drop, and to observe important trends such as improvement in separation efficiency when the load is increased, including agglomeration model is crucial.

Importance of sub-grid drag model
Studying the importance of using a heterogamous drag model instead of a homogeneous one was done by simulating three cases once with the presented sub-grid drag model (the revised Sarkar drag model [42]) and once using a homogeneous drag model (the Wen-Yu drag model [100]). The results are compared in Fig. 7 and Fig. 8. It should be noted that in both cases the agglomeration model is included. Although apparently the overall pressure drop predictions were not significantly affected by using either drag models (see Fig. 7), when the  Wen-Yu drag model is used the separation efficiency was over predicted (greater than 96.5% in all cases) due to over prediction of drag force as shown in Fig. 8. Using Wen-Yu drag model led to deviation of predicted separation efficiency from the measured value up to 6%, while by using the revised Sarkar drag model this deviation did not exceed 3%. In addition to this over prediction, the Wen-Yu drag model was also inferior to the revised Sarkar model in capturing the trend of improvement in separation efficiency due to an increase in solid load (see Fig. 8); the measurements reported a 4.6% increase in separation efficiency when the solid load is increased from 0.3 kg solid/kg air in case-4 to 0.6 kg solid/kg air in case-6, while the revised Sarkar model and Wen-Yu drag model predicted an 1.8% and an 0.8% increase in separation efficiency respectively. Although none of the models captured the exact trend, the revised Sarkar model performed better.
Therefore, although pressure drop prediction can be achieved with the same level of accuracy by using the homogeneous Wen-Yu drag model, for accurate prediction of separation efficiency, and consequently better-observing of improvement in separation efficiency when the load is increased, including a sub-grid drag model (the revised Sarkar drag model) is helpful.
In Fig. 9, a snapshot of a converged solution of case-2 showing particles colored with heterogeneity index H d ¼ β Sarkar β Wen−Yu is presented; the lower value of heterogeneity index, the higher over prediction of the drag force by the Wen-Yu drag model compared to a more realistic prediction of the revised Sarkar drag model. Having many particles with heterogeneity index values smaller than one shows an extent of drag over prediction by the Wen-Yu drag model and this illustrate the importance of including a sub-grid drag model. Simulations and results presented in the following sections were obtained by including agglomeration model and the revised the revised Sarkar drag model. It should be noted that these models were also included in the mesh independence study.

Pressure drop
In Fig. 10, pressure drop obtained from simulations are compared with the measurements for all seven cases; the model predicts the overall pressure drops in a reasonable agreement with experimental measurements, as the predicted values did not deviate more than 12% from the measurements.
Experimental measurements showed that the overall (riser + cyclone) pressure drop is more sensitive towards gas flow rate than solid flow rate, and the same was observed in the simulations. Considering case1-3 that have identical gas flow rates, varying the solid load did not affect the overall pressure drop significantly; measurements showed a slight reduction as the solid load was increased while the simulations predicted a slight increase. The same trend was observed in case4-6. The reason that the overall pressure drop is not affected significantly by the solid load is that the pressure drop in the riser and in the cyclone change differently with increasing the solid load so that they are cancelled out to some level. As the solid load is increased, the pressure drop of the cyclone is reduced due to the dampening of swirling motion (the way this dampening leads to a reduction in pressure drop which seems counter-intuitive is discussed in details in the literature [57]), while the pressure drop in the riser is increased as the gas needs to carry more particles leading to more loss of energy. The model can picture these well-known increase of pressure drop in the riser and reduction of pressure drop in the cyclone due to an increase in the solid load (see calculated pressure drops in the riser and the cyclone in case1-3 and case4-6 in Fig. 10).
Overall, considering all cases, as observed by both the measurements and the simulations, an increase in gas flow rate leads to greater pressure drop. To see just the effect of gas flow rate case-2 and case-6 that have an identical solid load (0.6 kg solid/kg air) can be compared. An increase in the gas flow rate from 8.2 m 3 /s in case-2 to 9.4 m 3 /s in case-6 led to a considerable increase in the overall and also both the riser and the cyclone pressure drop; Higher gas flow rate leads to more dissipation of flow energy due to stronger friction on the walls just like in normal pipe flow.

Separation efficiency
Overall separation efficiencies achieved from simulations are compared with the measurements for applied operating conditions in Fig. 11. Clearly observed from Fig. 11, the model estimates the overall separation efficiency quite accurately as the predicted values matched well with experimental measurements with small deviations (the deviations did not exceed 3%).
Unlike pressure drop, the separation efficiency seems to be more sensitive towards solid flow rate than gas flow rate. Such an observation was seen both in experimental measurements and simulations when comparing different operating conditions. Measurements showed that with a constant gas flow rate (case1-3 and case4-6), an increase in solid load leads to a measurable improvement in separations efficiency; such improvement is also widely reported in the literature. The presented model is well able to predict such improvement with comparable trends as can be observed in Fig. 11.
An increase in the solid load dampens the swirling motion in the cyclone so that particles experience less centrifugal force towards the wall and consequently are more likely to escape from the vortex finder to the outlet. However, the impact of the increase in particle agglomeration due to an increase in particle load is apparently significant enough to compensate for this dampening. Higher particle load leads to more particle collisions and consequently, more agglomeration and overall larger agglomerates that experience greater centrifugal force compared to overall smaller agglomerates in lower load although the swirling velocity is dampened. This will lead to improvement in the overall separation efficiency despite dampening of swirling motion. In addition, an increase in the solid load attenuates the turbulence (turbulent kinetic energy is reduced) that might improve the separation of particles as particles will be less dispersed by the turbulence. To see such a trend, accurate simulation of flow pattern and accurate modeling of agglomeration phenomenon are required, and the current model is accurate enough to do so.  Comparing case-2 with case-6 that have a similar solid load (0.6 kg solid/kg air) shows that increasing the gas flow rate while the solid load is maintained constant (by increasing solid flow rate) will not lead to a measurable change in separation efficiency of the cyclone. This means that the separation efficiency is more governed by solid load than gas flow rate. Such a trend was observed both in simulations and measurements. Fig. 12 presenting contours of gas tangential velocity, axial velocity, pressure and solid volume fraction averaged over 5 s after convergence for one of operating conditions (case-4) shows typical flow characteristics expected in cyclones; the rotational features of the cyclone flow ( Fig. 12(a)), ascending flow in the inner vortex and descending in the outer vortex ( Fig. 12(b)), drop in the pressure in the direction from the wall to the center (Fig. 12(c)), and concentration of particles in near wall region and formation of particle strand (Fig. 12(d)).

Flow pattern
To provide a picture of flow pattern in the cyclone and how it is affected by operating conditions, contours of tangential velocity and axial velocities, are presented in Fig. 13 and Fig. 14 respectively for three operating conditions; Case-2 (G8.2-S4.7), (b) Case-4 (G9.4-S3.2), and (c) Case-6 (G9.4-S5.8) (case-2 and case-6 have same solid loading ratio but different gas flow rate, while case-4 and case-6 have same gas flow rate but different solid loading ratio). In addition, radial profiles at a specific height (2 m below cyclone roof, the location is shown with yellow line in Fig. 12) are provided to better show the variation of gas flow field inside the cyclone. Although local measurements were not conducted for the studied pilot-scale cyclone, the provided flow pattern is in qualitative agreement with what is expected from cyclones.
The expected swirling motion in cyclone can be observed in Fig. 13. The radial distribution of tangential velocity is well-matched with expectations; increase in the magnitude of tangential velocity as going far away from the center up to a region close to the wall where a sharp reduction in the magnitude of tangential velocity is observed. It can be observed that as expected an increase in the solid load from 0.3 kg solid/kg air in case-4 to 0.6 kg solid/kg air in case-6 while the gas flow rate is kept constant leads to a reduction in swirling motion of the gas phase. This dampening of swirling motion is more severe in the lower region of the cyclone (the apex) where the particles are concentrated (see Fig. 12(d)). In cyclones, it is generally believed that such reduction in swirling motion due to an increase in solid load will lead to a reduction in pressure drop, as observed in the present simulation (reduction in cyclone pressure drop from case-4 to case-6 in Fig. 10).
As shown in Fig. 12(c), high-pressure regions are observed close to the cyclone wall which reduce as moving towards the center, a typical feature of a swirling motion. A lower radial gradient of static pressure is observed when the solid load is increased due the dampening of swirling motion. Fig. 14 shows that the axial velocity profile is negative (descending) near the walls, meaning that the gas flow is swirling down (outer vortex) carrying particles concentrated in the near-wall region towards the dust bin as expected. On the other hand, far from the walls, around the center, the axial gas velocity is positive (ascending), meaning that the gas flow is swirling up (inner vortex) in the center as expected. It was observed that the width of the inner vortex is not sensitive to solid load or gas flow rate. As it can be observed in Fig. 14(d), the axial gas velocity increases in the direction from the wall to the center and reaches its maximum at the center. It can be observed that the high axial velocity region shrinks and move upwards as the solid load is increased. Fig. 12(d) shows that the particles are concentrated in the near-wall region as expected from cyclones; higher particle concentration in the direction from the center to the wall reaching the maximum on the wall. A concentration of particles was also observed at the conical end of the cyclone due to having a very thin dipleg.
In Fig. 15, magnitudes of tangential and axial velocity in the case of a particle-free cyclone and in the case of a highly-loaded cyclone (case-2, SL = 0.6 kg solid/kg air) are compared. In particle-free case, the magnitude of the tangential velocity is significantly greater than axial velocity, and therefore the fluctuations in tangential direction are considerably greater than the other directions. This means that in particle-free case the turbulence is anisotropic [101]. Therefore, in particle-free cases (and also in ordinary cyclone with very low solid load) assuming isotropic eddy-viscosity hypothesis (which is the base of EVMs turbulence model) might be erroneous and the turbulence model used in the present study might not be sufficient. However, as it can be observed in Fig. 15 (b), in the case of a highly loaded cyclone the magnitudes of tangential and axial velocity are comparable, so the fluctuations in different directions are comparable as well, therefore isotropic eddy-viscosity hypothesis assumed in k − ω sst turbulence model used in the present work is credible. In other words, although in case of cyclones with low solid load turbulence model leaving the isotropic eddy-viscosity hypothesis such as RSM or LES are required, in case of a highly loaded cyclone where the swilling motion is dampened so that velocity components are in the same order, using turbulence models assuming the isotropic eddy-viscosity hypothesis such as k − ω sst turbulence model used in the present work is sound and sufficient.

Grade efficiency
The grade efficiency curves obtained from CFD simulations for three operating conditions are present in Fig. 16. The "fish-hook" behavior is observed in all cases as there are at least one minima in each grade efficiency curve in the range of 1-10 μm. Such grade efficiencies (an overall increase of fractional efficiency as particle size become bigger with local minima) were reported in cyclones with high solid loading [102]. As expected and reported, an increase in solid load led to an upward shift in the grade efficiency curve. It should be noted that case-2 and case-6 that have identical solid loads are presenting similar grade efficiencies suggesting that the separation efficiency is more controlled by the solid load than gas flow rate in the present case and applied operating conditions. As it can be observed in Fig. 16, fractional efficiencies of particles finer than 15 μm were more sensitive to the solid load while the larger particles appeared to be less affected by the solid load level. Apparently, particles bigger than 50 μm were almost totally separated in the studied cyclone at the applied operating conditions. High values of fractional efficiency for sub-micron particles can be only justified by the effect of solid load and particle agglomeration as cyclones with low solid loads present very low values of fractional efficiencies for sub-micron particles (in some case 0% fractional   efficiency) since agglomeration does not play a role when the solid load is too low.
In addition to grade efficiency curves, to better see the impact of the formation of agglomerates, PSD of particles leaving the domain (after agglomeration) are presented for these three cases in Fig. 17 and are compared with PSD of injected particles. The significant difference between PSD of particles leaving the domain and PSD of particles injected into the domain shows the importance of the agglomeration phenomenon. As the solid load is increased, overall larger particles are formed as the curve is shifted to the right that is explained by the fact that higher load leads to more particle collision and agglomeration and larger clusters. Although around 30% of injected particles are finer than 6 μm, almost all particle leaving the domain were bigger than 6 μm as they all had formed agglomerates inside the domain. In case-2 and case-6 that have similar solid loads, the PSDs of particles leaving the domain were similar suggesting that in the studied case the agglomeration is mostly governed by the solid load.
In Fig. 18, PSD of particle in six cross sections of the domain are shown to see how the PSD vary in the device and where agglomeration takes place more frequently. The most significant change in the PSD between consecutive planes occurs between the injection plane and plane 1 meaning that the agglomeration mostly occurs in the injection pipe. Comparing the PSD in plane 1 and 2 shows that agglomeration is still occurring up to plane 2 as the PSD varies considerably. It seems that the agglomeration is almost complete by plane 2 and not occurring in the volume between plane 2 and 3 (the riser) as the PSDs in these planes are quite the same. From plane 3 to the plane 4 and to the particle exit plane the PSD change slightly which can be either due to further agglomeration of particles and that more fine particles escape to the vortex finder as going down in the cyclone.

Conclusions
In the present study, a hybrid multiphase model coupled with an agglomeration model based on a stochastic Lagrangian inter-particle collision and a sub-grid drag modification based on the revised Sarkar et al.
drag model [42] was developed for simulation of industrial cyclones with high particle load. The model was implemented in ANSYS Fluent to simulate a pilot-scale cyclone operating under industrial cement plant conditions. The suggested model performs numerically stable converging in a reasonable duration using available resources making it applicable and affordable for simulation of industrial-scale cyclones. The model provides reasonable flow pattern, particle motion, and pressure drop and separation efficiency predictions. This model was validated against experimental measurements (separation efficiency and pressure drop) under different gas and solid flow rates. The measurements were done on a pilot plant cyclone with an internal diameter of 1.6 m and solid loading ratio 0.34-0.66 kg solid/ kg air. Reasonable agreement with experimental pressure drops and separation efficiencies were achieved proving the capability of the model in capturing important phenomena occurring inside such a device.
The following points are drawn from the present study as the main conclusions: • In addition to accurate prediction of pressure drop and separation efficiency, the presented model can also capture the influence of changes in operation conditions observed in our experiments. The model well predicts the improvement in separation efficiency, reduction in the pressure drop of the cyclone body and increase in the pressure drop of the riser when the solid load is increased. The model is also capable of predicting the increase in the pressure drop of both cyclone and the riser when the gas flow rate is increased. • Including the agglomeration model for simulation of highly loaded cyclones is vital if reasonable predictions of overall and fractional separation efficiencies and also pressure drop are desired, in particular in cases including fine particles (finer than 15 μm). Without such an agglomeration model, the separation efficiency is highly under predicted while the pressure drop is highly overpredicted by the model. In addition, the important trend of improvement in separation efficiency when the solid load is increased cannot be captured unless the agglomeration model is implemented. With the agglomeration model inclusion, the well-known "fish-hook" behavior (a minima in grade efficiency curve) in cyclones can be captured by the CFD simulations. • The Wen-Yu homogenous drag model over predicts the separation efficiency when using coarse grids in a large-scale cyclone showing the necessity to use a sub-grid drag modification for the simulation of large-scale cyclones. However, the Wen-Yu drag model might be sufficient for the prediction of pressure drop in cyclones as it provides similar predictions as a heterogeneous drag model. • The applied sub-grid drag model called the revised Sarkar et al. drag model [42] that has previously been tested successfully for a wide range of applications is now also verified and validated for simulation of highly loaded industrial-scale cement plant cyclones. The model provides reliable predictions of important performance parameters of the setup in agreement with experimental measurements. We believe this model can be used with confidence for the prediction of drag in such a device.
Since there is not a specific sub-grid drag model intended for the cyclones in the literature, as future work, a sensitivity analysis on other drag models (both heterogeneous and homogeneous) and also sensitivity analysis on other models and parameters used in the present work can be conducted to investigate how they could influence predictions of performance parameters in cyclones.

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.