A Model to Improve Granular Temperature in CFD-DEM Simulations

: CFD-DEM (computational ﬂuid dynamic-discrete element method) is a promising approach for simulating ﬂuid–solid ﬂows in ﬂuidized beds. This approach generally under-predicts the granular temperature due to the use of drag models for the average drag force. This work develops a simple model to improve the granular temperature through increasing the drag force ﬂuctuations on the particles. The increased drag force ﬂuctuations are designed to match those obtained from PR-DNSs (particle-resolved direct numerical simulations). The impacts of the present model on the granular temperatures are demonstrated by posteriori tests. The posteriori tests of tri-periodic gas–solid ﬂows show that simulations with the present model can obtain transient as well as steady-state granular temperature correctly. Moreover, the posteriori tests of ﬂuidized beds indicated that the present model could signiﬁcantly improve the granular temperature for the homogenous or slightly inhomogeneous systems, while it showed negligible improvement on the granular temperature for the signiﬁcantly inhomogeneous systems.


Introduction
Fluidized beds with fluid (gas/liquid)-solid flows are operating units adopted in chemical and energy processes. In the past decades, the computational fluid dynamic (CFD) simulations have become an effective tool in the study of fluid-solid flows in the fluidized beds. The predictability of CFD highly depends on the constitutive models. Accordingly, much attention has been focused on the development of models for CFD simulations with a higher level of accuracy.
In industrial applications, the most commonly used CFD simulation methods for fluid-solid flows are TFM (two-fluid model) and CFD-DEM (computational fluid dynamic-discrete element method). TFM treats both the fluid and solid phase as a continuum, thus the effective solid stresses due to particle collisions are required [1,2]. While in CFD-DEM, the solid phases are considered as discrete entities and the particle collisions are fully resolved [3,4]. Closure for fluid-particle interaction forces, such as the drag force, unsteady interaction force, and lift force, are required in both TFM and CFD-DEM. In general, the drag force plays a significant role in determining the hydrodynamics of the fluid-solid systems [5].
The closures for the drag force have been studied by many researchers. In addition to the widely used empirical correlations obtained from experimental studies [6][7][8], many correlations for the drag force in the homogenous fluid-particle system have been proposed through PR-DNSs (particle-resolved direct numerical simulations) over the last few decades [9][10][11][12]. Moreover, some researchers studied the effect of granular temperatures on the drag force. In recent years, the significant deviation of individual particle drag force from the average drag force has been pointed out. The mean and maximum relative deviation of the drag force can be as high as 25% and 40%, respectively, for randomly distributed stationary particles [13]. It is reported that the largest drag force on an individual particle would be twice the mean drag force calculated from traditional drag models [14].
The importance of capturing the particle velocity fluctuations was raised in predicting the core-annular structure observed in riser flows [7,15]. Yu et al. [16] recently reported that it is necessary to consider the fluctuations of the meso-scale drag force to attain correct granular temperature in coarse-grid CFD-DEM simulations. Without the consideration of drag fluctuations, traditional CFD-DEM simulations cannot correctly predict the magnitude of granular temperature even in the homogenous systems. Akiki et al. [17] found that the individual drag is dependent on the local structures of neighboring particles, and they proposed a pairwise interaction extended point-particle model for the prediction of the drag force. The accuracy of their model deteriorates as the solid volume fraction increases, and they used machine learning to predict the individual drag in moderate and dense flows [18]. Ma et al. [19] adopted a second-order structure tensor to measure the local anisotropic distribution of particles. However, their method only yields an averaged drag force for particles in a computational cell. Tenneti et al. [20] gave a fundamental understanding of the generation of the granular temperature and developed the Langevin model to correct the granular temperature of the gas-solid system. Esteghamatian et al. [21] adopted time-series analysis methods to consider the stochastic drag force. The integral time and deviation of the drag coefficient are required from PR-DNS simulations.
Overall, the research on improving granular temperatures is still limited. The implementation of most available models is complex or computationally expensive. In addition, it is noted that these models do not utilize the existent drag force difference in the same cell in CFD-DEM. Based on available PR-DNS data, we hypothesized that we could produce a correlation for the correct magnitude of the drag force difference. Then, the obtained drag force difference in the cells in CFD-DEM could be magnified to match the correct magnitude of drag force difference. Following this idea, in this paper, a simple model is proposed to improve the granular temperatures through enlarging the existent drag force difference on individual particles in the same computational cell.
The rest of the paper is organized as follows. In the next section, the model to enhance the drag force fluctuations is presented based on the magnitude of the drag fluctuations obtained from PR-DNSs. Then, posteriori validations of the proposed model on three different problems are implemented. Finally, the conclusions are summarized.

