Simulation and Optimization of Pressure Swing Adsorption Systems

Ind. Eng. Chem. Res. 2009, 48, 2327–2343

2327

Simulation and Optimization of Pressure Swing Adsorption Systems Using Reduced-Order Modeling Anshul Agarwal,†,‡ Lorenz T. Biegler,*,†,‡ and Stephen E. Zitney† Collaboratory for Process and Dynamic Systems Research, National Energy Technology Laboratory, P.O. Box 880, Morgantown, West Virginia 26507-0880, and Chemical Engineering Department, Carnegie Mellon UniVersity, Pittsburgh, PennsylVania 15213

Over the past three decades, pressure swing adsorption (PSA) processes have been widely used as energyefficient gas separation techniques, especially for high purity hydrogen purification from refinery gases. Models for PSA processes are multiple instances of partial differential equations (PDEs) in time and space with periodic boundary conditions that link the processing steps together. The solution of this coupled stiff PDE system is governed by steep fronts moving with time. As a result, the optimization of such systems represents a significant computational challenge to current differential algebraic equation (DAE) optimization techniques and nonlinear programming algorithms. Model reduction is one approach to generate cost-efficient low-order models which can be used as surrogate models in the optimization problems. This study develops a reducedorder model (ROM) based on proper orthogonal decomposition (POD), which is a low-dimensional approximation to a dynamic PDE-based model. The proposed method leads to a DAE system of significantly lower order, thus replacing the one obtained from spatial discretization and making the optimization problem computationally efficient. The method has been applied to the dynamic coupled PDE-based model of a twobed four-step PSA process for separation of hydrogen from methane. Separate ROMs have been developed for each operating step with different POD modes for each of them. A significant reduction in the order of the number of states has been achieved. The reduced-order model has been successfully used to maximize hydrogen recovery by manipulating operating pressures, step times and feed and regeneration velocities, while meeting product purity and tight bounds on these parameters. Current results indicate the proposed ROM methodology as a promising surrogate modeling technique for cost-effective optimization purposes. 1. Introduction Separation of gases accounts for a major fraction of the production cost in chemical, petrochemical, and related industries. There has been a growing demand for economical and energy efficient gas separation processes. The new generation of more selective adsorbents developed in recent years has enabled adsorption-based technologies to compete successfully with traditional gas separation techniques, such as cryogenic distillation and absorption. The last few decades have seen a considerable increase in the applications of adsorptive gas separation technologies, such as pressure swing adsorption (PSA); the applications range from bulk gas separations and drying to trace contaminant removal.1–3 In particular, from an economic and environmental point of view, hydrogen separation and carbon dioxide capture from flue gas streams are the most promising applications of PSA. For hydrogen separation, PSA has been found to be the most favorable process compared to other processes on the basis of flexibility, reliability and purityrecovery tradeoff.4–7 The essential feature of a PSA process is that the preferentially adsorbed species are removed and the bed is regenerated by reducing the total pressure. With extensive industrial applications, there is a significant interest for an efficient modeling, simulation, and optimization strategy. In H2 recovery with very high purity, PSA has been extensively studied to understand the dynamics of the processes.7–9 However, the design and optimization of PSA processes have largely remained an experimental effort because of the complex nature * To whom correspondence should be addressed. E-mail: lb01@ andrew.cmu.edu. † National Energy Technology Laboratory. ‡ Carnegie Mellon University.

of the mathematical models describing practical PSA processes.10 The mathematical models are described by coupled nonlinear partial differential equations (PDEs) distributed in space and time, with different initial and boundary conditions defining the steps of the process and high nonlinearities arising from nonisothermal effects. The computational effort required to solve such systems is usually quite expensive and prohibitively time-consuming. Besides this, stringent product specifications, required by many industrial processes, often lead to convergence failures of the optimizers. Sophisticated optimization strategies have been developed and applied to PSA systems with significant improvement in the performance of the process (for example, increased purity or recovery of the product gas). For the optimization of PSA and rapid PSA processes, Nilchan et al.11 proposed complete

Figure 1. A two-bed four-step PSA process.

10.1021/ie071416p CCC: $40.75  2009 American Chemical Society Published on Web 05/28/2008

2328 Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009 Table 1. Physical Properties and Basic Simulation Parameters parameter

value

bed length (L) bed porosity (b) bed radius (Rb) particle radius (Rp) particle porosity (p) diffusivity (Dx) particle density (Fp) bed density (Fb) thermal diffusivity (KL) heat capacity of solid (Cps) heat transfer coefficient (h) lumped mass transfer coefficient (k) heat of adsorption (∆H) gas viscosity (µ) R ambient temperature (Tw) feed temperature (Tfeed) feed composition feed pressure (Pfeed) purge pressure (Ppurge) pressurization time (tp) adsorption time (ta)

1m 0.404 0.25 m 5.41 × 10-3 m 0.546 1.3 × 10-5 m2/s 716.3 kg/m3 426.7 kg/m3 1.2 × 10-6 J/m/s/K 1046.7 J/kg/K 60 J/m2/s/K (0.136, 0.259) (CH4, H2) 1/s (24124, 8420) (CH4, H2) J/mol 3.73 × 10-8 kg/m/s 8.314 J/mol/K 300 K 310 K (0.7, 0.3) (CH4, H2) 6 × 105 Pa 1.5 × 105 Pa 5s 50 s

discretization of the PDAEs, which results in a huge nonlinear programming problem and can be solved with any state of the art NLP solver. Smith et al.12 considered the cyclic operating schedule of a PSA system as a mixed-integer nonlinear program (MINLP) for optimization. Ko et al.13,14 used an SQP-based approach for optimization of PSA and fractionated vacuum PSA (FVPSA) processes with cyclic steady state. Jiang et al.15 optimized a PSA system by using direct sensitivity and implementing a Newton-based approach to accelerate convergence for cyclic steady states. However, even the most efficient of these approaches can still be time-consuming for large systems. Jiang et al. 16 reported a CPU time of 50-200 h on a 2.4 GHz linux machine for a 5-bed 11-step PSA process optimization to maximize hydrogen recovery. A multiobjective optimization of a simple single bed air drying PSA process by Sankararao et al.17 took 720 h on a 2.99 GHz Pentium 4 machine. This gives a strong motivation to develop cost-efficient and robust optimization strategies for PSA processes. Moreover, in the case of flowsheet optimization, if dynamic PSA models are incorporated with other steady-state models in the flowsheet, then it will require much faster approaches for integrated optimization. Model reduction is a powerful tool that permits the systematic generation of cost-efficient representations of large-scale systems that, in particular, result from discretization of PDEs.18 Over the past decade, powerful model reduction approaches such as proper orthogonal decomposition (POD) have been developed that provide an accurate reduction of large spatially distributed models to much smaller models. POD (also known as Karhunen-Loe´ve approximation) based reduced-order modeling technique is ubiquitous and has been extensively used in obtaining low-dimensional models for efficient simulation19–21 and control22,23 in fluid dynamics. Yuan et al.24 developed reduced-order models (ROMs) for bubbling fluidized beds. Armaou et al.25 applied the concept of model reduction to the time-dependent parabolic PDEs and generated ROMs for diffusion-reaction process. Theodoropoulou et al.26 have extended the applications of ROM-based modeling to chemical vapor deposition. ROMs are derived from solutions of detailed distributed models through the representation of eigenfunction expansions. Using singular value decomposition, spatially distributed eigenfunctions can be derived empirically through a set of well-defined minimization problems. ROMs are then

Table 2. Model Equations of PSA component mass balance

