Towards better understanding the stiffness of nanocomposites via parametric study of an analytical model modeling parameters and experiments

Abstract

The stiffness of polymeric materials can be improved dramatically with the addition of nanoparticles. In theory, as the nanoparticle loading in the polymer increases, the nanocomposite becomes stiffer; however, experiments suggest that little or no stiffness improvement is observed beyond an optimal nanoparticle loading. The mismatch between the theoretical and experimental findings, particularly at high particle loadings, needs to be understood for the effective use of nanoparticles. In this respect, we have recently developed an analytical model to close the gap in the literature and predict elastic modulus of nanocomposites. The model is based on a three-phase Mori-Tanaka model coupled with the Monte-Carlo method, and satisfactorily captures the experimental results, even at high nanoparticle loadings. The developed model can also be used to study the effects of agglomeration in nanocomposites. In this paper, we use this model to study the effects of agglomeration and related model parameters on the stiffness of nanocomposites. In particular, the effects of particle orientation, critical distance, dispersion state and agglomerate property, and particle aspect ratio are investigated to demonstrate capabilities of the model and to observe how optimal particle loading changes with respect these parameters. The study shows that the critical distance defining agglomerates and the properties of agglomerates are the key design parameters at high particle loadings. These two parameters rule the optimal elastic modulus with respect to particle loading. The findings will allow researchers to form design curves and successfully predict the elastic moduli of nanocomposites without the exhaustive experimental undertakings.

Introduction

Nanoparticles have great potential to enhance the mechanical properties of polymeric materials. Various nanoparticles are shown to improve the stiffness of polymers: nano clay,^1,2 carbon nanotubes,^3,4 cellulose nanocrystals (CNC),^5,6 layered aluminosilicates. ⁷ This improvement is mostly due to the nanoparticles’ high surface-to-volume ratio,^8,9 and extraordinary properties. It is also shown that the addition of nanoparticles can change the polymer crystallinity ¹ and initiate specific interactions between polymer chains and nanoparticles, ¹⁰ also known as the polymer-particle interface.

The volume of the polymer-particle interaction is maximized when nanoparticles are uniformly dispersed; however, nanoparticles tend to agglomerate due to their high surface area and energy. Particularly, high nanoparticle loadings could result in agglomeration ¹¹ and cause overlaps of the polymer-particle interfaces. Beyond a certain particle concentration, agglomeration is shown to be a limiting or even detrimental factor to the targeted properties of polymers in many studies.^12–15 This adverse effect could cause mismatches between experimental findings and model predictions. Therefore, producing polymer nanocomposites with accurately predicted properties is still a challenge.^16,17

Predicting nanocomposites’ elastic properties at various particle loading levels is an important aspect of nanocomposite studies. Accurate predictions increase the efficient use of nanoparticles and lower the cost and time spent on the experimental work. In recent years, many numerical and analytical models were developed to predict the elastic modulus of nanocomposites.^18–23 Numerical models, especially molecular-level simulations, such as molecular dynamics, may produce accurate predictions in small length scales; however, they require high computational power to achieve macro-scale responses due to complex interactions between simulated elements. On the other hand, analytical micromechanical models may serve as a viable alternative for predicting the elastic modulus of polymer nanocomposites. They are, in general, easy to use, low cost, and reasonably accurate.^6,24

Short-fiber micromechanical models such as the Mori-Tanaka, ²⁵ have been used or modified to predict nanocomposites’ behaviour. For example, Jinsu et al. used the Mori-Tanaka micromechanics to predict elastic moduli of the silane functionalized ceramic nanocomposite. ²⁶ It was shown that the Mori-Tanaka model agrees well with the experimental data of 0.04 volume fraction of TiO₂ in acrylate matrix. ²⁶ The Halpin-Tsai micromechanics is another short-fiber composite model ²⁷ that was modified and implemented to estimate material properties of nanocomposites. Zhang et al. ⁴ included carbon nanotube distribution, waviness, and networks for accurate predictions. Their model agrees well with two existing literature data; however, they suggested that a more compressive model capable of considering carbon nanotube agglomeration is needed. Researchers also introduced new variables to micromechanical models or used multiscale approach for more accurate predictions. Arash et al. ²⁸ developed an effective interface model that considers the interface between the reinforcing element and matrix to predict elastic modulus of carbon nanotube polymer composites. The properties of the interfacial region were obtained based on molecular simulations, and these properties were utilized in a modified Mori-Tanaka model. Kim et al. ²⁹ investigated the mechanical properties of carbon nanotube modified carbon fiber reinforced epoxy composites by developing a multiscale composite model. The Halpin-Tsai model was applied to obtain carbon nanotube/Resin 2-phase system and then woven fiber micromechanics was used for the integration of fibers into the model. The model overestimates the experimental results. Kim et al. ²⁹ claimed that the assumption of perfect bonding increases the discrepancy between modeling and experimental results.

Many investigators included agglomeration as another variable into analytical micromechanical models for better predictions. Luo and Daniel ³⁰ proposed a three-phase analytical model based on a hybrid use of the Mori-Tanaka model and rule of mixtures to predict the elastic modulus of polymer clay nanocomposites. The developed model incorporates random orientation of clay layers and various exfoliation levels to capture the mechanical responses. In their experimental work, only low concentration levels (up to 1 w%) of clays were utilized, and higher concentration levels were not explored. The predictions agree well with the experimental findings when experimental parameters are accurately implemented. Shi et al. ⁴ studied the CNC waviness and agglomeration and their effect on the mechanical responses of the composites. The authors employed Eshelby’s inclusion model to predict the elastic modulus of the composite with a spherical inclusion/agglomeration assumption; however, they withheld the exact definition of agglomerates. Villoria and Miravete ³¹ focused on the effects of clusters on the stiffness of composites. They applied the Reuss model to predict the stiffness of clusters because it was assumed that fibers would stick side by side in agglomerated regions. Their study predicts composite stiffness as a function of clustered fibers; however, agglomerated regions’ definition needs further expansion.

Many researchers solely investigated the clustering/agglomeration effect on the mechanical properties of composites from different perspectives. Bhattacharyya and Lagoudas ³² studied different possible microstructures to calculate the clustering effect on the stiffness of composites. In the study, the clustering state is represented by probability distribution functions and a sensitivity study is conducted with respect to different clustering states. Although it is a great approach to studying the clustering effect, Bhattacharyya and Lagoudas did not discuss the exact definition of agglomerates in their study and the elastic modulus. Segurado et al. ³³ studied the effect of particle clustering in finite element analysis by dispersing particles randomly in the computer setting. They idealized the inhomogeneous structure with clustered particles in a spherical manner and created this structure with random dispersion. The dispersion state was limited to random dispersion and only 15 v% particle loading was studied. Hammerand et al. ³⁴ also used the finite element method to calculate the elastic modulus of carbon nanotube embedded epoxy for different cases of clustering. In this study, clustering is defined as a deviation state from uniformly dispersed particles; however, the exact definition of a cluster is not presented, and RVE is limited to only a few particles. Further, they only studied low particle loadings that is up to 1.0 v%.

