The FEM analysis of FGM piezoelectric semiconductor problems

Composite Structures 163 (2017) 13–20

Contents lists available at ScienceDirect

Composite Structures journal homepage: www.elsevier.com/locate/compstruct

The FEM analysis of FGM piezoelectric semiconductor problems J. Sladek a,⇑, V. Sladek a, H.H.-H. Lu b, D.L. Young b a b

Institute of Construction and Architecture, Slovak Academy of Sciences, 84503 Bratislava, Slovakia Department of Civil Engineering, National Taiwan University, Taipei 10617, Taiwan

a r t i c l e

i n f o

Article history: Received 6 October 2016 Revised 24 November 2016 Accepted 5 December 2016 Available online 9 December 2016 Keywords: Finite element method Functionally graded materials Piezoelectric semiconductors

a b s t r a c t The finite element method is developed for 3-D general boundary value problems for a piezoelectric semiconductor with functionally graded material properties. The electron density and electric current are additionally considered in the constitutive equations for piezoelectric semiconductors. A general variation of material properties with Cartesian coordinates is treated in the numerical analyses. The influence of material parameter gradation and initial electron density is investigated in the static case. Numerical results are presented for a 3-D beam under a static and impact mechanical load. Transversely isotropic material properties are considered in this study. Ó 2016 Elsevier Ltd. All rights reserved.

1. Introduction Piezoelectric materials (PE) are frequently utilized in smart structures and devices for many engineering applications. These materials have attracted considerable attention and are used in the production of transducers, sensors and actuators as they offer certain potential advantages in their performance over those using conventional materials. They have the capability to convert energy from one form to another (among electric and mechanical) [1–4]. Piezoelectric materials can be characterized as either dielectrics or semiconductors. In PE semiconductors (conducting PE) the induced electric field produces also the electric current. The acousto-electric effect in piezoelectric semiconductors is the interaction between the mechanical fields and the moving charges [5,6]. Yang and Zhou [7] showed that acoustic wave traveling can be amplified by an application of an initial or biased DC electric field in piezoelectric semiconductors. Amplification of acoustic waves is employed in acousto-electric receivers and transductors [8,9]. Other more complicated models to describe behavior of piezoelectric semiconductors have also been reported in the literature, see, e.g., [10,11]. Considerably more attention has, hitherto, been paid to investigate problems involving dielectric materials than with semiconductors, as the former are relatively simpler to treat than the latter. New advanced materials and structures often have to satisfy multiple functions. A state-of-the-art review and suggestions for ⇑ Corresponding author. E-mail address: [email protected] (J. Sladek). http://dx.doi.org/10.1016/j.compstruct.2016.12.019 0263-8223/Ó 2016 Elsevier Ltd. All rights reserved.

future directions to increase innovation in developing these materials and structures is given by Ferreira et al. [12]. Functionally graded materials (FGMs) can be considered as a special case of multifunctional materials, where the volume fractions of its constituents vary in a predominant direction [13]. At the macroscopic structural level, one can consider these materials as a continuum with graded material properties. These materials have been introduced to take advantage of the superior properties of its constituents. This notion has been extended to piezoelectric materials to achieve a low dielectric constant [14,15]. The solution of a general boundary value problem for coupled fields in continuously varying material properties in the medium requires advanced numerical methods due to the high mathematical complexity involved. The finite element method (FEM) and the boundary element method (BEM) have been frequently applied to solve boundary value problems. Allik and Hughes [16] first utilized the FEM to investigate the behavior of piezoelectric ultrasonic transducers, and many researchers have also delved into the related topics [17,18]. The BEM is an alternative, powerful computational tool if the fundamental solution is available. The BEM has been successfully applied to piezoelectricity by Lee [19], and Ding and Liang [20]. Wünsche et al. [21,22] have also applied the BEM to fracture problems in 2-D piezoelectric solids. Unfortunately, the fundamental solutions or Green’s functions are not available for 3-D problems in piezoelectric anisotropic materials or in a medium with continuously varying properties. Since governing equations for piezoelectric semiconductors are more complicated than for non-conducting piezoelectrics, the fundamental solution is not available here too. Then, it is possible to derive only boundary-

14

J. Sladek et al. / Composite Structures 163 (2017) 13–20

domain integral formulation with a simple test function in the weak form and advantages of a pure boundary integral formulations are lost. Meshless formulations have also become due to their high adaptivity and relatively low cost to prepare the input data for the numerical analysis. A variety of such methods has been applied to piezoelectric problems [23–25] and functionally graded material properties in non-conducting dielectric piezoelectrics have also been considered in the numerical analyses reported in [26]. To this end, recent reviews of meshless methods include those by Nguyen et al. [27] and by Sladek et al. [28]. There are only a few papers devoted to the solution of boundary value problems in piezoelectric semiconductor materials in the literature. They deal with the anti-plane crack problem in an unbounded domain with a semi-infinite crack [29] or a finite crack [30] under stationary conditions. The Fourier transform technique was applied to reduce the problem to a pair of dual integral equations. Yang et al. [31] have derived 2-D equations for multilayered shells of piezoelectric semiconductors to investigate the propagation of torsional waves. Sladek et al. [32–34] analyzed 2-D inplane and anti-plane crack problems in a finite domains under a mechanical, thermal and electric load. To the best of the authors’ knowledge, there is no paper on the FEM analysis of 3-D general piezoelectric semiconductor solids. In this paper, the governing equations with the corresponding boundary conditions for 3-D problems are derived from the variational principle, and from which the FEM formulation is subsequently developed and implemented. Numerical examples are presented in which a simply supported beam is analyzed in 3-D for the displacements under a mechanical load, electric potential and electron density for the coupled problem. The influence of material parameter gradation as well as the initial electron density on the responses is also investigated. 2. Governing equations for piezoelectric semiconductor

