Multiscale simulation of gas transport in mixed-matrix membranes with interfacial polymer rigidification

Journal Pre-proof Multiscale simulation of gas transport in mixed-matrix membranes with interfacial polymer rigidification Gloria M. Monsalve-Bravo, Ravi C. Dutta, Suresh K. Bhatia PII:

S1387-1811(19)30841-8

DOI:

https://doi.org/10.1016/j.micromeso.2019.109982

Reference:

MICMAT 109982

To appear in:

Microporous and Mesoporous Materials

Received Date: 20 November 2019 Revised Date:

20 December 2019

Accepted Date: 20 December 2019

Please cite this article as: G.M. Monsalve-Bravo, R.C. Dutta, S.K. Bhatia, Multiscale simulation of gas transport in mixed-matrix membranes with interfacial polymer rigidification, Microporous and Mesoporous Materials (2020), doi: https://doi.org/10.1016/j.micromeso.2019.109982. This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain. © 2019 Published by Elsevier Inc.

Credit Author Statement Gloria M. Monsalve-Bravo: Conceptualization, Methodology, COMSOL Simulations, Writing- Original draft preparation. Ravi C. Dutta: MD, methodology, Writing. Suresh K. Bhatia: Conceptualisation, Supervision, Writing- Reviewing and Editing,

Multiscale Simulation of Gas Transport in Mixed-Matrix Membranes with Interfacial Polymer Rigidification Gloria M. Monsalve-Bravo, Ravi C. Dutta, Suresh K. Bhatia* School of Chemical Engineering, The University of Queensland, Brisbane QLD 4072, Australia ABSTRACT We investigate permeation in non-ideal mixed-matrix membranes (MMMs) having interfacial defects based on a multiscale simulation approach, involving sequential nanoscale and macroscale simulations. Our approach combines insight from equilibrium molecular dynamics (EMD) simulation of the transport in the individual phases, and in the interfacial region, with macroscopic simulation of the transport through solution of coupled 3-D partial differential equations (PDEs) via a finite-element method (FEM). Thus, unlike existing adaptations of classical effective medium theory (EMT), our approach avoids the empirical fitting of filler-polymer interfacial properties, such as interfacial thickness and permeability, against experimental permeation data. While the proposed multiscale simulation approach is general and applicable to any MMM system, it is illustrated in the context of single CO 2 and CH 4 permeation in two different polyimide-based MMMs; with one using MFI-type zeolite as filler phase and the other using carbon molecular sieves (CMS). For both systems, our nanoscale simulations reveal an interfacial region with rigidification of the polymer and reduced permeability at the polymer-adsorbent interface. Comparison of permeability predictions of classical EMT-based models considering interfacial rigidification, with our macroscale simulations, demonstrates the former to be best at low particle loadings. Further, our macroscale simulations reveal the existence of an optimum filler particle size, never discussed before, which maximizes the permeability in non-ideal MMMs; this arises from the competing effects of decrease in permeability with increase in filler size and decrease of the relative interfacial resistance with increase of filler particle size. Moreover, the CO 2 CH 4 perm-selectivity is found strongly sensitive to the interfacial resistance for both finite uniform and non-uniform particle size distributions, an effect not captured by existing EMT models and overlooked in prior fundamental approaches. Keywords: Mixed-matrix membrane; non-ideal interface; polymer rigidification; membrane simulation; permeability. *Corresponding Author: E-mail: [email protected]

Page 1 of 41

1. Introduction Membranes for gas separation are increasingly used in industrial applications due to their great potential in energy and cost savings when compared to conventional gas separation techniques [1– 3]. Commonly, membranes are classified as organic (polymer) or inorganic based on their material [4], of which the performance of polymer membranes is often limited by a trade-off relation between the selectivity and permeability [5], while implementation of inorganic membranes in practical applications is generally limited by their fragility and higher cost over polymer membranes [4]. Mixed-matrix membranes (MMMs), consisting of a selective inorganic filler phase embedded in a continuous polymer matrix, have emerged as alternative to overcome these disadvantages [2]. These composite membranes combine the high intrinsic permeability and separation efficiency of advanced molecular sieving and/or nanoscale materials with the robust processing capabilities and mechanical properties of polymers [2,6]. However, success of MMMs is modest to date due to challenges associated with the membrane formation that lead to non-ideal membrane structures [1,3]. Over the past decade, much effort has been devoted to the screening and selection of suitable filler and polymer combination to achieve near ideal MMM structures [2–4,6,7]. However, the transport through both ideal and non-ideal MMMs remains much less investigated [8–10]. Existing permeation models, based on well-established effective medium theory (EMT) [11–13], always assume a uniform concentration field across the MMM; an assumption that has been recently shown inadequate for ideal MMMs operating at typical practice pressures [14–17]. Some efforts have also been devoted to advance EMT to accommodate non-ideal interfacial morphologies in the form of a rigidified polymer layer [18–22], partially blocked filler layer [23,24] and/or homogeneous voided layer [18,19,25,26]. However, these new models share the uniform field assumption of classical EMT, which influences the perceived effect of these non-ideal interfacial morphologies on the permeability [8]. Thus, while existing non-ideal permeation models are often shown to fit well single gas permeabilities [23,26–28], they lead to interface thicknesses of about ( 20 − 200 nm ) [21–

Page 2 of 41

23,27,29], that are one or two orders of magnitude greater than those recently found through equilibrium molecular dynamics (EMD) simulation of MMMs ( 1 − 2 nm ) [30–32]. As an alternative to predictions based on EMT models, there has been a growing interest in rigorously simulating the gas transport through MMMs [14–16,33–36]. This simulation approach is based on numerical solution of coupled three-dimensional (3-D) partial differential equations (PDEs) describing the transport through the composite membrane via the finite-element method (FEM), with recent advances and current challenges in the use of this simulation technique for predicting gas permeability in MMMs reviewed elsewhere [8]. FEM simulations have been shown to accommodate non-uniformities in the concentration field and gas permeability due to variation of intrinsic system properties such as finite particle size [14,34], particle shape [35,36], isotherm nonlinearity [15], and membrane geometry [16,33]. However, the approach has yet to be used to predict the effect non-ideal filler polymer interfacial morphologies on the gas permeability. While some progress in this direction has recently been made for MMMs with spherical particles [37], gas diffusion at the interface is still considered concentration-independent, with the interface thickness fitted upon comparison with experimental permeation data in MMMs while considering the pressure gradient across the MMM of about 0.2 − 0.3 bar . Under this conditions, nonlinear isotherm-related effects are insignificant and, therefore, similar to EMT models, current simulation predictions are only satisfactory for MMMs operating at low pressure, for which the Henry’s law limit for gas adsorption is well accepted. In the above fundamental approaches, the gas properties (e.g. permeability, diffusivity, and solubility) in the MMM individual phases are generally based on single-gas experimental permeation data in pure polymer and inorganic membranes [9,17,27,35,38]. However, gas permeabilities in filler-based inorganic membranes often differ from one study to another, due to challenges associated with such membrane fabrication, in which intergrowth of crystals is necessary to avoid large pores. Thus, for example, reported transport properties of several gases in ZIF-8 membranes vary significantly [39–42], and can differ from the intrinsic transport coefficient

Page 3 of 41

pertinent to the filler particles in MMMs; as recently demonstrated by using pulsed field gradient (PFG) NMR to measure light gas diffusivities inside ZIF-8 particle beds and ZIF-8-based MMM [43–45]. As an alternative, molecular dynamics (MD) simulations have recently been adopted to predict the gas transport through individual inorganic (e.g. MOFs and ZIFs) particles and polymer materials [30,46–49], and subsequently these MD-based permeabilities are used in EMT to calculate the gas permeability in MMMs, considering an ideal polymer-filler interface and Henry law isotherms (i.e. permeabilities are concentration-independent, as required by conventional EMT) [50,51]. Thus, the implementation of MD simulations in conjunction with rigorous macroscale FEM simulations, while considering interfacial non-ideality and accounting for concentration dependence of individual phase properties, including those in the filler-polymer interface, and their effect on the overall gas permeability in the MMM, is therefore a key advance, which forms the subject of this work. To the above end, we introduce a novel multiscale simulation approach to investigate gas transport through both ideal and non-ideal MMMs, integrating insight gained from EMD simulations with rigorous macroscale simulations of the transport through MMMs. Here, nanoscale Monte Carlo and EMD simulations are used to characterize gas adsorption isotherms and diffusivities in the MMM individual regions, including the interface region. These EMD-based region-specific properties are subsequently used to implement the FEM to simultaneously solve 3-D PDEs describing the transport in the MMM as a whole. Such full-scale membrane system, with finite thickness in the permeation direction, comprises spherical filler particles surrounded by a homogeneous rigidified layer, randomly located in a polymer matrix. In this way, our macroscale simulations explicitly consider effects of the filler-polymer interfacial resistance, finite filler particle size and its distribution together with isotherm nonlinearity in the MMM constituent phases, never investigated before together through a fundamental approach. As an immediate consequence, our FEM approach also allows us to examine the applicability of existing pseudo-two-phase models of MMMs that consider interfacial non-ideality. We note that while our multiscale approach enables

Page 4 of 41

first principles-based in silico design of MMMs, without necessitating macroscopic phase-specific permeation measurements, experimental microscopic diffusion measurements from PFG NMR [43– 45] may also be used in conjunction with the macroscopic simulation to explore the overall gas transport through the MMM. In this work, we illustrate the application of our multiscale simulation technique for the permeation of pure CO 2 and CH 4 in a MMM comprising of MFI-type zeolite as filler phase and polyimide (PI) as continuous phase. For this system, individual phase properties are fully based on EMD simulations, without the need of empirically fitting the interface properties (thickness and permeability), as commonly done in practice [21–23,27,52]. Further, we also implement our macroscale simulations to investigate a MMM comprising carbon molecular sieves (CMS) as filler phase and Matrimid®5218 as continuous phase, in which thickness of the interfacial region is calculated from EMD simulations and individual phase properties are calculated using available experimental data of the system. For both systems, we compare predictions of our simulations to those based on popular permeation models [18,20,25,53]. Our results demonstrate that existing pseudo-two-phase models of MMMs with interfacial defects are satisfactory only at low filler particle loadings (φ ≤ 0.2) , and that the thickness of the interfacial region based on fits using these models is considerably overpredicted. 2. Modeling permeation in non-ideal mixed-matrix membranes The concept of a non-ideal ideal MMM is often associated with a three-phase composite comprising a dispersed filler phase, a continuous polymer matrix and the interface region between these two phases [8,9,27]. To date, a myriad of models have been proposed based on this concept [18–20,23,24,26,54,55], with all derived either as (i ) the solution of the transport problem in a composite comprising spherical core-shell inclusions embedded in a continuous matrix [18,28] or

(ii) the assumption that the filler-interface system can be idealized as a pseudo-spherical particle embedded in the continuous matrix by applying classical EMT-based models several times [19,20,22,25,29,53,54,56]. Of these, the Felske model [18] of the first group and the pseudo-twoPage 5 of 41

