A detailed gas-solid fluidized bed comparison study on CFD-DEM coarse-graining techniques

Computational Fluid Dynamics - Discrete Element Method (CFD-DEM) is a numerical tool used for detailed ﬂuidized bed studies. However, CFD-DEM is computationally expensive, leading to a restriction regarding the number of simulated particles. Therefore, coarse-graining techniques have been developed to increase the simulation scale or reduce computational requirements for CFD-DEM simulations in ﬂu-idized beds. In this work

Lagrangian fashion.The latter indicates the major drawback of the CFD-DEM approach since the tracking of each individual particle is computationally expensive.Therefore, the simulation scale is generally restricted to 10 5 to 10 8 particles (Golshan et al., 2020).It is obvious that this limits the suitability of CFD-DEM to study large gas-solid fluidized beds.The second method, MP-PIC (Andrews and O'Rourke, 1996;Snider, 2001;O'Rourke and Snider, 2010), is less computationally expensive to study large systems.It uses a so-called parcel approach where each particle represents a number of real particles.The collisions between particles are represented via particle stresses.A point of concern is the reported relatively large disagreement with results obtained from more accurate CFD-DEM method (Verma and Padding, 2020).The third method is the TFM (Ding and Gidaspow, 1990;Kuipers et al., 1992;Verma et al., 2013) where both the gas and particle phase are treated as continuous media using the kinetic theory of granular flow (KTGF) to describe the particle-particle interactions.The TFM is fast but there are several drawbacks with respect to dense granular flows, treatment of polydispersity and tractable shapes of particles (Van Sint Annaland et al., 2009;Verma and Padding, 2020).
A coarse-grained CFD-DEM technique is an approach to circumvent the restriction with respect to the simulated number of particles in CFD-DEM avoiding the limitations of MP-PIC and the TFM.Coarse-graining is based on the substitution of the actual number of DEM particles by a smaller number of larger particles while applying scaling laws.In the literature, several coarse-graining scaling laws were presented as recently reviewed by Di Renzo et al. (2021).One of the most extensively researched methods was introduced by Sakai and Koshizuka (2009) for the usage of pneumatic conveying.The method was subsequently applied to fluidized beds, where Sakai et al. (2010) studied the effect of the coarse-graining method with respect to the bed height, bubble size and speedup.The method was validated by Sakai et al. (2014) by comparing simulation results with the experimentally obtained pressure drop and bed height.More recently, spout fluidized beds were studied (Takabatake et al., 2018;Mori et al., 2019).
Another scaling approach was recently proposed by Mu et al. (2020), applying the so-called exact scaling laws of Feng and Owen (2014).Mu et al. (2020) performed simulations of a pseudo-2D circulating fluidized bed and compared the obtained solids holdup, solids mass flux and solids cluster characteristics with experimental data of Carlos Varas et al. (2017).Benyahia and Galvin (2010) simulated shear and riser flow using a coarse-graining technique.The riser flow coarse-grained simulations were compared via time-averaged lateral profiles in solids volume fractions, gas and solids velocities and the granular temperature.Their results were considerably deviating from the non-coarse grained system results, leading to a proposed change in restitution coefficient.This method was the basis for the extensions to reactive systems (Lu et al., 2016), heat transfer (Lu et al., 2017b) and polydisperse gas-solid flow (Lu et al., 2018).Radl et al. (2011) used a similar approach for changing the collisional dissipation in a granular jet study.
In the literature, various other coarse-graining methodologies are available, see the review of Di Renzo et al. (2021) and recent contributions (Zhou and Zhao, 2021;Kanjilal and Schneiderbauer, 2021;Yu et al., 2021;Zhou et al., 2022;Chu et al., 2022).This raises the question which scaling methodology performs best with respect to the true system behavior.Another issue encountered in coarse-graining is the problematic two-way gas-particle coupling as pointed out by Di Renzo et al. (2021).An increase in particle size violates the CFD-DEM principle where the particle size should be smaller than the grid size (Peng et al., 2014).This issue was addressed by Radl and Sundaresan (2014, 2017a, 2020) where sub-grid corrections were applied or the grid size was scaled with the coarse-graining ratio.Grid size scaling can become problematic at larger coarse-graining ratios as the gas flow cannot be accurately resolved.Well before coarse-graining studies were commenced, two-way coupling complications in Euler-Lagrange models were

