Binary-Like Topology Optimization of Piezoelectric Metamaterial Plate with Interface Circuits Using Extended Plane Wave Expansion Method

Abstract

Piezoelectric metamaterial plate (PMP) is being investigated for structural vibration energy harvesting (SVEH), in which an interface circuit is often used. Thus, it is a challenge to perform bandgap optimization of such an elastic–electro–mechanical coupling system. This paper presents a binary-like topology optimization scheme by dividing the unit cell into identical pieces, where a {0, 1} matrix is optimized to indicate material distribution. Firstly, a unified motion equation is derived for the elastic plate and the piezoelectric patch, and an electromechanical coupling model is built for a self-powered synchronized charge extraction circuit. Then, an extended plane wave expansion method is presented to model the bandgap character of the PMP with interface circuits (PMPICs), and the numerical solution of the dispersion curves is derived based on the Bloch theorem. Next, an extended genetic algorithm is applied for the topology optimization of the PMPIC. In the end, numerical and finite element simulations are performed to validate the proposed method. The results demonstrate that both the structure and the circuit can be optimized simultaneously to obtain the maximum first-order bandgap at a given central frequency. Therefore, the proposed method should provide an effective solution for the topology optimization of a PMPIC for broadband SVEH.

1. Introduction

The goal of structural health monitoring (SHM) is to improve the safety, reliability, and/or ownership costs of engineering systems by autonomously monitoring their conditions and predicting incipient damages [1,2]. Most SHM systems are wired in a wide area or in a harsh environment, resulting in both economical and practical limitations. Owing to the development of MEMS technologies, low power consumption microsensors are designed to embed into engineering structures, and then condition-related signals can be sampled and transmitted wirelessly. Thus, embedded wireless sensor network (WSN) provides effective solutions for SHM [3]. Currently, however, a major concern on embedded WSN is how to provide long-term sources of power [4]. Nowadays, the most promising solution for extending the useful life of embedded WSN is to harvest environmental energy for generating electrical energy, which is commonly called energy harvesting [5].

As one of the promising energy harvesting technologies, piezoelectric vibration energy harvesting (PVEH) has been widely studied to realize self-powered WSNs due to their large energy density and ease of implementation [4,5,6,7,8]. In most studies, piezoelectric cantilever beams with a proof mass are used as the medium. As for SHM systems, however, a cantilever-like PVEH device has two main drawbacks. The first one is that it needs extra space to hold a bulky proof mass and clamping part. The second one is that it is difficult to effectively harvest low-frequency vibrations due to the length limitation. To cope with these drawbacks, piezoelectric patch-based vibration energy harvesters have been proposed to convert the vibration energy of host structures into electrical energy [9,10,11]. Compared with cantilever-like ones, patch-like PVEH devices can be easily integrated into their host structures and are especially feasible for thin structures used in the fields of aerospace, automotive, and marine applications. The key is how to effectively control vibration propagation. During the past two decades, artificial materials and structures with periodic modulations of physical properties have attracted significant attention, especially bandgaps, which are defined as those frequency bands where waves are completely blocked [12]. In particular, bandgaps just meet the challenges of low-frequency and broadband vibration energy harvesting in SHM systems, and artificial material with bandgaps is often called bandgap material. According to the generation mechanisms of bandgaps, periodic structural configurations can be used to design phononic crystals or metamaterials for energy harvesting. Chen et al. investigated one-dimension phononic piezoelectric cantilever beams to widen the resonant bandwidth of a harvester [13]. Wu et al. first used a phononic crystal to harvest acoustic energy [14]. However, the lattice constant of phononic crystal is on the order of wavelength so that its length scale is often too large for low-frequency bandgaps. Liu et al. firstly proposed the concept of acoustic metamaterial for solving the length-scale problem of bandgap materials [15]. Due to the local resonant behaviors of these auxiliary oscillators, bandgaps can be generated in a low-frequency regime. Therefore, metamaterial-like structures provide an effective way for harvesting structural vibrations, which are also called metastructures such as metasurfaces. Sugino et al. presented analytical and numerical investigations on a one-dimensional locally resonant piezoelectric metastructure with segmented electrodes under transverse vibrations [16]. Wang et al. proposed to harvest energy from vortex-induced vibration by using metasurface [17]. Kherraz et al. studied numerically and experimentally locally resonant bandgaps in a homogeneous piezoelectric plate covered by a 1D periodic array of thin electrodes connected to inductive shunts [18]. Chen et al. first investigated one kind of piezoelectric metamaterial plate (PMP) for vibration energy harvesting and analyzed the effects of geometric and material parameters on vibration bandgaps by finite element simulations [19]. Later, Chen et al. continued to build an elastic–electro–mechanical model of the PMP with interface circuits (PMPICs) [20]. Zhang et al. proposed an acoustic metamaterial that had dual functionality of sound insulation and energy harvesting simultaneously [21]. Existing studies have shown that PMPs provide a promising way of harvesting structural vibrations.

Bandgap characteristics of periodic structures are highly dependent on the topology of the unit cell. In recent years, many studies have been conducted to design and optimize the unit cell in order to obtain targeted bandgaps. In particular, topology optimization has been identified as a powerful computational tool for the design of bandgap material [22]. Sigmund et al. further improved the above topology optimization method by considering both scalar and coupled problems [23]. A 2D phononic crystal with an arbitrarily asymmetric lattice was also obtained by combining the genetic algorithm with the finite element model [24]. Zhong et al. maximized the bandgap width of a phononic crystal by using a GA combined with the fast plane wave expansion method [25]. Xie et al. proposed an improved fast plane wave expansion method for topology optimization of phononic crystals [26]. In 1993, Xie and Steven first attempted to carry out topology optimization of continuum structures using binary {0, 1} variables [27]. In 2018, Sivapuram and Picelli proposed the topology optimization of binary structures (TOBS) method, where the material layout is defined by a set of binary {0, 1} design variables [28]. By now, the TOBS method has attracted much attention due to the advantage of suiting well separate physics domains. The current challenge for structural topology optimization is to develop reliable techniques to account for different physics interactions. In 2020, Sivapuram and Picelli further extended the TOBS method to design topologies by considering fluid-structure interaction loads [29], pressure, and thermal loads simultaneously [30].

