Optimization of the Performance of Bimorph Piezoelectric Energy Harvesters with Honeycomb Metamaterials Using Artificial Neural Networks and NSGA-II

Mohamed Taha Mhiri^{1, 2}, Walid Larbi^2,*, Mnaouar Chouchane¹ and Mohamed Guerich³

¹Université de Monastir, Ecole Nationale d’Ingénieurs de Monastir, Laboratoire de Génie Mécanique, LR99ES32, 5019, Monastir, Tunisie
²Structural Mechanics and Coupled Systems Laboratory, Conservatoire National des Arts et Métiers, 292 rue Saint Martin, 75003, Paris, France
³De Vinci Higher Education, De Vinci Research Center, Paris, France
E-mail: walid.larbi@lecnam.net
*Corresponding Author

Received 30 September 2025; Accepted 15 December 2025

Abstract

Honeycomb-based metamaterials have recently attracted considerable interest for their potential in energy-harvesting applications. In this paper, we focus on two types of bimorph harvesters incorporating honeycomb substrates with positive and negative Poisson’s ratios (PPR and NPR). The substrates are coupled with piezoceramic (PZT) layers to enable electromechanical conversion. To improve energy-harvesting efficiency, we propose an optimization framework that integrates finite element (FE) simulations for data generation, a pretrained neural network for rapid performance prediction and the NSGA-II evolutionary algorithm for multiobjective optimization. The proposed strategy enables the identification of optimal geometric parameters of the honeycomb cells, particularly for cantilever resonators where the fundamental vibration modes are critical and must be tuned to specific eigenfrequencies for engineering applications. The optimization is applied to an initially unoptimized bimorph harvester with a honeycomb substrate. The simulation results indicate that an optimized harvester with an NPR metamaterial substrate can increase the power-to-mass ratio at a resonance frequency of 160 Hz by approximately 14.72% compared to the unoptimized honeycomb harvester. Moreover, compared to a harvester with a solid substrate, the improvement reaches 5%. Furthermore, a life cycle analysis was conducted, showing that the honeycomb substrate can significantly increase the operational lifetime of the piezoceramic layer. This improvement arises because the honeycomb substrate reduces the equivalent stress on the active layer. At high acceleration levels, the predicted operational lifetime is enhanced by more than tenfold compared to the solid substrate, since the PZT material is brittle and lacks flexibility. This improvement is notable because it surpasses the typical lifespan of both standard button batteries and rechargeable batteries highlighting a strong potential of this design for future applications.

Keywords: Piezoelectric vibration energy harvesting, metamaterial, artificial neural network, NSGA-II.

1 Introduction

In the Internet of Things (IoT) society, it is essential to build wireless sensor networks that can efficiently transmit and receive data collected by sensors through wireless communication [1]. Consequently, a large number of wireless sensor nodes are deployed, making their power supply a critical concern. Most of these sensors are powered by conventional solutions such as button batteries, which can typically support devices requiring 1 mW or less [2]. However, the limited lifetime and maintenance requirements of such batteries have driven research toward alternative energy sources. Over the past decade, piezoelectric energy harvesting (PEH) has attracted considerable attention due to its capability to provide a sustainable power supply for low-power sensors.

One of the most commonly employed PEH configurations is the cantilever beam, which consists of a piezoelectric layer bonded to a solid elastic substrate. Typically, there are two types of cantilever harvesters. If the harvester uses a single piezoelectric layer it is known as an unimorph configuration [3]. In contrast, a bimorph configuration consists of two piezoelectric layers bonded to a solid elastic substrate [4].

To improve the performance of PEHs and bring them closer to conventional batteries, researchers have focused on boosting their output power. This is particularly critical because non-rechargeable button batteries, still widely used, offer a relatively high power-to-mass ratio of about 220 mW/g, despite a lifespan of less than 5 years [5]. Therefore, effectively replacing these batteries demands significant enhancements in PEH efficiency, both in terms of power output and device longevity. Hence, optimized structures via metaheuristic algorithms are proposed to enhance the power/mass ratio [6–8]. Additionally, there is a large amount of research about the topology of the cantilever beam harvester [9, 10]. In recent years, another notable class of structural design is the honeycomb structure inspired by natural patterns. Due to their lightweight and impressive mechanical performance, these structures can enhance the output power of PEHs [11].

The physical characteristics of honeycomb materials are obtained by the geometry of their unit cells. Their elastic modulus, Poisson’s ratio and bulk modulus can be flexibly varied while maintaining the same occupied volume. Due to their advantageous properties, honeycomb structures have been employed in sensors [12], biomedical engineering [13] and energy absorbing applications [14]. The parameters of the unit cells can be tuned to achieve either a positive or negative Poisson’s ratio. Taking advantage of the negative Poisson’s ratio (NPR) of the metamaterial structural substrate referred to the “auxetic structure,” longitudinal and lateral average stresses of the piezoelectric layer bonded to the substrate have the same sign. As a result, the output power of PEHs increases. Recently, auxetic structural properties have been applied in the energy harvesting field to substantially increase the PEH power-to-mass ratio [14, 15]. The auxetic effect is particularly noticeable when using a polyvinylidene fluoride (PVDF) active layer, as its rigidity is typically lower than that of the metamaterial substrate. This amplifies auxetic behaviour, making the PEH more effective than the one with a conventional substrate. In this context, Ichige et al. [16] introduced a PEH that employed a metamaterial as the substrate. Compared with a conventional substrate PEH, their proposed harvester achieved 48% reduction in the resonant frequency, and 3.2 times increase in output power. Ravanbod and Ebrahimi-Nejad [17] proposed an auxetic meta-structure booster (AS-II) which exhibited remarkable improvement in PEH performance. The AS-II harvester achieved 14.61 times more harvested power and 5.8 times less structural weight compared to the solid substrate harvester. Additionally, He et al. [18] proposed cantilever harvesters with honeycomb substrates of positive, zero, and negative Poisson’s ratios using PVDF actives layers. They employed an optimization method combining finite element simulations, neural networks, and an NSGA-II optimization algorithm to identify the best geometric parameters of the substrate of a harvester having a resonance frequency of 75 Hz. Their results showed an increase of about 30% of the produced power compared to the harvester with a conventional substrate. Nevertheless, using a piezoelectric ceramic (PZT) active layer does not offer significant advantages in terms of harvested power, since the auxetic behaviour is not pronounced because of the high stiffness of the active layer compared to the metamaterial substrates [18]. In the literature, it has been reported that honeycomb structures, whether PPR or NPR, do not achieve higher harvested power than solid substrates when the substrate dimensions are kept constant. However, when a tip mass is added, Chandrasekharan and Thompson [19] showed that a bimorph beam with a honeycomb substrate and PZT layers provides a power-to-mass ratio about 25% higher than that of a solid substrate. The advantage of using a PZT layer instead of PVDF is that the PZT offers much higher piezoelectric coupling coefficients and power output, making it the preferred choice for energy harvesters [20]. The literature does not provide a definitive recommendation for combining piezoceramic layers with metamaterial beam substrates, particularly in configurations without a tip mass. Consequently, there is a need for a careful optimization of the metamaterial unit cell geometry to enhance the power-to-mass ratio of the harvester.

This paper proposes the use of a pretrained neural network (NN), trained through backpropagation algorithms, to predict the electromechanical response of PEHs with both PPR and NPR configurations, thereby reducing reliance on computationally expensive FE simulations. The trained metamodel is further integrated with an evolutionary algorithm to optimize the physical attributes of the harvesters for enhanced performance. The NN is trained on a dataset generated from an FE model, which was rigorously validated against a reference model from the literature. In addition, a lifetime analysis of honeycomb-based PPR and NPR structures is carried out using FE simulations combined with Basquin’s law, allowing an approximate estimation of the durability of these harvesters. The results demonstrate that such honeycomb-based PEHs can sustain high accelerations over an extended lifetime, thereby offering a viable pathway to narrow the gap between piezoelectric harvesters and conventional button batteries.