Roman letters m i
Particle mass (kg) F contact;i Particle contact force (N)  Kitagawa et al. (2001).Discontinuous two-way coupling or grid refinement in combination with local coupling (such as volume-weighted techniques (Van der Hoef et al., 2006)) can lead to large unphysical deviating pressures and velocities.Furthermore, it should be noted that upon grid refinement griddependent solutions are resulting.In order to overcome this dependency, Link et al. (2005) used a porous cube approach where the cube volume is dependent on the particle diameter.Besides, smooth mapping functions such as a Gaussian (Peskin, 1977), sine-wave (Kitagawa et al., 2001) or polynomial (Deen et al., 2004) were proposed in the literature.The fourth-order polynomial of Deen et al. (2004) is a cheap computational two-way coupling method including a dependency on the particle diameter and not on the grid cell size.This method was previously applied for liquid-gas Euler-Lagrange simulations, see Lau et al. (2014Lau et al. ( , 2020)).The aim of this work is twofold, we first critically compare the coarse-graining scaling laws of Sakai and Koshizuka (2009) and Mu et al. (2020).The idea of this contribution is to illustrate the necessity of an in-depth analysis in order to adequately use coarsegraining techniques to gas-solid fluidized beds without losing the gas-solid contacting mechanisms of the real granular material.Secondly, we investigate the usage of a grid-independent fourthorder polynomial proposed by Deen et al. (2004) to avoid the aforementioned discussed two-way coupling complications in coarsegraining simulations.In this work, we use a newly developed CFD-DEM model by combining the CFD framework FoxBerry (Kamath et al., 2020) with the DEM code MercuryDPM (Weinhart et al., 2020).The model validation can be found in the supplementary material.
Our paper is organized as follows.In the next section, the coarse-graining CFD-DEM modeling methods are described.Thereafter, the simulation setup and results of the comparison study are presented in Section 3. In Section 4 we discuss the usage of the smooth fourth-order polynomial function for coarse-graining purposes.A final conclusion regarding both the scaling law comparison study and the usage of the polynomial two-way coupling method is given in Section 5.

General equations
The gas phase is described by the continuity equation and the volume-averaged Navier-Stokes equations given by: where the gas phase density is calculated according to the ideal gas law.s g is the gas phase stress tensor using the general Newtonian form: S represents the momentum source term for the gas-solids interaction, given by Eq. 4, where D stands for the usage of a volumeweighted or polynomial distribution function using the particle position r i relative to the computational cell r.For more detailed information about the Euler-Lagrange coupling methods, see Deen et al. (2004Deen et al. ( , 2006)).The momentum source term contribution is implemented via a semi-implicit method (Kuipers et al., 1993;Verma et al., 2013;Patil et al., 2014).
The interface momentum exchange coefficient, b, describes the drag of the gas phase acting on a particle, which can be calculated via a drag force closure.In this work, the drag correlation proposed by (Beetstra et al., 2007) is used: with The translational and rotational movement of each individual particle can be obtained by solving Newton's equations captured respectively Eq. 7 and 8: where v i is the particle velocity.The forces on the right-hand side for the translational motion (Eq.7) are respectively due to pressure gradient, drag, gravity and contact forces due to particle-particle or particle-wall collisions.The contact forces are captured by a softsphere model originally developed by Cundall and Strack (1979).
T i is the torque, I i the moment of inertia and x i the rotational velocity.
In both scaling models, the total particle mass and volume of the coarse-grained system are equal to the original system (coarse-graining ratio equal to 1).A coarse-grained particle represents l 3 original particles, where l stands for the coarse-graining ratio.This means that the particle diameter is scaled with a factor l. Therefore, the number of particles in a coarse-grained system will reduce by a factor l À3 .In the subsequent subsections, the mathematical description of the two different scaling laws is given.
2.2.Scaling law of Sakai and Koshizuka (2009) In this scaling law, the translational motion of the coarsegrained particle (indicated by c; j) is assumed to be an average of the original particles by multiplying all contributions by a factor l 3 as presented in Eq. 9.By using a coarse-graining ratio equal to 1, it can easily be shown that the original balance, given in Eq. 7, is obtained.The scaling with l is such that the kinematics of the coarse-grained particle is similar to that of the original (noncoarse grained) particle.This means that both the velocity and the acceleration remain the same irrespective of l.Since m c;j scales with l 3 all forces on the right-hand side should scale with l 3 as well.
m c;j dv c;j dt The applied Beetstra drag model uses the Stokes drag consisting of the particle diameter which scales with l.Therefore, an additional factor l 2 is applied in order to scale the drag force with a factor l 3 .
The Reynolds number inside the drag correlation is unchanged by using the original particle diameter: The contact force in Eq. 9 uses the linear spring-dashpot model, where the normal component is scaled with a factor l 3 .k n ; d n and g n are respectively the normal stiffness, the overlap and the damping coefficient.In coarse-grained systems, the original physical properties such as the stiffness and damping coefficient can be used.
In a similar fashion, the tangential component of the contact force is given by: where l 0 is the friction coefficient, k t the tangential stiffness, d t the overlap and g t the damping coefficient.Using the rotational motion shown in Eq. 8, the coarse-grained rotational motion will become: This scaling is such that the surface velocity of the particle remains the same independent of l.This is also important for the ratio between normal and tangential contact forces.

