Optimization of a Pressure Swing Adsorption Process for Nitrogen

Article pubs.acs.org/IECR

Optimization of a Pressure Swing Adsorption Process for Nitrogen Rejection from Natural Gas Surya Effendy, Chen Xu, and Shamsuzzaman Farooq* Department of Chemical & Biomolecular Engineering, National University of Singapore, 4 Engineering Drive, 117585 Singapore S Supporting Information *

ABSTRACT: Nitrogen rejection from natural gas is a capital intensive process. The most widely used technology, cryogenic distillation, is uneconomical for gas wells in remote areas or with small production rate. In this work, a detailed simulation and optimization is carried out to explore the feasibility of pressure-vacuum swing adsorption (PVSA) technology for onsite nitrogen rejection at remote/low-throughput gas wells. The designed PVSA unit possesses high mobility at a cost that scales well with low throughput. A nonisothermal and nonisobaric simulation model including detailed description of both micropore and barrier resistances in a kinetically selective carbon molecular sieve adsorbent has been used for process optimization. The process parameters constitute the decision variables of the optimization problem, with recovery and productivity as the objective functions, and industrial pipeline purity specification as an inequality constraint. Starting with a 6-step PVSA cycle, it is shown that a modified 9-step PVSA cycle promises to improve further the process performance. The process simulation model is the most comprehensive model for kinetically controlled PSA processes to date and the optimum results exceed the performances expected from the trends revealed in parametric studies published so far. Preliminary costing of the optimized process reveals that the modified PSA cycle has the potential to produce the desired product at a lower cost than the competing alternatives. (such as MOLECULAR GATE and Nitrogen Sponge).5 Among these technologies, PSA separation has attracted much attention due to its insensitivity to low production rate. PSA has also been studied in other CH4 enrichment applications from CH4/N2 mixture, such as coal-bed gas,6 coal-mine ventilation gas and LNG vent gas.7 Here, we report the findings from a detailed simulation and optimization study of the PSA technology using a Takeda carbon molecular sieve (CMS) for upgrading N2-contaminated natural gas to the required pipeline specification of CH4 purity for transmission to the consumers. For a PSA-based natural gas upgrading process to be competitive, the adsorbent should be N2-selective, so as to obtain CH4 as the high-pressure raffinate product. This reduces the work needed to pressurize CH4 back to pipeline pressure. Among the known commercial adsorbents, CMS possesses the requisite kinetic selectivity for N2 over CH4 to meet the process requirement. Since the kinetic diameter of N2 is smaller than CH4, CH4 can be produced as the high-pressure raffinate product. Several new adsorbents in the literature likewise rely on the size difference between N2 and CH4 to create kinetic selectivity.8−11 However, despite the promising kinetic selectivity of various materials, very few process studies were able to produce CH4 product stream satisfying both the strict

1. INTRODUCTION Natural gas is the cleanest fossil fuel available. Its major component is methane (CH4), which produces only water vapor and carbon dioxide (CO2) upon combustion. It generates less atmospheric pollutants, e.g., nitrogen and sulfur oxides, and green-house emission per unit amount of extractable energy. The production and consumption of this clean energy source has been increasing at a rate of 2.5% and 2.4% per annum respectively over the last ten years.1 This motivates exploration and utilization of gas wells that were considered uneconomical in the past due to their high contamination level or low throughput. The major contaminants in natural gas are higher hydrocarbons, moisture, nitrogen (N2), CO2, hydrogen sulfide (H2S) and other sulfur components.2 N2 rejection is particularly challenging due to its close chemical and physical properties with CH4. Cryogenic distillation has been utilized as the N2 removal process in most industrial natural gas treatment plants. However, the cost for nitrogen removal by cryogenic distillation per unit volume of natural gas increases rapidly with decrease in gas well capacity.3 This imposes a significant penalty on the overall treatment cost of natural gas from wells located in remote areas, or with low production rate. Consequently, alternative technologies capable of reducing N2 level in raw natural gas (some as high as 20 mol %) to a pipeline level of less than 5 mol % for such wells are being proposed. These technologies include membrane separation,4 solvent or chelating chemical absorption, cryogenic absorption using lean oil and pressure swing adsorption by molecular sieves © XXXX American Chemical Society

Received: February 6, 2017 Revised: April 16, 2017 Accepted: April 17, 2017

A

DOI: 10.1021/acs.iecr.7b00513 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

between the bulk phase and macropore gas offered better agreement with the experimental results than the simplified bilinear driving force (bi-LDF) model used by Grande and Rodrigues16 for the same separation. Here, we extend the PSA simulation model of Bhadra and Farooq12 for kinetic separation of a binary mixture in CMS adsorbents to make it nonisothermal and nonisobaric. These two features are important for large scale industrial process design. The simulation model is equally applicable to adsorbents that do not show barrier resistance at the micropore mouth or the stronger-than-expected concentration dependence. This is done by appropriately assigning zero/large values to the relevant parameters. 2.1. Assumptions. The feed is assumed to be an ideal gas mixture of CH4 and N2 with known composition. The feed composition will be varied in order to find the threshold N2 fraction up to which CH4/N2 separation by PSA could meet the industrial requirements. Several major assumptions are made during the construction of the model. The fluid flow through the packed column is represented by an axially dispersed plug flow model. All the properties of the adsorbent and fluid in the packed column are radially isotropic. Prandtl number and heat capacity ratio are assumed to be constant. Fluid viscosity of the gas phase varies with temperature, but not composition, while heat capacity varies with composition, but not temperature. The solid phase is in thermal equilibrium with the surrounding fluid throughout the PSA process. The drag coefficient in the external piping is assumed to be constant. The performance of the PSA nitrogen rejection unit can be adequately captured by monitoring the product purity, product recovery and productivity at cyclic steady-state. 2.2. PSA Cycle. We begin with a 2-column, 6-step cycle comprising feed pressurization (FP), high-pressure adsorption (HPA), depressurizing pressure equalization (DPE), countercurrent blowdown (CnBn) to atmospheric pressure, countercurrent evacuation (CnEv), and pressurizing pressure equalization (PPE). The cycle is schematically shown in Figure 1. During DPE/PPE, the product ends of the two columns are connected, causing material flow from the high- to the lowpressure column. The columns eventually reach some

