Geometrical Investigation of Piezoelectric Patches for Broadband Energy Harvesting in Non-Deterministic Composite Plates

Abstract

Mechanical energy is the most ubiquitous form of energy that can be harvested and converted into useful electrical power. For this reason, the piezoelectric energy harvesters (PEHs), with their inherent electromechanical coupling and high-power density, have been widely incorporated in many applications to generate power from ambient mechanical vibrations. However, one of the main challenges to the wider adoption of PEHs is how to optimize their design for maximum energy harvesting. In this paper, an investigation was conducted on the energy harvesting from seven piezoelectric patch shapes (differing in the number of edges) when attached to a non-deterministic laminated composite (single/double lamina) plate subjected to change in fiber orientation. The performance of the PEHs was examined through a coupled-field finite element (FE) model. The plate was simply supported, and its dynamics were randomized by attaching randomly distributed point masses on the plate surface in addition to applying randomly located time-harmonic point forces. The randomization of point masses and point force location on a thin plate produce non-deterministic response. The design optimization was performed by employing the ensemble-responses of the electrical potential developed across the electrodes of the piezoelectric patches. The results present the optimal fiber orientation and patch shape for maximum energy harvesting in the case of single and double lamina composite plates. The results show that the performance is optimal at 0° or 90° fiber orientation for single-lamina, and at 0°/0° and 0°/90° fiber orientations for double-lamina composites. For frequencies below 25 Hz, patches with a low number of edges exhibited a higher harvesting performance (triangular for single-lamina/quadrilateral for double-lamina). As for the broadband frequencies (above 25 Hz), the performance was optimal for the patches with a higher number of edges (dodecagonal for single-lamina/octagonal for double-lamina).

1. Introduction

Ambient mechanical vibration is ubiquitous; harvesting this energy to generate electricity has been considered a pivotal point for researchers in fulfilling the desire for renewable green energy [1,2,3]. Various studies have proposed vibration energy harvesting technologies that can efficiently harness ambient vibration energy and provide sustainable power sources. These technologies demonstrate potential use in powering electronics in a variety of applications, such as micro-powered electronic devices, embedded sensors in structures, medical devices, and wireless sensor networks [4,5,6,7,8]. Moreover, the piezoelectric materials have been used as sensors and/or actuators in the form of layers or patches embedded and/or surface bonded on structural plates for vibration control Bodaghi et al. [9] The energy harvesting mechanism based on piezoelectric energy harvesting has received the utmost interest from many researchers. Piezoelectric energy harvesters (PEHs) are the prominent harvesting systems owing to their high energy conversion capabilities and simple structure [10]. The versatility of PEHs has facilitated their incorporation in various energy harvesting techniques such as fluid-based energy harvesting (Truitt and Mahmoodi, 2013) [11]. Due to the growing interest in PEHs, many investigations have been conducted on the performance optimization of PEHs for optimum power generation. These studies were focused on optimizing the performance of various PEH types, such as the cantilever-type PEH, MEMS-based PEH, and sandwich-type PEH.

Park et al. [12] performed an optimization study on a cantilever-type PEH in which the free tip is excited by the rotary motion of a mechanical device and maximized the power output. Wein et al. [13] performed topology optimization of a cantilever-type PEH using stress norm constraints to optimize the electrical power for a sinusoidal vibrational excitation. Zhu et al. [14] conducted a theoretical and experimental analysis to study the effect of geometrical dimension on the energy harvesting performance of a unimorph cantilever-type PEH with fixed resonance frequency. Kim and Le [15] used a finite element analysis-based topology optimization to calculate the topologies of a substrate plate and a piezoelectric layer in a vibrating cantilever plate application. The designs with optimal topologies were proposed for three different piezoelectric materials, and their voltage output was compared to investigate the most suitable material for the maximum harvesting capacity. Sordo et al. [16] proposed optimizing a MEMS-based PEH by maximizing the number of resonant modes within the narrow operating bandwidth. Song et al. [17] proposed an optimization strategy to maximize the power output density of a PVDF-based cantilever-type energy harvester. This strategy showed that the power output mainly depends on the maximum allowable stress within the beam structure and the device’s working frequency, which can be obtained by adjusting the geometry of the piezoelectric layers. Izadgoshasb et al. [18] explored the effect of the orientation of the cantilever beam structure on the power generation efficiency of the PEH. This type of PEH was proved to have higher voltage output than the conventional PEHs due to the selection of materials and geometry of the core and the metal layers. Kim et al. [19] investigated the performance optimization of cantilever-based PEH by adjusting the position of the piezoelectric generators based on the changes in the phase angle of the substrate. Ichige et al. [20] proposed a PEH with mechanical metamaterials for the elastic layer optimized for high power output and low resonance frequency. Peralta et al. [21] investigated the performance of a bimorph cantilever-type PEH subjected to base excitation to optimize its performance. They performed a parametric study to investigate the effect of shape perturbation on the fundamental frequency, amplitude, and frequency response function.

Most of these studies have investigated the performance optimization of the cantilever-type PEHs at the dominant resonant modes, which are predominantly evident at the low-frequency range. However, investigations on the performance of PEHs at the high-frequency ranges are seldom found in research. Modelling the vibration response of structures at high-frequency ranges is more challenging than the low-frequency range [22]. The vibration response can be divided into three main frequency ranges: low, mid, and high-frequency. The response exhibits dominant resonant peaks at the low-frequency range, and the structure response tends to have a high variance. In contrast, the high-frequency peaks are no longer distinct, and the response variance becomes smaller. The mid-frequency range holds the characteristics of both the low and the high-frequency ranges [23].

Structures vibrating at the low-frequency range are subjected to long-wavelength deformations and are referred to as deterministic subsystems (DSs); these structures are insensitive to structural uncertainties and can be analyzed using deterministic modelling techniques such as the finite element analysis (FEA). The structures vibrating at high-frequency ranges are subjected to short-wavelength deformations and are sensitive to uncertainties. Such structures are referred to as non-deterministic subsystems (Non-Ds) and analyzed using statistical modelling techniques such as statistical energy analysis (SEA) [24]. The mid-frequency is the range where some components in the system have distinct modal responses and others that behave statistically all in the same frequency range. Therefore, a single modelling technique may not be sufficient to predict the response of such structures [25]. Although the computational resources for FEA can handle smaller elements to capture short wavelengths, it would not be efficient due to computational expenses. Therefore, at high-frequency ranges, SEA is performed to get the ensemble average and predict the response of such vibrations. Hence, SEA is more favourable to use to analyze the dynamics of high-frequency vibrations. In SEA, the modal overlapping factor ($\mathrm{MOF}$) is used to quantify the degree of overlap in the modal response. MOF is defined as the ratio between the modal bandwidth at the half-amplitude and the average modal spacing, and it can be used to establish the three main frequency ranges: low, mid, and high-frequency. For orthotropic plates, the average modal frequency spacing is one of the parameters used in SEA analysis; which can be analytically defined in terms of mechanical properties of the plate using geometric mean of the wave velocities in two principal directions $x$ and $y$ [26]:

(1)${\overline{\delta f}}_{ortho}=\frac{h}{\sqrt{3}A}\sqrt{c_{Lx}^{’}{c}_{Ly}^{’}}$

(2)$c_{Lx}^{’}=\sqrt{E_1/\rho (1-{\upsilon }_{12}^2}), c_{Ly}^{’}=\sqrt{E_2/\rho \left(1-{\upsilon }_{21}^2\right)}$

where $h$ is the thickness, $A$ is the surface area, while $E_1$ and $E_2$ are modulus of elasticity, ${\upsilon }_{12}$ and ${\upsilon }_{21}$ are Poisson’s ratios in two principal directions, and $\rho$ is the mass density. ($c_{Lx}^{’}$ and $c_{Ly}^{’}$) are the longitudinal wave speed in the $x$ and $y$ directions, respectively (as shown in Figure 1). Equation (1) is only valid for wavelength much greater than the thickness of the plate. The $\mathrm{MOF}$ of flexural modes in thin orthotropic plates can be analytically written as [27]:

