Comparison of sub-grid drag laws for modeling ﬂuidized beds with the coarse grain DEM–CFD approach

Fluidized particulate systems can be well described by coupling the discrete element method (DEM) with computational ﬂuid dynamics (CFD). However, the simulations are computationally very demanding. The computational demand is drastically reduced by applying the coarse grain (CG) approach, where several particles are summarized into larger grains. Scaling rules are applied to the dominant forces to obtain precise solutions. However, with growing grain size, an adequate representation of the interaction forces and, thus, representation of sub-grid effects such as bubble and cluster formation in the ﬂuidized particulate system becomes challenging. As a result, particle drag can be overestimated, leading to an increase in average particle height. In this work, limitations of the system-to-grain ratio are identiﬁed but also a dependency on system width. To address this issue, sub-grid drag models are often applied to increase the accuracy of simulations. Nonetheless, the sub-grid models tend to have an ad hoc ﬁtting, and thorough testing of the system conﬁgurations is often missing. Here, ﬁve different sub-grid drag models are compared and tested on ﬂuidized bed systems with different Geldart group particles, ﬂuidization velocity, and system-to-grain diameter ratios.


Introduction
Fluidized beds play an essential role in the chemical and process industry [1].An upward-directed gas stream fluidizes a bulk of particles in a fluidized bed.Their good heat and mass transfer characteristics make it an efficient system for multiple chemical or physical applications [2,3].However, such systems are inherently unstable, and effects occur on multiple time and length scales.On a microscopic level, particle-particle and particle-wall interactions, as well as the interaction of particles with the fluid, are dominant.On the mesoscale level, bubble-and cluster formation dominate due to local fluid flow effects.These strongly affect the macroscale, e.g., the fluidization regime and mixing characteristics [4].
Many approaches have been proposed to model particlefluid multiphase systems to account for the different time and length scale phenomena [5].Direct numerical simulations (DNS) can be applied on a microscopic scale.Here, fluid dynamics is solved by computational fluid dynamics methods such as Lattice Boltzmann (LBM) or finite volume method (FVM).The particle surface can be resolved using the immersed-boundary method [6][7][8].Particle movement is solved with the discrete element method.The methods are computationally very demanding and limited to small-scale systems [9].Mesoscale effects can be resolved using either the two-fluidmodel (TFM) or coupling DEM and particle-unresolved CFD (DEM-CFD).TFM is an Euler-Euler approach in which particle-particle, particle-wall, and particle-fluid interactions are modeled by closure conditions such as the kinetic theory of granular flow [10].Much research is focused on developing accurate closure models, as shown by [11][12][13].DEM-CFD is a Lagrange-Euler method where particle movement is tracked by solving Newton's second law.Particle-particle and particle-wall interactions are described using contact models such as the Hertzian model or the linearspring-dashpot model [14].Closures are only necessary to model the fluid-particle interaction with appropriate drag laws.However, the well-resolved TFM needs a grid resolution with a cell size less than 10 particle diameters [15], which is not feasible for industrial applications.Also, the DEM-CFD approach is computationally very demanding and not feasible for industrial systems with particle numbers ranging from 10 8 to 10 12 [5].
A common way to decrease the computational demand is to reduce the cell number by using a coarse mesh.This can be done for both TFM and DEM-CFD.In the TFM, both the description of the solid and fluid phases is affected.In the DEM-CFD only the fluid phase is implicated by a mesh coarsening.In the DEM, the computational demand with regard to the solid phase can be decreased by reducing the number of particles following the coarse grain (CG) approach [16].Here, several original particles are lumped into a larger grain under the assumption that the individual movement of a particle is not crucial to the fluidized bed mesoscale characteristics as shown by [17].The CG approach gives the possibility to simulate large-scale reactors in a realistic time frame.Some large-scale examples are given in [18][19][20].However, if the cells as part of a TFM or CG-CFD-DEM are too coarse in combination with a too strong coarse graining for the latter, it can lead to an inadequate representation of the mesoscale characteristics.A particle moving in a particle cluster experiences less drag than individual particles [1].Neglecting this phenomenon in coarse simulations, thus assuming the drag in a grain is similar to the sum of the particles in the grain, where each particle exerts drag as a single particle, can lead to an overestimation of drag force and, thus, an erroneous estimation of bed expansion [21][22][23].To address the aforementioned issue, much research has been conducted on so-called sub-grid drag models that account for the mesoscale effects while still using a coarse mesh or coarse mesh and coarse grains.Sub-grid drag models can be classified into empirical, structure-based, energy minimization-based, filtered, and artificial neural networkbased models [24].The sub-grid drag model can be related to empirical correlations, for example, the O'Brien and Syamlal drag law, which is based on an air-FCC system with solid circulation fluxes [25].Wang et al. [26] developed a model based on the assumption that for Geldart B and D particles, interphase drag, e.g., drag between the solid and fluid phase, only occurs in a particle-dense phase.In contrast, the dilute phase, e.g., bubbles, remains free of particles and free of interphase drag.This leads to a two-phase structure of the system for the drag estimation, which belongs to the category of a structure-based drag law.
A large and well-investigated drag law group is based on the concept of energy minimization of particle transport and suspension [27,28].The energy minimization multi-scale (EMMS) model has been expanded to account for bubbly and circulating fluidized beds recently [29,30].Lastly, filtered drag correlations are often applied.Here, a correction factor is introduced, obtained by regressing results of finely resolved TFM or DEM-CFD simulations [31][32][33].
The filtered drag models can tend to ad hoc regressed results, so applicability to other system configurations always needs evaluation.A recent approach in the development of filtered drag laws involves leveraging artificial neural networks for regression purposes [34,35].Lu et al. [35] performed a comparative analysis between filtered drag laws derived with artificial neural network (ANN) and the nonlinear regression technique.In this case, the filtered drag law generated through nonlinear regression showed more accurate results.Moreover, attempts to couple the ANN method to the EMMS model were successfully made to extend EMMS to multiple macroscale markers such as operating conditions, e.g., gas inlet velocity, system width, or even material properties that are averaged over the cross-sectional area of the fluidized bed system [36,37].
The heterogeneous drag laws are, in most cases, developed for the application of coarse TFM simulations.Few researches have been conducted on the application of heterogeneous drag laws in the coarse grain DEM approach framework.Lu et al. [38] demonstrated the general applicability of the EMMS model with coarse grains and showed good agreement with experimental data.Moreover, Jurtz et al. [39] investigated the EMMS approach for a bubbling fluidized bed with Geldart class A particles, which improved the results compared to a classical homogeneous model.The EMMS model has also been applied to DEM simulations [40] and for an industrial example of a polymerization process using coarse grains [41], which was also in good agreement with experimental data.Radl et al. [31] extended their filtered drag law to the CG approach by introducing a factor reducing the drag on the grains with increasing grain size.However, they stated that their model might be ad hoc regressed, and it is only tested for coarse grain factors up to 8.
As can be seen, research about the applicability of such models to the CG-DEM-CFD approach has been demonstrated for EMMS and filtered drag models.However, a detailed analysis of the applicability of sub-grid drag models to CG-DEM-CFD simulations is still missing.This work compares popular drag models using the CG-DEM-CFD approach for fluidized beds.Often the drag models are limited to specific operating conditions of the system.Therefore, limitations of the drag laws concerning fluidization regime, system width, and particle types (Geldart class B and D) are analyzed while changing the coarse grain factor, thus, the number of grains in the system.Geldart class A particles play an essential role in fluidized bed engineering and will be the focus of future work.