rij ðx; sÞ ¼ cijkl ðxÞekl ðx; sÞ  ekij ðxÞEk ðx; sÞ ð1Þ

J i ðx; sÞ ¼ qM0 lij Ej ðx; sÞ  qdij M;j ðx; sÞ where rij , Dj , and q are stress tensor, electric displacement field, and the electric charge of electron, respectively. Symbols M and Ji are used for the electron density field variable and the density of electric current, respectively. Material parameters cijkl ðxÞ, eijk ðxÞ, hij ðxÞ, lij ðxÞ and dij ðxÞ are the elastic, piezoelectric, dielectric, electron mobility and carrier diffusion material coefficients, respectively. For FGMs, these coefficients are dependent on the spatial coordinates. The strain tensor eij and the electric field vector Ej are related to the displacements ui and the electric potential / by

1 2

eij ¼ ðui;j þ uj;i Þ; Ej ¼ /;j

r11 3 2 c11 6r 7 6c 6 22 7 6 12 6 7 6 6 r33 7 6 c13 6 7 6 6r 7 ¼ 6 0 6 23 7 6 6 7 6 4 r13 5 4 0 r12

c12 c11 c13 0 0

c13 c13 c33 0 0

0 0 0 c44 0

0

0

0

0

0 0 0 0 c55

0 0 0 0 0

3 2 e11 0 76 e 7 6 0 76 22 7 6 76 7 6 76 e33 7 6 0 76 7 6 76 2e 7  6 0 76 23 7 6 76 7 6 54 2e13 5 4 e15 32

2e12

0 c66

0

0 0 0 e15 0 0

2 3 3 e11 e31 6 7 7 2 3 e31 72 3 6 e22 7 E1 6 7 E1 7 6 e33 7 e33 76 7 7 74 E2 5 ¼¼ C6 7  L6 4 E2 5; 6 2e 7 0 7 6 23 7 7 E3 6 7 E3 7 4 2e13 5 0 5 2e12

0

ð3Þ

2

e11

2

3

e11

3

6e 7 7 32 3 2 3 36 6 e22 7 2 6 22 7 0 0 0 0 e15 0 6 E1 D1 E1 h11 0 0 6 7 7 6e 7 e 33 7 6 7 6 6 7 6 7 7 76 7 þ 4 0 h11 0 54 E2 5 ¼¼ G6 33 7 þ H6 4 D2 5 ¼ 4 0 0 0 e15 0 0 56 4 E2 5; 6 2e 7 6 2e 7 6 32 7 6 32 7 e31 e31 e33 0 0 0 6 D3 E3 E3 0 0 h33 6 7 7 4 2e31 5 4 2e31 5 2

3

2

2e12

2e12 ð4Þ

2 3 2 l11 J1 6 7 6 4 J 2 5 ¼ qM0 4 0 0 J3

0

0

32

E1

3

2

d11 0

0

32

M;1

3

2

E1

3

2

M;1

3

6 7 6 76 7 6 7 6 7 l22 0 7 54 E2 5  q4 0 d22 0 54 M;2 5 ¼¼ qM0 A4 E2 5  qF4 M;2 5; 0 l33 M;3 M;3 E3 E3 0 0 d33 ð5Þ T

1 ðc11 2

where c66 ¼  c12 Þ and G ¼ L . The quasi-static approximation of the first Maxwell equation is assumed, since the frequency of external loadings is assumed to be significantly lower than the characteristic frequency of electromagnetic fields changes. Consider a piezoelectric solid of domain V with boundary C. The electro-elastic Gibbs free energy density function U is given as

U¼

1 1 1 rij eij þ Dk /;k þ Jk M;k 2 2 2

ð6Þ

Then, the variation of free energy is equal to R R dU ¼ V ðrij deij þ Dk d/;k þ J k dM ;k ÞdV ¼ V ðrij dui;j þ Dk d/;k þ J k dM;k ÞdV R R ¼  V ðrij;j dui þ Dk;k d/ þ Jk;k dMÞdV þ C ðnj rij dui þ nk Dk d/ þ nk J k dMÞdC R R ¼  V ½rij;j dui þ Dk;k d/ þ J k;k dMdV þ C ft i dui þ Q d/ þ SdMgdC ð7Þ

We consider an one-carrier piezoelectric semiconductor of ntype with continuous, non-homogeneous material properties and with electron density M 0 in the unloaded state with vanishing initial electric field E0 . A general variation of material properties with Cartesian coordinates is considered. The constitutive equations represent the coupling of the mechanical and electrical fields and electric current fields [5,34]

Dj ðx; sÞ ¼ ejkl ðxÞekl ðx; sÞ þ hjk ðxÞEk ðx; sÞ

