CFD Modeling of Solid Suspension in a Stirred Tank: Effect of Drag Models and Turbulent Dispersion on Cloud Height

Many chemical engineering processes involve the suspension of solid particles in a liquid. In dense systems, agitation leads to the formation of a clear liquid layer above a solid cloud. Cloud height, deﬁned as the location of the clear liquid interface, is a critical measure of process performance. In this study, solid-liquid mixing experiments were conducted and cloud height was measured as a function operating conditions and stirred tank conﬁguration. Computational ﬂuid dynamics simulations were then performed using an Eulerian-Granular multiphase model. The e ﬀ ects of hindered and unhindered drag models and turbulent dispersion force on cloud height were investigated. A comparison of the experimental and computational data showed excellent agreement over the full range of conditions tested.


Introduction
The suspension of solid particles in a liquid is a key requirement of many industrial processes. Examples relevant to this work include crystallisation, dissolution, and adsorption processes. Each of these can be characterized as complex multiphase processes that are often facilitated by mechanically stirred tanks. Maximizing contact between the solid and liquid phases facilitates mass transfer and reaction, therefore, assessing the ability of the process to suspend particles is a key objective. The quality of suspension can be quantified by three parameters: just suspension velocity, solids distribution, and cloud height. The latter is defined as the location of the interface between solids-rich liquid and clear liquid regions. Cloud height is critical because of limited mixing between the suspended solids and the upper clear liquid layer. The current study aims to develop and validate a computational model of a solids suspension process using experimental measurements of cloud height.
The hydrodynamics in stirred vessels are complex, threedimensional, and turbulent. The interplay between solid and liquid materials, vessel design, impeller type and location, and level of agitation all determine the efficiency of the mixing process. Numerous empirical models have been developed to relate tank performance to operating conditions and geometry [1,2]. However, these models lack wider application, owing to the complexity of the flow field. This has helped to motivate the need for more detailed analysis.
Computational fluid dynamics (CFD) methods are powerful tools for improving our understanding of mixing in stirred tanks. In the last few years, encouraging results have been obtained in the simulation of solid suspension systems [3]. The solids distribution has been traditionally simulated using either Lagrangian [4][5][6] or Eulerian multiphase modeling approaches [7][8][9][10][11][12][13][14]. The Lagrangian model treats solids as individual particles subject to hydrodynamic forces that are approximated using single-particle empirical models. The Eulerian model, on the other hand, treats the solids as an additional continuous phase. Interphase transport models are added to account for the interaction between the solids and the fluid phase.
The Granular model is an extension of the Eulerian model wherein the solids are modeled as a pseudo-fluid. The pseudo-fluid physical properties, such as viscosity, solids pressure and stress, are derived from the kinetic theory of granular flow. According to Montante et al. [15], the Granular model has been found to be more appropriate for modeling processes exhibiting high local solid concentrations.
Turbulence modeling is also critical to the success of CFD analysis in solid-liquid systems. Turbulence plays an important role in developing single-phase flow structures, and it is generally accepted that it will also have some influence on the solids distribution in stirred tanks. Previous CFD investigations have highlighted the importance of turbulence-assisted solids dispersion in stirred tanks [16][17][18][19][20]. However, a quantitative evaluation of the effect of turbulent dispersion on cloud height has not been considered in detail.
To address the above, this study aims to develop a CFD modeling strategy for predicting cloud height in a solid suspension process. An Eulerian-Granular multiphase model is used to model the particle suspension process. Experimental data is used to justify the selection of model inputs, specifically the drag model and turbulent dispersion force. Each of these were evaluated independently to illustrate their relative importance. Finally, a comparison of the model results with corresponding experimental data establishes the validity of the modeling approach across a wide range of experimental conditions.

