A Computational Transport Model for Wall-Cooled Catalytic Reactor

Mar 22, 2008 - A computational transport model is proposed for simulating the concentration, temperature, and velocity distributions in a wall-cooled ...
2 downloads 0 Views 352KB Size
Download PDF
2656

Ind. Eng. Chem. Res. 2008, 47, 2656-2665

A Computational Transport Model for Wall-Cooled Catalytic Reactor G. B. Liu, K. T. Yu, X. G. Yuan,* and C. J. Liu State Key Laboratory for Chemical Engineering (Tianjin UniVersity), Chemical Engineering Research Center and School of Chemical Engineering and Technology, Tianjin UniVersity, Tianjin, 300072, China

A computational transport model is proposed for simulating the concentration, temperature, and velocity distributions in a wall-cooled fixed-bed catalytic reactor producing vinyl acetate from acetic acid and acetylene without assuming empirical Sct and Prt numbers or using the experimental coefficient of dispersion. The model consists of the basic differential equation of mass transfer with the recently developed c2 - c model for its closure and the accompanied computational fluid dynamics and heat transfer models. With the present model, the radial and axial profiles of turbulent mass-transfer diffusivity Dt, turbulent thermal diffusivity Rt, and turbulent kinematic viscosity νt, as well as their distributions in the reactor are obtained. It is found that the local ratios of νt/Dt and νt/Rt, or the local Prandtl (Prt) and Schmidt (Sct) numbers, are varying in the reactor, although the shapes of the νt, Rt, and Dt profiles show somewhat similarity. The simulated results are found to be in agreement with the experimental measurements reported by Valstar et al.33 Introduction The fixed-bed reactors are the most commonly used equipment for conducting industrial heterogeneous catalytic reactions, such as the carbon monoxide conversion and ammonia synthesis in the basic chemical industry, the ethylene oxide and vinyl acetate synthesis in the petrochemical industry, and the catalytic reforming and polymerization in the petroleum refining process.1 The studies concerning the performances of such kinds of reactors have been extensively reported.1-4 The 1D plug-flow model is first used for reactor design and analysis, where the concentration and temperature gradients are assumed only to occur in the axial direction.5 Later, the 1D flow model with axial mixing is introduced to take into account the mixing effect that is influential to the temperature and concentration gradients as well as the reactor performances.6 It was reported7-17 that the porosity in fixed-bed displays uneven distribution, especially when the ratio of reactor tube diameter to particle diameter is less than 8. Further studies18-21 showed that the momentum, heat, and mass transports can be modeled correctly only if the uneven porosity was considered. With the increase of computer power, modeling works22,23 on the fluid flow and heat transfer in the fixed bed have been done using the interstitial CFD approach, and the detailed velocity and temperature profiles were given and validated with the experiments. Among these, simulation was done so far under the condition that at most 44 balls were structurally arranged into a tube. However, this is quite different from the industrial reactors where a large number of catalyst or specific packing particles are packed, and an interstitial CFD approach may lead to a considerable computation burden even with powerful computation facilities available. At the same time, some investigations used the 2D pseudo-homogeneous model with the consideration of the radial distribution of porosity20,21 and gave good simulated results. The advancement of applying a pseudo-homogeneous CFD model to the reactor design has enabled us to calculate the velocity profile, whereas the temperature and concentration distributions were obtained using either the empirical Prandtl (Prt) and Schmidt (Sct) numbers or * To whom correspondence should be addressed. Tel.: +86 22 27404732. Fax: +86 22 27404496. E-mail: [email protected].

the experimental coefficient of dispersion for substituting the effective diffusivities of heat and mass transfer.1-3 Nevertheless, the empirical dispersion correlations are often obtained using the inert tracer technique. In fact, such empirical correlations, even if available, are not always applicable to different reactors concerned because the effective diffusivities are dependent on many factors, such as bed structure, operating conditions, and many others.24,25 In this article, a computational transport model, which is based on the computational mass transfer methodology recently developed,26-28 is proposed for predicting the velocity, temperature and concentration distributions in fixed-bed reactor, without using the empirical Sct and Prt numbers or the experimental dispersion coefficient. The proposed model consists of the basic differential equation of mass transfer with the c2 - c equations26,29 for its closure and the accompanied equations of computational fluid dynamics (CFD) and computational heat transfer (CHT). The closures of CFD and CHT differential equations are made by employing the k - 30,31 and t2 - t32 equations. The simulation in this work is made to a wall-cooled catalytic fixed-bed reactor for the synthesis of vinyl acetate from acetic acid and acetylene as reported by Valstar et al.33 Computation was implemented by using commercial software FLUENT 6.2. The simulated results on velocity, temperature, and concentration distributions in axial and radial directions as well as the conversion of the catalytic reaction are compared with the experimental data for testing the model reliability. Model Formulations The processes undertaking simultaneously in a catalytic reactor include the reactive mass transfer, the momentum transfer, and heat effect from the chemical reaction. The following assumptions are applied for the present simulation: (1) The reactor operates in a steady state and is in pseudohomogeneous flow; (2) The fluid is incompressible, and the flow is axially symmetrical; (3) The catalyst temperature is equal to the fluid temperature, and the coolant temperature at the wall is constant; (4) The activity of the catalyst remains unchanged. The modeling equations should involve both the differential equations of mass transfer with chemical reaction and the

