CFD for the Design and Optimization of Slurry Bubble Column Reactors CFD for the Design and Optimization of Slurry Bubble Column Reactors

Despite the notion that computational fluid dynamics (CFD) models are considered com- plicated, expensive, time-consuming and difficult to formulate, their implementation offers an advanced prospect to move beyond empirical models, which inherit severe limitations in terms of flexibility, scale-up, and optimization of slurry bubble column reactors (SBCRs). This is because complex hydrodynamics coupled with chemical reac- tions in such reactors increase the uncertainty in using empirical models, leading to significant startup delays and overruns. Recent work by Basha et al. has shown that properly validated CFD models provide an exceptional opportunity to gain detailed temporal and spatial information about the local hydrodynamics and overall behavior as well as performance of SBCRs. This chapter provides a comprehensive overview of different CFD frameworks which could be used to model SBCRs, namely the multi-Eulerian, direct numerical simulations (DNS) and large Eddy simulation (LES). The steps required in developing CFD models and the optimization of different sub-models, such as inter- phase interactions, solid-phase representation, bubble population balance, bubble-induced turbulence, mass transfer and reaction kinetics are highlighted. Different convergence criteria for meshing, solution stability and techniques for maximizing the CFD model scale without compromising accuracy are addressed. An example of using CFD multi-Eulerian frameworks to describe the local hydrodynamics in a pilot-scale SBCR (0.3-m ID, 3-m height) operating under the Fisher-Tropsch (F-T) synthesis process are also provided.


Introduction
Slurry bubble column reactors (SBCRs) are multitubular columns, operating with gas-liquidsolid 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-toliquid (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 S are the concentrations of the reactant in the bulk gas and inside the catalyst particle, respectively; k L , k G and k S are the liquid-side, gas-side and liquid-solid mass transfer coefficients, respectively; He is the Henry's Law constant, a and a S are the gas-liquid and liquid-particle interfacial areas; k o is the pseudo-kinetic rate and η 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.

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

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); u ! k = (u, v, w) is the velocity of phase k; _ m kn is the mass transfer rate between phases k and n; u ! kn is the relative velocity between phases k and n; α is the volume fraction of each phase; M ! k is the overall interfacial momentum; p f is the fluid pressure and τ is the stress tensors.
The stress tensor τ ! k is represented as: μ eff , is the effective viscosity, which is typically the sum of the molecular and turbulent viscosities.
For the liquid-phase, the effective viscosity, accounting for the three contributions (1) molecular viscosity, μ L, M , (2) shear induced turbulence viscosity, μ L, T and (3) additional term due to the bubble induced turbulence, μ L, B is represented as:  Figure 2. Schematic of the model components and their couplings.
The turbulent viscosity is expressed as: For the gas-phase, the effective viscosity is represented as the sum of the molecular viscosity μ G, M and the turbulence-induced viscosity: where the turbulence-induced viscosity for the gas-phase is represented as:

Turbulence model
As the mass and momentum balances are obtained through the ensemble averaging formulations, the terms u k and u k 0 represent the mean and fluctuating components of the velocity, thus the unclosed terms in the stress tensor, Eq. (4), should be modeled in order to accurately account for mixing and turbulence effects. Generally, liquid mixing in SBCRs is a resultant of three major contributing mechanisms [21,22]: (1) global convective circulation of the liquidphase induced by the non-uniform radial gas holdup distribution; (2) turbulent diffusion due to the presence of large and small eddies generated by the rising gas bubbles and (3) molecular diffusion, which is negligible when compared with the other diffusion mechanisms. Thus, to predict the effects of turbulence, CFD models primarily concentrate on the methods which make use of turbulence models. These methods are typically based on scaling laws and are specifically developed to account for the effects of turbulence without recourse to a prohibitively fine mesh and direct numerical simulation (DNS). Most of the conclusions and observations regarding turbulence are based on the order-of-magnitude estimates, which follow from logical applications of scaling laws and dimensional analysis.
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 μ T is the turbulent or eddy viscosity, which unlike molecular viscosity, is not a fluid property, but depends on the local state of flow or turbulence; μ T is a scalar quantity, which varies significantly within the flow domain; and 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 G k is the generation of turbulence kinetic energy due to the mean velocity gradients; G b is the generation of turbulence kinetic energy due to buoyancy; Y M is the contribution of the fluctuating dilation in compressible turbulence to the overall dissipation rate; C ε1 , C ε2 , C ε3 are constants which are varied based on the application. Pr À1 eff , k and Pr À1 eff , ε are the inverse effective Prandtl numbers for k and ε, respectively derived analytically using the Random-Number Generator (RNG) theory. In addition, S k and S ε are user-defined source term. The equations used to calculate some of the above parameters are: Pr À1 eff À 1:3929 α 0 À 1:3929 where α 0 ¼ 1, ( μ mol ) is the molecular viscosity. In the high Reynolds number limit 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].

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

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 liquidphase, 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

