Use of Precise Area Fraction Model for Fine Grid DEM Simulation of ICFB with Large Particles

: The heterogeneous structures in a gas–solid ﬂuidized bed can be resolved in discrete element simulation so long as the grid is ﬁne enough. In order to conveniently calculate mean porosity in ﬁne grid simulations, a precise area fraction model is given for two-dimensional simulations. The proposed area fraction model is validated by the discrete element simulation test on a small-scale internal circulation ﬂuidized system of large particles, using a ﬁne grid size of two particle diameters. Simulations show that the discrete element method can perform well in modelling time-varying waveforms for the physical quantities in an internal circulating ﬂuidized bed, employing the precise gas area fraction model. This thought of precise calculation can be generalized to construct a volume fraction porosity model for three-dimensional simulation by use of the similar symmetry of a rectangular grid. Moreover, to construct these area and volume fraction models is to enrich and perfect the underlying model of ﬁne grid simulation.


Introduction
Internal circulating fluidized beds (ICFB) can find lots of applications in industrial fields such as energy, environment, chemical engineering, etc. [1]. The non-intrusive access of an ICFB is too difficult to hamper the detailed experimental study of granulation. As a complement to experimental study, numerical methods are being used more and more frequently for simulating gas-solid two-phase flow behaviors and transfer characteristics which are significant for the design of an ICFB. Among these methods, the discrete element method (DEM) [2][3][4][5] deals with the solid phase as discrete particles, and thus is advantageous in providing the particle level information. DEM has become a powerful tool in fluidized bed researches and has shown good prospects.
However, the reliable applications of DEM simulations for various fluidized operation types need to be further validated in many ways. Now DEM can find its mature applications in simulations of fundamental fluidization phenomena and analysis of qualitative tendencies, rather than in quantitative researches. Since parameters for both gas and solid phases are calculated within grid circumstance in DEM simulation, results have much to do with the grid size chosen. For example, if one uses a coarse grid in a common DEM simulation code, the distinct heterogeneous particle distribution within the grid will cause large deviation in the calculated drag force [6,7]. On the contrary, fine grid simulation, although it has high spatial resolution, will lead to the complexity of parameter calculation.
Nowadays more and more researchers tend to use a fine grid in their DEM simulations [8][9][10][11]. It can be easily imagined that a fine grid is advantageous in reducing the amount of smoothing of the flow field. It is suggested that the fine grid size of 2-4 particle diameters be sufficient to resolve the complex structures in fluidized beds of various types of particles [12]. Moreover, the main issue of