Scaling law of Mu et al. (2020)
This scaling method aims to reduce the number of discrete particles while maintaining constant dimensionless groups between the original and the coarse-grained system.Mu et al. (2020) used the translational motion (see Eq. 7) of the particles in the vertical direction (i.e.z-direction) as a starting point.The forces due to the pressure gradient and the gravity can be combined, leading to: Besides, Eq. 7 can be rewritten by splitting the contact force into the normal and tangential contributions and dividing by m i g.This results in the following balance, where each group should remain constant in order to maintain the same hydrodynamic behavior. 1 The first term on the right-hand side implies that the ratio between the pressure and gravity force remains constant, leading to: The second group implies that the ratio between the drag and gravity force should remain constant.This ratio was expressed in terms of the particle Reynolds number, Re (see Eq. 6), the Archimedes number Ar and N AR : where: By keeping the Reynolds number, the voidage, N AR and b (which is in principle a function of the Reynolds number and the voidage) constant, the correct scaling is achieved.The third term consists of a constant ratio between the normal contact and gravity force where the overlap and the normal velocity are normalized using respectively the particle diameter and impact velocity (v 0 ): Besides, the third term implies that the ratio between the normal and tangential contact forces remains constant.This is guaranteed in the linear spring-dashpot model.Summarizing, dimensionless groups are obtained that should remain constant in the case of coarse-grained systems.The dimensionless groups are rewritten into a more pragmatic form: In these dimensionless groups, the following parameters will be kept constant: ; g ; v 0 ; q p ; q g , while the seven other parameters will change in order to ensure proper scaling: Using the dimensionless groups presented in Eq. 20, the parameters are scaled according to:

Scaling law comparison
In order to study the different scaling methods, a 3D fluidized bed column with dimensions 0.3 Â 0.06 Â 0.036 m (height Â width Â depth) was simulated.In the base case, 720 000 glass beads were used to maintain a 3D configuration when a large coarse-graining ratio is applied.The scaling laws were studied using coarsegraining ratios 1.5, 2, 2.5, 3 and 3.5.The resulting numbers of particles were equal to 213 333, 90 000, 46 080, 26 667 and 16 793 respectively.The simulation parameters are listed in Table 1.See Deen et al. (2007) for more details about calculating the spring stiffness and damping coefficient.The scaled simulation parameter values are presented in the supplementary material Sections 2 and 3.In order to avoid the aforementioned two-way coupling limitations, a relatively large grid size was applied.A schematic representation of the setup is shown in Fig. 1.In the initialization step, the particles were placed in a lattice structure.The boundary conditions for the gas flow were set as a uniform superficial air velocity equal of 0.50 m/s at the bottom, no-slip boundary conditions for the side walls and a fixed pressure boundary condition of 1 atm at the top of the domain.An isothermal temperature of 20 C was set that was used for calculating the air density via the ideal gas law.A volume-weighted mapping method was employed for the two-way gas-particle coupling.
As discussed in Section 2, the Beetstra drag model was applied in the simulations (Beetstra et al., 2007).However, Sakai and Koshizuka (2009) originally designed the scaling methodology using the Gidaspow (1994) drag model.Therefore, both drag models were investigated in this study but negligible differences between the drag models regarding the scaling methods were obtained.For sake of clarity, therefore only the results using the Beetstra et al. (2007) drag model are shown in this work.

Instantaneous particle configurations
Fig. 2 presents instantaneous snapshots of particles near the center plane of the column for the base case simulation (i.e.no applied scaling).The base case simulation shows fluidizing behavior including the well-known bubble formation, propagation and eruption.The bubble eruption causes clearly visible chaotic solids movement at the top of the bed.Figs. 3 and 4 show the instantaneous particle configurations of the coarse-graining ratios equal to 2 and 3 taken at 2.5, 5 and 7 s of simulation time for respectively the scaling laws of Mu et al. (2020) and Sakai and Koshizuka (2009).As a result of coarse-graining, fewer but bigger particles are present in the systems with larger coarse-grain ratios.Both scaling laws feature fluidizing behavior as seen in the instantaneous particle configurations.It is observed that the flow patterns of the original system are captured in more detail compared to the coarse-grained systems.This is due to the coarse-grained representation of N number of original particles.The lost details are especially noticable in chaotic regions of the bed such as the eruption zone.It should be noted that the coarse-graining scaling law of Mu et al. (2020) also results in a remarkably large bed expansion as shown in Fig. 3B.The shown instantaneous image profiles could differ since bubble formation, propagation and eruption is a very dynamic process.However, this large bed expansion is not observed in the base case simulations.Therefore, this needs to be investigated further.The obtained particle configurations using the method of Sakai and Koshizuka (2009) are quite similar to the base case.However, it is hard to draw quantitative conclusions from individual snapshots and consequently, a more in-depth analysis is needed.

Time-averaged pressure drop
In order to avoid startup effects, the first second is discarded for the subsequent analyses.The pressure drop is a key parameter for fluidized beds.Once the minimum fluidization velocity is attained, the pressure drop equals the bed weight.However, this is dependent on the gravitational acceleration acting on the bed.The gravitational acceleration is scaled in the method of Mu et al. (2020) as shown in Eq. 21.Therefore, the pressure drop is made dimensionless using the column cross-sectional area, A, the applied gravitational acceleration, g, and the solids mass.This results in an average dimensionless pressure drop (DPA/(mg)) presented in Fig. 5A.The dimensionless pressure drop is accurately captured by the coarse-grained simulations and a maximum deviation equal to 1.2% compared to the base case is found.
In fluidized beds, the pressure drop standard deviation is strongly correlated with the energy dissipation due to the particle collisions as shown by Li andKuipers (2007, 2013).Typically, an increase in energy dissipation leads to more heterogeneity and thus a larger pressure drop standard deviation.Fig. 5B presents the dimensionless pressure drop standard deviation for both scaling laws.It can clearly be seen that the standard deviation remains constant until the coarse-graining ratio is equal to 2 in the case of the method of Mu et al. ( 2020), while it slightly decreases for the  value 2.5 where-after it remains constant again, resulting in a maximum deviation equal to 10%.The methodology of Sakai and Koshizuka (2009) results in a reduced standard deviation up to 52% of the original value.An explanation for the reduced level of pressure fluctuations is the coarse-grained representation of N number of original particles.In the original system, the motion, collisions and the gas-solids interaction are described for each individual particle, while in a coarse-graining system these components are averaged into a representation of the original system.This averaged representation is by construction more homogeneous, explaining the lower pressure drop standard deviation.However, the direct influence on the resulting gas-solids interaction should be investigated in further detail by analyzing the solids volume fraction profiles and bubble characteristics and will be reported subsequently.

