Modeling of a Tubular Reactor Producing Epichlorohydrin with

Sep 30, 2010 - and pertinent kinetic parameters are estimated by a genetic algorithm combined with the Levenberg-Marquardt algorithm. The time-varying...
0 downloads 0 Views 3MB Size
Download PDF
Ind. Eng. Chem. Res. 2011, 50, 1187–1195

1187

Modeling of a Tubular Reactor Producing Epichlorohydrin with Consideration of Reaction Kinetics and Deactivation of Titanium Silicate-1 Catalyst Woohyun Kim,† Choamun Yun,†,⊥ Young Kim,‡ Jeongho Park,† Sunwon Park,*,† Ki Taek Jung,§ Yong Hwa Lee,§ and Sae Heon Kim§

Downloaded via KAOHSIUNG MEDICAL UNIV on November 10, 2018 at 23:00:13 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Department of Chemical & Biomolecular Engineering, Korea AdVanced Institute of Science and Technology, 335 Gwahak-ro, Yuseong-gu, Daejeon, 305-701, South Korea, Korea Institute of Machinery and Materials, 104 Sinseongno, Yuseong-gu, Daejeon, 305-343, South Korea, and Hanwha Chemical R&D Center, 6 Shinseong-dong, Yuseong-gu, Daejeon, 305-804, South Korea

A mathematical model for a tubular reactor producing epichlorohydrin (ECH) over titanium silicate-1 (TS-1) catalyst is proposed in this work. The Eley-Rideal mechanism is selected to describe the catalytic reaction, and pertinent kinetic parameters are estimated by a genetic algorithm combined with the Levenberg-Marquardt algorithm. The time-varying catalytic activity and effectiveness factor of a commercial TS-1 catalyst pellet are additionally included in the two-dimensional dynamic mass and energy transport model for the packedbed tubular reactor. In a case study, the optimal length of the tubular reactor is determined by simulations using the developed model to maintain both conversion and selectivity above 96% during 90 days. The proposed model will contribute to effective design of industrial ECH production processes. 1. Introduction Epichlorohydrin (ECH) is a major raw material for epoxy resins and synthetic glycerin. With steady growth in the demand for epoxy resins for electronics, automobiles, and windmill generators, the production of ECH has increased significantly. The conventional process for industrial production of ECH, which originates from the 1950s, generates a large amount of wastewater containing lime and chlorinated organic compounds such as allyl chloride, hypochlorite, and dichlorohydrines.1-3 Although several alternative processes have been developed for the production of ECH, the efficiency with these approaches is still low or the supply of raw material is unstable.4-10 However, recent studies on epoxidation catalyzed by titanium silicate (TS1) catalyst in aqueous methanol indicate the possibility of commercializing a new ECH production process.3,10-15 This process shows high conversion and selectivity for ECH and emits little wastewater.3 [TS-1] + MeOH a [TS-1 · MeOH]

2. Reaction Kinetics As the first step, a mechanistic kinetic model is developed for the reaction of allyl chloride and hydrogen peroxide over TS-1 catalyst. In previous studies, the catalytic reaction mechanism, shown in reactions R1-R4, was determined to follow the Eley-Rideal mechanism. The rate law corresponding to this mechanism is as follows:13-15 Rmain )

k0 dCECH kk1CHPCALC ) MTS-1 dt 1 + k0 1 + k1CHP + k2CALC

(1)

(R1)

k ) Ae-Ea/RT

(R2)

In addition to the main reaction synthesizing ECH, two major side reactions consume ECH, as depicted in Figure 2. The rates

[TS-1 · MeOH] + H2O2 a [TS-1 · MeOH · H2O2] [TS-1] + allyl chloride a [TS-1 · allyl chloride]

characteristics of industrial TS-1 catalysts is established. Finally, the mass and energy transport models of the reactor are simulated to identify the optimal reactor length for long-term operation of a reactor producing ECH above the desired quality level.

(R3)

[TS-1 · MeOH · H2O2] + allyl chloride a [TS-1 · MeOH] + epichlorohydrin + water (R4) In this work, a dynamic model of a packed-bed tubular reactor is developed to facilitate commercialization of the ECH production process over TS-1 by determining the optimal design of the reactor. The following steps, depicted in Figure 1, explain the hierarchy of this work. First, the kinetics of relevant reactions is modeled, and second, a dynamic reactor model reflecting the * To whom correspondence should be addressed. Tel.: 82-42-3503920. Fax: 82-42-862-5961. E-mail: [email protected]. † Korea Advanced Institute of Science and Technology. ‡ Korea Institute of Machinery and Materials. § Hanwha Chemical R&D Center. ⊥ Present address: Doosan Heavy Industries & Construction, Daejeon, South Korea.

