Necessity and Feasibility of 3D Simulations of Steam Cracking

Article pubs.acs.org/IECR

Necessity and Feasibility of 3D Simulations of Steam Cracking Reactors Pieter A. Reyniers, Carl M. Schietekat, David J. Van Cauwenberge, Laurien A. Vandewalle, Kevin M. Van Geem,* and Guy B. Marin Ghent University, Laboratory for Chemical Technology, Technologiepark 914, 9052 Gent, Belgium S Supporting Information *

ABSTRACT: Using detailed kinetic models in computational fluid dynamics (CFD) simulations is extremely challenging because of the large number of species that need to be considered and the stiffness of the associated set of differential equations. The high computational cost associated with using a detailed kinetic network in CFD simulations is why one-dimensional simulations are still used, although this leads to substantial differences compared to reference three-dimensional simulations. Therefore, a methodology was developed that allows one to use detailed single-event microkinetic models in CFD simulations by on the fly application of the pseudo-steady-state assumption to the radical reaction intermediates. Depending on the reaction model size, a speedup factor of more than 50 was obtained compared to the standard ANSYS Fluent routines without losing accuracy. As proof of concept, propane steam cracking in a conventional bare reactor and a helicoidally finned reactor was simulated using a reaction model containing 85 species: 41 radicals and 44 transported species. Next to a drastic speedup of the simulations due to the kinetic network reduction technique, significant differences were observed between the bare and the finned reactor in the three-dimensional simulations. In particular, the ethene selectivity is reduced by 0.20% by application of the helicoidally finned reactor. The one-dimensional simulations were not able to correctly predict the selectivity effect of the different reactor geometry.

1. INTRODUCTION Accurate simulation of reactive fluid flow requires firstprinciples models for both the chemistry and the fluid flow dynamics. For the former only comprehensive kinetic models consisting of elementary reaction steps will yield valid results over a wide range of conditions. In the field of combustion, oxidation, and pyrolysis, the fundamental modeling of gasphase chemistry has thrived over the last two decades driven by an improved knowledge on the occurring reaction families, the availability of rate coefficients,1 and the automation of kinetic model generation through computer codes.2−10 This increase in modeling accuracy is accompanied by a steep increase in the size of kinetic models, with state-of-the-art combustion and pyrolysis kinetic models easily containing several hundreds of species and thousands of reactions11−14 to account for all reaction pathways and intermediate species. Figure 1 displays the size of selected kinetic models for thermal decomposition, oxidation, and combustion processes over the last two decades. It is clear that the number of species and reactions in these kinetic models has increased drastically. Notwithstanding the developments in accurate chemistry modeling, implementation of these detailed kinetic models in computational fluid dynamics (CFD) simulations remains a daunting task because of two specific characteristics of these kinetic models. First, the simulation time increases considerably with increasing size of the kinetic model as each transported species adds an extra continuity equation to the system of differential equations that is to be integrated. Second, these kinetic models show a large spread in interacting time scales associated with species and reactions.15 This spread originates from the distinct difference in reactivity of radical and © XXXX American Chemical Society

Figure 1. Number of reactions as a function of the number of species of selected reaction models for oxidation, pyrolysis, and combustion (after Lu and Law) over the period 1996−2012. Reproduced (adapted) with permission from ref 54. Copyright 2008 Elsevier Ltd.

molecular species and from reactions in partial equilibrium. These two characteristics of the kinetic models, i.e., large size and stiffness, render the combination of multidimensional CFD simulations with these kinetic models very costly in terms of computation, often beyond the capabilities of (super)computers available in most research centers.16 This Received: July 8, 2015 Revised: October 30, 2015 Accepted: November 6, 2015

A

DOI: 10.1021/acs.iecr.5b02477 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research has resulted in two main categories of reactive flow simulations.16 On one hand, zero-dimensional or one-dimensional models using detailed kinetic models are often used for the assessment of the accuracy of (automatically generated) kinetic models and for the simulation of labscale setups or even for industrial reactors if simplification of the fluid flow is possible. On the other hand, advanced engineering tools relying on CFD methods used for the design of complex, industrially relevant geometries usually model the occurring chemistry rather approximately by only a few reactions.17,18 Merging these two categories, i.e., unification of detailed kinetic models and CFD simulations, would allow more accurate design and optimization of chemical reactors in the broad sense. This challenge has been widely recognized in the combustion engineering community,19 and several methods for reducing the computational cost associated with detailed kinetic models have been proposed. These can be roughly divided into three categories: kinetic model reduction methods, storage/ retrieval-based methods, and adaptive chemistry methods. Model reduction decreases the number of reactions and/or species or reduces the number of independent variables. Storage/retrieval methods are based on tabulation strategies for the species rate of formation. Finally, adaptive chemistry methods use several reduced kinetic models, each valid for a subset of simulation space. As the model reduction technique based on the pseudo-steady-state assumption (PSSA) has proven its worth in speeding up the calculation of the chemical source terms,15 it is applied in this work. Also, for steam cracking of hydrocarbons there is strong interest to shift from conventional one-dimensional modeling, the current state of the art, toward more fundamental, multidimensional reactive flow modeling 20 to simulate enhanced three-dimensional reactor technologies and account for nonuniform radial temperature and species profiles. In the present work, steam cracking reactors are modeled using detailed kinetic models of approximately 100 species. To limit the computational cost, the pseudo-steady-state assumption is applied to the free radicals both a priori and on the fly. The methodology is validated with reference simulations, and the obtained reduction in computational cost is quantified. Finally, the application of an optimized helicoidally finned reactor in an industrial steam cracking unit is simulated, and the effect of the reactor geometry on product yields is analyzed. A comparison of the results of the three-dimensional simulations with onedimensional simulations indicates the necessity of 3D simulations for industrial reactor geometries.

2.1.2. Three-Dimensional Reactor Model. The steady-state governing equations of a compressible, reactive, single-phase fluid flow are summarized in Table 1. To be able to explicitly Table 1. Steady-State Governing Equations for Compressible Reactive Flow with Conjugate Heat Transfer continuity equation momentum equations gas-phase energy equation solid-phase energy equation species transport equation

∇·(ρu ̅ ) = 0

eq 1

∇·(ρuu ̅ ̅ ) = −∇p + ∇·τ

eq 2

∇·(u ̅ (ρE + p)) = ∇·(λf ,eff ∇T −

∑ hjJj̅ ) + Sh

eq 3

j

∇·(λs∇T ) = 0

eq 4

∇·(uc̅ j) = −∇·Jj ̅ + R j , ∀ j = 1, nspec − 1

eq 5

account for heat transfer from the metal reactor outer wall to the process gas (conjugate heat transfer), the Laplace conduction equation in the metal surrounding the process gas, eq 4 in Table 1, is simultaneously solved with the governing equations for the fluid flow. In the gas-phase energy equation, E = h − (p/ρ) + (|u|̅ 2/2) with the sensible enthalpy of the ideal gas mixture calculated as T h = ∑jnspec = 1 xjhj and hj = ∫ Tref cp,j(T)dT with cp,j defined as a NASA polynomial as a function of temperature with coefficients for the different species estimated with RMG’s ThermoDataEstimator.5 The effective conductivity of the process gas λf,eff is the sum of the laminar and turbulent thermal conductivity λf,eff = λf,l + λf,t with λf,l calculated from the kinetic theory of gases and λf,t = (cpμt)/(Prt) with Prt = 0.85. Jj̅ is the diffusion flux of species j, including contributions from both the laminar and turbulent diffusivity, i.e., Jj̅ = −(1/MMj)(ρDj,m + (μt/Sct))∇cj with Sct = 0.7. Turbulent diffusion overwhelms laminar and thermal or Soret diffusion in steam cracking reactors.20 Hence, laminar diffusion is modeled in a rather approximate fashion using the kinetic theory of gases and Soret diffusion is neglected. The final term Sh in the gas-phase energy equation is the enthalpy of reaction. The metal thermal conductivity λs is modeled as a temperature polynomial. Rj is the rate of formation of species j, and nspec is the number of species. The renormalization group (RNG) k−ε model was used for the bare tubes, while the Reynolds stress model was used for the helicoidally finned tube as it showed better agreement to experimental data in our previous work.20 The state of the art in the computational fluid dynamics simulation involves implementation of so-called eddy resolving techniques in which turbulent eddies are partially (large eddy simulation, LES) or fully (direct numerical simulation, DNS) resolved, contrary to the most widely used method in which turbulence is fully modeled (Reynoldsaveraged Navier−Stokes, RANS). Explicitly resolving turbulent eddies, albeit partially or fully, is required to adequately model secondary flow structures and their influence on the reactor operation; however, LES and DNS are still prohibitively expensive in terms of computational cost and therefore not easily applicable for systems with complex chemistry such as pyrolysis.26,27 Several applications in combustion28−31 have proven that the combination of eddy resolving techniques with (detailed) chemical kinetics is indeed possible, although the restrictions on the size of the computational domain and the governing Reynolds numbers are still quite severe. To be able to model an industrial-scale reactor at its high-Reynolds

