Modeling Permeation through Mixed-Matrix Membranes: A Review

Abstract

Over the past three decades, mixed-matrix membranes (MMMs), comprising an inorganic filler phase embedded in a polymer matrix, have emerged as a promising alternative to overcome limitations of conventional polymer and inorganic membranes. However, while much effort has been devoted to MMMs in practice, their modeling is largely based on early theories for transport in composites. These theories consider uniform transport properties and driving force, and thus models for the permeability in MMMs often perform unsatisfactorily when compared to experimental permeation data. In this work, we review existing theories for permeation in MMMs and discuss their fundamental assumptions and limitations with the aim of providing future directions permitting new models to consider realistic MMM operating conditions. Furthermore, we compare predictions of popular permeation models against available experimental and simulation-based permeation data, and discuss the suitability of these models for predicting MMM permeability under typical operating conditions.

1. Introduction

In the last few decades, membrane technologies have attracted increasing attention to be used in a variety of industrial applications, which include gas separation [1,2,3], water desalination [4,5,6], food processing [7], pervaporation [8,9], membrane contactors [10,11], and membrane reactors [12,13]. In many of these applications, membrane technologies are preferred over conventional separation techniques (e.g., distillation, absorption, and adsorption) due to their superior features such as [14,15,16,17]: (i) stable production with high separation efficiency [18]; (ii) low energy consumption with no phase change requirements [15,19]; (iii) simple operation with convenient modular scale-up [20,21,22]; and (iv) small environmental footprint [18,23]. Nevertheless, implementation of membrane technologies in practical applications has been limited by challenges with the engineering of robust materials able to be effective under a variety of operating conditions and environments [21,24], with only polymer membranes currently available in large-scale applications; yet failing to overcome Robeson’s [25,26] trade-off curves between the selectivity and permeability [25,26,27,28].

Different alternatives have been explored to enhance polymer membranes to perform beyond Robeson’s upper bound, including surface modification [29], facilitated transport [30,31], polymer blends [32], and mixed-matrix membranes (MMMs) [33]. Amongst these alternatives, devoted attention to the synthesis of MMMs has been intensified over the last three-decades [21,34,35], with a myriad of studies focusing on novel materials to increase efficiency of CO₂ capture [18,35,36], natural gas purification [37,38,39], water purification [40,41], and olefin/paraffin separation [42,43,44]. Thus, much effort has been devoted to the optimization of MMMs synthesis [34,42,45,46,47,48]; with a number of works even reporting fabrication of defect-free MMMs [37,49,50,51].

Ideally, a mixed-matrix membrane (MMM) consists of a selective inorganic filler phase embedded to continuous polymer matrix [21,35]. In this way, an MMM combines high intrinsic permeability and separation efficiency of advanced molecular sieving materials (e.g., zeolites, carbons, metal-organic frameworks) or nanoscale materials (e.g., carbon nanosheets or nanotubes) with robust processing capabilities and mechanical properties of glassy polymers [23,52]. Consequently, MMMs are expected to have higher efficiency than those based on their polymer counterpart, thus exceeding the trade-off between the permeability and selectivity [15,21,53].

MMMs are commonly prepared either with symmetric or asymmetric structure [17,34,54]. Symmetric MMMs consist of a uniform dense composite film of thickness $20\mathsf{\mu }\mathrm{m}\le \ell \le 100\mathsf{\mu }\mathrm{m}$ [41,46] while asymmetric MMMs comprise a thin selective composite skin layer of thickness $2\mathsf{\mu }\mathrm{m}\le \ell \le 5\mathsf{\mu }\mathrm{m}$ coated on a highly porous non-selective core layer of thickness $50\mathsf{\mu }\mathrm{m}\le \ell \le 300\mathsf{\mu }\mathrm{m}$ [3,55,56,57]. In this way, thicknesses of MMMs are large enough to disregard effects of interfacial entrance and exit barriers on the transport; these barriers which have been shown to significantly decrease the permeant diffusivity in nanoporous materials only when the overall system thickness is $\ell \le 0.1\mathsf{\mu }\mathrm{m}$ [58,59,60]. Such barriers may nevertheless be important at the interfaces of nano-sized fillers and zeolites in nanocomposites, where potential of mean force calculations demonstrate their significance [61]. However, filler size in MMMs is of the order of 0.1–1 μm or larger [34,41], and such barriers are insignificant relative to the internal resistance in the filler particles.

Current models for permeation in MMMs are adaptations of early theories for the transport in heterogeneous media, either following the effective medium approach (EMA) [62,63,64,65] or resistance model approach (RMA) [66,67,68,69]. In such approaches, the effective permeability is usually based on the permeabilities and volume fractions of the MMM constituent phases [21,34,70,71,72], with these phase-specific properties largely assumed constant. Thus, applicability of early RMA/EMA models [62,63,64,66] and later adaptations [65,68,73,74] is often limited to narrow MMM operating conditions (Henry’s law region) and ideal polymer-particle morphologies. Furthermore, although experimental studies on MMMs have shown that deviations from Henry’s law are common under usual operating pressures (1–4 bar) [15,36,46,75], isotherm nonlinearity is incorporated into the permeation models through the Darken or free volume theories while assuming a uniform field based on the mean permeant concentration [71,76,77,78,79], an assumption that needs to be relaxed for further progress.

Over the last decade, increased efforts have been devoted to advance permeation models to integrate effects of filler morphology (e.g., particle size, shape, and agglomeration) [65,80,81] and defects at the particle-polymer interface in the form of a rigidified polymeric [76,77,78,79,82,83,84], and void [69,77,79,85] or pore-blocked [84] regions. However, while considering such non-idealities, these models share the uniform field assumption inherent to the EMA and RMA [69,76,77,86,87]. Thus, effects of the filler morphology together with isotherm nonlinearity are often embedded in a single empirical morphology-related parameter, such as in the Pal [65], Lewis-Nielsen [73], and Higuchi [88] models. Thus, extending the predictive capability of existing models through more rigorous approaches able to be valid over a variety of conditions and systems remains a possibility [15,21,36].

In this work, we review existing approaches for engineering models of permeation through MMMs. To do so, we first introduce the concepts of permeability and selectivity in the context of MMMs and discuss how their mathematical formulation is integrated to existing models for the transport in composite media. Here, we also classify permeation models by approach (i.e., RMA and EMA) and discuss their range of applicability based on their fundamental assumptions. Finally, we compare RMA/EMA models to simulation-based and experimental permeation data for various gases (e.g., $C{O}_2$, $C{H}_4$, $H_2$, $O_2$, $N_2$) in several MMM systems, and provide future directions on how to progress existing models to undertake realistic MMM operating conditions.

2. Gas Transport through Mixed-Matrix Membranes

Gas separation through membranes can take place by different mechanisms [89,90]. Three main diffusion mechanisms have been well-accepted to describe gas transport through membranes [91,92]: (i) Knudsen diffusion; (ii) molecular sieving (molecular diffusion); and (iii) solution-diffusion (sorption-diffusion), with detailed discussion of these transport mechanisms available elsewhere [91]. In general, the diffusion mechanism is assumed to change from solution-diffusion to Knudsen diffusion with increase of the pore size in the membrane material [18,36]. Based on this consideration, transport through inorganic porous membranes has been largely associated with the Knudsen diffusion [35,91], that in membranes based on nanomaterials such carbon molecular sieves (CMSs) [46], zeolitic imidazolate frameworks (ZIFs) [93,94] and metal organic frameworks (MOFs) [95], has been associated with the molecular diffusion, and that through glassy polymers has been associated with the solution-diffusion [30,89,90,91,96,97,98]. While MMMs combine transport principles of both polymer and inorganic membranes, diffusion through them is understood via the solution-diffusion mechanism [21,34]. This mechanism assumes that permeant molecules dissolve (adsorb) on one side of the membrane, diffuse across the membrane and then are released (desorbed) at the other side [14,19,90,99], as depicted in Figure 1 for the $C{O}_2/C{H}_4$ separation.

