Two-dimensional transient analysis of wave propagation in functionally graded piezoelectric slabs using the transform method

International Journal of Solids and Structures 52 (2015) 72–82

Contents lists available at ScienceDirect

International Journal of Solids and Structures journal homepage: www.elsevier.com/locate/ijsolstr

Two-dimensional transient analysis of wave propagation in functionally graded piezoelectric slabs using the transform method Yi-Hsien Lin a, Yi-Shyong Ing b, Chien-Ching Ma a,⇑ a b

Department of Mechanical Engineering, National Taiwan University, Taipei 10617, Taiwan, ROC Department of Aerospace Engineering, Tamkang University, Tamsui 25137, Taiwan, ROC

a r t i c l e

i n f o

Article history: Received 3 March 2014 Received in revised form 31 August 2014 Available online 6 October 2014 Keywords: FGPM Stress waves Durbin method Numerical Laplace inversion Transient response

a b s t r a c t In this study, the transient responses of a functionally graded piezoelectric slab were analyzed using the transform technique. The slab was subjected to a dynamic anti-plane concentrated force and an in-plane point electric displacement on the top surface, and the bottom surface was assumed to be open-circuit or grounded. The analytical solutions were obtained in the Laplace transform domain, and a numerical inversion was performed using the Durbin method. The numerical results showed that an applied point electric displacement may instantaneously induce shear stress waves due to the piezoelectricity for both open-circuit and grounded cases. When the bottom surface was grounded, a visible electric wave was induced by the grounded boundary and the influence of this wave was significant when the gradient coefficient of the functionally graded piezoelectric materials (FGPMs) was large enough. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Piezoelectric materials possess the important property of coupling between mechanical and electrical fields, which renders them useful in many areas of modern technology. These materials have therefore been widely used for a long time as electromechanical transducers, filters, sensors and actuators, to mention only a few applications. In recent years, various applications of piezoelectric materials have been implemented in non-destructive evaluation, ultrasonic medical imaging, smart structures, and the active control of sound and vibration. To increase the lifetime and reliability of advanced piezoelectric structures, functionally graded material has been considered as a way to improve interface problems. This new type of piezoelectric material is referred to as a functionally graded piezoelectric material. A functionally graded material (FGM) can be prepared by continuously changing the constituents of multi-phase materials in a predetermined volume as a fraction of the constituent material (Khor et al., 1997; Kwon and Crimp, 1997; Nogata, 1997). Most researchers have analyzed the composition of FGPMs with three types of functions, power-law, polynomial, and exponential, which are widely used because these functions provide advantages in theoretical investigation. For the power-law case, Wu et al. (2003) considered the material constants of a FGPM

⇑ Corresponding author. Tel.: +886 2 23659996; fax: +886 2 23631755. E-mail address: [email protected] (C.-C. Ma). http://dx.doi.org/10.1016/j.ijsolstr.2014.09.021 0020-7683/Ó 2014 Elsevier Ltd. All rights reserved.

cylindrical shell as a power-law functions in the radial direction. When an axisymmetric thermal or mechanical loading is applied on the cylindrical shell, an exact solution is obtained through the power series expansion method together with the Fourier series expansion method. For the polynomial FGPM, Yu et al. (2007, 2009) used the Legendre orthogonal polynomial series expansion to determine the wave characteristics in spherically curved plates or hollow cylinders composed of FGPMs with an open-circuit condition. For the FGPM with an exponential function variation, Li et al. (2004) investigated the dispersion relations of Love waves in a layered functionally graded piezoelectric structure for electric open and short cases. Zhong and Shang (2003) presented an exact threedimensional solution of the FGPM rectangular plate for the simply supported and grounded boundary. Zhong and Yu (2008) proposed a two-dimensional general solution for FGPM beams with arbitrary graded functions, and the numerical calculation was based on the cantilever beam with exponential variation. Time-harmonic response of a vertically graded transversely isotropic, linearly elastic half-space is analytically determined by Eskandari-Ghadi and Amiri-Hezaveh (2014) by introducing a new set of potential functions. The potential functions are set in such a way that the governing equations be simple and with physical meaning as well. The study of wave propagation in piezoelectric materials or FGPMs is a rather involved problem. The situation is even more formidable when non-homogeneity has to be considered. It is therefore not surprising that only scant information regarding transient wave propagation problems has been presented. As in

73

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

the case of the Stoneley wave, whose mechanical displacements are in the sagittal plane, the amplitude of this wave decreases with the distance away from the interface into both media (Stoneley, 1924). Bleustein (1968) and Gulayaev (1969) simultaneously discovered that there exists a shear horizontal (SH) electro-acoustic surface mode in a class of transversely isotropic piezoelectric media, which is known today as the BG wave. The BG wave is a unique result in the repertoire of surface acoustic wave (SAW) theory because it has no counterpart in purely elastic solids. As a matter of fact, since its discovery, the BG wave theory has become one of the cornerstones of modern electro-acoustic technology. It  was shown that a BG wave can exist in cubic crystals of 43m and  1 0Þ plane and their 23 classes, along [1 1 0] direction on the ð1 equivalent orientations. The velocity equations for piezoelectric and elastic surface waves were derived and their characteristics were discussed by Tseng (1970). A pure shear elastic surface wave (MT wave) can propagate along the interface of two identical crystals, in class 6 mm, when the z-axes of these crystals, both parallel to the interface and perpendicular to the propagation direction, are in opposite directions (Maerfeld and Tournois, 1971). The general equations and the fundamental piezoelectric matrix derived for the anti-plane wave motion and Floquet theory were applied to obtain the passing and stopping bands in a periodically layered infinite space by Honein and Herrman (1992). Taking into account both the optical effect as well as the contribution from the rotational part of an electric field, the obtained solutions not only were valid for any wave speed range but also provided accurate formulas to evaluate the acousto-optic interaction due to piezoelectricity. As the wave speed is much less than the speed of light, the solution degenerates to the well-known BG wave or MT wave (Li, 1996). The surface acoustic wave (SAW) can be excited and detected efficiently by using an interdigital transducer (IDT) placed on a piezoelectric substrate. Therefore, a vast amount of effort was invested in the research and development of SAW devices for military and communication applications, such as delay lines and filters for radar. The propagation mode in most devices is the Rayleigh wave on a free surface of a piezoelectric substrate. Recently, Ma et al. (2007) used Cagniard’s method to construct the full-field transient solutions of piezoelectric bi-materials. Each term represents a physical transient wave. The existence condition of the surface wave of piezoelectric bi-materials is restricted to the situation where the shear wave speeds of the two piezoelectric materials are very close. Laplace transform and the inversion of Laplace transform are usually used for analyzing the dynamic wave-propagation problem. However, using the analytical inversion of Laplace transform for analysis renders the mathematical composition excessively complex and difficult; therefore, it is only suitable for relatively simple geometric structures. The numerical inversion of Laplace transform is more practical for calculating complex problems. Lin and Ma (2011) used the numerical Laplace inversion and the Durbin method to obtain the transient stress responses for a 20-layered elastic structure. Ma et al. (2012) used the Durbin method to analyze the long-time responses for a functionally graded slab. Lin and Ma (2012) found that the transient responses of the continuously distributed multilayered media can be simulated by the effective material concept and are applicable to analyze the wave propagation problem in a functionally graded slab. Ing et al. (2013) presented an extended Durbin method for a two-sided Laplace inversion, and evaluated the transient stress and electric displacement of a two-layered piezoelectric medium. In this study, the transient wave propagation problem of an FGPM slab was investigated. The slab was subjected to dynamic point loading on the top surface. The transient response was analyzed by employing the Laplace transform method. The interface