The remainder of this paper is organized as follows. Section 2 describes the design of the harvester. Section 3 presents the finite element model along with its validation. Section 4 discusses the structure and performance of the NN models, while Section 5 details the optimization process. Section 6 presents the results and their interpretation. Finally, Section 7 concludes the paper.

2 Piezoelectric Harvester with Metamaterial Substrate

In this section, the proposed structure of the PEH is presented in detail. Additionally, the different types and dimensions of the microstructures employed in the optimization study are also introduced.

2.1 Harvester Design

The harvester considered in this paper is based on a bimorph configuration with a base excitation. This structure consists of two piezoelectric layers bonded onto a microstructure substrate as presented in Figure 1. The device operates in the $d_{31}$ mode with polarization of the PZT layers oriented in opposite directions meaning that the layers are connected in series to the load resistance R.

Figure 1 Two-dimensional representation of the bimorph energy harvester device.

2.2 Microstructure Geometries

Instead of using a conventional plain substrate, microstructure substrates are used in this research work. Typically, there are two types of microstructures: one with a positive Poisson’s ratio (PPR) and another with a negative Poisson’s ratio (NPR). The cell optimized in this work for the PPR case is the conventional honeycomb structure, as presented in Figure 2(a). For the NPR case, the optimized cell is the re-entrant hexagonal structure, also called the auxetic honeycomb, as shown in Figure 2(b).

Figure 2 Unit cell representation of (a) conventional and (b) auxetic honeycomb.

The PPR honeycomb structure is inspired from natural biological forms [21], while the NPR honeycomb is probably the first auxetic material proposed in the literature for energy harvesting applications [22]. The unit cell geometry is defined by the following parameters: $L_u$ (the unit cell length), $H_u$ (the unit cell height), t (the wall thickness) and $\theta$ (the cell angle). The internal angles of the two honeycomb types have opposite orientations. Furthermore, l and h denote the inclined and vertical lengths, respectively. To estimate the effective equivalent mechanical properties of the microstructure substrate, homogenization techniques can be employed. In this study, the Gibson formulation [23] is used, which is particularly suitable under the assumption that the ratio (t/l) is very small as the effective properties are determined using Bernoulli beam theory. For higher (t/l) ratios, axial and shear deformations must also be considered in addition to bending deformation. The inclined length (l) and the vertical length (h) (shown in Figure 2) can be determined using the following trigonometric relations.

$$l=\frac{L_u}{4\cos \theta }$$	(1)
$$h=\frac{H_u}{2}-l\sin \theta$$	(2)

The fabrication of metallic honeycomb structures, particularly those exhibiting NPR and PPR behaviours, has become increasingly feasible due to advancements in manufacturing techniques. Additive manufacturing (AM) methods such as selective laser melting and electron beam melting now enable the direct fabrication of complex honeycomb structures in metals like stainless steel and titanium, with high precision [24]. In parallel, traditional manufacturing methods using techniques such as Computer Numerical Control (CNC) machining or micro-end milling remain useful and widely used to produce brass and aluminium honeycomb cores. For example, Qiu et al. [25] investigated the machining characteristics of aluminium honeycomb cores, focusing on cutting forces, entrance angle, and cutting speed. By considering these parameters, honeycomb structure with high precision in the millimetre scale can be fabricated.

3 Finite Element Formulation for a Bimorph Piezoelectric Energy Harvester

3.1 Local Equations of the Coupled System

In this study, PEH is modelled as a structure composed of an elastic domain, denoted by ${\mathrm{\Omega }}_E$, and one or more piezoelectric domains, denoted by ${\mathrm{\Omega }}^{(p)}$, where $p\in \{1,2\}$ refers to the number of piezoelectric patches bonded to the elastic substrate. The set (${\mathrm{\Omega }}_E,{\mathrm{\Omega }}^{(1)},{\mathrm{\Omega }}^{(2)}$) constitutes a partition of the overall domain ${\mathrm{\Omega }}_s$. The structure is subjected to a prescribed displacement $w_i^d$ on a section ${\mathrm{\Gamma }}_w$ and a surface force density $F_i^d$ on the complementary boundary section ${\mathrm{\Gamma }}_T$.

The electric boundary conditions consist of a prescribed electric potential ${\psi }^d$ on ${\mathrm{\Gamma }}_{\psi }$ boundaries and a surface charge density $Q^d$ applied on the remaining boundaries ${\mathrm{\Gamma }}_q$. The components of the linearized strain tensor are denoted by ${\epsilon }_{ij}$ and the corresponding stress tensor components are represented by $T_{ij}$. Additionally, $D_i$ denotes the $i$th component of the electric displacement vector and $E_i$ represents the components of the electric field vector. The parameter ${\rho }_s$ indicates the mass density of the structure and $n_i$ is the unit outward normal vector to the boundary of the domain ${\mathrm{\Omega }}_s$. The local governing equations for the coupled piezoelectric–structure interaction problem, as illustrated in Figure 1, are presented along with the appropriate boundary conditions [6].

$$T_{ij,j}={\rho }_s\frac{{\partial }^2{w}_i}{\partial t^2}\mathrm{in}{\mathrm{\Omega }}_s$$	(3a)
$$T_{ij}{n}_j=F_i^d\mathrm{on}{\mathrm{\Gamma }}_T$$	(3b)
$$w_i=w_i^d\mathrm{on}{\mathrm{\Gamma }}_w$$	(3c)

$$D_{i,i}=0\mathrm{in}{\mathrm{\Omega }}_s$$	(4a)
$$D_i{n}_i=-Q^d\mathrm{on}{\mathrm{\Gamma }}_q$$	(4b)
$$\psi ={\psi }^d\mathrm{on}{\mathrm{\Gamma }}_{\psi }$$	(4c)

Equation (3a) represents the elastodynamic equation in the case of negligeable body forces. Equations (3b) and (3c) define the mechanical boundary conditions. Equation (4a) corresponds to the electric charge equation for a dielectric medium, while Equations (4b) and (4c) specify the electric boundary conditions.

In the piezoelectric material, the stress tensor $T_{ij}$ and the electric displacement $D_i$ are linearly related to the strain tensor ${\epsilon }_{kl}$ and the electric field $E_k$ via the piezoelectric constitutive equations [26, 27].

$$T_{ij}=c_{ijkl}{\epsilon }_{kl}-e_{kij}{E}_K$$	(5)
$$D_i=e_{ikl}{\epsilon }_{kl}+{\in }_{ik}{E}_k$$	(6)

where $c_{ijkl},e_{kij}$ and ${\in }_{ik}$ represent respectively the elastic modulus at constant electrical field, the piezoelectric constant coefficient and dielectric permittivity at constant strain. Additionally, the following two gradient relationships which relate the linearized strain tensor ${\epsilon }_{kl}$ to the displacement $w_k$ and the electric field $E_k$ to the electric potential $\psi$ are used [6]:

$${\epsilon }_{kl}=\frac{1}{2}(w_{k,l}+w_{l,k})$$	(7)
$$E_k=-{\psi }_{,k}=-\frac{V^{(p)}}{h_p}$$	(8)

where $V^{(p)}$ is defined as the applied potential difference (voltage) between the electrode surfaces of the piezoelectric element and $h_p$ is the thickness of the piezoelectric layer.

3.2 Finite Element Formulation

The variational formulation in terms of $(w_i,V^{(p)})$ presented in [28] is adopted to derive the finite element formulation for elastic structures incorporating piezoelectric elements. The resulting finite element equations governing the damped electromechanical response under harmonic excitation at an angular frequency $\omega$ in steady-state conditions are given by:

$$-{\omega }^2{\text{M}}_m\text{U}+j\omega {\text{C}}_m\text{U}+{\text{K}}_m\text{U}+{\text{K}}_c\text{V}=\text{F}$$	(9)
$${\text{K}}_c^T\text{U}-{\text{K}}_e\text{V}=\text{Q}$$	(10)

