A Nonlinear Free Vibration Analysis of Functionally Graded Beams Using a Mixed Finite Element Method and a Comparative Artificial Neural Network

Abstract

Based on the Hamilton principle combined with the Timoshenko beam theory, the authors developed a mixed finite element (FE) method for the nonlinear free vibration analysis of functionally graded (FG) beams under combinations of simply supported, free, and clamped edge conditions. The material properties of the FG beam gradually and smoothly varied through the thickness direction according to the power-law distributions of the volume fractions in the constituents, and the effective material properties of the FG beam were estimated using the rule of mixtures. The von Kármán geometrical nonlinearity was considered. The FE solutions of the amplitude-frequency relations of the FG beam were obtained using an iterative process. Implementing the mixed FE method showed that its solutions converged rapidly and that the convergent solutions closely agreed with the accurate solutions reported in the literature. A multilayer perceptron (MP) back propagation neural network (BPNN) was also developed to predict the nonlinear free vibration behavior of the FG beam. After appropriate training, the prediction of the MP BPNN’s amplitude-frequency relations was entirely accurate compared to those obtained using the mixed FE method, and its central processing unit time was less time-consuming than that of the mixed FE method.

1. Introduction

In recent years, functionally graded (FG) material, which is composed of a two- or multi-phase material according to some spatial distributions of the volume fractions of these constituents, has gradually become an important industrial material and has been used to form a variety of beam-, plate- and shell-like structures in advanced engineering due to the flexibility of its application and the continuous distributions of its material properties over the physical domain of the structures in which it is used [1,2]. FG structures could be designed to improve their specific structural performances by changing the volume fractions of the constituents. In addition, FG structures could replace conventional fiber-reinforced composite (FRC) structures to prevent delamination failure, which often occurs at the interfaces between the adjacent layers of FRC structures due to the material properties that suddenly change at such places. Some review articles discussing the development, manufacture, and application of FG and FRC structures and their corresponding mechanical analyses can be found in the literature [3,4]. Among these, the review conducted in this work focuses on the literature related to the linear and nonlinear free vibration analyses of FG beams with different boundary conditions.

The linear structural analyses of FG beams have been presented. Based on the principle of virtual displacements (PVD) combined with the first-order shear deformation theory (FSDT), Chakraborty et al. [5,6] developed a new beam element for the static bending, thermo-elastic, free vibration, and wave propagation analyses of FG beams. The differences between the above-mentioned structural behaviors of the metal-ceramic two-phase FG beam and those of the pure metal and ceramic beams were highlighted. Aydogdu and Taskin [7] studied the free vibration responses of FGM beams using various equivalent single-layered theories (ESLTs), such as the classical beam theory, the FSDT, and the higher-order shear deformation theory (HSDT). Simsek [8,9] investigated the static bending, free vibration behaviors of FG beams under various boundary conditions using the HSDT, in which some effects on the deflection and frequency parameters of the FG beams were examined, such as the effects of slender ratios, material-property gradient indices, and different formulations. The HSDT was extended to the thermo-mechanical vibration analysis of sandwich beams with FG carbon nanotube-reinforced composite (CNTRC) face sheets by Ebrahimi and Farazmandnia [10]. Using Carrera’s unified formulation (CUF) [11], Pietro et al. [12] developed a hierarchical one-dimensional finite element method to study the thermo-elastic response of FG beams. Coskun et al. [13] developed a third-order plate theory for the static bending, free vibration, and static buckling analyses of FG porous microplates, where three different porosity distributions through the thickness direction of the plate were considered. Tham et al. [14] investigated the free vibration behavior of laminated FG CNTRC doubly curved shallow shells using a new four-variable refined theory. Wu and Xu [15] presented strong and weak formulations of a mixed higher-order formulation for the static bending analysis of FG beams when subjected to thermo-mechanical loads. Based on the mixed higher-order shear deformation theory, Wu and Li [16] examined the thermal buckling behavior of FG beams, and they also developed the multi-objective optimization of FG beams to maximize critical temperature change parameters and minimize the total mass of the FG beam using a non-dominated sorting-based genetic algorithm.

ESLTs combined with the von Kármán geometrical nonlinearity (VKGN) were also applied as extensions to FG beams’ geometrical nonlinear static and free vibration analyses. Ma and Lee [17] presented an exact, closed-form solution for the nonlinear static responses of single-layered FG beams under different boundary conditions using the FSDT and accounting for the VKGN. Based on the HSDT and using a physically neutral surface and the VKGN, Zhang [18] examined the nonlinear bending behavior of single-layered FG beams using the Ritz method. The HSDT was also extended to a bending analysis of sandwich beams with CNTRC face sheets by Salami [19] and a nonlinear bending analysis of FG CNTRC plates in thermal environments by Shen [20]. Based on the FSDT, Ghayesh [21] investigated the nonlinear forced vibration behavior of axially FG Timoshenko tapered beams. Shen et al. [22] presented the nonlinear vibration analysis of FG graphene-reinforced composite laminated beams when resting on elastic foundations in thermal environments using the HSDT. Ding et al. [23] applied the Euler-Bernoulli beam theory to the nonlinear vibration analysis of FG beams, in which the effects of the rotary inertia of the cross-section and the neural surface position were considered. Eltaher et al. [24] carried out the post-buckling and nonlinear vibration analyses of beams resting on a nonlinear foundation. Based on Hamilton’s principle, Chaudhari and Lal [25] carried out a nonlinear free vibration analysis of FG CNTRC beams, which were subjected to thermal loading, in which the HSDT and the VKGN were used to derive a weak-form formulation, from which a Lagrangian C⁰ element with four degrees-of-freedom per node was developed. Based on the FSDT and combined with the VKGN, Mirzaei and Kiani [26] analyzed the nonlinear free vibration characteristics of temperature-dependent sandwich beams with CNTRC face sheets. Based on the third-order shear deformation theory, Babaei et al. [27] investigated the small- and large-amplitude-free vibration behaviors of FG beams resting on an elastic foundation consisting of Winkler, shear, and nonlinear springs.

Over the past few decades, artificial neural networks (ANNs) have emerged as useful mathematical tools with a self-learning ability and have been used in various advanced engineering applications. The ANN approach was inspired by the biological neural system, such that its algorithm was designed as a multilayer structure in which each layer consists of some interconnected neurons. The ANN was trained using some given data sources, and the learning experience was stored in the weight numbers and biases, which were used to connect the relations between the outputs of the last layer and the inputs of each neuron in the current layer. Due to its superior abilities to handle massively parallel computation and self-learning, the ANN can be conventionally used in function approximations, pattern recognition, classification, etc. [28]. Recently, ANNs have further shown their diverse applications, such as machine learning [29], the modeling of structural and material behavior [30], and optimal design [31,32].

Chakraverty et al. [33] developed an ANN to estimate the vibration characteristics of plate structures, in which initial weight numbers were generated using regression analysis. Reddy et al. [34] used an ANN to predict the natural frequencies of laminated composite plates under clamped boundary conditions, in which a multilayer perceptron (MP) back propagation neural network (BPNN) was used, and the data sources obtained using the linear shell elements were provided for training. Jodaei et al. [35] developed a state space differential quadrature (SSDQ) method for the free vibration analysis of functionally graded (FG) annular plates under various boundary conditions, and these SSDQ solutions were used to train an ANN. Fetene et al. [36] developed a FEM-based neural network for the inverse prediction of bending a cantilever beam. The comparative modeling of the static and buckling behavior of laminated composite structures with the aid of ANNs has also been proposed by Subramani and Sharmila [37] and Liu et al. [38], respectively.

