Modeling and Numerical Simulation of Ammonia Synthesis Reactors Using Compositional Approach

Diogo Silva Sanches Jorqueira^1*, Antonio Marinho Barbosa Neto², Maria Teresa Moreira Rodrigues¹

¹Department of Chemical Systems Engineering, Campinas University, Campinas, Brazil.
²Simworx, Engineering, Research & Development, Campinas, Brazil.

DOI: 10.4236/aces.2018.83009                

Abstract

Ammonia synthesis reactors operate in conditions of high pressure and high temperature. Consequently, the flow inside these reactors always presents interaction between components in the feed mixture. A modeling accounts these interactions with pressure, temperature and the molar fraction is essential to converter simulation more realistic. The compositional approach based on cubic equations of state provides the influences of the component of a gas mixture using mixing rules and binary interaction parameters. This multicomponent description makes the model more robust and reliable for properties mixture prediction. In this work, two models of ammonia synthesis reactors were simulated: adiabatic and autothermal. The fitted expression of Singh and Saraf was used. The adiabatic reactor model presented a maximum relative error of 1.6% in temperature and 11.4% in conversion while the autothermal reactor model presents a maximum error of 2.7% in temperature, when compared to plant data. Furthermore, a sensitivity analysis in input variables of both converter models was performed to predict operational limits and performance of the Models for Ammonia Reactor Simulation (MARS).

Keywords

Cubic Equation of State, Adiabatic Converter, Autothermal Converter, Temkin-Pyzhev Expression

1. Introduction

In ammonia synthesis, fixed bed reactors are always used. Only one reaction occurs in the gas phase using iron catalyst particles [1] [2] , as described in Equation (1).

${\text{N}}_{\text{2}}+{\text{3H}}_{\text{2}}\leftrightarrow {\text{2NH}}_{\text{3}}\text{,Δ}{H}_r^o=-\text{45}\text{.94}\text{kJ}/\text{mol}\text{of}{\text{NH}}_{\text{3}}$ (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 [1] or chemical activities [2] in formulation. One of the first models was presented by the authors Emmett and Kummer [3] . They used Temkin and Pyzhev expression [1] 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 [4] [5] [6] [7] . However, in 1964, Nielsen and collaborators proposed modifications in reaction rate [8] . They verified that chemical activities could be included in the rate model.

Dyson and Simon [2] 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 [9] - [15] .

However, in previous publications, the fugacity coefficients were derived from a correlation that they varied with temperature and pressure, but not composition [9] - [15] . 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 NH₃) 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 [16] and PR [17] cubic EoS are the two most widely used cubic equations in the process industry. According to Ahmed [18] , the general form of two-parameter cubic EoS is shown at Equation (2).

$p=\frac{RT}{v-b}-\frac{a}{\left(v+{\delta }_{1}b\right)\left(v+{\delta }_{2}b\right)}$ (2)

For the results of this work, only PR-EoS was used, with the parameters ${\delta }_1=1+2^{0.5}$ and ${\delta }_2=1+2^{0.5}$ . Defining the following Equations (3) to (5) below:

$A=\frac{ap}{{\left(RT\right)}^2}$ (3)

$b=\frac{bp}{RT}$ (4)

$Z=\frac{pv}{RT}$ (5)

and substituting into Equation (2), as described at Barbosa Neto [19] , the implicit form of the cubic EoS, in the compressibility factor Z, is obtained:

$\begin{array}{l}{Z}^3+\left[\left({\delta }_1+{\delta }_2-1\right)B-1\right]Z^2+\left[A+{\delta }_1{\delta }_2{B}^2-\left({\delta }_1+{\delta }_2\right)B\left(B+1\right)\right]Z\\ -\left[AB+{\delta }_1{\delta }_2{B}^2\left(B+1\right)\right]=0\end{array}$ (6)

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.

$A={\displaystyle \underset{i=1}{\overset{nc}{\sum }}{\displaystyle \underset{j=1}{\overset{nc}{\sum }}{y}_i{y}_j{A}_{ij}}}$ (7)

$B={\displaystyle \underset{i=1}{\overset{nc}{\sum }}{y}_i{B}_i}$ (8)

The parameter A_ij, in the Equation (7), is calculated from Equation (9).

$A_{ij}=\sqrt{A_i{A}_j}\left(1-k_{ij}\right)$ (9)