10.1021/ie070737y CCC: $40.75 © 2008 American Chemical Society Published on Web 03/22/2008

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2657

accompanied equations of momentum and heat transfer as shown in the following sections. The Equation of Mass Transfer with Chemical Reaction. The mass-transfer equation set consists of the differential

((

(

)

(

)

where Mi is the molar weight of component i, Rs is the apparent reaction rate, which will be given in the subsequent section, and Fb is the bulk density of catalyst. In eq 2, the negative and the plus signs refer to the reactant and product components, respectively. The Di,eff in eq 1 is the effective diffusivity of mass transfer, represented as,

Di,eff ) Di,L + Di,t

(3)

where Di,L is the molecular diffusivity and Di,t is the masstransfer diffusivity induced by the turbulent flow. The turbulent diffusivity Di,t can be calculated on the basis of the Boussinesq hypothesis from the fluctuating concentration variance c2i and its dissipation rate i,c:26,27,29

( )

2 k ci Di,t ) Cc0k  i,c

(4)

- 2γFi,c

Di,L +

2

∂Ci ∂Ci + ∂x ∂r

(7)

Di,t ∂(γFi,c) σ c ∂x

Di,L +

2

i,c ci2

Di,t ∂(γFi,c) r σ c ∂r

- Cc2γF

i,c2

c2 i,c (8) Cc3γF k

The constants in eqs 4, 7, and 8 are34 Cc0 ) 0.11, Cc1 ) 1.8, Cc2 ) 2.2, Cc3 ) 0.8, σc ) 1.0, σc ) 1.0. The Fluid Dynamics Equations. The fluid dynamic equations representing the momentum transfer include the continuity equation, momentum equations, and the standard k- model equations based on the axial symmetrical assumption. The continuity equation:

∂(γFu) 1 ∂(γFrV) + )0 ∂x r ∂r

(9)

The axial and radial momentum equations:

∂(γFuu) 1 ∂(rγFuV) ∂ ∂u 2 + γµeff 2 - (∇‚u) ∂x r ∂r ∂x ∂x 3 1 ∂ ∂u ∂V rγµeff + r ∂r ∂r ∂x

[ (

)] [ (

)]

∂p - γ(R0u + Fg) ∂x

) -γ

(10)

∂(γFuV) 1 ∂(rγFVV) 1 ∂ ∂V 2 rγµeff 2 - (∇‚u) + ∂x r ∂r r ∂r ∂r 3 ∂ ∂p ∂u ∂V V 2γµeff γµeff ) -γ - 2γµeff 2 + + (∇‚u) ∂x ∂r ∂x ∂r 3r r γR0V (11)

[

[ (

1/2

)

2

∂Ci ∂Ci + ∂x ∂r

1 ∂ r ∂r

) Cc1γFDi,t

(2)

)

[( ) ( ) ] (( ) ) (( ) ) [( ) ( ) ]

∂Ci 1 ∂ γrFDi,eff + Si (1) r ∂r ∂r

Si ) ( MiRsFb

((

2

) 2γFDi,t

∂(γFui,c) 1 ∂(γFrVi,c) ∂ + ∂x r ∂r ∂x

where Ci is the Reynolds-average concentration for component i, Si is the source term representing the mass of component i generated by the chemical reaction and can be calculated from the reaction rate,

)

Di,t ∂(γFci2) 1 ∂ Di,L + r r ∂r σc ∂r

equation of mass transfer for component i and the proposed c2 - c two-equation model for the closure. On the basis of the axial symmetrical assumption, the mass-transfer equation can be written as,

∂(γFuCi) 1 ∂(γrFVCi) ∂Ci ∂ γFDi,eff + + ) ∂x r ∂r ∂x ∂x

)

Di,t ∂(γFci2) ∂(γFuci2) 1 ∂(γFrVci2) ∂ Di,L + + ∂x r ∂r ∂x σc ∂x

(

)]

)]

where where c2i and i,c are defined as follows:

ci2 t cici i,c t Di,L

( ) ∂ci ∂ci ∂xj ∂xj

∇‚u ) (5) (6)

The two-equation c2i - i,c model for the closure of eq 1, consisting of the equations of fluctuating concentration variance and its dissipation rate, are given below:

∂u 1 ∂(rV) + ∂x r ∂r

µeff ) µ + µt µt ) FCµ

k2 

(12) (13) (14)

In the equations, F is the fluid density; u is the interstitial velocity vector; u and V are the axial and radial velocities respectively in u; µ, µt, and µeff represent the molecular, turbulent, and effective viscosities of the fluid respectively, and γ is the void fraction that is related to the dimensionless distance z from the center to the reactor wall as given by de Klerk’s,17

2658

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

γ ) 2.14z2 - 2.53z + 1, z e 0.637

(15)

γ ) γ∞ + 0.29 exp(- 0.6z) cos(2.3π(z - 0.16)) + 0.15 exp(- 0.9z), z > 0.637 (16) where R0 denotes the resistant coefficient of porous media composed of a catalyst, which can be calculated by the modified Ergun equation:35

R0 ) 150µ

(1 - γ) 2

2

γ de

2

+ 1.75F

(1 - γ) γ2de

|u|

(17)

The kinetic energy k equation:

(( ) )

