Numerical simulation of elastic wave propagation in functionally graded cylinders using time-domain spectral finite element method

Abstract

Functionally graded material is a new class of composite materials with gradually varying material contents. In this article, a two-dimensional axisymmetric time-domain spectral finite element method is proposed for the analysis of elastic wave propagation in functionally graded cylinders. The stiffness and mass matrix of an arbitrary quadrilateral axisymmetric spectral element is developed. The element nodes are collocated at Gauss–Lobatto–Legendre points. High-order Lagrangian polynomials are used as shape functions and the Gauss–Lobatto–Legendre quadrature rules are applied. They provide the capability to model spatial varying material properties and lead to a diagonal mass matrix. The effectiveness of the proposed model is verified by comparing with published data. Elastic wave propagation in a functionally graded cylinder subjected to an impulsive loading is studied detailed. The effect of material grading pattern is analyzed. The results demonstrate the effectiveness of the present method for the analysis of elastic wave propagation in functionally graded solids with axial symmetry and the material composition variation has an important effect on structural wave propagation behavior.

Introduction

Functionally graded material (FGM) is a new kind of composite materials with gradually varying material properties tailored to specific applications. Compared with the common layered composite material, one of the superior properties of FGM is that the continuous gradation in material content can overcome the interfacial problem arisen from sudden change in material properties at the interfaces of different phases. ¹ FGMs have been increasingly used in aerospace, nuclear, automotive and defense industries. In many of these applications, FGM structures usually serve as critical parts in the whole assemble. However, on the other hand, these FGM structures are often subjected to extreme loadings and prone to becoming damaged. To effectively detect damages in FGM structures which is critical to enable system reliability and safety, many damage detection methods have been proposed in recent years, such as vibration-based methods, and ultrasonic-based methods. Among these current methods, the guided wave damage detection method has attracted much attention due to its capability for online and long range interrogation of the structures. This method is based on the active or passive elastic wave stress propagation in the structures. To better implement the method, it is important to analyze and understand the wave propagation phenomenon in various structures. Therefore, the analysis of wave propagation in FGM structures is issued in this article.

A number of analytical studies have been carried out to study elastic wave propagation in functionally graded (FG) structures.^2–7 An infinite FG plate subjected to a point impact was studied by Sun and Luo. ² High-order shear deformation plate theory (HSDT) and Hamilton principle were used to obtain approximate analytical solutions in their work. Dispersion and power flow solutions for a FGM plate were given using the extension of Legendre polynomial approach by Lefebvre et al. ³ Lamb wave propagation in a FGM plate was studied based on the power series technique by Cao et al. ⁴ A semi-analytical method, which combines the state space method and the one-dimensional (1D) differential quadrature method, was developed to obtain the dynamic responses of FG circular plates by Nie and Zhong ⁵ A simplified 1D model is proposed to develop insight into the stress wave propagation in FGMs by Bruck ⁶ The free and forced vibrations of a FGM annular sector plate with simply supported radial edges and arbitrary circular edges were analyzed by Nie and Zhong ⁷ Although these analytical works can give some insight into the elastic wave propagation in FGM structures, they are limited to very simple boundary and loading conditions.

In recent years, many works on the analysis of wave propagation in complex FG structures by the conventional finite element method have been reported.^8–11 For example, a beam element based on first-order shear deformation theory was proposed to study wave propagation behaviors of FGMs by Chakraborty et al. ⁸ The analysis of a FG hollow cylinder under dynamic impact was carried out in Shakeri and colleagues,^9,10 where the FG cylinder was assumed to be made of many sub-cylinders. Each sub-cylinder was considered as an isotropic layer and the material properties in each layer were either constant ⁹ or assumed to vary linearly in the thickness direction. ¹⁰ Thus, FGM properties were resulted by suitable arrangement of layers in multilayer cylinder. The axisymmetric finite element model was applied to study the wave propagation in a thick hollow cylinder with finite length made of two-dimensional (2D) FGM subjected to impact internal pressure by Asgari et al. ¹¹ Although it is validated that the conventional finite element method (FEM) has the ability to simulate the elastic wave propagation in FGM structures with arbitrary geometry and boundary conditions, the computational cost is often very high. The reason is that at high frequency, extremely precise finite element discretization in time and space domain is required to yield a converged dynamics solution. Hence, alternatives to the conventional FEM have been considered.

Recently, the time-domain spectral finite element method (SFEM), which was first proposed by Patera ¹² in fluid dynamics, has been extended to analyze the problems of structural wave propagation.^13–19 This method can be viewed as a special type of finite element methods and bears some similarity with the p-version FEM. The key idea of time-domain SFEM is the adoption of specific high-order shape functions. Local element nodes are distributed non-uniformly at the Gauss–Lobatto–Legendre (GLL) points. The Lagrangian interpolation polynomials pass through the GLL points. The GLL quadrature rules are applied to deduce element mass and stiffness matrices. As a result, high interpolation accuracy is achieved. Besides, the element mass matrix is diagonal due to the orthogonality of approximation functions. Hence, the time-domain SFEM is more effective compared to the conventional FEM for wave propagation problems. It is also noted that using high-order shape functions in SFEM can give better approximation of gradually varying material properties. Therefore, the time-domain SFEM is particularly suitable for the wave propagation analysis of complex FGM structures. However, the related works are few except that in Hedayatrasa et al., ²⁰ a 2D plane strain time-domain SFEM was proposed to analyze the lamb wave propagation in FGMs. In addition, it is noted that beside the above time-domain spectral element methods, the transformed-domain SFEM methods have also been proposed to effectively solve the problems of elastic wave propagation.^21,22

