Static Analysis of Functionally Graded Piezoelectric Beams under Thermo-Electro-Mechanical Loads

Abstract

This paper presents the analysis of static bending of beams made of functionally graded piezoelectric materials (FGPMs) under a combined thermo-electro-mechanical load. The Euler Bernoulli theory (EBT), first-order shear deformation theory (FSDT) and third-order shear deformation theory (TSDT) were employed to compare the accuracy and the reliability of each theory in applications. The material properties vary continuously through the thickness direction. The material compositions were selected from the PZT family. The governing equations were derived from Hamilton's principle and solved using the finite element method and Fourier series method. Cubic Hermit interpolation shape function was used for estimating the transverse deflection, and the linear interpolation function was used for the axial displacement and the shear rotation as well. Fourier series expansion, based on the boundary conditions, were employed to solve the governing equations analytically. The accuracy of the method was validated by comparing the results with the previous studies. Finite element results were compared with the analytical results presented in this paper. A comprehensive parametric study is conducted to show the influence of the voltage, shear deformation, material composition, end supports, and the slenderness ratio on the thermo-electro-mechanical characteristic.

1. Introduction

The concept of the well-known functionally graded materials (FGMs) was extended to a piezoelectric material to improve its reliability and performance, and a new kind of material, such as a functionally graded piezoelectric material (FGPM), was developed [1]. Monomorph and bimorph piezoelectric actuators are used for applications which require large displacement output. Since these actuators are manufactured by bonding one or two piezoelectric layers to a metal shim layer, the bonding interface has low strength. On the other hand, large stress discontinuity and concentration are initiated in the bonding interface due to the difference of material properties between the shim and piezoelectric layer [2]. To solve this problem, a novel class of actuators were developed in 1995 in which one or more layers were made from FGPMs with thermo-electromechanical properties varying smoothly along the layer thickness [3]. The FGPMs are designed to withstand against high thermal disturbances specially in the cases with shock nature, and it has been proven, both theoretically and experimentally, that the durability, reliability, and actuation performance of actuators can be greatly improved by the use of FGPMs. A number of fabrication techniques for the FGPM structures have been proposed and successfully implemented [2, 4–8].

Many structural systems are routinely modelled as beams. Hence, the static behavior of rectangular FGPM beams in terms of the transverse deflection, transverse shear strain, and stresses may be examined. The design and analysis of these types of smart structures may require the use of the finite element method. Detailed reviews of the finite element (FE) modelling for adaptive structural elements consisting of piezoelectric hybrid plates have been reported by Benjeddou [9] and Mackerle [10]. Based on the layerwise theory, Saravanos et al. [11] applied a numerical investigation for the dynamic behavior of piezoelectric laminates. Kapuria and Alam [12] used a finite element model for dynamic analysis of hybrid piezoelectric beams based on the zigzag theory. Heiliger and Reddy [13] presented the bending and vibration problems for FGM structures based on the TSDT. In another study, Fernandes and Pouget [14] proposed a new approach for governing the equations of the laminated plates with piezoelectric layers by expressing the electric potential as a function of the thickness coordinate. The shear correction factor of elastic displacements was considered in their study. The modelling and optimal design method for the piezoelectric microactuators was developed by DeVoe and Pisano [15]. They assumed that the piezoelectric coefficient varies linearly along the thickness direction and other material parameters assumed to be constant. Hauke et al. [16] performed an experimental study, and they showed that, compared to the homogenous beam, the deflection of an FGPM beam is only slightly smaller when beam is long; on the other hand, the internal mechanical stress is drastically reduced. Due to these advanced properties, the mechanical stresses in beams with FGPM sensors and actuators are significantly decreased [17]. Closed-form solutions for the FGPM cantilevers with a linear gradient for the mass density and with a linear gradient for the piezoelectric constant e₃₁ were performed by Zhifei and Shi [18, 19]. He used the two-dimensional theory of elasticity to develop the governing equations and employed the Airy stress function to solve them. Based on the EBT, Kruusing [20] presented an analytical solution to calculate the tip deflection and the blocking force of a cantilever actuator with a graded elastic constant when the actuator was subjected to a transverse point load at the tip. Joshi et al. [21, 22] obtained exact solution for the electrically induced bending [21] and thermoelectric bending [22] of a laminated beam consisted of an FGPM layer and an active piezoelectric layer. Yang and Xiang [23] employed generalized differential quadrature (GDQ) method to solve the governing equations of the FGPM beam, and they validate the accuracy of GDQ method comparing to the other numerical methods.

The literature survey gives a clear indication that studies on the static deflection and stresses under ambient loading conditions of FGPM beams, using the TSDT, have been rarely reported in the literature. In the present paper, the principle of stationary potential energy is used to obtain the static finite element equilibrium equations of the FGPM beams with uniformly distributed mechanical, electrical, and thermal loads. The finite element formulation follows the conventional displacement-based approach employed to derive the stiffness and force matrices. In addition, another approach to obtain the governing equations of the FGPM beam based on the EBT is used in the paper. The Fourier series expansion is employed to solve the equations analytically.