(3)$\mathrm{MOF}=\frac{\mathsf{\pi }}{2}\frac{f\eta }{\delta f_{ortho}}$

where $\eta$ is the damping loss factor and $f$ is the frequency. $\mathrm{MOF}$ can be used to indicate the high-frequency threshold value, i.e., $\mathrm{MOF}=2.5$ for plates. The effect of uncertainties due to mass variations, stiffness variations due to boundary conditions, and discontinues has a minimal influence on the low-frequency modes. However, at high-frequency, the structural response is sensitive to such uncertainties. Due to this, the response will be non-deterministic, and the ensemble average is taken. Literature involving energy harvesting have addressed such case; however, the change of the structure anisotropy has not been investigated in the literature.

In this paper, an investigation on the performance of the ceramic-type PEH attached to a non-DS composite plate was conducted using FEA. The investigation focuses on changing the piezoelectric shape and the plate’s anisotropy on the developed voltage. The study proposed seven different piezoelectric patch shapes holding the material properties of the PZT-5A attached to a simply supported laminated composite. Randomly located point masses and time-harmonic point forces were distributed on the host plate to generate uncertainties in the structure, producing the effect of an infinite plate that generates the effect of a non-DS. The investigation was performed by employing the ensemble average response. The electric potential across the piezoelectric patch was obtained, and the ensemble responses were analyzed.

This paper is organized as follows: Section 2 presents the methodology where the mathematical modelling and the FEA model construction are discussed including: composite laminate construction, geometry, material properties, meshing and boundary conditions. Section 3 complies with the ensemble response results obtained from the performed simulations and discusses the response behavior. The paper ends with a summary of the completed work and presents the finding of the research.

2. Methodology

The dynamics of the laminated composite plate with the piezoelectric patch attached to it were simulated using ANSYS Mechanical APDL^® (ANSYS Inc., Canonsburg, PA, USA). In this section, the mathematical modelling and the FE model construction are presented.

2.1. Theory of Laminated Composite Plates

2.1.1. Strain-Displacement Relationship

This section considers a fiber-reinforced laminated composite with distributed piezoelectric patches as sensors bonded on the surface. Figure 1 represents the coordinates of the laminated composite plate. The $x$-$y$ plane coincides with the mid-plane, whereas the $z$-axis is normal to the mid-plane.

The displacements field at any point within the laminated composite is given based on the first-order shear deformation plate theory (FSDT) (see Appendix A) by [28]:

(4)$\begin{array}{c}u\left(x,y,z,t\right)=u_0\left(x,y,t\right)+z{\varphi }_x\left(x,y,z\right)\\ v\left(x,y,z,t\right)=v_0\left(x,y,t\right)+z{\varphi }_y\left(x,y,z\right)\\ w\left(x,y,z,t\right)=w_0\left(x,y,t\right)\end{array}$

where $u$, $v$ and $w$ are the displacements in the $x$, $y$, and $z$ directions, and ($u_0$, $v_0$,$w_0$) denote the displacements on the plane $z=0$. While ${\varphi }_x$ and ${\varphi }_y$ are the rotations of a transverse normal about the $y$- and $x$-axes, respectively (they do not follow the right-hand rule):

(5)$\frac{\partial u}{\partial z}={\varphi }_x, \frac{\partial v}{\partial z}={\varphi }_y$

In FSDT, the straight lines originally normal to the mid-plane will not remain straight during bending. They will rotate with angular displacements ${\varphi }_x$ and ${\varphi }_y$ independent of the transverse displacement and its derivatives. This allows the development of stress and shear strains during bending. Additionally, the transverse shear strain is assumed to be constant across the thickness. For thin-plates, the rotation functions ${\varphi }_x$ and ${\varphi }_y$ approach the respective slopes of the transverse deflections:

(6)${\varphi }_x=-\frac{\partial w_0}{\partial x}, {\varphi }_y=-\frac{\partial w_0}{\partial y}$

The nonlinear strains associated with the displacement field in Equation (1) are [29]:

(7)$\left\{\begin{array}{l}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\gamma }_{yz}\\ {\gamma }_{xz}\\ {\gamma }_{xy}\end{array}\right\}=\left\{\begin{array}{l}{\epsilon }_{xx}^{(0)}\\ {\epsilon }_{yy}^{(0)}\\ {\gamma }_{yz}^{(0)}\\ {\gamma }_{xz}^{(0)}\\ {\gamma }_{xy}^{(0)}\end{array}\right\}+z\left\{\begin{array}{c}{\epsilon }_{xx}^{(1)}\\ {\epsilon }_{yy}^{(1)}\\ {\gamma }_{yz}^{(1)}\\ {\gamma }_{xz}^{(1)}\\ {\gamma }_{xy}^{(1)}\end{array}\right\}=\left\{\begin{array}{c}\frac{\mathrm{\partial }{u}_0}{\mathrm{\partial }x}+\frac{1}{2}{\left(\frac{\mathrm{\partial }{w}_0}{\mathrm{\partial }x}\right)}^2\\ \frac{\mathrm{\partial }{v}_0}{\mathrm{\partial }y}+\frac{1}{2}{\left(\frac{\mathrm{\partial }{w}_0}{\mathrm{\partial }y}\right)}^2\\ \frac{\mathrm{\partial }{w}_0}{\mathrm{\partial }y}+{\varphi }_y\\ \frac{\mathrm{\partial }{w}_0}{\mathrm{\partial }x}+{\varphi }_x\\ \frac{\mathrm{\partial }{u}_0}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{v}_0}{\mathrm{\partial }x}\frac{\mathrm{\partial }{w}_0}{\mathrm{\partial }x}\frac{\mathrm{\partial }{w}_0}{\mathrm{\partial }y}\end{array}\right\}+z\left\{\begin{array}{c}\frac{\mathrm{\partial }{\mathrm{\varphi }}_x}{\mathrm{\partial }x}\\ \frac{\mathrm{\partial }{\mathrm{\varphi }}_y}{\mathrm{\partial }y}\\ 0\\ 0\\ \frac{\mathrm{\partial }{\mathrm{\varphi }}_x}{\mathrm{\partial }y}+\frac{\mathrm{\partial }{\mathrm{\varphi }}_y}{\mathrm{\partial }x}\end{array}\right\}$

where (${\epsilon }_{xx}^{\left(0\right)}$,${\epsilon }_{yy}^{\left(0\right)}$,${\gamma }_{yz}^{\left(0\right)}$,${\gamma }_{xz}^{\left(0\right)}$,${\gamma }_{xy}^{\left(0\right)}$) are the membrane strains, and (${\epsilon }_{xx}^{\left(1\right)}$,${\epsilon }_{yy}^{\left(1\right)}$,${\gamma }_{yz}^{\left(1\right)}$,${\gamma }_{xz}^{\left(1\right)}$,${\gamma }_{xy}^{\left(1\right)}$) are the flexural (bending) strains.

2.1.2. Constitutive Equations of a Lamina

In formulating the constitutive equations that represent the mechanical behaviour of a typical lamina of a composite laminate, two assumptions are taken into consideration: (1) no gaps or empty regions are present in the lamina (i.e., lamina is a continuum), and (2) the lamina has a linear elastic behaviour. The first assumption considers the macromechanical/micromechanical behaviour of the lamina, while the latter implies the validation of the generalized Hooke’s law. The laminated plate in Figure 1 has a total thickness $h$ and is composed of a total of $N$ orthotropic layers and the orientation is at angle ${\theta }_k$ to the laminate coordinate $x$. Certain layers can be for the sensing purpose, i.e., piezoelectric layers. The linear constitutive equations for the $kth$ orthotropic-piezoelectric lamina can be written as [29]:

