## Abstract

The transient, coupled thermo-piezoelectric response of a functionally graded, radially polarized hollow cylinder under dynamic axisymmetric loadings was investigated in the present paper. To take into account the simultaneous coupling of displacement, temperature and electric fields as well as non-Fourier heat conduction effect, the Chandrasekharaiah theory of generalized thermo-piezoelectricity was employed. Except thermal relaxation time which was taken to be constant, profiles of all other material properties follow a volume-fraction-based rule with different non-homogeneity indices for each property. To solve three governing coupled partial differential equations, the Galerkin finite-element method was used in the Laplace domain. To restore time, a numerical scheme was employed for the Laplace inversion. When the cylinder was exposed to a highly transient thermal loading, effects of the non-homogeneity index and thermal relaxation time on the results were investigated.

## 1. Introduction

Functionally graded materials (FGMs) are composite media that have continuously changing material properties. The smooth variation of properties within FGMs results in lower stress concentration and intensity factors, higher fracture toughness and improved residual stress distribution compared with traditional laminated composites. FGMs have been widely used as thermal shields, wear-resistant linings, heat exchanger tubes, heat engine components and even prostheses (Birman & Byrd 2007).

In contrast, piezoelectrics are media in which elastic and electric fields are coupled, i.e. when a piezoelectric is deformed, a voltage is generated and vice versa. When a piezoelectric device operates at an elevated temperature, coupling of the temperature with the elastic and electric fields has also to be taken into account (Parton & Kudryavtsev 1988). Piezoelectrics have found numerous applications in the smart materials industry as sensors and actuators (Prasad *et al.* 2005).

If a piezoelectric material can be tailored in a way that its properties vary continuously in a certain direction, it can be used to improve the performance of smart structures such as surface acoustic wave sensors (Liu *et al.* 2003*a*). With advancements in the metallurgical science, it is possible to fabricate such non-homogeneous structures. They are called functionally graded piezoelectric materials (FGPMs). For instance, Takagi *et al.* (2002) fabricated a functionally graded piezoelectric bimorph actuator with an optimal material property profile leading to very high deflections with low stress values within the actuator.

To mathematically model the behaviour of piezoelectric materials exposed to different loading conditions, a proper thermo-piezoelectricity theory has to be chosen. Research has shown that the classical, uncoupled theories of thermo-piezoelectricity were unable to explain some physical observations such as the wave-like behaviour of temperature and the effect of strain and electric field rates on the temperature distribution (Hetnarski & Eslami 2008). In response, the generalized theories of thermo-piezoelectricity were proposed to account for the aforementioned phenomena. For instance, Chandrasekharaiah (1988) extended the generalized Lord & Shulman (1967) theory of thermoelasticity to a thermo-piezoelectricity theory.

There is a lack of literature on the generalized, coupled thermo-piezoelectric problems. A handful of papers dealt with one-dimensional problems owing to their relative simplicity. El-Karamany & Ezzat (2005) solved the problem of a semi-infinite, one-dimensional piezoelectric medium using a free energy function method and a state-space approach in the Laplace domain. They compared their results using different classical and generalized thermo-piezoelectricity theories. Majhi (1995) solved the problem of a semi-infinite piezoelectric rod heated over a portion of its length using the Chandrasekharaiah theory. He *et al.* (2002) used the state-space method to solve the problem of a semi-infinite rod subjected to a sudden temperature change at the rod end. More recently, Babaei & Chen (2009) investigated a one-dimensional problem of a finite thermo-piezoelectric material excited by a moving heat source using the Chandrasekharaiah theory.

To the authors’ best knowledge, there is no paper in the literature investigating the behaviour of FGPM structures using the generalized, coupled thermo-piezoelectricity theories. In the present paper, the Chandrasekharaiah theory is employed to investigate the transient response of an FGPM cylinder. The cylinder is polarized in the radial direction and subjected to an axisymmetric thermal shock. Except the thermal relaxation time, which is taken to be constant, all other material properties vary in the radial direction. The profile of these material properties is chosen on the basis of the volume fraction rule with different non-homogeneity indices for each one of the properties. This is known to be one of the most general forms of the property profile considered for FGMs. The Galerkin finite-element method with linear elements and test functions is employed to solve the governing, linear differential equations in the Laplace domain. The solutions for the displacement, electric and temperature fields are then inverted to the time domain by a numerical method. Finally, effects of the non-homogeneity indices and thermal relaxation time on the results are analysed and shown graphically.

