A coupled CFD simulation approach for investigating the pyrolysis process in industrial naphtha thermal cracking furnaces

Mohsen Rezaeimanesh,Ali Asghar Ghoreyshi,,S.M.Peyghambarzadeh,Seyed Hassan Hashemabadi

1 Chemical Engineering Department,Babol Noshirvani University of Technology,Shariati Ave.,47148-71167 Babol,Iran

2 Department of Chemical Engineering,Mahshahr Branch,Islamic Azad University,Mahshahr,Iran

3 Computational Fluid Dynamics (CFD) Research Laboratory,School of Chemical,Petroleum and Gas Engineering,Iran University of Science and Technology (IUST),Tehran,Iran

Keywords:Naphtha thermal cracking CFD modeling Systems engineering Residence time Product yield

ABSTRACT In the steam thermal cracking of naphtha,the hydrocarbon stream flows inside tubular reactors and is exposed to flames of a series of burners in the firebox.In this paper,a full three-dimensional computational fluid dynamics (CFD) model was developed to investigate the process variables in the firebox and reactor coil of an industrial naphtha furnace.This comprehensive CFD model consists of a standard k-ε turbulence model accompanied by a molecular kinetic reaction for cracking,detailed combustion model,and radiative properties.In order to improve the steam cracking performance,the model is solved using a proposed iterative algorithm.With respect to temperature,product yield and specially propylene-toethylene ratio (P/E),the simulation results agreed well with industrial data obtained from a mega olefin plant of a petrochemical complex.The deviation of P/E results from industrial data was less than 2%.The obtained velocity,temperature,and concentration profiles were used to investigate the residence time,coking rate,coke concentration,and some other findings.The coke concentration at coil exit was 1.9×10-3%(mass)and the residence time is calculated to be 0.29 s.The results can be used as a scientific guide for process engineers.

1.Introduction

As the main petrochemical process for production of light olefins,the steam thermal cracking of hydrocarbons has a high level of importance for chemical,in special polymer,industries.Although the catalytic cracking process has recently received some attentions,but the main process for production of these olefins is still the thermal steam cracking.Production of olefins is normally carried out through a series of endothermic reactions taking place in long tubes reactor,i.e.the steam cracking coils,embedded in a large long rectangular furnace.In the firebox,the required energy for cracking reactions is provided by a number of gas fired long flame burners in the floor and/or gas-fired radiation burners distributed on the side walls of the furnace.Complexity in modeling of a naphtha cracking furnace arises from the nature of many complicated physical and chemical processes occurring simultaneously such as momentum,heat and mass transfer,chain reactions,and combustion [1].All these processes are directly influenced by the endothermic reactions as well as formation of coke in the reactor coil.In this situation,computational fluid dynamics (CFD) models can be useful but more challenging and computationally expensive[2].

The simulation of these furnaces with its all complexities has always been a big challenge for process engineers and numerous attempts were made in the last two decades.For the simulation of the cracking reactions inside the reactor coils a molecular reaction scheme has been developed for cracking of ethane,propane,butanes and their mixture by Sundaram and Froment [3,4].In 1985,Kumar and Kunzru [5] had presented a molecular reaction scheme for cracking of naphtha based on experimental results.This model was able to predict the product distribution for heavier hydrocarbon (HC) feedstocks such as kerosene,and gas oil [5–7].Although some detailed kinetic radical reaction schemes have been used for better prediction of reaction mechanism in onedimensional (1-D) or two-dimensional (2-D) simulations [2,8–13],due to high computational cost,some studies on CFD simulations of reactor coils used the molecular schemes [14,15].

In the last years of 20th century,full CFD simulations were applied to investigate fireboxes for the first time by de Saegheret al.[16].An attempt for full CFD simulation of both reaction and radiation sides was carried out by Detemmerman and Froment[17].They used zone method of Hottel for the furnace simulation[18].Heynderickxet al.[8] and Oprinset al.[9] have used a twodimensional (3-D) model for a firebox and 1-D numerical model for the coil side of a naphtha thermal cracking furnace with some improvements in radiation models.The turbulent combustion was studied for modeling of firebox.It was shown that detailed reaction kinetics (DRK),such as EBU-SRK (Eddy Break Up (EBU)/finite rate formulation with Simplified Reaction Kinetics (SRK))and EDC-DRK(Eddy Dissipation Concept (EDC) model and makes use of Detailed Reaction Kinetics (DRK)) models,could be used for the furnaces under normal firing conditions [19].Huet al.[2]performed a coupled CFD simulation for an industrial naphtha cracking furnace with SL-II coil type,long-flame and radiation burners.The temperature,concentration,and velocity profiles were obtained and product yields were investigated.Some of the recent studies try to simplified cracking models and predict product distributions rapidly with black boxes [20].Fanget al.[13]showed that a 1-D reactor model with a 2-D data-based recirculation zonal method can predict the behavior of an industrial naphtha furnace better and so faster than a 1-D reactor model with simple and classical CFD simulation of the firebox.But,the detail information about the furnace can help the designers/engineers to configure their plans/plants.