µt ∂k ∂(γFuk) 1 ∂(rγFVk) ∂ γ µ+ + ∂x r ∂r ∂x σk ∂x

(( ))

µt ∂k 1 ∂ γr µ + r ∂r σk ∂r

{ [(∂u∂x) + (∂V∂r) + (Vr) ] + (∂u∂r + ∂V∂x) } - γF

) γµt 2

2

2

2

2

(18)

keff ) kfγks(1.0-γ)

The kf includes the contributions by molecular conductivity kL ) FcpR and turbulent conductivity kt ) FcpRt,

kf ) Fcp(R + Rt)

(( ))

µt ∂ 1 ∂ γr µ + r ∂r σ ∂r

{ [( ) ( ) ( ) ] (

)}

 ∂u 2 ∂V 2 V2 ∂u ∂V 2 + 2 + + + k ∂x ∂r r ∂r ∂x c2γF

keff ) [cpF(R + Rt)]γ ks(1.0-γ)

2 (19) k

books.1-3 The heat transfer equation and its closure t2-t equations are given below. The heat transfer equation in axial symmetry:

1 ∂ ∂ ∂T ∂ (γFcpuT) + (γrFcpVT) ) k + ∂x r ∂r ∂x eff ∂x 1 ∂ ∂T rk + Q (20) r ∂r eff ∂r

)

(

)

where T is the Reynolds-average temperature and Q is the heat source coming from the exothermic heat of the chemical reaction:

Q ) -FbRs∆Hr

(24)

The turbulent thermal diffusivity Rt, which is an important parameter in computing the heat transfer in the exothermic

Rt ) CD0k

() k t2  t

(25)

The t2-t model by Jones and Musonge32 is composed of the fluctuating temperature variance t2 and its dissipation rate t as shown below:

The constants in eqs 14, 18, and 19 are taken to be36 Cµ ) 0.09, σk ) 1.0, σ ) 1.3, c1 ) 1.44, and c2)1.92. The Heat Transfer Equations. This part of equations includes the heat transfer equation expressed in time-averaged temperature and its closure equations. Because the synthesis of vinyl acetate from acetic acid and acetylene is a strong exothermic reaction, and the excess heat should be removed quickly by the coolant flowing outside of the reactor, which leads to a large temperature gradient along the radial directions. The heat transfer in the reactor with the wall-cooled jacket has been studied and reported in a number of articles 4,37-39 and

(

(23)

where R and Rt are regarded as molecular and turbulent thermal diffusivities, respectively. From eq 22, the effective thermal conductivity can be expressed by

µt ∂ ∂(γFu) 1 ∂(rγFV) ∂ γ µ+ + ∂x r ∂r ∂x σ ∂x

) c1γµt

(22)

reactor, is determined by the t2-t model, in which Rt is defined by:

The kinetic energy dissipation rate  equation:

(( ) )

geometry and the fluid characteristics. The calculated result in the present case shows that better simulation can be achieved if the weighted geometric mean of the thermal conductivities ks and kf39,41 is used:

(21)

In eq 20, keff is the effective thermal conductivity, which involves the contributions by the thermal conductivities of the catalyst ks and the fluid kf. The effective thermal conductivity of a porous media saturated with flowing fluid depends on its

The fluctuating temperature Variance t2 equation:

(( ) ) [( ) ( ) ]

Rt ∂t2 1 ∂ ∂ ∂ γF R + + (γFut2) + (rγFVt2) ) ∂x r ∂r ∂x σt ∂x

( ( ))

Rt ∂t2 1 ∂ ∂T 2 ∂T rγF R + + 2γFRt + r ∂r σt ∂r ∂x ∂r

2

- 2γFt (26)

The fluctuating temperature Variance dissipation rate t equation:

(( ) )

Rt ∂t ∂ 1 ∂ ∂ (γFut) + (rγFVt) ) γF R + + ∂x r ∂r ∂x σt ∂x

( ( ) ) [( )

Rt ∂t t2 t 1 ∂ rγF R + - CD1γF - CD2γF r ∂r σt ∂r k t2 t ∂ui ∂T 2 ∂T 2 CD3γFRt + - CD4 uiuj (27) ∂x ∂r k ∂xj

( )]

The constants in eqs 25-27 are taken to be32 CD0 ) 0.10, CD1 ) 2.0, CD2 ) 0.9, CD3 ) 1.7, CD4 ) 1.4, and σt ) 1.0. It should be mentioned that the foregoing modeling formulation is isotropic. Strictly speaking, the flow as well as other transport properties in some chemical equipment, like catalytic reactor, are isotropic. Nevertheless, the use of an anisotropic model (i.e., anisotropic eddy-diffusivity) needs to add more complexity into the differential equations for CFD and other parts, which means much more computer work with the risk of not being convergent. As an approximation, in this article the conventional isotropic transport models were adopted for simplifying the computation.

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2659 Table 1. Properties of Reaction Mixture case. no

t0 °C

molar ratio

1 2 3 4

176.1 176.0 186.4 176.1

1.5 1.5 1.5 4.0

µ cp kL aver. mol. 10-5 J kg-1 Jm-1 weight G kg/kmol Kgm-2s-1 Nsm-2 K-1 K-1 s-1 39.6 39.6 39.6 32.8

0.242 0.186 0.242 0.200

1.372 1.369 1.376 1.375

