Optimal Design of a Two-Stage Membrane System for Hydrogen Separation in Refining Processes

This paper fits into the process system engineering field by addressing the optimization of a two-stage membrane system for H2 separation in refinery processes. To this end, a nonlinear mathematical programming (NLP) model is developed to simultaneously optimize the size of each membrane stage (membrane area, heat transfer area, and installed power for compressors and vacuum pumps) and operating conditions (flow rates, pressures, temperatures, and compositions) to achieve desired target levels of H2 product purity and H2 recovery at a minimum total annual cost. Optimal configuration and process design are obtained from a model which embeds different operating modes and process configurations. For instance, the following candidate ways to create the driving force across the membrane are embedded: (a) compression of both feed and/or permeate streams, or (b) vacuum application in permeate streams, or (c) a combination of (a) and (b). In addition, the potential selection of an expansion turbine to recover energy from the retentate stream (energy recovery system) is also embedded. For a H2 product purity of 0.90 and H2 recovery of 90%, a minimum total annual cost of 1.764 M$·year−1 was obtained for treating 100 kmol·h−1 with 0.18, 0.16, 0.62, and 0.04 mole fraction of H2, CO, N2, CO2, respectively. The optimal solution selected a combination of compression and vacuum to create the driving force and removed the expansion turbine. Afterwards, this optimal solution was compared in terms of costs, process-unit sizes, and operating conditions to the following two sub-optimal solutions: (i) no vacuum in permeate stream is applied, and (ii) the expansion turbine is included into the process. The comparison showed that the latter (ii) has the highest total annual cost (TAC) value, which is around 7% higher than the former (i) and 24% higher than the found optimal solution. Finally, a sensitivity analysis to investigate the influence of the desired H2 product purity and H2 recovery is presented. Opposite cost-based trade-offs between total membrane area and total electric power were observed with the variations of these two model parameters. This paper contributes a valuable decision-support tool in the process system engineering field for designing, simulating, and optimizing membrane-based systems for H2 separation in a particular industrial case; and the presented optimization results provide useful guidelines to assist in selecting the optimal configuration and operating mode.


