Modeling the Coke Formation in the Convection ... - ACS Publications

May 7, 2010 - The presence of liquid hydrocarbon droplets in the convection section of a steam cracker may cause serious fouling problems due to coke ...
0 downloads 0 Views 4MB Size
5752

Ind. Eng. Chem. Res. 2010, 49, 5752–5764

Modeling the Coke Formation in the Convection Section Tubes of a Steam Cracker Sandra C. K. De Schepper, Geraldine J. Heynderickx,* and Guy B. Marin Laboratory for Chemical Technology, Ghent UniVersity, Krijgslaan 281 (S5), B9000 Gent, Belgium

The presence of liquid hydrocarbon droplets in the convection section of a steam cracker may cause serious fouling problems due to coke formation, especially in high temperature zones. In order to investigate these fouling problems, a model has been developed and implemented in a CFD code to accurately simulate the behavior of a hydrocarbon droplet impinging on a hot surface. The impact energy of the droplet and the hot surface temperature are found to determine the impact behavior. On the basis of the newly developed model, the positions where droplets preferentially collide with the convection section tube walls and liquid material is deposited are determined. Furthermore, a kinetic model describing coke formation out of liquid hydrocarbon droplets in the temperature range of 450-700 K has been developed to calculate the rate of coke formation in each zone of the convection section tube. The calculated coke layer thickness on the most vulnerable tube locations and the industrially available values correspond remarkably well. 1. Introduction Steam cracking of hydrocarbons is the most important source of olefins and aromatics, the main feedstock in the petrochemical industry. A steam cracking unit consists of two main sections. The radiation section is fired with side-wall radiation burners and/or long-flame floor burners to supply the necessary heat for the cracking reactions taking place in the reactor tubes suspended in the middle of the steam cracker radiation section. Up to 95% of the heat transferred from the flue gas section to the reactor tubes is due to radiation. However, the temperature of the flue gas leaving the radiation section is still on the order of 1500 K. The remaining flue gas energy is used in the convection section of the steam cracker to evaporate the liquid hydrocarbon feed, to produce overheated steam, and to overheat the hydrocarbon-steam mixture before introducing it into the reactor coils suspended in the radiation section. The heat transfer to the heat exchangers in the convection section is mainly due to convection. Finally, hardly 5% of the energy fed to the radiation section of the steam cracker is lost via the convection section outlet chimneys to the environment.34 Comparable to the coking problems in the tubes of the radiation section,19,20 fouling problems can occur in the tubes of the convection section as well, especially when heavy liquid hydrocarbon feeds are used. This can result in an unexpected shut-down of the entire steam cracker. The damage can result in the need to replace the convection section tubes. The economic losses are considerable. In order to prevent these shutdowns, the origin of the coke formation in the convection section tubes has to be studied. This is done in the present study by modeling the coke formation in the convection section tubes. This way, the positions in the convection section tubes that are vulnerable for coke formation are determined. On the basis of that knowledge, the operating conditions of the steam cracker are to be adapted in order to avoid or reduce the coke formation problem at those positions. Results of previous studies,10,34 performed by the authors, have shown that the liquid hydrocarbon feed is not completely evaporated when leaving the vaporizer, located at the top of the convection section. The twophase vapor-liquid hydrocarbon flow flowing out of the * To whom correspondence should be addressed. E-mail: [email protected].

vaporizer is mixed with overheated steam. The resulting hydrocarbon-steam mixture flow is of the spray flow type: the fraction of not yet evaporated liquid hydrocarbons is present as small droplets. Finally, in the mixture overheater located at the bottom of the convection section, the feed is completely vaporized and further heated. The liquid hydrocarbon droplets that enter the lower heat exchanger are expected to be at the origin of the coke formation inside these convection section tubes. These droplets containing the heavy feedstock hydrocarbons are found to stick easily onto the tube walls. Coke is then expected to be formed out of these sticking droplets through dehydrogenation and polymerization reactions. In order to obtain further insight into coke formation in the convection section, it is very important to determine where droplets preferentially collide with the tube walls. The study of the coke formation thus requires the modeling of the spray flow behavior on the one hand and the application of a coke formation model on the other hand. 2. Modeling 2.1. Spray Flow Modeling. For the modeling of liquid droplets in a continuous gas phase, a Eulerian-Lagrangian approach is used.1,2 For the continuous gas phase, the well-known set of Reynolds-Averaged Navier-Stokes equations (RANS) is solved (see Table 1).1,3 Due to the averaging, the so-called Reynolds stress terms appear in the equations. These are unclosed terms that need to be modeled by introducing a turbulence model. The Reynolds Stress Model (RSM)4 is used in the present study to simulate the turbulence effects in the mixture overheater. It is the most elaborate turbulence model available in the commercial software package Fluent.5 It closes the RANS equations by solving the transport equation for each of the six Reynolds stresses (see Table 2) appearing in the momentum equations due to the averaging, together with a transport equation for the dissipation rate of the turbulent kinetic energy. For details on these equations, reference is made to Ranade2 and the Fluent User Guide.5 The RSM is chosen because it accounts for the effects of streamline curvature, swirl, rotation, and rapid changes in the strain rate in a more rigorous manner than the well-known standard k-ε closure model.6 Moreover, it is extensively used

10.1021/ie100091e  2010 American Chemical Society Published on Web 05/07/2010

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

The droplet size is found by solving (for the x direction in Cartesian coordinates):

Table 1. Reynolds-Averaged Navier-Stokes eqs 1-3 3

