Impact of Contact Scaling and Drag Calculation on the Accuracy of Coarse-Grained Discrete Element Method

The accuracy of coarse-grained discrete element method (CGDEM) relies on appropriate scaling rules for contact and fluid-particle interaction forces. For fluidized bed applications, different scaling rules are used and compared with DEM results. The results indicated that in terms of averaged values as mean particle position and voidage profile, the coupling of computational fluid dynamics and CGDEM leads to accurate results for low scaling factors. Regarding the particle dynamics, the approach leads to an underestimation of RMS values of particle position indicating a loss of particle dynamics in the system due to coarse graining. The impact of cell cluster size on drag force calculation is studied. The use of energy minimization multiscale drag correction is investigated, and a reduced mesh dependency and good accuracy are observed.


Introduction
Fluidized particulate systems are widely used in the chemical and process industry. Common applications are heterogeneous catalysis, drying processes, polymerization, and coal combustion or gasification. Those applications are most often conducted in fluidized beds. This reactor type is characterized by physical phenomena acting on different length scales. The microscopic scale is determined by collisions between individual particles and between particles and the wall, as well as interactions between particles and the fluid phase. The microscale interactions significantly control energy, momentum, and mass transport and lead to local mesoscale phenomena like cluster and bubble formation and, in the case of nonspherical particles, a preferred particle orientation [1]. On the macroscopic level, those phenomena significantly affect mixing characteristics and flow regime [2,3].
Several numerical methods have been developed to describe this complex multiscale system, predictively. A numerical approach that describes the system based on microscale phenomena is the coupling of direct numerical simulation (DNS) to calculate the fluid dynamics and discrete element method (DEM) for the description of particle dynamics. Here, the interface between particles and fluid phase is spatially resolved using either the immersed-boundary method [4][5][6], interpolation techniques [7], or body-fitted meshes. The latter approach is cutting-edge technology for the particle-resolved simulation of fixed-bed reactors, where particle interactions are negligible [8]. The benefit of DNS compared to other numerical methods is that no additional closures are needed, e.g., the drag coeffi-cient is not an input parameter that is needed for the simulation but can be extracted from the simulation results. Although, recently, DNS was used for the simulation of fluidized beds [9][10][11][12], its applicability is limited since the method is numerically very demanding, and only a few thousand particles can be simulated today. Therefore, it is often used to develop closures for other models or to generate reference and benchmark cases for comparison against other not fully resolved methods [13].
Two well-established methods to describe the system on meso-scale level is the kinetic theory of granular flow (KTGF) and the coupling of computational fluid dynamics (CFD) with DEM (CFD-DEM). The former method describes the continuous and dispersed phase as separate Eulerian phases and only resolves scales bigger than individual particles. Interactions between fluid and particles, e.g., drag, lift, and turbulent dispersion, are incorporated via additional source terms in the momentum equation of each phase. For the modeling of particle-particle interactions, two new properties are introduced: the granular pressure and granular viscosity. Both properties are related to granular temperature [14]. Below a critical packing density, particle interactions are characterized by collisions between particles. The probability of particles to collide can be estimated by distribution functions that are used to calculate granular temperature, and with this, granular pressure and viscosity [15][16][17].
Above the critical packing density, particle dynamics are characterized by friction. Granular pressure and viscosity can be described by relations developed by Johnson and Jackson [18] or Schaeffer [19]. Overall, KTGF relies on multiple closures for the aforementioned properties and also the interphase momentum and possible heat and mass transfer. Its development is still a vital research topic [20][21][22][23]. Although the method is widely used [24][25][26][27][28], it has some constraints that limit its application on industrial-scale problems. To resolve all microand mesoscale phenomena, an adequate mesh size Δx 1) is needed. For Geldart Group A and B particles, a mesh size of Δx 10d p and Δx 16d p (d p meaning the particle diameter) was reported, respectively [29,30]. In many industrial applications, this constraint is not always satisfied, which leads to significant deviations regarding particle dynamics and bed height [31].
CFD-DEM is a numerical method that couples CFD and DEM, whereas the latter method was developed by [32]. Contrary to KTGF, only a small number of closures (see details below) are needed. The trajectory of each individual particle is calculated using Newton's law of motion. Particle-particle and particle-wall interactions are described by contact models [33,34]. Hence, these interactions are resolved on the microscale, while model-based closures are needed to describe fluidsolid interactions, e.g., drag force and possible heat and mass transfer [14]. Because of the high numerical effort to track all particles plus their contacts and long simulation times that are needed to get statistically independent results, the applicability of CFD-DEM is currently limited to simulations with few million particles [35]. To resolve micro-and mesoscale effects, a relatively fine mesh with a size of Δx 1.5d p is needed [36]. Therefore, not only the number of particles, but also the number of cells could be the bottleneck.
Although modern CFD and DEM codes are highly parallelized and can additionally take advantage of graphics processing unit (GPU) computing, this is not enough to apply KTGF or CFD-DEM on industrial-scale cases. A reduction of cell and/or particle count is needed, which leads to a reduced spatially resolution of the system (''coarsening''). In the framework of KTGF, a reduction of cell count leads to non-resolved mesoscale phenomena, e.g., cluster formation, which results in an overprediction of drag forces [29,37,38]. To increase accuracy of the drag calculation with the lower cell count, sub-grid models are needed to account for not resolved effects.
Since an accurate prediction of drag force is essential to describe fluidized particle systems, several research groups developed filtered drag models. A description and comparison of different approaches can be found in Schneiderbauer et al. [39]. The listed models predict significantly different results under the same conditions. The development of a general model is still an ongoing research topic. Nevertheless, it is obvious that the number of needed closures will further increase.
An alternative approach for the simulation of large-scale fluidized particle systems is the multiphase particle-in-cell method (MP-PIC). Like CFD-DEM, it is a combination of Eulerian and Lagrangian framework [40]. The flow field of the continuous phase is calculated by solving the Navier-Stokes equations. Additionally, a transport equation for a probability density function (PDF) is solved that describes particle dynamics and particle-particle interactions. The PDF is discretized by Lagrangian parcels that represent a specified number of original particles. By integrating the PDF, the phase fraction of the dispersed phase and the average particle density and velocity are determined. The values get mapped from the mesh to the Lagrangian parcels, which trajectories get updated.
To avoid unphysically high packing densities at high time step sizes, an additional model is used that applies additional stresses, either based on particle volume fraction only [41], or based on particles' coefficient of restitution, average particle density, granular temperature, and phase fraction of the dispersed phase [42]. The fundamental concept of this model requires a statistical representative number of particles within each cell; therefore, a grid size of Δx = 10d p -100d p is suggested. Inherent with the needed coarsening of cells, important flow features may not be resolved and might be compensated by sub-grid modeling of drag, as shown by [43,44]. The number of closures needed for MP-PIC is in the order of magnitude as for KTGF. For complex reactor configurations, including internals like heat exchanger, baffles, and nozzles, the mesh requirement cannot always be fulfilled, leading to instabilities and unphysically high packing densities.
A promising and relatively new modeling approach that can be applied for industrial-scale simulations is the coupling of CFD and coarse-grained DEM (CFD-CGDEM). The main concept is the same as for CFD-DEM, but not all individual particles are tracked. Instead, a specified number of particles gets lumped into a representative CG-parcel, as illustrated in Fig. 1, which reduced the numerical effort significantly. Trajectories