2.1. Permeability and Selectivity in MMMs

In the solution-diffusion mechanism, the permeant transport is driven by the chemical potential gradient $(\nabla \mu )$ across the membrane, in which $\mu$ is only dependent on the concentration gradient $(\nabla c)$ while the fugacity $(f)$ is assumed uniform across the membrane (cf. Figure 1) [100]. Under these considerations, the permeant flux $(J)$ can be defined as:

$$J=\overline{D}\overline{S}(-\Delta f)/\ell =P(-\Delta f)/\ell ,\Delta f=f_2-f_1$$

(1)

with detailed derivation of Equation (1) found elsewhere [100,101]. Here, $\overline{D}$ and $\overline{S}$ are the concentration-averaged diffusivity and solubility, respectively; and $\Delta f=f_2-f_1$ is the fugacity difference between the retentate $(f_1)$ and permeate $(f_2)$ sides of the membrane (cf. Figure 1), respectively [102]. Further, Equation (1) is usually rearranged as [46]:

$$P=J\ell /(-\Delta f),\Delta f=f_2-f_1$$

(2)

as the permeant flux $(J)$, membrane thickness $(\ell )$, and fugacity difference $(\Delta f)$ can be measured in practice [42,99,103,104]. In Equation (2), the permeability $(P)$ is calculated in $\mathrm{Barrer}$, with [90,105]:

$$1\mathrm{Barrer}=1\times 10^{-10}\frac{{\mathrm{cm}}^3\left(\mathrm{STP}\right)\mathrm{cm}}{{\mathrm{cm}}^2\mathrm{s}\mathrm{cmHg}}=3.348\times 10^{-16}\frac{\mathrm{mol}\mathrm{m}}{{\mathrm{m}}^2\mathrm{s}\mathrm{Pa}}.$$

(3)

For gas mixtures, the permselectivity is used to characterize the MMM separation efficiency, with the permselectivity of species $A$ relative to species $B$ $({\alpha }_{AB}^{*})$ defined as [91,102,103,106]:

$${\alpha }_{AB}^{*}=P_A/P_B=(D_A{S}_A)/(D_B{S}_B)$$

(4)

where ${\alpha }_{AB}^{*}$ is also known as the ideal selectivity. In Equation (4), the permeability of the slower permeant is usually placed in the denominator, and thus ${\alpha }_{AB}^{*}>1$ [99].

2.2. Diffusion and Sorption in MMMs

Based on the permeant flux definition in Equation (1), the permeability can be defined as the product of a kinetic $(\overline{D})$ and a thermodynamic $(\overline{S})$ contribution, with [91,99]:

$$P=\overline{D}\overline{S}$$

(5)

where, assuming the retentate side fugacity to be negligibly small (cf. Figure 1), the mean diffusivity $(\overline{D})$ and solubility $(\overline{S})$ can be expressed as [107,108,109]:

$$\overline{D}=(1/c_1){\displaystyle {\int }_{c_2\approx 0}^{c_1}D(c)dc}$$

(6)

$$\overline{S}=(1/f_1){\displaystyle {\int }_{c_2\approx 0}^{c_1}dc}=c_1/f_1.$$

(7)

respectively. In Equations (6) and (7), $c$ is the permeant adsorbed concentration and $D$ is the Fickian diffusivity, well-accepted to follow the Darken relation as $D=D_o(d\ln f/d\ln c)$ [110,111]. Here, $D_o$ is the permeant mobility, assumed concentration-independent [112,113].

At low pressures, Henry’s law is often adequate to express the gas concentration, with $c=K_{H}f$ [89,109]. Thus, by following Equations (6) and (7), both diffusivity and solubility are concentration- independent, with the permeability in the MMM constituent phases given by [109,110]:

$$P_f=D_{of}{K}_{Hf}$$

(8)

$$P_c=D_{oc}{K}_{Hc}$$

(9)

with subscripts $f$ and $c$ denoting the filler and continuous phases, respectively.

At moderate/high pressures, the solubility in the filler phase is commonly assumed to follow the Langmuir isotherm, leading to:

$${\overline{S}}_f=K_f{c}_f^s/(1+K_f{f}_1)$$

(10)

which often describes well the sorption equilibrium in various porous fillers [42,48,75,114,115,116]. Here, $K_f$ is the affinity constant and $c_f^s$ the capacity in the filler phase. Then, by following Equation (5) while assuming the diffusivity to be concentration-independent, with $\overline{D}=D$ ($\overline{D}\ne D_o$) based on Equation (6) [117], the permeability in the filler phase is given by [99]:

$$P_f=D_f{K}_f{c}_f^s/(1+K_f{f}_1)$$

(11)

where $D_f$ is the permeant diffusivity in the filler. Similarly, a combination of Henry’s law and the Langmuir model is often used to describe the mean polymer solubility, following [75,106]:

$${\overline{S}}_c=K_h+K_c{c}_c^s/(1+K_c{f}_1)$$

(12)

with $K_h$ being the Henry’s law constant, $K_c$ the affinity constant and $c_c^s$ the capacity in the matrix. Thus, by following Equation (5), the permeability in the polymer is given by [117,118]:

$$P_c=D_h{K}_h+D_c{K}_c{c}_c^s/(1+K_c{f}_1)$$

(13)

where $D_h$ and $D_c$ are the diffusivities in the Henry’s law and Langmuir sites, respectively. Equation (13) is known as the dual-mode/partial immobilization model, as the penetrant is assumed fully mobile in the Henry environment and partially mobile in the Langmuir environment [118,119,120,121]. In the next section, the permeant permeabilities in the filler $(P_f)$ and continuous $(P_c)$ phases are largely assumed concentration-independent following Equations (8) and (9) [86,122,123] or based on concentration-averaged solubility and diffusivity via Equations (11) and (13) [15,36,99].

3. Models for Gas Permeation in Mixed-Matrix Membranes

While a universal description of the gas transport through MMM is a complex problem [21,71,87,112], the modeling of permeation through MMMs is largely based on early theories for thermal/electrical conduction of heterogeneous media [62,65,124,125,126]. These theories were extended to MMMs on the basis of the analogy between the thermal/electrical conductivity and the permeability of composite materials [66,81] in the presence of linear flux laws. Based on this analogy, two main approaches have been extensively used to predict the permeability in MMMs: the resistance model approach (RMA) and the effective medium approach (EMA), described in Section 3.1 and Section 3.2, respectively. Besides these early approaches, simulation-based rigorous modeling of MMMs has attracted increased attention in recent years, and thus this latter approach is described in Section 3.3. In what follows, subscripts $c$, $f$, $i$, and $eff$ refer to the continuous phase (polymer matrix), filler phase (selective phase), polymer-filler interface and MMM as a whole, respectively.

3.1. Resistance Model Approach

The resistance model approach (RMA) relies on analogy between the current flow through a series-parallel array of resistors (Ohm’s law) and the permeation rate through a composite membrane (Fick’s law) [66,67,125,127]. Under this consideration, the MMM permeability is inversely proportional to the overall transport resistance, with the fundamental equation of the RMA being [87,128,129,130,131,132]:

$$F=(-\Delta f)/R_{eff},\Delta f=f_2-f_1$$

(14)

where $R_{eff}$ is the overall transport (permeation) resistance and $F=A_{cross}J$ is the permeant flow rate, with $A_{cross}$ the cross-sectional area in the flow direction. Consequently, the permeant flux $(J)$ through the membrane can also be expressed as:

$$J=(-\Delta f)/(A_{cross}{R}_{eff}),\Delta f=f_2-f_1.$$

(15)

On comparing flux definitions in Equations (1) and (15), the resistance $(R_{eff})$ can be equated as [66,68]:

$$R_{eff}=\ell /(P_{eff}{A}_{cross}).$$

(16)

In this way, if an expression for the equivalent resistance $(R_{eff})$ is known, the MMM permeability $(P_{eff})$ can be calculated via Equation (16) [125,132].

