Modeling and Simulation of an Industrial Top-Fired Methane Steam Reforming Unit

,


INTRODUCTION
The mathematical modeling of a methane steam reforming unit (MSRU) is a challenging task.This process involves complex heat phenomena because radiation, convection, and conduction mechanisms are responsible for the simultaneous heat transfer between several control volumes within the equipment.In parallel, a network of chemical reactions creates radial and axial temperature and concentration gradients in the tube and its surroundings.The precise determination of all phenomena occurring in this system requires then the representation of highly complex and interconnected systems so that simplifications are usually made.Excessive simplification may, however, limit the adherence of the model to the actual system.There is no universal model applicable in all cases, but one might be the most appropriate to an application.Several studies about modeling and simulation of MSRUs are reported in the literature, using different approaches, from the simplest to the most complex, depending on the purpose for which the model was developed.Such models can be distinguished by the representation of the catalytic bed equation as well as by the treatment given to model the radiant furnace.
−7 Usually, the models are stationary, one-dimensional, and pseudohomogeneous.The heterogeneous model evaluates the temperature and concentration profiles between the gas and solid phases as well as the gradients inside the solid phase so that the model might present convergence problems due to the steep gradients in the catalyst particle. 8The pseudohomogeneous models, on the other hand, implicitly consider the transport limitations in the interior and on the surface of the catalyst particles in the kinetic expressions. 8−19 The precise quantification of the flaming radiation or hot gases is very difficult to obtain, mainly due to the complexity of this phenomenon, the nonlinear dependence on temperature, and the participation of a gaseous medium in the emission and absorption of radiation.Lobo and Evans 20 propose the first model for the calculation of radiant heat exchange in furnaces, which is the basis for most of the following improvements, such as the models from Roesler 14 and Hottel and Sarofim, 13 whose approaches are frequently used by many authors.Among these two approaches, Hottel and Sarofim 13 are the most used due to their simplicity when estimating the heat exchange and exhaust temperature of the flue gases in the furnace.
There are still some works that study both the catalytic tubular reactors and the radiant furnace, where the combustion reactions occur.−25 Most of these models are stationary and used to monitor, simulate, and optimize the temperature profiles along the length of the MSRU.These works consider that the combustible gas is burnt completely and isothermally at the entrance of the furnace at the feed temperature, assuming (i) a profile for the release of the combustion heat through the axial coordinate of the MSRU, 23−26 (ii) a profile for the temperature of the outer wall of the reforming tube, 2,27 or (iii) even a profile for the temperature of the furnace. 28The heat release profile is defined according to empirical correlations that generally have as input the flame length and the total heat of combustion, as suggested by Roesler 14 and Hyde et al. 29 Roesler 14 proposes a parabolic heat release profile, whose parameters are adjusted to industrial data. 24,25,30Hyde 29 represents the heat profile as an exponential function, which allows us to obtain a maximum flow of heat at one-third from the top of the furnace.These profiles are usually developed by extrapolating the scarcely available experimental data. 12,23,28None of these approaches considers the effects of the kinetics of the combustion reactions along the length of the furnace.
The combustion reactions greatly influence the composition and temperature profiles inside the furnace, reforming reactor, tube walls, and refractory.Therefore, a more precise representation of the combustion phenomena would be of great value for the process monitoring.Despite the extensive literature review, up to the author's knowledge, the models that attempt to describe the combustion phenomena in the MSRU are scarce.Lee et al. 31 and Pret et al. 5 consider the combustion reactions when modeling the steady-state reforming process and Tran et al. 19 when approaching the dynamic model.The three studies describe the combustion phenomenon based on the model developed by Magnussen and Hjertager, 32 which relates the rate of combustion to the rate of turbulent energy dissipation.The reaction rate is limited by the mixing effect between fuel, oxygen, and products, allowing the prediction of a slower temperature profile, as observed in the process operation.These authors, 5,19,31 however, propose a more complex model based on computational fluid dynamics (CFD), which might not be applied for control and optimization studies.Authors that proposed optimal control applications on MSRU 33,34 emphasized the computational difficulties in combining online realtime calculations with such CFD studies.Some simplifications are commonly adopted in the mathematical modeling by neglecting the following effects: the energy transferred by convection at the furnace, the radiation reflected by the reforming tubes, the participation of a gaseous medium in the process of radiation absorption, shading effects due to the presence of other tubes, and finally, the dependence of the gases properties on temperature and pressure. 1,2,10,21,22,25,35−38 Since temperature variations of up to 1000 K occur throughout the length of the MSRU, the physical properties of the gases can vary considerably.To ensure the accuracy of the model, therefore, it is essential to consider the temperature and pressure dependency when computing the physical properties of the gas mixture, such as viscosity, thermal conductivity, and heat capacity.
The most detailed models found in the literature basically use three control volumes (the CFT, tube wall, and furnace) for modeling of the MSRU.The refractory wall is usually neglected.Large temperature gradients, however, might compromise the structure of the refractory material due to thermal shock and cause energy losses to the surrounding environment.These losses imply a greater consumption of combustible gas and, therefore, a higher operational cost.The description of the heat transfer phenomena at the refractory wall within the furnace then would be of great value.The knowledge of temperature gradient allows taking corrective actions and then avoiding compromising the refractory material.Moreover, such a model could be useful for investigating the influence of the temperature gradient on the thermomechanical properties of the refractory material (mechanical strength, conductivity, and porosity) and hence on service life.
The main goal of this paper is to develop a model of an industrial MSRU, considering four control volumes, the CFT, tube wall, furnace, and refractory, using a more simplified approach than the one commonly reported in the literature, by the computational fluid dynamics approach.The mass, momentum, and energy balances, together with other constitutive relations, make up the phenomenological model, which is described by a set of differential algebraic equation (DAE).This stationary model will be used for practical purposes, helping operators evaluate and monitor key variables for the reforming process as temperature and concentration profiles along the four control volumes considered.In addition, it will serve as a basis for the development of a dynamic model for the purposes of process control and optimization in future studies.
The combustion reaction rates are described by a kinetic model so that the fuel gas is burned and its composition changes along the length of the furnace.Most studies in the literature consider that the fuel gas is burnt completely and isothermally at the entrance of the furnace.As a result, the usual literature approach assumes a flat flue gas composition profile.We explicitly consider the dependency of the gas physical properties on temperature and discuss the ideal gas assumption due to the high content of hydrogen in the gas mixture.The model further considers more detailed heat transfer effects as the convection at the furnace and the energy reflected by the refractory and by the tube wall.Last but not least, the refractory is additionally modeled, allowing us to predict, for example, production losses and material failures, minimizing unplanned downtime.
The paper is organized as follows.Section 2 describes a case study of an industrial MSRU.Section 3 presents the steady-state mathematical model of this equipment, which describes in detail each control volume investigated.The next section deals with the implementation and numerical aspects of this model.The

