The Barrier Properties of Flake-Filled Composites with Precise Control of Flake Orientation

Additive manufacturing, especially in the form of 3D printing, offers the exciting possibility of generating heterogeneous articles with precisely controlled internal microstructure. One area in which this feature can be of significant advantage is in diffusion control, specifically in the design and fabrication of microstructures which optimize the rate of transport of a solute to and from a contained fluid. In this work we focus on the use of flakes as diffusion-control agents and study computationally and theoretically the effect of orientation on the barrier properties of flake-filled composites. We conducted over 1500 simulations in two-dimensional, doubly-periodic unit cells each containing up to 3000 individual flake cross-sections which are randomly placed and with their axes forming an angle ( π 2 θ − ) with the direction of macroscopic diffusion. We consider long-flake systems of aspect ratio (α ) 100 and 1000, from the dilute ( 0.01 αφ = ) and into the concentrated ( 40 αφ = ) regime. Based on the rotation properties of the diffusivity tensor, we derive a model which is capable of accurately reproducing all computational results ( 0.01 40 αφ < < and 0 π 2 θ < < ). The model requires as inputs the two principal diffusivities of the composite, normal and parallel to the flake axis. In this respect, we find the models of Lape et al. [1] and Nielsen [2] form an excellent combination. Both our model and our computational data predict that at 0 θ > the quadratic dependence of the Barrier Improvement Factor (BIF) on (αφ ) is lost, with the BIF approaching a plateau at higher values of (αφ ). This plateau is lower as (θ ) increases. We derive analytical estimates of this maximum achievable BIF at each level of misalignment; these are also shown to be in excellent agreement with the computational results. Finally we show that our computational results and model are in agreement with experimental evidence at small values of (θ ). How to cite this paper: Tsiantis, A. and Papathanasiou, T.D. (2017) The Barrier Properties of Flake-Filled Composites with Precise Control of Flake Orientation. Materials Sciences and Applications, 8, 234-246. https://doi.org/10.4236/msa.2017.83016 Received: January 7, 2017 Accepted: March 6, 2017 Published: March 9, 2017 Copyright © 2017 by authors and Scientific Research Publishing Inc. This work is licensed under the Creative Commons Attribution International License (CC BY 4.0). http://creativecommons.org/licenses/by/4.0/

the dilute ( 0.01 αϕ = ) and into the concentrated ( 40 αϕ = ) regime.Based on the rotation properties of the diffusivity tensor, we derive a model which is capable of accurately reproducing all computational results ( 0.01 40 αϕ < < and 0 π 2 θ < < ).The model requires as inputs the two principal diffusivities of the composite, normal and parallel to the flake axis.In this respect, we find the models of Lape et al. [1] and Nielsen [2] form an excellent combination.
Both our model and our computational data predict that at 0 θ > the qua- dratic dependence of the Barrier Improvement Factor (BIF) on ( αϕ ) is lost, with the BIF approaching a plateau at higher values of ( αϕ ).This plateau is lower as ( θ ) increases.We derive analytical estimates of this maximum achievable BIF at each level of misalignment; these are also shown to be in excellent agreement with the computational results.Finally we show that our computational results and model are in agreement with experimental evidence