Here, ${\mathit{M}}_m$ is the $(N\times N)$ global mass matrix, ${\mathit{C}}_m$ is the $(N\times N)$ global damping calculated using Rayleigh damping, ${\mathit{K}}_m$ is the global $(N\times N)$ stiffness matrix, ${\mathit{K}}_c$ is the global $(N\times P)$ electromechanical coupling matrix and ${\text{K}}_e(P\times P)$ is the diagonal global capacitance matrix. Additionally, the vector $\mathit{F}(t)=\mathit{F}{e}^{j\omega t}$ represents the global $(N\times 1)$ vector of mechanical forces which can be expressed as $\mathit{F}(t)=-{\omega }^2{\mathit{m}}^{\ast }{w}_b{e}^{j\omega t}$, where ${\mathit{m}}^{\ast }$ is the effective mass $(N\times 1)$ vector and $w_b$ is the base transversal displacement. $\mathit{Q}(t)=\mathit{Q}{e}^{j\omega t}$ is the global $(P\times 1)$ vector of electrical charge output, $\mathit{U}(t)=\mathit{U}{e}^{j\omega t}$ is the global $(N\times 1)$ vector of mechanical coordinates $w_i$ and $\mathit{V}(t)=\mathit{V}{e}^{j\omega t}$ is the global $(P\times 1)$ vector of output voltages across the piezoelectric elements. $N$ and $P$ are respectively the number of mechanical degrees of freedom and the number of piezoelectric elements. Since the piezoelectric patches have thin, highly conductive electrodes on their top and bottom surfaces, all finite elements within a patch generate the same voltage, so that each patch is represented by a single electrical degree of freedom.

Equation (9) represents the mechanical equation of motion, incorporating electrical coupling and a forcing vector $\mathit{F}$. In contrast, Equation (10) represents the electrical circuit equation with a mechanical coupling term. In this study, the harvested energy is dissipated through a resistive load $R$. Applying Ohm’s law, the following equation relates the voltage vector $\mathit{V}(t)$ to the charge vector $\mathit{Q}(t)$:

$$\mathit{V}(t)=R\dot{\mathit{Q}}(t)$$

(11)

Consequently, Equation (10) can be rewritten using Equation (11) under harmonic excitation and steady state conditions as follows:

$$\text{U}=\text{V}\left(\frac{\frac{1}{R}+j\omega {\text{K}}_e}{j\omega {\text{K}}_c^T}\right)$$

(12)

By combining Equations (9) and (12), a finite element matrix equation can be derived to calculate the voltage across the load resistance. The resulting voltage Frequency Response Function (FRF) of the harvester is then obtained as follows:

${\displaystyle \frac{\mathit{V}(t)}{-{\omega }^2{w}_b{e}^{j\omega t}}}$	$={\displaystyle \frac{j\omega {\mathit{K}}_C^T}{\left(\frac{1}{R}+j\omega {\mathit{K}}_e\right)}}$
	$\times {\left(-{\omega }^2{\mathit{M}}_m+j\omega {\mathit{C}}_m+{\mathit{K}}_m+{\displaystyle \frac{j\omega {\mathit{K}}_c{\mathit{K}}_c^T}{\left(\frac{1}{R}+j\omega {\mathit{K}}_e\right)}}\right)}^{-1}{\mathit{m}}^{\ast }$	(13)

The power FRFs can be derived from the voltage FRFs given by Equation (3.2) as follows:

${\displaystyle \frac{\mathit{P}(t)}{{(-{\omega }^2{w}_b{e}^{j\omega t})}^2}}$	$={\displaystyle \frac{1}{R}}({\displaystyle \frac{j\omega {\mathit{K}}_C^T}{\left(\frac{1}{R}+j\omega {\mathit{K}}_e\right)}}(-{\omega }^2{\mathit{M}}_m+j\omega {\mathit{C}}_m$
	${+{\mathit{K}}_m+{\displaystyle \frac{j\omega {\mathit{K}}_c{\mathit{K}}_c^T}{\left(\frac{1}{R}+j\omega {\mathit{K}}_e\right)}})}^{-1}{\mathit{m}}^{\ast }){}^2$	(14)

3.3 Validation of the FE Model

The validation data used in this section are derived from the paper of Chandrasekharan and Thompson [19]. The harvester considered in [19] consists of a cantilever beam with two PZT layers bonded to a brass honeycomb substrate. The global geometric parameters of the harvester are as follows: the overall length $L=66.62$ mm and the width $b=13.68$ mm. The thickness of the honeycomb substrate $(h_s)$ is 0.76 mm, while the thickness of the active layer, made of PZT-5H ($h_p$) is 0.26 mm.

The honeycomb cell parameters are as follows: $\theta ={30}^{\circ }$, $L_u=5.33$ mm, $H_u=4.56$ mm and $t=0.152$ mm. The substrate is modelled as a homogenized layer, and all the equivalent properties integrated into the finite element model are determined using the Gibson formulation. The complete effective elastic properties and the relative density of the brass honeycomb structure are presented in Table 1 [19].

Table 1 Effective elastic constants and density of hexagonal honeycomb structures [19]

Mechanical Properties	Honeycomb
$E_{S,1}^{\ast }$(GPa)	0.22
$E_{S,2}^{\ast }$(GPa)	0.22
$E_{S,3}^{\ast }$(GPa)	11.05
$G_{S,12}^{\ast }$(GPa)	0.05
$G_{S,13}^{\ast }$(GPa)	2.06
$G_{S,23}^{\ast }$(GPa)	2.06
${\nu }_{S,12}^{\ast }={\nu }_{S,21}^{\ast }$	1
${\nu }_{S,31}^{\ast }={\nu }_{S,32}^{\ast }$	0.34
${\nu }_{S,13}^{\ast }$	0.0008
${\nu }_{S,23}^{\ast }$	0.0008
${\rho }^{+}={\rho }_E^{\ast }/{\rho }_E$	0.114

$E_{S,1}^{\ast }$,$E_{S,2}^{\ast }$, $E_{S,2}^{\ast }$ are the elastic moduli in the $x_1$, $x_2$, and $x_3$ directions. Additionally, the in-plane shear modulus is represented by$G_{S,12}^{\ast }$and the out-of-plane shear modulus denoted by $G_{S,13}^{\ast }$ and $G_{S,23}^{\ast }$. Moreover, ${\rho }_E^{\ast }$ and ${\rho }_E$ are the densities of the cellular structure and solid material, respectively. The Poisson’s ratios ${\nu }_{S,ij}^{\ast }$ define the strain relationships between the orthogonal directions $i$ and $j$. The PZT-5H material has a density ${\rho }_p=7200$ kg/m³, a permittivity of 25.55e^-9 F/m, a piezoelectric strain constant of 320 pm/V, a Young’s modulus $E_P=62$ GPa and a Poisson’s ratio of 0.31.

In this work, proportional damping with ${\zeta }_1=0.019$ and ${\zeta }_2=0.033$ are applied to the electromechanical system to ensure consistency with [19, 29]. Additionally, the harvester is subjected to a harmonic base excitation with an amplitude of 1 g. Based on the equations presented in Section 3, a three-dimensional finite element model is built using the parameters listed above. A frequency-domain study is performed by sweeping the frequency from 100 to 260 Hz in 0.25 Hz increments. The model employs a tetrahedral mesh consisting of 6438 elements.

Figure 3 shows a comparison between the power output at the optimal load resistance of 20 k$\mathrm{\Omega }$ obtained from the current FE model and the results reported by Chandrasekharan and Thompson [19]. The good agreement between the two confirms the validity of the present model, with a maximum 3.45% discrepancy in power and 1.1% discrepancy in frequency.

Figure 3 Comparison of FRF power curves between the curves predicted by the model used in this paper and that of Chandrasekharan and Thompson [19].

4 Metamodel Based on Artificial Neural Networks