Time-averaged solids volume fraction distributions
The time-averaged solids volume fraction distribution versus the axial position is presented in Fig. 6.In a coarse-graining methodology, it is critical to maintain the solids holdup distribution of the original system as this is the main indicator of the bed height.The scaling method of Mu et al. (2020) shows a decrease in solids volume fraction as the coarse-graining ratio increases.As a result, the bed height is increased.Li and Kuipers (2007) investigated in detail the effect of particle-particle collisions and gas-particle interactions on the flow pattern in dense gas-fluidized beds.An increase in particle-particle dissipation while keeping the gas-particle interaction constant leads to a more heterogeneous flow structure including larger but less frequently appearing bubbles.On the other hand, the method of Sakai and Koshizuka (2009) shows more correspondence with the original system.However, a slight change is observed for the coarse-grain ratios 3 and 3.5, especially in the eruption zone a more step-wise solids volume fraction profile is observed.In this zone, the solids movement is very chaotic.By applying coarse-graining, this chaotic movement will be less pronounced and consequently a more homogeneous flow is obtained which results in a more step-wise solids volume fraction profile.However, it should be noted that the bed height is independent of the level of coarse-graining.Using the snapshots shown in Fig. 4 for illustrative purposes, we clearly observe that the chaotic behavior is indeed suppressed when one increases the coarse-graining ratio.Another reason could be the particle diameter which approaches the grid cell size at larger coarsegraining ratios.This could become problematic and will be discussed in more detail in Section 4.
A cross-sectional average solids volume fraction at a height of 6 cm is analyzed to further quantify the differences between both scaling methods.In fluidized beds, the cross-sectional average solids volume fraction follows typically a U-shaped profile where dense regions are formed near the walls while in the center more bubbles are observed leading to a lower average solids volume fraction.As discussed in the axial solids volume fraction profiles, the methodology of Mu et al. (2020) leads to an increased bed expansion due to an increasingly heterogeneous flow.The bubble size becomes larger and this is substantiated by Fig. 7A where the U-shaped profile is changed as a function of the coarsegraining ratios.The 3 and 3.5 ratios result in a flat profile in the center of the column that could even indicate bubble slug formation.The scaling law of Sakai and Koshizuka (2009) produces better correspondence with the original system as shown in Fig. 7B for the coarse-graining ratios up to 2.5 for both the wall and center zone.The ratios 3 and 3.5 clearly deviate from the original system.

Time-averaged solids vertical mass flux
The time-averaged vertical solids mass flux is obtained by multiplying the cell local average solids volume fraction at a height of 6 cm with the average obtained vertical particle velocities of all particles in the cell of interest.This is captured by: In a fluidized bed, an upward solids motion is present in the center of the bed as the particles are dragged along in the wake of rising bubbles, while a downward flow exists in the wall region leading to respectively a positive or negative solids vertical mass flux.Fig. 8 presents the solids mass flux for both methods.The original system is not well-presented by the scaling methodology of Mu et al. (2020) where the solids volume flux becomes flatter with increasing coarse-graining ratio.The obtained time-averaged solids volume fraction profiles contribute to this trend.The method of Sakai and Koshizuka (2009) represents the solids motion well especially for the ratios up to 2. At a ratio of 2.5, the correspondence with the original system becomes less satisfactory.This is similar to the trend observed for the quantities presented before.