A number of permeation models have been proposed for MMMs [68,69,72,76,82,87,133], based on Equations (15) and (16), as listed in

Table 1. Models based on the resistance model approach.

Model	Key Equations
Series [47]	$P_{eff}=\frac{P_c{P}_f}{P_c{\varphi }_f+P_f(1-{\varphi }_f)}$	(17)
Parallel [47]	$P_{eff}=P_f{\varphi }_f+P_c(1-{\varphi }_f)$	(18)
Te Hennepe [130]	$P_{eff}=P_c{\left[(1-{\varphi }_f^{\frac{1}{3}})+\frac{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{\varphi }_f^{\frac{1}{3}}{P}_c}{P_c(1-{\varphi }_f)+\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.P_f{\varphi }_f}\right]}^{-1}$	(19)
Cussler [72]	$P_{eff}=P_c{\left[(1-{\varphi }_f)+{\left(\frac{P_f}{{\varphi }_f{P}_c}+\frac{4(1-{\varphi }_f)}{{\lambda }_f^2{\varphi }_f^2}\right)}^{-1}\right]}^{-1}$	(20)
Ebneyamini [68]	$P_{eff}=\tau \left[(1-{\varphi }^{2/3})P_c+\frac{P_c{P}_f{\varphi }^{2/3}}{{\varphi }^{1/3}{P}_c+(1-{\varphi }^{1/3})P_f}\right]$	(21)
KJN [69]	$P_{eff}^{Oriented}=P_c{\left[\left(1-\frac{\cos \theta }{\cos \theta +{\lambda }_f\sin \theta }{\varphi }_f\right)+\frac{P_c}{P_f}\left(\frac{1}{\cos \theta +{\lambda }_f\sin \theta }\right){\varphi }_f\right]}^{-1}$	(22)
	$P_{eff}^{Random}=\frac{\pi }{2}{P}_c{\left[{\displaystyle {\int }_0^{\frac{\pi }{2}}\frac{P_c}{P_{eff}^{Oriented}(\theta )}d\theta }\right]}^{-1}$	(23)

, with the most popular models considering platelet or cubic filler particles [72,130,133] and few newly developed models considering tubular particles [69,76,82]. The simplest RMA models idealize the MMM as a two-phase laminated composite comprising multiple sheets of polymer and selective material alternated in series or in parallel to the flow direction [70,134], as depicted in Figure 2a,b, respectively. Figure 2 also depicts electrical circuit analogues used to calculate the overall transport resistance of both composites.

Following the electrical circuit analog in Figure 2 and the above definition for the permeation resistance in Equation (16), the permeability for the multilayer composite in series yields Equation (17) in Table 1, while that of a multilayer composite in parallel yields Equation (18) [47,70]. The series model in Equation (18) is assumed to provide the lower bound for the permeability of a given penetrant in an ideal MMM [135]. Alternatively, the parallel model in Equation (18) is assumed to provide to the upper bound for the effective permeability of a given penetrant in an ideal MMM [47,70].

In addition to the two-resistance based models in Equations (17) and (18), more complex models, including three-resistances or more, have been proposed following the RMA [68,69,72,82,87,130]. In these models, the additional diffusion resistances are intended to account for tortuosity effects in the permeant diffusion path when there are large differences in permeabilities amongst the filler and polymer phases [68,87] or defects in MMM structure [69,82,87]. Amongst existing RMA models, those of Te Hennepe [130] in Equation (19) and Cussler [72] in Equation (20), based on three-resistance circuit analogs, have widely been applied to zeolite-polymer MMMs [114,128,136].

Te Hennepe et al. [130] considered the one-dimensional transport in zeolite-rubber MMMs. They idealized the MMM as a lamella containing composite layers [131], in which each composite layer comprised two regions. The first region consisted of polymer and the second one of polymer and zeolite particles (mixed-region). In this model, the polymer region was assumed in series with parallel resistances of the second mixed-region [130,131], which led to Equation (19) in Table 1. Alternatively, Cussler [72] considered two-dimensional transport in the MMM. He assumed the resistance of the polymer region in series with that of a second mixed-region, similar to Te Hennepe et al. [130]. However, transport in the mixed-region was assumed to occur in the permeation direction through filler phase and perpendicular to the permeation direction through polymer phase [15,36]. This assumption led to Equation (20) in Table 1, in which ${\lambda }_f=w_f/{\ell }_f$ is the aspect ratio of the filler phase with $w_f$ and ${\ell }_f$ being the flake width and thickness [21,136].

Recently, Ebneyamini et al. [68] proposed a semi-empirical four-resistance model for ideal MMMs comprising cubical particles. To do so, an empirical correction factor $(\tau )$ was introduced to a one-dimensional four-resistance model. The final model is given by Equation (21) in Table 1, and referred here as the Ebneyamini model. In this model, $\tau$ was estimated via simulation of the 3D particle-polymer system and adjusted to follow Langmuir-type equations. Thus, $\tau$ is assumed to accommodate tortuosity effects arising from large differences amongst the MMM constituent phase permeabilities.

Modeling MMMs with tubular filler has received less attention than those having cubic or platelet fillers, with only few studies [69,87] developing RMA models for nanotube-MMMs. The first of these RMA models was proposed by Kang et al. [69], who accommodated the orientation of tubular fillers in the calculation of the overall transport resistance. The final model was named by authors as the Kang-Jones-Nair (KJN) model, and for uniformly oriented fillers is given by Equation (22) in Table 1. For randomly oriented fillers, the KJN model is rewritten as Equation (23) in Table 1, with $P_{eff}^{Oriented}(\theta )$ in Equation (23) following Equation (22). Here, $\theta \in [0,\pi /2]$ is the orientation of the tubular filler, measured with respect the permeation direction, and ${\lambda }_f=d_f/{\ell }_f$ is the aspect ratio of the tubular filler with $d_f$ and ${\ell }_f$ being the diameter and length of the cylindrical particle, respectively. The predictions of the KJN model are always lower than those based on the series model, with the KJN model simplifying to the series model when $\theta =0$. Figure 3 depicts a comparison of the permeability $(P_{eff})$ profiles based on models of Table 1 with ${\alpha }_{fc}=P_f/P_c=10$ in all models and ${\lambda }_f=0.25$ in the Cussler and KJN models. Figure 3 also depicts the MMM structure assumed by each model.

3.2. Effective Medium Approach

The crux of the effective medium approach (EMA) lies in the substitution of a given composite system by an equivalent effective homogeneous one having the properties of the composite [137,138,139]. The resulting effective composite properties are generally functions of the volume fraction and permeabilities of the composite constituent phases, similar to RMA-based models (cf. Section 3.1). However, EMA differs from RMA in the way the filler phase is considered within the composite. While most RMA models assume regular distributions (e.g., simple cubic or body centered cubic lattices) of platelet and/or cubic particles [47,68,131,140], EMA models consider random distributions of spherical inclusions [64,65,141,142,143]. For ease of analysis, we here classify EMA models in two main groups. The first group corresponds to EMA models following Maxwell’s theory and second to those following Bruggeman’s theory, with models associated with each theory described in Section 3.2.1 and Section 3.2.2, respectively.

3.2.1. Maxwell Theory

Maxwell [62] calculated the electrical conductivity of infinitely diluted cluster of particles embedded in an infinite matrix [138]. To do so, he assumed that the far field potential of this cluster was equivalent to that of a homogeneous sphere having the original composite volume [65,135,137], as depicted in Figure 4. Based on this assumption, Maxwell defined the composite conductivity as that of the homogeneous sphere [21,141,144,145], which led to Equation (24) in

Table 2. Models based on Maxwell’s theory.

