Axial and Shear Buckling Analysis of Multiscale FGM Carbon Nanotube Plates Using the MTSDT Model: A Numerical Approach

The present paper investigates the axial and shear buckling analysis of a carbon nanotube (CNT)-reinforced multiscale functionally graded material (FGM) plate. Modified third-order deformation theory (MTSDT) with transverse displacement variation is used. CNT materials are assumed to be uniformly distributed, and ceramic fibers are graded according to a power-law distribution of the volume fraction of the constituents. The effective material properties are obtained using the Halpin–Tsai equation and Voigt rule of the mixture approach. A MATLAB code is developed using nine noded iso-parametric elements containing 13 nodal unknowns at each node. The shear correction factor is eliminated in the present model, and top and bottom transverse shear stresses are imposed null to derive higher-order unknowns. Comparisons of the present results with those available in the literature confirm the accuracy of the existing model. The effects of material components, plate sizes, loading types, and boundary conditions on the critical buckling load are investigated. For the first time, the critical buckling loads of CNT-reinforced multiscale FGM rectangular plates with diverse boundary conditions are given, and they can be used as future references.


Introduction
In the analysis and design of all civil engineering structures, the buckling response of the CNT-reinforced FGM plate caught the attention of many researchers in recent years. Currently, critical buckling loads are obtained using the Corr and Jennings [1] simultaneous iteration technique. The critical buckling load is the maximum load in the elastic range of the material above which plates start to deflect laterally. If the material is stressed beyond the elastic range and into the non-linear (plastic) range, the buckling strength of a plate is smaller than the elastic buckling strength of a plate. When the load approaches the critical buckling load, the plate will bend significantly, and the material's stress-strain behavior will diverge from linear. In FGM-type composite material, properties of material constituents are varied according to the required performance. In this paper, the material constituents were a metal matrix, CNT reinforcement, and fiber. The final material was made in two phases. Here, we calculated the minimum edge compressive load in the form of the non-dimensional critical buckling load, which is required to initiate the instability of the plate structure. FGM is widely employed in many areas such as machine, construction, defense, electronic, chemical, pharma, energy sources, nuclear, automotive, and shipbuilding industries. Because of the expanding use of FGMs in a range of structural applications, detailed theoretical models are required to anticipate their behavior. the buckling behavior. A refined plate theory based on the secant function was used by Abdulrazzaq et al. [30] to study the thermal buckling stability of clamped nano-size FGM plates. From their study, it can be observed that the buckling behavior of clamped FGM nanoplates was very sensitive to various parameters such as aspect and side-to-thickness ratios, material graduation, thermal condition, etc. The study of the influence of small-scale parameters on the vibration and buckling behavior of CNT-reinforced FGM plates was done by Shahraki et al. [31]. The CNT-based FGM nanoplate was considered to rest on a Kerr elastic foundation. Costa and Loja [32] represented the static analysis of a dual-phase moderately thick FGM plate. The CNT reinforcements were assumed to be added to the matrix material in the first phase.
Even though various studies on the buckling of FGM plates have been conducted based on a range of plate theories, no studies on the buckling of multiscale FGM plates based on the MTSDT theory were found. The present MTSDT mathematical theory has been modified to represent the kinematics field that captures normal and transverse crosssection deformation modes. The assumed in-plane fields incorporate the cubic degree of thickness terms and quadratic degree of thickness terms for the transverse component. The C1 continuity requirement associated with third-order shear deformation theory is avoided by developing a C0 FE formulation by replacing the out-of-plane derivatives with independent field variables. The present study can be used for the design and analysis of various types of hybrid composite curve panels, which are used in various engineering fields. The design charts can be obtained by the present model, which may be useful for the designer. Material properties, such as Young's modulus, are supposed to change with plate thickness according to a power-law distribution of the volume percentage of the constituents. To the best of the authors' knowledge, no experimental results on the present work are available in the literature; hence, present model results were validated with the closed-form elasticity solution and numerical analysis results available in the literature. To study the influence of various parameters, the non-dimensional critical buckling load was calculated for numerical analysis.