1)
List of symbols at the end of the paper.
for the CG-parcel and fluid-particle interactions are calculated as in CFD-DEM; however, the calculated interaction forces need to be scaled to assure the same conditions for the original and the coarse-grained system. In comparison to filtered-KTGF, the CFD-CGDEM approach has the benefit that the history of each CG-parcel is tracked, which is especially of interest if a heterogeneous chemical reaction on particle surfaces play a role. This is also true for MP-PIC, but in comparison CFD-CGDEM relies on a significantly lower number of needed closures. Recent studies have also shown advantages regarding the performance of CFD-CGDEM in comparison to MP-PIC [45]. A rapidly increasing number of researchers use CFD-CGDEM. The investigated problems include fluidized and non-fluidized particle systems [46][47][48][49][50][51][52][53][54]. Additional physical phenomena like heat transfer [48], adhesive interparticle forces [55,56], chemical reactions [57], and particle size distributions [46,54,58] have also been incorporated. Nevertheless, many different scaling rules for particle-particle and fluid-particle interaction forces have been reported in the literature, including the incorporation of sub-grid drag models. To date, only few comparisons between the different scaling rules were made, and no statement is given if system-dependent constraints that limit the maximum applicable coarse-graining ratio l 3 (number of particles per CG-parcel) exist.
Based on the simulation of an uniaxial compression test of cohesive particles, Thakur et al. [59] investigated the scaling of cohesive forces and found that in terms of time evolution of the axial stress and porosity they should by scaled with a factor of l 2 instead of l 3 . The fluidization of cohesive particles was investigated by Tausendschön et al. [60]. The authors found that regarding the prediction of the average slip velocity, a stress-based scaling is superior to a Bond number-based scaling. The scaling of contact forces for wet and dry particles in a vertical mixer was examined Chan and Washino [61]. They compared different scaling rules (l 0 , l 2 , and l 3 ) and found for dry particles that the radial profile of the average particle velocity is invariant of the contact scaling method used. For wet particles, the authors show that l 2 -scaling of liquid bridge capillary and viscous forces lead to a good prediction of the original system. The question is whether the latter findings, regarding the invariance of contact force scaling for dry particles, can also be transferred to fluidized particle systems.
Therefore, the aim of this work is to compare different scaling rules proposed in the literature, and to compare their impact on pressure drop and particle dynamics in the field of fluidized beds. Industrial applications are characterized by a huge number of particles (> 10 11 ). To efficiently apply CFD-CGDEM for this kind of applications, scaling factors of l 3 10 5 are necessary to keep the number of parcels reasonable. To the authors' knowledge, the maximum scaling factors applied so far by other authors are limited to values below l 3 < 10 3 . Therefore, the limits of CFD-CGDEM are explored by conducting simulations with increasing values for l 3 . In the first part of this study, moderate scaling factors of l 3 125 are used, while in the latter part a factor of l 3 = 300 000 is applied. The results are compared with CFD-DEM results and experimental data.
Since the use of CFD-CGDEM inherently leads to a reduced spatially resolution, it cannot be ensured that all important local flow structures are resolved. This can cause an overprediction of the calculated drag forces, and with this, an overestimation of bed height. In analogy to filtered-KTGF, the loss of spatial information might be compensated by using filtered sub-grid drag models as the energy minimization multiscale method (EMMS). The impact of EMMS-based drag correction on the calculated static pressure and its standard deviation is investigated within the scope of this work.