2

ð2Þ

Many piezoelectric ceramic materials exhibit transversely isotropic elastic behavior with hexagonal symmetry of 6 mm class with the poling direction parallel to the x3 -axis and x1  x2 plane being the isotropic plane. The Voigt notation is used here for material coefficients in the matrix form of the constitutive equations:

The variation of the work done by external loads is given by

Z

dW ¼

dWdV Z Z Z Z _ t i dui dC þ þ ¼  qMd/dV þ qMdMdV V

V

Z þ

CS

V

Ct

CQ

SdMdC

 d/dC Q ð8Þ

_ is the body where –(qM) represents the body electric charge; qM electric current; and prescribed quantities are denoted by the overbar. The variation of the kinetic energy is written as

Z

dK ¼

Z

dKdV ¼ V

V

qu_ i du_ i dV

ð9Þ

€ i is the acceleration. where q is the mass density and u From the Hamilton principle for deformable bodies, Rt ðdW þ dK  dUÞds ¼ 0, the following governing equations are 0 obtained,

rij;j ðx; sÞ ¼ qu€i ðx; sÞ Dk;k ðx; sÞ ¼ qMðx; sÞ _ sÞ þ J ðx; sÞ ¼ 0 qMðx;

ð10Þ

k;k

and the two kinds of boundary conditions (b.c.) can be stated as follows:

i ðx; sÞ on Cu ; Cu  C Essential b:c: : ui ðx; sÞ ¼ u  Mðx; sÞ ¼ Mðx; sÞ on CM ; CM  C  sÞ on C/ ; C/  C /ðx; sÞ ¼ /ðx;

ð11Þ

15

J. Sladek et al. / Composite Structures 163 (2017) 13–20

