Membrane separation study for methane-hydrogen gas mixtures by molecular simulations

Direct simulation results for stationary gas transport through pure silica zeolite membranes (MFI, LTA and DDR types) are presented using a hybrid, non-equilibrium molecular dynamics simulation methodology introduced recently. The intermolecular potential models for the investigated CH$_{4}$ and H$_{2}$ gases were taken from literature. For different zeolites, the same atomic (Si and O) interaction parameters were used, and the membranes were constructed according to their real (MFI, LTA, or DDR) crystal structures. A realistic nature of the applied potential parameters was tested by performing equilibrium adsorption simulations and by comparing the calculated results with the data of experimental adsorption isotherms. The results of transport simulations carried out at 25$^0$C and 125$^0$C, and at 2.5, 5 or 10 bar clearly show that the permeation selectivities of CH$_{4}$ are higher than the corresponding permeability ratios of pure components, and significantly differ from the equilibrium selectivities in mixture adsorptions. We experienced a transport selectivity in favor of CH$_{4}$ in only one case. A large discrepancy between different types of selectivity data can be attributed to dissimilar mobilities of the components in a membrane, their dependence on the loading of a membrane, and the unlike adsorption preferences of the gas molecules.


Introduction
Zeolites are made up of silicon, aluminum, and oxygen atoms linked together so that they form structurally well-defined pores. A high regularity of the structure distinguishes zeolites from other porous materials and makes their high selectivity in catalytic and separation processes possible [1,2]. A usual size of zeolite micropores is similar to that of many small molecules. Therefore, it is possible for some molecules to enter the zeolite and then get stuck in the pores, where they can even react with each other, while the other molecules can move through the zeolite channels faster. In separation processes, the advantage of zeolites over other types of adsorbents/membranes is that they typically offer good endurance to high pressure and temperature and can often tolerate harsh chemical environments. The Si-Al ratio of the zeolite framework is an important factor of the application. Zeolites with lower Si-Al ratios are more hydrophilic, while zeolite membranes with higher Si-Al ratios have fewer structural defects. High-silica zeolite membranes are greatly preferred in gas separation applications, and among them silicalite (MFI) is one of the most commonly studied.
Purification of methane from carbon dioxide (carbon dioxide is one of the major contaminants of natural gas) turns out to be an especially attractive gas separation process. Another important practical system that is involved in the process of purification of synthetic gas obtained from steam reforming of natural gas is the methane-hydrogen system. The available literature data, however, are scarce for all the systems that combine methane-hydrogen gas mixtures with the use of zeolite membranes [3,4].
The behavior of material systems, characterized at the atomic level, can be effectively studied using classical molecular simulations. The atomistic simulation results can explain or, in some cases, substitute the experimental results. Nowadays, these simulations play an important role in the description of processes that occur in crystalline adsorbents and membranes [5][6][7][8][9][10][11][12][13]. To simulate steady-state transport of molecules, we need a method that satisfies two criteria: it must ensure the real dynamics of a system, and it must maintain the steady-state driving force at the microscopic level. The dynamics of a system can be investigated microscopically by means of molecular dynamics (MD), dynamic Monte Carlo [8] and other direct or indirect simulation methods (e.g., a very fast NP+LEMC method [11]). To preserve a constant driving force, these methods are often linked with the other techniques, e.g., dual control volume or local equilibrium Monte Carlo (LEMC) techniques [12]. There are many composite methods having many limits of applicability and numerous advantages and disadvantages, such as the gradient relaxation MD [13] or other earlier gradient techniques (e.g., [14,15]), a self-adjusting plates technique [16], the external field MD [17] or its boundary-driven version [18], and the dual control volume grand canonical molecular dynamics (DCV-GCMD) [5]. In general, the main problem with simulating steady-state membrane transport at the atomic level is that we cannot simulate (at least virtually) bulk fluid phases of realistic size. Therefore, particle depletion/accumulation in the simulated bulk regions frequently occurs, if particle reinjection/removal steps are not applied. A sudden appearance and annihilation of molecules, however, can disturb the steady-state flux of transporting particles. Recently, we developed a simple atomistic simulation method for membrane transport that can maintain a driving force without significantly disturbing the previously-developed steady-state flux [19], while it properly mimics the common experimental situation in gas permeations measurements, where pressure is the main control parameter.
In this work, direct simulation results for steady-state gas transport through some relevant pure silica zeolite membranes are reported using our novel hybrid MD simulation method. In what follows, we briefly outline the applied transport simulation technique, and then specifications of the performed simulations and results for the separation of the technologically important CH 4 -H 2 mixture (that might also be relevant in the development of engine fuels with high hydrogen content) will be presented in more detail.