2. NUMERICAL MODELS 2.1. Governing Equations. 2.1.1. One-Dimensional Reactor Model. Instead of performing CPU-intensive simulations with a three-dimensional reactor model, plug flow reactor models are more routinely used to simulate steam cracking reactors. The common method to account for the effect of three-dimensional reactor technologies on product yields and maximum external tube metal temperatures is by multiplication of the friction factor and Nusselt number by a factor derived from nonreactive experiments or simulations as, e.g., reported by Albano et al.21 and Schietekat et al.22 The onedimensional plug flow model CHEMKIN23 using a proprietary stiff solver based on DASPK24,25 was used for comparison with the results of the three-dimensional simulations in case the pseudo-steady-state assumption is not applied. B


Article

Industrial & Engineering Chemistry Research

hydrogen abstractions, β scissions, and C−H scissions, which is described by the remaining 16 subnetworks. During generation of a μ network, the pseudo-steady-state assumption is applied to the μ radicals. In combination with the μ radical hypothesis, this results in a set of algebraic equations that can be solved analytically to the μ radicals’ concentrations. Hence, the concentration of each μ radical can be expressed as a function of the concentration of the starting species of the μ network. Substitution of these concentrations in the rate of formation terms makes the latter only dependent on the concentrations of molecular species and β(μ) radicals. This not only reduces the number of species continuity equations to be solved but also largely reduces the stiffness of the set of differential equations. Consequently, the complete network contains 85 species, of which 41 are radicals and 44 molecules. Validation of the reduced kinetic model for propane cracking was performed by comparison to the full model of Dijkmans et al.33 and is available in the Supporting Information. The thermodynamic data of the C4− molecules and β(μ) radicals was derived from first-principles CBS-QB3 calculations of Sabbe et al.36,37 Thermodynamic data of the μ species was estimated using RMG’s ThermoDataEstimator. The Arrhenius parameters necessary to calculate the reaction rate coefficients were calculated using a group additivity framework that is an extension of Benson’s group additivity.38−40 The reference parameters and group additive values used by Dijkmans et al.33 were also adopted here. Most are based on first-principles calculations, and none were fitted to experimental data. To avoid the computationally costly calculation of the reverse reaction rate coefficient by expressing thermodynamic consistency during the fluid dynamics simulation, the reversible β reactions were split into two irreversible reactions by a priori performing a linear regression of the reverse Arrhenius parameters to the values expressed by thermodynamic consistency in the relevant temperature range, i.e., between 700 and 1300 K. 2.2.2. Pseudo-Steady-State Assumption for β(μ) Radicals. Elimination of the μ radicals by expressing the PSSA during reaction model generation greatly reduces the number of species continuity equations and the system stiffness. However, the difference in time scales between molecules and β radicals is so large that the system remains stiff. This requires very strong under-relaxation of the continuity equations of the β radicals to avoid divergence, making the simulation of industrial reactors too costly in terms of computation. Therefore, the PSSA is applied to the β radicals as well. As the rate of formation of these β radicals contains terms that are second order in the PSS species, a priori analytic calculation of the β radicals concentration as done for the μ radicals is impossible. Lu and Law41 developed a method to obtain analytic solutions even when second-order terms in the PSS species concentrations are present. The nonlinear algebraic equations are first linearized and then analytically solved with a directed graph. A good agreement was seen for the ignition delay time. This methodology was not applied here as the nonlinear terms in the PSS equations are strongly influenced by the higher pressure drop induced by three-dimensional reactor technologies and contribute more than 10% to the rate of consumption for some PSS species, such as methyl and allyl radicals. Consequently, to accurately model the effect of the reactor geometry on product yields, accounting for these terms is necessary. Hence, the concentration of the β radical βj is calculated numerically on the fly, i.e., during the flow

number operating conditions, turbulence is here accounted for via RANS modeling. An adiabatic entrance zone with a length of 1.5 m was added upstream all reactors to ensure fully developed profiles for velocity and turbulence parameters at the reactor inlet. The computational grid was taken from our previous work.20 An explicit grid independence study was carried out for nonreactive flow. For the reactive flow, local grid refinement based on the prevailing gradients was imposed without significant effect on the simulation results. The results in the present work can hence be considered as grid independent. At the inlet boundary, the process gas temperature, mass flow rate, turbulence parameters k and ε, and composition of the process gas were imposed. The turbulence parameters were calculated for a turbulence intensity of 8% and a characteristic length scale of 10% of the tube hydraulic diameter. At the reactor outlet, a constant pressure boundary condition was set. The no-slip boundary condition was set at the tube inner walls. The enhanced wall treatment was used to “bridge” the solution variables in the near-wall cells. This model combines a twolayer model with enhanced wall functions by blending linear and logarithmic laws of the wall. Similar grids as in our previous work were used, ensuring that the near-wall cells are within the viscous sublayer, satisfying the y+ < 1 condition.20 The heat flux profile applied to the outer surface of the reactor was taken from a furnace simulation where the boundary condition applied to the reactor tubes was the industrially measured outer wall temperature profile. The adopted models in the furnace simulation were similar to those of Hu et al.32 The resulting heat flux profile as a function of the axial reactor coordinate is given in the Supporting Information. The commercial CFD code ANSYS Fluent 13.0 was adopted to solve the governing equations. 2.2. Calculation Rate of Production Term. The calculation routine for the rate of production terms Rj in the species continuity eqs (5) was tailored to steam cracking kinetics and implemented in a user-defined function (UDF). 2.2.1. Single-Event Microkinetic Model. The adopted singleevent microkinetic model was generated using the same methodology as applied by Dijkmans et al.33 As the details of the methodology have been explained previously,13,33−35 only a brief summary of the key aspects is given here. The reaction model consists of two parts: the β network and the μ network. The majority of the β network consists of the monomolecular and bimolecular reactions between species with 5 carbon atoms or less. The radicals involved in these reactions are called β radicals. Other reactions involving C6+ β radicals (e.g., the benzyl radical) and the βμ radicals (β character at low temperatures and μ at high temperatures, e.g., the 3-methyl-3pentene-2-yl radical) are also part of the β network. The adopted β network was obtained by reducing the β network of Dijkmans et al.33 to the species relevant for propane cracking, the feedstock studied in this work. This resulted in a β network of 758 reversible elementary reactions between 42 molecules and 41 radicals. The β network is appended with the μ network, which is a collection of independent subnetworks. In the present work, 31 subnetworks were included to account for the formation of aromatics. In 15 subnetworks, a β(μ) radical adds to an olefin to produce a C6+ radical. These radicals are assumed to only undergo monomolecular reactions, i.e., the socalled μ radical hypothesis. Unsaturated μ radicals can undergo ring closure reactions to form cyclic (di)olefins. These cyclic species can subsequently dehydrogenate to form aromatics via C


Article


Table 2. Chemical Lifetime of and the Relative Error on Species j When Applying the PSSA to Species j for a Selection of Radicals at z = 5 m 873 K |Δci

lifetime (s) •

H CH3• C2H3• C2H5• C3H5−1• C3H5−2• C3H5−3• C3H7−1• C3H7−2• C4H9−1• C4H9−2•

8.72 2.69 1.44 4.39 7.15 3.28 3.94 6.01 2.37 6.72 3.96

× × × × × × × × × × ×

1073 K

