Gas Permeation Model of Mixed-Matrix Membranes with Embedded Impermeable Cuboid Nanoparticles

In the packaging industry, the barrier property of packaging materials is of paramount importance. The enhancement of barrier properties of materials can be achieved by adding impermeable nanoparticles into thin polymeric films, known as mixed-matrix membranes (MMMs). Three-dimensional numerical simulations were performed to study the barrier property of these MMMs and to estimate the effective membrane gas permeability. Results show that horizontally-aligned thin cuboid nanoparticles offer far superior barrier properties than spherical nanoparticles for an identical solid volume fraction. Maxwell’s model predicts very well the relative permeability of spherical and cubic nanoparticles over a wide range of the solid volume fraction. However, Maxwell’s model shows an increasingly poor prediction of the relative permeability of MMM as the aspect ratio of cuboid nanoparticles tends to zero or infinity. An artificial neural network (ANN) model was developed successfully to predict the relative permeability of MMMs as a function of the relative thickness and the relative projected area of the embedded nanoparticles. However, since an ANN model does not provide an explicit form of the relation of the relative permeability with the physical characteristics of the MMM, a new model based on multivariable regression analysis is introduced to represent the relative permeability in a MMM with impermeable cuboid nanoparticles. The new model possesses a simple explicit form and can predict, very well, the relative permeability over an extensive range of the solid volume fraction and aspect ratio, compared with many existing models.


Mixed-Matrix Membranes (MMMs) as Barrier Materials
Barrier thin films are widely used in food and beverage packaging, coating, drug release, fuel cells, and batteries [1][2][3]. Many of these films are made of polymer materials, of which the barrier properties are constrained due to the limitation in their thicknesses. By adding impermeable inorganic fillers into the polymeric matrix of the thin films, the permeation of water and gas molecules in barrier materials is further prevented or delayed. The embedding of a homogenously-dispersed inorganic filler into the continuous phase of a polymeric matrix is typically referred to as a mixed-matrix membrane (MMM). The shape, volume fraction, and dimensions of the nanoparticles within the polymer are essential factors that determine the effective gas permeability of an MMM. Wolf et al. [3] performed an exhaustive review of the effects of the shape of the fillers on the barrier properties of polymer/impermeable nanocomposite membranes. Based on 1000 published data, their results showed a very wide variation of the relative permeability of MMMs. They nevertheless confirmed the superiority of layered nanoparticles to decrease the membrane relative permeability compared to isodimensional and elongated nanoparticles. Many researchers [4][5][6] also believe that fillers have high aspect ratios create a more tortuous path for diffusion and therefore improves the barrier properties by orders of magnitude even for small volume fractions of the fillers [6]. For this reason, flake-like or plate-like fillers such as mica and clay minerals (hectorite, saponite, and montomorillonite) are studied intensively due to their high aspect ratios [6][7][8]. The aspect ratio is commonly defined as the ratio of the mean diameter of a circle of the same area of the filler to the mean thickness of the filler. There is a clear need to better understand the effect of the shape of nanoparticles on the relative permeability of MMMs and to have reliable theoretical predictions.
A number of researchers have studied and proposed models to characterize the permeation for both ideal and non-ideal morphology of MMMs [9]. Ideal membrane morphology refers to a two-phase membrane system with good adherence between the continuous phase and the discontinuous phase. The majority of theoretical models to predict the effective permeability of MMMs were developed based on ideal membrane morphology, such as the models proposed by Maxwell [10], Bruggeman [11], Lewis-Nielsen [12,13], Pal [14], Cussler [15,16], and Bharadwajl [17]. Non-ideal morphology refers to MMMs for which a non-ideal polymer-particle interface exists due to poor polymer−particle adhesion, polymer-packing disruption near the dispersed particles, and repulsive forces between the two phases. MMMs are sometimes characterized using a three-phase system, having the interface voids as the third phase, as proposed in a few modified models [18]. However, in reality, the interfacial defects are the result of many factors and vary from experiment to experiment such that the accurate determination of the exact interface morphology is a challenge. In this study, we assume that the dispersed phase is nonporous (impermeable) and ideal MMM morphology exists with no defects and no distortion at the filler−polymer interface. It is also assumed that, for a given MMM, the fillers have identical shapes and sizes, and they are uniformly distributed within the polymeric membrane. The case where interfacial defects exist is briefly discussed at the end of this study.