A PMPIC is composed of a piezoelectric metamaterial plate and an interface circuit. The aim of PMPIC design is to obtain proper bandgaps to meet the needs of low frequency and broadband. Bandgaps are strongly related to elastic wave propagation and the interface circuit; therefore, elastic–electro–mechanical coupling effects should be considered in topology optimization. To the best of the author’s knowledge, this is the first work to optimize the topology of a PMPIC by considering elastic–electro–mechanical coupling behaviors. In order to optimize the topology of the PMPIC, a binary topology optimization process is first proposed based on an extended plane wave expansion (PWE) method in this paper. The remainder of this paper is organized as follows. In Section 2, the problem of topology optimization of a PMPIC is outlined, and a binary-like topology optimization scheme is presented. In Section 3, a unified motion equation is derived for the elastic plate and the piezoelectric patch, and an electromechanical coupling model is built. Then, the bandgap character of the PMPIC is modeled based on an extended plane wave expansion method and the Bloch theorem is used to solve the above model. In Section 4, an extended genetic algorithm is applied for the topology optimization of the PMPIC. Numerical and finite element simulations are carried out in Section 5. Finally, the conclusions are summarized in Section 6.

2. Problem Formulation of Topology Optimization of a PMPIC

2.1. Basic Configuration of a PMPIC

A schematic diagram of a PMPIC is shown in Figure 1. A thin rectangular plate is used as the thin plate, which can be decomposed into periodic unit cells. In each unit cell, a hole is cut off and then filled by a piezoelectric patch that supports one or multiple small square masses. The piezoelectric patch and the masses form a local resonator. The top and bottom electrodes of the piezoelectric patch are connected to an interface circuit for electric power conversion. Under excitations, transverse vibrations of the PMP will propagate along the plane of the thin plate. Due to the existence of local resonators, vibrations within specific frequency bands cannot propagate through the periodic cells, leading to locally resonant bandgaps. In this case, vibration energies will be localized in some unit cells, which then can be converted into electricity through the piezoelectric patch.

According to the bandgap theory, the band structure characteristics of the PMPIC are strongly dependent on the topology of the unit cell. Key components of each unit cell include the thin plate, the piezoelectric patch, the masses, and the interface circuit. Thus we need to optimize the topology of the unit cell to meet the band requirements of the PMPIC in practice, including material type, size, and shape, etc.

2.2. Binary Topology Optimization of the PMPIC

In this paper, the goal is to maximize the bandgap width so that more vibration energy can be harvested by the PMPIC. Generally speaking, a topology optimization problem always consists of distributing two material phases in a fixed design domain, and then an objective function is minimized or maximized under a number of constraints. The design domain is discretized into many small elements and material properties are assumed to be constant within each element. The design variable is the material density, which is equal to 0 or 1. If the design variable takes a zero value, we should have the first phase material in the element. Otherwise, we should have the second phase material in the element. In this problem, the physical size, shape, and connectivity of the unit cell are unknown in advance.

For structural topology optimization, a common method is called the solid isotropic material with penalization (SIMP) approach [31]. Material properties are modeled as the relative material density. In this case, design variables are the element relative densities. However, the drawback of this approach is that no physical material exists with properties described by power–law interpolation. In order to overcome this limitation, the levelset method has been recently introduced in the field of topology optimization [32]. Level set-based topology optimization methods can directly provide clear boundaries in optimal configurations that avoid the presence of intermediate densities. The level set method has been successfully used for topology optimization of metamaterials for the sake of controlling elastic wave propagation. However, the level set method is not feasible for optimizing the PMPIC due to the following two reasons: (i) the optimization function is difficult to be defined by including the interface circuit and the piezoelectric effect and (ii) the optimal result of the level set method may cause the piezoelectric patch to many independent pieces. In this case, the result may be inaccurate due to the piezoelectric effect, and PVEH is also inconvenient to carry out.

In this paper, a binary-like topology optimization scheme is proposed for the PMPIC, which is inspired by the idea of the SIMP approach. As shown in Figure 2, a square PMPIC unit cell is considered, where the interface circuit is not plotted. Then, the unit cell is divided into small elements and each element is marked by 0 or 1. If an element takes a zero value, we should choose thin plate material in the element. Otherwise, we should choose piezoelectric material in the element. In this way, the PMPIC unit cell can be represented by a 0–1 matrix. Then, the optimization problem is formulated to find the optimal topology of the 0–1 matrix for maximizing the bandgap width. The advantages of the above scheme include (i) compared with the SIMP approach, the above binary-like topology optimization method will not generate nonexistent materials with intermediate densities; (ii) compared with the level set approach, the above binary-like topology optimization method can easily control the connectivity of all “1” digits so that the piezoelectric patch forms a whole; (iii) the above binary-like topology optimization problem can be efficiently solved by common 0–1 programming algorithms or evolutionary optimization algorithms. In particular, it should be noted that each isolated host structure in the optimal result will be used as one mass, as shown in Figure 1, instead of mounting additional masses. That is to say, multiple masses may exist in the unit cell.