Attempts to predict elastic properties of nanocomposites are invaluable; however, current studies in literature suggest that there is a gap in knowledge as the available analytical models are, in general, limited in their ability to capture experimental trends at high nanoparticle loadings. Furthermore, exact definition of an agglomerate is missing in the current literature. Intuition and existing analytical model predictions suggest that improvement in stiffness increases as the particle loading increases. On the other hand, the reinforcement is limited to optimal particle loading. Additional particle loadings higher than the optimal particle loading become detrimental to the mechanical properties of nanocomposites^35–39. It can be claimed that the existing analytical models do not account for the complex nature of nanocomposites and agglomeration, and thus they fail to predict the elastic modulus of nanocomposites at high particle loadings.

Recently, we developed an analytical model that uses a three-phase Mori-Tanaka model and the Monte-Carlo method to predict the stiffness of nanocomposites. ⁴⁰ The model predictions and experimental findings match well. A comprehensive study of the model’s parameters can allow us to examine and understand the effect of agglomeration on nanocomposites to address the aforementioned gap in the literature.

Consequently, in this paper, we focus on and study the effect of the model’s parameters to establish a sensitivity analysis for the model. We investigate the effect of critical design variables defined in the model, such as the critical distance, agglomerates’ properties, aspect ratio, particle loading, and various dispersion states of nanoparticles. The predictions of the proposed model are cross-examined with experimental results from the previous study ⁴⁰ where polyamide 6 (PA6) is reinforced with cellulose nanocrystals (CNC).

Modelling

The current micromechanical models^41,42 assume uniformly dispersed particles; however, obtaining uniformly dispersed nanoparticles is unrealistic because of nanoparticles’ high surface energy and area. This non-uniform dispersion state of nanoparticles could result in mismatches between model predictions and experimental results, particularly at high particle loadings. We aim to lower these mismatches by introducing an agglomerate phase into the model. The model, in this study, combines the Monte-Carlo approach with analytical models to capture the effect of agglomeration on the elastic modulus of nanocomposites.

The previously developed model ⁴⁰ aims to offer a simple yet accurate method to predict elastic properties of nanocomposites with respect to particle loading. Its mainframe is introduced here for ease of following; however, interested readers can see the details in the original publication. ⁴⁰

The proposed model is developed based on six main steps: (1) particles are dispersed in a computational setting, (2) agglomerated particles are detected using a machine learning method, (3) agglomerates are averaged based on their volume fractions and particles concentrations, (4) the elastic modulus of the averaged agglomerate is calculated using analytical models, (5) a three-phase Mori-Tanaka model is applied for homogenization, and (6) the Monte-Carlo method is utilized for obtaining statistical information about the elastic modulus of composites. Each step has various parameters, and these parameters are explained in the subsections of this section. For the visualization purpose, the schematics of the homogenization approach are given in Figure 1.

Figure 1.

Schematics of the homogenization approach.

We assume that nanocomposites can be modelled using a three-phase Mori-Tanaka approach where the phases are the particle phase, the agglomerate phase, and the matrix phase. The application of this model requires the shape and aspect ratio of particles, as well as the elastic properties and volume fractions of the phases. The shape and aspect ratio of particles are assigned based on literature data; on the other hand, the elastic properties and volume fractions of the phases are determined using analytical and computational tools. The phases’ volume fractions are calculated based on a computational approach that we developed. A 2-dimensional space (1000 × 1000 nm) in a computational environment is set up for particle dispersion using MATLAB software. The number of particles is calculated based on a particle loading, and particles are dispersed in the two-dimensional space using the Mersenne-Twister algorithm, a pseudorandom number generator. Various dispersion states are covered by changing value of dispersion parameter, $\mu \left[d\right]$ value (detailed in Particle Dispersion section).

The dispersed particles are classified either as agglomerated particles or free particles using the hierarchical clustering method, a machine learning method. ⁴³ The implemented hierarchical clustering algorithm measures Euclidian Distances between particles to detect agglomerated particles. When the Euclidian Distance between two particles is shorter than the value of the critical distance parameter, $\gamma \left[d\right]$ (detailed in Critical distance), particles are considered to be agglomerated. In our approach, particles that are close enough (i.e. within the order of 10 nm) to each other are counted as agglomerated particles and matrix material can exist between them. These regions, i.e. agglomerated particles together with the matrix that is between agglomerated particles, are defined as agglomerated regions. A spherical border is drawn around agglomerated particles to define an agglomerate and to calculate particle concentration within each agglomerate. Once the particle concentrations and volume of each agglomerate are known, a volume-based averaging process is applied to obtain an averaged agglomerate. The final averaged agglomerate is taken as the agglomerate phase.

The volume fractions of the particle and matrix phases are calculated based on the detected agglomerated regions. The particle phase represents the free particles (non-agglomerated particles). Because the total number of agglomerated particles and initial particles are known, the number of free particles is easily calculated. The total volume of free particles is found by multiplying the volume of a particle by the number of free particles. In the end, the particle phase’s volume fraction is obtained by dividing the volume of free particles by the composite volume. Lastly, the matrix phase’s volume fraction is calculated by subtracting the volume fractions of the agglomerate and particle phases from one.

Once the volume fractions are calculated, the elastic properties of each phase are required to utilize the three-phase Mori-Tanaka method. The matrix material is tested for its properties whereas the particle phase’s elastic properties are retrieved from the literature because testing the nanoparticle’s properties is considerably challenging and requires another major study. For the agglomerate phase’s properties, we used two micromechanical models: Halpin-Tsai and Reuss model (detailed in section 2.5).

The calculated volume fractions and elastic properties of the phases are inserted into the Mori-Tanaka model to homogenize the composite system and predict the elastic modulus of the nanocomposite. Main modelling equations are presented here. Standard notation is used throughout the manuscript. The notations , t, T, and T are vector, scalar values, and second-order tensor, respectively. Subscript letters represent the phases of the system, such as m for matrix, p for particle, and a for agglomerate phase. The closed-form of the three-phase Mori-Tanaka model is given by equation (1) ⁴⁴

$\mathit{C}=\left(v_m{\mathit{C}}_{\mathit{m}}+v_f\left\{{\mathit{C}}_{\mathit{p}}{\mathit{A}}_{\mathit{p}}\right\}+v_a{\{\mathit{C}}_{\mathit{a}}{\mathit{A}}_{\mathit{a}}\}\right){\left(v_m\mathit{I}+v_p{\mathit{A}}_{\mathit{p}}+v_a{\mathit{A}}_{\mathit{a}}\right)}^{-1}$

(1)where $\mathit{C}$ is the stiffness tensor of the composite, ${\mathit{C}}_{\mathit{m}}$ is the stiffness tensor of the matrix, ${\mathit{C}}_{\mathit{p}}$ is the stiffness tensor of the particle, ${\mathit{C}}_{\mathit{a}}$ is the stiffness tensor of the agglomerate, $\mathit{I}$ is the identity tensor, ${\mathit{A}}_{\mathit{p}}$ is the strain concentration tensor of the particle, ${\mathit{A}}_a$ is the strain concentration tensor of the agglomerate, $v_m$ is the volume fraction of the matrix, $v_p$ is the volume fraction of the particles and $v_a$ is the volume fraction of the agglomerate, and the curly brackets “{}” stand for the indication of orientation averaging (detailed in section 4.3). The strain concentration tensor of $phase\hspace{0.17em}i$ is given by equation (2)