The terms A_i and B_i are defined by Equations (10) and (11), respectively. The factor m(ω_i) is given in the literature [17] . Values of Ω_a and Ω_b changes with EOS.

$A_i={\Omega }_a\left(\frac{p_{ri}}{T_{ri}^2}\right){\left[1+m\left({\omega }_i\right)\left(1-\sqrt{T_{ri}}\right)\right]}^2$ (10)

$B_i={\Omega }_b\left(\frac{p_{ri}}{T_{ri}}\right)$ (11)

The expression shown at Equation (12) for fugacity coefficients is obtained.

$\ln \left({\phi }_i\right)=\left(Z-1\right)\left(\frac{B_i}{B}\right)-\ln \left(Z-B\right)-\left(\frac{A}{\Delta B}\right)\left(\frac{2{\psi }_i}{A}-\frac{B_i}{B}\right)\ln \left[\frac{Z+{\delta }_{1}B}{Z+{\delta }_{2}B}\right]$ (12)

With $\Delta ={\delta }_1-{\delta }_2$ and

${\psi }_i={\displaystyle \underset{i=1}{\overset{nc}{\sum }}{A}_i{}_j{y}_j}$ (13)

The fugacity and activity are also important for reaction rate used. They are given in Equations (14) and (15).

$f_i={\phi }_i{y}_{i}P$ (14)

$a_i=\frac{{\phi }_i{y}_{i}P}{P_{ref}}=\frac{f_i}{P_{ref}}$ (15)

Other important properties are the heat capacities, which are given in Equations (16) and (17). More details are found in Tosun [20] .

$C_p=C_p^{ig}+C_p^{res}$ (16)

$C_v=C_v^{ig}+C_v^{res}$ (17)

2.2. Rate Expression

The rate used in ammonia reactors of this research was made by Singh and Saraf [9] , as given in Equation (18). We compute the chemical activities predicted by a suited thermodynamic model. The previous works presented fugacity coefficients given correlations [2] .

$r_{{\text{NH}}_{\text{3}}}=4.11\times {10}^{10}\exp \left(\frac{-163422}{RT}\right)\left[K_{eq}^2{a}_{{\text{N}}_{\text{2}}}{\left(\frac{a_{{\text{H}}_{\text{2}}}^3}{a_{{\text{NH}}_{\text{3}}}^2}\right)}^{\alpha }-{\left(\frac{a_{{\text{NH}}_{\text{3}}}^2}{a_{{\text{H}}_{\text{2}}}^3}\right)}^{1-\alpha }\right]$ (18)

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 [8] . Therefore, it was corrected by an effectiveness factor η for a heterogeneous reactor, as shown in Equation (19) and Table 1.

$\eta =b_0+b_{1}T+b_2{x}_{{\text{N}}_{\text{2}}}+b_3{T}^2+b_4{x}_{{\text{N}}_{\text{2}}}+b_5{T}^3+b_6{x}_{{\text{N}}_{\text{2}}}^3$ (19)

To summarize this section, the equilibrium constant (K_eq) in Equation (19) is found in correlation given by Gillespie and Beattie [21] .

Table 1. Coefficients for Equation (19) [2] .

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 [2] . The only difference between the following relations is the area of reactor. Therefore, the coordinates can be done in length L or volume V.

$\frac{\text{d}{x}_{{\text{N}}_{\text{2}}}}{\text{d}L}=\frac{A{r}_{{\text{NH}}_{\text{3}}}\eta }{2F_{{\text{N}}_{\text{2}}}^0}$ (20)

$\frac{\text{d}{x}_{{\text{N}}_{\text{2}}}}{\text{d}V}=\frac{r_{{\text{NH}}_{\text{3}}}\eta }{2F_{{\text{N}}_{\text{2}}}^0}$ (21)

Equations (22) and (23) express the energy balances. As the reactor is at high pressure, the estimation of C_p is essential. Moreover, the heat of reaction is computed according to correlation [21] .