Several radiative heat transfer models such as Adiabatic,P-1,Rosseland,discrete ordinates model (DOM),discrete transfer radiation model (DTRM),and surface-to-surface (S2S) were studied in the fireboxes [12,21,22].Habibiet al.[21] showed that Rosseland model is not valid in a steam cracking firebox and Huet al.[12]recommended the use of DOM for parallel run in FLUENT 14.0 between the others.

For coke formation in the naphtha cracking furnaces,a model was proposed following an experimental study performed by Kumar and Kunzru [23].Afterwards,some models for both light and heavy HC feedstocks were derived by Plehierset al.[24] and Reynierset al.[25].

In a pyrolysis process,an important parameter which can be resulted from the simulations is the exact residence time.On one hand at a very low residence time,the secondary reactions which are very important for production of the light olefins,cannot be propagated and the yield of the main products is not satisfactory.On the other hand,the high residence time leads to the increase of the coke formation or cracking of the heavier molecules to the lighter products.Therefore,a compromise should be made between these two contradicting features.Enniset al.[26] suggested a residence time of lower than 0.1 s for a millisecond naphtha pyrolysis.Heynderickxet al.[8] reported a residence time of 0.45 s for a 4/4/2/1/1/1 split naphtha cracker with 47 m coil length.

Although the rapid methods for prediction of product distributions like machine learning techniques are nowadays welcomed[20],full 3-D CFD simulation provides comprehensive information for engineers and designers[27],despite the fact that it is computationally too expensive.As indicated in the above-mentioned studies,the coil is considered as one-dimensional plug flow reactor,except for the study carried out by Huet al.[15] in which 3-D equations were used inside the SL-II furnace.In this work,the coupled reactor/furnace performance of an industrial SRT-III naphtha cracking furnace with a capacity of 120 kt?a-1was simulated by a 3-D CFD model for the first time.To the best of our knowledge,such a type of furnace has never been studied.The results obtained from this coupled simulation provides better insight into the behavior of this industrial steam cracking furnace and can be used to improve efficiency for the unit.A comprehensive data gathering,measurements,and analysis were accomplished for both firebox and reactor domains.Over 200 variables and parameters were studied during the first days of each furnace run-length,including several analyses of compositionsviaonline and offline gas chromatography (GC).The products yields,the amount of produced coke,and the residence time were investigated as the main variables.

2.Numerical Simulation Procedure

2.1.Governing equations

2.1.1.Flow model

Having robustness,accuracy,and reasonable economy,the Reynolds-averaged Navier–Stokes (RANS) equations were used to model a compressible,reactive,single-phase,turbulent flow in steady-state condition.The governing equations consist of the momentum,heat,and mass transfers,and two transport equations for turbulent kinetic energy(TKE),k,and its dissipation rate(DR),ε,(thek-ε turbulence model) are described as follows:

Mass Balance (continuity equation):

Momentum balance:

Energy balance:

Species transport equation:

k-equation:

ε-equation:

where σh(=Prt/Cp=1?0),σY(=ρSct=1?0),σk(=1.0),and σε(=1.3)are the Prandtl number of enthalpy equation,the Schmidt number,and the turbulent Prandtl numbers,respectively.C1ε(=1.44),andCμ(=0.09) in the above-mentioned expressions are thek-εmodel constants.SEdenotes the heat generation/consumption in the reaction which is radiated in the firebox and convected in the reactor coils.The perturbation number for the input flow is equal to 17%,which is a relatively large number.The ideal-gas mixing law is used for calculation of ρ(mixing density),CP(heat capacity),λ (thermal conductivity),and μ (viscosity).Also,D(mass diffusivity) is calculated from kinetic theory.Considering the adiabatic straight entrance length for the reactor coil,fully developed gas flow in radiation section is guaranteed.

2.1.2.Combustion model

In order to take into account of the turbulence–chemistry interactions,an eddy dissipation concept (EDC) model was adopted.Based on this model,it is assumed that reactions occur in fine scales or small turbulent structures.The source terms in species equation considered as reaction rates are modeled as:

where ρ denotes mean density.Yiandrepresent the mean and the fine-scale mass fraction of speciesi,respectively.Xis the reacting fraction of the fine scales and is assumed to be unity.ξ*and τ*are the length fraction and residence time of the fine scales,respectively,and can be calculated from the following equations using Kolmogorov scales theory [11]:

whereCξ(=2.1377) andCτ(=0.4083) are volume fraction constant and time scale constant,respectively.According to Zhanget al.[11] this value for the time scale constant works well with most interior reacting flows utilizing standardk-ε model and EDC model.