∑ ∂x∂ · (Fuj ) ) S

∂ (F) + ∂t

∂ (Fuj ) + ∂t i ∂p + ∂xi

i j

j

j)1

∑ ∂x∂

j

j)1

[( µ

∂uji ∂ujj 2 + - δij ∂xj ∂xi 3

3

∂uj1

∑ ∂x

1

1)1

)]

+



∑ ∂x∂ · (FΕuj ) ) ∑ ∑ ( ∂x∂ (τ 3

3

3

i)1

j)1

j

j)1

(2)

3

j)1

∂ (FE) + ∂t

(1)

M

dxp ) up dt

(9)

3

∑ ∂x∂ · (Fuj uj ) ) 3

-

i

i

j)1

5753

ij

j

∂ (-Fu'iu'j) + SF,i ∂xj

)

- Fu'iu'j)uji -

j

3

∑ ∂x∂ q + S j

j)1

E

(3)

j

Table 2. Transport Equations for Reynolds Stresses (eq 4) ∂ (Fu'iu'j) ∂t

+

(local time derivative)

∂ (Fu u' u' ) ) ∂xk k i j (Cij≡convection)

[

]

∂ ∂ ∂ µ (u' u' ) - [Fu'iu'ju'k + p(δkju'i + δiku'j)] + ∂xk ∂xk ∂xk i j

(

(DT,ij≡turbulent diffusion)

- F u'iu'k

(

)

(DL,ij≡molecular diffusion)

∂u'j ∂u'i + u'ju'k - Fβ(giu'jθ + gju'iθ) ∂xk ∂xk (Gij≡bouyancy production)

(Pij≡stress production)

)

(4)

∂u'i ∂u'i ∂u'j ∂u'j +p + - 2µ ∂xk ∂xi ∂xk ∂xk (φij≡pressure strain)

(εij≡dissipation)

when particle trajectories are to be calculated. As one of the goals in the present work is to follow a hydrocarbon droplet in a convection section tube, this model is selected to close the RANS equations for the continuous gas phase. The velocity of a discrete droplet is calculated by integrating the force balance on the droplet, written in a Cartesian reference frame as

(5)

dVp gy(Fp - F) 18µ CDRed,y (V - Vp) + ) × + Fy 2 dt 24 Fp Fpdp

(6)

∆mp m ˙ mp,0 p,0

(10)

Momentum exchange: SF,i )

SE )

dwp gz(Fp - F) 18µ CDRed,z (w - wp) + ) × + Fz (7) 2 dt 24 Fp Fpdp Notice that the instantaneous velocities u, V, and w of the continuous gas phase are used in the above expressions. The three contributions on the right-hand side (RHS) of eqs 5-7 represent the drag force, the gravity force, and additional forces (e.g., the virtual mass force, the thermophoretic force, the Brownian force, and the Saffman lift force).2 Red is the relative Reynolds number, which is defined as (for the x direction in Cartesian coordinates): Fdp |up - u| µ

SM )

(

18µCDRe 24Fpd2p

)

(up,i-ui) + Fother m ˙ p∆t

(11)

Heat exchange:

dup gx(Fp - F) 18µ CDRed,x (u - up) + ) + Fx × 2 dt 24 Fp Fpdp

Red,x ≡

Equations 5 and 9 are solved simultaneously to determine the velocity and the position of the droplet at any given time in the x direction. An analogous approach is followed for the y and z directions. As such, the trajectory of the droplet is determined. As the flow of the continuous gas phase is turbulent, the turbulent dispersion of the droplets is to be accounted for. This is done by applying a random walk model2 to determine the fluctuating component of the continuous phase velocity. Finally, coupling between the discrete and the continuous phase is required. To determine the required amount of coupling, a regime map is used (see Figure 1).2 When the droplet size is very small or the liquid volume fraction is very low, the gas phase flow field is not affected by the droplet trajectories, and the influence of the discrete phase on the continuous phase can be neglected: one-way coupling suffices. However, in the case of an increased discrete phase volume fraction, the discrete and continuous phases influence one another, and two-way coupling is required. Finally, if the number of droplets is very high, interaction between the droplets has to be accounted for as well. This is referred to as four-way coupling. On the basis of the data obtained by De Schepper et al.,10 a two-way coupling is required as the volume fraction of the discrete phase, R2, equals 10-4. While the droplet trajectory is calculated, the software keeps track of the mass, momentum, and heat gained or lost by the droplet. The corresponding equations are Mass exchange:

(8)

[

∆mp ∆mp c ∆T + (-hfg + mp,0 p p mp,0



Tp c Tref p

]

dT) m ˙ p,0 (12)

The above expressions are used as source terms in the conservation equations for the continuous gas phase (see Table 1). Two-way coupling implies alternating solving of the equations describing both phases until convergence is obtained.2 The latter implies that droplet-droplet interactions are not accounted for. The main reason is the low volume fraction of the discrete phase (R2 ) 10-4), resulting in a dilute suspension. 2.2. Modeling of the Droplet Impact Behavior. The droplets that impinge on the internal tube surface are expected to be at the origin of the coke formation. However, before the coke formation can be modeled and calculated, the impact behavior of the droplets on the tube surface has to be analyzed. Wachters and Westerling11 described different collision models where water droplets either stick, rebound, or splash on the surface (Figure 2). Mundo et al.12 confirmed this behavior for ethanol droplets and for water-ethanol-sucrose droplets. Gavaises et al.13 described the impact behavior of diesel droplets impinging on a hot surface, and the results are similar to those

