The optimal planform of a cantilever unimorph piezoelectric vibrating energy harvester (PVEH) with a device-layer edge block

The study considers the optimal planform of a cantilever piezoelectric vibrating energy harvester (PVEH) with an edge block that is patterned only in the device layer. The optimal response of a PVEH is achieved when the strain in the piezoelectric layer is uniform. It is possible to design the planform of the PVEH such that when it vibrates at a given frequency of excitation vibrations, the mode shape results in a uniform strain in the piezoelectric layer. Another design choice is the size of the edge block. The size of the edge block affects the natural frequency of the structure and hence may be adjusted such that the natural frequency of the PVEH matches the frequency of the excitation vibrations. Previous studies considered the optimal planform of a PVEH with a massive edge block that is patterned in both the device and handle layers of the wafer. For this case, it was found that the optimal planform of such a device is trapezoidal. However, there is another class of PVEH devices in which the edge block is patterned only in the device layer (i.e., a device-layer edge block). In the present study, we show that the optimal planform for such a PVEH is defined by Bessel functions, and we demonstrate the predictive capabilities of our analytic model by comparison to results of finite-element simulations.

1 Introduction

Piezoelectric vibrating energy harvester (PVEH) devices have been extensively studied because of their potential to power autonomous sensors. The design of such devices varies extensively depending on the frequency distribution of the environmental vibrations. The simplest case is when the frequency of the environmental vibrations is fixed and a priori known. In this case, it is possible to optimize a cantilever PVEH such that its natural frequency is compatible with the excitation vibrations. In the context of PVEHs that target a specific excitation frequency, it is often assumed that the optimal performance is achieved when the strain amplitude in the piezoelectric layer is uniform (Salman et al., 2024b, 2025).

In typical cantilever PVEH devices that are made using silicon-on-insulator (SOI) wafer technology (Cowen et al., 2014; Du et al., 2017; Jia and Seshia, 2016; Pillai et al., 2019), the cantilever is patterned only in the device layer, whereas the edge block is patterned in both the device and handle layers of the wafer. A massive edge block will increase the amplitude of vibrations and therefore may enhance the performance of the system. However, if the excitation vibrations are at a sufficiently high frequency, a massive edge block may be detrimental, and the edge block is patterned only in the device layer. This simplifies the fabrication process and may also simplify packaging because the cantilever can vibrate freely and it is not necessary to prepare extra space to allow for the vibrations of a handle-layer edge block.

In a recent paper (Salman et al., 2024a), we considered a cantilever PVEH with a massive edge block that extends into the handle layer. The inertia of this edge block dominates the vibration response of the structure. We showed that for such a massive edge block, the optimal planform of the PVEH is a trapeze. The same conclusion was presented in many previous studies which were based on experiments, simulations or a combination of the two (Ben Ayed et al., 2009, 2014; Baker et al., 2005; Benasciutti et al., 2010; Chen et al., 2009; Dietl and Garcia, 2010; Gallina and Benasciutti, 2013; Goldschmidtboeing and Woias, 2008; Halvorsen and Dong, 2008; Hosseini and Hamedi, 2016; Jia and Seshia, 2016; Lee et al., 2009; Mateu and Moll, 2005; Matova et al., 2013; Miller et al., 2008, 2011; Muthalif and Nordin, 2015; Park et al., 2012; Raju et al., 2018; Rosa and De Marqui Junior, 2014; Roundy et al., 2005; Salmani et al., 2015; Yang et al., 2009; Zhang et al., 2017). However, our work was the first rigorous analysis of the problem, and we presented an explicit functional form of the optimal trapeze planform (Salman et al., 2024a). That model was the first ever to offer predictive capabilities. In contrast to those devices, there is a different class of PVEHs in which the edge block is patterned only in the device layer (e.g., Fig. 1) and does not extend into the underlaying handle layer. In this case, the inertia of the beam is as significant as the inertia of the edge block.

In the present study, we show that the optimal planform of a cantilever PVEH with a device-layer edge block is not a trapeze with straight contours but rather a planform with curved contours described by Bessel functions. We validate our model and demonstrate its predictive capabilities using finite-element simulations.

2 Analysis

Figure 1 presents a schematic description of a cantilever with an edge block that is patterned only in the device layer. The beam has a length L, a uniform thickness hć and a rectangular cross section with varying width w(x). The edge block at the far edge has a mass m_L and a moment of inertia I_L. The moment of inertia I_L refers to the z axis, and it is given relative to the center point (i.e., $y=z=$ 0) of the edge cross section x=L. The distance from that point to the center of mass of the edge block is d in the x direction.

