Predictive Dynamic Model of a Small Pressure Swing Adsorption Air

Ind. Eng. Chem. Res. 1999, 38, 3761-3775

3761

Predictive Dynamic Model of a Small Pressure Swing Adsorption Air Separation Unit Kenneth G. Teague, Jr. and Thomas F. Edgar* Department of Chemical Engineering, University of Texas at Austin, Austin, Texas 78712

A predictive dynamic model of a small pressure swing adsorption (PSA) air separation process was developed for the purposes of evaluation, optimization, and control of oxygen generation systems on board military aircraft. A mathematical model of the adsorption beds was formulated by application of fundamental mass- and energy-transport modeling techniques. These equations were discretized using the Galerkin finite element technique. The resulting ODE systems were coupled with ODEs describing the rate of change of pressure in each bed and models of the feed and exhaust valves and purge orifice. The model was developed so that it is possible to predict the dynamic response of product oxygen composition and feed air consumption to step changes in feed pressure, product flow rate, and cycle time. A laboratory PSA unit similar in size to an on-board oxygen generation system (OBOGS) was constructed to validate the model. The laboratory unit was constructed so that step changes could be implemented and the responses observed for comparison with the model. All parameters in the model were estimated from literature sources with the exception of the feed/exhaust valve and purge orifice discharge coefficients. Excellent dynamic predictions of bed pressure, cycle-averaged feed flow rate, and cycle-averaged bed temperature vs time in response to step changes in all three input variables compared to the two-bed PSA data were achieved without additional parameter estimation from two-bed data. Introduction A small-scale application of pressure swing adsorption (PSA) is employed to provide oxygen-enriched air for the crew of a military aircraft. Such a unit is known as an on-board oxygen generation system (OBOGS). OBOGS is preferred over bottled oxygen because it is safer and it does not place a limitation on flight time as bottled oxygen does. The application of PSA to air separation on a military aircraft presents a need for a dynamic model of the process which can predict the response of product composition and feed consumption to step changes in feed pressure, feed composition, cycle time, and product flow rate. Ultimately, the military needs a tool that can be used to predict how a given contaminant or chemical warfare agent will affect the crew if a contaminant enters the feed air of the OBOGS. A dynamic model could also be used to evaluate and optimize prospective oxygen generation units. This process involves assessing the performance of the unit over a range of operating conditions. Presently, this is done experimentally by developing a large experimental "matrix" that relates process performance (oxygen purity and recovery) to operating variables (product flow, feed pressure, cycle time, and temperature). This information can be used to locate the optimal operating conditions or the optimal purge orifice characteristics for a given breathing flow requirement, and it can be used to respond to questions that may arise after the unit is placed into service. A third need arises from efforts to develop model-based control of on-board PSA units. The above three challenges require that a computer model of a small PSA air separation process be devel* To whom correspondence should be addressed. Phone: (512) 471-3080. E-mail:[email protected].

oped. Specifically, for a given set of process characteristics (bed dimensions and valve and purge orifice specifications), the computer model must (1) correlate the initial cyclic steady-state performance (product oxygen composition and oxygen recovery), given the product flow rate, cycle time, feed pressure (1-5 bar), and (2) predict the dynamic response of the unit from the initial cyclic steady state to changes in the same variables. The PSA cycle has been examined extensively in the published literature since Shendalman and Mitchell1 modeled the removal of trace CO2 from helium by a PSA unit. Much of this work has been reviewed more recently.2 A limited emphasis, however, has been placed on computing the transient response of process variables between cyclic steady states, in response to changes in the input variables (e.g., feed pressure, cycle time, and product draw rate). Such predictions require that the bed pressure (and dynamic response of the bed pressure) be computed as well. A few investigations have addressed either the dynamic response or computation of the bed pressure history. Ruthven has noted that the solution for the bed pressure history is problematic.2 It is more convenient to make assumptions about the bed pressure history. Beaman3 modeled a PSA unit as a complete system, accounting for the beds, valves, and a product plenum. His model of air separation is isothermal and handles the two major components (oxygen and nitrogen). This model allows specification of known input variables (feed and exhaust pressure and product flow). However, the solution is semianalytic and thus limited to the assumptions made. Experimental validation of bed pressure history predictions also was not considered. Chou and Huang4 developed a model of a four-bed pressure swing air separation unit. They showed that

10.1021/ie990181h CCC: $18.00 © 1999 American Chemical Society Published on Web 09/09/1999

3762

Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999

Pressurization and blowdown occur at the beginning of the applicable step and the product is drawn even while the producing bed is pressurizing. Thus, the product flow rate oscillates. A pressure regulator on the product end serves to dampen the fluctuation of the product flow rate. The following typical assumptions applied to the development of the bed model:

Figure 1. Simple two-bed PSA unit.

by modeling the valves, they could favorably predict pressure versus time. The transient case was not of interest in this study however. In the works of Jianyu and Zhenhua5 and Crittenden et al.,6,7 computation of the bed pressure history is pursued analogous to an electric circuit. To predict the dynamic response of a PSA unit to changes in any process variable, the scope of the model must include all beds which comprise the process, and the scope must extend beyond the adsorbent beds themselves. This approach is necessary so that the correct boundary conditions as a function of time are applied to the beds. At a minimum this scope must include the valves (or other pressure drops, such as the purge orifice in the case of OBOGS) that affect the flow rate of gas entering and leaving the beds. The purpose of this paper is to present a model which addresses computation of the bed pressure history and prediction of transient behavior. We develop a PSA model in which the typical mass and energy conservation equations are coupled (1) with an ODE which enforces overall mass conservation to compute the pressure history and (2) with models of the valves and purge orifice to compute the gas flow in and out of the beds. These flows are required by the overall balance equation. Computation of the bed pressure does not involve the Ergun equation because the pressure drop across the beds is negligible. This assumption is examined and justified. The PSA model has sufficient scope to predict the dynamic response of the unit to changes in the input variables. Mathematical Model A two-bed OBOGS is probably the simplest kind of PSA unit, a reflection of its application. See Figure 1. There are only two distinct steps in the operation of this unit, i.e., during the first half-cycle, production is on bed 1 while bed 2 is purged countercurrently. During the second half-cycle, the role of each bed is simply reversed.

1. Heat and mass transport in the radial direction is neglected and therefore transport only in the axial direction is considered. 2. Gas-phase axial dispersion can be lumped into a single axial dispersion coefficient. 3. The pressure swing air separation process is equilibrium-driven. Macropore diffusion provides the major resistance to mass transfer; under such conditions the simple linear driving force (LDF) assumption can be employed.2,8-10 4. Nitrogen is adsorbed in bulk quantity, and thus the change in the gas-phase velocity along the axial dimension is significant. The gas-phase material balance must be written accordingly.2,8-10 5. The gas phase is ideal. 6. The beds are adiabatic. 7. The molar heat capacity of the adsorbed phase is equal to that of the gas phase. 8. Pressure drop across the bed is negligible.11 An additional assumption was employed to simplify the derivation of the conservation equations: The mass and energy conservation equations were written to conserve molar quantities rather than mass quantities. The simplification is reasonable since the molecular weights of the conserved components (oxygen and nitrogen) are similar. The impact of this assumption on the mass-balance closure and on the results will be discussed. Justification for Negligible Pressure Drop Assumption. Sundaram and Wankat11 examined in detail the importance of pressure drop across the bed during pressurization and blowdown. By combining Darcy's Law with the continuity equation, they derived a dimensionless constant, τ, that varies directly with pressurization time, the adsorbent particle size, and initial column pressure, but varies inversely with gasphase viscosity and column length. For pressurization steps where the value of τ exceeds 0.5, pressure gradients can be assumed to be insignificant. Using the criteria established by Sundaram and Wankat, considering an OBOGS bed packed with a (16 × 40)-mesh 13X zeolite, having a 1/2-m length, and having a minimum pressurization time of 5 s (observed experimentally), τ g 140. Physically, this result indicates that the pressure drop across the beds is negligible compared to pressure drops that exist upstream and downstream from the beds, and thus, the bed pressure can safely be assumed to be uniform during the entire cycle. For this reason, conservation of momentum does not need to be considered in the model of the bed(s). Nevertheless, models of the upstream and downstream pressure drops are still needed because their characteristics determine the flow of gas through the process. Bed Mass Balances. The component balance for the gas phase is developed from a mole balance on a cylindrical differential element of the adsorbent bed:

Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3763