1680 1680 1710 1800

0.0333 0.0333 0.0344 0.0380

Table 2. Catalyst Characteristics average length, mm average diameter, mm effective diameter de, mm specific external surface, m2g-1 specific internal surface, m2 g-1 bed porosity bulk density, kg m-3 particle density, kg m-3 thermal conductivity, Jm-1 K-1 s-1

5.4 2.8 3.3 0.00217 350 0.36 570-600 910 0.184

Figure 1. Boundary conditions setup.

(170/Rg) atm-1 at a molar ratio of acetylene/acetic acid ) 4. The heat of reaction, ∆Hr is expressed as a function of temperature:

Table 3. Coefficients of Eq 31 for Calculating Component Enthalpy42 component

A

B

vinyl acetate acetic acid acetylene

-298.4 -417.9 228.0

-6.9870E-2 -5.8243E-2 1.5754E-3

C 3.9316E-5 3.3466E-5 -3.5319E-6

Application of the Proposed Model to a Wall-Cooled Fixed-Bed Catalytic Reactor. For testing the applicability of the proposed model to the reactor, a simulation was made to a wall-cooled fixed-bed catalytic reactor as reported by Valstar et al.33 for the synthesis of vinyl acetate from acetic acid and acetylene with zinc acetate on activated carbon as a catalyst. The internal and external diameters of the reactor are respectively 0.041 and 0.0449 m with a 1 m height of the catalyst bed, and the gaseous mixture flows upward from the bottom of the reactor. The reactor tube is surrounded by a second tube with an internal diameter of 0.0725 m. Oil is pumped with a linear velocity of 3.2 m/s through the annular space between tubes. The oil temperature change is controlled to within 0.5 °C, which is used to maintain the temperature in the annulus. The radial average conversions and the temperature profiles along the radial direction at different axial positions were measured. The properties of the reaction mixture for the four cases tested are listed in Table 1, and the catalyst characteristics are listed in Table 2. Other detailed information about the experiments can be found from the published article.33 The overall chemical reaction of vinyl acetate synthesis is as follows: Zn(Ac)2

CH3COOH + CHtCH 98 CH3COOCHdCH2 + ∆Hr (28) The apparent reaction rate of foregoing reaction was reported by Valstar et al.:33

Rs ) k∞ exp(- E/RgT)pAC 1 + exp(- ∆H1/RgT) ×exp(- ∆S1/Rg)pHAc + Cr pVA

(29)

The parameters of eq 29 were given as follows: k∞ ) 5100 kmol/kg s atm, E ) 85 000 kJ/kmol, ∆H1 ) 31 500 kJ/kmol, ∆S1 ) -71 000 kJ/kmol K, Cr ) 2.6 atm-1 at a molar ratio of acetylene/acetic acid ) 1.5, and Cr ) exp(-70 000/RgT) exp-

∆Hr ) ∆f HVA,m - ∆f HHAc,m - ∆f HAC,m

(30)

∆f Hi,m ) A + BT + CT2

(31)

where i can be VA, HAc, or AC; the coefficients for the enthalpy of different components are listed in Table 3. Boundary Conditions. Inlet Condition. At the reactor bottom, x ) 0, the boundary condition is given by Liu et al.:28 u ) uin, V ) 0, k ) 0.003uin2,  ) 0.09k3/2/de, CAC ) CAC,in, CHAc ) CHAc,in, T ) Tin ) Tw, c2i ) 0, ci ) 0, t 2 ) (0.082∆T)2, t ) 0.4 (/k) t 2, and ∆T is set to be 1 K at the start of the computation. Axis Condition. At the axis of reactor, r ) 0, all of the variables Φ have a zero gradient due to the assumption of axial symmetry, that is, (∂Φ/∂r) ) 0. Outlet Condition. Because the ratio of the reactor length L to the diameter of catalyst de is about 300, the fluid flow is close to the fully developed condition at the reactor outlet. Therefore, the outlet boundary condition of zero normal gradients is applied for all of the flow variables except pressure. Wall Condition. At the reactor wall, r ) R, the no slip condition is applied to the fluid flow, thus all of the variables related to flow are set to be zero, that is, u ) V ) k )  ) c2i ) i,c ) t2 ) t ) 0. For the mass-transfer equation, the zero flux condition at the wall is adopted, that is, (∂Ci/∂r) ) 0 at r ) R. The flow behavior in the near-wall region is approximated using the standard wall functions.36 At the wall, let T ) Tw as reported by a number of researchers.1,4,43 Numerical Algorithm. The governing equations of the proposed model as applied to the reactor concerned were solved numerically by using the commercial software FLUENT 6.2 with finite volume method. The SIMPLEC algorithm was used to solve the pressure-velocity coupling problem in the momentum equations. The second-order upwind spatial discretization scheme was used for all differential equations. The convergence criterion was that all of the variables’ residues were less than 1.0-6. As a higher number of nodes in the numerical computation results in better simulation, case 1 of Valstar’s experiments was taken as an example for finding the optimal resolution. The simulation with 1000 nodes at the axial direction and 100 nodes

2660

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 2. Simulated (line) and measured (circle) conversion profiles along the height of the reactor.

