Comprehensive Simulation and Optimization of an Ethylene

Article pubs.acs.org/IECR

Comprehensive Simulation and Optimization of an Ethylene Dichloride Cracker Based on the One-Dimensional Lobo−Evans Method and Computational Fluid Dynamics Chaochun Li, Guihua Hu, Weimin Zhong,* Hui Cheng, Wenli Du, and Feng Qian* Key Laboratory of Advanced Control and Optimization for Chemical Processes, Ministry of Education, East China University of Science and Technology, Shanghai 200237, China ABSTRACT: Coupled simulations of an ethylene dichloride (EDC) cracking furnace and reactor are conducted using onedimensional Lobo−Evans and computational fluid dynamics (CFD) models. Optimization is performed using the first model, in which the fuel gas allocation operator α is investigated to improve performance indices such as selectivity, conversion, and fuel gas consumption (per vinyl chloride monomer production). The optimum coil outlet temperature (COT) is suggested to make a good compromise among the performance indices. The CFD model is used to validate the optimized results. A standard k−ε two-equation model is applied to simulate turbulence, and a finite-rate/eddy dissipation model is used to model a premixed combustion of the sidewall burners. The discrete ordinate model is applied to simulate the radiative heat transfer of a furnace in a CFD simulation. The EDC cracking process in the reactor, as well as the flow, combustion, and radiative heat transfer in the furnace, is provided in the CFD model.

1. INTRODUCTION A vinyl chloride monomer (VCM) plant consists of oxychlorination, ethylene dichloride (EDC) rectification, EDC thermal cracking, VCM purification, and hydrogen chloride recovery. Thermal cracking is the core process and major energy consumer by converting EDC feedstock into VCM. An EDC cracker is equipped with furnace and reactor tubes. Fuel gas and air enter the furnace through burners, going through diffusion combustion and turbulent mixing. Heat profiles generated in the furnace are transferred to the process gas inside the reactor tubes, raising the process gas temperature and supporting endothermic cracking reactions along the reactor tubes. The furnace is mainly divided into convection and radiation sections. The EDC raw material in the reactor tubes is preheated in the convection section, whereas the EDC cracking reactions in the reactor tubes occur in the radiation section when the temperature of the EDC is heated up to approximately 670 K. The processes inside the furnace radiation section are strongly affected by the endothermic cracking reactions in the reactor coil and vice versa. The heat flux profile along the reactor tubes forms the connection between the process and flue gases. However, the experimental investigation of this profile is extremely difficult because of the long cycle periods, high cost, and limited measuring capabilities. Therefore, an accurate description of the heat transfer processes inside the furnace radiation section and a reasonable allocation of the energy within the furnace are crucial for the design and optimization of cracking furnaces. Over the past decades, significant progress in modeling hydrocarbon or chloride cracking furnaces has been made. Historically, these models include zero-, one-, and multidimensional models as well as computational fluid dynamics (CFD) models. In 1939, a zero-dimensional model called the Lobo− Evans model was developed by Lobo and Evans.1 The radiation section of the furnace is divided into the flue gas zone, heat © 2012 American Chemical Society

absorbing plane, and heat reflecting plane. The combustion is supposed to be perfectly mixed under constant temperature. However, this model is not accurate if the temperature of the combustion greatly varies along the height. Thus, Lobo2 developed a thermal gradient method based on the onedimensional Lobo−Evans method in 1974. The furnace is divided into several zones along the height, and each zone is calculated in the Lobo−Evans method. Afterward, multidimensional methods such as the zone,3 Monte Carlo,4 and heat flux methods5 were presented and applied in the furnace heat flux calculation. Froment,6 Vandamme and Froment,7 Huang et a1.,8 Xie et al.,9 and Xu et al.10 proposed several mathematical models. However, these simplified models cannot provide a detailed description of the processes inside the cracking furnace. CFD has recently become increasingly more popular in processing complicated heat and flow problems in furnaces. Similar to EDC crackers, ethylene crackers have been simulated using CFD. Taking into consideration the thermal cracking reactions in reactor coils, Detemmerman and Froment11 and Heynderickx et al.12 joined a zone model and CFD to perform coupled simulations for radiative heat transfer in an ethylene cracking furnace. Stefanidis et al.13−16 performed a detailed simulation of the combustion and radiation mechanisms in a furnace using CFD. More recently, Hu et al.17 used CFD to obtain detailed information on the velocity, temperature, and composition of flue gas in a naphtha cracking furnace. Only a few studies have focused on the simulation of ethylene chloride crackers using CFD. Park et al.18 used a twodimensional CFD model to calculate flow patterns and heat release based on equations of continuity, momentum, and Received: Revised: Accepted: Published: 645

September 10, 2012 December 1, 2012 December 13, 2012 December 13, 2012 dx.doi.org/10.1021/ie302436r | Ind. Eng. Chem. Res. 2013, 52, 645−657

Industrial & Engineering Chemistry Research

Article

⎧Q + σA [T 4 − T 4] − Q − ΔH − Q 0 gm2 1 gm1 r1 l1 ⎪ s1 ⎪ =0 ⎪ ⎪Q s2 + σA 0[Tgm14 − Tgm2 4] + σA 0[Tgm3 4 − Tgm2 4] ⎪ ⎪ − Q r 2 − ΔH2 − Q l 2 = 0 ⎪⋮ ⎪ ⎪Q + σA [T 4 4 0 gm(i − 1) − Tgmi ] + σ si ⎪ ⎪ 4 4 ⎨ A 0[Tgm(i + 1) − Tgmi ] − Q ri − ΔHi − Q li = 0 ⎪ ⎪⋮ ⎪ 4 4 ⎪Q s(N − 1) + σA 0[Tgm(N − 2) − Tgm(N − 1) ] + σ 4 4 ⎪ A [T 0 gm(N ) − Tgm(N − 1) ] − Q r(N − 1) − ΔH(N − 1) ⎪ ⎪ − Q l(N − 1) = 0 ⎪ 4 4 ⎪Q + σA [T 0 gm(N − 1) − Tgm(N ) ] − Q rN − ΔHN ⎪ sN ⎪ −Q =0 ⎩ lN