Model to Improve the Granular Temperature
In this section, the relative deviation of the drag force in the homogenous system is obtained from PR-DNSs performed by Huang et al. [13]. Then, a model to improve the granular temperature in CFD-DEM simulations is proposed.

Mean Relative Deviation of the Drag Force in Homogenous Systems
The data from PR-DNSs performed by Huang et al. [13] were used in the present study. In their simulations, the second-order accurate immersed boundary-lattice Boltzmann method (IB-LBM) proposed by Zhou and Fan [22] was adopted. Flows past the randomly distributed particles with various granular temperature were simulated. For the details of the PR-DNSs, the readers are referred to Huang et al. [13]. The magnitude of the drag fluctuation, measured by the mean relative deviation of the drag force, can be readily obtained from the PR-DNSs. The effects of φ (solid volume fraction), Re (Reynolds number), and Re T (granular temperature-based Reynolds numbers) on the mean relative deviation of drag force were considered. The mean relative deviation of the drag force, σ F d , was defined as where F d is the drag force on individual particles; the angular brackets denote the average over all particles in the computational domain. The Reynolds number, Re, is given as where v p is the particle velocity, u f is the fluid velocity, d p is the particle diameter, and ρ f and µ f are the fluid density and dynamic viscosity, respectively. The granular temperature-based Reynolds number is the ratio between fluid viscous and inertial force due to particle fluctuations, namely where T is granular temperature and is defined as In Huang et al.'s work [13], the PR-DNSs with six different solid volume fractions (φ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6) were performed for six granular temperature-based Reynolds numbers (Re T = 0, 5, 8.66, 17.3, 26 and 34.6). The Reynolds number Re was varied from practically zero to approximately 120. As Re T is much smaller than Re in the practical fluidized bed simulations [23], only the simulation setups with Re T /Re < 0.5 were selected in this work. Based on the PR-DNS results, a simple correlation for the mean relative deviation of the drag force, σ dns F d , was fitted, which is valid in the range 0.1 ≤ φ ≤ 0.6, 5 ≤ Re ≤ 120, 0 ≤ Re T /Re ≤ 0.5. The correlation for the mean relative deviation of the drag force is given as σ dns where σ dns F d ,Re T =0 is the mean relative deviation of the drag force when Re T = 0. The parameters k, c, y 0 , A, and t are functions of the solid volume fraction and are listed below: c = (0.859 + 9.26φ)/(1 + 5.366φ) (8) The PR-DNS results and corresponding fitting curves calculated by Equations (5)-(11) are shown in Figure 1. Figure 1a depicts that σ dns F d ,Re T =0 varies with Re at different φ. It can be seen that there is a notable increase of σ dns F d ,Re T =0 with the increase of φ. Moreover, in most cases, σ dns F d ,Re T =0 increases with increasing Re and it levels off at high Re. At φ = 0.6, we do observe a slight decrease of σ dns F d ,Re T =0 with increasing Re. The reasons for this phenomenon are not clear. The fact that the trends of σ dns changing with Re are different between φ = 0.5 and φ = 0.6 makes the fitting much more difficult. Hence, a relatively larger fitting error could be observed at φ = 0.5. increases with the increase of the solid volume fraction. This trend is consistent with the observations made in Figure 1a for significantly decreases with decreased solid volume fractions (see Figure 1a). The small denominator tends to magnify the scattering of the numerator.
is the ratio between the mean relative deviation of the drag force at arbitrary Re T and that at Re T = 0 under the same solid volume fraction and Reynolds number. The symbols are from PR-DNS simulations, and the fits of the PR-DNS results are represented by the corresponding solid lines computed by Equation (5). Figure 1b. It shows that σ dns increases with the increase of Re T /Re and the increasing rate increases gradually. Moreover, σ dns first decreases with the increase of the solid volume fraction and then increases slightly when the solid volume fraction rises from 0.5 to 0.6. It should be noted that σ dns F d also increases with the increase of the solid volume fraction. This trend is consistent with the observations made in Figure 1a for σ dns F d ,Re T =0 . Close observation shows that the scattering of PR-DNS results of σ dns F d /σ dns F d ,Re T =0 generally increases with decreased solid volume fractions. This is understandable since the denominator σ dns F d ,Re T =0 significantly decreases with decreased solid volume fractions (see Figure 1a). The small denominator tends to magnify the scattering of the numerator.

