Abstract
This paper presents an efficient approach to simulate frequency and wave propagations in functionally graded material (FGM) bars by using novel weak form quadrature element method (WQEM). Based on MindlinHerrmann rod theory, a time domain Nnode quadrature bar element is proposed. Detailed formulations are given. Dynamic problems of FGM bars are investigated by using the proposed weak form quadrature bar element. Comparisons are made with results obtained by using frequency domain spectral element method (SEM) and by using strong form differential quadrature method (DQM) to verify the developed quadrature bar element. It is shown that one 21node bar element can yield very accurate frequencies and that the proposed element can efficiently simulate the wave propagation in FGM bars. Compared to results based on the simple rod theory, the results based on MindlinHerrmann rod theory should be more reliable, especially for the wave forms and group velocity.
1. Introduction
The material properties of functionally graded materials (FGMs) can be varied per desired spatial coordinates in such a way as to achieve desired mechanical, thermal, or electrical properties. This capability can help improve the strength and toughness of a structure. Therefore, FGMs have recently received a great attention in fields of aerospace, automobile, as well as electronics [113]. The problems of FGMs have been studied theoretically [16], approximately [7], and numerically [813].
Due to the complicated mathematical structure of the differential equations, it is not an easy task to obtain analytical solutions for general cases [2]. Therefore, various numerical methods, such as spectral element method (SEM) [10], differential quadrature method (DQM) [11], and finite element method (FEM) [10, 12, 13], and meshless method [8, 9] have been used to study the static, buckling, as well as the dynamic response of FGM bars, beams, circular and rectangular plates, cylinders and shells.
Although FEM is a very powerful computational method to deal with complex structural engineering problems, however, very fine meshes are usually needed to obtain accurate solutions, especially at the high frequency region [10]. The weak form quadrature element method (WQEM) [14], formulated based on the principle of minimum potential energy and differential quadrature rule, seems more flexible than the strong form differential quadrature method (DQM) [15]. It is shown that high accuracy can be achieved by the WQEM with smaller number of nodal points [14]. Besides explicit formulas for computing various derivatives of the displacement with respect to space variables at integration points are available, therefore, it is extremely simpler to formulate the stiffness and mass matrix of weak form quadrature elements with different number of nodes.
The purposes of this paper are two folds. One is to develop a weak form quadrature element model which has the capability to accurately predict the dynamics and wave propagations in FGM bars. The other is to investigate the dynamics and wave characteristics of FGM bars, namely, free vibration and wave propagation of FGM bars with radial variation of material properties. Based on MindlinHerrmann rod theory [16], a time domain Nnode quadrature bar element is proposed. Detailed formulations are given. To validate the formulations, results are compared with existing data as well as solutions obtained by using the DQM. Based on the results reported herein conclusions are drawn.
2. Formulation of an Nnode quadrature bar element
2.1. Expressions of strain energy and kinetic energy
For completeness considerations, the MindlinHerrmann theory [10, 16, 17] of a FGM bar with circular cross section shown in Fig. 1 is reviewed. The displacement fields is given by:
where $u(x,t)$ is the axial displacement along the central axis of the bar and $\psi (x,t)$ is the radial contraction.
Fig. 1Sketch of a circular threelayer cross section of the FGM bar
For comparisons, the circular crosssection of the FGM bar schematically shown in Fig. 1 consists of three layers with different materials, namely, the core ($0\le r\le {r}_{c}$), the inner layer (${r}_{c}\le r\le {r}_{i}$), and the outer layer (${r}_{i}\le r\le {r}_{o}$), where ${r}_{o}$ is the outer radius of the circular cross section.
With Eq. (1), strains in the FGM bar can be expressed as:
where the prime denotes the first order derivative with respective to $x$, subscripts $x$, $r$ and $\theta $ denote the axial, radial and circumferential directions, respectively.
To match with three dimensional exact theory, two adjustment coefficients, denoted by ${\kappa}_{\text{1}}$, ${\kappa}_{\text{2}}$ are introduced, namely [10, 16, 17]:
where $\nu $ is Poisson’s ratio of the bar material and is assumed the same for the entire cross section of the bar for simplicity. Note that the values of these two adjustment parameters were selected based on matching the limiting value of the first mode wave and setting a common tangent point for dispersion curves of different Poisson’s ratios [17], thus are approximate values. More details on the two adjustment coefficients may be found in [10, 16, 17].
For linear elastic materials, the strain energy of the FGM bar element with length $L$ is given by [10]:
where $E{A}_{1}$, $E{A}_{2}$, $K$, and $GJ$ are defined by:
In Eq. (5), $\lambda \left(r\right)$ and $\mu \left(r\right)$ are the wellknown Lamé constants defined by:
where $E\left(r\right)$ is the Young’s modulus of the FGM.
The kinetic energy of the FGM bar element with length $L$ is given by:
In Eq. (7), $\rho A$ and $\rho J$ are computed by:
where $\rho \left(r\right)$ is the mass density of the FGM.
Let subscripts $c$, $i$, and $o$ represent the core, inner layer, and the outer layer. Assumption is made that elasticity moduli and mass densities in the core and the outer layer are uniform and that those in the inner layer are varying in the radial direction according to the power law given by [10]:
where $n$ is the nonnegative power law exponent, and ${F}_{i}\left(r\right)$ represents either the elasticity modulus ${E}_{i}\left(r\right)$ or the mass density ${\rho}_{i}\left(r\right)$ of the inner layer.
With the help of symbolic computing software Maple@, explicit formulas of $E{A}_{1}$, $E{A}_{2}$, $K$, $GJ$, in Eq. (4) and $\rho A$, $\rho J$ in Eq. (7) can be easily obtained by substituting Eq. (9) into Eq. (5) and Eq. (8). Detailed expressions can be found in [10].
2.2. Nnode quadrature bar element
Consider an Nnode weak form quadrature bar element, each node has two degrees of freedom, i.e., ${u}_{j}$, ${\psi}_{j}$ ($j=\text{1,}\text{}\text{2,...,}\text{}N$). A 9node weak form quadrature bar element is skematically shown in Fig. 2.
Either Gauss quadrature or GaussLobattoLegendre (GLL) quadrature can be used for the numerical integration. After performing the quadrature, the strain energy of the bar element can be written as:
$+\frac{LE{A}_{2}}{4}\sum _{k=1}^{N}{H}_{k}\left({\psi}_{k}\mathrm{*}{\psi}_{k}\right)+\frac{LGJ{\kappa}_{1}}{4}\sum _{k=1}^{N}{H}_{k}\left(\sum _{i=1}^{N}{A}_{ki}{\psi}_{i}\sum _{j=1}^{N}{A}_{kj}{\psi}_{j}\right),$
where ${A}_{ki}$ are weighting coefficient of the first order derivatives with respect to $x$ at the integration point ${x}_{k}$, and ${H}_{k}$ is the weight of the quadrature. Explicit formulas to compute ${A}_{ki}$ are available. Details may be found in [14].
Fig. 2Sketch of a 9node weak form quadrature bar element
Employing the row summation technique [18] and performing the quadrature, the kinetic energy can be expressed as:
where ${l}_{ik}$ are the value of the Lagrange interpolation function ${l}_{i}\left(x\right)$ at the integration point ${x}_{k}$. If the integration points are also used as the nodal points of the element, one has ${l}_{ik}={\delta}_{ik}$.
The quadrature element equation of motion can be symbolically written as:
where the double over dots denote the second order derivative with respect to time, and [$m$] and [$k$] are mass matrix and stiffness matrix, respectively.
If the nodal displacement vector is defined as:
Then the elements in stiffness matrix [$k$] are:
And the diagonal terms in mass matrix [$m$] are:
For free vibration analysis, $\left\{f\left(t\right)\right\}=\left\{0\right\}$ and $\left\{q\right\}=\left\{Q\right\}\mathrm{s}\mathrm{i}\mathrm{n}\omega t\text{,}$ where $\omega $ is the circular frequency. Thus Eq. (12) can be simplified as:
Solving the above equation by an eigenvalue solver yields the frequencies and their corresponding mode shapes.
For wave propagations, Eq. (12) can be effectively integrated by using the central finite difference method, since the mass matrix is in diagonal form.
3. Results and discussion
To verify the proposed Nnode quadrature bar element, the problem presented in [10], i.e., free vibration of a cantilever FGM bar is reanalyzed. The bar is axially uniform and made of Alumina (Al_{2}O_{3}) and mild Steel. The elasticity modulus and mass density are $E=$390 GPa and $\rho =$3950 kg/m^{3} for the Alumina, and $E=$210 GPa and $\rho =$7800 kg/m^{3} for the mild Steel. Poisson’s ratio is 0.3 for the entire bar. Geometries are $L=$ 2 m, ${r}_{o}$, ${r}_{i}$ and ${r}_{c}$ are 0.01 m, 0.009 m, and 0.001 m, respectively.
For comparisons, two types of FGM axial bars are considered. One type, denoted by ‘${E}_{c}>{E}_{o}$’, is that the core and the outer layer are made of Alumina and mild Steel. The other, denoted by ‘${E}_{c}<{E}_{o}$’, is that the core and the outer layer are made of mild Steel and Alumina. The bar is clamped at $x=$ 0 and free at $x=L$. Only one Nnode quadrature bar element is used in the free vibration analysis. Thus the qudrature element mesh is similar to the one shown in Fig. 2, the only difference from the one shown in Fig. 2 is the number of nodes.
To show the convergence, different $N$ is used. Results and comparisons are listed in Table 1 and Table 2 for the two types of FGM.
Table 1Comparison of the natural frequencies (Ec>Eo, n= 1)
Mode  1  2  3  5  10 
Present ($N=$ 15) Present ($N=$ 19) Present ($N=$ 21)  789.544 789.544 789.544  2368.59 2368.59 2368.59  3947.50 3947.50 3947.50  7104.59 7104.59 7104.59  15137.4 14987.0 14988.9 
DQM ($N=$ 21) DQM($N=$ 101)  789.544 789.544  2368.59 2368.59  3947.50 3947.50  7104.59 7104.59  14981.9 14988.9 
SEM [10]  789.955  2369.84  3949.65  7108.82  15001.8 
In Table 1 and Table 2, the data cited from [10] are obtained by using one frequency domain spectral element. It is seen that the results obtained by using the DQM with the same number of grid points are similar to the ones by the weak form quadrature element method. Only the tenth mode frequency is sightly different. To improve the solution accuracy, 101 grid points are used in the DQ analysis. The DQ results with $N=$ 101 should be accurate to be servered as the reference data. As expected, the data are not changed except for the the tenth mode frequency. It is seen that one 21node weak form quadrature bar element can yield very accurate frequencies.
Table 2Comparison of the natural frequencies (Ec<Eo, n= 1)
Mode  1  2  3  5  10 
Present ($N=$ 15) Present ($N=$19) Present ($N=$ 21)  1011.92 1011.92 1011.92  3035.73 3035.73 30335.73  5059.43 5059.43 5059.43  9106.21 9106.21 9106.21  19407.3 19213.7 19216.2 
DQM ($N=$ 21) DQM($N=$ 101)  1011.92 1011.92  3035.73 3035.73  5059.43 5059.43  9106.21 9106.21  19207.0 19216.1 
SEM [10]  1012.50  3037.52  5062.58  9112.87  19240.7 
When Poisson’s ratio is zero, then MindlinHerrmann rod theory is reduced to the simple bar theory to some extent. The results are listed in Table 3 and Table 4. It is seen that one 21node quadrature bar element can yield exactly the same results as the ones obtained by the DQM with 101 points. Comparisons of the results to the corresponding ones listed in Table 1 and Table 2 reveal that the difference caused by the two theories is very small and can be neglected in engineering applications.
Table 3The natural frequencies based on simple theory (Ec>Eo, n= 1)
Mode  1  2  3  5  10 
Present ($N=$ 21)  789.546  2368.64  3947.73  7105.91  15001.3 
DQM ($N=$ 101)  789.546  2368.64  3947.73  7105.91  15001.4 
Table 4The natural frequencies based on simple theory (Ec<Eo, n= 1)
Mode  1  2  3  5  10 
Present ($N=$ 21)  1011.93  3035.78  5059.63  9107.33  19226.5 
DQM ($N=$ 101)  1011.93  3035.78  5059.63  9107.33  19226.5 
The first two mode shapes are shown in Fig. 3 for the cases of ${E}_{c}>{E}_{o}$ ($n=$ 1) and Fig. 4 for ${E}_{c}<{E}_{o}$ ($n=$ 1). The shapes are similar for the two cases. The variation of the radial contraction, $\psi \left(x\right)$, can be clearly seen. At the free end, $\psi $ is zero and at the fixed end, $\psi $ reaches its maximum value. Therefore, care should be taken to apply the fixed end boundary conditions, only $u=$0 is applied and $\psi \left(0\right)$ is not zero.
Next the axial wave propagations in two types of FGM bars are investigated by using the proposed quadrature bar element. Both ends of the bar are free. The geometries and material properties are the same as the ones in the free vibration analysis. The excitation signal is a fivepeak toneburst wave with a center frequency of 100 kHz [18] and is applied at $x=$0. The axial waves are measured at the middle point of the bar ($x=L/$2).
Fig. 3The first and second mode shapes for the case of Ec>Eo (n= 1)
a)
b)
Fig. 4The first and second mode shapes for the case of Ec<Eo (n= 1)
a)
b)
In the numerical simulations, 60 equal length 9node bar elements are used. Thus the model contains 481 nodes and 962 degrees of freedom. Since the mass matrix is a diagonal matrix, thus the equation of motion is explicitly integrated step by step by using the central difference method. A time increment of 5×10^{8} s is used in the numerical integration. Differnt power law exponents are adopted in the investigations. Simulations are plotted in Fig. 5 and Fig. 6. Normalized amplitudes are shown.
Fig. 5 displays the axial waves measured at $x=L/$2 with type 1 FGM (${E}_{c}>{E}_{o}$). The axial waves are propagating in the FGM bars and then being reflected repeatedly from two free ends. The group velocity of the axial wave has the smallest value when $n=$ 0.25. The group velocity of the axial wave increases with the increase of $n$. Thus Fig. 5(c) shows that the axial wave when $n=$4 is propagating faster than those when $n=$ 1 and $n=$ 0.25. Severe dispersion is observed.
Fig. 6 shows the axial waves measured at $x=L/$2 with type 2 FGM (${E}_{c}<{E}_{o}$). The axial waves are propagating in the FGM bars and then being reflected repeatedly from two free ends. The group velocity of the axial wave has the largest value when $n=$0.25. The group velocity of the axial wave decreases with the increase of $n$, different from the previous cases. Thus Fig. 6(a) shows that the axial wave when $n=$ 0.25 is propagating faster than those when $n=$ 1 and $n=$ 4. Besides the dispersion is not as severe as the one shown in Fig. 5.
Fig. 5Axial waves in the FGM axial bars with different n (Ec>Eo)
a)
b)
c)
It has been shown that the effect on the natural frequencies is negligible whether MindlinHerrmann rod theory or simple rod theory is used. To investigate the effect of the two different theories on the characteristic of wave propagation in the FGM bar, the axial waves measured at the middle point of the bar ($n=$1) are shown in Fig. 7 for ${E}_{c}>{E}_{o}$ and Fig. 8 for ${E}_{c}<{E}_{o}$. Simple rod theory is used.
It is seen that the wave forms based on simple rod theory, whether ${E}_{c}>{E}_{o}$ or ${E}_{c}<{E}_{o}$, are not dispersive. Comparisons of Fig. 7 and Fig. 8 with Fig. 5(b) and Fig. 6(b) reveal that the group velocity based on the simple rod theory is faster than the one based on MindlinHerrmann rod theory. Therefore, the results based on MindlinHerrmann rod theory should be more reliable.
Fig. 6Axial waves in the FGM axial bars with different n (Ec<Eo)
a)
b)
c)
Fig. 7Axial waves in the FGM axial bars based on simple rod theory (Ec>Eo)
Fig. 8Axial waves in the FGM axial bars based on simple rod theory (Ec<Eo)
To give an explanation why the wave form for ${E}_{c}>{E}_{o}$ is more dispersive than the one for ${E}_{c}<{E}_{o}$, the axial waves measured at the middle point of the bar are plotted in Fig. 9 for bar with material of Alumina (Al_{2}O_{3}) and Fig. 10 for bar with material of mild Steel. Poisson’s ratio is 0.3 and MindlinHerrmann rod theory is used.
Fig. 9Axial waves in the Alumina (Al2O3) bar
Fig. 10Axial waves in the mild Steel bar
From Fig. 9 and Fig. 10, it is seen that the predicted wave forms based on MindlinHerrmann rod theory show dispersive even for bars with isotropic materials. The wave form in the bar with material of mild Steel is more dispersive than the one in the bar with material of Alumina. This is the reason why the wave form in the FGM bar is more dispersive for ${E}_{c}>{E}_{o}$ than that for ${E}_{c}<{E}_{o}$, since the FGM bar contains more material of mild Steel.
It should be mentioned that the waves of $\psi (L/2,t)$ are similar to the corresponding axial waves shown in Figs. 5 and 6, thus are ommitted for the space limitations.
4. Conclusions
Based on MindlinHerrmann rod theory, a novel weak form quadrature element method is proposed in this paper to analyze the dynamic behavior of FGM bars. A time domain Nnode weak form quadrature bar element is formulated in detail. Dynamic problems of FGM bars are investigated by using the proposed quadrature bar element. For verifications, comparisons are made to results obtained by using frequency domain spectral element method and strong form differential quadrature method with large number of grid points. It is shown that one 21node bar element can yield very accurate frequencies and that the proposed element can accurately simulate the wave propagation in FGM bars. Compared to results based on the simple rod theory, the difference on the nature frequencies corresponding to the axial vibration is negligible. However, the wave forms and group velocity are quite different. Therefore, MindlinHerrmann rod theory should be more reliable and useful in practice for wave propagation analysis.
References