Industrial & Engineering Chemistry Research
experimental data necessary for its solution are explained in Section 5.In Section 6, the temperature and composition profiles predicted by the reformer model are evaluated against experimental data from industry and from the literature.The results indicate then that the developed mathematical model is able to explain the experimental behavior of the industrial MSRU.This paper ends with the final conclusions and suggestions for future works in Section 7.

PROCESS DESCRIPTION
The case study is an industrial MSRU in the Northeastern Petrochemical Company−COPENOR (Camacari, Brazil), depicted in Figure 1.The MSRU is composed of 72 CFTs (14.6 m high), which are fed by a mixture of the natural gas, steam, and CO 2 , producing the reformed gas.The combustion in the 21 burners located at the top of the furnace provides heat to highly endothermic reactions.The heat is distributed by convection and radiation mechanisms to the reforming tubes and to the refractory surface.Part of this energy is reflected, and the other part is transferred by conduction into the tubes and to the environment.

MATHEMATICAL MODEL
Figure 2 shows the four control volumes considered in the stationary MSRU model: (i) the catalyst-filled tube (CFT), (ii) tube wall, (iii) furnace, and (iv) refractory.This division of the reformer into four control volumes is considered enough to describe the stationary behavior of the equipment, thus not needing to subdivide it into more control volumes according to the spatial position, as is often seen in computational fluid dynamics approaches.In this type of approach, the furnace is divided into a finite number of isothermal zones of volume and surface, 13 which requires laborious calculations to determine the radiation heat exchange areas between these zones in the furnace.This would increase the complexity of the model, which is not desirable, since it is intended to obtain a model for practical purposes to help the operator in monitoring the temperature and composition profiles of the industrial MSRU.Despite its simplified formulation, when compared to CFD models, the model is quite comprehensive in relation to the study of different heat transfer mechanisms and its influence on the performance of the industrial reformer.
Figure 2 also illustrates the heat transfer phenomena considered as well as the main variables of the model.The heat generated by the combustion reactions inside the furnace is transferred by radiation and convection to the tube wall (Q rad f→w and Q conv f→w ) and to the refractory (Q rad f→ref and Q conv f→ref ).Part of the energy transferred to the tube is conducted across the tube wall (Q cond w ) and the other part is irradiated (Q rad w ) to the furnace (a f × Q rad w ) and to the refractory, (1 − a f ) × Q rad w , where a f is the gas absorptivity.Similarly, part of the energy transferred to the refractory is conducted across the refractory wall (Q cond ref ) and the other part is irradiated (Q rad ref ) to the furnace (a f × Q rad ref ) and to the tube wall, (1 In the interior of the tube, the heat is transferred by radiation and convection, Q rad w→g and Q conv w→g , respectively.Heat is lost to the environment from the refractory wall by convection (Q conv ref→env ).The heat transferred by convection and radiation, illustrated in Figure 2, are defined as (1)  Industrial & Engineering Chemistry Research (8) where the furnace side area (A L f ) is computed as a rectangular parallelepiped using the dimensions of the equipment (length, a = 9.36 m; width, b = 4.86 m), r in is the inner radius equal to 0.053 m, s w = 0.01 m is the tube thickness, Δz is the length of control volumes, and β 1 and β 2 are parameters equal to 0.58 and 1.0 used to correct the amount of heat transferred by radiation from the refractory and tube wall, respectively.These parameters consider, among others, the effects of irradiation and shading (shadow region caused by a tube in the adjacent ones and vice versa), caused by the geometry of the furnace and the way in which the tubes are spaced inside it.To determine the furnace gas emissivity (ε f ), several models have been used in the literature, from the simplest ones, where the furnace gas is considered as a gray gas, to the most complex ones, calculating the gas emissivity according to its composition, pressure, temperature, and furnace geometry, which influences the radiation wavelength.In the gray gas model, it is assumed that the absorption coefficient does not depend on the wavelength; however, it may be a function of the temperature and the concentration of the chemical species present in the furnace.Due to its simplicity, the gray gas hypothesis is widely used for the study of radiation in furnaces 39−42 and used in this work.To reduce the complexity of the proposed model, considering that the objective is to use it for practical purposes, the emissivity values of the gas furnace (ε f ), tubes (ε w ), and the refractory (ε ref ) are considered constant parameters and equal to 0.3758, 28 0.85, 6,25,26,28,43 and 0.60, 25,26,28 respectively.
Beek's correlation 44 is used to calculate the convective coefficient on the tube side, (h t ), 1,3,36,45,46 and Dittus−Boelter's equation 47 is used to calculate the convective coefficient on the furnace side, (h f ) 22,24,28 which assumes a turbulent flow for the gas in the furnace.The hydraulic mean diameter (D h f ) is calculated through the relationship between the cross-sectional area (free area for flow) (A) and the wet perimeter of the cross section (P). 48For the configuration of the industrial MSRU, in which the tubes are arranged in a line, this correlation is given as where p represents the tube pitch at the right angle to the direction of flow, that is, the distance between the centers of two adjacent tubes, r in is the inner radius of the tube, and s w is its thickness.The introduction of eq 11 in the proposed model allows us to study different arrangements of tubes inside the furnace, modifying only the values of A and P. Thus, the effect of this geometry on thermal exchanges and, consequently, on the temperature and composition profiles can be investigated.The convective coefficient to the external environment is a constant Ideal mixing rules are considered to compute the physical properties of the process gases in all control volumes.The specific heat (cp , i ) and specific enthalpy (h ̂i) of pure components are determined from NASA correlations. 49The dynamic viscosity (μ i ) and thermal conductivity (λ i ) are calculated according to Yaws 50 and these properties in the gas mixture (μ, λ) are determined by the Wilke method and Wassiljewa equation, respectively. 51.1.Catalyst-Filled Tube (CFT).To capture the main aspects of steam reforming, three reactions are studied: the steam reforming reaction with methane (R1), water-gas shift reaction (R2), and reverse methanation reaction (R3) Other reactions, such as carbon formation, occur less frequently and are often neglected since they would make the kinetic model too complex for practical purposes. 3,6,28,52,53everal kinetic models have been developed to represent the kinetic rates of reforming reactions, such as those reported by Akers and Camp, 54 Moe and Gerhard, 55 Bodrov et al., 56−58 Ross and Steel, 59 Al-Ubaid, 60 Allen et al., 61 Sing and Saraf, 21 De Deken et al., 62 Numaguchi and Kikuchi, 63 Xu and Froment, 52 and Hou and Hughes. 373][4][5]38,53,64 The kinetic model from Xu and Froment 52 considers the hypotheses formulated by Langmuir, 65 which describes the adsorption of the reactants on the surface of the catalyst. The kineic constants (K c, kt g ) are computed by the Arrhenius law and the adsorption constants (K a, i g ), by the Van't Hoff expression, using the parameters given by Xu and Froment, 52 according to The equilibrium constants (K e, kt g ) are calculated according to Twigg 66 g g (16)   where the parameter c in eq 14 is a conversion factor for the pressure from atom to Pa.
The gaseous mixture fed to the CFT contains different hydrocarbons, mainly methane.Higher-molecular-weight hydrocarbons, such as ethane and propane, have a relatively lower composition so that it is assumed that all the higher alkanes are irreversibly hydrocracked as soon as they enter the CFT, according to Latham 26 Assuming that the hydrocracking reaction (R4) is instantaneous and isobaric, mass and energy balances are developed to compute the composition and temperature of the equivalent stream, which feeds the CFTs, as summarized in Table 1.The new equivalent stream after the hydrocraking at the tube inlet contains six species: CH 4 , CO 2 , CO, H 2 , H 2 O, and N 2 .Figure 2 depicts the inlet variables to the hydrocracking (y i,0 g , F 0 g , P 0 g , T 0 g ) and equivalent stream (y i, in g , F in g , P 0 g , T in g ).The CFT model is formulated based on an infinitesimal control volume, assuming the following hypotheses: • The inlet stream is uniformly distributed in the CFT.
• The process gas is an ideal mixture; enthalpy variations resulting from the mixing of the components are neglected.• All the reforming tubes are identical and subject to the same temperature and pressure profiles.
• The gas inside the CFT possesses a plug flow velocity profile.
• The deposition of carbon on the catalyst surface is not considered.• There are axial temperature and concentration gradients in the CFT.• There is no distinction between the gas phase and the catalyst particles so that the temperature of the catalyst particles is uniform and equal to that of the process gas.• The catalyst particle distribution is uniform throughout the axial and radial coordinates.• The fluid flows without the existence of preferred paths.
• The physical properties of the process gas (thermal capacity, viscosity, and thermal conductivity) are functions of its temperature.• The gas flow is compressible.
• The kinetic and potential energies of all supply and output streams are negligible.The material and energy balances, which describe the CFT, are shown in Table 2.The density and flow velocity of the process gas are functions of the axial coordinate and are derived from the ideal gas equation (eq 19) and the definition of the mass flow (eq 20), respectively.The (Ni − 1) differential equations (eq 21) and the equation of the sum of the mass compositions (eq 22) define the compositions of all six components present in the gas mixture (CH 4 , CO 2 , CO, H 2 , H 2 O, and N 2 ).According to the pseudohomogeneous approach, effectiveness factors (η), indicated in Table 2, are constants multiplied to the rates of reforming reactions to account for transport limitations inside and on the surface of the catalyst particles. 10,24,25,67This simplification is acceptable in most models presented in the literature and has the reduction of the computational effort as an advantage, as stressed by Wesenberg and Svendsen. 68Ergun's equation, 69 eq 23, gives the pressure along the CFT, and the energy balance (eq 24) gives the temperature along the reactor.The density (ρ cat ), porosity (ϕ), and diameter (D cat ) of the catalyst particle are equal to 2355.
Table 2. System of Equations Describing the CFT boundary condition: (0) ∑ Industrial & Engineering Chemistry Research • m −3 , 3,23,70 0.519 e, 1,3,26 and 0.0054 m, 71 respectively.The enthalpy of the reaction kt, ΔH kt , is a function of T g and is calculated by the specific enthalpies of the products and reagents (h ̂i g ), which are obtained from their respective specific heat (cp , i g ).The effect of deactivating the catalyst is not considered in the proposed model.Its consideration would imply a deeper study of the parallel and undesirable reactions of carbon formation, among others, which is not part of the scope of this work.
3.2.Tube Wall.The energy balance at the tube wall is developed based on the following assumptions: • The properties of the material that make up the tube (thermal conductivity and emissivity) are independent of temperature.• Heat conduction through the tube wall occurs in the radial direction (r), while heat conduction in the axial direction (z) is neglected.The application of these hypotheses results in the following differential equation with boundary conditions derived from the heat transfer phenomena described by Figure where k w is the thermal conductivity of the tube, 28. 3 and a f is the absorptivity of gas in the furnace equal to 0.698.3.3.Furnace.The rigorous modeling of the energy exchange in the furnace is a complex process and requires a complete description of its geometry, tube arrangement, and burners.In addition, it is necessary to investigate the entire mechanism involved in the combustion reactions, which means evaluating the kinetic rate of numerous chemical reactions.Although the development of such a model is possible, it would require a great amount of computational effort, hindering its industrial application.Therefore, some hypotheses are made when developing the material and energy balances for an infinitesimal element at the furnace: • The heat released in the combustion reactions is transferred to the CFTs and to the refractory surface through convection and radiation.• The gas inside the furnace behaves like an ideal gas.
• The kinetic and potential energies of all supply and output streams are negligible.• The temperature gradient exists only in the axial direction of the furnace.• The pressure drop along the length of the furnace is negligible.• The physical properties of furnace gases (density, heat capacity, viscosity, and thermal conductivity) are functions of its temperature.• The combustion air and fuel gas are mixed before feeding the combustion chamber.• The furnace gas is modeled as one gray gas.
The mixture that feeds the furnace is composed by hydrocarbons, mainly methane.As carried out for the reforming tube, all the higher alkanes are assumed to be completely burned, according to the general reaction Mass and energy balances are then developed to calculate the composition and temperature of the equivalent stream fed to the furnace burners.The gaseous feed then contains seven species: CH 4 , CO 2 , CO, H 2 O, H 2 , O 2 , and N 2 .Table 3 summarizes the equation system to compute the new equivalent stream after the precombustion at the furnace inlet.Figure 2 depicts the inlet variables to the precombustion (y i,0 f , F 0 f , P 0 f , T 0 f ) and the equivalent stream (y i, in f , F in f , P 0 f , T in f ).This equivalent stream is used as input data to the furnace model, which is summarized in Table 4.The density and flow velocity of the gases in the furnace change along the axial coordinate and are derived from the ideal gas equation (eq 30) and the definition of the mass flow (eq 31), respectively.Similar to the CFT, the component mass balance is described by Ni − 1 differential equations (eq 32) and one algebraic equation (eq 33), which sums up the mass fractions.The energy balance expressed in eq 34 provides the furnace temperature profile, which is significantly influenced by the kinetic model used to describe the combustion reactions.The enthalpy of the combustion reaction kf, ΔH kf ,is the function of T f and is calculated by the specific enthalpies of the products and reagents (h ̂i f ), which are obtained from their respective specific heat (cp , i f ), determined from the NASA correlations. 49he literature usually employs empirical expressions for the distribution of heat along the furnace and assumes a homogeneous composition in the furnace.Unlike this usual approach, the proposed model incorporates phenomenological combustion kinetics and heat transfer phenomena to predict composition and temperature profiles along the furnace.The combustion reactions of CH 4 , CO, and H 2 are described according to Tran et al. 19 + Tran et al., 19 when developing a model of a MSRU based on CFD, suggest the following kinetic model for the combustion reactions Table 3. Model of Precombustion to Compute the Equivalent Stream for the Furnace Industrial & Engineering Chemistry Research where K c, kf f are the kinetic constants of reactions kf = 6, 7, and 8, calculated by the Arrhenius equation; C i f are the concentrations of the species i involved in the reactions, and α j , where j = 1, ...,7 are the reaction orders of each of these species.The authors consider in the mass and energy balances the mixing effect between fuel, oxygen, and products by relating reaction rates with the rate of turbulent energy dissipation since they are concerned with a CFD model.Our approach, however, does not consider the mixing effect since this would increase the complexity of the proposed model and the required computational effort.
The kinetic parameters (activation energy, pre-exponential factor, and order of reaction) are qualitatively adjusted to validate the temperature profile observed in industrial MSRUs, which presents a peak temperature around 1/3 from the top of the furnace.The rigorous estimation of such parameters is difficult due to the scarcity of experimental and literature data of the temperature profile along the length of the furnace.In the MSRU investigated by Latham, 26 for example, only two external wall temperature measurements of the tube are available, one located at 3.57 m and the other at 8.56 m from the top of the reformer, which correspond to averages of different tubes.Furthermore, such measurements have a high degree of uncertainty associated with the infrared devices used to measure the temperature at the radiant furnaces.The availability of temperature measurements at different points along the length of the equipment would allow a more rigorous and quantitative estimation of the kinetic parameters.
3.4.Refractory Wall.The modeling of the refractory wall may contribute to an increase in the on-service life of the refractory material.Therefore, unlike most studies in the literature, the proposed model computes the temperature gradient in the refractory wall.The following hypotheses are considered: • The heat flow through the refractory wall occurs by conduction and from there to the environment by convection and radiation.• The properties of the refractory (thermal conductivity and emissivity) are independent of temperature.• Thermal energy is exchange by radiation and convection with CFT and furnace gases.• Heat conduction through the refractory wall occurs in the radial direction (r), while heat conduction in the axial direction (z) is neglected.The following energy balance for the refractory wall results in with boundary conditions derived from the heat transfer phenomena described by Figure 1 where k ref is the thermal conductivity of the tube, 2.6 W•m −1 • K −1 . 19It is convenient to remember the influence of this parameter on the insulation of the furnace.The higher the thermal conductivity of the refractory material, the greater the heat losses to the environment.These losses demand a higher consumption of fuel gas to meet the thermal load necessary for the reforming reactions.Therefore, it is a parameter that can be used for purposes of analyzing the efficiency of the MSRU.In a simplistic way, this efficiency can be estimated based on the knowledge of heat losses through the furnace and the energy released by burning fuel gas.This heat released in combustion Table 4. System of Equations Describing the Furnace boundary condition: (0) Industrial & Engineering Chemistry Research can be determined, in turn, by the product between the amount of fuel used and its calorific value.Another possibility of estimating efficiency would be by determining the fraction of this released energy that is transferred to the process gas inside the reform tubes.