5754

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

Bai and Gosman15 defined the stick behavior as “a wall impaction in which droplets with little energy impact on a relatively cool surface and remain there in a nearly spherical state.” The rebound impact behavior was described by Kandlikar and Steinke14 as a situation where the droplets experience an (in)elastic reflection from the wall. The experimental results of Wachters and Westerling11 (see Figure 4), concerning this impact behavior, were curve-fitted by Mundo et al.,12 resulting in the following equation to calculate the normal velocity of the rebounded droplet: Weout ) 0.6785Wein exp(-0.04415Wein) Figure 1. Coupling between phases in dispersed multiphase flows.2

(14)

The tangential velocity remains unchanged. Splashing of a droplet occurs when the impact energy (i.e., the Weber number) of the impinging droplet is high. The droplet breaks into several smaller droplets, which, in turn, are rebounded in different directions. Only a small part of the droplet mass sticks to the tube surface. As the simulation of this model is very time-consuming, a simplified model has been developed by Grover and Assanis,16 limiting the number of rebounded, smaller droplets to three, next to a single droplet sticking to the surface: min ) ms + mwf ) m1 + m2 + m3 + mwf

(15)

The total mass of the three splashed droplets, ms, is estimated from empirical correlations: 0.0 We < 80 ms ) 0.5 80 < We < 600 min 0.75 We > 600

(16)

Upon splashing, the tangential momentum of the impinging droplet is divided over the three newly formed droplets. In the model of Grover and Assanis,16 30% of the tangential moment is transferred to the wall droplet. The resulting motion of the droplet over the surface is however not considered in this work. According to Grover and Assanis,16 the tangential momentum conservation is thus written by:

Figure 2. Droplet impact behavior regimes.

0.7minVT,in ) m1VT,1 + m2VT,2 + m3VT,3

(17)

0.3minVT,in ) mwfVT,wf

(18)

Finally, an energy conservation equation is required, accounting for the kinetic energy, the surface energy, and the dissipated energy. The energy of the impinging droplet is equated to the energy of the droplets that are formed: Figure 3. Decision chart for determination of droplet impact behavior.

of Wachters and Westerling11 and Mundo et al.12 The impact behavior of heavy hydrocarbon droplets has not yet been described in the literature, but it is reasonable to suppose the behavior is similar to that of diesel droplets. The variables determining the droplet impact behavior are the wall temperature, compared to the boiling temperature of the droplet, and the impact energy of the droplet, expressed by the Weber number: FpV2ndp We ) σ

(13)

A decision chart on the impact behavior of a droplet is shown in Figure 3.

( )

We0.83 dmax 1 in minV2in + πd2inσ ) 0.2 πd2inσ 0.33 2 din Re (Esurf,in) in (Evisc)

(Eke,in) 3

+

∑ 21 m V

2 j j

j)1 (Eke,s)

2

1 2 + mwfVwf 2 (Eke,wf)

3

+

∑ πd σ 2 j

j)1 (Esurf,s)

(19) An expression for the energy dissipated upon impact, Evisc, has been developed by Mao et al.17 According to Pasandideh-Fard et al.,18 the maximum extension diameter for the droplet can be calculated from: dmax ) 0.5 dinRe0.25 in

(20)

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

5755

Figure 4. Relationship between the normal velocity components at arrival and departure of a droplet impinging on a hot wall.11

For the calculation of the energy of the newly formed droplets, reference is made to Grover and Assanis.16 Next to the spray flow modeling, a coke formation model is required to calculate the amount of coke formed out of the liquid material deposited on the tube surface. 2.3. Coke Formation Modeling. Coke formation is inherent to the processing of hydrocarbons. In steam cracking, the coke deposit on the wall of the reactor coils affects the cracking process by a pressure drop increase over the reactor coil, a heat flux decrease from flue gas to process gas, hot spots on the tube wall, corrosion due to carbonization of the reactor coils, and so forth. The major part of the research on coke formation in the steam cracking process therefore concentrates on the process conditions in the reactor coils.19,20 Temperatures in the convection section tubes are however considerably lower,34 and it is thus not justified to use the kinetic coke formation models, applied in the reactor coils, in the convection section tubes. A literature review has been done to find information on coke formation in convection section tubes and corresponding kinetic coke formation models in the temperature range of the mixture overheater tube wall temperatures, 450 to 700 K. Specific information on coke formation in the convection section tubes of a steam cracker unit is found in a patent.21 A new kind of nozzle for the vaporization of a heavy hydrocarbon feedstock with steam is described. The patent describes how the normally used nozzles easily give rise to coke formation in the flow line, especially in that part of the line where steam is introduced for the final vaporization of the hydrocarbon feedstock. The following phenomena are assumed to occur. When the partially vaporized hydrocarbon feed is brought in contact with superheated steam to complete the vaporization, the temperature of the droplets will immediately rise to temperatures close to the overheated steam temperature. The remaining part of the lighter components in the liquid droplets will immediately vaporize, preferentially at the droplet surface. As a consequence, an external layer of even heavier liquid hydrocarbons is formed after evaporation of these lighter hydrocarbons. Further evaporation of the liquid hydrocarbons in the center of the droplet is limited by the slow diffusion of the lighter hydrocarbons through the outer layer of heavy liquid hydrocarbons. This physical inability of the droplets to completely evaporate can lead to coke formation when the droplet temperatures remain high. It is thus clearly described in the patent that the phenomenon of coke formation in the convection section tubes is related to liquid hydrocarbon droplets sticking on the hot walls. However, no information concerning a kinetic coke formation model is given.

