The model was derived for a steady-state, packed-bed tubular reactor with variable cross section area and is based on the work of Magnanelli et al. (2019). The main assumptions of the model are: (i) plug flow, (ii) negligible temperature difference between gas bulk and catalyst, and (iii) validity of Ergun’s equation to compute pressure drop.

In a control volume π r 2 d l along the reactor axis the mass balance is: d y i n ̇ d l = ν i   R   π r 2 (1) d n ̇ d l = ∑ k   ∈   i ν k   R   π   r 2 (2) where y i is the mole fraction of species i = { H 2 , N 2 , N H 3 , A r , C H 4 } , n ̇ is the total mole flow rate, r is the reactor radius, and ν i is the stoichiometric coefficient of species i in the reaction: 1 / 2 N 2 + 3 / 2 H 2 ⇌ N H 3 (3)

The reaction rate of ammonia formation, R , is given by the expression (Nadiri et al., 2024): R = K b K c 2 α K a 2 f N 2 f H 2 3 f N H 3 2 α − f N H 3 2 f H 2 3 1 − α (4) where α is the kinetic exponent (0.5 ≤ α ≤ 0.75), K b and K c are kinetic coefficients, K a is the thermodynamic equilibrium constant, and f i the fugacity of species i in the mixture.

The kinetic coefficients, Equations 5 and 6, follow Arrhenius-type kinetic expressions (Nielsen et al., 1964): K b = K b 0   exp − E b R T (5) K c = K c 0 ⁡ exp − E c R T (6) where R is the ideal gas constant and T is the gas temperature. The values for the pre-exponential factors, K b 0 , K c 0 , and for the activation energies, E b and E c , are reported in Table 1.

The thermodynamic equilibrium constant, K a is calculated as a function of temperature following the expression (Gillespie and Beattie, 1930): log ⁢ K a = − 2.691122 ⁡ log ⁡ T − 5.519265 × 1 0 − 5 ⁢ T + 2.6899 + 1.848863 × 1 0 − 7 ⁢ T 2 + 2001.6 T (7) The fugacity of the species i in the gas mixture at temperature, T , pressure, p , and composition, y (mole fraction), is calculated as the product of the fugacity coefficient, ϕ i and the pressure of i in the mixture, f i T , p , y = ϕ i T , p , y   y i   p (8) The fugacity coefficient ϕ i is derived from the Beatty-Bridgeman equation of state (Månson and Andresen, 1986): R T ⁡ ln ϕ i = β i + D i R   T   p (9) where β i and D i terms are calculated according to the expressions: β i = R   T   B i − A i   R   C i T 2 (10) and D i = A i 1 / 2 − ∑ j y j A j 1 / 2 2 + 1 T 2 C i 1 / 2 − ∑ j y j C j 1 / 2 2 − 3 4 R T B i 1 / 3 − ∑ j y j B j 1 / 3 B i 2 / 3 − ∑ j y j B j 2 / 3 (11) The values of parameters A i , B i and C i are provided in Table 2.

The energy balance for a differential section of the reactor d l , neglecting potential and kinetic energy contributions, is the following: d n ̇ h d l = J q   2   π   r   C (12) where n ̇ is the total molar flow, h is the molar enthalpy of the mixture, J q is the local heat flux, and C is the local curvature of the heat transfer surface, C = 1 + U 2 (13) where U represents the local radius variation in the reactor, d r d l = U (14) The enthalpy of the mixture is calculated using the residual approach: h T , p , y = h I G T , y + h R T , p , y (15) The ideal gas contribution to enthalpy, h I G , taking as reference state 298.15 K and 1 atm, is calculated from the standard enthalpy of formation, Δ H f , i 0 and the ideal gas heat capacity, C p i 0 , h I G T , y = ∑ i y i Δ H f , i 0 + ∫ 1 T C p i 0 d T (16) The values for Δ H f , i 0 and C p i 0 are taken from the data base of Poling et al. (2001). The residual enthalpy of the mixture, h R , is calculated from the Beatty-Bridgeman equation of state, h R T , p , y = B m − 2 A m R T − 4 C m T 3 p (17) where A m , B m and C m are calculated from the following the mixing rules: A m = ∑ i y i A i 1 / 2 2 (18) B m = ∑ i ∑ j y i y j B i 1 / 3 + B j 1 / 3 3 (19) C m = ∑ i y i C i 1 / 2 2 (20) where parameters A , B and C are reported in Table 2.

