Simulation of the Evaporation Tube Banks in the Convection Section

Subscriber access provided by Warwick University Library

Article

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 Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/acs.iecr.7b02806 • Publication Date (Web): 05 Sep 2017 Downloaded from http://pubs.acs.org on September 10, 2017

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

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

Page 1 of 45

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

Industrial & Engineering Chemistry Research

For Table of Contents Only 547x220mm (96 x 96 DPI)

ACS Paragon Plus Environment


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

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*, 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 *.Corresponding Author: Tel: +86-13636576934, E-mail address: [email protected] **.Corresponding Author: Tel: +86-21-64252060, E-mail address: [email protected]

Abstract Simulation of the evaporation process of a naphtha feedstock in the convection section of a steam cracking furnace is reported in this paper. A two-step evaporation model was developed to determine component changes of the vapor and liquid phases. Firstly, a feedstock-simplifying method, based on Eckert’s work, was developed in which the original hydrocarbon mixture is represented by a system of real components. Secondly, component distribution was determined by vapor-liquid equilibrium theory, using an activity coefficient method. The group contribution method (GCM) was 1


Page 2 of 45

Page 3 of 45

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


adopted to calculate the interaction parameters in the activity coefficient method. Using a mechanistic heat transfer method, a comparison between the evaporation model developed in this paper and a simplified evaporation model was made. Simulation results are in good agreement with industrial data. Some parameter profiles in the heat transfer process are provided, as are the flow patterns of the naphtha feedstock. Keywords: Convection section, Evaporation model, Real component, Vapor-liquid equilibrium, Flow pattern map.

1. Introduction Steam cracking furnace is the key unit to produce ethylene in chemical industry. It mainly contains two parts: radiation section and 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 bundle, preheating and vaporizing the hydrocarbon feedstock inside the tube 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 have usually focused on the radiation section.1-9 With an increasing demand for ethylene and propylene, and the decreasing supply of light feedstocks, inexpensive but heavier feedstocks, such as heavy naphtha and heavy vacuum gas oil (HVGO) are being cracked more often. Apparently, this change has the potential to cause more fouling problems within the tubes in the convection section. In addition, 2



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

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. In order to analyze and solve these problems, an accurate description of the evaporation process in the convection section is necessary. Over the past few decades, significant progress in the modeling of the convection sections of steam crackers has been made, mostly based on the computational fluid dynamics (CFD) method. De Schepper10-12 conducted a CFD study on the feedstock temperature field inside tubes, and the flue-gas flow velocity and temperature fields in the convection chamber. Mahulkar13-14 and De Schepper15 analyzed the coking mechanism and simulated the coke-forming process of heavy feedstocks in the convection section. Some suggestions on how to reduce coke formation were also provided based on these studies. Hu16 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 that include: (1) simulation by the CFD technique is extremely time-consuming and, as a consequence, has limited application in optimization and control of steam cracking furnace; 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, such as thermodynamic ( , , and ) and transport properties (kinematic viscosity (ν), thermal conductivity (λ), specific heat capacity ( ) and vapor pressure ( )), which are very important for 3


Page 4 of 45

Page 5 of 45

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


heat transfer calculations in flow boiling processes. Thus, a fast and accurate evaporation model of the convection section is necessary for broad application in chemical industries. Many researchers have focused on the heat transfer during flow boiling.17-19 Among previous studies, a mechanistic model was developed by Wojtan.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 vapor quality. To calculate accurate physical properties, an evaporation model that considers component changes is necessary. Verhees22 developed a 1-D model for the coupled simulation of a steam cracker convection section based on Wojtan’s correlations. In his model, a set of pseudo-components was adopted for the characterization of a multicomponent hydrocarbon feedstock. In addition, an evaporation model, based on the assumption that all pseudo-components evaporate one by one, was used. The evaporation order of these pseudo-components 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 cannot be ignored. Firstly, the widely accepted pseudo-component method23-24 relies on pseudo-components which are primarily defined by their (pseudo) boiling points and some additional parameters, mostly density, molecular weight or viscosity. However, for the determination of other thermodynamic and transport properties, some predictive methods are needed. In 4



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

Page 6 of 45