Model to Enhance Granular Temperature in CFD-DEM Simulations
In CFD-DEM simulations, particles located in the same fluid grid experience different drag forces due to the difference of the particle velocities, the local solid volume fractions, and local gas velocities. The drag force fluctuation of the particle in i direction, F d,i , is where F d,i is the average drag component in i direction (x, y or z direction) on the particles located in the same computational grid. The overbar denotes the average in a computational cell. The drag force fluctuation F d,i obtained directly from CFD-DEM simulations is generally much lower than that computed from PR-DNSs. To make the magnitude of F d,i approach that from PR-DNSs, the modified drag force fluctuation is calculated as follows In Equation (13) (13). The determination of h can be achieved by simple iterative methods introduced in textbooks of numerical analysis, such as the shooting method. The solution is believed to be converged when the relative difference between the mean relative deviation calculated based on F d,i and σ exp F d is less than 1%. Setting the initial value of h as unity, the shooting method usually converges rapidly in only several iterations. The converged h usually attains the value ranging from 2 to 10, making the magnitude of F d,i generally much larger than that of F d,i . With the final value of h determined in the iteration, the modified particle drag force, can still be obtained even when the used coordinates have an arbitrary angle with the original framework.

Determination of the Expected Mean Relative Deviation of the Drag Force
The abovementioned expected drag deviation σ exp F d is calculated through the following expression: where cor is the correction coefficient. To obtain the correction coefficient, gas-solid flows suspended in a tri-periodic domain are simulated using MFIX-DEM (https://mfix.netl.doe.gov/) [24]. The computational domain is a cube with side lengths of L = 9d p . In these simulations, only one computational cell is used, and the particle drag forces are given by Equation (15). The steady-state granular temperature is obtained after the simulations reach thier statistical steady state. It should be noted that all the particles in these simulations are fully elastic and hence the restitution coefficients for colliding particles are set unity. Gas-solid flows with 100 ≤ ρ s /ρ f ≤ 2000, 0.1 ≤ φ ≤ 0.4 and 10 ≤ Re ≤ 100 are simulated. For a given set of ρ s /ρ f , φ and Re, the correction coefficient, cor, is adjusted to ensure the steady-state granular temperature attains the same value as that calculated by the formula proposed by Tenneti et al. [20]. Consequently, the correlation for the correction coefficient with ρ s /ρ f , φ and Re is fitted and is given as g(φ) = −1.039 exp(−φ/0.066) + 1.071 (19) h(Re) = 0.272 exp(−Re/21.126) + 0.938 (20)

Posteriori Validations
This section is to assess the performance of the proposed model in improving the granular temperature. In the MFIX code, the particle drag forces and the correction coefficient of the mean relative deviation of drag force are given by Equations (15) and (17), respectively. Posteriori validations on three different types of flows are performed: gas-solid flows in a tri-periodic domain, 3-dimensional liquid-solid fluidized beds, and 3-dimensional gas-solid fluidized bed.

Gas-Solid Flows in a Tri-Periodic Domain
CFD-DEM simulations of gas-solid flows suspended in a tri-periodic domain were performed over a wide range of solid volume fractions (0.1 ≤ φ ≤ 0.4), Reynolds numbers (10 ≤ Re ≤ 100), and density ratios (100 ≤ ρ s /ρ f ≤ 2000). The average drag model proposed by Tenneti et al. [11] was used. Note that these PR-DNS results were used to attain the coefficient in Equation (16) in Section 2.3. The simulations run for 100τ p to achieve the statistical steady state and then run for another 100τ p to gain statistical values. τ p is the particle time defined as τ p = d p / v p − u f , where v p and u f are the mean particle velocity and mean fluid velocity in the periodic domain, respectively.
The steady-state granular temperatures obtained from simulations are compared with those calculated by the formula proposed by Tenneti et al. [20] in Figure 2. It can be seen that the steady-state granular temperatures obtained from simulations agree with their formula very well. This demonstrates the effectiveness of the proposed model in improving the prediction of granular temperature in CFD-DEM simulations. Besides, the results show that the granular temperature scaled by the slip velocity decreases with the increase of the solid-fluid density ratio, Reynolds number, and solid-volume fraction. Close observation shows that the deviation between the proposed model and CFD-DEM increased with decreased density ratio. This may be because the proposed formula of the mean relative deviation of the drag force (Equations (5) and (6)) came from simulations of fixed particles in which the Stokes number (S t = Reρ s /18ρ f ) was high. The observed deviation could be explained considering that the Stokes number decreases with the decrease of density ratio. It is also observed that the deviation between the proposed model and CFD-DEM slightly increased with a decreased solid volume fraction. The physical reason for this is not clear. One possible reason is that the proposed correction coefficients (Equation (17)) are not well fitted at small solid volume fractions.
The impact of the present model on the evolution behavior of granular temperature was also evaluated. The evolution of granular temperatures under different initial granular temperatures were compared with the transient PR-DNS results obtained in Tenneti et al. [20], and are shown in Figure 3. The data from the case with ρ s /ρ f = 100, φ= 0.1 and Re= 10 were selected as an illustrative example. It can be seen that the evolution of granular temperatures from the CFD-DEM simulations with the present model agrees favorably well with that from the PR-DNSs by Tenneti et al. [20]. Moreover, for cases with different initial granular temperatures, the converged steady-state granular temperatures from the present model reach almost the same value, which is very close to that obtained from the PR-DNSs. It is noted that the granular temperature from traditional CFD-DEM simulation (without the present model) decays to zero. Overall, this test showed that CFD-DEM with the present model has the capacity to correctly produce the transient as well as the steady-state granular temperature in small computational domains.
Energies 2020, 13, 4730 7 of 12 CFD-DEM increased with decreased density ratio. This may be because the proposed formula of the mean relative deviation of the drag force (Equations (5) and (6)) came from simulations of fixed particles in which the Stokes number (

Re 18
ρ ρ = t s f S ) was high. The observed deviation could be explained considering that the Stokes number decreases with the decrease of density ratio. It is also observed that the deviation between the proposed model and CFD-DEM slightly increased with a decreased solid volume fraction. The physical reason for this is not clear. One possible reason is that the proposed correction coefficients (Equation (17)) are not well fitted at small solid volume fractions. The impact of the present model on the evolution behavior of granular temperature was also evaluated. The evolution of granular temperatures under different initial granular temperatures were compared with the transient PR-DNS results obtained in Tenneti et al. [20], and are shown in Figure 3.
The data from the case with were selected as an illustrative example. It can be seen that the evolution of granular temperatures from the CFD-DEM simulations with the present model agrees favorably well with that from the PR-DNSs by Tenneti et al. [20]. Moreover, for cases with different initial granular temperatures, the converged steady-state granular temperatures from the present model reach almost the same value, which is very close to that obtained from the PR-DNSs. It is noted that the granular temperature from traditional CFD-DEM simulation (without the present model) decays to zero. Overall, this test showed that CFD-DEM with  [20]. The granular temperature is scaled by the mean fluid-particle slip velocity W = v p − u f , where the angular brackets denote the average over the system. The solid-fluid density ratio for (a), (b), and (c) are 100, 500, and 1000, respectively. As the agreements between scatters and lines were similar for the density ratios 2000 and 1000, the comparisons are only shown for density ratio at 1000.    Re T−initial is the initial granular temperature-based Reynolds number imposed on the particles. For details of imposing initial granular temperature on the particles, please refer to Huang et al. [13]. It should be noted that the scatters from CFD-DEM without the new model for Re T−initial = 0.11 are not shown, as the corresponding granular temperatures were too small.

Liquid-Solid Fluidized Beds
CFD-DEM simulations were performed for liquid-solid flows in bi-periodic fluidized beds with periodic boundary condition imposed both in the x and y directions. In these simulations, Energies 2020, 13, 4730 8 of 12 the vertical direction is the z direction and the lateral direction is the x or y direction. The size of fluidized bed was 15d p in the vertical direction and 8d p in the both lateral directions. The liquids were uniformly injected from the bottom of the beds with inlet velocity U in , while the top of the beds was imposed by the constant pressure boundary condition. The simulation parameters are listed in Table 1. The computational grid size of 2d p was selected in the following simulations. Following the suggestion of Esteghamatian et al. [5], model B was adopted in the present CFD-DEM simulations. The details of the difference between models A and model B in the CFD-DEM simulations was not the focus of the present study and hence will not be explained here. The interested reader can refer to Zhou et al. [25]. The bounding cube (BC) scheme [26] was used to map information from particles to the computational grid and vice-versa. The average drag model proposed by Beetstra et al. [9] was used in the following simulations. Table 1. Simulation parameters of liquid-solid fluidized beds. Note: µ f is the dynamic viscosity of liquid. U m f is the minimum fluidized fluid velocity derived from Gidaspow [7].

Simulation Parameters Values
The evolution of granular temperature in the main and lateral flow directions of the CFD-DEM simulations were compared with the PR-DNS results from Esteghamatian et al. [21]. The granular temperature in the main and lateral directions are defined as where the main flow direction (vertical direction) is along the positive z axis. The flows in this fluidized bed are statistically the same in the x and y directions, so only the granular temperature in the x direction was selected. v p,z and v p,x are particle velocities in the z and x directions, respectively. The evolutions of granular temperatures under various inlet Reynolds numbers are shown in Figure 4. It can be seen that the granular temperature from CFD-DEM without the proposed model decays to very small values after the systems reach their steady states, while the granular temperatures from CFD-DEM with the proposed model are significantly improved. Close observation shows that, compared with the PR-DNSs, CFD-DEM simulations with the proposed model overpredict/underpredict the granular temperature at low/high inlet Reynolds numbers.
that the ratio of the steady-state granular temperature in the vertical direction to that in the lateral direction increases with increasing Reynolds number. This indicates that the inhomogeneity of the flows would increase with the increase of Reynolds number. A smaller computational grid size is required to resolve the existing inhomogeneity. The computational grid size of 2 p d used in this work was not sufficient for resolving the inhomogeneity. Thus, the unresolved inhomogeneity leads to the under-prediction of the granular temperature at relatively high Reynolds numbers.
Energies 2020, 13, x FOR PEER REVIEW 10 of 13  [21] are also provided. The granular temperatures are scaled by the square of the inlet fluid velocity.

Gas-Solid Fluidized Beds
CFD-DEM simulations were performed for gas-solid flows in bi-periodic fluidized beds with periodic boundary conditions imposed both in the x and y directions. In these simulations, the vertical direction was the z direction and the lateral direction was the x or y direction. The size of the fluidized bed was 50 p d in the vertical direction and 10 p d in the both lateral directions. The gas was At relatively small Reynolds numbers (Re = 4), the proposed model over-predicts the granular temperature. This may be because Equation (16) over-predicts the expected drag deviations. After all, the proposed drag deviations are developed based on simulation results from nearly homogeneous systems. They are expected to work well only for the homogenous gas-solid systems in Section 2.3. Nevertheless, compared to the traditional CFD-DEM, the CFD-DEM with the proposed model significantly improves the prediction of granular temperature. At relatively high Reynolds numbers, the proposed model under-predicts the granular temperature. This is because the inhomogeneity becomes important at relatively high Reynolds numbers. In Figure 4, it is obvious that the ratio of the steady-state granular temperature in the vertical direction to that in the lateral direction increases with increasing Reynolds number. This indicates that the inhomogeneity of the flows would increase with the increase of Reynolds number. A smaller computational grid size is required to resolve the existing inhomogeneity. The computational grid size of 2d p used in this work was not sufficient for resolving the inhomogeneity. Thus, the unresolved inhomogeneity leads to the under-prediction of the granular temperature at relatively high Reynolds numbers.

Gas-Solid Fluidized Beds
CFD-DEM simulations were performed for gas-solid flows in bi-periodic fluidized beds with periodic boundary conditions imposed both in the x and y directions. In these simulations, the vertical direction was the z direction and the lateral direction was the x or y direction. The size of the fluidized bed was 50d p in the vertical direction and 10d p in the both lateral directions. The gas was uniformly injected from the bottom of the beds with inlet velocity U in , while the top of the beds was imposed by the constant pressure boundary condition. The computational grid size was 2d p and the bounding cube (BC) scheme was used. The simulation parameters are listed in Table 2. The simulations ran for 200τ p to reach a steady-state and ran another 500τ p to obtain statistical results. In Table 3, the steady-state granular temperatures of the CFD-DEM simulations are compared with those of the PR-DNS results [21]. It can be seen that the magnitude of the granular temperature without the present model would be under-predicted in both directions. Moreover, with the present model, the granular temperature was under-predicted in the vertical direction while it was over-predicted in the lateral direction. The present model only showed better performance in terms of the average granular temperature of three directions. On the whole, the CFD-DEM simulation with the present model performs as well as but not better than that without the present model. This may be because of the existence of significant inhomogeneity in the gas-solid flows, which was mentioned in Esteghamatian et al.'s work [21]. The PR-DNS results in Table 3 show the granular temperature in the vertical direction is around three times that in the lateral direction. This also indicates the significant inhomogeneity in the flows. As explained in Section 3.2, the unresolved inhomogeneity would lead to under-prediction of the granular temperature. Table 3. Comparisons of the steady-state granular temperature between CFD-DEM simulations and PR-DNSs [21].  [21]. The values in the brackets are the relative errors compared with the PR-DNS results.

Conclusions
CFD-DEM generally under-predicts granular temperature due to the use of the average drag force. To improve the prediction of granular temperature, the present work proposed a model to enhance the drag force fluctuations on particles (see Equation (15)). The mean relative deviations of the drag force from PR-DNS simulations were calculated and used as a reference to determine the modified drag force.
The present model was implemented in the MFIX-DEM software. To assess the impact of the present model on the CFD-DEM simulations, three types of fluid-solid systems were tested. In the systems of gas-solid flows suspended in a tri-periodic domain, the steady-state granular temperatures calculated by CFD-DEM agree favorably well with the formula proposed by Tenneti et al. [20]. Besides, the evolution of the granular temperature also matches the transient solution from PR-DNSs [20]. For the liquid-solid systems, the CFD-DEM without the present model significantly under-predicts the granular temperatures, and there is an appreciable improvement of the granular temperatures when the present model is used. It was also found that the inhomogeneity increases with an increasing Reynolds number and the unresolved inhomogeneity would lead to the under-prediction of the granular temperature. Though the agreement with PR-DNS results [21] was not as good as that in the first validation case, the effectiveness of the present model in improving the prediction of granular temperature was also demonstrated. For the gas-solid systems, CFD-DEM simulation with the present model was unable to a give better performance than that without the present model. This may also be because of the existence of the significant inhomogeneity. In conclusion, CFD-DEM with the proposed model is able to consider the effect of the unresolved drag force fluctuations on the granular temperature, while the effect of the unresolved inhomogeneity on the granular temperature cannot be considered.
Overall, the present study shows that granular temperature calculated from CFD-DEM simulations could be profoundly improved through simply magnifying the existent drag force fluctuation on different particles in the same computational cell. Future studies should be focusing on further improvement of the proposed model. For example, the estimation of the coefficients in Equation