Precise Area Fraction Model
In general, two-dimensional grid porosity, i.e., the gas area fraction, can be expressed as where ε 2p is the solid area fraction, f i is the grid area occupied by particle i and N is the total number of particles overlapped with the grid. For fine grid simulation, a detailed check of overlap should be made and thus the precise expression of f i is complicated. In the following, we will give a precise area fraction (PAF) model. Figure 1 shows all five cases of particle-grid overlap. For convenience of editing a model code, the five cases when the particle center is relatively close to a certain grid vertex should be separately analyzed. Figure 2 shows the region divisions of the particle-grid overlap. If the particle center is in the total region noticed, then f i > 0. The lower-left one fourth of the total region determined by the grid center (x, y) is also noticed. If the particle center is in this region, it is relatively close to the lower-left grid vertex. Moreover, the lower-left region is further divided into five sub-regions which correspond to the five cases of grid-particle overlap.   According to Figure 2, when the particle center is relatively close to the lower-left grid vertex, there are five cases of grid-particle overlap. Here it is assumed that a square grid is used, and the grid side is at least two times longer than the particle diameter. As will be seen, the model can be easily generalized to a rectangular grid because of similar symmetry. Figure 3 gives the first two cases. For Case 1, the particle center ( , ) satisfies − < < and − < < . Then is calculated as For Case 2, it satisfies that − < < , − < < and ( − ) + ( − ) < .
Then is calculated as Figure 2. Region division of grid-particle overlap.
According to Figure 2, when the particle center is relatively close to the lower-left grid vertex, there are five cases of grid-particle overlap. Here it is assumed that a square grid is used, and the grid side is at least two times longer than the particle diameter. As will be seen, the model can be easily generalized to a rectangular grid because of similar symmetry. Figure 3 gives the first two cases. For Case 1, the particle center (x, y) satisfies x 1 − r < X < x and y 1 − r < Y < y. Then f i is calculated as For Case 2, it satisfies that Then f i is calculated as  . Cases 1 and 2 when particle center is relatively close to lower-left grid vertex. Figure 4 shows Case 3 when the particle center is relatively close to the lower-left grid vertex. For this case, there are four variations for which the particle center satisfies ( − ) + ( − ) < . However, the expressions of for the four variations are formally consistent. According to Figure  4, . Cases 1 and 2 when particle center is relatively close to lower-left grid vertex. Figure 4 shows Case 3 when the particle center is relatively close to the lower-left grid vertex.
For this case, there are four variations for which the particle center satisfies (X − x 1 ) 2 + (Y − y 1 ) 2 < r. However, the expressions of f i for the four variations are formally consistent. According to Figure 4, Then f i is calculated as Symmetry 2020, 12, x FOR PEER REVIEW 5 of 15 For the second variation the particle center satisfies + < < , < < + . Then the expression of is  Figure 5 shows Case 4 when the particle center is relatively close to the lower-left grid vertex. For this case, there are two variations. For the first variation the particle center satisfies x 1 < X < x 1 + r and y 1 + r < Y < y. Then the expression of f i is For the second variation the particle center satisfies x 1 + r < X < x, y 1 < Y < y 1 + r. Then the expression of f i is   Figure 6 shows Case 5 when the particle center is closer to the lower-left grid vertex. For this case, there are also two variations. For the first variation, the particle center satisfies < < , − < < . The expression of is For the second variation, the particle satisfies − < < , < < . Then the expression of is For the above five cases, is calculated provided that the lower-left grid vertex is relatively close to the considered particle. If the considered particle is relatively close to another grid vertex,  Figure 6 shows Case 5 when the particle center is closer to the lower-left grid vertex. For this case, there are also two variations. For the first variation, the particle center satisfies Symmetry 2020, 12, x FOR PEER REVIEW 6 of 15  Figure 6 shows Case 5 when the particle center is closer to the lower-left grid vertex. For this case, there are also two variations. For the first variation, the particle center satisfies For the second variation, the particle satisfies − < < , < < . Then the expression of is For the above five cases, is calculated provided that the lower-left grid vertex is relatively close to the considered particle. If the considered particle is relatively close to another grid vertex, For the second variation, the particle satisfies For the above five cases, f i is calculated provided that the lower-left grid vertex is relatively close to the considered particle. If the considered particle is relatively close to another grid vertex, f i can be similarly calculated. Put f i into Equation (1) and calculate the precise ε 2p and ε 2D , then we construct the PAF model.

Extension to Three-Dimensional Simulation
We have constructed a three-dimensional porosity model which considers all the cases of particle-grid overlap. Let us take just a particle to be a sphere and a grid to be a cube. Figure 7 gives the schematic of all nine cases of particle-grid overlap. For Case 2 to Case 4, the solid volume can be obtained by precisely calculating the volume of the sphere minus the volume of one, two or three spherical segments, respectively. However, generally the solid volume is complicated and cannot be precisely calculated for the rest of the five cases. can be similarly calculated. Put into Equation (1) and calculate the precise ε and ε , then we construct the PAF model.