particular, the group contribution method, which requires detailed molecular-structure information to estimate interaction parameters, cannot be used for calculation of vapor-liquid equilibrium in this context. In 2001, Eckert made a comparison of molecular weight, critical temperature, critical pressure and the acentric factor for 98 pure components and studied the reliability of the estimation methods of physical properties in the pseudo-component method.25 The results showed that the estimation methods mostly led to larger errors. Secondly, due to interaction of the molecules, the species in the actual mixture cannot evaporate one by one, as described in Verhees’s model. Thus, the evaporation mechanism in Verhees’s work, which is referred to as the "simplified model" in this paper, causes certain errors in the simulation results. Vapor-liquid equilibrium theory (VLE) 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

feedstock-simplifying model is developed based on Eckert’s work.29 In the second step, vapor-liquid equilibrium theory, with some predicted interaction parameters, is used to calculate component distribution during evaporation. On the basis of coupled simulations between the tube section and the convection chamber, the flow-pattern oriented heat transfer method proposed by Wojtan was adopted to calculate heat 5


Page 7 of 45

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


transfer rate and predict flow patterns during flow boiling. Consequently, comparative testing of the evaporation mechanism using Verhees’s model and that described in this paper was undertaken. Vapor qualities, heat transfer 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.

2. Mathematical model 2.1 Heat balance The heat transferred to a fluid in a tube section of ∆ is calculated by:

Q = htotal P∆L(Tg − T )

(1)

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

 1 d d 1 d ln + = +  h  tp 2hw D h f D 

−1

(2)

For the single vapor/liquid phase flow, the increase of fluid temperature over the tube section ∆ is caused by single vapor/liquid phase forced convection: Q = FC p (Tout − Tin )

(3)

When evaporation occurs, the vapor phase forms inside the tube. The compositions in the two phases vary during evaporation. The heat conservation equation Eq. (4) was used to calculate the outlet temperature during flow boiling in a section of tube, as follows:

6



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

（ Fein ∑ yi C Pv ,i + F (1 − ein ) ∑ xi C P l ,i )(Tout − Tin )=Q + F (eout − ein ) ∆H latent (4) i

i

2.2 Heat transfer model 2.2.1 Internal heat transfer

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

Figure1 depicts several typical flow regimes encountered during the flow boiling process in a horizontal tube, in which distinct heat transfer mechanisms are needed in order to address different flow patterns. The heat transfer correlations developed by Wojtan et al. are shown in Table 1. In Wojtan’s work, 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 heat transfer coefficients, resulting in Eqs. (5)-(9). Under various flow regimes, the dry angle, θ, is defined in different forms as follows:

θ= 0, for slug flow, intermittent flow and annular flow;

θ= θ , for stratified flow;

θ =

θ = *

*

!

!"#$#

+,

%

&.()

!

!"#$#

%

θ , for stratified-wavy flow;

&.()

θ , for slug/stratified-wavy flow. 7


Page 8 of 45

Page 9 of 45

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


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 only transferred 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, the heat transfer coefficient decreases dramatically and then levels out at the value for mist flow (M). Hence, the heat transfer coefficients can be calculated from Eq. (12), by linear interpolation, where ℎ -./0 1 is the heat transfer coefficient of the vapor quality at the dryout inception point, ./0 , calculated by Eq.(5), and ℎ203 -./4 1 is the heat transfer coefficient of the vapor quality at the mist inception point, ./4 , calculated by Eq.(12). Table 1. Tube interior heat-transfer correlation equations. Model Wojtan et al.

Equation S, SW, Slug, Slug + SW, I, and A flow regimes:

h tp =

θ dry hv + (2π − θ dry ) hwet 2π

0.4 hv = 0.023Re0.8 v Prv

(5)

kv D

(6)

hwet = [hnb3 + hcb3 ]1/3

(7)

hcb = 0.0133Reδ0.69 PrL0.4

kL

δ

hnb = 55(Pr )0.12 (− log Pr )−0.55 M −0.5q0.67

8


(8) (9)


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

1/3    3π  1/3 1/3 π(1−ε) +  (1−2(1−ε) +(1−ε) −ε )  2 θstrat = 2π −2  − 1−ε ε 1− 2(1−ε) 1+ 4 1−ε 2 +ε 2 / 200 ( ) ( ) ( )  

))