Numerical methodology
The DEM is used to solve the movement of particles.Within this method, the soft sphere approach is employed with the linear spring-dashpot model to model the interaction between particles.For fluid movement, computational fluid dynamics (CFD) is used based on fluid phase continuity and momentum conservation.The DEM and CFD are coupled and solved as part of the software MFiX [42,43].As MFiX does not provide for the sub-grid drag laws considered here, MFiX is extended using the usr_drag module provided by MFiX.
Regarding the DEM-CFD, a full description of the model and applied equations is given in Sect.A.1.

Coarse grain DEM
Coarse graining is performed by summarizing several original particles into larger grains.The coarse grain number n cg = n p /n g is defined as the number of particles lumped into a grain.The coarse grain factor, also known as grain-toparticle size ratio, is defined as: where d g is the grain diameter and d p is the particle diameter.Contact forces in CG-DEM are calculated by assuming compact grains, so the total volume and the total mass are constant.A rescaling of the interaction forces is needed to represent the original particles properly.Several scaling approaches have been developed recently.For our simulation system, the scaling method by [44] combined with [45] was proven one of the most accurate looking at multiple output parameters in our previous work [46] and it is applied here.The stiffness components in normal and tangential directions are scaled with the coarse grain number: whereas, for the damping components, the ratio of normal to tangential damping is kept constant.The restitution coefficient is modified to represent intra-grain collisions such as: Lastly, when sliding occurs, dissipation is conservated by scaling the friction coefficient μ as follows As hydrodynamic forces, pressure gradient force F ∇ p p and drag force F d p are considered.The pressure gradient force acting on the grain is assumed to be similar to the sum of particles inside the grain: with V g as the volume of the grain.In our previous work [46], it was found that the scaling of drag force in a similar way gives the most accurate results, such as: where ε f is the void fraction, u f and v g are the fluid and grain velocity, respectively, and β is the interphase momentum transfer coefficient.Drag models depending on slip velocities and void fractions are applied to calculate this coefficient.In this study, the drag law developed by Gidaspow [47] is used, which belongs to the class of homogeneous drag models.It is a combination of the Ergun drag law [48] for dense regions and Wen-Yu drag law [49] for dilute regions: where μ f is the fluid viscosity, ρ f is the fluid density and v p is the particle velocity.The drag coefficient C d is calculated as: where Re p is the particle Reynolds number, which is defined as: However, for coarse simulations, either by applying a larger coarse grain factor or by applying a larger grid size, the accurate representation of sub-grid behavior becomes challenging.Therefore, modification of the interphase momentum transfer coefficient is performed by introducing a heterogeneity index that accounts for the mesoscale effects: Here, β Wen-Yu is the interphase momentum transfer coefficient obtained by the Wen-Yu drag law.β Het is the interphase momentum transfer coefficient obtained by heterogeneous sub-grid drag laws.The sub-grid drag laws will be explained in the following section.

Sub-grid drag laws
This subsection elaborates on the sub-grid drag laws used in this study.All necessary equations used to implement the drag laws in the software MFiX are given in Sect.A.2.The structure-based drag law is explained first, followed by the EMMS drag law group and the two filtered drag laws.
Structure-based drag law: As a structure-based drag law, the drag law by Wang et al. [26] is applied.The modification by Wang et al. to the homogeneous drag model is based on the assumption that the drag between solid and fluid only occurs in a particle-dense phase.For Geldart class B and D beds, bubbles remain primarily free of particles, which justifies the assumption for the drag for these cases.In the model, the drag force in each CFD cell is adjusted to only account for the fraction of the cell occupied by the particle-dense phase.The fraction of the particle-dense phase is calculated based on a Reynolds-Archimedes correlation.The interested reader is referred to Wang et al. [26] for detailed information.The model was initially designed for the TFM.The modified model improved results compared to the unmodified drag law looking at a fluidized bed system at the pilot and at the industrial scale [26].For the CG-DEM-CFD approach, the model has not been tested yet.Nevertheless, the underlying assumptions are also valid for a Lagrange-Euler framework, which is why this work investigates the model.