The main objective of this article is studying the elastic wave propagation in functionally grade cylinders by an efficient time-domain SFEM. An axisymmetric solid time-domain spectral finite element model is developed. The proposed SFEM, which combines the geometrical flexibility of the conventional FEM and the rapid convergence of spectral method, is capable to accurately analyze axial and radial oscillations with high-order approximation functions. This article is organized as follows. In section “Theoretical background,” the related theoretical background is introduced. In section “Spectral finite element modeling,” the axisymmetric solid spectral finite element formulations for the wave propagation in FG cylinders are deduced. The main purpose of section “Model validation” is to validate the proposed method by comparing with published data. In section “Numerical example and discussions,” elastic wave propagation in a FGM cylinder is studied in detail. In the last section, the conclusions are addressed.

Theoretical background

Basic equations

Consider a FG cylinder as shown in Figure 1. A cylindrical coordinate system (r, θ, and z) is applied to describe the cylinder geometry and deformation. The thickness of the cylinder is L and the radius is R. The cylinder is made of 1D FGM. The geometry of the cylinder and the loading considered here are not the functions of circumferential coordinate θ. It is known that there are an infinite number of guided wave modes in a cylindrical waveguide. Torsional waves can also be created under axisymmetric conditions. In this article, however, only axisymmetric wave modes are considered. In this case, the deformation becomes symmetrical with respect to the z-axis and the wave propagation problem can be considered as a plane axisymmetric problem. Thus, by neglecting the body force, the equilibrium differential equations written in 2D cylindrical coordinate (r, z) are ²³

$$\frac{\partial {\sigma }_{rr}}{\partial r}+\frac{\partial {\tau }_{rz}}{\partial z}+\frac{{\sigma }_{rr}-{\sigma }_{\theta \theta }}{r}=\rho \frac{{\partial }^{2}u}{\partial t^2}$$

(1)

and

$$\frac{\partial {\tau }_{rz}}{\partial r}+\frac{\partial {\sigma }_{zz}}{\partial z}+\frac{{\tau }_{rz}}{r}=\rho \frac{{\partial }^{2}w}{\partial t^2}$$

(2)

where t is time; $u=u(r,z)$ and $w=w(r,z)$ denote the displacement components along the radial and axial directions, respectively; ${\sigma }_{rr}{\sigma }_{zz}{\sigma }_{\theta \theta }$ , and ${\tau }_{rz}$ are the stress components. The deformations of the cylinder are assumed to be small. The linear strain–displacement equations in the cylindrical coordinate system are given by

$$\mathit{\epsilon }=\left[\begin{array}{c}\hfill \begin{array}{c}\hfill {\epsilon }_{rr}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{\theta \theta }}\end{array}\hfill \\ \multicolumn{1}{c}{{\epsilon }_{zz}}\\ \multicolumn{1}{c}{{\gamma }_{rz}}\end{array}\right]=\left[\begin{array}{c}\hfill \begin{array}{c}\hfill \frac{\partial u}{\partial r}\hfill \\ \multicolumn{1}{c}{\frac{u}{r}}\end{array}\hfill \\ \multicolumn{1}{c}{\frac{\partial w}{\partial z}}\\ \multicolumn{1}{c}{\frac{\partial u}{\partial z}+\frac{\partial w}{\partial r}}\end{array}\right]$$

(3)

where ${\epsilon }_{rr},{\epsilon }_{\theta \theta },{\epsilon }_{zz}$ , and ${\gamma }_{rz}$ are the strain components. According to Hooke’s law, the constitutive equations for a linearly elastic material in matrix form are

$$\mathit{\sigma }=\left[\begin{array}{c}\hfill \begin{array}{c}\hfill {\sigma }_{rr}\hfill \\ \multicolumn{1}{c}{{\sigma }_{\theta \theta }}\end{array}\hfill \\ \multicolumn{1}{c}{{\sigma }_{zz}}\\ \multicolumn{1}{c}{{\sigma }_{rz}}\end{array}\right]=\mathit{D}\mathit{\epsilon }$$

(4)

where the coefficients of elasticity is

$$\mathit{D}=\frac{E}{(1+\mathrm{\mu })(1-2\mathrm{\mu })}\left[\begin{array}{c}\hfill \begin{array}{c}\hfill \begin{array}{cc}\hfill 1-\mathrm{\mu }\hfill & \hfill \mathrm{\mu }\hfill \end{array}\begin{array}{cc}\hfill \mathrm{\mu }\hfill & \hfill 0\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\hfill \mathrm{\mu }\hfill & \hfill 1-\mathrm{\mu }\hfill \end{array}\begin{array}{cc}\hfill \mathrm{\mu }\hfill & \hfill 0\hfill \end{array}}\end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\hfill \mathrm{\mu }\hfill & \hfill \mathrm{\mu }\hfill \end{array}\begin{array}{cc}\hfill 1-\mathrm{\mu }\hfill & \hfill 0\hfill \end{array}}\\ \multicolumn{1}{c}{\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \end{array}\begin{array}{cc}\hfill 0\hfill & \hfill \frac{1-2\mathrm{\mu }}{2}\hfill \end{array}}\end{array}\right]$$