Figure 3. Simulated and experimental measurement of temperature profiles along the radial direction at different axial positions of the reactor. Solid black line, present model with 50 × 550 grid; solid green line, present model with 100 × 1000 grid; dashed red line, model by Valstar et al. (Valstar et al., 1975); circles, measured data at x ) 0.15 m, x ) 0.31 m, x ) 0.45 m, x ) 0.61 m, and x ) 0.79 m.

at the radial direction was compared with the other choice of 550 nodes at the axial direction and 50 nodes at the radial direction. The results show that no substantial difference on the

radial temperature profiles between them was found (part a of Figure 3, denoted by the green lines for 1000 × 100 grids and black lines for 550 × 50 grids). Therefore, the latter mentioned

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2661

Figure 5. Axial relative velocity along the radial direction at x ) 0.45 m (case 1).

Figure 4. Temperature profile of the reactor (case 1).

grid format was used for the present simulation with a total number of about 27 500 quadrilateral cells. Simulated Results and Discussion The Radial Average HAc Conversion along the Height of Reactor. The simulated radial concentration distribution was averaged at a different height of reactor to find the average conversion along the axial direction. The simulated and measured HAc conversions by Valstar et al.33 versus the height of the reactor under different operating conditions are plotted in Figure 2. Satisfactory agreement between them is seen. The Radial Temperature Profile. The simulated radial temperature profiles at different reactor heights and under different operating conditions are compared with the experimental measurements by Valstar et al. and their prediction33 as shown in Figure 3. It is seen that the simulated temperature profiles by the present model are closer to the measurements and better than the simulation by Valstar et al. The temperature profile of the whole reactor for case 1 is displayed in Figure 4. The Axial Relative Velocity Profile. The simulated axial velocity profile along the radial direction is shown in Figure 5 in the form of relative velocity that refers to gas flow. The profile demonstrates that the axial velocity falls quickly to zero at the wall surface due to the no slip boundary, reaches the maximum at the near-wall region where the porosity is highest and oscillates in a strongly damped fashion toward the center of the fixed bed. The axial velocity displays relatively less fluctuation at about 2-3dp from the wall. The shape of the axial velocity distribution is in agreement with published experimental investigations.18,20,44,45 The Turbulent Mass-Transfer Diffusivity Profile. As a feature of applying the present model, the turbulent diffusivity for mass transfer Di,t and the concentration distribution can be obtained. Taking case 1 for example, the calculated DAC,t profile

along the height of reactor is shown in part a of Figure 6, in which the DAC,t profile is seen to reach the fully developed condition after the reactant traveling a short distance about 20de from the entrance. The variation of diffusivity along the radial direction is shown in part b of Figure 6, where the DAC,t is seen to be nearly constant only at the central region, then reaches to a maximum at some distance near the wall region, and sharply drops down to zero at the wall surface. All of these trends are attributed to the nonuniform distributions of velocity, temperature, and concentration. The shape of radial DAC,t profiles at different axial positions are consistent with the research results28 and the experiments24,25,38 by using gas trace technique. It is also shown from parts a and c of Figure 6 that the diffusivity profiles of DAC,t and DHAc,t along the reactor are in close similarity, as they are in the same velocity and temperature fields although their concentrations are a little different. The Turbulent Thermal Diffusivity Rt Profile. The turbulent thermal diffusivity Rt can be calculated according to the accompanied t2-t model. Still, taking case 1 for example, the calculated Rt profiles are plotted in Figure 7, in which, like the diffusivity Di,t, Rt quickly reaches to the fully developed state after flowing about 50de in distance from the reactor entrance (part a of Figure 7) and sharply decreased to the minimum at some distance near the wall region (part b of Figure 7). The volume-weighted average value of Rt based on the superficial velocity for case 1 was found to be on the order of magnitude of about 1.0 × 10-4 m2/s by the present simulation. Olbrick and Potter46 found experimentally that the fluid phase effective radial Peclet number in the fixed-bed was about 7-8 under the condition of R/dp > 5, which means that the Rt for the present case is about 1.14 × 10-4 m2/s. Therefore, Rt predicted by the proposed model is on the same order of magnitude as the experiment. The Effective Thermal Conductivity Profile. The effective thermal conductivity profiles at different axial positions are shown in Figure 8. It can be seen that keff sharply drops down at the near-wall region, which is in agreement with the experimental observation by Kwong and Smith.37 The volumeweighted average effective thermal conductivity keff computed by the proposed model in case 1 is about 0.21 J m-1 s-1 K-1, which is less than that of 0.43 J m-1 s-1 K-1 as reported by Valstar et al.33 obtained by using the plug-flow model. Nevertheless, according to the study on catalytic reactor,20 the effective thermal conductivity obtained by using more accurate flow model, which took into account the uneven porosity profile, turned to be about 20-40% lower than that using the plugflow model. Furthermore, they also found that keff by the

2662

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

Figure 6. (a,b) Turbulent mass-transfer diffusivity profiles DAC,t and (c) DHAc,t of the reactor (case 1).

Figure 7. Turbulent thermal diffusivity profiles Rt along the reactor (case 1).

Figure 8. Effective thermal conductivity profiles along the radial direction at different axial positions (case 1).

more accurate flow model did not show a length effect, which is consistent with the prediction by the proposed model (Figure 9). The Turbulent Kinematic Viscosity νt Profile. As a result of the present simulation, the νt can be obtained. Taking case 1 for example, the νt profile along the height of reactor is shown in Figure 10, in which νt is fully developed after the reactant entering the reactor with uneven distribution along the radial