EMMS group:
Two EMMS drag laws, the EMMS/Matrix [27,28], and the EMMS/steady-state model [29,30], are investigated in this work.The EMMS/Matrix drag law, initially developed for flow in risers, accounts for the heterogeneous behavior by establishing transport equations of the particle flow in the dense, particle-rich phase, the dilute, fluid-rich phase, and the interaction of dense and dilute phase of the fluidized bed.This set of equations leads to an unclosed system of equations.Therefore, a stability criterion is introduced to minimize the energy consumption for suspending and transporting the particles in the flow [27,28].The solution to the minimization problem is an interphase momentum transfer coefficient for a specific range of slip velocities and void fractions.From this, the heterogeneity index can be calculated with Eq.The filtered drag law encompasses both fluid coarsening, which refers to the coarsening of the CFD cell size, and coarse graining.A correction factor is introduced for coarse grain DEM to decrease the interphase momentum transfer coefficient for increasing grain size.However, their model is only tested for l < 8, so inaccuracies may occur for larger grain diameters and different .The original drag law was developed for TFM simulations.Now, the revision allows using the drag law for DEM-CFD simulations, which makes it more applicable for CG-DEM-CFD.Moreover, by introducing a macroscale marker for the gas inlet velocity, the model is applicable to systems with a wider range of gas inlet velocities.It was shown that coarse grain simulations performed with filtered drag laws derived from DEM-CFD simulations yield more accurate results than filtered drag models derived from TFM simulations [35,51].Similarly to the Radl drag law, the heterogeneity index is calculated with regression parameters that differ with void fraction and slip velocity.More information is available in [35].
A comparison of the heterogeneity index as a function of the void fraction for the different sub-grid drag laws is given in Fig. 1 for a constant slip velocity of 0.98 m/s for Geldart B particles and a gas inlet velocity of v in = 1.7 m/s.For the regression based models, the ratio of the filter size to characteristic cluster length scale f /L c equals 1.08.The Wang drag law [26] decreases drag compared to the Gidaspow drag law until a void fraction of ε f = 0.55, leading to a heterogeneity index below 1.From there, the heterogeneity index reaches values up to 4. The EMMS/Matrix model, as originally established for fast risers, lowers the drag for the void fraction area 0.42 < ε f < 0.99 between 1 and 2 orders of magnitude.No correction is performed for void fractions below 0.42 and above 0.99.The EMMS/steady-state model reduces the drag for 0.45 < ε f < 0.85 up to one order of magnitude, but a drag correction is not conducted for dense and dilute areas.The drag law by Radl et al. shows a reduction of He = 1 to He = 0.7 up to a void fraction ε f = 0.6.Subsequently, the heterogeneity index slightly increases again but

Simulation design
The simulations are designed to assess the drag laws' application range thoroughly.For this, a cylindrical system configuration is chosen with three different system widths to assess whether a limitation on the coarse grain factor or the ratio of system diameter to particle/grain diameter N exists.Moreover, two gas inlet velocities are chosen to account for two fluidization states.Monodisperse particle beds are generated with Geldart B or Geldart D particles to assess the effectiveness of the sub-grid drag laws for different particle types.Simulations are performed with the open-source software MFiX [42,43].The properties for the unscaled benchmark particle systems are given in Table 1.
For the particle properties, the properties of the particleparticle interactions are given.The properties of particle-wall interactions are the same, except for the friction coefficient, which is reduced from 0.1 to 0.09.Rolling friction is not considered because appropriate scaling rules of rolling friction for the coarse grain approach are still a subject of current research [16].The benchmark simulations are then repeated with coarse grain factors l = 2, 4, 6, 8 and 10, applying the scaling rules explained in Sect.This leads to a constant system-to-particle ratio N p of 100 to 50 and 25, respectively.A constant initial bed height is ensured by reducing the particle number from 4 × 10 5 to 2.25 × 10 5 and 1 × 10 5 particles, respectively.The spring constant is adjusted for Geldart B and Geldart D type particles to keep the maximum overlap between particles or particles and wall below 1% as recommended by [52] but ensure a reasonable DEM time step.The ratio between particle collision time and DEM time step is set to 25.The CFD time step is set constant to 1 × 10 −4 s.The CFD cell size is set to dx = 2d p for unscaled simulations and dx = 2d g for coarse grain simulations in each direction.The cells are cut into smaller cells at the boundaries to resolve the cylindrical geometry of the domain using the cut-cell method.Therefore, the divided particle volume method (DPVM) is used for mapping Lagrangian information on the Eulerian mesh to avoid numerical instabilities.Traditionally, void fraction and other interphase quantities are assigned to the CFD cell, in which the centroid of the particle is, which is called the particle centroid method [53].For DPVM, these quantities are also distributed to the neighboring CFD cells based on a distributing weight function, which allows using smaller cell-to-particle ratios [43,53].To map the gas velocity back to the particle's position, two steps are performed in MFiX.At first, the gas velocities are mapped to an interpolation grid.
At second, a second-order accurate Lagrange polynomial is Fig. 2 Geometries for the three different system widths and boundary conditions used to interpolate the gas velocity from the interpolation grid to the particle's position [54].The gas flows through the domain using a velocity boundary condition at the bottom and a pressure outlet at the top, as illustrated in Fig. 2.
No-slip boundaries are set for the walls of the fluidized bed.
The physical time simulated is 20 s.The gas inlet velocity is ramped in the first second of the simulation to ensure smooth fluidization.
For evaluation of the simulation accuracy, the normalized average particle or grain height is calculated.Therefore, at first, the average particle or grain height h < p/g> for each time frame is calculated using the formula by [55]: where h k is the axial position of the center of cell k and V k is the volume of the cell k.At second, the average particle or grain height is then normalized by dividing the initial average particle height h < p,init> of the static bed given in Table 1.
The general simulation procedure is as follows: In the first step, unscaled DEM-CFD simulations with the Gidaspow drag law are performed for each system configuration, which is used as a reference.Three different system diameters are investigated to analyze the limitation of the system-toparticle/grain diameter ratio N .In the second step, each configuration is repeatedly simulated using coarse grain factors between l = 2-10 with the homogeneous Gidaspow drag law.The simulations with a system-to-particle ratio of N p = 50 are only repeated until l = 8 to ensure at least 3 CFD cells in each direction and maintain a cell size of dx = 2d p/g .The limitations of resolving sub-grid phenom-ena as a function of the system-to-grain diameter ratio or the coarse grain factor are analyzed and determined.In the third step, the simulations with large coarse grain factors of l = 6−10 and l = 6−8 for N p = 50 are repeated with the heterogeneous drag laws defined in Sect.2.2.Most subgrid drag models are optimized on a specific particle type and fluidization regime.Geldart group B and Geldart group D particles are analyzed to assess applicability for different particle types.The simulations are conducted in bubbling and slugging to turbulent regimes to cover two different types of fluidization.

