Simulation of the Evaporation Tube Banks in the Convection Section

Sep 5, 2017 - Simulation of the evaporation process of a naphtha feedstock in the convection section of a steam cracking furnace is reported in this a...
0 downloads 10 Views 4MB Size
Article pubs.acs.org/IECR

Simulation of the Evaporation Tube Banks in the Convection Section of a Steam Cracking Furnace Using an Evaporation Model Benfeng Yuan, Yu Zhang, Xin Peng, Wenli Du, Guihua Hu,* and Feng Qian* State Key Laboratory of Chemical Engineering, East China University of Science and Technology, 130 Meilong Road, Shanghai 200237, China. Key Laboratory of Advanced Control and Optimization for Chemical Processes, Ministry of Education, East China University of Science and Technology, Shanghai, China, 200237 ABSTRACT: Simulation of the evaporation process of a naphtha feedstock in the convection section of a steam cracking furnace is reported in this article. A two-step evaporation model was developed to determine component changes in the vapor and liquid phases. First, a feedstock-simplifying method, based on the work of Eckert and Vaněk (Comput. Chem. Eng. 2005, 30, 343−356), was developed in which the original hydrocarbon mixture is represented by a system of real components. Second, the component distribution was determined by vapor− liquid equilibrium theory, using an activity coefficient method. The group contribution method (GCM) was employed to calculate the interaction parameters in the activity coefficient method. Using a mechanistic heat-transfer method, a comparison was made between the evaporation model developed in this work and a simplified evaporation model. The simulation results are in good agreement with industrial data. Some parameter profiles for the heat-transfer process are provided, as are the flow patterns of the naphtha feedstock. flue-gas flow velocity and temperature fields in the convection chamber. Mahulkar et al.13,14 and De Schepper et al.15 analyzed the coking mechanism and simulated the coke-forming process of heavy feedstocks in the convection section. Some suggestions for methods to reduce coke formation were also made on the basis of these studies. Hu et al.16 took the temperature change of the feedstock into account when determining the physical properties of both the vapor and liquid phases. However, the above-mentioned simulation methods have drawbacks including the facts that (1) simulations using the CFD technique are extremely time-consuming and, as a consequence, have limited application in the optimization and control of steam cracking furnaces and (2) variations in the components of the vapor and liquid phases cannot be determined by the CFD technique. These variations have a significant impact on the physical properties of the vapor and liquid phases, including both thermodynamic properties [critical temperature (Tc), critical pressure (Pc), critical volume (Vc), and acentric factor (ωac)] and transport properties [kinematic viscosity (ν), thermal conductivity (λ), specific heat capacity (Cp), and vapor pressure (Pvapor)], which are very important for heattransfer calculations in flow boiling processes. Thus, a fast and accurate evaporation model of the convection section is necessary for broad application in chemical industries.

1. INTRODUCTION The steam cracking furnace is the key unit for the production of ethylene in the chemical industry. It mainly consists of two parts: a radiation section and a convection section. In the steam cracker, the flue gas that leaves the radiation section enters the convection section, which contains several horizontal tube bundles. The flue gas flows through the tube bundles, preheating and vaporizing the hydrocarbon feedstock inside the tubes so as to recover the remaining energy of the flue gas. Considering that the efficiency of the cracking process is crucial for the product yield and the economics of the steam cracker, previous efforts to improve steam cracking efficiency have usually focused on the radiation section.1−9 With the increasing demands for ethylene and propylene and the decreasing supplies of light feedstocks, inexpensive but heavier feedstocks, such as heavy naphtha and heavy vacuum gas oil (HVGO), are being cracked more frequently. Apparently, this change has the potential to cause more fouling problems within the tubes in the convection section. In addition, more heat is needed to ensure that the heavier feedstocks evaporate completely and are able to reach the reaction temperature at the end of the convection section. To analyze and solve these problems, researchers need an accurate description of the evaporation process in the convection section. Over the past few decades, significant progress has been made in the modeling of the convection sections of steam crackers, mostly based on the computational fluid dynamics (CFD) method. De Schepper et al.10−12 conducted a CFD study on the feedstock temperature field inside tubes and the © 2017 American Chemical Society

Received: Revised: Accepted: Published: 10813

July 9, 2017 September 3, 2017 September 5, 2017 September 5, 2017 DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research Many researchers have focused on the heat transfer during flow boiling.17−19 Among previous studies, a mechanistic model was developed by Wojtan et al.20,21 In this model, correlations were constructed between flow patterns and corresponding heat-transfer coefficients. The correlations for the flow-pattern map depend on the physical properties of the vapor and liquid phases and the vapor quality. To calculate accurate physical properties, an evaporation model that considers component changes is necessary. Verhees et al.22 developed a one-dimensional model for the coupled simulation of a steam cracker convection section based on Wojtan et al.’s correlations.20,21 In this model, a set of pseudocomponents was employed for the characterization of a multicomponent hydrocarbon feedstock. In addition, an evaporation model, based on the assumption that the pseudocomponents evaporate one by one, was used. The evaporation order of these pseudocomponents was determined by boiling points higher than the saturation temperature at a given vapor pressure. Although component changes and flow patterns were taken into account, some disadvantages of this model cannot be ignored. First, the widely accepted pseudocomponent method23,24 relies on pseudocomponents that are primarily defined by their (pseudo) boiling points and some additional parameters, mostly density, molecular weight, and/or viscosity. However, for the determination of other thermodynamic and transport properties, some predictive methods are needed. In particular, the group contribution method, which requires detailed molecular-structure information to estimate interaction parameters, cannot be used for the calculation of vapor−liquid equilibria in this context. In 2001, Eckert made a comparison of molecular weights, critical temperatures, critical pressures, and acentric factors of 98 pure components and studied the reliability of the methods for estimating physical properties in the pseudocomponent method.25 The results showed that the estimation methods mostly led to larger errors. Second, because of interactions among the molecules, the species in an actual mixture cannot evaporate one by one, as assumed in Verhees et al.’s model.22 Thus, the evaporation mechanism in Verhees et al.’s work,22 which is referred to as the “simplified model” in this article, causes certain errors in the simulation results. Vapor−liquid equilibrium (VLE) theory is the common evaporation mechanism used to describe multicomponent evaporation and is widely used in distillation processes.26−28 In this work, a two-step evaporation model is developed for the simulation of the flow boiling of naphtha feedstocks in the convection section. In the first step, the naphtha feedstock is characterized by several real components, and a feedstocksimplifying model is developed based on Eckert and Vaněk’s approach.29 In the second step, vapor−liquid equilibrium theory, with some predicted interaction parameters, is used to calculate the component distribution during evaporation. On the basis of coupled simulations between the tube section and the convection chamber, the flow-pattern-based heattransfer method proposed by Wojtan et al. was employed to calculate the heat-transfer rate and predict flow patterns during flow boiling. Consequently, comparative testing of the evaporation mechanism using Verhees et al.’s model22 and that described in this article was undertaken. Vapor qualities, heattransfer coefficients, feed-fluid temperatures, external tube-skin temperatures, flue-gas temperatures, and external tube heat fluxes in the evaporation bank of the convection section are presented. Moreover, flow-pattern maps for the naphtha feedstock are given.