Numerical Methods
The basis of the CFD-CGDEM model applied in this work is a soft-sphere CFD-DEM approach based on the model developed by [32]. A full description of the model used is given in the Supporting Information, Sect. S1.

Coarse-Grained DEM
The scaling rules presented in the literature can be categorized into two approaches. The first approach is based on the Buckingham Π-theorem and assumes that minimum fluidization velocity u mf , Archimedes number Ar = (Δrgd p 3 r g )/m g 2 , and the particle-based Reynolds number Re = (u r d p r g )/m g are kept constant during coarse graining (CG) [62,63]. With l = d p,CG /d p,0 being the ratio of the diameter of CG-parcel to that of the original particle, fluid dynamic viscosity and particle density are scaled as follows: where Eq. (1) is the result of keeping minimum fluidization velocity and Reynolds number constant, and Eq. (2) is derived from the invariance of Archimedes number. A direct scaling of acting forces is not done in this case, rather the forces are scaled only indirectly via changed material properties. This method has the disadvantage that changing material properties can lead to difficulties or very complex scaling rules if other transport phenomena, such as heat and mass transfer and chemical reactions, are included, as these are based on the original material properties. Most researchers use a different approach that is based on energy conservation. This approach assumes that all particles within a CG-parcel share the same translational velocity that matches the velocity of the original particle and rotate with equal mean velocity. Based on these assumptions, scaling rules for all acting forces can be postulated.
For the scaling of contact forces, different rules are proposed in the literature. Most researchers argue that the total energy has to be equal between CG-parcel and the original particle system, and, therefore, propose a scaling of F c,CG = l 3 F c,0 [49, 51-53, 57, 64-66]. In contrast, other authors use a scaling of F c,CG = l 2 F c,0 and explain it with either an invariant dimensionless normal overlap during the collision [67], the invariance of stresses [68], or the invariance of total energy, and the assumption that the contact time scales linearly with l [54]. This scaling approach leads to an invariance of the coefficient of restitution as shown by [60].
Other authors consider the voidage within the CG-parcel e CGP and argue with the invariance of kinetic energy and energy dissipation during the collision process [69,70]. Instead of an explicit scaling of the calculated contact forces, they scale the restitution coefficient by . A combination of scaled restitution coefficient and l 3 -scaled contact force can also be found in the literature [48].
To ensure the same rotational energy between both systems, the rotational speed of the CG-parcel should be scaled in the form w CGP = l -1 w 0 [64]. However, this is most often neglected [71]. Since pressure gradient and gravitational force are directly proportional to the particle mass, a consensus can be found in the literature that these forces should be scaled as F Δp,CG = l 3 F Δp,0 and F grav,CG = l 3 F grav,0 .
There are several inconsistent hypotheses of the scaling rules for the drag force in literature. Most researchers propose a scaling rule of F d,CG = l 3 F d,0 . They derive that rule from the necessity that the drag force acting on the CG-parcel should be equal to the sum of drag forces that act on all original particles within the parcel [51-53, 57, 58, 64-66, 72]. A scaling rule of F d,CG = lF d,0 is proposed by [71] and justified by the argument that the surface where the drag force acts on should be invariant. For the Stokes flow regime (Re < 1), a drag scaling of F d,CG = l 2 F d,0 is proposed by [68].
In the scope of this work, only energy conservation-based scaling rules are used. All simulations are done with the commercial simulation platform Simcenter STAR-CCM+ by Siemens Digital Industries Software. If not otherwise stated, the original implementation of CFD-CGDEM in Simcenter STAR-CCM+ is employed and, therefore, the following scaling rules are applied:

EMMS Drag Correction
Fluid-particle interaction forces are calculated based on the phase fraction of the dispersed phase within each cell. Therefore, for CFD-DEM as well as for CFD-CGDEM, the cell size should be bigger than the particle or parcel size to avoid unphysical phase fractions that lead to singularities in the equation system. The use of CFD-CGDEM inherently includes a coarsening of the mesh. This can result in unresolved flow features, which leads to an inaccurate drag prediction [73]. In analogy to MP-PIC and filtered-KTGF, the use of sub-grid drag models may improve the accuracy of CFD-CGDEM results in those cases.
In this work, the application of EMMS-based drag correction on simulation results is investigated. The EMMS model was developed by Li et al. [74]. The concept is based on the idea that the formation of local particle clusters is the reason for heterogeneity inside fluidized beds. It is assumed that particles are evenly distributed inside and outside a cluster. Several constraints and correlation-based closure conditions lead to an equation system that is still not fully closed. Therefore, a cost function is formulated that describes the energy consumed by the suspension and the transport of particles. The minimization of this cost function is the final closure needed. A full description of the model, including its equations, can be found in [75].
The output of the model is a heterogeneity index He d = b E /b that accounts for unresolved structures and describes the deviation from the applied drag model, which in most cases is a combination of the Ergun equation [76] for high packing fractions and the model of Wen and Yu [77] in the dilute regime. This model is often called the Gidaspow drag model [17]. Since the perpetual computing of the EMMS-based drag coefficient b E during a coarse grid simulation is numerically highly demanding, it is common practice to derive correlations from the solution of the optimization problem, which describe the heterogeneity index as a function of Reynolds number and voidage.
In the scope of this work, the correlation derived by Herbert and Reh [78] is used whenever an EMMS-based drag correction is applied. Its full formulation and parameters are given in the Supporting Information, Sect. S2.