Results and discussion
As described above, this section is divided into two subsections to analyze and discuss the impact of drag laws on the accuracy of the coarse grain approach.Firstly, the cases with the homogeneous Gidaspow drag law are analyzed for various coarse grain factors.Secondly, the cases showing significant deviation from the unscaled benchmark cases are then analyzed using the different sub-grid drag laws.The goal is an improvement of the results for all tested configurations.

Effect of coarse graining using the Gidaspow drag law
In Fig. 3, the normalized average particle or grain height for Geldart class B (Fig. 3a) and D (Fig. 3b) particles is shown by using the Gidaspow drag law under varying system-toparticle/grain diameter ratios N. The ratio N is varied by changing the system diameter and the particle/grain diameter.Here, the normalized average particle or grain height is temporally averaged over the time frame of t = 5−20 s, which results in one value for each simulation.
Looking at the Geldart class B systems, almost no impact of coarse graining with respect to the normalized particle height can be observed until a ratio of N < 12.5 for both inlet velocities, except for the system with a ratio of N p = 50, where the bed height is increasing only at a ratio lower than N < 9.4.A ratio of N = 12.5 corresponds to a maximum coarse grain factor of l = 8 in the wide system and l = 4 in the narrow system, as seen in the color bar, representing the coarse grain factor l. Below N < 12.5 and N < 9.4, a relatively sharp increase in particle height is observable.Interestingly also is the influence of absolute system width.For similar ratios N , the normalized particle height in the wide system only increases from 4.04 to 4.27 even for coarse grain factors of l = 10 (N = 10), whereas for the width of N p = 50 and N p = 75, the bed height rises sharply from 4.04 and 4.17  For Geldart D class systems, at the large inlet velocity of v in = 5.2u mf , the normalized particle height is constant for ratios N until N = 25, with one exception for N p = 75.There, at N = 37.5, the value is already overpredicted, meaning the deviation to the benchmark case is more than 5%.Similar, to the Geldart class B, there is also an absolute influence of the system width as a narrower system width is leading to larger normalized particle heights at the same ratio N .For the low gas inlet velocity of v in = 2.6u mf , the normalized particle height is first underestimated and then increased for larger grain sizes.Therefore, no limitation regarding the coarse grain factor or ratio N can be concluded for the low gas inlet velocity for Geldart class D particle systems.
As shown in Fig. 3, an overprediction of the particle height appears already at larger N for both particle types when using a larger gas inlet velocity.The Geldart class D particles are less elevated because the fluidization regimes differ depending on particle type.However, the deviation with decreasing ratio N is slightly larger for Geldart D particles.For Geldart class D, the lowest normalized particle heights are reached at around 2.6−2.8 for benchmark cases, but reach values up to 4 for low ratios N .For Geldart class B, the increase is from 4 for benchmark cases to 5 for low ratios N .The reason for the larger error using Geldart class D particles might be explainable by two effects.Firstly, at the large inlet velocity, different fluidization regimes are reached for Geldart B particles and Geldart D particles.A good indicator for the fluidization regime is the intensity of pressure fluctuations in the fluidized bed [1,56].Therefore, a fast Fourier transformation (FFT) on the pressure drop signal is performed.The pressure drop is calculated between the inlet and the outlet of the fluidized bed system.A sampling rate of 100 Hz considering data taken from the time frame of t = 5−20 s is chosen.In Fig. 4, the results of the FFT are shown in the form of the power spectral density (PSD) plotted against frequency for the unscaled benchmark cases at the two velocities.Exemplary, the results from the widest system (N p = 100) are shown.For the simulations with the Geldart class B particles, the power spectral density always stays below 1500 Pa 2 /Hz.The main peaks at the low gas inlet velocity are between 5 and 15 Hz, and the PSD goes up to 250 Pa 2 /Hz.For the gas inlet velocity at 14u mf , the dominant peak is at 5 Hz with a much higher PSD but narrower width, which is typical for the slugging regime [1,56].Looking at the simulations with Geldart D particles, the power spectral density is one magnitude higher.Interestingly, the simulation at lower gas inlet velocity shows much higher power spectral densities at specific frequencies than at high gas inlet velocity.A decrease in the PSD is typically seen when going from a bubbling or slugging regime to a turbulent regime [1,56].Therefore, the frequency analysis indicates that the system with Geldart D particles at the high gas velocity is reaching the turbulent regime.In contrast, the system with Geldart B particles is reaching the slugging regime.For the turbulent regime, a higher heterogeneity of the fluidized bed is expected.This might lead to overprediction in particle height as the coarse systems cannot resolve the heterogeneous characteristics of the fluidized bed.Secondly, the difference in absolute gas inlet velocity can lead to the different behavior of Geldart B to Geldart D grains.The fluid velocity to achieve the higher fluidization regime is u in = 2.75 m/s for Geldart class D particles, and u in = 1.7 m/s for Geldart class B. This could result in higher slip velocities being attained by Geldart D particles, potentially Moreover, it is also observable that the coarse grain results in the narrow system N p = 50 have a higher deviation from the benchmark than in the wider systems.Therefore, there might also be other effects causing the deviation in the narrow systems for both Geldart types.There can be four reasons why the effect cannot be corrected.Firstly, the D/d-ratio N can be too narrow to distribute the particles in the radial direction leading to larger distribution in the axial direction.Secondly, the portion of particle-wall interaction compared to particleparticle interaction increases significantly and might lead to a deterioration of the results.Thirdly, the large grains distort the typical fluid flow pattern in narrow configurations.Lastly, due to the inherent mesh coarsening, the adequate resolution of fluid flow and interaction quantities for drag, such as porosity and slip velocities, become challenging.For example, for l = 8 and N p = 50, there are only 9 grid points in the radial direction.In addition, the coarse grid might lead to erroneous interpolation of the velocities between the no-slip boundary condition at the cylinder wall to the fluid velocity at the position of the particle depending on the interpolation scheme applied.A second-order accurate Lagrangian polynomial is already used.Switching to a higher order, in this case a fourth-order accurate Lagrangian polynomial did not enhance the results and lead to an even higher overestimation of the normalized particle height.Interpolation techniques exist to reduce the CFD cell-to-particle ratio and might lead to improved results.However, interpolation errors can appear and must be assessed for such methods, so they were not applied in these cases.
Figure 5 shows snapshots of the wide simulation system (N p = 100) taken at 10 s to give a visual example of the change of the fluidized bed dynamics while applying coarse graining.Exemplary, the simulation results for particle velocity for Geldart B particles with an inlet velocity of u in = 14.2u mf are shown.For the benchmark case, l = 1, we can see most of the particles concentrating in the lower part of the system.Some particles move into the upper part, indicating a dense region and a dilute region.The separation between the dense and the dilute phases is still observable for a coarse grain factor of l = 4, which reduces the particle number from 400,000 to 6250.Although particles are not elevated as high as in the unscaled benchmark case, the bed average void fraction profile is shifted to a lower void fraction with increased bed height.This effect is even intensified for the coarse grain factor of l = 10, as also already shown in Fig. 3, the normalized particle height is increased.
The snapshot shows that particles are elevated, but differentiation between the dense and dilute regions is not possible anymore, as one grain holds 1000 particles.A similar effect is also observable by looking at the contour plots of the void fraction in Fig. 6.Here, exemplary shown for Geldart class D particles in the wide system (N p = 100) at a gas inlet velocity of u in = 5.2u mf .For l = 1 and l = 2, several bubbles and small particle clusters are visible.For l = 4 and l = 6, the bubbles turn into larger voids.A differentiation between dense and dilute regions becomes challenging for large coarse grain factors of l = 8 and l = 10.
To sum up the subsection, the simulations performed with the Gidaspow drag model show no overprediction of particle height for low coarse grain factors and thus, larger D/d-ratios Fig. 6 Representative contour plots for void fraction profiles taken at 5 s of the simulation for Geldart D type particles at fast fluidization u in = 5.2u mf and a system width of N p = 100.The black line represents a void fraction of ε = 0.9.The decreasing resolution of the simulation leads to erroneous bubble determination and a loss of the ability to determine particle clusters, which is observable looking at the slices going from l = 4 to l = 6 N .When using larger coarse grain factors, so low ratios N , the bed height is overpredicted, probably due to an overestimation of the drag force.This effect is intensified when using higher gas inlet velocities.At low gas inlet velocities, most particles concentrate at the bottom of the system and the influence of drag force on the particles is rather small, which leads to accurate results or even underprediction of particle height.However, for larger gas inlet velocities, the heterogeneous behavior of the fluidized bed can be intensified, which might be less resolved with the coarse simulations.As the deviation is dependent on the gas inlet velocity, the overestimation might increase with large gas velocities, looking, for example, at circulating fluidized beds or risers.When comparing the fluidized bed simulations containing Geldart D particles with the simulations containing Geldart B particles, the tendencies are similar, especially for the low inlet velocities.For the large inlet velocities, the deviation for Geldart D particles is more exaggerated than with the Geldart B particles.