The output power of a PEH that integrates a microstructure substrate with PZT layers is generally lower than that of a harvester using a solid continuous beam for the same dimensions of the harvester. However, one of the main objectives of incorporating microstructures into the harvester is to reduce its weight, which is a critical constraint in many aeronautic and aerospace applications. As mentioned earlier, the research work presented in this paper aims to optimize the geometry of the unit cells presented in Figure 2. The focus is solely on the first vibration mode, as it typically yields the highest output power. Moreover, the first natural frequency can generally be tuned to match specific target frequencies in engineering applications, which are often low (typically below 200 Hz) in real-world environments. The initial geometry of the cells is presented in reference [19] which details the validation step for their presented results. In this section, the neural network model used to predict the performance of PEHs is presented and validated. A metamodel is employed in this study to reduce the computational cost of evolutionary algorithms in solving the optimization problem. To reduce this computational demand, an artificial neural network (ANN) is employed as an efficient alternative, significantly decreasing the computation time. Furthermore, the ANN model is trained on three-dimensional finite element simulations to ensure higher accuracy compared to analytical models, especially since strong anisotropy effects arise when the honeycomb structure deviates from conventional designs. Additionally, the nondominated sorting genetic algorithm II (NSGA-II) is then used to optimize the geometric parameters of the NRP and PPR unit cells. NSGA-II is particularly suitable because it is a derivative-free, multi-objective optimization algorithm, capable of handling complex, non-differentiable objective functions such as power output and structural mass. ANN evaluates the performance of each candidate design in the NSGA-II population, enabling rapid identification of the Pareto-optimal front. This coupling of ANN and NSGA-II allows for an effective exploration of the design space, balancing competing objectives such as maximizing power output while minimizing weight, and provides a robust framework for optimizing microstructured PEHs in practical engineering applications.

4.1 Artificial Neural Network Architecture for PEH Metamodeling

A metamodel is a simplified representation of a complex simulation, designed to save time and computational resources during tasks such as optimization and sensitivity analysis. Traditional techniques, such as Response Surface Models are commonly used in this context but often struggle when dealing with highly nonlinear systems or interactions between multiple variables. To overcome these limitations, more advanced modelling approaches like ANN have been adopted. Inspired by the structure and functioning of the human brain, ANNs are powerful tools for modelling complex, nonlinear relationships and are widely used for prediction tasks.

In this paper, a Multilayer Perceptron (MLP) is employed as a metamodel to predict the performance of PEH. Two separate MLPs were trained using the backpropagation algorithm: one for harvesters with PPR honeycomb and another for those with NPR honeycomb. The backpropagation algorithm adjusts the weights of the connections between neurons by propagating the error gradient backward through the network. This training process aims to minimize the difference between the predicted outputs and the actual target values, allowing the network to accurately learn the underlying patterns within the data.

To evaluate the effectiveness of the harvester, especially when using a microstructure substrate, it is essential to consider both the peak power at the first resonance frequency and the total mass of the device. Consequently, the ANN architectures developed for the two cases (PPR and PPR cases) each has four input parameters corresponding to the unit cell geometry, and two output parameters: the peak power and the total mass of the harvester. The design of the ANN is illustrated in Figure 4.

Figure 4 Structure of the ANN model.

Choosing the number of hidden layers and the number of nodes per layer is a critical step in developing effective neural network models. Detailed design guidelines can be found in [30]. Although increasing the number of layers and nodes may enhance the model’s capacity, it can also reduce its generalization ability and increase the risk of overfitting, which is often undesirable. In this study, a balanced architecture is selected, consisting of two hidden layer with 16 neurons. The ReLU activation function is adopted to effectively perform the regression tasks in both cases (PPR and NPR). The key parameters of the two constructed neural networks are summarized in Table 2.

Table 2 MLP parameters

Parameter	Value
Number of nodes in the input layer	4
Number of hidden layers	2
Number of nodes in each hidden layer	16
Learning rate	0.01
Number of nodes in the output layer	2
Size of the total data set (for each case)	400

A total of 400 random samples were generated within the specified ranges defined in the Equations (15)–(16) given below for each design variable with responses calculated using the previously described finite element model. Among these, 288 simulations were used for training, 72 for validation, and 40 for testing. This data split is used for both designs. Each simulation is performed using a 22 k$\mathrm{\Omega }$ load resistance, which is selected as the optimal value for all cases, since the resonance frequencies of the simulated systems are all within the vicinity of 179 Hz.

The optimization problem seeks to determine the optimal values for the geometric parameters ($L_u$, $H_u$, $t$) and the cell angle ($\theta$). These parameters are subject to the following bounds, which are chosen to ensure realistic and physically feasible designs, as supported by [16]:

$$\begin{array}{c}3.5<L_u<7\mathrm{mm}\hfill \\ 3<H_u<5\mathrm{mm}\hfill \\ 0.1<t<0.25\mathrm{mm}\hfill \end{array}$$

(15)

The cell angle ($\theta$) is subject to specific constraints that define the two distinct optimization scenarios for: (PPR)

$${20}^{\circ }<\theta <50$$

(NPR)

$$-{50}^{\circ }<\theta <-{20}^{\circ }$$

(16)

Optimization aims to maximize energy harvesting performance while minimizing structural mass. This multi-objective problem is formally stated as minimizing the weighted vector of the negative harvested power and the structural mass:

Objective: Min (-Power, Mass)

(17)

The training performance of the neural networks is quantified using the coefficient of determination $(R^2)$, a standard metric in regression analysis that measures the proportion of variance in the target variable explained by the model predictions.

$$R^2=1-\frac{{\sum }_{i=1}^n{(y_i-{\widehat{y}}_i)}^2}{{\sum }_{i=1}^n{(y_i-\overline{y})}^2}$$

(18)

Here, $y_i$ represents the observed values, ${\widehat{y}}_i$ denotes the predicted values, $\overline{y}$ is the average of the observed values, and $n$ is the total number of data points. The $R^2$ value varies between 0 and 1, where 1 signifies a perfect fit to the data.

In addition to $R^2$, we incorporated the Mean Squared Error (MSE) into the loss function to further evaluate the model’s performance. MSE calculates the average of the squared differences between the predicted values, ${\widehat{y}}_i$ and their corresponding observed values $y_i$ [31]:

$$\mathit{MSE}=\frac{1}{n}\sum _{i=1}^n{(y_i-\widehat{y_i})}^2$$

(19)

Since ANN predicts multiple outputs with different units (peak power in watts and mass in kilograms), the combined MSE reported here corresponds to the average of normalized squared errors across outputs and is therefore dimensionless. MSE will be used to evaluate the accuracy of the prediction. The outputs are normalized individually, typically by subtracting the mean and dividing by the standard deviation of each output.

4.2 Neural Networks Training Accuracy

After training the ANN model, it is necessary to evaluate the performance of the model and its capability to mimic the behaviour of the harvester in terms of peak power and mass. Figure 5 shows the variation of the MSE for both training and validation sets for the PPR and NPR configurations. Figure 5(a), which represents the PPR configuration, presents excellent convergence behaviour and high prediction accuracy, with the MSE loss reaching 0.0012 for the training set and 0.0019 for the validation set, both at 300 epochs. Figure 5(b), which represents the NPR configuration, shows MSE losses of 0.0020 for the training set and 0.0025 for the validation set at 300 epochs.

Figure 5 Mean squared error for training and validation data vs. the number of epochs: (a) PPR and (b) NPR.

Figure 6 shows the regression diagrams using the test set. Figure 6(a) represents the power predictions of the PPR honeycomb, while Figure 6(b) represents the power predictions of the PPR honeycomb harvester. In both cases, the regression correlation coefficient, $R^2$, as defined in Equation (18), is above 0.98. In addition, Figures 6(c) and 6(d) show the corresponding error distributions. The prediction errors remain very small and centred around zero, indicating that the ANN does not exhibit noticeable bias and maintains stable performance across the dataset. Overall, these results demonstrate that ANN provides highly reliable predictions for both power and mass.

Figure 6 Regression fit and error distributions for power and mass of the PRP honeycomb harvester.

Figure 7 shows the predictive performance of the regression models for both power and mass, along with the corresponding error distributions for the NPR honeycomb harvester. Figures 7(a) and 7(b) compare the predicted values to the true data for power and mass, respectively. In both cases, the points align closely with the identity line, demonstrating excellent agreement and very high correlation coefficients $(R>0.98)$. Figures 7(c) and 7(d) show the histograms of the prediction errors. The power error distribution is centred near zero with a few isolated outliers, while the mass error distribution is even more concentrated, indicating highly stable and unbiased predictions. Overall, the figure highlights the accuracy and robustness of the developed models.