phase (PTP) Maxwell [25], Bruggeman [20] and Chiew-Glandt [53] models of the second group remain popular for estimating the interface properties [19,22,26,27]. Table 1 summarizes the key equations of these models, in which P and φ denote the permeability and volume fraction, respectively, and superscripts/subscripts f , c , i , m and g denote the filler phase, continuous phase, interface, MMM, and combined filler phase and interface composite, respectively. Table 1. Summary of popular models for permeation in non-ideal mixed-matrix membranes. Model Key equations  2(1 − φgc ) + (1 + 2φgc )(η γ )  Pm = Pc  Felske model [18]   (2 + φgc ) + (1 − φgc )(η γ ) 

PTP Maxwell model [25]

 P + 2 Pc − 2φgc ( Pc − Pg )  Pm = Pc  g   Pg + 2 Pc + φgc ( Pc − Pg ) 

(2)

 P + 2 Pi − 2φ fi ( Pi − Pf )  Pg = Pi  f   Pf + 2 Pi + φ fi ( Pi − Pf ) 

(3)

1 3

PTP Bruggeman model [20]

PTP Chiew-Glandt [53]

(1)

−1  Pm   ( Pg Pc ) − 1   = 1 − φgc      Pc   ( Pg Pc ) − ( Pm Pc ) 

(4)

1

−1  Pg  3  ( Pf Pi ) − 1   = 1 − φ fi      Pi   ( Pf Pi ) − ( Pg Pi ) 

(5)

1 + 2β gcφgc + (κ 2gc − 3β gc2 )φgc2  Pm = Pc   1 − β gcφgc  

(6)

1 + 2β fiφ fi + (κ 2fi − 3β f2i )φ fi2  Pg = Pi   1 − β fiφ fi  

(7)

Among the models in Table 1, that of Felske in Eq. (1) corresponds to the exact first order solution of the transport problem in a dispersed composite comprising non-interacting core-shell particles surrounded by a continuous matrix [8,9]. The Felske model is based on the Laplace solution for the potential of a dispersion similar to the Maxwell model [11]. Thus, Eq. (1) is only applicable to diluted suspensions (φ gc ≤ 0.2) [28,56], in which φ gc = φ f + φi is the volume fraction of total dispersed phase (i.e. filler particle and interfacial layer) following [28,53]

φgc = φ f (1 + l i ro )3

(8)

where l i is the interface thickness and ro is the particle radius while η and γ are given by [18,28]

η =  2 + (1 + l i ro )3  α fc − 2 1 − (1 + l i ro )3  α ic

(9)

Page 6 of 41

γ = 1 + 2(1 + l i ro )3  − 1 − (1 + l i ro )3  α fi

(10)

with α fc = Pf Pc , α ic = Pi Pc and α fi = Pf Pi . The Felske model simplifies to the Maxwell model when the thickness of the interfacial layer is negligible (l i → 0) [28,55]. Further, despite the permeability being a function of particle size in the Felske model, the model disregards inhomogeneities in the concentration field arising from nonlinear adsorption and finite particle size in the MMM, shown by us to significantly affect its permeability [8,14–17]. Consequently, the Felske model, in addition to being only applicable to low filler concentrations, is limited to MMMs having negligible filler particle sizes operating at low pressures. Only under these conditions, the concentration gradient across the MMM may be considered uniform[14,15]. The PTP Maxwell model in Table 1 assumes that the non-ideal MMM can be idealized as a pseudo-two-phase composite [23,57], in which the polymer matrix is one phase and the combined filler phase and interface is the other one (i.e. pseudo dispersed filler phase) [9,56]. To do so, the Maxwell model is first used to calculate the permeability in a composite comprising the filler particles and the interface via Eq. (3) [57], in which φ fi is the volume fraction of the filler phase in the filler-interface composite system following [23,25,26]

φ fi = φ f (φ f + φi ) = 1 (1 + l i ro )3

(11)

where ro and l i have the same connotation as in the Felske model. Then, the Maxwell model is used for a second time to calculate the permeability of the pseudo-two-phase system via Eq. (2), in which φ gc also has the same connotation as in the Felske model. Recently, this model has been shown consistent with the Laplace solution of the transport problem in the composite upon disregarding discontinuity-related effects at the interface in the Felske derivation [58]. As expected, the PTP Maxwell model leads to very similar predictions as the Felske model [8,9], with the model also simplifying to the Maxwell model when l i → 0 and thus is only applicable to diluted suspensions (φ gc ≤ 0.2) [23,25,26].

Page 7 of 41

While the above technique of using the Maxwell model multiple times to calculate the permeability in non-ideal MMMs was first introduced by Mahajan and Koros [25] to capture the effect of a voided interface, the approach has since been used in conjunction with other EMT models to describe common non-ideal filler-polymer morphologies [8,9,19,23,24,27,57], including the PTP Bruggeman [20] and Chiew-Glandt [53] models in Table 1. These adaptations of the Bruggeman [13] and Chiew-Glandt [59] models consider the MMM as a pseudo-two-phase composite by applying the EMT twice in the same fashion as the PTP Maxwell model [20,27], and thus φ gc and φ fi in Eqs. (4)-(7) have the same connotation as in the PTP Maxwell model. Here,

β gc = (α gc − 1) (α gc + 2) with α gc = Pg Pc and β fi = (α fi − 1) (α fi + 2) with α fi = Pf Pi in Eqs. (6) and (7), in which κ 2k is given by [60]

κ 2k = ak + bkφk

(12)

ak = −0.002254 − 0.123112 β k + 2.93656 β k2 + 1.6904 β k3

(13)

bk = 0.0039298 − 0.803494 β k − 2.16207 β k2 + 6.48296 β k3 + 5.27196 β k4

(14)

3 2

where k = gc in Eq. (6) and k = fi in Eq. (7). Further, while both Bruggeman and Chiew-Glandt models have been shown applicable to relatively concentrated systems (0 ≤ φ ≤ 0.4) [8,9,14– 16,27,60], their adaptations in Table 1 lack physical foundation in the multiple application of EMT. This is because despite EMT assumes that one of the filler phase is finely dispersed in the other one, so that the composite can be considered macroscopically homogeneous [14,61], such a condition is not met in Eqs. (5) and (7) because the filler-interface composite is not an actual dispersion. Commonly, permeation models in Table 1 are used in conjunction with experimental permeation data in the MMM to fit the interface region properties (i.e. thickness and permeability) [19,22– 27,53]. To do so, the interfacial permeability ( Pi ) is fixed to a given value while the interface thickness (l i ) is calculated upon error minimization between the experimental data and predictions based on a given permeation model, by assuming the filler particle size is known [8,27,57]. When polymer rigidification is considered at the filler particle surface, the interfacial permeability is often

Page 8 of 41

assumed a fraction of that in the bulk polymer following Pi = Pc γ , where γ , referred to as the chain immobilization factor, is an empirically fixed parameter used to differentiate the permeability of the rigid polymer layer from that of the bulk polymer matrix [23,27,53,54,57]. In practice, γ is often assumed between 3 − 4 [19,22,56,62], which leads to filler-polymer interface thicknesses of about 20 − 200 nm [22,23,27,57]. However, based on our atomistic simulations, we show here that this thickness is much lower than these values. 3. Simulation technique for gas permeation in non-ideal mixed-matrix membranes In this section, we first describe our approach to investigate gas adsorption and diffusion in the individual MMM phases and the polymer-filler interface region via nanoscale simulations. Subsequently, we describe our approach to characterize the gas transport in the MMM as a whole via macroscale simulations. 3.1. Equilibrium molecular dynamics simulations 3.1.1.

Model system and simulation details

A model MMM consisting of PI polymer and the inorganic filler is considered in our EMD simulations, as illustrated in Figure 1(a) for the MFI-PI composite. Further, PI polymer is used to represent the polymer matrix in both MFI-PI MMMs and CMS-Matrimid®5218. The neat PI polymer is constructed by considering 15 PI polymer chains, each having 10 monomers of biphenyltetracarboxylic dianhydride (BPDA) and 1,3-bis(4-aminophenoxy) benzene (APB). The initial polymer chains at very low density are generated using Packmol [63], and the final polymer structure is prepared by following a 21 -step compression/relaxation scheme, as detailed elsewhere [64]. The inter and intra-molecular non-bonded interactions in the polymer chains are described using the 9 − 6 form of Lennard-Jones (LJ) potential in accordance with the polymer consistent force field (PCFF) [65]; this established force field is extensively employed to represent the structure of PI polymer and PI-gas interactions [30,66,67]. In the MFI-PI system, all-silica type MFI zeolite, consisting of a 3-D pore network, is considered as filler. The potential parameters used to represent MFI are provided in the Supplementary data Page 9 of 41

(c.f. Table S-1) [68]; this force field is successfully employed to represent the gas adsorption and diffusion characteristics in several all-silica type zeolites including MFI [69]. The MFI surface is cleaved at the (1 0 0) plane, and all surface silica and oxygen atoms are saturated with − OH and −H groups, respectively. In the CMS-PI system, CMS is represented by a porous graphene having

nanopores of 5 − 6 Å diameter. The unsaturated carbon atoms on the rim of a nanopore are passivated with hydrogen as discussed in detail elsewhere [70,71]. In both system, the atoms on the surface layers are relaxed using a conjugate gradient method without optimizing the cell dimensions as well as bulk atoms using VASP software [72–74]. The resulting structures were subsequently used in the MD simulations without considering the flexibility of the framework.

Figure 1. (a) Structure of the PI-MFI hybrid system, and (b) Variation of PI polymer density with the distance from the surface in the presence of MFI and porous graphene. Here, ρ o is the bulk polymer density. In each composite, the thickness of the interface region is determined by computing the PI density near the inorganic surface, by dividing the simulation box into bins of 1 Å each in the direction normal to the inorganic surface ( z − direction in Figure 1(a)). Figure 1(b) depicts the relative PI density ( ρ ρo ) profiles for both composites as a function of position from the inorganic surface. Here, the polymer in the interface region is up to 30% and 130% denser than the bulk polymer in the presence of MFI-zeolite and porous graphene, respectively. This is associated with the relatively strong interaction between PI and porous graphene compared to that between PI and MFI-zeolite, and suggests that despite the interface thickness being nearly the same in both systems, the resistance offered to gas transport in the interface region is significantly larger in the CMS-PI than Page 10 of 41

