Ind. Eng. Chem. Res. 2007, 46, 4755-4763
4755
Kinetic Modeling of Acid Catalyzed Hydrocracking of Heavy Molecules: Application to Squalane Eric Vale´ ry,†,‡ Denis Guillaume,*,† Karine Surla,† Pierre Galtier,† Jan Verstraete,† and Daniel Schweich§ Institut Franc¸ ais du Pe´ trole, BP 3, 69390 Vernaison, France, and LGPC, ESCPE-CNRS, 43 bd du 11 noVembre 1918, BP 2077, 69616 Villeurbanne, France
A kinetic model for hydroisomerization of long chain paraffins has been developed by extending the single events modeling methodology to very large reaction networks. The lumped coefficients of the reaction network can now be directly calculated based on a decomposition strategy of the various intervening species into lateral chains and active centers. This approach advantageously avoids the time-consuming generation of the complete exhaustive list of molecules and single events reactions. The properties of lateral chains are computed through recursive formulas. The identified model is confronted to experimental data of squalane hydrocracking. Introduction Using the concept of single events is an efficient way of modeling acid-catalyzed reactions in refining without any reductive assumptions on the composition of the feed.1,2 Good results have been published for the modeling of hydrocracking at lab scale3-6 and at pilot plant scale.7 Extrapolation to industrial units was also successfully achieved in modeling the isomerization of butanes.8,9 Application of the methodology to heavier feeds and especially to molecules with high carbon numbers still remains a bottleneck of the approach, as the computer generated reaction network exponentially increases in size10,11 (both in terms of species and reactions), which leads to untractable networks and long computing times. In order to solve this difficulty, the equations were inspected and reformulated yielding two alternative methods: one based upon the concept of structural classes6 and the other based upon the concept of lateral chains.12,13 Both methods circumvent the network generation. The method based on structural classes, proposed by Martens and Marin,6 consists essentially in counting structures that have the same global symmetry and branching number in order to determine the frequency of each structure by carbon number and degree of branching. This method does therefore not need any network generation. However, for a given network, the structural classes have to be enumerated by the user. Hence, this method seems to need a classification of the structural classes, which has to be manually updated when the maximum degree of branching of the network is increased. The method developed in this work is based on a lateral chain decomposition approach. The computations are entirely automated with the highest carbon number and the highest degree of branching in the network as the only input. This technique was subsequently used to derive a kinetic model for the hydrocracking of squalane, a C30 molecule with 6 methyl substituents. Confrontation to experimental data is also presented in the present paper. Relumping of Large Networks A posteriori lumping of species by carbon number and degree of branching in large single event networks was introduced by * To whom correspondence should be addressed. † Institut Franc ¸ ais du Pe´trole. ‡ Current address: Novasep SAS, BP 50, 54340 Pompey, France. § LGPC, ESCPE-CNRS.
Vynckier and Froment.11 This development is based on the assumption that species inside a lump are in thermodynamic equilibrium. This lumping scheme was not simply devised to handle the lack of component-based analytical data in the modeling of hydrocracking of paraffins, but its underlying assumptions are chemically justified and were reasonably verified by Schweitzer et al.7 and Vale´ry12 for the hydrocracking of n-hexadecane. As usual, the main assumptions are the following: • Hydrogenation/dehydrogenation, protonation/deprotonation are at equilibrium. • Intermediate species (alkenes, carbenium ions) are at pseudo-steady state. • The physical adsorption constants (Langmuir) depend only on the number of carbon atoms of the molecules. • The amount of adsorbed olefins is negligible with respect to the amount of paraffins. We assume further that gas-liquid equilibrium prevails. This enables expressing the rate laws using the partial pressures instead of the concentrations of the dissolved species.6 Let Fo, Fp, and Fq be lumps, i.e., groups of paraffins. In each lump, all hydrocarbons have the same number of carbon atoms and they are in chemical equilibrium. Then, let “reac” be an elementary step, such as a protonated cyclopropane (PCP) or protonated cyclobutane (PCB) isomerization step, Fo f Fp or a β-scission cracking step Fo f Fp + Fq. The lumped rate expression for the disappearance of lump Fo is given by
R
{
reac(m,u) FofFp(+Fq)
[
})
ikl∈
{
∑
reac(m,u) FofFp(+Fq)
kpr
kdep
}
]
(m,Oij)KDHijkreac(m,u)yi
CH+CsatbFoPFo
PH2DENg + DENgacid
(1)
where Csat is the physical adsorption capacity of the zeolite (assumed to be lump independent), Ct is the total concentration of acidic sites, bFo is the physical adsorption constant of lump Fo and PFo is its partial pressure, and “m” and “u” are the type (secondary or tertiary) of the reactant and product carbenium ions.11 The indices “i” refers to paraffin Pi in Fo, “ij” to the jth olefin Oij resulting from dehydrogenation of paraffin Pi, “ik” to the kth carbenium ion Rik+ resulting from protonation of olefin Oij, and “ikl” to the lth activated complex #ikl involved in the reaction that consumes Rik+. Paraffin i, olefin ij, ion ik, and
10.1021/ie061559w CCC: $37.00 © 2007 American Chemical Society Published on Web 06/13/2007
4756
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007
Figure 1. Examples of the two types of activated complexes.
Figure 2. Possible active zones (between square brackets) for β-scission (left) and pcp-isomerization (right). L1 to L5 can be -H, -CH3, or -C2H5 (see text for exceptions).
complex ikl have the same number of carbon atoms within lump Fo. The other parameters are given by
kreac(m,u,ikl) )
σRik+ k˜ (m,u), σikl reac σOij k˜pr(m) kpr (m,Oij) ) kdep σRik+ k˜dep(m,Oij) /
/
k˜dep(m,Oij) ) k˜dep(m,Oref)e(-[∆Gf (Oij)-∆Gf (Oref)]/RT), /
e(-∆Gf(Oij)/RT) )
e(-∆Gf (Oij)/RT) (2) σOij
where k˜X and kX (X ) reac, pr, dep, ...) are respectively the single event kinetic constant and the corresponding elementary step constant, σY (Y ) Oij, Rik+, #ikl, ...) the symmetry number of species Y, ∆Gf(Y) the Gibbs free energy of formation of species Y, ∆G/f its intrinsic Gibbs free energy of formation (rotation symmetries and chirality excluded). Finally,
Figure 3. Several pcp-isomerization steps in a given lump may involve a common activated complex unless the external symmetry of the reacting ion causes the steps to be identical.
Rearranging eq 1 yields the following:14
R
{
reac(m,u) FofFp(+Fq)
(-∆G/f (Pi)/RT)
yi )
e
∑ e(-∆G (P )/RT) P ∈F / f
i
k˜pr(m)
k˜reac(m,u)PFo
i
with
DENg ) 1 +
)
LCFreac(m,u) ofFp(+Fq)
PH2DENg + DENgacid k˜de´ pr(m,Oref)
, KDHij ) e(-[∆Gf(Oij)+∆Gf(H2)-∆Gf(Pi)]/RT),
0
DENgacid
})
σOref
∑ CsatbiKDH K*(Oij T Oref)σ ∀ion ij
Rik+
∑q bF PF q
/
q
(-[∆Gf (Oref(nc))+∆Gf(H2)]/RT) LCFreac(m,u) ) CtCsatbliq × nc e ofFp(+Fq)
1
k˜pr(m) P Pi
k˜dep(m,Oref)
(3) where yi is the mole fraction of paraffin Pi in lump Fo, KDHij the equilibrium constant for dehydrogenation of Pi to Oij, DENg the inhibition term due to physical adsorption, and DENgacid the inhibition term due to “chemisorption” on the acidic sites, i.e., protonation of olefins.
∑e
(-∆G(Pi)/RT)
Pi∈Fo
({
∑
reac(m,u) ikl∈ FofFp(+Fq)
1
}
σ#ikl
)
(4)
In this equation, the energy terms are obtained from Benson’s group contribution method.15 It is important to note that in our approach the last factor in the lumping coefficient (i.e., the sum over the symmetry numbers of the activated complexes) is independent of the produced ions Rqr+. This factor only depends on the reacting ion Rik+ and on the various activated complexes
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007 4757
that result from it. This means that calculating a lumping coefficient requires only knowledge of the ions involved, and thus of the resulting activated complexes, instead of a list of the elementary steps. However, the sum extends over all the steps, and some may involve the very same activated complex #ikl. Historically, the sums were based on the detailed network description. With heavy molecules, the network becomes very large, thus the computation time and the memory requirements become tremendous if still feasible. To solve this problem, Martens and Marin6 proposed a method based on structural classes. The main idea is to generate the few structures that are involved in hydrocarbons that have the same intrinsic Gibbs energy of formation and symmetry number: such a structure is a class. Then, the sums in (4) reduce to sums over the classes weighted by the number of reactions that involve each species in a class. This method is again an enumeration technique, since it requires a listing of the classes and the members of the class; however, it is much faster than the exhaustive generation of all the species. Introducing the intermediate activated complex allows the derivation of an alternate recursive approach to predict the properties of these complexes and especially the sum of the reciprocal of the symmetry numbers.12 We will see that this recursive method also applies for evaluating the sum of energy terms. Hence, the computation of the lumping coefficients is based on an efficient computer algorithm that only needs to know the maximum number of carbon atoms and the maximum branching number of the species in the reaction network. Description of an Activated Complex. In this section, we only consider β-scission and pcp-isomerization. We further assume that lateral branches are methyl or ethyl at most. Every activated complex is composed of three parts (Figure 1): (1) an “activated zone”, ZA, that carries the positive charge and is characteristic of the elementary step considered; (2) a “left” lateral chain, A, and (3) a “right” lateral chain, B; both lateral chains are independent of the reaction type. Based on the structure of the active zone postulated in Figure 1, the generalization of the active zones can be inferred, as shown in Figure 2. The substituents L1 to L5 can be methyl or ethyl or simply hydrogen. Nevertheless, in pcp-isomerization, L1 cannot be ethyl otherwise a propyl branch would be generated, and L3 and L4 cannot be both H in β-scission otherwise a primary carbenium ion would be produced. The global symmetry number σ# of the activated complex # is given by
σ ) σAσZAσBσext(A,ZA,B)
(5)
where σA, σB, and σZA are the partial symmetry numbers of each part and σext(A,ZA,B) accounts for external symmetries that occur for instance when lateral chains are identical and L2 ) L3 in a pcp-isomerization. Recursive Method for Evaluating the Sum of Reciprocals of Symmetry Numbers. The sum of symmetry numbers now becomes
ikl∈
{
∑
reac(m,u) FofFp(+Fq)
1
}
σ#ikl
)
∑ A-ZA-B∈R σ
ncorr(A,ZA,B) AσZAσBσext(A,ZA,B)
reactions of type (m,u) that involve complex A-ZA-B. In β-scissions, a given activated complex is only involved in a single cracking step because the complex is an “intermediate hybrid” between the reacting ion and the resulting ion and olefin (see the left-hand side of Figure 2). As a consequence ncorr ) 1 for any β-scission step. Conversely, pcp-isomerizations may involve the same complex for four reversible isomerization steps as shown in Figure 3. When the lateral chains A and B are different while the total number of carbon atoms is the same in all ions, the corresponding starting paraffins (on the left-hand side of Figure 3) belong to the same lump. There are thus only 4 (s,t) isomerization steps remaining that occur between different lumps and involve the same activated complex, and ncorr ) 4. When A ) B and L1 ) L2, there is a single isomerization step (ncorr ) 1) that involves an activated complex with an external symmetry of 2. We will now focus on β-scission; the reader is referred to Vale´ry12 for details concerning pcp-isomerization. The evaluation of (6) is first based on the following decomposition: • A sum over the various active zones ZA; there are up to 90 different active zones in β-scission. In a given lump, choosing an active zone of nCZA carbon atoms and bZA branches defines the numbers of carbon atoms and of lateral branches still available for all the possible lateral chains A and B. • A subsum over lateral chains A that have a given number of carbon atoms, nCA, and a given number of lateral branches, bA. We will speak of “equivalent lateral chains” for this set of lateral chains defined by nCA and bA. As soon as the class of equivalent lateral chains A and the active zone are given, the class of equivalent lateral chains B is fixed because the total numbers of carbons and of branches are fixed by the lump. Let us remark that in β-scission: (1) ncorr ) 1 (see above); (2) the equivalence class of lateral chains are defined by the lumps of the product ion and olefin; (3) knowing that the active zone has a rigid structure (2-electron 3-center bond), there is no external symmetry and σext ) 1. Hence, eq 6 becomes the following:
(6)
where R is the set of the activated complexes, A-ZA-B, involved in reactions of type (m,u) going from lump Fo to lump Fp (and Fq in a β-scission) and ncorr(A,ZA,B) is the number of
Figure 4 illustrates “homologous” activated complexes involved in the (s,s) β-scission of a tribranched C16 hydrocarbon to a monobranched C9 hydrocarbon (olefin) and a dibranched C7 hydrocarbon (ion). Activated complexes are said to be homologous when they have the same active zone and when the lateral chains A (or B, equivalently) belongs to the same equivalence class. Homologous activated complexes are the individual contributions to each term of the sum over ZA in (7). In the case of Figure 4, there are 32 homologous activated complexes. The sums between square brackets in (7) are calculated by a recursive method. Let nP be the number of carbon atoms of the “main chain” of a lateral chain. The main chain is easily identified, as it is the longest linear chain of the lateral chain. The lateral chain is therefore made of a terminal CH3 group and successive elementary patterns chosen among -CH2-, -CHMe-, -CHEt-, -CMeMe-, -CMeEt-, and -CEtEt-. Each of these patterns contributes to the length of the main
4758
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007
This formula is valid for nP > 3. For low b or nC, the second or third argument of the function US may become negative; in this case, the contribution to US must be set to 0. Ignoring chirality, the symmetry number (1, 3, or 9) of the added pattern is defined by the structure of the allowed lateral branches (-H, -CH3, -C2H5 presently). Chirality (division by 2) is defined as soon as the chain on the left of the inserted pattern is more than 3 carbon atoms long (i.e., different from -H, -CH3, and -C2H5): the central carbon is chiral when its two remaining substituents are different. Further details and formulas for nP e 3 are given in the appendix. We finally define the symmetry property of lateral chain, SPLC(nC,b), by the following: np)nc
SPLC(nc,b) )
∑ US(np,nc,b) ) A∈ ∑ n )n -2b p
{
c
nc carbons b branchings
1 σA
(10)
}
This is exactly the sum within square brackets in eq 7. The starting value of index nP in this sum is obtained when all the branches are ethyl groups, the maximum value when the lateral chain is linear. Equation 7, which was derived for β-scissions, now becomes
∑ A-ZA-B∈Rσ
1
)
AσZAσB
1
∑ ZA σ
SPLC(nCA,bA)SPLC(nCB,bB)
ZA
(11) Figure 4. Example of homologous activated complexes and equivalent lateral chains.
chain. Figure 5 shows the six possibilities for increasing the length of a lateral chain starting from CH3-. If σL is the symmetry number of the lateral chain with nP carbons along the main chain, then the symmetry number after inserting a pattern is σL times the symmetry number of that pattern as shown on the right of Figure 5. Let US(nP,nC,b) be defined by the following:
US(np,nc,b) )
{
A∈
∑ Length of main chain)n
1 p
Number of carbon atoms)nc Number of branches)b
}
σA
(8) σZA )
The recursion formula then becomes
US(np,nc,nb) ) US(np - 1,nc - 1,nb) × 1 - CH2- insertion (σCH2 ) 1) +US(np - 1,nc - 2,nb - 1) ×
2 3
-CHMe- insertion
(σ
CHMe
)
3 2
)
1 +US(np - 1,nc - 3,nb - 2) × 9 -CMeMe- insertion (σCMeMe ) 9) +US(np - 1,nc - 4,nb - 2) ×
2 9
-CEtMe- insertion +US(np - 1,nc - 5,nb - 2) ×
1 9
(σ
CEtMe
In this expression, the symmetry number of the active zone still remains to be determined. It is deduced from simple properties of the activated complex of Figure 2 left. Let nLi ) 0 when Li is a hydrogen substituent and nLi ) 1 otherwise. The internal nLi symmetry number for each Li group is then σint Li ) 3 . Chirality is born by carbons C1, C2, or C3. They are chiral when their substituents are different. Remark that when the lateral chains A and B are more than 2 carbon atoms long, the chirality of the active zone no longer depends on these chains because L1 to L5 are ethyl branchings at most. Finally, the symmetry number of the active zone is given by
)
9 2
)
-CEtEt- insertion (σCEtEt ) 9) (9)
3nL1+nL2+nL3+nL4+nL5 2nch
(12)
where nch is the number of chiral carbon atoms among C1, C2, and C3. Combining this expression with eq 11 therefore allows calculating the last factor in the expression of the lumping coefficient in eq 4 easily by a fast method that does not require the exhaustive generation of the reaction network. Recursive Method for Evaluating the Sum of Energy Terms. The sum of energy terms is calculated by /
∑
Pi∈F0
e
(-∆Gf(Pi)/RT)
)
∑
Pi∈F0
e(-∆Gf (Pi)/RT) σ Pi
(13)
where ∆G/f (Pi) is the intrinsic Gibbs energy of formation of paraffin Pi that can be calculated using Benson’s group contribution method. Here again, a recursive method can be used. For the sake of simplicity, we will ignore the gauche interactions and the external symmetry of the paraffins. However, the method can account for these symmetrical molecules, and the reader is referred to Vale´ry12 for further details. The method is based on the concept of lateral chains as defined above: a given paraffin can be decomposed into two zones, identical to an activated complex without active zone.
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007 4759
Figure 5. Principle of the recursive method for calculating the sum of the reciprocal of symmetry numbers.
A paraffin with nP carbon atoms along the main chain is thus the concatenation of a lateral chain A with nPA carbon atoms and a lateral chain B with nP - nPA carbon atoms along the respective main chains. The symmetry number of the concatenated chains is the product of the symmetry numbers of the chains (external symmetry is ignored!); the intrinsic Gibbs energy is the sum of the individual Gibbs energies (gauche interaction is ignored!); the exponential of the normalized intrinsic Gibbs energy is then the product of the individual exponentials. Hence, eq 13 reduces to
using the “added pattern” method discussed in the previous section. Let /
UT(np,nc,b) )
∑
Pi∈F0 Pi)A-B
σ Pi
)
TPLC(nc,b) )
{
∑
(
∑
/
/
∑
equivalent B chains such that A-B∈F0
)
e(-∆Gf (B)/RT)
}
σB
The factor 1/2 results from the fact that chains A and B can be interchanged, i.e., a paraffin can be considered as either A-B or B-A. The length, nPA, of chain A is arbitrary and remains to be chosen. For reasons of chirality explained in the Appendix, we set nPA ) 3. Then, the properties of chain B can be calculated
)
np)nc np)nc-2b
(14)
(15)
σA
∑
carbons along the main chain
σA
e(-∆Gf (A)/RT)
nc carbon atoms b lateral branchings
2A:chains with nPA e(-∆Gf (A)/RT)
}
σA
/
/
1
∑
A∈ nplength of main chain nc carbon atoms b lateral branchings
A∈
e(-∆Gf (Pi)/RT)
{
e(-∆Gf (A)/RT)
UT(np,nc,b) (16)
where TPLC(nc,b) is the “thermodynamic property of lateral chains”. With nPA ) 3, eq 14 becomes the following: /
∑
Pi∈F0
e(-∆Gf (Pi)/RT) σ Pi
bB)Min(4,b) nCB)Min(nc-1,9)
)
1 2
∑
bB)0
∑
nCB)3
UT(3, nCB, bB)TPLC(nC - nCB,b - bB) (17) Finally, based on Figure 5 and eq 9, UT is recursively given by the following:
4760
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007 /
UT(np,nc,nb) ) UT(np - 1,nc - 1,nb) × 1 × e-∆Gf (CH2)/RT -CH22 +UT(np - 1,nc - 2,nb - 1) × × 3 e
-[∆G/f (CH3)+∆G/f (CH)]/RT
Table 1. Final Parameter Estimates for the Hydrocracking of Squalane kinetic parameter elementary reaction
best value
confidence interval
kpcp(s,s) kpcp(s,t) kcr(s,s) kcr(s,t)a kcr(t,s) kcr(t,t)
0.652 × 102 1.273 × 103 0.407 × 103 0 0.404 × 104 2.056 × 105
(0.654 × 100 (0.464 × 102 (0.359 × 102
-CHMe-
2 +UT(np - 1,nc - 3,nb - 1) × × 3 /
/
/
e-[∆Gf (CH3)+∆Gf (CH2)+∆Gf (CH)]/RT -CHEt-
a
(0.181 × 104 (0.295 × 105
Numerically nonsignificant.
1 +UT(np - 1,nc - 3,nb - 2) × × 9 /
/
e -[2∆Gf (CH3)+∆Gf (C)]/RT -CMeMe2 +UT(np - 1,nc - 4,nb - 2) × × 9 /
/
/
e-[2∆Gf (CH3)+∆Gf (CH2)+∆Gf (C)]/RT -CEtMe1 +UT(np - 1,nc - 5,nb - 2) × × 9 /
/
/
e-[2∆Gf (CH3)+2∆Gf (CH2)+∆Gf (C)]/RT -CEtEt- (18) Further details and the initial UT(3,nC,b) are given in the Appendix. As a conclusion, eqs 9-12 and 15-18 enable calculating the lumping coefficients much faster than generating the complete reaction network and determining the lumping coefficients from the complete list of single events in a generated reaction network. Extension of the Method. The method described above is fairly flexible. It can easily be extended to more general situations involving n-propyl, i-propyl, n-butyl, i-butyl, t-butyl, ... lateral branches. It is sufficient to consider new insertion patterns, other than those considered in Figure 5 and to adapt recursion formulas 9 and 18 together with the initial values (see the Appendix). Equations 10-12, 16, and 17 remain unchanged. The recursive method allows calculating properties other than those presented above. Let
Figure 6. Molar repartition by carbon number: (points) experimental; (lines) simulation.
Figure 7. Mean branching number by carbon atom number.
UX(np,nc,nb) ) UX(np - 1,nc - 1,nb) × X(CH2) -CH2- insertion +UX(np - 1,nc - 2,nb - 1) × X(CHMe) -CHMe- insertion +UX(np - 1,nc - 3,nb - 1) × X(CHEt) -CHEt- insertion +UX(np - 1,nc - 3,nb - 2) × X(CMeMe) -CMeMe- insertion +UX(np - 1,nc - 4,nb - 2) × X(CEtMe) -CEtMe- insertion +UX(np - 1,nc - 5,nb - 2) × X(CEtEt) -CEtEt- insertion np)nc
XPLC(nc,nb) )
∑
np)nc-2nb
UX(np,nc,nb)
(19)
where X(Y) is any property of the inserted pattern Y. In the previous sections, the properties were the reciprocal of the symmetry numbers or the exponential of the normalized Gibbs free energies. If we set X ) 1, then U1 is the number of lateral chains that have nC carbon atoms, b lateral branches, and a main
chain nP carbons long; 1 PLC(nC,b) is thus the number of lateral chains that have nC carbon atoms and b lateral branches. Since a paraffin made of nC + 1 carbon atoms and having b lateral branches is made of a terminal methyl concatenated with a lateral chain of nC carbons and b branches, then 1PLC(nC,b) is the number of paraffins in the lump with nC - 1 carbon atoms and b lateral branches. When X is the square of the reciprocal of a symmetry number or the square of a thermodynamic term involved in (13), one defines series “S2PLC” and “T2PLC” that are involved in the corrections for external symmetries (see the work of Vale´ry12). Up to now, chirality and external symmetry were considered to be independent from each other, although it is wrong for some unfrequent molecules. For instance, 15-methyl-nonacosane is a nonchiral and nonsymmetric molecule. This paraffin is accounted for in (17) by the term UT(3,3,0) × TPLC(27,1). In TPLC(27,1), 15-methyl-nonacosane is accounted for by a n-C14 lateral chain (nP ) 14 > 3) to which the prochiral elementary pattern -CHMe- is added. The resulting 15th carbon is thus considered to be chiral according to the method presented above and summarized in Figure 5. However, this 15th carbon is not chiral in the paraffin since it has two identical substituents, namely, n-C14-. We have no simple solution for this “symmetry coupling” that is presently neglected. This problem is not encountered in SPLC series because the active zone is respon-
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007 4761
Figure 8. Elementary patterns of a lateral chain: nonchiral branched carbon and chiral branched carbon.
sible for a lack of symmetry between the lateral chain A and the complementary ZA-B set.
the problem size and correlation between parameters, the protonation, deprotonation, and reaction rate coefficients were reparameterised as follows:
Model vs Experimentation Experimentation. Experiments were performed in a fixed bed reactor at a pressure of 150 bar and temperature of 375 °C. The reactor was loaded with 80 g of a modified Y-zeolite catalyst. The catalyst was sulfided with a hydrogenated light cycle oil (LCO). The pretreating feed was injected at a liquid hourly space velocity (LHSV) of 1 h-1, with a volumetric H2/HC ratio of 1000 liter at standard temperature and pressure per liter. The temperature is increased to 350 °C in 12 h. The sulfiding is continued until the H2S composition in the effluent remains constant. Before carrying out the kinetic experiments, a line-out period of 150 h at 375 °C with the sulfiding feed was applied in order to reach a stable catalyst performance. As feedstock, pure squalane (2, 6, 10, 15, 19, 23-hexamethyltetracosane) was chosen (purity > 99%). Kinetic data were collected at various LHSVs ranging from 0.5 to 1.5 h-1. The reactor effluent was separated in a flash drum at ambient temperature. The gas fraction was analyzed online, while the liquid effluent was sampled and analyzed on a GC equipped with a polydimethyl siloxane HP PONA column and a flame ionization detector (FID) in order to obtain the carbon number distribution of the products. To assess the isomerization degree as a function of carbon number, the liquid effluent was distilled into narrow cuts for which the degree of branching was determined by 13C NMR. Parameter Identification. The single event kinetic model was incorporated in a plug flow reactor model. Gas and liquid were supposed to be at equilibrium inside the reactor. Phase compositions were calculated by means of a Grayson-Streed based flash module. The differential equations of the reactor were solved with LSODE.16 Parameters were identified by nobs (yi - yˆi)2 minimizing the sum of the squared deviations ∑i)1 with a Levenberg-Marquardt algorithm.17,18 In order to reduce
kreac(m,u) )
k˜pr(m)
k˜reac(m,u)
k˜depr(m,Oref)
(20)
The values of the final parameter estimates are listed in Table 1. The secondary-tertiary cracking rate coefficient was found to be numerically nonsignificant according to the estimation algorithm and was therefore put to zero. This is most certainly due to coupling effects or to indiscernability between (s,t) and (t,s) cracking on the lumps. The performances of the data fit are summarized in the following graphs. Figure 6 shows the comparison between the experimental and the calculated molar fractions of each carbon number fraction during the cracking of a squalane at an LHSV of 0.5 h-1. At a cracking conversion of 86.4%, a good agreement is obtained between the experimental and the calculated cracking conversion as well as for the cracked product distribution (Figure 6). The distribution of the various lumps within each carbon number fraction can be characterized by its degree of branching. Figure 7 illustrates that the model correctly predicts the degree of branching of each carbon number fraction, and hence its composition in terms of the various lumps. Conclusion The generation and relumping of detailed single event reaction network was greatly simplified by introducing a novel approach that allows the direct calculation of the lumping coefficients, thus avoiding the generation of the complete exhaustive list of molecules and single events reactions. The method is based on a decomposition strategy of the various intervening species into lateral chains and active centers and leads to a set of recursive
4762
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007
Table 2. Starting Lateral Chains and Initial Values of US and UT chain involved
US(nP,nC,b)
CH3-
US(1,1,0) ) 1/3
CH3-CH2-
US(2,2,0) ) 1/3
CH3-CHMe-
US(2,3,1) ) 1/9
CH3-CMeMe-
US(2,4,2) ) 1/81
CH3-CH2-CH2-
US(3,3,0) ) 1/3
CH3-CHMe-CH2CH3-CH2-CHMe-
}
US(3,4,1) ) 1/9 + 2/9
CH3-CH2-CHEt-
US(3,5,1) ) 1/9
}
CH3-CH2-CMeMeCH3-CHMe-CHMeCH3-CMeMe-CH2-
}
CH3-CH2-CEtMeCH3-CHMe-CHEtCH3-CHMe-CMeMeCH3-CMeMe-CHMeCH3-CH2-CEtEtCH3-CHMe-CEtMeCH3-CMeMe-CHEt-
US(3,5,2) ) 1/27 + 2/27 + 1/81
}
}
US(3,6,2) ) 1/27 + 2/27 US(3,6,3) ) 1/81 + 2/243 US(3,7,2) ) 1/81 US(3,7,3) ) 2/81 + 2/243
CH3-CMeMesCMeMe-
US(3,7,4) ) 1/729
CH3-CHMe-CEtEt-
US(3,8,3) ) 1/81
CH3-CMeMe-CEtMe-
US(3,8,4) ) 2/729
CH3-CMeMe-CEtEt-
US(3,9,4) ) 1/729
UT(nP,nC,b) 1 UT(1,1,0) ) e 3 / / 1 UT(2,2,0) ) e-[∆Gf (CH3)+∆Gf (CH2)]/RT 3 / / 1 UT(2,3,1) ) e-[2∆Gf (CH3)+∆Gf (CH)]/RT 9 / / 1 UT(2,4,2) ) e-[3∆Gf (CH3)+∆Gf (C)]/RT 81 / / 1 UT(3,3,0) ) e-[∆Gf (CH3)+2∆Gf (CH2)]/RT 3 / / / / / / 1 2 UT(3,4,1) ) e-[2∆Gf (CH3)+∆Gf (CH2)+∆Gf (CH)]/RT + e-[2∆Gf (CH3)+∆Gf (CH2)+∆Gf (CH)]/RT 9 9 / / / 1 UT(3,5,1) ) e-[2∆Gf (CH3)+∆Gf (CH2)+∆Gf (CH)]/RT 9 / / / / / 1 2 UT(3,5,2) ) e-[3∆Gf (CH3)+∆Gf (CH2)+∆Gf (C)]/RT + e-[3∆Gf (CH3)+2∆Gf (CH)]/RT + 27/ 27 1 -[3∆Gf (CH3)+∆G/f (CH2)+∆G/f (C)]/RT e 81 / / / / / / 1 2 UT(3,6,2) ) e-[3∆Gf (CH3)+2∆Gf (CH2)+∆Gf (C)]/RT + e-[3∆Gf (CH3)+∆Gf (CH2)+2∆Gf (CH)]/RT 27 27 / / / 1 2 -[4∆G/f (CH3)+∆G/f (CH)+∆G/f (C)]/RT e UT(3,6,3) ) e-[4∆Gf (CH3)+∆Gf (CH)+∆Gf (C)]/RT + 81 243 1 -[3∆G/f (CH3)+3∆G/f (CH2)+∆G/f (C)]/RT UT(3,7,2) ) e 81 / / / / 2 UT(3,7,3) ) e-[4∆Gf (CH3)+∆Gf (CH2)+∆Gf (CH)+∆Gf (C)]/RT + 81 2 -[4∆G/f (CH3)+∆G/f (CH2)+∆G/f (CH)+∆G/f (C)]/RT e 243 1 -[5∆G/f (CH3)+2G/f (C)]/RT e UT(3,7,4) ) 729 / / / / 1 UT(3,8,3) ) e-[4∆Gf (CH3)+2∆Gf (CH2)+∆Gf (CH)+∆Gf (C)]/RT 81 2 -[5∆G/f (CH3)+3∆G/f (CH2)+2∆G/f (C)]/RT e UT(3,8,4) ) 729 1 -[5∆Gf/(CH3)+2∆Gf/(CH2)+2∆Gf/(C)]/RT UT(3,9,4) ) e 729
formulas that allow the calculation of the properties of the lateral chains. The new technique therefore allows the extension of the single event modeling methodology to very large reaction networks. A kinetic model for hydroisomerisation and hydrocracking of squalane was developed via this technique. After parameter estimation, the identified model was confronted to experimental data. Overall, a good agreement between the predicted and measured product distributions was obtained. The method described above is fairly flexible. It has been shown that it can easily be extended to more general lateral branches. It is possible to introduce new insertion patterns, even new calculated properties, and to adapt recursion formulas consequently. Appendix The smallest lateral chain is CH3-, and the reciprocal of its symmetry number corresponds to US(1,1,0) ) 1/σCH3 ) 1/3. For larger lateral chains, the terms are obtained by adding the elementary patterns given at the top of Figure 8. The problem is to account for chirality. The examples given in Figure 8 show that adding a prochiral pattern next to the third carbon atom of the main chain adds systematically a chiral carbon. Conversely, for the second and third carbon, a prochiral pattern is not always responsible for a chiral carbon (second row of Figure 8). Thus, US(1,1,0) to US(3,9,4) and UT(1,1,0) to UT(3,9,4) require specific calculations. The results for these calculations are given in Table 2.
-[∆G/f (CH3)]/RT
Nomenclature Roman bFo ) physisorption constant for lump Fo Csat ) saturation concentration inside the zeolite CH+ ) concentration of free acid sites Ct ) total concentration of acid sites DENg ) physisorption term (Langmuir type) DENgacid ) adsorption term (Langmuir type) kdep ) deprotonation rate coefficient kpr ) protonation rate coefficient kpcp ) PCP branching rate coefficient kcr ) cracking rate coefficient nb ) number of branches nc ) number of carbon atoms Oref ) reference olefin Pi ) paraffin i PFo ) partial pressure of lump Fo PH2 ) partial pressure of hydrogen R{reac(m,u)FofFp} ) rate of reaction reac, of type (m,u), between lump Fo and lump Fp R ) ideal gas constant T ) temperature yi ) experimental mole fraction of observable i yˆi ) calculated mole fraction of observable i Greek ∆Gf ) gibbs energy σ ) symmetry number
Ind. Eng. Chem. Res., Vol. 46, No. 14, 2007 4763
Literature Cited (1) Froment, G. F. Kinetic modeling of acid-catalyzed oil refining processes. Catal. Today 1999, 52, 153-163. (2) Froment, G. F. Single event kinetic modeling of complex catalytic processes. Catal. ReV. 2005, 47, 83-124. (3) Baltanas, M. A.; Van, Raemsdonck, K. K.; Froment, G. F.; Mohedas, S. R. Fundamental kinetic modeling of hydroisomerization and hydrocracking on noble-metal-loaded faujasites. 1. Rate parameters for hydroisomerization. Ind. Eng. Chem. Res. 1989, 28, 899-910. (4) Svoboda, G. D.; Vynckier, E.; Debrabandere, B.; Froment, G. F. Single event rate parameters for paraffin hydrocracking on a Pt/US-Y zeolite. Ind. Eng. Chem. Res. 1995, 34, 3793-3800. (5) Martens, G. G.; Froment, G. F. Kinetic modeling of paraffins hydrocracking based upon elementary steps and the single event concept. In Reaction kinetics and the deVelopment of catalytic processes; Froment, G. F., Waugh, K. C., Eds.; Elsevier Science BV: New York, 1999. (6) Martens, G. G.; Marin, G. B. Kinetics for hydrocracking based on structural classes: model development and application. AIChE J. 2001, 47 (7), 1607-1622. (7) Schweitzer, J. M.; Galtier, P.; Schweich, D. A single events kinetic model for hydrocracking of paraffins in a three phase reactor. Chem. Eng. Sci. 1999, 54, 2441-2452. (8) Guillaume, D.; Surla, K.; Galtier, P. From Single Events theory to molecular kinetics-application to industrial process modeling. Chem. Eng. Sci. 2003, 58 (21), 4861-4869. (9) Surla, K.; Vleeming, H.; Guillaume, D.; Galtier, P. A single events kinetic model: n-butane isomerization. Chem. Eng. Sci. 2004, 59 (2223), 4773-4779. (10) Baltanas, M. A.; Froment, G. F. Computer generation of reaction networks and calculation of product distributions in the hydroisomerization
and hydrocracking of paraffins on Pt-containing bifunctional catlysts. Comput. Chem. Eng. 1985, 9, 71-81. (11) Vynckier, E.; Froment G. F. Modeling of the kinetics of complex processes upon elementary steps. In Kinetic and thermodynamic lumping of multicomponent mixtures; Astarita, G., Sandler, S.I., Eds.; Elsevier Science Publishers BV: Amsterdam, 1991; pp 131-161 (12) Vale´ry, E. Application de la the´orie des e´ve´nements constitufs a` l’hydrocraquage de paraffines lourdes. Ph.D. Thesis, Universite´ Claude Bernard Lyon I, Lyon, France, 2002. (13) Guillaume, D.; Vale´ry, E.; Surla, K.; Galtier, P.; Verstraete, J. Single Events Modeling-extension to large networks. Communication at ECCE4, Grenade, Espagne, September 21-25, 2003. (14) Vale´ry, E. Mode´ lisation cine´ tique des re´ actions catalytiques d’hydrocraquage par la the´ orie des e´ Ve´ nements constitutifs. Internal IFP Report, 2000. (15) Benson, S. W. Thermochemical kinetics, second ed.; Wiley & Sons: New York, 1976. (16) Hindmarsch, A. C. Odepack, a systematized collection of ODE solvers. In Scientific Computing, Stepleman, R.S., et al., Eds.; NorthHolland: Amsterdam, 1983; pp 55-64. (17) Levenberg, K. A Method for the Solution of Certain Problems in Least Squares. Q. Appl. Math. 1944, 2, 164-168. (18) Marquardt, D. An Algorithm for Least-Squares Estimation of Nonlinear Parameters. SIAM. J. Appl. Math. 1963, 11, 431-441.
ReceiVed for reView December 5, 2006 Accepted February 8, 2007 IE061559W