Figure 1. Hierarchy of modeling, simulation, and optimization.

10.1021/ie101071m  2011 American Chemical Society Published on Web 09/30/2010

(2)

1188

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

η)

actual overall rate of reaction rate of reaction resulting if entire interior surface × exposed to external pellet surface conditions

(7)

Figure 2. Side reactions consuming epichlorohydrin.

of the side reactions are described in power-law kinetic models rather than in detailed mechanistic kinetics, because the mechanisms of the side reactions have not yet been exactly verified and these rates are very slow compared to the main reaction. For side reaction 1, depicted in Figure 2, we assume that a small amount of hydrogen chloride from allyl chloride is the only chlorine source for the production of 1,3-dichloro-2propanol (1,3DC2PO). It is also assumed that hydrogen peroxide acts as a catalyst. Therefore, the rate of side reaction 1 depends on temperature and the concentrations of ECH and allyl chloride, respectively. On the other hand, the rate of side reaction 2 (see Figure 2), producing 1-chloro-3-methoxy-2-propanol (1C3M2PO), is a function of temperature and the concentration of ECH; note that the amount of methanol does not affect the reaction rate, because methanol is assumed to always be abundant. For the kinetic modeling of the side reactions, it is also assumed that both 1,3DC2PO and 1C3M2PO are exclusively produced by the two side reactions. r1,3DC2PO )

dC1,3DC2PO ) ks1CECHR1CALCβ1 dt ks1 ) A1e-Ea1/RT

For practical reasons, for powdered catalysts it is often assumed that their entire active sites are exposed to bulk fluid flow. Thus, the effectiveness factor of the powder is 1 while that of a pellet ranges from 0 to 1. The rate law of the main reaction is directly multiplied by the effectiveness factor in order to reflect the decrease of the reaction rate caused by diffusion limitations. The effectiveness factor is assumed to be applied to the main reaction only although side reaction 2 (see Figure 2) is a catalytic reaction. Because the rate of side reaction 2 is normally around 1/100 of the main reaction rate, the limitation of the reaction rate caused by diffusion resistance will be negligible. Rmain

Rmain ) Φ

dC1C3M2PO ) ks2CECHR2 MTS-1 dt

ks2 ) A2e-Ea2/RT

(5)

(6)

The decomposition of hydrogen peroxide releasing oxygen might be considered as another side reaction. However, the decomposition of hydrogen peroxide is negligible below 40 °C, where the rate of the decomposition is very low.16 In this 4-mlong packed-bed tubular reactor, the hot spots, representing the hottest temperatures, are managed below 40 °C (or 35 °C at some hot spots), while the space time of the process flow is less than 20 s. Furthermore, laboratory-scale experimental results showed that most of the hydrogen peroxide is consumed in the main reaction. 3. Reactor Modeling 3.1. Effectiveness Factor. If the diffusion of reactants to the surface of a catalyst pellet is not sufficiently fast, the concentration of reactants near the surface starts to decrease so rapidly that the reaction rate is limited as a result of the low concentration of reactants. Accordingly, the actual overall rate of reaction is not as fast as it could be if the interior surface of the pellet were exposed to the external pellet surface conditions. The decrease in the overall rate of reaction is measured and represented by the effectiveness factor of the pellet.17

k0 ηAe-Ea/RTk1CHPCALC 1 + k0 1 + k1CHP + k2CALC

(9)

where Φ)

r1C3M2PO )

(8)

3.2. Catalyst Activity Function. The catalytic activity of solid catalysts in heterogeneous catalysis normally decays with time. The catalyst activity function (Φ), which describes the freshness of the catalyst, is employed to present the actual rate of reactions. While the function can be either separable or nonseparable from the rate law, most of the reactions following the Eley-Rideal mechanism are known to involve a separable catalytic activity model.18

(3)

(4)

k0 ηAe-Ea/RTk1CHPCALC ) 1 + k0 1 + k1CHP + k2CALC

current conversion conversion without catalyst deactivation

