Reactor Modeling Part 4 (Petroleum Refining)

When a high-purity hydrogen stream without gas recycle is used, such as in the case of some laboratory and bench-scale HDT reactors, or when the gas recycle has been subject to the purification process in commercial units, values of partial pressure (pG) and liquid molar concentrations (CL) of H-S, NH3, and LHC at the entrance of the catalytic bed (z = 0 and z = LB for co – current and countercurrent operation, respectively) are equal or very close to zero. For commercial HDT reactors without the high purification of a gas recycle stream, values of partial pressures (pG) and liquid molar concentrations (CL) of H2S, NH3, and LHC at the entrance of the catalytic bed (z = 0) differ from 0 (Mederos et al., 2006).

The axial and radial dispersion terms of mass and heat result in a second-order differential equation for all phases; consequently, two boundary conditions are necessary. According to Danckwerts (Wehner and Wilhelm, 1956), the Danckwerts’ boundary condition at z = 0 is

which can be simplified to

These boundary conditions are true because the axial dispersion of mass and heat is relatively small and the concentration and temperature gradients at the reactor inlet are quite flat (Chen et al., 2001). The Danckwerts’ boundary condition at z = LB is is the real condition at the outlet of the reactor, but as it is an infinite type of boundary condition, it is difficult to apply under most numerical methods.

The boundary condition

TABLE 2.14. Boundary Conditions (t > 0) of Generalized Mass and Heat Balance Equations

		Gas Phase		Liquid Phase
			No		No
Condition	Operation Mode	Dispersion	Dispersion	Dispersion	Dispersion
	Co – current
	Countercurrent
	Co – current				-
	Countercurrent
	Co – current Countercurrent
	Co – current/ countercurrent
	Co – current/ countercurrent

Because of this, many researchers prefer to use Danckwerts’ boundary conditions, where the temperature gradients at the exit of the reactor are zero. It must be pointed out that unless the reactor is infinitely long, this exit condition will not be true because the fluid will not reach equilibrium with its surroundings. This lack of equilibrium at the exit would be even more pronounced at the high temperatures frequently experienced in TBRs. The boundary conditions proposed by Young and Finlayson (1973) avoid this abnormality. They derived inlet and outlet conditions similar to those of Danckwerts for the case in which both axial and radial dispersion is present. However, the outlet condition was reported with a nonzero gradient, which is reduced to Danckwerts form when radial dispersion is neglected. The boundary condition dTS / dz = 0 for z > Lb isa Danckwerts type also, but as the solid phase is commonly nonexistent for z > LB, it is generally regarded as a valid exit condition. For the same reason, as the reaction does not take place without the solid catalyst, the concentration gradients can be assumed to be zero at the exit. Therefore, the Danckwerts boundary conditions are used when there is no mass (or heat) dispersion outside the reactor (Chao and Chang, 1987- Feyo De Azevedo et al., 1990).

			Solid Phase
Dispersion	No Dispersion	Surface (i=All Compounds)	Interior (i = All Compounds)

The boundary condition atis sometimes considered to be

(Shokri and Zarrinpashne, 2006; Shokri et al., 2007). According to Chen et al. (2001), the following boundary conditions will be employed for heat transfer in the thermowell:

The general energy heat balance for the liquid phase is

Assuming no temperature gradient within the catalyst particles (isothermal catalyst), the general heat balance for the solid phase is

The resulting set of PDEs coupled with the respective initial and boundary conditions are then solved simultaneously using an appropriate numerical method (Melis et al., 2004 ).

Example of Simplification of the Generalized Model Sometimes, to simplify heat transfer modeling in HDT reactors, the three processes involved (heat transfer in the solid, liquid, and gas phases) can be lumped into only one equation with a pseudohomogeneous heat balance. In the following, an example of how to obtain a pseudohomogeneous heat balance from the generalized heat balance presented in Table 2.12 is shown.