Extension to Three-Dimensional Simulation
We have constructed a three-dimensional porosity model which considers all the cases of particle-grid overlap. Let us take just a particle to be a sphere and a grid to be a cube. Figure 7 gives the schematic of all nine cases of particle-grid overlap. For Case 2 to Case 4, the solid volume can be obtained by precisely calculating the volume of the sphere minus the volume of one, two or three spherical segments, respectively. However, generally the solid volume is complicated and cannot be precisely calculated for the rest of the five cases. The intersected geometries in Case 5 and 9 are called quasi semi-segment and quasi quartersegment, respectively. Then the intersected solid volume for all cases can be uniformly expressed as (11) where k and n are equal to 0 or 1, l and m are equal to 0, 1, 2 or 3. As is known, Vsphere and Vsegment can be easily and precisely calculated. In fact, one can also directly calculate Vquasi semi-segment and Vquasi quartersegment by use of a single mathematical formula, employing the double definite integral and compound numerical integral. See [13] for more details.

Simulation Methods
Gas motion over grid k is determined by the Navier-Stokes equations as where is the gas density, is the gas velocity, is the gas pressure, is the velocity of particle , is the viscous stress tensor, is the momentum exchange coefficient [14], and is the number of particles within the grid.
The last term in Equation (13) is the coupling term. This way of coupling, although less frequently adopted by researchers, seems to be more capable of predicting bed expansion [7,15]. Xu and Yu [4] argued that the coupling between gas-solid two phases should follow Newton's third law of motion. That is to say, drag force on individual particles should react on gas from individual particles. Although the strictly coupling ways is widely used, according to Feng and Yu [16], the difference between the simulations is less significant unless a large particle layer keeps fluidized in the flow region close to the inlet. As the gas-solid interaction scale is thought to be larger than the fine grid scale, we forgo the use of the frequently-adopted coupling way. The intersected geometries in Case 5 and 9 are called quasi semi-segment and quasi quarter-segment, respectively. Then the intersected solid volume for all cases can be uniformly expressed as (11) where k and n are equal to 0 or 1, l and m are equal to 0, 1, 2 or 3. As is known, V sphere and V segment can be easily and precisely calculated. In fact, one can also directly calculate V quasi semi-segment and V quasi quarter-segment by use of a single mathematical formula, employing the double definite integral and compound numerical integral. See [13] for more details.

Simulation Methods
Gas motion over grid k is determined by the Navier-Stokes equations as where ρ g is the gas density, u is the gas velocity, p is the gas pressure, v i is the velocity of particle i, τ g is the viscous stress tensor, β is the momentum exchange coefficient [14], and N k is the number of particles within the grid. The last term in Equation (13) is the coupling term. This way of coupling, although less frequently adopted by researchers, seems to be more capable of predicting bed expansion [7,15]. Xu and Yu [4] argued that the coupling between gas-solid two phases should follow Newton's third law of motion. That is to say, drag force on individual particles should react on gas from individual particles. Although the strictly coupling ways is widely used, according to Feng and Yu [16], the difference between the simulations is less significant unless a large particle layer keeps fluidized in the flow region close to the inlet. As the gas-solid interaction scale is thought to be larger than the fine grid scale, we forgo the use of the frequently-adopted coupling way. The gas area fraction ε 2D of a grid is calculated by the proposed PAF model. Note that drag correlation is based on the real three-dimensional experimental results. The grid porosity in the Navier-Stokes equations should be three-dimensional porosity ε 3D which is ε g in Equations (13) and (14). Then ε 2D is transformed into ε 3D [4] as Equations (12) and (13) are discretized by the finite volume method based on the collocated grid. Uniform velocity inlet, pressure outlet and impermeable wall are defined as the boundary conditions. The Navier-Stokes equations are solved according to the method in [17].
The soft-sphere model is used to deal with particle-particle and particle-wall collisions. Presently the time step for the gas motion is 0.00005 s and the DEM time step is 0.00001 s. This DEM time step is much smaller than the required maximum time resolution [2] defined by the spring constant.
The translational motion of each particle i is determined by where ρ p is particle density, V p is particle volume, F Di is drag force, F Ci is collision force and p i is local gas pressure. The rotational motion of particle i is described as where ω i is the particle angular velocity, I is the inertia moment of the particle as spherical and T Ci is the torque due to collision. The drag force F Di is computed by the Wen and Yu correlation [14], using a so-called complete circumstance-independent model to calculate local porosity [7], as where u i is the local gas velocity which is calculated by area weighted averaging, and the standard drag coefficient C Di is expressed as The Reynolds number of particle i in Equation (18) is In the present simulations, the parameters involved are given in Table 1, with the geometry, operation conditions and material properties basically similar to those in the experimental research of van Wachem et al. [18].  Figure 8 demonstrates the particle distributions in the present simulations and previous experiments [18]. The simulated bubble size is in agreement with that in the experiments. The inconsistency of bubble shape may indicate that the real complicated inlet effect cannot be properly reflected by just setting the simple inlet boundary condition. In the simulations the bubble wake effect is stronger than in the experiments. Thus at 0.42 s, the bubble seems bullet-bottom; while at 0.5 s, there are much more particles thrown into the big bubble.  Figure 8 demonstrates the particle distributions in the present simulations and previous experiments [18]. The simulated bubble size is in agreement with that in the experiments. The inconsistency of bubble shape may indicate that the real complicated inlet effect cannot be properly reflected by just setting the simple inlet boundary condition. In the simulations the bubble wake effect is stronger than in the experiments. Thus at 0.42 s, the bubble seems bullet-bottom; while at 0.5 s, there are much more particles thrown into the big bubble.