and boundary conditions were used to construct the system of equations for determining the field vectors in the FGPM slab. Subsequently, the Durbin method was used to implement the double inversions for the transient response of wave propagation in the FGPM slab. 2. Statement of the problem For a linear FGPM medium, the constitutive equations can be expressed as

rij ¼ cijkl Skl  ekij Ek ; Di ¼ eikl Skl þ eik Ek ;

ð1Þ ð2Þ

where rij, Skl, Ek and Di indicate the stress tensor, strain tensor, electric field vector, and electric displacement vector, respectively. In the absence of body forces, the governing equations can be expressed as

rij;j ¼ qu€i ;

ð3Þ

Di;i ¼ 0;

ð4Þ

where q is the mass density, ui is the displacement vector, and the superposed dot indicates material differentiation with respect to time. The material property parameters of the FGPM are assumed to vary exponentially along the y-direction according to the following exponential law

c44 ¼ c440 eby ;

e15 ¼ e150 eby ;

e11 ¼ e110 eby ; q ¼ q0 eby ;

ð5Þ

where c44, e11, and e15 denote the elastic modulus, dielectric permittivity, and piezoelectric constant, respectively. b represents the gradient coefficient of the FGPM. To reduce the complexity of mathematical analysis, all of the physical properties were assumed to vary in the same way. Consider the following anti-plane displacement and electric potential fields:

u1 ¼ u2 ¼ 0;

u3 ¼ wðx; y; tÞ;

/ ¼ /ðx; y; tÞ:

ð6Þ

From the constitutive relations of FGPM poled in the z-direction, the nontrivial components of stresses and electric displacements are

@w @/ þ e150 eby ; @x @x @w @/ syz ¼ c440 eby þ e150 eby ; @y @y by @w by @/  e110 e ; Dx ¼ e150 e @x @x @w @/  e110 eby : Dy ¼ e150 eby @y @y

sxz ¼ c440 eby

ð7Þ ð8Þ ð9Þ ð10Þ

Substituting Eqs. (5)–(10) into Eqs. (3) and (4), the governing equations can be obtained as follows:

@w @/ € þ be150 ¼ q0 w; @y @y @w @/  be110 ¼ 0; e150 r2 w  e110 r2 / þ be150 @y @y

c440 r2 w þ e150 r2 / þ bc440

ð11Þ ð12Þ

where

r2 ¼

@2 @2 þ 2 2 @x @y

ð13Þ

is the two-dimensional Laplacian. It is worth noting that Eqs. (13) and (12) are two coupled partial differential equations. By introducing a function

w¼/

e150

e110

w;

ð14Þ

then the constitutive equations are reduced to the following forms:

74

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

@w @w þ e150 eby ; @x @x

ð15Þ

@w @w ¼ ~c440 eby þ e150 eby ; @y @y

ð16Þ

sxz ¼ ~c440 eby syz

nj ¼

@w ; Dx ¼ e110 e @x

ð17Þ

@w ; @y

ð18Þ

by

Dy ¼ e110 eby where

e2150

~c440 ¼ c440 þ

e110

ð19Þ

;

and Eqs. (11) and (12) can be split into two uncoupled sets of equations

@w ¼ 0; @y @w 2 € r2 w þ b ¼ b w; @y

r2 w þ b

ð20Þ ð21Þ

where

rffiffiffiffiffiffiffiffi b¼

q0

ð22Þ

~c440

is the slowness of the shear wave. Note that the shear wave speed is cs = 1/b. The one-sided Laplace transform with respect to time and the two-sided Laplace transform with respect to x are defined by

f ðx; y; sÞ ¼

Z

1

f ðx; y; tÞest dt;

ð23Þ

0

^f ðk; y; sÞ ¼

Z

1

mj ¼

b þ ð1Þjþ1

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 b2 þ 4ðb  k2 Þs2

2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi jþ1 b þ ð1Þ b2 þ 4ðe2  k2 Þs2

ð24Þ

j ¼ 1; 2;

;

3 ^ y; sÞ 2 wðk; M 11 7 6 ^  y; sÞ 7 6 M 6 /ðk; 7 6 21 6 7 6 6 ^ 6 sxz ðk; y; sÞ 7 6 M 31 7¼6 6 ^ 7 6 6s 6 yz ðk; y; sÞ 7 6 M 41 7 6 6^ 6 Dx ðk; y; sÞ 7 4 M 51 5 4 M 61 ^ ðk; y; sÞ D y 2

M 14

3

M 12

M 13

M 22

M 23

M 32

M 33

M 42

M 43

M 52

M 53

2 3 M 24 7 7 A1 76 7 M 34 76 A2 7 76 7; M 44 7 74 B1 5 7 M 54 5 B2

M 62

M 63

M 64

M11 ¼ em1 y ; M21 ¼

e150

e110

M 12 ¼ em2 y ;

em1 y ;

M 22 ¼

M 13 ¼ 0;

e150

e110

em2 y ;

M14 ¼ 0;

M23 ¼ en1 y ;

M31 ¼ ~c440 skeðbþm1 Þy ;

M 32 ¼ ~c440 skeðbþm2 Þy ;

M33 ¼ e150 skeðbþn1 Þy ;

M 34 ¼ e150 skeðbþn2 Þy ;

Application of the one-sided and two-sided Laplace transforms to Eqs. (14)–(18) yields

M51 ¼ 0;

M 53 ¼ e110 skeðbþn1 Þy ;

^ y; sÞ ¼ wðk; ^ y; sÞ þ e150 wðk; ^ y; sÞ; /ðk;

M61 ¼ 0;

e110

^ y; sÞ; ^ y; sÞ þ e150 e skwðk; s^xz ðk; y; sÞ ¼ ~c440 e skwðk; ^ y; sÞ ^ y; sÞ @ wðk; @ wðk; s^yz ðk; y; sÞ ¼ ~c440 eby þ e150 eby ; @y @y ^ ðk; y; sÞ ¼ e eby skwðk; ^ y; sÞ; D

ð27Þ