Introduction
Additive manufacturing, especially in the form of 3D printing, offers the exciting possibility of generating articles with precisely controlled internal microstructure.One area in which this feature can be of significant advantage is in diffusion control, specifically in the design and fabrication of microstructures which allow for minimization of the transport of a solute to/from a contained fluid.
Flake-filled polymeric composites, incorporating mica, glass or metallic flakes have found many uses in this direction, as they offer significant processing and property advantages, namely high dimensional stability and low warpage in molding, uniform in-plane mechanical properties, corrosion protection, sound insulation as well as appearance and color control [3] [4] [5] [6] [7].Micron-sized flakes of inorganic materials such as mica, nano-scale platelets of clay minerals such as hectrite, saponite and montmorillonite and more recently graphene-oxide platelets of aspect ratios well over 1000, have been used for this purpose [8].It has been demonstrated that incorporation of such fillers aligned perpendicular to the direction of macroscopic diffusion can be very effective in increasing the tortuosity of the diffusion path of the diffusing species.When the flakes are randomly placed, as would be the case in a flake composite manufactured from the melt, the predicted improvement in barrier efficiency ranges from being ~( αϕ ) in dilute systems, where ( α ) is the aspect ratio and ( ϕ ) the volume fraction of the flakes, to being ~( ) 2 αϕ in more concentrated disper- sions [1] [9]- [16].
One notable disadvantage of traditional processing methods vis-à-vis flakefilled composites is the fact that flake orientation cannot be precisely controlled.In such operations (extrusion, compression or injection molding, thermoforming and others) flake orientation is achieved due to the propensity of the flakes to orient in accordance to the prevailing flow field-either in the main direction of flow when the flow is shear or transverse to it when the flow is extensional [17].An additional shortcoming of traditional flow-processing routes is the inability to utilize high flake loadings since, in that case, the viscosity of the resulting melt becomes prohibitively high.Given the capability afforded by 3D printing to fully control flake orientation as well as generate articles with flake loadings approaching those at maximum packing, it is desirable to predict the effective diffusion coefficient (or its inverse, the barrier improvement factor, BIF) as an explicit function of the flake orientation angle and for very high, previously untenable, concentrations.
The two main approaches which have been used in the literature to-date for this purpose are (i) an ad-hoc incorporation of orientation metrics in existing models for the BIF [8] [18] and (ii) derivation of BIF models from diffusion path calculations [19] [20] [21].In both cases, the proposed models have been derived for low or very-low flake concentrations and have not been extensively tested in the moderate to high-concentration regime, which will be of importance in any 3D printing application.In addition, by not respecting the rotational properties of the diffusivity tensor, these models are not grounded on a sound theoretical footing.This paper addresses the above issues both computationally and theoretically, by proposing a model based on the two principal diffusivities of a flake composite.We also show that the implications of our theoretical model are fully supported by extensive computational results.

Computational
We carry out steady-state diffusion computations in doubly-periodic unit cells containing up to 3000 individual flake cross-sections.These are added in the domain sequentially, using a Random Sequential Addition (RSA) procedure.
Specifically, at each flake placement attempt, two random numbers are used to assign the coordinates of the flake center while its orientation angle ( θ ) is fixed and the same for all flakes.If, after placement, no overlap with other flakes is detected, the process continues with the next flake, until the desired number of flakes has been placed, or, until no flake can be placed after 50,000 attempts; in this case no geometry is generated.In order to enable subsequent meshing of the computational domain, a minimum allowable distance between flakes is imposed; this is ( 2t ) where ( t ) is the thickness of the flake.Since in this work we have dealt with flakes of high aspect ratio ( 100 α = and 1000 α = ), this minimum distance requirement is deemed reasonable so as to not result in excessive local mesh refinement.In a rectangular unit cell of dimensions (H) and (L) containing (N) flakes of dimensions ( , t α ), the flake area fraction ( ϕ ) is We have looked at systems in which 0.01 40 αϕ ≤ ≤ .At higher values of ( αϕ ) it be- comes impossible to fill the space with non-overlapping flake crossections.This not-withstanding, the present study is to our knowledge the first to investigate systems of such large concentration.
In multi-particle simulations, use of doubly-periodic cells is essential when dealing with elongated particles so as to eliminate artifacts of oriented (or, depleted) layers which would otherwise appear adjacent to cell boundaries [11] [16].
The effect of the RVE dimensions on the computed effective diffusivity is also eliminated when using periodic unit cells.A sample unit cell, showing flakes oriented at an angle θ = 0.8 rad with respect to the horizontal (x) axis (extended slightly outside the limits of the unit cell to show the doubly-periodic geometry) is shown in Figure 1.These meshes are then used by Open Foam to solve the steady-state diffusion equation 2 0 C ∇ = , (C) being the solute concentration, and provide its distribu- tion in the domain of interest.The assumption of an isotropic matrix material is also made.The solution also supplies the value of C ∂ ∂n at each point on the upper (or lower) boundary.As a result, the mass flux along this boundary can be calculated as where the subscript (n) indicates numerically computed value, n is the outward unit vector and (L) is the width of the unit cell.Because of impermeable flakes crossing boundaries, which results in sudden local changes of the flux, care must be taken in performing this integration.In this work, we used adaptive intervals and only accepted values of the integral when these were convergent with refinement.Equating this flux with the one obtained from Fick's law in a macroscopic equivalent cell (whose diffusivity is D eff ) we obtain where C ∆ is the macroscopically imposed concentration difference and (H) the height of the unit cell.These effective diffusivities will be presented and discussed for various values of ( θ ), ( α ) and ( ϕ ) in the following sections.