Transport simulation method
The idea behind our method is based on the fact that in experimental gas permeation arrangements, pressure is the property that can be controlled relatively simply: the partial pressure for each component on the input side of the membrane and the total pressure on the permeate side. The essential point of our atomistic simulation scheme (called pressure-tuned, boundary driven molecular dynamics technique, PBD-MD) is that pressure is controlled by adjusting the particle numbers on the two sides of the membrane indirectly [19].
We started from the DCV-GCMD approach [5] containing two equilibrium control cells distinguished by unequal intensive thermodynamic parameters (such as chemical potential or pressure) on the two sides of the membrane, where random particle insertion/deletion moves are applied to maintain bulk phase properties in the control cells and thus to uphold the desired (constant) driving force through the membrane. However, to avoid such non-physical particle moves in the vicinity of the membrane affecting any eventual sorption layers on the membrane surfaces and interfering with earlier stabilized flux of the permeating particles, adjusting the particle numbers is restricted to zones far from the membrane region. This means that while all properties calculated and monitored in the control cells predominantly stem from the movements and interactions of particles in these regions, the artificial perturbation of the system is only present in the boundary regions of the simulation cell [18]. Our earlier test clearly showed that the particles inserted into the system have no "memories" of their initial velocities before they reach the interaction range of the membrane [19].
Here, we consider the system as being at a constant temperature T with a fixed number of particles N and a box volume V. To attain the target pressure p (or partial pressures in the case of mixtures), regular perturbations in the number of particles are applied close to the boundaries of the total simulation box. In this way, as the chemical potentials of the control cells used in the traditional DCV-GCMD technique are connected to the corresponding partial pressures, the pressure can be controlled effectively in the control cells. In our trial-and-error type pressure-tuning approach, the particle insertion and deletion steps are performed randomly (i.e., randomly chosen particles and positions) at both ends of the simulation box, in the direction of the transport, far enough from the membrane region (in this respect, this technique is similar to the Particle Counting (PACO) method of Berti et al. [20]). One particle insertion (+1 case) or deletion (−1 case) is executed periodically if the following general criterion is satisfied: where N control cell is the number of particles in the control cell, and p control cell and p target are the calculated pressure and the designated pressure, respectively, in the same control cell. Here, p denotes partial pressure on the feed side or total pressure on the permeate side, and N control cell denotes the number of particles of the individual components on the feed side or the total number of particles on the permeate side. The attempts are always accepted, except for those random insertion steps which result in particle overlap (in such cases, the insertion attempt is repeated).

