Verification and Validation of Spatially Averaged Models for Fluidized Gas‐Particle Suspensions

Two different methods for closure modeling of the unresolved terms appearing in the filtered two-fluid model are discussed and compared. The spatially averaged two-fluid model is based on generalizing the concepts of large eddy simulation to gas-particle flows. In the approximate deconvolution method-two-fluid model approach, the unresolved terms are modeled by an approximate deconvolution method, where an approximation of the unfiltered solution is obtained by repeated filtering. Finally, these models are applied to a lab-scale and a pilot-scale fluidized bed. Both approaches yield fairly good agreement with a highly resolved reference simulation as well as with experimental data. Additionally, both methods deliver reasonable grid-independent solutions up to a grid resolution of 2 cm in the case of Geldart type A particles.


Introduction
Fluidized beds are widely used in a variety of industrially important processes. During the last decades, the analysis of the hydrodynamics or the efficiency of fluidized beds through numerical simulations has become increasingly common [1][2][3][4][5][6][7][8][9][10], where the two-fluid model (TFM) approach has proven to provide fairly good predictions of the hydrodynamics of gassolid flows [1].
However, due to computational limitations, a highly resolved simulation of industrial-scale reactors is still unfeasible [11,12]. Therefore, it is straightforward to apply coarse grids, which allow a considerable reduction of the computational effort. Nevertheless, such a procedure inevitably neglects the small (unresolved) scales [13,14]. Consequently, many sub-grid drag modifications have been proposed to account for the effect of these small mesoscale structures on the resolved macroscales [2,12,[15][16][17][18][19][20][21].
In this context, filtered two-fluid models (fTFMs) [12,[15][16][17][18][19]22] as well as heterogeneity-based two-fluid models [17,20,21,23] have been developed. While the heterogeneitybased sub-grid models, such as the energy minimization multiscale (EMMS) model [23,24], make some specific assumptions about the form of the unresolved clusters and streamers, the filtered approach aims to find an appropriate constitutive relation for the unresolved terms appearing due to filtering of the balance equations. Closure models are commonly deduced from highly resolved simulations, either Euler-Euler or Euler-Lagrange [25,26], which are filtered using filters of different sizes. Various markers, such as solids volume fraction and slip velocity, are then employed to classify the sub-filter scale state and averaged to obtain statistics of the filtered quantities. Commonly, based on these markers, functional fits are provided that the fTFM can be applied in coarse-grid simulations [12,[15][16][17][18][19]22].
In previous studies [27][28][29], a different approach for deriving constitutive relations for the unresolved terms appearing in fTFM has been presented. On the one hand, a spatially averaged two-fluid model (SA-TFM) was advanced, which is based on the concepts of turbulence modeling [30]. For example, the residual correlation for the filtered drag force can be expressed by the covariance of the solid-phase volume fraction and gasphase velocity [31]. Thus, the filtered drag force can be computed from the variance of the solid volume fraction and the turbulent kinetic energy of the gas phase. Closure models for these quantities have been derived from first principles [27,32].
On the other hand, an approximate deconvolution method (ADM) [33,34] for the large eddy simulation (LES) of turbulent gas-solid suspensions has been presented recently [31]. With such an approach, an approximation of the unfiltered solution is obtained by repeated filtering, allowing the determination of the unclosed terms of the filtered equations directly.
The main advantage of such so-called structural models is their purely mathematical nature without requiring physical assumptions; thus, these are established with no a-priori knowledge of the physics of the interactions between the resolved and sub-grid scales (SGSs) [35].
In this paper, the SA-TFM and approximate deconvolution method-two-fluid model (ADM-TFM) are discussed and applied to two different fluidized beds. First, these coarse-grid models are employed for a lab-scale fluidized bed of Geldart type B particles. The results are analyzed with respect to timeaveraged solids volume fraction and Reynolds stress and subsequently compared with fine grid reference data [36]. Second, the predictions of both turbulence models are applied to a pilot-scale fluidized bed of Geldart type A particles [37] and opposed to experimental data. The latter fluidized bed is operated in bubbling and turbulent regimes.

Filtered Two-Fluid Models
In this study, a standard TFM formulation is employed to discuss the application of spatially averaged models to moderately dense gas-particle flows. It is started from the continuity and momentum equations for the gas and solid phases [1,6]: where f 1) denotes the solids volume fraction, r q is the density of the gas phase (q ≡ g) and the solid phase (q ≡ s), and v i and u i are the ith component of the gas-and solid-phase velocities, respectively. Furthermore, p is the gas-phase pressure, Σ g,ik denotes the viscous gas-phase stress, S kc ik is the mesoscale kinetic theory-based particle-phase stress, S fr ik is the frictional stress, b means the mesoscale drag coefficient, and g i is the standard acceleration due to gravity.
Constitutive models employed for the stress tensors can be found in Tabs. 3 and 4 in the previous study [2]. The drag coefficient b can be expressed by the particle relaxation time t p , i.e., b = fr s /t p , where t p is given by the Wen and Yu model [38].
In the case of coarse grids, not all relevant flow features are resolved by Eqs. (1)-(4). Thus, balance equations for the macroscale flow can be found by applying a spatial filter operation to these equations, which gives their filtered (discretized or resolved) counterpart. The filtered complement of a continuous space-time variable g(x,t) is given by: where G denotes the filter operator defined by the convolution of g(x,t) with a weighting function G(x,y,Δ fi ) satisfying ∫GdV = 1. Thus, g x; t ð Þ is determined by: where Δ fi denotes the filter width and G G? ð Þ. By applying Eq. (6) to the continuity and momentum equations, the fTFM equations [27-29, 31, 36] are obtained: where Favre averages were employed: The Reynolds stress-like contribution stemming from the convective terms reads: d D g;ik denotes the rate-of-strain tensor of the gas phase evaluated from the filtered gas-phase velocity.
For the derivation of Eqs. (8) and (10) it has been assumed that the unresolved part of the filter gas-phase pressure is negligible [18,27]. Furthermore, the contribution of the filtered mesoscale molecular stress, S kc ik , may be neglected in many  moderately dense gas-solid flows, i.e., f < 0.5, at sufficiently large filter sizes. A budget analysis [27] shows that for the grid sizes employed here the molecular stress is insignificant compared to the Reynolds stress-like contribution. Nevertheless, the contribution from the frictional contacts S fr ik becomes significant close to the maximum packing and is, therefore, retained in Eq. (10).
It has to be noted that sub-grid effects of the frictional stress are not included in this study, i.e., S fr ik c S fr ik , where c ð Þ indicates that S fr ik is evaluated from filtered quantities. Models for S fr ik will be addressed in future publications. Following previous studies [31,36], a good approximation of the filtered drag force is as follows: where the drift velocity u d,i = v i sv i g is introduced, which requires additional modeling. Furthermore, b t p is evaluated from resolved quantities.
To sum up, to finally close the fTFM Eqs. (7)-(10), constitutive relations for the Reynolds stress terms (Eqs. (12) and (13)) as well as for the drift velocity (Eq. (14)) have to be introduced. These are discussed in the following.