^ y; sÞ ^ ðk; y; sÞ ¼ e eby @ wðk; D : y 110 @y

ð29Þ

x

by

110

ð26Þ

ð28Þ

Similarly, the governing Eqs. (20) and (21) can be transformed to ordinary differential equations as follows:

^ y; sÞ ^ y; sÞ @ 2 wðk; @ wðk; 2 ^ y; sÞ ¼ 0; þb  ðb  k2 Þs2 wðk; 2 @y @y ^ y; sÞ ^ y; sÞ @ 2 wðk; @ wðk; ^ y; sÞ ¼ 0; þ b  ðe2  k2 Þs2 wðk; @y2 @y

1

where

2

M 52 ¼ 0;

M54 ¼ e110 skeðbþn2 Þy ; M 63 ¼ e110 n1 eðbþn1 Þy ;

M 62 ¼ 0;

M64 ¼ e110 n2 eðbþn2 Þy : We define the global field vector c, the response vector b, and the phase-related receiver matrix Rcv as follows:

3 ^ w 6 ^ 7 6 / 7 7 6 7 6s 6 ^xz 7 7 b¼6 ^yz 7; 6s 7 6 6^ 7 6 Dx 7 5 4 ^ D 2

2

3 A1 6A 7 6 27 c ¼ 6 7; 4 B1 5 B2

y

2

3

M 11

M 12

M 13

M 14

6M 6 21 6 6 M 31 R cv ¼ 6 6M 6 41 6 4 M 51 M 61

M 22

M 23

M 32

M 33

M 42

M 43

M 52 M 62

M 53 M 63

M 24 7 7 7 M 33 7 7: M 44 7 7 7 M 54 5 M 64 ð37Þ

ð30Þ ð31Þ

where e is an auxiliary positive real perturbation parameter. It is noted that the final expressions involved are always evaluated at e = 0 in the end of the manipulation. This technique has been widely used in applying transform methods to solve partial differential equations. The general solutions to Eqs. (30) and (31) in the double transformed domain have the following form:

^ y; sÞ ¼ A1 em1 y þ A2 em2 y ; wðk; ^ y; sÞ ¼ B en1 y þ B en2 y ; wðk;

M24 ¼ en2 y ;

M42 ¼ ~c440 m2 eðbþm2 Þy ;

;

M44 ¼ e150 n2 eðbþn2 Þy ;

by

ð36Þ

where Mij are phase-related receiver elements and can be expressed as

M43 ¼ e150 n1 eðbþn1 Þy ;

ð25Þ

ð35Þ

> 0. The four undetermined functions A1, A2, B1, and B2 can be resolved explicitly using the related boundary conditions. Following the process of a global matrix method proposed by Lee and Ma (1999), the following relations in matrix form can be obtained from Eqs. (25)–(29), (32) and (33):

M41 ¼ ~c440 m1 e

1

ð34Þ

2 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2 2 in which Re b2 þ 4ðe2  k2 Þs2 b þ 4ðb  k2 Þs2 > 0 and Re

ðbþm1 Þy

f ðx; y; sÞeskx dx:

j ¼ 1; 2;

;

ð32Þ ð33Þ

Therefore, Eq. (36) can be rewritten into the following matrix form:

b ¼ Rcv c:

ð38Þ

Once the global field vector c is determined by boundary conditions, the solutions in the double transformed domain (i.e., the response vector b) can be obtained after relating c with the phase-related receiver matrix Rcv. Next, we will derive the solutions for a FGPM slab with different types of boundary conditions, which will be used to determine the unknown coefficients in Eqs. (32) and (33). 3. Determination of unknown coefficients In the following Sections 3.1 and 3.2, two different types of electrical boundary conditions on the bottom surface (i.e., open-circuit and grounded conditions) will be considered, respectively.

75

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

U1 ¼



M 43 ðk; 0; sÞ M 44 ðk; 0; sÞ



M 63 ðk; 0; sÞ M 64 ðk; 0; sÞ

;

L2 ¼





M 41 ðk; h; sÞ M42 ðk; h; sÞ ; M 61 ðk; h; sÞ M62 ðk; h; sÞ ð51Þ

and

2

M 41 ðk; 0; sÞ M42 ðk; 0; sÞ

0

0

3

7 6 7 6 M 61 ðk; 0; sÞ M62 ðk; 0; sÞ 0 0 7 6 D¼6 7; 6 0 0 M 43 ðk; h; sÞ M 44 ðk; h; sÞ 7 5 4 M 63 ðk; h; sÞ M 64 ðk; h; sÞ 3 0 0 M 43 ðk; 0; sÞ M44 ðk; 0; sÞ 7 6 6 0 0 M 63 ðk; 0; sÞ M64 ðk; 0; sÞ 7 7 6 U¼6 7; 7 60 0 0 0 5 4 0

0

2

Fig. 1. Geometrical configuration for FGPM slab.

3.1. The open-circuit boundary conditions Consider an initially undisturbed and infinitely long (in the x-direction) FGPM plate as shown in Fig. 1. The behavior of the plate is confined to the x–y plane. At time t = 0, an anti-plane concentrated force of Pfm(t)d(x) and an in-plane point electric displacement of D0fe(t)d(x) are suddenly applied to the top surface of y = 0. The boundary conditions can be written as

syz ðx; 0; tÞ ¼ Pfm ðtÞdðxÞ;

ð39Þ

Dy ðx; 0; tÞ ¼ D0 f e ðtÞdðxÞ;

ð40Þ

where P and D0 indicate the magnitudes of the anti-plane force and electric displacement, respectively. The bottom surface of y = h is regarded as traction-free and open-circuit, thus the boundary conditions are

0 0 2

0 0

0 0

0 0

3

7 6 6 0 0 0 07 7 6 L¼6 7: 6 M 41 ðk; h; sÞ M42 ðk; h; sÞ 0 0 7 5 4 M 61 ðk; h; sÞ M62 ðk; h; sÞ 0 0 ð52Þ From Eq. (49), the global field vector c is given by

c ¼ M1^t:

ð53Þ

syz ðx; h; tÞ ¼ 0;

ð41Þ

The inverse of the coefficient matrix M is represented by extracting the diagonal block matrix D out of the expression as

Dy ðx; h; tÞ ¼ 0:

ð42Þ

M ¼ DðI  RÞ;

Applying the double Laplace transforms to Eqs. (39)–(42) gives

where

s^yz ðk; 0; sÞ ¼ P^f m ðsÞ; ^ ðk; 0; sÞ ¼ D ^f ðsÞ; D

ð43Þ

R ¼ D1 ðL þ UÞ:

0 e

ð44Þ

Therefore, the global field vector c is represented by

s^yz ðk; h; sÞ ¼ 0;

ð45Þ

c ¼ ðI  RÞ1 s;

^ ðk; h; sÞ ¼ 0; D y

ð46Þ

where the source vector s is

y

