Ind. Eng. Chem. Res. 1992, 31,5849
58
Simulation of Transport in Laminar, Tubular Reactors and Application to Ethane Pyrolysis John C.Fagley, Jr. BP Research, 4440 Warrensville Center Rd, Cleveland, Ohio 44128
Simulation of laminar flow reactors has been an area of interest for over three decades. Previous workers in this area have employed a wide variety of simplifying assumptions, e.g., irreversible first-order reaction, constant physical properties, and assumed parabolic velocity profile. One motivation for making these simplifying assumptions has been to make solution of the governing equation sets tractable. With the advent of more powerful machines and more robust numerical methods, these assumptions are no longer requisite to posing a solvable problem. Rigorous solution of the full set of governing equations will now allow us to accurately simulate (and interpret experiments and scale up) a wider variety of conditions. This is demonstrated through application of the Phoenics simulator to laminar ethane pyrolysis. The ability to predict natural convection in the system with little additional effort is also demonstrated.
Introduction Several articles have been written on laminar flow reactor simulation a t different levels of complexity. Early work by Cleland and Wilhelm (1956) and Hsu (1965) dealt with the accuracy of the plug flow assumption by solution of the governing two-dimensional concentration equation for an isothermal system with first-order, irreversible reaction. The results of these works included plots of radial concentration profiles and a criterion for determining when the plug flow assumption could be used. (Note that this criterion applies only to plug flow in an isothermal system.) Rothenberg and Smith (1966) were among the first to couple an energy balance along with the species conservation equation, and results showed the impact of heat of reaction on temperature profiles and heat-transfer coefficients. Assumptions incorporated into their model included fmt-order, irreversible reaction, parabolic velocity profile, isothermal wall, and constant physical properties. Trombetta and Happel (1965) and Andersen and Coull (1970) included descriptions for isothermal or adiabatic wall, nonequimolar reactions, and, in the latter work, temperature dependence of molecular diffusivity. In both cases, temperature variation of specific heat and thermal conductivity were neglected, and parabolic velocity profiles were assumed. Results were presented as concentration and temperature profiles, heat-transfer-coefficient correction factors, and factors for correction of reactor length for laminar as compared to plug flow. Sundaram and Froment (1979) presented an equation for Nusselt number correction factors in laminar flow using similar assumptions: parabolic velocity profile, constant physical properties, isothermal wall, and fmt-order irreversible reaction. These previous works demonstrated the ability to solve the governing mass and energy balance equations for simplified systems, i.e., systems with simplified kinetics, only two or three species of importance, etc. Also, the published results are useful for determining when the plug flow assumption is warranted and give ways of correcting for laminar flow in reaction systems where the simplifying assumptions apply. Unfortunately, the assumptions used in these earlier works do not adequately apply to many specific reaction systems. The assumption of constant wall temperature or an adiabatic wall is frequently not useful in simulation of homogeneous, laminar flow systems, because the reaction of interest is a thermal decomposition being carried out in a tubular reactor with a heating length, reaction zone, and cooling length. In turn, the radial and longitudinal variation in temperature and composition realized in such a system may have significant impact on 0888-5885f 92f 2631-QQ58$Q3.QQ fQ
the physical properties of interest: molecular diffusivity, specific heat, thermal conductivity, and viscosity. Even when the tube temperature in the reaction zone itself is isothermal, results of the earlier works mentioned above may not apply due to (1) radial temperatures gradients existing in the gas at the beginning of the reaction zone as a result of heat transfer in the heat-up zone and (2) significant conversion of reactants near the wall in the heat-up zone, as well as near the centerline in the cooling zone. In cases with severe radial temperature gradients, even the assumption of parabolic velocity profile may break down due to rapid heating (and lower density) of gas near the wall, as well as temperature dependence of viscosity and nonequimolar reaction. In terms of describing the kinetics of the system, the assumption of a single irreversible, finborder reaction may be too approximate. For instance, we may want to use the simulation to help determine or verify the details of a kinetic mechanism. When serial/parallel reaction mechanisms are involved, the local composition of the products formed may become important in physical property evaluation, and the heat of reaction of individual steps may be necessary for accurate description of the temperature profile. In this article, the ability to solve more detailed descriptions of the governing momentum, mass, and energy balance equations for laminar flow of gas within a tubular reador is demonstrated. The model allows for temperature and composition variation of density, viscosity, thermal conductivity, specific heat, and molecular diffusivity. Multiple reaction steps, both serial and parallel, can be accommodated. Each reaction step is assigned its own heat of reaction which may vary with temperature. Molecular diffusion of each species is allowed, taking into account that some species in the system may diffuse significantly faster than others. Wall temperature is allowed to vary with length in the reactor, giving us the ability to simulate heat-up, reaction, and cool-down zones. Also, the velocity profile is not assumed to be parabolic, but rather is determined from solution of the momentum balance. As mentioned above, this may be particularly important in situations where large radial temperature gradients exist, or when a large change in the number of moles occurs during reaction. The ability to extend this simulation from two (r,z ) dimensions to three (r,z, 0) is also demonstrated, with prediction of natural convection in the system. One disadvantage of the approach described here is that so many parameters are introduced, it is not possible to generate generalized charts for correction factors that may 0 1992 American Chemical Society
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 69 Table I. Comparison of Experimental Weight Percent of Species in Effluent with Simplified Kinetic Model Results Using the Plug Flow (One-Dimensional) Assumption" 1198 K,0.05 s 1198 K,0.03 s 1245 K,0.013 s 1245 K,0.023 s exptl model exptl model exptl model exptl model 26 26 10 10 20 20 6.0 5.8 C2H6 4.7 5.5 8.9 6.7 5.8 7.0 10.7 8.2 CH4 C2H4 60 59 66 66 62 62 66 67 CBH6 5.7 5.0 10 12 7.4 5.8 11 13 The contact times shown are the kinetically averaged reaction times.
be applied to simplified calculations for laminar flow reaction systems. Instead, this article presents a technique for detailed modeling of specific systems. This technique should be used when the simplifying assumptions presented by earlier workers in this area are inapplicable. The advantage of the approach described here is that it makes possible the detailed description of relatively complex systems, where radial variation of temperatures and of concentrations of multiple species, as well as multiple reaction steps, have impact on scaleup and performance of the system. One area of interest in laminar flow studies is hydrocarbon pyrolysis. Commercial units invariably operate in the turbulent flow regime. At the same time, it is desirable to study these systems in the laboratory where the scale of the system results in laminar flow. For example, in ethane pyrolysis, commercial coils that are 5-8 cm in diameter and 50-100 m in length may operate at Reynolds numbers of 100000-200000 (see, e.g., Robertson and Hanesian (1975)). On the other hand, a laboratory reactor with 0.5-cm i.d. and a 0.4-m-long reaction zone will have Reynolds numbers of the order of 100-400 (depending upon temperature and degree of conversion). The simulation described here is capable of simulating laboratory-scale ethane pyrolysis reactors. The ultimate goal of this exercise is to quantitatively describe the impact of radial concentration and temperature profiles on the chemistry of the system and to be able to then extrapolate the results to a commercial-scale unit. Results show that transport plays an important role in laminar flow ethane pyrolysis, especially at high reaction zone temperatures, e.g., 850 "C and above. Another purpose of this article is to demonstrate our ability to rigorously simulate transport in a laminar flow reactor, by application of the modeling technique to laboratory-scale ethane pyrolysis. Simulation Technique A laminar flow, tubular, ethane pyrolysis reaction system was simulated by use of the Phoenics package available through CHAM of North America (see Rosten and Spalding (1987)). The governing equations solved are conservation equations written for totaJ fluid flow, species flow, momentum and energy, of the form accumulation - net inflow by convection net inflow by diffusion = net generation (1) The exact form of the equations used is given in the Appendix. In the most general case, the Phoenics simulation is capable of solving the governing equations in un-steadystate mode in three dimensions, using Cartesian, cylindrical, or body-fitted coordinates. As applied here, two dimensions (r, z ) are used without time dependence. The governing PDEs (partial differential equations) are solved by using finite volume cella within Phoenics. This is accomplished by integrating the PDEs over a finite volume computational cell and representing the resulting volume averages by way of interpolating schemes. This
technique has similarities to both finite difference and finite element techniques. The result is a set of algebraic equations that are solved in either full-field or "slab-wise" mode. Species mass fractions, enthalpies, and pressures are represented at the center of the cells, while velocities are arranged in a staggered fashion at the center of cell faces. Boundary conditions are treated as source terms in the applicable finite volume cells. For instance, heat/momentum losses at the walls are treated as source terms in the cells adjacent to the walls. Interested readers are directed to Rosten and Spalding (1987), with further background on the finite volume approach given by Patankar (1980). Simulation Details for Ethane Pyrolysis Because use of the Phoenics simulator to study chemical reaction systems is a somewhat novel use of this simulator, and because the execution time of simulation is very sensitive to the representation of physical properties within the system, an outline of the details of simulation of the ethane pyrolysis system will be given here. For simulation of the tubular, laminar flow, ethane pyrolysis system, five chemical species were chosen. These are C2H6, CzH4, Hz, C6H6,and CHI. In this case, ethane, methane, and hydrogen are actual chemical species. C2H4 is a pseudospecies, used to represent both ethylene and acetylene, while c,& is a pseudospecies used to represent C3 and heavier species. The following simplified kinetic mechanism with Arrhenius temperature dependence was used to represent the reaction system: CzH6 2CH4 - Hz ko = 11 X 1013 l / s E, = 74 kcal/g-mol +
CzH6
-
-
ko = 1.6 X 1013 l / s ko = 9.6
CzH4 + Hz E, = 64 kcal/g-mol
CZH4 (1/3)C6& + "2 X 10IOl/s E, = 56 kcal/g-mol
Reactions were assumed to be irreversible and first order. The frequency factors and activation energies for this mechanism were regressed from experiments conducted at reaction zone tube temperatures of 900-1000 "C, with ethane conversions ranging from 74 to 94%. The experiments were carried out in 3- and 5-mm4.d. silicon carbide reactors, 67 cm in length. The reactors where heated by placing the center third of the length inside electric furnaces. The experiments were conducted at pressures just above atmospheric, with no dilution of feed ethane. Additional information on the experimental procedure is given by Velenyi et al. (1991). The kinetic parameters were regressed by use of the plug flow assumption, and the plug flow model results showed an average error of 10.3% for weight fraction of species in the effluent. A comparison of kinetic model and experimental results are given in Table I. Note that this simplified kinetic mechanism will be only approximate for two reasons. First, the actual mechanism of ethane pyrolysis is a complex set of serial and parallel
60 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
reactions involving hundreds of free radical and molecular species (e.g., see McConnell and Head (1983)). Second, the kinetic parameters used in this simplified mechanism have been regressed by use of the plug flow assumption, an assumption that will later be shown to be only approximate under these conditions. However, this simplified kinetic mechanism is sufficient to explore the impact of transport on the system. With the simplified mechanism, the reactions will be occurring on approximately the correct time scale with the correct type of temperature dependence, product distribution, and heat of reaction. Comments on use of more detailed kinetic descriptions will be made later. The heats of reaction for these three reactions were taken as -15.5, +32.7, and -5.9 kcal/g-mol reaction as written. As described in the Appendix, these are the heats of reaction at the standard state of 298 K and 1atm between gaseous species. Variation of heats of reaction with temperature are taken into account in the way that enthalpy and temperature are related. Density of the mixture was determined by using the ideal gas law, which is a very reasonable assumption for the temperatures and pressures under consideration. Because physical property calculations can account for a large portion of the computation time in simulation of chemical systems (e.g., see Fagley (1984)),shortcut methods were used for calculating viscosity, thermal conductivity, and molecular diffusivity. First, viscosity and thermal conductivity calculations were performed over the temperature range of 25-1000 OC using three representative mixtures corresponding to 0%,40%, and 80% ethane conversion. Values of these transport properties and mixing rules were taken from Perry and Chilton (1973). The following correlations, found to give less than 5% error over the temperature and concentration ranges of interest, were then used in the simulation: k = -0.023 + ~ ~ ~ ( 0 . 0 63.88 6 X 10-5T) 1.48 X 10-4T (2) p =
2.7567 X
+
+
lo4 + 2.4178 X
10-8T
(3) with Tin kelvin, k in J/(ms.K), and p in kg/(ms). Species specific heats were assumed to vary linearly with temperature. This facilitated determination of temperature from enthalpy, as compared to use of second-order temperature dependence of specific heat. Binary molecular diffusivities for each pair of species in the mixture were determined according to the method of Fuller, Schettler, and Gidding (see Fuller et al. (1966)). Subsequently, the expression shown in eq A6 of the Appendix was used to determine molecular diffusivity of each species into the mixture. During execution of the simulation, it was suspected that these mixture diffusivity calculations were consuming a large portion of the total computation time. In an effort to reduce execution time, a simulation run was made during which values of several sets of Dim where printed out along with temperature. On the basis of these values, simplified expressions for Dim were determined:
Dim = A,T1.75/P (4) where Dim is in m2/s, T is in kelvin, P is in atm, and A, is 1.6 X 2.5 X 8.3 X 10-lo,1.0 x and 5.1 x lo-'' for CH4,H2, C2H6, C2H4 and c6&, respectively. Two runs were then compared, one with the rigorous diffusivity calculation and another with shortcut diffusivity calculation. While results from the shortcut diffusivity calculations were very comparable to the more rigorous run (mass fractions differing by less than 3% over the length of the reactor and temperatures by less than 4 K), the shortcut method realized a 69% reduction in execution time! The
shortcut method for diffusivity determination was used on all subsequent runs. Within Phoenics, certain boundary conditions (in the form of net generation source terms) for eq 1 can be specified through the user input file. For example, the momentum sink term associated with a wall can be specified in the input file by setting the z velocity to zero at the outer face of cells adjacent to the wall. Similarly, inlet velocities, pressures, and species mass fractions may be set in the user input file. While there is allowance within Phoenics to specify a constant wall temperature, there is currently no mechanism for specifying a wall temperature profile in the user input file. For this reason, a source term had to be added to the enthalpy balance at cells along the wall in order to simulate the type of tube temperature profiles realized in ethane pyrolysis. This source term appears in a user-supplied subroutine, and has the form (5) where Shis the enthalpy source term in cells adjacent to the wall, J/(m3-s),k is the thermal conductivity in the cell, A, is the area of the cell face at the wall, d is the distance from the cell center to the wall (which may be corrected for differences between conduction in Cartesian and cylindrical coordinates), V, is the cell volume, Tw is the user-specified wall temperature at the z coordinate of the cell center, and T is the gas temperature in the cell. Also, Phoenica has no hilt-in feature to take into account the enthalpy associated with species that diffuse from one temperature to another as described in the Appendix. A user-supplied subroutine was also included in the simulation that adds an enthalpy source term determined from diffusion fluxes and the species' enthalpies. In the past, Phoenica has been used primarily to study fluid flow and heat transfer, and in some cases combustion processes. Using this simulator to study kinetically controlled reactors is somewhat novel and requires the addition of approximately 500 lines of FORTRAN code in the form of subroutines called by the simulator. To summarize, the user-supplied subroutines provide the following calculations: (1)physical property evaluation of density, thermal conductivity, molecular diffusivity, specific heat, and viscosity that vary with temperature and composition; (2) enthalpy balance source terms for representation of wall temperatures, molecular diffusion contributions, and heats of reaction; (3) species mass balance source terms for representation of chemical reaction; and (4) determination of gas temperature given enthalpy and composition. Having defined physical properties, kinetics, and the balance source terms described above, the user need specify only a few quantities to execute a simulation run: (1)tube inside diameter, length, and temperature profile; (2) inlet flow rate, composition, temperature, and pressure; and (3) numerical solution parameters such as the number of radial and longitudinal grids to be used, relaxation parameters, and convergence tolerances. The simulation then determines radial and longitudinal velocities, compositions, enthalpies, temperatures, and pressures for each finite volume cell in the system.
Simulation Results Velocity, Temperature, and Concentration Profiles. The model was used to simulate several experiments, as described in Velenyi et al. (1991). Reactor conditions for one set of runs are given as follows: tube i.d. = 5.0 mm; tube length = 67 cm; tube temperature profile varying linearly between 0 cm at 227 "C, 25 cm at 1007 OC, 50 cm
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 61
P 1:
”g
I
a,
-t
4 1
l v
i
1 L
B
DISTANCI,
1.1
CY
Figure 1. Temperature profiles from simulation of 5-mm4.d. tubular reactor with 72% conversion of feed ethane. Temperatures (“C): 0,tube; 0,gas at wall; A, mixing cup; V, centerline.
at 1007 “C, and 67 cm a t 127 “C; pure ethane feed rate = 29.7 cm3/s (STP); system pressure = 1 atm. The simulation was carried out with 50 longitudinal grids and 20 radial grids. The simulation results gave a final product with 72% conversion of feed ethane, 3.9% yield of methane, 4.8% yield of hydrogen, 57.6% yield of ethene, and 5.7% yield of benzene. Here, yield is used to denote grams of product per gram of ethane fed. The simulation results are very close to the actual product yields measured in the laboratory, with relative errors less than 5%. Figure 1 shows a plot of temperatures vs length along the tube for this simulation run. “Tube temperature” is the user-specified tube wall temperature. The high-temperature plateau is representative of the near-constant tube temperature realized within the furnace. The “gas at w d ” temperatures are the simulated gas temperatures in the outermost radial grid, whose grid center is located 0.125 mm from the wall. “Mixing-cup temperature” is the mixing-cup temperature of the fluid, i.e., the volumetricflow-averaged temperature of the fluid. The “centerline temperature” is the predicted temperature at the centerline of the reactor. This figure shows the large radial temperature gradients that can result when ethane pyrolysis experiments are run at elevated temperatures. Note that, at the beginning of the reaction zone, the centerline temperature is more than 300 “C lower than the gas near the wall! Further down the reactor, near the end of the reaction zone, the temperature difference drops to about 100 “C. Note that considerable reaction is still occurring in the first 5-10 cm of the cooling zone, as indicated by the centerline temperature. On the basis of Figure 1alone, we can see the difficulties that would arise in attempting to simulate the temperatures in the reactor with a one-dimensional model. If we take the simplified kinetic set used in the model to be the true kinetics for a moment, we can determine the level of conversion that would be obtained by using the temperature profiles shown in Figure 1. The results of this exercise are summarized: “actual” two-dimensional conversion = 72% ; one-dimensional conversion (tube temperature) = 99% ; one-dimensionalconversion (mixing-cup temperature) = 47 % ; one-dimensional conversion (centerline temperature) = 28%. These results highlight the difficulty of finding an appropriate temperature profile for use in one-dimensional modeling of the system. Laboratory-measured tube tem-
1.2
I
:
:
1.4
1.8
8.8
:
:
:
1.1
1.2
1.4
RADIUS,
: 1.1
:
: 1.1
: 1.1
f
: 2.2
2.4
2.0
YY
Figure 2. Radial profiles of temperature and longitudinal velocities at 27-em length in 5-mm4.d. reactor with 72% final ethane conversion. High-temperature case. (-0-1 Gas temperature; (101)z velocity. 1.17
Y @.It
i ‘i
A
R 0.6
RADIUS,
YY
Figure 3. Radial mass fraction profiles at 27-em length in 5-mm4.d. reactor with 72% final ethane Conversion. High-temperature case. (-O-) Methane X 10; (=O=)ethane; (--A.-)ethylene.
peratures or centerline temperatures will not give an adequate description of the system temperature. Even the mixing-cup temperature is not a good approximation of the kinetic, flow-averaged temperature of the gas because it does not adequately represent the slow-moving, rapidly reacting gas near the wall. Figure 2 shows temperature and longitudinal velocity profiles in the radial direction at a longitudinal distance of 27 cm, near the beginning of the reaction zone. The gas velocity is very close to the parabolic profile associated with nonreacting, isothermal laminar flow. Small deviations from exact parabolic flow occur because of nonequimolar reaction, as well as density and viscosity variation with temperature. Note the huge variation in temperature (and therefore reaction rates) with radius. The slowly moving material near the wall is near 980 “C, and reacting very rapidly. At the same time, the slowly moving material near the core is only at 660 “C, far below the temperature where ethane would react on this time scale. Figure 3 shows radial mass fraction profiles for ethane, methane, and ethylene at the same 27-cm length used in Figure 2. Note that at the centerline, where temperatures
4
62 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
"T
E C R
b 198 0.0
1 0.2
1 0,4
:
:
0.8
e,*
:
: 1.0
: 1.2
:
:
1.4
1.8
: 1.1
:
m
2.0
2.2
2.4
2.0
I
0
RADIUS.
Figure 4. Radial profiles of temperature and longitudinal velocities at 47-cm length in 5-mm4.d. reactor with 72% final ethane conversion. High-temperature case. (-04 Gas temperature; (=O=) z velocity.
::q
0.0
0.2
0.4
0.8
S.8
1.a
1.2
1.4
RADIUS,
6
0
YU
1.8
1.8
2.0
2.2
2.4
IS
10
20
U
26
DISTANCE,
40
IS
46
68
SS
CM
Figure 6. Temperature profiles from simulation of 6.5-mm4.d. reactor with 22% final conversion of ethane. Low-temperature case. (-O-) Gas at wall; (EO=) gas at midradius; (-.A...) gas at centerline; (- - 0 - -) wall temperature.
2.8
YY
Figure 5. Radial mass fraction profiles at 47-cm length in 5-mm-i.d. reactor with 72% final ethane conversion. High-temperature case. (-0-) Methane X 10; (=O=) ethane; (.-.A-) ethylene.
are far too cold for ethane to react on this time scale, the mass fraction of ethane is 0.96, indicating a local ethane conversion of 4%! This nonzero conversion of ethane is due to molecular diffusion of reactant away from the centerline toward the wall and diffusion of products away from the wall toward the centerline. This degree of diffusion is made possible by the diffusivity of these species at these elevated temperatures combined with the small distance (2.5 mm) from the centerline to the wall. Figure 4 shows temperature and longitudinal velocity profiles in the radial direction a t a longitudinal distance of 47 cm, near the end of the hot tube temperature plateau. Here the temperature profile is much less severe. The velocity at the centerline is much higher than that shown in Figure 2, due primarily to heating of the gas at the centerline. Figure 5 shows the radial mass fraction profiles at the same distance. Once again,the radial concentration gradients are made much less steep by the diffusion of reactant and products. On the basis of these figures, we can conclude that the largest source of errors in using the plug flow, one-dimensional assumption will be in treating the temperature
a.O 0.m
:
0.2
:
a.4
: 0.0
:
0.8
:
1.1
:
1.2
:
1.4
: 1.0
RADIUS,
:
1.8
:
2.S
;
2.2
:
2.4
:
1.0
:
2.1
:
8.0
zd 1.1
YY
Figure 7. Radial profiles of mass fraction, temperature, and longitudinal velocity for 6.4-mm4.d. reactor with 22% final ethane conversion. Low-temperature case. (-0-1 Gas temperature; ( = O r ) velocity, m/s; (.-O.-) mass fraction of C2Hs.
as a function of longitudinal distance only. While some additional error will occur due to radial variations in concentration, these will be of a much smaller magnitude. Figure 6 summarizes results for a lower temperature simulation run. Here, the maximum tube wall temperature is 860 "C. Note that the difference between the tube and centerline gas temperatures is much smaller here than in the higher temperature case shown in Figure 1. Because the radial gradients are much smaller at this lower operating temperature (e.g., 20-40 O C from centerline to tube wall), it will be much easier to represent temperatures in a one-dimensionalsimulation. Figure 7 shows temperature, ethane, and longitudinal velocity profiles in the radial direction midway through the reaction zone. Note that velocities are about an order of magnitude lower than in the previous case. With lower velocities, the rate of radial heat transfer becomes faster relative to the rate at which material moves through the system, resulting in smaller radial temperature gradients. Note that radial diffusion
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 63 Table 11. Impact of Tube Diameter, Conversion, and Operating Temperature on Temperature Gradient and One-Dimensional Simulations for Ethane Pyrolysisa I-D conv, % ~~
tube T , O C 775
tube i.d., mm 3 5
max T, - Tcl,O 10 5 4 13
C
11
850
3 5
1007
3 5
8 28 16 11 61 38 19 381 129 106 544 206 156
2-D conv, % 51.0 69.0 95.9 49.3 71.2 93.4 46.5 73.8 95.5 36.8 67.3 93.5 41.0 72.9 97.4 31.2 70.9 86.0
tube T 51.4 69.2 96.1 52.3 73.5 93.6 53.5 76.3 96.3 50.6 76.3 95.6 82.4 92.5 99.8 86.9 98.6 99.7
TR/2
47.3 68.0 95.5 47.2 70.9 92.4 42.0 68.1 95.1 24.9 56.8 90.9 17.0 47.2 95.8 4.1 41.5 67.2
Tcl 46.2 67.0 95.4 45.4 69.7 92.1 38.9 65.8 95.0 20.4 51.9 89.5 9.8 38.4 94.1 2.0 27.2 57.6
"Tube T" is the tube temperature in the hot zone, while "tube i.d." is the inside diameter of the tube. The column headed "mas T, - Tc," is the difference between the wall temperature and the centerline gas temperature a t the beginning of the reaction zone. "2-D conv" is the "correct" conversion given by the full two-dimensional model. "I-D conv" are the approximate, one-dimensional conversions obtained using the same kinetics and three different temperature profiles obtained from the two-dimensional run: "tube T" uses the tube wall temperature uses the temperature profile midway between the centerline and wall, and "TcY uses the centerline gas temperature. profile, "TR12"
of material is quite fast in this case relative to the rate at which material moves through the system, resulting in very small radial ethane concentration gradients. As seen from these two simulation cases, the impact of transport on the system will vary with temperature of operation. At high temperatures, very high velocities are required to obtain a given level of conversion. At high velocities, material is moving through the system faster than the rate of radial heat transfer, and cold gas can exist near the center of the tube well into the reaction zone. At lower temperature, lower gas velocities are required to obtain a given level of conversion. At low gas velocities, material is moving through the system more slowly than the rate of radial heat transfer, and the gas near the centerline is close in temperature to the gas at the wall. This means that characteristic temperatures for one-dimensional simulation will be easier to obtain at lower temperature operation, and that the one-dimensional simulation will be a more accurate representation of the physical system. A study of the effects of operating temperature, tube diameter, and degree of conversion on radial gradients is given in the next section. Model Results: Impact of Operating Parameters. A series of simulation runs were carried out to determine the impact of various operating parameters. Two tube diameters were used, 3- and 5-mm i.d. These are the smallest diameters that may easily be used, while larger diameter tubes will give larger radial gradients (up to the point where natural convection may start-see following section). Three different maximum tube temperatures were used: 1007,850, and 775 "C. The tube temperature profiles are similar to the ones used in the previous section-a 70-cm-long reactor with the middle third held at a constant temperature. In addition, inlet velocity was adjusted in each case to obtain three different conversion levels: one below 50%,one near 70%, and one near 90%. Five quantities are reported for each simulation run in Table 11. The first is the difference between the tube wall temperature and the centerline gas temperature at a distance of 28 cm, approximately onefifth of the way through the hottest section of the tube. The lower this temperature difference, the less radial temperature profiles will affect the accuracy of a one-dimensional simulation of the system.
The second quantity reported is the ethane conversion level obtained from the two-dimensional simulation. The third through fifth quantities reported are the ethane conversion levels obtained with a one-dimensional simulation, and the same kinetic model, using three different temperature profiles obtained from the two-dimensional simulation: the tube temperature, the gas temperature at the centerline, and the gas temperature midway between the centerline and the wall. The differences between the two-dimensional simulated conversion and the three one-dimensional simulated conversions are an indication of the type of inaccuracies to be expected when a onedimensional, plug flow treatment of the kinetic data is used. Several conclusions may be drawn from this table. The radial temperature gradient in the reaction zone increases dramatically with increasing reaction temperature. At higher operating temperatures, the gas velocity through the reactor must be higher to achieve a given level of conversion. As the gas temperature velocity becomes higher, the rate at which the material moves through the reactor becomes faster than the rate at which heat is radially transferred, resulting in colder gas near the centerline. Similarly, for a given tube temperature, conversion decreases with increasing gas velocity. Lower conversions give rise to larger radial temperature gradients as a result. Changing from a 5-mm4.d. tube to a 3-mm-i.d. tube helps decrease the radial temperature gradients, but not completely. In those cases where large radial temperature gradients exist, the results of the one-dimensional simulations yield noticeably different predicted ethane conversions. At a tube temperature of 1007 "C, even the 3-mm4.d. results show the plug flow assumption to be inadequate. These are the cases where the plug flow assumption is unacceptable, and using the one-dimensional assumption will give erroneous kinetic parameters. Note that conversion numbers near 100% can be misleading. For example, predicted conversions of 97.4% and 99.8% may seem to be fairly close. However, one represents 2.6% unreacted, while the other represents 0.2% unreacted. For a first-order decomposition, the reaction rate needs to be increased by roughly 70% to go from a fiial predicted conversion of 97.4% to 99.8%! This means
64 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 Table 111. Longitudinal Velocity Profiles in the Heating Zone for a 3-mm4.d. Tube with a Maximum Tube Temwrature of 1123 K" radial grid u 8 ~m/s , %par, m/s 20 0.062 0.061 0.180 0.181 19 0.292 0.293 18 0.397 0.398 17 0.497 0.498 16 0.590 0.591 15 0.677 0.678 14 0.758 0.759 13 0.833 0.833 12 0.901 0.902 11 0.963 0.964 10 1.02 1.02 9 1.07 1.07 8 1.11 1.11 7 1.15 1.15 6 1.18 1.18 5 1.21 1.21 4 1.22 1.22 3 1.237 1.237 2 1.243 1.243 1
Table IV. Longitudinal Velocity Profiles in the Heating Zone for a 6-mm4.d.Tube with a Maximum Tube Temperature of 1223 K"
"uZcim are simulation results, while u , , is ~ a parabolic velocity having the same uz,l as the simulation results. Tz0= 818.5 K,TI = 792.0 K.
uZrimare simulation results, while u , , is ~ ~a parabolic velocity ~ simulation results. T,,= 509.5 K,TI = having the same u ~as, the 317.5 K.
that difficulty will be experienced in determining accurate kinetic parameters from measured temperature profiles and the one-dimensional assumption at lo00 "C and high conversion in a 3-mm-i.d. tube. On the basis of the results of Table 11, we can define limits where plug flow treatment of kinetic data are acceptable. For example, if we require that radial temperatures gradient must be less than 20 "C, then we should not operate a 3-mm4.d. tube at temperatures of 850 "C with conversions less than 70%. As shown in the table, none of the results with tube temperatures of 1007 "C should be handled using the plug flow assumption. Velocity Profiles: Deviation from Parabolic. Two runs of the simulator were made to determine sensitivity of velocity profiles to reactor conditions. In the first run, a reactor tube of 3-mm4.d. by 0.67-m length was simulated with a maximum tube temperature of 1123 K extending over the central 0.25 m of the tube. The tube temperature at the inlet was taken to be 500 K. The average inlet gas velocity was set at 0.3 m/s at a temperature of 390 K. Table I11 shows the z velocity profile determined from the simulation in the heat-up zone at z = 0.34 m. Also shown in this table are the calculated parabolic velocities, based on a match of ul,sim. In fact, for these conditions the simulated velocity profile remains very close to parabolic throughout the reactor. The close agreement between the simulation results and the calculated parabolic profile indicates that the parabolic assumption would be an excellent approximation under these reactor conditions. Note that there is only a 26.5 K temperature difference between the fluid near the wall and that at the centerline. For this run, simulated radial velocity profiles are zero to within the accuracy of the simulation. In the second run, a reactor tube of 5-mm i.d. by 0.762-m length was simulated with a maximum tube temperature of 1223 K extending over the central 0.28 m of the tube. The tube temperature at the inlet was taken to be 373 K. The average inlet gas velocity was set to 1.08 m/s, at a temperature of 300 K. Table IV shows the z velocity profile determined by the simulator at distance z = 0.048 m. Comparison with z , , reveals ~ ~ that this profile is significantly different from parabolic. The fluid near the wall is nearly 200 K hotter than the fluid at the centerline.
There is a corresponding radial variation in fluid viscosity and density resulting from this large radial temperature gradient. Also, at the z location shown in Table IV,there is an appreciable negative radial velocity (averaging 0.003 m/s). Note that the effect of the temperature gradient on the velocity profiles is somewhat complex. The temperature dependency of viscosity will impact the shape of the longitudinal velocity profile. Lower density at the wall indicates gas expanding more rapidly at the wall than near the centerline, which will tend to increase velocity along the wall and push gas inward toward the centerline. The data shown in Table IV are the cumulative result of these effects. For this second run, the velocity becomes nearly parabolic toward the end of the reaction zone as the radial temperatures begin to equilibrate. Then in the cooling zone, the opposite effect occurs: fluid near the wall is significantly cooler than fluid near the centerline, resulting in velocities at the wall being lower than those predicted by a parabolic profile. Correspondingly, there is a positive radial velocity in this region. Note that accurate description of u, near the wall is very important for accurate representation of reaction in the system. This is because a good portion of the reaction occulg near the wall where the slowest moving, hottest fluid exists. These results show the importance of solving the momentum balances to determine the velocity profiles, instead of making the parabolic assumption, as tube diameters become larger and tube temperatures hotter. Larger diameters give larger radial temperature gradients for the obvious reason that heat must conduct over a longer distance to reach the fluid near the centerline. Higher tube temperatures give larger radial temperature gradients for two reasons. First, with hotter tube temperatures, the cold inlet gas must be heated to a higher temperature in the reaction zone. Second, and most important, higher reaction temperatures require higher gas velocities to achieve a given level of conversion in a given length reactor. This means that, at higher tube temperatures, we are exposing the gas to a more rapid tube temperature rise as the gas flows down the tube, resulting in larger radial temperature gradients.
radial grid 20 19 18 17 16 15 14 13 12 11 10 9 8 7 6 5 4 3 2 1
m/s 0.170 0.482 0.765 1.02 1.25 1.45 1.62 1.77 1.89 1.99 2.07 2.13 2.18 2.21 2.24 2.26 2.28 2.29 2.296 2.299
Uzairn,
uz,par* m/s
0.114 0.332 0.539 0.735 0.919 1.09 1.25 1.40 1.54 1.67 1.78 1.88 1.98 2.06 2.13 2.18 2.23 2.26 2.288 2.299
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 65
Impact of Natural Convection in Horizontal Reactors. Because of the large radial temperature gradients occurring in laminar ethane pyrolysis experimenb at elevated temperatures, the possibility of natural convection within the system arises. At the beginning of the reaction zone, cold, dense gas exists near the centerline of the reactor, while hot, less dense gas is present at the wall. For a horizontal tube, two natural convection cells will be set up, symmetric about a vertical plane bisecting the reactor along its length. The cold, dense gas near the center will be drawn downward, while the hot, less dense gas near the wall will rise upward along the wall. As a result, streamlines will become helical. Looking down the reactor from inlet to outlet, two vortices will be set up in the heating section. The vortex on the left will be moving clockwise, while the vortex on the right will be moving counterclockwise. Suppose that, near the end of the reaction zone, only small radial temperature gradients exist. As the gas enters the cooling section, the colder, more dense gas will be near the wall and the hotter, less dense gas near the centerline. This will cause the vortices to move in the opposite direction. (Note that the flow profiles in the cooling zone resulting from natural convection will become more complex for the case where some of the gas near the centerline is still appreciably cooler than the gas at the wall near the end of the reaction zone.) The vortices resulting from natural convection will have a characteristic time constant, i.e., the time required for the vortex to complete one “revolution”. If the characteristic time the gas is in the reactor is considerably shorter than the time scale of a vortex revolution (e.g., 2 or more orders of magnitude smaller), then natural convection will have virtually no influence on the system and the twodimensional representation of the system described above will be adequate. However, as the time the gas is in the reactor approaches or exceeds the time scale of a vortex revolution, natural convection will play a very important role in transfer of mass and heat within the system. Under these circumstances, a two-dimensional representation of the system will be inadequate, and a three-dimensional representation that takes into account gravity forces on the system will be required. In the beginning of the reaction zone, the natural convection vortices will tend to take hot, more fully reacted gas away from the walls and cause it to swirl down into the colder core. Also, cold, less fully reacted gas from the center will move down toward the bottom of the reactor and then swirl up along the sides of the reactor. The net effect will be increased heat transfer from the wall to the gas and smaller radial gradients in both temperatures and concentrations. The amount of natural convection within the system will be affected by several factors. First, the more viscous the gas, the smaller the tendency for natural convection. In the limit, if the gas were “infinitely viscousn,natural convection could not occur. Second, as radial temperature gradients become more pronounced, the difference in density and therefore gas buoyancy becomes larger, resulting in increased tendency for natural convection. Last, as the tube becomes larger in diameter, the drag along the wall per cross section of gas becomes smaller, and natural convection becomes more likely. Suppose we are studying an ethane pyrolysis system at constant conversion, say 70%. As we increase the wall temperature to conduct experiments at higher and higher temperatures, larger longitudinal gas velocities are required to maintain our 70% conversion level. As shown in Table 11, as tube temperatures become hotter, radial thermal gradients become larger. Thus at high temperature op-
Table V. Effect of Natural Convection in a 5-mm4.d. Tube” final 2-D final 3-D conversion, conversion, tube temp, K % 5% max T diff, K 1123 93.5 93.5 4 1123 67.3 67.3 4 1280 31.3 31.2 16 1280 71.0 71.0 15 “The “max T diff“ is the maximum difference between two-dimensional and three-dimensional simulated temperatures. This maximum difference occurred near the beginning of the reaction zone in each case. Further down the reaction zone, temperatures were different by less than 1 K in each case. Table VI. Importance of Natural Convection in a 41-mm4.d.Tube with a Hot Zone Tube Temperature of 1173 Ka tube i.d., inlet vel, angular T mm mls outlet C2H6 variation, K 41 0.04 0.069-0.570 146 41 0.10 0.309-0.649 280 ““Inlet vel” is the inlet velocity of the gas at 386 K. “Outlet CzHsngives the mass fraction of unreacted ethane at the outlet. The variation of unreacted ethane across the cross section is the result of both radial heat conduction and natural convection “swirling” of flow profiles. ”Angular T variation” is the angular variation of temperature at a constant radial coordinate equal to half the inner radius of the reactor, at a reactor length corresponding to the beginning of the reaction zone.
eration, natural convection vortices will have shorter revolution times. At the same time, for high-temperature operation, gas will be moving through the system faster in order to maintain 70% conversion. For this reason, it becomes necessary to test the importance of natural convection over a temperature range. Similarly, it is necessary to test the importance of natural convection over a range of conversions, because lower conversions have larger radial temperature gradients but higher longitudinal gas velocities at a given tube temperature. Determining the impact of tube diameter on natural convection is not as difficult because as tube diameter increases, both thermal gradients and the volume-to-surface ratio of the system increase. Each of these factors contributes toward more natural convection. Thus, once we determine the tube diameter at which natural convection just begins to become important, we will know that all smaller diameters are unaffected, while all larger diameters are affected. Toward the end of determining the importance of natural convection (and therefore the accuracy of the twodimensional model described above), a three-dimensional model of tubular ethane pyrolysis was developed with the use of the Phoenics simulator. In addition to the radial and longitudinal coordinates r and z employed before, an angular coordinate 0 was introduced. A gravity source term was added to the mechanical energy equation to account for a horizontal tube. Because of symmetry considerations, only ?r radians of the tube were simulated. Several tube temperatures, geometries, and conversion levels were simulated. The results of these runs are summarized below in Tables V and VI. Table V shows that the impact of natural convection in the 5-mm-i.d. tube is minimal. While some angular variation in temperature exists between the two-dimensional case without natural convection and the three-dimensional case with natural convection, these temperature variations are localized near the beginning of the reaction zone, and the overall impact on conversion and selectivity is minimal. These simulation runs validate the use of the two-dimen-
66 Ind. Eng. Chem. Res., Vol. 31, No.1, 1992
-
W
4
: L . I M - 0 1 ./I.
PYROLYSIS OF C2H6
- M T CONU - SONG E*
I PHOENICS
Figure 8. Velocity field at constant z for three-dimensional simulation of a 41-mm4.d. tube. Inlet velocity (386 K) is 0.10 m/s. Maximum tube temperature is 1173 K.
I PHOENICS - VAT CONU - SONG E* Figure 9. Mass fraction of unreacted ethane at constant z for three-dimensional simulation of a 41-mm4.d. tube. Inlet velocity (386 K) is 0.10 m/s. Maximum tube temperature is 1173 K.
1
PYROLYSIS OF C2H6
sional simulation studies described earlier in this report. Table VI shows that natural convection is significant in the 41-mm-i.d. reactor, assuming laminar flow. This is indicated by the large angular temperature variations shown. Note that, in this diameter tube, the outlet mass fraction of ethane varies significantly across the cross section of the reactor. This is the result of radial conduction of heat accompanied by natural convection “swirling” of the flow profiles. The large gradients in both temperature and species concentrations existing in these large diameter tubes make them less suitable for kinetic studies than the smaller diameter tubes. Figure 8 is a plot of velocity vectors in the beginning of the cooling zone in the 41-mm4.d. reaction tube. The cold, more dense gas near the centerline is moving downward, as the hotter, less dense gas near the wall rises. Figures 9 and 10 show how the natural convection vortices alter concentration and temperature profiles within the reactor. The coldest and least reacted gas has drifted downward from the centerline of the tube due to natural convection. Comments and Conclusions 1. Many previous workers have developed simulations to describe transport of mass and energy in reacting, laminar flow systems. The present work is an important extension of the previous work in that momentum transport is also considered. As shown in this paper, full so-
x
- MRT CM(U - SONG EW
1 PHOENICS Figure 10. Gas temperatures at constant z for three-dimensional simulation of a 41-mm4.d. tube. Inlet velocity (386 K) is 0.10 m/s. Maximum tube temperature is 1173 K. PWOLKIS OF C2H6
lution of the governing momentum balance equations, in place of assuming a parabolic velocity profile, is required for accurate simulation at elevated temperatures for tubular ethane pyrolysis. In addition, the current work shows that simulation is still possible when many of the simplifying assumptions made by previous workers are removed. Thew simplifying assumptions include having only two or three species of interest, a single reaction, constant tube wall temperatures, and, in many cases, constant physical properties. 2. This work demonstrates that rigorous simulation of transport in laminar flow systems is now possible without prohibitive computing expense. A typical two-dimensional simulation described here required less than 3 CPU min of execution time on a VAX 8800. Extending the simulation to three dimensions in order to study natural convection in the system was also possible with minimal manpower, and these simulations required less than 15 CPU min each. 3. The possibility of employing serial/parallel reaction mechanisms involving multiple species within the framework of the Phoenics simulator, with rigorous solution of the governing mass, energy, and momentum balance equations has been demonstrated. Simulation with detailed kinetics, involving dozens or even hundreds of reactions and species, could also be accomplished in the following manner: Beginning with a simplified kinetic set, the simulation can be used to determine velocity and temperature profiles. Subsequently, these profiles can be fixed, and used with a more detailed kinetic set in order to facilitate more rigorous chemical description without prohibitive computing expense. This technique will yield accurate results, provided that the simplified kinetic description gives compositions sufficiently accurate to characterize physical properties and heats of reaction of individual steps. 4. The parametric studies described here show the importance of transport in laminar ethane pyrolysis and the need for two-dimensional simulation including solution of the momentum balances. This is particularly true at elevated temperature, the direction being pursued by the new “millisecond furnace” technology of Kellogg (see,e.g., Ross (1983)). In future work, the Phoenica simulator would be especially useful in describing temperatures, velocities, and compositions near the wall surface during carbon formation and deposition and in scaling up into the turbulent regime.
Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992 67 Table VII. Identity of 4 and Source Terms within Phoenics for Various Conservation Equations conservation eauation variable d source terms S, total continuity unity mass i formed by reaction per mass species i x,, i = 1-n volume-time momentum ui, i = 1, 2 , 3 pressure gradient term (and gravity term) energy h viscous heating, pressure variation with time ~~~
Nomenclature Ai = parameter for diffusion coefficient correlation,component i, atmm2/(s.K1.75) A, = area of cell face at wall, m2 C,, = specific heat, J/(kgK) Dij = binary diffusivity of species i in species j , m2/s Dim= diffusivity of species i in mixture, m2/s d = distance from cell center to wall, m E, = activation energy, kcal/g-mol gi = total body force per unit mass of component i, m/s2 h = specific enthalpy, J/kg hi = partial mass enthalpy of species i, J / k g i AHrxn,J = heat of reaction j, J/kg of rxn j ji = mass flux vector of species i relative to mass-average velocity, kg/ (m2d k = thermal conductivity, J/(ms.K) ko = reaction rate frequency factor, l / s P = absolute pressure, atm R = radius of tubular reactor, m Rj = rate of reaction j, kg of rxn j/(m3.s) r = radial coordinate, m S, = source term for variable 4, units of 4 times kg/(m3.s) Sh= enthalpy source term, J/(m3.s) T = absolute temperature, K T,, = centerline temperature Tg= gas temperature T,,= temperature at radial grid n TRI2= temperature at r = R / 2 T, = tube wall temperature t = time, s V , = volume of cell, m3 v = mass-average velocity vector, m/s ui = component of velocity vector in direction i u, = component of velocity vector in direction z x i = mass fraction of species i, kg of i/kg of total yi = mole fraction of species i, kg-mol of i/kg-mol of total z = longitudinal coordinate, m re = diffusion term coefficient for variable 4, kg/(ms) 0 = angular coordinate, rad p = viscosity, kg/(ms) p = gas density, kg/m3 pi = mass concentration of species i, kg of i/m3 4 = balance equation variable, v for momentum balance, h for enthalpy, xi for mass balance, unity for continuity 7 = viscous stress tensor, N/m2 v. = divergence operator, l/m v = gradient operator, l / m Appendix: Derivation of Conservation Equations The governing conservation equations used within Phoenics are of the form (see Rosten and Spalding (1987))
Term by term, this equation states that accumulation minus net flow in by forced convection minus net flow in by diffusion is equal to net source generation. Table VI1 shows the identity of 4 for the various conservation equations used in the simulation.
Values of the coefficient for the diffusive term, r are also dependent upon which conservation equation is geing considered. For the energy equation, r, is k/C,. For the momentum equation, r, is viscosity I.L. For the species i balance equation, I', becomes pDim,where Dimis the diffusivity of component i in the mixture. In Phoenics, the source terms shown in Table VI1 for pressure gradient, viscous heating, and pressure variation with time are generated automatically. The generation terms for formation of species by reaction are supplied by the user in the form of FORTRAN subroutines. In order to put the governing conservation equations into this form, several assumptions are required. Some of these assumptions have been deemed invalid for a multicomponent gas system. As noted below, inappropriate assumptions have been overridden in the current work by setting coefficients r, to zero and by introduction of extra source terms S,. Total Continuity. According to Bird et al. (1960), the total continuity equation in a multicomponent system may be written as aP
- + (v.pv) = 0 at
Referring to eq Al, by setting 4 to unity and S , to zero, the form of the total continuity equation used within Phoenics is shown to be correct for multicomponent systems. Species Continuity. From Bird et al. (1960)) the continuity equation for species i may be written as
DPi - - - - p i ( ~ . ~- )(v.ji)+ ri Dt The quantity ji is the mass flux of species i relative to the mass average velocity, v, and is related to the gradient of the mass fraction according to Fick's f i s t law through the effective diffusivity of the mixture (Bird et al., 1960): ji = -pDimvxi (A4) Combining eqs A2, A3, and A4, and noting that pi is equal to pxi, yields
Comparing eqs A5 and A1 shows that the Phoenics form is correct, by setting @ equal to mass fraction xi, I?, equal to pDim,and source term S, equal to ri. The effective diffusivity is difficult to determine in the general case. In this work, the following approximation is used: for summation over all j unequal to i. As stated in section 18.4 of Bird et al. (1960),this approximation becomes exact when either all Dij are equal or when all species except i move at the same velocity. Momentum Balance. The momentum balance for a multicomponent mixture is given by eq 18.3-2 in Bird et al. (1960) as
For a two-dimensional (r, z ) simulation in a horizontal tube, the gravity term is neglected. Referring to eq A l , if we set 4 equal to v, and source term S , equal to vp, the Phoenics form is correct provided that V.7 =
-v.(pvv)
(A81
68 Ind. Eng. Chem. Res., Vol. 31, No. 1, 1992
According to Bird et al. (1960), this is the assumption required to develop the Naviel-Stokes equations from the general momentum balance. Equation A8 is exactly correct under the following circumstances; Newtonian fluid, constant viscosity, constant (V-v)and negligible “second viscosity coefficient”. Note that constant (v-v) is satisfied if density is constant along a streamline, or if flow is incompressible. There are also cases in which eq A8 is approximately correct. According to Cham of North American (Rawnsley, 1990), the approximation shown in eq A8 is commonly used in fluid flow studies in cases where flow is dominated by convective, as opposed to viscous, effects. Cham has performed a series of studies using the full momentum equation with VT and the “linearized” They have found that form replacing V - rwith -v.(~vv). the linearized form yields results that are virtually identical with the full form when the cell Peclet number for momentum exceeds 10. The cell Peclet number is the ratio of convective forces to viscous forces within a cell, and for momentum is defined as ( p u , L / p ) , where L is the cell length in the z direction. For the ethane pyrolysis studies conducted here, cell Peclet numbers were calculated at a variety of cell locations. In general, the Peclet numbers will be smallest for cells near the wall, where u, is small. For the case of 20 radial grid spacings, the smallest cell Peclet numbers for momentum were calculated as 50, while cell Peclet numbers for momentum at the centerline of the tube were found to be of the order of 1200-1500. On the basis of these results, the assumption shown in eq A8 has been assumed to be applicable for the simulations described in this paper. Energy Balance. According to Bird et al. (1960), the energy balance for a multicomponent system can be expressed as
do + v.(pvh) + v.q = r:vv + c(ji-gi) at
(~9)
Note that the energy balance given by Bird et al. contains no heat of reaction term. This is because their enthalpies are based on elements, not molecules, at a standard state, such that heats of reaction are “automatically” accounted for in the calculation of temperature from enthalpy. In the current work, enthalpies are based on molecules at the standard state of gaseous species at 298 K and 1atm. For this reason, a source term of the form CRjAHr,m,jfor reactions j = 1,2, and 3 was added to the Phoenics energy balance, where R j is the kilograms of reaction occurring per cubic meter per second, and A“,,, is the heat of reaction at the standard state in terms oijoules per kilogram of reaction j occurring. The energy flux q shown in eq A9 is the sum of energy fluxes due to heat conduction, interdiffusion of molecular species, and the Dufour effect. The Dufour effect is typically unimportant (Bird et al., 1960), in which case V.q v * ( - k ~ T x h i ( - ~ D i ~ ~ ~ i (A10) )) For the energy balance equation, the Phoenics simulator does ”automatically” take into account the double dot product of r and VV, which is the viscous heating term. However, referring to eq A1 with 4 set to h, the Phoenics form does have two shortcomings for proper representation of the energy balance for multicomponent systems. First, there is no “automatic” source term to represent interdiffusion of molecular species. This is the second term on the right-hand side of eq A10. A user-supplied source term of this form was added to the Phoenics simulation, and its contribution was found to be significant. In a test run of the simulation for tubular ethane pyrolysis, temperatures within the reaction zone varied by as much as 12 K
+
when this source term was turned off and on. The second difficulty of the Phoenics form of the energy equation for representation of multicomponent systems is in attempting to represent the temperature gradient term with the enthalpy gradient divided by C,. The enthalpy of the system is a function of temperature, pressure, and mass fraction of species present. The gradient of the enthalpy can be represented as
Vh(TQ,X1,XZ,...,2,) = (dh/dP)~P+ (dh/dT)vT
+ C ( d h / & i ) ~ ~( Ai l l )
For the simulations described here, pressure gradients are small, and enthalpy is only a weak function of pressure, so that the first term on the right-hand side of eq A l l may be omitted. By definition, dh/dT at constant pressure and composition is the specific heat, C,. However, kvT cannot be set equal to (k/C,)vh because the last term on the right-hand side of eq A l l is appreciable. For this reason, rhwas set to zero for this equation within Phoenics and a user-supplied source term of v-(-kvT) was added to the energy equation. Summary. For the simulations described here, the total and species continuity equations within Phoenics were used in the original Phoenics form, with diffusivity of species i into the mixture determined by eq A6 and a user-specified reaction term added. The momentum balance was also taken as the original form, noting that the approximation shown in eq A8 is employed. The Phoenics energy balance equation was altered in three ways. First, the term involving the gradient of enthalpy was turned off. Second, the two source terms shown in eq A10 were added. Third, a heat of reaction term was included. Registry NO.C2&, 74-84-0; CHI, 74-82-8; CzH4,7485-1; C& 71-43-2.
Literature Cited hdersen, T. S.;C o d , J. Evaluation of Models for Tubular, Laminar Flow Reactors. AIChE J . 1090,16, 542-553. Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 1960; pp 495-625. Cleland, F. A.; Wilhelm, R. H. Diffusion and Reaction in ViscousFlow Tubular Reactors. AICHE J. 1956,2, 489-497. Fagley, J. C. Flexibility and Efficiency in Modular, Dynamic, Chemical Plant Simulation. Ph.D. Dissertation, The University of Michigan, Ann Arbor, 1984. Fuller, E. N.; Schettler, P. D.; Gidding, J. C. A New Method for Prediction of Binary Gas-PhaseDiffusion Coefficients. Ind. Eng. Chem. 1966,58,19. Hsu, C. A Method of Solution for Mass Transfer with Chemical Reaction Under Conditions of Viscous Flow in a Tubular Reactor. AIChE J. 1965,Il, 938-940. McConnell, C. F.; Head, B. D. Pyrolysis of Ethane and Propane. In Pyrolysis: Theory and Industrial Practice; Albright, L. F., Crynes, B. L.. Corcoran,W. H., Eds.: Academic Press: New York, 1983; Chapter 2. Patankar. S. V. Numerical Heat Transfer and Fluid Flow: HemiSphere’Publishing: New York, 1980.’ Perry, R. H.; Chilton, C. H. Chemical Engineers’ Handbook; McGraw-Hilh New York 1973; pp 3.210-3.215. Rawnsley, S. M. of Cham North America, Huntsville, AL, personal communication, 1990. Robertson, R. W. J.; Hanesian, D. An Optimization Study of the Pyrolysis of Ethane in a Tubular Reactor. Ind. Eng. Chem. Process Des. Dev. 1975, 14, 259. Ross, L. L. Pyrolysis Furnace Design: Conventionaland Novel. In Pyrolysis: Theory and Industrial Practice; Albright, L. F., Crvnes. B. L.. Corcoran.W. H.., Eds.:, Academic Press: New York. 19i)3; Chap& 13. Rosten, H. I.; Spalding, D. B. The Phoenics Beginner’s Guide; Cham Ltd.: London, 1987: Chapters 1 and 2. Rothenberg, R. I.; Smith, J. M. Heat Transfer and Reaction in Laminar Tube Flow. AIChE J. 1966, 12,213-220.
Ind. Eng. Chem. Res. 1992,31,69-75 Sundaram, K. M.; Froment, G. F. A Comparison of Simulation Models for Empty Tubular Reactors. Chem. Eng. Sci. 1979,34, 117-124. Trombetta, M. L.; Happel, J. Analysis and Design of Gas Flow Reactors with Applications to Hydrocarbon Pyrolysis. AZChE J. 1970,11, 1041-1050.
69
Velenyi, L. J.; Song, Y.; Fagley, J. C. Carbon Deposition in Ethane Pyrolsyis Reactors. Submitted for publication in Znd. Chem. Eng. Res., 1991.
Received for review August 1, 1990 Accepted February 14,1991
KINETICS AND CATALYSIS Reaction Pathways in Lubricant Degradation. 3. Reaction Model for n -Hexadecane Autoxidation Steven Blaine and Phillip E. Savage* Department of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109-2136
We have developed a reaction model that describes the autoxidation of n-hexadecane under severe conditions. The model rate equations were solved numerically concurrent with optimization of the reaction rate constants to yield the solution that provided the best correlation with experimental data. This numerical analysis showed that hydroperoxides were formed in the reaction step with the smallest reaction rate constants, and they decomposed in the steps with the largest rate constants. Thus, hydroperoxides were key reaction products. A comparison of the kinetics determined in this study with those reported in the literature showed that the rate constant for hydroperoxide formation during the initial stages of the autoxidation compared favorably with that determined by others. The reaction rate constants for the secondary reactions determined for the overall autoxidation under severe conditions, however, differed from those appearing in the literature for the initial stages of paraffin autoxidation. Thus, the kinetics of paraffin autoxidation appeared to depend on the extent of oxidation. An illustration of the implications of these findings to lubricant degradation is also presented.
Introduction The degradation of a liquid lubricant is frequently accompanied by deleterious changes in its physical and chemical properties that can adversely affect ita performance. Oxidation of the lubricant base oil has been identified as the primary agent of degradation (Fenske et al., 1941; Korcek et al., 1986; Gunsel et al., 1988; Naidu et al., 1984,1986), and this has motivated many investigations of base oil oxidation (Colclough, 1987; Jette and Shaffer, 1988; Tseregounis et al., 1987; Naidu et al., 1984, 1986; Korcek and Jensen, 1976; Dialmond et al., 1952; Spearot, 1974; Hsu et al., 1986). Petroleum base oils are complex mixtures of hydrocarbons, however, and heteroatom-containing compounds that can function as prooxidants or antioxidants are also often present. This complexity has frustrated resolution of the reaction fundamentals for the oxidation. One alternate approach for resolving the fundamentals of base oil oxidation is to study the oxidation of a simpler yet chemically relevant model reactant. Paraffinic hydrocarbons, which have C-C and C-H bonds that mimic those in petroleum base oils, are excellent model reactants. Many studies have focused on the oxidation of simple hydrocarbons to elucidate reaction pathways, kinetics, and mechanisms (e.g., Garcia-Ochoa et al., 1989; Lee et al., 1987; Brown and Fish, 1969; Van Sickle, 1972; Van Sickle et al., 1973; Suresh et al., 1988; Mill et al., 1972; Korcek et al., 1972; Jensen et al., 1979,1981;Hombek et al., 1989 W d y et al., 1988; Hazlett et al., 1977). Much of this work, 0888-5885/92/2631-0069$03.00/0
however, was limited to the initial stages of the oxidation, where the hydrocarbon concentration remains essentially constant and hydroperoxides are the major oxidation product. These studies of the initial stages of hydrocarbon oxidation have led to the detailed resolution of reaction fundamentals for the formation of primary (Jensen et al., 1979; Emanuel, 1965; Brown and Fish, 1969; Benson, 1981; Denisov et aL, 1977; Van Sickle et al., 1973; Mill et al., 1972; Boss and Hazlett, 1975; Garcia-Ochoa et al., 1989) and some secondary (Jensen et al., 1981; Brown and Fish, 1969; Boss and Hazlett, 1975; Garcia-Ochoa et al., 1989) oxidation producta. Under more severe reaction conditions (i.e., longer reaction times and higher temperatures), however, the oxidation chemistry becomes more complex, as evidenced by the formation of a broader spectrum of oxidation products (Brown and Fish, 1969). Several researchers (Boss and Hazlett, 1969; Hazlett et al., 1977; Reddy et al., 1988) have investigated the oxidation of hydrocarbons under severe conditions, but much of this work involved only a qualitative or semiquantitative analysis of the reaction products. These studies provided some insight to the reaction pathways, but the lack of a unified protocol for the complete quantitative chemical analysis of hydrocarbon oxidation products may have prevented a more quantitative mathematical treatment of the reaction kinetics. In the first paper of this series (Blaine and Savage, 1991a), we presented an analytical protocol for quantifying the total concentrations of the different oxygen-containing 0 1992 American Chemical Society