Figure 1. Flow regimes during evaporation in a horizontal tube.

2. MATHEMATICAL MODEL 2.1. Heat Balance. The heat transferred to a fluid in a tube section of length ΔL is calculated as Q = htotalP ΔL(Tg − T )

(1)

where htotal is the total heat-transfer coefficient, which is defined as htotal

−1 ⎛ 1 d d 1 d ⎞⎟ ⎜ =⎜ + ln + 2hw D hf D ⎟⎠ ⎝ htp

(2)

For the flow of a single vapor/liquid phase, the increase in the fluid temperature over the tube section of length ΔL is caused by the forced convection of the single vapor/liquid phase

Q = FCp̅ (Tout − Tin)

(3)

When evaporation occurs, a vapor phase forms inside the tube. The compositions of the two phases vary during evaporation. The following heat conservation equation was used to calculate the outlet temperature during flow boiling in a section of tube [Fe in ∑ yC + F(1 − e in) ∑ xiC Pl, i](Tout − Tin) i Pv, i i

i

= Q + F(eout − e in)ΔHlatent

(4)

2.2. Heat-Transfer Model. 2.2.1. Internal Heat Transfer. Figure1 depicts several typical flow regimes encountered during the flow boiling process in a horizontal tube, where distinct heat-transfer mechanisms are needed to address different flow patterns. The heat-transfer correlations developed by Wojtan et al.20,21 are listed in Table 1. In Wojtan et al.’s work,20,21 the fact that part of the heat is directly transferred into the vapor phase on the dry upper perimeter of the tube during stratified flow (S), as well as stratified-wavy flow [including slug flow, slug/stratified-wavy flow, and stratified-wavy (SW) flow], is taken into account. Therefore, the average heat-transfer coefficient is a combination of the vapor-phase and liquid-phase heattransfer coefficients, resulting in eqs 5−9. Under various flow regimes, the dry angle, θdry, has different forms as follows: For slug flow, intermittent flow, and annular flow, θdry = 0. For stratified flow, θdry = θstrat. ⎛ Gwavy − G ⎞0.61 For stratified-wavy flow, θdry = ⎜ G − G ⎟ θstrat . ⎝ wvy strat ⎠ 0.61 x ⎛ Gwavy − G ⎞ For slug/stratified-wavy flow, θdry = x ⎜ G − G ⎟ θstrat IA ⎝ wavy strat ⎠ For slug flow, intermittent flow, and annular flow, the entire perimeter of the tube is saturated with the liquid flow. In this case, heat is transferred only into the liquid phase, and the equations that calculate the heat-transfer coefficients for the three flow regimes are identical. During annular flow with partial dryout (D), heat is transferred directly into the vapor phase in the dry region, and the heat-transfer coefficient decreases dramatically and then levels out at the value for mist (M) flow. Hence, the heat-transfer 10814

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research Table 1. Tube Interior Heat-Transfer Correlation Equations for the Model of Wojtan et al.20,21 S, SW, Slug, Slug + SW, I, and A Flow Regimesa

h tp =

θdryhv + (2π − θdry )hwet

(5)



hv = 0.023Re v 0.8Prv 0.4 3

kv D

(6)

3 1/3

hwet = (hnb + hcb )

hcb = 0.0133Reδ 0.69PrL 0.4 hnb = 55Pr

θstrat

(7) kL δ

−0.55

0.12

( −log Pr)

(8) M

−0.5 0.67

q

(9)

⎧ ⎛ 3π ⎞1/3 = 2π − 2⎨π(1 − ε) + ⎜ ⎟ [1 − 2(1 − ε) + (1 − ε)1/3 − ε1/3] ⎝ 2 ⎠ ⎩ ⎪



⎫ − (1 − ε)ε[1 − 2(1 − ε)]{1 + 4[(1 − ε)2 + ε 2]}/200⎬ ⎭





ε=

(10)

⎧ ⎛x x⎪ 1 − x⎞ ⎜⎜ + ⎟⎟ ⎨ + − [1 0.12(1 x )] ρv ⎪ ρL ⎠ ⎝ ρv ⎩ 1.18(1 − x)[gσ(ρL − ρv )]0.25 ⎫ ⎪ ⎬ + 0.5 ⎪ GρL ⎭

−1

(11)

Mist-Flow Regime

kv (12) D Annular Flow with Partial Dryout x − xdi = h tp(xdi) − [h tp(xdi) − hmist(xde)] xde − xdi hmist = 0.0117ReH 0.79PrV1.06Y −1.83

hdryout a

(13)

A, annular; I, intermittent; S, stratified; SW, stratified-wavy.

is provided by the equation31

coefficients can be calculated from eq 12, by linear interpolation, where htp(xdi) is the heat-transfer coefficient of the vapor quality at the dryout inception point, xdi, calculated by eq 5, and hmist(xde) is the heat-transfer coefficient of the vapor quality at the mist inception point, xde, calculated by eq 12. 2.2.2. External Heat Transfer. For the coupled simulation of evaporation tube banks in the convection section, the heat transfer in the convection chamber also needs to be considered. Because the evaporation tube banks are located on the top of the convection section, where the flue-gas temperature is low at about 373−573 K, radiative heat transfer can be neglected. The mass velocity of the flue gas increases significantly when the flue gas passes through the tube bundle. Reference velocities for calculations of flow and heat transfer, for both in-line and staggered arrangements, are given by the equation30 ⎛ φ ⎞ φT T Umax = max⎜⎜ Uapp , Uapp⎟⎟ φD − 1 ⎝ φT − 1 ⎠

hf = 0.1378

ρdUmax μ

(16)

2.3. Evaporation Model. 2.3.1. Characterization of Feedstock Using Real Components. The characterization of a heavy feedstock using real components is similar to the approach of the pseudocomponent method. Unlike the physical properties of pseudocomponents, the physical properties of real components, such as thermodynamic and transport properties, can be obtained easily. Normally, the definition of a system of real components for a complex mixture is a two-step process: First, a set of substitute components is selected, and then the composition of this mixture is determined. The substitute components are first selected using characterization curves, in particular, true-boiling-point (TBP) curves. As only ASTM D86 distillation curves are usually available in the industry, the TBP curve of the naphtha feedstock in this work was calculated by transforming the ASTM D86 curve using an empirical correlation.32 As in the pseudocomponent method, the boiling point range of the TBP curve must be cut to obtain nonoverlapping temperature intervals (Ti, Ti+1, for i = 1, ..., I) that cover the entire temperature range. For each temperature interval, a set of candidate real components is selected from a database on the basis of their boiling points within the temperature interval. Hence, a sufficiently large naphtha molecular database that contains as many components of a naphtha feedstock as possible is required for accurate calculations. However, the detailed composition of a naphtha

