Gas Barrier Properties of Multilayer Polymer–Clay Nanocomposite Films: A Multiscale Simulation Approach

: The paper discusses the development of a multiscale computational model for predicting the permeability of multilayer protective ﬁlms consisting of multiple polymeric and hybrid layers containing clay minerals as ﬁllers. The presented approach combines three levels of computation: continuous, full atomic, and quantitative structure–property correlations (QSPR). Oxygen and water are chosen as penetrant molecules. The main predictions are made using the continuum model, which takes into account the real scales of ﬁlms and nanoparticles. It is shown that reliable predictions of the permeability coefﬁcients can be obtained for oxygen molecules, which is not always possible for water. The latter requires the reﬁnement of existing QSPR methods and interatomic interaction potentials for the atomistic level of calculations. Nevertheless, we show that the maximum effect on permeability reduction from the addition of clay ﬁllers to the hybrid layer can be achieved by using nanoparticles with large aspect ratios and a high degree of orientational order. In addition, the use of the hybrid layer should be combined with the use of polymer layers with minimal oxygen and water permeability. The constructed model can be used to improve the properties of protective coatings for food and drug storage and to regulate the gas permeability of polymeric materials.


Introduction
Gas permeability is one of the most important properties in assessing the operational performance of packaging and protective materials designed to ensure long-term storage of food and pharmaceuticals [1][2][3][4][5][6][7][8][9][10][11][12][13][14].The choice of polymers for these applications is explained by their high manufacturability and design flexibility, while products based on them have a lower weight that cannot be offered by glass and metal.
One of the main problems associated with the preservation of foods and drugs is their relatively rapid degradation when in contact with oxygen-containing molecules in the environment [5,8,10,[15][16][17].This is due to the fact that all organic materials (oils, lipids, sugars, and proteins) are susceptible to oxidation and many of them can bind to water.In turn, medicines contain biologically active molecules (enzymes, hormones, vitamins, antibiotics, etc.) that have many chemically active functional groups (-S-, -OH, -C(O)OH, -C(O)O-, -NH-, -NH 2 , -NO 2 , etc.).Therefore, the penetration of even a small amount of oxygen or water molecules into the packaging can significantly reduce the consumer properties of the products and gradually lead to their complete deterioration.
Recently, flexible multilayer films have been increasingly used in the production of packaging for fresh vegetables and meat, as well as bags and containers for the storage of young wine, drugs, medicines, etc. [4,7,9,11,13,18,19].Commercial multilayer films typically consist of three to nine layers of different polymers with or without an aluminum interlayer.The outer layers consist of polymers with suitable mechanical properties, such as polyethylene terephthalate (PET), polypropylene (PP), and high-density polyethylene.The inner layer is the sealing layer, which provides a hermetic seal to protect the product.Copolymers consisting of ethylene are often used as sealants due to their low melting points [20,21].These layers are typically bonded together with various types of adhesives.To ensure the long-term preservation of some drugs, the used packaging materials must have extremely low water vapor and oxygen permeabilities of less than 10 −1 g/(m 2 •day) and 10 −1 cm 3 /(m 2 •day), respectively.Furthermore, flexible coatings have to withstand repeated bending loads.Typically, thin aluminum interlayers are introduced to achieve the required ultra-low permeability [5,22].
Today, some biopolymer films are also being adapted with synthetic polymers to meet the current trends to reduce polymer waste and facilitate recycling [13][14][15]23,24].At the same time, the use of aluminum layers increases the environmental impact and reduces the recyclability of such protective coatings.
To fulfill the requirements of environmental safety on the one hand and to ensure the barrier properties required on the other, it is proposed to use hybrid multilayer films that include pure polymers, clay-based polymer nanocomposites (Figure S1a), and continuous inorganic layers (see Figure S1b; note that "S" indicates any reference to the Supporting Materials file) [3,[25][26][27][28].In the first case, the permeability of the hybrid layer depends on the filler properties, its concentration, and the permeability of the matrix.In the second case, permeability is mainly defined by defects in the inorganic layer [29,30], such as pinholes (see Figure S1b).These can be caused by impurity particles on the substrate surface on which the barrier films are deposited.Experimental data show that the scale of these pinholes is (sub)micron, and their density for ultrabarrier films (whose permeability is approximately 10 −6 cm 3 /(m 2 •day)) should not exceed 100 pcs/cm 2 [30].In addition, defects in pure polymers and hybrid films can be induced by the crystallization of polymer molecules.This process causes the formation of grain boundaries and leads to a significant increase in the diffusion rate [3,31,32].Inorganic coatings based on nanolaminates, in which the amorphous phase is stabilized, can be used as a way to suppress crystallization [33,34].In addition, additional defects in inorganic layers can be caused by cracks when flexible coatings are bent [35,36] (see Figure S2).Therefore, hybrid protective layers based on polymer nanocomposites, as shown in Figure S1a, are more suitable for flexible barrier coatings because they are less subject to mechanical degradation.These coatings are chosen as the objects of our study.
Two types of clay-polymer nanocomposite films are possible, as shown in Figure 1: nanocomposites with a more or less random orientation of the clay fillers (Figure 1a), and systems with an enforced orientation order of the clay particles (Figure 1b).The first structure can be obtained by mixing polymer/clay in solution or melting followed by film formation.The second can be obtained, for example, by the layer-by-layer deposition method, where the substrate is sequentially treated with anionic and cationic aqueous solutions of polymers and clay particles [37], which eventually leads to the formation of the nanolaminate structure [38].Experimental studies have shown that such materials have very low permeability.In particular, less than 0.005 cm 3 /(m 2 •day) for oxygen has been demonstrated [37].
Therefore, further improvement of packaging materials requires research aimed at the simultaneous optimization of the structures of the hybrid protective layer and the protective multilayer film as a whole.In addition, there is a separate problem of selecting a suitable polymer for each layer of the film.These issues require long-term experimental studies that combine the efforts of large teams of chemists and engineers, as well as the use of various laboratory equipment.In this situation, the use of theoretical models describing mass transfer through multilayer hybrid polymeric media [6,[39][40][41][42] can significantly reduce the overall duration of the development process.However, building adequate theoretical models is a nontrivial problem, since it is necessary to consider the structure of each polymer, the inorganic fillers, the entire barrier material, and the properties of the interactions of all components with the molecules of penetrating gases.In this case, the use of computer simulation is a compromise solution, since it can work with models that combine several theoretical methods.At the same time, for laboratory use conditions, ready software solutions are in great demand which allow predicting the properties of materials in automatic mode.Typically, available software packages allow for the study of monolayer films only [43][44][45][46].Therefore, further improvement of packaging materials requires research aimed at the simultaneous optimization of the structures of the hybrid protective layer and the protective multilayer film as a whole.In addition, there is a separate problem of selecting a suitable polymer for each layer of the film.These issues require long-term experimental studies that combine the efforts of large teams of chemists and engineers, as well as the use of various laboratory equipment.In this situation, the use of theoretical models describing mass transfer through multilayer hybrid polymeric media [6,[39][40][41][42] can significantly reduce the overall duration of the development process.However, building adequate theoretical models is a nontrivial problem, since it is necessary to consider the structure of each polymer, the inorganic fillers, the entire barrier material, and the properties of the interactions of all components with the molecules of penetrating gases.In this case, the use of computer simulation is a compromise solution, since it can work with models that combine several theoretical methods.At the same time, for laboratory use conditions, ready software solutions are in great demand which allow predicting the properties of materials in automatic mode.Typically, available software packages allow for the study of monolayer films only [43][44][45][46].
The goal of the present study was to develop a computational technique for the rapid quantitative evaluation of oxygen and water vapor permeation through multilayer organic-inorganic barrier films.It is based on a multiscale modeling approach that combines several computational methods [47,48].We developed a continuum model of a multilayer coating film consisting of several polymer layers and one polymer/clay nanocomposite layer.The model is parameterized using the quantitative structure-property correlation (QSPR) methods of Bicerano [49] and Askadskii [50].For additional verification of these two methods, full-atom modeling methods based on molecular dynamics (MD) and Monte Carlo (MC) methods are used.They allow the study of the gas permeability of polymers, based only on the structure of the chemical components and a composition of the material with minimal external parametrization [51].However, MD and MC are limited to relatively small systems on the nanometer scale.This fact reduces their applicability to the modeling of multilayer polymer films and nanocomposites.Therefore, we have paid special attention to exploring the possibility of using these methods to make predictions of water and oxygen permeability when using material samples of small size.
The article is organized as follows: Section 2 describes the hierarchy of levels of the model to determine the water vapor and oxygen permeability of the multilayer protective film.We give a detailed description of each level of modeling: continuum, QSPR, and full atomic based on MD and MC methods.In Section 3, we describe the results of calculating the transport properties of selected polymers using the full atomic approach and QSPR.The goal of the present study was to develop a computational technique for the rapid quantitative evaluation of oxygen and water vapor permeation through multilayer organicinorganic barrier films.It is based on a multiscale modeling approach that combines several computational methods [47,48].We developed a continuum model of a multilayer coating film consisting of several polymer layers and one polymer/clay nanocomposite layer.The model is parameterized using the quantitative structure-property correlation (QSPR) methods of Bicerano [49] and Askadskii [50].For additional verification of these two methods, full-atom modeling methods based on molecular dynamics (MD) and Monte Carlo (MC) methods are used.They allow the study of the gas permeability of polymers, based only on the structure of the chemical components and a composition of the material with minimal external parametrization [51].However, MD and MC are limited to relatively small systems on the nanometer scale.This fact reduces their applicability to the modeling of multilayer polymer films and nanocomposites.Therefore, we have paid special attention to exploring the possibility of using these methods to make predictions of water and oxygen permeability when using material samples of small size.
The article is organized as follows: Section 2 describes the hierarchy of levels of the model to determine the water vapor and oxygen permeability of the multilayer protective film.We give a detailed description of each level of modeling: continuum, QSPR, and full atomic based on MD and MC methods.In Section 3, we describe the results of calculating the transport properties of selected polymers using the full atomic approach and QSPR.Section 4 discusses the results obtained using the continuum model.Section 5 discusses the obtained results and the problems associated with the use of the selected methods.In the final, Section 6, a summary of the work and conclusions on the results obtained are given.

Multiscale Model Description
As mentioned above, a typical flexible barrier coating consists of several polymer layers, each providing specific functionality.The outer layer protects against mechanical damage and has hydrophobic properties.The substrate layer, which contains materials such as polyethylene terephthalate (PET), provides mechanical strength.All parts must be firmly bonded together, and special adhesive layers are used for this purpose.
The fastest predictions of the properties of multilayer polymer films can be made using continuum models, where each layer is considered as a nonporous membrane or as a continuous medium.The concept of dissolution and diffusion, which describes the phenomena of gas transport in each layer, was first formulated in [52][53][54][55].According to this concept, the penetration of small molecules through a membrane consists of three stages: (1) gas adsorption on the surface of the material, (2) diffusion of the sorbed gas molecules through the volume of the l layer, and (3) desorption from the opposite surface (see Figure S1a).Thus, the transport properties of a multilayer barrier are determined by the sum of the reverse transport resistances of each layer.To use this model, it is necessary to know the properties of all the polymeric materials in the system, primarily their water vapor and oxygen permeabilities.In the case of a nanocomposite polymer layer with fillers, the permeability is determined by the distribution of filler particles in the layer and the diffusion of penetrant molecules through the polymer matrix.Therefore, we developed a continuum model that solves the diffusion problem of small gas molecules in the case of an arbitrary distribution of clay particles in a polymer matrix.
However, in some cases, data on the transport properties of penetrant molecules in a polymer layer may not be available, e.g., when new polymers are synthesized.To solve parameterization problems, we consider using the well-established quantitative structureproperty relationship (QSPR) regression models of Bicerano and Askadskii [49,50].These models allow us to make quick estimates based only on the chemical structure of the repeated polymer unit.They are also convenient for screening polymers by designing new variants of chemical structures.In addition, they are well described and have an open parameterization, which makes them portable and verifiable.
If the chemical structure of selected polymers is too different from that of the types of polymeric materials used to parameterize the QSPR models, these models may produce predictions with a large error.In such cases, atomistic methods based on molecular dynamics and Monte Carlo methods can be used to test the predictions for selected polymers.For our multilayer polymer coatings model, we used the implementation of these methods in our MULTICOMP package [48].
Thus, our model is a multiscale computational scheme based on the use of a hierarchy of three approaches: continuum, QSPR, and atomistic.Let us describe each of the modeling levels in detail.

Continuous Level
Methods to model the transport properties of multilayer hybrid films at the continuum level typically assume that permeability is determined by the product of solubility and diffusion (see Section S2).Inorganic additives to polymers are assumed to be impermeable barriers to permeating molecules that increase their diffusion path, i.e., decrease the permeability of the film.
The effect of diffusion barriers is determined by both the volume fraction of inorganic fillers and their shape and orientation.Thus, for the simple case of an impermeable layered filler orientated perpendicular to the flow, the effective diffusion coefficient decreases as (Nielsen's model [56]): where α and ϕ are the aspect ratio and volume fraction of the filler particles, and D 0 and D are the diffusion coefficients in the pure polymer matrix and nanocomposite, respectively.In Nielsen's model, permeability decreases linearly with increasing filler aspect ratio, which explains the widespread use of plate-like fillers, such as clay particles or graphene, in nanocomposite barrier films [57,58].The Nielsen approximation is accurate in the dilute regime but may be inaccurate in the semi-dilute case where the volume fraction of the filler is small ϕ << 1 but its overlap is large αϕ >> 1, as suggested by Cussler et al., who derived the following approximation for 2D geometry [59]: where σ is the aspect ratio of a slit between clay particles (ratio of the width of the slit to slit thickness).The stronger dependence of the effective diffusion coefficient on the overlap factor compared to that of Nielsen's model is related to the reduced cross-sectional area between the clay plates.More complex expressions have been obtained for certain configurations of 2D layered fillers [59,60], but a general analytical solution cannot be obtained for an arbitrary configuration.Therefore, numerical simulation methods are used as a good compromise to analyze transport in more realistic configurations of inorganic fillers.
There are two main approaches for the numerical modeling of the transport properties of multilayer hybrid coatings at the continuum level: the deterministic method, based on the solution of partial differential equations, and the stochastic method, based on the Monte Carlo technique.The use of the first method is limited by the presence of fillers with a large aspect ratio or a large difference in the characteristic sizes in the system (for example, between the hole size in the inorganic layer and the thickness of the polymer layer).This leads to high computational costs and the need to use complex adaptive grids to solve the diffusion equation.Therefore, Monte Carlo methods are often used to analyze transport in complex configurations of inorganic layers.For example, they have been used to calculate diffusion in polymers with layered fillers [61].Since this approach is well established, we use it in our continuum model to describe gas transport.
In our model, the protective coating film consists of alternating polymer and hybrid polymer/leather layers (see Figure 1).Each polymer material is characterized by its diffusion coefficient D of penetrant molecules and their solubility S, as well as the thickness of the layer.In addition, it is assumed that different filler particles do not intersect and do not stick out from the layer.Now let us discuss the implementation of the chosen model.It is based on tracking the random motion of a molecule, taking into account the presence of impenetrable barriers.In the initial state, the molecule is located in the first polymer layer near the outer surface of the film, and then the molecule makes microsteps, i.e., random displacements whose value depends on the distance d to the nearest barrier: The minimum step size d min depends on the minimum geometric size in the system (e.g., the size of the hole or the thickness of the layered filler).After displacement, the total diffusion time increases by the value which is determined by the diffusion coefficient D i of a molecule in an i-th layer.
The probability of a molecule passing from one polymer layer to another is determined by the ratio of its solubilities: where S i is the solubility of the current layer and S i+1 is the solubility of the layer into which the molecule is trying to penetrate.The model gives the diffusion time τ diff of the molecule throughout the thickness of the multilayer coating.This time is then averaged over the number of attempts for a set of molecules.By comparing the average diffusion times in a multilayer system with inorganic fillers, <τ diff >, and without them, <τ diff > free , we can estimate the change in the permeability of the coating as: where P free is the permeability of the multilayer system without fillers and inorganic layers: and P i free = S i •D i is the permeability of the i-th layer of the multilayer system without fillers.
Therefore, to predict and optimize the protective properties of the barrier coating, it is necessary to determine the parameters of the transport properties (S i and D i ) for each polymer layer.To solve this problem, we used quantitative structure-property correlations (QSPR) and atomistic modeling methods discussed in the following sections.

The Quantitative Structure-Property Correlations
For rapid evaluations of polymer permeability, we used the QSPR regression models of Bicerano and Askadskii [49,50].These models are well known and tested and allow us to calculate many properties of polymeric materials (such as gas transport, Young's modulus, etc.).Previously, we implemented Bicerano's models to predict polymer properties in our MULTICOMP software package [48].However, the Bicerano approach lacks a model corresponding to the permeability of polymers to water P(H 2 O).Therefore, we use the Bicerano model only to estimate the oxygen permeability P(O 2 ).To estimate P(H 2 O), we selected the Askadskii model [62,63].
The structural formulas of the monomers are used as input information for the Askadskii and Bicerano models (see Figure S3).To prepare input data, the structures of the monomers are decomposed into predefined atoms (Askadskii) or structural fragments (Bicerano), and then their number is counted.The water vapor permeability in the Askadskii method is determined by the following expression: where p 0 is a constant value equal to 3.002 (in Barrer units) for all polymers, ∆E is the activation energy, R is the gas constant, V Vdw is the van der Waals volume, and T is the temperature.The value ∆E is calculated as the sum of the atomic and group contributions ∆E j : For the parameterization of the atomic and group contributions, we use the values presented in Table S1.Analysis of this table leads to the conclusion that the lowest permeability value may be obtained for aliphatic polymers (with a small positive contribution to the activation energy) and for polymers with the maximum number of chlorine atoms, which gives a negative contribution to the activation energy.
The Bicerano model [49] for predicting oxygen permeability includes the following expressions: These equations contain the following parameters: cohesion energy (E coh ), molecular volume (V), number of rotating bonds (N rot ), number of non-hydrogen atoms (N), number of acyclic carbon-carbon double bonds (N C=C ), number of ester groups in the main polymer chain (N bb,ester ), number of substituents in aromatic cycles in the main chain (X 4 ), the sum of Cl and Br atoms attached to sp 3 carbon atoms (N hheq,σ ), the sum of cyanide groups attached to sp 3 carbon atoms (N cyanideeq,σ ), and number of hydroxyl hydrogen atoms and aromatic cycles with hydrogen bonds (N hb,ar ).For the compounds studied, the calculation of P(O 2 ) requires knowledge of the values of the parameters E coh , V, V VdW , N rot , N, X 4 , and N hheq,σ .The other parameters are zero.The parameters E coh , V, and V VdW are estimated using the Bicerano method (implemented by us in the MULTICOMP package [48]).N rot , N, X 4 , and N hheq,σ are estimated from 2D structural formulas (see description in Ref. [49]).We used atomistic modeling as a computational experiment to verify the results of predictions made using regression QSPR models.
It should be noted that these models are useful for analyzing polymer databases to screen for materials with optimal properties.This analysis can be extensive, using databases of known polymers, or intensive, based on virtually constructed polymer structures.We considered the second possibility (see Section S4) when selecting polymeric materials to parameterize our continual model.We wrote a special program for this purpose and took into account the resulting ambiguities when determining the uniqueness of the generated polymer structures using the InChIkey identifier (see Section S4.1) [64,65].

Full Atomic Level
One of the most widely used theoretical approaches to estimate the transport properties of polymeric materials is atomistic modeling based on a combination of molecular dynamics and Grand Canonical Monte Carlo (GCMC) [66][67][68][69][70].The implementation of this scheme is shown in Figure 2. volume (V), number of rotating bonds (Nrot), number of non-hydrogen atoms (N), number of acyclic carbon-carbon double bonds (NC=C), number of ester groups in the main polymer chain (Nbb,ester), number of substituents in aromatic cycles in the main chain (X4′), the sum of Cl and Br atoms attached to sp 3 carbon atoms (Nhheq,σ), the sum of cyanide groups attached to sp 3 carbon atoms (Ncyanideeq,σ), and number of hydroxyl hydrogen atoms and aromatic cycles with hydrogen bonds (Nhb,ar).For the compounds studied, the calculation of P(O2) requires knowledge of the values of the parameters Ecoh, V, VVdW, Nrot, N, X4′, and Nhheq,σ.The other parameters are zero.The parameters Ecoh, V, and VVdW are estimated using the Bicerano method (implemented by us in the MULTICOMP package [48]).Nrot, N, X4′, and Nhheq,σ are estimated from 2D structural formulas (see description in Ref. [49]).We used atomistic modeling as a computational experiment to verify the results of predictions made using regression QSPR models.
It should be noted that these models are useful for analyzing polymer databases to screen for materials with optimal properties.This analysis can be extensive, using databases of known polymers, or intensive, based on virtually constructed polymer structures.We considered the second possibility (see Section S4) when selecting polymeric materials to parameterize our continual model.We wrote a special program for this purpose and took into account the resulting ambiguities when determining the uniqueness of the generated polymer structures using the InChIkey identifier (see Section S4.1) [64,65].

Full Atomic Level
One of the most widely used theoretical approaches to estimate the transport properties of polymeric materials is atomistic modeling based on a combination of molecular dynamics and Grand Canonical Monte Carlo (GCMC) [66][67][68][69][70].The implementation of this scheme is shown in Figure 2.  In the first stage, samples of the polymeric materials to be studied are generated and optimized.The MD-GCMC procedure is then performed for the chosen number of N GCMCMD steps to obtain the optimized structure with penetrant molecules (in this study, for O 2 and H 2 O).
In the GCMC method, the solubility coefficient is estimated by calculating an adsorption isotherm curve at fixed temperature and volume.The number of molecules in the ensemble can vary (only the chemical potential µ is fixed).The solubility coefficient is obtained by approximating the adsorption curve with a straight line (as the pressure tends to zero) and calculating the slope of the constructed straight line N cell (p), where N cell is the calculated equilibrium number of penetrant molecules in the simulation cell for a given injection pressure p.
For the resulting matrix + penetrant model, MD simulation is performed to calculate the mean square displacement (MSD) of the penetrating molecule.MSD allows the calculation of the diffusion coefficient D [66,67] using the Einstein formula: where r i (t) is the position of the center of mass of a selected molecule at time t and N is the number of diffusing molecules.If the root-mean-square displacement curve can be approximated with a straight line with slope k, the equation given above can be simplified as follows: The resulting permeability coefficient P of the polymeric material is calculated as The described combination of MD+GCMC and MD simulations is implemented using the LAMMPS code as part of a new "Permeability" module (developed specifically for this study) and included in our MULTICOMP package [48].This module performs the two-step atomistic simulation (see Figure 2) that allows estimating the solubility coefficient S and the diffusion coefficient D for the selected structure of the solvent and the selected penetrant molecule, and uses them to calculate the permeability coefficient P using Equation (15).Unlike previous similar studies using the MD-GCMC method [66,[71][72][73][74], the atoms in the polymer matrix and the penetrant molecules are not fixed during the MD-GCMC or MD simulation runs.This allows us to describe the structural changes due to the addition of penetrant molecules to the matrix.

Results of Full Atomistic Calculations
As a first step in the implementation of the designed multilevel simulation of protective films, we selected and performed a comparative study of the transport properties of several polymers at the atomistic level.The diffusion, solubility, and permeability coefficients were calculated.Oxygen (O 2 ) and water (H 2 O) molecules were used as penetrants.The SPC/E model was used for water molecules [75].Furthermore, the permeability coefficients were predicted using QSPR methods.
To select suitable polymers for our continuous model, we generated a series of polymers and screened them using QSPR methods, showing that organochlorine materials have the best barrier properties (see Section S4).However, such polymers do not meet environmental requirements, as they release toxic gases during degradation [76].Other obtained polymers have higher vapor permeability values, which limits the possibility of improving the barrier properties of polymer films by modifying the chemical composition of the matrices.In this situation, it seems more promising to improve the barrier properties of polymer films by using hybrid layers with inorganic fillers, such as clays and standard polymer matrices.Therefore, the choice of polymeric materials was based on data on their availability and use in the composition of protective materials (see Section S4.2).As a result, we selected common plastics such as polyethylene terephthalate (PET), polyethylene (PE), and two fluorinated polymers: polyvinylidene fluoride (PVDF) and polytetrafluoroethylene (PTFE).A large amount of experimental data is also available for these materials.
Therefore, the calculation of the transport properties of PET, PE, PVDF, and PTFE serves on the one hand to parameterize our continual model, and on the other hand, the comparison with the experiment is used to verify the adequacy of the chosen methods.

Generation of Polymer Matrices
Polymer matrices were prepared using constructor modules (integrated into our previously developed MULTICOMP package [48]) and the Class II Polymer Consistent Force Field (PCFF) [77].Polymer chains were constructed using the "Polymer Chain" module.The degree of polymerization of the chains was set to 33 for PE, PTFE, PVDF, and 11 for PET.The length of the polymer chains was chosen so that each chain contains approximately 200 atoms, as this value is sufficient to predict the transport properties of polymers [78].The samples of matrices were then obtained from 10 polymer chains using the "Polymer/Composite Constructor" module.They were subjected to geometry and relaxation optimization to bring them to an equilibrium state (T = 300 K), and an additional 100 ps annealing step was used in the NVT ensemble with an Andersen/Behrendsen thermostat [79,80].The models of the polymer chains and sample matrices are shown in Figure 3. Density is an important parameter that can affect the accuracy of polymer permeability results.In our simulations, the density of each sample is assumed to be fixed.Therefore, when we build polymer samples, their final sizes correspond to the experimental density values.The densities and fractions of free volume in the obtained samples, calculated using the "Structure Properties" module of the MULTICOMP package, are given in Table 1.It can be seen from the table that among the selected polymers, PET has the smallest free volume, which should lead to low gas permeability.Density is an important parameter that can affect the accuracy of polymer permeability results.In our simulations, the density of each sample is assumed to be fixed.Therefore, when we build polymer samples, their final sizes correspond to the experimental density values.The densities and fractions of free volume in the obtained samples, calculated using the "Structure Properties" module of the MULTICOMP package, are given in Table 1.It can be seen from the table that among the selected polymers, PET has the smallest free volume, which should lead to low gas permeability.

Gas Permeability of Selected Polymers
The gas barrier properties of PET, PE, PVDF, and PTFE were calculated using the "Permeability" module integrated into the MULTICOMP package.All calculations were performed using the PCFF force field [77].As a separate issue, we considered the sensitivity of our results to the choice of the valence force field.For this purpose, the solubility of oxygen and water molecules was calculated using VFF COMPASS [81] and DREIDING [82] (see Section S5).

Calculations of the Diffusion Coefficients
The diffusion coefficients for oxygen and water molecules were estimated from the slope of the lines approximating their MSDs calculated from the MD simulations (see Table 2).A total of 10-50 ns MD runs were performed for each material.The MSD dependencies obtained, shown in Figure 4, have strong fluctuations that can be related to the relatively small number of diffusion trajectories realized for a single molecule.This is limited by the size of the chosen polymer matrix structure (~2000 atoms in the simulation cell).However, we can state that all the dependencies reached the asymptotic regime within the computational time.

Calculations of the Diffusion Coefficients
The diffusion coefficients for oxygen and water molecules were estimated from the slope of the lines approximating their MSDs calculated from the MD simulations (see Table 2).A total of 10-50 ns MD runs were performed for each material.The MSD dependencies obtained, shown in Figure 4, have strong fluctuations that can be related to the relatively small number of diffusion trajectories realized for a single molecule.This is limited by the size of the chosen polymer matrix structure (~2000 atoms in the simulation cell).However, we can state that all the dependencies reached the asymptotic regime within the computational time.The observed trend (PTFE > PVDF > PE > PET) of the diffusion coefficients of both O2 and H2O molecules is determined by the characteristics of the interaction of dissolved molecules with polymers and the free volume values in the created samples.The results in Table 2 also allow us to conclude that our results are qualitatively in agreement with  2 also allow us to conclude that our results are qualitatively in agreement with the experiment.In terms of quantitative agreement, the values obtained are in agreement with the experimental data within one order of magnitude.In general, all diffusion coefficients are higher, which can also be attributed to the limited number of available diffusion trajectories [83].However, it should be noted that accuracy within one or two orders of magnitude is characteristic of experimental measurement of diffusion coefficients in polymers [84].Comparison with previous theoretical studies also shows that our results are within orders of magnitude [85].The small difference observed between the diffusion coefficients for the H 2 O and O 2 molecules has also been confirmed experimentally [86] and theoretically [87,88].Therefore, it can be concluded that the atomistic approach provides reliable qualitative and, within an order of magnitude, quantitative estimates for the diffusion coefficients of small molecules in polymers.

Calculations of the Solubility Coefficients
The combined GCMC/MD approach is used to estimate the solubility coefficient S. The MC method is used to model the H 2 O and O 2 adsorption processes in the polymer matrix samples described above (see Section 3.1).Solute molecules can be created and destroyed with random probabilities, and the well-known Metropolis algorithm is used to accept or reject these configurational moves of a penetrant molecule [89].The GCMC simulation is invoked every 100 times.The motion of the inserted molecules is described in the subsequent MD run, which is performed to obtain the optimized configuration of the matrix/penetrant solution.These steps are used to reach the equilibrium state at a given temperature and chemical potential.When the penetrant insertion pressure p is changed, and thus the chemical potential of the system, we can obtain the dependence for the equilibrium number of penetrant molecules in the polymer N cell (p) and extract the solubility coefficient by linear approximation.The calculated N cell (p) for oxygen and water molecules are shown in Figure 5, with each point averaged over 10 ns MD-GCMC runs.The experimental values [84] and estimated values of the solubility coefficients are given in Table 3.
within an order of magnitude, quantitative estimates for the diffusion coefficients of small molecules in polymers.

Calculations of the Solubility Coefficients
The combined GCMC/MD approach is used to estimate the solubility coefficient S. The MC method is used to model the H2O and O2 adsorption processes in the polymer matrix samples described above (see Section 3.1).Solute molecules can be created and destroyed with random probabilities, and the well-known Metropolis algorithm is used to accept or reject these configurational moves of a penetrant molecule [89].The GCMC simulation is invoked every 100 times.The motion of the inserted molecules is described in the subsequent MD run, which is performed to obtain the optimized configuration of the matrix/penetrant solution.These steps are used to reach the equilibrium state at a given temperature and chemical potential.When the penetrant insertion pressure p is changed, and thus the chemical potential of the system, we can obtain the dependence for the equilibrium number of penetrant molecules in the polymer Ncell(p) and extract the solubility coefficient by linear approximation.The calculated Ncell(p) for oxygen and water molecules are shown in Figure 5, with each point averaged over 10 ns MD-GCMC runs.The experimental values [84] and estimated values of the solubility coefficients are given in Table 3.It should be noted that in our calculations we consider relatively small matrices consisting of ~2000 atoms.This was done to minimize the required computational resources.For the system sizes used, the solubility of the order of S ~10 −7 (cm 3 (GC)/(cm 3 Pa)) corresponds to only 1-10 molecules of the penetrant gas in the simulation cell (for the interval considered p).In the present study, we considered p ~10-1000 atm, because for this pressure interval, for the chosen system size, the dependence of N cell (p) is linear so that Henry's law can be applied for solubility estimation.At lower pressures, gas molecules occupy only a small number of the most accessible positions for dissolution ("adsorption sites"), so a strong penetrant-matrix interaction can lead to a large overestimation of the coefficient S. Furthermore, for pressures p below 1 atm, for the prepared material samples, only one or zero penetrant molecules are present in the simulation cells during the productive runs, making the calculated solubility values unreliable.As we increase the insertion pressure p above 10 3 -10 4 atm, the molecules penetrate tiny pores with poor connectivity.Moreover, because of the decrease in distance between the penetrants, strong interactions occur between the dissolved molecules.
In the case of water, this can even lead to the formation of water clusters inside the polymer matrix at penetrant penetration pressures p ~500-1000 atm.Such phenomena are often observed in atomistic simulations with water molecules and can be attributed to the consideration of hydrogen bonding [89] and dipoles [90,91], which is expected for water molecules [92].The formation of large water clusters can also affect the structure of the polymer matrix by inducing local polymer swelling.These two effects, among others [93], lead to noticeable nonlinearities in the dependence of N cell (p) and can strongly bias the value of the solubility coefficient.Thus, the values of the solubility coefficient can be very sensitive to many parameters, especially those defining the interactions between the penetrant molecules and the polymer matrix.A detailed discussion of the sensitivity of the calculated solubility coefficients to the parameters of the chosen interatomic potentials is given in Section S5.
In summary, our calculations show that the obtained solubility coefficients for oxygen are almost one order of magnitude lower than the experimental values (except for PTFE, which is in very good agreement with the experiment; see Table 3).Overall, the resulting trend in the value of the solubility coefficient for fluorine-containing materials (obtained in Section 3.2.1,i.e., PTFE > PVDF > PE > PET) is qualitatively in agreement with both experiments.At the same time, as we can see, the solubility coefficients for water are typically higher than the experimental ones.We attribute the latter fact to the effects of the water-water interaction that leads to cluster formation and the observed overestimation of the solubility coefficient.Nevertheless, as for the diffusion coefficient estimates, we can conclude that the atomistic approach provides reliable qualitative and, within an order of magnitude, quantitative estimates for the solubility coefficients of small molecules in the selected polymerics.

Calculations of Permeability Coefficients
The diffusion and solubility coefficients obtained for the materials considered allow us to calculate the permeability coefficients using Equation (15).Their product and comparison with the experimental and QSPR results are summarized in Table 4.As can be seen, PET has the lowest value of the permeability coefficients, which correlates well with its smallest free volume (see Table 1).All permeability coefficients for oxygen molecules are in good agreement with the experimental data and qualitatively agree with the results of the Bicerano model.At the same time, the coefficients P(H 2 O) do not agree with the Askadskii model.Furthermore, the Askadskii model does not follow the trend observed experimentally for the selected materials.Such a strong discrepancy in our results may be related to the difficulty of correctly describing water/polymer interactions.It is generally accepted that relatively small changes in the chemical nature of the polymer chain, which affect the abundance of available polar substituents, can have marked effects on water transport [94].Although hydrocarbon polymers exhibit low water uptake and permeability, polar substituents on polymer chains can significantly increase water uptake by forming specific interactions with water molecules [95].Unlike water, gases such as O 2 interact very little or not at all with the absorbent, so the matrix structure does not undergo swelling strain or other rearrangements.Since the interactions between O 2 molecules and O 2 /polymer are less complex, the coefficients of solubility, diffusion, and permeation are much easier to calculate with atomistic simulations [22].A more accurate selection of the modeling parameters (system size, valence force field parameters, etc.) and the extension of the set of statistics will significantly increase the accuracy and, accordingly, the predictive capabilities of the developed model.
On the basis of these conclusions, the predictions of the QSPR models for oxygen and water permeability were selected for the parameterization of the continuous model.The QSPR results appear to be more reliable because the models were trained on a fairly large representative sample of polymers.Nevertheless, the full-atom approach proved to be a useful tool for semi-quantitative predictions, especially in the absence of representative experimental data.

Results Based on the Continuous Model
This section presents the results of the analysis of the barrier properties of multilayer protective coatings using the model described in Section 2.2.
In the previous section, four polymers (PET, PE, PVDF, and PTFE) were selected and tested as components of a multilayer coating.As can be seen from the results obtained, they have relatively high water permeability values.Other polymer matrices (see Table S2) have slightly different vapor permeability values, which limits the ability to improve the barrier properties of polymer films by modifying their chemical structures.Therefore, it can be suggested that the main factor in improving the barrier properties of polymer films is the use of hybrid layers with inorganic fillers.
A four-layer stack, shown in Figure 1, is considered as a model of the flexible protective coating.It consists of a 30 µm polytetrafluoroethylene (PTFE) top layer, a 10 µm adhesive layer, a 60 µm polymer layer such as polyethylene (PE) with inorganic fillers, and a 60 µm substrate layer (PET) [96].Polymer parameters calculated with the Bicerano and Askadskii methods (see Table 4) were used to parameterize this model.To investigate the effect of inorganic fillers on the permeability of hybrid coatings, calculations were performed with different filler volume contents and degrees of their orientational ordering.Inorganic fillers were given in the form of thin disks and were characterized by the aspect ratio α = 2R/H, where R is the radius and H is the thickness of the particle.The orientation of the fillers was characterized by the direction of the normal to the disk surface and was set to be uniform in a cone with an angle θ max relative to the normal to the film surface (the isotropic distribution of fillers corresponds to the angle θ max = π/2).
The lateral size of the system was chosen so that the inorganic layer contained at least hundreds of filler particles.The number of particles over which the diffusion time was averaged in the Monte Carlo method was 500, ensuring an error of less than 10% in the permeability calculations for a given filler distribution.Diffusion time calculations are averaged over three different filler configurations at a given filler volume fraction.
First, the effect of the orientational order of the clay particles on the permeability of the hybrid organic-inorganic layer was investigated.Figure 6 shows the dependence of the ratio of the permeability of the filled polymer layer to the permeability of the pure polymer on the filler volume fraction, ϕ, at different values of the angle θ max = 12.5 • , 25 • , and 45 • with a fixed aspect ratio of the clay particle α = 60.As can be seen in Figure 6, the decrease in layer permeability is about 50% for ϕ = 5%, which is in good agreement with a large amount of experimental data (see, for example, the review [22]).A larger change in the permeability of the polymer layer containing disordered clay particles in some experiments [22] can be attributed to a change in the degree of crystallinity of the polymer matrix [97], which is not considered in this study.It can be seen that the calculated effect of a decrease in permeability has a linear dependence on low filler contents, which is consistent with the Nielsen model [56].It should be noted that the value of the overlap parameter αϕ, in this case, can be greater than one (αϕ = 3 for ϕ = 5%), but the Nielsen linear model describes the numerical results satisfactorily.The difference from the nonlinear dependence in the analytical model of Cassler et al. [56] for 2D layered structures may be due to the three-dimensional nature of the nanocomposite structure and the shape of the filler particles in the form of disks used in our calculations.matrix [97], which is not considered in this study.It can be seen that the calculated effect of a decrease in permeability has a linear dependence on low filler contents, which is consistent with the Nielsen model [56].It should be noted that the value of the overlap parameter αφ, in this case, can be greater than one (αφ = 3 for φ = 5%), but the Nielsen linear model describes the numerical results satisfactorily.The difference from the nonlinear dependence in the analytical model of Cassler et al. [56] for 2D layered structures may be due to the three-dimensional nature of the nanocomposite structure and the shape of the filler particles in the form of disks used in our calculations.Furthermore, Figure 6 shows that increasing the degree of orientational order of the layered filler significantly reduces the permeability of the nanocomposite and also allows the use of composites with a high filler volume fraction.Therefore, from a practical point of view, it is important to use a method to obtain a nanocomposite that allows one to Furthermore, Figure 6 shows that increasing the degree of orientational order of the layered filler significantly reduces the permeability of the nanocomposite and also allows the use of composites with a high filler volume fraction.Therefore, from a practical point of view, it is important to use a method to obtain a nanocomposite that allows one to obtain a high degree of orientational order of the filler.As a possible technique, layer-by-layer deposition of clay-polymer laminates can be used, as in [37,98,99].
Next, we investigate the variant of the protective coating that contains the hybrid layer of the clay-polymer laminate.The model structure used is shown in Figure 1b.It includes a 10 nm laminate layer and a 100 nm polymer PE layer.
The calculated permeability of a barrier layer with a clay/polymer laminate as a function of the number of bilayers is shown in Figure 7 for different sizes of clay particles (with R of 1, 3, 10, and 30 µm).The permeability of the barrier layer without the laminate was approximately 2.5 × 10 −1 (g/m 2 •day).As can be seen, a layer of 30 bilayers of laminate for clay particles with R of 30 µm leads to a decrease in vapor permeability by more than 500 times, which is sufficient for use in protective coatings for the long-term storage of food and medicine.In addition, clay/polymer laminates can be used as organic/inorganic barrier layers between polymer layers.From the results obtained, it can be estimated that 10 polymer/clay bilayers (with 10 µm clay particles) provide a vapor permeability of 10 −2 g/(m 2 •day).This is in qualitative agreement with experimental data on the decrease in permeability of nanolaminate hybrid films [37].
Thus, to obtain the maximum permeability reduction effect from the addition of inorganic fillers to the polymer matrix, it is necessary to use filler particles with the highest aspect ratio and degree of their orientational order.The product of the aspect ratio of the disc-shaped filler and its volume fraction must exceed 10 to achieve a significant permeability reduction effect.This should be combined with the use of a polymer matrix with minimal oxygen and water permeability.
Thus, to obtain the maximum permeability reduction effect from the addition of inorganic fillers to the polymer matrix, it is necessary to use filler particles with the highest aspect ratio and degree of their orientational order.The product of the aspect ratio of the disc-shaped filler and its volume fraction must exceed 10 to achieve a significant permeability reduction effect.This should be combined with the use of a polymer matrix with minimal oxygen and water permeability.

Discussion
Our study shows that the full atomistic methods allow us to semi-quantitatively describe the permeability of the polymers for oxygen molecules, which interact weakly with the polymer matrix.All obtained results are in reasonable agreement with the experiment and with the predictions of the Bicerano model.However, in the case of the water vapor permeability of polymers, our results show a strong dependence on the water/polymer model of interatomic potentials, which explains the observed discrepancy.
At the same time, the predictions of water permeability coefficients using the Askadskii method are not in good agreement with the experimental data.This can be explained both by the use of a limited amount of experimental data to fit the parameters of this regression scheme and by the presence of parameters that are not taken into account.These are the density and degree of crystallinity of the matrix, which depends on the preparation method, and the history of the material samples.The latter is difficult to account for in

Discussion
Our study shows that the full atomistic methods allow us to semi-quantitatively describe the permeability of the polymers for oxygen molecules, which interact weakly with the polymer matrix.All obtained results are in reasonable agreement with the experiment and with the predictions of the Bicerano model.However, in the case of the water vapor permeability of polymers, our results show a strong dependence on the water/polymer model of interatomic potentials, which explains the observed discrepancy.
At the same time, the predictions of water permeability coefficients using the Askadskii method are not in good agreement with the experimental data.This can be explained both by the use of a limited amount of experimental data to fit the parameters of this regression scheme and by the presence of parameters that are not taken into account.These are the density and degree of crystallinity of the matrix, which depends on the preparation method, and the history of the material samples.The latter is difficult to account for in QSPR models since they correspond to an idealized (highly averaged) structure of polymer materials.Nevertheless, such models are in great demand for the explicit prediction of the properties of pure polymers, which is necessary for the parameterization of the various theoretical models.They can also be useful for rapid quantitative evaluations of the permeability properties of components of designed protective coatings in laboratory studies.
The analysis of the transport properties of hybrid films with layered inorganic fillers using our continuous model has shown that, in the case of the use of orientationally disordered fillers, their effect makes it possible to reduce the permeability of the polymer film several times.However, it is difficult to achieve better results because of the mutual overlapping of the filler particles.Therefore, it is more promising to use nanolaminate materials with orientationally ordered clay particles.In this case, our model predicts a strong decrease in the permeability of the protective multilayer film (up to two or three orders of magnitude), especially when nanoparticles with a high aspect ratio are used (of the order of 1000).This theoretical result is in reasonable agreement with the experimental data for nanolaminate systems [37].

Conclusions
In this work, a multiscale computational model was developed to predict the barrier properties of multilayer protective films consisting of alternating polymer and hybrid layers containing clay minerals as fillers.Oxygen and water molecules were used as penetrant molecules.The model combines three levels of calculation, namely continuum, all-atom, and QSPR.The use of the continuum level makes it possible to take into account both the actual scale of the structure of the protective film and the dimensions and aspect ratio of the clay particles, which is very difficult to do with all-atom MD and MC simulations.To quickly parameterize the continuum model, the possibilities of using calculation schemes based on QSPR and full-atom models were considered.They allow the estimation of the transport properties of polymer materials for gas molecules (i.e., diffusion, solubility, and permeability coefficients).In testing the full atomic calculation scheme that combines MD and GCMC methods, we focused on considering relatively small volumes of matter, because of the need to minimize computational costs in laboratory studies.
We can conclude that the full atomic level of simulations allows us to make reliable, qualitative, and, within an order of magnitude, quantitative estimates of the diffusion and solubility coefficients of oxygen for relatively small polymeric samples.This allows good predictions of the oxygen permeability of pure polymers in agreement with the results of the Bicerano model and experimental studies.At the same time, we obtained evidence that existing QSPR schemes need to be improved to better predict water vapor permeability.This can be achieved by extending and increasing the base of polymers used to train the QSPR methods.In the case of all-atom simulations, it is necessary to make a more accurate choice of the sizes of the polymer samples and to adjust the parameters of the valence force field responsible for describing the water/polymer interactions.This requires separate studies that were not included in the objectives of this work.
Nevertheless, using the parameterization obtained and the continuum model, we were able to demonstrate that the maximum effect on the reduction in permeability from the addition of clay fillers to the polymers is obtained when using minerals in which the filler particles have high aspect ratios.It is also desirable to ensure a high degree of orientational order.We have shown that the product of the aspect ratio of the disk-shaped filler and its volume fraction must be greater than 10 to obtain a noticeable effect in reducing permeability.In addition, this should be combined with the use of polymer matrices with minimal oxygen and water permeability.
The advantage of our computational approach is its relative universality.This applies to the use of any low-molecular-weight penetrants, both for the development of barrier coatings and for membranes.At the same time, it is a closed-loop approach, since it can work without the use of external parameterizations as all the necessary coefficients (diffusion, solubility, and permeability) can be obtained through express analysis based on QSPR methods, and direct full-atom molecular dynamics and Monte Carlo simulations.The latter is important when using newly developed polymeric materials.
The proposed computational scheme can be used to improve the properties of protective coatings for food and drug storage and may also be useful in regulating the gas permeability of new polymeric materials.

Figure 1 .
Figure 1.Models of a barrier layer based on hybrid organic-inorganic multilayer films: (a) orientationally disordered filler particles and (b) laminate structure with orientationally ordered filler particles.Different colors of the layers schematically show the use of different materials for the formation of the multilayer protective films.

Figure 1 .
Figure 1.Models of a barrier layer based on hybrid organic-inorganic multilayer films: (a) orientationally disordered filler particles and (b) laminate structure with orientationally ordered filler particles.Different colors of the layers schematically show the use of different materials for the formation of the multilayer protective films.

Figure 2 .
Figure 2. General scheme for calculating permeability using combined MD and GCMC modeling.Figure 2. General scheme for calculating permeability using combined MD and GCMC modeling.

Figure 2 .
Figure 2. General scheme for calculating permeability using combined MD and GCMC modeling.Figure 2. General scheme for calculating permeability using combined MD and GCMC modeling.

Figure 3 .
Figure 3. Chain models of selected polymers (a) polyethylene (PE), (b) polyvinylidene fluoride (PVDF), (c) polytetrafluoroethylene (PTFE), (d) polyethylene terephthalate (PET), and material samples built on their basis.The grey and blue colors in the simulation cells show the outer and inner sides of the Connolly surfaces constructed by rolling a test ball of radius 1A.

Figure 4 .
Figure 4. Calculated MSD curves obtained for (a) oxygen and (b) water molecules for PET, PE, PVDF, and PTFE at 300 K.Each curve is the result of averaging 10 MD runs.Dashed lines show the value of the diffusion coefficient estimated by a linear approximation.

Figure 4 .
Figure 4. Calculated MSD curves obtained for (a) oxygen and (b) water molecules for PET, PE, PVDF, and PTFE at 300 K.Each curve is the result of averaging 10 MD runs.Dashed lines show the value of the diffusion coefficient estimated by a linear approximation.The observed trend (PTFE > PVDF > PE > PET) of the diffusion coefficients of both O 2 and H 2 O molecules is determined by the characteristics of the interaction of dissolved molecules with polymers and the free volume values in the created samples.The results in Table 2 also allow us to conclude that our results are qualitatively in agreement with the experiment.In terms of quantitative agreement, the values obtained are in agreement with the experimental data within one order of magnitude.In general, all diffusion coefficients are higher, which can also be attributed to the limited number of available diffusion trajectories[83].However, it should be noted that accuracy within one or two orders of magnitude is characteristic of experimental measurement of diffusion coefficients in polymers[84].Comparison with previous theoretical studies also shows that our results are within orders of magnitude[85].

Figure 5 .
Figure 5. Calculated dependencies of the number of (a) oxygen and (b) water molecules in the simulation cell N cell (p) for the selected polymer materials PET, PE, PVDF, and PTFE obtained using the PCFF force field.

Figure 6 .
Figure 6.The calculated dependence of the ratio of the permeability of the polymer layer with fillers P to the permeability of a pure polymer layer Pp on the filler volume fraction with the aspect ratio of α = 60 at different values of the angle θmax = 12.5°, 25°, and 45°.

Figure 6 .
Figure 6.The calculated dependence of the ratio of the permeability of the polymer layer with fillers P to the permeability of a pure polymer layer P p on the filler volume fraction with the aspect ratio of α = 60 at different values of the angle θ max = 12.5 • , 25 • , and 45 • .

Figure 7 .
Figure 7. Vapor permeability of the modified protective coating with a clay-polymer laminate layer as a function of the number of laminate bilayers for clay particle sizes of 1, 3, 10, and 30 µm.

Figure 7 .
Figure 7. Vapor permeability of the modified protective coating with a clay-polymer laminate layer as a function of the number of laminate bilayers for clay particle sizes of 1, 3, 10, and 30 µm.

Table 1 .
Densities and free volume fractions for PET, PE, PVDF, and PTFE.

Table 2 .
Diffusion coefficients of oxygen and water molecules for PET, PE, PVDF, and PTFE.

Table 3 .
Solubility coefficients of oxygen and water molecules for PET, PE, PVDF, and PTFE.S(O

Table 4 .
Permeability coefficients of oxygen and water molecules for PET, PE, PVDF, and PTFE.