in the MFI-zeolite/PI composite. Further, ρ ρ o steeply increases in the region 3 Å ≤ z ≤ 6 Å for both systems, which indicates polymer rigidification near both surfaces. Here, the first ρ ρ o peak in the presence of porous graphene is closer to the surface than that of PI near MFI-zeolite, which suggests the absence of cavities near the interface of the porous graphene, unlike the PI-MFI composite in which cavities arise due to the presence of end silanol groups near the interface. The length and proximity of first ρ ρ o peak to the inorganic surfaces suggest that the interface thickness, only consisting of a few polymer layers, is largely independent of the system, with this tendency differing from empirical results based on popular EMT models, in which the interface thickness has been shown strongly dependent on the filler-polymer pair combination [21–23,27,29]. Thus, the rigidified polymer region has a thickness of l i =1.2 nm ; however, small void like region at the external surface of MFI due surface terminal groups, having a thickness of 0.26 nm are considered to be part of the MFI-PI interface region. Therefore, the molecules located within this region of overall thickness of l i = 1.46 nm are then considered to compute gas adsorption and transport at the interface region for this system. In this study, CO 2 and CH 4 are chosen as adsorbates and represented by EPM2 [75] and 5 -site atomistic [76] models, respectively. The gas molecules are also treated as rigid in the entire simulation. A hybrid potential is used to model the non-bonded PI-MFI and PI-gas interactions, as detailed in the Supplementary data (c.f. Table S-2). For this we use modified interaction parameters for the polymer, corresponding to a 12 − 6 potential, by choosing the position at which the potential is zero to be the same as that of the 9 − 6 potential. This provides the 12 − 6 parameter σ = rm 1.5 3 1

, where rm is the position of the potential minimum of the 9 − 6 potential. Further, LorentzBerthelot rules are used to obtain the corresponding interaction parameters. Finally, the EMD simulations with periodic boundary conditions in all three dimensions are implemented using the LAMMPS package [77]. Temperature and pressure in the system are maintained by applying Nosé-Hoover thermostat and Berendsen barostat, respectively. Here, a cut-

Page 11 of 41

off distance of 1.4 nm is used for both potential energy and electrostatic interaction calculations. Further, long-range electrostatic interactions are incorporated using the Ewald summation method, and Verlet method with a time step of 1 fs used to integrate the particle equations of motion. Every simulation is run for 60 ns in the isothermal-isobaric (NPT) ensemble, with 10 ns allowed for equilibration. We note here that only polymer atoms are allowed to re-scale in the z-direction to control the pressure, while the position of filler atoms is unchanged. Then, corrected diffusivities are calculated by averaging over several independent atomistic simulation runs, each starting from a different initial configuration. In addition, the reported error bars are based on the standard deviation in the results determined by dividing the total simulation run into four equal parts. 3.1.2.

Adsorption, diffusion and permeation properties in the polymer, filler and interfacial region

First, The CO 2 and CH 4 adsorption properties in the neat polymer, filler as well as in the interface region are investigated by implementing GCMC simulations with insertion, deletion, rotational and translational moves having equal probability for approximately 2 ×107 steps, including the initial equilibration steps of 5 ×106 using DL_MONTE package [78]. The swelling behavior of polymer upon gas adsorption is accounted by implementing NPT-EMD simulations after GCMC simulations, as detailed elsewhere [79]. This procedure is repeated 10 times until a constant polymer density is achieved during the last 3 runs of the EMD simulation, indicating no further swelling of polymer upon gas adsorption. Here, the permeant capacity and associated error are computed by considering the last 3 GCMC runs. Then, the simulation results are fitted to Langmuir model in the filler and dual-mode model in the filler-polymer interface and bulk polymer, with [30,79] c f (Cb ) = k f c sf Cb (1 + k f Cb )

(15)

ci (Cb ) = k hi Cb + ki cis Cb (1 + ki Cb )

(16)

cc (Cb ) = khc Cb + kc ccs Cb (1 + kc Cb )

(17)

Page 12 of 41

where Cb = f g RT is the (ideal) gas pseudo-bulk concentration, with f g , T and R being the gas fugacity (or pressure), gas temperature and ideal gas constant, respectively, and k f , ki kc the permeant affinity constant in the filler phase, interface and bulk polymer phase, respectively, and c sf , cis and ccs are the gas saturation capacity in the filler, interface and bulk polymer, respectively.

Further, khi and khc are the permeant’s Henry law constants in the interface region and bulk polymer, respectively. Here, the application of the dual-mode model to adsorption in the polymer is empirical, as the effect of associated polymer swelling is embedded in the fitted isotherm parameters. Second, to describe the gas transport in the neat polymer, filler as well as in the interface region from the EMD simulations, corrected diffusivities ( Do ) of the gas molecule are determined. We note that for pure gases, Do is equivalent to Maxwell-Stefan (M-S) diffusivity [17,80]. In the neat filler and polymer materials, the Do is calculated by tracking the center of mass motion of all adsorbed molecules in the system using an Einstein relationship, as [81,82]

Do =

1 1 lim t →∞ 6N t

r r ∑ [ r (t ) − r (0)] N

k =1

k

2

k

(18)

with rk (t ) being the center of mass of molecule k at a time t . On the other hand, Do in the interface region of length l i , is determined by defining a collective coordinate n following [83]

dn = ∑ k ∈ interface(t ) [ dzk l i ]

(19)

where dzk is displacement of a gas molecule k in the z − direction during time dt in the interfacial region. The quantity n(t ) can be determined upon integrating Eq. (19) using the stored trajectory data from a given EMD simulation run. Since the interface region comprises an open system, the identities of the gas molecules in this region change with time and n(t ) executes a random walk, for which mean square displacement obeys the Einstein relation over sufficiently long time, leading to63

Dn = n2 (t ) 2t

(20)

where Dn , defined as the diffusion coefficient of n, can be related to Do , as [30]

Page 13 of 41

Do = Dn l i ( ρ Ac ) = Dn l 2i

where

N mol

(21)

N mol

is the average number of adsorbate molecules in the interface region, ρ is gas

density in the interface region and Ac is the area of the interface region. Third, we define the permeability in each phase of the MMMs by relating the transport (Fickian) diffusivity with Do through the Darken relation [84], using Eqs. (15)-(17) to describe the adsorption isotherms. Thus, the permeabilities in the filler, interface region and bulk polymer are given by: Cb 2 ≈0 1 Dof (Cb )  k f c sf (1 + k f Cb )  dCb ∫ RT (−∆Cb ) Cb1

(22)

Pi (Cb1 ) = −

Cb 2 ≈ 0 1 Dof (Cb )  khi + ki cis (1 + ki Cb )  dCb ∫ C RT (−∆Cb ) b1

(23)

Pc (Cb1 ) = −

Cb 2 ≈ 0 1 Doc (Cb )  khc + kc ccs (1 + kcCb )  dCb ∫ C RT (−∆Cb ) b1

(24)

Pf (Cb1 ) = −

where ∆Cb = Cb 2 − Cb1 and Dof , Doi and Doc are the corrected diffusivities in the filler, interface region, and polymer, respectively. Further, isotherm parameters (e.g. c sf , cis , ccs , k f , ki , kc khi , and khc ) have the same connotation as in Eq. (15)-(17). A comparison of the EMD simulation-based permeabilities in the MMM individual phases is presented in Section 4.1 (c.f. Figure 7). Isotherm parameters and corrected diffusivities are reported in the Supplementary data (c.f. Table S-3 and Table S-4) for both composite systems. For MFI-PI system, all individual phase parameters are calculated from the EMD simulation adsorption isotherms and corrected diffusivities. For the CMS-Matrimid®5218, isotherm parameters correspond to those reported in the literature [27,62,85] in the neat filler and polymer phases, and corrected diffusivities are back-calculated via Eqs. (22)(24) using literature reported neat phase permeabilities [27,62,85] because the gas diffusion coefficients and adsorption isotherms in the interfacial region are not accessible through EMD simulations due to strong PI polymer densification near the porous graphene (c.f. Figure 1(b)). 3.2. Full-scale simulation of the transport through MMMs

Page 14 of 41

We implemented the FEM to solve the coupled 3-D PDEs describing the transport in both ideal and non-ideal MMMs using the COMSOL Multiphysics® software package with LiveLinkTM for MATLAB®. In what follows, we limit ourselves to the description of our macroscale simulation procedure to calculate the permeability in non-ideal MMMs, as that for the ideal is detailed elsewhere [14,15]. First, 3-D MMMs are constructed with randomly distributed non-overlapping spherical fillers using a Force-Biased algorithm (FBA) [86] considering three types of particle size distributions. For every packing fraction and particle size distribution, we generated 5 independent configurations and averaged all simulation results amongst configurations. Figure 2(a) depicts representative simulation boxes with ro l = 0.020 , φ = 0.4 , in which ro

is the mean particle radius and l

the MMM thickness in the x -direction, while Figure 2(b) depicts corresponding particle size probability histograms. Here, the polydisperse distribution corresponds to a normal probability distribution with

ro

l = 0.020 and standard deviation of σ = 0.0040 . In all simulations, the

thickness (l i ) of the interfacial layer surrounding each particle (depicted in red color in Figure 2(a)) is considered uniform and independent of particle size, as discussed in Section 3.1.1.

Figure 2. Representative filler particle packing distributions with ro l = 0.020 and φ = 0.4 : (a) 3-D view of sphere packing structures, and (b) particle size probability histograms. Our full-scale simulations consider steady-state transport through 3-D MMMs following:

(∇ ⋅ J ) = 0

(25)

Page 15 of 41

with J being the steady-state flux at location ( x, y, z ) in the MMM, in which the Fick’s law is used to describe the flux in the filler ( J f ) , rigidified polymer ( J i ) , and bulk polymer ( J c ) phases, as [14–16]

where Cb = f g

J f = D f (Cb )( −∇Cb )

(26)

J i = Di (Cb )( −∇Cb )

(27)

J c = Dc (Cb )(−∇Cb ) (28) RT is the pseudo-bulk concentration at location ( x, y, z ) in the MMM, with f g

being the gas fugacity (or pressure) at location ( x, y, z ) , and T and R having the same connotation as defined in Section 3.1.2. Here, D f (Cb ) , Di (Cb ) and Dc (Cb ) are the local transport (Fickian) diffusivities in each phase, which are concentration-dependent through the Darken relation following:

D f (Cb ) = Dof (Cb )  k f c sf (1 + k f Cb ) 

(29)

Di (Cb ) = Doi (Cb )  khi + ki cis (1 + ki Cb ) 

(30)

Dc (Cb ) = Doc (Cb ) khc + kc ccs (1 + kcCb ) 

(31)

where all parameters in Eq. (29)-(31) have the same connotation as defined in Section 3.1.2. Further, we define the following boundary conditions: x = 0, Cb (0, y, z ) = Cb1

(32)

x = l, Cb (l, y, z ) = Cb 2

(33)

to solve Eqs. (25)-(31). Here, Cb1 and Cb 2 are the retentate and permeate pseudo-bulk concentrations of species i in the MMM, respectively. Further, we consider PBC at the membrane ends in the y and z directions and continuity of flux and pseudo-bulk concentrations across fillerpolymer boundaries (interface region). Figure 3 depicts the boundary conditions in a representative 2-D view of a single core-shell particle embedded in the polymer matrix. Further, we use adjustable tetrahedral meshes and a stationary fully coupled parallel sparse direct solver (MUMPS) to obtain all numerical solutions, as described elsewhere [14].

Page 16 of 41