energy in an EDC cracking furnace, with the P1 model as the radiation heat model. However, the flow heat flux and temperature of the combustion and flue gas in the furnace were not provided. Generally, CFD simulation is commonly used in industrial equipment, especially in crackers. However, the computational load is an important factor in CFD simulation, making this method unsuitable for optimization. In this work, an industrial MISTUI EDC cracker with 100 kt/a capacity is examined. The cracker is simulated using onedimensional Lobo−Evans and CFD models. In one hand, optimization is performed using the first model for its timesaving property. First, fuel gas allocation operator α is investigated to improve performance indices such as selectivity, conversion, and unit fuel gas consumption (per VCM production). Latter, the optimum coil outlet temperature (COT) is also suggested to make good compromise among these performance indices by adjusting the total fuel gas flow rate. In the other hand, CFD simulation is used to validate the accuracy of the Lobo−Evans model and provide more details on the optimization results. A premixed combustion mechanism is adopted to account for the sidewall firing condition. Then, the discrete ordinate (DO) model is used to calculate the radiative heat transfer in the furnace. The weighted-sum-ofgray-gases model (WSGGM) is used to calculate the absorption coefficient of the mixture of gases. Furthermore, the tube skintemperature profile calculated in the one-dimensional platform is used to initialize that of the CFD model, thus accelerating the convergence speed of the coupled CFD simulation between the furnace and the reactor. As a result, the flow, combustion, heat transfer, and thermal reaction are provided by the CFD simulation. The temperature, velocity, and concentration distributions of the flue gas in the furnace and the heat flux profile on the reactor tube skin are obtained.

with five parts in consideration: the heat released by the fuel gas Qsi, the net radiation heat brought in/out from the axial neighboring zone, the heat absorbed by the reactor tubes Qri, the net heat of flue gas enthalpy change ΔHi, and the furnace wall heat dissipation Qli. 2.2. Reactor Model. A one-dimensional plug flow reactor model is assumed based on a high Reynolds number and a low viscosity of the reaction side stream. The simulation of the reactor integrates a set of continuity equations for the process gas species, along with the energy and momentum equations together. These equations originally proposed by Plehiers et al.19 can be written as dFj dz

2. MATHEMATICAL MODEL IN ONE-DIMENSIONAL LOBO−EVANS METHOD

i

πd t 2 4

πd t 2 dT ≡ Q (z ) π d t + dz 4

(4)

∑ rri(−ΔH )i i

(5)

⎛ 1 Pt ⎞ dpt ⎞ d⎛ 1 ⎞ 1 ⎜⎛ 1 dT − = + Fr ⎟ ⎜ ⎟ ⎜ ⎟+ 2 ⎠ dz ⎝ M m ⎠ M m ⎝ T dz η·G RT ⎠ dz ⎝ MmPt

2.1. Furnace Model. On the basis of the clear temperature gradient along the height inside the EDC cracker furnace, the radiation section of the furnace is divided into several zones along the height of the furnace. It is assumed that the fuel gas sprayed out of each burner row is totally consumed in its local zone.

(6)

with the friction factor: Fr = 0.092

ζ Re−0.2 + πR b dt

(7)

For the tube bends, the factor ζ is used for the pressure drop supplement given by Nekrasov.20

(1)

⎛ dt ⎞ Λ ⎞⎛ ⎟⎜0.051 + 0.19 ζ = ⎜0.7 + 0.35 ⎟ ⎝ 90° ⎠⎝ Rb ⎠

where Qsi is the heat release of fuel gas in the ith zone, Gsi is the fuel gas mass flow sprayed in the ith zone, and Qsl is the low heat value of the fuel gas. The heat profile absorbed by the reactor tubes in any furnace zone consists of radiation and convective heat transfers as defined below. Q ri = σaDAcpi Fi(Tgmi 4 − Twi 4) + hRcARi (Tgmi − Twi)

= (∑ sijrri)

∑ Fjcpj j

Q si = GsiQ sl

(3)

(8)

where Rb and Λ represent the tube bend radius and bend angle, respectively. The radical chain mechanism of the thermal cracking is wellknown in principle. Experimental studies and rigorous mathematical models of the complex EDC cracking reaction can be found in the literature.21−26 However, it is complicated and time-consuming to realize in coupled EDC cracker simulation. In this paper, a first-order successive molecular reaction model simplified from these radical chain mechanisms

(2)

According to one-dimensional Lobo−Evans method, the energy balance equation of the furnace could be written as 646

dx.doi.org/10.1021/ie302436r | Ind. Eng. Chem. Res. 2013, 52, 645−657


Article

is considered to describe the thermal cracking kinetics mechanism of EDC. k1

k2

EDC → VCM → byproduct

According to the thermal cracking kinetics presented above, two serried reactions are found in the tube with the following components: EDC, VCM, hydrogen chloride (HCl), and the byproduct represented by acetylene (C2H2). The kinetic equation is described as ki = Ai exp(−Ei/(RT)). The preexponential factors A1 and A2 are 5.507 × 107 and 5.580 × 107 s−1, respectively, and the reaction activation energy factors E1 and E2 are 123.86 and 137.78 kJ/mol, respectively, which are all regressed by the industrial data. 2.3. Coupled Numerical Simulation between the Furnace and the Reactor. The heat is strongly coupled between the furnace and the reactor. Hence, the heat flux transferred from the furnace to the reactor must be balanced with that absorbed by the process gas along the reactor length. In this work, the coupled numerical solution of onedimensional Lobo−Evans method is realized in Matlab(2009a), and the flue gas temperature and the heat flux profiles are considered as the core variables. The boundary conditions of the inlets in the furnace and the reactor are determined according to Table 1. The details of the flow diagram are illustrated in Figure 1.

Figure 1. Flow diagram of the coupled solution algorithm of the onedimensional Lobo−Evans model.

Table 1. Furnace Dimension and Operating Conditions Furnace Segment length (m) (x-direction) width (m) (y-direction) height (m) (z-direction) number of wall burners Firing Conditions fuel gas flow rate (kg/h) oxygen excess (vol%) fuel composition (wt%) C4H8 C4H6 C4H10 C3H8 H2,C1,C2 Material Properties emissivity of the furnace wall emissivity of tube skin Reactor Coils number of reactor row number of stright tubes (per row) number of bends (per row) inlet tube diameter (10−3 m) outlet tube diameter (10−3 m) thickness of tube (10−3 m) EDC feed rate (t/h) coil inlet temperature (K) coil outlet pressure (kPa)

3.1.1. Flow model. In the steady state, the governing equations in the furnace are established based on the Reynoldsaveraged Navier−Stokes (RANS) equations. The continuity, momentum, and energy transport equations in the furnace are expressed as

20.898 1.900 6.700 136 999.3 3

∂ (ρUi) = 0 ∂xi

76.0% 20.0% 2.6% 0.7% 0.7%

(9)

⎡ ⎤ ∂Uj ∂peff 2 ∂U ⎞ ∂ ⎢ ⎛⎜ ∂Ui ∂ − δij l ⎟⎟⎥ (ρUU + + μeff ⎜ i j) = ∂xi 3 ∂xl ⎠⎥⎦ ∂xi ∂xj ⎢⎣ ⎝ ∂xj ∂xj + SMi ⎛ ∂ ∂ ⎜ ∂T [Ui(ρE + p)] = − keff ∂xi ∂xj ⎜⎝ ∂xj

0.75 0.85 2 20 19 101.3 114.3 6.5 42 533.15 2355

(10)

∑ hjJj ⃗ j

