Impact of Radiation Models in Coupled Simulations of Steam

Subscriber access provided by UNIVERSITY OF MICHIGAN LIBRARY

Article

Impact of radiation models in coupled simulations of steam cracking furnaces and reactors Guihua Hu, Carl Schietekat, Yu Zhang, Feng QIAN, Geraldine J. Heynderickx, Kevin Marcel Van Geem, and Guy B Marin Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie5042337 • Publication Date (Web): 13 Feb 2015 Downloaded from http://pubs.acs.org on February 18, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Impact of radiation models in coupled simulations of steam cracking furnaces and reactors Guihua Hu1, Carl M.Schietekat2, Yu Zhang1,2, FengQian1,*, Geraldine Heynderickx2, Kevin M. Van Geem2,*,Guy B. Marin2 1

Key Laboratory of Advanced Control and Optimization for Chemical Processes of Ministry of

Education, East China University of Science and Technology, Shanghai 200237, China. 2

Laboratory for Chemical Technology, Ghent University, Technologiepark 914, B-9052Ghent,

Belgium. KEYWORDS: radiation model; P-1; DOM; DTRM; steam cracking; furnace; combustion; PDF ABSTRACT: As large floor-fired furnaces have many applications in refinery and (petro-) chemical units and about 80% of heat transfer in these furnaces is by radiation, the accurate description of radiative heat transfer is of the most importance for accurate design and optimization. However, the impact of using different radiation models in coupled furnace/reactor simulations has never been evaluated before. Therefore coupled furnace/reactor simulations of an industrial naphtha cracking furnace with a 130kt/a capacity have been conducted. Computational fluid dynamics simulations were performed for the furnace side, while the onedimensional reactor model COILSIM1D was used for the reactor simulations. The Adiabatic, P1, discrete ordinates model (DOM) and discrete transfer radiation model (DTRM) were evaluated for modeling the radiative heat transfer. The results with DOM and the DTRM radiation model are very similar both on the furnace and the reactor side. The flue gas temperature using DOM is higher than when using the P-1 radiation model, resulting in higher incident radiation. Comparing the simulated results of all radiation models to the industrial product yields and run lengths shows that DOM and DTRM outperform the others. As DOM has a broader application range than DTRM, and because the current implementation of DTRM in

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 49

FLUENT/14.0 cannot be run in parallel yet, DOM is the recommended radiation model for run length simulations of steam cracking furnaces. 1. INTRODUCTION Large floor-fired furnaces have many applications in refinery and (petro-) chemical units. One of the most important applications is in the hot section of steam cracking units for the production of olefins and aromatics. In this case several tubular reactors are suspended in the furnace and the heat released by the burners in the furnace is transferred to the reactor tubes by convection, radiation and conduction. This heat is on the one hand required for the production of steam and evaporation of the liquid feed. More importantly the heat is needed to drive the endothermic thermal cracking reactions and to overcome the insulating effect of the forming cokes layer on the Fe-Cr-Ni heat resistant steel reactors.1, 2The total heat absorbed by the reactor coils in the radiation section of a steam cracking furnace is about 42-47% for 100% floor firing.3 This value, which is also known as furnace thermal efficiency, consists of both radiative heat transfer and convective heat transfer. Because of the high temperature in the furnace, i.e. above 1300 K, radiation dominates the heat transfer process.4, 5Hence, an accurate prediction of radiative heat transfer is a prerequisite for a correct simulation of these furnaces. In the past decades, many researchers have used the Lobo-Evans method,6 the Belokon’s method7 and zone methods4,

8, 9

for simulating industrial steam cracking furnaces. The fuel

combustion was not simulated rigorously, instead a predefined heat release rate was imposed to estimate the composition and temperature of the flue gas. Furthermore, convective heat transfer to the reactor tubes was often ignored. These simplifications obviously cause a certain error, which may result in inaccurate design optimization. By virtue of the development of accurate models and continuously growing computational power, computational fluid dynamics (CFD) has steadily grown to become an important and indispensable simulation tool for the chemical industry. More particularly for the simulation of 2 ACS Paragon Plus Environment

Page 3 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


steam cracking furnaces, different CFD models have been evaluated by many researchers over the last two decades. Wang and Zhang10 used the P-1 radiation model to calculate radiative heat transfer. Zhou and Jia11 adopted the discrete ordinates model (DOM), but they introduced empirical formulas for the calculation of the flue gas radiative properties, which could introduce large errors in the results. Coelho12 and Stefanidis et al.13 assessed the influence of adopting nongray radiative properties of the flue gas mixture. More recently, Hu et al. ,14, 15 Yang et al.16 and Hassan et al.17 performed coupled simulations of the furnace and the reactor tubes, in which DOM was applied in the furnace simulation. However, in all these studies only a single radiation model was applied and the results can therefore not be used for the comparison of the performance of various radiation models. Keramida et al.18 compared the discrete ordinates and six-flux radiation model for a natural gas diffusion flame and concluded that the two models performed similarly, both showing good agreement with the experimental data. Li et al.19 compared different radiation models for heat transfer in a vertical pipe. Mendes et al.20 adopted DOM and Rosseland model for the onedimensional simulation of a premixed flame in inert porous media. Zhang et al.21 looked at the influence of the flue gas radiative properties and burner geometry on the flame front in the firebox, the heat transfer to the coils and the product selectivities has been investigated. To the authors’ knowledge only Habibi et al.22, 23 evaluated different radiation models in steam cracking furnaces. They assessed the impact of the radiation model on the temperature distribution in the furnace, the heat flux to the reactors, the species profiles in the furnace and the structure of the burner flames. In the first paper,22 a simplified reaction network and the Finite Rate/EddyDissipation model were used, while in the second paper23 a more complex reaction network and the eddy dissipation concept (EDC) model were used. These model differences resulted in different heat flux profiles for the different radiation models, that is, in the first paper a lower heat flux to the reactors was simulated with the P-1 model compared to DOM while in the



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 49

second paper the contrary was simulated. The main conclusion of both papers was that the differences between using DOM and P-1 are negligible. However, given the broader application range of DOM, use of this model was recommended, notwithstanding the model’s higher CPU time requirements compared to the P-1 model. Habibi et al. did not perform a coupled simulation of the furnace and the reactor tubes, instead a fixed process gas temperature profile was imposed from which the heat flux profile to the reactors and the tube external wall temperature profile was calculated. As the heat transfer to the reactors influences the heating of the process gas and the cracking reactions, the process gas temperature profile will change drastically when altering the radiation model, which will have an effect on the furnace simulation. The latter can only be addressed via a series of coupled furnace/reactor simulations. Furthermore, in the work of Habibi et al.,22, 23 it is not even assured that the heat balance over the entire system, i.e. furnace and reactors is closed. Therefore in this work coupled furnace/reactor simulations have been performed for an industrial naphtha cracking furnace with 130kt/a ethene capacity. In the furnace model, the compressible formulation of the Reynolds-averaged Navier-Stokes (RANS) equations is adopted to simulate the flue gas flow. The standard k-ε model is used for closure. As non-premixed burners are used, the probability density function (PDF) model24 is used for the simulation of the fuel combustion. The Adiabatic, P-1, DOM and discrete transfer radiation model (DTRM) radiation models are used for modeling of the radiative heat transfer. A weighted-sum-of-graygases model (WSGGM)14, 15 is used to calculate the radiative characteristics of the flue gas. For the reactor simulations, a plug flow model combined with a single-event microkinetic freeradical model is adopted. The influence of the adopted radiation model on the flue gas flow, radiative properties, product yields, heat flux profiles to the reactor tubes, external tube metal and process gas temperature profiles and run length is analyzed. Accordingly, the best radiation model for the simulation of steam cracking furnaces is determined.