Table 3. Nine-Step Global Kinetic Model27 nine-step global kinetic model Bulk Reactions: O2 + Fuel f ROOH

(21)

ROOH + Fuel f Solubles

(22)

ROOH + Fs f P

(23)

P + Fuel f Solubles

(24)

ROOH + Fuel f Dbulk

(25)

Dbulk + Fuel f 2Dbulk

(26)

Wall Reactions O2 + Fuel f D

(27)

PfD

(28)

Dbulk f D

(29)

The coke formation out of fuel droplets impinging on a hot surface can be related to the thermal instability of the fuel. A review of the thermal instability of fuel was written by Katta et al.22 Nevertheless, the coking mechanism is still largely unknown. Most of the modeling is done using so-called global chemistry models, which try to postulate a kinetic rate expression based on non-elementary reaction steps. From a one-step model,23 to a two-step mechanism,24 and a three-step mechanism accounting for mass transfer and diffusion effects24 further improved by Krazinski et al. in 1992,26 the most elaborate model is the nine-step global mechanism proposed by Katta et al.,27 consisting of six bulk reactions and three wall reactions (see Table 3). For the work considered here, one of the wall reactions, eq 28 (see Table 3) was selected from the global model. As coke formation takes place on the wall, the bulk reactions are not considered. Furthermore, the hydrocarbon/steam mixture is oxygen-free, while coke formation from a bulk component like Dbulk cannot be brought into agreement with the droplet impact model that is used. In the remaining reaction, the coke precursor P is transformed into a coke deposit D, with a reaction rate given by

5756

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

r ) A e-Ea/RT[P] (kg/m2s)

(30)

where A ) 260 m/s Ea ) 71.18 kJ/mol The values for the kinetic parameters were obtained by fitting the experimental results obtained by Ervin and Williams28 in the temperature range of 298 to 546 K, closely corresponding to the range of the wall temperatures of the convection section tubes. The validity of the model and the parameters was demonstrated by Spadaccini et al.29 for several fuels. More recently, Kuprowicz et al.30 developed an 18-elementary-step reaction mechanism replacing the six bulk reactions of the model of Katta et al.27 (see Table 3). They were however not able to determine elementary reaction steps for the wall reactions. The coke precursor P is defined as the hydrocarbon component(s) with the highest tendency to form a coke deposit. Given the vaporization of the feed and the results obtained by Kopinke et al.,31,32 the aromatic components in the hydrocarbon feed are selected. The kinetic model of Ervin and Williams28 (eq 30) is implemented in an extensive User Defined Function (UDF) added to the commercial software package Fluent5 to calculate the coke formation, in combination with the droplet trajectory calculation. The main remark to be made is that the kinetic model development was done in the presence of oxygen, a component that is not found in the hydrocarbon-steam mixture in the convection section tubes. 3. Solution Techniques The calculations are performed making use of a second-order upwind calculation scheme1 for the solution of the conservation equations and the RSM model equations. The PISO algorithm33 is used to couple pressure and velocity. Convergence is assumed to be reached when the scaled residuals have lowered in order of magnitude by a factor of 4, and by 6 orders in magnitude for energy. The number of grid cells required to obtain grid-size independent simulation results has been determined by performing simulations on different grids. The grid is characterized by the definition of three layers of small cells near the wall to correctly capture the phenomenon of droplet-wall interactions. 4. Convection Section Geometry and Operating Conditions For an overview and discussion of the complete convection section geometry, reference is made to De Schepper et al.10,34 In this work, the results of a Eulerian-Lagrangian simulation for the eight tubes of mixture overheater 1 (see Figure 5) with a geometry as shown in Figure 6 are presented. The tubes of mixture overheater 1 each make three horizontal passes through the convection section. The passes are connected by two (isolated) U bends situated outside the convection section. As the final goal of the research is to have an overview of the exact position where a droplet is deposited and an exact knowledge of the amount of coke formed at that position, the heat exchanger tubes are divided in several zones and subzones. Smaller zones are defined at those positions where more droplet impingement and thus more coke formation is to be expected. Along a three-pass tube, 79 zones are defined, each having six subzones.

Figure 5. Cross-section of the convection section with indication of the heat exchangers and the tube positions.

5. Results and discussion First, a simulation is performed to calculate the flow and temperature fields of the gas phase without droplets in the heat exchanger tubes.10,24 Next, discrete droplets are injected at the inlet of a tube of mixture overheater 1. The droplet conditions to be defined are the inlet velocity, diameter, and temperature. A random introduction of the droplets over the cross section of the tube is required. A solid cone injection is applied: velocity and direction of the droplets are randomly generated within a defined cone at the tube inlet. Notice that the use of the commercial software package Fluent5 implies that neither a boiling point trajectory for the hydrocarbon droplets nor a temperature dependence of the physical properties of the hydrocarbon droplets can be taken into account. A coupled simulation of the entire steam cracker convection section, that is, flue gas side and heat exchanger tubes, has been performed, accounting for these temperature dependencies.34 In a later stage, both studies will be coupled. Table 4 gives an overview of the different simulations that have been performed next to the base case simulation, in order to determine the influence of the injected droplet diameter and the tube position in the mixture overheater 1 of the convection section cross section. Finally, a case study is performed where a liquid spray with a droplet diameter distribution is injected into the different tubes. Results for the base case simulation using 1 µm droplets, during which it is observed that only the stick impact behavior of the droplets is observed (see Figure 3), are shown in Figure 7. Most of the 1 µm droplet impact is calculated for the inlet tube zones 4-7, in combination with smaller peaks for zones 28 and 54, the first and second U bends of the three-pass tube, respectively. It must thus be concluded that, in the inlet bend and in the U bends, the (small) droplets are not able to completely follow the flow field of the continuous gas phase.

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

