CFD simulations of stirred-tank reactors for gas-liquid and gas-liquid-solid systems using OpenFOAM^®

2.1 Governing equations

The Eulerian multi-fluid method (Drew 1983) is used to model the systems in this work. All phases are treated as interpenetrating continua and represented by their volume fractions. The governing equations are summarized below.

The continuity equations are written as

where αk, ρk and U_k are the volume fraction, density and velocity of the gas (k = g), liquid (k = l) and solid phase (k = s), respectively.

The momentum equation is written as

where p is the shared pressure, τ_k is the stress tensor of phase k, g is the acceleration due to gravity and M_jk represents the momentum exchange between phases j and k. The momentum exchange between phases can include contributions from the following forces: drag, lift, virtual mass, turbulent dispersion and wall lubrication. In this work, drag, virtual mass and turbulent-dispersion forces are included. The buoyancy force results in the shared-pressure term in Eq. (2). The turbulent-dispersion force is required to ensure that the multi-fluid model is hyperbolic when a separate bubble pressure is not employed. Due to the intense mixing provided by the impellers, other exchange forces such as lift were neglected.

2.2 Turbulence models

In the present work, large-eddy simulations (LES) are used to account for the turbulence in the reactor for each phase. The sub-grid scale (SGS) viscosity for phase k is calculated using the Smagorinsky sub-grid model (Smagorinsky 1963):

where the Smagorinsky constant C_S = 0.168 in this work. Δ is the filter width calculated as the cubic root of the volume of each cell, and D_k is the rate-of-strain tensor of phase k. The stress tensor of phase k, τ_k in Eq. (2), is expressed as

where μ_{k, lam} is the dynamic viscosity of phase k, and I is the unit tensor. The SGS turbulent kinetic energy of phase k can be found by assuming local equilibrium for k_k, expressed as

where the SGS stress tensor Bk=23ρkkkI−2μk,SGS[Dk−13tr(Dk)I] and the constant C_E = 1.048 in this work.

2.3 Multiple reference frame (MRF) model

The rotation of the impellers is treated by using the MRF model (Luo, Issa, and Gosman 1994), which divides the computational region into a rotating region and a stationary region. In the stationary region, governing equations are solved in a rotating reference frame. At the interface of the two regions, a reference frame transformation is performed.

The fluid velocities can be transformed from the stationary reference frame to the rotating reference frame by calculating the relative velocities as

where U_{r, k} is the relative velocity of phase k viewed from the rotating reference frame, U_k is the absolute velocity viewed from the stationary reference frame, ω is the angular velocity of the rotating coordinate system relative to the stationary reference frame, and r is the position vector of a point in the rotating region with respect to the origin of the rotating reference frame. When the governing equations are solved in the rotating reference frame, the continuity equation is written as

The momentum equation in the rotating reference frame is expressed as

where the extra term αkρk(ω×Uk) represents the Coriolis and centripetal accelerations.

3.1 Description of reactor 1 and operating conditions

Reactor 1, as described in Ford et al. (2008), is a cylindrical vessel with a flat bottom and four baffles. A six-blade Rushton turbine (RT) is mounted on the shaft near the bottom of the reactor. Air at 300 K is injected uniformly from a ring sparger located at the bottom of the tank into water (300 K). Details on the tank and impeller geometries are listed in Table 1.

The commercial mesh generation software Pointwise^® is used to generate 420,624 three-dimensional hexahedral and 2436 polyhedral cells with 1,308,248 faces for Reactor 1, as shown in Figure 1. The grid convergence study is performed and given in Table S1 in Supplementary Materials. The gas-liquid flow is found from the two-fluid solver reactingTwoPhaseEulerFoam in OpenFOAM^® (OpenFOAM 2015). The initial water level in the reactor is 0.21 m. Air is injected uniformly from the sparger at gas flow rate 1.5 × 10⁻⁴ m³/s (9 LPM in Ford et al. 2008). Three impeller rotation speeds are considered 200, 350 and 700 rpm. No-slip boundary conditions are used for both the gas and liquid phases at the tank wall and baffles. Constant gas inlet velocity and atmospheric pressure are used at the gas inlet and outlet, respectively. At the shaft and impeller surfaces, zero relative velocity (U_{r, k} = 0) is assumed for both phases. The simulated flow time for all cases of Reactor 1 is 65 s, and results are time averaged over the last 60 s.