pipeline composition and industrial CH4 recovery requirements.12 One possible reason for this is that all the adsorption cycles analyzed in these work were far from optimal in both process parameters and operating conditions. In this work, we focus on the feasibility of N2 rejection using PSA on CMS by optimizing the operating process parameters. We begin with simulating and optimizing a standard 6-step PSA cycle using commercial Takeda CMS for N2/CH4 separation. Both CH4 recovery and adsorbent productivity are optimized for given CH4 purity requirements. Detailed simulation and optimization methodology are described along with the new features of the simulation model that distinguishes it as the most comprehensive model for kinetically controlled PSA processes to date. The impact of CH4 purity requirement and decision variable values on column performance is discussed. An alternative PSA cycle is proposed and its performance is compared with the 6-step cycle. A cost analysis of the proposed PSA process is also presented and compared with the competing alternative technologies.

2. FEATURES OF THE SIMULATION MODEL AND SUPPORTING EXPERIMENTAL EVIDENCE The CMS adsorbent particles have bidispersed pore structure comprising macropores and micropores. Transport of gases in macropores and micropores of CMS has been discussed and modeled in several papers. Qinglin et al.13 proposed and verified a dual-resistance model to describe transport behavior of four gases (N2, O2, CH4, and CO2) in three CMS samples (one from Bergbau-Forschung and two from Takeda). The single component uptake in the micropores was observed to be controlled by a barrier resistance at the pore mouth followed by a distributed pore diffusional resistance inside the micropores. Transport parameters were obtained by fitting the model to the experimental uptake data. These parameters were observed to have much stronger dependence on adsorbent loading than that expected from chemical potential gradient as the driving force for diffusion. An empirical correlation between transport parameters and adsorbent loading was proposed to account for the stronger concentration dependence and validated with experimental data. Qinglin et al.14 later extended this model to the adsorption kinetics of binary and ternary gas mixtures of N2, O2, CH4, and CO2, and validated the extended model with experimental integral uptake data on two of the adsorbents. The multisite Langmuir isotherm model was found to effectively describe unary equilibrium and predict the mixture behavior in both CMS samples. Bhadra and Farooq12 incorporated this detailed dual resistance model developed for the micropores in CMS adsorbents (while also accounting for the presence of the macropores) into an isothermal and isobaric PSA process simulation model, and applied it for CH4−N2 adsorption in two CMS adsorbents. Although Takeda CMS was able to achieve a methane purity of 96%, the recovery was only 70%, which was too low for industrial application. The PSA column performance was found to be sensitive to the process parameters, thus highlighting the importance of constrained process optimization in order to find the best operating conditions while meeting the process requirements. A nonisothermal, isobaric PSA simulation model using a concentration dependent micropore diffusion model with chemical potential gradient as the driving force was simulated and validated by Khalighi et al. (2011) for propylene/propane separation on 4A zeolite.15 The concentration dependent micropore diffusion model and equilibrium assumption

Figure 1. 6-step cycle comprising feed pressurization (FP), highpressure adsorption (HPA), depressurizing pressure equalization (DPE), counter-current blowdown to atmospheric pressure (CnBn), counter-current evacuation (CnEv), and pressurizing pressure equalization (PPE). B


Article


Here qi̅ is the adsorbed amount of component i, normalized against its saturation loading. The terms Qi,1 and Q2b have the following forms:

intermediate pressure. Because DPE and PPE steps are coupled, the durations of DPE and PPE must be equal. Moreover, the durations of the FP and HPA steps must be equal to the durations of the CnBn and CnEv steps in order to synchronize the occurrences of the DPE and PPE steps. Because the durations of FP, HPA, CnBn, and CnEv are not constrained by the equality requirement for a given set of decision variables, it may be necessary to insert an idle (IDL) step in the cycle to fulfill this requirement. In this study, we have opted to insert IDL between CnEv and PPE or PPE and FP. 2.3. Model Equations. The simulated PSA column consists of four interlinked subsystems: external fluid-phase, column wall, macropores, and micropores. Each subsystem consists of a set of nondimensionalized mass and energy balances, and they are linked to each other via correlations or constitutive equations. For brevity, only the macropore and micropore balances are shown in the present work, since they contain features which have not been reported in the literature. The external fluid-phase and wall balances are presented in the Supporting Information. 2.3.1. Macropore Balances. Transport in the macropores is assumed to be molecular diffusion-controlled: ∂yp, i ∂t

yp, i ∂P

+

=−

P ∂t 1 − εp

(qp, j

εp

−

⎡ 2 ⎤ 1 Dm, i ⎢ ∂ yp, i 2 ∂yp, i ⎥ − qj|r = rc ) + + τ R p2 ⎢⎣ ∂R2 R ∂R ⎥⎦

P

r = rc+

⎛ Ep, i ⎞ ∂ ⎛ 2 ∂qi̅ exp + r 2Q 2p(Q i ,1 − 1) ⎜− ⎟ ⎜⎜r Q 2p 2 ∂r ⎝ RT ⎠ ∂r ⎝ r ∂qj̅ ⎞ ⎟ ∂r ⎟⎠

(6)

The right-hand side of eq 6 describes the micropore flux. Note that

(1)

Q 2p = 1 +

∑i βp, iqi̅ 1 − ∑i qi̅

(7)

where eq 7 is a semiempirical correction to the crystal diffusivity, analogous to eq 4.13 2.4. Initial and Boundary Conditions. The columns are initialized with feed at room temperature and pressure, and the solid-phase at equilibrium with the fluid-phase. The boundary conditions for the external fluid-phase balances vary according to the step being simulated. The complete set of boundary conditions is provided in the Supporting Information. The boundary conditions for eqs 1 and 6 are obtained by equating fluxes at the surface of the relevant particles, and by assuming smooth symmetric profiles at the center of spherical particles. The boundary conditions for macropore transport are