B?senhoferet al.[28]compared several EDC models with different conversion ratio of the mass exchange rates and showed that the plug flow reactors (PFRs) fine structure treatment with the EDC,can give reasonable results for classical combustion conditions.

The used fuel gas is consisted of CH4/C2H4/C2H6/C3H8/H2with compositions of 85.52%,0.85%,4.94%,1.48%,and 7.21% (mol),respectively.It is burned in up to 10% (vol) of excess air.The C2H4content in the fuel gas was negligible and its heating value is close toTherefore,due to simplicity,ethylene was not explicitly included and methane is considered instead,with a concentration of 86.37% (mol).The side effect of this assumption on the fuel gas combustion is not considerable.The reaction scheme derived by Dryer and Glassman and Westbrook and Dryer is used[29,30]:

The combustion of H2is very fast.Regular measurements of the fuel and exhaust gas flow rate,temperature,pressure and composition is done during two run-length of the furnace.

The coupled ordinary differential equations of rate equations are stiff and conventional numerical integration is computationally expensive in CFD.Pope introduced theIn-situadaptive tabulation(ISAT)algorithm for numerically integration of reaction rate equations which resulted in 1000 times faster calculations compared to the conventional numerical integration referred to as direct integration [31].Kumar and Mazumder [32] presented an adoptive ISAT procedure in which a remarkable reduction in computational costs was achieved by two orders of magnitude acceleration in chemistry calculations.Therefore,in order to evaluateSY,the Arrhenius rates of reaction were integrated numerically using the ISAT approach.

2.1.3.Radiation model

The radiation is the main heat transfer mechanism in the furnace.In this study,according to recommendation given by Huet al.for the use of the discrete ordinate method(DOM)in parallel computing,this method was adapted to solve the radiative transfer equation (RTE) [12,21,33]:

The density of the flue gas is calculated in each cell from the local temperature,pressure,and mass fraction of species with the assumption of ideal-gas mixture [34]:

Temperature and concentration dependence of other thermophysical properties of the gas mixture are expressed as follows:

in which μiand λiare the viscosity and thermal conductivity of the pure componenti,respectively.φijdenotes the binary interaction parameter and is calculated from:

in which ‘‘f”can be either the viscosity (by Wilke method) or the thermal conductivity(by Mason and Saxena method)[35].The viscosity and thermal conductivity of the pure low pressure gas components were calculated based on method of Chunget al.[35].

Finally,the heat capacity of the mixture is computed based on a mass fraction average of the heat capacities of pure species as:

whereAi,Bi,Ci,andDiare the constant coefficients for the pure speciesi[36].

Due to the importance of the radiative heat transfer of flue gas in the furnace,the weighted-sum-of-gray gases model (WSGGM)was applied to calculate the absorption coefficient of the flue gas.In this method,the emissivity of the real gas is defined as the weighted sum of the emissivity coefficients of a number of gray gases and a transparent gas as follows:

According to the technical datasheets,the emissivity coefficients of tube skin and firebox walls are 0.85 and 0.75,respectively.

2.1.4.Reactor model

A full CFD simulation of a reactor coil with detailed kinetics is still a big challenge for the simulators[2,37].In this study,a molecular reaction set was used and a full CFD simulation is carried out for investigation of flow patterns and material compositions in the coils.Similar to firebox,a steady-state,compressible,reactive,single-phase,turbulent fluid flow based on the RANS equations with the standard k-ε turbulence model was used for simulation of the fluid flow in the reactor coils.

2.1.4.1.Feedstock and kinetic reactions.Naphtha is composed of several ingredients mainly paraffins.For the industrial furnace under study,the detailed composition of the feed entering the reactor was determinedviaGC mass analytical technique.Some other characteristics such as the specific gravity,the average molecular weight of the mixture,the PIONA weight fractions,and the boiling point distillation curve were also measured and tabulated in Table 1.Since the cracking of naphtha in the reactor coils is carried out at high temperature and low pressure,the process gas can be considered as an ideal gas mixture.Therefore,the physical properties of the process gas were calculated using Eqs.(27)–(31) [15].

In order to study of replacing the current feed with a little heavier one,having the indicated properties of last column in Table 1,the validated model is used to predict the product distribution.

In the current simulation,the molecular reaction scheme for cracking of naphtha was based on the experimental work of Kumar and Kunzru[5].The naphtha used in the current study is different from that used in Kumar and Kunzru’s study [5].As suggested by them,the initial selectivities should be adjusted for other naphthato achieve agreement between the model and measurements.Thus,the stoichiometry coefficient of primary reaction was tuned regarding the industrial measurements.The reactions and their parameter of the rate constants are listed in Table 2.Although Kumar and Kunzru did not report the parameters of rate coefficients for the reverse reactions (20)–(23),some studies did[38,39].Due to similarities,the parameters of rate coefficients for the reversible reactions and the reaction (24) were obtained from study of Seifzadeh Haghighiet al.[38].It is worth mentioning that,these reactions have very low conversion in the reverse path due to the low rate coefficients.