Page 5 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


2. NUMERICAL MODELS 2.1. Flow Model. As the flow of the flue gas is turbulent and complex combustion and heat transfer processes occur in the furnace, the mass, momentum and energy transport equations need to be solved. The RANS equations are adopted to describe the three-dimensional compressible turbulent flow. The equations are closed using the standard k-ε model.25 The general form of the resulting transport equations can be written as:

∂( ρU jφ ) ∂x j

=

∂  ∂φ  Γφ + Sφ ∂x j  ∂x j 

(1)

2.2. Combustion Model. In the simulated furnace, only floor burners are used, in which the fuel gas and air enter the furnace separately, i.e. non-premixed combustion. Therefore, the probability density function (PDF) model is used to describe the turbulent combustion process. The PDF model assumes that the chemical reaction rate is infinitely fast and irreversible. Fuel and oxidant never coexist in space as chemical reactions occur immediately at contact and complete in a one-step conversion to the final products, i.e. carbon dioxide and steam. Therefore, the instantaneous thermochemical state of the flue gas in the furnace can be expressed by a mixture fraction:

f =

Z i − Z i ,ox

(2)

Z i , fuel − Z i ,ox

where Zi is the elemental mass fraction for element i. The subscript ox denotes the value at the oxidizer stream inlet and the subscript fuel denotes the value at the fuel stream inlet. As the assumption of equal diffusivities is generally acceptable for turbulent flows where the turbulent convection overwhelms molecular diffusion,26 the species transport equations can be replaced by a transport equation of the mean mixture fraction: ∂ ∂ µt ∂ f ( ρU i f ) = ( ) + Sm ∂xi ∂xi σ t ∂xi

(3)

In addition the mean mixture fraction variance obeys the following conservation equation: 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

∂f ∂ ∂ µt ∂ f ′2 ε ( ρU i f ′ 2 ) = ( ) + C g µt ( ) 2 − Cd ρ f ′2 ∂x i ∂x i σ t ∂x i ∂x i k

Page 6 of 49

(4)

where f ′ = f − f .The constants σ t , C g and Cd are 0.85, 2.86 and 2.0 respectively. By solving the conservation equations of mean mixture fraction and its variance, the mass fraction of each species in the flue gas, e.g. carbon monoxide, can be obtained based on a predefined PDF table. This predefined PDF table was generated based on the assumption that chemical equilibrium is reached as soon as the fuel and oxidant are mixed. The following table parameters were used: the number of mean mixture fraction points (21), the number of mixture fraction variance points (11), the maximum number of species (20) and the number of mean enthalpy points (45). The main advantage of the PDF approach is that it reduces the number of transport equations from the number of species that are involved in the chemical reaction to two, i.e. only the transport equations for mean mixture fraction and mixture variance need to be solved. Hence, the simulations are greatly accelerated. Nevertheless, the predicted combustion rate with the PDF approach could be higher than the actual combustion rate, as was shown by Stefanidis et al.27 However, the aim of this work is to evaluate the impact of different radiation models on heat transfer and reaction rates in a coupled steam cracking furnace/reactor simulation. The impact of using detailed combustion schemes is in this case minor based on the work of Stefanidis et al.,27 and hence, does not need to be accounted for in this work. For future research, the Eddy Dissipation Concept (EDC) model is suggested as a more accurate combustion model. It assumes that the molecular mixing and the subsequent reactions are taken place in the small turbulent structures that are in the order of Kolmogorov scale, in which kinetic energy is dissipated into heat. These so-called fine structures are treated locally as adiabatic Perfectly Stirred Reactors (PSR) under constant pressure, in which the reaction rates are calculated using finite-rate chemical kinetics as well as the local temperature and species concentration. By introducing this more sophisticated and reasonable turbulence-chemistry interaction assumption, the model


Page 7 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


avoids overprediction of combustion rate when the mixed-is burnt approximation is not satisfied.27, 28 2.3. Radiation Model. As stated before, it is necessary to adopt an accurate radiation model, because radiative heat transfer dominates the heat flux to the reactor tubes. For a specific problem, one radiation model may be more appropriate than another depending on the optical thickness of the simulated medium. Optical thickness is a dimensionless value given by:

τ = αL

(5)

where α is the medium absorption coefficient and L is a characteristic length for the simulated geometry. Several radiation models have been developed and implemented in CFD packages over the last decades including the P-1 model, Rosseland model, discrete-ordinates model (DOM), discrete transfer radiation model (DTRM) and surface-to-surface (S2S) radiation model. The main assumption of the S2S model is that any absorption, emission, or scattering of radiation can be ignored. Consequently, only "surface-to-surface'' radiation needs to be considered. The S2S model is suitable for enclosure radiative transfer with non-participating media and is thus not applicable to the flue gas medium of the cracking furnace. The Rosseland model is valid for an optical thickness larger than 3 whereas the optical thickness in steam cracking furnaces is below 1 as will be shown below. Hence, the Rosseland model is not valid for the simulation of steam cracking furnaces as was confirmed by Habibi et al.22 Therefore DOM, DTRM and P-1 model were selected for evaluation in this work. A simulation without radiative heat transfer was performed, referred to as the Adiabatic model. As excess air is provided and combustion is complete, almost no soot particles are formed and the scattering coefficient can be assumed to be negligible compared to the absorption coefficient. Hence, the considered radiation models are all applied to the radiative transfer equation (RTE) for an absorbing, non-scattering and emitting medium at position in the direction given by:



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

r r r r σT 4 dI ( r , s ) + αI ( r , s ) = α ds π

Page 8 of 49

(6)

Assuming that the refractive index is equal to unity. 2.3.1.Discrete Ordinates Model. DOM solves the RTE for a finite number of discrete solid

r angles, each associated with a vector direction s fixed in the global Cartesian system.DOM transforms equation (6) into a transport equation for radiation intensity in the spatial Cartesian coordinates(x, y, z). As many transport equations as there are directions are solved.The solution method of these equations is identical to that used for the flow and energy equations.DOM

r considers the RTE in the direction s as a field equation and can be written as: r r r r r σT 4 ∇ ⋅ ( I (r , s ) s ) + αI ( r , s ) = α

(7)

π

The wall surface temperature is calculated from the following equation,29 while imposing qout : q out = (1 − ε w ) q in + ε w σ T w4

(8)

where 4π r r qin = ∫r r I in s ⋅ ndΩ

(9)

s ⋅n > 0

r The boundary intensity for all outgoing directions s at the wall is given by: I0 =

