A Study on the Effect of Carbon Nanotubes’ Distribution and Agglomeration in the Free Vibration of Nanocomposite Plates

The present work aimed to characterize the free vibrations’ behaviour of nanocomposite plates obtained by incorporating graded distributions of carbon nanotubes (CNTs) in a polymeric matrix, considering the carbon nanotubes’ agglomeration effect. This effect is known to degrade material properties, therefore being important to predict the consequences it may bring to structures’ mechanical performance. To this purpose, the elastic properties’ estimation is performed according to the two-parameter agglomeration model based on the Eshelby–Mori–Tanaka approach for randomly dispersed nano-inclusions. This approach is implemented in association with the finite element method to determine the natural frequencies and corresponding mode shapes. Three main agglomeration cases were considered, namely, agglomeration absence, complete agglomeration, and partial agglomeration. The results show that the agglomeration effect has a negative impact on the natural frequencies of the plates, regardless the CNTs’ distribution considered. For the corresponding vibrations’ mode shapes, the agglomeration effect was shown in most cases not to have a significant impact, except for two of the cases studied: for a square plate and a rectangular plate with symmetrical and unsymmetrical CNTs’ distribution, respectively. Globally, the results confirm that not accounting for the nanotubes’ agglomeration effect may lead to less accurate elastic properties and less structures’ performance predictions.


Introduction
Carbon nanotubes' (CNTs) remarkable characteristics makes them well suited for acting as a reinforcing phase in composite materials where high strength and low density are required or in applications where superior physical, mechanical, thermal, and electrical properties have to be guaranteed [1,2]. Considering these nanoparticles' potential use in a number of engineering fields, it is thus important to characterize the influence that their distribution and eventual agglomeration may have in the behaviour of structures, namely in their free vibrations' behaviour.
A number of researchers have developed very relevant work in the field of composite materials reinforced with nanoparticles. Among them, one can refer the work developed by Sobhani Aragh et al. [3], where an equivalent continuum model based on the Eshelby-Mori-Tanaka (EMT) approach was employed to estimate the elastic properties of functionally graded carbon nanotube-reinforced C 2020, 6, 79 2 of 32 composite (FG-CNTRC) cylindrical panels to study their free vibration response. These authors found out that CNTs' volume fraction can be used for the management of vibrational behaviour of structures. A similar study was developed by Sobhadi Aragh et al. [4], this time using the third-order shear deformation theory to investigate the dynamic behaviour of FG-CNTRCs cylindrical shells. These authors concluded that symmetric distributions of the CNTs' volume fraction, over the thickness direction, provided a greater improvement of the dynamic behaviour when compared to uniform or unsymmetrical volume fraction distributions. Yas and Heshmati [5] studied the influence of various CNTs distributions and orientations, among other parameters on the dynamic behaviour of a functionally graded nanocomposite Timoshenko beam.
Due to their high aspect ratio, low bending stiffness, and van der Waals forces, CNTs tend to agglomerate into bundles [6,7], with this phenomenon interfering directly with the material properties of a CNTRC (carbon nanotube-reinforced composite). To estimate the CNTRC's properties more realistically, it is necessary to take into account this agglomeration effect. However, in most of the studies including CNTRC, the most common homogenization scheme employed is the extended rule of mixtures (ROM) not taking into account this effect.
To overcome this problem, Shi et al. [8] proposed a two-parameter agglomeration model based on the EMT for CNTRC property estimation with randomly oriented straight CNTs to investigate the effect of nanotube waviness and agglomeration on the mechanical properties of CNTRCs.
Van der Waals forces play an important role on CNT micromechanical models namely, to model the interphase between the resin and the CNT for nanocomposites, for predicting their mechanical properties as proposed in the studies led by Shokrieh [9,10]. There is another important consideration about the van der Waals forces when assessing the free vibration behaviour of MWCNT (multi-walled carbon nanotubes), which is the interaction of the van der Waals forces between two different SWCNT (single-walled carbon nanotube) layers and their effect on the natural frequencies and modes' shapes. The study developed by He et al. [11] indicated that when considering van der Waals forces, the modes' shapes change when varying the radii in MWCNT from very small to very large, and the that these forces have a significant impact in higher natural frequencies. Another study that considers the influence of the van der Waals forces between two different SWCNT layers is the study developed by Strozzi and Pellicano [12] on the linear vibrations of TWCNTs (triple-walled carbon nanotubes), where these interactions are modelled by means of a radius-dependent function that relates the displacements of two adjacent layers.
Since the EMT agglomeration model was proposed, a number of studies have emerged using this technique for property estimation. Among them, the dynamic behaviour of an FG-CNTRC Timoshenko beam resting on a Pasternak foundation was studied under the influence of CNT agglomeration, showing that the natural frequencies of the beam are highly influenced by the agglomeration effect [13]. Another study performed using beam models was developed by Heshmati and Yas [14], where the dynamic behaviour of FG nanocomposite beams was evaluated under the influence of CNT agglomeration, and it was also shown that the agglomeration exerted a significant weakening effect in CNTRC.
Kamarian et al. [15] studied the dynamic behaviour under free vibrations of CNTRC conical shells under the influence of agglomerated CNT using the EMT approach and they concluded that the agglomeration significantly affects the natural frequency of the structure. Tornabene et al. [16] performed a study on the effect of agglomerated CNTs on the natural frequencies of FG-CNTRC doubly curved shells, where several parametric studies were developed varying the CNT volume fraction, the CNTs distributions along the thickness, as well as the agglomeration parameters, proving their influence on the dynamic behaviour of the structure. Tornabene et al. [17] published another study now on the static response of nanocomposite plates and shells reinforced by agglomerated CNTs, where the static response of these structures was shown to be considerably affected by the agglomeration of CNTs.
Eringen's nonlocal elasticity theory was employed to develop bending and buckling analyses of agglomerated CNTRC nanoplates resting on a Pasternak foundation by Daghigh et al. [18]. With this C 2020, 6, 79 3 of 32 study, it was concluded that the presence of the foundation beneath the nanoplates reduced the effects of CNT agglomeration. The same author led another study on the dynamic behaviour of size-and temperature-dependent agglomerated CNTRC nanoplates resting on a visco-Pasternak foundation [19], where it was concluded that ignoring the agglomeration effect causes an overestimation of the elastic properties and the presence of the Pasternak foundation clearly increased the nanoplate natural frequency.
More recently, Bisheh et al. [20] studied the dynamics of wave propagation in agglomerated CNTR piezoelectric composite cylindrical shells. The developed analytical model showed a capacity to find the effects of nanotube agglomeration on wave propagation characteristics of these smart nanocomposite shells for the applications of structural health monitoring and energy harvesting. The performance of CNTR porous composite plates bonded with piezoceramic material in the lower and upper faces under static mechanical and electrical loading was evaluated by Moradi-Dastjerdi et al. [21], concluding that embedding pores up to around 90% of the volume of the middle plate leads to an increase of around 11.3% and 2.2% in mechanical and electrical deflections, respectively. The inclusion of CNTs in the middle plate was demonstrated to be beneficial; however, this was to limited values of the CNT volume fraction due to the agglomeration effect.
A modelling paradigm for improving the piezoelectric performance for piezoelectric matrix-inclusion composites based on lead-free ceramics through improved matrices and optimal polycrystallinity in the piezoelectric inclusions was developed by Krishnaswamy et al. [22]. These authors demonstrated for this case that CNT agglomeration near nanotube percolation and the clustering of nanotubes can lead to better matrix hardening and higher permittivity, leading to an improvement exceeding 30% in piezoelectric response.
As mentioned, CNTs have remarkable characteristics that are important for many engineering applications; however, CNTs tend to agglomerate even for low volume fraction's distributions [23]. Hence, assessing the influence of the agglomeration effect on the free vibration behaviour of nanocomposites is fundamental to better design structures built from these materials. The present work aimed to characterize the free vibration behaviour of FG-CNTRC square and rectangular plates concerning the effect of CNTs' distributions and agglomeration, regarding its influence on the natural frequencies and modes' shapes of these plate-type structures.