Figure 1 Schematic description of a PVEH device, which is patterned only in the device layer. The PVEH is constructed from a cantilever and an edge block. The cantilever has a length L, a uniform thickness h and a rectangular cross section with a width w(x) which changes along the x axis. The edge block has a length L_EB and a width w_EB, and it is connected at x=L.

It is assumed that the Euler–Bernoulli beam theory is applicable, where the governing equation that determines the out-of-plane displacement y(t,x) in free vibrations, and is given by

$$\begin{array}{}\text{(1)}& \mathit{\rho }hw\left(x\right){\displaystyle \frac{{\partial }^{\mathrm{2}}y(t,x)}{\partial t^{\mathrm{2}}}}=-E{\stackrel{\mathrm{̃}}{I}}_{yy}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial x^{\mathrm{2}}}}\left(w\left(x\right){\displaystyle \frac{{\partial }^{\mathrm{2}}y(t,x)}{\partial x^{\mathrm{2}}}}\right).\end{array}$$

Here, ρ is the density of the elastic material and E is the Young modulus. The term $w\left(x\right)E{\stackrel{\mathrm{̃}}{I}}_{yy}$ is the effective bending rigidity of the beam cross section, where ${\stackrel{\mathrm{̃}}{I}}_{yy}$ is given by

$$\begin{array}{}\text{(2)}& {\stackrel{\mathrm{̃}}{I}}_{yy}=\underset{y=-h/\mathrm{2}}{\overset{h/\mathrm{2}}{\int }}{y}^{\mathrm{2}}\mathrm{d}y={\displaystyle \frac{\mathrm{1}}{\mathrm{12}}}{h}^{\mathrm{3}}.\end{array}$$

We assume that the vibrations of the cantilever can be described as the product of a time-periodic function T(t) and a deflection mode Y(x):

$$\begin{array}{}\text{(3)}& y(t,x)=T\left(t\right)\cdot Y\left(x\right).\end{array}$$

The amplitude of the axial strain at the top surface of the beam is therefore given by

$$\begin{array}{}\text{(4)}& {\mathit{\epsilon }}_{xx}\left(x\right)=-{\displaystyle \frac{d^{\mathrm{2}}Y\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}}{\displaystyle \frac{h}{\mathrm{2}}}.\end{array}$$

The aim here is to determine the functional form of w(x) such that the strain on the top surface is uniform (${\mathit{\epsilon }}_{xx}\left(x\right)=\stackrel{\mathrm{‾}}{\mathit{\epsilon }}$). Substituting this into Eq. (4) yields

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{‾}}{\mathit{\epsilon }}=-{\displaystyle \frac{d^{\mathrm{2}}Y\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}}{\displaystyle \frac{h}{\mathrm{2}}}.\end{array}$$

From Eq. (5), we can deduce that the curvature along the beam must be uniform. It follows that

$$\begin{array}{}\text{(6)}& Y\left(x\right)=-{\displaystyle \frac{\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}{x}^{\mathrm{2}}{\displaystyle \frac{\mathrm{d}Y\left(x\right)}{\mathrm{d}x}}=-{\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}x{\displaystyle \frac{d^{\mathrm{2}}Y\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}}=-{\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}},\end{array}$$

where the boundary conditions Y(0) = 0 and Y'(0) = 0 have been imposed.

Substituting Eqs. (3) and (6) into Eq. (1) yields

$$\begin{array}{}\text{(7)}& \mathit{\rho }w\left(x\right)\stackrel{\mathrm{‾}}{\mathit{\epsilon }}{x}^{\mathrm{2}}{\displaystyle \frac{d^{\mathrm{2}}T\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}=-E{\stackrel{\mathrm{̃}}{I}}_{yy}{\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}T\left(t\right){\displaystyle \frac{d^{\mathrm{2}}w\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}},\end{array}$$

which may be rearranged in the form

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{\mathrm{1}}{T\left(t\right)}}{\displaystyle \frac{d^{\mathrm{2}}T\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}=-E{\stackrel{\mathrm{̃}}{I}}_{yy}{\displaystyle \frac{\mathrm{2}}{\mathit{\rho }h}}{\displaystyle \frac{\mathrm{1}}{w\left(x\right)\cdot x^{\mathrm{2}}}}{\displaystyle \frac{d^{\mathrm{2}}w\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}}=-{\mathit{\omega }}^{\mathrm{2}}.\end{array}$$

It follows that the time harmonic solution must satisfy

$$\begin{array}{}\text{(9)}& {\displaystyle \frac{d^{\mathrm{2}}T\left(t\right)}{\mathrm{d}{t}^{\mathrm{2}}}}+{\mathit{\omega }}^{\mathrm{2}}T\left(t\right)=\mathrm{0},\end{array}$$

and without loss of generality, we may consider the solution

$$\begin{array}{}\text{(10)}& T\left(t\right)=\cos \left(\mathit{\omega }t\right).\end{array}$$

Accordingly, the cantilever width w(x) must satisfy the following differential equation:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{d^{\mathrm{2}}w\left(x\right)}{\mathrm{d}{x}^{\mathrm{2}}}}-{\displaystyle \frac{\mathit{\rho }h}{\mathrm{2}E{\stackrel{\mathrm{̃}}{I}}_{yy}}}{\mathit{\omega }}^{\mathrm{2}}{x}^{\mathrm{2}}w\left(x\right)=\mathrm{0}.\end{array}$$