IMPLEMENTATION AND NUMERICAL ASPECTS
First, the equivalent streams that feed the CFTs and the furnace (y i, in g , y i, in f , T in g , T in f ) are computed, as indicated in Tables 1 and 3.The nonlinear systems are solved by the BzzNonLinearSystem numerical solver, available from the BzzMath library developed by Buzzi-Ferraris. 72he mathematical model developed for the MSRU results in a differential algebraic equation (DAE) system, which involves control volumes that show different dynamics.On the one hand, the reforming reactions proceed slowly.On the other hand, the combustion reaction in the furnace occurs very quickly, resulting in highly steep temperature and concentration profiles.Thus, the model is represented by a system of stiff ordinary differential equations coupled with algebraic equations, which requires accurate, robust, and efficient numerical methods.Different numerical methods have been proposed in the literature to solve these stiff DAE systems: LSODI, 73 DASSL, 74 DASPK, 75 and BzzDae. 76−80 The greater robustness of the BzzDae method allows the user to overcome the faults from other codes and to reduce the critical selection of relative and absolute tolerances.In addition, the efficiency of the BzzDae method in terms of computing time is significantly enhanced mainly due to better Jacobian evaluation criteria.

MODEL VALIDATION
The model is validated against experimental data from the literature 26 and from our industrial partner, which are summarized in Table 5.In the following tables and figures, cases 1 to 3 refer to the experimental data provided by Latham 26 and Cases 4 and 5 to the data obtained by our industrial partner.