Results and Discussion
In the following we present the results of a comprehensive computational study of diffusion across doubly-periodic unit cells, each containing up to N = 3000 randomly placed impermeable flakes of rectangular cross-section.In such a system, the orientation of each flake is defined by the orientation angle (θ) formed between the vertical axis (y), which is taken to be the direction of macroscopic diffusion, and the outward normal vector on the flake surface.The horizontal axis is indicated as ( x ).Parametric studies have shown that for N > 200 the ob-tained D eff were practically indistinguishable; this conclusion extended for ( αϕ ) as large as 40; therefore most of our computations have been carried out in RVEs containing 500 flake cross sections.We look at systems ranging from dilute to concentrated and in which the fiber orientation (θ) changes between zero (flake orientation perpendicular to the direction of diffusion) to π 2 (fibers oriented along the direction of macroscopic diffusion).We have carried out computations in unit cells similar to those of Figure 1 for 100 α =

Effect of Flake Misalignment on Effective Diffusivity
Representative results of the distribution of (C), also showing the corresponding flake distributions, are shown in Figure 2 and Figure 3.
We define as D 11 the diffusivity of such a system when θ = 0˚ (all flakes oriented perpendicular to the direction of macroscopic diffusion) and D 22 the diffusivity when θ = 90˚ (all flakes oriented parallel to the direction of diffusion).D 11 and D 22 are the principal values of the two-dimensional diffusivity tensor, D.
The diffusivity tensor ′ D corresponding to a system in which the flakes assume an orientation angle θ (counterclockwise with respect to the x-axis) can be determined according to the relation Hence Therefore, the effective diffusivity of this system in the direction (y) forming  The distribution of flakes is also visible.N = 500.
an angle ( π 2 θ − ) with the axis of the flakes will be ( ) ( ) ( ) We will investigate the use of Equation (2) to determine ( ) eff D θ , provided the principal permeabilities D 11 and D 22 are known.By comparing its predictions to our computational results we will identify which models for D 11 and D 22 give the best agreement with computation.
In the first instance we have compared the computational results for dilute cases ( 0.01 αϕ = and 0.1 αϕ = ) with the predictions of Equation ( 2 ).For progressively higher of ( αϕ ) there is a growing discrepancy.
It is of course possible to use diffusivity models for D 11 and D 22 more suitable for concentrated suspensions.A review and evaluation of available models has been carried out by Chen and Papathanasiou [11].Of the models discussed there, we single out those of Cussler and co-workers [1] [9] mainly because of their relevance to the systems modeled here (randomly placed flakes) as well as due to the small number of adjustable parameters needed in their implementation and their extensive use in the literature.Lape et al. [1] proposed that for diffusion across arrays of unidirectional randomly placed flakes it is ( ) In deriving this model, the tortuosity factor was taken to be 1 + αφ/3 and it was further assumed that the ratio of the areas available for diffusion is Implicit in the above derivation is the assumption that the diffusion paths around a flake form straight lines; therefore it is not unreasonable to treat the factor "3" in the expression above as a geometrical parameter that may be adjusted if so suggested by the data.Since that was found to be the case in analyzing our data, we use the model of Lape et al. [1] in the form: in which ( λ ) is an adjustable geometrical parameter.Another model suitable for concentrated aligned flake systems [9] reads ( ) where ( β ) is also an adjustable geometric factor.The following Figure 4 gives a comparison between the computational results, for flakes with 100 α = and for 1 αϕ = and 10 αϕ = , in unit cells similar to those of Figure 3 and the predic- tions of Equation ( 2), in which D 11 is taken from [1] [9] and D 22 from [2].It is evident that use of models for D 11 more suitable for concentrated systems results in significantly improved predictions of D eff for all (θ).The model of Lape et al. [1] gives an excellent agreement with the computational results for 2.5 λ = even for αϕ as low as 0.01 (especially away from ~π 2 θ ) with a slight adjustment of λ to 2.7 at