( (

  x 1− x   (1 + 0.12 (1 − x ) )  +   ρL    ρv x  ε=  ρv 1.18 (1 − x ) ( gσ ( ρ L − ρ v ) )0.25  +  0.5   G ρ   L

Page 10 of 45

(10)

−1

(11)

Mist flow regime: 1.06 −1.83 hmist = 0.0117Re0.79 H PrV Y

kv D

(12)

Annular flow with partial dryout:

hdryout = htp ( xdi ) −

x − xdi [htp ( xdi ) − hmist ( xde )] xde − xdi

(13)

2.2.2 External heat transfer For coupled simulation of evaporation tube banks in the convection section, the heat transfer in the convection chamber also needs to be considered. Since the evaporation tube banks are located on the top of the convection section, where the flue gas temperature is low about 373K to 573K, 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 Eq. (14) 30:  ϕ  ϕ U max = max  T U app , T U app  ϕD −1  ϕT − 1 

(14)

where ϕ D = ϕ L 2 + (ϕ L /2 ) 2 is the dimensionless diagonal pitch in the case of the 9


Page 11 of 45

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


staggered arrangement. The Reynolds number for any arrangement is defined by: Re D =

ρ dU max µ

(15)

Since finned annulus tubes are used in the evaporation banks, a correlation for calculating the heat transfer in the finned tube banks is provided by Eq. (16) 31:

d  k h f =0.1378 G ReG0.718 PrG1/8  p  Db  l 

0.296

(16)

2.3 Evaporation model 2.3.1 Characterization of feedstock using real components The characterization of the heavy feedstock using real components is similar to that of the pseudo-component method. Unlike pseudo-components, 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, i.e. a set of substitute components is firstly selected and then the composition of this mixture is determined. The substitute components are firstly selected using characterization curves, in particular true boiling point (TBP) curves. As only ASTM D86 distillation curves are usually available in industry, the TBP curve of the naphtha feedstock in this work was calculated by transforming the ASTM D86 curve using an empirical correlation.32 Like the pseudo-component method, the boiling point range of the TBP curve must be cut in order to obtain non-overlapping temperature intervals (0 , 06) , for 7 = 1... 8) that cover the entire temperature range. For each temperature interval, a set of 10



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

candidate real components is selected from a database by their boiling point within the temperature interval. Hence, a sufficiently large naphtha molecular database, which contains all the components of a naphtha feedstock as much as possible, is required for accurate calculations. However, the detailed composition of a naphtha feedstock cannot be obtained straightforwardly. Thus in this work, the SIMCO software is adopted to predict the detailed composition of a naphtha feedstock by application of the molecular reconstruction method.33-34 It has a power molecular library which is obtained from many real naphthas by experimental methods. Some industrial indices of the naphtha feedstock are taken as input of this software, as listed in Table 2. Table 2. Industrial indices of the naphtha feedstock (IBP, initial boiling point, NBP, final boiling point). Specific density (20 ℃) (g/cm3) Specific density (15 °C) (g/cm3) ASTM D86 Boiling Points (K) IBP 10% 50% 90% 95% FBP PIONA（ wt%） n-Paraffins i-Paraffins Olefins Napthenes 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

After the candidate components are 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 11


Page 12 of 45

Page 13 of 45

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


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:

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

K

∑ω

k

k =1

ξ m , k ,c

→ min c

(17)

In Eckert’s work,29 9: is the weighting factor that reflects the precision of the measurements. ;2,:, and ;,:, 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, with the specific density curve used occasionally. Unfortunately, other physical properties, such as thermodynamic and transport properties, which are crucial to the heat transfer calculation, cannot be easily measured. Hence, simply comparing the common properties is not enough and the small relative differences of the common properties calculated by Eq. (17) may not represent the small differences of other properties. To address this issue, another similar optimization task is established that can be defined by the objective function: l

ξ m , l − ξ r ,l

l =1

ξ m ,l

∑

→ min

(18)

In this equation, ;2, and ;, are the measured thermodynamic and transport properties and those of the selected components, respectively. Before calculating Eq. (18), the representative component for each temperature interval should be selected firstly. Thus, Eq. (17) is still needed. However, if the optimization procedure needs to 12



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

be implemented for each characterization, it is complicated and time-consuming. If the tendency that the small differences of the common properties or part of these properties between the temperature interval and the selected components correspond to the small differences of the thermodynamic and transport properties can be proved, then only the common properties need to be considered and the calculation process can be greatly simplified. To prove this ideal, 9: in Eq. (17) is redefined as the weighting factor and reflects the influence of the common properties on thermodynamic and transport properties. Here, 9: is in the range of 0 – 1 and ∑ 9: = 1. In this work, only the TBP and specific density curves were obtained. Hence, only two weighting factors, 9 >?@ and 93/ , need to be determined. For the calculation accuracy, the best choice is using the thermodynamic and transport properties measured by the experiment to compare with those of the selected components. However, as it is difficult to measure the thermodynamic and transport properties, some predictive methods for the basic properties of the pseudo-components were adopted.35 Although this is a compromise solution, compared with the pseudo-component method, the thermodynamic and transport properties of the real component are easily obtained which further simplifies the calculation process and gets similar results. The results show that when 9>?@ exceeds 0.9 and 93/ is within 0.1, the thermodynamic and transport properties of the selected real component are the closest to the predicted values. In other words, the normal boiling point of the selected component is closer to the mean temperature of the temperature interval, the thermodynamic and transport 13


Page 14 of 45

Page 15 of 45

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


properties of the selected component are closer to those of the pseudo component which is represented for the temperature interval. Thus, selecting the real component by simply comparing the boiling point is precise enough for these calculations, and the influence of 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, while some thermodynamic properties (critical pressure, , and critical temperature, ) of the selected real components are depicted in Figure 3.

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

14



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

Page 16 of 45

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

After the representative components are determined, the composition distribution of these components needs to be determined. This paper adopted the optimization method proposed by Eckert.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 are generally unable to 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 follows: I +1

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

ν < ν iL+1 , i = 1,..., I L i

(19)

ν i −1 < ν iL < ν i +1 , i = 2,..., I In this equation, A0B and A0C denote the starting and end distillation fractions. To determine the relative proportion of the selected component, 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 shown in Table 3.

15


Page 17 of 45

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


Figure 4. Construction of distilled-fraction intervals with real components. Table 3. The composition of the selected components by optimization The selected component

Mass Fraction (%)

Neopentane

4.13

Isopentane

1.22

Cyclopentane

4.72

n-hexane

18.48

2,3-Di-methyl heptane

14.74

2,4-Di-methyl octane

8.14

2,2,4-Tri-methyl nonane

16.82

n-Nonane

7.17

3-methyl decane

12.91

n-Undecane

10.82

n-Dodecane

0.85

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, due to the molecular interaction, the 16



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

Page 18 of 45

real system cannot be taken as an ideal system. Thus, two parameters, the vapor phase fugacity (D0 ) and liquid activity (E0 1 coefficients, are introduced. A well-known vapor-liquid equilibrium equation with the two parameters was adopted:

ϕiv yi P = xiγ i Pi s

(20)

In this equation, D0 denotes the pressure correction of the true gas to the ideal gas, E0 denotes the mole fraction correction of the true solution to the ideal solution. The vapor-phase fugacity coefficient is usually computed from the equations of state (EoS). However, for the mixture, accurately quantifying the interactions between each pair of molecules is difficult. To avoid such a fastidious work, an alternative solution is adopted using a predictive method, which can estimate the interactions from mere knowledge of the structure of molecules. Thus, PPR78 model36 is adopted, which combines the Peng–Robinson(PR) EoS37 and a group contribution method (GCM).36, 38

This approach can be seen as a GCM to estimate 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), D0 , which is temperature dependent, can be predicted. Table 4. Equations of the PPR78 model

17


Page 19 of 45

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


PR78 model ln ϕ iv =

bi A  2 ( Z m − 1) − ln( Z m − B ) −  bm 2 2 B  am

N

∑x

j

j =1

ai a j −

bi bm

  Z m + ( 2 + 1) B   ln     Z m − ( 2 − 1) B 

 R = 8.314472 Jmol −1 K − 1   b = 0.0777960739 RTc ,i  i Pc ,i  2  R 2Tc2,i  T  = + − a 0.457235529 1 m (1 )    i i Pc ,i  Tc ,i    2  if ω i ≤ 0.491 mi = 0.37464 + 1.5 422ω i − 0.26992ω i    if ω i > 0.491 mi = 0.379642 + 1.48503ω i − 0.164423ω i2 + 0.016666ω i3

(21)

Mixing rule N N   a m = ∑ ∑ yi y j i =1 j =1  N   b m = ∑ y i bi  i =1   A = am P  m (RT )2   B = P bm  m RT

a i a j 1 − K ij (T ) 

(22)

GCM

Kij (T ) =

Bkl N N a j (T )  1 g g 298.15 ( Akl −1)   ai (T )  − ∑∑(αik − α jk )(αil − α jl ) Akl × ( ) −  −  2  k =1 l =1 T /K bj    bi 

2

(23)

ai (T )a j (T ) bb i j

In Eq. (20), the activity coefficient of the liquid phase is also required. Several models have been proposed to calculate the activity coefficient, such as the Wilson, universal quasi chemical (UNIQUAC), and non-random two liquid (NRTL) models.39 As was the case for the calculation of BIP, some interaction parameters in these methods need to be fitted using vapor-liquid equilibrium data. Unfortunately, the experimental data are hard to obtain. This paper adopts the UNIFAC model to estimate activity coefficients.40 In the UNIFAC model, the activity coefficient is 18



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

Page 20 of 45

calculated as the sum of combinatorial and residual contributions, described as follows:

ln γ i = ln γ ic + ln γ iR

(24)

where ln γ ic = ln

ϕi xi

+1−

ϕi xi

−

Z 2

 ϕi ϕi   ln + 1 −  θi   θi

(25)

and ln γ ic = ∑ vki ln Γ k − ln Γik 

(26)

k

The molecular volume fraction (D0 ) and the surface area fraction (F0 ) are defined in terms of the segment numbers G and constant parameter numbers H as:

ϕi =

xi ri xq and θi = i i ∑ x j rj ∑ xjqj j

(27)

j

In Eq. (26), Γ: is the activity coefficient of group J at the mixture composition, and

Γ:0 is the activity coefficient of group J at a group composition corresponding to pure component 7 . Γ: and Γ:0 are given by:

 z  θτ ln Γ k = Qk  1 − ln(∑ θ mτ mk ) − ∑ i ki 2  m i ∑ θ jτ ji j 

    

(28)

To calculate the detailed compositions of the two phases, a combination of the phase-equilibrium equation and the mass-balance equation is necessary. The mass-balance equation is given in Eq. (29) and the calculation process is implemented in AspenPlus.

19


Page 21 of 45

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


Fzi = Vyi + Lxi

(29)

3. Simulations 3.1 Tube geometry and operating conditions

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

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-Ι) bank. In addition to the FPH and HTC-I tube banks, the economizer tube bank (ECO) is also taken into account for a continuous temperature 20



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