⎡⎛ ∂U ⎤⎞ ∂Ui ⎞ ∂U j ⎟⎟ − 2 l δij ⎥⎟ + Sh + + Uiμeff ⎢⎜⎜ ⎢⎣⎝ ∂xi ∂xj ⎠ 3 ∂xl ⎥⎦⎟⎠

(11)

where the source term SMi in the momentum equation includes gravitation (SMx = SMy = 0 and SMz = −ρg), and the source term Sh in the energy equation includes the heat of chemical reaction and any other volumetric heat sources. For the compressible turbulent flow that exits in the furnace, two standard k−ε equations with a high Reynolds number in the turbulence model are adopted to close the governing eqs 9−11. The turbulent kinetic energy and its dissipation rate are expressed as

3. MATHEMATICAL MODELS IN THE CFD METHOD 3.1. Furnace Model. Compared with the one-dimensional Lobo−Evans model, the heat transfer model, combustion mechanism, turbulence, and flue gas flow are considered in the CFD model.

μ ⎞ ∂k ⎤ ∂ ∂ ⎡⎢⎛ (ρkUi) = ⎜μ + t ⎟ ⎥ + Gk − ρε ∂xi ∂xj ⎢⎣⎝ σk ⎠ ∂xj ⎥⎦ 647

(12)


Industrial & Engineering Chemistry Research μ ⎞ ∂ε ⎤ ∂ ∂ ⎡⎢⎛ ε (ρεUi) = ⎜μ + t ⎟ ⎥ + (c1Gk − c 2ρε) ∂xi ∂xj ⎢⎣⎝ σε ⎠ ∂xj ⎥⎦ k

Article

2H 2 + O2 → 2H 2O CO + 0.5O2 → CO2

(13)

The mixture is assumed to obey the ideal gas law. The viscosity, thermal conductivity, and specific heat of the mixture are computed from the properties of individual species, which are all functions of temperature. 3.1.3. Radiation Model. Radiation is the predominant mode of heat transfer in the cracking furnace because of the high temperature of the flue gas and that of the furnace wall. Hence, an appropriate radiation model is strongly needed. Several wellknown radiation models, such as P-1, P−N, Rosseland, discrete transfer, and discrete ordinates (DO), are usually used to simulate radiation heat transfer. The DO model is probably the best suited for radiative transfer solutions with localized heating sources, especially in a cracking furnace. It considers the radiation effects of nongray gases and works across the full range of optical thicknesses. Moreover, the scattering and particle effects between gas and particulates can be finely processed in the DO model. This model also has a broader application range and is easily solved simultaneously with the coupled calculation of the fluid flow and energy equations. In this study, the radiative transfer equation is discretized in the DO model for a finite number of solid angles associated with the vector direction s ⃗ fixed in the global Cartesian coordinate system (x,y,z). The finite-volume scheme and its extension are used to solve unstructured mesh problems. The general DO model is expressed as

In addition, the turbulence model equations are invalid for wall-bounded flow regions and close to walls. Thus, the standard wall functions based on the proposal of Launder and Spalding27 to the wall boundary conditions of the furnace and reactor are applied to bridge the viscosity-affected region between the wall and the fully turbulent region. 3.1.2. Combustion Model. In this study, the fuel and air are perfectly premixed and sprayed out to combustion. On the basis of the work of Magnussen and Hjertager,28 finite-rate chemistry combining the eddy dissipation concept in the turbulence-chemistry interaction model is applied to describe the combustion in the EDC furnace. The chemical source terms are computed in the finite-rate model using Arrhenius expressions without considering the turbulent fluctuations. The Arrhenius molar creation/destruction rate R̂ i,r of species i in reaction r is expressed as N

R̂ i , r = Γ(νi″, r − νi′, r )(k f , r ∏ [Cj , r ](ηj′,r + ηj″,r ) ) j=1

(14)

Given that all the flames in the furnace are premixed, the mixture of fuel gas and air will burn as soon as they enter the computational domain. The minimum of the Arrhenius and the eddy dissipation reaction rates are considered in the calculation of the combustion and turbulence models. In the eddy dissipation model, the eddy dissipation reaction rate that also indicates the production net rate of species i caused by reaction r, Ri,r, is given by the minimum in following equation: R i,r

dI ( r ⃗ , s ⃗ ) n 2 σT 4 = −(K a + K s)I( r ⃗ , s ⃗) + K a ds π 4π Ks + I( r ⃗ , Ω′) Φ(Ω, Ω′) dΩ′ 4π 0

∫

⎛ ⎛ Y ⎞ ε R ⎟⎟ , = min⎜vi′, rM w , iAρ min⎜⎜ ⎜ k ⎝ vR′ , rM w , R ⎠ ⎝

vi′, rM w , iABρ

⎞ ε ∑P YP ⎟ k ∑Nj v″j , rM w , j ⎟⎠

(16)

The furnace is filled with a large amount of flue gas components (mostly CO2 and H2O), which should participate in the radiation heat transfer. Given the low scattering effect of flue gas, the scattering coefficient Ks usually is naturally set at zero. Only the absorption coefficient of the flue gas Ka must be calculated to close the radiation transfer equation. The components of the flue gas have different absorption bands. Hence, WSGGM is used to describe the radiation characteristics of flue gas mixtures. The emissivity of the flue gas is basically assumed as the weighted sum of the emissivities of a number of gray gases, and a transparent gas over the distance s is expressed as30

(15)

To simulate precisely the combustion in an EDC furnace, the use of the full reaction mechanism, such as the radical mechanism, is indispensable for considering all the intermediates. However, a practical simulation including a threedimensional flow with the full reaction mechanism is beyond the capability of present computers. Therefore, the most realistic solution is to adopt a simplified reaction mechanism. In this study, the two-step reaction mechanism of Westbrook and Dryer29 with CO as the intermediate is used because it is one of the most commonly used in CFD applications. Considering the fuel used in the trial whose composition is listed in Table 1, the C1, C2, and H2 contents in the fuel gas are relatively small and not explicitly accounted for, which slightly affects the fuel gas combustion. In our simulations, these contents in the fuel gas are assumed to be replaced by H2. The simplified two-step mechanism is described as

I

ε=

∑ αε ,i(T )(1 − e−k ps) i

(17)

i=0

The absorption coefficient is derived from following equation when s > 10−4 m. Ka = −

ln(1 − ε) s

(18)

For the DO radiation model, the furnace wall is considered as the diffusion reflection plane. The boundary intensity of the furnace wall for all outgoing directions s ⃗ at the wall is given by

C3H8 + 3.5O2 → 3CO + 4H 2O C4 H6 + 3.5O2 → 4CO + 3H 2O

Im,w =

C4 H8 + 4O2 → 4CO + 4H 2O

(1 − εw ) π

4π

∫s ⃗·n⃗ >0

Iins ⃗ ·nd⃗ Ω + εw n2