The local pressure drop is estimated from Ergun’s equation, d p d l = − 150   μ D p 2 1 − ϵ 2 ϵ 3 + 1.75   ρ   v D p 1 − ϵ ϵ 3 v (21) where μ is the viscosity of the mixture, D p is the catalytic pellet diameter, ϵ the void fraction, ρ , v are the gas density and velocity, respectively. The values of the parameters used in Equation 21 are given in Table 3.

The local gas flow velocity, v , is obtained using Equation 22: v = m ̇ ρ   π   r 2 (22) where m ̇ is the gas mass flow rate, r is the reactor radius, and ρ is the density of the mixture calculated from Equation 23: ρ = p z   R   T (23) The compressibility factor, z , is obtained from the Beattie-Bridgeman equation of state, Equations 24, 25: z = 1 + β m R   T 2   p (24) and β m = R   T   B m − A m − R   C m T 2 (25) where A m , B m and C m are calculated from the mixing rules presented in Equation 18.

The entropy production in the reactor, ( d S / d t ) irr , is the integral of the local entropy production per unit length, σ ′ , along the reactor length, Equation 26 d S d t irr = ∫ 0 l σ ′   d l (26)

The local entropy production per unit length, Equation 27 is obtained using irreversible thermodynamics as the product of fluxes and thermodynamic forces (Kjelstrup et al., 2010), σ ′ = R   − Δ G r T + v − 1 T d p d l   π   r 2 (27)

The first term on the right-hand-side of Equation 27 is the reaction rate, R , multiplied by minus the reaction Gibbs energy, Δ G r , divided by the local temperature, T . The second term is the product of the local velocity of the gas, v , by the pressure drop per unit length, d p / d l , divided by the temperature.

The thermodynamic reaction force is calculated from Equation 28 (Nummedal et al., 2003): − Δ G r T = − R   ln a N H 3 a H 2 3 / 2 a N 2 1 / 2 K a (28) where R is the ideal gas constant, a i is the activity of species i , calculated as the fugacity of component i in the mixture, f i ( T , p , y ) , divided by the reference pressure, p 0 = 1 atm, while the equilibrium constant, K a , is calculated from Equation 7.

As Base Case (BC) we have taken an ideal reactor of 10 m length with a constant radius of 0.1 m. By ideal we mean that the reactor has an optimal temperature profile that maximizes the mole fraction of the ammonia outlet, y N H 3 ( L = 10   m ) . This is a relevant choice, consistent with the purpose of chemical reactors; to produce chemicals. The maximum production of ammonia is found from the optimal local heat flux, J q , at an inlet temperature of the reactive mixture, T ( l = 0 ) . The optimization problem formulation, Equations 29–31, for this case is as follows: max J q l , T l = 0 y N H 3 L = 10   m (29) s.t.

Equations 1-22 T min ≤ T l = 0 ≤ T max (30) J q , m i n ≤ J q l ≤ J q , m a x (31) where Equations 1–22 are equality constraints corresponding to the differential-algebraic system of the reactor model, including its entropy production. The lower and upper bounds for the inlet temperature and for heat flux were set to T min = 300 K, T max = 950 K, J q , m i n = − 5 MW/ m 2 , and J q , m a x = 5 MW/ m 2 , respectively. The high heat flux constraint allows exploring optimal solutions beyond the current achievable limits in chemical reactor technology. Sensitivity of the BC to the heat flux constraint is addressed in Section 4.3.