∂Cm ∂qm ∂2Cm ∂(uCm) (1) + Fs(1 - ) )+ DL ∂t ∂t ∂z ∂z2

q/1 ) K1p1 + mc,2 mc,1

∑ ∑ j)0 i)0

The total component balance is

∂(uC) ∂q ∂C + Fs(1 - ) )∂t ∂t ∂z

(

(

mc,2 mc,1(K

∑ ∑ j)0 i)0

M

∑ qm

(3)

m)1

P RT

Substitution and expansion of derivative terms yields

∂u u ∂T 1 ∂q dP ∂T + Fs(1 - ) )+ (5) P dt T ∂t C ∂t ∂z T ∂z The solid- (adsorbed-) phase component mole balance is

(6)

It has been assumed that the rate of exchange of mass between the gas and solid phase can be computed with a linear driving force rate expression in which km is the mass-transfer coefficient for component m.12 Bed Energy Balance. An energy balance on the cylindrical differential element of the adsorbent bed yields the following energy conservation statement:

[CCvg + Fs(1 - ) (Cpa + qCvg)]

(∑ M

Fs(1 - )

m)1

qdiff m

)

∂qm ∂t

Adsorption Isotherm. The equilibrium uptake in molecules/cavity versus pressure of a pure gas adsorbing on zeolite crystals can be modeled by the SSTM equation:13

(10)

for i + j g 2, iβ1 + jβ2 e v and where

-Q0,m RT

(11)

This is the multicomponent form of the SSTM, which is needed to calculate the uptake of each component from a multicomponent gas mixture. The parameters that are unique for each gas-zeolite crystal combination are K0, Q0, β, and mc. They are determined from pure component equilibrium data using the pure component form of the SSTM (eq 8). Commercial zeolites, like those used in this work, are zeolite crystals that are pelletized with a binder.15 The mass fraction of any binder must be taken into account when converting the uptake in molecules/cavity computed from eq 10 to the desired value in mol/g:

q0m ) q0(1 - bf)q/m

(12)

where bf is the mass fraction of the binder in the pelletized form. Overall Material Balance. The bed pressure history can be computed from an overall mole balance, given the molar flow of gas in to and out of the beds:

∂T

∂t ∂T ∂2 T ∂u ) -uCCvg -P + λg ∂z ∂z ∂z2 (7)

)

- iβ1/v - jβ2/v)(i+j)

Km ) K0,m exp (4)

∂qm ) km(q0m - qm) ∂t

i j 1p1) (K2p2) (1

i!j!

If the gas phase is assumed to be ideal, then

C)

(i -1)!j!

1 + K1p1 + K2p2 +

(2)

where

q)

)/

(K1p1)i(K2p2)j(1 - iβ1/v - jβ2/v)(i+j)

l ∂q dz + AtFs(1 - )∫0 dz ) FIN - FOUT ∫0l ∂C ∂t ∂t B

At

B

(13) This equation expresses the fact that the net change in gas- and adsorbed-phase inventory, integrated over the bed, is equal to the net flow into the bed. Introducing (4), this equation, in terms of pressure, is

At

dz + ∫0l RT l C ∂T ∂q At ∫0 (+ Fs(1 - ) ) dz ) T ∂t ∂t

dP dt

B

B

q* ) Kp + (Kp)2(1 - 2β/v)2 + ‚‚‚ + 1 + Kp +

(Kp)mc (1 - mcβ/v)mc (mc - 1)!

(Kp)2 (Kp)mc (1 - 2β/v)2 + ‚‚‚ + (1 - mcβ/v)mc 2! mc! (8)

where

K ) K0 exp

-Q0 RT

(9)

For a two-component system, the equilibrium uptake of component 1 is14

FIN - FOUT (14)

In eq 14, the derivative of pressure is with respect to time only since the pressure drop across the bed is assumed negligible. Pressure Drop Model. The inlet and outlet gas flows for each bed are determined by upstream and downstream pressure drops, respectively. These pressure drops are either the feed/exhaust valves or the purge orifice. For simplicity, a single pressure drop model was sought for all pressure drops. An empirical model to which pressure drop data could be fit was derived from an expression describing flow through an isothermal tube:

3764


π W) 8

x ()

2 2 d5 (Pu - Pd ) fl RT

6

10

(15)

in which W is the mass flow rate through the valve, d and l are the diameter and length of the tube, respectively, and f is a friction factor which is a complicated function of Pu and Pd. An empirical form, which can be used to compute the molar flow rate and can be fitted to experimental data, is obtained by dividing both sides by the gas molecular weight MW and defining an empirical discharge coefficient

CD )

x( ) d5 fl

(16)

to obtain

x

π F ) CD 8

10

(Pu2 - Pd2)

6

(17)

MWRT

The discharge coefficient CD is simply fit to the appropriate molar flow versus pressure drop data to model the feed/exhaust valve or purge orifice. The pressures Pu and Pd are the upstream and downstream pressures respectively that apply for the pressure drop. The upstream and downstream locations are based upon the assumed sign convention for positive flow. Thus, if Pu < Pd, then the following model applies:

x

π F ) -CD 8

(Pd2 - Pu2) 106 MWRT

(18)

System Model. The remaining unknowns are the (2) molar gas flows in to and out of each bed (F(1) IN, FIN, (1) (2) FOUT, and FOUT). The computation of these values is cycle-dependent. For the first half-cycle,

tively. By convention, positive flow is in the direction from feed end to product end of either bed. For computation of the purge flow, the upstream pressure is supplied by the regenerating bed and the producing bed supplies the downstream pressure. Thus, the purge flow typically has a negative sign, except for a short period of time following a cycle switch when the regenerating bed has a higher pressure than that of the producing bed. To limit the solution times, the two-component form of the model was considered for this work (M ) 2). Equations 1 (for m ) 1), 5, 6 (for m ) 1, 2), 7, and 14 represent 2M + 2 equations in the unknown state variables Cm (m ) 1), u, qm (m ) 1, 2), T, and P. Thus, for the two-component system, we solve one PDE for each of C1, u, q1, q2, and T and a single ODE (eq 14) for P. Both integrals in eq 14 can be evaluated since all variables appearing in the integrals are known functions of the state variables T(z,t), qm(z,t), and P(t). Equations 3, 4, 10, 11, and 12 are used to evaluate the appropriate terms in the equation system. This system must be solved for both beds simultaneously. Equations 19-23, or 24-28 are used to evaluate the inlet and outlet gas flows. The resulting system is fully specified, given the feed pressure (PSUP), exhaust pressure (PEX), product demand flow rate (FPROD), and appropriate initial and boundary conditions for the beds, which will be discussed next. Initial Conditions. For each bed, the following initial conditions apply:

P(0) ) Pinit Cinit ) Cm,init ) ym,initCinit

Pinit RTinit

(29) (30)

(m ) 1, ..., M - 1) (31)

u(z,0) ) 0

(32)