Table 1 The properties of feedstock used/planned for the industrial plant

2.1.4.2.Coke formation.As an undesired byproduct of thermal cracking processes of HCs,coke causes several severe problems such as conductive heat transfer resistance over the coil body,making hot spots and rapture on the coils,decreasing reactor diameter,volume,and therefore capacity,residence time and,consequently,changing the product distribution.Hence,the rate of coke formation and coke concentration along the reactor should be calculated.The coke formation precursors and reactions are determined based on the study of Kumar and Kunzru [23].These reactions and their rate expressions and parameters are described as(R24)–(R28)in Table 2.The reaction(R24)is considered regarding studies of Sundaramet al.[40]and Kumar and Kunzru[23]for pyrolytic coke formation.The density,thermal conductivity,and specific heat of the coke were considered to be 1600 kg?m-3,5.5 W?m-1?K-1,and 720 J?kg-1?K-1,respectively.

In the cracking process,the steam not only reduce the partial pressure of the HC feedstock,but also led to partial removal of the coke laid on the coils wall.This phenomenon takes place in the form of the following reactions and CO content on the coil outlet can be considered as a removal index.The reactions (R29) in Table 2 presents the removal effect of the steam and carbon dioxide.Thus,in the current simulation,29 reactions have been considered over 22 major species for the reactor side.

2.2.Boundary conditions

A coupled velocity inlet/pressure outlet condition was applied for firebox.The furnace works in natural draft conditions at the outlet pressure of -60 Pa gage.The flue gas composition,flowrate/velocity and temperature are presented in Table 3.

In the reactor side,mass flow inlet boundary was used for the gas enters radiation tubes reactor and all variables were prescribed at the inlet.The pressure is 80 kPa gage at the exit.The other conditions are described in the following section (Table 3).

2.3.Geometry characteristics and data gathering

2.3.1.Process conditions and reactor configuration

After preheating the feed up to the temperature of starting chain reaction (about 900 K) in the convection section,it enters the radiation section (or firebox).Fig.1(a) shows the scheme of the cracking furnace.There are six reactor coils in the simulated furnace which are configured on SRT-III type arrangement.This coil configuration and flow paths is demonstrated in Fig.1(b) and the details of the reactor configurations are previously summarized in Table 3.There are 12 long flame burners at the floor (6 burners in two columns) and 72 radiation burners positioned in four rows on both sides of the furnace wall.As convection section provides a fully developed profiles for velocity and turbulence parameters at the reactor inlet,an adiabatic entrance zone was considered at upstream of the reactor.In fact,the isolated tubes passing the furnace walls between convection and radiation section indicate the adiabatic zone.Due to symmetrical feature of the furnace,a geometry with 1 reactor coil was selected for the reactor side simulation.

Table 2 The reactor side reactions and their rate constants [5]

2.3.2.Industrial sampling,analysis,and data gathering

In order to measure,analyze,and process the industrial data,which contained over 200 variables,a regular sampling and recording from both firebox and reactor domains were conducted during several cycles of furnace run-length.Five sources of data have been considered:

1.Sampling and analysis of fuel gas,liquid feed,and process gas by an offline Agilent?3800 Varian GC.All samplings were performed regularly during the furnace run-length.For the liquid feed sampling,an ice-jacket bottle was utilized to avoid the outflow of lighter materials.Furthermore,since the temperature of cracking products leaving the transfer line heat exchanger(TLE)is still above 550 K,using a water-cooled sample cooler was necessary.Also,in order to separate the condensates from the gas sample,an empty impinger was used before collecting the process gas into the sampling bomb or ball.

2.Analysis of the flue gas by an online portable analyzer of Testo?350 M/Xl.The compositions of the resulted O2,H2,NO,CO,CO2,SO2,NO2,NOxand CxHyfrom the simulation is then compared to the industrial stack.

3.Recording the operating conditions from both of the distributed control system (DCS) and local indicators.

4.Measurement of coil surface temperature and furnace wall temperatures by a calibrated pyrometer of Raytek?Raynger 3i Plus in the 40 points of each coil and the facing furnace wall.

5.The coke concentration in the operation was determined by measuring the total amount of produced coke in one furnace run-length.Two forms of carried and deposited coke were identified.One of those are particles carried by the process gas and finally collected from the bottom of the fractionator during a normal operation.The second one,deposited on the inner reactor coil skin and TLEs and were drawn out during the first decoking process.The corresponding measurements based on different deposition sources are described as below:

(1) A sampling and analysis of flue gas during decoking process was performed.The weighted average concentration of CO2was calculated.

(2) Spalled coke layers moving to the next process unit operation and collected in the coke pot was weighed.

(3) The water-jet washed coke from the TLEs was weighed.

2.4.Computational grids