The general heat balance for the gas phase is When modeling commercial HDT reactors it is generally accepted that axial dispersion can be omitted; then HA3 = HB3 = HC3 = 0. If adiabatic operation is also assumed, then HA4 = HB4 = HC4 = 0 (isothermal reactor in the radial direction) and HA8 = HB8 = 0 (no heat transfer between the fluid phase and the reactor wall).

The pseudohomogeneous model is based on the fact that the temperature difference among the gas, liquid, and catalyst at any particular axial position of the reactor is negligible. Hence, TG = TI = TL = TSS = T, and in consequence the following terms can be neglected: HA5 = HA7 = HB5 = HB7 = HC9 = HC10 = 0. The temperature gradients TG – TI and TI – TL in the HA6 and HB6 terms are also neglected. The final heat balance equation for the gas, liquid, and solid phases is then

Equation (2.137) in terms of model parameters for co-current operation gives

The final simplified pseudohomogeneous heat balance is

If the catalyst wetting efficiency is also assumed to be complete, fw = 1 (e.g., when properly designed liquid distributors and high liquid velocities in commercial reactors are used). Therefore, reactions occur only in the liquid phase; then Eq. (2.139) is rewritten as

According to Feyo De Azevedo et al. (1990), some researchers have included the radiation effects in nonisothermal reactors in the radial direction (HA4 = HB4 = HC4 ^ 0) by means of the effective radial thermal conductivity

whereis the conductive plus convective contributions to the effective radial conductivity,the radiant contribution,the Stephan-Boltzmann constant, anda radiant transfer factor defined as

where e is the particle emissivity. Therefore, rearranging Eq. (2.151) as