Figure 1:

(a) Illustration of the computational mesh for reactor 1. (b) Enlarged view of the mesh near the impeller.

3.2 Description of reactor 2 and operating conditions

The stirred-tank reactor, studied by (Shewale and Pandit 2006) and referred to as Reactor 2, is a cylindrical vessel with a flat bottom and four baffles. Three impellers are mounted on the shaft and rotate synchronously. The bottom impeller is a six-blade RT, and the middle and top impellers are pitched-blade downflow turbines (PBTD). Details on the tank and impeller geometries are listed in Table 2. In total 669,888 three-dimensional hexahedral cells are generated for Reactor 2 using Pointwise^®, as shown in Figure 2. The grid convergence study for Reactor 2 is given in Table S2 in Supplementary Materials.

Figure 2:

(a) Illustration of the computational mesh for reactor 2. (b, c) Enlarged views of the mesh near the RT impeller (bottom) and PBTD impeller (top).

The OpenFOAM^® solver reactingTwoPhaseEulerFoam (OpenFOAM 2015) is used to solve for the gas-liquid (air-water) flow in the reactor. Initially the water level in the tank is 0.87 m. A gas-liquid-solid flow is also simulated with OpenFOAM^® solver reactingMultiphaseEulerFoam (OpenFOAM 2015). The three phases are air, water and polymethyl methacrylate (PMMA) particles, respectively. The initial mixture level is 0.81 m, with solid and liquid volume fractions 0.2 and 0.8, respectively. For both two-phase and three-phase systems, air at 300 K is injected uniformly from the sparger located below the RT into a water-PMMA mixture at 300 K with superficial velocity 0.01 m/s. A uniform bubble diameter (0.5 cm) is used for all simulations of Reactor 2. For the gas-liquid-solid simulations, the particle diameter is assumed constant at 150 µm with density 1190 kg/m³. The impellers are rotating at speeds of 3.75 and 5.08 rps in the down-pumping direction for the gas-liquid simulations, and a rotation speed of 5.08 rps for the gas-liquid-solid simulation. At the tank wall and baffles, a no-slip boundary condition is used for all phases. Constant gas inlet velocity and atmospheric pressure are used at the gas inlet and outlet, respectively. At the shaft and impeller surfaces, zero relative velocity (U_{r, k} = 0) is used for all phases. Simulated flow time for the gas-liquid system is 65 s, and results are time averaged over the last 60 s. For the gas-liquid-solid system, results are time averaged over 40 s due to the heavy computational cost.

4.1 Gas-liquid simulations of reactor 1

The overall gas holdup in the vessel αOverall can be calculated as

where H₀ is the initial water level in the tank (0.21 m). Table 3 lists the averaged mixture height and overall gas holdup in the tank at different impeller rotation speeds and for three bubble diameters of 0.5, 1.5 and 2.5 mm. The comparison of other bubble diameters (3.5 and 5 mm) with the experimental gas holdups are given in the Supplementary Materials (Table S3). The simulation results are compared to experimental data reported in Ford et al. 2008, where αCT and αGlobal are the values obtained using different techniques. Overall, the agreement between the simulations and experiments is good and indicate that a bubble diameter of approximately 0.5 mm can capture the dependence of the overall gas holdup on all rotational speed. Thus, we will look next at the predictions of the local gas holdup at different heights and radial positions in the reactor.

Figure 3 compares the experimental data for the gas holdup to the simulation results as a function of the height and radial position in the reactor for three bubble sizes. In Figure 3 (a, c and e) the gas volume fractions along a line (x-axis in Ford et al. 2008) at height 0.065 m is plotted for bubble diameters of 0.5, 1.5 and 2.5 mm. From the experimental results the gas volume fraction increases in the impeller region. Although a similar trend is observed from simulations, in general, the simulations underestimate the gas volume fraction, and only the 0.5 mm bubble diameter at 700 rpm mixing rate agrees well with experiments. On the other hand, as shown in Table 3 the average gas hold-up values are well correlated with the experimental results where αCT and αGlobal are the averaged gas holdup in the CT imaging region and global value in the tank, respectively. This agreement indicates that, since CT imaging is a time averaged technique and cannot capture the cavities behind the impellers (especially at high mixing rates due to blade rotation), relatively small discrepancies between experiments and simulations are obtained.