σTw 4 π

(19)

3.2. Reactor Model. In the CFD simulation, similar to that used for the furnace, a reactor requires the integration of a set

C4 H10 + 4.5O2 → 4CO + 5H 2O 648



Article

flow rate requirement of each burner, not all the burners are turned on. The turning on/off status of the burners is also shown in Figure 3, where the burners marked with black circles are turned off and the switch status of the burners in the neighboring rows is mutually staggered. The CFD simulation of an EDC cracker is divided into the thermally coupled reactor and furnace models. In the reactor model, hexahedral cells are used to describe the tube wall of the straight parts, whereas tetrahedral cells are used to describe the bends. In the furnace model, the hexahedral cells are used to discretize the furnace zone, whereas the tetrahedral cells are used to discretize the physical domain including the burner and tube skin zones. A total of 263 540 and 1 645 382 grid cells are present in the reactor tubes and furnace, respectively. According to the operating parameters in Table 1, the inlet boundary conditions of the fuel gas, air in the burners, and the EDC raw material of the reactor inlet are determined as the mass flow inlet boundary. Thus, the mass flow, temperature, and mixture composition of the fuel gas and air are given, whereas the mass flow and temperature of the reactor tube inlet are also determined. The outlet boundary conditions of the furnace and reactor tube are all determined as pressure outlet boundaries. The former is set to −50 Pa gauge at the outlet of the furnace, and the latter is defined as the coil outlet pressure. The heat boundary on the furnace walls is calculated on the basis of the industrial heat loss value, which is 3% of the total heat supply. The CFD simulation of EDC cracker is realized in ansys14.0.0(fluent). The coupled calculation algorithm is shown in Figure 4. In the simultaneous numerical solution,

of continuity equations for the product yields, along with energy and momentum equations, which are calculated from RANS equations with the k-ε model for turbulence as the closure model. If a complicated reaction mechanism similar to that in the radical model described by Song et al.26 is used, thousands of variables would be involved in the overall dimension reactor calculation. Currently, this calculation is too computational to be accepted in the CFD simulation. In addition, the simulations between the one-dimensional Lobo− Evans and CFD platforms should be coherent. Therefore, the reaction mechanism described as the first-order successive molecular reaction model in the one-dimensional Lobo−Evans model is adopted as the furnace model. A total of 11 balance equations should be solved to obtain the velocity, temperature fields, pressure fields, and product yield profiles along the reactor length. 3.3. Grid Division, Boundary Conditions, and Numerical Method. Only half of the industrial EDC cracker is simulated because of the geometric similarity and symmetry along the width. A schematic representation of the simulated industrial EDC cracker is illustrated in Figure 2. The furnace is

Figure 2. Front view of the considered half segment of an EDC cracker.

heated with 68 radiation burners lined in four rows on each side of the long wall. The arrangement of the burners and the closeup of a burner are shown in Figure 3. According to the fuel gas

Figure 4. Flowchart of the coupled solution algorithm for the mathematical CFD models in the furnace and the reactor.

several iterations are needed between the two models to reach convergence. First, to hasten the convergence, an external tube skin-temperature profile calculated from the Lobo−Evans platform is assigned as the tube skin boundary of the furnace model using a user defined function (UDF). Second, the combustion details and heat transfer in the furnace simulation, as well as the first heat flux profile boundary condition for the reactor tubes, are obtained. Next, the transfer and reaction models of the reactor tube are calculated using a new reactor

Figure 3. The arrangement of the burners and the close-up of a burner. 649



Article

tube skin-temperature profile. The coupled calculation is reiterated until the maximal difference of tube skin temperature between two successive iterations reaches a predefined threshold value, which is set at 1 K in this coupled simulation. The governing equations in the furnace and reactor tubes are solved sequentially using a control volume-based method. The convection items are discretized implicitly through a secondorder upwind scheme. The semi-implicit method for the pressure-linked equation algorithm is adopted to solve the coupled momentum, energy, and species transport equations.

4. RESULTS AND DISCUSSION The coupled EDC cracker is simulated based on the simplified one-dimensional Lobo−Evans method and CFD. The former is relatively simple for calculation. Thus, the optimization problem is solved in this platform. However, it cannot provide the combustion and flue gas details in the furnace, so the optimized results are crosschecked using CFD. Although the calculation is much more time-consuming in CFD, it provided the combustion reaction mechanism and turbulence, temperature, and mass concentration of the flue gas in the furnace. 4.1. Fuel Gas Allocation Optimization Result and Analysis. The EDC cracker has a series of performance indices, including EDC cracking conversion, selectivity, and fuel gas consumption (per unit VCM production). The operational conditions include the air/fuel ratio, furnace pressure, and fuel gas allocation. In the EDC cracking industry, the fuel gas and the air are premixed with a fixed air/fuel ratio and the furnace pressure is controlled at −50 Pa gauge. These conditions ensure that the fuel gas is fully combusted and that no CO emission exists at the furnace outlet. Only the fuel gas allocation at all the burners has not been discussed yet. Given that the reactor and the furnace are seriously thermally coupled, the temperature of the fuel gas, the heat flux absorbed by the reactor tube, and the temperature, cracking conversion, and selectivity of the process gas are strongly affected by the fuel gas allocation. Therefore, finding the best solution to the fuel gas allocation is profitable. With the MISTUI EDC cracker discussed in this paper, the fuel gas allocated in the four rows is controlled by two valves. The fuel gas mass flow in the fourth row at the low position of the sidewall is controlled by one valve, whereas the fuel gas mass flow of the rest is controlled by another. Thus, the fuel gas allocation operator α is defined as the ratio of the fuel gas mass flow in the fourth flow path to the total fuel gas mass flow. In this section, the influences of the fuel gas allocation operator α are discussed below. In practice, fuel gas allocation operator is limited in the range of [0.1, 0.5], avoiding some risks such as flameout, damages of firebricks and tubes. However, to enrich the analysis data, the range is released to [0.1, 0.8] in the simulation. The average flue gas temperature curves of the different fuel gas allocations in the furnace are demonstrated in Figure 5. When the fuel gas allocation operator α expands from 0.1 to 0.8, the temperature of the bottom of the furnace increases from 927 to 1134 K, whereas that at the top of the furnace decreases from 959 to 837 K. The temperature of the flue gas is more uniform when α is 0.3 (Figure 5). When α is below 0.3, the temperature increases at the bottom part of the furnace and decreases at the top along the furnace height. By contrast, the temperature decreases monotonically when α is above 0.3. Additionally, Figure 5 shows that all the curves at the different fuel gas allocations almost intersect at one point at the second furnace zone. The fuel gas allocation scheme, limitation of the one-

Figure 5. Fuel gas temperature profiles along the furnace partition.

