Ind. Eng. Chem. Res. 1988,27, 2015-2022
2015
Simulation and Optimization of an Industrial Ammonia Reactor Said S. Elnashaie,* Mohamed E. Abashar, and Abdulrahman S. Al-Ubaid Chemical Reaction Engineering Group (CREG),Chemical Engineering Department, College of Engineering, King Saud University, P.O. Box 800, Riyadh 11421, Saudi Arabia
The heterogeneous model for ammonia converter is used to simulate an industrial reactor having three adiabatic beds with interstage cooling. The two-point boundary value differential equation for the catalyst particles is solved by using the orthogonal collocation method. The model is used to investigate the optimal behavior of the reactor, and comparison between different approaches to the optimal operation of the reactor is carried out. I t is found that the heterogeneous model with optimization based on the actual rate of reaction gives the best performance, with an improvement of 6 7 % in % NH, at the exit compared with the suboptimal three-bed reactor with interstage cooling. The temperature-constrained optimization of the reactor is also investigated in order to limit the temperature to the limits that the catalyst can withstand, and it was found that the appreciable improvement in % NH3 achieved with unconstrained optimization is still maintained. Mathematical models for simulation and optimization studies have been built up by many workers for various types of ammonia synthesis reactors (Annable, 1952; Van Heerden, 1953; Kjaer, 1963; Baddour et al., 1965; Shah, 1967; Gaines, 1977; Singh and Saraf, 1979; Reddy and Husain, 1982; Mahfouz, 1985; Elnashaie et al., 1988; Mansson and Andresen, 1986). In this work, a heterogeneous model is developed for a reador having three adiabatic beds with interstage cooling. The effectiveness factor is calculated by the empirical relationship of Dyson and Simon (1968) as well as through the solution of the two-point boundary value differential equation of diffusion and reaction in the porous catalyst pellets. The two-point boundary value differential equation is solved by using the orthogonal collocation method. The results of the heterogeneous model are compared with the results of an industrial reactor and calculated values presented by Singh and Saraf (1979). They used the shooting technique to solve the two-point boundary value differential equation to calculate the effectiveness factor. This shooting technique took a large number of iterations and did not satisfy the boundary conditions at the center due to the divergence of the iterative path (Singh, 1985) while the collocation method used in this investigation did not given any of these difficulties. The heterogeneous model developed is used to predict the optimal temperature and ammonia concentration profiles. Optimization investigations are carried out based on the intrinsic rate of reaction as well as the actual rate of reaction, i.e., the product of the rate of reaction and the effectiveness factor. The heterogeneous model is also used to obtain the corresponding ammonia concentration profile for the optimal temperature profile obtained from the pseudohomogeneous model. Constrained optimization is also investigated to obtain the suboptimal policy and performance that does not damage the catalyst.
Rate Expression The intrinsic modified form of the Temkin rate expression (Dyson and Simon, 1968) is used in the developed model, as follows: (1) RNHs = kz[K2fN,(fH,3/fNH,2)a- (fNH,2/fH,3)1-"l where RNHsis the reaction rate in kmol of NH3/ (h.m3 of catalyst bed), k2 is the velocity constant for the reverse reaction in kmol/(h-m3),and K, is the equilibrium constant of the reaction
* Author
(2) f/,N, + Y,H2 + NH, where f N 2 , f ~and ~ fNH3 , are the fugacities of nitrogen, hydrogen, and ammonia, respectively and a is a constant. The velocity constant, kz,is estimated by the Arrhenius relation of the form (Singh and Saraf, 1979) 122 = kzo exp(-Ez/RJ') (3)
The respective values for the Montecatini Edison catalyst are a = 0.55, Ez = 1.635 X lo8 J/kmol, and log k,, = 14.7102, and RG is the universal gas constant and T is absolute temperature in K. The equilibrium constant, K,, is calculated from (Gaines, 1979) log K , = -2.691122 log T - 5.519265 X 10-5T + 1.848863 X 10-'T2 + 2001.6/T + 2.6899 (4) The fugacity of component i is given by definition as fi = r$iXiP (5) where r$i is the fugacity coefficient of component i, X iis the mole fraction of component i, and P is the total pressure. The fugacity coefficients are calculated by means of equations given by Mahfouz (1985). A drawback of eq 1 is that it is obviously not valid for very low ammonia concentrations since the first term diverges. Model Development A reactor having three adiabatic beds with Montecatini Edison catalyst and interstage cooling is used to produce ammonia according to the exothermic reaction given by eq 2. A schematic diagram of the reactor is given in Figure 1. 1. T h e Overall Reactor Model. A heterogeneous one-dimensional model is developed for the reactor, and the following assumptions are made: 1. The reactor is operating at steady state. 2. There is negligible heat-transfer resistance between the pellets and the gas (Vancini, 1971). 3. Axial dispersion is negligible due to the high gas velocities (Reddy and Husain, 1982). 4. The thermal and concentration gradients in the radial direction are negligible (Shah, 1967). 5. The pressure drop along the reactor is negligible (Baddour et al., 1965). 1.1. Material Balance on t h e Bulk Gas. A molar differential balance for nitrogen in the catalyst bed gives
to whom all correspondence should be addressed.
0888-5885f 88f 2627-2015$01.50f 0 0 1988 American Chemical Society
2016 Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988 7 = bo
+ blT + b2Z + b3T2+ b4Z2 + b5T3 + b6Z3
(9)
A more rigorous approach for the treatment of the effectiveness factor problem is through formulating the diffusion-reaction equations, solving them, and calculating the effectiveness factor. The catalyst pellet equations are developed using the well-established assumptions presented by previous investigators (Mahfouz, 1985; Dyson and Simon, 1968; Singh and Saraf, 1979). A molar differential balance for component i inside the catalyst particle gives
where Ni is the molar flux of the ith component in the r direction, E is the void fraction of the packed bed and is equal to 0.46 (Dyson and Simon, 1968), and y i is the stoichiometric coefficient of the ith component in the reaction scheme given in eq 2. Following the ordinary convention, yi is positive for products, negative for reactants, and zero for inerts. Nitrogen, hydrogen, ammonia, methane, and argon have been designated as Components 1,2, 3,4, and 5, respectively. The boundary conditions for eq 10 are at r = 0.0 Ni = 0.0 dXi/dr = 0.0
A
at r = RP xi OUTLET d.PRODUCT
'4'
lNLEl
FEED
Figure 1. Schematic diagram of a three-bed reactor.
where Z is the fractional conversion of nitrogen, V is the volume of catalyst bed in m3, 7 is the effectiveness factor, and FN: is the initial molar flow rate of nitrogen in kmol/ h. The fractional conversion of nitrogen at any cross section of the bed can be written as Z = (molar flow of Nz at inlet - molar flow of N2 at cross section)/molar flow of N2 at inlet (7)
=
xi,
(11)
where RP is the radius of the sphere equivalent to an inm for dustrial sized particle and is equal to 2.85 X industrial particles of 6-10 mm (Dyson and Simon, 1968) and Xi, is the mole fraction at the surface of the catalyst particle and this will be assumed to be the same as in the gas phase. Following some lengthy but straightforward manipulation of the equations (Mahfouz, 1985),we end up with the following dimensionless equation for the catalyst pellets:
All mole fractions can be expressed in terms of the feed mole fractions and the fractional conversion of nitrogen. 1.2. Energy Balance on the Bulk Gas. The energy balance for a differential element of the catalyst bed gives dT dV -
( - ~ R ) ' I R N H , ( Z~TP)
(8)
&CPmix
where AHR is the heat of reaction in J/kmol of NH3, m is the total mass flow rate in kg/h, and C p is the specific heat of reacting gas mixture in J/(kgKT The heat of reaction has been calculated by using the empirical correlation presented by Strelzoff (1981). The molar specific heat of different components is calculated from the equations given by Shah (1967), and the average molar specific heat of the mixture varies with the degree of conversion. 2. Effectiveness Factor and the Catalyst Particle Equations. The effect of diffusional resistance inside the catalyst which is important for the industrial catalyst particles (6-12 mm) is expressed in terms of the effectiveness factor, 7. Dyson and Simon (1968) have derived an empirical correlation to obtain the effectiveness factor for ammonia synthesis reaction as a function of temperature and conversion: ,
and the boundary conditions in terms of w are at w = 0.0 dXi/dw = 0.0 at
w
= 1.0
xi = xi,
(13)
where w = r/Rp. The total concentration, C, is found from
5fi e = -i-= l
RT where fi is the fugacity of component i. The effective diffusion coefficient is calculated from (Wheeler, 1955) 1 D , =-OD, (15) le 2 where 6 is the intraparticle porosity and Di is the bulk diffusion coefficient of component i. The intraparticle porosity is approximately 0.52 according to Mahfouz
Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988 2017 (1985). The diffusion coefficients at 0 "C and 1 atm are given by 1 - xig D: = n i = 1-3 (16)
The diffusion coefficients calculated from eq 16 are then corrected to the temperature and pressure at the surface of the catalyst pellet by
where P is in atmospheres. Equation 12 is a two-point boundary value differential equation and is solved in this work by the orthogonal collocation method (Villadsen and Stewart, 1967; Finlayson, 1972), and the effectiveness factor is evaluated with high accuracy by using the summation formula of the orthogonal collocation method.
Optimization Several investigators carried out optimization studies on the ammonia reactors by using the simplified pseudohomogeneous models (Gmelin Handbook of Inorganic Chemistry, 1936; Mansson and Andersen, 1986);of course the use of such models which do not take into account all physical and chemical processes occurring inside the catalyst particles does not describe the system accurately enough. In this study, the developed heterogeneous model which is checked against industrial results is used to predict the optimal temperature and concentration profiles for adiabatic ammonia reactor. The new optimization using the heterogeneous model is based on two different bases, the intrinsic rate expression as well as the actual rate of reaction. This is the first realistic optimization study using a heterogeneous model to be reported in the literature for the ammonia reactor. The optimal temperature and concentration profiles using the intrinsic rate expression can be found by numerically solving the optimality condition aRNHs(X,T,P)/aT = 0.0
(18)
at each point along the reactor simultaneously with the mass balance equation of the bulk fluid (eq 6) to determine the optimal temperature and concentration profiles along the length of the reactor. The optimization based on the maximization of the intrinsic rate of reaction does not give the absolute maximum. The absolute maximum for the heterogeneous model can be obtained by considering the optimization problem based on the actual rate of reaction. The optimality condition a(TRN,,(x,T,p))/aT = 0.0
(19)
combined with the material balance equation (eq 6) is used to predict the optimal temperature and concentration profiles along the reactor. At the entrance of the reactor, the optimal temperature is usually very high, causing a very high rate of reaction; accordingly the effectiveness factor is very small, and this is the pore diffusion control regime. In this region the collocation method gives accurate results only with a large number of collocation points (7-10 points), while the shooting method drastically fails. However a good approximate value of the effectiveness factor can be obtained for this region by using the effective reaction zone method (Gavalas, 1968; Petersen, 1965; Cresswell, 1970). The ap-
plication of the method to our case proved to give excellent evaluation of the values of 9 in this region.
Computational Algorithms Nitrogen is considered as the base component for calculations. The initial feed composition, temperature, flow rate, pressure inside the reactor, and the volume of the catalyst for each bed are taken from the work of Singh and Saraf (1979). For the simulation studies, the pseudohomogeneous model can be solved by direct integration of the mass and heat balance initial value differential equations. However, for the heterogeneous model, the effectiveness factor is computed by using the orthogonal collocation method at each point during the integration of the bulk phase mass and heat balance initial value differential equations. Automatic step size in the integration of the initial value differential equations was used to ensure accuracy. The following algorithms are used for solving the optimization problems: A. Optimization Based on the Intrinsic Rate of Reaction for the Pseudohomogeneous Model (Opt. Hom. Modl.). The optimal temperature profile for the homogeneous model is obtained in the same way as for the heterogeneous model, maximizing the intrinsic rate expression described in section B except now the effectiveness factor is equal to unity in the material balance equation all along the length of the reactor. B. Optimization of the Heterogeneous Model Based on the Intrinsic Rate of Reaction (Opt. Het. Modl. Intr.). 1. The optimality condition (eq 18) is solved to find the optimal temperature at known conversions at each point along the length of the reactor. Since the intrinsic rate expression is quite complex because of the fugacity coefficients involved, numerical methods for finding the optimum temperature must be used. Fibonacci search is used to determine the temperature at which the rate of reaction (eq 1) is maximum at known conversion at a certain depth within the catalyst bed. In the case under investigation, the conversion is considered zero at the inlet of the first bed. 2. After the optimal temperature at the entrance of the reactor is computed, the effectiveness factor can be calculated by using either the subroutine that calculates the effectivenessfactor by the orthogonal collocation technique or the subroutine that calculates the effectiveness factor for the transport-limited regime depending on the value of the optimal temperature. 3. From the computed temperature, effectiveness factor, and knowledge of the concentration at the entrance point (feed condition), the material balance equation (eq 6) is integrated by the Runge-Kutta subroutine with automatic step size to determine the conversion for the next step. 4. Steps 1, 2, and 3 are repeated for the whole length of the reactor to give the optimal temperature and concentration profiles along the length of the reactor. C. Simulation Using the Optimal Temperature Profile of the Pseudohomogeneous Model in the Heterogeneous Model (Hom. Het. Modl.). In some cases, for the sake of comparison, the optimal temperature profile obtained from the pseudohomogeneous model is used to compute the change of concentration along the length of reactor using heterogeneous model, in order to check the implications of using the optimal control policy obtained from the pseudohomogeneous model on the actual reactor. D. Optimization of the Heterogeneous Model Based on the Actual Rate of Reaction (Opt. Het. Modl. Act.). 1. The optimality condition (eq 19) is solved to find the
2018 Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988
optimal temperature at known conversion at each point along the length of the reactor. Fibonacci search is used to determine the temperature at which the actual rate of is maximum at a known converion reaction (qRRNHI(X,T,P)) at each point along the length of the reactor. In this case, 7 is evaluated a very large number of times to find the maximum of the actual rate of reaction at each point along the length of the reactor, and the orthogonal collocation method makes this possible because of its high efficiency. The computation is initiated at zero conversion at the inlet of the reactor, and the optimal temperature is found at that conversion. 2. From the known value of (qRNHs(X,T,P))at the optimal temperature at the entrance of the reactor, the material balance equation (eq 6) is integrated by the Runge-Kutta subroutine with automatic step size to determine the conversion for the next step size. 3. Steps 1 and 2 are repeated for each step size to predict the temperature and conversion profiles along the length of the reactor. E. Equilibrium Compositions and Temperatures. If the rate of reaction is equated to zero at all temperatures along the length of the reactor, the equilibrium concentration profile along the length of the reactor can be obtained. The equilibrium temperature profile can be obtained by equating the rate of reaction to zero at all conversions along the length of the reactor. In order to solve the equilibrium equations, the Secant method is used (Gerald, 1970).
Results and Discussion The results of reactor modeling and simulation using the heterogeneous model and the orthogonal collocation method for the solution of the pellet equations are presented here with the pellet effectiveness factor profiles and discussion of the reasons for its nonmonotonic behavior. This is the first time that these profiles have been published. Modeling and Simulation. Table I gives the results of different models tested. Good agreement was obtained between the computed results using the heterogeneous models and the reactor performance; this is in contrast to the pseudohomogeneous model which deviates largely from the reactor performance, as shown in Table I. The three-collocation method gives the best results for ammonia mole percent and temperature prediction, with much less computational effort compared with the shooting method. The accuracy of the collocation method can be improved numerically as much as it is desired by increasing the number of collocation points without increasing the computational time considerably as shown in Table I for one, two, and three collocation points. The collocation method is a very effective method, and it uses much less computation time, which make it suitable for the optimization study that follows. Figure 2 shows the effectiveness factor profile along the reactor. The effectiveness factor profile of the first bed shows a nonmonotonic behavior along the length of the bed. The effectiveness factor is a measure of the effect of the diffusional limitations on the rate of reaction. For negligible diffusional resistances, the effectiveness factor should be equal to 1.0. For constant rate of reaction, increasing the diffusional resistances causes the effectiveness factor to decrease. For constant diffusional resistance, the effectiveness factor increases as the rate of reaction decreases. In the present case the intrinsic rate of reaction, diffusional resistances, and the equilibrium constant are changing along the length of the reactor. A similar effectiveness factor profile with a maximum
0
!
0
'
,
'
I
3 95
J
'
~
'
"
190 VOL
'
'
"
"
"
"
11 85
~ ~ 15 80
' 19 75
OF C A T B E D ( m 3 )
Figure 2. Effectiveness factor profile (using three collocation points). 9 00
300 0
-1
I
3 95
7 90
11 8 5
V O L . OF. C A T . B E D
15 8 0
19 7 5
(173')
Figure 3. Optimal temperature profile (Opt. Hom. Modl).
(nonmonotonic) for steam reforming of natural gas was shown by DeDeken et al. (1982). This nonmonotonic behavior of 7 with reactor length is due to complex interaction between rate of reaction, diffusional resistances, and equilibrium constants of reversible reactions. This phenomenon is discussed elsewhere (Elnashaie et al., 1988).
Optimization of the Reactor Performance The optimal temperature profile that gives the highest conversion to ammonia is obtained for different models by using algorithms described in sections A, B and D. Such an optimization investigation using the heterogeneous model and using two different bases for optimization, Le., maximizing the intrinsic rate and maximizing the actual rat of reaction, is made possible because of the high efficiency of the orthogonal collocation method. The optimal results are compared with the heterogeneous model results rather than the experimental results for the sake of consistency of the comparison. The optimization investigation presented in this paper does not deal with the very important practical problem of the implementation of the optimal policy. This problem is at present being investigated in this department both theoretically and experimentally. A heterogeneous model using three internal collocation points is used throughout this optimization study. Optimization Based on the Intrinsic Rate of Reaction for Pseudohomogeneous Model (Opt. Hom. Modl.). The optimal temperature and ammonia concentration profiles obtained by using the pseudohomogeneous model based on maximizing the intrinsic rate expression are shown in Figures 3 and 4, respectively. The equilib-
Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988 2019
I
I
-
"
0 3.95
0
7.90
11.85
15.80
19.75
'
,
0
'
" 3.95
VOL. OF. CAT. BED ( m3)
Bed-I
~
Bid
'
~
'
"
'
'
7.90
"'"""4 11.85
15.80
19.75
V O L . OF. C A T . BED ("1
Figure 4. Optimal concentration profiles (Opt. Hom. Modl. and Hom. Het. Modl).
-
Bed -111
B i d I1
Bed-I
- I1
8rd
Figure 6. Optimal concentration profiles (Opt. Het. Modl. Intr. and Opt. Het. Modl. Act.).
""I
- 111
o c I u ( ~ Irate 01 reacllon
inlrinric rate el reaction U
3
nom. hct.
0.60
In
v)
W W z
?
I-
W
420
-
300
LL W
"
0
~
~
' . 3.95
"
'
~ ' ' 7.90
'
'
L ' 11.85
l
l
i l l 15.80
l
l
' 19.75
VOL. OF. CAT. BED ( m 3 i
Figure 5. Optimal temperature profiles (Opt. Het. Modl. Intr. and Opt. Het. Modl. Act.).
rium curve shown in Figure 3 is the temperature at which the optimal reaction mixture would be in equilibrium, whereas the equilibrium curve in Figure 4 is an equilibrium ammonia concentration at the optimal temperatures. The optimal temperature profile is very high at the inlet of the reactor (Figure 3) where the ammonia mole percent is very low (Figure 4), resulting in the highest rate of reaction. The first bed plays an important role due to the high temperature zone that prevails in it, especially at the entrance of the reactor, resulting in high ammonia mole percent as shown in Figure 4. This bed gives 15.75 mol % ammonia, which is equivalent to 24.51 % improvement with respect to the pseudohomogeneousmodel simulation results. This means that the first bed is the backbone for the improvement obtained by the optimization technique. Optimization of the Heterogeneous Model Based on the Intrinsic Rate of Reaction (Opt. Het. Modl. Intr.). Optimization in this part is performed on the heterogeneous model based on the maximization of the intrinsic rate of reaction (eq 18). The optimal temperature and ammonia concentration profiles obtained from such policy are presented in Figures 5 and 6, respectively. The optimal temperature is quite high (Figure 5) at the entrance of the reactor, giving a high rate of intrinsic reaction but low values of the effectiveness factor, as shown in Figure 7. This means that there is a large resistance for diffusion, causing a very steep concentration profile inside the catalyst pellet. The solution of the pellet equations in this region has been discussed earlier. The effectiveness factor for the optimal temperature profile in this case ranges from 0.058 at the inlet of the first
0
3.95
7.90
11.85
15.80
19.75
V O L . O F . C A T . BED (m3)
Figure 7. Optimal effectiveness factor profiles for all optimization cases.
bed to 0.86 at the exit from the third bed (Figure 7). This large variation of the effectiveness factor demonstrates clearly the importance of the accurate determination of its values. It also explains the unrealistic prediction that can be obtained from the pseudohomogeneous model for reactor performance where diffusional resistances are not negligible. Table I1 gives the results for the heterogeneous model with the optimization based on the intrinsic rate of reaction for the case of unconstrained optimization as well as the case where the optimal temperature profile is constrained with the maximum temperature that the iron catalyst can withstand, which is about 800 K accoridng to Mansson and Andresen (1986). For unconstrained optimization, the first bed gives 11.83 mol % ammonia, which is equivalent to 11.18% improvement in ammonia concentration. The second and third beds give 17.23 and 20.48 mol % of ammonia, respectively, which are equivalent to improvements of 9.89% and 6.06% over the adiabatic reactor with interstage cooling. For the constrained optimization, the improvement over the simulation of the practical system is higher than the case of optimization of the heterogeneous model based on the intrinsic rate of reaction (see Table 11). This is due to the fact that in this case the unconstrained optimization does not take into consideration the harmful effect of the excessively high temperature on the effectiveness factor at the entrance of the reactor, while the constrained case limits this harmful effect to some extent by limiting the temperature at the entrance of the reactor.
2020 Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988 Table I. Comparison between Experimental Data and Calculated Results (Total Feed Flow, 242 160 (N m8)/h; Pressure, 226 atm) composition, mol 70 N2 H2 NH3 CH, Ar temp, "C conversion Bed I: Volume of the Catalyst, 4.75 m3 Inlet 22.19 67.03 2.76 5.46 2.56 385.000 0.0000
Outlet exptl calcd (shooting) Hom. Modl. Het. Modl. (Pol. Dyson) Het. Modl.: 1 col. pt 2 col. pts 3 col. pts
20.10 20.03 19.51 20.32
61.00 60.58 59.05 61.44
10.50 10.74 12.65 9.68
5.70 5.89 5.99 5.83
2.70 2.76 2.81 2.73
507.000 504.000 535.268 493.311
0.1578 0.1624 0.1978 0.1422
19.92 20.04 20.06
60.26 60.61 60.67
11.15 10.71 10.64
5.91 5.88 5.88
2.77 2.76 2.76
513.892 507.530 506.470
0.1701 0.1618 0.1604
Bed 11: Volume of the Catalyst, 7.2 m3 Inlet exptl calcd (shooting) Hom. Modl. Het. Modl. (Pol. Dyson) Het. Modl.: 1 col. pt 2 col. pts 3 col. pts
20.10 20.03 19.51 20.32
61.00 60.58 59.05 61.44
10.50 10.74 12.65 9.68
5.70 5.89 5.99 5.83
2.70 2.76 2.81 2.73
433.000 433.000 433.000 433.000
0.1578 0.1624 0.1978 0.1422
19.92 20.04 20.06
60.26 60.61 60.67
11.15 10.71 10.64
5.91 5.88 5.88
2.77 2.76 2.76
433.000 433.000 433.000
0.1701 0.1618 0.1604
18.20 18.58 18.28 18.83
57.10 56.25 55.35 56.33
15.90 16.11 17.22 15.19
6.10 6.17 6.23 6.12
2.70 2.89 2.92 2.87
502.000 505.000 495.777 510.884
0.2555 0.2591 0.2780 0.2431
18.60 18.68 18.69
56.33 56.56 56.60
16.01 15.72 15.68
6.16 6.15 6.15
2.89 2.88 2.88
501.079 503.670 504.148
0.2574 0.2524 0.2517
Outlet exptl calcd (shooting) Hom. Modl. Het. Modl. (Pol. Dyson) Het. Modl.: 1 col. pt 2 col. pts 3 col. pts
Bed 111: Volume of the Catalyst, 7.8 m3 Inlet exptl calcd (shooting) Hom. Modl. Het. Modl. (Pol. Dyson) Het. Modl.: 1 col. pt 2 col. pts 3 col. pts
18.20 18.58 18.28 18.83
57.10 56.25 55.35 56.99
15.90 16.11 17.22 15.19
6.10 6.17 6.23 6.12
2.70 2.89 2.92 2.87
415.000 415.000 415.000 415.000
0.2555 0.2591 0.2780 0.2431
18.60 18.68 18.69
56.33 56.56 56.60
16.01 15.72 15.68
6.16 6.15 6.15
2.89 2.88 2.88
415.000 415.000 415.000
0.2574 0.2524 0.2517
exptl. calcd (shooting) Hom. Modl. Het. Modl. (Pol. Dyson) Het. Modl.: 1 col. pt 2 col. pts 3 col. pts
17.80 17.55 17.35 17.77
53.90 53.16 52.58 53.85
19.10 19.93 20.66 19.08
6.30 6.37 6.41 6.33
2.90 2.99 3.01 2.97
455.000 463.000 458.559 465.859
0.3091 0.3226 0.3343 0.3088
17.65 17.70 17.71
53.50 53.65 53.67
19.52 19.33 19.31
6.35 6.34 6.34
2.98 2.97 2.97
460.397 461.868 462.105
0.3160 0.3129 0.3126
Outlet
Later in this paper, this harmful effect at the entrance of the reactor can be rigorously minimized by optimization of the heterogeneous model based on the actual rate of reaction. Simulation Using the Optimal Temperature Profile of the Pseudohomogeneous Model in the Heterogeneous Model (Hom. Het. Modl.). The optimization results obtained from the pseudohomogeneous model do not represent the actual situation because it neglects masstransfer resistances in the catalyst pellets which do exist in the reactor. However, the results of the pseudohomogeneous model are important because they give higher limit for production when the diffusion resistances are decreased by using smaller catalyst particles (which causes high pressure drop). In this section, the results of reactor simulation using the heterogeneous model and the optimal temperature profile obtained from the pseudohomogeneous
model are presented. This simulation is important because it gives the behavior of the actual reactor if the optimal temperature profile applied to it is the one obtained from the pseudohomogeneous model. Using the optimal temperature profile of the pseudohomogeneous model (Figure 3) in the heterogeneous model gives the corresponding optimal concentration profile shown in Figure 4. From this figure, the optimal concentration profile shows similar behavior to the pseudohomogeneous, one but it obviously shifted to lower values. The results of this approach and constrained (limited) temperature policy are shown in Table 111. The constrained temperature policy shows almost similar results with a slight difference in the first and the third beds (Table 111) in comparison with the unconstrained case. It is interesting to notice that application of the optimal temperature profile obtained from the pseudohomogeneous model on the heterogeneous
Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988 2021 Table 11. Comparison between the Heterogeneous Model and Optimization of the Heterogeneous Model Based on the Intrinsic Rate of Reaction (Opt. Het. Modl. Intr.) Opt. Het. Modl. Intr. Het. Modl. limited temo Bed I NH3, mol % 10.64 11.83 12.15 0.1604 0.1828 conversn 0.1887 temp, O C 506.740 492.906 489.865
NH,, mol % conversn temp, O C
Bed I1 15.68 17.23 0.2517 0.2781 504.148 451.395
NH,, mol % conversn temo, OC
Bed I11 19.31 20.48 0.3126 0.3314 462.105 430.166
17.38 0.2807 450.374
20.57 0.3328 429.532 % improvement of NH3 concn with respect to the heterogeneous model bed I bed I1 bed I11 Opt. Het. Modl. Intr. 11.18 9.89 6.06 Opt. Het. Modl. Intr. (limited temp) 14.19 10.84 6.53 Table 111. Comparison between the Heterogeneous Model and Simulation of the Heterogeneous Model Using the Optimal Temperature Profile Obtained from the Homogeneous Model (Hom. Het. Modl.) Hom. Het. Modl. Het. Modl. limited temp Bed I NH3, mol % 10.64 12.26 12.25 conversn 0.1604 0.1907 0.1905
NH3, mol % conversn
Bed I1 15.68 17.44 0.2517 0.2817
17.44 0.2817
Bed I11 19.31 20.58 20.59 0.3126 0.3330 0.3332 % improvement of NH3 concn with respect to the heterogeneous model bed I bed I1 bed I11 Hom. Het. Modl. 15.23 11.22 6.58 Hom. Het. Modl. (limited temp) 15.13 11.22 6.63 NH3, mol % conversn
model gives results which are better than that obtained from the optimal temperature policy obtained from the heterogeneous model based on the maximization of the intrinsic rate of reaction, showing clearly that optimization for the heterogeneous model based on the intrinsic rate of reaction is not really optimal. Optimization of the Heterogeneous Model Based on Actual Rate of Reaction (Opt. Het. Modl. Act.). The optimization of the heterogeneous model based on the maximization of the intrinsic rate expression does not give an actual maximum because the optimal temperature is not predicted from the actual rate of reaction, especially at very low values of the effectiveness factor at which diffusional resistance has a great effect on the value of the actual rate of reaction. Optimization based on the actual takes more computational rate of reaction (T&'~~(X,T,P)) effort because it requires the evolution of the 9 many times at each point of the reactor to give at each point the temperature that maximizes qRNH>(X,T,P) and not RNH3(X,T,P) only. Such an optimization computation is carried out to obtain the actual optimal temperature and con-
Table IV. Comparison between the Heterogeneous Model and Optimization of the Heterogeneous Model Based on the Actual Rate of Reaction (Opt. Het. Modl. Act.) Opt. Het. Modl. Act. Het. Modl. limited temD Bed I NH3, mol % 10.64 12.38 12.34 0.1604 0.1929 0.1922 conversn temp. " C 506.470 472.082 472.25
NH,, mol % conversn temp, O C
Bed I1 15.68 17.57 0.2517 0.2838 504.148 442.192
17.55 0.2835 442.240
Bed I11 19.31 20.71 20.70 0.3126 0.3351 0.3349 462.105 424.551 424.74 % improvement of NH3 concn with respect to the heterogeneous model bed I bed I1 bed I11 Opt. Het. Modl. Act. 16.35 12.05 7.25 Opt. Het. Modl. Act. (limited temp) 15.98 11.93 7.20 NH,, mol % conversn temo, " C
~~
centration profiles for the ammonia reactor. The optimal temperature and ammonia concentration obtained by this method are shown in Figures 5 and 6, respectively. The actual optimal temperature profile (Figure 5) is lower than that obtained by optimization based on the maximization of the intrinsic rate expression especially at the first bed, resulting in higher optimal concentration profile as shown in Figure 6. The optimal effectiveness factor is also shown in Figure 7. It lies between the optimal effectiveness factor profiles of the optimization of the heterogeneous model based on the intrinsic rate of reaction and the one based on the insertion of the optimal temperature profile of the pseudohomogeneous model into the heterogeneous model. The optimal effectiveness factor from optimization based on actual rate of reaction is higher than that obtained from optimization based on the intrinsic rate expression. From Table IV, the optimization based on actual rate of reaction gives 12.38, 17.57, and 20.71 mol % ammonia from the first, second, and the third beds, respectively and these are equivalent to 16.35%, 12.05%, and 7.25% improvement of ammonia concentration over the adiabatic operation with interstage cooling. The limited temperature policy is also implemented here, giving slightly lower performance (Table IV) than that of unconstrained optimization based on the actual rate of reaction. This limited temperature gives 15.98%, 11.93%, and 7.2% improvement of ammonia concentration from the first, second, and third beds, respectively, over that obtained from the adiabatic beds with interstage cooling. All optimal policies show that the overall increase of ammonia concentration due to these policies ranges from 6.06% to 7.25% improvement in comparison with the performance of adiabatic reactor with interstage cooling predicted from the accurate heterogeneous model. These improvements are very appreciable. To get an idea about the dollar value of such improvements (6-7% for the three beds, 9-12% for two beds, and 11-16% for one bed), a simple calculation is carried out for the dollar value of only 1% improvement in the % NH3 in the outlet of the reador for the Saudi ammonia production which is about 3600 MTPD. On the basis of the assumption of 330 working days a year and $100.0 for a metric ton of ammonia, this 1% corresponds to a dollar value of about $1.2 X 106/year; for Canada (production rate of about 11000 MTPD) the
2022 Ind. Eng. Chem. Res., Vol. 27, No. 11, 1988
dollar value of 1% improvement i n 70 NH:, is equivalent to $4.0 X 106/year. Therefore, for a 6-7 % improvement, these figures will become $7.9-9.2 X 106/year for Saudi Arabia and $24-28 X 106/year for Canada. Another important point to notice is that the percent improvement from the optimization results is much higher for small beds than for larger beds. Therefore, when the economics of the process dictate the use of smaller beds, the optimal performance will give improvements which are m u c h higher than the figures given.
Conclusions The collocation method is found to be suitable and efficient for the solution of the catalyst pellets two-points boundary value differential equation for the ammonia synthesis reaction. The heterogeneous model is used to find the o p t i m a l
temperature profile that maximizes the production of ammonia, and it is shown that the optimization based on the actual rate of reaction gives the true optimal policy compared with the policy obtained from optimization of the heterogeneous model based on the intrinsic rate of reaction. The improvement of the reactor performance using the t r u e optimal policy is found to be quite appreciable, amounting to 6-7% increase in percentage ammonia at the exit of the three beds. The constrained optimization using a maximum temperature limit on the reactor temperature gives % NH, values which are very close to those obtained from the unconstrained optimization, and the improvement in performance discussed earlier is maintained. The heterogeneous model used in this study is at the present time being perfected to match exactly the practical results for a large n u m b e r of reactor configurations, adiabatic and nonadiabatic, and the practical implementation of the optimal control policy is being investigated for different reactor configurations.
Nomenclature C = total concentration Cpw = specific heat of reacting gas mixture, J / ( k g K ) D i= diffusion coefficient of component i D: = diffusion coefficient of component i at 0 “C and 1a t m Dj? = diffusion coefficient of component j in component i E, = activation energy for ammonia decomposition, J/kmol Die = effective diffusion coefficient of component i f i = fugacity of component i f N 2 , fH2, fNHs = fugacity of nitrogen, hydrogen, and ammonia, respectively F N = molar flow rate of nitrogen a t reactor inlet, kmol/h A d R = heat of reaction, J / k m o l of NH, k 2 = velocity constant of reverse reaction k,, = frequency factor in Arrhenius equation for k 2 K , = equilibrium constant of reaction 2 m = total mass flow rate, kg/h N , = molar flux of component i in r direction P = pressure, a t m r = radial coordinate of spherical catalyst particle R N H B = rate of ammonia formation, kmol of NH3/(h.m3 of catalyst bed) RG = universal gas constant RP = radius of spherical particle T = temperature, K V = volume of catalyst bed, m3 X I = mole fraction of component i X , , = mole fraction of component i in gas (bulk phase) 2 = conversion based on nitrogen
”
Greek Symbols cy = kinetic parameter y 1 = stoichiometric coefficient of component i $L= fugacity coefficient of component i 9 = effectiveness factor 6 = void fraction of packed bed w = dimensionless radial coordinate of a spherical particle % = intraparticle porosity Registry No. H2, 1333-74-0;N2,7727-37-9; NH3, 7664-41-7.
Literature Cited Annable, D. “Application of the Temkin Kinetic Equation to Ammonia Synthesis in Large Scale Reactors”. Chem. Eng. Sci. 1952, I, 145-154. Baddour, R. F.; Brain, P. L. T.; Logeais, B. A.; Eymery, J. P. “Steady-State Simulation of an Ammonia Synthesis Converter”. Chem. Eng. Sci. 1965,20, 281-292. Cresswell, D. L. “On the Uniqueness of the Steady State of a Catalyst Pellet Involving Both Intraphase and Interphase Transport”. Chem. Eng. Sci. 1970,25, 267-275. Dyson, D. C.; Simon, J. M. “A Kinetic Expression with Diffusion Correction for Ammonia Synthesis on Industrial Catalyst”. Ind. Eng. Chem. Fundam. 1968, 7,605-610. DeDeken, J. C.; Devos, E. F.; Froment, G. F. “Steam Reforming of Natural Gas: Intrinsic Kinetics, Diffusional Influences, and Reactor Design”. React. Eng. ACS Symp. Ser. 1982,196, 181-197. Elnashaie, S. S.; Mahfouz, A. T.; Elshishini, S. S. “Digital Simulation of an Industrial Ammonia Reactor”. Chem. Eng. Process. 1988, 23(3), 165-177. Finlayson, A. ”The Method of Weighted Residuals and Variational Principles”. In Mathematics in Science and Engineering; Bellman, R., Ed.; Academic: New York, 1972; Vol. 87. Gaines, L. D. “Optimal Temperature for Ammonia Synthesis Converters”. Ind. Eng. Chem. Process Des. Dev. 1977, 16, 381-388. Gaines, L. D. “Ammonia Synthesis Loop Variables Investigated by Steady-State Simulation”. Chem. Eng. Sci. 1979, 34, 37-50. Gerald, C. F. Applied Numerical Analysis; Addison-Wesley Series in Mathematics; Varga, R. S.,Ed.; Addison-Wesley: Reading, MA, 1970. Gavalas, G. R. Nonlinear Differential Equations of Chemically Reacting System; Springer-Verlag: New York, 1968. Gmelin Handboek of Inorganic Chemistry Verlag Chemie: Berlin, 1936. Kjaer, J. Calculation of Ammonia Converters on an Electronic Digital Computer; Akademisk Forlag: Copenhagen, 1963. Mahfouz, A. T. “Steady-State Modeling and Simulation of Industrial Ammonia Synthesis Reactor”. PbD. Dissertation, The University of Cairo at Cairo, Egypt, 1985. Mansson, B., Andresen, B. “Optimal Temperature Profile for an Ammonia Reactor”. Ind. Eng. Chem. Process Des. Dev. 1986,25, 59-65. Petersen, E. E. Chemical Reaction Analysis; Prentice Hall: Englewood Cliffs, NJ, 1965. Reddy, K. V.; Husain, A. “Modelingand Simulation of an Ammonia Synthesis Loop”. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 359-366. Shah, M. J. “Control Simulation of Ammonia Production”. Znd. Eng. Chem. 1967,59, 72-83. Singh, C. P. P. “A Generalized Method for Calculating Effectiveness Factor”. Chem. Eng. Sci. 1985, 40, 2375-2378. Singh, C. P. P.; Saraf, D. N. “Simulation of Ammonia Synthesis Reactors”. Ind. Eng. Chem. Process Des. Dev. 1979,18,364-370. Strelzoff, S. Technology and Manufacture of Ammonia; Wiley: New York, 1981. Villadsen, J. V.; Stewart, W. E. ”Solution of Boundary-Value Problems by Orthogonal Collocation”. Chem. Eng. Sci. 1967, 22, 1483-1501. Vancini, C. A. Synthesis of Ammonia; The MacMillan Press: London, 1971. Van Heerden, C. “Autothermic Process Properties and Reactors Design”. Znd. Eng. Chem. 1953, 45, 1242-1247. Wheeler, A. Catalysis; Emmett, P. H., Ed.; Reinhold: New York, 1955. Received for review November 25, 1987 Revised manuscript received May 17, 1988 Accepted June 12, 1988
