Isogeometric Analysis for Free Vibration of Functionally Graded Plates Using a New Quasi-3D Spectral Displacement Formulation

: This paper presents a novel quasi-three-dimensional shear deformation theory called the spectral displacement formulation (SDF) for analyzing the free vibration of functionally graded plates. The SDF expresses the unknown displacement ﬁeld as a unique form of the Chebyshev series in the thickness direction. By increasing the truncation number in the Chebyshev series, the bending analysis results can approach the three-dimensional elasticity solution and satisfy the traction-free boundary conditions without requiring a shear correction factor. The SDF is an extension of the classical plate theory, thereby naturally avoiding the shear-locking phenomenon. These characteristics enable the SDF to apply to plates of arbitrary thickness while maintaining accuracy. The nonuniform rational B-spline-based isogeometric approach is employed to enhance the applicability of this theory to free vibration analysis of functionally graded plates with complex geometries and different boundary conditions. Numerical examples are presented to demonstrate the accuracy and reliability of the proposed method in analyzing the free vibration of functionally graded plates.


Introduction
Functionally graded materials (FGMs) are multiphase composite materials characterized by the variation in composition and structure over volume, resulting in their thermodynamic and physical properties also exhibiting gradient changes across the volume.This unique characteristic enables FGMs to possess tailored and optimized properties to meet specific engineering requirements.Due to their exceptional mechanical strength, thermal shock resistance, and high-temperature durability, among other qualities, FGMs are widely applicable in diverse fields, including aerospace equipment, electronic devices, automotive engines, chemical components, and medical equipment [1,2].
The mechanical analysis of functionally graded material (FGM) plates has attracted significant attention from researchers [3][4][5].The classical plate theory (CPT) [4], based on the Kirchhoff-Love assumptions, is commonly employed to analyze thin FGM plates.Nevertheless, its applicability to moderately thick FGM plates may produce less accurate results.The first-order shear deformation theory (FSDT) [4], which takes the shear effect into account, is therefore developed to improve the CPT.However, a shear correction factor (SCF) is necessary to modify the calculation of shear stresses, and determining the appropriate value of SCF can be challenging as it varies between different problems.Additionally, FSDT is susceptible to the shear-locking problem when the thickness to length ratio becomes very small.Numerous types of higher-order shear deformation theories (HSDTs) [4] have been developed to elucidate the nonlinear, parabolic changes in transverse shear stresses throughout the thickness, thereby obviating the need to calculate the SCF.HSDTs can be categorized into models based on polynomial shape functions, such as the third-order shear deformation theory (TSDT) [6], the fifth-order shear deformation theory [7], and the nthorder shear deformation theory [8].Alternatively, there are models based on nonpolynomial shape functions, such as trigonometric HSDTs [9][10][11][12], exponential HSDTs [13,14], hyperbolic HSDTs [15][16][17][18], and hybrid HSDTs [19][20][21].Traction-free boundary conditions on the top and bottom surfaces are automatically satisfied in these HSDTs, and the displacement and stress results obtained from HSDTs are also more accurate than those obtained from CPT and FSDT.
Carrera et al. [22] demonstrated the significance of the thickness stretching effect on the response analysis of FGM plates and shells.Through multiple comparative studies, they found that considering the thickness stretching effect is crucial in the computation of moderately thick FGM plates and shells.Furthermore, they discovered that in the case of sandwich structures, including thin plates and shells, the thickness stretching effect must be taken into account.Quasi-three-dimensional (3D) theories are HSDTs with high-order variations throughout the thickness for the transverse displacement [4].By accounting for higher-order variations in both in-plane and transverse displacements, the quasi-3D theory takes into account the shear deformation effect and the thickness stretching effect.Generally, Carrera's unified formulation (CUF) proposed by Carrera [23] and extended by Demasi [24] is commonly used for implementing quasi-3D theories.Numerous quasi-3D theories have been proposed in the literature.Zenkour [25] developed a quasi-3D theory for FGM plates using trigonometric functions.Matsunaga [26,27] constructed a quasi-3D theory by expanding in-plane and transverse displacements using power series.Ferreira et al. [28] employed CUF to develop a sinusoidal quasi-3D shear deformation theory.Neves and Ferreira et al. [29,30] conducted similar work using hyperbolic or hybrid functions instead of sinusoidal functions.Mantari and Soares et al. [19] introduced a general formulation that allows for deriving various quasi-3D theories using polynomials, trigonometric, or hybrid functions.They [31,32] optimized the number of unknowns in tangential-based quasi-3D shear deformation theories by separating transverse displacement into bending and shear compohnents, requiring only four unknowns to describe the unknown displacement field model.Thai and Kim et al. [33] developed a quasi-3D theory with five unknowns using sinusoidal functions, where the transverse displacement is divided into bending, shear, and stretching parts.Following a similar approach, Thai et al. [16], Hebali et al. [17], Bessaim et al. [34], and Bennoun et al. [35] employed hyperbolic functions to construct quasi-3D theories.At the same time, Belabed et al. [36] and Mantari et al. [20,21] used combined exponential and hyperbolic functions and combined exponential and trigonometric functions to develop quasi-3D theories.Quasi-3D theories, as modifications of HSDTs, typically yield more accurate results in response analysis compared to conventional HSDTs.
The free vibration analysis of FGM plates based on quasi-3D shear deformation theories has been extensively discussed.Talha and Singh [37] designed a nine-node C 0 continuous isoparametric element with 13 degrees of freedom per node based on quasi-3D theory.They investigated the influence of aspect ratio, thickness ratio, volume fraction index, and boundary conditions on the bending and free vibration responses of FGM plates.Matsunaga [26,27] considered the effects of transverse shear and normal stresses, as well as rotational inertia on the natural frequencies and buckling stress of FGM plates and shallow shells.Thai and Choi [38] studied the free vibration of FGM plates on a Pasternak foundation using a new quasi-3D shear deformation theory with only four unknowns.Hebali et al. [17] analyzed the bending and free vibration of FGM plates using a novel hyperbolic quasi-3D theory.Abualnour et al. [39] developed a new trigonometric quasi-3D theory to analyze the vibration problems of FGM plates, where the displacement field includes undetermined integral terms.Mantari et al. [20,21] investigated the free vibration analysis of FGM plates resting on elastic foundations by using a generalized quasi-3D hybrid-type higher-order shear deformation theory (HSDT).The significant feature of this formulation is that it considers only five unknowns instead of six unknowns in the well-known trigonometric plate theory (TPT).Zenkour and Sobhy [40] utilized quasi-3D theory to analyze the statics, buckling, and free vibration of simply supported sandwich plates.Li et al. [41] employed a new quasi-3D theory to analyze the free vibration of simply supported FGM plates.This new quasi-3D theory introduces two reformulated transverse shear strain functions to describe the distribution of transverse shear strains along the plate thickness.Additionally, using an undetermined integral format of displacement field enhances computational efficiency.Shahsavari et al. [42] presented a novel quasi-3D hyperbolic theory for the free vibration analysis of FGM porous plates resting on elastic foundations.A comprehensive parametric study is carried out to consider the effects of porosity fraction index, stiffness of foundation parameters, and volume fraction index on the frequencies of imperfect FGM plates.Hai Van et al. [43] investigated the free vibration of nonuniform-thickness 2D-FGM sandwich porous plates using a quasi-3D theory.The findings revealed that the mechanical behaviors of the porous sandwich plates are significantly affected by the hygrothermal environment, the boundary conditions, and the patterns of porosity.Jafari and Kiani [44] obtained a free vibration response of composite plates which are reinforced with graphene platelets using a quasi-3D plate model.Pham et al. [45] proposed the free vibration response from bi-directional FGM rectangular plates in the fluid medium based on a quasi-3D plate theory.Kaddari et al. [46] investigated a new quasi-3D hyperbolic shear deformation theory to discuss the statics and free vibration of FGM porous plates resting on elastic foundations.The influences of the porosity parameter, power-law index, aspect ratio, thickness ratio, and the foundation parameters on the vibration of porous FGM plate are considered.Daikh and Zenkour [47] studied the free vibration of porous FGM sandwich plates using a new and simple higher-order shear deformation theory.The effect of porosity, sandwich plate geometry, and heterogeneity parameters on the FGM sandwich plate's free vibration was investigated.
As previously stated, considerable focus has been on analyzing the free vibration of FGM plates.It is worth noting that in numerical analysis, HSDTs generally require C 1 continuity, which can be easily achieved through the isogeometric method.Isogeometric analysis (IGA), introduced by Hughes et al. [48], aims to seamlessly integrate finite element analysis (FEA) with computer-aided design (CAD).By directly employing the widely used nonuniform rational B-splines (NURBS) basis functions from CAD software in FEA, IGA methods exhibit notable advantages, primarily due to their ability to accurately represent geometric shapes.NURBS basis functions can precisely describe the geometry of the research object, even with relatively sparse mesh discretization, thereby ensuring high numerical accuracy.Moreover, NURBS-based IGA offers flexibility in satisfying arbitrary order continuity requirements, which proves challenging in conventional FEA.Consequently, IGA methods have found extensive application in analyzing FGM plates.
This paper presents a novel quasi-3D theory called the spectral deformation formulation (SDF) [49,50], where the unknown displacement field is expanded as a unique form using the Chebyshev series.Similar to CUF, this theory exhibits hierarchical refinement capability, enabling a perfect match with the 3D elasticity solution.Additionally, the SDF is an extension of CPT, thus naturally avoiding the shear-locking problem.Consequently, this theory is applicable to both thin and thick plates while maintaining high accuracy [50].This study analyzes the free vibration of FGM plates by incorporating this theory with IGA methods.Numerical results are compared with existing literature, demonstrating the reliability and accuracy of this quasi-3D theory in the analysis of free vibration in FGM plates.
The paper is organized as follows.Section 2 describes the FGM plate model and the quasi-3D deformation theory based on the SDF.The governing equations are derived using the d'Alembert principle and the principle of virtual work, followed by discretization using NURBS.Section 3 presents relevant numerical examples and discussions.Section 4 provides the corresponding conclusions.

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coordinates x and y represent the in-plane directions, and z denotes the thickness direction (see Figure 1).

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coordinates  and  represent the in-plane directions, and  denotes the thickness direction (see Figure 1).The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic material, while the bottom surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a power-law distribution [3].The relationship between the volume fractions of the ceramic and metal is defined as  m +  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the metal, respectively, and  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of the FGM plate can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's ratio, and mass density of the ceramic and the metal, respectively., , and  are the equivalent elastic modulus, Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence, allowing for high-precision approximation of target functions with few truncation terms [52].In this study, the first kind of Chebyshev polynomials is employed as the spectral basis to expand the displacement field of the FGM plate in the thickness direction as the following spectral series: where  = [, , ] T and  ,  , and  are displacements in the  ,  , and  directions, respectively. and  are the truncation numbers for the transverse and the in-plane displacement fields, respectively, with the condition that  ≤ .{  } =0  are the generalized displacements defined as follows: The top surface of the plate (z = h/2) consists of ceramic material, while the bottom surface (z = −h/2) transitions from ceramic to metal using a power-law distribution [3].The relationship between the volume fractions of the ceramic and metal is defined as

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coordinates  and  represent the in-plane directions, and  denotes the thickness direction (see Figure 1).The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic material, while the bottom surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a power-law distribution [3].The relationship between the volume fractions of the ceramic and metal is defined as  m +  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the metal, respectively, and  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of the FGM plate can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's ratio, and mass density of the ceramic and the metal, respectively., , and  are the equivalent elastic modulus, Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence, allowing for high-precision approximation of target functions with few truncation terms [52].In this study, the first kind of Chebyshev polynomials is employed as the spectral basis to expand the displacement field of the FGM plate in the thickness direction as the following spectral series: where  = [, , ] T and  ,  , and  are displacements in the  ,  , and  directions, respectively. and  are the truncation numbers for the transverse and the in-plane displacement fields, respectively, with the condition that  ≤ .{  } =0  are the generalized displacements defined as follows:

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coordinates  an  represent the in-plane directions, and  denotes the thickness direction (see Figure 1).The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic material, while the bottom surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a power-law distribution [3 The relationship between the volume fractions of the ceramic and metal is defined as  m  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the metal, respectively, an  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of the FGM pla can be computed using the following expression:

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence, allowing fo high-precision approximation of target functions with few truncation terms [52].In this stud the first kind of Chebyshev polynomials is employed as the spectral basis to expand the di placement field of the FGM plate in the thickness direction as the following spectral series: where  = [, , ] T and  ,  , and  are displacements in the  ,  , and  directions, r spectively. and  are the truncation numbers for the transverse and the in-plane di placement fields, respectively, with the condition that  ≤ .{  } =0  are the generalize displacements defined as follows: Mathematics 2023, 11, x FOR PEER REVIEW

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coor  represent the in-plane directions, and  denotes the thickness direction (s The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic material, wh surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a power-law di The relationship between the volume fractions of the ceramic and metal is de  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the metal, res  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's ratio, and of the ceramic and the metal, respectively., , and  are the equivalent ela Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence high-precision approximation of target functions with few truncation terms [52 the first kind of Chebyshev polynomials is employed as the spectral basis to e placement field of the FGM plate in the thickness direction as the following spe where  = [, , ] T and  ,  , and  are displacements in the  ,  , and  spectively. and  are the truncation numbers for the transverse and the placement fields, respectively, with the condition that  ≤ .{  } =0  are th displacements defined as follows: is expressed as: Mathematics 2023, 11, x FOR PEER REVIEW

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where th  represent the in-plane directions, and  denotes the thickness direc ⁄ ) transitions from ceramic to metal using a power-The relationship between the volume fractions of the ceramic and met  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the me  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material proper can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's rat of the ceramic and the metal, respectively., , and  are the equival Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral conve high-precision approximation of target functions with few truncation ter the first kind of Chebyshev polynomials is employed as the spectral ba placement field of the FGM plate in the thickness direction as the follow where  = [, , ] T and  ,  , and  are displacements in the  ,  , spectively. and  are the truncation numbers for the transverse a placement fields, respectively, with the condition that  ≤ .{  } =0  displacements defined as follows: where Mathematics 2023, 11, x FOR PEER REVIEW 4 of

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coordinates  a  represent the in-plane directions, and  denotes the thickness direction (see Figure 1) The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic material, while the botto surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a power-law distribution [ The relationship between the volume fractions of the ceramic and metal is defined as  m  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the metal, respectively, a  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of the FGM pla can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's ratio, and mass dens of the ceramic and the metal, respectively., , and  are the equivalent elastic modulu Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence, allowing f high-precision approximation of target functions with few truncation terms [52].In this stud the first kind of Chebyshev polynomials is employed as the spectral basis to expand the d placement field of the FGM plate in the thickness direction as the following spectral series: where  = [, , ] T and  ,  , and  are displacements in the  ,  , and  directions, spectively. and  are the truncation numbers for the transverse and the in-plane d placement fields, respectively, with the condition that  ≤ .{  } =0  are the generaliz displacements defined as follows:

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, where the coordinates  and  represent the in-plane directions, and  denotes the thickness direction (see Figure 1).The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic material, while the bottom surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a power-law distribution [3].The relationship between the volume fractions of the ceramic and metal is defined as  m +  c = 1, with  c is expressed as: where  c and  m are the volume fractions of the ceramic and the metal, respectively, and  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of the FGM plate can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's ratio, and mass density of the ceramic and the metal, respectively., , and  are the equivalent elastic modulus, Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence, allowing for high-precision approximation of target functions with few truncation terms [52].In this study, the first kind of Chebyshev polynomials is employed as the spectral basis to expand the displacement field of the FGM plate in the thickness direction as the following spectral series: where  = [, , ] T and  ,  , and  are displacements in the  ,  , and  directions, respectively. and  are the truncation numbers for the transverse and the in-plane displacement fields, respectively, with the condition that  ≤ .{  } =0  are the generalized displacements defined as follows: are the volume fractions of the ceramic and the metal, respectively, and g is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material properties of the FGM plate can be computed using the following expression:

Functionally Graded Plate
Consider an FGM plate composed of ceramic a  represent the in-plane directions, and  denotes t where  c and  m are the volume fractions of the ce  is the gradient index [51].
Based on the rule of mixtures [2], the equivalen can be computed using the following expression: where  c ,  c ,  c and  m ,  m ,  m are the elastic modu of the ceramic and the metal, respectively., , an Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantag high-precision approximation of target functions with the first kind of Chebyshev polynomials is employed placement field of the FGM plate in the thickness dire where  = [, , ] T and  ,  , and  are displacem spectively. and  are the truncation numbers fo placement fields, respectively, with the condition displacements defined as follows:

Functionally Graded Plate
Consider an FGM plate composed of cerami  represent the in-plane directions, and  denote where  c and  m are the volume fractions of the  is the gradient index [51].
Based on the rule of mixtures [2], the equiva can be computed using the following expression where  c ,  c ,  c and  m ,  m ,  m are the elastic mo of the ceramic and the metal, respectively., , Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advant high-precision approximation of target functions w the first kind of Chebyshev polynomials is employ placement field of the FGM plate in the thickness d where  = [, , ] T and  ,  , and  are displac spectively. and  are the truncation numbers placement fields, respectively, with the conditio displacements defined as follows: where E c , ν c , ρ c and E m , ν m , ρ m are the elastic modulus, Poisson's ratio, and mass density of the ceramic and the metal, respectively.E, ν, and ρ are the equivalent elastic modulus, Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectral convergence, allowing for high-precision approximation of target functions with few truncation terms [52].In this study, the first kind of Chebyshev polynomials is employed as the spectral basis to expand the displacement field of the FGM plate in the thickness direction as the following spectral series: where u = [u, v, w] T and u, v, and w are displacements in the x, y, and z directions, respectively.M and N are the truncation numbers for the transverse and the in-plane displacement fields, respectively, with the condition that M ≤ N. {u i } N i=0 are the generalized displacements defined as follows: where {u i } N i=0 , {v i } N i=0 , and {w i } M−1 i=0 represent the basic unknowns, which are functions of the in-plane coordinates (x, y).{κ i (z)} N i=1 are the spectral basis excluding the first constant term [53].( ) represents the first-order derivative, ( ) ,x and ( ) ,y indicate partial derivatives with respect to x and y, respectively.
According to the d'Alembert principle and the principle of virtual work, the integral governing equations for the free vibration of the FGM plate are written as: By combining Equations ( 3)- (7), Equation ( 8) can be rearranged into the following two-dimensional (2D) format: where D and Λ represent the modulus matrix and inertia matrix, respectively, ω denotes the angular frequency, and

Discrete Governing Equations
The NURBS-based isogeometric method fulfills the C 1 -continuity requirement for discretizing the generalized displacements.In the case of a surface geometry model of the FGM plate, the model is discretized using 2D NURBS basis functions [48] as described below: where i = (i 1 , i 2 ), with i 1 and i 2 representing the sequence numbers of the ζ and η dimen- sions of the NURBS basis functions, respectively, i denotes the lexicographical order associated with i [55].x i represent the geometric control points, Q is the number of NURBS basis functions, and Ω is the mapping of the plate surface Ω in the 2D parameter space.R i correspond to the 2D NURBS basis functions, as described in Appendix B. Their further application in the discretization of the in-plane displacement fields is as follows: where q = [u 0 , . . . ,u N , v 0 , . . . ,v N , w 0 , . . . ,w M−1 ] T is a vector that contains all the basic unknowns.q i are the control vectors, and the displacement and the strain fields can be expressed as: where U i and E i are sparse block coefficient matrices, as described in Appendix C for the detailed expressions.
By combining Equation (13) with Equation ( 9) and rearranging, the discretized governing equations for the free vibration analysis of the FGM plate can be derived as follows: where J a represents the Jacobian determinant [56], K i j represent the block stiffness matrices, and M i j represent the block mass matrices.Equation ( 14) can be reformulated in a global form as follows: where K, M, and d represent the global stiffness matrix, global mass matrix, and displacement vector, respectively, which are assembled from K i j , M i j , and q j .The local support property of the NURBS basis functions ensures the sparsity of the stiffness matrix K and the mass matrix M [48].

Boundary Conditions and Multiple Patches Coupling
NURBS basis functions generally exhibit noninterpolatory property, which limits their direct applicability for enforcing boundary conditions and coupling adjacent NURBS patches [57].A general approach using penalty methods is employed in this paper to impose boundary conditions and couple multiple patches [57].For an FGM plate with midplane geometry represented by 2D NURBS basis functions, as described in Equation ( 11), the boundary is typically composed of segmented isoparametric curves.As shown in Figure 2, the constrained boundary ∂S is denoted as follows: where ς represents the parameter coordinate corresponding to the boundary, [ς 1 , ς 2 ] denotes its range of variation.The ` signifies the parameter functions with fixed parameter coordinate values.
where  ⌈⌉ are sparse block matrices of size 3 × ( + 2 + 2), and the definition of their nonzero elements can be obtained by combining Equation (10) and referring to Equation (A7) in Appendix C. The global matrix  is assembled from the block matrices  ⌈⌉ in lexicographical order.Substituting the parameter coordinates value corresponding to the boundary into Equation ( 17) yields the following expression: By applying the Galerkin method, the homogeneous boundary conditions can be discretized to yield the following expression: where  c = ‖d̀d ⁄ ‖ represents the arc length coefficient, and ‖ ‖ denotes the Euclidean norm [56]. is a symmetric sparse matrix and it has the same size as the stiffness matrix in the governing equations.By multiplying Equation ( 19) by a penalty coefficient [58,59] and directly substituting it into Equation ( 15), the boundary conditions can be enforced.The appropriate value of the penalty coefficient should be determined through numerical testing.
The penalty method can be extended to couple adjacent NURBS patches.As shown in Figure 3, adjacent NURBS patches are coupled by satisfying the following conditions: Due to the linear independence of the spectral basis in the thickness direction, as shown in Equation ( 3), constraining the displacement field is equivalent to imposing constraints on the generalized displacement field {u i } N i=0 individually.The discretized constrained displacement field can be expressed as: where L i are sparse block matrices of size 3 × (M + 2N + 2), and the definition of their nonzero elements can be obtained by combining Equation (10) and referring to Equation (A7) in Appendix C. The global matrix L is assembled from the block matrices L i in lexicographical order.Substituting the parameter coordinates value corresponding to the boundary into Equation ( 17) yields the following expression: By applying the Galerkin method, the homogeneous boundary conditions can be discretized to yield the following expression: where J c = dx/dς represents the arc length coefficient, and denotes the Euclidean norm [56].Ω is a symmetric sparse matrix and it has the same size as the stiffness matrix in the governing equations.
By multiplying Equation ( 19) by a penalty coefficient [58,59] and directly substituting it into Equation ( 15), the boundary conditions can be enforced.The appropriate value of the penalty coefficient should be determined through numerical testing.
The penalty method can be extended to couple adjacent NURBS patches.As shown in Figure 3, adjacent NURBS patches are coupled by satisfying the following conditions: where ∂S a and ∂S b represent the coupling boundaries of patch a and patch b, ς a and ς b are the parameter coordinates at ∂S a and ∂S b , respectively.Assuming that patch a and patch b have the same truncation order as the SDF, { ùa i } N i=0 and ùb  By combining Equations ( 3)-( 6) and Equation (10), and referring to Appendix C, the discretization of Equation (20b,c) follows a similar form as Equation ( 19).Subsequently, an appropriate penalty coefficient is chosen [58,59] to enable the coupling of NURBS patches.

Numerical Results and Discussion
The NURBS midplane depicted in Figure 4 includes a square plate, a circular plate, a square plate with a complicated cutout, and an L-shaped plate with multiple holes.The square plate modeled uses a single NURBS patch, the circular plate is composed of 5 NURBS patches, the square plate with a complicated cutout consists of 8 NURBS patches, and the L-shaped plate with multiple holes comprises 18 NURBS patches.The grid lines forming the mesh are represented by solid red lines, and the control points are denoted by solid blue dots.

Numerical Results and Discussion
The NURBS midplane depicted in Figure 4 includes a square plate, a circular plate, a square plate with a complicated cutout, and an L-shaped plate with multiple holes.The square plate modeled uses a single NURBS patch, the circular plate is composed of 5 NURBS patches, the square plate with a complicated cutout consists of 8 NURBS patches, and the L-shaped plate with multiple holes comprises 18 NURBS patches.The grid lines forming the mesh are represented by solid red lines, and the control points are denoted by solid blue dots.
Using the aforementioned four models as examples, this section compares the obtained numerical results with existing literature to validate the accuracy and generality of the proposed SDF-based IGA for the free vibration analysis of FGM plates.Numerical tests indicate that fourth-order NURBS basis functions and truncation numbers M = 5, N = 6 in the SDF-based IGA provide converged response results.In all cases, (p + 1) × (q + 1) Gaussian integration points are used for the in-plane elements [48].The number of integration points in the thickness direction is determined according to the value of the gradient index.The penalty coefficient is chosen as 10 4 times the maximum element value of the relevant stiffness matrix.The material properties of the components involved in the examples are listed in Table 1.The code is compiled using MATLAB 9.10 and executed on a computer with an Intel(R) Core (TM) i7-9750H CPU (2.60 GHz), Windows 10 64-bit operating system, 32 GB RAM, and 12 threads.The equipment was sourced from Lenovo, headquartered in Beijing, China.
The NURBS midplane depicted in Figure 4 includes a square plate, a circular plate, a square plate with a complicated cutout, and an L-shaped plate with multiple holes.The square plate modeled uses a single NURBS patch, the circular plate is composed of 5 NURBS patches, the square plate with a complicated cutout consists of 8 NURBS patches, and the L-shaped plate with multiple holes comprises 18 NURBS patches.The grid lines forming the mesh are represented by solid red lines, and the control points are denoted by solid blue dots.Using the aforementioned four models as examples, this section compares the obtained numerical results with existing literature to validate the accuracy and generality of the proposed SDF-based IGA for the free vibration analysis of FGM plates.Numerical tests indicate that fourth-order NURBS basis functions and truncation numbers  = 5， = 6 in the SDF-based IGA provide converged response results.In all cases, ( + 1) × ( + 1) Gaussian integration points are used for the in-plane elements [48].The number of integration points in the thickness direction is determined according to the value of the gradient index.The penalty coefficient is chosen as 10 4 times the maximum element value of the relevant stiffness matrix.The material properties of the components involved in the examples are listed in Table 1.The code is compiled using MATLAB 9.10 and executed on a computer with an Intel(R) Core (TM) i7-9750H CPU (2.60 GHz), Windows 10 64-bit operating system, 32 GB RAM, and 12 threads.The equipment was sourced from Lenovo, headquartered in Beijing, China.Consider a simply supported FGM square plate consisting of Al/Al 2 O 3 with a powerlaw distribution, as shown in Figure 4a.The material properties of the ceramic and metal are provided in Table 1, and the side length of the plate is 1.For comparison purposes, the dimensionless frequency is defined as ω = ωh ρ c /E c , where ω represents the angular frequency, h is the thickness of the plate, ρ c is the density of the ceramic material, and E c is the elastic modulus of the ceramic material.
As shown in Table 2, the dimensionless natural frequency results for Al/Al 2 O 3 square plates obtained using the SDF-based IGA are compared with the results from existing literature for different aspect ratios (a/h) and gradient indices (g).The quasi-3D solutions provided by Matsunaga [26] are obtained based on Carrera's unified formulation.The HSDT results given by Thai [60] are obtained based on a simple higher-order shear defor-mation theory.The FSDT results [61] correspond to analytical solutions based on first-order shear deformation theory.It can be observed that the SDF-based IGA results conform very well with the quasi-3D results [26].The FSDT results [61] slightly overestimate the frequencies of FGM plates at a high value of gradient index due to the use of a constant SCF for any values of gradient index.Conversely, the HSDT results [60] slightly underestimate frequencies since the HSDT neglects the thickness stretching effect.This highlights the necessity for considering the thickness stretching effect in the analysis of plate bending vibrations.When the thickness decreases, the first three bending frequencies significantly decrease.Additionally, when the gradient index increases, the frequencies also decrease.However, the influence of the gradient index on the frequencies is not as significant as the Influence of thickness.Figure 5 illustrates the corresponding bending vibration modes of the Al/Al 2 O 3 square plate for the aspect ratio a/h = 10 and the gradient index g = 10.The 3D-mode shapes provide a clear visual representation of the stretching effect in the plate's thickness.deformation theory.The FSDT results [61] correspond to analytical solutions based on first-order shear deformation theory.It can be observed that the SDF-based IGA results conform very well with the quasi-3D results [26].The FSDT results [61] slightly overestimate the frequencies of FGM plates at a high value of gradient index due to the use of a constant SCF for any values of gradient index.Conversely, the HSDT results [60] slightly underestimate frequencies since the HSDT neglects the thickness stretching effect.This highlights the necessity for considering the thickness stretching effect in the analysis of plate bending vibrations.When the thickness decreases, the first three bending frequencies significantly decrease.Additionally, when the gradient index increases, the frequencies also decrease.However, the influence of the gradient index on the frequencies is not as significant as the influence of thickness.Figure 5 illustrates the corresponding bending vibration modes of the Al/Al2O3 square plate for the aspect ratio  ℎ ⁄ = 10 and the gradient index  = 10.The 3D-mode shapes provide a clear visual representation of the stretching effect in the plate's thickness.

Circular Plate (Al/Al2O3)
Consider a clamped FGM circular plate composed of aluminum (Al) and aluminum oxide (Al2O3), as shown in Figure 4b.For the purpose of comparison, a power-law distribution is adopted in this example, given by  m = (1 2 ⁄ −  ℎ ⁄ )  ,  c = 1 −  m and the radius is 1.The dimensionless frequency is defined as  ̅ = 100ℎ√ c  c ⁄ , where ℎ is the thickness of the plate.
In Table 3, the dimensionless natural frequency results for Al/Al2O3 circular plates are listed for various thickness-to-radius ratios with a gradient index of 1.The converged results obtained using the SDF-IGA method are compared with the semianalytical solution based on the first-order shear deformation theory (FSDT) [62], the Abaqus results based

Functionally Graded Plate
Consider an FGM plate composed of ceramic and metal, w  represent the in-plane directions, and  denotes the thicknes where  c and  m are the volume fractions of the ceramic and  is the gradient index [51].
Based on the rule of mixtures [2], the equivalent material p can be computed using the following expression: (, , ) = ( c ,  c ,  c ) c + ( m ,  m ,  m where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisso of the ceramic and the metal, respectively., , and  are the e Poisson's ratio, and mass density, respectively.

Spectral Displacement Formulation
Chebyshev polynomials possess the advantage of spectra high-precision approximation of target functions with few trunca the first kind of Chebyshev polynomials is employed as the spec placement field of the FGM plate in the thickness direction as the where  c and  m are the volume  is the gradient index [51].
Based on the rule of mixture can be computed using the follow (, , ) where  c ,  c ,  c and  m ,  m ,  m ar of the ceramic and the metal, resp Poisson's ratio, and mass density

Spectral Displacement Formula
Chebyshev polynomials poss high-precision approximation of ta the first kind of Chebyshev polyno placement field of the FGM plate in and the radius is 1.The dimensionless frequency is defined as ω = 100ωh ρ c /E c , where h is the thickness of the plate.
In Table 3, the dimensionless natural frequency results for Al/Al 2 O 3 circular plates are listed for various thickness-to-radius ratios with a gradient index of 1.The converged results obtained using the SDF-IGA method are compared with the semianalytical solution based on the first-order shear deformation theory (FSDT) [62], the Abaqus results based on 3D elasticity theory (Abaqus) [63], together with the IGA solutions based on the modified firstorder shear deformation theory (sFSDT-IGA) [64] and the higher-order shear deformation theory (HSDT-IGA) [65].These comparisons demonstrate good agreement between the SDF-IGA results and the reference solutions for both thick and thin FGM plates.Reducing the thickness of the circular plate leads to a significant decrease in the plate's natural frequencies.The first six mode shapes of the Al/Al 2 O 3 circular plates with a thickness-toradius ratio h/r = 0.1 are shown in Figure 6.The sixth bending mode shape exhibits an irregular shape.on 3D elasticity theory (Abaqus) [63], together with the IGA solutions based on the modified first-order shear deformation theory (sFSDT-IGA) [64] and the higher-order shear deformation theory (HSDT-IGA) [65].These comparisons demonstrate good agreement between the SDF-IGA results and the reference solutions for both thick and thin FGM plates.Reducing the thickness of the circular plate leads to a significant decrease in the plate's natural frequencies.The first six mode shapes of the Al/Al2O3 circular plates with a thickness-to-radius ratio ℎ  ⁄ = 0.1 are shown in Figure 6.The sixth bending mode shape exhibits an irregular shape.A symmetric Al/ZrO 2 square plate with a complicated cutout is considered to verify the accuracy of the proposed SDF-IGA method in the analysis of free vibrations of irregularshaped structures.The geometric model is shown in Figure 4c, where the outer edges are either simply supported or clamped.The length-to-thickness ratio is chosen as h/L = 0.05, and the other detailed original geometric information can be found in the paper by Nguyen and Nguyen-Xuan [66].For observation and comparison, the dimensionless frequency is defined as ω = ωL 2 ρ c /E c /h .The first four natural frequency results for the square plates with multiple holes under thickness h = 0.1 and h = 0.2 are listed in Tables 6 and 7, respectively.The reference results are obtained by ANSYS software, utilizing Shell181 elements based on the FSDT and employing 100 uniform layers in the thickness direction.It can be observed that the variation in the gradient index significantly affects the dynamic characteristics of the Al/SiC L-shaped plates with multiple holes, with an increase in the gradient index resulting in a notable decrease in the fundamental frequency of the plates.Some slight differences between the results obtained in this study and those from ANSYS can be attributed to the use of different shear deformation theories.Moreover, as the gradient index increases, the frequencies decrease.Significant differences exist in the frequencies of homogeneous plates (g = 0) and FGM plates (g = 0).Figure 8 displays the first four mode shapes of the plate for the case of h = 0.1 and g = 4.The mode shapes exhibit the expected shape characteristics of symmetrical structures.
sults from the SDF-based IGA show good agreement with the reference solutions for different gradient indices.Increasing the gradient index decreases the natural frequencies of the plates, but the effect is insignificant.Additionally, when the boundary conditions change from simply supported to clamped, the natural frequencies of the plates significantly increase.The first six mode shapes of the square plate with a complicated cutout under the clamped boundary condition, with a gradient index equal to 1, are displayed in Figure 7.The mode shapes exhibit noticeable symmetry.Table 6.The dimensionless frequencies ω of the Al/SiC L-shaped plates with multiple holes (h = 0.2).ing in a notable decrease in the fundamental frequency of the plates.Some slight differences between the results obtained in this study and those from ANSYS can be attributed to the use of different shear deformation theories.Moreover, as the gradient index increases, the frequencies decrease.Significant differences exist in the frequencies of homogeneous plates ( = 0 ) and FGM plates ( ≠ 0 ). Figure 8 displays the first four mode shapes of the plate for the case of ℎ = 0.1 and  = 4.The mode shapes exhibit the expected shape characteristics of symmetrical structures.

Conclusions
In the classical plate theory (CPT), there is no independent variable for the angles of rotation, which naturally avoids the shear-locking phenomenon.The spectral displacement formulation (SDF) is an extension of the CPT and naturally possesses the advantage of avoiding shear locking.Chebyshev polynomials offer several advantages, such as orthogonality and spectral convergence, which make them well suited for approximating complex functions with high accuracy using a small number of terms.Combining the Chebyshev polynomials and complete 3D elasticity constitutive equations ensures the high precision of the SDF.Isogeometric analysis (IGA) is incorporated to analyze the free vibrations of functionally graded material (FGM) plates with complex geometries and different boundary conditions.
Numerical examples involving various types of plates, such as square plates, circular plates, square plates with complicated cutouts, and L-shaped plates with multiple holes, have been utilized to validate the effectiveness of the SDF-based IGA method in analyzing FGM plates.The results demonstrate that the proposed method accurately calculates the free vibrations of both thick and thin plates.The investigations in this study have also highlighted the significant influences of thickness and boundary conditions on the natural frequencies of the plates.Additionally, the gradient index of the FGM plates has been shown to affect the natural frequencies.These findings further emphasize the importance of considering these factors in the design and analysis of FGM plates.
With its demonstrated accuracy, the proposed SDF-based IGA method offers a practical and efficient approach for predicting the natural frequencies and mode shapes of FGM plates.Future research can expand upon this method by considering additional factors such as porosity, voids, and microstructural defects, as they may impact the vibrational behavior of FGM plates.Furthermore, generalizing this method to analyze nonlinear, forced vibrations, and other dynamic behaviors of FGM plates can provide valuable insights for designing and optimizing structures in various engineering applications.
The NURBS surface defined in the 2D parameter space, using the NURBS basis functions R p 1 ,p 2 i 1 ,i 2 (ζ, η), can be expressed as [48]: where S i 1 ,i 2 represents the geometric control points, and where R i (ζ, η) represent the 2D NURBS basis functions rearranged in lexicographical order [55].

Appendix C. Matrices U i and E i
U i and E i are block coefficient matrices, and the nonzero elements of these matrices are listed as follows: and (A8)

(
, , ) = ( c ,  c ,  c ) c + ( m ,  m ,  m ) m ( where  c ,  c ,  c and  m ,  m ,  m are the elastic modulus, Poisson's ratio, and mass densit of the ceramic and the metal, respectively., , and  are the equivalent elastic modulu Poisson's ratio, and mass density, respectively.

Figure 1 .
Figure 1.Functionally graded plate.The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic materi surface ( = −ℎ 2⁄ ) transitions from ceramic to metal using a power-The relationship between the volume fractions of the ceramic and met  c = 1, with  c is expressed as:

and
Mathematics 2023, 11, x FOR PEER REVIEW 4 of 20

Figure 1 .
Figure 1.Functionally graded plate.The top surface of the plate ( = ℎ 2 ⁄ ) consists surface ( = −ℎ 2 ⁄ ) transitions from ceramic to me The relationship between the volume fractions of th  c = 1, with  c is expressed as:

Figure 1 .
Figure 1.Functionally graded plate.The top surface of the plate ( = ℎ 2 ⁄ ) consi surface ( = −ℎ 2 ⁄ ) transitions from ceramic to m The relationship between the volume fractions of  c = 1, with  c is expressed as:

Figure 2 .
Figure 2. Illustration of the NURBS boundary.Due to the linear independence of the spectral basis in the thickness direction, as shown in Equation (3), constraining the displacement field is equivalent to imposing constraints on the generalized displacement field {  } =0  individually.The discretized constrained displacement field can be expressed as: () = ∑  ⌈⌉ () ⌈⌉  ⌈⌉=1

Figure 2 .
Figure 2. Illustration of the NURBS boundary.

20 {
of the generalized displacement field at ∂S a and ∂S b , respectively.Similarly, { ὲa i } N i=0 and ὲb i N i=0 are the values of the generalized strain field at ∂S a and ∂S b , respectively.Mathematics 2023, 11, x FOR PEER REVIEW 8 of  a ∋ ̀a( a ) = ̀b( b ) ∈  b ̀ a ( a ) = ̀ b ( b ), ∀ ∈ [0, ] ̀ a ( a ) = ̀ b ( b ), ∀ ∈ [0, ] , ∀ a =  b ∈ [̅ 1 , ̅ 2 ] where  a and  b represent the coupling boundaries of patch a and patch b,  a and  b are the parameter coordinates at  a and  b , respectively.Assuming that patch a and patch b have the same truncation order as the SDF, {̀ a } =0  and {̀ b } =0  are the values of the generalized displacement field at  a and  b , respectively.Similarly, {̀ a } =0  and {̀ b } =0  are the values of the generalized strain field at  a and  b , respectively.

Figure 3 .
Figure 3. Illustration of the NURBS patches coupling.

Figure 3 .
Figure 3. Illustration of the NURBS patches coupling.

3. 2 .
Circular Plate (Al/Al 2 O 3 ) Consider a clamped FGM circular plate composed of aluminum (Al) and aluminum oxide (Al 2 O 3 ), as shown in Figure 4b.For the purpose of comparison, a power-law distribution is adopted in this example, given by Mathematics 2023, 11, x FOR PEER REVIEW

Figure 1 .
Figure 1.Functionally graded plate.The top surface of the plate ( = ℎ 2 ⁄ ) consists of ceramic surface ( = −ℎ 2 ⁄ ) transitions from ceramic to metal using a The relationship between the volume fractions of the ceramic a  c = 1, with  c is expressed as:

Table 4 .
The dimensionless frequencies  ̅ of the Al/ZrO2 square plates with a complicated cutout under a simply supported boundary condition (ℎ/ = 0.05).

Table 1 .
Material properties of ceramics and metals.
3.1.Square Plate (Al/Al2O3) Consider a simply supported FGM square plate consisting of Al/Al2O3 with a power-

Table 1 .
Material properties of ceramics and metals.

Table 2 .
The dimensionless frequencies ω of the Al/Al 2 O 3 square plates.

Table 2 .
The dimensionless frequencies  ̅ of the Al/Al2O3 square plates.

Table 5 .
The dimensionless frequencies ω of the Al/ZrO 2 square plates with a complicated cutout under a clamped boundary condition (h/L = 0.05).

Table 7 .
The dimensionless frequencies ω of the Al/SiC L-shaped plates with multiple holes (h = 0.1).