Materials and Methods
The present section is subdivided into two sub-sections, wherein the most relevant aspects, techniques, and methodologies that support the present study are summarized.

Properties' Estimation through the Eshelby-Mori-Tanaka Approach
In FG-CNTs, CNTs tend to agglomerate in a polymer matrix, due to its high aspect ratio, low bending stiffness (small diameter and lower elastic modulus in its radial direction), and to van der Waals forces. Shi et al. [8] developed a two-parameter model to describe the effect of agglomeration on the elastic properties of randomly oriented CNTs.
The present Eshelby-Mori-Tanaka property estimation model [8] is based on the equivalent fibre concept. The equivalent fibre properties are estimated by multiscale finite element method (FEM) analysis or through molecular dynamics (MD) simulations. Through these simulations, one determines the properties of the composite, and then the properties of the equivalent fibre might be calculated using the rule of mixtures (ROM) presented in Equation (1). A virtual-equivalent fibre consists of a straight CNT embedded in a polymeric resin and its interphase. For the SWCNT with a chiral index of C 2020, 6, 79 4 of 32 (10,10), the equivalent fibre would be a solid cylinder with a 1.424-nm diameter [10]. The properties of the equivalent fibre for the present studies are listed in Table 1.
where E LEF , E TEF , G EF , and ν LEF are the longitudinal elastic modulus, transverse elastic modulus, shear modulus, and Poisson's ratio, respectively, of the equivalent fiber. E LC , E TC , G C , and ν C are the longitudinal elastic modulus, transverse elastic modulus, shear modulus, and Poisson's ratio, respectively, of the composite, which were determined using multiscale FEM or MD simulations. E M , G M , and ν C are the elastic modulus, shear modulus, and Poisson's ratio, respectively, of the matrix. Finally, V EF and V M are the volume fraction of the equivalent fibre and the volume fraction of the matrix, respectively. This model considers that a part of the CNTs' reinforcement is dispersed throughout the matrix, while another part appears in a cluster form due to the agglomeration effect. This model is based on the Eshelby inclusion model, where it is assumed that the inclusions appear in spherical shapes and they are represented in Figure 1 [13]. The total volume of CNT reinforcement in the representative volume element (RVE) is denoted by V r and is divided in V inclusion r , which is the volume of CNTs in the agglomerated inclusions, and V m r , which is the volume of CNTs dispersed in the matrix. Note that from now on, the subscripts r and m will stand for the reinforcing phase and for the matrix, respectively. The employed homogenization scheme here described can be originally found in Section 4.1 of the reference article [8].
, , , and are the longitudinal elastic modulus, transverse elastic modulus, shear modulus, and Poisson's ratio, respectively, of the composite, which were determined using multiscale FEM or MD simulations. , , and are the elastic modulus, shear modulus, and Poisson's ratio, respectively, of the matrix. Finally, and are the volume fraction of the equivalent fibre and the volume fraction of the matrix, respectively.
This model considers that a part of the CNTs' reinforcement is dispersed throughout the matrix, while another part appears in a cluster form due to the agglomeration effect. This model is based on the Eshelby inclusion model, where it is assumed that the inclusions appear in spherical shapes and they are represented in Figure 1 [13]. The total volume of CNT reinforcement in the representative volume element (RVE) is denoted by and is divided in , which is the volume of CNTs in the agglomerated inclusions, and , which is the volume of CNTs dispersed in the matrix. Note that from now on, the subscripts and will stand for the reinforcing phase and for the matrix, respectively. The employed homogenization scheme here described can be originally found in Section 4.1 of the reference article [8].   The two parameters used to describe the agglomeration are µ and η, where µ denotes the volume fraction of the inclusions with respect to the volume V of the RVE, η is the volume ratio of CNTs that are dispersed in agglomerated inclusions with respect to the total volume of CNTs V r , and where V inclusion is the total volume of the inclusions within the RVE: when µ = 1, all the CNTs are uniformly dispersed in the matrix. With the decrease of µ, the agglomeration degree is more severe. If η = 1, all the CNTs are located within the agglomerated inclusions. If η = µ, then the volume fraction of CNTs inside the inclusions is the same as the volume fraction CNTs outside, which means the CNTs are uniformly distributed. When η > µ and the bigger the value of η, the spatial distribution of CNTs throughout the matrix is more heterogeneous.
The average volume fraction f r of CNTs in the composite is described by: Although assuming the nanotubes are transversely isotropic, considering that the nanotubes are randomly dispersed within the inclusions as well as within the matrix, one might consider that the material behaves as an isotropic material inside and outside the inclusions.
Using the Mori-Tanaka method for composite properties' estimation for randomly oriented CNTs and considering the agglomeration effect, one has: n r l r l r 0 0 0 l r k r + m r k r − m r 0 0 0 l r k r − m r k r + m r 0 0 0 0 0 0 where n r , l r , k r , m r , and p r are Hill's elastic moduli calculated by the inverse of the compliance matrix of the equivalent fibre. Note that E L , E T , E Z , G LT , G TZ , G TZ , and ν LT are its properties, which can be determined using ROM, but first the properties of the composite must be determined using a multiscale FEM or an MD simulation analysis [13]. The Hill's elastic moduli are used to determine the effective bulk moduli of the composite inside the inclusions K in and outside the inclusions K out , and the same goes for the shear moduli inside and outside the inclusions (G in and G out , respectively) [8]: where: C 2020, 6, 79 6 of 32 The bulk and shear moduli of the matrix are K m and G m , respectively. Finally, the effective bulk modulus K and the effective shear modulus G of the composite can be determined using the following expressions: where: where ν out is the Poisson's ratio outside the inclusions. Finally, one can determine the effective Young modulus E and Poisson's ratio ν of the composite using: Besides using the EMT agglomeration model to estimate the elastic properties of the FG-CNTRC, the composite mass density is calculated using the Voigt's rule [24,25]: where f m and ρ m are the matrix volume fraction and mass density, respectively. The matrix volume fraction is given by To calculate the mentioned properties, one must define how the CNTs are distributed along the matrix in the FG-CNTRC. For this purpose, three different CNTs volume fractions distributions were considered in this work, these volume fraction distributions are represented in Figure 2. The first is the uniform distribution (UD): An unsymmetrical functionally graded distribution (USFG) of the CNTs: C 2020, 6, 79 7 of 32 A symmetrical functionally graded distribution (SFG) of the CNTs: Knowing that z = − h 2 , h 2 , where h is the FG-CNTRC thickness and in which f * r is defined by [5]: where w r is the mass fraction of the carbon nanotubes.
C 2020, 6, x 7 of 33 An unsymmetrical functionally graded distribution (USFG) of the CNTs: A symmetrical functionally graded distribution (SFG) of the CNTs: Knowing that = − , , where ℎ is the FG-CNTRC thickness and in which * is defined by [5]: where is the mass fraction of the carbon nanotubes.

