Ind. Eng. Chem. Res. 1994,33, 1126-1135
1126
The Relevance of Thermokinetic Interactions and Numerical Modeling to the Homogeneous Partial Oxidation of Methane Brian F. Gray’ Department of Mathematics and Statistics, The University of Sydney, Sydney, New South Wales, Australia
John F. Griffiths School of Chemistry, University of Leeds, Leeds LS2 9541: England
Gary A. Foulds Department of Molecular Sciences, Division of Chemical Engineering & Industrial Chemistry, James Cook University, Queensland, Townsville, Queensland 4811, Australia
Brian G. Charlton and G. Stewart Walker CSIRO Division of Coal and Energy Technology, PMB 7, Menai, New South Wales 2234, Australia
Both a two-dimensional thermokinetic model and a detailed multidimensional one, only tractable by numerical computation, are applied to the high-pressure, low-temperature, noncatalytic, direct partial oxidation of methane in a continuous stirred tank reactor (CSTR). It is shown that all the phenomena encountered experimentally in this system so far can be accounted for qualitatively by the two-dimensional model and semiquantitatively by the numerical approach. These phenomena include the parametric dependencies of methanol selectivity and yield and the occurrence and significance of hysteresis, bifurcation, oscillation, and extinction, of the direct partial oxidation reaction.
Introduction The direct conversion of methane to methanol by partial oxidation would represent a technological breakthrough, provided that necessary conversion and selectivity criteria are met. To date, the most promising experimental results have been reported for the homogeneous gas-phase conversion of methane to methanol (Yarlagadda et al., 1988; Gesser et al., 1987). However, there is considerable variation in the literature, with experimentally determined methanol selectivities varying from less than 10% (Fukuoka et al., 1989) to over 80% (Yarlagadda et al., 1988; Gesser et al., 1987). In addition, multiple steady states (Foulds et al., 1993) and oscillations (Yarlagadda et al., 1990) have also been observed under certain process conditions. A range in conversion and selectivity is also manifested in the reaction models that have been developed thus far for the homogeneous gas-phase reaction (Danen et ai., 1991; Durante et al., 1989; Onsager et al., 1989; Vedeneev et al., 1988a-d; Vardanyan et al., 1981). This, coupled with the fact that the process economics are very sensitive tomethanol selectivity, initiated the present investigation, from a theoretical point of view, of the nonlinear nature of this highly exothermic reaction system. The aim is to obtain a qualitative and semiquantitative understanding of the thermokinetic interactions, and a secondary objective is to obtain a workable quantitative numerical model. These possibilities arise because a large amount of kinetic data for methane oxidation has been obtained over the last 20 years or so, particularly in the range of “low temperature oxidation”, relevant to the partial oxidation to methanol. Comprehensive rate constant data, involving some 20 odd chemical species, are now available. For example, it has been used by the Russian group (Vedeneev et al., 1988a) to simulate the consumption of 02 and the production of methanol in an 0888-588519412633-1126$04.50/0
isothermal experiment at 614 K and 10 MPa at an 02 content of 12.5 % . The success of this simulation led these authors (Vedeneev et al., 1988b)to believe that the kinetic scheme was useful in a predictive sense, and they used the kinetic schemeto investigate details of the autoacceleration present in this reaction (also found in our own computations). They studied the sensitivity of the reaction behavior for small variations of the rate constants using standard techniques, i.e., calculating the local sensitivity coefficients, which are defined as
the logarithmic form being used in order to make dimensionless. p i j measures the sensitivity of the concentration of species i to variation of reaction rate constant kj at time t. These quantities are useful for investigating numerical difficulties, but they are of little physical significance since it is practically impossible to vary a single rate constant experimentally for a gas-phase reaction. The only control parameters at the disposal of the experimentalist are the equivalence ratio of fuel to 02,the surface to volume ratio, the total pressure, the ambient temperature, and the residence time in a continuous stirred tank reactor (CSTR). We have not therefore pursued the question of local sensitivity coefficients, rather we have directed our attention to the sensitivity, in the large, of the kinetic behavior. The sensitivity considered is with respect to some of the control parameters, in particular, total pressure, equivalence ratio, ambient temperature, and residence time. In a continuation of work in this area, Vedeneev et al. (1988~)consider factors influencing the yields (not the selectivities) of products, in particular, CH30H. Further indication that the full kinetic scheme for methane oxidation has some predictive value comes from another 0 1994 American Chemical Society
Ind. Eng. Chem. Res., Vol. 33, No. 5,1994 1127 paper by Vedeneev et al. (1988d), where ignition delays are calculated for rich methane-air mixtures. Our interest in this work is that it represents the only work of this group which is not isothermal, Le., it includes the effects of exothermicity of the reactions on the temperature of the reacting mixture. This is done by supplementing the kinetics (or mass conservation equations) with the energy conservation equation, assuming Newtonian heat loss from the reacting mixture to the surrounding heat bath. We assume that the reacting mixture is well-stirred, i.e., there are no temperature gradients within it. In this case there is a temperature discontinuity at the vessel wall and the above assumption is that the heat loss is proportional to it, the proportionality constant being the surface to volume ratio of the vessel times a heat-transfer coefficient, x , dependent on the physical nature of the wall. Of course this calculation only has meaning when the system is beyond the self-ignition limit and self-ignition is actually being obtained. In such cases much higher temperatures occur than in the low-temperature regime of interest here, and the kinetics are completely different. Experimental work highly relevant to the present problem is exemplified by Yarlagadda et al. (1990). These workers used a flow tube reactor, and their objective was to find optimum methanol production conditions. However, in and near such regions, they observed temperature oscillations in the reacting gas, with frequencies of the order of 5 min-l. The amplitudes of the oscillations are unknown, since the thermocouple used was far too large to have adequate time response, but they are very probably of the order of those observed in cool flame oscillations, well characterized at lower pressures for many hydrocarbons and their derivatives (Lignola and Reverchon 1987). At lower pressures, the phenomenology of cool flames is reasonably well understood in terms of the chain-thermal or thermokinetic theory of Gray and Yang (1965,1969ad). Use of these models in the low pressure area has given useful insights into the links between oscillating cool flames, negative temperature coefficients of reaction rate, hysteresis phenomena, and bifurcation behavior, where the systemjumps discontinuously from one stateto another on varying a bifurcation parameter, such as ambient temperature. The whole philosophy used in this approach is to identify the system control parameters and then study the variation of the behavior of the system as one or more of these is varied. The benefits of this approach are that modeling is kept in close contact with experimentally feasible procedures and experimental results are presented, and indeed looked for, in a way which is theoretically sound. With systems which are very nonlinear (here largely due to the Arrhenius temperature dependence of the reaction rate constants), this is essential as,for example, the system can be in two different states at the same parameter values, depending on its history. In such circumstances, clear distinction between system variables (concentrations, reaction mixture temperature) and control parameters (ambient temperature, flow rate, etc.) is absolutely essential and has not always been present in previous work, particularly in the high-pressure region. It is expected that if cool flames and associated nonlinear phenomena occur in the high-pressure region, they would be smoothly linked with the low-pressure phenomena and possess similar characteristics. With this assumption we would immediately expect to find a negative temperature coefficient of reaction rate, discontinuities in the concentrations as functions of ambient temperature, and hysteresis phenomena.
Two complementary approaches have been used in the modeling. Firstly, a low-dimensional model, already reasonably well understood, is analyzed qualitatively and semiquantitatively to obtain the general structure of the behavior with respect to the bifurcations that occur in the system. In subsequent computational work on the higher dimensional model these phenomena, when they occur, are easily recognizable and understandable, which would not be the case in a purely numerical approach. Part of this work has been presented at the Symposium of Natural Gas Upgrading, Division of Petroleum Chemistry, Inc., American Chemical Society (Foulds et al., 1992a).
The Two-Dimensional Model This model was the first used to discuss low-temperature oxidation phenomena and was introduced by Gray and Yang (1965, 1969a-d) as a generic model capable of describing the puzzling nonlinear phenomena which had been reported over many years, such as reaction mixture temperature oscillations, visible oscillatory cool flames, hysteresis phenomena, and multistability and its associated jump (bifurcation) phenomena. The formal kinetics in the two-dimensional scheme are as follows:
ki
fuel
-
X
kti
X
product 1
kta
X
product 2
x + ...
kb
2x + ...
X may be regarded as an autocatalytic species or a radical intermediate partaking in a chain branching reaction (hence the notation kb for branching reaction, ktl for termination reaction to product 1,etc., ki then being the initiation reaction). We shall restrict ourselves to the case where the activation energies satisfy the inequality
this producing the kind of behavior shown by the partial oxidation reaction. For the detailed discussion of this restriction see Gray and Yang (1969a). This inequality is the only requirement, that the overall heat release rate exhibit a region of negative temperature coefficient, Le., decrease with increasing reaction mixture temperature. All the other aspects that we wish to describe such as hysteresis, multistability, and thermally coupled oscillations flow from this. It is worth noting a t this point that the richness of behavior shown by a single exothermic reaction in a CSTR is both remarkable and well understood (Golubitaky and Schaeffer 1985). It is therefore not surprising that the more complex reaction scheme considered here can show such complex nonlinear behavior. However, the single exothermic reaction does not show a negative temperature coefficient, and we believe that the two-dimensional scheme discussed here is the simplest form capable of describing all the behavior that we are interested in.
1128 Ind. Eng. Chem. Res., Vol. 33, No.5,1994
The chemical part of this two-dimensional model consists of
ki
fuel
kb
(initiation)
X
h.
(branching)
x - 2 x ku
ha
2 X +stable products
-
(termination)
hu
2X
stable products
(termination)
rrr
Here X is an autocatalytic chain carrier concentration; the symbols k and h are the reaction rate constants and (minus) the enthalpies of reaction, respectively. We assume that the consumption of the original fuel is negligible so that kj is not time dependent. With those assumptions,and that of Newtonian heat loss, we can write the kinetic and energy equations as
@ dt = ki(kfl + k,- k , ) X C p V e = kihi(kf,hfI+ ktJzt2 + k&X dt
(1)
= U T - Td) (2)
in dimensional form. C, p, and V are respectively the mixture heat capacity, density, and volume. L is the reactor surface area multiplied by an empirical heattransfer coefficient, x , dependent on the physical and chemical nature of the wall. Eqs 1and 2 can be written in the dimensionless form as follows dx dt = fJu) + W dt
+
) X
= fz(u) 8(u)x - l(u - u,)
(3) (4)
where u is dimensionless temperature, u, is dimensionless ambient temperature, and x is the dimensionless chain carrier concentration. 1 is the dimensionless form of L. The dimensionless form of would simply be the net branchingfactor,i.e., the coefficientofX in (1). Similarly the dimensionless form of 8 would be the coefficient of X in (2). Equations 3 and 4 are highly nonlinear (due to the Arrhenius terms), although they are linear in x. Thus, x can be eliminated from the steady-state forms of eqs 3 and 4 to obtain the transcendental equation
+
the first two terms representing the chemical heat release rate and the third term representing loss of heat by Newtonian cooling and outflow of hot reacting mixture. Note that 4(u) < 0 in the range where we wish to solve eq 3. is the steady-state temperature. In the combustion literature, this equation is usually treated graphically by plotting the heat release term and the loss term on the same diagram and interpreting the intersections of the two curves as solutions of eq 3. Reference to Figure 1clarifiesthe discussion considerably, since eq 3 may have up to five solutions at a point in parameter space. Local stability analysis shows that the lowest steady state is always stable and that the second and fourth
U.
U
Figure 1. Heat production rate and heat loss ratecurvea YB dimensionlesa ambient temperature. The effect of varying dimensionlesa ambient temperature.
-- -
solutions are always saddle points. Every saddle point and a pair has a pair of curves which finish at it as t of curves which originate at it (and finish at it as t --), and these are known as separatrices. The former pair separate integral curves in the u-u plane which go upwards past the saddle from those which go downwards past the saddle, these two groups finishing up at completely different steady states. This particular pair is called "unstable separatrices" and in this sense act as a "watershed". The third steady state can show variable stability. It can be locally stable or locally unstable, with or without a stable limit cycle surrounding it in the phase (u-x) plane. Of course, the latter condition represents sustained thermokinetic oscillations in the system. We can obtain considerable qualitative insight into the behavior of eqs 1and 2 as the parameter u. is very slowly varied, first upwards, then downwards. Initially, Figure l a describes the system at low u., the lowest stable steady state of the system (i) being populated. On increasing ua, with other parameters fixed, which is a very typical experimentalsituation,atangency between the twocurves occurs and steady states i and ii merge in a saddle-node bifurcation, Figure lb. At this critical ambient temperature, a jump in the reactant temperature occurs, and the system populates steady state iii when that is stable. This intermediate state is accessible on normal starting only between ~ ~ , ~ ~ a n d To u ' anticipate ~ , ~ , . thenumerical results for the CHd02 system, this intermediate steady state represents the one needed to execute direct partial oxidation. At steady state (i) virtually no reaction is taking place, and beyond steady state (iv) is the stable steady state representing full-fledged combustion to largely CO and COz. The above comments are true for increasing un3normal procedure after startup. However, Figure ICreveals that the required direct partial oxidation state (iii) exists over a much wider range of u. when this is being decreased. Extinction does not occur until ua,exis reached, and this is considerably lower
Ind. Eng. Chem. Res., Vol. 33, No. 5, 1994 1129
8
U.
U-
e -b Y
8 r a *
ua
e
UII
U
U,
bsn
U
Figure 2. Heat production rate and heat loss rate curves v8 dimensionless ambient temperature. The effect of varying gas flow rate (residence time).
than u , , ~ ~Of. course, this type of ignition-extinction hysteresis is well known in the combustion literature, but its importance does not seem to have been clearly recognized in this case of “partial ignition” or “cool-flame ignition”. In the next section we shall see that numericaltreatment of a realistic oxidation scheme for CH4 gives first such behavior and sometimes gives higher yield and/or selectivity for methanol in the “upwardly inaccessible” region between Ua,cr and Ua,exReference to Figure Id shows a situation where the middle direct partial oxidation state is completely inaccessible, simply because the heat losses are too small and Ua,er corresponds to full-fledged ignition, with no stable steady state existing in the low-temperature regime. This undesirable condition could arise due to excessive thermal lagging of the reaction vessel, too low a flow rate (long residence time), or too large a reaction vessel (surface to volume ratio too small). Similarly, hysteresis effects occur in this model if the flow rate (or residence time) is varied. This is described very simply by variation of the slope of the straight line at fixed intersection (u,) with the u axis, as shown in Figure 2. Again, the direct partial oxidation state is accessible over a considerably wider range as the slope increases from emaller values, i.e., the flow rate increases. These qualitative parametric dependencies are extremely useful in giving a perspective in what would otherwise be a complicated numerical problem, even in the two-variable case. That they occur in a similar way in a multivariable kinetic scheme, which is a reasonable representation of low temperature methane oxidation, is a great help in interpreting and planning such numerical work and ultimately in the experimental work described elsewhere (Foulds et al., 1993). To conclude this section, it is pertinent to discuss some of the bifurcation diagrams which can occur in the two-
Ua
m
Figure 3. Heat production rate and heat loss rate curves vs dimensionless ambient temperature. The Occurrence and location of bifurcations.
dimensionalmodel, particularly those covering oscillatory reaction states. These have not yet been fully tabulated, and 20 or 30 of them are expected. Some have been obtained (Gray and Gonda, 1983), and these indicate that oscillatory states can arise through Hopf bifurcations and/ or through homoclinic bifurcations. Figure 3 shows some examples. In Figure 3a oscillations commence with finite amplitude on increasing u, at the homoclinic bifurcation (hm) and gradually decrease in amplitude until they vanish at the Hopf bifurcation (h). Similar effects occur in Figure 3b, except the homoclinic bifurcation occurs between the limit cycle and the upper saddle point. Figure 3c shows a bifurcation diagram in a parameter range without oscillations, while Figure 3d shows the case where the lower tangency condition (sn = saddle-node bifurcation) ceases to exist at a hysteresis bifurcation.
The Multidimensional Model The multidimensional CSTR model utilizes a set of ordinary differential equations, comprising the species conservation equations and the energy conservation equation: dxi dt = fi(X1...Xn,U) + kf($ - x i )
du = R(x, ... xn,u) + k f ( U a - u ) - I(u - 24,) dt where
xi
i = 1... n are dimensionless concentrations,
fi
(7)
represent
1130 Ind. Eng. Chem. Res., Vol. 33, No. 5, 1994 Table 1. Kinetic Scheme Showing Enthalpies, Preexponential Factors, and Activation Energies (All in SI Units) for Relevant Gas-Phase Reactions Used in the Modeling -AHg')
ZG)
-110450 115400 -115400 371200 27300 338800 102400 194040
1.00 X IO' 4.00 X 106 1.5 X 10'3 2.50 X 10' 1.00 X 106 1.00 X 106 1.00 X 10' 4.0 X 10' 1.0 x 106 2.0 X 10" 8.5 X 106 5.0 X 10' 5.0 X 10' 6.0X 106 1.00 X 106 1.12 X lo7 2.51 X lo' 2.51 X lo' 1.12 X lo' 2.00 X 106 4.79 X 10' 1.50 X 1 O a 1.00 X 106 2.00 X lo2 1.00 X lo4 3.00 X 106 2.00XlOa 5.00 X 109 1.95 X lo7 6.56 X lo7 2.30X lo7 4.00X lo6 6.30 X 106 2.80 X 106
00000 -184700 62700 63500 -167500 132600 loo00 85770 14640 268620 339750 56900 -120510 -70400 125700 533640 265263 104200 -66900 197300 64000 -8247 70660 8247 72279 502544
EIR
reaction
-AHG)
Zoo00 CH4 + 02 = CH3 + HOz
430266 217587 503026 438483 176800 159707 236170 238651 303576 -210300 122600 60280 15750
CH3 + 02 CHaOz CH3Oz CH3 + 02 CH3 + CH3 = CzHe CH&z + CH302 = CH3O + CHsO + 02 CH3O + CH3O CHsOH + CHzO CH30 + 02 CH20 + HOz CHaOH + 02 = CHzO + HOz CH30 + OH CHzO + HzO CH302H = CH3O + OH CHI + OH = CH3 + Hz0 CHzO + 0 = HCO + OH CHzO + 02 = HCO + HOz CHzO + OH = HCO + HzO CHzO + HOz = HCO + HzOz CzHe + OH = CzHs + HzO CzHe + 0 = CzHs + OH CzHs + 0 = CzH4 + OH CzHs + OH = Cz& + HzO CzH4 + OH = CHzO + CH3 CH3 + 02 = CH3O + 0 HCO + M = H + CO + M HCO + 02 = CO + HOz CO + 0 + M = COz + M CO + HOz = COz + OH CO + OH = COz + H H+Oz=OH+O H + 02 + M = HOz + M OH + Hz = HzO + H 0 + Hz = OH + H OH+O=Oz+H OH+H=Hp+O 550 OH + OH = HzO + 0 0 H + OH + M HzO + M
0 14574 0 0 0 3000 25718 0 20447 1230 2313 19000 0 4022 1340 3200 3200 1233 483 14594 9553 3007 0 5032 302 8298 -654 2850 5337 0 3500
the source (sink) terms for each species arising from the chemistry, u is the dimensionless temperature of the reacting mixture, kf is the feed rate, u4 is the feed and reactor wall temperature, R represents the net thermal effect of the chemistry, 1is again a heat-transfer coefficient incorporating the surface to volume ratio of the reactor and the flow rate, and the superscript zero indicates feed concentration. In this work all but two (CH4and 0 2 ) of the xi0 are zero. For example the species conservation equation for methane is
where k l = 21 exp(-El/RT) etc. and kf is a flow rate. [CH4I0 would be the feed concentration of methane. Similarly for the methoxy radical
k,[CH,OI[OHl + k,[CH,O,HI + k,[CH30,12 + k,[CH,l 1021 - k,[Ml [CH,OI - k1o[CH@l [CHJ k11[CH301[CHZOI - k,2[CH301 [CH,O,I k1,[CH3OI[H021- kh[CH@lo - [CH301] (9) where here [CHsOIO = 0 since methoxy radicals are not present in the feed. In this system we have a n + 1 dimensionless phase space (XI,...,xn,u)and a five-dimensional parameter space, (kf,u4, 1, xio, xz0). We expect the phase space to be essentially two dimensional in so far as the experimental system we are describing has been observed to show only
o
-21700 38000 -43000 0 -505682 -441475 -503026 -1057500 117790 109860 36960 44050 -16240 338600 3000 -66900 -282100 -236100 -110400
-58600
20')
EIR
reaction
1.26 X l(r 0 H+O+M=OH+M 9.10 X lo2 -2650 OH + OH + M = HzOz + M 0 O+O+M=Oz+M 5.00X 102 1.03 X 109 0 H+H+M=Hz+M 1.8OX 106 0 HOz + HOz = HzOz 02 7.00 X 108 950 H O z + H = O H + O H 3.00 X 10s 2000 H O z + H z H z O + O 2.50 X 10' 350 HOz+H=Hz+Oz 6.00 X 106 0 HOz OH HnO + 02 3.10 X 10" 23830 HzOz M OH + OH + M 3.60X 106 0 HzOz + OH a HOz + HzO 1.00X 106 3030 CHzOH + 02 = CHzO HOz 1.00X 106 5190 C&Oz CHzO = CH30zH + HCO 5.00 x 107 10500 CH3O M = CHaO H M 5.00X 106 25000 Cz& + 02 C ~ H S HOz 1.00x106 2500 CzHa + 02 = CzH4 HOz 1.00x 106 9ooo CzHe + HOz CzHa HzO2 8.00X106 3000 CH3 CHzO = C& + CHO 2.2 X 10'0 52920 HzO M OH H M 2.2 X 106 48850 H z + H = H + H + H 1.2 X 106 54250 0 2+ M = 0 0 M 1.0 X lo8 4oo00 C O z + M = C O + O + M 5.20 X 106 2593 CH3OH H = CH3 + HzO 4.00x106 1010 CH3OH + OH = CHzOH + Hz0 1.70X 106 1157 CHaOH 0 = CHzOH + OH 1.8OX 106 4951 CH3OH CHa = CHaOH CH4 6.30 X 106 9780 CHsOH HOz CHzOH + Hz02 1.00X 106 0 CHzOH + CHzOH = CHzO + CHsOH 2 . 0 0 ~105 3550 CHsO CH( CHaOH CH3 4.0Xl0' 1000 CH5O CHzO CHsOH HCO 4.00 X 106 0 CHaO + CH30z = CHBOOH+ CH3O 1.0 X loe 0 CH30 + HzOz CHsOH 02 2.0 X lo' 541 CH3 + HOz = CHsO OH 1.ox 106 1126 CH3Oz C& = CH302H + CH3
+ +
+
+ +
+ + + +
+
+
+ +
+
+ +
+ + + +
+
+
+
+
+
+
+
+
two dimensional behavior, i.e., periodic solutions in addition to an asymptotic approach to singularities which are attractors and represent physically stable steady states of the system satisfying
R(xl ... xn,u) + k&u4-E) - l(u - u4) = 0 ji(X,
...x n , g + k,(x? - X i ) = 0
(10) (11)
where the underscore indicates a steady-state value. The chemical kinetic input is shown in Table 1. The data have been derived from various sources in the literature (Baulch et al., 1988;Griffiths and Sykes, 1989; Griffiths and Phillips, 1990). An average heat capacity and density are assumed for the mixture in the computation of time dependent behavior. In this work Twd,the wall or ambient temperature, has been chosen as the preferred parameter, and flow rate (residence time) and surface to volume ratio have been chosen as secondary parameters, along with 02 concentration and total pressure. Typical calculated results are shown in Figure 4, which represents steady state situations at 9.5% 0 2 , 5.0 MPa total pressure, and 5 s residence time. The arrows represent the direction of variation of Twd on the solid curve. Thus, at approximately670K,an upward transition occurs (corresponding to the lower tangency in Figure lb) and significant self-heating and methanol production commences. The upper branch of the solid curve represents the degree of self heating (AT = T - Twd)in the direct partial oxidation state as Twdis decreased. Here, state (iii)is decreasing toward ita extinction value as shown in Figure IC. This occurs just below 590 K at the highest methanol concentration obtainable under these parametric conditions (3.4% of carbon products). Note that this
Ind. Eng. Chem. Res., Vol. 33, No. 5, 1994 1131 300
160
I
T
I
P
- 8 0 -
l-
4
80-
40-
20-
T wall (K)
Twall (K)
Figure 4. Simulation of AT and methanol concentration v8 T d , using 9.5% 0 2 in the feed, 5.0 MPa pressure, and a residence time of 5 8.
Figure 6. Simulation of AT and methanol concentration vs T d , using 9.5% 0 2 in the feed, 5.0 MPa pressure, and a residence time of 20 8.
I_
T
200-
-
g'"
k Q imW -
-------
,
____ --__
20
i
T / I
11
2
T
, "..".. , .....,...__ CHgH ..__ , '..__ ...__
P
I.. II
+
-...___
lo
$u
NO HYSTERESIS
F
TMII (K) Figure 6. Simulation of AT and methanol concentration vs Twd, using 9.5% 02 in the feed, 5.0 MPa pressure, and a residence time of 10 8.
methanol production is not attainable at all on upwards only sweeps of TWd. Figure 5 shows a similar set of runs for a 10-5residence time, and the interpretation in terms of Figure 1 is once more evident. Again extinction on the downward sweep occurs at the maximum methanol concentration (13.8% of carbon products). Figure 6,at a 20-5residence time, begins to exhibit some quantitatively different features, the hysteresis loop having decreased from -80 K for a 5-5residence time to only -10 K in this case. This variation is easily understood by reference to Figure 2. Also in this case, there is a local maximum in methanol concentration, corresponding to 25.2% of the carbon products, before the extinction point at 640K. This indicates that a local maximum in methanol concentration exists beyond the extinction points in the earlier cases, if such areas were accessible. Figure 7 completes this particular set of results with a 30-5residence time. In this case, the straight line in Figure 2 is now so steep that there is no hysteresis in the lowtemperature range, but a region of sustained oscillation has appeared with a methanol maximum of 26.2% of carbon products, a t a temperature just above the extinction of oscillations at a bifurcation. The variation of calculated methanol selectivity and methane conversion, under these conditions, with varying residence time, is shown in Figure 8. Initially, methanol selectivity is found to increase with increasing residence time, before leveling off at about 30%. Methane con-
Figure 7. Simulation of AT and methanol concentration vs T d , using 9.5%0 2 in the feed, 5.0 MPa pressure, and a residence time of 30 a. 4ol
01
I
1
10
I
20
30
40
! 50
Residence Time (secs)
Figure 8. Simulation of methane conversion at maximummethanol selectivity vs residence time, using 9.5%02 in the feed and 5.0 MPa pressure.
version is found to decrease rather weakly as the residence time is increased, leveling off at 6% conversion. In Figure 9 the variation of methanol selectivity with wall temperature is shown for various 02 concentrations. In agreement with experimental results (Foulds et al., 1992b,1993;Ford 1992;Rytz andBaiker, 1990;Yarlagadda et al., 1988), the selectivity to methanol increases with decreasing 02 concentration, with a local maximum at a wall temperature of 445 "C of 50% for 2.5% 02. Higher
1132 Ind. Eng. Chem. Res., Vol. 33, No. 5, 1994
Table 2. Simulation of the Effect of the Overall Heat-Transfer Coefficient on Hysteresis, Conversion, and Selectivity at Maximum Methanol Yield, Using 9.5% Oa in the Feed, a Pressure of 3.0 MPa, and a Residence Time of 10 s L (W rn4 K-l) 5500
CHI Twd (K) T (K) AT (K) conv (%) 148 6.3 674 822
6.3
5000
667
823
156
4500
658
826
168
6.5
4000
648
828
180
6.4
3500
640
830
190
6.6
3000
623
828.5
205.5
6.5
mol m4
% total C in products
molm9 % total C in products
molm" % total C in products molm4 % total C in products molm4 % total C in products rnolm" % total C in products 1
Bo,
380
1
1
1
400
420
,
I
4 4 0 4 8 0
480
5Ml
Twall ("c) Figure 9. Simulation of methanol selectivityvs 2'411 using an overall heat-transfer coefficient of L = 4500 W m4 K-I, a pressure of 3.0 MPa, and feed oxygen concentrations of 2.5,5.0, 7.5,and 9.5%.
' 3uI
CHI CHaOH CH2O CO COz C& 3.6 0.16 19.7 7.0 0.61 483.3 62.5 22.5 1.9 11.4 0.5 3.5 0.15 19.8 7.1 0.70 483.0 11.0 0.5 62.4 22.4 2.21 3.2 0.13 20.1 7.2 0.71 482.4 63.0 22.7 2.22 10.1 0.4 482.6 3.3 0.14 20.0 7.2 0.75 10.3 0.4 62.6 22.5 2.35 481.6 3.0 0.12 20.1 7.45 1.03 9.1 0.36 62.7 22.5 3.11 3.3 0.14 19.9 7.24 0.81 482.3 10.3 0.44 62.1 22.6 2.52
380
s
sm
380
Twall (C) Figure 10. Simulation of the effect of the heat-transfer coefficient, L, on the position and magnitude of hysteresis using 9.5% 0 2 in the feed, 3.0 MPa pressure, and a residence time of 20 s (A) L = 5000 W m4 K-I; (B) L = 6OOO W m-S K-1; (C) L = 6500 W m-S K-l; (I)) L = 7500 W m-9 K-1).
selectivities could be obtained by allowing the 02 concentration to decrease further. Variation of the heat transfer coefficient (Le., the slope of the straight line in Figure 2) clearly affects the existence and size of the hysteresis effect. A particular example is shown in Figure 10 where clearly the hysteresis loop decreases as the heat-transfer coefficient is increased, until it finally disappears. Further data on the chemical effects of varying the heattransfer coefficient are listed in Table 2. The selectivity for methanol, at its optimum, does not seem to be a sensitivefunction of the heat-transfer coefficient. Clearly, this is because the reactant temperature at maximum
OZ 3.00 2.91 2.61 2.56 2.55 2.58
Table 3. Simulation of the Effect of Residence Time and Overall Heat-Transfer Coefficient on the Hysteresis Effect, Using 9.5% 0 2 in the Feed and a Pressure of 3.0 MPa L (Wm-S K-9 Twa~,rxt('C) T w a ~ i g("C) loop ("C) T~ (8) 393.6 32.6 10 3500 361.0 371.0 394.2 23.2 4000 395.0 4500 380.0 15.0 388.0 395.8 7.8 5000 391.8 5500 394.0 3.8 376.0 381.0 5.0 20 3000 390.0 390.0 0.0 4000 399.6 6000 399.6 0.0 Table 4. Simulation of the Effect of Increasing the Preexponential Factor (Z) of the Initiation Reaction on Hysteresis and Methanol Yield. All Determined Using L = 6000 W m-s K-1, a 20-8 Residence Time, 9.5% 0, in the Feed, and a Pressure of 5.0 MPa
1 x 107
/ I
C2H4 0.40 1.30 1.48 1.51 0.50 1.57 0.56 1.75 0.077 2.32 0.64 2.00
373
369
12.8 (492,373)' 2 x 107 360 353 15.9 (475,357)' 3 x 107 350 338 19.0 (470,349)' 4 x 107 341 334 20.6 (460,341)' 5 x 107 339 328 21.9 (456,337)' 1 x 108 324 313 26.9 (452,332)' Gas temperature and wall temperaturerespectively at maximum CH30H yield.
methanol yield is very nearly independent of the heattransfer coefficient, although the wall temperature is very different. Of course, the data in Table 2 are produced by scanning !fwdfor the methanol optimum at each given heat-transfer coefficient, and it is not surprising that the answer occurs at the same reactant temperature, as we are observing here what is essentially a reactor-independent property. Clearly, in a flow reactor which is nonadiabatic, the straight line loss curve in Figure 2 has two components, one inversely proportional to the heat transfer coefficient and the second inversely proportional to the residence time. Both of these components will affect the presence and nature of hysteresis. The results of varying both are shown in Table 3. As expected, the hysteresis loop decreases as the heat transfer coefficient is increased, as it does when the residence time is increased. Individual variation of rate constants as parameters has not been carried out in this study as, in general, there is no physical parallel for such variation. Neither has a sensitivity analysis of the system of differential equations
Ind. Eng. Chem. Res., Vol. 33, No. 5, 1994 1133 with respect to individual rate constants been carried out, although this has been done by other workers (Mackie, 1991). The effect of increasing the preexponential factor of the initiation reaction has also been simulated, as this is a possible physical variation. This was done using 9.5% 02 in the feed a t a total pressure of 5.0 MPa. The results are listed in Table 4. As the preexponential factor increases, the hysteresis range increases and the maximum methanol yield doubles for a 10-fold increase in rate constant. The temperature a t which maximum yield occurs also drops considerably. Simulations were also carried out at 3.0 MPa, and similar comments apply. As far as the oscillatory regions of reaction are concerned, in these regions the methanol concentration varies over many orders of magnitude during each oscillation. The period of these oscillations is of the order of a minute so that, in principle, some sort of methanol enrichment by means of a microprocessor-operated control system is not out of the question. In Figure 11an example, using a wall temperature of 652 K, 9.5% 02 in the feed, 5.0 MPa pressure, and a residence time of 30 s, is shown. The period of oscillation is about 90 s, and the temperature amplitude is about 140 "C. The methanol oscillation is biperiodic for reasons not understood at the moment, and its concentration varies over the range 1.5-9.5 mol ma during each oscillation, the minimum corresponding with the peak in the singly periodic 02 oscillation. In Figure 12the 02oscillations are shown using the same conditions. Interestingly the 02 concentration is almost completely depleted during each cycle. The corresponding CO concentrations are depressingly high (from the point of view of selectivity for methanol) even at its minimum value of 18 mol m3. However, there are instants during the oscillation when selectivity for methanol is very much better than the average over one period, which is usually similar to the steady-state value near the oscillatory region. In Figure 13 the selectivity for methanol is shown as a function of time over one oscillation period under real
16
652K
650K
400
600
654K
1000
800
1200
656K
1400
1800
1600
TIME (sec)
Figure 11. Simulation of methanol concentration vs time in an oscillatory region using 9.5% 02 in the feed, 5.0 MPa pressure, and a residence time of 30 8. TwaiR being varied from 650 to 656 K. 6of
AVERAGE VALUE = 16.0%
10-
k! W
*
5-
0~
I
I
I
I
I
I
I
I
I
I
1134 Ind. Eng. Chem. Res., Vol. 33, No. 5, 1994 conditions. It is well over 20% for a significant fraction of the cycle compared with the steady-statevalue of around
16%.
There is clearly much scope for investigating the phase relationships during oscillation between the desired product, in this case methanol, and byproducts, with the possibility of significant increases in selectivity. Conclusions
Detailed numerical modeling of the direct partial oxidation of methane under fuel-rich conditions at high pressures, including heat transfer effects, gives results which are in semiquantitative agreement with experimental CSTR results. The results are also in qualitative agreementwith asingle two-dimensionalmodel which gives a very convenientmode of interpretation of otherwise complicated nonlinear phenomena. Observed phenomena such as hysteresis, bifurcation, sustained oscillation, and homoclinic extinction are all accounted for semiquantitatively b y the numerical model and qualitatively by the two-dimensional model. Parametric dependencies of methanol selectivity and methane conversion on wall temperature, residence time, % 02,total pressure, and heat-transfer characteristicsare also accounted for in semiquantitative fashion by the numerical approach. Observations of higher methanol yields in downward sweeps of wall temperature, comparedwith upward sweeps, are accounted for by the model. Acknowledgment
The authors wish to thank the Australian Research Council for the award of an Australian Senior Research Fellowship (B.F.G.) and B.H.P Co. Ltd. for financial assistance. Nomenclature Bij = local sensitivity coefficient = dimensionless concentration of component i kj = reaction rate constant t = time ki = rate constant for initiation reaction kb = rate constant for branching reaction ktz = rate constant for termination reaction for product 1 kt2 = rate constant for termination reacton for product 2 hi = negative enthalpy of reaction for initiation reaction hb = negative enthalpy of reaction for branching reaction htl = negative enthalpy of reaction for termination reaction for product 1 ht2 = negative enthalpy of reacton for termination reacton for product 2 AH = heat of reaction (J mol-') u = dimensionless temperature 5 = dimensionless steady-state temperature u4 = dimensionless ambient or wall temperature u . , ~=~critical dimensionless ambient temperature u ~= extinction , ~ ~ dimensionless ambient temperature 1 = dimensionless heat-transfer coefficient incorporating the surface to volume ratio of the reactor hm = homoclinic bifurcation h = Hopf bifurcation sn = stable node k, = feed rate to CSTR R = net thermal effect of chemistry the source (sink) terms for each species arising from the fi chemistry TWd= wall or ambient temperature (K or "C) xi
T = reaction mixture temperature (Kor "C) Twd,ext = wall temperature at extinction (K or "C) Twdd= wall temperature at ignition (Kor "C) T = residence time ( 8 ) Zi = preexponential factor (mol-' m3 8-1) x = empirical heat-transfer coefficient dependent on physical characteristics of reactor vessel (Wm-2 K-1) X = autocatalytic chain carrier concentration (mol dma) S = reactor surface area (m2) V = volume of reactor (m3) L = x (S/V): heat transfer-coefficient incorporating the surface to volume ratio of the reactor (Wm a K-1) C = heat capacity of gas mixture (J kg-1 K-1) p = density of gas mixture (kg m-3) Ei = activation energy (J mol-' K-1) R = universal gas constant
Literature Cited Baulch, D. L.; Griffiths, J. F.; Pappin, A. J.; Sykes, A. F. Stationary State and Oscillatory Combustion of Hydrogen in a CSTR. Combust. and Flame 1988, 73, 163-185. Danen, W. C.; Ferris, M. J.; Lyman, J. L.; Oldenborg, R. C.; Rofer, C. K.; Steit, G. E. Methane to Methanol Conversion by Direct Partial Oxidation. Preprints, Symposiumon Methane Upgrading, 201st National Meeting of the AmericanChemicalSociety,Division of Petroleum Chemistry, Atlanta, GA, April 14-19,1991; American Chemical Society: Washington, DC, 1991; pp 166-171. Durante, V. A.; Walker, D. W.; Seitzer, W. H.; Lyons, J. E. Vapor Phase Hydroxylation of Methane. Preprints of 3B Symposium on Methane Activation, Conversion and Utilization, 1989 International Chemical Congress of Pacific Basin Societies, Honolulu, HI, 1989; V0148, pp 23-26. Ford, M. J. The Noncatalytic Partial Oxidation of Natural Gas to Methanol. Preprints, Symposium on Natural Gas Upgrading 11, 203rd National Meeting of the American Chemical Society, Division of Petroleum Chemistry,San Francisco, CA, April 5-10, 1992;American Chemical Society: Washington, DC, 1992;pp 3440. Foulds, G. A.; Gray, B. F.; Griffiths, F.; Walker, G. S. Theoretical Modelling of the Gas Phase Partial Oxidation of Methane. Preprints, Symposium on Natural Gas Upgrading 11, 203rd National Meeting of the American Chemical Society, Division of Petroleum Chemistry, San Francisco, CA, April 5-10, 1992a; American Chemical Society: Washington, DC, 1992; pp 51-60. Foulds,G. A.; Miller, S. A.; Walker G. S. Gas Phase Partial Oxidation of Methane. Preprints, Symposiumon Natural Gas Upgrading11, American Chemical Society, Division of Petroleum Chemistry, San Francisco Meeting, April 5-10, 1992b; pp 26-33. Foulds, G. A.; Gray, B. F.; Miller, S. A.; Walker, G. S. Homogeneous Gas-PhaseOxidation of Methane Using Oxygen as Oxidant in an Annular Reactor. Ind. Eng. Chem. Res. 1993,32, 780-787. Fukuoka, N.; Omata, K.; Fujimoto,K. Effect of Additives on Partial Oxidation of Methane. Preprints of 3B Symposiumon Methane Activation, Conversion and Utilization, 1989InternationalChemical Congress of Pacific Basin Societies, Honolulu, HI, 1989; Vol. 135, pp 106-107. Gesser, H. D.; Hunter, N. R.; Morton, L. A.; Yarlagadda,P. S.; Fung, D. P. C. The Direct Conversion of Methane to Methanolby a High PressurePartial Oxidation Reaction. Prepr.Pap.-Am.Chem.SOC., Diu. Fuel Chem. 1987,32 (3), 255-259. Golubitsky, M.; Schaeffer, D. G. Singularities and Groups in Bifurcation Theory VI; Springer-Verlag: New York, 1985. Gray, B. F.; Yang, C. H. On the Unification of the Thermal and Chain Theories of Explosion Limits. J. Phys. Chem. 1965, 69, 2747-2750. Gray, B. F.; Yang, C. H. Unified Theory of Explosions, Cool Flames and Two-Stage Ignition, Part 1. Trans. Faraday SOC.1969a, 65, 1603-1613. Gray, B. F.; Yang, C. H. Unified Theory of Explosions, Cool Flames and Two-Stage Ignition, Part 2. Trans. Faraday Soc. 1969b,65, 1614-1622. Gray, B. F.; Yang, C. H. Unified Theory of Explosions, Cool Flames and Two-Stage Ignition, Part 3. Trans. Faraday SOC.1969c, 65, 2133-2141. Gray, B. F.; Yang, C. H. On the Slow Oxidation of Hydrocarbonsand Cool Flames. J. Phye. Chem. 1969d, 73,3395-3406.
Ind. Eng. Chem. Res., Vol. 33, No.5, 1994 1135 Gray, B. F.; Gonda, I. The Unified Thermal and Chain Branching Model of Hydrocarbon Oxidation in a C.S.T.R. Roc. R. SOC.1983, A398,133-152. Griffitha, J. F.; Sykes, A. F. Numerical Studies of a Thermokinetic Model of Ethanal Oxidation in a C.S.T.R. R o c . R. SOC.London 1989, A422,289-310. Griffitha, J. F.; Phillips, C. H. Combustion of Ditertiary Butyl Peroxide. Combust. and Flame 1990,81,304-316. Lignola, P. G.; Reverchon, E. Cool Flames. Bog. Energy Combust. Sci. 1987, 13, 75-96. Mackie, J. C. Partial Oxidation of Methane: The Role of the Gas Phase Reactions. Catal. Rev. Sci. Eng. 1991 33 (1,2), 169-240. Onsager, 0. T.; Soraker, P.; Lodeng, R. Experimental Investigation and Computer Simulation of the Homogeneous Gas Phase Oxidation of Methane to Methanol. Preprinta of 3B Symposium on Methane Activation, Conversion and Utilization, 1989 International Chemical Congress of Pacific Basin Societies, Honolulu, HI, 1989; Vol. 135, pp 113-114. Rytz, D. W.; Baiker, A. Partial Oxidation of Methane to Methanol in a Flow Reactor a t Elevated Pressure. Znd. Eng. Chem. Res. 1991,30,2287-2292. Vardanyan, I. A.; Yan, S.; Nalbandyan, A. B. Mechanism of the Thermal Oxidation of Methane 1. Mathematical Modeling of the Reaction. Kinet. Katal. 1981,22,845-851. Vedeneev,V. I.; Gol'denberg, M. Ya.; Gorban', N. I.; Teitel'boim, M. A. Quantitative Model of the Oxidation of Methane a t High Pressures. I. Description of Model. Kinet. Katal. 1988a, 29,7-13.
Vedeneev, V. I.; Gol'denberg, M. Ya.; Gorban', N. I.; Teitel'boim, M. A. Quantitative Model of the Oxidation of Methane at High Pressures. 11.Mechanismof Autoacceleration.Kinet. Katal. 198813, 29,14-20. Vedeneev,V. I.;Gol'denberg, M. Ya.; Gorban',N. I.;Teitel'boim,M. A. Quantitative Model of the Oxidation of Methane at High Pressures. 111. Mechanism of Formation of Reaction Products. Kinet. Katal. 1988c, 29,1291-1296. Vedeneev, V. I.; Gol'denberg, M. Ya.; Gorban', N. I.; Karnaukh, A. A.; Teitel'boim, M.A. Quantitative Model of the Oxidation of Methane at High Pressures. IV. Ignition Delays. Kinet. Katal. 1988d, 29, 1297-1304. Yarlagadda, P. S.; Morton, L. A.; Hunter, N. R.;Gesser, H. D. Direct Conversion of Methane to Methanol in a Flow Reactor. Ind. Eng. Chem. Res. 1988,27,252-256. Yarlagadda, P. S.; Morton, L. A.; Hunter, N. R.; Gesser, H. D. Temperature Oscillations During the High Pressure Partial Oxidation of Methane in a Tubular Flow Reactor. Combust.Flame 1990, 79,216218.
Receiued for review September 2, 1993 Revised manuscript received January 3, 1994 Accepted January 7, 1994. Abstract published in Advance ACS Abstracts, March 1, 1994.
