Computational fluid dynamics (CFD)-discrete element method (DEM) simulations are designed to model a pseudo-two-dimensional (2D) fluidized bed, in which bed thickness is minimal compared to height and length. Predicted bed behavior varies as the simulations are conducted on increasingly refined computational grids. Pseudo-2D simulation results, in which a single computational cell spans the bed thickness, are compared against fully-three-dimensional (3D) simulations results. Both pseudo-2D and fully-3D simulations exhibit high accuracy when sufficiently refined. Indicators of bed behavior, such as bed height, bed height fluctuation, bubble generation frequency, and segregation, do not appear to converge as the cell size is reduced. The Koch-Hill and Gidaspow drag laws are alternately employed in the simulations, resulting in different trends of results with computational grid refinement. Grid refinement studies are used to quantify the change in results with grid refinement for both three-dimensional, uniform refinement, and for two-dimensional refinement on pseudo-2D computational grids. Grid refinement study results indicate the total drag converges as the computational grid is refined, for both 3D and pseudo-2D approaches. The grid refinement study results are also used to distinguish the relatively grid-independent results using the Koch-Hill drag law from the highly grid-dependent Gidaspow drag law results. Computational cell size has a significant impact on CFD-DEM results for fluidized beds, but the grid refinement study method can be used to quantify the resulting numerical error.

## Introduction

Fluidized bed research has become more widespread due to increased usage in industries such as pharmaceutical manufacturing and food processing, as well as newly developed applications like energy production from biomass [1]. Fluidized beds are excellent for coating, granulation, and drying applications due to high mixing and heat transfer rates [2,3]. Research in this area includes many studies for pseudo-two-dimensional (2D) fluidized beds, in which one dimension of the experimental bed is reduced in size sufficiently to support the assumption that gradients in the reduced dimension are negligible compared to gradients in the remaining two dimensions. Pseudo-2D beds are commonly employed due to ease of experimental data collection [4,5], and the low computational cost of modeling [2,3,6]. Pseudo-2D fluidized bed experiments have brought about a large increase in knowledge [7], but comparison of computational modeling results for pseudo-2D beds with experimental results has revealed large discrepancies [4,8].

There has been increasing discussion of if, and when, 2D simulations can be used to predict phenomena that occur in three-dimensional (3D) experiments. Xie et al. [3,9] found that the 3D experiments can be represented by 2D simulations when the forces are predominantly in only two directions. Pseudo-2D fluidized beds fit this description with gradients predominantly in only two directions, and therefore it is possible that a 2D model will be appropriate [9].

A significant factor in simulation accuracy is the computational cell size. Peng et al. [10] were the first to quantify the upper and lower critical cell bounds for cell size. Similarly, while simulating a bubble column using a 2D Eulerian–Eulerian mathematical model, Picardi et al. [6] reported that simulation accuracy did not necessarily increase with computational grid refinement. We have reported reductions in simulation accuracy for both over- and under-refined computational grids for computational fluid dynamics-discrete element method (CFD-DEM) simulation of flow through a fixed-particle bed [11].

We have investigated how simulation accuracy varies within this critical range, as well as distinct simulation solution behaviors that occur if the cell size is outside of the critical bounds. We found cell size had a significant effect on simulation accuracy for flow through fixed-particle beds, but this effect could be quantified using a grid refinement study under certain conditions [11,12]. To the author's knowledge, changes in predicted bed behavior, as the computational cell size is systematically modified, have not been reported for fluidized bed simulations.

In the present study, we vary computational grid refinement in all three dimensions to examine how grid refinement and the pseudo-2D assumption affect accuracy of the simulation results [13] for a fluidized bed. We compare our CFD-DEM simulation results with experimental data published by Goldschmidt et al. [14], for a fluidized bed of binary-sized particles. Bed-height fluctuation frequency, segregation of the fine and coarse particle components, and height fluctuation of the two components are compared with experimental data. The average number and size of bubbles are examined for each computational grid, as the extent of segregation is known to be highly dependent on bubbling characteristics [15]. Finally, we quantify expected numerical error in the simulations using a grid refinement study.

We begin with an overview of the CFD-DEM mathematical model, including our specific choices for void fraction and fluid–particle interaction calculations. Next, the Computational Model section details the physical system being modeled and the construction of computational grids. Results are analyzed with emphasis on their dependence on the computational grid size and drag law. We conclude this paper with our contributions and several suggestions for future work.

## Mathematical Model

The CFD-DEM mathematical model consists of CFD equations for the fluid phase, DEM equations for dynamics of the solid particles, and a coupling routine. The fluid-phase equations include the conservation of mass and conservation of momentum, written as
$∂ε∂t+▽·(εuf)=0$
(1)
$∂(εuf)∂t+▽·(εufuf)=−ε▽pρf−Rpf+▽·τ$
(2)

Turbulence is included in the mathematical model with the standard kε turbulence model for incompressible flow [2]. We compare results with this turbulence model to simulation results using the realizable kε turbulence model, and to results using the laminar assumption.

