Ind. Eng. Chem. Res. 1998, 37, 1435-1443
1435
Modeling and Optimization of Continuous Adsorption in a Dessicant Rotor Frantisˇ ek S ˇ teˇ pa´ nek,†,‡ Milan Kubı´cˇ ek,‡,§ and Milosˇ Marek*,† Department of Chemical Engineering, Department of Mathematics, and Center for Nonlinear Dynamics of Chemical and Biological Systems, Prague Institute of Chemical Technology, Technicka´ 5, 166 28 Praha 6, Czech Republic
A mathematical description of a continuous adsorption process in the so-called “dessicant rotor” based on the assumption of a locally adiabatic operation is developed. The steady state of the cyclic process is formulated as a mixed boundary value problem for the governing spatially onedimensional parabolic PDE’s and solved by the shooting method. The efficiency of several quasiNewton schemes is compared with that of the fixed-point iteration method. An algorithm for the determination of operational parameters satisfying given constraints and minimizing the amount of heat supplied to the process is demonstrated as an example. Introduction Separation processes involving fluid-solid interactions face the problem of limited mobility of the solid phase and are therefore usually operated in a periodic, unsteady-state manner. Pressure and temperature swing adsorptions (PSA and TSA) can serve as typical examples of these processes (see, e.g., Ruthven et al., 1994). If such a unit operation is to be connected within a larger system, it is often desirable that it operate continuously. This can be accomplished by periodically switching between two or more nonstationary units, simulating a countercurrently moving bed (Tonkovich and Carr, 1996), or designing an apparatus which operates truly continuously. Such arrangements for affinity extraction of proteins have been studied by Rodrigues et al. (1992). Bridges and Baker (1992) have considered continuous liquid chromatography separations, and Byers and Williams (1996) investigated the feasibility of continuous ion exchange. Another possible implementation of a continuous adsorption process is the so-called sorption wheel, schematically shown in Figure 1. According to Keller (1995), the process enjoys widespread industrial utilization, with its major applications being gas purification (e.g., VOC removal) and air-drying (cf. Munters (1992), the developers of the process). The latter application is used as an example in this work. To the best of our knowledge, no comprehensive mathematical treatment of this process has appeared in the open literature to date, with the exception of Konrad and Eigenberger (1994), who studied a stepwise operated rotary adsorber. We present a basic, spatially one-dimensional model which is structurally equivalent to the description of adsorption in a fixed-bed column regenerated by temperature swing. The main difference is that, unlike in most simulation works in the field of regenerative adsorption processes (cf., e.g., Liu and Ritter, 1996), we * Author to whom correspondence is addressed. E-mail:
[email protected]. Phone: +4202 2435 3104. Fax: +4202 311 7335. † Department of Chemical Engineering. ‡ Center for Nonlinear Dynamics of Chemical and Biological Systems. § Department of Mathematics.
Figure 1. Outline of the desiccant rotor and the corresponding process scheme (HX is a heat exchanger).
use a specific formulation, similar to the situation in practice, where limits are imposed upon the composition of the outlet streams while some of the process parameters (e.g., the desorption gas temperature or the cycle time) are unknown. This specific formulation requires a new computational approach. The solution procedure described here is also applicable to the determination of key operational parameters of the classical TSA and PSA processes. The Principle of Operation A slowly rotating monolithic wheel (Figure 1) passes through two separate zones in which streams of air flow
S0888-5885(97)00478-8 CCC: $15.00 © 1998 American Chemical Society Published on Web 02/24/1998
1436 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
across the wheel in opposite directions (a similar configuration for a monolithic reactor was used by Bakker et al., 1996). In the upper zone, adsorption takes place and a dried air stream leaves the wheel. Desorption occurs in the lower zone in which hot purge air regenerates the adsorbent. The desorption zone typically comprises less than half of the wheel. Either a part of the dried air or the humid atmospheric air (or a mixture of the two) can be used for purge purposes. For given conditions of the humid inlet air (its temperature, relative humidity, and volumetric flowrate), the operation of the wheel is determined by four parameters: the period of rotation and the composition, temperature, and flowrate of the desorption stream. The performance of the process is characterized by the amount of moisture removed from the air and the cost (mainly energy consumption) per unit of dry air leaving the apparatus at steady state. The process is said to be at its steady state when the concentration and temperature profiles at any angular position within the wheel do not change with time. A flowsheet representation of the process along with the notation used for variables in the text is displayed in Figure 1. Model Formulation The rotor is formed by a monolith with a honeycomblike internal structure, containing a large number of small channels. The walls of the channels are covered by a moisture-adsorbing agent (e.g., silica gel in the case of air-drying), which is integrated directly into the matrix of the monolith. We therefore adopt a pseudohomogeneous model which characterizes the wheel by a macroscopic porosity and assumes three state variablessconcentration in the gas phase, denoted c, concentration in the solid phase q, and temperature T (assuming thermal equilibrium between the fluid and the solid phase). A rigorous description of the wheel would lead to partial differential equations in three spatial dimensionssa computationally rather demanding problem. Assuming locally adiabatic operation of the wheel (i.e., neglecting temperature gradients in both radial and angular directions) allows us to reduce the model into a spatially one-dimensional form, which permits the parametric studies and optimization to be performed within reasonable CPU times. The effect of thermal coupling in the angular direction will be treated in a separate study. We chose the time t that it takes for a volume element of the wheel to rotate from an angular position φ ) 0 (set to be at the start of the adsorption compartment) to some angular position φ as the independent variable. If τp is the total period of rotation, we have
φ ) 2π(t/τp)
(1)
Considering further the usual assumptions (cf. Tondeur, 1995) such as mass transfer described by the linear driving force, plug flow with linear velocity unaffected by trace component removal, zero pressure drop, and temperature independence of solid phase density and heat capacities, we arrive at the following three partial differential equations, structurally similar to those used to describe periodic fixed-bed processes:
∂(vc) 1 - ∂q ∂c )∂t ∂z ∂t
∂q ) k(q* - q) ∂t
(F C g
p,g
+
∂(FgvT) 1- ∂T FsCp,s ) -Cp,g + ∂t ∂z 1- ∂q (-∆H) (2) ∂t
)
The interstitial gas velocity, v, depends on temperature according to
v ) v0(T/T0)
(3)
as follows from the continuity condition, the assumptions of constant total pressure P, and the ideal gas law. Similarly, the temperature dependence of the gas-phase density Fg is
Fg ) Fg,0(T0/T)
(4)
In the above equations, the subscript 0 refers to conditions at inlet. Other quantities, namely, the mean heat capacities Cp,g and Cp,s, the mean density of the solid phase Fs, and the overall mass transfer coefficient k, are assumed to be temperature independent. Equilibrium and mass-transfer data for the model were those used by Park and Knaebel (1995), who measured breakthrough curves for the water-air-silica system. They found the adsorption equilibrium for this system to be well-described by the Dubinin-Ashtakhov isotherm, which is type IV according to Brunauer’s classification and has the form
q* ) q1e-(Ea/E1) + q2e-(Ea/E2) 2
2
(5)
where q1 and q2 (assumed constant) are capacities of the adsorbent attributed to multilayer sorption and pore filling, respectively, and E1 and E2 are the respective energies of binding. Ea is the adsorption potential energy given by
Ea ) -RT ln(φ)
(6)
where φ is the relative humidity defined by sat φ ) pH2O/pH 2O
(7)
The saturated vapor pressure depends on temperature according to the Antoine equation sat ))Aln(pH 2O
B C+T
(8)
The isosteric heat of adsorption, which is the sum of an adsorption-specific term and the condensation heat of water, is given by
∆H ) -Ea -
BRT2 (C + T)2
(9)
as derived and verified on experimental data by Park and Knaebel (1995). The values and units of all parameters appearing in the equations above are listed in Table 1. A modification was made to the total adsorbent capacity, qm ) q1 + q2. We reduced qm by a factor of 10 to account for the fact that the monolith of
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1437 Table 1. Values of Fixed Parameters Used in the Model Fg Fs Cp,g Cp,s k
1.164 kg/m3 1.35 × 103 kg/m3 1.02 × 103 J/kg‚K 7.94 × 102 J/kg‚K 0.018 s-1
q1 q2 E1 E2
32.123 × 102 mol/m3 8.5193 × 102 mol/m3 3.4434 × 103 J/mol 10.931 × 103 J/mol 0.67
the wheel does not entirely consist of an active adsorbent material. The wheel rotates with a certain speed; we denote τp as the period of rotation. The time an element of the wheel spends in the adsorption zone is τa and in the desorption zone τd; hence, it holds that τa ) ra/dτp and τd ) (1 - ra/d)τp, where ra/d is a design parameter determining the relative sizes of the adsorption and desorption compartments. When τa is known (determined by a method described in the next section), the desorption period τd can then be calculated:
1 - ra/d τd ) τa ra/d
(10)
The solution of system (2) must fulfill periodically switched boundary conditions in space:
(t mod τp) ∈ 〈0; τa)
(t mod τp) ∈ 〈τa; τp)
{ {
c(t,0) ) ca T(t,0) ) Ta v|z)0 ) va c(t,L) ) cd T(t,L) ) Td v|z)L ) -vd
(11)
(12)
In addition, we want the 1-periodic solution (cyclic steady state) of the set of PDE’s (2) to fulfill the following boundary conditions in time:
q(τd,z)0) ) rprgq*(ca)
T(0,z) ) T(τp,z)
(13)
In the current formulation, τp is a priori unknown; it is determined adaptively for each set of operating parameters so that a certain level of drying of the inlet air be achieved in the adsorption compartment:
∫0 y(t,L) dt ) rbrkya
1 τa
τa
(14)
We have denoted ya and yj the average mole fractions of water vapor in the humid air and the product air, respectively, rbrk the drying parameter, and L the depth of the wheel. The flow (linear velocity va) and properties (ca, Ta) of the inlet air are assumed to be fixed; ca is calculated from the relative humidity φa as it is a more common quantity to express the content of moisture in air. The properties of the desorption air stream are determined as follows. Let the desorption air be a mixture of a part of the dry air and the humid atmospheric air. Its composition is then given by
cd ) [rmixyj + (1 - rmix)ya]
P RTd
(15)
where rmix ∈ 〈0; 1〉 is the ratio of the molar flow of the dry air fraction to the total molar flow of the desorption
(16)
where rprg is a parameter expressing the degree of regeneration of adsorbent in the desorption subperiod. As will be shown below, the choice of rprg significantly affects the performance of the process. Let us consider a situation where the goal is to determine the values of variable process parameters, given an inlet stream of known properties and an existing device (implying fixed design parameters) in such a way that certain requirements on the properties of the outlet stream be fulfilled and direct operating costs be minimized at the same time. For a given flowrate of the inlet stream, the direct operating costs are composed of three principal itemssdirect labor, power consumption associated with rotating the wheel, and heat supplied to an exchanger in which the desorption air is heated. While the first two items are more or less fixed, the last one gives a degree of freedom for optimization with respect to some variables appearing in the model formulated above, namely, Td. We therefore define the following dimensionless criterion CF to which the operating costs are proportional:
CF )
q(0,z) ) q(τp,z)
23.478 3.9849 × 103 K -39.724 K 0.75 0.35 m
stream. The flow of the desorption air is again determined adaptively in order to achieve a chosen level of purging of the wheel in the desorption compartment. We formulate the condition for determining vd so that
c(0,z) ) c(τp,z)
yj )
A B C ra/d L
nd Td - Ta n0 Ta
(17)
The molar flows nd and n0 (the net dry air stream leaving the process) are related to the respective linear velocities by
P nd ) S(1 - ra/d)vd RTd
(18)
P n0 ) Sra/dv0 RTa
(19)
and
where S is the cross-sectional area of the rotor and v0 is related to va and vd through rmix:
Ta 1 - ra/d v0 ) va - rmixvd Td ra/d
(20)
Substituting nd and n0 into eq 17 from the above equations leads to
CF )
vd 1 - ra/d Td - Ta v0 ra/d Td
(21)
which is the final form of the objective function to be minimized. For the sake of completeness, we also mention the relations between the partial pressure p, mole fraction y, relative humidity φ, and concentration c, because they are used throughout the computations:
1438 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 sat p ) yP ) φpH ) cRT 2O
(22)
They follow directly from the ideal gas law. Methods of Solution The direct determination of cyclic steady states of periodically operated fixed-bed processes was outlined by Croft and LeVan (1994a,b), who posed the problem as a periodic boundary value problem in time and solved it by the method of lines (MOL) and the shooting method, utilizing the Newton method for the resulting algebraic system with the Jacobi matrix obtained by integrating variational equations. Conceptually similar but employing the finite elements method instead of the finite difference method was the approach of Salinger and Eigenberger (1996), who studied a reverse-flow afterburner. In both of these cases, arc-length continuation (e.g., Kubı´cˇek and Marek, 1983) was performed to trace the dependence of the periodic solution and its stability on a parameter. If the only goal is to find a stable and unique 1-periodic solution, a computationally more efficient approach is to use one of the quasiNewton methods, as was done by Smith and Westerberg (1993) or Kvamsdal and Herzberg (1997). It has to be noted, however, that like the Newton method itself, the quasi-Newton methods can generally converge also to an unstable solution. When the method of lines is used for the solution of the PDE’s, a set of ordinary differential equations with periodically switched righthand sides results. The dynamics of such systems was formalized in the work of Klı´cˇ and Pokorny´ (1996) and studied in the case of a reverse-flow reactor by R ˇ eha´cˇek et al. (1998). The efficiency of possible computational approaches to such systems is discussed in Unger et al. (1997). Before we describe the overall solution procedure of the optimization problem, let us discuss the numerical methods that were selected for the individual subproblems. The core of the computational efforts is the solution of the periodic boundary value problem represented by eqs 2, 11, 12, and 13. We have used the method of lines; i.e., the spatial derivatives in eq 2 were replaced by finite difference formulas. If we arrange the state variables at all m + 1 equidistant spatial mesh points into a common vector x ) [c0, q0, T0, c1, q1, T1, ..., cm, qm, Tm]T ∈ Rn, n ) 3(m + 1), we can regard the integration of eqs 2 over one time cycle τp to define a discrete mapping P(x) on Rn, called the period map. The cyclic steady state (1-periodic solution) then corresponds to a fixed point of the period map. We thus have to solve the algebraic problem
F(x) ) P(x) - x ) 0
(23)
Let x* be a fixed point of the period map P (i.e., it satisfies F(x*) ) 0). Stability of x* as well as convergence properties of the iteration sequence
xk+1 ) P(xk)
k ) 0, 1, 2, ...
(24)
are determined by the eigenvalues of the linearization matrix of P(x) enumerated at x*, called the monodromy matrix and denoted M:
M ) ∂P(x)/∂x
(25)
If all eigenvalues of M lie within the unit circle in the
complex plane, x* is an asymptotically stable fixed pointsan attractorsand it holds that limkf∞ xk ) x* provided that x0 is chosen in the domain of attraction of x*. We will only consider this case from now on. The rate of convergence toward x* is determined by the socalled spectral radius F(M) ) max|λi|, where λi’s are the eigenvalues of M. The closer F(M) is to unity, the slower the convergence. As the evaluation of M is rather timeconsuming, it has not been explicitly used during our calculations. Instead, the stability of found solutions was tested in selected cases by performing a few iterations of the direct substitution method (see below). For the solution of eq 23, we have tested the fixedpoint iteration method (alternative names: direct substitution or Picard method), the Wegstein method and three quasi-Newton methodssnamely, those of Davidon-Fletcher-Powell, Broyden, and Zontendijk. For the direct substitution method, which is in fact a dynamical simulation of the original system the iteration sequence is represented by eq 24. This method has the advantage of being very robust; it always converges to x* if the initial estimate is chosen within its domain of attraction. The disadvantage of this method is that it is of first order of convergence which may imply quite a large number of iterations before a given precision is reached. As a convergence criterion, we have used the weighted root-mean-square (RMS) norm defined by
k )
x ( ) 1
n
∑
ni)1
2
xk+1 - xki i
(26)
wki
because the components of x differ in orders of magnitude (recall that x contains c, q, and T). The scaling factor wki is chosen to be proportional to xki . Convergence of the Picard method can be accelerated by the Wegstein method, known especially from flowsheeting simulations with recycle streams. Let xk, xk-1, and xk-2 be three consecutive iterations obtained by the direct substitution method. Then the new, extrapolated value of xk+1 is obtained according to
[
]
|xk - xk-1| xk+1 ) xk-1 + 1 - sk k-1 |x - xk-2|
-1
(xk - xk-1) (27)
Here sk ∈ 〈-1, 1〉 is a parameter determining the orientation and possible damping of the extrapolation in the kth step. The component-by-component extrapolation of the Wegstein method makes an implicit assumption that there is not a very strong coupling between the individual components of x (i.e., that M is close to a diagonal matrix, which is not the case here). Therefore, the Wegstein method only works well when the concentration and temperature profiles do not differ qualitatively from iteration to iteration. This is valid in the vicinity of the solution; otherwise, the Wegstein method leads to overshooting and divergence. Finally, three quasi-Newton methods were considered (e.g., Edgar and Himmelblau, 1988). Quasi-Newton methods differ from the classical Newton-Raphson iteration in that they do not require the Jacobi matrix JF of F(x) (equal to M - I) in an analytical form, but they build an approximation of the inverted Jacobi matrix from the knowledge of the history of the convergence process. This is their biggest advantage because obtaining M by integrating the variational equations is
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1439
Figure 2. Comparison of the rate of convergence of the solution methods tested, starting from a good initial estimate. The first initializating iteration is always done by the Picard method.
computationally rather expensive. For instance, if the dimension n of the original system were 300, we would have to integrate additional 300 × 300 ) 900 000 ordinary differential equations and then solve a system of 300 linear algebraic equations in order to take one Newton step. This would only pay off if one step of the Newton method were at least as efficient as 300 steps of the direct substitution method, which is not the case for the system considered here, as can be seen in Figure 2, curve 5 (cf. also the comparison made by Kvamsdal and Herzberg, 1997, p 828). A new iteration xk+1 is obtained by a quasi-Newton method according to the formula
xk+1 ) xk - RHkF(xk)
(28)
R ∈ (0, 1〉 is a damping factor and H is an approximation of the inverted Jacobi matrix JF-1 updated either according to the Davidon-Fletcher-Powell scheme
Hk+1 ) Hk +
∆xk(∆xk)T (∆xk)T∆Fk
-
Hk∆Fk(∆Fk)THk (∆Fk)THk∆Fk
(29)
the Broyden update rule
Hk+1 ) Hk +
(∆xk - Hk∆Fk)(∆xk - Hk∆Fk)T (∆xk - Hk∆Fk)T∆Fk
(30)
or the Zontendijk formula
Hk+1 ) Hk -
(Hk∆Fk)(Hk∆Fk)T ∆FTkHk∆Fk
(31)
In the above update rules, the ∆’s are defined as follows:
∆xk ) xk+1 - xk and ∆Fk ) F(xk+1) - F(xk) (32) In our simulations, we took H0 to be the unit matrix which ensures that the first step of the quasi-Newton methods is in fact equivalent to a direct substitution step if R ) 1. This has a positive effect on the robustness of the methods because, in adsorption processes with adaptive switching periods, the first few cycles are typically sufficient to generate temperature
and concentration profiles topologically close to the final solution. In the case when as the starting profile we take the end profile obtained for similar values of parameters, a single Picard iteration is sufficient to bring the state vector safely deep into the domain of attraction of the quasi-Newton methods. The rate of convergence of all the methods just described is compared in Figure 2, for a typical situation encountered during our optimization calculations, i.e., for a case where the iteration starts from a good initial estimate. In parametric studies, such an initial estimate is typically available from previous iteration runs performed with similar values of parameters. The starting value for the very first calculation can be obtained by the direct substitution method if the quasiNewton methods diverge. During the tests, we found the Zontendijk method to be superior in terms of both speed and stability of convergence (it converged even in some situations where the other quasi-Newton methods failed). That is why it was chosen as the solution method for the subsequent optimization calculations. The Optimization Procedure For given properties of the inlet air stream, characterized by its composition (expressed in terms of relative humidity φa), temperature Ta, and flowrate (expressed by the interstitial velocity va), we want to determine all process parameters (the period of rotation τp, the ratio of mixing of the dried air into the purge air rmix, and the purge air temperature Td and flowrate vd) which ensure the chosen product quality (characterized by the drying ratio rbrk) and minimize the objective function CF defined by eq 21. As already mentioned, instead of using vd as an optimization variable directly, we have introduced the more general purge factor rprg; cf. eq 16. This leaves three variable parameters to be determined: Td, rmix, and rprg. As will be shown below, both τp and vd are obtained as “byproducts” during enumeration of the objective function. For optimization with respect to rprg and rmix, we have used the systematic search method as we were interested not only in the actual position of the optimum but also in the influence of these two parameters on the shape of the objective function. As an alternative approach, a global optimum with respect to all three variables (Td, rprg, rmix) might first be sought for using a rigorous optimization method (e.g., Nelder-Mead) and then the optimum is perturbed in order to obtain information about the shape of the objective function. Due to the fact that CF does not exist in analytical form, this sensitivity study would still have to be made by systematically changing rprg and rmix in different directions around the optimum. Moreover, the BVP defined by eq 16 proved to be rather sensitive to the initial estimate of vd, a fact that favors the use of systematic search, as it permits a certain degree of “manual control” over the calculations. For a selected combination of rprg and rmix, the objective function CF was minimized with respect to Td by the quadratic search method. This method was chosen in order to minimize the number of necessary enumerations of the objective function as this is a lengthy iterative procedure, consisting of the following steps: 1. Starting from an initial state-variable profile x0, integrate (numerically) eqs 2 in time with boundary conditions (11). At each time step, evaluate the average
1440 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
mole fraction of water vapor in the outlet stream according to eq 14. Stop when the required value of rprgya is exceeded. 2. To determine τa exactly, return to the last time step (let t ) r and yj ) yjr) and integrate a transformed system with yj as an independent variable and an augmented vector [x, t] as the dependent variables from yjr to rbrkya. 3. Calculate τd from τa according to eq 10. Store τp ) τa + τd. 4. Determine vd in the following way (shooting method): (a) Choose an initial approximation vd0. (b) Integrate eqs 2 with boundary conditions (12) from τa to τp. (c) Update vd by the secant method. (d) Compare q(τp,z)0) with rprgq(ca). If the difference is greater than vd, return to step b. 5. Calculate the residual ||x(τp) - x0||. If it exceeds cs, update x0 by a quasi-Newton step and return to step 1. 6. The cyclic steady state is reached; calculate the objective function CF. In a typical case, there were 5-6 iterations of the vd loop (step 4), 3-4 iterations of the quasi-Newton method (step 5), and about 6 enumerations (recall that three are necessary just for the initial bracketing during the quadratic search method) of CF needed in order to determine the optimum Td with a relative error of 10-4 (four significant digits). Usual CPU times on an HP 700 series workstation were thus on the order of 50 min for one optimization run with respect to Td (i.e., with rmix and rprg fixed).
Figure 3. Dependence of the relative cost function on desorption air temperature for several values of inlet air relative humidity. The values of CF were scaled by φa to put them on a comparable basis (CF′ ) 100CF/φa).
Results and Discussion To demonstrate the functionality of the abovedescribed procedures, we have performed a series of optimization calculations and parametric studies. All results reported below were obtained for fixed values of design parameters (ra/d ) 0.75 and L ) 0.35 m), and flowrate and temperature of the inlet stream (Ta ) 298 K, va ) 0.225 m/s). The content of moisture in the inlet air was varied in the range of 40-80% in order to simulate real-life conditions. First, the effect of the desorption air temperature on the objective function CF was investigated for fixed values of rmix ) 1.0 (i.e., purge with dry air) and rprg ) 0.5. As can be seen from Figure 3 which displays results for a 95% air moisture removal, CF is a unimodal function of Td. For lower desorption temperatures, the total heat supplied per unit of dry air is higher despite the lower temperature. The reason for this is that higher flowrates (i.e., higher values of vd) are required when Td is low (recall that the heat flow is proportional to product of the temperature difference and the masssand hence volumetric flowrate). For the same reasons, when Td is too high, vd is smaller but its decrease is not sufficient to offset the effect of the increased temperature difference. The optimum then lies in the region where the combined influence of these two reverse effects is minimal. A trend that can be deduced from Figure 3 is that the desorption temperature corresponding to the minimum of CF increases with increasing φa. There are two reasons for this fact. First, the period of rotation is shorter for higher contents of moisture in the inlet air because the adsorbent’s capacity gets saturated earlier. This leaves less time for desorption, and a bigger driving force is necessary if a certain degree of adsorbent
Figure 4. Optimum desorption gas velocity vd* as a function of inlet air relative humidity φa for different levels of required drying (rmix ) 1.0 and rprg ) 0.5).
regeneration is to be achieved in the shorter desorption time. A higher driving force is induced by higher desorption temperature. The second reason is that not only is the desorption time shorter but more moisture (in absolute terms) has to be removed from the silica gel within this shorter period. The total result of these two effects is that there is a monotonic increase of the optimal Td with increasing φa. A confirmation of the above-stated qualitative observations was made by a parametric study whereby the values of operational parameters at optimum were followed as a function of φa. The study was made for three different levels of required drying, namely, rbrk ) 0.02, 0.05, and 0.10, which corresponds to the removal of 98, 95, and 90% of the moisture originally contained in the process air, respectively. Figure 4 reveals the results for vd. Similar upward-sloping trends also hold for Td as well as for the objective function CF; however, the rotation period decreases with increasing φa. An interesting phenomenon can be observed on the 0.10 branch. At φa equal to approximately 54%, a second local minimum emerges and attains a lower value of the objective function for φa ) 57%. The original local optimum ceases to exist when φa reaches
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1441
Figure 5. Transition between two local minima of CF as a function of Td with increasing φa (rbrk ) 0.10, rmix ) 1.0, and rprg ) 0.3).
about 61%, with the “new” optimum becoming unique from then on. While the optimum Td* was increasing with increasing φa below the 57% threshold, a sudden jump occurs on the transition to the new minimum, with Td* increasing from 117 to 195 °C. As φa continues to increase, Td* decreases and levels off at a value of 168 °C for φa greater than 70%. This is a remarkable phenomenon because in all other situations Td* normally increases with φa. On the vd* vs φa graph, the situation just described is visible as a “break” on the rbrk ) 0.10 curve in the specified φa range. The genesis of the transition between the two local optima is illustrated in detail in Figure 5. So far, rprg and rmix were kept fixed. Next we investigated the effect of the purge air composition and the degree of adsorbent regeneration on operating costs. Both effects are shown in Figure 6 for the case of 95%
moisture removal from inlet air with relative humidity φa ) 60%. Although the optimum desorption temperature Td* increases as rmix varies from 1.0 to 0.0 (i.e., from purge with dry air to wet air), the cost function moves in the opposite direction. For example, in the case mentioned above and rprg ) 0.5, CF is higher by20.5% when purge is made with dry air compared to purge with atmospheric air. The reason for this behavior is that, although sorbent regeneration with dry air means a higher driving force for mass transfer and could be thus perceived as “better”, CF was defined (recall from eq 21) to be proportional to heat supplied per unit of net product. When part of the product dry air is used for desorption, we arrive at higher relative costs. Thus, the conclusion is that, in this range of parameters and product quality requirements, it does not pay off to use the dry air for purge.
1442 Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998
Figure 6. Optimum desorption temperature Td* as a function of rprg for different compositions of the desorption air, determined by rmix (φa ) 60% and rbrk ) 0.05).
Next, the effect of the purge ratio on operating costs was investigated. It was shown by Rousˇar and Ditl (1993) that complete purge in the desorption step of a PSA cycle does not necessarily lead to maximum product purity. In our thermal regeneration air-drying example, operating costs are considered and it is also evident that neither complete purge nor an excessive amount of adsorptive left in the wheel leads to optimum performance. By definition, complete purge would imply infinitely large vd which means infinite CF, as follows from eq 16 and the linear driving force approximation for mass transfer. On the other hand, large values of rprg imply short rotation periods. If at least partial regeneration of the adsorbent is to be made within this short τd, vd and/or Td must be high, again leading to higher values of the objective function. It can be seen in Figure 6 that, for all values of rmix, the dependence Td* vs rprg goes through a single minimum and so does the objective function. With the inclusion of rmix and rprg, it can be stated that all operating parameters were covered and a global optimum was found in the parameter space (rmix ) 0 and rprg ) 0.52 in this particular case). The situation is somewhat different when a higher level of drying is required (e.g., rbrk ) 0.02). In the rotor configuration we have considered, there is no cooling step (Davis and LeVan, 1989) connected in the cycle and this results in early breakthrough due to thermal carryover. As the requirement is imposed upon the average composition of the product air, the early breakthrough is offset by dryer air coming out of the rotor in the central part of the adsorption compartment. However, for higher humidities of the inlet air the thermal carryover effect is too strong to be offset and the only means of achieving the required level of drying is by using the dry air for purge. Mathematically, it is possible to view this situation as a constrained optimization problem where, with respect to rmix, the optimum lies at the boundary of the feasible region. Concluding Remarks The optimization procedure described above serves as a tool for determining process parameters that fulfill given requirements on product quality and at the same time minimize direct operating costs. As the calculation
of a global optimum still requires larger CPU times than what is feasible for process control in real time, we see the main use of the algorithm to be in precalculating the optimum process parameters (rotation times and purge air temperature and flowrate) for a series of inlet air properties and flowrates and several values of the drying ratio rbrk and using such data in the form of graphs or tables for setting control parameters of an adsorption wheel. Finally, it has to be mentioned that adsorption in a sorption wheel as a process lends itself to further quantitative investigation. For example, testing the validity of the locally adiabatic model used in this work against a rigorous 2D or 3D description of the sorption monolith would be an interesting effort. Dynamic control (based on the Pontrjagin maximum principle) of a transient behavior of the rotor when inlet conditions change is another area yet to be explored. A whole class of new problems could also be represented by evaluating novel designs of the sorption wheel featuring, e.g., multiple compartments that would allow for the inclusion of a cooling step or the withdrawal of product gas fractions of different quality. Acknowledgment This work was partially supported by Grant No. 104/ 97/0602, Czech Grant Agency, and Grant No. VS96073, Czech Ministry of Education. Literature Cited Bakker, W. J.; Vriesendrop, M.; Kapteijn, F.; Moulijn, J. A. Sorbent Development for Continuous Regenerative H2S Removal in a Rotating Monolith Reactor. Can. J. Chem. Eng. 1996, 74, 713. Bridges, S.; Baker, P. E. Modelling Continuous Chromatographic Separations. Chem. Eng. Sci. 1992, 47, 1299. Byers, Ch. H.; Williams, D. F. Efficient Recovery of Lanthanides by Continuous Ion Exchange. Ind. Eng. Chem. Res. 1996, 35, 993. Croft, D. T.; LeVan, M. D. Periodic States of Adsorption CyclessI. Direct Determination and Stability. Chem. Eng. Sci. 1994a, 49, 1821. Croft, D. T.; LeVan, M. D. Periodic States of Adsorption CyclessII. Solution Spaces and Multiplicity. Chem. Eng. Sci. 1994b, 49, 1831. Davis, M. M.; LeVan, M. D. Experiments on Optimization of Thermal Swing Adsorption. Ind. Eng. Chem. Res. 1989, 28, 778. Edgar, T. F.; Himmelblau, D. M. Optimization of Chemical Processes; McGraw-Hill: Singapore, 1988. Keller, G. E. Adsorption: Building Upon a Solid Foundation. Chem. Eng. Prog. 1995, 11, 56. Klı´cˇ, A.; Pokorny´, P. On Dynamical Systems Generated by Two Alternating Vector Fields. Int. J. Bif. Chaos. 1996, 6, 2015. Konrad, G.; Eigenberger, G. Rotating adsorber for air purification and solvent recovery (Rotoradsorber zur Abluftreinigung und Lo¨sunsmittel-ru¨ckgewinnung). Chem.-Ing.-Tech. 1994, 66, 321. Kubı´cˇek, M.; Marek, M. Computational Methods in Bifurcation Theory and Dissipative Structures; Springer-Verlag: New York, 1983. Kvamsdal, M. H.; Herzberg, T. Optimization of PSA Systemss Studies on Cyclic Steady-State Convergence. Comput. Chem. Eng. 1997, 21, 819. Liu, Y.; Ritter, J. A. Pressure Swing AdsorptionsSolvent Vapor Recovery: Dynamics and Parametric Study. Ind. Eng. Chem. Res. 1996, 35, 7. Munters GmbH. Desiccant Rotors; Company product information catalogue: Hamburg, 1992. Park, I.; Knaebel, K. S. Adsorption Breakthrough Behavior: Unusual Effects and Possible Causes. AIChE J. 1992, 38, 660.
Ind. Eng. Chem. Res., Vol. 37, No. 4, 1998 1443 R ˇ eha´cˇek, J.; Kubı´cˇek, M.; Marek, M. Periodic, Quasiperiodic and Chaotic Spatiotemporal Patterns in a Tubular Catalytic Reactor with Periodic Flow Reversal. Comput. Chem. Eng. 1998, 22, 283. Rodrigues, M. I.; Zagor, C. A.; Maguer, F.; Asenjo, J. A. Dynamic Modelling, Simulation and Control of Continuous Adsorption Recycle Extraction. Chem. Eng. Sci. 1992, 47, 263. Rousar I.; Ditl, P. Pressure Swing Adsorption: Analytical Solution for Optimum Purge. Chem. Eng. Sci. 1993, 48, 723. Ruthven, D. M.; Farooq, S.; Knaebel, K. Pressure Swing Adsorption; VCH Publishers: New York, 1994. Salinger, A. W.; Eigenberger, G. The Direct Calculation of Periodic States of the Reversed Flow ReactorsI. Methodology and Propane Combustion Results. Chem. Eng. Sci. 1996, 51, 4903. Smith, O. J.; Westerberg, A. W. Acceleration of Cyclic Steady-State Convergence for Pressure Swing Adsorption Models. Ind. Eng. Chem. Res. 1992, 31, 1569.
Tondeur, D. Paradigms and Paradoxes in Modeling Adsorption and Chromatographic Separations. Ind. Eng. Chem. Res. 1995, 34, 2782. Tonkovich, A. L. Y.; Carr, R. W. Experimental Evaluation of Designs for the Simulated Countercurrent Moving Bed Separator. AIChE J. 1996, 42, 683. Unger, J.; Kolios, G.; Eigenberger, G. On the Efficient Simulation and Analysis of Regenerative Processes in Cyclic Operation. Comput. Chem. Eng. 1997, 21, S167.
Received for review July 11, 1997 Revised manuscript received December 12, 1997 Accepted December 18, 1997 IE970478K