(14)

where φD = φL 2 + (φL/2)2 is the dimensionless diagonal pitch in the case of the staggered arrangement. The Reynolds number for any arrangement is defined by ReD =

⎛ d p ⎞0.296 kG ReG 0.718PrG1/8⎜ ⎟ Db ⎝ l ⎠

(15)

Because finned annulus tubes are used in evaporation banks, a correlation for calculating the heat transfer in finned tube banks 10815

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research

Before eq 18 is used, the representative component for each temperature interval should be selected. Thus, eq 17 is still needed. However, if the optimization procedure needs to be implemented for each characterization, it is complicated and time-consuming. If one can confirm that the tendency that the small differences in the common properties or some of these properties between the temperature interval and the selected components correspond to the small differences in the thermodynamic and transport properties, then only the common properties need to be considered, and the calculation process can be greatly simplified. To demonstrate this ideal, ωk in eq 17 is redefined as the weighting factor and reflects the influence of the common properties on the thermodynamic and transport properties. Here, ωk is in the range of 0−1 and ∑ωk = 1. In this work, only the TBP and specific density (sd) curves were obtained. Hence, only two weighting factors, ωTBP and ωsd, need to be determined. For calculation accuracy, the best choice is to use thermodynamic and transport properties measured by experiment for comparison with those of the selected components. However, because it is difficult to measure thermodynamic and transport properties, some predictive methods for the basic properties of the pseudocomponents were employed.35 Although this is a compromise solution, compared with the pseudocomponent method, this approach allows the thermodynamic and transport properties of the real component to be easily obtained, which further simplifies the calculation process and gives similar results. The results show that, when ωTBP exceeds 0.9 and ωsd is within 0.1, the thermodynamic and transport properties of the selected real component are the closest to the predicted values. In other words, when the normal boiling point of the selected component is closer to the mean temperature of the temperature interval, the thermodynamic and transport properties of the selected component are closer to those of the pseudocomponent that is represented for the temperature interval. Thus, selecting the real component based on a simple comparison of the boiling point is sufficiently precise for these calculations, and the influence of the specific density can be ignored. The measured and theoretical (defined by real components) TBP and specific-density curves of the naphtha feedstock are compared in Figure 2, and some thermodynamic properties (i.e., critical pressure, Pc, and critical temperature, Tc) of the selected real components are depicted in Figure 3. After the representative components have been established, the composition distribution of these components needs to be determined. This work used the optimization method proposed by Eckert and Vaněk.29 As shown in Figure 4, each real component represents a distilled fraction interval, and its normal boiling point is the mean temperature of that interval. However, the intervals generally cannot be chosen such that they cover the entire range without gaps and/or overlaps. An optimization task is therefore established where the objective function can be defined as

feedstock cannot be obtained in a straightforward manner. Thus, in this work, SIMCO software was used to predict the detailed composition of a naphtha feedstock through the application of the molecular reconstruction method.33,34 SIMCO has a power molecular library that was obtained from many real naphthas by experimental methods. Some industrial indexes of the naphtha feedstock are taken as input of this software, as listed in Table 2. Table 2. Industrial Indexes of the Naphtha Feedstocka index

value

specific density (g/cm3) 20 °C 15 °C ASTM D86 boiling points (K) IBP 10% 50% 90% 95% FBP PIONAb characterization (wt %) n-paraffins i-paraffins olefins naphthenes aromatics

0.715 0.7196 318.15 342.15 365.15 434.15 442.15 464.15 30.16 41.56 1.37 16.13 10.78

IBP, initial boiling point; FBP, final boiling point. bPIONA, paraffin, isoparaffin, olefins, naphthenes, and aromatics. a

After the candidate components have been determined, one component should be selected from the candidate components for the representation of each temperature interval; as a consequence, some comparisons are needed. To accurately capture the characteristics of the original mixture, the physical properties of the selected component should be as close as possible to those of the corresponding temperature interval. Thus, a weighted sum of the relative differences of the various properties was used for selection. The criterion is defined by K

∑ ωk

|ξr, k ,c − ξm, k ,c| ξm, k ,c

k=1

→ minc

(17)

In Eckert and Vaněk’s work, ωk is the weighting factor that reflects the precision of the measurements. ξm,k,c and ξr,k,c are the measured properties and the properties of the candidate components, respectively. The most common characterization curve used in the industry is the TBP curve, although the specific density curve is used occasionally. Unfortunately, other physical properties, such as thermodynamic and transport properties, that are crucial for heat-transfer calculations cannot be easily measured. Hence, simply comparing the common properties is not sufficient, and the small relative differences in the common properties calculated using eq 17 might not represent the small differences in other properties. To address this issue, another similar optimization task is established that can be defined by the objective function 29

l

∑ l=1

|ξm, l − ξr, l| ξm, l

I+1

F (ν L ) =

∑ (νiR− 1 − νiL)2 → min i=2

L

νi < νi L+ 1 ,

→ min (18)

i = 1, ..., I

νi − 1 < νi L < νi + 1 ,

In this equation, ξm,l and ξr,l are the thermodynamic and transport properties that have been measured and the thermodynamic and transport properties of the selected components, respectively.

i = 2, ..., I

(19)

In this equation, vLi and vRi denote the starting and ending distillation fractions, respectivelyy. To determine the relative 10816

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research Table 3. Compositions of Selected Components as Determined by Optimization

Figure 2. Measured (squares) and theoretical (solid lines) TBP (left) and specific density (right) curves. The theoretical curves were obtained using the substitute mixtures.

mass fraction (%) 4.13 1.22 4.72 18.48 14.74 8.14 16.82 7.17 12.91 10.82 0.85

To avoid such a demanding task, an alternative solution is employed using a predictive method that can estimate the interactions from mere knowledge of the structures of molecules. In particular, the PPR78 model36 was employed, which combines the Peng−Robinson(PR) EoS37 and a group contribution method (GCM).36,38 This approach can be seen as a GCM for estimating the binary interaction parameter (BIP) for the classical mixing rule of the PR78 EoS, as shown Table 4. Through the use of eqs 21−23, φvi , which is temperaturedependent, can be predicted. In eq 20, the activity coefficient of the liquid phase is also required. Several models for calculating activity coefficients have been proposed, such as the Wilson, universal quasichemical (UNIQUAC), and nonrandom two liquid (NRTL) models.39 As was the case for the calculation of BIPs, some interaction parameters in these methods need to be fitted using vapor−liquid equilibrium data. Unfortunately, experimental data are hard to obtain. This work used the universal functional activity coefficient (UNIFAC) model to estimate activity coefficients.40 In the UNIFAC model, the activity coefficient is calculated as the sum of combinatorial and residual contributions, described by

Figure 3. Measured (squares) and theoretical (solid lines) critical temperature (Tc) and critical pressure (Pc) as functions of TBP. The theoretical curves were obtained using the substitute mixtures.