Simulation details
In the transport simulations, standard MD was used employing the leap-frog algorithm as an integrator with the time step of 2 fs. The Berendsen thermostat with the thermal coupling parameter of 200 fs (very weak coupling) was applied for controlling the temperature. The temperature was defined by subtracting off the streaming velocity from the x component of the particle velocities, where +x direction is the direction of the transport. The simulation box was confined by impenetrable, soft repulsive walls in the direction of the transport, while periodic boundary conditions were applied in the other two directions. (Note that the use of periodic boundary conditions in the direction of the transport is incompatible with the applied particle number adjustment procedure.) The geometry of this simulation box is depicted in figure 1. The boundary regions were considerably wider than the range of the repulsive walls on the two sides of the box. Random initial velocities were assigned to the inserted particles according to the Maxwell-Boltzmann distribution at the prescribed temperature. One particle insertion or deletion step was performed periodically in both boundary (wall) regions, outside the range of the repulsive walls. The length of the simulation period between two particle number modifications was taken as 2000 consecutive MD steps and N control cell and p control cell were collected as averages for these periods. The average pressure values were calculated by the virial expression. The membrane transport processes were simulated for at least 30 million MD steps, but in several cases twice as long runs were needed to collect a sufficient amount of transfer events (at least several hundred).
We also executed equilibrium adsorption simulations in the grand canonical ensemble (fixed chemical potential, volume and temperature) using the standard grand canonical Monte Carlo (GCMC) technique. The pressure/partial pressure of the adsorbate molecules in the gas phase was given indirectly by specifying the chemical potential of the component. The simulation box size was equal to the size of the investigated zeolite crystal and periodic boundary conditions were applied in all three spatial directions.

23002-3
The external surface adsorption, therefore, was not taken into account. The adsorption simulations were conducted for 200 million Monte Carlo steps.
We present results for the adsorption (equilibrium state) and transport (steady-state) of gases on allsilica MFI, LTA and DDR zeolites at T = 298.15 and 398.15 K, and at p = 250 and 500 (or 1000) kPa. These are built up from SiO 4 tetrahedrons, with linear and/or zig-zag (sinusoidal) channels [21]. In the MFI zeolite [figure 2 (a), (b)], straight channels of elliptical cross-section extending along one direction are cross-linked by zig-zag channels of nearly circular cross-section extending in the other perpendicular direction (in both cases the channel diameters are about 0.7 nm). The crystal lattice parameters of its orthorhombic unit cell are a = 2.0090 nm, b = 1.9738 nm and c = 1.3142 nm, with all the lattice angles being 90°. The pure silicon type of the LTA zeolite [figure 2 (c), (d)] has a cubic structure (all the lattice angles are 90°) with nearly spherical sodalite cages and straight channels with a diameter of about 1.1 nm between them. The crystal lattice parameter of its unit cell is a = 1.1919 nm. The pure silica DDR zeolite [figure 2 (e), (f)] has a two-dimensional pore network with channels of slightly elliptical cross-section and with the characteristic channel diameter of about 0.8 nm. The crystal lattice parameters of its hexagonal unit cell are a = 1.
The membrane structures were constructed according to the available crystallographic information from the IZA database [21] and taken to be defect-free. For the adsorption simulations, we built the model structures in the 1 × 1 × 2, 2 × 2 × 2 and 2 sin γ × 2 × 1 arrangements of the unit cells with MFI, LTA and DDR adsorbents, respectively. In the transport simulations, 1 or 2 unit-cell thick and 10-40 unit-cell wide (tall) membranes were used with the intention that the membrane sizes of the different zeolites should be as similar as possible. Two crystal orientations were realized in the membranes to place the (1) straight or (2) zig-zag channels of the zeolites in parallel with the x direction [in LTA membranes the straight channels were rotated by 45°for case (2)]. Thicknesses and surface sizes of these membranes containing 5−8000 Si and O atoms are summarized in table 1. The particles of the zeolite lattices were fixed at the crystallographic positions, so the frameworks were kept rigid [22]. In the steady-state transport simulations, the box size in the direction of the transport was set to 80 nm. To describe the molecular interactions, the shifted and cut Lennard-Jones (LJ) pair potential was employed together with the Lorentz-Berthelot combining rule. The effective interaction potential for a pair of particles (α and β) was calculated as where is the standard 12-6 LJ potential, u αβ c is the value of the potential at the cut-off distance r αβ c , and σ and are the size and energy parameters of the LJ potential, respectively. The spherical cut-off distance was set to r αβ c = 3.5σ αβ leading to u αβ c = −0.00217478 αβ . All particles interacted with the walls according to the Weeks-Chandler-Anderson potential. This is also a shifted and cut LJ potential but it has a short cut-off 2 1/6 σ wall , and thus its interaction can only be repulsive. This particle-wall potential was chosen because it is convenient to be used in MD simulations. Its choice does not influence the results in the membrane region.
Single-site (i.e., not atomically detailed) models were used for the gas molecules, and the interaction parameters were taken from the literature [23,24]. In addition, the same atomic interaction parameters were applied for different zeolites (see table 2): literature σ and parameters [23,25] of the zeolitic Si and O atoms were slightly adjusted to better reproduce the available equilibrium adsorption data for MFI.