Functional Closures (SA-TFM)
A closure model for the drift velocity can be found from rewriting its definition as follows [32,39]: Greek symbols, i.e., κ, indicate that there is no summation over identical indices. Following the previous studies [27,31,32], the covariance fv k À fv k can be approximated by [31,32]: with k g,κ being the anisotropic components of the turbulent kinetic energy of the gas phase k g,κ = 0.5( v κ v κ gv κ g v κ g ) and x gf,κ being the correlation coefficient between gas-phase velocity v κ and the solids volume fraction j. Furthermore, f 2 denotes the sub-grid variance of the solids volume fraction. It has to be emphasized that Eq. (16) accounts for the anisotropic nature of the drift velocity due to k g,κ and x gf,κ , which is especially important in fluidized beds [22,39]. k g,κ can be directly obtained from Eq. (56) in [27] by disregarding the summation over the second index appearing in the shear production term. Finally, it has been shown recently [27] that the correlation coefficients x gf,κ are approximately This, in turn, implies that a non-zero drift velocity yields a drag reduction.
The Reynolds stress contributions are given by [27]: (18) where k q,i accounts for the anisotropic macroscale normal stress. Nevertheless, here an isotropic model for the macroscale viscosity m q,t is used since the anisotropy of the off-diagonal components of R q,ik appears less important than for the diagonal components [31], i.e., k q,i . It has to be emphasized that the constitutive relations for k g,i , k s,i , m q,t , and f 2 have already been discussed in detail previously and, thus, they are not considered here but can be found in previous works [27,28,32].