qout

(10)

π

2.3.2. P-1 Model. The P-1 radiation model is the simplest form of the more general P-N model, which is based on the expansion of the radiation intensity I into an orthogonal series of spherical harmonics.30 If only four terms in the series are used, the following equation is obtained for the radiative flux qr:

q r = − Γ ∇G

(11)

where


Page 9 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Γ=

1 3α

(12)

Hence, the problem is simplified to finding a solution for the incident radiation G. The transport equation for G is: ∇ ⋅ (Γ ∇ G ) − α G + 4ασ T 4 = S G

(13)

Combining equations (11) (11)and (13), the following equation is obtained:

− ∇ ⋅ q r = αG − 4ασT 4

(14)

Therefore, this expression for − ∇ ⋅ qr can be directly substituted into the energy equation (1) to account for heat sources due to radiation. 2.3.3. Discrete Transfer Radiation Model. The main assumption of the discrete transfer radiation model (DTRM) is that radiation leaving a surface in a certain range of solid angles can be approximated by a single ray. The change of radiant intensity, dI, along a path, ds, can be written as:

ασT 4 dI + αI = ds π

(15)

The DTRM integrates equation (15) along a series of rays emanating from the boundary walls. If α is constant along the ray, then I(s) can be estimated as

I (s ) =

σT 4 ( 1 − e −αs ) + I 0 e −αs π

(16)

Therefore, the energy source due to radiation in equation (1) (1)is computed by summing the change in intensity along the path of each ray that is traced through the fluid control volume. 2.4. Reactor and Feedstock Model. As the reaction families dictating the chemistry in steam cracking reactors have been well-known for many years, fundamental process simulation tools for the simulation of steam cracking reactors have been used extensively since the pioneering work of Dente et al.31 in the late seventies. Although several authors have used the twodimensional32,

33

and three-dimensional15,

24, 34-37

reactor models for the simulation of steam 9

ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 49

cracking reactors, the one-dimensional plug flow model for steam cracking is mostly adopted. In this work, the one-dimensional plug flow reactor model COILSIM1D was used to simulate the reactors.38 As steam cracking is a non-isobaric and non-isothermal process, the species continuity equations are solved together with the energy and momentum equation. A singleevent microkinetic model containing thousands of reactions between over 700 species is used to simulate the steam cracking free-radical chemistry. A more detailed description of the reactor model and kinetic network can be found in previous works.39, 40 As the adopted reactor model uses a detailed single-event microkinetic model, a detailed molecular feedstock specification is necessary. However, such a detailed feedstock analysis, e.g. obtained by comprehensive twodimensional gas chromatography,41 is often lacking for industrial units. Typically only a limited number of feedstock characteristics, so-called commercial indices, obtained via standardized ASTM procedures are available. For the feedstock of the simulated industrial unit, the density, average molecular mass, global PONA analysis and ASTM 2887 boiling point analysis were available. To reconstruct the detailed molecular feedstock from these characteristics, the simulation package SimCo was used which uses an artificial neural network (ANN) trained for a set of 300 reference naphthas and a Shannan entropy method.42 As the considered naphtha lies outside the application range of the ANN, the Shannon entropy method was adopted. 2.5. Thermophysical and Radiative Properties. The ideal gas equation of state is used. The flue gas density is calculated in each cell from the local pressure, temperature and species mass fraction:

ρ=

p RT ∑

Yi Mi

(17)

The dynamic viscosity and thermal conductivity of the flue gas mixture are calculated from the local temperature as:24


Page 11 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


µ = 7.61×10−6 + 3.26 ×10−8 × T

(18)

λ = 7.67 ×10 −3 + 5.88 ×10−5 × T

(19)

The heat capacity of the flue gas mixture is computed as the mass-weighted average of the heat capacities of the pure species as: c p = ∑ Yi c p ,i

(20)

i

where cp,i is the heat capacity of pure species i calculated from a polynomial function of temperature: c p ,i = Ai + BiT + C iT 2 + DiT 3

(21)

where Ai, Bi, Ci and Di are constant coefficients for the pure species i. Domain-based weighted-sum of gray gas model (WSGGM) is used for the calculation of the emissivity of the flue gas mixture. The basic assumption of the WSGGM is that the total emissivity over a path length scan be calculated as the weighted sum of a number of gray gasses and a transparent gas: I

ε = ∑ aε ,i (T )(1 − e −κ ps ) i

(22)

i =0

The weighting factor of the transparent gas can be calculated: I

aε , 0 = 1 − ∑ aε , i

(23)

i =1

The temperature dependence of aε ,i can be approximated by any function, but the most common approximation is a temperature polynomial: J

aε ,i = ∑ bε ,i , j T

j −1

(24)

j =1

The absorption coefficient of the flue gas can be obtained from the following two formulae: when > 10 :



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 49

ln(1 − ε ) s

(25)

α = ∑ aε , i κ i p

(26)

α =− when 10 : I

i =0

3. SIMULATION SET-UP 3.1. Furnace Geometry and Operating Conditions. Figure 1 shows one-eighth of the considered naphtha cracking furnace. The dimensions of the complete furnace and the operating conditions are summarized in Table 1. Fuel gas is uniformly distributed over 64 floor burners and per burner divided over two types of inlets. The primary inlets lie inside the air inlet, while the secondary burners are outside of the air inlet as shown in Figure 1c. All reactors comprise two passes: an inlet and an outlet pass. The inlet passes enter the furnace from the top and bend alternating towards both walls through an S-type bend close to the furnace floor. U-type bends connect each inlet pass to its outlet pass realigning all reactors in a single row.

Figure 1.Schematic of the naphtha cracking furnace: (a) overall geometry of the furnace; (b) detail of the reactor tube bends; (c) details of a burner.


Page 13 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1.Furnace dimensions and operating conditions

Furnace specifications Length [m] (x-direction) 24.89 Width [m] (y-direction) 2.964 Height [m] (z-direction) 11.609 Number of reactor tubes [-] 176 Number of bottom burners [-] 64 Firing conditions Fuel gas ﬂow rate in bottom [g/s] 2221.9 Excess air [vol%] 1.1 Fuel composition [mol%] CH4 88.55 H2 11.14 CO 0.17 C 2 H4 0.14 Material properties Emissivity of the furnace wall [-] 0.75 Emissivity of the reactoralloy [-] 0.85 Reactor coils Inlet tube inner diameter [m] 0.045 Outlet tube inner diameter [m] 0.051 Inlet tube outer diameter [m] 0.0556 Outlet tube outer diameter [m] 0.0666 Feed rate [g/s] 16088.9 Steam dilution ratio [kg/kg] 0.5 Coil inlet temperature [K] 853.15 Coil outlet pressure [MPa] 0.077 Feedstock composition [wt.%] Paraffins 58.5 Olefins 0.0 Naphthenics 31.5 Aromatics 10.0 3.2. Boundary Conditions. In the furnace simulation, the inlet boundary conditions of the fuel gas and air from the burners are set as mass flow inlet boundaries. The outlet boundary condition of the furnace is set as a pressure outlet boundary with a value of-24.5 Pa to account for the under-pressure created by the fan in the furnace stack. The no-slip boundary condition is imposed to all walls. The energy boundary condition on the furnace refractory walls is set as a uniform heat flux based on the heat loss value calculated from a heat balance over the furnace using the plant design data. Because of structural similarity, only a representative segment of one-eighth of the industrial cracking furnace is simulated. Symmetrical boundary conditions are 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 49