Results
In the adsorption simulations, the equilibrium selectivity was defined as follows: Here, q denotes the loading of the zeolite (adsorption inside the adsorbent in mol · kg −1 , relative to the mass of the adsorbent). In the transport simulations, the permeation selectivity (or dynamical selectivity) Atom/molecule σ / nm ( /k B ) / K O [23] 0.270 130.0 Si [25] 0.070 20.0 CH 4 [23] 0.373 147.9 H 2 [24] 0.296 36.7 Soft repulsive wall 0.300 120.0 was calculated from where j is the component flux. For comparison, the permeation ratio (R P ) was also evaluated; R P is formally similar to S P but here the pure component fluxes are used, and thus it can be regarded as the idealized limiting case. First, the applied size and energy parameters of the potential models were tested with equilibrium adsorption calculations for CH 4 and H 2 . Single-component adsorption isotherms were determined by GCMC simulations to select/verify the potential parameters for the pure gas components and the zeolite atoms. Here, we present the results obtained with the final parameters shown in table 2.
In the case of MFI, as expected, the calculated loadings with CH 4 agree quite well with the experimental adsorption isotherm data [26] in the investigated temperature and pressure ranges [ figure 3 (a)]. At lower pressures, a relatively good agreement with the available experimental data [27] can be also observed with H 2 [ figure 3 (b)]. There are some inaccessible cavities in this zeolite for the gas molecules, which turned out to be effectively inaccessible in the equilibrium simulations. For the pure silica LTA zeolite, we found only one experimental adsorption isotherm at room temperature or above, and only with CH 4 , and its data [28] were systematically overestimated in the simulations [figure 3 (c), (d)]. However, when we artificially prevented the creations of CH 4 molecules inside the sodalite cages, an acceptable matching between simulation and measurement could be detected. In such a way we took into account the physical diffusion pathways in the zeolite, since the CH 4 molecules are incapable of passing through the small windows of the sodalite cage. For the pure silica DDR zeolite, experimental data with CH 4 were available at two temperatures [29] [figure 3 (e), (f)] and we did not find such data with pure H 2 . As it is seen, the simulated adsorption isotherms with CH 4 deviate considerably from their experimental counterparts. In our grand canonical simulations, however, this zeolite contained a large quantity of occupiable pores or cavities, which have, according to the literature [30,31], too small windows to allow molecule access, and so these simulations also overestimate the proper results (in these simulations, due to the complicated geometry, we could not preselect the regions where artificial molecule creations are forbidden). Accessible volume data from the literature make possible a rough estimation of an average (temperature, pressure and composition-independent) correction factor for the calculated adsorption loading of DDR: this amounts to 0.6 or 0.7 with CH 4 or H 2 , respectively. Notwithstanding all these, we must acknowledge that we had to make a compromise at this point to keep the transferability of the pair potential parameters and highlight the differences in the simulated results that originate from the structural alterations between the zeolites.
Mixture adsorption was studied using GCMC simulations with an equimolar CH 4 -H 2 gas system at temperatures of 298.15 and 398.15 K, and at the total pressures of 250 and 1000 kPa (more accurately, the CH 4 -H 2 mixture was a nearly equimolar mixture, as we set the partial pressures of the components to be equal). The obtained equilibrium loading and selectivity values are shown in figure 4. The S E values are between ∼4 and ∼80, so the degree of adsorption is noticeably higher for CH 4 than for H 2 . We found the absolute loadings with H 2 to be rather low (not shown), and this suggests some exclusion effect against the smaller gas molecules. The relatively larger S E values for the model DDR zeolite are due to its higher loadings with CH 4 (even if we consider the necessary corrections with the ratio of its inaccessible cavities). Generally, the calculated selectivity decreases with an increasing temperature and pressure.
Our simulations using the PBD-MD technique delivered results for binary gas transport through the zeolite membranes with the CH 4   feed side pressures of 250 and 1000 kPa. The permeate side pressure was always equal to 100 kPa (approximately ambient condition). As a reference, single-component transport simulations were also performed, where the feed side pressure was equal to either the partial pressure of the component or to the total feed side pressure of the corresponding mixture simulation. Before we proceed with the permeation results, we give a validation of our PBD-MD approach for the present conditions. The concentration profile of the transporting particles is a good indicator whether the control cells can be considered as the ones representing bulk phases or not. For a few simulations  . It should be noted, furthermore, that the calculated total pressure on the permeate side was always correct in the simulations (the target pressures were generally reproduced within 1%), even if the partial pressures and component concentrations looked different here [e.g., figure 5 (c)]. This corresponds to the selectivity-related properties of the membrane. The calculated permeation selectivities and permeation ratios are depicted in figure 6. It is seen that the S P values are mostly less than 1, which means that, in their majority, the model membranes are selective for the transport of H 2 . Sufficiently high or low selectivity data for practical use are absent in the figure. Separation of CH 4 from H 2 on the permeate side is only reliable with the zig-zag model MFI membrane at a lower temperature and higher pressure, where the dynamical selectivity is almost equal to 2. On the other hand, separation of H 2 from CH 4 is preferred with the model DDR membranes, especially at a higher temperature and pressure (S P < 0.25). In this respect, the model LTA membranes are of intermediate type. The calculated S P values are consistently larger than the corresponding R P values, indicating the presence of momentum coupling between the unlike gas molecules during the mixture transport and, more specifically, the effect of the prevalent adsorption of CH 4    the most remarkable fact is that the permeation selectivity values are much lower than the equilibrium selectivity values. It means that the adsorption preference for CH 4 observed in all equilibrium cases does not predictably go hand in hand with transport selectivity in favor of CH 4 . This remains true even if we know that the feed side pressure in the dynamic simulations cannot be unequivocally compared to the equilibrium pressure of the adsorption simulations. In the light of the pressure gradient applied in the transport simulations across the membranes, however, it is somewhat logical that the equilibrium adsorption loadings are generally larger than the dynamic loadings of the corresponding membranes with p feed = p adsorption (compare data of table 3).
The S P values predominantly exhibit some decrease with an increasing temperature, which is attributed to a smaller decrease of the flux of H 2 with temperature (cf. table 3). At the same time, the general declining pressure dependence of S P observed here is due to a more intense increase of the flux of H 2 with pressure. The only exception to this latter tendency is the transport through the model MFI membranes (with both straight and zig-zag channels) at a lower temperature, where the growth in the flux of CH 4 is larger. This behavior cannot be explained by anomalous loadings of the MFI membrane in dynamic situations. As data of table 3 show, the CH 4 to H 2 loading ratios of the membranes in the dynamic simulations roughly follow the trends of the S E values obtained from the equilibrium adsorption simulations. However, the much milder temperature and pressure dependences of the R P values than those of the S P values (see figure 6) underline the importance of the momentum coupling between different components connected to some adsorption preferences.