The variables of entropy production minimization problem include the radius change, reactor length, heat flux and gas inlet temperature. The model is reformulated in terms of a dimensionless length, 0 ≤ l * ≤ 1 , such that the reactor length, l , is calculated as the product of the dimensionless length and the total reactor length, L tot , Equation 32, l = L tot   l * (32) The differential operator, d / d l , applied to any variable in the model is reformulated as according to Equation 33: d d l = 1 L tot d d l * (33) Once these changes are introduced in the model, the optimization problem is formulated as according to Equations 34–41: min J q l * ,   U l * ,   L tot ,   T l * = 0 ∫ 0 1 σ ′   L tot   d l * (34) s.t.

Scaled Equations 1-22 T min ≤ T l * = 0 ≤ T max (35) J q , m i n ≤ J q l * ≤ J q , m a x (36) U min ≤ U l * ≤ U max (37) L min ≤ L tot ≤ L max (38) r min ≤ r l * ≤ r max (39) m cat l * = 1 ≤ m cat,max (40) y N H 3 l * = 1 ≤ y N H 3 , m i n (41)

The equality constraints apply to the scaled model equations using dimensionless length l * . Equations 36–38 are path constraints corresponding to the optimization controls: heat flux, J q , local radius change, U , and total reactor length, L tot . Equation 39 defines a lower and upper bound for the radius, while the terminal constraints given by Equations 40–41 ensure that the amount of catalyst be less or equal to the catalyst mass spent in the Base Case (BC), and that the ammonia mole fraction be at least equal to the one obtained in the BC. The upper and lower bounds for this optimization problem are shown in Table 4.

The model was implemented in the optimization modeling language Pyomo (Bynum et al., 2021). The differential algebraic model was discretized using orthogonal collocation on finite elements using the pyomo. dae extension (Nicholson et al., 2017). In this case the discretization consisted of 200 finite elements with 3 Radau collocations points. The discretized model was solved in the large scale nonlinear optimization solver IPOPT 3.13.2.

The Base Case (BC) is a reactor with maximum ammonia produced, given a heat flux and gas inlet temperature. The optimal gas temperature profile shown in Figure 1a implies a high inlet gas temperature (870 K) followed by cooling along the reactor. The main determinant of this behavior is the kinetic expression, Equation 4, which takes into account two antagonistic contributions: on one hand, the term K b / K c 2 α increases with temperature, on the other hand, the equilibrium constant of the reaction, K a , decreases as the temperature increases since the reaction is exothermic.

The reactor temperature profile (Figure 1a) is high at start (above 800 K) and declines sharply below 800K (Figure 1c). In the high temperature zone, the system is far from chemical equilibrium. The reaction rate experiences a very sharp decrease (Figure 1c). The high temperature zone corresponds to the first 0.5 m of the reactor (5% of the total reactor length) and is the most efficient reaction zone generating 40% of the total ammonia produced in the reactor.

An important point not handled here is catalyst deactivation. In terms of catalyst utilization, the high-temperature zone (up to 0.5 m along the reactor length) produces 0.62 mol/s of ammonia per kilogram of catalyst, compared to 0.048 mol/s per kilogram in the low-temperature zone (beyond 0.5 m along the reactor). However, catalysts exposed to temperatures above 800 K can rapidly deactivate, while below this threshold, catalyst deactivation remains low (Liu, 2014; Chonggen et al., 2011). Under these conditions, 95% of the total catalyst mass contributes to 60% of the total ammonia production.

The optimized reactor produces 22.3 mol of ammonia per second (32.7 metric ton/day) with an effluent composed of 22.01 (mole %) ammonia (Figure 1b). This optimal ammonia mole fraction is close to the 23.45 (mole %) obtained by Månson and Andresen (1986) who optimized an ammonia synthesis reactor with similar inlet composition and pressure, but without taking into account the pressure drop. The main difference between our and their results lies in the higher gas inlet temperature of about 950 K, which leads to a higher ammonia mole fraction profile. Nonetheless, the optimal temperature and ammonia profiles reported by Månson and Andresen (1986) are in good agreement with those obtained in this work (see Figures 1a,b). This makes it now interesting to move one more step, and see if the energy efficiency also can be improved in terms of having a smaller entropy production.