Model	Key Equations
Maxwell [62]	$P_{eff}=P_c\left[\frac{1+2{\beta }_{fc}{\varphi }_f}{1-{\beta }_{fc}{\varphi }_f}\right]=P_c\left[\frac{P_f+2P_c-2{\varphi }_f(P_c-P_f)}{P_f+2P_c+{\varphi }_f(P_c-P_f)}\right]$	(24)
Maxwell-Wagner-Sillars model [81]	$P_{eff}=P_c\frac{{\lambda }_f{P}_f+(1-{\lambda }_f)P_c-(1-{\lambda }_f){\varphi }_f(P_c-P_f)}{{\lambda }_f{P}_f+(1-{\lambda }_f)P_c+{\lambda }_f{\varphi }_f(P_c-P_f)}$	(25)
Chiew-Glandt [141]	$P_{eff}=P_c\left[\frac{1+2{\beta }_{fc}{\varphi }_f+(K_2-3{\beta }_{fc}^2){\varphi }_f^2}{1-{\beta }_{fc}{\varphi }_f}\right]$	(26)
	$K_2=a+b{\varphi }_f^{\frac{3}{2}}$	(27)
	$a=-0.002254-0.123112{\beta }_{fc}+2.93656{\beta }_{fc}^2+1.6904{\beta }_{fc}^3$	(28)
	$b=0.0039298-0.803494{\beta }_{fc}-2.16207{\beta }_{fc}^2+6.48296{\beta }_{fc}^3+5.27196{\beta }_{fc}^4$	(29)
MB-B model [152]	$\frac{1}{R^{i-1}}\frac{d}{dR}\left[R^{i-1}{P}_m(R)\frac{d{C}_m}{dR}\right]=0$	(30)
	${\overline{\varphi }}_f(R)=\frac{3}{2r_o^3}{\displaystyle {\int }_0^{r_o}{\displaystyle {\int }_0^{\pi }{\varphi }_f\left(R\prime (R,r,\theta )\right)r^2\sin \theta d\theta dr}}$	(31)
Lewis-Nielsen [73]	$P_r=\frac{1+2{\beta }_{fc}{\varphi }_f}{1-\beta {\varphi }_{fc}{\psi }_m}$	(32)
	${\psi }_m=1+\left[(1-{\varphi }_m)/{\varphi }_m^2\right]{\varphi }_f$	(33)
Higuchi [88]	$P_{eff}=P_c\left[1+\frac{3{\varphi }_f{\beta }_{fc}}{1-{\varphi }_f{\beta }_{fc}-K_H(1-{\varphi }_f){\beta }_{fc}^2}\right]$	(34)
Felske model [77]	$P_{eff}=P_c\left[\frac{2(1-{\varphi }_{fi})+(1+2{\varphi }_t)(\eta /\gamma )}{(2+{\varphi }_{fi})+(1-{\varphi }_t)(\eta /\gamma )}\right]$	(35)
	$\eta =(2+{\delta }^3){\alpha }_{fc}-2(1-{\delta }^3){\alpha }_{ic}$	(36)
	$\gamma =(1+{\delta }^3)-(1-{\delta }^3){\alpha }_{fi}$	(37)
	${\varphi }_{fi}=\frac{{\varphi }_f^N}{{\varphi }_f^N+[(1-{\varphi }_f^N)/{\delta }^3]}$	(38)
Pseudo-two-phase Maxwell [85]	$P_{eff}=P_c\left[\frac{P_{eff}^f+2P_c-2{\varphi }_{fi}(P_c-P_{eff}^f)}{P_{eff}^f+2P_c+{\varphi }_{fi}(P_c-P_{eff}^f)}\right]$	(39)
	$P_{eff}^f=P_i\left[\frac{P_f+2P_i-2{\varphi }_s(P_i-P_f)}{P_f+2P_i+{\varphi }_s(P_i-P_f)}\right]$	(40)
	${\varphi }_s=1/{(1+{\ell }_i/r_o)}^3$	(41)

, with ${\beta }_{fc}=({\alpha }_{fc}-1)/({\alpha }_{fc}+2)$ and ${\alpha }_{fc}=P_f/P_c$. Further, because the Maxwell model disregards particle interaction, it is only applicable to dilute suspensions $({\varphi }_f\le 0.2)$.

The Maxwell model was later extended to spheroids by Wagner and Sillars [126]. To do so, the conduction problem was reformulated with oriented spheroidal inclusion, with these spheroids oriented along the axis of the potential difference [81,126]. This consideration led to Equation (25) in Table 2, referred as the Maxwell-Wagner-Sillars model [81]. Here, ${\lambda }_f\in [0,1]$ is the particle shape factor, and for prolate spheroids ${\lambda }_f\in [0,1/3]$ while for oblate spheroids ${\lambda }_f\in [1/3,1]$. In the limits, when ${\lambda }_f=0$, ${\lambda }_f=1/3$, and ${\lambda }_f=1$, Equation (25) reduces to the parallel model, series model (cf. Table 1) and Maxwell model, respectively. The Maxwell-Wagner-Sillars model is only applicable to diluted ellipsoid dispersions $({\varphi }_f\le 0.2)$, similar to the Maxwell model [3,21,41,146].

Several attempts have been made to extend Maxwell’s equation to concentrated composites $({\varphi }_f>0.2)$ [64,65,73,141,142], of which, the one with the most significant results is that of Jeffrey [142,147]. He showed that relative conductivity $(P_r)$of a dispersion can be expressed in series of ${\varphi }_f$ following the form $P_r=P_{eff}/P_c=1+K_1{\varphi }_f+K_2{\varphi }_f^2+\cdots$ [141,148], in which each $K$-term in the series $(K_1,K_2,\dots ,K_n)$ accommodates the interaction of successively larger sets of particles [141,149]. Further, Jeffrey [142,147] demonstrated that only the first-order term of the Maxwell model is exact, where Equation (24) yields $P_r=P_{eff}/P_c=1+3{\beta }_{fc}{\varphi }_f+\mathcal{O}({\varphi }_f^2)$ with $K_1=3{\beta }_{fc}$ when $\varphi \to 0$. This limiting equation was later used by Bruggeman [139] and Pal [65], as described in Section 3.2.2.

Although Jeffrey [142,147] hypothesized that $K_2$ in the series was dependent on ${\beta }_{fc}$ and ${\varphi }_f$, his final equation corresponds to the low-density limit of interacting spheres [141]. Thus, Jeffrey’s model provides very similar predictions to the Maxwell model [150]. Later, Chiew and Glandt [141] estimated $K_2$ in the series using pair-correlation functions of hard-sphere fluid simulations, which led to Equation (26) in Table 2. In this work, the resulting values of $K_2$ were tabulated as function of ${\alpha }_{fc}$ and ${\varphi }_f$ [141,149]. Later, Gonzo et al. [145] fitted Chiew and Glandt’s results for $K_2$, which led to Equation (27) in Table 2. Further, the Chiew-Glandt model corresponds to the exact solution to the series truncated after $K_2$-term [141,149], with the model applicable to moderate/high filler loadings $({\varphi }_f\le 0.645)$ [141,145,151].

Recently, Monsalve-Bravo and Bhatia [151,152,153] used the Chiew-Glandt model in conjunction with the one-dimensional transport equation for the permeant to describe the permeability for various gases in flat and hollow fiber MMMs. The semi-analytical model is given by Equation (30) in Table 2, referred to as the MB-B model. Here, $i=1$ for a flat MMM and $i=2$ for a hollow fiber MMM while $C_m$ is the position-dependent pseudo-bulk concentration in the MMM, with boundary conditions $C_m=C_{m1}=f_1/R_g{T}_g$ at $R=R_1$ and $C_m=C_{m2}=f_2/R_g{T}_g$ at $R=R_1$. In Equation (30), $P_m(R)$ is the effective local MMM permeability, estimated using the Chiew-Glandt model in Equation (26), in which the constituent phase permeabilities ($P_f$ and $P_c$) are concentration-dependent via the Darken model [154,155]. Further, the locally averaged filler volume fraction $({\overline{\varphi }}_f)$ in Equation (31) averages the filler volume fraction $({\varphi }_f)$ over the particle volume at a given membrane location, with $r_o$ being the filler particle radius and $R\prime (R,r,\theta )$ the location in the particle relative to the position $R$ in the MMM. Thus, the MB-B model incorporates effects of particle size [151], membrane geometry [152] and isotherm nonlinearity [153] in the calculation of the MMM permeability. The MB-B model reduces to the Chiew-Glandt model when $r_o/\ell \to 0$.