RESULTS AND DISCUSSION
Despite the high operating temperature of the process gases, the presence of hydrogen in the gas mixture and the high operating pressure of the CFT raise the question of whether the hypothesis of ideal gas is valid.Attempting to investigate this ideal behavior, as considered by most of the works reported in the literature, some tests are performed in the separation and phase equilibrium calculations (SPECS) simulator, from Technical University of Denmark, using the Soave−Redlich−Kwong equation, with classical mixing rule and binary interaction parameter k ij = 0.0. 81The compressibility factors are computed at the inlet and outlet of the tube and furnace as well as at the point with the highest temperature along the furnace, referred to as a peak position.Two data sets are selected: case 2 (Latham 26 ) and case 4 from the industrial partner.Table 6 shows that the compressibility factors (Z) for each operating condition are very close to unity, confirming that the ideal gas hypothesis holds true.
Table 7 shows the relative deviations (RD) obtained for the flue gas temperature (T out f ) and for the temperature (T out g ), pressure (P out g ), and composition (y i, out g ) of the reformed gas at the outlet of the MSRU.The experimental data were taken from Latham. 26For each case, the RDs from the model by Latham 26 are presented together with the RD calculated in this work.The calculated pressure at the outlet of the CFT is not shown by Latham, 26 and then Table 7 does not present the corresponding RD.The proposed model can well reproduce the experimental data.When comparing our model predictions with the simulated data from Latham, 26 similar RDs are observed for all variables, except methane composition, whose RD was significantly different from that obtained by the proposed model.
Figure 3 compares the temperature and composition profiles of both models.The results for cases 1 and 2 are similar to case 3; therefore, only case 3, which presented smallest deviations (Table 7), is illustrated in Figure 3 for the sake of simplicity.The  Industrial & Engineering Chemistry Research industrial data reported by Latham 26 are also illustrated in Figure 3: the temperature of the flue gas at the outlet (T out f ), temperature (T out g ), and composition (y i, out g ) of the reformed gas at the outlet as well as the outer wall temperature (T ext w ) at two intermediate points along the axial coordinate.At the first section of the reactor, the profiles differ significantly due to the different modeling approach of the furnace.Unlike Latham, 26 who provided the heat distribution profile of the furnace, our model employs the kinetic rates to describe the combustion reactions, which proceed more slowly.The proposed model presents an abrupt decrease of nearly 100 K in T g just after the reactor inlet because the reforming reactions occur very quickly due to the high methane concentration at the inlet.This indicates the endothermic net effect of the parallel reforming reactions.Then, the reforming reactions occur slowly and a decrease of 20 K is observed for T g .Around 1.6 m, the furnace temperature abruptly increases, providing the necessary heat for increasing the reaction rates and the temperature inside the tube.Due to the scarcity of experimental data along the axial coordinate, it is not possible to affirm which of the simulated profiles for the composition and temperature are closer to the industrial unit.From 3 m onward, both models behave similarly and can predict the experimental data with accuracy since small relative deviations are obtained, as Figure 3 illustrates.As expected, H 2 is formed to the extent that steam and methane are consumed.Moreover, there is an amount of CO higher than that of CO 2 in the reactor outlet, resulting from the reversal of watergas shift reaction (R2) at high temperatures due to its exothermic nature.
It was not possible to directly validate the temperature profiles on the inner and outer surfaces of the refractory since the experimental data of these variables are not available, neither in the literature nor in industry.Such temperature profiles are more closely correlated with the compositions and temperature profiles of the tube and furnace (Figure 3).As these profiles are consistent, as discussed previously, we assume that the equations describing the temperature gradients at the inner and outer surfaces of the refractory wall are valid.
Table 8 presents the deviations of the model predictions against the experimental data from our industrial partner, cases 4 and 5.The model can predict the data with accuracy since deviations are below 6%, except for the furnace temperature, which is around 10%.Even this higher deviation is acceptable because the reliability of the experimental data cannot be guaranteed.In the industrial MSRU, flue gas temperaturemeasuring instruments are considered critical since they need periodic calibrations.
6.1.Sensitivity Analysis on Performance.The proposed model considers some phenomena often neglected in the literature: (i) the energy irradiated by the refractory (Q rad ref , eq 8) and by the tube wall (Q rad w , eq 9), (ii) the heat transfer by convection from the furnace to the tube wall (Q conv f→w , eq 2) and to the refractory (Q conv f→ref , eq 3), and (iii) the heat transfer by radiation from the tube wall to the process gas (Q rad w→g , eq 7).Attempting to understand how these phenomena influence the model prediction, sensitivity analysis is carried out, taking the case 3 as the base case as it presented smallest deviations.