Tutuncu N., Temel B. A novel approach to stress analysis of pressurized FGM cylinders, disks and spheres. Composite Structures, Vol. 91, 2009, p. 385390.

Murín J., Aminbaghai M., Kuti V. Exact solution of the bending vibration problem of FGM beams with variation of material properties. Engineering Structures, Vol. 32, 2010, p. 16311640.

Arslan E., Eraslan A. N. Bending of graded curved bars at elastic limits and beyond. International Journal of Solids and Structures, Vol. 50, 2013, p. 806814.

Huang Y., Li X. F. Buckling analysis of nonuniform and axially graded columns with varying flexural rigidity. Journal of Engineering Mechanics, Vol. 137, Issue 1, 2011, p. 7381.

Torabi J., Kiani Y., Eslami M. R. Linear thermal buckling analysis of truncated hybrid FGM conical shells. Composites Part B: Engineering, Vol. 50, 2013, p. 265272.

AkgÖz B., Civalek Ö. Longitudinal vibration analysis of strain gradient bars made of functionally graded materials (FGM). Composites Part B: Engineering, Vol. 55, 2013, p. 263268.

Pradhan K. K., Chakraverty S. Free vibration of Euler and Timoshenko functionally graded beams by RayleighRitz method. Composites Part B: Engineering, Vol. 51, 2013, p. 175184.