(5)

In equation (5), E is Young’s modulus and $\mathrm{\mu }$ is Poisson’s ratio. For FGM, both are the functions of space coordinates. The detailed material model for FGM will be given in next section.

Figure 1.

Geometry of a functionally graded cylinder.

Material properties

The cylinder is made of a combined ceramic–metal material and the material composition is graded continuously along either its thickness direction (z-direction) or its radial direction (r-direction). Let V_m and V_c denote the volume fraction of the metal phase and the ceramic phase, respectively. According to the homogenization scheme in Reddy and Chin, ²⁴ V_m and V_c are related as

$$V_m+V_c=1$$

(6)

For the case of the cylinder’s material is graded through the z-direction, it is assumed that the top surface of the cylinder is pure ceramic and the bottom surface is pure metal. The variation of the volume faction of the ceramic can be defined by the well-known power law

$$V_c={\left(\frac{z}{L}\right)}^n$$

(7)

where z denotes the z-coordinate of the cylinder and n is a non-negative constant. Using the classical rule of mixture, the distribution of the material properties can be expressed as

$$P\left(z\right)=(P_c-P_m){\left(\frac{z}{L}\right)}^n+P_m$$

(8)

where P is the material properties such as Young’s modulus, Poisson’s ratio, and the mass density $\rho$ .

Equation (8) indicates that at the top surface boundary (z = L) of the cylinder, a typical material property takes on the value P = P_c. Similarly, at the bottom surface boundary (z = 0) of the cylinder, a typical material property takes on the value P = P_m. Figure 2 shows the influence of power law index n on volume fraction. It is noted that in the case of n = 1, the material property is linear variation along the thickness direction.

Figure 2.

The influence of compositional gradient exponents on volume fraction variation.

For the case of the cylinder’s material is graded through the r-direction, the grading pattern can be modeled similarly as equation (8), in which the variable z is substituted by variable r and the thickness L is substituted by the cylinder’s radius, R.

Spectral finite element modeling

Due to axisymmetry, the FG cylinder is approximated by finite elements in the plane of revolution (r-z plane). A 2D arbitrary quadrilateral spectral element is proposed to mesh the structure. The element is based on the high-order Lagrange polynomial interpolation. Figure 3 shows a typical 2D spectral element of fourth order in local and global coordinate system. It is noted that unlike the conventional finite element, the nodes of the spectral element are non-uniformly spaced. It is very important to avoid Runge’s phenomenon, which refers to a problem of oscillation near the edges of an interval that occurs when using polynomial interpolation with polynomials of high degree over a set of nodes spaced in a uniform grid. ²⁵ For the proposed spectral element, the local nodes are located at the GLL points, which are the roots of the following equations

$$(1-{\xi }^2)P{\prime }_n\left(\xi \right)=0\mathrm{and}(1-{\eta }^2)P{\prime }_n\left(\eta \right)=0$$

(9)

where $P{\prime }_n(\xi )$ and $P{\prime }_n(\eta )$ represent the first derivative of the nth order Legendre polynomial in the $\xi$ direction and $\eta$ direction, respectively. It is assumed that the interpolation order $n_{\xi }$ in the $\xi$ direction is the same as the interpolation order $n_{\eta }$ in the $\eta$ direction. By solving equation (9), the local coordinates of all nodes can be determined.

Figure 3.

Typical 2D spectral element of fourth order in local (a) and global (b) coordinate system.

The interpolation of the displacement component u in the local coordinates $(\xi ,\eta )$ is defined as

$$u=\sum _{i=1}^{n_{\xi }+1}\sum _{j=1}^{n_{\eta }+1}{N}_i\left(\xi \right)N_j\left(\eta \right)u_{ij}$$

(10)

where $N_i$ is the nth order Lagrange interpolating polynomials in 1D. For the 2D spectral element of order n, the total node number is m = (n + 1)². The corresponding 2D shape function can be written as

$$N_k(\xi ,\eta )=N_i\left(\xi \right)N_j\left(\eta \right)$$

(11)

where k = 1,…, m; i, j = 1,…, n + 1. Therefore, the element displacement field can be written in the form of matrix as

$$\begin{array}{cc}\multicolumn{1}{c}{\mathit{u}=\left\{\begin{array}{c}\hfill u\hfill \\ \multicolumn{1}{c}{w}\end{array}\right\}=}& [\begin{array}{cc}\hfill N_1(\xi ,\eta )\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill N_1(\xi ,\eta )\hfill \end{array}\dots \begin{array}{cc}\hfill N_m(\xi ,\eta )\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{0}& \hfill N_m(\xi ,\eta )\hfill \end{array}]\hfill \\ \left[\begin{array}{c}\hfill \begin{array}{c}\hfill u_1\hfill \\ \multicolumn{1}{c}{w_1}\\ \multicolumn{1}{c}{u_2}\end{array}\hfill \\ \multicolumn{1}{c}{w_2}\\ \multicolumn{1}{c}{⋮}\\ \multicolumn{1}{c}{u_m}\\ \multicolumn{1}{c}{w_m}\end{array}\right]=\mathit{Nq}\hfill \end{array}$$

(12)

with N denotes the element matrix of the shape functions and q is the element nodal displacement vector.

Substituting equation (12) in equation (3), the element strain matrix can be obtained and defined as