Solid Volume Fraction
It has been reported that traditional coarse grid DEM simulations may capture some features which deviate from experiments [18]. Drag calculation was improved and better simulations were obtained in our coarse grid simulation work [6,7]. Those previous simulation results should be referred to, but are not given in the following. Only the present simulation results and the pertinent experimental results are given for comparison to validate the present fine grid simulations. Figure 9 demonstrates the measured and calculated solid volume fraction fluctuations which are averaged in the horizontal plane at 45 mm above the distributor. Although small amplitude fluctuations can be seen in both the experiments and the present simulations, there are more of them in the simulations. Thus, the simulated fluctuations seem a little faster than the measured. The former traverse about 9 cycles in 5-7 s while the latter traverse about 7.5 cycles in 2-4 s. The simulated solid volume fraction fluctuation slightly deviates from the waveform shown in the experiments. This deviation may partly result from the calculation method for the solid volume fraction in the simulations. The solid volume fraction is computed as one plus ε . However, the transformation of ε to ε is adopted, which to some extend may deviate from the real value of porosity, although this transformation formula is widely used. On the whole, the simulated waveform for the solid volume fraction is as good as that in [7], which is much more consistent with the experiment data than that in [6,18].

Solid Volume Fraction
It has been reported that traditional coarse grid DEM simulations may capture some features which deviate from experiments [18]. Drag calculation was improved and better simulations were obtained in our coarse grid simulation work [6,7]. Those previous simulation results should be referred to, but are not given in the following. Only the present simulation results and the pertinent experimental results are given for comparison to validate the present fine grid simulations. Figure 9 demonstrates the measured and calculated solid volume fraction fluctuations which are averaged in the horizontal plane at 45 mm above the distributor. Although small amplitude fluctuations can be seen in both the experiments and the present simulations, there are more of them in the simulations. Thus, the simulated fluctuations seem a little faster than the measured. The former traverse about 9 cycles in 5-7 s while the latter traverse about 7.5 cycles in 2-4 s. The simulated solid volume fraction fluctuation slightly deviates from the waveform shown in the experiments. This deviation may partly result from the calculation method for the solid volume fraction in the simulations. The solid volume fraction is computed as one plus ε 3D . However, the transformation of ε 2D to ε 3D is adopted, which to some extend may deviate from the real value of porosity, although this transformation formula is widely used. On the whole, the simulated waveform for the solid volume fraction is as good as that in [7], which is much more consistent with the experiment data than that in [6,18].