αϕ =
, while the model of Cussler et al.Additional comparisons for 1000 α = and higher values of ( αϕ ) are shown in Figure 5, in terms of the BIF.There are 50 computational data-points at each value of ( αϕ ) and those at the same ( θ ) almost completely overlap.This has been shown before [11], namely that spatial randomness has a very small effect on the diffusivity of such systems.
In summary, our computational results and the comparisons presented above have shown that the effective diffusivity D eff of a system of randomly placed flakes oriented at an angle ( π 2 θ − ) with the direction of macroscopic diffusion can be predicted by  2.5 λ = .As explained above, this model is in excellent agreement with the computational data for the entire range of ( αϕ ) and (θ) studied.In addition, we compare the predictions of our model to an experimental result [22] [23], namely that for small values of the misalignment angle (θ) it is in which 0 θ = corresponds to a composite fully aligned normal to the direc- tion of diffusion.If D 11 is the diffusivity of the fully-aligned system, the BIF implied by the above statement will be From Equation (2) it can be seen that the BIF implied by our model, (setting, without loss of generality or relevance, D 22 ~D0 ~1) is As shown in Figure 6, at each value of ( θ ) the predictions of Equation (7)   approach asymptotically those of Equation ( 6) albeit at progressively higher values of D 11 (that is, for more dilute systems) as ( θ ) increases.However, the li- miting behavior of Equation (7) in the concentrated regime (small D 11 ) suggests a qualitatively different behavior for the BIF.Our computational results support this prediction, as will be elaborated upon in the following section.With reference to Figure 6, if the model of Lape et al. [1] is adopted for D 11 , a value of Therefore, our model is consistent with Equation ( 6) well into the semi-concentrated regime, for small misalignment angles.

The Effect of Flake Concentration
In aligned systems, it is known [1] [9] [10] that the BIF scales with ( ) 2 αϕ at higher concentrations, and linearly with ( ) αϕ in the dilute regime.No such definitive information is available when deviations from perfect alignment occur.Figure 7 shows all our computational results for 1000 α = . It is clear that while the quadratic rise with ( αϕ ) is indeed observed in aligned systems ( 0 θ = ), this asymptotic behavior is lost as (θ) increase and the BIF approaches a   5); as in Figure 4 and Figure 5 the agreement between the two is excellent.

Limiting Behavior of the BIF at Very High (αφ).
In light of the excellent agreement between computational results and Equation ( 5) it is possible to use the latter to obtain analytical estimates of the leveling-off values of the BIF at each (θ), by observing that the first term of Equation ( 5) becomes negligible at high ( αϕ ), leaving Figure 8 compares our computational results to the predictions of Equation ( 8) as well as the approach to that limit based on Equation (5).A conclusion is obvious-the quadratic rise of the BIF with ( αϕ ) is lost when 0 θ > .For a mi- salignment as small as 5.7˚ (0.1 rad) the upper limit on the achievable BIF from Equation ( 13) is 104-a three-fold decrease from the theoretical BIF of a perfectly aligned composite with 40 αϕ = and a multi-fold decrease from an aligned composite of even higher ( αϕ ).In fact, for such concentrated systems the departure from the theoretical BIF can be very rapid at small misalignment angles, as can be inferred from Equation ( 8).This we show in Figure 9 in which we plot the predictions of Equation ( 8) along with our computational results for 1000 α = and 40 αϕ = .
The above comments and results are particularly pertinent to high aspect ratio Figure 8.The approach to the BIF limit (as predicted by Equation ( 8), dotted lines) for θ = 0.1, 0.2 and 0.4 (in rad) as well as the predictions of Equation ( 5) (solid lines).Points are computational results.α = 1000.8).The rate of decline in barrier performance with even a slight misalignment is very significant at small (θ), when (αφ) is large.
flakes, such as found in exfoliated nanoclay or graphene composites, for which even at low ( ϕ ) a high ( αϕ ) value can be achieved; in our simulations in which 1000 α = , the maximum αϕ of 40 translates into 4% ϕ = .Evidently, Εquation (8) in that case says that the limiting BIF is only a function of the misalignment angle-and our computations are in complete agreement with this prediction.At higher loadings, Equation (8) predicts that the limiting BIF will increase for larger values of ( αϕ ).

