Model of Coal Pyrolysis. 2. Quantitative Formulation and Results

122

Ind. Eng. Chem. Fundam. 1981, 20, 122-132

Stein, S.; Golden, D.; Benson, S.; Shaw, R., paper presented at the Coal Chemistry Workshop, Stanford Research Institute, Aug 1976. Suuberg, E. M. sc.D. Thesis, Massachusetts Institute of Technology, 1977. Suuberg, E. M.; Peters, W. A.; Howard, J. B. Ind. Eng. Chem. Process Des. D ~ v1978, . 17. 37-46. Szwarc, M. J. Chem. fhys. 1048, 16, 128-136. Tingey, G. L.; Morrey, J. R. Batteiie Energy Program Report, Battelle, Pacific Northwest Laboratories, 1972. van Krevelen, D. W. f o @ m r 1975, 76, 615-620. Walling, C.; Lepley, A. R. Inf. J. Chem. Kinet. 1971, 3 , 97-104.

Wachowska, H. Fuel1979, 58, 00-108. Whltehurst, D. D. ACS Synp. Ser. 1978, No. 71. Williams, G. H. "Homolytic Aromatic Substhution”, P w ~ m W York, 1960. Yokohama. S.; BOdiiy, D. M.; Wiscw. W. H. FUel1979, 58, 162-170. Zavksas, A. A.; Meliklan, A. A. J. Am. Chem. Soc. 1975, 97. 2757-2763.

Received for review May 27, 1980 Accepted January 5,1981

Model of Coal Pyrolysis. 2. Quantitative Formulation and Results George R. Gavalas,’ Ravl Jaln, and Paul How-Kel Cheong Division of Chemistry and Chemical Engineering, California Institute of Technology, Pasadena, Callfomla 91 125

The reaction rates are expressed in terms of concentrations of reactive conflgurations which are calculated by randomly distributing the functional groups among the a carbons. The probability of formation of a free unit (tar molecule) following a bridge dissociation is ca!culated by a similar random placement technique. The rates of change of the state variables and the rates of product formation are expressed in termq of the various reaction rates and the rate of loss of free units. Computer simulations are carried out for an hvc bituminous coal and a subbituminous coal and the results are compared with limited experimental data. A sensitivity analysis is carried out to study the importance of various structural and kinetic parameters relative to the yield of various products.

Introduction In part 1, we proposed a representation of coal as a collection of functional groups and suggested a set of free radical reactions to describe chemical change during pyrolysis. The mechanism of formation of gases, tar,and char was explained in terms of these reactions. In this part, we formulate a mathematical model for quantitative description of the kinetics of coal pyrolysis. Following the general formulation, simulation results will be presented for the pyrolysis of two coals. Compared to modeling other free radical reaction systems (e.g., pyrolysis of hydrocarbons), the modeling of coal pyrolysis involves three additional complications: (i) rate expressions for the reactions are complicated because the reactants are not individual molecules but combinations of functional groups on the coal matrix (Table 11, part 1); (ii) the concentrations of functional groups change due to reactions and due to loss with departing volatile products; (iii) complications arise due to the fact that reactions occur in the condensed phase; for example, products of dissociation may undergo recombination before escaping the condensed phase. In the following sections we present methods developed to deal with the above difficulties. Then we proceed to formulate the differential equations and obtain numerical results for two coals. Reactive Configurations and Reaction Rates In part 1we defined a set of state variables representing the concentrations of various functional groups and certain reIated quantities and developed a procedure for calculating the initial values of the state variables from analytical data. In this section we first illustrate the difficulty in relating the state variables to reaction rates and then proceed to develop techniques to overcome this difficulty. Consider the methyl group in the two molecules PhCH2CH3,Ph-CHCH,. Only the first molecule can dissociate to produce a CH3-radical (reaction 2 in part 1). The 0196-4313/81/1020-0122$01.25/0