^ ^ where f m and f e are double transformed functions of fm and fe, t is defined as respectively. The external load vector ^

2

b ¼ Rcv ðI  RÞ1 s; ð47Þ

y

By using Eq. (36), Eqs. (43)–(46) can be represented as

3 2 32 3 ~c440 m1 ~c440 m2 A1 e150 n1 e150 n2 P^f ðsÞ 6 m 7 7 6 6 7 6 ^ 0 0 e110 n1 e110 n2 7 6 76 A2 7 6 D0 f e ðsÞ 7; 7 4 ~c440 m1 em1 h ~c440 m2 em2 h e150 n1 en1 h e150 n2 en2 h 54 B1 5 ¼ 6 4 0 5 B2 0 0 e110 n1 en1 h e110 n2 en2 h 0 2

ð48Þ

or in compact form as

Mc ¼ ^t;

ð49Þ

where the coefficient matrix M is a 4  4 matrix and can be expressed as follows:

M¼DþLþU¼



D1

U1

L2

D2



ð50Þ

;

where   M41 ðk; 0; sÞ M42 ðk; 0; sÞ ; D1 ¼ M61 ðk; 0; sÞ M62 ðk; 0; sÞ

 D2 ¼

ð55Þ

ð56Þ

ð57Þ

By substituting Eq. (56) into Eq. (38), we can obtain

3

3 2 ^  7 6 P f m ðsÞ 7 6 ^ 6 D ðk; 0; sÞ 7 6 ^ y 7 6 D0 f e ðsÞ 7 ^t ¼ 6 7: 7¼ 6 7 ^yz ðk; h; sÞ 7 6 6s 5 4 0 5 4 ^ ðk; h; sÞ 0 D

s^yz ðk; 0; sÞ

s ¼ D1^t:

ð54Þ

 M43 ðk; h; sÞ M44 ðk; h; sÞ ; M63 ðk; h; sÞ M64 ðk; h; sÞ

ð58Þ

where b is the response vector of the FGPM slab in a double transformed domain. With the transformed solution at hand, the inverse Laplace transforms were performed to obtain the transient solution in time domain. It is noted that the matrix formulation proposed in this study is effective and can be extended to solve complex structures such as multilayered FGPM problems. For the problem of a FGPM slab, the four unknown coefficients can be derived explicitly. The process is tricky, and therefore the detailed derivation is omitted. Finally, we obtained A1, A2, B1, and B2 as follows:

  P^f m ðsÞ 1 ;  ð1 þ aÞ m1 ~c440 1  eðm1 m2 Þh  ðm m Þh  P^f m ðsÞ e 1 2 A2 ¼  ;  ð1 þ aÞ m2 ~c440 1  eðm1 m2 Þh   D0 ^f e ðsÞ 1 ; B1 ¼   n1 e110 1  eðn1 n2 Þh  ðn n Þh  ^ D0 f e ðsÞ e 1 2 B2 ¼ ;  n2 e110 1  eðn1 n2 Þh

A1 ¼

ð59aÞ

ð59bÞ

ð59cÞ

ð59dÞ

76

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

ð60Þ

n B1 ¼ P^f m ðsÞ  ð1 þ aÞeðm1 þm2 n1 n2 Þh f31 þ D0 ^f e ðsÞ    f32 eðm1 n1 Þh  f33 eðm2 n1 Þh  D1 ;

ð70cÞ

Substituting Eqs. (59a)–(59d) into Eq. (36), the analytical solutions for shear stresses and electric displacements in the double transformed domain were found to be

n ^ ^ B2 ¼ Pf m ðsÞ  ð1 þ aÞeðm1 þm2 n1 n2 Þh f41 þ D0 f e ðsÞ    f42 eðm1 n2 Þh  f43 eðm2 n2 Þh  D1 ;

ð70dÞ

where

a¼

D0 e150 : Pe110

    1þa eðbþm1 Þy  s^xz ðk; y; sÞ ¼ P^f m ðsÞsk  m1 1  eðm1 m2 Þh    ðm m Þhþðbþm Þy      2 1þa e 1 2 a eðbþn1 Þy     m2 n1 1  eðm1 m2 Þh 1  eðn1 n2 Þh    ðn n Þhþðbþn Þy  1 2 2 a e  ; ð61Þ þ n2 1  eðn1 n2 Þh   ðbþm Þy  1 e  eðm1 m2 Þhþðbþm2 Þy s^yz ðk; y; sÞ ¼ P^f m ðsÞ  ð1 þ aÞ 1  eðm1 m2 Þh  ðbþn Þy  ðn1 n2 Þhþðbþn2 Þy 1 e e ; a 1  eðn1 n2 Þh    eðbþn1 Þy ^ ðk; y; sÞ ¼ D ^f ðsÞsk 1 D x 0 e ðn n Þh n1 1e 1 2   ðn n Þhþðbþn Þy  2 1 e 1 2 ; þ n2 eðn1 n2 Þh  1  ðbþn Þy  ðn n Þhþðbþn2 Þy ^ ðk; y; sÞ ¼ D ^f ðsÞ e 1  e 1 2 : D y 0 e 1  eðn1 n2 Þh

where

f11 ¼ f13 ¼ f21 ¼

ð62Þ

f23 ¼ f31 ¼

ð63Þ

f33 ¼

f41 ¼

n1 n2 e150 m2 n1  ; e 440 e110 e150 C m2 ðn1  n2 Þ

e110

e110

f43 ¼

e110

3.2. The grounded boundary conditions

;

1

m1 n2 e150

e110

e110

n1 ðm2  m1 Þ

e110 1

e110

f22 ¼

m1 n2 n1 n2 e150  ; e 440 e110 e150 C

;

n2 ðm1  m2 Þ

ð64Þ

m2 n2 n1 n2 e150  ; e 440 e110 e150 C

;

n1 n2 e150 m1 n1  ; e 440 e110 e150 C m1 ðn1  n2 Þ

f12 ¼

;

f32 ¼

1

m2 n2 e150

e110

e110

! e 440 m1 m2 C  ; e150

!



e 440 m1 m2 C ; e150

f42 ¼

1

e110

! e 440 n1 m2 e150 m1 m2 C  ; e150 e110 !

e 440 m1 n1 e150 m1 m2 C  ; e150 e110

D ¼ g1 eðm2 n1 Þh þ g2 eðm2 n2 Þh  g3 eðm1 n1 Þh  eðm1 n2 Þh ; Consider an anti-plane concentrated force Pfm(t)d(x) and an inplane point electric-displacement D0fe(t)d(x) applied on the top surface of the FGPM plate, and the bottom surface is regarded as traction-free and grounded. The boundary conditions can be expressed as follows:

syz ðx; 0; tÞ ¼ Pfm ðtÞdðxÞ;

ð65Þ

Dy ðx; 0; tÞ ¼ D0 f e ðtÞdðxÞ;

ð66Þ

