Combined effect of carbon nanotubes distribution and orientation on functionally graded nanocomposite beams using finite element analysis

Structural tailoring can provide a promising performance for Functionally Graded (FG) components in engineering. Moreover, utilizing advanced Carbon Nanotube (CNT) as embedded reinforcement in nanocomposite structures, excellent mechanical properties can be tailored and designed to meet requirements. This research addressed the issue of a particular effect for CNT orientation and gradation distribution on static and free vibration analysis of Functionally Graded CNT-Reinforced Composite (FG-CNTRC) beams. First, an efficient finite beam element capable of controlling both parameters was derived based on the Timoshenko beam theory. Single-Walled CNT (SWCNT) was used as primary reinforcement and graded through-thickness. Then, an extensive parametric study was done for model convergence, static, and dynamic analysis. The proposed model offers unique shape function depends on material properties and cross-section geometry, high-accuracy, and expanded to cover both orientations and grading exponents. This expansion allows passive-control of the beam stiffness and strength without any increment in structural weight. Wherein constituent materials quantities and volume fractions were not changed. Finally, obtained findings concerned about orientation angle and power-law exponent, which showed that they significantly affect the structural response, and therefore offer a practical approach of structure tailoring for applied loads, required response, and specific weight limitations.


CNT m 12
Poisson's ratio of CNT  Young's modulus for CNT and matrix f t Transverse distributed load function along the length G 12 ,G 13 ,G 23 Composite shear moduli G G ,

Introduction
The first observation for CNT was in 1991 by Iijima [1]. From that moment, numerous investigations have been done. It concluded significant physical and mechanical properties of the new carbon form. CNTs stiffness and strength are higher than any current materials [2]. In addition, low specific weight introduced it to be marvelous reinforcement in nanocomposite structures, and to be a replacement for conventional fiber composites on aircraft structures, space vehicles, and other engineering applications [2][3][4][5]. CNTs were extensively studied regarding the functionally graded distribution concept, which significantly enhanced the mechanical behavior of CNT reinforced composite structures, during these years [6]. On the other hand, the main benefit of using Functionally Graded Material (FGM) instead of traditional is that the internal composition of their component materials can be tailored to satisfy the requirements of a given structure. Therefore, it significantly enhanced the mechanical properties of these proposed structures. FGM structures were extensively investigated in many kinds of research such as [7][8][9][10][11][12].
In this paper, the mathematical and finite element derivation was performed for analysis of FG-CNTRC beams based on the Timoshenko Beam Theory (TBT) and under the consideration of three different distribution patterns. Contrary to the previous studies [19, 22-25, 27, 30, 31], the new added value of the presented model was that both CNTs orientation angles and the gradient power-law factor were variable. The modified rule of mixture was employed to obtain the mechanical properties. The governing differential equations were derived after applying Hamilton's principle for energy formulas. Then, displacement fields were assumed, and the corresponding shape function was derived, in turn, stiffness and mass matrices. Convergence checks were performed for different responses and inputs to ensure the efficiency of the proposed model. Moreover, extensive parametric studies were carried out for static and free vibration response, including results validation with published ones. Result analysis illustrated the effects of CNT orientation angles and power-law index on the beam responses and showed that the proposed gradation and orientation adjustment allow structural tailoring of beam structures for enhancement of mechanical properties such as stiffness and strength.

FG-CNTRCs approach
Three types of aligned CNT, namely, Uniformly Distributed (UD-CNTRC), FG type-A (FGA-CNTRC), and FG type-X (FGX-CNTRC) for SWCNT inside the polymer matrix, were considered as shown in figure 1. The investigations focused on FGA with variable power-law exponent n, where the other two types UD and FGX, were used for validation with literature. Geometric characteristics and coordinate systems of the beam were also demonstrated in figure 2.  -For FGX-CNT  Figure 3 illustrates the V CNT distribution through the thickness as equations (1), (2), and (3). The composite properties were found using the extended rule of mixture [23]. The young's modulus, shear modulus, composite density, and Poisson's ratio were obtained as equation (4).

Theoretical and finite element formulations
Based on the TBT, equations of motion formulation were developed. The element was selected as a two-node shown in figure 4 based on [60].

