In order to predict the performance of the SMR process, it is necessary to simulate the tube side process and provide a detailed profile of the heat flux, gas composition, carbon forming potential and the pressure inside of the reformer tubes incrementally. The calculations involve solving material and energy balance equations along with reaction kinetic expressions for the nickel catalyst.
The general overall reaction for the steam reforming of any hydrocarbons can be defined as Equation (i) (1, 2):
In this work, steam reforming of the natural gas is described with the following equations, as the methane is the major constituent of the natural gas. Equation (ii) (1, 2):
In parallel with this SMR equilibrium, the water gas shift (WGS) reaction proceeds according to Equation (iii) (1, 2):
Minette et al. (17) in their work stated that the second SMR reaction is often not accounted for assuming it follows directly from combining Equations (i) and (ii). However, the work of Xu and Froment (7–8) showed that the second SMR reaction expressed by Equation (iv) follows an independent reaction path and must be accounted for in combination with Equations (i) and (ii), as confirmed by the measurements of Minette et al. (17):
As mentioned, the described reactions proceed in indirectly heated reformer tubes filled with nickel-containing reforming catalyst and are controlled to achieve only a partial methane conversion. In a top fired reformer usually up to 65% to 68% conversion based on methane feed can be accomplished, leaving around 10 mol% to 14 mol% methane per dry basis (1, 2).
The overall SMR reaction of methane is endothermic and proceeds with an increase of volume at the elevated pressure of 20 bar to 40 bar and temperatures from 800°C to 1200°C at the exit of the reformer tubes in the presence of metallic nickel catalyst as an active component. Besides pressure and temperature, the S/NG molar ratio has a beneficial effect on the equilibrium methane concentration (18).
Another reason for applying the appropriate (higher) S/NG molar ratio is to prevent carbon deposition on the reforming catalyst. The side effect of carbon deposition is a higher pressure drop and the reduction of catalyst activity. As the rate of endothermic reaction is lowered, this can cause local overheating of the reformer tubes (hot spots and bands) and the premature failure of the tube walls. The carbon formation may occur via Boudouard reaction, methane cracking and carbon monoxide and carbon dioxide reduction. These reactions are reversible with dynamic equilibrium between carbon formation and removal. Under typical steam reforming conditions, Boudouard reaction and carbon monoxide and carbon dioxide reduction cause carbon removal, whilst methane cracking leads to carbon formation in the upper part of the reformer tube (19). Greenfield SMR units based on natural gas regularly use a S/NG molar ratio of around 3.0, while older installations are in the range from 3.5 to 4.0 (1). From the theoretical point of view any S/NG molar ratio which is slightly over 1.0 will prevent cracking, because the rate of carbon removing reactions is faster than the rate of carbon deposition reactions. However, from the practical point of view (catalyst limitations and sufficient quantity of steam for the downstream process step of WGS conversion), the minimum molar ratio which applies at the industrial level is 2.5. To account for all these facts, the model was validated for S/NG molar ratios in the range from 2.0 to 6.0.
The nickel content in relation to the composition and structure of the support differs considerably from one catalyst supplier to another. This is the reason why it is difficult to relate data from industrial plants to laboratory experiments. Reformer simulations frequently use a numerical approach in which the experimental data serves for reaction rate calculations which are described by closed analytical expressions. From the reaction rates perspective, it is possible to calculate the equilibrium gas composition for a given pressure and S/NG molar ratio at different temperatures. On top of this, the equilibrium curve which is defined by the corresponding enthalpy changes versus temperature also presents a useful parameter in the estimation of the catalyst performance. The comparison of the mentioned equilibrium curves with the working curves (working point) and the subsequent operator’s adjustments of the influencing process parameters according to the evaluated recommendations seem a useful tool to improve the catalyst performance.
In order to describe the kinetic conditions which are necessary for the determination of equilibrium methane molar concentration (a measure for the theoretical conversion) and enthalpy change over different nickel catalysts in relation with temperature at different S/NG molar ratios and reforming pressures, the model uses the following reaction rates for the equilibrium Equations (ii) to (iv) (7–8, 20), Equations (v)–(viii):
where r presents the reaction rates for methane, carbon monoxide and carbon dioxide in kmol m–3 s–1; p stands for the species partial pressures (in atm); T is the temperature (in K); while R is the gas constant (in kJ kmol–1 K–1).
Kinetic rate constants ki are given by the general Arrhenius relationship, Equation (ix) (7–8, 20), where i denotes the number of reactions from Equation (i) to (iii):
The units of k2 and k4 (Equation (ii) and (iv)) are kmol bar0.5 kg–1cat h–1), while the unit of k3 (Equation (iii) is kmol bar–1 kg–1cat h–1).
Table I (20) gives the parameters for the activation energies, Ei, and for the pre-exponential factors, Ai, used in the model, valid for most of the commercial nickel catalysts with either MgAl2O4 or CaAl12O19 support.
Table I
Parameters for the Activation Energies, E i, and for the Pre-Exponential Factors, Ai
| Equilibrium reaction |
Activation energy, Ei
|
Pre-exponential factor, Ai
|
| Unit |
Value |
Unit |
Value |
| Reaction no. 2 |
kJ mol–1 |
–240.100 |
kmol bar0.5 kg–1cat h–1 |
4.22 × 1015 |
| Reaction no. 3 |
kJ mol–1 |
–67.130 |
kmol bar–1 kg–1cat h–1 |
1.96 × 106 |
| Reaction no. 4 |
kJ mol–1 |
–243.900 |
kmol bar0.5 kg–1cat h–1 |
1.02 × 1015 |
Apparent adsorption equilibrium constants Ki in Equation (x) are defined by the general expression given in (7–8, 20), where i denotes the species in Equations (i), (ii) and (iii) or methane, water, hydrogen and carbon monoxide:
Bi is the pre-exponential factor expressed in bar-1 or unitless, while ΔHi is the absorption enthalpy change expressed in kJ mol–1.
Table II presents the pre-exponential factors and the absorption enthalpy changes for species given in Equation (x), and the same is also valid for most of the commercial nickel catalysts with either MgAl2O4 or CaAl12O19 support.
From Equations (v) to (vii) it can be concluded that the concentration of hydrogen cannot be zero, because dividing with zero would make calculated reaction rates infinite. So, according to this, it is necessary to ensure the minimum content of hydrogen in the natural gas stream to ensure applicability of these equations in the model. From the process side, hydrogen is necessary for two reasons. Firstly, it is important for the removal of organic sulfur compounds present in the natural gas by the cobalt-molybdenum catalyst, as sulfur is a poison for the nickel catalyst (reaction between organic sulfur compounds and hydrogen to give hydrogen sulfide which is subsequently absorbed by zinc oxide bed). Secondly, hydrogen will always keep the nickel catalyst in the reduced state of metallic nickel and hence maintain adequate catalyst activity in the reformer tubes.
Table II
Parameters for the Pre-Exponential Factor, Bi, and for the Absorption Enthalpy Changes ΔHi
| Species |
Pre-exponential factor, Bi
|
Absorption enthalpy change, ΔHi
|
| Unit |
Value |
Unit |
Value |
| Methane |
bar–1 |
6.65 × 10–4 |
kJ mol–1 |
38.280 |
| Water |
– |
1.77 × 105 |
kJ mol–1 |
–88.680 |
| Hydrogen |
bar–1 |
6.12 × 10–9 |
kJ mol–1 |
82.900 |
| Carbon monoxide |
bar–1 |
8.23 × 10–5 |
kJ mol–1 |
70.650 |
From the general stoichiometry and according to defined reaction rates, the model can calculate the molar flow rates of species i in kmol h–1 in the presence of an adequate quantity of nickel catalyst with the ultimate result of methane and water conversions. The relations used to determine the methane and water conversions are as follows (21, 22), Equations (xi)–(xii):
A denotes the catalyst tube cross-sectional area in m2; ρB represents the catalyst bed density in kg m–3; Fi is the molar flow rate of the species methane and water in kmol h–1; while ηi is the effectiveness factor for methane and water.
To account for the variations in reaction rate throughout the catalyst pellet, a parameter called effectiveness factor, η, is defined. This is the ratio of the overall reaction rate in the catalyst pellet and the reaction rate at the external surface of the catalyst pellet. Effectiveness factor is a function of Thiele modulus, Φ, which is related to the catalyst volume and the external surface area of the catalyst pellets. Taking into account reaction rates given by Equations (v)–(vii) and following the mechanism given by Xu and Froment (7, 8), the effectiveness factor can be calculated from Equation (xiii):
where p is the partial pressure of the species in bar; r presents the reaction rates for methane, carbon monoxide and carbon dioxide in kmol m–3 s–1; while ξ is the dimensionless intracatalyst coordinate.
Effectiveness factor profiles along the length of the reformer tube are calculated for all key species given in Equations (ii) to (iv) by solving two-point boundary differential equations for the catalyst pellets with the help of scripts and functions in the form of m-files, which was reconciled with the data from the simulator flowsheet.
The algorithm uses the following relationship for calculation of species concentration profiles inside the catalyst layer under reconciled conditions (17), Equations (xiv)–(xv):
with the corresponding boundary conditions, Equations (xvi)–(xvii):
where ξ is the dimensionless intracatalyst coordinate; De,A is the species effective diffusivity in m3fluid m–1catalyst s–1; p denotes the partial pressure of species in bar; R is the universal gas constant in kJ kmol–1 K–1; T is the bulk fluid temperature in K; h is catalytic layer thickness in m and ρs is the active solid density in kgcatalyst m–3catalyst.
The interfacial (gas-solid) mass and heat transfer limitations are negligible and were not accounted for, because the high volume flow velocity and sufficient turbulence have been assumed which reflects the operation conditions inside of the reformer tubes.
Due to model simplification and minimisation of the computational time the simplest geometry of a slab of catalyst has been assumed, which is a satisfactory assumption for the computational routine required for industrial application. The model has been tested with coating thickness in the range from 10 μm to 50 μm and the best fit with the actual process data was achieved with the catalyst coating of 10 μm.
The species effective diffusivity is determined by Equation (xviii):
where ɛs is the internal void fraction or porosity of the catalyst in m3fluid m–3catalyst; τ denotes the catalyst tortuosity and 
is the average diffusivity of species A.
The average diffusivity of species is determined by Equation (xix):
where DA is the diffusivity of the reacting species A given by Equation (xx) and S(rp,i) is the void fraction taken by the pores with radii ranging from rp,i to rp,i +1:
where DkA is the Knudsen diffusivity in m3fluid m–1catalyst s–1.
In order to have an appropriate computational speed of effectiveness factor (which is performed by m-file), the actxserver command is used for the interconnection through the COM automation server that controls the simulator. The COM interface establishes a two-way communication between the simulator and MATLAB® through shared memory block, which is built as level-2S-function. The approximation of the catalyst effectiveness factor is determined by correlating the kinetic model results with the plant process data, and the model is validated to get maximum alignment with the actual process data.
Conversions of methane and water are calculated by Equations (xxi)–(xxii) (22):
The Ergun equation for the determination of the pressure drop across the plug flow reactor (PFR) is used and solved as an ordinary differential equation (23–31), Equation (xxiii):
where ρ denotes the pressure in bar; ρ is the fluid density in kg m–3; v is the fluid velocity in m s–1; dp is the catalyst particle diameter in m; ∈ is the catalyst void fraction and Re is the particle Reynolds number.
The temperature variation of the reacting mixture (natural gas and steam) along the reformer tube is calculated according to the following relationship, Equation (xxiv):
where G is the reacting mixture flow rate in kg h–1; 
denotes average specific heat of the gas mixture in kJ kg–1 K–1; U is the overall heat transfer coefficient between the reformer tubes and their surrounding in m2 h K kJ–1; Tt,0 is the temperature of the furnace that surrounds the reformer tubes; ΔHi is the enthalpy change in kJ kmol–1; ρB represents the catalyst bed density in kg m–3; ηi is the effectiveness factor for each of the species in reacting mixture and ri is the reaction rates in kmol m–3 s–1.
The reformer catalyst tubes are simulated as PFR in which the flow field is modelled as plug flow, implying that the stream is radially isotropic (without mass or energy gradients). According to this, axial mixing is negligible. As the reactants flow the length of the reformer tube, they are continually consumed, hence, there is an axial variation in the concentration. Since reaction rate is a function of concentration, the reaction rate varies axially. To get the solution for the PFR (axial profiles of compositions, temperature and so forth) the reformer tubes are divided into several sub-volumes. Within each sub-volume, the reaction rate is spatially uniform. A mole balance executes routine calculation procedure in each sub-volume j according to Equation (xxv) (28, 29):
Because the reaction rate is spatially uniform in each sub-volume, the third term reduces to rjdV and at steady state, the above expression reduces to Equation (xxvi):
The firing side (furnace combustion model) was simulated according to the previous work of Zečević and Bolf (32) which is able to calculate adiabatic and real flame temperatures, quality and quantity composition of the waste gases, according to the known composition of the fuel gas and inlet temperatures of fuel and combustion air, with possibility to control all critical process parameters by implementation of proposed gain-scheduled model predictive control.
The basic input requirements for the model are:
-
Integration information: number of reformer tube segments, minimum step fraction, minimum step length
-
Tube dimensions: total volume, length and internal diameter of the reformer tube, number of tubes, wall thickness
-
Tube packing: void fraction
-
Catalyst data: diameter, sphericity, solid density, solid heat capacity, number of holes, tortuosity, mean pore radius, catalyst characteristic length, catalyst support
-
Inlet process composition: flow rate, natural gas composition, pressure, temperature
-
Outside tube wall temperature: measured values
-
Heat transfer coefficient
-
Activity coefficient.