(8)$\begin{array}{c}{\left\{\begin{array}{c}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_4\\ {\sigma }_5\\ {\sigma }_6\end{array}\right\}}^{\left(k\right)}={\left[\begin{array}{ccccc}{Q}_{11}& Q_{12}& 0& 0& 0\\ Q_{12}& Q_{22}& 0& 0& 0\\ 0& 0& Q_{44}& 0& 0\\ 0& 0& 0& Q_{55}& 0\\ 0& 0& 0& 0& Q_{66}\end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{\epsilon }_1-{\alpha }_1\Delta T\\ {\epsilon }_2-{\alpha }_2\Delta T\\ {\epsilon }_4\\ {\epsilon }_5\\ {\epsilon }_6\end{array}\right\}\\ -{\left[\begin{array}{ccc}0& 0& e_{31}\\ 0& 0& e_{32}\\ 0& e_{24}& 0\\ e_{15}& 0& 0\\ 0& 0& 0\end{array}\right]}^{\left(k\right)}\left\{\begin{array}{c}{E}_1\\ E_2\\ E_3\end{array}\right\}\end{array}$

where $Q_{ij}^{\left(k\right)}$ are the elastic stiffnesses, $e_{ij}^{\left(k\right)}$ are the piezoelectric coefficients of the $kth$ lamina, and ${\sigma }_i$, ${\epsilon }_i$, and $E_i$ are the stress, strain, and electric field components, respectively. The thermal expansions are denoted by ${\alpha }_1$ and ${\alpha }_2$, and $\Delta T$ is represents the temperature increment from a reference state. It should be noted that for the layers that do not contain piezoelectric materials, the part of the piezoelectric coefficients $e_{ij}^{\left(k\right)}$ in Equation (5) is omitted. The coefficients $Q_{ij}^{\left(k\right)}$ are given in terms of engineering constants of the $kth$ layer as:

(9)$\begin{array}{c}{Q}_{11}=\frac{E_1}{1-{\upsilon }_{12}{\upsilon }_{12}}, Q_{12}=\frac{{\upsilon }_{12}{E}_2}{1-{\upsilon }_{12}{\upsilon }_{12}}, Q_{22}=\frac{E_2}{1-{\upsilon }_{12}{\upsilon }_{12}},\\ Q_{66}=G_{12}, Q_{22}=G_{23}, Q_{55}=G_{13}\end{array}$

The laminated composite is made of several orthotropic layers, and their material axes are oriented arbitrarily with respect to the laminate coordinate. Thus, the constitutive equations of each layer must be transformed to the laminate coordinate ($x$,$y$,$z$). Consequently, the stress-strain relations in the laminated coordinates are given as:

(10)$\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{yz}\\ {\sigma }_{xz}\\ {\sigma }_{xy}\end{array}\right\}=\left[\begin{array}{ccccc}{\overline{Q}}_{11}& {\overline{Q}}_{12}& 0& 0& {\overline{Q}}_{16}\\ {\overline{Q}}_{12}& {\overline{Q}}_{22}& 0& 0& {\overline{Q}}_{26}\\ 0& 0& {\overline{Q}}_{44}& {\overline{Q}}_{45}& 0\\ 0& 0& {\overline{Q}}_{45}& {\overline{Q}}_{55}& 0\\ {\overline{Q}}_{16}& {\overline{Q}}_{26}& 0& 0& {\overline{Q}}_{66}\end{array}\right](\left\{\begin{array}{c}{\epsilon }_{xx}\\ {\epsilon }_{yy}\\ {\epsilon }_{yz}\\ {\epsilon }_{xz}\\ {\gamma }_{xy}\end{array}\right\}-\left\{\begin{array}{c}{\overline{\alpha }}_1\\ {\overline{\alpha }}_2\\ 0\\ 0\\ {\overline{\alpha }}_6\end{array}\right\}\mathrm{\Delta }T)-\left[\begin{array}{c}{\overline{e}}_{31}\\ {\overline{e}}_{31}\\ 0\\ 0\\ {\overline{e}}_{36}\end{array}\right]E_z$

where ${\overline{Q}}_{ij}$, (${\overline{\alpha }}_1$,${\overline{\alpha }}_2$,${\overline{\alpha }}_3$), and ${\overline{e}}_{ij}$ are the transformed elastic stiffnesses, thermal coefficients of expansion, and piezoelectric coefficients, respectively:

(11)$\begin{array}{c}{\overline{Q}}_{11}=Q_{11}{\cos }^4\theta +2\left(Q_{12}+2Q_{66}\right){\sin }^2\theta {\cos }^2\theta +Q_{22}{\sin }^4\theta \\ {\overline{Q}}_{12}=\left(Q_{11}+Q_{22}-4Q_{66}\right){\sin }^2\theta {\cos }^2\theta +Q_{12}({\sin }^4\theta +{\cos }^4\theta )\\ {\overline{Q}}_{22}=Q_{11}{\sin }^4\theta +2\left(Q_{12}+2Q_{66}\right){\sin }^2\theta {\cos }^2\theta +Q_{22}{\cos }^4\theta \\ {\overline{Q}}_{16}=\left(Q_{11}-Q_{12}-2Q_{66}\right)\sin \theta {\cos }^3+(Q_{12}-Q_{22}+2Q_{66}) {\sin }^3\theta \cos \theta \\ {\overline{Q}}_{16}=\left(Q_{11}-Q_{12}-2Q_{66}\right)\sin \theta \cos \theta +(Q_{12}-Q_{22}+2Q_{66}) \sin \theta {\cos }^3\theta \\ {\overline{Q}}_{66}=\left(Q_{11}+Q_{22}-2Q_{12}-2Q_{66}\right){\sin }^2\theta {\cos }^2\theta +Q_{66}({\sin }^4\theta +{\cos }^4\theta )\\ {\overline{Q}}_{44}=Q_{44}{\cos }^2\theta +Q_{55}{\sin }^2\theta \\ {\overline{Q}}_{45}=\left(Q_{55}-Q_{44}\right) \cos \theta \sin \theta {\overline{Q}}_{55}=Q_{55}{\cos }^2\theta +Q_{44}{\sin }^2\theta \end{array}$

(12)$\begin{array}{c}{\overline{\alpha }}_1={\alpha }_1{\cos }^2\theta +{\alpha }_2{\sin }^2\theta \\ {\overline{\alpha }}_2={\alpha }_1{\sin }^2\theta +{\alpha }_2{\cos }^2\theta \\ {\overline{\alpha }}_1=2\left({\alpha }_1-{\alpha }_2\right)\sin \theta \cos \theta \end{array}$

(13)$\begin{array}{c}{\overline{e}}_{31}=e_{31}{\cos }^2\theta +e_{32}{\sin }^2\theta \\ {\overline{e}}_{32}=e_{31}{\sin }^2\theta +e_{32}{\cos }^2\theta \\ {\overline{e}}_{36}=\left(e_{31}-e_{32}\right)\sin \theta \cos \theta \end{array}$

2.2. Model Construction

2.2.1. Composite Laminated Plate

The non-DS system proposed in this study consisted of a laminated composite plate with a piezoelectric patch attached to its surface. The laminated composite is an assembly of layers (plies) of fibrous composite materials. The fibers are joined together to provide specific engineering properties, including in-plane stiffness, bending stiffness, and bending strength. The properties of the composite laminate can be significantly altered by changing the orientation of the fibers in each ply. In this study, a rectangular epoxy-glass laminated composite is used as the host plate. Each ply of the epoxy-glass laminate had a thickness of 1 mm. This study proposed two different laminated composite plate models: (i) a single-lamina composite oriented at 0, 45, and 90-degrees. (ii) a double-lamina composite stacked by five different stacking sequences. Table 1 summarizes the fiber orientations of the proposed laminated composite models.

The properties and dimensions of the laminated composite plate in Figure 1 are compiled in Table 2 below.

2.2.2. Piezoelectric Patches

PEHs utilize the piezoelectric effect in which an electric charge is induced in response to the applied mechanical strain (i.e., to change in the material polarization). This is referred to as the direct piezoelectric effect. When mechanical stress is applied on the piezoelectric materials, the geometry of their atomic structure changes due to the movement of the positive and negative ions resulting in polarization. Figure 2a presents an illustration of the direct piezoelectric effect.

