How Computational Grid Refinement in Three Dimensions Affects Computational Fluid Dynamics-Discrete Element Method Results for Psuedo-Two-Dimensional Fluidized Gas–Solid Beds

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.

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.

where F₀ and F₃ are given by van Buijtenen et al. [8].

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

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.

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.

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.

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.

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.

where x_j is the mass fraction, and d_j 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 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).

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.

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

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.

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.

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 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.

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.

The average height of the particle bed center of mass, $⟨hbed⟩$⁠, and the amount of variation of this height, RMS_bed, 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.

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.

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.

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.75d_p,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.

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.

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.

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.

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.

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.

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.

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

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

a =

quadratic regression coefficient

b =

linear regression coefficient

c =

constant regression coefficient

C_d =

coefficient of drag

d_p =

particle diameter

d_S =

Sauter-mean particle diameter

$dsmall,dlarge$ =

small and large particle diameters, respectively

f =

frequency of data collection

f_f,i =

fluid force on particle i

f_grav,i =

gravitational force on particle i

F_t,ij =

tangential contact force exerted on particle i by particle j

$Fijc$ =

contact force exerted on particle i by particle j

F₀, F₃ =

fitted functions for Koch-Hill drag law

h =

volume-average cell size

h_fit(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

I_i =

moment of inertia of particle i

m_i =

mass of particle i

N =

total number of computational cells

p =

pressure

R =

nondimensional cell size

R_i =

radius of particle i

Re_p =

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

R_pf =

momentum exchange between fluid and particle

s =

amount of segregation

t =

time

t_max =

time instant of final data collection

t_min =

time instant of initial data collection

u_f =

fluid velocity

v_i =

velocity of particle i

V_i =

volume of particle i

w_i =

rotational velocity of particle i

x_j =

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

β =

coefficient of inter-phase momentum transfer

ε =

void fraction

μ_f =

fluid absolute viscosity

ρ_f =

fluid density

τ =

fluid phase stress tensor

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.