Experimental Setup and Procedure
The experimental work was carried out at BHR Group as part of the fluid mixing processes (FMPs) industrial research program. The tests were performed in a transparent Perspex tank with a diameter (T) of 0.61 m, equipped with four standard baffles and a torispherical base, as shown in Figure 1. The impellers were centrally mounted pitched blade turbines (PBTs) with a diameter (D) of 0.2 m and a blade width of 0.055 m. One set of tests was carried out with a single-impeller located at a distance from the tank bottom (C 1 ) of 0.15 m. For the dual impeller tests, a second impeller was installed on the same shaft at an off-bottom distance (C 2 ) of 0.39 m. The liquid height (H) was kept constant at 1.5T in all experiments.
The solid phase was comprised of sand particles with a nominal average diameter of 180 μm and a density of 2,630 kg/m 3 . The liquid phase was water with a density of 1,000 kg/m 3 and a viscosity of 0.001 Pa-s. Solids concentrations of 10% (w/w) and 15% (w/w) were tested. These values correspond to 4% and 6.2% by volume, respectively.
The solids cloud height was obtained by measuring the distance from the bottom of the tank to the cloud surface. Measurements were determined visually from outside the tank. Given the turbulent nature of the flow, the cloud surface was not flat. Therefore, the minimum and maximum positions of the surface at a fixed radial position on the wall of the tank were recorded. The average cloud height (H c ) could then be calculated for each experimental condition. Measurements were carried out at impeller speeds (N) ranging from 150 to 450 RPM in 50 RPM increments.
For both impeller systems, the minimum impeller speed for just suspension (N JS ) was also determined at each solids loading. This is the minimum speed at which the impeller(s) operate to ensure that the solids are suspended, that is, not settled at the vessel bottom. N JS was determined by placing a mirror underneath the tank and observing the bed of solids at the bottom while the impeller speed was increased in small increments. The criterion used to determine whether the solids were suspended followed that of Zwietering [1], who defined N JS as "the impeller speed at which no material remains stationary on the vessel base for more than 1-2 seconds." Figure 2 shows the results for the fractional average cloud height (H c /H) versus impeller speed normalized by N JS . At impeller speeds below N JS (N/N JS < 1), the cloud height is considerably higher with the dual impeller system than with the single impeller system (up to about 70% higher). The difference becomes less significant at higher impeller speeds.
With the single-impeller system, there is no significant difference between the two different solid concentrations.
There is a small difference at high impeller speeds with the dual-impeller system. However, this could be due to experimental error (typically 5 to 10%). interpenetrating continua. The continuity and momentum equations are solved for all phases and the coupling between phases is obtained through pressure and interphase exchange coefficients. The conservation equation for each phase q is [21]:

Modeling and Numerical Simulations
where terms with the subscript " f s" account for exchange between the fluid and solid phases. One fluid phase and one solid phase were modeled in this work. They are denoted individually with subscripts f and s, respectively, in the rest of this section.
In the Eulerian-Granular model, the solid phase momentum equation includes an additional solid pressure term. The solids pressure and solids stress tensor terms (τ q ) are determined from the kinetic theory of granular flow as described by Ding and Gidaspow [22]. The Euler-Granular model also includes an additional transport equation for the granular temperature.
In (2), F q accounts for the effects of lift force, virtual mass force and any other external forces acting on the solid phase.
The term R f s represents fluid solid interphase momentum exchange which is expressed as: where K f s is the interphase exchange coefficient, which is defined in terms of the drag model proposed by Ding and Gidaspow [22]: for α f > 0.8, For turbulent flows, the drag force (C D ) is a function of average velocity. A correction term is included in the interphase momentum exchange equation to account for the effect of turbulent dispersion due to instantaneous fluctuations, as proposed by Deutsch and Simonin [23]. The resulting interphase momentum exchange term is then given by: Here U represents phasic average velocities and v dr is the drift velocity: In this expression, D f = D s = D t, f s are the fluid-solid turbulent dispersion coefficients, and σ f s is the dispersion Prandtl number (which is set equal to 0.75 in this study).

Geometry and Mesh.
The geometric models and meshes were constructed in the commercial software, GAM-BIT 2.4 (ANSYS, Inc., Canonsburg, PA, USA). Full hexahedral meshes were constructed for each of the two tank configurations (Figure 3). The computational mesh for the single impeller tank had 47, 72, and 122 cells in the radial, angular, and axial directions, respectively; while the mesh for the dual impeller configuration had 47, 72, and 177 cells in the respective directions. These settings resulted in 402k hexahedral elements for the single-impeller tank and 576k hexahedral elements for the dual-impeller tank. (Note that multiplying the mesh sizings in each direction will result in a mesh count that is higher than the final mesh count quoted in the text. This is because there is no mesh at impeller and shaft locations and because there are small regions of paved mesh which are not regularly spaced.) The mesh spacing was more or less uniform in each direction, with mesh refinement on wall surfaces and in the impeller zones.   A grid-dependency study was carried out to evaluate mesh suitability for the single impeller configuration. The procedure involved comparing results from the baseline mesh to a refined grid with 3.2 million hexahedral cells, representing an almost 8-fold increase over the baseline mesh. The effects on one flow parameter-impeller power number-and one multiphase model output-cloud height-were extracted and compared (see Table 1). The Gidaspow drag model and turbulent dispersion force were included in all cases. The maximum difference in the impeller power numbers was less than 4%, and the maximum difference in the predicted cloud height was less than 2%. Therefore, the resolution of the baseline mesh was deemed satisfactory.