and temperature, respectively. Note that, in general, ∑i yp,i is not equal to unity. In effect, macropore concentration is being nondimensionalized using external fluid-phase pressure and temperature. Note that the two terms on the right-hand side of eq 1 correspond to micropore uptake and macropore flux, respectively. The term

(5)

′ i Dc0,

j

where P and T are the external fluid-phase pressure

qi ∂ciim ciim ∂qj

=

∑

external fluid-phase balances directly (see Section S1 in the Supporting Information). We thus introduce a new variable C p, iRT

(4)

2.3.2. Micropore Balances. Transport in micropores is given

The terms on the left-hand side of eq 1 possess an unusual definition. First, note that the macropore need not necessarily have the same pressure as the external fluid-phase. Because transport in the macropore is driven by molecular diffusion, there is no reason to model macropore pressure as a separate intrinsic variable. Consequently, the appropriate variable to model should be the macropore concentration (Cp,i). However, it is convenient to be able to calculate (yi − yp, i |R̅ = 1 ) in the

yp, i ≡

1 − ∑i qi̅

by

∂t

j

∑i βb, iqi̅

qp,̅ i ⎛ ΔU ⎞ Pyp, i = b0, i exp⎜ − i ⎟ ⎝ RT ⎠ RT (1 − ∑i qp,̅ i)ai

T ∂t

∑

(3)

where eq 4 is a semiempirical correction to the barrier resistance based on local crystal loading. A more detailed discussion on this semiempirical correction may be found in a previous work by our research group.13 Equation 1 requires evaluation of qp,i, which is obtained from the multisite Langmuir isotherm:

∂qi̅

⎛ E b, i ⎞ RT ′ i exp⎜ − k b0, ⎟Q | ⎝ RT ⎠ 2b r = rc P

1 − ∑i qi̅

Q 2b = 1 +

yp, i ∂T qi ∂ciim ciim ∂qj

aiqi̅

Q i ,1 = 1 +

∂yp, i ∂R

= R=R p

k f, i τ (y − yp, i |R = R p ) Dm, i εp i

(8)

is the magnitude of the

∂yp, i

contribution of the adsorbed-phase concentration of component j on the chemical potential gradient of component i. This term may be evaluated explicitly for a multisite Langmuir isotherm: ⎧ i = j ⇒ Q i ,1 ⎫ qi̅ ∂ciim ⎪ ⎪ ⎨ ⎬ = im ⎪ ⎪ ci ∂qj̅ ⎩ i ≠ j ⇒ Q i ,1 − 1⎭

∂R̅

=0 R=0

(9)

The boundary conditions for micropore transport follow the same form, although there are additional complexities due to chemical potential gradient and pore size distribution. As mentioned in the previous section, these are accounted for via the terms Qi,1, Q2b, and Q2p:

(2) C


Article

Industrial & Engineering Chemistry Research ∂qi ∂r

+ (Q i ,1|r = rc − 1)∑

=

∂r

r = rc

⎛ E b, i − Ep, i ⎞ Q 2b ′ i rc k b0, exp⎜ − ⎟ ′ i ⎝ RT 3 Dc0, ⎠ Q 2p

∑ j

∂qi ∂r

j

r = rc

qi ∂ciim ciim ∂qj

composite error of all numerical methods, i.e., finite volume method, orthogonal collocation, ode15s, CSS criterion and other nonlinear equation solvers is estimated by decreasing the number of cells in the external fluid-phase to 15, decreasing the number of macropore and micropore nodes to 1 and 3, and doubling all tolerances for the aforementioned cases. The composite error is found to be acceptable at ca. 1% for product purity, 0.5% for recovery, and 10−5 molCH4/(kgads·s) for productivity. The simulations were run on Dell Precision Tower 7910 XCTO with Intel Xeon processor E5-2699 v3 (18C, 2.3 GHz, Turbo HT45M) and 32 GB (4 × 8GB) 2133 MHz DDR4 running on Windows 10. Following the algorithm outlined above, a vast majority of cases reach cyclic steady-state within 90 cycles. The average simulation time required to reach cyclic steady-state was 15 min, and almost never exceeded 45 min. 2.6. Model Validation. The proposed model is validated against the experimental data published by Qinglin. 20 Validation is done for runs 9 through 20 in Table 6.2 of the aforementioned work. In the interest of brevity, two of the validations are shown in Figure 2. The remainder of the validation may be found in the Supporting Information. A comparison between the present model and that reported by Bhadra and Farooq,12 together with a comprehensive case study showing the non-negligible impact of the isothermal and isobaric assumption, can also be found in the Supporting Information. 2.7. Optimization Method. Product recovery and productivity are selected as objective functions for optimization, subject to an inequality constraint on product purity. The inequality constraint is derived from the energy content requirement for pipeline natural gas,21 and is found to be between 93 to 97% depending on prevailing legislations. The decision variables and their upper and lower bounds during optimization are listed in Table 1. Energy consumption is not considered as an objective function since the energy cost is negligible compared to capital and other operating costs. Nondominated Sorting Genetic Algorithm-II (NSGA-II), available from MATLAB Central, was used for optimization. Deb et al.22 recommended the number of candidate solutions to be 10 × the number of decision variables. There are 5 decision variables, leading to 50 candidate solutions. However, a conservative value of 60 candidate solutions was used instead. Optimization was found to be complete within 40 generations, but was run for a further 20 generations. Using 12 parallel processors, the optimization was generally complete within 3 to 4 days on the aforementioned workstation.

∂qj

r = rc

(qp, j − qj|r = rc ) r = rc+

(10)

=0 (11)

r = rc