$$\mathit{\epsilon }=\mathit{Bq}$$

(13)

where $B=[B_1{B}_2\cdots B_k\cdots B_m]$ is the strain–displacement matrix and sub-matrix B _k is

$$\left[B_k\right]=\left[\begin{array}{c}\hfill \begin{array}{cc}\hfill \frac{\partial N_k}{\partial r}\hfill & \hfill 0\hfill \\ \multicolumn{1}{c}{\frac{N_k}{r}}& \hfill 0\hfill \end{array}\hfill \\ \multicolumn{1}{c}{\begin{array}{cc}\hfill 0\hfill & \hfill \frac{\partial N_k}{\partial z}\hfill \\ \multicolumn{1}{c}{\frac{\partial N_k}{\partial z}}& \hfill \frac{\partial N_k}{\partial r}\hfill \end{array}}\end{array}\right]$$

(14)

The element stress matrix can be deduced by combining equations (4) and (13)

$$\mathit{\sigma }=\mathit{DBq}$$

(15)

The resulted element mass and stiffness matrix M^e and K^e can be achieved by applying Hamilton principles for each element

$${\mathit{M}}^e=2\pi \times \underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}\rho (\xi ,\eta ){\mathit{N}}^T(\xi ,\eta )\mathit{N}(\xi ,\eta )r\det \left(\mathit{J}\right)d\xi d\eta$$

(16)

$${\mathit{K}}^e=2\pi \underset{-1}{\overset{1}{\int }}\underset{-1}{\overset{1}{\int }}{\mathit{B}}^T\mathit{DB}r\det \left(\mathit{J}\right)d\xi d\eta$$

(17)

where $\rho$ and D represent the mass density and the elastic material property matrix, respectively, which are the functions of space coordinates; J is Jacobian matrix that denotes the mapping from the global coordinates (r, z) to the reference coordinates ( $\xi ,\eta$ ), defined as

$$\left[\mathit{J}\right]=\left[\begin{array}{cc}\hfill \frac{\partial r}{\partial \xi }\hfill & \hfill \frac{\partial z}{\partial \xi }\hfill \\ \multicolumn{1}{c}{\frac{\partial r}{\partial \eta }}& \hfill \frac{\partial z}{\partial \eta }\hfill \end{array}\right]$$

(18)

The element matrices $M^e$ and $K^e$ can be calculated numerically using the GLL integration rules as follows

$$M^e=2\pi \sum _{i=1}^n{\omega }_i\sum _{j=1}^n{\omega }_j\rho (\xi ,\eta ){\mathit{N}}^T(\xi ,\eta )\mathit{N}(\xi ,\eta )rdet\left[J(\xi ,\eta )\right]$$

(19)

$$\begin{array}{c}{K}^e=2\pi {\displaystyle \sum }_{i=1}^n{\omega }_i{\displaystyle \sum }_{j=1}^n\hfill \\ {\omega }_j\mathit{B}{(\xi ,\eta )}^T\mathit{D}(\xi ,\eta )\mathit{B}(\xi ,\eta )rdet\left(J\right)rdet\left[J(\xi ,\eta )\right]\hfill \end{array}$$

(20)

where ${\omega }_i$ and ${\omega }_j$ represent the quadrature weights in two directions $\xi ,\eta$ , respectively. The relationship between the weights and shape function is

