Modeling and Numerical Simulation of Ammonia Synthesis Reactors Using Compositional Approach

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).


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) 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 highpressure deviates gas from ideal treatment, approximating gas molecules.
As most of substances are non-polar (except for NH 3 ) 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.

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).

( )( )
For the results of this work, only PR-EoS was used, with the parameters ( ) 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: 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 A ij , in the Equation (7), is calculated from Equation ( 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.
( ) ( ) 11) Advances in Chemical Engineering and Science The expression shown at Equation (12) for fugacity coefficients is obtained.
( ) ( ) ( ) With 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 [20].

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].
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.
To summarize this section, the equilibrium constant (K eq ) in Equation ( 19) is found in correlation given by Gillespie and Beattie [21].

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.
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].
( ) ( ) 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].

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.

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

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.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.

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.

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.   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 Advances in Chemical Engineering and Science

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.

Sensitivity Analysis
In the sensitivity analysis the same data set of validation was used to fit α. 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.

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.15K) 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 2 can be converted, giving a more accentuated profile (T in = 643.15K).
On the other hand, high inlet temperatures present high reaction rates, converting more N 2 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.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 2 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.15K and 663.15K is smaller than 663.15K 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.

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 2 •K is 753 K, while for U = 650 W/m 2 •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 2 •K, the final T value is 724 K and U = 650 W/m 2 •K is 687 K.It takes place because a large U removes more energy in the reactor, decreasing temperature more rapidly.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 2 •K, there are errors computing η (higher than 1).

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.15K 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 2 •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.Step size along the reactor in numerical method Gas viscosity

List of Greek Letters
. Defining the following Equations (3) to (5) below:

Figure 4
Figure 4 presents the N 2 conversion variation in an adiabatic converter with α alterations.

Figure 4 .
Figure 4. Conversion profile in first bed using α variations at MARS (EoS approach).
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.

Figure 5
Figure5shows 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 .
Figure 5. Plant data and MARS EoS model temperature profile (adiabatic reactor).

Figure 6 .
Figure 6.Plant data and MARS EoS model conversion profile (adiabatic reactor).

Figure 7 .
Figure 7.Comparison between the MARS EoS model and plant data (autothermal reactor).

Figure 14 presents
Figure14presents 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.

1 δ 2 δΩΩψ
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) a Constant for cubic EoS (PR or SRK) b Constant for cubic EoS (PR or SRK) i Mixture factor in cubic EoS for i component
Another challenge is the interdependence of variables in differential equations.Therefore, the code was made using a modular structure in Wolfram Ma-thematica® Programming Language.Moreover, our algorithm is called MARS (Models for Ammonia Reactor Simulation).
4.75 Advances in Chemical Engineering and Science

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.

Table 7
presents the relative errors at each point of the autothermal reactor.It D. S. S. Jorqueira et al.DOI: 10.4236/aces.2018.83009134 Advances in Chemical Engineering and Science