I
S. M. ROBERTS and C. G. LASPE
TRW Computers Co., Beverly Hills, Calif.
On-line Computer Control of Thermal Cracking Computer simulation reveals the shape of temperature, pressure, and conversion profiles along the tube length. From the structure of these profiles, simple conversion equations were developed for com puter control
AN development of
is presented for the conversion expressions for on-line computer control of a thermal cracking reaction. Through the simulation on a digital computer of a tubular reactor, certain structures in the temperature, pressure, and conversion profiles were revealed. Based on this, the conversion can be described adequately for control in terms of effluent temperature and pressure, total flow rate, and tube diameter. The utilization of readily measured variables rather than an entire profile to characterize reactor performance provides a simple and effective control tool. Two conversion expressions were developed, one based on a multiple regression, the other, on theory. Both expressions may be evaluated in seconds from plant measurements and may be used as predictors. APPROACH
Program Design
Reaction.
The reaction considered is
kl
A @ B + C kz
The compound A , initially a liquid, is heated to its boiling point, vaporized, and superheated until it is thermally cracked to compounds B and C. The reaction is an endothermic first order system. Simulation Program. The cracking reaction for the tubular reactor was simulated by a step-by-step calculation on a digital computer. The program computes as a function of tube length 0 0 0 0
Pressure Fluid temperature Wall temperature Conversion
The computer program considers the reactor to be made up of three stages: Stage 1, preheating; Stage 2, vaporization; Stage 3, superheating and thermal cracking.
For all three stages, the heat flux distribution is specified as a step function of tube length. For the operating plant to be simulated, the information for computing heat transfer coefficients, especially for the scale and outside film coefficients, was not available. T o obviate this difficulty, a lumped heat transfer coefficient was arbitrarily chosen to account for the heat resistance of the inside scale, the pipe wall, and the outside film coefficient. Computer results showed that the tube wall temperature was reasonable, based on this lumped heat transfer coefficient. Inside film coefficient and over-all heat transfer coefficient, U, were computed for each element of length. The tube wall temperature was calculated as a function of the heat flux, the fluid temperature, and the over-all heat transfer coefficient. The size of the incremental length A1 was chosen small enough so that a log mean temperature difference driving force was not required. The friction factor, f,was taken as one fourth of the value given by Brown ( I ) . The program takes into account the pressure drop loss across the return bends. All the fluid properties are evaluated at the upstream side of the incremental length. Discussion of Computer Results
A set of computer runs was made to determine the effect of the size of AI and AX,,,, on the speed of convergence and on the accuracy of the calculated conversion. I n general, the smaller the Al, the less sensitive the calculated conversion is to the size of the AXtest. For a fixed size of Al, as the AXtest becomes smaller, the number of iterations required for convergence increases, as would be expected. I n all cases considered, the number of iterations per A1 is small, averaging around two iterations. For many cases, a AI of approximately 3 feet and a AX of about 0.01 mole are reasonable compromises between accuracy and speed.
A set of computer runs was made to determine the Arrhenius constant. Various values of k ~ were , chosen and tested by comparing the resultant calculated conversion, temperature distribution, and over-all pressure drop with the observed plant data. The selected value of ko occurred when the calculated conversion was within 2% of the measured value, when calculated pressure drop across the reactor was within 2 p.s.i. of the measured value, and when there were only minor differences in the temperature profile. Typical computer results are plotted in Figures 1, 2, and 3. To protect proprietary information, the scales of the ordinates have been adjusted. The scale of the abscissa is true so that the trends, shape, and slopes of the various curves have been preserved. From the typical results, certain structures in the solutions are revealed for the various stages in the reactor, and for the temperature, pressure, and conversion profiles. As a general rule, the tube length is allocated so that the preheating section takes about 4oy0 of the tube length, the vaporization and superheating parts of the reaction zone each take about lo%, and the cracking reaction occurs in the remaining 40%. In regard to the temperature profile, the temperature increases linearly with length in the preheating zone. Any slight curvature near the vaporization zone is due to the change in heat flux. I n the two-phase vaporization zone, the temperature drops slightly, perhaps as much as 3' F., while in the superheating portion of the reaction zone, the temperature increases linearly rapidly. When the conversion reaches 2 or 3y0,there is a knee in the temperature profile. From about the 5% conversion point, the temperature is once again linear with length, but with much smaller slope than in the superheating region. T h e structure of the pressure profile shows that the pressure is approximately VOL. 53, NO. 5
M A Y 1961
343
ON-LINE COMPUTER CONTROL
if ;
N
0
z 0 a
(pi
-
a C.
0
~
AI
3 W
(pi
4
c
E
LL
I
L
Lo 0
T
0
z
0 a
&
s
3
t
Y
n
a
Lz I
B
d
z W -1
m
z
IL
w
L
W
I I :
e
0
c @
k W -
a -II 0
l
z a
+
=0
~
3 P
v
(I)
m
344
INDUSTRIAL A N D ENGINEERING CHEMISTRY
w 0
(pi
i
Stage I-Preheating In the preheating stage, the liquid A is heated from its initial temperature and pressure to its boiling point. The process can be described b y a pressure balance and a heat balance. Using finite difference approximations, the equations are: 2fG2 A -P = _ _ A1 pd,g
(1)
relation (2), which is stored as a table in the program. The equations for Stage 2 are quite similar to those of Stage 1. The wall temperature, however,
evaluated on mole fraction basis As derived in Derivation of Kinetic Equations, the combined material balance and kinetic equations yield
is computed b y
The trial and error calculation is made on AX, the quantity o f compound A reacting in the length A/. Briefly, a trial value of AX is substituted into Equation 10, which is solved for AT, and TPz. After evaluating the fluid properties and the pressure drop, Equation 13 is solved for AX. The adequacy of the trial AX is found b y testing whether 1 AXtrial - AXoaloulatedI is less than AX,,. The AXGost is a prescribed number. Details of the calculation are given in the flowsheet. The iteration scheme was chosen originally for its simplicity. It also proved to be extremely rapid. Hand calculations showed that the best choice for the trial value of AX at the beginning of Stage 3 was AX 0. In addition, the best trial value of AX for the next iteration was the previously calculated value of AX. This was true for the iterations at a given A/ and also for the initial trial value of AX for the next A/. Computer Runs. Runs were made to determine: 8 Most efficient size of A/ and Ax,,, o Value of the Arrhenius constant, ko,that fits plant data best o Effect of tube diameter on conversion Efiect of heat flux per pound of reactant on conversion 0 Effect of heat flux distribution at constant total heat flux per pound of reactant on conversion 0 Effect of inlet pressure and flow rateonconversion
- TI (8) and the Equation 2 is replaced ATw = T w
(9)
(3) 1~.- 1
1
(4)
u-h,+x h _ c,G
-_ -
0.023
-
(5)
+ TFZ
(6)
R~O.*PI-~/~
A T w = Tw -
(AP)zr =
TFI
($)
2
lu
(7)
The computational procedure for Stage 1 is straightforward and no trial and error calculation is required. The test for exiting from Stage 1 to Stage 2 is the comparison of the temperature and pressure at the downstream side of A/ with the fluid vapor pressure at the corresponding TF2. See the flowsheet for details.
Stage 2-Vaporization In Stage 2, since the liquid boils as it flows, a two-phase fluid flow and heat transfer problem must be solved. The process can be described by material, pressure, and heat balances. An estimate of the two-phase heat transfer coefficient is taken as three times the coefficient of the liquid, if it were flowing alone. The data of McAdams (3) indicate that this is a reasonable assumption. The two-phase pressure drop is estimated by the Lockhart and Martinelli cor-
The calculations for Stage 2 are straightforward and no trial and error calculations are required. The exit from Stage 2 is made when the quantity of vaporized material exceeds the total feed rate. Suitable interpolation is then made within the program to guarantee a material balance.
Stage 3-Superheating Reaction
and
Entrance into Stage 3 may be made either directly after loading the data or in sequence following Stages 1 and 2. In addition to the material, pressure, and heat balance equations, common to Stages l and 2, Stage 3 calculations include kinetic equations. ln contrast to Stages 1 and 2, the Stage 3 does involve a trial and error calculation. The heat balance, Equation 2, is amended to take into account the heat of reaction
The reaction velocity constants are given by kl = koe
- 1.987 ( T pAE* A V g+ 460)
k. = K ,
(11) (12)
kz
The properties of the mixed gases flowing in Stage 3 are
=
VOL. 53, NO. 5
MAY 1961
345
I II 2.44 2.34 544 544
TEMPERATURE
50
-
200%
REACTOR TUBE LENGTH,
Figure 1.
L
Effect of tube diameter on reactor profiles
I
II
2,480 1 9,000 19,000 22,125 496
2,480 14,900 1 2,800 17,700 446
I
II
111
6,400
6,400
Curve
&?/aAQ, in B.t.u./hr.
5
$500
100 -&400 Q w
300
[
-t
IO0
0
0
I
0
I 250
a 200
O.IL
I 02L
0.3L
04L
0 5L
06i
0.7L
0.8L
0.9L
I.OL
0
1100
1 t
- 10
_c---
-! 50
I100
800 900
LL.
700
5
600
sq. ft. from 0.0001 to
-
0.400L 0.400L to 0.700L 0.700L to 1.OOOL
500
400
6,400
1 8 , 4 0 0 ~1 6 , 4 0 0 2 0 , 4 0 0 18,400 20,400
16,400
Flow rate, Ib./hr. 22,125 22,125 22,125 Conversion, 51 .14 50.30 51 .82 H e a t input, B.t.u./lb. feed 544 544 544
300
7 0
200 100
0
Curve
AQ/AAQ, B.t.u./hr.
c 100
sq. ft. from 0.OOOL to 0.400L 0.400L to 0.8201 0.82OL to 1 .OOOL Flow rate, Ib./hr. Heat inDut, B.t.u./lb. f e e d
0
O.IL
0.2L
0.3L
0.41
0.5L
0.6L
0.7L
OBI
0.9L
1,OL
REICTOR TUBE LENGTH,L
Figure 3.
Effect of heat flux distribution at constant heat flux on reactor profiles
linear with tube length in the preheating zone, the vaporization zone and the superheating zone, with the slopes increasing in each zone. I n the reaction zone, the pressure profile is slightly concave. The conversion distribution is characterized by insignificant amount of conversion until the temperature reaches about 850" F. T h e conversion is curvilinear with length until about 5% conversion after which it becomes linear with length. The effect of coking in the reactor, causing smaller tube diameter, is shown in Figure 1. With smaller tube diameter the pressure drop increases, the temperature in the reaction zone increases, and the conversion drops.
346
The effect of total heat flux per pound of reactant on conversion is illustrated typically in Figure 2. AS might be expected, the higher the total heat flux per pound of reactant, the higher the conversion, with all other conditions constant. The effect of heat flux distribution for constant heat flux per pound of reactant is shown in Figure 3. The highest conversion is obtained (Curve 111) when the largest heat flux is concentrated a t the vaporization zone. T h e lowest conversion is obtained (Curve 11) when the lowest heat flux is concentrated a t the vaporization zone. Uniform heat flux distribution (Curve I) gives intermediate conversion. There-
INDUSTRIAL AND ENGINEERING CHEMISTRY
fore, judicious distribution of thermal energy can yield higher conversion at no extra cost. The information presented here confirmed suspicions of the plant operators and resulted in their changing the firing pattern. The study revealed that higher conversions are favored by lower flow rates and also by higher pressures when the thermal energy per pound of reactant is held constant. Although for this reaction thermodynamics indicate that higher pressures favor lower conversion, the reaction kinetics drive the system toward higher conversion at higher pressures. Whether there is an optimum pressure and flow rate for this reaction has not been investigated as yet.
ON-LINE COMPUTER CONTROL
Development of On-line Computer Control Equation
Fx
on-line computer control of a chemical reactor, it is necessary to compute or measure conversion conveniently and then to adjust the set points of various instruments to bring the existing conversion back to the desired conversion level. To analyze rapidly and accurately the reactor effluent stream may be a difficult instrumentation problem. Conditions such as high pressures and temperatures and corrosive chemicals are not uncommon. In addition, there is the problem of obtaining a representative sample. Despite these difficulties instrument manufacturers, through spectography and chromatography, are often able to produce offthe-shelf items that meet severe industrial environments. When instrumentation is available, it is advantageous to compute conversion from the measured process variables as a second check. When the instrumentation is not available, it is a necessity. I t is desired, therefore, to develop from the plant data or from a computer simulation such as this, a convenient expression relating conversion to the process variables. Since the performance of the reactor is characterized by pressure, temperature, and conversion profiles (a distributed system), there is difficulty in “capsulizing” the distributions. Some questions that arise are: 0 Should an “average” or “effective” temperature and pressure be used to describe the conversion? 0 How are the terms “average” or “effective” defined? 0 How would this be done in terms of measurements commonly made in an operating plant?
Two ways of expressing conversion in terms of conveniently measured variables were discovered from the simulation. The first method was based on curve fitting the computer data. The second method has a theoretical basis. Data from many runs showed that a simple relationship existed at each
The information r e q u i r e d t o study t h e cracking reaction i s not necessarily t h e same as t h e informution n e e d e d f o r con trol. Control equations must contain v a r i a b l e s easy t o measure a n d m a n i p u l a t e in t h e p l a n t
conversion level among the temperature, pressure, flow rate, and tube diameter. At constant heat flux and heat flux distribution the conversion could be expressed as:
C
=
+
G I T c ~ P- CBNQ f ~4d. - cs
Conversion Correlation from Theoretical Basis. From an analysis of an isothermal and isobaric reactor, a functional group, called w , the severity factor, was developed which characterized the reactor performance :
(14)
The relationship covered a range of conversions from 40 to 5570, a pressure range of 225 p.s.i., a temperature range of 200” F., and 2 to 1 range of flow rates. The conversion as determined by this correlation was within 2% of the computer results. With this insight into the nature of the conversion dependence, a multiple regression analysis of the data obtained at the reactor outlet was then made. The conversion was correlated with the reactor effluent temperature and pressure, which can be conveniently measured. The multiple regression analysis also yielded a conversion expression of the same form as Equation 14. As expected, it gave approximately the same accuracy as the curve fitting method and was applicable in the same ranges. This analysis answers the questions raised above in regard to characterizing the reactor performance conveniently. The use of reactor outlet condition for this particular problem satisfied the requirements of accuracy in the correlation and ease of measurement in the plant. While it might be argued that a regression analysis should have been done first, it was not apparent from the maze of computer data how to describe this distributed system. The correlating of variables without any insight into what variables or any insight into the form of the equation is generally a time consuming, expensive, and often fruitless undertaking. The correlating of variables, without an appreciation of the process, can result in expressions that are physically and chemically false but which are statistically satisfactory.
The computer results were found to correlate as a semilog plot of severity referred to a base severity C = ai
+
u z loglo w / w b
(16)
The conversion calculated from Equation 16 is within 2% of the value obtained from the step by step computer solution. With an on-line process control computer either correlation gives in a few seconds approximately the same results as given by 11/2 minutes of IBM 704 calculations. Utilization of the Conversion Expressions. For steady-state conditions the effluent conversion is described by Equations 14 and 16 in terms of the total feed rate, the effective tube diameter, effluent pressure and temperature. The information required to study thoroughly the cracking reaction is not necessarily the information that is required to control the reactor. The control expressions have extracted from the reactor study the information required for control in terms of the variables that can be measured and manipulated for control. The forms of the conversion expressions are the important things and not necessarily the numerical values of the constants. In an operating plant, the relative accuracy of the conversion expression will vary with changes in the impurities in the feed stock, the heat flux distribution, the local coking, plugging, and the like. This is, of course, a reflection of the shifting profiles within the reactor. T o maintain either equation as an effective control device, it is necessary to reevaluate the constants in the equations periodically. This requires measuring VOL. 53, NO. 5
MAY 1961
347
the flow rate, the effluent composition, the outlet pressure and temperature. The effective rube diameter may be lumped in with constants as they are determined. One of the principal uses of the conversion expression for on-line control is, of course, to calculate conversion level as a function of current operating conditions. Where needed, compensating adjustments can be made by the computer. For example, the outlet pressure may vary with changes in the pressure level of the purification section following the reactor. The on-line computer can compensate for the pressure change b y a change in the outlet temperature or the flow rate according to the conversion expression to maintain the desired conversion. The conversion expression also serves as a predictor. For example, the plant operator wishes to run the thermal cracker over a period of time to make the required product rate at maximum profit or minimum cost. This requires knowing how to vary the process variables with time in order to shut down the reactor a t the optimum time due to coking. The conversion expression, a coking rate expression, heat and material balance equations, as well as economic factors are required. Methods already described (4, 5) will handle this predictive problem. Derivation of Kinetic Equations
Through a section of a tubular reactor the volumetric flow rate per unit time is
Using the definition of conversion and simplifying, Equation A-5 may be rewritten:
=
= = =
3Pf11n(l -pC2)2P2
c* -- - 3 -c 2P
P
Substituting the definitions of the CA, CB, and CC and re-arranging, the Equation A3 may be written as A-4 shown below. Equating Equations A-2 and A-4, and rearranging gives Equation A-5 below. The Equation 13 is the finite difference approximation of Equation A-5. For an isothermal and isobaric reactor, the Equation A-5 may be integrated to give conversion as a function of severity factor.
= = = = = = = =
(A-7)
w =
kA (2 ) I = severity factor z NoRT 2
Nomenclature
P
a], a2 = constants in Equation 16
A,
= cross sectional area of tube, sq.
AAO
=
CA
ft . incremental outside heat transfer area, sq. ft. = conversion, fraction of compound A that reacts = concentration of Compound A
CB
cu. ft. = concentration of Compound B =
moles/hr.
(i f XX ) RT’ moles/hr.
= concentration of Compound C
( -) X
moles/hr. No X RT’ cu. ft. = specific heat, B.t.u./lb.-” F. CP 6 1 , CZ, c3, c4, c j = constants in Equation 14 d = tube diameter, fr. AE* = energy of activation, B.t.u.;lb., mole-’ F. f = friction factor g = constant = 4.173 X 108, ft./ hr.2 G = mass rate of flow = W/A,, Ib./hr. sq. ft. 6 = inside film convection heat transfer coefficient, B.t.u./hr. sq. f t . - O F. h, = combined inside scale, pipe wall, and outside film heat transfer coefficient, B.t.u./hr. sa. ft.-” F. H, = heat of reaction, B.t.u./lb.mole k = thermal conductivity, B.t.u./ (hr.) (sq. ft.) (” F./ft.) =
+
E..
’ F.
product rate, moles/hr. upper limit on the absolute difference between the trial X and calculated X = compressiblity factor = increment = density, Ib./cu. ft. = latent heat of vaporization, B.t.u./lb. = severity factor, Equation A-9 or 15 =
st
A
N o -x
+
= tube wall temperature, = mass flow rate, Ib./hr.
(A-9)
C
stant, cu. ft./hr. reverse reaction velocity constant (cu. ft.2/mole) kiikz tube length, ft. equivalent length of straight pipe to 180° bend, ft. molecular weight feed rate to reactor, moles/hr. total flow rate in reactor = No X,moles/hr. pressure, p.s.i. vapor pressure, p.s.i. Prandtl number heat transferred, B.t.u./hr. gas constant Reynolds number temperature of flowing fluid, O
@ = l +kz- -P kl R T
cu. ft.
The rate of reaction is expressed as:
= =
where the integration is taken over the length of the reaction zone. For an isothermal and isobaric reaction, the Equation A-6 may be integrated analytically:
C,
The average residence time is
= Arrhenius constant = forward reaction velocity con-
x
w
=
Subscripts 1 = upstream 2 = downstream BV = boiling vapor MG = mixed gases max = maximum b = base C = critical R = reduced U = U-bend 2 = inside e = effective = average avg = two phase tP Acknowledgment
The authors acknowledge the assistance of David Roberts of Space Technology Laboratories, Inc., for his skillfull coding of the problem. literature Cited
(1) Brown, G. G. and others, “Unit Operations,” p. 140, John Wiley, New York, 1950. (2) Lockhart, R. W., Martinelli, R. C., Chem. Eng. Progr. 45, 39-48 (1948). (3) McAdams; W. A., “Heat Transmission,” 3rd ed., p. 396-8, McGrawHill, New York. 1954. (4) Roberts, S. M., Advances in Computational and Mathematical Techniques in Chemical Engineering, Chern. Eng. Progr. Symflosium Ser. 5 6 , No. 31 103-10 (1960). (5) Roberts, S. M., Proceedings of the “Optimization Techniques in Chemical Engineering Symposium,” p. 169-88, sponsored by New York University, AIChE, ORSA, New York University, New York, May 18, 1960.
RECEIVED for review August 12, 1959 RESUBMITTED August 18, 1960 ACCEPTEDFebruary 21, 1961
348
INDUSTRIAL A N D ENGINEERING CHEMISTRY
