## 1. Introduction

Slurry bubble column reactors (SBCRs) are multitubular columns, operating with gas-liquid-solid systems mostly under isothermal conditions as the heat of reaction is removed through carefully designed cooling tubes or pipes with a large surface area for heat transfer [1, 2]. The power needed for mixing and suspending the solid (catalyst) in SBCRs is primarily provided by the gas flowing upward at high superficial velocity within the churn-turbulent flow regime. The gas is often sparged from the bottom of the reactor through specially designed distributors, with controlled pressure drop and nozzle sizes.

SBCRs are used in numerous important industrial applications, such as: (1) catalytic hydrogenation of carbon monoxide to syncrude via Fischer-Tropsch (F-T) synthesis: The Oryx gas-to-liquid (GTL) process (34,000 bpd) built in Qatar by Sasol [3]; (2) hydrogenation of glucose to sugar alcohols (D-glucose, sorbitol, D-xylose, mannitol, etc.): 60,000 metric tons annual capacity plant built in China by Global Bio-Chem [4]; (3) slurry phase hydrocrackers [5]; (4) catalytic wet air oxidation treatment of wastewater [6]; (5) catalytic synthesis of organic chemicals and polyolefins [7]; (6) polymerization of ethylene in a slurry of cyclohexane and a solid catalyst (Chromium, Ziegler Natta), the Solvay process [8]; and (7) the ALFOL process for synthesis of fatty alcohols developed by Conoco [7, 9]. These reactors are often preferred over multitubular fixed-bed reactors due to their numerous advantages, including better temperature control, easier construction, ability of online catalyst addition and withdrawal, higher effectiveness factor, higher gas and liquid holdups, lower pressure drop and more reasonable interphase mass transfer rates with low energy input [1]. SBCRs, however, inherit some disadvantages, such as high degree of liquid back-mixing, strong catalyst attrition, increased side reactions, problematic catalyst separation from the molten products, complex hydrodynamics and ultimately difficult scale up.

In SBCRs, the reactants in the gas-phase have to travel from the gas-phase bulk through the gas-film, gas/liquid interface, liquid-film, liquid-phase bulk, liquid/solid film to reach the solid (catalyst) surface. The reactants then have to diffuse inside the catalyst pores until they reach its active sites, where the chemical reaction takes place. This series of events is illustrated in **Figure 1**. Once the reaction products are formed, they start diffusing all the way back from the active sites to the catalyst to the gas-phase bulk.

The overall reaction rate, Eq. (1), is shown in terms of the resistances encountered in the reaction process:

In this equation, *C _{G}* and

*C*are the concentrations of the reactant in the bulk gas and inside the catalyst particle, respectively;

_{S}*k*and

_{L}, k_{G}*k*are the liquid-side, gas-side and liquid-solid mass transfer coefficients, respectively;

_{S}*He*is the Henry’s Law constant,

*a*and

*a*are the gas-liquid and liquid-particle interfacial areas;

_{S}*k*is the pseudo-kinetic rate and

_{o}*η*is the catalyst effectiveness.

Since the diameter of the catalyst particles frequently used in SBCRs is micron-sized, the interfacial area between the liquid and the catalyst surface is considerably large, and as such, the mass transfer resistance between the liquid and solid phases could be neglected. Moreover, since the vapor pressure of the products (wax) in the gas-phase is negligible, the resistance to mass transfer in the gas-film could be neglected. Thus, the two remaining resistances which have to be considered are the reaction kinetic and the volumetric liquid-side mass transfer. SBCRs typically operate in the churn-turbulent flow regime to guarantee complete suspension of the fine catalyst particles at high solid loading. The gas holdup in SBCRs could reach 50% by volume and consequently the mass transfer behavior is primarily controlled by the gas-liquid interfacial area (*a*).

Modeling of SBCRs is a complex task which requires, among others, detailed description of the hydrodynamics, reaction kinetics and mass as well as heat transfer parameters. Over the past century, extensive efforts have been made to accurately model the behavior and performance of SBCRs with the aim of enhancing the understanding of the complex hydrodynamics and their effects on the reactor performance in order to enable better design, operation and troubleshooting [10]. Earlier studies focused mainly on experimentally examining the macroscopic fluid dynamic behavior of three-phase fluidized beds and developing empirical correlation. Similarly, empirical one-dimensional (1-D) models have been proposed for SBCRs [11, 12, 13], such as the axial dispersion models and multiple cell circulation models, which provide valuable information and predictions of the overall reactor performance. Two-dimensional (2-D) models were also proposed to account for dispersion coefficients, which could be calculated from the first principles, but were empirically obtained. The flow structure and internal recirculation zones in SBCRs, however, were ignored in both 1-D and 2-D models.