Bubble-induced turbulence model
The bubble-induced turbulence is represented by introducing two source terms, S k and S ε , into the k-ε Eq. [36,37] as: where q LG is the covariance of the liquid-phase and the dispersed gas-phase velocities; u r and u 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 (C D and C VM ) appear in the formulation. The equations for the relative and drift velocities are: The drift velocity (u ! d ) is a statistical quantity due to the conditional averaging. It may not be negligible as it accounts for the dispersion effects due to the bubbles transport by the turbulent fluid motion. D . For the purpose of practical computations [38], the dispersion tensor is simplified to its diagonal form [37,[39][40][41] as: where I is a 3 Â 3 identity matrix

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 P J R J i represents the net change in the number density distribution due to the coalescence and breakup of bubbles; B 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 (f BV Þ: where d I and d II represent the diameters of the daughter bubbles in the binary breakage of a parent bubble of diameter d and volume v. The breakup rate of the bubbles can be represented as: where ξ min ¼ λ min =d and ξ ¼ λ=d, represent the size ratio between an eddy and a bubble in the inertial subrange and β is a constant derived from fundamental consideration and it equals 2. Subsequently, the birth and death rates due to breakup are represented as: 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: u T is the turbulent velocity in the inertial subrange of isotropic turbulence and can be estimated as: [46] 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 h o and h f represent the initial bubble film thickness and the critical film thickness. Their values were reported to be 10 À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: where χ ij is coalescence rate of bubbles (m À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.

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: λ S is the solids bulk viscosity, which describes the resistance of the particle suspension against compression and is written 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].
g 0 is the radial distribution function, accounting for the increase in the probability of collisions when the particle density increases; and can be expressed as [26,48]: where α S, max ¼ 0:62, and beyond this value, the radial distribution function goes to infinity.
Θ is the granular temperature, which is a measure of the kinetic energy contained in the fluctuating velocity for the solid particles and it is defined using the algebraic model by Ding and Gidaspow [48], which helps minimize the computational load by avoiding to solve an additional differential equation along with its closure models: μ S is the solid-phase shear viscosity, which is an elastic force, arising from the solids in response to shear, compression and extension stresses exerted on it by the continuous liquidphase. It should not be confused with the viscous forces, which arise during the fluid flow. The fluid viscosity is proportional to the rate of deformation over time, whereas, the solid viscosity is proportional to the amount of shear deformation. The following model by Gidaspow [14] is used to describe the solid-phase viscosity:

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 masstransfer 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 k L a is the volumetric liquid-side mass transfer coefficient, which is represented using an empirical correlation, Eq.
Moreover, when accounting for chemical reactions, an additional species conservation equation has to be considered as follows: where y i k represents the mass fraction of species i in phase k. Also, r i k and r ij kn represent the F-T reaction kinetics rate and the rate of chemical absorption, respectively.

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]:

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 Δ is the mesh spacing and α = 2, 3 for 2 and 3 dimensions.

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 (φ x; t ð ÞÞ is the related to the liquid volume fraction field function as shown in Eq. (43): where H c φ À Á is referred to as the Heaviside function, is represented as [52]: where, ε is a parameter in the order of the mesh size (typically 1.5Δ). It is important to note that a re-distancing algorithm is required to regularize the function, as the numerical errors accumulate. Also, an extra re-distancing function, preserving ∇φ ¼ 1 around the zero level of φ, is required.

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.

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 G x À x 0 j j ð Þis an appropriate spatial filter function: 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 ϕ x 0 ; t ð Þ can thus be described using Eq. (49), where g ϕ x 0 ; t ð Þ is the filtered resolvable component and ϕ x 0 ; t ð Þ is the sub-grid scale (SGS) component, which accounts for the unresolved spatial variations at lengths smaller than the filter width.
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 τ k 00 is decomposed into both resolved and SGS components as shown in Eq. (53).
The parameter τ k 00 ij is an unknown SGS tensor, representing the effects of SGS motion on the resolved fields of the LES simulation, and has to be represented using an SGS model. The prime objective of the SGS models is to accurately capture small dissipative scales by representing kinetic energy losses due to the viscous forces by accounting for their effect in a statistical manner. Numerous SGS models were developed and most of them are based on the Boussinesq hypothesis [56][57][58][59]. A brief description of the basic and dynamic SGS models is given in the following section.

Basic SGS model
In this model, the unresolved stress tensor (τ k 00 ij Þ, can be expressed according to the singlephase counterpart described in Eq. (54).
where μ SGS is the SGS eddy viscosity and Š ij is the strain rate of the large-scale resolved field expressed using Eq. (55).
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 is an empirical constant (Smagorinsky constant), which depends on the type of flow. It equals 0.18 for isotropic turbulence and 0.10 for flow near solid walls. For the multi-fluid approach, the value of the Smagorinsky constant for the continuous phase typically equals (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].

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 G x À x 0 j j ð Þand a test filter G x À x 0 j j ð Þ, where the latter has a larger width than the former. The component-weighted volume-averaging for the multi-fluid approach using the test filter can be defined using Eq. (57) for the continuous phase. When this test filter is applied to the instantaneous momentum equation, it results in Eq. (58).
where L is difference between the grid and test level stress tensors represented by Eq. (59): The test filter stress (T k 00 ) is expressed using Eq. (60) 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 d can be found elsewhere [64][65][66].

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 Re 3 L .
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].

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: 1. 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].

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

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

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]. Figure 8. Snapshots of local turbulence intensity contours at different times [2].
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.

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.