Due to the edge block accelerations, a resultant transverse force and a resultant moment are applied at the beam's far edge, as illustrated in Fig. 2.

Figure 2 The loads applied by the device-layer edge block on the far edge of the cantilever beam.

From Eq. (6), it follows that at the far edge of the beam, the transverse and angular accelerations of the edge block are given by

$$\begin{array}{}\text{(12)}& {\displaystyle }{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}y}{\partial t^{\mathrm{2}}}}\right|}_{x=L}={\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}{\mathit{\omega }}^{\mathrm{2}}{\displaystyle \frac{L^{\mathrm{2}}}{\mathrm{2}}}\cos \left(\mathit{\omega }t\right),\text{(13)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial t^{\mathrm{2}}}}\left({\left.{\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}x}}\right|}_{x=L}\right)={\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}{\mathit{\omega }}^{\mathrm{2}}L\cos \left(\mathit{\omega }t\right).\end{array}$$

The resultant force and moment at the far edge x=L that are applied by the edge block and induced due to the edge block transverse and angular accelerations are given by

$$\begin{array}{}\text{(14)}& {\displaystyle }{V}_L=-m_L\left[{\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}y}{\partial t^{\mathrm{2}}}}\right|}_{x=L}+d\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial t^{\mathrm{2}}}}\left({\left.{\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}x}}\right|}_{x=L}\right)\right],\text{(15)}& {\displaystyle }{M}_L=-\left[I_L{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial t^{\mathrm{2}}}}\left({\left.{\displaystyle \frac{\mathrm{d}Y}{\mathrm{d}x}}\right|}_{x=L}\right)+m_L{\left.d\cdot {\displaystyle \frac{{\partial }^{\mathrm{2}}y}{\partial t^{\mathrm{2}}}}\right|}_{x=L}\right].\end{array}$$

At the far edge, the resultant shear force and the resultant bending moment are given by

$$\begin{array}{}\text{(16)}& {\displaystyle }\begin{array}{rl}{V}_L& =-E{\stackrel{\mathrm{̃}}{I}}_{yy}{\left.{\displaystyle \frac{\partial }{\partial x}}\left(w\left(x\right){\displaystyle \frac{{\partial }^{\mathrm{2}}y(t,x)}{\partial x^{\mathrm{2}}}}\right)\right|}_{x=L}\\ & =E{\stackrel{\mathrm{̃}}{I}}_{yy}{\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}\cos \left(\mathit{\omega }t\right){\left.{\displaystyle \frac{\partial w\left(x\right)}{\partial x}}\right|}_{x=L},\end{array}\text{(17)}& {\displaystyle }\begin{array}{rl}{M}_L& =E{\stackrel{\mathrm{̃}}{I}}_{yy}{\left.\left(w\left(x\right){\displaystyle \frac{{\partial }^{\mathrm{2}}y(t,x)}{\partial x^{\mathrm{2}}}}\right)\right|}_{x=L}\\ & =-E{\stackrel{\mathrm{̃}}{I}}_{yy}{\displaystyle \frac{\mathrm{2}\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}{h}}\cos \left(\mathit{\omega }t\right)w\left(L\right).\end{array}\end{array}$$

After substituting Eqs. (12) and (13) into Eq. (14) and comparing to Eq. (16), it follows that