−8

10 10−5 10−6 10−5 10−7 10−6 10−3 10−7 10−5 10−7 10−7

PSSA

2.01 1.57 −2.96 −3.38 −3.88 −4.90 −3.94 −1.54 1.38 9.30 −2.30

× × × × × × × × × × ×

/ci| 10 10−4 10−5 10−4 10−5 10−5 10−2 10−5 10−4 10−4 10−5

(6)

j

41

Solving eq 6 numerically to cβ can be unstable, while a robust solution method is required as divergence of this algebraic solver would terminate the time-consuming flow simulation. A tailored algebraic solver was used that updates the concentration of the β radical βj in iteration step n + 1 via c βnj+ 1

=

c βnj

R βp j

R βcj

, n = 1...nmax (7)

|Rpβj|

|Rcβj|

The concentration will increase if > and decrease if the opposite holds. Convergence is reached when |Rpβj| = |Rcβj|, i.e., eq 6 expressing the PSSA for β radical βj, holds. The recursive relation is used until the relative change in concentration drops below a certain user-specified threshold value ε for all PSS species or until the maximum number of solver iteration steps nmax is surpassed ⎛ c n+1 − c n βj βj max⎜⎜ n ∀j c βj ⎝

⎞ ⎟ < ε or n = nmax ⎟ ⎠

9.93 2.46 3.09 6.93 2.61 3.45 4.39 2.04 3.68 3.14 1.29

× × × × × × × × × × ×

−9

10 10−7 10−7 10−7 10−8 10−7 10−5 10−8 10−7 10−8 10−8

|ΔciPSSA/ci| 3.21 3.35 −5.05 −2.93 −3.14 −2.15 2.51 −3.37 −3.33 −3.38 −3.38

× × × × × × × × × × ×

10−3 10−3 10−4 10−3 10−3 10−3 10−3 10−3 10−3 10−3 10−3

applying the PSSA for this species, denoted as ΔcPSSA (z), and j the chemical lifetime of each species j, defined as cj(z)/ ∑ivijri(z) at the axial coordinate z. Details regarding the procedure utilized in KINALC to calculate these properties can be found elsewhere.43,44 Two isothermal reactor simulations were performed at 873.15 and 1073.15 K, respectively. In both simulations, the total mass flow rate, steam dilution, and pressure were set to 43.66 × 10−3 kg s−1, 0.325 kg kg−1, and 200 kPa abs, respectively. The reactor inner diameter was 3.02 × 10−2 m. Table 2 shows the results for a selection of important radicals during propane pyrolysis for both simulations. The relative error on the concentration of a radical when applying the PSSA is found to be negligible, below 0.5%. The lifetime of the majority of radicals is below 10−4 s. The lifetime of the allyl radical (C3H5-3•) is about 2−3 orders of magnitude larger than the lifetime of most other β radicals at 873.15 K. This is caused by the resonance stabilization of this radical and results in a larger error when applying the PSSA to the allyl radical. However, at 1073.15 K the effect of the resonance stabilization is lower and the relative error is of the same order of magnitude as for the other β radicals. These results only quantify the error made on a specific radical when applying the PSSA to that radical. In section 3.2, the error on both molecules and radicals is assessed when applying the PSSA on all radicals by direct comparison to simulations without application of the PSSA. 2.2.3. Computation Cost Minimization. As the reaction rates and rates of formation are calculated multiple times during each flow iteration in the solver, the total simulation time depends strongly on the efficiency of these calculations. Therefore, the generation of the UDF file was automated with a Python script. Typically, reaction coefficients are stored in so-called stoichiometric matrices, and calculation of the rates of formation requires several do-loops over the number of species and number of reactions. These large do-loops were eliminated by hard coding the values of all variables needed for the calculation of the reaction rates and the rates of formation. To further reduce the computational cost, the calculations are processed using Maple’s CodeGeneration package to generate an optimal calculation procedure for the rates of formation. This package detects common subexpressions, calculates them once, and stores them for subsequent use to minimize the number of mathematical operations. An obvious example is to only calculate the denominator (RT)−1 once and to store its value instead of calculating it again for every reaction rate coefficient.

simulations by expressing the rate for consumption Rcβj to be equal to the rate of production Rpβj −R βcj(cM , cβ) = R βp(cM , cβ)

lifetime (s)

−5

(8)

A value of 1% for the threshold value ε and 100 for the maximum number of iterations nmax was found to be sufficient. In the first iteration steps of the flow solver, nmax is the limiting factor. However, after several flow iterations, the PSS species concentration profiles change only slightly and as the PSS concentrations of the previous iteration are taken as an initial guess for cβj, convergence is reached rapidly. Applying the PSSA to the β radicals further reduces the number of continuity equations to be integrated. As shown by Ren and Pope, the reactive system is regarded solely to consist of the non-PSS species to guarantee element conservation and to satisfy the first and second law of thermodynamics.42 To quantify the error made by application of the PSSA on the β radicals, the pyrolysis of propane was simulated with the one-dimensional plug flow model implemented in CHEMKIN.23 CHEMKIN uses a proprietary stiff solver based on DASPK,24,25 and the PSSA is not applied to any species. Subsequently, KINALC 43,44 was used to calculate the instantaneous error on the concentration of species j when D


Article


Table 3. Comparison Simulation Results of Three-Dimensional Simulations of Schietekat et al.20 and One-Dimensional Plug Flow Simulations with the 3D Mixing-Cup Temperature and Pressure Profile (T,p) and with the 3D Mixing-Cup Pressure Profile and Heat Flux Profile (Q,p) Scaled to the Same Conversion coil outlet temperature [K] pressure drop [kPa] propane conversion [-] product yields [wt %] H2 CH4 C2H2 C2H4 C2H6 C3H4 (PD) C3H6 C3H8 1,3-C4H6 1-C4H8 2-C4H8 n-C4H10 product mass selectivity [%] H2 CH4 C2H2 C2H4 C2H6 C3H4 (PD) C3H6 1,3-C4H6 1-C4H8 2-C4H8 n-C4H10 valuable olefinsa a

Bare_3D

Bare_1D T, p

Bare_1D Q, p

SmallFins_3D

SmallFins_1D T, p

SmallFins_1D Q, p

1190.9 29.13 84.55

1190.9 29.13 81.99

1192.6 29.13 84.55

1192.7 48.42 85.16

1192.7 48.42 83.50

1192.6 48.42 85.16

1.52 18.75 1.03 38.31 3.65 1.20 18.24 15.45 1.23 0.54 0.02 0.02

1.43 18.13 0.78 35.89 3.21 1.01 19.01 18.01 1.79 0.62 0.04 0.02

1.47 18.90 0.92 37.45 3.34 1.16 18.56 15.45 2.06 0.59 0.04 0.02

1.49 19.25 0.98 38.11 3.46 1.22 18.47 14.84 1.53 0.56 0.03 0.02

1.44 18.69 0.84 36.66 3.28 1.08 18.85 16.50 1.95 0.61 0.04 0.02

1.47 19.21 0.93 37.66 3.35 1.19 18.52 14.84 2.14 0.59 0.04 0.02

1.80 22.18 1.21 45.31 4.31 1.42 21.57 1.46 0.64 0.03 0.02 68.34

1.75 22.11 0.96 43.78 3.92 1.23 23.19 2.19 0.76 0.04 0.02 69.15

1.74 22.35 1.09 44.29 3.95 1.37 21.95 2.43 0.70 0.04 0.02 68.67

1.75 22.61 1.15 44.75 4.07 1.43 21.69 1.79 0.65 0.03 0.02 68.24

1.73 22.38 1.00 43.90 3.92 1.30 22.57 2.33 0.73 0.05 0.02 68.81

1.72 22.56 1.09 44.23 3.94 1.39 21.75 2.51 0.69 0.05 0.02 68.49

Valuable light olefins is defined as the sum of ethene, propene, and 1,3-butadiene.