Free Vibration Analysis
As previously mentioned, this work intends to characterize the free vibration behaviour of quadrilateral nanocomposite plates where the effects of different CNT distributions and agglomeration are to be assessed. Finite element models were implemented using FSDT. The quadrilateral plate is represented in Figure 3.

Free Vibration Analysis
As previously mentioned, this work intends to characterize the free vibration behaviour of quadrilateral nanocomposite plates where the effects of different CNT distributions and agglomeration are to be assessed. Finite element models were implemented using FSDT. The quadrilateral plate is represented in Figure 3. An unsymmetrical functionally graded distribution (USFG) of the CNTs: A symmetrical functionally graded distribution (SFG) of the CNTs: Knowing that = − , , where ℎ is the FG-CNTRC thickness and in which * is defined by [5]: where is the mass fraction of the carbon nanotubes.

Free Vibration Analysis
As previously mentioned, this work intends to characterize the free vibration behaviour of quadrilateral nanocomposite plates where the effects of different CNT distributions and agglomeration are to be assessed. Finite element models were implemented using FSDT. The quadrilateral plate is represented in Figure 3.  The plates are considered to be simply supported in the four edges, with a = b for the square plate and b/a = 3 for the rectangular plate. The aspect ratio considered between the x-edge and thickness is a/h = 10.

Constitutive Equations
Considering the FSDT, the displacement field is given by [25][26][27]: where θ x and θ y are the first-order rotations and u 0 , v 0 , and w 0 are the mid-surface displacements in the (x, y, z) directions, respectively, and t denotes time. The superscript 0 is going to be used to refer to the mid-surface of the plate. The present theory provides five degrees of freedom. The strain field is given by [25][26][27]: where xx , yy , and γ xy are the in-plane normal and shear strain components; and γ yz and γ xz are transverse shear strains components.
In the context of a small strains approach, which is the one considered in the present work, the strain field relates with the displacements in the following manner: Yielding to: where σ xx , σ yy , and τ xy are the in-plane normal and shear stress components; τ yz and τ xz are transverse shear stress components; and K c is the shear correction factor associated with the transverse shear stress components, which for this study is assumed to have the typical value of 5/6. Considering the random orientation of CNTs, the corresponding nanocomposites can be considered as an isotropic material, thus the elastic stiffness coefficients associated to [Q] are given as:

Governing Equations of Motion
Hamilton's principle can be generically stated as [25]: where T is the kinetic energy, U is the elastic strain energy, and W is the work performed by applied forces. Considering the aim of the present work, the last term will be null, while the kinetic energy and the elastic strain energy are given as: By carrying out the necessary mathematical manipulations and considering free harmonic vibrations, one achieves for the whole discretized domain the following equilibrium Equation: is the global mass matrix, ω i is the natural frequencies, and {u i } denotes the vibration mode shapes' vectors associated to the corresponding natural frequencies ω i [25][26][27].
The global stiffness matrix can be expressed as: where [K mm ] e , [K bb ] e , [K mb ] e , and [K ss ] e are stiffness matrices associated to individualized effects, namely due to membrane, bending, membrane-bending coupling, and transverse shear, respectively. The index e refers to a generic element and N is the total number of elements considered. The stiffness matrices components are determined by the expressions below [26,27]: where [A], [B], and [D] are the extensional stiffness, the bending extensional coupling stiffness, and the bending stiffness matrices, where its components are given by: These matrices consider that the generalized displacements vector {u} for each node within each element is ordered as follows: where N i is the interpolating functions. In the initial phase of the present work, two plate elements were considered: the bi-linear (Q4) and bi-quadratic (Q9) quadrilateral plate elements with 4 and 9 nodes, respectively [26,27], thus leading to the corresponding interpolating functions: The Jacobian matrix [J e ] relates the local coordinates (ξ, η) of a generic element e with the global coordinates (x,y): The global mass matrix cab be written as [27]: The matrix [M] e is the mass matrix associated with a generic element e and can be divided into each of the degrees of freedom contributions, namely: C 2020, 6, 79 11 of 32 where I 0 and I 2 are the inertias associated with the translational and the rotational degrees of freedom, respectively, and they are determined by the following equation, where ρ c is the mass density of the nanocomposite [25,27]:

Results
This section is divided into two sub-sections: a first one devoted to the verification of the developed codes and a subsequent one that addresses the characterization of the free vibrations behaviour and natural modes' shapes of square and rectangular plates.

On the Use of Eshelby-Mori-Tanaka Agglomeration Model
For the present verification study, which concerns the mechanical properties' estimation for the FG-CNTRC, the EMT agglomeration model proposed by Shi et al. [8] was used. The Hill's moduli of the CNTs considered for this purpose are listed in Table 2. The material of the matrix considered in this case has the following elastic properties: E m = 0.85 GPa, ν m = 0.3. This study was performed considering the experimental results obtained by Odegard et al. [29], where the material properties of the present materials were tested for different values of volume fraction of the reinforcement.
In Figure 4, we observe that the experimental results obtained by Odegard et al. [29] are close to the agglomerated model, where µ = 0.4 as predicted by Barai et al. [13,30] when considering a complete agglomerated model where η = 1, meaning that all CNTs are located in inclusions, showing that the employed agglomeration model based on the EMT method is in good agreement with the previous literature (note that the considered points of the experiment led by Odegard et al. were attained through graphic reading and might have slight differences).
In Figure 4, we observe that the experimental results obtained by Odegard et al. [29] are close to the agglomerated model, where µ=0.4 as predicted by Barai et al. [13,30] when considering a complete agglomerated model where η=1, meaning that all CNTs are located in inclusions, showing that the employed agglomeration model based on the EMT method is in good agreement with the previous literature (note that the considered points of the experiment led by Odegard et al. were attained through graphic reading and might have slight differences). Taking the fully dispersed case as a reference, µ=1, where the Young's modulus has the higher increase in function of the volume fraction, it is possible to observe that when µ decreases, the increase in the CNT volume fraction does not correspond to the expected increase of mechanical performance due to the severity of the agglomeration effect.
However, partial CNT agglomeration is a more common situation, with this meaning that both parameters η and µ are necessary to describe the agglomeration state [8]. The following graphics represented in Figure 5 show how the Young's modulus and the Poisson's ratio behave with the variation of these agglomeration parameters. For this illustration purpose, one considered the CNT volume fractions of 0.1%, 0.5%, and 1%, because they coincide with the points of the experiment led by Odegard et al.
One can observe that the highest values of the Young' modulus appear when µ = η. It is also possible to observe that the variation of the parameter µ has a higher impact on the elastic properties, when compared to η. Taking the fully dispersed case as a reference, µ=1, where the Young's modulus has the higher increase in function of the volume fraction, it is possible to observe that when µ decreases, the increase in the CNT volume fraction does not correspond to the expected increase of mechanical performance due to the severity of the agglomeration effect.
However, partial CNT agglomeration is a more common situation, with this meaning that both parameters η and µ are necessary to describe the agglomeration state [8]. The following graphics represented in Figure 5 show how the Young's modulus and the Poisson's ratio behave with the variation of these agglomeration parameters. For this illustration purpose, one considered the CNT volume fractions of 0.1%, 0.5%, and 1%, because they coincide with the points of the experiment led by Odegard et al.
One can observe that the highest values of the Young' modulus appear when µ = η. It is also possible to observe that the variation of the parameter µ has a higher impact on the elastic properties, when compared to η.

On the Use of FSDT for Free Vibration Analysis
The following verification studies consider a unit length square plate as schematically illustrated in Figure 3 with opposite hinged and simply supported edges in the directions of the plate plane (x,y). Two edge-to-thickness ratios were considered, namely a/h = 10 and a/h = 100.
The material considered for the studies is considered as an isotropic material, and its Poisson's ratio assumes the value of ν = 0.3, and the Young's modulus and the specific mass E and ρ are considered in their SI units [31].
The results are compared in Table 3 with an analytical solution obtained through the Rayleigh-Ritz method for vibration analysis of Mindlin plates [31]. The frequencies are presented in a non-dimensional form using the multiplier: C 2020, 6, 79 13 of 32 C 2020, 6, x 13 of 33

On the Use of FSDT for Free Vibration Analysis
The following verification studies consider a unit length square plate as schematically illustrated in Figure 3 with opposite hinged and simply supported edges in the directions of the plate plane (x,y). Two edge-to-thickness ratios were considered, namely a/h = 10 and a/h = 100.
The material considered for the studies is considered as an isotropic material, and its Poisson's ratio assumes the value of = 0.3, and the Young's modulus and the specific mass and are considered in their SI units [31].  A convergence study is also presented in the same table.
From Table 3, it is possible to conclude that for both elements, the results are very approximate, if not coincident to the analytic solution, which was expected since the analytic solution presented is for a Mindlin plate. However, one can observe that the element Q9 gives better/faster results even with a smaller number of elements, although the number of degrees of freedom is considerably higher.
Even when comparing the same number of degrees of freedom for the meshes Q4 10 × 10 and Q9 5 × 5 or Q4 20 × 20 and Q9 10 × 10, for example, one can observe a faster convergence for the element Q9, with this happening due the quadratic degree of the interpolating functions that provides a better approximation and thus a better performance for the Q9 plate model.

