Finite Element Modeling of Multiscale Diffusion in Intercalated Nanocomposites

This work is aimed to study the diffusion in 3D nanocomposites obtained with stacks of lamellar nanofillers characterized by the presence of permeable galleries, by finite element (FE) analysis. To this purpose, a geometricmodel, based on a randomdistribution of noninterpenetrating stacks, with each one beingmade of regularly spaced lamellae, was developed.Thedevelopedmodel is able to account for diffusion between stacks (interstack diffusion) as well as diffusion inside stacks (intrastack diffusion). Simulation results showed that intrastack diffusion, related to flow inside galleries, can be quite relevant, particularly at high values of gallery thickness. Comparison of the simulation results with literature models shows that when intrastack diffusion is not taken into account, the diffusion behavior in intercalated nanocomposites is not well predicted. Therefore, intrastack permeability of nanofillers such as organic modified clays cannot be neglected. Such intrastack diffusivity is shown to depend on the morphological features of the nanofiller requiring the development of a proper mathematical model.


Introduction
The substantial decrease of permeability to mass transport brought about by the addition of nanofiller is a major advantage of polymer matrix nanocomposites.The improvement of the barrier properties is mainly attributed to the increase of the tortuous path that diffusing molecules must travel inside the diffusing medium [1,2].On the other hand, it is generally accepted that such increase depends on different parameters, among which are volume fraction, orientation, and aspect ratio of the nanofiller [3,4].
The mass diffusion process in polymer nanocomposites is very often studied by finite elements simulation and/or mathematical modeling accounting for the effects of different factors, some of which are very difficult to quantify and measure (i.e., orientation of nanofiller).Due to the difficulties in characterizing the structure of polymer nanocomposites by wide angle X-ray diffraction (WAXD) [5] and transmission electron microscopy (TEM), [6], simulation and modeling are very often used, together with experimental data, to obtain indirect measurements of the morphologic features of nanocomposites.Some mathematical models developed in the literature are aimed to predict the tortuous path of diffusing particles in polymer nanocomposites.Different models were developed, which are able to account for the high aspect ratio of the nanofillers and are valid in 2D or in 3D problems [1,7].All the developed models assume that the nanofiller is perfectly impermeable to the diffusing species.Only recently, the effect of increased volume of impermeable phase, strictly correlated to this assumption, has been introduced [8].Consequently, 2D [9] and 3D [10] analytical models have been developed, which are able to account for diffusion inside the galleries.
Alternatively and, in some cases, complementarily, the diffusion process in polymer nanocomposites has been studied by finite element analysis.In many cases, simulations were run on perfectly aligned particles, perpendicular or parallel to the direction of diffusion, modeling the nanofiller as discs [11,12] or prisms [13].However, perfect alignment of nanoparticles is an ideal situation not commonly found in real nanocomposites system.Only recently, the effect of random distribution of particles [14], or of particles oriented with a certain degree with respect to the direction of diffusion [15], has been introduced.In addition, 2D [9] and 3D [10] FEM models, able to account for diffusion inside the galleries, have been recently developed.
Accounting for diffusion inside stacks is of primary importance in order to simulate and model the mass transport behavior in intercalated nanocomposite, which retain the basal order of the nanofiller.The structure of the nanocomposite is characterized by the presence of lamellar stacks, whose aspect ratio depends on both the number of platelets which are, on the average, present in each stack and the distance between individual platelets.The two different morphological features of the nanofiller are defined as degree of dispersion and degree of intercalation, respectively.A strongly intercalated, weakly dispersed nanocomposite is characterized by the presence of thick galleries, with each stack being composed of many platelets; in contrast, a strongly dispersed, weakly intercalated nanocomposite is characterized by thin lamellar galleries, but the stacks are composed of a few number of platelets.
In such situation, diffusion should be considered at two different scales; the first one associated with diffusion among stacks (interstack diffusion), and the second one associated with diffusion inside stacks (intrastack diffusion) [16].The dual scale approach has been recently introduced in order to model the mechanical behavior of polymer nanocomposites [17,18].The reported results show the relevance of the degree of dispersion on the mechanical behavior of polymer nanocomposites.The dual scale diffusion problem is similar to the problem of polymer flow during fiber impregnation in infusion processes used in composite materials fabrication.In this case, two mechanisms of impregnation have been clearly identified: the first one, occurring at lower time scales, characterized by a higher value of permeability, as defined by the Darcy law, is associated with interbundle flow, whereas the second one, occurring at higher time scale, and therefore characterized by a lower value of permeability, is associated with intrabundle flow [19].The first mechanism involves flow of the molten polymer between fiber bundles (macroimpregnation), whereas the second one involves polymer flow inside each bundle (microimpregnation).
The aim of this paper is to study the 3D diffusion behavior of nanocomposites in the presence of permeable lamellar stacks by finite elements (FE) analysis.To this purpose, a geometric model, made of randomly and uniformly distributed stacks, was developed by the use of a Matlab routine.To the best of our knowledge, this is the first time that such a multiscale approach is introduced to study 3D diffusion in polymer nanocomposites.The results from simulation are compared with different analytical models, showing the relevance of intrastack diffusion.Finally, a proper model, which is able to predict the intrastack diffusion coefficient, is developed.