In linear piezoelectricity, the equations of elasticity are coupled to the charge equation of electrostatics by the mean of piezoelectric electric constants (APDL^® adopts IEEE standard form) and the piezoelectric effect can be described in the strain-charge form as:

(14)$\left\{S\right\}=\left[s^E\right]\left\{T\right\}-\left[d^T\right]\left\{E\right\}$

(15)$\left\{D\right\}=\left[d\right]\left\{T\right\}+\left[{\epsilon }^T\right]\left\{E\right\}$

where, $\left\{S\right\}$ is the elastic strain vector, $\left\{T\right\}$ is the stress vector, $\left\{E\right\}$ is the electric field intensity vector and $\left\{D\right\}$ is the electric flux density vector. The matrix $\left[s^E\right]$ is the elastic compliance matrix, $\left[d\right]$ is the piezoelectric strain matrix and $\left[{\epsilon }^T\right]$ is the dielectric permittivity matrix at constant stress. For this piezoelectric model, the direction of the positive polarization coincides with the z-axis of the rectangular system (as shown in Figure 2b). The piezoelectric charge constants involved in this model are as follows:

• (i). Compression mode $d_{33}$: indicates the polarization generated in the 3-direction per unit of mechanical compression stress applied in the 3-direction to the piezoelectric. Or induced strain in the 3-direction per unit electric field applied in the 3-direction [30,31].

• (ii). Transverse Mode $d_{31}$: indicates the polarization in the 3-direction per unit stress applied in the 1-direction. Or, induced strain in the 1-direction per unit electric field applied in the 3-direction [32]

• (iii). Shear Mode $d_{15}$: indicates the polarization developed in the 1-direction per unit shear stress 5 applied (shear around the 2-direction) when there are no other external stresses.

This study proposed six different piezoelectric patch shapes (varying by the number of edges). The thickness of the piezoelectric patches is $0.5 \mathrm{mm},$ and the material used in this study is the ceramic-type PZT-5A. The location in which the piezoelectric patches is positioned on the plate was selected to ensure optimum deformation of the piezoelectric patches. It was essential to guarantee that the piezoelectric materials could capture the optimum number of bending modes resulting from the host plate’s deformation. For this reason, the patches were positioned at $3/4$of both the plate’s length and width. Figure 3 shows the shapes and locations of the piezoelectric patches. The material properties and dimensions of the piezoelectric models are compiled in Table 3 below.

2.2.3. FE Model Components and Randomization

The FE model consisted of a simply-supported laminated composite (epoxy-glass) with a ceramic-type piezoelectric patch holding the material properties of PZT-5A. Two different models were proposed throughout this study: a single-lamina and a double-lamina composite. The fiber orientation of the single-lamina composite was varied by $0^{°}$, ${45}^{°}$, and ${90}^{°}$. The double-lamina composite had two layers with fiber orientations of: $0^{°}/0^{°}$, $0^{°}/{45}^{°}$, $0^{°}/{90}^{°}$, ${45}^{°}/{90}^{°}$, and ${45}^{°}/{45}^{°}$. Six different piezoelectric patch shapes were proposed, differing in the number of edges. Figure 4 presents the benchmark model adopted in this study.

2.2.4. Randomization

To promote a dynamic randomization condition to the structure, 30 point masses were randomly attached to the top surface of the plate. The plate is then excited by 20 randomly located time-harmonic point forces, as shown in Figure 4. This randomized condition generated structural uncertainties which in turn produced the effect of an infinite plate. Hence, the real-life response found in non-DS was approached; where the structural response at the high-frequency became sensitive to such uncertainties. Ultimately, the ensemble responses were obtained, and the average response was used for the investigation.

2.2.5. Elements and Mesh

ANSYS was used for the meshing and simulation, as shown in Figure 5. The FE model was meshed using four different element types: SOLID185 was used for the host plate, SOLID5 for the piezoelectric model, MASS21 for the point mass, and CIRCU94 for the resistive load. SOLID185 is used for 3-D modelling of solid structures. It is defined by eight nodes having three degrees of freedom (DOF) at each node (translations in the nodal x, y, and z directions). SOLID5 is a 3-D couple-field solid element and is used for piezoelectric and structural couple-field analysis. This element has eight nodes with up to six DOFs at each node (three translations, one temperature, one voltage, and one magnetic potential). MASS21 is a point element having up to six DOF (three translations and three rotations) and is defined by a single node. As for the CIRCU94 element, it is a circuit element for use in piezoelectric-circuit analyses. This element can interface with SOLID5 to form a harvesting circuit. Information about the selected elements and their capabilities can be found in ANSYS^® theory reference [33]

2.2.6. Boundary Condition

The model in Figure 5 has two different boundary conditions: (1) mechanical boundary on the host plate, (2) electrical boundary on the piezoelectric model. The 3-D element used to model the plate has only three translations DOF, and no rotations, the modelling of a simply supported boundary condition can be complicated. Therefore, when considering the plate’s element, the translations DOF of the plate’s lower edges were constrained to approach the effect of a simply-supported boundary condition. As for the electrical boundary condition, the nodes on the upper surface of the piezoelectric patch were coupled to have the same electrical potential throughout the surface; consequently, the surface nodes would act as one. Similarly, the nodes on the bottom surface (electrical ground) were coupled and constrained to zero voltage. Ultimately, the upper electrode of the patch was connected to the bottom electrode (ground) through a very high resistive load. Hence, the electrical potential could be measured at the top piezoelectric surface.

3. Results and Analysis

This section presents the results of the ensemble-average responses of the different PEH shapes attached with the non-DS composite laminate. Two different analyses were performed to analyze the response of the harvesting system: (1) modal analysis, and (2) harmonic analysis. The non-DS plate is typically subjected to uncertainties; therefore, one must ensure that a proper frequency range is selected whereby the non-deterministic behaviour can be observed, which is typically noticed at the high-frequency range. A 5% loss factor $\left(\eta \right)$ was selected; this corresponds to a damping ratio $\left(\zeta \right)$ of 2.5% since $\eta =2\zeta$. The harmonic analysis was set to a frequency range from 0–800 $\mathrm{Hz}$ (800 substeps) with the excitation of the 20 time-harmonic randomly located point forces each of 5N. For the harmonic analysis, 10 runs were performed for each piezoelectric patch shape. The location of the point masses point forces was changed in each run. The 10 responses were then averaged to obtain the ensemble average response. The modal analysis was performed prior to the harmonic analysis to analyze the laminated composite’s natural frequencies and mode shapes. Such analysis was crucial to ensure that the PEH are experiencing the maximum strain at the specified location on the plate; thus, effectively generating maximum power capacity. Additionally, the modal analysis was performed to identify the number of bending modes present within the specified frequency range. Table 4 compiles the natural frequencies and the mode shapes of the composite plate obtained from the modal analysis by the FE method. Full harmonic analysis was performed due to its superior efficiency at high frequencies.

Cumulative frequency analysis was performed to get more insight into the performance of the six patch shapes within different frequency ranges. The cumulative frequency average voltage ($\mathrm{CFAV}$) clearly represents the modes as the value peaks at the region where the resonant modes occur. As the value drops, it means that not many modes are occurring within the respective frequency range. $\mathrm{CFAV}$ was obtained by calculating the integral of the voltage response $\left(V\right)$ over the studied frequency range $\left(f\right)$ as:

(16)${\mathrm{CFAV}}_{\mathrm{i}}=\frac{1}{{\mathrm{f}}_{\mathrm{i}}}\underset{0}{\overset{{\mathrm{f}}_{\mathrm{i}}}{{\displaystyle \int }}}\mathrm{V}\left(\mathrm{f}\right) \mathrm{df} ,{ \mathrm{where} \mathrm{f}}_{\mathrm{i}}=1,2,3\dots 800$

3.1. Single-Lamina Composite

3.1.1. Ensemble-Average Voltage Response of Single-Lamina