$\frac{\text{d}T}{\text{d}L}=\frac{A{r}_{{\text{NH}}_{\text{3}}}\eta \left(-\Delta H_r\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}$ (22)

$\frac{\text{d}T}{\text{d}V}=\frac{r_{{\text{NH}}_{\text{3}}}\eta \left(-\Delta H_r\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}$ (23)

Besides, as reactor contains catalyst particles, there is a pressure loss along the converter. It is estimated using the Ergun Equation [22] , 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 [9] .

$\frac{\text{d}P}{\text{d}L}=\frac{150{\left(1-\epsilon \right)}^2\mu u}{{\epsilon }^3{d}_p^2}-\frac{1.75\left(1-\epsilon \right)\rho u^2}{{\epsilon }^3{d}_p}$ (24)

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 [12] .

Figure 2. Scheme of an autothermal reactor [23] .

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.

$\frac{\text{d}T}{\text{d}L}=\frac{A{r}_{{\text{NH}}_{\text{3}}}\eta \left(-\Delta H_r\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}-\frac{U{A}^{\prime }\left(T-T_g\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}$ (25)

$\frac{\text{d}{T}_g}{\text{d}L}=-\frac{U{A}^{\prime }\left(T-T_g\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}$ (26)

$\frac{\text{d}T}{\text{d}V}=\frac{r_{{\text{NH}}_{\text{3}}}\eta \left(-\Delta H_r\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}-\frac{U{a}^{\prime }\left(T-T_g\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}$ (27)

$\frac{\text{d}{T}_g}{\text{d}V}=-\frac{U{a}^{\prime }\left(T-T_g\right)}{\stackrel{\dot{}}{m}{C}_{p_{mix}}}$ (28)

3. Methodology

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 [24] .

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 [9] . Only the first bed is calculated in temperature and conversion. Table 2 shows the parameters for this reactor.

Figure 4 presents the N₂ 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 [9] .

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 [9] 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.

${\text{rel}}_{\text{error}}=100\frac{\left|{\text{plant}}_{\text{data}}-{\text{simulation}}_{\text{data}}\right|}{{\text{plant}}_{\text{data}}}$ (29)

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 [9] .

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 4^th 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) [9] . The maximum error of temperature reached by other authors was less than 6% (Singh and Saraf, 1979) [9] .

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 3^rd 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% [9] .

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 [9] .

Table 7. Relative errors in reactant gas temperature in autothermal reactor [9] .

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 [9] .

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 [9] [12] .

In adiabatic operation, only the inlet temperature of the 1^st reactor is varied. The temperatures T_in2 and T_in3 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 (T_in = 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 N₂ can be converted, giving a more accentuated profile (T_in = 643.15 K).

On the other hand, high inlet temperatures present high reaction rates, converting more N₂ 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 T_in 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 (T_in = 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 N₂ 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 T_in 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/m²∙K is 753 K, while for U = 650 W/m²∙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/m²∙K, the final T value is 724 K and U = 650 W/m²∙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 (T_in = 643.15 K, 663.15 K and 683.15 K).

Figure 11. Effectiveness factor profiles in beds (T_in = 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/m²∙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/m²∙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/m²∙K).

Figure 14. Conversion profiles in autothermal reactor (U = 450, 650 and 850 W/m²∙K).

Figure 15. Effectiveness factor profiles in autothermal reactor (U = 450, 650 and 850 W/m²∙K).

5. Conclusions

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 T_in = 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/m²∙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.

Acknowledgements

The authors acknowledge CNPq (Process Number 131744/2016-0) and Unicamp for providing financial support.

Nomenclature

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

N₂ Nitrogen

H₂ Hydrogen

NH₃ Ammonia

CH₄ Methane

Ar Argon

List of Roman Letters

$r_{{\text{NH}}_{\text{3}}}\text{kmol}/\left({\text{m}}^3\cdot \text{s}\right)$ Ammonia production rate

P (Pa) System pressure in reactor

T (K) System temperature

$R\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Ideal gas constant

$K_{eq}$ Equilibrium constant for ammonia rate

$y_i$ Molar fraction of substance i in reactor

$v{\text{m}}^{\text{3}}/\text{mol}$ Molar volume of gaseous system

$b{\text{m}}^{\text{3}}/\text{mol}$ Covolume term for cubic EoS

$a\left(\text{Pa}\cdot {\text{m}}^{\text{6}}\right)/{\text{mol}}^{\text{2}}$ Function of temperature in cubic EoS

A Mixture term for cubic EoS computing interactions

B Mixture term for cubic EoS computing interactions

$A_{ij}$ Mixture term for cubic EoS computing interactions

Z Mixture compressibility factor

$k_{ij}$ Binary interaction parameter

$T_{r_i}$ Reduced temperature of component i

$P_{r_i}$ Reduced pressure of component i

$P_{ref}\text{Pa}$ Reference pressure

$f_i\text{Pa}$ Fugacity of i component in gaseous mixture

$a_i$ Chemical activity of i component in gaseous mixture

$g_e\text{J}/\text{mol}$ Gibbs excess energy for mixture

$C_{p_{mix}}\text{J}/\left(\text{kg}\cdot \text{K}\right)$ Mixture heat capacity in mass units

$C_{p_{mix}}\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Real heat capacity at constant pressure of mixture

$C_p^{res}\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Residual heat capacity at constant pressure of mixture

$C_p^{ig}\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Ideal gas heat capacity at constant pressure of mixture

$C_v\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Real heat capacity at constant volume of mixture

$C_v^{res}\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Residual heat capacity at constant volume of mixture

$C_v^{ig}\text{J}/\left(\text{mol}\cdot \text{K}\right)$ Ideal gas heat capacity at constant volume of mixture

$b_0\text{-}{b}_6$ Experimental coefficients for $\eta$ calculation

$x_{{\text{N}}_{\text{2}}}$ Nitrogen conversion

$L\text{m}$ Reactor length

$V{\text{m}}^3$ Reactor volume

$y_i$ Molar fraction of substance i in reactor

$\Delta H_r\text{J}/\text{mol}$ Heat of reaction for ammonia synthesis

$F_{{\text{N}}_{\text{2}}}^o\text{mol}/\text{s}$ Initial molar flow rate of nitrogen in reactor

$\stackrel{\dot{}}{m}\text{kg}/\text{s}$ Mass flow inside the reactor

$A{\text{m}}^2$ Sectional area of the reactor

$T\text{K}$ Reactant gas temperature in reactor models

$T_g\text{K}$ Cooling gas temperature in autothermal reactor

$U\text{W}/\left({\text{m}}^{\text{2}}\cdot \text{K}\right)$ Overall heat transfer coefficient in autothermal converter

$A^{\prime }{\text{m}}^{\text{2}}/\text{m}$ Specific heat exchange area in length

$a^{\prime }{\text{m}}^{\text{2}}/{\text{m}}^{\text{3}}$ Specific heat exchange area in volume

$u\text{m}/\text{s}$ Superficial velocity of gas in reactor

$d_p\text{m}$ Particle diameter

$h\text{m}$ or m³ Step size along the reactor in numerical method

List of Greek Letters

$\rho \text{mol}/{\text{m}}^3$ Molar density of gaseous system

$\Delta$ Constant for cubic EoS (PR or SRK)

${\delta }_1$ Constant for cubic EoS (PR or SRK)

${\delta }_2$ Constant for cubic EoS (PR or SRK)

${\Omega }_a$ Constant for cubic EoS (PR or SRK)

${\Omega }_b$ Constant for cubic EoS (PR or SRK)

${\psi }_i$ Mixture factor in cubic EoS for i component

${\omega }_i$ Acentric factor for i component

$m\left({\omega }_i\right)$ Acentric factor function

${\phi }_i$ Fugacity coefficient in gas phase

$\eta$ Effectiveness factor inside the catalyst pellet

$\epsilon$ Bed porosity

$\mu$ (Pa∙s) Gas viscosity

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Temkin, M. and Pyzhev, V. (1939) Kinetics of the Synthesis of Ammonia on Promoted Iron Catalysts. Journal of Physical Chemistry, 13, 851.
[2]	Dyson, D.C. and Simon, J.C. (1968) A Kinetic Expression with Diffusion Correction for Ammonia Synthesis on Industrial Catalyst. Industrial and Engineering Chemistry Fundamentals, 7, 605-610. s://doi.org/10.1021/i160028a013
[3]	Emmett, P.H. and Kummer, J.T. (1943) Kinetics of Ammonia Synthesis. Industrial and Engineering Chemistry, 35, 677-683. s://doi.org/10.1021/ie50402a012
[4]	Annable, D. (1952) Application of the Temkin Kinetic Equation to Ammonia Synthesis in Large-Scale Reactors. Chemical Engineering Science, 1, 145-154. s://doi.org/10.1016/0009-2509(52)87011-3
[5]	Baddour, R.F., Brian, P.L.T., Logeais, B.A. and Eymery, J.P. (1965) Steady-State Simulation of an Ammonia Synthesis Converter. Chemical Engineering Science, 20, 281-292. s://doi.org/10.1016/0009-2509(65)85017-5
[6]	Brian, P.L.T., Baddour, R.F. and Eymery, J.P. (1965) Transient Behaviour of an Ammonia Synthesis Reactor. Chemical Engineering Science, 20, 297-310. s://doi.org/10.1016/0009-2509(65)85019-9
[7]	Murase, A., Roberts, H.L. and Converse, A.O. (1970) Optimal Thermal Design of an Autothermal Ammonia Synthesis Reactor. Industrial & Engineering Chemistry Process Design and Development, 9, 503-513. s://doi.org/10.1021/i260036a003
[8]	Nielsen, A., Kjaer, J. and Hansen, B. (1964) Rate Equation and Mechanism of Ammonia Synthesis at Industrial Conditions. Journal of Catalysis, 3, 68-79. s://doi.org/10.1016/0021-9517(64)90094-6
[9]	Singh, C.P.P. and Saraf, D.N. (1979) Simulation of Ammonia Synthesis Reactors. Industrial & Engineering Chemistry Process Design and Development, 18, 364-370. s://doi.org/10.1021/i260071a002
[10]	Elnashaie, S.S., Abashar, M.E. and Al-Ubaid, A.S. (1988) Simulation and Optimization of an Industrial Ammonia Reactor. Industrial and Engineering Chemistry Research, 27, 2015-2022. s://doi.org/10.1021/ie00083a010
[11]	Dashti, A., Khorsand, K., Marvast, M.A. and Kakavand, M. (2006) Modeling and Simulation of Ammonia Synthesis Reactor. Petroleum & Coal, 48, 15-23.
[12]	Azarhoosh, M.J., Farivar, F. and Ale Ebrahim, H. (2014) Simulation and Optimization of a Horizontal Ammonia Synthesis Reactor Using Genetic Algorithm. RSC Advances, 4, 13419-13429. s://doi.org/10.1039/C3RA45410J
[13]	Farivar, F. and Ale Ebrahim, H. (2014) Modeling, Simulation, and Configuration Improvement of Horizontal Ammonia Synthesis Reactor. Chemical Product and Process Modeling, 9, 89-95. s://doi.org/10.1515/cppm-2013-0037
[14]	Akpa, J.G. and Raphael, N.R. (2014) Optimization of an Ammonia Synthesis Converter. World Journal of Engineering and Technology, 2, 305-313.
[15]	Khademi, M.H. and Sabbaghi, R.S. (2017) Comparison between Three Types of Ammonia Synthesis Reactor Configurations in Terms of Cooling Methods. Chemical Engineering Research and Design, 128, 306-317. s://doi.org/10.1016/j.cherd.2017.10.021
[16]	Soave, G. (1972) Equilibrium Constants from a Modified Redlich-Kwong Equation of State. Chemical Engineering Science, 27, 1197-1203. s://doi.org/10.1016/0009-2509(72)80096-4
[17]	Peng, D.Y. and Robinson, D.B. (1976) A New Two-Constant Equation of State. Industrial & Eng. Chemistry Fundamentals, 15, 59-64. s://doi.org/10.1021/i160057a011
[18]	Ahmed, T. (2016) Equations of State and PVT Analysis: Applications for Improved Reservoir Modeling. 2nd Edition, Elsevier, New York.
[19]	Barbosa Neto, A.M. (2015) Desenvolvimento de um Simulador PVT Composicional para Fluidos de Petróleo. Masters Dissertation, University of Campinas, Campinas.
[20]	Tosun, I. (2013) The Thermodynamics of Phase and Reaction Equilibria. Elsevier, New York.
[21]	Gillespie, L.J. and Beattie, J.A. (1930) The Thermodynamic Treatment of Chemical Equilibria in Systems Composed of Real Gases. A Relation of Heat of Reaction Applied to Ammonia Synthesis. The Energy and Entropy Constants for Ammonia. Physical Review, 36, 1008-1013. s://doi.org/10.1103/PhysRev.36.1008
[22]	Ergun, S. (1952) Fluid Flow through Packed Columns. Chemical Engineering Progress, 48, 89-94.
[23]	Edgar, T.F., Himmelblau, D.M. and Lasdon, L.S. (2001) Optimization of Chemical Processes. McGraw Hill Chemical Engineering Series, New York.
[24]	Chapra, S.C. and Canale, R.P. (2009) Métodos Numéricos para Engenharia. McGraw Hill, New York.