Figure 3. Boundary conditions for the solution of Eqs. (25)-(31). After simultaneously solving Eqs. (25)-(31), we calculate the perm-selectivity (α P ) , as:

α P = PA PB

(34)

with the effective permeability of each permeant in the MMM estimated as follows [8] P = J x l [ RT ( − ∆Cb ) ] , ∆Cbi = Cb 2 − Cb1

(35)

where J x is the steady-state flux through the membrane, based on the mean value of J ( x, y, z ) at any location x , following [15,16]

ω ( x) = ∫∫ ω ( x, y, z )dydz yz

∫∫

yz

(36)

dydz

where ω ( x, y, z ) is the local value of a given property at location ( x, y, z ) and ω ( x) corresponds to the mean value of ω ( x, y, z ) over the yz plane at location x in the MMM. We note that while our macroscale simulations are implemented using spherical particles, this technique can accommodate other particle shapes (e.g. disks, cubes, cylinders). To do so, Eqs. (25)-(31) can be solved in the 3-D MMM containing such non-spherical fillers using a suitable hard particle packing simulation algorithm [87,88], in place of the FBA [86] used here to pack spherical fillers in the MMM. 4. Results and discussion We adopt the following conventions in the subsequent sections. We refer to CO 2 as species A and CH 4 as species B . In all figures, percentage deviation amongst theoretical profiles and experimental or simulation-based profiles accompany every theoretical profile, with this percentage deviation following color conventions of the figure legend. Percentage deviations are calculated as the percentage normalized root-mean-squared error, using the mean value of the simulation or experimental data to normalize all root-mean-squared errors.

Page 17 of 41

4.1. Adsorption, diffusion, and permeation in MFI-PI mixed matrix membranes As discussed in Section 3.1.2, the adsorption isotherms of pure CO 2 and CH 4 in each phase of the MFI-PI MMM are determined through GCMC simulations coupled with EMD simulations. Figure 4(a) depicts a comparison of the CO 2 capacity (c A ) in the polymer, filler and interface region in the pressure range of 0 − 30 atm at T = 300K , in which isotherms based on Eqs. (15)-(17) are shown (including the percentage deviation) using lines. Further, a similar plot for adsorption of CH 4 (cB ) is depicted in Figure 4(b). Here, it is seen that adsorption capacities of both gases in the MFI-zeolite at any given pressure are much higher than those in the neat PI polymer due to the availability of large free volume in MFI-zeolite. Further, a decrease in gas adsorption capacity in the rigidified polymer region is expected due to reduction in the free volume in this region compared to that in the bulk polymer region [30]. On the other hand, an overall enhancement in the gas adsorption capacity in the interface region is observed in Figure 4. This can be attributed to the cavities present near the PI-MFI surface arising from the presence of terminal functional groups on the MFI surface that act as strong sorption sites. A similar increase in gas adsorption capacity in the presence of grain boundaries has also been reported in ZIF-8 crystals [89]. Further, we note that gas adsorption isotherms in PI polymer and MFI filler are in good agreement with experimental data from the literature [90–92] , as shown in the Supplementary data (c.f. Figure S-1).

Figure 4. Adsorption isotherms of (a) CO 2 , and (b) CH 4 in individual MMM phases at T = 300 K . Lines correspond to fits of Eq. (15) in the filler phase, Eq. (16) in the interface and Eq. (16) in bulk polymer phase of the EMD simulation data.

Page 18 of 41

To investigate this increase in gas adsorption capacity in the interface region, the gas concentration is extracted as a function of position from the MFI surface (c.f. Figure 1(a)). Figure 5 depicts the position dependence of CO 2 concentration in the interface region. Here, a sharp increase in gas adsorption capacity compared to the bulk polymer is evident near the MFI surface (  3 Å from the surface), leading to overall increase in gas adsorption capacity in the interface region compared to the bulk polymer. However, it is seen that gas adsorption capacity in the core region of rigidified region is much lower compared to the bulk polymer. Further, we note that similar behavior is observed for methane adsorption capacity at the polymer-filler interface.

Figure 5. Variation of CO 2 adsorption capacity in the interface region with position, measured from the MFI surface. Dashed profile accompanying symbols are drawn as a guide for the eye. Subsequently, the pressure-dependence of the corrected diffusivities of CO 2 ( DoA ) and CH 4 ( DoB ) in all phases at T = 300 K were extracted and are depicted in Figure 6. It is seen that DoA and DoB are higher at any given pressure in the MFI-zeolite than those in the PI polymer. On the other hand, DoA and DoB are about 60% lower in the interface region compared to bulk PI polymer. This can be attributed to decrease in free volume and absence of large pores due to rigidification of the polymer near the zeolite surface [30].

Page 19 of 41

Figure 6. Corrected diffusivities of (a) CO 2 , and (b) CH 4 in individual MMM phases at T = 300 K. Lines correspond interpolation/extrapolation curves of the EMD simulation data. We note that the lines in Figure 6 correspond to estimations based on interpolation of the data based on Piecewise Cubic Hermite Interpolating Polynomials in the pressure range of 1 − 10 atm and linear extrapolation functions in the pressure range below 1 atm . Further, our CO 2 and CH 4 self-, corrected and transport diffusivities in MFI zeolite, reported with increase of loading in our previous work [30], are also comparable to PFG NMR [93] and quasi-elastic neutron scattering [94,95] measurements in various zeolites. Thus, corrected diffusivities based on these experimental techniques provide a means of validating EMD simulations. Finally, the permeabilities of CO 2 and CH 4 ( PA and PB ) in all the individual phases of MMM were computed using Eqs. (22)-(24) and are depicted in Figure 7(a) and (b), respectively. It is seen that PA and PB in MFI-zeolite are significantly higher than in the neat PI membrane. This can be attributed to higher gas diffusion and adsorption in MFI- zeolite compared to that of neat PI polymer membrane. On the other hand, PA and PB are about 10% and 50% lower in the PI- MFI interface than in PI, respectively. This is due to significant decrease in the gas diffusivity in the interface region than increase in gas adsorption capacity.

Page 20 of 41

Figure 7. Permeabilities of (a) CO 2 , and (b) CH 4 in individual MMM phases at T = 300 K using Eqs. (22)-(24) to calculate the pressure-dependence of the permeability. 4.2. Permeation in MFI-PI mixed matrix membranes for negligibly small particle size Given that EMT models are more appropriate when particle size is negligibly small in the MMMs (ro l → 0) , Figure 8 depicts a comparison of the gas permeability based on the models reported in Table 1 with those of our macroscale simulations with ro l = 0.004 . Here, Figure 8(a) depicts the CO 2 permeability ( PA ) profiles with increase of the filler volume fraction ( φ ) at T = 300 K and ∆pi = 5 atm while Figure 8(b) depicts those of CH 4 ( PB ) . In each case, the inset depicts PA and PB

profiles in equivalent ideal MMMs (i.e. same filler and polymer materials but ideal interface), in which lines correspond to the Maxwell, Bruggeman, and Chiew-Glandt models. We note that phase permeabilities in the theoretical models in Figure 8 are based on Eqs. (22)-(24), which average the transport diffusivity dependence on the pseudo-bulk concentration (pressure) over the membrane [17].

Page 21 of 41

Figure 8. Comparison of the (a) CO 2 and (b) CH 4 permeability profiles based on simulation and EMT models with increase of the mean filler volume fraction in ideal and non-ideal MMMs with ro l = 0.004 at T = 300 K and ∆pi = 5 atm . In the insets of Figure 8, the Maxwell model under-predicts the simulation permeabilities at high loadings while that of Bruggeman suggests an opposite tendency. The discrepancy between these models and our simulation is associated with the way both models consider particle-particle interactions. The Maxwell model, on one hand, assumes particles are too diluted to interact, and the Bruggeman model, on the other, considers interaction only amongst newly added particles with those already in the composite (new particles are too diluted to interact). Further, close agreement between the Chiew-Glandt model and simulation is seen in the insets of Figure 8, which is associated with the way particle-particle interactions are captured through pair-correlation functions of hard-sphere simulations in this model [60]. In Figure 8, it is also seen that simulation-based PA and PB are about 8 − 10% lower in the nonideal MMMs than in the ideal MMMs, which suggests that both CO 2 and CH 4 diffusion through the MMM is strongly sensitive to the interfacial resistance. Consequently, we compare in Figure 9 the position-dependent transport diffusivity profiles ( DmA , DmB ) across simulated ideal and nonideal MMMs with ro l = 0.004 and φ = 0.45 . Here, DmA and DmB are calculated through Eq. (36) and their profiles are depicted in dimensionless form based on the gas corrected diffusivities in MFI-zeolite ( DofA DofB ) at 5 atm (c.f. Table S-3).

Figure 9. Local simulation-based Fickian diffusivity (permeability) profiles of (a) CO 2 , and (b)

Page 22 of 41

CH 4 across an ideal and a non-ideal MMM with ro l = 0.004 at T = 300 K and ∆pi = 5 atm . Dotted profiles accompanying symbols are drawn as a guide for the eye. In Figure 9, the decrease of DmA and DmB in the regions neighboring the membrane ends is associated with depletion of the MFI-zeolite concentration in these regions [15], due to the consideration of a finite filler particle size relative to the MMM thickness in the macroscale simulations (c.f. Section 3.2). As a result of the finite particle size, there can be no particle centers in the regions 0 ≤ x / l < ro / l and (1 − ro / l) < x / l ≤ 1 , and the filler volume fraction profile is nonuniform in the regions 0 ≤ x / l < 2ro / l and (1 − 2ro / l) < x / l ≤ 1 . Both DmA and DmB are weakly sensitive to the interfacial resistance in the end regions, as demonstrated by the overlap of DmA and DmB profiles of both ideal and non-ideal MMMs. Alternatively, a steep decrease of DmA and DmB

profiles is seen in the region 2 ro l ≤ x l ≤ (1 − 2 ro l) for the non-ideal MMM, which is associated with increase of the interfacial resistance due to increase of the interface region volume fraction with increase of the filler volume fraction in this inner region. Further, the curvature variation of DmA and DmB profiles in this inner region 2 ro l ≤ x l ≤ (1 − 2 ro l) is due to strong sensitivity of

the gas transport to the combined effect of nonlinear gas adsorption and diffusion in the MMM individual phases (c.f. Section 4.1). It is also seen in Figure 9 that the CH 4 transport diffusivity is more sensitive to the interfacial resistance than that of CO 2 , as decrease in magnitude of DmB at 2 ro l ≤ x l ≤ (1 − 2 ro l) for the non-ideal MMM is greater than that of DmA in the same region. This tendency is associated with an interfacial size-exclusion effect, as CH 4 diffusion through the interface region is further limited by its larger size (kinetic diameter) in comparison to that of CO 2 [30].Further, this interfacial sizeexclusion effect leads to overall improvement of the separation efficiency, which is shown in Figure 10 depicting the perm-selectivity (α P ) profiles against the CO 2 permeability ( PA ) for both ideal and non-ideal MMMs with ro l = 0.004 and 0 ≤ φ ≤ 0.45 . Here, α P increases about 13% in the