When the catalyst is fresh, Φ ) 1, and the value decreases from 1 to 0 as the catalyst is used for increasingly long time. The retarded rate of the reaction due to catalyst deactivation is predicted by the catalyst activity function. Among several chemical and mechanical factors inducing catalyst deactivation,17,19,20 fouling and coking by bulky byproducts, are known to be the major factors for the reaction of hydrogen peroxide over titanium silicate catalysts.21 Thermal degradation of the catalyst is assumed to be negligible, as the reactions occur at mild temperature below 50 °C. Therefore, the catalyst activity model for this reaction can be defined as a function of the rate of coke formation. As the rate of coking on the active sites is not readily measurable, the catalyst activity function is often defined as a function of time or the cumulative quantity of the product.22 3.3. Two-Dimensional Transport Phenomenon Model. The rates of reaction and catalyst deactivation are functions of temperature and concentrations in the reactor. To predict and optimize the overall reactor performance, two-dimensional mass and energy models are developed to describe the transport phenomena in the reactor described in Figure 3. The mass balance model of the unit reactor we have developed includes reaction kinetic models, the effectiveness factor, and the catalyst activity function. The model also considers dispersion combining the two-dimensional diffusion equation and the convective term in order to describe mixing driven by a concentration gradient in the flow.17,23-25

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

1189

(iii) r ) 0 (m), ∂T/∂r ) 0; no radial heat transfer at the center of the reactor tube (iv) r ) R (m), ke ∂T/∂rr)R ) kw(Two(z) - T(z,R)) ) Qw, Qw ) ho(Tc(z) - Two(z)) To describe heat transfer between the reactor and the cooling jacket, Qw, we also assumed that heat flux to the inner wall of the tube is the same as that through the wall, since the thickness of the wall is merely 2 mm and the tube is made of steel, which has high thermal conductivity. At the same time, the heat flux through the wall is the same as that from the outer surface of the tube to the coolant. The energy balance for the coolant is expressed below. Figure 3. Configuration of the 1-m-long unit reactor for modeling twodimensional transport phenomena.

-Qwa - FcCp,c

mass balance: De ∂Ci ∂Ci ∂2Ci ∂2Ci ∂Ci + D e 2 - Uz + Ri ) De 2 + r ∂r ∂z ∂t ∂r ∂z

a)

(10) Since the first, second, and third terms of eq 10 represent all the mass transport phenomena driven by a concentration gradient, the model can describe dispersion in the flow. The effective diffusivity, De, is introduced to present the average diffusivity to all directions, since the precise diffusion paths in a packed-bed reactor are not easily predictable.17 Diφpσc De ) τ˜

energy balance: ke ∂ r ∂T ∂2T + ke 2 + ∆HRx,iRi r ∂r ∂r ∂z ∂Ci ∂T + UzCi Cp )( -De ∂z ∂z

[ ( )]

(

)

(13)

4 πdL ) di πdi2L/4

(14)

The coolant is assumed to have a negligible temperature gradient in the radial direction considering that the flow rate and the velocity of the coolant are large and sufficiently fast. For the calculation of the heat flux from the outer wall of the tube to the coolant, the convective heat transfer coefficient is given as follows.29

( ) ( )

hodo doG ) 0.287 kc µc

(11)

The boundary conditions for the mass balance model are defined as follows: (i) z ) 0 (m), Ci(z,r) ) Ci,0 (ii) z ) 1 (m), ∂Ci/∂z ) 0; no axial diffusion at the end of the reactor tube (iii) r ) 0 and R (m), ∂Ci/∂r ) 0; no radial diffusion at the center and the inner wall of the reactor tube In the proposed model, channeling and back mixing are assumed to be negligible as the model satisfies the following criteria suggested in previous works:26,27 (i) ratio of catalyst bed height to catalyst pellet size (L/dp) g 50 (ii) ratio of inner diameter of the reactor to catalyst pellet size (di/dp) g 10 L/dp and di/dp of the unit reactor are 100 and 10.2, respectively.

∂Tc ∂Tc ) CwCp,c ∂V ∂t

0.61

Cp,cµc kc

0.33

(15)

Fa

4. Parameter Estimation To estimate the kinetic parameters of each reaction, objective functions to be minimized are given as the mean square error of the prediction of the reaction rate. For the main reaction minimize: f(K) ) N k0 Ae-Ea/RTk1CHP,iCALC,i 1 Ri,main,expt N i)1 1 + k0 1 + k1CHP,i + k2CALC,i

∑

(

(

))

2

(16)

where K ) [k0, k1, k2, A, Ea]T ) kinetic parameters N ) number of experimental data For side reaction 1 (see Figure 2) minimize: f(K) ) N 1 (R - (A1e-Ea1/RTCECH,iR1CALC,iβ1))2 (17) N i)1 i,side 1,expt

∑

∑ C )C i

p

∂T ∂t

(12)

