Computational Screening of Diffusive Transport in Nanoplatelet-Filled Composites: Use of Graphene to Enhance Polymer Barrier Properties

: Motivated by the substantial interest in various ﬁ llers to enhance the barrier properties of polymeric ﬁ lms, especially graphene derivatives, we perform a computational screening of obstructed di ﬀ usion to explore the design parameter space of nanoplatelet-ﬁ lled composites synthesized in silico. As a model for the nanoplatelets, we use circular and elliptical nonoverlapping and impermeable ﬂ at disks, and di ﬀ usion is stochastically simulated using a random-walk model, from which the e ﬀ ective di ﬀ usivity is calculated. On the basis of ∼ 1000 generated structures and di ﬀ usion simulations, we systematically investigate the impact of di ﬀ erent nanoplatelet characteristics such as orientation, layering, size, polydispersity, shape, and amount. We conclude that the orientation, size, and amount of nanoplatelets are the most important parameters and show that using nanoplatelets oriented perpendicular to the di ﬀ usion direction, under reasonable assumptions, with approximately 0.2% (w/w) graphene, we can reach 90% reduction and, with approximately 1% (w/w) graphene, we can reach 99% reduction in di ﬀ usivity, purely because of geometrical e ﬀ ects, in a defect-free matrix with perfect compatibility. Additionally, our results suggest that the existing analytical models have some di ﬃ culty with extremely large aspect ratio (extremely ﬂ at) nanoplatelets, which calls for further development.


INTRODUCTION
There is substantial interest in nanoplatelet-filled (bio)polymeric composites because of their barrier properties for obstructing the transport of gas, vapor, and liquid. We are concerned in particular with graphene and graphene derivatives for their potential of enhancing barrier properties, which some of the authors currently investigate experimentally. 1−3 Graphene, a 2D carbon monolayer forming a hexagonal lattice, possesses exceptional mechanical, thermal, and optical properties, high crystal and electronic quality, and extremely high surface area. 4,5 Graphene and its many derivatives have emerged as some of the most highly promising material classes of the future, with applications in energy storage, 6 electronics and optoelectronics, 7 biological and chemical sensors, 8 environmental decontamination and water desalination, 9 and many others. 10 There is a rather comprehensive literature on graphene/polymer nanocomposites and their permeability. This covers many different types of polymers, e.g., poly(lactic acid), poly(ethylene terephthalate), poly(vinyl chloride), polystyrene, cellulose, poly(vinyl alcohol), and poly-(ethylenimine), and many different types of graphene derivatives, e.g., graphene, various forms of (reduced) graphene oxide, and exfoliated graphite. A plethora of processing conditions and weight/volume fractions lead to reported results on a reduction in the (gas) permeability ranging from a few percent to above 99.94%. The amount of reduction depends on the chemistry as such but also on purely geometrical characteristics such as the heterogeneity and orientation of the graphene-based obstacles in the polymer matrix, determining the degree of tortuosity, i.e., the lengths of the diffusive pathways through the material. 11 −32 In Figure 1, examples of polymer nanocomposite morphologies studied by means of a digital scanning electron microscope (Carl Zeiss DSM 940, Carl Zeiss AG, Oberkochen, Germany) are shown. The polymeric matrix consists of a commercial low-density polyethylene (LDPE; M w = 92 kg/mol, PI = 7.6, and T m = 111°C; Borealis AB, Stenungsund, Sweden). The nanofillers are two commercial types of graphite nanoplatelets (GnPs; XG Sciences, Lansing, MI) with 5 and 25 μm mean diameters. The nanocomposites are manufactured via extrusion processing, resulting in highly oriented nanofillers in the extrusion flow direction. [1][2][3]33 In connection with these experiments, we are interested in gaining an understanding of the effect of the different material properties and processing conditions by means of simulation.
Some molecular-dynamics-based simulation studies on graphene oxide membranes, nanoporous graphene, and stacked layers of graphene sheets and their molecular-level interactions with a permeating species have been performed. 34−39 At the mesoscale level, more relevant for our work, purely geometrical obstruction effects on the diffusion/permeability have been studied in both 2D and 3D. These studies use finite-element and grid-based methods as well as theoretical methods to solve the diffusion and Laplace equations for the local chemical potential. There is, for the most part, a focus on round platelets with aspect ratios α ranging from 3 to 1000, the volume fraction, orientational distributions, random and ordered configurations, multiscale approaches to account for diffusion inside lamellar obstacles, and the impact of interaction between the polymer and filler. 31,40−56 In this work, we perform a computational screening of obstructed diffusion to explore the design parameter space of nanoplatelet-filled composites synthesized in silico. As a model for the nanoplatelets, we use circular and elliptical nonoverlapping and impermeable (with solubility 0 and without defects) flat disks (with an infinite aspect ratio, i.e., infinitely flat). This provides a simple model of graphene-based nanoplatelets dispersed in a homogeneous, isotropic, polymer matrix under the assumption of perfect compatibility between the matrix and filler, i.e., that the proximity to a filler particle does not influence the properties of the matrix through interactions or nucleation of crystal structures (implying that the diffusivity controls the permeability entirely). It is obvious that inhomogeneities in the matrix can impact the diffusivity, 31,56 but we focus on purely the geometrical effects of nanoplatelets on diffusion in this work. Diffusive transport of point particles is stochastically simulated using a random-walk model from which the effective diffusivity is calculated. On the basis of ∼1000 simulated structures and their corresponding simulated effective diffusivities, we systematically investigate the impact of different nanoplatelet characteristics, such as the angular orientation, layered structures, size, polydispersity, shape, and total amount. The aim of this computational screening paradigm is to compare the relative impact of these