{

-Dx

∂2y1 ∂x

2

(T1 )( ∂x )(∂T∂x ) + 2(P1 )( ∂x )(∂P∂x )} + ∂t + u ∂x + ∂y1

-2

∂y1

∂y1

∂y1

(

RT 1 - εb ∂q1 - y1 F P εb p ∂t

2

∑ ∂t

∂qi

i)1

)

) 0 (1)

y2 ) 1 - y1 1 ) CH4 2 ) H2

(2)

Ergun equation -

( )

2 1.75Fg 1 - εb 2 ∂P 150µ (1 - εb) ) u+ u 2 3 ∂x 4R 2Rp ε εb3 p b

(3)

LDF equation ∂qi ) ki(qi* - qi) i ) 1, 2 ∂t

(4)

energy balance

∑(

)

2

∂T

∂T

(εtFgCpg + FbCps) ∂t + FgCpgεbu ∂x - KL

∂qi ∂2T - Fb ∆Hi + ∂t ∂x2 i)1 2h (T - Tw) ) 0 (5) Rb

Langmur isotherm aiyiP

qi* )

2

1+

∑byP

ai ) R1ieR2iT bi ) β1ieβ2iT i ) 1, 2

(6)

i i

i)1

valve equation u ) uL(x/L) + u0(L - x)/L

uL, u0 )

{

CvPH

CvPH

1 - (PL/PH)2 MT

(7)

if PL/PH > 0.53

1 - (1.893)-2 else MT

cyclic steady state z(t0) ) z(tcycle) z : yi, qi, T ∀ i

(8)

Table 3. Isotherm Parameters for H2 and CH4 on Activated Carbon

methane hydrogen

R1

R2

β1

β2

0.0086 -0.0000379

-0.2155 -0.01815

0.0004066 2.2998

-0.010604 -0.05993

formulated through the substitution of the eigenfunction expansion into the PDE model using Galerkin projection. Truncation of those modes that have no significant contribution to the solution profile then leads to significant reduction in the number of states which eventually leads to a much smaller optimization problem. ROM-based optimization has been used in fluid flow optimal control problems. Fahl27 applied a trust region proper orthogonal decomposition algorithm for the optimal boundary control of a cavity flow. Bergmann et al.28 also applied the same algorithm for the optimal control of the circular cylinder wake flow considered in the laminar regime. Weickum et al.29 developed extended ROMs for optimal design problems, in

Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009 2329 Table 4. Mole Flux Variables at Each Operating Step exhaust flux ) Ej(t) ) (exhaust)j/(available area) outlet flux ) Oj(t) ) (outlet)j/(available area)

feed flux ) Fj(t) ) (feed)j/(available area) purge flux ) Pgj(t) ) (purge)j/(available area) j ) CH4, H2 pressurization

adsorption

depressurization

desorption

Fj(t) ) uPyj/RT|x)0 Oj(t) ) 0 Ej(t) ) 0 Pgj(t) ) 0

Fj(t) ) uPyj/RT|x)0 Oj(t) ) uPyj/RT}|x)L Ej(t) ) 0 Pgj(t) ) 0

Fj(t) ) 0 Oj(t) ) 0 Ej(t) ) -uPyj/RT|x)0 Pgj(t) ) 0

Fj(t) ) 0 Oj(t) ) 0 Ej(t) ) -uPyj/RT|x)0 Pgj(t) ) -uPyj/RT|x)L

Table 5. Boundary Conditions for Each Operating Step pressurization

adsorption

yj|x)0 ) yf, j yj|x)0 ) yf, j (∂yj/∂x)|x)L ) 0 (∂yj/∂x)|x)L ) 0 P|x)0 ) Phigh P|x)0 ) Phigh T|x)0 ) Tfeed T|x)0 ) Tfeed ∂T/∂x|x)L ) 0 ∂T/∂x|x)L ) 0 u|x) 0 ) ufeed u|x)0 ) ufeed u|x)L ) 0.125ufeed u|x)L ) 0

depressurization (∂yj/∂x)|x)0 ) 0 (∂yj/∂x)|x)L ) 0 P|x)0 ) Plow (∂T/∂x)|x)0 ) 0 ∂T/∂x|x)L ) 0 u|x)0 ) -ureg u|x)L ) 0

desorption (∂yj/∂x)|x)0 ) 0 yj|x)L ) yj|x)L of other bed P|x)0 ) Plow (∂T/∂x)|x)0 ) 0 ∂T/∂x|x)L ) 0 u|x)0 ) -ureg u|x)L ) -0.6ureg

which ROMs are developed for the whole design space before optimization and then a trust region based strategy is used for globalization. Bui-Thanh et al.30 provided a goal-oriented framework for generating optimization oriented ROMs by solving a model-constrained optimization problem. However, until now ROM-based optimization has been used only for small-scale optimal control or dynamic optimization problems, and its use for large-scale dynamic optimization problems, especially PSA, has not been explored. In this paper, we develop a POD-based reduced-order model for the large scale optimization of the PSA process. A dynamic PDE-based model of a bench-scale PSA system for separation of hydrogen from a methane-rich mixture is considered and some preliminary ROM-based optimization results for maximizing hydrogen recovery have been generated. The feed composition is 30% hydrogen and 70% methane. The next section presents the process description and the governing equations for the detailed PDE system with the optimization problem. Section 3 describes the proper orthogonal decomposition method for model reduction, while section 4 describes the methodology to generate the reduced-order model and its application to the PSA process. It also describes the solution strategy for solving the detailed PDE model. Section 5 summarizes the ROM-based simulation and optimization results. The concluding remarks and future work have been outlined in the final section. 2. PSA Process and Optimization Problem 2.1. PSA Process. Fixed-bed adsorption systems are typically operated in a cyclic manner with each bed repeatedly undergoing a sequence of steps. In this study, a coupled two-bed system, with beds packed with activated carbon as an adsorbent, has been considered for a PSA cycle operated with four steps: pressurization, adsorption, depressurization (or blowdown), and purge (or desorption), as shown in Figure 1. During the pressurization step for bed 1, high-pressure feed gas consisting of 30% hydrogen and 70% methane at a specified temperature is supplied to the bottom of bed 1, while bed 2 undergoes a depressurization step in which the pressure is reduced to purge pressure and methane is recovered in the exhaust stream. Next, in the adsorption step for bed 1, the high-pressure feed gas is continued at the bottom of the bed 1, where methane gets adsorbed and hydrogen is taken out as a product at the top of the bed. During this period, bed 2 undergoes a desorption step in which a small portion of hydrogen product is drawn as a product purge stream to the bed 2 at low pressure to purge the

