HUNGARIAN JOURNAL OF INDUSTRIAL CHEMISTRY VESZPREM Vol. 29. pp. 95 -104 (2001) DESIGN OF A PRESSURE SWING ADSORPTION MODULE BASED ON CARBON NANOTUBES AS ADSORBENT- A MOLECULAR MODELING APPROACH A. HEYDEN, T. DOREN, M. KOLKOWSKI and F. J. KEIL (Chair of Chemical Reaction Engineering, Technical University of Ham.burg-Harburg, Eissendorfer Str. 38, D-21073 Hamburg, GERMANY) Received: March 21,2001 The prese~t paper describes the first sequence of the Skarstrom cycle for bed 1 to analyze the breakthrough curves and study the I~rtance of. the mass ~ansfer resistance for the performance of the PSA process which separates a feed stre~ ~onsisti.ng of eqmmolar fractions o~ meth~e ~d ethane. The beds are filled with carbon nanotubes (CNT) of a certam_mner diameter (3 nm). The adsorptiOn bed 1s simulated by a one-dimensional, time-dependent dispersion model. A special feature of the present paper are the adsorption isotherms of the CHJC2H6 mixtures and the diffusivities which were calculated by molecular simulation methods. Introduction All adsorption separation processes involve two essential steps: a) adsorption, during which the preferentially adsorbed species are picked up from the feed gas, and b) desorption, during which these species are removed from the adsorbent, and is thus regenerated for the next cycle. The essential feature of the pressure swing adsorption (PSA) process is the removal of the adsorbed species by reducing the total pressure, rather than by raising the temperature or purging with a sweep ~as. The process operates under approximately Isothermal conditions. As the pressure can be changed much more rapidly than the temperature, it is possible to operate a PSA process on a much faster cycle, and thereby the throughput per unit of adsorbent bed volume is considerably increased. The major limitation is that PSA processes are restricted to components that are not too strongly . adso_rbed. Air separation. air drying, hydrogen punfication, carbon dioxide recovery and natu~al ~as purification are at present the most important applt~at10ns of the PSA process. Detailed descriptions of this process and its applications can be found in books by Ruthven et al. [1]. Yang [2], and Thomas and Crittenden [3]. An industrially employed PSA cycle is the Skarstrom cycle which in its basic form utilizes two pack~ adsorbent beds (see Fig.J). The Skarstrom cycle conststs of four steps: a) pressurization~ b) adsorption, c) countercurrent blowdown. and d) countercurrent purge. Contact information: E-mail: keil@tu-harburg.de In step 1, bed 1 is pressurized to the higher operatino pressure with feed from the feed end (open valves: 1,3; closed valves: 2,4-9) while isolated from bed 2. During this step bed 2 is blown down from high pressure to atmospheric pressure in the opposite direction. These steps are relatively fast in most industrial processes and often have no cignificant influence on the performance of the entire process. In step 2 (open valves: 1,3,5,6,8~ closed valves: 2,4,7,9), high pressure feed flows through bed 1. The product stream is separated from the stronger adsorbing species. The concentration front of the stronger adsorbing species reaches the exit of the bed at a far lower concentration compared to the weaker adsorbing component during this step. A fraction of the effluent stream is withdrawn as product, while the remaining gas is used to purge bed 2 at the low operating pressure. The direction of the purge flow is opposite to that of the feed flow. The purge to feed volume ratio must be at least one in order to sweep bed 2. Steps 3 and 4 follow the same sequence with the bed interchanged. The present paper describes lhe first sequence of the Skarstrom cycle for bed 1 to analyze the breakthrough curves and study the importance of the mass transfer resistance for the performance of the PSA process which separates a feed stream consisling of equimolar fractions of methane and ethane. The beds are filled with carbon nanotubes (CNT) of a certain inner diameter (3 nm). The adsorption bed is simulated by a one-dimensional, time~dependent dispersion model. A 96 RAFFINATE PRODUCT Bed 1 Bed 2 FEED Fig. I Schematic diagram of the basic two-bed pressure swing adsorption system. special feature of the present paper is the adsorption isotherms of the CHJC2~ mixtures and the diffusivities which were calculated by molecular simulation methods. Theory 1 Grand Canonical Monte Carlo Simulations of Adsorption Isotherms A Monte Carlo simulation (MC) generates configurations of a molecular system by making random changes to the positions of the species present (together with their orientations and conformations where appropriate). The potential energy of each configuration of the system, together with the values of other properties, can be calculated from the positions of molecules in space. The MC method thus samples from a 3N-dimensional space of the positions of N particles. In adsorption studies one would like to know the amount of material adsorbed as a function of pressure and temperature of the reservoir with which the porous material is in contact. A proper choice for adsorption studies is the grand canonical ensemble (or J.1 VT ensemble). In this ensemble, the chemical potential J.t, the volume V and the temperature T are fixed. This corresponds to an experimental setup where the adsorbed gas is in equilibrium with the gas in the reservoir. The equilibrium conditions are that the temperature and chemical potential of the gas inside and outside the adsorbent must be equal. The gas that is in contact with the adsorbent can be considered as a reservoir that imposes a temperature and chemical potential on the adsorbed gas. One has to know only the temperature and chemical potential of this reservoir to determine the equilibrium concentration inside the adsorbent. This is done with the grand-canonical ensemble. The temperature and chemical {X)tential are set. fixed and the number of particles is allowed to fluctuate during the simulation. One should keep in mind that the pressure is related to the chemical potential via an equation of state. Inside microporous materials the pressure is not defined. The partition function for the grand-canonical ensemble is given by ( 4-6]: Q(Ji.,V,T) = L exp(f3J.1N;q:VN J dsN exp[- fJU(sN)](l) N=O N. and the corresponding probability density N (sN N) oc exp(fJ!.JN)!N q: exp~ fJU(sN)] (2) ~~ ' N! with (3) q: · VN =kinetic energy contribution ofN molecules in the system. ( 4) In a grand-canonical simulation, one has to sample the distribution eq. (2) by doing the following trial moves: a) Displacement of particles. A particle is selected randomly and given a new configuration. This trial move is accepted with a probability: acc(s -7 s') = min~,exp{- fJ[U(s'N) -U(sN)n) (5) b) Insertion and removal of particles. A particle is inserted at a random position or a randomly selected particle is removed. The creation of a particle is accepted with a probability acc(N -7 N + 1)= =min[l.~·exp{B(JL-U(N +l)+U(N)]}l (6) (N +1) j and the removal of a particle is accepted with a probability acc(N -7 N -1) = = min[l.~{- {3[u+U(N -1)-U(N)]}l (?) Vqt j Details of the implementation of eqs. (5-7) may be found in [5]. The grand-canonical MC technique can be applied to spherical molecules and simple models of nonspherical molecules. A common model to represent molecules in simulation is the united atom model. Methane is represented by a single sphere whereas chain molecules are represented by several joined spheres each denoting a methyl group of the chain molecule (e.g. ethane is represented by two joined CH3-spheres). In more complicated cases biased sampling techniques have to be used [7]. The main feature of this more sophisticated MC method is that the trial moves are no longer completely random, instead the moves are biased in such a way that the molecule to be e.g. inserted has an enhanced probability to 4'fit" into the existing configuration. This approach was used for ethane in the present paper. The method to increase the speed of MC simulations for chain molecules by MC moves is called Configuration-Bias Monte Carlo (CBMC) and was developed by Siepmann and Frenkel [7]. Further details may be found in the book by Frenkel and Smit [5] and in articles by Vlugt et al. [8-10]. The main feature of this more sophisticated MC approach is that the trial moves are no longer completely random, instead they are biased in such a way that the molecule to be, for example, inserted has an enhanced probability to "fit" into the existing configuration. No information about the present configuration of the system is used in the generation of unbiased MC trial moves. This information is used only in the acceptance criteria. The CBMC approach considerably improves the conformational sampling for molecules with complicated structures. To satisfy detailed balance, one has to change the acceptance rules. In the CBMC scheme it is common to split the total potential energy of a trial site into two parts. The first part consists of an internal, bonded, intra-molecular potential uint which is used for the generation of trial orientations. The second part of the potential, the external potential uexr, is used to bias the selection of a site from the set of trial sites. It should be noted that the split of the potential into uint and uext is arbitrary and can be optimized for particular applications. For a new configuration n, a randomly chosen molecule of length N is regrown segment by segment. First, f trial sites for the first bead (the first atom or united atom) of the molecule are placed at random positions in the simulation box. In this work/ is equal to one. This value was chosen as in the present paper only small molecules (methane and ethane) are investigated. For larger molecules f can be larger than 1. The so-called Rosenbluth weight of this segment is then [5, 8-10]: f W1 (n) = L exp[- u:;r [3] (8) j=l The interaction energy u Ii includes all interactions of atom i with other atoms in the system. One of these trial sites is now selected with probability P.~ez(b-)= exp[-utd3] (9) h I ( ) W1 n For all other segments l of the molecule, k trial orientation bi are generated according to the Boltzmann weight of the internal potential of that segment P.?en (b _)db = exp[- u:t [3}ih lz I f r int/3LL expl-u1 po (10) In the present paper k is set equal to six. This value can be optimized [ 11]. The set of k different trial segments is denoted by (ll) Out of these k trial orientations one is chosen according to the Boltzmann weight of its external potential 97 (12) This procedure is repeated until the entire chain of length N has been grown. The Rosenbluth weight W(n) of the new configuration is then defined as The old configuration (a) is retraced in a similar way, except that for each segment only k-1 (f-1 for the first bead) trial orientations are generated with a probability given in eq. (10). The k-th lf-th) trial orientation is the old orientation. The Rosenbluth weight of the old configuration W(o) is then defined as w1 ( o) fi[± exp(- uif" fj ]l W(o) = l=Z J=l J (14) f·kN-1 where w 1(a) is the Rosenbluth weight of the first segment of the old configuration. To satisfy detailed balance, the new configuration is accepted with the probability acc(o -7 n) = mm 1,--. ( W(n)) W(a) (15) In the present paper the following acceptance rules for creation and deletion were employed: creation: deletion: . { N (o) 1 I acc(o-7n)=mm 1,-k_. __ Jt ftf3V W(o) (16) (17) In a trial move a molecule is displaced, and rotated in case it is not symmetrical in all directions. The distances and the angels of rotation are chosen randomly but are limited by a maximum displacement in x, y and z direction and maximum rotation angle. A good sampling of configuration space is achieved with an acceptance rate for the trial move of approximately 50%. The maximum limits of displacement are adjusted during the simulation according to this value. It should be kept in mind that a trial move in MC simulation is different from physical motion in that respect that no check for collision on the path between old and new position is done. Therefore, the maximum displacement of a molecule should always be small enough to prevent molecules from moving through the pore wall. There is no biasing-like algorithm for the trial move. Consequently, the acceptance criteria for a trial move is: 98 acc(O --7 n) z min{l,expf- {3(um (n) -um (o )) ]} (18) The new coordinates of each bead of the molecule after a translational move can be calculated as (19) In the rotational case an equivalent equation looks like able to change the internal orientation of the beads of one molecule the partial trial regrowth is used. Since the first bead is the same in the old and the new configuration the Rosenbluth weight of the first segment is in the old and new configuration the same and can be set equal to one: (27) The rotational matrix A is given by: (20) With this change the same acceptance rule as for a total regrowth trial is employed: The matrix A is described by a quaternion of unit norm. Such a quaternion may be thought of as a unit vector in a four~dimensional space: Qz(qo,ql,q2,q3) with q~+q~+qi+qi=l (22) · If the reference frame in the center of mass of the molecule is placed in such a way that the Eulerian angles are all zero for the orientation of the old configuration, then the quaternion is: Q = (1,0,0,0). To determine the quaternion of the new orientation and therewith the rotational matrix the following scheme is used: 1. Generate two random numbers , 1, ' 2 between -1 and + 1 which satisfy the condition: (23) 2. Generate a second set of two random numbers {3, G between -1 and + 1 which satisfy the same condition: S2 ={i +si ~I (24) 3. The quaternion for the new orientation is: (25) 4. To satisfy the condition that 5Q% of the rotational trials are accepted, the following constraint applies: (26) whereby "value" must be adjusted such that the 50%~rule is fulfilled [5]. For molecules which consist of more than one united·atom, trials that change the internal orientation of the beads of one molecule are necessary to speed-up the simulation. There are two trials that change this internal orientation, a total regrowth trial in which the coordinates of an beads of the molecules change~ and a partial regrowth trial where aU beads except the first one change their coordinates. For a total regrowth the conventional CBMC scheme described above was used. Jn dense systems the insertion of the first bead in a total regrowth trial is unlikely. consequently the success of the entire trial becomes improbable. To be nevertheless acc(o~n)=mm 1,--. ( W(n)) W(o) (28) In dense systems the creation and deletion events occur, despite the use of CBMC, only with a very low probability. Therefore, so-called swap trials were introduced where the identity of a molecule changes. This trial is equivalent of deleting a molecule of species i and replacing it by a molecule j. The following acceptance rule for the swap trials was employed: acc(o-7n)==min[l, N;(o) -- 1 --W(n)· 11 J (29) ,h W(o) N 1(o)+l After 150,000 MC trials equilibrium was reached at all pressures. Actually, we used 250,000 equilibration steps. By doing 20 simulations of the same system, one observes that 37,500 MC steps for each pressure are enough to sample the configurational space. 2 Molecular Dynamics In order to calculate the diffusion coefficients molecular dynamics (MD) calculations were done. The MD approach is a solution of the N-body problem of classical mechanics by integrating Newton's equations of motion for each atom k: (30) After the starting configuration is selected, the force on each united atom is calculated and the equations of motion are integrated. The integration is done by means of the Verlet-algorithm [5,6]. First, 250,000 equilibration steps are calculated and then the time steps follow to sample equilibrium to get the desired self- diffusion coefficients by employing the Einstein- relation: (31) with (32) The transition from self-diffusion coefficients, detennined by means of MD calculations, to transport diffusion coefficients was done by a Darken factor: D (alnfk) Dr,k = s,k dlnp k pi (33) To conclude, first MC simulations in the grand- canonical ensemble were done. Hereby, a specific equilibrium concentration of each component in the pore is determined at certain conditions in the bulk phase. Then a molecular dynamics calculation at the equilibrium pore concentrations was done to calculate the self-diffusion coefficients. The Darken factor can be obtain~d from equilibrium adsorption isotherm data. Therefore, it is necessary to calculate the derivative of the component fugacity with respect to the amount adsorbed of one component by leaving constant the amount adsorbed of the other components. In practice this is not possible since the input data of our simulations are total bulk pressure and gas fraction composition and it is hardly possible to obtain enough data at constant amount adsorbed of one component to calculate an accurate derivative. As a result Jacobian Transformations were used [12]. With the help of Jacobian Transformations, eq. (34), the derivative at constant amount adsorbed of one component can be replaced by multiple derivatives with respect to the input data. This facilitates greatly the calculation of accurate derivatives. a(lnfk ,In p 1 ) aijnpk ,lnpJ= (34) In case of non-polar molecules like methane and ethane, only external interactions have to be considered. In the present paper the Lennard-Jones (L-J) potential was employed. The interaction between two Lennard-Jones sites i andj is given by the potential [5,6]: (35) The parameters aii and eii are obtained from the Lorentz- Berthelot mixing rules: (36) {37) The L-J potential is an effective pairwise additive potential. The interaction energy between every L-J site with all other L-J sites must be summed up: 99 U(r)= LLUij(rii) (38) irc (39) Periodic boundary conditions and the minimum image convention [5,6] were employed along the pore axis. 3 Model of the Adsorber Bed For the PSA module a one-dimensional adsorber model was used. The mass balance for each component can be written at all operating conditions of the Skarstrom cycle as: dck ={ dck du]_l-e . dqk +D d2ck (40) u ':\ +ck ':\ Ps ':\ ax ':\ 2 , dt az az e at uz Assuming that the axial pressure drop is negligible, the sum of the gas concentrations of an components is constant and no function of the axial coordinate: (41) Under this condition, the overall mass balance during the high-pressure adsorption and low-pressure purge steps, where the total pressure is approximately constant, looks like: 3u 1-e "'dq1 C·-+-·Ps,L.,;-==0 oz e 1 dt (42) Augmented by an appropriate mass transfer model, eqs. (40-42) form a set of equations· which describe the time evolution of the concentration front of each component. For the mass transfer two models were used. First, for comparison, the rather unrealistic equilibrium model, which assumes that the average amount adsorbed per unit volume pellet is equal to the equilibrium amount adsorbed at the concentration of the surrounding gas phase: 3if1 aq; aq; ac;; a;=a;= dc 1 ·a; (43) Second, a linear driving force model was employed. which gives a quite realistic description of the adsorption kinetics: 100 Table 1 Operating conditions and design parameters of the PSA process simulations. Molar ratio of CHJC2Ilt; in feed Feed mole flow Adsorption pressure Desorption pressure Temperature Purge to feed volume ratio Bed length Bed diameter Bed porosity Duration of high pressure feed Duration of low pressure purge Representative pore diameter of the pellets Pellet density Pellet diameter (lqk f* -) a;-=kk\qk -qk where the Glueckauf constant kk [13] 15Dek kk=--'- R2 p 50150 4mol!s lObar 1bar 300K 2 2m 0.4m 0.4 140 s 140 s ·3nm 1500kg!m3 0.004m (44) (45) was introduced. The equilibrium adsorption data, q;, were obtained from configurational bias MC calculations and the diffusivities from MD and MC calculations (see below). The axial dispersion coefficient, Dax• was estimated by a correlation from Wen and Fan [14]: 1 Peax,p 0.3 (46} for 0.008 ~ 1.2 .g., '0 .8 0.8 l5 ~ 0.6 8 Bulk mole fraction of Ett1ane _j L------------- Fig.4 Like Fig.3 for 1 h