syz ðx; h; tÞ ¼ 0;

ð67Þ

/ðx; h; tÞ ¼ 0:

ð68Þ

Applying the double Laplace transforms to the boundary conditions in Eqs. (65)–(68), the boundary value problem is represented as a matrix form:

2

~c440 m1 ~c440 m2 6 0 0 6 4 ~c440 m1 em1 h ~c440 m2 em2 h e150 em1 h e150 em2 h 3 2 P^f ðsÞ 7 6 ^m 7 6  ¼ 6 D0 f e ðsÞ 7: 4 0 5 0

e150 n1 e110 n1 e150 n1 en1 h e110 en1 h

32 3 A1 e150 n2 6 A2 7 e110 n2 7 7 6 7 e150 n2 en2 h 54 B1 5 n2 h B2 e110 e

g1 ¼ g3 ¼

m1 n1 n2 e150

e110 m2 n1 n2 e150

e110



e 440 m1 m2 n1 C ; e150

g2 ¼

e 440 m1 n1 n2 e150 m1 m2 n2 C  ; e150 e110



e 440 m1 m2 n1 C ; e150

g4 ¼

e 440 m2 n1 n2 e150 m1 m2 n2 C  e150 e110

It is noted that the boundary conditions are similar to that of the open-circuit case as presented in Eqs. (35)–(42). The only difference is the electric boundary condition in the bottom surface (i.e., Eqs. (42) and (68)). However, the solutions are quite different and are much more complicated than the open-circuit case. From Eq. (36), the analytical solutions in the double transformed domain for stresses and electric-displacements are expressed as follows:

s^xy ðk; y; sÞ ¼ A1 ~c440 skeðbþm1 Þy þ A2 ~c440 skeðbþm2 Þy þ B1 e150 skeðbþn1 Þy þ B2 e150 skeðbþn2 Þy ;

ð71Þ

s^yz ðk; y; sÞ ¼ A1 ~c440 m1 eðbþm1 Þy þ A2 ~c440 m2 eðbþm2 Þy þ B1 e150 n1 eðbþn1 Þy þ B2 e150 n2 eðbþn2 Þy ;

ð72Þ

^ ðk; y; sÞ ¼ B e skeðbþn1 Þy  B e skeðbþn2 Þy ; D x 1 110 2 110

ð73Þ

The four unknowns in Eq. (69) can be solved explicitly and were obtained:

^ ðk; y; sÞ ¼ B e n eðbþn1 Þy  B e n eðbþn2 Þy : D y 1 110 1 2 110 2

ð74Þ

^ ^ A1 ¼ Pf m ðsÞ  ð1 þ aÞeðm2 n1 Þh f11 þ Pf m ðsÞ  ð1 þ aÞeðm2 n2 Þh f12 i ð70aÞ þ D0^f e ðsÞf13  D1 ;

It is worth noting that the analytical solutions presented in Eqs. (61)–(64) and (71)–(74) are in the double transformed domain. In this study, the Durbin method for the one-sided Laplace transform inversion and the extended Durbin method for the two-sided Laplace transform inversion were used to implement the numerical inversions. Furthermore, the present results can degenerate correctly into the case of slabs made of homogeneous piezoelectric materials.

ð69Þ

h

h A2 ¼ P^f m ðsÞ  ð1 þ aÞeðm1 n1 Þh f21  P^f m ðsÞ i ð1 þ aÞeðm1 n2 Þh f22  D0 ^f e ðsÞf23  D1 ;

ð70bÞ

77

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

4. Numerical Laplace inversion Transient behavior can be obtained via the Laplace transform and numerical inversion proposed by different researchers in the literature, including Papoulis (1957),Miller and Guy (1966), Durbin (1974), Narayanan and Beskos (1982). In this paper, we used the method proposed by Durbin, which is an accurate and efficient method for numerically inverting Laplace-transformed functions, in a combination of finite Fourier sine and cosine series. In the Durbin method, the inverse Laplace transform of a function f ðsÞ is expressed as the following series:

f ðtÞ ¼

Z cþi1 at 1   1 f ðsÞest ds ¼ 2e  Re f ðaÞ 2pi ci1 2 T      N X 2kp 2kpt cos Re f a þ i þ T T k¼0      2k p 2k p t  Im f a þ i sin : T T

ð75Þ

Note that the infinite series involved can only be summed up to N terms, and the transform parameter s is composed of real part a and imaginary part 2kTp:

2p s ¼ a þ ik T

for k ¼ 0; 1; 2; 3; . . . ; N;

ð76Þ

in which T is the total time interval of interest, and the number of equidistant points, N, is a finite positive integer for computing f(t). ^ The formula of the double Laplace inversions of a function f ðk; sÞ

Fig. 2. The variation of material parameters for different gradient coefficient b in FGPM.

2p for k1 ¼ 0; 1; 2; . . . ; N1 ; Tx k3 ¼ 0; 1; 2; . . . ; N3 :

k ¼ ax þ ikj

k2 ¼ 0; 1; 2; . . . ; N2 ; ð79Þ

It was found that atTt = 5 to 10 and axTx = 0.01 to 0.1 can be used to obtain accurate numerical results (Ing et al., 2013). By substituting Eqs. (61)–(64) and (71)–(74) into the formula of the double Laplace inversion, Eq. (77), the transient responses in time domain can be obtained numerically.

can be expressed as follows (Ing et al., 2013):

 at t

 k¼ax  2eax x 1 2e 1  Re ^f ðs; kÞ  Re s¼at 2 2 Tx Tt # ( "   N1 k¼ax X 2k1 pt cos Re ^f ðs; kÞ þ 2k p Tt s¼at þi T1 t k1 ¼0 " #  ))# k¼ax 2k1 pt sin  Im ^f ðs; kÞ 2k p Tt s¼at þi T1 t ( " #

 N k¼ax þi2kT3 p 3 X 2eat t 1 x ^   Re f ðs; kÞ Re þ 2 Tt s¼at k3 ¼0 # ( "   2k p N2 k¼ax þi T3 X 2k2 pt x cos Re ^f ðs; kÞ þ 2k p Tt s¼at þi T2 t k2 ¼0 " #  ))#   k¼ax þi2kT3 p 2k2 pt 2k3 px x sin cos  Im ^f ðs; kÞ 2k2 p s¼at þi T Tt Tx t " #  at t ( 2k3 p 2e 1 k¼ax þi T x  Im  Re ^f ðs; kÞ 2 Tt s¼at # ( "   2k p N k¼ax þi T3 2 X 2k2 pt x cos þ Re ^f ðs; kÞ 2k p s¼at þi T2 Tt t k2 ¼0 " # ))#    )) k¼ax þi2kT3 p 2k2 pt 2k3 px x ^  Im f ðs; kÞ sin sin : 2k p s¼at þi T2 Tt Tx t

f ðt; xÞ ¼