5757

Figure 6. Definition of the (sub)zones in a convection section tube. Table 4. Overview of the Performed Simulations and the Applied Conditions

droplet diameter (µm) droplet density (kg/m3) droplet boiling point (K) droplet surface tension (N/m) tube position

base case

influence droplet diameter

influence tube position

case study

1 718 622 0.018 1

0.1, 1, 10, 50, 100 718 622 0.018 1

0.1, 1, 10, 50, 100 718 622 0.018 1f8

diameter distribution 0.1 f 100 µm 718 622 0.018 1f8

Furthermore, only the stick impact behavior takes place due to the low values of the Weber number (see Figure 3), that is, due to the low impact energy of a small 1 µm droplet. Figure 8 shows the coke formation rates calculated using eq 30 for tube pass 1 and part of tube pass 2. For the following tube zones, this coke formation rate drops to values of almost zero. In Figure 8, it is seen that the coke formation rate reaches the highest values in the bottom part of the tube and the lowest values in the upper part of the tube. This corresponds with the results presented by De Schepper et al.10 There, it was calculated that the temperature of the flue gas that flows from the bottom to the top of the convection section, over the horizontally suspended tubes, has its highest value at the bottom of the tube and its lowest value at the top of the tube. Consequently, the same observation was made for the tube wall temperatures. In a set of simulations, the influence of the droplet diameter has been investigated. Figure 7 presents the results for the behavior of the droplets with different diameters in the tube passes 1 to 3, studied for each diameter separately. In Figure 7a, the droplet behavior in tube pass 1 is studied. As expected, it becomes more difficult for the droplet to follow the gas flow with growing droplet diameter. For the largest droplets (50 and 100 µm), the splash model is often applied for droplets that impact on the zones of the inlet bend, contrary to the observation for a small droplet where the stick model has to be applied as discussed above. Especially for the 100 µm droplets, the impact energy is high, resulting in a high value of the Weber number (see Figure 3). More liquid is splashed back into the gas flow as compared to the 50 µm droplets. This explains why, in the inlet zone, somewhat more liquid originating from the 50 µm droplets is deposited on these zones as compared to the 100 µm droplets, contrary to what should be expected at first. Due to the splashing, the deposition of liquid originating from the 100 µm droplets is spread out over the inlet zone as the splashed, smaller diameter droplets originating from these large droplets

are deposited at a further position along the tube. The smaller the droplet diameter, the better the droplet can follow the gas flow. As a result, the fraction of small droplets deposited in the inlet bend of the tube is small as compared to the large droplets. The fraction of the liquid originating from droplets with (originally) different diameters that is deposited in the first U bend (zone 28) is equal for each initial droplet diameter. Note that the deposition of liquid mass corresponding with the initially largest droplets is in fact corresponding with the deposition of liquid from smaller droplets formed out of the splashing of initially large droplets. In Figure 7b, it is seen that at the end of tube pass 1 all liquid originating from initially large droplets has already been deposited. Some (originally) small droplets are still able to follow the gas flow. From Figure 7c, it is determined that none of the injected liquid, not even when injected as small droplets, is able to leave the tube, as no liquid is deposited in the last tube zone 79. Note that, for each droplet diameter, the sum of the deposited liquid fractions over all tube zones equals 100%. Coke formation rate profiles for the different droplet diameters are presented as well. Figure 9 presents the coke formation rate profiles for the first tube pass. The coke formation rate for the largest droplets is extremely high in the inlet bend but drops with rising distance in the tube due to the splashing of the large droplets into smaller droplets. Along the tube pass, the coke formation rate of the smaller droplets becomes even higher than that of the (originally) large droplets because the small droplets that hit the tube wall will not splash but stick to the wall and are turned into coke. When analogous curves are studied for tube pass 2, the coke formation rate for the (originally) largest droplets (10-50 µm) falls to zero in the zones following the tube bend (that is zone 29 and following). For the smallest droplets, coke formation rates are calculated all over tube passes

5758

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

Figure 7. Fraction of liquid deposited in each tube zone of tube pass 1 (a), tube pass 2 (b), and tube pass 3 (c) for the five droplet diameters

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

5759

Figure 8. Coke formation rate (kg/m2s) for the base case, calculated with the kinetic modelseq 30.28

Figure 9. Coke formation rate (kg/m2s) for different droplet diameters, calculated with the kinetic modelseq 30.28

2 and 3, but the absolute values drop by a factor of 10 when comparing the values for pass 2 and pass 3 with the values for pass 1. In a next set of simulations (see Table 4), the influence of the position of the tube in a horizontal tube row is studied. However, no significant influence is observed. This result was to be expected, accounting for the results of De Schepper et al.,10,34 where it was shown that a changing tube position results in a slightly different temperature of the inner tube wall and a slightly changing amount of liquid input. Finally, a case study is performed where a hydrocarbon spray flow with droplets of different diameters is injected in the tube. The droplet diameter distribution is defined in Figure 10. The simulation results for the eight tubes on a horizontal row in the convection section of a steam cracker are transformed into data presenting the coke layer thickness in each tube when 30 days of operation have passed. Note that this is an estimated value based on the initial coking rate. An exact value requires a run length simulation of the convection section, taking into account the heat transfer resistance of the growing coke layer.35,36 In Figure 11, the results for the tubes in the different positions are compared. The results are shown for tube zones 1-35, as the