2. Material Properties

Figure 1 shows the configuration of FGPM beam made of PZT-4 and PZT-5 with the length L, width b, and thickness h. The material properties vary continuously along the z-direction. Top surface consisted of pure PZT-4, and the bottom surface is composed of pure PZT-5. Different mathematical models are developed to simulate the distribution of the material constituent in the FGMs [24–26]. One frequently used model for the FGMs is the power law model proposed by Reddy [26], which has the privilege of leading to simplified equations to be solved for analytical solutions. In this formulation, the volume fraction of material composition are assumed to obey the power law distribution with the simple law of constituent mixture as

$\begin{array}{cc}\hfill P_e\left(z\right)& =\left(P_t-P_b\right){\left(\frac{z}{h}+0.5\right)}^n\hfill \\ \hfill & \hspace{1em}+P_b\left(-\frac{h}{2}\le z\le +\frac{h}{2},\hspace{0.17em}0\le n\le \infty \right),\hfill \end{array}$

(1)

where P _e is the effective material property of the FGPM beam, P _t is top surface (PZT-4) property, and P _b is the bottom surface (PZT-5) property. In this paper, z=+h/2 leads to P _e = P _t that denotes that the top surface of the FGPM beam consists of PZT-4 and z=-h/2 leads P _e = P _b which means that, at the bottom surface of the FGPM beam, PZT-5 is the only material constituent. In the power law formulation, n is the volume fraction index of FGPM beam. For n=0, we have P _e = P _t , and for n=∞, we have P _e = P _b . The two cases correspond to the isotropic PZT-4 and PZT-5, respectively. The material constants are presented in Table 1.

Table 1:

Material properties of FGPM constituent.

		PZT-4 Yang and Xiang [23]	PZT-5 Yang and Xiang [23]
Elastic modulus (GPa)	C 11	81.3	60.6
	C 55	25.6	23.0
Piezoelectric constant (C/m2)	e 31	−5.2	−6.5
	e 15	12.7	17
Dielectric constant (10− 11 F/m)	∊11	1300.9	2770
	∊33	1150	3009
Pyroelectric constant (10− 5 C/m2K)	P	5.48	2.5
Density (Kg/m3)	ρ	7500	7500
Thermal conductivity (W/Km)	λ	2.1	1.5
Thermal expansion (10− 61/K)	α	2	10
Specific heat (J/km)	c P	339.45	776.53

Figure 1:

The geometry of the FGPM beam and the possible variation of material constituent through the thickness.

3. Governing Equations

3.1. Finite Element Method

The geometry of deformation for a beam based on different beam theories is presented in Figure 2. The rotation, $\partial w/\partial z$ , can be expressed as the sum of bending rotation ψ _x (or rotation of the normal) and rotation due to shear γ. The displacement field based on the TSDT, as stated in Heiliger and Reddy [13], is written as

$\begin{array}{c}u\left(x,z\right)=u_0\left(x\right)+z{\psi }_x+c_0{z}^3\left[{\psi }_x+w_{0,x}\right],\\ w\left(x,z\right)=w_0\left(x\right),\end{array}$

(2)

where c₀ = −4h²/3. Subscript comma denotes the differentiation, u₀ is the midplane axial displacement, and w₀ is the midplane transverse deflection. This displacement field accounts for traction-free condition on the top and bottom surfaces of the beam. A Cartesian coordinate system (x,z) is introduced in the present analysis. We consider an FGPM beam with length L and depth h, as shown in Figure 1. The beam is polarized in z-direction and has a graded profile in this direction. The potential field is approximated linearly through the thickness, as given by

$\Phi \left(z,x\right)=\frac{2z}{h}\varphi \left(x\right).$

(3)

Figure 2:

Kinematics of deformation of a beam in various beam theories.

Let ϵ _{i j} ,σ _{i j} ,D _i and E _i denote the strain tensor, stress tensor, dielectric displacement, and the induced electric field, respectively. The constitutive equations for a piezoelectric beam can be written as [27]

$\begin{array}{cc}\hfill \left\{\begin{array}{c}{\sigma }_x\\ {\tau }_{xz}\end{array}\right\}& =\left[\begin{array}{cc}{C}_{11}& 0\\ 0& C_{55}\end{array}\right]\left\{\begin{array}{c}{ϵ}_x\\ {ϵ}_{xz}\end{array}\right\}-\left[\begin{array}{cc}0& e_{11}\\ e_{15}& 0\end{array}\right]\left\{\begin{array}{c}0\\ E_z\end{array}\right\}-\left\{\begin{array}{c}\beta \\ 0\end{array}\right\}T,\hfill \\ \hfill \left\{\begin{array}{c}{D}_x\\ D_z\end{array}\right\}& \begin{array}{c}=\left[\begin{array}{cc}0& e_{15}\\ e_{31}& 0\end{array}\right]\left\{\begin{array}{c}{ϵ}_x\\ {ϵ}_{xz}\end{array}\right\}-\left[\begin{array}{cc}{∊}_{11}& 0\\ 0& {∊}_{33}\end{array}\right]\left\{\begin{array}{c}0\\ E_z\end{array}\right\}-\left\{\begin{array}{c}{P}_1\\ P_2\end{array}\right\}T\end{array}\hfill \end{array}$

