ARTICLE pubs.acs.org/IECR
Computational Modeling of Trickle Bed Reactors Zeljko V. Kuzeljevic*,† and Milorad P. Dudukovic Chemical Reaction Engineering Laboratory, Department of Energy, Environmental and Chemical Engineering, Campus Box 1180, One Brookings Drive, Washington University in St. Louis, St. Louis, Missouri 63130, United States ABSTRACT: In this study, a three-dimensional Eulerian computational fluid dynamics (CFD) reactive flow model is developed for the trickle flow regime in a packed bed. The experimental study of porosity distribution in a bed of extrudates has been performed and then used to implement the local porosity values in the CFD computational grid. The reactive flow study addressed gas and liquid reactant limited systems. For each case a closed form approach of coupling the bed and particle scale solution within the CFD framework is presented. Model predictions achieved satisfactory agreement with experimental data for liquid holdup, wetting efficiency, and reactant conversion.
’ INTRODUCTION Trickle bed reactors (TBRs) are packed beds of catalyst particles in which reactants in the gas and liquid phase are introduced at the top of the column and flow cocurrently down. They are used in petroleum and refinery processes such as hydrodesulfurization and hydrogenation, but also find application in oxidation of organics in wastewater effluents, volatile organic compound abatement in air pollution control, and enzymatic reactions.1 The use of TBRs offers many advantages: they operate at conditions that are close to plug flow, have high catalyst to liquid ratio which is useful in abating the homogeneous side reactions, high throughput range, and there is no danger of flooding. On the other hand, TBRs also have serious drawbacks. They are prone to liquid mal-distribution and incomplete catalyst wetting. This reduces the extent of catalyst utilization and, for the case of highly exothermic reactions, can lead to hot spots and reactor runaway.2,3 Also, to reduce the pressure drop, catalyst particle diameters are usually in the range of a couple of millimeters. Thus, in TBRs intraparticle diffusion effects can play a significant role.4 Currently, Eulerian computational fluid dynamics modeling (CFD) is employed as an important tool for performace assessment, development, and scale-up of TBRs. The underlying equations, obtained by volume averaging of momentum and mass conservation equations, have a very attractive form which does not require detailed geometry of the system as an input.5 However, the system of equations needs to be closed with descriptions of interphase momentum exchanges, i.e., the closure equations.6 Usually, these are grouped into phase interaction closures (i.e., gasliquid, gassolid, and liquidsolid interactions) and capillary pressure closure. The first group gives an estimate of momentum transfer from/to each phase due to drag exerted at the interfacial areas. The second gives the pressure difference between the nonwetting (i.e., gas) and wetting (i.e., liquid) phase, which arises due to the nonzero value of curvature of the interfacial area between the two phases.7 It is evident that a physical picture of the system introduced through closures and a proper description of the pore space are crucial for proper modeling. Studies811 showed a great potential for CFD-based modeling of TBR. However, they also showed the effect and importance of the choice of interaction and capillary closures and of the r 2012 American Chemical Society
implementation of a porosity map on the computational domain.12 For example, introducing the capillary term leads to the smearing of the liquid phase front as it travels down the bed and to the “diffusion” of liquid phase from the region of higher porosity into the region of lower porosity. The phase interaction closures reduce the speed at which the liquid front travels and increase the liquid saturation of the front.10,12 Gunjal and Ranade13 studied the hydrodesulfurization in a pilot plant TBR using CFD, thus expanding the use of the model to the reacting systems. Alopaeus et al.14 based their interactions on the slit model and developed the CFD approach that is applicable to Newtonian and non-Newtonian fluids, high and low pressure operation, and co- and countercurrent operation. Other representative CFD studies of TBR have been performed.1521 Comprehensive reviews of TBR phenomenological and CFD modeling are available in Kuipers and Van Swaaij22 and Carbonell.23 The details of the derivation of the volume averaged governing equations are given in Anderson and Jackson,24 Drew,6 and Ishii.5 In this study, we focus on developing a robust hydrodynamic Eulerian CFD model, and the use of the hydrodynamic solution predicted by the model together with species balances to yield the TBR performance assessment. In the Experimental section the measurements for statistical characterization of the porosity map and hydrodynamic parameters (liquid holdup) are presented. The results obtained are then used to implement the local porosity values in the CFD grid. The reactive flow study addresses gas and liquid limited systems and for each case a closed form approach of coupling bed and particle scale solution within CFD framework is presented and assessed against available experimental data.
’ EXPERIMENTAL Figure 1 is a schematic of the high pressure trickle bed reactor setup. It consists of a 16.3-cm inner diameter stainless steel Special Issue: Nigam Issue Received: April 7, 2011 Accepted: December 20, 2011 Revised: December 14, 2011 Published: January 18, 2012 1663
dx.doi.org/10.1021/ie2007449 | Ind. Eng. Chem. Res. 2012, 51, 1663–1671
Industrial & Engineering Chemistry Research
ARTICLE
Table 2. CT Scan Axial Positions position
distance from the packing
distance from the top
support mesh, cm
of the packing, cm
top
65
middle
35
bottom
2.5
5.0 35 67.5
Table 3. Governing Equations of the Eulerian CFD Model ∂ ðεα Fα Þ þ ∇ 3 ðεα Fα uBα Þ ¼ 0 ∂t
(1)
∂ ðεα Fα uBα Þ þ ∇ 3 ðεα Fα uBα uBα Þ ¼ εα ∇Pα þ εα Fα gB þ ∇ 3 τ α þ B F α (2) ∂t ~ Fα ¼
Figure 1. Experimental setup.
τ α ¼ εα μα ∇uBα þ ð∇uBα ÞT
Table 1. Experimental Setup and Operating Conditions column diameter, m
0.163
bed height, m
0.70
fluids used
water and air
packing
alumina extrudates (trilobes); 1.9 mm
bed porosity system pressure, barg
0.39 07
equivalent diameter
gas velocity, mm/s
30200
gas mass flux, kg/m2s
0.0361.68
liquid velocity, mm/s
1.99
liquid hourly space velocity, 1/hr
1046
reactor column, gas and liquid delivery system, and gasliquid separator. The setup enables the measurement of the overall liquid holdup and pressure drop. With the aid of gamma-ray computed tomography (CT), the cross-sectional distribution of porosity can be obtained at the designated axial positions. Liquid is stored in a pressurized tank and supplied via a pump to the distributor at the top of the column. The liquid and gas flow rates are controlled using the needle valves and measured using the flow meters placed in the inlet lines. After the gasliquid separation in a bypass separator the liquid is directed back to the tank. The column was set on the Arlyn weight scale thus enabling the measurement of the overall liquid holdup using the weighting method. Cross-sectional porosity maps are obtained using gamma-ray CT. For this purpose, the reactor column is placed between the gamma-ray source (Cs137) and 11 NaI scintillation detectors. During measurements, the source and the detectors rotate around the column and at each position the attenuation of gamma-rays detected by the 11 detectors is recorded with the aid of a data acquisition system connected to a PC. Then, a reconstruction procedure based on the expectation-maximization (EM) algorithm is used for image reconstruction. More details on experimental setup and operating procedure for the high pressure TBR is given in Kuzeljevic,25 and Kuzeljevic et al.26 CT hardware, data acquisition system, and reconstruction procedure are discussed in Kumar,27 Kumar and Dudukovic,28 Roy,29 Chen et al.,30 and Kuzeljevic.25 The principles of EM algorithm are given in Dempster et al.,31 while the application of the algorithm for the transmission tomography is discussed in Lange and Carson.32
n
∑ Kβα ðuBβ uBα Þ β¼1
(3)
(4)
Table 1 summarizes the experimental conditions. The bed was packed with 1.9-mm porous alumina extrudates (trilobes) to a height of 0.70 m. The resulting average porosity of the bed was 0.39. The fluids used were water and air. CT unit was employed to examine the cross-sectional distribution of porosity at three axial positions specified in Table 2. In the remainder of the text, we refer to these positions according to the values given in the second column of Table 2, i.e., as the distance from the packing support mesh. The experimental procedure was as follows. The bed was flooded and left overnight to ensure that catalyst particle pores were filled with water. Afterward, the bed was drained until only static liquid holdup remained and the system was pressurized. Weight scale was zeroed and gas flow initiated and set to the operating value. To minimize the hysteresis effects on the measurements,33 the Kan-liquid prewetting procedure was used. Thus, first, the liquid flow rate was increased until the pulsing flow regime was reached. Afterward, the liquid flow rate was set to the operating value. The measurements of liquid holdup (weight scale reading) were performed once the steady state flow conditions were established.
’ MODEL DESCRIPTION Governing Equations. Governing equations of the Eulerian CFD model are the volume averaged mass and momentum conservation equations given in Table 3. Closures. Phase interaction closures based on the model of Attou et al.34 were used (Table 4). Capillary closure is given by Grosser et al.35 model; i.e., correlation of Leverett J-function36 for the conditions of trickle flow. The resulting expression is equation (8).
2 0 3 pffiffiffi 1=2 2σ 1 εB 1=3 4 1 3 5 1 þ @ Þ Pc ¼ dp 1 εG 2 π
ð8Þ
Implementation of Experimental Porosity Distribution. The need for implementation of porosity profile on the computational domain has been recognized first in the studies of packed beds with single phase flow.37 Such approach has been further developed 1664
dx.doi.org/10.1021/ie2007449 |Ind. Eng. Chem. Res. 2012, 51, 1663–1671
Industrial & Engineering Chemistry Research
ARTICLE
Table 4. Phase Interaction Closures34 momentum exchange coefficient, kg/m3s
phases gasliquid
gassolid
liquidsolid
KGL ¼
E1 ð1 εG Þ2 1 εB 2=3 E2 ð1 εG Þ 1 εB 1=3 μ þ FG uBG uBL G dp εG dp2 1 εG 1 εG
(5)
KGS ¼
E1 ð1 εG Þ2 1 εB 2=3 E2 ð1 εG Þ 1 εB 1=3 μ þ FG uBG G dp εG dp2 1 εG 1 εG
(6)
KLS ¼
E1 ð1 εB Þ2 E2 ð1 εB Þ μL þ FG uBL dp εL dp2
and incorporated into the CFD models of TBRs.810,12,13,1518 We extend this previously proposed methodology to the threedimensional (3D) computational grid used in this study. In general, the observed porosity distribution is dependent upon the length scale (resolution) at which the measurements are performed. Typically the high-resolution (finer than the particle size) imaging yields the bimodal distribution, i.e., each pixel in the image is either occupied by solid phase or represents a void (see images obtained in a study of Baldwin et al.38 or van der Merwe et al.39). Axial averaging of such data (or equivalently, the use of techniques with coarser resolution) yields the typical oscillatory pattern of porosity in the near wall region that extends about 510 particle diameters toward the interior of the bed. For example, the experimental data of Roblee et al.,40 and Benenati and Brosilow41 have served as a basis for a number of correlations4244 predicting such damped oscillatory radial porosity profile. Also, exponential type correlations have been proposed.45,46 They predict the increase of porosity value toward the wall but without any oscillatory pattern. Judging based on available experimental data for the porosity profile40,41,47,48 a damped oscillatory pattern seems to exist in the near wall region. On the meso-scale, as shown in the experimental results below, overall bed porosity distribution is pseudo-Gaussian. Hence, the following methodology was proposed9,12,13,15,17 to assign porosity map in the CFD grid so that it follows Gaussian distribution. Mueller44 correlation (eq 9) is used to generate the longitudinally averaged radial profile of porosity.
εB ðrÞ ¼ εB þ ð1 εB ÞJ0 ðar Þexpð br Þ ða, bÞ ¼ f ðDc , dp Þ r ¼ r=dp
ð9Þ
Resulting profile is then integrated with respect to radial position to yield the averaged porosity in 10 radial sections (eq 10). Z 1 ÆεBi æ ¼ end init ri ri
riend riinit
εB ðrÞdr
ð10Þ
In each radial section, Gaussian distribution is imposed around the calculated mean with the standard deviation value as obtained based on experimental results of CT imaging given below. Note that, in principle, similarly to the use of Muller44 voidage profile, an exponential correlation45,46 could have been integrated to yield the average porosity values in each radial section of the CFD gird. Clearly, the entire procedure for porosity implementation would still hold. The use of the correlation that captures damped oscillatory behavior brings the model in line with the majority of the experimental results for porosity variation in the near wall region which may be important to capture if
(7)
finer meshing is done. At the resolution used in this study use of an exponential profile is also adequate. Species Balance and Particle Scale Models. The general expression for species balance adapted to our Eulerian CFD model is given by eq 11. The last term on the RHS represents the creation/disappearance of species i in chemical reactions. Its proper representation, in general, requires consideration of interphase and intraparticle mass transfer of reacting species. We here introduce the particle scale equations to account for mass transfer and partial wetting in the reactive flow CFD model. ∂ ðF Yi, α Þ þ ∇ 3 ðFα uB uBα Yi, α Þ ∂t α ¼ ∇ 3 ðFα Di, m ∇Yi, α Þ þ
NR
∑ Ri, r r¼1
ð11Þ
The two general cases considered here are (i) liquid reactant limited reaction with nonvolatile liquid phase, and (ii) gas reactant limited system. Both systems are considered isothermal; i.e., interphase and intraparticle thermal gradients are considered negligible. The particle scale implementation for case (i) is based on the exact solution for the partially wetted, infinitely long cylinder particle49 as given in Table 5. The hydrodynamic CFD solution is used to prescribe the local flow field, and values of wetting efficiency and rates of interphase mass transfer (see Table 6). The particle scale implementation for case (ii) is based on the approximate solution as proposed by Beaudry et al.50 To make analytical or numerical solution feasible, the existence of dry, wetted, and half-wetted catalyst particles is assumed. The solution is then sought for each case and overall catalyst efficiency is obtained as the weighted average of the three contributions as given in eq 25. η ¼ ð1 ηCE Þηd þ 2ð1 ηCE ÞηCE ηhw þ η2CE ηw
ð25Þ
The needed estimate of mass transfer parameters was obtained using correlations given in Table 6; however, experimental values are used where available. The CFD approach used in this study does not yield the value of wetting efficiency as the direct outcome of simulation since no particle scale geometry is assumed and gasliquid interface is not resolved. Hence, it is only possible to predict wetting efficiency indirectly.9 We utilized the correlations given in Table 6 that predict wetting efficiency based on the values of Reynolds number, Galileo number, and pressure drop gradient. In turn, dimensionless numbers and pressure drop are part of simulation predictions and in that sense our simulation and model affect the predicted values of wetting efficiency. Wetting efficiency and porosity are meso- and macro-scale parameters, but are here applied in each of the computational cells in the grid. However, it should be noted that all the 1665
dx.doi.org/10.1021/ie2007449 |Ind. Eng. Chem. Res. 2012, 51, 1663–1671
Industrial & Engineering Chemistry Research
ARTICLE
Table 5. Solution to the Particle Scale Model for Liquid Limited Reaction for with Nonvolatile Liquid Phase49 ∞
∑ Ai, j ej j¼1
(12)
¼ Gi
uðξ, kÞ ¼
e1 I0 ðϕξÞ þ ϕ I1 ðϕÞ
∞
∑ n¼2
en cos½ðn 1Þθ In þ 1 ðϕÞ þ n1 ϕ In ðϕξÞ
(13)
A11 ¼ π º1 þ ηCE ðt1 2 1Þß Ai1 ¼
(14)
sin ºði 1ÞπηCE ß ðti tj 1Þ i1
(15)
sin ºði 1ÞπηCE ß ðti tj 1Þ, 2 e i e ∞ i1 " # 1 sin½ði þ j 2ÞπηCE sin½ði jÞπηCE þ ðti tj 1Þ, Aij ¼ 2 i þ j2 ij
(16)
A1i ¼
Aii ¼
π 1 sin½2ði 1ÞπηCE 2 þ πηCE þ ðti 1Þ, 2 2 2ði 1Þ
G1 ¼ π ηCE t1
Gi ¼ ti
sin ºði 1ÞπηCE ß , i1
Ii ðϕÞ þ i 1, ti 1 ¼ ϕ Ii1 ðϕÞ
(17)
2eie∞
(18)
2eie∞
(19)
2eie∞
(20)
1eie∞
Table 6. Correlations Used in the Reactive Flow CFD Model kLS ¼ 4:25 Di, m ¼
Di, m εp 0:48 0:33 Re ScL dp ηCE L
7:4•10 10 ð2:6MB Þ1=2 T μB VA0:6
GaL 0:071 ηCE ¼ 1:617 Re0:146 L ηCE
Tan and Smith58
(21)
Wilke and Chang57
(22)
El-Hisnawi56
1=9 1=3 1 þ ΔP=Δz=ðFL gÞ ¼ 1:104ReL GaL
(23) 51
Al-Dahhan and Dudukovic
(24)
CT data were used to calculate the values of sectional porosity using the codes written in Matlab. Values of porosity were then patched into the solver using user defined functions (UDFs) written in C programming language. Phase interaction closures have been specified using UDFs. Capillary closure has been introduced as the source term in the momentum equations and specified using UDFs. Particle scale and interphase mass transfer models were specified using UDFs within species transport module of Fluent 6.2 Eulerian model. Initial conditions were as given in eq 26. On the cylindrical column wall the no-slip boundary condition is used. By using three different cell sizes in the computational grid it was verified that 40,000 cells on the domain are sufficient to achieve grid independent results. Figure 2. CFD model solution strategy.
equations are set in the meso-scale domain and, at the same time, all needed parameters are calculated using the meso-scale framework. Solution Strategy. An unstructured, three-dimensional computational grid representing the TBR column geometry has been created in GAMBIT. Unsteady state model is solved using Fluent (Figure 2) until steady state is reached. Mueller44 correlation and
εL ¼ ε0L εs ¼ 1 εB ¼ const uG ¼ 106 m=s uL ¼ 106 m=s uS ¼ 0 ¼ const
ð26Þ
Inlet boundary condition is set to mimic the top liquid distributor. The distributor is a 2.5-cm (1 in.) high box directly 1666
dx.doi.org/10.1021/ie2007449 |Ind. Eng. Chem. Res. 2012, 51, 1663–1671
Industrial & Engineering Chemistry Research
ARTICLE
Table 7. Porosity Distribution Parameters standard deviation
overall porosity
of porosity on the cross-sectional domain, %
(weight method)
position, cm
average porosity
65
0.40
14.8
35
0.41
16.1
2.5
0.40
13.9
0.39
Figure 3. Implementation of liquid velocity at the inlet boundary. Colorbars represent liquid velocity in m/s.
Figure 5. (a) Representative radial section (0.90 e r/(Dc/2) e 1.0) in the porosity map on the CFD computational domain. In each radial section, Gaussian distribution is imposed around the mean obtained by the integration of the Mueller44 correlation. (b) Resulting porosity distribution map for the entire computational domain.
’ RESULTS AND DISCUSSION Figure 4. Experimental results of porosity distribution for middle axial position (z = 35 cm): (a) cross-sectional porosity map; (b) radial profile of porosity; (c) porosity distribution.
connected to the liquid injection tube. The liquid injected in the distributor exits through 240 holes of 0.9-mm diameter. The air injected at the top of the column flows through 19 tubes (6.5-mm I.D) passed axially though the distributor. Hence, the inlet BC is set as shown in Figure 3.
Porosity Distribution. Porosity (voidage) distribution of the packed bed was characterized at three axial positions (see Table 2) using CT experimental results. The resolution was on the order of couple of particle diameters and the observed porosity distribution resembles a Gaussian one. The results are presented in Figure 4 and Table 7. Figure 4a presents the cross-sectional distribution of porosity at the middle axial location. Figure 4b shows the radial distribution of porosity at the middle axial location, while Figure 4c displays the values in a histogram. The results for the standard 1667
dx.doi.org/10.1021/ie2007449 |Ind. Eng. Chem. Res. 2012, 51, 1663–1671
Industrial & Engineering Chemistry Research
ARTICLE
Figure 6. Simulated and experimental values of liquid holdup in the bed of extrudates. (a) P = 4 barg, UG = 70 mm/s, (b) P = 1 barg, UL = 4.53 mm/s.
deviation of porosity in the three cross-sections (Table 7) are in agreement with the above-mentioned studies. Also, the results do not indicate a significant dependence of porosity distribution and its parameters on the axial position. Average value of the cross-sectional porosity obtained using CT is in good agreement (deviation