$$\begin{array}{}\text{(18)}& {\left.{\displaystyle \frac{\mathrm{d}w\left(x\right)}{\mathrm{d}x}}\right|}_{x=L}=-{\displaystyle \frac{m_L{\mathit{\omega }}^{\mathrm{2}}}{E{\stackrel{\mathrm{̃}}{I}}_{yy}}}\left[{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{L}^{\mathrm{2}}+d\cdot L\right].\end{array}$$

Substituting Eqs. (12) and (13) into Eq. (15) and comparing to Eq. (17) yields

$$\begin{array}{}\text{(19)}& w\left(L\right)={\displaystyle \frac{{\mathit{\omega }}^{\mathrm{2}}L}{E{\stackrel{\mathrm{̃}}{I}}_{yy}}}\left[I_L+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{m}_L\cdot d\cdot L\right].\end{array}$$

Now, Eq. (11) can be solved analytically from x=L back to x= 0 beginning with the two edge conditions for w(x) given by Eqs. (18) and (19). The analytic result of this is given by

$$\begin{array}{}\text{(20)}& \begin{array}{rl}w\left(x\right)& =a_{\mathrm{1}}\sqrt{x}\cdot I_{\mathrm{0.25}}\left(x^{\mathrm{2}}{\displaystyle \frac{\mathit{\omega }}{\mathrm{2}}}\sqrt{{\displaystyle \frac{\mathit{\rho }h}{\mathrm{2}E{\stackrel{\mathrm{̃}}{I}}_{yy}}}}\right)+a_{\mathrm{2}}\sqrt{x}\\ & \cdot K_{\mathrm{0.25}}\left(x^{\mathrm{2}}{\displaystyle \frac{\mathit{\omega }}{\mathrm{2}}}\sqrt{{\displaystyle \frac{\mathit{\rho }h}{\mathrm{2}E{\stackrel{\mathrm{̃}}{I}}_{yy}}}}\right).\end{array}\end{array}$$

Here, I_0.25 and K_0.25 are modified Bessel functions, and the constants a₁ and a₂ are determined from the boundary conditions (18) and (19).

3 Simulations

We consider a typical cantilever PVEH with a device-layer edge block. The edge block has a width of w_EB= 200 µm and various lengths, L_EB= 100, 200, 300 and 400 µm. We target a vibrations frequency ω= 30 000 rad s⁻¹ (i.e., f= 4.774 kHz), the cantilever length is L= 2000 µm, and the width w(x) varies along the beam. The considered PVEH devices are patterned only in the device layer, which has a uniform thickness h= 10 µm. For mechanical properties, we consider devices that are made using SOI wafer technology from single-crystalline silicon (SCS) in the (110) orientation, where the density is ρ= 2330 kg m⁻³ and the material moduli are E_xx= 169.7 GPa, G_xy= 80 GPa and ν_xz= 0.0606. For each of the edge block lengths L_EB, the constants a₁ and a₂ were computed by applying the boundary conditions in Eqs. (18)–(19), using MATLAB (R2023a). These constants were substituted into Eq. (20) to compute the beam width w(x). An example of an optimal Bessel tapered planform is presented in Fig. 3.

Figure 3 The optimal Bessel tapered planform for an edge block length of L_EB= 400 µm. The blue lines mark the beam width w(x), and the orange lines mark the edge block.

To validate our model, we simulate the eigenfrequency for each PVEH using COMSOL^™ 6.0 finite-element code (COMSOL, 2022). Each simulation includes a Bessel tapered planform with a matching device-layer edge block. Figure 4 presents a typical geometry and mesh. In the simulations, we used quadratic serendipity quadrilateral elements, with a restriction on each element size (i.e., minimum and maximum length). It was verified that the simulation results converge with mesh refinement.

Figure 4 A typical geometry and mesh used in the finite-element simulations. The zoomed-in view focuses on the connection of the edge block to the beam.

The aim of our work is to optimize the planform of the PVEH such that the strain distribution over the top surface of the cantilever is uniform. To this end, we define the measure of nonuniformity $S_{{\mathit{\epsilon }}_{xx}}$ of the axial strain over the top surface of the cantilever ε_xx by

$$\begin{array}{}\text{(21)}& S_{{\mathit{\epsilon }}_{xx}}={\displaystyle \frac{\sqrt{\frac{\underset{x=\mathrm{0}}{\overset{L}{\int }}{\mathit{\epsilon }}_{xx}^{\mathrm{2}}w\left(x\right)\mathrm{d}x}{\underset{x=\mathrm{0}}{\overset{L}{\int }}w\left(x\right)\mathrm{d}x}}-\frac{\underset{x=\mathrm{0}}{\overset{L}{\int }}\sqrt{{\mathit{\epsilon }}_{xx}^{\mathrm{2}}}w\left(x\right)\mathrm{d}x}{\underset{x=\mathrm{0}}{\overset{L}{\int }}w\left(x\right)\mathrm{d}x}}{\frac{\underset{x=\mathrm{0}}{\overset{L}{\int }}\sqrt{{\mathit{\epsilon }}_{xx}^{\mathrm{2}}}w\left(x\right)\mathrm{d}x}{\underset{x=\mathrm{0}}{\overset{L}{\int }}w\left(x\right)\mathrm{d}x}}}.\end{array}$$

The measure of nonuniformity in Eq. (21) is identically zero (i.e., $S_{{\mathit{\epsilon }}_{xx}}=$ 0) only when the axial strain over the top surface of the cantilever ε_xx is uniform (i.e., ε_xx is not dependent on x or z).

For each length of edge block L_EB, Table 1 presents the width of the beam at the clamped edge w(0) and the far edge w(L), the simulated eigenfrequency ω_n and the simulated nonuniformity of axial strain $S_{{\mathit{\epsilon }}_{xx}}$.

Table 1 The optimal parameters for different device-layer edge blocks and the simulated axial strain nonuniformity.

The values of the simulated frequency ω_n in Table 1 validate our model, which considered the target frequency of ω= 30 000 rad s⁻¹. The relative errors in ω_n suggest that there is a discrepancy between the analysis and the simulations. Our analysis considers the Euler–Bernoulli beam theory, which neglects shear deformations that are not neglected in the finite-element simulations, resulting in a discrepancy of 1 %.

The low values of the nonuniformity $S_{{\mathit{\epsilon }}_{xx}}$ in Table 1 suggest that the strain is indeed rather uniform. To examine the scale of these low-value nonuniformities, we compare the optimal Bessel tapered planforms to non-optimal planforms. As a test case, we consider the third entry in Table 1, with L_EB= 300 µm.

Figure 5 shows the strain distribution over the top surface of the optimal Bessel tapered beam for the third entry in Table 1. For this optimal Bessel tapered beam, the related nonuniformity is $S_{{\mathit{\epsilon }}_{xx}}=$ 3.037 × 10⁻⁴, and it is clear from Fig. 5 that the strain over the top surface is rather uniform.

Figure 5 The axial strain distribution of the optimal Bessel tapered planform for edge block length L_EB= 300 µm. The uniform color suggests that the strain is rather uniform. The reddish corners at the clamped edge suggest that some strain concentration occurs at these corners (zoomed in). The color scale is between red (maximal strain) and blue (minimal strain) with green at the average of maximal and minimal strains.

For comparison, we consider a non-optimal cantilever beam for the same edge block, in which the Bessel edges are replaced with straight lines. This beam consists of the same widths at the clamped and the far edge (i.e., same w(0) and w(L)), but this is a width that is linearly tapered.

Figure 6 shows the axial strain distribution over the top surface of the trapezoidal beam. For this non-optimal trapezoidal beam, the nonuniformity is $S_{{\mathit{\epsilon }}_{xx}}=$ 1.643 × 10⁻², which is 50 times larger than the value obtained for the optimal planform.

Figure 6 The axial strain distribution of a non-optimal trapezoidal planform for edge block length L_EB= 300 µm. The color scale is between red (maximal strain) and blue (minimal strain) with green at the average of maximal and minimal strains.

For the other cases of L_EB= 100, 200 and 400 µm, the strains of the Bessel tapered planforms were more uniform relative to the related trapezoidal planforms by a factor of ∼ 580, ∼ 175 and ∼ 18, respectively.

4 Discussion and summary

To simplify our model, we assume that the device-layer edge block responds like a rigid body. This facilitates the calculation of the resultant shear force and resultant moment at the edge cross section, where the edge block is connected to the cantilever. If, however, the edge block length and width are significantly large such that deformation of the edge block is substantial, then deriving the relevant resultant shear force and resultant moment may require more rigorous analysis.

For the sufficiently small edge blocks considered in the present study, our simulations confirm that the edge block deformation is marginal, and therefore the beam planform predicted by our model results in very small simulated strain nonuniformities in the piezoelectric layer.

In any case, the strain nonuniformities predicted by our model are small, and the difference between a trapezoidal planform and a Bessel function planform may not be large. However, in previous studies it was identified, through simulations, that in some cases with a nondominant edge block, the performance of a PVEH with curved contours is superior (Ben Ayed et al., 2014; Dietl and Garcia, 2010; Park et al., 2012; Salmani et al., 2015). The present study provides an analysis that explains why curved contours may be superior, and provides a model with predictive capabilities.

This study presents a derivation of the optimal planform of a PVEH with a device-layer edge block. It is shown that the optimal planform of such a device is defined by Bessel functions. Using our model, it is possible to design a PVEH with an optimal performance such that the amplitude of the strain in the piezoelectric layer is uniform. Finite-element simulations show good agreement with the analytic derivation, validating the presented model.

The relevance of this model is that it provides a tool with predictive capabilities for designing PVEH devices with optimal performance.