Figure 6 compiles the ensemble electric potential responses for the 10 randomized structures (for each of the seven piezoelectric patch shapes). The response of each shape was investigated at the three different fiber orientation of the laminated composite: $0^{°}$,${45}^{°}$ and ${90}^{°}$. The responses of the PEH hold the characteristics of the typical response of the non-DS structure. The peaks were clearly seen occurring the natural frequencies of the structure. In the low-frequency, the peaks were distinct, well-spaces, and high response variance could be observed. On the other hand, the peaks were no longer distinct at the high-frequency range, and the response demonstrated low response variance. At high-frequencies, the length of the bending wavelength deformation becomes comparable to the distances between the random point masses, and thus, the structure’s response became more sensitive to structural uncertainties, which in turn increases the variance. Compared to previous studies, it is important to note that the piezoelectric voltage did not exceed the breakdown voltage at which the piezoelectric property is lost [34].

The ensemble-average responses are compiled and plotted together for each fiber orientation in Figure 7 below. It could be clearly observed that the ensemble-average response was smoothed out as the system vibrates at the high-frequency, which justifies the purpose of using the SEA method to predict/analyze the response of the non-DS at higher-frequencies.

3.1.2. CFAV of Single-Lamina for Optimal Fiber Orientation

This analysis was performed to identify the optimal orientation of fibers used within the laminated composite, as shown in Figure 8. The CFAV was calculated over the entire frequency range for each fiber orientation for each piezoelectric patch shape. It was deduced that the performance differs from one orientation to the other; with the $0^{°}$ and ${90}^{°}$fiber orientations having a higher performance than that of the ${45}^{°}$. However, it could be clearly noticed that the harvesting performance at both the $0^{°}$ and ${90}^{°}$fiber orientations are homogeneous (for every piezoelectric patch shape).

3.1.3. CFAV of Single-Lamina for Optimal Piezoelectric Shape

This section calculated the CFAV over the entire frequency range for each of the seven piezoelectric shapes. This analysis was performed to compare the performance of the piezoelectric patch shapes at different fiber orientations. It was deduced that the performance varies at different frequency ranges. Figure 9 below depicts the harvesting performance of the PEHs at different frequency ranges. It was noticed that the piezoelectric patches with a lower number of edges (i.e., triangular, quadrilateral and hexagonal) patches contribute to a higher power generation than those with a higher number of edges. At frequencies below 25 $\mathrm{Hz}$, the triangular patch possessed the highest power generation capacity, followed by the quadrilateral patch. A similar outcome could be realized at different fiber orientations within the low-frequency range. At higher-frequency ranges, it was noticed that the harvesting performance of PEHs with a higher number of edges (i.e., dodecagonal and octagonal) started to increase. The dodecagonal patch contributed to the highest overall $\mathrm{CFAV}$.

3.2. Double-Layer Lamina

3.2.1. Ensemble-Average Voltage Response of Double-Lamina

The ensemble-responses of the randomized structures (PEHs attached to the non-DS laminated composite) are presented in Figure 10. The response of each piezoelectric patch shapes was obtained at each stacking sequence with fiber orientation of: $0^{°}/0^{°}$,$0^{°}/{45}^{°}$,$0^{°}/{90}^{°}$,${45}^{°}/{90}^{°}$, and ${45}^{°}/{45}^{°}$. The responses exhibited similar characteristics to those found in Figure 6. At the low-frequency, the peaks distinct and contribute to higher amplitudes; whereas, the peaks became broader and indistinct at the high-frequency. In addition, the difference in magnitude between responses at the higher frequency range was not significant. The ensemble-average responses are compiled and plotted together for each fiber orientation in Figure 11 below.

3.2.2. CFAV of Double-Lamina for Optimal FIBER orientation

In Figure 12, the $\mathrm{CFAV}$ was calculated over the entire frequency range to have an insight into the optimal orientation of fibers in the case of double-lamina. It was noticed that the performance of the PEHs when the fibers are orientated at $0^{°}/0^{°}$ and $0^{°}/{90}^{°}$ hold the highest harvesting capacities as compared to the other orientations, followed by the fiber orientations of $0^{°}/{45}^{°}$ and ${45}^{°}/{90}^{°}$. It was found that the performance in the case of the ${45}^{°}/{45}^{°}$ orientations were the lowest. The performance of the ceramic-type piezoelectric harvesters can be depicted as the patch is operating in the compression $d_{33}$ or the transverse mode $d_{31}$. While the $d_{15}$ is mostly not efficiently used in energy harvesting because its relation is based on shear stress [35]. In the compression mode, the charge constant $d_{33}$ is identified, where the the polarization generated in the 3-direction per unit of mechanical compression stress applied in the 3-direction to the piezoelectric (as shown in Figure 2b). As for the transverse mode, the constant $d_{31}$ indicates the polarization in the 3-direction per unit stress applied in the 1-direction. Both cases are depicted for $0^{°}$ and ${90}^{°}$ fibers orientations where the performance of the PEH is optimal. Hence, it is expected that the performance at ${45}^{°}/{45}^{°}$ orientation is lower.

3.2.3. CFAV of Double-Lamina for Optimal Piezoelectric Shape

In this section, the total $\mathrm{CFAV}$ was calculated for each piezoelectric patch shape and compared at each fiber orientation. It could be noticed that the performance of the PEH is altered when the structure is vibrating in the low and at the high-frequency range. At almost every fiber orientation, it could be deduced the quadrilateral patch possessed the highest harvesting performance as the structure was vibrating at frequencies below 25 $\mathrm{Hz},$ as shown in Figure 13. Patches with a higher number of edges contributed to lower energy harvesting at the low-frequency range. As the structure was vibrating at higher frequency, the performance of the patches with higher number of edges (i.e., octagonal, and dodecagonal) surpassed the quadrilateral patch in the total $\mathrm{CFAV}$. However, an increase in the performance of the triangular patch was detected at the high-frequency range. Overall, it was observed that the octagonal patch contributed to the highest overall $\mathrm{CFAV}$.

Figure 14 represents the plots for the normalized maximum $\mathrm{CFAV}$ for the single and double-lamina composites, where the variation in performance at different frequency ranges can be clearly depicted. For a single layer composite, it is observed that the performance of the triangular patch is higher at lower frequency bands, followed by the quadrilateral patch with almost 30% drop in the performance. The pentagonal patch had the lowest harvesting performance with 75% drop for frequencies below 25 Hz. The dodecagonal patch had the highest performance at intermediate and broadband frequencies, followed by the octagonal patch with around 20% lower harvesting performance.

4. Conclusions

In this study, an investigation of the energy harvesting from seven piezoelectric patch shapes (differing in the number of edges) was conducted. PEH models were attached to a non-DS laminated composite plate (single-lamina and double-lamina), and the performance was examined through a coupled-field FE model. The structures’ dynamics were randomized by subjecting the host plate to randomly distributed point masses and randomly located harmonic point forces. The randomized condition generated structural uncertainties required to achieve the effect of an infinite plate structure with the response of a non-DS system. The analysis was performed by employing the ensemble-average responses of the developed electrical voltage across the piezoelectric models. In the case of the single-lamina composite, it was observed that the performance was optimal when the fibers were oriented at either $0^{°}$ or ${90}^{°}$. At frequencies below 25$\mathrm{Hz}$, the quadrilateral and triangular patches possessed higher harvesting performance. The latter had a higher power generation capacity; however, the patches had a higher number of edges (i.e., dodecagonal, octagonal and hexagonal) exhibited better performance at higher frequency ranges. The dodecagonal patch contributed to the highest overall $\mathrm{CFAV}$ in the case of single-lamina composite. As for the double-lamina composite, the performance was found optimal at $0^{°}/0^{°}$ and $0^{°}/{90}^{°}$ fiber orientations. The performance of the quadrilateral patch was highest within the low-frequency range but decreased drastically at higher frequency ranges. It was observed that the octagonal patch contributed to the highest overall $\mathrm{CFAV}$ for the double-lamina composite. Figure 13 represents the normalized maximum cumulative voltages for single and double-lamina composites.