Main Prediction Models of the Relative Permeability of MMMs
Many models, both analytical and numerical, for predicting the relative permeability of MMMs have been proposed over the years. Six of these analytical predictive models, considered most pertinent in this investigation, are presented in Table 1. These analytical models were proposed by Maxwell [10], Bruggeman [11], Lewis-Nielsen [12,13], Pal [14], Cussler [15,16], and Bharadwajl [17]. These models have been developed to describe the permeation of a species through MMMs and to estimate their relative permeability (P r ): where P eff is the effective permeability of the MMM and P c is the permeability of the continuous phase, i.e., that of a neat polymeric membrane. In all these models, P r depends on the filler volume fraction φ, the permeability of the continuous phase P c , and the permeability of the dispersed phase P d . For a MMM with impermeable particles, P d = 0, the mathematical expressions of Table 1 are significantly reduced. The first model used to predict the permeation properties of MMM was the Maxwell's model, a model initially proposed to estimate the dielectric properties of polymer composites. Among all models listed, the Maxwell's model is the most commonly used. Besides the permeability coefficients of the continuous phase (the polymer) and the dispersed phase (the fillers), the Maxwell's model uses only the volume fraction φ as a model parameter, regardless of the particle shape, size distribution, and particle dispersion. The model is only applicable to a dilute two-phase MMM system containing spherical particles with φ < 0.2. Bruggeman [11], also attempting to predict the dielectric properties of composite materials with randomly dispersed particles, modified the Maxwell's model to predict the relative permeability over a larger range of the volume fraction. Lewis and Nielsen [12,13] studied the mechanical properties of composite materials by experimentally exploring the relationship between the relative elastic modulus of composite materials and the volume fraction of its spherical fillers. The equation introduced an additional function ψ taking into account the maximum packing fraction φ m . The model was often applied to predicting the P r of MMMs. Pal's model was originally developed to determine the effective thermal conductivity of composite materials using the differential effective medium approach. This model also introduced an additional parameter, the maximum packing fraction φ m . The Bruggeman's model is a special case of Pal's model when φ m = 1.0. Most models in Table 1 idealized the filler geometry and investigated only the effect of the filler volume fraction φ. To understand the impact of the filler's geometry, models such as the one proposed by Cussler's group relate the relative permeability to an aspect ratio and the volume fraction φ [15,16]. They extended the Maxwell's model by studying regular and random arrays of high aspect ratio impermeable particles, such as flakes and lamellae. They finally verified their model by performing experiments with the measurement of electrical resistance of salt solutions. Another well-known model is the two-dimensional model proposed by Bharadwajl [17]. Bharadwaj modified Nielsen's model and performed a theoretical study on the effects of filler sheet length, concentration, orientation, and degree of delamination on the relative permeability of MMMs. He concluded that dispersing a long sheet of inorganic filler in a polymer matrix was particularly beneficial for barrier properties. Table 1. Predictive models for the relative permeability (P r ) of a migrating species in a mixed-matrix membrane (MMM) with nanoparticles [19].
Bharadwaj (BDW) [17] , µ = 0 for randomly dispersed fillers, µ = 1 for fillers perfectly aligned perpendicular to the gas flux (12) As previously stated, based on the Wolf et al. [3] review, layered impermeable nanoparticles are more effective in reducing the permeability of MMMs than their isodimensional and elongated counterparts. Yet most of the available analytical models, including the majority of those listed in Table 1, are based on spherical nanoparticles. Wolf et al. [3] clearly demonstrated that the relative permeability of MMMs with impermeable fillers as a function of the volume fraction often behaves differently than expected, which gives rise to a series of hypotheses to explain the different trends observed. To be able properly determine the reasons for non-ideal behavior, it is essential to compare the relative permeability obtained experimentally with the one expected for an ideal MMM, thereby serving as a reference relative permeability. Therefore, the objective of this paper is two-fold. First, we would like to show that the models listed in Table 1 cannot be used to accurately predict P r of MMMs with dispersed layered nanoparticles over a wide range of volume fractions. Secondly, to develop a simple, yet general, mathematical model to predict P r of MMMs for different geometries and dimensions of the filler particles in a 3D setting. Considering that a general and easily identifiable analytical model, which covers all geometries of nanoparticles, may or may not exist, as an intermediate step between the numerical solution and an analytical model to predict P r , we also demonstrate the applicability of an artificial neural network model for the prediction of P r .