From the above literature survey, the authors found that the analyses mentioned above were rarely based on the weak-form formulation compared with those based on the strong-form formulation, and the corresponding equilibrium and motion equations derived using the Reissner mixed variation theorem (RMVT) were even fewer than those derived using the PVD, even though RMVT-based models were concluded to be superior to PVD-based models for the various analyses of plates and shells [11]. Hence, developing an RMVT-based weak-form formulation for analyzing FG beams is important in academic research and practical applications. In addition, as mentioned above, ANNs can be used to predict the structural behavior of FG beams by feeding them suitable data sources; thus, the development of a competitive ANN algorithm is necessary for future advanced studies with a lot of time-consuming characteristics, such as nonlinear modeling, optimal design, system recognition, etc.

In this work, the authors derived a weak-form formulation of an RMVT-based Timoshenko beam theory (TBT) for the large-amplitude free vibration analysis of FG beams under combinations of simply supported, free, and clamped edge conditions, in which the VKGN effect was considered. Nonlinear equilibrium equations of the finite element (FE) method based on the mixed TBT were derived using a variational approach, and an iterative process was used to obtain the numerical solutions for the issue of interest. The material properties were assumed to obey the power-law distributions that varied through the thickness direction of the FG beam according to the volume fractions of the constituents, and the effective material properties were estimated using the rules of mixtures. Four different boundary conditions, simple-simple (S-S), clamped-clamped (C-C), clamped-simple (C-S), and clamped-free (C-F), were considered. A parametric study of the critical effects on the amplitude-frequency relations in significant amplitude-free vibration cases was conducted, including geometrical nonlinearity, boundary conditions, aspect ratios, and material-property gradient indices. In addition, an MP BPNN was developed to predict the amplitude-frequency relations of FG beams, in which the Levenberg–Marquardt algorithm [28] was adopted to speed up the convergence rate of the MP BPNN. Some critical effects on the performance of the MP BPNN were closely examined, such as the number of layers used in the ANN and the number of neurons in each layer.

2. Weak-Form Formulation of the Nonlinear Mixed TBT

Based on the mixed TBT, a weak-form formulation for a large amplitude free vibration analysis of moderately thick FG rectangular beams under various boundary conditions was developed in this section. The symbols L, h, and b denote the FG beam’s length, thickness, and width. In addition, a set of Cartesian coordinates (x, y, z) for the kinematics description of the FG beam was located at the mid-surface of its left end.

The structural behavior of the FG beam was described using a mixed TBT [39,40,41,42], in which the shear deformation effects were considered to be a constant through the thickness direction of the FG beam, and the related displacement field was given, as follows:

(1)$u_x\left(x, z\right)=u\left(x\right)-z\varphi \left(x\right),$

(2)$u_y\left(x, z\right)=0,$

(3)$u_z\left(x, z\right)=w\left(x\right),$

where $u_i\left(x, z\right)$ (i = x, y, and z) denote the displacement components of the FG beam in the x, y, and z directions, respectively. $u\left(x\right)$ and $w\left(x\right)$ stand for the mid-plane displacement components of the FG beam in the x and z directions and $\varphi \left(x\right)$ is the total rotation in the x–z plane.

The strain-displacement relations of the FG beam considering the VKGN effect were given by:

(4)${\epsilon }_x=u{,}_x-z \varphi {,}_x+\left(1/2\right){\left(w{,}_x\right)}^2,$

(5)${\gamma }_{xz}=-\varphi +w{,}_x,$

(6)${\epsilon }_y={\epsilon }_z={\gamma }_{yz}={\gamma }_{xy}=0,$

where ${\epsilon }_x, {\epsilon }_y, {\epsilon }_z, {\gamma }_{xz}, {\gamma }_{yz} \mathrm{and} {\gamma }_{xy}$ are the strain components of the FG beam, and $g{,}_x=\partial g/\partial x$, in which $g=u, w, \mathrm{and} \varphi$.

The nonzero stress components of the FG beam were given by:

(7)${\sigma }_x=E\left[u{,}_x-z \varphi {,}_x+\left(1/2\right){\left(w{,}_x\right)}^2\right],$

(8)${\tau }_{xz}= k_c G\left(-\varphi +w{,}_x\right),$

where $G$ and $E$ are defined as the shear and Young’s moduli of the FG beam, respectively, and in general, these are specific functions of z, i.e., $G=G\left(z\right)$ and $E=E\left(z\right)$. $k_c$ denotes the shear stress correction factor of the FG beam, which was taken as 5/6 in this work.

The generalized force resultants, namely the axial force $N$, moment $M$, and shear force $Q$, of the FG beam, were defined by:

(9)$N=A_{xx} u{,}_x-B_{xx} \varphi {,}_x+(A_{xx}/2){\left(w{,}_x\right)}^2,$

(10)$M=-B_{xx} u{,}_x+D_{xx} \varphi {,}_x-(B_{xx}/2){\left(w{,}_x\right)}^2,$

(11)$Q=-k_c A_{xz}\left(-\varphi +w{,}_x\right),$

where for an $N_l-\mathrm{layered}$ FG beam:

$A_{xx}={\displaystyle \sum _{m=1}^{N_l} }b {\displaystyle {\int }_{z_{m-1}}^{z_m} E\left(z\right) dz}, B_{xx}={\displaystyle \sum _{m=1}^{N_l} }b {\displaystyle {\int }_{z_{m-1}}^{z_m} E\left(z\right) z dz}, D_{xx}={\displaystyle \sum _{m=1}^{N_l} }b {\displaystyle {\int }_{z_{m-1}}^{z_m} E\left(z\right) z^2 dz}, A_{xz}={\displaystyle \sum _{m=1}^{N_l} }b {\displaystyle {\int }_{z_{m-1}}^{z_m} G\left(z\right) dz},$

where $N_l$ denotes the total number of layers constituting the multi-layered beam, and z_m and z_{m − 1} are the thickness coordinates of the top and bottom surfaces of the m^th-layer when measured from the mid-surface of the FG beam.

Four different boundary conditions of the FG beam were considered as follows:

Case 1. S-S supported:

(12)$u=w=M=0 \mathrm{at} x=0 \mathrm{and} x=L;$

Case 2. C-S supported:

(13)$\begin{array}{cc}\hfill u=w=\varphi =0 & \mathrm{at} x=0,\hfill \\ \hfill \mathrm{and} u=w=M=0 & \mathrm{at} x=\mathrm{L};\hfill \end{array}$

Case 3. C-C supported:

(14)$u=w=\varphi =0 \mathrm{at} x=0 \mathrm{and} x=L;$

Case 4. C-F supported:

(15)$\begin{array}{cc}\hfill u=w=\varphi =0& \mathrm{at} x=0,\hfill \\ \hfill \mathrm{and} N=M=Q=0& \mathrm{at} x=L.\hfill \end{array}$

Based on Hamilton’s principle, a weak-form formulation was derived using a variational approach, in which the RMVT combined with the TBT and VKGN kinematics was used. The energy functional of the FG beam could be written in the form of:

(16)$I={\displaystyle {\int }_{t_1}^{t_2}\left(T-{\Pi }_R\right)} dt,$

where T and ${\Pi }_R$ represent the kinetic energy and Reissner’s potential energy of the FG beam, respectively, and can be given by:

(17)$\begin{array}{l}T={\displaystyle \sum _{e=1}^{N_e}}{\displaystyle {\int }_{x_e}^{x_{e+1}}\left[\left(1/2\right)m_0 {\left(\partial u^{(e)}/\partial t\right)}^2+\left(1/2\right)m_2{\left(\partial {\varphi }^{(e)}/\partial t\right)}^2\right. }\\ \left.+\left(1/2\right)m_0{\left(\partial w^{(e)}/\partial t\right)}^2-m_1\left(\partial u^{(e)}/\partial t\right)\left(\partial {\varphi }^{(e)}/\partial t\right)\right]dx,\end{array}$