3. Bandgap Algorithm of the PMPIC by Extended Plane Wave Expansion Method

In recent years, piezoelectric phononic crystals with shunting circuits were studied to adjust and control bandgap characteristics for vibration suppression [33]. Lian et al. proposed to solve piezoelectric phononic crystal connected with resonant shunting circuits using an enhanced PWE method [34]. Inspired by the idea, the authors extended the PWE method to model and solve the PMPIC in this section.

3.1. Transverse Motion Equation of the PMPIC

Based on the Kirchhoff plate theory [35], transverse shear strains of a thin plate can be neglected. Then, its middle layer of the thin plate is looked at as the middle surface. In this case, the relationship between the strains (${\epsilon }_x$, ${\epsilon }_y$, ${\gamma }_{xy}$) and transverse displacement ($w=w\left(x,y,t\right)$) can be represented as Equation (1).

$$\{\begin{array}{l}{\epsilon }_x^{}=-z\frac{{\partial }^{2}w}{\partial x^2}\\ {\epsilon }_y^{}=-z\frac{{\partial }^{2}w}{\partial y^2}\\ {\gamma }_{xy}^{}=-2z\frac{{\partial }^{2}w}{\partial x\partial y}\end{array}$$

(1)

According to Newton’s method, the transverse vibration equation of a thin plate is given by

$$\begin{array}{l}\rho h\ddot{w}=\frac{\partial }{\partial x}\left(\frac{\partial H_{11}}{\partial x}+\frac{\partial H_{12}}{\partial y}\right)+\frac{\partial }{\partial y}\left(\frac{\partial H_{21}}{\partial x}+\frac{\partial H_{22}}{\partial y}\right)\\ H_{11}={\displaystyle \int {\sigma }_{x}zdz,}{H}_{22}={\displaystyle \int {\sigma }_{y}zdz,}{H}_{12}=H_{21}={\displaystyle \int {\tau }_{xy}zdz}\end{array}$$

(2)

where $H_{ij}$ is the bending moment, $\rho$, $h$ are the density and the thickness of the thin plate, respectively, and ${\sigma }_x$, ${\sigma }_y$, ${\tau }_{xy}$ are the stresses.

(i) For the host elastic plate: As a transversely isotropic elastic material, the relationship between the strains and the stresses is as follows [36]:

$$\left[\begin{array}{l}{\sigma }_x^s\\ {\sigma }_y^s\\ {\tau }_{xy}^s\end{array}\right]=\frac{E_s}{1-{\mu }_s^2}\left[\begin{array}{ccc}1& {\mu }_s& 0\\ {\mu }_s& 1& 0\\ 0& 0& \left(1-{\mu }_s\right)/2\end{array}\right]\left[\begin{array}{c}{\epsilon }_x^s\\ {\epsilon }_y^s\\ {\gamma }_{xy}^s\end{array}\right]$$

(3)

Then, its dynamic equation of transverse motion can be expressed as Equation (4) based on Equations (1)–(3).

$$-{\varphi }_s\frac{{\partial }^{2}w}{\partial t^2}=\frac{{\partial }^2}{\partial x^2}\left[D_s\frac{{\partial }^{2}w}{\partial x^2}+{\upsilon }_s\frac{{\partial }^{2}w}{\partial y^2}\right]+\frac{{\partial }^2}{\partial y^2}\left[D_s\frac{{\partial }^{2}w}{\partial y^2}+{\upsilon }_s\frac{{\partial }^{2}w}{\partial x^2}\right]+2\frac{{\partial }^2}{\partial x\partial y}\left[{\chi }_s\frac{{\partial }^{2}w}{\partial x\partial y}\right]$$

(4)

where ${\varphi }_s={\rho }_s{h}_s$, $D_s=E_s{h}_s^3/12\left(1-{\mu }_s^2\right)$, ${\upsilon }_s=D_s{\mu }_s$, ${\chi }_s=D_s\left(1-{\mu }_s\right)$, ${\rho }_s$, $E_s$, ${\mu }_s$, $h_s$ are the density, Young’s modulus, the Poisson’s ratio, and the thickness of the thin plate, respectively [35].

(ii) For the piezoelectric patch: Furthermore, a piezoelectric patch can be assumed to be a transversely isotropic material. Then, the piezoelectric patch is also looked at as a thin plate, and its constitutive equation can be reduced to a 2D form as Equation (5).

$$\left[\begin{array}{c}{\epsilon }_x^p\\ {\epsilon }_y^p\\ {\gamma }_{xy}^p\\ D_3\end{array}\right]=\left[\begin{array}{cccc}{s}_{11}^E& s_{12}^E& 0& d_{31}\\ s_{12}^E& s_{11}^E& 0& d_{31}\\ 0& 0& s_{66}^E& 0\\ d_{31}& d_{31}& 0& e_{33}^T\end{array}\right]\left[\begin{array}{c}{\sigma }_x^p\\ {\sigma }_y^p\\ {\tau }_{xy}^p\\ E_3\end{array}\right]$$

(5)

where $s_{11}^E,s_{12}^E,s_{66}^E=2\left(s_{11}^E-s_{12}^E\right)$ are elastic compliance parameters at the constant electric field, $d_{31}$ is the piezoelectric constant in the 3-1direction, ${\epsilon }_{33}^T$ is the dielectric permittivity atconstant stress, and $D_3$, $E_3$ are the electric displacement and field along the z-axis, respectively. Moreover, Equation (5) can be rewritten as follows:

$$\left[\begin{array}{c}{\sigma }_x^p\\ {\sigma }_y^p\\ {\tau }_{xy}^p\\ D_3\end{array}\right]=\left[\begin{array}{cccc}{\overline{c}}_{11}& {\overline{c}}_{12}& 0& -{\overline{e}}_{31}\\ {\overline{c}}_{12}& {\overline{c}}_{11}& 0& -{\overline{e}}_{31}\\ 0& 0& {\overline{c}}_{66}& 0\\ {\overline{e}}_{31}& {\overline{e}}_{31}& 0& {\overline{\epsilon }}_{33}^s\end{array}\right]\left[\begin{array}{c}{\epsilon }_x^p\\ {\epsilon }_y^p\\ {\gamma }_{xy}^p\\ E_3\end{array}\right]$$

(6)

where, ${\overline{c}}_{11}=s_{11}^E/\left[{(s_{11}^E)}^2-{(s_{12}^E)}^2\right]$, ${\overline{c}}_{12}=-s_{12}^E/\left[{(s_{11}^E)}^2-{(s_{12}^E)}^2\right]$, ${\overline{c}}_{66}=1/s_{66}^E$, ${\overline{e}}_{31}=d_{31}/\left[s_{11}^E+s_{12}^E\right]$, ${\overline{\epsilon }}_{33}^s=e_{33}^T-2d_{31}^2/\left[s_{11}^E+s_{12}^E\right]$.

The piezoelectric patch is assumed to be polarized along the direction of z-axis, and then the electric field in the z-axis can be approximated as Equation (7) [10],

$$E_3=-V/h_p$$

(7)

where $V$ is the output voltage across the two electrodes on the piezoelectric patch, and $h_p$ is the thickness of the piezoelectric patch.

Based on Equation (6), the bending moments can be calculated as follows:

$$\begin{array}{l}{H}_{11}={\displaystyle {\int }_{h_s/2-h_p}^{h_s/2}{\sigma }_x^{p}zdz}=-D_p\left(\frac{{\partial }^{2}w}{\partial x^2}+{\mu }_p\frac{{\partial }^{2}w}{\partial y^2}\right)+{\lambda }_{p}V\\ H_{22}={\displaystyle {\int }_{h_s/2-h_p}^{h_s/2}{\sigma }_y^{p}zdz}=-D_p\left(\frac{{\partial }^{2}w}{\partial y^2}+{\mu }_p\frac{{\partial }^{2}w}{\partial x^2}\right)+{\lambda }_{p}V\\ H_{12}=H_{21}={\displaystyle {\int }_{h_s/2-h_p}^{h_s/2}{\tau }_{xy}^{p}zdz}=-D_p\left(1-{\mu }_p\right)\frac{{\partial }^{2}w}{\partial x\partial y}\end{array}$$

(8)

where, $D_p={\overline{c}}_{11}\left[{\left(h_s/2\right)}^3-{\left(h_s/2-h_p\right)}^3\right]/3$, ${\mu }_p={\overline{c}}_{12}/{\overline{c}}_{11}$, ${\lambda }_p=\left(h_s-h_p\right){\overline{e}}_{31}/2$.

By combining Equations (2), (7), and (8), we can obtain the transverse motion equation of the piezoelectric patch as

$$-{\rho }_p{h}_p\frac{{\partial }^{2}w}{\partial t^2}=\frac{{\partial }^2}{\partial x^2}\left[D_p\left(\frac{{\partial }^{2}w}{\partial x^2}+{\mu }_p\frac{{\partial }^{2}w}{\partial y^2}\right)-{\lambda }_{p}V\right]+\frac{{\partial }^2}{\partial y^2}\left[D_p\left(\frac{{\partial }^{2}w}{\partial y^2}+{\mu }_p\frac{{\partial }^{2}w}{\partial x^2}\right)-{\lambda }_{p}V\right]+2\frac{{\partial }^2}{\partial x\partial y}\left[D_p\left(1-{\mu }_p\right)\frac{{\partial }^{2}w}{\partial x\partial y}\right]$$

(9)

Referring to the idea in Ref. [34], an infinitesimal number ${\lambda }_s$ multiplying the voltage is introduced into Equation (4) in order to unify the forms of Equations (4) and (9). Then, the unified motion equation can be written as

$$-{\varphi }_{*}\frac{{\partial }^{2}w}{\partial t^2}=\frac{{\partial }^2}{\partial x^2}\left[D_{*}\frac{{\partial }^{2}w}{\partial x^2}+{\upsilon }_{*}\frac{{\partial }^{2}w}{\partial y^2}-{\lambda }_{*}V\right]+\frac{{\partial }^2}{\partial y^2}\left[D_{*}\frac{{\partial }^{2}w}{\partial y^2}+{\upsilon }_{*}\frac{{\partial }^{2}w}{\partial x^2}-{\lambda }_{*}V\right]+2\frac{{\partial }^2}{\partial x\partial y}\left[{\chi }_{*}\frac{{\partial }^{2}w}{\partial x\partial y}\right]$$

(10)

where the subscript (*) can be “s” for the thin plate or “p” for the piezoelectric patch, ${\upsilon }_{*}=D_{*}{\mu }_{*}$ and ${\chi }_{*}=D_{*}\left(1-{\mu }_{*}\right)$.

(iii) For the interface circuit: This electric current flowing to the interface circuit can be calculated as Equation (11) using the long-wave approximation assumption [37].

$$I\left(s\right)=D_3{A}_{e}s$$

(11)

where $s$ is the Laplace operator, and $A_e$ is the surface area of the electrode or the piezoelectric patch. On the other hand, the electric current can be written as

$$I\left(s\right)=\frac{V\left(s\right)}{Z\left(s\right)}$$

(12)

where $Z\left(s\right)$ denotes the impedance of the interface circuit.

According to previous studies [20,38], a synchronized charge extraction circuit (SCEC) with a resonant inductor is chosen as the interface circuit, which is shown in Figure 3. In most times, the switch is off, and the capacitor $C_1$ is much less than the clamped capacitor $C_p$ of the piezoelectric patch; hence, we can only consider the resonant inductor as the electric load. In this case, the complex impedance of the interface circuit is $j\omega L$.

By combing Equations (11) and (12), the following formulation can be obtained in the complex-frequency domain.

$$\left(1-\frac{A_e{\epsilon }_{33}^{s}L}{h_p}{\omega }^2\right)V={\omega }^2{\overline{e}}_{31}{A}_{e}L\left({\epsilon }_x^p+{\epsilon }_y^p\right)$$

(13)

Furthermore, Equation (13) can be rewritten as

$$\left({\alpha }_1+{\alpha }_2{\omega }^2\right)V={\omega }^2{\alpha }_3\left(\frac{{\partial }^{2}w}{\partial x^2}+\frac{{\partial }^{2}w}{\partial y^2}\right)$$

(14)

where ${\alpha }_1=1$, ${\alpha }_2=-L{C}_p$, ${\alpha }_3=A_e{\overline{e}}_{31}L{h}_s/2$, $C_p=A_e{\epsilon }_{33}^s/h_p$.

3.2. Bandgap Structures of the PMPIC by an Extended Plane Wave Expansion Method

For the complex 2D wave motion equation as Equation (10), it is difficult to obtainits analytical solution. Several numerical methods have been proposed to investigate 2D phononic crystals, such as the finite difference time domain method, finite element method, and PWE method. Due to its convenience, the PWE method is one of the most widely used methods to calculate bandgap structures of a 2D phononic crystal. As for piezoelectric coupling effects, an extended PWE method is presented here.

Due to a periodic structure, material properties (${\varphi }_{*}$, $D_{*}$, ${\upsilon }_{*}$, ${\lambda }_{*}$ and ${\chi }_{*}$) of the PMPIC are periodic functions of the spatial coordinates ($r$). In Figure 2, it is assumed that ${\Theta }_p$ denotes the set of small elements filled with piezoelectric patch and the center element belongs to ${\Theta }_p$. Then, material properties of the center element can be represented as Equation (15) according to its Fourier series expansion

$$g_0\left(r\right)={\displaystyle \sum _G{g}_0\left(G\right)\exp \left(j\mathrm{Gr}\right)}$$

(15)

where $g_0\left(G\right)$ is the Fourier coefficient, $g_0$ can be any of ${\varphi }_{*}$, $D_{*}$, ${\upsilon }_{*}$, ${\lambda }_{*}$ or ${\chi }_{*}$, $G=n_1{b}_1+n_2{b}_2$ is the reciprocal lattice vector and $b_1,b_2$ are the x-axis and y-axis lattice vectors in the reciprocal space, respectively, and $r=n_1{a}_1+n_2{a}_2$ is the primitive space vector. Here, $b_1=2\pi \left(1,0\right)/a,b_2=2\pi \left(0,1\right)/a$, and $a$ is the lattice constant of the unit cell. Furthermore, $g_0\left(G\right)$ can be defined as

$$g_0\left(G\right)=\frac{1}{S}{\displaystyle {\iint }_S{g}_0\left(r\right)\exp (-j\mathrm{Gr})}{d}^{2}r$$

(16)

where $S$ is the area of the unit cell. Furthermore, the integration of Equation (16) can be calculated as

$$g_0\left(G\right)=\{\begin{array}{l}{g}_{p}f+g_s\left(1-f\right),G=0\\ \left(g_p-g_s\right)P\left(G\right),G\ne 0\end{array}$$

(17)

where $g_p$ and $g_s$ denote material properties (${\varphi }_{*}$, $D_{*}$, ${\upsilon }_{*}$, ${\lambda }_{*}$ or ${\chi }_{*}$) of the piezoelectric patch and host structure, respectively. $f=1/M^2$ is the filling ratio of each small element in the unit cell. $P\left(G\right)$ is called a structural function, which is defined as

$$P\left(G\right)=\frac{1}{S}{\displaystyle \underset{S_p}{\int }{e}^{-j\mathrm{Gr}}{d}^{2}r}$$

(18)

where $S_p$ is the area of the small element. In this paper, the piezoelectric patch in any small element is square, so we will have

$$P\left(G\right)=f\mathrm{sinc}\left(\frac{G_{x}a}{2M}\right)\mathrm{sinc}\left(\frac{G_{y}a}{2M}\right)$$

(19)

where $G_x$ and $G_y$ are x and y components of the reciprocal vector $G$, and sinc(●) is the cardinal sine function.

According to the basic properties of Fourier transform, any other small element belonging to ${\Theta }_p$ will have $g_r\left(G\right)=g_0\left(G\right)\exp \left(j\mathrm{Gr}\right)$. Then, Fourier transform of material constants of the piezoelectric patch in the unit cell can be expressed as

$$g\left(G\right)={\displaystyle \sum _r{g}_r\left(G\right)\delta \left(r\right)}=g_0\left(G\right){\displaystyle \sum _r\exp \left(j\mathrm{Gr}\right)\delta \left(r\right)}$$

(20)

where $r={\left[r_1,\dots ,r_N\right]}^T$, $r_i$ denotes the center position of the ith small element and $N=M^2$. When the ith small element belongs to ${\Theta }_p$, $\delta \left(r_i\right)=1$, otherwise $\delta \left(r_i\right)=0$. Thus $\delta \left(r\right)$ can be looked at as a distribution function of piezoelectric material in the unit cell, which is just the optimization target of this paper.

Based on the Bloch theorem, the displacement and output voltage in Equation (10) can be defined as

$$\begin{array}{l}w\left(r,t\right)={\mathrm{e}}^{j\left(\mathrm{kr}-\omega t\right)}{w}_k\left(r\right)\\ V\left(r,t\right)={\mathrm{e}}^{j\left(\mathrm{kr}-\omega t\right)}{V}_k\left(r\right)\end{array}$$

(21)

where $k$ is the Bloch wave vector which is restricted to the first irreducible Brillouin zone. $w_k\left(r\right)$ is also a periodic function with the same period as that of material properties, so it can be expanded as

$$\begin{array}{l}{w}_k\left(r\right)={\displaystyle \sum _{G^{\prime }}{w}_k\left(G^{\prime }\right)}{e}^{j{G}^{\prime }r}\\ V_k\left(r\right)={\displaystyle \sum _{G^{\prime }}V\left(G^{\prime }\right)}{e}^{j{G}^{\prime }r}\end{array}$$

(22)

Based on Equations (21) and (22), we will have

$$\begin{array}{l}w\left(r,t\right)=e^{-j\omega t}{\displaystyle \sum _{G\prime }{w}_k\left(G^{\prime }\right)}{e}^{j\left(k+G^{\prime }\right)r}\\ V\left(r,t\right)=e^{-j\omega t}{\displaystyle \sum _{G^{\prime }}V\left(G^{\prime }\right)}{e}^{j\left(k+G^{\prime }\right)r}\end{array}$$

(23)

Substituting Equations (15) and (23) into Equation (10), we will have

$$\begin{array}{l}{\omega }^2{\displaystyle \sum _{G^{\prime }}{\varphi }_{*}\left(G^{″}-G^{\prime }\right)}{w}_k\left(G^{\prime }\right)={\displaystyle \sum _{G^{\prime }}{D}_{*}\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{\prime }\right)}_x^2{\left(k+G^{″}\right)}_x^2+{\left(k+G^{\prime }\right)}_y^2{\left(k+G^{″}\right)}_y^2\right]w_k\left(G^{\prime }\right)+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\displaystyle \sum _{G^{\prime }}{\upsilon }_{*}\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{\prime }\right)}_y^2{\left(k+G^{″}\right)}_x^2+{\left(k+G^{\prime }\right)}_x^2{\left(k+G^{″}\right)}_y^2\right]w_k\left(G^{\prime }\right)+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\displaystyle \sum _{G^{\prime }}{\lambda }_{*}\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{″}\right)}_x^2+{\left(k+G^{″}\right)}_y^2\right]V\left(G^{\prime }\right)+\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}2{\displaystyle \sum _{G^{\prime }}{\chi }_{*}\left(G^{″}-G^{\prime }\right)}{\left(k+G^{\prime }\right)}_x{\left(k+G^{\prime }\right)}_y{\left(k+G^{″}\right)}_x{\left(k+G^{″}\right)}_y{w}_k\left(G^{\prime }\right)\end{array}$$