For the energy balance model, we assume that all the heat flux caused by the difference of temperature is described as conduction with effective thermal conductivity, ke.23,28 This value is estimated by measuring the temperature change in the radial direction. The temperature change along the radial direction can be measured by varying the depth of a thermometer at a single location of the tubular reactor. For the unsteadystate energy balance model, the boundary conditions are defined as follows: (i) z ) 0 (m), T(z,r) ) Ti (ii) z ) 1 (m), ∂T/∂z ) 0; no axial heat transfer at the end of the reactor tube

where K ) [A1, Ea1, R1, β1]T ) kinetic parameters N ) number of experimental data For side reaction 2 (see Figure 2) minimize: f(K) )

1 N

N

∑ (R

i,side 2,expt

- (A2e-Ea2/RTCECH,iR2))2

i)1

(18) where K ) [A2, Ea2, R2]T ) kinetic parameters N ) number of experimental data The main reaction rate of the ith experiment, Ri,main,expt, is measured in laboratory-scale experiments for various concentrations of the reactants and temperatures. For the side reactions

1190

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

Figure 4. Application of a novel optimization technique combining a genetic algorithm and the Levenberg-Marquardt algorithm.

(see Figure 2), it should be noted that both 1,3-dichloro-2propanol (1,3DC2PO) and 1-chloro-3-methoxy-2-propanol (1C3M2PO) are generated in the presence of TS-1; in contrast, only 1,3DC2PO is produced without TS-1 catalyst. Therefore, kinetic parameters for side reaction 1 (see Figure 2), i.e., the production of 1,3DC2PO, are determined first from experimental results obtained under the condition that TS-1 catalyst is not used. Subsequently, the kinetic parameters for side reaction 2 (see Figure 2), i.e., the production of 1C3M2PO, are estimated from experimental results with TS-1 considering that both side reactions occur at the same time. Equations 16-18 are nonlinear multimodal problems that can be solved by optimization methods. In this work, a novel optimization technique that combines a genetic algorithm (GA) and the Levenberg-Marquardt algorithm (LMA) is proposed to estimate the parameters in complex kinetic models. A GA is a heuristic algorithm based on the theory of evolution in a natural system. The algorithm works with a set of solutions called a population, and each solution is called an individual. After an initial population is created, the objective function value for each individual is computed. Subsequently, individuals with superior objective function values are selected to generate the next generation of the population through genetic operators: crossover and mutation of the selected individuals. In general applications of a GA, this regeneration procedure is repeated until the solution satisfies the termination criteria.30,31 Although GAs have been applied to kinetic parameter estimation for various catalytic reactions,32-34 they suffer from the drawback that the computational load dramatically increases for more complex problems. When handling multidimensional problems, a significantly larger population is required, and the overall number of calculations exponentially grows with an increase of the population. In the proposed approach, this issue is resolved by locally employing a fast optimization method, the LMA.35 The LMA is mathematically derived by a linear approximation, and the searching method is based on calculation of a gradient.36-38

f(K) ≈ f(Ki) + ∇Tf(K)∆Ki

(19)

where ∆Ki ) Ki+1 - Ki, and Ki is the kinetic parameter of the ith iteration. To find the minimum of eq 19, the equation is differentiated with respect to ∆Ki and the resulting equation equals zero as follows. ∇f(K) ) ∇f(Ki) + H(Ki)∆Ki ) 0

(20)

∆Ki ) Ki+1 - Ki ) -[H(Ki)]-1∇f(Ki)

(21)

or

where H(Ki) is the Hessian matrix of f(x) in eqs 20 and 21. The search direction for the minimum of the objective function is defined as eq 21 and a solution is found by iterative calculations until ∆Ki is less than a specified tolerance. However, if the Hessian matrix of f(x) is a negative-definite matrix, the solution cannot converge to the optimum. The LMA modifies the Hessian matrix at each stage of the search to be positivedefinite, and the negative-definite Hessian matrix is replaced ˜ (Ki), in eq 22. with a modified Hessian matrix, H ˜ (Ki) ) [H(Ki) + βI] H

(22)

where I is a unit matrix and β is a positive constant, which is ˜ (Ki) is positive-definite. sufficiently large to ensure that H The initial point of the LMA determines the rate of convergence, since the number of iterations of eqs 21 and 22 depends on this variable. Therefore, the selection of the initial point is the most critical factor for the LMA when the solution space is large. The overall procedure of the combined method is described in Figure 4. A large solution space is searched by a GA without concern for the initial points. After a given number of generations, the GA provides the initial point of the LMA with the best solution found thus far. Since the initial point is already close to the optimum, the LMA converges rapidly to find the

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

1191

Figure 5. Absolute percent errors of predicted main reaction rates using kinetic parameters estimated by the GA combined with the LMA (blue bars), Newton’s method (red bars), and the conjugate gradient method (green bars). Table 1. Comparison of Prediction Results Estimated by the GA Combined with the LMA, Newton’s Method, and Conjugate Gradient Method