Comparison of sub-grid drag models
Heterogeneous drag models tend to be ad hoc modified, meaning they can be restricted to the simulation system they were generated with.Hence, the following requirements are established.The modified drag model should decrease the observed deviations of normalized particle height for coarse grain simulations.The model should work for different particle types (Geldart class B and D particles are tested), different fluidization regimes (two inlet velocities are investigated), and different system-to-grain diameter ratios (three different system diameters are analyzed).

Effect of particle type
Figure 7 shows the normalized particle height characteristics for the heterogeneous drag laws compared to the DEM benchmark for Geldart class D and B particles using the large inlet velocity and a coarse grain factor of l = 10 for the wide system of N p = 100.The unscaled benchmark results are shown on the left for comparison.The distribution of average normalized particle height over time (t = 5−20 s) is summarized in a boxplot, hence, showing the minimum, maximum, median, and quartiles of the normalized particle height in a concise way.Outliers are shown as red symbols.An outlier is a single data point that is larger than 1.5 times the value of the upper quartile (Q3) or 1.5 times lower than the value of the lower quartile (Q1).Therefore, the boxplots not only give information about average results but also about the dynamics of the fluidized bed system.For the benchmark, the distribution of particle height shows a symmetrical trend with a similar distribution of values below and above the median.For Geldart D, more outliers are observed as for the Geldart B simulation.These tendencies are also observed using the Gidaspow model with l = 10 but with a general trend of overestimating normalized particle height for both Geldart classes.In contrast, the normalized particle height is underestimated by the heterogeneous drag models.The structure-based drag law by Wang et al. reduced the median for both Geldart types to a similar height, which yields in a large error for Geldart B and good results for Geldart D particles.However, the distribution of particle height is much larger when using the drag law by Wang et al. compared to the unscaled benchmark simulation.
Table 2 Deviation of the median of the normalized particle heights compared to the benchmark cases at the high gas inlet velocity (for Geldart class D: u in = 5.2u mf and for class B: u in = 14.2u mf ) in the wide system N p = 100 Drag model Gidaspow [47] (%) Wang et al. [26] (%) EMMS/Matrix [27,28] (%) EMMS/s-s [29,30] 2, where the relative deviation of the median value of the normalized particle height compared to the unscaled benchmark simulation as reference is given.The relative deviation is defined as where x refers to the analysis variable, in this case, the median of the normalized bed height.The effect of the drag models on the fluidization of the particles is also exemplarily depicted in Fig. 8, where snapshots of the particle distribution at a simulation time of 10 s are presented for Geldart class D. The simulation with the Wang law shows only a slight elevation of the particles.The EMMS group has a similar tendency as well as the filtered drag laws.Moreover, the time-averaged axial void fraction profiles are shown in Fig. 9.The void fraction profiles for Geldart class D particles are given in Fig. 9a.Here, it is also visible that the heterogeneous drag models all reduce drag a similar amount, which yields in a slight underprediction of the bed height, whereas the homogeneous drag model by Gidaspow leads to overprediction.For Geldart B, shown in Fig. 9b, the

Effect of fluidization regime
In fluidized beds, the fluidization regime plays an important role.While in the bubbling regime, gas bubbles are formed regularly, accompanied by large fluctuations in pressure and bed height.Higher slip velocities are reached in the slugging regime, and bubbles transform into large voids, leading to even larger pressure fluctuations [1,56].In fluidized beds, the transition from one fluidization regime to another is often non-smooth and difficult to determine before conducting experiments or simulations.Thus, the heterogeneous drag model must apply to various regimes.
In Fig. 10, the normalized particle heights are shown for the sub-grid drag models and Gidaspow drag law, here exem- plary shown for the Geldart B type particle using a coarse grain factor of 10.The benchmark results are also depicted for comparison.The results of the simulations for the high inlet velocity of 14u mf show a large underestimation of normalized particle height except for the drag law by  3, where the deviation of the median compared to the unscaled benchmark simulations is shown.

