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

Abstract

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 field 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.

1. 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 nth-order 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⁰ 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¹ 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.

2. The Model and Discretization

2.1. 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).

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 ${𝓋}_{\mathrm{m}}+{𝓋}_{\mathrm{c}}=1$, with ${𝓋}_{\mathrm{c}}$ is expressed as:

$$\begin{array}{c}\begin{array}{c}{𝓋}_{\mathrm{c}}={\left(1/2+z/h\right)}^g\end{array}, z\in \left[-h/2,h/2\right] \end{array}$$

(1)

where ${𝓋}_{\mathrm{c}}$ and ${𝓋}_{\mathrm{m}}$ 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:

$$\begin{array}{c}\begin{array}{c}\left(E,\nu ,\rho \right)=\left(E_{\mathrm{c}},{\nu }_{\mathrm{c}},{\rho }_{\mathrm{c}}\right){𝓋}_{\mathrm{c}}+\left(E_{\mathrm{m}},{\nu }_{\mathrm{m}},{\rho }_{\mathrm{m}}\right){𝓋}_{\mathrm{m}}\end{array} \end{array}$$

(2)

where $E_{\mathrm{c}}$, ${\nu }_{\mathrm{c}}$, ${\rho }_{\mathrm{c}}$ and $E_{\mathrm{m}}$, ${\nu }_{\mathrm{m}}$, ${\rho }_{\mathrm{m}}$ are the elastic modulus, Poisson’s ratio, and mass density of the ceramic and the metal, respectively. $E$, $\nu$, and $\rho$ are the equivalent elastic modulus, Poisson’s ratio, and mass density, respectively.

2.2. 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:

$$\begin{array}{c}u=u_0+{{\displaystyle \sum }}_{i=1}^M{A}_i\left(z\right)u_i+{{\displaystyle \sum }}_{i=M+1}^N{A}_i\left(z\right)u_i\end{array}$$

(3)