average absolute percent error [%] root-mean-square error [(f(K))1/2]

a

main reaction side reaction 1a side reaction 2a main reaction [mol/min kg of catal] side reaction 1a [mol/m3 min] side reaction 2a [mol/min kg of catal]

GA combined with LMA

Newton’s method

conjugate gradient method

8.725 11.082 15.941 0.2229 0.2004 0.0013

21.364 19.328 17.543 0.6960 0.3023 0.0013

21.842 27.863 17.543 0.5611 0.4092 0.0013

See Figure 2.

Figure 6. Parameter estimation for the catalyst activity function of the TS-1 catalyst.

optimal solution. The results presented in Figure 5 and Table 1 indicate that the combined method consistently shows small errors, whereas other deterministic methods show inferior performance. After estimation of all the kinetic parameters, the effectiveness factor of the TS-1 pellets for the main reaction is determined. To determine the effectiveness factor of the TS-1 pellets, the

rate law of the main reaction without consideration of the effectiveness factor, eq 5, is applied to a model for the 1-mlong unit reactor of this process. The reaction rates calculated by simulation of the model are then compared with actual pilotscale experimental data. The ratio of the calculated values to the experimental data becomes the effectiveness factor of the TS-1 pellets used in this work.

1192

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

Figure 7. Configuration of the pilot-scale reactor.

As the last step of the parameter estimation, a catalyst activity model is defined. A catalyst activity model is usually an empirical equation and can be shown in many different types of functions. In this work, the model is derived as a function of the cumulative quantity of ECH production, as shown in eq 23. The three parameters of the model are determined by fitting the graph of the activity change of the 1-m-long unit reactor versus the cumulative quantity of produced ECH, as illustrated in Figure 6. Φ(CUMECH) )

a + b · CUMECH 1 + c · CUMECH

(23)

where a, b, and c are constants. 5. Application The pilot-scale reactor is composed of four unit reactors, as illustrated in Figure 7. Several thermal transmitters are installed to measure the temperature changes along the axial direction at the radial centers of the reactor. The pilot-scale reactor is charged with a mixture of methanol/water solvent, allyl chloride, and hydrogen peroxide, where the allyl chloride to hydrogen peroxide ratio is 4 to 1 on a molar basis. The temperature of each cooling jacket is independently adjusted to enhance the reactor performance. The reactor model is built and simulated in Aspen Custom Modeler.39 All the physical properties of the reactants and the products are calculated with the UNIQUAC method in Aspen Properties.40 A central finite difference method is selected to solve the differential equations of the developed model. Simulation of the two-dimensional reactor model generates various profiles of conversion, selectivity, and temperature. In this section, first, the reactor model is validated with pilotscale experimental data. Subsequently, the optimal length of the pilot-scale reactor is determined by dynamic simulations to maintain the reactor performance at or above the desired level during long-term operation. 5.1. Pilot-Scale Experiments and Dynamic Simulations of the Model. In the first set of pilot-scale experiments, the temperature of the feed is maintained at 35 °C and the flow rate is 4.84 L/h. The reactor is operated for 80 h. The temperature of each cooling jacket is dynamically adjusted to maintain the highest temperature of each unit reactor at roughly 40 °C during operation of the reactor. Up to 20 h of operation, it is observed that the measured temperatures at the hot spot of the first unit reactor are slightly higher than the simulated values. However, the temperature profiles of the simulation results are quite similar to the measured temperatures at most of the points, as shown in Figure 8. In the two-dimensional temperature profiles obtained by simulation, shown in Figure 9, the locations and the movements of the hot spots can also be predicted.

Figure 8. Center temperatures of the pilot-scale reactor after 20 and 60 h.

The experimental data of conversion and selectivity are obtained at each end of the first, second, and fourth unit reactors, respectively. In the dynamic simulation results, the conversion gradually decreases while the selectivity does not vary as much as the conversion does. This tendency is also observed in the pilot-scale experiments. Table 2 shows that the absolute percent errors of the simulation results are generally less than 5%. In the second set of pilot-scale experiments, the temperatures of the feed and the coolant in the first unit reactor are changed and the reactor is operated for 15 h. By changing the operating conditions, the changes in the reactor performance are examined and the accuracy of the developed model, which should reflect the dynamic behaviors of the pilot-scale reactor, can be evaluated. In this case, the feed temperature is maintained at 30 °C and the temperature at the hot spot of the first reactor is controlled to be lower than 35 °C. The temperature profiles obtained by dynamic simulations are shown in Figure 10. It should be noted that the differences between the measured and the simulated temperature profiles are generally less than 2 °C. Table 3 indicates that the conversion and the selectivity calculated by the simulations are close to the experimental data. The absolute percent errors of the simulation results are within 10%. For selectivity, in particular, all the percent errors are less than 1%. 5.2. Case Study: Determining the Reactor Length. When designing a catalytic tubular reactor, there are two major

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