Modeling of Diffusion in Nanocomposite
The first paper about modeling of diffusion in multiphase systems, most commonly quoted, is by Maxwell, who solved the problem of a suspension of spheres sparsely distributed in the continuum and therefore characterized by a reasonably small volume fraction, V  .The model accounts for the presence of permeable particles, as reported in the scheme of Figure 1(a), characterized by a diffusivity   , dispersed in a permeable matrix, characterized by a diffusivity   [20]: in which  norm , the diffusivity of the multiphase system, is defined as the diffusion coefficient of a homogeneous medium, exhibiting the same steady-state behavior as the two-phase composite.Different modifications of the Maxwell model have been proposed in order to account for the shape of dispersed particles.Among these, the model developed by Fricke [21] is valid for oblate spheroids, characterized by two axes,  and , with  >  and  being the rotation axis of the ellipse: where with AR the aspect ratio of the nanofiller, /.On the other hand, modeling the diffusion in polymer nanocomposites must be able to account for the high aspect ratio of the nanofiller.The most commonly quoted approach was developed by Bharadwaj [1], whose model is able to account for the orientation angle of the nanofiller and, in the case of a random distribution, can be written as The Bharadwaj model assumes that nanofiller is perfectly impermeable to the diffusing gas, as reported in the scheme of Figure 1(b).Nevertheless, a perfect dispersion of the nanofiller is very difficult to attain.Single platelets of nanofiller can only be found in perfectly exfoliated nanocomposites.In most cases, corresponding to intercalated nanocomposites, stacks of platelets characterize the structure of the nanocomposite, as reported in the scheme of Figure 2, which represents a 2D scheme of the nanofiller.In 3D, the rectangle becomes a cylinder of diameter .
The morphological features of the stacks are the width or diameter of the platelets , the average number of platelets in each stack,   , the thickness of platelets   , the number of galleries   =   − 1, and the thickness of galleries,   .The number of platelets in each stack defines the degree of dispersion, whereas the gallery thickness defines the degree of intercalation.The total thickness of the stack   is given by According to the scheme of Figure 2, it is possible to define two different values for the aspect ratio.The first one is the aspect ratio of platelets, AR  , defined as whereas the second one is the aspect ratio of a stack, AR  = /  , which, also accounting for (5), is This equation highlights that the aspect ratio of the stack is dependent on both the degree of intercalation, through the parameter   /  , and the degree of dispersion, through the parameter   [9].Further, two different values of volume fraction can be defined.The first one, V  , is defined as the ratio between the volume of the nanofiller and the total volume of the nanocomposite.If   is the number of stacks per unit volume of the nanocomposite, then the volume fraction V  is The second one, V  , is defined as the ratio between the total volume of lamellar stacks, including the galleries volume, and the total volume of the nanocomposite, which is Taking the ratio between ( 8) and ( 9), the following correlation can be obtained between V  and V  : Recently, Nazarenko et al. [22] recognized that, in assuming impermeability of the nanofiller, the volume fraction of inorganic phase must be substituted by the total volume fraction of impermeable phase, also including lamellar galleries.In such case, the Bharadwaj model is corrected to where, according to ( 7)-( 10), Therefore, the Nazarenko model can be considered an extension of the Bharadwaj model, in which the volume fraction and aspect ratio of impermeable phase are considered, as reported in the scheme of Figure 3(a).Nazarenko model is based on the assumption that, at relatively low values of   /  , diffusion inside the nanofiller can be neglected.Nevertheless, as reported in the scheme of Figure 3(b), diffusion can occur either between stacks (interstack diffusion) or inside stacks (intrastack diffusion).This is a typical dual scale problem, which has been solved to yield the following model for the diffusivity [10]: Comparison between (11) and (13) shows that the latter model allows to account for the interstack diffusion (through the term V  AR  /  , homologous to that found in the Nazarenko model) and intrastack diffusion (through the term

Simulation of Dual Scale Diffusion
Finite element analysis Comsol Multiphysics software (version 3.5, Comsol AB, Sweden) was used to study the 3D diffusion in nanocomposites, characterized by the presence of intercalated lamellar nanofillers characterized by permeable stacks.The nanocomposite was modeled as an array of stacks of impermeable lamellae uniformly and randomly dispersed in the integration domain.The position and orientation of stacks inside the cubic domain were randomized by means of a properly developed subroutine of Matlab 6.5, which is based on a Monte Carlo stochastic approach and includes noninterpenetrating conditions between stacks [9].As reported in a previous work [10], the average distance between the centers of the stacks is in the range of 0.5-1 times the diameter of the stacks, depending on volume fraction, number of platelets in each stack, and aspect ratio of the platelets.
Each stack is composed of a regular arrangement of parallel and equidistant platelets, represented as cylinders.A sketch of the domains obtained is reported in Figure 4.Each stack can intercept the domain boundary surfaces parallel to the direction of diffusion.In such case, the actual volume of the stack which falls inside the domain is considered.
The problem is solved by means of the mass transport module, diffusion submodule, in stationary conditions, assuming that the diffusivity inside the galleries is equal to the diffusivity of the matrix.
This assumption does not take into account the anomalous effects which could arise at low values of   .At least two factors can contribute to such anomalous diffusion: (a) The existence of a rigid amorphous fraction [23] at the lamellar sheet interface, which can drastically reduce the diffusivity inside galleries.
(b) The presence of low molecular weight organic modifier inside galleries, which can increase the diffusivity.
However, changes of diffusivities inside the galleries can be readily implemented in the simulation, as will be discussed in the following section.The software solves the mass balance equation: coupled with the proper boundary conditions at each platelet surface with normal ⃗  (impermeability), at the boundaries of the domain (transparent boundaries in Figure 4) parallel to the direction of diffusion ( axis), at the upper boundary of the domain perpendicular to the direction of diffusion (upper gray boundary  1 in Figure 4), and at the lower boundary of the domain perpendicular to the direction of diffusion (lower gray boundary  2 in Figure 4).The concentration  1 was chosen in order to have, for all the simulations, a concentration gradient  1 / 0 = 10 3 mol/m 4 , with  0 being the length of the integration domain along the direction of diffusion ( axis in Figure 4).Some simulations were also run using 2D geometry.In such case, each platelet was represented as rectangle  *  dimensions.Randomization of the position was obtained using the same procedure described for 3D simulations.
The Matlab routine also includes a subroutine for the calculation of the average path length of massless particles diffusing due to concentration gradient.The subroutine is based on the scheme reported in Figure 5 for a two-dimensional problem.
The diffusion path is made by a series of points (  ,   ,   ), and its slope can be calculated, according to Fick's law according to the scheme of Figure 5, as where  is the flux and  is concentration.Therefore, the position of each point can be calculated as Finally, the diffusion path length is calculated as The normalized path length,  norm , was obtained as the ratio between the average path length of diffusing particles inside the nanocomposite and  0 , which corresponds to the diffusion length in the absence of lamellae: The normalized diffusion coefficient,  norm , was obtained as the ratio between the normal flow rate on the two boundaries perpendicular to the direction of diffusion ( 1 and  2 ) and the flow rate in the neat matrix: In order to account for the morphology of the nanofiller, simulations were run adopting three different nanofiller volume fractions (0.015, 0.0225, and 0.03), a number of lamellae in each stack changing from 3 to 5, and three different gallery thicknesses (1.5, 2.5 and 3.5 times the lamellar thickness).The aspect ratio of platelets was held at a constant value of 50.The values chosen for simulation are in the typical range of montmorillonite nanofilled polymers, which are characterized by a lamellar thickness of about 0.883 nm [24] and separated by lamellar galleries of thickness between 1.5 and 4 nm [25].Each result reported is obtained as the average of 15 simulations; each one runs in a domain comprising 10 stacks.

Results of Numerical Analysis
The results of 2D simulations are reported in Figure 6, showing the diffusion paths calculated according to (20).The blue ellipses evidence some of the intrastack diffusion zones, characterized by the fact that the diffusing particle travels inside the lamellar galleries.
The normalized coefficient of diffusion obtained by 3D FE simulation at different values of gallery thickness and V  = 0.015 and   = 3 is reported in Figure 7.The two sets of data are relative to simulation with impermeable stacks and simulation with permeable stacks (diffusion inside galleries).For comparison purposes, the coefficient of diffusion corresponding to null gallery thickness,   = 0, is also reported.Both models (either accounting or neglecting intrastack diffusion) converge at the same value when   = 0, that is, impermeable galleries.On the other hand, neglecting intrastack diffusion yields a significant underestimation of the coefficient of diffusion.The error associated with the impermeable stack approximation increases with increasing gallery thickness.
Initially, a sensitivity analysis was performed, running simulations using a variable intrastack diffusivity.The plot of Figure 8 reports the coefficient of diffusion of the nanocomposite as a function of the ratio between intrastack and interstack diffusivity, for V  = 0.015,   /  = 1.5, and   = 3. One-way analysis of variance (ANOVA) was performed using the intrastack diffusivity as the factor level.Considering the set of 9 data ranging between  intrastack / interstack = 0.6 and  intrastack / interstack = 1.4,an  value equal to 1.84 was calculated.This value was compared with the critical  CV (8, 126) = 2.01, with 8 being the degree of freedom of the factor level and 126 = 9 * 14 the degree of freedom associated to the error.The estimated  value is lower than the critical  value, indicating that the different simulations provide an almost constant value of the coefficient of diffusion.On the other hand, when the ANOVA is performed on the set or 11 data ranging between  intrastack / interstack = 0.5 and  intrastack /  interstack = 1.5, the calculated  = 2.63 is higher than the critical  CV (10, 154) = 1.89, which indicates that the intrastack diffusivity has a significant influence on the simulated coefficient of diffusion.Therefore, the ANOVA analysis indicates that the assumption of intrastack diffusivity equal to interstack diffusivity can be considered valid in the range 0.6 ≤  intrastack / interstack ≤ 1.4.
The normalized coefficient of diffusion is reported as a function of the gallery thickness for different values of the nanofiller volume fraction, at constant value of   = 3, in Figure 9.For comparison purposes, the normalized coefficient of diffusion is also reported for a null gallery thickness; that is,   = 0.
Due to the quite high error bars, one-way analysis of variance (ANOVA) was initially performed on the simulated coefficient of diffusion.Considering 3 different gallery thicknesses as factor levels and 15 simulations on each level, a critical value of the  ratio ( CV (2, 42)) equal to 3.22 at a  11) Model eq. ( 9) Simulation Model eq. ( 11) Model eq. ( 9) Simulation Model eq. ( 11) Model eq. ( 9) significance level  = 0.05 is obtained.The  ratio values calculated for the coefficient of diffusion simulated at V  = 0.015, V  = 0.0225, and V  = 0.03 are 6.85, 7.35, and 5.67, respectively.Therefore, in any case,  ratio is much higher than  CV , which indicates that the differences between the results of the simulations obtained at different gallery thickness are relevant at the 5% significance level.For comparison purposes, the simulated normalized diffusion coefficient is reported for a gallery thickness equal to 0 (corresponding to impermeable stack with aspect ratio equal to 16.7 according to (7)).The normalized coefficient of diffusion decreases with increasing lamellar spacing, which is in agreement with the generally accepted concept and many experimental results indicating that improved barrier properties are achieved when a higher amount of polymer chain segments is intercalated.As expected, for each value of gallery thickness, the coefficient of diffusion decreases with increasing volume fraction of the nanofiller.The predictions according to the Nazarenko model, (11), also reported in Figure 9, are not able to capture the dependence of the coefficient of diffusion on the gallery thickness.In contrast, model prediction according to (13), accounting for intrastack diffusion, shows an excellent agreement with FE simulation results, as confirmed by the coefficient of determination,  2 = 0.972.
The normalized coefficient of diffusion is reported in Figure 10 as a function of the number of platelets in each stack, at different values of the volume fractions, keeping constant   /  = 2.5.
Even in this case, in order to account for the effect of the wide error bars, ANOVA analysis was performed, considering the number of lamellae in the stacks as the analysis factor.An  value equal to 130, 56, and 32 was found for the three different volume fractions, which again indicates that the differences between simulations run at different  11) Model eq. ( 9) Simulation Model eq. ( 11) Model eq. ( 9) Simulation Model eq. ( 11) Model eq. ( 9) number of lamellae in each stack are relevant at the 5% significance level.The simulation results are in good agreement with the expected behavior of nanocomposites, since a decrease of the number of lamellae in each stack, which is equivalent to an increase in the aspect ratio and hence to an improvement of the degree of dispersion, involves a decrease of the diffusion coefficient of nanocomposite.In this case, the model of (11), despite being able to capture the dependence of the coefficient of diffusion on the number of platelets, provides a very poor fitting of the simulation data.In contrast, model of ( 13) provides an excellent fitting of simulation data, being  2 = 0.988.The data in Figure 11 show the normalized coefficient of diffusion as a function of the volume fraction of the stacks, given by (10), for constant values of the aspect ratio of the stacks, given by (7).Therefore, the Fricke model, ( 2)-( 3), was used in order to calculate the normalized intrastack diffusivity,   /  .According to the Fricke model,   represents the diffusivity of the stacks, or, more correctly, the equivalent coefficient of diffusion of a homogeneous medium exhibiting the same transport properties of the stack.The excellent fitting reported in Figure 11 by the continuous line is evidenced by the value of  2 = 0.961.
The values calculated for   /  is reported in Figure 12 as a function of the aspect ratio of the stacks.As it can be observed, the intrastack diffusivity increases with decreasing aspect ratio of the stacks.In order to explain this behavior, the domain surrounding each stack can be considered, as that reported in Figure 13(a).The elliptic cylinder is tangent to the stack and is therefore characterized by the dimensions reported in Figure 13(a).The angle  is the angle formed between the direction normal to stack surface and  the direction of diffusion.For a random orientation of stacks, as that investigated in this work, in 3D  = cos −1 (1/ √ 3) and in 2D  = /4.The 2D visualization of the domain is reported in Figure 13(b).In the domains of Figures 13(a) and 13(b), diffusion can only occur inside the stacks.Therefore, due to the impermeability of the single platelets, it can be assumed that, in  1 and  2 , whereas in the volume  3 , the average path length of a diffusing particle is [10]  Therefore, due to a concentration difference between the two horizontal surfaces Δ, the average gradient in the direction parallel to the platelet surface, , is given as and, with  being the angle between the  and the  direction, the concentration gradient in the  direction is The flow rate is given by the concentration gradient multiplied by diffusivity and by the area crossed by diffusing particles, which is (  − 1)  : The same domain of Figure 13(a), made of stacks of impermeable platelets, can be schematized as reported in Figure 14(a), which is made of a single permeable inclusion, which has the same geometric characteristic (diameter and total thickness) of the stack and has a diffusivity equal to   .In such situation, for   ≪   , the diffusion paths follow parallel to the stack thickness.Further, the concentration gradient in the matrix is negligible compared to the concentration gradient in the stack.This is shown in Figure 14(b), where the concentration gradient and the diffusion path are shown in the corresponding 2D simulation.Therefore, the upper surface of the stack can be assumed to be at a constant concentration, equal to that of the upper boundary of the domain, whereas the lower surface of the stack can be assumed to be at a constant concentration, equal to the concentration of the lower boundary of the domain.Therefore, for a 3D problem, the concentration gradient inside the stack is parallel to the direction  (perpendicular to lamellar stack surface) and is given as The gradient in the  direction is therefore given as And the total flow rate in the  direction can be obtained by multiplying the gradient by the diffusivity of the stack and the flow passage,  2 /4: By combining (28) and (31), The values of stack diffusivity calculated according to (32) are also reported in Figure 12, showing a very good agreement with the data calculated from FE simulation.According to (32), the diffusivity of the stack is mainly dependent on the aspect ratio of the stack itself.For a fixed value of the aspect ratio of the stack, the diffusivity also depends on the number of platelets and the thickness of the galleries.