Solid Volume Fraction
It has been reported that traditional coarse grid DEM simulations may capture some features which deviate from experiments [18]. Drag calculation was improved and better simulations were obtained in our coarse grid simulation work [6,7]. Those previous simulation results should be referred to, but are not given in the following. Only the present simulation results and the pertinent experimental results are given for comparison to validate the present fine grid simulations. Figure 9 demonstrates the measured and calculated solid volume fraction fluctuations which are averaged in the horizontal plane at 45 mm above the distributor. Although small amplitude fluctuations can be seen in both the experiments and the present simulations, there are more of them in the simulations. Thus, the simulated fluctuations seem a little faster than the measured. The former traverse about 9 cycles in 5-7 s while the latter traverse about 7.5 cycles in 2-4 s. The simulated solid volume fraction fluctuation slightly deviates from the waveform shown in the experiments. This deviation may partly result from the calculation method for the solid volume fraction in the simulations. The solid volume fraction is computed as one plus ε . However, the transformation of ε to ε is adopted, which to some extend may deviate from the real value of porosity, although this transformation formula is widely used. On the whole, the simulated waveform for the solid volume fraction is as good as that in [7], which is much more consistent with the experiment data than that in [6,18].   Figure 10 demonstrates the measured and presently-calculated relative pressure fluctuations which are averaged in the horizontal plane at 45 mm above the distributor. Both the simulated and measured relative pressure values are between -200-200 Pa. As is noticed, small amplitude fluctuations of relative pressure occur in both the experiments (from 3 to 3.75 s) and the present simulations (from 5.5 to 6.25 s). This reflects the random and instable features of the real fluidization. The unconformity of the two frames in the figure lies in the waveforms. The simulated waveform reveals weaker periodicity. Pressure is a basic variable of the Navier-Stokes equations. It is hard to be accurately calculated, which mainly lies not in calculation error, but in model error. The momentum exchange source term in Equation (13) represents the gas-particle interaction which is calculated by the widely used correlation formula proposed by Wen and Yu [14]. This model error may be reduced in two ways. One way is to develop structure-dependent drag force in a coarse grid simulation [6]. Another way is to use a fine grid so that the heterogeneous structure can be fully resolved, which is also adopted in this work. Without consideration of the small amplitude fluctuations during the similar period of 0.75 s, the calculated fluctuations traverse about 6.5 cycles while the measured fluctuations traverse about 5.5 cycles. In this sense, the simulated fluctuation frequency is better than that in [7].  Figure 10 demonstrates the measured and presently-calculated relative pressure fluctuations which are averaged in the horizontal plane at 45 mm above the distributor. Both the simulated and measured relative pressure values are between -200-200 Pa. As is noticed, small amplitude fluctuations of relative pressure occur in both the experiments (from 3 to 3.75 s) and the present simulations (from 5.5 to 6.25 s). This reflects the random and instable features of the real fluidization. The unconformity of the two frames in the figure lies in the waveforms. The simulated waveform reveals weaker periodicity. Pressure is a basic variable of the Navier-Stokes equations. It is hard to be accurately calculated, which mainly lies not in calculation error, but in model error. The momentum exchange source term in Equation (13) represents the gas-particle interaction which is calculated by the widely used correlation formula proposed by Wen and Yu [14]. This model error may be reduced in two ways. One way is to develop structure-dependent drag force in a coarse grid simulation [6]. Another way is to use a fine grid so that the heterogeneous structure can be fully resolved, which is also adopted in this work. Without consideration of the small amplitude fluctuations during the similar period of 0.75 s, the calculated fluctuations traverse about 6.5 cycles while the measured fluctuations traverse about 5.5 cycles. In this sense, the simulated fluctuation frequency is better than that in [7].