ln γi = ln γic + ln γi R

Figure 4. Construction of distilled-fraction intervals with real components.

(24)

where

proportions of the selected components, the gap and/or overlap between the adjacent intervals needs to be as small as possible, and this is achieved through optimization of the interval distribution. The composition distribution of the selected components is reported in Table 3. 2.3.2. Calculation of Vapor−Liquid Equilibrium. Vapor− liquid equilibrium theory was applied to calculate the variations in the components of the two phases. For the mixture, because of molecular interactions, the real system cannot be taken as an ideal system. Thus, two parameters, namely, the vapor-phase fugacity (φvi ) and liquid activity (γi) coefficients, are introduced. A well-known vapor−liquid equilibrium equation with the two parameters was employed, namely φi v yP = xiγiPis i

component neopentane isopentane cyclopentane n-hexane 2,3-dimethylheptane 2,4-dimethyloctane 2,2,4-trimethylnonane n-nonane 3-methyldecane n-undecane n-dodecane

ln γic = ln

φi xi

+1−

φi xi



φ⎤ Z ⎡ φi + 1 − i⎥ ⎢ln θi ⎦ 2 ⎣ θi

(25)

and ln γic =

∑ vki[ln Γk − ln Γik]

(26)

k

The molecular volume fraction (φi) and the surface-area fraction (θi) are defined in terms of the segment numbers r and the constant parameter numbers q as xiqi xiri φi = and θi = ∑j xjrj ∑j xjqj (27) In eq 26, Γk is the activity coefficient of group k at the mixture composition, and Γik is the activity coefficient of group k at a group composition corresponding to pure component i. Γk and Γik are given by:

(20)

In this equation, φvi denotes the pressure correction of the true gas to the ideal gas and γi denotes the mole fraction correction of the true solution to the ideal solution. The vapor-phase fugacity coefficient is usually computed from an equation of state (EoS). However, for the mixture, accurately quantifying the interactions between each pair of molecules is difficult.

ln Γk = 10817

⎡ z ⎢ Q k 1 − ln(∑ θmτmk) − 2 ⎢⎣ m

∑ i

⎤ θτ i ki ⎥ ∑j θτ j ji ⎥ ⎦

(28)

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research Table 4. Equations of the PPR78 Model PR78 Model

⎛ bi A ⎜ 2 (Zm − 1) − ln(Zm − B) − 2 2 B ⎜⎝ am bm ⎡ Z + ( 2 + 1)B ⎤ ln⎢ m ⎥ ⎣ Zm − ( 2 − 1)B ⎦

ln φi v =

N

∑ xj

aiaj −

j=1

⎞ bi ⎟ bm ⎟⎠

(21)

⎧ R = 8.314472 J mol−1 K−1 ⎪ RT ⎪bi = 0.0777960739 c, i ⎪ Pc, i ⎪ 2 2 ⎛ ⎞⎤ ⎨ R Tc, i 2 ⎡⎢ where T ⎟⎥ 1 + mi⎜⎜1 − ⎪ ai = 0.457235529 Pc, i ⎢⎣ Tc, i ⎟⎠⎥⎦ ⎪ ⎝ ⎪ if ω ≤ 0.491, m = 0.37464 + 1.5422ω − 0.26992ω 2 i i i ⎪ i 2 3 ⎩ if ωi > 0.491, mi = 0.379642 + 1.48503ωi − 0.164423ωi + 0.016666ωi Mixing Rule N N ⎧ ⎪ am = ∑ ∑ yy aiaj [1 − K ij(T )] i j ⎪ i=1 j=1 ⎪ N ⎪ ⎪bm = ∑ ybi ⎪ i ⎨ i=1 ⎪ ⎪ A = a mP ⎪ m (RT )2 ⎪ Pbm ⎪ ⎪ Bm = ⎩ RT

(22)

Group Contribution Method N N ⎧ ⎪ 1⎡ g g ⎛ 298.15 ⎞(Bkl − 1/ Akl)⎤ ⎥ ⎟ K ij = ⎨ − ⎢∑ ∑ (αik − αjk)(αil − αjl)Akl ⎜ ⎥⎦ ⎝ T /K ⎠ ⎪ 2 ⎢⎣ k = 1 k = 1 ⎩

⎛ ai(T ) −⎜ − ⎜ bi ⎝

aj(T ) ⎞⎫ ⎟⎪ ⎬ bj ⎟⎠⎪ ⎭

⎡ a (T )a (T ) j ⎢2 i ⎢ b b i j ⎣

(23)

of the flue gas. In this work, it is reasonable to assume that the flue-gas temperature in each tube layer is uniform. Hence, only one tube is selected to represent all of the tubes in each layer. 3.2. Solution Strategy. To explain the calculation algorithm of the simulation for each tube bank, this article takes the HTC-I bank as an example. As shown in Figure 6, each tube pass is divided into 10 subzones. Considering that thermal insulation ensures that heat losses are small at the elbows, it is acceptable that the tube bends are not taken into account in the heat-transfer calculations. At the start of the simulation, the flue-gas temperatures between each pair of adjacent layers are first initialized, and the flue-gas temperature in each tube layer is taken to be the average value of the flue-gas temperatures of the two adjacent sides. For example, the flue-gas temperature used to calculate the external heat-transfer coefficient for pass 2 is equal to the average value of Tg5 (the temperature of the flue gas entering pass 2) and Tg4 (the temperature of the flue gas exiting pass 2) (Figure 6). As a consequence, the heat transfer in each subzone is calculated according to eq 1. Finally, the total heat absorbed by each tube pass is obtained by summing the transferred heat in each subzone, and a new set of flue-gas temperatures between the different layers can be calculated from the heat balance of the flue gas. This iteration continues until the maximum differences in the flue-gas temperatures between two iterations are all below the predefined threshold value of 1 K. The algorithm used for this calculation is shown in Figure 7.

To calculate the detailed compositions of the two phases, a combination of the phase-equilibrium equation and the massbalance equation is necessary. The mass-balance equation is given by Fzi = Vyi + Lxi

⎤ ⎥ ⎥ ⎦

(29)

and the calculation process was implemented in AspenPlus.

3. SIMULATIONS 3.1. Tube Geometry and Operating Conditions. The configuration of the convection section of a steam cracker is shown in Figure 5. This section is composed of eight tube banks. The evaporation of the naphtha feedstock occurs predominantly in the feedstock-preheating (FPH) bank and the high-temperature-coil (HTC-I) bank. In addition to the FPH and HTC-I tube banks, the economizer (ECO) tube bank is also taken into account for a continuous temperature drop of the flue gas and a full heat-exchange network. The temperature and mass flow of the flue gas leaving the HTC-II tube bank, which are used from the design values, are imposed as the inlet conditions of the flue-gas side. The tubular structure and operating conditions are listed in Table 5. The CFD simulation results reported by Hu et al.16 indicated that the temperature distribution of the flue gas tends to be uniform throughout the width and length directions at higher elevations, because the staggered arrangement of the tube bundles enhances the mixing 10818

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research Table 5. Tube Geometry and Operating Conditions Convection Chamber Specifications length (z direction) (m) 10.368 FPH number of tube passes 20 number of tubes 6 tube outer diameter (m) 0.0603 thickness of tube (m) 0.00554 transverse pitch (m) 0.102