${\mathit{A}}_i=\left[\mathit{I}+{\mathit{S}}_{\mathit{i}}{\left({\mathit{C}}_{\mathit{m}}\right)}^{-1}{\left({{\mathit{C}}_{\mathit{i}}-\mathit{C}}_{\mathit{m}}\right)}^{-1}\right]$

(2)where ${\mathit{S}}_i$ is the Eshelby tensor of $phase\hspace{0.17em}i$ and its closed form can be found in Mura’s book. ⁴⁵

In the last step of the model, the Monte-Carlo method is utilized to acquire comprehensive data about the elastic modulus of nanocomposites. The Monte-Carlo method uses repeated random sampling to predict an outcome range of an uncertain event or problem. ⁴⁶ The uncertain problem, in this study, is the lack of knowledge about the exact dispersion state of particles. The dispersion state is the main foundation of stiffness prediction. In our model, random sampling corresponds to determining the particles’ locations randomly in the defined 1000 × 1000 nm² space. As soon as particles are randomly dispersed, the aforementioned modelling steps are completed to calculate the volume fractions of the three phases. Based on the calculated volume fractions, the elastic modulus of the composite is found using the three-phase Mori-Tanaka model. For accurate and reliable predictions, this process is repeated one hundred times for each particle loading. We kept the repetitions at one hundred to ensure the timely completion of Monte-Carlo and reach the near-infinite number of possible predictions. The results for each particle loading are illustrated in boxplots.

Particle dispersion

In the developed model, we introduce a parameter $\mu \left[d\right]$ that controls the dispersion state of the particles in the computational setting. Experimentally, the dispersion state is usually quantified by measuring the distance between neighbouring particles that are observed on microscope images.^47–49 The investigators^47–49 measure and record neighbouring distances and present them as histogram diagrams. The obtained histograms are fitted into lognormal probability density functions to discuss the dispersion state of particles. Here, instead of evaluating microscope images to obtain lognormal distributions, we utilize lognormal probability density functions to generate synthetic dispersion states. A lognormal probability density function is formed based on two variables: the mean value ( $\mu$ ) and standard deviation ( $\sigma$ ) of the variable’s natural logarithm.^48,50 The dispersion parameter $\mu \left[d\right]$ is defined as the coefficient of particle diameter. The value of the dispersion parameter $\mu \left[d\right]$ is multiplied by the particle diameter to find the mean of logarithmic values ( $\mu$ ) that are used to establish lognormal probability density functions.

Particles’ locations are determined based on random numbers that are generated from the established lognormal distribution. The location of the first particle is selected from the uniform distribution. For the second particle, a random angle and a random distance are chosen from the uniform probability distribution and lognormal probability distribution, respectively. The chosen random values for the angle and distance are utilized to locate the second particle with respect to the location of the first particle. Each particle is located in the same manner – the angle and distance are selected to locate the new particle’s location with respect to the previous particle’s location. The dispersion of particles is continued until all particles are located. It is important to note that particles cannot occupy the same location in the model. Thus, in the case of particle intersection, a new location is assigned to the newest particle based on the random location selection from all possible locations.

In addition to the lognormal distribution, a uniform probability distribution is used to disperse particles. In the case of uniform distributions, the constant probability distribution function is used to assign the location of a particle. Once a particle is located, the next one is located using the uniform distribution. In the case of an overlap, a new location is assigned to the newest particle. Here, we aim to capture various dispersion states, from uniformly dispersed particles to agglomerated particles. Although different states of dispersion are established, agglomerates are needed to be defined and understood well to investigate their effects on the stiffness of nanocomposites.

Critical distance

The critical distance parameter $\gamma \left[d\right]$ is introduced to differentiate agglomerated particles from non-agglomerated particles as the distinction between them has been vague in the composite literature. Existing agglomerate quantification studies analyze and compare different microscope images to show relative agglomeration states of nanoparticles instead of focusing on individual agglomerates. ⁵¹ The distinction between agglomerated and non-agglomerated particles can only be accommodated if individual agglomerates are defined and quantified. Here, the critical distance parameter $\gamma \left[d\right]$ is proposed to detect agglomerated particles, define an agglomerate, and differentiate one agglomerate from another. The critical distance parameter is inspired by the “cut-off” argument defined in the “cluster” function that is used in the hierarchical clustering method in MATLAB.

The hierarchical clustering, a machine learning method, is used to group data sets into a cluster tree where the tree represents a hierarchy of clusters. Two hierarchical clustering methods exist: agglomerative (bottom-up approach) and divisive (top-down approach). In this study, we utilize the agglomerative hierarchical clustering method because it is a built-in tool in MATLAB Statistics and Machine Learning Toolbox, and it is easier to comprehend intuitively. In the agglomerative hierarchical clustering method, particles start with their own cluster and then combine into bigger clusters based on the distance between them. This method can be described in three main steps: (1) the distance between particles is calculated using a “distance” metric, (2) particles are linked and grouped with a “linkage” method, and (3) the number of clusters (agglomerates) is determined using the value of “cutoff” argument.

In this study, the agglomerative hierarchical clustering method is adopted as a solution to detect agglomerates based on the location data set. Firstly, the distance between particles is calculated using the “Euclidean distance” metric. The “Euclidean distance” is the length of a line between two points that is calculated using the Pythagorean theorem. Secondly, the particles are linked together based on the “single” linkage method. The “single” method takes the shortest distance between particles to link and group them. Herein, the “single” method is used because we assume that particles that are closest to each other should belong to the same agglomerate. Thirdly, each agglomerate is classified, and the number of agglomerates is found using the critical distance parameter. The critical distance parameter is used as the “Cut-off” argument that determines the number of agglomerates in the system. If the critical distance is higher than the shortest distance between two particles or agglomerates, they are grouped together and become members of the same agglomerate. The critical distance parameter is defined as a function of the particle diameter and the value assigned to $\gamma \left[d\right]$ represents the coefficient of the particle diameter. For example, $\gamma \left[d\right]=2$ means that the critical distance is 10 nm when the diameter of the particle is 5 nm. If the center-to-center distance of any particles is less than 10 nm, then they are counted as agglomerated particles, and particles belong to the same agglomerate. Thus, the value of $\gamma \left[d\right]$ is crucial for detecting agglomerated particles and differentiating one agglomerate from another.

The value of $\gamma \left[d\right]$ may depend on the constituents of the composite; however, data in the literature guides us to estimate a range. In a finite element analysis study, from Sheng et al. study, strain fields intersect when the distance between fillers is less than 4 times the diameter of the fillers. In another study, Liu et al. studied nanoparticle dispersion with molecular dynamics and discussed the interaction between nanoparticles. Liu et al. claimed that particles do not have to touch each other to be counted as agglomerates. A short separation of fillers by polymer chains is defined as the state of “local bridging of fillers”. The short separation distance is in the range of 1–2 times the particle diameter. In this study, we follow a similar range (1–4 times particle diameter) to examine the effects of critical distance on the state of agglomerates.

