Investigation of Mechanical Behaviors of Functionally Graded CNT-Reinforced Composite Plates

In this paper, the mechanical behavior of a functionally graded carbon nanotube-reinforced composite (FG-CNTRC) plate is numerically investigated. According to the concept of a hierarchical model, the displacement is decomposed into the in-field functions and the assumed thickness-wise monomial. The former is defined on the plate midsurface and is approximated by the 2-D meshfree natural element method (NEM). The FG-CNTRC plate is modeled as a homogenized orthotropic body, and its effective elastic properties are determined by referring to MD simulation and the linear rule of mixtures. Regarding the thickness-wise distribution of CNTs, one uniform and three functionally gradient distributions are taken. Through comparative numerical experiments, the reliability of the numerical method is justified with the maximum relative difference of 6.12%. The effects of the volume fraction and vertical distribution of CNTs, the plate width-thickness and aspect ratios, and the boundary conditions on the bending, free vibration, and buckling behaviors of FG-CNTRC plates are examined. It is found that the mechanical behavior of FG-CNTRC plates is significantly dependent of these major parameters.


Introduction
Carbon nanotubes (CNTs) have been spotlighted as an innovative material for the 21st century due to their outstanding thermo-mechanical properties [1]. They are fabricated by rolling graphene sheets into the form of a cylinder and are divided into single-and multi-walled carbon nanotubes according to the number of graphene sheets. The excellent properties not only promise CNTs as next-generation multi-functional reinforcements for polymer composites [2] but also remarkably extend the conventional polymer composites to a variety of engineering fields. As a representative application, carbon nanotube-reinforced composites (CNTRCs) were recently introduced, and extensive studies have focused on their fabrication methodology and thermo-mechanical behavior [3]. A simple technique for fabricating the aligned arrays of CNTs was presented by Ajayan et al. [4] by utilizing cutting thin slices. The elastic properties of CNTRC were predicted by Hu et al. [5] through structural deformation analyses for various loading conditions and were also evaluated by Han and Elliott [6] through molecular dynamics (MD) simulation.
However, the success in applying CNTRCs to engineering applications relies on the in-depth investigation of their thermo-mechanical behaviors because CNTRCs can have various forms of structural elements under various loading and boundary conditions. Wuite and Adali [7] presented the mechanical responses of a CNTRC beam that were predicted by multiscale simulation using laminated beam theory. Formica et al. [8] presented an equivalent continuum theory to investigate the vibration behavior of CNTRCs, and Arani et al. [9] investigated the buckling behavior of laminated CNTRC plates by combining analytical and FE methods. Wang and Shen [10] applied the higher-order shear Figure 1a depicts a rectangular polymer plate in which SWCNTs are inserted uniformly through the thickness, where CNTs are parallel to the x− direction and the dimensions of the plate are length a, width b, and thickness d. Figure 2b represents three different functionally graded (FG) CNT distributions: FG-V, FG-O, and FG-X. Since CNTRC plates can be viewed as a two-phase composite, their equivalent elastic properties can be estimated using the rules of mixtures or the Mori-Tanaka method [12]. In this study, a modified linear rule of mixtures in which the CNT efficiency parameters η j (j = 1, 2, 3) are introduced is employed. Based on this method, the equivalent elastic and shear moduli of the FG-CNTRC plate are calculated as follows [13]: are introduced is employed. Based on this method, the equivalent elastic and shear moduli of the FG-CNTRC plate are calculated as follows [13]: indicate the volume fractions of CNTs and the polymer matrix, respectively. In Equations (1)-(3), the orthotropic elastic properties of CNT are labeled as CNT while those of the isotropic polymer matrix are denoted by m . The scaledependence of the equivalent elastic properties of CNTRCs is reflected by the CNT efficiency parameters j  which were calculated by equating the equivalent elastic properties predicted by MD simulation with those obtained by the rule of mixtures [28]. The volume fraction CNT V of carbon nanotubes in the UD-and FG-CNTRC plates through the thickness is given by are introduced is employed. Based on this method, the equivalent elastic and shear moduli of the FG-CNTRC plate are calculated as follows [13]: indicate the volume fractions of CNTs and the polymer matrix, respectively. In Equations (1)-(3), the orthotropic elastic properties of CNT are labeled as CNT while those of the isotropic polymer matrix are denoted by m . The scaledependence of the equivalent elastic properties of CNTRCs is reflected by the CNT efficiency parameters j  which were calculated by equating the equivalent elastic properties predicted by MD simulation with those obtained by the rule of mixtures [28]. The volume fraction CNT V of carbon nanotubes in the UD-and FG-CNTRC plates through the thickness is given by V CNT and V m = 1 − V CNT indicate the volume fractions of CNTs and the polymer matrix, respectively. In Equations (1)-(3), the orthotropic elastic properties of CNT are labeled as CNT while those of the isotropic polymer matrix are denoted by m. The scaledependence of the equivalent elastic properties of CNTRCs is reflected by the CNT efficiency parameters η j which were calculated by equating the equivalent elastic properties predicted by MD simulation with those obtained by the rule of mixtures [28].
The volume fraction V CNT of carbon nanotubes in the UD-and FG-CNTRC plates through the thickness is given by Here, w CNT denotes the mass fraction of CNTs within the CNTRC plate, and ρ m and ρ CNT are the densities of the matrix and CNT, respectively.
In a similar manner, the equivalent Poisson's ratio ν 12 and the equivalent density ρ of the CNTRC plate are calculated as follows:

Analysis of Bending, Free-Vibration, and Buckling
The displacement u = u x , u y , u z T in the bending, free-vibration, and buckling of the FG-CNTRC plate is expressed by a (1,1,0)* hierarchical model [26], which is equivalent to the first-order SDT.
Here, Θ m α (x, y) and (2z/d) m are 2-D in-plane vector functions and 1-D assumed thickness monomials, respectively. Figure 2a demonstrates a uniform NEM grid consisting of N grid points (nodes) and M Delaunay triangles. The NEM grid is generated on the plate midsurface ω (see Figure 1a). For each grid point, the Laplace interpolation (L/I) function ϕ J (x) illustrated in Figure 2b is assigned with the help of Delaunay triangulation and the Voronoi diagram [27,29,30]. With a finite number of L/I functions, both actual and virtual displacements u h and v h are approximated as Here, q x , q y , q z = (1, 1, 0) and U m α,J are the nodal displacements corresponding to the thickness monomial (2z/d) m of u α at node J.
Under the assumption of q x = q y = q z = q for the convenience of concise expression, the virtual strain vector ε v h and the actual stress vector σ u h are denoted as with the (6 × 6) orthotropic elastic constant matrix D [31]. The partial differential operator B J is defined by with ϕ J,α = ∂ϕ j /∂α. Then, the principle of virtual work leads to the coupled simultaneous equations defined by KU = F (16) to compute the nodal displacements {U} m β,J = U m x.J , U m y,J , U m z,J T . In Equation (16), the stiffness matrix K and the load vector F are defined by The subscript RI in Equation (17) denotes the selectively reduced integration (SRI) [32,33] using one Gauss point to avoid shear locking for the bending-dominated thin structure. The numerical integration in Equation (17) is carried out by combining the thickness-wise analytical integration and the in-plane Gauss numerical integration. For the SRI integration, the elastic constant matrix D should be divided as follows: where G 23 = G 23 /κ and G 31 = G 31 /κ with κ = 6/5 being the modified shear correction factor. For the free-vibration analysis of the FG-CNTRC plate, the harmonic response u(x; t) = u(x) · e jωt is considered, together with the mass matrix [M] m αβ,I J defined by with the (3 × 3) identity matrix I. Then, the weak formulation of the eigenvalue problem provides the modal equations given by for computing the natural frequencies ω J N J=1 and natural mode s U J N J=1 . Next, the following virtual work principle is considered to derive the discretized equations of the linear buckling problem: Here, the virtual strain energy δU and the virtual work δW ex are given as follows: Substituting the previous Equations (11)- (14) into Equations (23) and (24) results in the eigenvalue equations given by to calculate the buckling loads λ J N J=1 and the associated buckling modes U J N J=1 . Here, the geometric stiffness matrix K G is defined by Here, (γ 1 , γ 2 ) = (1, 0) indicates the uniaxial buckling while (γ 1 , γ 2 ) = (1, 1) denotes the biaxial buckling. The lowest eigenvalue λ 1 becomes the critical buckling load N cr , and the in-plane loads N 0 x and N 0 y are γ 1 N cr and γ 2 N cr , respectively. Figure 3a shows a simply supported FG-CNTRC plate subject to a uniform distributed load q 0 equal to 1.0 × 10 5 N/m 2 . The width a is 0.1 m and the thickness d is a set variable for the sake of parametric investigation. The matrix is PmPV [9] and its isotropic elastic properties are E m = 2.1 GPa, ν m = 0.34, and ρ m = 1150 kg/m 3 . CNTs have the orthotropic material properties given by E CNT
The CNT efficiency parameter j  for the numerical simulations was chosen by referring to Zhu et al. [20] and is recorded in Table 1. Figure 3b shows a 21 21  uniform NEM grid generated on the midsurface of the plate. The stiffness and mass matrices and load vector were integrated by making use of the 2-D in-plane Gaussian integration using 7 points and the thickness-wise trapezoidal rule using 40 equal segments. The simply supported condition is implemented as follows:   The CNT efficiency parameter η j for the numerical simulations was chosen by referring to Zhu et al. [20] and is recorded in Table 1. Figure 3b shows a 21 × 21 uniform NEM grid generated on the midsurface of the plate. The stiffness and mass matrices and load vector were integrated by making use of the 2-D in-plane Gaussian integration using 7 points and the thickness-wise trapezoidal rule using 40 equal segments. The simply supported condition is implemented as follows: U  The bending analyses were carried out by changing the value of V * CNT , the widththickness ratio a/d, and the CNT distributions. The calibrated central deflections u c z /d are recorded in Table 2 for comparison with those of Zhu et al. [19]. First of all, one can see that the present method is in good agreement with Zhu et al., with the peak relative difference of 1.01% at FG-X for V * CNT = 0.11 and a/d = 10. In addition, the difference between the two methods becomes smaller and is proportional to a/d. This agrees with the fact that all the plate theories approach a certain limit theory (for example, the Kirchhoff theory for isotropic plates) as the thickness decreases [26].   Figure 4a comparatively represents the magnitude of u c z /d for three different values of V * CNT , where the case without reinforcement with CNTs is included. First, it is apparent that the central deflection is significantly reduced by reinforcing only a small amount of CNTs (more than seven times at a/d = 50 when compared with the non-reinforced plate). Second, the magnitude of u c z /d becomes smaller in proportion to V * CNT . Figure 4b represents the effect of the vertical CNT distribution on the magnitude of u c z /d, where FG-X produces the lowest value while FG-O leads to the highest one. FG-V and UD show the second and third highest levels, respectively. This relative variation of u c z /d can be understood from the fact that the bending stiffness is influenced by the vertical distribution pattern of CNTs. FG-X has the highest bending stiffness because CNTs are concentrated on the bottom and top sides, and vice versa for FG-O. This result informs us that the desired mechanical behavior of the FG-CNTRC plate can be obtained when the vertical CNT distribution is suitably designed. d / u z , where FG-X produces the lowest value while FG-O leads to the highest one. FG-V and UD show the second and third highest levels, respectively. This relative variation of d / u c z can be understood from the fact that the bending stiffness is influenced by the vertical distribution pattern of CNTs. FG-X has the highest bending stiffness because CNTs are concentrated on the bottom and top sides, and vice versa for FG-O. This result informs us that the desired mechanical behavior of the FG-CNTRC plate can be obtained when the vertical CNT distribution is suitably designed.    Figure 5a comparatively shows the thickness-wise distributions of calibrated normal stress σ xx = σ xx d 2 / q 0 a 2 to the vertical CNT distribution pattern. The width-thickness ratio a/d and V * CNT are 50 and 0.17, and the stresses were measured along the thickness at point O. UD shows a typical anti-symmetric linear variation, but three different FGs exhibit curved nonlinear variations. It is because the (1,1,0)* hierarchical model leads to liner axial stress distribution through the thickness of the homogeneous material. However, for the inhomogeneous material, the vertical distribution of the axial stress is dependent on the vertical distribution of CNTs. FG-O and FG-X produce anti-symmetric variation, but FG-V shows a leaning parabolic variation. This is caused by the thickness-wise distribution of CNTs, which characterizes the bending stiffness along with the thickness. Figure 5b comparatively shows the thickness-wise distributions of shear stress τ zx = τ zx d 2 / q 0 a 2 at point O. As in the previous Figure 5b, FG-O represents the highest value and FG-X provides the lowest value owing to the difference in the bending stiffness among the four CNTRC plates. Differing from UD, FG-O, and FG-X, FG-V exhibits a leaning shear stress distribution owing to the non-symmetric vertical distribution of CNTs, as in σ xx in Figure 5a. Therefore, this result confirms once again that the behavior of bending deformation and stress is strongly influenced by the thickness-wise distribution of the CNTs.  Figure 5a. Therefore, this result confirms once again that the behavior of bending deformation and stress is strongly influenced by the thickness-wise distribution of the CNTs. , and the CNT distribution pattern. A simply supported square FG-CNTRC plate shown in Figure 4a was taken, using the same material properties for the matrix and CNTs. However, the NEM grid density was changed from Next, the free-vibration behavior of the CNTRC plate was examined under the value of V * CNT , the width-thickness ratio a/d, and the CNT distribution pattern. A simply supported square FG-CNTRC plate shown in Figure 4a was taken, using the same material properties for the matrix and CNTs. However, the NEM grid density was changed from 21 × 21 to 41 × 41 in order to secure the modal analysis accuracy. A total of twenty-one modes were extracted from the Lanczos and Jacobi iterations. In Table 3, the six lowest calibrated natural frequencies ω = ω a 2 /d ρ m /E m are compared with those of Zhu et al. [20]. First, it can be observed that the natural frequencies of modes (2,1) and (2,2) are higher than those of modes (1,3) and (1,4), regardless of V * CNT and the thickness-wise CNT distribution. This mode sequence of FG-CNTRC plates is not the same as that of isotropic plates [34], which is due to the difference in the elastic properties between the x-and y-axes. The elastic properties of the y-axis are smaller than those of the x-axis because CNTs are parallel to the x-direction, as depicted in Figure 1a. Second, the comparison of detailed natural frequencies between the two methods informs us that the peak relative difference equal to 6.07% occurs at V * CNT = 0.14 of FG-V. Thus, the accuracy of the present method has been verified.  Figure 6a comparatively represents the variation in the calibrated fundamental frequencies (CFFs) ω (1,1) to the volume fraction V * CNT . It is found that CFFs increase in proportion to V * CNT because the stiffness increase is superior to the mass increase. This explanation can be justified from the fact that the non-CNTRC plate (i.e., V * CNT = 0) provides a much lower CFF. Figure 6b compares CFFs with respect to the thickness-wise CNT distributions. By referring to the previous Figure 5b for the calibrated central deflection, one can realize that the trend is completely reversed. That is, FG-X exhibits the highest CFF while FG-O provides the lowest one. UD and FG-V show the second and third highest levels, respectively. Thus, it has been justified once again that the vertical distribution of CNTs is important, and FG-X is the stiffest while FG-O is the weakest. The results of Figure 7a,b inform us that the vibration characteristic of the CNTRC plate is markedly influenced by the volume fraction V * CNT and the thickness-wise CNT distribution.
butions. By referring to the previous Figure 5b for the calibrated central deflection, one can realize that the trend is completely reversed. That is, FG-X exhibits the highest CFF while FG-O provides the lowest one. UD and FG-V show the second and third highest levels, respectively. Thus, it has been justified once again that the vertical distribution of CNTs is important, and FG-X is the stiffest while FG-O is the weakest. The results of Figure 7a,b inform us that the vibration characteristic of the CNTRC plate is markedly influ- Next, the critical buckling loads of the FG-CNTRC plate were calculated and compared with the analytical solutions of Chesmeh et al. [24], as given in Table 4. The plate was simply supported, its width-thickness ratio d / a was 100, and three different vertical distributions and volume fractions of CNTs were taken. The critical buckling loads were calibrated as The midsurface of the plate was uniformly discretized by 15 15  grid points, and the lowest six buckling load parameters and buckling modes were solved by Lanczos transformation and Jacobi iteration. It was found that the calibrated buckling loads obtained by the present NEM were slightly higher than the  becomes larger. Regarding the vertical CNT distribution, the highest and lowest levels occur at FG-X and FG-O and the second and third highest levels occur at UD and FG-V, respectively. This trend is almost the same as one of the previous free vibrations, and it is because the buckling load increases in proportion to the plate stiffness, which is influenced by the vertical CNT distribution. This trend is also observed in Figure 7b, which represents the influence of the CNT volume fraction  Next, the critical buckling loads of the FG-CNTRC plate were calculated and compared with the analytical solutions of Chesmeh et al. [24], as given in Table 4. The plate was simply supported, its width-thickness ratio a/d was 100, and three different vertical distributions and volume fractions of CNTs were taken. The critical buckling loads were calibrated as N cr = N cr a 2 /E m d 3 . The midsurface of the plate was uniformly discretized by 15 × 15 grid points, and the lowest six buckling load parameters and buckling modes were solved by Lanczos transformation and Jacobi iteration. It was found that the calibrated buckling loads obtained by the present NEM were slightly higher than the analytical solutions by Chesmeh et al. [24]. This is consistent with the fact that the numerical approximation provides buckling loads that are higher than the analytical solutions. The maximum relative error occurred at the uniaxial for V * CNT = 0.17 of FG-O, and its magnitude was 6.12%. Therefore, the accuracy of the present method has been justified once again.  Figure 7a compares the variations of N cr to the width-thickness ratio a/d for various distributions of CNTs. First of all, one can see that the calibrated buckling loads increase and saturate as a/d becomes larger. Regarding the vertical CNT distribution, the highest and lowest levels occur at FG-X and FG-O and the second and third highest levels occur at UD and FG-V, respectively. This trend is almost the same as one of the previous free vibrations, and it is because the buckling load increases in proportion to the plate stiffness, which is influenced by the vertical CNT distribution. This trend is also observed in Figure 7b, which represents the influence of the CNT volume fraction V * CNT on the calibrated buckling load N cr . It can be seen that N cr uniformly increases in proportion to the value of V * CNT . Therefore, the buckling could be effectively suppressed by increasing V * CNT and/or adopting FG-X. Figure 8a represents the influence of boundary conditions on the variation of N cr to the width-thickness ratio a/d, where S and C denote simply supported and clamped. A combination of four characteristics denotes a combination of boundary conditions that is enforced to sides 1 , 2 , 3 , and 4 shown in Figure 3. The highest and lowest levels are seen at CCCC and SSSS, while SCSC and CSCS provide the curves close to CCCC and SSSS, respectively. It is because two sides ( 2 and 4 ) play an important role in the boundary condition for the FG-CNTRC plates in which CNTs are parallel to the x− direction, and these two sides are clamped at SCSC while simply supported at CSCS. , where S and C denote simply supported and clamped. A combination of four characteristics denotes a combination of boundary conditions that is enforced to sides , , , and  shown in Figure 3. The highest and lowest levels are seen at CCCC and SSSS, while SCSC and CSCS provide the curves close to CCCC and SSSS, respectively. It is because two sides ( and ) play an important role in the boundary condition for the FG-CNTRC plates in which CNTs are parallel to the − x direction, and these two sides are clamped at SCSC while simply supported at CSCS.     One can see that the buckling load increases and saturates as the aspect ratio a/b becomes larger because the plate becomes stiffer as the constrained plate becomes narrower with the side length a kept the same. The difference in N cr between boundary conditions becomes insensitive as the aspect ratio increases. This is because the width of the two dominant sides 2 and 4 becomes shorter as the aspect ratio increases so that the type of boundary condition specified for these sides does not lead to a marked difference.