width (x direction) (m) Tube Bank fin thickness (m) fin number fin height (m)

mass flow rate of naphtha (kg/s) longitudinal pitch (m) 0.08833 inlet temperature (K) ECO Tube Bank number of tube passes 18 fin thickness (m) number of tubes 6 fin number tube outer diameter (m) 0.0603 fin height (m) thickness of tube (m) 0.00874 transverse pitch (m) 0.102 mass flow rate of boiler feedwater (kg/s) longitudinal pitch (m) 0.102 tube inlet temperature (K) HTC-I Tube Bank number of tube passes 6 thickness of tube (m) number of tubes 6 fin number tube outer diameter (m) 0.1143 fin height (m) transverse pitch (m) 0.201 longitudinal pitch (m) 0.203 mass flow rate of dilute steam (kg/s) fin thickness (m) 0.00127 inlet temperature of dilute steam (K) Flue-Gas Side mass flow rate (kg/s) 25.14 inlet temperature of the flue gas (K) flue-gas composition (mol %) O2 1.98 CO2 N2 71.07 Ar H2O 12.26

Figure 5. Schematic diagram of the convection section of a steam cracker.

4. RESULTS AND DISCUSSION 4.1. Comparison of VLE and Simplified Evaporation Mechanisms. To analyze the influence of the VLE model and the simplified model described by Verhees et al.22 on heat transfer, simulations were compared using a tube with a length of 300 m and a diameter of 0.0603 m. The naphtha feedstock enters the tube at an inlet temperature of 333.15 K and an inlet pressure of 5.8 bar. A heat flux of 15.836 kW/m2 is used so that the feedstock can be completely evaporated at the end. Figure 8 depicts the vapor quality profiles inside the tube as predicted by the VLE and simplified models. Using the simplified model, the evaporation process occurs earlier but has a lower evaporation rate than that predicted by the VLE model. Thus, the vapor quality profiles along the tube axis predicted by the two models are very different. The different profiles of the vapor quality also lead to differences in the mixture velocity profiles. As shown in Figure 9, because of an earlier evaporation, the mixture velocity predicted by the simplified model is obviously greater than that predicted by the VLE model in the first half of the tube. Because the vapor phase has a higher velocity than the liquid phase, this promotes the mixture velocity inside the tubes. In the latter half of the tube, the mixture velocities predicted by the two models are comparable, because the velocity of the vapor phase plays a dominant role in the mixture velocity. In general, the higher mixture velocity predicted by the simplified model will lead to a larger pressure drop than that predicted by the VLE model. Detailed flowpattern maps are presented in Figure 10. Vertical lines (x = xIA) separating the intermittent and annular flows at vapor qualities of 0.5 and 0.98 are included, along with a horizontal line labeled G495 that represents a mass flux of 495 kg m−2 s−1. It can be observed that, when the feedstock is determined, the

1.22 0.00127 157 0.01905 8.486 333.15 0.00127 157 0.01905 10.48 418.15 0.00602 0.01905 0.85 458.15

837.15

13.47 1.22

flow-pattern map is also determined. Thus, the flow patterns in the evaporation process mainly depend on the mass velocity. Moreover, because of variations in the physical properties of the vapor and liquid phases in the evaporation process, the vertical line that separates the intermittent and annular flow regimes changes with vapor quality. It can be observed that slug flow, intermittent flow, annular flow, and mist flow are observed in both evaporation processes. Although the flow-pattern maps produced by the two models are almost identical, it should be noted that the different vapor quality profiles result in different heat-transfer-coefficient distributions, as shown in Figure 11. Using the simplified model, the feedstock evaporates earlier, resulting in higher heat-transfer coefficients near the tube entrance than predicted by the VLE model. When the fluid enters the mist-flow region, both models predict an abrupt drop in heat-transfer coefficient, followed by a region in which the heat-transfer coefficient tends to be constant. Because the simplified model predicts a longer mist-flow region than the VLE model, it also predicts lower heat-transfer coefficients in the second half of the tube. These significant differences in heattransfer-coefficient profiles cannot be ignored, especially for the calculation of the external tube-skin temperature. Hence, the use of the VLE model to simulate the evaporation process in the evaporation bank of the steam cracker is necessary. 4.2. Application of the Evaporation Model to the Evaporation Bank. Utilizing the evaporation model developed in this work, the evaporation process of the naphtha feedstock in the evaporation tube banks was simulated. 10819

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research

Figure 8. Vapor quality profiles predicted by the VLE and simplified models.

Figure 6. Schematic diagram of the tube subzones of the HTC-I tube section.

Figure 9. Mixture velocity profiles predicted by the VLE and simplified models.

naphtha feedstock. Thus, the vapor quality at the inlet of the HTC-I bank is higher than that at the outlet of the FPH bank. Different evaporation flow patterns occur as a result of variations in vapor quality. Figure 13a shows the flow-pattern map of the FPH bank predicted at G650 (mass velocity of 650 kg m−2 s−1), which represents the mass velocity in the FPH bank. Such a high mass velocity and low vapor quality result in flow patterns that are located in the slug and intermittent flow regions. The flow-pattern map of the HTC-I bank is illustrated in Figure 13b. As some dilute steam is added to the feed fluid, the vapor quality increases and reaches 0.37. Moreover, because of a larger tube diameter, the mass velocity in this bank is only 189 kg m−2 s−1, which is much smaller than that in the FPH bank. Stratified-wavy flow, annular flow with partial dryout, and mist flow are predicted. Because of the low mass velocity, the flow pattern at the beginning of the tube is stratified-wavy flow. When the vapor quality exceeds xdi (the dryout inception point), annular flow with partial dryout occurs, resulting in an abrupt drop in the heat-transfer coefficient, as shown in Figure 14b. When the vapor quality exceeds xde (the mist-flow inception point), mist flow occurs, resulting in a constant heat-transfer coefficient. Heat-transfer coefficient profiles are depicted in Figure 14. In the FPH bank, the internal heat-transfer coefficient increases gradually along the tube axis, with a dramatic rise in zone 175, where the naphtha feedstock begins to evaporate. This is explained by the fact that, compared with the heat transfer in the single liquid phase, boiling heat transfer leads to a significant increase in the heat transfer inside the tube. Moreover, the

Figure 7. Flowchart of the algorithm used for the coupled simulation of the tubes and the flue-gas side.