ð77Þ The finite series in Eq. (77) can also be summed up to N1, N2, and N3 terms, and the time transform parameter s is composed of real part at and imaginary part 2kTjtp:

2p s ¼ at þ ikj for k1 ¼ 0; 1; 2; . . . ; N1 ; Tt k2 ¼ 0; 1; 2; . . . ; N2 ; k3 ¼ 0; 1; 2; . . . ; N3 ;

ð78Þ

and the space transform parameter k is composed of real part ax and 2k p imaginary part T jx :

5. Numerical results and discussion In this section, we performed the numerical calculation of the transient response for the FGPM slab. As shown in Fig. 2, the material constant was assumed to be a function of y in the region h 6 y 6 0, where PZT4 was the starting value. The material parameters underwent an exponential change along the negative y direction. The material constants for PZT4 are listed below.

c¼440 2:56  1010 N=m2 ;

e¼150 12:7 C=m2 ;

e¼110 6:46  109 F=m2 ;

b ¼ 3:852  104 s=m: When the top surface of the FGPM was subjected to a concentrated force and a point electric displacement load with a heaviside step function, i.e., fe(t) = fm(t) = H(t), detailed discussions for the bottom surface considered as the open-circuit and grounded cases are presented in Sections 5.1 and 5.2, respectively. The relevant parameter settings for the numerical calculation were atTt = 5, axTx = 0.01, N1 = 1000, N2 = 3000, and N3 = 3000. The static solution for the open-circuit boundary case (Wang, 2003) was compared to verify the accuracy of the transient solutions in this study. Employing the solution proposed by Wang (2003) in this problem results in the following static solution:

sSyz ¼

 P sin ph y 

p 

 : 2h cosh  h x  cos ph y

ð80Þ

5.1. Numerical results for the open-circuit boundary condition First, we considered the case where the gradient coefficient b = 0 in Eqs. (62) and (64) for the numerical calculation. When an anti-plane concentrated force was applied at the origin and the point electric-displacement was neglected, the non-dimensional transient shear stresses for the receivers below the load (0, 0.5h) and not below the load (h, 0.5h) were determined, as shown in Fig. 3. It can be seen that the incident shear wave reached

78

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

Fig. 3. Transient shear stress syz response and static solution for PZT4 (b = 0) for different receiver when the point electric-displacement is neglected (only concentrated force applied) and the bottom surface is open-circuit.

Fig. 5. Shear stress responses in the PZT4 (b = 0) slab for different electricdisplacement loading when the receiver is set at (0, 0.5h) and the bottom surface is open-circuit.

the receiver below the load (0, 0.5h) when the non-dimensional time t/bh = 0.5 (i.e., the arrival time of the incident shear wave to (0, 0.5h) is t = 0.5h/cs = 0.5bh) and the reflected shear wave from the bottom surface arrived at t/bh = 1.5 (i.e., the arrival time of the first reflected shear wave from the bottom surface to (0, 0.5h) is t = 1.5h/cs = 1.5bh). Subsequently, positive and negative values of reflected shear waves appeared alternately at t/bh = 2.5, 3.5, and 4.5. When the receiver was set away from the loading (i.e., (h, 0.5h)), the arrival times would delay either for the source wave (t/bh  1.12) or reflections (t/bh  1.80, 2.69, 3.64, and 4.61). Fig. 3 also shows that the transient solution oscillated near the static solution and rapidly tended to the static value after the arrival of each shear wave. When a point electric-displacement was applied at the origin and the mechanical concentrated force was absent, the transient responses of PZT4 (b = 0) for different receivers, (0, 0.5h) and (h, 0.5h), are shown in Fig. 4. Because of the piezoelectricity, the applied point electric-displacement at the origin instantaneously induced both the acoustic SH and electric waves in the PZT4 slab. Fig. 4 shows that the incident shear wave arrived at the receiver (0, 0.5h) at non-dimensional time t/bh = 0.5, and the first four

reflected shear waves successively reached the receiver at nondimensional time t/bh = 1.5, 2.5, 3.5, and 4.5, respectively. We can see from Fig. 4 that based on the hypothesis that the speed of the electric waves is infinite, at the instant t = 0, the significant stress values are induced. Similarly, when the receiver was set at the location (h, 0.5h), the incident shear stress wave induced by the electric displacement load arrived at non-dimensional time t/ bh = 1.12. Then, the first four reflections arrived at t/bh  1.80, 2.69, 3.64, and 4.61, respectively. In addition, because only an electric displacement load was applied, the static solution would approach zero when time reached to infinity. Fig. 5 shows the transient shear stresses at (0, 0.5h) when an anti-plane concentrated force and various electric displacement loads were applied to the origin (0, 0). This figure shows that the electric displacement effects induced significant stress values at the instant t = 0. It also demonstrated that the sign of the starting values was opposite the sign of values a. It is interesting to see from Fig. 5 that the transient response for a = 1 jumped to the static value at the instant when the loads were applied. This was because one part of the effects from the electric displacement load cancelled the response from the mechanical load. It can be verified that if we set a = 1 in Eq. (62), only the last term remained. Fig. 6

Fig. 4. Transient shear stress syz response and static solution for PZT4 (b = 0) for different receiver when the concentrated force is neglected (only point electricdisplacement applied) and the bottom surface is open-circuit.

Fig. 6. Shear stress syz responses at (0, 0.5h) for different gradient coefficient b when only unit concentrated force applied at the origin and the bottom surface is open-circuit.

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

Fig. 7. Electric-displacement responses of PZT4 (b = 0) for different receiver when the bottom surface is open-circuit.

Fig. 8. Electric-displacement responses at (0, 0.9h) for different gradient coefficient b when the bottom surface is open-circuit.

Fig. 9. Shear stress syz responses at (0, 0.5h) for different electric displacement loading a for PZT4 (b = 0) when the bottom surface is grounded.

79

Fig. 10. Shear stress responses of PZT4 (b = 0) for different receiver when unit electric-displacement loading ^f e ¼ 1=s applied at the origin and the bottom surface is grounded.

Fig. 11. Shear stress responses at (h, 0.5h) for different electric-displacement loading a for PZT4 (b = 0) when the bottom surface is grounded.

Fig. 12. Shear stress responses at (h, 0.5h) for different electric-displacement loading a for FGPM (b = 2) when the bottom surface is grounded.

80

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

Fig. 13. Shear stress syz responses at for different gradient coefficient b when unit concentrated force applied at the origin and the bottom surface is grounded.