where $u={\left[u,v,w\right]}^{\mathrm{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\le N$. ${\left\{u_i\right\}}_{i=0}^N$ are the generalized displacements defined as follows:

$$\begin{array}{l}{u}_0={\left[u_0,v_0,0\right]}^{\mathrm{T}}, A_i\left(z\right)=\mathrm{diag}\left[{\kappa }_i\left(z\right),{\kappa }_i\left(z\right),{\kappa }_i^{\prime }\left(z\right)\right]\hfill \\ u_i=\{\begin{array}{l}\begin{array}{c}{[u_i-w_{i-1,x},v_i-w_{i-1,y},w_{i-1}]}^{\mathrm{T}}\end{array}, 1\le i\le M\hfill \\ {\left[u_i,v_i,0\right]}^{\mathrm{T}}, \left(M+1\right)\le i\le N\hfill \end{array}\hfill \end{array}$$

(4)

where ${\left\{u_i\right\}}_{i=0}^N$, ${\left\{v_i\right\}}_{i=0}^N$, and ${\left\{w_i\right\}}_{i=0}^{M-1}$ represent the basic unknowns, which are functions of the in-plane coordinates $\left(x,y\right)$. ${\left\{{\kappa }_i\left(z\right)\right\}}_{i=1}^N$are the spectral basis excluding the first constant term [53]. ${( )}^{\prime }$ represents the first-order derivative, ${( )}_{,x}$ and ${( )}_{,y}$ indicate partial derivatives with respect to $x$ and $y$, respectively.

The relationship between strains and displacements can be expressed as [54]:

$$\begin{array}{c}\epsilon ={\epsilon }_0+{{\displaystyle \sum }}_{i=1}^M{B}_i\left(z\right){\epsilon }_i+{{\displaystyle \sum }}_{i=M+1}^N{B}_i\left(z\right){\epsilon }_i\end{array}$$

(5)

where $\epsilon ={\left[{\epsilon }_{xx},{\epsilon }_{yy},{\epsilon }_{zz},{\gamma }_{yz},{\gamma }_{xz},{\gamma }_{xy}\right]}^{\mathrm{T}}$, ${\left\{{\epsilon }_i\right\}}_{i=0}^N$ are generalized strains in the following form:

$$\begin{array}{l}{\epsilon }_0={\left[u_{0,x},v_{0,y},0,0,0,u_{0,y}+v_{0,x}\right]}^{\mathrm{T}}\hfill \\ B_i\left(z\right)=\mathrm{diag}\left[{\kappa }_i\left(z\right),{\kappa }_i\left(z\right),{\kappa }_i^{″}\left(z\right),{\kappa }_i^{\prime }\left(z\right),{\kappa }_i^{\prime }\left(z\right),{\kappa }_i\left(z\right)\right]\hfill \\ {\epsilon }_i=\{\begin{array}{l}\begin{array}{l}[u_{i,x}-w_{i-1,xx},v_{i,y}-w_{i-1,yy},w_{i-1},\hfill \\ { v_i,u_i,u_{i,y}+v_{i,x}-2w_{i-1,xy}]}^{\mathrm{T}}\hfill \end{array}, 1\le i\le M\hfill \\ {\left[u_{i,x},v_{i,y},0,v_i,u_i,u_{i,y}+v_{i,x}\right]}^{\mathrm{T}}, \left(M+1\right)\le i\le N\hfill \end{array}\hfill \end{array}$$

(6)

where ${\left(\right)}^{″}$ represents the second-order derivative; ${\kappa }_i^{\prime }\left(z\right)$ and ${\kappa }_i^{″}\left(z\right)$ can be referenced in Appendix A.

Considering the complete 3D elasticity constitutive [54], the stress field can be expressed based on Hooke’s law as follows:

$$\begin{array}{l}\begin{array}{c}\begin{array}{c}\sigma =D\epsilon =D{\epsilon }_0+{{\displaystyle \sum }}_{i=1}^{M}D{B}_i\left(z\right){\epsilon }_i+{{\displaystyle \sum }}_{i=M+1}^{N}D{B}_i\left(z\right){\epsilon }_i\end{array}\end{array}\hfill \\ D=\left[\begin{array}{ll}{D}_{\mathrm{n}}\hfill & \mathbf{0}\hfill \\ \mathbf{0}\hfill & D_{\mathrm{s}}\hfill \end{array}\right], {\left(D_{\mathrm{n}}\right)}_{ij}=\lambda +2\mu {\delta }_{ij}, {\left(D_{\mathrm{s}}\right)}_{ij}=\mu {\delta }_{ij}\hfill \end{array}$$

(7)

where$\sigma ={\left[{\sigma }_{xx},{\sigma }_{yy},{\sigma }_{zz},{\sigma }_{yz},{\sigma }_{xz},{\sigma }_{xy}\right]}^{\mathrm{T}}$, $\lambda =E\nu /\left(1+\nu \right)/\left(1-2\nu \right)$, and $\mu =E/\left(1+\nu \right)/2$ are Lamé parameters. ${\delta }_{ij}$is the Kronecker symbol.

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:

$$\begin{array}{c}{{\displaystyle \int }}_{\mathcal{V}}\delta {\epsilon }^{\mathrm{T}}\sigma \mathrm{d}v={{\displaystyle \int }}_{\mathcal{V}}\delta u^{\mathrm{T}}\left(-\rho \ddot{u}\right)\mathrm{d}v\end{array}$$

(8)

By combining Equations (3)–(7), Equation (8) can be rearranged into the following two-dimensional (2D) format:

$${{\displaystyle \int }}_{\Omega }\delta {\tilde{\epsilon }}^{\mathrm{T}}\tilde{D}\tilde{\epsilon }\mathrm{d}a={\omega }^2{{\displaystyle \int }}_{\Omega }\delta {\tilde{u}}^{\mathrm{T}}\Lambda \tilde{u}\mathrm{d}a$$

(9)

where $\tilde{D}$ and $\Lambda$ represent the modulus matrix and inertia matrix, respectively, $\omega$ denotes the angular frequency, and

$$\begin{array}{l}\tilde{u}=\left\{\begin{array}{c}{u}_0\\ ⋮\\ u_N\end{array}\right\}, \begin{array}{c}\begin{array}{c}\tilde{\epsilon }=\left\{\begin{array}{c}{\epsilon }_0\\ ⋮\\ {\epsilon }_N\end{array}\right\}\end{array}\end{array}, \Lambda =\left[\begin{array}{ccc}{\Lambda }_{00}& \dots & {\Lambda }_{0N}\\ ⋮& \ddots & ⋮\\ {\Lambda }_{N0}& \dots & {\Lambda }_{NN}\end{array}\right], \tilde{D}=\left[\begin{array}{ccc}{\tilde{D}}_{00}& \cdots & {\tilde{D}}_{0N}\\ ⋮& \ddots & ⋮\\ {\tilde{D}}_{N0}& \dots & {\tilde{D}}_{NN}\end{array}\right]\hfill \\ \left({\Lambda }_{00},{\Lambda }_{0j},{\Lambda }_{ij}\right)={{\displaystyle \int }}_{-h/2}^{h/2}\rho \left(I_{3\times 3},A_j\left(z\right),A_i\left(z\right)A_j\left(z\right)\right)\mathrm{d}z, 1\le i,j\le N\hfill \\ \left({\tilde{D}}_{00},{\tilde{D}}_{0j},{\tilde{D}}_{ij}\right)={{\displaystyle \int }}_{-h/2}^{h/2}\left(D,D{B}_j\left(z\right),B_i\left(z\right)D{B}_j\left(z\right)\right)\mathrm{d}z, 1\le i,j\le N\hfill \\ \{\begin{array}{l}{\Lambda }_{ji}={\left({\Lambda }_{ij}\right)}^{\mathrm{T}}\hfill \\ {\tilde{D}}_{ji}={\left({\tilde{D}}_{ij}\right)}^{\mathrm{T}}\hfill \end{array}, \forall i,j\in \left[0,N\right];\hfill \end{array}$$

(10)

2.3. Discrete Governing Equations

The NURBS-based isogeometric method fulfills the C¹-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:

$$\Omega =\left\{x=\left(x,y\right)|x\left(\zeta \right)={{\displaystyle \sum }}_{\lceil i\rceil =1}^Q{x}_{\lceil i\rceil }{ℛ}_{\lceil i\rceil }\left(\zeta \right), \zeta =\left(\zeta ,\eta \right)\in \widehat{\Omega }\right\}$$

(11)

where $i=\left(i_1,i_2\right)$, with $i_1$ and $i_2$ representing the sequence numbers of the $\zeta$and $\eta$ dimensions of the NURBS basis functions, respectively, $\lceil i\rceil$ denotes the lexicographical order associated with $i$ [55]. $x_{\lceil i\rceil }$ represent the geometric control points, $Q$ is the number of NURBS basis functions, and $\widehat{\Omega }$ is the mapping of the plate surface $\Omega$ in the 2D parameter space. ${ℛ}_{\lceil i\rceil }$ 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:

$$q\left(\zeta \right)={{\displaystyle \sum }}_{\lceil i\rceil =1}^Q{q}_{\lceil i\rceil }{ℛ}_{\lceil i\rceil }\left(\zeta \right)$$

(12)

where $q={\left[u_0,\dots ,u_N,v_0,\dots ,v_N,w_0,\dots ,w_{M-1}\right]}^{\mathrm{T}}$ is a vector that contains all the basic unknowns. $q_{\lceil i\rceil }$ are the control vectors, and the displacement and the strain fields can be expressed as:

$$\begin{array}{l}\begin{array}{c}\begin{array}{c}\tilde{u}\left(\zeta \right)={{\displaystyle \sum }}_{\lceil i\rceil =1}^Q{U}_{\lceil i\rceil }\left(\zeta \right)q_{\lceil i\rceil }\end{array}\end{array}\hfill \\ \tilde{\epsilon }\left(\zeta \right)={{\displaystyle \sum }}_{\lceil i\rceil =1}^Q{{\rm E}}_{\lceil i\rceil }\left(\zeta \right)q_{\lceil i\rceil }\hfill \end{array}$$

(13)

where $U_{\lceil i\rceil }$ and $E_{\lceil i\rceil }$ 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:

$$\begin{array}{c}\begin{array}{l}{{\displaystyle \sum }}_{\lceil j\rceil =1}^Q{K}_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }{q}_{\lceil j\rceil }={\omega }^2{{\displaystyle \sum }}_{\lceil j\rceil =1}^Q{M}_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }{q}_{\lceil j\rceil }, \forall i\in \left[1,Q\right]\hfill \\ K_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }={{\displaystyle \int }}_{\widehat{\Omega }}{{\rm E}}_{\lceil i\rceil }^{\mathrm{T}}\tilde{D}{{\rm E}}_{\lceil j\rceil }{J}_{\mathrm{a}}\mathrm{d}\widehat{a}\hfill \\ M_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }={{\displaystyle \int }}_{\widehat{\Omega }}{U}_{\lceil i\rceil }^{\mathrm{T}}\Lambda U_{\lceil j\rceil }{J}_{\mathrm{a}}\mathrm{d}\widehat{a}\hfill \end{array}\end{array}$$

(14)

where $J_{\mathrm{a}}$ represents the Jacobian determinant [56], $K_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }$ represent the block stiffness matrices, and $M_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }$ represent the block mass matrices. Equation (14) can be reformulated in a global form as follows:

$$Kd={\omega }^{2}Md$$

(15)

where$K$, $M$, and $d$ represent the global stiffness matrix, global mass matrix, and displacement vector, respectively, which are assembled from $K_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }$, $M_{\lceil \mathbf{i}\rceil \lceil \mathbf{j}\rceil }$, and $q_{\lceil j\rceil }$. The local support property of the NURBS basis functions ensures the sparsity of the stiffness matrix $K$ and the mass matrix $M$ [48].

2.4. 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:

$$\begin{array}{c}𝜕S=\left\{\grave{x}=\left(\grave{x},\grave{y}\right)|\grave{x}\left(\varsigma \right)={{\displaystyle \sum }}_{\lceil i\rceil =1}^Q{x}_{\lceil i\rceil }{\grave{ℛ}}_{\lceil i\rceil }\left(\varsigma \right), \varsigma \in \left[{\overline{\varsigma }}_1,{\overline{\varsigma }}_2\right]\right\}\end{array}$$

(16)

where $\varsigma$ represents the parameter coordinate corresponding to the boundary, $\left[{\overline{\varsigma }}_1,{\overline{\varsigma }}_2\right]$ denotes its range of variation. The $\left({ }^{̀}\right)$signifies the parameter functions with fixed parameter coordinate values.

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 ${\left\{u_i\right\}}_{i=0}^N$ individually. The discretized constrained displacement field can be expressed as:

$$c\left(\zeta \right)={{\displaystyle \sum }}_{\lceil i\rceil =1}^Q{L}_{\lceil i\rceil }\left(\zeta \right)q_{\lceil i\rceil }=L\left(\zeta \right)d$$

(17)

where $L_{\lceil i\rceil }$ are sparse block matrices of size $3\times \left(M+2N+2\right)$, 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_{\lceil i\rceil }$ in lexicographical order.

Substituting the parameter coordinates value corresponding to the boundary into Equation (17) yields the following expression:

$$\begin{array}{c}\begin{array}{c}\begin{array}{c}\grave{c}\left(\varsigma \right)=\grave{L}\left(\varsigma \right)d\end{array}\end{array}\end{array}, \mathrm{on}𝜕S$$

(18)

By applying the Galerkin method, the homogeneous boundary conditions can be discretized to yield the following expression:

$$\Omega d=\mathbf{0}, \Omega =\begin{array}{ll}{{\displaystyle \int }}_{{\overline{\varsigma }}_1}^{{\overline{\varsigma }}_2}{\grave{L}}^{\mathrm{T}}\grave{L}{J}_{\mathrm{c}}\mathrm{d}\varsigma \hfill & \mathrm{Clamped}\hfill \\ {{\displaystyle \int }}_{{\overline{\varsigma }}_1}^{{\overline{\varsigma }}_2}{\grave{L}}^{\mathrm{T}}{𝕔}_{\mathrm{s}}^{\mathrm{T}}{𝕔}_{\mathrm{s}}\grave{L}{J}_{\mathrm{c}}\mathrm{d}\varsigma \hfill & \mathrm{Simply}\mathrm{supported}\hfill \end{array}$$

(19)

where $J_{\mathrm{c}}=\parallel \mathrm{d}\grave{x}/\mathrm{d}\varsigma \parallel$ represents the arc length coefficient, and $\parallel \parallel$ denotes the Euclidean norm [56]. $\Omega$ 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:

$$\{\begin{array}{ll}𝜕S^{\mathrm{a}}\ni {\grave{x}}^{\mathrm{a}}\left({\varsigma }^{\mathrm{a}}\right)={\grave{x}}^{\mathrm{b}}\left({\varsigma }^{\mathrm{b}}\right)\in 𝜕S^{\mathrm{b}}\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\left(20\mathrm{a}\right)\\ {\grave{u}}_i^{\mathrm{a}}\left({\varsigma }^{\mathrm{a}}\right)={\grave{u}}_i^{\mathrm{b}}\left({\varsigma }^{\mathrm{b}}\right), \forall i\in \left[0,N\right],\forall {\varsigma }^{\mathrm{a}}={\varsigma }^{\mathrm{b}}\in \left[{\overline{\varsigma }}_1,{\overline{\varsigma }}_2\right]\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\left(20\mathrm{b}\right)\\ {\grave{\epsilon }}_i^{\mathrm{a}}\left({\varsigma }^{\mathrm{a}}\right)={\grave{\epsilon }}_i^{\mathrm{b}}\left({\varsigma }^{\mathrm{b}}\right), \forall i\in \left[0,N\right]\hfill & \hfill \hspace{1em}\hspace{1em}\hspace{1em}\left(20\mathrm{c}\right)\end{array}$$

where $𝜕S^{\mathrm{a}}$ and $𝜕S^{\mathrm{b}}$ represent the coupling boundaries of patch a and patch b, ${\varsigma }^{\mathrm{a}}$ and ${\varsigma }^{\mathrm{b}}$ are the parameter coordinates at $𝜕S^{\mathrm{a}}$ and $𝜕S^{\mathrm{b}}$, respectively. Assuming that patch a and patch b have the same truncation order as the SDF, ${\left\{{\grave{u}}_i^{\mathrm{a}}\right\}}_{i=0}^N$ and ${\left\{{\grave{u}}_i^{\mathrm{b}}\right\}}_{i=0}^N$ are the values of the generalized displacement field at $𝜕S^{\mathrm{a}}$ and $𝜕S^{\mathrm{b}}$, respectively. Similarly, ${\left\{{\grave{\epsilon }}_i^{\mathrm{a}}\right\}}_{i=0}^N$ and ${\left\{{\grave{\epsilon }}_i^{\mathrm{b}}\right\}}_{i=0}^N$ are the values of the generalized strain field at $𝜕S^{\mathrm{a}}$ and $𝜕S^{\mathrm{b}}$, respectively.

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.

3. 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, $\left(p+1\right)\times \left(q+1\right)$ 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⁴ 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. Material properties of ceramics and metals.

	Young’s Modulus (Gpa)	Poisson’s Ratio	Density (kg/m3)
Al	70	0.3	2707
Al2O3	380	0.3	3800
ZrO2	270	0.3	5700
SiC	427	0.17	3210

. 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.

3.1. Square Plate (Al/Al₂O₃)

Consider a simply supported FGM square plate consisting of Al/Al₂O₃ with a power-law 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 $\overline{\omega }=\omega h\sqrt{{\rho }_{\mathrm{c}}/E_{\mathrm{c}}}$, where $\omega$ represents the angular frequency, $h$ is the thickness of the plate, ${\rho }_{\mathrm{c}}$ is the density of the ceramic material, and $E_{\mathrm{c}}$ is the elastic modulus of the ceramic material.

As shown in

Table 2. The dimensionless frequencies $\overline{\omega }$ of the Al/Al₂O₃ square plates.

$\mathit{g}$	Method	${\overline{\mathit{\omega }}}_1$(1, 1)		${\overline{\mathit{\omega }}}_2$(1, 2)		${\overline{\mathit{\omega }}}_3$(2, 2)
		$\mathit{a}\mathbf{/}\mathit{h}\mathbf{=}\mathbf{5}$	$\mathit{a}\mathbf{/}\mathit{h}\mathbf{=}\mathbf{10}$	$\mathit{a}\mathbf{/}\mathit{h}\mathbf{=}\mathbf{5}$	$\mathit{a}\mathbf{/}\mathit{h}\mathbf{=}\mathbf{10}$	$\mathit{a}\mathbf{/}\mathit{h}\mathbf{=}\mathbf{5}$	$\mathit{a}\mathbf{/}\mathit{h}\mathbf{=}\mathbf{10}$
1	Quasi-3D [26]	0.1640	0.0443	0.3644	0.1063	0.5444	0.1640
	HSDT [60]	0.1631	0.0442	0.3607	0.1059	0.5254	0.1631
	FSDT [61]	0.1631	0.0442	0.3604	0.1059	0.5245	0.1631
	Present	0.1640	0.0443	0.3643	0.1063	0.5324	0.1640
10	Quasi-3D [26]	0.1306	0.0364	0.2790	0.0859	0.3981	0.1306
	HSDT [60]	0.1301	0.0364	0.2771	0.0856	0.3948	0.1301
	FSDT [61]	0.1324	0.0366	0.2856	0.0867	0.4097	0.1324
	Present	0.1305	0.0364	0.2789	0.0859	0.3980	0.1305

, the dimensionless natural frequency results for Al/Al₂O₃ 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 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/Al₂O₃ 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.

3.2. Circular Plate (Al/Al₂O₃)

Consider a clamped FGM circular plate composed of aluminum (Al) and aluminum oxide (Al₂O₃), as shown in Figure 4b. For the purpose of comparison, a power-law distribution is adopted in this example, given by ${𝓋}_{\mathrm{m}}={\left(1/2-z/h\right)}^g,{𝓋}_{\mathrm{c}}=1-{𝓋}_{\mathrm{m}}$ and the radius is 1. The dimensionless frequency is defined as $\overline{\omega }=100\omega h\sqrt{{\rho }_{\mathrm{c}}/E_{\mathrm{c}}}$, where $h$ is the thickness of the plate.

In

Table 3. The dimensionless frequencies $\overline{\omega }$ of the Al/Al₂O₃ circular plates ($g=1$).

$\mathit{h}/\mathit{r}$	Method	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
0.2	FSDT [62]	8.6535	16.7666	25.6486	28.7574	34.0756	35.0981
	Abaqus [63]	8.6403	16.7890	25.7661	28.9152	34.1893	35.3618
	sFSDT-IGA [64]	8.6486	17.0016	26.2512	28.7691	34.1216	36.1557
	HSDT-IGA [65]	8.6787	16.8595	25.8479	29.0092	34.0581	35.4875
	Present	8.6916	16.8787	25.8735	29.0259	34.1437	35.4900
0.1	FSDT [62]	2.3053	4.6934	7.5146	8.5181	10.7128	12.6197
	Abaqus [63]	2.2888	4.6661	7.4808	8.4829	10.6776	12.5877
	sFSDT-IGA [64]	2.3040	4.7137	7.5773	8.5244	10.8524	12.7017
	HSDT-IGA [65]	2.3076	4.7005	7.5318	8.5380	10.7483	12.6636
	Present	2.3079	4.7021	7.5341	8.5421	10.7507	12.6709
0.01	FSDT [62]	0.0236	0.0491	0.0805	0.0918	0.1178	0.1404
	Abaqus [63]	0.0234	0.0486	0.0798	0.0909	0.1167	0.1391
	sFSDT-IGA [64]	0.0236	0.0491	0.0805	0.0919	0.1180	0.1408
	HSDT-IGA [65]	0.0236	0.0492	0.0807	0.0924	0.1191	0.1431
	Present	0.0236	0.0491	0.0805	0.0918	0.1178	0.1404

, the dimensionless natural frequency results for Al/Al₂O₃ 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 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/Al₂O₃ circular plates with a thickness-to-radius ratio $h/r=0.1$ are shown in Figure 6. The sixth bending mode shape exhibits an irregular shape.

3.3. Square Plate with a Complicated Cutout (Al/ZrO₂)

A symmetric Al/ZrO₂ 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 irregular-shaped 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 $\begin{array}{c}\overline{\omega }=\omega L^2\sqrt{{\rho }_c/E_c}/h\end{array}$.

The dimensionless natural frequency results for the square plates with a complicated cutout under simply supported or clamped boundary conditions are listed in

Table 4. The dimensionless frequencies $\overline{\omega }$ of the Al/ZrO₂ square plates with a complicated cutout under a simply supported boundary condition ($h/L=0.05$).

Mode.	Method	$\mathit{g}=0.5$	$\mathit{g}=1$	$\mathit{g}=5$	$\mathit{g}=20$	$\mathit{g}=50$
1	3D-IGA [66]	6.75	6.58	6.71	6.46	6.19
	3D-IGSBFEM [67]	6.7340	6.5896	6.7108	6.5198	6.2596
	Present	6.7695	6.5825	6.7086	6.5565	6.2618
2	3D-IGA [66]	11.00	10.73	10.88	10.48	10.07
	3D-IGSBFEM [67]	10.9747	10.7355	10.8690	10.5745	10.1763
	Present	11.0038	10.7032	10.8487	10.6100	10.3118
3	3D-IGA [66]	12.36	12.06	12.24	11.79	11.32
	3D-IGSBFEM [67]	12.3640	12.0954	12.2590	11.9235	11.4698
	Present	12.4138	12.0769	12.2524	11.9814	11.6415
4	3D-IGA [66]	19.84	19.35	19.60	18.89	18.15
	3D-IGSBFEM [67]	19.8134	19.3808	19.6133	19.0241	18.3701
	Present	19.8995	19.3596	19.6091	19.0445	18.6461
5	3D-IGA [66]	21.28	20.77	19.73	19.05	18.81
	3D-IGSBFEM [67]	21.2322	20.7737	19.6994	19.0857	18.8038
	Present	21.0078	20.7304	19.6922	19.1810	18.8638
6	3D-IGA [66]	21.33	20.92	21.00	20.25	19.48
	3D-IGSBFEM [67]	21.2879	20.8954	20.9776	20.4249	19.6844
	Present	21.3682	20.8862	20.9538	20.5002	19.9453
7	3D-IGA [66]	30.20	29.45	29.72	28.68	27.59
	3D-IGSBFEM [67]	30.1039	29.4393	29.6770	28.9066	27.8665
	Present	30.2406	29.4271	29.6927	29.0596	28.2820
8	3D-IGA [66]	31.76	30.97	31.28	30.18	29.03
	3D-IGSBFEM [67]	31.6935	30.9945	31.2640	30.4492	29.3458
	Present	31.7979	30.9716	31.2701	30.6023	29.7779
9	3D-IGA [66]	35.63	34.74	34.74	33.53	32.55
	3D-IGSBFEM [67]	35.6922	34.9013	34.6652	33.4726	33.0217
	Present	35.7754	34.8132	34.6536	33.5100	33.1911
10	3D-IGA [66]	37.49	36.80	35.03	33.82	33.10
	3D-IGSBFEM [67]	37.4047	36.7308	35.1432	33.9928	33.0836
	Present	37.0479	36.7255	35.0856	33.9272	33.4448

and

Table 5. The dimensionless frequencies $\overline{\omega }$ of the Al/ZrO₂ square plates with a complicated cutout under a clamped boundary condition ($h/L=0.05$).

Mode	Method	$\mathit{g}=0.5$	$\mathit{g}=1$	$\mathit{g}=5$	$\mathit{g}=20$	$\mathit{g}=50$
1	3D-IGA [66]	14.98	14.62	14.79	14.41	13.84
	3D-IGSBFEM [67]	14.9928	14.6649	14.8201	14.4257	13.8929
	Present	15.0061	14.5983	14.7673	14.4459	14.0490
2	3D-IGA [66]	25.79	25.17	25.38	24.74	23.79
	3D-IGSBFEM [67]	25.8817	25.3108	25.5007	24.8407	23.9522
	Present	25.9223	25.2352	25.4412	24.8977	24.2389
3	3D-IGA [66]	25.95	25.32	25.54	24.9	23.93
	3D-IGSBFEM [67]	26.0072	25.4339	25.6269	24.9629	24.0689
	Present	26.0312	25.3438	25.5532	25.0065	24.3441
4	3D-IGA [66]	31.44	30.68	30.83	30.07	28.95
	3D-IGSBFEM [67]	31.4801	30.7790	30.9001	30.1256	29.0898
	Present	31.5092	30.6864	30.8316	30.1862	29.4175
5	3D-IGA [66]	32.45	31.67	31.8	31.02	29.87
	3D-IGSBFEM [67]	32.4235	31.6999	31.8004	31.0085	29.9515
	Present	32.4472	31.5919	31.7224	31.0604	30.2748
6	3D-IGA [66]	39.04	38.1	38.16	37.23	35.9
	3D-IGSBFEM [67]	39.0165	38.1390	38.1562	37.2304	36.0014
	Present	39.0812	38.0654	38.1203	37.3388	36.4234
7	3D-IGA [66]	42.44	41.42	41.46	40.46	39.02
	3D-IGSBFEM [67]	42.5408	41.5822	41.5650	40.5645	39.2392
	Present	42.4714	41.5087	41.5359	40.6890	39.7011
8	3D-IGA [66]	46.02	44.91	44.87	43.8	42.27
	3D-IGSBFEM [67]	46.3903	45.3397	45.2177	44.1488	42.7473
	Present	46.0243	45.2794	45.2382	44.3260	43.2688
9	3D-IGA [66]	52.6	51.33	51.2	49.99	48.28
	3D-IGSBFEM [67]	52.5765	51.3800	51.1881	49.9993	48.4299
	Present	52.5841	51.2439	51.1134	50.0930	48.9238
10	3D-IGA [66]	53.99	52.69	52.5	51.26	49.53
	3D-IGSBFEM [67]	53.9887	52.7576	52.5155	51.3034	49.7107
	Present	54.0316	52.6317	52.4478	51.4045	50.2191

. These results are compared with the converged IGA results based on the 3D elasticity theory (3D-IGA) [66] and the semianalytical results obtained using the IGA and scaled boundary element method based on the 3D elasticity theory (3D-IGSBFEM) [67]. The results 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.

3.4. L-Shaped Plate with Multiple Holes (Al/SiC)

Consider a symmetric L-shaped plate with multiple holes made of aluminum (Al) and silicon carbide (SiC), as shown in Figure 4d. The geometric parameters of the L-shaped plate are defined as follows: $a=1$, $b=1.3$, and $r=\varphi =0.3$. The outer edges of the plate are clamped. For observation and comparison, the dimensionless frequency is defined as $\begin{array}{c}\stackrel{𝓋}{\omega }=\omega L^2\sqrt{{\rho }_c/E_c}/h\end{array}$, where $L$ is the length of the outer edges.

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

Table 6. The dimensionless frequencies $\overline{\omega }$ of the Al/SiC L-shaped plates with multiple holes ($h=0.2$).

Mode	Method	$\mathit{g}=0$	$\mathit{g}=2$	$\mathit{g}=4$	$\mathit{g}=6$	$\mathit{g}=8$	$\mathit{g}=10$
1	Ansys	11.1813	6.9654	6.4343	6.2463	6.1241	6.0254
	Present	11.2659	7.2427	6.7304	6.5612	6.4437	6.3450
2	Ansys	12.0884	7.4777	6.8902	6.6787	6.5424	6.4343
	Present	12.1871	7.7832	7.2145	7.0171	6.8902	6.7821
3	Ansys	16.5205	10.0533	9.1979	8.8924	8.6950	8.5493
	Present	16.6615	10.5421	9.7149	9.4282	9.2449	9.0992
4	Ansys	21.4461	12.7840	11.6231	11.1907	10.9322	10.7395
	Present	21.6858	13.6018	12.4738	12.0696	11.8205	11.6278

and

Table 7. The dimensionless frequencies $\overline{\omega }$ of the Al/SiC L-shaped plates with multiple holes ($h=0.1$).

Mode	Method	$\mathit{g}=0$	$\mathit{g}=2$	$\mathit{g}=4$	$\mathit{g}=6$	$\mathit{g}=8$	$\mathit{g}=10$
1	Ansys	11.4680	7.3884	6.8996	6.7398	6.6270	6.5236
	Present	11.5432	7.4636	6.9842	6.8244	6.7116	6.6176
2	Ansys	12.5302	7.9994	7.4636	7.2850	7.1628	7.0500
	Present	12.6054	8.0934	7.5576	7.379	7.2568	7.1440
3	Ansys	17.296	10.9792	10.1990	9.9358	9.7572	9.6068
	Present	17.4370	11.1202	10.3494	10.0956	9.9170	9.7666
4	Ansys	22.8514	14.5512	13.3104	16.6380	12.7088	12.5114
	Present	23.0864	14.6170	13.5642	16.9388	12.9814	12.7746

, 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\ne 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.

4. 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.