Impact of Contact Scaling
When it comes to scaling of contact forces, two levels of scaling can be distinguished. The first one is the choice of length scale that is used for temporal discretization of the collision process, while the second one is the scaling of the calculated contact forces itself. In DEM, the recommended time step size is, among other parameters, limited by certain criteria, that are proportional to the particle diameter as a length scale. With increasing length scale, the time step size can be raised, which is beneficial for the performance of the simulation. Therefore, the question arises which length scale should be used for CFD-CGDEM simulations: the diameter of the original particle or that of the CG-parcel diameter? If the latter is possible without loss of accuracy, this would further increase the performance boost due to coarse graining, since not only a reduced number of particles need to be tracked, but also the time step size for the simulation could be enlarged.
As discussed earlier, different scaling rules for contact forces are proposed in the literature. It will be demonstrated how differences in contact scaling affect parameters like mean and RMS value of particle position and the axial void fraction profile along the reactor.
To evaluate the accuracy of CFD-CGDEM, a reference CFD-DEM simulation of a rectangular fluidized bed was conducted. The fluid phase was treated as incompressible gas. At the gas

Research Article
inlet, a constant mass flow rate was applied, while ambient pressure was specified at the outlet. All walls have no-slip boundary conditions. The most important conditions and properties of the CFD-DEM simulation are given in Fig. 2. According to [79], values of C fs = 0.15, C fr = 0.001, and e 0 = 0.97 are used for friction and restitution coefficients. Density and fluid viscosity were set to r g = 1.184 kg m -3 and m g = 0.01855 mPa s, respectively. A constant time step size of 50 ms was selected, and turbulent effects were considered using the realizable k-e turbulence model [80]. For temporal discretization of CFD-CGDEM simulations, two different approaches were employed. For the particle-based coarse graining, the time step is determined by the original particle diameter, while for the parcel-based coarse graining (PBCG), the parcel diameter is used. Additionally, two different contact scaling laws were applied. The contact forces are either scaled with l 3 (L3) or l 2 (L2). The remaining forces are scaled according to Eqs. (4)- (6). Because of CG, the basic DEM requirement that cells need to be bigger than parcels can get violated. This would result in unphysical local void fractions and cause divergence of the simulation. To circumvent this problem, momentum source smoothing using cell clusters (in the following abbreviated as source smoothing) is employed to evaluate the local void fraction. The cluster size for all simulations was set to 12 mm. A detailed description of the source smoothing approach is given later in the manuscript.
The average particle position h p is calculated by evaluating the solid phase fraction e s in each cell as: Where h k is the axial position of the cell centroid and V k denotes the volume of the respective cell. The mean average particle position is calculated as: (8) and its RMS value is given by: All simulations were run for 20 s physical time, while the time frame of 5-20 s was used to extract the simulation results.
Simulations with inlet velocities of u in = 1.2 m s -1 and u in = 1.6 m s -1 were conducted using the scaling laws and temporal discretization described above. CG and PBCG with l 3 contact scaling were conducted for scaling factors of l = 2-5. PBCG simulations using l 2 contact scaling were done for l = 2-3. The evaluated mean average particle position and its RMS values for both inlet velocities are given in Fig. 3. The results indicate that neither the scaling law nor the temporal discretization seems to have a significant impact on the mean average particle position or its RMS value. For scaling factors l 3, a very good agreement with CFD-DEM results (l = 1) can be found in terms of particle position. Nevertheless, for higher values, the mean average position of the particles gets slightly overpredicted, but still show a good to reasonable agreement.
The RMS of particle position characterizes the dynamic behavior of the particles. For most of the cases with increasing scaling factors, the RMS value gets more and more underpredicted in comparison to CFD-DEM results. Above findings are also presented in an animation provided as Video S1 in the Supporting Information. For PBCG temporal discretization and L3 contact scaling the particle movement for different scaling factors is illustrated. While the average bed height is predicted quite accurately for all configurations, the particle dynamics are predicted wrong for higher scaling factors. Even a slug-flow like flow regime is predicted for l = 5, which is not in accordance with the CFD-DEM results.
The impact of different temporal discretization and scaling factors on the time-averaged axial voidage profile is displayed in Fig. 4. It is obvious that almost no difference can be seen for PBCG and CG time discretization. A good agreement with CFD-DEM results can be observed for low scaling factors of l 3, while the deviations significantly increase for higher values. The results show that a parcel-based temporal discretization can be used without any loss of accuracy. This leads to a further speedup of the simulations since the DEM time step size can be significantly increased. The impact of the contact scaling rule on the axial voidage profile is indicated in Fig. 5. For the lower inlet velocity, almost no difference can be seen between L3-and L2-scaling. Only slight deviations can be observed if the inlet velocity is increased. The results demonstrate that contact scaling might not be as relevant as expected in that particular case and that a correct scaling of drag force might be of major importance. However, this might be related to the special configuration of the reactor, where,  Research Article due to the low reactor-to-parcel diameter ratio, wall effects may become dominant. The results of this study indicate that in terms of averaged values as mean particle position and voidage profile, the CFD-CGDEM approach leads to accurate results for low scaling factors. With increasing scaling factors, deviations get higher. Regarding the particle dynamics, the approach leads to an underestimation of RMS values of particle position indicating a loss of particle dynamics in the system due to CG. Furthermore, the results reveal that temporal discretization and contact scaling rule is not as important as expected.
This might be a surprising finding. However, Li and Kuipers [81] investigated the impact of particle restitution coefficient on flow structure. The authors found an impact of particle collisions on pattern formation, energy budget, and pressure fluctuations if particles are almost elastic (0.85 e 0 1). Below this limit, the effect of the coefficient of restitution vanished. Bed height and time-averaged pressure drop proved not to be affected by choice of the coefficient of restitution at all, which is in agreement with our findings. Nevertheless, the discussion demonstrates that further investigations are necessary to understand the impact of contact scaling laws on accuracy on CFD-CGDEM.
One possible explanation for increasing deviations at higher scaling factors might be related to the reactor-to-parcel diameter ratio, which decreases if l is raised. For the investigated design, a very low diameter ratio of 20 results for a scaling factor of l = 5. Recent investigations by Lu and Benyahia [82] indicate that larger scaling factors can be taken for large-scale reactors. Industrial applications are characterized by a huge number of particles, on the order of magnitude of 10 11 -10 13 . Therefore, to make CFD-CGDEM appropriate for this kind of applications, scaling factors of 10 5 are needed. This corresponds to a parcel-to-particle diameter ratio of l 46. To the authors' knowledge, CFD-CGDEM with scaling factors that high have not been conducted so far. In the following section, CFD-CGDEM will be applied to a semi-circular lab-scale unit that allows the use of very large scaling factors without running into a possible limitation concerning the reactor-to-parcel diameter ratio.