With the advent of increasing computational power, the use of computational fluid dynamics (CFD) has gained considerable attention in modeling purposes. Over the past decade, significant advances have been made in numerical modeling of two-phase, gas-solid and gas-liquid flow systems. However, understanding the behavior of the three-phase flows is still limited because of the complex phenomena underlying interactions among the phases, including the particle-bubble interaction and the liquid interstitial effect during particle-particle collision. Recently, several CFD models have been reported to simulate three-phase systems, such as fluidized beds and SBCRs [14, 15, 16, 17]. This chapter focuses mainly on the use of CFD modeling in the design of SBCRs.

## 2. CFD modeling of SBCRs

The presence of three (gas-liquid-solid) coexisting phases in SBCRs, exhibiting significant relative motions, where the exchange of mass, energy and momentum across their interfaces are a dynamic process, represent a challenging modeling task. In this situation, the two most commonly used CFD modeling strategies [18] are the Eulerian-Eulerian or multi-fluid Eulerian [19] and the Eulerian-Lagrangian frameworks [20]. The Eulerian-Eulerian approach assumes that the dispersed and continuous phases are separate interpenetrating quasi-continuous phases, which interact with each other within the computational domain in a fixed Eulerian frame, and a set of Navier-Stokes equations is solved for each of the phases. In this approach, the properties of fluid flow, momentum and heat transfer are represented at a macroscopic level rather than dealing with the individual constituents at a microscopic level. The coupling between the motion of the dispersed and continuous phases is achieved by implementing momentum exchange terms into the momentum balance equations of the respective phase. The Eulerian-Lagrangian approach, on the other hand, assumes that the dispersed phase consists of representative particles transported with the continuous phase and a set of Navier-Stokes equations, which include coupling between the continuous and dispersed phases, is solved only for the continuous phase, whereas the dispersed-phase particles are tracked by solving individual motion equations for each particle. It should be emphasized that when using CFD for modeling large-scale applications, it is critical to strike a balance between the model accuracy and usability and the computational time.

The preceding discussion indicates that the multi-fluid Eulerian modeling is the best suited for large-scale applications with multi-phase turbulent flows, such as SBCRs. Also, the Eulerian-Lagrangian modeling is best suited for fundamental studies, such as bubble-bubble and/or bubble-liquid interactions, and thus it is limited to small superficial gas velocities and phase holdups, unlike those in SBCRs. Therefore, the bulk of this chapter is dedicated to the multi-fluid Eulerian approach for modeling industrial-scale SBCRs.

Earlier, CFD works relied heavily on model simplifications to lower the computational cost, such as using a 2-D Cartesian geometry with axisymmetry and isothermal flow or lumping the gas and liquid phases into one pseudo-homogeneous phase. Recent modeling efforts [2, 17] have been more focused on accurately capturing the interphase interactions of the three phases in detailed 3-D geometries. Basha et al. [10] have comprehensively summarized CFD modeling efforts of three-phase reactors. They pointed out that one of the major drawbacks in these CFD modeling efforts of SBCRs is the lack of model validation. It is important to emphasize that careful and rigorous validation of CFD models is a critical step in the development process, which cannot be overlooked if the model is to be used for the design and optimization purposes. Moreover, validations need to be conducted keeping in mind that data obtained with air-water-solid systems in small-scale reactors at ambient pressures and temperatures cannot satisfactorily be used in the validation of SBCRs, which are designed for high-pressure, high-temperature and high superficial gas velocity applications. This is because such validation will bring into question the scaling up protocol to large-scale industrial SBCRs.

### 2.1. Multi-fluid Eulerian modeling approach

A schematic of the CFD model components based on a multi-fluid Eulerian approach for modeling SBCRs is shown in **Figure 2**. It consists of three main components: (i) the core hydrodynamics model, consisting of the Navier-Stokes equations; (ii) the multiple-fluid model based on an analog of Boussinesq approximation to represent the natural convection and (iii) the population balance equations to describe the size distribution of the dispersed gas-phase.