Approximate Deconvolution Method (ADM-TFM)
In contrast to the functional closure models, which are based on both mathematical derivation and physical assumptions, ADM allows a purely mathematical reconstruction of the unresolved contributions. Since the filtered momentum equations (Eqs. (8) and (10)) were obtained by applying a low-pass filter (Eq. (6)) to the macroscale balance equations (Eqs. (2) and (4)), the unresolved terms may be reconstructed by an inverse filter operation: and, therefore, the filtered drag force and the Reynolds stress contribution may be computed directly from the reconstructed f, v i , and u i . However, in most interesting cases, G is not invertible [34]. Higher-order reconstruction of G À1 can be achieved by the iterative deconvolution method of van Cittert [40], where the unfiltered quantities can be derived by a series of successive filtering operations (G) applied to the filtered quantities with [31,36]: In Eq. (20), I is the identity tensor and Q v is the approximation of the inverse filter G À1 . The truncation order of the expansion v determines the level of deconvolution [34]. Numerically, it is more efficient to rewrite Eq. (20) as an iterative expression [36]: For example, the closure model for the drift velocity (Eq. (14)) consequently reads [36]: Approximations for the Reynolds stress contributions can be found similarly [36]. In this paper, v = 2 is used, which is a good compromise between accuracy and numerical efficiency [36]. Finally, the rheology-based regularization proposed in previous work [36] is applied to all ADM simulations. Regula-rization accounts for the sub-grid scales of the Reynolds stress term, which are not recoverable by deconvolution.

Numerical Simulations
To evaluate the applicability of the sub-grid closure models presented in Sects. 2.1 and 2.2 at different scales (lab scale and pilot scale), the gas-solid fluidized beds of Geldart A and B particles were investigated. Similar to previous studies [17,36,41], a simple case with a superficial vertical gas velocity, W g in = 0.63 m s -1 , was studied (Fig. 1a). The pilot-scale bubbling fluidized bed consisted of a cylinder with an inner diameter of 0.267 m (Fig. 1b). This reactor was simulated for comparison with detailed experimental data reported by Zhu et al. [37]. The height of the reactor was 2.464 m with an added freeboard region expanding to a height of approximately 4.2 m. The freeboard had an inner diameter of 0.667 m to stop excessive particle entrainment out of the bed. The freeboard region was included in the simulation domain to accurately account for the large degree of bed expansion observed in some of the simulations conducted. Gas was injected through a velocity inlet on the bottom face of the reactor.
Different velocities were employed, namely, W g in = 0.4 m s -1 and 0.9 m s -1 . Gas exited at the top of the reactor through a pressure outlet at 0 Pa gauge pressure. In both cases, a no-slip boundary condition was applied for the gas phase and a free slip boundary condition for the solid phase at the side walls. This is contrast to the previous study [29], where no-slip boundary conditions were applied for both phases. Here, a free slip condition for the solid phase inhibits the too high production of k s near the wall due to high shear. The physical parameters are given in Tab. 1.
The open-source CFD code OpenFOAM 6 [42] was employed for numerical solution of the governing equations. Particularly, the OpenFOAM solver twoPhaseEulerFoam was modified to account for Eqs. (7)- (10). Time advancement is achieved by a variable time-step procedure, where the time step is limited by a maximum Courant number of 0.25. Pressure-velocity coupling was based on the PIMPLE algorithm [43]. Similar to [11], the SuperBee flux-limiter was used for all variable extrapolation. These code modifications can be downloaded at [44].