accumulated methane in the exhaust stream. In the next two steps the previous two steps are repeated with bed 1 and bed 2 interchanging roles. The target process for this work is benchscale as described in Ko et al.13,14 The design specifications and simulation conditions are listed in Table 1. 2.2. Model Equations. For modeling the above operation, the following assumptions have been taken to avoid complexities in the mathematical model: 1. All of the gases follow the ideal gas law. 2. There are no radial variations in temperature, pressure, and concentrations of the gases in the solid and the gas phase. 3. The gas and the solid phases are in thermal equilibrium and bulk density of the solid phase remains constant. 4. Pressure drop along the bed is calculated by the Ergun equation. 5. The adsorption behaviors are described by the Langmuir isotherm. It has been reported that this monolayer-based isotherm is preferable because of its ability to accurately capture linear behavior at low pressures.2 6. The adsorption rate is approximated by the linear driving force (LDF) expression. Sircar and Hufton31 demonstrated that the LDF model works because the estimation of the separation performance of an adsorptive process requires several sets of averaging of kinetic properties and the effect of local characteristics are lumped during integration. 7. A linear profile has been assumed for superficial gas velocity for all the steps. Cruz et al.32 suggested that this kind of assumption can be taken for PSA processes of bench-scale and there is no need to solve overall mass balance. Valve equations have been used to evaluate velocities on the boundaries. Based on the above assumptions, the mathematical model for the PSA process is listed in Table 2. Here, a lumped mass transfer coefficient has been used for the LDF model. The temperature-dependent adsorption isotherm parameters for hydrogen and methane on activated carbon (R1, R2 and β1, β2) have been taken from Knaebel et al.33 and are shown in Table 3. As suggested by eqs 1 and 2, the PDE for component mass balance has been solved only for methane and the mole fraction of hydrogen has been evaluated by ensuring that the mole fractions sum up to 1. This summation has been enforced because overall mass balance, which implicitly ensures this summation, has not been taken into account for superficial velocity calculation. Since the boundary conditions at the ends of the bed change abruptly and periodically during the operating cycle, PSA processes never reach steady state. Instead, after a startup time, the system approaches cyclic steady state (CSS), where the conditions at the start and end of each cycle in each bed are identical and lead to normal production. The CSS condition (eq 8) shows that the initial conditions are parametric and should be solved simultaneously with the PDE system. We denote the model in Table 2 as the rigorous model for which the ROM is developed. Tables 4 and 5 show the equations for mole flux variables, to calculate the performance (purities and recoveries), and the

2330 Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009

Figure 2. Conditions that need to be matched at the temporal boundaries of the steps, at the beginning and the end of the cycle and at spatial boundaries of both beds during adsorption/desorption.

boundary conditions for each operating step, respectively. Based on the mole flux variables, purities and recoveries of hydrogen and methane are given by

∫O

H2(t)

purityH2 )

dt (9)

n

∫ ∑ O (t) dt j

j)1

∫ E (t) dt ∫ E (t) dt ∫ O (t) dt - ∫ Pg recovery ) ∫ F (t) dt ∫ E (t) dt recovery ) ∫ F (t) dt purityCH4 )

CH4

(10)

n

j)1

j

H2(t)

H2

dt

H2

(11)

H2

CH4

CH4

(12)

the algebraic equations while g represents the design constraints which can be purity, pressure or product rate requirements. tfi is the final time of ith step. Φ represents the objective function and can maximize component recovery or profit at desired purity. The PSA model is governed by nonlinear conservation laws and consists of nonlinear convective-diffusion PDEs. The solution of such systems is characterized by steep profiles in the spatial dimension and hence requires a large number of spatially discretized nodes. This discretization leads to a large set of ordinary differential equations (ODEs) which, after combining for all the steps, results in a large scale optimization problem described by a large number of state variables. In this study, a reduced-order model has been sought for the PSA model described in Table 2, which can be described by a much smaller set of state variables, thus leading to a much smaller optimization problem and reducing the computational complexity. The proper orthogonal decomposition (POD) scheme is used to derive the reduced-order model.

CH4

2.3. PSA Optimization Problem. The PDE-constrained optimization problem for the PSA process can be solved by discretizing them in space, which leads to a set of ODEs. The dynamic optimization problem with CSS conditions, thus obtained after this discretization, can be expressed as: max Φ(y(t), z(t), y0, q) dyi s.t. ) fi(yi(t), zi(t), q) i ) 1, ... , C dt g(y(t), z(t), q) e 0 h(y(t), z(t), q) ) 0 y1(0) ) y0 yi(0) ) yi-1(tfi ) y0 ) yC(tfC)

(13)

i ) 2, ... , C

LB e (y0, q) e UB Here, fi represents the right-hand side of the discretized bed model, obtained after spatial discretization, for each operating step i. C represents the total number of operating steps, which is four for our PSA process. y is the combined vector of the differential state variables, while z is that of algebraic state variables. y0 are the parametric initial conditions for the differential state variables and q are the decision variables. Decision variables can be design parameters such as bed length, diameter or adsorbent packing, or process parameters such as flow rates, step times, operating pressures, and velocities. They are subject to lower (LB) and upper bounds (UB). h represents

3. Proper Orthogonal Decomposition Based Model Reduction Proper orthogonal decomposition (POD) is a model reduction technique which has been used in a wide variety of disciplines such as turbulence, image processing, signal analysis, data compression, oceanography, and process identification and control. The essential idea of POD is to generate an orthonormal subspace of basis functions which are optimal in the sense that they capture the maximum energy of the system and thus allow the state variables of the model to be represented with these (fewer) basis functions. In other words, given an ensemble of observations defined in a vector space, we seek to find a subspace, much smaller than the original vector space, such that the projection of all the observations on to that is maximal. The attractiveness of the POD technique lies in the fact that the basis functions are derived from the numerical solutions or experimental measurements of the system, thus capturing its physics. The method of snapshots has been used in this work to generate POD basis functions.34 POD-based model reduction begins with the collection of snapshot sets which consist of solutions of the PDEs at several time instants during the evolution of the system. These snapshot sets are obtained by solving a rigorous, very large-dimensional system obtained after the spatial discretization (and temporal also in some cases) of the PDEs. The determination of these sets is crucial to the effectiveness of POD-based reduced-order modeling. Hence, they must contain sufficient information to accurately represent

Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009 2331

the dynamics of the system. One then uses the set of snapshots to determine a POD basis set which can accurately capture the information contained in the snapshots using a much smaller set of basis functions. Given the snapshot set YSNAP ) {y1, ... , yNt}

Table 6. ROM for the PSA Process M

Nt

M

i

i)1

j

j 2

j)1

i)1

(yi, φj)φj|2

((

-Dx

j

j

M

d2y0i dx

+

2

∑

ayij

dx

i ∈ {1, ... , Nt}

M

1 - b b

Ergun equation

(∑ M

-

aPj

dφPj

j)1

dx

j

R

d 2 dx TR

dy0i + dx

daqij dt

2

)

, φyik

((

energy balance M

t

∑

FRg CRpg(φTk,φTj)

j)1

∑a

∑∑ r)1 j)1

M

∑ j)1

(

yij

j)1

) ))

dφyij , φyik + dx

TRyiRφqri R

P

) )

, φyik

) )

daqrj + dt

dφyij ayij , φyik ) 0 ∀ k ∈ [1, M] (18) dx

( )

M

∑ (φ

( ( ∑ ) ) ∑( ∑ ) ( ∑ )

dT0 daT + buFRg CRpg + dt dx 2

r)1

j)1

Langmuir isotherm

aiRyiRPR

1+

∑

biRyiRPR

M

aTj

(20)

j)1

dφTj , φTk + dx

M

∆Hr

M

(

* * qij, φqik)aqij - aqik - (q0i, φqik)) ∀ k ∈ [1, M]

j)1

daTk - Fb FbCps dt

* aqik )

(17)

M

M

1.75Fg 1 - b 2 150µ (1 - b)2 (u, φPk) + (u , φPk) (19) 2 3 2Rp 4Rp b b3

(φqrj, φTk)

j)1

d2φTj

d2T0 daqrj -KL + dt dx2

2h R (T - Tw, φTk) ) 0 Rb

∀ k ∈ [1, M] (21)

* , φqik aiR ) R1ieR2iT , biR ) β1ieR2iT

∀ k ∈ [1, M] (22)

aTj

M

j

)

daqik ) ki( dt

(16)

where σi denotes the singular values of YSNAP. The state variables are then expanded in terms of spatial POD basis functions and unknown temporal coefficients.

∑ a (t)φ (x)

T φqij

, φPk )

i)M+1