dimensional Lobo−Evans model, and inadequate furnace-zone division are believed to be responsible for this phenomenon. The heat flux profiles of tube skin along the reactor tube length are shown in Figure 6a. When α increases, which indicates that less fuel gas is released in the upper three rows and more fuel gas is released in the last row, the heat flux in the first half pass of the reactor deceases, whereas that in the last half pass of the reactor increases. The heat flux along the reactor length can be partitioned into four parts, and the heat flux profiles in each part decreases monotonously (Figure 6a). This phenomenon is caused by the limitation of the onedimensional Lobo−Evans model, where the furnace is divided into four parts. In this model, the temperature of each part is supposed to be uniform, and the temperature of the process gas increases monotonously along the reactor length. As a result, the temperature gradient between the flue and process gases becomes smaller along the reactor length. The heat flux profile is the most uniform when α = 0.3 (Figure 6a). The process gas absorbs the heat flux transfer from the tube wall to increase the temperature for the endothermic reaction. The temperatures of the reactor tube skin and process gas are shown in Figure 6 panels b and c, respectively. When α is small, the temperatures of the tube skin and process gas rapidly increase in the front half pass. In the last half pass, their temperatures increase relatively slowly and even decrease. The opposite condition occurs when the operator α is large. The COT increases from 737.5 to 777.6 K as the operator α varies from 0.1 to 0.8 Figure 6b), although the process gases in all the simulations are supposed to be preheated up to 533.15 K at the reactor inlet. The process gas composition varies with the cracking reaction along the reactor tube length (Figure 6d). The EDC cracking reaction begins when the temperature of the process gas reaches 670 K. When the operator α is small, the ratio of the fuel gas mass flow located at the top of the furnace is large. During this period, the temperature of the process gas increases rapidly, and the EDC cracking reaction happens earlier and more violently during the first half pass. As the reaction proceeds, the EDC concentration decreases along the reactor tube length, and the increase in the EDC cracking rate is limited by the concentration of EDC. These phenomena are more obvious when α is large. Several critical performance indices in the EDC cracker are taken into consideration in the industry. The variations in EDC cracking conversion, EDC cracking selectivity, and fuel gas 650



Article

Figure 7. (a) EDC cracking conversion, (b) EDC cracking selectivity, (c) fuel gas consumption per unit VCM production.

monotonously increases from 0.5518 to 0.5791. The EDC cracking selectivity of the process gas outlet with operator α variation is illustrated in Figure 7b. As α increases from 0.1 to 0.8, the selectivity decreases from 0.9511 to 0.9427. The fuel gas consumption (per unit VCM production) variation is shown in Figure 7c. When α increases from 0.1 to 0.8, the fuel gas consumption (per unit VCM production) initially increases slowly from 71.82 to 71.98 and then decreases to 69.03. The first minor increase in Figure 7c is believed to be caused by the limitation of the simplified one-dimensional Lobo−Evans model. Generally, the fuel gas consumption (per unit VCM production) should decrease monotonously with increasing α.

Figure 6. (a) Heat flux profiles, (b) tube skin-temperature profiles, (c) process gas temperature profiles, and (d) process gas composition along the reactor tube length.

consumption (per unit VCM production) with the fuel gas allocation operator α are shown in Figure 7. The EDC cracking conversion varies in the process gas outlet with the operator α (Figure 7a). When α increases from 0.1 to 0.8, the conversion 651



Article

With the increase in the operator α, the conversion of EDC increases, but the selectivity of EDC and the fuel gas consumption (per unit VCM production) decrease. As previously mentioned, the flue gas temperature and the heat flux profiles along the reactor length are more uniform when α is close to 0.3. This condition is believed to prolong the run length of the EDC cracker. Moreover, the minimum selectivity is suggested to be 0.948 to prevent excessive EDC loss and high energy consumption in the downstream. Considering all the constraints mentioned above and the field experience, the fuel gas allocation operator α is most profitable at 0.365 to make good balance of conversion, selectivity, fuel gas consumption (per unit VCM production). 4.2. COT Optimization and Analysis. Given that the original patent data of the MISTUI cracker is invalid after the addition of heat-exchange drum for EDC feed vaporization using the heat of process gas, no process gas component analysis instrument is implemented in situ. The COT of the reactor is the only measured and controlled variable of the EDC cracking performance in the industry. Thus, quantitative analysis of the relationship between COT and the cracking performance indices is particularly necessary. When the optimized fuel gas allocation operator α is implemented, the COT changes along with the total fuel gas mass flow rate. The simulation result on the changes in COT with varying fuel gas consumption is illustrated in Figure 8. As the fuel gas mass flow increases from 510 kg/h to 1500 kg/h, the COT almost linearly increases from 710 to 796 K.

Figure 8. Changes in the coil outlet temperature with varying fuel gas consumption.

The changes in the outlet process gas composition, conversion, selectivity, and fuel gas consumption (per unit VCM production) with COT variation are shown in Figure 9. The cracking severity is deepened with increasing COT (Figure 9a). In the coil outlet, the molar flow rate of residual EDC decreases from 158.5 kmol/h to 48.7 kmol/h, whereas those of VCM, HCl, and the byproduct increase from 52.8, 54.4, and 0.8 kmol/h to 146.1, 180.7, and 17.1 kmol/h, respectively. The EDC cracking conversion monotonously increases from 0.2528 to 0.7703 with a 0.066 increase in a 1 K increase in the COT (Figure 9b). By contrast, the EDC cracking selectivity monotonously decreases from 0.9848 to 0.8941 with a 0.09 decrease in a 1 K increase in the COT (Figure 9c). The fuel gas

Figure 9. Changes with COT variation: (a) process gas outlet composition, (b) EDC cracking conversion, (c) EDC cracking selectivity, and (d) fuel gas consumption per unit VCM production.

consumption (per unit VCM production) as a parabola approximation curve with the increase in COT is depicted in Figure 9d. First, it decreases from 79.11 to 71.075 kg/T and then increases to 81.27 kg/T. The minimum point of the fuel 652



Article