Page 22 of 45

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 adopted 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 indicate that the temperature distribution of the flue gas tends to be uniform throughout the width and the length directions at higher elevation, because the staggered arrangement of the tube bundles enhances the mixing of the flue gas. In this work, it is reasonable to assume that the flue gas temperature in each tube layer is identical. Hence, only one tube is selected to represent all the tubes in each layer. Table 5. Tube geometry and operating conditions. Convection chamber specifications Length (z-direction) (m) 10.368

Width (x-direction) (m)

1.22

FPH tube bank Number of tube passes Number of tubes Tube outer diameter (m) Thickness of tube (m) Transverse pitch (m) Longitudinal pitch (m)

20 6 0.0603 0.00554 0.102 0.08833

Fin thickness (m) Fin number Fin height (m)

0.00127 157 0.01905

Mass flow rate of naphtha (kg/s) Inlet temperature (K)

8.486 333.15

ECO tube bank Number of tube passes Number of tubes Tube outer diameter (m) Thickness of tube (m) Transverse pitch (m) Longitudinal pitch (m)

18 6 0.0603 0.00874 0.102 0.102

Fin thickness (m) Fin number Fin height (m)

0.00127 157 0.01905

Mass flow rate of boiler feed water (kg/s) Tube inlet temperature (K)