carbon-carbon bond in the second molecule cannot dissociate because of the unpaired electron on the a carbon. Dissociation is also impossible when the a carbon po&sesses a double bond. To calculate the rate of production of methyl radicals by dissociation, we need to determine the fraction of methyl groups which are attached to a carbons which have no radicals or double bonds. The total concentration Yz of methyl groups cannot be used directly in the rate expression. Similarly, to compute the rate of double bond formation, reaction 7, the concentration of a carbons with a methyl group and an unpaired electron has to he computed. Yet another example is bridge dissociation, reaction 5, where we need the number of combinations of two a carbms, each of which participates in an ethylene bridge but has no radical or double bond. Evidently, a scheme is needed to distribute functional groups, including radicals and double bonds, to a carbons in order to compute the concentrations of various combinations of functional groups. These combinations, termed reactive configurations, will be the “molecules” participating in reactions. The scheme developed here is based on random distribution. In the absence of information on any naturally preferred arrangement of functional groups, random distribution is a reasonable approach. Moreover, the initial distribution derived from the coalification process will tend to become randomized by successive dissociation, recombination, and addition reactions. Under the scheme proposed, the functional groups are first classified into five categories and then randomly distributed in groups of two’s or three’s to a carbons. The five categories are represented by the following auxiliary variables: MR,concentration of a radicals; MD,concentration of double bonds-double bonds on chains and double bonds on bridges; Mp,concentration of methylene bridges; MH,concentration of hydroaromatic structures, including those bearing radicals; and MG,sum of concentrations of a hydrogens, methyl and ethyl groups (including 0 1981 American Chemical Societv

Ind. Eng. Chem. Fundam., Vol. 20, No. 2, 1981

123

their /3 radicals), and ethylene bridges. The variables MR, carrying one radical, one hydrogen, and one methyl group. ..., MG can be directly expressed in terms of the state The concentrations of these configurations are p?p2 K m and plp2K m . The rates of the reactions are kg?p2 Km variables. The distribution of the auxiliary variables among a and k7pG2KRGG. Consider next the reaction carbons is subject to the following five constraints, but I I + Ph'. otherwise completely random: (i) at most one double bond Ph-C-Ph' Ph-C. per a carbon; (ii) a t most one radical per a carbon; (iii) a I I radical and a double bond cannot coexist on the same a carbon; (iv) at most one methylene bridge per a carbon; The reactive configuration here consists of an a carbon (v) at most one hydroaromatic structure per a carbon. The carrying a methylene bridge but no radicals or double first three are obvioua chemical constraints, while the last bonds. The rate is given by k6[(1 - p5)KHp~+ (1 two are introduced to simplify the mathematics by rep5)'KpGGI. The rates of all other reactions in which the ducing the number of possible combinations. As a conreactive configuration involves a single a carbon can be sequence of the fourth constraint, an a carbon can be calculated in a similar fashion. connected to at most two aromatic carbons. Such asRates of those reactions in which two a carbons are sumptions are unlikely to affect the results in any signifinvolved cannot be computed in terms of the K's alone. icant way. To illustrate this, let us consider reaction 5 Since each a carbon carries three substituents, or two if it carries a double bond, the possible combinations of the five auxiliary variables are I: KDH, KDP,K,; II: K ~ G , K R H KR~G, ~ , K m ; I11 K H ~ G KHm; , nT: KpGG, K W Each K represents the concentration of a carbons carrying the The reactive configuration for this reaction consists of a substituents indicated by the subscripts. The subscripts pair of a carbons, each carrying an ethylene bridge but no have the same meaning as in the M's; for example, KRHp radicals or double bonds. Therefore, the concentration of represents a carbons carrying a radical, a hydroaromatic this configuration cannot be expressed in terms of the ICs structure, and a methylene bridge. The order of the which represent single a carbons only. To compute the subscripts is immaterial, i.e., KRHp,KRpH, KmR, etc. denote concentration of this pair, we first divide the a carbons the same quantity. The eleven K's satisfy the following carrying ethylene bridges into classes Z and X, with and five balances without a radical or a double bond, respectively. These two classes induce the following combinations: (i) comKDH+ KDP+ KDG= MD (1) binations in which both a carbons are free of radicals and KRHG+ KRHP+ KRPG+ KRGG=MR (2) double bonds: L;, (ii) combinations in which only one a carbon is free of radicals and double bonds: Lxz; and KDH+ KRHG+ KRHP+ KHPG+ KHGG= MH (3) (iii) combinations in which both a carbons contain a radical KDP+ KRHP+ KRPG+ KHPG+ KPGG= MP (4) or a double bond Lzz. It should be noted that an a carbon may carry more than one ethylene bridge; therefore, it may KDG+ KRHG+ KRPG+ ~ K R G+GKHPG+ ~ K H G+G be involved in more than one combination. If XI, X2, X3 ~ K P G+G~ K G G=GMG (5) represent the concentrations of X-type a carbons carrying The sum of all the K's equals the total concentration of 1 , 2 , 3 ethylene bridges and Z1and Z2represent the cona carbons, Y13 This leads to the overall balance already centrations of Z-type a carbons carrying 1 and 2 ethylene satisfied by the M's bridges, the concentrations X1, X2, X3 and Z,, Z2can be easily computed in terms of the K's. The L's are to be MR + 2MD Mp + MH + MG = 3Y13 (6) calculated by a combinatorial placement scheme analogous The calculation of the K's subject to the constraints to the one used for the K's. In this case the constraints (1)-(5) is a combinatorial placement problem which is are solved in Appendix A. A relationship is thus established 2Lxx + Lxz = X I + 2x2 + 3x3 (7) between the concentrations of reactive configurations,ICs, 2Lzz + Lxz = z1 + 222 (8) which represent particular combinations of functional groups undergoing reactions, and the state variables, which The solution of the placement problem leads to explicit represent the total concentrations of the functional groups. expressions for the L's, (Jain, 1979; p 123). The ICs in categories I and I1 are zero before reaction due With Lxx being the concentration of the reactive conto the absence of radicals and double bonds. figurations for reaction 5, the corresponding rate is given To distinguish between the various components of by k5 (1 - ml - m,J2L*. The quantities ml and m2 category G, we define a set of probabilities p1through p6 characterize X and are defined in Appendix B. Similarly, in Appendix A. For example, p1is the probability that G for reaction 11 is an a hydrogen, p 2 the probability that it is a methyl H group, etc. Also defined in Appendix A are probabilities II PH1 and P H 2 distinguishing between hydroaromatic strucPh-6-C-Ph' PhC=C-Ph' + H. tures with and without a /3 radical. I I I t We are now in a position to calculate the rates of the the reactive configuration is Lxz, and its rate is given by reactions considered earlier. The calculation is illustrated kll (clc2Lxz).The quantities c1 and cz are also defined in by reactions 2 and 7 of Table 11, part 1. Appendix B. Ph-CH,CH, +Ph-CH, + CH,. With the help of the K's and the L's, we can calculate the rates of all reactions of the carbon-hydrogen structure Ph-CHCH, +Ph-CH=CH, + H. listed in Table I1 of part 1. A detailed listing of all the rates can be found in Jain (1979), p 139. I t is important The reactive configuration in the first reaction is an a carbon carrying two hydrogens and one methyl group, to note that despite their complex appearance, the reaction while in the second reaction it consists of an a carbon rates are of the standard form rj = kjfj(Yi), with the

-

+

-

124


functions f .( Yi) given by well defined transformations involved in the random placement of groups on a carbons. The rate of the phenolic condensation reaction is assumed to have a second-order dependence on phenolic concentration, while the rates of reactions of carboxylic and carbonyl groups are assumed to be first order. The kinetic parameters of all reactions were listed in part 1, Table 111. Rate of Generation of Free Units In this section we formulate the differential equations governing the rate of change of the state variables during pyrolysis. Recalling that Y15is the volume of the condensed phase, the equations are of the form (d/dt)(Y,Yi,) = Y15 (change due to reactions - loss with free units)

randomly among N1and N2 The distribution is carried out using the technique adopted earlier for the ICs and leads to the following combinations. Nl: N b , N b b , Nbbb; N $ Nb-, Nbbnn, Nbbbb.. For example, Nbbn represents the concentration of nuclei which bear three CY carbons, two of which carry at least one bridge with the third carrying no bridge. Similarly, Nbbnn stands for concentration of nuclei carrying four CY carbons, two with bridges and two without. It should be noted that N,, and N,,, are not acceptable combinations because each nucleus must bear at least one connection with the coal matrix. We now have

for i = 1-14. The first term on the right side can be readily computed in terms of the various reaction rates. The second term is zero for the state variables representing bridges, Le., Y4,Ye,and Yll, because as discussed below only single-unit fragments are considered. To obtain the second term for the other variables, two quantities need to be calculated: the rate of escape of free units and the concentration of functional groups on free units. The dissociation of a bridge can lead to the formation of a volatile fragment which could be a monomer, dimer, or perhaps trimer, i.e., consist of one, two, or three structural units. These products are collectively described as tar. The probability of formation of dimers and trimers has been shown by Cheong (1976) to be significantly lower than that of monomers. Moreover, dimers and trimers have significantly lower volatility and diffusion coefficients than monomers. On this account, we will consider dimers and trimers as remaining in the condensed phase and attribute tar production solely to the generation of single unit fragments. A unit becomes free following the dissociation of a bridge possessed by a precursor unit. A precursor unit is one that carries exactly one “breakable” bridge. To determine the rate at which bridge dissociation results in free units, we need the probabilities p g ,plo,pll that a unit becomes free following the dissociation of an ethylene, a methylene, or an ether bridge, respectively. To determine the probabilities p9 and plowe must at first distribute YEa carbon connections to YI4 nuclei. This distribution is complicated by the fact that different nuclei may carry different numbers of a carbons. Although Y8/Y14 is typically between 2 and 5 for subbituminous and bituminous coals, a given coal will contain a mixture of nuclei carrying 1, 2, 3 etc. a carbons. Accounting for all different possibilities involves unnecessary complexity. To simplify the situation, we assume that a nucleus carries 3 or 4 substituent a carbons. The number 3 or 4 represents a degree of substitution typical of the coals simulated in later sections. If N 1 and N 2 represent the concentrations of nuclei carrying 3 and 4 substituents a carbons N1 + N2 = YI4 3N1 + 4N2 = Ys N1 = 4Y1,- Ys N2 = Ys - 3Y14 The substituent a carbons can also be divided into two parts: (i) those which involve (Y carbons participating in at least one bridge, concentration Yab, and (ii) those which involve a carbons participating in no bridge, concentration Yan. These two types of N carbons must be distributed

straints follows the principles of Appendix A and is given in Jain (1979), p 129. Having determined the N’s we need two probabilities, p7and pa,that a given bridge is a breakable ethylene or a breakable methylene bridge, respectively. The computation of these two probabilities as functions of the Rs is given in Jain (1979), p 130. We now have the following expressions for pg,plo

Nbnn

+ 2 N b b n + 3 N b b b + Nbnnn + 2Nbbnn + 3Nbbbn + 4Nbbbb

=

Y8b

+ N b b n + 3Nbnnn + 2Nbbnn + Nbbbn = The determination of the N’s subject to these two con2Nbnn

Finally, the probability pll that a unit is free of ether linkages can be calculated by distributing the two oxygen groups, Ylo and Yll, among Y14units. This distribution depends on the ratio (Ylo + Y11)/Y14which represents the oxygen content of the coal. Depending on the value of this ratio, three cases of distribution have been considered in Jain (1979), p 131. In each case the calculation of the fraction of units lacking ether bonds is similar to that of the N’s above. The three probabilities, pg,pl0,pll, along with the rates of bridge dissociation, are sufficient for calculating the rate of generation of free units. Of the free units generated, some will reattach to sites on the coal matrix and will thus be unable to escape the coal particle. Reattachment can occur by radical-radical recombination or by radical addition to aromatic rings. The probability of reattachment depends on the rates of these reactions relative to the rate of diffusion. As mentioned in the qualitative discussion of part 1, the transport of free units in a coal particle involves transport through the condensed phase (the coal matrix including the micropores) and through the pore phase (the transitional pores and the macropores). The transport through the condensed phase takes place by slow activated diffusion. A simple mathematical analysis of activated transport and simultaneous reattachment of free clusters in the coal matrix is presented in Jain (1979), p 156. However, the quantitative application of such analysis is subject to considerable uncertainty regarding the rates of activated diffusion and in any case it would complicate the mathematics significantly by introducing spatial variation of the concentrations. We have adopted a simple empirical approach to simulate reattachment during transport through the condensed phase. It is assumed that only those free units which are generated in a fraction X of the condensed phase near the surface of the transitional pores and the macropores escape to the pore space. The bridge dissociation


reactions in the remaining fraction 1 - X do not lead to formation of free units. The effect of this assumption is to multiply the rate of bridge dissociation reactions by a factor X < 1. The fraction X is treated as an adjustable parameter. Once in the pore space, free units can still reattach to the pore surface. While not as likely as the reattachment in the condensed phase, this pore phase reattachment can also reduce the escape of units as evidenced by the role of particle size and external pressure on tar yield (Anthony and Howard, 1976; Suuberg, 1978). An analysis of transport limitation in the pore space was presented by Gavalas and Wilks (1980). Again, incorporation of transport processes would unduly complicate the chemical model. Consequently, in the present model we assume that once in the pore space, free units suffer no additional reattachment. This assumption is reasonably accurate at low pressures and small particle sizes. We are now ready to calculate the rate of escape of free units. Consider reaction 5, the dissociation of ethylene bridges

I I I I

Ph-C-C-Ph'

-

Ph-C.

I + .C-Ph' I I I

Following dissociation, either one or both of the units could be free, hence (rate of escape of units by reaction 5 ) = xr5[P9Pll(l - PS) + PCI2Pll(l - P11) + P$ll(l - PS) + Pg2pll(l - P11) + ~PSz~l121 = 2P9Pllxr5 The first four terms on the right-hand side account for one of the two units becoming free, while the last term represents the probability that both units become free. Similarly for other bridge dissociation reactions reaction 6 18 19 36, 37, 38

rate of escape of units 2P 10 P 1 1 x r 6 2 p 9 P 1 1 Xr1, 2PtoP11Xr19 2PZp,JJ,,X(r, + r37 +

r38)

(overall rate of escape of units due to bridge dissociation) = x[2P9Pll(r5 + r8) + 2Plfill(r6 + r19 + r36 + r37 + r38)l (9) Rates of Change of Functional Groups The escaping free units carry all the functional groups attached on their peripheral aromatic carbons. The rate of loss of a functional group with escapting units can be calculated by considering each bridge dissociation reaction separately. For each reaction the rate of escape of free units is multiplied by a probability function representing the concentration of a given functional group on the generated unit. To determine such a probability function, we first calculate the number of a carbons carried by units freed as a result of a reaction. Free units generated by reactions 5, 18, 37, and 38 carry either 3 or 4 a carbons, while half of the units produced by reactions 6,19, and 36 carry only two or three a carbons. Next, the probability of the presence of a functional group is computed in terms of the K s for each a carbon separately. Some functional groups like a hydrogen, methyl groups, and ethyl groups, can occupy one, two, or three substituent positions on an a carbon, while others like hydroaromatic structures and a radicals can occupy only one position. All possible combinations are accounted for while calculating the probability.

125

This detailed calculation leads to complicated expressions for the overall rates of loss of functional groups with escaping free units, as illustrated for the state variable Yl, in Table I. Simpler approaches may be adopted if care is taken to ensure that the state variables continue to satisfy the following balance (discussed in part 1) at all times during integration. Y1 + Yz Y3 + Y4 Y5 + (10) = 4Y13 The rate of loss of phenolic groups is proportional to the state variable Ylo and is computed directly. The Differential Equations (i) For small radicals, H-, CH3., C2H5. (rate of change) = (net production by reactions)

+

+

For example

d -[H.] dt

= Y15[rl

+ r7 + r9 + rlo + rll + r15 - r23 - r26al3r30 - a13r33- r361

where al = 1 - p4- p5 - p6. (ii) For the total moles of units (rate of change) = -(rate of escape of free units)

(iii) For the total moles of ethylene bridges (rate of change) = (net production by reactions) d z(y4y15) = 2Ydr22 - Xr5 - rll - r12- r13 - r14 - XrlJ (iv) For the total moles of double bond bridges and ether bridges (rate of change) = (net production by reactions) d -(y6y15) = 2Y15(rll + r12+ r13 + r13 dt d = zY1d39 dt (v) For the state variables Yl, Y2, Y3, Y5, Y8, Ys, Ylo, Y12, -(y11y15)

y13

(rate of change) = (net production by reactions) (rate of loss with escaping units) For example

+

+ +

Y15[a13(3r31 + 2~32)+ a2(3r3, + 2r35) 3r3, 2r38r1 - rll - r15 - r23 - r24 - r25 - ~ 2 -9 (a3+ 2a4 2p13)x (r30 + r31 + r32) - (a5 + 2 ~ ~ ~+ )r34( +r d ~1~Y15 (rate of loss of Yl with escaping units)

where a2 = 1 - p3 - p5

- p6;

a3 = 3p1pd2; a4= 3Pl2p&a5 = 2PlPd; Pd = - P1 - P4 - P5 - P6 The rate of change of the volume, dY15/dt, which appears a t the left-hand side of all differential equations, is yet to be evaluated. Since very little consideration has been given so far to physical changes brought about by chemical reactions, it will be assumed that the volume of the condensed phase changes solely due to the loss of units. In other words, the concentration of units in the condensed phase remains constant, dY14/dt = 0. This assumption

126

Ind. Eng. Chem. Fundam., Vol. 20, No. 2. 1981

is consistent with certain previously imposed constraints, for example, that the number of LY carbons per unit is either 3 or 4. We now have d ,(y14y15)

=

dY15 Y ' 4 d t

Using the differential equation for the units we obtain dY15 -- dt y15

-2Xpii-[pg(r5 y14 "

+ ria) + p10(r6 + r19 + r36 + r37 + Qe)I

DD

1-

+

e "

a,

0

4

+

.

2 "

DD e

h

Y

3B L

+ cL

+ ;R Y

v)

+ "

a,

0

4

3. h

Y

w 0

With the help of the last relationship, the differential equations can be reduced to the standard form in terms of the derivatives dYi/dt. Overall, a set of 18 simultaneous differential equations is obtained for various condensed phase species and the gas phase radicals. The solution of these equations also provides the instantaneous production of the gas phase species H20, H2, CHI, CzH4,C2HB.Although the above differential equations are in principle similar to those for any multi-species, multi-reaction system, the functions involved at the right sides are very complex because of the various statistical factors present and the lengthy expressions describing the loss of functional groups with escaping units. To reduce the numerical difficulties associated with the stiff nature of the system, the steady-state approximation can be applied to the gas-phase radicals H., CH3., CzH5-. Eliminating the concentrations of these three radicals, we obtain a system of 15 simultaneous differential equations in a form suitable for numerical integration. Computer Simulation The differential equations were integrated numerically on the computer to simulate pyrolysis of two bituminous coals. After a brief discussion of the computer code, the simulation results for each coal will be presented and compared with limited experimental data. The sensitivity of the results to the values of the adjustable parameters will also be discussed. The differential equations formulated above constitute a stiff set as a result of large differences in the time constants of the various elementary reactions. A package of FORTRAN subroutines called EPISODE, developed at the Lawrence Livermore Laboratory (Hindmarsh and Byrne, 1975,available from Argonne National Laboratory) for stiff systes, was used for the numerical integration. This package is an improved version of the earlier package GEAR (Hindmarsh, 1974),based on Gear's method (Gear, 1971). The main difficulty encountered in the integration was in the calculation of the Jacobian matrix. Due to the complexity of the functions involved, the derivatives could not be computed analytically and had to be approximated by difference quotients. The increments in Y,,automatically chosen by the program to calculate the difference quotients, were sometimes so large that the Y's were unable to satisfy such basic constraints as eq 10, which is vital to the various statistical distribution schemes. The continual updating of the Jacobian matrix accounted for a major fraction of the consumed computer time. Results of Simulation The parameters needed for the numerical calculations included the initial values for the state variables and the kinetic parameters for the reactions. All calculations presented below involve isothermal conditions. As presented in part 1, the calculation of the initial values of the state variables is based on elemental com-


Table 11. Analytical Data for Two Coals subbituhvc bituminous minous elemental compn, (PSOC (Hamilton, % 212) W. Kentucky) mineral matter/ash 3.0 7.9 77.0 73.5 carbon hydrogen 5.0 5.1 oxygen 12.55 8.7 1.55 1.5 nitrogen sulfur 0.90 3.3 1.3 1.3 density, kg/L aromaticity 0.70 0.69 hydrogen types, % 32 28.5 30 32 32 35 ~henolic 6 4.5 metky lene/eth ylene 0.27 0.25 bridges concn of units, 5.83 5.10 g-mol/L

position, content of various types of carbon and hydrogen, and certain structural parameters whose values have to be assumed. The information available on the two coals under consideration is listed in Table 11. The elemental composition of the Hamilton coal was determined by commercial laboratories, while that of PSOC-212 was taken from Solomon (1977). The helium density was determined from a graph of density-carbon content given by Gan et al. (1972). The aromaticity fa was estimated from a graph presented by Suuberg (1977) to be approximately 0.7 for both coals. A few trial runs were made with f a between 0.68 and 0.72, and it was found that 0.69 and 0.70 gave the best fit with the experimental data for the Hamilton coal and PSOC-212, respectively. The content of various hydrogen types in the Hamilton coal was estimated from the NMR spectra of pyrolysates reported by Gavalas and Oka (1978). Some adjustments were made in the fractions of (Y and 0 hydrogen to account for structural differences between the pyrolysates and the parent coal. Similarly, values slightly modified from those reported by Solomon (1979) are listed for hydrogen types in PSOC-212. The fraction of phenolic hydrogen was estimated from the amount of water produced by pyrolysis. The values of two of the structural parameters, the ratio of methylene to ethylene bridges, and the concentration of units had to be assumed. In selecting a value for the ratio of bridges we note that methylene bridges possess higher dissociation energy; therefore increasing the ratio has an adverse effect on tar production. Similarly, the concentration of units, Y14,affects tar production by determining the number of bridges per unit, which in turn determines the probabilities p g and plo.The values of Y14 are subject to some other constraints discussed in part 1. The values assigned to the two parameters are listed in Table 11. The structural information in Table I leads to the initial values of state variables listed in Table 111. The two coals are seen to be very different in their content of (Y hydrogens, ethylene bridges, and phenolic oxygen. State variables missing from Table I11 were initially zero. The parameter X was assumed to be 0.15 for the Hamilton coal and 0.2 for the PSOC-212. The parameter values for the computer simulation were chosen as follows. The kinetic parameters were estimated from the literature or chosen more or less arbitrarily. The structural parameters such as YI4(O), X , initial ratio of methylene-to-ethylene bridges were adjusted to provide reasonable agreement with experimental data. The ad-

127

Table 111. Initial Values of the State Variableso for Two Coals concn, gmol/L hvc bituminous subbituminous (Hamilton, (PSOC212) W.Kentucky) 19.5 21.2 4.04 4.51 1.56 1.74 21.4 17.6 0.501 0.449 19.5 18.1 3.90 3.00 16.6 15.9 volume, cm3/g 0.77 0.77

Yl y* y3 y 4

y, Y 8

YlO y,3 y,5

Unlisted variables are zero. Table IV. Values of Kinetic Parameters values used in best estimated simulation values reaction 1 2

3 4 5 6

7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39

logA

E , kcal/ g-mol

15.5 15.3 14.9 14.9 14.4 14.4 12.8 12.8 12.8 12.8 12.0 12.0 12.0 12.0 12.3 12.3 12.6 13.6 13.6 12.4 12.4 8.0 11.1 8.0 8.0 11.1 8.0 8.0 8.0 11.0 7.0 7.0 11.0 7.0 7.0 11.0 7.0 7.0 8.5

84.0 65.0 63.0 63.0 52.0 68.5 50.0 48.0 50.0 50.0 38.0 37.0 37.0 38.0 38.0 38.0 35.0 35.0 35.0 35.0 36.0 10.0 10.0 10.0 10.0 10.0 10.0 10.0 20.0 15.0 15.0 15.0 10.0 10.0 10.0 10.0 10.0 10.0 35.0

loaA 14.9 15.3 15.4 15.4 13.9 14.3 15.1 14.4 15.1 15.1 12.8 12.1 12.1 12.1 12.8 12.1 12.1 13.0 13.0 12.1 14.2

E , kcal/ g-mol 81.3 68.4 65.0 65.0 50.4 76.7 55.7 49.0 56.0 55.7 54.2 47.5 47.4 47.4 34.3 23.6 20.8 7.0 57.6 20.9 13.5

10.0 7.5 7.0 10.3 7.8 7.3

2.3 8.0 8.9 9.7 10.8 13.4

10.4 7.8 7.3 10.4 7.8 7.3 10.4 7.8 7.3

2.0 7.0 9.0 2.0 7.0 9.0 2.0 7.0 9.0

justment of these parameters was not systematic; therefore the set of kinetic and structural parameters used in the simulations is by no means optimal. Table IV lists the parameter values used in the simulation as well as the best estimated values. The two sets are different because many of the values used in the simulation were assigned rather arbitrarily before it was realized that they could be estimated by group additivity methods. The differences are largest for the radical dissociation reactions 7-20 and the addition reactions of the

128


-

5 0 7 I +

I

Exp. H 2 0 , S I O T

25

$40-

0

z

w

moot W5lOT

T600"C

I

_

W ' wt T

+ x

7

LOSS

Tor Yield Exp W i LOSS 5 I 0 T Exp Tar 510T

510'C

HZ

IO

20

3b Time

is)

4b

sb

60

Figure 3. Simulated and experimental yields of water vapor and hydrogen from the pyrolysis of a bituminous coal (Hamilton).

o

i//

~

Time

Figure 2. Simulated and experimental yields of hydrocarbon gases from the pyrolysis of a bituminous coal (Hamilton).

hydrogen atom. The differences are smallest for the primary dissociations 1-6 and hydrogen abstractions 23-29. In future applications it is recommended that the estimated values be used, at least as a starting point. In this case the structural parameters will also have to be adjusted, since kinetic and structural parameters jointly determine the distribution of products. The estimated parameters given in Table IV were obtained by adjusting the values given in Table I11 of part 1 as follows. (i) The activation energies of reactions 1, 2, 3, 4, 6, are reduced by 4 kcal/g-mol to account for the additional stabilization of the a carbon radical (benzyl radical) due to multiring nuclei and activating phenolic substituents. The activation energy of reaction 5 was likewise increased by 6 to account for the production of two benzyl radicals. The activation energies of reactions 7-14 were increased by 4 to account for the increased stability of the benzyl radical due to multiring nuclei and phenolic substituents. (ii) The parameters of Table 111, part 1, correspond to gas-phase kinetics. It is assumed that the same parameters apply to the condensed phase for all reactions except 22 and 29, which are diffusion limited. The parameters of these two reactions cannot be estimated a t present. The simulation results are shown in Figures 1, 2, and 3 for the Hamilton coal, and in Figures 4, 5, and 6 for PSOC-212. We will discuss the results for Hamilton coal in the context of the mechanism of wrolvsis detailed in " part 1. The results for PSOC-212 can be explained in a similar fashion. Consider the various results for 510 "C. The rates of weight loss and tar production are very high initially but diminish relatively rapidly as pyrolysis progresses. In the first second, the weight loss is almost exclusively due to tar production. The rates of production of hydrocarbon gases and hydrogen change very slowly, resulting in a nearly linear cumulative production curve. Ultimately, the cumulative amounts reach asymptotic

(5)

Figure 4. Simulated weight loss and tar yield from the pyrolysis of a subbituminous coal (PSOC 212). 3.01

-8 -

2.0Ctiq

565'C

__

Time ( s i

Figure 6. Simulated yield of water vapor and hydrogen from the pyrolysis of a subbituminous coal (PSOC212).

values, but this requires more than 60 s. Most of the water vapor is produced during the first few seconds of reaction. Tar is the major component of the volatiles, followed by water vapor, methane, ethane, hydrogen, and ethylene.


Table V. Sensitivity of Tar Yield to Structural and Kinetic Parameters perturbed parameter base value value

129

~

X

bridge-to-cluster ratio methylene-to-ethylene bridge ratio E4 E39

E,, E,, E, Table VI.

0.15 2.16 0.25 52 35 38, 37, 37 65, 63, 63

% perturbation

0.10 2.249 0.30 51 36 39, 38, 38 66, 64, 64

-33.3 4.07 20.0 -1.9 2.9 2.7 1.6

% change in yield

-11.6 -9.9 -8.1 +16.2 +5.2 +5.2 t 0.72

35 24 3 -41 -853 179 193 45

Sensitivity of Hydrocarbon Gases to Structural and Kinetic Parameters %

parameter

X bridge-to-cluster ratio methylene-to-ethylene bridge ratio

E4 E, Eia, E113 E12 E , , E33 E ,

perturbation -33.3 4.07 20.0 -1.9 2.9 2.7 1.6

These trends regarding the time variation of the rates of tar and gases are in agreement with a number of previous experimental studies. Figure 7 showing the evolution of bridges and double bonds is helpful to a mechanistic interpretation of the results. The rate of production of free units (tar) is the product of the bridge dissociation rate and the probability that a unit becomes free following dissociation. The first factor is proportional to the concentration of "breakable" bridges, while the second is determined by the total concentration of bridges. With the concentration of total bridges slightly increasing with time and that of breakable bridges decreasing with time, both factors determining the rate of tar production decrease as pyrolysis progresses. The difference between total and breakable bridges consists of the permanent bridges which include ether bridges and bridges with double bonds. Since the concentration of ether bridges in this case does not exceed 1.5 mol/L (the initial phenolic concentration was taken as 3 mol/L), the largest part of the permanent briges is due to double bond formation. As a result of the increase in total bridges and bridges with double bonds, tar formation essentially ceases after 20 s of pyrolysis time. The rate of hydrocarbon gas production initially lags behind that of tar because the activation energy for chain dissociation is higher than for bridge dissociation. The rate at which chains are lost with free clusters exceeds the rate of formation of small radicals by reactions 1, 2, and 3. However, as reaction proceeds, the production of tar slows down, while the production of small radicals continues by indirect elimination, reactions 7-13 and 15-17. These radicals subsequently abstract hydrogen atoms to form gases. The indirect elimination helps to maintain gas production for a long time, although eventually all volatile formation halts due to the depletion of chains and the formation of double bonds on the remaining chains. At 600 "C, volatile production proceeds much more rapidly, with tar and gas production ceasing after about 2 s and 10 s, respectively. The total amount of tar produced is almost the same as that at 510 "C. Clearly, temperature affects the rate rather than the total amount of tar produced. The latter depends mainly on structural limitations such as the bridge to unit ratio. The amount of water vapor formed at 600 "C is slightly less than at 510 "C,due to an increased loss of phenolic groups with escaping tar,before water forming condensation can occur. The other gases show a substantial increase at 600 "C. Since chain dissociation reactions have an activation energy

% change in alkanes -3.27 + 3.5 +0.7 +9.7 -1.9 -27.6 -12.6

% change in alkenes + 20.6 t i9.8 +5.4 + 10.3 t 39.4 -22.6 -4.0 I

s m e s

-6 2 486 27 -542 1358 - 837 -250

a0

t

T:SIO"C

> e

Saenes

9.8 86 3.5 -510 -66 -1022 -788

Total

10

ooy

Ib

A 3b 4b Time ( s )

sb

'

60

Figure 7. Simulated evolution of the concentration of bridges during the pyrolysis of a bituminous coal (Hamilton).

higher than that of bridge dissociation, higher temperature favors gas formation relative to tar formation. A certain fraction of chains which at low temperatures are lost with escaping free units, at higher temperatures dissociate to form small radicals and subsequently gases. The shift in selectivity from tar to gases with increasing temperature has been amply established experimentally. The results of the calculations for the PSOC-212 coal are qualitatively similar to those of the Hamilton coal. Quantitative differences arise mainly due to the higher phenolic content of PSOC-212 resulting in more ether linkages and decreased tar production. Comparisons of the calculated results with limited experimental data for the Hamilton coal are shown in Figures 1, 2, and 3. The experimental results were obtained by pyrolyzing a "captive" sample of coal particles, placed in the folds of an electrically heated wire mesh inside a reactor. The experimental technique has been described previously, Jain (1979) and Gavalas and W h (1980). The agreement in Figure 1 is reasonable, considering that complete tar recovery has proved very difficult. In Figure 2, the predicted yields of hydrocarbon gases are in very good agreement with the experimental yields. However, the comparison for water vapor in Figure 3 is not so favorable. Reliable determination of water vapor by gas chromatographic methods was difficult due to inaccurate integration of the broad and skewed peak. Parameter Sensitivity Studies Given the uncertainty in the values of some of the structural and kinetic parameters, it is important to determine the sensitivity of the predicted results to these

130


parameters. Sensitivity calculations were carried out for seven parameters, three structural and four kinetic, relative to the pyrolysis of the Hamilton coal at 510 " C for 30 s. The results are shown in Tables V and VI, where the sensitivity S is defined as S = 100 X (% change in predicted quantity/ % change in parameter). 1. A 33% reduction in X reduced the total amount of tar by almost 1270, that of alkanes by 3%, but increased the total amount of ethylene by 21 5%. The reduction in tar is expected because of the decrease in the rate of bridge dissociation and the resulting reduction in free unit formation. The change in X was found to have a greater effect on the tar-time variation than on the ultimate yield of tar. After 1 s, the difference in the cumulative yield was 25%, after 10 s 17%,while after 60 s it was only 9%. The yield of alkanes depends on the production of small radicals, which as mentioned earlier takes place by two mechanisms: direct dissociation of chains and indirect elimination during double bond formation. A reduction in X results in a lower rate of bridge dissociation and hence in lower concentration of (Y radicals. The production of alkanes by the second mechanism will accordingly decrease. However, this decrease is compensated to some extent by an increase in the direct dissociation, due to increased availability of chains (less chains lost with the tar produced). The net effect is a slight reduction in the amount of alkane gases. The mechanism of ethylene production is more complex and difficult to interpret. 2. A 4% increase in the bridge to unit ratio was implemented by reducing the concentration of units from 5.1 mol/L to 4.9 mol/L. A higher value of this ratio implies more cross-linking among units and a lower probability of forming free units following bridge dissociation. The yield of tar is thereby reduced substantially. The bridge to unit ratio does not affect gas formation by the indirect mechanism, but gases produced by the direct mechanism are increased due to the reduction in chains lost with the departing tar molecules. As a result of these effects the bridge to unit ratio influences the yields of all products. 3. An increase in the methylene to ethylene bridge ratio has a dual effect on product yields. First, since the activation energies of the two types of bridges are 68 and 52 kcal/g-mol, respectively, increasing the ratio increases the average activation energy for bridge dissociation and thereby lowers the tar yield and enhances the gas yield. Secondly, increasing the methylene to ethylene bridge ratio results in a higher total bridge concentration as shown in the calculation of initial conditions, part 1. The higher bridge to unit ratio also results in higher gas and lower tar yields. 4. The activation energy for ethylene bridge dissociation is the most important kinetic parameter. A reduction by 1kcal results in a 16% increase in the yield of tar and an approximately 5% increase in the yield of hydrocarbon gases. The value of 52 kcal used is very close to the estimated value, which in any case is subject to at least f2 uncertainty. 5. Reduction in the rate of phenolic condensation increases tar production by increasing the probability that a unit is free of ether bridges. The increase in the number of chains lost with the tar implies a smaller yield of gases. Thus, phenolic condensation serves to shift the selectivity to favor the gases at the expense of tar. This may be the main reason for the low tar yield of coals containing high amounts of oxygen such as subbituminous and lignites. 6. The inverse relationship between tar and gas formation is also clearly demonstrated by the effect of increasing the activation energies for formation of double

bond bridges. While the amount of tar formed is increased due to a reduction in the concentration of unbreakable bridges, the amount of gases formed is reduced due to the lower rates of indirect elimination reactions. 7. Reduction in the rates of chain dissociation results, as expected, in a higher gas yield. At the same time the yield of tar is also increased by the concomitant reduction in the concentration of (Y radicals. By comparing the sensitivity of gas yield to the rate constants of chain dissociation and double bond formation on bridges, it is clear that the indirect mechanism of gas formation is more effective than the direct mechanism.

Conclusions In the proposed model, coal is represented as a collection of 14 functional groups whose concentrations vary with the type and origin of the coal. Available analytical data were not sufficient to determine all the concentrations and had to be supplemented by the values of four structural parameters that were assigned on the basis of assumptions or treated as adjustable parameters. A set of elementary reactions was chosen guided by organic chemical theory and model compound studies. The associated rate parameters were estimated by the methods of thermochemical kinetics or from published model compound studies. A few rate constants could not be estimated and were treated as adjustable parameters. The formulation of reaction rate expressions was based on the concentrations of reactive configurations calculated from the concentrations of functional groups by combinatorial placement techniques. This procedure resulted in lengthy rate expressions. Computer calculations were carried out for two coals. Several kinetic and structural parameters used in these calculations were adjusted until reasonable agreement with experimental data was obtained. This adjustment was carried out informally and the set of parameters obtained is by no means optimal. Most of the kinetic parameters were not adjusted but were estimated by group additivity methods or taken from the literature. Several kinetic parameters were assigned arbitrarily and differ from the best estimated values. In view of the several adjustable parameters employed, the comparison given between computer calculations and limited experimental data is not intended to show the validity of the structural and mechanistic assumptions used in the model; instead, its ability to describe the experimental trends is shown and is believed to be significant. In particular, the model describes the correlation and interdependence between the yields of tar and gases which was not possible with earlier phenomenological models. The calculations were performed for isothermal conditions to allow comparison with isothermal pyrolysis yields; however, the model can be used equally well with variable temperature. A certain limitation exists with regard to the maximum temperature. Above about 700 "C ring condensation accompaned by the evolution of hydrogen and carbon monoxide become important. Since these reactions were not considered, the model cannot predict the hightemperature evolution of these gases and the ultimate weight loss. However, it would still describe the yields of tar and hydrocarbon gases which are formed relatively rapidly and thus precede significant ring condensation. Mass transfer in the coal phase has been treated empirically by introducing the adjustable parameter X. Mass transfer in the pore space has not been considered; hence the model cannot describe particle size and pressure effects. In its present state the model should not be considered as final and ready to apply but rather as a source of


mechanistic and kinetic information and of various mathematical techniques that may be variously combined in more applied and specialized models. Further work should be directed toward the following improvements: addition of two or three reactions such as ether bond dissociation and radical addition to double bonds; more extensive computations and parameter adjustment by comparison with structural and kinetic data; derivation of simpler rate expressions by using simplified procedures in the calculation of various probability factors.

of the ICs which can be computed, in principle, from the above probabilities. This calculation is, of course, impractical owing to the inability to compute the sum in the denominator. As in statistical mechanics, however, we can use an extremely accurate approximation wherein the mean of the ICs is set equal to their most probable value, attained by maximizing il subject to the constraints (1)-(5). The maximization of Q is most conveniently carried out by using Lagrange multipliers AD,..., XG, corresponding to (1)-(5), and maximizing the quantity J = In Q hD(KDH + KDp + KDG - MD) + ... XG(KDG+ ... + 3 K w -~ MG)

+

Acknowledgment This work was supported by the National Science Foundation under Grant No. AER 74-12161. Appendix A. Calculation of Concentrations of Reactive Configurations We first consider a fixed set of values for the ICs. The subscript of each K represents a certain type of a carbon; e.g., DH is an a carbon carrying a double bond and a hydroaromatic group. The eleven ICs can be considered as eleven types of “balls” of different colors DH, ..., GGG with KDH balls of the first color, etc. Next, we count the number of arrangements of these balls into the Y I 3 “cells”-the a carbons. Each ball has some structure, however, and can be placed in a cell in a number of ways. The DH, DP, DG, RGG, HGG, PGG balls can be placed in three ways each, the RHG, RHP, RPG, HPG balls in six ways each, and the GGG ball in only one way. The number of arrangements of the balls into the cells is therefore given by

+

Using Stirling’s approximation in the form In (K!)= K

In K - K and differentiating J with respect to each K there is obtained KDH = 3pDpH KDP = 3pDpP KDG = 3pDpG KRHG= ~PRPHPG KRHP= ~PRPHPP KHPG= ~ P H P P P G KRGG = 3p@G2 KHPG= ~PHPPPG KHGG = 3pHpG2 KPGG = 3pPpG2

WDH,KGGG)= ~ D H ~ D P~ G G G where

131

KGGG= PG’ where pi

= exi

i = D, R, H, P, G

Introducing (Al) into eq 1-5 we obtain five nonlinear algebraic equations for the five unknowns pD, ...,p G Although these equations can be solved exactly, the numerical calculations were performed with an approximate solution technique described in Jain (1979), p 177. Once the pihave been determined, eq A1 immediately provide the ICs. So far we have not differentiated between the various components of G, namely -H, -CH3, -CHz, -C2Hs, -CHzCHz, and ethylene bridges. To distinguish between the various components of G, we define the following probabilities

P1 = yl/MG PZ = Y z U - Ai - A z ) / M G p3 = y3(1- Ai - A z ) / M G P4 = y4/MG P6

Expanding the binomial coefficients we obtain ~ W D H...,, KGGG)= y13!

KDH!KDp!

x

... KHPG!KGGG!

~ ( K D H...+ +K P D ~ ) ~ ( K ~ G + . . . + K w G )

The probability of the given set of values KDH,...is given by ~(KDH, ...,K m ) / C Q , where the summation is over all sets of Rs that are consistent with the constraints (7)-(11). The desired concentrations are given by the mean values

= Y2AZ/MG

P6 = Y 3 A 2 / M G where Al = Y,/(Y2+ Y3+ Y5)and A2 = Ylz/(Yz+ Y3+ Ys).Thus, P~KRHGis the concentration of a carbons carrying a radical, a hydroaromatic structure, and a hydrogen atom. Similarly, p3pd(pw is the concentration of a carbons carrying a methylene bridge, a (-CzH5) group, and a hydroaromatic (-CH,) group. Thus, probabilities p1 through p6 will be used to characterize the group G. Similarly, we define p H 1 = Y5(1- Al - A 2 ) / M and ~ PH2 = Y J 2 / M H . These two probabilities distinguish between a

Ind. Eng. Chem. Fundam. 1981, 20, 132-137

132

hydroaromatic structure and its /3 radical. Appendix B. Probabilities Relating to Classes X a n d Z of a Carbons To characterize classes x and z we define the following probabilities ml = pr(X possesses a -CH2) =

-

@P$4(KHGG

+ KPGG) + 3P@5(2

-P4

-

P5)KGGG)/[Xl

m2 = pr{X possesses a hydroaromatic structure carrying a P radical) = b $ H 2 [ ( 2 - P 4 - P 5 - ~PG)KHGG + KHPGI)/[Xl c1 = pr(X possesses an a hydrogen) = { 2 p l p 4 [ K H G G + KPGG + 3(2 - P d K G G G I ) / [ X l c2 = pr{Z possesses an a radical) = b 4 [ K R H G + KRPG + (2 - P4)KRGGI)/[Zl

Literature Cited Anthony, D. B.; Howard, J. B. AIChE J. 1976, 79, 625-656. Cheong, P. H. W.D. Thesis, Californla Institute of Technology, 1977. Gan. H.; Nandi, S. P.; Walker, P. L., Jr. Fuel1972, 57, 272-277. Gavalas, G. R.; Oka, M. Fuel 1978, 57, 285-28s. Gavalas, G. R.; Wilks, K. A. AIChE J. 1980, 28, 201-212. Gear, C. W. "Numerical Initial Value Problems in Ordinary Differential Equations", PrenticaHall: Englewood Cliffs, NJ, 1971. Hindmarsh, A. C., Lawrence Livermore Laboratory, Report No. UCID-30001, 1974. Hindmarsh, A. C.; Byme, G. D., Lawrence Livermore Laboratory, Report No. UCID-30112. 1975. Jain, R. Ph.D. Thesis, California Institute of Technology, 1979. Solomon, P. R. "The Evolution of Pollutants During the RapM Devdatilization of Coal", United Technologies Research Center Report No. R76-9525882, 1977. Solomon, P. R. Prepr. Div. Fuel Chem., Am. Chem. Soc. 1979, 24(3), 154- 159.

Received for review May 27, 1980 Accepted January 5, 1981

Derivation of Generalized Darcy Equations for Creeping Flow in Porous Media Ronald G. Larson' Department of Chemical Engineering and MateriSls Science, Unlversi@of Minnesota, Minneapolis, Minnesota 55455

Generalized Darcy equations, that is, equations which express a proportionality between the pressure drop and the volumetric flow rate raised to some power, are herein derived from the creeping flow equations for single or multiphase flow of any fluids for which stress depends homogeneously upon fluid velocity, e.g., Newtonian or power-law fluids. The derivation does not rely upon volumeaveraging techniques. It is found that the simple proportionality between pressure drop and the flow rate raised to a power stems from the flow rate independence of the velocity streamlines.

Introd uction When a single phase flows through a porous medium at low Reynolds number, that is under creeping flow conditions, Darcy's law frequently holds; that is, the superficial velocity, q / A , is proportional to the average gradient of mechanical potential, MIL, and inversely proportional to the fluid viscosity, p (Darcy, 1856). Symbolically _4 -- k A P (1) A CLL where k, the permeability of the porous medium, is the constant of proportionality. The mechanical potential, P, equals the sum of the pressure, p , the gravitational potential, and, if present, the centrifugal potential. Originally proposed by Darcy as an empirical relation, Hubbert in 1956, and others since then (Whitaker, 1969; Slattery, 1969; Lehner, 1979) have derived it by volume averaging the Stokes equations. In the volume averaging approach, a continuum form of Darcy's law is obtained k p = -- v p Ir

(2)

or, for anisotropic porous media (3) *Bell Laboratories,600M o u n t a i n Ave., M u r r a y Hili, NJ 07974. 0196-4313/81/ 1020-0132$01.25/0

where k is the symmetric permeability tensor. v and P in these equations are averages over a volume large compared to the individual pores or cavities, but small compared to the gross dimensions of the porous medium, and h or Zr is position dependent. The gross form of Darcy's law, eq 1, is, for incompressible isothermal flow, an immediate consequence of the continuum form, eq 2. Although it is certainly a more powerful mathematical statement than the gross Darcy equation, the validity of the continuum equation is subject to suspicion. It is clear, for example, that one can define both a (scalar) overall volumetric flow rate through a porous medium and a continuum Navier-Stokes velocity which varies within each pore. But the validity of the volume-averaging approach depends upon the existence of a third velocity defined on an intermediate scale, a scale much larger than pore dimensions, but much smaller than the porous sample. To be well defined, this velocity, 0, must be insensitive to the size and shape of the averaging volume (Whitaker, 1969). The volume-averagedquantities are in practice virtually unmeasurable, and may, in many cases, be undefinable. For example, a packing of spheres in a thin tube whose diameter is only a few sphere diameters is a porous medium for which no volume large compared to pore dimensions yet small compared to tube dimensions exists. More importantly, in naturally occurring porous samples, such as sandstones, large variations in pore structure are usually present on all scales, because of shale streaks and 0 1981 American Chemical Society

Model of Coal Pyrolysis. 2. Quantitative Formulation and Results

Recommend Documents