(4)

in which β is defined by

$\beta =C_{11}\alpha$

(5)

and C _{i j} ,e _{i j} ,P _i ,α as well as ∊ _{i j} are the components of the effective elastic, piezoelectric, pyroelectric, thermal expansion, and dielectric constants of the piezoelectric material, respectively. The straindisplacement and the electric field potential relations are

$\begin{array}{c}{ϵ}_x=u_{x,x},\hspace{1em}\hspace{1em}{ϵ}_z=w_{,z}=0,\\ {ϵ}_{xz}=\frac{1}{2}\left(u_{x,z}+u_{z,x}\right),\\ E_z=-{\Phi }_{,z}.\end{array}$

(6)

With these assumptions, the governing equations can be derived using Hamilton's principle as

$\delta {{\displaystyle \int }}_{t_0}^{t_1}\left(K-H\right)dt+\delta {{\displaystyle \int }}_{t_0}^{t_1}R\hspace{0.17em}dt=0,$

(7)

where the variation of electric enthalpy density δH, the variation of kinetic energy δK, and the virtual work δR due to the mechanical load q(X,t) are

$\begin{array}{cc}\hfill \delta H& =b{{\displaystyle \int }}_{t_0}^{t_1}{{\displaystyle \int }}_{-h/2}^{h/2}\left({\sigma }_x\delta {ϵ}_x+{\tau }_{xz}\delta {\gamma }_{xz}+\delta \chi -D_z\delta E_z\right)dz\hspace{0.17em}dx,\hfill \\ \hfill \delta K& =b{{\displaystyle \int }}_{t_0}^{t_1}{{\displaystyle \int }}_{-h/2}^{h/2}\left(\dot{U}\delta \dot{U}+\dot{W}\delta \dot{W}\right)dz\hspace{0.17em}dx,\hfill \\ \hfill \delta R& =b{{\displaystyle \int }}_v{f}_{bi}\delta u_{i}dv+b{{\displaystyle \int }}_s\left(f_{si}\delta u_i+q_e\varphi +q_t\theta \right)ds+f_{ci}\delta u_i.\hfill \end{array}$

(8)

Here, f _{b i} ,f _{c i} ,f _{s i} ,q _e and q _t represent the body forces, concentrated loads, the specified traction, electrical surface charge, and heat flux, respectively. Using (4) and substituting (8) into the Hamilton principle (7) and integrating over the beam thickness gives

$\begin{array}{cc}\hfill & \delta {{\displaystyle \int }}_{t_0}^{t_1}{{\displaystyle \int }}_v(-\rho {\ddot{u}}_i\delta u_i-C_{ij}{ϵ}_i\delta {ϵ}_j+e_{ij}{E}_j\delta {ϵ}_i+{\beta }_i\theta \delta {ϵ}_i\hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}+e_{ij}{ϵ}_j\delta {\varphi }_{,i}+{ϵ}_{ij}{E}_j\delta {\varphi }_{,i}-P_i\theta \delta {\varphi }_{,i}+{\beta }_i{ϵ}_i\delta \theta \hfill \\ \hfill & \hspace{1em}\hspace{1em}\hspace{1em}\hspace{0.17em}-P_i{E}_i\delta \theta +\alpha \theta \delta \theta \phantom{-\rho {\ddot{u}}_i\delta u_i-C_{ij}{ϵ}_i\delta {ϵ}_j+e_{ij}{E}_j\delta {ϵ}_i+{\beta }_i\theta \delta {ϵ}_i})dv\hspace{0.17em}dt\hfill \\ \hfill & \hspace{1em}\hspace{1em}+{{\displaystyle \int }}_{t_0}^{t_1}{{\displaystyle \int }}_v{f}_{bi}\delta u_{i}dv\hspace{0.17em}dt+{{\displaystyle \int }}_{t_0}^t{f}_{ci}\delta u_{i}dt\hfill \\ \hfill & \hspace{1em}\hspace{1em}+{{\displaystyle \int }}_{t_0}^t{{\displaystyle \int }}_s\left(f_{si}\delta u_i+q\delta \varphi +q_e\theta \right)ds\hspace{0.17em}dt=0.\hfill \end{array}$

(9)

The degrees of freedom of the base element (e) are $(u_0,w_0,\partial w_0/\partial x,{\psi }_x)$ . To meet the convergence requirements, the interpolation functions for $(u_0,w_0,\partial w_0/\partial x,{\psi }_x)$ must be continuous at the element boundaries. Hence, $(w_0,\partial w_0/\partial x)$ were expanded using the cubic Hermite interpolation in terms of the nodal values of (w₀,w_0,x); respectively, and a linear interpolation is used for (u₀,ψ _x ). Thus, at the element level, each node has four degrees of freedom (u₀,w₀,w_0,x,ψ _x ), for the displacements. For the static case, by substituting the strain-displacement matrix and element degrees of freedom in the strain energy expression, (9), and using the principle of stationary potential energy as well as the first variation of Π, that is, δΠ=0, the finite element dynamic equations are derived as

