# A Coupling Electromechanical Inhomogeneous Cell-Based Smoothed Finite Element Method for Dynamic Analysis of Functionally Graded Piezoelectric Beams

Authors: Bin Cai; Liming Zhou
Journal: Advances in Materials Science and Engineering
 (2019)

---

## Abstract

To accurately simulate the continuous property change of functionally graded piezoelectric materials (FGPMs) and overcome the overstiffness of the finite element method (FEM), we present an electromechanical inhomogeneous cell-based smoothed FEM (ISFEM) of FGPMs. Firstly, ISFEM formulations were derived to calculate the transient response of FGPMs, and then, a modified Wilson-θ method was deduced to solve the integration of the FGPM system. The true parameters at the Gaussian integration point in FGPMs were adopted directly to replace the homogenization parameters in an element. ISFEM provides a close-to-exact stiffness of the continuous system, which could automatically and more easily generate for complicated domains and thus significantly decrease numerical errors. The accuracy and trustworthiness of ISFEM were verified as higher than the standard FEM by several numerical examples.


---

## Body

## 1. Introduction

Because of their outstanding electromechanical properties, easy fabricability, and preparation flexibility, piezoelectric materials are extensively applied as sensors and actuators to monitor and modulate the response of structures [1, 2]. Piezoelectric actuators and sensors are innovative for microscopic electromechanical systems and intelligent material systems, particularly in aerospace and medical fields [3]. Conventional piezoelectric sensors and actuators comprise multiple layers of various piezoelectric materials [4–7]. Moreover, piezoelectric layers with uniform material properties are limited by large bending displacement, stress concentration, creeping at high temperature, and failure from interfacial unbounding. All these phenomena are induced by mechanical or electric loading at layer interfaces [8].To overcome the above limitations, Zhu et al. introduced and fabricated functionally graded piezoelectric material (FGPM) sensors and actuators [9, 10]. FGPMs change nonstop in one or more directions without generating internal stress concentration despite the production of large displacements. Takagi et al. fabricated FGPM bimorph actuators by using a mixed system of lead-zirconate-titanate (PZT) and Pt [11]. Nowadays, FGPMs are widely used intelligent materials for sensors and actuators in microstructural engineering. Many efforts have been made to analyze the behaviors and static/dynamic responses of FGPMs (e.g., shells, beams, and plates), such as the wave propagation study of FGPM plates based on the laminate theory [12]. Moreover, exact 3D analysis of FGPM rectangular plates was conducted by using a state-space approach [13] and to investigate the natural frequencies and mode shapes after being poled perpendicular to the middle plane [14]. The above method was also applied to explore the free vibration of rectangular FGPM plates [15]. The semianalytical finite element method (FEM) was used to investigate the static response of anisotropic and linear functionally graded magneto-electro-elastic plates [16]. Lezgy-Nazargah et al. [17] carried out static and dynamic analyses on piezoelectric beams by using a refined sinus model. The results have been found in good agreement with the reference solutions for various electrical and mechanical constrained conditions. Meanwhile, the 3D exact state-space solution [18] and Peano series solution [19] were developed for the cylindrical bending vibration of the FGPM laminates, respectively. Layerwise FEM was adopted to investigate the displacement and stress responses of an FGPM bimorph actuator [20]. Qiu et al. inhibited the vibration of a smart flexible clamped plate by using piezoelectric ceramic patch sensors and actuators [21]. The Timoshenko beam theory was used to analyze the static and dynamic responses of FGPM actuators to thermo-electro-mechanical loading [22]. The first-order shear deformation theory was used to study the static bending, free vibration, and dynamic responses of FGPM plates under electromechanical loading [23]. The free and induced vibrations of FGPM beams under thermo-electro-mechanical loading were characterized using the 3-order shear deformation beam theory [24]. A high-order theory for FGPM shells was proposed based on the generalized Hamilton’s principle [8]. Although the FEM (h-version) is adequate for low-frequency vibration analysis, it is not well suited to the vibration analysis of medium- or high-frequency regimes [25]. The spectral finite element method (SFEM) [26, 27] and the weak form quadrature method (QEM) [28, 29] are developed for the dynamic analysis of FGPM beams and structures.Though FEM is the most widely used and effective numerical approach in practical issues in research and engineering (including mechanics of vibration), it is not necessarily fully perfect or cannot be further improved. For example, the probable overestimation of stiffness in solid structures may lead to locking behavior and inaccurate stress-solving [30]. By adding strain smoothing into FEM[31], Liu et al. established a series of cell-based [32–36], node-based [37, 38], edge-based [39–41], or face-based [42, 43] smoothed FEMs (S-FEMs) and their combinations [44–47]. These S-FEMs with different properties can be used to get desired solutions for a variety of benchmarks and practical mechanic issues [48–50]. The strain-smoothing operations can reduce or alleviate the overstiffness of standard FEM, significantly improving the accuracy of both primal and dual quantities [51]. Moreover, owing to absence of parametric mapping, the shape function derivatives and S-FEM models established in elasticity are not required to be insensitive to mesh distortion [52]. S-FEMs have been successfully extended to analyze the dynamic control of piezoelectric sensors and actuators, topological optimization of linear piezoelectric micromotors, statics, frequency, or defects of smart materials [53–63]. Zheng et al. [64] utilized the cell-based smoothed finite element method with the asymptotic homogenization method to analyze the dynamic issues on micromechanics of piezoelectric composite materials. Zhou et al. [65, 66] deduced the linear and nonlinear cell-based smoothed finite element method of functionally graded magneto-electro-elastic (MEE) structures and further examined the transient responses of MEE sensors or energy harvest structures considering the damping factors. However, there is little literature reported concerning the dynamic response of FGPMs using the electromechanical inhomogeneous cell-based smoothed finite element method. Because of versatility, S-FEMs become convenient and efficient numerical approaches to address different physical issues.Given the continuous change of the gradient of material properties along the thicknessx3 direction and with cell-based gradient smoothing, we deduced the basic formula of ISFEM and a modified Wilson-θ method to solve the integral solution of the FGPM system. The displacements and potentials of FGPM cantilever beams under sine wave load, cosine wave load, step wave load, and triangular wave load were analyzed in comparison with FEM.