Bubble frequency and size
The bubble frequency is obtained via a fast Fourier transform (FFT) on the pressure drop signals.The dominant frequencies can be obtained from the highest values in the power spectrum.However the dominant frequency cannot be easily captured by one value, therefore a Lorentz fit was applied after which the dominant frequency is obtained from the Lorentz fit peak value.Previously the effect of an increase in dissipation with increasing coarsegraining ratio was discussed for the scaling of Mu et al. (2020).This increase results in a reduced bubble frequency up to 44% of the original base case value for the scaling law presented by Mu et al. (2020) as shown in Fig. 9.The scaling method of Sakai and Koshizuka (2009) leads to bubble frequencies that are practically independent of the level of coarse-graining since a maximum deviation equal to 0.2 Hz compared to the base case is found.
The bubble size is an important parameter in fluidized beds that needs to be well represented also in a coarse-grained simulation.The equivalent bubble diameter is obtained from a postprocessing algorithm where the bubble regions were determined via a threshold porosity value ( g ) of 0.75 where the freeboard region was excluded from the analysis.In addition, a void should have a minimum size of 2.5 times the grid cell size (i.e. 5 mm) to be classified as a bubble.The equivalent bubble diameter is obtained from: Fig. 10 displays the equivalent bubble diameter for both scaling laws.The data is cut-off at the height where the eruption zone starts for the original system.In the results presented so far, the methodology of Mu et al. (2020) showed a clear trend towards a more heterogeneous flow structure that is also observed for the equivalent bubble diameter.In the bottom section, the bubble growth is well captured but a clear discrepancy up to 1.5 times the local bubble size of the base case is observed above a height of 4 cm.At a height of 7.5 cm, the equivalent bubble diameter shows better agreement.However, it should be noted that this region can be indicated as the start of the eruption zone.The equivalent bubble diameter also substantiates the relatively flat Ushaped profiles as the bubble diameter becomes extremely large compared to the original system.Fig. 10B demonstrates a good coincidence with the original system using the method of Sakai and Koshizuka (2009) as also reported by Sakai et al. (2010).However, the equivalent bubble size is underpredicted in the bottom regions.Besides, it becomes clear that the coarse-graining ratios of 3 and 3.5 produce larger deviations as found earlier.At these ratios, bubble size becomes relatively smaller (64% of the original local base case value) and it seems that the eruption zone starts already above a height of 7 cm.This is also observed in Fig. 6B and it could also be one of the reasons for the relatively flat vertical solids mass flux presented in Fig. 8B.