The mass and momentum conservation equations are:

where k indicates the phase (G for gas, L for liquid and S for solid);

The stress tensor

For the liquid-phase, the effective viscosity, accounting for the three contributions (1) molecular viscosity

The turbulent viscosity is expressed as:

For the gas-phase, the effective viscosity is represented as the sum of the molecular viscosity

where the turbulence-induced viscosity for the gas-phase is represented as:

#### 2.1.1. Turbulence model

As the mass and momentum balances are obtained through the ensemble averaging formulations, the terms

Generally, the turbulent length scale (eddy sizes) is a physical quantity related to the size of large eddies, containing the energy in the turbulent flow. A variety of length scales exists within the turbulent flow, where the largest eddies in the flow account for most of the momentum and energy transport, and their sizes are only constrained by the physical boundaries of the flow. The sizes of the smallest eddies, on the other hand, are determined by the viscosity. Therefore, the effect of viscosity increases with decreasing length scales. The smallest length scales are those where the kinetic energy is dissipated into heat. The turbulence eddies are visualized as molecules, constantly colliding, exchanging momentum and obeying laws similar to the kinetic theory of gases. Most models of Reynolds stress, using an eddy viscosity hypothesis, based on an analogy between the molecular and turbulent motions, are described as follows:

where *k* is the turbulent kinetic energy.

Although there are numerous turbulence models as outlined elsewhere [21, 22], the most widely used is the two-equation (*k-ε*) model in which the turbulent kinetic energy (*k*) and the turbulent energy dissipation (ε) are calculated based on the following governing equations, respectively:

where *k* and *ε*, respectively derived analytically using the Random-Number Generator (RNG) theory. In addition,

where

It is important to note that the selection of the most appropriate turbulence model for SBCRs will depend on the reactor geometry and the superficial velocities of the gas and liquid phases. Also, different models have to be tested and validated to determine which is the most suitable for the given SBCR design and operating conditions [18, 24].

#### 2.1.2. Momentum exchange terms and interphase coefficient correlations

In a SBCR, the solid particles are fluidized in the liquid-phase due to the momentum transfer to it from the gas-phase. Therefore, capturing accurately the interaction dynamics between the three phases is critical in developing the CFD model. The momentum exchange term in the momentum balance Eq. (3), which describes the interface forces between the phases is described as:

The right-hand side (RHS) terms of Eq. (15) represent the interphase drag, lift, virtual mass, wall force, and turbulence dispersion, respectively. The drag force is by far the most important force in describing the momentum transfer between phases, whereas other forces, such as lift, buoyancy virtual mass and turbulent dispersion have more a tuning effect, which can be optimized during the validation process.

Significant efforts have been reported investigating each of these momentum terms separately, and numerous terms have been summarized elsewhere [10]. While there is extensive work on two-phase flow systems, studies investigating three-phase hydrodynamics are rather limited. Moreover, the gas-liquid and liquid-solid drag models available were obtained mostly using air-water or air-water-glass beads under ambient conditions, which will bring into question about their applicability for modeling large-scale industrial SBCRs, operating with organic media under high pressures and temperatures. Several investigators [18, 25] have underscored the need for better and more representative drag correlations for CFD modeling of the hydrodynamics in three-phase reactors. Therefore, the selection of the appropriate terms and coefficients not only highly dependent on the system conditions, but also require extensive validation to determine the best combination that most accurately capture the behavior of the SBCR. Basha et al. [24] have recently identified the momentum interphase terms (**Table 1**) as well and the interphase coefficient expressions (**Table 2**) to represent SBCRs operating in the churn-turbulent flow regime for the Fischer-Tropsch synthesis process.

Term | Expression |
---|---|

Drag [26] | |

Lift [27] | |

Virtual Mass [28] | |

Lateral [29] | |

Turbulence dispersion [30] |

Term | Model |
---|---|

Gas-liquid drag [31] | |

Gas-solid drag [32] | |

Liquid–solid drag [33] | |

Liquid-gas lift [34] |

#### 2.1.3. Bubbles representation