applied on the side faces of the segment basically neglecting the influence of the furnace side walls on the fluid dynamics and heat transfer. The tube external wall temperature of the reactor tubes is obtained from the reactor simulations and is assigned by a user-defined function (UDF). The inlet turbulence properties are specified using the turbulence intensity and the turbulence length scale. The turbulence length scale is defined as 10% of the inlet hydraulic diameter and the turbulence intensity is calculated from the inlet Reynolds number using the following correlation:26 I = 0.16 Re −1 / 8

(27)

In the reactor simulations, the mass flow rate, composition and temperature at the inlet are set. At the outlet the coil-outlet-pressure (COP) is set and at the reactor tube external wall the heat flux profile obtained from the furnace simulation is imposed. 3.3. Numerical Solution. The governing equations of the furnace are solved numerically using the commercial CFD package FLUENT/14.0 which uses the finite volume method. The SemiImplicit Method for Pressure-Linked Equations (SIMPLE) is used for pressure-velocity coupling. As lightly more conservative under-relaxation factor of 0.7 is used for pressure, momentum, k and ε, and of 0.8 for temperature and radiative intensity. The convergence was judged by monitoring residuals, the overall mass and energy balances and the flue gas outlet temperature. The convergence criterion is that all scaled residuals are below 10−3. The exceptions are the energy and the radiation equations, where the criterion is 10−6. The net mass and energy imbalance should be less than 0.2% of the flux through the domain. The flue gas outlet temperature was monitored and changed less than 0.05 K over the final 100 iterations. Tetrahedral cells were used to discretize the computational domain around the burners and the reactor tubes. Hexahedral cells were used to discretize the rest of the computational domain. The total number of grid cells was 4,392,849. The adopted mesh is the same as used by Hu et al.14, 15 where mesh independence was verified.


Page 15 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.4. Coupled Furnace-Reactor Simulation. As heat to the reactors is supplied by the furnace burners and the temperature difference between the process gas and the reactor outer wall is the driving force for heat transfer, the fireside and the process gas side of a steam cracking furnace are strongly coupled. A coupled reactor-furnace simulation is necessary to obtain a closed heat balance over the entire unit as a change in heat flux to the reactors has an influence on the cracking reactions in the coil and the process gas temperature profile. As such, only coupled simulation of furnace and reactor can accurately predict the effect of the radiation model on key process parameters such as furnace efficiency, reactor external wall temperature, heat flux profiles to the reactors and run length. Furthermore, the iterative simulation procedure for coupled simulations of Hu et al.14, 15 is applied. In the first iteration, the industrially measured reactor external wall temperature profile is adopted in the furnace simulation. The simulated heat flux profile to the reactor tubes is then applied in the reactor simulation, yielding updated reactor external wall temperature profiles for the furnace simulation. This procedure is repeated until the maximum difference of the tube external wall temperature between two successive iterations is below a specified predefined threshold value of 1 K. 4. RESULTS AND DISCUSSION 4.1. Furnace Simulation. Figure 2 compares the flue gas z-velocity along the width of the furnace at x=5.377 m(see Figure 1) and at different heights for the four radiation models. At a height of 1 m, the jet zones near the furnace bottom at the burner locations are clearly seen. These velocity peaks decrease with increasing z-coordinate. As can be seen in Figure 2a and 2d, at a height of 1 m in the furnace the velocity varies near the burners from 40 to 50 m/s, while at a height of 10 m, i.e. near the convection zone the velocity range is 4-6 m/s and 1-3 m/s respectively, corresponding to the furnace outlet side and side far away from the furnace outlet. The reason of different profile is due to the asymmetrical structure of the furnace. More importantly for this study it can be seen 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 49

that the impact of the radiation model on the flue gas velocity profile is minor. Only for the adiabatic case clear differences can be observed with the other radiation models. It is clearly observed from Figure 2a-d that along the y direction of the furnace zero and negative z-velocities can be observed near the reactor tubes. The former is a result of imposing the non-slip boundary condition at the reactor tubes, while the latter is caused by the recirculation of the flue gas on both sides of the reactor tubes. The negative z-velocity close to the reactor tubes increases with increasing z-coordinate and confirms the flue gas recirculation typically encountered in floor-fired steam cracking furnaces.14,

15

This recirculation promotes

mixing of the flue gas, resulting in a more uniform temperature profile and a more uniform heat flux profile to the reactors.27 The small discontinuities seen in all graphs of Figure 2 are the locations of the reactor tubes where the flue gas z-velocity is obviously nonexistent. Owing to the bends in the reactors, the position of the discontinuities changes along the furnace height. For all heights the Adiabatic model overpredicts the velocity compared to the other models. DOM, P-1 and DTRM give similar results. The overprediction in the adiabatic simulation can be attributed to the higher flue gas temperature resulting in a lower flue gas density and higher z-velocities as discussed further.


Page 17 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2.Flue gas z-velocity component along the width of the furnace in the plane x=5.377 m at different heights;(a)1 m; (b)2 m; (c)4 m; (d)10 m for different radiation models: DOM; and Adiabatic;

,DTRM;

,P-1;

.

Figure 3 shows the flue gas temperature in a vertical cross section along the furnace length at y=2.75 m. The flue gas temperature of the adiabatic simulation is unrealistically high due to the absence of radiative heat transfer. For DOM, DTRM and P-1 model clear flame shapes are simulated. The results of DOM and DTRM simulations are very similar. However, slightly lower flue gas temperatures are simulated with the P-1 model, consistent with the results of Habibi et al.22, 23



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 49

Figure 3.Flue gas temperature [K] in a vertical cross section at y=2.75 m for different radiation models; from left to right:DOM,DTRM,P-1 and Adiabatic.

Figure 4 shows the flue gas temperature mixing-cup averaged over a horizontal plane along the furnace height for the different radiation models. The average flue gas temperature of the Adiabatic model is the highest, in agreement with Figure 3. The predictions of DOM and DTRM model are close with a maximum difference of 15 K. The P-1 model gives flue gas temperatures that are on average about 30 K lower than DOM. Maximum values are predicted at a height of about 5 m by DOM, DTRM and P-1 models, while the maximal flue gas temperature occurs at a height of 9 m for the Adiabatic model. In the outlet of furnace, the flue gas temperatures of DOM, DTRM, P-1 and Adiabatic model are respectively 1357K, 1361K, 1322K and 1993K, while the design data is 1364K (this is seen in Table 2 as mentioned below). Compared the different radiation models with the design data, the results of DOM and DTRM model are close to the design data, while those of P-1 and Adiabatic model are respectively lower and higher than the design data.


Page 19 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4.Mixing-cup averaged flue gas temperature [K] along the height of the furnace for different radiation models: DOM; ,DTRM; ,P-1; and Adiabatic; .

Figure 5 shows the incident radiation in the plane x=2.3 m for DOM and P-1 radiation models. Note that different ranges are used in Figure 5a and Figure 5b as the differences are too large to use a single range. It represents the rate of incident radiation per unit area originating from all directions over all wavelengths. It can be seen that the incident radiation for DOM is much higher than for the P-1 model. This corresponds with the flue gas temperature being higher with DOM compared to P-1, which is in agreement with Habibi et al.22



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 49

Figure 5.Incident radiation [W/m2] in a vertical cross-section at x = 2.3 m fordifferent radiation models: (a)DOM; (b) P-1.

Figure 6 shows the flue gas absorption coefficient in the plane x=2.3 m for the different radiation models. The maximum value is around 0.2 which is a typical value for combustion processes.43 The highest values are calculated where carbon dioxide and steam concentrations are high as these are the main absorbing species at the prevailing temperatures. Comparing the results for the three different radiation models, highly similar results are simulated.


Page 21 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 6.Absorption coefficient in a vertical cross-section at x =2.3 m for different radiation models (a)DOM;(b)DTRM; (C) P-1.

Figure 7 shows the computed contours of the optical thickness in the plane x=2.3 m for the different radiation models. In light of the definition of Habibi et al.,22 the optical thickness is defined as the product of the local mean absorption coefficient and the local integral turbulent length scale (k3/2/ε). For all the models, the optical thickness in the jet core zones of the burners is rather low due to the locally low absorption coefficient. Outside the jet core, larger values of the optical thickness are simulated due to the abundance of H2O and CO2. At the upper left side close to the reactors, low velocities and accompanying low kinetic energy levels explain the low optical thickness. A similar reasoning holds for the two recirculation zones in the outlet. As the optical thickness is below 3 throughout the entire domain, the Rosseland model is not applicable for steam cracking furnaces as shown by Habibi et al.22 The P-1 model is applicable at optical thicknesses above 1.22 Hence, it is clear that the P-1 model is used outside its applicability range 21 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 49

in some regions in the furnace. Even though the CPU and memory requirements are higher compared to the use of P-1, only DOM and DTRM should be used for the simulation of steam cracking furnaces.

Figure 7.Optical thickness in a vertical cross-section at x =2.3 m for different radiation models (a)DOM; (b)DTRM; (C) P-1.

Figure 8 shows the CO mass fraction in a vertical cross-section along the furnace length at y = 0.214 m, i.e. the center of the front wall burners. In the jet core zones of the burners, incomplete combustion of the fuel gas leads to high CO values. The flame height is calculated based on the CO mass fraction. In this case it is defined as the distance following the center line of the burner between the burner tip and the height at which the CO mass fraction is zero.14 Outside the jet core zones, complete combustion conditions lead to conversion of CO to CO2. It also shows that the flame height of DOM, DTRM and P-1 models is shorter than that of the Adiabatic model. As the temperature in the adiabatic simulation is higher, the density is lower resulting in higher velocities, leading to longer flames. 22 ACS Paragon Plus Environment

Page 23 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 8.CO mass fraction [wt%] in a vertical cross-section at y=0.214; from left to right: DOM, DTRM, P-1 and Adiabatic.

4.2. Reactor Simulation. Figure 9 shows the heat flux to the reactors per internal surface area, the external tube metal temperature and the process gas temperature along the reactor axial position respectively at start of run (SOR) conditions. In the adiabatic simulation, a very low heat flux to the reactors is simulated because of the neglect of the radiative heat transfer, resulting in unrealistically low external wall and process gas temperatures. Hence, the reactor results of the adiabatic simulation are not shown. As seen from Figure 9a, the heat flux profiles for DOM and DTRM are very similar and both show two maxima, i.e. one in each pass. With the P-1 model a higher heat input to the first pass of the reactors is simulated compared to DOM and DTRM, while a lower heat input to the second pass of the reactors is simulated. Habibi et al.23 simulated a higher heat input with the P-1 model compared to DOM along the entire reactor. Habibi et al.22, 23

did not evaluate the performance of the DTRM. The higher heat flux to the first pass with the

P-1 model leads to a higher process gas temperature at the inlet of the second pass of the reactor compared to DOM and DTRM. This leads to a lower driving force for heat transfer from the furnace to the second pass, resulting in the lower heat flux to the second pass with the P-1 model compared to DOM and DTRM. Habibi et al.22,

23

did not perform a coupled furnace-reactor

simulation but imposed a predefined process gas temperature. As a result the interaction between 23 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 49

the process gas temperature on one hand and the heat flux on the other hand is not accounted for. Hence, the heat fluxes to the reactor tubes predicted by Habibi et al. using the P-1 model are too high. This consideration shows the importance of a coupled furnace-reactor simulation for the evaluation of different radiation models as performed in this work. The reactor external wall temperature profiles depicted in Figure 9b basically show the same trends as the heat flux profile. DOM and DTRM show very similar results because of the similar heat flux profiles. With P-1 higher tube external wall temperatures are simulated in the first pass, whereas the profile flattens out in the second pass because of the lower heat input. In Figure 9c, the process gas temperature profile for the P-1 model increases faster than for DOM and DTRM because of the higher heat flux to the first pass. However, after 15 m the increase in temperature is highly similar because of the similar heat input. The similar increase in temperature after 15 m for P-1 compared to DOM and DTRM notwithstanding the lower heat input is explained by the higher cracking severity with P-1. As the cracking severity is higher, the heat of reaction is lower because of the higher reaction rates of exothermic reactions forming aromatics.

Figure 9. Heat flux per tube internal surface area[kW/m² internal] (a)SOR external wall temperature [K] (b) and Process gas temperature [K] (c) as a function of axial position with different radiation models: DOM; ,DTRM; and P-1; .

Figure 10 and Figure 11 show the ethene and propene yield along the reactor axial position respectively. For all radiation models, the ethene yield increases monotonically with the reactor axial position and the propene yield has a maximum. Because of the higher process gas temperature for P-1, cracking starts faster resulting in a higher cracking severity. This results in a


Page 25 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


higher ethene yield and a shift of the propene maximum upstream of the reactor because of the secondary reactions converting propene to ethene more upstream the reactor.

Figure 10.Ethene yield [wt%] profile as a function of axial position with different radiation models: DOM; , DTRM; and P-1; .

Figure 11.Propene yield [wt%] profile as a function of axial position with different radiation models: DOM; , DTRM; and P-1; .

In order to quantify the effect of the different radiation models on the predicted run length, three run length simulations were performed. Run length simulations require that reduction of the cross-sectional flow area due to a growing coke layer is accounted for. COILSIM1D combines the model of Plehiers et al.8 for light feedstocks and the model of Reyniers et al.44 for 25 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 49

