Modeling, Simulation, and Operation of a Sabatier Reactor Peter J. Lunde The Center for the Environment and Man. Inc.. Hartford, Connecticut 06120
+
+
+
The Sabatier reaction for the methanation of carbon dioxide, C02 4H2 2H20 CHI heat, is an important reaction in closed-cycle space life-support systems in which metabolic COn is used to provide product water for subsequent electrolysis into breathing oxygen. The basic kinetics of the rutheniumcatalyzed gas-phase reaction, which have been reported in an earlier paper, have been applied to a dynamic model of the gas-cooled reactor which was solved on a large digital computer. Separate zones for catalyst, reacting gas, wall, coolant, and outside insulation were modeled. Each zone varied axially but not radially in temperature, permitting simulation of axial conductivity in wall and bed with wall temperature as a dependent variable. Good experimental correlations of reactor temperature distribution and conversion efficiency were achieved at the steady state and for a series of step changes in feed composition which included a 15-min "off" period and subsequent recovery.
Introduction When a space station is placed in long-term orbit the life support system will recover carbon dioxide from the cabin for reprocessing into breathing oxygen. This paper is concerned with a key element of this oxygen generation system, the chemical reactor which combines carbon dioxide that has been extracted from the cabin atmosphere with hydrogen to produce water vapor which will be electrolyzed to provide cabin oxygen. The desired reaction, called the Sabatier reaction after the Belgian chemist who first investigated the hydrogenation of carbon compounds, proceeds with a proper catalyst a t 400-700°F according to
+
+
+
4H2 CO, == 2H,O CHI heat (1 1 The forward equilibrium is favored a t low temperatures. A simplified flow diagram for the entire oxygen reclamation system to be used on board a Space Station is shown in Figure 1. The water balance shown is the key to appreciation of the stoichiometry of the chemical reactor. Extra hydrogen for the reactor must be supplied from storage, and is conveniently made available in the form of water by using ordinary frozen foods instead of the dehydrated fare in current spaceflight use. Methane can then be dumped to space instead of the stoichiometrically more desirable acetylene, which would present far more difficult processing problems. Since supply of hydrogen by water storage also provides by-product oxygen, reduction of all the available carbon dioxide is unnecessary. Thus the feed mixture to the Sabatier reactor can be quite rich in COz, with a choice of H2:C02 ratio from 2.5 up to the stoichiometric ratio of 4.0. The Sabatier reaction is also important in the manufacture of methane from synthesis gas, for which nickel is the customary catalyst. When high one-pass yields are desired the equilibrium thermodynamics call for the low temperatures a t which the Sabatier reaction is predominant. The first step in the development of a design procedure for the space life-support application of the Sabatier reaction was the development of the basic reaction kinetics for the very active ruthenium catalyst. This has been reported in an earlier paper (Lunde and Kester, 1974). The present paper describes the modeling, simulation, and operation of a Sabatier reactor. The work is unique in its continuity from kinetic data to successful laboratory confirmation of the mathematical simulation.
Previous Work There is considerable background literature regarding 226
Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 3, 1974
formulation and solution of the equations describing exothermic reactions in tubular reactors. Solution of all the equations which can be written for a realistic situation overwhelms even the largest computers; therefore investigators in the field always define their models carefully to consider those problems regarded as most important. The easiest and most common simplification is to eliminate time dependence and be content with the steadystate solution to the problem. This approach is satisfactory for the design of many industrial reactors but unsatisfying for process control applications and stability investigations. In the present study the steady-state solution was of little value because a time-dependent simulation was clearly needed to predict performance of an intermittently active reactor in orbit around the earth. Dassau and Wolfgang (1963) attacked the general exothermic tubular reactor problem with an analog computer, simplifying the reactor to a string of homogeneous noninteractive reaction nodes so that ordinary differential equations could provide a solution. In 1966 Marek and HlavkEek introduced the use of digital computation to the problem, permitting a more accurate description of the axial temperature profile by partial differential equations, but limited their model to the less general adiabatic reactor. Gilles (1968) added cooling through the reactor walls and noted the extreme sensitivity of the predicted temperature profile to such cooling. However, only a simple hypothetical reaction was considered because of computational restrictions. Wicke and Padberg (1968) solved a similar model and described the characteristic reaction zone of the nonadiabatic reactor, emphasizing the movement with time of this zone as heat and mass transfer relationships were changed. Marek and HlavhEek (1968) added the radial variation of temperature across the catalyst bed, but had to eliminate time dependence to solve their equations. McGreavy and Cresswell (1969) further discussed the computational problem and proposed an analytical simplification which unfortunately holds only for the adiabatic reactor. In 1970 Andersen and Coull discussed the inadequacies of the plug flow models most investigators use and detailed several two-dimensional steady-state models. Heat transfer through or along the reactor wall was not considered. Paynter, et al. (1970), again attacked the problem of simulating the radial dimension, introducing the concept of reducing the size of the axial distance increments in zones of high reaction rate in order to reduce computation time. Sinai and Foss (1970) presented a lumped-parameter reactor model with accompanying liquid phase experi-
10Hi0
ELECTROLYSIS 10HzO+lOH~+5
0,
502
FOOD 2 C i H z
IO Hi
-
SABATIER
fCoz
OUMP TO
10Hz + $ COzSHzO+)CH,
4co2 METABOLISM Z C z H z + 5 02-
-
4COz +2HzO
t
(An A p p r o x i m a t e C o m p o s i t i o n )
Figure 1. Simplified molar process flows in a space station oxygen reclamation system.
mental data. The model itself is similar to the present work although it does not consider axial heat transfer along the reactor wall. Vortmeyer and Jahnel (1972) discussed the dynamics of the characteristic reactor temperature profile in detail, stressing the importance of axial conduction within the reactor bed on the temperature profile. DeWasch and Froment (1971) illustrated with a computer study the importance of both axial and radial heat transfer on the reactor performance. Eigenberger (1972) detailed results from a similar dynamic model, emphasizing the importance of axial heat conduction within the catalyst bed. The reaction kinetics were as simple as possible, and while heat transfer to the wall was considered, the wall temperature was invariant. His solutions were obtained by use of an automatic length increment size adjustment routine. None of the references cited consider the importance of axial conduction along the reactor wall to an accurate simulation, although axial conduction within the bed is often discussed. The wall conduction is especially important when the wall temperature is to be a dependent variable because the heat transmission of even a thin section of metal is often higher than a bed of the ceramic catalyst supports now widely used because their high surface area allows a high specific reaction rate. While the effects of conduction in the wall and bed are similar, the reactor designer must consider both and be prepared to modify bed and wall conductivity to meet his needs. Mathematical Model It is convenient to divide the mathematical model into a thermal model and a chemical model. These are linked by the heat generation rate, which is determined in the chemical model according to the temperatures of the thermal model. The thermal model chosen for this work does not consider radial temperature gradients in a continuous fashion, but instead assumes a constant (or linearly variable) temperature in each of five zones-the catalyst, reactant gas, wall, coolant gas, and insulation. These temperatures are continuously variable in length and interact realistically, thus modeling axial conduction in bed and wall and permitting the wall cooling rate to be a dependent variable. These are important assumptions which consume computer time, but are necessary to an accurate simulation of a practical reactor. The chemical model utilizes ordinary material balances and a kinetic expression derived in previous work from this laboratory. The key assumption necessary to solve the
chemical equations is that the exponential rise in reaction rate is limited at high temperatures by otherwise unmodeled phenomena. Otherwise the predicted reaction rate becomes so fast that simulation is extremely difficult. The Thermal Model. Figure 2 is a drawing of the reactor which was chosen for modeling. It is a single-pass tubular packed-bed reactor cooled by an annular cocurrent flow of coolant fluid, with the entire shell covered by insulation. Inlet gases pass through an inlet section and react in the packed reactor section. The model assumes that heat is generated on the catalyst surface and transfers by convection to the bulk of the reacting gas, which moves in plug flow. This gas heats the reactor wall by convection, and heat is transferred by conduction across the reactor wall and then by convection to the coolant stream. The coolant stream loses some of its heat to the ambient through the insulation. The catalyst reactant gas and coolant temperatures were considered to be functions of time and axial distance while the reactor wall, coolant wall, and insulation temperatures were considered radially variable as well, according to Fourier’s law. This allowed a simpler simulation by the use of experimentally derived overall heat transfer coefficients. Separate radial and axial thermal conductivities were used in the reactor bed and reactor wall, permitting lumping of the radial convection to the reactor wall with the radial conduction through i t and, if desired, a thermally anisotropic reactor bed. Since axial heat transfer was considered the key to a successful simulation, axial temperature gradients were continuously variable with length in both the reactor bed and the reactor wall as well as in the gas streams. The following differential equations express heat transfer among the five thermal zones. 1.Reactor (catalyst) heat balance
2. Reactant gas heat balance
3. Reactor wall heat balance
4. Coolant heat balance
where
c‘, =
_ I ‘0
h, r,
+
1 1 -r, k,
ro ln -
r,
T r , Tp, Qr,T,, and T, are functions of reactor length and time. Reactant gas composition and temperature determine cpg. Boundary conditions for these equations will always be satisfied by the initial temperature distribution for a particular simulation. Since the heat storage in the gaseous streams was small, m, and mp were set equal to zero prior to solution of the equations. One additional parameter is needed for solution of these equations-the Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 3, 1974
227
P1 anned Cool a n t Gas Entrance
t
Actual Coolant Gas Entrance
Coolant Gas Exit
I
n
-
Product Gas
Feed + Gas
Scale in Inches r ’
I
1
1
0
I
I
2
4
6
Figure 2. The laboratory reactor. For detailed dimensions, see Table II. heat generation rate Qr, a function of reactor length and time determined by the temperature, reaction kinetics, and material balance of the chemical model. Chemical Model. A given reactor location has an incoming gas composition, a rate of reaction, a corresponding heat of reaction, and an outlet gas composition. The incoming gas composition and the local temperature of the reactor location determine the other three parameters. The rate of reaction was calculated by using the rate equation for the Sabatier reaction given by Lunde and Kester (1974) for the yz% ruthenium catalyst later used in laboratory scale tests
where
K,(T)
=
exp[(1.0/1.987)(56,000/Tk2 f 34633/Tk 16.4 In Tk 0.00557Tk) + 33.1651 ( 7 )
+
E,, k , and n were determined experimentally. The partial pressures are functions of distance and time and are expressed in atmospheres, with time in hours. Simple stoichiometric relationships establish the value of P H ~ PcH,, , and P H ~as Oa function of PcoZ according to the basic Sabatier reaction, with the feed composition an initial condition a t the reactor entrance. The heat of reaction is proportional to the reaction rate
where
(9) a relationship derived from the specific heat of the gases as a function of temperature and the 298°K heat of reaction of -39,433 cal/mol.
Simulation of the Reactor Using a Digital Computer The simultaneous solution of the heat, material, and reaction equations was carried out by reducing them to difference form and selecting appropriate Ax and A t intervals with which to proceed down the reactor. Since digital computation inherently requires this lumping of a continuous process into discrete steps, oversize increments can produce inaccuracies and a point of instability can be reached in which the approximation 228
Ind. Eng. Chem., Process Des. Develop., Vol. 13,No. 3,1974
produces impossible conditions; e.g., sufficient heat may be conducted out of an element to lower its temperature below that of the element receiving the heat if the time increment is too large. Such conditions were noted by inspection of the difference equations and used as criteria for selection of time and distance increment sizes. The heat generation rate is not at all constant along the reactor length. The reaction tends to proceed very rapidly once it has been preheated to an ignition temperature which is determined by the gas composition and pressure and the catalyst activity. It then becomes “autothermal,” proceeding rapidly to a high-temperature equilibrium (as high as 1500°F for the Sabatier reaction) since the heat generated increases the rate of reaction exponentially. Further reaction beyond the autothermal zone occurs only when the gas mixture is cooled while still in the presence of the catalyst. Final temperatures as low as 400-500”F are necessary to achieve the nearly complete conversion desired for the space life-support application of the reaction. As the temperature drops, the slower reaction rate makes a proper simulation more important than in the autothermal reaction zone, for most of the catalyst volume is devoted to this zone of controlled cooling. The autothermal operation cannot generally be controlled and is difficult to simulate because the required time and distance increments become so small that the solution overloads even modern computer systems. The peak temperature builds up even more rapidly in the simulation than in the real reactor because the incremental solution lets the heat accumulate during one time increment, conducting it out of an element a t a rate reflecting the consequent temperature rise only in the succeeding time increment. To reduce the computer time for simulation of the autothermal zone three techniques were used. (1)The distance increment size was made variable and was adjusted during each calculational pass through the reactor. (2) The time interval was made small by executing several time increments sequentially within a given distance increment when time increments of less than 1 sec were required. (3) An upper bound was placed on the reaction rate by limiting i t to the rate a t a certain maximum temperature (830°F was used). This preserved the kinetics in the important lower temperature zones of the reactor and reduced computation times in t h e high-temperature zone. There is probably no effect on the accuracy of the simulation, because the autothermal zone is so short that a large factor of error in its size would not be significant. With these techniques, the simulation used about 1 min of UNIVAC 1108 computer time for each 4.5 min of simulated time.
r '
I
--
i n l e t Section
--
. .-
~
- -_-
.- .
-- , -
-
Reactor Section
1
Outlet
~
Experimental 1100
-
~
--
c
Cooled S e c t i o n
c
1000 L rSteady S t a t e Experirrental Datal "-
0
Temperatures
A
Conversions I
"""1pULCr
-
700
Simulated Temperature Distributior
-
i
'\\\\
> ,\
Y
E
a
2
600 500
-
200
-
1
LJ A
0
2
L
4
U
6
-
8
10
-
12
I
14
I
16
1
18
1
1
20
1
1
22
1
'
24
1
-
L
0
----.__1
,I
26
D l s t a n c e From F r o n t o f R e a c t o r ( I n c h e s )
Figure 3. Comparison of simulated and experimental steady-state temperature and composition distributions.
Experimental Operation Laboratory Reactor. To evaluate the mathematical model, a laboratory reactor was built. After initial tests were run to establish heat transfer coefficients, actual reactor temperatures and composition distributions were determined a t steady state and as a function of time during startup and during severe fluctuations in the feed supPly. Figure 2 is a drawing of the laboratory reactor. It is a l%-in. diameter tube 25 in. long, of which 16.5 in. was packed with catalyst. Four hundred forty grams of catalyst was used, a scaleup by a factor of 100 over the original kinetics data (Lunde and Kester, 1974). The tube was instrumented with eight center line thermocouples and six gas sampling ports. Front half cooling was passive, consisting of heat loss through the 1 in. thick fiberglass insulation plus losses by axial conduction to the cooler areas. Cooling of the rear half of the reactor was by cocurrent nitrogen flow. Reacting gases were preheated externally in order to start the reaction. At lower flow rates, the preheater was unnecessary once the reactor had ignited due to axial conduction along the catalyst bed and wall. Thermal Parameter Tests. Heat transfer coefficients were determined empirically by comparing simulated data with experimental data from the nonreacting reactor a t steady state when fed with a preheated CO:! feed gas under three different conditions of flow and cooling. Table I summarizes values found for the coefficients, which were within the range of analytical estimates. Slightly different values were used for the final simulation. These gave better results in matching the active temperature distributions which were several hundred degrees higher. Steady State Tests. The test reactor was operated with a variety of feed flow rates (from 4.1 to 24.5 ft3/hr at 1 atm) and inlet molar flow ratios (Hz:CO:! from 2.1 to 4) a t two pressures (0.3 and 1.0 atm) with varying cooling flows (0 to 1.2 ft3/min a t 1 atm). Data from only one of these tests are reported here. (Further data are available from the author upon request.) In the course of these tests it
was found that with cocurrent cooling, which can be modeled by the computer program, the reaction was quenched unless the cooling flow was introduced downstream of the high-temperature region of the reactor. It was also found that above a certain feed rate the high temperature reaction zone would slowly proceed downstream and "blow" out of the reactor. Just below this critical feed rate the reaction was stable, and a t this condition steady-state and dynamic test data were accumulated to test the mathematical model. Steady-state data were well-matched by the simulation (Figure 3) after the heat transfer coefficients which had been approximated in thermal parameter testing were trimmed to the "Final Simulation" values listed in Table I. Experimental and calculated reactor conversions as a function of distance through the reactor are also shown in Figure 3. The experimental outlet conversion agrees quite well with the model prediction. Intermediate conversions are also in good agreement, considering that there was up to 6% CO in these samples and both the experimental conversion calculation and the mathematical model ignored the possibility of CO formation. Table I1 lists all the parameters used for the steadystate simulation. With the exception of manipulated variables, these were kept constant during the subsequent simulation of the dynamic reactor operation. Dynamic Tests. With a suitable experimental steadystate temperature and composition profile established, a definitive test of the mathematical model was made by comparison of the simulation with actual performance during severe upsets to the reactor feed. The upsets lasted about 15 min each and were followed by recovery periods of about 15 min. First the CO:! and then the H:! feed was reduced to one-half the steady-state value, then the feed hydrogen was stopped entirely. This final upset was followed by a long recovery period which was monitored for nearly 1 hr. The top graph of Figure 4 shows the upsets which were imposed on the reactor in the form of a plot of time us. COz and H:! feed flow. Crucial times just before a perturInd. Eng. Chem., Process Des. Develop., Vol. 13,No. 3,1974
229
Table 11. Program Input Parameters D
_
l
A
C
B
D
E
F
H
G
Y
Total length of reactor inlet section Total length of reactor catalyst bed section Total length of reactor outlet section Radius of catalyst bed Outer radius of reactor wall (bed radius reactor wall thickness a ) Extreme outer radius of assembly, including insulation Distance from front of reactor inlet section a t which cooling begins Distance from front of reactor inlet section a t which cooling ends Catalyst coefficient, nc Activation energy Rate constant Maximum reaction rate equal t o that at temperature Axial and radial wall conductivity (stainless steel) Wall specific heat Wall density Coolant specific heat Coolant molecular weight Catalyst particle diameter Catalyst thermal conductivity Catalyst specific heat Catalyst density Coolant flow rate Coolant inlet temperature Reactor pressure Ambient temperature Heat transfer coefficients, Btu/hr ft2 O F Catalyst-to-gas (h,) Gas-to-wall (h,,.) Wall-to-coolant gas (h,) Coolant gas through insulation to ambient (Ui) Feed gas Flow (base case) Temperature (at reactor inlet) Inlet composition, base case (mole fractions)
+
o EXPERIIVENT-
,
SIMULATION
( P o i n t s have been moved 2 m i n u t e s t o l e f t t o acCoLint for a l o n g sample
-5
4.
line)
A 3 20 40 60 80
109
120
,I
Test Time (Minutes)
Figure 4. Experimental feed perturbations (top) and the resultant experimental and simulated reactor outlet conversions (bottom). Letters A-H key time to temperature distributions in Figures 5-8, and do not necessarily correspond with experimental analysis points. T a b l e I. Overall Heat Transfer Coefficients (Btu/hr f t 2 O F ) Heat transfer test no. 1
2
3
Catalyst to gas, h, 10 10 10 Gas to outer coolant wall, h, 10 7 6.5 Coolant wall to coolant, h, 30 16 No flow Coolant to ambient, U1 1 . 5 1.5 0.9 (through insulation) Q
Final simulation 25 14
29 1.25b
a The value for catalyst-to-gas coefficient does not affect the fit when there is no reaction, so a value of 10 was arbitrarily chosen. * See footnote b, Table 11.
bation and just prior to recovery have been given letter designations so that experimental temperature distributions can be conveniently identified and compared with corresponding computer-predicted data. Temperature data are easily obtained and tell more about reactor performance than gas composition information. Figures 5-8 show the predicted and experimental reactor temperature distributions a t designated times A to H which are keyed to Figure 4. Figures 5-7 show the temperature distribution just prior to an upset and just prior to recovery, while Figure 8 shows two distributions during the recovery from the upset produced by the hydrogen flow stoppage. The bottom graph in Figure 4 shows the predicted value of the outlet composition as per cent conversion of the stoichiometrically lean constituent in the feed (usually H2) and the experimental conversions as calculated from chromatographic analyses taken every 11min. Predicted conversion accuracy is excellent except for samples a t 10 and 91 min. The -1O-min sample, experimentally recorded as 100% conversion, was probably less because the chromatograph gain was held fixed throughout the tests and not adjusted to record properly the small (0.05 in.) COz peak which would have been appropriate for 98% conversion. The 91-min sample indicates that actual recovery after the 15-min hydrogen flow stoppage was faster than that predicted analytically. This effect will be discussed later and is believed related to the operation of reactor (and simulation) very near the maximum possible reactor flow. The general correspondence of experimental and simulated reactor temperature distributions as shown in Fig230
Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 3, 1974
6.5 in.
16.5 in. 2.5 in. 0.72 in. 0.77a in. 1.5bin. 16.125 in. 25.5 in. 0.25 -31,000 Btu/lb-mol of COa 0.2339 X 1 O l o hr -1 a t m 4.25 830°F 11 Btu hr f t
O F
0.15 Btu,lb O F 488 lb/ft3 0.25 Btu,/lb O F 29 lb/lb-mol 0.12 in.d 0.001 Btu/hr ft O F 0.278 Btu/lb O F 68 lb/ft3 1.24 ft3/min 90 O F 14.7 Ib/in.2 (absolute) 80°F 25 14 29 1.256 0.638 ft3,‘min 362°F \COz 0.3 ] H ? 0.7
Actual wall, thickness was 0.035 in. This was multiplied by 1.5 to approximate the dynamics of the unmodeled but equally thick outer coolant wall. The true extreme outer radius was 2.0 in. Had this value been used a trial value of U1 = 0.94 would have produced the same simulations. These values were derived in an early correlation of the data presented by Lunde and Kester (1974). The catalyst was actually 0.125 X 0.125-in. cylinders. ures 5-8 was excellent. Specific discussions of each upset are given below. “Low COz” Perturbation (Times A-B). The match a t time A is very good, At time B the model is somewhat high in temperature but the correspondence is reasonable. “Low Hz” Perturbation (Times C-D). This test is more severe than the “Low COz” perturbation since Hz is the limiting feed component. Correspondence between model and experimental data is good a t beginning (C) but somewhat off a t the end (D). There is no explanation a t present for the unusual temperature inflection in the computer simulation (time D, 8-16 in.). “Off Hz” Perturbation (Times E-F). Simulation of the recovery from the previous perturbation (time E) is not up to previous standards, being 60 to 140°F low in the 18-23 in. cooling zone. Careful examination of the com-
1200 r
-
Inlet Section
-- --
2
6
Outlet
I
1
1
1000 ‘loo
0
4
8
10
12
14
16
18
20
22
24
O
26
2
4
6
D i s t a n c e From F r o n t o f R e a c t o r ( I n c h e s )
Figure 5. Center temperature distributions for the “low COz” perturbation. Original near-steady state (time = 0): (A) experimental, (A’) computer simulation. After 15 rnin of low COz (time = 15 min); (B) experimental, (B’) computer simulation. I
L ~
1200
r
Inlet Section
Outlet Section
i -
8
10
12
14
16
18
20
22
24
26
D i s t a n c e From Front o f R e a c t o r ( I n c h e s )
-1
Figure 7. Center temperature distributions for the “off Hz” perturbation. After recovery from “low Hz” (time = 64 min): (E) experimental, (E’) computer simulation, After 15 min with Hz off plus 1 min of recovery (time = 80 min): (F)experimental, (F’) computer simulation.
~
II-1200 \
I
Cooled S e c t i o n
1100 1000
E
Outlet
- I1 --
Reactor Section
1
1
900 p
Y
1
Inlet Section
800
X
-
l
I 700
a
.!
600
1 I
C1,D’
C
2
4
6
8
10
l
I
1
.
1
I
l
1
12
14
16
18
20
22
24
26
1
400
f
300
r
---.
~~
0
D i s t a n c e From F r o n t o f R e a c t o r ( I n c h e s )
2
4
E
8
10
12
14
16
18
20
22
24
26
D i s t a n c e Frow F r o n t o f R e a c t o r (:nciies)
Figure 6. Center temperature distributions for the “low Hz” perturbation. After recovery from “low COz” (time = 34 m i d : (C) experimental, (C’) computer simulation. After 15 min of low Hz (time = 48 min): (D) experimental,(D’) computer simulation. puter output between times D and E shows that after time D temperatures in the cooling section continued to drop for 9 min and that a recovery was still in progress a t time E. The point at 20.25 in., for instance, which was 534°F at 48 min computer time (time D), fell as low as 480°F (at 57 min) and had risen to only 502°F at 64 min (time E). Thus the thermal time constants of the model seem a little longer than actually existed. The match a t time F after 15 rnin with the hydrogen feed completely off is somewhat better and although the same effects are evident the model is perfectly adequate for most purposes. Recovery from “Off Hz” Perturbation (Times G and H). Examination of both computer and experimental data at times G and H indicates that the reactor will soon “blow out;” i.e., the reaction peak temperature is moving down the reactor and the reaction will probably eventually cease entirely. This effect has been noted in the laboratory in unrecorded runs at too high a flow rate and is now shown to be a hazard of feed interruptions. The computer fit is not as good as in less sensitive operating situations, but the presence of an operating problem
Figure 8. Center temperature distribution during recovery from the “off Hz” perturbation. After 21 rnin (time = 100 rnin): (G) experimental, (G’) computer simulation. After 51 rnin (time = 130 min): (H) experimental, (H’)computer simulation, is clear and the problem is somewhat worse in the computer model than in practice. This conclusion is confirmed by examining the outlet conversions beginning at time F (plotted in Figure 4) noting the peaking, decline, and final rise of the computed conversion to the experimental value. (A similar effect is noted in comparing the experimental temperature distribution at time G us. time
H.1 Reactor Control. In one important way, neither the simulation nor the laboratory reactor accurately represented the performance of a working Sabatier reactor. The missing element is reactor control, which was not modeled because it would lengthen already considerable computation time requirements. An ideal reactor control might sense outlet conversion and then cause appropriate adjustments in coolant flow to give the desired chemical conversion. The practical difficulties of such a direct feedback scheme, however, are almost insurmountable because of the time delay between coolant change and measurement of the outlet composition change. Ind. Eng. Chem., Process Des. Develop., Vol. 13, No. 3, 1974
231
A more practical control system would adjust the reactor outlet temperature to a specified value desirable to achieve good conversion (500°F would be appropriate). Automatic or manual control of this setpoint would then maximize conversion for a given reactor design and feed conditions. Such a scheme was tested manually in the laboratory during the reactor tests. It was found to be responsive and generally desirable when the coolant flow was countercurrent to the reactor flow, and sluggish although adequate when the cocurrent coolant flow represented by the model was used. Accuracy of the Simulation. If the anticipated automatic control of the reactor is hypothesized, the importance of correctly determined heat transfer coefficients diminishes because the control will adjust coolant flow to give a correctly simulated reactor outlet temperature. This is an important point, for temperatures in the mathematical model (as well as the experimental reactor) were quite sensitive to the cooling gas flow rate and the heat loss direct to ambient from the hot reaction zone at the front of the reactor, which is difficult to predict analytically. With this single exception of heat transfer coefficients, the mathematical model accurately predicted reactor performance exclusively from fundamental relationships which incorporated a number of simplifying assumptions. It is especially gratifying that good accuracy was achieved without modeling the radial temperature gradients within the reactor, for this complication could well have driven the simulation off the largest digital computer with hybrid simulation the only recourse.
Conclusions 1. The reaction kinetics experimentally determined by Lunde and Kester (1974) for the ruthenium-catalyzed Sabatier reaction have been applied to a (dynamic) mathe1H, CO, 2H,O + CH, + heat matical reactor model and the results confirmed by laboratory testing of an actual reactor. 2. The mathematical reactor model assumed axial and radial heat transfer with continuous axial temperature distributions, but lumped radial temperatures into one zone each for the catalyst, gas, wall, cocurrent annular coolant flow, and outer insulation. These assumptions permitted accurate simulation of actual reactor performance. 3. Simulation difficulties were encountered in the initial reaction or autothermal zone where the rate of reaction was high enough to rapidly raise the reaction temperature and hence the reaction rate. Three techniques were used to reduce the computer time necessary to simulate the autothermal zone. (a) Length increment size was varied with distance after each calculational pass through the reactor so that short increments were used only in zones of rapid reaction. (b) When required time increments were small, many successive time steps were executed a t a single distance increment before moving to the next distance increment. (c) The reaction rate as determined from the kinetics of Lunde and Kester (1974) was limited to that calculated at 830°F to reduce the need for very small time and distance increments. These techniques (especially the third) reduced UNTVAC 1108 computer time to 1 min for every 4.5 min of dynamic simulation time. 4. Accurate steady-state simulation was quite sensitive to the values used for heat transfer coefficients. After minor adjustment of these coefficients, the dynamic simulation accurately represented reactor outlet compositions and temperature distributions before, during, and after
+
232
*
Ind. Eng. Chern., Process Des. Develop., Vol. 13,No. 3,1974
severe feed disturbances, including a 15-min stoppage of the hydrogen feed flow. 5. Control of the Sabatier reactor is most conveniently achieved through regulation of the coolant flow rate. Response is best when a countercurrent flow of coolant gas is used but satisfactory control is achieved with the cocurrent coolant flow which was simulated. When coolant flow control is to be used, the need for accurate prediction of heat transfer coefficients diminishes because changes in the coolant flow will hold the desired reactor outlet temperature. 6. When the reactor suffered severe feed disturbances after reaching a steady state near the maximum feed flow rates, conditions were predicted by the model and experimentally confirmed which caused the relatively short zone of high reaction to move down the reactor, leading to poor conversion and eventually to a shutdown of the reactor. Therefore, it cannot be assumed that a reactor, once started, will necessarily assume a particular steady state condition under all circumstances.
Acknowledgments The work reported in this paper was performed at the Hamilton Standard Division of United Aircraft Corporation under National Aeronautics and Space Administration Contract 9-9844. Thanks are due to Charles Verostko a t the Manned Spacecraft Center for his guidance throughout the program, and to R. A. Baum, who contributed to the mathematical modeling and did most of the computer programming. Nomenclature Ab, = catalyst bed superficial peripheral surface area (27rr,),ft2/ft of bed Abx = catalyst bed superficial cross-sectional area (nr12), ft2
A, = outside area of insulation, ft2/ft of bed A,, = wall outside peripheral area (27rr,,), ft2/ft of bed A,, = wall cross-sectional area [n(ro2- r,2)],ft2 ( A / & = surface-to-volume ratio for the catalyst pellets, ft- 1 E, = activation energy, btu/lb-mol of COz AH, = heat of reaction, btu/lb-mol of COz Ke = equilibrium constant, atm-2 [Pt] = partial pressure of subscripted species i, atm Qr = rate of heat generation, btu/hr f t of bed R = gas constant, 0.730 atm ft3/lb-mol "R 5" = temperature, O R Tk = temperature, "K U = overall heat transfer coefficient, btu/hr ft2O R W = gas flowrate, lb/hr cp = specific heat, btu/lb O R h = local heat transfer coefficient, btu/hr ft2 OR k = thermal conductivity, btu/hr ft O R k = reaction rate constant, hr-1 atm-o.25 m = weight per unit bed length, lb/ft of bed n = catalyst coefficient, dimensionless r = radius of cylindrical wall surface, f t t = time, hr x = axial distance or length along bed, f t
Subscripts a = ambient b = bed c = coolant gas g = reactantgas i = inside 1 = losstoambient o = outside p = peripheral r = reaction w = wall x = cross-sectional
Literature Cited Anderson,T. S., Coull, J.,A.l.Ch.E. J.. 16, 542 (1970). Dassau, W.J., Wolfgang, G. H.. Chem. Eng. Progr., 59,4, 43 (1963). DeWasch, A. P., Froment, G. P., Chem. Eng. Sci.. 26, 629 (1971). Eigenberger, G.. Chem. Eng. Sci., 27, 1909 (1972). Gilles. E. D.. Chem.-lno.-Techn.,40. 469 (1968). Lunde, P. J.', Kester, F. L.. lnd. Eng. Chem.. Process Des. Develop., 13, 27 (1974). McGreavy, C.. Cresswell, D. L.. Can. J. Chem. Eng., 4 7 , 6 , 583 (1969).
Marek, M., HlavaEek, V., Chem. Eng. Sci., 21,(193 (1966). Marek, M., HlavaEek, V., Chem.-lng.-Techn., 40,1086 (1968). Paynter, J. D., Dranoff, J. S.,Bankoff, S. G., Ind. Eng. Chem., Process Des. Develop., 9,303 (1970). Sinai, J., Foss, A.S..A.l.Ch.E. J., 16, 658 (1970). Chem. Eng. Sci., 27, 1485 (1972). Vortmeyer, D., Jahnel, W.. Wicke. E., Padberg, G., Chem.-Ing.-Techn., 40,1033 (1968).
Received for review M a y 9, 1973 Accepted M a r c h 6,1974
Nonisothermal Kinetics Studies of the Hydrodesulfurization of Coal Alfred L. Yergey,* Frederick W. Lampe,' Marvin L. Vestal,2 Alan G. Day, Gordon J. Fergusson, William H. Johnston,3 Judith S. Snyderman, Robert H. E ~ s e n h i g hand , ~ James E. Hudson Scientific Research lnstruments Corporation, Baltimore. Maryland 27207
The application of nonisothermal kinetic methods to the hydrodesulfurization of coal is described. The results indicate that with only a few exceptions, the hydrogen-induced release of sulfur as hydrogen sulfide can be described adequately, in a practical kinetic sense, by five processes. We have designated the five processes in terms of the source of sulfur as Organic I , Organic II, Pyrite, Sulfide, and Organic I I I. Preexponential factors and activation energies are determined for each. The importance of the back-reaction of H2S with partially desulfurized coal is discussed and it is shown that the most important of the back-reactions in the inhibition of H2S release is reaction with iron. It is also shown that Organic I I I is most likely formed by back-reaction of H2S with carbon in the partially desulfurized coal.
Introduction The current recognition of sulfur oxides as a major air pollutant has generated a renewed interest in an old industrial problem, namely, the desulfurization of coal and coke prior to their combustion. It is generally accepted that sulfur occurs in coal in three forms: (1) in organic chemical combinations and, therefore, as an integral part of the coal substance itself; ( 2 ) as pyrites and/or marcasites; and (3) as sulfates, which generally are of calcium and iron. In view of the nature and abundance of the organic sulfur in coal, it is readily apparent that chemical methods of sulfur removal, in which the coal substance itself is chemically altered, present the only practical means of efficient desulfurization prior to combustion. Thus knowledge of the chemical kinetics of this very complex chemical transformation is desirable from the practical-process point of view and because of any fundamental information concerning structure that may be derived. The desulfurization of coal and coke has been well-studied over the past half-century beginning with the investigations of Wibaut and Stoffel (1919) in Europe and Powell (1920) in the United States. An excellent review of work done on this subject prior to 1932 is given by Snow (1932), while a summary of work done to about 1945 is presented by Thiessen (1945). Reports of desulfurization studies have continued to appear to the present time (Batchelor, et al., 1960; Blayden, 1958; Brewer, 1949; Chapman, 1955; Chowlhury, et al., 1952; Curran, et al., 1958; Fuchs, 1951; Gray, et al., 1970; Mahmoud, et al., 1969; Mason, 1959; Zielke, et al., 1954, 1970). The general picture of hydrodesulfurization emerging
' *
Department of Chemistry, Pennsylvania State University, University Park, Pa. 16802. Department of Chemistry, University of Utah, Salt Lake City, Utah 84112. Deceased. Department of Fuel Sciences, Pennsylvania State University, University Park, Pa. 16902.
from these studies is that the organically. bound sulfur and the sulfur contained in pyrites are evolved from the coal or coke principally as hydrogen sulfide at processes that begin to be of significance at about 300°C. Presumably ferrous sulfide is a by-product of the pyrite decompositions and this sulfide is thought to decompose in the range of 400-600°C. Sulfates are thought to be reduced to sulfides a t temperatures below 500°C. Apparently a significant fraction of the original organically bound sulfur is retained in the solid product as chemically bound sulfur that is at least partially removable by hydrogen treatment a t temperatures of 700-1000°C. The reaction of the hydrogen sulfide with the partially desulfurized coal is rapid and unless this back-reaction is somehow inhibited practical desulfurizations are precluded. This back-reaction of the hydrogen sulfide may also be an important source of the strongly bound sulfur in the solid partially desulfurized product. A common feature of all the early investigations of coal desulfurization is the use of the classical, time-honored method of chemical kinetics in which the extent of conversion is measured as a function of time for a series of experiments, each series at a different but constant temperature. Unfortunately, this method is severely limited when applied to such a complex conversion as the hydrodesulfurization of coal because of the occurrence of a number of parallel chemical reactions, both at the desired reaction temperature and during the time the sample is being brought to the reaction temperature. The result is that one observes, at each temperature measurement, the overall end result of a superposition of reactions. A method in which the velocity is measured at a constant and known heating rate of the solid sample, and which circumvents this difficulty, has been developed recently by Juntgen and coworkers (Hanbaba, et al., 1968; Juntgen, 1964; Juntgen and Traenckner, 1964; Juntgen and van Heek, 1968; Peters and Juntgen, 1965; van Heek, et al., 1967a,b). The method permits evaluation of the usual kiInd. Eng. Chern., Process Des. Develop., Vol. 13,No. 3, 1974
233
