Discrete Simulation Model of Industrial Natural Gas Primary Reformer in Ammonia Production and Related Evaluation of the Catalyst Performance

The process described herein is based on the Kellogg Inc catalytic high-pressure reforming method for producing ammonia starting with natural gas feed. An ammonia plant steam reforming unit can produce 1360 tonnes per day of liquid ammonia. Figure 2 presents the steady-state flow sheet of the SMR unit build in UniSim^® Design R470 with the main process flow designated with the red line.

Fig. 2.

Natural gas feed at a pressure of about 32 bar enters the natural gas knock-out drum 120-F for elimination of entrained liquid. The outlet line of 120-F feeds the one-stage centrifugal natural gas feed compressor 102-J driven by back-pressure (40/4 bar) steam turbine 102-JT. Outlet pressure of natural gas is at the level of 42 bar. Hydrogen required for desulfurisation of the natural gas is injected into the paralleled natural gas stream entering the natural gas fired heater 103-B. The outlet temperature of 103-B is 400°C. The heated natural gas stream flows through two reactors in series. The first is the hydrogenator 101-D, which contains a single bed of cobalt-molybdenum catalyst. It converts the organic sulfur compounds to hydrogen sulfide in the presence of the hydrogen injected upstream of 103-B. The natural gas stream next passes into the desulfuriser reactor 102-D, which contains a single bed of zinc-oxide catalyst. In this reactor the hydrogen sulfide is converted to zinc sulfide which remains in the catalyst.

The desulfurised natural gas, plus residual hydrogen, leaves 102-D with a sulfur content of 0.25 ppm and at a temperature of 370°C. The natural gas plus residual hydrogen stream is joined by the process steam in the mixer. The process steam is at a pressure of about 40 bar and a temperature of 392°C. The steam flow is controlled with the steam-to-natural gas (S/NG) molar ratio controller.

The SMR feed gas flows to the mixed feed coil, which is located in the convection section of the SMR furnace. In this coil, the SMR feed is heated to about 510°C. After heating, the SMR feed flows down through ten rows of reformer tubes that are suspended in the radiant box of primary reformer 101-B. Eleven rows of forced draught down fired burners are located in parallel rows to the catalyst tubes, in total 198 burners. They raise the feed temperature to about 790°C at the outlet of the catalyst tubes. In addition, 11 tunnel burners are used to heat the waste gases passing from the radiant to the convection part of the SMR furnace. 520 catalyst tubes with a total length of 10 m and inside diameter of 0.0857 m contain 30 m³ of nickel reformer catalyst. The reformed gas (syngas) then flows to the secondary reformer for further processing.

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 k_i 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 k₂ and k₄ (Equation (ii) and (iv)) are kmol bar^0.5 kg^–1_cat h^–1), while the unit of k₃ (Equation (iii) is kmol bar^–1 kg^–1_cat h^–1).

Table I (20) gives the parameters for the activation energies, E_i, and for the pre-exponential factors, A_i, used in the model, valid for most of the commercial nickel catalysts with either MgAl₂O₄ or CaAl₁₂O₁₉ support.

Apparent adsorption equilibrium constants K_i 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:

B_i is the pre-exponential factor expressed in bar^-1 or unitless, while ΔH_i 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 MgAl₂O₄ or CaAl₁₂O₁₉ 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.

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 m²; ρ_B represents the catalyst bed density in kg m^–3; F_i 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; D_e,A is the species effective diffusivity in m³_fluid m^–1_catalyst 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 kg_catalyst m^–3_catalyst.

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 m³_fluid m^–3_catalyst; τ denotes the catalyst tortuosity and is the average diffusivity of species A.

The average diffusivity of species is determined by Equation (xix):

where D_A is the diffusivity of the reacting species A given by Equation (xx) and S(r_p,i) is the void fraction taken by the pores with radii ranging from r_p,i to r_{p,i +1}:

where D_kA is the Knudsen diffusivity in m³_fluid m^–1_catalyst 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; d_p 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 m² h K kJ^–1; T_t,0 is the temperature of the furnace that surrounds the reformer tubes; ΔH_i 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 r_i 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 r_jdV 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: