Ind. Eng. Chem. Res. 2001, 40, 3217-3224
3217
SEPARATIONS Fast Finite-Volume Method for PSA/VSA Cycle SimulationsExperimental Validation Richard S. Todd, Jianming He, Paul A. Webley,* Christopher Beh, Simon Wilson, and Matthew A. Lloyd Department of Chemical Engineering, P.O. Box 36, Monash University, Victoria, 3800 Australia
In a previous study (Webley, P. A.; He, J. Fast Solution-Adaptive Finite Volume Method for PSA/VSA Cycle Simulation. 1. Single Step Simulation. Comput. Chem. Eng. 2000, 23, 1701), a fast, solution-adaptive finite-volume technique for the simulation of a single adsorption step for a variety of boundary conditions was described. In this study, we apply this method to the simulation of nonisothermal PSA/VSA cycles of general complexity. Using successive substitution, stage-wise node refinement, and control algorithms, rapid cyclic steady state can be achieved even for problems with very long dynamics. We illustrate the advantages of node refinement as a useful tool for accelerating convergence to cyclic steady state. Simultaneous incorporation of control techniques into the successive substitution framework allows for rapid convergence of the PSA process to design specifications. We compare the results of our simulator to experimental data for a two-bed VSA process and find good agreement in pressures, flows, and temperatures. Introduction Adsorption has become a common and efficient unit operation for the separation of multicomponent gas streams in the chemical industry. Common applications for adsorption include oxygen and nitrogen enrichment from air, H2 purification, CO2 removal, and natural gas sweetening and drying.2 Given the large number of process variables associated with a PSA-based process for gas separation, setting up a comprehensive experimental program to determine the optimal operating conditions is economically unfavorable and time-consuming. For this reason, the use of numerical simulators based on the governing conservation equations for PSA systems appears to be the logical approach. However, the complexity of the governing equations for such processes usually results in computation times that are not far removed from those required experimentally, even when many limiting assumptions are made. Under these conditions, the advantages of a numerical approach over an experimental approach are not clear. For a numerical approach to prove worthwhile, an accurate model covering all fundamental principles of operation must be used, and computation times significantly faster than those achievable experimentally must be attained. Over the course of a real PSA cycle, the state (mole fractions, temperature, and pressure) that exists at each position within the adsorbent bed at any time is a complex function of the boundary conditions and the physical processes occurring simultaneously within the adsorbent bed. Bulk gas motion between each adsorbent particle, diffusive gas motion through randomly oriented porous channels of the adsorbent particle, and nonisothermal effects associated with the adsorption of a gas * Author to whom correspondence should be addressed. E-mail:
[email protected]. Fax: +61 3 9905 5686. Tel.:+61 3 9905 1874.
onto the adsorbent surface are just a small collection of the important physical processes to be considered. Equilibrium models, although they provide useful insight into the important process and design parameters, are not sufficiently accurate to allow for PSA design. Numerical methods employed to solve the resulting set of conservation equations for the adsorbent bed include orthogonal collocation,3 the method of lines,4 and finite-volume analysis,5 among others. These models discretize the spatial domain of the adsorbent bed and then solve the resulting time derivatives through an appropriate ODE solver. This procedure involves successive substitution, where bed conditions from the previous cycle become initial conditions for the next simulated cycle, continuing in this manner until cyclic steady state (CSS) is achieved. To overcome the slow convergence inherent in a successive substitution method, various acceleration procedures employing direct determination have been proposed. Most recently, Ding and LeVan6 demonstrated acceleration toward CSS for PSA and TSA simulations using hybrid Newton-Broyden methods. Smith and Westerberg7 considered a shooting method that explicitly combines the CSS condition within the process model. An alternative to successive substitution is direct determination. Harriott8 and Nilchan and Pantelides9 proposed complete discretization of both the temporal and spatial domains. By enforcing periodicity, the problem is reduced to a set of nonlinear equations that is solved through appropriate (frontal) methods. Although these methods might be favorable computationally for simple problems, direct determination has yet to be applied to a complex nonisothermal simulation. Most likely, a combination of successive substitution and direct determination could yield substantial overall performance improvement and robustness. A common assumption with small-diameter metalwalled columns of an experimental apparatus is that
10.1021/ie0008070 CCC: $20.00 © 2001 American Chemical Society Published on Web 06/09/2001
3218
Ind. Eng. Chem. Res., Vol. 40, No. 14, 2001
near-isothermal behavior prevails. Despite the damped thermal effect provided by small-diameter, metal-walled columns, neglecting thermal effects can significantly alter the final CSS solution. This is evident in the slow evolution to CSS, on the order of hundreds to thousands of cycles, for processes where nonisothermal conditions prevail.10 Although heat leaks might appear small, thermal convection over a cycle is often comparable to small heat leaks, and accumulation over thousands of cycles can lead to appreciably different final temperature profiles if heat leaks are not included. An in-depth analysis of the various simulators proposed by several authors reveals a lack of generality in approach to system modeling. The importance of a nonisothermal model becomes critical when a valid assessment between laboratory-scale units and an industrial-sized plant is to be made. A general model should capture not only operating characteristics associated with pressure, temperature, and composition changes, but also those associated with the operating performance of ancillary equipment such as feed compressor(s), vacuum pump(s), and valves located around each bed; thus, a large variety of boundary conditions must be provided. Another factor is the importance of pressure drop in the model, without which it is difficult to assess the effect of cycle time, especially for RPSA. The recent abundance of patents detailing different adsorbent structures such as monoliths and laminates represents a direct effort to reduce the pressure drop, which primarily leads to excessive power requirements. This paper extends the numerical model developed earlier1 for a single adsorption step to now simulate a general cyclic adsorption process. The use of a cyclicsteady-state acceleration procedure, termed node refinement, is also discussed within the framework of our cycle simulator. Our adsorption simulator is called MINSA for Monash Integrator for Numerical Simulation of Adsorption. Adsorption Model and Cycle Simulation Adsorption Model. The adsorption model used here is similar to that presented previously,1 and we do not repeat the related equations here. We include the energy balance in a manner identical to that described in Kumar et al.,11 with the assumption that the local gas and solid temperatures are equal. Conditions under which this assumption is true have been identified12 and cover the region of interest in our work. The isosteric heat of adsorption is calculated using the ClausiusClapeyron equation and the adsorption isotherm. To avoid computation of this derivative for each time step and each spatial node, a preprocessing calculation is made over a wide range of temperatures and partial pressures. Appropriate regression is then used to express the isosteric heat as a function of loading and temperature as the two independent variables. Pressure drop is modeled explicitly using the Ergun equation, whereas mass transfer is modeled using a linear driving force model based on solid loading as first established by Glueckauf and Coates.13 We use either a linear, single-site or a dual-site Langmuir model (as needed to represent the data) to describe the adsorption isotherm. Heat transfer to and from the bed and wall is important in accurately representing the thermal profiles in the bed, especially for bulk-gas separation processes, and we use the model of Li and Finlayson14 to represent the heat transfer coefficient. A more complex wall model (appropriately discretized) should
be used for systems in which the heat transfer up/down and through the wall represents an appreciable fraction of the total thermal convection in the beds. Our conservation equations are discretized in space using a solution-adaptive finite-volume method described earlier.1 In addition to solution of the simultaneous ODEs generated from the finite-volume discretization of the conservation equations, we add ODEs to the system to represent the mass and energy flows entering and leaving the bed. This allows us to calculate accurately the overall energy and mass balance for the cycle, parameters that we use to help us evaluate cyclic steady state (CSS). The simultaneous ODEs are solved in the time domain using LSODA or DVODE, a stiff, first-order, backward difference integrator package.15 Cycle Simulation. From a mathematical viewpoint, there is no need to simulate each bed in a multibed PSA system separately unless accurate dynamic data are needed (e.g., PSA start-up and short-time upsets) as each bed undergoes the same cycle separated in time. We therefore simulate a single bed through the sequence of steps required by the cycle. Bed-bed interactions, such as pressure equalization and purge, are dealt with by temporary storage of the conditions of the connecting stream, as done in previous studies.11 At the end of a cycle, the bed conditions are used as the initial conditions for the next cycle so that a successive substitution approach is used. Although this approach converges linearly at best, it does have advantages in elucidating the approach to cyclic steady state as well as providing useful information on criteria used to establish whether cyclic steady state has been achieved. It is important to note that our simulator calculates the bed pressure profiles directly from the imposed boundary conditions, e.g., valve settings, velocities, molar flow rates, etc. This has the benefit of mimicking actual the operation of a real PSA and providing bed pressure profiles as output for direct comparison to experimental values. Our simulator is designed to model any one of four boundary conditions, as described previously in detail.1 The four boundary conditions are flow through a valve of fixed valve coefficient (Cv), fixed velocity flow, fixed mole rate, and fixed pressure. These options allow us to simulate a wide variety of experimental devices such as pumps, mass flow controllers, etc. Cyclic Simulation and Steady-State Determination Successive Substitution and Node Refinement. As mentioned earlier, the approach used to converge to cyclic steady state is successive substitution. This method is slow to converge, especially when adiabatic conditions are simulated and convectively driven thermal waves are moving through the system. Because of the temperature dependence of the loading of species i, the slow temperature evolution means that bed profiles also follow this slow transient until the onset of cyclic steady state. As a result, our simulator must progressively iterate toward a steady-state solution through the successive substitution of very small perturbations in the state vector for each cycle. Neglect of this slow temperature transient results in a cyclic-steady-state solution well removed from that obtained experimentally. Despite these drawbacks, the successive substitution method is extremely robust and provides information that is used to determine whether cyclic steady state has been achieved.
Ind. Eng. Chem. Res., Vol. 40, No. 14, 2001 3219 Table 1. Set of Criteria Used to Determine Cyclic Steady State notation Emax
Enorm
physical meaning
calculation
maximum difference in each value of state vector X between cycles k and k - 1; variables are f1 ) P, f2 ) y, f3 ) n1, f4 ) n2, and f5 ) T
RMS error between cycles k and k - 1; M is the number of spatial points of comparison between the two cycles (for comparison, at the end of each cycle for N finite volumes, M ) N)
( max|f ki - f k-1 | i i X ) {f1, f2, ..., f5}
where
{
}
1/2
M
∑(f
k i
- f k-1 )2 i
i)1
M
E1 E2 Eheat
relative mass balance error for component 1 over one cycle relative mass balance error for component 2 over one cycle relative enthalpy balance error over one cycle
(moles of 1 out)/(moles of 1 in) (moles of 2 out)/(moles of 2 in) (enthalpy out)/(enthalpy in)
Ediff
average gradient of Eheat from cycle k - r to cycle k
Ekheat - Ek-r heat r
If the successive substitution is to be retained, a method for accelerating convergence to cyclic steady state becomes highly desirable. The approach adopted with our simulator is one of successive node refinement. One of the advantages of employing a high-order, finitevolume approach to discretization of the conservation equations is that relatively stable solutions to highly convective flow regimes can be achieved with a small number of node points.1 Also, the conservative nature of the numerical method means that mass and energy balance conservation can be achieved within machine precision, independent of the level of discretization. This means that we can run our simulation with a low number of nodes to cyclic steady state, as determined by closure of the mass and energy balances. The simulator commences iterations from an initial condition with only a small number of internal node points (usually less than 8) and runs until cyclic steady state is achieved. Cyclic steady state is determined by comparing bed profiles for temperature, pressure, and composition from one cycle to the next (using error norms), along with mass and energy balance closure, against predefined tolerances (see discussion below for more details). These tolerances can be quite generous at this stage of the simulation. The low number of nodes during this stage means that cycle simulation is extremely rapid (∼1-2 s/cycle on a DEC Alpha 500au workstation), and cyclic steady state is achieved within minutes. In addition to these criteria, we also require that design specifications be met to within a predefined tolerance at cyclic steady state (see Control Techniques below for a discussion on achieving design specifications). Once cyclic steady state is established, performance parameters, such as production rate, product purity, recovery, etc., are calculated. The next stage involves a finer level of discretization (typically 12 nodes) and interpolation of the nodal values for all state variables given in the profile from the previous node stage to the new node stage. Here, cubic-spline interpolation is used over each profile to determine new nodal values. Cyclic steady state for the previous node scheme then becomes the initial conditions for the new node procedure. From here, the simulation continues with the increased number of nodes until, again, cyclic steady state is achieved within a predefined tolerance for bed profiles and mass and energy balance closure. This process is repeated as many times as necessary to reduce the stage-to-stage variation in process performance to below a prescribed value. It should be emphasized that the appropriate number of nodes (or elements) required to
capture the physics of the adsorption problem is not known a priori and is likely to be different for different operating conditions, isotherms, etc. Simply assuming that a set number of nodes (e.g., 30) is sufficient does not guarantee that the correct solution will always be obtained. The node-refinement technique described above, however, ensures that an appropriate number of nodes will always be automatically used for any desired level of accuracy. Cyclic-Steady-State Criteria. The state vector X of dependent variables for the system contains pressure, mole fraction of one component, loadings of each component, and temperature. At cyclic steady state, at the end of a cycle, X returns identically to its value at the start of the cycle. This condition is mathematically sufficient for determining CSS but, from a practical standpoint, requires very tight tolerances with correspondingly small ODE tolerances (less than 10-9), resulting in excessive CPU time. The same result can be achieved more quickly if overall energy and mass balance errors are used as appropriate criteria (in addition to closure of X) with a reasonable relative ODE tolerance, e.g., 10-6. In summary, Table 1 shows the set of criteria that we use to establish automatically whether cyclic steady state has been achieved. We use the usual error quantities Emax and Enorm, as well as closure of the total energy and mass balances. The last criterion, Ediff, is required to avoid picking up the false zeros in the oscillating energy balance error evolution. Thus, Ediff enforces detection of a cyclic steady state only when the energy balance oscillation has “died out”. It is important that the simulator automatically detects the presence of CSS because evaluation of additional cycles is costly in CPU time. Achieving Design Specifications During operation of a typical PSA plant, it is important that certain specifications on system performance be controlled to within predefined tolerances. These constraints become important in the simulation of a real industrial process, where control schemes are often applied to particular operating variables for control on system performance. Some typical quantities controlled within a PSA process include purity of the product, the lowest pressure achieved during the blowdown step, the highest pressure achieved during the feed step, and intermediate pressures at a specific times. In the control of these (and other) parameters, boundary conditions, such as mole rate to a bed, inlet/outlet velocity or volumetric flow rate to a compressor, and valve positions
3220
Ind. Eng. Chem. Res., Vol. 40, No. 14, 2001
Figure 1. Schematic of experimental VSA system.
between the bed and tank(s), become the manipulated parameters. Control Techniques. In practice, valve openings, flow rates, etc., are manipulated via PID algorithms to ensure that target purities, pressures, and flow rates are achieved at cyclic steady state. Mathematically, there is no particular reason that similar schemes should be used to achieve desired cyclic steady states in a model of the process. Indeed, the use of standard iterative schemes such as modified Newton’s methods or Broyden methods should yield improved performance if the only goal is to achieve cyclic steady state.16 However, our experience with these iterative schemes has indicated that the strong nonlinearity of the process leads to problems with convergence and stability, particularly in purity control. Also, for each iteration, a cyclic-steady-state solution is required, thus increasing the time for each iteration. In contrast, schemes that adjust control variables on each cycle (such as PID loops) are slower to converge but are able to act simultaneously with the approach to cyclic steady state. Overall, we have experienced far more success with control loop implementation than with other iterative schemes. To implement mathematical PID control, the relationship between the controlled variables and the manipulated variables must be knownsin our case, it is our numerical solution of the PSA model equations. The control algorithm implemented within our simulator takes the form
ck+1 ) cki + Kc[(eki - ek-1 ) + Tieki + Td(cki - ck-1 )] (1) i i i where Kc, KcTi, and KcTd are proportional, integral, and derivative constants, respectively, and updates in the manipulated variables ci are made from cycle to cycle. The errors, ei, are differences between target and set points for each controlled variable. The disadvantage of the above scheme is that a mapping must be made between manipulated and controlled variables before the start of the cycle simulation, which requires some knowledge or experience of the adsorption system
Table 2. CPU Time Comparison with and without Node Refinement case
N ) 10
N ) 20
N ) 30
relative CPU time
node refinement no node refinement
0.04 -
0.35 -
0.61 2.58
1.0 2.58
behavior. The control parameters Kc, Ti, and Td are userspecified but can be changed during cycle simulation. As part of the criteria for cyclic steady state, we require that the target variables be within a specified tolerance of the setpoints. Illustration of Acceleration from Node-Refinement Technique We illustrate our simulation scheme with a four-step VSA cycle for the production of oxygen-enriched air over 5A zeolite. Each bed undergoes only four distinct phases: (a) high-pressure feed step with product takeoff, 25 s; (b) countercurrent evacuation of the bed through the feed end, 20 s; (c) countercurrent receipt of purge gas from another bed, 10 s; and (d) receipt of repressurization gas from product tank, 5 s. We first simulate the system without using the noderefinement procedure. Thirty internal node points are set through the bed, and the simulation is allowed to proceed to CSS. We then use node refinement beginning with 10 nodes and then stepping to 20 and finally 30 nodes. Table 2 shows a comparison of the relative CPU times for these two cases. We see that the node-refinement scheme results in an improvement of more than a factor of 2 in overall CPU time to produce identical results. For systems exhibiting very slow evolution to cyclic steady state, much of the work done by the simulator in getting to CSS can be done with a coarse node discretization, resulting in substantial savings in CPU time. Similar performance gains are not expected for rapidly converging systems, but it is hard to conceive of a system for which the node-refinement scheme would perform worse than no node refinement. The only additional calcula-
Ind. Eng. Chem. Res., Vol. 40, No. 14, 2001 3221 Table 4. Physical Properties of NaX Zeolite Bed and Adsorbent Properties bed diameter ) 0.1023 m (measured) dp ) 1.97 mm (equivalent diameter since pellets are 1/ 16 in. × 2.5 mm) b ) 0.375 (Haefner and Thodos, 1986) Fb ) 647 kg/m3 (measured) Lcolumn ) 1.95 m (measured) mass of adsorbent (2 beds) ) 20.74 kg (measured) adsorbent heat capacity ) 920 J/(kg K) (Haefner and Thodos, 1986)
Adsorbent Isotherm UOP APG 13X pellets, dual Langmuir isotherm neq i )
m1bipyi 2
1+
∑ j)1
bjpyj
+
m2dipyi 2
1+
∑d py j
j
j)1
where bi ) b°1 exp[∆q1i /RT] and di ) d°i exp[∆q2i /RT] m1 ) 2.99 gmol/kg, m2 ) 4.00 gmol/kg bN2° ) 8.3198 × 10-5 bar-1, bO2° ) 2.2206 × 10-4 bar-1 dN2° ) 2.1218 × 10-6 bar-1, dO2° ) 1.194 × 10-4 bar-1 ∆qN21 ) 17 889 J/gmol, ∆qO21 ) 11 322 J/mol ∆qN22 ) 17 441 J/gmol, ∆qO22 ) 11 163 J/mol rigorous heat of adsorption model used Figure 2. Eight-step VSA cycle used to determine experimental data.
Mass Transfer Parameters mean pore diameter (measured) ) 0.3 µm, macropore model
Table 3. Experimental Conditions and Control Variables for High- and Low-Purity Runs
CSS Tolerances mass and energy Balance errors to be within 0.05% between two concurrent cycles T, y, and p profiles to be within 0.05% between two concurrent cycles
Process Conditions for Low-Purity Run end of step pressures step 1 (not all controlled) step 2 step 3 step 4 step 5 step 6 step 7 feed temperature 18 °C product purity 80.2% O2
Process Conditions for High-Purity Run end of step pressures step 1 (not all controlled) step 2 step 3 step 4 step 5 step 6 step 7 feed temperature 23 °C product purity 94.9% O2
126.9 137.3 106.1 96.2 37.2 51.5 66.0
132.9 142.0 109.2 93.8 34.5 50.2 68.5
Control Conditions step 1 time ) 12 s inlet flow determined by valve coefficient equal to that in following step step 2 time ) 10 s inlet flow determined by valve coefficient; end of step pressure controlled to 137.2 kPa by adjusting valve step 3 time ) 5 s fixed valve provide purge from bed A to B (top-top) step 4 time ) 3 s fixed valve provide pressure equalization from bed A to B (top-top) step 5 time ) 12 s bed A on pump down with fixed valve coefficient step 6 time ) 10 s bed A on pump down with fixed valve coefficient controlled so that end of pump-down pressure is required value step 7 time ) 5 s bed A on receive purge step 8 time ) 3 s bed A on receive pressure equalization
tions done during node refinement involve spline fitting and interpolation, both of which consume very little CPU time. Application of the Simulator After validating our simulator against the predictions of the binary linear isotherm model for PSA, as derived by Knaebel and Hill,17 we proceeded to examine the ability of our model to match experimental data on a pilot-scale VSA unit. We model the ternary air system
(nitrogen, oxygen, argon) as a binary system assuming that argon behaves similarly to oxygen.18 Thus, the feed air to our simulator is nitrogen (78%) and oxygen (22%). Experimental Apparatus. The predictions of our simulator were compared to experimental data generated by our two-bed pilot-scale vacuum swing adsorption unit. A flow diagram of the apparatus is shown in Figure 1. The system consists of two beds constructed from 3-mm-thick PVC of internal diameter 0.104 m. Each bed is packed with the zeolite(s) of interest such that the packed bed height is approximately 1.8 m. A valve manifold is connected to each the top and the bottom of the beds. The top manifold consists of a separate product valve for each bed, as well as a separate switch valve, which leads to a control valve to allow controlled gas exchange between the beds. The lower manifold contains separate feed and vacuum valves on each bed. Also, the feed line is equipped with a control valve as is the vacuum line and product line. Thus, control of feed, vacuum, and product amount is possible over a wide range. A feed and product tank is also used to buffer swings in supply and product pressure and flow. The feed air line also contains an adsorbent drier in which water and carbon dioxide are removed. The drier also has the effect of slightly enriching the feed air in oxygen so that we feed 22% oxygen + 1% argon to our VSA system. This is accounted for in our simulation by increasing the oxygen content to 23%. Each bed is equipped with several thermocouples located axially from the entrance to the exit. Pitot tube flow meters are located on the feed and vacuum lines, each with its own pressure and temperature instrumentation to allow for calculation of molar flow rates. The product flowmeter is of the orifice type. The product composition is measured with a Servomex 1440C analyzer, which is accurate to within 0.5% O2. All digital and analog input/output is controlled and logged through
3222
Ind. Eng. Chem. Res., Vol. 40, No. 14, 2001
Table 5. Comparison of Experimental and Simulator Data run 1 80% purity quantity
exp
calc
run 2 95% purity % error
4490 3876 580 80.2 47.2 21.40
4322 3846 476 80.2 39.4 17.60
step 1. feed repressurization step 2. feed + product step 3. provide purge step 4. provide PE gas step 5. evacuation step 6. receive purge gas step 7. receive PE gas
126.9 137.3 106.1 96.2 37.2 51.5 66
end of step pressures (kPa) 126.8 136.7 106.4 96.7 37.6 51.8 -0.6 65.1 1.4
step 1. feed repressurization step 2. feed + product step 5. evacuation step 6. receive purge gas mass balance error (%)
end of step feed and vacuum flows (mmol/bed) 1410 1503 -6.6 1444 835 658 21.2 545 1661 1647 0.8 1547 277 276.0 0.4 264 0.8