## 2. Basic Equations for Piezoelectric Materials

### 2.1. Geometry and Coordinate System

Each beam has a rectangular uniform cross section and is made ofNi layers either completely or partially composed of FGPM beams. The Cartesian coordinate system (x1, x2, x3) and geometric parameters are illustrated in Figure 1.Figure 1
FGPM beams: Cartesian coordinate system and geometric parameters: length (L), width (b), and height (h).


### 2.2. Constitutive Equations

At thek-th layer, 3D linear constitutive equations are polarized along its global coordinates as follows:(1)σ11σ22σ33σ23σ13σ12k=c11c12c13000c12c22c23000c13c23c33000000c44000000c55000000c66kε11ε22ε332ε232ε132ε12k−00e3100e3200e330e240e1500000kE1E2E3k,(2)D1D2D3k=0000e150000e2400e31e32e33000kε11ε22ε332ε232ε132ε12k+χ11000χ22000χ33kE1E2E3k,where σij and εkl are the stress tensor and infinitesimal strain tensor, respectively, Ei and Di are electric field and electric displacement vector components, respectively; eik, ckl, and Xij are the piezoelectric, elastic, and dielectric material constants, respectively. Different from homogeneous piezoelectric materials, the three constants are dependent on coordinate x3. We assume that the material properties along the thickness direction are arbitrarily distributed as follows:(3)cijklx3=cijkl0fx3,eiklx3=eikl0fx3,χikx3=χik0fx3,i,j,k,l==1,2,3,where fx3 is an arbitrary function and ckl0, ekl0, and χkl0 are values at the plane x3 = 0.In an FGPM beam, equations (1) and (2) are reduced to(4)σ¯=C¯ε¯−e¯E¯,D¯=e¯Tε¯+χ¯E¯,where(5)σ¯=σ11σ33σ13,C¯=c11c130c31c33000c55,e¯=0e310e33e150,E¯=E1E3,ε¯=ε11ε33ε13,D¯=D1D3,χ¯=χ1100χ33.


### 2.3. Weak Formulation

The principle of virtual work for a piezoelectric medium of volume Ω and regular boundary surface Γ can be written as(6)δΠ=δU−δW=−∫Ωδε¯Tσ¯dΩ+∫ΓδuTFsdΓ+∫ΩδuTFvdΩ−∫ΩρδuTu¨dΩ+∫ΩδE¯TD¯dΩ+∫ΓQ¯δφdΓ−∫Ωq¯δφdΩ=0,where Fs, Fv, u, and φ are the vectors of surface force, mechanical body force, node displacements, and node electrical potentials, respectively; q¯, Q¯, and ρ are the electrical body charge, surface charge, and mass density, respectively; and δ is the virtual quantity.


## 2.1. Geometry and Coordinate System

Each beam has a rectangular uniform cross section and is made ofNi layers either completely or partially composed of FGPM beams. The Cartesian coordinate system (x1, x2, x3) and geometric parameters are illustrated in Figure 1.Figure 1
FGPM beams: Cartesian coordinate system and geometric parameters: length (L), width (b), and height (h).


## 2.2. Constitutive Equations

At thek-th layer, 3D linear constitutive equations are polarized along its global coordinates as follows:(1)σ11σ22σ33σ23σ13σ12k=c11c12c13000c12c22c23000c13c23c33000000c44000000c55000000c66kε11ε22ε332ε232ε132ε12k−00e3100e3200e330e240e1500000kE1E2E3k,(2)D1D2D3k=0000e150000e2400e31e32e33000kε11ε22ε332ε232ε132ε12k+χ11000χ22000χ33kE1E2E3k,where σij and εkl are the stress tensor and infinitesimal strain tensor, respectively, Ei and Di are electric field and electric displacement vector components, respectively; eik, ckl, and Xij are the piezoelectric, elastic, and dielectric material constants, respectively. Different from homogeneous piezoelectric materials, the three constants are dependent on coordinate x3. We assume that the material properties along the thickness direction are arbitrarily distributed as follows:(3)cijklx3=cijkl0fx3,eiklx3=eikl0fx3,χikx3=χik0fx3,i,j,k,l==1,2,3,where fx3 is an arbitrary function and ckl0, ekl0, and χkl0 are values at the plane x3 = 0.In an FGPM beam, equations (1) and (2) are reduced to(4)σ¯=C¯ε¯−e¯E¯,D¯=e¯Tε¯+χ¯E¯,where(5)σ¯=σ11σ33σ13,C¯=c11c130c31c33000c55,e¯=0e310e33e150,E¯=E1E3,ε¯=ε11ε33ε13,D¯=D1D3,χ¯=χ1100χ33.