Conclusion
In this study we proposed a model for the Barrier Improvement Factor (BIF) of misaligned flake composites which is valid up to very high flake concentrations, as could be found in composites fabricated by 3D printing.The model requires as inputs the two principal diffusivities of the composite, normal and parallel to the flake axis.In this respect, we find that the models of Lape et al. [1] and Nielsen [2] form an excellent combination.The simple algebraic form of the proposed model makes it usable without recourse to special computing facilities.This model was tested exhaustively by comparing to predictions of 2D computer simulations which included up to 3000 randomly placed but uniformly oriented flake cross-sections in each RVE.Each cross-section forms an angle ( π 2 θ − ) with the direction of macroscopic diffusion.Over 1500 simulations were carried out and upon comparison the model was found in agreement with computational results for all misalignment angles and for values of ( αϕ ) up to 40.Both our model and our computational data predict that at 0 θ > the quadratic depen- dence of the BIF on ( αϕ ) is lost, with the BIF approaching a plateau at higher values of ( αϕ ).This plateau is lower as (θ) increases.We derive analytical esti- The boundary conditions are cyclic on the right and left boundaries, namely the top and bottom boundaries, fixed values of concentration are prescribed.On the surface of each flake we impose 0 C ∂ ∂ = nindicating that the flakes are impermeable.At each level of ( α ) and ( ϕ ) we generate ~10 different realizations, each containing up to 3000 randomly

Figure 2 .
Figure 2. Distribution of concentration in geometries with θ = 0 and
), in which Nielsen's [2], model has been used for D 11 and D 22 , namely 22 shown that predictions of Equation (2) based on Nielsen's model for D 11 and D 22 are close to the computational results only for the very dilute regime ( ~0.01 αϕ

2 θ
The latter model Equation (4) can also be used at lower ( αϕ ) values with proper adjustment of the parameter ( β ); at almost parallel to the direction of diffusion) the numerical results are in very close agreement with Equation (2) for all concentrations.Since at ~π the term containing D 22 dominates, this shows that Nielsen's model for diffusion parallel to the flakes is a reliable one, even for (αφ)as high as 10.

Figure 5 .
Figure 5.Comparison of computational results (points) with predictions of Equation (2) for 20 αϕ = and 40 αϕ = .The legend refers to the model used in place of D 11 .For D 22 Nielsen's model [2] was used.In all cases 1000 α = where a value 11 0.01 D = will give 27 αϕ = .

Figure 6 .
Figure 6.Predictions of Equation (7) (broken lines) showing its asymptotic approach to the experimental result represented by Equation (6) (solid line).Larger values of D 11 correspond to more dilute systems.

Figure 7 .
Figure 7. Summary of computational results (circles) for 1000 α = and θ = 0, 0.1, 0.2, 0.4, 0.6, 0.8, and 1.0, all in rad.We observe the anticipated quadratic rise of the BIF with ( αϕ ) for higher values of ( αϕ ) at 0 θ = and also a progressively lower plateau reached at increasing values of the misalignment angle ( θ ).The predictions of Equation (5) are also shown as solid lines.Total of 1295 data points.