(24)

Similarly, Equation (14) can also be derived as Equation (25).

$$\left[{\displaystyle \sum _{G^{\prime }}{\alpha }_1\left(G^{″}-G^{\prime }\right)}+{\omega }^2{\displaystyle \sum _{G^{\prime }}{\alpha }_2\left(G^{″}-G^{\prime }\right)}\right]V\left(G^{\prime }\right)={\omega }^2{\displaystyle \sum _{G^{\prime }}{\alpha }_3\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{\prime }\right)}_x^2+{\left(k+G^{\prime }\right)}_y^2\right]w_k\left(G^{\prime }\right)$$

(25)

Finally, Equations (24) and (25) can be rearranged as the following standard eigenvalue equation.

$$\left\{\left[\begin{array}{cc}E& F\\ 0& R\end{array}\right]-{\omega }^2\left[\begin{array}{cc}A& 0\\ B& -C\end{array}\right]\right\}\left[\begin{array}{c}{w}_k\left(G^{\prime }\right)\\ V\left(G^{\prime }\right)\end{array}\right]=0$$

(26)

where $\begin{array}{l}E={\displaystyle \sum _{G^{\prime }}{D}_{*}\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{\prime }\right)}_x^2{\left(k+G^{″}\right)}_x^2+{\left(k+G^{\prime }\right)}_y^2{\left(k+G^{″}\right)}_y^2\right]+{\displaystyle \sum _{G^{\prime }}{\upsilon }_{*}\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{\prime }\right)}_y^2{\left(k+G^{″}\right)}_x^2+{\left(k+G^{\prime }\right)}_x^2{\left(k+G^{″}\right)}_y^2\right]\\ +2{\displaystyle \sum _{G^{\prime }}{\chi }_{*}\left(G^{″}-G^{\prime }\right)}{\left(k+G^{\prime }\right)}_x{\left(k+G^{\prime }\right)}_y{\left(k+G^{″}\right)}_x{\left(k+G^{″}\right)}_y\end{array}$, $F={\lambda }_{*}\left(G^{″}-G^{\prime }\right)\left[{\left(k+G^{″}\right)}_x^2+{\left(k+G^{″}\right)}_y^2\right]$, $R={\displaystyle \sum _{G^{\prime }}{\alpha }_1\left(G^{″}-G^{\prime }\right)}$, $A={\displaystyle \sum _{G^{\prime }}{\varphi }_{*}\left(G^{″}-G^{\prime }\right)}$, $C={\displaystyle \sum _{G^{\prime }}{\alpha }_2\left(G^{″}-G^{\prime }\right)}$, $B={\displaystyle \sum _{G^{\prime }}{\alpha }_3\left(G^{″}-G^{\prime }\right)}\left[{\left(k+G^{\prime }\right)}_x^2+{\left(k+G^{\prime }\right)}_y^2\right]$, $G^{\prime }$, $G^{″}$ denote the reciprocal lattice vectors.