The total cooling load required to carry out the process is 2.14 MW. As much as 67.3% of it (1.44 MW) is transferred in the first 2 m of the reactor (Figure 1d), while the last 2 m represents only 4.6% (0.1 MW) of the total transfer. One reason for this behavior, in addition to the reaction rate and thermodynamic equilibrium considered in the beginning of this section, is that the control of the heat released by the reaction is mainly located in the first part of the reactor.

The total pressure drop for the Base Case was 0.432 MPa (Figure 2a), corresponding to 1.65% of the inlet pressure (26.462 MPa), this pressure drop is constant per reactor length (0.0432 MPa/m) due to slight variations in flow velocity, between 1.8 and 2.4 m/s, (Figure 2b). In practical terms a lower pressure drop has a positive effect on reactor performance since ammonia production is thermodynamically and kinetically favored by high pressures, additionally higher reactor outlet pressures reduce the cooling duty necessary to recover ammonia by condensation before recirculation of unconverted reactants which would reduce econmic cost.

Entropy production minimization was carried out keeping the same ammonia production (22.3 mol/s) and composition (22.01% molar) as well as total mass of catalyst (290 kg), using local radius change U , reactor length, L tot , gas inlet temperature, T ( l * = 0 ) , and local heat flux, J q , as controls.

The optimal geometry (Figure 2c) results in a reactor with a length of 2.9 m and a radius distribution that increases linearly from 0.14 m in the inlet to 0.22 m in the outlet. The geometry diminishes the pressure drop by 96% with respect to the Base Case. The main reason for this lies in the in gas velocity reduction which has a quadratic contribution to the pressure drop expression as given by Ergun’s equation, Equation 21.

The gas inlet temperature in the reactor with minimum entropy production reactor is, however, 80 K higher than in the Base Case reactor (Figure 2d). Higher inlet temperatures increase the initial reaction rate of the minimum entropy production reactor followed by a steeper decrease than in the Base Case reactor (Figure 2e). Temperature profiles show that the length of hot reaction zone, T ≥ 800 K, is the same in both cases (0.5 m); however, the temperature in the hot zone tends to be higher in the minimum entropy production reactor. Besides, the temperature in the cold zone of the reactor with minimum entropy production is lower than this zone in the Base Case reactor.

The heat load distribution (Figure 2f) in the Base Case and the in the minimum entropy production design follow similar patterns. In both cases the greatest heat load is located in the first part of the reactor. The total cooling load in the reactor with minimum entropy production (2.49 MW) is 16% greater than in the Base Case. The main differences in heat load is from the first 0.5 m to the end of the reactor length. A way of using the cooling heat to increase thermal efficiency in the process is by generating steam (Flórez-Orrego and de Oliveira Junior, 2017b), which could increase the thermal efficiency of the minimum entropy production reactor with respect to the BC reactor.

The entropy production in the Base Case is 0.139 kW/K where 0.102 kW/K, 73% of the total entropy production (Figure 3a), is due to the reaction term in the entropy production equation. The reactor with minimum entropy production, Figure 3b, reduces the total entropy production by 57% with respect to the Base Case. A salient aspect of the thermodynamically optimal reactor is that the contribution of the reaction and pressure terms in the entropy production is 57% and 99% lower with respect to the Base Case, leading to a reactor design where almost all the irreversibility is due to reaction, since this term represents 98% the total entropy production. These are significant changes.

It is interesting to note that the local entropy production profiles show that entropy production is more evenly distributed in the EM case (Figure 3d) than in the BC (Figure 3c). The result give some support to the principle of equipartition of entropy production which identifies systems with minimum entropy as those that distribute entropy production evenly throughout the system (Tondeur and Kvaalen, 1987). However, exact equipartition is not expected for constrained optimization.

