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

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


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 sortingbased 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 singlelayered 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 0 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), clampedsimple (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.

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: The strain-displacement relations of the FG beam considering the VKGN effect were given by: x xz w, where are the strain components of the FG beam, and , in which , , and g u w φ = . The nonzero stress components of the FG beam were given by: ( ) 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.,

( )
notes 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: where for an layered − where l N denotes the total number of layers constituting the multi-layered beam, and zm and zm − 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: Case 2. C-S supported: Case 3. C-C supported: Case 4. C-F supported: 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: where T and R Π represent the kinetic energy and Reissner's potential energy of the FG beam, respectively, and can be given by: where the superscript e denotes a typical beam element, Ne denotes the total number of beam elements constituting the FG beam; 0 m , 1 m , and 2 m are the mass of the cross-sectional area and its moment of inertia, and are given by ; when the material properties constituting the FG beam were symmetric with respect to the mid-surface, Q 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., in which t represents the time variable, and ω 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:   ( ) e ij g 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 −5 .
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.

Feeding Forward Process
An MP BPNN was developed and briefly described in this section. For illustrative purposes, a schematic for the structure of the Nhlayer network is shown in Figure 1 A symbol of R-S 1 -S 2 in this work was thus used to represent a twolayer network algorithm with R inputs, S 1 and S 2 neurons in the first and second layers, and S 2 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:  In the first hidden layer, the neurons received external inputs as follows: , 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:

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: 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:  1  1  1  1  11  1  11  1 , , , In order to speed up the convergence rate, the Levenberg-Marquardt algorithm [26] was adopted for network training.
When the input q P was fed into the network and the corresponding network output h N q a was computed, the Levenberg-Marquardt BP algorithm could be initialized with: The total Marquardt sensitivity matrices m S  could then be obtained in the following order: m = Nh − 1, Nh − 2, …, 1, using the BP network, and this could be written as follows: where ( )( ) According to the Levenberg-Marquardt BP algorithm, all of the parameters could be modified using an iteration process as follows: where k denotes the number of iterations, and the initial value of μ , i.e., (0) μ , was taken to be a significant value (e.g., (0) μ = 10,000). If a step yielded a more significant value for F(x), then the step was repeated with ( ) k μ multiplied by a factor ϑ > 1 (e.g., ϑ = 2 in this work), and if a step yielded a smaller value for F(x), then the step was repeated with ( ) k μ divided by a factor ϑ > 1 (e.g., ϑ = 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 Re was less than 10 −5 , i.e., ( ) ( ) ( )  t a t a t t (28) (2) When the number of iterations was greater than 20,000.

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 ( linear nl ω ω / ) 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. . 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 , as well as , 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 powerlaw distributions through the thickness direction according to the volume fractions of the constituents, which were , in which p κ represents the material-property gradient index, and the effective material properties of this were estimated using the rule of mixture, as follows: where P represents E, υ , and ρ .
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 max

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 × 10 cm, and their material properties were considered the same as those used in Table 2. The dimensionless nonlinear frequency parameter was defined as After training, a set of optimal parameters, including the weight numbers , the initial values of which were randomly generated between -1 and 1. There were 100 sets of the initial values of the np 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 Re 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 × 10 cm, and their material properties were considered the same as those used in Table  2. The dimensionless nonlinear frequency parameter was defined as ( ) A I w / / max = 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 p n parameters for each MP BPNN, in which ( ) ( ) ( ) , the initial values of which were randomly generated between −1 and 1. There were 100 sets of the initial values of the np 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 Re 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 Re 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 Figures 3 and 4. Furthermore, it can be seen in each curve of Figures 3 and 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 Figures 3 and 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. 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 Figures 2 and 3. [ p κ , 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.

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.