reactor, respectively. In the plug flow reactor model a species’ rate of formation is calculated using the cross-sectional mixingcup temperature and concentrations. This value is not equal to the mixing-cup average rate of formation, i.e., Rj(C̅ , T̅ ) ≠ R j(C , T ). Consequently, depending on the magnitude of the radial gradient, an error on the predicted product yields is made. Hence, the higher conversion in the 3D simulations is explained by accounting for the high-temperature film near the reactor inner wall. As this high-temperature zone is less distinct in the SmallFins reactor, the agreement between the 1D and the 3D simulations is better. By matching the conversion, the total heat input is overpredicted by 1.07% and 2.59% in the 1D simulations. This shows that when a coupled furnace−reactor simulation45,46 is performed where the fuel flow rate to the furnace burners is adjusted to match a desired cracking severity, e.g., conversion, the simulated fuel flow rate will be overestimated when adopting a one-dimensional reactor model. Furthermore, the effect of the SmallFins reactor on product selectivities compared to the Bare reactor is not well captured by the 1D simulations at the same conversion of the 3D simulations (Q,p cases). For ethene, the 1D simulations show a decrease in selectivity of 0.06 wt % while the 3D simulations show a decrease of 0.56 wt %. In contrast to the increase of 0.12 wt % in propene selectivity simulated with the 3D model, a decrease of 0.20 wt % is seen from the 1D simulations. Finally, the increase of 1,3-butadiene selectivity is underestimated: 0.08 wt

3. RESULTS AND DISCUSSION 3.1. One-Dimensional vs Three-Dimensional Reactor Model. In this section, the 3D simulation results of Schietekat et al.20 are compared to one-dimensional plug flow simulation results to study the necessity of using a three-dimensional reactor model for steam cracking reactors. As the PSSA was not applied in these 3D simulations, the plug flow reactor model implemented in CHEMKIN was used.23 The two reactor configurations referred to as Bare and SmallFins reactor are discussed here. The reactor configuration details are given in the Supporting Information. For both reactors, two onedimensional simulations were performed. In the first type of 1D simulations, referred to with T, the process gas mixing-cup temperature is imposed according to the corresponding profile in the 3D simulation. These simulations correspond to the theoretical case where one would have perfect correlations for the friction and Nusselt numbers resulting in perfect agreement between one-dimensional and three-dimensional temperature and pressure profiles. In the second type of 1D simulation, referred to with Q, the heat flux profile was applied to the reactor inner wall and scaled to obtain the same conversion as in the corresponding 3D simulation. Both simulation types use the mixing-cup pressure profile from the 3D simulations. Table 3 compares the 3D and 1D results and shows that in the 1D simulations with the imposed temperature profile the conversion is lower than in the corresponding 3D simulations. The difference is 2.56 and 1.66 wt % for the Bare and SmallFins E


Article


application of the PSSA on the radicals. Second, the calculation of the reaction rates and rates of formation is optimized. Third, the stiffness of the system is reduced. The optimized calculations and the reduced stiffness allow for more favorable settings of the under-relaxation coefficients in the simulation, further reducing the time to convergence. The numerical contribution to the speedup factor is of the same importance as the PSSA contribution in the case of a small network with minor stiffness. With larger networks and higher stiffness, the numerical contribution to the speedup factor reduces to about 30% of the PSSA contribution to the speedup factor. The speedup factors reported in Table 4 are the result of the combined effect of the PSSA application and the more favorable numerical settings. A detailed explanation on the speedup factor is provided in the Supporting Information. Figure 2 shows the mixing-cup conversion and the mixingcup yields of the most abundant molecules as a function of the axial coordinate in the butane case. A good agreement between the PSSA and the non-PSSA cases is seen for all molecules. Figure 3 shows the mixing-cup yields of several radicals as a function of the axial coordinate in the butane case. Some differences are seen here. Most notably, the radicals that can be formed directly from the feedstock molecule butane have a nonzero value at the reactor inlet, which is most clearly seen for the methyl and ethyl radical. In general, a small underestimation is seen toward the end of the reactor as most clearly seen for the methyl radical. Nonetheless, the agreement is satisfactory. As the species concentrations vary significantly from the reactor centerline to the reactor inner wall20 and as the concentrations near the wall are important for an accurate calculation of the coke formation rate when using a fundamental coke formation model,47−49 the concentrations in the butane case as a function of radial coordinate at an axial position of 9.5 m are compared in Figures 4 and 5 for molecules and radicals, respectively. When applying the PSSA on the β radicals, the conversion near the reactor inner wall at an axial coordinate of 9.5 m is 0.27 wt % higher than in the non-PSSA case. Consequently, the yield of the different products is also slightly higher. A good agreement between the PSSA and the non-PSSA case is seen for the radical concentrations at the center of the reactor. Near the reactor inner wall, the radical concentrations are overpredicted by applying the PSSA. The average relative error is about 10%. For the most important radical coke precursors, i.e., hydrogen, methyl, and ethyl, the relative error is 9%, 15%, and 9%, respectively. In the fundamental coking model of Wauters and Marin,47,48 an increase of the concentration of these radicals by 90% resulted in an increase of the simulated coking rate by

% with the 1D compared to 0.33 wt % with the 3D model. Although these differences are small, they are significant to the economics of an industrial cracker. An estimation for a worldscale propane cracking unit with an annual ethene production rate of 106 metric tons indicates that an error on the predicted yields of 0.1 wt % results in a difference on predicted production of 1000 t/year. Ethene is valued at an average price of around 1000 $/t, so the difference in financial terms equals 106 $/year which is considerable, certainly when taking the lifetime of an industrial unit into account. One-dimensional simulations are not well suited to simulate the selectivity effect of an enhanced reactor geometry as radially averaged species concentrations and temperature do not provide sufficient information to correctly calculate the species’ rate of formation. Hence, an accurate prediction of the effect of a 3D reactor technology on product yields based on a 3D reactor model is necessary upon evaluation of the installation of such a technology. Nonetheless, only detailed single-event microkinetic models guarantee the desired accuracy. As stated above, the implementation of these models in CFD simulations calls for methods that reduce the computational cost such as the PSSA applied here. 3.2. Validation and Speedup. To validate the application of the PSSA on the β radicals and to quantify the obtained speedup, three cases with different β networks were studied. For all cases, a simulation was performed with the developed UDF and with the standard ANSYS Fluent routines for calculating the rates of formation. Simulations were performed for a bare tube on a 2D grid of 280 000 cells. An industrially measured external tube wall temperature profile was imposed on the reactor outer wall. The simulated reactor had a length of 10 m, an internal diameter of 0.032 m, and a wall thickness of 6 mm. The total mass flow rate was set to 0.01 kg s−1 and the inlet temperature to 873.15 K. The steam dilution was 0.333 kg kg−1. The details of the used reaction models and the obtained speedup factor by application of the developed UDF compared to the standard ANSYS Fluent routines are shown in Table 4. Table 4. Specifications of the Reaction Models and the Obtained Speedup Factor in the 3D Validation Simulations feedstock

Nr molecules

Nr radicals

Nr reactions

speedup factor

ethane butane propane

6 8 13

3 6 11

9 27 100

7.4 13.0 54.2

The speedup factor is seen to increase with the model size up to a value of around 50. The reason for the large speedup factor is 3-fold. First, less continuity equations are to be solved by the

Figure 2. Yields of molecules and butane conversion as a function of axial coordinate [m] in 3D validation simulation butane: symbols without PSSA, full lines with PSSA; (A) (■) butane conversion, (○) propene; (B) (●) methane, (□) ethene; (C) (▲) hydrogen, (▽) ethane. F


Article


Figure 3. Yields of radicals as a function of axial coordinate [m] in 3D validation simulation butane: symbols without PSSA, full lines with PSSA; (A) (■) methyl, (○) ethyl; (B) (●) 1-butyl, (□) 2-butyl; (C) (▲) hydrogen radical, (▽) allyl.

Figure 4. Butane conversion and product yields as a function of radial coordinate [×10−3 m] at an axial position of 9.5 m in 3D validation simulation butane: symbols without PSSA, full lines with PSSA; (A) (■) butane conversion, (○) propene; (B) (●) methane, (□) ethene; (C) (▲) hydrogen, (▽) ethane.

Figure 5. Yields of radicals as a function of radial coordinate [×10−3 m] at an axial position of 9.5 m in 3D validation simulation butane: symbols without PSSA, full lines with PSSA; (A) (■) methyl, (○) ethyl; (B) (●) 1-butyl, (□) 2-butyl; (C) (▲) hydrogen radical; (▽) allyl.

