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

Abstract

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, defined 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 configuration. Computational fluid dynamics simulations were then performed using an Eulerian-Granular multiphase model. The effects 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.

1. 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, three-dimensional, 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–6] or Eulerian multiphase modeling approaches [7–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–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.

2. 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 (𝑇) 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 (𝐷) 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 (𝐶₁) of 0.15 m. For the dual impeller tests, a second impeller was installed on the same shaft at an off-bottom distance (𝐶₂) of 0.39 m. The liquid height (𝐻) was kept constant at 1.5𝑇 in all experiments.

Figure 1 

Schematic representation of the stirred tank reactor.

The solid phase was comprised of sand particles with a nominal average diameter of 180 𝜇m and a density of 2,630 kg/m³. The liquid phase was water with a density of 1,000 kg/m³ 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 (𝐻_𝑐) could then be calculated for each experimental condition. Measurements were carried out at impeller speeds (𝑁) ranging from 150 to 450 RPM in 50 RPM increments.

For both impeller systems, the minimum impeller speed for just suspension (𝑁_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. 𝑁_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 𝑁_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 (𝐻_𝑐/𝐻) versus impeller speed normalized by 𝑁_JS. At impeller speeds below 𝑁_JS (𝑁/𝑁_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.

Figure 2 

Fractional average cloud height (𝐻𝑐/𝐻) as a function of impeller speed normalized by 𝑁JS for the single- and dual-impeller configurations for 10% (w/w) and 15% (w/w) solids.

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

3. Modeling and Numerical Simulations

3.1. Euler-Granular Theory

The simulations were performed using the Euler-Granular multiphase model available in the commercial CFD software, ANSYS Fluent 13.0 (ANSYS Inc., Canonsburg, PA, USA). This model assumes that each phase coexists at every point in space in the form of 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 𝑞 is [21]: 𝜕𝜕𝑡𝛼^𝑞𝜌_𝑞+∇⋅𝛼_𝑞𝜌_𝑞→𝑢_𝑞𝜕=0(1)𝜕𝑡𝛼^𝑞𝜌_𝑞→𝑢_𝑞+∇⋅𝛼^𝑞𝜌_𝑞→𝑢_𝑞→𝑢_𝑞=−𝛼^𝑞∇𝑝+∇⋅𝜏_𝑞+^𝑛_𝑞=1𝑅^𝑓𝑠+𝑚̇_→𝑓𝑠𝑢_𝑓𝑠+^→𝐹_𝑞,(2) where terms with the subscript “𝑓𝑠” 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 𝑓 and 𝑠, 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 (𝜏_𝑞) 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), 𝐹^→_𝑞 accounts for the effects of lift force, virtual mass force and any other external forces acting on the solid phase. The term 𝑅_𝑓𝑠 represents fluid solid interphase momentum exchange which is expressed as: 𝑅_𝑓𝑠=𝐾_𝑓𝑠^→𝑢_𝑓−^→𝑢_𝑠,(3) where 𝐾_𝑓𝑠 is the interphase exchange coefficient, which is defined in terms of the drag model proposed by Ding and Gidaspow [22]: 𝐾_𝑓𝑠=34𝐶^𝐷𝛼_𝑠𝛼_𝑓𝜌_𝑓|||^→𝑢_𝑓−^→𝑢_𝑠|||𝑑_𝑠𝐾𝛼_𝑓−2.65for𝛼_𝑓>0.8,_𝑓𝑠=150𝛼^𝑠1−𝛼_𝑓𝜇_𝑓𝛼_𝑓𝑑_𝑠2+1.75𝛼^𝑠𝜌_𝑓|||^→𝑢_𝑓−^→𝑢_𝑠|||𝑑_𝑠𝐶for𝛼_𝑓≤0.8_𝐷=24Re𝛼_𝑓1+0.15Re𝛼^𝑓^0.687,Re=𝜌^𝑓|||^→𝑢_𝑓−^→𝑢_𝑠|||𝑑_𝑠𝜇_𝑓.(4) For turbulent flows, the drag force (𝐶_𝐷) 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: 𝑅_𝑓𝑠=𝐾_𝑓𝑠^→𝑈_𝑓−𝑈^→_𝑠−𝐾^→𝑓𝑠𝑣_dr.(5) Here 𝑈^→ represents phasic average velocities and ^→𝑣_dr is the drift velocity:^→𝑣_dr=−𝐷^𝑓𝜎_𝑓𝑠𝛼_𝑓∇𝛼_𝑓−𝐷^𝑠𝜎_𝑓𝑠𝛼_𝑠∇𝛼_𝑠.(6) In this expression, 𝐷_𝑓=𝐷_𝑠=𝐷_{𝑡,𝑓𝑠} are the fluid-solid turbulent dispersion coefficients, and 𝜎_𝑓𝑠 is the dispersion Prandtl number (which is set equal to 0.75 in this study).

3.2. Simulation Setup

3.2.1. Geometry and Mesh

The geometric models and meshes were constructed in the commercial software, GAMBIT 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.

Figure 3 

CFD meshes for the single- and dual-impeller stirred tanks.

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.

Table 1 

Grid dependency study for the single-impeller tank. (A) 10% (w/w) solids loading and 300 RPM, (B) 15% (w/w) solids loading and 250 RPM.

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

4. Results

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

(a)

(b)

Figure 4 

Contour plots for 200 RPM, 10% (w/w) loading for single- and dual-impeller configurations on a midplane. (a) Velocity vector of liquid phase. (b) Volume fraction of solid phase.

(a)

(b)

Figure 5 

Contour plots for 400 RPM, 10% (w/w) loading for single- and dual-impeller configurations on a midplane. (a) Velocity vector of liquid phase. (b) Volume fraction of solid phase.

4.2. 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 𝑧-coordinate of an isosurface of solid phase volume fraction (note that the shaft axis is aligned with 𝑧-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.

Figure 6 

Effect of drag model. Volume rendered plots with isosurface for single- and dual-impeller configurations, 15% (w/w) loading, 300 RPM. The double-headed arrows marked “EXP” represent the location of experimental observation with arrow thickness representing the observed 8% spread in the experimental data.

Figure 7 

Effect of drag model. Plot of normalized cloud height against normalized impeller speed for 15% (w/w) loading and single-impeller.

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

Figure 8 

Effect of turbulent dispersion force (TDF). Volume rendered plots with isosurface for (a) single-impeller, 10% (w/w) loading, 400 RPM and (b) dual-impeller, 15% (w/w) loading, 150 RPM. The double-headed arrows marked “EXP” represent the location of experimental observation with arrow thickness representing the observed 8% spread in the experimental data. Gidaspow drag model.

4.4. 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 single- and 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. This may be due to the fact that the single impeller may be not able to fully suspend solids at low RPM. Modifications to the TDF model to include effect of local solid concentration may overcome this difficulty. This analysis is left for future study.

(a)

(b)

Figure 9 

(a) Plot of normalized cloud height against normalized impeller speed for 10% (w/w) loading; (b) plot of normalized cloud height against normalized impeller speed for 15% (w/w) loading.

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

Nomenclature

𝐶1:	Distance of lower impeller from the tank bottom, m
𝐶2:	Distance of upper impeller from the tank bottom, m
𝐶𝐷:	Drag coefficient
𝐷:	Impeller diameter, m
𝑑𝑠:	Particle diameter, m
𝐷𝑡,𝑓𝑠:	fluid-solid turbulent dispersion coefficient,  m2/s
𝐹→𝑞:	Lift, virtual mass, and any other external forces acting on the solid phase,  N/m3
𝐻:	Liquid fill level, m
𝐻𝑐:	Cloud height, m
𝐾𝑓𝑠:	Interphase exchange coefficient,  kg/m3-s
𝑚̇𝑓𝑠:	Fluid-solid interphase mass exchange,  kg/m3-s
𝑁:	Impeller rotational speed, rpm
𝑁JS:	Just suspension speed, rpm
𝑝:	Pressure shared by all phases,  N/m2
Re:	Relative Reynolds number
𝑅𝑓𝑠:	Fluid-solid interphase momentum exchange,  N/m3
𝑇:	Tank diameter, m
→𝑢𝑓:	Velocity of fluid phase, m/s
→𝑢𝑓𝑠:	Interphase velocity, m/s
→𝑢𝑞:	Velocity of 𝑞th phase, m/s
→𝑢𝑠:	Velocity of solid phase, m/s
𝑈→𝑓:	Phasic average velocity of fluid phase, m/s
𝑈→𝑠:	Phasic average velocity of solid phase, m/s
→𝑣dr:	Drift velocity, m/s.

Greek letters

𝛼𝑞:	Volume fraction of phase 𝑞
𝜏𝑞:	Stress-strain tensor for phase 𝑞,  N/m2
𝜇𝑓:	Kinematic viscosity of liquid phase, kg/m-s
𝜎𝑓𝑠:	Dispersion Prandtl number
𝜌𝑞:	Density of phase 𝑞,  kg/m3.

Subscripts

𝑓:	Fluid phase
𝑓𝑠:	Exchange between fluid and solid phases
𝑞:	Phase number
𝑠:	Solid phase
𝑡:	Turbulent.

Acknowledgments

The authors gratefully acknowledge the invaluable help provided by Prem Andrade and Mohan Srinivasa (ANSYS, Pune) for their guidance during the CFD simulations.