Figure 10. Droplet diameter distribution for the case study.

coke formation in the following zones is negligible. It can be seen that the profiles of the calculated curves are quite similar for the tubes in the different positions. This is explained by the fact that the fraction of the injected liquid deposited on corresponding zones of the different tubes is comparable. Furthermore, it is seen that, for all tubes, independent of their

5760

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

Figure 11. Continued.

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

Figure 11. Continued.

5761

5762

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

Figure 11. Coke layer thickness in the tubes of mixture overheater 1 after 30 days: (a) position 1 to (h) position 8.

position in the horizontal row, the thickest coke layer is formed on the bottom part of the tube, in particular for zone 3, which is the zone following the inlet bend of the tube. Thus, a first conclusion is that the location in the tube following the tube inlet of the tubes of mixture overheater 1 will suffer the most from the formation of coke. The tube zones following tube zone 3 also suffer from a high coke formation, but the cokes is seen to be formed all around the tube circumference and not only on the bottom part of the tube. Next, it is observed that the zones in the (first) U bend (zone 28) and following this U bend (zone 29) are critical zones for coke formation as well. The differences in coke layer thickness for the horizontal tube positions 1 (Figure 11a) to 8 (Figure 11h) can be explained by the differences in tube wall temperatures, on one the hand, and the different values for the total inlet mass flow of the droplets, on the other hand. The inlet mass flow of liquid droplets in the different tubes changes, as the fraction of evaporation of the feed in the vaporizer and the temperature of the overheated steam varies as well.10,36 The mass flow rates of the liquid

Table 5. Mass Flow Rate of Liquid Hydrocarbon Feed at the Entrance of the Mixture Overheater 1 Tubes at Different Positions liquid mass flow rate (kg/s) tube tube tube tube tube tube tube tube

position position position position position position position position

1 2 3 4 5 6 7 8

0.001125 0.001750 0.002125 0.002875 0.003125 0.003375 0.002750 0.002500

hydrocarbons at the entrance of the tubes at the different horizontal positions are given in Table 5. For the tubes in the outer positions, the tube wall temperatures are high. This results in a higher rate of coke formation than for the central tubes. However, the amount of liquid introduced in these outer tubes is low compared to that of the central tubes (see Table 5), resulting in a relatively lower liquid deposition on the walls of the outer tubes. Two opposite effects are thus observed. When

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

comparing the results for the tubes in different horizontal positions, it is seen that the thickest coke layer obtained after 30 days is calculated for the tube located in position 4 (see Figure 11d). For this tube, intermediate values for the liquid inlet mass flow rate as well as for the tube wall temperatures are calculated. After 30 days, a coke layer with a thickness of 0.6 mm is formed on the bottom part of zone 3 of the tube in position 4 (Figure 11d). Due to the low liquid inlet mass flow of droplets in the tube located in position 1 (see Table 5), the coke layer in this tube is the thinnest (Figure 11a). The amount of liquid introduced is the highest (see Table 5) for the tube in position 6. However, the tube wall temperatures are the lowest, explaining the lower amount of coke (see Figure 11f) formed in this tube as compared to the tube in position 4 (see Figure 11d). It must be remarked that the calculated coke layer thickness values correspond remarkably well with values supplied by the industry. The obtained models can be used in a coupled simulation of the steam cracker convection section34 to predict the exact growth of the coke layer as a function of time.35,36 Furthermore, the model can be used to adapt the operating conditions and/or the geometry of the convection section of a steam cracker to avoid or reduce the problem of coke formation in the convection section tubes. 6. Conclusions The formation of coke out of liquid droplets that stick on the wall of the convection section tubes of a steam cracker has been studied and modeled. Two models are required: a model that describes the impact behavior of the droplet and a model that describes the coke formation rate. These models are combined with one another and with a commercial CFD software package to simulate the vapor-liquid two-phase flow patterns using a Eulerian-Lagrangian approach. Depending on the boiling temperature of the hydrocarbon droplets and their impact energy on the tube wall, expressed by the Weber number, three different droplet impact behaviors are observed and modeled: stick, splash, and rebound. These models allow determination of the locations in a tube that are vulnerable for liquid deposition and thus for coke formation. The coke formation rate is described by using a kinetic model valid in the heat exchanger wall temperature range of 450 to 700 K. Using all these models, it is found that the droplet diameter has a large influence on the droplet behavior and on the location in the tube where the droplet is deposited. The calculated coke layer thickness corresponds remarkably well with the industrially observed values. Acknowledgment Sandra De Schepper would like to thank SABIC Europe for supporting this research study. The CFD research in thermal cracking (G.J. Heynderickx) is supported by the Fonds voor Wetenschappelijk Onderzoek-Vlaanderen (FWO-Vlaanderen) (project G.0022.09). Nomenclature Roman Symbols A ) pre-exponential factor, m/s cp ) heat capacity, J/kgK din ) diameter of the impinging droplet, m dj ) diameter of the splashed droplet, m dmax ) maximum extension diameter, m

5763