The interaction terms for barrier resistance on the right-hand side of eq 10 was neglected in the earlier work published from this laboratory.17 The contributions of the interaction terms were found to be negligible in this study. eq 10 consumed the bulk of computational time, and may lead to a prohibitively long simulation depending on the choice of the numerical method. The derivation of eq 10, as well as various simplifications and solution algorithms to reduce the computational time, are discussed in the Supporting Information. 2.5. Numerical Method. The model equations are solved using finite volume method with forward 3-point weighted essentially nonoscillating (WENO) flux limiter in the external fluid-phase, and the method of orthogonal collocation in the macropore and micropore subsystems. The external fluid-phase is discretized into 30 cells following the recommendation made in a previous work.18 The macropore concentration is found to be essentially independent of the radial location at all time steps in all simulations. This is to be expected, since mass transfer resistances of CH4 and N2 are strongly dominated by micropore diffusion. Consequently, the appropriate number of macropore nodes should be 1. However, we have opted for a conservative approach and have chosen 2 internal nodes instead. Micropore concentration is found to oscillate in the radial direction when more than 3 internal nodes are used. The degree of oscillation increases with increased number of internal nodes. Trial-and-error suggests that 4 internal nodes is able to provide a robust simulation capable of capturing concentration profile while avoiding unphysical oscillation. Orthogonal collocation is done using even series of Legendre polynomials. The internal points of the collocation scheme is selected to be the roots of the nth order (where n = N2 or N3) weighted Jacobi polynomial.19 The spatially discretized set of equations is then solved using MATLAB 2011b in-built function ode15s with the Vectorized and JPattern options. The simulation is run with the final state of one step being used as the initial state of the next step until cyclic steady-state is achieved. The cyclic steady-state is detected using the relative L1-norm error: E ≡ ||

f p+1 − f p δ + f p+1

3. RESULTS AND DISCUSSION 3.1. CH4 Purity, Recovery, and Productivity. The effect of varying N2 contamination in the feed on the recoveryproductivity Pareto front for a fixed product purity of 95% is shown in Figure 3. As expected, there is a trade-off between product recovery and productivity, and increasing N 2 concentration causes the Pareto curves to move in an unfavorable direction. The maximum achievable product recovery at 10%, 15%, and 20% N2 contamination in feed are 97.9%, 97.3%, and 95.8% respectively, while the maximum productivities are 0.00546, 0.00286, and 0.00182 molCH4/(kgads· s), respectively. The performance of the Molecular Gate SPEC PLANT for a feed containing 10% N2 and delivering CH4 product at a purity of 96% is included in Figure 3. The data is estimated using the

|| 1

(12)

Here fp is a vector containing all simulation variables at cycle “p” and δ is a small number, set to be 10−4. Let N4 be the cycle number at which E first decreases below 0.01. The cyclic steady-state is assumed to have been reached at 3 × max{N4, 30} cycles. This algorithm is tested against several randomly selected cases, and is found to be very conservative. The D


Article


Table 1. Lower and Upper Bound of Decision Variables Name of Decision Variables (units)

Lower Bound

Upper Bound

Adsorption pressure (bar) Evacuation pressure (bar) Adsorption velocity (m/s) Evacuation velocity (m/s) Adsorption duration (s)

3 0.1 0.005 0.02 800

16 0.8 0.055 0.2 5400

Figure 3. Productivity-recovery Pareto fronts at various levels of N2 contamination. The Molecular Gate process result is for 10% N2 in the feed inlet yielding 96% CH4 purity.

with a pore size of 3.7 Å. The recovery-productivity Pareto for the 6-step cycle at 96% CH4 purity is also included in Figure 3 for direct comparison. A 4-fold increase in productivity over Molecular Gate SPEC PLANT at equal product purity and recovery is achievable using the 6-step cycle and Takeda CMS. It should be noted that the column used in the Molecular Gate SPEC PLANT was 2.4 m long, versus the 1 m long column used to generate the Pareto curves in Figure 3. A separate optimization was run using a 2.4 m long column. The resultant Pareto front at 10% N2 contamination in feed and 96% product purity was found to shift slightly in an unfavorable direction. The implication of this result is discussed further in Section 3.4. Figure 4 shows the scatter plots of product recovery and productivity against selected decision variable values corresponding to the Pareto points in Figure 3. The decision variables, except the adsorption pressure, are not hitting the assigned upper and lower bounds. Based on the trend of recovery and productivity with respect to adsorption pressure, it is reasonable to deduce that there would be a marginal increase in recovery by allowing adsorption at atmospheric pressure. On the other hand, it is possible to double or even triple productivity by performing adsorption at wellhead pressure. For a Pareto front maximizing two objectives, the correlation between the first objective and any decision variable must be opposite to that of the second objective and the same decision variable. The strength of the correlation indicates the degree to which optimization is complete, as well as the sensitivity of the objective functions with respect to the decision variable. The strengths of the correlation for the scatter plots of adsorption pressure, velocity and duration suggest that the optimization is essentially complete. On the other hand, recovery and productivity do not show strong correlation with evacuation pressure, preferring to cluster around 28 000, 33 000, and 50 000 Pa for 10%, 15%, and 20% N2 contamination in feed, respectively. This implies that, within the domain of optimization for evacuation pressure, for each level of N2 contamination, there is a small region where both recovery

Figure 2. Validation for two cases whereby the column is initially saturated with 22% O2 and 78% N2. (a, b) Column is purged with pure O2. (c, d) Column is purged with pure N2. Note that the data corresponds to run 13 of the work.20

correlation provided on the company Web site.23 Note that the Molecular Gate SPEC PLANT is designed for operation in a wide range of gas fields, and may not be optimized for the specified feed conditions. The Molecular Gate SPEC PLANT, a mobile commercial PSA unit, is used to enrich CH4 from natural gas, biogas, landfill gas and coal mine gas. This technology utilizes a patented titanium silicate molecular sieve E


Article


incomplete without considering the remaining decision variables. For example, it is intuitively expected that recovery increases with adsorption velocity. This is not what is observed in Figure 4. This arises from the fact that adsorption velocity is negatively correlated with adsorption duration, and as adsorption duration decreases, recovery decreases. Increase in recovery due to adsorption velocity is therefore offset by a decrease in adsorption duration. Figure 5 shows the change in the duration of various steps (FP+HPA and CnEv) along the 15% N2 Pareto front in Figure