The MATLAB code assigns a number to each agglomerate so that we can detect individual agglomerates and differentiate one from another. The number assigned to the agglomerate is also associated with particles within that agglomerate. For example, if we consider a case where there are two agglomerates, one with 3 particles and the other one with 5 particles. Here, the one with 3 particles is associated with number 1 and the one with 5 particles is associated with number 2. Hence, one can tell which particle belongs to which agglomerate. Once the agglomerates and the particles within each agglomerate are numbered, we establish a spherical border around the agglomerated particles to define agglomerated regions. Figure 2 shows the schematics of an agglomerate. The green filled circles represent the fillers, and the white region between them represents the matrix. Because agglomerated particles do not have to be in contact, an agglomerate can include particles and matrix within its border. The border is fit around the agglomerated particles such that the longest distance between two particles is assumed to be the diameter of the agglomerate. With the knowledge of the agglomerate boundary, the volume fraction of particles within the agglomerate can be calculated and monitored for various critical distance values.

Figure 2.

Schematics of an agglomerate for better understanding of a center-to-center distance.

Particle orientation

Nanocomposites are stiffer in the direction of oriented fibers; however, representing and utilizing the orientation of fibers in modelling may vary and can be challenging. One possible representation of fibers’ orientation is a probability distribution function. According to Advani and Tucker, ⁵² this probability distribution function can be replaced by an orientation tensor for ease of calculation in continuum equations. The orientation of the fiber is required to establish the orientation tensor. Figure 3 describes the orientation of a fiber using spherical coordinate angles ( $\theta ,Ф$ ) under three main assumptions: fibers are rigid cylindrical, fibers’ length and diameter are uniform, and the number of fibers per unit volume is uniform.

Figure 3.

The definitions of $\theta ,Ф,$ and $fiber\hspace{0.17em}\mathit{f}$ in a Cartesian coordinate system.

The probability of having a fiber between angles ${\theta }_1$ and ${\theta }_1+d\theta$ and ${Ф}_1$ and ${Ф}_1+dФ$ is given by Equation (3)

$P\left({\theta }_1\le \theta \le {\theta }_1+d\theta ,\hspace{0.17em}{Ф}_1\le Ф\le {Ф}_1+dФ\right)=\phi \hspace{0.17em}\left({\theta }_1,{Ф}_1\right)\mathit{sin}{\theta }_{1}d\theta dФ$

(3)where $\phi \hspace{0.17em}\left({\theta }_1,{Ф}_1\right)$ is a probability distribution function. The orientation of a fiber can also be represented as a vector $\mathit{f}$ , and the components of the vector $\mathit{f}$ in a Cartesian coordinate system can be written by equations (4), (5), and (6)

$f_1=\mathit{sin}\theta \cos Ф$

(4)

$f_2=\mathit{sin}\theta \mathit{sin}\mathrm{Ф}$

(5)

$f_3=\cos \theta$

(6)Because the orientation can be represented by $\mathit{f}$ , the distribution function can be written as $\phi \left(\mathit{f}\right)$ . If it is assumed that all possible directions of a fiber correspond to the unit sphere and a fiber is a unit vector, the surface integral of the unit sphere can be calculated by equation (7)

${\int }_{\theta =0}^{\pi }{\int }_{\mathrm{Ф}=0}^{2\pi }\phi \hspace{0.17em}\left(\theta ,\mathrm{Ф}\right)\mathit{sin}\theta d\theta d\mathrm{Ф}=\oint \phi \left(\mathit{f}\right)d\mathit{f}=1$

(7)

Advani and Tucker ⁵² showed that the probability distribution function can be represented by tensors, and they provided a widely used orientation tensor that represents the orientation state of fibers. The fourth-order orientation tensor that is defined by dyadic products of the vector is given by equation (8)

$T_{ijkl}=\oint f_i{f}_j{f}_k{f}_l\phi \left(\mathit{f}\right)d\mathit{f}$

(8)

The integration of the product including all possible directions with $\phi$ as the weighting function can be used for the averaging process. The orientation averaging can be represented by $\{\hspace{0.17em}\}$ and expressed by equation (9)

$\overline{{\mathit{T}}_{ijkl}}=\left\{{\mathit{T}}_{ijkl}\right\}={\int }_0^{\pi }{\int }_0^{2\pi }{\mathit{T}}_{ijkl}\left(\theta ,\mathrm{Ф}\right)\sin \theta d\theta d\mathrm{Ф}$

(9)where ${\mathit{T}}_{ijkl}\left(\theta ,Ф\right)$ is the orientation distribution in tensor format, $\overline{{\mathit{T}}_{ijkl}}$ is the averaged orientation tensor and $\{\hspace{0.17em}\}$ represent averaging process. ⁵³

Aspect ratio of particles

The aspect ratio, defined as the length to diameter ratio, is another key design parameter in composite modelling and manufacturing. Some micromechanical models, such as Shear-log, Halpin-Tsai, Mori-Tanaka take the aspect ratio as a variable. In this study, we investigate the effect of the aspect ratio using the Mori-Tanaka model.

In the three-phase Mori-Tanaka model, the aspect ratio of particle phase and agglomerate phase is required for the elastic modulus calculations. The aspect ratio of the agglomerate phase is taken as one because it is assumed that particles agglomerate in a spherical form due to the minimum surface energy requirement and mechanical stability.^4,54 On the other hand, we investigated the effect of the aspect ratio of the free particles on stiffness with three different values: 15, 30 and 60. We selected that range because the aspect ratio of reinforcement used in our experiments, CNC, changes from 10 to 70 depending on the source of CNC. ⁵⁵

Material properties

Material properties of the particle, agglomerate and matrix phases are needed for the model implementation. The matrix (PA6) is tested for its elastic modulus, and the particle’s (CNC) elastic modulus is retrieved from the literature data ⁵⁵ for calculations. The agglomerate phases’ properties are more complicated than matrix and particle phase’s properties, because of the lack of the agglomerate definition and studies in the literature. Although many studies examine the effect of agglomeration on the properties of composites, ⁵¹ there is no consensus about the individual agglomerate’s elastic modulus. Some studies ⁵⁶ assume that agglomerates behave like a void, whereas others calculate their properties based on the Reuss model ⁵⁷ or the Halpin-Tsai model. ⁵⁸ For examples, Khodayari et al. ⁵⁹ used Reuss model to predict elastic modulus of cellulose nanofibrils and generalized the findings for multiphase nanocomposites. Esmizadeh et al. ⁵⁸ used the Reuss model and Halpin-Tsai model to predict elastic modulus of orgonoclay nanocomposite . In this study, either the Reuss model or Halpin-Tsai model are utilized to calculate the elastic modulus of agglomerated regions as they were also applied in various nanocomposite studies. ⁵⁹ These models would provide a good range of possible outcomes, and the effects of stiffness of agglomerates on the stiffness of composite can be observed. The Reuss model is given by equation (10)

${E_a}^{-1}=v_{am}{E_m}^{-1}+v_{ap}{E_p}^{-1}$

(10)where $E_a$ is the elastic modulus of the agglomerate, $E_m$ is the elastic modulus of the matrix, and $E_p$ is the elastic modulus of the particle; $v_{am}$ is the volume fraction of the matrix within the agglomerate, and $v_{ap}$ is the volume fraction of the particles within the agglomerate. Material properties of the particle and matrix for the agglomerate property calculation are kept the same as their individual properties.