Gas Transport in an MMM
A schematic diagram of a MMM is shown in Figure 1. It is assumed that the nanoparticles are uniformly distributed throughout the whole membrane, such that it is possible to define an elementary cubic unit comprised of a single nanoparticle located at the center of the polymer cube. The membrane, therefore, consists of multiple identical elementary units juxtaposed to form the complete membrane. It can be easily shown that the permeability of a permeating species is identical for both a single elementary unit and the whole membrane [20,21]. Consequently, to estimate the permeability of the entire MMM membrane, it is sufficient to solve Fick's second law of diffusion numerically for a single elementary unit. In this investigation, the polymer phase was assumed to be isotropic, and the nanoparticle was located at the center of an elementary unit was assumed impermeable to the permeating gas. Two different types of impermeable nanoparticles often suggested as fillers in barrier films were used in this investigation. These were spherical and cuboid nanoparticles. An example of the former is TiO 2 nanoparticles, and an example of the latter is montmorillonite (MMT) clay nanoparticles. The dimensions of each cubic elementary unit were L x , L y , and L z (Figure 1), whereas the spherical nanoparticle was defined with its diameter d p , and the cuboid nanoparticle dimensions were x p , y p , and z p . Figure 1. Schematic diagram of a MMM and one elementary unit with a nanoparticle located at its center; y is the direction of the gas permeation. The dimensions of the cuboid nanoparticle are x p , y p , and z p , and the spherical particle is d p . The dimensions of the polymer elementary unit are L x , L y , and L z .
Gas permeation across a membrane follows three steps: sorption, diffusion, and desorption. The sorption and desorption processes follow Henry's law. For a membrane of constant solubility S, the concentration of the migrating species at each surface of the membrane is proportional to the surrounding pressure as described in Equation (13). It was assumed that the concentrations of gas molecules in the gas phase and within the membrane at the two gas-membrane interfaces were in instantaneous equilibrium. The diffusion process follows the Fickian's first and second laws, as expressed in Equations (14) and (15). In a neat polymeric membrane, the diffusion occurs solely in the direction (y) perpendicular to the membrane. When permeable or impermeable particles were embedded within the polymeric matrix of the membrane, the diffusion occurred in the x, y, and z directions in Cartesian coordinates.
where p is the gas pressure adjacent to a membrane surface, and C is the concentration of the migrating species within the MMM, → J is the time-varying local permeation flux (a vector), t is the time, and x, y, and z are the Cartesian coordinates to account for the three-dimensional diffusion within the membrane where y is the main direction of permeation. In this investigation, it was assumed that the membrane (continuous phase) had a constant and isotropic diffusivity D whereas the nanoparticles (dispersed phase) was impermeable and had a nil diffusivity. It was further assumed that the diffusivity was independent of concentration. In this investigation, we were, therefore, considering an ideal MMM with impermeable nanoparticles.

