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

An open-source CFD software OpenFOAM® is used to simulate two multiphase stirred-tank reactors relevant to industrial processes such as slurry polymerization and fuel production. Gas-liquid simulations are first performed in a single-impeller stirred-tank reactor, studied experimentally by Ford, J. J., T. J. Heindel, T. C. Jensen, and J. B. Drake. 2008. “X-Ray Computed Tomography of a Gas-Sparged Stirred-Tank Reactor.” Chemical Engineering Science 63: 2075–85. Three impeller rotation speeds (200, 350 and 700 rpm) with three different bubble diameters (0.5, 1.5 and 2.5 mm) are investigated. Flow patterns compared qualitatively to those from experiments. Compared to the experimental data, the simulations are in relatively good agreement for gas holdup in the reactor. The second multiphase system is a multiimpeller stirred-tank reactor, studied experimentally by Shewale, S. D., and A. B. Pandit. 2006. “Studies in Multiple Impeller Agitated Gas-Liquid Contractors.” Chemical Engineering Science 61: 486–504. Gasliquid simulations are performed at two impeller rotation speeds (3.75 and 5.08 RPS). The simulated flow patterns agree with published pictures from the experiments. Gas-liquid-solid simulations of the multiimpeller stirred-tank reactor are also carried out at impeller rotation speed 5.08 RPS. The addition of solid particles with a volume fraction characteristic of slurry reactors changes the flow pattern significantly. The bottom Rushton turbine becomes flooded, while the upper pitched-blade downflow turbines present a radial-pumping flow pattern instead of down-pumping. Nonetheless, the solid phase has a similar flow pattern to the liquid phase, indicating that the particles modify the effective density of the fluid.