At low superficial gas velocities, bubble interaction can be neglected and the bubble diameter is not significantly influenced by coalescence and breakup processes. However, when modeling SBCRs, which typically operate in the churn-turbulent flow regime, the effect of bubbles interaction becomes significant and the effects of coalescence and breakup must be taken into account. This will allow accurate representation of the range of bubble sizes, which exists within the flow, and better description of the local interfacial area concentration, which is a critical parameter in the determination of the inter-phase mass, momentum and heat transfer.

Moreover, bubble size variation is intimately related to bubble-particle collisions, since the presence of solids was reported to significantly affect the bubble sizes within the SBCR. Similarly, the presence of bubbles is critical in the suspension of solid particles, which can be transported in the wake of rising bubbles or their settling velocity can be reduced due to bubbles holding up the solid particles. Furthermore, the presence of both solid particles and gas bubbles has a significant effect on the turbulent velocity fluctuations within the liquid-phase, where solid particles and large gas bubbles greater than the eddy length scale tend to enhance the turbulent velocity fluctuations, whereas solid particles and gas bubbles smaller than the eddy length scale will dampen the velocity fluctuations.

Generally, bubble-particle collisions can yield two different consequences: the particle is ejected from the bubble surface, or the particle penetrates the bubble leading to either bubble breakage or non-breakage. Bubble-particle collisions generate perturbations on the bubble surface. After the bubble-particle collision, three factors become crucial in determining the coalescence or breakage of the bubble [26]: (1) shear stress, which depends on the liquid velocity gradient and the relative bubble-particle impact speed, and tends to break the gas bubble; (2) surface tension force, which tends to stabilize the gas bubble and drive it to recover its original shape; and (3) viscous force, which slows the growth rate of the surface perturbation, and tends to stabilize the bubble. Accordingly, the gas bubble coalescence and breakage models are classified into four main categories [35]: (1) turbulent fluctuation and collision; (2) viscous shear stress; (3) shearing-off process and (4) interfacial stability. The bubbles induced turbulence and bubble population balance are discussed below.

### 2.1.3.1. Bubble-induced turbulence model

The bubble-induced turbulence is represented by introducing two source terms, *k-ε* Eq. [36, 37] as:

where _{d} are the relative and drift velocities, respectively; and C_{ɛ3} = 1.2.