Results and Discussion
A previous study [27] suggests that the grid resolution for kinetic theory-based TFM should be in the order of the characteristic length scale, L ch = (u t 2 /g)Fr p -2/3 , where Fr p = u t 2 /(gd s ) is the particle-based Froude number. In particular, this length scale is a good estimate for the size of the smallest clusters, i.e., where the energy of the clusters dissipates into ''molecular'' fluctuations. This grid size requirement has also been confirmed recently by other authors [11,45]. For the group B particles used for the lab-scale fluidized bed, the characteristic length scale is L ch = 1.3 mm = 8.55 d s , which is approximately by a factor of 6-12 smaller than the grid spacings of the coarse grid. For the group A particles employed in the pilot-scale fluidized bed, the characteristic length scale is 265 mm, which is by factor 40-80 smaller than the coarse grid spacing. A previous study [17] indicates that for coarse grid spacings of about 8L ch , the contribution stemming from the interparticle collisions is insignificant compared to the ''turbulent'' Reynolds stress. These profiles are further compared with fine grid reference simulation. For more details about the fine grid simulations, the reader is referred to [2,36].

Lab-Scale Fluidized Bed
It is observed that applying either functional closures (SA-TFM) or ADM to the unresolved terms in the filtered momentum equations (Eqs. (8) and (10)) yields good agreement with the fine grid reference simulation even though much coarser grids were employed. Furthermore, both closure ap-   proaches do not unveil a considerable dependence on the grid resolution for the grid sizes employed here, i.e., Δ f coarse = 6L ch to 12L ch . However, it has to be highlighted that in the case of Δ f coarse = 6L ch , neglecting the kinetic theory stresses may not be appropriate, which is indicated by the slight deviation of the profile of the time-averaged solids volume fraction in the case of SA-TFM. Fig. 2b depicts the vertical profile of the variance of the solids volume fraction. In the case of the fine grid the variance is computed as follows, where £ t indicates time-averaging. However, in the case of the coarse grid simulation, a second unresolved contribution becomes effective, which can be either determined from the constitutive model for f 2 in the case of SA-TFM or by deconvolution in the case of ADM-TFM [36]. Here, the subscripts r and ur denote the resolved and unresolved contributions to f 2 t . Consequently, the total variance is given by the sum of resolved and unresolved contributions. Fig. 2b clearly reveals that both closure approaches are also able to recover the variance of the solids volume fraction appropriately, which is important for mixing as well as heat and mass transfer. The profiles of f 2 give evidence of a pronounced bubbling and clustering within the fluidized bed even though using a coarse grid, which is also confirmed by comparing snapshots of the spatial distribution of the solids volume fraction (Fig. 3).
Figs. 2 and 3 demonstrate that in the case of coarse grids the standard TFM considerably overpredicts the expansion of the bed, i.e., TFM underestimates the mean solids concentration with the fluidized bed. This is because TFM neglects the contribution from the unresolved heterogeneous structures, which is also generally agreed in literature [2,16,20,21,46]. Furthermore, since no constitutive relation accounts for f 2 , f 2 t,ur is zero and this, in turn, implies that f 2 t solely consists of the resolved part f 2 t,r . Thus, the variance of the solids volume fraction is greatly undervalued as well. The snapshots in Figs. 3i and 3j further support this observation showing much smoother structures (i.e., less corrugated) than depicted in Figs. 3a-3h.
Recent studies [31,47,48] demonstrated that in fluidized gas-solid flows the macroscale solid stress, i.e., the Reynolds stress-like contribution, R s,ik , can be highly anisotropic. Cloete et al. [47] suggested that employing isotropic closures for the normal components of R s,ik , i.e., for k s,i , may lead to an overestimation of the bed expansion due to the overprediction of the horizontal stress (R s,11 ≡ k s,1 and R s,22 ≡ k s,2 ), which pushes the solids away from the bounding walls. However, this theory has not been thoroughly proven in coarse-grid simulations yet [47] and, therefore, this will be addressed in future publications. Fig. 4 demonstrates that both SA-TFM and ADM-TFM appropriately predict macroscale stress in vertical direction. Furthermore, the results indicated that the horizontal macroscale stress is much smaller than the vertical component (not shown here). Nevertheless, the SA-TFM approach slightly underpredicts R s,33 in the upper part of the bed, which might be attributed to the assumption of constant correlation coefficients [27]. However, in future work, it is intended to employ a dynamic procedure following Germano and Lilly [49,50], which may considerably improve the predictions of the SA-TFM approach due to locally adjusting the correlation coefficient x jg,κ (Eq. (16)) by test filtering [32].
Finally, the posteriori simulations using ADM-TFM have to be further discussed. In particular, the filtering operator (Eqs. (5) and (6)) divides the flow into so-called resolved and sub-filter scale (SFS) motions [51]. When Eqs. (7)-(10) are solved on a discrete grid, a discretization operator is applied to the equations as well, which further divides the turbulent flow field; these are the resolved SFS (RSFS) and unresolved SFS regions, where the latter is commonly called sub-grid scale. However, ADM is solely able to recover the RSFS contribution, while the SFS region is technically lost [31]. Thus, the contribution from the SGS region has to be modeled; this is commonly called regularization [51].
To account for the energy dissipation within the SGS region, an additional dissipative term on the right-hand side of Eqs. (8) and (10) has to be added. A previous study [36] clearly shows that employing an appropriate rheology model yields such a dissipative contribution. Nevertheless, in the case of much coarser grids than used for the lab-scale problem, the closure of the drift velocity (Eq. (22)), which appears in the closure for the filtered drag term, requires additional regularization as well. This will be discussed in the next section. In contrast, functional models such as the SA-TFM approach do not distinguish RSFS and SGS contributions. Constitutive relations are supposed to cover both contributions appropriately.