kinetic network in order to reduce the computational cost typically associated with industrial-scale reactor simulations. Figure 6 shows the external and internal tube metal temperature as a function of axial coordinate. Maximum values

about 640%, 165%, and 10%, respectively. Hence, it is clear that the error made on the radical concentrations near the wall by application of the PSSA is small but can have a significant influence on the simulated coking rates. 3.3. Simulation of an Industrial Propane Cracking Reactor. An industrial propane cracking reactor was simulated. The reactor is of the millisecond type, i.e., a straight tube with rather small length and diameter compared to other reactor designs. The reactor length is 10.55 m, and the tube crosssectional flow area is 715.72 mm2. A pure propane feedstock is assumed. The steam dilution is 0.325 kg kg−1, and the inlet temperature is 903.7 K. The outlet pressure is 205.7 kPa abs. A heat flux profile as a function of axial coordinate is applied to the reactor outer wall as described in section 2.1.2. The total heat input is adjusted in order to reach the same conversion in the two cases. In this work a conventional bare reactor (Bare) is compared to an optimized finned reactor (SmallFins). The reactor details are given in the Supporting Information. As in our previous work, the cross-sectional flow area, reactor length, and flow rate are the same for the two simulations.20 The pseudo-steady-state assumption was applied to the adopted

Figure 6. External and internal wall temperature [K] as a function of axial coordinate [m] in the 3D simulations of the industrial propane cracking reactor: (red solid line) bare external wall temperature; (red dashed line) bare internal wall temperature: (green dotted line) SmallFins external wall temperature; (green dot dash line) SmallFins internal wall temperature. G


Article

Industrial & Engineering Chemistry Research occur where the heat flux to the reactor is maximal. The maximum external tube temperature is 44 K lower in the SmallFins case compared to the Bare case. Obviously, this will result in a significant increase of the run length. Table 5 shows the 3D simulation results in the columns denoted with Bare 3D and SmallFins 3D. Note that the difference in selectivity in this case is lower as compared to the difference reported in Table 3

at a conversion of 85%. Higher conversions imply a higher heat input to the same reactor and thus stronger radial nonuniformities and higher absolute temperatures, which translate into more substantial differences in selectivity. Figure 7 shows the azimuthally averaged mixing-cup process gas temperature as a function of the radial coordinate at an axial

Table 5. Comparison Simulation Results for Bare and SmallFins Reactor Using a Three-Dimensional (3D) Model and a One-Dimensional Model (1D)

coil outlet temperature [K] pressure drop [kPa] P/E ratio [wt %/wt %] propane conversion [-] product yields [wt %] C4− [C5+-benzene] [benzene-naphthalene] H2 CH4 C2H2 C2H4 C2H6 C3H4 C3H6 C3H8 1,3-C4H6 1-C4H8 2-C4H8 i-C4H8 n-C4H10 1,3-cyclopentadiene Me-1,3-cyclopentadiene benzene toluene styrene naphthalene product mass selectivity [%] C4− [C5+-benzene] [benzene-naphthalene] H2 CH4 C2H2 C2H4 C2H6 C3H4 C3H6 1,3-C4H6 1-C4H8 2-C4H8 i-C4H8 n-C4H10 1,3-cyclopentadiene Me-1,3-cyclopentadiene benzene toluene styrene naphthalene

Bare 3D

Bare 1D

SmallFins 3D

SmallFins 1D

1164 22.9 0.73 74.2

1169 26.7 0.75 74.2

1163 46.8 0.74 74.2

1168 46.7 0.75 74.2

93.83 2.62 3.56 2.03 14.67 0.66 26.72 1.15 0.78 19.46 25.81 1.42 0.95 0.13 0.06 0.01 0.19 0.36 2.32 0.37 0.09 0.78

93.82 2.45 3.74 2.02 14.85 0.66 26.35 1.13 0.77 19.75 25.81 1.40 0.88 0.13 0.06 0.01 1.44 0.32 2.67 0.33 0.08 0.67

93.84 2.60 3.56 2.01 14.77 0.65 26.57 1.13 0.77 19.57 25.80 1.41 0.95 0.13 0.06 0.01 0.19 0.36 2.34 0.36 0.09 0.77

93.82 2.43 3.74 2.01 14.89 0.65 26.30 1.15 0.76 19.80 25.81 1.39 0.88 0.13 0.06 0.01 1.43 0.32 2.68 0.33 0.08 0.66

91.68 3.53 4.79 2.73 19.77 0.89 36.01 1.55 1.05 26.22 1.91 1.28 0.17 0.09 0.01 0.26 0.48 3.12 0.50 0.12 1.05

91.66 3.30 5.04 2.72 20.01 0.89 35.52 1.52 1.04 26.62 1.88 1.18 0.17 0.08 0.01 1.95 0.43 3.59 0.44 0.11 0.90

91.70 3.51 4.79 2.72 19.91 0.87 35.81 1.53 1.03 26.38 1.89 1.28 0.17 0.09 0.01 0.26 0.48 3.15 0.49 0.12 1.04

91.67 3.28 5.05 2.71 20.07 0.87 35.45 1.55 1.02 26.68 1.87 1.18 0.17 0.08 0.01 1.92 0.42 3.61 0.44 0.10 0.89

Figure 7. Azimuthally averaged mixing-cup process gas temperature [K] as a function of the radial coordinate [10−3 m] at an axial coordinate of 10.5 m in the 3D simulations of the industrial propane cracking reactor: (red solid line) Bare; (green dashed line) SmallFins.

coordinate of 10.5 m for the two reactors. The temperature difference between the center and the inner wall is 91 K for the Bare reactor. This is reduced to 65 K by application of the SmallFins reactor. A small bump is seen at a radial coordinate of 15 mm for the SmallFins reactor, which can be attributed to the higher temperature in the wake of the fins.20 Two phenomena influence the product selectivities in the SmallFins reactor compared to the Bare reactor: the larger reactor pressure drop and the more uniform cross-sectional temperature profile. Depending on the reaction path through which a certain product is formed, one of the two phenomena might have a decisive influence. The less uniform crosssectional temperature profile in the Bare reactor will enhance the rate of reactions with a high activation energy. On the other hand, the higher reactor pressure drop in the SmallFins will favor bimolecular over monomolecular reactions. The mixing-cup averaged propane conversion as a function of axial coordinate is shown in Figure 8. The conversion

Figure 8. Mixing-cup averaged propane conversion [%] as a function of axial coordinate [m] in the 3D simulations of the industrial propane cracking reactor: (red solid line) Bare; (green dashed line) SmallFins.

initially increases faster in the Bare reactor as the induction period for heating of the process gas is shorter. Indeed, the high-temperature zone near the reactor inner wall already induces C−C scissions more upstream, resulting in the formation of radicals to start reactions as shown in the inset of Figure 8. Further downstream in the reactor, the mixing-cup averaged temperature profile of the two reactors is very similar and the higher pressure in the SmallFins reactor results in H


Article


A lower ethene selectivity is simulated by application of the SmallFins reactor compared to the Bare reactor. Ethene is mainly formed by the β scission of the 1-propyl radical. This 1propyl radical is formed by hydrogen addition to propene. As shown in Figure 10, the hydrogen radical yield is higher in the

higher reaction rates explaining the higher conversion in the SmallFins reactor. As seen from Table 5, the selectivity toward methane is higher in the SmallFins reactor. Methane is mainly formed by hydrogen abstractions on propane by the methyl radical. Van Geem et al.49 simulated a higher methane yield using a twodimensional simulation model compared to a one-dimensional simulation model during cracking of ethane. This was attributed to the higher concentration of the methyl radical because of the high-temperature zone near the reactor inner wall in the 2D model. Following this reasoning, a higher methane selectivity would be expected in the Bare reactor. During cracking of ethane, the methyl radical is formed by C−C scission of ethane, which has a high activation energy, i.e., about 377 kJ mol−1.50 However, during cracking of propane, methyl is mainly formed by β scission of the 1-propyl radical. This reaction has a much lower activation energy (137 kJ mol−1), and hence, a negligible difference in mixing-cup average concentration of the methyl radical is simulated. As formation of methane from methyl mainly proceeds via bimolecular abstraction, the higher pressure drop in the SmallFins reactor causes the higher methane selectivity.