Ferreira A. J. M., Batra R. C., Roque C. M. C., et al. Natural frequencies of functionally graded plates by a meshless method. Composite Structures, Vol. 75, Issue 14, 2006, p. 593600.

Neves A. M. A., Ferreira A. J. M., Carrera E., et al. A quasi3D sinusoidal shear deformation theory for the static and free vibration analysis of functionally graded plates. Composites Part B: Engineering, Vol. 43, Issue 2, 2012, p. 711725.

Hong M., Park I., Lee U. Dynamics and waves characteristics of the FGM axial bars by using spectral element method. Composite Structures, Vol. 107, 2014, p. 585593.

Xiang H. J., Yang J. Free and forced vibration of a laminated FGM Timoshenko beam of variable thickness under heat conduction. Composites Part B, Vol. 39, Issue 2, 2008, p. 292303.

Yu Z., Chu F. Identification of crack in functionally graded material beams using the pversion of finite element method. Journal of Sound and Vibration, Vol. 325, Issue 12, 2009, p. 6984.

Alshorbagy A. E., Eltaher M. A., Mahmoud F. F. Free vibration characteristics of a functionally graded beam by finite element method. Applied Mathematical Modelling, Vol. 35, Issue 1, 2011, p. 41225.

Jin C. H., Wang X., Ge L. Y. Novel weak form quadrature element method with expanded Chebyshev nodes. Applied Mathematics Letters, Vol. 34, 2014, p. 5159.

Wang X., Wang Y., Yuan Z. Accurate vibration analysis of skew plates by the new version of the differential quadrature method. Applied Mathematical Modelling, Vol. 38, 2014, p. 926937.

Mindlin R., Herrmann G. A onedimensional theory of compressional waves in an elastic rod. Proceedings of the First US national congress of applied mechanics. Chicago, Illinois, 1951.

Yu C. P., Roesset J. M. Dynamic stiffness matrices for linear members with distributed mass. Tamkang Journal of Science and Engineering, Vol. 4, Issue 4, 2001, p. 253264.

Ge L. Y., Wang X., Wang F. Accurate modeling of PZTinduced Lamb wave propagation in structures by using a novel spectral finite element method. Smart Materials and Structures, Vol. 23, 2014, p. 09501814.
About this article
The project is partially supported by Funding of Jiangsu Innovation Program for Graduate Education (KYLX_0219), the Fundamental Research Funds for the Central Universities, and the Priority Academic Program Development of Jiangsu Higher Education Institutions.