Industrial & Engineering Chemistry Research
To evaluate the energy irradiated by the refractory and by the tube wall (i), the corrective factors β 1 and β 2 , respectively, are investigated.Figure 4 shows the deviations against the experimental data from Latham 26 for the output variables when these parameters are varied from 0 to 1.When the energy irradiated by the refractory and tube wall is neglected (β 1 = 0 or β 2 = 0), the model loses the prediction capacity as high relative deviations are observed, mainly for methane and carbon monoxide compositions.Furthermore, in such cases, the refractory wall temperature profile is inconsistent, with an excessively high temperature peak (above 10,000 K).When the energy irradiated by the refractory is considered, β 1 > 0, the deviations of tube and furnace outlet variables do not change significantly and the temperature profiles are similar to those obtained in the validation procedure (Figure 3, case 3); therefore, the model adopts β 1 = 1.Regarding the energy irradiated by the tube wall, smaller deviations are observed when 0.5 ≤ β 2 ≤ 0.75; then, after fine-tuning inside this range, β 2 = 0.58 was chosen to validate the proposed model.The temperature profiles of the furnace, T f , the inner refractory wall, T in ref , the process gas inside the CFT, T g , and the outer wall of the tube, T ext w , for different values of β 2 are illustrated in Figure 5a−d, respectively.As tube irradiation increases (higher β 2 ), higher furnace and refractory wall temperatures are reached, while process gas and tube wall temperatures decrease.This is due to the lower energy availability for processing of the reforming reactions, which leads to a larger amount of methane  Industrial & Engineering Chemistry Research in the reformed gas.Therefore, the energy irradiated by the tube wall plays a key role in describing the process with accuracy.
If the heat transfer by convection from the furnace to the tube wall (Q conv f→w , eq 2) and to the refractory (Q conv f→ref , eq 3) is neglected (ii), no significant variation in the responses predicted by the model are observed, as Figure 6a illustrates.This might be associated with the dominant effect of radiation when compared with the convection heat transport mechanism, as observed in Hottel and Sarofim, 13 Roesler, 14 Selcuk et al., 15,16 Keramida et al., 17 Ebrahimi et al., 18 and Tran et al. 19 Therefore, the convective term can be neglected without damage to the accuracy of the model.The heat transfer by radiation from the tube wall to the process gas, Q rad w→g (iii), on the other hand, should be considered to improve the model fit.As Figure 6b illustrates, the methane composition reaches relative deviations of around 20% if this term is neglected, although the other variables show deviations up to 10%.
Due to the dominant effect of the radiation mechanism, sensibility analysis was carried out to investigate the influence of the emissivity of the tube (ε w ) and the refractory (ε ref ) on the performance of the proposed model.For a tube emissivity range from 0.75 to 0.95, the sensitivity analysis indicates that the model is not sensible to this parameter.Regarding the emissivity of the refractory, some values are reported in the literature, among which can be cited: 0.6, 3,25,26,28 0.7, 23,43 0.75, 26 0.83, 82 and 0.9. 83Figure 7 illustrates that these emissivity values have a considerable impact on the temperatures of the refractory inner and outer walls.When ε ref = 0.9, the maximum temperatures of the refractory inner and outer surfaces decrease around 9% with respect to the values found for the base case (ε ref = 0.6).The other output variables were insensitive to these variations.Such results can be attributed to the gray gas model used to investigate radiation heat transfer.This approach proved to be unable to capture the effect of changing refractory emissivity on the temperature profiles and composition of the tubular reactor.These results agree with the study developed by Adams and Olver, 83 who investigated the effects of different coatings in a process furnace using the computational fluid dynamics approach.These authors changed the refractory emissivity from 0.4 to 0.9 and negligible differences were observed in the process variables when using only one gray gas model.On the other hand, when they performed the same change using a spectral gas line model (spectral-line-weighted, SLW), the efficiency of the radiant oven increased by 2% and the arc temperature decreased by 30 °C.Therefore, despite the observed results in Figure 7, it is not possible to conclude that the quality of the reformed gas is not influenced by this parameter.Due to its simplicity, the gray gas model causes greater deviations from the real behavior of gases in relation to the radiation absorption spectrum.Further studies can be conducted using a more rigorous model that describes, in more detail, the transfer of radiant heat in the furnace gases.