Methodology
The three-dimensional form of the Fick's second law of diffusion (Equation (15)) can be discretized using finite differences. Equation (16) represents the explicit discretization of Equation (15) for an interior mesh point (i, j, k) that allows calculating the concentration at the next time step knowing the current concentration at a given mesh point within the membrane. The change of concentration at mesh point (i, j, k), as a function of time, depends on the current concentration at point (i, j, k) and the concentrations at the six neighboring mesh points, as illustrated in Figure 2.
It is important to reiterate that it is only necessary to solve for a single elementary unit to estimate the effective permeability of the MMM. The elementary unit consists of a polymer cuboid (L x by L y by L z ) with a single nanoparticle at its center. To solve numerically for the temporal variation of the concentration at all points within the membrane, it was necessary to define the initial and boundary conditions. Concerning the initial condition, i.e., before the concentration step-change on one side of the membrane at time t = 0, a nil concentration was assumed throughout the membrane (Equation (17)). Concerning the boundary conditions, twelve relations were required to define the problem: Six boundary conditions at the periphery of the polymeric elementary unit and six boundary conditions at the polymer-solid interfaces. Equation (18) provided the boundary conditions on both sides of the membrane (y-axis was the permeation direction). It was assumed that at the onset of permeation, a step-change in the gas pressure was applied to the upstream side of the membrane (y = 0) whereas the gas pressure in the downstream side of the membrane (y = L y ) was kept under perfect vacuum. The resulting concentrations in the membrane at both surfaces were simply the product of the neighboring pressure and the solubility S. Because all elementary units were identical, symmetry conditions prevailed at the other four faces of the polymer parallelepiped as expressed in Equation (19) for BC [3][4][5][6] . Since the nanoparticle in the center of the elementary unit was impermeable, the mass flux at each of the six faces of the nanoparticle, assuming a cuboid nanoparticle, was zero (see BC [7][8][9][10][11][12] in Equation (20)). The latter boundary conditions imply that the concentration of the migrating species inside the nanoparticle was zero.
IC : The upstream and downstream fluxes were obtained from the concentration gradients at the upstream and downstream gas-membrane interfaces using Equations (21) and (22). Convergence to the steady permeation rate was assumed to be reached when the fluxes at the upstream and downstream surfaces were the same, within 0.001%. The effective permeability (P eff ) was calculated from the steady flux, the thickness L y of the elementary unit and the feed pressure difference between the two sides of the elementary unit of the membrane (∆p), as expressed in Equation (23).
In Equations (21) and (22), J y=0 is the upstream permeation flux and J y=Ly is the downstream permeation flux.

Numerical Simulation Results
This section presents the results for the series of 3D simulated permeation experiments of gas molecules across a polymeric elementary unit with an impermeable nanoparticle at its center. The impermeable nanoparticle at the center of the elementary unit was either a sphere or a cuboid, representing the shape of the TiO 2 and the MMT nanoparticles, respectively. For a spherical nanoparticle, the diameter was varied whereas, for a cuboid nanoparticle, the relative thickness, y p /L y , and the aspect ratio q were varied in order to assess the effect of the size and the shape of nanoparticles on the relative permeability (P r ) of MMMs. Similar to the Cussler's definition but in a 3D setting, the aspect ratio (q) is defined in a dimensionless form using Equation (24): In the current investigation, we assumed that the x-z plane of a cuboid nanoparticle was parallel to the surface of the membrane, i.e., q was the square root of the projected area divided by the thickness of the nanoparticle. Figure 3 presents the relative permeability of MMMs obtained from the simulated permeation experiments with spherical nanoparticles for a filler volume fraction φ ranging from 0 to 0.52 and with cuboid nanoparticles for a φ ranging from 0 to 0.77. Since q has a significant impact on P r , the results for the cuboid nanoparticles were grouped over a series of narrow ranges of the aspect ratio: from 0.51-0.56 to 7.00-7.29. Many more numerical results were obtained than those presented in Figure 3, and they will all be used in developing a model to predict P r . When q is large, the cuboid nanoparticle is a thin flat sheet with a large projected area. The migrating species must circumvent to reach the permeate side of the membrane, thereby reducing the permeability of the membrane. On the other hand, when q is small, the cuboid has a small projected area with a large thickness (y p ).
As expected, the relative permeability decreases with an increase in the volume fraction occupied by the nanoparticles within each q. It is evident that the aspect ratio has a major impact on the effective permeability of the membrane. Results of Figure 3 show that P r as a function of φ for MMMs with an impermeable sphere and an impermeable cube (having an aspect ratio of unity) was nearly identical. For impermeable cuboid nanoparticles with q smaller than unity, the relative permeability is greater than the relative permeability of MMMs with spherical and cubic nanoparticles. In contrast, when q is larger than unity, the situation is exactly the opposite. For a membrane with cuboid nanoparticles with the same φ, the greater the value of q, the smaller the value of P r . For all results presented in this investigation, the permeability of the continuous phase was assumed to equal to 5 × 10 −12 m 2 /s (5.0 × 10 −11 m 2 /s diffusivity and 0.10 solubility). However, it is important to note that the relative permeability was independent of the permeability of the continuous polymeric phase since the filler was impermeable. A different polymer permeability would only affect the time required to reach steady-state permeation and the steady-state permeation flux. With the numerous estimations of P r obtained numerically in this investigation, it was desired to verify if any correlation available in the literature could adequately predict the relative permeability of cuboids with large values of the q. Figure 4 compares the relative permeability obtained numerically with the one obtained using the six correlations of Table 1  Pal (PAL)) accurately predict P r over the whole range of φ. Only the models proposed by Cussler [15,16] and Bharadwajl [17], which were not developed for spheres, show more significant deviations. Based on the results of Figure 3, it is not surprising that the same four models were equally able to predict the relative permeability of the cuboid nanoparticle with an aspect ratio of unity (results not shown). For the prediction of P r for cuboid nanoparticles with a q value other than unity, Cussler's model predicts better than the four models for the two larger q values [5.0, 5.2] and [7.0, 7.3]; however, the deviations are still significant. The model proposed by Bharadwajl [17], which includes some geometrical parameters, can reasonably represent the data for small φ. However, this model becomes less accurate in predicting P r for higher φ. All the other models over-predicted P r of MMMs with cuboid nanoparticles when q was greater than unity and under-predicted P r when q was smaller than unity. It is evident that the available analytical models fail to predict accurately the relative permeability of MMMs with dispersed impermeable nano-cuboids with q > 1, in particular when φ becomes larger. The greater the value of q, the more significant the deviation between the simulated and predicted relative permeability is even at very small φ. For example, for q = 5, the deviation between the simulated P r and the one predicted by the best existing model (BWD) becomes evident at φ~0.02.  Table 1 as a function of the filler volume fraction (φ). The comparison is made for a polymeric elementary unit containing an impermeable nanoparticle: A sphere and cuboids with four different ranges of the aspect ratio (q) ( These results show that it is imperative to develop a model that would predict the permeation behavior of MMMs embedding impermeable nano-cuboids with q different from unity.