The DEM calculations are performed for each particle individually. The linear-momentum equation
$midvidt=∑jFijc+fgrav,i+ff,i$
(3)
calculates the particle acceleration from the sum of contact forces ($Fijc$), gravitational forces (fgrav,i), and fluid–solid interaction force (ff,i). The angular-momentum equation
$Iidwidt=∑j(Ri×Ft,ij)$
(4)
is based on the sum of tangential contact forces (Ft,ij), and the particle's moment of inertia (Ii). The normal and tangential contact forces arising from each collision are calculated using the Hertz contact model [16], which is employed to model both particle–particle and particle–wall collisions.
We employ an unresolved coupling routine, which treats the fluid–solid mixture in each cell as a continuum [17]. The coupling routine first calculates the void fraction in each of the CFD computational cells using the divided volume fraction method [18]. The particle–fluid interaction force,
$ff,i=Viβ(1−ε)[uf−vi]$
(5)
typically called the drag force, is calculated from the volume of the particle (Vi), the void fraction in each cell (ε), and the difference between the particle velocity (vi) and surrounding fluid velocity(uf). The coefficient of interphase momentum transfer, β in Eq. (5), is based on the chosen drag law.
The Koch-Hill drag law was developed from lattice-Boltzmann simulations of fluidized particles [19,20], and is therefore appropriate for this application [21]. Unless otherwise noted, the models presented use the Koch-Hill drag law, where β is given as
$β=18μfε2(1−ε)dp2[F0(ε)+12F3(ε)Rep]$
(6)

where F0 and F3 are given by van Buijtenen et al. [8].