Page 23 of 41

range 0 ≤ φ ≤ 0.45 for the non-ideal MMMs and slightly decreases about 0.1% for the ideal MMMs.

Figure 10. Variation of simulation-based perm-selectivity profiles with increase of the mean filler volume fraction in ideal and non-ideal MMMs with ro l = 0.004 at T = 300 K and ∆pi = 5 atm . 4.3. Effect of particle size on gas permeation in MFI-PI mixed matrix membranes We have also investigated the effect of the resistance offered by the interface region on the gas permeability with increase of particle size in the MMM. Thus, Figure 11 depicts a comparison of the simulation-based CO 2 and CH 4 permeability profiles, considering 0.004 ≤ ro l ≤ 0.160 in both ideal (closed symbols) and non-ideal (open symbols) MMMs, with those based on the Felske model (lines) at φ = 0.25, 0.35, 0.45 . The other EMT models of Table 1 are not considered here because their deviations are significantly larger (20 − 50%) than those of the Felske model in Figure 11, except predictions based on the PTP Maxwell model, which closely match those of Felske. As discussed in Section 3.1.1, we consider the thickness of interface region independent of particle size in both theoretical and simulation calculations, as effects of filler-polymer interactions only extend 12 − 15 Å from the MFI-PI interface (c.f. Figure 1(b)), and our prior simulations have shown this thickness to be invariant with filler particle size even for nano-sized crystals [30].

Page 24 of 41

Figure 11. Comparison of variation of (a) CO 2 and (b) CH 4 permeability with particle size based on the Felske model (solid and dashed lines) and on simulation, in ideal (closed symbols) and nonideal (open symbols) MMM at T = 300 K and ∆pi = 5 atm . Dotted profiles accompanying symbols are drawn as a guide for the eye. In Figure 11, PA and PB decrease with increase of particle size in the ideal MMMs, with this tendency associated with decrease of filler specific surface with increase of particle size for the same filler volume fraction [14,15]. On the other hand, PA and PB in the non-ideal MMMs increase with particle size in the region 0.004 < ro l ≤ 0.040 and closely overlap those of ideal MMMs in the region 0.040 < ro l ≤ 0.160 , which suggests that the permeability is sensitive to the interfacial resistance when particle size is small in the MMM. Here, the Felske model under-predicts the simulation-based PA and PB for every φ , with percentage deviations lying between 5 − 20% and increasing with increase of φ . Further, the theoretical PA and PB profiles are near constant with increase of particle size, with the model unable to mimic the steep increase of the simulation-based results with increase in particle size at 0.004 < ro l ≤ 0.040 , indicating that it is unable to capture the effect of particle size together with that of the interfacial resistance on the gas permeability. While a negligibly small particle size has been shown to lead to the highest permeability at fixed volume fraction in several ideal MMM systems [14,15,17], with this tendency also shown in Figure 11 for ideal MFI-PI MMMs and in the Supplementary data for ideal CMS-Matrimid®5218 MMMs (c.f. Figure S-2), the non-ideal simulation-based PA and PB profiles suggest that an optimum particle size can exists in non-ideal MMMs, and this does not correspond to the smallest particle

Page 25 of 41

size in the system. Here, PA and PB profiles exhibit a maximum at ro l = 0.040 for φ ≥ 0.35 , which is associated with the combination of two opposing finite-size-related effects on the effective permeability. On one hand, the permeability decreases with increase of particle size at fixed filler volume fraction, due to decrease of the lower resistance filler phase near the MMM ends (c.f. Figure 9), and due to increase in intra-particle resistance, with increase of particle size. On the other hand, the relative interfacial resistance decreases with increase of particle size at fixed filler volume fraction because the ratio of the interface thickness to particle size (l i ro ) is smaller for a larger particle. Further, this sensitivity of PA and PB profiles to the interfacial resistance leads to an increase of the perm-selectivity (α P ) with decrease of filler particle size. This is illustrated in Figure 12, which depicts a comparison of the variation of α P with ro l for φ = 0.45 based on the Felske model with simulation for both ideal and non-ideal MMMs. We note that while we found ro l = 0.040 to lead to a maximum in the gas permeability (flux) for this system (c.f. Figure 11), the optimum effective permeability in non-ideal MMMs is also dependent on the permeability in the interfacial region, Pi . This is because the interfacial resistance is proportional to l i ro and inversely proportional to Pi , and thus if Pi is decreased at fixed l i , the interfacial resistance increases; this means that only a larger particle size in the system can compensate for decrease of the MMM permeability due to increase of the resistance offered by the interface region. Further, this effect of a significantly low Pi (i.e. about 0.1 − 0.3% of Pc ) on the gas permeability is presented in the Supplementary data for the CMS-Matrimid®5218 MMMs. In this system, both gas permeability (c.f. Figure S-2) and perm-selectivity (c.f. Figure S-3) increase monotonically with particle size for non-ideal MMMs, and the largest possible particle size is preferred. Here, the Felske model provides better predictions of the CO 2 permeabilities than in Figure 11, while those of CH 4 have similar percentage deviations as in Figure 11. This tendency is

Page 26 of 41

expected, and is associated with the CO 2 having essentially linear adsorption isotherm in the polymer for this system in the pressure range of 0 − 50 psia .

Figure 12. Comparison of variation CO 2 CH 4 perm-selectivity with particle size based on the Felske model (solid line) and on simulation, in ideal (closed symbols) and non-ideal (open symbols) MMM with φ = 0.45 at T = 300 K and ∆pi = 5 atm . Dotted profiles accompanying symbols are drawn as a guide for the eye. Besides the above, a key outcome of our macroscale simulation technique is the consideration of the non-uniformity of the pseudo-bulk concentration (pressure) in the calculation of the gas permeability. This is a feature disregarded in classical EMT as well as its adaptations presented in Table 1, and which is strongly sensitive to both interfacial resistance and finite filler particle size, as illustrated in Figure 13 depicting a comparison of the CO 2 pseudo-bulk concentration profiles (Cb Cb1 ) of a yz -slice of a portion of (a) an ideal, and (b) non-ideal MMM. Here, the radius of the large particle size is 4 times that of the small one, and the interface region has the same thickness in both particles. In Figure 13(a), the interface region is only depicted as a reference, as transport properties in this region correspond to those in bulk PI. Arrows and streamlines (solid lines) in Figure 13 indicate the flux direction, with arrow thickness and number of streamlines being proportional to the flux magnitude. In Figure 13, it is seen that the gas transport through the large particle is only weakly sensitive to the resistance offered by the interface region. Here, the Cb Cb1 flow arrows and streamlines are nearly the same in the large particle in both Figure 13(a) and (b). This tendency agrees with the permeability profiles in Figure 11 closely overlapping with increase of particle size ( Page 27 of 41

0.040 < ro l ≤ 0.160 ). On the other hand, the permeant transport through the small particle is visibly sensitive to the interfacial resistance, as demonstrated by the Cb Cb1 color difference inside the small particle, arrows thicknesses, and number of streamlines when comparing Figure 13(a) and (b). This tendency is also in agreement with the permeability profiles in Figure 11 and confirms that the interfacial resistance has a significant effect on the gas transport when particle size is small due to increase of the interface region volume fraction in the MMM. Similar tendencies have been found recently for CO2 transport in ZIF-8/6FDA-DAM MMM through IR microimaging [96], in which CO2 was found to accumulate at the filler-polymer interfacial voided region.

Figure 13. Comparison of local dimensionless pseudo-bulk concentration (Cb Cb1 ) profiles of a yz -slice of a portion of a MMM with (a) ideal and (b) non-ideal filler-polymer interface. 4.4. Permeation in CMS-Matrimid®5218 mixed matrix membranes In this section, the permeation of CO2 and CH 4 in CMS-Matrimid®5218 MMMs is investigated. As discussed in Section 3.1.1, the thickness of the interface region is l i = 1.2 nm in this system. To estimate CO2 and CH 4 interfacial permeabilities, our macroscale simulation results were fitted to those reported in the literature [27,62,85] upon error minimization of the permeability predictions considering different chain immobilization factors (γ ) . To do so, we define the transport diffusivity in the interface region as Di = Decff γ , in which Deffc = Pco RT corresponds to the experimental gas permeabilities in the polymer (c.f. Table S-4). Figure 14 depicts a comparison of the experimental CO2 permeability ( PA ) and CH 4 ( PB ) with those based on simulation with increase of the mean volume fraction for three different chain immobilization factors. Here, the simulation profiles in an

Page 28 of 41

equivalent ideal MMM are also shown as a reference. In all simulations, we consider ro l = 0.020 , as experimental MMMs were reported to contain filler of mean particle size ro ≈ 0.5 − 1 µm , and had thickness 30 µm ≤ l ≤ 60 µm (i.e. ro l ≈ 0.0083 − 0.033 ).

Figure 14. Comparison of (a) CO 2 , and (b) CH 4 experimental permeabilities for the PI-CMS based MMM [27,62] with those based on simulation, for ro l = 0.020 , considering different chain immobilization factors (γ ) at T = 35 °C and ∆pi = 50 psia . Dotted profiles accompanying experimental data are drawn as a guide for the eye. In Figure 14, the simulation-based profile of γ A = 150 fits the experimental CO 2 permeabilities while that of γ B = 400 better fits those of CH 4 , with percentage deviations of 4.4% and 6.8% , respectively; this is consistent with the steep increase of the PI density for PI near the porous graphene (c.f. Figure 1(b)). Further, despite these percentage deviations being small, it is seen a better fit of the CO 2 experimental permeabilities than those of CH 4 and related to scatter of experimental result due to the smaller permeability values of CH 4 in comparison to those of CO 2 . Further, PA is more sensitive to the interfacial resistance than PB , which leads to decrease of about 10 − 40% of the perm-selectivity, also shown in the Supplementary data with increase of filler

particle sizes at φ = 0.40 (c.f. Figure S-3) and indicating that strong polymer densification at the particle surface has a significant effect on the overall permeant transport. As commonly done in practice [22,23,27,54], we back-calculated the thickness of the interface with the models in Table 1 considering γ A = 150 and γ B = 400 to estimate CO 2 and CH 4 interfacial permeabilities, respectively. Table 2 summarizes values of the interfacial thickness (l i ) Page 29 of 41

based on each one of these models along with the percentage deviation from the experimental CO 2 and CH 4 permeabilities. Here, the last two columns of Table 2 list the mean estimated thickness and its deviation from that predicted by the EMD simulations (l i = 1.2 nm) . Further, we consider the experimental permeabilities in the range 0 ≤ φ ≤ 0.25 to estimate the interfacial thickness in all cases, as EMT models in Table 1 were shown to only fit well the simulation-based permeabilities at low loadings (0 ≤ φ ≤ 0.25) in Figure 8. Table 2. Prediction of the interface thickness based on EMT models using ro l = 0.020 . CO 2 Model Felske model [18] PTP Maxwell model [25] PTP Bruggeman model [20] PTP Chiew-Glandt [53]