Pilot-Scale Fluidized Bed
As discussed above, the grid spacing used for the pilot-scale case is a factor of about 40-80 coarser than the grid size required for standard TFM. At such coarse grids, the SGS contribution gets non-negligible and, thus, the drag term necessitates regularization. A closure for the drag regularization can be found from Eq. (15) by using a gradient hypothesis [52], that is: where a constitutive relation for the ''turbulent'' viscosity m g,t is required. Following previous work [27], m g,t may be approximated by: where C mg 0.4. Eq. (24) is quite similar to the Smagorinsky model used for single-phase LES, but with an additional term accounting for the production of turbulent kinetic energy due to the interfacial work [32,53].
In Fig. 5, the time-averaged axial pressure gradient for the two different superficial gas velocities is depicted. Here, the time-averaging was employed over 20 s of simulation time. This time-averaging interval window was large enough to obtain time-independent mean values. The figure clearly reveals that employing kinetic theory-based (standard) TFM considerably underestimates the axial pressure gradient as the gas-solid drag force does not account for the unresolved heterogeneous structures [29,54]. In contrast, the SA-TFM model appropriately predicts the pressure gradient for both superficial gas velocities even though the grid spacing is nearly two orders of magnitude larger than the grid resolution required for TFM. These results are also quite similar to the results obtained using an earlier implementation of the SA-TFM approach based on the commercial CFD solver FLUENT [29].
In the present study, all results are achieved based on an OpenFOAM implementation. Furthermore, the ADM-TFM ap-  proach yields an even better estimate of the time-averaged pressure gradient. However, additional regularization for the filtered drag term is required. Neglecting the SGS contribution of the drag force (Eq. (23)) reveals a much smaller pressure gradient (dotted line in Fig. 5a). This, in turn, implies that at such coarse grids the SGS contribution is much larger than the RSFS part. Fig. 6 presents the time-averaged radial solids volume fraction at z = 0.6 m. The figure reveals that in the case of bubbling, i.e., W g in = 0.4 m s -1 , the coarse grid models correctly predict the radial profile of the solids volume fraction. This is observed for both grid resolutions. Furthermore, the degree of the segregation of the solid phase in the vessel is in fairly good agreement with the experimental data. This is in contrast to the previous study [27], where an isotropic simplification of the SA-TFM approach was applied; in this case, the solid segregation is considerably overestimated.
It has to be further noted that the particular experimental case with which the simulations are compared exhibited a large degree of non-symmetry in the flow, even after 30 s of averaging, due to a spiraling bubble motion in the bed [54]. This can be seen from the three different time-averaged radial measurements (R1, R2, and R3), where one of them differed significantly from the others. This asymmetry shown in Fig. 6 could not be accurately captured by the presented coarse grid approaches. Since the spiraling bubble motion causing this asymmetry is a fairly isolated phenomenon [54], a precise simulation match is therefore not pursued.
In the case of W g in = 0.9 m s -1 , the radial profiles are symmetric and both coarse grid approaches are in good agreement with the experimental data (Fig. 6). As already discussed along with Fig. 5, the standard TFM considerably underestimates the solids holdup for both superficial gas velocities.
Finally, Fig. 7 presents snapshots of the solids volume fraction. The figure clearly unveils that the SA-TFM as well as the ADM-TFM approach appear to be insensitive to the grid resolution with respect to the pressure gradient and solids holdup in the dense region of the bed. However, in the case of ADM-TFM, minor discrepancies between 1-and 2-cm grids can be observed in transition to the freeboard, which is much more distinct for the 1-cm grid. However, due to grid-coarsening, detailed features of the bubbles are inevitably lost, which have to be accounted for by the present closure modeling in the framework of the fTFM. The figure further indicates that even though both methods yield quite similar time-averaged properties (Figs. 5 and 6), the form of the bubbles becomes different. Especially in the case of Δ f coarse = 1 cm, much smaller bubbles are received from the ADM-TFM method, while in the case of SA-TFM the fluidized bed shows considerable slugging. This difference may be attributed to the underlying nature of both approaches. The idea of ADM is to reconstruct sub-filter information from the resolved field, whereas in the case of SA-TFM the sub-filter models include different physical assumptions. For example, constant correlation coefficients are used to close sub-filter variances, e.g., Eq. (16). Applying a dynamic adaption procedure may considerably improve the local predictions of the SA-TFM approach [32].
Finally, it has to be noted that for the present cases, ADM-TFM and SA-TFM required about 20 % more computational resources compared to standard kinetic theory-based TFM due to the larger number of model equations at the same grid resolution. However, since ADM-TFM and SA-TFM allow much coarser grid resolutions, their computational demand is at least three orders of magnitude lower than corresponding fine grid simulations.