CONCLUSIONS
This paper presented a rigorous simulation of an industrial MSRU, deeply describing heat transfer phenomena as well as the combustion kinetics through a phenomenological model.The model considers some commonly neglected heat transfer mechanisms, such as the intraphase exchange of radiant energy, the absorption of part of this energy by the furnace gases, and the  convection mechanisms in the furnace and radiation inside the tube.The results show the importance of some heat transfer mechanisms that are usually neglected (e.g., radiation inside the tubes) while revealing that other typically considered mechanisms may be neglected (heat transfer by convection inside the furnace).The detailed modeling of the combustion kinetics in the furnace allows the prediction of the concentration and temperature profiles without the need to previously fix such profiles through empirical equations.To validate the model, simulations were performed using data from the literature and an industrial partner.The proposed model compares favorably to the results reported by Latham, who used a more detailed approach to describe the radiation in the furnace.It was observed that the temperature profiles predicted by both models are similar.Additionally, the proposed model also reproduced the outlet conditions of the CFT and furnace with reasonable accuracy.The formulation of the proposed model through four control volumes proved to be adequate for the purpose of monitoring the stationary behavior of the industrial MSRU, so it is not necessary to subdivide the equipment into more control volumes according to the spatial position.In this context, it is possible to conclude that the simplifications adopted in the MSRU model are valid and do not produce errors in the simulation.The applicable range of the model proposed is not limited as a result of adopting these simplifications, since its validation was conducted from the adjustment of parameters related to radiation heat transfer, such as β 1 and β 2 .Therefore, for practical purposes, these parameters can be re-estimated when necessary.In fact, almost all the models presented in the literature that describe radiation in the furnace using more detailed approaches require some sort of data fitting.Such models often consider strong simplifying hypothesis, with heuristic profiles, which are functions of parameters determined by data fitting.This shows that parameter re-estimation is a common, valid strategy, which supports the method presented in this work.Due to the use of only one gray gas model, the effect of the variation of the refractory emissivity on the output variables could not be visualized in the proposed model, recommending the use of a more detailed model for the study of the radiation heat distribution in the industrial MSRU.Therefore, the proposed model offers the advantage of better understanding how the heat transfer mechanisms influence the concentration and temperature profiles throughout the equipment.We suggest, as a continuation of this work, to verify the influence of different kinetic models of the reforming reactions in the prediction of the model proposed for the MSRU.Additionally, this industrial unit must be investigated using a time-dependent model to be able to predict the composition and temperature profiles during specific equipment operation (start-up, shutdown, perturbations, etc.).This dynamic model can serve as a useful tool for investigating control systems, for example.