Influence of system width
At last, the influence of system width is analyzed.As demonstrated in Sect.4.1, there is a dependency on the accuracy in predicting normalized particle height on N .At lower ratios, a larger deviation to the unscaled simulation is found.The subgrid drag models must reduce the drag accordingly within the narrower systems.
Figure 11 shows the normalized particle heights for the three system widths with Geldart type B particles at the large inlet velocity for a coarse grain factor of l = 8.Note that l = 8 is the maximum coarse grain factor to ensure a minimum of 3 CFD cells in radial direction in the narrow system N p = 50.
For the benchmark cases, the behavior in the three systems shows no significant effect.For l = 8 and the Gidaspow model, the median and the distribution of average particle height increases from N p = 100 to N p = 75.From N p = 75 to N p = 50, the median stays constant.However, the distribution of particle height is much wider for N p = 50, reaching outliers that are larger than nine times the initial bed height.This indicates different behavior of the particle movement in the narrow system with large coarse grain factors.The drag law by Wang et al. reduces the drag strongly for the system width of N p = 100 and N p = 75.For the narrow width of N p = 50, the median of the normalized particle height is closer to the benchmark system, but a wider distribution of the bed height is observed.The median of the normalized particle height is reduced with both EMMS models for the narrower system width.However, bed height is largely underestimated for both laws.The EMMS/Matrix model predicts a similar normalized particle height for N p = 75 and N p = 50, indicating that drag reduction is stronger in the narrower system.The model by Lu et al. also gives good agreement for the systems with a system diameter of N p = 75 and N p = 100.All drag models fail to depict the behavior in the narrow system.Explana-  tions on why the narrow system could act differently were given in Sect.4.1, and they can also be the reason that the behavior cannot be corrected with sub-grid drag laws and other methods might be more applicable, such as interpolation methods to reduce the cell-to-particle ratio.The results are summarized in Table 4, where the deviation of the median compared to the unscaled benchmark simulations is given.

Section summary
The effect of the influence of Geldart type, fluidization regime and system width on the normalized particle height demonstrates clearly the different behavior of the sub-grid drag laws applied here, which for most cases for Geldart class B yields in a strong underprediction of normalized particle height.For the analyzed cases, the filtered drag model by The underprediction can also be explained, when comparing the length scales of the simulation systems, the models were established or verified with.Table 7 in Sect.A.3 summarizes the ratio of the CFD cell length to the characteristic cluster length scale L c .In the simulation system used here, the ratio reaches maximum values of 1.36 for Geldart class B and 0.48 for Geldart class D, while for the other studies, the ratios are larger.This is an indication that the models are more applicable when larger length scale ratios are used, such as for larger simulation systems or for smaller particle types.

Influence of coarse grain number
In Looking at the deviations of the median and standard deviation, the results emphasize that the behavior of the drag models also depends on the coarse grain factor.For a coarse Table 4 Deviation of the median of the normalized particle heights compared to the benchmark cases in the three system width for Geldart class B at the large inlet velocity of u in = 14.2u mf for l = 8 Drag model Gidaspow [47] (%) Wang et al. [26] (%) EMMS/Matrix [27,28] (%) EMMS/s-s [29,30]  grain factor of l = 6, the error for using the Gidaspow drag law is generally low, so applying sub-grid drag models, in this case, leads to a larger error.However, for l = 8 and l = 10, the error can be reduced with some sub-grid laws.

Conclusions
In this work, different heterogeneous drag models were analyzed toward their applicability to incorporate sub-grid behavior.Firstly, the results for the drag law by Gidaspow were analyzed.Different behavior of Geldart class D and B was observed, when predicting the normalized particle height with decreasing ratio N .The deviation to the benchmark case is larger at the large gas inlet velocity for both Geldart types, but more pronounced with the class D particles.The increase of particle height depends not only on the ratio N , but also on the absolute value of system width as the investigated narrow width (N p = 50) showed a larger deviation at similar ratios N .
For the simulations with large coarse grain factors of l = 6−10 and subsequently low ratios N , sub-grid drag laws were applied.Their performances in predicting particle height using configurations of different particles type, gas inlet velocity, and system width were analyzed.Looking at Geldart types, the simulations the sub-grid drag models obtained different results, which for the Geldart class D the error was less pronounced and for class B particles, the normalized height strongly underestimated except for the sub-grid law by Lu et al.Looking at the two inlet veloc-ities, a lower error at the low inlet velocity of u in = 4.2u mf , except for the drag law by Lu et al. was observed compared to results at the large inlet velocity of u in = 14.2u mf .The simulation results obtained with the narrowest system were not improved by any drag law, either due to a limiting ratio N or due to an insufficient number of CFD cells.The latter can be improved by using interpolation mapping methods to calculate void fractions and slip velocities, which will need further evaluation.Lastly, the results of the median and standard deviation of the normalized particle height were analyzed for the sub-grid drag laws with the lowest deviations for different coarse grain factors, which emphasized the different behavior of the sub-grid drag laws with different coarse grain factors.
The predictions not only vary depending on the Geldart type, inlet velocity, and system width, but also on the coarse grain factor.If a drag law is chosen, that does not fit the system's configurations, it can lead to an overcorrection of the problem and might increase the error.This work also shows that there is not one drag law that can be applied for all fluid regimes, Geldart types, and system configurations.Depending on the type of simulation, the sub-grid drag laws are a useful tool within the CG-DEM-CFD approach but always need some evaluation and comparison to other drag models.As mentioned by [57], the physical phenomena causing the micro-, meso-, and macroscale effects in fluidized beds are not fully understood yet, so a general drag formulation needs far more research.
where F h p (t) is the total hydrodynamic force acting on the particle, and F c p (t) is the contact force resulting from particle-particle or wall-particle contacts.T stands for the sum of all torques acting on the particle.In this study, the contact force is calculated using the linearspring-dashpot model.The contact force is the sum of a normal force and a tangential force, and it is calculated for each particle pair or particle-wall pair i j such as with k n as the normal stiffness, k t as the tangential stiffness.δ i j is the overlap between the particle-particle or particle-wall pair.n i j and t i j are the normal and tangential unit vectors, γ n and γ t are normal and tangential damping coefficients.v n,i j and v t,i j is the relative normal and relative tangential velocities.The hydrodynamic force consists of the pressure gradient force and the drag force: The fluid flow is solved with the gas-phase continuity and momentum conservation as follows: and with p f as the gas-phase pressure, τ f as the gas-phase stress tensor, and I f p as the momentum transfer term between the gas and the particles or solid phase p.