Author Contributions

Conceptualization, A.G.A.M.; methodology, A.G.A.M. and A.A.; software, A.A. and J.R.; validation, A.G.A.M., J.R. and A.N.W.; formal analysis, A.A., K.A.M.N. and N.H.D.N.; investigation, A.G.A.M. and A.A.; resources, A.G.A.M. and A.N.W.; data curation, A.A., K.A.M.N. and N.H.D.N.; writing—original draft preparation, A.A.; review and editing, A.G.A.M., J.R.; visualization, A.G.A.M., J.R. and A.N.W.; supervision, A.G.A.M.; project administration, A.G.A.M.; funding acquisition, A.G.A.M.; All authors have read and agreed to the published version of the manuscript.

Funding

The authors gratefully acknowledge the financial support by Qatar University-International Research Collaboration under grant number IRCC-2020-017, and the Qatar National Research Fund through the National Priorities Research Program under grant number NPRP 11S-1220-170112.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on reasonable request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest in preparing this article.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Coordinates of laminated composite plate.

Figure 2. Piezoelectric model. (a) Direct piezoelectric effect; (b) Directions of forces affecting a piezoelectric element.

Figure 3. Shapes and location of the piezoelectric patches on the plate.

Figure 4. Benchmark model.

Figure 5. FE Model mesh and boundary condition.

Figure 6. Ensemble-average voltage response for single-lamina.

Figure 6. Ensemble-average voltage response for single-lamina.

Figure 7. Compiled average response for single-lamina.

Figure 8. CFAV of the piezoelectric patches for optimal fiber orientation for single-lamina.

Figure 9. (a) frequencies below 25 Hz range (b) broadband frequency ranges for single-lamina.

Figure 10. Ensemble-average voltage response for double-lamina.

Figure 10. Ensemble-average voltage response for double-lamina.

Figure 11. Compiled average-response for double-lamina.

Figure 12. CFAV of the piezoelectric patches for optimal fiber orientation for double-lamina.

Figure 13. Performance (a) frequencies below 25 Hz range (b) broadband frequency ranges for double-lamina.

Figure 14. Normalized maximum cumulative voltage (a) single-lamina and (b) double-lamina composites.

Table 1

Proposed laminated composite models.

Laminated Composite Plate
Model	Fiber Orientation
Single-lamina composite	[Image omitted. Please see PDF.]
Double-lamina composite	Model A	$$0^{°}/0^{°}$$	[Image omitted. Please see PDF.]
	Model B	$$0^{°}/{45}^{°}$$	[Image omitted. Please see PDF.]
	Model C	$$0^{°}/{90}^{°}$$	[Image omitted. Please see PDF.]
	Model D	$${45}^{°}/{90}^{°}$$	[Image omitted. Please see PDF.]
	Model E	$${45}^{°}/{45}^{°}$$	[Image omitted. Please see PDF.]

Table 2

Material properties and dimensions of the laminated composite.

Length		$$\mathbf{100} cm$$
Width		$$\mathbf{120} cm$$
Ply Thickness		$$\mathbf{0.1} cm$$
Material		Epoxy-Glass
Modulus of elasticity	$$E_x$$	$$54 \mathrm{GPa}$$
	$$E_y$$	$$54 \mathrm{GPa}$$
	$$E_z$$	$$4.8 \mathrm{GPa}$$
Shear modulus	$$G_{xy}$$	$$3.16 \mathrm{GPa}$$
	$$G_{yz}$$	$$1.78 \mathrm{GPa}$$
	$$G_{xz}$$	$$1.78 \mathrm{GPa}$$
Poisson’s ration	$${\upsilon }_{xy}$$	$$0.06$$
	$${\upsilon }_{yz}$$	$$0.313$$
	$${\upsilon }_{xz}$$	$$0.313$$
Density		$$2000$$

Table 3

Material properties and dimensions of the piezoelectric models.

Material Property	Coefficient/Shape	Value
Compliance Matrix($${10}^{-12} {\mathrm{m}}^2$$/N)	$$s_{11}=s_{22}$$ $$s_{12}$$ $$s_{13}=s_{23}$$ $$s_{33}$$ $$s_{44}=s_{55}$$ $$s_{66}$$	$$16.4$$ $$-5.74$$ $$-7.22$$ $$18.8$$ $$47.5$$ $$44.3$$
Relative Permittivity Matrix($${\mathsf{\epsilon }}_{\mathrm{o}}=8.854·{10}^{-12}$$F/m)	$${\epsilon }_{11}={\epsilon }_{22}$$ $${\epsilon }_{33}$$	$$1730$$ $$1700$$
Piezoelectric Strain Matrix ($${10}^{-12}$$ C/N)	$$d_{15}=d_{24}$$ $$d_{31}=d_{32}$$ $$d_{33}$$	$$584$$ $$-171$$ $$374$$
Density (kg/m3)	$$\rho$$	$$7750$$
Side Length ($$\mathrm{cm}$$)	Triangular	$$6.8$$
	Quadrilateral	$$4.5$$
	Pentagonal	$$3.4$$
	Hexagonal	$$2.8$$
	Octagonal	$$2$$
	Dodecagonal	$$1.3$$
	Circular	$$2.8 \mathrm{radius}$$

Table 4

Natural frequencies (Hz) and mode shapes of the composite plate.

	Orientation	Mode 1	Mode 2	Mode 3	Mode 4
Single	$$0^{°}$$	7.8076	16.281	21.992	26.242
	$${45}^{°}$$	8.2785	17.486	22.536	29.368
	$${90}^{°}$$	7.8076	16.281	21.922	26.242
Double	$$0^{°}/0^{°}$$	12.265	23.161	30.253	37.537
	$$0^{°}/{45}^{°}$$	13.764	25.899	31.451	43.243
	$$0^{°}/{90}^{°}$$	12.265	23.161	30.253	37.537
	$${45}^{°}/{90}^{°}$$	12.019	23.711	30.196	40.927
	$${45}^{°}/{45}^{°}$$	13.263	26.130	31.791	45.171
Mode Shape		[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]
		[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]	[Image omitted. Please see PDF.]

Appendix A

The governing equation of motion of the FSDT is obtained using the dynamic version of the principle of virtual displacements:(A1)$$\frac{\partial N_{xx}}{\partial x}+\frac{\partial Nxy}{\partial y}=I_0\frac{{\partial }^2{u}_0}{\partial t^2}+I_1\frac{{\partial }^2{\varphi }_x}{\partial t^2}$$ (A2)$$\frac{\partial N_{xy}}{\partial x}+\frac{\partial Nyy}{\partial y}=I_0\frac{{\partial }^2{v}_0}{\partial t^2}+I_1\frac{{\partial }^2{\varphi }_y}{\partial t^2}$$ (A3)$$\frac{\partial Q_x}{\partial x}+\frac{\partial Q_y}{\partial y}+N+q=I_0\frac{{\partial }^2{w}_0}{\partial t^2}$$ (A4)$$\frac{\partial M_{xx}}{\partial x}+\frac{\partial Mxy}{\partial y}-Q_x=I_2\frac{{\partial }^2{\varphi }_x}{\partial t^2}+I_1\frac{{\partial }^2{u}_0}{\partial t^2}$$ (A5)$$\frac{\partial M_{xy}}{\partial x}+\frac{\partial Myy}{\partial y}-Q_y=I_2\frac{{\partial }^2{\varphi }_y}{\partial t^2}+I_1\frac{{\partial }^2{v}_0}{\partial t^2}$$