$\left[M_{uu}\right]\left\{\ddot{u}\right\}+\left[K_{uu}\right]\left\{u\right\}+\left[K_{u\varphi }\right]\left\{\varphi \right\}=\left\{F_m\right\}+\left\{F_{u\theta }\right\},$

(10)

$\left[K_{\varphi u}\right]\left\{u\right\}-\left[K_{\varphi \varphi }\right]\left\{\varphi \right\}=\left\{F_q\right\}-\left\{F_{\varphi \theta }\right\}.$

(11)

The elements of the mass matrix [M _{u u} ], membrane stiffness matrix [K _{u u} ], the bending stiffness matrix [K _{u ϕ} ], the piezoelectric stiffness matrices [K _{ϕ ϕ} ], and the force matrices [F _{u θ} ],[F _{ϕ θ} ],[F _q ^e ],[F _m ] are

$\begin{array}{c}\left[M_{uu}\right]={{\displaystyle \int }}_v{\left[N\right]}^T{\left[H\right]}^T\left[\rho \right]\left[H\right]\left[N\right]dv,\\ \begin{array}{c}\left[K_{uu}\right]={{\displaystyle \int }}_v\phantom{\frac{1}{2}}{\left[Bu\right]}^T{\left[C\right]}^T\left[Bu\right]dv\end{array},\\ \begin{array}{c}\left[K_{u\varphi }\right]={{\displaystyle \int }}_v\phantom{\frac{1}{2}}{\left[Bu\right]}^T{\left[e\right]}^T\left[B_{\varphi }\right]dv\end{array},\\ \begin{array}{c}\left[K_{\varphi \varphi }\right]={{\displaystyle \int }}_v\phantom{\frac{1}{2}}{\left[B_{\varphi }\right]}^T{\left[∊\right]}^T\left[B_{\varphi }\right]dv\end{array},\\ \begin{array}{c}\left[F_{u\theta }\right]={{\displaystyle \int }}_v\phantom{\frac{1}{2}}{\left[Bu\right]}^T{\left[\beta \right]}^T\theta dv\end{array},\\ \begin{array}{c}\left[F_{\varphi \theta }\right]={{\displaystyle \int }}_v\phantom{\frac{1}{2}}{\left[B_{\varphi }\right]}^T{\left[P\right]}^T\theta dv\end{array},\\ \begin{array}{c}\left[F_q^e\right]={{\displaystyle \int }}_s\phantom{\frac{1}{2}}{\left[N\right]}^T\left[q_e\right]ds\end{array},\\ \begin{array}{c}\left[F_m\right]={{\displaystyle \int }}_v\phantom{\frac{1}{2}}{\left[N\right]}^T{\left[H\right]}^T\left\{f_b\right\}dv+{{\displaystyle \int }}_s\phantom{\frac{1}{2}}{\left[N\right]}^T{\left[H\right]}^T\left\{f_s\right\}ds.\end{array}\end{array}$

(12)

When electric load is applied, (11) is omitted and $[K_{u\varphi }]\{\varphi \}$ is transferred to the right side of (10) and is treated as a force. For static analysis of the FGPM beam under electrical load, the following equation should be solved:

$\left[K_{uu}\right]\left\{u\right\}=-\left[K_{u\varphi }\right]\left\{\varphi \right\}+\left\{F_m\right\}+\left\{F_{u\theta }\right\}.$

(13)

In the present analysis, thermal load is considered as well. Constant surface temperature is imposed at the top and bottom surfaces. The variation in temperature is assumed to occur along the thickness direction only. Thus, the steady-state heat transfer equation is reduced to the one-dimensional equation as

$\frac{\partial }{\partial z}\left(K\left(z\right)\frac{\partial T}{\partial z}\right)=0$

(14)

with the boundary conditions such that T=T _t at z=h/2 and T=T _b at z=-h/2. Here, T is the temperature; T _t and T _b are the applied temperatures on the upper and lower surfaces, respectively. Term K is the thermal conductivity that may be assumed to follow the power law form, indicated by (1). The temperature field through the thickness of the beam is obtained by solving (14).

3.2. Fourier Series Solution

Consider an FGPM beam as shown in Figure 1. The material properties change gradually between the upper and lower surfaces according to (1). In the absence of the thermal load and temperature change, the constitutive equations for the stresses and the electric displacements for the steady-stated FGPM beam based on the EBT reduce to

$\begin{array}{cc}\hfill {\sigma }_x& =C_{11}{ϵ}_x+e_{11}\Phi \left(z,x\right),\hfill \\ \hfill D_z& =e_{31}{ϵ}_x+{∊}_{33}\Phi \left(z,x\right).\hfill \end{array}$

(15)

For mathematical simplification, the potential Φ(z,x) is assumed linearly across the FGPM beam thickness as [28]

