In ammonia synthesis, fixed bed reactors are always used. Only one reaction occurs in the gas phase using iron catalyst particles   , as described in Equation (1).
The reaction is exothermic in pressures from 150 to 300 atm, due to the decrease in mole numbers. Moreover, the temperature range is 600 to 800 K. After all, even with an exothermic reaction, the reaction rate in converters must be high. The High Pressure and High Temperature (HPHT) conditions in gas phase make properties in reactive fluid temperature vary too much along the converter. Some examples are density, fugacity, chemical activity and heat capacity.
The mathematical modeling of ammonia converters uses experimental rate expressions. They generally contain partial pressures  or chemical activities  in formulation. One of the first models was presented by the authors Emmett and Kummer  . They used Temkin and Pyzhev expression  to obtain experimental and predicted data. Moreover, they proposed the formulation of a reaction rate using fugacity and not partial pressures.
After a few years, many researchers continued to use Temkin and Pyzhev rate in their calculations     . However, in 1964, Nielsen and collaborators proposed modifications in reaction rate  . They verified that chemical activities could be included in the rate model.
Dyson and Simon  also used chemical activities in their work. They also proposed a kinetic correction factor to a heterogeneous reaction. Many studies continued to use Dyson and Simon expression, with small modifications. Many presented good validation with plant data and were focused on optimization analysis  -  .
However, in previous publications, the fugacity coefficients were derived from a correlation that they varied with temperature and pressure, but not composition  -  . The implications are that the chemical activity does not change with molar and mass fraction variations inside the reactor. Moreover, the mass, energy and momentum balances provided remarkable composition variations, even with only one reaction. This effect is more pronounced when fluid presents differences in molar weight (which is the case of ammonia production).
Therefore, a more sophisticated model can be used to provide a multicomponent description. This method is the principle of the compositional approach: the reactant fluid present changes of properties influenced not only by pressure and temperature but also by composition and intermolecular forces. The last and more important part is computed due to mixing rules in cubic EoS calculations. Furthermore, two facts reinforce the use of a compositional approach in ammonia synthesis: 1) the molecules in the gas phase are smaller and 2) the high-pressure deviates gas from ideal treatment, approximating gas molecules.
As most of substances are non-polar (except for NH3) and gases are at HPHT conditions, the Peng-Robinson (PR) and Soave-Redlich-Kwong (SRK) cubic EoS are chosen to model the system. Therefore, the EoS approach will replace the previous activities used in reaction rate. As discussed above, these models present a more reliable theory and a molecular formulation.
The main objectives of this work are: validate of two ammonia reactor models with plant data (adiabatic and autothermal); use a fitted reaction rate using a compositional approach; solve mass, energy and momentum balances in PFR model in steady-state and present a study of parametric sensitivity in input variables of ammonia converters.
2. Mathematical Background
2.1. Thermodynamic Modeling
The SRK  and PR  cubic EoS are the two most widely used cubic equations in the process industry. According to Ahmed  , the general form of two-parameter cubic EoS is shown at Equation (2).
For the results of this work, only PR-EoS was used, with the parameters and . Defining the following Equations (3) to (5) below:
and substituting into Equation (2), as described at Barbosa Neto  , the implicit form of the cubic EoS, in the compressibility factor Z, is obtained:
The van der Waals one-fluid mixing rules are used for the energy A, and for the volume B, parameters of the EoS, are expressed by Equations (7) and (8), respectively.
The parameter Aij, in the Equation (7), is calculated from Equation (9).
The terms Ai and Bi are defined by Equations (10) and (11), respectively. The factor m(ωi) is given in the literature  . Values of Ωa and Ωb changes with EOS.
The expression shown at Equation (12) for fugacity coefficients is obtained.
The fugacity and activity are also important for reaction rate used. They are given in Equations (14) and (15).
Other important properties are the heat capacities, which are given in Equations (16) and (17). More details are found in Tosun  .
2.2. Rate Expression
The rate used in ammonia reactors of this research was made by Singh and Saraf  , as given in Equation (18). We compute the chemical activities predicted by a suited thermodynamic model. The previous works presented fugacity coefficients given correlations  .
The rate given in Equation (18) is pseudo-homogeneous. Nielsen and collaborators provided a suitable range for ammonia rate expressions: 640 to 770 K and 150 to 310 atm  . Therefore, it was corrected by an effectiveness factor η for a heterogeneous reactor, as shown in Equation (19) and Table 1.
To summarize this section, the equilibrium constant (Keq) in Equation (19) is found in correlation given by Gillespie and Beattie  .
Table 1. Coefficients for Equation (19)  .
2.3. Mass, Momentum and Energy Balances (Adiabatic Reactors)
In the adiabatic operation of ammonia reactors, the volume of control increases its temperature using only the heat of reaction. However, this operation is not possible in one converter. The reactor volume must be separated into several reactors, as given in Figure 1.
The mass balance in one reactor is expressed only regarding nitrogen conversion (Equations (20) and (21)) because only one reaction takes place in converter  . The only difference between the following relations is the area of reactor. Therefore, the coordinates can be done in length L or volume V.
Equations (22) and (23) express the energy balances. As the reactor is at high pressure, the estimation of Cp is essential. Moreover, the heat of reaction is computed according to correlation  .
Besides, as reactor contains catalyst particles, there is a pressure loss along the converter. It is estimated using the Ergun Equation  , as expressed in Equation (24). However, this relation is only for dP computations in length L. For dV computations, the pressure loss dP/dV is set to zero  .
2.4. Mass, Momentum and Energy Balances (Autothermal Reactors)
The autothermal reactor uses the energy released by a reaction to heat the reactant gas. It operates with countercurrent flow, and it is divided into several tubes, as given in Figure 2. The catalyst is inside the tubes, while the cooling gas increases its temperature along the converter.
Figure 1. Indirect cooling adiabatic ammonia reactor  .
Figure 2. Scheme of an autothermal reactor  .
The mass balance and pressure loss in the autothermal converter do not change compared to the adiabatic model. The difference is in the temperature: now we have the reactant gas (that increases its temperature by the reaction and is cooling by cooling gas) and the cooling gas (that always increases its temperature). The energy balances are given in Equations (25) and (26) (for length variations) and Equations (27) and (28) (for volume variations). The minus sign in Equations (26) and (28) is related to countercurrent operation.
3.1. Solution of ODEs System
The adiabatic and autothermal reactors must be solved using a numerical method, once an analytical solution is complex. The method used was the Runge-Kutta-Fehlberg(RKF) using an error control strategy, described by Chapra and Canale  .
Another challenge is the interdependence of variables in differential equations. Therefore, the code was made using a modular structure in Wolfram Mathematica® Programming Language. Moreover, our algorithm is called MARS (Models for Ammonia Reactor Simulation).
The variables of interest in both models are the outlet temperature, pressure, and composition. The thermodynamic module computes the thermodynamic and transport properties; the kinetic module calculates the reaction rate and another module joins all the previous in a numerical method to solve the balances. As the problem is solved in one dimension, the stopping-criteria is the end of the reactor, as expressed in Figure 3.
3.2. Variation of α in Reaction Rate
As discussed before, the chemical activities originally are computed using a correlation―Singh and Saraf Rate (Equation (18)). However, when using the multicomponent approach, it is expected that chemical activity decreases, due to compositional interactions. Therefore, the reaction rate computed also decreases its value. So, the kinetic factor α is fitted in MARS. The fit is made according to an adiabatic reactor in the literature  . Only the first bed is calculated in temperature and conversion. Table 2 shows the parameters for this reactor.
Figure 4 presents the N2 conversion variation in an adiabatic converter with α alterations.
Figure 3. Flowchart for calculation of reactor.
Table 2. Input parameters and plant data for the 1st reactor  .
Figure 4. Conversion profile in first bed using α variations at MARS (EoS approach).
The outlet conversion increases when α rises because the reaction rate is higher. Moreover, at α = 0.570, the model does well in predicting the outlet conversion. Therefore, the original value of α = 0.550 in Singh and Saraf rate  should be replaced by α = 0.570 when using the compositional approach from now on. The conversion was chosen because it is the variables which have more effect on output composition and its error is more significant than temperature.
4. Results and Discussions
For validations in adiabatic and autothermal models, the relative error was calculated as expressed in Equation (29). In the case of multiple reactors, the maximum relative error will be taken.
Even with a suitable numerical method and a robust code, validations of reactors models are necessary. The RKF method and α = 0.570 are selected for Singh and Saraf modified rate. Furthermore, calculates an adiabatic reactor containing three fixed beds in series and an autothermal converter. Both models are reliable compared to plant data.
4.1. Model Validation
4.1.1. Adiabatic Reactor
In adiabatic case, we have three reactors in series. The first bed is computed in a variation of α. Therefore, all inlets parameters remain the same. The only additional information for simulation is the inlet temperature of the second and third reactors and their respective volume, as given in Table 3.
Figure 5 shows the temperature profile obtained. In this, the highest errors in temperature are noted in the first reactor. Moreover, the first converter is the
Figure 5. Plant data and MARS EoS model temperature profile (adiabatic reactor).
Table 3. Plant data for three adiabatic beds in series  .
place where the rate has its highest values. Therefore, it can provide more errors compared to the others. However, in the second and third converters, the simulated temperature gives good results compared to plant data. In all simulations of the adiabatic arrangement, 74 iterations are required with RKF. It is a good result compared if the 4th Order Runge-Kutta with fixed step size was used (with 50 iterations at each reactor, for example).
In addition, Table 4 presents the relative errors of temperature. It proves a good comparison between MARS and plant data (Singh and Saraf, 1979)  . The maximum error of temperature reached by other authors was less than 6% (Singh and Saraf, 1979)  .
Even with good results in temperature predictions, the conversion is another crucial variable. The composition of the reactor depends on conversion, after all. Furthermore, it is more sensitive to variations than temperature. Figure 6 shows the conversion profile computed by MARS.
In Figure 6, the highest error in conversion occurred in the third reactor. The first and second converters presented a good agreement with plant data. The main error in the third reactor is the previous errors inherited by the first and second reactors (the error was propagated). Moreover, the final converter is where the reaction rate presents the smallest value. The relative errors in conversion are summarized in Table 5. To summarize, even with the difference in conversion in the 3rd reactor, the MARS model is reliable for adiabatic reactor
Figure 6. Plant data and MARS EoS model conversion profile (adiabatic reactor).
Table 4. Comparison between outlet temperature in plant data and the MARS model.
Table 5. Comparison between outlet conversion in plant data and the MARS model.
simulation, because the temperature is usually the control variable in ammonia reactors. Errors in literature reached less than 0.5%  .
4.1.2. Autothermal Reactor
The autothermal converter has the reactor parameters given in Table 6. In this type of converter, the temperature is measured in several points of the reactor.
In Figure 7, differences are noted between simulation and plant data for reactant gas temperature inside reactor. First, only at the end of the autothermal reactor are the differences significant. These occur due to a change of dynamics in the ODE system.
However, the maximum temperature point is well predicted, which reinforces the method’s effectiveness. The maximum relative error was 2.7% in the end of reactor, due to high nonlinearity of equations.
Table 7 presents the relative errors at each point of the autothermal reactor. It
Figure 7. Comparison between the MARS EoS model and plant data (autothermal reactor).
Table 6. Input parameters for autothermal reactor  .
Table 7. Relative errors in reactant gas temperature in autothermal reactor  .
proves the reliability of the MARS model for the autothermal reactor. Other authors reached 3.2% of maximum relative difference in temperature, which was an acceptable value  .
4.2. Sensitivity Analysis
In the sensitivity analysis the same data set of validation was used to fit α. In the adiabatic model, the analysis is realized in three reactors in series. The inlet temperature was varied. Then, the effects on temperature, conversion and effectiveness factor profiles were detected in each case. While, in the autothermal model, the heat exchange coefficient was varied, with the same detections on output of reactor. The variations were computed in length L.
4.2.1. Adiabatic Reactor
Figure 8 presents the adiabatic reactors flowchart and dimensions. In all the analysis, the reactor consists of 3 beds, with indirect cooling. The values of volume and length were estimated using literature   .
In adiabatic operation, only the inlet temperature of the 1st reactor is varied. The temperatures Tin2 and Tin3 maintain the same values (as given in Table 8). Moreover, the RKF method with variable step size is used for numerical simulation.
The inlet temperature effect along converter temperature is given in Figure 9. Analyzing the Figure 9, we can note that the lowest inlet temperature (Tin = 643.15 K) gave the highest profiles along second and third reactors. After all, smaller temperatures provide smaller rates, which give minor conversions. Therefore, after first reactor, more N2 can be converted, giving a more accentuated profile (Tin = 643.15 K).
On the other hand, high inlet temperatures present high reaction rates, converting more N2 in the first bed. However, it becomes more difficult to react in the second and third converters. Therefore, the temperature rise is not so significant as in smaller Tin values, even with high rates. In the three cases, the number of iterations during the converging process is similar.
Figure 8. Flowchart for adiabatic beds parametric sensitivity.
Table 8. Parameters for inlet temperature variation in adiabatic model (643.15 K, 663.15 K and 683.15 K).
Figure 9. Temperature profiles in beds (Tin = 643.15 K, 663.15 K and 683.15 K).
Another essential measure inside our reactor is the composition. As there is only one reaction, the conversion profile gives the other molar fractions indirectly. For the inlet temperature variation case, Figure 10 shows the conversion profile.
In Figure 10, it is noted that larger inlet temperatures give higher N2 conversions. However, a high Tin does not guarantee to elevate conversions, because the reaction is exothermic. After all, the reaction is reversible, and high temperatures provide low equilibrium conversion. Moreover, the highest rise in conversions occurs in the first bed. It can be seen that the growth of conversion between interval of 683.15 K and 663.15 K is smaller than 663.15 K and 643.15 K. Therefore, there is a limit to the Tin value.
Figure 11 shows the effectiveness factor profiles. As inlet temperature increases, η factor also increases. This occurs due to a higher reaction rate. Moreover, the first reactor presents the smaller values of η. After all, the overall conversion is lower in first reactor.
4.2.2. Autothermal Reactor
Figure 12 presents the autothermal reactor flowchart and dimensions.
The parameters for simulation are given in Table 9.
Figure 13 presents the reactant gas temperatures profiles. In this, we note that the length in the reactor which presents the largest temperature does not change. However, the maximum temperature value for U = 450 W/m2∙K is 753 K, while for U = 650 W/m2∙K, is 739 K. Moreover, another important value is the final value of temperature for the reactant gas at the end of the reactor. For U = 450 W/m2∙K, the final T value is 724 K and U = 650 W/m2∙K is 687 K. It takes place because a large U removes more energy in the reactor, decreasing temperature more rapidly.
Figure 10. Conversion profiles in beds (Tin = 643.15 K, 663.15 K and 683.15 K).
Figure 11. Effectiveness factor profiles in beds (Tin = 643.15 K, 663.15 K and 683.15 K).
Table 9. Parameters for inlet temperature variation in autothermal model (450, 650 and 850 W/m2∙K).
Figure 14 presents the conversion profile. When U increases, the final conversion also increases. It occurs because more energy is removed by the reaction, decreasing the speed of the reverse reaction. However, U cannot be raised too much, otherwise the rate would be very low.
The effectiveness factor η profiles are given in Figure 15. As U increases, the temperature decreases in the reactor, raising equilibrium conversion and η values. However, in U = 850 W/m2∙K, there are errors computing η (higher than 1).
Figure 12. Adiabatic reactor flowchart and dimensions.
Figure 13. Reactant gas temperature profiles in autothermal reactor (U = 450, 650 and 850 W/m2∙K).
Figure 14. Conversion profiles in autothermal reactor (U = 450, 650 and 850 W/m2∙K).
Figure 15. Effectiveness factor profiles in autothermal reactor (U = 450, 650 and 850 W/m2∙K).
The compositional modeling using PR cubic EoS was essential to calculate the properties in reactive streams at HPHT conditions during the ammonia synthesis reactors simulation. The adiabatic reactor model presented a maximum relative error of 1.6% in temperature and 11.4% in conversion when compared to plant data. In sensitivity study, a value of Tin = 683.15 K gave the highest conversion of 25.16%. The autothermal reactor model presented a maximum error of 2.7% in temperature when compared to experimental points. In parametrical sensitivity, the highest conversion of 37.22% was provided by a value of U = 850 W/m2∙K. Therefore, both models proved reliable in simulating ammonia reactors for the set of data used.
Another improvement of the ammonia synthesis reactor models can be achieved using intraparticle diffusional approach. It computes the effectiveness factor without using an experimental correlation.
The authors acknowledge CNPq (Process Number 131744/2016-0) and Unicamp for providing financial support.
List of Abbreviations and Acronyms
HPHT High Pressure and High Temperature
PR Peng Robinson
SRK Soave Redlich Kwong
BWR Benedict Webb Rubin
EoS Equation of State
PFR Plug Flow Reactor
i Substance i
List of Symbols
List of Roman Letters
Ammonia production rate
P (Pa) System pressure in reactor
T (K) System temperature
Ideal gas constant
Equilibrium constant for ammonia rate
Molar fraction of substance i in reactor
Molar volume of gaseous system
Covolume term for cubic EoS
Function of temperature in cubic EoS
A Mixture term for cubic EoS computing interactions
B Mixture term for cubic EoS computing interactions
Mixture term for cubic EoS computing interactions
Z Mixture compressibility factor
Binary interaction parameter
Reduced temperature of component i
Reduced pressure of component i
Fugacity of i component in gaseous mixture
Chemical activity of i component in gaseous mixture
Gibbs excess energy for mixture
Mixture heat capacity in mass units
Real heat capacity at constant pressure of mixture
Residual heat capacity at constant pressure of mixture
Ideal gas heat capacity at constant pressure of mixture
Real heat capacity at constant volume of mixture
Residual heat capacity at constant volume of mixture
Ideal gas heat capacity at constant volume of mixture
Experimental coefficients for calculation
Molar fraction of substance i in reactor
Heat of reaction for ammonia synthesis
Initial molar flow rate of nitrogen in reactor
Mass flow inside the reactor
Sectional area of the reactor
Reactant gas temperature in reactor models
Cooling gas temperature in autothermal reactor
Overall heat transfer coefficient in autothermal converter
Specific heat exchange area in length
Specific heat exchange area in volume
Superficial velocity of gas in reactor
or m3 Step size along the reactor in numerical method
List of Greek Letters
Molar density of gaseous system
Constant for cubic EoS (PR or SRK)
Constant for cubic EoS (PR or SRK)
Constant for cubic EoS (PR or SRK)
Constant for cubic EoS (PR or SRK)
Constant for cubic EoS (PR or SRK)
Mixture factor in cubic EoS for i component
Acentric factor for i component
Acentric factor function
Fugacity coefficient in gas phase
Effectiveness factor inside the catalyst pellet
(Pa∙s) Gas viscosity