heavy feedstocks. As coke formation is much slower than the gas-phase cracking reactions, pseudo-steady state can be assumed for coke formation, i.e. the run time is increased in a stepwise manner. In each time step, the coking rate is calculated using the prevalent temperature and concentrations. By multiplication of this coking rate with the time step and the coke density, the coke layer growth in the current time step is determined. Coupled furnace-reactor run length simulations require a coupled furnace-reactor simulation at each time step and an iteration loop within each time step to adjust the fuel flow rate to maintain the desired cracking severity.8, 45 This is computationally not possible for the detailed furnace simulations of this work but requires simplification to a coarser grid and the assumption of a predefined combustion profile within the furnace.8, 45 In this work, the shape of the heat flux profile to the reactors is assumed to be constant over time and the heat flux profile is scaled in each time step to maintain the startof-run propene over ethene ratio. Three simulations were performed with the heat flux profiles shown in Figure 9a, i.e. for DOM, DTRM and P-1. No simulation was performed for the adiabatic case as the low heat input yields unrealistic results. The desired propene over ethene ratio for the three simulations was set to the SOR P/E ratio, i.e. 0.52, 0.52 and 0.42 for DOM, DTRM and P-1 respectively. The decoke criterion for this furnace is the external tube metal temperature reaching 1398.15 K. Figure 12a and b show the maximum tube metal temperature and the inlet pressure as a function of run time for the three simulations. The decoke criterion is calculated to be reached after 53, 51 and 56 days of operation when using DOM, DTRM and P-1 respectively. These are all slightly lower than the design run length of 60 days. Despite the higher cracking severity in the P-1 simulation compared to DOM and DTRM simulations, i.e. 0.42 compared to 0.52, a longer run length is simulated in the P-1 case. This shows the importance of a correct shape of the heat flux profile for an accurate run length prediction, justifying the detailed and computationally intensive simulations performed, analyzed and presented in this work. Figure 13 shows the coking rate as a function of reactor axial position at


Page 27 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


start of run conditions. The maximum coking rate is similar in all three cases and occurs at the axial position with maximum tube metal temperature (TMT), i.e. around 18 m. More upstream the reactor, the coking rate is higher when using the P-1 model but this does not affect the increase in the maximum TMT as seen from Figure 12a because the maximum TMT is determined by the maximum coking rate and not by the average coking rate. On the contrary, the reactor pressure drop increases faster with P-1 compared to DOM and DTRM as seen from Figure 12b because of the higher average coking rate resulting in a lower cross-sectional flow area along the entire reactor.

Figure 12.Maximum tube metal temperature [K] (a) and reactor pressure drop [kPa] (b) as a function of run time [days]: DOM; , DTRM; and P-1; .



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 49

Figure 13.SOR coking rate [10-6 kg/m²/s] as a function of axial position with different radiation models: DOM; , DTRM; and P-1; .

Table 2 compares the most important simulation results obtained using the different radiation models with the plant design data. Due to the higher total heat input, the coil-outlet-temperature (COT) is higher and propene to ethene ratio is lower with the P-1 model compared to DOM and DTRM. However, the maximum tube external wall temperature is similar as the heat flux to the second pass, i.e. where the maximum external tube wall temperature is reached, is lower with the P-1 model. The radiation efficiency reported in Table 2 is defined as the ratio of the radiative heat transfer to the reactors to the total heat transfer to the coils. As can be seen from Table 2 the dominant heat transfer mechanism is radiation in line with previous work.4 The radiation efficiency when applying the DOM, DTRM and P-1 models is higher than 90%.The higher total heat transfer to the reactors with the P-1 model results in a higher thermal efficiency of the furnace. The thermal efficiency is defined as the ratio of the heat transfer to the reactors to the total heat input to the furnace, i.e. both enthalpy inflow from the air and fuel inlets and the heat of combustion. The SOR thermal efficiency of the furnace with DOM and the DTRM is within 0.6% of the thermal efficiency of the plant design data. However, DOM and DTRM simulated coil-outlet-temperature and the maximum start-of-run external wall temperature are about 8-9 28 ACS Paragon Plus Environment

Page 29 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and 10-15 K lower compared to the design data respectively. This results in a lower simulated cracking severity, i.e. a higher propene to ethene ratio. Even though a higher propene to ethene ratio is simulated, the simulated ethene yield for DOM and DTRM is higher than the design ethene yield which might cast doubt on the kinetic model that was adopted to determine the design data. Comparing the simulated results of all radiation models to the design data, DOM and DTRM perform best. However, DOM has a broader application range than the DTRM model.26 Furthermore, contrary to DOM, the current implementation of the DTRM in FLUENT/14.0 is not compatible with parallel processing. As a result, the simulation clock time increases by a factor of eleven compared to DOM simulation, i.e. about two weeks for a DTRM simulation. Therefore DOM is the recommended radiation model to be used for the simulation of steam cracking furnaces. Table 2.Comparison of simulation results with different radiation models and design data

Items Flue gas outlet temperature [K] Coil-outlet-temperature [K] Propene to ethene ratio [wt.%/wt.%] Maximum SOR tube external wall temperature [K] Radiation efficiency [%] Thermal efficiency of the furnace [%] Simulation clock time of one furnace simulation [ks] Run length [days] Yields [wt%] Ethene Propene

Design data 1364 1161.2 0.44

DOM

DTRM

P-1

Adiabatic

1357 1153.2 0.52

1361 1152.2 0.52

1322 1178.2 0.42

1993 959.2 0.89

1268.2

1253.2

1258.2

1255.2

1004.2

N.A. 45.45

95.67 46.00

95.54 45.73

96.88 48.00

0 10.16

N.A.

110

1210

90

72

60

53

51

56

N.A.

28.03 12.34

29.63 15.41

29.59 15.49

30.66 13.00

2.83 2.51

5. CONCLUSIONS Coupled furnace/reactor simulations of an industrial naphtha cracking furnace with 130kt/a capacity were performed. In the furnace model, the compressible formulation of the ReynoldsAveraged Navier-Stokes (RANS) equations was adopted to simulate the flue gas flow. For the reactor simulations a plug flow model combined with a single-event microkinetic free-radical



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 49

model was adopted. For the first time different radiation models were evaluated in coupled furnace/reactor simulations to predict run lengths. The four evaluated radiation models are the Adiabatic, P-1, discrete ordinates and discrete transfer radiation model. Adopting the Adiabatic model resulted in unrealistically high flue gas temperatures, due to the neglect of radiative heat transfer which is the main mode of heat transfer in the furnace. The results with DOM and the DTRM radiation model are very similar. The flue gas temperature with DOM is higher than with the P-1 radiation model resulting in higher incident radiation. Calculation of the optical thickness shows that the optical thickness is below 1in some regions of the furnace. The P-1 model is only valid for an optical thickness above 1, explaining the discrepancies with the DTRM and DOM. With the P-1 model a higher total heat input to the reactors is simulated compared to DOM and DTRM. Also the heat flux profile along the reactor axial length is different with P-1 compared to DOM and DTRM. Comparing the simulated results of all radiation models to the plant design data, DOM and DTRM perform best. As DOM has a broader application range than the DTRM, DOM is the recommended radiation model to be used for the simulation of steam cracking furnaces, in particular when run length simulations are carried out. AUTHOR INFORMATION Corresponding author: Feng Qian, No. 130 of Meilong Road, Shanghai 200237, Email: [email protected] Kevin M. Van Geem, Laboratory for Chemical Technology, Ghent University, Technologiepark 914, B-9052 Ghent, Belgium, Email:[email protected] ACKNOWLEDGMENT This work is supported by Major State Basic Research Development Program of China (2012CB720500), National Natural Science Foundation of China (Key Program: U1162202), National Science Fund for Outstanding Young Scholars (61222303), National Natural Science