Drag law by Wang et al. [26]
The drag occurring in the dense phase is calculated by using the homogeneous drag model with the dense phase parameters: with ε f ,d as the void fraction and ε p,d the solid volume fraction in the dense phase.The Reynolds number is defined as: with U d as the superficial gas velocity of the dense phase, the is void fraction defined as: With these, the drag force within the emulsion phase can be calculated with: For any computational grid cell, a varying part of it is occupied by the dense phase, and the drag force within the cell amounts to: where ε p ε p,d is the emulsion phase fraction in the cell.Now, only ε p,d (ε f ,d = 1 − ε p,d ) and U d need to be determined.They are derived from experimental correlations.ε p,d can be calculated with the following equations: with and where ε mf is the voidage at minimum fluidization.With ε f ,d , U d can be calculated with the following relation: with v t being the particle terminal velocity and with the exponent n being: The drag coefficient for each computational cell can now be calculated with:

EMMS/Matrix [27, 28]
The EMMS/Matrix heterogeneity index is read into the software by bilinear interpolation.The heterogeneity index as a function of void fraction and slip velocity is given in Fig. 13.

EMMS/steady-state [29, 30]
In the following, the equations used for the EMMS/steadystate model are stated.Equation 31 is used for simulations with Geldart D particles at an inlet velocity of u in = 2.5u mf , Eq. 32 for simulations with Geldart D particles at an inlet velocity of u in = 5.2u mf .Equation 33 is used for simulations with Geldart B particles at an inlet velocity of u in = 4.2u mf , Eq. 34 for simulations with Geldart B particles at an inlet velocity of u in = 14.2u mf . with with

Drag law by Radl et al. [31]
β β is calculated as: , ε p h(ε p ) c corr (l, ε p ) (35) with and c corr = exp(0.05(l− 1)) (37) with a(ε p ) and h(ε p ) being piecewise continuous algebraic expressions presented in Tables 5 and 6.The function c corr was introduced to account for particle coarsening and has been validated for coarse grain factors up to l = 8.
The characteristic cluster length scale L c is given by: Fr n p (38) with v t being the terminal particle velocity and Fr being the particle Froude number.The exponent n was determined to be − 2  (46) and f as the filter size: where V cell is the volume of the fluid cell.The regressed coefficients are:

A.3 Comparison of length scales in the models
The applied models in this study were all verified based on simulations of fluidized beds and compared either to experimental or fine grid simulation data.For the coarse simulations in these studies, the CFD cell length dx and the characteristic cluster length scale L c from Eq. 38 are summarized in Table 7.For comparison, the length scales used in this study for both particle classes (Geldart class B and class D) for a coarse grain factor of l = 10 are also shown.

Fig. 1
Fig. 1 Heterogeneity index trends with changing void fraction at a constant slip velocity of u slip = 0.98 m/s for the applied drag laws

4 . 1 .
The system width is reduced from D = 100 to D = 75 mm and D = 50 mm for Geldart D particles and from D = 50 to D = 37.5 mm and D = 25 mm for Geldart B particles.

Fig. 3
Fig. 3 Normalized particle height for coarse grain factors between 2 and 10 for Geldart B particles (a) and Geldart D particles (b) at the two gas inlet velocities for the three system widths.The symbols mark the different system widths, unfilled at the low gas inlet velocity and filled at the high gas inlet velocity.The different colors of the symbols represent l

Fig. 4
Fig. 4 Power spectral density plots for the unscaled benchmark cases in the wide system (N p = 100) for Geldart class B (a) and D (b) particles.The blue line shows the result at the low gas inlet velocity, and the black line at the high gas inlet velocity

Fig. 5
Fig. 5 Snapshots taken at 10 s of a Geldart B type simulation with an inlet velocity of 14u mf in the wide system of N p = 100.a shows the benchmark case (l = 1), b for l = 4, c for l = 10.On the right-hand side, the time-averaged void fraction profiles for these cases over the system height are shown

Fig. 7
Fig. 7 Boxplots of the normalized particle height for Geldart B and Geldart D particles.Here, exemplary shown at the high gas inlet velocity (for Geldart class D: u in = 5.2u mf and for class B: u in = 14.2u mf ) in the wide system (N p = 100) for a coarse grain factor of l = 10

Fig. 8
Fig. 8 Snapshots taken at 10 s of Geldart class D simulations with an inlet velocity of 5.2u mf in the wide system of N p = 100 for the applied drag models at coarse grain factor of l = 10, a shows reference results for the unscaled benchmark case, b for the homogeneous drag law by Gidaspow, c for the EMMS/steady-state model and d for the EMMS/Matrix model.The snapshots for regression-based models by Radl, Lu, and the structure-based law by Wang are shown in (e), (f), and (g), respectively

Fig. 9 Fig. 10
Fig. 9 Time-averaged void fraction profiles along the relative system height for Geldart class D and B particles in the wide system (N p = 100) at the large inlet velocity (D: u in = 5.2u mf , class B: u in = 14.2u mf ) for a coarse grain factor of l = 10.The unscaled benchmark is given as a reference