The heat transfer and geometry controls used in the EM case concerned two key reaction aspects, viz, the reaction temperature and pressure. As previously discussed, the first part of the EM reactor, where most ammonia is produced, has a low pressure drop with respect to the BC. The higher pressure present in the EM thermodynamically favors the production of ammonia. Geometry optimization is key to reduce the compression work required for ammonia condensation once it leaves the reactor, and to increase the reaction performance.

In Section 4.1, it was assumed that the heat flux J q could be freely adjusted to achieve optimal operating conditions. However, in practical applications, heat transfer is often limited by material and design constraints. Such limitations may significantly influence the reactor’s thermodynamic performance by increasing entropy production. We propose that a reduced maximum heat flux hinders operation close to reversible conditions, thereby amplifying irreversibilities. Consequently, we expect that this will affect critical operational and design parameters such as the temperature profile, catalyst use, pressure drop, and reactor radius.

The heat flux bounds for the EM case were − 5 ≤ J q ≤ 5 MW/ m 2 . Now the effect of reducing the available heat flux on the entropy production is explored. To this end we define a variable Equation 42 κ q such that: κ q = J q , m a x 5  MW / m 2 (42)

From this expression we define variable lower and upper bounds of the J q constraint in terms of κ q , such that: − 5 ⋅ κ q ≤ J q ≤ 5 ⋅ κ q MW/ m 2 .

The optimization problem given by Equations 36–38 is solved for each value of κ q (0.1 ≤ κ q ≤ 1) in the constraint related to the allowed heat flux, Equation 36.

Figure 4a presents a Pareto front illustrating the relationship between the maximum allowable heat flux and the entropy production rate. It is observed that as the heat flux J q decreases, the entropy production rate in the reactor increases. This suggests that a larger heat flux enables the reactor to operate in a more reversible manner, thereby reducing thermodynamic losses. When heat transfer is limited, the reactor experiences greater irreversibilities associated with chemical reaction and pressure drop. By examining Figure 4a, a point can be identified between the limits of 2 and 3 MW/ m 2 , beyond which, even with a higher heat flux allowed (e.g., up to 5 MW/ m 2 ), the additional reduction in generated entropy is comparatively small. Moreover, increasing the system’s maximum heat flux beyond this point may not be economically justifiable. From a practical perspective, it may be more feasible and cost-effective to design a cooling system that operates near the point of diminishing returns, around 0.5 MW/ m 2 (Figure 4b). It is not necessary to oversize the cooling system beyond the point where the benefits in terms of entropy production reduction are marginal.

Figure 4b shows the heat flux profiles along the reactor for different values of J q , m a x . For higher J q , m a x values, the heat flux profile exhibits more pronounced variations along the reactor, indicating the need for more intense heat transfer in the first part of the reactor to maintain the optimal temperature profile. As J q , m a x decreases, the heat flux profile becomes more uniform, and the absolute values of the heat flux are lower. Horizontal segments in some regions mean that the heat flux has reached its imposed upper or lower bound. Although the local heat fluxes are lower when heat transfer limitations are increased, to maintain the same conversion, the reactor length must be longer, as seen in Figure 4b. The total heat load removed will then be lower, as shown in Figure 4c.

As shown in Figure 4d, the gas temperature at the reactor inlet tends to decrease as the maximum allowed heat flux decreases. This occurs because, the optimization procedure aims to prevent an excessive increase in temperature in the initial section of the reactor, where the exothermic ammonia synthesis reaction is faster and releases a significant amount of heat with a lower heat removal capacity. By starting with a lower temperature, the heat generated along the reactor can be better managed under limited cooling capacity. This is consistent with observations showing that, toward the end of the reactor, temperatures tend to rise as the maximum heat flux decreases. The presence of “hot spots” in the temperature profile becomes more evident as the maximum allowed heat flux is reduced. This is because the limited capacity to remove the heat generated by the reaction leads to local accumulations of energy. When the reactor’s cooling capacity is limited, it becomes impossible to maintain the optimal temperatures associated with higher heat flux, which prevents the minimization of entropy production. From a local perspective, to initiate the process at a lower temperature may reduce the entropy produced by the reaction; however, this also significantly lowers the reaction rate at the reactor inlet (Saffari et al., 2021). To achieve the same ammonia conversion under these conditions, the reactor would need to operate at higher temperatures in other zones and may additionally require an increased length to compensate for the reduced thermal efficiency, resulting in greater overall entropy production.