As an alternative to the exact second order solution of Chiew and Glandt [141], several empirical modifications of the Maxwell model [80,156,157] have been proposed for concentrated composites. These models empirically embed packing-related effects [158] and variation in filler properties, such as particle agglomeration, size, and shape [71,80,86], within a single parameter in the model [80,86,153], with the Lewis-Nielsen model in Equation (32) being one of most popular of these models. In this model, packing-related effects are accommodated via the maximum filler volume fraction $({\varphi }_m)$ [124,159], with the Lewis-Nielsen model simplifying to that of Maxwell when ${\varphi }_m\to 1$. Alternatively, Higuchi [88] introduced empirical parameter, $K_H$, to the Maxwell model to account for particle-particle interactions and asphericity effects arising from particle shape variation, with the Higuchi model given by Equation (34) [70,85] and $0\le K_H\le 0.78$ for spherical particles [160]. Thus, Equation (34) reduces to the Maxwell model when $K_H=0$ [74]. Figure 5 compares the permeability $(P_{eff})$ profiles based on ideal models in Table 2, with ${\alpha }_{fc}=P_f/P_c=10$ and specific model parameters values listed in the legend.

Based on Maxwell’s theory, several models have been proposed to account for non-ideal polymer-particle morphologies [77,78,84,86,160], either by: (i) solving the transport problem in an MMM comprising spherical core-shell inclusions [77,86] or (ii) assuming the particle-interface system as a pseudo-phase dispersed in the polymer matrix [71,85]. In the first group of non-ideal models are found that of Felske [77], and later Pal’s adaptation [86]. In the second group are found the pseudo-two-phase Maxwell model [85], and its adaptations [84,161]. The Felske model in Equation (35) [77] is the exact first order solution to the transport problem through a dispersed composite comprising non-interacting core-shell particles [160]. In this model, ${\varphi }_{fi}={\varphi }_f+{\varphi }_i$ is the volume fraction of the total dispersed phase following Equation (38), with $\delta =({\ell }_i+r_o)/r_o$ and ${\varphi }_f^N$ the nominal filler volume fraction [70,86,160]. Further, ${\alpha }_{fc}=P_f/P_c$, ${\alpha }_{ic}=P_i/P_c$, and ${\alpha }_{fi}=P_f/P_i$ in Equations (36) and (37). The Felske model is only applicable to dilute suspensions $({\varphi }_{fi}\le 0.2)$ [71,160,162].

In the pseudo-two-phase Maxwell model in Equation (39) [85], it is assumed that the three-phase composite can be idealized as pseudo two-phase composite [163,164], with the polymer matrix being one phase and the combined filler-interface system being the other phase (i.e., pseudo dispersed filler phase) [21,71]. In Equation (39), ${\varphi }_{fi}$ is the volume fraction of total dispersed phase, given by Equation (38) [70]. Further, $P_{eff}^f$ is the permeability of the combined filler-interface system, given by Equation (40). Here, ${\varphi }_s$ is the volume fraction of the filler in the filler-interface composite, and following Equation (41) [70,85,163]. Further, $P_i$ is assumed well described by Knudsen diffusion mechanism when interfacial voidage is considered [70,163]. Alternatively, when polymer rigidification is considered at the filler surface $P_i=P_c/\sigma$, with $\sigma \in [3,4]$ being the chain immobilization factor; an empirical parameter used to differentiate the rigid polymer permeability from that in the bulk polymer [21,163,164]. Later, Li et al. [84] modified the pseudo-two-phase Maxwell model to accommodate both partial pore blockage and polymer rigidification. In this later work, the Maxwell model is accordingly applied three times to include both effects.

3.2.2. Bruggeman’s Theory

The Bruggeman’s theory proceeds from the premise that the field of neighboring particles can be taken into account by randomly adding the dispersed phase incrementally while considering the surrounding medium as the existing composite with effective transport properties at each stage [64,65,139]. This concept of incremental homogenization is illustrated in Figure 6, in which $P_{eff}^1$, $P_{eff}^2$ and $P_{eff}$ are the MMM permeabilities at different stages of the homogenization process.

Bruggeman [64] based his expression for the dielectric permeability of dispersed composites on the assumption that the Maxwell model described well the composite permeability in the diluted limit [65,139], as depicted in Figure 6. To do so, He considered that the differential increment in the dielectric permeability $(dP)$ resulting from addition of new particles was well described by $P_r=P_{eff}/P_c=1+3{\beta }_{fc}{\varphi }_f$, that upon substitution of $P_c\to P$, $P_m\to P+dP$, and ${\varphi }_f\to d{\varphi }_f/(1-{\varphi }_f)$ leads to Equation (42) in

Table 3. Models based on Bruggeman’s theory.

Model	Key Equations
Bruggeman [64]	${\left[\frac{P_{eff}}{P_c}\right]}^{\frac{1}{3}}\left[\frac{{\alpha }_{fc}-1}{{\alpha }_{fc}-(P_{eff}/P_c)}\right]={\left[1-{\varphi }_{{}_f}\right]}^{-1}$	(42)
Pal [65]	${\left[\frac{P_{eff}}{P_c}\right]}^{\frac{1}{3}}\left[\frac{{\alpha }_f-1}{{\alpha }_f-(P_{eff}/P_c)}\right]={\left[1-\frac{{\varphi }_f}{{\varphi }_m}\right]}^{-{\varphi }_m}$	(43)
Pseudo-two phase Bruggeman model [78]	${\left[\frac{P_{eff}}{P_c}\right]}^{\frac{1}{3}}\left[\frac{(P_{eff}^f/P_c)-1}{(P_{eff}^f/P_c)-(P_{eff}/P_c)}\right]={\left[1-{\varphi }_{fi}\right]}^{-1}$	(44)
	${\left[P_{eff}^f/P_i\right]}^{\frac{2}{3}}=1-{\varphi }_s,P_f\ll P_i$	(45)
	${\left[\frac{P_{eff}^f}{P_i}\right]}^{\frac{1}{3}}\left[\frac{(P_f/P_i)-1}{(P_f/P_i)-(P_{eff}^f/P_i)}\right]={\left[1-{\varphi }_s\right]}^{-1},P_f\ge P_i$	(46)
Pseudo-two phase Pal model [71]	${\left[\frac{P_{eff}}{P_c}\right]}^{\frac{1}{3}}\left[\frac{(P_{eff}^f/P_c)-1}{(P_{eff}^f/P_c)-(P_{eff}/P_c)}\right]={\left[1-\frac{{\varphi }_{fi}}{{\varphi }_m}\right]}^{-{\varphi }_m}$	(47)
	${\left[\frac{P_{eff}^f}{P_i}\right]}^{\frac{1}{3}}\left[\frac{(P_f/P_i)-1}{(P_f/P_i)-(P_{eff}^f/P_i)}\right]={\left[1-\frac{{\varphi }_s}{{\varphi }_m}\right]}^{-{\varphi }_m}$	(48)

after integration [138,165]. The Bruggeman model is assumed to be applicable to moderate filler loadings $(0\le {\varphi }_f\le 0.35)$ [135,139,146].

Similar to Bruggeman [64], Pal [65] assumed that the increment $dP$ resulting from the addition of new particles was well-described by $P_r=P_{eff}/P_c=1+3{\beta }_{fc}{\varphi }_f$, with $P_c\to P$ and $P_m\to P+dP$; however, ${\varphi }_f\to d{\varphi }_f/(1-{\varphi }_f/{\varphi }_m)$ was used instead of ${\varphi }_f\to d{\varphi }_f/(1-{\varphi }_f)$, which leads to Equation (43) in Table 3 after integration [86,145]. In this model, ${\varphi }_m$ has the same connotation as in the Lewis-Nielsen model (cf. Section 3.2.1) [80,86,151,153]. The Pal model reduces to Bruggeman’s result when ${\varphi }_m=1$. Further, the Pal model always predicts higher values for the permeability $(P_{eff})$ than that of Bruggeman for a given system.