gas consumption (per unit VCM production) is found at the COT of 742 K. To satisfy the selectivity constraints mentioned above and guarantee a high conversion and a low fuel gas consumption (per unit VCM production), the COT and fuel gas should be 756 K and 993.5 kg/h, respectively. 4.3. Details and Discussion in the CFD Simulation. To validate the accuracy of the one-dimensional Lobo−Evans model and provide details on the combustion, turbulence, temperature, and mass concentration of the flue gas, the optimized results are resimulated in the CFD platform again. The details of the reactor are shown in Figure 10, in which Figure 10a depicts the profiles of the external tube skin-, internal tube skin- and process gas temperatures along the reactor length. During the first half pass of the reactor, the process gas temperature rapidly increases because the absorbed heat is mainly used to increase the process gas temperature. The EDC cracking reaction occurs when the process gas temperature reaches 670 K. The process gas temperature slowly increases during the last half pass because most of the absorbed heat is used for the endothermic reaction. Furthermore, Figure 10a also illustrates that the temperature of the external tube skin is approximately 10 K higher than that of the internal tube skin, which is approximately 40 K higher than the temperature of the process gas. The profile of the process gas composition molar flow rate along the reactor length is depicted in Figure 10b. Given the low temperature of the process gas in the first 140 m length of the reactor, almost no cracking reaction happens at this point. When the process gas temperature reaches 670 K, the main EDC cracking and side reactions are both accelerated. In addition, the molar flow rate of EDC decreases, whereas those of VCM, HCl, and the byproduct relatively increase. The variation in process gas velocity along the reactor tube length is illustrated in Figure 10c. During the first half pass, the volume flow rate is hastened because the volume of the process gas expands with increasing temperature. During the last half pass, the volume flow rate of the process gas increases more obviously than that during the first half pass. The reason is that the cracking reaction is the decomposition reaction and the total molar flow rate of all process gas species tremendously increases along the reactor length. The process gas pressure distribution along the reactor length is shown in Figure 10d. The coil inlet pressure of the process gas is 2904 kPa, and it continuously drops along the reactor length until the coil outlet pressure reaches 2413 kPa. The simulated pressure drop coincides with that of the industry. The z-velocity distribution of the flue gas at different planes in the furnace is illustrated in Figure 11. The z-velocity distribution of the fuel gas at different planes (i.e., x = 1, 5, 10, 15, and 19 m) is shown in Figure 11 panels a, b, c, d, and e, respectively. The premixed fuel gas/air mixture sprays out from the sidewall burners and swiftly combusts to become flue gas. As in the close-up of a burner depicted in Figure 3, the nozzles are tilted toward the sidewall, so the high-speed flue gas jets are ejected to the sidewall at high speed. Afterward, the downward and upward flue gas jets collide and diffusively migrate toward the reactor plane. Hence, the z-velocity close to the reactor plane is much slower than that close to the sidewall, resulting in a more uniform high-temperature distribution and more residence time of the flue gas in the furnace. The flue gas velocity close to the furnace radiation outlet is higher than the average velocity in the furnace because the space size of the furnace radiation section outlet becomes smaller (Figure 11a− e). The z-velocity distribution of the flue gas at plane y =

Figure 10. (a) Process gas composition profile along the reactor length, (b) temperature profile along the reactor length, (c) process gas velocity along the reactor length, and (d) process gas pressure along the reactor length.

0.0485 m is illustrated in Figure 11f. Considering the high fuel gas allocation ratio in the fourth row, the velocity of flue gas 653



Article

Figure 11. Z-velocity distribution at different planes: (a) x = 1 m, (b) x = 5 m, (c) x = 10 m, (d) x = 15 m, (e) x = 19 m, and (f) y = 0.0485 m.

close to the fourth row is higher than that close to the first three rows. The complexity of the flue gas velocity distribution determines the flue gas temperature distribution in the furnace. The flue gas temperature contours at different sections along the furnace lengths of 1, 5, 10.5, 15, and 18.5 m are shown in Figure 12 panels a, b, c, d, and e, respectively. A clear temperature gradient is found along the height of the furnace. The temperature of the flue gas close to the furnace bottom is high, whereas that close to the furnace top is low because more fuel gas is allocated to the burners at the last row. The flue gas temperature contours at different sections along the furnace widths of 0.0485, 0.097, and 0.71285 m are depicted in Figure 12f−h, respectively. The premixed fuel gas/air mixture is ejected from the sidewall burner, fast turns to high-temperature combustion and collides with the furnace wall. The reactor is mainly heated by the radiation from the high-temperature wall. Thus, a clear temperature gradient can be found along the width of the furnace as it is shown in Figure 12f−h. Besides that, as y moves, the flue gas mixes better, the high-temperature spots become lower, and the temperature of flue gas becomes more evenly distributed (Figure 12f−h). The O2 mass concentration contours at different sections along the furnace lengths of 1.5, 6, 10.5, 15, and 18.5 m are illustrated in Figure 13panels a, b, c, d, and e, respectively. Given that the fuel gas and air are premixed, the O2 mass fraction decreases because of the combustion of the fuel gas (Figure 13a−e). As a result, the O2 mass fraction is almost equal to 0.03 in most areas of the furnace, except at the jet zones. The O2 mass concentration contours at the plane y = 0.0485 m, which is close to the sidewall, are demonstrated in Figure 13f. Although the fuel gas and air are premixed, the combustion fails to be fully immediately conducted as the

Figure 12. Flue gas temperature contour at different planes in the furnace: (a) x = 1.5 m, (b) x = 5 m, (c) x = 10.5 m, (d) x = 15 m, (e) x = 18.5 m, (f) y = 0.0485 m, (g) y = 0.097 m, and (h) y = 0.71285 m.

mixtures spray out of the burners because of the high velocity and low temperature of the fuel gas/air mixture. Therefore, high O2 concentration spots are formed close to the sidewall burners. The fuel gas is mainly composed by C4s, which has a poor thermal stability and can be thermally decomposed and oxidized in a very short distance. Given that the combustion of the fuel gas in the furnace is a limited diffusion combustion, an incomplete combustion of the fuel gas generates CO. The CO mass concentration contours at different sections along the furnace lengths and widths of x = 1.5 m, x = 5 m, x = 9.5 m, x = 15 m, x = 18.5 m, and y = 0.0485 m are shown in Figure 14 panels a, b, c, d, e, and f, respectively. The CO mass concentration rapidly decreases to zero because ambient excess air is joined to the combustion to guarantee that the fuel gas and CO are completely consumed. The termination position of the combustion flame is known to be the end of the oxidation reaction. Hence, the flame area can be estimated by the distribution of CO contents. The CO mass concentration at plane y = 0.0485 m is shown in Figure 14f, which is also the plane close to the sidewall burners. The high CO mass concentration spots depict the flame area. The CO2 mass concentration contours at different sections along the furnace lengths and widths of x = 1.5 m, x = 5 m, x = 9.5 m, x = 15 m, x = 18.5 m, and y = 0.0485 m are depicted in Figure 15 panels a, b, c, d, e, and f, respectively. The premixed 654



Article

Figure 13. O2 mass concentration contours at different planes: (a) x = 1.5 m, (b) x = 6 m, (c) x = 10.5 m, (d) x = 15 m, (e) x = 18.5 m, and (f) y = 0.0485 m.

Figure 15. CO2 mass concentration contours at different planes: (a) x = 1.5 m, (b) x = 5 m, (c) x = 9.5 m, (d) x = 15 m, (e) x = 18.5 m, and (f) y = 0.0485 m.