On the Verification of the FSDT Model for Free Vibration Analysis of an FG-CNTRC Beam Estimated with the EMT Agglomeration Model
The beam was discretized using the Q4 and Q9 plate element models, considering that in the y direction, there is only one element. The natural frequencies convergence study was developed by increasing the number of elements in the x direction.
To allow for the intended verification purpose, the selected beam was the one studied by Yas and Heshmati [5], and it is illustrated in Figure 6. The frequencies are presented in their dimensionless form by considering the multiplier: where A is the area and I is the second moment of area of the cross-section of the beam. The reinforcement material considered in the next verification studies was the SWCNT with a chiral index of (10,10) presented in Table 1, for the matrix phase, the material properties are = 10 , = 0.3, and = 1150 / . The length-to-thickness aspect ratio considered was L/h = 20 and * = 0.075. The length-to-thickness relation, L/b, is equal to 50. The other authors considered randomly oriented straight CNTs in the matrix, without considering the agglomeration effect. For this purpose, the material properties of the composite were calculated using the EMT agglomeration model using its parameters µ = η, which ensures a uniformly dispersed state of the carbon nanotubes.
To ensure the results of the analysis of the natural frequencies, a convergence study was made considering the distribution UD and using a clamped-clamped (C-C) boundary condition (B.C.). Table 4 presents the results for the bi-linear element Q4 and for the bi-quadric element Q9. The The reinforcement material considered in the next verification studies was the SWCNT with a chiral index of (10,10) presented in Table 1, for the matrix phase, the material properties are E = 10 GPa, ν = 0.3, and ρ = 1150 Kg/m 3 . The length-to-thickness aspect ratio considered was L/h = 20 and f * r = 0.075. The length-to-thickness relation, L/b, is equal to 50.
The other authors considered randomly oriented straight CNTs in the matrix, without considering the agglomeration effect. For this purpose, the material properties of the composite were calculated using the EMT agglomeration model using its parameters µ = η, which ensures a uniformly dispersed state of the carbon nanotubes.
To ensure the results of the analysis of the natural frequencies, a convergence study was made considering the distribution UD and using a clamped-clamped (C-C) boundary condition (B.C.). Table 4 presents the results for the bi-linear element Q4 and for the bi-quadric element Q9. The deviations (dev) were calculated as: In this convergence study, it is observed that using the bi-linear element Q4, the convergence is slower than when using the bi-quadratic element Q9, for the same number of elements. However, the element Q9 has a higher computational cost due to its higher number of nodes per element (nine against four).
The results of the convergence studies seem to be in reasonable agreement with the results of the considered reference (see [5]) mainly for the first natural frequency. The deviations increase for higher frequencies. These results are expected as one is using 2-D elements, which provide a greater stiffness to the computational model. Additionally, it will be expected that for higher order frequencies, modes that are not expected to appear using Timoshenko beam theory will appear. It is the case of the torsional mode shown in Figure 7.
Additionally, to verify the three CNT volume fraction distributions (UD, USFG, and SFG) and other boundary conditions for the present problem, one performed another study considering the present models to be compared with the Euler Bernoulli beam element and the Timoshenko beam element. In Tables 5 and 6, where the first three natural frequencies are presented, where one reads H-H, C-F, and C-H, it must be understood as hinged-hinged, clamped-free, and clamped-hinged, respectively.
Despite the convergence study considered until 100 elements, which was the mesh considered by the reference author [5], one has considered a mesh of 50 elements for the next verification with the Q9 element as it provides lower deviations, and also because from 50 to 100 elements, the deviations' improvement was minimal.
The results presented seem to be in reasonable agreement with the reference, considering that different assumptions are involved when considering a plate and a beam model. The deviations were calculated considering the Timoshenko beam theory as a reference. As expected, they are lower when compared with the Euler-Bernoulli beam results, since the Euler-Bernoulli beam theory does not consider transverse shear components. The present results are higher when compared with the Timoshenko beam theory results, since the 1-D model used to implement it provides a less stiff structure when compared to the 2-D approach considered in the present study. Additionally, to verify the three CNT volume fraction distributions (UD, USFG, and SFG) and other boundary conditions for the present problem, one performed another study considering the present models to be compared with the Euler Bernoulli beam element and the Timoshenko beam element. In Tables 5 and 6, where the first three natural frequencies are presented, where one reads H-H, C-F, and C-H, it must be understood as hinged-hinged, clamped-free, and clamped-hinged, respectively.
Despite the convergence study considered until 100 elements, which was the mesh considered by the reference author [5], one has considered a mesh of 50 elements for the next verification with the Q9 element as it provides lower deviations, and also because from 50 to 100 elements, the deviations' improvement was minimal.   In this last case, an important trend was identified by Yas and Heshmati [5], in which the frequencies obtained when using the SFG distribution were higher than when using the USFG and the UD, and it was stated that this phenomenon was caused by a better use of the CNTs' distribution in higher bending stress regions, where they are more concentrated, thus concluding that its bending stiffness was larger than for USFG and UD distributions. In the present study, we can observe the same phenomenon. In this case study, it was also found that there was one more natural frequency between the second and the third natural frequencies considered by Yas and Heshmati [5], which is associated with the natural modes of vibration of a plate. The first four modes of vibration for a UD and (C-C) using the Q9 element type are illustrated in the Figure 7. As expected, the use of a plate model to characterize beams' natural frequencies is able to provide more information about its modes of vibration.

Effect of the Agglomeration on FG-CNTRCs Square Plates' Free Vibrations Behaviour
In the following subsubsections, the free vibrations behaviour of a simply supported square plate as illustrated in the Figure 3 was evaluated using the element Q4 with 20 × 20 elements and Q9 with 15 × 15 elements, where finer meshes were not considered due to its high computational cost and because satisfactory results were obtained with coarser meshes as concluded in Section 3.1.2. The plate's aspect ratio a/h is equal to 10.
The material of the matrix considered in this study has the following elastic properties: E m = 2.1 GPa, ν m = 0.34, and ρ m = 1150 kg/m 3 , and the material properties of the reinforcement are listed in Table 2. A value of f * r = 0.075 was considered, as it was found that for a 7.5% concentration, a large amount of CNTs are concentrated in inclusions [23]. The reinforcement distributions considered were the UD, USFD, and the SFD with different levels of agglomeration tested. The dimensionless frequencies to be presented were obtained using the following expression:

Free Vibration Analysis without the Agglomeration Effect
The present study does not consider an agglomerated state, µ = η, and its results serve as a reference for the next ones that consider agglomerated states. For the UD distribution, the results are presented in Table 7. It is possible to observe that when compared to the other two distributions, the SFG provides the best vibrational characteristics, since its natural frequencies assume higher values. This behaviour is attained because the CNTs are in a higher concentration distributed to higher stress regions.
It is also possible to observe from this table when comparing the results from the Q4 element with the Q9 element, we can see that the frequencies tend to be higher for the Q4 element, which correspond to the expected behaviour.
The next figures present the plots corresponding to the results obtained using only the Q9 plate element, since the results with the Q4 element were similar and the slight differences in the natural frequencies do not interfere with the modes' shape. Figure 8 shows the first vibrational modes for the UD distribution of CNTs without the agglomeration effect, noting that the third mode was omitted since it is symmetrical with the second. It is possible to observe that when compared to the other two distributions, the SFG provides the best vibrational characteristics, since its natural frequencies assume higher values. This behaviour is attained because the CNTs are in a higher concentration distributed to higher stress regions.
It is also possible to observe from this table when comparing the results from the Q4 element with the Q9 element, we can see that the frequencies tend to be higher for the Q4 element, which correspond to the expected behaviour.
The next figures present the plots corresponding to the results obtained using only the Q9 plate element, since the results with the Q4 element were similar and the slight differences in the natural frequencies do not interfere with the modes' shape. Figure 8 shows the first vibrational modes for the UD distribution of CNTs without the agglomeration effect, noting that the third mode was omitted since it is symmetrical with the second.  Table 7, one can observe that this distribution has multiple modes, which is due to the symmetries involved in the analysed plate.
In Figure 9, one can see the first five modes for the USFG distribution. Although some frequencies are close, there is a slight difference that may be due to the CNT asymmetry distribution.
In Figure 10, the first modes for the SFG distribution are illustrated, where the third mode was not displayed since it is symmetrical with the second mode.
Considering the better performance of the biquadratic (Q9) plate model when compared to the bilinear one (Q4), we will only consider Q9 in the next studies.  Table 7, one can observe that this distribution has multiple modes, which is due to the symmetries involved in the analysed plate.
In Figure 9, one can see the first five modes for the USFG distribution. Although some frequencies are close, there is a slight difference that may be due to the CNT asymmetry distribution. In Figure 10, the first modes for the SFG distribution are illustrated, where the third mode was not displayed since it is symmetrical with the second mode.  In Figure 10, the first modes for the SFG distribution are illustrated, where the third mode was not displayed since it is symmetrical with the second mode.

Free Vibration Analysis for Complete Agglomerated States (η = 1)
This study comprehends a complete state of agglomeration, when η = 1 for three different values of µ = {0.25, 0.5, 0.75}. In these agglomerated states, it is considered that all the CNTs are aggregated in the spherically shaped inclusions [8].
The modes' shapes for agglomerated states will just be displayed if an observable change on them was identified for any of these three CNT volume fraction distributions.
The results presented in Table 8 are associated to the UD distribution. From the three states of agglomeration considered, it is possible to observe that the worse vibrational behaviour is precisely when µ = 0.25, which corresponds to the most heterogeneous distribution of CNTs and where the level of agglomeration is more severe.
As it was foreseen in Section 3.1.1, the bigger the difference between the values of the agglomeration parameters, the more its elastic properties would be affected by the agglomeration of the CNTs.
The differences between the natural frequencies are very significant, when comparing these two states of complete agglomeration. When comparing these same frequencies with the ones obtained for a non-agglomerated state in Section 3.2.1, the difference is even higher. Although, the agglomeration effect seems to not interfere with the modes' shape for this distribution of CNTs along the thickness.
The natural frequencies obtained for the three cases of complete agglomeration considering the USFG distribution are listed in Table 9. From this table, one can conclude that for all cases of complete agglomeration studied, the USFG is the CNT distribution that has the worse dynamic behaviour, when comparing it with the same states of agglomeration for the other CNT distributions. Similarly to the previous case, the complete agglomerated situations for the USFG CNT volume fraction distribution seem not to interfere significantly with the natural modes' shape for this distribution.
For the SFG distribution (Table 10), one can observe that, as predicted, the lower the value of the parameter µ, the lower the natural frequencies, just like for the previous distributions. The modes' shapes for the SFG distribution are presented in Figure 11 for the severest case of complete agglomeration and in Figure 12 for the least severe case of complete agglomeration studied. The intermediate case is not displayed, since it has the same behaviour as the other cases. For these three cases of complete agglomeration, it is possible to observe a significant change of the fifth mode's shape, when compared with the non-agglomerated case. For the other three modes' shapes, no significant changes were observed.
agglomeration. λ η = 1 μ = 0.25 η = 1 μ = 0.5 η = 1 μ = 0.75 The modes' shapes for the SFG distribution are presented in Figure 11 for the severest case of complete agglomeration and in Figure 12 for the least severe case of complete agglomeration studied. The intermediate case is not displayed, since it has the same behaviour as the other cases. For these three cases of complete agglomeration, it is possible to observe a significant change of the fifth mode's shape, when compared with the non-agglomerated case. For the other three modes' shapes, no significant changes were observed.
Globally, one can say that for a completely agglomerated case, the more heterogeneous the distribution, the lower the natural frequencies will be for the three CNT distributions studied.
It is also possible to conclude that besides the level of agglomeration, the SFG distribution shows a better vibrational behaviour due to its distribution of CNTs in higher bending stress areas; however, the lesser the value of µ, the lower the differences of the natural frequencies become between distributions.   Globally, one can say that for a completely agglomerated case, the more heterogeneous the distribution, the lower the natural frequencies will be for the three CNT distributions studied.
It is also possible to conclude that besides the level of agglomeration, the SFG distribution shows a better vibrational behaviour due to its distribution of CNTs in higher bending stress areas; however, the lesser the value of µ, the lower the differences of the natural frequencies become between distributions.

Free Vibration Analysis for Partially Agglomerated States
In a more general situation, one has situations that are distinct from the completely nonagglomerated state or from the completely agglomerated state, hence being very important to achieve a plausible description of the degree of agglomeration through the parameters µ and η [8].
To investigate the free vibrations' behaviour of simply supported FG-CNTRC square plates, two different partially agglomerated situations were evaluated, one with η<µ and the other with η > µ, when µ = 0.5 for both situations. Table 11 presents the first five natural frequencies for the UD distribution.