Fig. 11
Fig. 11 Normalized particle height for Geldart B type particles at the three different system widths for l = 8 at the large inlet velocity (u in = 14u mf ) this subsection, the drag law by Lu et al. is investigated on its accuracy with varying coarse grain factors of l = 6, 8, 10 for Geldart class B and the drag laws by Lu et al., Wang et al., and the EMMS/Matrix drag law are investigated for Geldart class D and compared to the Gidaspow drag law.As an example, the deviations of the median and standard deviation of the normalized particle height in the system width of N p = 100 are analyzed.The wider system is chosen to avoid the error that is caused by the narrower system widths and not the reduced ratio N .In Fig.12a and b, the results for the deviation of the median particle height for both Geldart types are shown.For Geldart type D, the deviation for the drag law by Gidaspow increases with increasing coarse grain number.The heterogeneous drag laws show the opposite trend, leading to a smaller deviation with a larger coarse grain number, which results in reduced deviation for coarse grain factors of l = 8 and l = 10.For Geldart type B, the trends are similar.When looking at the deviation of the standard deviation of the particle height signal shown in Fig.12c, the results for the homogeneous Gidaspow drag law increase until l = 8 and then decrease, similar behavior is observed for EMMS/Matrix model and the drag law by Lu et al.The drag law by Wang shows a large deviation for all cases for the standard deviation of particle height.For the Geldart B particles shown in Fig. 12d, the deviations are increasing with increasing coarse grain number for the homogeneous model.The drag law by Lu et al. shows a reduced error except for l = 8.

1 Fig. 12
Fig. 12 Absolute values of the relative deviation of the median and standard deviation of particle height to the unscaled Benchmark case for Geldart class B particles (a, b) and Geldart class D (c, d) at the large gas inlet velocity (for Geldart class D: u in = 5.2u mf and for class B: u in = 14.2u mf ) for different coarse grain factors l in the wide system of N p = 100

Table 1
Simulation parameters 5 /2.25 × 10 5 /1 × 10 5 4 × 10 5 /2.25 × 10 5 /1 × 10 5 is always below 1.The drag law by Lu et al. decreases drag significantly for a void fraction of ε f = 0.42 and ε f = 0.99, but only minor drag reduction is conducted for a void fraction range of ε f = 0.6 − 0.8.No correction is performed for void fractions below 0.42 and above 0.99, similarly to the EMMS/Matrix model.Note that the heterogeneity index calculated with the drag laws by Radl et al., Lu et al., and the EMMS/Matrix law also depend on the slip velocity.Thus, the actual behavior can be different based on the reached slip velocities in the system.
[35]Radl et al.[31](%) Lu et al.[35](%)The EMMS/Matrix law reduces the normalized particle height for both particle types.The prediction of the median for Geldart class D particles is similar to the benchmark case.The prediction for Geldart class B particles is much lower and reaches only twice the initial bed height.Looking at the EMMS/steady-state model, for Geldart type B, the bed height is reduced from four to two times the initial bed height.For Geldart type D, the effect is less pronounced.However, less outliers on the top are observed compared to the EMMS/Matrix model, which is less in agreement with the benchmark case.Depending on the Geldart class, the model uses two different correlations.The correlation for Geldart D type represents the particle dynamics better in this case.The regression-based Radl model also underestimates the results for both Geldart classes.The distribution for Geldart D particles gets wider, whereas for Geldart B particles, more outliers above the upper whisker are observed.The regression-based model by Lu et al. represents the results well, but the distribution becomes slightly wider.The wider distribution of particles in the bed can be observed for the homogeneous Gidaspow drag law and the heterogeneous drag laws by Radl et al., Lu et al. and Wang et al., which can be the statistical influence because the total number of particles is decreased significantly.For the EMMS models with Geldart class B particles, the movement of particles is generally much lower as the models generate the strongest drag reduction.The results are also summarized in Table Lu et al.,as already discussed in Sect.4.2.1.Looking at the low inlet velocity of 4.2u mf , the Gidaspow model overpredicts the bed height only slightly, indicating that the homogeneous drag law is also applicable in low-fluidization regimes.However, the filtered drag models must incorporate this and only apply slight modifications to the drag model.The drag law by Wang et al. leads to a large overprediction of bed height in this case.The EMMS/Matrix drag law and drag law byRadl et al.show almost no movement of the particles, a narrow distribution, with only a slight increase in height compared to the initial bed height.The EMMS/steady-state model gives a slight decrease compared to the Gidaspow model at low inlet velocities, indicating that it is suitable for lower gas velocities for Geldart class B in the bubbling regime.The results are also summarized in Table

Table 3
Deviation of the median of the normalized particle heights compared to the benchmark cases at the two gas inlet velocities for Geldart class B in the wide system N [35]30] Drag model Gidaspow[47](%) Wang et al.[26](%) EMMS/Matrix[27,28](%) EMMS/s-s[29,30](%) Radl et al.[31](%) Lu et al.[35](%) Lu et al. performedwell except for the narrow system width N p = 50, where the deviation could not be corrected.The drag law by Wang et al. gave much wider distributions of normalized bed height.However, it could reduce the deviation for the median looking at the Geldart class D example.Also, the EMMS/Matrix drag law predicted the Geldart D class normalized particle height well.For the low inlet velocity, a drag correction is not necessary.The EMMS/steady-state model is developed for continuous fluidized bed (CFB) and bubbling regime, where it performed well at the low inlet velocity.Similar reasons explain the deviations from the drag laws byRadl etal.and Wang et al.Moreover, the model by Wang et al. does not incorporate macroscale markers, e.g., gas inlet velocity or system size.

Table 7
CFD cell size and characteristic cluster length scale of the particle systems in literature.These simulations were used to verify the applied drag models.For comparison, the length scales for the system in this study are also shown for Geldart B and Geldart D particles for a coarse grain factor of l = 10.The terminal velocity is calculated with Eq. 46, the characteristic cluster length scale with Eq. 38 Wang et al. [26] TFM Experimental data based on Werther & Wein 4.53 5.00 × 10 −2 4.94 × 10 −3 10.12 Wang et al. [26] TFM Experimental data based on Johnson et al. 24.68 2.26 × 10 −1 2.69 × 10 −2 8.42 + a 2 ε p + a 3 ε 2 p + a 4 ε 3 p + a 5 ε 4 p )(1 − e a 19 ε p ) 1 + e a 20 (ε p −0.55) P = a 14 + a 15 ε p + a 16 ε 2 X = a 21 + a 22 + a 23 f * + a 24 ε p + a 25 u * with the filtered slip velocity u slip = |u f − v p | and the terminal particle velocity v t : v t = gd 2 p (ρ g − ρ p ) 18 • μ g 25 = [4.802706,12.753237, 67.982539,