## 2.3. Weak Formulation

The principle of virtual work for a piezoelectric medium of volume Ω and regular boundary surface Γ can be written as(6)δΠ=δU−δW=−∫Ωδε¯Tσ¯dΩ+∫ΓδuTFsdΓ+∫ΩδuTFvdΩ−∫ΩρδuTu¨dΩ+∫ΩδE¯TD¯dΩ+∫ΓQ¯δφdΓ−∫Ωq¯δφdΩ=0,where Fs, Fv, u, and φ are the vectors of surface force, mechanical body force, node displacements, and node electrical potentials, respectively; q¯, Q¯, and ρ are the electrical body charge, surface charge, and mass density, respectively; and δ is the virtual quantity.


## 3. Electromechanical ISFEM

The solving domain Ω is discretized intonp elements, which contain Nn nodes; the approximation displacement u¯ and the approximation electrical potential φ¯ for the FGPM problem can be expressed as(7)u¯=∑i=1npNiuui=Nuu,φ¯=∑i=1npNiφφi=Nφφ,where Nu and Nφ are the ISFEM displacement shape function and electrical potential shape function, respectively.Four-node element is divided into four smoothing subdomains. Field nodes, edge smoothing nodes, center smoothing nodes and edge Gaussian points, the outer normal vector distribution, and the shape function values are shown in Figure2.Figure 2
Smoothing subcells and the values of shape functions (reproduced from Zheng et al. [64] under the Creative Commons Attribution License/public domain).At any pointxk in the smoothing subdomains Ωik, the smoothed strain ε¯xk and the smoothed electric field E¯xk are(8)ε¯xk=∫ΩikεxΦx−xkdΩ,(9)E¯xk=∫ΩikExΦx−xkdΩ,where εx and Ex are the strain and electric field in FEM, respectively, and Φx−xk is the constant function:(10)Φx−xk=1Aik,x∈Ωik,0,x∉Ωik,where(11)Aik=∫ΩikdΩ.Substituting equation (10) into equations (8) and (9), then we have(12)ε¯xk=1Aik∫ΓiknukudΓ,(13)E¯xk=1Aik∫ΓiknφkφdΓ,where Γik is the boundary of Ωik and nuk and nφk are the outer normal vector matrixes of the smoothing domain boundary(14)nuk=nx1k00nx3knx3knx1k,nφk=nx1knx3k.Equations (12) and (13) can be rewritten as(15)ε¯xk=∑i=1neB¯uixkui,E¯xk=−∑i=1neB¯φixkφi,where ne is the number of smoothing elements.(16)B¯uixk=1Aik∫ΓkNiunx1k00Niunx3kNiunx3kNiunx1kdΓ,(17)B¯φixk=1Aik∫ΓkNiφnx1kNiφnx3kdΓ.At the Gaussian pointxbG, equations (16) and (17) are (18)B¯uixk=1Aik∑b=1nbNiuxbGnx1k00NiuxbGnx3kNiuxbGnx3kNiuxbGnx1klbk,B¯φixk=1Aik∑b=1nbNiφxbGnx1kNiφxbGnx3klbk,where xbG and lbk are the Gaussian point and the length of the smoothing boundary, respectively, and nb is the total number of boundaries for each smoothing subdomain. As the shape function is linearly changed along each side of the smoothing subdomain, one Gauss point is sufficient for accurate boundary integration [30].The essential distinction between ISFEM and FEM is that FEM needs to construct the shape function matrix of the element, while ISFEM only needs to use the shape function at the Gaussian point of the smoothing element boundary and does not require to involve the shape function derivatives. The above can reduce the continuity requirement of the shape function, and therefore, the accuracy and convergence of the method are improved.The dynamic model of the FGPM electromechanical system can be derived from the Hamilton principle in the following form:(19)M¯q¨+K¯q=F¯,where(20)M¯=Muu000,q¨=u¨φ¨,K¯=KuuKuφKuφTKφφ,q=uφ,F¯=FQ,Muu=∑eMuue,Muue=diagm1,m1,m2,m2,m3,m3,m4,m4,where mi=piT, Aik (i = 1, 2, 3, 4) is the mass of the i-th smoothing element corresponding to node i, T is the smoothing element thickness, and pi is the density of Gaussian integration point of the i-th smoothing subdomain:(21)Kuu=∑i=1ncB¯uiTC¯B¯uiAik,Kuφ=∑i=1ncB¯uiTe¯B¯φiAik,Kφφ=−∑i=1ncB¯φiTχ¯B¯φiAik,F=∫ΩNuTFvdΩ−∫ΓqNuTFsdΓ,Q=∫ΩNφTQ¯dΩ+∫ΓqNφTq¯dΓ,where nc=np×ne.The application of the inhomogeneous smoothing element is to calculate stiffness matrix of the element. The parameters of four smoothing subdomainsAik (i = 1, 2, 3, 4) are various in the elements, so the actual parameters at the Gaussian integration point are taken directly in order to reflect the changes of material property in each element.