1193

Figure 9. Two-dimensional temperature profile after 20 h. Table 2. Simulation Errors of the Developed Modela

Table 3. Comparison between Experimental Data of the Second Pilot-Scale Experiment and Simulated Valuesa

b

error (%) 2m 4m operation time conversion selectivity conversion selectivity conversion selectivity

errorb (%)

1m

15 min 20 h 70 h 80 h

-4.21 -5.24 9.09 -4.04

0.43 0.32 0.47 0.46

15.72 1.62 -

0.11 0.94 -

5.51 1.55 1.35

0.23 0.73 0.88

a -, not measured. b Error (%) ) 100[(simulated value - measured value)/measured value].

1m 2m 4m operation time conversion selectivity conversion selectivity conversion selectivity 2h 4h 10 h 15 h

-1.49

-0.18

9.01 3.94 0.51 -5.08

-0.26 0.07 0.23 -0.07

5.61 5.34 5.29 4.77

0.97 -0.48 -0.03 -0.32

a -, not measured. b Error (%) ) 100[(simulated value - measured value)/measured value].

Figure 11. Predicted performance change of the pilot-scale reactor during 90-day operation.

Figure 10. Center temperature of the pilot-scale reactor after 2 and 14 h.

concerns: catalyst deactivation and the trade-off relationship between conversion and selectivity. Due to the catalyst deactivation in the reactor, decrease of the conversion is inevitable during long-term operation, as shown in Figure 11. A longer tubular reactor could raise the conversion level. However, a longer tubular reactor also increases the residence time of ECH and, accordingly, the side reactions consuming ECH can be more stimulated. Moreover, from an economics point of view, the length of the tubular reactor should be minimized, because it is the most important factor regarding the investment cost of an industrial-scale reactor.

In this case study, the objective is to identify the minimum number of serially connected unit reactors, where both conversion and selectivity should be maintained at or above 0.96 for 90 days. The flow rate and the composition of the feed are set to be the same as those in the first pilot-scale experiment. The overall procedure is described in Figure 12. Initially, the change of the reactor performance during 90-day operation is predicted using the model developed for the 4-mlong tubular reactor. If the conversion and the selectivity do not satisfy the desired levels, one unit reactor is added and the model is modified to reflect the corresponding changes in the design. A dynamic simulation of the updated model is then executed again to evaluate the reactor performance. This procedure is iteratively followed until the minimum number of unit reactors to achieve the desired conversion and selectivity is found. The results of the case study are presented in Table 4. The minimum length of the reactor is determined within several iterations.

1194

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011

Figure 12. Iterative procedure to determine the optimal length of the reactor. Table 4. Results of the Case Study number of unit reactors reactor performance after 90 days

5

6

7

8

conversion selectivity

0.853 0.964

0.906 0.962

0.942 0.960

0.961 0.961

The model will significantly contribute to the development of an industrial-scale reactor by providing required information for design and scale-up of the reactor. Acknowledgment This work was supported by Hanwha Chemical R&D Center.

6. Conclusions A mathematical model for a tubular reactor producing ECH has been developed. The model includes the mechanistic reaction kinetics as well as the characteristics of the TS-1 catalyst, i.e., the deactivation and the effectiveness factor of the catalyst pellet. The predicted temperatures, conversion, and selectivity deviate only slightly from the experimental results, and the model is capable of reflecting dynamic behaviors of the reactor under different operating conditions. Furthermore, the design problem of the catalytic tubular reactor has been discussed in relation to a case study. After iterative dynamic simulations, the optimal length of the reactor was determined such that both conversion and selectivity are maintained at or above the desired levels during long-term operation.

Nomenclature A, A1, and A2 ) preexponential factors CALC ) concentration of allyl chloride CECH ) concentration of epichlorohydrin CHP ) concentration of hydrogen peroxide Ci ) concentration of component i Cp ) heat capacity of mixture Cp,c ) heat capacity of coolant Cw ) number of moles of coolant per unit volume C1,3DC2PO ) concentration of 1,3-dichloro-2-propanol C1C3M2PO ) concentration of 1-chloro-3-methoxy-2-propanol CUMECH ) cumulative moles of produced epichlorohydrin De ) effective diffusivity Di ) diffusivity of component i