10.48 418.15

HTC-I tube bank Number of tube passes Number of tubes Tube outer diameter (m) Transverse pitch (m)

6 6 0.1143 0.201

Thickness of tube (m) Fin number Fin height (m)

0.00602

21


0.01905

Page 23 of 45

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


Longitudinal pitch (m) Fin thickness (m) Flue gas side Mass flow rate (kg/s) Flue gas composition (mol %) O2 N2 H2O

0.203 0.00127

Mass flow rate of dilute steam (kg/s) Inlet temperature of dilute steam (K)

0.85 458.15

25.14

Inlet temperature of the flue gas (K)

837.15

1.98 71.07 12.26

CO2 Ar

13.47 1.22

3.2 Solution strategy In order to explain the calculation algorithm of the simulation for each tube bank, this paper will take the HTC-Ι bank as an example. As is shown in Figure 6, each tube pass is divided into ten 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 heat transfer calculation. At the commencement of the simulation, the flue gas temperatures between each of the two adjacent layers are firstly 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 by Eq. (1).

22



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

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

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 by heat balance of the flue gas. This iteration continues until the maximum differences of 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.

23


Page 24 of 45

Page 25 of 45

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


start

Construct geometrical model of Construction of geometrical the of tubes model tubes Construct geometrical model of convention chamber