Artificial Neural Network Model for the Prediction of the Relative Permeability
If the purpose is to find rapidly an accurate model that could be used for predicting the dependent variable, such as the relative permeability (P r ), an artificial neural network can be used. Artificial neural networks are now commonly used for a myriad of engineering applications. The high degree of plasticity of its structure is the main reason for its ability to efficiently represent the underlying causal relationship between input and output data. In this investigation, a three-layer feedforward neural network (FFNN) was used to predict P r as a function of some input variables. Cybenko [22] showed that a three-layer FFNN was sufficient to encapsulate any input-output relationship if a sufficient number of neurons are used.
In this investigation, the FFNN consisted of an input layer, a hidden layer, and an output layer, as shown in Figure 5. The input layer contains many neurons corresponding to the number of independent variables plus the bias neuron. The input layer transfers each set of independent variables to each functional neuron of the hidden layer. Each functional neuron of the hidden layer performs a weighted summation of all process inputs and a nonlinear transformation of the weighted summation to generate the output of each neuron of the hidden layer. The outputs of the hidden layer, including the bias neuron of the hidden layer, are then sent to the output neuron. The output neuron performs the same task as the neurons of the second layer to generate the final output of the FFNN. A sigmoid function (Equation (25)) was used as the transfer function for the hidden and output neurons. In this investigation, the inputs and outputs of the FFNN of Figure 5 are, by definition, already scaled between 0 and 1, so that it was not necessary to scale them as it is usually done. Figure 5. Feedforward neural network for the prediction of the relative permeability (P r ) for impermeable cuboid nanoparticles as a function of the normalized projected area (x p z p /L x L z ) and the relative thickness (y p /L y ).
Various input data vectors were explored in order to find the simplest neural network structure to accurately predict P r . The two simplest structures of the FFNN were the ones that used only the relative dimensions of the nanoparticles within an elementary unit. In one case, the three relative lengths, x p /L x , y p /L y and z p /L z , were used. In the other case, the relative projected area (x p z p /L x L z ) and the relative thickness (y p /L y ) as shown in Figure 5, were used. To determine the neural model for the prediction of P r , 359 data points obtained by solving numerically the governing partial differential equation were divided equally into a training and a validation data set. The quasi-Newton nonlinear regression algorithm was used to adjust the weights of the FFNN that minimize the sum of squares of the training data set. At each iteration, the sum of squares of the validation data set was also evaluated, and the set of weights that minimize the sum of squares of the validation data set was retained. The coefficient of regression for the FFNN with six hidden neurons (including the bias) was 0.9998 for both neural network structures mentioned above. The parity plot based on the FFNN of Figure 5 is presented in Figure 6. The FFNN was able to predict very accurately P r of MMMs containing impermeable cuboid nanoparticles. The accuracy of the neural model was excellent over the entire range of x p z p /L x L z and y p /L y as the parity plot of Figure 6 shows. The most significant deviation was observed for thin cuboids covering nearly the entire x-z area of the elementary unit, thereby associated with low relative permeability. The latter condition was, however, extreme and will not be encountered in reality. The thinner the cuboid, the larger the deviation. It is believed that the FFNN can be used with confidence to predict P r . The mathematical model of the FFNN is given in Equations (26)- (28). The results obtained with the neural network suggest that a strong relationship exists between the relative projected area and the relative thickness. It is, therefore, hopeful of developing an analytical model with the two geometrical parameters. It differs from the traditional modelling approach mostly based on the volume fraction φ [23].