Analogous to the pseudo-two-phase Maxwell [85], both Bruggeman and Pal models have been extended to describe non-ideal MMMs. The former of these adaptations is known as the pseudo-two phase Bruggeman model in Equation (44) [78], with the model describing the effect on the MMM permeability of a voided or rigidified interfacial region at the filler particle surface [79]. In this model, ${\varphi }_{fi}$ is given by Equation (38) [70] and ${\varphi }_s$ given by Equation (41) [70,85,163]. Further, $P_{eff}^f$ is given by Equation (45) for voided interphase$(P_f<<P_i)$ [81] and by Equation (46) for rigidified interface [166]. Similarly, the pseudo-two phase Pal model in Equation (48) [71] accounts for the effect of polymer rigidification on the MMM permeability [79,85,160,164]. In this way, ${\varphi }_{fi}$ and ${\varphi }_s$ are given by Equations (38) and (41), respectively [70]. Here, $P_{eff}^f$ and $P_i$ have the same connotation as in the pseudo-two-phase Maxwell model (cf. Section 3.2.1) [79]. Recently, Idris et al. [79] modified the pseudo-two-phase Bruggeman model to account for both effects of polymer rigidification and presence of voids at the particle surface. In this later work, the Bruggeman model is accordingly applied three times to include both effects.

Finally, Bruggeman also developed a symmetric theory for the transport in composites [126,137]. However, models such as those of Landauer [63] and Böttcher [167], associated with this symmetric theory, are not discussed here. This is because the symmetric Bruggeman’s theory considers the composite to be a random assembly of spherical particles of different materials [138,168], in which all components in the composite are continuous in the medium [141]. This condition is not met in MMMs, where one phase is preferentially assumed dispersed in the other.

3.3. Simulation-Based Rigorous Modeling Approach

The simulation-based rigorous modeling approach (SMA) is based on numerical solution of coupled partial differential equations (PDEs) describing the transport through the MMM via the finite-element method (FEM) [20,151,169], or any other suitable method (e.g., finite volume method or boundary element method), to discretize the 3D computational system [123]. Thus, the SMA is based on the assumption that the steady-state transport through the MMM is well-described by the continuity equation [151,170], as:

$$\nabla \cdot J=0$$

(49)

where $J$ is the permeant steady-state flux, with the Fick’s law describing the flux in both dispersed $(J_f)$ and continuous $(J_c)$ phases, as [153]:

$$J_f=D_f(C_f)(-\nabla C_f)$$

(50)

$$J_c=D_c(C_c)(-\nabla C_c).$$

(51)

respectively. In Equations (50) and (51), ∇C and D are the concentration gradient and diffusivity in a given phase, respectively. In this way, Equations (49)–(51) are solved in the 3D MMM [151], in which the resulting steady-state flux is used to calculate the MMM permeability via Equation (2) [152].

The SMA offers several benefits in comparison to earlier RMA and EMA approaches, amongst which the main advantage is that this approach can easily incorporate the permeability dependence on the concentration field [153] and finite-system size effects [151], largely disregarded in the former approaches. In this way, the SMA can accommodate effects on the permeability of intrinsic system properties such filler particle size [151] and shape [171,172], isotherm nonlinearity [153], and membrane geometry [152,169]. These effects are empirically treated in most RMA/EMA models [20], such in the Ebneyamini [68], Higuchi [88], and Pal [65] models.

The effect of filler size on the MMM permeability has been the focus of several studies [20,151,169], with Singh et al. [20] being the first to investigate such an effect for spherical particles and MMMs operating in the Henry’s law region. In this work, the relative permeability $(P_{eff}/P_c)$ was found independent of the particle size in the range of sizes investigated. However, later studies [151,169] suggested that the MMM permeability decreased with increase of particle size, and associated this depletion in the permeability with decrease of the specific polymer-filler interfacial area with increase of particle size at fixed volume fraction [151,173]. In these latter studies, the discrepancy amongst studies on the effect of particle size on the MMM permeability was associated with poor mesh quality in the FEM implementation in the former study [169], where the permeability was also found sensitive to changes of the product ratio $K_f{D}_f/K_c{D}_c$; inconsistent with definition in Equation (5) (cf. Section 2.2) [174]. Further, for MMMs operating at low pressure, other investigations have been focused on tubular [172] and layered fillers [171], with these investigations indicating sensitivity of the MMM permeability to filler shape, orientation, and packing.

Recently, the effect of isotherm nonlinearity was evaluated in MMMs with spherical fillers [153]. In this study, the MMM permeability was found more sensitive to isotherm nonlinearity in the filler phase than in the continuous phase. Similar to earlier studies in the Henry’s law region [151,169], increase of filler particle size was found to decrease the MMM permeability. Further, comparison of the simulation predictions to EMA models suggested that the Chiew-Glandt model in Equation (26) describes well the effective permeability of ideal MMMs when the relative particle size is negligible $(r_o/\ell \to 0)$, as in such case the effects of isotherm nonlinearity and finite-system size effects are negligible in the system. Figure 7 depicts a comparison of several EMA models to simulation results [151] for small filler particle size, and showing the Chiew-Glandt model to match well the simulation results, with the smallest percentage deviation $(4.7\%)$.

Following the SMA, the effect of the membrane geometry has been evaluated [153,169]. Yang et al. [169] were the first to evaluate such an effect, by comparing the permeability of full-scale hollow fiber and flat-dense MMMs. By considering randomly distributed spherical inclusion and constant filler and polymer diffusivities, they found the hollow-fiber MMMs to have higher permeabilities than flat-dense MMMs. The increased efficiency of hollow fiber MMMs was associated with decrease in the length of the permeant diffusion pathways in the hollow fiber relative to those in a flat MMM. However, a more recent study found opposite behavior [152], by comparing MMMs having the same filler particle size, thickness, and volume, with the permeability in the hollow fiber found lower than in the flat MMM. Further, this tendency was associated with asymmetric filler phase distribution at the MMM ends, and arising from curvature change in the hollow fiber. The discrepancy amongst the different studies was associated with inaccuracies in the FEM implementation, which also led to the incorrect finding that the effective permeability was sensitive to variation of the ratio $K_f/K_c$ and $D_{of}/D_{oc}$ at constant $K_f{D}_{of}/K_c{D}_{oc}$ [149].

4. Predicting the Effective Permeability in Mixed-Matrix Membranes

Although there are a myriad of models for permeation in MMMs [62,68,80,85,147,151,175], with the popular ones described in Section 3, only a few of these match experimental permeation data [120,145,153,176]. While ideal models are often used as reference in practice [45,93,99], they are commonly shown to poorly describe experimental MMM permeabilities [71,128,166]. Alternatively, models incorporating interfacial non-idealities have been largely shown to match experimental data [70,71,78,79,161]. However, properties of the interfacial layer are always empirically fitted in these models [160,166,177]. When compared to experimental data, RMA and EMA models are used to either estimate the permeability in: (i) the MMM or (ii) any of its constituent phases by fitting MMM permeability data. In the first case, filler and continuous phase permeabilities are experimentally measured, and EMA/RMA are used to estimate the permeability in the MMM [99,107,110]. In the second case, the MMM and continuous phase permeabilities are experimentally measured, and EMA/RMA models are used to estimate the filler phase permeability [42,81], by fitting MMM permeation data. Here, we discuss the use of EMA/RMA models to predict filler permeabilities in Section 4.1 and MMM permeabilities in Section 4.2.

4.1. Estimation of the Filler Phase Permeability through EMA and RMA Models

