Article pubs.acs.org/IECR
Effective Kinetic Modeling for Homogeneous Reaction Networks with Branches from Loops Frederick Chern, Tai-Shang Chen, and Jia-Ming Chern* Department of Chemical Engineering, Tatung University, 40 Chungshan North Road, Third Sec. Taipei, 104, Taiwan ABSTRACT: The algorithm of Helfferich to derive the general yield ratio and rate equations for multistep homogeneous reactions with simple nodes and loops is extended to cope with reaction networks with branches from nodes and nonsimple networks. Systematically using the Bodenstein approximation of quasi-stationary behavior of reaction intermediates and network reduction technique, the yield ratio and rate equations are readily compiled for networks of arbitrary configuration and complexity. The applications of the algorithm are illustrated by some homogeneous reactions.
■
INTRODUCTION Reaction and separation are two central themes in chemical engineering education and chemical industries. Although chemical reactor units, compared with separation units, usually represent a minor portion in process flow diagrams, they are of vital importance. Adequate chemical reactor design and control can save the efforts of separating products from reactants. The fundamentals and design aspects of chemical reactors have been well documented.1−5 For design, scale-up, optimization, online control, and trouble shooting of chemical reactors, rate equations for reactants and products as functions of the reactant and/or product concentrations and temperature are indispensable.6,7 Since many chemical reactions of industrial importance involve complex reaction pathways and networks, the determination of the reaction rate laws becomes difficult. Therefore, power-law type rate equations are usually adopted to correlate the experimental data in laboratory-scale and pilotplant-scale reactors. These power-law type rate equations are empirical in nature and might not be used confidently in scaleup and optimization of chemical reactors.8 Safe scale-up of reactors by large factors calls for reaction-kinetic modeling that correctly reflects the events at the molecular level, that is, mathematics based on reaction networks or pathways.7 A methodology of systematically using the Bodenstein approximation of quasi-stationary behavior of reaction intermediates and a network reduction technique were applied to derive the general rate and yield ratio equations of multistep homogeneous reactions with simple nodes and loops.9 The same methodology was extended to derive general yield ratio and rate equations for multistep homogeneous reactions with snowflake-type networks.10 The algorithm of the general rate equation was applied to study some homogeneous reactions of environmental interest.11−13 When applying the algorithm to study more homogeneous reactions, we found that the previous algorithm was not applicable to many reaction networks consisting of branches from loops or nonsimple networks with reaction steps higher than first order in intermediates. Therefore, a new algorithm is developed to allow fast derivation of the rate equations for such reaction networks and illustrated by some multistep homogeneous reactions. © 2014 American Chemical Society
2. METHODOLOGY As shown in the previous publications,7,9,10 the Bodenstein approximation of quasi-stationary behavior of intermediates permits any “simple” network (that is, with all intermediates at trace level and no steps of higher order in intermediates) to be reduced to one with only pseudosingle steps between nodes and end points. Specifically, the net rate contribution of a multistep, reversible simple network segment between nodes Xj and Xk X j ↔ X j + 1 ↔ ... ↔ X j + n ↔ ... ↔ X k
(1)
(coreactants and coproducts not shown) is rj → k = Λjk[X j] − Λkj[X k]
(2)
where the “segment coefficients” Λ are given by k−1
∏i = j λi , i + 1
Λjk =
Djk
k−1
and Λkj =
∏i = j λi + 1, i Djk
(3)
with k
∑
Djk =
(
i−1
k−1
∏
λm , m − 1 ∏ λm , m + 1)
i=j+1 m=j+1
m=i
(4)
(products Π = 1 if lower index exceeds upper). More easily, Djk is generated as the sum of the products of the rows of the square matrix of order n with elements 1 along the diagonal, with rate coefficients of mth step in the mth column, and with forward rate coefficients above and reverse rate coefficients below the diagonal: Received: Revised: Accepted: Published: 12983
May 14, 2014 July 24, 2014 July 28, 2014 July 28, 2014 dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research ⎡ 1 ⎢ ⎢ λj + 1, j ⎢ ⎢ λj + 1, j ⎢ ⎢ ... ⎢ λj + 1, j ⎢ ⎢⎣ λj + 1, j
λj + 1, j + 2 λj + 2, j + 3 ... λk − 2, k − 1 λk − 1, k ⎤ ⎥ 1 λj + 2, j + 3 ... λk − 2, k − 1 λk − 1, k ⎥ ⎥ λj + 2, j + 1 1 ... λk − 2, k − 1 λk − 1, k ⎥ ⎥ ... ... ... ... ... ⎥ λj + 2, j + 1 λj + 3, j + 2 ... 1 λk − 1, k ⎥ ⎥ λj + 2, j + 1 λj + 3, j + 2 ... λk − 1, k − 2 1 ⎥⎦
Article
steps Xi → J, with J being the immediate neighbors of Xi in the network (intermediates or end members).
3. THREE-NODE REACTION NETWORK Consider the following three-node reaction network with reduced pseudosingle steps between the nodes and end members shown in Figure 1a. (5)
For example, for j = 0 and k = 4 (three intermediates): D04 = λ12λ 23λ34 + λ10λ 23λ34 + λ10λ 21λ34 + λ10λ 21λ32
(6)
The λ coefficients are pseudo-first-order rate coefficients of quasi-single molecular steps and are the products of the actual rate coefficients and the concentrations of any coreactants of the respective steps. For example, for the step X0 + A → X1, λ01 = k01CA; Xi + Xj → Xk, λij = kijCXj = kij′. The Bodenstein approximation of quasi-stationary behavior of intermediates lies on the fact that the intermediate concentrations are much smaller compared with those of end-member reactants or products. As the half-life of the intermediates in the reaction mixture is exceedingly short. They must therefore be assumed to disappear at nearly the same rate as the one with which they are formed, the respective concentrations therefore being stationary in time.14 Since the intermediate concentration CXj is assumed to be constant in time, the product of this constant with another constant, the rate coefficient kij becomes a constant in time. Therefore, the Xi + Xj → Xk step can be viewed as a reduced step XI → Xk with a pseudo-first-order rate coefficient of kij′. Therefore, the algorithm developed in this study is not restricted to networks containing steps involving two or more molecules of intermediates as reactants. With the introduction of the pseudo-first-order rate coefficients, the intermediate concentrations can be obtained from the matrix operation as shown in the Appendix. The rate equations for the end members can thus be obtained by the developed algorithm. However, the algorithm developed in this study is applicable only to homogeneous reactions. For heterogeneous catalysis or homogeneous catalysis, total site or enzyme balance must be taken into account to obtain the intermediate concentrations. The resultant formulas for the intermediate concentrations are quite different; therefore, the developed algorithm cannot be applied in heterogeneous catalysis or homogeneous catalysis. Following a similar methodology derived by Chern and Helfferich,10 the rate law of formation of any product Pj (or of consumption of any reactant Pj) in a multinode network with cycle configuration can be compiled with the following rules: (i) The rate law rpj is composed of additive contributions from all pi ↔ pj (i ≠ j), each expressing the net rate at which Pj is formed from (or, if negative, reverts to) the respective Pi. (ii) Each contribution pi ↔ pj consists of a factor (ΠijCpi − ΠjiCpj, a factor |Mij|, and a factor 1/|det M|. Here, Πij is the product of all forward Λ coefficients on the path from Pi to Pj; Πji is the product of all reverse Λ coefficients on that path; M is a matrix of order n, with elements −Si (i = 1, 2, ..., n) along its diagonal, elements Λkm in row m and column k for all Xk and Xm, directly connected with one another, and elements zero in all other positions; Mij = Mji is a matrix similarly compiled, but with mth rows and columns omitted, where the m’s are the indices of all intermediates on the path from Pi to Pj, and the node coefficient Si is the sum of the Λ coefficients of the reduced
Figure 1. Three-node reaction network with pseudo-single steps between nodes and end members.
The rate law rp2 is composed of five contributions of net reaction rates between P2 and the other three end members P1, P3, and P4. As an example of the compilation of these contributions, for pathway 1 in Figure 1b: rP1↔ X1↔ X 2 ↔ P2 =
(Λ P11Λ12Λ 2P2C P1 − Λ1P1Λ 21Λ P22C P2)|M12| |det M| (7)
Because the pathway does not pass intermediate X3, |M12| = S3 = Λ31 + Λ32 + Λ3P4. For pathway 2 in Figure 1c: rP1↔ X1↔ X3↔ X 2 ↔ P2 =
Λ P11Λ13Λ32Λ 2P2C P1 − Λ1P1Λ31Λ 23Λ P22C P2 |det M| (8)
For pathway 3 in Figure 1d: rP4 ↔ X3↔ X1↔ X 2 ↔ P2 =
Λ P43Λ31Λ12Λ 2P2C P4 − Λ3P4Λ13Λ 21Λ P22C P2 (9)
|det M|
For pathway 4 in Figure 1e: rP4 ↔ X3↔ X 2 ↔ P2 =
(Λ P43Λ32Λ 2P2C P4 − Λ3P4Λ 23Λ P22C P2)|M32| |det M| (10)
Because the pathway does not pass intermediate X1, |M32| = S1 = Λ12 + Λ13 + Λ1P1. For pathway 5 in Figure 1f: rP3 ↔ X 2 ↔ P2 =
(Λ P32Λ 2P2C P3 − Λ 2P3Λ P22C P2)|M 22| |det M|
(11)
Because the pathway does not pass intermediate X1 and X3, |M22| = S1S3 − Λ31Λ13. 12984
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
The reaction rate for P2, rp2 is the sum of the above five pathway rates. The common denominator in all of the pathway rates is the absolute value of the following determinant:
Pathway 5 as shown by Figure 2e: rP4 ↔ X3↔ X1↔ P1 =
−S1 Λ 21 Λ31 det M =
|det M| (18)
Because the pathway does not pass intermediate X2, |M13| = S2 = Λ21 + Λ23 + Λ2P2 + Λ2P3. Pathway 6 as shown by Figure 2f:
Λ12 −S2 Λ32 Λ13 Λ 23 −S3
(12)
where
rP4 ↔ X3↔ X 2 ↔ X1↔ P1 =
S1 = Λ12 + Λ13 + Λ1P1
Λ P43Λ32Λ 21Λ1P1C P4 − Λ3P4Λ 23Λ12Λ P11C P1 |det M| (19)
S2 = Λ 21 + Λ 23 + Λ 2P2 + Λ 2P3 S3 = Λ31 + Λ32 + Λ3P4
In some cases, the rate equations of the coreactants and/or coproducts in the reaction steps between the nodes are of interest. In this case, the reaction rates between the nodes can be calculated. For example, rX1→X2 = Λ12CX1 − Λ21CX2, rX1→X3 = Λ13CX1 − Λ31CX3, and rX2→X3 = Λ23CX2 − Λ32CX3. In order to calculate the reaction rates between the nodes, the intermediate concentrations can be expressed in terms of the end member concentrations. (see Appendix). The developed algorithm can be extended to chemical reaction networks with an arbitrary number of nodes and branches and will be presented in future publications. Using this algorithm, the yield ratio equations of the end members can be easily obtained.7,9,10
(13)
The above algorithm was derived by solving the following matrix equation for the intermediate node concentrations obtained from the Bodenstein approximation and using the end-member rate definition rp2 = Λ2P2CX2 − ΛP22CP2 (see Appendix). Using the same algorithm, the rate equations for other end members can be readily obtained; the rate law for P3 is symmetrical to that of P2, while the rate laws for P1 and P4 are also symmetrical. For example the rate law for P1 is composed of six contributions from six pathways. Pathway 1 as shown by Figure 2a: rP2 ↔ X 2 ↔ X1↔ P1 =
(Λ P43Λ31Λ1P1C P4 − Λ3P4Λ13Λ P11C P1)|M13|
4. APPLICATIONS It is very easy to apply the developed methodology to highly complex reaction networks. Several steps are involved in
(Λ P22Λ 21Λ1P1C P2 − Λ 2P2Λ12Λ P11C P1)|M12| |det M| (14)
Table 1. Hypothetical Mechanisms of Two Parallel Reactions: A + B → P and A + C → Q mechanism I k 01, k10
A + cat ←⎯⎯⎯→ X1 k12
A + cat ←⎯⎯⎯→ X1
X1 ⎯→ ⎯ X2
k13
k13
X1 → X3 k 24
k46
X4 ⎯→ ⎯ P + cat k35
X3 + C ⎯→ ⎯ X5
X 5 ⎯→ ⎯ Q + cat
k 23
X 2 ⎯→ ⎯ X3 k 24
X 2 + B ⎯→ ⎯ X4 k46
X4 ⎯→ ⎯ P + cat k35
X3 + C ⎯→ ⎯ X5 k57
rP2 ↔ X 2 ↔ X3↔ X1↔ P1 =
k 01, k10
A + cat ←⎯⎯⎯→ X1 k12
X1 ⎯→ ⎯ X2 X1 → X3
X 5 ⎯→ ⎯ Q + cat
Pathway 2 as shown by Figure 2b:
mechanism III
k13
X1 → X3
X 2 + B ⎯→ ⎯ X4
k57
k 01, k10
k12
X1 ⎯→ ⎯ X2
Figure 2. Three-node reaction network with six pathways for rp1.
mechanism II
k32
X3 ⎯→ ⎯ X2 k 24
X 2 + B ⎯→ ⎯ X4 k46
X4 ⎯→ ⎯ P + cat k35
X3 + C ⎯→ ⎯ X5 k57
X 5 ⎯→ ⎯ Q + cat
Λ P22Λ 23Λ31Λ1P1C P2 − Λ 2P2Λ32Λ13Λ P11C P1 |det M| (15)
Pathway 3 as shown by Figure 2c: rP3 ↔ X 2 ↔ X1↔ P1 =
(Λ P32Λ 21Λ1P1C P3 − Λ 2P3Λ12Λ P11C P1)|M12|
Figure 3. Reaction networks of two parallel reactions with different mechanisms.
|det M| (16)
Pathway 4 as shown by Figure 2d: rP3 ↔ X 2 ↔ X3↔ X1↔ P1 =
applying the general rate equations. First of all, a proposed reaction mechanism should be appropriately arranged as a cyclic reaction network. Second, the general rate equation of the corresponding network in terms of the node coefficients, S,
Λ P32Λ 23Λ31Λ1P1C P3 − Λ 2P3Λ32Λ13Λ P11C P1 |det M| (17) 12985
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
Table 2. Rate Equations of A + B → P and A + C → Q with Different Mechanisms mechanism I
rate equation
ΛA1 = k 01Ccat ,
Λ1A = k10 ,
lumped parameters
Λ1P = k12 , Λ1Q = k13
ka =
k 01k12Ccat k10 + k12 + k13
kb =
k 01k13Ccat k10 + k12 + k13
kc =
ka kb
S1 = Λ1A + Λ1P + Λ1Q − rB = rP = − rC = rQ =
k 01Ccatk12CA ΛA1Λ1PCA = = kaCA S1 k10 + k12 + k13 ΛA1Λ1Q CA S1
=
k 01Ccatk13CA = k bCA k10 + k12 + k13
− rA = rP + rQ = (ka + k b)CA RP/Q ≡ II
rP = kc rQ
ΛA1 = k 01Ccat , Λ 23 = k 23 ,
Λ1A = k10 ,
Λ 2P = k 24C B ,
S1 = Λ1A + Λ12 + Λ13, − S1 |det M| =
0
Λ12 ‐S2
Λ12 = k12 ,
ka = k 01Ccat
Λ13 = k13
Λ3Q = k 35CC
S2 = Λ 2P + Λ 23,
=
kc =
k10 k12
kd =
k13 k12
= S1S2S3
Λ13 Λ 23 ‐S3 − rB = rP =
k 24 k 23
S3 = Λ3Q
0 0
kb =
ΛA1Λ12Λ 2PCAS3 Λ Λ Λ C = A1 12 2P A S1S2S3 S1S2
k 01Ccatk12k 24C BCA kak bCAC B = (k10 + k12 + k13)(k 24C B + k 23) (1 + kc + kd)(1 + k bC B)
− rC = rQ =
ΛA1Λ13Λ3Q CAS2 + ΛA1Λ12Λ 23Λ3Q CA S1S2S3
=
k 01Ccatk13k 35CCCA(k 24C B + k 23) + k 01Ccatk12k 23k 35CCCA (k10 + k12 + k13)(k 24C B + k 23)k 35CC
=
k 01Ccatk13CA(k 24C B + k 23) + k 01Ccatk12k 23CA (k10 + k12 + k13)(k 24C B + k 23)
=
kakdCA(1 + k bC B) + kaCA (1 + kc + kd)(1 + k bC B)
− rA = rP + rQ = RP/Q ≡
kak bCAC B + kakdCA(1 + k bC B) + kaCA (1 + kc + kd)(1 + k bC B)
kak bCAC B k bC B rP = = rQ kakdCA(1 + k bC B) + kaCA kd(1 + k bC B) + 1
12986
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
Table 2. continued mechanism III
rate equation
ΛA1 = k 01Ccat , Λ32 = k 32 ,
Λ1A = k10 ,
Λ 2P = k 24C B ,
S1 = Λ1A + Λ12 + Λ13, − S1 |det M| =
0
Λ12 − S2
lumped parameters
Λ12 = k12 ,
Λ13 = k13
ka = k 01Ccat
Λ3Q = k 35CC
S2 = Λ 2P,
k 01Ccatk12k 24C BCA(k 35CC + k 32) + k 01Ccatk13k 32k 24C BCA (k10 + k12 + k13)(k 24C B)(k 35CC + k 32)
=
k 01Ccatk12CA(k 35CC + k 32) + k 01Ccatk13k 32CA (k10 + k12 + k13)(k 35CC + k 32)
=
kaCA(1 + k bCC) + kakdCA (1 + kc + kd)(1 + k bCC) ΛA1Λ13Λ3Q CAS2 S1S2S3
=
kd =
k13 k12
ΛA1Λ13Λ3Q CA
=
k 01Ccatk13k 35CCCA (k10 + k12 + k13)(k 35CC + k 32)
=
kak bkdCACC (1 + kc + kd)(1 + k bCC)
− rA = rP + rQ = RP/Q ≡
k10 k12
ΛA1Λ12Λ 2PCAS3 + ΛA1Λ13Λ32Λ 2PCA S1S2S3
=
− rC = rQ =
kc = = S1S2S3
Λ13 Λ 23 − S3 − rB = rP =
k 35 k 32
S3 = Λ3Q + Λ32
0 0
kb =
S1S3
kaCA(1 + k bCC) + kakdCA + kak bkdCACC (1 + kc + kd)(1 + k bCC)
k C (1 + k bCC) + kakdCA (1 + k bCC) + kd rP = a A = rQ kak bkdCACC k bkdCC
Table 3. Reaction Mechanism of Ethane Pyrolysis reaction k1 •
•
C2H6 → CH3 + CH3 k2
•
CH3 + C2H6 → •C2H5 + CH4
•
k3
C2H5 → •H + C2H4
•
Figure 4. Yield ratios versus initial B concentration plots for different mechanisms.
k4
H + C2H6 → •C2H5 + H 2
T3−2 T3−3 T3−4
•
k5
T3−5
•
k6
T3−6
C2H5 + •C2H5 → C4 H10
C2H5 + •C2H5 → C2H6 + C2H4
and the Λ coefficients is applied. Finally, the S and the Λ coefficients are replaced by the products of ordinary rate coefficients and concentrations of respective coreactants if any. The following examples will demonstrate this application procedure. 4.1. Parallel Reactions with Rival Mechanisms. Consider two parallel homogeneous reactions, A + B → P and A + C → Q, with hypothetical reaction mechanisms shown in Table 1. In each hypothetical mechanism, the reaction is initiated by the reaction of the reactant A with a coreactant cat to form an intermediate X1, followed by sequential or branched reactions, and cat is formed as a coproduct in the last step. It is important to note that species X1 to X5 are not adsorbed species in heterogeneous catalysis. Instead, they are intermediate species with a very short half-life, compared with the reactants A, B, and C and the products P and Q, in the homogeneous reactions.
step T3−1
•
k7
T3−7
k8
T3−8
CH3 + H 2 → •H + CH4
•
•
H + C2H4 → •C2H5 k9
C2H5 + C2H4 → •CH3 + C3H6
T3−9
In the traditional method, eight ordinary differential equations are needed to analyze batch reaction data of the three reactants A, B, and C and the five intermediate products X1 to X5. Eight rate coefficients are needed for mechanism I, and nine rate coefficients are needed for mechanisms II and III. To apply the methodology developed in this study, the reaction mechanisms are first rearranged into the reaction networks shown in Figure 3. As is shown by Figure 3a, the reaction network for mechanism I belongs to the snowflake type; the methodology 12987
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
Figure 6. Reaction network of ethane thermal cracking at high ethane conversion.
needed to model the reaction kinetics for the two parallel reactions with mechanism I; three ordinary differential equations and four lumped parameters are needed for the reactions with mechanisms II and III. This approach without assuming that any reaction step is a rate-determining step will adequately describe the reaction kinetics of the two parallel reactions and greatly reduce the efforts for the kinetic data analysis and parameter estimation. Before detailed kinetic analysis and parameter estimation is carried out for given overall reactions, the rival mechanisms should be discriminated by comparing with the experimental data. Usually, the initial yield ratios obtained with different initial reactant concentrations are compared with those predicted by all rival mechanisms. For example, the predicted yield ratios of P over Q versus initial B concentration can be plotted schematically as in Figure 4 to discriminate the rival mechanisms. As is shown by Table 2, the yield ratio RP/Q for mechanism I is a constant, independent of the initial concentrations of B and C. Therefore, R P/Q vs C B plots at varying initial C concentrations should be overlapping straight lines as shown in Figure 4a. Also shown by Table 2, the yield ratio RP/Q for mechanism II, independent of the initial C concentration, linearly increases with the initial concentration of B for lower B concentrations but approaches a plateau for very high B concentrations. Therefore, RP/Q vs CB plots at varying initial C concentrations should be overlapping curves as shown in Figure 4b. Also shown by Table 2, the yield ratio RP/Q for mechanism III, independent of the initial B concentration, decreases with increasing initial concentration of C as shown in Figure 4c. Once a possible reaction mechanism is identified using the initial rate data, some of the lumped parameters can be obtained from suitable linearized plots of the initial rate data. For mechanism I, the slope of the rP versus CA plot gives the lumped parameter ka while that of rQ versus CA plot gives the lumped parameter kb. For mechanism II, the lumped parameters kb and kd can be determined from the slope and intercept of the CB/(RP/Q) versus CB plot: CB/RP/Q = ((1 + kd)/kb) + kdCB. For mechanism III, the lumped parameters kb and kd can be determined from the slope and intercept of the RP/Q versus 1/ CC plot: RP/Q = (1/kd) + ((1 + kd)/kbkd)(1/CC). 4.2. Ethane Thermal Cracking. Thermal cracking of ethane at temperatures in the range of 800 to 1200 K and ambient or lower pressures yields mainly ethene and hydrogen, but small amounts of methane, butane, and propene are also formed. The reaction is generally held to be essentially first order in ethane.7 According to Rice and Herzfeld,15 the reaction mechanism consists of the reaction steps T3−1 to T3−6 at low ethane conversion as listed in Table 3.
Figure 5. Reaction network of ethane thermal cracking at low ethane conversion.
Table 4. Rate Equations of Ethane Pyrolysis at Low Conversion ΛC2H6 • C2H5 =
k1k 2CC2H6 k 2CC2H6
= k1, Λ•C2H5• H = k 3 , Λ•H • C2H5 = k4CC2H6
Λ•C2H5C4H10 = k5C•C2H5 = k5′, Λ•C2H5C2H4 = k6C•C2H5 = k6′ S•C2H5 = k 3 + k5′ + k6′, S•H = k4CC2H6
|det M| =
C•C2H5 = C•H =
− (k 3 + k5′ + k6′)
k4CC2H6
k3
− k4CC2H6
ΛC2H6 • C2H5CC2H6S•H |det M|
=
ΛC2H6 • C2H5Λ•C2H5• HCC2H6 |det M|
= (k5′ + k6′)k4CC2H6
k1CC2H6k4CC2H6 (k5′ + k6′)k4CC2H6 =
=
k1k 3CC2H6 (k5′ + k6′)k4CC2H6
k1CC2H6 k5′ + k6′ =
k1k 3 (k5′ + k6′)k4
rC2H6 = − 2ΛC2H6 • C2H5CC2H6 − k4CC2H6C•H + k6′C•C2H5 =− 2k1CC2H6 − k4CC2H6 =− 2k1CC2H6 −
k1CC2H6 k1k 3 + k6′ k5′ + k6′ (k5′ + k6′)k4
k1k 3CC2H6
+
k1k6′CC2H6
k5′ + k6′ (k5′ + k6′) ⎛ 2k1k5′ + k1k6′ + k1k 3 ⎞ = −⎜ ⎟CC2H6 k5′ + k6′ ⎝ ⎠ ∝ C C2H 6
rCH4 = ΛC2H6 • C2H5CC2H6 = k1CC2H6 ∝ CC2H6
rC2H4 = k 3C•C2H5 + k6′C•C2H5 = rC4H10 = k5′C•C2H5 = rH2 = k4CC2H6C•H =
k1k5′CC2H6 k5′ + k6′ k1k 3CC2H6 k5′ + k6′
k1(k 3 + k6′)CC2H6 k5′ + k6′
∝ C C2H 6
∝ C C2H 6 ∝ C C2H 6
developed by Chern and Helfferich10 can be directly applied for kinetic analysis. Because of the existence of the steps between intermediates X2 and X3, the reaction networks for mechanisms II and III are not snowflake-type networks and the previous methodology cannot be applied. Using the previous and new methodologies, the rate equations can be readily derived and shown in Table 2. It is important to note that only three ordinary differential equations and three lumped parameters are 12988
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
Table 5. Intermediate Radical Concentrations of Ethane Pyrolysis at High Conversion Λ•C2H5C4H10 = k5C•C2H5 = k5′,
Λ•C2H5C2H4 = k6C•C2H5 = k6′
S•CH3 = k 2CC2H6 + k 7C H2 S•C2H5 = k 3 + k5′ + k6′ + k 9CC2H4 S•H = k4CC2H6 + k 8CC2H4
− (k 2CC2H6 + k 7C H2)
k 9CC2H4
0
k 2CC2H6
− (k 3 + k5′ + k6′ + k 9CC2H4)
k4CC2H6 + k 8
k3
− (k4CC2H6 + k 8
|det M| = k 7C H2
C C2H 4
C C2 H 4) =(k 2CC2H6 + k 7CH2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′) k1CC2H6
− (k 3 + k5′ + k6′ + k 9CC2H4)
k4CC2H6 + k 8CC2H4
k3
− (k4CC2H6 + k 8CC2H4)
C•CH3 = = =
|det M|
k1CC2H6(k4CC2H6 + k 8CC2H4)(k5′ + k6′ + k 9CC2H4) (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′) k1CC2H6(k5′ + k6′ + k 9CC2H4) (k 2CC2H6 + k 7C H2)(k5′ + k6′)
C•C2H5 =
k1k 2CC2H6 2S•H |det M|
+
k1k 7CC2H6C H2k4CC2H6 |det M|
+
k1k 7CC2H6C H2k 8CC2H4 |det M|
2
k1k 2CC2H6 (k4CC2H6 + k 8CC2H4) + k1k 7CC2H6C H2(k4CC2H6 + k 8CC2H4)
=
= =
(k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′) k1k 2CC2H6 2 + k1k 7CC2H6C H2 (k 2CC2H6 + k 7C H2)(k5′ + k6′) k1CC2H6(k 2CC2H6 + k 7C H2)
=
(k 2CC2H6 + k 7C H2)(k5′ + k6′) k1CC2H6 k5′ + k6′
C•H = =
=
k1k 2k 3CC2H6 2 |det M|
+
k1k 7CC2H6C H2S•C2H5 |det M|
k1k 2k 3CC2H6 2 + k1k 7CC2H6C H2(k 3 + k5′ + k6′ + k 9CC2H4) (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′) k1CC2H6[k 2k 3CC2H6 + k 7C H2(k 3 + k5′ + k6′ + k 9CC2H4)] (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′)
constants k5 and k6 are constants and can be defined as new rate coefficients k5′ and k6′. Ethane is the only source of the intermediates •C2H5 and •H.
These reaction steps T3−1 to T3−6 can be rearranged as a reaction network shown in Figure 5a. k1
k2
The two reaction steps C2H6→•CH3→•C2H5 can be reduced
Λ C2H6 • C2H5
The pathway leading to •C2H5 is C2H6⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→•C2H5. Because this pathway does not pass through •H, C•C2H5 is calculated by C•C2H5 = ΛC2H6•C2H5CC2H6S•H/|det M|. The pathway
Λ C2H6 • C2H5
to a single step, C2H6⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→•C2H5, with ΛC2H6•C2H5 = k1k2CC2H6/k2CC2H6 = k1, and the reduced network is shown in Figure 5b. In the reduced network, there are only two intermediates whose concentrations can be obtained from the algorithm developed in this study. The rate equations of different products can be readily calculated. Table 4 summarizes the calculation procedure and results using the developed algorithm. In Table 4, two pseudo rate coefficients k5′ and k6′ are defined. According to the Bodenstein approximation, the intermediate concentration C•C2H5 does not vary with time and thus is a constant. Therefore, its products with the rate
Λ C2H6 • C2H5
Λ•C2H5• H
leading to •H is C2H6⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→•C2H5⎯⎯⎯⎯⎯⎯⎯⎯→•H. Because this pathway passes through all the intermediates, C•H is calculated by C•H = ΛC2H6•C2H5Λ•C2H5•HCC2H6/|det M|. After substitution, cancellation and collection of the terms, the intermediate concentrations and the rate equations for the end members can be obtained, and the results are shown in Table 4. It is important to note that at low ethane conversion, the ethane cracking rate and all of the product reaction rates are first order with respect to ethane. 12989
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
Table 6. Rate Equations of Ethane Pyrolysis at High Conversion − rC2H6 = k1CC2H6 + k 2CC2H6C•CH3 + k4CC2H6C•H − k6′C•C2H5 k1k 2CC2H6 2(k5′ + k6′ + k 9CC2H4)
=k1CC2H6 +
(k 2CC2H6 + k 7C H2)(k5′ + k6′)
+
k1k4CC2H6 2[k 2k 3CC2H6 + k 7C H2(k 3 + k5′ + k6′ + k 9CC2H4)] (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′) =
k1k5′CC2H6 k5′ + k6′
+
k1k 2CC2 2H6(k5′ + k6′ + k 9CC2H4) (k 2CC2H6 + k 7C H2)(k5′ + k6′)
−
k1k6′CC2H6 k5′ + k6′
+
k1k4CC2H6 2[k 2k 3CC2H6 + k 7C H2(k 3 + k5′ + k6′ + k 9CC2H4)] (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′)
rCH4 = k 2CC2H6C•CH3 + k 7C H2C•CH3 = =
k1CC2H6(k5′ + k6′ + k 9CC2H4)(k 2CC2H6 + k 7C H2) (k 2CC2H6 + k 7C H2)(k5′ + k6′) k1CC2H6(k5′ + k6′ + k 9CC2H4) k5′ + k6′
rC2H4 = k 3C•C2H5 + k6′C•C2H5 − k 8CC2H4C•H − k 9CC2H4C•C2H5 =
k1(k 3 + k6′ − k 9CC2H4)CC2H6 −
k5′ + k6′ k1k 8CC2H6CC2H4[k 2k 3CC2H6 + k 7C H2(k 3 + k5′ + k6′ + k 9CC2H4)] (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′)
rC4H10 = k5′C•C2H5 =
k1k5′CC2H6 k5′ + k6′
rC3H6 = k 9CC2H4C•C2H5 =
k1k 9CC2H6CC2H4 k5′ + k6′
rH2 = k4CC2H6C•H − k 7C H2C•CH3 =
k1k4CC2H6 2[k 2k 3CC2H6 + k 7C H2(k 3 + k5′ + k6′ + k 9CC2H4)] (k 2CC2H6 + k 7C H2)(k4CC2H6 + k 8CC2H4)(k5′ + k6′)
−
k1k 7CC2H6C H2(k5′ + k6′ + k 9CC2H4) (k 2CC2H6 + k 7C H2)(k5′ + k6′)
study is not to discuss the detailed kinetic modeling of ethane thermal cracking; the simplified ethane thermal cracking kinetics serves as an example to illustrate how the rate equations of a chemical reaction network with branches from loops can be readily obtained without assuming which single step is the rate-determining step.
At high ethane conversion when sufficient ethylene and hydrogen are formed, they might react with the free radicals • C2H5 and •H to produce more side products as shown by steps T3−7 to T3−9 in Table 3. In this case, all of the reaction steps in Table 3 are rearranged to a new reaction network as shown in Figure 6. In the network shown in Figure 6, the only source of the intermediates is ethane. Using the same algorithm, the intermediate concentrations can be obtained and the results are shown in Table 5. After the intermediate radical concentrations are obtained, the rate equations for the end members can be readily formulated, and the results are shown in Table 6. In Table 6, if the rate coefficients k7, k8, and k9 are set to be zero, the rate equations for all the end members become identical with those shown in Table 4. A comparison of Table 4 and Table 6 reveals that the ethane thermal cracking kinetics is not as simple as first-order one. Many more products are produced in ethane cracking; the actual reaction kinetics is thus much more complicated than those shown in Table 6. For example, a reaction mechanism consisted of 224 species, and 1722 reaction steps were used to predict the ethane thermal pyrolysis at high conversion.16 The primary purpose of this
5. CONCLUSIONS A general algorithm was developed to derive the rate equations for reactants and products of homogeneous reaction networks with branches from loops. It saves time in deriving the equations for hypothetical networks. Using the algorithm, the pathway rate equations and the intermediate concentrations can be directly obtained; the rate equations for the reactants and products can thus be derived by simple substitutions, cancellations, and collection of terms. For irreversible reactions, the algorithm allows some of the remaining rate equations to be replaced by simpler yield-ratio equations with which rival mechanisms can be discriminated. Moreover, a lower number of ordinary differential equations and lumped kinetic parameters are needed to model the reaction kinetics and thus greatly reduce the efforts for the kinetic data analysis and parameter estimation. With the introduction of pseudo rate 12990
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research
Article
Denote the absolute value of the common denominator −Δ as |det M|; the numerators in eqs A-3 to A-5 should be multiplied by −1 to become positive. For example, the positive numerator for CX1 becomes
coefficients, i.e., the products of the rate coefficients and the intermediate concentrations, the algorithm developed in this study enables us to derive the rate equations for nonsimple networks such as those involved in free radical reactions.
■
APPENDIX Applying the Bodenstein approximation of quasi-stationary behavior of reaction intermediates leads to the following algebraic equations: ⎧−S1C X + Λ 21C X + Λ31C X = −Λ P 1C P 1 2 3 1 1 ⎪ ⎪ ⎪ Λ12C X − S2C X + Λ32C X 1 2 3 ⎨ ⎪ = −Λ P22C P2 − Λ P32C P3 ⎪ ⎪ Λ13C X + Λ 23C X − S3C X = −Λ P 3C P ⎩ 1 2 3 4 4
Λ P11C P1
Λ P43C P4 = Λ P11C P1
(A-1)
(A-2)
Λ P11C P1
=
Δ1 −Δ1 = Δ −Δ
is needed to multiply with the Λ coefficient and end member concentration CP1 leading to X1. Similarly, (i) ΛP22CP2Λ21S3 is the contribution of pathway P2 → X2 → X1. Because this pathway does not pass through intermediate X3, |M12| or S3 is needed to multiply with all of the Λ coefficients and end member concentration CP2 leading to X1. (ii) ΛP22CP2Λ23Λ31 is the contribution of P2 → X2 → X3 → X1 consisting of all of the Λ coefficients and end member concentration CP2 leading to
(A-3)
Λ31
Λ12 −Λ P22C P2 − Λ P32C P3 Λ32 CX2 =
−Λ P43C P4
Λ13
−S3
−S1 Λ 21 Λ31
=
X1. (iii) ΛP32CP3Λ21S3 is the contribution of P3 → X2 → X1. Because this pathway does not pass through intermediate X3, |M12| or S3 is needed to multiply with all of the Λ coefficients and end member concentration CP3 leading to X1. (iv)
Δ2 −Δ2 = Δ −Δ
Λ12 −S2 Λ32
ΛP32CP3Λ23Λ31 is the contribution of P3 → X2 → X3 → X1 consisting of all of the Λ coefficients and end member concentration CP3 leading to X1. (v) ΛP43CP4Λ31S2 is the contribution of P4 → X3 → X1. Because this pathway does not pass through intermediate X2, |M13| or S2 is needed to multiply with all of the Λ coefficients and end member concentration CP4 leading to X1. (vi) ΛP43CP4Λ32Λ21 is the contribution of P4 → X3 → X2 → X1 consisting of all of the Λ coefficients and end member concentration CP4 leading to X1. Using the above rules, the positive numerators for CX2 and
Λ13 Λ 23 −S3 (A-4)
−S1 Λ 21
−Λ P11C P1
Λ12 −S2 −Λ P22C P2 − Λ P32C P3 C X3 =
Λ13 Λ 23
−Λ P43C P4
−S1 Λ 21 Λ31
Λ 23 −S3
Λ 23 −S3
Λ13 Λ 23 −S3 −Λ P11C P1
−S2 Λ32
−S2 Λ32
Λ12 −S2 Λ32
−S1
(A-6)
is the contribution of pathway P1 → X1. Because this pathway does not pass through the other two intermediate X2 and X3, |M11| or
Λ 21 Λ31
−S1 Λ 21 Λ31
+ Λ P22C P2(Λ 21S3 + Λ 23Λ31)
According to eq A-6, the intermediate concentration CX1 can be expressed in terms of the end member concentrations with contributions from all possible pathways:
−Λ P22C P2 − Λ P32C P3 −S2 Λ32 C X1 =
Λ 23 −S3
(Λ31S2 + Λ32Λ 21)
where S1 = Λ1P1 + Λ12 + Λ13 ,S2 = Λ2P2 + Λ2P3 + Λ21 + Λ23, and S3 = Λ3P4 + Λ31 + Λ32. The intermediate concentrations can be solved by Cramer’s rule:
Λ 23 −S3
−S2 Λ32
Λ 23 −S3
+ Λ P32C P3(Λ 21S3 + Λ 23Λ31) + Λ P43C P4
⎤ −Λ P11C P1 ⎡−S1 Λ 21 Λ31 ⎤⎡ C X1 ⎤ ⎡ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎢ Λ12 −S2 Λ32 ⎥⎢C X 2 ⎥ = ⎢−Λ P22C P2 − Λ P32C P3 ⎥ ⎥ ⎢ ⎥⎢ ⎥ ⎢ ⎥⎦ −Λ P43C P4 ⎣ Λ13 Λ 23 −S3 ⎦⎢⎣ C X3 ⎥⎦ ⎢⎣
−Λ P43C P4
Λ P22C P2 + Λ P32C P3 −S2 Λ32
−Δ1 =
or in matrix form:
−Λ P11C P1
Λ 21 Λ31
=
Δ3 −Δ3 = Δ −Δ
Λ12 −S2 Λ32 Λ13 Λ 23 −S3
CX3 become
(A-5) 12991
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
Industrial & Engineering Chemistry Research −S1 −Δ2 =
Λ P11C P1
Article
Λ31
Λ12 Λ P22C P2 + Λ P32C P3 Λ32 Λ P43C P4
Λ13
−S3
= Λ P11C P1(Λ12S3 + Λ13Λ32) + Λ P22C P2 + Λ P32C P3
−S1 Λ31 Λ13 −S3
−S1 Λ31 Λ13 −S3
+ Λ P43C P4(Λ32S1 + Λ31Λ12) (A-7)
−S1 Λ 21 −Δ3 =
Λ P11C P1
Λ12 −S2 Λ P22C P2 + Λ P32C P3
■
Λ P43C P4
Λ13 Λ 23
= Λ P11C P1(Λ13S2 + Λ12Λ 23) + Λ P22C P2 −S1 Λ 21 Λ12 −S2
(A-8)
The rate equations for the end members can be obtained by substituting the intermediate concentrations to the rate equations. For example, rP2 = Λ2P2CX2 − ΛP22CP2 rP2 =
(Λ P11Λ12Λ 2P2C P1 − Λ1P1Λ 21Λ P22C P2)S3 + + +
|det M| Λ P11Λ13Λ32Λ 2P2C P1 − Λ1P1Λ31Λ 23Λ P22C P2 |det M| (Λ P43Λ32Λ 2P2C P4 − Λ3P4Λ 23Λ P22C P2)S1 |det M| Λ P43Λ31Λ12Λ 2P2C P4 − Λ3P4Λ13Λ 21Λ P22C P2 |det M| (Λ P32Λ 2P2C P3 − Λ 2P3Λ P22C P2)
+
−S1 Λ31 Λ13 −S3
|det M|
(A-9)
As shown by eq A-9, the rate equation for P2 consists of five pathway rates, as described in the text.
■
REFERENCES
(1) Hill, C. G. An Introduction to Chemical Engineering Kinetics and Reactor Design; Wiley: New York, 1977. (2) Holland, C. D.; Anthony, R. G. Fundamentals of Chemical Reaction Engineering, 2nd ed.; Prentice-Hall: New York, 1989. (3) Levenspiel, O. Chemical Reaction Engineering, 3rd ed.; Wiley: New York, 1999. (4) Fogler, H. S. Essentials of Chemical Reaction Engineering; PrenticeHall: New York, 2011. (5) Froment, G. F.; Bischoff, K. B.; De Wilde, J. Chemical Reactor Analysis and Design, 3rd ed.; Wiley: New York, 2011. (6) Helfferich, F. G. Kinetics of Homogeneous Multistep Reactions; Elsevier: New York, 2001. (7) Helfferich, F. G. Kinetics of Multistep Reactions, 2nd ed.; Elsevier: New York, 2004. (8) Chern, J.-M. General Rate Equations and Their Applications for Cyclic Reaction Networks: Single-Cycle Systems. Ind. Eng. Chem. Res. 2000, 39, 4100. (9) Helfferich, F. G. Systematic Approach to Elucidation of Multistep Reaction Networks. J. Phys. Chem. 1989, 93, 6676. (10) Chern, J.-M.; Helfferich, F. G. Effective Kinetic Modeling of Multistep Homogeneous Reactions. AIChE J. 1990, 36, 1200. (11) Chang, M.-W.; Chen, T.-S.; Chern, J.-M. Initial Degradation Rate of p-Nitrophenol in Aqueous Solution by Fenton Reaction. Ind. Eng. Chem. Res. 2008, 47, 8533. (12) Chang, M.-W.; Chern, J.-M. Decolorization of Peach Red Azo Dye, HF6 by Fenton Reaction: Initial Rate Analysis. J. Taiwan Inst. Chem. Eng. 2010, 41, 221. (13) Chang, M.-W.; Chung, C.-C.; Chern, J.-M.; Chen, T.-S. Dye Decomposition Kinetics by UV/H2O2: Initial Rate Analysis by Effective Kinetic Modelling Methodology. Chem. Eng. Sci. 2010, 65, 135. (14) Christiansen, J. A. The Elucidation of Reaction Mechanisms by the Method of Intermediates in Quasi-Stationary Concentrations. Adv. Catal. 1953, 6, 311. (15) Rice, F. O.; Herzfeld, K. F. The Thermal Decomposition of Organic Compounds from the Standpoint of Free Radicals. VI. The Mechanism of Some Chain Reactions. J. Am. Chem. Soc. 1934, 56, 284. (16) Xu, C.; Al Shoaibi, A. S.; Wang, C.; Carstensen, H.-H.; Dean, A. M. Kinetic Modeling of Ethane Pyrolysis at High Conversion. J. Phys. Chem. A 2011, 115, 10470.
(Λ 23S1 + Λ 21Λ13) + Λ P32C P3(Λ 23S1 + Λ 21Λ13) + Λ P43C P4
Djk = denominator in rate equation (eq 3) for simple pathway or segment j ↔ k, dimension depends on pathway kij = rate coefficient of step i → j, dimension depends on molecularity k, ki, ki′ = lumped coefficients, dimension depends on network M, Mij = matrices for arbitrary branched networks Pi = network end members (reactants and products) ri = net rate of formation of i, mol m3 s−1 ri→j = rate contribution of step or segment i → j, mol m3 s−1 Rij ≡ ri/rj = yield ratio of products i and j Si = sum of Λ coefficient of all segments originating from node i in branched networks, s−1 Xi = trace-level intermediates λij = pseudo-first order coefficient of step i → j, s−1 Λij = segment coefficient of segment i → j, s−1 Πij = product of Λ coefficients along path i → j, si‑j
AUTHOR INFORMATION
Corresponding Author
*Phone: 886-2-77364674. Fax: 886-2-25861939. E-mail:
[email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The financial support from the National Science Council of Taiwan under grant NSC102-2221-E-036-023 is highly appreciated.
■
NOMENCLATURE A, B, C, P, Q = coreactants and coproducts Ci = concentration of i, mol m−3 12992
dx.doi.org/10.1021/ie501971a | Ind. Eng. Chem. Res. 2014, 53, 12983−12992