l i [nm] 1.95 1.95 23.55 9.88

CH 4 PA Error [%] 8.01 8.01 8.01 8.14

l i [nm] 3.17 3.17 20.81 13.67

PB Error [%] 10.90 10.90 10.99 11.03

l i ± σ l i [nm] 2.56 ± 0.86 2.56 ± 0.86 22.18 ± 1.94 11.78 ± 2.68

li

Error [%] 113 113 1748 882

In Table 2, despite percentage deviations between the experimental permeabilities and prediction of the EMT being relatively small ( 8 − 11% ), these models largely over-estimate the interface thickness when permeation data in the MMM are fitted, with the Felske and PTP Maxwell models having the smallest deviation of about 100%. This tendency is associated with the inability of these models to capture finite-size effects, which are masked by the fitting of the interface thickness, evident in the comparison presented in the Supplementary data (c.f. Figure S-2) of the Felske model to our macroscale simulations considering l i = 1.2 nm with increase of particle size. Further, while the Bruggeman and Chiew-Glandt models have been shown to closely match the permeability in various near ideal MMMs [9,14,27,60], including the ideal MFI-zeolite/PI MMMs investigated here (c.f. Figure 8), their adaptions largely over-predict the interface thickness in Table 2 (deviations between 800 − 1700% ). Thus, the discrepancy between theoretical predictions in Table 2 suggests that it is not appropriate to apply the effective medium approximation to the particle-interface composite by considering it as a random dispersion (c.f. Section 2), and consequently further

Page 30 of 41

theoretical developments are needed to accurately capture the effect of the interfacial resistance together with particle size through the EMT. 4.5. Effect of particle size distribution on gas permeation in CMS-Matrimid®5218 mixed matrix membranes Because the pseudo-bulk concentration profiles in Figure 13 were shown strongly sensitive to the interfacial resistance with decrease of particle size, we also study the effect of particle size distribution on the gas permeability and perm-selectivity in CMS-Matrimid®5218 MMMs. Figure 15 depicts a comparison of the experimental CO 2 permeability ( PA ) and CO 2 CH 4 permselectivity (α P ) with those based on simulation for three different particle size distributions with mean particle size of ro l = 0.020 (c.f. Figure 2). For both gases, closed symbols correspond to ideal MMMs while open symbols to those of the non-ideal MMMs. In Figure 15, both PA and α P profiles of ideal MMMs overlap for all distributions, with simulation profiles having similar deviations from the experimental data for all distributions of about 31% and 20% , respectively. This tendency indicates that permeability and perm-selectivity in ideal MMMs are only weakly sensitive to the effects of non-uniform particle size distribution. On the other hand, while both PA and α P profiles of the non-ideal MMMs with monodisperse and polydisperse distributions overlap in Figure 15, those of the bidisperse distribution fall above them. This tendency suggests that the permeability and perm-selectivity in non-ideal MMM are sensitive to non-uniform particle size distribution, especially when there is a large difference between the particle sizes in the system (c.f. Figure 2), and is associated with decrease of the interfacial resistance with increase of the number of large particles in the system (c.f. Figure 13). While this effect is also present in any polydisperse distribution, the probability of having large particles in the polydisperse distribution of Figure 15 is less than 1% in comparison to that of the bidisperse distribution of about 30% (c.f. Figure 2). This tendency suggests that wide particle size distributions lead to greater deviations of the permeability from the predictions based on the mean particle size in

Page 31 of 41

non-ideal MMMs; in agreement with the CO 2 and CH 4 permeability profiles presented in the Supplementary data (c.f. Figure S-2) for this system, exhibiting a permeability increase with increase of particle size.

Figure 15. Comparison of simulation-based and experimental (a) CO 2 permeability, and (b) CO 2 CH 4 perm-selectivity for the PI-CMS based MMM[27,62], for various particle size distributions, with ro l = 0.020 at T = 35 °C and ∆pi = 50 psia . Dotted profiles accompanying experimental data are drawn as a guide for the eye. 5. Conclusions We introduce a multiscale simulation technique to predict the effective gas permeability in a nonideal Mixed matrix membrane (MMM). Our approach integrates equilibrium molecular dynamics simulation predictions of the adsorption isotherms and diffusion coefficients with full-scale simulations solving the coupled 3-D partial differential equations at the macroscale via the finiteelement method. Based on our results, existing adaptations of classical EMT models are found to deviate significantly from those of simulation, as they fail to capture effects of particle-particle interactions together with the interfacial resistance on the gas permeability, and lead to very large fitted values of the filler-polymer interfacial thickness (with a maximum deviation of 1800% ). Furthermore, our macroscale simulations reveal the existence of an optimum particle size in nonideal MMMs, that leads to optimum permeabilities at any filler volume fraction, which arises due to the competing effects of decrease in the MMM permeability with increase of particle size and decrease of the interfacial resistance with increase of particle size.

Page 32 of 41

By using our multiscale simulation technique, we have investigated the transport of CO 2 and CH 4 in MFI-zeolite/PI and CMS-PI MMMs. In both systems, the interfacial permeability is found to be lower than that in the bulk polymer, with this tendency associated with reduced gas diffusion in the interface region due increase of the polymer density near the filler surface. Further, we find that although the thickness of the interfacial layer is only about 1.2 − 1.5 nm , the MMM permselectivity is strongly sensitive to the interfacial resistance. In both PI-based MMMs, the effective permeability in non-ideal MMMs is much lower than that of the ideal MMM when the particle size is small with such sensitivity to the interfacial resistance disappearing with increase of particle size. Finally, we have also studied the effect of the particle size distribution on the gas permeability in both ideal and non-ideal MMMs, never considered before in existing fundamental approaches. We found that the permeability is weakly sensitive to non-uniform particle size distributions in ideal MMMs, for which case the MMM permeability can be based in the mean particle size of the system. On the other hand, we found that the permeability is sensitive to non-uniform particle size distributions in non-ideal MMMs, and deviates from that based on the mean particle size when the MMM has a wide particle size distribution. Here, the MMMs having bidisperse distribution improved gas diffusion due to significant decrease of the interfacial resistance in the large particles. Acknowledgments This research has been supported by a grant (No. DP150101824) from the Australian Research Council, through the Discovery scheme. This research was undertaken with the assistance of the computational resources provided at the NCI National Facility systems at the Australian National University (ANU), and those at the Pawsey Supercomputing Centre in Western Australia, through their National Computational Merit Allocation Schemes supported by the Australian Government and the Government of Western Australia. Appendix A. Supplementary data Supplementary data to this article can be found online.

Page 33 of 41

Page 34 of 41

References [1]

M. Galizia, W.S. Chi, Z.P. Smith, T.C. Merkel, R.W. Baker, B.D. Freeman, 50th Anniversary Perspective: Polymers and Mixed Matrix Membranes for Gas and Vapor Separation: A Review and Prospective Opportunities, Macromolecules. (2017). doi:10.1021/acs.macromol.7b01718.

[2]

Y. Cheng, Z. Wang, D. Zhao, Mixed Matrix Membranes for Natural Gas Upgrading: Current Status and Opportunities, Ind. Eng. Chem. Res. 57 (2018) 4139–4169. doi:10.1021/acs.iecr.7b04796.

[3]

H.B. Park, J. Kamcev, L.M. Robeson, M. Elimelech, B.D. Freeman, Maximizing the right stuff: The trade-off between membrane permeability and selectivity, Science (80-. ). 356 (2017). doi:10.1126/science.aab0530.

[4]

M. Vinoba, M. Bhagiyalakshmi, Y. Alqaheem, A.A. Alomair, A. Pérez, M.S. Rana, Recent progress of fillers in mixed matrix membranes for CO2 separation: A review, Sep. Purif. Technol. 188 (2017) 431–450. doi:10.1016/j.seppur.2017.07.051.

[5]

L.M. Robeson, The upper bound revisited, J. Memb. Sci. 320 (2008) 390–400. doi:10.1016/j.memsci.2008.04.030.

[6]

M. Rezakazemi, A. Ebadi Amooghin, M.M. Montazer-Rahmati, A.F. Ismail, T. Matsuura, State-of-the-art membrane based CO2 separation using mixed matrix membranes (MMMs): An overview on current status and future directions, Prog. Polym. Sci. 39 (2014) 817–861. doi:10.1016/j.progpolymsci.2014.01.003.

[7]

T.-S. Chung, L.Y. Jiang, Y. Li, S. Kulprathipanja, Mixed matrix membranes (MMMs) comprising organic polymers with dispersed inorganic fillers for gas separation, Prog. Polym. Sci. 32 (2007) 483–507. doi:10.1016/j.progpolymsci.2007.01.008.

[8]

G.M. Monsalve-Bravo, S.K. Bhatia, Modeling Permeation through Mixed-Matrix Membranes: A Review, Processes. 6 (2018) 172. doi:10.3390/pr6090172.

[9]

H. Vinh-Thang, S. Kaliaguine, Predictive Models for Mixed-Matrix Membrane Performance: A Review, Chem. Rev. 113 (2013) 4980–5028. doi:10.1021/cr3003888.

[10] B. Zornoza, C. Tellez, J. Coronas, J. Gascon, F. Kapteijn, Metal organic framework based mixed matrix membranes: An increasingly important field of research with a large application potential, Microporous Mesoporous Mater. 166 (2013) 67–78. doi:10.1016/j.micromeso.2012.03.012. [11] J.C. Maxwell, A treatise on electricity and magnetism, Clarendon Press, Oxford, UK, 1873. [12] Y.C. Chiew, E.D. Glandt, Effective conductivity of dispersions: The effect of resistance at the particle surfaces, Chem. Eng. Sci. 42 (1987) 2677–2685. doi:10.1016/00092509(87)87018-5. [13] D.A.G. Bruggeman, Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen, Ann. Phys. 24 (1935) 636. doi:10.1002/andp.19354160705. [14] G.M. Monsalve-Bravo, S.K. Bhatia, Extending Effective Medium Theory to Finite Size Systems: Theory and Simulation for Permeation in Mixed-Matrix Membranes, J. Memb. Sci. 531 (2017) 148–159. doi:10.1016/j.memsci.2017.02.029. [15] G.M. Monsalve-Bravo, S.K. Bhatia, Concentration-dependent transport in finite sized

Page 35 of 41