Figure 7 Regression fit and error distributions for power and mass of the PRP honeycomb harvester.

The results presented in Figures 6 and 7 confirm the accuracy of the ANN model. The next section will present the integration of the trained models with an evolutionary algorithm to optimize the dimensional parameters of the PEHs.

5 Optimization of the Unit Cell of the Microstructure Substrate

To mimic the behaviour of the harvester using ANN can be highly effective, especially when applied to an optimization problem, because this approach significantly reduces computational time. For instance, performing full FE simulations over the excitation range from 100 Hz to 260 Hz with a 0.25 Hz increment requires approximately 22 minutes using an Intel^® Core i7-8665 processor and 8 GB of RAM, whereas ANN can predict the same responses almost instantaneously. This demonstrates the value of using a metamodel, which enables faster evaluations and more efficient optimization. To further enhance the optimization process, the NSGA-II algorithm is employed, effectively handling multi-objective optimization tasks.

NSGA-II is a multi-objective optimization algorithm developed to solve problems involving multiple conflicting objectives [32]. It belongs to the family of evolutionary algorithms that mimic natural selection and genetic processes to efficiently explore and identify multiple optimal solutions for a given problem. Figure 8 presents the flowchart of the NSGA-II algorithm. For a detailed explanation of the algorithm, the readers may want to consult [33].

Figure 8 Optimization process for NSGA-II.

The trained backpropagation neural network models are integrated into the fitness function of the NSGA-II algorithm to optimize key performance attributes of the substrate cells of the PEHs, as defined in Equations (15) and (16), with the objective of the optimization to maximize the peak power and to minimize the overall mass of the harvester (Equation (17)). The upper and lower limits of the harvester variables are set based on fabrication technology constraints. The challenge of the fabrication is to make a thin thickness; however, values of the thickness close to the lower limit can be fabricated according to Hexcel Corporation [34]. The optimization process is carried out using a population size of 100 and 500 generations. The resulting optimized designs and a comparative analysis of various PEH configurations will be presented in the following sections.

6 Results and Discussion

The convergence of the NSGA-II algorithm was evaluated using the hypervolume (HV) indicator, computed at each generation. HV measures the portion of the objective space dominated by the evolving Pareto front relative to a reference point set 10% worse than the worst observed objectives. As it does not require knowledge of the true Pareto-optimal front, it reliably reflects both convergence and diversity. In this study, the hypervolume rises rapidly during the initial generations and then stabilizes, as shown in Figure 9(a) for the PRP harvester and Figure 9(b) for the NPR harvester, indicating that the algorithm has reached a steady-state Pareto front.

Figure 9 Convergence analysis for (a) PPR and (b) NPR.

These figures show that the hypervolume grows quickly in the first 10 generations, meaning the initial iterations improve the solutions most. After that, growth slows and levels off around generations 20–25, with small fluctuations between 760–780 for NSGA-II PRP optimization and 1850–1860 for the NPR harvester problem. This shows that the algorithm converges well, with most improvements happening in the first half of the run. A sensitivity analysis of the population size, mutation rate, and crossover probability indicates that a population size of 100 provides a good balance between solution quality and computational time, while a mutation rate of 0.15 and a crossover probability of 0.7 yield the best overall performance.

Building on this convergence analysis, Figures 10(a) and 10(b) show the Pareto fronts of the optimization results for the initial geometry presented previously in Section 2.3 for the PPR and NPR cells.

Figure 10 Pareto front for (a) PPR and (b) NPR.

These plots illustrate the trade-off between peak power output and mass of the harvester. Each point on the fronts represents a non-dominated solution, highlighting the inherent conflict between maximizing power and minimizing mass. As expected, achieving higher peak power typically requires increased mass, while reducing the harvester mass leads to a reduced power.

The shape and spread of the front for PPR in Figure 10(a) shows a relatively linear trend, suggesting a steady trade-off and a predictable design space. In contrast, the NPR front in Figure 10(b) spans a wider range of values, indicating greater flexibility in performance but also a steeper performance-to-mass compromise. This highlights the fundamental conflict between the two objectives: improving output power generally requires heavier or more complex structures, limiting the ability to reduce simultaneously the mass of the harvester.

The distribution and continuity of the Pareto fronts validate the quality of the optimization process, reflecting both convergence and diversity of the non-dominated solutions. The selection of a final design point depends on application-specific constraints: mass-sensitive systems, such as lightweight electronics [1], may favour low-mass solutions (e.g., PPR-optimized designs), while performance-critical applications, such as aerospace systems [2], might prefer configurations that deliver higher absolute power, even at the cost of increased mass (e.g., NPR-optimized solutions). This trade-off is particularly relevant in aerospace, where power output is often prioritized over minimal mass to satisfy the durability and robustness of the harvester [3].

In this study, priority is given to maximizing the output power per unit mass of the PEHs. Consequently, the optimal solutions, which provide the highest power-to-mass ratio for the PPR PEH and NPR PEH are (1951.51 $\mu$W, 5.03 g) and (2498.34 $\mu$W, 6.01 g), respectively, as illustrated in Figure 10 (corresponding to the top-right point). The optimal design for each microstructured substrate is shown in Figure 11.

Figure 11 Optimized designs associated with the optimal solutions: (a) PPR and (b) NPR.

Figure 12 illustrates the relationship between each geometric design parameter of the microstructure’s unit cell ($L_U$, $H_U$, $\theta$, and $t$) and the output power based on the optimized solutions from the Pareto front. Figure 12(a) corresponds to the PPR harvester, while Figure 12(b) shows the results for the NPR harvester.

The parameters $L_U$ and $H_U$ (the horizontal and vertical dimensions of the unit cell) show a negative correlation with power output. This indicates that smaller cell sizes enhance deformation transmission to the piezoelectric layer, thus increasing energy conversion efficiency. However, it is important to note that the optimized cell sizes are relatively small but not the smallest possible to fabricate. Moreover, increasing the absolute value of the angle ($\theta$) within the range of 40^∘ to 50^∘ appears to enhance the mechanical behaviour of the structure, thereby promoting higher energy harvesting performance. Similarly, the thickness parameter (t) also correlates positively with power output. A greater wall thickness improves the stiffness of the structure and enhances stress transfer, ultimately boosting power generation.

Figure 12 Evolution of the Pareto front as a function of the power output: (a) PPR and (b) NPR.

Figure 13 compares the FEM simulation results of the power output per acceleration (g $=$ 9.81 m/s²) in the vicinity of the first resonant frequency for each optimized PEH, the non-PPR design from [19], and the continuous solid brass model from [19], using the optimal electrical load resistance specific to each configuration. The non-PPR PEH delivers 1395.5 $\mu$W at a mass of 4.1 g, while the continuous solid brass PEH achieves 4426.83 $\mu$W at 9.5 g.

Figure 13 Harvested power responses of PEHs for different harvester substrate designs.

These results demonstrate that the optimized designs perform as intended. We introduce a significant performance metric, power-per-mass ratio ($P_m$), to better evaluate and compare the efficiency of each configuration.The optimized honeycomb and optimized auxetic configurations show an increase in $P_m$ of 13.8% and 22%, respectively, compared to the unoptimized honeycomb design. However, the solid continuous brass configuration still demonstrates the highest performance, with $P_m$ value 10% higher than that of the auxetic configuration and 19% higher than that of the optimized honeycomb configuration.

To enable a fair and meaningful comparison between the different energy harvester configurations, it is essential that their first natural frequencies are tuned to match the same input excitation frequency. Indeed, any deviation from the resonant frequency can lead to a substantial reduction in power output due to the high sensitivity of piezoelectric harvesters to resonance conditions. In this study, the input driving frequency is assumed to be approximately 160 Hz. To match this frequency with the first bending mode of each configuration, the length of each configuration is adjusted accordingly. Table 3 shows the power output of each configuration based on the validated FE model.