Table 6 compares the simulation results with the design data. The simulated outlet temperatures of the three banks are in good agreement with the design data. To analyze the heat transfer in the evaporation process, some parameter profiles in the FPH and HTC-I tube banks are discussed based on Wojtan et al.’s heat-transfer correlations.20,21 Figure 12 illustrates the flow conditions in the evaporation bank. Evaporation mainly occurs in the HTC-I bank, with less than 10% of the naphtha feedstock evaporating in the FPH bank. Some dilute steam, whose mass flow accounts for 10% of the mass flow of the naphtha feedstock, is added to the naphtha feedstock in the HTC-I bank to promote the evaporation of the 10820

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research

Table 6. Comparison of Industrial Data and Simulation Results item

design data (K)

simulation data (K)

430 523.15 609

425.78 530.17 598.37

385.15 463.15 632.15

390.02 460.89 639.45

fluid FPH ECO HTC-I flue gas FPH ECO HTC-I

boundary layer are enhanced, resulting in a continuous increase in the internal heat-transfer coefficient. However, the internal heat-transfer coefficient is observed to drop abruptly in tube zone 15, after which it tends to be constant. This is explained by the fact that annular flow with partial dryout starts in this zone, as can be seen in Figure 13b and is finally transformed into mist flow. In addition, there is also an abrupt rise in zone 25, where the naphtha feedstock is evaporated completely and the flow regime is transformed from mist flow into single-phase vapor flow. As shown in Figure 14, the internal heat-transfer coefficients are much larger than the corresponding external heat-transfer coefficients. Accordingly, the total heat-transfer coefficients are mainly dependent on the external heat-transfer coefficients. In particular, with the evaporation of the naphtha feedstock, the internal heat-transfer coefficient increases dramatically, resulting in an increased difference between the internal and external heat-transfer coefficients. In the latter stage of the evaporation process in the HTC-I bank, the internal heat-transfer coefficient decreases. Thus, part of the heat resistance is gradually transferred to the inside tube, and the heat transfer inside the tube begins to impact the total heat transfer more. Moreover, the total heat-transfer coefficient profiles are similar to those for internal heat transfer in the two banks. In particular, as can be seen in Figure 14b, an abrupt drop in the total and internal heat-transfer coefficients is observed. Figure 15 presents the temperature profiles of the feed fluid, external tube skin, and flue gas. In Figure 15a, both the feedfluid and external tube-skin temperatures in the FPH tube bank increase gradually along the tube axis. Because more heat is used for evaporation, the increase in the feed-fluid temperature levels off in zone 175, and a reduced temperature difference between the feed-fluid temperature and the external tube-skin temperature is also observed. Figure 15b shows the feed-fluid and external tube temperature profiles in the HTC-I bank. In this section, some dilute steam is added to the feedstock to promote evaporation of the feedstock. Moreover, because of addition of some dilute steam, a lower partial pressure for each pure component is obtained which further promotes evaporation. Hence, more heat is needed for evaporation, which leads to an inlet temperature lower than the outlet temperature of the FPH bank. On the other hand, because of the higher flue-gas temperature, a higher external tube heat flux is obtained, resulting in a rapid rise in the feed-fluid temperature. As can be seen in Figure 15b, a slight increase in the already-rising feedfluid temperature is observed in zone 25, where the naphtha feedstock is fully evaporated. For the external tube-skin temperature, an abrupt rise is observed in zone 15 and can be explained by the partial dryout of the feedstock flow that occurs inside the tube, resulting in a rapid drop in the

Figure 10. Comparison of the flow-pattern maps predicted by the VLE and simplified models: (a) VLE model and (b) simplified model.

Figure 11. Heat-transfer coefficient profiles predicted by the VLE and simplified models.

vapor phase exhibits a significant increase in the mixture velocity, which leads to an increased Reynolds number, thereby enhancing the heat transfer inside the tube as well. Figure 14b shows the heat-transfer coefficient profiles in the HTC-I bank. As can be seen, the internal heat-transfer coefficient at the HTC-I bank inlet is higher than that predicted at the FPH bank outlet. This is because the higher temperature of the flue gas in the HTC-I bank leads to a higher heat flux, as shown in Figure 15. With the intensive evaporation and increased vapor quality, the mixture velocity and the disturbances in the 10821

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research

Figure 12. Vapor quality profiles along the tube axis in the (a) FPH and (b) HTC-I banks.

Figure 13. Flow-pattern maps in the (a) FPH and (b) HTC-I banks as predicted by Wojtan et al.’s method.20,21

Figure 14. Heat-transfer coefficient profiles along the tube axis in the (a) FPH and (b) HTC-I banks (IHTC, internal heat-transfer coefficient; ETHC, external heat-transfer coefficient; THTC, total heat-transfer coefficient).

fluid in the HTC-I tube bank leads to a larger heat flux than that in the FPH tube bank, as shown in Figure 16. Figure 16a depicts the external tube heat-flux profile along the tube axis in the FPH bank. The heat flux decreases along the tube axis in each pass and is a result of the external tubeskin temperature increasing along the axial direction, heated by the flue gas, which leads to a decline in the temperature difference between the external tube skin and the flue gas. Moreover, the heat flux in the first zone of each pass is higher than that in

heat-transfer coefficient. Once the naphtha feedstock has completely evaporated, the external tube-skin temperature falls abruptly in zone 25 due to a corresponding rise in the heattransfer coefficient. Because of the identical temperature of the flue gas in each layer, the flue-gas temperature has a stepped increase both in the FPH and HTC-I tube banks which leads to a stepped increase in the external heat-transfer coefficient. From Figure 15, the larger temperature difference between the flue gas and feed 10822

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research

Figure 15. Temperature profiles of the feed fluid, external tube skin, and flue gas along the tube axis in the (a) FPH and (b) HTC-I banks (t, feed-fluid temperature; tw, external-skin temperature; tg, flue-gas temperature).

Figure 16. External tube heat-flux profiles along the tube axis in the FPH and HTC-I banks.

two-phase vapor−liquid system during flow boiling, a flowpattern-based method was used in the simulation. The significant differences between the simulation results predicted by the simplified model and the VLE model indicate that the use of the VLE model to simulate the evaporation banks of the convection section is necessary. The simulation results are in good agreement with the design data. The profiles of vapor qualities, heat-transfer coefficients, feed-fluid temperatures, external tube-skin temperatures, flue-gas temperatures, and external tube heat fluxes are evaluated. Because of the high mass velocity and low vapor quality, slug and intermittent flows are predicted in the FPH bank, which enhance heat transfer inside the tube. In the HTC-I bank, the flow regime is mainly located in the stratified-wavy flow region because if the low mass flow. In the latter stage of evaporation, the occurrence of annular flow with partial dryout causes an abrupt decrease in heat transfer. The heat resistance is mainly concentrated on the flue-gas side, where the external heat-transfer coefficient is much smaller than the internal heat-transfer coefficient, with the exception that in the latter stage of evaporation and during single-phase vapor flow. This work has mainly focused on the evaporation of a naphtha feedstock in the evaporation banks. To obtain an even more accurate simulation of the flow boiling process in the convection section, the techniques presented in this work should