Introduction
Currently, membranes are playing an important role in a wide range of industrial applications [1].If compared with other separation processes (absorption, adsorption, and cryogenic processes), membrane-based separation processes require simple pieces of equipment which are easy to operate (compressors and/or vacuum pumps, membrane modules) and imply low investment and operating costs [2].However, despite these advantages, research efforts are still needed to enhance separation efficiency of membrane-based processes [3,4].
This paper fits into the process system engineering (PSE) field, which is, together with material science, an important research area towards the enhancement of membrane-based separation processes.In PSE, mathematical modeling and simulation constitute a valuable tool, not only for understanding the systems' behavior, but also for both finding or developing novel processes and optimizing them by considering different criteria (efficiency, cost, cost-effectiveness, etc.).Certainly, modeling and simulation of membrane modules and processes are of great importance to ensure a successful transfer into industrial applications.
In order to show the vast variety of mathematical model types, solution strategies, and computational tools, we here include a brief overview of some articles that provided contributions to membrane-based separation in the process system engineering field, regardless of any particular membrane material, geometry and flow pattern, and type of the treated gas mixture.
As regards commercial and in-house process simulators, Ahmad et al. [5] developed a Visual Basic sub-routine of a single-stage membrane model which was included into process simulator Aspen HYSYS with the aim of simulating different multi-stage membrane configurations.In a similar way, He and Hägg [6] combined Aspen HYSYS with an in-house membrane program (ChemBrane) to evaluate the techno-economic feasibility of different membrane systems consisting of hollow fiber carbon membranes and fixed site carrier membranes for CO 2 removal from natural gas.Low et al. [7] presented a parametric simulation study of a single-stage membrane process for capturing CO 2 from a flue gas, considering real gas behavior and using process simulator PRO/II.They simulated a cross-flow configuration and used Peng-Robinson equation of state as thermodynamic model.The authors investigated the influence of feed compression and/or permeate vacuum, moisture level and temperature on the CO 2 separation degree, membrane area requirement, and energy consumption.A result showed that the higher the membrane CO 2 permeance, the lower the membrane area requirements; while the higher CO 2 /N 2 selectivity, the higher CO 2 purity and the lower energy requirement.Process simulator ChemCAD has been applied to study membrane-based processes in several applications, such as production of oxygen-enriched air [8], CO 2 capture in power plants [9,10], syngas purification [11], and methane enrichment process [12].For instance, Merkel et al. [9] performed a techno-economic comparison of two membrane process designs for 90% CO 2 capture from a 600 MW coal-fired power plant.The authors concluded that a two-step counter-current sweep design requires less power and membrane area than a two-step/two-stage design.Bounaceur et al. [13] employed a CAPE-OPEN compliant simulation software MEMSIC to simulate two case studies: natural gas treatment (CO 2 /CH 4 separation through a glassy polymer) and recovery of volatile organic compounds (acetone/N 2 separation through a rubbery polymer).MEMSIC allows simulating membrane modules, considering the possibility of selecting different flow configurations and conditions (ideal mixing, co-current and counter-current plug flow patterns) and different transport mechanisms of molecular species through a membrane (constant permeability, Henry law, etc.).
With respect to numerical methods, Makaruk and Harasek [14] presented an iterative method based on finite difference Gauß-Seidel method to calculate the permeation of a gas mixture consisting of CH 4 , CO 2 , and O 2 in a study of membrane systems involving hollow fiber modules.Co-current, counter-current, and cross-flow configurations were simulated.The authors considered a relaxation factor to avoid convergence problems.Kundu et al. [15] developed a solution technique based on Gear's backward differentiation formulae (BDF) method to solve the ordinary differential equations as an initial value problem in two successive steps to simulate asymmetric hollow fiber membranes for air separation.They concluded that the procedure leads to minimum computational time with improved solution stability and that computational complexity does not depend on the number of components.
When designing cost-efficient multi-stage membrane systems, which are required to simultaneously achieve high product purity and high recovery levels [16][17][18], all trade-offs among the model variables must be investigated.The number of trade-offs rises along with the number of stages.Therefore, a process-level analysis is highly suggested by several authors [19][20][21].Thus, the application of the aforementioned process simulators and sequential methods to optimize multi-stage processes may be computationally expensive (time-consuming and large number of iterations) due to the presence of recycle streams [14].Moreover, these methods may also present convergence problems for optimization purposes if no good initial values are guessed in the iterative solution procedures.
In relation to optimization methods to deal with these trade-offs, genetic algorithms (GA) have been recently applied to optimize membrane-based separation processes [16,22,23].Yuan et al. [22] addressed the study of multi-stage membrane processes employing GA and using the function 'gamultiobj' supported in MATLAB to solve a two-criteria optimization problem (minimization of energy and total membrane area).In line with [22], Gabrielli et al. [16] addressed the optimal design of membrane-based gas separation processes to treat a CO 2 /N 2 binary mixture typical of the flue gas of a coal-fired power plant after drying by proposing a multi-objective optimization problem implemented in MATLAB and solved with GA.As a result, the layout of the process, operating variables, and membrane material were obtained.Attainable regions for purity-recovery specifications were presented through Pareto sets in term of energy and membrane area.In general, no specific knowledge about the problem is required by GAs, which is one of the main advantages of this approach.Also, GAs can be easily parallelized in computer clusters.They are derivate-free and well suited for high complexity problems with discontinuous models and for parameter estimation problems.The solutions obtained by GAs are strongly dependent on the required parameters; for instance, the number of generations, population, crossover rate, and mutation rate, among others.
There are software tools for designing experiments (DoE), such as Modde ® , that are used by researchers for membrane process optimization.Since they are not formally based on first principles and no extensive knowledge of the field is required, they can be easily adopted.Some application examples dealing with membrane systems can be found in [24,25].
Mathematical programming tools and robust optimization techniques can be used to tackle combinatorial, large, and highly nonlinear optimization problems and to deal with the existing trade-offs in any process and system.Few articles addressing methods for the simultaneous optimization of multi-stage membrane systems can be found in the literature [16,23,[26][27][28][29][30].Qi and Henson [26] and Scholz et al. [27] developed mixed-integer nonlinear programming (MINLP) models for the optimization of multi-stage membrane systems using GAMS (General Algebraic Modeling System) as implementation tool.Discrete decisions are related to the selection of stages and are modeled by using integer variables.Qi and Henson [26] employed their model to study several case studies such as natural gas sweetening and enhanced oil recovery, considering multi-component separations and using DICOPT (DIscrete and Continuous OPTimizer) as optimization solver.They highlighted the robustness of the proposed MINLP model.Scholz et al. [27] used a branch-and-bound type global optimization solver (Branch-And-Reduce Optimization Navigator BARON) to optimize a biogas upgrading process.The used solver guarantees finding the global optimum solution, which is an important aspect of PSE.It should be noted that local search optimization algorithms (CONOPT, MINOS) cannot guarantee the attained solution to be a global one.Recently, Zarca et al. [30] concluded that a two-stage membrane process can be successfully applied for H 2 recovery by using a polymeric membrane, and for syngas recovery by using a polymer-ionic liquid composite membrane.Their conclusion is supported by a techno-economic analysis based on a nonlinear mathematical programming (NLP) model developed for this purpose, which was implemented and solved in GAMS/CONOPT environment.The net present value was proposed as the objective function to be minimized.They assumed that retentate streams are at atmospheric pressure and they did not consider the presence of an expander for power recovery.The authors highlighted the usefulness of the mathematical modeling and optimization tools in the PSE field in general and their use for improving membrane system performance in particular.
In this research line, the motivation of this paper is to extend the optimization problem proposed in Zarca et al. [30] by additionally considering the possibility of both including an expansion turbine to recover power and using vacuum to create the driving force for permeation, aiming at minimizing the total annual cost of the process while meeting desired H 2 product purity and H 2 recovery target levels.The inclusion of an expansion turbine and vacuum raises the degrees of freedom of the resulting optimization problem and, consequently, the trade-offs exiting between the model variables, as it is clearly explained in the next section devoted to the process description.

Process Description
A schematic of the studied two-stage membrane system for hydrogen separation is illustrated in Figure 1.As shown, the feed stream F 0 passes through the compressor C1 and the heat exchanger HEX1 to increase the pressure and to reach the operating temperature, respectively.Afterward, F 0 can be optionally mixed with a fraction of the retentate stream obtained in the first membrane stage MS1 (RR MS1 ) and/or the obtained in the second membrane stage MS2 (RR MS2_MS1 ).Then, the resulting stream F MS1 is fed to the membrane stage MS1, where it is separated into two streams: a stream that permeates through the membrane (permeate P MS1 ) and a stream that is retained by the membrane (retentate R MS1 ).
In this study, the driving force for selective permeation can optionally be created by different ways: (a) compressors C1 and/or C2, (b) vacuum pumps VP1 and/or VP2, and (c) a combination of compressors and vacuum pumps.In other words, the ways to create the driving force in each membrane stage will be a result of the optimization problem, specifically through p H (retentate pressure) and p L (permeate pressure) which are optimization variables of the process model.(Of course, the driving force also depends on the mole fractions of the components at both membrane sides).For instance, if no vacuum pumps are selected, p L will be the atmospheric pressure.On the other hand, if no compressors are selected, the driving force will be created by vacuum pumps and p H will be the atmospheric pressure.Thus, the formulated mathematical model includes the possibility to select (feed and permeate) compression and/or application of vacuum in the permeate streams.
Processes 2018, 6, x FOR PEER REVIEW 4 of 23 atmospheric pressure and they did not consider the presence of an expander for power recovery.The authors highlighted the usefulness of the mathematical modeling and optimization tools in the PSE field in general and their use for improving membrane system performance in particular.In this research line, the motivation of this paper is to extend the optimization problem proposed in Zarca et al. [30] by additionally considering the possibility of both including an expansion turbine to recover power and using vacuum to create the driving force for permeation, aiming at minimizing the total annual cost of the process while meeting desired H2 product purity and H2 recovery target levels.The inclusion of an expansion turbine and vacuum raises the degrees of freedom of the resulting optimization problem and, consequently, the trade-offs exiting between the model variables, as it is clearly explained in the next section devoted to the process description.

Process Description
A schematic of the studied two-stage membrane system for hydrogen separation is illustrated in Figure 1.As shown, the feed stream F0 passes through the compressor C1 and the heat exchanger HEX1 to increase the pressure and to reach the operating temperature, respectively.Afterward, F0 can be optionally mixed with a fraction of the retentate stream obtained in the first membrane stage MS1 (RRMS1) and/or the obtained in the second membrane stage MS2 (RRMS2_MS1).Then, the resulting stream FMS1 is fed to the membrane stage MS1, where it is separated into two streams: a stream that permeates through the membrane (permeate PMS1) and a stream that is retained by the membrane (retentate RMS1).
In this study, the driving force for selective permeation can optionally be created by different ways: (a) compressors C1 and/or C2, (b) vacuum pumps VP1 and/or VP2, and (c) a combination of compressors and vacuum pumps.In other words, the ways to create the driving force in each membrane stage will be a result of the optimization problem, specifically through p H (retentate pressure) and p L (permeate pressure) which are optimization variables of the process model.(Of course, the driving force also depends on the mole fractions of the components at both membrane sides).For instance, if no vacuum pumps are selected, p L will be the atmospheric pressure.On the other hand, if no compressors are selected, the driving force will be created by vacuum pumps and p H will be the atmospheric pressure.Thus, the formulated mathematical model includes the possibility to select (feed and permeate) compression and/or application of vacuum in the permeate streams.In addition, if the optimal solution selects the inclusion of the compressor C1, the mathematical model has the option to include the expander EXP to recover mechanical power from the retentate stream R MS1 .The model also includes the possibility of recycling part of R MS2 back to MS2 (RR MS2 ) and/or to MS1 (RR MS2_MS1 ).
As it can be noticed, there are different trade-offs between the pressure ratio (p H /p L ), the electric power requirement to run compressors and vacuum pumps, and the membrane area.For instance, the higher the pressure ratio, the higher the H 2 purity in the permeate stream and the lower the membrane area but the higher the electric power requirement to run the compressor.The optimal pressure ratio values depend on the relationships between investment and operating costs.However, these trade-offs can show opposite trends if the expander EXP is included into the analysis since the higher the retentate pressure, the higher the amount of power recovered by the expander and, therefore, the lower the total net amount of electric power required to operate the process.
Based on this, it is clear that the finding of optimal solutions could be a difficult and time-consuming task if parametric optimization is employed instead of simultaneous optimization as is proposed in this work.Thus, the relevance of this paper is to simultaneously optimize the trade-offs that are difficult to distinguish at first glance.As a result of the model, the following process design decisions that minimize the total annual cost (capital and operating expenditures) are determined by:

•
The way in which the driving force for permeation is created in each membrane stage, i.e., the solution will indicate if the driving force is created by (a) employing compression without applying vacuum, or (b) applying vacuum in the permeate streams without employing compression, or (c) combining both compression and vacuum.

•
The inclusion or not of the expansion turbine for power recovery.

•
The inclusion or not of (one or more) recycle streams.

•
The optimal values of operation pressure, temperature, composition, and flow rate of each process stream.

•
The optimal sizes of the process units (membrane areas, heat transfer areas, power capacity of compressors and vacuum pumps).
The specifications of the gas mixture feed and the desired H 2 product purity and H 2 recovery are assumed as model parameters.In this study, H 2 product purity and H 2 recovery target levels of 0.9 and 90%, respectively, are specified, which are typical values for fluid catalytic cracking (FCC) off-gases and other refinery purge streams [31].As a first approximation, the same polymeric material is considered for both membrane modules and the same permeance values are used in both cases.

Main Model Assumptions
The following assumptions are considered for modeling the membrane units [20]:

•
All the mixture components can permeate through the membrane.

•
The operating pressure does not affect the component permeability.

•
No pressure drop is considered in the retentate and permeate sides.

•
The feed and retentate streams are at the same pressure.

•
Plug flow pattern is considered at both sides of the membrane unit.

•
Isothermal condition is assumed within the membrane module.

•
Fick's first law is used.
Next, the mathematical model that describes the complete process, including the cost model to calculate the total annual cost in terms of the capital and operating expenditures, is presented.

Mathematical Model
Figure 2 shows a schematic of the whole process and the membrane module, and the nomenclature used to introduce the mathematical model.

Mathematical Model
Figure 2 shows a schematic of the whole process and the membrane module, and the nomenclature used to introduce the mathematical model.

Mass Balances
By applying the backward finite difference method (BFDM) to the membrane module MS1 illustrated in Figure 2b, the following set of algebraic equations is obtained: (F indicates flow rate, Q heat load, T temperature, p pressure, x composition in feed/retentate streams, y composition in permeate streams, W power, A area of membrane and area of heat exchanger.The subscript i represents a mixture component, subscript j a discretization point along the membrane length, and subscript 0 the fresh feed).

Mass Balances
By applying the backward finite difference method (BFDM) to the membrane module MS1 illustrated in Figure 2b, the following set of algebraic equations is obtained: (J−1) where ξ i and A MS1 refer to the permeance of component i and membrane surface area, respectively.p H and p L MS1 are the operating pressures in the retentate and permeate sides, respectively.The index j refers to a discretization point which varies from 0 to 19 (J = 19, i.e., 20 discretization points is considered).
A similar set of constraints is formulated for the second membrane stage MS2.The mass balances in splitters SP1 and SP2 are expressed as follows:
Finally, the power W EXP recovered in the expander (isothermal expansion) is given by Equation ( 17): As mentioned earlier, the model was developed in such a way that the expander can be selected or removed from the optimal solution, depending on the optimal flow rate values in the splitter SP1 (Equation ( 7)).For instance, if the optimal value of R MS1 ## in Equation ( 7) is equal to 0, then the expander is removed and no electric power is recovered.

Energy Balances and Transfer Areas of Heat Exchangers
The energy balances in the heat exchangers and the corresponding heat transfer areas are calculated by Equations ( 18)-( 26): LMTD HEX1 = LMTD HEX2 = LMTD HEX3 = The parameter U refers to the overall heat transfer coefficient, which is assumed to be 277.7 W•m −2 •K −1 for all heat exchangers.

Connecting Constraints
The following constraints are used to relate model variables defined inside and outside of the membrane module MS1: x MS1,i = x MS1,i,j = 0 (28) x MS1,i,j Similar constraints are considered for the second membrane stage MS2.

Performance Variables
The total membrane area (TMA) and the total heat transfer area (THTA) are calculated by Equations ( 33) and (34), respectively: The total and net power requirements (TW and TNW, respectively) are calculated as follows: TNW = TW − W EXP (36)
The numerical values of the model parameters used in this optimization study are taken from Zarca et al. [30], which are listed in Table 1.
The resulting NLP optimization model was implemented in GAMS and solved with CONOPT.

Feed specification
Flow rate, kmol

Results and Discussion
Figure 3 illustrates the optimal process configuration obtained by solving the stated optimization problem, which is hereafter named solution 'OPT'.Figure 4 shows the percentage contribution of the annualized CAPEX (annCAPEX) and OPEX to the TAC value (Figure 4a), the percentage contribution of each process unit to the total equipment acquisition investment C INV (Figure 4b), and the percentage contribution of the electric cost, cooling water cost, and membrane replacement cost to the total utility cost (Figure 4c).As shown in Figure 3, the optimal configuration involves two compressors (C1 and C2) and one vacuum pump (VP1).As the optimal value of R EMS1 is equal to 0, the expander was excluded from the optimal configuration.Figure 4a shows that the minimum TAC value is 1.764 M$•year −1 , where the annCAPEX contribution is more significant than the OPEX contribution (62% vs. 38%).
Figure 3 2).In the first membrane stage, the feed compressor C1 requires 0.197 MW and the vacuum pump VP1 requires 0.048 MW, leading to an operating pressure ratio of 29.9 (0.598/0.020).In the second membrane stage, the electric power required in the compressor C2 for compressing the permeate obtained in the first stage is 0.053 MW.Thus, the required total electric power capacity to run the process is 0.298 MW, which represents 75.9% of the C INV (1.087 of 1.431 M$, see Table 2).
Regarding heat transfer requirements, the heat exchangers require in total 13.1 m 2 of heat transfer area, which represent only 2.9% of the C INV (0.041 of 1.430 M$, see Table 2).The results presented in Figure 4b show that the compressor C1 is the largest contributor to the C INV , followed by the compressor C2 and the membrane area A MS1 of the first stage, contributing with 48.5%, 22.1%, and 18.8%, respectively.Clearly, in each membrane stage, the compressor is significantly more expensive than the membrane itself (0.694 vs. 0.269 M$ in the first stage, and 0.316 vs. 0.034 M$ in the second one, see Table 2).Regarding the optimal distribution of the total utility cost CRM, it can be seen in Figure 4c that the cost for electric power demand CEP is the largest contributor to CRM with 90.8%, followed by the cost for membrane replacement CMR with 7.3%.The contribution of the cost associated to the cooling water CCW is insignificant as it represents only 1.7% of the CRM value.
(a) (b) Regarding the optimal distribution of the total utility cost C RM , it can be seen in Figure 4c that the cost for electric power demand C EP is the largest contributor to C RM with 90.8%, followed by the cost for membrane replacement C MR with 7.3%.The contribution of the cost associated to the cooling water C CW is insignificant as it represents only 1.7% of the C RM value.created by combining compression and vacuum and no expander is included.
Regarding the optimal distribution of the total utility cost CRM, it can be seen in Figure 4c that the cost for electric power demand CEP is the largest contributor to CRM with 90.8%, followed by the cost for membrane replacement CMR with 7.3%.The contribution of the cost associated to the cooling water CCW is insignificant as it represents only 1.7% of the CRM value.Figure 5 shows the optimal composition profiles of all the components in the retentate and permeate streams along the membrane length (which is represented by the discretization points) of the permeation stages MS1 and MS2.In this separation, the permeate profile is the most interesting profile to analyze because it concentrates H2, which is the desired component in the product stream.Figure 5 shows the optimal composition profiles of all the components in the retentate and permeate streams along the membrane length (which is represented by the discretization points) of the permeation stages MS1 and MS2.In this separation, the permeate profile is the most interesting profile to analyze because it concentrates H 2 , which is the desired component in the product stream.
It can be observed in Figure 5a that the H 2 mole fraction in the permeate stream in the first membrane increases from 0.370 at the discretization point j = 19 (at the feed exit) to 0.710 at the discretization point j = 0 (where the permeate stream leaves the membrane module MS1).Similarly, Figure 5b shows how the H 2 mole fraction in the permeate stream in the second membrane increases from 0.813 at the discretization point j = 19 to 0.90 at j = 0 (where the desired design specifications of 0.90 H 2 purity and 90% H 2 recovery are achieved).Figure 5c,d show the mole fraction profiles of the components in the retentate streams in MS1 and MS2, respectively, which flow in opposite direction to the permeate streams.The H 2 mole fraction in the retentate stream in MS2 (Figure 5d) reaches the lowest value of 0.363 at the discretization point j = 19 i.e., where the H 2 mole fraction in the permeate stream also reaches its lowest value (0.813).Thus, as the retentate flows inside the membrane module, its H 2 mole fraction decreases from the highest value of 0.71 to the lowest value of 0.363.Finally, it can be clearly observed in Figure 5c,d that the CO 2 , CO, and N 2 mole fractions increase as the H 2 mole fraction in the retentate stream decreases.
to the permeate streams.The H2 mole fraction in the retentate stream in MS2 (Figure 5d) reaches the lowest value of 0.363 at the discretization point j = 19 i.e., where the H2 mole fraction in the permeate stream also reaches its lowest value (0.813).Thus, as the retentate flows inside the membrane module, its H2 mole fraction decreases from the highest value of 0.71 to the lowest value of 0.363.Finally, it can be clearly observed in Figure 5c,d that the CO2, CO, and N2 mole fractions increase as the H2 mole fraction in the retentate stream decreases.

