1424
Ind. Eng. Chem. Res. 1987, 26, 1424-1434
Corti% Corberin, V.; Corma, A.; KremeniE, G. Ind. Eng. Chem. Prod. Res. Deu. 1985, 24, 62. Crozart, M.; Germain, J. E. Bull. SOC.Chim. Fr. 1973, 1973, 2498. Cullis, C. F.; Hucknall, D. J. Catalysis; Bond, G. C., Webb, G., Eds.; The Royal Society of Chemistry: London, 1982; Vol. 5, Chapter 7. Exner, 0. Collect. Czech. Chem. Commun. 1966, 31, 3222. Fisher, T. Statistical Methods for Chemists; Wiley: New York, 1955. Hinshelwood, C. N. The Kinetics of Chemical Change; Oxford University Press: New York, 1940; p 207. Hucknall, D. L. Selective Oxidation of Hydrocarbons; Academic: London, 1974. Isaev, 0. V.; Margolis, L. Y. Kinet. Katal. 1960, 1, 237. Jaswal, I. S.; Mann, R. F.; Junsola, J. A.; Downie, Y. Can. J. Chem. Eng. 1969, 47, 284. Keulks, G. W.; Yu, Z.; Krenzke, L. D. J . Catal. 1983, 84, 38. Keulks, G. W.; Rosynek, M. P.; Daniel, Ch. Ind. Eng. Chem. Prod. Res. Deu. 1971, 10, 138. Khiteeva, V. M.; Rzakulieva, Sh. M. Russ J . Phys. Chem. 1981,55, 1202. Klissurski, D. G. Proceedings of the Climax Third International Conference on the Chemistry and Uses of Molybdenum; Barry,
H. F., Mitchell, P. C. H., Eds.; Climax Molybdenum Company: Ann Arbor, MI, 1979; p 123. Krenzke, L. D.; Keulks, G . W. J . Catal. 1980, 64, 295. Mars, P.; van Krevelen, D. W. Chem. Eng. Sci. (Spec. Suppl.) 1954, 3, 41. Matsuura, I. Proceedings of the 6th International Congress on Catalysis; Bond, G. C., Wells, P. B., Tompkins, F. C., Eds.; The Chemical Society: London, 1977; p 819; Vol. 2. Minachev, Kh. M.; Antoshin, G. V.; Guin, N. K.; Klissurski, D. G.; Abadzhijeva, N. Ts. Chem. Phys. Lett. 1977,48, 515. Monnier, J. R.; Keulks, G. W. J . Catal. 1981, 68, 51. Ono, T.; Kubokawa, Y. Bull. Chem. Soc. Jpn. 1982, 55, 1748. Takasu, Y.; Matsui, M.; Tamura, H.; Kawamura, S.; Matsuda, Y.; Toyoshima, I. J . Catal. 1981, 69, 51. Tasc6n, J. M. D.; Cort6s CorberBn, V.; KremeniC, G.; Gonzglez Tejuca, L. J . Colloid Interf. Sei. 1985, 106, 269. Uda, T.; Lin, T. T.; Keulks, G. W. J . Catal. 1980, 62, 26. Viswanath, V. J.; Chanda, M. J . Appl. Chem. Biotechnol. 1978,28, 1.
Received for review September 19, 1986 Accepted March 6, 1987
Modeling of Exothermic Solid-Solid Noncatalytic Reactions Jan Puszynski, Jan Degreve, and Vladimir Hlavacek* Department of Chemical Engineering, S t a t e University of N e w York at B u f f a l o , Buffalo, New York 14260
A certain class of refractory materials can be synthesized in the sintering or in the self-propagating regime by a direct reaction between the metal powder and carbon, boron, or silicon. A detailed analysis and modeling of such processes are given. Criteria are derived for stable, unstable, and degenerated combustion as well as approximate relations for the velocity of a constant-pattern profile, as a function of the heat-transfer coefficient. The approximate formulas are verified by a detailed numerical study. Different types of combustion (constant-pattern profiles, standing waves, and complex oscillatory behavior) have been observed. There exists an important class of heterogeneous noncatalytic reactions of the type solid-gas and solid-solid which are accompanied by a substantial generation of heat. The most important are oxidation of metals and nonmetallic compounds, direct reaction of nitrogen or hydrogen with certain metals, synthesis of carbides, borides, sulfides, or silicides of metals from pure components, reduction of metallic oxides by aluminum, synthesis of intermetallic compounds (e.g., Ni-AI), combustion of solid propellants, and many others. In several of the examples mentioned above, the process takes place as follows: heat of reaction or externally transferred heat evaporates the material, the reaction occurs in the gas phase, and the product is obtained as a solid or a liquid. In other cases, the exothermic chemical reaction results in an excessive gas-phase formation (oxidation of propellants). In this paper, we will address reactions occurring without gasphase formation. Such reactions are sometimes referred to as gasless combustion reactions. Many of these reactions are of major technological importance. Reduction or roasting of ores, thermite reactions in ordnance applications, and synthesis of new refractory materials from transition or rare earth metals and boron or carbon are a few examples. The exothermic noncatalytic reaction can proceed through the reacting mixture in a uniform way, and the concentration of the components at any instant is almost the same everywhere (the kinetic or “sintering” regime). However, if the system was locally ignited by external 0S8S-5885/87/2626-1424$01.50/0
means, a steep reaction front can propagate through the mixture. Evidently, the latter phenomenon is closely related to flame propagation in gaseous mixtures and is frequently referred to as the self-propagating regime. Synthesis of borides and carbides of transition and rare earth metals can be performed both in sintering and self-propagating regimes. The synthesis in the sintering regime requires a substantial energy input, a high reaction temperature, and a long reaction time. On the other hand, the self-propagating synthesis calls only for a short-range initial energy input. The reaction is self-sustaining owing to the high values of the activation energy and heat of reaction (Merzhanov, 1974). Recent experimental observations by Novikov et al. (1974) and Maksimov et al. (1981) and asymptotic analysis by Margolis (1983) reveal that there may exist several types of propagating waves (constant-pattern profiles, oscillatory waves, and spinning and fingering fronts). However, in this paper we will analyze only the one-dimensionalsituation. Breaking of symmetry in the transverse direction will be the topic of another paper. There is extensive combustion-oriented literature on propagation of fronts in reacting systems. Many of the results can be applied to the theory of noncatalytic exothermic reactions in solid-solid and solid-gas systems. In addition, there is a recent extensive treatment of combustion systems based on activation energy asymptotics. In the present paper, we make an attempt to amalgamate the dispersed information on propagation phenomena in 0 1987 American Chemical Society
Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987 1425 exothermic combustion systems into a compact theory supported by detailed numerical calculations. This work represents a starting point for the description of more complex effects including breaking of symmetry, melting and solidification, diffusion and gas filtration, and free convection phenomena.
Model Formulation In this chapter model equations will be formulated which describe exothermic solid-solid noncatalytic heterogeneous reaction systems. An assumption is made that a heterogeneous reaction occurs between two solid materials, S and M, S(s) + M(s) P(s)
-
7
> 0:
u=O 7
= 0:
8=Bi
O < U < ~
at
1 aT d2T -d2T + + + dr2 r ar a9
+
r2 a& M - W C , exp(-k/RT
(2)
-ac,- - -koC, exp(-E/RT)
(3) at We will assume, in our further analysis, that there is no symmetry breaking. This assumption reduces eq 2 to a one-dimensional case.
The symmetry-breaking situations will be analyzed elsewhere (Puszynski and Hlavacek, 1987). The boundary conditions can be formulated in the form aT t > 0: z = O T=Ti 2-m -=0(5) az
The initial conditions are
t = 0: O l z < m T=To (6) We can easily render these equations dimensionless.
8 ~ 8 0
(9) (10)
.[ 21
112
u =
-
au
Here we have denoted
and that the reaction rate can be described by integer power law kinetics: To develop the model, the following physical assumptions have been made: (1)A heterogeneous mixture of solid powders, S and M, behaves as an isotropic homogeneous system. (2) The temperature dependence of the reaction rate constant can be expressed in the Arrhenius form. (3) Heat conduction in the solid phase can be described in terms of the Fourier law. (4) Mass diffusion of the solid reactants or products does not occur. (5) All physical properties (density, heat capacity, and effective thermal conductivity) are assumed constant. (6) The reaction is not accompanied by melting effects. (7) Radiation effects inside the porous layer do not play an important role. Heat loss from the system can be described by an effective linear driving force. (8) The reactant, M, is in excess so that the reaction process can be considered of the first order with respect to S. Gas-solid noncatalytic reactions of the type S(s) + A(g) P(s) can be handled in a similar way if diffusion effects are unimportant. This is satisfied for high-pressure systems and thin layers of S. The exothermic reaction is supposed to occur in a porous cylinder initially composed of the powders of S and M. At time t = 0, the initial temperature Tois suddenly increased at the end of the cylinder to Ti. This temperature is high enough to ignite the mixture. Heat and mass balances in the porous cylindrical body can be described by
a0 -=0
U-m
( p C,) R T,2 exp(E/ R T,)
t* =
Eko(-AmCso R T,
PC,,RT*~
Y = C,&(-AH)
O=
E
For the nonadiabatic case, eq 7 reads
Here the term 4 8 - 8 0 ) represents the heat loss from the surface and 4KHt* a=-
dPC,
Numerical Solution of the Model Equations The governing equations describing propagation of a reaction front in the solid material are represented by a nonlinear parabolic partial differential equation (eq 11) and an ordinary differential equation (eq 8). In our particular problem, because of high values of both heat of reaction and activation energy, very steep reaction fronts may occur. Within the reaction front, the solution varies rapidly over a very thin region. The "boundary layers" on concentration and temperature profiles represent a very difficult numerical problem. In addition, for the oscillating reaction fronts, the governing equations become stiff. To calculate the reaction-front propagation, we have used two types of finite difference procedures: (a) highorder compact finite difference scheme; (b) transformation of the space interval with reaction front velocity iteration. Both methods will be briefly described below. High-Order Compact Finite Difference Scheme This approach makes use of the well-known StormerNumerov formula for solving boundary value problems of the type Y" = f ( x , y )
(12)
with boundary conditions x=xo x=x1
y=yo y=y1
(13)
The Stormer-Numerov approximation reads (Ciment et al., 1978)
1426
Ind. Eng. Chem. Res., Vol. 26, No. 7 , 1987
Yn-1 - 2Yn
+ Yn+l
=
h2
+ lOf(xnt y n ) + f ( X n + I , . ~ n + l ) l + 10(h4) (14)
zV(xn-19
Yn-1)
A quasi-linear parabolic partial differential equation of the type given by eq 11 can be rewritten in the form
x = tanh ( a [ )
The reduced coordinate x is now x E (-1, 1). Equations 18 and 19 now read a0 = [a(l - x 2 ) ] 2
ae a2e + a(1 - x2)(2a+ + V,) + ax
ar2
a7
Making use of approximation 14, we can write
with boundary conditions x=-1 0=0, x = l
Evidently, the procedure retains the tridiagonal structure of the problem. This scheme requires almost the same computational effort as the standard Crank-Nicolson approach; however, the accuracy is O(h4,k 2 ) . The nonlinear finite-difference equations were solved by quasi-linearization and iteration. An automatic timestep adjustment has been adapted in the following way: (i) Extrapolate (from the previous two time profiles, j and j - 1) to estimate the first approximation of the dependent variables at the new time-level, j + l. The extrapolated values are used for calculation of the nonlinear source term and its derivatives in the quasilinearized scheme. (ii) Solve a set of linear finite-difference equations. (iii) Calculate eyi = yp+’lP1 - yplj+l. If for each leYil < e*, then the particular time profile is considered to be resolved. If for any leril > e* and iteration counter p greater than the maximum number of iterations p., then the present time step, k , is reduced, 1212 k , new extrapolation is performed, and the calculation starts again. After five successfully resolved time profiles are found during which the time step was not reduced, the time step is doubled, 2k k . Here e , = 1 X lo-* and p. = 6. The number of grid points was changed from 200 to 1200 according to the complexity of the combustion process in question.
-
-
Transformation of the Space Interval with Reaction Front Velocity Iteration When the reaction process occurs in an infinite system, it is appropriate to use a new system of coordinates which is moving with the reaction front: t = u - VET (17)
0 = 1
(22)
8=00
Equations 20 and 21 with boundary conditions 22 have been solved by using central and upwind approximation of second and first derivatives, respectively, and fully implicit time integration. In order to calculate a new profile, it is convenient to fix one point on the temperature or conversion profile (e.g., we have used the point for which the conversion is equal to 0.5) as the origin of the coordinate system. In such a case, the velocity of the propagating front V,has to be found in an iterative way to satisfy = 0 at the origin. eq 21 with
Analysis of the Governing Equations Differential eq 7 and 8 have been extensively analyzed in the chemical engineering literature for bodies of a finite dimension (Liu and Amundson, 1962). It has been shown that standing waves usually occur in finite “diffusionexothermic reaction” systems (Hlavacek and Marek, 1971). For large-system dimensions, propagating waves have been also demonstrated. Merzhanov and Averson (1971) showed that, for zero-order reaction and weakly supercritical conditions, degenerate standing-propagating waves exist, while for strongly supercritical conditions a propagating wave results. Qualitative behavior of “diffusion-reaction” system, specified by eq 7-10, is given by the rate of heat conduction and/or heat generation. Consider the case of a very slow chemical reaction. The initial temperature 4 gives rise to a traveling temperature wave which is approximately described by the solution of
7
7
By use of transformation 17, eq 7 and 8 can be rewritten: 7
O 0:
u=o
u--”
0=0,
%=Bi 0=0,
(24)
Since the heat conduction process is much faster than the exothermic chemical reaction, the heat generation term in eq 7 can be neglected. The solution of eq 23 is given by 0 = Bo
-
with boundary conditions (-03
-m
B
= 0i
0=0,
For numerical purposes, the interval of t E (-a,+a) must be replaced by a finite domain. An appropriate transformation can “compress” $, in a finite interval. We have considered, for instance,
+ (0, - 0,)
erfc
(5)
(25)
Evidently, there are two stages of the process: (i) preheating of the solid which is much faster than the chemical reaction, and (ii) slow chemical reaction taking place at 0 = Oi. This situation describes the “sintering” regime of the reaction between the solids. The conversion profile can be easily calculated from eq 8 for 0 = 0i. The temperature, Bi, is too low to ignite a self-propagating regime. Temperature and conversion profiles for such a degener-
Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987 1427 ated situation can be calculated from equations describing a batch process: do dr
- = (1 - 7 ) exp
1961). In the following analysis, these equations will be considered in their steady-state form. After substitution of eq 19 in eq 18, the resulting equation can be integrated from f to +m to obtain
Vf17
do
df + Vf(0 - eo) - Y =0
d17 dr
- = y ( l - 17) exp
For higher ignition temperature, the reaction can be self-sustaining. After the transient start-up period has passed, a constant-pattern profile regime develops (sometimes called a planar wave, creeping profile, or traveling wave) and propagates through the system with constant velocity. Constant-pattern profiles have been predicted, analyzed, and experimentally observed in chemical exothermic reacting systems by many authors: in homogeneous gas systems by Frank-Kamenetskij (1969), Zeldovich and Barenblatt (1959), and Bush and Fendell (1974); in heterogeneous catalytic systems by Fieguth and Wicke (1971), Puszynski and Hlavacek (1984), and Rhee et al. (1973); and in heterogeneous nonctalytic systems by Merzhanov et al. (1972,1973) and Hardt and Phung (1973). The condition for existence of constant-pattern profiles can be inferred from relations developed for a homogeneous explosion
For large values of activation energy, the reaction occurs in a very narrow layer and the reaction temperatures can be approximated by the adiabatic temperature. For such a condition, the second term in eq 29 can be replaced by
Vf (30) Y in the reaction zone. Here l / y represents a dimensionless adiabatic temperature rise. Based on this assumption, eq 29 yields Vf(0 - 0,)
Since the reaction zone is infinitesimally narrow, the temperature change within the reaction zone can be neglected. Evidently, the second term in eq 18 may be omitted. Setting dO/df = q and using the Frank-Kamenetskij approximation, exp[o/(l + DO) N exp 0 , eq 18 can be integrated: q+ - q- =
where 2 ocr = $0 + 1 - 2p For y C ycr,the self-propagating combustion front exists.
A degenerate type of “constant-pattern profiles” may also exist. This situation occurs if the rate of the exothermic reaction in the cold part of the system is comparable with the rate of the reaction front propagation. A traveling wave is created by the hot source at u = 0; however, the already preheated “cold” part of the system gives rise to a finite reaction rate. As a result, a mixed travelingstanding wave is observed. In “asymmetrical”reacting systems (adiabatic systems with Le # 1 or nonadiabatic systems), instability of reaction fronts was observed by Zeldovich and Barenblatt (1959), Hlavaceck and Votruba (1978), and Matkowsky and Sivashinsky (1978). The traveling wave may lose stability via a Hopf bifurcation. The resulting temperature and conversion profiles feature oscillatory behavior. Under certain conditions, the three-dimensional model may become unstable to a two- or three-dimensional perturbation and a variety of reaction front phenomena may result, for instance, pulsating fronts with traveling waves along the front or pulsating fronts with standing waves along the front. These phenomena were experimentally observed (Maksimov et al., 1981) and also recently analyzed by methods of activation energy asymptotics (Margolis, 1983). Our detailed nonlinear analysis will be published elsewhere. In the next two sections, we will develop analytical formulas for the speed of the constant-pattern profile in adiabatic and nonadiabatic systems and for the onset of instability via the Hopf bifurcation.
Propagating Front Velocity Adiabatic Case. The velocity of a propagating constant-pattern profile can be approximately determined from the asymptotic behavior of eq 18 and 19 (Novozhilov,
-
Y -$ Vf
8+ 0-
exp 0 do
The integration limits can be found from eq 31. Assuming that the reaction occurs only inside a narrow layer (0- = 0, $+ = -a), the velocity of the propagating front reads
Vf = Y (33) This relation does not give very accurate results for finite values of the activation energy. A better agreement with accurate results can be obtained from eq 34
v,= y ( l - p + 202)
(34)
The exact value of V , can be calculated numerically from eq 35
This equation was solved by using the Gear’s integration subroutine with the following boundary conditions: 0 = $0 + E: q = o
o = oi:
q = l
Here t is a very small positive number typically t = 0.01 + 0.05. For t = 0, the right-hand side in eq 35 becomes singular and a numerical solution cannot be obtained. Nonadiabatic Case. In this section, relations between the velocity of a propagating front, the hot-spot temperature, and the critical heat-transfer coefficient will be developed. We also report the critical values of these parameters for which an extinction process will occur. The development of the above-mentioned relations is based on an assumption that the width of the reaction zone is narrow in comparison with the preheating and cooling zones. The energy balance of the preheating and cooling zones reads
1428 Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987
d2T dx2
The relation for a dimensionless critical value of the heat transfer coefficient can be found from eq 41
~ K H -(TTi) = 0 d
dT dx
X -- p C p V - -
The boundary conditions are (a) preheating zone
T=To
x=-w
Numerical calculations revealed that the agreement between the numerically computed and analytically calculated values of critical parameters is very good.
T=T,
x=O
(b) cooling zone
T=Tm
x=O
Stability Analysis of Gasless Combustion Systematic studies of stability in gasless systems in which a steady propagating combustion front becomes unstable were done by Maksimov and Shkadinsky (1971) and later by Matkowsky and Sivashinsky (1978). The analysis is based on the assumption that the reaction zone is infinitely narrow. The analytically developed stability criteria for gasless combustion read for constant-pattern profile p < pCr
T=Ti
x-+m
The solution of eq 36 is
p
The heat balance for the reaction zone can be written
for oscillatory behavior
> pc,
where pcr =
4 + 2(5)112
(47)
and After substitution of the first derivative in eq 40, we get -
m
or in the dimensionless form '
L""d(pCpV)2J
For the nonadiabatic case, the velocity of the propagating combustion front is
(
V = C1exp -Tm:2)
Fmc(Pc, ~ c =)
where C1 is a constant. A simultaneous solution of eq 41 and 42 yields the velocity, V , and the hot-spot temperature, T,, for a given value of the heat-transfer coefficient, KH. The exponential function, eq 42, may intersect the function, eq 41, at more than one point. However, for certain parameter values, there is no intersection. Evidently, for the limiting case, the two intersections collapse into one. This condition can be easily calculated from the conditions f(Tm)= f'(Tm)= 0
Shkadinsky et al. (1971) correlated the condition of stability for the first-order reaction model, based on numerical calculations as follows
(43)
where f ( T...J =
,
"
pmc> pmc
2, an increase in the reaction heat results in a more stable combustion, while, for T,/To < 2, an increase in the reaction heat gives rise to a less stable combustion process. Figure 1 indicates in the parametric plane "0-y" the regions of oscillatory combustion, constant-pattern profiles, and degenerate combustion. This figure was constructed by plotting the relations 28 and 50. The accuracy and validity of these criteria are discussed later in this text, based on our detailed numerical study embracing a large domain of parameters. The effect of the heat loss parameter, a , on stability of gasless combustion was investigated by Skhadinsky and
Ind. Eng. Chem. Res.,-Vol. 26, No. 7 , 1987 1429 0.3
Table I. Adiabatic Temperature Rise, Heat of Reaction, and Melting Temperature for Certain Exothermic Noncatalytic Reactions of the Type s-s and s-g compd HfB, TiB TaB, ZrC NbC TaC TiSi TaSi2 HfN ZrN TiN NbN
r 0.2
t
10.1
I OSC I LLATORV
I l l 1 1
REGIME
AT,, K 3520 3350 2700 3400 2800 2700 2000 1800 5100 4900 4900 3500
,3“-
T,, K
kJ/kg 1676.0 2710.9 1034.9 1948.4 1340.8 754.2 1751.4 490.2 1919.0 3477.7 5447.0 2216.5
3520 2500 3370 3690 3750 4270 2190 2470 3580 3220 3220 2740
Table 11. Magnitudes of Activation Energy and Dimensionless Parameters (3, y, and Bi for Certain Noncatalytic Reactions reaction E , kJ/mol p 0.0 0.1
0.2
c
0.3
Figure 1. Phase plane of instabilities in adiabatic systems.
-- -
Moo3 + 3B MOB + BzO3 Ti 2B TiB, Zr 2B ZrB; Hf + 2B HfBz Hf + 2B + O.5HfB2 4 HfB2 Hf + 2B + 0.65HfBz -.+ HfB2
+ +
-+
20
113.1 318.4 310.1 398.1 398.1 398.1
0.21 0.091 0.098 0.08 0.056 0.045
Bo
0.29 -4.4 0.10 -10.2 0.11 -9.4 0.087 -11.6 0.063 -15.83 0.051 -19.50
40
I
C
20
LO
Figure 3. Time-space conversion and temperature profiles ( p = 0.1, = 0.2, o0 = -5, ei = 0, a = 0). l: (A) 0.5, (B) 45.5, (c)93.5, (D) 141.5, (E) 189.5.
10.‘ 0.1
0.2
-
r
0.3
Figure 2. Phase plane of instabilities in nonadiabatic systems for (3 0.
Khaikin (1972). The criterion for the Hopf bifurcation point is 4
+ 2[ 4 + (1 +
5)]1’2
a
(51) .- ,
_1 -Ycr
1+-
Ycr
2
A steadily propagating combustion front is stable for y > ycr. For < ycr,the linear stability analysis predicts that small perturbations grow with time and the combustion process becomes unstable; for an unstable front, temperature and velocity oscillate. For y cycT. Both numerical and approximate calculations are in good agreement for different values of /3 and y. Occurrence of Oscillating Fronts. For supercritical conditions (Shkadinsky et al., 1971),oscillatory fronts can be observed. A typical behavior in this region is shown in Figure 9. The temperature can overshoot the adiabatic temperature, and at these high temperatures, the solid material is completely burned. This is the onset of the cooling cycle. The temperature drops well below the adiabatic temperature, and the fresh unreacted material is slowly preheated. After the temperature of the fresh material has reached certain critical value, a new ignition process occurs and the temperature in the reaction front exceeds the adiabatic temperature again, etc. With increasing values of the activation energy, the profiles become very steep (see Figure 10). An extinction of the oscillatory regime may occur in the nonadiabatic case. For low values of the heat loss parameter, oscillations show behavior similar to that for the adiabatic situation (see Figure 11). However, for the nonadiabatic case, the maximum temperature in the reaction front can be significantly higher. This can be easily explained since there is a lower conversion in the preheating zone so that at the moment of ignition more solid material can react in the very narrow reaction front. For
-8. 100
Figure 11. Time-space conversion and temperature profiles (6 = = 0.5, A~ = io. 0.03, = 0.08, eo = 9.5, ei = 3, a = o.oooi).
higher values of the cooling parameter, a , temperature and conversion profiles become steeper and the hot-spot temperature in the reaction front increases. For a > cy,,, extinctions of the reaction front takes place (see Figure 12). Oscillation and extinction phenomena are highly undesirable in practical applications because of nonuniform product properties resulting from melting effects and/or thermal stresses. To avoid these phenomena, the solid mixture should be preheated and ignited. Figure 13 shows concentration and temperature profiles for a higher ignition temperature Bo (compare Figure 10). These profiles were calculated using 1200 equidistant grid points and a numerical scheme of the type O(h4, k2). For a lower number of grid points, numerical oscillations at temperature profiles were observed. The computer time required for a calculation of wave propagation in an oscillatory regime was 2 orders of magnitude or more longer (it depends on complexity) than for a constant-pattern profile regime and varied from 60 to 3600 s on a CYBER 730. The effect of the reaction heat was discussed earlier. We have shown, based on the analysis done by Matkowsky and Sivashinsky (1978), that two different regimes exist, depending on the ratio of adiabatic and initial temperatures.
1432 Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987
I
n 1 .a Y 0.5
0.5
0.0 100
50
U
150
3.0
80
I
I
4.0
0.0 20
-4.0
60
40
L
80
Figure 14. Time-space conversion and temperature profiles-low heat of reaction (0= 0.0469, y = 0.5, Bo = -8, Bi = -6, a = 0). T, =
-8.0
0.5, A 7 = 500. 100
50
U
150
Figure 12. Time-space conversion and temperature profiles (6 = 0.03, y = 0.08, 00 = -9.5, Bi = 3, a = 0.001). 7: (A) 20, (B) 30, ('2) 40, (D) 450, (E) 480, (F)510, (G) 540, (H)1000.
2.0 Yf
1.0
1.0 17
0.0
0. 0 20.
40.
u "f 5.0
10.
5. 0.0 2
1 15.0
Q 20.
40.
u
"f
Figure 13. Time-space conversion and temperature profiles (6 = 0.03,y = 0.08, 60' = -3, 8, = 3, CY = 0.0). 7 1 = 0.5, AT = 0.5.
For Ta/ To< 2, an increase of the combustion temperature results in an oscillatory behavior. Such an experiment can be carried out by dilution of the reactants by products. Numerical simulation reveals that the stability criteria cannot be applied for both very high and very low values of the reaction heat. For small values, of the heat of reaction, constant-pattern profiles or oscillatory combustion cannot occur and the extent of the reaction is limited by heat conduction from the ignition source (cf. Figure 14). In this case, a forced burning of the substance with continuous heat supply from the source takes place. The self-propagating stage is absent. Figure 15 shows the behavior of systems with identical values of the activation energy and initial temperature but with different values of the heat of reaction. An increase of the reaction heat results in a stable combustion characterized by a constant velocity of the propagating front (cf. Figure 15a). A further increase of -AHleads to oscillations (cf. Figure 15b,c) and finally again to a stable combustion. Table I11 reports a comparison of numerically calculated regions of stable and unstable combustion with those predicted from eq 47 and 50. When a finite cylindrical sample is ignited at both ends, two reaction fronts move in opposite directions toward each other. If the combustion takes place in a subcritical regime, the waves initially move as constant-pattern pro-
10.0
5.0
0.5
1.0
Figure 15. Dependence of dimensionless combustion velocity vs. time for different values of reaction heat (/3 = 0.0469, 6, = -8, a = 0). (a, top) (1) y = 0.125, Oi = 0; (2) y = 0.075, Bi = 5.33. (b, middle) y = 0.05, Bi = 12. (c, bottom) y = 0.04, Bi = 17. Table 111. Comparison of Results of Numerical Calculations with Those Predicted by Equations 47 a n d 50 (6'0
= -8, 6. = 0.0469) 7
e;
0.500 0.125 0.075 0.050 0.040 0.035 0.030
-6.00
0.00 5.33 12.00 17.00 20.55 25.33
L1
ump
3.87 8.00 8.53 8.18 7.74 7.402 6.96
2.265 1.020 0.920 0.928 0.965 0.9985 1.050
numerical simulation slow chemical reaction regime constant-pattern profile oscillatory behavior oscillatory behavior oscillatory behavior oscillatory behavior constant-pattern profile
files. Later, however, since the unreacted material is heated by conduction from both reaction fronts, the tem-
Ind. Eng. Chem. Res., Vol. 26, No. 7, 1987 1433
25
50
U
50
U
7
5.0
0.0
-5.0 25
F i g u r e 16. Time-space conversion and temperature profiles (twosided ignition) (0= 0.1, y = 0.2, Bo = -5, Bi = 0 , a = 0). r1 = 0.5,AT = 24. 1.0 ll
0.0
1
0
7 10.0 5.0
0.0
-5.0
25
50
75
100
U
F i g u r e 17. Time-space conversion and temperature profiles (twosided ignition) (@= 0.03, y = 0.08, 8, = -9.5, Bi = 3, a = 0). il = 0.5, AT = 11.
perature in the middle part of the sample gradually increases. The collision of these two propagating waves results in a temperature in the middle of the sample which is significantly higher than the adiabatic temperature (cf. Figure 16). More pronounced effects are observed for a supercritical regime (cf. Figure 17). Conclusions Solid-solid noncatalytic reactions featuring high values of both activation energy and reaction heat can occur in different regimes, depending on values of these parameters and operating conditions. Numerical analysis reveals that the following combustion regimes can be expected: (i) degenerate combustion; (ii) constant-pattern profile regime; (iii) simple oscillatory behavior; and (iv) multipleperiod oscillatory combustion. For high values of the above-mentioned parameters, the numerical solution of governing equations gives rise to steep temperature and concentration profiles and very narrow reaction fronts. Since the reaction takes place in a thin zone, the numerical solution of the problem represents a formidable task. Sophisticated schemes have been devised which take advantage of automatic time-step adjustment, high-order approximation of the space derivatives, and nonuniform distribution of grid points in space. Analytical conditions have been developed describing a transition between degenerate and constant-pattern profile combustion. Approximate formulas for the velocity of the front of the constant-pattern profile regime have
been derived, both for the adiabatic and nonadiabatic conditions. A numerical technique to determine an exact value of the front velocity has been proposed and tested. A comparison of approximate and exact numerical values of the velocity reveals good accuracy of the analytical expressions. For a nonadiabatic case, explicit conditions in an analytical form have been found that guarantee constantpattern type of propagation of the reaction front. If these conditions are violated, an extinction process takes place. For very high values of the governing parameters, a Hopf bifurcation can occur. As a result, oscillations in the system take place. Analytical formulas, based on the activation-energy asymptotic approximation, have been scrutinized and the results compared with the exact numerical values. The investigation reveals that the asymptotic analysis predicts in a qualitative way the transition, although, the exact values must be determined numerically. From our numerical simulation of a nonadiabatic oscillatory case, it can be inferred that the temperature profiles are much steeper than in the adiabatic case and the hot-spot temperature significantly exceeds the adiabatic combustion temperature. Acknowledgment This work was supported by NSF Grant CPE83-10627. Nomenclature C = concentration, kg m-3 C, = heat capacity, J kg-' K-l d = diameter, m E = activation energy, J kmol-' AH = heat of reaction, J kg-' ko = rate constant K H = heat-transfer coefficient, J m-2 s-' K-' n, m = reaction orders r = radial coordinate, m R = universal gas constant, J kmol-' K-l t = time, s T = temperature, K u = dimensionless coordinate V = velocity of combustion front, ms-' V, = dimensionless velocity of combustion front x , z = coordinates, m Greek Symbols a = dimensionless heat-transfer coefficient p = dimensionless activation energy y = dimensionless heat of reaction 7 = conversion X = heat conductivity, J m-' s-' K-' k , F~~ = bifurcation paramters, eq 47 and 50 v = frequency of oscillations, s-' 8 = dimensionless temperature p = density, kg m-3 4 = angular coordinate T = dimensionless time Subscripts a = adiabatic cr = critical i = ignition 0 = initial m = maximum s = solid * = scale Literature Cited Averson, A. E.; Barzykin, V. V.; Merzhanov, A. G. Inz.-Fiz. Zh. 1965, 9, 244 (in Russian). Borovinskaja, I. P.; Merzhanov, A. G.; Novikov, N. P.; Filolenko, A. K. Fiz. Goreniya Vzryua 1971, 7, 454 (in Russian). Bush, W. B.; Fendell, F. E. Combust. Sci. Technol. 1974,17,1499.
1434
I n d . Eng. Chem. Res. 1987, 26, 1434-1441
Ciment, M.; Leventhal, S.; Weinberg, B. J . Comp. Phys. 1978, 28, 135. Fieguth, P.; Wicke, E. Chem.-1ng.-Tech. 1971,43,604 (in German). Frank-Kamenetskij, D. A. Diffusion and Heat Transfer i n Chemical Kinetics, 2nd ed.; Plenum: New York-London, 1969; p 422. Hardt, A. P.; Phung, P. V. Combust. Flame 1973, 21, 77. Hlavacek, V.; Marek, M. Proceeding of the Fourth European Symposium, Brussels, 1968; Pergamon: Oxford-New York, 1971; Section 11, Paper 5. Hlavacek, V.; Votruba, J. Ado. Catal. 1978, 27, 59. Liu, S. L.; Amundson, N. R. Znd. Eng. Chem. Fundam. 1962, I , 200. Maksimov, E. I.; Shkadinsky, K. G. Fir. Goreniya Vzryua 1971, 3, 459. Maksimov, Yu. M.; Merzhanov, A. G.; Pak, A. T.; Kuchkin, M. N. Combust. Explos. Shock Waves (Engl. Transl.) 1981, 17, 393. Margolis, S. B. SIAM J . Appl. Math. 1983, 43, 351. Matkowsky, B. Y.; Sivashinsky, G. I. SIAM J . Appl. Math. 1978,35, 465. Merzhanov, A. G. Arch. Combust. Process. 1974,5, 17 (in Russian). Merzhanov, A. G.; Averson, A. E. Combust. Flame 1971, 26, 89.
Merzhanov, A. G.; Filolenko, A. K.; Volodin, Yu E. Dokl. Akad. Nauk SSSR 1972, 206, 905 (in Russian). Merzhanov, A. G.; Filolenko, A. K.; Borovinskaja, I. P. Dokl. Akad. Nauk SSSR 1973,208, 892 (in Russian). Novikov, N. P.; Borovinskaja, I. P.; Merzhanov, A. G. Fir. Goreniya Vzryua 1974, 2, 201 (in Russian). Novozhilov, B. V. Dokl. Akad. Nauk SSSR 1961, 142, 151 (in Russian). Puszynski, J.; Hlavacek, V., unpublished results, 1987. Puszynski, J.; Hlavacek, V. Chem. Eng. Sei. 1984, 39, 681. Rhee, H. K.; Foley, D.; Amundson, N. R. Chem. Eng. Sei. 1973,28, 607. Shkadinsky, K. G.; Khaikin, B. I.; Merzhanov, A. G. Combust. Explos. Shock Waues (Engl. Transl.) 1971, 7, 15. Shkadinsky, K. G.; Khaikin, B. I. Symposium on Combustion Processes, Moskva, 1972 Zeldovich, Yu. B.; Barenblatt, G. I. Combust. Flame 1959, 2, 61.
Received f o r review September 13, 1985 Accepted March 2, 1987
An Experimental Investigation of Oxygen Enrichment in a Silicone Capillary Permeator with Permeate Recycle Sudipto Majumdar, Lawrence B. Heit,+Amitava Sengupta, and Kamalesh K. Sirkar* Department of Chemistry and Chemical Engineering, Stevens Institute o f Technology, Hoboken, N e w Jersey 07030
Experimental investigations have been carried out on oxygen enrichment of air by permeation through silicone-capillary membranes in a permeator, with partial or total recycle of the oxygen-enriched permeate to the fresh-air feed. A vacuum recycle mode was utilized, and permeate pressures were in the range from 1.013 X lo3to 4.053 X lo3 P a (0.01-0.04 atm) for feed air pressure varying between 1.023 X lo5 and 1.114 X lo5 P a (1.01-1.10 atm). Permeator performances are determined for the recycle ratios of 0.0, 0,1990, 0.6067, 0.7462, 0.9039, and 1.0. A model incorporating axial pressure drop and capillary deformation is proposed for the recycle permeator. Model simulations agree with the experimental data quite well for all recycle ratios. For high recycle ratios, the recycle permeator displays a maximum in oxygen enrichment a t a finite net stage cut. This maximum increases as the recycle ratio is increased. A recycle permeator is a convenient way to produce a richer permeate than can be obtained from a permeator operating in the conventional mode. However, the compressor power and membrane area required per unit net permeated product increases with recycle. Membrane gas separation is being increasinly adopted by the chemical process industry (MacLean and Graham, 1980). The gas-separating membranes used in such processes do not usually possess the desirable combination of both high selectivity and high permeability. Therefore, considerable attention has been focused on the development of efficient separation schemes which use the available membranes (Matson et al., 1983; Sengupta and Sirkar, 1986). Such separation schemes may utilize two different membranes in the same permeator, with reversed membrane selectivities for different feed gas components (Kimura et al., 1973; Sirkar, 1980; Stern et al., 1984; Sengupta and Sirkar, 1984; Stern and Perrin, 1985; Sengupta and Sirkar, 1987). Alternately, a higher overall separation factor can be achieved if one can effectively increase the concentration of the more permeable species in the high-pressure feed stream. In one such scheme, a fraction of the permeate is recycled and mixed with the high-pressure feed gas. The experimental investigation and simulation of such a scheme, widely followed in industry (Schell and Houston, 1983; Kilgour and Hilliard, 1983),has
* To whom
correspondence should be addressed. Ciba-Geigy Ltd., Summit, N J 07901.
0888-5885187 /2626-1434$01.50/0
been studied here. Oxygen enrichment of air, known to be very useful in improving combustion efficiency, has been chosen as the gas-separation system for study. Silicone rubber membranes are used often to obtain 02-enriched air due to their high permeability to O2 and moderate permselectivity for O2over N2 (Ward et al., 1976; Hedman, 1982). A silicone rubber membrane permeator operated conventionally with a fresh-air feed has, however, an upper limit of O2 enrichment. Higher d2enrichment of air is possible if a fraction of the permeated 02-enriched air is recycled and mixed with the fresh-air feed. For a given fresh-air feed flow rate, the net rate of production of 02-enriched air is decreased by such partial permeate recycle. The net stage cut (defined as the ratio of the net production rate of the 02-enrichedair to the fresh-air feed flow rate) with recycle is lower than that without recycle. How the net stage cut varies with the improved composition of the permeate is of considerable interest. Of particular interest is an experimental verification of a maximum in the permeate mole fraction of O2 when plotted against the net stage cut for a given recycle ratio (defined as the fraction of the total permeate recycled to mix with the fresh feed). The numerical simulations of Teslik and Sirkar (1981),Teslik (1983),McCandless (1985), and Stern et al. (1984) suggest such behavior. This is in 1987 American Chemical Society