Each side were discretized using tetrahedral cells for near burners,fittings,curvatures,and etc.and hexahedral for rest of thedomains.A mesh independency study was carried out for each domain and the adopted mesh was verified.The details of the adoption are presented in Table 4.In this study,the minimum Yplus (Y+) on the walls was considered as discretization criteria for the flow field.An adopted mesh with Y+value greater than 300 is not appropriate for standard wall functions [41].Therefore,‘‘Large scale”has been rejected.Due to lower relative error in product yields,the third adapted mesh,denoted as ‘‘Adapted 2”,was selected.Although finer grids of‘‘Fine scale”were capable of more accurate calculations with lower relative errors,the computational costs drastically increased.Finally,the computational domains contained 2,858,771 and 2,298,265 volumes for firebox and reactor coil,respectively.The adopted mesh had skewness and orthogonality of 0.893 and 0.121,respectively.

Table 3 Furnace configuration and operating conditions

Table 4 The grid refinement for reactor coil simulation

2.5.Solution procedure for coupled furnace-reactor simulation

In order to simulate the whole furnace,the heat flux profile on the coil walls as the common boundary condition between the two flow domains is required.Therefore,a tube wall temperature and heat flux profile were pre-defined by a user-defined functions(UDFs) according to the industrial measurements and the simple 2-D or 3-D simulations,respectively.An iterative algorithm was used to simulate the coupled domains and determine the exact heat flux profile.In this algorithm,a firebox simulation is conducted for updating the tube wall temperatures.Then,a coil simulation is carried out to determine the heat fluxes.This path is repeated until the profiles have been converged and the results were comparable to the industrial data.The algorithm and sequence of solving the governing equations is shown in Fig.2.Each equation was solved using the previously updated values of the calculated variables.A second-order upwind scheme is used for discretizing the non-linear equations.The linearized set of equations were solved numerically based on finite volume discretization technique in Fluent?software.The Semi-Implicit Method for Pressure Linked Equations (SIMPLE) algorithm was applied to couple the momentum and continuity equations.The under-relaxation values for pressure,momentum,energy,and species material balance equations are 0.3,0.6,0.9 and 0.9,respectively.The converged simulation results were compared to relevant industrial data.

Fig.1.(a) The steam cracking furnace configuration (b) the SRT-III type of reactor coils.

Fig.2.The proposed algorithm for solving the model governing equations.

2.6.The model summarization

The presented model for the both domains of the radiation section,i.e.reactor coil and firebox,consists of several main equations considering the momentum,heat and mass transfers,and also rate equations for the reactions.In this work,the total number of partial differential equations(PDEs)that are simultaneously solved for the reactor and firebox sides are 28 and 17,respectively.Hereby,for each side of the furnace,the continuity and momentum equations give the pressure and velocities of each computational cell (4 PDEs).The TKE and its dissipation rate are calculated fromk-ε equations (2 PDEs).The energy equation gives the temperature(1 PDE).Also,some PDEs are required for solving the radiative heat transfer in firebox.There are 22 species for the reactor and 8 species for the firebox side.The concentrations of (N-1) species are determined by the (N-1) PDE equations utilizing ordinary differential equations (ODEs) of reaction rates as the source terms,noting that there are 29 ODEs for the reactor and 5 for the firebox.Thenth component concentration is calculated from following algebraic equation:

Also,there is a solid-phase energy equation:

The iterative algorithm of Fig.2 was applied to solve the set of equations of each section with conjugating boundary conditions until being converged.Depending on the initial estimates and number of iterations required,the computations for a coupled simulation took an average of 180 hours CPU time on a computer with 128 computational nodes in parallel.

The results are validated using percentage of absolute relative difference (ARD) defined as follows:

where sim.denotes simulation results and ind.indicates the industrial plant data.

3.Results and Discussion

3.1.Firebox simulation

Although the amount of excess air can be determined by the analysis of stack gas,determination of the proper/required amount of excess air for the industrial furnace is essential.Also,it can be used as a validation criterion for the model.For this purpose,as a part of the present study,the air was mixed with the fuel gas in several excess air levels up to 30 % (vol).Effects of these mixtures on furnace conditions have been investigated at three %(volume fraction) of excess air,i.e.0,10 %,and 30 %(vol) in order to make a comparison with industrial data.Fig.3 shows the temperature profiles when the furnace operated with different amounts of excess air.These temperature profiles are depicted at a typical width of the furnace (z=1280 mm) including a column of wall burners.The flue gas composition is shown in Table 5.It is obvious that the present of excess air guarantees the complete combustion and also increases NO production.Also,it has a cooling effect which can decrease the firebox temperature due to the lower temperature of the air.The best match between the industrial and simulated data in terms of both temperature profile and stack gas composition including O2,H2,NO,CO,CO2,and CxHyis obtained at the 10%excess air.On one hand,trying an excess air lower than the required or proper value,for example equal to stoichiometric ratio,caused the temperature of the combustion gases to become higher than the practice.On the other hand,it can be observed from Fig.3(a) that limited entry of the air into the furnace may results in incomplete combustion of the fuel,but there would be any coolant effect of the air and consequently an increment in the maximum temperature of the combustion gases and radiated energy.Hence,10% excess air is the better value for this furnace.At 10% excess air condition,almost uniform temperature profile was obtained due to the lower temperature of the unreacted air(Fig.3(b)).