y(x, t) )

R

P

* Cpgj(t)φCpgj(x)

LDF equation

Nx

σi2

) ( ( )(

, φyik +

j)1

∑a j)1

dy0i dayik + u + dt dx

From a physical point of view, these first M basis functions capture more energy of the snapshot field than any other set of M orthonormal spatial basis functions. Here, |.| denotes the L2 -norm and (,) denotes the Euclidean inner product. The POD basis functions can be computed from the singular value decomposition (SVD) of the snapshot matrix. Let YSNAP ) UDVT denote the SVD of YSNAP; then the M POD basis functions {φ1,..., φM} are given by the first M left singular vectors, i.e., φi ) ui for i ) 1,..., M, where ui denotes the ith column of U. By construction, POD basis functions are orthonormal. Since UD ) YSNAPV, it follows that the basis vectors are linear combinations of the snapshots. Thus, the PODbased reduced-order modeling has advantages over an empirical model reduction because the basis functions capture the physics of the problem. The aforementioned error of the truncated basis representation is given by

∑

2

( )( ∑ (

RFp

j)1

j)1

εPOD(M) )

∑

d2φyij

Fgj(t)φFgj(x)

M

CRpg(x, t) )

aTj(t)φTj(x)

component mass balance

j)M+1

∑a j)1

M

TR(x, t) ) T0(x) +

Nx

i

M

FRg (x, t) )

qij(t)φqij(x)

j)1

where

∑ (y , φ )φ ,

∑a

j)1

(15)

yi )

* * qij(t)φqij(x)

j)1

M

qiR(x, t) ) q0i(x) +

Nx

∑ | y - ∑ (y , φ )φ | ) ∑ | ∑ i

∑a

qi*R(x, t) )

yij

(14)

with the fields y ) y(x,ti), where Nt is the number of snapshots and Nx is the number of spatial discretization nodes, the POD procedure computes an orthonormal set of basis functions {φ1,..., φNx} which minimizes the error in the truncated basis representation of first M e Nx basis functions, i.e., Nt

yij

j)1

i

εPOD(M) )

M

∑ a (t)φ (x)

yiR(x, t) ) y0i(x) +

dx2

, φTj +

)

R

R

j)1

The PDE system is then projected using a Galerkin-type projection on to the subspace spanned by these basis functions to yield a much smaller M-dimensional ODE system for the unknown expansion coefficients. This system of ODEs (or differential-algebraic equations if there are algebraic variables) represents the reduced-order model (ROM) for the rigorous model. The singular values represent the amount of energy captured by each basis function. They are arranged in descending order in D and it is commonly observed that the first few singular values are significantly larger than the subsequent ones, thus representing most of the captured projection of the system. Therefore, in POD-based ROM, the basis functions corresponding to those smaller singular values are dropped, which leads

cyclic steady state

(∑ M

az,k(0)|pres )

j)1

)

az,j(tcycle)|desφz,j|des, φz,k|pres ∀ k ∈ [1, M], z : yi, qi, T

(23)

to a much smaller subspace spanned by very few basis functions. Hence, a significant model reduction is achieved since generally the value of M is much smaller compared to the value of Nx. M can be less than 10 whereas Nx can be on the order of 100s. 4. ROM-Based Dynamic Optimization Problem for PSA As described in the previous section, the first step in the development of the reduced-order model is the collection of


Figure 3. First six POD basis functions of methane mole fraction for adsorption.

Figure 4. Singular values of gas-phase mole fraction of methane and temperature.

the snapshots of the state variables at various time intervals from the rigorous model. A solution strategy for the PDE-based PSA model has been outlined as below. 4.1. Rigorous Model Simulation Strategy. As mentioned in section 2.1, both PSA beds interact over the cycle and this interaction needs to be considered during the simulation, as conditions at the end of the beds need to be matched. Figure 2

shows the spatial domain and step sequence in time for both beds. Both beds undergo the same steps but in a staggered sequence so that a continuous product flow can be maintained. It shows that since the product stream from one bed is used as a purge stream for the other bed, the conditions at the end of the beds need to be matched during the adsorption/desorption step. Moreover, the final conditions of a step serve as the initial


Figure 5. Comparison of methane mole fraction for all the steps after CSS.

conditions for the next step and CSS condition requires that the conditions be matched at the beginning and at the end of the cycle for both beds. The boundaries where conditions need to be matched are shown in bold.

In our work, we use the method of lines to spatially discretize partial differential equations into sets of differential algebraic equations. The solution of the adsorption system PDEs is characterized by steep moving fronts in spatial domain for gas


Figure 6. Comparison of temperature profiles for all the steps after CSS.

concentrations and temperature. At high pressure, steepness of profiles increases even further, thus increasing computational complexity. In order to ensure the accurate approximation of the adsorption models, the selection of an appropriate discretization method is essential. Because first- and second-order finite difference or finite element methods often introduce physically

unrealistic oscillations (dispersion) and unwanted numerical smearing (damping) near steep adsorption fronts, it is important to use high-resolution methods to mitigate this numerical noise.35 In addition, an accurate solution needs to be addressed with conservative methods, which preserve the mass and energy balance in the spatial direction. Therefore, an upwind-based


Figure 7. Comparison of the profile of gas-phase mole fraction of methane for adsorption step obtained after integrating in MATLAB and solving simultaneously in AMPL Table 7. Comparison of Rigorous Model and ROM Based on the Performance Variables performance variables

rigorous model

ROM

H2 purity H2 recovery CH4 purity CH4 recovery

0.9987 0.1095 0.9421 0.2091

0.9987 0.1094 0.9425 0.2094

Table 8. Optimization Variables and their Values at which ROM was built variable

value

high operating bed pressure (PH) low operating bed pressure (PL) pressurization step time (tp) adsorption step time (ta) feed velocity (ufeed) regeneration velocity (ureg)

5 × 105 Pa 1.5 × 105 Pa 5s 50 s 0.1 m/s -0.05 m/s

Table 9. Optimization Results for Case I problem size and computational time no. of variables no. of constraints no. of function evaluations total no. of iterations total CPU seconds optimal parameters high operating bed pressure (PH) low operating bed pressure (PL) pressurization step time (tp) adsorption step time (ta) performance variables from ROM optimization H2 purity H2 recovery CH4 purity CH4 recovery performance variables from rigorous model simulation at optimum H2 purity H2 recovery CH4 purity CH4 recovery

42760 42756 16 15 195.44 5.2 × 105 Pa 1.3 × 105 Pa 3s 53 s 0.9988 0.1628 0.9541 0.1771 0.9991 0.1629 0.9491 0.1769

finite volume method has been applied since it is particularly suited for modeling conservative properties.36 Here, we apply the van Leer flux limiter with a smoothing function near the origin as the interpolation scheme; this is required for spatially discretized nodes at boundaries of each finite volume, and the scheme is applied for both mole fraction and temperature. This approach was also used in Jiang et al.,15 to avoid oscillations and ensure that solution profiles remains positive.