New Analytical Model for Pr of MMMs with Impermeable Cuboid Nanoparticles
Although the neural network model presented in the previous section provides excellent predictions of the relative permeability of cuboid nanoparticles for a wide range of geometrical parameters, it does not provide a clear insight into the impact of the input variables on P r as the analytical models listed in Table 1 do. On the other hand, these models, except for the Bharadwajl's correlation [17], fail to accurately predict P r for cuboid nanoparticles with q different from unity. Therefore, it is desired to propose a new analytical model that will be valid for the widest possible ranges of φ and the geometrical parameters of the cuboid nanoparticles.
To clearly elucidate the relationship between the shape and the relative dimensions of nanoparticles and P r , multivariate covariance analysis was performed to assess the underlying correlation between all possible geometrical factors and P r . For this analysis, the Pearson correlation coefficient (PCC), defined in Equations (29) to (31), was used [24]. Cov where σ is the standard deviation of the variable E, Cov(E A , E B ) is the covariance of E A and E B , E A and E B are the average values of E A and E B , and PCC(E A , E B ) is the Pearson correlation coefficient between E A and E B . In this analysis, E A corresponded to one geometrical variable to be tested, and E B was P r determined numerically. The same 359 data points used for developing the neural network were used for this analysis. The results of this statistical analysis for five potential geometrical factors are presented in Table 2. The results show that x p z p /L x L z had, in agreement with the results of the previous section, the highest negative PCC with P r . The x p /L x ranked second because it was equivalent to the square root of the x p z p /L x L z . The volume fraction φ also had a significant correlation factor with P r . The aspect ratio and the relative thickness also had some impact on P r but to a lesser degree. It is important to note that the five selected geometrical variables are not all mutually independent. To better comprehend the effect of these geometrical variables on P r , one needs to examine the permeation process. Gas molecules entering a membrane will diffuse freely through the polymer matrix in the absence of impermeable nanoparticles. In that case, based on an elementary unit, the entire surface area of the polymer L x L z is available for diffusion. When a nanoparticle is introduced into an elementary unit, the gas molecules have to adopt a more tortuous path to diffuse around the impermeable nanoparticle. Figure 7 presents a plot of the isoconcentration lines within the polymer where the concentration within the cuboid particle is zero. Since the diffusion streamlines of gas molecules in the presence of an impermeable cuboid nanoparticle run perpendicular to the isoconcentration lines, one can easily imagine the diffusion path of these gas molecules. Indeed, the diffusion streamlines are perpendicular to the x-z plane of the elementary unit at the gas-membrane interface and deviate progressively as they approach the impermeable nanoparticle. The density of the isoconcentration lines is indicative of the intensity of the local flux. The flux increases when the diffusion channel narrows down and decreases when it widens up.
Considering the dependence of the local permeation flux on the size of the permeation channel for a given y-position, a dimensionless parameter A*, defined in Equation (32), was introduced. This dimensionless parameter was simply the ratio of the projected area that is available for diffusion (L x L z − x p z p ) and the maximum surface area for diffusion (L x L z ). It is logical to postulate that an increase in A* will lead to an increase in P r of the MMM, and vice versa. Figure 8 clearly shows the strong relationship between P r of a MMM with cuboid nanoparticles and the dimensionless parameter A*.
where x p z p is the projected area of the nanoparticle, while L x L z is the total permeation area of an elementary membrane unit. In addition to the influence of A* on P r , Figure 8 illustrates the impact of the relative thickness y p /L y over the range spanning from 0.0500 to 0.9833. For a nanoparticle with a fixed y p /L y , a decrease in x p z p leads to an increase of A* and an increase in P r . For a fixed A*, P r decreases when y p /L y increases. The plot of P r versus A* at the largest y p /L y is approaching the 45 • line. In other words, P r approaches A* when y p /L y tends to 1.0. It is clear that the model to be developed must include the strong linear relationship of P r with A* and a nonlinear component to account for the effect of y p /L y .  To characterize and better understand the contribution of the nonlinear portion (NLP) of the relative permeability curve, Figure 9 is presented to highlight the breakdown of P r in Figure 8 into its linear (P r = A*) and NLP as a function of A*. This data exploration is essential in searching for the right form of the analytical equation in the development of an accurate predictive model. The contribution of the NLP that is related to y p /L y is clearly illustrated in Figure 9. The nonlinear term first increases rapidly with A* up to a maximum value before decreasing more gently to zero as A* tends to unity. The location and the magnitude of the maximum are a function of the relative thickness y p /L y . The maximum value is located in the interval of A* between 0.16 to 0.38. The contribution of the nonlinear part can be as high as 0.37 of the value of P r at y p /L y = 0.05 and A* = 0.18. Based on the insight of the above information, a new analytical model for the prediction of P r is now presented in Equations (33) to (35). The model depends strictly on two simple geometrical parameters: A* and y p /L y . Equation (33) is the sum of the linear and nonlinear parts of the estimated P r . The parameters a and b in Equation (33), given by Equations (34) and (35), respectively, are a function of only the relative thickness. It is interesting to note that φ is not used explicitly in Equation (33). On the other hand, both A* and y p /L y indirectly determine the value of φ. Figure 10 shows the comparison between the relative permeability obtained by solving Fick's second law of diffusion numerically and the one predicted with the proposed model. To assess the accuracy of the model prediction, the average prediction error was calculated using Equation (36).
where ε A is the average prediction error, n is the number of data for the average error calculation, and P r andP r are the relative permeability obtained numerically and the one calculated by the proposed model, respectively. Table 3 summarizes the average prediction error for the combination of the geometrical parameters presented in Figure 10. In this investigation, 359 numerically predicted values of the relative permeability were used to fit the model by minimizing Equation (36) using the Levenberg-Marquardt algorithm. Figure 10. Plots of the P r of MMMs embedding cuboid nanoparticles as a function of A* for five different values of y p /L y . Table 3. Average prediction errors of the P r of Figure 10 using the proposed model evaluated over five values of the y p /L y . Given its simplicity, the proposed model predicts P r accurately over a wide range of A* and y p /L y . Table 3 shows that for each group, the relative particle thickness, ε A is always below 0.03. The proposed model resorting to two simple geometrical parameters can be used with confidence to predict the permeability of MMMs with impermeable cuboid nanoparticles.