dp ) droplet diameter, m E ) total energy per unit mass, J/kg Ea ) activation energy, kJ/mol Eke,in ) kinetic energy of the impinging droplet, kg m2/s2 Eke,s ) kinetic energy of the splashing droplets, kg m2/s2 Eke,wf ) kinetic energy of the wallfilm droplet, kg m2/s2 Esurf,in ) surface energy of the impinging droplet, kg m2/s2 Esurf,s ) surface energy of the splashing droplets, kg m2/s2 Evisc ) energy of the droplet dissipated upon impact, kg m2/s2 F ) external body forces, kg/m2s2 Fi ) forces acting on the droplet in the i direction with i ∈ (x,y,z), N Fx ) forces acting on the droplet in the x direction, N Fy ) forces acting on the droplet in the y direction, N Fz ) forces acting on the droplet in the z direction, N g ) gravitational acceleration, g ) 9.81, m/s2 gx ) gravitational acceleration in the x direction, m/s2 gy ) gravitational acceleration in the y direction, m/s2 gz ) gravitational acceleration in the z direction, m/s2 hfg ) latent heat, J/kg m ) mass, kg min ) mass of the impinging droplet, kg mp ) mass of the droplet, kg ms ) total mass of the three splashing droplets, kg mwf ) mass of the wallfilm droplet, kg m1 ) mass of splashing droplet 1, kg m2 ) mass of splashing droplet 2, kg m3 ) mass of splashing droplet 3, kg mp,0 ) original mass of the droplet, kg m ˙ p ) mass flow rate of the droplets, kg/s P ) pressure, Pa [P] ) mass concentration of the coke precursor, kg/m3 p ) pressure, Pa qj ) conductive heat flux in the xj direction, J/m2s Red,x ) Reynolds number based on droplet diameter and relative velocity ) (F dp|up - u|)/µ (for the x direction) r ) reaction rate, kg/m2s SE ) source term in the energy equation, J/m3s SF ) source term in the momentum conservation equation, kg/m2s2 SM ) source term in the mass conservation equation, kg/m3s SRk ) source term for volume fraction of phase k, 1/s T ) temperature, K Tboiling ) boiling temperature (Figure 2), K Tp ) droplet temperature, K Tw ) wall temperature, K t ) time, s u ) velocity of the continuous phase in the x direction, m/s ui ) velocity component in the xi direction, m/s ui′ ) Reynolds average of ui, m/s ui′ ) Reynolds fluctuation of ui, m/s uj ) velocity component in the xj direction, m/s uj′ ) Reynolds average of uj, m/s uj′ ) Reynolds fluctuation of uj, m/s uk ) velocity of phase k, m/s up ) droplet velocity in the x direction, m/s V ) velocity of the continuous phase in the y direction, m/s Vj ) velocity of splashed droplet, m/s Vn ) droplet velocity component normal to the surface, m/s Vp ) droplet velocity in the y direction, m/s VT ) tangential droplet velocity component, m/s VT,in ) tangential velocity component of the impinging droplet, m/s VT,wf ) tangential velocity component of the wallfilm droplet, m/s VT,1 ) tangential velocity component of the splashing droplet 1, m/s

5764

Ind. Eng. Chem. Res., Vol. 49, No. 12, 2010

VT,2 ) tangential velocity component of the splashing droplet 2, m/s VT,3 ) tangential velocity component of the splashing droplet 3, m/s We ) Weber number ) (FpV2n dp)/σ w ) velocity of the continuous phase in the z direction, m/s w ) dimensionless coke deposit, wt % wp ) droplet velocity in the z direction, m/s x ) Cartesian coordinate, m xi ) xi direction or coordinate, with i ∈ (x, y, z), m xj ) xj direction or coordinate, with j ∈ (x, y, z), m y ) Cartesian coordinate, m z ) Cartesian coordinate, m Greek symbols Rk ) volume fraction of phase k δij ) Kronecker function µ ) dynamic viscosity, Pa · s µk ) dynamic viscosity of phase k, Pa · s F ) fluid density, kg/m3 Fp ) density of the droplet, kg/m3 σ ) gas-liquid surface tension, N/m τij ) viscous stresses with (i,j) ∈ (x,y,z), kg/ms2

Literature Cited (1) Versteeg, H. K.; Malalasekera, W. An Introduction to Computational Fluid Dynamics, The Finite Volume Method; Longman Scientific and Technical, Harlow: Essex, U. K., 1995. (2) Ranade, V. V. Computational Flow Modeling for Chemical Engineering; Academic Press: San Diego, CA, 2002; Vol. 5. (3) Anderson, J. D. Computational Fluid Dynamics: The Basics with Applications; McGraw-Hill: New York, 2005. (4) Gibson, M. M.; Launder, B. E. Ground Effects on Pressure Fluctuations in the Atmospheric Boundary Layer. J. Fluid Mech. 1978, 86, 491. (5) Fluent 6.3 User’s Guide; Fluent Inc.: Lebanon, NH, 2006. (6) Launder, B. E.; Spalding, D. B. Lectures in Mathematical Models of Turbulence; Academic Press: London, 1972 (7) Cortes, C.; Gil, A. Modeling the gas and particle flow inside cyclone separators. Prog. Energy Combust. Sci. 2007, 33, 409. (8) Hidayat, M.; Rasmuson, A. Numerical assessment of gas-solid flow in a U-bend. Chem. Eng. Res. Des. 2004, 82, 332. (9) Udaya Bhaskar, K.; Rama Murthy, Y.; Ravi Raju, M.; Sumit, T.; Srivastava, J. K.; Ramakrishnan, N. CFD simulation and experimental validation studies on hydrocyclone. Miner. Eng. 2007, 20, 60. (10) 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. (11) Wachters, H. J.; Westerling, N. A. J. The heat transfer from a hot wall to impinging water drops in a spheroidal state. Chem. Eng. Sci. 1966, 21, 1047. (12) Mundo, C.; Sommerfeld, M.; Tropea, C. Droplet-Wall Collisions: Experimental Studies of the Deformation and Breakup Process. Int. J. Multiphase Flow 1995, 21, 151. (13) Gavaises, M.; Theodorakakos, A.; Bergeles, G. Modeling wall impaction of diesel sprays. Int. J. Heat Fluid Flow 1996, 17, 130. (14) Kandlikar, S. G.; Steinke, M. E. Contact angles and interface behavior during rapid evaporation of a liquid on a heated surface. Int. J. Heat Mass Transfer 2002, 45, 3771. (15) Bai, C.; Gosman, A. D. Development of methodology for spray impingement simulation. SAE Technical Paper 950283; SAE International: Warrendale, PA, 1995.

