SIMULATION AND OPTIMIZATION OF AN INDUSTRIAL PSA UNIT
C.Barg^{1}, J.M.P.Ferreira^{2}, J.O.Trierweiler^{1 }and A.R.Secchi^{1* }^{1}Departamento de Engenharia Química, Universidade Federal do Rio Grande do Sul,
Rua Sarmento Leite, 288/24, 90050170, Porto Alegre  RS, Brazil
^{2}COPESUL, Companhia Petroquímica do Sul, III Pólo Petroquímico, Triunfo  RS, Brazil
(Received:November 26, 1999 ; Accepted: April 6, 2000)
]]>
Abstract  The Pressure Swing Adsorption (PSA) units have been used as a low cost alternative to the usual gas separation processes. Its largest commercial application is for hydrogen purification systems. Several studies have been made about the simulation of pressure swing adsorption units, but there are only few reports on the optimization of such processes. The objective of this study is to simulate and optimize an industrial PSA unit for hydrogen purification. This unit consists of six beds, each of them have three layers of different kinds of adsorbents. The main impurities are methane, carbon monoxide and sulfidric gas. The product stream has 99.99% purity in hydrogen, and the recovery is around 90%. A mathematical model for a commercial PSA unit is developed. The cycle time and the pressure swing steps are optimized. All the features concerning with complex commercial processes are considered.
Keywords: Pressure swing adsorption, simulation, optimization, gas separation.
INTRODUCTION
]]>The pressure swing adsorption (PSA) systems have received increasing attention from the specialized literature. Several articles have been published in the past years, with focus on the numerical simulation of such units and on experimental results. In general, these studied units are very simple, and have too many simplifications when compared to industrial units. This fact reflects the required computational complexity when working with real industrial systems, mainly when the objective is to optimize some process parameters.
In the past years, some works on optimization of pressure swing adsorption units have been done. Almost all those works are related to simple experimental units. Kvamsdal and Hertzberg (1995a) have worked with three systems of trace separation to investigate the effect of mass transfer during the blowdown step. They presented another work (Kvamsdal and Hertzberg, 1995b) on optimization of a trace separation system, where the importance of model simplifications to save time was emphasized, and the adsorption pressure, purge rate and column length were optimized. Chlendi et al. (1995) have proposed a method to characterize the behavior of pressure swing adsorption cycles based on experiments design and polynomial fittings to represent the system. The objective of this work was to make the problem easier to be optimized. There are some studies (Park et al., 1998; Kumar, 1994) on the optimal height of different kinds of adsorbents to separate multicomponent mixtures. Nilchan (1997) presented a more systematic approach to solve optimization of periodic adsorption process. It is proposed the discretization of the time domain with the inclusion of cyclic steady state conditions. The main disadvantage of this method is the increase in the computer memory requirement.
In the present work, modeling and simulation of an industrial hydrogen purification unit is made. The system is fed with a stream containing about 95% hydrogen, almost 5% methane, and traces of carbon monoxide, humidity and benzene. The product has about 99.99% purity in hydrogen, and the unit recovers about 90% of the hydrogen fed. The beds have three adsorbent layers, at feed entrance there is one composed of alumina, followed by an activated carbon layer, and after that, there is a zeolite layer. The alumina is used as a guard bed, in order to prevent irreversible contamination of the other adsorbents with water or benzene. This occurs because of the strong interaction of these heavier compounds with the adsorbents. As alumina has weaker adsorption strength, they can be dessorbed by pressure reduction, which would be impossible for the activated carbon or for the zeolite.
This work does not consider the presence of water or benzene in the feed, and the adsorption of carbon monoxide or methane in alumina is not considered too. Thus the alumina bed have no influence in the model results. Because of the weak interaction of hydrogen with all adsorbents, it is treated as an inert.
First of all, a model to describe an industrial pressure swing adsorption unit for hydrogen purification was made. All the process characteristics are considered for the model conception, in order to describe the system behavior correctly. The model was then implemented in gPROMS^{â} , and some real industrial cases were simulated. After verifying that the model is in accordance with plant data, the influence of some process parameters were studied, enabling the identification of optimal points. The main focus is the influence of the pressure swing steps on the process performance, which is represented by the relation of two characteristic times, T1 and T2. In the industry, for this particular system, the cycle time is manipulated to control the product purity, but the relation of the pressure swing steps is set to constant value. The objectives of this work are to study the influence of variations on the relation of the pressure swing steps and to identify its optimal value.
PROCESS DESCRIPTION
]]> The hydrogen purification plant has six vessels and operates in twelve steps. Each vessel is filled with three different adsorbents, alumina, activated carbon and zeolite. The sequencing of the different steps follows a predefined rule, which is best understood in Fig. 1.
The steps performed for each bed in a complete cycle are the following: pressure equalization 1 (provide), hold, pressure equalization 2 (provide), pressure equalization 3 (provide), provide purge, blowdown, purge, pressure equalization 3 (receive), pressure equalization 2 (receive), pressure equalization 1 (receive), repressurization, and adsorption. Figure 1 shows the duration of each step, which is a composition of two characteristic times, T1 and T2. The pressure equalization 1 has a fixed duration of 25 seconds, thus the hold step take (T1 – 25) seconds, and the repressurization step take (T2 + T1 – 25) seconds. A typical pressure variation with time is shown in Fig. 2.
]]> The purge is the main regenerating step, which enables the dessorption of the adsorbed material. The other steps are used basically to reduce the bed pressure to the purge pressure, and, after the purge step, to increase it to the adsorption pressure. The equalization steps are used to improve recovery, because they use the gas that is living a bed at a higher pressure and has to be depressurized to fill a bed that is increasing pressure to enter adsorption step. It is obvious that in the pressure equalization step the receiving pressure vessel cannot reach the adsorption pressure, once the final pressure of that bed must be lower or equal to the final pressure of the providing pressure vessel. In order to reach the adsorption pressure, the bed receives a fraction of the product gas in the repressurization step.
Figure 3 represents a brief description of the unit. The valves marked are open. This diagram shows the way the beds are connected one to another. The valves have to open and close synchronized to change the system from one step to another, in accordance with the necessary connections between the beds. The different bed connection schemes can be known by the analysis of Fig. 1.
The pressure variation in the variable pressure steps has a linear dependency with time. There is a control system that controls the interconnection valves in order to follow a predefined linear set point. In the diagram described in Fig. 3 the bed 1 is equalizing pressure with bed 3, bed 6 is equalizing pressure with bed 4, bed 5 is in blowdown, and bed 2 is in adsorption step. This configuration corresponds to the first step in the Fig. 1.
MATHEMATICAL MODEL AND OPTIMIZATION STRATEGY
]]>The PSA is a distributed system that can be modeled using a set of partial differential and algebraic equations. Although it has an intrinsically dynamic behavior, it can achieve a cyclic steady state after a certain number of cycles. A mathematical model for a complex commercial PSA unit is developed. The assumptions adopted are summarized below.
(1) The flow model adopted is the plugflow with axial dispersion. No radial gradients are considered.
(2)The system is nonisothermal, with thermal axial dispersion. Is considered local thermal equilibrium between the gas and the solid phases.
(3)To maintain the model generality, the system is considered a bulk separation process, then the change of velocity due to adsorption (or desorption) is taken into account by the overall material balance.
(4)The multicomponent adsorption equilibrium is computed adopting the extended Langmuir model. The isotherm constants are taken from literature, obtained from experiments with a single component (Park et al., 1998). No consideration was made concerning interaction among different molecules.
(5) The adsorption of hydrogen is considered negligible.
(6)A linear driving force model is adopted to compute the mass transfer dynamics, with constant overall mass transfer coefficient.
(7) Darcy´s Law is used to compute the pressure drop across the bed.
(8)The ideal gas law is assumed.
]]> (9) A linear time dependence of the pressure in the pressure equalization, provide purge, and repressurization steps is known from plant data.(10)All transport parameters, as well as physical properties of gas and solid phases are taken from the literature or estimated from empirical correlations.
(11)Heat transfer to the surroundings is negligible.
The different layers obey the same balance equations, but with different sets of physical properties and equilibrium parameters. Thus a balance equation must be done for each adsorbent layer.
The adsorption isotherm parameters are taken from Park et al. (1998).
Balance Equations
With the preceding assumptions, the balance equations to model the system are given below.
]]>a) Overall material balance applied to the gas phase, where the latest term on the right hand side of the equation takes into account the adsorbed quantity:
(1) 
where r is the gas density, v is the gas interstitial velocity, e is the bed porosity, q_{i} is the amount of component i adsorbed, z is the axial dimension, and t is the time.
b) Component material balance, for the gas phase:
(2) 
where n is the number of components, D_{z} is the axial dispersion coefficient, and C_{i} is the molar concentration of component i in the gas phase.
]]>c) Energy balance:
(3) 
where Cp_{g} is the mean heat capacity of the gas, Cp_{s} is the heat capacity of the adsorbent, K_{z} is the thermal axial dispersion coefficient, T is the temperature, and DH_{ads,i} is the heat of adsorption of component i.
d) Adsorption isotherm for adsorbed components:
(4a) 
]]> 
(4b) 
(4c) 
where a_{1,i}, a_{2,i}, b_{0,i}, b_{1,i} are the adsorption isotherm parameters for the component i in each adsorbent material, P_{i} is the partial pressure of component i in the gas phase, and q_{i,eq} is the amount of i adsorbed in equilibrium with the gas phase partial pressure of i at system temperature in that point.
e) Linear driven force equation, to model the mass transport between gas and solid:
(5) 
f) Darcy's equation:
(6) 
where P is the pressure, m is the viscosity, and dp is the mean particle diameter.
g) Ideal gas equation:
(7) 
Initial conditions
The solution of the previous equations needs some initial conditions to solve the equations with temporal derivatives. The bed is considered initially clean, and the pressure is equal to the initial pressure of the first step of the bed. The temperature is assumed to be equal to the feed temperature.
The following initial conditions are used to solve the set of differential and algebraic equations.