Figure 9. Effective thermal conductivity profile along axial direction (case 1).

direction and falls rapidly to zero at the wall surface. As seen from part b of Figure 10, part b of Figure 6, and part b of Figure 7, the radial νt, DAC,t, and Rt along the radial direction at different

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2663

Figure 10. (a) Turbulent kinematic viscosity νt profile in the whole reactor, (b) turbulent kinematic viscosity νt profiles along the radial direction (case 1).

reactor heights are varying, and the local ratios of νt/Dt and νt/Rt, or the local Sct and Prt numbers, are changing in the whole reactor. Furthermore, according to the presented model, the computed volume-averaged DAC,t and Rt are 2.8 × 10-3 m2/s and 3.59 × 10-4 m2/s, respectively, which is different from the calculated values of DAC,t ) 5.13 × 10-4 m2/s and Rt )4.22 × 10-4 m2/s by simply using an empirical Sct number (usually take Sct ) νt/Dt ) 0.7) and Prt number (usually take Prt ) νt/Rt ) 0.85) together with volume-averaged turbulent kinematic viscosity νt ) 4.48 × 10-4 m2/s by CFD simulation. It demonstrates that the conventional method of using a fixed Sct number and Prt number for modeling the turbulent mass and heat transfer process is only a rough approximation. Conclusions (1) A computational transport model is proposed, which consists of computational fluid dynamics, computational heat transfer, and the recently developed computational mass-transfer models for simulating a wall-cooled catalytic reactor. The proposed model enables us to obtain the radial and axial distributions of velocity, temperature, and concentration at once without assuming certain empirical Prt and Sct numbers or using the experimental coefficient of dispersion. A feature of the proposed model is introducing the c2-c model to close the turbulent mass-transfer equation to obtain the unknown turbulent mass-transfer diffusivity. The predicted velocity, temperature, and concentration distributions are in agreement with the experimental data reported from the literature. (2) The simulation results show that the turbulent diffusivity of mass transfer, the turbulent thermal diffusivity of heat transfer, and the turbulent kinematic viscosity of momentum transfer are changing and nonuniformly distributed in the reactor, although their distributions show somewhat similarity. (3) From the profiles of νt, DAC,t, and Rt along the radial direction at different reactor heights, it demonstrates that the local ratios of νt/Dt and νt/Rt, or the local Sct and Prt numbers, are varying throughout the whole reactor. Acknowledgment The authors acknowledge the financial support by the National Natural Science Foundation of China (contract no. 20736005) and the assistance from the staff in the State Key Laboratories for Chemical Engineering (Tianjin University).

Notation c2 ) fluctuating concentration variance Ci ) average component concentration, mass fraction Cµ, c1, c2 ) constant in k- model equations Cc0, Cc1, Cc2, Cc3 ) constant in c2-c model equations CD0, CD1, CD2, CD3, CD4 ) constant in t2-t model equations cp ) specific heat, J kg-1 K-1 de ) effective diameter of catalyst, m dp ) nominal diameter of catalyst, m Di,eff ) component effective diffusivity for mass transfer, m2 s-1 Di,t ) component turbulent diffusivity for mass transfer, m2 s-1 G ) gas-phase flow rate per unit cross-section area, kg m-2 s-1 ∆Hr ) heat of reaction, kJ mol-1 k ) turbulent kinetic energy, m2 s-2 keff ) effective thermal conductivity, J m-1 K-1 s-1 k0e ) effective thermal conductivity of stagnant fixed beds, J m-1 K-1 s-1 kf ) fluid thermal conductivity, J m-1 K-1 s-1 kL ) molecular thermal conductivity of fluid, J m-1 K-1 s-1 ks ) thermal conductivity of catalyst, J m-1 K-1 s-1 kt ) turbulent conductivity of fluid based on superficial velocity, J m-1 K-1 s-1 L ) reactor length, m Mi ) molecular weight of component, kg kmol-1 Pe ) Peclet number Prt ) turbulent Prandtl number pi ) partial pressure of component, atm Q ) thermal source term of temperature equation, J m-3 s-1 R ) position in radial direction, m R ) radius of the column, m R0 ) the resistant coefficient of porous media Rg ) gas constant, kJ kmol-1 K-1 Rs ) apparent reaction rate, kmol kg-1 (cat) s-1 Sct ) turbulent Schmidt number Si ) source term, kg m-3 s-1 t0 ) inlet temperature, °C t2 ) temperature variance, K2 T ) temperature, K u ) interstitial velocity vector, m s-1 U ) superficial velocity, m s-1 x ) axial position, m

2664

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008