Conclusions
In this work, a finite elements' model was developed, in order to simulate the 3D diffusion into polymer nanocomposites with nanofillers made of lamellar stacks characterized by permeable galleries, as can occur when organic modified montmorillonites are dispersed in a polymer matrix.The developed method is able to introduce a random distribution of nonoverlapping stacks.Each stack is composed of an arbitrary number of platelets, separated by galleries of arbitrary thickness.To the best of our knowledge, this is the first time that this dual scale approach has been used to simulate the diffusion in polymer nanocomposites, occurring either between lamellar stacks or within lamellar stacks.The results obtained from numerical simulation show that intrastack diffusion, that is, diffusion inside galleries, can be quite relevant, particularly at the higher thickness of the galleries.Therefore, the assumption usually made in modeling and simulating diffusion in polymer nanocomposites, based on the impermeability of the nanofiller, cannot be considered valid.
As a consequence of this, the literature models, which do not account for the intrastack diffusivity, are not able to correctly predict the coefficient of diffusion of intercalated nanocomposites.In contrast, literature models (the proposed one), which account for intrastack diffusion, show an excellent agreement with the simulation data.
The dependence of the coefficient of diffusion on the volume fraction of the stack has been used to estimate an equivalent coefficient of diffusion of a homogeneous medium exhibiting the same transport properties of the stack.This equivalent diffusivity was shown to be significantly dependent on the aspect ratio of the stack.An analytical model which is able to represent such property of the stack was proposed.This model shows a very good agreement with the data obtained from simulation.