Free Vibration Analysis for Partially Agglomerated States
In a more general situation, one has situations that are distinct from the completely non-agglomerated state or from the completely agglomerated state, hence being very important to achieve a plausible description of the degree of agglomeration through the parameters µ and η [8].
To investigate the free vibrations' behaviour of simply supported FG-CNTRC square plates, two different partially agglomerated situations were evaluated, one with η < µ and the other with η > µ, when µ = 0.5 for both situations. Table 11 presents the first five natural frequencies for the UD distribution. From this table, when comparing the two states of partial agglomeration, the best free vibration behaviour is achieved with µ = 0.5 and η = 0.25, which corresponds to a state where there is less volume fraction of CNTs in agglomerated inclusions, when compared with µ = 0.5 and η = 0.75. Although significant differences were found for the natural frequencies when comparing both states of agglomeration, and when comparing these agglomerated states with the results obtained in Section 3.2.1 for a non-agglomerated state, the natural modes' shapes do not suffer much influence from the agglomeration effect, for these two partially agglomerated states. In Table 12, one presents the first five natural frequencies for the USFG distribution. In this case the highest natural frequencies appear for η = 0.25, and once again the lesser the volume of CNTs inside the agglomerated inclusions, the better the dynamic free vibrations' behaviour obtained in the CNTRC. It was observed that the natural modes' shapes remained unchanged despite the agglomeration effect when compared to the results obtained for the non-agglomerated state in Section 3.2.1, for the USFG distribution. Table 13 presents the first five natural frequencies for the SFG distribution. It was also observed in this case that the highest natural frequencies appear for η = 0.25, which translates into a better vibrational behaviour. For this last distribution, it was observed that the modes' shape remained unchanged despite the agglomeration effect, for these two cases of partial CNTs agglomeration. The fifth mode change, when in complete agglomeration, may indicate that the more severe the agglomeration, the more prone the interference in the modes' shapes.

Effect of the Agglomeration on FG-CNTRCs Rectangular Plates' Free Vibration Behaviour
In the present case study, the free vibration behaviour of a simply supported rectangular plate was explored under the influence of the agglomeration of CNTs for three different distributions along the thickness direction. The plate is illustrated in Figure 3, with the aspect ratios of a/h = 10 and a/b = 3. The dimensionless results were obtained according to Equation (54), as done in Section 3.2.

Free Vibration Analysis without the Agglomeration Effect
The present study does not consider an agglomerated state, µ = η, in order to consider the achieved results as a reference for the next one, which will consider agglomerated states. Since here it is considered that there is no influence of the agglomeration effect, this study illustrates the best free vibrations' performance of these rectangular FG-CNTRC plates.
The natural frequencies are presented in Table 14. As expected, one observes that there are no multiple modes for any of the distributions considered, with this happening due to the symmetry loss. It is also possible to say that for the rectangular plate, the tendency of the SFG distribution having the higher natural frequencies is maintained followed by the UD distribution and at last the USFG distribution with the poorest dynamic behaviour. This results from the distribution of CNTs along the thickness direction, as the higher the concentration of CNTs in high bending stress areas, the higher the natural frequencies. Figure 13 shows the natural modes' shape without the influence of the CNTs' agglomeration for the UD distribution. It is also possible to say that for the rectangular plate, the tendency of the SFG distribution having the higher natural frequencies is maintained followed by the UD distribution and at last the USFG distribution with the poorest dynamic behaviour. This results from the distribution of CNTs along the thickness direction, as the higher the concentration of CNTs in high bending stress areas, the higher the natural frequencies. Figure 13 shows the natural modes' shape without the influence of the CNTs' agglomeration for the UD distribution. Figure 14 shows the natural modes' shape without the influence of the CNTs' agglomeration for the USFG distribution. Figure 15 shows the natural modes' shape without the influence of the CNTs' agglomeration for the SFG distribution.   Figure 14 shows the natural modes' shape without the influence of the CNTs' agglomeration for the USFG distribution.     Figure 15 shows the natural modes' shape without the influence of the CNTs' agglomeration for the SFG distribution.   Similarly to the squared plate, for the rectangular plate, the same three states of complete agglomeration, when η = 1 for three different values of µ = {0.25, 0.5, 0.75}, were also considered.
The results listed in Table 15 correspond to the UD distribution. The free vibrations' behaviour of the rectangular plate tends to worsen with the increase of the severity of the complete agglomerated state, showing lower natural frequencies with the decrease of the parameter µ.  The results listed in Table 15 correspond to the UD distribution. The free vibrations' behaviour of the rectangular plate tends to worsen with the increase of the severity of the complete agglomerated state, showing lower natural frequencies with the decrease of the parameter µ. For this case, one can observe that the agglomeration effect exerts influence in the values of the natural frequencies, although the natural modes' shape seems to be unchanged when comparing to the shapes presented in Section 3.3.1 for the non-agglomerated state for this distribution.
For the USFG distribution, the results for the natural frequencies obtained for the three cases of complete agglomeration can be observed in Table 16. From this table, one can conclude that for all cases of complete agglomeration studied, the USFG is the distribution that has the worse dynamic behaviour for the rectangular plate, when compared with the same states of complete agglomeration for the other CNT distributions. The modes' shapes for the USFG distribution are presented in Figure 16 for the most severe case of complete agglomeration and in Figure 17 for the intermediate case of complete agglomeration studied. For this distribution, it is possible to see that the agglomeration effect, besides decreasing the natural frequencies, also changes the shapes of the natural modes in some cases, namely the fourth and the fifth mode of the USFG using η = 1 and µ = 0.25, and in the fifth mode when η = 1 and µ = 0.5, when compared with the non-agglomerated state presented in Section 3.3.1. These differences are apparently intensified with the agglomeration state, since for the least severe case, no differences were observed; for the intermediate situation, a difference in the fifth mode was observed; and at last, for the severest situation of complete agglomeration, differences in the fourth and fifth mode were observed.
when compared with the non-agglomerated state presented in Section 3.3.1. These differences are apparently intensified with the agglomeration state, since for the least severe case, no differences were observed; for the intermediate situation, a difference in the fifth mode was observed; and at last, for the severest situation of complete agglomeration, differences in the fourth and fifth mode were observed.  when compared with the non-agglomerated state presented in Section 3.3.1. These differences are apparently intensified with the agglomeration state, since for the least severe case, no differences were observed; for the intermediate situation, a difference in the fifth mode was observed; and at last, for the severest situation of complete agglomeration, differences in the fourth and fifth mode were observed.  For the SFG distribution (Table 17), one can observe that, as predicted, the lower the value of the parameter µ, the lower are the natural frequencies, just like for the previous distributions. Once again, the best dynamic behaviour seems to prevail for this distribution regardless of the agglomeration effect when compared with the same agglomeration states for the USFG and UD distributions.   For the SFG distribution (Table 17), one can observe that, as predicted, the lower the value of the parameter µ, the lower are the natural frequencies, just like for the previous distributions. Once again, the best dynamic behaviour seems to prevail for this distribution regardless of the agglomeration effect when compared with the same agglomeration states for the USFG and UD distributions. For this distribution, one can observe that the agglomeration effect has a high impact on the natural frequencies, just like for the other distributions. However, the natural modes' shape did not present significant differences when compared to the non-agglomerated situation.