composites: Modified effective medium theory, J. Memb. Sci. 550 (2018) 110–125. doi:10.1016/j.memsci.2017.12.058. [16] G.M. Monsalve-Bravo, S.K. Bhatia, Comparison of Hollow Fiber and Flat Mixed-Matrix Membranes: Theory and Simulation, Chem. Eng. Sci. 187 (2018) 174–188. doi:10.1016/j.ces.2018.04.037. [17] G.M. Monsalve-Bravo, S. Smart, S.K. Bhatia, Simulation of Multicomponent Gas Transport through Mixed-Matrix Membranes, J. Memb. Sci. 577 (2019) 219–234. doi:10.1016/j.memsci.2019.02.013. [18] J.D. Felske, Effective thermal conductivity of composite spheres in a continuous medium with contact resistance, Int. J. Heat Mass Transf. 47 (2004) 3453–3461. doi:10.1016/j.ijheatmasstransfer.2004.01.013. [19] A. Idris, Z. Man, S.M. Abdulhalim, F. Uddin, Modified Bruggeman models for prediction of CO2 permeance in polycarbonate/silica nanocomposite membranes, Can. J. Chem. Eng. 95 (2017) 2398–2409. doi:10.1002/cjce.22933. [20] A. Shariati, M. Omidkhah, M.Z. Pedram, New permeation models for nanocomposite polymeric membranes filled with nonporous particles, Chem. Eng. Res. Des. 90 (2012) 563– 575. doi:10.1016/j.cherd.2011.08.010. [21] E. Chehrazi, A. Sharif, M. Omidkhah, M. Karimi, Modeling the Effects of Interfacial Characteristics on Gas Permeation Behavior of Nanotube–Mixed Matrix Membranes, ACS Appl. Mater. Interfaces. 9 (2017) 37321–37331. doi:10.1021/acsami.7b11545. [22] Y. Li, H.-M. Guan, T.-S. Chung, S. Kulprathipanja, Effects of novel silane modification of zeolite surface on polymer chain rigidification and partial pore blockage in polyethersulfone (PES)–zeolite A mixed matrix membranes, J. Memb. Sci. 275 (2006) 17–28. doi:10.1016/j.memsci.2005.08.015. [23] T.T. Moore, W.J. Koros, Non-ideal effects in organic–inorganic materials for gas separation membranes, J. Mol. Struct. 739 (2005) 87–98. doi:10.1016/j.molstruc.2004.05.043. [24] Z. Sadeghi, M. Omidkhah, M.E. Masoumi, R. Abedini, Modification of existing permeation models of mixed matrix membranes filled with porous particles for gas separation, Can. J. Chem. Eng. 94 (2016) 547–555. doi:10.1002/cjce.22414. [25] R. Mahajan, W.J. Koros, Mixed matrix membrane materials with glassy polymers. Part 1, Polym. Eng. Sci. 42 (2002) 1420–1431. doi:10.1002/pen.11041. [26] Y. Shen, A.C. Lua, Theoretical and experimental studies on the gas transport properties of mixed matrix membranes based on polyvinylidene fluoride, AIChE J. 59 (2013) 4715–4726. doi:10.1002/aic.14186. [27] D.Q. Vu, W.J. Koros, S.J. Miller, Mixed matrix membranes using carbon molecular sieves II. Modeling permeation behavior, J. Memb. Sci. 211 (2003) 335–348. doi:10.1016/S03767388(02)00425-8. [28] R. Pal, Permeation models for mixed matrix membranes, J. Colloid Interface Sci. 317 (2008) 191–198. doi:10.1016/j.jcis.2007.09.032. [29] F.P. Di Maio, A. Santaniello, A. Di Renzo, G. Golemme, Description of gas transport in perfluoropolymer/SAPO-34 mixed matrix membranes using four-resistance model, Sep. Purif. Technol. 185 (2017) 160–174. doi:10.1016/j.seppur.2017.05.024.

Page 36 of 41

[30] R.C. Dutta, S.K. Bhatia, Structure and Gas Transport at the Polymer–Zeolite Interface: Insights from Molecular Dynamics Simulations, ACS Appl. Mater. Interfaces. 10 (2018) 5992–6005. doi:10.1021/acsami.7b17470. [31] M. Benzaqui, R. Semino, N. Menguy, F. Carn, T. Kundu, J.-M. Guigner, N.B. McKeown, K.J. Msayib, M. Carta, R. Malpass-Evans, C. Le Guillouzer, G. Clet, N.A. Ramsahye, C. Serre, G. Maurin, N. Steunou, Toward an Understanding of the Microstructure and Interfacial Properties of PIMs/ZIF-8 Mixed Matrix Membranes, ACS Appl. Mater. Interfaces. 8 (2016) 27311–27321. doi:10.1021/acsami.6b08954. [32] R. Semino, N.A. Ramsahye, A. Ghoufi, G. Maurin, Microscopic Model of the Metal– Organic Framework/Polymer Interface: A First Step toward Understanding the Compatibility in Mixed Matrix Membranes, ACS Appl. Mater. Interfaces. 8 (2016) 809–819. doi:10.1021/acsami.5b10150. [33] A.-C. Yang, C.-H. Liu, D.-Y. Kang, Estimations of effective diffusivity of hollow fiber mixed matrix membranes, J. Memb. Sci. 495 (2015) 269–275. doi:10.1016/j.memsci.2015.08.030. [34] T. Singh, D.-Y. Kang, S. Nair, Rigorous calculations of permeation in mixed-matrix membranes: Evaluation of interfacial equilibrium effects and permeability-based models, J. Memb. Sci. 448 (2013) 160–169. doi:10.1016/j.memsci.2013.08.010. [35] T. Wang, D.-Y. Kang, Predictions of effective diffusivity of mixed matrix membranes with tubular fillers, J. Memb. Sci. 485 (2015) 123–131. doi:10.1016/j.memsci.2015.03.028. [36] T. Wang, D.-Y. Kang, Highly selective mixed-matrix membranes with layered fillers for molecular separation, J. Memb. Sci. 497 (2015) 394–401. doi:10.1016/j.memsci.2015.09.057. [37] M.M.M. Sharifzadeh, M.Z. Pedram, A.E. Amooghin, A new permeation model in porous filler–based mixed matrix membranes for CO2 separation, Greenh. Gases Sci. Technol. 9 (2019) 719–742. doi:10.1002/ghg.1891. [38] C. Zhang, Y. Dai, J.R. Johnson, O. Karvan, W.J. Koros, High performance ZIF-8/6FDADAM mixed matrix membrane for propylene/propane separations, J. Memb. Sci. 389 (2012) 34–42. doi:10.1016/j.memsci.2011.10.003. [39] C. Zhang, R.P. Lively, K. Zhang, J.R. Johnson, O. Karvan, W.J. Koros, Unexpected Molecular Sieving Properties of Zeolitic Imidazolate Framework-8, J. Phys. Chem. Lett. 3 (2012) 2130–2134. doi:10.1021/jz300855a. [40] Y. Pan, W. Liu, Y. Zhao, C. Wang, Z. Lai, Improved ZIF-8 membrane: Effect of activation procedure and determination of diffusivities of light hydrocarbons, J. Memb. Sci. 493 (2015) 88–96. doi:https://doi.org/10.1016/j.memsci.2015.06.019. [41] D. Liu, X. Ma, H. Xi, Y.S. Lin, Gas transport properties and propylene/propane separation characteristics of ZIF-8 membranes, J. Memb. Sci. 451 (2014) 85–93. doi:10.1016/j.memsci.2013.09.029. [42] Q. Song, S.K. Nataraj, M. V Roussenova, J.C. Tan, D.J. Hughes, W. Li, P. Bourgoin, M.A. Alam, A.K. Cheetham, S.A. Al-Muhtaseb, E. Sivaniah, Zeolitic imidazolate framework (ZIF8) based polymer nanocomposite membranes for gas separation, Energy Environ. Sci. 5 (2012) 8359–8369. doi:10.1039/C2EE21996D. [43] R. Mueller, V. Hariharan, C. Zhang, R. Lively, S. Vasenkov, Relationship between mixed Page 37 of 41

and pure gas self-diffusion for ethane and ethene in ZIF-8/6FDA-DAM mixed-matrix membrane by pulsed field gradient NMR, J. Memb. Sci. 499 (2016) 12–19. doi:10.1016/j.memsci.2015.10.036. [44] R. Mueller, S. Zhang, C. Zhang, R. Lively, S. Vasenkov, Relationship between long-range diffusion and diffusion in the ZIF-8 and polymer phases of a mixed-matrix membrane by high field NMR diffusometry, J. Memb. Sci. 477 (2015) 123–130. doi:10.1016/j.memsci.2014.12.015. [45] E.M. Forman, A. Baniani, L. Fan, K.J. Ziegler, E. Zhou, F. Zhang, R.P. Lively, S. Vasenkov, Ethylene diffusion in crystals of zeolitic imidazole Framework-11 embedded in polymers to form mixed-matrix membranes, Microporous Mesoporous Mater. 274 (2019) 163–170. doi:https://doi.org/10.1016/j.micromeso.2018.07.044. [46] I. Erucar, S. Keskin, Computational Methods for MOF/Polymer Membranes, Chem. Rec. 16 (2016) 703–718. doi:10.1002/tcr.201500275. [47] S. Keskin, D.S. Sholl, Selecting metal organic frameworks as enabling materials in mixed matrix membranes for high efficiency natural gas purification, Energy Environ. Sci. 3 (2010) 343–351. doi:10.1039/B923980B. [48] I. Erucar, S. Keskin, Computational screening of metal organic frameworks for mixed matrix membrane applications, J. Memb. Sci. 407–408 (2012) 221–230. doi:https://doi.org/10.1016/j.memsci.2012.03.050. [49] G. Yilmaz, S. Keskin, Molecular modeling of MOF and ZIF-filled MMMs for CO2/N2 separations, J. Memb. Sci. 454 (2014) 407–417. doi:https://doi.org/10.1016/j.memsci.2013.12.045. [50] I. Erucar, S. Keskin, Molecular Modeling of MOF-based Mixed Matrix Membranes, Curr. Org. Chem. 18 (2014) 2364–2380. https://www.ingentaconnect.com/content/ben/coc/2014/00000018/00000018/art00007. [51] C. Altintas, S. Keskin, Molecular Simulations of MOF Membranes and Performance Predictions of MOF/Polymer Mixed Matrix Membranes for CO2/CH4 Separations, ACS Sustain. Chem. Eng. 7 (2019) 2739–2750. doi:10.1021/acssuschemeng.8b05832. [52] M. Tatlier, Ç. Atalay-Oral, Estimation of interphase properties of various mixed matrix membranes, Compos. Interfaces. (2018) 1–13. doi:10.1080/09276440.2018.1535105. [53] S. Maghami, M. Sadeghi, A. Mehrabani-Zeinabad, Recognition of polymer-particle interfacial morphology in mixed matrix membranes through ideal permeation predictive models, Polym. Test. 63 (2017) 25–37. doi:10.1016/j.polymertesting.2017.07.021. [54] Y. Li, T.-S. Chung, C. Cao, S. Kulprathipanja, The effects of polymer chain rigidification, zeolite pore size and pore blockage on polyethersulfone (PES)-zeolite A mixed matrix membranes, J. Memb. Sci. 260 (2005) 45–55. doi:10.1016/j.memsci.2005.03.019. [55] A.J. Petsi, V.N. Burganos, Interphase layer effects on transport in mixed matrix membranes, J. Memb. Sci. 421 (2012) 247–257. doi:10.1016/j.memsci.2012.07.021. [56] B. Shimekit, H. Mukhtar, T. Murugesan, Prediction of the relative permeability of gases in mixed matrix membranes, J. Memb. Sci. 373 (2011) 152–159. doi:10.1016/j.memsci.2011.02.038. [57] T.T. Moore, R. Mahajan, D.Q. Vu, W.J. Koros, Hybrid membrane materials comprising