Figure 5. Plot of productivity and recovery against the duration of various steps (open symbols represent recovery and filled symbols represent productivity).

4 (i.e., the diamond symbols ◊⧫). The total duration of pressurization and adsorption steps is initially lower than that of blowdown and evacuation steps at 97% recovery, but eventually becomes higher at 91% recovery. The durations are equal at a product recovery of approximately 92.7%. Following the definition in Section 2.2, the idle time vanishes at a product recovery of 92.7%. This suggests that, by attempting to maximize productivity, the optimizer minimizes idle time. Although idle steps are deliberately inserted in locations where they are unlikely to harm product purity or recovery, optimization disfavors the idle step in general since they are unproductive. This validates a rule of thumb of minimizing idle steps, frequently used to design PSA cycles. The total duration of the FP and HPA steps fall between the time scales of for CH4 uptake (37 000 s) and N2 uptake (140 s), the time scales being defined as ti ≡ Figure 4. Scatter plots of product recovery and productivity against (a) adsorption pressure, (b) adsorption velocity, (c) adsorption duration, and (d) evacuation pressure (square symbols □■ for 10% N2, diamond symbols ◊⧫ for 15% N2, triangle symbols Δ▲ for 20% N2, the open symbols represent recovery and the filled symbols represent productivity).

rc 2

(

E

′ i exp − p,i 15Dc0, RT

amb

)

(13)

Note that the time scale defined in eq 13 is equivalent to the reciprocal of the constant of the commonly used linear driving force approximation,24 evaluated at the feed temperature and at zero loading. This appears to be a valuable rule of thumb for choosing an appropriate HPA duration for kinetically controlled adsorptive separation. 3.2. An Improved PSA cycle. To improve further the PSA separation performance, a new 9-step cycle, schematically shown in Figure 6, was developed and analyzed. Like the 6-step cycle, the 9-step cycle contains feed pressurization (FP) and high-pressure adsorption (HPA) steps, followed by depressurizing pressure equalization (DPE). The regeneration steps in the 6-step cycle, consisting of counter-current blowdown (CnBn) and counter-current evacuation (CnEv), have been divided into four smaller steps: cocurrent blowdown (CoBn), CnBn, cocurrent evacuation (CoEv), and CnEv. The improve-

and productivity are maximized. This region corresponds to a balance between excessively low evacuation pressures, where slow-desorbing CH4 is lost in the blowdown gas, and excessively high evacuation pressures, where N2 does not desorb, and thus the bed is not properly regenerated. It is worth noting that in Figure 4 the variations in objective function values are caused by simultaneous change in all five decision variables. Any attempt to explain changes in an objective function value due to one decision variable is F


Article


Figure 6. Proposed 9-step PSA cycle.

The optimization results for the 6-step and 9-step cycles are compared in Figure 7. It is obvious that the 9-step cycle

ment in the 9-step cycle is achieved by recycling the streams collected from the CoBn and CoEv steps, and using them to pressurize the other column undergoing a recycle pressurization (RP) step. The regenerated column then undergoes PPE, and receives RP stream from another column. The cycle then repeats itself. The operation of the CoBn and CoEv steps involve some control measures. The CoBn step is terminated when the effluent meets a certain composition criterion, or when the column reaches atmospheric pressure. If the column remains above atmospheric pressure after the CoBn step, CnBn and CnEv steps are carried out to regenerate the column. If the CoBn step is terminated because the column reaches atmospheric pressure, additional CH4 is recovered at subatmospheric pressure in a similar fashion via the CoEv step. The column is then regenerated via the CnEv step. The control logic implies that the CnBn and CoEv steps are mutually exclusive−the former is immediately followed by subatmospheric regeneration (CnEv), whereas the latter cannot be preceded by the CnBn step because the product end of the column would be exhausted of recycle stream of desired composition. The 9-step cycle is expected to improve CH4 recovery at the expense of additional energy consumption. The CH4 recovered from the CoBn and CoEv steps must be recompressed to high pressure, which increases energy consumption. Nevertheless, the annualized cost of electricity required to run the 9-step cycle remains negligible (approximately 7%) compared to capital and other operating costs, and therefore the recovery advantage is expected to dominate. Note that because the DPE and PPE steps are coupled, the durations of these steps must be equal, and the duration in between these steps must also be equal. An IDL step is inserted in between the CnEv and PPE steps or the RP and FP steps to enforce the latter equality. Six decision variables are used to optimize the 9-step cycle. Five of these decision variables are identical to those used to optimize the 6-step cycle (see Table 1). The sixth decision variable is the product-end composition at which the CoBn and CoEv steps are terminated, henceforth referred to as composition cutoff. The lower bound for composition cutoff is set at −15% of CH4 composition in feed, while the upper bound is set as pure methane. Note that the 9-step cycle reduces to the 6-step cycle as the composition cutoff approaches 100% methane.

Figure 7. Performance improvement of the 9-step cycle over the 6step cycle.

outperforms the 6-step cycle, especially with respect to the maximum achievable recovery. At all levels of N2 contamination in the feed, the maximum CH4 recovery reaches up to 99%. However, it should be noted that excessively high recovery will result in an incombustible blowdown gas upon mixing with air. This is detailed in a later section. In general, the observed trends in scatter plots of product recovery and productivity against decision variable values corresponding to the Pareto points in Figure 7 for the 9-step cycle are identical to those in Figure 4. For a given recovery, adsorption pressure, adsorption velocity and evacuation velocity of the 9-step cycle are higher than the values in the 6-step cycle, whereas the adsorption step duration decreases sharply. Figure 8 shows the sixth variable: the composition cutoff along the Pareto front. As intuitively expected, lower composition cutoff results in more CH4 recovered, but at the expense of reduced productivity. A closer examination reveals that the increased recovery arises due to additional recycle during the RP step, whereas the reduced productivity arises due to longer RP step and shorter HPA duration. Note that the RP step is slow relative to the FP step because the pressure ratio between the RP and CoEv steps varies between 16.8 and 49.2 along the Pareto front. As shown in Figure 9, the product of adsorption pressure, adsorption step duration and adsorption velocity along the Pareto fronts overlap for the 6- and 9-step cycles. Because this G