Impact of Drag Calculation
In the case of fluidized particle systems, it is widely accepted that an accurate prediction of drag force is of major importance to describe particle and fluid dynamics right. As already discussed, coarsening of particles and cells can lead to unresolved flow features, and, therefore, a misprediction of drag. Based on simulations done for a semi-circular unit with a diameter of 304.8 mm, the impact of large scaling factors, drag calculation, Chem. Eng. Technol. 2020, 43, No. 10, 1959-1970 2020 The Authors. Published by Wiley-VCH GmbH www.cet-journal.com  and cell coarsening on the results in terms of pressure drop and pressure fluctuations is investigated. The results are compared with experimental data. Experiments were conducted in the lower fluidized-bed region of a semi-circular circulating fluidized-bed reactor at Particulate Solid Research, Inc. The unit has a width of 0.2858 m and a maximum depth of 0.1517 m. Compressed air with a superficial velocity of u in = 0.107 m s-1 is supplied from the bottom to the plenum and distributed with a perforated plate. Two pressure taps were used in this study, one at a position slightly above the gas distributor, the other at a height of 0.9144 m (36''). A detailed description of the experimental setup can be found in [83].
The reactor was filled with fluid catalytic cracking (FCC) catalyst particles (Geldart Group A) that have a particle density of r p = 1200 kg m -3 and a Sauter mean diameter of d p,32 = 64 μm until an initial bed height of approximately 1.66 m was reached. The underlying particle size distribution is displayed in Fig. 6.
For the numerical simulations, the computer-aided design (CAD) description of the reactor was meshed using hexahedral dominated trimmed cells with an almost uniform edge length of Δx 11 mm. For the calculation of the flow field, the gas was treated as constant density fluid with a density of r g = 1.184 kg m -3 , and a dynamic viscosity of m g = 0.01855 mPa s. Turbulence effects were considered by using the realizable k-e turbulence model. A no-slip wall boundary condition was used for the reactor wall. At the inlet the mass flow is specified, based on the respective inlet velocity. Ambient pressure conditions are specified at the outlet.
For the particulate phase, the PBCG-L3 approach with a scaling factor of l 3 = 300 000 was used. This corresponds to a parcel-to-particle diameter ratio of 67. The total number of parcels in the domain was 682 447. The Hertz-Mindlin contact model was applied to calculate parcel-parcel and parcel-wall interactions, using a static friction coefficient of C fs = 0.6, a rolling resistance of C fr = 0.001, and a restitution coefficient of e 0 = 0.5. The time step size was set to a value of 50 μs for all simulations.
Inherent with coarsening the particles is a coarsening of the volume cells, since the fundamental requirement of CFD-DEM concerning the cell size in relation to particle size still needs to be fulfilled. However, for complex geometries like reactors which include internals, a coarse geometrical discretization is not always possible. One approach to circumvent this problem is the aggregation of adjacent cells into cell clusters that are bigger than the original cells. The cluster is used to calculate the voidage based on the number of particles within the cluster. This value is later mapped to all cells that build the cluster to calculate the drags forces based on the averaged voidage and the still cell-based particle slip velocity. To investigate the effect of source smoothing on the results, all simulations were performed with or without that method. If source smoothing is used, eight volume cells are aggregated into one cluster.
As in the previously investigated cases, the Gidaspow drag model [17] was employed. However, this time, also the impact of EMMS-based drag correction is investigated by performing the simulations with or without applied drag correction.
The progression of pressure drop inside the bed evaluated between a height of 0 and 0.9144 m is depicted in Fig. 7. The results show that pressure drop is significantly underpredicted if the uncorrected Gidaspow drag law is used and source smoothing is applied. A relative deviation of approximately 50 % can be observed. By excluding the source smoothing method, the calculated pressure drop increases significantly by around 1500 Pa and only slightly underpredicts the experimentally determined value, showing a relative deviation of less than 14 %. It can also be seen that pressure fluctuations get damped because of source smoothing. A possible reason for source smoothing significantly impacting results could be related to  the rather large ratio of the cluster-to-parcel size of 2.8 which was used.
If EMMS-based drag correction is incorporated, the calculated pressure drop significantly increases, which is expected since the corrected drag coefficient is lowered, which leads to a reduction of bed height. When source smoothing is used, the calculated pressure drop slightly overpredicts the experimental value, showing a relative deviation of approximately 16 %. Without source smoothing, the predicted value increases slightly, leading to a relative deviation of 29 %. The simulation results indicate that EMMS-based drag correction minimizes the impact of source smoothing since it accounts for not resolved flow structures.
Differences between EMMS-based simulation and experimental results might be related to the heterogeneity index correlation that was used, which was determined for particles (r p = 2500 kg m -3 and d p = 300 mm) that were only roughly similar to the ones used in experiments and simulations. Better suited correlations may improve accuracy significantly. Videos S2 and S3 show the simulated particle dynamics for the investigated cases and video material taken at the experimental facility. Bed height is significantly higher if the uncorrected Gidaspow drag law is applied. Qualitatively the EMMS-corrected results are much closer to what was visually observed in experiments. The animation also illus-trates the reduced impact of source smoothing if EMMS-based drag correction is employed.
The simulation results regarding time-averaged pressure drop and standard deviation of its fluctuations are summarized in Fig. 8. The static pressure at h = 0 , which corresponds to the overall pressure drop, is predicted accurately for all configurations. A maximum relative deviation of 2.5 % is found if the uncorrected Gidaspow drag law with source smoothing is used. However, the good prediction of total pressure drop is not surprising since it is an integral value that only depends on total particle mass within the system. Nevertheless, this indicates that the total pressure drop is an insufficient quantity for validation. Regarding the standard deviation of pressure fluctuations, using source smoothing leads to a reduction of pressure fluctuations for uncorrected and corrected drag laws. For h = 36 ,  best accordance is achieved if no source smoothing is employed. In this case, a relative deviation of 5 and 23 % is observed for EMMS-corrected and uncorrected simulations, respectively. At h = 0'', pressure fluctuations are mispredicted by almost all simulations. This could possibly be related to unresolved effects from the distributor plate since the plenum and gas distributor were not explicitly modeled in the simulations. However, to some extent, differences regarding bed fluctuations will be unavoidably related to CFD-CGDEM since the loss of spatial resolution will inherently lead to a loss of heterogeneity that causes fluctuations.