Page 31 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Foundation of China (21276078, 61174118), Shanghai Key Technologies R&D Program (12dz1125100) and Shanghai Leading Academic Discipline Project (B504). The financial support from the BOF Bilateral Scientific Cooperation (ECUST/LCT) and the Long Term Structural Methusalem Funding by the Flemish Government (No.BOF09/01M00409) and ‘111_ Project by the Chinese Government (No. B08021) are acknowledged. CMS acknowledges financial support from a doctoral fellowship from the Fund for Scientific Research Flanders (FWO). NOMENCLATURE

Roman

bε ,i , j

=

emissivity gas temperature polynomial coefficients,

-

C=

linear-anisotropic phase function coefficient,

-

c1,c2=

constants for each component,

-

Cp=

heat capacity of mixture,

J/kg/K

f=

mixture fraction,

-

mean mixture fraction,

-

mixture fraction variance,

-

Fi=

molar flow rate of component i,

kmol/s

G=

incident radiation,

W/m2

I=

radiation intensity,

J/m2/s

k=

turbulent kinetic energy, absorption coefficients,

m2/s2, -

kn =

specific absorption coefficients of the nth gray gas,

Pa-1m-1

L=

characteristic length of the computational domain,

m

LM =

mean beam length of flue gas,

m

Mi =

molecular weight of species i,

g /mol

r n=

normal pointing out of Domain,

-

f= f ′2 =



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 49

p=

static pressure, sum of the partial pressures of all absorbing gases,

Pa

qr=

radiation flux,

W/m2

qin=

intensity of the incoming ray,

W/m2

qout=

net radiative heat flux from the surface,

W/m2

r r =

position vector,

-

R=

ideal gas constant,

R=8.314 J/mol/K

S=

furnace wall area,

m2

r s=

direction vector,

-

r s′ =

scattering direction vector,

-

s

path length,

m

Sh =

source term in the energy equation,

J/m3/s

SG=

user-defined radiation source,

W/m2

Sφ=

source term,

-

T =

local temperature and process-gas temperature,

K

Tw=

surface temperature of the point P on the furface,

K

=

Ui, Uj,Ul= velocity component in the ith, jth or lth direction,

m/s

V=

total gas volume,

m3

xi,xj, xl=

coordinate direction in the ith, jth or lth direction,

m

Yj =

mass fraction of species j,

-

absorption coefficient,

1/m

α ε ,i =

emissivity weighting factors for the ith fictitious gray gas,

-

α g ,n =

emissivity weighting factors for the nth fictitious gray gas,

-

ρ =

gas density,

kg/m3

µ =

viscosity of gas molecules,

kg/m/s

Greek

α

=


Page 33 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


µt =

turbulent viscosity,

kg/m/s

ε=

dissipation rate of turbulent kinetic energy, emissivity,

m2/s3,-

εg =

emissivity of the real gas,

-

εw

wall emissivity,

-

∆ε =

correction factor of emissivity,

-

τ =

optical thickness,

-

Γφ =

generalized diffusion coefficient,

-

σ=

Stefan-Boltzmann constant,

σ =5.672×10-8W/m2 K4

σS =

scattering coefficient,

1/m

λ=

thermal conductivity,

W/m/K

φ=

dependent variable,

-

Φ=

phase function,

-

Ω=

hemispherical solid angle,

degrees

Ω′ =

solid angle,

degrees

i=

symbol for spatial coordinates, symbol for pure species,

-

j=

symbol for spatial coordinates, symbol for pure species,

-

l=

symbol for spatial coordinates,

-

in=

at tube inlet

-

out=

at tube outlet

-

ox=

value at the oxidizer stream inlet,

-

fuel=

value at the fuel stream inlet,

-

=

Subscript

REFERENCES (1) Muñoz Gandarillas, A. E.; Van Geem, K. M.; Reyniers, M. F.; Marin, G. B. Coking Resistance of Specialized Coil Materials during Steam Cracking of Sulfur-Free Naphtha. Ind. Eng. Chem. Res. 2014, 53 (35), 13644-13655.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 49

(2) Muñoz Gandarillas, A. E.; Van Geem., K. M.; Reyniers, M. F.; Marin, G. B. Influence of the reactor material composition on coke formation during ethane steam cracking. Ind. Eng. Chem. Res. 2014, 53 (15), 6358-6371. (3) Heinz Zimmermann; Roland Walzl. Ethylene. in ULLMANN'S Encyclopedia of Industrial Chemistry; John Wiley & Sons, Inc., 2012; Vol.13, pp 465-529. Online: DOI: 10.1002/14356007.a10_045.pub3 (4) Ramana Rao, M. V.; Plehiers, P. M.; Froment, G. F. The coupled simulation of heat transfer and reaction in a pyrolysis furnace. Chem. Eng. Sci. 1988, 43 (6), 1223-1229. (5) Xiang, J. F. Progress in Research of Heat Transfer Model for Cracking Tubes in Radiation Section of Furnace. Shiyou Huagong 2007, 36 (4), 407-411. (6) Hou, K. H.; Wang, Z. X. Modeling calculation for radiant chamber of a furnace for light hydrocarbon pyrolysis with steam by using Lobo-Evance's sectional method. Daqing Shiyou Xueyuan Xuebao 1989, 13 (4), 55-63. (7) Zhu, H. D.; Zhang, L. Simulation of dichloroethane pyrolysis process. Huaxue Fanying Gongcheng Yu Gongyi 1995, 11 (3), 251-258. (8) Plehiers, M. P.; Reyniers, G. C.; Froment, G. F. Simulation of the run length of an ethane cracking furnace. Ind. Eng. Chem. Res. 1990, 29 (4), 636-641. (9) Qian, J. L. Tubular heating furnace; China Petrochemical Press: Beijing, 2003. (10) Wang, G. Q.; Zhang, L. J. Numerical simulation of gas flow and heat transfer in cracking heater. Petrochem. Technol. 2005, 34 (7), 652-655. (11) Zhou, H. Z.; Jia, Z. G. Three dimensional numerical simulation of flow and combustion in firebox of ethylene cracking furnace. Petrochem. Technol. 2007, 36 (6), 584-590. (12) Coelho, P. J. Numerical simulation of radiative heat transferfrom non-gray gases in threedimensional enclosures. J. Quant. Spectrosc. Radiat. Transf. 2002, 74 (3), 307-328. (13) Stefanidis, G. D.; Merci, B.; Heynderickx, G. J.; Marin, G. B. Gray/nongray gas radiation modeling in steam cracker CFD calculations. AIChE J. 2007, 53 (7), 1658-1669. (14) Hu, G. H.; Wang, H. G.; Qian, F.; Van Geem, K. M.; Schietekat, C. M.; Marin, G. B. Coupled simulation of anindustrial naphtha cracking furnace equipped with long-flame and radiation burners. Compu. Chem. Eng. 2012, 38, 24-34. (15) Hu, G. H.; Wang, H. G.; Qian, F.; Zhang, Y.; Li, J. L.; 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 (24), 13672-13685. (16) Yang, J. W.; Tai, N. N.; Wang, L. J.; Xiao, J. Z.; Yang, C. H. Numerical simulation of the flue gas and process side of coking furnaces. Ind. Eng. Chem. Res. 2012, 51, 15440-15447. (17) Hassan, G.; Pourkashanian, M.; Ingham, D.; Ma, L.; Newman, P.; Odedra, A. Predictions of CO and NOx emissions from steam cracking furnaces using GR12.11 detailed reaction mechanism-A CFD investigation. Comput. Chem. Eng. 2013, 58, 68-83. 34 ACS Paragon Plus Environment