Article


containing methane), which is not achievable in this study. We apply the lower limit to determine the maximum allowable methane recovery which will keep the blowdown gas combustible upon mixing with air (for torching, which is the simplest legal disposal method for the blowdown gas). A more complete discussion on maximum allowable recovery can be found in the Supporting Information. The relationship between maximum allowable recovery, degree of N2 contamination in the feed and product purity is shown in Figure 10.

Figure 8. Composition cutoff along the Pareto front.

Figure 10. Upper limit of methane recover at different N 2 contamination levels and two different CH4 product purities above which the blowdown gas will become incombustible after mixing with air. Figure 9. Scatter plot of PHPA × vHPA × tHPA against recovery at the Pareto front.

3.4. Costing of Proposed PVSA Unit. N2 rejection via cryogenic distillation requires a high feed throughput to be economical, and is therefore unsuitable for remote gas reservoirs with limited production rate. A few alternative processes have been commercialized for use as mobile N2 rejection units. For example, the Molecular Gate SPEC PLANT is an adsorption process that can be mounted on a truck, capable of processing feed with flow rate up to 10 MMscfd (Million Standard Cubic Feet per Day).23 The total capital and operating cost of the Molecular Gate SPEC PLANT decreases from $0.50/Mscf feed for a plant processing 2 MMscfd feed to $0.35/Mscf feed when the plant capacity increases to 10 MMscfd. Lokhandwala et al. (2010) reported the cost for membrane separation of N2 from CH4 as $0.46/Mscf product when the inlet N2 concentration is 15% and the feed flow rate is 10 MMscfd, with a corresponding CH4 recovery of 86%.26 The overall cost of nitrogen rejection by PVSA process is estimated in this study and compared with these two commercial alternatives. 3.4.1. Overview of Costing. The cost of the proposed PVSA unit may be broken down into six parts. The first part is the costs associated with the purchase, installation, and maintenance of the adsorption columns. The column length of the Molecular Gate SPEC PLANT is 2.4 m. Consequently, to ensure parity, the Pareto fronts reported for the proposed PVSA unit (see Figure 7) are recalculated for a column length of 2.4 m. Because the column cost benefits from economies of scale, it is sensible to design for an adsorption column with relatively large aspect ratio (length to diameter ratio). Here, the aspect ratio is kept constant at 3, leading to a column diameter of 0.8 m. The second and third parts concern the costs associated with the purchase, installation, and maintenance of the CoBn/CoEv and the CnEv evacuation pumps, respectively. Two separate pumps for the cocurrent and the counter-current evacuation steps are necessary because the blowdown gas from the CoBn/