Table 5 Industrial composition of the flue gas in comparison with the simulation results at different excess air amounts

Table 6 Comparison of simulation resulting from original and modified kinetic model with industrial and design data

Fig.4 shows the temperature contour around the floor burners.As can be seen,the temperature is lower around the burner inlet due to the lower fuel inlet temperature.Hence,it is observed that the maximum temperature occurs almost near the burners.Fig.5 shows the velocity and pressure contours around the wall burners.From the velocity contour,it can deduce that the gas velocity has considerably increased in the upper half of the burner due to the accumulation of combustion gases,and this trend continues until the products outlet.Similar to the study of Huet al.[2],two recirculation zones are seen in the firebox.These circulations cause enhancement in convective heat transfer due to the increase in flue gas residence time.The first zone occurs in the floor towards the middle of the firebox due to the high-velocity flow of the longflame burners where high-temperature flue gas returns to the bottom of reactor tubes (Fig.4).The second one is resulted from the flow of flue gas towards the convection section which affects both temperature and pressure profiles of the exhaust gas (Fig.5(a)).Although the mean velocity is less than 20 m?s-1in the firebox,the velocity rises to 50 m?s-1at the outlet.The pressure contour shows a 10 Pa pressure drop along the furnace which causes the natural draft.

Fig.3.The temperature (in K) profiles at z=1280 mm (a) no excess air,(b) 10 %(vol) excess air,and (c) 30 %(vol) excess air.

Fig.4.Temperature contour around the floor burners.

Fig.5.The (a) velocity and (b) pressure contours of the furnace around the wall burners.

The temperature profile on the furnace wall is shown in Fig.6(a).The maximum temperature is observed at the height of 8000 to 9300 mm in the middle of the radiation burners.This region is also affected by the high temperature combustion gasses inside the lower recirculation area of the firebox.The maximum temperature observed in this region is therefore caused by the conjugate effect of long-flame and radiation burner.Fig.6(a)shows the counters of simulation results for the furnace wall.The industrial measurements of the wall temperatures at different heights (h) and lengths(z)of the firebox and the corresponding simulation results are compared in Fig.6(b).This comparison shows that the simulation can well predict the practical temperature profiles with the relative error lower than 2%.

3.2.Reactor simulation

Fig.6.The temperature profiles on the furnace wall,(a) the simulation results and (b) the simulation results versus industrial measurements in the selected points.

The coil effluent which contains a variety of materials,was impinged and the dry gas was analyzed by GC.A comparison between design data,industrial data,and the simulation results is shown in Table 6.Considering suitable threshold,the simulation results are almost in a good agreement with the industrial data and shown to be reliable,especially the propylene-to-ethylene ratio (P/E) as a key factor in the cracking furnaces,is found to be nearly the same for the industrial data and simulation results.Hence,the temperature,pressure,velocity,and concentration distributions and the other variables in the reactor coil are investigated using simulation results.

The modification of the stoichiometric coefficients in the primary reaction caused the average of ARDs of product yields to decrease from 38.6%to 12.1%,which is mainly due to the deviation of the remaining propane.

It should be noted that,there is a bunch of uncertainties that makes hard the exact validation of modelling and simulation results against the designed and controlled operations in the labs.In the current work,although many accurate methods were applied for sampling,analysis,and measurement,deviations from the industrial data were observed in some cases which can be attributed to two main factors;the first was the instant variations in real conditions of the unit operation,especially in flow rates and compositions.For example,non-uniform splitting the feed into the four inlet tubes,or flame variations due to fuel gas perturbation or mixing with the air.The second was some physical deviations from the model assumptions,such as tube arching,burner failure or malfunction,which may customarily happen in an industrial unit.