Figure 10. Hydrogen radical yield [10−7 wt %] as a function of axial coordinate [m] in the 3D simulations of the industrial propane cracking reactor: (red solid line) Bare; (green dashed line) SmallFins.

Bare than in the SmallFins reactor. This higher hydrogen radical yield results in a higher rate of formation of 1-propyl and subsequently a higher ethene production. As this pathway consumes propene, it also explains the higher propene selectivity in the SmallFins reactor. The selectivity toward 1,3-butadiene, 1-butene, and 2-butene is similar in the SmallFins reactor compared to the Bare reactor. Figure 11 shows a scheme of the most important reactions

The selectivity toward ethane is lower in the SmallFins reactor than in the Bare reactor. Ethane is mainly formed by hydrogen abstraction by the ethyl radical and by recombination of two methyl radicals. Figure 9 shows the rate of formation of

Figure 11. Reaction pathway to C4-olefins.

leading to the formation of these components. 1-Butene is formed via the recombination of a methyl and an allyl radical. Subsequent hydrogen abstraction yields the 1-buten-3-yl radical, which can form 2-butene or 1,3-butadiene. As the formation of 1-butene is bimolecular, this reaction is favored at high pressure, explaining the higher yields of the C4-olefins. It is noted that during steam cracking of heavier paraffins, e.g., present in naphtha feedstocks, 1-butene is mainly formed via β scissions, e.g., via β scission of the 3-hexyl radical forming 1butene and ethyl during thermal cracking of hexane.51 When cracking heavier feedstocks such as naphtha these reactions will only by marginally affected when a 3D reactor technology is used. Hence, the yields of C4 olefins are not expected to be drastically affected when using a 3D reactor technology for heavier feedstocks. . The difference in 1,3-butadiene selectivity by application of the SmallFins reactor is smaller than in our previous work.20 This is caused by part of 1,3-butadiene reacting with vinyl radicals to yield benzene as shown by Sabbe et al.52 and Dijkmans et al.33 Another main route to benzene is the addition of allyl to allene and propyne. These pathways are accounted for in the β network by reactions involving benzyl and cyclopentadienyl radicals taken from the high-temperature

Figure 9. Rate of formation [mol m−3 s−1] of ethane in the SmallFins reactor as a function of axial coordinate [m]: (red solid line) by ethyl hydrogen abstraction; (green dashed line) by methyl−methyl recombination.

ethane by these two contributions. The methyl recombination is bimolecular and hence favored at higher pressure. The concentration of the ethyl radical in the Bare reactor is higher than in the SmallFins reactors because of the high-temperature zone near the reactor inner wall. Ethyl radicals are mainly formed via C−C scission of propane having a high activation energy, i.e., about 370 kJ mol−1.50 As the lion’s share of the ethane formation goes through ethyl radicals as shown in Figure 9, the higher ethane formation via the ethyl radical in the Bare reactor compensates for the high-pressure effect on the bimolecular recombination reaction in the SmallFins reactor.

I


Article


Figure 12. Mixing-cup average yield of vinyl radical [10−4 wt %] (A), allene [10−4 wt %] (B), and propadiene [10−4 wt %] (C) as a function of axial coordinate [m] in the 3D simulations of the industrial propane cracking reactor: (red solid line) Bare; (green dashed line) SmallFins.

evaluation of the application of a 3D reactor design in a specific cracker. Table 5 also displays results of one-dimensional simulation in which the heat flux was adjusted to obtain the same conversion. Details about the simulation setup can be found in the Supporting Information. The coil outlet temperature is overpredicted by 5 K for both the Bare case and the SmallFins case. The pressure drop for the Bare case is overpredicted, while the pressure drop for the SmallFins case agrees well with the three-dimensional simulation. The one-dimensional simulations overpredict the selectivity for methane, propene, 1,3cyclopentadiene, and benzene, while they underpredict the selectivity for hydrogen, ethene, and 1,3-butadiene. The differences between the Bare reactor and the SmallFins reactor in terms of selectivity are not well captured by the onedimensional simulations. To mitigate this problem, the onedimensional simulation was extended by accounting for the radial profiles of temperature and radical species. However, this proved also insufficient to capture the geometry effect of the SmallFins reactor compared to the Bare reactor. The results and discussion on these one-dimensional simulations can be found in the Supporting Information.

combustion kinetic model developed by the CRECK modeling group.53 As both routes involve a bimolecular first reaction, they are favored at higher pressure, i.e., in the SmallFins reactor. On the basis of these considerations solely, a much higher increase in benzene selectivity would be expected in the SmallFins reactor. However, both routes to benzene involve typical high-temperature species, i.e., vinyl and allene or propadiene, respectively. Vinyl is formed by hydrogen abstraction from ethene, mainly by methyl and hydrogen radicals, which are hydrogen abstractions with high activation energies (84.4 and 72.1 kJ mol−1, respectively). Allene and propadiene are formed by a hydrogen abstraction from propene to form propen-2-yl and a subsequent β scission. Like the formation of vinyl, the formation of propen-2-yl has a high activation energy. These considerations mean that the vinyl, allene, and propadiene yields are higher in the Bare reactor as shown in Figure 12. Hence, the rate of formation of benzene in the SmallFins reactor in comparison to the Bare reactor is a balance between an increase caused by a higher pressure and a decrease caused by a greater cross-sectional temperature uniformity. This statement is corroborated by Figure 13,

4. CONCLUSIONS Applying the pseudo-steady-state assumption for radical species when using detailed reaction models in CFD simulations allows one to drastically reduce the computational cost by decreasing the number of species continuity equations and by mitigating the stiffness of the resulting mathematical problem. Considering the PSSA also results in significant speedup in 3D simulations. Depending on the kinetic model’s size, a speedup factor from 7.4 to 54.2 was obtained compared to the standard ANSYS Fluent routines. A good agreement was seen for the species yields compared to the conventional 3D simulations. The deviations are largest for the radicals weight fraction near the reactor inner wall with an average relative error of 10%. Applying the methodology for a kinetic model of 85 species to a bare and helicoidally finned propane cracking reactor shows that the finned reactor has a lower coking tendency because the maximum tube temperature was reduced by 44 K. However, at the same time a selectivity loss for ethene of 0.20% was observed, while propene and 1,3-butadiene selectivity changed by 0.16% and −0.02%, respectively. Reaction path analysis on the detailed kinetic model reveals that the changes in ethene and propene selectivity are correlated as less propene is converted to ethene in the SmallFins reactor because of the lower radial temperature gradients and consequently lower hydrogen radical concentrations. The small change in 1,3-

Figure 13. Rate of formation of benzene [mol m−3 s−1] as a function of radial coordinate [10−3 m] at an axial coordinate of 7 m in the 3D simulations of the industrial propane cracking reactor: (red solid line) Bare; (green dashed line) SmallFins.

showing the benzene rate of formation as a function of radial coordinate at an axial coordinate of 7 m. In the reactor center, a higher rate of formation is seen in the SmallFins reactor caused by the higher pressure. However, near the reactor inner wall the rate of formation in the SmallFins reactor is lower as caused by the lower temperature. For the here adopted reactor designs, feedstock, and process conditions, the balance between higher pressure and lower temperature results in similar selectivities. However, for other designs, conditions, and/or feedstocks, one effect might dominate over the other, showing the need for these highly detailed three-dimensional simulations upon J


Article

Industrial & Engineering Chemistry Research nspec Nu p Pr q Q r R ri Sc T u̅ x y y+ z

butadiene and benzene is caused by the higher pressure drop in the finned reactor, favoring bimolecular recombination reactions such as those of methyl and allyl.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.5b02477. Validation of the reduced reaction network against the full single-event microkinetic network from the work of Dijkmans et al.33; description of the simulation details of the industrial reactor simulation; averaging procedures used throughout this work; comparison between the simulation results and reference data from industry; speedup of the 3D simulations when implementing the PSSA; stiffness reduction as a consequence of applying the pseudo-steady-state assumption is verified; details on the computational setup of the 1D simulations implementing PSSA as done with the software package COILSIM1D, also providing an overview and discussion of the results (PDF)