shows the transient shear stresses at (0, 0.5h) for the various gradient parameters b. When only an anti-plane concentrated force was applied to the origin (0, 0), the transient shear stresses for parameter b = 2, 1, 0, 1, and 2 were represented by red, blue, purple, green, and black lines, respectively. It can be seen from Fig. 6 that the greater transient shear stress was obtained for a smaller parameter b (or the larger variation of material constants). Fig. 7 shows the transient electric displacement responses of PZT4 (b = 0) for different receivers. It can be seen from Fig. 7 that the electric displacement responses jumped to the corresponding static value instantaneously. The electric displacement responses were not influenced by the arrivals of shear wave fronts, and they became smaller the further the receiver was away from the loading. Regardless of whether an anti-plane concentrated force was applied on the top surface of PZT4, the transient electric displacement responses have the same results. In other words, a concentrated force does not affect the transient behavior of the electric displacement under the open-circuit boundary condition. It can be verified that only the magnitude of the electric displacement D0 was present in Eq. (64). When an in-plane electric displacement load was applied on the top surface of the FGPM slab, the transient electric displacement at (0, 0.9h) for a different gradient coefficient b is shown in Fig. 8. It can be seen from Fig. 8 that the greater electric displacement value was obtained for a smaller parameter b. This was because all of the material properties (including the dielectric permittivities and piezoelectric constants) increase more rapidly from the top surface to the bottom surface for the case with a smaller b. The same phenomenon can also be seen in Fig. 2. 5.2. Numerical results for the grounded boundary condition

Fig. 14. Electric-displacement responses for different receiver and different electric-displacement loading when a unit concentrated force applied at the origin and the bottom surface is grounded.

Fig. 15. Electric-displacement responses at (0, 0.9h) for different gradient coefficient b when the bottom surface is grounded.

In this subsection, the lower surface of the FGPM slab was considered as a grounded boundary for numerical calculation. The transient responses for shear stresses and electric displacements are shown in Figs. 9–15. Here, the gradient parameter b = 0 was considered for numerical calculation first. When an anti-plane concentrated force was applied to the origin, the non-dimensional transient shear stresses syz at (0, 0.5h) for different electric displacement loads are shown in Fig. 9. The results show that the first four shear waves reached the receiver at non-dimensional time t/bh = 0.5, 1.5, 2.5, 3.5, and 4.5, respectively. Similarly, the electric displacement effects induced significant stress values at the instant t = 0. In addition, the larger the value a, the smaller the produced stress at t = 0. The transient behavior in Fig. 9 had similar features as that indicated in Fig. 5, except for the following two phenomena: (1) The transient responses in Fig. 9 had slight changes at time t/bh = 1 and t/bh = 3. This was because when the incident shear stress wave arrived at the bottom surface, it would instantaneously generate the electric wave, which influenced the shear stress value. (2) Under the different electric displacement loads, the stress values for the grounded boundary were larger than that for the open-circuit boundary. If the gradient coefficient b = 0 and the anti-plane concentrated force was absent, the transient shear stresses syz for different receivers, (0, 0.5h) and (h, 0.5h), are shown in Fig. 10. Fig. 10 shows that the first five shear waves reached the observation point (0, 0.5h) at non-dimensional time t/bh = 0.5, 1.5, 2.5, 3.5, and 4.5, respectively. When the receiver was moved to the right side with a distance h (i.e., (h, 0.5h)), the arrival time of the incident wave qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 equals h þ ð0:5hÞ =cs  1:12bh. Similarly, the arrival time of the first reflected wave from the bottom surface equals qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 h þ ð1:5hÞ =cs  1:80bh, the arrival time of the second reflected qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 wave from the top surface equals h þ ð2:5hÞ =cs  2:69bh, and so on. However, Fig. 10 shows that some significant waves reached

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

the receiver at t/bh = 0.5, 1.5, 2.5, 3.5, and 4.5. This was because the ground boundary induced a shear wave at the instant t = 0 under the electric displacement load. Then the shear wave propagated back and forth between the upper and bottom surfaces. When an anti-plane concentrated force was applied at the origin, the transient shear stresses under various electric displacement loads at the receiver (h, 0.5h) for the PZT4 (b = 0) and FGPM (b = 2) slabs are shown in Figs. 11 and 12, respectively. We found from Figs. 11 and 12 that adding the functionally graded effect did not influence the arrival times of shear wave fronts. This was because when all of the functionally graded parameters only underwent similar changes along the y direction, the speed of the shear wave remained unchanged. However, this does affect the value of the transient shear stress and possibly lowers its value. When an anti-plane concentrated force was applied to the origin and the electric displacement load was neglected, the transient shear stresses at the receiver (0, 0.5h) for various gradient parameters (b = 2, 1, 0, 1, and 2) are shown in Fig. 13. The transient behavior in Fig. 13 had similar features as that indicated in Fig. 6, except for the effects induced by the electric wave at time t/ bh = 1 and t/bh = 3. When b = 0 and an anti-plane concentrated force was applied to the origin, the electric displacement responses for various receivers are shown in Fig. 14. Because of the existence of electric potential difference (grounded boundary), the transient electric displacements had significant changes at non-dimensional time t/bh = 1.0 and t/bh = 3.0. When the electric displacement load was absent (a = 0), Fig. 14 shows the phases of transient electric displacement were irrespective of the observation points. This implies that when the shear stress waves arrived at the bottom surface of the PZT4 slab, the effects of the electric displacement would be induced at the same time. We also found that a negative electric displacement change (negative slope) was generated at t/bh = 1.0 and t/bh = 3.0 when the receiver lay vertically below the point load. However, a positive change was generated when the receiver was not vertically below the point load. The black line shown in Fig. 14 indicates the transient electric displacement for the observation point (0, 0.5h) under the electric displacement load a = 1. The results show that because of the electric effects, the numerical value of the transient electric displacement was highly visible from the beginning. The values of the transient electric displacement clearly varied when the shear stress waves (t/bh = 1.0 and t/bh = 3.0) reached the receiver. Moreover, some significant changes also occurred at non-dimensional time t/bh = 2.0 and t/bh = 4.0. This was because when the electric wave coming from the electric displacement load arrived instantaneously at the bottom surface, it would generate the shear stress wave that propagated back and forth in the PZT4 slab. These repeatedly reflected shear stress waves then affected the transient electric displacement and caused slight changes at t/bh = 2.0, 4.0, and so forth. When an electric displacement load was applied to the origin and the mechanical concentrated force was absent, the transient electric displacements at (0, 0.9h) for various gradient parameters (b = 2, 1, 0, 1, and 2) are shown in Fig. 15. When the gradient parameter b was decreased, Fig. 15 shows that the sudden jumps caused by the arrivals of shear waves would become more obvious, and the corresponding transient electric displacement value would be larger. By comparing Fig. 8 with Fig. 15, the transient electric displacements have significant differences between the open-circuit and grounded boundary conditions when a point electric-displacement is applied on the top surface of the FGPM. The influences of shear waves on the electric displacement response can be seen distinctly for the grounded boundary. In other words, an external electric displacement load would induce the shear waves and affect the electric displacement responses due to the existence of an electric potential difference.

81