, the pseudohomogeneous radial heat dispersion term (HA4 + HB4 + HC4 withis expressed as

where the last term on the right side of Eq. (2.144) represents the radioactive heat transfer, and its exclusion may produce errors in estimation of some other heat transfer parameters.

Estimation of Model Parameters

To solve the set of ordinary differential equations (ODEs) (for the steady-state regime) or the set of PDEs (for the dynamic regime), it is necessary to evaluate several parameters and chemical properties of the system. Those parameters can be estimated with existing correlations, whose accuracy is of great importance for the entire state of robustness of the reactor model.

Effective Diffusivity The description of steady-state diffusion and reaction of a multicomponent liquid mixture in a porous catalyst particle requires appropriate definition of the species fluxes to the particle. Fick’s law for equi-molar counterdiffusion through an ideal cylindrical pore is given by

where Df is the effective diffusivity, by means of which the structure (porosity and tortuosity) of the pore network inside the particle is taken into consideration in the modeling. Table 2.15 shows Bosanquet’s formula to estimate the effective diffusivity inside the catalyst particle (Bosanquet, 1944), which consists of two diffusion contributions: Knudsen diffusivity (DfK) and molecular diffusivity (DfMi). Both Bosanquet’s formula and Fick’s law can also be applied with sufficient accuracy to cases involving a narrow unimodal pore-size distribution and very dilute mixtures. The restrictive factor F(Xg) accounts for additional friction between the solute and the pore walls. The exponent Z is 4 for Xg < 0.2 (Iliuta et al., 2006). The tortuosity factor of the pore network, t, is used in the calculation of Dfei because the pores are not oriented along the normal direction from the surface to the center of the catalyst particle, and its value generally varies between 3 and 7, but for HDT process it is commonly assumed to be equal to 4. It is also possible to estimate the tortuosity factor assuming valid the upper bound correlation given by Weissberg (1963) for a packing of random spheres as shown in Table 2.15.

Effectiveness Factor The effectiveness factor of independent reactions can be defined as the ratio of the volumetric average of the reaction rate into the particle to the reaction rate at the surface of the particle as proposed by Thiele (1939) and Zeldovich (1939) :

Analytical solutions for Eqs. (2.146) and (2.147) are possible only for single reactions and for zero- and first-order rate expressions. The various correlations used in the literature to estimate the catalyst effectiveness factor for isothermal and irreversible reactions are shown in Table 2.16.

For kinetic models other than the power-law approach, such as Langmuir-Hinshelwood-Hougen-Watson (LHHW)-type kinetic expressions, there is no analytical solution of Eq. (2.146) – Therefore, an alternative method to avoid the numerical integration of Eq. (2.146) is the Bischoff generalized modulus approach (Bischoff, 1965), which enables an analytical solution to any type of rate equation and single reaction.

TABLE 2.15. Estimation of Effective Diffusivity

Parameter	Gas Phase Liquid Phase
Molecular diffusivity coefficient
Knudsen diffusivity coefficient
Effective diffusivity
Restrictive factor
Ratio of radius of gyration over pore radius
Mean pore radius
Tortuosity factor
Binary diffusion coefficient (P in atm)
Dynamic liquid viscosity
Molar volume of solute (i) in liquid phase and liquid solvent (L)
Solvent critical specific volume

As mentioned previously, the effectiveness factor in commercial HDS catalysts has been reported to be in the range 0.4 to 0.8.Expressions for isothermal first-order reactions with irregularly shaped catalysts lead under steady-state conditions to acceptable results with errors not exceeding 20% (Aris, 1975; Dudukovic, 1977). Bischoff (1965) proposed a general modulus to predict the effectiveness factor for any reaction type within a relatively narrow region. If reactions of order less than one-half are excluded, the spread between all the various curves is about 15%. The mean deviation of values of n calculated from the empirical correlation proposed by Papayannakos and Georgiou (1988) is less than 2.4% from those predicted with the normalized modulus for simple order reactions proposed by Froment and Bischoff (1990) in the range 0.05 < n < 0.99.

TABLE 2.16. Estimation of Catalyst Effectiveness Factor

Kinetic Model	Reaction Order	Shape	Thiele Modulus	Effectiveness Factor
Power law		Spheres and crushed
		Any geometry
		Any geometry
		Any geometry
LHHW		Spheres and crushed
		Pellet, cylinder, 2-, 3 – lobe, etc.

Global Gas-Liquid Mass Transfer The gas-liquid interphase mass transfer flux is described in terms of the simple two-film theory:

The overall external resistance to mass transfer (KLi) is composed by the resistance to mass transfer in the gas ( kG) and liquid ( kt) films. Estimates of the compressibility factor Z sometimes give values close to 1; therefore, ideal gas law could be used (Mejdell et al., 2001).

For slightly soluble gases such as H2 , the value of Henry’ s constant (Hi ) exceeds unity, and then mass transfer resistance in the gas film can be neglected (Zhukova et al., 1990). Therefore, the total mass transfer is approximately equal to the liquid-side mass transfer coefficient:

The liquid film mass transfer coefficient ( k;L) is calculated using the correlations reported in Table 2.17 .

Gas-Liquid Equilibrium The gas-liquid equilibrium along the catalyst bed is represented in the mass balance equations by the Henry’s law constant. The constants related to this law for different chemical species available in the system may be defined in two ways, described below:

1. Solubility coefficients. Employing solubility coefficients, the following expression is used to estimate Henry’s constant:

where stands for the component – solubility and vN is the molar volume under normal conditions. Using this expression implies knowledge of the gaseous component solubility in the liquid phase considering the process temperature effect. Korsten and Hoffmann (1996) have reported the next correlations to evaluate this parameter only for H2 and H2S: for hydrogen sulfide. However, this last correlation may not be no adequate to evaluate the solubility of H2S in oil fractions at the complete temperature range used in commercial units. In that case, it is possible to use an EoS to estimate this Henry’s constant.

TABLE 2.17. Correlations to Estimate Model Parameters

Parameter	Symbol	References	Parameter	Symbol	References
Holdup gas		Calculated from £b = £l + £g	Density of the liquid		Ahmed (1989)
			phase
Holdup liquid		Charpentier and Favier (1975), Satterfield (1969),	Catalyst bulk density		ASTM (2003)
		Specchia and Baldi (1977), Ellman et al. (1990)
Bed void fraction		Haughey and Beveridge (1969), Carberry and	Dynamic viscosity of		Glaso (1980) , Ahmed (1989) ,
		Varma (1987), Froment and Bischoff (1990)	liquid		Brule and Starling (1984)
Catalyst wetting efficiency (or		Ring and Missen (1991), Al-Dahhan and	Dynamic viscosity of		Ahmed (1989), Brule and Starling
contacting effectiveness)		Dudukovic (1995)	gas		(1984)
Two – phase pressure drop		Larkins et al. (1961) , Ellman et al. (1988)	Diffusion coefficients		Wilke and Chang (1955)
Gas-liquid interfacial area		Iliuta et al. (1999)	for gases Diffusion coefficients		Tyn and Calus (1975)
			for liquids
Gas-solid interfacial area		Puranik and Vogelpohl (1974), Onda et al. (1967)	Specific heat of gas		Lee and Kesler (1975), Perry et al.
			phase		(2004)
Mass transfer gas-solid		Petrovic and Thodos (1968) , Dwivedi and	Specific heat of		Lee and Kesler (1975, 1976),
coefficient		Upadhyay (1977)	liquid phase		Perry et al. (2004)
Mass transfer liquid-solid		Dudukovic et al. (2002), Dwivedi and Upadhyay	Specific heat of solid		Perry et al. (2004)
coefficient		(1977), Bird et al. (2002), Evans and Gerald	phase
		(1953) , Wilson and Geankoplis (1966) , Goto
		and Smith (1975), Satterfield et al. (1978),
		Specchia et al. (1974)
Mass transfer coefficient of		Goto and Smith (1975)	Gas-liquid heat		Marroqum de la Rosa et al.
liquid side at G-L interface			transfer coefficient		(2002), Chilton and Colburn
					(1939)
Mass transfer coefficient of		Goto and Smith (1975), Reiss (1967),Yaici et al.	Liquid-solid heat		Chilton and Colburn (1939)
gas side at G-L interface		(1988)	transfer coefficient
Axial dispersion of gas		Hochman and Effron (1969), Sater and	Chilton-Colburn		Froment and Bischoff (1990) , Hill
		Levenspiel (1966) , Demaria and White (1960)	‘-factor for energy		(1977), Bird et al. (2002), Gupta
			transfer		et al. (1974)
Axial dispersion of liquid		Gierman (1988) , Hochman and Effron (1969) ,	Effective thermal		Hashimoto et al. (1976)
		Sater and Levenspiel (1966), Tsamatsoulis and	conductivity radial
		Papayannakos (1998)
Radial dispersion of gas		Fahien and Smith (1955)	Effective thermal		Tarhan (1983) , Dixon (1985)
			conductivity axial
Radial dispersion of liquid		Fahien and Smith (1955), De Ligny (1970),	Thermal conductivity		API (1997) , Chung et al. (1988)
		Herskowitz and Smith (1978a,b)	of f phase
Binary interaction parameter		Moysan et al. (1983), Ronze et al. (2002), Riazi	Heat of reaction		Tarhan (1983)
for H2-oil using PR EoS		(2005), Lal et al. (1999), Magoulas and Tassios
		(1990)
Binary interaction parameter		Feng and Mather (1993a,b), Carroll and Mather	Heat of vaporization/		Soave (1972) , Peng and Robinson
for H2S – oil using PR EoS		(1995)	condensation		(1976)
Binary interaction parameter		API (1997)	Thiele modulus		Bischoff (1965)
for NH3- oil using SRK EoS
Density of the gas phase		Soave (1972) , Peng and Robinson (1976)	Effectiveness factor		Froment and Bischoff (1990) , Aris (1975)

 

for hydrogen and

2. Equation ofstate. It is also possible to obtain Henry-s constant assuming local equilibrium at the liquid-gas interface:

where the equilibrium constant (Ke q,i ) is calculated using an adequate EoS (e.g., Peng-Robinson (PR), Soave-Redlich-Kwong (SRK), t-van der Waals, Grayson-Streed). For this expression, it was necessary to calculate a two-phase thermodynamic equilibrium previously at each local point along the catalytic bed to estimate the interface temperature and the molar compositions in liquid and gas phases. The main advantage of this expression is that it takes into account the volatility of the feedstock; however, it also increases the computing time too greatly.

Another way to calculate the Henry’s constant of gaseous solute in a solvent is to use the next thermodynamic assumption:

where yt is the fugacity coefficient of a gaseous compound i (solute) in the liquid phase (solvent), and its calculation using EoS is addressed below. The SRK and PR equations of state are the most widely used in HDT process modeling, being defined by the generalized expression

For pure compounds, values of parameters a and b are given by

The generalized temperature function a (Tf) was proposed by Soave (1972) to be an equation of the form with

with

When applied to mixtures, the classical mixing rules may be considered to evaluate parameters a and b:

where ki k are the binary interaction parameters, which may be obtained from the references reported in Table 2.17. It is important to point out that the quality of Henry’s constant calculation depends enormously on the accuracy of these interaction parameters. The liquid-phase fugacity coefficient can be derived from Eq. (2.155) to give

with

TABLE 2.18. Parameters of Soave-Redlich-Kwong and Peng-Robinson Equations of State

Parameter	SRK	EoS PR
	1
	0
	0.42748
	0.08664
	0.48
	1.574
	-0.176
	1
	-1
	-AB

where ZL is the compressibility factor of the liquid phase at saturation obtained from the solution of Eq. (2.155) expressed in its cubic compressibility factor form:

To evaluate the liquid-solid mass transfer coefficients, for example, Eq. (2.168) must be expressed as

The values of the universal parameters and DZ) are given in Table 2.18 .