(1) F(1) IN ) fv(PSUP,P )

(19)

(2) F(2) IN ) fv(PEX,P )

(20)

F(1) OUT ) FPROD - FPURGE

(21)

T(z,0) ) Tinit

(34)

F(2) OUT ) FPURGE

(22)

Tw(z,0) ) Tw,init

(35)

FPURGE ) fpurge(P(2),P(1))

(23)

The exception is the initial uptake. However, the adsorbed phase is initially in equilibrium with the gas phase. Therefore, the initial uptake, qm, is computed with the equilibrium isotherm. The initial gas-phase composition is assumed to be that of air. Boundary Conditions. The boundary conditions on each bed depend on the mode. The two distinct modes of operation imposed on a bed during an entire cycle are pressurization/production (first half-cycle) and blowdown/purge (second half-cycle). But at the beginning of the first half-cycle, gas enters the pressurizing bed at both ends. This part of the half-cycle will be called pressurization. Eventually, the pressurizing bed pressure exceeds that of the purging bed, and the gas flow at the product end reverses, so that flow through the entire bed is in the same direction throughout (from feed to product end). This part of the half-cycle will be called production. Similarly, at the beginning of the second

For the second half-cycle, (1) F(1) IN ) fv(PEX,P )

(24)

(2) F(2) IN ) fv(PSUP,P )

(25)

F(1) OUT ) FPURGE

(26)

F(2) OUT ) FPROD - FPURGE

(27)

FPURGE ) fpurge(P(1),P(2))

(28)

where fv(Pu,Pd) is computed with (17) or (18), depending upon whether Pu - Pd is positive or negative, respec-

qm(z,0) ) q0m(y1,initP, ..., yM,initP)

(m ) 1, ..., M) (33)


half-cycle, gas exits the bed at both ends. This step is called blowdown. Finally, when the pressurizing bed pressure exceeds that of the purging bed, gas flow through the purging bed is in the same direction throughout (from product to feed end). This step is called purge. Mathematically, the pressurization step is distinguished from the production step by the sign of the flow at the product end. Similarly, the blowdown step is distinguished from the purge step by the sign of the flow at the product end. It is necessary to identify all four separate modes of operation so that the proper boundary conditions are applied to the beds. Pressurization and blowdown boundary conditions are applied immediately at the start of a new half-cycle. The boundary conditions are changed to that of production and purge, respectively, when FOUT on the appropriate bed changes sign.

Pressurization (FOUT < 0): -D -D

|

∂Cm ∂z

|

∂Cm ∂z

z)lB

z)0

) u(0,t)(ym,0-C(0,t) - Cm(0,t)) (36)

) -u(lB,t)(Cm(lB,t) - ym,l+BC(lB,t)) (37) u(0,t) )

|

FINRT AtP

(38)

∂T ) uCCpg(T0- - T(0,t)) -λg ∂z z)0

(39)

∂T -λg ) -uCCpg(T(lB,t) - Tl+B) ∂z z)lB

(40)

|

Production (FOUT > 0):

|

∂Cm -D ∂z

z)0

) u(0,t)(ym,0-C(0,t) - Cm(0,t)) (41)

|

∂Cm ∂z

u(0,t) )

|

)0

(42)

FINRT AtP

(43)

z)l

∂T -λg ) uCCpg(T0- - T(0,t)) ∂z z)0

(44)

∂T )0 ∂z z)lB

(45)

|

Blowdown (FOUT > 0): ∂Cm ∂z

| |

∂Cm ∂z

(46)

)0

(47)

FOUTRT AtP

(48)

z)l

u(l,t) )

| |

)0

z)0

Purge (FOUT < 0): ∂Cm ∂z -D

|

∂Cm ∂z

z)lB

|

z)0

(51)

) -u(lB,t)(Cm(lB,t) - ym,l+BC(lB,t)) (52) u(l,t) )

-λg

)0

|

FOUTRT AtP

∂T )0 ∂z z)0

|

∂T ) -uCCpg(T(lB,t) - Tl+B) ∂z z)lB

(53) (54) (55)

Note that the quantities u(0,t) and u(lB,t) are both computed during the pressurization step, but the latter is only used to apply the boundary condition for Cm and T at lB. Equation 5 requires specification of u at exactly one end at all times, regardless of the bed’s mode of operation. The values of ym,0- and T0- are determined by the composition and temperature of the feed and only apply to pressurization/production. The values of ym,l+ and Tl+ are determined by the composition and temperature of gas exiting the alternate bed, i.e., (1) (2) ym,0 - ) ym,0- ) ym,FEED

(56)

(2) (1) ym,l + ) ym (l,t)

(57)

(2) (1) ym,l + ) ym (l,t)

(58)

T0(1)- ) T0(1)- ) TFEED

(59)

(2) + ) T (lB,t) Tl(1) B

(60)

(1) + ) T (lB,t) Tl(2) B

(61)

Numerical Solution Dimensionless variables were defined for the independent variables t and z and for the state variables Cm (m ) 1), u, qm (m ) 1, 2), T, and P. These variables were used to express eqs 1, 5, 6, 7, and 14 in dimensionless form. The axially varying state variables Cm, u, qm, and T were discretized in the axial dimension using the Galerkin finite element technique,16 coded in FORTRAN. Further details have been documented.17 Integrals arising from this formulation were evaluated with quadrature rules. The integrals in eq 14 were evaluated piecewise over each finite element, also using quadrature rules. Thus, eq 14 is a single ODE in dP/dt which must be solved simultaneously with the ODE system arising from the finite element discretization of Cm, u, qm, and T. A uniform mesh with 100 elements was employed for the discretization. The discretized model can now be expressed in the residual form:

∂T )0 ∂z z)0

(49)

G(Y, Y4 ,t) ) 0

∂T )0 ∂z z)lB

(50)

In addition to the ODE system (62), the ODE solver requires that the vectors Y and Y4 be initialized with values that satisfy (62). Thus, at t ) 0 it is desirable to

(62)

3766


specify boundary conditions for which Y4 ) 0 to avoid having to root-solve (62) for Y4 at t ) 0. Therefore, in the numerical solution, feed pressure, exhaust pressure, and product flow rate are changed from initial values satisfying Y4 ) 0 to the desired values over a short finite time after t ) 0. These initial values are PSUP ) PEX ) Pinit and FPROD ) 0. Thus,

PEX ) PEX,spec - (PEX,spec - Pinit)e-At

(63)

PSUP ) PSUP,spec - (PSUP,spec - Pinit)e-At

(64)

FPROD ) FPROD,spec - (FPROD,spec - FPROD,init)e-At (65) where typically FPROD,init ) 0 and PSUP,spec, PEX,spec, and FPROD,spec are the specified values for the supply pressure, exhaust pressure and product draw rate, respectively. A is a time constant that determines the duration over which these values are changed from their initial values to the specified values. A value of A ) 0.5 s was found to be practical since smaller values resulted in longer integration times with no impact on the simulation results. The ODE solver also requires a routine which computes the Jacobian J of G at a given time t. The ADIFOR system developed by Bischof et al.18 at Argonne National Laboratory and Rice University was employed for this computation. The ADIFOR system takes as input the FORTRAN subroutine describing the function G and outputs FORTRAN code that computes gradients of G. The output code is derived through symbolic chain differentiation of the input code. These gradients are then used to compute the Jacobian:

Jij )

∂Gi ∂Gi +γ ∂Yj ∂Y˙

(66)

j