Figure 2 .
Figure 2. Representation of the control volumes, the heat transfer mechanisms, and main variables of the industrial reformer model.

Figure 3 .
Figure 3. Temperature and composition profiles predicted by the model proposed in this work and the one by Latham: 26 case 3.

Figure 5 .
Figure 5. Simulated temperature profiles for changing values of β 2 : (a) temperature of the flue gas (T f ), (b) refractory inner surface (T in ref ), (c) process gas (T g ), and (d) tube outer surface (T ext w ).Orange dots represent experimental measurements from Latham. 26

Figure 6 .
Figure 6.Effect of the convective heat transfer in the furnace side (Q conv f→w and Q conv f→ref) and the radiation in the tube side (Q rad w→g ) on the model accuracy.

Figure
Figure Predicted temperature profiles for refractory inner (T in ref ) and outer surface (T ext ref ) when ε ref is varied.

1 )
factor of the radiation emitted by the refractory β 2 corrective factor of the radiation emitted by the tube η kt effectiveness factor of the reforming reactions kt = 1, ...,3 K e, kt g equilibrium constant of the reformation reactions kt = 1, ...,3 ξ kr extent of hydrocracking reaction kr ξ kc extent of precombustion reaction kc T env external environment temperature (K) v f flow velocity of the furnace gas (m • s −1 ) v g flow velocity of the process gas (m • s −1 ) k f flue gas thermal conductivity (W • m −2 • K −1 ) μ f flue gas viscosity (Pa • s −to the furnace before precombustion (K) Q rad ref heat reflected by the refractory (J • s −1 ) Q rad w heat reflected by the tube wall (J • s −1 ) Q cond ref heat transferred by conduction through the refractory wall (J • s −1 ) Q cond w heat transferred by conduction through the tube wall (J • s −1 ) Q conv f→ref heat transferred by convection from the furnace to the refractory (J • s −1 ) Q conv f→w heat transferred by convection from the furnace to the tube wall (J • s −1 ) Q conv ref→env heat transferred by convection from the refractory to the external environment (J • s −1 ) Q conv w→g heat transferred by convection from the tube wall to the process gas (J • s −1 ) Q rad f→ref heat transferred by radiation from the furnace to the refractory (J • s −1 ) Q rad f→w heat transferred by radiation from the furnace to the tube wall (J • s −1 ) Q rad w→g heat transferred by radiation from the tube wall to the process gas (J • s −1 ) D h f hydraulic mean diameter (m) R ideal gas constant (m 3 • Pa • kmol −1 • K −1 ) r in inner radius of CFT (m) K c, kf f kinetic constant of the combustion reactions kf = 6, ...,8 Δz length of each control volume element (m) w i f mass fraction of component i in the furnace gas w i, in f mass fraction of component i in the gas fed to the furnace after precombustion w i g mass fraction of component i in the process gas w i, in g mass fraction of component i in the process gas fed to tube after hydrocracking h ̃i, in f molar enthalpy of component i in the furnace gas after precombustion (J • kmol −1 ) h ̃i,0 f molar enthalpy of component i in the furnace gas before precombustion (J • kmol −1 ) h ̃i,0 g molar enthalpy of component i in the process before hydrocracking (J • kmol −1 ) h ̃i, in g molar enthalpy of component i in the process gas after hydrocracking (J • kmol −1 ) F in f molar flow of furnace gas after precombustion (kmol • s −1 ) F 0 f molar flow of furnace gas before precombustion (kmol • s −1 ) F in g molar flow of process gas fed to CFT after hydrocracking (kmol • s −1 ) F 0 g molar flow of process gas fed to CFT before hydrocracking (kmol • s −1 ) y i, in f molar fraction of component i in the furnace gas after precombustion y i,0 f molar fraction of component i in the furnace gas before precombustion y i, in g molar fraction of component i in the process gas fed to the tube after hydrocracking y i,0 g molar fraction of component i in the process gas fed to the tube before hydrocracking MW f molecular weight of furnace gas (kg • kmol −1 ) MW in f molecular weight of furnace gas after precombustion (kg • kmol −1 ) MW g molecular weight of process gas (kg • kmol −1 ) MW in g molecular weight of process gas after hydrocracking (kg • kmol −1 ) MW i f molecular weight of species i in the gas furnace (kg • kmol −1 ) MW i g molecular weight of species i in the process gas (kg • kmol −1) Nkr number of hydrocracking reactions Nkc number of precombustion reactions Nkf number of reactions in the furnace Nkt number of reactions inside the CFT α j order of reaction with respect to component i, where i ∈ {CH 4 , CO, CO 2 , O 2 , H 2 O} and j = 1, ..., 7 ϕ porosity of the catalytic bed A a, i g pre-exponential adsorption factor of component i ∈ {CH 4 , CO, H 2 O, H 2 } A c, kt g pre-exponential factor of reforming reactions kt = 1, ...,3 ρ g process gas density (kg • m −3 ) P gprocess gas pressure (Pa) P 0 g process gas pressure fed to the CFT (Pa) T g process gas temperature (K) T in g process gas temperature after hydrocracking (K) μ g process gas viscosity (Pa • s −1 ) T 0 g process gas temperature before hydrocracking (K) r radial coordinate (m) E a, kt g reaction activation energy of reforming reactions kt = 1, ...,3 (J • kmol −1 ) ΔH kf reaction enthalpy of combustion reactions kf = 6, ...,8 (J • kmol −1 ) ΔH kt reaction enthalpy of reforming reactions kt = 1, ...,3 (J • kmol −1 ) ε ref refractory emissivity s ref refractory thickness (m) h ̂i f specific enthalpy of component i in the furnace gas (J • kg −1 ) h ̂i g specific enthalpy of component i in the process gas (J• kg −1 ) ĉp , i f specific heat of component i in the furnace gas (J • kg −1 • K −1 ) ĉp , i g specific heat of component i in the process gas (J • kg −1 • K −1 ) cp f specific heat of the furnace gas (J • kg −1 • K −1 ) cp gspecific heat of the process gas (J• kg −1 • K −1 ) ΔH i 0 standard chemical adsorption energy of component i ∈ {CH 4 , CO, H 2 O, H 2 } (J • kmol −1 ) σ Stefan−Boltzmann constant (s −1 • m −2 • K −4 ) υ i, kff stoichiometric coefficient of combustion reactions kf = 6, ...,8 υ i, kr g stoichiometric coefficient of hydrocracking reaction kr υ i, kc f stoichiometric coefficient of precombustion reaction kc υ i, kt g stoichiometric coefficient of reforming reactions kt = 1, ...,3 T in ref temperature in the inner refractory wall (K) Industrial & Engineering Chemistry Research T in w temperature in the inner wall of the tube (K) T out ref temperature in the outer refractory wall (K) T out w temperature in the outer wall of the tube (K) k ref thermal conductivity of the refractory (W • m −1 • K −1 ) ṁf total furnace gas mass flow (kg • s −1 ) ṁg total process gas mass flow (kg • s −1 ) A t tube cross-sectional area (m 2 ) ε w tube emissivity k w tube thermal conductivity (W • m −1 • K −1 ) fraction in the process gas before precombustion P wet perimeter of the cross section (m)

Table 1 .
2 kg Model of Hydrocracking to Compute the Equivalent Stream for the CFT

Table 5 .
Tube-Side and Furnace-Side Inlet Conditions 71a aValues associated to chemical species are molar fractions (in %).

Table 6 .
Compressibility Factor (Z) Computed at Different Operating Conditions

Table 7 .
Relative Deviations (in %) Obtained from the Model by Latham 26 and from the Model Proposed in This Work (Experimental Data from Latham 26 )

Table 8 .
Relative Deviations (in %) between the Results Predicted by the Proposed Model and the Industrial MSRU under Study The authors would like to thank COPENOR, especially John Kennedy Fernandes and Alan Rocha dos Santos Pinho Costa for providing data and knowledge about the industrial plant under study and for authorizing the publication of this information in this article; the Coordination for the Improvement of Higher Education Personnel (CAPES), the National Council for Scientific and Technological Development (CNPQ), and the CARIPLO foundation for the financial support; and Professor Eliseo Ranzi and Dr. Ing.Andrea Bassani for their suggestions and support regarding the modeling approach.