The randomly orientated Halpin-Tsai model is calculated based on the particle volume fraction within agglomerated regions by equation (11)

$E_a=\raisebox{1ex}{$5$}\!\left/ \!\raisebox{-1ex}{$8$}\right.E_L+\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$8$}\right.E_T$

(11)where $E_L$ is the longitudinal modulus and $E_T$ is the transversal modulus. The longitudinal modulus $E_L$ is calculated by equation (12) ⁶⁰

$E_L=E_m\frac{1+2\beta v_{ap}{ƞ}_l}{1-{ƞ}_l}$

(12)where β and ${ƞ}_l$ are geometrical parameters. β converges to 2 for transverse modulus or converges to $2\alpha$ for longitudinal modulus, where $\alpha$ is the aspect ratio of the particles defined as ratio of length (l) to diameter (d) of particles. The transversal modulus $E_T$ is calculated by equation (13)

$E_T=E_m\frac{1+2\beta v_{ap}{ƞ}_t}{1-{ƞ}_t}$

(13)where ${ƞ}_t$ is a geometrical parameter for transversal modulus, and ${ƞ}_l$ and ${ƞ}_t$ are given by Eq. (14) and Eq. (15), respectively

${ƞ}_l=\frac{\frac{E_P}{E_m}-1}{\frac{E_p}{E_m}+\beta }$

(14)

${ƞ}_t=\frac{\frac{E_P}{E_m}-1}{\frac{E_p}{E_m}+2}$

(15)

Experimental methodology

The experimental work to validate the model was presented in a previous study. ⁴⁰ Herein, we summarize the procedure, and readers are encouraged to see the previous study for detailed information.

PA6 was dissolved in formic acid, and then CNC was added to the suspension based on the target concentration. The prepared suspension was kept under agitation for approximately an hour. The prepared suspension was sonicated for 45 min before the spinning process. 2 mL of the suspension was cast on a rectangular glass substrate, and the glass substrate was spun at 2000 r/min for 15 s and at 3000 r/min for 30 s.

CNC-PA6 nanocomposite samples were tested using TA Instrument ElectroForce 3200 with 10 N load cells. The elastic modulus of samples was calculated from the linear region of the stress-strain curve and reported for comparison of the model predictions and experimental results.

Model implications and discussion

Effect of critical distance and particle orientation

The elastic modulus of randomly and aligned fiber composites in the longitudinal direction are calculated with the critical distance parameter $\gamma \left[d\right]$ for the values of 1, 2, 3, and 4. The model outputs for the aligned and randomly oriented particles are given in Figure 4(a) and (b), respectively. In the model predictions, the elastic modulus of agglomerates is calculated based on the Reuss model, the aspect ratio is taken as 30 and particles are dispersed uniformly. The matrix material is assumed to be isotropic material with the elastic modulus of 911 MPa according to experimental results retrieved. ⁴⁰ Reinforcing material, CNC, is also assumed to be isotropic material with the elastic modulus of 150 GPa. The predictions are shown in box plots as a hundred elastic modulus values for each particle loading are predicted during the Monte-Carlo method. The red crosses in the graphs show outliners, and the dashed lines are used to connect average values.

Figure 4.

Elastic modulus predictions of cellulose nanocrystals with polyamide 6 (CNC-PA6) composite as a function of CNC loading for (a) aligned and (b) randomly oriented particles at various $\gamma \left[d\right]$ values.

The trends of predictions in Figure 4(a) and (b) are similar at the same values of $\gamma \left[d\right]$ ; however, the trends vary among different $\gamma \left[d\right]$ values. The critical distance parameter $\gamma \left[d\right]$ is used to detect agglomerated and free particles. If the critical distance is higher than the shortest distance between two particles, those particles are counted as agglomerated particles. When $\gamma \left[d\right]=1$ , the model predicts continuous enhancement with respect to particle loading because the three-phase Mori-Tanaka converges to the two-phase Mori-Tanaka: the particle and matrix phases without an agglomerate phase. Selection of $\gamma \left[d\right]=1$ (critical distance of 5 nm) eliminates the chance of agglomerate formation because the closest distance between particles’ centers cannot be lower than 5 nm. As a result, each particle behaves as a free particle, and an approximate linear trend is observed in the case of $\gamma \left[d\right]=1$ . On the other hand, the trend is highly different with the selection of $\gamma \left[d\right]=2$ compared to the selection of $\gamma \left[d\right]=1$ . The predictions increase up to an optimal point of particle loading and then start to drop with particle over loadings. This prediction trend is also valid for $\gamma \left[d\right]=3$ and $\gamma \left[d\right]=4$ , although the predicted values are different. There are 3 main reasons for this type of trend: (1) the high probability of agglomerate formation, (2) the low elastic modulus of agglomerates and (3) the low aspect ratio of agglomerates.

Particles are expected to agglomerate when $\gamma \left[d\right]>1$ because the shortest distance between particles can only be lower than the critical distance at values of $\gamma \left[d\right]>1$ . The probability of agglomerate formation increases with increasing particle loadings and higher $\gamma \left[d\right]$ values. High particle loadings decrease the shortest distance between particles, and high $\gamma \left[d\right]$ values result in higher critical distances. Thus, increase in both parameters leads to the higher probability of agglomeration. The higher probability of agglomeration accelerates the agglomerate phase domination in overall model response, and the properties of agglomerates become major contributors to the model predictions. Because the Reuss model predicts lower than the two-phase Mori-Tanaka model and is used for predicting the elastic modus of agglomerates, the three-phase model predicts lower values after an optimal particle loading. The lowest elastic moduli are observed for $\gamma \left[d\right]$ value of 4 because the critical distance reaches its highest value among the four values. This high value of critical distance yield easier particle agglomeration and low elastic modulus.

Another reason for lower model predictions after the optimal particle loadings at $\gamma \left[d\right]>1$ is the low aspect ratio of agglomerates. We assume that the agglomerate phase is in spherical form and has an aspect ratio of one, which is much lower than the aspect ratio of CNC. This decrease in the aspect ratio, according to the Mori-Tanaka model calculations, also lowers the model predictions. The decrease in the aspect ratio of phases, the probability of agglomerate formation and the low elastic modulus of agglomerates result in differences between the prediction trends of various $\gamma \left[d\right]$ values.

The prediction values for the aligned particles are higher than randomly oriented ones, although the trends are similar at the same value of $\gamma \left[d\right]$ . The dramatic influence of alignment on the elastic modulus can be seen when in Figure 4(a) and (b) are compared at $\gamma \left[d\right]=2$ . According to the model predictions, the maximum elastic modulus for aligned particles is 2800 MPa, whereas it is approximately 1300 MPa for the randomly oriented particles. The particle alignment provides approximately 2.5 times stiffer nanocomposite than randomly oriented particles at 5.0 v% CNC loading. Stiffer nanocomposites in the longitudinal direction are obtained with aligned particle reinforced composite unless particle loadings are not higher than 7.5 v%. Although particle alignment provides stiffer nanocomposites in the longitudinal direction, obtaining aligned particles in nanocomposites could require a specialized manufacturing method such as electrospinning. The electric field can be used to align particles in electrospinning method; however, most of the engineering manufacturing methods such as injection molding could result in random orientation of particles. Manufacturing method should be taken into account while applying the model for the elastic model prediction of nanocomposites.