■

Greek Symbols

α correction factor to account for 3D reactor design in 1D model [-] ε threshold value [-] θ azimuthal coordinate [rad] λ thermal conductivity [W m−1 K−1] μ dynamic viscosity [kg m−1 s−1] ρ density [kg m−3] τ stress tensor [Pa]

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

Sub- and Superscripts

Notes

avg c center eff

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS C.M.S. and P.A.R. acknowledge financial support from a doctoral fellowship from the Fund for Scientific Research Flanders (FWO). D.V.C. acknowledges financial support from the combined IWT and BASF Antwerp N.V. Baekeland program. The authors also acknowledge the financial support from the Long Term Structural Methusalem Funding by the Flemish Government, grant number BOF09/01M00409. The computational work was carried out using the STEVIN Supercomputer Infrastructure at Ghent University, funded by Ghent University, the Flemish Supercomputer Center (VSC), the Hercules Foundation, and the Flemish Government− department EWI.

■

f f fr h i in j l m M max p s t β

NOMENCLATURE

Variables

A c Ci+ Ci− cp d E f F h h J̅ L ṁ MM nmax nreac

number of species [-] Nusselt number [-] pressure [Pa] Prandtl number [-] heat flux [W m−2] −1 heat of reaction = ∑j nspec m−3] = 1 Rjhj [J mol radial coordinate [m] rate of formation [mol m−3 s−1] reaction rate of reaction i [mol m−3 s−1] Schmidt number [-] temperature [K] velocity vector [m s−1] mole fraction [molj molf −1] mass fraction [kgj kgf −1] dimensionless wall distance [-] axial coordinate [m]

cross-sectional area [m2] molar concentration [mol m−3] species with i or more carbon atoms [-] species with i or less carbon atoms [-] molar heat capacity [J mol−1 K−1] diameter [m] specific total energy [J kg−1] Fanning friction coefficient [-] molar flow rate [mol s−1] convective heat transfer coefficient [W m−2 K−1] specific enthalpy [J mol−1] diffusion flux [mol m−2 s−1] reactor length [m] mass flux [kg m−2 s−1] molecular mass [kg mol−1] maximum number of iterations [-] number of reactions [-]

mixing-cup average consumption at the centerline of the reactor effective, i.e., sum of laminar and turbulent contributions of the fluid formation friction heat transfer reaction inner tube wall species laminar of the mixture molecule maximum in cross section production of the solid turbulent beta radical

Abbreviations

COT COP PSS PSSA TMT UDF

■

coil outlet temperature [K] coil outlet pressure [Pa] pseudo-steady-state [-] pseudo-steady-state assumption [-] tube metal temperature [K] user-defined function [-]

REFERENCES

(1) Vinu, R.; Broadbelt, L. J.; Prausnitz, J. M. Unraveling Reaction Pathways and Specifying Reaction Kinetics for Complex Systems. Annu. Rev. Chem. Biomol. Eng. 2012, 3, 29−54. (2) Ranzi, E.; Faravelli, T.; Gaffuri, P.; Sogaro, A. Low-temperature combustion: Automatic generation of primary oxidation reactions and lumping procedures. Combust. Flame 1995, 102, 179−192. (3) Côme, G. M.; Warth, V.; Glaude, P. A.; Fournet, R.; BattinLeclerc, F.; Scacchi, G. Computer-aided design of gas-phase oxidation

K


Article

Industrial & Engineering Chemistry Research mechanismsApplication to the modeling of n-heptane and isooctane oxidation. Symp. (Int.) Combust., [Proc.] 1996, 26, 755−762. (4) Broadbelt, L. J.; Stark, S. M.; Klein, M. T. Computer-generated pyrolysis modeling - On the fly generation of species, reactions and rates. Ind. Eng. Chem. Res. 1994, 33, 790−799. (5) Green, W. H.; Allen, J. W.; Ashcraft, R. W.; Beran, G. J.; Goldsmith, C. F.; Harper, M. R.; Jalan, A.; Magoon, G. R.; Matheu, D. M.; Petway, S.; Raman, S.; Sharma, S.; Van Geem, K. M.; Song, J.; Wen, J.; West, R. H.; Wong, A.; Wong, H.-W.; Yelvington, P. E.; Yu, J. RMG-Reaction Mechanism Generator, v3.3. http://rmg.sourceforge.net/ . (6) Vandewiele, N. M.; Van Geem, K. M.; Reyniers, M. F.; Marin, G. B. Genesys: Kinetic model construction using chemo-informatics. Chem. Eng. J. 2012, 207-208, 526−538. (7) Blurock, E. S. Detailed mechanism generation. 1. Generalized reactive properties as reaction class substructures. J. Chem. Inf. Model. 2004, 44, 1336−1347. (8) Blurock, E. S. Detailed mechanism generation. 2. Aldehydes, ketones, and olefins. J. Chem. Inf. Model. 2004, 44, 1348−1357. (9) Rangarajan, S.; Bhan, A.; Daoutidis, P. Language-oriented rulebased reaction network generation and analysis: Description of RING. Comput. Chem. Eng. 2012, 45, 114−123. (10) Rangarajan, S.; Bhan, A.; Daoutidis, P. Language-oriented rulebased reaction network generation and analysis: Applications of RING. Comput. Chem. Eng. 2012, 46, 141−152. (11) Herbinet, O.; Pitz, W. J.; Westbrook, C. K. Detailed chemical kinetic oxidation mechanism for a biodiesel surrogate. Combust. Flame 2008, 154, 507−528. (12) Glaude, P. A.; Herbinet, O.; Bax, S.; Biet, J.; Warth, V.; BattinLeclerc, F. Modeling of the oxidation of methyl esters-Validation for methyl hexanoate, methyl heptanoate, and methyl decanoate in a jetstirred reactor. Combust. Flame 2010, 157, 2035−2050. (13) Van Geem, K. M.; Reyniers, M. F.; Marin, G. B. Challenges of modeling steam cracking of heavy feedstocks. Oil Gas Sci. Technol. 2008, 63, 79−94. (14) van Goethem, M. W. M.; Kleinendorst, F. I.; van Leeuwen, C.; van Velzen, N. Equation-based SPYRO® model and solver for the simulation of the steam cracking process. Comput. Chem. Eng. 2001, 25, 905−911. (15) Dijkmans, T.; Schietekat, C. M.; Van Geem, K. M.; Marin, G. B. GPU based simulation of reactive mixtures with detailed chemistry in combination with tabulation and an analytical Jacobian. Comput. Chem. Eng. 2014, 71, 521−531. (16) Tosatto, L.; Bennett, B. A. V.; Smooke, M. D. A transport-fluxbased directed relation graph method for the spatially inhomogeneous instantaneous reduction of chemical kinetic mechanisms. Combust. Flame 2011, 158, 820−835. (17) Boileau, M.; Staffelbach, G.; Cuenot, B.; Poinsot, T.; Berat, C. LES of an ignition sequence in a gas turbine engine. Combust. Flame 2008, 154, 2−22. (18) Hu, G.; Wang, H.; Qian, F.; Zhang, Y.; Li, J.; Van Geem, K. M.; Marin, G. B. Comprehensive CFD Simulation of Product Yields and Coking Rates for a Floor- and Wall-Fired Naphtha Cracking Furnace. Ind. Eng. Chem. Res. 2011, 50, 13672−13685. (19) Liang, L.; Stevens, J. G.; Farrell, J. T. A Dynamic Multi-Zone Partitioning Scheme for Solving Detailed Chemical Kinetics in Reactive Flow Computations. Combust. Sci. Technol. 2009, 181, 1345−1371. (20) Schietekat, C. M.; Van Cauwenberge, D. J.; Van Geem, K. M.; Marin, G. B. Computational Fluid Dynamics-Based Design of Finned Steam Cracking Reactors. AIChE J. 2014, 60, 794−808. (21) Albano, J. V.; Sundaram, K. M.; Maddock, M. J. Applications of Extended Surfaces in Pyrolysis Coils. Energy Prog. 1988, 8, 160−168. (22) Schietekat, C. M.; van Goethem, M. W. M.; Van Geem, K. M.; Marin, G. B. Swirl flow tube reactor technology: An experimental and computational fluid dynamics study. Chem. Eng. J. 2014, 238, 56−65. (23) Kee, R. J.; Rupley, F. M.; Miller, J. A.; Coltrin, M. E.; Grcar, J. F.; Meeks, E.; Moffat, H. K.; Lutz, G.; Dixon-Lewis, A. E.; Smooke, M. D.; Warnatz, J.; Evans, G. H.; Larson, R. S.; Mitchell, R. E.; Petzhold, L.

R.; Reynolds, W. C.; Caracotsios, M.; Stewart, W. E.; Glarborg, P.; Wang, C.; Adigun, O.; Houf, W. G.; Chou, C. P.; Miller, S. F.; Ho, P.; Young, D. J. CHEMKIN, Release 4.1.1; Reaction Design, Inc.: San Diego, CAA, 2007. (24) Li, S.; Petzold, L. Design of New DASPK for Sensitivity Analysis, UCSB Technical report. 1999. (25) Shengtai, L.; Linda, P. Design of new Daspk for Sensitivity Analysis; University of California at Santa Barbara: Santa Barbara, 1999. (26) Van Cauwenberge, D. J.; Schietekat, C. M.; Floré, J.; Van Geem, K. M.; Marin, G. B. CFD-based design of 3D pyrolysis reactors: RANS vs. LES. Chem. Eng. J. 2015, 282, 66. (27) Zhu, M. Large eddy simulation of thermal cracking in petroleum industry; Institut National Polytechnique de Toulouse: Toulouse, 2015. (28) Mahesh, K.; Constantinescu, G.; Apte, S.; Iaccarino, G.; Ham, F.; Moin, P. Large-eddy simulation of reacting turbulent flows in complex geometries. J. Appl. Mech. 2006, 73, 374−381. (29) Walter, M.; Kornev, N.; Hassel, E. Large Eddy Simulation of Turbulent Reactive Mixing at High Schmidt and Reynolds Numbers. Chem. Eng. Technol. 2015, 38, 1608−1616. (30) Coelho, P. J. Approximate solutions of the filtered radiative transfer equation in large eddy simulations of turbulent reactive flows. Combust. Flame 2009, 156, 1099−1110. (31) Di Sarli, V.; Di Benedetto, A. Sensitivity to the Presence of the Combustion Submodel for Large Eddy Simulation of Transient Premixed Flame-Vortex Interactions. Ind. Eng. Chem. Res. 2012, 51, 7704−7712. (32) Hu, G.; Wang, H.; Qian, F.; Van Geem, K. M.; Schietekat, C. M.; Marin, G. B. Coupled simulation of an industrial naphtha cracking furnace equipped with long-flame and radiation burners. Comput. Chem. Eng. 2012, 38, 24−34. (33) Dijkmans, T.; Pyl, S. P.; Reyniers, M.-F.; Abhari, R.; Van Geem, K. M.; Marin, G. B. Production of bio-ethene and propene: alternatives for bulk chemicals and polymers. Green Chem. 2013, 15, 3064−3076. (34) Clymans, P. J.; Froment, G. F. Computer-generation of reaction paths and rate-equations in the thermal-cracking of normal and branched paraffins. Comput. Chem. Eng. 1984, 8, 137−142. (35) Hillewaert, L. P.; Dierickx, J. L.; Froment, G. F. ComputerGeneration of Reaction Schemes and Rate-Equations for Thermal Cracking. AIChE J. 1988, 34, 17−24. (36) Sabbe, M. K.; De Vleeschouwer, F.; Reyniers, M. F.; Waroquier, M.; Marin, G. B. First Principles Based Group Additive Values for the Gas Phase Standard Entropy and Heat Capacity of Hydrocarbons and Hydrocarbon Radicals. J. Phys. Chem. A 2008, 112, 12235−12251. (37) Sabbe, M. K.; Saeys, M.; Reyniers, M. F.; Marin, G. B.; Van Speybroeck, V.; Waroquier, M. Group additive values for the gas phase standard enthalpy of formation of hydrocarbons and hydrocarbon radicals. J. Phys. Chem. A 2005, 109, 7466−7480. (38) Sabbe, M. K.; Reyniers, M. F.; Van Speybroeck, V.; Waroquier, M.; Marin, G. B. Carbon-centered radical addition and beta-scission reactions: modeling of activation energies and pre-exponential factors. ChemPhysChem 2008, 9, 124−40. (39) Sabbe, M. K.; Reyniers, M. F.; Waroquier, M.; Marin, G. B. Hydrogen radical additions to unsaturated hydrocarbons and the reverse beta-scission reactions: modeling of activation energies and pre-exponential factors. ChemPhysChem 2010, 11, 195−210. (40) Sabbe, M. K.; Vandeputte, A. G.; Reyniers, M. F.; Waroquier, M.; Marin, G. B. Modeling the influence of resonance stabilization on the kinetics of hydrogen abstractions. Phys. Chem. Chem. Phys. 2010, 12, 1278−98. (41) Lu; Law, C. K. Systematic Approach To Obtain Analytic Solutions of Quasi Steady State Species in Reduced Mechanisms. J. Phys. Chem. A 2006, 110, 13202−13208. (42) Ren, Z. Y.; Pope, S. B. Entropy production and element conservation in the quasi-steady-state approximation. Combust. Flame 2004, 137, 251−254. L


Article

Industrial & Engineering Chemistry Research (43) Kovacs, T.; Zsely, I. G.; Kramarics, A.; Turanyi, T. Kinetic analysis of mechanisms of complex pyrolytic reactions. J. Anal. Appl. Pyrolysis 2007, 79, 252−258. (44) Turanyi, T.; Tomlin, A. S.; Pilling, M. J. On the error of the quasi-steady-state approximation. J. Phys. Chem. 1993, 97, 163−172. (45) Plehiers, P. M.; Reyniers, G. C.; Froment, G. F. Simulation of the run length of an ethane cracking furnace. Ind. Eng. Chem. Res. 1990, 29, 636−641. (46) Heynderickx, G. J.; Froment, G. F. Simulation and Comparison of the Run Length of an Ethane Cracking Furnace with Reactor Tubes of Circular and Elliptical Cross Sections. Ind. Eng. Chem. Res. 1998, 37, 914−922. (47) Wauters, S.; Marin, G. B. Kinetic Modeling of Coke Formation during Steam Cracking. Ind. Eng. Chem. Res. 2002, 41, 2379−2391. (48) Wauters, S.; Marin, G. B. Computer generation of a network of elementary steps for coke formation during the thermal cracking of hydrocarbons. Chem. Eng. J. 2001, 82, 267−279. (49) Van Geem, K. M.; Heynderickx, G. J.; Marin, G. B. Effect of radial temperature profiles on yields in steam cracking. AIChE J. 2004, 50, 173−183. (50) Luo, Y.-R. Handbook of bond dissociation energies in organic compounds; CRC Press: Boca Raton, FL, 2003. (51) Van Geem, K. M.; Reyniers, M. F.; Marin, G. B.; Song, J.; Green, W. H.; Matheu, D. M. Automatic reaction network generation using RMG for steam cracking of n-hexane. AIChE J. 2006, 52, 718−730. (52) Sabbe, M. K.; Van Geem, K. M.; Reyniers, M.-F.; Marin, G. B. First principle-based simulation of ethane steam cracking. AIChE J. 2011, 57, 482−496. (53) Cavallotti, C.; Polino, D.; Frassoldati, A.; Ranzi, E. Analysis of Some Reaction Pathways Active during Cyclopentadiene Pyrolysis. J. Phys. Chem. A 2012, 116, 3313−3324. (54) Lu, T.; Law, C. K. Toward accommodating realistic fuel chemistry in large-scale computations. Prog. Energy Combust. Sci. 2009, 35, 192−215.

■

NOTE ADDED AFTER ASAP PUBLICATION After this paper was published ASAP November 23, 2015, several corrections were made to the text. The corrected version was reposted December 1, 2015.

M


Necessity and Feasibility of 3D Simulations of Steam Cracking

Recommend Documents