Free Vibration Analysis for Partially Agglomerated States
Similarly to the study developed for the square plate, the free vibrations' behaviour of simply supported FG-CNTRC rectangular plates, under different partially agglomerated situations were evaluated, one with η < µ and the other with η > µ, when µ = 0.5 for both situations. As already mentioned, partial agglomeration is more common than complete agglomeration in CNTRC.
In Table 18, we present the first five natural frequencies for the UD distribution. For this CNTs' distribution, when comparing the two states of partial agglomeration, the best behaviour found is achieved with µ = 0.5 and η = 0.25, which corresponds to a state where there is less volume fraction of CNTs located in agglomerated inclusions, when compared with µ = 0.5 and η = 0.75.
Significant differences are presented for the natural frequencies when comparing both states of agglomeration, and when comparing these agglomerated states with the results obtained in Section 3.3.1 for a non-agglomerated state. However, in terms of the natural modes' shapes, no significant differences were found for this partially agglomerated situation when compared with the non-agglomerated state. Table 19 presents the first five natural frequencies for the USFG distribution. For the rectangular plate with the USFG distribution, the tendency for the highest natural frequencies to appear for η = 0.25 is maintained, and once again, the lesser the volume of CNTs inside the agglomerated inclusions, the better the dynamic behaviour of the plate. For these two situations of partial agglomeration, when comparing the modes' shapes results with the ones for a non-agglomerated state (Section 3.3.1), no significant differences were observed for this CNTs unsymmetrical distribution along the thickness direction, different to the previous results for the complete agglomeration (Section 3.3.2), which is another indicator that a more severe state of agglomeration can be more prone to influencing the natural modes.
Finally, in Table 20, we present the first five natural frequencies for the SFG distribution. For this last vibration analysis, it was observed that this distribution provides the highest natural frequencies, which translates into a better vibrational behaviour, when compared to the other distributions due to its enhanced bending stiffness properties. When comparing both states of agglomeration, one can say that the best dynamic behaviour is reached for η = 0.25, when there is a lesser volume of agglomerated CNTs' inclusions.
For this last distribution, one can say that the shape of the natural modes remained unchanged despite the agglomeration effect when compared to the results obtained for the non-agglomerated state in Section 3.3.1; however, as already mentioned, there are significant differences in the values of the natural frequencies.

Discussion
From the case studies considered, the element Q9 tends to provide lower natural frequencies when compared to the element Q4, due to its greater ability to describe the plate's shape kinematics, yielding also to more expressive differences for higher order frequencies.
Considering the results of the free vibrations' behaviour of the FG-CNTRC quadrilateral plates, one can say:

•
Since CNTs tend to agglomerate for relatively low CNT volume fractions, this leads to an important conclusion, namely the one that not taking into account the effect of the CNTs' agglomeration might mislead into an overestimation of the elastic properties of the CNTRC, thus to a less accurate prediction of structure behaviour. • For all CNTs' distribution in non-agglomerated situations, the natural frequencies of these structures were always higher when compared with the results for the same CNT distribution if under an agglomerated condition whether it would be a complete or a partial one. • For all cases considered, i.e., no agglomeration, three different states of complete agglomeration and the two different states of partial agglomeration, the CNTs' SFG distribution along the thickness direction provided higher natural frequencies, when comparing to the other two distributions considered for the same state of agglomeration. This happens due to the higher concentration of CNTs in high bending stress regions.

•
The results of the complete agglomeration studies also demonstrate that the decrease of the agglomeration parameter µ deteriorates the free vibrations' behaviour of these structures, leading to lower natural frequencies for all the distributions considered.

•
When considering partial agglomeration of the CNTs in the FG-CNTRC plates, it can be concluded that the higher the agglomeration parameter η, the lower the natural frequencies for the CNTs' distributions considered. With respect to the natural modes' shape associated to the different natural frequencies: • On the FG-CNTRC square plates, no significant differences in the modes' shapes due to the influence of the agglomeration effect were found for the majority of the cases. However, for the SFG CNTs distribution, the fifth mode was affected for complete agglomeration.

•
On the FG-CNTRC rectangular plates, for the SFG and UD CNTs distributions, the natural modes' shapes seem not to be affected with the agglomeration of CNTs; however, for the USFG distribution, differences in the fourth and fifth modes for the complete agglomerated CNTs were found. Table 21 summarizes the observed differences in the modes' shapes considering the parametric studies performed. For the majority of the agglomeration cases studied, the agglomeration effect did not show an influence on the modes' shape. Table 21. Summary of the agglomeration effect in the natural modes shapes of the quadrilateral plates.

Square Plate
CNTs distribution η = 1 µ = 0. For the cases where differences were observed, they appear for higher modes (fourth and fifth) and for situations of complete agglomeration of the CNTs. As for the CNTs' volume fraction distributions, for the square plate, the observed differences occur in the SFG distribution, while for the rectangular plate, the observed differences appeared in the USFG distribution. This allows the conclusion that depending on the agglomeration state, its influence on the natural modes' shapes also depend on the geometry of the plate and on the CNTs' distribution along the thickness.