## 4. Modified Wilson-θ Method

The modified Wilson-θ method is an important scheme and an implicit integral way to solve the dynamic system equations [63]. If θ > 1.37, the solution is unconditionally stable. The detailed procedures are showed as follows:


### 4.1. Initial Calculation

(1)
Formulate generalized stiffness matrixK¯, mass matrix M¯, and damping matrix C¯
(2)
Calculate initial values ofq, q˙, q¨
(3)
Select the time step Δt and the integral constant θ (θ = 1.4)(22)a0=6θΔt2,a1=3θΔt,a2=2a1,a3=θΔt2,a4=a0θ,a5=−a2θ,a6=1−3θ,a7=Δt2,a8=Δt26.
(4)
Formulate an effective generalized stiffness matrixK˜: K˜=K¯+a0M¯+a1C¯


### 4.2. For Each Time Step

(1)
Calculate the payload at timet + θΔt:(23)F˜t+θΔt∗=F¯t+θF¯t+Δt−F¯t+M¯a0qt+a2q˙t+2q¨t+C¯a1qt+2q˙t+a3q¨t.
(2)
Calculate the generalized displacement at timet + θΔt:(24)K˜qt+θΔt=F˜t+θΔt∗.
(3)
Calculate the generalized acceleration, generalized speed, and generalized displacement at timet + Δt:(25)q¨t+Δt=a4qt+θΔt−qt+a5q˙t+a6q¨t,q˙t+Δt=q˙t+a7q¨t+Δt+q¨t,qt+Δt=qt+Δtq˙t+a8q¨t+Δt+2q¨t.


## 4.1. Initial Calculation

(1)
Formulate generalized stiffness matrixK¯, mass matrix M¯, and damping matrix C¯
(2)
Calculate initial values ofq, q˙, q¨
(3)
Select the time step Δt and the integral constant θ (θ = 1.4)(22)a0=6θΔt2,a1=3θΔt,a2=2a1,a3=θΔt2,a4=a0θ,a5=−a2θ,a6=1−3θ,a7=Δt2,a8=Δt26.
(4)
Formulate an effective generalized stiffness matrixK˜: K˜=K¯+a0M¯+a1C¯


## 4.2. For Each Time Step

(1)
Calculate the payload at timet + θΔt:(23)F˜t+θΔt∗=F¯t+θF¯t+Δt−F¯t+M¯a0qt+a2q˙t+2q¨t+C¯a1qt+2q˙t+a3q¨t.
(2)
Calculate the generalized displacement at timet + θΔt:(24)K˜qt+θΔt=F˜t+θΔt∗.
(3)
Calculate the generalized acceleration, generalized speed, and generalized displacement at timet + Δt:(25)q¨t+Δt=a4qt+θΔt−qt+a5q˙t+a6q¨t,q˙t+Δt=q˙t+a7q¨t+Δt+q¨t,qt+Δt=qt+Δtq˙t+a8q¨t+Δt+2q¨t.


## 5. Numerical Examples

Four numerical examples were conducted under sine wave load, cosine wave load, step wave load and triangular wave load, respectively. FGPM cantilever beams of the same dimensions (lengthL = 40 mm, width h = 5 mm and thickness b = 1 mm) were subjected to forced vibration (Figure 3). The material constants are shown in Table 1. And initial conditions were q=0 and q˙=0 at t = 0 moment. The FGPM beams were made of PZT-4 or PZT-5H on basis of exponentially graded piezoelectric materials with the following material properties:(26)cklx3=ckl0fx3,eklx3=ekl0fx3,χklx3=χkl0fx3,where fx3=enx3/h and n is the gradient parameter.Figure 3
Geometry of the cantilever beam.Table 1
Material constants.
Material
c
11 (GPa)
c
12 (GPa)
c
13 (GPa)
c
33 (GPa)
c
44 (GPa)
ρ (kg/m³)
e
31 (C/m2)
e
33 (C/m2)
e
15 (C/m2)
χ
11/χ0
χ
33/χ0
χ
0 (C/N·m2)
PZT-4
139.0
77.80
74.30
115.0
25.60
7600
−5.20
15.10
12.70
1476.0
1301.0
8.854 × 10−12
PZT-5H
126.0
79.60
84.10
114.0
23.30
7500
−6.50
17.44
23.30
1468.9
1698.3
8.854 × 10−12


### 5.1. Sine Wave Load