finding may be attributed to the inadequate combustion of the fuel gas caused by the high-velocity and low-temperature jet of the burners, such that a substantial amount of intermediate CO is generated instead of CO2. 4.4. Comparison among the Two Platforms and the Industry. Through the simplified one-dimensional Lobo− Evans simulation, the fuel gas allocation is optimized, and the most profitable COT is suggested at full capacity. Meanwhile, the optimized result is validated with the CFD method. Using this method, the details on the combustion in the furnace are obtained, which cannot be obtained either in the onedimensional Lobo−Evans simulation or in the industry. Finally, the optimized fuel gas allocation operator a = 0.365 and the suggested COT 756 K are implemented to the MISTUI EDC cracker in the industry. On the basis of Table 2, the simulation results in the one-dimensional Lobo−Evans platform and in the CFD platform show high agreement with the industrial data.

5. CONCLUSION A coupled reactor/furnace simulation of an industrial MISTUI EDC cracker with 100 kt/a capacity is performed for the first time using the one-dimensional Lobo−Evans and CFD approaches. Optimization is performed using the one-dimensional Lobo− Evans approach because it is relatively simple and time saving. The heat distribution on the furnace is seriously influenced by the fuel gas allocation in the furnace, which also affects the crucial performance indices, such as EDC cracking conversion, selectivity, and fuel gas consumption (per VCM production). The fuel gas allocation operator α is investigated to improve these indices, and the best solution is found at α = 0.365. Furthermore, with the fuel gas allocation operation factor α settled at 0.365, the coil outlet temperature is optimized at 756

Figure 14. CO mass concentration contours at different planes: (a) x = 1.5 m, (b) x = 5 m, (c) x = 9.5 m, (d) x = 15 m, (e) x = 18.5 m, and (f) y = 0.0485 m.

fuel gas/air mixture is consumed in the combustion, and a large amount of CO2 is generated and diffused in the furnace (Figure 15a−e). The CO2 mass concentration has a small gradient and a uniform distribution in most parts of the furnace (i.e., close to 0.17), except for that in the burner’s outlet. A low CO2 mass concentration is found at the burner’s outlet (Figure 15f). This 655



Article

Acpi = cold plane, m2 aD = angle factor ARi = tubes’ skin surface in the ith zone, m2 B = empirical constant c1,c2 = standard k − ε model constants, c1 = 1.44 and c2 = 1.92 Cpj = heat capacity of the jth species, J/mol/K D = mass diffusion coefficient, m2/s dt = tube diameter, m E = total energy per unit mass, J/kg E1E2 = reaction activation energy, kJ/mol E(L) = equivalent conversion coefficients in the location L of the tube Fi = radiation heat transfer coefficient Fj = molar flow rate, mol/s Fr = friction factor G = the mass flow rate of the process gas Gk = generation of turbulent kinetic energy, J/m3/s Gsi = fuel gas mass flow in the ith zone, kg/h hrc = convective heat-transfer coefficient ΔHi = net heat of flue gas enthalpy change W I = radiation intensity, J/m2/s Iin = intensity of the incoming ray, J/m2/s Jj⃗ = diffusion flux of species j, kg/m2/s K = heat transfer coefficient Ka = absorption coefficient, 1/m Ks = scattering coefficient, 1/m keff = effective conductivity, W/m/K k1,k2 = pre-exponential factors, 1/s ki = absorption coefficient of the ith gray gas, 1/m L = reactor length, m Mm = average molecular weight of all the process gas species, kg/mol Mw,i = molecular weight of species i, kg/mol n = refractive index n⃗ = normal pointing out of the domain N = number of chemical species in the system Ni = molar flow rate of the ith component p = pressure and the sum of the partial pressures of all absorbing gases, Pa qir = radiation heat transfer, J/m3 Qli = the furnace wall heat dissipation in the ith zone W Qsl = low heating value of the flue gas, J/kg Qsi = heat release of fuel gas in the local zone, W Qri = heat absorbed by the reactor tubes located in ith zone, W pt = total pressure, Pa p′ = fluctuating pressure, Pa peff = effective pressure, Pa rri = reaction rate, mol/m3/s r ⃗ = position vector R = ideal gas constant, R = 8.314 J/mol/K Rb = radius of the reactor bend, m Ri = net rate of production of species i by chemical reaction, mol/m3/s Re = Reynolds number s = path length s ⃗ = direction vector sij = stochiometry factor S = cross sectional area of the reactor tube m2 Sh = source term in the energy equation, J/m3/s T = local temperature and process gas temperature, K Tgmi = average temperature of combustion in the ith zone, K

Table 2. Comparison of Simulation Results and Industrial Data item temp of flue gas at outlet (K) temp of process gas at outlet (K) coil inlet pressure (kPa) coil pressure drop (kPa) performance indices conversion selectivity fuel gas consumption (per unit VCM)/ kg/t

industrial data

optimized data in the Lobo−Evans platform

simulation data in the CFD Platform

900.0

913

894.5

756.0

756.0

753.6

2904

2953

2904

549

598

491

0.550 0.950 72.0

0.554 0.949 71.7

0.547 0.945 72.5