z ) dimensionless distance, z ) (R - r)/de Zn(Ac)2 ) zinc acetate Greek Letters R, Rt ) molecular and turbulent thermal diffusivities, m2 s-1  ) turbulent dissipation rate, m2 s-3 i,c ) turbulent dissipation rate of fluctuating concentration variance, s-1 t ) turbulent dissipation rate of fluctuating temperature variance, s-1 Φ ) variable γ ) porosity distribution of the random packing bed. γ∞ ) porosity in an unbounded packing µ, µt, µeff ) molecular, turbulent and effective viscosity, kg m-1 s-1 F ) density, kg m-3 Fb ) bulk density of catalyst, kg m-3 νt ) turbulent kinematic viscosity, νt )µt/F, m2 s-1 σc, σc ) model parameters in c2-c model equations σt ) model parameter in t2-t model equations σk, σ ) model parameters in k- model equations Subscripts AC ) acetylene eff ) effective f ) fluid HAc ) acetic acid i ) AC, HAc, VA or 1, 2, 3 in ) inlet L ) laminar s ) solid t ) turbulent VA ) vinyl acetate w ) wall Literature Cited (1) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design; John Wiley & Sons, Inc.: New York, 1990. (2) Wakao, N.; kaguei, S. Heat and Mass Transfer in Packed Beds; Gordon and Breach Science Publishers: New York, 1982. (3) Tarhan, M.O. Catalytic Reactor Design; McGraw-Hill: New York, 1983. (4) Winterberg, M.; Tsotsas, E.; Krischke, A.; Vortmeyer, D. A Simple and Coherent Set of Coefficients for Modelling of Heat and Mass Transport with and without Chemical Reaction in Tubes Filled with Spheres. Chem. Eng. Sci. 2000, 55 (5), 967-979. (5) Krishnaswamy, S.; Gunn, R. D.; Agarwal, P. K. Low-Temperature Oxidation of Coal. 2. An Experimental and Modelling Investigation Using a Fixed-Bed Isothermal Flow Reactor. Fuel 1996, 75 (3), 344-352. (6) Hong, R.; Li, X.; Li, H.; Yuan, W. Modeling and Simulation of SO2 Oxidation in a Fixed-Bed Reactor with Periodic Flow Reversal. Catal. Today 1997, 38 (1), 47-58. (7) Roblee, L. H. S.; Baird, R. M. Tierney, J. W. Radial Porosity Variations in Packed Beds. AIChE. J. 1958, 4 (4), 460-464. (8) Benenati, R. F.; Brosilow, C. B. Void Fraction Distribution in Beds of Spheres. AIChE. J. 1962, 8 (3), 359-361. (9) Scott, G. D. Radial Distribution of the Random Close Packing of Equal Spheres. Nature 1962, 194, 956-957. (10) Ridgway, K.; Tarbuck, K. Radial Voidage Variation in Randomly Packed Beds of Spheres of Different Sizes. J. Pharm. Pharmacol. 1966, 18, 168-175. (11) Thadani, M. C.; Peebles, F. N. Variation of Local Void Fraction in Randomly Packed Beds of Equal Spheres. Ind. Eng. Chem. Proc. Des. DeV. 1966, 5 (3), 265-268. (12) Ridgway, K.; Tarbuck, K. J. Voidage Fluctuations in RandomlyPacked Beds of Spheres Adjacent to a Containing Wall. Chem. Eng. Sci. 1968, 23 (9), 1147-1155.