Geometrical Configuration and Effective Material Properties
A multiscale FGM plate of length a, width b, and thickness h, as shown in Figure 1 was considered. In the buckling response of the plate, the rectangular Cartesian platform coordinates × and y were used. The co-ordinate planes × = 0, a and y = 0, b define the boundaries of the plate. The reference surface is the middle surface of the plate, defined by z = 0, where z is the thickness co-ordinate measured from the un-deformed middle surface of the plate. The performance of these FGM plates might be improved by using a multiscale hierarchical FGM as shown in Figure 2, which is made possible by combining the continuous fiber phase, the metal matrix, and CNT reinforcement. In such circumstances, the overall homogenization process can be divided into two phases: in the first phase, the dispersion of CNT in the metallic matrix yields a nanocomposite, and in the second phase, this nanocomposite receives ceramic inclusions in a graded manner, resulting in a CNT-reinforced multiscale composite. Since the CNTs are expected to be evenly distributed and randomly oriented throughout the matrix, the final mixture is considered an isotropic mixture. It is also expected that the bonding between CNT and matrix and dispersion of CNT in the matrix are perfect. Each CNT is assumed to be straight and has the same aspect ratio and mechanical properties. The matrix material is considered void-free, and the bonding between the matrix and fiber is excellent. To evaluate the effective elastic properties of the material, a suitable approach should be adopted. A combination of the Halpin-Tsai equation [33] and homogenization scheme can be adapted to predict the effective material properties of a three-phase multiscale FGM plate. The Halpin-Tsai equation is an empirical formula, known to be fit for calculating effective material properties of the mixture of the matrix and low fraction of the CNT reinforcement. The elastic properties of an anisotropic mixture of CNT and the matrix can be expressed as follows: The volume fraction of carbon nanotube V CN and Poisson's ratio of the nanocomposite ν MNC are calculated as [34].
The volume fraction of dispersed fiber constituents is expressed as follows: where h and Z are the respective total thickness and thickness coordinate in the transverse direction, having an origin on the middle surface of the plate. The exponential power n permits the ceramic fiber to fluctuate in the thickness direction. The effective material characteristics of the final material fluctuate continuously according to Equation (5). In this paper, effective elastic material properties are calculated using a homogenization approach based on the Voigt rule of the mixture. as shown below: Because of the dispersion of carbon nanotubes in the metal matrix, the effective Young's modulus of the nanocomposite phase may be used instead of Young's modulus of the matrix phase in the preceding equation. In this work, we assume the dispersion of carbon nanotubes in metal; therefore, we must first compute the effective material properties of the nanocomposite.

Governing Equation
The governing equation for buckling analysis is derived by using the MTSDT mathematical model. A rectangular plate of size (a × b) is assumed to be perfect in geometry.

Displacement Equation
The in-plane displacement (u and v) and transverse displacement (w), which is based on the MTSDT, are represented as follows: In the above matrix [A], all higher-order terms are determined by eliminating the transverse shear (τ xz = τ yz = 0) at the outer surface of the plate; then, the modified in-plane displacement field is as follows: where f 1 (z) = z − f 2 (z) and f 2 (z) = 4z 3 3h 2 . The final expression for the in-plane displacement and transverse displacement fields: To replace the C 1 continuity with C 0 continuity to assure the field variables are continuous within the element, the out-of-plane derivatives are replaced by the following relation: However, due to the above substitution, there is an introduction of additional nodal unknowns that impose extra constraints, which are enforced variationally through a penalty approach as follows: Hence, in the present formulation, the displacement variables are as follows:

Strain Displacement Relationship
The linear strain corresponding to the displacement fields is expressed as follows: Further incorporation of the final expression for the displacement field (Equation (9)) into the above equation leads to the following expression: The above strain equation can be generalized into the following expression: where {ε} 6X1 = ε xx ε yy ε zz γ xy γ xz γ yz and The relationship between the strain vector {ε} and displacement vector {d} can express by the following relationship:

Element Description and Shape Function
A nine-noded iso-parametric element (shown in Figure 3) was employed for the present finite element model with 13 unknown variables at each node. The nodal unknowns at any point within the domain were expressed in terms of the shape function. At each element, the displacement field and the element geometry are defined as follows: The shape function N i is the function of the natural coordinate system used in the finite element modeling, and it is expressed as follows:

Constitutive Relationship
In this study, we considered that the multiscale composite material is an isotropic material at each point of its domain, and the constitutive relationship between stress and strain is as follows: (19) where the constitutive matrix is expressed as [35];

Buckling Analysis
The strain energy of the plate may be written as By putting the value of Equation (19) in the above Equation (21), we obtain The global stiffness matrix of the multiscale composite plate is obtained by equating the total energy of the system to zero. (23) To derive the membrane stiffness matrix [K m ], the membrane strain associated with the deflection can be calculated as [36] follows:

[K] = [B] T [D][B]dxdy
Or it can be written as The matrix {θ} and strain vector {ε m } can be related as 3h 2 , and f 6 (z) = − z 3 3 . By using Equation (26) and the strain displacement relationship, the stress stiffness matrix [K m ] can be written as follows: dz and the stress matrix [S] in terms of plane stress N x , N y , and N xy can be expressed as follows: The governing equation for calculating the critical buckling is expressed as follows:

Computation of the Critical Buckling Load
In this analysis, the governing equation for buckling analysis [K]{δ} = λ[K m ]{δ} was solved by the simultaneous iteration technique of Corr and Jennings [1] for the computation of eigenvalues and eigenvectors. In this method, [K] is positive definite and can be decomposed into Cholesky factors as The governing equation for buckling analysis characterizes the standard eigenvalue problem, and these have been solved to extract the eigenvalues and the eigenvectors. In this equation, 1/λ is the eigenvalue. Therefore, the eigenvalue corresponding to the lowest buckling loads is obtained using the simultaneous iteration technique. The methodology is explained as follows: 1.
Set a trial vector [U] and ortho-normalize.

2.
Back The numerical results are calculated in the form of non-dimensional critical buckling as shown below: The rectangular plates shown in Figure 4 are subjected to in-plane loading in two different directions. In the given Figure 4

Numerical Results
To calculate the critical buckling load, the eigenvalue problem was determined. Various numerical results were used to obtain the mechanical buckling behavior of CNTreinforced multiscale FGM rectangular plates using the proposed 9-noded isoparametric elements. The finite element code was developed in Matlab to perform the numerical simulation. The numerical values were calculated for 3 × 3 gauss integration points. The material components adopted in this study are listed in Table 1. Table 1. Material properties of constituents.

Constituents' Material Material Properties
Matrix

Comparison and Convergence
To determine the best suitable mesh size for the present numerical analysis, the given plate was divided into various mesh sizes in the ×and y-directions. This convergence study was carried out for different volume fractions of ceramic fiber and with a side-tothickness ratio a/h = 10, as shown in Table 2. The non-dimensional critical buckling load was determined for mesh sizes varying from 2 × 2 to 6 × 6. It was observed that the critical buckling load converged for the mesh size 5 × 5. Therefore, a 5 × 5 mesh size was adopted for the complete numerical analysis.
To validate the present MTSDT theory, the non-dimensional critical buckling load was calculated for a different side-to-thickness ratio of simply supported square plates under uniaxial and biaxial compressive loadings. The numerical values in Table 3 represent the critical buckling load for the Al/Al 2 O 3 plate with 0% weight fraction of CNT reinforcement. The presented numerical results were compared with a previous numerical study [37] and are in good agreement with the reference. The mode shape of a simply supported plate for the first three nodes is presented in Figures 5 and 6 for n = 0 and n = 1, respectively.