The Jacobian evaluation using this approach is faster by 1 or more orders of magnitude and is exact to within machine precision, resulting in a more stable integration compared to an integration which employs a finitedifference evaluation of the Jacobian. Solutions to G were obtained by integration with DASSL19 in half-cycle intervals. At the beginning of each half-cycle, a valve state logical flag was toggled. The pressure upstream from each feed/exhaust valve and the boundary conditions at the feed end of both beds were assigned applicable values throughout the duration of the half-cycle, depending upon the valve state flag value. Step changes in any of the specified variables, PSUP, PEX, and FPROD, if applicable, were implemented at this time also. Such step changes were implemented in a manner analogous to eqs 63-65:

PEX ) PEX,newspec - (PEX,newspec - Pprevspec)e-A(t-t0) (67) PSUP ) PSUP,newspec - (PSUP,newspec - Pprevspec)e-A(t-t0) (68) FPROD ) FPROD,newspec (FPROD,newspec - FPROD,prevspec)e-A(t-t0) (69) where t0 is the time at which the specification change is implemented and typically coincides with the begin-

ning of a new half-cycle. In this way, the step change is implemented over a finite time after the start of the half-cycle. We remark that the values of PSUP and PEX obtained from eqs 67 and 68 are applied to the appropriate feed/exhaust valve, depending upon the valve state. Thus, starting with the second half-cycle, a pressure discontinuity is applied to the feed/exhaust valve equations nonetheless. The magnitude of this change is |PSUP - PEX|. This particular discontinuity was not found to hinder the integrator at sufficiently small step sizes. Boundary conditions at the product end of both beds are dependent upon the sign of FOUT, as shown in eqs 36-55. Since u ) 0 when FOUT changes sign, the boundary conditions applied at the product end of the bed remain continuous. The average CPU time on a IBM RS/6000 running AIX per simulated PSA cycle was 13 min. Following step changes in any of the specified variables, 10-30 cycles were required for a new cyclic steady state to be reached. Thus, on average, simulation of dynamic responses of the PSA unit to all process variables took 4.25 h of CPU time. Experimental Section A two-bed PSA unit similar in size to that of an onboard oxygen generation system was constructed to test the predictive capabilities of the dynamic model. In particular, the PSA unit and supporting equipment were constructed so that it was possible to implement step changes in the following input variables: product flow, feed pressure, and cycle time. By design, nearly true steps in the variables were possible. The product oxygen composition, feed consumption, temperature of the bed at the axial midpoint (radial center), and the input variables (feed pressure and temperature, product flow rate, and cycle time) were all recorded as a function of time. The two-bed PSA unit which was constructed is illustrated in Figure 2. The beds were stainless steel pipe, 24.5-in. long with a 2-in. inside diameter and 16gauge wall thickness. The beds were loaded with UOP Oxysiv-5 zeolite (type 13X). The zeolite was loaded from a 50-lb airtight container in which the zeolite is shipped, and thus activation of the zeolite was not pursued. The dew point of the feed air was maintained at -40 °F with a PSA air-drying system employing activated alumina as the adsorbent. The beds, solenoid valves, and purge orifice were all plumbed with 1/4-in. National Pipe Thread (NPT) stainless fittings. The purge orifice was microbore tubing (0.028-in. i.d.; 1/8-in. o.d.) having a length of 6 in. The 1/ -in. feed piping accommodated feed flows up to 150 2 standard liters per minute (SLPM) of air. Because the product flow was less than 20 SLPM, the piping beyond the product-end solenoid valves was 1/4-in. stainless steel plumbed with 1/4-in. Swagelok. To implement step changes in feed pressure, the process flow scheme from the compressed air supply to the feed point was duplicated, as shown in Figure 2. Two solenoid spool valves were activated to implement a switch from one feed source (at the desired initial pressure) to the second feed source (at the desired final pressure). Step changes in the product flow rate were implemented by manually changing the setpoint on the product flow control valve. The PSA cycle was driven


Figure 2. Two-bed PSA unit. Table 1. Specifications for Pressure and Flow Measurement Devices (units of bar or SLPM) model no.

measured variable

manufacturer

input

accuracy

linearity

BHL-4202-01-05 FM3925V FC260

feed and bed 1 pressure

TransInstruments Tylan Tylan

0-6 0-1000 0-5