(13) Goodling, J. S.; Vachon, R. I.; Stelpflug, W. S.; Ying, S. J.; Khader, M. S. Radial Porosity Distribution in Cylindrical Beds Packed with Spheres. Powder Technol. 1983, 35 (1), 23-29. (14) Toye, D.; Marchot, P.; Crine, M.; Pelsser, A. M.; L’Homme, G. Local Measurements of Void Fraction and Liquid Holdup in Packed Columns Using X-ray Computed Tomography. Chem. Eng. Process. 1998, 37 (6), 511-520. (15) Kufner, R.; Hofmann, H. Implementation of Radial Porosity and Velocity Distribution in a Reactor Model for Heterogeneous Catalytic Gas Phase Reactions (TORUS-Model). Chem. Eng. Sci. 1990, 45 (8), 21412146. (16) Mueller, G. E. Radial Void Fraction Distributions in Randomly Packed Fixed Beds of Uniformly Sized Spheres in Cylindrical Containers. Powder Technol. 1992, 72 (3), 269-275. (17) de Klerk, A. Voidage Variation in Packed Beds at Small Column to Particle Diameter Ratio. AIChE. J. 2003, 49 (8), 2022-2029. (18) Lerou, J. J.; Froment, G. F. Velocity, Temperature and Conversion Profiles in Fixed-Bed Catalytic Reactors. Chem. Eng. Sci. 1977, 32 (8), 853-861. (19) Kalthoff, O.; Vortmeyer, D. Ignition/Extinction Phenomena in a Wall Cooled Fixed Bed Reactor: Experiments and Model Calculations Including Radial Porosity and Velocity Distributions. Chem. Eng. Sci. 1980, 35 (7), 1637-1643. (20) Daszkowski, T.; Eigenberger, G. A Reevaluation of Fluid Flow, Heat Transfer and Chemical Reaction in Catalyst Filled Tubes. Chem. Eng. Sci. 1992, 47 (9-11), 2245-2250. (21) Papageorgiou, J. N.; Froment, G. F. Simulation Models Accounting for Radial Voidage Profiles in Fixed-Bed Reactors. Chem. Eng. Sci. 1995, 50 (19), 3043-3056. (22) Nijemeisland, M.; Dixon, A. G. Comparison of CFD Simulations to Experiment for Convective Heat Transfer in a Gas-Solid Fixed Bed. Chem. Eng. J. 2001, 82 (1-3), 231-246. (23) Dixon, A. G.; Nijemeisland, M.; Stitt, E. H. Packed Tubular Reactor Modeling and Catalyst Design Using Computational Fluid Dynamics. In AdVances in Chemical Engineering; Marin, G. B., Ed.; Academic Press, Elsevier: New York, 2006; Vol. 31. (24) Fahien, R. W.; Smith, J. M. Mass Transfer in Packed Beds. AIChE. J. 1955, 1 (1), 28-37. (25) Dorweiler, V. P.; Fahien, R. W. Mass Transfer at Low Flow Rates in a Packed Column. AIChE. J. 1959, 5 (2), 139-144. (26) Sun, Z. M.; Liu, B. T.; Yuan, X. G.; Liu, C. J.; Yu, K. T. New Turbulent Model for Computational Mass Transfer and its applicatIon to a Commercial-Scale Distillation Column. Ind. Eng. Chem. Res. 2005, 44 (12), 4427-4434. (27) Liu, G. B.; Yu, K. T.; Yuan, X. G.; Liu, C. J. New Model for Turbulent Mass Transfer and its Application to the Simulations of a PilotScale Randomly Packed Column for CO2-NaOH Chemical Absorption. Ind. Eng. Chem. Res. 2006, 45 (9), 3220-3229. (28) Liu, G. B.; Yu, K. T.; Yuan, X. G.; Liu, C. J.; Guo, Q. C. Simulations of Chemical Absorption in Pilot-Scale and Industrial-Scale Packed Columns by Computational Mass Transfer. Chem. Eng. Sci. 2006, 61 (19), 6511-6529. (29) Liu, B. T. Study of a New Mass Transfer Model of CFD and Its Application on Distillation Tray, Ph.D. Thesis, Tianjin University, Tianjin, China, 2003. (30) Yin, F. H.; Sun, C. G.; Afacan, A.; Nandakumar, K.; Chuang, K. T. CFD Modeling of Mass-Transfer Processes in Randomly Packed Distillation Columns. Ind. Eng. Chem. Res. 2000, 39, 1369-1380. (31) Guo, B. Y.; Yu, A. B.; Wright, B.; Zulli, P. Comparison of SeVeral Turbulence Models Applied to the Simulation of Gas Flow in a Packed Bed, Third International Conference on CFD in the Minerals and Process Industries, Melbourne, Australia, 2003; pp 509-514. (32) Jones, W. P.; Musonge, P. Closure of the Reynolds Stress and Scalar Flux Equations. Phys. Fluids 1988, 31 (12), 3589-3604. (33) Valstar, J. M.; Berg, P. J.; Oyserman, J. Comparison between Two Dimensional Fixed Bed Reactor Calculations and Measurements. Chem. Eng. Sci. 1975, 30 (7), 723-728. (34) Zhang, Y. S. Turbulence; National Defense Industry Press: Beijing, China, 2002. (35) Yin, F. H.; Sun, C.; Afacan, G. A.; Nandakumar, K.; Chuang, K. T. CFD Modeling of Mass-Transfer Processes in Randomly Packed Distillation Columns. Ind. Eng. Chem. Res. 2000, 39, 1369-1380. (36) Launder, B. E.; Spalding, D. B. Lectures in Mathematical Models of Turbulence; Academic Press: London, 1972. (37) Kwong, S. S.; Smith, J. M. Radial Heat Transfer in Packed Beds. Ind. Eng. Chem. 1957, 49 (5), 894-903.

Ind. Eng. Chem. Res., Vol. 47, No. 8, 2008 2665 (38) Schertz, W. W.; Bischoff, K. B. Thermal and Material Transport in Non-Isothermal Packed Beds. AIChE. J. 1969, 15 (4), 597-604. (39) Nield, D. A. Estimation of the Stagnant Thermal Conductivity of Saturated Porous Media. Int. J. of Heat Mass Transfer 1991, 34, 15751576. (40) Nijemeisland, M.; Dixon, A. G. CFD Study of Fluid Flow and Wall Heat Transfer in a Fixed Bed of Spheres. AIChE. J. 2004, 50 (5), 906921. (41) Bejan, A.; Kraus, A. D. Heat Transfer Handbook; John Wiley & Sons, Inc.: New Jersey, 2003. (42) Yaws, C. L. Yaws’ Handbook of Thermodynamic and Physical Properties of Chemical Compounds; Knovel: New York, 2003. (43) Tsotsas, E.; Schlu¨nder, E. U. Heat Transfer in Packed Beds with Fluid Flow: Remarks on the Meaning and the Calculation of a Heat Transfer Coefficient at the Wall. Chem. Eng. Sci. 1990, 45 (4), 819-837.

(44) Maximo, M.; Spinn, C. W.; Smith, J. M. Velocities and Effective Thermal Conductivities in Packed Beds. Ind. Eng. Chem. 1951, 43 (1), 225232. (45) Giese, M.; Rottscha¨fer, K.; Vortmeyer, D. Measured and Modeled Superficial Flow Profiles in Packed Beds with Liquid Flow. AIChE J. 1998, 44 (2), 484-490. (46) Olbrich, W. E.; Potter, O. E. Heat Transfer in Small Diameter Packed Beds. Chem. Eng. Sci. 1972, 27 (9), 1723-1732.

ReceiVed for reView May 24, 2007 ReVised manuscript receiVed January 6, 2008 Accepted January 30, 2008 IE070737Y