Figure 3:

Comparison of simulated and experimental gas volume fractions as a function of height and radial position. Black lines with circles, red lines with triangles and blue lines with rectangles represent the experimental results for 700, 350 and 200 rpm mixing rates respectively, and corresponding colored symbols represent simulation results: (a, c, e) Gas volume fraction along a line at height 0.065 m for 0.5, 1.5 and 2.5 mm bubble diameter. (b, d, f) Averaged gas volume fractions at different heights for 0.1, 1.5 and 2.5 mm bubble diameter at different heights.

These results are more obvious in Figure 3(b, d, f) which provides the comparison of gas volume fractions along the reactor in the axial direction. As given in Table 1, the distance of the end of blades from the bottom is 0.076 m and up to this height the experimental and simulation results are different yet in closeness, proximity. Away from the impeller blades, the simulation results are almost the same as the experimental measurements. In comparison of Figure 3(a, b) at high mixing rate (700 rpm – dispersed regime) although the gas holdup is estimated well in the impeller region with 0.5 mm bubble diameter, better agreement is obtained for higher bubble diameters, 1.5 and 2.5 mm, outside of this region. Therefore, it can be concluded that the bubble size distribution is most likely not uniform in the reactor and break-up of bubbles happens within the impeller region due to applied shear and coalescence is the dominant mechanism above the blades. Furthermore, in the loaded regime (350 rpm) good agreement is obtained with 1.5 mm bubble diameter (Figure 3(d)) since a more uniform velocity profile (hence shear applied on the bubbles) can be expected in this regime. This conclusion is supported by the gas and liquid phase velocity magnitude profiles shown in Figure 4 (and in Supplementary MaterialsFigure S4).

Figure 4:

Gas velocity magnitude profiles for bubble diameter d = 1.5 mm at different mixing rates and different heights (0.065, 0.087 and 0.11 m): (a–c) 200 rpm, (d–f) 350 rpm, (g–i) 700 rpm.

In Figure 4(a–c) contour plots of the gas-phase velocity magnitudes for the bubble diameter of 1.5 mm are given at impeller rotation speed 200 rpm at three heights 0.065, 0.087 and 0.11 m, corresponding to z = 0.8, 3.0 and 5.3 cm in Figure 5 of (Ford et al. 2008), respectively. (The gas-phase velocity profiles for d = 0.5 and 2.5 mm are given in Supplementary MaterialsFigures S1 and S2). Comparing the three different bubble diameters, although the velocity profiles are similar, the velocity magnitudes increase proportionally with the bubble size. These results clearly indicate that the impeller is flooded at this rotation speed and gas flow rate. The bubbles concentrate near the impeller and shaft and are barely dispersed away from the impeller region (Paglianti 2002) and bubbles spread radially due to the decreased pressure and bubble rise velocity (Ford et al. 2008). Compared to (Ford et al. 2008), similar patterns are observed, i.e., gas holdup is high near the blades, and bubbles spread radially with increasing vertical distance. The gas volume fraction is very low near the tank wall. (In the Supplementary Materials, the flow pattern of water at impeller rotation speed 200 rpm is shown in Figure S4(a–c) at the same heights with the gas phase velocities. The flow pattern of water for 0.5 and 2.5 mm bubble diameters are also provided in Figures S3 and S5). For all bubble diameters no radial-pumping flow pattern can be observed because the impeller is flooded.

Figure 5:

Contour plots of gas volume fraction at 3.75 rps for gas-liquid system:

(a) Axial cross section. (b) Radial cross section at bottom RT (0.15 m). (c) Radial cross section at middle PBTD (0.45 m). (d) Radial cross section at top PBTD (0.75 m).