(0.006 (30 (0.05

(0.006 (10 (0.025

flow

by a Transistor-Transistor Logic (TTL) compliant square wave signal supplied to the solenoid valve driver by a function generator. Step changes in cycle time were implemented by manually changing the frequency of the function generator. All step changes were implemented at a time coincident with the beginning of a new halfcycle. Figure 2 illustrates the process variables that were measured and recorded as a function of time. These variables included the flow rate and temperature of each feed, the temperature and pressure at the feed point, product flow rate and temperature, the pressure of bed 1, and the composition of the product gas. Temperature measurements also were obtained at the center radial position at distances of 5 in. (front), 1215/16 in. (middle), and 2015/16 in. (rear) from the feed end of bed 1. Specifications for the pressure transducers, feed flow-

meters, and product flow controllers are shown in Table 1. Temperature measurements were obtained with Omega K-type thermocouples. The composition of the product gas was measured with a Perkin-Elmer medical gas analyzer. All measurement signals were analogfiltered so that frequencies above 20 Hz were attenuated. By design of the data acquisition and control system, a typical data set consists of process data from the first 21/2 cycles of the initial cyclic steady state, followed by dynamic data associated with the process variable step change until the new cyclic steady state was reached. A suitable operating region was identified experimentally to determine an operating region (cycle time, product flow, and feed pressure) in which a significant response of the product oxygen composition to a change in product flow, feed pressure, and cycle time could be

3768


Table 2. Summary of Experimental PSA Data cyclic run steady no. state

cycle time (s)

feed feed product pres. temp. flow (bar) (°C) (mol/s)

Table 4. Model Parameters

average product composition O2

N2

Ar

1

initial 10.6329 2.89 final 10.6329 2.87

24.8 24.6

0.00056 0.956 0.000 0.044 0.00267 0.496 0.479 0.025

2

initial 10.4364 1.02 final 10.4364 4.48

24.0 24.0

0.0 0.210 0.780 0.010 0.00178 0.706 0.263 0.031

3

initial 10.4254 4.49 final 10.4254 3.55

24.0 24.8

0.00178 0.706 0.263 0.031 0.00178 0.589 0.384 0.026

4

initial 9.5933 4.40 final 13.3205 4.47

24.8 25.0

0.00178 0.608 0.365 0.027 0.00178 0.808 0.156 0.036

5

initial 17.4333 4.48 final 21.1433 4.48

25.2 25.3

0.00187 0.880 0.082 0.039 0.00179 0.952 0.006 0.042

6

initial 35.5790 1.02 final 35.5790 4.69

23.5 23.5

0.0 0.210 0.780 0.010 0.00249 0.883 0.077 0.040

7

initial 40.6667 4.66 final 40.6633 4.66

25.9 25.9

0.00224 0.917 0.042 0.041 0.00105 0.958 0.000 0.042

Table 3. SSTM Parameters

gas

m

β (Å3)

K0 (molecules/ (cavity bar))

nitrogen oxygen argon

14 32 38

82.1 35.2 29.8

4.41 × 10-4 1.23 × 10-3 1.74 × 10-3

Q0 (kcal/mol) 4.51 3.14 2.91

observed. This operating region was found to correspond to a feed pressure of 50 psig, product flow rate of 0-5 SLPM, and a cycle-time of 10-40 s. Summary of PSA Runs. Approximately 25 experimental runs were conducted for which feed pressure, product flow, and cycle time were stepped. The experimental conditions selected were modest deviations within the identified operating region. From these runs, seven cases were selected for which the most significant unit dynamics were observed. These cases were used to test the model and are summarized in Table 2. As shown in Table 2, two of the experimental runs involve start-up of the PSA unit from idle conditions to a cyclic steady state. The remaining five runs involve step changes of one of the specified variables, PSUP and FPROD, and cycle time. Process variables on the unit change from an initial cyclic steady state to a final cyclic steady state. Parameter Estimation Isotherm Parameters. Equilibrium data for components of air adsorbing on type 13X zeolite were determined by Miller.20 Pure and multicomponent equilibrium data were obtained using the volumetric method. Table 3 shows the parameter values for 13X which Miller determined from the pure component equilibrium data. The pure component data were determined over the temperature range of -70-50 °C and over a pressure range of 0-5 bar. The capacity factor q0 for type 13X zeolite crystals was found by Miller20 to be 0.000556 (mol/g)/(molecules/ cavity). Miller also found that 5A and 13X type zeolites manufactured by UOP have a binder mass fraction of 20-24%.20-21 A binder mass fraction of 20% was assumed in this work. Differential Adsorption Heats. Assuming the gas phase is ideal, the differential heat of adsorption, qdiff, is related to the isosteric heat of adsorption, Q0, by the following equation:22

parameter packed bed length, lB (cm) packed bed radius, Rb (cm) zeolite solid density, Fs (g/cm3) zeolite void fraction, zeolite capacity factor, q0 ((mol/g)/(molecules/ cavity)) zeolite binder mass fraction, bf zeolite cage volume, v (Å3) zeolite pellet diameter, Rp (mm) axial dispersion coefficient, DL (cm2/s) axial thermal conductivity, λL (cal/(cm s °C)) O2 mass transfer coefficient, kO2 N2 mass transfer coefficient, kN2 O2 molecular weight, MWO2 N2 molecular weight, MWN2 O2 heat capacity, Cp,O2 N2 heat capacity, Cp,N2 bed heat capacity, Cpa (cal/(g °C) diff (cal/mol) differential adsorption heat of N2, qN 2 diff differential adsorption heat of O2, qO2 (cal/mol)

qdiff ) Q0 - RT

value 62.87 2.347 1.53 0.53 0.000556 0.2 1150 1 20 1.66 × 10-3 122 46 32.0 28.0 6.96 7.03 0.19 3918 2548

(70)

Using the pure component isosteric heats of adsorption, Q0,m, based on the fit of the isotherm data to the SSTM equation (see Table 3), the differential heats of adsorption for nitrogen and oxygen at 25 °C are shown in Table 4. Bed Physical Parameters. The solid particle density, Fs, of 13X molecular sieve crystals was measured by Miller.20 Munkvold experimentally determined the void fraction, , of Oxysiv-5.10 The adsorbent bulk density is given by

Fb ) Fs(1 - )

(71)

The measured length of the packed section, lB, and the bed radius, Rb, yield the packed volume, VB. The mass of the packed section can thus be computed from Fb. It was found in this work that the computed mass of the packed section in each bed agreed to within 2% of the measured values, providing confidence in the parameter values used, which are shown in Table 4. Gas-Phase Parameters. The gas phase is assumed by the model to be ideal, and therefore the only required parameters are the molecular weight and ideal gas heat capacity. The heat capacities are further assumed by the model to be constant. The values reported by Felder and Rousseau23 for 0 °C were used. Molecular Diffusivity. Literature correlations were used to estimate the effective axial diffusivity, the effective axial thermal conductivity, and the masstransfer coefficient for each component. The best model predictions were obtained with constant values for each of these parameters, computed at the purge conditions of 1 bar and 25 °C. The molecular diffusivity for the pair oxygen-nitrogen was needed for each of the computations, for which the theory of Chapman and Enskog, described by Bird et al.,24 was employed. Axial Dispersion Coefficient. The effective axial diffusivity, D, in eq 1 can be estimated using a correlation described by Ruthven:12

D ) γ1Dm + γ22Rpu

(72)

where γ1 and γ2 are typically 0.7 and 0.5, respectively. Rp is the adsorbent particle radius and u is the superficial velocity (see eq 1). On the basis of an average


velocity of 14 cm/s (determined by simulation) and adsorbent pellets having a diameter of 1 mm, an axial diffusion coefficient of 1 cm2 /s is obtained. This value, however, was revised so that the slope of an experimentally observed nitrogen breakthrough curve matched the simulated curve. The simulation employed mass-transfer coefficients computed in a manner discussed in the next section, and the breakthrough experiments were conducted at this average velocity. A diffusion coefficient of 20 cm2/s was found to model the nitrogen breakthrough curve most accurately. Details about the breakthrough study are discussed in a subsequent section. Mass-Transfer Coefficients. The rate of mass exchange between the gas and adsorbed phases in a zeolite-based air separation process is controlled by the diffusivity of each gas in the macropores of the zeolite pellets, as noted by Ruthven and co-workers.2 For this case a linear driving force rate expression can be applied, for which the mass-transfer coefficient is given by

km )

ΩDe,m Rp2

Figure 3. Fit of CD to the feed/exhaust valve pressure drop data.

(73)

When used in PSA modeling, the value of Ω is 15,25 except for PSA cycles having a cycle time much shorter than those used in this work. The effective diffusivity, De,m, is given by

De,m )

Dp,m dq0m 1 + (1 - p) dCm

(74) Figure 4. Fit of CD to the purge orifice pressure drop data.

in which the particle diffusivity, Dp,m, can be estimated by12

Dp,m )

Dm τ

(75)

In the study of PSA air separation by Farooq et al.8 using 0.7-mm 5A pellets, it was assumed that p ) 0.33 and τ ) 3. In this study, 1-mm 13X pellets ((16 × 40)mesh Oxysiv-5) were assumed. Considering that both crystal types are pelletized by similar processes and that the pellet diameters are similar, it was assumed that the same values of p and τ could be used. Axial Thermal Conductivity. An estimate for the gas-phase axial thermal conductivity follows from heatand mass-transfer analogy:12

λg ) CCvgD

(76)

Using this equation, an order of magnitude estimate of 1 × 10-3 cal/(cm s °C) was inferred by assuming an average bed velocity of 14 cm/s at purge conditions, based on simulation results. The estimate is 2 orders of magnitude greater than the thermal conductivity of the zeolite, thus supporting the assumption that solidphase conduction is negligible compared to that in the gas phase. Bed Heat Capacity. Literature data for the heat capacity of Oxysiv-5 in particular could not be located; however, a value of 0.19 cal/(g °C) was used by Chou and Huang4 for type 5A zeolite. This value was found to be a reasonable estimate for Oxysiv-5.

Feed/Exhaust Valve Parameters. As discussed, the feed/exhaust valves were modeled with eq 17. The molar flow of air through one valve was measured over a range of feed pressures ranging from 1 to 5 bar. The discharge coefficient, CD, was fitted to these data using the nonlinear parameter estimation program GREG.26 A discharge coefficient of CD,v ) 0.01 was identified. Figure 3 illustrates the fit of eq 17 to the pressure drop data. Purge Orifice Coefficient. In a manner similar to that of the feed/exhaust valves, the molar flow of air through the purge orifice was measured over a range of feed pressures. A fit of the pressure drop data suggested a value of CD,purge ) 0.018, as shown in Figure 4. The fit is not as favorable as that for the feed/exhaust valve, and there is a dependence on the flow direction. The data from both flow directions were used to obtain an average fit. Alternative pressure drop models were considered but not found to yield better overall fit of both data sets. Subsequent simulations of the initial cyclic steady state of all selected experimental runs with this coefficient value yielded predictions of the product oxygen/argon composition which were consistently low by 0.1-0.2 mole fraction. However, adjustment of the purge orifice coefficient was found to yield agreement between the predicted and experimentally observed cycle-averaged product oxygen/argon composition at the initial cyclic steady state of each experimental run. This observation could indicate that the installed characteristics may differ from the uninstalled characteristics. The value of the purge orifice coefficient which was used for each run is shown in Table 5. Parameter Evaluation by Breakthrough Simulation. To verify that physically meaningful parameters