(16) Grover, R. O.; Assanis, D. N. A Spray Wall Impingement Model Based Upon Conservation Principles. Fifth International Symposium on Diagnostics and Modeling of Combustion in Internal Combustion Engines, 2001; p 551. (17) Mao, T.; Kuhn, D.; Tran, H. Spread and Rebound of Liquid Droplets upon Impact on Flat Surfaces. AIChE J. 1997, 43, 2169. (18) Pasandideh-Fard, M.; Qiao, Y. M.; Chandra, S.; Mostaghimi, J. Capillary effects during droplet impact on a solid surface. Phys. Fluids 1996, 8 (3), 650. (19) Cai, H.; Krzywicki, A.; Oballa, M. C. Coke formation in steam crackers for ethylene production. Chem. Eng. Process. 2002, 41, 199. (20) Wauters, S.; Marin, G. B. Kinetic modeling of coke formation during steam cracking. Ind. Eng. Chem. Res. 2002, 41, 2379. (21) Chandrasekharan, K.; Kloth, A. G. J.; Van Westrenen, J. Apparatus and process for Vaporizing a heaVy hydrocarbon feedstock with steam. International application published under the patent cooperation treaty. International Publication Number WO 01/90277 A1, 2001. Patent assigned to Shell Internationale Research Maatschappij B.V. (22) Katta, V.; Jones, E.; Roquemore, W. Modeling of the deposition process in liquid fuels. Combust. Sci. Technol. 1998, 139, 75. (23) Chin, J.; Lefebvre, A. Influence of flow conditions on deposits from heated hydrocarbon fuels, ASME paper 92-GT-114; Presented at the International Gas Turbine and Aeroengine Congress and Exposition, Cologne, Germany, 1995. (24) Giovanetti, A.; Szetela, E. Long term deposit formation in aViation turbine fuel at eleVated temperature, AIAA Paper 86-0525; 24th Aerospace Sciences Meeting, Reno, NV, 1986. (25) Deshpande, G.; Micheal, A.; Solomon, P.; Malhotra, R. Modeling of the thermal instability of aViation fuels; Presented at the 198th ACS National Meeting,Symposium on the chemical aspects of hypersonic propulsion, Miami, FL, 1989. (26) Krazinski, J.; Vanka, S.; Pearce, J.; Roquemore, W. A computational fluid dynamics and chemistry model for jet fuel thermal stability. J. Eng. Gas Turbines Power 1992, 114, 104. (27) Katta, V.; Jones, E.; Roquemore, W. DeVelopment of global chemistry model for jet-fuel thermal stability based on obserVations from static and flowing experiments, AGARD-CP-536, Paper 19, NATO Research and Technology Organisation, 1993. (28) Ervin, J.; Williams, T. Global kinetic modeling of aviation fuel fouling in cooled regions in a flowing system. Ind. Eng. Chem. Res. 1996, 35, 4028. (29) Spadaccini, L. J.; Sobel, D. R.; Huang, H. Deposit formation and mitigation in aircraft fuels. J. Eng. Gas Turbines PowersTrans. ASME 2001, 123 (4), 741. (30) Kuprowicz, N.; Zabarnick, S.; West, Z.; Ervin, J. Use of measured species class concentrations with chemical kinetic modeling for the prediction of autoxidation and deposition of jet fuels. Energy Fuels 2007, 21, 530. (31) Kopinke, F.; Zimmermann, G.; Reyniers, G.; Froment, G. Relative rates of coke formation from hydrocarbons in steam cracking of naphtha. 2. Paraffins, naphtenes, mono-, di-, and cycloolefins, and acetylenes. Ind. Eng. Chem. Res. 1993, 32, 56. (32) Kopinke, F.; Zimmermann, G.; Reyniers, G.; Froment, G. Relative rates of coke formation from hydrocarbons in steam cracking of naphtha. 3. Aromatic hydrocarbons. Ind. Eng. Chem. Res. 1993, 32, 2620. (33) Issa, R. I. Solution of the Implicitly Discretized Fluid Flow Equations by Operator Splitting. J. Comput. Phys. 1986, 62, 40. (34) 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. (35) Plehiers, P. M.; Reyniers, G. C.; Froment, G. F. Simulation of the run length of an ethane cracking furnace. Ind. Eng. Chem. Res. 1990, 29, 636. (36) Heynderickx, G. J.; Froment, G. F. Simulation and Comparison of the run length of an ethane cracking furnace with reactor tubes of circular and elliptical cross sections. Ind. Eng. Chem. Res. 1998, 37 (3), 914.

ReceiVed for reView January 14, 2010 ReVised manuscript receiVed April 6, 2010 Accepted April 12, 2010 IE100091E