Ind. Eng. Chem. Res., Vol. 50, No. 3, 2011 di ) tube inner diameter do ) tube outer diameter Ea, Ea1, and Ea2 ) activation energy Fa ) arrangement factor Fc ) coolant flow rate G ) coolant mass flux k0, k1, and k2 ) kinetic parameters ho ) convective heat transfer coefficient ∆HRx,i ) Heat of reaction kc ) conductivity of coolant ke ) effective conductivity kw ) conductivity of tube R ) radius of tube R ) ideal gas constant r ) radial direction Ri ) rate of formation of the consumption of component i Rmain ) rate of main reaction r1,3DC2PO ) rate law for formation of 1,3-dichloro-2-propanol r1C3M2PO ) rate law for formation of 1-chloro-3-methoxy-2propanopl T ) temperature of process flow Tc ) coolant temperature Ti ) feed temperature Two ) temperature at the outer surface of the tube [TS-1] ) mass of catalyst Uz ) velocity of fluid in axial direction Qw ) heat flux to outer surface of tube z ) axial direction Greek Symbols R1, β1, and R2 ) kinetic parameters η ) effectiveness factor Φ ) catalyst activity function φp ) porosity τ˜ ) tortuosity σ ) constriction factor µc ) coolant viscosity

Literature Cited (1) Nagato, N.; Mori, H.; Maki, K.; Ishioka, R. (Kawaski). U.S. Patent 4,634,784, 1987. (2) Bijsterbosch, J. W.; Das, A.; Kerkhof, F. P. J. M. Clean Technology in the Production of Epichlorohydrin. J. Cleaner Prod. 1994, 2, 181. (3) Pandey, R. K.; Kumar, R. Eco-friendly Synthesis of Epichlorohydrin Catalyzed by Titanium Silicate (TS-1) Molecular Sieve and Hydrogen Peroxide. Catal. Commun. 2007, 8, 379. (4) Milas, N. A.; Sussman, S. The Hydroxylation of the Double Bond. J. Am. Chem. Soc. 1936, 58, 1302. (5) Milas, N. A.; Sussman, S. The Hydroxylation of Unsaturated Substances. IV. The Catalytic Hydroxylation of Unsaturated Hydrocarbons. J. Am. Chem. Soc. 1937, 59, 2345. (6) Milas, N. A.; Sussman, S.; Mason, H. S. The Hydroxylation of Unsaturated Substances. V. The Catalytic Hydroxylation of Certain Unsaturated Substances with Functional Groups. J. Am. Chem. Soc. 1939, 61, 1844. (7) Milas, N. A.; Maloney, L. S. The Hydroxylation of Unsaturated Substances. VI. The Catalytic Hydroxylation of Cyclopentadiene. J. Am. Chem. Soc. 1940, 62, 1841. (8) Payne, G. B.; Smith, C. W. Reactions of Hydrogen Peroxide. III. Tungstic Acid Catalyzed Hydroxylation of Cyclohexene in Nonaqueous Media. J. Org. Chem. 1957, 22, 1682. (9) Sheldon, R. A. Synthetic and Mechanistic Aspects of Metal-catalysed Epoxidations with Hydroperoxides. J. Mol. Catal. 1980, 7, 107. (10) Clerici, M. G.; Romano, U. U.S. Patent 4,824,976, 1989. (11) Neri, C.; Anfossi, B.; Esposito, A.; Buonomo, F. U.S. Patent 4,833,260, 1989. (12) Escrig, P. D. F.; Martin, J. M. C. U.S. Patent 6,160,138, 2000.

1195