ACS Applied Nano Materials
Article parameters. The ambition is to explore by a computational screening the effect of varying different geometrical parameters independently and to discover and understand the generic design rules for graphene−polymer nanocomposites and the tailoring of mass-transport properties. To our knowledge, this joint study of many different parameters has not been done before for these materials, and this effort will aid in guiding future experimental work.

METHODS
2.1. Structure Generation. Random configurations of nonoverlapping elliptical flat disks (including circular disks as a special case) are generated in a cubic simulation domain with side L = 100 μm (the algorithms are all scale-independent, but a scale on the order of 100 μm is a realistic setting for our problem) using a hard-particle Markov chain Monte Carlo (MCMC)-type algorithm. First, disks are placed randomly (and possibly overlapping), either uniformly distributed in the whole simulation domain with random orientations or subject to some constraints, see below. Second, the configurations are relaxed, iteratively performing random translations and rotations of all particles until all overlaps have been removed. Finally, the configurations are equilibrated, performing a large number of random translations and rotations of all particles ensuring a distribution in the location and orientation that is as uniform as possible. The overlap criterion is based on the Perram−Wertheim overlap criterion 57 for two ellipsoids of arbitrary orientation, reduced to the "degenerate" case of two ellipses considered as flat ellipsoids with one semiaxis equal to zero (the resulting overlap criterion is well-defined except in the case of two coplanar ellipses for which the intersection is a single point, but this is immaterial for simulation purposes). The algorithm is implemented in Julia (www.julialang.org), 58 and the code is available in a Github repository (https://github.com/roding/whitefish_ generation, version 0.1). On a dual Intel Xeon E5-2699 v4 setup, the execution time is, on average, ∼1 min (single thread).
We generate several series of structure data sets to study the impact of different structural parameters, i.e., orientation, layering, number of disks, polydispersity, shape, and total surface area (by which we mean the sum of the surface areas of the disks, not counting both sides). Some structures are anisotropic; we have defined the z axis as the direction "through" the material, and even though we perform 3D diffusion simulations, the effective diffusivity will later be calculated in the z direction. We study the following (examples of structures from the different data sets are shown in Figure 2): (i) Orientation: For a total surface area of 10 5 μm 2 and 500 disks (with radius ∼8 μm), the maximum angular deviation relative to the z axis is varied between 0 (all disks orthogonal to the axis; the angular constraint is with respect to the normal vector of the disk) and π/2 (free orientation).
(ii) Layering: For a total surface area of 10 5 μm 2 and 500 disks (with radius ∼8 μm), the disks are compressed into a layer centered in the simulation domain with the thickness along the z axis varied between 25 and 100 μm (the latter meaning no constraint).
(iii) Number: For a total surface area of 10 5 μm 2 , the number of disks is varied between 100 and 1000, hence distributing the same total surface area differently and changing the radius from ∼5.5 to ∼18 μm.
(iv) Polydispersity: For a total surface area of 10 5 μm 2 and 500 disks, the surface areas of the disks are log-normal-distributed with a coefficient of variation (i.e., ratio of the standard deviation and the mean) between 0 and 1 (the latter being a very broad distribution). Because randomly sampling areas will create a random variation in the total surface area, we normalize the total surface area to 10 5 μm 2 in order to isolate the effect of polydispersity.
(v) Shape: For a total surface area of 10 5 μm 2 and 500 disks, the semiaxis ratio is varied between 1 and 10, ranging from circular disks to very elongated elliptical disks.
(vi) Total surface area: For a constant radius of ∼8 μm as above, the total surface area is varied from 5 × 10 4 to 5 × 10 5 μm 2 by varying the number of disks from 250 to 2500.
Finally, we also study the combined effect of some of these parameters with a discussion toward practical feasibility and usefulness.
2.2. Diffusion Simulation. Diffusion in the generated structures is simulated using a "random-walk particle tracking" technique. 59,60 An ensemble of N = 4 × 10 6 diffusing point particles is initially uniformly distributed in the simulation domain. Stochastic particle motion is simulated as a Gaussian random walk with time resolution δt = 0.0025 arbitrary units (a.u.) and diffusion coefficient D 0 = 1 a.u., hence adding random normal distributed displacements to the position in each time step with zero mean and standard deviation δ D t 2 0 (assuming that D 0 is constant corresponds to assuming a perfect compatibility between the matrix and filler, i.e., that the proximity to a filler particle does not influence the properties of the matrix). The position is recorded at each major time step Δt = 0.25 a.u. The simulation proceeds up to t max = 2000 a.u. (for all structures except the layering data set) or 5000 a.u. (the layering data set). Single rejection boundary conditions 61 are used; a proposed displacement is only accepted if it does not pass through any disk(s) (this is equivalent to zero solubility within the disks). A time-dependent effective diffusion coefficient (i.e., obstruction factor) in the z direction is computed as where z n (t) is the z position of the nth particle at time t. The effective diffusivity is obtained as the asymptotic value D ∞ /D 0 ≤ 1. The finiteness of δt implies that the probability of accepting a displacement is smaller than 1 and hence that, effectively, D 0 < 1 in practice (or, more precisely, D δt < 1). We compensate for this by dividing all estimated effective diffusion coefficients with the corresponding "empirical" D 0 (D δt ) as obtained in that structure (the impact on the results is, in most cases, negligible, however). An advantage of this simulation technique is that the nanoplatelets can be exactly represented without discretization and approximation, and depending on the choice of time resolution δt, the diffusive transport can be simulated with arbitrary precision. The algorithm is provided in a parallel implementation in Julia (www.julialang.org), 58 and the code is available in a Github repository (https://github.com/roding/ whitefish_diffusion, version 0.1). On a dual Intel Xeon E5-2699 v4 setup, the execution time is, on average, ∼2.5 h (88 threads). In Figure  3, some examples of computed effective diffusivity curves are shown.

RESULTS AND DISCUSSION
First, we investigate the effect of the distribution of angles relative to the z direction. We impose the constraint that the maximum angular deviation relative to the z axis is varied between 0 (all disks perpendicular to the axis; the angular constraint is with respect to the normal vector of the disk) and

ACS Applied Nano Materials
Article π/2 (free orientation). As shown in Figure 4, the effective diffusivity increases as the angular constraint is relaxed. The enhancement in the barrier properties (i.e., the reduction from the effective diffusivity equal to unity) is ∼15% for the randomly oriented case and ∼45% for the perfectly perpendicular case. Hence, these results are in line with the predictions of Fredrickson and Bicerano, 42 stating a 3-fold improvement when moving from random to perfectly aligned configurations. Second, we investigate the effect of having a layered or a nonlayered structure. The disks are compressed into a layer centered in the simulation domain with the thickness along the z axis varied between 25 and 100 μm (the latter meaning no constraint, i.e., that the disks are completely uniformly distributed). As shown in Figure 5, the effect is quite small, with a slightly decreasing effective diffusivity as the disks are increasingly uniformly distributed. This small effect could be because of increased ordering in the structure when the disks are more tightly packed; indeed, whereas the mean angular deviation from the z axis is constant (∼π/2), the standard deviation of angular deviations increases from ∼0.45 to ∼0.55 rad when the thickness goes from 25 to 100 μm, indicative of a less angular ordering for a less compressed layer.
Third, we investigate the effect of the number of disks, varied between 100 and 1000 for a constant total surface area. Hence, the surface area is distributed over a different number of disks, thereby varying the radius from ∼5.5 to ∼18 μm. As shown in Figure 6, the effect is quite pronounced, with lower effective diffusivities for a few large disks than for many small ones. Apparently, large obstacles more efficiently block diffusion, whereas with small obstacles, more possible diffusive pathways are available.
Fourth, we investigate the effect of the polydispersity or size distribution of the disks, with the surface areas of the disks being log-normal-distributed with a coefficient of variation (i.e., ratio of the standard deviation and the mean) between 0 and 1 (the latter being a fairly broad distribution). By normalization to a fixed mean area and total surface area, the effect of the polydispersity can be studied completely independently of the mean disk area. As shown in Figure 7, the effect is rather small, with slightly decreasing effective diffusivity for increasing polydispersity. We attribute this to a few large disks efficiently blocking the diffusion. The spread in the effective diffusivity is also increasing. This is an effect of the larger differences between statistically equal polydisperse configurations than between monodisperse ones (the latter differing only in the localization and orientation, not in where the large and small disks are). In conclusion, the polydispersity is, to some extent, giving the same effect as increasing the mean disk area, although the latter has more effect. This is in line with work by Lape et al., 62 which predicts a decrease in the diffusivity with increasing polydispersity.

ACS Applied Nano Materials
Article Fifth, we investigate the effect of the shape by varying the semiaxis ratios of the disks between 1 and 10, ranging from circular disks to very elongated elliptical disks. As shown in Figure 8, more elongated elliptical disks are less efficient as barriers than circular disks. This actually came as a surprise; because elongated disks are more "entangled", we believed at first that they should give more tortuous diffusive pathways. However, in reality, we observe the opposite effect, which can be understood by considering that it is easier to "diffuse around" particles that are very small in one dimension. This is thereby an effect related to the fact that smaller particles in larger numbers block diffusion less efficiently as concluded earlier.
Sixth, we investigate the effect of the total surface area. For a constant radius of ∼8 μm, as in most investigations above, the total surface area is varied from 5 × 10 4 to 5 × 10 5 μm 2 by varying the number of disks from 250 to 2500. As shown in Figure 9, the effect is quite pronounced, as would be obvious. The effect of the total surface area will be further investigated below.
After investigating these six parameters and their impact on the overall effective diffusivity, we proceed to a more in-depth analysis of a special, "best" case where we focus on circular, monodisperse disks, aligned perpendicular to the diffusion direction and with no layering; i.e., the disks are uniformly distributed in the simulation domain. For three different particle radii, 7.5, 12.5, and 17.5 μm, we perform simulations for the total surface areas of 10 5 , 2 × 10 5 , ..., 15 × 10 5 μm 2 (these simulations are more computationally demanding than before because of the increased number of disks, and, therefore, we perform a rather small number of them in this final simulation series). We study three different radii because the previous investigation made it clear that the radius is an important factor; the three chosen values are realistic values for graphene-based fillers (see Figure 1 and the description). In the simulations, the nanoplatelets are infinitely thin and, hence, have zero volume fraction; the obstruction effects are completely due to the surface area and not the volume. In order to make more sense of these final results in a physical context, assume that we have a generic polymer matrix with a density of 1 g/cm 3 and that we have an average of 10 layers of graphene in each nanoplatelet (the density of the graphene is 0.77 mg/m 2 , and it is noteworthy that <10 layers is typically considered "graphene" and >10 layers is typically considered GnPs). The results can be plotted as a function of the weight fraction, shown in Figure 10 and demonstrating as expected that the larger particles have superior barrier properties. In previous investigations, we only reached rather moderate obstruction effects (D ∞ /D 0 > 0.55), but now we see that the synergetic effect of the angular alignment and larger total surface area can provide more than 99% reduction in diffusivity purely because geometrical effects. We see that, for 17.5-μmradius disks and with approximately 0.2% (w/w) graphene, we can reach 90% reduction and, with approximately 1% (w/w) graphene, we can reach 99% reduction in diffusivity, purely because of geometrical effects. Of course, the fewer layers of graphene on average, the smaller the weight fraction of the filler necessary to obtain this level of obstruction. Furthermore, some comparison with analytical models for diffusivity (or rather permeability) is of interest. We considered three models for platelets/flakes aligned perpendicularly to the direction of diffusion and incorporating volume fraction ϕ and aspect ratio (diameter-to-thickness ratio) α, namely, the Nielsen model, 63 the Lape−Nuxoll−Cussler model, 62 Figure 8. Effective diffusivity as a function of the semiaxis ratio for elliptical disks, ranging from 1 (circular disks) to 10 (very elongated elliptical disks).
Assuming that we have 10 layers of graphene, each 0.335 nm thick, in each nanoplatelet, the thickness is 3.35 nm and, hence, the aspect ratios for the three different particle diameters are between 4500 and 10500. As a function of the weight fraction, the predictions of these three models are also plotted in Figure  10. None of them fit particularly well. However, if the aspect ratio α is treated as a fitting parameter that can be tweaked rather than as a physical parameter that is a known constant, a near-perfect fit can be found for the Lape−Nuxoll−Cussler model by choosing a value of α equal to ∼0.32 times its true value. The fact that the Lape−Nuxoll−Cussler model fits well only with a tweaked parameter and that none of the models fit well without tweaking may be indicative of a difficulty to capture the effects of very large aspect ratios in analytical models.

CONCLUSION
We have performed computational screening of the effective diffusivity and barrier properties in nanoplatelet-filled composites synthesized in silico. As a model for the nanoplatelets, we use circular and elliptical nonoverlapping and impermeable flat disks dispersed in a homogeneous, isotropic (polymer) matrix with constant solubility, assuming perfect compatibility between the matrix and filler. Exploring the design space of this model using ∼1000 simulated structures and effective diffusivity simulations, we assessed the importance of several geometrical parameters independently such as the orientation, layering, size, polydispersity, shape, and amount of nanoplatelets. We found that the most crucial parameters are, not very surprisingly, the angular orientation/alignment, size, and amount of nanoplatelets. Further investigation into these parameters demonstrated that, under reasonable assumptions, with approximately 0.2% (w/w) graphene, we can reach 90% reduction and, with approximately 1% (w/w) graphene, we can reach 99% reduction in diffusivity, purely because of geometrical effects, not relying on, e.g., crystal nucleation. Of course, the fewer layers of graphene on average, the smaller the weight fraction of the filler necessary to obtain this level of obstruction. The investigated model constitutes a rather simplified model of a polymer−graphene nanocomposite; nevertheless, generic design rules and their implications on the effective diffusivity can be understood through the use of this model and facilitate the tailoring of mass-transport properties of barrier materials, thus guiding future experimental work. There are promising directions to pursue in further work, such as characterization of the materials by means of spatial statistics for a further understanding of the critical geometrical features, e.g., tortuosity and correlation functions, and assessing the effect of an inhomogeneous matrix and defects in the nanoplatelets. Additionally, our results suggest that existing analytical models have some difficulty with extremely large aspect ratio (extremely flat) nanoplatelets, which calls for further development.