$\Phi \left(z,x\right)=\frac{z}{h}\left({\varphi }^{+}\left(x\right)-{\varphi }^{-}\left(x\right)\right),$

(16)

in which ϕ⁺ and ϕ⁻ are the electric potentials on the upper and lower surfaces of the FGPM beam, respectively. For the static analysis of the beam, the virtual work done by the electromechanical internal forces in the FGPM beam is

$\delta H={{\displaystyle \int }}_v\left({\sigma }_x\delta {ϵ}_x-D_z\delta E_z\right)dv.$

(17)

Using (15) and the electric field relation of (16), and substituting into (17) and using the virtual work by carrying the variation formulation, the governing equations are obtained as

$\begin{array}{c}\begin{array}{c}\delta u_0:{\acute{A}}_1{u}_{0,xx}+{\acute{A}}_2{w}_{0,xxx}+{\acute{A}}_3{\varphi }_{,x}^{+}+{\acute{A}}_4{\varphi }_{,x}^{-}=0,\\ \delta w_0:{\acute{B}}_1{u}_{0,xxx}+{\acute{B}}_2{w}_{0,xxxx}+{\acute{B}}_3{\varphi }_{,xx}^{+}+{\acute{B}}_4{\varphi }_{,xx}^{-}={\acute{B}}_5,\\ \delta {\varphi }^{+}:{\acute{C}}_1{u}_{0,x}+{\acute{C}}_2{w}_{0,xx}+{\acute{C}}_3{\varphi }^{+}+{\acute{C}}_4{\varphi }^{-}=0,\\ \delta {\varphi }^{-}:{\acute{D}}_1{u}_{0,x}+{\acute{D}}_2{w}_{0,xx}+{\acute{D}}_3{\varphi }^{+}+{\acute{D}}_4{\varphi }^{-}=0,\end{array}\end{array}$

(18)

where $\acute{A},\acute{B},\acute{C},\acute{D}$ are given in the Appendix. To solve the governing equations simultaneously, dimensionless values are defined as

$\begin{array}{c}{\bar{u}}_0=\frac{u}{L},\hspace{1em}\hspace{1em}\bar{x}=\frac{x}{L},\hspace{1em}\hspace{1em}{\bar{w}}_0=\frac{w_0}{L},\hspace{1em}\hspace{1em}{\bar{\varphi }}^{+}=\frac{e_{31t}}{C_{11t}L}{\varphi }^{+},\\ {\bar{\varphi }}^{-}=\frac{e_{31t}}{C_{11t}L}{\varphi }^{+},\end{array}$

(19)

where C_11t and e_31t are Young's modulus and piezoelectric constant of the upper surface of the FGPM beam, respectively. The sign (⁻) indicates dimensionless value. Using the dimensionless parameters, the governing (18) is given as

$\begin{array}{c}{\acute{a}}_1{\bar{u}}_{0,\bar{x}\hspace{0.17em}\bar{x}}+{\acute{a}}_2{\bar{w}}_{0,\bar{x}\hspace{0.17em}\bar{x}\hspace{0.17em}\bar{x}}+{\acute{a}}_3{\varphi }_{,\bar{x}}^{+}+{\acute{a}}_4{\varphi }_{,\bar{x}}^{-}=0,\\ {\acute{b}}_1{u}_{0,\bar{x}\hspace{0.17em}\bar{x}\hspace{0.17em}\bar{x}}+{\acute{b}}_2{w}_{0,\bar{x}\hspace{0.17em}\bar{x}\hspace{0.17em}\bar{x}\hspace{0.17em}\bar{x}}+{\acute{b}}_3{\varphi }_{,\bar{x}\hspace{0.17em}\bar{x}}^{+}+{\acute{b}}_4{\varphi }_{,\bar{x}\hspace{0.17em}\bar{x}}^{-}={\acute{b}}_5,\\ {\acute{c}}_1{u}_{0,\bar{x}}+{\acute{c}}_2{w}_{0,\bar{x}\hspace{0.17em}\bar{x}}+{\acute{c}}_3{\varphi }^{+}+{\acute{c}}_4{\varphi }^{-}=0,\\ {\acute{d}}_1{u}_{0,\bar{x}}+{\acute{d}}_2{w}_{0,\bar{x}\hspace{0.17em}\bar{x}}+{\acute{d}}_3{\varphi }^{+}+{\acute{d}}_4{\varphi }^{-}=0,\end{array}$

(20)

where ${\acute{a}}_s,{\acute{b}}_s,{\acute{c}}_s$ , and ${\acute{d}}_s$ are dimensionless constants. In the present analysis, an analytical solution is obtained for the simply supported FGPM beam with the following boundary conditions as:

$\begin{array}{cc}\hfill {\sigma }_x=0\hspace{1em}\hspace{1em}⟶{\bar{u}}_{,\bar{x}}=0\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}& \bar{x}=\mathrm{0,1},\hfill \\ \hfill \bar{w}=0\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}& \bar{x}=\mathrm{0,1},\hfill \\ \hfill \bar{\varphi }=0\hspace{1em}⟶{\bar{\varphi }}^{+}={\bar{\varphi }}^{+}=0\hspace{1em}\hspace{1em}\hspace{1em}& \bar{x}=\mathrm{0,1}.\hfill \end{array}$

