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.

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 [47]. 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 [3236], node-based [37, 38], edge-based [3941], or face-based [42, 43] smoothed FEMs (S-FEMs) and their combinations [4447]. These S-FEMs with different properties can be used to get desired solutions for a variety of benchmarks and practical mechanic issues [4850]. 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 [5363]. 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 thickness x3 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 of Ni layers either completely or partially composed of FGPM beams. The Cartesian coordinate system (x1, x2, x3) and geometric parameters are illustrated in Figure 1.

2.2. Constitutive Equations

At the k-th layer, 3D linear constitutive equations are polarized along its global coordinates as follows:where and are the stress tensor and infinitesimal strain tensor, respectively, and are electric field and electric displacement vector components, respectively; , , and are the piezoelectric, elastic, and dielectric material constants, respectively. Different from homogeneous piezoelectric materials, the three constants are dependent on coordinate . We assume that the material properties along the thickness direction are arbitrarily distributed as follows:where is an arbitrary function and , , and are values at the plane x3 = 0.

In an FGPM beam, equations (1) and (2) are reduced towhere

2.3. Weak Formulation

The principle of virtual work for a piezoelectric medium of volume Ω and regular boundary surface Γ can be written aswhere Fs, Fv, u, and φ are the vectors of surface force, mechanical body force, node displacements, and node electrical potentials, respectively; , , 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 into np elements, which contain Nn nodes; the approximation displacement and the approximation electrical potential for the FGPM problem can be expressed aswhere and 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 Figure 2.

At any point in the smoothing subdomains , the smoothed strain and the smoothed electric field arewhere and are the strain and electric field in FEM, respectively, and is the constant function:where

Substituting equation (10) into equations (8) and (9), then we havewhere is the boundary of and and are the outer normal vector matrixes of the smoothing domain boundary

Equations (12) and (13) can be rewritten aswhere is the number of smoothing elements.

At the Gaussian point , equations (16) and (17) are where and are the Gaussian point and the length of the smoothing boundary, respectively, and 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:wherewhere , (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 is the density of Gaussian integration point of the i-th smoothing subdomain:where .

The application of the inhomogeneous smoothing element is to calculate stiffness matrix of the element. The parameters of four smoothing subdomains (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 matrix , mass matrix , and damping matrix (2)Calculate initial values of , , (3)Select the time step Δt and the integral constant θ (θ = 1.4)(4)Formulate an effective generalized stiffness matrix :

4.2. For Each Time Step

(1)Calculate the payload at time t + θΔt:(2)Calculate the generalized displacement at time t + θΔt:(3)Calculate the generalized acceleration, generalized speed, and generalized displacement at time t + Δ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 (length L = 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 and 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:where and n is the gradient parameter.

5.1. Sine Wave Load

The load F 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 8 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.

The variations of u3 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.

The variations of u3 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.

5.2. Cosine Wave Load

The cosine load F 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.

5.3. Step Wave Load

The step load F 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.

5.4. Triangular Wave Load

The triangular load F 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.

6. Conclusions

An electromechanical ISFEM was proposed given the continuous changes of the gradient of material properties along the thickness x3 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 accuracy

This study indicates a possibility to select suitable grading controlled by the power law index according to the application.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare no conflicts of interest.

Authors’ Contributions

Liming Zhou and Bin Cai performed the simulations. Bin Cai contributed to the writing of the manuscript.


Special thanks are due to Professor Guirong Liu for the S-FEM Source Code in http://www.ase.uc.edu/∼liugr/software.html. This work was financially supported by the National Key R&D Program of China (grant no. 2018YFF01012401-06), Jilin Provincial Department of Science and Technology Fund Project (grant no. 20170101043JC), Jilin Provincial Department of Education (grant nos. JJKH20180084KJ and JJKH20170788KJ), Fundamental Research Funds for the Central Universities, Science and Technological Planning Project of Ministry of Housing and Urban–Rural Development of the People's Republic of China (2017-K9-047) and Graduate Innovation Fund of Jilin University (grant no. 101832018C184).