To prevent catalyst deactivation, it is necessary to limit the temperature to 800 K, (Chonggen et al., 2011; Liu, 2014). It is observed that under low heat transfer constraints ( J q , m a x between 5.0 and 1.0 MW/m²), hot spots with temperatures exceeding 800 K appear in the initial sections of the reactor, covering between 15% and 25% of the reactor length. These conditions may limit reactor operation, at least until new developments in thermally more stable catalysts are achieved. However, the thermal profile obtained under the highest heat transfer restriction that is ( J q , m a x = 0.5 MW/ m 2 ) exhibits satisfactory temperature levels that allow the use of existing catalysts and also enable the implementation of standard heat exchange systems. This profile is similar to those observed under industrial conditions in reactors with integrated heat exchange, represented by circle black markers in Figure 4d (Baddour et al., 1965). Under these heat transfer conditions, the entropy production is 0.17 kW/K, which is close to the value found for the Base Case (0.14 kW/K, Figure 3a). However, for the Base Case, the required heat transfer rate is 5.0 MW/ m 2 , which prevents its implementation using conventional heat transfer systems.

Figure 4e illustrates the relationship between pressure drop across the reactor and the maximum allowable heat flux for different design constraints. It is evident that as the maximum heat flux increases, the pressure drop across the reactor tends to decrease, leading to lower energy losses due to pressure drop and thereby minimizing this source of irreversibility. This reduction in pressure drop is achieved at the expense of an increase in reactor diameter, as shown in Figure 4f. The trend observed in Figure 4f suggests that when higher maximum heat flux per unit area is allowed (i.e., higher J q , m a x ), the optimal reactor design favors a geometry with a larger radius along its length. Conversely, when heat transfer is more restricted, the reactor radius decreases. Under these conditions, heat transfer is maximized through an increase in the lateral surface area-to-volume ratio (Mahapatra and Dasappa, 2014). This phenomenon also explains the marked reduction in the reactor diameter in the initial sections, followed by an increase in diameter toward the middle and end of the reactor. Since hot spots occur in these regions, a reduction in diameter enhances the heat transfer area per unit volume, thereby compensating for the limitations imposed on heat transfer. For the simulated capacity of 32.7 metric ton/day, this configuration features (Figure 4f line red) a variable-radius reactor geometry and a temperature profile below 800 K, ensuring catalyst stability.

As observed, a tubular reactor with constant radius (Base Case) results in high entropy production (0.14 kW/K), whereas a reactor with a monotonically increasing variable radius, achieving the same ammonia production and operating under the same allowable heat flux (5 MW/ m 2 ), produces approximately 0.06 kW/K (only 43% of the entropy production with respect to the Base Case). However, this heat flux is excessive for current heat transfer systems. Therefore, a more realistic design, corresponding to a maximum heat flux of 0.5 MW/ m 2 , leads to a geometry like the one shown in Figure 4f (red line) and Figure 5, with an entropy production of 0.17 kW/K, similar to that of a tubular reactor operating at 5 MW/ m 2 maximum heat flux.

From a mechanical design perspective, variations in radius along the reactor axis, especially if pronounced, would lead to a non-uniform distribution of stresses along the reactor walls due to internal pressure. To withstand these variable stresses and the high operating pressures, high-strength materials would be necessary, along with more sophisticated manufacturing processes, such as additive manufacturing (3D printing), instead of using a simple constant-radius tube. This would require further detailed studies, both at the simulation level and through experimental prototypes.

It is well known in chemical engineering that more reversible conditions can be approached by increasing the contact area for heat exchange. Similar systematic studies on decrease in the entropy production of the chemical reactor with varying geometric variables, has not been done. That geometric variables are important is, however, beyond doubt. This opens up for new important applications of our well-proven engineering mathematical procedures.