(18)$\begin{array}{c}{\Pi }_R={\displaystyle \sum _{e=1}^{N_e}}{\displaystyle {\int }_{x_e}^{x_{e+1}}\left\{ N^{(e)}\left[u^{(e)}{,}_x+\left(1/2\right){\left(w^{(e)}{,}_x\right)}^2\right]+M^{(e)} \left({\varphi }^{(e)}{,}_x\right)\right.}+Q^{(e)}\left({\varphi }^{(e)}-w^{(e)}{,}_x\right)-\left(S_{Axx}/2\right){\left(N^{(e)}\right)}^2-\left(S_{Dxx}/2\right){\left(M^{(e)}\right)}^2\\ \left.-S_{Bxx}\left(N^{(e)}\right)\left(M^{(e)}\right)-\left[1/\left(2k_c{A}_{xz}\right)\right]{\left(Q^{(e)}\right)}^2\right\}dx+{\overline{N}}_0 u_0-{\overline{N}}_L u_L+{\overline{M}}_0{\varphi }_0-{\overline{M}}_L{\varphi }_L-{\overline{V}}_0{w}_0+{\overline{V}}_L{w}_L,\end{array}$

where the superscript e denotes a typical beam element, N_e denotes the total number of beam elements constituting the FG beam; $m_0$, $m_1$, and $m_2$ are the mass of the cross-sectional area and its moment of inertia, and are given by $m_0={\displaystyle {\int }_A \rho dA}$, $m_1={\displaystyle {\int }_A \rho z dA}$, and $m_2={\displaystyle {\int }_A \rho z^2 dA}$; ${\delta }_{kl}$ is the Kronecker delta symbol, in which ${\delta }_{kl}$ = 1 when $k=l$, while ${\delta }_{kl}$ = 0 when $k\ne l$; $S_{Axx}=D_{xx}/\left(A_{xx}{D}_{xx}-B_{xx}{B}_{xx}\right)$, $S_{Bxx}=B_{xx}/\left(A_{xx}{D}_{xx}-B_{xx}{B}_{xx}\right)$, and $S_{Dxx}=A_{xx}/\left(A_{xx}{D}_{xx}-B_{xx}{B}_{xx}\right)$; when the material properties constituting the FG beam were symmetric with respect to the mid-surface, $B_{xx}=0$; otherwise, $B_{xx}\ne 0$. In addition, the variables, ${\overline{N}}_0,$ ${\overline{N}}_L,$ ${\overline{M}}_0, {\overline{M}}_L,$ ${\overline{V}}_0 , \mathrm{and} {\overline{V}}_L$, represented the applied axial forces, moments, and shear forces of the FG beam at the edges.

Based on the RMVT-based TBT, the generalized displacement variables $u^{(e)}$, $w^{(e)}, \mathrm{and} {\varphi }^{(e)}$ and the generalized force resultant variables $N^{(e)}, M^{(e)}, \mathrm{and} Q^{(e)}$ of each element were selected as the primary variables subject to variation, and the time function of each field variable was assumed to be the harmonic function, i.e., $F^{(e)}\left(x, z, t\right)={\tilde{F}}^{(e)}\left(x, z\right) e^{i \omega t}$, in which t represents the time variable, and $\omega$ is the natural frequency of the FG beam. Then, by applying Hamilton’s principle and imposing the continuity conditions at the nodes of the beam element, the authors finally obtained the motion equations of the loaded FG beam, as follows:

(19)${\displaystyle \sum _{e=1}^{N_e} }\left({\left[\begin{array}{cccccc}0& 0& 0& k_{14}^{(e)}& 0& 0\\ 0& 0& 0& 0& 0& k_{26}^{(e)}\\ 0& 0& 0& 0& k_{35}^{(e)}& k_{36}^{(e)}\\ k_{41}^{(e)}& 0& 0& k_{44}^{(e)}& k_{45}^{(e)}& 0\\ 0& 0& 0& k_{54}^{(e)}& k_{55}^{(e)}& 0\\ 0& k_{62}^{(e)}& k_{63}^{(e)}& 0& 0& k_{66}^{(e)}\end{array}\right]}^{(m)} +{\left[\begin{array}{cccccc}0& 0& 0& 0& 0& 0\\ 0& g_{22}^{(e)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& g_{42}^{(e)}& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]}^{(m-1)}\left.-{\omega }^2 {\left[\begin{array}{cccccc}{m}_{11}^{(e)}& 0& m_{13}^{(e)}& 0& 0& 0\\ 0& m_{22}^{(e)}& 0& 0& 0& 0\\ m_{31}^{(e)}& 0& m_{33}^{(e)}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]}^{(m)}\right){\left\{\begin{array}{c}{u}_j^{(e)}\\ w_j^{(e)}\\ {\varphi }_j^{(e)}\\ N_j^{(e)}\\ M_j^{(e)}\\ Q_j^{(e)}\end{array}\right\}}_{ i=j=1,\dots , n_d}^{(m)}={\displaystyle \sum _{e=1}^{N_e} }{\left({\left\{\begin{array}{c}0\\ 0\\ 0\\ 0\\ 0\\ 0\end{array}\right\}}^{(m)}\right)}_{i=1,\dots , n_d},\right.$

where $k_{ ij}^{(e)}, g_{ ij}^{(e)} \mathrm{and} m_{ ij}^{(e)}$ are relevant coefficients and the superscripts m and (m − 1) refer to the m^th and (m − 1)^th iterative processes.

Based on the linear vibration theory without considering the VKGN, the linear natural frequencies of the FE beam, which were independent of the modal variables, could be obtained by letting the geometrical stiffness coefficients be zero (i.e., $g_{ ij}^{(e)}$ = 0). Then, the corresponding eigenvectors were scaled up by allowing the maximum transverse displacement of the FG beam to be identical to a given vibration amplitude ${\overline{w}}_{\max }$, and $g_{ ij}^{(e)}$ could be determined using these scaled-up eigenvectors. The nonlinear natural frequencies of the FG beam, which were dependent upon the modal variables, could thus be obtained using the updated eigenvalue system equations (Equation (19)). The iteration process ended when the relative error between the nonlinear frequencies of m^th and (m − 1)^th iterations were less than 10⁻⁵.

Using Equation (19) and the above-mentioned iterative process, the authors could investigate the large amplitude free vibration characteristics of moderately thick FG beams under different boundary conditions.

3. The MP BPNN

3.1. Feeding Forward Process

An MP BPNN was developed and briefly described in this section. For illustrative purposes, a schematic for the structure of the N_h- layer network is shown in Figure 1, in which N_h = 2; $p_j \left(j=1,\dots , R\right)$ denotes the R inputs; $w_{ij}^m$ is the weight representing the connection to the i^th-neuron of the m^th-layer from the j^th-neuron of the (m − 1)^th-layer, such that m = 1,…, N_h, i = 1,…, S^m and j = 1,…, S^{(m − 1)}, as well as S^m, which is the total number of neurons used in the m^th-layer; $b_j^m$ is the bias of the j^th-neuro of the m^th-layer; $n_i^m$ denotes the net input of the i^th-neuron of the m^th-layer; $n_i^m={\displaystyle \sum _{j=1}^{S^{(m-1)}}{w}_{ij}^m}{a}_j^{(m-1)}+b_j^m$; $a_i^m\left(i=1,\dots , S^m\right)$ is the i^th-output of the m^th-layer, and $a_i^m =f^m\left(n_i^m\right)$, where $f^m$ is the transfer function used in the m^th-layer, and $a_j^0 =p_j$. A symbol of R-S¹-S² in this work was thus used to represent a two-layer network algorithm with R inputs, S¹ and S² neurons in the first and second layers, and S² outputs, in which the number of neurons used for the last layer was the same as the number of final outputs.

For a multilayer network, the output of one layer became the input to the following layer, the operation of which as thus given as:

(20)$a^{m+1}=f^{m+1}\left(n^{m+1}\right)=f^{m+1}\left(W^{m+1}{a}^m+b^{m+1}\right) \mathrm{for} m=N_h,\dots , 1, 0,$

where $b^{m+1}={\left\{\begin{array}{cccc}{b}_{ 1}^{m+1}& b_{ 2}^{m+1}& \cdots & b_{S^{m+1}}^{m+1}\end{array}\right\}}^T$; $n^{m+1}={\left\{\begin{array}{cccc}{n}_{ 1}^{m+1}& n_{ 2}^{m+1}& \cdots & n_{S^{m+1}}^{m+1}\end{array}\right\}}^T$. In this work, the transfer function for the hidden layers of m = 1,…, (N_h − 1) was selected as the log-sigmoid function, while the last output layer of m = N_h was selected as the linear transfer function, such that $f^{m+1}={\left\{\begin{array}{cccc}{f}_{ 1}^{m+1}& f_{ 2}^{m+1}& \cdots & f_{S^{m+1}}^{m+1}\end{array}\right\}}^T$, in which $f_{ i}^{m+1}=1/\left(1+e^{-n_{ i}^{m+1}}\right)$ and $f_i^{N_h}=n_i^{N_h}$.

In the first hidden layer, the neurons received external inputs as follows:

(21)$a^0=p,$

where $p=\left\{\begin{array}{cccc}{p}_1& p_2& \cdots & p_R\end{array}\right\}$, which provided the starting point for Equation (21).

The outputs of the neurons in the last layer were considered as the network outputs and defined as:

(22)$a=a^{N_h}.$

3.2. Backpropagation Process

In order to operate the algorithm, there were appropriate input data sets and their corresponding target output sets, which were provided and given as:

(23)$\left(p_1, t_1\right), \left(p_2, t_2\right), \dots \left(p_q, t_q\right), \dots , \left(p_Q, t_Q\right).$

The learning rule used in this ANN algorithm is called the least mean square algorithm, which means that the network parameters could be adjusted to minimize the mean square error as follows:

(24)$\begin{array}{ll}\mathrm{Min} \widehat{F}\left(x\right)& ={\displaystyle \sum _{q=1}^Q}{\left(t_q-a_q\right)}^T\left(t_q-a_q\right)\\ & ={\displaystyle \sum _{q=1}^Q{\left(e_q\right)}^T}{e}_q={\displaystyle \sum _{q=1}^Q{\displaystyle \sum _{j=1}^{S^{N_h}}{\left(e_{jq}\right)}^2}}\\ & ={\displaystyle \sum _{i=1}^{n_e}{\left(v_i\right)}^2},\end{array}$

where $x=\left\{w_{11}^1,\cdots , w_{S^1 R}^1, b_1^1,\cdots , b_{S^1}^1,\cdots , w_{11}^{N_h},\cdots , w_{S^{N_h} S^{N_h-1}}^{N_h},b_{ 1}^{N_h},\cdots , b_{S^{N_h}}^{N_h} \right\}$, in which the total number (n_p) for the parameters, including all the weights and biases in this MP BPNN algorithm, is $n_p=\left(R+1\right)S^1+\left(S^1+1\right)S^2+\cdots \left(S^{N_h-1}+1\right)S^{N_h}$. $v^T=\left\{v_1 v_2 \cdots v_{n_e} \right\}=\left\{e_{11} e_{21} \cdots e_{S^{N_h} 1} e_{12} e_{22} \cdots e_{S^{N_h} 2} \cdots e_{S^{N_h} Q}\right\}$, in which the total number (n_e) of the square error terms is $n_e=Q S^{N_h}$.

In order to speed up the convergence rate, the Levenberg–Marquardt algorithm [26] was adopted for network training.

When the input $P_q$ was fed into the network and the corresponding network output $a_q^{N_h}$ was computed, the Levenberg–Marquardt BP algorithm could be initialized with:

(25)${\tilde{S}}^{ N_h}=\left[\begin{array}{cccc}{\tilde{S}}_{ 1}^{ N_h}& {\tilde{S}}_{ 2}^{ N_h}& \cdots & {\tilde{S}}_{ Q}^{ N_h}\end{array}\right]$

where ${\tilde{S}}_{ q}^{N_h}=-{\dot{F}}^{N_h}\left(n_{ q}^{N_h}\right)$, and the dimension of ${\dot{F}}^m\left(n^m\right)$ is $S^m\mathrm{x}{S}^m$ and that of ${\tilde{S}}^m$ is $S^m\times \left(Q S^m\right)$ for m = 1,…, N_h.

The total Marquardt sensitivity matrices ${\tilde{S}}^{ m}$ could then be obtained in the following order: m = N_h − 1, N_h − 2, …, 1, using the BP network, and this could be written as follows:

(26)${\tilde{S}}^{ m}=\left[\begin{array}{cccc}{\tilde{S}}_{ 1}^{ m}& {\tilde{S}}_{ 2}^{ m}& \cdots & {\tilde{S}}_{ Q}^{ m}\end{array}\right] \mathrm{for} m=N_h-1, N_h-2,\dots , 1,$

where ${\tilde{S}}_q^{ m}={\dot{F}}^{ m}\left(n_{ q}^m\right) {\left(W^{m+1}\right)}^T {\tilde{S}}_{ q}^{ m+1}$ for q = 1,…, Q.

According to the Levenberg–Marquardt BP algorithm, all of the parameters could be modified using an iteration process as follows:

(27)$X^{(k+1)}=X^{(k)}-{\left[J^T\left(x^{(k)}\right)J\left(x^{(k)}\right)+{\mu }^{(k)} I\right]}^{-1} J^T\left(x^{(k)}\right)v\left(x^{(k)}\right),$

where k denotes the number of iterations, and the initial value of $\mu$, i.e., ${\mu }^{(0)}$, was taken to be a significant value (e.g., ${\mu }^{(0)}$ = 10,000). If a step yielded a more significant value for F(x), then the step was repeated with ${\mu }^{(k)}$ multiplied by a factor $\vartheta$ > 1 (e.g., $\vartheta$ = 2 in this work), and if a step yielded a smaller value for F(x), then the step was repeated with ${\mu }^{(k)}$ divided by a factor $\vartheta$ > 1 (e.g., $\vartheta$ = 3 in this work). Eventually, the iteration process converged. In this work, the stop criteria of the iteration process were considered as follows:

• (1). When the root of the mean square relative error R_e was less than 10⁻⁵, i.e.,

  (28)$R_e=\sqrt{{\displaystyle \sum _{q=1}^Q\left[{\left(t_q-a_q\right)}^T\left(t_q-a_q\right)/\left(t_q^T{t}_q\right)\right]}}/(Q S^{N_h})\le {10}^{-5}$

• (2). When the number of iterations was greater than 20,000.

4. Illustrative Examples

4.1. Large Amplitude-Free Vibration Analysis Using the Mixed FE Method

In this section, the authors investigate the large amplitude-free vibration behavior of homogeneous isotropic and FG isotropic beams with combinations of simply supported, free, and clamped boundary conditions using the mixed FE method.

Table 1 shows the results of the convergence study for the mixed FE solutions of the nonlinear-to-linear natural frequency ratios (${\omega }_{nl}/{\omega }_{linear}$) of homogeneous isotropic beams under S-S boundary conditions. This issue has also been studied using various FE methods, such as the Lagrange-type FE method [43] and the Galerkin FE method [44]. The accuracy and convergence rate of the RMVT-based FE methods with different orders used to expand the primary variables in the element domain could thus be validated by comparing their solutions with those available in the literature.

The geometric parameters of the beams were L/h = 20 and h = 0.1 m. The material properties were $E=322.3 \mathrm{GPa}$, and $\upsilon =0.24$. A maximum dimensionless amplitude was defined as obtaining the nonlinear natural frequencies of the beams.

It can be seen in Table 1 that the RMVT-based (or mixed) FE solutions converged rapidly and that convergent solutions were obtained when the 16 linear, eight quadratic or four cubic elements were used. Moreover, the convergent solutions were in excellent agreement with those obtained using various PVD-based FE methods [43,44].

Table 2 compares the mixed FE solutions for the natural frequency ratios of FG beams under C-C boundary conditions with the solutions reported in the literature, in which 16 cubic elements were used. In addition, mixed FE solutions for the S-S, C-S, and C-F boundary conditions were also presented. The FG beam was made of a ceramic-metal two-phase material, comprised Silicon nitride Si3N4 (ceramic) and a stainless steel SuS304 (metal), and the material properties of the Si3N4 and SuS304 materials were $E_c=322.3 \mathrm{GPa}$, ${\upsilon }_c=0.24,$ and ${\rho }_c=2370 \mathrm{K}\mathrm{g}/{\mathrm{m}}^3$, as well as $E_m=207.8 \mathrm{GPa}$, ${\upsilon }_m=0.3178 ,$ and ${\rho }_m=8166 \mathrm{K}\mathrm{g}/{\mathrm{m}}^3$, in which the subscripts c and m denoted the ceramic and metal materials, respectively. The material properties of the FG beam were assumed to obey the power-law distributions through the thickness direction according to the volume fractions of the constituents, which were ${\Gamma }_c\left(z\right)={\left[1/2+\left(z/h\right)\right]}^{{\kappa }_p}$ and ${\Gamma }_m\left(z\right)=1-{\Gamma }_c\left(z\right)$, in which ${\kappa }_p$ represents the material-property gradient index, and the effective material properties of this were estimated using the rule of mixture, as follows:

(29)$P\left(z\right)=P_c{\Gamma }_c\left(z\right)+P_m{\Gamma }_m\left(z\right)$

where P represents E, $\upsilon$, and $\rho$.

It can be seen in Table 2 that the nonlinear-to-linear natural frequency ratios increased when the maximum dimensionless amplitude became greater for C-C, C-S, and S-S boundary conditions. Again, it can be seen in Table 2 that for the C-C boundary conditions, the mixed FE solutions closely agreed with the solutions obtained by Elmaguiri et al. [45] and Ke et al. [46] with Euler-Bernoulli’s beam theory accounting for the VKGN effect. These results also show that the magnitudes of the nonlinear-to-linear natural frequency ratio at increments for different boundary conditions were arranged in descending order: S-S boundary conditions > C-S boundary conditions> C-C boundary conditions.

On the other hand, for C-F boundary conditions, the nonlinear-to-linear natural frequency ratios slightly decreased when the maximum dimensionless amplitude became greater. This was because the geometrical nonlinearity effect and the boundary constraints resulted in an axial tensional force for the S-S, C-S, and C-C boundary conditions, enhancing the overall stiffness of the FG beams. In contrast, this resulted in a very small compressive force for the C-F boundary conditions, which slightly weakened the overall stiffness of the FG beams. The results also showed that the influence of the geometrical nonlinearity effect on the natural frequencies of the FG beams for different boundary conditions was arranged in descending order: S-S > C-S > C-C > C-F boundary conditions.

The modal shapes of the FG beam under different boundary conditions are shown in Figure 2, in which the maximum dimensionless amplitude and the material-property gradient index were taken to be ${\overline{w}}_{\max }=w_{\max }/\sqrt{I/A}$ = 1 and ${\kappa }_p$ = 0.3, respectively.

4.2. Large Amplitude–Free Vibration Analysis Using the Developed MP BPNN

In this section, the large amplitude free vibration analysis of the FG rectangular beams mentioned above was undertaken again using the developed MP BPNN algorithm, which was trained using the mixed FE solutions. The boundary conditions of the FG beams were considered to be clamped-clamped edge conditions. The cross-section of the FG beam was 10 cm $\times$ 10 cm, and their material properties were considered the same as those used in Table 2. The dimensionless nonlinear frequency parameter was defined as ${\overline{\omega }}_{nl}={\omega }_{nl} \sqrt{{\rho }_m L^2/E_m}$.

After training, a set of optimal parameters, including the weight numbers $w_{ij}^m$ and biases $b_j^m$, where i = 1,…, S^m, j = 1,…, S^{(m − 1)}, and m = 1,…, N_h, could be obtained, and the trained MP BPNN could replace the mixed FE method to predict the nonlinear natural frequencies of the FG beams significantly faster compared to the case with the mixed FE method. In the training of this MP BPNN, R was taken as three and represented the length-to-thickness ratio (L/h), the material-property gradient index (${\kappa }_p$), and the maximum dimensionless amplitude ($w_{\max }/ \left(\sqrt{I/A}\right)$), in which L/h = ${10}^{1+0.05i}$, i = 0,…, 20; ${\kappa }_p$ = ${10}^{-2+0.2i}$, i = 0,…, 20; $w_{\max }/ \left(\sqrt{I/A}\right)$ = 1 + i, i = 0,…, 3, such that the total number of the training data sets was 21 × 21 × 4 = 1764. The number of outputs was taken as two, which were the first and second lowest natural frequencies of the FG beams. Thus, there were 1764 training data sets for the first lowest and second lowest nonlinear natural frequencies of the FG beams, which were obtained using the mixed FE methods and used to train the $n_p$ parameters for each MP BPNN, in which $n_p=\left(R+1\right)S^1+\left(S^1+1\right)S^2+\cdots \left(S^{N_h-1}+1\right)S^{N_h}$, the initial values of which were randomly generated between −1 and 1. There were 100 sets of the initial values of the n_p parameters used for each MP BPNN algorithm in this analysis, for which the first three best results are shown in Table 3 and were estimated using the R_e of the outputs compared to the other 200 selected testing data.

In this section, the large amplitude free vibration analysis of FG rectangular beams mentioned above was undertaken again using the developed MP BPNN algorithm, which was trained using the mixed FE solutions. The FG beams’ boundary conditions were considered clamped-clamped edge conditions. The cross-section of the FG beam was 10 cm $\times$ 10 cm, and their material properties were considered the same as those used in Table 2. The dimensionless nonlinear frequency parameter was defined as ${\overline{\omega }}_{nl}={\omega }_{nl} \sqrt{{\rho }_m L^2/E_m}$.

After training, a set of optimal parameters, including the weight numbers $w_{ij}^m$ and biases $b_j^m$, where i = 1,…, S^m, j = 1,…, S^(m − 1), and m = 1,…, N_h, could be obtained, and the trained MP BPNN replaced the mixed FE method to predict the nonlinear natural frequencies of the FG beams significantly faster compared to the case with the mixed FE method. In the training of this MP BPNN, R was taken as three, and represented the length-to-thickness ratio (L/h), the material-property gradient index (${\kappa }_p$), and the maximum dimensionless amplitude ($w_{\max }/ \left(\sqrt{I/A}\right)$), in which L/h = ${10}^{1+0.05i}$, i = 0,…, 20; ${\kappa }_p$ = ${10}^{-2+0.2i}$, i = 0,…, 20; $w_{\max }/ \left(\sqrt{I/A}\right)$ = 1 + i, i = 0,…, 3, such that the total number of the training data sets was 21 × 21 × 4 = 1764. The number of outputs was taken as two, which were the first and second lowest natural frequencies of the FG beams. Thus, there were 1764 training data sets for the first lowest and second lowest nonlinear natural frequencies of the FG beams obtained using the mixed FE methods, and these were used to train the $n_p$ parameters for each MP BPNN, in which $n_p=\left(R+1\right)S^1+\left(S^1+1\right)S^2+\cdots \left(S^{N_h-1}+1\right)S^{N_h}$, the initial values of which were randomly generated between −1 and 1. There were 100 sets of the initial values of the n_p parameters used for each MP BPNN algorithm in this analysis, for which the first three best results are shown in Table 3; these were estimated using the R_e of outputs as compared with another 200 selected testing data.

It can be seen in Table 3 that the 3-6-6-2 and 3-4-4-4-2 BPNN algorithms yielded optimal results, for which the values of R_e were 0.0006% and 0.0008%, respectively. The 3-6-6-2 BPNN algorithm was thus used to predict the first lowest and second lowest nonlinear natural frequencies of the FG beams, as shown in Figure 3 and Figure 4. Furthermore, it can be seen in each curve of Figure 3 and Figure 4 that the solutions obtained using the mixed FE method, drawn using the solid lines and the developed MP BPNN algorithm, and drawn using the discrete symbols closely agreed with each other for a randomly selected group of 50 sets of input data.

Table 4 compares the computer execution time required for the developed MP BPNN and the mixed FE method when randomly selecting 50 input data sets, as shown in Figure 3 and Figure 4. It was shown that the central processing unit (CPU) time required for the mixed FE method was about 2800 times greater than that needed for the developed MP BPNN algorithm. Thus, the trained MP BPNN algorithm was superior to the mixed FE method because it was less time-consuming.

The results also showed that the nonlinear frequency parameters decreased when the length-to-thickness ratio increased. The nonlinear frequency parameters depended on the maximum dimensionless amplitude, and these values increased when the maximum dimensionless amplitude was greater.

5. Conclusions

In this article, based on the RMVT and TBT, a mixed FE method was developed for the nonlinear free vibration analysis of FG beams under the combinations of simply supported, free, and clamped edge conditions. Implementing the mixed FE method shows that its convergent solutions were obtained when 16 linear, eight quadratic, or four cubic elements were used and that these convergent solutions were in excellent agreement with the accuracy solutions reported in the literature. These results also indicated that the nonlinear natural frequency-to-linear natural frequency ratios increased when the maximum dimensionless amplitude increased. On the other hand, for C-F boundary conditions, the nonlinear-to-linear natural frequency ratios decreased when the maximum dimensionless amplitude became greater. The influence of the geometrical nonlinearity effect on the natural frequencies of the FG beams for different boundary conditions was arranged in descending order: S-S > C-S > C-C > C-F boundary conditions.

An MP BPNN algorithm was also developed, in which the Levenberg–Marquardt algorithm was used to speed up MP BPNN’s convergence rate. After training using the mixed FE solutions, it was shown that the trained MP BPNN algorithm was superior to the mixed FE method because it was approximately 1/2800 times less time-consuming. It is worth mentioning that the above comparison did not account for the cost of generating FEM data and training the BPNN. It is expected that the developed MP BPNN algorithm could be extended to research subjects, including inverse problems, optimal design, structural control, etc., when the traditional analytical or numerical methods used in such studies are very time-consuming. Therefore, applying the MP BPNN to the research mentioned above is ongoing.

Author Contributions

Conceptualization, C.-P.W.; methodology, C.-P.W.; software, S.-T.Y. and J.-H.L.; validation, C.-P.W., S.-T.Y. and J.-H.L.; formal analysis, C.-P.W.; investigation, C.-P.W.; resources, C.-P.W.; data curation, S.-T.Y. and J.-H.L.; writing—original draft preparation, C.-P.W.; writing—review and editing, C.-P.W.; visualization, C.-P.W.; supervision, C.-P.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data supporting reported results can be found at https://drive.google.com/drive/folders/1TmbVHAdo6sZNMGytFrmZPrLPL4_Ye4aN?usp=share_link (accessed on 1 May 2023).

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Schematic of the structure of a two-layer network.

Figure 2. The 1st, 2nd, and 3rd mode shapes of the FG beams under (a) S−S, (b) C−S, (c) C−C, and (d) C−F boundary conditions.

Figure 2. The 1st, 2nd, and 3rd mode shapes of the FG beams under (a) S−S, (b) C−S, (c) C−C, and (d) C−F boundary conditions.

Figure 3. Variations in the nonlinear frequency parameters of FG beams with the aspect ratio using the mixed FE method and the developed MP BPNN algorithm for [Forumla omitted. See PDF.] = 1 and 4, and [Forumla omitted. See PDF.] = 4; (a) The first lowest frequency parameters, (b) The second lowest frequency parameter.

Figure 4. Variations in the nonlinear frequency parameters of FG beams with the aspect ratio using the mixed FE method and the proposed MP BPNN algorithm for [Forumla omitted. See PDF.] = 0.1, 1, and 10, and [Forumla omitted. See PDF.] = 2; (a) The first lowest frequency parameters, (b) The second lowest frequency parameters.

Table 1

Results of the convergence study for the mixed FE solutions of the nonlinear-to-linear natural frequency ratios of isotropic beams under S-S boundary conditions.

Theories	${\mathit{w}}_{\mathbf{max}}/\sqrt{\mathit{I}/\mathit{A}}$
	1.0	2.0	3.0	4.0
Current four linear elements	1.1196	1.4193	1.8117	2.2488
Current eight linear elements	1.1187	1.4164	1.8065	2.2412
Current 16 linear elements	1.1187	1.4163	1.8062	2.2407
Current four quadratic elements	1.1173	1.4124	1.8006	2.2342
Current eight quadratic elements	1.1183	1.4155	1.8053	2.2398
Current 16 quadratic elements	1.1186	1.4161	1.8061	2.2406
Current four cubic elements	1.1187	1.4164	1.8064	2.2410
Current eight cubic elements	1.1187	1.4162	1.8062	2.2407
Sarma and Varadan [43]	1.1180	1.4142	1.8028	2.2361
Bhashyam and Prathap [44]	1.1180	1.4141	1.8027	2.2359

Table 2

Comparisons between the mixed FE solutions of the nonlinear-to-linear natural frequency ratios (${\omega }_{nl}/{\omega }_{linear}$) of FG beams under C-C boundary conditions and the solutions reported in the literature, in which the mixed FE solutions for the S-S, C-S, and C-F boundary conditions are also presented.

Boundary Conditions	$${\mathsf{\kappa }}_{\mathit{p}}$$	Theories	$${\overline{\mathit{W}}}_{\mathit{max}}={\mathit{W}}_{\mathit{max}}/\sqrt{\mathit{I}/\mathit{A}}$$
			1	2	3	4
C-C	0.3	Mixed FE method	1.0305	1.1164	1.2445	1.4017
		Elmaguiri et al. [45]	1.0227	1.0875	1.1869	1.3121
		Ke et al. [46]	1.0220	1.0852	1.1831	1.3079
	1	Mixed FE method	1.0304	1.1158	1.2434	1.4000
		Elmaguiri et al. [45]	1.0225	1.0871	1.1860	1.3106
		Ke et al. [46]	1.0025	1.0873	1.1874	1.3149
	2	Mixed FE method	1.0296	1.1128	1.2374	1.3906
		Elmaguiri et al. [45]	1.0215	1.0832	1.1780	1.2980
		Ke et al. [46]	1.0232	1.0900	1.1929	1.3237
S-S	0.3	Mixed FE method	1.1492	1.4666	1.8681	2.3102
	1	Mixed FE method	1.1704	1.4964	1.8993	2.3397
	2	Mixed FE method	1.1695	1.4899	1.8855	2.3182
C-S	0.3	Mixed FE method	1.0696	1.2379	1.4667	1.7299
	1	Mixed FE method	1.0764	1.2481	1.4779	1.7409
	2	Mixed FE method	1.0755	1.2438	1.4691	1.7272
C-F	0.3	Mixed FE method	1.00004	0.99998	0.99979	0.99948
	1	Mixed FE method	1.00008	1.00004	0.99989	0.99961
	2	Mixed FE method	1.00006	1.00002	0.99985	0.99956

Table 3

General information regarding the developed MP BPNNs and their performance.

No. Neuronsof Each Hidden Layer	No. Hidden Layers (i.e., $({\mathit{N}}_{\mathit{h}}-\mathbf{1})$)
	1			2			3			4
	No.Parameters	TrainingTime (s)	Re ofOutputs	No. Parameters	TrainingTime (s)	Re ofOutputs	No. Parameters	TrainingTime (s)	Re ofOutputs	No. Parameters	TrainingTime (s)	Re ofOutputs
2	14	40.99	0.2328%	20	116.50	0.1913%	26	83.46	0.0341%	32	314.25	0.0296%
		31.98	0.2331%		170.48	0.1917%		52.81	0.1315%		525.97	0.0569%
		32.50	0.2331%		103.02	0.2206%		56.06	0.1426%		107.98	0.0572%
4	26	45.13	0.1393%	46	70.94	0.0048%	66	435.87	0.0008%
		25.45	0.1421%		38.68	0.0167%		188.26	0.0043%
		35.50	0.2676%		124.84	0.0168%		78.46	0.0128%
6	38	163.32	0.0237%	80	173.84	0.0006%
		55.70	0.0338%		252.95	0.0007%
		55.11	0.0467%		533.02	0.0008%
8	50	82.67	0.0228%
		86.69	0.0266%
		122.42	0.0359%

Table 4

A comparison of the CPU time required for the proposed MP BPNN and the mixed FEM when randomly selecting a group of 50 input data sets is shown in Figure 2 and Figure 3.

$$[{\mathsf{\kappa }}_{\mathit{p}}, {\mathit{w}}_{\mathbf{max}}/\sqrt{\mathit{I}/\mathit{A}}]$$	Time (s)
	ANN	FEM
[4, 1]	0.0027	7.1167
[4, 4]	0.0025	7.5926
[0.1, 2]	0.0026	7.4780
[1, 2]	0.0026	7.0521
CPU: Central processing unit.

References

1. Koizumi, M. Recent progress of functionally graded materials in Japan. Ceram. Eng. Sci. Proc.; 1992; 13, pp. 333-347.

2. Koizumi, M. FGM activities in Japan. Compos. Part B; 1997; 28, pp. 1-4. [DOI: https://dx.doi.org/10.1016/S1359-8368(96)00016-9]

3. Jha, D.K.; Kant, T.; Singh, R.K. A critical review of recent research on functionally graded plates. Compos. Struct.; 2013; 96, pp. 833-849. [DOI: https://dx.doi.org/10.1016/j.compstruct.2012.09.001]

4. Shen, H.S. Functionally Graded Materials: Nonlinear Analysis of Plates and Shells; CRC Press: Boca Raton, FL, USA, 2009.

5. Chakraborty, A.; Mahapatra, D.R.; Gopalakrishnan, S. Finite element analysis of free vibration and wave propagation in asymmetric composite beams with structural discontinuities. Compos. Struct.; 2002; 55, pp. 23-36. [DOI: https://dx.doi.org/10.1016/S0263-8223(01)00130-1]

6. Chakraborty, A.; Gopalakrishnan, S.; Reddy, J.N. A new beam finite element for analyzing functionally graded materials. Int. J. Mech. Sci.; 2003; 45, pp. 519-539. [DOI: https://dx.doi.org/10.1016/S0020-7403(03)00058-4]

7. Aydogdu, M.; Taskin, V. Free vibration analysis of functionally graded beams with simply supported edges. Mater. Des.; 2007; 28, pp. 1651-1656. [DOI: https://dx.doi.org/10.1016/j.matdes.2006.02.007]

8. Simsek, M. Static analysis of a functionally graded beam under a uniformly distributed load by Ritz method. Int. J. Eng. Appl. Sci.; 2009; 1, pp. 1-11.

9. Simsek, M. Fundamental frequency analysis of functionally graded beams using different higher-order beam theories. Nucl. Eng. Des.; 2010; 240, pp. 697-705. [DOI: https://dx.doi.org/10.1016/j.nucengdes.2009.12.013]

10. Ebrahimi, F.; Farazmandnia, N. Thermo-mechanical vibration analysis of sandwich beams with functionally graded carbon nanotube-reinforced composite face sheets based on a higher-order shear deformation beam theory. Mech. Adv. Mater. Struct.; 2017; 24, pp. 820-829. [DOI: https://dx.doi.org/10.1080/15376494.2016.1196786]

11. Carrera, E. Theories and finite elements for multilayered plates and shells: A unified compact formulation with numerical assessment and benchmarking. Arch. Comput. Methods Eng.; 2003; 10, pp. 215-296. [DOI: https://dx.doi.org/10.1007/BF02736224]

12. De Pietro, G.; Hui, Y.; Giunta, G.; Belouettar, S.; Carrera, E.; Hu, H. Hierarchical one-dimensional finite elements for the thermal stress analysis of three-dimensional functionally graded beams. Compos. Struct.; 2016; 153, pp. 514-528. [DOI: https://dx.doi.org/10.1016/j.compstruct.2016.06.012]

13. Coskun, S.; Kim, J.; Toutanji, H. Bending, free vibration, and buckling analysis of functionally graded porous micro-plates using a general third-order plate theory. J. Compos. Sci.; 2019; 3, 15. [DOI: https://dx.doi.org/10.3390/jcs3010015]

14. Tham, V.V.; Quoc, T.H.; Tu, T.M. Free vibration analysis of laminated functionally graded carbon nanotube-reinforced composite doubly curved shallow shell panels using a new four-variable refined theory. J. Compos. Sci.; 2019; 3, 104. [DOI: https://dx.doi.org/10.3390/jcs3040104]

15. Wu, C.P.; Xu, Z.R. Strong and weak formulations of a mixed higher-order shear deformation theory for the static analysis of functionally graded beams under thermo-mechanical loads. J. Compos. Sci.; 2020; 4, 158. [DOI: https://dx.doi.org/10.3390/jcs4040158]

16. Wu, C.P.; Li, K.W. Multi-objective optimization of functionally graded beams using a genetic algorithm with non-dominated sorting. J. Compos. Sci.; 2021; 5, 92. [DOI: https://dx.doi.org/10.3390/jcs5040092]

17. Ma, L.S.; Lee, D.W. Exact solutions for nonlinear static responses of a shear deformable FGM beam under an in-plane thermal loading. Eur. J. Mech. A/Solids; 2012; 31, pp. 13-20. [DOI: https://dx.doi.org/10.1016/j.euromechsol.2011.06.016]

18. Zhang, D.G. Nonlinear bending analysis of FGM beams based on physical neutral surface and high order shear deformation theory. Compos. Struct.; 2013; 100, pp. 121-126. [DOI: https://dx.doi.org/10.1016/j.compstruct.2012.12.024]

19. Salami, S.J. Extended high order sandwich panel theory for bending analysis of sandwich beams with carbon nanotube reinforced face sheets. Phys. E; 2016; 76, pp. 187-197. [DOI: https://dx.doi.org/10.1016/j.physe.2015.10.015]

20. Shen, H.S. Nonlinear bending of functionally graded carbon nanotube-reinforced composite plates in thermal environments. Compos. Struct.; 2009; 91, pp. 9-19. [DOI: https://dx.doi.org/10.1016/j.compstruct.2009.04.026]

21. Ghayesh, M.H. Nonlinear vibrations of axially functionally graded Timoshenko tapered beams. J. Comput. Nonlinear Dyn.; 2018; 13, 041002. [DOI: https://dx.doi.org/10.1115/1.4039191]

22. Shen, H.S.; Lin, F.; Xiang, Y. Nonlinear vibration of functionally graded graphene-reinforced composite laminated beams resting on elastic foundations in thermal environments. Nonlinear Dyn.; 2017; 90, pp. 899-914. [DOI: https://dx.doi.org/10.1007/s11071-017-3701-0]

23. Ding, J.; Chu, L.; Xin, L.; Dui, G. Nonlinear vibration analysis of functionally graded beams considering the influences of the rotary inertia of the cross section and neural surface position. Mech. Based Des. Struct. Mach.; 2018; 46, pp. 225-237. [DOI: https://dx.doi.org/10.1080/15397734.2017.1329020]

24. Eltaher, M.A.; Mohamed, N.; Mohamed, S.A.; Seddek, L.F. Periodic and nonperiodic modes of postbuckling and nonlinear vibration of beams attached to nonlinear foundations. Appl. Math. Model.; 2019; 75, pp. 414-445. [DOI: https://dx.doi.org/10.1016/j.apm.2019.05.026]

25. Chaudhari, V.K.; Lal, A. Nonlinear free vibration analysis of elastically supported nanotube-reinforced composite beam in thermal environment. Procedia Eng.; 2016; 144, pp. 928-935. [DOI: https://dx.doi.org/10.1016/j.proeng.2016.05.119]

26. Mirzaei, M.; Kiani, Y. Nonlinear free vibration of temperature-dependent sandwich beams with carbon nanotube-reinforced face sheets. Acta Mech.; 2016; 227, pp. 1869-1884. [DOI: https://dx.doi.org/10.1007/s00707-016-1593-6]

27. Babaei, H.; Kiani, Y.; Eslami, M.R. Large amplitude free vibrations of FGM beams on nonlinear elastic foundation in thermal field based on neutral/mid-plane formulations. Iran. J. Sci. Technol. Trans. Mech. Eng.; 2021; 45, pp. 611-630. [DOI: https://dx.doi.org/10.1007/s40997-020-00389-y]

28. Hagan, M.T.; Demuth, H.B.; Beale, M. Neural Network Design; PWS Pub. Comp.: New York, NY, USA, 1995.

29. Yagawa, G.; Okuda, H. Neural networks in computational mechanics. Arch. Comput. Methods Eng.; 1996; 3, pp. 435-512. [DOI: https://dx.doi.org/10.1007/BF02818935]

30. Zhang, Z.; Friedrich, K. Artificial neural networks applied to polymer composites: A review. Compos. Sci. Technol.; 2003; 63, pp. 2029-2044. [DOI: https://dx.doi.org/10.1016/S0266-3538(03)00106-4]

31. Ootao, Y.; Kawamura, R.; Tanigawa, Y.; Nakamura, T. Neural network optimization of material composition of a functionally graded material plate at arbitrary temperature range and temperature rise. Arch. Appl. Mech.; 1998; 68, pp. 662-676. [DOI: https://dx.doi.org/10.1007/s004190050195]

32. Ootao, Y.; Kawamura, R.; Tanigawa, Y.; Imamura, R. Optimization of material composition of nonhomogeneous hollow circular cylinder for thermal stress relaxation making use of neural network. J. Therm. Stress.; 1999; 22, pp. 1-22.

33. Chakraverty, S.; Singh, V.P.; Sharma, R.K. Regression based weight generation algorithm in neural network for estimation of frequencies of vibrating plates. Comput. Methods Appl. Mech. Eng.; 2006; 195, pp. 4194-4202. [DOI: https://dx.doi.org/10.1016/j.cma.2005.08.008]

34. Reddy, M.R.S.; Reddy, B.S.; Reddy, V.N.; Sreenivasulu, S. Prediction of natural frequency of laminated composite plates using artificial neural networks. Engineering; 2012; 4, pp. 329-337. [DOI: https://dx.doi.org/10.4236/eng.2012.46043]

35. Jodaei, A.; Jalal, M.; Yas, M.H. Free vibration analysis of functionally graded annular plates by state-space based differential quadrature method and comparative modeling by ANN. Compos. Part B; 2012; 43, pp. 340-353. [DOI: https://dx.doi.org/10.1016/j.compositesb.2011.08.052]

36. Fetene, B.N.; Shufen, R.; Dixit, U.S. FEM-based neural network modeling of laser-assisted bending. Neural Comput. Appl.; 2018; 29, pp. 69-82. [DOI: https://dx.doi.org/10.1007/s00521-016-2544-9]

37. Subramani, T.; Sharmila, S. Prediction of deflection and stresses of laminated composite plate with artificial network aid. Int. J. Modern Eng. Res.; 2014; 4, pp. 51-58.

38. Liu, X.; Du, X.; Niu, X. A neural network method applied in prediction eigenvalue buckling for sandwich plates. Inform. Technol. J.; 2013; 12, pp. 8129-8134.

39. Reissner, E. On a certain mixed variational theorem and a proposed application. Int. J. Numer. Methods Eng.; 1984; 20, pp. 1366-1368. [DOI: https://dx.doi.org/10.1002/nme.1620200714]

40. Carrera, E. An assessment of mixed and classical theories on global and local response of multilayered orthotropic plates. Compos. Struct.; 2000; 50, pp. 183-198. [DOI: https://dx.doi.org/10.1016/S0263-8223(00)00099-4]

41. Wu, C.P.; Lai, W.W. Free vibration of an embedded single-walled carbon nanotube with various boundary conditions using the RMVT-based nonlocal Timoshenko beam theory and DQ method. Phys. E; 2015; 68, pp. 8-21. [DOI: https://dx.doi.org/10.1016/j.physe.2014.12.004]

42. Wu, C.P.; Lai, W.W. Reissner’s mixed variational theorem-based nonlocal Timoshenko beam theory for a single-walled carbon nanotube embedded in an elastic medium and with various boundary conditions. Compos. Struct.; 2015; 122, pp. 390-404. [DOI: https://dx.doi.org/10.1016/j.compstruct.2014.11.073]

43. Sarma, B.S.; Varadan, T.K. Lagrange-type formulation for finite element analysis of non-linear beam vibrations. J. Sound Vib.; 1983; 86, pp. 61-70. [DOI: https://dx.doi.org/10.1016/0022-460X(83)90943-4]

44. Bhashyam, G.R.; Prathap, G. Galerkin finite element method for nonlinear beam vibrations. J. Sound Vib.; 1980; 72, pp. 191-203. [DOI: https://dx.doi.org/10.1016/0022-460X(80)90652-5]

45. Elmaguiri, M.N.; Haterbouch, M.; Bouayad, A.; Oussouaddi, O. Geometrically nonlinear free vibration of functionally graded beams. J. Mater. Environ. Sci.; 2015; 6, pp. 3620-3633.

46. Ke, L.L.; Yang, J.; Kitipornchai, S. An analytical study on the nonlinear vibration of functionally graded beams. Meccanica; 2010; 45, pp. 743-752. [DOI: https://dx.doi.org/10.1007/s11012-009-9276-1]