Dynamic analysis of functionally graded material bars by using novel weak form quadrature element method

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 [1-13]. The problems of FGMs have been studied theoretically [1-6], approximately [7], and numerically [8-13].

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 Mindlin-Herrmann rod theory [16], a time domain N-node 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 N-node quadrature bar element

2.1. Expressions of strain energy and kinetic energy

For completeness considerations, the Mindlin-Herrmann theory [10, 16, 17] of a FGM bar with circular cross section shown in Fig. 1 is reviewed. The displacement fields is given by:

1

$u\mathrm{*}\left(x,r,t\right)=u\left(x,t\right),\mathrm{}\mathrm{}\mathrm{}v\mathrm{*}\left(x,r,t\right)=r\psi \left(x,t\right),$

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 three-layer cross section of the FGM bar

For comparisons, the circular cross-section 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:

2

${\epsilon }_{xx}=u^{\mathrm{\text{'}}},\mathrm{}\mathrm{}\mathrm{}{\epsilon }_{rr}={\epsilon }_{\theta \theta }=\psi ,\mathrm{}\mathrm{}\mathrm{}{\gamma }_{xr}=r{\psi }^{\mathrm{\text{'}}},\mathrm{}\mathrm{}\mathrm{}{\gamma }_{x\theta }={\gamma }_{r\theta }=0,$

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

3

${\kappa }_2=\frac{3.383(1+\nu {)}^2}{(2.014+\nu )(0.891+\nu )},\mathrm{}\mathrm{}\mathrm{}{\kappa }_1={\kappa }_2{\left(\frac{0.862+1.14\nu }{1+\nu }\right)}^2,$

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

4