Theoretical formulation
First, displacement and strain fields were introduced. The axial and transverse displacement fields were expressed as equation (5). The Strain-displacement field was deduced as equation (6).
Where i, j=1:6 and k, l=1:3 The plane stress assumptions were applied by setting σ 33 and eliminate the corresponding strain ε 33 from the constitutive relation. Equation (7) as one layer was expressed by equation (8). For FG-CNTRC material According to CNT orientation, the transformation from principle axes (1,2, and 3) to geometrical axes (x, y, and z) was applied. Thus, the transformed stress-strain relationship in equation (8) was written as equation (10).¯¯( Then, substituting by σ yy =σ yz =σ xy =γ yz =γ xy =0 and ε yy ≠0 as [62] y-direction was stress-free and the plane stress assumptions were used. After reduction, equation (10) was rewritten in equation (11).
The reduced transformed stiffness components took the form of equation (12).¯¯˜¯( Introducing the following stiffness and mass coefficients as [7] in equation (13) and (14). Third, using the above expressions, an energy formulation was deduced. The virtual kinetic-energy of the beam structure was stated by [63]. Virtual internal strain-energy was given by [64]. The final form of virtual energies were derived as equation (15) and (16).
Finally, zero external work was assumed, and Hamilton's principle was applied to obtain the governing differential equations in terms of three degrees of freedom u 0 , w 0 and f [11] as followed in equation (17).̈( )

Finite element formulations
Based on the static part solution of the differential equations (17), displacements were assumed as polynomial functions. For u 0 and f, the order chooses to be from the second rank and third for w 0 . The exact solution for the displacement has ten constants, and only six independent were needed. The four dependent can be found as a function of the independent six. So, another coupling stiffness coefficients were introduced as equation (18). And the displacements were stated as equation (19). Equation (19) was rewritten in matrix form in equation (20).
The exact shape function [ù(x)] was given in appendix A. Stiffness matrix was derived as equation (21).
From the dynamic part of the governing differential equations, the mass matrix was deduced as the sum of the following four sub-matrices as equation (22).
Distributed load as a force element was expressed as equation (23).
Stiffness, mass matrix, and distributed load vector were listed in appendix A, after performing matrices multiplications, summations and integration stated previously. Finally, obtained matrices were substituted in the governing equation of motion in equation (24).

Numerical results
The presented element was validated by comparing the results with two previously published data based on nonlinear studies and high order theories. A parametric study for FG-CNTRC beams was performed. The study includes model convergence, natural frequencies, static deformation, and stress distributions analysis. A MATLAB code was conducted to perform the following analysis using the presented finite element model. The boundary conditions of the beams were considered as followed. The geometric characteristics were stated in each separate case with unit width and specific length to thickness ratio (L/h). The shear correction factor was taken as k s =5/6. Material properties were chosen as [22,65] and listed in table 1. The efficiency parameters for CNT/matrix were obtained by comparing the elastic properties result of molecular dynamics simulations with the ones calculated by extended rule of mixture for specific CNT total volume fraction V tcnt . Table 2 represents the selected values and their corresponding efficiency parameters as [22,66,67].

Model convergence
The convergence of the presented model was checked for both static deformation and natural frequency. The geometric characteristics were L/h=25, CNT orientation angles θ=0°,90°, and V tcnt =0.12.
For static convergence tests, cantilever CNTRC beams were subjected to uniformly distributed load of intensity 1N/m, and the normalized tip deflection was given as [60] in equation (25). S-S CNTRC beams were chosen for the natural frequency convergence test, where, the non-dimensional natural frequency was calculated using equation (26).¯( It can be seen from figure 5 that the model prediction was converged at only one element with super-convergence property for all proposed CNT distributions. This feature occurs because the stiffness matrix was exact for point loading as the displacement field was exactly derived from the governing static differential equations homogeneous part. In contrast, the convergence of the eigenvalues started at around n≈10, as in figure 6, because the mass distribution was approximated and was not the exact solution for the dynamic part of governing partial differential equations. It was clear from previous that the present model was converged in a suitable number of elements and free of shear-locking or divergence.

Dynamic analysis
Numerical results were tabulated to investigate the effect of CNT volume fraction, CNT orientation angle, CNT distributions, boundary conditions, and slenderness ratios on the non-dimensional natural frequencies of FG-CNTRC beams. Moreover, validation of free vibration was performed by comparing the results with published ones. For validation, first, three non-dimensional natural frequencies for UD, FGA, and FGX-CNTRC beams were compared with the results obtained from [22] and [68] in table 3. Ansari et al [22] performed a nonlinear study based on TBT, and Shen et al [68] also represented a nonlinear study based on high-order beam theories. It can be seen that results have an acceptable agreement with the literature depends on nonlinear and high order theories analysis. Table 4 also compared with the results obtained by [22]. It also studies the effect of boundary condition and slenderness ratio change on the accuracy of the present model. Noted that the non-dimensional natural frequencies of all the following results were calculated as equation (27).   To study the effect of CNT orientation angle change on the natural frequencies, table 6 shows that increasing of CNT angle lowered the natural frequency values. It can also be noticed, huge reducing jump from θ=0°to θ=45°in comparison with the reduction from θ=45°to θ=90°. This observation strictly affects the tailoring process using the CNT orientation angles.
For FGA-CNTRC beams, a suitable range of power-law gradation index n from 0.3 to 3 as [11] has been chosen to avoid over single-material domination. Figure 7 shows the effect of the variation of power-law exponent on the first three non-dimensional natural frequencies for the S-S FGA-CNTRC beam. By increasing n, it decreased with exceptional values. This result indicates an essential application for using FGA-CNTRC beams in response control by tailoring the gradation index.

Static analysis
Static response and stress behavior were parametrically studied for the presented three types of CNT distributions. The analysis starts with an investigation of the boundary conditions and CNT orientation angles' effect on the static response. Then, the effect of change of power-law exponent n on the static response for the FGA-CNTRC beam was performed. Finally, extensive stress analysis for all CNT distributions was conducted. Table 7 tabulated the normalized transverse deflection of CNTRC beams according to equation (26) due to the applied unit distributed load. Various boundary conditions and CNT orientation angles were used. The normalized deflections were calculated at the mid-point of each case. The overall lower deflection values appear on FGX, then UD and FGA, respectively. For all proposed CNT distributions, the deflection increased by changing the boundary conditions from C-C to C-S, then S-S, and C-F, respectively. The most critical notice was about the CNT orientation angle effect on the value of transverse deflection. It can see from the results, the significant change of deflection values exists between θ=0°to θ=45°also, the relatively tiny values change between θ=45°to θ=90°. The same notice was observed in natural frequency analysis. This notice could open new ways to response control in addition to dynamic properties tailoring, especially on response constraint design problems. The normalized mid-point and tip deflections for the C-F FGA beam were drawn in figure 8 against the power-law exponent change. The deflection increased by the increment of exponent value with logic slope difference between mid-point and tip deflection. The percentage of deflection change reaches around 100% increment between n=0.3 to n=3. This evidence supports the ability of response control in FGA-CNTRC beams using power exponent control, with the same total CNT volume fraction and growth weight.
For stress analysis, figure 9(a) shows normal stress configurations for UD, FGA, and FGX-CNTRC cantilever beams. The UD represents the most straightforward linear conventional normal stress distribution, the same upper and lower value with an inverse sign to indicate compression and tension sides. FGA with n=1 obtains non-conventional normal stress with the highest tension values among all other distributions at the lower surface, and the smallest compression one at the upper side. In between, quadratic distribution appears to influence the neutral plan position, which moves downward in the case of n=1. FGX did not influence the neutral plan position, but cubic distribution appears and smoothes UD response with larger maximum stress values at both compression/tension sides. Figure 9(b) represents shear stress for the proposed CNT distributions. TBT assumption of constant shear stress along thickness appears in the UD case. The corrected shape had been drawn for the UD case to show the real shear stress shape and zero shears at upper and lower surfaces, to explain the represented error due to this constant shear stress assumption. However, it cannot be precisely predicted in FGA, and FGX cases, the obtained  results were acceptable to understand the stress manner across the thickness. Both were not straight lines with the same shape as CNT volume fraction distribution through the normalized thickness figure 3. Stress behavior was changed from the CNT orientation point of view. Figure 10(a) illustrates the change in normal stress shape due to the CNT orientation angle change. It can be seen, a relatively small change appears between θ=45°to θ=90°in comparison with the quadratic distribution θ=0°in addition to inverse increasing and decreasing values at maximum compression/tension on the upper/lower surfaces, respectively. However, shear stress configures a different manner. Figure 10(b) represents shear stress values variation with θ. As logic results, there was a very tiny difference in shear at θ=0°and θ=90°, but the difference only appears at the range between them.

Results analysis
CNT orientations θ and power-law exponent n were proved to be the significant factors affecting the beam response in the previous section 4.
Further investigations were performed to study the individual and mutual effect on the FGA-CNTRC beam static and free vibration response. Other parameters were chosen to be fixed as the following, C-F beam, f t =1 N/M, L/h=25, and V tcnt =0.12. The power-law index margin was selected from 0.3 to 3 and CNT angles from 0°to 90 o .  First, the effect on normalized first natural frequency was studied where the normalization followed equation (27) and multiplied by factor equals to 100. Figure 11 shows that the orientation angle had greater influence than grading exponent on the natural frequency, and the combination of both parameters capable of tailoring the frequency around 400% from the lowest value. Second, the effect on the normalized deflection was investigated at x=0 or the mid-point of the cantilever beam and the normalization as equation (25). Figure 12 illustrates the same notice about orientation angle but increasing of tailoring ability in case of deflection control to 700%.
Next, absolute normal stresses at the top-most and bottom-most point at x=−L/2 were explored according to the same criteria. Figure 13(a) represented normal stress at the upper surface with a higher effect of θ than n and control-ability reached 800%. Nevertheless, figure 13(b) showed the great influence of both parameters on the lower surface normal stress and variation equal to 300%.
Finally, shear stresses also graphically investigated in figure 14. Mid, upper, and lower points were highly distorted due to the error obtained from the constant shear stress assumption of TBT. Therefore, intermediate point values at h=25% & h=−25% were chosen. The figure presented a similar and upset shape for the two cases with the great effect of both orientation and grading index. The change in shear stress value showed in figure 14(a) equals to 500%, and in figure 14(b) equals 54%.

Conclusion
The motivation of this research was to provide the idea of tailoring FG-CNTRC beams for static and free vibration responses using controlled-orientation and grading values. Mathematical and finite element modeling was  performed, followed by the governing equation of motions, shape function, stiffness, and mass matrices derivations. Parametric studies were conducted and hoped that corresponding results could be useful in future studies. The endsection illustrated graphically the real contributions of major parameters θ and n on different responses obtained.
The proposed model offers a unique shape function, high-accuracy, expanded to cover CNT orientation angle, and gradation index. It has a super-convergence property for static analysis and good predicting convergence for   free vibration studies. The obtained shape function depends on both material properties and cross-section geometry. High accuracy was achieved after validation by comparing the results to published ones. Moreover, founded results accuracy reaches high order shear deformation theories and nonlinear studies findings.
The parametric study investigated including boundary conditions, CNT volume fraction, various CNT distribution, power-law exponent for FGA-CNTRC, slenderness ratio, and CNT orientation angles. The results indicated that the CNT orientation and power-law index were the main critical factors controlling the different responses. This observation strictly affects the tailoring process of the beam stiffness and strength. It could open new ways to response control in addition to properties tailoring without any increment in the structure weight. For example, natural frequency, static deflection, normal stress, and shear stress can be controlled by at least 400%, 700%, 300%, and 54%, respectively. Future work will include axially graded properties for CNT. Also, the tailoring process will be studied in a real design problem to conceptualize the whole approach.

Data Availability
The raw/processed data required to reproduce these findings cannot be shared at this time as the data also forms part of an ongoing study.

Appendix
Exact shape function   Mass matrix elements  L  I  I L  IL  M  M  L  I  I L  IL   M  M  L  I L  I L  IL  I L  I L  I   M  M  L  IL  IL  IL  IL  IL  I  I  I   M  L  IL  IL  IL  IL  IL  I   M  M  L  IL  IL  IL  IL  IL  I  I 2  3  ,  2  48  3  12  ,  24  3  12  13  42  84  42  294  1680  35   11  126  147  21  231  1260  1260  1260  210  3  3  28  56  28  84  560