Conclusions
In this article, the previously presented SA-TFM [27,32] and ADM-TFM [31,36] approaches were applied to different fluidized-bed regimes, such as bubbling and turbulent, utilizing group A and B particles. The SA-TFM approach is based on the turbulent behavior of the heterogeneous gas-solid struc-tures. These appear as additional sources for the turbulent kinetic energies of both phases and for the sub-filter variance of the solids volume fraction due to the interfacial work. Those quantities, in turn, characterize the sub-grid heterogeneity and, therefore, the unresolved terms in the filtered balance equations can be determined. In the ADM-TFM method, the unresolved terms are closed by an ADM. With such an approach, an approximation of the unfiltered solution is obtained by repeated filtering allowing the determination of the RSFS contribution of unclosed terms of the filtered equations directly.
The main findings are: Neglecting the contribution of the unresolved terms fails to predict the hydrodynamics of fluidized beds using a coarse mesh.
Both approaches, SA-TFM and ADM-TFM, are able to predict the hydrodynamics of bubbling and turbulent fluidized beds correction even in the case of very coarse meshes.
Furthermore, both approaches appear to be nearly insensitive to the grid spacing up to 80 times the grid resolution required for standard kinetic theory-based TFM.
However, in the case of the ADM-TFM method, additional regularization of the drag term is required for very coarse meshes, where the SGS contribution gets significant compared to the RSFS part.
To conclude, this study indicates that the ADM-TFM approach is more accurate than the SA-TFM method when relatively fine grids can be afforded computationally. This can be explained by the lower degree of modeling involved in ADM. However, the accuracy of ADM-TFM decreases with increasing grid resolution since ADM by itself does not account for the SGS contribution. Thus, additional regularization methods for the Reynolds stress and the unresolved part of the drag force are required, which necessitates additional modeling, e.g., Eq. (25). Thus, the results of this study reveal that at coarser grid resolutions the application of SA-TFM is preferable with respect to accuracy and computational costs, while in the case of ADM-TFM additional regularization decreases the numerical performance of the method.
However, several tasks remain. First, these methods have to be further tested in the fast fluidization regime. Especially, the regularization method for the filtered drag force should be evaluated in more detail. Second, the performance and accuracy of these approaches have to be studied using non-orthogonal and non-structured meshes as well. Finally, both approaches have to be augmented by heat and mass transfer, which is inevitable for their industrial application. These issues will be addressed in future studies.

Acknowledgment
The financial support by the Austrian Federal Ministry for Digital and Economic Affairs and the National Foundation for Research, Technology and Development as well as from the K1MET Center for Metallurgical Research in Austria (www.k1-met.com) is gratefully acknowledged.
The author has declared no conflict of interest.