Simulation results using the Gidaspow drag law (when specifically noted) are also included to examine how results utilizing different drag laws behave differently with grid refinement. The Gidaspow drag law [22], where β is given as
$β={150(1−ε)2μfεdp2+1.75(1−ε)ρfdp|uf−vi|for (ε<0.8)34Cd(1−ε)ρfdp|uf−vi|ε−2.65for (ε≥0.8)$
(7)
is a combination of the Ergun drag law (for ε < 0.8) and the Wen and Yu drag law (for ε ≥ 0.8).

Simulations are performed using the open-source CFDEM project [18], a coupling between OpenFOAM and LIGGGHTS.

## Computational Model

The computational model is designed to represent the pseudo-2D fluidized bed experiment of Goldschmidt et al. [14]. First, a description of the experiment is given, with a listing of all input quantities. An overview of the computational grids, simulation cases, and postprocessing is also provided.

### Experimental Description.

For their detailed experimental publication, Goldschmidt et al. [14] fluidized a binary-sized mixture of particles in a pseudo-2D bed, and reported quantitative data for bed expansion, bed oscillation frequency, and segregation [14]. Input properties for our simulations, including bed dimensions, particle diameters, number of particles, and particle properties, are based on their reported data, as shown in Table 1.

Table 1

Input properties (Goldschmidt et al. [14] and Beetstra et al. [23])

Particle properties
Lg. Particle diameter (dlarge)2.49 mm
Sm. particle diameter (dsmall)1.52 mm
Number Sm. particles27,720
Number Lg. particles17,940
Density (ρs)2523 kgm–3
Coeff. of restitution0.97
Coeff. of friction0.1
Fraction of fines (xf)25%
Sauter mean diameter (dS)2.14 mm
Fluid properties
Density (ρf)1.2041 kgm–3
Viscosity (μf)1.98 × 10–5 kgm–1 s–1
Inlet velocity1.3 ms−1
Domain size
Thickness (T)15 mm
Length (L)150 mm
Height (H)450 mm
Particle properties
Lg. Particle diameter (dlarge)2.49 mm
Sm. particle diameter (dsmall)1.52 mm
Number Sm. particles27,720
Number Lg. particles17,940
Density (ρs)2523 kgm–3
Coeff. of restitution0.97
Coeff. of friction0.1
Fraction of fines (xf)25%
Sauter mean diameter (dS)2.14 mm
Fluid properties
Density (ρf)1.2041 kgm–3
Viscosity (μf)1.98 × 10–5 kgm–1 s–1
Inlet velocity1.3 ms−1
Domain size
Thickness (T)15 mm
Length (L)150 mm
Height (H)450 mm

Fluid properties were not reported by Goldschmidt et al. [14]. To remain consistent with previously published models of this experiment, we used the fluid properties from Ref. [23], as shown in Table 1.

### Computational Grids.

Simulations are performed on multiple computational grids, to investigate the effect of computational cell size, and the pseudo-2D assumption on the results. For each computational grid, a unique descriptor, number of cells in each direction, and cell size in each direction are listed in Table 2.

Table 2

Table of cases

Number of cellsCell size (mm)
(Direction)(Direction)
Grid(T)(L)(H)(T)(L)(H)R
T1-L101103015.015.015.07.0
T1-L201206015.07.57.54.4
T1-L301309015.05.05.03.4
T1-L4014012015.03.83.82.8
T2-L10210307.515.015.05.6
T2-L20220607.57.57.53.5
T2-L30230907.55.05.02.7
T2-L402401207.53.83.82.2
T3-L10310305.015.015.04.9
T3-L20320605.07.57.53.1
T3-L30330905.05.05.02.3
T3-L403401205.03.83.81.9
T4-L20420603.87.57.52.8
T4-L404401203.83.83.81.8
Number of cellsCell size (mm)
(Direction)(Direction)
Grid(T)(L)(H)(T)(L)(H)R
T1-L101103015.015.015.07.0
T1-L201206015.07.57.54.4
T1-L301309015.05.05.03.4
T1-L4014012015.03.83.82.8
T2-L10210307.515.015.05.6
T2-L20220607.57.57.53.5
T2-L30230907.55.05.02.7
T2-L402401207.53.83.82.2
T3-L10310305.015.015.04.9
T3-L20320605.07.57.53.1
T3-L30330905.05.05.02.3
T3-L403401205.03.83.81.9
T4-L20420603.87.57.52.8
T4-L404401203.83.83.81.8
The level of refinement of each computational grid is quantified by a nondimensional cell size, R, which describes the grid size relative to the average size of the particle mixture [10]. Calculation of the nondimensional cell size first requires calculation of the volume-average cell size and the Sauter-mean particle diameter. The volume-average cell size for each grid can be calculated as
$h=[1N∑i=1N△∀i]1/3$
(8)

where N is the total number of computational cells, and $∀i$ is the volume of each individual cell [24]. For grids with cubic cells, the volume-average cell size is equivalent to the side length of a computational cell.

The Sauter-mean diameter is a good description of the average particle size in the mixture, and is calculated as
$dS=1∑jxjdj$
(9)

where xj is the mass fraction, and dj is the particle diameter of each component of the binary mixture. The Sauter-mean diameter of the 25% fines particle mixture is 2.14 mm, as listed in Table 1.

The ratio of volume-averaged cell size (h) to Sauter-mean particle diameter (dS)
$R=hdS$
(10)
serves as the nondimensional cell size (R), and is recorded for each computational grid in Table 2.

The simulation is performed on a total of 14 computational grids, listed in Table 2. The naming convention is based on the number of cells in each direction. For example, computational grid T2-L30 has two cells in the thickness direction and 30 cells in the length direction. Cells' sizes in the length and height direction are always equivalent for this study, as gradients in both of these directions are expected to be large due to the presence of bubbles.

In Table 2, the computational grids are grouped by the number of cells in the thickness direction, with the first group containing all pseudo-2D computational grids. The grids in each group represent a 2D refinement with a constant number of cells in the thickness direction. Grids in bold font have cubic (equilateral) cells, and the four bold grids represent a uniform, 3D refinement. The 3D refinement is performed with the smallest possible refinement intervals (based on the change in cell length in the thickness direction). These refinement intervals are used for all other refinement types for consistent comparison. A combination of all grids with L10 forms a one-dimensional (1D) refinement, and the same can be said for the grids with L20, L30, and L40, individually. Therefore, varying combinations of the 14 grids result in one 3D refinement, three 2D refinements, and four 1D refinements. The eight refinement types and their associated grids, listed from least to most refined, are listed in Table 3. Note that each refinement group is labeled by the type of refinement (1D, 2D or 3D) followed by a description of the level of refinement in the directions not being actively refined (coarse, med, fine, finest).

Table 3

Table of refinement combinations

LabelLeast—→Most refined
3DT1-L10T2-L20T3-L30T4-L40
2D_coarseT1-L10T1-L20T1-L30T1-L40
2D_medT2-L10T2-L20T2-L30T2-L40
2D_fineT3-L10T3-L20T3-L30T3-L40
1D_coarseT1-L10T2-L10T3-L10
1D_medT1-L20T2-L20T3-L20T4-L20
1D_fineT1-L30T2-L30T3-L30
1D_finestT1-L40T2-L40T3-L40T4-L40
LabelLeast—→Most refined
3DT1-L10T2-L20T3-L30T4-L40
2D_coarseT1-L10T1-L20T1-L30T1-L40
2D_medT2-L10T2-L20T2-L30T2-L40
2D_fineT3-L10T3-L20T3-L30T3-L40
1D_coarseT1-L10T2-L10T3-L10
1D_medT1-L20T2-L20T3-L20T4-L20
1D_fineT1-L30T2-L30T3-L30
1D_finestT1-L40T2-L40T3-L40T4-L40

Additionally, we have conducted 1D refinement in the height direction, and details of computational grid configuration and corresponding simulation results are included in the Appendix. These simulation results are used to provide additional detail for determining the placement of critical cell sizes. The inclusion of these simulation results is noted in Secs. 4.3 and 4.4.

The boundary conditions are consistent for simulations across all of the computational grids, with CFD boundary conditions shown on a sketch of the domain in Fig. 1. In the CFD calculation, the inlet is maintained at a uniform velocity, while the outlet is assigned as constant atmospheric pressure. The walls are treated as slip boundaries, and the effects of this choice are analyzed later. In the DEM calculation, the inlet and sides are assigned as walls to prevent simulated particles from exiting the domain.

Fig. 1
Fig. 1
Close modal

### Postprocessing.

Results from each simulation are postprocessed to compare with the results of Goldschmidt et al. [14]. First, the vertical height of the center of mass of all of the bed material is captured at a frequency of 25 Hz. The average center of mass height over the entire simulated time is calculated as
$⟨hbed⟩=∑tmintmaxhbed(t)fΔt$
(11)

where f is the frequency of data collection in Hz, and Δt is the total time for which data are collected.

The variation of the bed center of mass height from the average is reported as a root-mean-square (RMS) of the difference at each time, calculated as
$RMS(hbed)=∑tmintmax(hbed(t)−⟨hbed⟩)2fΔt$
(12)

The bed RMS indicates the relative amount of movement of the particle mixture. Additionally, a fast Fourier transform is applied to the bed center of mass height with time, to determine the dominant oscillation frequencies.

Next, the center of mass height of the small and large particles (hsmall, hlarge) is captured at a frequency of 25 Hz, and the center of mass height is averaged over one second intervals. The data are analyzed using a regression to obtain a best-fit second-order polynomial of the form
$hsmallfit(t)=at2+bt+c$
(13)
shown for the small particles, where a, b, and c are regression coefficients.
For each particle type, the RMS variation in center of mass height against the best-fit polynomial is calculated as
$RMS(hsmall)=∑tmintmax(hsmall(t)−hsmallfit(t))2fΔt$
(14)
again shown for the small particles.
The degree of segregation over the first ten seconds, for a mixture with 25% weight fraction fines, is calculated from the large and small particle average center of mass heights [14] as
$s=34[⟨hsmall⟩⟨hlarge⟩−1]$
(15)

where the individual particle components' average center of mass height, $⟨hlarge⟩$ and $⟨hsmall⟩$, are calculated using Eq. (11). The coefficient of 3/4 is based on Goldschmidt et al.'s [14] segregation equation and the 25% weight fraction of small particles in the simulated mixture. The segregation degree in Eq. (15) is the average degree of segregation over the first ten seconds, while the instantaneous degree of segregation at the end of the simulation is expected to be greater.

The bed behavior is defined by the intensity of bubbling. The simulations are performed for conditions within the bubbling fluidization regime, so the average number and size of the bubbles are indicative of the bubbling intensity.

For uniformity of postprocessing across all of the computational grids, the bed is treated in a pseudo-2D manner. This is consistent with the digital image analysis used by Goldschmidt et al. [14]. Viewed from the front of the bed, the computational domain is divided into a 2D postprocessing grid consisting of the same number of cells in the length and height direction as the computational grid. For each 2D postprocessing cell within the dense-phase particle bed, if any computational grid cell in the direction of bed thickness has a void fraction greater than 0.8, a typical cutoff for bubble detection [7], the postprocessing cell is identified as residing within a bubble. Once this analysis is completed on all computational cells, the postprocessing cells residing in bubbles are grouped to identify bubbles spanning multiple cells. The number, diameter, and overall average diameter of bubbles are reported.

#### Grid Refinement Study.

We perform a grid refinement study to quantify how the solution changes as the cell size is reduced. The calculation method is described by Celik et al. [24], along with a discussion of the significance of the calculation quantities. We confine our discussion of the grid refinement results to the observed order of accuracy (p) and the grid convergence index (GCI). The observed order of accuracy is the apparent rate at which the solution converges to a single value as the cell size is reduced. The GCI is the expected amount of numerical error in the solution calculated on the finest grid, based on the observed solution trend. A very small value of GCI is desired, as this indicates minimal error in the solution originating from the size of the computational cells.

## Results

Results are reported with the eight refinement methods based on simulation results on 14 total computational grids. We begin by examining the effect of boundary conditions and turbulence model on the accuracy of our model. After using these results to justify our model choices, we analyze how results are affected by the level, and type, of grid refinement. The bed expansion, bed oscillation frequency, segregation index, and bubbling intensity are examined for all eight refinement methods, and major trends are discussed. Results from simulations utilizing the Gidaspow drag law are compared with simulations using the Koch-Hill drag law. Finally, a grid refinement study is used to quantify the results' changes with grid refinement, and to further justify our choice of drag law. The results are discussed in terms of an optimal cell size range in which the results remain within an acceptable level of error.

### Effect of Boundary Conditions and Turbulence Model.

We begin by analyzing our choices of boundary conditions and turbulence model. Regarding boundary conditions, Beetstra et al. [23] noted that neither slip nor nonslip boundary conditions accurately represent the fluid–wall interaction for simulations of pseduo-2D fluidized beds. Substituting nonslip boundary conditions as an alternate to the chosen slip condition on the short sidewalls, we find minimal changes to simulation results, except for simulations on pseudo-2D grids. When the sidewalls are designated nonslip and the third dimension is represented by only a single cell, the bubbling is almost completely suppressed. The same suppression of bubbling occurs when all walls are designated nonslip (sidewalls and front/rear walls), regardless of the number of cells in the thickness direction. We choose slip boundary conditions on all walls, since this best represents physical reality across all of the computational grids.

We use the standard k-epsilon model [2] to model turbulence, but we also analyze use of the realizable k-epsilon model. Simulation results utilizing the two models vary minimally, with a 0.4% change in average drag force and a 0.05% change in average bed height. The results from the realizable k-epsilon model predict greater segregation with a 1.7% higher segregation index.

We also compare our simulation results with the standard k-epsilon model to simulations results employing the laminar assumption. Simulation results with the laminar assumption also varied negligibly from those of the standard k-epsilon model, with a 0.6% change in average drag force, and a 1.2% change in bed height. The laminar assumption resulted in a segregation index 2.7% higher than the simulation results utilizing the standard k-epsilon model. We choose the standard k-epsilon model for computational efficiency, as these results indicate the choice of turbulence model does not significantly affect model accuracy.

### Bed Expansion and Oscillation.

The average height of the particle bed center of mass, $⟨hbed⟩$, and the amount of variation of this height, RMSbed, are plotted against the computational cell size in Figs. 2 and 3, respectively. Both the bed height and RMS generally increase with decreasing cell size for 2D and 3D refinement methods. The 3D (dashed) line represents refinement in all three directions, starting with a coarse, pseudo-2D grid and terminating in a highly refined, fully 3D grid. The 2D (solid) lines represent refinement in the width and height direction, while maintaining a constant grid size in the thickness (pseudo-2D) direction. Solution result trends for bed height and RMS are similar for 2D and 3D refinement in Figs. 2 and 3. Refinement only in the direction of thickness (1D, dotted lines) has minimal effect on the output bed height and RMS, regardless of level of refinement in the other two directions.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

Increases in bed height and RMS with increased grid refinement are due to increases in bubble intensity, which is described in detail in Sec. 4.4. Increase in the area of voids in the bed increases the bed height, and the amount of bed height fluctuation (RMS) as the voids move through the bed.

A Fourier transform is performed on the bed center of mass height with time. All simulations over-predict the dominant frequency, reporting values between 2.1 and 2.7 Hz, compared to Goldchmidt et al.'s reported 1.9 Hz [14]. Secondary and tertiary dominant frequencies are recorded mainly between 2.0 and 2.8 Hz, with several reporting a secondary frequency of 0.1 Hz. The predicted frequency results show no obvious trend with grid refinement or grid refinement type.

### Segregation.

Goldschmidt et al. [14] report a segregation degree between 14.1% and 17.4%. Their reported segregation degree from three experiments varies over 3%, so segregation degree slightly outside this range is not unexpected. Simulation results from all grids are within, or very close to this range of expected segregation, as shown in Fig. 4. There is a generally increasing trend of segregation with decreasing cell size, but the results fluctuate considerably.

Fig. 4
Fig. 4
Close modal

Additionally, Goldschmidt et al. [14] report average height and RMS of the large and small particles. Simulations consistently underpredict the small and large average particle height, possibly due to a difference in the number of particles (Goldschmidt et al. [14] only give approximate numbers). We normalize the RMS of each particle species by the species average bed height. Results for the large and small particle normalized RMS are shown in Figs. 5 and 6, respectively. The upper and lower normalized RMS from Goldschmidt et al. [14] are included in both figures. The simulation results remain within, or very close to, the expected normalized RMS range, except for the finest grids, namely those included in the 1D_Finest refinement. All computational grids in the 1D_Finest refinement have cell dimensions of 1.75dp,sauter in the length and height directions. Simulations performed on these finest grids overestimate the amount of movement of the small particles. Additional simulations (see the Appendix) included in Fig. 6 are consistent with this observation, with overpredictions of the small particle movement when the cells dimensions are less than two Sauter-mean diameters in any direction with non-negligible gradients (height and length for this case). This reduction in accuracy on highly refined grids is consistent with previous publications of a lower critical grid size (minimum acceptable cell size) [6,10,11]. The overestimated movement of smaller particles, along with the increasing movement of the large particles, on the finest grids, may be due to the previously reported increase in standard deviation of porosity as the cell size is reduced [12]. Large variation in the porosity profile will cause larger variation in the calculated drag force, overestimating the local fluid drag force on the particles for the finest grids. For future simulations, we recommend cell lengths in all direction to be two Sauter-mean particle diameters or longer.

Fig. 5
Fig. 5
Close modal
Fig. 6
Fig. 6
Close modal

### Bubbling Intensity.

The bubbling intensity is determined by the average number of bubbles present in the bed at any time and the average diameter of the bubbles. The average number of bubbles and average bubble diameter are plotted against the computational cell size in Figs. 7 and 8, respectively. Additional simulation results based on 1D refinement in the thickness direction (see the Appendix) are included in Fig. 7. It is clear in Fig. 7 that the number of bubbles is highly dependent on the level of grid refinement in the flow direction (2D and 3D refinements). The predicted number of bubbles increases as the cell size is decreased, which is expected, as calculations on finer grids can detect small bubbles missed by coarser grids. Alternately, the predicted number of bubbles appears to be almost independent of the level of refinement in the thickness direction (1D refinement). We know from Goldschmidt et al. [14] that the physical conditions defining our computational model correspond to an actively bubbling fluidized bed. We expect an average of at least one bubble in the domain, at minimum. When simulations predict lower than this limit, we can assume that the computational grid is too coarse. At a nondimensional computational cell size of 4.4, one of these additional simulations predicts an average number of bubbles less than one, while the other two simulations conducted at this level of refinement predict an average number of bubbles greater than one. All simulations with nondimensional cell size less than 4.4 predict an average number of bubbles greater than one. Therefore, the critical cell size is around 4.4 Sauter-mean particle diameters, and we recommend the use of nondimensional cell sizes less than this value.

Fig. 7
Fig. 7
Close modal
Fig. 8
Fig. 8
Close modal

The average bubble diameter generally increases as the computational cell size is decreased in Fig. 8. The one-directional refinements with coarse and medium grids report bubble diameters that appear to be independent of cell size. However, the bubble diameters from one-directional refinements with fine and finest grids are highly dependent on cell size. Refinement in the thickness direction on these finer grids allows for capturing of bubbles that only span a fraction of the full bed thickness, resulting in a rapid increase in predicted bubble diameter with refinement. The 2D and 3D refinements exhibit a consistent increase in predicted bubble diameter with decreasing cell size, except for the 2D_Coarse (pseudo-2D) refinement. The pseudo-2D grids may be unable to resolve the full outline of each bubble as the bubbles do not always span the full bed thickness.

Overall, the predicted amount of bubbling intensity increases with grid refinement, especially with refinement in the flow direction. This trend is due to better resolution with grid refinement, and may also be due to nonlinearity of the drag function with void fraction.

### Effect of Drag Law.

Additional simulations are performed on the computational grids in the 3D refinement using the Gidaspow drag law. Simulations using the two drag laws report a similar number of bubbles in the bed. Simulations using the Gidaspow drag law report both a higher void fraction in the expanded bed, and larger average bubble diameters, resulting in a significantly higher predicted average bed height.

The main difference in results with the Gidaspow drag law is the segregation index, shown in Fig. 9. Simulations using the Koch-Hill drag law report segregation index values close to the experimentally reported values as the computational grid are refined. The simulation on the coarsest grid using the Gidaspow drag law reports a segregation index slightly above the reported experimental range. In contrast, the Gidaspow simulations on refined grids all report segregation index values far below the experimental range. Our choice of the Koch-Hill drag law is justified by the ability to capture the known segregation phenomenon at all levels of grid refinement.

Fig. 9
Fig. 9
Close modal

### Grid Refinement Study.

We perform a grid refinement study using the drag force value across uniformly refined grids. We have two sets of uniformly refined grids: 3D refinement and the set of pseudo-2D grids (2D_Coarse). The drag force profile with computational cell size is shown in Fig. 10 for the two refinement sets. Each set consists of four grids, which can be combined in four different ways for processing by the grid refinement study procedure.

Fig. 10
Fig. 10
Close modal

Grid refinement studies for the 3D refined grids report an observed order of accuracy between 2.5 and 3.5 and GCI between 0.05% and 0.5%. For the uniformly refined pseudo-2D grids, the order of accuracy is between 2 and 4.5, with a slightly increased GCI between 0.15% and 1.6%. The grid refinement study for both of these refinement types indicates quick solution convergence and minimal numerical error. An additional grid refinement study is performed for the 3D refined grids with simulations utilizing the Gidaspow drag law in place of the Koch-Hill drag law. The grid refinement study using the Gidaspow simulations reports the observed order of accuracy between 0.15 and 0.85. The GCI values range from 3% to 26%, significantly higher than the GCI calculated for simulations utilizing the Koch-Hill drag law. These results agree with previous statements [15] that the Koch-Hill drag law is better suited for this application, as it was derived from LB simulations of similar conditions.

The CFD numerical solution is expected to converge with a spatial order of accuracy of two. However, the definition of a formal order of accuracy for the CFD-DEM computational scheme is hampered by the interaction between fluid and solids. The grid refinement studies indicate that despite this influence, the result for total drag force converges slightly faster than expected when the mathematical model is appropriate (Koch-Hill). Other model behavior parameters do not show the same convergence behavior with grid refinement as the total drag force, such as the rapid increase in the number of bubbles with grid refinement (Fig. 7). The apparent divergence of some parameters as the grid is refined stems from the complex interaction of numerical convergence with increased spatial resolution of the fluid–particle interaction.

## Conclusion

We perform simulations of a pseudo-2D fluidized bed on 14 total computational grids with different levels and types of refinement, including four pseudo-2D grids. The 14 computational grid solutions are combined in eight different ways to analyze the effect of 1D, 2D, and 3D refinement. Resultant variables such as bed height, bed oscillation, segregation index, small and large particle oscillation, number of bubbles, and average bubble diameter are plotted, with refinement type, against computational cell size.

Reduction of the grid to pseudo-2D (one computational cell spans the bed thickness) does not greatly affect results. However, the refinement in the remaining two dimensions significantly affects the results. One-dimensional refinement in the thickness direction appears to have a negligible effect on results, compared to 2D and 3D refinement. This is consistent with grid generation theory for CFD; refinement should be conducted in the directions of the largest gradients to enhance resolution of the flow. Based on our overall findings, we recommend a volume-average computational cell size less 4.4 particle diameters, and a minimum cell length of two particle diameters in any direction with non-negligible gradients. In the future, simulation of larger 3D fluidized beds should be conducted with varying computational cell size, in order to more precisely determine the limits of the critical cell size range.

Grid refinement studies converge for both fully 3D refined grids and for the set of pseudo-2D grids. The grid refinement study is able to capture the convergence behavior of the total drag force with grid refinement, and divergence behavior when the mathematical model was not appropriate due to use of the Gidaspow drag law.

For future simulations of pseudo-2D fluidized beds, pseudo-2D simulations are shown to be able to accurately capture the system behavior, given sufficient refinement in the two remaining dimensions. Studies of a similar nature should be conducted to assess suitability of pseudo-2D computational grids to similar systems. Theoretical analysis targeting the effect of computational cell size on the CFD-DEM coupling routine should be performed in order to explain the large variation of predicted bed behavior with cell size in this study.

To the authors' knowledge, this is the first application of the grid refinement procedure to CFD-DEM simulations of fluidized beds. For future studies, the use of a grid refinement study is suggested to determine numerical error due to the relatively large cell sizes needed for CFD-DEM computations.

## Acknowledgment

Additional funding for this work was provided by Procter and Gamble and the UC Simulation Center.

## Funding Data

• National Science Foundation Graduate Research Fellowship Grant No. (1102690).

## Nomenclature

• a =

•
• b =

linear regression coefficient

•
• c =

constant regression coefficient

•
• Cd =

coefficient of drag

•
• dp =

particle diameter

•
• dS =

Sauter-mean particle diameter

•
• $dsmall,dlarge$ =

small and large particle diameters, respectively

•
• f =

frequency of data collection

•
• ff,i =

fluid force on particle i

•
• fgrav,i =

gravitational force on particle i

•
• Ft,ij =

tangential contact force exerted on particle i by particle j

•
• $Fijc$ =

contact force exerted on particle i by particle j

•
• F0, F3 =

fitted functions for Koch-Hill drag law

•
• h =

volume-average cell size

•
• hfit(t) =

fitting equation for center of mass of each particle type

•
• $⟨hlarge⟩$ =

vertical center of mass of large particles

•
• $⟨hsmall⟩$ =

vertical center of mass of small particles

•
• Ii =

moment of inertia of particle i

•
• mi =

mass of particle i

•
• N =

total number of computational cells

•
• p =

pressure

•
• R =

nondimensional cell size

•
• Ri =

•
• Rep =

particle Reynolds number, $Rep=ερfdp|uf−vi|μf$

•
• Rpf =

momentum exchange between fluid and particle

•
• s =

amount of segregation

•
• t =

time

•
• tmax =

time instant of final data collection

•
• tmin =

time instant of initial data collection

•
• uf =

fluid velocity

•
• vi =

velocity of particle i

•
• Vi =

volume of particle i

•
• wi =

rotational velocity of particle i

•
• xj =

weight fraction of particle type j

•
• $∀i$ =

volume of computational cell i

•
• H =

computational domain height

•
• L =

computational domain length

•
• T =

computational domain thickness

### Greek Symbols

Greek Symbols

• β =

coefficient of inter-phase momentum transfer

•
• ε =

void fraction

•
• μf =

fluid absolute viscosity

•
• ρf =

fluid density

•
• τ =

fluid phase stress tensor

### Appendix

Several simulations were conducted in order to assess the effect of computational grid refinement in the height direction on the solution. These simulation results are labeled “Addtl. Sims” when plotted. Computational grid setup and corresponding results for average number of bubbles and small particle normalized RMS are included in Table 4.

Table 4

Table of cases and results for investigation of computational grid refinement in height direction

Number of cells (direction)
Grid(T)(L)(H)RNumber of bubblesRMSsmall/$⟨hsmall⟩$
T1-L10-H601100603.40.99.1%
T1-L10-H90110902.71.49.4%
T1-L10-H1201101202.11.810.3%
T2-L20-H30220305.60.78.6%
T2-L20-H90220904.92.08.8%
T2-L20-H1202201204.43.010.5%
T3-L30-H30330304.41.28.7%
T3-L30-H60330603.11.89.0%
T3-L30-H1203301202.85.011.3%
Number of cells (direction)
Grid(T)(L)(H)RNumber of bubblesRMSsmall/$⟨hsmall⟩$
T1-L10-H601100603.40.99.1%
T1-L10-H90110902.71.49.4%
T1-L10-H1201101202.11.810.3%
T2-L20-H30220305.60.78.6%
T2-L20-H90220904.92.08.8%
T2-L20-H1202201204.43.010.5%
T3-L30-H30330304.41.28.7%
T3-L30-H60330603.11.89.0%
T3-L30-H1203301202.85.011.3%

## References

1.
Teaters
,
L. C.
, and
Battaglia
,
F.
,
2014
, “
On the Computational Modeling of Unfluidized and Fluidized Bed Dynamics
,”
ASME J. Fluids Eng.
,
136
(
10
), p.
104501
.
2.
Luo
,
K.
,
Wu
,
F.
,
Yang
,
S.
, and
Fan
,
J.
,
2015
, “
CFD-DEM Study of Mixing and Dispersion Behaviors of Solid Phase in a Bubbling Fluidized Bed
,”
Powder Technol.
,
274
, pp.
482
493
.
3.
Xie
,
N.
,
Battaglia
,
F.
, and
Pannala
,
S.
,
2008
, “
Effects of Using Two- Versus Three-Dimensional Computational Modeling of Fluidized Beds: Part I, Hydrodynamics
,”
Powder Technol.
,
182
(
1
), pp.
1
13
.
4.
Busciglio
,
A.
,
Vella
,
G.
,
Micale
,
G.
, and
Rizzuti
,
L.
,
2009
, “
Analysis of the Bubbling Behaviour of 2D Gas Solid Fluidized Beds—Part II: Comparison Between Experiments and Numerical Simulations Via Digital Image Analysis Technique
,”
Chem. Eng. J.
,
148
(
1
), pp.
145
163
.
5.
Taghipour
,
F.
,
Ellis
,
N.
, and
Wong
,
C.
,
2005
, “
Experimental and Computational Study of Gas-Solid Fluidized Bed Hydrodynamics
,”
Chem. Eng. Sci.
,
60
(
24
), pp.
6857
6867
.
6.
Picardi
,
R.
,
Zhao
,
L.
, and
Battaglia
,
F.
,
2016
, “
On the Ideal Grid Resolution for Two-Dimensional Eulerian Modeling of Gas-Liquid Flows
,”
ASME J. Fluids Eng.
,
138
(
11
), p.
114503
.
7.
Li
,
T.
,
Gopalakrishnan
,
P.
,
Garg
,
R.
, and
Shahnam
,
M.
,
2012
, “
CFD-DEM Study of Effect of Bed Thickness for Bubbling Fluidized Beds
,”
Particuology
,
10
(
5
), pp.
532
541
.
8.
Van Buijtenen
,
M. S.
,
Van Dijk
,
W. J.
,
Deen
,
N. G.
,
Kuipers
,
J. A. M.
,
,
T.
, and
Parker
,
D.
,
2011
, “
Numerical and Experimental Study on Multiple-Spout Fluidized Beds
,”
Chem. Eng. Sci.
,
66
(
11
), pp.
2368
2376
.
9.
Xie
,
N.
,
Battaglia
,
F.
, and
Pannala
,
S.
,
2008
, “
Effects of Using Two- Versus Three-Dimensional Computational Modeling of Fluidized Beds: Part II, Budget Analysis
,”
Powder Technol.
,
182
(
1
), pp.
14
24
.
10.
Peng
,
Z.
,
Doroodchi
,
E.
,
Luo
,
C.
, and
,
B.
,
2014
, “
Influence of Void Fraction Calculation on Fidelity of CFD-DEM Simulation of Gas-Solid Bubbling Fluidized Beds
,”
AIChE J.
,
60
(
6
), pp.
2000
2018
.
11.
Volk
,
A.
,
Ghia
,
U.
, and
Stoltz
,
C.
,
2017
, “
Effect of Grid Type and Refinement Method on CFD-DEM Solution Trend With Grid Size
,”
Powder Technol.
,
311
, pp.
137
146
.
12.
Volk
,
A.
,
Ghia
,
U.
, and
Liu
,
G. R.
,
2018
, “
Assessment of CFD-DEM Solution Error against Computational Cell Size for Flows Through a Fixed-Bed of Binary-Sized Particles
,”
Powder Technol.
,
325
, pp. 519–529.
13.
Volk
,
A.
, and
Ghia
,
U.
,
2017
, “
How Computational Grid Refinement in Three Dimensions Affects CFD-DEM Results for Psuedo-2D Fluidized Gas-Solid Beds
,”
ASME
Paper No.
FEDSM2017-69222.
14.
Goldschmidt
,
M. J. V.
,
,
J. M.
,
Mellema
,
S.
, and
Kuipers
,
J. A. M.
,
2003
, “
Digital Image Analysis Measurements of Bed Expansion and Segregation Dynamics in Dense Gas-Fluidised Beds
,”
Powder Technol.
,
138
(
2–3
), pp.
135
159
.
15.
Bokkers
,
G. A.
,
van Sint Annaland
,
M.
, and
Kuipers
,
J. A. M.
,
2004
, “
Mixing and Segregation in a Bidisperse Gas-Solid Fluidised Bed: A Numerical and Experimental Study
,”
Powder Technol.
,
140
(
3
), pp.
176
186
.
16.
Zhu
,
H. P.
,
Zhou
,
Z. Y.
,
Yang
,
R. Y.
, and
Yu
,
A. B.
,
2007
, “
Discrete Particle Simulation of Particulate Systems: Theoretical Developments
,”
Chem. Eng. Sci.
,
62
(
13
), pp.
3378
3396
.
17.
Hager
,
A.
,
Kloss
,
C.
,
Pirker
,
S.
, and
Goniva
,
C.
,
2012
, “
Parallel Open Source CFD-DEM for Resolved Particle-Fluid Interaction
,”
Ninth International Conference on CFD in the Minerals and Process Industries
CSIRO
, Melbourne, Australia, Dec. 10–12, pp.
1
6
.http://www.cfd.com.au/cfd_conf12/PDFs/068HAG.pdf
18.
Goniva
,
C.
,
Kloss
,
C.
,
Deen
,
N. G.
,
Kuipers
,
J. A. M.
, and
Pirker
,
S.
,
2012
, “
Influence of Rolling Friction on Single Spout Fluidized Bed Simulation
,”
Particuology
,
10
(
5
), pp.
582
591
.
19.
Hill
,
R. J.
,
Koch
,
D. L.
, and
,
A. J. C.
,
2001
, “
Moderate-Reynolds-Number Flows in Ordered and Random Arrays of Spheres
,”
J. Fluid Mech.
,
448
, pp.
243
278
.
20.
Hill
,
R. J.
,
Koch
,
D. L.
, and
,
A. J. C.
,
2001
, “
The First Effects of Fluid Inertia on Flows in Ordered and Random Arrays of Spheres
,”
J. Fluid Mech.
,
448
, pp.
213
241
.
21.
Li
,
L.
,
Li
,
B.
, and
Liu
,
Z.
,
2017
, “
Modeling of Spout-Fluidized Beds and Investigation of Drag Closures Using OpenFOAM
,”
Powder Technol.
,
305
, pp.
364
376
.
22.
Gidaspow
,
D.
,
1994
,
Multiphase Flow and Fluidization: Continuum and Kinetic Theory Descriptions
,
, San Diego, CA.
23.
Beetstra
,
R.
,
Van der Hoef
,
M. A.
, and
Kuipers
,
J. A. M.
,
2007
, “
Numerical Study of Segregation Using a New Drag Force Correlation for Polydisperse Systems Derived From Lattice-Boltzmann Simulations
,”
Chem. Eng. Sci.
,
62
(
1–2
), pp.
246
255
.
24.
Celik
,
I. B.
,
Ghia
,
U.
,
Roache
,
P. J.
,
Freitas
,
C. J.
,
Coleman
,
H.
, and