Heat Transfer Coefficients Correlations employed for mass transfer can be used to calculate the parameters for energy transfer between phases by employing the Chilton and Colburn (1939) analogy. The Chilton-Colburn ‘-factor for mass transfer (jD) is given by

where g(Re/) is a correlation that is a function of the Reynolds number of the respective phase f to be evaluated; some of these correlations are shown in Table 2.17 – If one needs to estimate the liquid-solid heat transfer coefficient, for example, Eq. (2.169) is rewritten in the expression

There are various correlations to estimate the mass transfer coefficients, but the Chilton-Colburn analogy is not usually employed to evaluate them, whereas due to the lack of correlations to estimate the heat transfer coefficients in the gas or liquid film side at the gas-liquid interface and in the liquid film at the liquid-solid interface, it is common to use the Chilton-Colburn analogy, by equating Eqs. (2.168) and (2.170) (jD = jH), to estimate these coefficients. The physical and geometrical properties involved in the dimensionless numbers must be evaluated at conditions of reaction and for each phase in a heterogeneous reactor.

To use the Chilton-Colburn analogy, it is necessary to consider the following conditions (Bird et al., 2002):

• Constant physical properties

• Small net mass transfer rates

• No chemical reaction

• No viscous dissipation heating

• No absorption or emission of radian energy