Input the initial feedstock conditions

Initialize the flue gas temperature of each layer

Solve equations for heat and mass transfer

Output the absorbed heat of each tube pass

N

Input the initial flue gas conditions

Solve equations for heat and mass transfer

Output new flue gas temperature profiles

Convergence? Y End

Figure 7. Flow chart of the algorithm used for the coupled simulation of the tubes and the flue gas side.

4. Results and discussion 4.1 Comparison of the VLE and simplified evaporation mechanisms To analyze the influence of the VLE model and the simplified model described by Verhees 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 24



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

temperature of 333.15 K and 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．．Vapor quality profiles predicted by the VLE and simplified models.

25


Page 26 of 45

Page 27 of 45

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


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

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 the difference in the mixture velocity profiles. As is shown in Figure 9, due to an earlier evaporation, the mixture velocity predicted by the simplified model is obviously larger than that predicted by the VLE model at the first half of the tube. Because the vapor phase has a larger velocity than the liquid phase which promotes the mixture velocity inside the tubes. When at the latter half of the tube, the mixture velocities predicted by the two model are comparable, this is because the velocity of the vapor phase plays a 26



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

dominant role in the mixture velocity. In general, the larger mixture velocity predicted by the simplified model will lead to a larger pressure drop than that predicted by the VLE model. Detailed flow pattern maps are presented in Figure 10. The vertical lines (X = XIA) separating intermittent and annular flows at vapor qualities of 0.5 and 0.98 are given, the horizontal line labeled "GMNO " represents a mass flux of 495 kg/m2s. It can be observed that when the feedstock is determined, the flow pattern map is also determined. Thus, the flow patterns in the evaporation process mainly depend on the mass velocity. Moreover, due to variations in the physical properties of the vapor and liquid phases in the evaporation process, the vertical line that separates intermittent and annular flows 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 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 larger heat transfer coefficients near the tube entrance than that 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 it tends to be constant. Since the simplified model predicts a longer mist flow region than that predicted by the VLE model, lower heat transfer coefficients in the second half of the tube are also predicted. These significant differences in heat transfer coefficient profiles cannot be ignored, especially for the 27


Page 28 of 45

Page 29 of 45

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


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.

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

28



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

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

4.2 Application of the evaporation model to the evaporation bank Table 6. Comparison of industrial data and simulation results. Item

Design data

Simulation data

Fluid (FPH)/K

430

425.78

Fluid (ECO)/K

523.15

530.17

Fluid (HTC-I)/K

609

598.37

Flue gas (FPH)/K

385.15

390.02

Flue gas (ECO)/K

463.15

460.89

Flue gas (HTC-I)/K

632.15

639.45

Utilizing the evaporation model developed in this work, the evaporation process of the naphtha feedstock in the evaporation tube banks was simulated. Table 6 compares the simulation results with the design data. The simulated outlet temperatures of the 29


Page 30 of 45

Page 31 of 45

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


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’s heat transfer correlations.

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

Figure 13. Flow pattern maps predicted by Wojtan’s method in the FPH and HTC-I banks.

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 that of the naphtha feedstock, is added to the naphtha feedstock in the HTC-I bank 30



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

to promote the evaporation of the naphtha feedstock. Thus, the vapor quality at the inlet of the HTC-I bank is larger than that at the outlet of the FPH-bank. Different evaporation flow patterns occur due to variations in vapor quality. Figure 13 (a) shows the flow pattern map of the FPH bank predicted at G650 (mass velocity of 650 kgm-2s-1), which represents the mass velocity in the FPH bank. Such a high mass velocity and low vapor quality result in flow patterns which are located in the slug and intermittent flow regions. The flow pattern map of the HTC-Ι bank is illustrated in Figure 13 (b). As some dilute steam is added to the feed fluid, the vapor quality increases and reaches 0.37. Moreover, due to a larger tube diameter, the mass velocity in this bank is only 189 kgm-2s-1 which is much smaller than that of the FPH bank. Stratified-wavy flow, annular flow with partial dryout, and mist flow are predicted. Due to the low mass velocity, the flow pattern at the beginning of the tube is strait-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 14 (b). When the vapor quality exceeds Xde (the mist-flow inception point), mist flow occurs resulting in a constant heat transfer coefficient.

