Ind. Eng. Chem. Fundam. 1980, 79, Hoyaux, M. F., “Arc Physics”, p 172, Springer-Verlag, Berlin, 1968. “J.A.N.A.F. Thermochemical Tables”, 2nd ed,NSRDS-NBS37 National Bureau of Standards, Washington, D.C., 1970. Jordan, G. R., Bowman, B., Wakelam, D., J. Phys. D. Appl. Phys., 3 , 1089 (1970). Langmuir, I., Phys. Rev., 34, 401 (1912). Lincoln, K. A., Covington, M. A., Int. J. Mass Spectrom. Ion Phys., 18, 191 (1975). Lundell, J. H., Dickey, R. R. A . I . A . A . J., 11, 216 (1973). Lundell, J. H., Dickey, R. R., Prog. Astronaut. Aeronaut., 56, 405 (1977). Matsunawa, A., et al. Trans. J. WRI(Jpn), 6, 161 (1977). Mentel, J. Appl. Phys., 14, 269, 361 (1977a, b). Mentel, J. Appl. Phys., 15, 179 (1978). Meyer, J., Brit. J. Appl. Phys., Ser. 2 , 2 , 221 (1969). Nishiguchi, K., Matsunawa, A. Techno/. Rep. Osaka Univ., 19, No. 897 (1969) (Figures 15 and 16); 21, No. 987 (1971).
243-250
243
Ozisik, M. N. “Boundary Value Problems of Heat Conduction”, pp 301, 326, International Textbook Co., Scranton, Pa., 1968. Patankar. S. V., Spaldins, D. B., “Heat and Mass Transfer in Boundary Layers”, 2nd ed, p 71, Intertext, 1970. Schwabe, W. E. J. Met., 24, 65 (1972). Sparrow, E. M., Gregg, J. L., Trans. A . S . M . E . ,78, 1823 (1956). Vargattik, N. B., Vasilevskaya, Yu, D., translation of Teplof. Vys. Temp.. 7 , 913 (1969). Ward,-d. A:, Ph.D. Dissertation, University of Canterbury, N.Z., 1974. Whittaker, A. G., Kintner, P. L., Carbon, 14, 257 (1976). Wiles, P. G., Abrahamson, J., Carbon, 16, 341 (1978). Wiles, P. G., Ph.D. Dissertation, University of Canterbury, N.Z., 1979.
Received for review June 12, 1979 Accepted D e c e m b e r 27, 1979
Diffusion and Reaction in a Stagnant Boundary Layer about a Carbon Particle. 4. The Dynamical Behavior Eduardo Mon University of Minnesota, Minneapolis, Minnesota 55455
Neal R. Amundson’ University of Houston, Houston, Texas 77004
A transient model of single-particle carbon combustion is presented based on the assumption that the surrounding gas-phase boundary layer is in a quasi-steady state. The analysis is subdivided according to the particle size into a one-film and a continuous flame model. I n the latter, more general case, the model may be formulated as an initial value problem at the boundary of a nonlinear ordinary differential equation. The numerical solution reveals the combustion history of the particle from ignition to complete burn-off which displays a number of unsuspected and occasionally anomalous features of behavior. Some of these include multiple bifurcations of the quasi-steady state, transitions between such states, and nonuniqueness of the distributed transient model.
Introduction The physical system considered in this paper is that of an impervious carbon particle surrounded by a stagnant boundary layer through which oxygen diffuses from the ambient gas and reacts with the carbon surface to produce mainly carbon monoxide in the temperature range of interest. The CO so obtained is oxidized homogeneously to carbon dioxide in the boundary layer whereas the resulting COz diffuses in part to the carbon surface where it is reduced to carbon monoxide. The latter reaction is endothermic, the first two exothermic, so that in the mathematical model we must consider multicomponent diffusion-reaction coupled with heat conduction across the boundary layer. The quasi-steady state behavior of this system was considered by Caram and Amundson (1977) and Mon and Amundson (1978). Both studies showed that multiple steady states result for certain parameter combinations, their stability character having later been examined by Mon and Amundson (1979). In the present study we consider the actual dynamics of the system-the history of the particle and combustion process from ignition to complete burn-off. This work necessarily draws to a great extent on the results of the previously mentioned papers, which should be consulted before proceeding with this one. Preliminary Development In the following treatment it is assumed, given the relatively large difference in thermal time constants be0196-4313/80/1019-0243$01.00/0
Table I.
Notation
in-
spe-
dex
cies
i = 1 CO, i = 2 CO i = 3 0,
reaction
rate expression
2C + 0, = 2CO R, = k,e-Ei/RTc3 C + CO, = 2CO R , = k,e-E2/RTc, CO + 1/20, = CO, R = k , e - E 3 ’ R T ~ , ~ , ” Z
tween the solid particle and the gas-phase boundary layer, that the gas is in a quasi-steady state whether the system of particle and boundary layer is in a quasi-steady state or not. Some discussion of this is in order for, in fact, the burning of a particle is always a transient process in which the particle shrinks in size and eventually burns out. However, we can envisage as in the previous papers that a t any instant the particle is in a steady state with its surroundings and compute what the particle temperature will be. As we will show in this paper, there is a locus of steady states for fixed ambient conditions as the particle size decreases, this locus having, in general, multiple branches. It is to be stressed that any point on the locus would be a steady state for a particle of that size, if the size could be instantaneously fixed. But suppose a particle with a given size and some other temperature were placed in the same environment,.that particle then would not, in general, be in a quasi-steady state. However, it is assumed here that the boundary layer assumes a steady state instantaneously consistent with the particle temperature and the surrounding conditions. The state of such a system 0 1980 American Chemical
Society
244
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980
Table 11. Phvsical Data Used in the Calculations k , =4.016 X
kg-mol/m3 = 1 . 7 5 x 104 J/kg-mol K D =2X m2/s E, = 17966 K E ? = 29790 K E:, = 15098 K f = 0.98 k , = 3.007 x l o 5 m/s
k , = 8.107 X 10' m3'2/kgmol''2s k = 6.7 X lo-' J/m s K Q1 = 2.24 X l o g J/kg-mol Q, = -1.70 x 10' J/kgmol 9, = 2.82 X 10' J/kgmol p c = 0.8 X l o 3 kg/m3 T M = 1400 K
c,
Carrying out the differentiation in (6) and using (8)gives a more explicit form for the boundary condition on the energy flux d T, (12) elr=a(t) = RsHc - cs - qR
log m/s
c = 1.21 x
dt
+
will not be on the locus and the system will pass through a transient, the particle temperature changing, and with the boundary layer always assuming a steady state but changing as the particle temperature changes. The difference then in the structure of the solution will be that for points on the locus the time derivative of the particle temperature will be zero while for points not on the locus it will not be zero. The reactions that occur at the surface of the particle and in the boundary layer are summarized in Table I, along with the species present and rate expressions that will be used. The temperature of the particle will be taken to be uniform. We neglect the Stefan flow and assume the physical properties to be constant. Specific values used in the calculations are given in Table 11. The mass and energy balances in the time dependent domain a ( t ) C r C b ( t ) have the form I d --(rzNJ = uiR i = 1, 2, 3 r2 dr I d =0 --(?e) r2 dr The boundary conditions associated with the mass balance at the particle surface, r = a ( t ) are N1 = -Rz (3) Nz = 2Rz 2R1 (4)
+
N3 = -R1 (5) whereas the energy boundary condition at the surface takes the form
A[ dt
% 3 z 3 ( t ) ~ E c +] 4m2(t)(e+ QR) = 0
where R, is the combustion rate, 2R1 R2. Notice that when dT,/dt = 0, eq 12 reduces to the boundary condition of the quasi-steady state model, this being the major difference in structure between the present and the previous model. The moving boundary is conveniently immobilized by taking the boundary layer thickness to be linearly proportional to the particle radius and by using the coordinate transformation f = a(t)/r. It is then possible to show that the aid of the stoichiometric coupling of the molar fluxes and the Stefan-Maxwell equations, as done in our first paper, that the following relations hold in the fixed domain a / b C f C 1
The energy flux with respect to stationary coordinates was similarly shown to be
e
( ):
-k-d T - N1Q3 + f 2 R , H2- dr
(15)
The energy equation (2), when integrated with the aid of the boundary condition (12) results in
Substitution of (15) into (16), use of the proper flux expression for C 0 2 ,and setting dr = -(a/t2)df gives
where Q1 1 -AHl = H3/2 + H, - H2, a3 = cDQ3/k, and L = pcC,/3kM,. The integration of (17) provides a relation between the C 0 2 and temperature profiles
(6)
subject to
T, = Ts0 t = 0 Similarly, the carbon balance is
(7)
It should also be noted that (17) when evaluated at the particle surface and combined with (3) yields a more useful formulation of the energy boundary condition
with t=O (9) At the external edge of the boundary layer, r = b ( t ) ,we have xj = Xib i = 1, 2, 3 (10) T = Tb a=ao
The radiant energy flux q R , is q R = g ~ ( T 2 TM') 7
(11)
where TMis a model parameter. We refer to the condition TM= Tbas an isolated or single particle and the condition TM# Tb as an interacting particle.
The problem to be solved then consists of the above equation and its corresponding initial condition (7), as well as the shrinkage eq 8 and 9. At any given time t , evaluation of the right-hand side of these equations requires the computation of the temperature profile in the boundary layer to evaluate (dT/df) and the mole fractions at the surface XI,,xZs,x3,. To ogtain these quantities, we must solve simultaneously, and in conjunction with the above initial value problem, the distributed conservations (1)and (2). However, with the foregoing development it is possible to reiterate the model in a simpler form. For instance, we have in eq 13,14, and 18 three relations between four of
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980 245
the unknown profiles so that it is only necessary to solve one of the four (1)and (2). Furthermore, the model is additionally simplified when the nonlinear rate term R vanishes. This motivates us to subdivide further development of the model according to whether or not a CO flame develops in the boundary layer. Part I. A One-Film Model Small particles, in the pulverized fuel range, may be treated by a one-film model since generally the CO oxidation reaction does not take place in the boundary layer. Under these circumstances, the mass balances can be integrated directly to give, as in the steady state case
24 L
c
SS1-c
lo IO 12 14 AMBIENT TEMPERATURE (K)XlO-'
where y1 = 1, y2 = y3 = y. The mole fractions a t the surface of the particle are obtained by writing eq 20 at [ = 1and use of the boundary conditions (3)-(5). This yields
6k2 Yi D where Pki is the stoichiometric coefficient of the ith component in the kth heterogeneous reaction and 6 = 1- a / b. Equation 18, when evaluated a t the surface, becomes P2i
-a
- exp(-c2/T,)xl, = 0 (21)
The values of xls and xh appearing above can be obtained in terms of the particle size and temperature from (21)
X38
=
X3b
a bki 1 + -- exp(-t,/T,)
(24)
YD Substitution of (23) and (24) into (22) and introduction of the dimensionless particle temperature y = T,/Tso, and the fractional decrease in radius cy = 1 - (a(t)/ao),allow us to formulate the model as dY - = f(Y,d (25) dt dtr -- = g(y,cy) dt with initial conditions y(0) = 0; cy(0) = 0 (27) The Quasi-Steady State. The quasi-steady state values of the particle temperature a t a given particle size are determined by the roots of the nonlinear algebraic equation f(YP) = 0
(28)
In our previous work, the steady-state results were presented by considering the particle temperature to be a function of the ambient temperature at a fixed value of the particle size. In the case of the one-film model it is possible to solve for Tb in terms of T, to obtain a locus determined by
0
20
40
60
80
100
5 DECREASE IN RADIUS
Figure 1. Some loci of quasi-steady state solutions for a particle burning in 15% 02.
The dynamic behavior of the particle, however, is best presented in the plane of particle temperature-particle radius or percent decrease in radius at a fixed value of the ambient temperature. It is therefore of interest to first present the locus of quasi-steady state solutions in such a plane. In this study we will choose a particle with a radius of 50 pm exposed to an atmosphere containing 15% O2and use the expression for the radiant heat flux, qR, for an interacting particle. For our present purposes, there will be no CO or COz in the ambient which means a single reaction takes place at the surface. Treatment of the case where the C-C02 reaction also takes place is by no means more difficult but neither are new qualitative features of behavior introduced. The boundary layer thickness will be fixed at two particle diameters. The resulting locus of solutions in the T,-Tb plane is shown in Figure 1A. The shape of the locus in the T, vs. % decrease in radius plane will depend on the value chosen for the ambient temperature, and the range of possibilities is better illustrated by selecting the three different values of T b indicated by the dotted lines in Figure 1A. They correspond to 1000,1200, and 1300 K. The points numbered in the Ts-Tblocus will thus lie on the ordinate of the T,- % decrease in radius plane. Notice that at this initial value of the particle size there are three different steady states that correspond to the 1000 and 1200 K values of the ambient temperature. At 1300 K, however, there is a unique steady state. The prediction of three steady states with a one-film model has previously been pointed out by Vulis (1961), Field et al. (1967), and Ubhayakar (1976). The T, vs. % decrease in radius plot must be obtained from (28) but we may avoid the solution of this nonlinear equation by generating a series of T,-Tb loci (29) at closely spaced values of the particle size and linearly interpolating the value of the particle temperature that corresponds to the ambient temperature and particle size of interest. The result is shown in Figure 1B for the three chosen Tb values. Along the first steady state (SS1-a), corresponding to 1000 K, T, decreases almost linearly toward Tb. The second steady state (SS2-a) increases in particle temperature at a small rate initially, then more sharply near 40% decrease t o disappear at (42%, 1960 K) as it joins the third steady state, thus forming an upper loop. Past 42% decrease only the lower state is possible. A similar behavior is displayed at Tb = 1200 K (curves labeled b) but now the range of particle size over which all three steady states coexist is enlarged as the upper loop does not disappear until the initial radius has decreased by 82%. Finally, for the T b = 1300 K case (curves labeled c) only a highly ignited upper steady state is possible until the
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980
246 24,
I
-A-
I
% DECREASE IN
RADIUS
% DECREASE IN RADIUS
Figure 2. A select number of transient temperature histories for the particle with the quasi-steady state structures (b) and (a) of Figure
1B.
radius has decreased by 41%. However, for this given size a bifurcation emerges at T, = 1450 K thus giving rise to an additional two steady states. These three states coexist until the radius has shrunk by 88% after which only the lowest state is possible. The most striking feature is that the number of quasi-steady states is not constant throughout the history of the particle, but rather they vary from three to one or from one to three and back to one depending on the ambient conditions. As to what the actual behavior of the particle temperature will be for a given set of initial conditions, and particularly as the particle burns along a given steady state that disappears, can only be answered by the solution of the dynamic model. Transient Behavior. Equations 25-27 were found to be stiff. However, the semi-implicit third-order RungeKutta method with variable step size described by Michelsen (1976) and Villadsen and Michelsen (1978) provided excellent convergence. A set of transients with varying initial particle temperatures is shown in Figure 2A for the 50 pm radius particle with the ambient at 1200 K and 15% 02. When the particle is initially heated to the temperature of the upper state, it burns essentially next to it up to the point where this steady state disappears (transient 2A-a). The temperature then drops from 2323 K at 75% decrease to 1210 K a t 87% in 0.06 s as there is a shift to the lower steady state. Although the lifetime is 1.24 s, the radius decreases by 75% in the first 0.07 s of combustion. With an initial temperature of 1485 K, near the intermediate steady state, combustion proceeds initially parallel to this state but gradually moves away toward the upper state (transient 2A-c). The temperature history then coincides with that of transient 2A-a, burning along the upper state but shifting to the lower state as the upper disappears. The particle lifetime is then 1.30 s. Decreasing the starting temperature by only 1' (Tso= 1484 K) causes the particle to burn initially along the intermediate state but as it drifts away the shift is to the lower state (transient 2A-d). As expected, any initial temperature above 1485 K will cause the particle to approach the upper state, whereas any initial temperature below 1485 K will cause it to approach the lower, of which transients 2A-b and 2A-e are two examples. These results agree with the conclusion that the intermediate steady state is unstable, the upper and lower stable as obtained in our previous stability analysis (Mon and Amundson, 1979). In that analysis we capitalized on the difference in time constants between the gas and the solid phase to show that the stability of the system could be established by observing the dissipative behavior of a perturbation in the gas phase. It is thus interesting to note that the stability of the system translates into the present model where
0
20 %
40 60 SO DECREASE IN RADIUS
0
Figure 3. Effect of the initial particle size on the temperature history for the conditions corresponding to case (a) of Figure 1B. The initial temperature in all cases is 300 K.
we take the gas phase to be in a steady state. Figure 2B shows several transients when the same particle burns at a lower ambient temperature, 1000 K. Transient 2B-a shows the particle initially heated to the temperature of the upper state. The behavior is similar to that of transient 2A-a except that the shift to the lower state does not occur as fast, causing some overshoot of the loop formed by the upper states. When the initial temperature is 1700 K (transient 2B-b), the temperature increases toward the upper steady state but the approach is not fast enough for the particle to reach that state before it disappears. That the intermediate state is unstable is again indicated by transient 2B-c. In case 2B-d (Td = 1400 K) we see that there is only a small decrease in radius in the approach to the lower steady state as a result of the lower combustion rate a t such temperatures. The high ambient temperature case is studied in Figure 3. Transient 3a shows the temperature history of a particle initially at room temperature, 300 K. It reaches 1500 K in 0.06 s and the unique, high-temperature steady state in 0.10 s. Combustion then continues along this state but shifts to the lower as in the previous cases when the upper disappears. Any other particle in the size range 33.5 pm I a. I 50 pm shows similar behavior when started from 300 K (transient 3b). However, at a. = 33 pm, transient 3c, the particle temperature is deflected toward the intermediate steady state, coincides with it temporarily, but eventually drifts to the lower state. Since this particle burns at lower temperatures its lifetime is 1.25 s, whereas the lifetime of the 50-pm particle is less, 0.38 s. This indicates that given two different sized particles (50 and 33 pm) burning under identical conditions, the larger particle can be consumed before the smaller one-a pathological result that appears to defy common sense. Another important feature of these results is the prediction of a transition from the upper to the lower steady state at a given size of the particle. At lower ambient temperatures than 1000 K, the lower steady state is often an unignited steady state. In that case the above transition would correspond to a spontaneous extinction of the burning particle. The phenomenon of spontaneous extinction at a critical size of the particle has been documented experimentally by Ubhayakar and Williams (1976). We must finally warn that this extinction does not necessarily coincide instantaneously with the disappearance of the upper ignited steady state as the lag displayed in transients 2B-a and 2B-b indicates. Part 11. A Continuous Flame Model For larger sized particles, we must consider the effect of the CO oxidation in the boundary layer. In this case
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980
the mass balances cannot be directly integrated as in the previous section. As stated in the preliminary development, to evaluate the right-hand side of (19) and (8) requires the solution of one of the distributed conservation equations (1-2). For numerical convenience we opt for the energy equation. Substitution of the expression for the energy flux (15) in (2) leads, after some straightforward manipulations to d2T kt4a2Q3k3c3/2 e x p ( - ~ ~ / T ) x ~=x0~ ~ (30) I~ dt2 The mole fractions x 2 , x3 can be obtained in terms of the temperature by solving for the C 0 2profile in (18) and using the result in (13) and (14). This gives
+
x2
(T - T b )
+
247
temperature and size history of the particle we must solve the problem defined by (19) and (7)-(9). Using (36) in (8) allows us to write these equations formally as (37)
with y,(O) = 1; a(0) = 0 (39) This initial value problem is coupled through the gradient of the temperature at the surface to the nonlinear ODE (33)-(35), or
= X2b - Ya3
Y(0) = 1; Y(1) = TS/Tb (41) A more convenient form of (40) is obtained by elimination of dy,/dt with (37). The result is
To summarize, the model therefore consists of the simultaneous solution of the following set of equations However, the mole fractions at the surface that appear in the heterogeneous reaction rate expressions above are still unknown. Suppose we evaluate eq 13, 14, and 18 a t the surface of the particle. We then obtain three equations in four unknowns (xis, xZs,x3s, and dT,/dt). As opposed to the steady-state case, the surface temperature does not form part of the unknowns since it is specified by the initial condition (7). The resulting equations at [ = 1 are linear in the mole fractions and we can solve explicity for their values in terms of dTs/dt. Formally
i = 1, 2, 3
(33)
The specific functional dependence above is identical to eq 44-46 of Mon and Amundson (1978) if 7, is redefined as
(34) for a single particle and as
for an interacting particle. It is more useful to rewrite (33) by eliminating dT,/dt with (19)
(43) (44)
subject to the initial conditions y,(O) = 1; a ( 0 ) = 0 (46) and the boundary conditions Y(0) = 1; Y(1) = Tb/T,o (47) Any integration procedure of the initial value problem (43)-(44) requires repeated evaluations of 4 and 6. For instance, the semi-implicit third-order Runge-Kutta procedure described previously requires eleven function evaluations per time step. In turn, each function evaluation requires the solution of the distributed equation (45) that incorporates all the undesirable numerical properties of its corresponding steady-state form discussed in our first paper. The complexity of any computational solution scheme should thus be evident. Most of the results to be presented here were obtained by the use of orthogonal collocation in terms of the ordinates, yi (i = 1, ...,N) at the N zeroes of a suitable Jacobi polynomial (Michelsen and Villadsen, 1978). This converts (45) into a system of nonlinear algebraic equations in N unknowns N+ 1
N+l
C B , j ~ j+ ]
=o
X ( ~ , Y , , Y ~ , ~ L ,
E Aoj.Yj) = 0
1=0
(48)
j = 1, ..., N
which shows that in the transient case, the surface mole fractions are a function of not only the surface temperature but also the gradient of the temperature at the surface. It will be convenient to normalize the domain by introduction of the variable u = (1- [)/(1 - ( a / b ) )and thus formulate the model in dimensionless form. To obtain the
or in matrix form BY + x(~,Y,,u,Y,&) = 0 (49) A given guess of the temperature profile, 9, connecting the known end values y(0) and y ( l ) , allows us to calculate an improved profile by the Newton-Raphson method. Con-
248
Ind. Eng.
Chem. Fundam., Vol. 19, No. 3, 1980
vergence of this iterative scheme leads to a value of (dy/du),=, that can be used in a given evaluation of the right-hand side of (43) and (44). Although this procedure is repeated eleven times in any given time step, convergence is generally very fast due to the slight difference between the profiles within a given step. However, since the profile that results at a specific time step must serve as initial guess for the next step, the length of the time step must be regulated so that the difference between the guessed and actual solution does not cause the NewtonRaphson search to diverge. Notice that this method does not make use of the concept of a feasibility region as used in our previous papers for the steady-state case. Convergence of the solution of the distributed equation generally depends on the steepness of the solution profile and the accuracy of the guessed profile. While the collocation method worked in most cases, it failed for certain very steep profiles obtained when the flame intensity was such as to consume all the oxygen before it diffused to the surface of the particle. Under these circumstances, results were obtained by application of a finite difference method to the distributed problem. This method also converts (45) into an algebraic system and the solution procedure is similar. The time requirements for the finite difference method, however, were substantially greater than for global collocation. T h e Quasi-Steady State. While discussing the quasi-steady state we must note a correction of the T,-Tb loci for the particle represented in Figures 3A and 4A of Mon and Amundson (1978). Using a much finer finite difference grid, it was possible to detect an additional branch of the locus that coincides everywhere in the temperature range presented with the zero oxygen line. Following a similar pattern, we postulate the existence of a similar branch for the particle burning in pure oxygen, Figure 2B, although a numerical solution could not be obtained in the latter case due to the steepness of the temperature profile. As in the case of the one-film model, it is desired to present the locus of quasi-steady state solutions in the plane of particle temperature vsi percent decrease in radius (aX lo2). We w ill choose to study the particle whose T,-Tb locus was previously presented in Figure 4B of our first paper and is reproduced in Figure 4A. This locus corresponds to a particle of 500 pm radius burning in 21% O2 with no CO or C 0 2 in the ambient and a boundary layer thickness of one particle radius. The radiation term is of the single particle type. The three dotted lines in Figure 4A indicate the ambient temperature selected, and each yields a T, vs. CY X lo2 locus. Choosing Tb = 900 K in Figure 4A yields the locus of Figure 4B. Here the points marked 3 to 6 correspond to those similarly marked in Figure 4A. The results have been extended to larger particle sizes so that the initial radius is 1000 pm and the 500-pm particle of Figure 4A corresponds to 50% decrease in Figure 4B. It should be understood that we are omitting the additional unignited state below 1000 K. As may be seen, up to 35% decrease there are only two ignited steady states. Both of these states display significant CO oxidation with a pronounced temperature maximum in the boundary layer. At 35% decrease a bifurcation point appears, giving rise to an additional two steady states. The upper two states, however, converge and disappear at CY = 0.78, whereas the lower two remain until CY = 0.98. The number of steady states thus varies as 3-5-3-1. At Tb = 700 K in Figure 4C we obtain a situation where the loops formed by each set of steady states do not coexist at any given particle size and the variation in the total
I
1
Ye DECREASE IN R A W S
0
20
40
60
80
% DECREASE IN RADIUS
100
0
20 40 60 X DECREASE IN RADIUS
80
IW
Figure 4. Some loci of quasi-steady state solutions for a particle burning in air.
number of states is then 3-1-3-1. Finally, Figure 4D shows the case of T b = 1100 K. The points marked 7 , 8 , and 9 correspond to the steady states of the 500-pm particle in Figure 4A but the results have been extended to a larger size, 750 pm, in order to show that initially there can be a single ignited steady state. Two different bifurcation points appear in the course of burnoff, one at CY = 0.05 and another a t CY = 0.70 with the result that the variation in the total number of steady states is 1-3-5-3-1. Transient Behavior. A significant result that arises in connection with the structure of the mathematical model is the non-uniqueness of the solution to the distributed problem (45). Numerical solution of eq 45 by the methods already described indicates that a maximum of three different solutions are possible in the boundary layer for a certain range of the parameters. Physically, this means that there is more than one possible set of composition and temperature profiles that result in the same unsteady-state value of the solid temperature. When such multiplicity occurs, we have three different values of dy/dul,,, to choose from in evaluating and 8. Depending on which is chosen, a different temperature history is obtained emanating from a single point y,(O), (~(0). Each of the transient solutions bears a qualitative resemblance to one of the three upper steady states. We will thus refer to them in the course of the discussion at first, second, or third transient depending on which of the three upper steady states they resemble (in order of increasing particle temperature), respectively. It should also be noted that although each transient solution shares the same initial particle temperature, the flux of heat at the surface differs in each case and therefore so do the values of the mole fractions at the surface. Convergence of the numerical scheme to a given transient depends on the character of the guessed temperature profiles. The total number of transient solutions at a given size of the particle varies with the value of the particle temperature given by the initial condition (7). This variation is shown in Figure 5 where we have plotted the value of dT,/dt as a function of T, for two different sized particles, thus illustrating the multivalued character of 4 and 8 in-
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980 249
/ \ \ \ h ; b:l s 0 LL
?
4
5
15
x o
0
W
0
_ _ - + -
20
40
60
80
800
20
40
60
80
00
% D E C R E I S E IN RADIUS
%DECREASE IN RADIUS
z + -4
_.__.---.
15 16 17 IS PARTICLE TEMPERATURE (K) i 1 6 *
Figure 5. Illustration of the multiplicity predicted by eq 45 for two different sized particles burning under the conditions of Figure 4B.
troduced by the multiple solutions of eq 45. For the 4Wpm particle of Figure 5, multiplicity occurs in the range 1610 K I T, I 1750 K. The points where the curve crosses the abscissa correspond to the ignited steady-state solutions represented in Figure 4B a t CY = 0.60. Also shown in Figure 5 is the locus for a 700 pm radius particle. This size corresponds to a 30% decrease in radius of the particle with the steady-state behavior of Figure 4B. At this size there are only two ignited steady-state solutions, yet three transient solutions are possible in the range 1490 K IT, d 1580 K. It should also be noted that the range of multiplicity decreases for larger sized particles, so that by a = 750 pm there is a unique solution. However, past 78% decrease in Figure 4B, when the loop formed by the upper states disappears, there is again a unique transient solution a t any value of the particle temperature. For the other ambient temperatures considered, Figures 4C and 4D, there is multiplicity in the case of the 500-pm particle a t Tb = 700 K for 1575 K IT, I 1735 K but a unique solution in the case of the 750-pm particle at Tb = 1100 K, although no new features of behavior are added. As previously mentioned, the multiple solutions of the distributed equation result in “branching” of the temperature history. That is, given an initial particle temperature and size, we may march along the integration of eq 43 and 44 using one of the possible solutions of (45), each leading to a different temperature history. A typical set of results is shown in Figure 6A for the 400-pm particle in a 900 K ambient environment. Depending on which of the three transients we choose to march along, we obtain the branches labeled c, b, and a, respectively in Figure 6A. A given transient approaches the steady state which more closely resembles it. This is better understood if we examine the magnitude and sign of dT,ldt as shown in Figure 5. For instance, at T, = 1740 K the lower two transients have a negative value of dT,/dt, which explains why the particle temperature initially drops from 1740 K in both cases. For the upper transient the opposite is true. In all three cases burning proceeds next to the steady state approached, except for the upper two branches that shift to the next lower state as the loop formed by the upper steady states disappears. Branches d and e represent the first and second transients that originate out of T, = 1610 K. They approach the same steady states as when the initial particle temperature was 1740 K, except in this case the approach is from below. Branch f originates out of a steady state with no significant CO oxidation and its behavior is similar to some of the results obtained in the one-film model. The behavior of several particles burning at the ambient temperature of 700 K (locus of Figure 4C) is shown in
Figure 6. A select number of temperature histories for the ambient conditions that result in the quasi-steady state structures of Figures 4B and 4C.
Figure 6B. As previously stated a t a, = 500 pm, multiplicity of the transient solution occurs in the range 1575 K I T, I 1735 K. Thus, near the upper bound of this range (1720 K) it is possible to obtain a branch of the first type (labeled a) that quickly drops toward the lowest steady state, ignoring its initial proximity to the upper two. On the other hand, a second transient out of the intermediate state a t T,, = 1620 K generates a branch that burns next to that state up to 43% decrease. At this point, as that state disappears, the CO flame becomes quenched and the particle temperature drops by about 100 K but it eventually increases as it begins to approach the upper state of the second loop. The remaining history does not differ from the one-film case. Another two branches (of the first type) out of T,, = 1750 K and CY = 0.41, CY = 0.42 are presented to illustrate the domain of attraction of the upper state. That branch out of CY = 0.42 also drops by 100 K before its slope changes sign. However, the branch out of CY = 0.41 is able to break free from that state, although it joins the intermediate state temporarily, soon drifting to the lowest state. We also note upon comparing branches a and b that although both correspond to equally sized particles, the one with the higher initial temperature, branch a, results in a greater lifetime than the particle of branch b. Similarly, branch c that represents a particle smaller and at a higher initial temperature than that of branch b, nevertheless possesses a longer lifetime-a result analogous to that illustrated in Figure 3. We must finally emphasize that due to the assumption of a quasi-steady-state boundary layer, and the resulting multiplicity, the actual stability character of the system does not necessarily translate into the present continuous flame model as it did in the case of one-film model. This is why we observe, for instance, that branch b in Figure 6B predicts stable burning along a steady state that our previous stability analysis suggested to be unstable. In fact, this branch may not exist in practice, if unstable. A model that includes transient effects in the boundary layer would therefore be necessary to detect such stability nature. Summary a n d Conclusions The quasi-steady state combustion history of a carbon particle was presented by means of one-film and a continuous flame model. This history indicates that depending on the specific set of ambient conditions, the total number of quasi-steady states may undergo multiple bifurcations during burn-off. Of the cases considered here two of the most general variation patterns were found to be 3-1-3-1 and 1-3-5-3-1. Solution of the developed transient models shows that the dynamic behavior of the particle cannot be generalized
250
Ind. Eng. Chem. Fundam., Vol. 19, No. 3, 1980
from the quasi-steady-state results. According to the initial size and temperature, considerable burning can occur away from any of the quasi-steady states. Furthermore, significant and successive drops in temperature and combustion rate may occur as transitions take place between upper and lower steady states. Anomalous behavior was found to result from some quasi-steady-state structures that produced a shorter lifetime for the larger of two particles burning under the same set of constant ambient conditions. Of particular interest is the fact that due to the assumption of a quasi-steady-state boundary layer, there can be multiple composition and temperature profiles that result for a unique value of the particle temperature specified as an initial condition. A maximum of three such transient solutions are obtainable in a certain range of initial particle temperatures, each generating widely different temperature histories approaching one of the three upper quasi-steady states. Thus it was possible for a particle a t high temperature to become unignited or for one at relatively low temperature to become readily ignited.
Nomenclature A = matrix defined by eq 49 a = particle radius, m B = defined by eq 49 b = value of the radial coordinate at the external edge of the boundary layer, m c = molar concentration of the gas mixture, kg-mol/m3 C, = heat capacity of the carbon particle, J/kg-mol K D = diffusivity, m2/s e = energy flux with respect to stationary coordinates, J/m2 8
E, = internal energy of carbon, J/kg-mol (= H,, enthalpy) E; = activation energy of the ith reaction, kcal/kg-mol f = defined by eq 25 g = defined by eq 26 H , = enthalpy of carbon, J/kg-mol Hi = partial molal enthalpy of the ith species, J/kg-mol k = thermal conductivity of the gas mixture, J / m s K ki = frequency factor for the ith reaction, m/s L = defined by eq 17 M , = molecular weight of carbon, kg/kg-mol Ni = molar flux relative to stationary coordinates, kg-mol/m2 S
q R = radiant energy flux,J/m2 s
Qi = heat of reaction of the ith reaction, J/kg-mol of reactant r = radial coordinate, m R1 = rate of the C + 0, reaction, kg-mol/m2 s R2 = rate of the C + COz reaction, kg-mol/m2 s R = rate of the CO + O2 reaction, kg-mol/m3 s R, = rate of combustion (= 2R1 + R,) kg-mol/m2 s t = time, s T = temperature, K
u = dimensionless distance variable [ = ( 1 - €)/(1- (a/b))] ii = defined by eq 49 xi
= mole fraction of the ith component
y = dimensionless temperature (= T/T,) 7 = defined by eq 49
Creek Letters a = fractional decrease a3 = defined by eq 17
in particle radius (= 1 - & ) / a o )
= stoichiometric coefficient of the jth species in the ith heterogeneous reaction yi = defined by eq 20 y = dimensionless parameter r = defined by eq 40 6 = defined by eq 21 q = activation energy of the ith reaction divided by the universal gas constant, K t = emissivity of carbon 0 = defined by eq 38 p = defined by eq 36 p c = density of carbon, kg/m3 u = Stefan-Boltzmann constant, J/s m2 K4 7 = lifetime of the carbon particle, s 7 , = defined by eq 34 and 35 t = dimensionless distance (= a / r ) vi = stoichiometric coefficient of the ith component in the homogeneous reaction $ = defined by eq 37 x, ii = defined by eq 48 and eq 49 il = defined by eq 45 Subscripts s = particle surface b = external edge of the boundary layer 0 = initial value @ij
Literature Cited Caram, H. S., Amundson, N. R., Ind. f n g . Chem. Fundem., 18, 171 (1977). Field, M. A., Gill, D. W., Morgan, 8. E., Hawksiey. P. G. W., "Combustion of Pulverized Coal", p 192, BCURA, Cheney and Sons, Ltd., Banbury, England, 1987. Michelsen. M. L.. AICM J . , 22. 594 (1976). Mon, E., AmuMon, N. R., Ind. fng. Chsm. Fundem., 18, 182 (1979). Mon. E., Amundson, N. R., Ind. Eng. Chem. Fundem., 17, 313 (1978). Ubhayakar, S. K.. Combust. Flame, 28, 23 (1976). Ubhayakar, S. K., Williams, J. F.. J . Electrochem. Soc.,123, 747 (1976). Villadsen, J. Michelsen, M. L. "Solution of Differential Equation Models by Polynominai Approximation", p 305, hentice-Hall, Englewood Cliffs, N.J., 1978. Vulis, L. A. "Thermal Regimes of Combustion", McGraw-Hili, New York, N.Y., 1961.
Received f o r review June 26, 1979 Accepted February 11, 1980 The authors are grateful to the University of Houston for providing financial support during the course of this research. For the past year this work was funded by Department of Energy Grant ET-78-G-01-3020.