The quantities ($$N_{xx}$$,$$N_{yy}$$,$$N_{xy}$$) are the in-plane force resultants, ($$M_{xx}$$,$$M_{yy}$$,$$M_{xy}$$) are the moment resultants, and ($$I_0$$,$$I_1$$,$$I_2$$) are the mass moment of inertia:(A6)$$\begin{array}{c}\left\{\begin{array}{c}{N}_{xx}\\ N_{yy}\\ N_{xy}\end{array}\right\}={\int }_{-h/2}^{h/2}\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}dz,\hspace{1em}\left\{\begin{array}{c}{M}_{xx}\\ M_{yy}\\ M_{xy}\end{array}\right\}={\int }_{-h/2}^{h/2}\left\{\begin{array}{c}{\sigma }_{xx}\\ {\sigma }_{yy}\\ {\sigma }_{xy}\end{array}\right\}zdz,\hspace{1em}\left\{\begin{array}{c}{I}_0\\ I_1\\ I_2\end{array}\right\}=\\ {\int }_{-h/2}^{h/2}\left\{\begin{array}{c}1\\ z\\ z^2\end{array}\right\}{\rho }_{0}dz\end{array}$$ (A7)$$=\frac{\partial }{\partial x}\left(N_{xx}\frac{\partial w_0}{\partial x}+N_{xy}\frac{\partial w_0}{\partial y}\right)+\frac{\partial }{\partial y}\left(N_{xy}\frac{\partial w_0}{\partial x}+N_{yy}\frac{\partial w_0}{\partial y}\right)$$

The transverse force resultants are denoted by ($$Q_x$$,$$Q_y$$) and are given by:(A8)$$\left\{\begin{array}{c}{Q}_x\\ Q_y\end{array}\right\}=K{{\displaystyle \int }}_{-h/2}^{h/2}\left\{\begin{array}{c}{\sigma }_{xz}\\ {\sigma }_{yz}\end{array}\right\}dz=K\left[\begin{array}{cc}{A}_{44}& A_{45}\\ A_{45}& A_{55}\end{array}\right]\left\{\begin{array}{c}\begin{array}{c}{\gamma }_{yz}^{\left(0\right)}\\ {\gamma }_{xz}^{\left(0\right)}\end{array}\end{array}\right\}$$ where $$A_{44}$$, $$A_{45}$$, and $$A_{55}$$ are the transverse shear stiffnesses and $$K$$ denotes the shear correction factor; which is used to account for the discrepancy between the actual shear forces and the ones due to the contact stress state predicted by the first-order theory.

The elasticity matrix specifies the stiffness $$\left[{\mathit{c}}_{\mathit{i}\mathit{j}}\right]$$ or the compliance coefficients $$\left[{\mathit{s}}_{\mathit{i}\mathit{j}}\right]$$. The form used for the elasticity matrix utilizes the compliance form and is input as:(A9)$$\left[\mathit{s}\right]=\left[\begin{array}{cccccc}{\mathit{s}}_{\mathbf{11}}& {\mathit{s}}_{\mathbf{12}}& {\mathit{s}}_{\mathbf{13}}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ {\mathit{s}}_{\mathbf{21}}& {\mathit{s}}_{\mathbf{22}}& {\mathit{s}}_{\mathbf{33}}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ {\mathit{s}}_{\mathbf{31}}& {\mathit{s}}_{\mathbf{32}}& {\mathit{s}}_{\mathbf{33}}& \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{s}}_{\mathbf{44}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{s}}_{\mathbf{55}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{s}}_{\mathbf{66}}\end{array}\right]$$

The piezoelectric matrix can be input as either a stress matrix $$\left[\mathit{e}\right]$$ or a strain matrix $$\left[\mathit{d}\right]$$. The stress matrix is associated with the input of the anisotropic elasticity in the stiffness form $$\left[\mathit{c}\right]$$, while the strain matrix corresponds to the input of the elasticity in the compliance form $$\left[\mathit{s}\right]$$. The piezoelectric matrix is input in the strain from $$\left[\mathit{d}\right]$$ as:(A10)$$\left[\mathit{d}\right]=\left[\begin{array}{cccccc}\mathbf{0}& \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{d}}_{\mathbf{15}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}& {\mathit{d}}_{\mathbf{24}}& \mathbf{0}& \mathbf{0}\\ {\mathit{d}}_{\mathbf{31}}& {\mathit{d}}_{\mathbf{32}}& {\mathit{d}}_{\mathbf{33}}& \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right]$$

Ansys can convert the piezoelectric strain matrix $$\left[\mathit{d}\right]$$ to a piezoelectric stress matrix $$\left[\mathit{e}\right]$$ using the elasticity matrix $$\left[\mathit{c}\right]$$ by:(A11)$$\left[\mathit{e}\right]=\left[\mathit{c}\right]\left[\mathit{d}\right]$$

As for the dielectric matrix, it uses the electric permittivity and is input as a dielectric permittivity matrix at constant stress $$\left[{\mathbf{\epsilon }}^{\mathit{T}}\right]$$:(A12)$$\left[{\mathbf{\epsilon }}^{\mathit{T}}\right]=\left[\begin{array}{ccc}{\mathbf{\epsilon }}_{\mathbf{11}}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& {\mathbf{\epsilon }}_{\mathbf{22}}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& {\mathbf{\epsilon }}_{\mathbf{33}}\end{array}\right]$$

ANSYS can transform the dielectric matrix at constant stress $$\left[{\mathbf{\epsilon }}^{\mathit{T}}\right]$$ to a dielectric matrix at constant strain $$\left[{\mathbf{\epsilon }}^{\mathit{S}}\right]$$ by:(A13)$$\left[{\mathbf{\epsilon }}^{\mathit{S}}]=[{\mathbf{\epsilon }}^{\mathit{T}}\right]-{\left[\mathit{e}\right]}^{\mathit{T}}\left[\mathit{d}\right]$$

References