Figure 1 :
Figure 1: Sketches for modeling of diffusion in multiphase systems: (a) permeable spherical particles and (b) impermeable high aspect ratio platelets.

Figure 3 :
Figure 3: Scheme for modeling of diffusion in intercalated polymer nanocomposites: (a) impermeable stacks and (b) permeable stacks.

Figure 4 :
Figure 4: FEM simulation domain in 3D for modeling of diffusion in intercalated nanocomposites.

Figure 5 :
Figure 5: Procedure used for calculation of the path length of diffusing particles.

Figure 7 :Figure 8 :
Figure 7: Normalized coefficient of diffusion as a function of gallery thickness at V  = 0.015 and   = 3 calculated form 3D simulations.

Figure 9 :
Figure 9: Normalized coefficient of diffusion as a function of gallery thickness at different V  and   = 3.

Figure 10 :
Figure 10: Normalized coefficient of diffusion as a function of the number of platelets at different V  and   /  = 2.5.

Figure 11 :
Figure 11: Normalized coefficient of diffusion as a function of volume fraction of stacks at different values of the aspect ratio of the stacks and Fricke model fitting.

Figure 12 :
Figure 12: Coefficient of diffusion of the stacks as a function of the stack aspect ratio.

Figure 13 :
Figure 13: Schemes for prediction of diffusive flow rate in the presence of stack with permeable galleries (a) 3D and (b) 2D.

Figure 14 :
Figure 14: Schemes for prediction of diffusive flow rate in the presence of stacks regarded as homogenous inclusions and diffusivity   (a) 3D and (b) 2D solution.