$$\sum _{i=1}^n{\omega }_i{N}_n{N}_m=\{\begin{array}{c}\hfill C_{n}ifn=m\hfill \\ \multicolumn{1}{c}{0ifn\ne m}\end{array}$$

(21)

This significant orthogonality can result in a diagonal mass matrix when implementing the GLL integration. It is also noted that using GLL integration can decrease the complexity and improve the computation efficiency, but this integration scheme will typically result in “under-integration,” which may underestimate values of computed eigenfrequencies of the analyzed systems, when compared with other more accurate integration scheme like the Gauss integration.

It is noted that the radial distance r in matrix B is assumed variable. The material properties are evaluated at each GLL integration points. This can be realized that first transform the integration point coordinate in the local element coordinate system to the global coordinate system and then use the spatial coordinates in the global system to calculate the corresponding radial distance and material properties. Hence, material properties can be obtained for each GLL integration point while the numerical integration is performed.

It can be found from equation (14) that when r becomes very small, the matrix B _k becomes singularity, which results in the difficulty of calculating the element matrices. In this article, two measures are presented to solve this problem.

For the stiffness matrix calculation, the singularity can be eliminated by imposing exact boundary conditions. In plane axisymmetric problem, the ideal boundary conditions can be expressed as

$$u=0,\mathrm{along}r=0$$

(22)

The boundary condition is introduced to the stiffness matrix during calculation. The rows and columns of stiffness matrix corresponding to the r = 0 boundary degrees of freedoms are imposed as zero value, except for elements on the diagonal of the matrix, which should be fixed as one. ¹³ By this means, the effect of r = 0 are removed and the singularity of stiffness matrix can be solved.

For the mass matrix calculation, according to equation (16), the elements of mass matrix relative to r = 0 are zeros. For example, an element whose left boundary coincides with the symmetry axis is shown in Figure 4. The elements of this element mass matrix which are related to the green integration points are zeros. This will bring singularity when using the mass matrix for the process of direct numerical integration. In this article, an average mass distribution method is proposed. As shown in Figure 4, consider the blue integration points near the green integration points. The half of the mass concentrated at the blue points is given to the green integration points. Thus, the mass at the green points are not zeros. This method can eliminate zero elements in the main diagonal elements of mass matrix.

Figure 4.

The illustration of the element mass distribution method.

Now by assembling the element matrixes, the global dynamic equilibrium equations for the FG cylinder can be obtained as

$$\mathit{M}\stackrel{··}{\mathit{q}}+\mathit{Kq}=\mathit{F}$$

(23)

where the matrix M and K represent the global mass matrix and the global stiffness matrix, respectively. F is the global vector of external force applied to the system. To solve the second-order ordinal-differential equation (23), the central difference method with suitable time step is applied in this article. In the case of zero initial conditions, that is, $\mathit{q}=0$ and $\stackrel{·}{\mathit{q}}=0$ at t = 0, the iterative equation is

$$\frac{1}{\mathrm{\Delta }{t}^2}\mathit{M}{\mathit{q}}_{\mathit{t}+\mathbf{\Delta }\mathit{t}}=\mathit{F}-(\mathit{K}-\frac{2}{\mathrm{\Delta }{t}^2}\mathit{M}){\mathit{q}}_t-\frac{1}{\mathrm{\Delta }{t}^2}\mathit{M}{\mathit{q}}_{\mathit{t}-\mathbf{\Delta }\mathit{t}}$$

(24)

where $\mathrm{\Delta }t$ denotes the time step of integration. The algorithm is only conditionally stable and the stable condition is

$$\mathrm{\Delta }t\le \mathrm{\Delta }{t}_{cr}=\frac{T_n}{\pi }$$

(25)

where $\mathrm{\Delta }{t}_{cr}$ denotes the critical stable time step, and $T_n$ represents the smallest period in the finite element assemblage.

Model validation

This section is to validate the proposed axisymmetric time-domain spectral element model. Transient response in a FGM cylinder calculated by the conventional finite element method by Hosseini et al. ¹⁰ is used as a model problem. As shown in Figure 5, a thick hollow cylinder with material gradually changed along the radial direction is subjected to dynamic inner pressure. The inner radius of the cylinder is $r_i=0.25\mathrm{m}$ and the outer radius is $r_o=0.50\mathrm{m}$ . The height of the cylinder is h = 0.01 m. The pressure loading function is defined as

$$P\left(t\right)=\{\begin{array}{cc}\hfill 4\times t\mathrm{GPa},\hfill & t\le 0.005\mathrm{s}\\ \multicolumn{1}{c}{0,}& t>0.005\mathrm{s}\end{array}$$

(26)

Figure 5.

The functionally graded hollow cylinder subjected to inner pressure.

It can be seen from equation (26) that the pressure increases until t = 0.005 s and transient wave propagation will occur after the end of the loading. In this example, the inner and outer wall of the cylinder is pure ceramic and pure metal, respectively. The properties of each material used in the FGM cylinder are listed in Table 1.

Table 1.

Material properties of ceramic and metal phase.

Material	Material properties
	Young’s modulus E (GPa)	Density $\rho$ (kg/m3)	Poisson’s ratio $\nu$
Ceramic	380	3800	0.3
Metal	70	2707	0.3

In Hosseini et al., ¹⁰ the FGM cylinder is assumed to be made of many sub-cylinders. Each sub-cylinder is considered as an isotropic layer. The material properties within a layer vary linearly in the layer thickness direction. The FGM cylinder is assumed to be in plane strain condition and is subjected to axisymmetric dynamic loading. Therefore, only 1D wave responses were solved by the conventional finite element method.

The time history and the frequency spectrum of the excitation loading are given in Figure 6. It can be seen that the excitation load has a frequency content of around 3 kHz. This frequency content is low. Generally, for the finite element method, about 20 elements are required per wavelength of the highest effective.^25,26 In this validation example, the highest effective frequency results in the maximum element size of about 0.080 m. In this article, the cylinder is meshed by the proposed spectral finite elements with $1\times 20$ elements. Each element is a fourth-order spectral element with $5\times 5$ nodes and the element size is $0.01\mathrm{m}\times 0.0125\mathrm{m}$ .The time step of integration is chosen as $1\times {10}^{-8}\mathrm{s}$ . Therefore, the mesh used here has the capability to simulate the dynamic behaviors.

Figure 6.

The time history (a) and the frequency spectrum (b) of the excitation loading.

Figure 7 shows the comparison of the time histories of the radial displacement component at the middle of the thickness for n = 0.5 obtained by the present method with those published data in Hosseini et al. ¹⁰ A very good agreement is observed. This validates the effectiveness of the proposed spectral element method. The differences between the current results with those published in Hosseini et al. ¹⁰ are mainly from the 2D axisymmetric analysis is performed in this paper, while only 1D wave propagation equation was solved in Hosseini et al. ¹⁰ The effect of the Poisson ratio can be considered by the proposed method. Besides, the present high-order spectral element can give better approximation of the varying material properties than low-order linear finite element. Therefore, it is believed that the proposed axisymmetric time-domain spectral element model can give better results for the analysis of elastic wave propagation in FG cylinders.

Figure 7.

The comparison of the time histories of the radial displacement at the middle point.

Numerical example and discussions

Model description

As shown in Figure 8, a FG cylinder with thickness L = 0.05 m and the radius r = 0.036 m is considered. Its top surface is pure ceramic (Al₂O₃) and the bottom is pure metal (Ni). The properties of two material phases are listed in Table 2.

Figure 8.

Functionally graded cylinder under the axial impact.

Table 2.

Material properties of top and bottom surfaces.

Material	Material properties
	Young’s modulus E ( ${10}^9\mathrm{Pa}$ )	Density $\rho$ (kg/m3)	Poisson’s ratio $\nu$
Al2O3	310	4400	0.24
Ni	210	8800	0.31

To study the elastic wave propagation in the FGM cylinder, an impulsive loading is applied to the region with radius $r_p=0.006\mathrm{m}$ at the center of the top surface and the pressure loading is defined as

$$f(r,t)=-\psi \left(t\right)P_0(H\left(r\right)-H(r-r_p))$$

(27)

where $P_0=1\times {10}^6\mathrm{Pa}$ , and H represents the Heaviside step function. The time-varying function $\psi (t)$ is a triangular pulse function with the peak load becomes at t₁ = 0.5 µs and lose its effect at t₂ = 1.0 µs. Although the loading and unloading duration is 1 µs, the problems were solved until 20 µs aimed to obtain the response of the cylinder for a sufficiently long time interval.

To investigate wave propagation in the FG cylinder with a reasonable accuracy and computational cost, convergence study was carried out. It was observed that using sixth-order spectral element and meshing the structure with 6 elements in the radial direction and 16 elements in the axial direction can give results with a reasonable accuracy and acceptable time and source costs. The time step of integration is chosen as $1\times {10}^{-8}\mathrm{s}$ to satisfy the stable condition. To study the effect of the compositional gradient exponents on wave propagation behavior, the analysis was repeated for the metal-rich (n = 10), balanced (n = 1), and ceramic-rich (n = 0.1) compositions, respectively.

As shown in Figure 9, three evaluation points are chosen to output displacement and stress time-history results. Point K is close to the loading area and near the ceramic surface (z = 0.05 m). Point L is located at the middle of the cylinder, and Point M is near the metal surface (z = 0 m).

Figure 9.

The coordinates of the chosen observation points and boundary points.

Results and discussions

Figures 10 and 11 show the distribution of the axial displacement component $w(r,z)$ and the radial displacement component $u(r,z)$ for the compositional gradient exponent n = 0.1, 1 and 10 at $t=5\mathrm{\mu }s$ , $t=7.5\mathrm{\mu }s$ , and $t=10\mathrm{\mu }s$ , respectively. It can be seen that after the load is applied to the top surface of the FG cylinder, the longitudinal and transverse waves are produced in the cylinder. The longitudinal waves travel faster than the transverse waves. Multiple reflections are observed when the waves travel to the geometry boundaries of the cylinder.

Figure 10.

Distribution of the axial displacement component $w(r,z)$ (m): (a) t = 5 µs, (b) t = 7.5 µs, and (c) t = 10 µs.

Figure 11.

Distribution of the radial displacement u(r, z) (m): (a) t = 5 µs, (b) t = 7.5 µs, and (c) t = 10 µs.

In Figure 10(a), the distributions of the axial displacement component $w(r,z)$ at $t=5\mathrm{\mu }s$ for n = 0.1, n = 1, and n = 10 are given. At this time, the maximum of the axial displacement occurs at DC edge in the negative direction. The profiles of the waves for n = 0.1, n = 1, and n = 10 are almost similar. However, it is shown that for three compositional gradient components, the waves travel with different speeds. The waves travel the fastest for the ceramic-rich composition (n = 0.1) and the slowest for the metal-rich composition (n = 10). Differences in axial displacement level can also be observed. The axial displacements for the ceramic-rich composition (n = 0.1) and balanced (n = 1) exhibit lower levels in comparison with those for the metal-rich composition (n = 10). The distributions of the axial displacement component at $t=7.5\mathrm{\mu }s$ for n = 0.1, n = 1, and n = 10 are plotted in Figure 10(b). It can be seen that the longitudinal waves in the cylinder with the ceramic-rich composition (n = 0.1) have traveled to the bottom surface and the reflections took place. Since the speed of transverse wave is lower than the speed of longitudinal wave, the distance between two wave fronts increases with the time. Figure 10(c) shows the axial displacement distributions at $t=10\mathrm{\mu }s$ . At this time, for the cylinder with the ceramic-rich composition (n = 0.1), the axial displacements are concentrated near the bottom right corner due to reflections. The maximum displacements appear in the negative direction near BC edge. For the cylinder with the metal-rich and balanced composition, the longitudinal waves have reached to the bottom surface and the reflections take place.

Figure 11(a) shows the distribution of the radial displacement component u(r, z) for n = 0.1, n = 1, and n = 10 at $t=5\mathrm{\mu }s$ . It can be seen that the radial displacements are almost positive, which means the structure in the state of radial expansion. Due to axisymmetry, the radial displacement component at the symmetry axis DA is zero. A wave front toward the BC edge can be found. The wave speeds for different material compositions are different. The displacement magnitudes increase with changing material composition from ceramic-rich (n = 0.1) to metal-rich (n = 10). This is due to the decreasing stiffness of the functionally grade cylinder with increasing through-thickness ceramic phase. Figure 11(b) shows the distributions of the radial displacement component at $t=7.5\mathrm{\mu }s$ . Since the wave speed in the cylinder with ceramic-rich composition (n = 0.1) is the fastest, wave reflections on the AB edge can be found. It is noted that the maximum displacement appears in the negative direction near the CD edge for the cylinder with metal-rich composition (n = 10). The radial displacement distributions at $t=10\mathrm{\mu }s$ are shown in Figure 11(c). Multiple reflections on AB and BC edge can be observed.

The distributions of the radial stress component ${\sigma }_{rr}$ are given in Figure 12 for the material gradient exponents n = 0.1, n = 1, and n = 10 at $t=5\mathrm{\mu }s$ , $t=7.5\mathrm{\mu }s$ , and $t=10\mathrm{\mu }s$ , respectively. Figure 12(a) shows the radial stress distributions in the cylinder for n = 0.1, n = 1, and n = 10 at $t=5\mathrm{\mu }s$ . It can be seen that the compositional gradient exponent n does not affect obviously the profile of radial stress and their critical locations. However, the maximum of the axial stress in the cylinder with the metal-rich composition and balanced composition are larger than those in the cylinder with the ceramic-rich composition. It can be found that the peak axial stress value increases as the material composition changed from the ceramic-rich to metal-rich. The distributions of axial stress at $t=7.5\mathrm{\mu }s$ are shown in Figure 12(b). For the cylinder with ceramic-rich composition, the radial stress values decrease due to reflections. At this time, the waves have not traveled to the bottom surface for the cylinder with balanced composition and metal-rich composition. Figure 12(c) shows the distributions of axial stress at $t=10\mathrm{\mu }s$ . Some stress concentrated locations can be found due to multiple wave reflections.

Figure 12.

Distribution of the radial stress ${\sigma }_{rr}$ (r, z) for n = 0.1, n = 1, and n = 10 (Pa): (a) t = 5 µs, (b) t = 7.5 µs, and (c) t = 10 µs.

Figure 13 shows the distributions of the axial stress component ${\sigma }_{zz}$ at $t=5\mathrm{\mu }s$ , $t=7.5\mathrm{\mu }s$ , and $t=10\mathrm{\mu }s$ for the material gradient exponents n = 0.1, n = 1, and n = 10, respectively. It can be found that the values of axial stress component are nearly twice higher than the other stress components. Figure 13(a) shows the distribution of the axial stress in the cylinder at $t=5\mathrm{\mu }s$ for n = 0.1, 1, and 10. It is noted that the stress profiles are almost same for n = 0.1, 1, and 10. After the loading is applied to the top surface, two stress wave packets are produced in the cylinder. The faster is longitudinal wave, which induces compressional deformation in the structure. This can be identified by two neighboring stress regions, where one is in compression and the other is in tension. Due to the difference in the material composition, the maximum stress value is different at this time. The maximum axial stress value is found in the FG cylinder for the metal-rich composition. The distributions of the axial stress in the cylinder at $t=7.5\mathrm{\mu }s$ for n = 0.1, 1, and 10 are plotted in Figure 13(b). At this time, for the case of n = 0.1, the longitudinal waves have been reflected by the bottom surface. Because free-surface boundary condition is considered here, it can be found that the longitudinal compressional wave is changed to the tension wave through the boundary reflection. In the case of n = 1, it is just at the beginning of the reflection. It is seen that the maximum axial stress is in tension near the bottom left corner. The distribution of the axial stress in the cylinder at $t=10\mathrm{\mu }s$ for n = 0.1, 1, and 10 is plotted in Figure 13(c). It can be observed that most areas of the cylinder are in the state of tension because of reflections.

Figure 13.

Distribution of the axial stress ${\sigma }_{zz}$ (r, z) at different times for n = 0.1, n = 1, and n = 10 (Pa): (a) t = 5 µs, (b) t = 7.5 µs, and (c) t = 10 µs.

Figure 14 shows the distributions of the shear stress component ${\sigma }_{rz}$ at $t=5\mathrm{\mu }s$ , $t=7.5\mathrm{\mu }s$ , and $t=10\mathrm{\mu }s$ for the material gradient exponents n = 0.1, n = 1, and n = 10, respectively. The levels of the shear stress are lower in comparison with other stress components. At the time of $t=5\mathrm{\mu }s$ , it can be found that the maximum shear stress is associated with the transverse wave induced by the applied impulsive loading. It is due to the shear deformation produced by the transverse wave. It can also be validated that the transverse waves have lower speed than the longitudinal waves.

Figure 14.

Distribution of the shear stress ${\sigma }_{rz}$ (r, z) at different times for n = 0.1, n = 1, and n = 10 (Pa): (a) t = 5 µs, (b) t = 7.5 µs, and (c) t = 10 µs.

Figure 15 shows the distributions of the hoop stress component ${\sigma }_{\theta \theta }$ at $t=5\mathrm{\mu }s$ , $t=7.5\mathrm{\mu }s$ , and $t=10\mathrm{\mu }s$ for the material gradient exponents n = 0.1, n = 1, and n = 10, respectively. In the FG cylinder, the hoop stress ${\sigma }_{\theta \theta }$ exhibits similar distribution but at lower levels compared to the axial stress ${\sigma }_{zz}$ . The values of the hoop stress increase when the material composition changes from ceramic-rich (n = 0.1) to metal-rich (n = 10). Figure 15(b) shows the hoop stress distribution at $t=7.5\mathrm{\mu }s$ . It can be found that the stress values decrease for the gradient exponent n = 0.1 when reflections occur. The distributions of hoop stress at $t=10\mathrm{\mu }s$ are shown in Figure 15(c). Multiple reflections can be observed.

Figure 15.

Distribution of the hoop stress ${\sigma }_{\theta \theta }$ (r, z) at different times for n = 0.1, n = 1, and n = 10 (Pa): (a) t = 5 µs, (b) t = 7.5 µs, and (c) t = 10 µs.

The time histories of the radial u and axial w displacement components at the evaluation points K, L, and M are illustrated in Figures 16 and 17. The time history of the radial displacement at point K is plotted in Figure 16(a). After the loading is applied, the wave reaches to the point K at about $t=1.1\mathrm{\mu }s$ and peaks at $t=1.6\mathrm{\mu }s$ . The radial displacement component at point K approximately equals to zero when the reflected waves not arrived at the point K from $t=5\mathrm{\mu }s$ . It can be seen that the wave speed for n = 10 is slower than it for n = 0.1 and n = 1. Figure 16(b) shows the time histories of the radial displacement at point L. The loading starts an effect at the point L at about $t=3\mathrm{\mu }s$ . The radial displacement becomes maximum in negative direction at about $t=7\mathrm{\mu }s$ for n = 0.1 and n = 1, whereas the peak appears at about $t=10\mathrm{\mu }s$ for n = 10 due to different wave speeds. Figure 16(c) illustrates the time histories of the radial displacement at point M. Although the profiles of radial displacement are similar, the difference of wave speed caused by different compositional gradient exponents is obvious. Figure 17 shows the time histories of the axial displacement at the reference points (K, L, and M) for three gradient exponents. Obviously, the axial displacement is greater than radial displacement. Figure 17(b) shows that the loading gets to affect the point L in the z-direction at about $t=3\mathrm{\mu }s$ for n = 0.1 and n = 1. In the case of n = 10, the wave speed is distinctly slower than other conditions. Overall, it can be found that the compositional gradient exponents n do not have the distinct effect on the displacement profiles, but work on the maximum displacement level and wave propagation speed.

Figure 16.

Time histories of the radial displacement at points K, L, and M: (a) K(0.006,0.04), (b) L(0.018,0.025), and (c) M(0.03,0.01).

Figure 17.

Time histories of the axial displacement at points K, L, and M: (a) K(0.006,0.04), (b) L(0.018,0.025), and (c) M(0.03,0.01).

The time histories of ${\sigma }_{rr},{\sigma }_{zz}$ , ${\sigma }_{rz}$ , and ${\sigma }_{\theta \theta }$ stress components at the evaluation points K, L, and M are plotted in Figures 18–20. Figure 18 shows the stresses components at point K. After the loading applied, the radial stress ${\sigma }_{rr}$ appears at $t=1.1\mathrm{\mu }s$ and peak at $t=1.6\mathrm{\mu }s$ . The stress value equals about zero during $t=5\mathrm{\mu }s$ to $t=7\mathrm{\mu }s$ because the reflected waves have not reached the point. It is obvious that the curves of stress components oscillate with low amplitude after $t=10\mathrm{\mu }s$ due to reflections. Since point K is close to the loading area, the peak stresses at point K are higher than those at points L and M. Figure 19 shows the time histories of the stress components at point L for n = 0.1, n = 1, and n = 10. It can be found that the wave speed is the fastest for n = 0.1 compared to the other two cases. The time histories of the stress components at the point M are illustrated in Figure 20. Because the point M is far from the loading surface, the arrival time differences between the longitudinal wave packet and the transverse wave packet are obvious. For n = 10 (metal-rich), the longitudinal wave speed can be calculated by the distance and the arrival time of the axial stress wave packets. It is estimated as c_L = 5412.7 m/s, which is very close to the theoretical longitudinal wave speed for pure metal (Ni), $c_{L,m}=5751.3\mathrm{m}/\mathrm{s}$ . The transverse wave speed can be calculated by the distance and the arrival time of the shear stress wave packets. It is estimated as $c_S=3752.3\mathrm{m}/\mathrm{s}$ , which is close to the theoretical transverse wave speed for pure metal (Ni), $c_{S,m}=3018.0\mathrm{m}/\mathrm{s}$ .

Figure 18.

Time histories of the stress components at point K.

Figure 19.

Time histories of the stress components at point L.

Figure 20.

Time histories of the stress components at point M.

Conclusion

The axisymmetric solid time-domain spectral element method is proposed to study the elastic wave propagation in FG cylinders. The cylinder’s material is graded either in the radial direction or in the thickness direction. The time-domain SFEM based on the high-order interpolation functions and the GLL node collocation and integration is attractive because it can avoid Runge’s phenomenon and needs lower computational costs compared to the conventional finite element method. In addition, using high-order interpolation function can give better approximation of gradually varying material properties. The validation of the proposed method is carried out by comparing with published data.

A FG cylinder with material graded in the thickness direction subjected to an impulsive shock loading is studied detailed by the present method. The displacement and stress distributions and the time histories of displacement and stress components are obtained for various grading patterns. The results show that the compositional gradient exponent n does not have a distinct effect on the displacement or stress profiles but has an obvious influence on the wave traveling speed as well as the displacement and stress levels. The proposed method can be used to other types of axisymmetric FG structures like thick truncated cones and hollow cylinders, and also various loading and boundary conditions can be considered to these problems.