## 2. Problem statement

As shown in figure 1, a long, hollow, FGPM cylinder is considered with its axis lying along the *z*-direction. The cylinder is poled and graded in the radial direction. The inner and outer radii of the cylinder are *a* and *b*, respectively. The temperature change *θ*, electric potential *Φ* and pressure *P* are specified on the surfaces of the cylinder, as shown in figure 1.

The profile of material properties follows the volume-fraction-based gradient rule as follows:
2.1
where *χ* is any property of the cylinder except the thermal relaxation time; the subscripts *a* and *b* denote the value of the property at the inner and outer surfaces, respectively, and *n*_{χ} is the non-homogeneity index of the corresponding material property, *χ*.

Figure 2 shows how a material property varies in the radial direction with different non-homogeneity indices. When *n*_{χ}=0, property *χ* is constant at the *χ*_{b} value, as shown in figure 2. It is worth noting that the dependence of material properties on temperature has not been sought in the present paper.

## 3. Fundamental relations and governing equations

The constitutive relations for a linear piezoelectric material are (Chandrasekharaiah 1988):
3.1
in which *σ*_{ij}, *S*, *D*_{i} and *q*_{i} are the stress, entropy, electric displacement and heat flux; *C*_{ijkl}, *e*_{ijk}, *β*_{ij}, *ρ*, *c* and *γ*_{i}, are, respectively, the elastic and piezoelectric coefficients, thermal moduli, density, specific heat and pyroelectric coefficients; ∈_{ij}, *τ* and *K*_{ij} are the dielectric coefficients, thermal relaxation time and coefficient of thermal conductivity; a subscript comma denotes the partial differentiation with respect to what follows it; *θ* is the temperature change with respect to the initial temperature, *T*_{0}, i.e. *θ*=*T*−*T*_{0}; *ε*_{kl} and *E*_{k} represent the strain and electric fields.

In the absence of electric current, free charge and body force, the equations of energy and motion, as well as the Coulomb equation for the conservation of electric charge, are:
3.2
where *u*_{i} is the displacement component.

The linear strain–displacement relation is
3.3
and for a quasi-stationary electric field, the following relation holds:
3.4
in which *Φ* is the electric potential.

For the current axisymmetric, plane strain problem
3.5
where *θ* is the azimuthal direction and (.) can be any variable.

Considering equation (3.5), the non-zero components of strain and electric field, equations (3.3) and (3.4), are
3.6
in which *u*=*u*_{r} for convenience.

So, the entropy and non-zero components of the stress, electric displacement and heat flux (equation (3.1)) are
3.7
where *c*_{αδ}=*c*_{ijkl} (*i*,*j*,*k*,*l*=1,2,3; *α*,*δ*=1,2,…,6), *e*_{mα}=*e*_{mij} (*m*=1,2,3; *α*=1,2,…,6) and *β*_{1}=*β*_{11},*β*_{3}=*β*_{33}.

Equation (3.2) is also reduced to the following set of equations:
3.8
We now introduce normalized parameters as follows:
3.9
where an overbar denotes a normalized value; is the speed of the elastic wave propagation in a homogeneous, linearly elastic isotropic medium and *η**=*ρ*_{a}*c*_{a}/*K*_{11a} is the reciprocal of the thermal diffusivity of the inner surface of the cylinder.

By substituting equation (3.7) into equation (3.8), the normalized governing equations of the problem are as follows:
3.10
in which *υ*=*e*_{33a}*b**c***η**/∈_{33a} and ∂( )/∂*t*=( )_{,t}, and the overbar is also dropped for convenience.

## 4. Solution procedure