31


Page 32 of 45

Page 33 of 45

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


Figure 14. Heat transfer coefficient profiles along the tube axis in the FPH and HTC-I banks (IHTC, internal heat transfer coefficient; ETHC, external heat transfer coefficient; THTC, total 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, the boiling heat transfer leads to a significant increase in the heat transfer inside the tube. Moreover, the vapor phase has a significant increase in the mixture velocity which leads to an increased Reynolds number enhancing the heat transfer inside the tube as well. Figure 14 (b) shows the heat transfer coefficient profiles in the HTC-Ι bank, 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 higher heat flux, as shown in Figure 15. With the intensive evaporation and increased vapor quality, the mixture velocity and the disturbances in the boundary layer are enhanced which result in a continuous increase in the internal heat transfer coefficient. However, the internal heat transfer coefficient is observed to drop abruptly 32



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

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 13 (b), and is finally transformed into mist flow. In addition, there is also an abrupt rise in zone 25 where the naphtha feedstock becomes totally evaporated and the flow regime is transformed from mist flow into single-phase vapor flow. As is 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 naphtha feedstock evaporation, the internal heat transfer coefficient increases dramatically resulting in an increased difference between the internal and external heat transfer coefficients. As for the latter stage of the evaporation process in the HTC-I bank, the internal heat transfer coefficient decreases. Thus, part of heat resistance is gradually transferred to the inside tube and the heat transfer inside the tube begins to impact on 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 seen in Figure 14 (b), an abrupt drop in the total and internal heat transfer coefficients is observed.

33


Page 34 of 45

Page 35 of 45

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


Figure 15. Temperature profiles of the feed fluid, external tube skin and flue gas along the

tube axis in the FPH and HTC-I banks (t, feed fluid temperature; tw, external skin temperature; tg, flue gas temperature).

Figure 15 depicts the temperature profiles of the feed fluid, external tube skin and flue gas. In Figure 15(a), both the feed fluid and external tube skin temperatures in the FPH tube bank increase gradually along the tube axis. Due to more heat used for evaporation, the increase of 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 15 (b) shows the feed fluid and external tube temperature profiles in the HTC-Ι bank. In this section, some dilute steam is added to the feedstock to promote evaporation of the feedstock. Moreover, due to 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, due to 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 34



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

Figure 15 (b), a slight increase in the already rising feed fluid 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 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 heat transfer coefficient. Due to 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 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 16. External tube heat flux profiles along the tube axis in the FPH and HTC-I banks.

35


Page 36 of 45

Page 37 of 45

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


Figure 16 (a) 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 tube skin 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 at 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 opposite direction to the flue gas flow. This is caused by the decline in the temperature difference between the external skin and the flue gas. However, due to 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 16 (b) depicts the external tube heat flux profile for the HTC-Ι 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 due to 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 36



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

feedstock-simplifying method and vapor-liquid equilibrium theory was adopted to simulate the evaporation process. To calculate the heat transfer of the vapor-liquid two phases during flow boiling, a flow-pattern oriented 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. Due to the high mass velocity and low vapor quality, the 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 strait-wavy flow region due to 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 the naphtha feedstock in the evaporation banks. In order to obtain an even more accurate simulation of the flow boiling process in the convection section, the techniques presented in this work should be further extended to the coupled simulations of the entire convection section. This 37


Page 38 of 45

Page 39 of 45

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


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.

38



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

Nomenclature a

A

R: b

T:

attractive parameter of EoS external tube surface area, m2 group-contribution parameters in the PPR78 model molar covolume of EoS group-contribution parameters in the PPR78 model

UUU

mean heat capacity of the vapor/liquid phase, J/kg/K

CWX,Y, /CW[,Y,

heat capacity of component 7 in the vapor/liquid phase, J/kg/K

\

internal tube diameter, m

d^

fin pitch, m

_

external tube diameter, m

eYa /ebc d

e0f

inlet/outlet vapor quality in each zone mass velocity of the feedstock, kg/m2s binary-interaction parameter

g

fin height, m

h

mass velocity of the liquid plus vapor, kg/m2s

G

stratified flow transition mass velocity, kg/m2s

GiX

wavy flow transition mass velocity, kg/m2s

hkl

nucleate boiling heat transfer coefficient, W/m2/K

hbc

dryout heat transfer coefficient, W/m2/K

hm

external heat transfer coefficient, W/m2/K

hal

convection heat transfer coefficient, W/m2/K

hnY

mass velocity of the liquid phase, kg/m2s

mist flow heat transfer coefficient, W/m2/K

hb[

total heat transfer coefficient, W/m2/K

h^

internal heat transfer coefficient, W/m2/K

hX

vapor convective heat transfer coefficient, W/m2/K

hi

tube wall thermal conductivity, W/m/K

∆p4q

latent heat of the mixture, J/kg

hio

wet perimeter heat transfer coefficient, W/m2/K

39


Page 40 of 45

Page 41 of 45

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


K

thermal conductivity of the flue gas, W/m2/K

u

molecular weight, kg/kmol

Nw

number of different groups defined by the group contribution method

Pressure, Pa

K s /K t

Pr

Prt /Prs P

thermal conductivity of the vapor/liquid phase, W/m2/K

Prandtl number of the flue gas side Prandtl number of the vapor/liquid phase reduced pressure [P/Pcrit]

q

internal heat flux, W/m2

{

total heat transferred from flue gas in tube zone, W

q0 rY

Re}

Re~ Re

ReX /Re

area parameter of component 7 in the UNIFAC model volume parameter of component 7 in the UNIFAC model Reynolds number of the flue gas homogenous Reynolds number in mist flow Reynolds number of the flue gas side Reynolds number of the vapor/liquid phase temperature, K

0q /

inlet/outlet temperature of the mixture in each zone, K

Tw

temperature of the flue gas, K

Un*

maximum velocity in the minimum flow area, m/s

v

molar volume, m3/mol

.

./4 ./0 .0 Y

0

mass velocity of the vapor phase, kg/m2s vapor quality dryout completion quality dryout inception quality mass (mole) fraction of ith species in the liquid phase multiplication factor mass (mole) fraction of ith species in the vapor phase

Greek symbols αY γY δ ε

fraction of molecule i occupied by group k

liquid activity coefficient of component i liquid film thickness, m cross-sectional vapor void fraction 40



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

ξ

ρ

symbol for property density, kg/m3

̅

density of the vapor-liquid mixture, kg/m3

ρX /ρt

vapor/liquid density, kg/m3

θ

dry angle of the tube perimeter, rad

θ

dry angle of the tube perimeter in stratified flow, rad

λ

μ

friction coefficient dynamic viscosity, Ns/m2

τ

interaction parameter between two groups in the UNIFAC model

φ}

dimensionless diagonal pitch of the staggered arrangement

υY φt

number of group k in component i in the UNIFAC model dimensionless longitudinal pitch

φ

dimensionless transverse pitch

AYt /AY

distilled volume or mass fraction of the left/right of interval i

φXY Ω

vapor phase fugacity of component i

internal cross sectional surface area, m2

41


Page 42 of 45

Page 43 of 45

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


Acknowledgments This work was supported by the National Natural Science Foundation of China (61333010, 61590923, 61422303).

References (1) Oprins, A. J. M.; Heynderickx, G. J. Calculation of three-dimensional 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. Three-dimensional asymmetric flow and temperature fields in cracking furnaces. Ind. Eng. Chem. Res. 2001, 40, 5087-5094. (3) 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, 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) Hu, G.; Wang, H.; Qian, F. 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

42



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

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.; Sage, R. Crude split figured by computer. Hydrocarbon Process. 1969, 48, 143-&. (25) Eckert, E. Do we need pseudocomponents? Chem. Listy, Prague 1997, 95. (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. China 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 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 High-Boiling-Point Petroleum Fractions (350°C+). Oil Gas Sci. Technol. 2012, 67, 451-477.

43


Page 44 of 45

Page 45 of 45

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


(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.; Brana-Mulero, 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.

44


Simulation of the Evaporation Tube Banks in the Convection Section

Recommend Documents