the last zone of the previous pass, which is caused by higher flue-gas temperature in the next pass. The entire curve shows a trend in which the external tube heat flux decreases gradually in the direction opposite to the flue-gas flow. This is caused by the decrease in the temperature difference between the external skin and the flue gas. However, because of an abruptly rising total heat-transfer coefficient and the extra heat required for evaporation, the external tube heat flux increases only slightly as evaporation begins. Figure 16b depicts the external tube heat-flux profile for the HTC-I bank. As was observed for the FPH bank, the heat flux decreases in each pass and falls sharply and abruptly in tube zone 20, where the vapor quality is about 0.8. This is because, when the vapor quality reaches a certain level, partial dryout of the tube occurs, leading to a decrease in the heat-transfer coefficient. When the naphtha feedstock is fully evaporated, the heat flux increases because of an increase in the total heat-transfer coefficient.

5. CONCLUSIONS The evaporation process of a complex naphtha feedstock in the convection section of a steam cracking furnace was simulated. An evaporation model based on a feedstock-simplifying method and vapor−liquid equilibrium theory was employed to simulate the evaporation process. To calculate the heat transfer of the 10823

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research

Q = total heat transferred from the flue gas in the tube zone, W ri = volume parameter of component i in the UNIFAC model ReD = Reynolds number of the flue gas ReH = homogeneous Reynolds number of mist flow ReG = Reynolds number of the flue-gas side Rev/Reσ = Reynolds number of the vapor/liquid phase T = temperature, K Tin/Tout = inlet/outlet temperature of the mixture in each zone, K Tg = temperature of the flue gas, K Umax = maximum velocity in the minimum flow area, m/s V = mass velocity of the vapor phase, kg m−2 s−1 v = molar volume, m3/mol x = vapor quality xde = dryout completion quality xdi = dryout inception quality xi = mass (mole) fraction of the ith species in the liquid phase xIA = vapor quality at transition from intermittent to annular flow Y = multiplication factor yi = mass (mole) fraction of the ith species in the vapor phase

be further extended to the coupled simulations of the entire convection section. This modeling will also provide useful results for the radiant section. Moreover, combined simulations of the convection and radiation sections will allow the cracking furnace to be optimized in an integrated fashion.



AUTHOR INFORMATION

Corresponding Authors

*Tel.: +86-13636576934. E-mail: [email protected] *Tel.: +86-21-64252060. E-mail: [email protected] ORCID

Feng Qian: 0000-0003-2781-332X Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS This work was supported by the National Natural Science Foundation of China (61333010, 61590923, 61422303). NOMENCLATURE a = attractive parameter of EoS A = external tube surface area, m2 Akl = group-contribution parameters in the PPR78 model b = molar covolume of EoS Bkl = group-contribution parameters in the PPR78 model C̅ p = mean heat capacity of the vapor/liquid phase, J kg−1 K−1 CPv,i,/CPl,i, = heat capacity of component i in the vapor/liquid phase, J kg−1 K−1 d = internal tube diameter, m dp = fin pitch, m D = external tube diameter, m ein/eout = inlet/outlet vapor quality in each zone F = mass velocity of the feedstock, kg m−2 s−1 Kij = binary-interaction parameter l = fin height, m L = mass velocity of the liquid phase, kg m−2 s−1 G = mass velocity of the liquid plus vapor, kg m−2 s−1 Gstrat = stratified flow transition mass velocity, kg m−2 s−1 Gwavy = wavy flow transition mass velocity, kg m−2 s−1 hcb = nucleate boiling heat-transfer coefficient, W m−2 K−1 hdryout = dryout heat-transfer coefficient, W m−2 K−1 hf = external heat-transfer coefficient, W m−2 K−1 hmist = mist-flow heat-transfer coefficient, W m−2 K−1 hnb = convection heat-transfer coefficient, W m−2 K−1 htotal = total heat-transfer coefficient, W m−2 K−1 htp = internal heat-transfer coefficient, W m−2 K−1 hv = vapor convective heat-transfer coefficient, W m−2 K−1 hw = tube wall thermal conductivity, W m−1 K−1 hwet = wet perimeter heat-transfer coefficient, W m−2 K−1 KG = thermal conductivity of the flue gas, W m−2 K−1 Kv/KL = thermal conductivity of the vapor/liquid phase, W m−2 K−1 M = molecular weight, kg/kmol Ng = number of different groups defined by the group contribution method P = pressure, Pa Pr = reduced pressure = P/Pc PrG = Prandtl number of the flue-gas side PrL/PrV = Prandtl number of the vapor/liquid phase q = internal heat flux, W/m2 qi = area parameter of component i in the UNIFAC model

Greek Symbols



αik = fraction of molecule i occupied by group k γi = liquid activity coefficient of component i Γk = activity coefficient of group k at the mixture composition, and Γik is the activity coefficient of group k at a group composition corresponding to pure component i Γik = activity coefficient of group k at a group composition corresponding to pure component i δ = liquid film thickness, m ΔHlatent = latent heat of the mixture, J/kg ε = cross-sectional vapor void fraction θdry = dry angle of the tube perimeter, rad θstrat = dry angle of the tube perimeter in stratified flow, rad λ = friction coefficient μ = dynamic viscosity, N s/m2 υki = number of group k in component i in the UNIFAC model ξ = symbol for property ρ = density, kg/m3 ρ̅ = density of the vapor−liquid mixture, kg/m3 ρv/ρL = vapor/liquid density, kg/m3 τ = interaction parameter between two groups in the UNIFAC model φD = diagonal pitch of the staggered arrangement φL = longitudinal pitch φT = transverse pitch φvi = vapor-phase fugacity of component i vLi /vRi = distilled volume or mass fraction of the left/right of interval i Ω = internal cross-sectional surface area, m2

REFERENCES

(1) Oprins, A. J. M.; Heynderickx, G. J. Calculation of threedimensional flow and pressure fields in cracking furnaces. Chem. Eng. Sci. 2003, 58, 4883−4893. (2) Oprins, A. J. M.; Heynderickx, G. J.; Marin, G. B. Threedimensional asymmetric flow and temperature fields in cracking furnaces. Ind. Eng. Chem. Res. 2001, 40, 5087−5094.

10824

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825

Article