• No pressure diffusion, thermal diffusion, or forced diffusion

Theoretical Calculations of Some Parameters Relative to the Catalyst

Bed Dilution of the catalyst bed with inert material is a common practice in experimental HDT reactors (Sie, 1996). The following simple formula is employed to calculate the dilution factor:

where Vc is the catalyst volume and V is the volume of inert particles, both obtained experimentally.

An equivalent particle diameter (dpe), defined as the diameter of a sphere that has the same external surface (or volume) as the actual catalyst particle,

On the other hand, the Chilton-Colburn /-factor for energy transfer (jH) is given by

is an important particle characteristic that depends on particle size and shape. For fixed beds with catalyst extrudates of commercial size, the equivalent particle diameter can be calculated according to Cooper et al. (1986):

Bed void fraction (or bed porosity) for undiluted catalyst bed can be calculated with the following correlation reported by Froment and Bischoff (1990), Haughey and Beveridge (1969), and Carberry and Varma (1987):

This correlation was developed for undiluted packed beds of spheres; however, if the equivalent particle diameter concept is used, it can also be employed for nonsphere particles. Once the bed void fraction is determined, the particle density can be calculated as follows (Tarhan, 1983):

Since the continuous models are based on the volume-average form of the transport equations for multiphase systems, equations expressing conservation of volume are (Whitaker, 1973)