The loadF applied to the free end and the load waveform is demonstrated in Figure 4. A convergence investigation with respect to meshes was first carried out. Four smoothing subcells were used for electromechanical ISFEM with ∆t = 1 × 10−3 s. The variations of displacement u3 and electrical potential φ at the loading point of the PZT-4-based FGPM beam combined with respect to time are shown in Figure 5. The results at n = −5 with the element number of 480, 800, 1200, or 1680 were compared with the reference solution [67]. The variations of u3 and φ at the loading point combined with respect to time at n = −1, 0, 1, and 5 in comparison with the reference solution are shown in Figure 6 [67]. Figure 7 illustrates the total energy norm Err versus the mesh density at t = 0.002 s and t = 0.01 s. The simulation results are well consistent among different numbers of meshes, which demonstrate the high convergence of ISFEM.Figure 4
Sine wave load.Figure 5
Variations of the displacementu3 and electrical potential φ at the loading point with time when the gradient parameter was n = −5. (a) Displacement. (b) Electrical potential.
(a)
(b)Figure 6
Variations of the displacementu3 and electrical potential φ at the loading point with time when the gradient parameter was n = −1, 0, 1, and 5. (a) Displacement (n = −1). (b) Electrical potential (n = −1). (c) Displacement (n = 0). (d) Electrical potential (n = 0). (e) Displacement (n = 1). (f) Electrical potential (n = 1). (g) Displacement (n = 5). (h) Electrical potential (n = 5).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 7
Convergence rate of ISFEM.Figure8 shows a comparison of calculation time between ISFEM and FEM at Intel ® Xeon ® CPU E3-1220 v3 @ 3.10 GHz, 16GB RAM. The bandwidth of the system matrices for the FEM and ISFEM is identical. However, in the ISFEM, only the values of shape functions (not the derivatives) at the quadrature points are needed and the requirement of traditional coordinate transform procedure is not necessary to perform the numerical integration. Therefore, the ISFEM generally needs less computational cost than the FEM for handling the dynamic analysis problems.Figure 8
Comparison of computational efficiency.The variations ofu3 and φ at the loading point combined with t = 0.0001 s, 0.0004 s, 0.0025 s, and 0.004 s in the PZT-4-based FGPM cantilever beam are shown in Table 2. The 80 × 10 meshes of ISFEM at n = −5, −1, 0, 1, and 5 are shown in Figure 9, and FEM is considered 160 × 20 elements. Clearly, the results of ISFEM with 80 × 10 elements are the same as the calculated results of FEM using 160 × 20 elements, suggesting ISFEM has higher accuracy.Table 2
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-4.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−10 m
12.852
12.658
2.760
2.751
1.707
1.702
1.006
1.004
0.0756
0.0754
φ × 10−2 V
11.848
11.912
0.442
0.427
0.212
0.207
0.0961
0.0944
0.00112
0.00110
0.0004
u
3 × 10−9 m
28.517
28.367
4.735
4.715
2.855
2.842
1.738
1.731
0.187
0.187
φ × 10−1 V
5.851
5.747
0.597
0.582
0.349
0.339
0.193
0.188
0.0805
0.782
0.0025
u
3 × 10−8 m
13.345
13.191
1.538
1.528
0.917
0.913
0.564
0.561
0.0882
0.0873
φ (V)
2.447
2.391
0.198
0.191
0.112
0.109
0.0621
0.0605
0.00386
0.00372
0.004
U
3 × 10−9 m
77.799
77.682
10.154
10.146
6.165
6.123
3.727
3.725
0.514
0.516
φ (V)
1.434
1.406
0.129
0.125
0.0756
0.0740
0.0414
0.0401
0.00224
0.00221Figure 9
Mesh of the FGPM beam.The variations ofu3 and φ at the loading point combined with time of the PZT-5H-based FGPM cantilever beam are listed in Table 3. Clearly, when n changes from −5 to 5, the maximum u3 and φ decrease, which is consistent with the PZT-4-based FGPM cantilever beam. Furthermore, the results of ISFEM with 80 × 10 elements are the same as the calculated results of FEM using 160 × 20 elements, suggesting ISFEM has higher accuracy.Table 3
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-5H.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−10 m
13.691
13.424
2.951
2.942
1.826
1.821
1.076
1.074
0.0802
0.0800
φ × 10−2 V
12.098
12.493
0.413
0.403
0.190
0.194
0.0825
0.0843
0.000624
0.000649
0.0004
u
3 × 10−9 m
31.394
31.196
5.583
5.561
3.381
3.367
2.050
2.042
0.207
0.206
φ × 10−1 V
5.901
6.057
0.599
0.607
0.349
0.350
0.192
0.190
0.0736
0.0735
0.0025
u
3 × 10−8 m
16.122
16.006
2.017
2.003
1.118
1.117
0.741
0.736
0.107
0.106
φ (V)
2.522
2.584
0.219
0.219
0.122
0.123
0.0691
0.0689
0.00392
0.00393
0.004
u
3 × 10−9 m
76.019
75.052
10.339
10.264
6.264
6.261
3.795
3.768
0.501
0.496
φ (V)
1.402
1.437
0.113
0.113
0.0649
0.0653
0.0351
0.0350
0.00179
0.00178