Comparison of Optimal and Sub-Optimal Solutions
In order to investigate how the optimal solution changes if other ways to create the driving force for permeation are applied, and if the expander is forced to be part of the process configuration, the same process optimization model is solved for the following cases:

Comparison of Optimal and Sub-Optimal Solutions
In order to investigate how the optimal solution changes if other ways to create the driving force for permeation are applied, and if the expander is forced to be part of the process configuration, the same process optimization model is solved for the following cases: 1.
The driving forces are created by compression of the feed and permeate streams only.To do this, the (low) pressures in the permeate streams were set to the atmospheric pressure.The optimal solution obtained for this case is shown in Figures 6 and 7, and is hereafter named 'SUBOPT1'.

2.
The expander is included in the process configuration.To do this, the flow rate R MS1 ## in Equation ( 7) (mass balance in the splitter SP1) is set to zero.The optimal solution obtained for this case is illustrated in Figures 8 and 9, and is hereafter named 'SUBOPT2'.
Both solutions SUBOPT1 and SUBOPT2 are sub-optimal solutions compared to the optimal solution OPT.Table 2 compares the costs obtained for these three solutions.The comparison shows that the TAC value obtained for SUBOP1 is 15.5% higher than that obtained for OPT (2.038 vs. 1.764M$•year −1 ) as consequence of the increase of both OPEX by 15.2% (from 1.095 to 1.262 M$•year −1 ) and the annCAPEX by 16.1% (from 0.669 to 0.776 M$•year −1 ).Also, it can be seen in Figure 6 that the total membrane area required in SUBOP1 is 24.6% lower than that required in OPT (4298.07 vs. 5701.67m 2 ) but contrarily the electric power required in the former is 41.0% higher than that required in the latter (0.420 vs. 0.298 MW).Thus, the increase of the annCAPEX value in SUBOPT1 is because of the fact that the increase of the investment associated to the compressors is more important than the decrease of the investment associated to the membrane areas.Regarding the OPEX, a 0.122 MW increase of the total electric power in SUBOPT1 with respect to OPT leads to a 5.77 × 10 −2 M$•year −1 increase of the cost for electricity required to run the compressors (from 0.141 to 0.198 M$•year −1 ).In addition, as in SUBOPT1 the driving force required in the membrane stages is created only by compressors, the operation pressure and temperature of the compressed stream in SUBOPT1 are higher than those obtained in OPT (1.252 vs. 0.598 MPa and 642.3 vs. 520.1 K, respectively).Then, in order to reach the operating temperature in both membrane stages (313.15K), the cost associated with the cooling water CW in SUBOPT1 increased by 1.11 × 10 −3 M$•year −1 with respect to that in OPT (from 2.79 × 10 −3 to 3.90 × 10 −3 M$•year −1 ).Finally, the cost associated to the membrane replacement in SUBOPT1 is 24.6% lower than that in OPT due to the required total membrane area decreases by around 1403.60 m 2 .Figure 7 shows the component composition profiles in the retentate and permeate streams in the membrane modules MS1 and MS2 corresponding to SUBOPT1.The comparison of results presented in Table 2 shows that the sub-optimal solution SUBOPT2 has the highest TAC value, which is 23.7% and 7.1% higher than those obtained in OPT and SUBOPT1, respectively.This is a confirmation that the expander was correctly removed (in terms of costs) from the optimal OPT and the sub-optimal SUBOPT1 solutions.In SUBOPT2, the percentage contributions of OPEX and annCAPEX to the TAC are 58.3% and 41.7%, respectively.The net electric power demand in SUBOPT2 is 45.2% and 61.2% lower than those in OPT and SUBOPT1, respectively, but the total membrane area in SUBOPT2 is 26.6% and 67.9% higher than those obtained in OPT and SUBOPT1, respectively.Then, the investment required by the expander for power recovery in SUBOPT2 (which represents 25% of the CINV) is more important than the decrease of the cost for electric power i.e., in OPEX savings.
Also, it can be seen in Table 2 that the inclusion of the expander modifies the relevance order of the process-unit contributions to the total equipment acquisition investment CINV, with the compressor C1 still being the largest contributor but now followed by the membrane area of the first stage, instead of the compressor C2 like in OPT and SUBOPT1.Figure 9 shows the component composition profiles in the retentate and permeate streams in the membrane modules MS1 and MS2 corresponding to SUBOPT2.The comparison of results presented in Table 2 shows that the sub-optimal solution SUBOPT2 has the highest TAC value, which is 23.7% and 7.1% higher than those obtained in OPT and SUBOPT1, respectively.This is a confirmation that the expander was correctly removed (in terms of costs) from the optimal OPT and the sub-optimal SUBOPT1 solutions.In SUBOPT2, the percentage contributions of OPEX and annCAPEX to the TAC are 58.3% and 41.7%, respectively.The net electric power demand in SUBOPT2 is 45.2% and 61.2% lower than those in OPT and SUBOPT1, respectively, but the total membrane area in SUBOPT2 is 26.6% and 67.9% higher than those obtained in OPT and SUBOPT1, respectively.Then, the investment required by the expander for power recovery in SUBOPT2 (which represents 25% of the C INV ) is more important than the decrease of the cost for electric power i.e., in OPEX savings.
Also, it can be seen in Table 2 that the inclusion of the expander modifies the relevance order of the process-unit contributions to the total equipment acquisition investment C INV , with the compressor C1 still being the largest contributor but now followed by the membrane area of the first stage, instead of the compressor C2 like in OPT and SUBOPT1.Figure 9 shows the component composition profiles in the retentate and permeate streams in the membrane modules MS1 and MS2 corresponding to SUBOPT2.The comparison of results presented in Table 2 shows that the sub-optimal solution SUBOPT2 has the highest TAC value, which is 23.7% and 7.1% higher than those obtained in OPT and SUBOPT1, respectively.This is a confirmation that the expander was correctly removed (in terms of costs) from the optimal OPT and the sub-optimal SUBOPT1 solutions.In SUBOPT2, the percentage contributions of OPEX and annCAPEX to the TAC are 58.3% and 41.7%, respectively.The net electric power demand in SUBOPT2 is 45.2% and 61.2% lower than those in OPT and SUBOPT1, respectively, but the total membrane area in SUBOPT2 is 26.6% and 67.9% higher than those obtained in OPT and SUBOPT1, respectively.Then, the investment required by the expander for power recovery in SUBOPT2 (which represents 25% of the CINV) is more important than the decrease of the cost for electric power i.e., in OPEX savings.
Also, it can be seen in Table 2 that the inclusion of the expander modifies the relevance order of the process-unit contributions to the total equipment acquisition investment CINV, with the compressor C1 still being the largest contributor but now followed by the membrane area of the first stage, instead of the compressor C2 like in OPT and SUBOPT1.Figure 9 shows the component composition profiles in the retentate and permeate streams in the membrane modules MS1 and MS2 corresponding to SUBOPT2.