Amongst the studies estimating the filler permeability through EMA/RMA models [42,81,128], that of Bouma et al. [81] estimated the relative permeability of $O_2$ in the filler phase $({\alpha }_{fc}=P_f/P_c)$ in an MMM comprising co-polyvinylidene fluoride-hexafluoropropene (co-PVDF) as continuous phase and a nematic liquid crystalline mixture (E7) as dispersed phase. In this work, Bouma et al. [81] found the Maxwell model to lead to ${\alpha }_{fc}=P_f/P_c=15.0\pm 8.1$ while that of Bruggeman to ${\alpha }_{fc}=P_f/P_c=9.2\pm 3.8$. Based on the relative error from each model fit, they concluded that the Bruggeman model better fitted the MMM permeabilities at higher filler concentrations.

Later, Gonzo et al. [145] found that the Chiew-Glandt model better fitted Bouma’s experimental permeation data [81] as well as in other gases (e.g., $C{O}_2$, $He$, $H_2$, $O_2$, $N_2$) for various MMMs upon comparison of the Maxwell, Bruggeman, Chiew-Glandt, Higuchi and Landauer models (cf. Section 3.2). Furthermore, they found the Chiew-Glandt model to describe well the MMM permeability at both low and moderate filler volume fractions $(0\le {\varphi }_f\le 0.4)$ [145]. Gonzo et al. [145] associated the overall success of the Chiew-Glandt model with the way particle-particle interactions are accommodated in this model [141,149], as the asymmetric integration technique used by Bruggeman [64] considers that newly added spheres are too dilute to interact among themselves, and interact only with the previously homogenized system [135,139,158]. This assumption may be inaccurate at high filler loadings.

In a similar fashion to Bouma et al. [81], Erdem-şenatalar et al. [128] estimated the filler phase permeability $(P_f)$ by fitting experimental permeation data of various gases (e.g., $C{O}_2$, $O_2$, $N_2$, $C{H}_4$, $n-C_4{H}_{10}$, and $i-C_4{H}_{10}$) in Zeolite/Polymer MMMs (e.g., Zeolite 4A, 5A, 13X, and Silicalite-1) to the Series, Parallel, Te Hennepe, Maxwell, and Landauer models. They found large differences in $P_f$ estimations amongst models for all gases and associated poor $P_f$ predictions with the limited capability of RMA/EMA models to account for interfacial defects. However, here we associate these differences in $P_f$ with the use of RMA/EMA models beyond their range of applicability, as most of these models are only suitable at low filler loadings.

Figure 8 depicts a comparison of the experimental effective permeability and that predicted by RMA/EMA models for one of the experimental data sets used by Erdem-şenatalar et al. [128]. However, we here calculated $P_f$ via nonlinear regression of each model with the experimental permeation data in the range $0\le {\varphi }_f\le 0.2$, as only at low volume fractions RMA/EMA models provide comparable predictions of $P_{eff}$ for a given ${\alpha }_{fc}=P_f/P_c$ (cf. Figure 3 and Figure 5). In Figure 8, $P_f$ predictions and percentage deviations between RMA/EMA models and experimental permeabilities are shown in the legend. Here, the experimental permeation data were originally reported by Duval et al. [178] for $C{O}_2$ in silicalite-1/EPDM MMMs, in which size of zeolite particles ranged between $1-5\mathsf{\mu }\mathrm{m}$, membrane thicknesses between $50-200\mathsf{\mu }\mathrm{m}$, and no morphological defects were reported in the MMMs structure.

In Figure 8, we used the Chiew-Glandt model in place of the Landauer model, as the Landauer model has been shown not suitable for MMMs [81,141]. Here, the series model $(P_f\to \infty )$ under-predicts the experimental permeabilities, with deviations of about $25\%$. This tendency is expected, as the series model corresponds to the lower bound for the permeability in ideal MMMs (cf. Section 3.1). Further, the filler phase permeability varies from 298 < P_f < 896 Barrer, with the Chiew-Glandt model providing the best match to the experimental MMM permeabilities (deviation of $6.7\%$). This suggests that the Chiew-Glandt models describes well the MMM permeability when membrane structure has no associated defects, and particle size is negligible in comparison to the membrane thickness.

4.2. Prediction of the MMM Permeability through EMA and RMA Models

Among the studies comparing simulation-based [20,151,152,169,171,172] or experimental permeabilities of MMMs to EMA/RMA models [45,120,145,153,166,179], that of Monsalve-Bravo and Bhatia [151] compared simulation-based MMM permeabilities to various EMA models. In this work, they showed the Chiew-Glandt model to describe well the simulation-based MMM permeabilities when the ratio of the particle size (radius) to the membrane thickness is small (cf. Figure 7), and that the MB-B model (cf. Section 3.2.1) matched well the simulation-based MMM permeabilities when relative filler particles sizes were in the range$0.004\le r_o/\ell \le 0.16$.

Figure 9 depicts a comparison of the simulations-based MMM permeabilities reported in [151] and the MB-B, Pal, Bruggeman, Lewis-Nielsen, Chiew-Glandt, and Maxwell models, with ${\alpha }_{fc}=P_f/P_c=100$ in all models, $r_o/\ell =0.040$ in the MB-B model, and ${\varphi }_m=0.645$ in the Pal and Lewis-Nielsen models. Further, percentage deviations between EMA models and simulation results are shown in the legend. In Figure 9, the MB-B model matches the simulation-based permeabilities, with a percentage deviation of $2.0\%$. Further, Monsalve-Bravo and Bhatia [151] showed that increase of filler particle size decreased the effective permeability of MMMs and associated this behavior with both depletion of the filler volume fraction at the membrane ends (due to the finite character of the system) and decrease of total available filler phase surface area, at a given volume fraction, with increase of the particle size.

Monsalve-Bravo and Bhatia [153] later extended their EMA model to accommodate isotherm nonlinearity using a self-consistent approach to calculate the filler phase permeability. In this work, they compared their model predictions to experimental permeation data of Vu et al. [46], who fabricated CMS-based MMMs for separation of $C{O}_2/C{H}_4$ and $O_2/N_2$ using Matrimid 5218 and Ultem 1000 as matrices. For the CMS/Ultem 1000 system, Monsalve-Bravo and Bhatia [153] showed that their modified effective medium theory (EMT) better predicted the permeability of all gases ($O_2$, $N_2$, $C{O}_2$, and $C{H}_4$) at different MMM feeding pressures and concluded that isotherm nonlinearity in the filler phase has a strong effect on the MMM permeability, with this effect becoming less significant with decrease of the particle size in the membrane system.

Figure 10 depicts a comparison of the experimental permeabilities reported in [46] for $O_2/N_2$ and the MB-B, Pal, Bruggeman, Lewis-Nielsen, Chiew-Glandt, and Maxwell models, with $r_o/\ell =0.040$ in MB-B model and ${\varphi }_m=0.645$ in the Pal and Lewis-Nielsen models. Here, percentage deviations between the experimental permeabilities and EMA models are shown in the legend of each plot. In both Figure 10a,b, the MB-B models is in good agreement with the experimental permeabilities, having the smallest deviation amongst considered EMA models.

While Monsalve-Bravo and Bhatia [153] studied the effect of isotherm nonlinearity for one of the membrane systems (CMS/Ultem 1000) reported by Vu et al. [46], the experimental permeabilities of O₂, N₂, CO₂, and CH₄ in the other system (CMS/Matrimid 5218) have been fitted to several non-ideal models [71,120,166]. This is because Vu et al. [46] hypothesized that polymer rigidification at the particle surface was the cause of the decreased permeabilities of all gases in the CMS/Matrimid 5218 MMMs. Among these works, Vu et al. [120] fitted the experimental permeability of $C{O}_2$ and $C{H}_4$ of [46] to the pseudo-two phase Maxwell model (cf. Section 3.2.1). They estimated the thickness of the rigidified interfacial layer as ${\ell }_i=0.075\mathsf{\mu }\mathrm{m}$ by assuming filler particle diameter equal to $d_o=1\mathsf{\mu }\mathrm{m}$ and interfacial permeability to $P_i=P_c/\sigma$, with $\sigma =3$. In this work, different values of $\sigma$ ranging between $1$ and $4$ were also used to adjust the interfacial thickness at fixed particle size $(d_o=1\mathsf{\mu }\mathrm{m})$, which led to the conclusion that increase of $\sigma$ yields a decrease in the interfacial thickness. A tendency later corroborated by Moore et al. [164] for zeolite-based MMMs.