Polynomial mapping function for two-way coupling
In the previous section, the scaling methodology of Sakai and Koshizuka (2009) produces a very good correspondence with the original system except for the largest two coarse-graining ratios of 3 and 3.5.Two main reasons were given for the observed differences.The first one considers the averaging of the original system, leading to a more homogeneous flow behavior by increasing the coarse-graining ratio.The second one is related to the particle to grid cell size limitation in the case of the applied volumeweighting technique.
Therefore, additional simulations were carried out where a smoothing function instead of a volume-weighing method was applied.In the case of volume-weighing, the grid size determines the length scale of the smoothing.When using a smoothing function, the length scale of the smoothing can in principle be chosen independent of the grid size.The fourth-order polynomial function of Deen et al. (2004) is applied with a mapping width of 3d p .This width is in accordance with the porous cube approach of Link et al. (2005).Please note that the mapping width linearly increases upon coarse-graining ratio via 3d c;j = 3ld p;i .A slightly smaller 3D fluidized bed column with dimensions 0.24 Â 0.06 Â 0.036 m (height Â width Â depth) was simulated.In order to show that this twoway coupling method can be applied for particle sizes exceeding the grid cell size, we reduced the Eulerian grid size and increased the maximum researched coarse-graining ratio.The grid cell size reduction could result in more resolved flow structures compared to the larger grid cell size applied in the previous section.However, this is not the focus of this work.The polynomial two-way coupling was studied for ratios 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5 and 5.The resulted number of particles were 720 000, 213 333, 90 000, 46 080, 26 667, 16 793, 11 250, 7 901 and 5 760 respectively.The simulation parameters are listed in Table 2.
Similar to the results reported in the previous section, the scaling law of Sakai and Koshizuka (2009) describes the original sys- (B) Sakai and Koshizuka (2009) Fig. 8. Solids mass flux profile at a height of 6 cm.The method of Sakai and Koshizuka (2009) represents the original data well (except for large coarse-graining ratios) while the method of Mu et al. (2020) is not able to.et al. (2020) shows a clear diminishing bubble frequency up to 44% of the original base case value while the method of Sakai and Koshizuka (2009) results in a maximum bubble frequency deviation equal to 0.2 Hz compared to the base case value.
tem well.Therefore, we will discuss only the solids volume fraction, solids mass flux profile and bubble size data as a function of coarse-grain ratio.Comparable time-averaged solids volume fractions were found as in the previous section.However, the particle size is clearly larger than the grid cell size.This was not feasible with the volume-weighted function used so far, as the grid cell size determined the length scale of the smoothing and thus the maximum allowable coarse-graining ratio.However, Fig. 11 still results in deviating time-averaged solids volume fractions for large coarse-graining ratios.Therefore, we conclude that the dominant reason for the observed discrepancies was caused by the averaging nature of the original system by applying coarse-graining scaling laws.Besides, the coarse-graining ratios equal to 4.5 and 5 indicate that a maximum permissible coarse-graining threshold is reached for the system of interest.In these two cases, the column depth to particle diameter becomes too small (depth/d p < 15), leading to a more dominant influence of the column walls.This results in pseudo-2D fluidized bed behavior.Similar to the volume-weighted case, the solids mass flux profiles presented in Fig. 12 show flattening at large coarse-graining ratios.Analogous to the discussion for the time-averaged solids volume fractions, we conclude that this behavior is mainly caused by the averaging nature of coarse-graining and not due to the grid dependency of the volume-weighting technique.Besides, we want to point out that the results for larger ratios could differ as a result of the particle diameter increase.This leads to a larger, and possibly dominant influence of the walls.The equivalent bubble diameter has been calculated using the same procedure as described in the previous section.In Fig. 10, a slightly different pattern was found compared to the results shown in Fig. 13.As all cases show these differences, it is expected that this is related to the reduced grid cell size leading to a better resolved flow field which is important to capture the bubble characteristics.Furthermore, it is clearly visible that the equivalent bubble diameter is well predicted at the employed coarse-graining ratios (differences for the coarsegraining ratios up to 3.5 are found to be a maximum of 30%), whereas in the volume-weighted two-way coupling method, the larger coarse-graining ratios (3 and 3.5) showed often large discrepancies (maximum value up to 36%).Therefore, we conclude here that the smoothing function indeed will improve the results   for larger coarse-graining ratios.This could also be substantiated by our previous discussion (see the introduction and Kitagawa et al. (2001Kitagawa et al. ( , 2004)) about two-way coupling functions and the resulting unphysical deviating pressures and velocities for discontinuous functions when the particle diameter approaches the grid cell size.However, it is still observed that the bubble size is also underpredicted at a coarse-graining ratio equal to 5 up to values of 48% of the original local bubble size.This is consistent with the previously shown figures.Therefore, as discussed, this may indicate a ratio limitation for the particular fluidized bed system of interest.
The scaling law of Sakai and Koshizuka (2009) is able to well predict the original system.However, at large coarse-graining ratios, differences were observed due to the averaging nature and the particle to grid cell size limitation in the case of a volumeweighting technique.By applying the fourth-order polynomial function, the smoothing is independent of the grid size.Therefore, the observed differences are only caused by the averaging nature of coarse-graining.This clearly improves the equivalent bubble diameter profiles at higher coarse-graining ratios.Furthermore, the particle diameter can be clearly larger than the grid cell size by applying the fourth-order polynomial function, allowing one to simulate larger coarse-graining ratios while keeping the original grid cell size.

Conclusions
In this work, the CFD-DEM coarse-graining laws of Mu et al. (2020) and Sakai and Koshizuka (2009) were critically compared in a gas-solid fluidized bed.An extensive analysis via pressure drop, solidity, solids mass flux, bubble frequency and bubble size showed that the method of Mu et al. (2020) leads to large discrepancies.On the other hand, the scaling law of Sakai and Koshizuka (2009) is able to well predict the original system.Special care needs to be taken as coarse-graining representation simulates N number of original particles.This means that coarse-graining is averaging the behavior of the original system, leading to a suppressed eruption zone, resulting in sharper solidity fronts.However, this averaging does not lead to changing bubble characteristics, which are found to be independent of the level of coarse-graining.This is very important as the bubble characteristics mainly determine the gas-solids contacting mechanisms.
In the second part of this work, we addressed the particle diameter to grid cell size ratio.By increasing the coarse-graining level, this particle diameter to grid cell size ratio becomes problematic for a conventional volume-weighted coupling.In order to maintain the flow field details, it is required to not enlarge the grid cell size.Therefore, we proposed a smooth and well-established fourthorder polynomial of Deen et al. (2004) as a solution.As shown, the particle diameter was well-exceeding the grid cell size for higher coarse-graining ratios while the scaling law of Sakai and Koshizuka (2009) was able to predict the original system for larger coarse-graining ratios.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.Fig. 13.Equivalent bubble diameter as a function of axial position using the method of Sakai and Koshizuka (2009) and the two-way coupling mapping function of Deen et al. (2004).

Fig. 1 .
Fig. 1.Schematic representation of the 3D fluidized bed column with dimensions 0.3 Â 0.06 Â 0.036 m (height Â width Â depth) including the initial particle lattice configuration.The boundary conditions for the gas flow were set as a uniform superficial air velocity equal of 0.50 m/s at the bottom, no-slip boundary conditions for the side walls and a fixed pressure boundary condition of 1 atm at the top of the domain.

Fig. 2 .
Fig.2.Instantaneous particle configurations taken at 2, 5, 7 and 8 s in the center of the column (0.5 depth) for the base case simulation (i.e.coarse-graining ratio equal to 1).

Fig. 3 .
Fig.3.Instantaneous particle configurations of coarse-graining ratios, l, equal to 2 (A, B, C) and 3 (D, E, F) at 2.5, 5 and 7.5 s in the center of the column (0.5 depth) for the scaling law developed byMu et al. (2020).

Fig. 4 .
Fig.4.Instantaneous particle configurations of coarse-graining ratios, l, equal to 2 (A, B, C) and 3 (D, E, F) at 2.5, 5 and 7.5 s in the center of the column (0.5 depth) for the scaling law developed bySakai and Koshizuka (2009).

Fig. 5 .Fig. 6 .
Fig.5.Mean and standard deviation dimensionless pressure drop using the column cross-sectional area, A, the applied gravitational acceleration, g, and the solids mass.

Fig. 7 .
Fig. 7. Average solids volume fraction at a height of 6 cm.Corresponding with earlier observations, the method of Mu et al. (2020) does not agree with the original system while the method ofSakai and Koshizuka (2009) shows good agreement except for large coarse-graining ratios.

Fig. 9 .
Fig.9.Dominant bubble frequency for the two scaling laws.The method ofMu et al. (2020) shows a clear diminishing bubble frequency up to 44% of the original base case value while the method ofSakai and Koshizuka (2009) results in a maximum bubble frequency deviation equal to 0.2 Hz compared to the base case value.

Fig. 10 .
Fig.10.Equivalent bubble diameter as a function of axial position.The method ofSakai and Koshizuka (2009) shows a better agreement with the original case compared to the scaling method ofMu et al. (2020).

Fig. 11 .
Fig. 11.Average solids volume fraction versus the axial position using the fourthorder polynomial.

Fig. 12 .
Fig. 12. Solids mass flux profiles at a height of 6 cm.Comparable to the previous section, the solids mass flux reduces in the center.

Table 1
Parameters used for the comparison simulations.The bold values indicate that proper scaling needs to be applied in either one of the two or both scaling laws to the shown base case value.The scaled simulation parameter values are presented in the supplementary material Sections 2 and 3.

Table 2
Sakai and Koshizuka (2009)ng the fourth-order polynomial for the two-way coupling.The bold values indicate that proper scaling needs to be applied in the method ofSakai and Koshizuka (2009).The scaled simulation parameter values are presented in the supplementary material.