Free Vibration Responses of Functionally Graded CNT-Reinforced Composite Conical Shell Panels

Functionally graded CNT (carbon nanotube)-reinforced composites (FG-CNTRCs) are intensively studied because the mechanical behaviors of conventional composites can be dramatically improved. Only a small amount of CNTs are appropriately distributed through the thickness. However, the studies on conical shell panels have been poorly reported when compared with beams, plates and cylindrical shells, even though more parameters are associated with the mechanical behaviors of conical shell panels. In this context, this study intends to profoundly investigate the free vibration of FG-CNTRC conical shell panels by developing an effective and reliable 2-D (two-dimensional) numerical method. The displacement field is expressed using the first-order shear deformation shell theory, and it is approximated by the 2-D planar natural element method (NEM). The conical shell surface is transformed into the 2-D planar NEM grid, and the approach for MITC3+shell element is employed to suppress the shear locking. The developed numerical method is validated through the benchmark experiments, and the free vibration responses of FG-CNTRC conical shell panels are investigated with respect to all the associated parameters. It is found from the numerical results that the free vibration of FG-CNTRC conical shell panels is significantly influenced by the volume fraction and distribution pattern of CNTs, the geometry parameters of the conical shell, and the boundary condition.


Introduction
Carbon nanotube (CNT) has attracted great attention as a promising material for the 21st century thanks to its excellent physical and chemical properties, such as its high mass-strength ratio. Since its introduction, it has been widely adopted as a next-generation pillar for polymer composites because the addition of a small amount of CNTs dramatically increases the stiffness of CNT-reinforced composites called CNTRC [1,2]. For use as mechanical and structural members, CNTRCs are usually developed in the forms of beams, plates and shells. Later, the thickness-wise distribution pattern of CNTs was purposefully assumed [1,3] according to the concept of functionally graded materials (FGMs) [4]. In FGMs, the thickness-wise volume fractions of base materials vary continuously and may be optimally tailored so as to maximize the target function [5]. The CNTRCs with functionally graded CNT distributions through the thickness are called FG-CNTRCs, and several representative ones, such as FG-U, FG-V, FG-O, FG-X and FG-Λ, were introduced. The purpose of these FG-CNTRCs was to allow one to choose one among them which is most suitable for the target performance, and it requires a profound investigation of the mechanical characteristics of each FG-CNTRC. In this context, extensive research efforts have focused on the investigation of mechanical responses such as static bending, free vibration and buckling. The research works on CNTRCs and FG-CNTRCs would be referred to the references [6][7][8].
troublesome shear locking for the bending-dominated thin structures [28]. The developed NEM-based numerical method is verified through the benchmark experiments, and the free vibration responses of FG-CNTRC conical shell panels are investigated with respect to the CNT-associated parameters, the conical shell parameters, the sandwich core thickness, and the boundary conditions. This paper is organized as follows. Following the introduction, the FG-CNTRC conical shell panel and its displacement, strain and stress fields are addressed in Section 2. The natural element approximation of free vibration of FG-CNTRC conical shell panels is explained in Section 3. The numerical results of benchmark and parametric experiments are presented and discussed in Section 4, and the final conclusion is given in Section. Figure 1a represents a circular conical shell panel reinforced with single-walled carbon nanotubes (SWCNTs). A coordinate system (x, θ, z) is positioned on the mid-surface of the panel, and the geometric configuration of the conical shell panel is characterized by semi-vertex angle α, sub-tended angle θ 0 , length , small radius R 0 and thickness h. The radius R(x) of conical shell panel at any point along its length is determined by