product is approximately proportional to amount of gas fed to the column during the HPA step, we may conclude that the amount of N2 (and CH4) fed during the HPA step is practically equal in both cycles at a given recovery. In other words, N2 working capacity remains unchanged in the two cycles. This results from the use of higher adsorption pressure, which tends to increase the working capacity, but neutralizes the effect by feeding more N2 during the RP step. At any point along the Pareto fronts of the 6- and 9-step cycles in Figures 3 and 7, the electricity extractable from combusting the blowdown and evacuation gas stream from the regeneration steps is sufficient to power the PSA unit. To give an illustration, at the recovery-productivity Pareto front of the 9-step cycle for a product containing 95% CH4 and a feed containing 15% N2, the maximum achievable recovery (98.4%) allows production of 4340 J of electricity per mole of CH4 produced, assuming a conservative power generation efficiency of 30%, versus an energy use of 1830 J per mole of CH4 produced. This also illustrates why energy use is not a significant performance measure of the proposed PVSA unit. However, to extract electricity, a separate generator unit will have to be installed. 3.3. Maximum Allowable CH4 Recovery. The performances of the PSA processes in this work are quantified via product recovery and productivity at a given product purity and feed composition. However, at a given product purity and feed composition, there exists a range of product recoveries that are undesirable. The upper limit of this range corresponds to legal requirement for methane emission for discharging the blowdown gas to the atmosphere (500 ppmv proposed by US EPA).25 The lower limit corresponds to the lower flammability limit of the blowdown gas upon blending with air. The upper limit of the recovery range is of the order 99+% (to keep below the legal limit of direct emission of the blowdown gas H


Article


train (Ntrain) needed to effect continuous production, continuous operation of vacuum pumps at a constant volumetric flow rate, and proper alignment of DPE/PPE steps is

CoEv step is recycled as feed for the RP step, and thus cannot be allowed to mix with the blowdown gas from the CnEv step. The fourth part concerns the costs associated with the purchase, installation, and maintenance of the HPA compressor. The HPA compressor is needed to compress the product to pipeline pressure. The fifth and sixth parts relate to the cost of purchasing adsorbents, and the cost of electricity needed to run the process. To improve unit mobility, an electricity generator utilizing the blowdown gas is not installed. Consequently, the cost of electricity incurred is the cost of purchasing electricity from the grid as an industrial consumer. Several general observations should be made regarding the proposed comparison with the two commercial alternatives. One, the reported total capital and operating costs of the Molecular Gate SPEC PLANT and the membrane system exclude the cost of recompression to pipeline pressure.26,27 Two, at the reported total costs, the Molecular Gate SPEC PLANT and the membrane system produce a methaneenriched stream at 15 and 13.8 bar, respectively. Three, at the reported total costs, these two alternatives process a feed containing 10% and 15% N2 respectively. Four, the desired pipeline pressure is not obvious, since it depends on the point of entry into the pipeline. Based on these observations, it is argued that the costing algorithm of the proposed PVSA unit should be dependent on the alternative being compared with. For comparison between the proposed PVSA unit and the Molecular Gate SPEC PLANT, a feed containing 10% N2 is used and the product is compressed to 15 bar. The cost of a HPA compressor is also included for the proposed PVSA unit. For comparison with the membrane system, a feed containing 15% N2 is used and the product is recompressed to 13.8 bar. Once again, the cost of a HPA compressor is included for the proposed PVSA unit. The costing procedure follows the correlations of Turton et al.28 The total capital and operating cost of Molecular Gate SPEC PLANT quoted prior was for 2009, whereas that of the membrane system was for 2010. These are adjusted to 2016 prices using the respective Chemical Engineering Plant Cost Indices (CEPCI). The units are costed assuming a depreciation period of 7.5 years. Operating conditions with methane recovery in excess of the maximum allowable recovery (resulting in an incombustible blowdown gas) is not considered for costing. 3.4.2. Scheduling. The PVSA unit operation is made continuous by grouping several columns to make a train. This process entails selecting sets of operating conditions that lead to a continuous production, as well as operation of CoBn/ CoEv and CnEv pumps. Furthermore, as discussed in Sections 2.2 and 3.2, the DPE step of one column in the train must align with the PPE step of another column in the train. Formally, suppose t HPA a = HPA tcyc bHPA

(14a)

aCoBn/CoEv tCoBn + tCoEv = tcyc bCoBn/CoEv

(14b)

tCnEv a = CnEv tcyc bCnEv

(14c)

Ntrain = LCM{bHPA , bCoBn/CoEv , bCnEv , 2}

(14d)

Here LCM is the least common multiple operator. To illustrate the use of eqs 14a−14d, consider a hypothetical scheduling problem wherein the HPA, CoBn/CoEv and CnEv steps occupy 1/3rd, 1/12th, and 1/6th of the total cycle time, respectively. Then the number of columns per train needed to make the operation continuous is 12. The resultant schedule is shown in Figure 11.

Figure 11. A 12-column schedule corresponding to a cycle wherein HPA, CoBn/CoEv and CnEv steps occupy 1/3rd, 1/6th, and 1/12th of the total cycle time, respectively. At any point in time, 4 columns undergo HPA, 1 column undergoes CoBn/CoEv, and 2 columns undergo CnEv. Each DPE step is coupled to a PPE step.

To prevent the operation of the unit from becoming excessively complex, we constrain the number of columns per train to be less than or equal to 10:

Ntrain ≤ 10

(15)

Equations 14a−14d are applied to all points simulated by NSGA-II. The number of columns calculated from eqs 14a−14d is found to be greater than the limit imposed by eq 15, i.e., there is no simulation point for which Ntrain is less than or equal to 10. This arises due to the fraction of total cycle time occupied by CoBn and CoEv steps. An analysis of the points simulated by NSGA-II indicates that this fraction varies within a narrow range of 0.00542 and 0.0432. Following eq 14b, aCoBn/CoEv < 0.043 × bCoBn/CoEv

(16)

For aCoBn/CoEv to be a positive integer, bCoBn/CoEv must be greater than or equal to 24. The requirement for continuous CoBn/CoEv pump operation thus bounds Ntrain below by 24. Since, from eq 15, the number of columns per train is also bounded above by 10, the feasible set is empty. The requirement for a continuous CoBn/CoEv pump operation thus imposes an undesirable constraint that needs to be removed. Constraint removal is achieved by allowing the evacuation pumps to do work against a valve, as illustrated in Figure 12. valve valve Formally, if tVP,CoBn + CoEv and tVP,CnEv are the durations per cycle for which the CoBn/CoEv and CnEv pumps are doing work against a valve, then

Here tcyc is the total cycle time, and a and b are coprime positive integers. Then the minimum number of columns per I


Article


Figure 13. By allowing work done against valves, the number of columns required per train is reduced to 6. The time during which a valve acts against the CoBn/CoEv vacuum pump corresponds to the area enclosed in the dotted lines.

Figure 12. Illustration of CnEv pump scheduling for one pump serving two columns. Work is done against a valve whenever CnEv step is not active. The CnEv pump operation is now continuous. The same logic may be applied to the CoBn/CoEv pump. valve t VP,CoBn + CoEv

tcyc

valve t VP,CnEv

tcyc

⎛t ⎞ + tCoEv = ceil⎜⎜ CoBn × Ntrain⎟⎟ tcyc ⎝ ⎠ tCoBn + tCoEv − × Ntrain tcyc

⎛t ⎞ t = ceil⎜⎜ CnEv × Ntrain⎟⎟ − CnEv × Ntrain tcyc ⎝ tcyc ⎠

constraints are reflected in eqs 14a−14d. The requirement a continuous CoBn/CoEv pump operation is found to be restrictive, and is removed by allowing the evacuation pumps to do work against a valve (see Figure 12). However, the energy use of the evacuation pumps must then be adjusted to account for the work done against the valve. This idea is encapsulated in eqs 17a and 17b. 3.4.3. Number of Trains. For each properly scheduled point, eq 14d determines the number of columns per train, which in turn determines the amount of feed that can be processed per train:

(17a)

(17b)

Here ceil() is the ceiling operator. Moving back to the sample problem, suppose that there is a cycle wherein HPA, CoBn/ CoEv, and CnEv steps occupy 1/3rd, 1/12th, and 1/6th of the total cycle time, respectively. By allowing CoBn/CoEv and CnEv pumps to do work against a valve, the vacuum pump constraints, i.e., eqs 14b and 14c are automatically satisfied. The minimum number of columns per train may therefore be reduced to 6, as illustrated in Figure 13. The graphical interpretation of eqs 17a and 17b is also shown in the figure. Simulation points fulfilling the remaining constraints (eqs 14a, 14d, and 15) are said to be properly scheduled, and are used for costing. Note that the cost of electricity must be adjusted to account for the work done against the valve; this has been achieved by scaling the energy use during the CoBn/ CoEv and CnEv steps up by some factor f i calculated as follows: fi = 1 +

Ftrain =

(1 − yN ,F )P

(18a)

2

Here V is the volume of the column in m3, yN2,F is the mole fraction of N2 in the feed, and P and Π are the recovery and productivity for the properly scheduled point under consideration. The amount of feed that can be processed in a mobile PVSA unit is given by F = Ftrainntrain

(18b)

Here ntrain is the number of trains in the PVSA unit. This introduces an unusual discontinuity in the costing problem because not all feed throughputs are achievable. To illustrate the point, supposing that Ftrain is equal to 1.2 MMscfd, the possible feed throughputs are 1.2 MMscfd (corresponding to 1 train), 2.4 MMscfd (corresponding to 2 trains), 3.6 MMscfd (corresponding to 3 trains), and etc. In this work, we have opted to perform costing for each possible feed throughput corresponding to each properly scheduled point. More details on optimization of the discontinuous costing problem are discussed in the subsequent sections. 3.4.4. Adsorption Column. The base and the bare module costs (Cc and CB,c respectively) are calculated from the correlations:

valve t VP, i

ti × Ntrain

(1 − εb)Vρs × Ntrain × Π

(17c)

Here i runs across CoBn/CoEv and CnEv steps. In reference to Figure 13, the energy use of the CoBn/CoEv pump is thus increased by a factor equal to the ratio of the area labeled tvalve VP,CoBn to the area corresponding to the CoBn steps. To summarize, four constraints on the scheduling of the proposed PVSA unit was originally considered (continuous production, continuous operation of CoBn/CoEv and CnEv pumps, and proper coupling of DPE and PPE steps). These

log Cc = 3.4974 + 0.4485log V + 0.1074(logV )2 J

(19a)


Article

Industrial & Engineering Chemistry Research C B,c =

ICEP,2016 ICEP,2001

Cc × (2.52 + 1.82FMFP)

enriched stream. The bare module costs corresponding to eq 21a and eq 21b are

(19b)

C B,HPA =

FM and FP are factors taking into account the material and the operating pressure of the column respectively, ICEP,2016 is the CEPCI for 2016, and ICEP,2001 is the CEPCI for 2001. For an adsorption column, the factor FM takes on the value 1, whereas FP is calculated from (P + 1)D ⎧ ⎫ ⎪ ⎪ 2 × [850 − 1.5(P + 1)] + 0.00315 , 1.25⎬ FP = max⎨ 0.0063 ⎪ ⎪ ⎭ ⎩

C B,d =

(19c)

(20a)

Here Ql is the pumping capacity of the evacuation pump corresponding to step l in m3/h, which may be either CoBn/ CoEv or CnEv. Note that significant cost reduction may be achieved by sharing the evacuation pump across the trains. Also, depending on the selected schedule, the number of columns (or valves) subjected to the CoBn/CoEv or the CnEv steps may vary. Consequently, the pumping capacity is calculated as

Cl × (1.89 + 1.35FMFP)

(21d)

(22a)

CGR = C T,B + 0.35(ntrainNtrain × C B,c + C B,CoBn/CoEv + C B,CnEv + C B,HPA + C B,d)

(22b)

The grassroots fee accounts for auxiliary facilities, which in this case is the vehicle on which the unit is mounted. The total grassroots cost is then amortized over a period of 7.5 years with a 10% interest rate (taken from the Molecular Gate SPEC PLANT cost report23) to give f, and combined with the adsorbent and electricity costs to give the total treatment cost (CT):

(20c)

Here FM and FP are equal to unity. Note that 2014 is the year in which eq 20a was developed. 3.4.6. HPA Compressor. The base cost of the HPA compressor without a drive (CHPA) and the base cost of the drive (Cd) are calculated from the correlations: log C HPA = 2.9945 + 0.9542 log(ηẆHPA)

C HPA × 1.5

The total grassroots cost (CGR) is obtained by adding similar grassroots (35%) fee:

Here vF,l is the evacuation velocity. Equation 20a is only valid for pumping capacities below 26 000 m3/h; however, this limit is not reached for the maximum feed throughput investigated in this study, i.e., 10 MMscfd. Note that the evacuation velocity is one of the decision variables of the optimization problem. The corresponding bare module costs of the evacuation pumps are ICEP,2014

ICEP,2001

(21c)

+ C B,HPA + C B,d)