A successive substitution method has been used to determine CSS, in which the PSA cycle is started with a specified initial condition and a series of cycles are simulated after that, with the initial condition of each cycle taken from the final condition of the previous cycle, until the bed conditions do not change from cycle to cycle. Since a bench-scale PSA model is considered in this work,13 cyclic steady state is achieved after 100-120 cycles. 4.2. ROM Formulation Strategy. The snapshots obtained from the rigorous model simulation were used to generate the POD basis functions and the model equations of PSA were then projected on to the subspace spanned by them, which eventually gave a set of differential algebraic equations (DAEs). Table 6 shows the DAE system obtained after the Galerkin projection of the model, shown in Table 2, on to the basis functions. The state variables are expressed as a summation of the spatial basis functions and substituted in the PDEs. The inner product with each basis function then results in the DAE system. In this way, a square system of DAEs is obtained. The boundary conditions were handled very differently in the case of the reduced-order model. The snapshots taken for the POD basis functions comprised only the interior data of the bed, and the boundary data are not considered here. Hence the ROM was developed and solved only for the interior of the bed. Unlike rigorous model simulation, the boundary values were evaluated after the interior had been obtained. For Dirichlet boundary conditions, the boundary values were retained as they were in the rigorous model. For Neumann boundary conditions, the boundary values were evaluated using the interior values obtained after solving the reduced-order model for the unknown coefficients. Since the snapshots were obtained after the rigorous model reached cyclic steady state, the snapshots for both beds for the corresponding steps were identical, i.e., the snapshots of bed 1 for the pressurization step are identical to the snapshots of the pressurization step for bed 2, and so on for all the steps. Therefore, the reduced-order model could be built for only one bed and the other bed was ignored in the ROM. Coupling for the adsorption/desorption steps of the two beds was ensured by the adsorption and desorption steps of the same bed. Hence, a greater model reduction was achieved since the DAEs of the reduced-order model were simulated only for one bed containing expansions for the state variables of only that bed, while the rigorous model DAEs were simulated for both beds simulta-


Figure 8. Comparison of mole fraction of methane after optimizing case I.

neously. The number of DAEs for both the rigorous model and the ROM is compared in section 5.1. 4.2. ROM-Based Optimization. The reduced-order model can be used in place of the spatially discretized rigorous model for the optimization problem (13). Since the set of DAEs in

the ROM is much smaller compared to the set of DAEs obtained after spatial discretization of the rigorous model, this results in a much smaller optimization problem. One of the advantages of ROM-based optimization is the ability to use complete spatial and temporal discretization for

Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009 2337 Table 10. Perturbation Results perturbed variable

PH

PL

tp

ta

optimal value perturbed value H2 recovery (after perturbation)a

5.2 × 105 5.3 × 105 0.1675

1.3 × 105 1.25 × 105 0.1723

3 2.5 0.1643

53 55 0.1632

a

Optimal hydrogen recovery ) 0.1628.

optimization which is much more difficult with rigorous modelbased optimization, as many more discretization nodes are required to capture the complex dynamics of the PSA process governed by steep fronts. Therefore, a complete discretization for all the operating steps leads to a huge set of algebraic equations, which becomes prohibitively time-consuming and memory intensive for the optimization, as reported by Jiang et al.15 and Ko et al.13 On the other hand, as a much smaller set of DAEs is considered in the ROM-based optimization problem, temporal discretization leads to a much smaller nonlinear programming problem (NLP). This is a significant advantage since efficient, state-of-the-art NLP solvers, such as IPOPT,37 can now be used for ROM-based optimization. The ROM-based optimization problem based on (13) is outlined below. Here the trapezoidal rule has been used in this work to discretize ODEs in the temporal dimension. max Φ(a, a0, q) k + 0.5h(Fi(ak-1 s.t. aki ) ak-1 i i , q) + Fi(ai , q)) i ) 1, ... , C, k ) 0, ... , Ti g(a, q) e 0 h(a, q) ) 0 a01 ) a0 Ti-1 a0i ) (ai-1 φi-1, φi)

i ) 2, ... , C

a0 ) (aTCCφC, φ1) LB e (a0, q) e UB (24) Here, Fi represents the ROM counterpart of the corresponding rigorous model function fi described by eq 13. a is the combined vector of the unknown coefficients for all the state variables for all the steps C and φ are the POD basis functions. Ti represents the number of temporal discretization nodes for each operating step i. a0 is the reduced set of initial conditions for the unknown coefficients in the expansion for the first operating step. The initial condition and cyclic steady state condition changes and a reduced set of equations based on the inner product with POD basis functions is used for reduced set of state variables a0i . The inequalities, algebraic equations, and the decision variables in this optimization problem remain the same and depend upon the state variables a in the ROM. A reduced-order model is accurate for the values of the parameters at which it is constructed. However, it loses its accuracy at a different point in the decision variable space since the snapshots do not capture the spatial behavior and dynamics of the system at that new point. The error in the solution given by the ROM increases as we go further away from the point where it was constructed. Hence, tighter variable bounds are needed in the ROM-based optimization problem, so that the algorithm remains close to the point where ROM was developed. This property suggests that an adaptive scheme is required for ROM-based optimization, so that the ROM will be repeatedly updated based on the error from the previous ROM optimization. This adaptive scheme is beyond the scope of this work, where we focus on optimization problems within a neighborhood of the constructed ROM. However, we note in the last section that

this optimization scheme forms a key step in an adaptive optimization that will be developed in the future. 5. Case Study and Computational Results 5.1. Comparison of Rigorous Model and ROM at CSS. As mentioned in section 4.1, the PDE-based PSA model was discretized in space to obtain a set of DAEs which were then integrated. For our study, 35 spatial nodes were considered for discretization of the PDEs and the resulting ODE system with algebraic equations was integrated using BDF method (ode15s in MATLAB38). Cyclic steady state was achieved up to desired tolerance after simulating the model repeatedly for 120 cycles. Snapshots were taken after CSS was reached and separate POD basis functions were derived for gas phase mole fractions, solid state loadings, temperature, gas density, specific heat and equilibrium solid concentrations, respectively. Moreover, separate POD basis functions were generated for pressurization, adsorption, depressurization and desorption steps. Figure 3 shows POD basis functions of gas-phase methane mole fraction for adsorption step. Since these functions are derived from data snapshots, their shapes are different for all four operating steps. Figure 4 shows singular values for gas-phase mole fraction of methane and temperature for all the operating steps. The slope of the curve shows that the singular values decrease sharply and the first 5-6 values account for around 99.9% of the energy of the system. Therefore, all the state variables for all the operating steps could be represented using only 5-6 spatial POD basis functions compared to 35 spatially discretized nodes. With 35 spatial discretization points, a total of 2800 DAEs were obtained after applying method of lines on the rigorous model. This includes equations for all the four operating steps and for both beds. However, using only 5 basis functions, the reducedorder model obtained after Galerkin projection comprised only 200 DAEs which is 1/14th of the rigorous system. As mentioned in section 4.2, the ROM was developed only for one bed which led to significant model reduction. The ODEs of the reduced-order model for all the operating steps were discretized in time, giving a large set of algebraic equations, and solved simultaneously in AMPL39 using IPOPT. Ti ) 30 nodes were used for temporal discretization. Instead of solving this model repeatedly for CSS, the initial conditions were taken as variables and reduced CSS conditions (shown in Table 6) were solved simultaneously with the model equations in AMPL. Figure 5 compares the profiles of gas-phase methane mole fraction obtained from rigorous model simulation and the ROM model with 5 POD basis functions. It can be seen that 5-basis approximation captures the dynamics of the problem very well and the profiles are nearly identical. Similarly, the temperature profile in Figure 6 shows the accuracy of the approximation for all the steps. The profiles were plotted after CSS was achieved. It can be observed that the profiles of gasphase mole fraction of methane are steep in the spatial dimension, especially for adsorption and depressurization steps, and the steepness has been effectively captured by the reducedorder model. Moreover, the behavior of the model in the temporal dimension has also been adequately captured by the ROM. To check for the accuracy of complete discretization of ROM in AMPL, the DAE system of the ROM was integrated in MATLAB and solved in AMPL without the CSS conditions. Figure 7 compares the profile of the gas-phase mole fraction of methane for the adsorption step for both approaches. The profiles compare very well and show that the additional errors introduced by a predetermined temporal discretization in AMPL were