Sensitivity Analysis
This section presents a sensitivity analysis to see how the optimal solution discussed in the previous section is affected when the model parameters related to the design specifications i.e., H2 product purity and H2 recovery levels are varied one at the time by keeping constant the values of the remaining parameters.

Sensitivity of the Optimal Solution to the H2 Product Purity Level
The optimization results obtained by varying the H2 product purity in −0.01 (−1.11%) and +0.01 (+1.11%) with respect to the value used in the previous section (0.90) are compared in Table 3.

Sensitivity Analysis
This section presents a sensitivity analysis to see how the optimal solution discussed in the previous section is affected when the model parameters related to the design specifications i.e., H 2 product purity and H 2 recovery levels are varied one at the time by keeping constant the values of the remaining parameters.

Sensitivity of the Optimal Solution to the H 2 Product Purity Level
The optimization results obtained by varying the H 2 product purity in −0.01 (−1.11%) and +0.01 (+1.11%) with respect to the value used in the previous section (0.90) are compared in Table 3.
According to Table 3, the TAC value increases by 2.11% (from 1.764 M$•year −1 to 1.802 M$•year −1 ) when the H 2 product purity increases in 0.01 with respect to the reference case (0.90 H 2 purity).This increase is due to the increase of both the annCAPEX and OPEX values by 2.24% and 2.04%, respectively.It is interesting to note that, compared to the reference case, the electric power required by compressors and vacuum pump increases in total 0.015 MW while the optimal total membrane area decreases 134.4 m 2 .The increased H 2 product purity level (0.91) implies to increase the H 2 concentration in the permeate stream of the first stage to reach the H 2 purity in the permeate stream of the second stage (the target level of 0.91).Then, the optimal cost-based trade-offs indicate that it is more beneficial to increase the total electrical power (from 0.598 MPa to 0.621 MPa) to operate the process at a higher operating pressure value p H rather than to increase the total membrane area.
On the other hand, the TAC value decreases by 1.32% (from 1.764 M$•year −1 to 1.741 M$•year −1 ) when the H 2 product purity decreases in 0.01 with respect to the reference level.In this case, the annCAPEX and OPEX values decrease by 1.36% and 1.30%, respectively.
By comparing the percentage variation in Table 3, it can be concluded that the variations of TAC, annCAPEX, and OPEX are nonlinear.Also, it can be seen that the order of importance of the contributions of the process units is not influenced by the changes made on the H 2 product purity.
Finally, it is worth of mention that the (optimal) trends of the total electric power and membrane area are opposite, independently of the H 2 purity target level.This means that the higher the desired H 2 purity, the higher the demanded total electric power and the lower the required total membrane area; while the lower the H 2 purity, the lower the electric power and the higher the membrane area.4 shows the influence of the H 2 recovery variation on the optimal solution.The minimum TAC value obtained for a H 2 recovery of 89% is 1.738 M$•year −1 , which represents a reduction of 1.49% with respect to the obtained for 90%.An increment in the H 2 recovery of 0.01 leads to an increase in the TAC of 0.029 M$•year −1 which represents an increase of 1.64%.
Unlike the optimal trends observed for H 2 product purity variations, Table 4 shows that the requirements of electric power and membrane area follow the same behaviors, i.e., they both decrease with decreasing H 2 recovery levels, and vice versa.In addition, the variations of TAC, annCAPEX, and OPEX shown in Table 4 are not as nonlinear as those shown in Table 3.
The application of the here proposed model-based optimization approach has been successfully applied by the authors in other processes such as single-purpose seawater desalination plants [34][35][36], dual-purpose desalination plants [37,38], combined cycle power plants [39,40], amine-based CO 2 capture processes [41,42], biological wastewater treatment plants [43,44], real-time optimization in energy systems [45], liquid biofuel processors for H 2 production coupled to stationary fuel cells [46,47] and its associated heat exchanger network for optimal energy integration [48].For instance, based on process models, this approach allowed to find a new cost-effective configuration for multi-stage flash desalination processes [35,36], which differ from the conventional configuration in the locations of distillate, recycle, and discharge streams.In combined cycle power plants [39,40], the simultaneous optimization approach allowed to find a new layout for the heat recovery steam generator (HRSG), which differs from the conventional configuration in the way in which heat exchangers (economizers, evaporators, and superheaters) are interconnected between them.