Industrial & Engineering Chemistry Research (3) Heynderickx, G. J.; Oprins, A. J. M.; Marin, G. B.; Dick, E. Threedimensional flow patterns in cracking furnaces with long-flame burners. AIChE J. 2001, 47, 388−400. (4) 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, 103−113. (5) 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, 635−649. (6) 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. (7) Guihua, H.; Honggang, W.; Feng, Q. Numerical simulation on flow, combustion and heat transfer of ethylene cracking furnaces. Chem. Eng. Sci. 2011, 66, 1600−1611. (8) Hu, G.; Schietekat, C. M.; Zhang, Y.; Qian, F.; Heynderickx, G.; Van Geem, K. M.; Marin, G. B. Impact of radiation models in coupled simulations of steam cracking furnaces and reactors. Ind. Eng. Chem. Res. 2015, 54, 2453−2465. (9) Van Geem, K. M.; Ž ajdlík, R.; Reyniers, M.-F.; Marin, G. B. Dimensional analysis for scaling up and down steam cracking coils. Chem. Eng. J. 2007, 134, 3−10. (10) De Schepper, S. C. K.; Heynderickx, G. J.; Marin, G. B. CFD modeling of all gas−liquid and vapor−liquid flow regimes predicted by the Baker chart. Chem. Eng. J. 2008, 138, 349−357. (11) De Schepper, S. C. K.; Heynderickx, G. J.; Marin, G. B. Modeling the evaporation of a hydrocarbon feedstock in the convection section of a steam cracker. Comput. Chem. Eng. 2009, 33, 122−132. (12) De Schepper, S. C. K.; Heynderickx, G. J.; Marin, G. B. Coupled simulation of the flue gas and process gas side of a steam cracker convection section. AIChE J. 2009, 55, 2773−2787. (13) Mahulkar, A. V.; Heynderickx, G. J.; Marin, G. B. Simulation of coking in convection section of steam cracker. Chem. Eng. Trans. 2012, 29, 1375−1380. (14) Mahulkar, A. V.; Heynderickx, G. J.; Marin, G. B. Simulation of the coking phenomenon in the superheater of a steam cracker. Chem. Eng. Sci. 2014, 110, 31−43. (15) De Schepper, S. C. K.; Heynderickx, G. J.; Marin, G. B. Modeling the coke Formation in the convection section tubes of a steam cracker. Ind. Eng. Chem. Res. 2010, 49, 5752−5764. (16) Hu, G.; Yuan, B.; Zhang, L.; Li, J.; Du, W.; Qian, F. Coupled simulation of convection section with dual stage steam feed mixing of an industrial ethylene cracking furnace. Chem. Eng. J. 2016, 286, 436− 446. (17) Chen, J. C. Correlation for boiling heat transfer to saturated fluids in convective flow. Ind. Eng. Chem. Process Des. Dev. 1966, 5, 322−329. (18) Gungor, K. E.; Winterton, R. H. S. A general correlation for flow boiling in tubes and annuli. Int. J. Heat Mass Transfer 1986, 29, 351− 358. (19) Liu, Z.; Winterton, R. H. S. A general correlation for saturated and subcooled flow boiling in tubes and annuli, based on a nucleate pool boiling equation. Int. J. Heat Mass Transfer 1991, 34, 2759−2766. (20) Wojtan, L.; Ursenbacher, T.; Thome, J. R. Investigation of flow boiling in horizontal tubes: Part IA new diabatic two-phase flow pattern map. Int. J. Heat Mass Transfer 2005, 48, 2955−2969. (21) Wojtan, L.; Ursenbacher, T.; Thome, J. R. Investigation of flow boiling in horizontal tubes: Part IIDevelopment of a new heat transfer model for stratified-wavy, dryout and mist flow regimes. Int. J. Heat Mass Transfer 2005, 48, 2970−2985. (22) Verhees, P.; Amghizar, I.; Goemare, J.; Akhras, A. R.; Marin, G. B.; Van Geem, K. M.; Heynderickx, G. J. 1D Model for Coupled Simulation of Steam Cracker Convection Section with Improved Evaporation Model. Chem. Ing. Tech. 2016, 88, 1650−1664. (23) Edmister, W. C. Improved integral technique for petroleum distillation calculations. Ind. Eng. Chem. 1955, 47, 1685−1690.

(24) Hariu, O. H.; Sage, R. C. Crude split figured by computer. Hydrocarbon Process 1969, 48, 143−148. (25) Eckert, E. Do we need pseudocomponents? Chem. Listy 2001, 95, 368−373 (in Czech). (26) Steyer, F.; Sundmacher, K. VLE and LLE Data for the System Cyclohexane + Cyclohexene + Water + Cyclohexanol. J. Chem. Eng. Data 2004, 49, 1675−1681. (27) Hoffmaster, W. R.; Hauan, S. Using feasible regions to design and optimize reactive distillation columns with ideal VLE. AIChE J. 2006, 52, 1744−1753. (28) Modla, G.; Lang, P. Separation of an Acetone−Methanol Mixture by Pressure-Swing Batch Distillation in a Double-Column System with and without Thermal Integration. Ind. Eng. Chem. Res. 2010, 49, 3785−3793. (29) Eckert, E.; Vaněk, T. New approach to the characterisation of petroleum mixtures used in the modelling of separation processes. Comput. Chem. Eng. 2005, 30, 343−356. (30) Khan, W. A.; Culham, J. R.; Yovanovich, M. M. Convection heat transfer from tube banks in crossflow: Analytical approach. Int. J. Heat Mass Transfer 2006, 49, 4831−4838. (31) Qian, J. L. Tubular Heating Furnace; Petrochemical Press: Beijing, 2003. (32) De Saegher, J. J. Modellering van Stroming, Warmtetransport en Reactie in Reactoren voor Thermisch Kraken van Koolwaterstoffen. Ph.D. Thesis, Ghent University, Ghent, Belgium, 1994. (33) 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, 1020−1034. (34) 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, 3174−3188. (35) Plazas Tovar, L.; Wolf Maciel, M. R.; Maciel Filho, R.; Batistella, C. B.; Celis Ariza, O. J.; Medina, L. C. Overview and Computational Approach for Studying the Physicochemical Characterization of HighBoiling-Point Petroleum Fractions (350°C+). Oil Gas Sci. Technol. 2012, 67, 451−477. (36) Jaubert, J. N.; Privat, R.; Mutelet, F. Predicting the phase equilibria of synthetic petroleum fluids with the PPR78 approach. AIChE J. 2010, 56, 3225−3235. (37) Peng, D.-Y.; Robinson, D. B. A new two-constant equation of state. Ind. Eng. Chem. Fundam. 1976, 15, 59−64. (38) Xu, X.; Jaubert, J.-N.; Privat, R.; Duchet-Suchaux, P.; BranaMulero, F. Predicting binary-interaction parameters of cubic equations of state for petroleum fluids containing pseudo-components. Ind. Eng. Chem. Res. 2015, 54, 2816−2824. (39) Asgarpour Khansary, M.; Hallaji Sani, A. Using genetic algorithm (GA) and particle swarm optimization (PSO) methods for determination of interaction parameters in multicomponent systems of liquid−liquid equilibria. Fluid Phase Equilib. 2014, 365, 141−145. (40) Larsen, B. L.; Rasmussen, P.; Fredenslund, A. A modified UNIFAC group-contribution model for prediction of phase equilibria and heats of mixing. Ind. Eng. Chem. Res. 1987, 26, 2274−2286.

10825

DOI: 10.1021/acs.iecr.7b02806 Ind. Eng. Chem. Res. 2017, 56, 10813−10825