Natural b:c: : t i ðx; sÞ ¼ ti ðx; sÞ on Ct ; Ct [ Cu ¼ C; Ct \ Cu ¼ £ Sðx; sÞ ¼ Sðx; sÞ on CS ; CS [ CM ¼ C; CS \ CM ¼ £  ðx; sÞ on CQ ; CQ [ C/ ¼ C; CQ \ C/ ¼ £ Q ðx; sÞ ¼ Q

ð12Þ The initial conditions for the primary fields are assumed as

ui ðx; sÞjs¼0 ¼ ui ðx; 0Þ; u_ i ðx; sÞjs¼0 ¼ u_ i ðx; 0Þ; /ðx; sÞjs¼0 ¼ /ðx; 0Þ; Mðx; sÞjs¼0 ¼ Mðx; 0Þ in domainV: 3. The finite element method formulation In the variational formulation of the FEM, the weak-form of a boundary value problem in piezoelectricity can be derived from the Hamilton variational principle discussed in the previous section. It can be written as

Z t (Z 0

V

€ i dui þ qMd/ ½rij dui;j þ Dk d/;k þ J k dM;k þ qu

¼

Ct

0

V

Z

tdui dC þ

ð17Þ

In view of (16) and (17), the variational formulation (14) can be rewritten as

Z h

CQ

 Cþ Qd/d

CS

CQ

SdMdC ds

ð13Þ

T

fJg ¼ ðJ1 J2 J 3 Þ

T

T

r22 r33 r23 r13 r12 Þ ; feg ¼ ð e11 e22 e33 2e23 2e13 2e12 Þ

Then, the integral statement (130 ) can be written as

^ t fugdV ¼ ½fdegT frg  fdEgT fDg þ fdKgT fJg þ fdugT ½D V Z Z Z fdugT fTgdC þ fdugT fQgdC þ fdugT fSgdC ¼

fSg ¼ ð 0

0

2

q @t@ 2

B B 0 B B ^ ½Dt  ¼ B 0 B B @ 0 0

0 2

q @t@ 2 0

t 3

0

0

0

0

q

@2 @t 2

0

0

0

T

0 Þ ; fQ g ¼ ð 0

0

0

 Q

0

0

0

0

0 q @t@

n X fqðeaÞ ðtÞgNa ðn1 ; n2 ; n3 Þ;

ð uðeaÞ 1 ðtÞ

ðeaÞ u2 ðtÞ

ðeaÞ u3 ðtÞ

ðeaÞ

ð19Þ

/

ðtÞ M

ðeaÞ

ðtÞ Þ

T

where the superscripts, (ea), denote the global number of the a-th node of V e , (a ¼ 1; 2; . . . ; n). The gradients in the Cartesian coordinate system can then be expressed in terms of the derivatives with respect to intrinsic coordinates as 1 0 1 0 1 @=@x1  @=@n1 @x1 =@n1 @x2 =@n1 @x3 =@n1  B C B C B C B C C B C e B e e 1 e B @=@x2 C ¼ ½Y B @=@n2 C; ½Y  ¼ ½J  ; ½J  ¼ B @x1 =@n2 @x2 =@n2 @x3 =@n2 C @ A @ A @ A   @=@x3 V e @=@n3 @x1 =@n3 @x2 =@n3 @x3 =@n3 V e 0 ðecÞ c 1 ðecÞ ðecÞ x1 @N =@n1 x2 @N c =@n1 x3 @N c =@n1 C n B X B ðecÞ c C B x @N =@n xðecÞ @N c =@n xðecÞ @N c =@n C ¼ 2 2 2C B 1 2 3 A c¼1 @ ðecÞ

ðecÞ

ðecÞ

x1 @N c =@n3 x2 @N c =@n3 x3 @N c =@n3 0 c 1 @N =@n1 n B C X B c C Þ: ¼ B @N =@n2 Cð x1ðecÞ x2ðecÞ xðecÞ 3 @ A c¼1

T

1

@N c =@n3

Hence,

 n X @f  ðeaÞ ea ¼ f bi ðn1 ; n2 ; n3 Þ;  @xi V e a¼1 @Na @Na @Na þ Y e12 þ Y e13 ; @n1 @n2 @n3 @Na @Na @Na ea þ Y e22 þ Y e23 ; b2 ðn1 ; n2 ; n3 Þ ¼ Y e21 @n1 @n2 @n3 @Na @Na @Na ea b3 ðn1 ; n2 ; n3 Þ ¼ Y e31 þ Y e32 þ Y e33 @n1 @n2 @n3 ea

Furthermore, from the constitutive equations, the energetically conjugated fields are given by

frg ¼ ½Cfeg  ½LfEg; fDg ¼ ½LT feg þ ½HfEg; fJg ¼ qM0 ½AfEg þ q½FfKg

ðtÞg ¼

0Þ ;

C 0 C C C : 0 C C C q A

0

ðeaÞ

ð14Þ

CS

fTg ¼ ð t 1 t 2 T 0 0 0 S Þ ,

fug  fuðx; tÞgjV e ¼

0

T

with

ð18Þ

The primary field variables in 3-D can be approximated over each standard finite element, V e , in terms of their respective nodal values and the corresponding shape function, Na, as

fq

T T fug ¼ ðu1 u2 u3 / M Þ ; fDg ¼ ðD1 D2 D3 Þ ; fKg ¼ ðM;1 M;2 M;3 Þ

CQ

CS

0

ð13 Þ

CS

Ct

CQ

a¼1

It is appropriate to introduce the following column vectors:

frg ¼ ð r11

^e   ½L½B ^ E Þ  fdEgT ð½LT ½B ^ e  þ ½H½B ^ E Þþ fdegT ð½C½B i ^ E  þ q½F½B ^ M Þ þ fdugT ½D ^ t  fugdV þ fdKgT ðqM 0 ½A½B Z Z Z ¼¼ fdugT fTgdC þ fdugT fQgdC þ fdugT fSgdC Ct

)

Z

_ € i dui þ qMd/  qMdMdV ½rij dui;j þ Dk d/;k þ J k dM;k þ qu Z Z Z  SdMdC tdui dC þ Qd/d ¼ Cþ Ct

Z

0 0 0 0 C C ^ b2 0 0 0 C 0 1 C ^1 0 0 0 0 b C ^ C B C 0 b3 0 0 C ^ e fug; fEg ¼ B ^ E fug; ^2 0 C Cfug ¼ ½B B0 0 0 b Cfug ¼ ½B C @ A ^1 0 0 C 0 b C ^3 0 0 0 0 b ^3 b ^2 0 0 C C b A ^1 0 0 0 ^2 b b 1 0 1 0 ^1 M ;1 0 0 0 0 b B C B C B ^i ¼ @ ^ M fug; with b C ^2 C fKg ¼ B Cfug ¼ ½B 0 0 0 0 b @xi @ M ;2 A ¼ B @ A ^3 M ;3 0 0 0 0 b

V

where we have utilized dui js¼0 ¼ 0 ¼ dui js¼t . From (13), it can be seen that the following weak form should be satisfied at any time instant

Z

^1 b B B B0 B B B0 B feg ¼ B B^ B b3 B B B0 @

) _  qMdMdV ds Z t (Z

1

0

ð16Þ

Since the secondary fields feg, fEg and fKg can be expressed in terms of the primary field variables fug as

b1 ðn1 ; n2 ; n3 Þ ¼ Y e11

ð20Þ

and dVjV e ¼ j det½J e jdn1 dn2 dn3 .Thus, from (17) and (20), the approximations of the secondary fields on V e are given by

16

J. Sladek et al. / Composite Structures 163 (2017) 13–20

0

1 ea b1 0 0 0 0 B 0 bea 0 0 0 C B C 2 B C ea  n X B C 0 0 0 0 b  3 ðeaÞ ^ e fug ¼ B C; feg  fegjV e ¼ ½B ½Bea g; ½Bea e fq e  ¼ B ea ea C e V 0 b 0 0 b B C 3 1 a¼1 B C ea ea @ 0 b3 b2 0 0 A ea ea b2 b1 0 0 0 0 1 ea 0 0 0 b1 0  n X  B C ea ðeaÞ ^ E fug ¼ ½Bea g; ½Bea fEg  fEgjV e ¼ ½B E fq E  ¼ @ 0 0 0 b2 0 A; Ve ea a¼1 0 0 0 b3 0 0 ea 1 0 0 0 0 b1  n X B ea C ðeaÞ ^M fug ¼ ½Bea g; ½Bea fKg  fKgjV e ¼ ½B M fq M  ¼ @ 0 0 0 0 b2 A e V

a¼1

0 0 0 0

ea b3

ð21Þ

In the case of static problems, the time derivative terms are vanishing and we have to deal with the system of the algebraic equations instead of the system of the ODE. 4. Numerical examples A piezoelectric beam, as shown in Fig. 1, with homogeneous material properties is analyzed in the first example. The following geometry parameters are considered: length l, width w, and height h are 0.1 m, 0.005 m, and 0.01 m, respectively. The material properties correspond to those of aluminum nitride (AlN) [35]: ð0Þ

ð0Þ

ð0Þ c13

ð0Þ c23

ð0Þ

c11 ¼ c22 ¼ 4:03  1011 Nm2 ; c12 ¼ 1:43  1011 Nm2 ; ¼

11

¼ 1:04  10

ð0Þ

Using the FEM approximations, the variations of the primary field variables are expressed in terms of the variations of the nodal values as

fdug  fdugjV e ¼

n X fdqðeaÞ gNa ðn1 ; n2 ; n3 Þ; fdqðeaÞ g

¼ ð duðeaÞ 1

ðeaÞ

ðeaÞ

du2

du3

d/ðeaÞ

dMðeaÞ Þ

Similarly, variations of the secondary fields are given as

fdeg  fdegjV e

n X

ð23Þ

ðeaÞ ½Bea g E fdq

ð24Þ

a¼1 ðeaÞ g fdKg  fdKgjV e ¼ sumna¼1 ½Bea M fdq

ð25Þ

Substituting equations (19)–(25) into the variational statement (18) yields the following:

Z

m X n X

1

Z

1

e¼1 a;c¼1

1

1

Z

1

fdqðeaÞ g

T

n

1

T T ec ½Bea E  ð½L ½Be 

T

ec ec ½Bea e  ð½C½Be   ½L½BE Þþ T ec ½Bea M  ðqM 0 ½A½BE 

Cet

e¼1 a¼1

CeS

ð0Þ

ð0Þ

ð0Þ

e15 ¼ 0:39 Cm2 ; e31 ¼ 0:66 Cm2 ; e33 ¼ 1:57 Cm2 ; ð0Þ

ð0Þ

h11 ¼ h22 ¼ h33 ¼ 8:092  1011 C ðmVÞ1 ;

q ¼ 3255 kgm3 ; q ¼ 1:602  1019 C; M0 ¼ 1:0  106 m3 where the superscripts, ð0Þ, denote the grading material property at x3 ¼ 0.The beam is discretized by 2000 (40  5  10) brick elements with a uniform distribution. A uniform static compressive load is applied at the top of the beam. The complete set of the prescribed boundary conditions is given as follows:

~ ~t3 ðxÞ ¼ 100Pa; ~t 1 ðx; tÞ ¼ ~t 2 ðxÞ ¼ QðxÞ ¼ ~SðxÞ ¼ 0; ~ ~ ~tðxÞ ¼ 0; /ðxÞ ¼ MðxÞ ¼ 0;

þ ½H½Bec þ q½F½Bec E Þ þ M Þþ o a c ^ e ðecÞ þN N ½Dt  fq gj det½J jdn1 dn2 dn3 Z Z Z m X n X T ¼¼ fdqðeaÞ g ð fTgNa dC þ fSgN a dC þ fQ gNa dCÞ 

ð0Þ

d11 ¼ d22 ¼ d33 ¼ 7:0  104 m2 s1 ;

a¼1

fdEg  fdEgjV e ¼

ð0Þ

l11 ¼ l22 ¼ l33 ¼ 3:0  102 m2 ðVsÞ1 ;

T

n X ðeaÞ ¼ ½Bea g e fdq

Nm ;

c33 ¼ 3:82  1011 Nm2 ; c44 ¼ c55 ¼ 1:20  1011 Nm2 ;

ð0Þ

ð22Þ

a¼1

2

CeQ

ð26Þ

at x 2 Cjx3 ¼0

~ ~tðxÞ ¼ 0; QðxÞ ¼ ~SðxÞ ¼ 0; atx 2 Cjx2 ¼0 [ Cjx2 ¼w ~ ~tðxÞ ¼ 0; QðxÞ ¼ ~SðxÞ ¼ 0; atx 2 Cjx1 ¼0 [ Cjx1 ¼l To preclude rigid body motion, and taking into account the ~ 3 ðxÞ ¼ 0 at symmetry, we assume u

where Cet ¼ Ce \ Ct ; CeS ¼ Ce \ CS ; CeQ ¼ Ce \ CQ . Denoting the global nodes xðeaÞ as xk , k 2 f1; 2; . . . ; Kg, the sumPm Pn ðeaÞ gðÞea in (26) can be replaced by mations e¼1 a¼1 fdq PK P ek a k k , where ek denotes the finite element V ek involvk¼1 ek fdq gðÞ ing the node xk , and ak is the local number of the node xk on V ek . Since (26) is valid for any arbitrary fdqk g satisfying the boundary conditions, the following system of ordinary differential equations is obtained: n Z XX

1

1

ek c¼1

Z

1

1

n T e c ½Beek ak  ð½C½Beek c   ½L½BEk Þ

e a T ½BEk k  ð½LT ½Beek c 

e c

e a T

e c

þ ½H½BEk Þ þ ½BMk k  ðqM0 ½A½BEk  o ek c ^ t  fqðek cÞ gj det½J ek jdn1 dn2 dn3 þ q½F½BM Þ þ þNak Nc ½D Z Z XZ ¼ ð e fTgNak dC þ e fSgNak dC þ e fQ gNak dCÞ; þ

ek

Ct k

ðk ¼ 1; 2; . . . ; KÞ

CSk

CQk

ð27Þ

at x 2 Cjx3 ¼h

Fig. 1. A simply supported piezoelectric beam.

17

J. Sladek et al. / Composite Structures 163 (2017) 13–20

x 2 Cjx1 ¼0 [ Cjx1 ¼l

~ 1 ðxÞ ¼ u ~ 2 ðxÞ ¼ 0 at and x3 ¼ h=2; u

x 2 Cjx1 ¼0;x3 ¼h=2 [ Cjx1 ¼l;x3 ¼h=2 andx2 ¼ 0: The variations of the beam deflection, the electric potential, and the electron density along the x1 -axis at x2 ¼ 0 and x3 ¼ h=2, are presented in Figs. 2–4, respectively. The accuracy of the present computational method is verified by results obtained using an alternative numerical method, namely, the local radial basis function collocation method (LRBFCM) [36]. Good agreement between the two sets of numerical results can be observed. The induced electric potential has the largest value at the center of beam, x1 ¼ l=2. The biggest difference between the FEM and LRBCFM results is observed at the ends of the beam. It can be due to a reduced support domain there. Nodes in a support domain are

0

2,00E+07

-5E-10

1,80E+07

-1E-09

1,60E+07

-1,5E-09

1,40E+07

LRBCFM

-2E-09

FEM

Μ [m-3]

u3 [m]

used for approximation of physical fields in the LRBCFM. In the FEM the approximation has a local character and there is no difference between the approximation at the center or at the end of the beam. This is similarly observed, in Fig. 4, for the variation of the density of electrons. After verifying the computer code for almost non-conducting piezoelectric material, the influence of the initial electron density on the beam deflection, induced electric potential and the electron density is analyzed. From Fig. 5 and 6, one can observe that there is virtually no influence of this quantity on the beam deflection and electric potential. However, the electron density is strongly affected by the initial electron density as can be seen in Fig. 7; it is significantly larger if the initial electron density is increased.

-2,5E-09 -3E-09

1,20E+07 1,00E+07 8,00E+06

-3,5E-09

6,00E+06

LRBCFM

-4E-09

4,00E+06

FEM

-4,5E-09

2,00E+06 0,00E+00

-5E-09 0

0,2

0,4

0,6

0,8

0

1

0,2

0,6

0,8

1

x1 [m]

x1 [m] Fig. 2. Variation of the transverse displacement along the x1 -axis at x2 ¼ 0 and x3 ¼ h=2 for static analysis.

Fig. 4. Variation of electron density along x1 -axis at x2 ¼ 0 and x3 ¼ h=2 for static analysis.

-0,1

0 -5E-10

-0,15

LRBCFM

-1E-09

FEM

-0,2

-1,5E-09 -0,25

u3 [m]

φ [V]

0,4

-0,3

M0 = 0.

-2E-09

M0 = 1.e7 [m-3]

-2,5E-09 -3E-09

-0,35

-3,5E-09

-0,4

-4E-09 -0,45

-4,5E-09

-0,5

-5E-09 0

0,2

0,4

0,6

0,8

1

x1 [m] Fig. 3. Variation of the electric potential along x1 -axis at x2 ¼ 0 and x3 ¼ h=2 for static analysis.

0

0,2

0,4

0,6

0,8

1

x1 [m] Fig. 5. Variation of the vertical displacements along the line ( x1 ; x2 ¼ 0; x3 ¼ h=2) for different values of M 0 in the beam under static load, t3 ¼ 100 Pa.

18

J. Sladek et al. / Composite Structures 163 (2017) 13–20

The influence of functionally graded material parameters on the beam deflection, electric potential and electron density is analyzed in the last example. The same material gradation is considered for all elastic, piezoelectric, and dielectric coefficients. A common symbol P is used to denote all material parameters. The exponential variation of material parameters is considered along the x3 direction

Pðx3 Þ ¼ Pð0Þ edx3 =h ; Where Pð0Þ , d, and h denote the material property at x3 ¼ 0, grading parameter, and beam height, respectively. The value Pð0Þ corresponds to the material used in the previous example for a homogeneous material. Basically, a general gradation of material parameters with different variations for individual parameters can be considered.

Variations of the beam deflection, electric potential and electron density along the x1 -axis for two different gradation exponents d and the initial electron density M0 ¼ 1:0  106 m3 are presented in Figs. 8–10. The elastic coefficient increases in value along the x1 -direction for positive gradation exponent d. The beam deflection is reduced if the elastic coefficients are larger. If the beam deflection (strain) is smaller, the induced electric potential and electron density are reduced as well. In the next example we analyze the same homogeneous beam under a pure mechanical load with Heaviside time variation t3 Hðs  0Þ. The same boundary conditions as those in the static case and the zero initial conditions including the initial electron

0 homogeneous

-5E-10

FGM: delta = ln(3/2)

-0,1

-1E-09

-0,15

M0 = 1.e7 [m-3]

u3 [m]

φ [V]

-1,5E-09

M0 = 0.

-0,2

delta = ln(2)

-0,25

-2E-09 -2,5E-09 -3E-09

-0,3

-3,5E-09

-0,35

-4E-09

-0,4

-4,5E-09 -5E-09

-0,45

0

0,2

0,4

-0,5 0

0,2

0,4

0,6

0,8

0,6

0,8

1

x1 [m]

1

x1 [m] Fig. 6. Variation of electric potential along the line (x1 ; x2 ¼ 0; x3 ¼ h=2) for different values of M 0 in the beam under static load, t 3 ¼ 100 Pa.

Fig. 8. Variation of the vertical displacements along the line ( x1 ; x2 ¼ 0; x3 ¼ h=2) for different values of the gradation parameter d in the FGM beam under static load, t3 ¼ 100 Pa.

0

2,50E+08

homogeneous

-0,05

FGM: delta = ln(3/2)

-0,1

2,00E+08

delta = ln(2)

1,50E+08

φ [V]

Μ [m-3]

-0,15 -0,2 -0,25 -0,3

1,00E+08 M0 = 1.e6 [m-3]

-0,35

M0 = 1.e7 [m-3]

-0,4

5,00E+07

-0,45 -0,5

0,00E+00 0

0,2

0,4

0,6

0,8

1

x1 [m] Fig. 7. Variation of density of electrons along the line ( x1 ; x2 ¼ 0; x3 ¼ h=2) for different values of M 0 in the beam under static load, t 3 ¼ 100 Pa.

0

0,2

0,4

0,6

0,8

1

x1 [m] Fig. 9. Variation of electric potential along the line ( x1 ; x2 ¼ 0; x3 ¼ h=2) for different values of the gradation parameter d in the FGM beam under static load, t3 ¼ 100 Pa.

19

J. Sladek et al. / Composite Structures 163 (2017) 13–20

2.5

2,00E+07 1,80E+07

2

1,60E+07

φ/φstat

M [m-3]

1,40E+07 1,20E+07 1,00E+07

1.5

1

8,00E+06 6,00E+06

non-conducting

0.5 homogeneous

4,00E+06

M0=1.e7[1/m**3]

FGM: delta = ln(3/2)

2,00E+06

0

delta = ln(2)

0

0.00005

0,00E+00 0

0,2

0,4

0,6

0,8

0.00015

0.0002

0.00025

time [sec]

1

x1 [m]

0.0001

Fig. 12. Temporal variation of the induced electric potential at the center ðx1 ¼ l=2; x2 ¼ 0; x3 ¼ h=2Þ under an impact load with Heaviside variation.

Fig. 10. Variation of the electron density along the line ( x1 ; x2 ¼ 0; x3 ¼ h=2) for different values of the gradation parameter d in the FGM beam under static load, t3 ¼ 100 Pa.

5. Conclusions

2.5 non-conducting

u3/u3stat

2

M0=1.e7[1/m**3]

1.5

1

0.5

0 0

0.00005

0.0001

0.00015

0.0002

0.00025

The finite element method is developed to analyze threedimensional boundary value problems under stationary boundary conditions and with functionally graded material properties in piezoelectric semiconductors. It follows from the results for a simply supported beam that the value of the initial electron density has no influence on the transverse displacement and the induced electric potential. However, the same cannot be said for the induced electron density. The influence of functionally graded material properties in the piezoelectric beam has also been investigated. If piezoelectric material coefficient increases in value along the beam height, the beam deflection, induced electric potential and electron density are reduced when compared to a beam with corresponding homogeneous properties. Therefore, the gradation of material properties can be utilized for optimal design of piezoelectric semiconductor structures. Acknowledgements

time [sec] Fig. 11. Temporal variation of the beam deflection at the center ðx1 ¼ l=2; x2 ¼ 0; x3 ¼ h=2Þ under an impact load with Heaviside variation.

density are considered. The time step Ds ¼ 1  105 ½sec is applied in numerical analyses. Influence of electric conductivity on the beam deflection and induced electric potential is investigated. Temporal variations of the beam deflection and induced electric potential at the point ðx1 ¼ l=2; x2 ¼ 0; x3 ¼ h=2Þ for various initial electron densities M 0 are presented in Figs. 11 and 12, respectively. Both quantities are normalized by the values corresponding to the static loading. One can observe that initial electron density has only a small influence on the beam deflection and induced electric potential if beam is under a pure mechanical load. A similar conclusion has been made for a static case. As can be seen in 11–12, the transient responses at the point

ðx1 ¼ l=2; x2 ¼ 0; x3 ¼ h=2Þ oscillate around the static values. Then, the peak values are about two times larger than static values.

The authors acknowledge the financial supports from the Slovak Science and Technology Assistance Agency registered under number APVV-14-0216. We are also grateful for the financial supports by the National Science Council of Taiwan under the grant nos. 103-2923-E-002-004-MY2 and 104-2811-E-002-012. References [1] Avellaneda M, Harshe G. Magnetoelectric effect in piezoelectric/magnetostrictive multilayer (2–2) composites. J Intell Mat Syst Struct 1994;5:501–13. [2] Berlingcourt DA et al. Piezoelectric and piezomagnetic materials and their function in transducers. Phys Acoust 1964;1:169–270. [3] Landau LD, Lifshitz EM. In: Lifshitz EM, Pitaevskii LP, editors. Electrodynamics of Continuous Media. New York: Pergamon Press; 1984. [4] Nan CW. Magnetoelectric effect in composites of piezoelectric and piezomagnetic phases. Phys Rev B 1994;50:6082–8. [5] Hutson AR, White DL. Elastic wave propagation in piezoelectric semiconductors. J Appl Phys 1962;33:40–7. [6] White DL. Amplification of ultrasonic waves in piezoelectric semiconductors. J Appl Phys 1962;33:2547–54. [7] Yang JS, Zhou HG. Amplification of acoustic waves in piezoelectric semiconductor plates. Int J Solids Struct 2005;42:3171–83. [8] Heyman JS. Phase insensitive acoustoelectric transducer. J Acoust Soc Am 1978;64:243–9.

20

J. Sladek et al. / Composite Structures 163 (2017) 13–20

[9] Busse LJ, Miller JG. Response characteristics of a finite aperture, phase insensitive ultrasonic receiver based upon the acoustoelectric effect. J Acoust Soc Am 1981;70:1370–6. [10] De Lorenzi HG, Tiersten HF. On the interaction of the electromagnetic field with heat conducting deformable semiconductors. J Math Phys 1975;16:938–57. [11] Maugin GA, Daher N. Phenomenological theory of elastic semiconductors. Int J Eng Sci 1986;24:703–31. [12] Ferreira ADBL et al. Multifunctional material systems: a state-of-the-art review. Compos Struct 2016;151:3–35. [13] Suresh S, Mortensen A. Fundamentals of functionally graded materials: The Institute of Materials, London; 1998. [14] Zhu X et al. A functionally gradient piezoelectric actuator prepared by metallurgical process in PMN-PZ-PT system. J Mat Sci Lett 1995;14:516–8. [15] Han F et al. Responses of piezoelectric, transversally isotropic, functionally graded and multilayered half spaces to uniform circular surface loading. CMES: Comp Model Eng Sci 2006;14:15–30. [16] Allik H, Hughes TJR. Finite element method for piezoelectric vibration. Int J Num Meth Eng 1970;2:151–7. [17] Benjeddou A. Advances in piezoelectric finite element modeling of adaptive structural elements: a survey. Comput Struct 2000;76:347–63. [18] Saravanos DA et al. Layerwise mechanics and finite element for the dynamic analysis of piezoelectric composite plates. Int J Solids Struct 1997;34:359–78. [19] Lee JS. Boundary element method for electroelastic interaction in piezoceramics. Engn Anal Bound Elem 1995;15:321–8. [20] Ding H, Liang J. The fundamental solutions for transversely isotropic piezoelectricity and boundary element method. Comput Struct 1999;71:447–55. [21] Wünsche M et al. A 2D time-domain collocation-Galerkin BEM for dynamic crack analysis in piezoelectric solids. Engn Anal Bound Elem 2010;34:377–87. [22] Wünsche M et al. Dynamic crack analysis in piezoelectric solids with nonlinear electrical and mechanical boundary conditions by a time-domain BEM. Comp Meth Appl Mech Eng 2011;200:2848–58.

[23] Ohs RR, Aluru NR. Meshless analysis of piezoelectric devices. Comput Mech 2001;27:23–36. [24] Sladek J et al. The MLPG for bending of electroelastic plates. CMES: Comp Model Eng Sci 2010;64:267–97. [25] Sladek J et al. Laminated elastic plates with piezoelectric sensors and actuators. CMES: Comp Model Eng Sci 2012;85:543–72. [26] Sladek J et al. Evaluation of fracture parameters in continuously nonhomogeneous piezoelectric solids. Int J Fract 2007;145:313–26. [27] Nguyen VP et al. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [28] Sladek J et al. Applications of the MLPG method in engineering & sciences: a review. CMES: Comp Model Eng Sci 2013;92:423–75. [29] Yang J. An anti-plane crack in a piezoelectric semiconductor. Int J Fract 2005;136:L27–32. [30] Hu Y et al. A mode III crack in a piezoelectric semiconductor of crystals with 6mm symmetry. Int J Solids Struct 2007;44:3928–38. [31] Yang J et al. Amplification of acoustic waves in piezoelectric semiconductor shells. J Intell Mat Syst Struct 2005;16:613–21. [32] Sladek J et al. Fracture analysis in piezoelectric semiconductors under a thermal load. Engn Fract Mech 2014;126:27–39. [33] Sladek J et al. Dynamic anti-plane crack analysis in functionally graded piezoelectric semiconductor crystal. CMES: Comp Model Engn Sci 2014;99:273–96. [34] Sladek J et al. Influence of electric conductivity on intensity factors for cracks in functionally graded piezoelectric semiconductors. Int J Solids Struct 2015;59:79–89. [35] Auld BA. Acoustic fields and waves in solids. New York: John Wiley and Sons; 1973. p. 357–82. [36] Lu HHH et al. Three-dimensional analysis for functionally graded piezoelectric semiconductors. J Intell Mater Syst Struct 2016. http://dx.doi.org/10.1177/ 1045389X16672566.