Conclusions
This paper presented the optimization results of two-stage membrane processes for H 2 separation from a multi-component gas mixture by applying nonlinear mathematical programming approach for cost minimization.
A set of nonlinear algebraic equations obtained by discretization of a set of nonlinear differential ordinary equations, applying backward finite difference method (BFDM), which describe gas permeation through a polymeric membrane in a counter-current flow pattern was implemented in GAMS.
Expanders, compressors, vacuum pumps, heat exchangers, splitters, mixers, and membrane modules were appropriately interconnected to represent different two-stage membrane process configurations (embedded configurations) to find the structure that determines the minimum total annual cost, while meeting target H 2 purity and H 2 recovery specifications in the product (permeate) stream.
Firstly, when the proposed optimization model was solved without imposing any structural constraints, no expander was selected and the driving force in the membrane stages was created by combining feed/permeate compression and retentate under vacuum (optimal solution OPT).For cost model and design specifications considered in this work, a minimum total annual cost of about 1.764 M$•year −1 was obtained.Then, when constraints were imposed to prevent the use of vacuum in the permeate streams, the expander was again eliminated from the sub-optimal solution SUBOPT1, resulting in a total annual cost of about 2.038 M$•year −1 .This result clearly indicates that no benefit is obtained by including an expander to decrease the net electric power requirement.In other words, the investment required by the expander does not compensate the resulting decrease in the electric power requirement.This was confirmed by solving a third optimization problem, in which a constraint forcing the presence of the expander was considered.The total annual cost increased by about 24% x MS#,i,j mole fraction of the component i in the retentate stream of the membrane stage MS# at the discretization point j, dimensionless.
x MS#,R,i mole fraction of the component i in the retentate stream leaving the membrane stage MS#, dimensionless.
y MS#,i mole fraction of the component i in the permeate stream leaving the membrane stage MS#, dimensionless.
y MS1,i,j mole fraction of the component i in the permeate stream of the membrane stage MS# at the discretization point j, dimensionless.