Prediction of Experimental Data with the Proposed Model
To validate the proposed model for the prediction of the relative permeability of MMMs with elongated impermeable nanoparticles, the data of a recent paper published by Zahid et al. [25] was used. This paper reports an investigation on the effect of the morphology of four types of layered graphene flakes, used as fillers, on the gas barrier properties of thermoplastic polyurethane (TPU) films. Fortunately, this paper provides the geometrical parameters of the nanofillers and the experimental permeability data for the oxygen transport across polyurethane films. The digitalized permeability data were converted into relative permeability data using Equation (1). Results in their paper were presented in terms of the mass fraction of the fillers such that, using the bulk density of the four types of layered graphene flake powders provided in their paper and a density of 1.08 g/cm 3 [26] for the specific thermoplastic polyurethane used, the filler volume fraction was calculated. A void fraction of the nanofiller powder flakes of 0.38 [27] was used in the conversion from mass fraction to volume fraction. The comparison of the experimental relative permeability P r obtained by Zahid et al. [25] and the P r predicted by the proposed model is presented in Figure 11. The geometrical parameters A * and y p /L y depend on the spacing between the particles for fixed particle dimensions; however, this information was not provided by Zahid et al. [25]. Therefore, for each set of experimental results, three values of A * were used in the model to obtain reasonable estimations. Correspondingly, y p /L y was decreased for an increase of A * to maintain constant the filler volume fraction φ. Results show that the experimental data fall among the predicted P r curves for the three values of A * . The proposed model assumed ideal MMMs whereas, in an actual membrane, the dimensions, orientation, and spatial distribution may show significant variability. Despite these potential discrepancies, the comparison of Figure 11 shows that the model can adequately represent the relative permeability of the four types of MMMs presented in the paper of Zahid et al. [25]. On the other hand, there are many studies in the literature where the relative permeability increases with the impermeable filler volume fraction and many studies where the trend varies significantly with the volume fraction [3]. Many reasons were postulated to explain these discrepancies. We strongly believe that the comparison of the relative permeability of the experimental data with the data of ideal MMMs, calculated with the proposed model, could help in finding the sources of the non-ideality. The ideal MMM serves as a benchmark or a point of reference to explain these non-idealities.