1. Wei, C.; Jing, X. A comprehensive review on vibration energy harvesting: Modelling and realization. Renew. Sustain. Energy Rev.; 2017; 74, pp. 1-18. [DOI: https://dx.doi.org/10.1016/j.rser.2017.01.073]

2. Liu, H.; Fu, H.; Sun, L.; Lee, C.; Yeatman, E.M. Hybrid energy harvesting technology: From materials, structural design, system integration to applications. Renew. Sustain. Energy Rev.; 2021; 137, 110473. [DOI: https://dx.doi.org/10.1016/j.rser.2020.110473]

3. Zou, H.-X.; Zhao, L.-C.; Gao, Q.-H.; Zuo, L.; Liu, M.; Tan, T.; Wei, K.-X.; Zhang, W.-M. Mechanical modulations for enhancing energy harvesting: Principles, methods and applications. Appl. Energy; 2019; 255, 113871. [DOI: https://dx.doi.org/10.1016/j.apenergy.2019.113871]

4. Wang, X. Piezoelectric nanogenerators—Harvesting ambient mechanical energy at the nanometer scale. Nano Energy; 2012; 1, pp. 13-24. [DOI: https://dx.doi.org/10.1016/j.nanoen.2011.09.001]

5. Chen, J.; Qiu, Q.; Han, Y.; Lau, D. Piezoelectric materials for sustainable building structures: Fundamentals and applications. Renew. Sustain. Energy Rev.; 2019; 101, pp. 14-25. [DOI: https://dx.doi.org/10.1016/j.rser.2018.09.038]

6. Ali, F.; Raza, W.; Li, X.; Gul, H.; Kim, K.-H. Piezoelectric energy harvesters for biomedical applications. Nano Energy; 2019; 57, pp. 879-902. [DOI: https://dx.doi.org/10.1016/j.nanoen.2019.01.012]

7. Roundy, S.; Wright, P.K. A piezoelectric vibration based generator for wireless electronics. Smart Mater. Struct.; 2004; 13, pp. 1131-1142. [DOI: https://dx.doi.org/10.1088/0964-1726/13/5/018]

8. Casciati, S.; Faravelli, L.; Chen, Z. Energy harvesting and power management of wireless sensors for structural control applications in civil engineering. Smart Struct. Syst.; 2012; 10, pp. 299-312. [DOI: https://dx.doi.org/10.12989/sss.2012.10.3.299]

9. Bodaghi, M.; Damanpack, A.; Aghdam, M.; Shakeri, M. Non-linear active control of FG beams in thermal environments subjected to blast loads with integrated FGP sensor/actuator layers. Compos. Struct.; 2012; 94, pp. 3612-3623. [DOI: https://dx.doi.org/10.1016/j.compstruct.2012.06.001]

10. Sezer, N.; Koç, M. A comprehensive review on the state-of-the-art of piezoelectric energy harvesting. Nano Energy; 2021; 80, 105567. [DOI: https://dx.doi.org/10.1016/j.nanoen.2020.105567]

11. Truitt, A.; Mahmoodi, S.N. A review on active wind energy harvesting designs. Int. J. Precis. Eng. Manuf.; 2013; 14, pp. 1667-1675. [DOI: https://dx.doi.org/10.1007/s12541-013-0226-4]

12. Park, J.; Lee, S.; Kwak, B.M. Design optimization of piezoelectric energy harvester subject to tip excitation. J. Mech. Sci. Technol.; 2012; 26, pp. 137-143. [DOI: https://dx.doi.org/10.1007/s12206-011-0910-1]

13. Wein, F.; Kaltenbacher, M.; Stingl, M. Topology optimization of a cantilevered piezoelectric energy harvester using stress norm constraints. Struct. Multidiscip. Optim.; 2013; 48, pp. 173-185. [DOI: https://dx.doi.org/10.1007/s00158-013-0889-6]

14. Liang, Z.; Xu, C.; Ren, B.; Di, W.; Li, L.; Luo, H.; Chen, Z.; Su, J. Optimization of cantilevered piezoelectric energy harvester with a fixed resonance frequency. Sci. China Ser. E Technol. Sci.; 2014; 57, pp. 1093-1100. [DOI: https://dx.doi.org/10.1007/s11431-014-5556-7]

15. Kim, C.; Lee, J. Topology optimum design of unimorph piezoelectric cantilevered Mindlin plates as a vibrating electric harvester. J. Mech. Sci. Technol.; 2014; 28, pp. 4131-4138. [DOI: https://dx.doi.org/10.1007/s12206-014-0925-5]

16. Sordo, G.; Serra, E.; Schmid, U.; Iannacci, J. Optimization method for designing multimodal piezoelectric MEMS energy harvesters. Microsyst. Technol.; 2016; 22, pp. 1811-1820. [DOI: https://dx.doi.org/10.1007/s00542-016-2848-9]

17. Song, J.; Zhao, G.; Li, B.; Wang, J. Design optimization of PVDF-based piezoelectric energy harvesters. Heliyon; 2017; 3, e00377. [DOI: https://dx.doi.org/10.1016/j.heliyon.2017.e00377]

18. Izadgoshasb, I.; Lim, Y.Y.; Lake, N.; Tang, L.; Padilla, R.V.; Kashiwao, T. Optimizing orientation of piezoelectric cantilever beam for harvesting energy from human walking. Energy Convers. Manag.; 2018; 161, pp. 66-73. [DOI: https://dx.doi.org/10.1016/j.enconman.2018.01.076]

19. Kim, T.; Ko, Y.; Yoo, C.; Choi, B.; Han, S.; Kim, N. Design optimisation of wide-band piezoelectric energy harvesters for self-powered devices. Energy Convers. Manag.; 2020; 225, 113443. [DOI: https://dx.doi.org/10.1016/j.enconman.2020.113443]

20. Ichige, R.; Kuriyama, N.; Umino, Y.; Tsukamoto, T.; Suzuki, T. Size optimization of metamaterial structure for elastic layer of a piezoelectric vibration energy harvester. Sens. Actuators A Phys.; 2021; 318, 112488. [DOI: https://dx.doi.org/10.1016/j.sna.2020.112488]

21. Peralta, P.; Ruiz, R.; Natarajan, S.; Atroshchenko, E. Parametric study and shape optimization of piezoelectric energy harvesters by isogeometric analysis and kriging metamodeling. J. Sound Vib.; 2020; 484, 115521. [DOI: https://dx.doi.org/10.1016/j.jsv.2020.115521]

22. Muthalif, A.G.A.; Wahid, A.N. Estimation and measurement of effective line mobility on a non-deterministic thin plate excited by a piezoelectric patch. J. Vibroengineering; 2020; 22, pp. 98-110. [DOI: https://dx.doi.org/10.21595/jve.2019.20865]

23. Muthalif, A.G.; Wahid, A.N.; Nor, K.A. Estimating ensemble average power delivered by a piezoelectric patch actuator to a non-deterministic subsystem. J. Sound Vib.; 2014; 333, pp. 1149-1162. [DOI: https://dx.doi.org/10.1016/j.jsv.2013.10.012]

24. Muthalif, A.G.A.; Langley, R.S. Active control of high-frequency vibration: Optimisation using the hybrid modelling method. J. Sound Vib.; 2012; 331, pp. 2969-2983. [DOI: https://dx.doi.org/10.1016/j.jsv.2012.02.012]

25. Cotoni, V.; Shorter, P.; Langley, R. Numerical and experimental validation of a hybrid finite element-statistical energy analysis method. J. Acoust. Soc. Am.; 2007; 122, pp. 259-270. [DOI: https://dx.doi.org/10.1121/1.2739420] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17614486]

26. Shastry, S.N.; Swamy, S.S. Study on high frequency vibrations using statistical energy analysis for different composite materials. IOP Conference Series: Materials Science and Engineering; IOP Publishing: Bristol, UK, 2020; 925, 012044.

27. Seçgin, A. Numerical determination of statistical energy analysis parameters of directly coupled composite plates using a modal-based approach. J. Sound Vib.; 2013; 332, pp. 361-377. [DOI: https://dx.doi.org/10.1016/j.jsv.2012.08.020]

28. Ngo-Cong, D.; Mai-Duy, N.; Karunasena, W.; Tran-Cong, T. Free vibration analysis of laminated composite plates based on FSDT using one-dimensional IRBFN method. Comput. Struct.; 2011; 89, pp. 1-13. [DOI: https://dx.doi.org/10.1016/j.compstruc.2010.07.012]

29. Soares, C.A.M.; Soares, C.M.M.; Freitas, M.J.M. Mechanics of Composite Materials and Structures; Springer: Dordrecht, The Netherlands, 1999; [DOI: https://dx.doi.org/10.1007/978-94-011-4489-6]

30. Benjeddou, A. Field-dependent nonlinear piezoelectricity: A focused review. Int. J. Smart Nano Mater.; 2018; 9, pp. 68-84. [DOI: https://dx.doi.org/10.1080/19475411.2018.1439850]

31. Fakhzan, M.A.; Muthalif, A.G.; Lumpur, K. Vibration based energy harvesting using piezoelectric material. Proceedings of the 2011 4th International Conference on Mechatronics: Integrated Engineering for Industrial and Societal Development, ICOM’11—Conference Proceedings; Kuala Lumpur, Malaysia, 17–19 May 2011; pp. 17-19.

32. Trindade, M.A.; Benjeddou, A. Finite element characterisation of multilayer d31 piezoelectric macro-fibre composites. Compos. Struct.; 2016; 151, pp. 47-57. [DOI: https://dx.doi.org/10.1016/j.compstruct.2015.10.011]

33. ANSYS® Academic Mechanical, Release 17.0, Help System, Element Reference, Ansys, Inc. Canonsburg, PA, USA, 2020.

34. Bodaghi, M.; Damanpack, A.; Aghdam, M.; Shakeri, M. Geometrically non-linear transient thermo-elastic response of FG beams integrated with a pair of FG piezoelectric sensors. Compos. Struct.; 2014; 107, pp. 48-59. [DOI: https://dx.doi.org/10.1016/j.compstruct.2013.07.045]

35. Chilabi, H.J.; Salleh, H.; Supeni, E.E.; As’Arry, A.B.; Rezali, K.A.M.; Atrah, A.B. Harvesting energy from planetary gear using piezoelectric material. Energies; 2020; 13, 223. [DOI: https://dx.doi.org/10.3390/en13010223]