(20b)

ICEP,2016

ICEP,2016

C HPA × 2.5

C T,B = 1.18(ntrainNtrain × C B,c + C B,CoBn/CoEv + C B,CnEv

⎛t ⎞ + tCoEv πD 2 vF, lεbntrain × ceil⎜⎜ CoBn × Ntrain⎟⎟ Ql = 4 tcyc ⎝ ⎠

C B, l =

ICEP,2001

3.4.7. Adsorbent and Electricity. The price of Takeda CMS is taken to be $2/kg-ads30 and it is assumed that the adsorbent will be replaced every year. The electricity required to run the process is taken to be the sum of that required to run the CoBn/CoEv and CnEv pumps, and to recompress the product stream to 13.8 bar. These are calculated assuming isentropic processes with constant efficiencies (0.72 for compressors, 0.6 for evacuation pumps) and a fractional engineering downtime of 0.01, whereas the cost of electricity is assumed to be 0.08$/ (kW·h). Although the isentropic efficiency of an evacuation pump is known to be a function of evacuation pressure, the small relative contribution of the electricity used by the vacuum pumps to the total treatment cost (between 1.9 and 5.1%) justifies the use of a constant efficiency. More details on the calculation of electricity use from isentropic compression equation may be found in.18 3.4.8. Total Cost. The total bare module cost (CT,B) is obtained by taking a sum of the bare module costs of the adsorption column (multiplied with the number of columns), the evacuation pumps and the HPA compressor, and adding to that sum contingency (15%) and engineering (3%) fees:

Here P is the maximum pressure in the column expressed in barg, and D is the diameter of the column in m. The default minimum value of 1.25 corresponds to a low pressure operation of

Optimization of a Pressure Swing Adsorption Process for Nitrogen

Recommend Documents