3770


Table 5. Model Predictions of the Average Product Oxygen/Argon Compositions, Average Feed Consumptions and Average Bed Temperatures at the Middle/Center Position for All Runs at Initial and Final Cyclic Steady States; Model vs Experiment cyclic steady state

purge orifice coeff.

1

initial final

2

run no.

average product O2/Ar compositions

average feed consumption (mol/s)

average bed temp. at middle/center position (°C)

exp.

model

exp.

model

exp.

model

0.004 0.004

1.000 0.514

0.919 0.437

0.0474 0.0475

0.0532 0.0536

25.8 26.0

28.3 28.4

initial final

0.0025 0.0025

0.220 0.733

0.220 0.735

0.0000 0.0866

0.0000 0.0871

22.5 27.3

22.5 29.4

3

initial final

0.0025 0.0025

0.733 0.612

0.735 0.608

0.0864 0.0598

0.0871 0.0689

27.4 25.6

29.4 27.8

4

initial final

0.0025 0.0025

0.637 0.842

0.676 0.863

0.0896 0.0740

0.0886 0.0729

27.6 27.3

30.9 30.7

5

initial final

0.0025 0.0025

0.919 0.994

0.942 0.978

0.0611 0.0525

0.0592 0.0508

27.5 27.4

31.1 31.3

6

initial final

0.0025 0.0025

0.220 0.921

0.220 0.933

0.0000 0.0362

0.0000 0.0376

23.5 27.3

23.5 29.5

7

initial final

0.0025 0.0025

0.958 1.000

0.926 0.998

0.0362 0.0329

0.0331 0.0316

26.4 26.0

31.9 32.4

were obtained, a brief study was conducted in which oxygen/nitrogen breakthrough experiments were conducted and compared to simulation results. One bed from the PSA unit was fed pure oxygen at 2.94 bar and 21.3 °C at a constant rate of 47.3 SLPM, until thermal equilibrium was reached. The feed was then switched to pure nitrogen at a constant feed temperature, feed pressure, and downstream valve position. A new steadystate was reached at 2.96 bar and 21.5 °C feed conditions and a constant flow of 44.4 SLPM. The feed was then switched back to pure oxygen until again the bed reached thermal equilibrium. Responses of the effluent composition, feed flow rate, bed pressure, and radial center temperature of the bed at front, middle, and rear axial positions were measured as functions of time. The model of this apparatus consisted of a one-bed model, a model of the upstream solenoid valve lying between the constant pressure source and the bed, and a model of the downstream flow control valve, beyond which the bed effluent was vented to atmospheric conditions. The discharge coefficients of these two valves were estimated from the experimentally observed pressure drop across each valve at final steady-state conditions. The model assumed constant feed pressure and temperature and constant valve characteristics. The flow rate of each feed source was measured upstream with 100-SLPM flowmeters manufactured by Sierra Instruments. The exhaust flow was measured with a 50-SLPM Tylan flowmeter. These flowmeters introduced negligible pressure drop into the system. As shown in Figure 5, the single-bed model was found to predict reasonably the nitrogen breakthrough time and the shape of the breakthrough curve. The net change in bed inventory was determined by integrating the difference between the feed consumption rate of nitrogen and exhaust flow rate from initial steady state to final steady state. This net inventory change agreed with the model-predicted value to within 2%. Similar agreement was found when the feed was switched back to pure oxygen. These observations suggest that the isotherm parameters and mass-transfer coefficients are reasonable. As discussed, a revised axial diffusion coefficient was identified. Figures 6 and 7 show the experimentally observed and computed responses of bed temperatures during nitrogen and oxygen breakthrough, respectively. The

Figure 5. Dynamic response of the product nitrogen compositions to a step change in the feed nitrogen mole fractions from 0 to 1 at 0 s; simulation vs experiment.

Figure 6. Dynamic response of the bed temperatures to a step change in feed nitrogen mole fractions from 0 to 1 at 0 s; simulation vs experiment.

magnitude, propagation rate, and dispersion of the thermal wave are all predicted reasonably by the model, establishing confidence in the computation of bed thermal behavior, bed heat capacity, and axial thermal conductivity. We note that the single-bed model also included a bed wall energy balance, with appropriate terms in both


Figure 7. Dynamic response of the bed temperatures to a step change in the feed oxygen mole fractions from 0 to 1 at 0 s; simulation vs experiment.

Figure 9. Dynamic response of the bed temperatures to a step change in the feed nitrogen mole fractions from 0 to 1 at 0 s; simulation vs experiment, constant qdiff ) 3233 cal/mol.

Figure 8. Computed temperature profiles at the end of the first half-cycle for cycles 1, 2, 6, 11, 12, 16, 17, 20, and 21 for run 1: CD,valve ) 0.1; CD,purge ) 0.003.

Figure 10. Dynamic response of the bed temperatures to a step change in the feed oxygen mole fraction from 0 to 1 at 0 s; simulation vs experiment; constant qdiff ) 3233 cal/mol.

energy balances needed to model the change in bed temperature due to heat loss to the surroundings. In Figure 9, the bed temperature decay over a period of 200 s (the typical length of a PSA simulation) is only 1 °C. Thus, we concluded that the PSA beds could be assumed to be adiabatic. This assumption is further supported by the PSA data which will be discussed. Energy-Balance Simplification. Predictions of the bed thermal behavior for the PSA unit with the model as formulated were unexpected, however. In particular, the model predicted that the temperature of the bed near the rear end decreases, presumably without bound, with each subsequent cycle. This prediction is illustrated in Figure 8. As discussed earlier, the conservation statements were formulated for simplicity to conserve molar quantities. Thus, the overall mass balance is not exactly satisfied. Consequently, the energy balance does not close exactly. A small error in the mass- and energybalance closure over the course of a single cycle may, however, compound over multiple cycles to yield the unexpected behavior. A slightly simplified energy balance was found to eliminate the unrealistic behavior:

where qdiff is an average differential adsorption heat value determined by averaging the pure component differential adsorption heat values. This average value is 3233 cal/mol. Figures 9 and 10 illustrate the impact of this simplification on the predictions of bed thermal behavior for the breakthrough simulations. Comparison of these figures with the more rigorous predictions suggested that the impact of this simplification on the predictions of PSA dynamic behavior would not be significant.

[CCvg + Fs(1 - ) (Cpa + qCvg)] Fs(1 - )qdiff

∂T ∂t

∂q ∂T ∂2T ∂u ) -uCCvg + λg 2 - P (77) ∂t ∂z ∂z ∂z

Results and Discussion Table 2 summarizes the dynamic data obtained with the PSA unit. Table 5 summarizes the model predictions of the average product oxygen/argon composition, average feed consumption, and average bed temperature at the middle/center position compared to the experimentally observed values for the initial and final cyclic steady state of all runs. Refer to Table 2 for the experimental conditions of each run. The model versus experiment comparisons in Figures 11-22 illustrate the process dynamics which can be modeled accurately with the model and those process dynamics which would require a much more rigorous model to predict. The comparison of simulation and experimental results generally revealed that the best agreement was with the cycle-averaged value of most

3772


Figure 11. Dynamic response of product O2/Ar compositions to the unit start-up; simulation vs experiment (run 2).

Figure 12. Dynamic response of the product O2/Ar compositions to the unit start-up; simulation vs experiment (run 6).