Based on Equation (26), we can obtain the corresponding value of $\omega$ once the wave vector $\left(k_x,k_y\right)$ is specified. Generally speaking, the wave vector is always chosen along the edges of the irreducible part of the first Brillouin zone. In this case, the dispersion curves will be obtained, and then vibration bandgaps can be calculated.

4. Binary-Like Topology Optimization of the Unit Cell by Extended Genetic Algorithm

As shown in Figure 2, the distribution of the piezoelectric patch can be represented as a $M\times M$ binary matrix. Thus, topology optimization of the unit cell is equivalent to find the optimal matrix to maximize the bandgaps. Moreover, the binary matrix can be further rewritten as a $M^2$-dimension vector, which can be looked at as encoding the matrix. In this case, the problem of the mentioned binary-like topology optimization will be solved by the genetic algorithm.

According to the previous work [20], the PVEH performance can be improved greatly once the $L-C_p$ resonant frequency is equal to the excitation frequency. For engineering applications, the central frequency of target structural vibration is assumed to be $f_0$, and then we can choose the inductor as follows:

$$L=\frac{1}{4{\pi }^2{C}_p{f}_0^2}$$

(27)

It is obvious that $L$ depends on $C_p$. Furthermore, $C_p$ is related to the distribution of the piezoelectric patch. Thus, it is valuable that both the piezoelectric patch and the interface circuit can be topologically optimized simultaneously to realize high-efficiency and broadband PVEH in structures. The detailed procedure is described as follows.