Conclusion and Outlook
In the scope of this work, the impact of different parameters on the accuracy of CFD-CGDEM was investigated. A comparison between CFD-DEM and CFD-CGDEM of Geldart Group D particles in a rectangular fluidized bed shows that scaling of contact forces seems to be of minor importance in the case of fluidized beds. It was also demonstrated that a parcel-based temporal discretization could be chosen without affecting the results. This leads to an additional speed-up of the simulations. For scaling factors l 3 > 3, increasing deviations were observed. Simulations based on an extended semi-circular unit with a diameter of 304.8 mm, however, have shown that CFD-CGDEM with huge scaling factors, in this case l 3 = 300 000 are possible, which makes the model applicable for industrial-scale simulations. This indicates that the maximum allowable scaling factor might be a function of the reactor-to-parcel dimension ratio. More studies in that direction are necessary to identify a general criterion for this. Furthermore, the results of the semicircular unit show that drag correction might be necessary if very high scaling factors are used. Otherwise, unresolved flow structures may lead to a significant misprediction of drag force, which results in an inaccurate local pressure and bed height.
Although the pressure drop was also predicted right using the uncorrected Gidaspow drag law, EMMS-based drag correction was identified as a suitable model since it additionally was able to capture pressure fluctuations within the particle bed. Furthermore, the EMMS-based results were qualitatively a bit closer to which was visually observed in the experiments. Yet, if this is also true for other reactor configurations or if other drag correction models are better suited should be further investigated. The EMMS model only accounts for the loss of information due to grid coarsening. A model that additionally considers the loss of information with respect to the particle velocity distribution, caused by the particle coarsening as it was recently presented by [84], might further enhance the accuracy of CFD-CGDEM.
It was demonstrated that source smoothing based on cell clustering could significantly impact the simulation results; however, in the present case this effect could be minimized by using EMMS-based drag correction. Nevertheless, this important aspect must be further investigated, since a working source smoothing scheme is essential for the simulation of industrial applications that include internals and complex geometries. Besides the source smoothing based on cell clustering, also other source smoothing methods, e.g., diffusive smoothing [85], simple and Laplacian LES filter, or mesh-adaptive coarse graining [86], need to be evaluated and compared against each other.
Overall, the results of this study illustrate that CFD-CGDEM is a promising method for the numerical description of industrial-scale fluidized particle systems. Nevertheless, more effort should be put in developing general guidelines concerning scaling rules, maximum scaling factors, and the impact of subgrid drag correction for a series of different particle systems (fluidized and non-fluidized) to (a) identify the level of accuracy that can be reached using CFD-CGDEM and (b) assure that this accuracy can be achieved by following a certain best practice. contact forces are scaled by l 3 MP-PIC multiphase particle in cell PBCG parcel-based coarse graining PDF probability density function