There are several Y-and U-shaped bends in the process gas path which causes occurrence of radial pressure gradient inside the coil,especially at the exit of these bends.As a result,the velocity profile is affected and some regions of high local velocities are emerged.One of such regions is indicated by a dashed oval in Fig.7.The figure shows the velocity of the fluid in the longitudinal section of the coil,in which the mentioned effects are seen.The gas velocity through the coil depends on the density of the fluid,which itself is a function of fluid temperature,pressure,and the number of the produced molecules (exentropicity).As the reactions proceed and the number of moles increases,the velocity increases.Fig.8 depicts the pressure,density and velocity of the process gas along the coil.The inlet flow was split into four tubes and flows with the velocity of about 100 m?s-1.In pass #1,due to the sharp decrease in molecular weight and density of the ideal gas which is caused by cracking of heavy molecules in the primary reaction,the gas velocity increases sharply.According to the continuity equation of compressible fluid,density change is the only factor affecting the velocity of the process gas in each pass.In passes #3 and #4,the progress of secondary reactions increases the production of heavy aromatics while increasing the production of lighter molecules such as ethylene and methane.As a result,the change in the average molecular weight will be negligible.On the other hand,temperature and pressure of the gas remain relatively constant during the last pass,which leaves the density almost unchanged.Consequently,the average velocity of the gas remains almost constant.The gas is discharged from the coil at an average velocity of 240 m?s-1.The CFD simulation showed the residence time is 0.29 s with an average gas velocity of 220 m?s-1.

Fig.7.The velocity magnitude of the coil (m?s-1).

Fig.8.Variation in pressure,density,and average velocity of the process gas along the coil.

Fig.9.The dry weight based yields of the products and the coke concentration along the coil length.

Fig.9 shows variations of different products yields along the coil.The graph shows that ethylene yield increases with increasing coil length,due to increase of the residence time.But the propylene production shows a decreasing trend after experiencing a maximum of 23% at the middle of the fourth pass.This reduction can be attributed to the higher conversion of the secondary reactions and further cracking to the lighter hydrocarbons.The trend of the other olefin products,i.e.butylene(s) and butadiene(s),is similar to propylene but deceasing trend for them occurred nearly at the coil exit.Most of the heavier HCs are approximately cracked before reaching this region,but the lighter paraffins still can be dehydrogenated to olefins and,simultaneously,ethane yield starts to decrease slightly at the coil exit.These trends emphasize the importance of the residence time on the P/E ratio and as well as the heavier olefins production.As shown in the figure,the two curves of propylene and ethylene are diverging from each other along the coil length.This can be interpreted as a decrease in the P/E ratio as a result of an increase in the residence time.As the temperature and concentration of precursors increase in the last pass,the coke production intensity increases and the coke concentration increases rapidly,reaching up to 0.0019% (mass) at the reactor outlet.

The effect of high local velocity of the fluid on the rate of reactions and concentration of products was studied by a plotting radial changes.The selectivity of propylene,as a major product of naphtha cracking process,highly depends on the residence time.The changes of the composition of propylene is shown in Fig.10 for pass #4 of the coil.As can be seen from this Fig.,more details of spatial level of the principal product could be provided by the CFD simulation.Under normal operating conditions,the possibility propylene decomposition to acetylene,methane,and ethylene is higher than propylene production.The secondary reactions of 11 and 16 with the fairly high pre-exponential factors of 7.284×1012and 7.638×1012,respectively,are in charge of the propylene consumption reaction.

The dashed and solid lines in Fig.10 represent the composition of propylene in two sides of the coil.Comparing the dashed regions in Fig.s 7 and 10 reveals that the composition of propylene is affected by the velocity profiles in the side whose velocity is higher.Therefore,the residence time is lower so that reactions have less time to progress rather than the other side.

Fig.10.The radial changes of the propylene mass fraction along the pass #4.

Fig.11.An illustration of the process gas and tube skin temperature profile.

The coking concentration is shown in Table 6.The simulation results demonstrated that the average concentration of coke at the coil exit is 0.0019%(mass)on dry basis.However,according to the aforementioned method in Section 2.3.2,the coke weight in the operation was determined by measuring the total amount of the produced coke in one furnace run-length.The coke concentration in the flow was estimated to be about 0.0022% (mass).Also,of the total produced coke,87%is found to deposit on the inner tube skin.The remaining coke was carried by the process gas,or left the coil in the form of spalling.The main reason for the difference between the simulation results and the industrial data can be attributed to the effect of di-methyl di-sulfide (DMDS) and other sulfuric components such as H2S and SO2existing in the feed,but have been neglected in the simulation.Some studies alleged that such components can affect the catalytic role of Ni on the coke formation reactions by producing sulfuric radicals like HS’[42].However,Wanget al.[43] explained that addition of more than 0.01%(mass) DMDS to feed can increase the chance of coke formation.In the current industrial plant,the feed included 0.01% to 0.05%(mass)DMDS.Therefore,the simulation results in the present study are more consistent with the hypothesis posed by Wanget al.

Fig.11 shows the temperature profile of external tube skin as well as the process gas temperature along the tube length.The fluid temperature curve indicates that the fluid enters the coil at a temperature of approximately 900 K.Despite the endothermic nature of the cracking reactions,the process gas temperature increases gradually according to the large heat input.The temperature of downward flow of the first and third passes begin to decrease at the end of these passes in which heat input decreases but the reactions progress.The tube skin temperature is affected by both inside and outside conditions and has a similar trend to the facing wall.Also,the resulted tube skin temperatures meet the measured industrial data.The difference between the temperature of tube skin and process gas shown in Fig.11,indicates the net heat flux radiated (or convected) to the process fluid.The maximum temperature difference is located at the middle of each pass,and it has the highest level at the first pass.