Equation (3.10) is a system of linear, partial differential equations with variable coefficients. To solve equation (3.10), we first conduct the Laplace transform to eliminate time and then solve the governing equations in the Laplace domain by employing the Galerkin finite-element method. This solution will be eventually transformed to the time domain by a numerical Laplace inversion approach. Details of the solution procedure are presented in the following sections.

### (a) The Laplace transform technique

The Laplace transform and its inversion are defined as follows:
4.1
where *α* is an arbitrary real number greater than all real parts of singularities of ; *s* is the Laplace variable; a tilda denotes a transformed function and *i* is the imaginary unit.

After applying the Laplace transform to equation (3.10) and considering zero initial conditions, the governing equations in the Laplace domain read:
4.2
where *ι*=(1+*τ**s*)*s*. In equation (4.2), we have dropped tildas for convenience.

### (b) The Galerkin finite-element method

We now employ finite-element method to solve equation (4.2). Linear elements and shape functions are used to discretize the governing equations. The shape of the elements, local coordinates and shape functions are shown in figure 3.

The unknowns of the problem are then approximated over the elements as follows: 4.3 Now, we apply the Galerkin method (Eslami 2003) to equation (4.2): 4.4 To lower the order of the second-order derivatives and to extract mixed boundary conditions, which are the stresses on the boundaries, we then perform integration by part to all of the second-order derivatives in equation (4.4). After the manipulations and integrations, we find the element matrices. The matrices are then assembled to construct a system of linear algebraic equations for the unknowns at each node. Then all the unknowns can be readily obtained.

### (c) Numerical inversion of the Laplace transform

To find the final solution, the finite-element results of the previous section have to be transformed to the time domain. To do this, we use the numerical method proposed by Durbin (1974):
4.5
in which
4.6
where Δ*t* is the time increment and *T*_{period} is the time period over which the inversion is performed. Re and Im denote the real and imaginary parts of their arguments, respectively; *L* and *N* are the arbitrary parameters that influence the accuracy of the solution. To minimize both the discretization and truncation errors, it is recommended to observe the following constraint (Durbin 1974):
4.7
In the present work, the following values were chosen for the above parameters: *α**T*=5 and *N**L*=4500.

## 5. Results

In this section, we analyse the effects of both the non-homogeneity index and thermal relaxation time on the response of the cylinder. All boundary conditions shown in figure 1 are taken to be zero except the non-dimensional temperature at the inner surface of the cylinder, which is *θ*_{a}=10^{4}*t*e^{−1000t} (*t* is non-dimensional). The non-dimensional values of the inner and outer radii of the cylinder are taken to be 0.01 and 0.02, respectively.

The inner and outer surfaces of the cylinder are made of lead zirconate titanate and cadmium selenide, whose properties are listed in table 1.

Figure 4*a*–*e* depicts the effect of the non-homogeneity index on the distribution of the displacement, electric potential, temperature change, and radial and hoop stresses, respectively, with *t*=0.0011 and *τ*=0.05. In this section, the non-homogeneity indices are taken to be the same for all varying material properties, i.e. *n*_{χ}=*n*. Three different values of the non-homogeneity index were used in the calculation, namely, *n*=0 (cadmium selenide), *n*=1 (linear variation) and *n*≫1 (lead zirconate titanate).

In figure 4, a non-zero thermal relaxation time results in non-Fourier heat conduction and consequently the appearance of thermal wavefronts in the solutions. These wavefronts travelling from inner surface to the outer surface in figure 4 serve as a messenger of the thermal disturbance along the thickness: in front of the wavefront, the solutions remain at their initial values. However, there is no wavefront in the electric potential distribution, as shown in figure 4*b*.

As shown in figure 4*a*,*c*–*e*, an increase in *n* leads to a contraction of the undisturbed portion of the results in the cylinder. Since the velocity of the thermal wave in a hyperbolic heat conduction problem is equal to , this contraction may be attributed to a higher average of this term when *n* is increased. It is worth noting that the contraction of the displacement, when *n* is increased from 1 to *n*≫1, is not significant, as depicted in figure 4*a*.

Figure 5*a*–*f* shows the effects of the non-homogeneity index on the history of the solutions at the middle point of the cylinder.