Table 3 highlights that the solid configuration delivers the highest peak power (3548.52 $\mu$W), but also exhibits the largest mass (9.1 g), resulting in a power-per-mass ratio of 393.8 $\mu$W/g. In contrast, the optimized auxetic harvester achieves a lower peak power (2357.21 $\mu$W) but with a significantly lower mass (5.7 g), yielding the highest power-per-mass value of 413 $\mu$W/g, indicating superior efficiency in terms of power output per unit mass. The optimized honeycomb configuration performs comparably to the solid substrate harvester, but with a slightly higher power-per-mass ratio (395.52 $\mu$W/g). The unoptimized honeycomb configuration still performs significantly worse than the others, highlighting the importance of structural optimization.

In this part of the paper, based on finite element method (FEM) simulations, we comprehensively examine the performance of the proposed configurations in terms of the durability and lifespan of the active layer of the harvester. According to the constitutive relations of the piezoelectric effect (strain–charge form) provided in Equations (5) and (6). These coupled equations indicate that the efficiency of the PEHs is directly related to the ability of the harvester to transmit mechanical stress and converts it into strain.

Table 3 Quantitative comparison of the performance of the four piezoelectric harvesters with equal resonant frequency of 160 Hz

Harvester	Mass (g)	Peak Power ($\mu$W)	$P_m$ ($\mu$W/g)	Length (mm)
Solid	9.01	3548.52	393.8	62.5
Unoptimized honeycomb	4.4	1588.38	360.1	71
Optimized honeycomb	5.03	1989.51	395.52	66.6
Optimized auxetic	5.7	2357.21	413	64.9

Figure 14 illustrates the distribution of the maximum von Mises stress in the piezoelectric layer of the solid, the optimized PPR, and the optimized NPR energy harvesters, evaluated at their first resonance frequency of 160 Hz. The results correspond to simulations conducted under a 1 g acceleration applied along the z-axis. The highest stress concentrations (indicated in red in Figure 14) are consistently located near the fixed boundary in all configurations. While the maximum von Mises stress values in the piezoelectric layers are comparable for the optimized PPR and NPR designs. The harvester with the solid substrate exhibits a noticeably higher stress level. The peak stress values are equal to 7.62 MPa for the solid harvester, 3.95 MPa for the PPR harvester, and 3.71 MPa for the NPR harvester. All of these stress values are well below the tensile strength of the piezoelectric material (30 MPa).

We conducted a simplified simulation to evaluate the influence of the amplitude of the base excitation acceleration (g) on the performance of the energy harvesters for each configuration. An increase in acceleration generally leads to higher mechanical stress in the active piezoelectric layer, which can enhance power output. However, this improvement is constrained by the fatigue limits of the piezoelectric material, as excessive stress can significantly shorten its operational lifespan. In this study, we aim to assess the trade-off between power generation and durability by analysing the piezoelectric layer’s life cycle under increasing acceleration levels. The objective is to identify an optimal acceleration value that maximizes power output while ensuring a reasonable life expectancy of the piezoelectric layer.

Figure 14 von Mises stress distributions in the piezoelectric patch for (a) solid, (b) NPR and (c) PPR harvesters.

Figure 15 presents the evolution of the maximum von Mises stress and power output as a function of the amplitude of the base acceleration for three different brass substrate types: solid, honeycomb, and auxetic. As expected, the power output increases nonlinearly with the amplitude of the acceleration. This trend is approximately quadratic ($\mathit{Power}\propto g^2)$. This is due to the quadratic dependence of the power on the induced voltage generated by the piezoelectric material.

Among the three configurations, the solid substrate generates the highest power output, especially at higher accelerations. This is attributed to its high stiffness and efficient mechanical energy transmission to the piezoelectric layer, and this can lead to larger strain, which induces high voltage. However, this also leads to significantly higher mechanical stress, which may affect the long-term reliability of the device. In contrast, the honeycomb and auxetic substrates generate lower power output but exhibit reduced stress levels. This trade-off seems advantageous, especially in high-acceleration applications.

Figure 15 Maximum power and maximum von Mises stress vs acceleration amplitude at resonant frequency.

In many practical applications, the ambient acceleration levels can be relatively high. For instance, in industrial environments and certain biomedical applications, accelerations can reach up to 10 g or even higher [35]. While higher acceleration levels may increase the energy harvested, they can also impose significant mechanical stress on the active piezoelectric layer, potentially causing damage and thereby reducing its operational lifespan. This trade-off between power output and durability is critical in the design of energy harvesters.

Figure 16 Maximum stress vs. estimated lifetime for different substrates at 160 Hz.

Figure 16 illustrates the relationship between the cycle life of the PZT 5H material and the equivalent stress applied to it. These results were derived using Basquin’s law, in conjunction with the S–N (stress–number of cycles) curves presented in [17]. The analysis compares three substrate configurations: solid, honeycomb, and auxetic under harmonic excitation at 160 Hz.

For an acceleration equal to 3 g, the corresponding electrical power outputs are 30.6 mW, 16 mW, and 12.8 mW for the solid, honeycomb, and auxetic substrates, respectively. These outputs are associated with mechanical stress levels of 22 MPa (solid), 11.3 MPa (honeycomb), and 11.9 MPa (auxetic). According to the fatigue model, a stress level of 22 MPa leads to a predicted lifetime of approximately 0.28 years for the PZT 5H when mounted on a solid brass substrate. In contrast, the honeycomb and auxetic designs result in significantly lower stress amplitudes, yielding estimated lifetimes exceeding 50 years, which can be considered functionally infinite. It should be noted that these estimates do not account for scatter, fatigue crack growth, or other operational uncertainties and should therefore be interpreted as indicative upper bounds rather than guaranteed lifetimes.

These findings demonstrate that while the solid substrate provides the highest power output, it drastically compromises the operational durability of the piezoelectric element. The microstructures substrates (honeycomb and auxetic) offer a more balanced trade-off between energy harvesting performance and fatigue life, making them more suitable for long-term operation in vibratory environments especially under high acceleration amplitudes. However, under low acceleration levels, where the induced stress remains below critical fatigue thresholds, the solid substrate may also achieve a long operational lifetime while maintaining superior power output.

In addition to the analysis of the active layer, the substrate layer must withstand all the induced stress to ensure the long-term functionality of the harvester. Figure 17 presents the von Mises stress distribution for the optimal design under an acceleration of 4 g and a resonant frequency of 160 Hz for: (a) PPR configuration and (b) NPR configuration. The maximum stress reaches 13.5 MPa in the PPR configuration and 14.7 MPa in the NPR configuration.

Figure 17 Stress distribution of the microstructure substrate: (a) PPR and (b) NPR.

To ensure that the honeycomb substrate can withstand the induced stress over time and maintain its structural integrity, its fatigue performance must be carefully evaluated. This typically involves the use of an S–N curve, which must be plotted by considering both the optimized cell geometry and the material properties. In the literature, investigations about this topic may be found. For instance, Palomba et al. [36] investigated the collapse mechanisms of aluminium honeycomb sandwich panels under fatigue bending through three-point bending tests. They experimentally plotted the S–N curves for commercial honeycomb panels with an overall thickness of approximately 11 mm. Their results showed that the endurance limits ranged between 39 MPa and 43 MPa. Similarly, Bianchi et al. [37] presented in plane shear fatigue behaviour for aluminium honeycomb cores and provided corresponding S–N curves for different loading directions. Nevertheless, these studies cannot be directly extrapolated to the current case due to the differences in material and geometry. Knowledge of the fatigue behaviour of the microstructures considered in this paper is relatively poor and should be the subject of further investigations. Furthermore, Neèemer et al. [38] showed that the decrease of the material density leads to a decrease of the endurance limits. As a result, the relative density is a critical parameter for fatigue performance and can be calculated as follows:

$${\rho }_E^{\ast }=\frac{\left(\frac{h}{l}+2\right)\frac{t}{l}}{2\cos \theta \left(\frac{h}{l}+\sin \theta \right)}$$

(20)

where $l$ and $h$ are calculated using Equations (1) and (2). Based on Equation (20) the proposed harvesters with cellular substrates densities have densities of 2561.1 kg/m³ and 4082.1 kg/m³, respectively, for the PPR and NPR substrates.