### 5.2. Cosine Wave Load

The cosine loadF applied to the free end and load waveform is shown in Figure 10. The variations of u3 and φ at the loading point combined with time of PZT-4- and PZT-5H-based FGPM cantilever beams are listed in Tables 4 and 5, respectively. The calculated results of ISFEM with 80 × 10 elements are the same as those of FEM using 160 × 20 elements, implying that ISFEM possesses higher accuracy.Figure 10
Cosine wave load.Table 4
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-4 under cosine load.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−9 m
20.427
20.159
4.387
4.382
2.713
2.710
1.600
1.599
0.120
0.120
φ × 10−2 V
188.326
189.821
7.029
6.896
3.369
3.269
1.528
1.5044
0.0178
0.0175
0.0004
u
3 × 10−8 m
21.208
21.234
2.738
2.743
1.610
1.611
1.006
1.004
0.141
0.142
φ × 10−1 V
27.889
27.949
3.315
3.216
0.193
0.187
1.113
1.098
0.0648
0.0639
0.0025
u
3 × 10−9 m
−33.487
−32.845
−12.058
−11.952
−4.445
−4.345
−4.436
−4.392
−0.225
−0.221
φ (V)
−0.138
−0.137
−0.132
−0.128
−0.0491
−0.0479
−0.0486
−0.0470
−0.00086
−0.00085
0.004
u
3 × 10−9 m
−11.111
−11.501
−0.984
−0.969
−5.549
−5.616
−0.351
−0.339
−0.0981
−0.0913
φ (V)
−1.408
−1.365
−0.0176
−0.0174
−0.0653
−0.0635
0.00234
0.00229
0.000596
0.000579Table 5
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-5H under cosine load.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−9 m
2.902
2.901
4.691
4.685
2.902
2.901
1.710
1.710
0.128
0.127
φ × 10−2 V
192.287
198.967
6.561
6.742
3.017
3.092
1.312
1.342
0.00992
0.0101
0.0004
u
3 × 10−8 m
23.941
23.948
3.486
3.493
2.069
2.064
1.281
1.284
0.160
0.159
φ × 10−1 V
28.178
29.094
3.5841
3.613
2.108
2.121
1.203
1.201
0.0621
0.0616
0.0025
u
3 × 10−9 m
117.958
122.357
−4.831
−4.723
−7.090
−6.942
−1.777
−1.715
0.795
0.8136
φ (V)
0.50863
0.5017
−0.0407
−0.0393
−0.0688
−0.0652
−0.0151
−0.0148
0.00338
0.00342
0.004
u
3 × 10−9 m
−183.26
−179.45
−27.847
−27.063
−9.636
−9.387
−10.233
−9.998
−1.287
−1.263
φ (V)
−2.273
−2.176
−0.284
−0.277
−0.0931
−0.0885
−0.0951
−0.0949
−0.00466
−0.00451


### 5.3. Step Wave Load

The step loadF applied to the free end and load waveform is indicated in Figure 11. The variations of u3 and φ at the loading point combined with time of PZT-4- and PZT-5H-based FGPM cantilever beams are illustrated in Figures 12 and 13, respectively. It is clearly shown that ISFEM possesses higher accuracy than FEM for the calculated results of ISFEM using 80 × 10 elements are the same as FEM using 160 × 20 elements.Figure 11
Step wave load.Figure 12
Variations of displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-4 when the gradient parameter was n = −5, 0 and 5 under step load. (a) Displacement. (b) Electrical potential.
(a)
(b)Figure 13
Variations of displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-5H when the gradient parameter was n = −5, 0 and 5 under step load. (a) Displacement. (b) Electrical potential.
(a)
(b)


### 5.4. Triangular Wave Load

The triangular loadF applied to the free end and load waveform is shown in Figure 14. The variations of u3 and φ at the loading point combined with time of PZT-4- and PZT-5H-based FGPM cantilever beams are shown in Figures 15 and 16, respectively. It shows that the solutions of CS-FEM with less elements are the same as the solutions of FEM using more elements.Figure 14
Triangular wave load.Figure 15
Variations of displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-4 when the gradient parameter was n = −5, 0, and 5 under triangular load. (a) Displacement. (b) Electrical potential.
(a)
(b)Figure 16
Variations of the displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-5H when the gradient parameter was n = −5, 0, and 5 under triangular load. (a) Displacement. (b) Electrical potential.
(a)
(b)


## 5.1. Sine Wave Load