(21)

Considering the assumed boundary conditions, to solve the system of (20), a finite Fourier transformation can be used as

$\begin{array}{cc}\hfill {\bar{u}}_{0m}& ={{\displaystyle \int }}_0^1{\bar{u}}_0\left(\bar{x}\right)\cos \left(m\pi \bar{x}\right)d\overline{x,}\hfill \\ \hfill {\bar{w}}_{0m}& ={{\displaystyle \int }}_0^1{\bar{w}}_0\left(\bar{x}\right)\sin \left(m\pi \bar{x}\right)d\bar{x},\hfill \\ \hfill {\bar{\varphi }}_m^{+}& ={{\displaystyle \int }}_0^1{\bar{\varphi }}^{+}\left(\bar{x}\right)\sin \left(m\pi \bar{x}\right)d\bar{x},\hfill \\ \hfill {\bar{\varphi }}_m^{-}& ={{\displaystyle \int }}_0^1{\bar{\varphi }}^{-}\left(\bar{x}\right)\sin \left(m\pi \bar{x}\right)d\bar{x}.\hfill \end{array}$

(22)

Formulas for the inverse of transformation are obtained by using relationship from the theory of Fourier series. To obtain the mth Fourier components of unknown variables, the system of (20) must be solved based on the choice of m.

4. Results and Discussion

Consider an FGPM beam with top surface made of monomorph PZT-4 and the bottom surface of PZT-5. The dimensionless deflection is defined as w₀/h.

The present analysis is validated in Table 2 and Figure 3 through comparisons with the existing solutions available in the literatures. The tip deflection of a standard cantilever monomorph beam with different slenderness ratios L/h and different volume fraction indices is presented by Yang and Xiang [23], in order to gain an insight into the effect of transverse shear deformation. Figure 3 shows the dimensionless w₀/h tip deflection of an FGPM cantilever monomorph actuator with L/h=6 and different values of volume fraction indices (n), that is subjected to a uniform pressure q=10 kN/m² and an applied voltage V=20 V. The tip deflection is increased significantly as the volume fraction index changed from 0→1.75 and then tends to remain constant in the EBT but decreased regularly in the TSDT and FSDT. The relative difference between the TSDT and EBT solutions remained around 2.7–3.0% as the volume fraction index (n) changes from 0→10.

Table 2:

Dimensionless tip deflections of the cantilever FGPM beam under electromechanical load (n=0.2).

	L/h
Theory	4	5	6	10	15	25
Yang and Xiang [23] (EBT)	0.603	1.380	2.758	20.121	100.062	764.572
Present (FSDT)	0.642	1.408	2.799	20.245	100.35	765.64
Present (TSDT)	0.622	1.417	2.81	20.285	100.36	765.78
Relative difference %	3.055	2.611	1.850	0.8084	0.297	0.158

Figure 3:

Tip deflection of the FGPM cantilever monomorph with a varying volume fraction index (L/h=6).

Table 2 presents the dimensionless tip deflection of an FGPM cantilever monomorph beam (n=0.2) with different values of L/h under a uniform pressure q=10 kN/m² and an applied voltage V=20 V. The deflection of beams with three different theories (TSDT, FSDT, and EBT) is listed for direct comparison. The relative difference between the TSDT and EBT solutions, which is defined as [(TSDT − EBT)/TSDT]×100% to show the effect of shear deformation is not negligible for L/h┼6 (thick beam), indicating that the TSDT beam theory is more accurate in this case. For the slender actuators with L/h>6 this difference is less than 1% and decreased as the value of L/h is increased.

Figures 4 and 5 compare the tip deflections of the FGPM beams (L/h=6) with different volume fraction indices when subjected to a uniform pressure q=10 kN/m² or an applied voltage V=100 V, using the FSDT and TSDT. The results are in good agreement with those reported by Yang and Xiang [23], where GDQ method is used for solving the governed equations. When only an external voltage is applied, the electrically induced bending deformation is primarily dependent upon the piezoelectric coefficient e₃₁ and the applied electric load. So, increasing the volume fraction index (n) leads to fully reached PZT-5 beam with higher value of e₃₁ and consequently less deflection.

Figure 4:

The effect of the volume fraction index on the tip deflection of FGPM beam under mechanical load.

Figure 5:

The effect of the volume fraction index on the tip deflection of FGPM beam under electrical load.

The accuracy of TSDT, FSDT, and EBT was investigated using the PZT4-PZT5 FGPM beam with different slenderness ratios under electromechanical load. Since the EBT dose not account for the effect of the shear deformation, this theory predicted less deflection compared with the FSDT, and TSDT. By increasing the proportion of L/h (slender beams), the effect of the shear deflection is reduced and the results of the TSDT, FSDT and EBT converge together. Figure 6 depicts the transverse deflection of the clamped-clamped beam subjected to electrical load V=10 V. As expected, for the thick beams (e.g., L/h=4), the EBT estimated inaccurate deflection. When the magnitude of L/h is approaching to 20, the effect of transverse shear became negligible, and it may be seen that for L/h=20, the difference of the midpoint deflection of the beam estimated by EBT and TSDT is less than 2%.