Figure 9. Oscillatory profiles of gas-phase methane mole fraction due to relaxed bounds. Table 11. Comparison of Hydrogen Performance for Relaxed Optimization performance variable

ROM

rigorous model

H2 recovery H2 purity

0.3482 0.9981

0.2763 0.8032

Table 12. Optimization Results for Case II problem size and computational time no. of variables no. of constraints no. of function evaluations total no. of iterations total CPU seconds optimal parameters high operating bed pressure (PH) low operating bed pressure (PL) pressurization step time (tp) adsorption step time (ta) feed velocity (ufeed) regeneration velocity (ureg) performance variables from ROM optimization H2 purity H2 recovery CH4 purity CH4 recovery performance variables from rigorous model simulation at optimum H2 purity H2 recovery CH4 purity CH4 recovery

43040 43034 20 13 184.493 5.1 × 105 Pa 1.4 × 105 Pa 4s 51 s 0.11 m/s -0.0495 m/s 0.9986 0.1511 0.9507 0.1827 0.9974 0.1509 0.9495 0.1826

negligible for this system of equations. Hence, as mentioned in section 4.3, the complete discretization approach is not only accurate for this set of equations but also allows simultaneous optimization with cyclic steady state conditions using a stateof-the-art NLP solver IPOPT. Table 7 compares the purity and recovery of hydrogen and methane obtained from the rigorous model and the reducedorder model at CSS. This can be seen as another basis to

compare the accuracy of ROM. The values were obtained from a 5-basis approximation and are accurate up to 3 decimal places. 5.2. Optimization with ROM. A general ROM-based optimization problem has been discussed in section 4.3. In this work, the ROM has been optimized to maximize recovery of hydrogen subject to the DAEs of ROM, CSS conditions and constraints on the operating variables. Table 8 shows the decision variables considered for this optimization problem and their values at which ROM was built. As mentioned before, the ROM is accurate at the parameter values at which it is constructed and the error in the solution given by ROM increases as we go further away from that point. Hence, tight bounds have been used for the parameters in the ROM-based optimization problem. Two separate cases for optimization were considered. In the first case, feed and regeneration velocities were not taken as the decision variables, although they were included as decision variables in the second case. The optmization results of both cases are discussed in the subsequent sections. 5.2.1. Case I: Optimization without Velocities. The reducedorder model was optimized in AMPL using IPOPT with operating pressures and step times as the decision variables. The feed and regeneration velocity were held constant in this case. Cyclic steady state conditions were included in the optimization problem. The optimization problem can be described as below: max RecoveryH2 s.t. eqs 18 - 23 4.8 × 105 e PH e 5.2 × 105 1.3 × 105 e PL e 1.7 × 105 3 e tp e 7 47 e ta e 53 purityH2 g 0.998

(25)


Figure 10. Comparison of mole fraction of methane after optimizing case II.

The optimal values of the parameters along with computational time are included in Table 9. The recovery of hydrogen could be increased up to 16.3% despite the tight bounds on decision variables. The optimum was achieved in only 195 CPU seconds since a much smaller model was used for optimization.

Moreover, the optimization was cheap since fewer function evaluations were needed to achieve the optimum. To check the accuracy of the optimization results, the rigorous model was simulated at the optimal values of the decision variables at cyclic steady state. Since the ROM is developed at


Figure 11. Profiles of mole fraction of methane after ROM optimization.

some other point and did not contain any information about this new point, it is bound to have some errors. Therefore, we analyzed this error by simulating the original rigorous model at this new optimum point. The purities and recoveries obtained from rigorous model simulation have also been listed in Table 9. These values are very close to the ones obtained from ROM optimization and it can be inferred that ROM is accurate at this optimum point too. The increase in hydrogen recovery is almost the same in both cases, which indicates that the error is not large at the optimum point in the ROM. Moreover, the hydrogen purity constraint is also satisfied by the rigorous model simulation results. To further consolidate this observation, gasphase mole fraction profiles of methane have been plotted in Figure 8 for all the steps. The methane mole fraction profiles were generated from rigorous model simulations at the optimum point and compared with profiles obtained from ROM after optimization. The profiles match reasonably accurately and show that the behavior of the ROM is close to the rigorous model at this optimum point. It is worth noting that all the decision variables were at their bounds after optimization. Therefore, to verify that (locally) optimal values for ROM are also optimal for the rigorous model, we evaluate the KKT conditions of the rigorous model at this point, x*. A general rigorous model based optimization problem with its KKT conditions is given by max f(x) s.t.

∇f(x*) + ∇ c(x)λ + µl - µu ) 0

c(x) ) 0

c(x*) ) 0,

aexeb

µl(x* - a) ) 0

µl, µu g 0

µu(x* - b) ) 0 Shifting the decision variables to the beginning of vector x and defining a null space basis matrix that satisfies ZT[I 0]T ) I and ZT∇c(x*) ) 0 allows the KKT conditions to be written as

ZT ∇ f(x*) - µu ) 0 : for decision variables at upper bound (µl ) 0) ZT ∇ f(x*) + µl ) 0 : for decision variables at lower bound (µu ) 0) Since µu, µl g 0, the KKT conditions simplify to (ZT∇f(x*))i e 0 if decision variable i is at its lower bound, and (ZT∇f(x*))i g 0 if decision variable i is at its upper bound. These properties were indeed verified by perturbing the bounded decision variables and evaluating the rigorous model for the perturbed decision variables. The results are shown in Table 10. It can be observed from the profiles obtained after ROM optimization that the error in the reduced-order model remains small at the optimum. This may be due to the tight bounds on the decision variables in optimization problem (25). An optimization was also solved with relaxed bounds on the decision variables as shown below: 3 × 105 e PH e 13 × 105 1 × 105 e PL e 2.5 × 105 0.5 e tp e 10 30 e ta e 80

(26)

Although, with the given hydrogen purity constraint, hydrogen recovery could be increased up to 34.8% in this case, the profiles were oscillatory and unrealistic. Figure 9 presents the methane gas-phase mole fraction profiles for optimization problem (26). For adsorption, the oscillations are quite big and they tend to increase as step time increases. In the case of depressurization, there is a jump in the profile before steep decrease along spatial dimension. Large oscillations in the profiles thus show large error in the reduced-order model at the optimum. Moreover, Table 11 shows that when the rigorous model was simulated to CSS at the optimal values then hydrogen purity dipped to 80%, compared to 99.8% given by ROM optimization. Hence, it can


be inferred that tighter restrictions are required on decision variables for ROM-based optimization. Recovery of hydrogen can be further improved by following this ROM-based optimization approach. For further optimization from this new optimal point, we can either retain the same ROM if the error does not seem too large or construct a new one by taking new snapshots from the rigorous model simulation. This process can be repeated and we can use ROM for optimization instead of the detailed rigorous model. An adaptive trust region framework can be used for ROM-based optimization.27,40,41 Here, the ROM updating strategy can be determined based on the ratio of actual to predicted improvement in the objective and the trust region radius can be changed based on the update rules. 5.2.2. Case II: Optimization with Velocities. In this case, besides operating pressures and step times, feed and regeneration velocities were also taken as optimization decision variables and recovery of hydrogen was maximized with cyclic steady state conditions. The optimization problem subject to tight bounds is described as below:

1.3 × 105 e PL e 1.7 × 105 3 e tp e 7 47 e ta e 53 0.08 e ufeed e 0.12 -0.053 e ureg e -0.047

(28)