Bed Layer Height
Since in the study conducted by van Wachem et al. [18], the expression of bed layer height was not given, we use a visual observation technique as in [7], rather than the on-line output technique as in [6]. Figure 11 shows the simulated particle locations at the time nodes of the locally highest and lowest bed layer height. Figure 12 shows the measured and presently-simulated bed layer height expansion. Again, it is noticed that the waveforms are the most distinct difference. The simulated time-varying seems like a piecewise linear function of time because it is approximated by just using local extremes. In [7], both the mean minimum and maximum values of the simulated are a little higher than those in the experiments. However, in the present simulations, both the mean minimum and maximum values of the simulated are very close to those in the experiment. The simulated fluctuations traverse about five cycles in 6-8 s while the measured fluctuations traverse about 5.5 cycles in 2-4 s. The simulated fluctuation frequency is only a little lower than the measured fluctuation frequency. These results show that the bed expansion in the present simulations is better than that in the previous simulations [7]. Moreover, the tendency of slow expansion and fast fallback, shown in the experiments, can be well described in the present simulations.

Bed Layer Height
Since in the study conducted by van Wachem et al. [18], the expression of bed layer height H was not given, we use a visual observation technique as in [7], rather than the on-line output technique as in [6]. Figure 11 shows the simulated particle locations at the time nodes of the locally highest and lowest bed layer height. Figure 12 shows the measured and presently-simulated bed layer height expansion. Again, it is noticed that the waveforms are the most distinct difference. The simulated time-varying H seems like a piecewise linear function of time because it is approximated by just using local extremes. In [7], both the mean minimum and maximum values of the simulated H are a little higher than those in the experiments. However, in the present simulations, both the mean minimum and maximum values of the simulated H are very close to those in the experiment. The simulated fluctuations traverse about five cycles in 6-8 s while the measured fluctuations traverse about 5.5 cycles in 2-4 s. The simulated fluctuation frequency is only a little lower than the measured fluctuation frequency. These results show that the bed expansion in the present simulations is better than that in the previous simulations [7]. Moreover, the tendency of slow expansion and fast fallback, shown in the experiments, can be well described in the present simulations.

Conclusions
This article suggests that the use of a fine grid is beneficial to resolve the complex heterogeneous structures in fluidized beds of various types of particles. The major complexity of fine grid simulation is thought to lie in the calculation of grid parameter. The article mainly presents a porosity model and validates it in two-dimensional simulation of an ICFB. The following conclusions can be obtained.
(1) The PAF model is given to precisely calculate the solid or gas area fraction, and theoretically to more properly calculate the grid porosity. (2) The simulated big bubble is generally consistent with the experiment results in shape and size.
The present simulation works well enough to model the basic bubbling phenomenon of an ICFB.

Conclusions
This article suggests that the use of a fine grid is beneficial to resolve the complex heterogeneous structures in fluidized beds of various types of particles. The major complexity of fine grid simulation is thought to lie in the calculation of grid parameter. The article mainly presents a porosity model and validates it in two-dimensional simulation of an ICFB. The following conclusions can be obtained.
(1) The PAF model is given to precisely calculate the solid or gas area fraction, and theoretically to more properly calculate the grid porosity. (2) The simulated big bubble is generally consistent with the experiment results in shape and size.
The present simulation works well enough to model the basic bubbling phenomenon of an ICFB.

Conclusions
This article suggests that the use of a fine grid is beneficial to resolve the complex heterogeneous structures in fluidized beds of various types of particles. The major complexity of fine grid simulation is thought to lie in the calculation of grid parameter. The article mainly presents a porosity model and validates it in two-dimensional simulation of an ICFB. The following conclusions can be obtained.
(1) The PAF model is given to precisely calculate the solid or gas area fraction, and theoretically to more properly calculate the grid porosity. (2) The simulated big bubble is generally consistent with the experiment results in shape and size.
The present simulation works well enough to model the basic bubbling phenomenon of an ICFB.