Effect of Boundary Conditions on Uniaxial and Biaxial Compression
The variation of the non-dimensional critical buckling load for various boundary conditions is represented in Table 4. The numerical values were calculated for 0%, 2.5%, and 5% weight fraction of CNT reinforcement under uniaxial and biaxial loading. From Table 4, the maximum value of the critical buckling load was obtained by clamped (CCCC) boundary conditions, whereas the CFCF boundary condition yielded the minimum value of the critical buckling load. The CCCC boundary condition indicated that the plates were fixed on all four sides, and the CFCF boundary condition indicated that the plates were fixed and free on adjacent sides. In the case of a 0% weight fraction of CNT, approximately (80-85)% difference in the critical buckling load was observed between the CCCC and CFCF boundary conditions, and a (60-63)% difference was observed between the CCCC and SSSS boundary conditions. However, if we assume a 5% weight fraction of CNT in the mixture, slightly less difference was observed between these boundary conditions. From Table 4, it was also observed that for all boundary conditions, the plate had a higher critical buckling load under uniaxial compression than under biaxial compression. The first three mode shapes of the plate under biaxial compressive and shear loading are presented in Figures 5 and 6, respectively. The mode shapes were drawn for simply supported and clamped-free boundary conditions. As seen from the mode shape of the plate, the essential boundary conditions were satisfied at the supports.

Effect of CNT and Volume Fraction Index (n) on the Critical Buckling Load
The numerical results for the critical buckling load at a different weight fraction of SWCNT and MWCNT are presented in Tables 5 and 6 under uniaxial and biaxial compressive loading, respectively. All numerical results were obtained for the 1st six modes, and plates were restrained with a simple supported condition. In this case, the Al/ZrO2 plate was assumed to be reinforced with SWCNT and MWCNT at 0%, 2.5%, and 5% weight fractions. From these tables, it was observed that SWCNT performed better than MWCNT under uniaxial and biaxial compression. At a 5% weight fraction of the CNT, SWCNT had a 17% higher critical buckling load for n = 0.5 and 37% higher critical buckling load for n = 10 than MWCNT. This happened because of the magnitude of the elastic properties of the SWCNT and MWCNT. Since at n = 0 only ZrO 2 fibers were present, no difference was observed. From n = 0.5 to 10, the proportion of ZrO 2 started to decrease and the proportion of the nanocomposite started to increase. Due to this increase in the nanocomposite proportion, a greater difference was observed at n = 10. By increasing the volume fraction index from n = 0 to n = 10, the amount of fiber in the mixture decreased, which led to a decrease in the stiffness of the plate. Therefore, the critical buckling load decreased as the volume fraction index increased. By increasing the weight fraction of the CNT up to 5%, the critical buckling load increased by 43% in the SWCNT case and 13% in the MWCNT case because of the stiffness of the plate increased by increasing the amount of CNT in the mixture. Figure 7 shows plots for freely supported and clamped boundary conditions for different SWCNT fractions and fiber volumes. At W_cnt = 5%, the plate had an approximately equal critical buckling load from n = 0 to n = 10. The plate with a multiscale phase behaved similar to the plate with only ceramic fiber at W_cnt = 5%. At W_cnt = 0% and 2.5%, the critical buckling load decreased with an increase in the volume fraction index (n), but a greater decline was observed at the 0% weight fraction of the CNT.

Effect of the Side-to-Thickness Ratio (a/h) and Aspect Ratio (b/a) of the Plates
The effect of the side-to-thickness ratio of the plate is presented in Figures 8 and 9. Figures are plotted for simply supported and clamped boundary conditions. In this case, all values were calculated for different volume fraction index (n) and 0% weight fraction values of the CNT. From Figures 8 and 9, it can be observed that under ceramicrich conditions, i.e., n = 0, the plate had the maximum critical buckling load, and in the case of nanocomposite-rich conditions, i.e., n = 10, the plate had the minimum critical buckling load. Figures 8 and 9 present the uniaxial compression and biaxial compression, respectively. Under uniaxial and biaxial compression, the critical buckling load for simply supported plates increased by increasing the a/h ratio up to 20; after that, no significant change in the critical buckling load was observed. This is because in a simply supported plate, the stiffness of the plate increases to only a/h = 20, and the same variation was observed by Reddy et al. [37]. In the case of a clamped supported plate, the critical buckling load increased to a/h = 100 because under clamped support conditions, the stiffness of the plate increased to a/h = 100.  The variation of the critical buckling load with the aspect ratio of the plate can be seen in Figures 10 and 11. For the given simply supported and clamped boundary conditions, numerical results were obtained for a constant side-to-thickness ratio of the plate, i.e., a/h = 10. In Figures 10 and 11, it is observed that the critical buckling load increases by increasing the b/a ratio of the plate for both types of loading conditions.

