Free vibration analysis of a functionally graded material beam: evaluation of the Haar wavelet method

The current study focuses on the evaluation of the Haar wavelet method, i.e. its comparison with widely used strong formulation based methods (FDM-finite difference method and DQM-differential quadrature method). A solid element 3D finite element model is developed and the numerical results obtained by using simplified approaches are confirmed.


INTRODUCTION
* Accuracy and complexity are two key factors characterizing any numerical method. The Haar wavelet method (HWM) considered in the current study was introduced by Chen and Hsiao in [1,2] almost 20 years ago and up to now it has been applied for solving a wide class of differential and integral equations covering engineering, economic, etc. problems [3][4][5][6][7]. An overview of the applications of the HWM is given in [8]. The wavelet techniques based on the use of an operational matrix of integration are developed for solving ordinal and partial differential equations in [1][2][3][4][5][6][7][8][9][10] and for integral equations in [11][12][13]. All these studies implement the strong formulation based approach of the HWM. The weak formulation based approach of the HWM was introduced in [14].
Most of the authors characterize the HWM as a simple and effective method [3][4][5][6][7][8][9]. These estimates cover mainly implementation of the HWM, less its * Corresponding author, Maarjus.Kirs@ttu.ee accuracy and convergence results, which are still under development. It is shown in [15] that in the case of function approximation with direct expansion into the Haar wavelet the convergence is of order one. However, according to the HWM approach considered, the highest order derivative included in the differential equation is expanded into a series of Haar functions. Thus, the estimate given in [15] holds good for estimating the accuracy of the highest order derivative, but not the solution of the differential equation. Recently, the convergence theorem of the HWM was proved in [16] for the nth order ordinal differential equations (ODEs) ( 2 n  ). It was stated that the order of convergence of the HWM is equal to two. In [17] the accuracy estimates for the extrapolated results in the case of the fourth order ODE are derived, and it is shown that the order of convergence of the extrapolated results is equal to four (Richardson extrapolation is applied).
The application area of new simple methods often includes problems with advanced material models, constitutive laws, etc., which are not yet (well) covered by commercial software.
A new trend in the development of wavelet methods can be outlined as solution of fractional differential and integral equations [15,[18][19][20][21][22][23][24], which is an area not yet well covered by commercial software (finite element method (FEM), etc.). It is observed in [21] that in the case of fractional ODE the order of convergence of the HWM is equal to two if higher order derivative α in the fractional differential equation exceeds one (α > 1). However, in the case of 0 < α < 1 the order of convergence of the HWM tends to the value 1 + α.
In [25,26] the HWM was adapted for the analysis of structures of functionally graded material (FGM). In the current study the vibration analysis of the FGM beams is performed and the results obtained by the HWM are compared with the corresponding results obtained by using the finite difference method (FDM) and the differential quadrature method (DQM). Selection of FDM and DQM for comparison of results was motivated by the fact that these methods are widely used numerical methods in engineering and are based on strong formulation (the complexity of implementation is similar). The methods considered are implemented by the authors in the MATLAB code.
In order to verify the obtained results and prepare solution procedures for structures with complex geometry and loading cases, the solid element 3D finite element model was developed.

BASICS OF HAAR WAVELETS
The Haar function is defined in [8,9] as Any function ( ) f x that is square integrable and finite in the interval ([A, B] can be expanded into Haar wavelets as The integrals of the Haar functions (1) of order n can be calculated analytically as [9] The integrals of the Haar functions determined by (5) are continuous functions in the interval   , A B .

FREE VIBRATION ANALYSIS OF THE FGM BEAM
In the following the free vibration analysis of the FGM beam is considered [27][28][29]. It is assumed that the material properties of the beam of length L vary axially. The governing differential equation of the beam can be written as The varying properties of the bending stiffness ( ) EI x and the distributed mass per unit length ( ) A x  are described by exponential functions as The reference values of the bending stiffness and distributed mass per unit length at x = 0 are denoted by (0) EI and (0), A  respectively. Relation (7) is used in a number of papers [27][28][29]. The volume fractions of the material corresponding to relation (7) can be derived as 2 2 2 The wavelet method approach considered can be applied for a wide range of functions describing properties of FGM. In the following a more general power law relation for describing FG materials is considered: Here k is the non-negative power-law exponent describing the material variation profile along the length of the beam and the indexes L and R stand for the values of the material properties on the left and right support of the beam, respectively. Relations (9)-(10) seem to be the most widely used relations for describing FGM properties found in the literature [30].
In the following the solution of the partial differential equation (6) is assumed in the form Considering Eqs (7) and (11), the governing differential equation (6) can be rewritten in a non-dimensional form as

.
A EI As a result, the vibration analysis problem of the FGM beam considered above is converted to solving the ordinal differential equation (12). The particular boundary conditions are introduced in Section 7.

THE HAAR WAVELET DISCRETIZATION METHOD
Herein the most commonly used approach of the HWM is employed. According to this method, the highest order derivative existing in a differential equation is expanded into Haar wavelets. Thus, Eq. (12) implies that the fourth order derivative should be expanded into Haar wavelets as (14) where N = 2M is the resolution used.
The solution of the differential governing equation (12) W(X) can be obtained by integrating the expansion (14) four times with respect to X as In (15) the operational matrix of integration (4) P is defined by formulas (5) and T a is a vector of coefficients. The integration constants c 0 ,…,c 3 can be determined for each particular boundary condition separately. Corresponding expressions of the integration constants are omitted for conciseness sake.
Inserting the solution of (15) in the differential equation (12) and assuming uniform grid points in the form one obtains a linear system of algebraic equations, which can be solved with respect to coefficient vector Finally, substituting the values of T a in (15) gives the solution of the posed problem in an analytical form.

CONVERGENCE AND ACCURACY ESTIMATES
The convergence theorem for the HWM is given in [16] for the nth order ODE ( 2 n  ) as THEOREM: Let us assume that is a continuous function on [0,1] and its first derivative is bounded Then the HWM, based on the approach in [1,2], will be convergent, i.e. M E will vanish as the number of collocation points approaches N infinity. The convergence is of the order two The proof of the theorem is given in [16]. Furthermore, the quadrate of the 2 L -norm of the error function can be estimated as 2 2 2 In the case of the considered problem the highest order derivative in differential equation equals four ( 4) n  and formula (19) reduces to Furthermore, it is proved in [17] that in the case of the general fourth-order ODE the accuracy of the results of the HWM can be improved from two to four by applying Richardson's extrapolation method. The theoretical estimates pointed out above are validated numerically in the following section.

FEM SIMULATION MODEL
Commercial analysis software Mechanical APDL 16.0 was used to develop a 3D finite element simulation model for free vibration analysis of an axially functionally graded beam. The FGM beam was partitioned through its length into a number of strips with constant material properties inside the strip (see Fig. 1). Figure 1 shows the mesh of the zoomed right-hand side of the beam corresponding to the third row of Table 1 (5 elements in the thickness and width directions and 500 elements in the length direction). The elements considered were cubical 3D 8-Node Homogeneous Structural Elements SOLID185. The detailed mesh values used are given in column 1 of Table 1.
The geometrical parameters of the beam considered are width (b), height (h), and length (L). The material properties of the steel and aluminium used in the FEM analysis are given in Table 2. The boundary conditions considered correspond to a cantilever beam. The results obtained from FEM analysis were originally in the  dimensional form, i.e. computed for a particular beam with the given geometry, rigidity, and mass per unit length values. In order to compare these results with the results of the FDM, DQM, and HWM, the frequency parameter was converted into the non-dimensional form using the following formula: where f stands for natural/dimensional frequency parameter value (in Hz). The FEM results are discussed in detail in the following section.

NUMERICAL RESULTS
In the following five different boundary conditions of the FGM beam are considered (see Fig. 2) and the results obtained by applying HWM, FDM, and DQM are compared (two symmetric and three non-symmetric conditions). The first two values of the fundamental frequency parameter  are presented in Tables 3 and 4 for a pinnedpinned beam, in Tables 5 and 6 for a clamped-clamped beam, in Tables 7 and 8 for a clamped-pinned beam, in Tables 9 and 10 for a pinned-clamped beam, and in Tables 11 and 12 for a clamped-free beam).
Note that in the FE model all supports with pinned boundary conditions (a, c, and d) have the ability to move in the horizontal direction ( 0 , 0 ) y z x u u u    . In Tables 1-10 the properties of the beam are considered to vary according to formula (7), i.e. by exponential functions.    In Tables 3-10 the value of the parameter  is taken equal to 2. The exact solutions computed based on transcendental algebraic equations derived in [27] are given in the headings of Tables 3-10. Obviously, the convergence of the HWM (also of the FDM and DQM) to the exact solution can be observed in all these tables.
The numerical rates of the convergence, computed for the solutions presented in Table 3, are presented in  Table 13.
In Table 13 the values of N start from 16 because each rate of convergence was computed on the basis of three consecutive values of the solution [16]. The rate of the convergence of the HWM and FDM obviously tends to two, but the DQM has an ultrafast rate for N ≤ 32 and a negative rate for N > 32 (loss of accuracy).
In Table 14, the convergence rates of the extrapolated results of the HWM are given for four different boundary conditions considered above. The Richardson extrapolation method was applied, and it can be seen from Table 14 that the order of the convergence of extrapolated results tends to four in the case of all boundary conditions considered.
Based on results given in Tables 3-12, it can be concluded that in the case of the posed problem the highest accuracy was achieved by applying the DQM, also in most cases the accuracy of the results obtained by the HWM is higher than that obtained by the FDM (there fundamental frequencies in Tables 3 and 7 are exceptions). Detailed analysis of DQM results shows that the maximum accuracy was achieved extremely quickly with N = 16 or N = 32; thereafter the accuracy of the solution decreased with increasing resolution. These results are in agreement with the theoretical concept of the DQM (it is based on the use of high order polynomials whose denominator vanishes for large N) and results found in the literature.
It can be seen from Fig. 3 that in the case of the parameter value k = 1.5 the functions of the elasticity modulus corresponding to the exponential and power law functions (7) and (9) are close (here a steel/aluminium cantilever beam with  = -0.549306 is considered).  In Tables 11 and 12 the FG material properties corresponding to steel/aluminium are considered with steel in the left and aluminium in the right support. The particular values of the material used are presented in Table 2.
The exponential model (7) does not include directly material properties at the right end of the beam. The required values of the material are obtained by determining the value of the parameter  ( = -0.549306).
The first four mode shapes for the above-considered FG steel/aluminium cantilever beam are depicted in Fig. 4. The corresponding mode shapes obtained by a FEM are shown in Figs 5-8.
The results given in Table 15 were obtained by applying the HWM with the general power law function (9)- (10).
The steel/aluminium FGM with properties given in Table 2 is considered and the value of the exponent is taken equal to 1.5 (k = 1.5). The boundary conditions for a clamped-clamped beam are applied.      The FEM results computed for all boundary conditions considered above are given in the last row of each table. The number of elements used is 10 × 10 × 1000.
Dimensionless coordinate x/L First four mode shapes of a FGM Beam Obviously, the values of the frequency parameters computed using 3D FEM analysis are in excellent agreement with those given in Tables 3-12, obtained by applying the HWM, DQM, and FDM.
An example of the results of more detailed FEM analysis is given in Table 15. In this table the first free frequencies for a pinned-pinned beam are presented. The convergence of the solution with increasing mesh can be observed. The results are in agreement with corresponding results obtained by applying the HWM, FDM, and DQM given in Tables 3 and 4.

CONCLUSIONS
Three strong formulation based numerical methods (HWM, FDM, and DQM) were applied for the analysis of the FGM beam and the obtained results were compared. The algorithms for all methods were coded by the authors in MATLAB. Good performance was observed in the case of all three methods used.
It can be concluded that in the case of the considered problem the accuracy of the solutions obtained by applying the HWM and FDM was in the same range. However, in most cases the accuracy of the results of the HWM outperformed that of the FDM. The accuracy of the DQM appears to be higher than that of the HWM and FDM. The convergence results presented in Table 13 confirm the accuracy of the HWM, FDM, and DQM. Similar accuracy was observed also for cylindrical shells in [16].
The obtained numerical results were validated with the solid element 3D finite element model developed for analysing more complex FGM structures. The results obtained with applying the 3D FEM and HWM were found to be in good (rather excellent) agreement.
Our future studies will focus on the application of the HWM for the analysis of nanostructures and solving fractional differential equations, which are not yet well covered by commercial software solutions. An interesting subtopic, whose research is underway, is adaption of global optimization methods and techniques, developed by the workgroup of design composite structures [31][32][33][34][35][36] to design nano-and graphene structures.