Introduction
Stirred-tank reactors (STRs) are widely used in industrial processes, such as in the chemical, biological, pharmaceutical, as well as many other industries. Usually, multiphase systems are of most interest in these processes, including twophase (gas-liquid/gas-solid/solid-liquid/liquid-liquid) systems and three-phase (gas-liquid-solid) systems (Achouri et al. 2012;You et al. 2014). Gas-liquid STRs with a single impeller have been studied extensively using experiments. Investigations of the effects of reactor parameters and operating conditions such as the geometry of the tank, impeller type, size and location, impeller rotation speed and gas flow rate on the reactor characteristics including flow regimes and patterns, gas hold-up, power-consumption, mixing time and mass-transfer coefficient between the phases have been carried out for decades (Ford et al. 2008;Joshi et al. 2011a;Kumaresan and Joshi 2006;Miller 1974;Nienow et al. 1985;Nienow 1998;Parasu-Veera, Patwardhan, and Joshi 2001;Rewatkar et al. 1993;Sawant and Joshi 1979). For large-scale, tall-thin, industrial gas-liquid STRs, multiple impellers are usually preferred over a single impeller because multiple impellers can have better gas utilization due to the high residence time of gas bubbles, maintain homogeneity in the reactor, provide high surface-to-volume ratio and offer lower shear than single-impeller STRs for shear-sensitive systems (Ahmed et al. 2010;Dutta and Pangarkar 1995;Himmelsbach et al. 2006;Kasat and Pandit 2004;Khopkar et al. 2006;Shewale and Pandit 2006;You et al. 2014). However, the complexity of the flow in STRs with multiple impellers increases with increasing number of impellers as any change in the reactor parameters, such as tank aspect ratio, number, type, size, location and configuration of impellers and operating conditions, may influence the reactor operating characteristics significantly. Therefore, many researchers have studied multiple-impeller STRs experimentally (Kasat and Pandit 2004;Shewale and Pandit 2006;Vr'abel et al., 1999Vr'abel et al., , 2000You et al. 2014) for gas-liquid two-phase flows. On the other hand, experiments of gas-liquid-solid STRs with single or multiple impellers are reported less in the literature than gas-liquid STRs (Boyer, Duquennt, and Wild 2002;Chen, Bao, and Gao 2009;Conway, Kyle, and Rielly 2002;Dohi et al., 2002Dohi et al., , 2004Dutta and Pangarkar 1995;Fishwick, Winterbottom, and Stitt 2003;Joshi et al. 2011a;Micale et al. 2000;Nienow and Bujalski 2002;Rewatkar and Joshi 1992), as the addition of a solid phase increases the flow complexity in the reactor and greatly limits the range of optical-based measurement techniques that can be applied.
Computational fluid dynamics (CFD) simulation is another approach to study mixing and fluid dynamics in STRs. While, many single-phase CFD simulations of STRs with single or multiple impellers has been reported in literature, for gas-liquid systems the main focus is devoted to single impeller STRs (Achouri et al. 2012;Arlov, Revstedt, and Fuchs 2008;Bakker and Van Den Akker 1994;Gosman et al. 1992;Jahoda, Tomaskove, and Moste 2009;Jahoda et al. 2007;Joshi et al. 2011a;K'alal, Jahoda, and Fort 2014;Lane, Schwarz, and Evans 2005;Luo, Issa, and Gosman 1994;Moilanen et al. 2008;Murthy and Joshi 2008;Petitti et al. 2013;Ranade and Van Den Akker 1994;Scargiali et al. 2007;Tyagi et al. 2007) and in these studies, the effects of reactor parameters and operating conditions, and different models such as turbulence models, drag correlations and bubble coalescence and breakage models have been studied. As mentioned above, due to complicated flow regimes and patterns in gasliquid STRs with multiple impellers, not as many works as for single-impeller STRs are reported in the literature (Ahmed et al. 2010;Joshi et al. 2011b;Kerdouss, Bannaari, and Proulx 2006;Khopkar and Tanguy 2008;Khopkar et al. 2006;Min et al. 2008;Tyagi et al. 2007). Even more, for gas-liquid-solid systems, very few studies are reported for single-impeller STRs (Murthy, Ghadge, and Joshi 2007;Panneerselvam, Savithri, and Surender 2009) while for multi-impeller STRs, works found in the literature are very limited and almost all of these studies are performed using commercial CFD software, such as ANSYS FLUENT ® and CFX ® .
In this work, an open-source CFD software, Open-FOAM ® , is used to numerically investigate two stirred-tank reactors studied experimentally (Ford et al. 2008;Shewale and Pandit 2006), but also of interest for industrial applications such as slurry reactors. Gas-liquid simulations are first performed in a single-impeller STR, reported in (Ford et al. 2008). Then a multi-impeller STR, studied experimentally by (Shewale and Pandit 2006), is simulated for both gas-liquid and gas-liquid-solid systems. However, it should be noted, in their study, Shewale and Pandit (2006) only studied gas-liquid two phase system. On the other hand, in this study, in order to observe the effect of the addition of solid phase to the gas and liquid phase individually, three-phase (gas-liquid-solid) simulations are also presented. In the remainder of this work, Section 2 describes modeling of the multiphase system. The equations used to model the system are outlined. Section 3 introduces the two simulated STRs, and the simulation conditions are described in detail. Section 4 discusses the simulation results, where Section 4.1 presents results obtained in the single-impeller STR and comparison to experiments, and Section 4.2 presents gas-liquid simulation results of the multi-impeller STR. Gas-liquid-solid simulations are discussed in Section 4.3. Conclusions are drawn in Section 5. Results for parameters outside the range reported in the main text are available in the Supplementary Materials.
2 Modeling of the multiphase system

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

Interphase drag correlation
The drag force between the gas and liquid, or the solid and liquid phases is expressed as where the dispersed phase can be the gas (k = g) or solid phase (k = s), and the coefficient K kl depends on the specific drag correlation. The drag correlation of (Tomiyama et al. 1998) is used for the drag force between the gas and liquid phases: where d g is the diameter of gas bubbles, and the drag coefficient C D is The bubble Reynolds number Re and Eötvös number Eo are defined as with μ l and σ being viscosity and surface tension, respectively. The Wen and Yu drag correlation (Wen and Hu 1966) is used for fluid-solid interactions: where

Re
(1 + 0.15Re)0.687 Re < 1000 and the slip Reynolds number Re is with d s being the representative particle diameter (taken to be constant).

Virtual mass force
A constant virtual-mass coefficient (C VM = 0.5) (Drew and Lahey 1987) is used to account for the gas-liquid virtual mass force, defined by where D Dt represents the material derivative based on each velocity field.

Turbulent-dispersion force
The gas-liquid turbulent-dispersion force can be expressed as where a constant turbulent-dispersion coefficient (C TD = 1.0) (Burns et al. 2004) is used, and k l is the turbulent kinetic energy of the liquid. This force is required to keep the multi-fluid model hyperbolic when the bubble pressure is neglected in the shared-pressure model. In summary, in Eq. (2), for the gas phase (k = g), and M ls = −M sl = −F D ls for the solid phase.

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 B k = 2 3 ρ k k k I − 2μ k, SGS [D k − 1 3 tr(D k )I] and the constant C E = 1.048 in this work.

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 (ω × U k ) represents the Coriolis and centripetal accelerations.

Simulated reactors and operating conditions
Two stirred-tank reactors are studied in the present work. The first reactor was investigated experimentally by (Ford et al. 2008) and is referred to as Reactor 1. The second reactor, referred to as Reactor 2, was investigated experimentally by (Shewale and Pandit 2006). Details about these reactors and the operating conditions are described next.

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 twofluid 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 −4 m 3 /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.

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. 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 twophase 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 3 . 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.

Results and discussion
In (Ford et al. 2008), X-ray computed tomography (CT) images of the local gas holdup obtained at different planes are presented. They also report local gas holdup along a line at a specific height and average gas holdup at different heights. Section 4.1 discusses simulation results of this work and compares simulation results to the experimental data. In Shewale and Pandit (2006), experimental results obtained for a gas-liquid flow are discussed, including flow patterns for different regimes and the relationship between the overall gas holdup and power consumption. Section 4.2 presents the gas-liquid simulation results of Reactor 2 and comparisons to experimental data. For the gas-liquid-solid system, no experiment results are reported in the literature, but this system is of interest for slurry polymerization and other industrial processes. Section 4.3 therefore discusses simulation results only.

Gas-liquid simulations of reactor 1
The overall gas holdup in the vessel α Overall can be calculated as where H 0 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.
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 rpmdispersed 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 Materials Figure S4).  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 gasphase velocity profiles for d = 0.5 and 2.5 mm are given in Supplementary Materials Figures 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)  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 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 welldispersed 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.
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 impeller rotation speed 700 rpm.) Like the air velocity profile at this rotation speed, a radial-pumping flow pattern for liquid phase is obtained. 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 (P G /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

Gas-liquid simulations of reactor 2
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 downpumping 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,   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.

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.
With the addition of a solid phase, the radialpumping flow pattern at the bottom RT seen in the gasliquid 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 radialpumping 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.

Conclusions
Multiphase CFD simulations were performed for two STRs. In Reactor-1, gas-liquid simulations were performed. Similar flow patterns with experiments in Ford et al. (2008) are found and the simulation gas holdup values are in good agreement with experimental overall gas holdup values. In Reactor-2, the flow patterns from gas-liquid simulations agree well with experiments (Shewale and Pandit 2006), however, calculated gassed power consumption is smaller than the experimental value. When solid particles added to the Reactor-2 simulations, the bottom RT tends to be flooded, while a radial-pumping flow pattern is observed at the upper PBTDs instead of downpumping flow.