Compared to the previous case, the hydrogen recovery can now be increased to 21.6%. However, this extra increment in the recovery was marred by oscillatory profiles for state variables obtained after optimization. Figure 11 shows the profiles of gasphase mole fraction of methane. It can be observed that profiles are oscillatory for adsorption and depressurization steps, similar to the previous case of without velocities. The profiles are physically unrealistic and ROM does not adequately describe the behavior of the process at the optimum. Hence, this case also suggests that tight restrictions are required on the decision variables and an adaptive trust region based strategy with appropriateROMupdationisrequiredforROM-basedoptimization. 6. Conclusions and Future work

max recoveryH2 s.t. eqs 18-23 4.9 × 105 e PH e 5.1 × 105 1.4 × 105 e PL e 1.6 × 105 4 e tp e 6 49 e ta e 51 0.09 e ufeed e 0.11 -0.0505 e ureg e -0.0495 purityH2 g 0.998

4.8 × 105 e PH e 5.2 × 105

(27)

The model was optimized in AMPL using IPOPT. The optimization results along with computational time have been shown in Table 12. The recovery of hydrogen could be increased up to 15.11% despite tight bounds on decision variables. As in the previous case, a fast optimization could be perfomed and the optimum was achieved in only 184 CPU seconds. However, this case was marginally more expensive as slightly more function evaluations were needed. As done in the previous case, the rigorous model was simulated at the optimal values to check the accuracy of the optimization results. The purities and recoveries obtained from rigorous model simulation have also been listed in Table 12. These values are very close to the ones obtained from ROM optimization, indicating that the error is not large at the optimum point in ROM. Moreover, the gasphase methane mole fraction profiles, obtained from rigorous model simulations and ROM optimization, have been compared in Figure 10. Barring small oscillations in adsorption profile, the profiles match reasonably accurately, thus showing that the behavior of the ROM is close to the rigorous model at the optimum point. Unlike the previous case without velocities, in this case the hydrogen purity constraint is violated slightly by the rigorous model simulation. The rigorous model gives a hydrogen purity of 99.74% which is slighly less than the specified lower bound 99.8%. Since the decision variables are at their bounds after optimization, the rigorous model has also achieved optimum in this case. On the other hand, the slight violation can be made acceptable by imposing a tighter purity constraint in the ROM optimization. Since tight bounds are used in the aforementioned optimization problem, the error in the ROM remains small. Instead, an optimization was also performed with slightly relaxed bounds, as shown below:

This work presents the development of a reduced-order model for the two-bed four-step PSA process for a hydrogen-methane system using proper orthogonal decomposition scheme. The method of snapshots was used to generate POD basis functions, and the PDEs of the process were projected to generate the ROM. It was shown that a significant model reduction was achieved using the POD scheme and 4-5 basis functions sufficently describe the dynamic model for all the operating steps. The reduced-order model based on 5-basis approximations captures the dynamics of the process very accurately. The profiles for the state variables derived from both the rigorous original model and the reduced-order model simulations were shown to be nearly identical. Moreover, both reduced-order model and rigorous model predicted approximately the same values for purity and recovery of hydrogen and methane. Thus, ROM was successfully able to mimic the behavior of the PSA system. With ROM-based optimization, we were able to use complete spatial and temporal discretization for optimization, which is very difficult to achieve with rigorous model based optimization, since it leads to a prohibitively large NLP. Since ROM-based problem consists of a small set of DAEs, it does not lead to a large NLP. The ROM was used to maximize recovery of hydrogen with operating pressures and step times as the decision variables. Tight bounds were used for parameters due to limited validity of the ROM in the vicinity of the point of construction. Extremely fast and cheap optimization could be performed with ROM and the optimum was achieved in approximately three CPU minutes. Moreover, the optimum attained match the optimum achieved with rigorous model optimization and the profiles of the state variables obtained after ROM-based optimization precisely matched the profiles of the rigorous model simulation at the optimum. In another case, feed and regeneration velocities were also included with other parameters and equivalently fast and cheap optimization was performed to maximize hydrogen recovery with tight bounds on parameters. In this case also the optimum achieved was an optimum for the rigorous model and the profiles obtained after ROM-based optimization matched the ones from rigorous model simulation. An optimization with relaxed bounds on parameters, for both cases, yielded oscillatory and physically inaccurate profiles for state variables, thus suggesting a need for tight restrictions on parameters for ROM-based optimization. Current results indicate


that a proposed ROM methodology is a promising surrogate modeling technique for cost-effective optimization purposes, but also that an adaptive trust region strategy is needed to determine when the ROM needs to be updated within the optimization cycle. In future work a rigorous trust-region based framework will be designed in which the problem will be constrained by the DAEs of the reduced-order model and optimization will be performed in a trust region such that the decision variables are bounded. The trust region update rules and sufficient decrease condition for the objective will be used to determine the size of the trust-region and based on the decrease and error in the ROM a strategy will be designed to decide whether to update the ROM or not.27,41 To ensure that the first-order consistency condition is met and the optimum obtained from ROM-based optimization corresponds to the optimum of the original problem, a scaling function, such as the one proposed by Alexandrov et al.,40 can be incorporated in the objective function. Numerical inconsistencies, such as oscillations, can also be treated with this error control mechanism. Finally, while this work considered a bench-scale PSA process, future work will deal with large industrial-scale systems addressed with our proposed ROM-based optimization. Acknowledgment This technical effort was performed in support of the National Energy Technology Laboratory’s ongoing research in Process and Dynamic Systems Research under the RDS contract DEAC26-04NT41817. Nomenclature Cpg ) average heat capacity of gas (J/kg/K) Cps ) heat capacity of adsorbent (J/kg/K) Dx ) gas diffusivity (m2/s) ∆H ) isosteric heat of adsorption (J/mol) h ) heat transfer coefficient (J/m2/s/K) KL ) effective axial thermal conductivity (J/m/s/K) kj ) lumped mass transfer coefficient for component j (1/s) L ) bed length (m) P ) total pressure (Pa) PH ) high operating pressure within bed (same as Phigh) (Pa) PL ) low operating pressure within bed (same as Plow) (Pa) Pfeed ) feed pressure (Pa) Ppurge ) exhaust pressure (Pa) qj ) solid-phase loading for component j (mol/kg) qj* ) equilibrium solid loading for component j (mol/kg) Rb ) bed radius (m) Rp ) particle radius (m) T ) gas temperature within the bed (K) Tfeed ) feed temperature (K) Tw ) column wall temperature(K) ta ) adsorption time (same as desorption time) (s) tp ) pressurization time (same as depressurization time) (s) u ) gas superficial velocity (m/s) ufeed ) feed velocity (m/s) ureg ) exhaust velocity (m/s) yj ) gas-phase mole fraction for component j Greek Letters Rkj ) temperature-dependent Langmuir constant for component j (k ) 1, 2) (mol/kg/Pa) βkj ) temperature-dependent Langmuir constantfor component j (k ) 1, 2) (1/Pa) b ) bed void fraction

p ) particle void fraction t ) total void fraction µ ) gas viscosity (kg/m/s) Fb ) bed density (kg/m 3) Fg ) gas density (kg/m 3) Fp ) particle density (kg/m 3) φ ) POD basis functions Subscripts des ) desorption pres ) pressurization i, k ) ith or kth POD basis function j ) jth component (1 ) CH4, 2 ) H2)