(8) 
where T_{F } is the temperature of the feed stream.
]]> Boundary ConditionsThe balance equations are applied to all the steps of the process. The differences from one step to another are accounted by the boundary conditions. On pressure variable steps, the pressure changes are assumed to have a linear dependency with time, here represented by P(t). This linear dependency is used as a boundary condition on the pressure. The boundary conditions used in the model are shown in Eqns. 9 to 14.
Fluid entering a bed:
(9) 
where X_{i} is the molar fraction of component i in the gas phase, X_{i,in} is the molar fraction of component i and T_{in} is the temperature in the stream that is entering the bed.
Fluid leaving bed and closed end:
]]> 
(10) 
For pressure and velocity, the boundary conditions are the following:
Adsorption
(11) 
where v_{f} is the interstitial velocity at feed entrance, and P_{ads} is the adsorption pressure.
Pressure equalizations, repressurization and provide purge:
]]> (12) 
where P(t) is the previously known pressure variation with time.
Blowdown:
(13) 
Purge:
(14) 
where P_{D} is the purge pressure, P_{out} is the pressure and v_{out} is the interstitial velocity in the exit of the bed which is providing purge gas.
]]>Intersection of the Layers
The intersection of the layers is modeled with a set of continuity equations. In order to avoid numerical problems when calculating the derivatives near the border between the two adsorbent layers, the axial dimension is divided in two parts, related to each adsorbent.
Then the derivatives near the borderline are calculated without using values of the other axial distribution domain, because the software calculates this derivatives with backward or forward finite difference formula, in exception of the centered finite difference approximation used on the points inside the interval. This way, there is a separation of the two axial distribution domains. They must be connected by a set of continuity equations (Eqns. 1516).
(15) 
(16) 
Optimization Strategy
Due to the mathematical complexity of the problem, the strategy selected to achieve the optimal point is the simplest possible. The complete discretization method, described in the work of Nilchan (1997) has required an unavailable amount of computer memory, because of the elevated number of variables generated by the time discretization.
The objective is to study the combined effect of the cycle time, here represented by the T1 + T2 parameter, and of the T1/T2 ratio on the system behavior. A mesh with different values for these parameters was generated, as can be seen in Table 1, resulting in twentyfive different simulation cases.
]]> With the results of product purity and recovery for these cases it is possible to find an optimal value for the two parameters for this specific case study.
Numerical Solution of the Model
To solve the set of partial differential and algebraic equations (PDAE) the gPROMS^{â} (Pantelides, 1996) software is used. This program is a generalpurpose process simulator that can work with models involving distributed systems.
The set of equations has to be changed when a bed goes from one step to another, due to the change on the boundary conditions. This work is also executed by the simulator. The user should just implement a routine that defines when the step changes should be made. This way, any modification to be done in the simulation can be made quickly.
The spatial discretization of the set of partial differential and algebraic equations (PDAE) results a set of differential and algebraic equations (DAE) which is integrated in relation to time using the SRADAU code.
The main advantage in using gPROMS^{â} is the reduction of the time spent to implement the model, once the simulator does several mathematical manipulations, as the spatial discretisation, the changing among different steps and the temporal integration of the resulting set of differential and algebraic equations.
The model was solved in a PC with Pentium II^{®} 350MHz processor. Each domain is discretized with the centered finite difference approximation, using fourth intervals. The simulations took about one fifth of the plant time, and have required about 75Mb of main memory. Thus to simulate one cycle of about three thousand seconds, it spent about six hundred seconds of computer time.
]]>
RESULTS AND DISCUSSION
The first objective of the work is to reproduce the plant behavior by numeric simulation, and then to find the optimal cycle time and optimal ratio of the characteristic times (T1 and T2). The operational parameters optimized in this study are the sum of the characteristic times T1 and T2, and its ratio, T1/T2.
The analysis of the plant behavior is made with two main parameters, the product purity and the hydrogen recovery.
The product purity is defined as the ratio between the molar amount of hydrogen collected in the product during the last cycle and the total molar amount of the product.
(17) 
where T_{cycle} is the cycle time, R_{i} is the adsorbent vessel internal radius, e is the bed porosity, V_{p} is the gas velocity at product end, and X_{H2, p} is the hydrogen molar fraction at product stream.
]]> The hydrogen recovery is defined as the ratio of the molar amount of hydrogen collected in the product during the last cycle and the molar amount of hydrogen fed to the bed in the last cycle.
(18) 
where X_{H2,f} is the hydrogen molar fraction at feed stream.
Thus the product purity and the hydrogen recovery are calculated at the beginning of a cycle, remaining constant until another cycle starts. Table 2 shows that the product purity and the hydrogen recovery in the simulation agree with plant data, for the same operational conditions.
]]> A basic case was first simulated. Figure 4 shows the product purity and hydrogen recovery of the system as a function of time. The plant is assumed to start with a hundred percent purity and recovery. This simulation is made with the cycle time equal to 5568 seconds. The system is considered to achieve the cyclic steady state (CSS) when the difference of product purity between one cycle and the proceeding is less then 10^{5}.
The cyclic steady state is achieved, and then used as initial condition to the other simulations. All perturbations on the values of T1 and T2 are done as step variations, after initializing the system with the cyclic steady state with T1 + T2 equal to 928 seconds and T1/T2 equal to 0.48. Then the values of T1 and T2 are 302 and 626 seconds respectively, and the cycle time is 5568 seconds. The variations of the cycle time are based on variations of ± 30%, and the variations on T1/T2 are chosen to cover a wide operational range.
Figure 5 shows the variation of hydrogen recovery and product purity with the parameter T1 + T2, which reflects the value of the cycle time. The results showed that product purity decreases with the increase of the cycle time, while the hydrogen recovery presented low reduction. This relationship among cycle time, recovery and purity is qualitatively well known. The utility of this kind of work is to identify the optimal cycle time in order to avoid product overspecification. The optimal cycle time would be the value that brings the product inside specification limits with the greatest possible recovery. The optimal cycle time may be related to different feed conditions.
]]>
In the industrial unit there is a feed forward control to automatically change the cycle time when the feed flow rate is changed, acting in a linear way. Such relationship assumes that, reducing the feed flow rate in fifty percent, the cycle time will be increased by two times. A very interesting study would be the investigation of the variation of the optimal cycle time with the feed flow rate, the feed temperature, and the adsorption pressure, in order to verify if the optimal value is really inversely proportional to the amount fed to the unit. Such study is going to be done as a continuation of this study, once the investigation of the optimal value for the ratio T1/T2 seems to be unnecessary.
Figure 6 shows the variation of the product purity and the hydrogen recovery with the ratio T1/T2 for two different values for the T1 + T2 parameter. The variations are verified only in the fourth decimal.
The variations of the purity and recovery showed a very low dependency with the variations made in the ratio T1/T2. It can be attested that, in this specific case, variations on the ratio T1/T2 cause no effect on the product purity or on the hydrogen recovery. In fact, the amount of gas used to purge the bed is the same, because the difference of the initial and the final pressure of the provide purge step is equal, only the rate of variation of the pressure with time is changed. The difference on the time to depressurize the bed causes different purge gas velocities, which presented no effect on the product purity or hydrogen recovery in this case. The choice of the ratio T1/T2 to be used is probably made based on mechanical aspects of the system.
A fact that is also observed is the simultaneous increase in hydrogen recovery and product purity when the ratio T1/T2 is reduced. Obviously, if the cycle time remains unchanged and the hydrogen purity is increased, then the additional hydrogen transferred to the product stream will reflect in reduction of the amount of hydrogen purged, increasing also the recovery.
]]> CONCLUSIONS
A model to simulate and optimize an industrial PSA unit was developed. All the features concerning with a complex industrial plant were considered. The characteristics of the process were well represented by a nonisothermal model using the extended Langmuir isotherm to represent the adsorption equilibrium and a linear driving force model to describe the adsorption dynamics.
The effects of the cycle time and the ratio T1/T2, two characteristic parameters of the process, were studied. The influence of the cycle time on product purity and hydrogen recovery was studied, showing that it is possible to foresee the optimal cycle time to achieve a defined product specification with maximum recovery. The effect of the ratio T1/T2 showed a negligible influence in product purity and hydrogen recovery.
REFERENCES
Chlendi M., Tondeur D., Rolland F. (1995), "A method to obtain a compact representation of process performances from a numerical simulator – example of pressure swing adsorption for pure hydrogen production", Gas Sep. Purif., vol. 9, 125135. [ Links ]
Kumar R. (1994), "Pressure swing adsorption process: Performance optimum and adsorbent selection," Ind. Eng. Chem. Res., vol. 33, 16001605. [ Links ]
Kvamsdal H. M., Hertzberg T. (1995), "Optimization of pressure swing adsorption systems – The effect of masstransfer", Chem. Eng. Sci., vol. 50, 12031212. [ Links ]
Kvamsdal H. M., Hertzberg T. (1995), "Pressure swing adsorption – Optimization of a trace separation system." Comp. Chem. Eng., vol. 19, S339S344. [ Links ]
Nilchan, S. (1997), "The Optimization of Periodic Adsorption Processes", PhD. Thesis. Imperial College. [ Links ]
Pantelides C.C. (1996), "An Advanced Tool for Process Modelling, Simulation and Optimisation," Presented at Chemputers Europe III, Frankfurt. [ Links ]
Park J., Kim J., Cho S., Kim J. and R. T. Yang. (1998), "Adsorber dynamics and optimal design of layered beds for multicomponent gas adsorption," Chem. Eng. Sci. , vol. 53, 39513963. [ Links ]
*To whom correspondence should be addressed
]]>