This model is rigorously derived by writing the equation of motion for a single bubble and rearranging it in terms of the fluid velocity, where the drag and mass coefficients (

The drift velocity (

where I is a 3 × 3 identity matrix

### 2.1.3.2. Bubble population balance

It is necessary to take into account the bubble breakup and coalescence phenomenon in the CFD model when a bubble column or a slurry bubble column operates in the heterogeneous flow regime. The typical approach is to use population balance models, which describe the variation in a given population property over the space and time within a velocity field. Assuming that each bubble class travels at the same mean algebraic velocity, the individual number density of a bubble class *i* can be expressed as [42]:

where _{C} and B_{B} represent the birth rate due to coalescence and breakage of bubbles; and D_{C} and D_{B} represent the death rate due to coalescence and breakage of bubbles, respectively.

There are numerous breakup and coalescence models reported in the literature, and available in commercial CFD packages. For modeling an F-T SBCR, the model proposed by Lou and Svendsen [43] for the breakup rate of bubbles was determined to be the most suitable. Their model was developed based on the assumption of bubble binary breakup (each bubble breaks up into two distinctly smaller bubbles) under isotropic turbulence. The “daughter” bubble sizes were accounted for using a dimensionless variable (

where *d _{I}* and

*d*represent the diameters of the daughter bubbles in the binary breakage of a parent bubble of diameter

_{II}*d*and volume

*v*. The breakup rate of the bubbles can be represented as:

where

On the other hand, the bubble coalescence occurring due to bubble collision is caused by three major forces: (1) wake entrainment, (2) buoyancy and (3) random turbulence. Wake entrainment has been widely accepted to be negligible [44]; and the effect of buoyancy is eliminated since all bubbles of the same class are assumed to travel at the same mean velocity. Therefore, the only remaining driving force is the random turbulence. The coalescence rate due to random turbulent collision from Prince and Blanch [45] was determined to be the most suitable for F-T SBCR modeling, as:

The contact time between two colliding bubbles with radii ri and rj is:

The time for two bubbles with radii ri and rj to coalesce is:

where ^{−4} and 10^{−8} m, respectively. rij is the equivalent radius which can be expressed as:

Then, the number density for individual bubble groups governed by birth and death due coalescence can be represented as:

(29) |

where ^{−3} s^{−1}), (Eq. (24)).

It is important to note that this bubble population balance model will be developed for a SBCR operating under typical pressures and temperatures of F-T synthesis; and it might be possible that other correlations/parameters could provide a better fit under different conditions.

### 2.1.3.3. Solid-phase representation

For the solid-phase, the effect of particle-particle interactions can be accounted for by introducing additional terms into the stress tensor. The kinetic theory of granular flow (KTGF) is used to represent the solids behavior. The stress tensor for the solid-phase is represented as:

P_{s} is the solids pressure, which represents the normal solid-phase forces due to particle-particle interaction and is expressed [47] as:

The first term in the RHS of the solid pressure equation is the kinetic contribution, which accounts for the momentum transferred through the system by particles moving across the imaginary shear layers in the flow. The second term in the RHS is the collisional contribution, which represents the momentum transferred by direct collisions.

ep is the restitution coefficient, representing the ratio of normal relative velocity after and before the collision. It equals 0.9 as proposed by Ding and Gidaspow [48].

where

### 2.1.3.4. Kinetics and mass transfer

In cases of absorption and boiling processes occurring in a large-scale flow system, significant heat and mass exchanges that occur across the interfaces separating the gas and liquid phases must be appropriately accounted for within the three-fluid model. The interphase mass-transfer rate depends on the mass-transfer coefficient, the interfacial area, concentration and the rate of chemical reaction. The mass-transfer coefficient is a function of the local hydrodynamics, which are influenced, on one hand, by the bubble shrinkage due to physical or chemical absorption and, on the other hand, by the change of the physical properties due to the heterogeneous distributions of the chemical species.

In calculating the gas-liquid mass transfer in the transient CFD simulation, it is not possible to obtain the mass transfer coefficients from the hydrodynamic data generated using multi-Eulerian simulation. This is primarily due to the limitations at the interface in the jump boundary conditions from the gas to liquid, which ideally require empirical mass transfer data or correlations to be incorporated in the CFD model.

Using CFD to model mass transfer and interfacial phenomena from the first principles is feasible at a very small scale, such a single droplet, using a Lagrangian or a volume of fluid (VOF) scheme. However, it should be noted that Gidaspow et al. [14, 48] used the granular temperature approach to derive the mass transfer coefficients in multiphase systems. Nonetheless, this approach would not be possible for a large-scale SBCR. Actually, attempting multi-fluid Eulerian to perform multiphase-multicomponent mass transfer would result in numerically unstable sources and sinks.

Therefore, the species mass transfer rate from the dispersed phase to the continuous phase per unit volume, which appeared in Eqs. (2) and (3), is defined as:

where

(38) |

Moreover, when accounting for chemical reactions, an additional species conservation equation has to be considered as follows:

where *i* in phase *k*. Also,

### 2.2. Direct numerical simulation (DNS) approach

In the case of DNS approach, the governing equations for the mass and momentum balances of each phase would be identical to Eqs. (1) and (2). However, the viscous shear stress (Eq. (3)) would not be resolved using turbulence models, but is resolved by tracking the evolution, position and flow structure near the interphase interfaces with no ensemble averaging. Typically, interface tracking or capturing algorithms are divided into three categories, front tracking, level set method and volume of fluid method. A brief outline of these methods is given below and more details can be found elsewhere [49, 50, 51, 52]:

#### 2.2.1. Front tracking method

In this method, a surface mesh within a control volume is used to track the interface front within the computational domain. The discretization of the balance equations, shown in **Figure 3**, is carried out on two separate sets of meshes to describe the solution: (1) a 3-D stage-structured non-adaptive Eulerian mesh and (2) a 2-D unstructured, adaptive and triangular Lagrangian front mesh. The Lagrangian coordinate of each marker on the interface is derived by integration from the initial position at time t = 0, as shown in Eq. (40), where the interface velocity is determined using Eq. (41).

where

#### 2.2.2. Level set method

In this method, a continuous function which is zero at the interface, positive on one side and negative on the other, is defined over the computational domain. The function (

where

where,

#### 2.2.3. Volume of fluid (VOF) method

The VOF method is an Eulerian treatment of the interface that requires accurate algorithms for the advection of the volume fraction function to preserve the mass conservation, which cannot be achieved by conventional differencing schemes due to numerical diffusion. The composition field must be advected by either the high-resolution interface capturing scheme to approximate the fluxes of volume fraction or the geometric reconstruction of the interface based on the simple line interface construction (SLIC). In this method, an indictor function, Eq. (45), is introduced to track the presence of one of the phases in the whole computational domain; and then it is used to evaluate the mixture physical properties. The advantage of this method is that it is easier to implement when compared with the other two methods and requires less computational effort. Moreover, the interface positions are not stored for each time-step, which allows for easier representation of large surface deformations and surface breakup and merging. However, the main disadvantage is that it is highly resolution dependent.

#### 2.2.4. Large Eddy simulation (LES) approach

In the large Eddy simulation (LES) approach, only large eddies are computed directly and thus a low-pass spatial filter is applied to the instantaneous Navier-Stokes conservation equations to formulate the 3-D unsteady governing equations for large-scale motions. This filtering method, also referred to as explicit filtering, allows for mesh independent results to be achieved since the truncation errors at the smallest resolved scales are not large. Different filtering methods can be found elsewhere [53]. Within the LES approach, the volume-averaging procedure can be derived by applying a spatial filter function, among others, available in the literature [54]. Also, in the LES approach, the filtering operation can be done according to Eq. (46), where

The different phases are then identified by defining a quantity reflecting the average volumetric fraction of each phase inside the computational mesh volume. Eq. (47) shows the average volumetric fraction of the *k*-th phase. A component weighted volume averaging process can further be derived [55] as shown in Eq. (48). The instantaneous property

After filtering, the mass and momentum equations, the unsteady multi-fluid flows can be expressed using Eqs. (50) and (51), respectively.

The filtered-averaged viscous shear stress component can be represented using Eq. (52).

The SGS stress tensor

The parameter

### 2.2.4.1. Basic SGS model

In this model, the unresolved stress tensor (

where

The basic model of Smagorinsky [59] assumes that the SGS eddy viscosity can be described using a length and velocity scale, as shown in Eq. (56).

where *C _{1}*)

^{0.5}. In some cases, a viscosity dampening function is used to accurately capture turbulent viscosity near solid surfaces. Details of this model are available elsewhere [54, 60, 61, 62].

### 2.2.4.2. Dynamic SGS model

The dynamic SGS model is based on the self-adaptive procedure developed by Germano et al. [63], which applies two different filters, a grid filter as described above

where L is difference between the grid and test level stress tensors represented by Eq. (59):

The test filter stress (

Consequently, the filtered and sub-test stresses are represented using Eqs. (61) and (62).

The value of *C _{d}* is a case-specific constant. The value and derivation of

*C*can be found elsewhere [64, 65, 66].

_{d}## 3. Advantages, limitations and prospects of DNS and LES

When modeling turbulence, it is important to realize that there is a wide spectrum of turbulent scales coexisting within the flow, which must be resolved. The larger eddies contain most of the turbulent kinetic energy, whereas the smaller eddies dissipate the energy that they receive from the larger eddies. Thus, the primary challenge in representing turbulence lies in accurately computing (as well as measuring) the contributions of all the scales within the spectrum. Typically, the size of the computational domain must be at least an order of magnitude larger than the scales characterizing the turbulence energy, while the computational mesh must be fine enough to resolve the smallest dynamically significant length-scale (the Kolmogorov micro-scale) for accurate simulation.

The advances of CFD modeling of turbulent multiphase flows have been primarily focused on the development and application of models to accurately capture flow behavior at higher spatial and temporal resolutions. However, in modeling of SBCRs for design and optimization purposes, the grid size should be significantly larger than the gas bubbles and solid particles. Therefore, there is a significant dependence on having accurate interphase closure relationships to ensure model validity, especially in the churn-turbulent flow regime. Nonetheless, higher resolution models, such as DNS and LES, which resolve the turbulences at significantly smaller length scales without being dependent on interphase closure relationships, present a promising prospect for future modeling of SBCRs.

In the DNS approach, the Navier-Stokes equations are numerically solved without any turbulence models or closure relationships. The conservation equations are derived by considering the microscopic instantaneous equations governing each phase. The smallest possible lengths and timescales, which are compatible with continuum formulation are considered and all the scales of motion in addition to interfacial configurations are captured or fully resolved. This approach requires model resolution on the scale of the largest as well as the smallest turbulent eddies. Also, the exact location of the interfaces separating the different phases that co-exist within the flow are determined using suitable micro-level evolutionary tracking methods (Front tracking, level-set, VOF). Therefore, the DNS approach is computationally very expensive and at present it can be applied only to low Reynolds number flows over simple geometries. In general, it has been recognized that the computational effort required for DNS simulations are proportional to

In the LES approach, the structure of the turbulent flow is viewed as the distinct transport of large- and small-scale motions. The large-scale motion is directly simulated on the scale of the underlying computational mesh; whereas the small-scale motion is represented using sub-grid scale (SGS) models. Since the large-scale motion is generally much more energetic and by far the most effective transporter of the conserved properties than the small-scale one, such an approach, which treats the large eddies as approximates of the small eddies makes perfect sense to be adopted for large-scale design and optimization scenarios. However, in comparison with the multi-fluid Eulerian approach, which has been developed to simulate large-scale macroscopic flow behavior with modeled microscopic behavior, the required computational resources for LES is comparatively very large, which currently prevents their wide scale application beyond laminar multi-phase systems. Additional details on the various numerical methods and approaches are discussed elsewhere [54, 67].

## 4. CFD modeling of a pilot-scale SBCR using the multi-fluid Eulerian approach

The multi-fluid Eulerian CFD formulation described above was used to model the local hydrodynamics (local phase holdups, phase velocity profiles and liquid turbulence intensity) in a pilot-scale SBCR (0.3-m ID, 3-m height), operating under the typical Fisher-Tropsch (F-T) process conditions. The multi-fluid Eulerian approach was implemented into ANSYS Fluent v 14.5, where the governing equations are solved using an Eulerian-multiphase segregated solver algorithm. The 3-D time-dependent simulations are conducted, both due to the nature of the geometry investigated and the bubble plume oscillations, which are characteristic of the churn-turbulent flow regime [68]. The RNG k-ε turbulence model was used, as it provides the best validation results as previously demonstrated elsewhere [69]. At the bottom of the reactor, Dirichlet velocity and volume fraction conditions for all phases were set, and a second-order spatially accurate QUICK scheme [70, 71] was employed to discretize all equations.

where β is a constant.

Moreover, a multiphase variant of the SIMPLE scheme was used for pressure-velocity coupling [72]. The first order implicit time stepping was then used to advance the solution in time. Before each simulation, mesh and time independence studies were carried out in order to optimize the solution and computational time. In all simulations, quasi-steady state numerical solutions were obtained. This means that at the end of the calculations, all variables exhibit small oscillations around steady-state values, indicating that the statistical averages were reached for all variables.

Although the inlet volumetric flow rates are known, the velocity distribution is not specified. The most widely used practice is to use the knowledge of fully developed flow in pipes to specify the inlet velocity distribution. Therefore, for laminar flow through a cylindrical inlet pipe, one can specify a parabolic velocity profile as an inlet boundary condition; however, if the feed pipes have complex shapes, which is typical in SBCRs, it is necessary to develop an additional model, which include appropriate boundary conditions.

The outlet boundary conditions were implemented taking into account the following [73, 74]: (1) F-T SBCRs are typically operated in a semi-batch mode where the liquid level may reach the top of the reactor [75] and (2) Due to the F-T reaction, the liquid (mainly hydrocarbon products) formed inside the reactor is continuously removed using appropriate filtration devices located inside the reactor. Therefore, the following outlet boundary conditions were executed to account for the two aforementioned considerations:

The real computation domain is selected to be taller than the initial height of the reactor, similar to the work by Troshko and Zdravistch [73].

The initial liquid height is set by initializing the liquid volume fraction to 1 in the zone up to a known initial liquid height, while the gas volume fraction is set to 1 in the region above that.

The ambient media in the computational zone above the reactor is set to be stagnant CO, such that the values for the pressure, backflow gas volume fraction, backflow CO gas species concentration and backflow turbulent parameters, are 1 atmosphere, 1, 1, and 0, respectively [73].

The implementation of these outlet boundary conditions allowed for modeling of the slurry-phase height expansion due to the gas holdup and heterogeneous reactions without being affected by the gas-phase backflow. Moreover, due to the strong non-linear characteristics of the model, relaxation coefficients (Patankar [67]) are introduced in the momentum conservation equations. The convergence criterion adopted from Patankar [67], based on the pressure, is given by:

where N is the mesh size

A picture of the SBCR is shown in **Figure 4**, and additional details can be found elsewhere [76]. The gas-phase used was N_{2}, the liquid-phase was an F-T reactor wax and the solid-phase used was an iron-based catalyst, and details of the physical properties of the system are available elsewhere [76]. The SBCR is provided with a spider-type gas sparger containing six arms, where each arm has 6 orifices of 0.005 m ID on each side and on the bottom, totalling 108 orifices on the gas sparger. It should be emphasized that there are no orifices at the top of the arms so that solid particles could not plug them and the generated gas jets should be able to lift any solid particles which might settle at the bottom flange of the reactor. A picture of the gas sparger described above is given in **Figure 5 (a)**. The gas sparger is screwed onto a 0.0254 m ID pipe to a height of 0.102 m from the bottom flange of the reactor. Also, its maximum height from the bottom of the reactor is about 0.152 m. The constructed geometry of this gas sparger is shown in **Figure 5 (b)**.

### 4.1. CFD prediction of the local gas holdup in the pilot-scale SBCR

The local gas-holdup is a critical parameter as it governs local pressure variations, liquid recirculation, overall mixing and heat and mass transfer within the SBCR. **Figure 6** shows the CFD model predictions as snapshots of gas holdup contours at different times of 0, 20, 40 and 60 s at an inlet superficial gas velocity = 0.20 m/s, temperature = 443 K and catalyst concentration = 11 vol%. It is clear from the figure that (1) there are significant axial and radial variations of the gas holdup within the reactor, (2) such variations are nonlinear and time-dependent and (3) the maximum gas holdup exhibited are located at the center of the reactor.

### 4.2. CFD prediction of the local liquid recirculation in the pilot-scale SBCR

**Figure 7** shows the CFD model predictions in a form of different snapshots of liquid velocity vectors in the SBCR. As can be seen in these shots smaller and faster liquid recirculation cells are present in the vicinity of the sparger from the startup until it reaches a steady-state. These smaller and faster liquid recirculations are primarily due to the geometry of the spider-type sparger used in this simulation, where all gas sparging orifices were located on the sides and bottom and none the top of each sparger arm. Moreover, only large and slow liquid recirculation cells appear to develop above the gas sparger region (around 1.2 diameter from the bottom flange). Furthermore, stronger liquid recirculations and backmixing are present near the walls, as the liquid rises upward with the gas bubbles at the center of the reactor and then flows downward near the wall. Thus, the maximum liquid velocity occurs at the center of the reactor.

### 4.3. CFD prediction of the turbulence intensity in the pilot-scale SBCR

The CFD model was also used to generate different snapshots of the local turbulence intensity contours at different times in the SBCR, as shown in **Figure 8**. The turbulence intensity is defined as the root-mean-square of the velocity fluctuations to the mean flow velocity, Eq. (65). Generally, a turbulence intensity of ≤1% is considered low and that of ≥10% is considered high [77]. As shown in **Figure 8**, higher liquid turbulence intensities are observed in the vicinity of the sparger up to 20 s, after that the turbulence intensity becomes more evenly distributed throughout the reactor, with higher turbulence intensities near the center of the reactor.

### 4.4. Sensitivity of CFD model predictions to its sub-models and parameters

When developing CFD models for the design and optimization purposes, it is critical on one hand to identify the relative importance of each of the incorporated sub-models and parameters, and on the other hand to eliminate any unnecessary terms that needlessly increase the computational time without an acceptable increase in the prediction accuracy. For SBCRs, Basha et al. [2], demonstrated that the interphase drag is the most dominant exchange term, whereas, the other terms could be eliminated without significant effect on model predications, as long as the model was carefully validated. In addition, eliminating the bubble population balance was found to significantly reduce the accuracy of the CFD model predictions and increase errors by up to 9%, despite decreasing computational time by up to 23%. Thus, the degree of complexity employed in CFD modeling is dependent on the required prediction accuracy. Obviously, increased accuracy significantly increases the required computational time.