Relationships between phase holdups inside the catalyst solid are given in the following expressions:

The external surface area of catalyst particles per unit of reactor volume for PBRs can be calculated as

Catalyst porosity (eS) may be calculated with the following equation from the experimental data for total pore volume (Vg):

In an extreme case where the experimental Sg parameter is not available in order to estimate the average pore radius, one can use the correlation proposed by Mace and Wei (1991):

where

Some parameters that account for bed characterization are experimentally measurable, others are experimental or can be obtained through simulations, and others are empirical. Of course, although it is better to obtain the local porosity experimentally, measurements required the use of advanced techniques. To do that, computational calculations are preferred.

Most of empirical correlations for predicting liquid saturation, pressure drop, and flow regimes are based on experiments performed at atmospheric conditions. Since industrial trickle-bed reactors are operated at high pressures and temperatures, the applicability of these correlations for such operating conditions needs to be investigated (Nguyen et al., 2006). The majority of correlations that have been developed on a laboratory scale may not work for large-scale reactors (operated at high pressure and temperature) due to significant changes in hydrodynamic characteristics with a reactor scale (Gunjal and Ranade, 2007). Although extensive studies in hydrodynamic correlations are available in the literature (Dudukovic et al., 2002) it seems that many works in modeling reactors for petroleum fractions still employ the classical correlations (i.e., those derived from reasonable assumptions with simple expressions). Some researchers have employed the same correlations of previous papers without checking the accuracy of these expressions or the range in which they are applicable. On the other hand, although the benefit of using correlations based on neural networks has been reported, many data are necessary to employ this approach. To overcome this feature, an alternative is to create a database with different types of crude oils and fractions in order to develop an online program able to determine the various parameters involved in modeling a HDT reactor. At present it is not available. The only valuable effort seems to be that of Larachi et al. (1999), who have developed a simulator based on neural networks for the prediction of some hydrodynamic parameters for trickle-bed reactors. Although neural networks are updated continuously, which favors its use for parameter predictions, it is necessary to develop a fundamental relationship that takes advantage of novel techniques for characterization of hydrodynamic parameters.

The limitations of the various models reported in the technical literature are closely related to the number of parameters involved and to the reliability of the data available. Therefore, appropriate correlations for mass and energy transfer should be employed in order to calculate each of the terms used in a model reactor (i.e., correlations developed under similar conditions, such as the same flow regime, pressure, and temperature, assuming similar liquid system and porous particles). Generalized correlations for the prediction of properties having a broad range of variation should be avoided when modeling a TBR in detail, because they could produce some miscalculations. Constant values assumed a priori, such as tortuosity factor, binary interaction parameters, heat of reaction, and specific heat, should be used as a reference when experimental data are available. To simplify a TBR model, some researchers have ignored the low heat of some reactions because its contribution is not significant in the energy balance, which seems to be a reasonable assumption.

Using different approaches to fluid dynamics, kinetics or thermodynamics can lead to very different conclusions in predictions of reactor performance. For example, Gunjal and Ranade (2007) have reported 15% more conversions considering all parameters to be constant and assuming only uniform porosity (i.e., no nonuniform porosity). Akgerman et al. (1985) has reported 24 to 38% higher conversions when considering volatiles with respect to nonvolatiles, and Inoue et al. (2000) has predicted the use of larger reactors employing an nth kinetic model when more accurate kinetic expression has been utilized. Another important finding is the consideration of a thermowell in the energy balance for a pilot plant, as pointed out by Chen et al. (2001). These are the reasons to account for detailed models which allow both for making accurate descriptions of the chemical phenomena and for reliable preliminary calculations when designing a TBR reactor.