process variables. Cycle averaging refers to integration of the process variable over the duration of one cycle divided by the cycle time. As revealed by Figures 17 and 18, the bed pressure history is correctly predicted at any point in time at both short and long cycle times. Although the actual purge rate cannot be ascertained, the good prediction of pressure history over the entire duration of the cycle suggests that a reasonable computation of the purge rate is obtained if the purge orifice coefficient can be correctly estimated. A reasonable prediction of the detailed product oxygen composition response is possible at sufficiently long cycle times (greater than 15 s), as revealed by Figure 12. At shorter cycle times, the measurement delay is a sufficiently significant percentage of the cycle time so that the prediction of the product oxygen composition at a specific point in time is somewhat in error, as shown in Figure 11. Figure 15 reveals that detailed prediction of the feed consumption versus time is not possible with the model as formulated; the predicted amplitude of the oscillation is considerably larger than that observed. Additional models of upstream pressure drops are necessary to obtain more acceptable prediction. From Figures 19 and 20 it is apparent that the model incorrectly predicts the amplitude of the bed temperature oscillation. The decrease in average bed temperature versus time, as shown in Figures 21 and 22, is very

Figure 13. Dynamic response of the cycle-averaged product O2/ Ar compositions to the unit start-up; simulation vs experiment (run 2).

Figure 14. Dynamic response of the cycle-averaged product O2/ Ar compositions to the unit start-up; simulation vs experiment (run 6).

Figure 15. Dynamic response of the feed flow rate to the unit start-up; simulation vs experiment (run 2).

small, indicating that heat loss to the surroundings does not account for the smaller observed amplitude in bed temperature swing relative to that predicted by the model. Rather, a more probable explanation is that there is a significant unmodeled temperature gradient within the zeolite pellets at the end of the pressurization/production step. This explanation is further supported by the fact that the observed amplitude is larger


Figure 16. Dynamic response of the cycle-averaged feed flow rate to a step increase in the cycle time; simulation vs experiment (run 5).

Figure 19. Dynamic response of the bed temperatures at the middle/center position to the unit start-up; simulation vs experiment (run 2).

Figure 17. Dynamic response of bed 1 pressure to the unit startup; simulation vs experiment (run 2).

Figure 20. Dynamic response of the bed temperatures at the middle/center position to the unit start-up, simulation vs experiment (run 6).

Figure 18. Dynamic response of the bed 1 pressure to the unit start-up; simulation vs experiment (run 6).

for the longer cycle time (Figure 20) compared to the short cycle time case (Figure 19). Consequently, predictions of the cycle-averaged bed temperature typically were higher, by 2-6 °C. Table 5 generally reveals that, with reasonable correlation of the cycle-averaged product oxygen composition at the initial cyclic steady state, subsequent predictions of the dynamic response of cycle-averaged product oxygen composition were generally excellent

Figure 21. Dynamic response of the cycle-averaged bed temperatures at the middle/center position to the unit start-up; simulation vs experiment (run 2).

(within 3 mol %). As discussed, for cycle times greater than 15 s, detailed predictions of the product oxygen composition vs time also were reasonably accurate. Predictions of the cycle-averaged feed flow rate versus time at cyclic steady state and the dynamic response of cycle-averaged feed flow rate versus time typically agreed with experimentally observed values to within 5%. However, two cases involved low-pressure cyclic

3774


Figure 22. Dynamic response of the cycle-averaged bed temperatures at the middle/center position to the unit start-up; simulation vs experiment (run 6).

steady states where the error was as high as 15%. (Error was greatest for the lowest pressure case, run 1.) In the development of the mathematical model, it was noted that the conservation equations were written to conserve molar quantities rather than mass quantities. Satisfaction of these statements was verified by checking that the simulation output for each PSA bed satisfied the following mole balance:

∆n )

∫0t FIN dt - ∫0t FOUT dt f

f

(78)

and the following statement of the first law of thermodynamics:

∆U )

∫0t Hˆ g,FEEDFIN dt - ∫0t Hˆ g,PRODFOUT dt f

f

(79)

These statements were found to be satisfied to less than 0.1%. In contrast, the overall mass balance

∆W )

∫0t WIN dt - ∫0t WOUT dt f

f

(80)

was found to be 0.5% in error after 535 s in run 6. While this error certainly is reasonable, it is a possible contributing factor to the unexpected predictions with the more rigorous energy balance.

flows were computed by models of the pressure drops lying upstream and downstream from the bed. A simplification to the bed energy balance, in which a constant differential heat of adsorption is employed, was necessary to obtain physically realistic predictions of the bed temperature. It is likely that this simplification is necessary since the extensive material balance fails to conserve mass by 0.5%. Literature correlations were used to estimate most parameters in the model and found to be physically reasonable by single-bed breakthrough experiments and simulation. The feed/exhaust valve discharge coefficient and an initial estimate for the purge orifice discharge coefficient were determined from pressure drop data. Predictions of cyclic steady-state bed pressure versus time and dynamic response of bed pressure versus time were in agreement with experimental data. Using the purge orifice coefficient estimated from pressure drop data, predictions of the cycle-averaged product oxygen composition versus time for a given cyclic steady state were low by 10-20 mol %. However, excellent correlation of these data were possible by adjusting the purge orifice coefficient. With reasonable correlation of the cycle-averaged product oxygen composition at cyclic steady state, subsequent predictions of the dynamic response of cycle-averaged product oxygen compositions were generally excellent (within 3 mol %). For a cycle time greater than 15 s, detailed prediction of the product oxygen composition versus time also was possible. Predictions of the cycle-averaged feed flow rate versus time at cyclic steady state and the dynamic response of cycle-averaged feed flow rate versus time typically agreed with experimentally observed values to within 5%, with the exception of two cases at low feed pressure where the error was as high as 15%. Only cycleaveraged values of the feed flow rate are useful. The detailed dynamics of the feed flow rate computed by the model do not reflect reality. The simulations suggested that the magnitude of the temperature swing over a cycle is significantly greater than that observed experimentally. Because heat loss to the surroundings cannot account for this difference, we postulated that a significant unmodeled temperature gradient within the zeolite pellets exists at the end of the pressurization/production step. Prediction of cycleaveraged bed temperature versus time was achieved to within ( 2 °C.

Conclusions A computer model of a small nonisothermal pressure swing air separation process was formulated, solved, and experimentally validated to address technical needs related to the design and operation of on-board oxygen generation systems (OBOGS). The scope of the model includes the two adsorbent beds, feed/exhaust solenoid valves, and the purge orifice. We have demonstrated that the model is capable of predicting the dynamic response of product composition, feed consumption, bed temperature, and bed pressure to step changes in input variables of interest. These input variables are feed pressure, product draw rate, and cycle time. The conservation equations describing each bed consist of M - 1 component balance PDEs, a total component balance PDE (which is solved for the velocity profile), an extensive material balance ODE which is solved for the bed pressure history, and a bed energy balance. The extensive material balance requires the molar gas flows entering and exiting the bed. These

Acknowledgment Financial support from The United States Air Force and McDonnell Aircraft Company are gratefully acknowledged. Nomenclature Roman Characters k ) mass-transfer coefficient of a single adsorbate (s-1) l ) length (cm) lB ) adsorbent bed length (cm) mc ) maximum number of molecules (an integer) that can be adsorbed by a zeolite cavity mc,i ) maximum number of component i molecules (an integer) that can be adsorbed by a zeolite cavity p ) partial pressure (bar) p0m ) partial pressure of component m corresponding to equilibrium conditions with uptake qm of component m (bar)