Critical Buckling Load for Various Types of Plates
The variation in the critical buckling load for various types of plates by considering the different volume fraction indexes is presented in Tables 7 and 8. All numerical values were calculated for the critical buckling load for the 1st six modes for simply supported plates. A plate made of a different type of metal and ceramic fibers behaves differently under uniaxial and biaxial compression. Here, we assumed that the 0% weight fraction of the CNT was used as a reinforcement. In the case of the Al/Al 2 O 3 plate, the critical buckling load increased by (60-74)% by increasing the volume fraction ratio under uniaxial and biaxial compression. In the case of the Ti-6Al-4V/ZrO 2 plate, the value for the critical buckling load increased by only (23-28)% by increasing the volume fraction index. Under uniaxial and biaxial compression, the Al/Al 2 O 3 plate had the highest critical buckling load value among all types of plates made of different metal matrix and fiber components. The differences in the critical buckling load values are due to the different elastic modulus values of the components. In the case of the Al/Al 2 O 3 plate, the difference in the elastic modulus of the Al matrix and Al 2 O 3 fiber was much larger. Due to this fact, greater variation in the critical buckling load was observed by increasing the volume fraction index.

Effect of Biaxial and Shear Loading of the Plate
The non-dimensional critical buckling load for a simply supported plate under various in-plane forces is presented in Tables 9 and 10. Numerical results were calculated for different fiber volume fractions and weight fractions of CNT reinforcement. Table 9 represents the variation in the critical buckling load for the 1st mode under various shear loading and constant biaxial loading values (N y /N x = 1). It is noted that by increasing the shear loading from 0 to 2, the non-dimensional critical buckling load decreased by 18% for all fractions of the CNT reinforcement. The reason for this is increased shear loading, reducing the stiffness of the plate.   Table 10 represents the variation in the critical buckling load for various in-plane compressive and shear forces. All values were calculated for a simply supported plate at different weight fractions of the CNT. From Table 10, it can be seen that by increasing the ratio of in-plane compressive forces in the y and × directions, the critical buckling load for all shear loads is reduced. Further, for all uniaxial and biaxial compressive forces, the critical buckling load decreases as the shear load increases. This is because as the compressive and shear loads increase, the buckling resistance of the plate decreases.

Conclusions
In this paper, an MTSDT mathematical theory was adopted to represent the kinematic field. The in-plane displacement fields integrate the cubic degree of thickness terms and quadratic degree of thickness terms for the out-of-plane displacement field. A nine-noded isoparametric element with 13 unknowns at each node was adopted for the finite element formulation. Effective elastic properties of the multiscale FGM material were predicted by using the Halpin-Tsai equation and the Voigt rule of mixture approach. The effect of various parameters on the critical buckling behavior of a multiscale FGM plate is presented, and the following conclusions were drawn from this numerical analysis:

•
The critical buckling load parameter was at a maximum under clamped boundary conditions. • By increasing the volume fraction index (n), the critical buckling is decreased due to less stiffness being obtained at a higher volume fraction index.

•
As the weight fraction of CNT increased, the critical buckling load increased because CNT imparted more stiffness to the material.

•
The side-to-thickness ratio (a/h) and aspect ratio (b/a) of the plates had a significant impact on the buckling behavior of the plate. Increasing the a/h ratio increased the critical buckling load, and increasing the b/a ratio decreased the critical buckling load. • Due to the given elastic properties of the Al and Al 2 O 3 , the Al/Al 2 O 3 plate yielded the maximum value of the critical buckling load among all plates.

•
For the same ratio of in-plane compression in the y-and x-direction, the critical buckling load decreased with increases in in-plane shear loading.

•
For all values of in-plane shear loading, the critical buckling load decreased with an increase in the ratio of in-plane compression in the y-and x-direction.
It was observed in the present study that CNT fibers and reinforcement play a very important role in the buckling response of a plate structure. The results presented in this study are new for the buckling behavior of multiscale FGM plates. Therefore, it is believed that the results obtained are very useful for the analysis and design of this type of plate structure.   In-plane and out-of-plane displacement on the midplane. ω 1 and ω 2 Rotation of the normal about the midplane on y-and x-axes.