Figure 4(d, e, f) gives contour plots of the gas-phase velocity magnitude profiles at impeller rotation speed 350 rpm. Compared to Figure 6 in Ford et al. 2008, similar patterns can be observed, and the impeller is loaded. The gas volume fraction is still high near the impeller (see Figure S9(e) in Supplementary Materials). However, a high gas volume fraction region can be found near the impeller tips, which indicates bubbles are beginning to be dispersed. As seen in Figure 6 of (Ford et al. 2008), a ring pattern can be observed in Figure 4, which means that gas is going radially and around the impeller. As is typical of the loaded regime, a well-dispersed gas phase can be expected, while the gas volume fraction is low below the impeller. Similar to air flow pattern, radial-pumping flow patterns is observed in the liquid phase (see Figure S4(d–f) in Supplementary Materials) as well.

Figure 6:

Flow patterns of gas phase at 3.75 rps for gas-liquid system: (a) 2D velocity vectors on the axial cross section. (b) Contour pot of magnitude of velocity on the axial cross section. (c) Contour plot of axial velocity on the axial cross section. (d–f) Contour plots of magnitude of velocity on the radial cross section at bottom RT (d), middle PBTD (e), and top PBTD (f). (g–i) contour plots of axial velocity on the radial cross section at bottom RT (g), middle PBTD (h), and top PBTD (i).

Contour plots of the gas-phase velocity magnitude profiles at impeller rotation speed 700 rpm are shown in Figure 4(g, h, i). Similar patterns to Figure 7 in (Ford et al. 2008) can be observed, and the flow is in the completely dispersed regime so that the gas is well dispersed throughout the reactor. (Figure S4(g, h, i) in the Supplementary Materials shows the flow pattern of water at impeller rotation speed 700 rpm.) Like the air velocity profile at this rotation speed, a radial-pumping flow pattern for liquid phase is obtained.

Figure 7:

Contour plots of gas volume fraction at 5.08 rps for gas-liquid system:

(a) Axial cross section. (b) Radial cross section at bottom RT (0.15 m). (c) Radial cross section at middle PBTD (0.45 m). (d) Radial cross section at top PBTD (0.75 m).

4.2 Gas-liquid simulations of reactor 2

Table 4 lists the gassed impeller power consumption P_G in the tank at two different impeller rotation speeds (3.75 and 5.08 rps) found from the simulations. The comparison of fractional gas holdup, ε_G, at the corresponding P_G both for experiments and simulations are also given in Table 4.

For the experimental fractional gas holdup values, in Shewale and Pandit (2006) a power-law type correlation as a function of (PG/V) and V_G is used, where P_G, V and V_G are the gassed impeller power consumption, bulk volume and superficial gas velocity, respectively, with fitting constants given in (2). Our simulation results overpredict the experimental fractional gas holdups (see Table 4) with approximately 30% percent error. The difference can be attributed to different power consumptions, calculated in experiments and in simulations. In Shewale and Pandit (2006), power consumption is calculated according to

where P_o and Q_G are the ungassed impeller power consumption and gas flow rate. Although a similar correlation is proposed in (Michael and Miller 1962) for single disc turbine, α and n parameters are added in Shewale and Pandit (2006) as fitting parameters for PBTD and are given in Table 2 of (Shewale and Pandit 2006).

On the other hand, in the simulations the gassed impeller power consumption in the vessel (P_G) is calculated by multiplying the magnitude of the torque (T) by the angular velocity ω:

where N_R is the impeller rotational speed. The torque is calculated by integrating the cross product of the total stress (σ) on the surfaces of the impellers and the position vector (r), written as

where n is the normal vector of cell face on the surface of shaft and impellers, and total stress (σ) is defined as

with stress tensors τ_g and τ_l being calculated using Eq. (10). It is unclear whether the results from Eq. (17) can be compared quantitatively to the correlation in Eq. (16) which is based on the total power consumption.

Figure 5 gives contour plots of the gas volume fraction in Reactor 2 at rotation speed 3.75 rps. Results indicate that the gas is dispersed effectively at the middle and top PBTDs, but not so well below and above the bottom RT, which agrees with the picture of Effective dispersion-Flooding (DE-F) regime shown in Shewale and Pandit (2006). A ring pattern can be found at the bottom RT, which means gas is going radially and around the impeller, but the impeller cannot disperse gas very effectively to the region away from the impeller. The flow pattern of air at impeller rotation speed 3.75 rps is shown in Figure 6. Air is mainly going upward from the sparger to the top of the tank. At the bottom RT, air is going up around the impeller. Upward and downward circulations are barely formed below and above the impeller. At the middle and top PBTDs, downward circulations can only be observed very close to the impellers.