Modeling of FG-CNTRC Conical Shell Panel
CNTs are aligned along the x−axis with a specific functionally graded distribution pattern through the thickness. Four representative distribution patterns of CNTs are represented in Figure 1b, where CNTs are uniformly distributed in FG-U while CNTs are rich at the mid-surface in FG-O, rich at the outer shell in FG-X and rich at the inner shell in FG-Λ, respectively. The volume fractions of CNTs and underlying matrix are denoted by V CNT (z) and V m (z) which are in the relation given by where the CNT volume fractions V CNT (z) corresponding to the four representative CNT distribution patterns are expressed as with Here, w CNT is the mass fraction of CNTs occupied within a unit volume of a conical shell panel, and ρ CNT and ρ m denote the CNT density and the matrix density. As a dual-phase composite material, FG-CNTRC structures are usually modeled as an orthotropic material, and their effective elastic properties are evaluated using the homogenization methods such as the liner rule of mixtures by introducing the directionwise CNT efficiency parameters η j (j = 1, 2, 3). The effective elastic and shear moduli of the FG-CNTRC conical shell panel are calculated according to [3]: Polymers 2023, 15, 1987 4 of 16 It is assumed that G 13 = G 12 and G 23 = 1.2G 12 . Meanwhile, the effective density ρ and the effective Poisson's ratio ν 12 through the thickness are assumed as It is assumed that 12   By adopting the first-order shear deformation shell theory, the displacement field being the vector of displacement components at the mid-surface of the shell panel. Then, the strain-displacement relations based on the co-  By adopting the first-order shear deformation shell theory, the displacement field u = {u, v, w} T is expressed as with d = u 0 , v 0 , w 0 , ϑ x , ϑ y T being the vector of displacement components at the midsurface of the shell panel. Then, the strain-displacement relations based on the coordinate system (x, θ, z) give Here, H and H s are (3 × 5) and (2 × 5) gradient-like matrices defined by Polymers 2023, 15, 1987 5 of 16 with H x = ∂/∂x, H θ = ∂/R∂θ, RS = sinα/R and RC = cosα/R. Then, the constitutive equations are expressed as with