Figure 6:

The effect of the slenderness ratio L/h on the deflection of clamped-clamped FGPM beam under electrical load (V=100 V).

A simply supported FGPM beam with the slenderness ratio L/h=6 and L=0.025 m subjected to uniform mechanical load q=10 kN/m² on the top surface is considered to compare the results of Fourier series solution with the finite element solution based on the EBT. Figure 7 presents the electric field E _z at the middle point, across the thickness direction of the FGPM beam, caused by mechanical load obtained from the Fourier series and finite element solution. Based on Figure 7, increasing the volume fraction index (n) in the range 0≤n≤2 leads to the growth of the electric field, and after that the electric field starts to decrease for higher values of n. The reason is the difference between two coefficients e₃₁ and ∊₃₃ for the PZT-4 and PZT-5. Comparing the curves in Figure 7, it is clear that, for 0≤n≤2 (rich PZT-4 in material composition), Fourier series solution and finite element solution coincide properly, but, for n≥2 (increasing the PZT-5 portion in material composition), the results based on the Fourier series solution show higher value of the electric field. Although in both solutions the electric field is assumed to vary linearly through the thickness, but the finite element solution based on the EBT leads to inaccurate results compared to the Fourier series solution. It is predicted that the reason is mainly due to using the approximated shape functions in the finite element solution.

Figure 7:

Comparison of the electric field versus volume fraction indices at the middle of simply supported FGPM beam under mechanical load obtained from Fourier series and finite element solution. (L/h=6,q=10 kN/m²).

Behavior of the simply supported FGPM beam under thermal loading is investigated using the EBT and TSDT. In the present analysis, constant surface temperature is imposed at the top and bottom surfaces as the thermal load.

The temperature of the top PZT-4 rich surface is fixed at 400°K and that of bottom PZT-5 -rich surface is kept constant at 350°K. It is assumed that the reference temperature is T₀ = 295°K. Homogeneous (n=0) and inhomogeneous (n=0.5,2,10) simply supported beams with slenderness ratio L/h=25 are considered for the analysis. Figure 8 compares the deflection of beams with various volume fraction indices. As expected, for homogeneous beam n=0, TSDT and EBT predicted deflection almost as the same, while, for the inhomogeneous beam, EBT is not capable to predict shear deformation, and, thus, it estimated less deflection. By increasing the volume fraction index (e.g., n=10) and making beam more homogeneous, the results of two theories converged again. Figure 9 depicts the upper surface electric field E _z of the FGPM simply supported beam (L/h=25) versus the volume fraction indices under thermal gradient (ΔT=50°C) using TSDT. As expected from (6) and (11), the electric field decreased by increasing the exponent (n) and consequently increasing the ∊₃₃, because the electric potential Φ(z,x) is proportional with $1/{∊}_{33}$ .

Figure 8:

Dimensionless deflection of the simply supported beam under thermal loading (ΔT=50°C) (dash line indicates EBT, stretch line indicates TSDT).

Figure 9:

The upper surface electric field of the FGPM simply supported beam versus volume fraction index under thermal gradient (ΔT=50°C).

The axial stress and transverse shear stress in the PZT4-PZT5 FGPM beam with slenderness ratio L/h=15 under a uniform pressure q=10 kN/m² and voltage V = 20 V were evaluated using the TSDT and FSDT. The results presented correspond to x=L/2. The stresses are nondimensional using the relation, $({\bar{\sigma }}_{xx},{\bar{\tau }}_{xz})=({\sigma }_{xx},{\tau }_{xz})h^2/(q{L}^2)$ . Figures 10 through 12 depict the nondimensional shear and axial stress of a clamped-clamped beam.

Figure 10:

Nondimensional shear stresses ${\bar{\tau }}_{xz}$ in clamped-clamped PZT4-PZT5 beam across the thickness using TSDT (L/h=15).

A thorough examination of these results reveals that for the shear stress, independent of the material combination, profiles are close to one another irrespective of the value of n in Figure 10. This nature of closeness of shear stress profiles may be due to small difference in Young's modulus of PZT4-PZT5 (81.3 GPa and 60.6 GPa, resp.) [23]. Figure 11 shows that, for homogenous beam (n=0,n=400), the FSDT predicted constant shear stresses, but TSDT predicted more accurate second-order polynomial function for shear stress through the thickness (Figure 10). Based on Figures 10 and 11, the condition of ${\bar{\sigma }}_{xx}=0$ is not satisfied using the FSDT. That may be accounted as a drawback for the FSDT compare with TSDT. Both TSDT and FSDT accounted the same axial stress through the thickness in Figure 12, which means that shear deformation does not influence the axial stress.

Figure 11:

Nondimensional shear stresses ${\bar{\tau }}_{xz}$ in clamped-clamped PZT4-PZT5 beam across the thickness using FSDT (L/h=15).

Figure 12:

Nondimensional axial stresses ${\bar{\sigma }}_{xx}$ in clamped-clamped PZT4-PZT5 beam across the thickness using TSDT and FSDT (L/h=15).

It is observed from Figure 10 that the shear stress across the beam thickness is not symmetric about the midplane. The relative location of the neutral plane is not the same for the FGPM beams with different volume fraction indices. The location of the neutral plane has been presented in Table 3 by Δδ through the various volume fraction indices. The slenderness ratio L/h=15 is considered with L=1.2 m. It is observed through Table 3 that the value of Δδ is greater when TSDT is employed. This occurred because for the FSDT shear stress remains constant through the thickness, wherein the rotation due to shear assumptions does not hold well in FSDT. For volume fraction indices other than homogeneous composition, the neutral plane shifted towards the stiffer material (PZT-4) face. The maximum shift from midplane occurred for n=2.

Table 3:

Distance of neutral plane from midplane Δδ (mm).

Distance of neutral plane from midplane
n	TSDT	FSDT
0	0	0
0.05	0.1662	0.1650
0.2	0.5747	0.5701
0.5	1.0820	1.0741
1	1.4350	1.4262
1.8	1.5298	1.5190
2	1.5191	1.5080
5	1.1242	1.1371
10	0.7484	0.7413
500	0.0203	0.0200

5. Conclusion

The thermo-electromechanical behavior of functionally graded piezoelectric beams is investigated in this study. A consistent finite element model as well as the Fourier series are presented for the static analysis of FGPM beams. Using Hamilton's principle, the governing differential equations are derived which took into account the effect of transverse shear deformation. The formulation is derived for both TSDT and FSDT. The FE formulation for static analysis is validated by comparing the results of simply-supported, clamped-clamped, and cantilever beams with the previous studies. The TSDT is effectively used to study extensively the static displacement field components, axial stress, and shear stress distribution in various FGPM beams with different slenderness ratios and volume fraction indices. It is found from the present study that the effect of shear deformation should be taken into consideration when L/h┼6 and the volume fraction index change the position of neutral plane. The error of the FSDT solution for the thick beams is large in predicting shear stress across the thickness and yields inaccurate results.

Appendix

We have:

$\begin{array}{cc}\hfill {\acute{A}}_1& ={{\displaystyle \int }}_{-h/2}^{+h/2}{C}_{11}\left(z\right)dz,\hfill \\ \hfill {\acute{A}}_2& ={{\displaystyle \int }}_{-h/2}^{+h/2}-z{C}_{11}\left(z\right)dz,\hfill \\ \hfill {\acute{A}}_3& ={{\displaystyle \int }}_{-h/2}^{+h/2}\frac{e_{31}\left(z\right)}{h}dz,\hfill \\ \hfill {\acute{B}}_1& ={{\displaystyle \int }}_{-h/2}^{+h/2}z{C}_{11}\left(z\right)dz,\hfill \\ \hfill {\acute{B}}_2& ={{\displaystyle \int }}_{-h/2}^{+h/2}{z}^2{C}_{11}\left(z\right)dz,\hfill \\ \hfill {\acute{B}}_3& ={{\displaystyle \int }}_{-h/2}^{+h/2}\frac{z{e}_{31}\left(z\right)}{h}dz,\hfill \\ \hfill {\acute{B}}_4& ={{\displaystyle \int }}_0^1{f}_3\sin \left(m\pi \bar{x}\right)dx,\hfill \\ \hfill {\acute{C}}_1& ={{\displaystyle \int }}_{-h/2}^{+h/2}{e}_{31}\left(z\right)dz,\hfill \\ \hfill {\acute{C}}_2& =-{{\displaystyle \int }}_{-h/2}^{+h/2}z{e}_{31}\left(z\right)dz,\hfill \\ \hfill {\acute{C}}_3& =-{{\displaystyle \int }}_{-h/2}^{+h/2}\frac{{∊}_{33}\left(z\right)}{h}dz,\hfill \\ \hfill {\acute{C}}_4& ={{\displaystyle \int }}_{-h/2}^{+h/2}\frac{{∊}_{33}\left(z\right)}{h}dz,\hfill \\ \hfill {\acute{D}}_1& =-{{\displaystyle \int }}_{-h/2}^{+h/2}{e}_{31}\left(z\right)dz,\hfill \\ \hfill {\acute{D}}_2& ={{\displaystyle \int }}_{-h/2}^{+h/2}z{e}_{31}\left(z\right)dz,\hfill \\ \hfill {\acute{D}}_3& ={{\displaystyle \int }}_{-h/2}^{+h/2}\frac{{∊}_{33}\left(z\right)}{h}dz,\hfill \\ \hfill {\acute{D}}_4& =-{{\displaystyle \int }}_{-h/2}^{+h/2}\frac{{∊}_{33}\left(z\right)}{h}dz.\hfill \end{array}$

(23)