Page 35 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(18) Keramida, E. P.; Liakos, H. H.;Founti, M. A.;Boudouvis, A. G.; Markatos, N. C. Radiative heat transfer in natural gas-fired furnaces. Int J Heat Mass Transf 2000, 43 (10), 1801-1809. (19) Li, T.; Li, W. L.; Yuan, Z. L. Different radiative models for heat and mass transfer characteristics in vertical pipe. Zhongguo Dianji Gongcheng Xuebao 2007, 27 (2), 92-98. (20) Mendes, M. A. A.; Pereira, J. M. C.; Pereira, J. C. F. A numerical study of the stability of one-dimensional laminar premixed flames in inert porous media. Combust. Flame 2008, 153 (4), 525-539. (21) Zhang, Y.; Schietekat, C. M.; Qian, F.; Van Geem, K. M.; Marin,G. B. Impact of flue gas radiative properties and burner geometry in furnace simulations. AIChE J. 2014, 10.1002/aic.14724. (22) Habibi, A.; Merci, B.; Heynderickx, G. J. Impact of radiation models in CFD simulations of steam cracking furnaces. Comput. Chem. Eng. 2007, 31, 1389-1406. (23) Habibi, A.; Merci, B.; Heynderickx, G. J. Multiscale modeling of turbulent combustion and NOx emission in steam crackers. AIChE J. 2007, 53 (9), 2384-2398. (24) Lan, X. Y.; Gao, J. S.; Xu, C. M.; Zhang, H. M. Numerical simulation of transfer and reaction processes in ethylene furnaces. Trans IChemE, Part A, Chem. Eng. Res. Des. 2007, 85 (A12), 1565-1579. (25) Tao, W. Q. Numerical Heat Transfer, Second Edition; Xi'an Jiaotong University Press: Xi'an, 2001. (26) ANSYS FLUENT Theory Guide; Release 14.0: ANSYS, Inc, 2011. (27) Stefanidis, G. D.; Merci, B.; Heynderickx, G. J.; Marin, G. B. CFD simulations of steam cracking furnaces using detailed combustion mechanisms. Comp. Chem. Eng. 2006, 30 (4), 635649. (28) Lysenko, D.; Ertesvåg, I.; Rian, K. Numerical Simulations of the Sandia Flame D Using the Eddy Dissipation Concept. Flow Turbul Combust 2014, 93 (4), 665-687. (29) Modest, M. F. Radiative Heat Transfer, Second Edition. Academic Press: New York, 2003. (30) Siegel, R.; Howell, J. R. Thermal Radiation Heat Transfer; Hemisphere Publishing Corporation: Washington,D.C., 1992. (31) Dente, M.; Ranzi, E.; Goossens, A. G. Detailed prediction of olefin yields from hydrocarbon pyrolysis through a fundamental simulation model (SPYRO). Comput. Chem. Eng. 1979, 3, 61-75. (32) Van Geem, K. M.; Heynderickx, G. J.; Marin, G. B. Effect of radial temperature profiles on yields in steam cracking. AIChE J. 2004, 50(1), 173-183. (33) Sundaram, K. M.; Froment, G. F. Two dimensional model for the simulation of tubular reactors for thermal cracking. Chem. Eng. Sci. 1980, 35, 364-371. 35 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 36 of 49

(34) Schietekat, C. M.; Van Cauwenberge, David J.; Van Geem, K. M.; Marin, G. B. Computational Fluid Dynamics-Based Design of Finned Steam Cracking Reactors. AIChE J. 2014, 60 (2), 794-808. (35) 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. (36) Detemmerman, T.; Froment, F. Three dimensional coupled simulation of furnaces and reactor tubes for the thermal cracking of hydrocarbons. Rev Inst Fr Pet 1998, 53, 181-194. (37) De Saegher, J. J.; Determmerman, T.; Froment, G. F. Three dimensional simulation of high severity internally finned cracking coils for olefins production. Rev Inst Fr Pet 1996, 51 (2), 245-260. (38) Van Geem, K. M.; Hudebine, D.; Reyniers, M. F.; Wahl, F.; Verstraete, J. J.; Marin, G. B. Molecular reconstruction of naphtha steam cracking feedstocks based on commercial indices. Comput. Chem. Eng. 2007, 31 (9), 1020-1034. (39) Van Geem, K. M.; Marin, G. B.; Hedebouin, N.; Grootjans, J. Energy efficiency of the cold train of an ethylene cracker. Oil Gas Eur Mag 2008, 34, 95-99. (40) 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. (41) Van Geem, K. M.; Pyl, S. P.;Reyniers, M. F.; Vercammen, J.;Beens,J.;Marin, G. B. Online analysis of complex hydrocarbon mixtures using comprehensive two-dimensional gas chromatography. J. Chromatogr. A 2010, 1217 (43), 6623-6633. (42) Pyl, S. P.; Van Geem, K. M.; Reyniers, M. F.; Marin, G. B. Molecular reconstruction of complex hydrocarbon mixtures: an application of principal component analysis. AIChE J. 2010, 56 (12), 3174-3188. (43) Viskanta, R.; Menguc, M. P. Radiation heat transfer in combustion systems. Prog. Energy Combust. Sci. 1987, 13 (2), 97-160. (44) Reyniers, G. C.; Froment, G. F.; Kopinke, F. D.; Zimmermann, G. Coke formation in the thermal cracking of hydrocarbons. 4. Modeling of coke formation in naphtha cracking. Ind. Eng. Chem. Res. 1994, 33 (11), 2584-2590. (45) 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.


Page 37 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


200x148mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

560x438mm (96 x 96 DPI)


Page 38 of 49

Page 39 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


374x168mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

107x80mm (96 x 96 DPI)


Page 40 of 49

Page 41 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


128x111mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

131x131mm (96 x 96 DPI)


Page 42 of 49

Page 43 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


133x121mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

365x162mm (96 x 96 DPI)


Page 44 of 49

Page 45 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


370x105mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

106x79mm (96 x 96 DPI)


Page 46 of 49

Page 47 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


112x84mm (96 x 96 DPI)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

165x75mm (96 x 96 DPI)


Page 48 of 49

Page 49 of 49

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


115x86mm (96 x 96 DPI)


Impact of Radiation Models in Coupled Simulations of Steam

Recommend Documents