Ind. Eng. Chem. Process Des. Dev. 1982, 21 359-367
359
Modeling and Simulation of an Ammonia Synthesis Loop K. V. Reddy and Asghar Husaln' Regional Research Laboratoty, Hyderabad 500 009, India
The paper describes modeling and simulation of the synthesis loop of an ammonia plant under operation. Models for the different process units are formulated using basic principles, and the unknown parameters are determined by matching their outputs with the plant data. The interconnected models are then used to simulate the loop performance under various conditions. Effects of H2/N, ratio in the recycle gas, loop pressure, recycle gas flow rate, and inerts concentration are studied on the ammonia production rate, fractional hydrogen conversion, and gross profitability. Among these, the H2/N, ratio is found to be most crucial and its optimum value lies around 2.50, which is further confirmed by plant implementation of the simulation results. At this value, one can easily
locate optimum of other parameters for a known plant configuration.
Introduction Ammonia is among the important chemicals produced the world over in high tonnage due to its wide use in the manufacture of fertilizers. A typical ammonia production process consists of (a) production of the synthesis gas, (b) compression of the gas to the required pressure, and (c) synthesis loop in which its conversion to ammonia takes place. Although the first two sections have their own importance, the converter which is part of the synthesis loop is crucial in the overall control strategy of the plant. Hence, the present work concentrates on the modeling and simulation of the synthesis loop of an ammonia plant in operation, which uses the petroleum feedstock as a source of hydrogen. Figure 1 shows a simplified diagram of the synthesis loop under study, which is operated under a pressure of 300 X 105-500 X lo5 Pa. The product gas from the converter at about 453 K passes through the waste heat boiler where it is cooled to about 413 K. This gas flows through the condenser and is further cooled to below its dew-point temperature. The gas-liquid mixture from the condenser then enters into the ejector where it is mixed with the fresh compressed gas, and the liquid ammonia is separated from the gas. This gas is recycled to the converter. The fresh feed is at about 308 K and consists of hydrogen and nitrogen in approximately stoichiometric proportions together with small amounts of argon and methane. The liquid ammonia separated in the ejector drum flows to a letdown column and drum where pressure is successively reduced to remove the dissolved gases (i.e., H2, N2,Ar, and CHI). The product from the letdown drum is 99+ 7%pure liquid ammonia. Various purge streams shown in Figure 1 are necessary to maintain inerts (Ar and CHI) in the system to a minimum level, of course, at the cost of letting off part of the hydrogen. In the present study (Reddy, 1980), mathematical models of the individual units in the synthesis loop (Figure 1) are developed based on the physicochemicalprinciples involved pertaining to the phenomena occurring in each of them. The unknown model parameters are generated by matching the model outputs with the consistent plant data collected over a period of time. Care has been exercised to generate such parameters which are physically consistent in all respects. The developed models for various units are then interconnected and used to simulate the loop performance under various operating conditions and optimizing it with respect to the ammonia production rate. 0196-4305/82/ 112 1-0359$01.25/0
It will be clear that the approach adopted in this study is based on writing specific computer programs for the individual loop units, rather than using any of the standard modular packages such as reviewed by Flower and Whitehead (1973) or Motard et al. (1975). That way it has been possible to incorporate much greater depth of detail with an aim to provide better insight into the process. Plant Data The plant data are collected in two different periods over a span of nearly two years, for consecutive days in each period and at regular intervals in the three shifts of operation per day. The loop being a high pressure system, no flow meters are provided in the plant except the one which measures the ammonia production rate (A3,Figure 1). In addition, the following measurements are available: (i) temperatures of streams F, F,, P, P,, C1, C2,W1, W2, A,, and the converter bed temperatures; (ii) pressures of streams F,, F, letdown column and drum; (iii) analyses of streams F, F1, P, and P1. The day averages are found and data sets 1-5 showing reasonable consistency with the flowsheet capacity, Figure 2, are used in the modeling and simulation study. Composition and temperature of stream A, are estimated by dew-point calculation. Then compositions and temperatures of the outlet streams of the letdown column and drum and their liquid/vapor ratios are successively calculated by flash vaporization. In subloop I (Figure l),nitrogen and methane balances provide flowrates of streams F and P1.These particular components are selected in view of the facts that (i) total number of moles reduces due to chemical reaction in the converter, (ii) component balances for nitrogen and hydrogen cannot be independent, and (iii) in the plant analysis only nitrogen, methane, and ammonia contents are determined. In the subloop I1 (Figure l), since no chemical reaction occurs, therefore, total balance and any single component balance provide flowrates of streams F1 and P. Converter In the converter the following catalytic reaction takes place Nz + 3H2 ;=t 2NH3 Since the forward reaction occurs at elevated temperatures releasing a large amount of heat, this heat has to be removed to obtain a reasonable overall conversion as well as to protect the life of the catalyst. At the same time, the 0 1982 American Chemical Society
360
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 3, 1982
I
Purge Gas,P2
Purge Gas, P)
i I
I 1
I
I
LETDOWN
~
I
I
i Liquid NH3 (Product)
Figure 1. Ammonia synthesis loop. ,Outer Annulus
j
/Innor Anrulus
Catalyst Basket
//
/Cooling
Ibes
410 LOG
390 380
i
370
,360
z
Z 350 340
330
- 320 LL
310
$ 2i
290
300
280
>
*
270 260
A m m o n 1 0 Pioduclion Rate, lonnes/ h w r
Figure 2. Plant material balances.
released energy can be utilized to heat the incoming feed gas to the proper reaction temperature. For this reason, an ammonia converter is self-supporting. As shown by van Heerden (19531, an autothermic process possesses multiple steady-state solutions; hence, its stable operation is quite crucial. The ammonia converters are available in various designs, as described by Walas (1959). Their modeling has received considerable attention during the past three decades. Annable (1952) is among the first few who studied modeling of the steady-state reactor by using the well known Temkin-Pyzhev (1939) rate equation, while placing major emphasis on the kinetic parameters. Kjaer (1958) formulated model equations for a Haeber-Bosh type of reactor by accounting for variations in the temperatures in both the longitudinal and radial directions. Zayarnyi (1962) also developed a steady-state model for the Haeber-Bosh type of converter. Baddour et al. (1965) reported a simple model for the TVA reactor, representing the plant temperatures and ammonia production rate within 1 5 2 0 % .
QS
Product
Gas
I Cold-Shot
Figure 3. Casale synthesis converter.
Shah et al. (1965) also reported modeling of the Haeber-Bosh type converter, their primary aim being to observe the effect of different variables on the operational stability of the converter. Later, Shah (1967) carried out a simulation study on a quench type converter with the same objective. Almasy et al. (1967) concentrated on the design aspects of the heat exchange portion of the NEC,TVA, and simple cocurrent tubular converters. Paavo Uronen (1971) reported modeling of the quench type converter, suggesting optimal purging control of the loop. Kubec et al. (1974) developed a model for a radial flow quench type converter. Gains (1977) showed substantial improvement in efficiency by maintaining optimal temperatures in a quench type
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 3, 1982 361
converter. Singh and Saraf (1979) considered diffusion effects within catalyst pores to describe the synthesis rate equation for catalysts of different makes; they found generally good agreement between plant data and simulation results for both the autothermic reactors and those having adiabatic catalyst beds with interstage cooling. In the present study a Casale type of converter is involved, which is more complex in construction as shown in Figure 3. It consists of two annular spaces (channels 1 and 4) surrounding the catalyst basket. In the catalyst bed (channel 3), 84 cooling tubes (channel 2) of 35 mm o.d., 30 mm i.d., are arranged in a triangular pitch of 85 mm. The promoted iron oxide catalyst particles ranging in size from 12 to 21 mm surround the cooling tubes. Some of the tubes contain electrical heaters for use during the start-up; the rest of the tubes have similar tube cores to get a uniform flow of gas in all the tubes. Inside and outside diameters of the catalyst basket are 928 and 948 mm, respectively. On the outside of the basket, rods of 8 mm diameter are welded to form helices, which act as turbulence promoters in the channel 4. Similar helices are present in the channel 1. In channels 1 and 2 the synthesis gas is heated during its flow. The heated gas mixed with the cold-shot (bypass containing portion of the cold synthesis gas) enters the catalyst bed (channel 3). After the reaction, the product gas flows through channel 4, transferring heat to the incoming synthesis gas. In writing the model equations, radial temperature gradients and corresponding concentration gradients are neglected in view of the results reported by Kjaer (1958). Any axial diffusion of enthalpy is similarly ignored in the light of results given by Emery (1964). Since the gas mixture flows at a sufficiently fast rate, its residence time being a few seconds, not much additional information of practical significance would have been available by making the model more complex in the form of partial differential equations; hence, time has not been included as a variable. The steady-state model so developed consists of the following ordinary differential equations which describe mass and heat transport within the four channels of the converter. dTi/dZ = [ U ~ I A ~ I-( T Ti) ~ - (& + ha)As(Tsk - Ts)l/SCFiCpi (1) i
dT2/dZ = - U343dT3 - T2)/SCFiCpi
(2)
I
dT3/dZ = [2/3(-Wt)lr3fA33- U32AdT3 - T2)U34A34(7'3 - T4)I/Cnicpi (3) i
dT4/dZ = - [U34A34(T3 - T4) - U41A41(T4 - Tl)I/Cn&pi (4) i
d X 3 / U = fr3A33/F3
(5)
dTD/dZ = (T3 - T3m)2 Equations 1 to 4 describe the change of temperature along the length of the converter in each of the channels 1-4. Equation 5 gives the rate of fractional conversion of hydrogen in the catalyst bed. Equation 6 is introduced for bringing a close match between the observed and computed bed temperatures. The details regarding solution of the above equations, known parameters, and generation of unknown parameters are given elsewhere (Reddy and Husain, 1978). In addition, abridged information is presented in the Appendix. Greater emphasis has been given to replicate within reasonable accuracy measurements from the plant before the converter model is in-
900
850
650
600
1
1
- Measured _____
Computed
>
550 0 500
02
01
Converler
03
04
Fraclmal
05
06
Length
07
OS
09
10
(Bottom lo Top)
Figure 4. Typical match of bed temperatures. Table I. Equilibrium Concentrations at 323 K
pressure, Pa(x l o 5 )
98.06 294.18 490.30
mol % ammonia in the vapor phase vendor supplied data computed
26.00 13.20 10.60
26.15 13.51 10.62
corporated in a simulation study of the loop. Figure 4 shows a typical match obtained of the catalyst bed temperatures.
VLE Vapor-liquid equilibrium relationships for ammonia in the presence of impurities (i-e.,N2, H2, Ar, and CHI) are needed in the modeling of other loop units namely condenser, ejedor, letdown column, and drum. Highly reliable VLE data are available to the authors from vendors of the plant in the range 273-323 K and 49 X 105-490 X lo5 Pa, with a vapor phase molar ratio of H2:N2:Ar:CH4equal to 3:1:0.19:0.46. These are suitably correlated thermodynamically by evaluating vapor phase fugacity coefficients using the Redlich-Kwong equation of state, with one of its parameters being temperature dependent for ammonia; mixing rules using parametric data of the pure individual components are applied. An algorithm for computing the equilibria is given elsewhere (Reddy and Husain, 1980). Three typical concentrations of the vapor phase in equilibrium with the liquid phase are given in Table I. Condenser The ammonia condenser in the synthesis loop is a vertical shell and tube type without any baffles. Condensation taking place inside the tubes at high pressure is a complex phenomenon because of the simultaneous heat and mass transfer as well as the solubility of other components in the condensed ammonia. In principle, for a rigorous modeling one has to account for the mass transfer fluxes of all the components present. However, inconsistency can arise in calculating mass transfer coefficients for the multicomponent system under consideration, using the relationships such as those given by Bird et al. (1960) or Treybal (1968), due to at least (n + 1) equations being available for determining n fluxes. At the same time, the computed coefficients will not account for the interaction of fluxes.
362 Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 3, 1982
As an alternative, Stefan-Maxwelldifferential equations (Hirschfelder et al., 1964) can be solved. But for a fivecomponents system their solution will be extremely laborious. Since the total solubility of noncondensables in ammonia is not more than 5%, a simplified approach is adopted here; it is based on considering condensation of ammonia alone and accounting for the solubility of noncondensables at the exit conditions only. Hence, considering the gas entering the condenser as a binary system of ammonia and noncondensables (i.e., N2, Ha, Ar, CHI) the following model equations are written for a differential length of the consenser, which are applicable as long as the interface temperature remains above the saturation temperature of the gas
/ I
/
- Equipment Parameters - Operating Data - c in Eq (10)
-~
_ - ~ -
G
a
t
L e
@and
1
(7)
-dT,- - -h,(Td - TJTDINt dZ
I
(8)
WCPW
The interface temperature Td is then calculated by ha(T,j - T,)= ht(To - Td)
J
(9)
Condensation starts as soon as the interface temperature falls below the saturation temperature of the vapor. In that case, the following equations are applicable
b Results
Figure 5. Procedure for condenser model parameter estimation.
producing a match with the plant data; the correlation is as follows
In eq lo,+ is a correction factor, known as Ackerman's correction (Schrodt, 1973) to take care of the sensible heat effects between the bulk gas phase and the interface. In order to determine Td,which will be below the saturation temperature in this case, a heat flux balance is made as f 0110ws h,(Td - T,) = ht(To - Td)$
+ NlMlh + D,
(13)
For the section in which condensate film is already existing
-+ha
> 2100, where Gw = w/s
De =
4(P:
- rDO2/4)
XDO Correlations for the coefficients ht and hl are available from Kern (1950) and McAdams (1954), respectively. Mass Flux. In order to calculate mass flux, N,, the following equation is used
hl
Fc,kflc,,(Tdi-' D, =
for (D,G,/p,)
- Tdi)
TD1NtA.Z
(15)
However, for the section in which condensationjust starts, D, = 0 and h, = h,. Heat Transfer Coefficients. In eq 14 the coolant side heat transfer coefficient, h,, is to be evaluated using the heat transfer coefficient between coolant and the interface, h;, and fouling and tube wall resistances. For a baffled condenser, a correlation is given by Kern (1950) for evaluating h i . But for a condenser without baffles, which is the present case, no separate correlation is available. Therefore, the same type of correlation can be used by introducing a correction factor C to be determined by
Sh = O.O23(
-) '.* ( DIGf Pm
?)I3
(22)
Thus, in the above model equations for the condenser the only unknown parameter is C in eq 16 which is determined by matching the model output with the plant data. The plant measurements of the bulk gas temperature, To, at the condenser exit are available. Therefore starting with an initial value of C, several trials are made to correct its value by producing a match between the
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 3, 1982
Calculate -Component Fugacity in Vapor Phase -Liquid N$ Fugacity ai Saturation Temperature 1
Calculate
'
7
1
I
- Activity Coetficients - Pressure Corrections - Henry's Constants
*
-
I 1
~
363
density, vapor pressures, and latent heat data from the International Critical Tables (1926) are fitted with suitable polynomials. For calculating viscosity and thermal conductivity of the gaseous mixture, the procedure given by Kjaer (1963) is applied; the critical constants are those from Perry (1973). The diffusion coefficient of ammonia in the noncondensables is obtained from Gilliland's (1934) correlation. This is multiplied by pm, being the density of the mixture as estimated from the RK equation of state, to give bulk phase diffusion coefficient, &, as required in eq 20 and 22. Ejector A single-stage ejector attached to a vertical drum mixes product from the condenser with the fresh synthesis gas, the latter serving as the primary fluid. The liquid portion of the condenser product is thus separated from the gas by falling into the vertical drum. A rigorous ejector model should consider mass, enthalpy, and momentum balances. However, preliminary calculations indicated that pressure drop over the loop is low enough to consider any momentum balances. Hence, the model consists of the following overall material balance F, F1 F = F, + F, (23)
+ +
component material balance
Results
FvYvi + F l ~ l + i FYfi = Fryi liquid/vapor ratio
+ FeXi
(24)
R = F,/F,
(25)
energy balance data set
measured
computed
1
320.0 323.0 316.0 321.0 318.0
322.9 324.7 315.0 318.6 314.8
2 3 4 5
computed and the measured Tovalues. For this purpose, model equations are integrated from top to bottom of the condenser using Euler's method according to the procedure shown in Figure 5. The best match obtained for the data sets 1-5, as given in Table 11, evaluates the proper value of C in eq 16.
Physical Properties Several properties are required for solving model equations of the converter as well as of the condenser. The relations for specific heats of the individual components and the heat of reaction based on ammonia are given by Shah et al. (1965). The specific heat of the gaseous mixture is evaluated by summing up molal contributions of the components. Molar volumes for nitrogen, hydrogen, and methane are available from Perry (1973), and for argon an equation from Reid and Sherwood (1966) is used. FOPliquid ammonia thermal conducitivity, viscosity, specific heat,
Fvepm(To - TJ + FlCpl(T0 - Te) + (F,Yvi + FY, - FrYi) = FCpf(T, - Tf) (26) The iterative procedure for solving eq 23 to 26 is shown in Figure 6, in which the vapor-liquid equilibrium relationship given by Reddy and Husain (1980) is applied. The average values of measured recycle gas compositions for various data sets are compared in Table I11 with those computed at the average values of pressure and ejector inlet conditions; it shows satisfactory agreement within the measurement accuracy. Letdown Column and Drum These involve the normal flash calculations, in which the VLE relations as given by Reddy and Husain (1980) are used. Loop Simulation Having all the models for the individual loop units formulated and their validity checked, they are now interconnected to simulate overall performance of the synthesis loop. Four independent parameters, namely (i) H2/N2ratio in the recycle gas, (ii) loop pressure, (iii) recycle flowrate, and (iv) inerts concentration in the recycle gas, are identified to study the effect of each of them, one at a time, on the following: (a) ammonia production rate,
Table 111. ComDarison of Recvcle Gas ComDositions recycle gas composition, mole fraction data measured computed set no. NH, N* H2 Ar CH, NH 1 N2 H* 0.1796 0.6306 1 0.0997 0.0257 0.0644 0.0990 0.1700 0.6409 2 0.0898 0.1886 0.6270 0.0672 0.1024 0.1756 0.0274 0.6280 3 0.0792 0.1932 0.6467 0.0273 0.0536 0.0924 0.1866 0.6413 4 0.0873 0.1866 0.6467 0.0273 0.0521 0.0956 0.1780 0.6485 5 0.0885 0.2023 0.6213 0.0285 0.0594 0.0889 0.1939 0.6239
~~
Ar 0.0258 0.0273 0.0270 0.0269 0.0286
CH, 0.0641 0.0665 0.0525 0.0508 0.0591
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 3, 1982
364
,!
-Recycle Conditions -Equipment Parameters -Synthesis Gas Flowrates -Coolant Flowrate --Initial Temperatures for Converter
I
Bed Inlet and Product Gas
0 31
-
115
-
115
030-
;
1 046
-1 -
d
Heat
3i
,
Condenser Model
042
2
- 1 040
-1
038
+
- I 036 1 034
r
I
i flonvergence Obtained o Recyle Conditions
046
1 044
-1
j 1 1
i
,
,&Fo(Mode!
\
l
E d i f y Initial I Temperatures
Boiler
I 058
- 1 056 - 1 054 - 1 052 - 1 050
I
onvergence Obtained on Converter Channels 1 8 2 . 3 Top Temperatures
060
-1
1
I
n
h
11 024 I
1 022
Column & Drum
Figure 7. Strategy for loop simulation. 0 35-
0 34114
K
1 .:7
115 r
I
-
3 33-
i 05
0 32-
1 34
031-
1
190
I
ILCE
03
u LCZ
F030-
-*
I*
0'
ij 0 2 9 -
1
i
028-
1 00
0 27-
G 99
0 26-
0 98
0 97
0 25--20
0
25
3
3 0 n 2 I N2 R o t c
3 5
3
Loop P r ~ S L u l t- 336 34 i IO5 pa, R l i y i l a R o l e - 90Zgm malrs/5sc ,lner*% :8 2'1.
Figure 8. Effect of Hz/N2ratio.
(b) fractional H2 conversion, and (c) profitability ratio. In each case, the gross profit is first calculated from gross profit/day = (sale revenue/day) (compression cost/day) - (reforming cost/day) (27) The reforming cost per ton of gas is provided by the plant management. The compreasion cost is obtained from a correlation of the plant data in tenas of the final pressure of Compression and flowrate of the gas. Then, taking gross profit per day of the conditions of data set 3 as one unit, the profitability ratio is calculated for all other cases. The strategy used for the overall loop simulation is shown in Figure 7, and the results are plotted in Figures 8-11. It is clear from Figure 8 that the Hz/N2ratio equal to 2.5 in the recycle gas is the best one, providing maximum ammonia production rate as well as profitability. The same result is obtained in the converter simulation
700
0
800
900
0 (Loop Pressure
IO00
Ricycl. i
1100
12W
Flowrotr
, qm
336 3 4 I l o s p a ,
n2/Nz
i
13W
14W
1500
0
moles / s c . 2 5,
lnerlr
i
8 2%)
Figure 10. Effect of recycle flowrate.
study already reported (Reddy and Husain, 1978), which also tallies with the results of Gains (1977,1979) for converters operated at lower pressures. The desirability of operating with an exceas of nitrogen over the stoichiometric ratio has been recognized in the plant operation as demonstrated in an article by Gremillion (1979). Figure 9 shows effect of the loop pressure at the optimum H2/N2ratio and the recycle rate and inerts concentration fixed around the conditions of data set 3. As expected, an increase in pressure is favorable to the reaction equilibrium; however, the fractional H2conversion becomes asymptotic while the two other curves clearly indicate an optimum pressure at the fixed conditions. Figure 10 points out an optimum recycle flowrate with respect to the
Ind. Eng. Chem. Process Des. Dev., Vol. 21, No. 3, 1982 365 - 1 13
-
I I2
-
I I!
-
! 10
4
-109 -103
-
! 07
-
106
1.8
- 105 -104
2
-103
2
7102
d