The impact of the particle concentration and $\gamma \left[d\right]$ value on agglomeration tendency is investigated by extracting the volume fractions of each constituent. Table 1 shows the volume fractions of each constituent at various particle loadings and $\gamma \left[d\right]$ levels. The initial volume fraction of particles ( $v_{ip}$ ) is given in the first column of Table 1. Based on the initial volume fractions, the number of particles is calculated and reported in the third column. As expected, when the initial volume fraction of particles increases from 0.35 to 0.85, the number of particles increases from 1800 to 4372. Based on the selected $\gamma \left[d\right]$ , the number of agglomerates and the volume fractions of each constitute are calculated. It can be seen that agglomerated particles, given in the fourth column, increase as $v_{ip}$ and $\gamma \left[d\right]$ increase. However, the number of agglomerates (# Agg) does not increase in the same ratio, and it even decreases when $v_{ip}=0.85$ .

Table 1.

Volume fractions and number of each constitute at various $\gamma \left[d\right]$ values and partcile $v_p$ %.

$v_{ip}$	$\gamma \left[d\right]$	# Particle	# Agg	# Agg. particles	$v_p$	Ave $v_a$	Ave $v_{ap}$	$v_m$
0.35	2	1800	367	980	0.016	0.084	0.226	0.9
0.35	3	1800	368	1463	0.007	0.21	0.135	0.783
0.85	2	4372	737	3803	0.011	0.392	0.188	0.597
0.85	3	4372	59	4318	0.001	0.579	0.145	0.42

The decrease in the number of agglomerates is related to the particle concentration and size of agglomerates. When $\gamma \left[d\right]$ is set to 2 for 0.85 $v_{ip}$ , high number of small agglomerates are formed (737 agglomerates with 0.392 total average volume of fraction agglomerates). On the other hand, when $\gamma \left[d\right]$ is set to 3 for 0.85 $v_{ip}$ , low number of big agglomerates are formed (59 agglomerates with 0.579 total average volume of fraction agglomerates). Thus, it can be suggested that big agglomerates are formed with less concentrated particles. The volume fractions of particles within agglomerates ( $v_{ap}$ ), in the eighth column, at 0.85 $v_{ip}$ can be compared to observe the effect of $\gamma \left[d\right]$ at values of 2 and 3. Agglomerated particle fraction $\left(v_{ap}\right)$ decreases from 0.188 to 0.145 when $\gamma \left[d\right]$ increases from 2 to 3 for the same initial particle fraction. It again suggests that big agglomerates are formed with less concentrated particles.

The volume fraction of free particles and matrix highly depend on $\gamma \left[d\right]$ . It is clear that free particle fraction ( $v_p$ ) drops with increasing $\gamma \left[d\right]$ for the same initial particle fractions. Similar trend is also observed for the matrix phase fractions. The matrix volume fraction ( $v_m$ ), in the ninth column, decreases as ( $v_{ip}$ ) and $\gamma \left[d\right]$ increases because of higher agglomerated matrix and particle fractions. Because the free particles are the main reinforcing elements, the model is expected to predict lower composite elastic modulus with increasing $\gamma \left[d\right]$ .

Effect of critical distance and agglomerate property

The effect of $\gamma \left[d\right]$ and agglomerate’s properties on the longitudinal elastic modulus of randomly oriented particle nanocomposites is shown in Figure 5. The values of 2 and 4 for the $\gamma \left[d\right]$ are investigated along with two different models, which are applied to calculate agglomerates’ properties: the Reuss and Halpin-Tsai models. While the $\gamma \left[d\right]$ values and properties of agglomerates are changed, the other parameters are fixed in the model calculations. Particles are assumed to be randomly oriented, and the aspect ratio ( $\alpha$ ) is taken as 30. Locations of particles are randomly selected from the uniform probability distribution function.

Figure 5.

The effect of $\gamma \left[d\right]$ and agglomerate’s property on the elastic modulus of randomly particle nanocomposites based on the Reuss and Halpin-Tsai (HT) used for agglomeration.

In the case of the Halpin-Tsai model, the elastic modulus of nanocomposites increases with increasing particle loading regardless of the $\gamma \left[d\right]$ value. The average elastic modulus improves from 1890 MPa to 2700 MPa when the particle loading increases from 8.5 v% to 11.5 v% with $\gamma \left[d\right]=4$ . The slope of this increase is 270 MPa/v%. In a similar manner, the average elastic modulus increases from 965 MPa to 1295 MPa when the particle loading increases from 0.5 v% to 3.5 v% when $\gamma \left[d\right]=4$ . The slope for the second case is 110.0 MPa/v%. It can be claimed agglomeration contributed to the elastic modulus positively with the Halpin-Tsai model. that both slopes are very close to each other. A similar positive slope is also observed for the model predictions when $\gamma \left[d\right]=2$ . The slope of the predictions between 8.5 and 11.5 is 118 MPa/v%, and the slope between 0.5 and 3.5 is 91 MPa/v%. The increase in the slope shows that agglomeration results in more enhancement in the stiffness of the nanocomposite, which confirms higher elastic modulus values when $\gamma \left[d\right]=4$ . This enhancement depends on particles fraction within agglomerates and remaining matrix fraction within composite. For example, the elastic modulus of agglomerates is calculated as 7969 MPa and 5955 MPa in one of the runs at 1.5 and 10.5 v% respectively when $\gamma \left[d\right]=2$ , and the Halpin-Tsai model is applied. The elastic modulus of agglomerates become 3495 MPa and 3411 MPa at 1.5 and 10.5 v% respectively when $\gamma \left[d\right]=4$ . The elastic modulus values decrease with increasing $\gamma \left[d\right]$ because while matrix fraction within agglomerate increases with increasing $\gamma \left[d\right]$ , matrix fraction in the composite decreases.

The prediction trends are quite different in the case of the Reuss model. We first observe an increase in the elastic modulus, and then a decrease regardless of the $\gamma \left[d\right]$ value. The increase is particularly observed at low particle loadings. For example, the average elastic modulus of the composite increases from 980 MPa to 1190 MPa when the particle loading is increased from 0.5 to 3.5 v% when $\gamma \left[d\right]=2$ . We see a similar improvement in the elastic modulus up to 1.5 v% particle loading, in the case of $\gamma \left[d\right]=4$ . The maximum average elastic modulus is obtained at 1.5 v% when $\gamma \left[d\right]=4$ and 7.5 v% when $\gamma \left[d\right]=2$ . The higher than optimal particle loadings result in a decrease in the elastic modulus because of the agglomerates’ properties. The agglomerates’ properties start to dominate the overall response with high particle loadings and high $\gamma \left[d\right]$ values. Because the Reuss model predicts a lower value than the Halpin-Tsai and Mori-Tanaka models, the agglomerate domination lowers the elastic modulus of the composite. This decrease starts earlier with $\gamma \left[d\right]=4$ due to the higher probability of agglomeration. Thus, it can be said that the impact of $\gamma \left[d\right]$ is apparent when the Reuss model is applied because of the low stiffness of agglomerates. As the agglomerate dominates in the composite material, overall response shifts to the Reuss model predictions at high particle loadings.

Predicting the mechanical behaviour of nanocomposites becomes difficult with increasing particle loadings because experimental findings vary substantially. For example, the elastic modulus values^59,58,57,56 reach a plateau value or decrease with particle over loading. Peng et al. ⁶¹ observed a plateau value at 3.5 w% of CNC in PA6 whereas Morelli et al. ⁶² observed that optimal particle loading is between 5 and 10 w% for unmodified CNC. These studies demonstrate that the optimal particle loading, and the amount of increase and decrease depend on experimental settings; however, the proposed model suggests that the combination of agglomerates’ property and $\gamma \left[d\right]$ determines the overall composite reaction at high concentration levels. The trends in the elastic modulus of composites in^37–39 validate that the Reuss model is more appropriate for the predictions of elastic modulus of agglomerates in the proposed model.

The critical distance parameter is the main parameter that controls the agglomerate formation, and higher $\gamma \left[d\right]$ values result in a higher agglomeration chance. In Figure 5, the values of 2 and 4 for $\gamma \left[d\right]$ are examined. The implementation of $\gamma \left[d\right]=4$ in the model dramatically lowers the predictions in the case of the Reuss model due to agglomerate domination. Thus, the appropriate value of $\gamma \left[d\right]$ is crucial for accurate predictions. The exact value of the critical distance $\gamma \left[d\right]$ may depend on the material system or interface chemistry. Sheng et al. ² modelled polymer/clay composite using finite element analysis. The particles are classified as “isolated”, “partially overlapped” and “completely overlapped” according to the distance and load transfer efficiency. The load transfer efficiency starts to decrease as the distance between nanoparticles is approximately below 20 nm, and they define particles that are closer than these values as overlapping particles. The value of 20 nm corresponds to approximately 4 distance/diameter ratio of clay in their study. Therefore, it can be claimed that the value of $\gamma \left[d\right]$ between 1 and 4 is a reasonable range.

Effect of orientation and aspect ratio of the particles

The effect of orientation and aspect ratio of the particles is shown in Figure 6. Uniform probability distribution is utilized for the particle dispersion, and $\gamma \left[d\right]$ is kept at the value of 2. The Reuss model is applied to calculate the modulus of agglomerates in this case.

Figure 6.

Predictions of the longitudinal elastic modulus of CNC-PA6 composites as a function of CNC concentrations for various orientation and aspect ratio ( $\alpha$ ) of the particles.

Figure 6 contains six predictions based on various aspect ratios and particle orientations. The highest three elastic modulus values are from the aligned composites. For a better understanding of the orientation effect, the aligned and randomly oriented CNC with the aspect ratio ( $\alpha$ ) value of 15 can be compared. The model predicts an average elastic modulus of 1750 MPa for aligned particles and 1125 MPa for randomly oriented particles at 7.5 v%, approximately 60% difference. This difference is higher when the aspect ratio is set to a higher value. For example, the model predicts an average elastic modulus of 3800 MPa for aligned particles and 1460 MPa for randomly oriented particles at 7.5 v% with the aspect ratio of 60, approximately 170% difference. This high difference suggests that particle orientation should be considered as a major contributor to the elastic modulus of nanocomposites. According to the 6 model predictions, it can be said that orientation carries more importance than the aspect ratio. The aligned particle composite with the aspect ratio of 15 provides approximately 50% higher elastic modulus than the randomly oriented particle composite with the aspect ratio of 60, which proves the importance of alignment.

The elastic modulus of nanocomposites is improved with an increasing aspect ratio according to the Mori-Tanaka model. This improvement is much higher when particles are aligned. For example, at 6.5 v% CNC, the increase in the elastic modulus is approximately 150% (from 1600 MPa to 3600 MPa) for the aligned particles when the aspect ratio changes from 15 to 60; however, it is approximately 40% (from 1125 MPa to 1460 MPa) for the randomly oriented particles at the same particle loading.

The maximum improvement with respect to neat PA6, ∼330% increase, is predicted to be at 7.5 v% with $\gamma \left[d\right]$ value of 2 and the aspect ratio of 60. However, the average CNC aspect ratio used in this study is approximately 30, and it is difficult to obtain fully aligned CNC reinforced composites. Thus, the maximum increase in elastic modulus would be expected to be lower than the 330% increase.

Effect of dispersion and critical distance

The effect of dispersion is explored by selecting the location of particles from the lognormal probability distribution instead of using the uniform probability distribution. The main reason to choose lognormal distribution is that distances between neighbour particles obey lognormal distribution according to dispersion quantification studies.^63,51,48

The lognormal probability distribution is defined by $\mu$ (mean of logarithmic values) and $\sigma$ (standard deviation of logarithmic values). It is challenging to know exact $\mu$ and $\sigma$ values because they vary in the studies^{68,67,66,65,64,51,48} depending on material system, production, and characterization methods. Instead of implementing an exact value, various $\mu$ values are used to study the dispersion effect on the composite’s modulus. The parameter $\mu \left[d\right]$ is defined as the coefficient of the particle diameter in this section i.e., $\mu \left[d\right]=0.5$ means that the mean of logarithmic values is 2.5. The more uniform particle dispersion is obtained with a higher $\mu \left[d\right]$ value.

The effects of the dispersion $\mu \left[d\right]$ and critical distance parameters $\gamma \left[d\right]$ on the elastic modulus of the composite at different particle loadings (1.0–15.0 w%) are presented with the surface plots in Figure 7. The surface plots are generated with 0.2 and 0.04 intervals of $\mu \left[d\right]\hspace{0.17em}and\hspace{0.17em}\gamma \left[d\right],\hspace{0.17em}\mathrm{r}\mathrm{e}\mathrm{s}\mathrm{p}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{v}\mathrm{e}\mathrm{l}\mathrm{y}.$ In these predictions, the agglomerate’s property is calculated based on the Reuss model, the aspect ratio is taken as 30, and free particles are assumed to be aligned. The black spheres in the figures represent the average experimental results from the previous study. ⁴⁰ The black spheres are located in the figures based on predictions' minimum mean squared error with respect to the obtained experimental data.

Figure 7.

Predictions of elastic modulus of nanocomposites as functions of $\gamma \left[d\right]$ and $\mu \left[d\right]$ at 1.0 (a), 2.5 (b), 5.0 (c), 7.5 (d), 10.0 (e) and 15.0 (f) w% of CNC and corresponding experimental findings from ⁴⁰ shown with black spheres.

The surface plots show that the model can predict experimental findings well for each CNC loading. The model calculates the highest outputs at $\gamma \left[d\right]=1$ because the model assumes no agglomeration in the system, which results in a two-phase Mori-Tanaka model. On the contrary, the model predicts lower elastic modulus values at high $\gamma \left[d\right]$ values and high particle loadings because of agglomeration. The transition between high and low ends is well captured by setting various levels of $\gamma \left[d\right]$ and $\mu \left[d\right]$ at any particle loadings.

Experimental results are relatively high compared to neat PA6 (911 MPa) at 1.0, 2.5 and 5.0 w% CNC and they are captured with low $\gamma \left[d\right]$ and high $\mu \left[d\right]$ . On the other hand, there is a dramatic drop in the elastic modulus values after 5.0 w% CNC loading. The drop and low elastic modulus values are captured with high $\gamma \left[d\right]$ and low $\mu \left[d\right]$ values. The high $\gamma \left[d\right]$ and low $\mu \left[d\right]$ values result in the domination of agglomerates in the model and since the Reuss model is used for agglomeration, the surface plots bend toward the lower end.

Higher values of $\mu \left[d\right]$ in the model generate a better dispersion and stiffer nanocomposites, and many experimental studies^70,69 suggested that well-dispersed particles provide stiffer nanocomposites. Azizi et al. ³⁹ observed an improvement in the elastic modulus to a certain extent and then a decrease in the modulus after optimal particle loading. According to Azizi, poor dispersion of CNC results in the drop of the elastic modulus. Similar to the experimental work, ³⁹ the poor dispersion results in decrease in elastic modulus in the model. This drop after an optimal particle loadings, is satisfactorily captured by the developed model. One can observe the effect of dispersion by examining Figure 7(b) at various values of $\mu \left[d\right]$ . The elastic modulus values increase with increasing $\mu \left[d\right]$ values (increasing dispersion) at a constant $\gamma \left[d\right]$ value. It is difficult to know exact value of $\mu \left[d\right]$ because the exact value may depend on the process, materials, and even on particle loading. A detailed study should be conducted to correlate the critical distance and the logarithmic mean of neighbour distance ( $\mu$ ) to the variables in the experiments. Each material system may require particular modeling inputs; however, it is shown here that the modeling outputs can capture a wide range of experimental results.

Another aim of this study is to cross-examine the model predictions with respect to particle loadings. Here, we fine-tune the values of the model parameters depending on the experimental setting. In the experimental study, ⁴⁰ the particles and the matrix were dissolved in a solvent before producing the samples. Various agitation methods were used to disperse particles uniformly in the suspension during the mixing process. Thus, it seems reasonable to use the uniform probability distribution to predict the locations of particles in the model. Based on TEM images of the previous study, ⁴⁰ free particles were mostly aligned, and the average aspect ratio was taken to 30. Here, we assume fully aligned particles with the aspect ratio of 30 for the calculations. Furthermore, based on the modeling predictions in this study and experimental results from the literature, the Reuss model seems acceptable to calculate the agglomerates’ elastic modulus. According to these assumptions, Figure 8 exhibits the predictions of elastic modulus of nanocomposites with respect to CNC particle loading along with some conventional analytical models.

Figure 8.

Predictions of elastic modulus of nanocomposites with respect to CNC concentration and experimental findings adapted from. ⁴⁰

A drop in elastic moduli, and high standard deviations, at higher CNC loadings are observed in Figure 8. This can be related to agglomeration state of the particles and their potential effects on homogeneity of specimen thickness. During the experiments, the agglomerated particles were observed with the naked eye at high CNC loadings. TEM analysis from our previous study ⁴⁰ also shows how agglomerates reach micron-size dimensions at 15.0 w% CNC loading where the thickness of the sample was approximately 5 μm. It can be discussed that some micron-sized agglomerates are partially embedded into the matrix. These defects may lead to variations in results of the experiments and could lower the elastic modulus of the nanocomposite. As can be seen, in high loadings, this could even result in nanocomposite’s elastic modulus to be lower than the matrix’s elastic modulus. Although the nature of the specimens result in variations in experimental results, our model captures the trend of the experimental results relatively well compared to existing conventional models.

It is still important to understand the deviations from experimental findings. In our predictions, we assumed aligned particles with a specific aspect ratio value. However, even if it is reasonable to assume particles are aligned, the complete particle alignment at the nanosize level is challenging to achieve. The alignment might also depend on the particle loading. In addition, particle size is likely to have distribution instead of a certain value; however, an average value is used in this study instead of distribution. These uncertainties are thought to be the main reasons for the difference between the model predictions and experimental results. In addition, the dispersion is assumed to be uniform for all concentrations in these predictions, which may not be completely valid. The particle dispersion may change with increasing particle concentration due to increased viscosity or particle-particle interactions. Kalfus et al. ⁷¹ discussed these physical relations and noted that a stable dispersion could be possible with favorable interactions between particle and matrix. Kalfus et al. ⁷¹ also state that temperature and particle loading affects the particle dispersion. Thus, the effect of temperature and particle loading can be related to dispersion parameter in the model. Moreover, Janjar et al. ⁷² reviews the physical relationships between nanoparticles and properties of polymer nanocomposites and shows how matrix type – amorphous versus crystalline can change mechanical response and dispersion state. The interactions between chains and nanoparticles surfaces based on molecular dynamic data and material type can be used when deciding the value of the critical distance parameter. Although variabilities exist in the model and the experimental setting, it can be said that the model satisfactorily captured the experimental results. The proposed model reflects the experimental trend well whereas common analytical models predict a continues increase with increasing particle loading. Linear relationship between the elastic modulus and particle loadings.

Conclusion

The proposed study focuses on a parametric study of our recently developed model that predicts the elastic modulus of nanocomposites based on the three-phase Mori-Tanaka model: free particles, agglomerated regions, and a matrix. The aim is to study the effect of agglomerations on nanocomposites at a broad range of reinforcement loadings. Since the agglomeration tendency of nanoparticles is a challenging problem to predict, the Monte-Carlo and hierarchical clustering method are proposed to capture the agglomerate formation. Along with agglomerate formation, a systematic study is performed to understand the effect of aspect ratio, critical distance, particle orientation, agglomerate property, and dispersion state of particles.

The critical distance parameter and elastic modulus of agglomerates are the key design parameters at high particle loadings. The critical distance parameter defines the agglomerated region that contains matrix and particles. Higher critical distance values result in easier agglomeration. As agglomerates are formed, the agglomerates’ properties begin to dominate the general response of nanocomposites. In that sense, predicting agglomerates’ properties become crucial for the model. Either the Halpin-Tsai or the Reuss models are assigned to predict agglomerates’ properties and examine agglomerates’ effect on the elastic modulus of nanocomposites. While the Halpin-Tsai overestimates the elastic modulus, the Reuss model provides more reasonable results.

The aspect ratio, orientation and dispersion of particles are also investigated to understand the model’s sensitivity. The drastic effect of aspect ratio and particle orientation on the elastic modulus of nanocomposites is observed at high particle loadings. At any aspect ratio, aligned fillers exhibits shows steeper slope than randomly oriented particles. It is concluded that alignment has more influence on the elastic modulus than the aspect ratio in the range of 15–60 for the case of uniform particle dispersion. Particle dispersion, another key parameter, is studied by locating the particles differently in the computational setting. The more uniform dispersion results in higher elastic modulus of nanocomposites regardless of other parameters. It is also observed that the agglomeration is inevitable at high particles even if particle locations are selected from uniform distribution function.

In this work, the capability of the previously developed model is analyzed with a systematic study. The proposed approach captures the experimental trend of elastic modulus of CNC reinforced PA6 samples relatively well compared to the conventional analytical models. As future of this work, the manufacturing methods and critical distance could be categorized based on the material system, and they could be correlated to the dispersion state of nanoparticles. With the existing data in the literature, we believe these correlations can be established via machine learning techniques.