Brass has a density of approximately 8700 kg/m³ and an endurance limit of about 150 MPa. The effective endurance limit of the cellular structure can be estimated by scaling this limit according to the relative density [39, 40]:

$${\sigma }_e^{\ast }\approx {\sigma }_{e,s}{\left(\frac{{\rho }_E^{\ast }}{{\rho }_E}\right)}^n$$

(21)

where ${\rho }_E^{\ast }$ and ${\rho }_E$ are the densities of the cellular structure and solid material, respectively. The parameter $n$ is an empirical exponent that depends on the geometry and deformation mechanism of the cellular architecture. According to Gibson et al. [23, 40, 41], the mechanical properties of cellular materials scale with $n=2$ for closed-cell foams and $n=1.5$ for open-cell foams. In our case, the PRP and NPR microstructures are similar to open-cell foams; therefore, choosing $n=1.5$ provides a reasonable first-order estimate of the endurance limit. Additionally, [39] suggested $n=1.2$ for tensile loading. Since tensile stresses are generally higher than bending stresses for the same material, selecting $n=1.5$ for bending appears reasonable and consistent with a first-order estimation.

By substituting the known values, the endurance limits of the PPR and NPR structures are estimated to be approximately 24 MPa and 48 MPa, respectively. These results highlight the significant influence of cellular density on fatigue performance: the higher density of the NPR substrate increases its effective endurance limit, leading to a potentially longer fatigue life.

Basquin’s law [42] is applied here to predict the high-cycle fatigue life of the cellular structures. Based on a 4 g excitation, the calculated bending stresses are approximately 13.5 MPa for the PPR structure and 14.7 MPa for the NPR structure. These stress levels fall within the high-cycle fatigue regime, where Basquin’s law provides a reliable relation between stress amplitude and number of cycles to failure. Therefore, using this law allows for an estimation of the fatigue life under the given operating conditions. Based on these stress values, the predicted fatigue life is roughly 8 years for the PPR harvester, whereas the NPR harvester exhibits an essentially infinite operational lifespan under the same conditions. This substantial difference is attributed to the higher material content in the auxetic cellular substrate, which increases its effective endurance limit and enhances overall durability, making it the superior honeycomb structure in terms of both power generation and long-term performance.

7 Conclusion

This paper presented a multi-objective optimization of a piezoelectric energy harvester (PEH) aimed to maximize the power-to-mass ratio. The harvester considered in this study is a bimorph in the d₃₁ configuration with a honeycomb substrate. This substrate is made of a lightweight material, which reduces the overall mass of the system. Such harvesters are primarily designed to replace batteries and can be deployed in various environments. Therefore, optimizing the unit cell of the substrate layer is essential to enhance the performance of cantilever-type harvesters.

A NN model is developed, and the NSGA-II algorithm is employed to identify the optimal unit cell geometry for the honeycomb positive ratio (PPR) and the honeycomb negative ratio (NPR). The finite element model used to generate the initial dataset for NN training and testing is validated against published literature results. The results demonstrate that the NN models accurately capture the complex behaviour of the harvester, making the problem more effectively solved.

The optimization results clearly demonstrate the critical role of unit cell design in enhancing the performance of honeycomb-based substrates. The optimized honeycomb harvesters outperform conventional solid substrate designs in both power-to-mass ratio and predicted durability, indicating an extended operational lifetime. Concerning the lifetime estimates, they are based on assumptions that provide a first-order approximation of stress and durability and do not capture local failure mechanisms or complex operational effects. Therefore, the results should be considered predictive rather than experimentally confirmed. Nevertheless, the integrated use of finite element models and neural-network-based metamodels allows efficient exploration of a wide design space and identification of high-performance configurations. These findings suggest that both PPR and NPR honeycomb harvesters are expected to achieve longer lifetimes than the solid substrate design, highlighting the potential of microstructural optimization to guide the development of efficient and durable PEH.

References

[1] F K Shaikh and S Zeadally. (2016). Energy harvesting in wireless sensor networks: A comprehensive review. Renewable and Sustainable Energy Reviews, 55:1041–1054. doi: 10.1016/j.rser.2015.11.010.

[2] P Harrop and R Das. (2009). Energy harvesting and storage for electronic devices 2009-2019. IDTechEx.

[3] Y Veli et al. (2023). Numerical analyses on the responses of unimorph and bimorph piezoelectric harvesters in time domain, APME/EMMD, 2022(1):1–7.

[4] M T Mhiri, M Chouchane, M Guerich and W Larbi. (2024). Modeling and analysis of a macro-fiber piezoelectric bimorph energy harvester operating in d33-Mode using Timoshenko theory. Mechanics of Advanced Materials and Structures, 31(30). doi: 10.1080/15376494.2024.2331597.

[5] Pkcell Power. (2025). Exploring the CR2032 battery: a complete guide. Available: https://pkcellpower.com/news/exploring-the-cr2032-battery-a-complete-guide.

[6] M T Mhiri, W Larbi, M Chouchane and M Guerich. (2025). Optimization of a macrofiber piezoelectric energy harvester using artificial neural networks. Composite Structures, 367:119269. doi: 10.1016/j.compstruct.2025.119269.

[7] S Nabavi and L Zhang. (2017). Design and optimization of piezoelectric MEMS vibration energy harvesters based on genetic algorithm. IEEE Sensors J, 17(22):7372–7382. doi: 10.1109/JSEN.2017.2756921.

[8] S Bagheri, N Wu and S Filizadeh. (2020). Application of artificial intelligence and evolutionary algorithms in simulation-based optimal design of a piezoelectric energy harvester. Smart Mater Struct, 29(10):105004. doi: 10.1088/1361-665X/ab9149.

[9] M He, X Zhang, L dos Santos Fernandez, A Molter, L Xia and T Shi. (2021). Multi-material topology optimization of piezoelectric composite structures for energy harvesting. Composite Structures, 265:113783. doi: 10.1016/j.compstruct.2021.113783.

[10] B V De Almeida, D C Cunha and R Pavanello. (2019). Topology optimization of bimorph piezoelectric energy harvesters considering variable electrode location. Smart Mater Struct, 28(8):085030. doi: 10.1088/1361-665X/ab2c3e.

[11] B Chen, Y Jia, F Narita, C Wang and Y Shi (2022). Multifunctional cellular sandwich structures with optimised core topologies for improved mechanical properties and energy harvesting performance. Composites Part B: Engineering, 238:109899. doi: 10.1016/j.compositesb.2022.109899.

[12] Y Jiang et al. (2018). Auxetic mechanical metamaterials to enhance sensitivity of stretchable strain sensors. Advanced Materials, 30(12):1706589. doi: 10.1002/adma.201706589.

[13] A Lazarus and P M Reis. (2015). Soft actuation of structured cylinders through auxetic behavior. Advanced Engineering Materials, 17(6):815–820. doi: 10.1002/adem.201400433.

[14] Q Gao, W-H Liao and L Wang. (2020). On the low-velocity impact responses of auxetic double arrowed honeycomb. Aerospace Science and Technology, 98:105698. doi: 10.1016/j.ast.2020.105698.

[15] A Moayedizadeh and D Younesian. (2022). Application of the meta-substrates for power amplification in rotary piezoelectric energy harvesting systems: Design, fabrication and testing. Energy Reports, 8:5653–5667. doi: 10.1016/j.egyr.2022.04.022.

[16] R Ichige, N Kuriyama, Y Umino, T Tsukamoto and T Suzuki. (2021). Size optimization of metamaterial structure for elastic layer of a piezoelectric vibration energy harvester. Sensors and Actuators A: Physical, 318:112488. doi: 10.1016/j.sna.2020.112488.

[17] M Ravanbod and S Ebrahimi-Nejad. (2024). Perforated auxetic honeycomb booster with re-entrant chirality: A new design for high-efficiency piezoelectric energy harvesting. Mechanics of Advanced Materials and Structures, 31(27):9857–9872. doi: 10.1080/15376494.2023.2280997.

[18] L He, H Kurita, Z Wang and F Narita. (2024). Structural optimization of PVDF cellular resonators for energy-harvesting enhancement based on backpropagation neural network and NSGA-II algorithm. Sensors and Actuators A: Physical, 376:115608. doi: 10.1016/j.sna.2024.115608.