Literature Cited (1) Keller, G. E. Gas Adsorption Processes: State of the Art. In Industrial Gas Separations; Whyte, T. E. Ed.; American Chemical Society: Washington, DC, 1983; ACS Symposium Series 223, Vol. 145. (2) Ruthven, D. M.; Farooq, S.; Knaebel, K. S. Pressure Swing Adsorption; VCH Publishers: New York, 1994. (3) Yang, R. T. Gas Separation by Adsorption Processes; Butterworth: Boston, MA, 1987. (4) Miller, G. Q.; Sto¨cker, J. Selection of a Hydrogen Separation Process; NPRA Annual Meeting, San Francisco, CA, 1989. (5) Peramanu, S.; Cox, B. G.; Pruden, B. B. Economics of hydrogen recovery processes for the purification of hydroprocessor purge and offgases. Int. J. Hydrogen Energy 1999, 24 (5), 405. (6) Voss, C. Applications of Pressure Swing Adsorption Technology. Adsorption 2005, 11, 527. (7) Sircar, S.; Golden, T. C. Purification of Hydrogen by Pressure Swing Adsorption. Sep. Sci. Technol. 2000, 35 (5), 667. (8) Cen, P.; Yang, R. T. Separation of a Five-Component Gas Mixture by Pressure Swing Adsorption. Sep. Sci. Technol. 1985, 20 (9, 10), 725. (9) Yang, J.; Lee, C. Adsorption Dynamics of a Layered Bed PSA for H2 Recovery from Coke Oven Gas. AIChE J. 1998, 44 (6), 1325. (10) Sircar, S. Pressure Swing Adsorption. Commentaries. Ind. Eng. Chem. Res. 2002, 41, 1389. (11) Nilchan, S.; Pantelides, C. C. On the Optimisation of Periodic Adsorption Processes. Adsorption 1998, 4, 113. (12) Smith, O. J.; Westerberg, A. W. The Optimal Design of Pressure Swing Adsorption Systems. Chem. Eng. Sci. 1991, 46, 2967. (13) Ko, D.; Siriwardane, R.; Biegler, L. T. Optimization of Pressure Swing Adsorption Process using Zeolite 13X for CO 2 Sequestration. Ind. Eng. Chem. Res. 2003, 42, 339. (14) Ko, D.; Siriwardane, R.; Biegler, L. T. Optimization of Pressure Swing Adsorption and Fractionated Vacuum Pressure Swing Adsorption Processes for CO2 Capture. Ind. Eng. Chem. Res. 2005, 44, 8084. (15) Jiang, L.; Biegler, L. T. Simulation and Optimization of PressureSwing Adsorption Systems for Air Separation. AIChE J. 2003, 49 (5), 1140. (16) Jiang, L.; Fox, V. G.; Biegler, L. T. Simulation and Optimal Design of Multiple-Bed Pressure-Swing Adsorption Systems. AIChE J. 2004, 50 (11), 2904. (17) Sankararao, B.; Gupta, S. K. Multi-Objective Optimization of Pressure Swing Adsorbers for Air Separation. Ind. Eng. Chem. Res. 2001, 46 (11), 3751. (18) Sorensen, D. C.; Antoulas, A. C. Approximation of large-scale dynamical systems: An overview. Int. J. Appl. Math. Comput. Sci. 2001, 11, 1093. (19) Berkooz, G.; Holmes, P.; Lumley, J. L. The Proper Orthogonal Decomposition in the Analysis of Turbulent Flows. Annu. ReV. Fluid Mech. 1993, 25, 539. (20) Kunisch, K.; Volkwein, S. Galerkin Proper Orthogonal Decomposition Methods for a General Equation in Fluid Dynamics. SIAM J. Numer. Anal. 2002, 40 (2), 492. (21) Rowley, C. W.; Colonius, T.; Murray, R. M. Model reduction for compressible flows using POD and Galerkin projection. Physica D 2004, 189, 115. (22) Kunisch, K.; Volkwein, S. Control of the Burgers Equation by a Reduced-Order Approach Using Proper Orthogonal Decomposition. J. Opt. Theory Appl. 1999, 102 (2), 345. (23) Ravindran, S. S. Proper Orthogonal Decomposition in Optimal Control of Fluids; NASA/TM-1999-209113, 1999. (24) Yuan, T.; Cizmas, P. G.; O’Brien, T. A reduced-order model for a bubbling fluidized bed based on proper orthogonal decomposition. Comput. Chem. Eng. 2005, 30, 243.

Ind. Eng. Chem. Res., Vol. 48, No. 5, 2009 2343 (25) Armaou, A.; Christofides, P. D. Dynamic optimization of dissipative PDE systems using nonlinear order reduction. Chem. Eng. Sci. 2002, 7, 5083. (26) Theodoropoulou, A.; Adomaitis, R. A.; Zafiriou, E. Model Reduction for Optimization of Rapid Thermal Chemical Vapor Deposition Systems. IEEE Trans. Semicond. Manuf. 1998, 11, 85. (27) Fahl, M. “Trust-region Methods for Flow Control based on Reduced Order Modelling”, Doctoral Dissertation, Trier University, 2000. (28) Bergmann, M.; Cordier, L.; Brancher, J.-P. Optimal rotary control of the cylinder wake using proper orthogonal decomposition reduced-order model. Phys. Fluids 2005, 17, 097101. (29) Weickum, G.; Eldred, M. S.; Maute, K. Multi-point Extended Reduced Order Modeling For Design Optimization and Uncertainty Analysis;Proc. 47th AIAA/ASME/ASCE/AHS/ASC Struct., Struct. Dyn., Mater. Conf., Newport, RI, 2006; AIAA Paper 2006-2145. (30) Bui-Thanh, T.; Willcox, K.; Ghattas, O.; Waanders, B. v. B. Goaloriented, model-constrained optimization for reduction of large-scale systems. J. Comput. Phys. 2007, 224 (2), 880. (31) Sircar, S.; Hufton, J. R. Why Does the Linear Driving Force Model for Adsorption Kinetics Work? Adsorption 2000, 6, 137. (32) Cruz, P.; Santos, J. C.; Magalhaes, F. D.; Mendes, A. Simulation of separation processes using finite volume method. Comput. Chem. Eng. 2005, 30, 83. (33) Knaebel, S. P.; Ko, D.; Biegler, L. T. Simulation and Optimization of a Pressure Swing Adsorption System: Recovering Hydrogen from Methane. Adsorption 2005, 11, 615.

(34) Sirovich, L. Turbulence and the dynamics of coherent structures. Q. Appl. Math 1987, 45 (3), 561. (35) LeVeque, R. J. Numerical Methods for ConserVatiVe Laws; Birkha¨user: Basel, Switzerland, 1992. (36) 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. (37) Wa¨chter, A.; Biegler, L. T. On the Implementation of an Interior Point Filter Line Search Algorithm for Large-Scale Nonlinear Programming. Math. Prog. 2006, 106 (1), 25-57. (38) MATLAB User’s Guide, The MathWorks, Inc., 1994-2005. (39) Fourer, R.; Gay, D. M.; Brian, W.; Kernighan, A. Modeling Language for Mathematical Programming. Manage. Sci. 1990, 36, 519. (40) Alexandrov, N. M.; Dennis, J. E. Jr.; Lewis, R. M.; Torczon, V. A trust-region framework for managing the use of approximation models in optimization. Struct. Multidisc. Optim. 1998, 15, 16. (41) Toint, P. L. Global convergence of a class of trust-region methods for nonconvex minimization in Hilbert space. IMA J. Numer. Anal. 1988, 8 (2), 231.

ReceiVed for reView October 18, 2007 ReVised manuscript receiVed January 24, 2008 Accepted January 28, 2008 IE071416P

Simulation and Optimization of Pressure Swing Adsorption Systems

Recommend Documents