K to make good compromise among all the performance indices by adjusting the fuel gas mass flow rate. The CFD simulation platform is used to crosscheck the simplified model. It also provides a proper theoretical basis for gaining better insight in the behavior of the industrial cracker. To accelerate the convergence speed of the coupled CFD simulation between the furnace and reactor, the tube-skin temperature calculated in the one-dimensional Lobo−Evans platform is used to initialize that of CFD. The flow, combustion, heat transfer, and thermal reaction are provided by the CFD simulation. The temperature, velocity, and concentration distributions of the flue gas in the furnace and the heat flux on the reactor tube skin are precisely obtained, which are useful for the optimal design of EDC cracking furnace/coils and for the positioning and turning on/off statuses of burners. However, the coupled CFD simulation requires a number of days of simulation time. Furthermore, the optimized fuel gas allocation operator α = 0.365 and suggested COT 756 K are implemented in the industry. After a steady long run time on full capacity in industry, the data obtained from the industry shows high agreement with the optimization results in the one-dimensional Lobo−Evans and CFD platforms.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: 86-21-64252060 (F.Q.), 86-21-64251250 (W.Z.). Fax: 86-21-64252060 (F.Q.), 86-21-64251250 (W.Z.). E-mail: [email protected] (F.Q.), [email protected] (W.Z.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is supported by Major State Basic Research Development Program of China (2012CB720500), National Natural Science Foundation of China (Key Program: U1162202), National Natural Science Foundation of China (General Program: 61174118), and Shanghai Leading Academic Discipline Project (B504).

■

NOMENCLATURE α = fuel gas allocation operator A = empirical constant A0 = effective axial radiation surface, m2 656



Article

(12) Heynderickx, G. J.; Oprins, A. J. M.; Marin, G. B.; Dick, E. Three-dimensional flow patterns in cracking furnaces with long-flame burners. AIChE J. 2001, 47 (2), 388−400. (13) Stefanidis, G. D.; Merci, B.; Heynderickx, G. J.; Marin, G. B. CFD simulations of steam cracking furnaces using detailed combustion mechanisms. Comput. Chem. Eng. 2006, 30 (4), 635−649. (14) Stefanidis, G. D.; Heynderickx, G. J.; Marin, G. B. Development of reduced combustion mechanisms for premixed flame modeling in steam cracking furnaces with emphasis on NO emission. Energy Fuels. 2006, 20 (1), 103−113. (15) 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. (16) Stefanidis, G. D.; Van Geem, K. M.; Heynderickx, G. J.; Marin, G. B. Evaluation of high-emissivity coatings in steam cracking furnaces using a non-grey gas radiation model. Chem. Eng. J. 2008, 137 (2), 411−421. (17) Hu, G. H.; Wang, H. G.; Qian, F. Numerical simulation on flow, combustion and heat transfer of ethylene cracking furnaces. Chem. Eng. Sci. 2011, 66, 1600−1611. (18) Park, Y.; Choi, B. S.; Yi, J. Simulation of imbalance reduction between two reactors in an ethylene dichloride cracker. Chem. Eng. Sci. 2005, 60, 1237−1249. (19) 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. (20) Nekrasov, B. B. Hydraulics; Peace Publishers: Moscow, 1969. (21) Barton, D. H. R. Kinetics of the Dehydrochlorination of substituted hydrocarbons. Part I. Induced dehydrochlorination. J. Chem. Soc. 1949, 148−155. (22) Barton, D. H. R.; Howlett, K. E. Kinetics of the dehydrochlorination of substituted hydrocarbons. Part II. The mechanism of the thermal decomposition of 1:2-dichloroethane. J. Chem. Soc. 1949, 156−164. (23) Borsa, A. G. Industrial Plant/Laboratory Investigation and Computer Modeling of 1,2-Dichloroethane Pyrolysis. Ph.D. Dissertation, Colorado School of Mines, Golden, CO, 1999. (24) Huybrechts, G.; Wouters, G. Mechanism of the pyrolysis of 1,2dichloroethane in the absence and presence of added chlorine. Int. J. Chem. Kinet. 2002, 34 (5), 316−321. (25) Ashmore, G. P.; Owen, A. J.; Robinson, P. J. Chlorine-catalysed pyrolysis of 1,2-dichloroethane. Part 1. Experimental results and proposed mechanism. J. Chem. Soc., Faraday Trans. 1982, 78, 657− 676. (26) Song, B. H.; Choi, B. S.; Kim, S. Y.; Yi, J. Kinetics and mechanism of pyrolysis of the 1,2-dichloroethane. J. Korean Inst. Chem. Eng. 2001, 99, 481−487. (27) Launder, B. E.; Spalding, D. B. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 1974, 3, 269−289. (28) Magnussen, B. F.; Hjertager, B. H. On mathematical models of turbulent combustion with special emphasis on soot formation and combustion. Proc. 16th Int. Combust. Symp. 1976, 719−729. (29) Westbrook, C. K.; Dryer, F. L. Simplified reaction-mechanisms for the oxidation of hydrocarbon fuels in flames. Combust. Sci. Technol. 1981, 27, 31−43. (30) Siegel, R.; Howell, J. R. Thermal Radiation Heat Transfer; Hemisphere publishing Corporation: Washington, DC, 2002.

Twi = average temperature of tube outer skin in the ith zone, K Ui,Uj,Ul = velocity component in the ith, jth, or lth direction, m/s v′ = fluctuating velocity, m/s V = volume flow velocity of the process gas, m3/s xi, xj, xl = coordinate direction in the ith, jth, or lth direction, m YP = mass fraction of product P YR = mass fraction of reactant R Greek Letters

αε,i = emissivity weighting factors for the ith fictitious gray gas β = convective heat-transfer coefficient of the process gas λw = thermal conductivity of the tube wall ε = dissipation rate of turbulent kinetic energy, m2/s3 and emissivity δij = Kronecker delta ρ = gas density, kg/m3 ρs = density of species s, kg/m3 Γ = net effect of third bodies on the reaction rate μt = turbulent viscosity, kg/m/s μeff = effective viscosity, kg/m/s νi,r′ = stoichiometric coefficient for reactant i in reaction r νi,r″ = stoichiometric coefficient for product i in reaction r η = unit conversion factor η′j,r = rate exponent for reactant species j in reaction r ηj,r″ = rate exponent for product species j in reaction r ζ = supplementary friction factor for the reactor bend Λ = tube bend angle, deg σ = Stefan−Boltzmann constant, σ = 5.672 × 10−8 W/m2 K4 σε = standard k−ε model constant σε = 1.3 Φ = phase function Ω = cross sectional area of local reactor tube/solid angle Ω′ = solid angle

■

REFERENCES

(1) Lobo, W. E.; Evans, J. E. Heat transfer in the radiant section of petroleum heaters. Trans. AIChE 1939, 35, 743. (2) Lobo, W. E. Design of furnaces with fuel gas temperature gradients. Chem. Eng. Progress 1974, 70 (1), 65−71. (3) Hottel, H. C. Radiation Heat Transmission, 3rd ed.; McGraw Hill Book Company: New York, 1959. (4) Hottel, H. C.; Sarofim, A. F. Radiative Transfer, 3rd ed.; McGraw Hill Book Company: New York, 1967. (5) Roesler, F. C. Theory of radiative heat transfer in co-current tube furnace. Chem. Eng. Sci. 1967, 22, 1325−1338. (6) Froment, G. F. Thermal cracking for olefins production fundamentals and their application to industrial problems. Chem. Eng. Sci. 1981, 36 (8), l271−1282. (7) Vandamme, P. S.; Froment, G. F. Scaling up of naphtha cracking coils. Ind. Eng. Chem. Process Des. Dev. 1981, 20, 366−376. (8) Huang, S. N.; Qian, F.; Yu, J. S. Zone method used for ethylene cracking furnace modelling. J. East China Inst. Chem. Technol. 1992, l 8 (suppl.), 160−166. (9) Xie, D. M.; Qian, F.; Yu, J. S. The mathematical model for naphtha pyrolysis in SRT-III cracking furnace I. The modelling and calculation of clean tube. Petrochem. Technol. 1993, 22 (12), 813−819. (10) Xu, Q.; Chen, B. Z.; He, X. R.; Zhang, L. Simulation for naphtha pyrolysis in clear radiation tube of SRT-IV cracking furnace. Comput. Appl. Chem. 2001, 18 (3), 223−228. (11) Detemmerman, T.; Froment, G. F. Three-dimensional coupled simulation of furnaces and reactor tubes for the thermal cracking of hydrocarbons. Oil, Gas, Sci. Technol,. Rev. IFP. 1998, 53 (2), 181−194. 657


Comprehensive Simulation and Optimization of an Ethylene

Recommend Documents