(13) Clerici, M. G.; Bellussi, G.; Romano, U. Synthesis of Propylene Oxide from Propylene and Hydrogen Peroxide Catalyzed by Titanium Silicate. J. Catal. 1991, 129, 159. (14) Bellussi, G.; Carati, A.; Clerici, M. G.; Maddinelli, G.; Millini, R. Reactions of Titanium Silicate with Protic Molecules and Hydrogen Peroxide. J. Catal. 1992, 133, 220. (15) Gao, H.; Lu, G.; Suo, J.; Li, S. Epoxidation of Allyl Chloride with Hydrogen Peroxide Catalyzed by Titanium Silicate 1. Appl. Catal., A: Gen. 1996, 138, 27. (16) Chen, L. Y.; Chuah, G. K.; Jaenicke, S. Propylene Epoxidation with Hydrogen Peroxide Catalyzed by Molecular Sieves Containing Framework Titanium. J. Mol. Catal. A: Chem. 1998, 132, 281. (17) Fogler, H. S. Elements of Chemical Reaction Engineering, 3rd ed.; Prentice Hall PTR: Upper Saddle River, NJ, 1999. (18) Lynch, D. T.; Emig, G. On the Separability of Catalyst Activity and Kinetic Behavior. Chem. Eng. Sci. 1989, 44, 1275. (19) Bartholomew, C. H. Mechanisms of Catalyst Deactivation. Appl. Catal., A: Gen. 2001, 212, 17. (20) Bartholomew, C. H.; Farrauto, R. J. Fundamentals of Industrial Catalytic Processes, 2nd ed.; Wiley: Hoboken, NJ, 2006. (21) Wang, Q.; Wang, L.; Chen, J.; Wu, Y.; Mi, Z. Deactivation and Regeneration of Titanium Silicate Catalyst for Epoxidation of Propylene. J. Mol. Catal. A: Chem. 2007, 273, 73. (22) Froment, G. F. Modeling of Catalyst Deactivation. Appl. Catal., A: Gen. 2001, 212, 117. (23) Peters, P. E.; Schiffino, R. S.; Harriott, P. Heat Transfer in Packedtube Reactors. Ind. Eng. Chem. Res. 1988, 27, 226. (24) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena, 2nd ed.; Wiley: New York, 2002. (25) Froment, G. F.; Bischoff, K. B. Chemical Reactor Analysis and Design, 2nd ed.; Wiley: New York, 1990. (26) Akpan, E.; Akande, A.; Aboudheir, A.; Ibrahim, H.; Idem, R. Experimental, Kinetic and 2-D Reactor Modeling for Simulation of the Production of Hydrogen by the Catalytic Reforming of Concentrated Crude Ethanol (CRCCE) over a Ni-based Commercial Catalyst in a Packed-bed Tubular reactor. Chem. Eng. Sci. 2007, 62, 3112. (27) Rase, H. F. Chemical Reactor Design for Process Plants; Wiley: New York, 1977. (28) Bahrami, M.; Yovanovich, M. M.; Culham, J. R. Effective Thermal Conductivity of Rough Spherical Packed Beds. Int. J. Heat Mass Transfer 2006, 49, 3691. (29) Stultz, S. C.; Kitto, J. B. SteamsIts Generation and Use, 40th ed.; Babcock and Wilcox: New York, 1992. (30) Holland, J. H.; Adaptation in Natural and Artificial Systems; University of Michigan Press: Ann Arbor, MI, 1975. (31) Goldberg, D. E. Genetic Algorithms in Search, Optimization, and Machine Learning, 1st ed.; Addison-Wesley Professional: Reading, MA, 1989. (32) Daneshpayeh, M.; Khodadadi, A.; Mostoufi, N.; Mortazavi, Y.; Sotudeh-Gharebagh, R.; Talebizadeh, A. Kinetic Modeling of Oxidative Coupling of Methane over Mn/Na2WO4/SiO2 Catalyst. Fuel Process. Technol. 2009, 90, 403. (33) Wolf, D.; Moros, R. Estimating Rate Constants of Heterogeneous Catalytic Reactions without Supposition of Rate Determining Surface Steps an Application of a Genetic Algorithm. Chem. Eng. Sci. 1997, 52, 1189. (34) Wolf, D.; Buyevskaya, O. V.; Baerns, M. An Evolutionary Approach in the Combinatorial Selection and Optimization of Catalytic Materials. Appl. Catal., A: Gen. 2000, 200, 63. (35) Venugopal, G.; Suryakant; Balaji, C.; Venkateshan, S. P. A Hybrid Optimization Technique for Developing Heat Transfer Correlations Based on Transient Experiments. Int. J. Heat Mass Transfer 2009, 52, 1954. (36) Levenberg, K. A Method for the Solution of Certain Nonlinear Problems in Least Squares. Q. Appl. Math. 1944, 2, 164. (37) Marquardt, D. W. An Algorithm for Least-squares Estimation of Nonlinear Parameters. J. Soc. Indust. Appl. Math. 1963, 11, 431. (38) Edgar, T. F.; Himmelblau, D. M.; Lasdon, L. S. Optimization of Chemical Process, 2nd ed.; McGraw-Hill: New York, 2001. (39) Web site for Aspen Custom Modeler: http://www.aspentech.com/ products/aspen-custom-modeler.cfm. (40) Web site for Aspen Properties: http://www.aspentech.com/products/ aspen-properties.cfm.

ReceiVed for reView May 10, 2010 ReVised manuscript receiVed August 26, 2010 Accepted September 7, 2010 IE101071M