Figure 11a depicts a comparison of the experimental permeabilities for from [46] and that predicted by Felske, pseudo-two-phase Maxwell and pseudo-two-phase Bruggeman models (cf. Section 3.2). Similar to Vu et al. [120], we here assumed $\sigma =3$ and $d_o=1\mathsf{\mu }\mathrm{m}$ to calculate ${\ell }_i$ in all models. Further, ${\ell }_i$ predictions and percentage deviation of each model from the experimental permeation data are shown in the legend of Figure 11a. Here, non-ideal models (overlapped) are in fair agreement with the experimental permeabilities, having percentage deviations of about $20\%$.

Alternatively, Figure 11b depicts a comparison the predicted interfacial thickness with increase of the chain immobilization factor $(\sigma )$ while assuming the filler phase particle size $d_o=1\mathsf{\mu }\mathrm{m}$. Here, an increase in the chain immobilization factor $(\sigma )$, decreases the thickness of the interfacial layer $({\ell }_i)$, similar to qualitative findings of Vu et al. [120] based on the pseudo-two phase Maxwell model. Further, the Felske and pseudo-two phase Maxwell model provide similar predictions of ${\ell }_i$ while the pseudo-two phase Bruggeman model predicts larger ${\ell }_i$ values. This behavior is expected because for a given MMM system the original Bruggeman model predicts higher permeabilities than the Maxwell model (cf. Figure 7 and Figure 10), and thus this effect is transferred to the interfacial thickness at fixed $P_f$, $P_i$, and $P_c$.

In practice, the chain immobilization factor is commonly assumed in the range $3\le \sigma \le 4$ [120,161,163,164]. This assumption often leads to interfacial layer thicknesses between the $0.05\mathsf{\mu }\mathrm{m}\le {\ell }_i\le 0.1\mathsf{\mu }\mathrm{m}$ [120], as depicted in Figure 11b. However, Dutta and Bhatia [180], via equilibrium molecular dynamics simulations, recently found that thickness of the rigidified layer to be ${\ell }_i\approx 0.001\mathsf{\mu }\mathrm{m}$ for polyimide near MFI (mordenite framework inverted) zeolite. While we recognize that thickness of the interfacial layer may vary between MMM systems, Dutta and Bhatia’s findings [180] suggests that interfacial permeability$(P_i)$ is much lower than that in the bulk polymer $(P_c)$, also evident in Figure 11b. Thus, further theoretical developments are needed to accurately assess the permeant properties in the interface, as current models rely on empirically fitting these interfacial properties.

Finally, while RMA is much less popular for describing permeation in MMMs. With the introduction of nano-flake [181,182,183] and nanotube [184,185,186] based MMMs, models such as those of Cussler [51] and KJN [48] have become more popular as reference to compare with experimental [76,114,187] or simulation-based permeation data [171,172]. Amongst these works, Wang et al. [171,172] performed simulations of 3D MMMs comprising ideal platelet and tubular particles, and compared simulation predictions to the Cussler and KJN models (cf. Section 3.1). For platelet fillers, they found the Cussler model to be in good agreement with their simulations only when isotropic diffusion was considered in the simulations [171]. For tubular particles, they found the KJN model to under-predict the simulation-based permeabilities [172]. While these simulation studies considered concentration-independent phase permeabilities (P_f and P_c), difference amongst simulation results and RMA models suggest that particle shape has a strong effect on the MMM permeability. Thus, new theoretical developments are required that overcome the limitation of current EMA/RMA models on embedding filler-related effects within empirical fitting parameters, such as in the Ebneyamini [68], Lewis-Nielsen [73], Higuchi [88], and Pal [65] models.

5. Conclusions

Existing models for transport through mixed-matrix membranes (MMMs) are adaptations of well-established theories for the thermal/electric conductivity of heterogeneous media, grounded either in the resistance model approach (RMA) or the effective medium approach (EMA). In these approaches, the MMM permeability is based on both volume fractions and permeabilities of the MMM constituent phases while assuming uniformity of the field across the MMM, and negligible filler phase particle sizes compared to membrane thickness. These considerations lead to concentration-independent permeabilities across the MMM, which mask effects of isotherm nonlinearity and finite filler size in the effective permeability calculation. Furthermore, it has been recently shown [153] that deviations to the Henry’s law, common under usual operating conditions of MMMs, lead to nonlinear concentration profiles and non-uniform permeabilities across the MMM. Thus, EMA/RMA models need to progress beyond these limitations to appropriately characterize the transport of species across the membrane.

Among existing permeation models, that of Maxwell and its adaptations (cf. Table 2), corresponding to the low-density limit of non-interacting dispersed spheres, have remain popular as reference in practical applications. However, this model has also been commonly shown to under-predict experimental MMM permeation data [45,51,81,145]. For concentrated composites, a number of analytical expressions have been proposed, including those corresponding to the Chiew-Glandt, Lewis-Nielsen, Higuchi (cf. Table 2), Bruggeman, and Pal models (cf. Table 3). Here, the Chiew-Glandt result is the exact second order solution of the transport problem through a dispersion. This model has been shown to fit simulation-based and experimental permeation data when filler particle size is negligible (cf. Figure 7 and Figure 8). Besides, an extended version of the Chiew-Glandt model (MB-B model in Table 2), accommodating both particle size and isotherm nonlinearity, has been shown to describe well the permeability for various gases in CMS-Ultem 1000 MMMs (cf. Figure 10). Alternatively, the Lewis-Nielsen, Higuchi, Bruggeman, and Pal models correspond to empirical modifications of the Maxwell or Bruggeman models to integrate particle interaction or packing-related effects; however, assuming negligible particle size in the MMM. Thus, while these models have been shown to match experimental permeation data [81], their associated empirical parameters are always fitted against the experimental MMM permeabilities, which mask effects of isotherm nonlinearity and finite filler size. Therefore, new theoretical models need to advance beyond this limitation and distinguish between effects of different mechanism to be able to achieve breakthroughs for the optimization of key separation processes.

Much effort has been devoted to mathematically represent effects of defective interfacial polymer-particle morphologies, in the form of a rigidified polymer [77,86], interfacial voids [163] or a blocked [161,179] layer around the filler particle surface, as well as particle shape (e.g., cubic, spheroidal, tubular) [69,72]. Nevertheless, while incorporating such morphological effects, existing models inherit the uniform field assumption of early RMA/EMA models, and thus disregard effects of isotherm nonlinearity and finite system size. While this drawback is addressed via simulation-based rigorous modeling of MMMs [152,171,172], EMA/RMA models need to progress beyond this limitation to appropriately characterize the transport of species across the membrane and its dependence on inherent matrix-filler properties. Furthermore, when non-ideal MMM morphologies are considered, the fitting of thickness of the interfacial layer based on an assumed permeability in this layer using MMM permeability data and EMA/RMA models is misleading (cf. Figure 11), and advances to independently characterize the interfacial layer thickness and its permeability are required.

Finally, while modeling of MMM through the RMA have received considerably less attention than EMA, with the introduction novel filler nanomaterials (e.g., carbon nanotubes and graphene nanosheets), models such that Te Hennepe, Cussler, and KJN have gained some popularity to be used as reference upon comparison with experimental and/or simulation-based permeation data. These models have been shown to successfully represent simulated MMM with flake-like and tubular particles operating at low pressures [171,172]. However, because they share the uniform concentration field assumption with EMA models, their applicability is limited to the Henry’s law region. With the upcoming advances on 3D printing for synthesis of membrane materials [6], these effects of filler particle shape and packing distribution together with isotherm nonlinearity on the MMM permeability will need to be appropriately incorporated in the new generation of models.