The loadF applied to the free end and the load waveform is demonstrated in Figure 4. A convergence investigation with respect to meshes was first carried out. Four smoothing subcells were used for electromechanical ISFEM with ∆t = 1 × 10−3 s. The variations of displacement u3 and electrical potential φ at the loading point of the PZT-4-based FGPM beam combined with respect to time are shown in Figure 5. The results at n = −5 with the element number of 480, 800, 1200, or 1680 were compared with the reference solution [67]. The variations of u3 and φ at the loading point combined with respect to time at n = −1, 0, 1, and 5 in comparison with the reference solution are shown in Figure 6 [67]. Figure 7 illustrates the total energy norm Err versus the mesh density at t = 0.002 s and t = 0.01 s. The simulation results are well consistent among different numbers of meshes, which demonstrate the high convergence of ISFEM.Figure 4
Sine wave load.Figure 5
Variations of the displacementu3 and electrical potential φ at the loading point with time when the gradient parameter was n = −5. (a) Displacement. (b) Electrical potential.
(a)
(b)Figure 6
Variations of the displacementu3 and electrical potential φ at the loading point with time when the gradient parameter was n = −1, 0, 1, and 5. (a) Displacement (n = −1). (b) Electrical potential (n = −1). (c) Displacement (n = 0). (d) Electrical potential (n = 0). (e) Displacement (n = 1). (f) Electrical potential (n = 1). (g) Displacement (n = 5). (h) Electrical potential (n = 5).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)Figure 7
Convergence rate of ISFEM.Figure8 shows a comparison of calculation time between ISFEM and FEM at Intel ® Xeon ® CPU E3-1220 v3 @ 3.10 GHz, 16GB RAM. The bandwidth of the system matrices for the FEM and ISFEM is identical. However, in the ISFEM, only the values of shape functions (not the derivatives) at the quadrature points are needed and the requirement of traditional coordinate transform procedure is not necessary to perform the numerical integration. Therefore, the ISFEM generally needs less computational cost than the FEM for handling the dynamic analysis problems.Figure 8
Comparison of computational efficiency.The variations ofu3 and φ at the loading point combined with t = 0.0001 s, 0.0004 s, 0.0025 s, and 0.004 s in the PZT-4-based FGPM cantilever beam are shown in Table 2. The 80 × 10 meshes of ISFEM at n = −5, −1, 0, 1, and 5 are shown in Figure 9, and FEM is considered 160 × 20 elements. Clearly, the results of ISFEM with 80 × 10 elements are the same as the calculated results of FEM using 160 × 20 elements, suggesting ISFEM has higher accuracy.Table 2
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-4.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−10 m
12.852
12.658
2.760
2.751
1.707
1.702
1.006
1.004
0.0756
0.0754
φ × 10−2 V
11.848
11.912
0.442
0.427
0.212
0.207
0.0961
0.0944
0.00112
0.00110
0.0004
u
3 × 10−9 m
28.517
28.367
4.735
4.715
2.855
2.842
1.738
1.731
0.187
0.187
φ × 10−1 V
5.851
5.747
0.597
0.582
0.349
0.339
0.193
0.188
0.0805
0.782
0.0025
u
3 × 10−8 m
13.345
13.191
1.538
1.528
0.917
0.913
0.564
0.561
0.0882
0.0873
φ (V)
2.447
2.391
0.198
0.191
0.112
0.109
0.0621
0.0605
0.00386
0.00372
0.004
U
3 × 10−9 m
77.799
77.682
10.154
10.146
6.165
6.123
3.727
3.725
0.514
0.516
φ (V)
1.434
1.406
0.129
0.125
0.0756
0.0740
0.0414
0.0401
0.00224
0.00221Figure 9
Mesh of the FGPM beam.The variations ofu3 and φ at the loading point combined with time of the PZT-5H-based FGPM cantilever beam are listed in Table 3. Clearly, when n changes from −5 to 5, the maximum u3 and φ decrease, which is consistent with the PZT-4-based FGPM cantilever beam. Furthermore, the results of ISFEM with 80 × 10 elements are the same as the calculated results of FEM using 160 × 20 elements, suggesting ISFEM has higher accuracy.Table 3
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-5H.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−10 m
13.691
13.424
2.951
2.942
1.826
1.821
1.076
1.074
0.0802
0.0800
φ × 10−2 V
12.098
12.493
0.413
0.403
0.190
0.194
0.0825
0.0843
0.000624
0.000649
0.0004
u
3 × 10−9 m
31.394
31.196
5.583
5.561
3.381
3.367
2.050
2.042
0.207
0.206
φ × 10−1 V
5.901
6.057
0.599
0.607
0.349
0.350
0.192
0.190
0.0736
0.0735
0.0025
u
3 × 10−8 m
16.122
16.006
2.017
2.003
1.118
1.117
0.741
0.736
0.107
0.106
φ (V)
2.522
2.584
0.219
0.219
0.122
0.123
0.0691
0.0689
0.00392
0.00393
0.004
u
3 × 10−9 m
76.019
75.052
10.339
10.264
6.264
6.261
3.795
3.768
0.501
0.496
φ (V)
1.402
1.437
0.113
0.113
0.0649
0.0653
0.0351
0.0350
0.00179
0.00178


## 5.2. Cosine Wave Load