Ind. Eng. Chem. Res., Vol. 38, No. 10, 1999 3775 q ) total adsorbate uptake (mol/g) qm ) uptake of component m (mol/g) q* ) total uptake when the gas and adsorbed phases are at equilibrium (molecules/cavity) q0 ) total uptake when the gas and adsorbed phases are at equilibrium (mol/g) q/m ) uptake of component m when the gas and adsorbed phases are at equilibrium (molecules/cavity) q0m ) uptake of component m when the gas and adsorbed phases are at equilibrium (mol/g) qdiff ) differential heat of adsorption (cal/mol) t ) time (s) u ) superficial velocity (cm/s) v ) volume of zeolite cage (Å) At ) bed total cross-sectional area (cm2) C ) total gas-phase concentration (mol/cm3) Cm ) gas-phase concentration of component m (mol/cm3) Cpg ) constant pressure molar heat capacity of the gas phase (cal/(mol K)) Cvg ) constant volume specific heat capacity of the gas phase (cal/(mol K)) Cpa ) constant pressure molar heat capacity of the adsorbent (cal/(g K)) CD ) valve discharge coefficient D ) axial diffusion coefficient (cm2/s) Dm ) molecular diffusivity of component m (cm2/s) De,m ) effective diffusivity of component m (cm2/s) Dp,m ) particle diffusivity of component m (cm2/s) F ) molar flow (mol/s) J ) Jacobian of the system of equations G K ) Henry's Law constant (cm3/g) Km ) Henry's Law constant for component m (cm3/g) M ) number of components MW ) molecular weight (g/mol) MW ) average molecular weight (g/mol) P ) pressure (bar) Pu ) upstream pressure (bar) Pd ) downstream pressure (bar) PSUP ) supply pressure (bar) PEX ) exhaust pressure (bar) P0 ) total pressure corresponding to equilibrium conditions at total adsorbate uptake q (bar) Q0 ) isosteric heat of adsorption in the SSTM equation (kcal/mol) Q0,m ) isosteric heat of adsorption of component m in the SSTM equation (kcal/mol) R ) universal gas constant, 83.14 (cm3 bar)/(mol K) Rb ) radius of adsorbent bed (cm) Rp ) adsorbent pellet radium (cm) T ) temperature (°C or K) W ) mass flow rate (g/s) Greek Letters β ) molecular volume of adsorbate (Å3) ) bed void fraction (dimensionless) p ) particle void fraction (dimensionless) λg ) thermal conductivity of gas phase in the axial dimension (cal/(cm s K)) Fb ) bulk density of zeolite (g/cm3) Fs ) zeolite solid density (g/cm3) σ ) diffusive flux (mol/(cm2 s)) σT ) heat flux due to conduction in the axial dimension through the bed (cal/(cm2 s)) τ ) tortuosity factor (dimensionless)

Literature Cited (1) Shendalman, L. H.; Mitchell, J. E. A Study of Heatless Adsorption in the Model System CO2 in He, I. Chem. Eng. Sci. 1972, 27, 1449. (2) Ruthven, D. M.; Farooq, M. S.; Knaebel, K. S. Pressure Swing Adsorption; VCH Publishers: New York, 1994. (3) Beaman, J. J. A Dynamic Model of a Pressure Swing Oxygen Generation System. J. Dyn. Syst., Meas., Control 1985, 107, 111. (4) Chou, C.-T.; Huang, W.-C. Simulation of a Four-Bed Pressure Swing Adsorption Process for Oxygen Enrichment. Ind. Eng. Chem. Res. 1994, 33, 1250. (5) Jianyu, G.; Zhenhua, Y. Analog Circuit for Simulation of Pressure Swing Adsorption. Chem. Eng. Sci. 1990, 45, 3063. (6) Crittenden, B. D.; Guan, J.; Ng, W. N.; Thomas, W. J. Dynamics of Pressurization and Depressurization During Pressure Swing Adsorption. Chem. Eng. Sci. 1994, 49, 2657. (7) Crittenden, B. D.; Guan, J.; Ng, W. N.; Thomas, W. J. Pressure, Concentration and Temperature Profiles in a 5A Zeolite Adsorbent Bed During Pressurisation and Depressurisation With Air. Chem. Eng. Sci. 1995, 50, 1417. (8) Farooq, S.; Ruthven, D. M.; Boniface, H. A. Numerical Simulation of a Pressure Swing Adsorption Oxygen Unit. Chem. Eng. Sci. 1989, 44, 2809. (9) Fernandez, G. F.; Kenney, C. N. Modeling of the Pressure Swing Air Separation Process. Chem. Eng. Sci. 1983, 38, 827. (10) Munkvold, G. D. A Multicomponent Predictive Model for Pressure Swing Adsorption Applied to Air Separation. Ph.D. Dissertation, University of Texas at Austin, Austin, TX, 1990. (11) Sundaram, N.; Wankat, P. C. Pressure Drop Effects in the Pressurization and Blowdown Steps of Pressure Swing Adsorption. Chem. Eng. Sci. 1988, 43, 123. (12) Ruthven, D. M. Principles of Adsorption and Adsorption Processes; Wiley: New York, 1984. (13) Ruthven, D. M. Sorption and Diffusion in Molecular Sieve Zeolites. Sep. Purif. Methods 1976, 5, 189. (14) Ruthven, D. M. Sorption of Oxygen, Nitrogen, Carbon Monoxide, Methane, and Binary Mixtures of These Gases in 5A Molecular Sieve. AIChE J. 1976, 22, 753. (15) Breck, D. W. Zeolite Molecular Sieves; Wiley: New York, 1973. (16) Becker, E. B.; Carey, G. F.; Oden, J. T. Finite Elements: An Introduction; Prentice Hall: Englewood Cliffs, NJ, 1981. (17) Teague, K. G., Jr. Predictive Dynamic Model of a Small Nonisothermal Pressure Swing Air Separation Process. Ph.D. Dissertation, University of Texas at Austin, Austin, TX, 1999. (18) Bischof, C.; Carle, A.; Corliss, G.; Griewank, A.; Hovland, P. ADIFORsGenerating Derivative Codes from Fortran Programs; Argonne Preprint MCS-P263-0991; Argonne National Laboratory: Argonne, IL, 1995. (19) Brenan, K. E.; Campbell, S. L.; Petzold, L. R. Numerical Solution of Initial-Value Problems in Differential-Algebraic Equations; North-Holland: New York, 1989. (20) Miller, G. W. Adsorption of Nitrogen, Oxygen, Argon, and Ternary Mixtures of These Gases in 13X Molecular Sieve. AIChE Symp. Ser. 1987, 83, 28. (21) Miller, G. W.; Knaebel, K. S.; Ikels, K. G. Equilibria of Nitrogen, Oxygen, Argon, and Air in Molecular Sieve 5A. AIChE J. 1987, 33, 194. (22) Ross, S.; Olivier, J. P. On Physical Adsorption; Wiley: New York, 1964. (23) Felder, R. M.; Rousseau, R. W. Elementary Principles of Chemical Processes; Wiley: New York, 1978. (24) Bird, R. B.; Stewart, W. E.; Lightfoot, E. N. Transport Phenomena; Wiley: New York, 1960. (25) Raghaven, N. S.; Hassan, M. M.; Ruthven, D. M. Numerical Simulation of a PSA System Using a Pore Diffusion Model. Chem. Eng. Sci. 1986, 41, 2787. (26) Caracotsios, M. Model Parametric Sensitivity Analysis and Nonlinear Parameter Estimation: Theory and Applications. Ph.D. Dissertation, University of Wisconsin, Madison, WI, 1986.

Received for review March 10, 1999 Revised manuscript received July 26, 1999 Accepted July 28, 1999 IE990181H

Predictive Dynamic Model of a Small Pressure Swing Adsorption Air

Recommend Documents