$V=\frac{1}{2}{\int }_0^L\left[E{A}_1(u\mathrm{\text{'}}{)}^2+2K\psi u\mathrm{\text{'}}+E{A}_2{\psi }^2+GJ{\kappa }_1(\psi \mathrm{\text{'}}{)}^2\right]dx,$

where $E{A}_1$, $E{A}_2$, $K$, and $GJ$ are defined by:

5

$\left\{\begin{array}{l}E{A}_1={\int }_A\left[\lambda \right(r)+2\mu (r\left)\right]dA,K=2{\int }_A\lambda \left(r\right)dA,\\ E{A}_2=4{\int }_A\left[\lambda \right(r)+\mu (r\left)\right]dA,GJ={\int }_A{r}^2\mu \left(r\right)dA.\end{array}\right.$

In Eq. (5), $\lambda \left(r\right)$ and $\mu \left(r\right)$ are the well-known Lamé constants defined by:

6

$\lambda \left(r\right)=\frac{E\left(r\right)\nu }{\left(1+\nu \right)\left(1-2\nu \right)},\mathrm{}\mathrm{}\mathrm{}\mu \left(r\right)=\frac{E\left(r\right)}{2\left(1+\nu \right)},$

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:

7

$T=\frac{1}{2}{\int }_0^L\left[\rho A(\dot{u}{)}^2+\rho J{\kappa }_2(\dot{\psi }{)}^2\right]dx.$

In Eq. (7), $\rho A$ and $\rho J$ are computed by:

8

$\rho A={\int }_A\rho \left(r\right)dA,\mathrm{}\mathrm{}\mathrm{}\rho J={\int }_A\rho \left(r\right)r^{2}dA,$

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

9

$F_i\left(r\right)=\left\{\begin{array}{l}{F}_c,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\left(0\le r\le r_c\right),\\ \left(F_o-F_c\right){\left(\frac{r-r_c}{r_i-r_c}\right)}^n+F_c,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\left(r_c<r<r_i\right),\\ F_0,\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\mathrm{}\left(r_i\le r\le r_o\right),\end{array}\right.$

where $n$ is the non-negative 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. N-node quadrature bar element

Consider an N-node weak form quadrature bar element, each node has two degrees of freedom, i.e., $u_j$, ${\psi }_j$ ($j=\text{1,}\text{2,...,}N$). A 9-node weak form quadrature bar element is skematically shown in Fig. 2.

Either Gauss quadrature or Gauss-Lobatto-Legendre (GLL) quadrature can be used for the numerical integration. After performing the quadrature, the strain energy of the bar element can be written as:

10

$V=\frac{LE{A}_1}{4}\sum _{k=1}^N{H}_k\left(\sum _{i=1}^N{A}_{ki}{u}_i\sum _{j=1}^N{A}_{kj}{u}_j\right)+\frac{LK}{2}\sum _{k=1}^N{H}_k\left({\psi }_k\sum _{j=1}^N{A}_{kj}{u}_j\right)$
$+\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 9-node weak form quadrature bar element

Employing the row summation technique [18] and performing the quadrature, the kinetic energy can be expressed as:

11

$T=\frac{\rho AL}{4}\sum _{i=1}^N\left(\sum _{k=1}^N{H}_k{l}_{ik}{\dot{u}}_i^2\right)+\frac{\rho JL{\kappa }_2}{4}\sum _{i=1}^N\left(\sum _{k=1}^N{H}_k{l}_{ik}{\dot{\psi }}_i^2\right),$

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:

12

$\left[m\right]\left\{\ddot{q}\right\}+\left[k\right]\left\{q\right\}=\left\{f\left(t\right)\right\},$

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:

13

${\left\{q\right\}}^T=⌊u_1,u_2,...,u_N,{\psi }_1,{\psi }_2,...,{\psi }_N⌋.$

Then the elements in stiffness matrix [$k$] are:

14

$\begin{array}{l}{k}_{ij}=\frac{{\partial }^{2}V}{\partial u_i\partial u_j},\mathrm{}\mathrm{}\mathrm{}{k}_{i\left(j+N\right)}=\frac{{\partial }^{2}V}{\partial u_i\partial {\psi }_j},\mathrm{}\mathrm{}\mathrm{}{k}_{\left(i+N\right)j}=\frac{{\partial }^{2}V}{\partial {\psi }_i\partial u_j},\mathrm{}\\ k_{(i+N)(j+N)}=\frac{{\partial }^{2}V}{\partial {\psi }_i\partial {\psi }_j},\left(i,j=1,\mathrm{}2,...,N\right).\end{array}$

And the diagonal terms in mass matrix [$m$] are:

15

$m_{ii}=\frac{{\partial }^{2}T}{\partial {\dot{u}}_i\partial {\dot{u}}_i},\mathrm{}\mathrm{}\mathrm{}{m}_{(i+N)(i+N)}=\frac{{\partial }^{2}T}{\partial {\dot{\psi }}_i\partial {\dot{\psi }}_i},\left(i=1,\mathrm{}2,...,N\right).$

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:

16

$\left[k\right]\left\{Q\right\}={\omega }^2\left[m\right]\left\{Q\right\}.$

Solving the above equation by an eigen-value 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 N-node quadrature bar element, the problem presented in [10], i.e., free vibration of a cantilever FGM bar is re-analyzed. The bar is axially uniform and made of Alumina (Al₂O₃) and mild Steel. The elasticity modulus and mass density are $E=$390 GPa and $\rho =$3950 kg/m³ for the Alumina, and $E=$210 GPa and $\rho =$7800 kg/m³ 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 N-node 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.

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 21-node weak form quadrature bar element can yield very accurate frequencies.

When Poisson’s ratio is zero, then Mindlin-Herrmann 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 21-node 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.

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 five-peak tone-burst 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 9-node 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 Mindlin-Herrmann 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 Mindlin-Herrmann rod theory. Therefore, the results based on Mindlin-Herrmann 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₂O₃) and Fig. 10 for bar with material of mild Steel. Poisson’s ratio is 0.3 and Mindlin-Herrmann 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 Mindlin-Herrmann 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 Mindlin-Herrmann 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 N-node 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 21-node 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, Mindlin-Herrmann rod theory should be more reliable and useful in practice for wave propagation analysis.