Physical Models and Boundary Conditions.
The realizable k-ε turbulence model was used to resolve the turbulent solid-liquid flow field. For comparison purposes, interphase momentum exchange terms were calculated using the Gidaspow (hindered) and Schiller-Naumann (unhindered) drag models. The effect of including turbulence-assisted solids dispersion in the interphase exchange term was also evaluated [23,24]. The multiple reference frame (MRF) approach was used to model impeller rotation. Although an approximation, this approach has been shown to produce satisfactory results,  especially when the interaction between baffles and impellers is weak. The free surface was assumed to be fixed, thus a symmetry boundary condition was prescribed. A no-slip boundary condition was applied for all other wall surfaces.
High-order discretization schemes were chosen to minimize numerical diffusion. The Green-Gauss Node Based Gradient option was used for gradient calculations at cell interfaces and the QUICK discretization scheme was used for momentum, volume fraction and turbulence equations. All cases were initialized with a uniformly distributed solid phase. The steady state solver was utilized to solve for all flow variables.

Flow and Solid Distributions.
Liquid velocities provide insight into the mixing processes occurring inside stirred tanks. Figures 4 and 5 show model results for liquid velocities and solid phase volume fraction on a tank midplane. Recirculation zones are evident near the tank bottom. These recirculation zones are predicted to be larger in the dual impeller configuration and increase in size with RPM. A distinct interface between clear liquid at the top of the tank and the cloudy region is observed at low RPM. The model also predicts that solids are more uniformly distributed throughout the entire reactor at higher RPM.

Effect of Drag.
Drag forces induce significant interphase momentum transfer and tend to dominate other interphase processes. Therefore, drag force model selection is critical to the accuracy of the solid-liquid suspension model. The Gidaspow and Schiller-Naumann drag models were investigated. The Gidaspow drag model has been specifically tailored to take into consideration particle induced hindrance. For this reason, it is expected to give good predictions for high solid volume fractions (as was observed under certain operating conditions studied here). Additional simulations with the Schiller-Naumann model were performed to help elucidate the importance of particle hindrance in this system. Figure 6 shows volume rendered (volume rendered postprocessing permits visualization of the solid-phase distribution throughout the entire domain by increasing the opacity of the foreground results) contour plots of solids concentration and cloud height for the single and the dual impeller configurations. The double-headed arrows marked "EXP" represent the location of the experimentally observed cloud height with the arrow thickness representing the observed ±8% variation in the experimental data. When post-processing the CFD model, the cloud height was defined as the maximum z-coordinate of an isosurface of solid phase volume fraction (note that the shaft axis is aligned with z-axis). Following Kasat et al. [12], the location of this isosurface is defined as the average solid phase volume fraction. It can be seen from Figure 6 that the Gidaspow model is a better predictor of cloud height versus the Schiller-Naumann model. Experimental results further substantiate the Gidaspow drag model as a better predictor of cloud height (see Figure 7). These results suggest that particle hindrance effects are an important aspect of the model, even at the relatively low solids loading studied here.

Effect of the Turbulent Dispersion Force (TDF).
The contribution of the turbulent dispersion force is significant when the size of the turbulent eddies is larger than the particle size [12]. Figure 8 shows how the inclusion of TDF substantially increases the cloud height. Moreover, the predicted cloud heights give better comparison with the experimental data when TDF is included. In an extreme example, it was observed that omitting TDF leads to a 250% difference in the cloud height prediction (see Figure 8(b)). Similar studies performed with the Schiller-Naumann drag model showed an identical trend, although the predicted cloud heights were observed to be slightly lower than that predicted by the Gidaspow. Therefore, TDF must be included in the model for the reactor configuration and conditions studied here.

Experimental Validation.
Results from the previous sections establish the model requirements based on a limited set of experiments. Applying the modeling strategy over the full range of experimental operating conditions establishes model validation. The experimental study included singleand dual-impeller tanks operating at six to seven speeds and for two different solids concentrations (10% (w/w) and 15% (w/w)), resulting in a total of 26 experimental measurements. Numerical simulations were performed for each of these experimental points. The profiles of the normalized cloud heights are plotted in Figures 9(a) and 9(b). It can be seen that the predicted cloud height for all operating conditions are in good agreement with the experimental data. For the single-impeller system, the simulations overpredict cloud height at lower RPM for both 10% (w/w) and 15% (w/w) solids loading.

Conclusion
A CFD modeling strategy was developed to predict cloud height in stirred tanks with single-and dual-axial PBT impellers. The CFD model was based on Eulerian-Granular multiphase theory. The effects of two different drag models, Gidaspow and Schiller-Naumann, and a drift velocity-based turbulent dispersion force were analyzed. The modeling strategy was developed for a limited set of experimental conditions and then applied across the full range solids loading, agitation rates, and reactor configurations. It was determined that the inclusion of turbulent dispersion force was critical to accurate prediction of cloud height for the conditions studied. The drag model was observed to have a much lower influence. The resulting modeling strategy was then compared to the full experimental range of cloud height measurements. The model was found to be predictive of cloud height across a broad range of experimental conditions, the exception being the single-impeller configuration at low agitation rates. It is the authors' opinion that solids concentration-based inclusion of the TDF can be studied further to improve the model prediction at lower agitation rates.