When *n*=0, the amplitudes of the dynamic displacement, stresses and electric displacement are the smallest of the results of all three *n* values considered, as shown in figure 5*a*,*d*–*f*. However, the largest amplitude of the electric potential occurs when *n*=0, at the smaller time values, approximately *t*<0.1, as shown in figure 5*b*. The smallest amplitude of the electric potential happens when *n*≫1, as shown in figure 5*b*. In addition, the amplitude of the electric displacement increases monotonically with increasing *n*, as shown in figure 5*f*.

When *n*=1, both the dynamic displacement and the temperature change exhibit the largest amplitudes of all three *n* values considered, as depicted in figure 5*a*,*c*, respectively.

Although not shown, it is worth noting that when *n*=2,3,…, the results are almost the same as those for *n*=1.

Figure 6*a*–*e* depicts the effects of thermal relaxation time, *τ*, on the results when *n*=0.

When *τ*=0, there does not exist any wavefront in the cylinder, i.e. the disturbance spreads throughout the cylinder right after the inner side is exposed to the thermal shock and there is no undisturbed portion in the cylinder, as shown in figure 6*a*,*c*–*e*. Moreover, a larger thermal relaxation time results in a larger undisturbed portion over which the field variables remain at their initial values, as shown in figure 6*a*,*c*–*e*.

When *τ*≠0, the electric potential distributes almost linearly in front of the wavefront of the temperature, as shown in figure 6*b*.

Figure 7*a*–*f* depicts the effects of the thermal relaxation time on the history of the field variables. It is seen that when *τ*=0, the amplitudes of all the variables reach their minimum. The amplitude of the displacement increases monotonically as the thermal relaxation time increases, as shown in figure 7*a*. Furthermore, by increasing the thermal relaxation time, the time needed for the temperature to become zero is also increased, as illustrated in figure 7*c*.

At the very beginning of the thermal disturbance (approx. for *t*<0.05), *τ*=0.05 results in higher amplitudes of the electric potential, stresses and electric displacement; however, when *τ*=0.5, the amplitudes of the results increase more quickly and overtake others at about *t*=0.05.

Furthermore, a homogeneous cylinder of lead zirconate titanate is considered in the calculation. Both solutions exhibit almost exactly the same trend for the electric potential, which corroborates the correctness of the derivations and solution of the current paper.

It is worth noting that commercial multiphysics analysis software such as Comsol and Ansys may also be employed to numerically solve the present problem; however, additional user-written codes are needed; these codes have to introduce hyperbolic heat conduction theory and take into account non-homogeneity of the medium.

## 6. Conclusions

The generalized coupled thermo-piezoelectricity theory of Chandrasekharaiah has been employed to investigate the multiphysics behaviour of a radially polarized FGPM hollow cylinder exposed to a thermal shock. The governing thermo-piezoelectricity equations of the non-homogeneous cylinder have been developed and solved using the Galerkin finite-element method in the Laplace domain, and the results were inverted into the time domain numerically. The effects of the non-homogeneity index and thermal relaxation time on the results were also investigated and shown graphically. The conclusions drawn are as follows.

— A non-zero thermal relaxation time results in non-Fourier heat conduction and wavefronts in the dynamic distribution of the displacement, temperature and stresses. However, the electric potential does not exhibit any wave-like feature, and the thermal shock applied on the inner surface is felt by the electric field throughout the whole thickness of the cylinder immediately.

— An increase in the non-homogeneity index,

*n*, leads to a contraction of the undisturbed portion of the cylinder over which the field variables are still at their initial values.— Of all three values of the non-homogeneity index, namely,

*n*=0,*n*=1 and*n*≫1, the amplitudes of the displacement, stresses and the electric displacement reach their smallest values at*n*=0. In contrast, the amplitudes of the displacement and temperature change reach their largest values at*n*=1.— A higher value of the thermal relaxation time,

*τ*, results in a larger undisturbed region in the thickness direction of the cylinder. Moreover, when*τ*=0, all field variables exhibit minimal variations.

## Footnotes

- Received October 14, 2009.
- Accepted November 2, 2009.

- © 2009 The Royal Society