Conclusions
The bending, free-vibration, and buckling behaviors of FG-CNTRC plates were investigated by 2-D NEM. Based on the hierarchical model, the displacement was split into the 2-D in-plane vector functions and the assumed thickness monomials. The in-plane displacement part was solved by the 2-D natural element method. The numerical experiments were performed to validate the introduced method and to explore the bending, free-vibration, and buckling characteristics of the FG-CNTRC plates. The following observations were drawn from the numerical results:

•
The accuracy of the present method is justified with the peak relative differences between the proposed method and the references, which are 1.01% for the central deflection, 6.07% for the natural frequencies, and 6.12% for the buckling loads.

•
The central deflection is significantly reduced by introducing only a small amount of CNTs, and its reduction is proportional to V * CNT . The natural frequencies and buckling loads show a completely reverse trend. • The (2,1) and (2,2) modes show higher frequencies than the (1,3) and (1,4) modes because CNTs are parallel to the x-direction. In other words, the stiffness in the x-direction is much higher than that in the y-direction. • UD, FG-O, and FG-X produce anti-symmetric normal stress distributions and symmetric shear stress distributions; FG-V does not show such distribution but rather a leaning parabolic distribution.

•
The buckling load is sensitive to the type of boundary condition, and it increases and saturates in proportion to the aspect ratio. The difference in the buckling load between boundary conditions becomes insensitive in proportion to the aspect ratio.

•
The comparison of the magnitudes of the central deflection, fundamental frequency, and buckling load shows that the order of stiffness of the FG-CNTRC plates is FG-X > UD > FG-V > FG-O.
The current study was conducted from the macroscopic point of view. However, the microscopic factor such as the microstructural evolution would be important, and this represents a topic that deserves future work.