Radial-pumping and down-pumping flow patterns of water are found at the bottom RT, and air is going up around the impeller (see Figure S10 in Supplementary Materials). Upward and downward circulations are barely formed below and above the impeller. At the middle and top PBTDs, downward circulations can only be observed very close to the impellers. Radial-pumping and down-pumping flow patterns of water at the bottom RT and upper PBTDs at 3.75 rps are more obvious than the flow patterns of air, shown in Figure 5. At the bottom RT, an upward and a downward circulation can be seen clearly below and above the impeller, and water is moving downward near the tank wall. At the middle PBTD, a downward circulation is formed, but still close to impeller region. An upward circulation is also formed between the bottom RT and middle PBTD. At the top PBTD, a downward circulation if formed, and expanded to regions away from the impeller.

Figure 7 shows contour plots of gas volume fraction at rotation speed 5.08 rps. Gas is dispersed effectively at the upper PBTDs, and dispersed better above the bottom RT, compared to gas dispersion at rotation speed 3.75 rps. At the bottom RT, a high gas volume fraction region can also be observed near the tank wall and baffles. Simulation results are similar to the photo shown in Shewale and Pandit (2006), described as the Effective-Loading (DE-L) regime.

The air flow pattern at impeller rotation speed 5.08 rps shown in Figure 8. Like the flow pattern of air at 3.75 rps, air is mainly going upward. A radial-pumping flow pattern can be observed at the bottom RT, although the upward and downward circulations are very subtle below and above the impeller. At the middle and top PBTDs, downward circulations are formed close to the impellers (Figure S11 in the Supplementary Materials shows the water flow pattern at rotation speed of 5.08 rps). A radial-pumping flow pattern is clear at the bottom RT. Below the bottom RT, an upward circulation is formed, while above the impeller, a downward circulation can be observed. At the upper PBTDs, downward circulations that are expanded to regions away from the impellers are generated.

Figure 8:

Flow patterns of gas phase at 5.08 rps for gas-liquid system: (a) 2D velocity vectors on the axial cross section. (b) Contour pot of magnitude of velocity on the axial cross section. (c) Contour plot of axial velocity on the axial cross section. (d–f) Contour plots of magnitude of velocity on the radial cross section at bottom RT (d), middle PBTD (e), and top PBTD (f). (g–i) Contour plots of axial velocity on the radial cross section at bottom RT (g), middle PBTD (h), and top PBTD (i).

4.3 Gas-liquid-solid simulations of reactor 2

Contour plots of gas, solid and liquid volume fractions at impeller rotation speed 5.08 rps are shown in Figure 9. The addition of solid particles changes the volume fraction profile significantly. Air is sucked into the mixture near the free surface. At the bottom RT, the ring pattern of gas volume fraction disappears, and the impeller tends to become flooded. A high gas volume fraction is observed near the shaft and impellers, while outside the impeller regions gas is not dispersed effectively. Accordingly, solid and water volume fractions are low near the shaft and impellers. For the solid phase, the particle concentration increases slightly with decreasing distance from the bottom of the vessel.

Figure 9:

Contour plots of gas (left), solid (middle) and liquid (right) volume fractions at 5.08 rps: (A–C) axial cross section. (D–F) radial cross section at bottom RT. (G–I) radial cross section at middle PBTD. (J–L) radial cross section at top PBTD. (A, D, G, J) gas phase. (B, E, H, K) solid phase. (C, F, I, L) liquid phase.

With the addition of a solid phase, the radial-pumping flow pattern at the bottom RT seen in the gas-liquid system disappears, and downward circulations at the upper PBTDs can barely be observed. Air mostly goes upward throughout the vessel, and only goes downward in a small region between the blades of the upper PBTDs. Compared to the gas-liquid system, the water flow pattern changes dramatically with the addition of solid particles (See Figure S12 in Supplementary Materials). The radial-pumping flow pattern at the bottom RT cannot be observed anymore. Instead, an upward circulation is formed near the impeller. The upper PBTDs present a radial-pumping more than a down-pumping flow pattern. Upward and downward circulations are formed below and above both impellers. Furthermore, solid particles follow the liquid, and therefore the solid phase has a similar flow pattern to the liquid phase.