[19] N Chandrasekharan and L L Thompson. (2016). Increased power to weight ratio of piezoelectric energy harvesters through integration of cellular honeycomb structures. Smart Mater Struct, 25(4):045019. doi: 10.1088/0964-1726/25/4/045019.

[20] S Mahboubizadeh, S T Dilamani and S Baghshahi. (2024). Piezoelectricity performance and $\beta$-phase analysis of PVDF composite fibers with BaTiO3 and PZT reinforcement. Heliyon, 10(3):e25021. doi: 10.1016/j.heliyon.2024.e25021.

[21] Q Zhang et al. Bioinspired engineering of honeycomb structure – Using nature to inspire human innovation. Progress in Materials Science, 74:332–400. doi: 10.1016/j.pmatsci.2015.05.001.

[22] Q Li, Y Kuang and M Zhu. (2017). Auxetic piezoelectric energy harvesters for increased electric power output. AIP Advances, 7(1):015104. doi: 10.1063/1.4974310.

[23] S Malek and L Gibson. (2015). Effective elastic properties of periodic hexagonal honeycombs. Mechanics of Materials, 91:226–240. doi: 10.1016/j.mechmat.2015.07.008.

[24] S P Katdare et al. (2025). A review on mechanical metamaterials and additive manufacturing techniques for biomedical applications. Mater Adv, 6(3):887–908. doi: 10.1039/D4MA00874J.

[25] K Qiu, W Ming, L Shen, Q An and M Chen. (2017). Study on the cutting force in machining of aluminum honeycomb core material. Composite Structures, 164:58–67. doi: 10.1016/j.compstruct.2016.12.060.

[26] R Aloui, W Larbi and M Chouchane. (2018). Finite element modeling and analysis of a bimorph piezoelectric energy harvester. Design and modeling of mechanical systems—III. M Haddar, F Chaari, A Benamara, M Chouchane, C Karra and N Aifaoui, Éds. Lecture notes in mechanical engineering. Cham: Springer International Publishing, 1205–1214. doi: 10.1007/978-3-319-66697-6_118.

[27] W Larbi, J-F Deü and R Ohayon. (2012). Finite element formulation of smart piezoelectric composite plates coupled with acoustic fluid. Composite Structures, 94(2):501–509. doi: 10.1016/j.compstruct.2011.08.010.

[28] L P D Silva, J-F Deü, W Larbi and M A Trindade. (2014). An efficient finite element approach for reduction of structural vibration and acoustic radiation by passive shunted piezoelectric systems. 10th World Congress on Computational Mechanics, 2302–2321. doi: 10.5151/meceng-wccm2012-18811.

[29] Y Liao and H A Sodano. (2008). Model of a single mode energy harvester and properties for optimal power generation. Smart Mater Struct, 17(6):065026. doi: 10.1088/0964-1726/17/6/065026.

[30] A K Jain, Jianchang Mao and K M Mohiuddin. (1996). Artificial neural networks: a tutorial. Computer, 29(3):31–44. doi: 10.1109/2.485891.

[31] N R Draper and H Smith. (1998). Applied regression analysis. Hoboken, NJ: Wiley. doi: 10.1002/9781118625590.

[32] Y Fan et al. (2024). Optimization method for the length of the outsourcing concrete working plane on the main arch rib of a rigid-frame arch bridge based on the NSGA-II algorithm. Structures, 59:105767. doi: 10.1016/j.istruc.2023.105767.

[33] K Deb, A Pratap, S Agarwal and T. Meyarivan. (2002). A fast and elitist multiobjective genetic algorithm: NSGA-II. IEEE Trans Evol Computat, 6(2):182–197. doi: 10.1109/4235.996017.

[34] M J Kochenderfer and T A Wheeler. (2019). Algorithms for optimization. Cambridge: The MIT Press.

[35] S Priya and D J Inman Éds. (2009). Energy harvesting technologies. Boston, MA: Springer US. doi: 10.1007/978-0-387-76464-1.

[36] G Palomba, V Crupi and G Epasto. (2019). Collapse modes of aluminium honeycomb sandwich structures under fatigue bending loading. Thin-Walled Structures, 145:106363. doi: 10.1016/j.tws.2019.106363.

[37] G Bianchi, G S Aglietti and G Richardson. (2012). Static and fatigue behaviour of hexagonal honeycomb cores under in-plane shear loads. Appl Compos Mater, 19(2):97–115. doi: 10.1007/s10443-010-9184-5.

[38] B Nečemer, M Vesenjak and S Glodež. (2019). Fatigue of cellular structures – a review. SV-JME, 65(9):525–536. doi: 10.5545/sv-jme.2019.6070.

[39] D Zenkert, A Shipsha and M Burman. (2006). Fatigue of closed cell foams. J Sandwich Structures & Materials, 8(6):517–538. doi: 10.1177/1099636206065886.

[40] L J Gibson. (1981). The elastic and plastic behaviour of cellular materials. University of Cambridge: Doctor of Philosophy. https://www.repository.cam.ac.uk/items/9725785d-9410-48f6-ae05-cc030dec4657.

[41] L J Gibson and M F Ashby. (1997). Cellular solids: structure and properties, 2nd ed. Cambridge: Cambridge University Press. doi: 10.1017/CBO9781139878326.

[42] F Kun, H A Carmona, J S Andrade and H J Herrmann. (2008). Universality behind Basquin’s Law of Fatigue. Phys Rev Lett, 100(9):094301. doi: 10.1103/PhysRevLett.100.094301.

Biographies

Mohamed Taha Mhiri is a PhD candidate in a dual doctoral program in mechanical engineering between the National Engineering School of Monastir (ENIM), Tunisia, and the National Conservatory of Arts and Crafts (CNAM), Paris, France. He earned his mechanical engineering degree from ENIM, Tunisia, in 2022. His research focuses on modelling and optimizing piezoelectric energy harvesters.

Walid Larbi is Full Professor of Mechanics and Civil Engineering at the Conservatoire National des Arts et Métiers (CNAM), Paris, and a member of the LMSSC laboratory. His research interests include structural dynamics, non-linear seismic analysis, reduced-order modelling, and vibro-acoustics. He has authored more than 100 publications and supervised over 12 PhD theses. He is head of the CNAM BTP “Structures” programs, training over 70 engineers per year, and leads CNAM’s international activities in Lebanon, Morocco, and Tunisia. In 2024, he was named Knight in the French Order of Academic Palms.

Mnaouar Chouchane is a Professor of Mechanical Engineering in the National School of Engineering in Monastir, Tunisia, where he has been a faculty member since 1990. He received an Engineering Degree from the University of Tunis in 1984 followed by a master’s degree in mechanical engineering from Washington University in Saint Louis, USA, in 1986 and a PhD from Georgia Institute of Technology, USA, in 1989. His research interests are in the area of vibration and acoustics of machines and structures and smart machines and structures with recent publications on piezoelectric harvesting and passive and active control of rotor vibration. He published more than 20 papers in peer reviewed journals and organized or co-organized the last five editions of the CMSM conference and co-edited the published conference proceedings.

Mohamed Guerich is a Professor of Mechanical Engineering in the Ecole Supérieure d’Ingénieurs Léonard de Vinci (ESILV), France, since 2001. He received a Civil Engineering diploma and a research master’s degree in mechanical engineering from the Ecole Nationale d’Ingénieurs de Tunis (ENIT), Tunisia, in 1991 then a PhD from Université de Technologie de Compiègne (UTC), France, in 1991 followed by an Habilitation to Supervise Research (HDR) from UTC in 2019. His research focuses on the Vibro-Acoustic modelling of structures (Reductions, Damping, Optimizations) and Homogenization of Composite Materials. His recent publications concern the modelling and the optimization of a piezoelectric energy harvester. He has published 14 papers in peer reviewed journals and more than 30 Communications in international conferences with peer review and proceedings.

European Journal of Computational Mechanics, Vol. 34_3&4, 325–362
doi: 10.13052/ejcm2642-2085.34346
© 2026 River Publishers