(1) Definition of the objective function and the fitness: In order to achieve broadband and low-frequency vibration energy harvesting, we hope to maximize the width of the first-order bandgap. Thus, the optimization objective function can be defined as

$$B\left(Q\right)=\underset{Q}{\min }{\omega }_2\left(Q,k\right)-\underset{Q}{\max }{\omega }_1\left(Q,k\right)$$

(28)

where ${\omega }_1,{\omega }_2$ denote the first-order and second-order dispersion curves calculated by Equation (26), and $Q$ is the $M\times M$ binary matrix representing the topology distribution of the piezoelectric patch. Then, the optimization variables are the elements of $Q$, which are defined as

$$Q_{ij}=\{\begin{array}{l}0,isthehoststructure\\ 1,isthepiezoelectricmaterial\end{array},\hspace{0.17em}\hspace{0.17em}i,j=1,\dots ,M$$

(29)

Different from previous studies on vibration isolations, all elements of “1” digits in $Q$ should be connected as a whole in order to ensure an integrated piezoelectric patch in this paper. Otherwise, Equation (9) does not hold again. By coincidence, this task can be easily overcome by the algorithm of finding the number of islands (NIs) in the {0, 1} matrix ($Q$). At the same time, the fitness of each $Q$ is defined as the width of the first-order bandgap ($B$).

(2) Initialization: Initialization parameters include geometric and material properties of the thin plate and the piezoelectric patch, the population size, the iterative times, the selection rate, the crossover probability, the mutation probability, and $M$. Here, it should be noted that each individual in the population must be selected to make its NIs be equal to 1.

(3) Selection: For each individual in the population, the first-order and second-order dispersion curves are simulated by Equation (26), and then the first-order band gap is calculated to define its fitness. If the first-order bandgap does not exist, a small number is assigned as its fitness to keep this individual. Based on the roulette wheel selection algorithm, the number (e.g., the population size × the selection rate) of individuals in the initialization parent population is selected to construct a new mating pool.

(4) Encoding and crossover: Firstly, we choose randomly two individuals in the mating pool as the parents and transform them into a $M^2$-dimension vector. Then, the selected two individuals are mated to create two offspring for replacing their parents according to a crossover probability ($P_{cr}$). Similarly, each offspring will be checked whether its NIs areequal to 1 or not. If not, the crossover operation needs to be carried out again. Inthis way, an offspring population will be created.

(5) Mutation and decoding: Each individual in the offspring population is mutated according to the mutation probability $P_{mu}$. The mutation position is chosen randomly, where the element value switches to its opposite. Additionally, the mutated individual will be checked whether its NIs are equal to 1 or not. If not, the mutation operation needs to be carried out again. In this way, a mutated offspring population will be created. Next, each individual in the mutated offspring population is rearranged as a $M\times M$ binary matrix. The resulted offspring population will be merged with the parent population to create the new parent population.

(6) Termination: Repeating steps (3) to (5) until the iterative times are reached. Then, we can obtain the optimal first-order bandgap width and the corresponding optimal binary matrix. Based on the binary matrix, the optimal topology of the unit cell is realized.

In summary, the binary-like topology optimization flowchart using the extended GA algorithm is plotted in Figure 4.

5. Simulations and Discussions

In order to validate the proposed method, numerical simulations are carried out to carry out the binary-like topology optimization procedure using the extended GA algorithm. Geometrical and material parameters in simulations are listed in

Table 1. Geometrical parameters of the unit cell.

$\mathit{a}$ (mm)	${\mathit{h}}_{\mathit{s}}$ (mm)	${\mathit{h}}_{\mathit{p}}$ (mm)
100	1	0.2

and

Table 2. Material parameters of the unit cell.

Parameters	PZT-5H	Al
Density $\rho$ (kg/m3)	7500	2700
Young’s modulus E (Pa)	63 × 109	70 × 109
Poisson’s ratio $\mu$	0.34	0.33
Elastic compliance $s_{11}$ (m2/N)	16.5 × 10−12	/
Elastic compliance $s_{12}$ (m2/N)	−4.78 × 10−12	/
Elastic compliance $s_{66}$ (m2/N)	42.6 × 10−12	/
Piezoelectric constant $d_{31}$ (C/N)	−274 × 10−12	/
Dielectric permittivity ${\epsilon }_{33}^T$ (F/m)	30.1 × 10−9	/

, respectively.

5.1. Feasibility of the Proposed Algorithm

The unit cell with configuration A, as shown in Figure 5a, is considered for the first numerical example and M is randomly selected as 10. Then, it is divided into 10 × 10 small elements. For this configuration, its finite element model (FEM) is built by the Comsol software and here the inductor is chosen as $L=1$ H. Then, the first and second dispersion curves are calculated by the proposed algorithm and the FEM simulation, respectively. The results are shown in Figure 6, and we can see that (i) numerical results are close to those of FEM simulations; hence, the proposed bandgap algorithm is testified to be effective; (ii) there is no full bandgaps under configuration A, while local bandgaps exist along the directions of $X\Gamma$ and $\Gamma M$; (iii) the bias between the numerical and FEM results may mainly come from the approximation error of the long-wave hypothesis in Equation (11).

5.2. Effects of the Interface Circuit (L)

Furthermore, the unit cell with configuration B as Figure 7 is considered for the second numerical example, and M is randomly selected as 12. Then, it is divided into 12 × 12 small elements.

Firstly, the inductor is chosen as $L=1$ H, and the first four disperse curves are calculated as Figure 8a. It can be seen that (i) besides local bandgaps, there is a full bandgap, which is called the first bandgap. The bandgap width is equal to 7.99; (ii) compared with configuration A, it may be easier to generate full bandgaps by using a nonrectangular configuration of piezoelectric material. Indeed, it provides the basis of topology optimization in this paper.

Then, in order to investigate the effects of the interface circuit, different values of the inductor ($L=1\mathrm{H},5\hspace{0.17em}\mathrm{H},10\mathrm{H}$) are chosen. The first-order bandgaps are calculated and shown in Figure 8b, respectively. It can be seen that (i)the first bandgap will shift to the low-frequency region with the increase of $L$; therefore, large $L$ will benefit low-frequency vibration energy harvesting; (ii) with the increase of $L$, the width of the first bandgap firstly increases and then decreases. It means that there may be an optimal $L$ to obtain the maximum width of the first bandgap, which also will indicate the feasibility of topology optimization.

5.3. Topology Optimization of the PMPIC

In order to validate the topology optimization algorithm, the procedure shown in Figure 4 is carried out by MATLAB Software. Firstly, the inductor is fixed as $L=1$ H, and the unit cell is divided into 10 × 10 small elements (e.g., $M=10$). Parameters of the GA are set as the population size, the selection rate, the crossover probability, and the mutation probability, which are equal to 50, 0.5, 0.7, and 0.001, respectively.

In order to show the structural topology clearly, two colors are applied for different materials. The Blue area denotes the host base structure, and the yellow area denotes the piezoelectric patch. The optimal topology configuration of the unit cell is shown in Figure 9a. Moreover, the first four dispersion curves are plotted in Figure 9b, and the maximum bandgap width is equal to 36.6891 Hz ([124.1, 160.8] Hz). Obviously, the maximum bandgap width is much larger than that of configuration B, and the structural topology of the piezoelectric patch is nonrectangular.

Secondly, it is desired to optimize the structural topology and the interface circuit simultaneously. Here, the central frequency is assumed to be 200 Hz, and the unit cell is also divided into 10 × 10 small elements (e.g., $M=10$). The optimal topology configuration of the unit cell is shown in Figure 10a, the first four dispersion curves are plotted in Figure 10b, and the maximum bandgap width is equal to 22.4 Hz ([192.3, 214.7]Hz). At the same time, the inductor is calculated as $L=1.221\mathrm{H}$. Obviously, $f_0=200\mathrm{Hz}$ is close to the center of the optimal bandgap. Four isolated host structures in the optimal result can be used for local masses.

In the end, it should be emphasized that M is a key parameter and the precision of material distribution increases with M. However, large M will cause the amount of optimization calculation to increase rapidly. Without loss of feasibility, this paper pays less attention on calculation efficiency. As for this problem, it may be overcome by using advanced calculation tools such as supercomputing platforms or fast optimization algorithms.

6. Conclusions

In this paper, a self-powered SCEC is used as the interface circuit connected to the PMP for harvesting structural vibration energy. In order to optimize the topology structure of the PMPIC, a binary-like topology optimization process is first presented. Then, simulations are conducted to validate the proposed method. The main highlights of this paper may include the following: (i) A binary-like topology optimization approach is first proposed for the PMPIC by dividing the unit cell into a {0, 1} matrix indicating material distribution so that the coupled structure and circuit can be optimized simultaneously; (ii) an elastic–electro–mechanical coupling model is built for the coupled PMP and self-powered SCEC with a resonant inductor; (iii) an extended plane wave expansion method is presented to model the bandgap character of the PMPIC, which is numerically solved by the Bloch theorem; (iv) GA is extended to perform topology optimization of the PMPIC by adding the restriction of NIs. In the future, a better approximation of Equation (11) can be studied to replace the long-wave hypothesis in order to increase its accuracy, and fast optimization algorithms need to be explored to reduce the computation cost. In addition, more experimental tests need to be performed to make sure theoretic results. For real-world experiment design, some uncertainties should be addressed and explored in the future, such as fabrication error, nonuniformity of piezoelectrical material, model approximation error, the time delay of the self-powered switch, etc.