The cosine loadF applied to the free end and load waveform is shown in Figure 10. The variations of u3 and φ at the loading point combined with time of PZT-4- and PZT-5H-based FGPM cantilever beams are listed in Tables 4 and 5, respectively. The calculated results of ISFEM with 80 × 10 elements are the same as those of FEM using 160 × 20 elements, implying that ISFEM possesses higher accuracy.Figure 10
Cosine wave load.Table 4
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-4 under cosine load.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−9 m
20.427
20.159
4.387
4.382
2.713
2.710
1.600
1.599
0.120
0.120
φ × 10−2 V
188.326
189.821
7.029
6.896
3.369
3.269
1.528
1.5044
0.0178
0.0175
0.0004
u
3 × 10−8 m
21.208
21.234
2.738
2.743
1.610
1.611
1.006
1.004
0.141
0.142
φ × 10−1 V
27.889
27.949
3.315
3.216
0.193
0.187
1.113
1.098
0.0648
0.0639
0.0025
u
3 × 10−9 m
−33.487
−32.845
−12.058
−11.952
−4.445
−4.345
−4.436
−4.392
−0.225
−0.221
φ (V)
−0.138
−0.137
−0.132
−0.128
−0.0491
−0.0479
−0.0486
−0.0470
−0.00086
−0.00085
0.004
u
3 × 10−9 m
−11.111
−11.501
−0.984
−0.969
−5.549
−5.616
−0.351
−0.339
−0.0981
−0.0913
φ (V)
−1.408
−1.365
−0.0176
−0.0174
−0.0653
−0.0635
0.00234
0.00229
0.000596
0.000579Table 5
Displacementu3 and electrical potential φ at the loading point with the gradient parameter n = −5, −1, 0, 1, and 5 of the FGPM beam based on PZT-5H under cosine load.
Time (s)
Field variables
n = −5
n = −1
n = 0
n = 1
n = 5
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
ISFEM
FEM
0.0001
u
3 × 10−9 m
2.902
2.901
4.691
4.685
2.902
2.901
1.710
1.710
0.128
0.127
φ × 10−2 V
192.287
198.967
6.561
6.742
3.017
3.092
1.312
1.342
0.00992
0.0101
0.0004
u
3 × 10−8 m
23.941
23.948
3.486
3.493
2.069
2.064
1.281
1.284
0.160
0.159
φ × 10−1 V
28.178
29.094
3.5841
3.613
2.108
2.121
1.203
1.201
0.0621
0.0616
0.0025
u
3 × 10−9 m
117.958
122.357
−4.831
−4.723
−7.090
−6.942
−1.777
−1.715
0.795
0.8136
φ (V)
0.50863
0.5017
−0.0407
−0.0393
−0.0688
−0.0652
−0.0151
−0.0148
0.00338
0.00342
0.004
u
3 × 10−9 m
−183.26
−179.45
−27.847
−27.063
−9.636
−9.387
−10.233
−9.998
−1.287
−1.263
φ (V)
−2.273
−2.176
−0.284
−0.277
−0.0931
−0.0885
−0.0951
−0.0949
−0.00466
−0.00451


## 5.3. Step Wave Load

The step loadF applied to the free end and load waveform is indicated in Figure 11. The variations of u3 and φ at the loading point combined with time of PZT-4- and PZT-5H-based FGPM cantilever beams are illustrated in Figures 12 and 13, respectively. It is clearly shown that ISFEM possesses higher accuracy than FEM for the calculated results of ISFEM using 80 × 10 elements are the same as FEM using 160 × 20 elements.Figure 11
Step wave load.Figure 12
Variations of displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-4 when the gradient parameter was n = −5, 0 and 5 under step load. (a) Displacement. (b) Electrical potential.
(a)
(b)Figure 13
Variations of displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-5H when the gradient parameter was n = −5, 0 and 5 under step load. (a) Displacement. (b) Electrical potential.
(a)
(b)


## 5.4. Triangular Wave Load

The triangular loadF applied to the free end and load waveform is shown in Figure 14. The variations of u3 and φ at the loading point combined with time of PZT-4- and PZT-5H-based FGPM cantilever beams are shown in Figures 15 and 16, respectively. It shows that the solutions of CS-FEM with less elements are the same as the solutions of FEM using more elements.Figure 14
Triangular wave load.Figure 15
Variations of displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-4 when the gradient parameter was n = −5, 0, and 5 under triangular load. (a) Displacement. (b) Electrical potential.
(a)
(b)Figure 16
Variations of the displacementu3 and electrical potential φ at the loading point with time of the FGPM beam based on PZT-5H when the gradient parameter was n = −5, 0, and 5 under triangular load. (a) Displacement. (b) Electrical potential.
(a)
(b)


## 6. Conclusions

An electromechanical ISFEM was proposed given the continuous changes of the gradient of material properties along the thicknessx3 direction and with cell-based gradient smoothing. The modified Wilson-θ method was deduced to solve the integral solution of the FGPM system. The displacements and potentials of cantilever beams combining with sine load, cosine load, step load, and triangular load were analyzed by ISFEM in comparison with FEM.(1)
ISFEM is correct and effective in solving the dynamic response of FGPM structures(2)
ISFEM can reduce the systematic stiffness of FEM and provides calculations closer to the true values(3)
ISFEM is more efficient than FEM and takes less computation time at the same accuracyThis study indicates a possibility to select suitable grading controlled by the power law index according to the application.



---

*Source: 2812748-2019-08-06.xml*