Page 38 of 41

organic polymers with rigid dispersed phases, AIChE J. 50 (2004) 311–321. doi:10.1002/aic.10029. [58] S. Maghami, M. Sadeghi, A. Mehrabani-Zeinabad, M. Zarabadi, B. Ghalei, The Role of Interfacial Morphology in the Gas Transport Behavior of Nanocomposite Membranes: A Mathematical Modeling Approach, Ind. Eng. Chem. Res. 58 (2019) 11022–11037. doi:10.1021/acs.iecr.9b01398. [59] Y.C. Chiew, E.D. Glandt, The effect of structure on the conductivity of a dispersion, J. Colloid Interface Sci. 94 (1983) 90–104. doi:10.1016/0021-9797(83)90238-2. [60] E. Gonzo, M. Parentis, J. Gottifredi, Estimating models for predicting effective permeability of mixed matrix membranes, J. Memb. Sci. 277 (2006) 46–54. doi:10.1016/j.memsci.2005.10.007. [61] R. Pal, New models for thermal conductivity of particulate composites, J. Reinf. Plast. Compos. 26 (2007) 643–651. doi:10.1177/0731684407075569. [62] D.Q. Vu, W.J. Koros, S.J. Miller, Mixed matrix membranes using carbon molecular sieves: I. Preparation and experimental results, J. Memb. Sci. 211 (2003) 311–334. doi:10.1016/S0376-7388(02)00429-5. [63] L. Martínez, R. Andrade, E.G. Birgin, J.M. Martínez, PACKMOL: A package for building initial configurations for molecular dynamics simulations, J. Comput. Chem. 30 (2009) 2157–2164. doi:10.1002/jcc.21224. [64] G.S. Larsen, P. Lin, K.E. Hart, C.M. Colina, Molecular Simulations of PIM-1-like Polymers of Intrinsic Microporosity, Macromolecules. 44 (2011) 6944–6951. doi:10.1021/ma200345v. [65] H. Sun, S.J. Mumby, J.R. Maple, A.T. Hagler, An ab Initio CFF93 All-Atom Force Field for Polycarbonates, J. Am. Chem. Soc. 116 (1994) 2978–2987. doi:10.1021/ja00086a030. [66] R.C. Dutta, S.K. Bhatia, Atomistic Investigation of Mixed-Gas Separation in a Fluorinated Polyimide Membrane, ACS Appl. Polym. Mater. 1 (2019) 1359–1371. doi:10.1021/acsapm.9b00146. [67] R.C. Dutta, S.K. Bhatia, Interfacial barriers to gas transport in zeolites: distinguishing internal and external resistances, Phys. Chem. Chem. Phys. 20 (2018) 26386–26395. doi:10.1039/C8CP05834B. [68] A. Goj, D.S. Sholl, E.D. Akten, D. Kohen, Atomistic Simulations of CO2 and N2 Adsorption in Silica Zeolites: The Impact of Pore Size and Shape, J. Phys. Chem. B. 106 (2002) 8367– 8375. doi:10.1021/jp025895b. [69] K. Makrodimitris, G.K. Papadopoulos, D.N. Theodorou, Prediction of Permeation Properties of CO2 and N2 through Silicalite via Molecular Simulations, J. Phys. Chem. B. 105 (2001) 777–788. doi:10.1021/jp002866x. [70] R.C. Dutta, S. Khan, J.K. Singh, Wetting transition of water on graphite and boron-nitride surfaces: A molecular dynamics study, Fluid Phase Equilib. 302 (2011) 310–315. doi:10.1016/j.fluid.2010.07.006. [71] M. Lalitha, S. Lakshmipathi, S.K. Bhatia, Defect-Mediated Reduction in Barrier for Helium Tunneling through Functionalized Graphene Nanopores, J. Phys. Chem. C. 119 (2015) 20940–20948. doi:10.1021/acs.jpcc.5b05567.

Page 39 of 41

[72] G. Kresse, J. Furthmüller, Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev. B. 54 (1996) 11169–11186. doi:10.1103/PhysRevB.54.11169. [73] G. Kresse, J. Furthmüller, Efficiency of ab-initio total energy calculations for metals and semiconductors using a plane-wave basis set, Comput. Mater. Sci. 6 (1996) 15–50. doi:10.1016/0927-0256(96)00008-0. [74] G. Kresse, J. Hafner, Ab initio molecular dynamics for liquid metals, Phys. Rev. B. 47 (1993) 558–561. doi:10.1103/PhysRevB.47.558. [75] J.G. Harris, K.H. Yung, Carbon Dioxide’s Liquid-Vapor Coexistence Curve And Critical Properties as Predicted by a Simple Molecular Model, J. Phys. Chem. 99 (1995) 12021– 12024. doi:10.1021/j100031a034. [76] Y. Sun, D. Spellmeyer, D.A. Pearlman, P. Kollman, Simulation of the solvation free energies for methane, ethane, and propane and corresponding amino acid dipeptides: a critical test of the bond-PMF correction, a new set of hydrocarbon parameters, and the gas phase-water hydrophobicity scale, J. Am. Chem. Soc. 114 (1992) 6798–6801. doi:10.1021/ja00043a027. [77] S. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J. Comput. Phys. 117 (1995) 1–19. doi:10.1006/jcph.1995.1039. [78] J.A. Purton, J.C. Crabtree, S.C. Parker, DL_MONTE: a general purpose program for parallel Monte Carlo simulation, Mol. Simul. 39 (2013) 1240–1252. doi:10.1080/08927022.2013.839871. [79] R.C. Dutta, S.K. Bhatia, Transport Diffusion of Light Gases in Polyethylene Using Atomistic Simulations, Langmuir. 33 (2017) 936–946. doi:10.1021/acs.langmuir.6b04037. [80] A.I. Skoulidas, D.S. Sholl, Self-Diffusion and Transport Diffusion of Light Gases in MetalOrganic Framework Materials Assessed Using Molecular Dynamics Simulations, J. Phys. Chem. B. 109 (2005) 15760–15768. doi:10.1021/jp051771y. [81] S.K. Bhatia, M.R. Bonilla, D. Nicholson, Molecular transport in nanopores: a theoretical perspective, Phys. Chem. Chem. Phys. 13 (2011) 15350–15383. doi:10.1039/C1CP21166H. [82] D.S. Sholl, Understanding Macroscopic Diffusion of Adsorbed Molecules in Crystalline Nanoporous Materials via Atomistic Simulations, Acc. Chem. Res. 39 (2006) 403–411. doi:10.1021/ar0402199. [83] F. Zhu, E. Tajkhorshid, K. Schulten, Collective Diffusion Model for Water Permeation through Microscopic Channels, Phys. Rev. Lett. 93 (2004) 224501. doi:10.1103/PhysRevLett.93.224501. [84] L.S. Darken, Diffusion, mobility and their interrelation through free energy in binary metallic systems, Trans. Am. Inst. Min. Engrs. 175 (1948) 184–194. [85] D.Q. Vu, Formation and characterization of asymmetric carbon molecular sieve and mixedmatrix membranes for natural gas purification, Ph.D. Dissertation, The University of Texas at Austin, TX, 2001. [86] V. Baranau, U. Tallarek, Random-close packing limits for monodisperse and polydisperse hard spheres, Soft Matter. 10 (2014) 3826–3841. doi:10.1039/C3SM52959B. [87] G.W. Delaney, P.W. Cleary, The packing properties of superellipsoids, EPL (Europhysics

Page 40 of 41

Lett. 89 (2010) 34002. [88] Y. Jiao, F.H. Stillinger, S. Torquato, Distinctive features arising in maximally random jammed packings of superballs, Phys. Rev. E. 81 (2010) 41304. doi:10.1103/PhysRevE.81.041304. [89] C. Chen, A. Ozcan, A.O. Yazaydin, B.P. Ladewig, Gas permeation through single-crystal ZIF-8 membranes, J. Memb. Sci. 575 (2019) 209–216. doi:10.1016/j.memsci.2019.01.027. [90] T.T. Moore, W.J. Koros, Gas sorption in polymers, molecular sieves, and mixed matrix membranes, J. Appl. Polym. Sci. 104 (2007) 4053–4059. doi:10.1002/app.25653. [91] R. Krishna, J.M. van Baten, E. García-Pérez, S. Calero, Diffusion of CH4 and CO2 in MFI, CHA and DDR zeolites, Chem. Phys. Lett. 429 (2006) 219–224. doi:10.1016/j.cplett.2006.08.015. [92] E. García-Pérez, J.B. Parra, C.O. Ania, A. García-Sánchez, J.M. van Baten, R. Krishna, D. Dubbeldam, S. Calero, A computational study of CO2, N2, and CH4 adsorption in zeolites, Adsorption. 13 (2007) 469–476. doi:10.1007/s10450-007-9039-z. [93] J. Kärger, H. Pfeifer, F. Stallmach, N.N. Feoktistova, S.P. Zhdanov, 129Xe and 13C PFG n.m.r. study of the intracrystalline self-diffusion of Xe, CO2, and CO, Zeolites. 13 (1993) 50–55. doi:https://doi.org/10.1016/0144-2449(93)90022-U. [94] G.K. Papadopoulos, H. Jobic, D.N. Theodorou, Transport Diffusivity of N2 and CO2 in Silicalite:  Coherent Quasielastic Neutron Scattering Measurements and Molecular Dynamics Simulations, J. Phys. Chem. B. 108 (2004) 12748–12756. doi:10.1021/jp049265g. [95] H. Jobic, W. Schmidt, C.B. Krause, J. Kärger, PFG NMR and QENS diffusion study of nalkane homologues in MFI-type zeolites, Microporous Mesoporous Mater. 90 (2006) 299– 306. doi:https://doi.org/10.1016/j.micromeso.2005.10.020. [96] S. Hwang, R. Semino, B. Seoane, M. Zahan, C. Chmelik, R. Valiullin, M. Bertmer, J. Haase, F. Kapteijn, J. Gascon, G. Maurin, J. Kärger, Revealing the Transient Concentration of CO2 in a Mixed-Matrix Membrane by IR Microimaging and Molecular Modeling, Angew. Chemie Int. Ed. 57 (2018) 5156–5160. doi:10.1002/anie.201713160.

Page 41 of 41

• • • • •

Prediction of gas permeability and selectivity in MMM using multiscale simulations MD simulations reveal interfacial densification at the filler particle surface Optimum filler particle size exists that maximizes permeability in non-ideal MMMs Selectivity is highly sensitive to interfacial resistance at filler particle surface Permeation models considering interfacial defects are best at low filler loading

Declaration of interests ☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. ☒The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

None