4.Conclusions

Thermal and hydrodynamic behavior of an industrial naphtha cracking SRT-III furnace was simulated using full CFD considering both of the firebox and reactor sides.The furnace under study was equipped with both radiation and long-flame burners.An algorithm was designed for solving a set of equations by aids of a coupled CFD simulation.The resulting product distributions agreed well with measured industrial yields,especially for the main productsi.e.ethylene and propylene.The velocity,temperature and concentration profiles were also investigated using the simulation.According to the simulation results,10 %(vol) excess air showed the best agreement with measurements,which guaranteed a complete combustion,a small amount of NO production,and also provided required temperature for the fuel combustion.The process gas velocity profile demonstrated that a sharp reduction of velocity occurs up to its half pass and then becomes nearly constant.The simulation resulted an average gas velocity of 220 m?s-1which led to the residence time of 0.29 s.Such a coupled full CFD simulation provides a foundation for achieving higher mass production and also producing heavier olefins,especially propylene,by manipulating the conditions of the naphtha furnaces.In spite of a timeconsuming computation,the consistency between the simulation results and plant experimental data is promising and shows that such a good accuracy can only be handled by a coupled full CFD simulation.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

The authors greatly acknowledge the support of Bandar-e-Imam petrochemical company(BIPC),Iran,which help the authors provide the industrial data and also their technical support to this research.

Nomenclature

Ai,Bi,Ci,Dithe coefficients of the heat capacity of pure speciesi

bεi,jemissivity gas temperature polynomial coefficients

Cimolar concentration of species i,mol?m-3

Cpiheat capacity of componenti,J?kg-1?K-1

Cξvolume fraction constant (=2.1377)

Cτ time scale constant (=0.4083)

Dmolecular diffusivity,m2?s-1

Etotal energy per unit mass,J?kg-1

Eaactivation energy,J?mol-1

Gkgeneration of turbulent kinetic energy,J?m-3?s-1

generation of dissipation of turbulent kinetic energy,J?m-3?s-1

ggravitational constant (=0.9806 m?s-2)

hmean static enthalpy per unit mass,J?kg-1

hienthalpy of formation for componentiper unit mass,J?kg-1

Iradiation intensity,J?m2?s-1

Ji,jdiffusion flux of speciesjin speciesi,kg?m2?s-1

kturbulent kinetic energy per unit mass,J?kg-1

kiArrhenius parameters for reaction rate coefficient in reactioni,s-1

krrate constant for reactionr,s-1

k0rate constant for reaction,s-1

Mimolecular weight,g?mol-1

Nnumber of chemical species in the system

nrefractive index

Prtturbulent Prandtl number

ppressure,Pa

Rideal gas constant (=8.314 J?mol-1?K-1)

rrivolumetric rate of reaction i,mol?m-3?s-1

r position vector

SEsource term in the energy equation,J?m-3?s-1

Sksource term in the turbulent kinetic energy equation,J?m-3?s-1

SMsource term in the momentum equation,kg?m-2?s-2

SYsource term in the species transport equation,J?m-3?s-1

source term in the dissipation of turbulent kinetic energy equation,J?m-3?s-1

Sctturbulent Schmidt number

s direction vector

s′scattering direction vector

Ttemperature,K

uvelocity component in the spatial directions,m?s-1

Xreacting fraction of the fine scales (=1)

xcoordinate direction,m

Ymass fraction of species

Yimean mass fraction on of speciesi

fine-scale mass fraction of speciesi

α absorption coefficient,m-1

emissivity weighting factors for theith fictitious gray gas

βi,rstoichiometric coefficient of speciesiin reactionr

δijKronecker delta

ε dissipation rate of turbulent kinetic energy,m2?s-3

εgemissivity of the real gas

κ absorption coefficient

λ thermal conductivity,W?m-1?K-1

μ molecular viscosity,kg?m-1?s-1

μtturbulent viscosity,kg?m-1?s-1

υ kinematic viscosity,m2?s-2

ξ*length fraction of the fine scales

ρ density,kg?m-3

σ Stefan–Boltzmann constant (=5.672×10-8J?s-1?m-2?K-4)

σhPrandtl number of enthalpy equation

σkturbulent Prandtl number

σsscattering coefficient,m-1

σYSchmidt number

σεturbulent Prandtl number

τ*residence time of the fine scales,s

φ phase function

φijinteraction parameter

Ω′solid angle,degrees

Subscripts

E energy

g gas

h thermal

isymbol for spatial coordinates,symbol for pure species

jsymbol for spatial coordinates

ksymbol for spatial coordinates

M momentum

msymbol for pure species

rreaction

ref reference

s solid

t turbulent