Analysis of Free-Vibration Using 2-D NEM
For the free vibration approximation of the FG-CNTRC conical shell panel by 2-D NEM, the mid-surface is discretized into a finite number of nodes and Delaunay triangles, as shown in Figure 2. A coordinate (x, s) is adopted to identify the position on the midsurface using the relation of s = Rθ, and Laplace interpolation functions ψ J (x, s) [25,26] are assigned to each node J. Letting the total number of nodes be N, the NEM approximation with d J = u 0 , v 0 , w 0 , ϑ x , ϑ y T J being the nodal vector of displacement components at node J.
Referring to the references [25,26], the definition of the Laplace interpolation function and its computation on the conical surface is complex and painstaking. To overcome this difficulty, the physical NEM grid C = [0, ] × [−s L , s R ] with M triangles on the conical surface is mapped to the computational NEM grid R = [0, ] × [−θ 0 /2, θ 0 /2] on the rectangular plane according to the inverse of the geometry transformation T C defined by where Then, Laplace interpolation functions ψ J (x, s) are mapped to ϕ J (ζ 1 , ζ 2 ), and one can obtain the inverse Jacobi matrix J −1 given by for the above geometry transformation. Moreover, the partial derivatives H x and H θ in Equations (13) and (14) are switched to according to the chain rule.  (22) according to the chain rule. By substituting Equations (21) and (22) into Equations (13) and (14) (23) according to the inverse transformation  By substituting Equations (21) and (22) into Equations (13) and (14), one can obtainĤ andĤ s in terms of H 1 and H 2 according to the inverse transformation T −1 C . Then, the NEM approximations of the bending-membrane strain ε in Equation (11) and the transverse shear strain γ in Equation (12) end up with The transverse shear strains in Equation (25), which are directly computed from the standard C 0 −displacement approximation, suffer from shear locking [24,28,29]. To suppress the shear locking, the transverse shear strains are interpolated by employing the concept of the 3-node triangular MITC3+shell finite element depicted in Figure 3. Each element e in the physical NEM grid C is to be mapped to the master triangular element . Moreover, the displacement field is re-interpolated within each element with d e K = u e 0 , v e 0 , w e 0 , ϑ e y , ϑ e y T K being the local nodal vector within e and N K (ξ, η) being the linear triangular FE shape functions [30]. Then, the element-wise transverse shear strains are interpolated aŝ Polymers 2023, 15, 1987 7 of 16 using the values at the typing points andĉ = γ xz . The analytical manipulation of Equations (27) and (28) using Equations (12) and (26), together with the chain rule between two coordinates (s, x) and (ξ, η), leads tô withB e being the (2 × 15) matrices in terms of ξ, η, z, α and R and d e = {d e 1 , d e 2 , d e 3 } being the (15 × 1) element-wise nodal vector.
using the values at the typing points and The ana nipulation of Equations (27) and (28) using Equations (12) and (26)  Meanwhile, for the free vibration analysis of FG-CNTRC conical shell standard Galerkin weak form can be derived from the dynamic form of energ [31] ( ) ( where m is the ( ) Meanwhile, for the free vibration analysis of FG-CNTRC conical shell panels, the standard Galerkin weak form can be derived from the dynamic form of energy principle [31] where m is the (5 × 5) symmetric matrix defined by with the (3 × 3) identity matrix I and m 2 = diag z 2 , z 2 . Substituting Equations (24) and (28), together with the constitutive relations Equations (15) and (16), into Equation (29) results in the modal equations given by . Here, the stiffness and mass matrices are defined by with the shear correction factor κ = 5/6, the length L e of the longest sides of Delaunay triangles, and a positive constant α called the shear stabilization parameter [21,24]. The value of 5/6 is taken for the shear correction factor because the displacement field is expressed by the FSDT. The value of α is recommended to be 0.1 [21] and taken through the preliminary experiment. This modification of shear modulus was proposed for the sake of further stabilization of the MITC3 element. The numerical integration of three matrices is performed triangle by triangle, and 7 Gauss integration points are used for K σ and M while 1 Gauss point is used for K e s . A flowchart for the free vibration analysis by the present 2-D NEM-based method is represented in Figure 4.   Meanwhile, three types of boundary conditions, simply-supported (S), clamped (c) and free, are given, where S and C are given as at side 1 shown in Figure 2, for example. It is noted that the simply-supported condition (37) implies that the shell panel is movable [32].

Verification
The numerical formulae given in Section 3 were coded in Fortran and integrated into the NEM program [2,26], which was previously developed for plate-like structures. The present method was compared with the other reference methods for two completely clamped isotropic and one FG-CNTRC conical shell panels, which are represented in Figure 5. The geometry dimensions of the first example are R 0 = 1.0 m, = 3.0 m, h = 0.03 m, α = π/6 and θ 0 = π/3, and the elastic modulus, Poisson's ratio, and the density are E = 2.1 GPa, ν = 0.3 and ρ = 1150 kg/m 3 , respectively. The natural frequencies were normalized asω = ω 2 ρh/D with D = Eh 3 /12 1 − ν 2 being the flexural rigidity.

Verification
The numerical formulae given in Section 3 were coded in Fortran and integrated into the NEM program [2,26], which was previously developed for plate-like structures. The present method was compared with the other reference methods for two completely clamped isotropic and one FG-CNTRC conical shell panels, which are represented in The lowest four natural frequencies were computed for different NEM grid densities and recorded in Table 1. Referring to Figure 2, the grid density n m × indicates that the total node numbers are m and m along the 2 ζ -and − 2 ζ axes. It is observed that the four normalized natural frequencies uniformly decrease and show a stable convergence with a relative difference % / 100 × = ω ω ∆ γ less than 0.3%. According to this convergence experiment, the NEM grids with equal or similar density to

21×
were used for the remaining numerical experiments.  The lowest four natural frequencies were computed for different NEM grid densities and recorded in Table 1. Referring to Figure 2, the grid density m × n indicates that the total node numbers are m and m along the ζ 2 -and ζ 2 −axes. It is observed that the four normalized natural frequencies uniformly decrease and show a stable convergence with a relative difference γ = |∆ω|/ω × 100% less than 0.3%. According to this convergence experiment, the NEM grids with equal or similar density to 21 × 21 were used for the remaining numerical experiments. In Table 2, the lowest four normalized natural frequencies at the grid density 21 × 21 are compared with the reference solutions obtained by the numerical methods. Where the term without locking suppression indicates the numerical results obtained by the present method in which the transverse shear stiffness matrix K e s in Equation (34) is obtained using the full integration without using MITC3+shell elements. One can see that big errors occur in the free vibration analysis when the locking is not suppressed, such that the four normalized natural frequencies are remarkably larger than the reference solutions. On the other hand, the present method using the MITC approach shows a good agreement with three reference solutions such that the maximum relative difference is 2.61% at mode 4. Here, Bardell et al. [33] used the h-p finite element method, Au and Cheung applied the finite strip method, while Xiang et al. [17] used the element-free kp-Ritz method, respectively. These three methods used curved meshes, but the present method uses a 2-D planar grid, which is easier to implement. The second example is the completely isotopic conical shell panel in which the semivertex and sub-tended angles α, θ 0 and the thickness h are smaller than those of the first example. In other words, the second example is close to a clamped very thin plate-like structure. The geometry dimensions are R 0 = 1.0 m, /s = 0.2, s/h = 1000,α = 7.5 0 and θ 0 = 20 0 . The material properties are the same as in the first example. The parameter b in the normalization factor is represented in Figure 1. It is clearly seen that the present method is in good agreement with the other three numerical methods, such that the maximum relative difference is 1.53% at mode 3. However, the case without locking suppression shows big errors such that the four normalized natural frequencies are larger than the reference solutions. However, the errors are shown to be smaller than those in Table 2 because the membrane deformation decreases as the sub-tended angle θ 0 becomes smaller. In other words, the shell panel becomes a plate-like structure in reverse proportion to θ 0 as represented in Figure 5b, so that the membrane-induced locking becomes less severe. Here, Lim and Liew [35] employed the Ritz method using the pb-2 shape functions to overcome the mathematical complexity stemming from the shell geometry, similar to the present method. However, differing from the present method, they integrated the whole shell as a single element using the pb-2 shape functions.
It has been justified from Tables 2 and 3 that the present method does not suffer from shear locking. Moreover, it is worth noting that the present results were obtained using 2-D planar, not 2-D curved, and coarse NEM grids. The detailed comparison showing the superiority of NEM against FEM at coarse grids may be referred to in the previous work [36]. The third example is FG-CNTRC conical shell panels manufactured with poly (methyl methacrylate) (PMMA) reinforced with (10,10) SWCNTs. The material properties of the two materials are presented in Table 4, and the efficiency parameters η j (j = 1, 2, 3) introduced in Equations (5)-(7) for PMMA/CNT are given in Table 5, respectively. These material properties and parameters are referred to the rule of mixture results at the temperature of T = 300K, which were reported by Han and Elliott [37] and Shen and Xiang [38]. The geometry dimensions are R 0 = 0.2 m, = 0.8 m, α = 45 0 and θ 0 = 120 0 , and the volume fraction V * cnt of CNTs is set by 0.12. The natural frequencies were normalized aŝ ω = ω 2 ρ m h/E m /(2πh) with ρ m and E m being the material properties of PMMA. The fundamental frequencies were computed for three different thickness ratios R 0 /h, four different CNT distribution patterns and two different boundary conditions. In Table 6, the results are compared with the reference solutions, which were obtained by Xiang et al. [17] using the element-free kp-Ritz method, where SSSS and CCCC indicate that all sides 1 ~4 of the shell panel are simply supported and clamped, respectively. It is found that both methods are in good agreement, such that the maximum relative difference is 4.75% at FG-Λ for R 0 /h = 25 and SSSS. From the detailed comparison, it is found that the normalized first frequencies obtained by the present method are, as a whole, higher than those of the Ritz method for SSSS and vice versa for CCCC. In addition, it is seen that the relative differences for SSSS are, as a whole larger than those for CCCC.

Parametric Investigation
Next, the free vibration responses of FG-CNTRC conical shell panels are parametrically investigated with respect to the major parameters. The example shown in Figure 5c is taken without changing the material properties, but the volume fraction V * cnt of CNTs and the geometry dimensions of the shell panel are changed depending on the simulation case. The effect of the thickness ratio R 0 /h onω 1 is firstly investigated. It is seen from Figure 6a,b thatω 1 uniformly increases in proportion to R 0 /h, regardless of the volume fraction and the distribution pattern of CNTs. It is entirely owing to the calibration with the thickness h, so, actually, ω 1 decreases with R 0 /h because the stiffness of the shell panel decreases as the thickness h becomes smaller. Meanwhile,ω 1 uniformly increases with increasing the value of V * cnt because the stiffness increase due to the increase of V * cnt gives rise to more effect onω 1 than the mass increase. This explanation is manifest from the fact thatω 1 of the non-CNTRC shell panel is much lower. It is found from Figure 6b thatω 1 is also influenced by the CNT distribution pattern such that the magnitude order ofω 1 is FG-X, FG-U, FG-Λ and FG-O. This is because the stiffness of the shell panel is also affected by the CNT distribution pattern.
rically investigated with respect to the major parameters. The example shown in Figure   5c is taken without changing the material properties, but the volume fraction   Figure 7a,b represent the variations ofω 1 to the semi-vertex angle α for different CNT distribution patterns and different boundary conditions, respectively. Where CSCS indicates that two sides 1 and 3 are clamped while the other two sides 2 and 4 are simply supported, by referring to Figure 2. It is observed thatω 1 uniformly decreases with increasing the value of the semi-vertex angle, regardless of the CNT distribution pattern and the boundary condition. The radius increase of the shell along the x−axis becomes more apparent as α increases, and accordingly, the stiffness of the shell panel decreases in proportion to the semi-vertex angle. Meanwhile, it is found from Figure 7b thatω 1 is remarkably influenced by the boundary condition such that the magnitude order ofω 1 is CCCC, CSCS, SSSS and CFCF. This order is consistent with the order of boundary constraint. Figure 8a,b represent the variations ofω 1 to the sub-tended angle θ 0 for different thickness ratios and different boundary conditions. All the plots in the two figures show noticeable fluctuations with respect to θ 0 . It was caused by the sensitivity of numerical natural frequency to the modification of the NEM grid owing to the change in the circular shell side length according to the change of θ 0 . Moreover, this fluctuation is also reported in the numerical results of Xiang et al. [17]. It is seen from Figure 8a that the variation of ω 1 is not remarkable until θ 0 decreases to near 30 0 , but thereafterω 1 abruptly goes up in reverse proportion to θ 0 , regardless of the thickness ratio R 0 /h. Moreover, the remarkable difference inω 1 between the thickness ratios until near θ 0 = 30 0 becomes negligible as θ 0 decreases. It is because conical shell panels approach narrow plate-like structures of zero curvature as θ 0 tends to zero, so the fundamental frequency increases owing to the increase of structural stiffness. And the calibrated fundamental frequencies of plate-like structures with different thicknesses approach the same limit as the thickness tends to zero [34,39]. The limit value depends on the boundary condition, as represented in Figure 8b, such that it becomes larger in proportion to the constraint strength of boundary condition. It is worth noting that the variation ofω 1 is not sensitive to θ 0 for CFCF, because two sides 2 and 4 which are relatively long are not constrained.
CSCS indicates that two sides ① and ③ are clamped while the other two sides ② and ④ are simply supported, by referring to Figure 2. It is observed that 1 ω uniformly decreases with increasing the value of the semi-vertex angle, regardless of the CNT distribution pattern and the boundary condition. The radius increase of the shell along the − x axis becomes more apparent as α increases, and accordingly, the stiffness of the shell panel decreases in proportion to the semi-vertex angle. Meanwhile, it is found from  θ tends to zero, so the fundamental frequency increases owing to the increase of structural stiffness. And the calibrated fundamental frequencies of plate-like structures with different thicknesses approach the same limit as the thickness tends to zero [34,39]. The limit value depends on the boundary condition, as represented in Figure 8b, such that it becomes larger in proportion to the constraint strength of boundary condition. It is worth noting that the variation of 1 ω is not sensitive to 0 θ for CFCF, because two sides ② and ④ which are relatively long are not constrained. The free vibration of the FG-CNTRC conical shell panel was investigated with respect to the shell radius 0 R . It is observed from Figure 9a,b that 1 ω shows a uniform decrease with increasing of the value of 0 R , regardless of the CNT distribution pattern and the panel length. It is because the pane stiffness decreases while its total mass increases as the shell radius increases. However, the decreased slope becomes lower in proportion to the shell radius. Regarding the CNT distribution pattern, the magnitude order of 1 ω is the same as in the previous simulation cases. Meanwhile, it is seen from  The free vibration of the FG-CNTRC conical shell panel was investigated with respect to the shell radius R 0 . It is observed from Figure 9a,b thatω 1 shows a uniform decrease with increasing of the value of R 0 , regardless of the CNT distribution pattern and the panel length. It is because the pane stiffness decreases while its total mass increases as the shell radius increases. However, the decreased slope becomes lower in proportion to the shell radius. Regarding the CNT distribution pattern, the magnitude order ofω 1 is the same as in the previous simulation cases. Meanwhile, it is seen from Figure 9b thatω 1 uniformly increases in proportion to , but which is entirely owing to the calibration with 2 . The non-calibrated fundamental frequency ω 1 decreases with increasing the value of because the mass increase is superior to the stiffness increase as the shell becomes longer.
proportion to the shell radius. Regarding the CNT distribution pattern, the magnitude order of 1 ω is the same as in the previous simulation cases. Meanwhile, it is seen from Figure 9b that 1 ω uniformly increases in proportion to  , but which is entirely owing to the calibration with 2  . The non-calibrated fundamental frequency 1 ω decreases with increasing the value of  because the mass increase is superior to the stiffness increase as the shell becomes longer.
(a) (b) Next, two types of three-layered symmetric sandwich FG-CNTRC conical shell panels are considered. One is composed of a PMMA homogeneous core and two FG-CNTRC faces, and the other consists of a FG-CNTRC core and two PMMA homoge- Next, two types of three-layered symmetric sandwich FG-CNTRC conical shell panels are considered. One is composed of a PMMA homogeneous core and two FG-CNTRC faces, and the other consists of a FG-CNTRC core and two PMMA homogeneous faces. The free vibration response was investigated with respect to the core thickness ratio h c /h. It is seen, from Figure 10a for the PMMA homogeneous core, thatω 1 slightly increases with the increase of h c /h, but thereafter, it drops remarkably as h c /h tends to unity. It is because the mass decrease is superior for the core increase up to a certain thickness, but a further increase in core thickness gives rise to a remarkable decrease in panel stiffness. This trend is shown to be severer as the semi-vertex angle α increases. On the other hand, the CNTRC core shows a different variation ofω 1 to the core thickness, as represented in Figure 10b. The normalized first frequencyω 1 uniformly increases in proportion to h c /h, and the slope of the increase becomes higher with the increase of h c /h and α. It is because the increase of CNTRC core gives rise to the increase of panel stiffness, and this effect becomes more significant as h c /h and α increase.

Conclusions
The free vibration of functionally graded CNT-reinforced composite (FG-CNTRC) conical shell panels was analyzed by a 2-D numerical method. The numerical method was developed in the framework of the 2-D planar natural element method (NEM) by introducing a geometry transformation between the shell surface and the planar NEM grid and by adopting the MITC3+shell element. The benchmark experiments were per-

Conclusions
The free vibration of functionally graded CNT-reinforced composite (FG-CNTRC) conical shell panels was analyzed by a 2-D numerical method. The numerical method was developed in the framework of the 2-D planar natural element method (NEM) by introducing a geometry transformation between the shell surface and the planar NEM grid and by adopting the MITC3+shell element. The benchmark experiments were performed to validate the developed 2-D NEM-based numerical method. Moreover, the free vibration responses of FG-CNTRC conical shell panels were profoundly investigated with respect to the associated parameters. The numerical results draw the following major observations:

•
The numerical method reliably and effectively solves the free vibration of FG-CNTRC conical shell panels without causing shear locking with the maximum relative difference of less than 5% to the reference solutions, even for coarse and 2-D planar NEM grids.

•
The normalized first frequencyω 1 increases proportional to the CNT volume fraction V * CNT , but it uniformly decreases with increasing the values of semi-vertex angle α and shell radius R 0 . • To the sub-tended angle θ 0 , theω 1 does not show remarkable change except for slight fluctuation until the decrease of θ 0 to a certain value. However, a further decrease of θ 0 causes the sudden increase ofω 1 , except for CFCF. The increased intensity is dependent on the boundary condition, but it is not sensitive to R 0 /h. • Theω 1 uniformly increases with thickness h and length owing to the calibration, but actually, its non-calibrated value ω 1 uniformly decreases proportional to these two parameters.

•
In the symmetric sandwich FG-CNTRC conical shell panels,ω 1 for the PMMA core slightly increases and then drops remarkably proportional to h c /h, and vice versa for the CNTRC core.