en a computer is used for control of a chemical w r o c e s s , a mathematical description of the process units is necessary for the computer to perform such tasks as feedforward and feedback control, as well as off-lme and on-line optimization. Often, historical plant data are used to obtain these relationships by means of regression analysis. Unfortunately, the reliability of regression models is limited by the system noise. Also, these models are valid only for the range of the values of variables in which the data was collected, and extrapolation of the relationships is risky. If a priori knowledge of the behavior of one or more critical units in the plant is available, the design of computer control strategy is greatly facilitated. Off-line simulation of the actual plant units by means of analog or digital computers will provide such a knowledge, provided that the mathematical model chosen for simulation describes the plant with reasonable accuracy. Once the accuracy of the simulation model in predicting the plant behavior is established, the same model can be reduced by various simplification techniques (8,73) to obtain mathematical relationships of the plant unit for use in on-line control and optimization. I n this paper, we will describe a simulation study of an ammonia synthesis reactor-one of the critical units in the ammonia plant. In Figure 1, we show a flow diagram of an ammonia process. One can see that aside from the synthesis reactor, the primary and secondary reformers, as well as the shift converter, are also of importance in the overall control strategy for the ammonia plant. However, a study of the reformers and CO converter units, which are involved in the production of synthesis gas, will logically be the subject of later investigation. It is recognized that a given ammonia plant may have operating constraints which must be considered in the ovepall optimization; for example, the compressor capacity may be limiting in some plants,
THIS CASE S T U D Y O F THE MATHEMATICAL SIMULATION OF STEADY-STATE BEHAVIOR IN A MULTIPLE BED, QUENCH-TYPE REACTOR VIVIDLY DEMONSTRATES BENEFITS OF MODERN SIMULATION TECHNIQUES
-
MANESH J. SHAH 72
INDUSTRIAL A N D ENGINEERING CHEMISTRY
whereas in others, capacity of one of the synthesis gas production units may impose restrictions in computer control. In what follows, however, the synthesis reactor is treated as a separate unit with the object of understanding its behavior and obtaining the key variables that will lead to its stable, sustained optimum operation. T o describe the reactor operating conditions as accurately as possible, the simulation model should take into consideration all the physical and chemical processes taking place in the reactor. I n doing this, unfortunately, we develop an extremely complex mathematical model whose solution presents forrnidable difficulties even on a large computer. The purpose of our simulation is not to replicate exactly the temperature and concentration profiles in the reactor (presuming that measurements of these profiles are available from an industrial reactor). We are interested mainly in predicting the reactor behavior when changes are made in the controllable variables of the reactor and specifically in studying the variation of ammonia yield as a result of the changes. Thus, a mathematical model which approximates the reactor to the point that it predicts the trends of the reactor output and reactor stability with reasonable accuracy will be adequate for our simulation. Simplifying assumptions can then be made in deriving the simulation equations to the point that the solution of the equations can be carried out in reasonable time on a large digital computer. The Synthesis Reactor
Ammonia is formed from synthesis gas consisting of Hs, N'L,and inerts (CH4 and argon) at elevated pressures (200 to 300 atm.) and temperatures (400' to 500" C.).
Since the reaction is exothermic, the reaction is carried out in various types of autothermic reactors or quench type reactors (79). We will consider in this paper an autothermic reactor with multiple catalytic beds. A
simplified diagram of the reactor is shown in Figure 2. We see that the reactor serves as its own heat exchanger to heat the incoming synthesis gas. The energy balance in the reactor is obtained by equating the heat of reaction released to the total sensible heat gain of reactants and products. As van Heerden (7) points out, stable operation of an autothermic reactor such as the one used in ammonia synthesis is quite critical. The stable operating point of the reactor depends on variables such as feed composition, feed temperature, heat exchange rate, catalyst activity, pressure, etc. This dependence is illustrated in Figure 3, in which the rate of heat generation and the rate of heat consumption by the gases to reach the reaction temperature are plotted as functions of temperature. The S-shaped curve represents the rate of heat generation due to the reaction while the family of curves represents rate of heat consumption as sensible heat and heat exchange. At points 1, 2 and 3, the heat produced equals heat removed. Point 1 is trivial, point 2 is a possible operating point, but is inherently unstable. Slack, Allgood, and Maune (75) have discussed the practical difficulty of maintaining a reactor at point 2. Point 3 is a stable steady-state operating point representing the highest percentage conversion of reactants. The difficulty of operating at point 2 of Figure 3 can be seen by considering the effect of a small perturbation from this point. For a small negative change in temperature, the heat production rate will be less than the heat consumption rate and the temperature will continue to fall unless proper corrective control action is taken. For a positive perturbation in temperature, the opposite effect is obtained as indicated by the arrow. Thus, point 2, while mathematically feasible, is a semistable point of reactor operation. Point 3 is an inherently stable point since a perturbation there will not lead to the trivial point 1, as indicated
VOL. 5 9
NO. 1 J A N U A R Y 1 9 6 7 73
Figure 7 . Flow dragram of an ammonia proccss
in Figure 3, and, in addition, it yields higher conversion to ammonia than point 2. Control of an autothermic reactor such as the ammonia reactor will involve manipulation of control variables in such a way as to keep the reactor at p i n t 3. As we shall see later, the optimum conversion in the reactor is obtained where points 2 and 3 coincide, as illustrated by p i n t a in Figure 3. If the feed rate is increased or the inert content of the feed gas increases, the heat production curves shift to the right. The feed rate may increase to the p i n t that the heat consumption curve becomes tangent to the heat production curve as at point b. This is known as the “blow-out’’ condition. The combination of low feed rate and high heat consumption rate may also produce a blow-out condition as shown by point a. Logeais (72), Eymery (6), Baddou r(Z), and van Heerden (7) indicate that in actual practice, for TVA-type reactors, the desirable stable point for maximum ammonia conversion is in the vicinity of the blow-out point. I t is clear from the above discussion that there is a definite range of feed rates for reactor operation. It is also reasonable to expect that the effects of changing over variables such as feed composition, pressure, etc., will be quite similar to the effect of changing temperature and feed rates. Since the optimum stable performance of the reactor is critically dependent on the values of these variables, the ammonia reactor is an ideal candidate for digital computer control. If one can acquire knowledge of the operating characteristics of the reactor, digital computer control can assist in obtaining optimum reactor performance continuously. Model Devdopmenl: Di(hnnHal Equation. and Eoundary Conditions
A simulation model for reactors of the type considered here consists of different equations describing mass, AUTHOR Manesh J . Shoh is a Senior Member of Control Systm Development Center of thc Data Processing Division, Znternational Business Machines Corp., Sun Jose, Calq. Thc author achwlcdges thc assistrmcc of Prof. J . H . Dufin and M I . Chcster James and thc suggestions of Dr. J . C. Fritz. Mr. E. Fenton, MI. D . Wyman, and Mr. W.M . Synpovidcd programming assistance. 74
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
energy, and momentum transport within the reactor. The most general equations would include time as a variable, so that both the dynamic and steady-state behavior of the system can be investigated. In the study of ammonia reactor dynamics, Eymery (6) encountered formidable difficulties in the solution of the resultant partial differential equations, so that the equations had to be linearized. The linearization was found by Eymery (6) to be valid only for temperature variation of f 5” C. because of the sensitive reactor behavior. The residence time for the gases in the reactor is of the order of 2 to 5 seconds so that the dynamics in gas phase are fairly fast. It is true that the dynamics of the solid catalyst will be governed by heat transfer, diffusion, and chemical reaction rates between the solid and gas. However, a simulation model which must be used in optimization scheme of the plant must be simple because of the limited capabilities of on-line digital computers. Thus, it is impractical to consider in a simulation model the dynamics of the catalyst. The recycle effects in a quench-type reactor (see Figure 2) considered in this paper will probably dampen the effects of the lag in solid dynamics and reduce the errors that may be encountered in ignoring the dynamics of the reactor. In any case, we will only consider the steady-state model, which, as we shall show later, can be used in a steady-state optimization scheme, off-lime as well as on-line (depending on the frequency of plant disturbances). The model has also been used to obtain information on interactions between the variables of an ammonia reactor. Figure 2 shows a simplified diagram of a typical multiple-bed quench-type ammonia reactor consisting of several sections-viz., a series of catalyst beds and a heat exchanger. The feed is divided into several parts; one portion goes to the heat exchanger and the others are divided into several “cold shots” for mixing with the gases between consecutive catalyst beds. Reactor designs of Uhde, Kellogg, and other companies (79) utilize more than two beds with direct or indirect cooling of reaction mixture between consecutive beds. The analysis presented in this paper can easily be extended to such types of reactors but for illustration purposes we will discuss only two-bed reactors.
w e WUI maxe cer model development in manageable form. Some of the assumptions are similar to those used by van Heerden (7). However, unlike other workers (2, 3, 6, 72), we will consider the nonideal behavior of the gases not only in the reaction rate expressions but also in the energy balance equation, where the specific heats and heat of reaction. are functions of tempercture and pressure. As we shall see in the solution of equations, at operational pressures and temperatures in the reactor, these nonidealities have an important bearing on the model solutions. The design of the converter is such that a thin layer of gas surrounds the outside of the catalyst bed wall. The gas acts as an insulator so that the bed is very nearly adiabatic. One can therefore assume that the radial temperature gradients within the bed are small. Furthermore, the ratio of bed diameter to the catalyst pellet diameter is so large that the velocity profile is essentially uniform along the bed diameter. Since, in simulation model, one is interested only in the variation of overall exit gas composition and temperature (integrated over radius), the assumption of radial uniformity in temperature, velocity, and concentration will be adequate for our steady-state model. For the reaction
‘/aNai’ / a Ha e “a a mass balance for hydrogen is taken over a differential reactor section leading to
Fidxi = riUd2
(1)
where the subscript 1 refers to hydrogen. The energy balance equation is obtained by equating the heat of reaction term to the sum of the sensible heat gain for the gases yielding
[F
nt(C,(T,p))tldT = -a/a[Affn(T,p)laFidxi
(2)
Note that the heat of reaction term is for conversion of 1 mole of Ha;n, in the left hand term in Equation 2 indicates quantity of each component and is a function of xi and hence Z. Both C, and A“, are functions of temperature and pressure. Equations 1 and 2 also hold for subsequent beds; we have avoided further subscripting Fi, n,, and xi for each bed to avoid confusion in notation. However, it must be remembered that these three quantities are different in each catalyst bed. In the heat-exchanger section of the reactor, no chemical reaction takes place and it consists simply of tubes with reactor gases on the inside flowing essentially countercurrent to feed gases on the outside. The energy balance equation for the gases inside the tubes is
[Ent’(C,(p,T)Jt]dT’ = -U’(T’ i
The term n’ in Equation 3 is independent of tube lengthz‘and, infact,foratwo-bedreactornt’ = nt(La)i.e., the gas composition at the exit of the second reactor bed. The energy balance equation for the gas on the outside (or the shell side) of the heat-exchanger tubes depends on the design of the reactor. In many ammonia reactors, as shown in Figure 2, the design is such that only part of the incoming cold gas feed is preheated in the heat-exchanger section, and the remainder bypassed. The bypassed gas is divided into several cold shots to quench the hot effluent gases from catalyst beds before entry to subsequent catalyst beds. If 6 is the fractional amount of feed gas entering the shell side of the heat exchanger, the energy balance for the gas on the shell side is
a[xF,,(C&,T))c]dT,’
= U’(T,’
- T’) O‘dd
(5)
- T#’)u’dz’ (3)
where the primes are used to denote the variables in the heat-exchanger section. T’is the gas temperature in the tubes, and To’ is the gas temperature in the shell ride. The heat-transfer coefficient U’may be expressed in terms of the height of a heat-transfer unit as
f tonpno-
TABLE I. HEAT RELATIONSHIPS FOR HEAT CAPACITIES OF VARIOUS GASESa
I
Gas
Hydrogen Nitrogen Ammonia
Methane
1
Argon
I
0.03875 X 10-9 Ta
Heat Capacity Equation
+ + + +
( C P ) l = 6.952 - 0.04576 X 10-ZT 0.09563 X 10-6T2 - 0.2079 X 10-QT3 (Cp).2 = 6.903 - 0.03753 X 10-2T 0.1930 X 10-6T2 - 0.6861 X 10-9T3 ( C p )= ~ 6.5846 - 0.61251 X 10-2T 0.23663 X 10-5T2 - 1.5981 X l O - Q T 3 [96.1678 0.067571~ (-0.2225 1.6847 X 10-4p)T (1.289 X IOF4 - 1.0095 X 10-7P)T21 (Cp)d = 4.750 1.200 X 10-2T f 0.3030 X 10-6 T2 - 2,630 X 10-9T3 (Cp)6 = 4.9675
+ +
+
+
K,(T,p)
=
+ 23.050]}
(1.7343 - 8.143 X 10-4p)
(5.714 X lo-’#
- 2.6714 X
Fi(1 -
nl
I n addition, since the gas stream exiting from each bed is mixed with cold shots, an expression is needed for the temperature of the mixed gas leaving the mixing zone to provide the boundary condition for the subsequent catalyst reactor sections. Thus, for a two-bed reactor
(7)
Equations 1, 2, 3, and 5 constitute the expressions for the simulation model. It still remains to obtain relationships for reaction rate, heat of reaction, heat capacities, and stoichiometric relations between n i and xl, before proceeding with the solution. The reaction rate is expressed by the modified TemkinPyzhev equation (78) =
3fk, 9.8068 X lo-‘ X
where k,, K , and K, are respectively the rate of reverse reaction, equilibrium constant (ideal gas) and fugacity correction term. These three terms are functions of temperature and pressure and were expressed as below using data of Annable ( 7 ) , Hougen (8) and other sources (79).
300 (I) exp[-24,092.2/T + 33.55661 0.63
kT =
K
=
exp(UT) - AFO
=
exp{i.
(9)
[-9184.0 -
7.2944 In T 76
+ 0.34996 X lo2 T + 0.016781 x
T2-
I N D U S T R I A L A N D E N G I N E E R I N G CHEMISTRY
T2 (11)
-
(12)
x1
-
3
F1
-
2
3 xiFi
(15)
where Fi = 6 ( F d , ) , since only part of the feed goes to the first bed. The expressions for n i for the second bed follow in the same manner, with appropriate change in x l ( 0 ) and F , = Fdi in the second bed. I t has been indicated (27) that catalyst activity may change through the reactor as much as 15 to 20% as the reaction mixture proceeds from top to bottom of the catalyst tube. A catalyst activity factor f is incorporated in Equation 8 not only to take the decay of catalyst activity into account but also to incorporate the change of catalyst activity along the reactor length; f was assumed to be a linear function of z. Next, an expression for pressure drop within the reactor is required. Typical pressure drop calculations for this type of reactor using data from references (4, 5, 70, and 79) indicated that the A p / p values were small. Therefore, p , is expressed in a linear form
p = p o - wz The values of heat capacities (CP)% and heat of reaction (AHRP).I are dependent on temperature and pressure. The specific heats are expressed by Hougen and Watson (8) as polynomials in T; the pressure dependence of (C,) was obtained by expressing the coefficients of the polynomials as linear functions of p from data at two pressures (200 and 600 atm.). Table 1 lists the heat capacity relationships in c.g.s. units for various gases; the pressure dependency for NP, HP, CH4, and Ar is negligible. The expression for heat of reaction was obtained in a similar manner using data at pressures of 200 and 1,000 atm. to determine the pressure dependency of coefficients. The resulting expression is (AHEP)3 = -91,840
T
+
XI)
Ftotai
+ (1 - 6 ) Fdifor all i
T
The details of the derivation of these equations have already been described in another publication by Shah, Duffin, and James (74). Stoichiometric and material balance relations lead to the expressions for nt as follows :
ne = Fz
[ni (0112 = [ n G ) 1 1
+
2.0 X
Temperature Range 500’-900” K Pressure Range 200-1000 atm.
Where TI, T2, and the subscripts 1 and 2 outside the large brackets refer to the gas temperatures and cornpositions in the catalyst bed for first and second reactors, respectively. T, ’(I,’)incidentally is the same temperature as the temperature for the feed gas. Also
(10)
10-2T2
- 7.2949T
+ 0.03356 X
+ 0.34996 X
10-5T3 - 0.11625 X 10-9T4 -
+ (14.3595 + 4.4552 X 104p) X T - T’ (8.3395 X lo4 + 1.928 X 104)) 51.21 + 0.14215p (17) (6329.3 - 3.1619p)
r
.
I
RfAO l N l l l A l CON011101~,fARAYtltlS
where again, T is in O K., p is in atmospheres and AH is in gram calories/gram mole ammonia. The simulation model thus consists of Equations 1 through 17. Equations 1 and 2 are repeated for the consecutive catalyst beds; however, the initial temperature, pressure, and composition for each bed (hence the initial boundary values) will be different for each bed. Numerical Solution
The simulation equations are coupled, and nonlinear, and must be solved numerically. The boundary conditions for a two-bed reactor are
TdO) = T.’(O)
(18)
Fzgure 4. Block dtogrmn of computer logrc, Method 1
for the first catalyst bed
(19) TdO) = f[T,’(L‘), Ti(L3, 61 for the second catalyst bed, as given by Equation 6 T’(0) = TdLd
(20)
for the heat-exchanger section and
To‘@’) = T(0)
(21)
which is the feed gas temperature (see Figure 2), for the heat-exchanger section. The flow rates are
b,(0)11 = 6
’
Fa.
(22)
for first bed and
[nc(O)h = nc(L1) -k (1.0
- 6)Fd.
(23)
for second bed. It is obvious that thia is a two-point boundary value problem, because the temperature at the entrance of the h s t bed, Tl(0)(Equation 18) is not known and must be assumed to start the numerical solution. If at the end of heat-exchanger calculations we find Equation 21 satisfied, the assumed value of Tl(0) is c o m t . Otherwise, a new value must be assumed and the solution repeated until Equation 21 is satihed. We a h note that when the assumed Tl(0) does not satisfy Equation 21, the :valuation of Tt(0) using T(0) instead of T,’(L’) is in m r , since T(O)#T,’(L’) in this case. yowever, if there is no great difficulty in convegence to the wrrect T1(0),one may be able to ignore this initial discrepancy. 3 n the other hand, we may assume T’(L’) and integrate >a&ards in the heat-exchanger section to obtain Tl(0) = T,’(O).Then, we can carry out the forward ntegration on the catalyst section and test for T&) = T‘(0). The equations thus may be solved using this .est and guessing T‘(L’). We note, however, that for ntegration in the heat-exchanger section, we assume d composition of the gas stream emerging from the second bed. This difficulty can be eliminated by successively replacing the assumed conversion with calculated conversions at the end of each iteration.
I n the present work, both of the above mentioned methods for solution were used. The use of a modified rcgulofdsi technique (9) of convergence helped to get the c o m t value of Tl(0)or T’(L’) in four or five iterations in almost every case. The solution to the simulation equations was carried out on an IBM 7090 using a DSL/90 language for programming. The details of the language and its usage may be found elsewhere (77). The program was divided into blocks as shown in Figure 4. We will describe it briefly here. In the first step, all calculations that have to be performed only once are carried out; C,’(T(O))and U’are calculated in this block. I n the second step, the integrations for the catalyst section are carried out with an initial guess of T,(O)= Tl(O),using a Milne predictorcorrector integration method for evaluating x1 and Tl. In the third step, the mixing Equations 6 and 7 are used for calculating the temperature and composition at the entrance of the second bed, and I and L are set equal to z m and Ls, respectively. The integrations in the second catalyst bed are carried out next. I n the fourth step, the T‘(0) is set equal to TI(&) and the heat-exchanger integrations are carried out. I n the fifth step, the value of T,’(L’) is compared with the value of feed temperature (called T(0)) and if the difference is larger than one degree, the modiied r e d o falsi technique provides a new estimate of T,(O)= Tl(0), and the program reverts to the second step. The Milne predictor-corrector integration routine used here for sohtion of the coupled dflerential equations chose its own step size, &, subject to the limitation of absolute error or relative error specified. Errom of less than 0.01’ K. in temperatures and 1 X 1 0 4 in fractional hydrogen conversion were allowed in o w calculations. For all reactor conditions for which a solution was obtained, two values of Ts(0) = T,(O) satisfied the required boundary conditions indicating two operating points of the reactor. V O L 5 9 NO. I
J A N U A R Y 1 9 6 7 77
I n the second method of solution as shown in Figure 5, the heat exchanger calculations are performed in step 2, M o r e the integrations are carried out in catalyst beds. I n the heat exchanger since To'@')is known and if one more temperature as well as conversion at the end of the second catalyst bed are wumed, the heat-exchanger calculations yield the other two unknown stream tempera-. It was observed that if average values of C, were used in Equations 3 and 5, the differential equations could be solved analytically as follows. If C,(T') and and C,(T,') are constant in Equations 3 and 5, Equation 3 can be integrated to yield :
T,'(O) =
. U'u'L' - S / a ) -+ T'(0)[l U' . (I'. L' [(@la) - e(l-81-) .
- T,'(L') . (1
B
(24)
I
,
'"Y,.unW
,
5
a =
C nr(La) . Cpd(T')ev ,-1
(25)
and
.c 5
B=8
i-1
Fd,
. C,'(T,')..
We note that To'@') is given and if T'(0) is assumed, Equation 24 enables us to calculate To'(0) which is equal to Tl(0). Equation 24, however, is an implicit equation since a and are also dependent on T,'(O) and T'(L') which are not known. It can be shown that T'(L') is given by
T'(L') = T'(0)
+ -aB (T,'(L')- To(0))
(26)
A trial and umr solution is involved in Equation 24 for which a simple successive approximation technique is adequate. The temperatures calculated by Equations
,
Figure 5. Block dagrm of cmnpufar logic, Method 2
Figwe 6. Tanpnohrra and cmwersion poofrs in reac&w bed 1
70
where
INDUSTRIAL AND ENGINEERING CHEMISTRY
Figure 7. Ammoma produc(imt and hydrogen dcpkrion pq rmtw bad 1
____ ..~
Figwe 8. Effectof vm'afimt in 6 m reactwpnfotmmtcs
24 and 26 agreed to within
of the temperatures resulting from the numerical solution of Equations 3 and I/*'%
5.
Discusdon of RaruN.
!+$pi" .> +.. r In fligures 6 and 7, we show temperature T,fraction
5), then, we assume T'(O)](which is the same as Tz(L31
H,conversion, and H g and NH,flow rate profiles in the reactor bed for a typical set of reactor conditions. Al-
and the composition of the gases coming out of the second catalyst bed and calculate T,'(O) [which is equal to T1(0)].In steps 3 and 4, the integrations of material and energy balance equations for the catalyst bed are carried out in the same manner as in the first method of solution. At the end of calculations in the second bed one checks if T&) is equal to the assumed value of T'(0). If the dserence IT&) - T'(0)I is larger than a specified value e, the calculations are repeated assuming a new value of T'(O), and replacing the assumed gas composition entering the heat exchanger with the calculated values exiting from the second bed. Thua, we revert to step 2 and repeat the integrations until IT&) TASl is less than the specified error limit.
though no comparhn of the profiles with actual reactor bed temperature8 has been made to date, the temperatures and composition of the gases agreed within 3% with the values between beds and at points where the measurements were available. These values were provided by the plant designer. We note also that equilibrium conditions are reached at the end of the firstbed. As indicated earlier, the p u v of this work is not to replicate exactly any particular reactor profiles; rather, it is to observe the performance of the reactor as the controllable and uncontrollable variables change their values. Several computer "UNI were, therefore, made varying some of the parameters such as feed temperature, pressure, composition, flow rate, feed split 8,
In step 2 of the second method of solution [see Figure
-
Figure 9. Effecfof prersure change 01
psion
-T
Figurr 10. Effat of um'afi'on i n inbf f m p m l w e on conversion
Figure 12. Ralafim.rhip befwca rsoctor skabilify and inkf feed ternpsrohrra V O L 5 9 NO. 1 J A N U A R Y 1 9 6 7
79
and heat transfer coefficient (directly related to tube fouling). I n Table 11, we have shown the values of these parameters for the runs and the resulting reactor exit gas ammonia production and temperature. The former is of importance in plant ammonia production whereas the latter is an indication of the heat recovery possible from the exit gas stream. For each set of parameter values, two solutions to model equations are indicated by two values of assumed T ~ ( L z )which satisfy the boundary conditions. The run numbers l a , 2a, etc., represent the condition of the upper stable point which is the desirable operating point for the reactor. Our computer results indicated that the reactor is not only bistable, but that beyond certain critical values of the various parameter combinations, no solution of the model equations could be obtained. The solution to the equations does not exist when the heat of reaction does not match the energy required to raise the feed gases to reaction temperatures, so that the reactor is unstable under these conditions, being outside the range determined by quench and blowing points shown in Figure 3. An illustration of the instability is indicated in Figures 8 through 12, where we have plotted H2 conversion and T A S on the ordinate and the various control parameters on the abscissa. Thus, if feed split is decreased from 1.0 to 0.75, the difference between the two assumed values of T2(L2) (called T A S ) decreases. I t decreases sharply from 0.75 to 0.722 and the decrease is not linear; at 6 = 0.72, there is no solution possible. Similar situation arises when the inerts are increased from 9.6 to 18%, the pressure decreased from 240 to 160 atm., feed increased and T ( 0 ) decreased. From a physical standpoint, the instability is indeed reasonable since, for example, an increase in inerts only wastes the heat of reaction in heating inerts together with the NBand Hz to the reaction temperature. As a result, beyond a certain inert concentration, the heat of reaction will not be sufficient to raise the reaction mixture to the reaction temperature. Thus, the effect of varying the parameters will not only show up in the hydrogen conversion but on the reactor stability as well. I t is also seen that the maximum ammonia is obtained closest to the point of reactor instability for T(O), and feed rate. When the total inerts are increased from 9.6 to 16.670 and higher, the ammonia production of the reactor is affected only slightly. However, the stability of the reactor is affected drastically. It was observed, however, that if the H2, Nz, NH3, and Ar feeds were maintained constant, and the methane rate were increased (see runs 1, 19, and ZO), ammonia production improved. These runs show the effect of increasing methane to argon ratio on the reactor behavior. At first glance this contradictory behavior appears to be an effect of increased space velocity. However, a closer look at the results reveals that the increase in ammonia conversion is due to a drop in Tl(O),with increased methane in feed. The decrease in temperature of gas entering the first catalyst bed, due to equilibrium considerations, gives higher conversion. Needless to say, as methane rate 80
INDUSTRIAL AND E N G I N E E R I N G C H E M I S T R Y
is increased beyond a certain point, the reactor will become unstable due to high percentage of total inerts. The effect of increase or decrease in per cent ammonia in feed is shown in the runs 23 and 24. An increase of 1.270 ammonia brings about decrease in conversion of about 1.570 and vice versa. The profit function for the ammonia plant will depend on the cost of increased refrigeration load to keep the ammonia in feed low so as to gain increased production. We see that the increase of ammonia in feed, moreover, leads towards instability of the reactor. Finally, the effect of variation of the (H2/N2) ratio on conversion of ammonia is shown in runs 1, 21, and 22. The ammonia production seems to vary only slightly in spite of the rather severe change in the ratio. In older design plants utilizing 300-atm. pressure for synthesis, the (Hz/N2) ratio was considered critical. I n the low pressure plants it appears that this ratio is not critical in the reactor. Although it is desirable to hold the ratio of (H2/N2) in gases entering the synthesis loop at first stage of compressor to 3/1 from final stoichiometric considerations, it appears that any variation in this ratio will not have any direct influence on conversion in synthesis reactor. There is an indirect effect, however, caused by gradual changes of pressure in the synthesis loop when the ratio of (H2/N2) in feed stream to synthesis loop is higher or lower than 3 for a long time. Logeais (72) and Eymery (75) indicate that the TVA reactor is difficult to maintain stable at the lower stability point. The upper stable point on the other hand, is closer to the critical blow-out point of the reactor in face of varying reactor conditions. For this reason, in normal practice, the temperature at the entrance of the first bed is kept 20" to 25" C. from the blow-out temperature or the temperature at the exit of the second bed is monitored closely. Our results indicate, however, that in spite of the safe temperature margin, an increase in the per cent inerts or decrease in 6 may also lead to reactor instability. Furthermore, Table I1 also indicates that the difference in ammonia production between optimum reactor performance and performance at one of the lower stable points can be considerable. For example, the mole fraction of ammonia in runs 1 and l a suggest that a l5Y0 increase in conversion can be achieved by operating at the upper stable point. Thus, the economic incentive for operating the reactor near the optimum conversion point by manipulating the variables affecting reactor performance and stability using a digital computer control is quite large. Conclusions
I t is evident, from previous discussion, that one cannot attain optimum stable reactor performance by monitoring feed temperature alone as in the normal industry practice. T o use our simulation results as a guide for on-line control, however, some additional measurements are necessary: composition of feed gases at reactor inlet, flow rates of feed and cold bypass, and feed temperature at reactor inlet and inlet to the catalyst beds. Depending
version, for the existing thermodynamic conditions, the heat-produced curve will drop off sharply as the temperature increases. This effect is shown in Figure 10. The position of the heat-consumed curve will be a function of the heat transfer and the feed temperature which, in turn, are functions of the controlled variables. The optimum point of operation would be the “ 0 ” the limit of stability. The position and shape of the heat-produced line are a function of feed rate, reaction kinetics, per cent inerts, etc. The simulation model was
upon a particular reactor design, the controllable variables will be feed flow rate, per cent inerts (indirectly controlled in the purge section), feed temperature (controlled in heat exchangers prior to reactor), and other variables characteristic of a particular industrial unit. These variables determine the relative positions of the heat production and heat consumption lines in Figure 3 and, consequently, the positions of the blow-out and quench points, In an actual production unit operating at a high temperature and near the equilibrium conTABLE 11.
RESULTS FROM SIMULATION MODEL “2
Run No.
P(O), Atm.
1 la 2 2a 3 3a 4 4a 5 5a 6 6a 7 7a 8 8a 9 9a
200 200 200 200 200 200 200 200 200 200 220
10 10a 11 lla 12 12a 13 13a 14 14a 15 15a 16 16a 17 17a 18 18a
Produced G. Moles/
Tz(L2) = TAS
T’(L‘)
XH 2
Sec .
744.6 806.8 740.4 794.6 772.8 772.7 754.15 814.5 760.9 820.2 735.7 816.8 729.0 823.5 755 2 793.9 765.9 781.8
551.8 576.9 556.9 578.6 559.5 572.1 545.9 568.1 548.8 560.5 518.8 580.9 545.9 583.3 556.5 571.6 560.3 566.6
0.1770 0.2050 0.1848 0.2084 0.1894 0.2032 0.1649 0.1930 0.1544 0.1835 0.1747 0.2087 0.1700 0.2130 0.1830 0.1994 0.1885 0.1952
386.7 446.9 403.8 455.5 414.0 444.0 360.3 421.8 337.5 401 . O 380.0 456.0 371.4 465.6 400.0 435.9 412.0 426.6
751.7 801.4 757.4 797.7 780.9 781.5 737.3 809.0 756.5 801 .O 760.5 801.2 770.2 799.8 734.0 818.7 756.8 794.6
547.9 567.4 544.5 560.2 542.2 542.3 548.2 578.0 557.4 575.0 558.3 574.1 560.6 573.6 547.0 580.3 557.8 572.8
0.1887 0.2097 0.1961 0.2130 0.2167 0.2168 0.1738 0.2037 0.1877 0.2061 0.1857 0.2044 0.1912 0.2028 0.1650 0.1961 0.1950 0.2095
412.36 458.3 428.5 465.4 473.6 473.8 341.8 400.7 451.3 495.5 487.0 536.1 543.0 576.0 378.4 449.7 412.1 442.8
No Lution possible when Inerts increased to 18 .5%b 19 678.2 Methane increased from 200 423.0 0.8 751.5 19a 714.6 9.9% t o 1 1 . 3 % 799.0 20 685.1 200 423.0 0.8 Methane increased to 760.6 20a 707.6 12.9% 790.3 21 200 666.9 423.0 0.8 Hz/N2 = 2 . 0 Total 736.6 21 a 717.5 feed rate same as 1 802.0 22 0.8 200 682.0 H ~ / N z= 4 . 0 423.0 756.0 22a 716.7 801 .O 23 423.0 0.8 200 692.3 Ammonia in feed incsd. 769.4 23a 707.0 f r o m 2 . 7 to 3.9oj, 566.8 24 0.8 200 661 . 2 423.0 Ammonia in feed de729.7 24a 726.4 creased to l . 9% 814.2 a Inerts 13.75%. b H,,Nz decrekied or increased to keep total feed rate same as in Run No. 9 .
555.2 574.0 559.5 570.8 549.0 574.8 556.3 574.1 559.8 566.8 547.8 581.7
0,1863 0.2090 0.1967 0.2117 0.1966 0.23 0.1712 0.1904 0.1925 0.2001 0.1699 0,2074
407.1 456.7 429.8 462.6 381.3 445.1 398.5 444.0 416.5 432.7 367.0 453.2
423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0 423.0
0.8 0.8 0.75 0.75 0.73 0.73 0.9 0.9 1. o 1 .o 0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8
200
410.0
0.8
200
400.0
0.8
200
380.0
0.8
200
423.0
0.8
200
423.0
0.8
200
423.0
0.8
200
423.0
0.8
200
423.0
0.8
200
423.0
0.8
240 180 170
Feed Rate
6-syit
673.05 721 23 678.5 721.5 683.5 707.8 664.15 707.7 653.5 693.7 666.4 728.9 660.9 733.7 681.5 711.3 689.5 701.9
Normala Normala Normal Normal Normal Normal Normal Normal Normal Normal Normal Normal Normal
I
Normal
No solution at P(0) = 160 atm. 676.2 Normal 714.4 678.6 Normal 709.5 Normal 692.8 693.0 667.1 Reduced by 10% 723.6 Increased by 10% 682.4 716.7 Increased by 20Oj, 685.3 716.3 Increased by 30y0 692.1 715.6 Inerts reduced to 9,6%b 664.0 730.6 Inerts increased to 682.6 711.6 16.6170 to 18.5%b
I
1
!
VOL. 5 9
NO. 1
JANUARY 1967
81
first used to define the operating stability limits of these variables and was then used in conjunction with an optimization program, where the variables yielding maximum throughput were determined. I n general, it is best to operate a converter near its stability limit-point “0” of Figure 11-to obtain low temperatures and high conversions. The price paid for this type of operation is a decrease in system stability. Converter control by variable manipulation should be such that operating conditions keep the converter near the stability limit with enough safety margin to prevent blow-out or quench. The simulation model has been of great help in defining these operating conditions and their limits and in devising control schemes for maintaining these conditions at optimum values. The reactor model was incorporated in a model of the entire ammonia synthesis loop which includes several heat exchangers and ammonia recovery system. The synthesis loop model was then used in an optimization program (76) to determine not only the values of independently controllable parameters within the synthesis loop but also the feed required in the synthesis loop to give maximum ammonia production and heat recovery within the given constraints of the various units of the synthesis loop. Some of the interactions of the different variables which are not evident from a one-variableat-a-time simulation were found to play important role in the optimization of the synthesis loop (20). The control strategy for the ammonia plant based on our simulation and optimization work in the synthesis as well as reforming section is still under development. However, the basic principle is as follows. If the plant is production limited, determine the conditions for optimum ammonia production and heat recovery in the synthesis section. The makeup feed required in the synthesis loop determined from this calculation leads to requirement of natural gas, steam, air, etc., in the “front end” of ammonia plant. The optimum performance of the front end is carried out next centered in these requirements and an infeasible condition in the front end will lead only to a change in the constraints (upper or lower limit) on the makeup stream. In an operational plant, on-line optimization calculation may be carried out, say, every ‘ / z to 1 hour depending on the frequency of the plant disturbances. Although the optimization is on a steady-state model, it is assumed that the plant is steady at least for the time of optimization calculations, which determine the new values of set points on the analog controller (or for direct digital control algorithms, if these are replacing the analog controllers). The problem of how best and how fast the set points be changed is, of course, determined by the dynamic characteristics of each control loop. If the dynamics of the process can be approximated by a second-order system with dead time, one can use “bangbang” control technique as suggested by Latour ( 7 7 ) to move the plant from one point to another. This method of moving the control valve from upper saturation point to lower saturation point at calculated time intervals is shown to give optimal control based on 82
INDUSTRIAL A N D ENGINEERING CHEMISTRY
Pontryagin’s maximum principle. However, it is suspected that plant operators are unlikely to favor violent upsetting of the control valves in this fashion and experience indicates that a trial and error procedure as suggested in a recent program ( 7 7 ) to achieve the change in set point will be more acceptable to the plant personnel. I n this method, time interval At is determined for the total set-point change and then, for example, 50y0 of the set point change is made at the end of optimization calculations, 75y0 change at time At/2 from the start and - 25% change at At from the start. Summarizing simulation results and their discussion, some recommendations for optimum reactor performance may be listed as follows: 1. The inert gas in feed should be at a safety limit, but methane/argon ratio can be as high as possible. 2. Keep the feed temperature as low as possible. 3. Adjust the feed rate, feed temperature, and other parameters so that the reactor remains within the stability limits, and yet gives optimum ammonia production as determined by simulation. When the simulation results are used for a particular plant, the reactor model must be reduced to reIate the ammonia yield to controllable and uncontrollable (but measurable) variables listed above, if the model is to be used for on-line optimization calculation. For our plant this was carried out by using a multiple stepwise regression technique on the simulation results to yield algebraic expressions from the differential equations. These expressions (generally nonlinear) were then used in an on-line control and optimization routine for a small digital computer (IBM 1800) normally employed in process control. The off-line optimization ( I S ) on a large computer (IBM 7090) with the original differential equation model determined the best operating conditions of the plant (far from design data). The online optimization program with the simplified model on the other hand presumes that the plant does not drift too far from the optimum over a period of 1/2 to 1 hour and hence, not only does the validity of the simple model hold, but also the optimization calculations with the simplified model can be performed in a very short computer time. NOMENCLATURE = cross-sectional area of reactor bed = heat-transfer area per unit length of heat-exchanger tube (Cn)& = molal heat capacity of component z AF, = standard free energy change = feed rate of component i to synthesis loop Fdl = feed rate of component i to the catalyst bed F, = catalyst activity f = height of a heat transfer unit defined by Equation 4, H‘ heat-exchanger section (AHRD)..= heat of reaction at any pressure p based on any component 2 . (Sote: z here is limited to one of components taking part in reaction) K = equilibrium constant KV = fugacity correction factor = reverse reaction constant kr L, = length of reactor bed i = length of heat-exchanger tubes L’ a a’
(3) Brian, P. L. T., Baddour, R. F., Eymery, J. P., Ibid., p. 297. (4) “Chemical Engineers’ Handbook,” J. H . Perry, ed., 4th ed., sec. 5, p. 51, McGraw-Hill, New York, 1963. (5) Comings, E., “High Pressure Technology,” p. 280, McGraw-Hill, New York, 1956. (6) Eymery, T. P., M I T Report ESL-R-191 (1964). (7) Heerden, C. van, IND. END.CHEM.45, 1242 (1953). (8) Hougen, 0. A. Watson, K. M., “Chemical Process Principles,” Part 111, p. 886, Wiley, New York, 1962. (9) Kjaer J. “Measurement and Calculation of Temperature and Conversion in Fixed B)ed batalytic Reactors,” Jull, G. Jerrerups, Forlag, Copenhagen, 1958. (10) La idus L “Digital Computation for Chemical Engineers,” McGraw-Hill, New ’19z2. (11) Latour, P. R., Koppel, L. B Coughanowr, D. R., “Time Optimal Control of Chemical Processes for Set Poznt Changes,” 60th A.1.Ch.E. National Meeting, Atlantic City, September 1966. (12) Logeais, B., M.S. thesis, MIT, 1959. (13) Shah, M. J Proc. S mp. on Application of Mathematical Models in Chem. Eng. June Res., 1965. Deegn and Sroduction, A.1.Ch.E.-1.Ch.E. Joint Meeting, London,
flow rate of component i at some point in reactor total system pressure gas constant rate of reaction of component i rz temperature of gas inside catalyst bed at any point T temperature of gas inside heat-exchanger tube at any T‘ point = temperature of gas outside heat-exchanger tubes at any ?Ig’ point = overall heat transfer coefficient defined by Equation 4, U’ heat-exchanger tube W = coefficient in pressure drop = fractional conversion of component i xi z = distance in flow direction within catalyst bed 2‘ = distance in flow direction within heat-exchanger tubes = fugacity coefficient of any component i V, = defined by Equation 25 cy, P Subscripts = component designating index 1-hydrogen, 2- nitrogen, 3-ammonia, 4-methane, 4-argon = number of catalyst bed fii
P R
= = = = = =
cork,
(14) Shah, M. J., Duffin, J. H., James, C., IBM/San Jose, Tech. R e p . TR.02.364 (1965). (15) Slack, A. V., Allgood, N. Y., Maune, H. E., Chem. Eng. Prog. 49,393 (1953). (16,) S$th H V “A Process Optimization Program for Non-Linear Optimization I B h Pro&am Information Dept., Hawthorne, N. Y . , Program No. IBM002i-7090-1965. (17) Syn W. M Wyman D. G “DSL/90 Digital Simulation User’s Guide,” IBM/