Discussion and conclusions
We used the novel PBD-MD simulation technique to study the stationary membrane transport at the molecular level. The method makes possible an accurate determination of steady-state fluxes of gases 23002-10 through microporous membranes. We investigated the permeation of CH 4 and H 2 gases across pure silica zeolite membranes (MFI, LTA and DDR types) at 25°C and 125°C and at 2.5 and 5 (or 10) bar.
The applied shifted and cut LJ potential is a less detailed interaction model for the simulations, and certainly, we could have better reproduced the available experimental data for the adsorption of pure CH 4 and H 2 gases by using more realistic all-atom potential models involving sophisticated inter-and intramolecular interactions. By this simple modelling, we rather preserved the transferability of the zeolite parameters and stressed the structural differences of the investigated zeolites. On the other hand, in such a way the total computational need of the executed transport simulations could be kept at a relatively low level (a couple of months with some dozens of CPU cores in a cluster computer).
We did not find experimental data for the investigated transport systems. At the same time, there would be many difficulties with the comparison of experimental and such atomistic simulation results. Experimental (synthesized) zeolite membranes typically consist of microcrystals with random orientations, and consequently, they exhibit an anisotropic pore geometry. The permeation properties of these membranes can be very different from those with oriented microcrystals [32] (these latter can be prepared in exceptional synthesis conditions). Moreover, real transport processes in zeolite membranes can occur in the intercrystalline pores. Other common problems include the presence of defects in the microcrystals and contaminations on the membrane surfaces, which may significantly affect the flux and selectivity of the membranes in experiments. Furthermore, an experimental zeolite membrane is built up from a thin zeolite layer and a macroporous supporting layer (this provides the membrane with the requisite mechanical strength); the support resistance may also influence the separation performance. From the simulation viewpoint, the main obstacle is that the extremely high demand of computational power today prevents the atomistic simulation of systems with real-life size (or nearly real-life size) membrane thicknesses (which is in the order of a micrometer). According to our former experience, the dependence of the transport properties on membrane thickness calculated with commonly applied simulation system sizes (with 1, 2, 3, . . . , unit-cell thick membranes) cannot be easily used to forecast the magnitude of the scale-up effect, since this dependence is not necessarily linear at the atomic scale. For these reasons, we consider our present calculations only as one step towards model studies of the interplay of adsorption and diffusion in the permeation processes with realistic microporous membrane structures.
The results for the investigated systems and conditions have two striking features. First, there are no profound differences between the three microporous zeolite structures of similar channel diameter concerning both their equilibrium and transport selectivities. Second, the most important trend that is revealed from the present molecular simulations is that these model zeolites suffer very low transport selectivities for CH 4 over H 2 as compared to the equilibrium ones. Higher equilibrium adsorption of CH 4 observed with these zeolites is more frequently linked with the permeation selectivity for H 2 over CH 4 . The contrast between the equilibrium and transport selectivity data can be ascribed to the quite different adsorption preferences of the gas molecules, the strongly dissimilar mobilities of the components in the membrane and their dependence on the loading of the membrane. Here, a lower mobility of one gas component is compensated by a more intense adsorption capability, and the reverse situation occurs with the other component, leading to very similar calculated component fluxes, and consequently, to permeation selectivity values not far from 1. The manifestation of this compensating effect is most remarkable in the case of the model DDR zeolite, where one of the largest equilibrium loadings of the adsorbent is combined with the significantly lower fluxes of the membranes. A certain role of the different degree of temporary sorption of molecules in dynamic conditions is evidenced by the fact that for each zeolite, somewhat higher permeation selectivities were obtained for the zig-zag channels of longer transfer route than for the straight channels, which have an otherwise similar diameter. Note that the applied pair of the LTA model membranes is a special case, where the transfer route of one of these membranes is artificially lengthened through rotating the straight channels by 45°. At the same time, the latter observation with the LTA membranes suggests that the use of larger membrane thicknesses would probably result in the increase of permeation selectivity. We obtained a truly higher permeation rate of CH 4 in only one case: for the model MFI membrane with zig-zag channels in the direction of the transport at a lower temperature and higher pressure. This prevalent transport of CH 4 through the membrane seems to be a favorable interplay between a relatively high equilibrium selectivity and a smaller channel diameter with an optimal channel structure/orientation, which latter impedes the diffusion of H 2 molecules to a greater extent. In the present systems, the transport through the membranes was more often easier for the  far faster moving and smaller H 2 molecules. To this end, however, there was one condition to be fulfilled: the absence of a significant steric hindrance from the adsorbed or slowly moving CH 4 molecules. In our case, apparently, the mobility of the H 2 molecules was not noticeably reduced due to such a blocking effect.