Discussion and Conclusions
In this paper, the permeation process of gas molecules through MMMs embedding impermeable spherical and cuboid nanoparticles of different aspect ratios was simulated using a three-dimensional finite-difference solution of Fick's second law of diffusion. These numerical simulations allowed assessing the effect of the shape and size of nanoparticles on the barrier properties of the resulting MMMs. At the same time, these simulations, which can be referred to as numerical experiments, allowed gaging the applicability of different analytical models to predict the relative permeability of these membranes.
Simulation results showed that, for the same volume fraction, the relative permeability of MMMs with cuboid nanoparticles with an aspect ratio larger than unity was lower than the relative permeability of membranes with spherical nanoparticles. The relative permeability decreases with an increase in the aspect ratio of cuboid nanoparticles. The Maxwell's and similar models predict the relative permeability of MMMs with spherical nanoparticles very well, within a wide (small and medium) particle size range. However, the Maxwell's model and other analytical models developed for spherical nanoparticles, showed a weak prediction of the relative permeability of MMMs with cuboid nanoparticles. Using an artificial neural network, we developed a feedforward neural network model, which very accurately predicts the relative permeability of MMMs with both spherical and cuboid nanoparticles for a wide range of sizes and the aspect ratios. However, although very accurate, the developed feedforward neural network model lacks analytical models' simplicity and explicability. To overcome this limitation of neural network-based models, we used the Pearson correlation coefficient and performed a multivariate covariance analysis to elucidate the relationship between the shape and the relative dimensions of nanoparticles and the relative permeability of MMMs. This led to an analytical model, which involves only two geometrical parameters, the ratio of the projected area available for diffusion and the maximum surface area for diffusion, which depends on the aspect ratio of the nanoparticle, and the relative thickness of the nanoparticle. Despite its simplicity, the new model accurately predicts the relative permeability of the MMMs with cuboid nanoparticles for a wide range of sizes and aspect ratios. To our best knowledge, the model developed in this study is the first one in the literature, which is applicable for MMMs with both spherical and layered nanoparticles.
The paper assumed an ideal case where the MMMs had ideal interfacial compatibility, no particle agglomeration, and particles were perfectly aligned. However, in reality, voids could exist due to poor adhesion between fillers and the polymer matrix; agglomeration of particles may create nanogaps between the particles and the polymer chains [9]. These non-ideal morphologies may increase or decrease the Pr such that the effects of non-ideality need to be further investigated even though it is very challenging to represent non-ideality in practice. One needs to idealize non-ideal morphology for modelling. This investigation on ideal MMMs is the first, but very important, step to develop an accurate model for MMMs with cuboid nanoparticles. Indeed, as noted by Vinh-Thang and Kaliaguine [9], good predictions are only achieved for MMMs in which morphology could be represented as a combination of voids with an ideal filler/polymer interface. In addition, as suggested by Hamid and Jeong [28,29], the embedding of nanosized particles within the polymer matrix will favor a uniform distribution of the nanoparticles and better wetting property between the filler phase and the polymer phase [30] compared to the distribution of micron-sized particles, which are often used in the fabrication of MMMs. The smaller size particles will allow the fabrication of thin membranes and will enhance filler-polymer interface contact.