6. Conclusions This study analyzed the transient responses of an FGPM slab. The slab was subjected to a dynamic anti-plane concentrated force and an in-plane electric displacement load at the top surface, and the bottom surface was assumed to be open-circuit or grounded. The transient responses for shear stresses and electric displacements were obtained explicitly in the Laplace transform domain. The Durbin method was then used to implement a numerical inversion. The Durbin method is suitable for calculating transient responses of FGPM slabs and we suggested preferable parameters (N1 = 1000, N2 = 3000, N3 = 3000, at Tt = 5, and axTx = 0.01) for employing it. From the numerical results, the following conclusions can be drawn: (1) The concentrated force has no influence on the transient behavior of the electric displacement under the open-circuit boundary condition. (2) Because one part of the effects from the electric displacement load cancelled the response from the mechanical load, the transient response of shear stresses for a = 1 under the open-circuit boundary condition jumped to the static value at the instant t = 0. (3) The ground boundary could induce a significant shear wave at the instant t = 0 under the electric displacement load. The induced shear stress wave had a considerable influence on the transient shear stress and electric displacement. (4) Either for the open-circuit or grounded boundary conditions, the smaller transient values of the shear stresses and electric displacements can be obtained for the larger gradient parameter b. Moreover, the sudden jumps caused by the arrivals of shear waves would become more obvious when the gradient parameter of the FGPMs was smaller.

Acknowledgments The authors gratefully acknowledge the Grant NSC 101-2923-E002-017-MY3 from the National Science Council, Republic of China, awarded to National Taiwan University. References Bleustein, J.L., 1968. A new surface wave in piezoelectric materials. Appl. Phys. Lett. 13, 412–413. Durbin, F., 1974. Numerical inversion of Laplace transforms: an efficient improvement to Dubner and Abate’s method. Computer J. 17, 371–376. Eskandari-Ghadi, M., Amiri-Hezaveh, A., 2014. Wave propagations in exponentially graded transversely isotropic half-space with potential function method. Mech. Mater. 68, 275–292. Gulyaev, Y.V., 1969. Electroacoustic surface waves in solids. Sov. Phys. JETP Lett. 9, 37–38. Honein, B., Herrman, G., 1992. Wave propagation in non-homogeneous piezoelectric materials. Act. Contr. Noise Vib. 38, 105–112. Ing, Y.S., Liao, H.F., Huang, K.S., 2013. The extended Durbin method and its application for piezoelectric wave propagation problems. Int. J. Solids Struct. 50, 4000–4009. Khor, K.A., Gu, Y.W., Dong, Z.L., 1997. Plasma spraying of functionally graded NiCoCrAlY/Yttria stabilized ZrO2 coating using composite powders. In: Srivatsan, T.S. et al. (Eds.), Composites and Functionally Graded Materials, vol. 80. American Society of Mechanical Engineerings, Materials Division (Publication) MD, pp. 89–105. Kwon, P., Crimp, M., 1997. Automating the design process and powder processing of functionally gradient materials. In: Srivatsan, T.S. et al. (Eds.), Composites and Functionally Graded Materials, vol. 80. American Society of Mechanical Engineerings, Materials Division (Publication) MD, pp. 73–88. Lee, G.S., Ma, C.C., 1999. Transient elastic waves propagating in a multi-layered medium subjected to in-plane dynamic loadings I. theory. Proc. R. Soc. London A. 456, 1355–1374. Li, S., 1996. The electromagneto-acoustic surface wave in a piezoelectric medium: the Bleustein–Gulyaev mode. J. Appl. Phys. 80, 5264–5269. Li, X.Y., Wang, Z.K., Huang, S.H., 2004. Love waves in functionally graded piezoelectric materials. Int. J. Solids Struct. 41, 7309–7328.

82

Y.-H. Lin et al. / International Journal of Solids and Structures 52 (2015) 72–82

Lin, Y.H., Ma, C.C., 2011. Transient analysis of elastic wave propagation in multilayered structures. CMC-Comput. Mat. Contin. 24 (1), 15–42. Lin, Y.H., Ma, C.C., 2012. The investigation of effective material concept for the transient wave propagation in multilayered media. J. Mech. 28 (2), 247–260. Ma, C.C., Lin, Y.H., Lin, S.H., 2012. Transient wave propagation in a functionally graded slab and multilayered medium subjected to dynamic loadings. CMCComput. Mat. Contin. 31 (1), 37–64. Ma, C.C., Chen, X.H., Ing, Y.S., 2007. Theoretical transient analysis and wave propagation of piezoelectric bi-materials. Int. J. Solids Struct. 44, 7110–7142. Maerfeld, C., Tournois, P., 1971. Pure shear elastic surface wave guided by the interface of two semi-infinite media. Appl. Phys. Lett. 19, 117–125. Miller, M.K., Guy, W.T., 1966. Numerical inversion of the Laplace transform by use of Jacobi polynomials. SIAM J. Numer. Anal. 3, 624–635. Narayanan, G.V., Beskos, D.E., 1982. Numerical operational methods for timedependent linear problems. Int. J. Numer. Methods. Eng. 18, 1829–1854. Nogata, F., 1997. Learning about design concepts from natural functionally graded materials. In: Srivatsan, T.S. et al. (Eds.), Composites and Functionally Graded Materials, vol. 80. American Society of Mechanical Engineerings, Materials Division (Publication) MD, pp. 11–18. Papoulis, A., 1957. A new method of inversion of the Laplace transform. Quart. Appl. Math. 14, 405–414.

Stoneley, R., 1924. Elastic waves at the surface of separation of two solids. Proc. R. Soc. London. A 106, 416–428. Tseng, C.C., 1970. Piezoelectric surface waves in cubic crystals. J. Appl. Phys. 41, 2270–2276. Wang, S.Y., 2003. The theoretical analysis of mechanical and electric field of piezoelectric materials for multi-layered, wedge shape and circular boundary problems. Master’s thesis, National Taiwan University, ROC. Wu, X.H., Shen, Y.P., Chen, C.Q., 2003. An exact solution for functionally graded piezothermoelastic cylindrical shell as sensors or actuators. Mater. Lett. 57, 3532–3542. Yu, J.G., Wu, B., Huo, H.L., He, C.F., 2007. Wave propagation in functionally graded piezoelectric spherically curved plates. Phys. Status Solidi (b) 244 (9), 3377– 3389. Yu, J.G., Wu, B., Chen, G.Q., 2009. Wave characteristics in functionally graded piezoelectric hollow cylinders. Arch. Appl. Mech. 79, 807–824. Zhong, Z., Shang, E.T., 2003. Three-dimensional exact analysis of a simply supported functionally gradient piezoelectric plate. Int. J. Solids Struct. 40, 5335–5352. Zhong, Z., Yu, T., 2008. Electroelastic analysis of functionally graded piezoelectric material beams. J. Intell. Mater. Syst. Struct. 19, 707–713.