Figure 1 .
Figure 1.Schematic of the studied two-stage membrane process for H2 separation.(MS represents membrane units, C compressors, HEX heat exchangers, VP vacuum pumps, EXP expansion turbine (expander), SP splitters, M mixers.PMS represents the permeate stream in MS, RMS retentate stream in MS, RRMS (internal) retentante recycle in MS, RRMS2_MS1 (external) retentate recycle from MS2 to MS1, R ## represents the fraction of R that is expanded in EXP and R # the fraction of R that is not expanded in EXP, and p H and p L indicate high and low pressure, respectively).

Figure 1 .
Figure 1.Schematic of the studied two-stage membrane process for H 2 separation.(MS represents membrane units, C compressors, HEX heat exchangers, VP vacuum pumps, EXP expansion turbine (expander), SP splitters, M mixers.P MS represents the permeate stream in MS, R MS retentate stream in MS, RR MS (internal) retentante recycle in MS, RR MS2_MS1 (external) retentate recycle from MS2 to MS1, R ## represents the fraction of R that is expanded in EXP and R # the fraction of R that is not expanded in EXP, and p H and p L indicate high and low pressure, respectively).

Figure 2 .
Figure 2. Schematic representation and nomenclature: (a) whole process; (b) membrane modules.(F indicates flow rate, Q heat load, T temperature, p pressure, x composition in feed/retentate streams, y composition in permeate streams, W power, A area of membrane and area of heat exchanger.The subscript i represents a mixture component, subscript j a discretization point along the membrane length, and subscript 0 the fresh feed).

Figure 2 .
Figure 2. Schematic representation and nomenclature: (a) whole process; (b) membrane modules.(Findicates flow rate, Q heat load, T temperature, p pressure, x composition in feed/retentate streams, y composition in permeate streams, W power, A area of membrane and area of heat exchanger.The subscript i represents a mixture component, subscript j a discretization point along the membrane length, and subscript 0 the fresh feed).
also shows that the optimal solution OPT requires 5701.7 m 2 of membrane area (5063.6 m 2 in the first stage and 638.1 m 2 in the second one) which, according to Figure 4b, represents 21.1% of the total equipment acquisition cost C INV (0.303 of 1.431 M$, see Table

Figure 3 .
Figure 3. Optimal process configuration obtained in solution (OPT), in which the driving force is created by combining compression and vacuum and no expander is included.

Figure 3 .
Figure 3. Optimal process configuration obtained in solution (OPT), in which the driving force is created by combining compression and vacuum and no expander is included.

Figure 4 .
Figure 4. Optimal cost distributions: (a) contributions of the total annualized capital expenditures (annCAPEX) and total operating expenditures (OPEX) to the total annual cost (TAC); (b) Equipment investment (CINV); (c) contributions of the electric power (EP) cost, membrane replacement (MR) cost, and cooling water (CW) cost to the total utility cost (CRM).

Figure 4 .
Figure 4. cost distributions: (a) contributions of the total annualized capital expenditures (annCAPEX) and total operating expenditures (OPEX) to the total annual cost (TAC); (b) Equipment investment (C INV ); (c) contributions of the electric power (EP) cost, membrane replacement (MR) cost, and cooling water (CW) cost to the total utility cost (C RM ).

Figure 5 .
Figure 5. Optimal mole fraction profiles of the components in the membrane modules obtained in solution OPT, in which the driving force for permeation is created by combining compression and vacuum and no expander is included for power recovery: (a) Permeate in MS1; (b) Permeate in MS2; (c) Retentate in MS1; (d) Retentate in MS2.

Figure 5 .
Figure 5. Optimal mole fraction profiles of the components in the membrane modules obtained in solution OPT, in which the driving force for permeation is created by combining compression and vacuum and no expander is included for power recovery: (a) Permeate in MS1; (b) Permeate in MS2; (c) Retentate in MS1; (d) Retentate in MS2.

Figure 6 .Figure 7 .
Figure 6.Sub-optimal solution SUBOPT1, in which the driving force for permeation is created by compression of the feed and permeate streams only (no vacuum) and no expander is included.

Figure 6 .
Figure 6.Sub-optimal solution SUBOPT1, in which the driving force for permeation is created by compression of the feed and permeate streams only (no vacuum) and no expander is included.

Figure 6 .Figure 7 .
Figure 6.Sub-optimal solution SUBOPT1, in which the driving force for permeation is created by compression of the feed and permeate streams only (no vacuum) and no expander is included.

Figure 7 .
Figure 7. Optimal mole fraction profiles of the components in the membrane modules obtained in solution SUBOPT1, in which the driving force for permeation is created by compression of the feed and permeate streams only (no vacuum) and no expander is included: (a) permeate in MS1; (b) permeate in MS2; (c) retentate in MS1; (d) retentate in MS2.

Figure 8 .
Figure 8. Sub-optimal solution SUBOPT2, in which the driving force for permeation is created by combining compression and vacuum and the expander is included for power recovery.

Figure 8 .
Figure 8. Sub-optimal solution SUBOPT2, in which the driving force for permeation is created by combining compression and vacuum and the expander is included for power recovery.

Figure 8 .
Figure 8. Sub-optimal solution SUBOPT2, in which the driving force for permeation is created by combining compression and vacuum and the expander is included for power recovery.

Figure 9 .
Figure 9. Optimal mole fraction profiles of the components in the membrane modules obtained in solution SUBOPT2, in which the driving force for permeation is created by combining compression and vacuum and the expander is included for power recovery: (a) permeate in MS1; (b) permeate in MS2; (c) retentate in MS1; (d) retentate in MS2.
the compressor C# associated with the membrane stage MS#, K. T MS# operating temperature in the membrane stage MS#, K. T out HEX# outlet temperature from the heat exchanger HEX#, K. W C# power required by the compressor C# associated with the membrane stage MS#, MW.W VP# power required by the vacuum pump VP# in the membrane stage MS#, MW.W EXP power recovered in the expander EXP, MW. x i,0 mole fraction of component i in the feed stream, dimensionless.x MS#,i inlet composition of the component i in the membrane stage MS#, dimensionless.

Table 2 .
Comparison of the optimal cost values of the three solutions OPT, SUBOPT1, and SUBOPT2.