e-mail: Fabian.Diaz@ecopetrol.com.co Jesus.Castro@ecopetrol.com.co
]]> (Received, April 15, 2009; Accepted May 26, 2010)* To whom correspondence may be addressed
ABSTRACT
A mathematical model for simulation of refinery furnaces is proposed. It consists of two different submodels, one for the process side and another for the flue gas side. The process side is appropriately modeled as a plug flow due to the high velocity of the fluid inside the tubes. The flue gas side is composed by a radiative chamber and a convective section both connected by a shield tube zone. Both models are connected by the tube surface temperature. As the flue gas side model uses this temperature as input data, the process side model recalculates this temperature. The procedure is executed until certain tolerance is achieved. This mathematical model has proved to be a useful tool for furnace analysis and simulation.
Palabras clave: furnace, simulation, mathematical models, oil refinery.
RESUMEN
En este trabajo se presenta el desarrollo de un modelo matemático para simulación de hornos de refinería, el cual consiste en dos sub-modelos diferenciados, uno para simular el lado proceso y otro para simular el gas de combustión. El lado proceso es apropiadamente modelado con un perfil de velocidad plano debido a la alta velocidad del fluido dentro de los tubos. El lado gas de combustión está compuesto por una cámara de radiación y una sección de convección, ambas unidas por una zona de tubos de choque. Los dos sub-modelos interactúan a través de la temperatura de superficie de los tubos, siendo esta un dato de entrada al sub-modelo del lado gas de combustión y es re-calculada por el sub-modelo del lado proceso. Este procedimiento es ejecutado en un ciclo iterativo hasta que cierta tolerancia es alcanzada. Este modelo matemático ha demostrado ser una herramienta muy útil para el análisis y simulación de hornos.
Key words: hornos, simulación, modelos matemáticos, refinería de petróleo.
INTRODUCTION
Furnaces for refinery services are considered a key piece of equipment as they provide the necessary heat to carry out distillation processes, cracking reactions and some other processes. They consume considerable amounts of energy and their cost usually ranges between 10% and 30% of the total investment. Given its importance and its cost inside the refination chain it's evidently the significance of a mathematical model for furnaces that can be implemented with the process simulators.
]]> Proper furnace simulation is necessary to analyze variations in regular process conditions such as changes in the feed charge, firing conditions, fuel composition, internal geometry, etc,. A simple mathematical model that takes into account these parameters is, therefore, a very useful simulation tool for refiners. Also considering the increasing necessity of processing heavier oils, this furnace simulator can predict the response of the equipment to changes in the design operating conditions.Many researchers have modeled furnaces using defined components in the feed especially for thermal cracking processes such as Niaei, Towfighi, Sadrameli & Karimzadeh (2004), Detemmerman & Froment (1998), Joo, Lee, Lee & Park (2000), Heynderickx, Oprins & Marin (2001) but few researchers have used pseudo-components in the feed which is necessary to simulate petroleum crude. Maciel & Sugaya (2001) used 16 pseudo-components to simulate thermal cracking inside refinery furnaces, however, in this work reactions are not considered. Commercial furnace simulators such as FRNC and HTRI Xfh use a separate source to create and calculate properties of pseudocomponents; the mathematical model proposed in this work does not need a separate source for that purpose.
DESCRIPTION OF THE MODEL
The process sub-model simulates in detail the behavior of the charge that flows through the furnace tubes. Physical properties, phase equilibrium, pressure drop and flow regimes are calculated for each integration step.
Petroleum feeds are characterized from the True Boiling Point (TBP) curve and the gravity curve. If the gravity curve is not available it can be extrapolated from the total gravity admitting the Watson characterization factor (KW) constant through the whole curve. The charge is divided into pseudo-components, with the gravity and the normal boiling point known their properties can be calculated using several correlations available in the literature. In this model API TDB (1997) and Aladwani & Riazi (2005) recommendations are used.
The flue gas side sub-model calculates the radiative and convective heat fluxes tube by tube. Tube surface temperatures and fin tip temperatures are calculated using API Standard 530 (2001) and ESCOA (2002) methods respectively. With the flue gas temperatures along the furnace the draft profile can be calculated.
The furnace model proposed manages to calculate the main operational variables for refinery services by combining successfully a process model and a flue gas side model. Data entry for the charge is very simple, only a distillation curve and a gravity curve (or total gravity) are needed. Pure light ends and steam can also be integrated to the charge in order to analyze the process behavior. This capability to analyze different charges that are typical in refinery services makes the model proposed a valuable tool for process engineers who can easily predict operational variables for the daily operation.
The model internally calculates, for every integration step, physical properties and phase equilibrium for petroleum fractions; also calculates heat transfer coefficients, holdup, flow regimes and pressure drop in the case of two phase flow. Pressure drop is an important parameter since it exponentially increases with the vapor fraction in the charge. All this calculations make the performance of the model superior when compared to other commercial software for furnace simulation.
MATHEMATICAL MODEL FOR FURNACE SIMULATION
The mathematical model proposed for furnace simulation consists of two different sub-models connected by the tube surface temperature and a pre-processor for petroleum feed characterization:
PETROLEUM FEED CHARACTERIZATION
A TBP curve and a density value must be provided to the simulator. In a first step, the TBP is completed from 0 to 100% of distillated volume using an interpolation polynomial. Then the extrapolated curve is uniformly divided in 30 cuts and the NBP of each cut is calculated using a composite Simpson's rule (See Figure 1). To calculate the density of the 30 pseudo-components, the Watson characterization factor is admitted constant through the whole curve. With the NBP and the SPGR known, molecular weight, critical properties and acentric factor are easily calculated with the correlations recommended by Aladwani & Riazi (2005) and the API TDB (1997), see Table 2.
Flue gas side sub-model
The Hottel one-gas-zone radiation method (Hottel, 1974) has proven to be an appropriate and effective model to simulate radiation chambers when highly detailed radiative conditions are not needed. The flue gas side model uses the mentioned radiation method to calculate the average and maximum heat flux to the tubes and the flue gas temperature. Hottel (1974) proposed that the bridgewall temperature can be considered as the following function:
Where QR is the heat supplied to the sink in the radiation chamber, QF is the supplied in fuel, TF is the so called adiabatic pseudo-flame temperature and d is a dimensionless constant. Based on the analysis of furnace simulation results and industrial data the value of 1,2 is recommended instead of the value of 1,33 suggested by Hottel (1974).
A different model is used in the shield tube zone, a model that includes radiation and forced convection in bare tube bundles. The radiative heat flux to the shield tubes is calculated using the recommendations published by Stehlík, Kohoutek & Jebácek (1996). The tube bundle is replaced by an equivalent plane surface since the first and second row of tubes intercept approximately all the radiation coming from the radiation chamber.
The convective heat fluxes in the bare and extended surface tube zones are calculated with the methods listed in Table 1.
]]> Pressure drop for the flue gas in convection sections is calculated using the ESCOA (2002) methods. API Standard 560 is used to calculate the draft profile.Figure 2 illustrates how radiative and convective heat fluxes are calculated from tube surface temperatures, flue gas temperatures and draft profile are also calculated in the flue gas side sub-model. Tube surface temperatures are re-calculated in the process sub-model using API Standard 530 (2001) methods.
Process sub-model
The process side is appropriately modeled as a plug flow due to the high velocity of the fluid inside the tubes (Maciel & Sugaya, 2001). Liquid or gaseous feeds can be used. Light oils vaporize as walk through the furnace tubes. In this case, two parallel plug flows (vapor and liquid) are modeled and several correlations regarding two-phase flow in tubes are used for a more realistic approach to the industrial process. For a detailed description of the process sub-model procedure see Diaz (2008) and Díaz, Wolf , Sugaya, Maciel & Costa (2008). The equations that describe the process sub-model are as follows:
Energy balance:
Pressure:
In the case of vertical flow, the contribution of the static head to the pressure drop is calculated using Equation 4; however, this static head is not totally recovered since vaporization takes place as flow proceeds downward through the furnace.
Phase equilibrium K-values are estimated using a modified Soave-Redlich-Kwong equation presented in the API TDB (1997); convergence is achieved by using a modified Rachford & Rice (1952) procedure. Since the equations of state are not appropriate to calculate the density of the liquid phase, different methods are used (See Table 3). To calculate the density of the vapor phase the compressibility factor found in the equilibrium calculations is used. Methods for physical properties are listed in Table 3:
If the equilibrium calculations determine the existence of a second phase, alternate procedures are activated to determine parameters for flow in two phases. Two-phase flow parameters are calculated using meth ods listed in Table 4. The whole procedure followed by the process sub-model is illustrated in Figure 3.
The flue gas side sub-model uses the tube surface temperature as input data, so the calculations start with an assumed tube surface temperature. Heat fluxes are calculated and then used in the process sub-model to heat up the feed charge. Since tube surface temperatures are highly dependent on bulk process conditions, they are recalculated using the API RP 530 methods. Results are compared and if certain tolerance is not achieved, the procedure starts all over again with the updated temperatures. The whole procedure is illustrated in Figure 4.
APPLICATION OF THE MODEL
The mathematical model for furnace simulation proposed has been tested with refinery plant data from the Barrancabermeja Refinery (Barrancabermeja, Colombia). The furnaces are typical box type configuration with burners located in the floor and refractory-backed radiant tubes in a single row (See Figure 5). The feed charges are crude oils in the range of 21,3 to 44,9 API gravity. The results of the simulations are summarized in Table 5. where good agreement between the calculated and the measured data is observed.
Results of the simulation of the furnace H2001 are analyzed to evaluate the furnace performance. Three different petroleum feeds are charged to the simulation as seen on Table 6. Tubes are numbered being 1 the entrance and 34 the exit of the feed.
Figures 6 and 7 show how the pressure drop increases exponentially with the vapor fraction. In the last four tubes of the furnace the tendency in the pressure drop changes because the tube diameter changes from 6 to 8 inches. It is also seen a higher pressure drop in the lighter feed since it produces more vapor.
Shield tubes (Tubes 13-16) receive an important heat flux load and consequently show high tube surface temperature (See Figure 8). However the maximum surface temperature appears in the last tubes where the feed is warmer and closer to the peak heat flux in the flame. The change in diameter from 6 to 8 inches decreases the heat transfer coefficient from
]]>
1700 to 1000 W/m2K (See Figure 10.) and the peak heat flux (See Figure 9.) make the tube 31 indicate the higher tube surface temperature. The heavier feed (F1) has the lower heat capacity and removes less heat from the flue gases. It is seen on Table 7, that F1 leaves the furnace 8K cooler than F3.
The higher tube surface temperatures observed in the simulations with the heavier feeds are caused mainly by the lower heat transfer coefficient (HTC) (Figure 10) since the heat fluxes remain practically unchanged for the three feeds (Figure 9). The difference among the HTC's in the feeds tends to decrease at the last tubes of the furnace because of two opposite effects: Higher temperature increases the HTC while higher vaporization decreases the HTC. Therefore, the heavier feed is cooler, less vaporized and the HTC tends to increase. HTC's were calculated with the API Standard 530 (2001) methods.
]]> CONCLUSIONS
ACKNOWLEDGMENTS
The authors want to thank the engineers Martha Yolanda Morales from Apoyo Técnico a la Producción (ATP) in the Barrancabermeja Refinery (Barrancabermeja, Colombia) and Jaqueline Saavedra Rueda from Instituto Colombiano de Petróleo (ICP) for their support in the acquisition of mechanical and processing data necessary to simulate the furnaces presented in this publication. Finally, the authors gratefully acknowledge Ecopetrol - ICP for the financial support that permitted the development of this work.
]]> REFERENCES
Aladwani, H. A., & Riazi, M. R. (2005). Some guidelines for choosing a characterization method for petroleum fractions in process simulators. Chem. Eng. Res. Des. 83, 160-162. [ Links ]
API Standard 530 (2001). Annex B, Calculation of maximum radiant section tube skin temperature, American Petroleum Institute. [ Links ]
API TDB (1997). Technical Data Book-Petroleum Refining. American Petroleum Institute, 6th edition. [ Links ]
Awad, M. M., & Muzychka Y. S. (2005). Bounds on two phase flow. Part 2. Void fraction on circular pipes, ASME Inter. Mechanical Eng. Congress and Exposition. Orlando, Fl., USA, 823-833. [ Links ]
Beattie, D. R., & Whalley, P. B. (1982). A simple two-phase frictional pressure drop calculation method. Int. J. Multiphase Flow, (8), 83-87. [ Links ]
Briggs, D. E., & Young, E. H. (1963). Convection heat transfer and pressure drop of air flowing across triangular pitch banks of finned tubes. Chem. Eng. Prog. Symp. Ser., 59: 1-10. [ Links ]
Bromley, L. A., & Wilke, C. R. (1951). Viscosity behavior of gases, Ind. Eng. Chem., 43, 1641-1648. [ Links ]
Chato, J. C., Yashar, D. A., Wilson, M. J., Kopke, H. R., Graham, D. M., & Newell, T. A. (2001). An Investigation of Refrigerant Void Fraction in Horizontal, Microfin Tubes. Inter. J. of HVAC& Research, 7(1), 67 - 82. [ Links ]
Coddington, P., & Macian, R. (2002). A study of the performance of void fraction correlations used in the context of drift-flux two-phase flow models. Nucl. Eng. Des., 149: 323 - 333. [ Links ]
Dean, D. E., & Stiel, L. I. (1965). The viscosity of nonpolar gas mixtures at moderate and high pressures. AIChE J., 11: 526-532. [ Links ]
Detemmerman, T., & Froment, F. (1998). Three dimensional coupled simulation of furnaces and reactor tubes for the thermal cracking of hydrocarbons. Revue de l'institut français du pétrole, 53 (2), 181 - 194. [ Links ]
Díaz, Fabian A. (2008). Desenvolvimento de Modelo Computacional para Craqueamento Térmico, M. Sc. Thesis dissertation, Universidade estadual de Campinas, UNICAMP, Campinas, Brazil, 138pp. [ Links ]
Díaz, Fabian A., Wolf Maciel M. R., Sugaya M. F., Maciel Filho R., & Costa C. B. B. (2008). Furnace modeling and simulation for the thermal cracking of heavy petroleum fractions. XVII Congresso brasileiro de engenharia quimica, COBEQ, Recife, Brazil. [ Links ]
ESCOA (2002). Engineering Manual, version 2.1. [ Links ]
García, F., García, J. M., García, R., & Joseph, D. D. (2007). Friction factor improved correlations for laminar and turbulent gas-liquid flow in horizontal pipelines. Int. J. Multiphase Flow, 33 (12), 1320-1336. [ Links ]
Gnielinski, V., Zukauskas, A., & Skrinska, A. (1983). Banks of plain and finned tubes. Heat Exchanger Design Handbook, Hemisphere, New York. [ Links ]
Grimson, E. D. (1938). Heat transfer and flow resistance of gases over tube banks. Transaction of ASME. (58), 381-392. [ Links ]
Hankinson, R., & Thomson, G. (1979). A new correlation for saturated densities of liquids and their mixtures. AIChE J. 25 (4), 653-663. [ Links ]
Heynderickx, Geraldine J., Oprins, Arno J. M., & and Marin, Guy B. (2001). Three-Dimensional Flow Patterns in Cracking Furnaces with Long-Flame Burners. AIChE J. 47 (2), 388 - 400. [ Links ]
Hottel, Hoyt C. (1974). First estimates of industrial furnace performance-the one-gas-zone model re-examined. Heat Transfer in Flames, 5-28. [ Links ]
Hughmark, G. A. (1962). Holdup in gas-liquid flow. Chem. Eng. Prog, 53 (4), 62 - 65. [ Links ]
Joo, E., Lee, K., Lee, M. & Park, S. (2000). CRACKER - A PC based simulator for Industrial cracking furnaces. Computers and Chemical Engineering. 24: 1523-1528. [ Links ]
Kendall, J. & Monroe, K. P. (1917). The viscosity of liquids. 2. The viscosity composition curve for ideal liquid mixtures. J. American Chemical Society, 39 (9), 1787-1802. [ Links ]
Khan, W. A., Culham, J.R. & Yovanovich, M. M. (2006). Convection heat transfer from tube banks in crossflow: Analytical approach. Int. J. of thermophysics and heat transfer, 49 (25-26), 4831-4838. [ Links ]
Kim, J. & Ghajar, A. J. (2006). A general heat transfer correlation for non-boiling gas-liquid flow with different flow patterns in horizontal pipes. Int. J. Multiphase Flow, 32, 447- 465. [ Links ]
Lee, B.I. & Kesler, M. G. (1975). A generalized thermodynamic correlation based on three-parameter corresponding states. AIChE J.,Vol. 21 (3), 510-527. [ Links ]
Maciel Filho, R. & Sugaya, M. F. (2001). A computer aided tool for heavy oil thermal cracking process simulation. Computers and Chemical Engineering. 25: 683 - 692. [ Links ]
Niaei, A., Towfighi, J., Sadrameli, S. M. & Karimzadeh, R. (2004) The combined simulation of heat transfer and pyrolysis reactions in industrial cracking furnaces. Applied Thermal Engineering. 24: 2251 - 2265. [ Links ]
Olujic, Z. (1985). Predicting two-phase-flow friction loss in horizontal pipes, Chem. Eng., 92 (13), 45 - 50. [ Links ]
Rachford, H. H. & Rice, J. D. (1952). Procedure for Use of Electrical Digital Computers in Calculating Flash Vaporization Hydrocarbon Equilibrium, J. of Petroleum Technology, 4, s.1, 19, s.2, 3. [ Links ]
Riazi, M. R. & Daubert, T. E. (1987). Characterization parameters for petroleum fractions. Ind. Eng. Chem. Res, 26 (4), 755-759. [ Links ]
Rouhani, Z. & Axelsson, E. (1970). Calculation of volume void fraction in the subcooled and quality region. Int. J. Heat Mass Transfer, 13, 383 - 393. [ Links ]
Spencer, C. F. & Danner, R. P.(1973). Prediction of Bubble-Point Density of Mixtures. J. Chem. Eng. Data, 18 (2): 230. [ Links ]
Stehlík, P., Kohoutek, J. & Jebácek, V. (1996) Simple mathematical model of furnaces and its possible applications. Comp. Chem. Eng., 20 (11), 1369 -1372. [ Links ]
Stiel, L. I., & Thodos, G. (1961). The viscosity of nonpolar gases at normal pressures. AIChE J., 7: 611- 615. [ Links ]
Twu, Chorng H. (1984). An internally consistent correlation for predicting the critical properties and molecular weights of petroleum and coal-tar liquids. Fluid phase equilibria, 16 (2), 137-150. [ Links ]
Twu, Chorng H. (1985). Internally consistent correlation for predicting liquid viscosities of petroleum fractions. Ind. Eng. Chem. Process. Des. Dev., 24: 1287-1293. [ Links ]
Weierman, C. (1976). Correlations ease the selection of finned tubes, Oil and Gas J., 74 (36), 94-100. [ Links ]
Woldesemayat, M. A. & Ghajar, A. J. (2007). Comparison of void fraction correlations for different flow patterns in horizontal and upward inclined pipes. Int. J. Multiphase Flow, 33, 347 - 370. [ Links ]
Zukauskas, A. (1972). Heat transfer from tubes in Cross Flow. Advances in heat transfer, 8, 93 - 160. [ Links ] ]]>