Towards advanced piezoelectric metamaterial design via combined topology and shape optimization

Abstract

Metamaterials open up a spectrum of artificially engineered properties otherwise unreachable in conventional bulk materials. For electromechanical energy conversion systems, lightweight materials with high hydrostatic piezoelectric coupling coefficients and negative Poisson’s ratio can be obtained. Thus, in this contribution, we explore the possibilities of piezoelectric metamaterials design by employing structural optimization. More specifically, we apply a sequential framework of topology and shape optimization to design piezoelectric metamaterials with negative Poisson’s ratio for electromechanical energy conversion under uniform pressure. Topology optimization is employed to generate the initial layout, whereas shape optimization fine tunes the design and improves durability and manufacturability of the structures with the help of a curvature constraint. An embedding domain discretization (EDD) method with adaptive domain and shape refinement is utilized for an efficient and accurate computation of the state problem in the shape optimization stage. Multiple case studies are conducted to determine the importance of desired stiffness characteristics, symmetry conditions and objective formulations on the design of piezoelectric metamaterials. Results show that the obtained designs are highly dependent on the desired stiffness characteristics. Moreover, the addition of the EDD-based shape optimization step introduces significant changes to the designs, confirming the usability of the sequential framework.

1 Introduction

Metamaterials have seen growing interest in recent years due to the their exceptional properties that are otherwise unreachable in natural or conventional materials. Such properties can be achieved via an engineered geometry on the microstructural level, which, when periodically repeated, can reproduce the effective properties on the macroscopic level. The broad spectrum of reachable physical properties includes, e.g. acoustic and piezoelectric acoustic metamaterials for wave manipulation and enhanced sound absorption and insulation (Ji and Huber 2022), photonic metamaterials for ultrahigh-resolution imaging or cloaking devices (Soukoulis and Wegener 2011), piezoelectric metamaterials for vibration suppression and energy harvesting at wide bandgaps (Hu et al. 2017), electromagnetic or chiral metamaterials with applications in remote aerospace devices, military equipment or medical diagnostic devices (Valipour et al. 2022).

Piezoelectrics are particularly interesting from the perspective of metamaterial design due to their capability of electromechanical energy conversion. Therefore, they find application in actuation, sensing (Shaikh and Zeadally 2016), hydrophones, energy harvesting (Wei and Jing 2017) and many other engineering and medical devices (Vijaya 2012). Some of the inherent limitations of ceramic piezoelectrics, like brittleness and high stiffness, can be overcome, for instance, by introducing porosity in the microstructure (Della and Shu 2008; Kar-Gupta and Venkatesh 2006; Zhang et al. 2017). Moreover, progress in additive manufacturing has a crucial contribution in the development of piezoelectric materials with arbitrary geometries on the microscale level, including lead-free piezoelectrics like PVDF nanocomposites (Bodkhe et al. 2017) or BCZT cellular materials (Köllner et al. 2022). Porous piezoelectrics exhibit lower stiffness, brittleness and higher sensitivity, which is highly beneficial for ultrasonic transduction (Smith 1989) or ceramics for biomechanical structures (Mancuso et al. 2021).

By altering the geometry of the microstructure, some mechanical properties like stiffness or density can be tailored for specific applications but most importantly, materials with unconventional characteristics like negative effective Poisson’s ratio, so-called auxetic materials, can be obtained (Alderson and Alderson 2007; Carneiro et al. 2013). The exceptional ability of auxetic materials to contract in transverse directions under compressive load opens up a range of possibilities for energy harvesting (Zhao et al. 2022). Auxetic cellular honeycomb (HC) structures are the most popular auxetic metamaterials, widely studied in the literature due to their relatively simple geometry (Iyer and Venkatesh 2011; Iyer et al. 2015; Khan et al. 2019). Stress build-up within auxetic HC was verified numerically and experimentally (Whitty et al. 2003) and various configurations of auxetic HC were compared for their superior electromechanical properties towards energy harvesting applications (Köllner et al. 2023).

In the last two decades, the role of structural optimization methods, density-based and level-set topology optimization in particular (Bendsoe and Sigmund 2003; Wang et al. 2003), has immensely increased in engineering design. Geometries generated using structural optimization reach beyond designers intuition and broaden our horizons in understanding the relationship between shape, connectivity and the resultant properties of the structure. Advances in additive manufacturing only strengthen these developments, since they allow manufacturing of topology optimized shapes without engineer’s interference with the obtained structures. These developments have been successfully extended towards the generation of new geometries for metamaterials. Besides the classical, re-entrant design of HC for auxetic metamaterials, alternative designs of periodic microstructures have been generated mathematically using topology optimization for linear elastic materials (Neves et al. 2000; Zhang et al. 2007) with the objective to generate structures with desired stiffness properties, in Andreassen et al. (2014) with consideration of manufacturability and 3D geometries.

In Paulino et al. (2009), the optimal design of functionally graded composites was shown, which benefits from the presence of grey transition regions within density-based topology optimization. Topology optimization of thermomechanical microstructures was addressed in Moshki et al. (2022) using an alternating active-phase algorithm to handle the presence of multiple materials. Optimal design of biodegradable metamaterials for medical applications was addressed in Zhang et al. (2021). Alternative optimization methods like the ground structure approach (Li et al. 2021) or machine learning frameworks (Challapalli et al. 2021; Fernández et al. 2022) have been applied to generate optimal lattice structures or density-based structures (Kollmann et al. 2020).

Optimal microstructure designs for piezoelectric materials were shown in Silva et al. (1997), Sigmund et al. (1998) and Silva and Kikuchi (1999) with the objective function to maximize the hydrostatic coupling coefficient $d_h$ or the figure of merit $d_h \times g_h$. These studies on the design of piezoelectric microstructures using topology optimization give a detailed overview of how the results are affected by the choice of different objective functions focusing on energy conversion performance. The designs obtained in these works resulted in auxetic behaviour, proving that $d_h$ is a good choice of objective function for hydrostatic applications. In Vatanabe et al. (2014), the design of piezoelectric microstructures was successfully explored using topology optimization for maximum phononic bandgaps. In Xu et al. (2017), the dynamic characteristics of piezoelectric microstructure were optimized using a bidirectional evolutionary structure optimization (BESO) method. The goal of this work was to maximize the damping dissipation velocity, which improves the active control performance. In Ding and Xu (2022), the design of piezoelectric microstructures for vibration suppression using the BESO method was introduced by considering the initial disturbance.

In addition to microstructural design for piezoelectric structures, recent advances in te optimal design of piezoelectric structures include reliability-based topology optimization with voltage uncertainty (Yang et al. 2022), which uses the PEMAP-P model that additionally considers the polarization direction as a design variable. In Donoso et al. (2023), a novel piezoelectric transducer design methodology is presented that focuses on optimizing electrode configurations for manufacturable and functional devices. In Wang et al. (2023), a velocity field level-set method is presented to determine the optimal material distribution of a piezoelectric layer on a plate with active vibration control. Compared to density-based topology optimization, the level-set method improved the accuracy of vibration control by avoiding intermediate densities.

Among various structural optimization methods the density-based approach to topology optimization (Bendsøe 1989), also referred to as the SIMP (Solid Isotropic Material with Penalization) method, is most widely studied due to its simple formulation and implementation (Andreassen et al. 2011; Aage et al. 2015), large design space (no initial guess of the design is necessary) and robustness. However, by using the material interpolation law and the lack of mechanical boundaries, it also poses certain challenges, in particular in imposing boundary-related constraints, e.g. curvature constraints. On the contrary, node-based shape optimization methods rely upon the direct update of the boundary of the structure (Ding 1986). The presence of the mechanical boundary guarantees an accurate solution of the state problem and enables simple incorporation of boundary-related constraints, like stress or curvature constraints. The main hurdle, however, is the limitation to small design updates and the necessity of a valid initial guess of the design. The standard approach to node-based shape optimization deals with an explicitly meshed structure and the design update concerns not only the boundary nodes (the factual design variables) but also the inner nodes (the dependent nodes), thereby greatly increasing the complexity of the regularization schemes. Some of these drawbacks are addressed in implicit boundary approaches, like, for instance, the level-set methods (Allaire et al. 2002), where larger design updates are not an issue due to the separation of the computational domain, which consists of a structured mesh, and the representation of the shape by a level-set function. Alternatively, node-based shape optimization using an embedding domain discretization (EDD) method was outlined in Riehl and Steinmann (2017), in which the variable shape is explicitly modelled by lower-dimensional segments, e.g. line segments in 2D; hence, the inner nodes of the variable shape are avoided. Moreover, a structured, adaptively refined background mesh is utilized as the computational domain. This approach alleviates the update of the inner elements of the variable shape and mitigates the complexity of regularization schemes, effectively enabling much larger design updates and at the same time facilitating more efficient computations with flexible refinement strategies for the background mesh.

To profit from the advantages of different structural optimization approaches, various combined methods were proposed (Sokolowski and Zochowski 1999; Norato et al. 2007; Eschenauer et al. 1994; Christiansen et al. 2014; Lian et al. 2017; Riehl and Steinmann 2015; Nguyen et al. 2020; Andreasen et al. 2020; Stankiewicz et al. 2021). In our work, we adapt the sequential approach as in Dev et al. (2022); Stankiewicz et al. (2022), in which a standard, density-based topology optimization is used to generate an initial structure and node-based shape optimization using EDD follows up to fine tune the design. This combination proved to be a good fit due to a number of reasons. First, the choice of EDD-based shape optimization requires generation of only the boundary of the shape (1D segments) after the topology optimization step, which greatly simplifies the transition step from topology to shape optimization. Secondly, a coarsely discretized domain for topology optimization is chosen, as we do not require highly accurate state problem computations, since it only provides the initial configuration of the shape optimization step. Next, we incorporate a shape curvature constraint as in Stankiewicz et al. (2022) to improve manufacturability and durability of the final design, which is otherwise a challenging task in topology optimization due to the lack of an explicit boundary. Finally, we exploit the advantages of both topology and shape optimization, as listed before, to provide a versatile and powerful structural optimization routine capable of designing complex geometries from a large design space, at the same time guaranteeing high accuracy of the state problem solutions and a variety of possible constraints.

The contribution of this work is twofold. Firstly, we develop a shape optimization methodology for the design of metamaterials by exploiting computational homogenization within the EDD framework. We investigate how shape optimization further improves the metamaterial after generating an initial design with topology optimization to verify the significance of such sequential approach. Secondly, we perform an advanced study on how the microstructural design of piezoelectric structures depends on the desired stiffness characteristics by running case studies with various constraint sets on both the “normal” stiffness coefficients as well as the “shear” stiffness coefficients. Moreover, a comparison of different objective functions for the generation of auxetic structures and of various symmetry conditions is performed.

The article is organized as follows. In Sect. 2 we introduce the theoretical background for the computation of effective materials properties of piezoelectric structures. We employ asymptotic expansion in order to provide analytical expressions for the effective materials tensors and for the computation of sensitivities. In Sect. 3 we outline the sequential optimization routine. Section 4 contains the results of the case studies and the discussion. The conclusions and future outlook are addressed in Sect. 5.

2 Theoretical background

2.1 Piezoelectricity

The direct piezoelectric effect is observed when electric charges accumulate on the surface of a non-conducting dielectric material due to mechanical strain. The converse piezoelectric effect occurs when mechanical stress is induced in response to an applied electric field.

The following theoretical formulations are outlined based on Ieee standard on piezoelectricity (1988). We adopt a standard tensor notation. For the formulation of the state problem, we consider a piezoelectric domain $\Omega$. At any point in the domain $\Omega$ , the mechanical balance of linear momentum has to be fulfiled, which, assuming there are no body forces acting on $\Omega$, reads

$$\begin{aligned} \text {div} \varvec{\sigma } = \varvec{0} \qquad \text {in} \ \Omega \end{aligned}$$

(1)

and where $\varvec{\sigma }$ is the Cauchy stress. Additionally, the electrostatic equation is formulated for an insulating material as follows;

$$\begin{aligned} \text {div} \varvec{D} = 0 \qquad \text {in} \ \Omega \end{aligned}$$

(2)

where $\varvec{D}$ is the electric displacement. Note that for the electric behaviour, we are only concerned with the quasistatic electric approximation; hence, no dynamic coupling with mechanical strains is considered. Mechanical strain is obtained as the symmetric gradient of the displacement field $\varvec{u}$

$$\begin{aligned} \varvec{\varepsilon }= \frac{1}{2}\left[ \nabla \varvec{u} + \nabla ^T \varvec{u} \right] \end{aligned}$$

(3)

and the electric field is derived from a scalar electric potential $\phi$

$$\begin{aligned} \varvec{E}= - \nabla \phi \end{aligned}$$

(4)

By considering the mechanical and electrical properties of the structure as well as the piezoelectric coupling between them, the enthalpy density for a piezoelectric structure takes the form

$$\begin{aligned} H = \frac{1}{2}\varvec{\varepsilon }: \varvec{C}: \varvec{\varepsilon }- \varvec{E}\cdot \varvec{e}: \varvec{\varepsilon }- \frac{1}{2}\varvec{E}\cdot \varvec{\epsilon }^S\cdot \varvec{E}\end{aligned}$$

(5)

where $\varvec{C}$ is the elasticity tensor, $\varvec{e}$ is the piezoelectric coupling tensor and $\varvec{\epsilon }^S$ is the permittivity tensor. The linear constitutive relations are then obtained by taking derivatives of H w.r.t the field gradients

$$\begin{aligned}{} & {} \varvec{\sigma }= \frac{\partial {H}}{\partial {\varvec{\varepsilon }}} = \varvec{C}: \varvec{\varepsilon }- \varvec{E}\cdot \varvec{e} \end{aligned}$$

(6)

$$\begin{aligned}{} & {} \varvec{D}= -\frac{\partial {H}}{\partial {\varvec{E}}} = \varvec{e}: \varvec{\varepsilon }+ \varvec{\epsilon }^S\cdot \varvec{E} \end{aligned}$$

(7)

Equations (6) and (7) are represented in a stress-charge form, in which the elasticity tensor corresponds to the material stiffness, the piezoelectric coupling tensor and the permittivity tensor define the relation at constant strain between the electric field and the Cauchy stress and the electric displacement, accordingly. Alternatively, the constitutive relations can be represented in a strain-charge form, namely

$$\begin{aligned}{} & {} \varvec{\varepsilon }= \varvec{S}: \varvec{\sigma }+ \varvec{E}\cdot \varvec{d} \end{aligned}$$

(8)

$$\begin{aligned}{} & {} \varvec{D}= \varvec{d}: \varvec{\sigma }+ \varvec{\epsilon }^T\cdot \varvec{E} \end{aligned}$$

(9)

where $\varvec{S}= \left[ \varvec{C}\right] ^{-1}$ is the compliance matrix, $\varvec{d}= \varvec{e}: \varvec{S}$ the piezoelectric tensor and $\varvec{\epsilon }^T$ the permittivity tensor in the strain-charge format

$$\begin{aligned} \varvec{\epsilon }^T= \varvec{e}: \varvec{S}: \varvec{e}+ \varvec{\epsilon }^S \end{aligned}$$

(10)

The boundary conditions for a piezoelectric problem are given as follows

$$\begin{aligned} \begin{aligned} \varvec{\sigma }\cdot \varvec{n} = \bar{\varvec{t}} \qquad&\text {on} \ \partial \Omega _N \\ \varvec{D}\cdot \varvec{n} = \bar{q} \qquad&\text {on} \ \partial \Omega _N \\ \varvec{u}= {\bar{\varvec{u}}} \qquad&\text {on} \ \partial \Omega _D \\ \phi = {\bar{\phi }} \qquad&\text {on} \ \partial \Omega _D \end{aligned} \end{aligned}$$

(11)

where $\bar{\varvec{t}}$ and $\bar{q}$ are prescribed tractions and free surface charges, respectively.

2.2 Homogenization for piezoelectrics

In the following, we briefly introduce the methodology to compute effective material properties of piezoelectric composites. In this work, the effective material properties are predicted numerically using the finite element method. For structural optimization purposes, the effective material properties require an analytical representation (Silva et al. 1997). Therefore, a first-order asymptotic expansion of the field variables (Chung et al. 2001) is employed according to the theory of homogenization in order to separate the micro- and macroscale:

$$\begin{aligned}{} & {} \varvec{u}^\epsilon = \varvec{u}_0 \left( \varvec{x} \right) + \epsilon \varvec{u}_1 \left( \varvec{x}, \varvec{y} \right) \end{aligned}$$

(12)

$$\begin{aligned}{} & {} \phi ^\epsilon = \phi _0 \left( \varvec{x} \right) + \epsilon \phi _1 \left( \varvec{x}, \varvec{y} \right) \end{aligned}$$

(13)

where $\varvec{y} = \varvec{x} / \epsilon$ ($\epsilon$ is a small parameter) represents an upscaling of the microstructure that captures the heterogeneous nature of the material. $\varvec{u}_0$ and $\phi _0$ are independent of $\varvec{y}$, since they describe the macroscopic (homogenized) part of the solution, while $\varvec{u}_1$ and $\phi _1$ represent the fluctuation part of the solution due to the material variation on the microscale. Hence, a portion of the microstructure is selected, the so-called unit cell or the representative volume element (RVE) $\text {Y}= [0,y_1] \times [0,y_2] \times [0,y_3]$ that best approximates the complex behaviour of the structure when periodically repeated in all directions. The fluctuations of the solution can be further formulated as follows (Gałka et al. 1992)

$$\begin{aligned}{} & {} \varvec{u}_1 = \varvec{\chi } \left( \varvec{x}, \varvec{y} \right) : \varvec{\varepsilon }\left( \varvec{u}_0\right) - \varvec{\zeta } \left( \varvec{x}, \varvec{y} \right) \cdot \varvec{E}\left( \phi _0 \right) \end{aligned}$$

(14)

$$\begin{aligned}{} & {} \phi _1 = \varvec{\psi } \left( \varvec{x}, \varvec{y} \right) : \varvec{\varepsilon }\left( \varvec{u}_0\right) - \varvec{\eta } \left( \varvec{x}, \varvec{y} \right) \cdot \varvec{E}\left( \phi _0 \right) \end{aligned}$$

(15)

where $\varvec{\chi }$ is a third-order tensor consisting of characteristic displacements $\varvec{\chi }^\mathrm{(ij)}$ (Torquato and Haslach Jr 2002), $\varvec{\zeta }$ and $\varvec{\psi }$ are second-order tensors, where $\varvec{\zeta }^\mathrm{(k)}$ and $\psi ^\mathrm{(ij)}$ are characteristic coupling functions and $\varvec{\eta }$ is a first-order tensor, where $\eta ^\mathrm{(k)}$ are characteristic potentials. Indexes $\mathrm{(ij)}$ and $\mathrm{(k)}$ indicate a relation to the corresponding unit tests, i.e. the $\mathrm{(ij)}$-th unit strain tests and $\mathrm{(k)}$-th unit electric field tests.

To obtain the respective microscopic strain and electric field tensors, we utilize the chain rule

$$\begin{aligned} \frac{\partial }{\partial \varvec{x}^\epsilon } = \frac{\partial }{\partial \varvec{x}} + \frac{1}{\epsilon } \frac{\partial }{\partial \varvec{y}} \end{aligned}$$

(16)

By assuming that $\epsilon$ is sufficiently small, we obtain the following kinematic relation for the partial derivative of the displacement field

$$\begin{aligned} \frac{\partial {\varvec{u}^\epsilon }}{\partial {\varvec{x}^\epsilon }} = \frac{\partial {\varvec{u}_0}}{\partial {\varvec{x}}} + \epsilon \frac{\partial {\varvec{u}_1}}{\partial {\varvec{x}}} + \frac{\partial {\varvec{u}_1}}{\partial {\varvec{y}}} \approx \frac{\partial {\varvec{u}_0}}{\partial {\varvec{x}}} + \frac{\partial {\varvec{u}_1}}{\partial {\varvec{y}}} \end{aligned}$$

(17)

and by employing Eqs. (14) and (15), we arrive at

$$\begin{aligned}{} & {} \varvec{\varepsilon }^{\epsilon } = \nabla _{\varvec{x}^\epsilon }^\textrm{sym} \varvec{u}^\epsilon \approx \left[ \varvec{I} + \nabla _{\varvec{y}}^\textrm{sym} \varvec{\chi } \right] : \varvec{\varepsilon }_0- \nabla _{\varvec{y}}^\textrm{sym} \varvec{\zeta } \cdot \varvec{E}_0 \end{aligned}$$

(18)

$$\begin{aligned}{} & {} \varvec{E}^{\epsilon } = - \frac{\partial {\phi ^\epsilon }}{\partial {\varvec{x}^\epsilon }} \approx \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \cdot \varvec{E}_0- \frac{\partial {\varvec{\psi }}}{\partial {\varvec{y}}}: \varvec{\varepsilon }_0 \end{aligned}$$

(19)

At this point, we can express the total strains and total electric fields for each of the (ij)-th and (k)-th unit tests as

$$\begin{aligned}{} & {} \varvec{\varepsilon }^\mathrm{(ij)} \approx \left[ \varvec{I} + \nabla _{\varvec{y}}^\textrm{sym} \varvec{\chi } \right] : \varvec{\varepsilon }_0^{(\mathrm ij)} \end{aligned}$$

(20)

$$\begin{aligned}{} & {} \varvec{\varepsilon }^\mathrm{(k)} \approx -\nabla _{\varvec{y}}^\textrm{sym} \varvec{\zeta } \cdot \varvec{E}_0^{(\mathrm k)} \end{aligned}$$

(21)

$$\begin{aligned}{} & {} \varvec{E}^\mathrm{(ij)} \approx - \frac{\partial {\varvec{\psi }}}{\partial {\varvec{y}}}: \varvec{\varepsilon }_0^{(\mathrm ij)} \end{aligned}$$

(22)

$$\begin{aligned}{} & {} \varvec{E}^\mathrm{(k)} \approx \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \cdot \varvec{E}_0^{(\mathrm k)} \end{aligned}$$

(23)

where $\varvec{\varepsilon }_0^\mathrm{(ij)}$ and $\varvec{E}_0^\mathrm{(k)}$ are the macroscopic (unit) strain and electric field due to $\mathrm{(ij)}$-th and $\mathrm{(k)}$-th unit tests. For a detailed derivation of the microscopic and macroscopic equations, we refer to (Gałka et al. 1992; Nelli Silva et al. 1999). In the following, we directly proceed to the homogenized tensors in the energy bilinear form. By utilizing the previously introduced quantities, they are given as follows.

Fig. 1

d21 (d31), d22 (d33) and the combined $d^h$ conversion on the example of the BTO crystal structure. The titanium atom (in red) is slightly off-centre in the vertical direction, indicating the polarization direction. When external loading is applied, the deformation of the structure causes the titanium atom to further displace, affecting the dipole moment and the resultant electric displacement $\varvec{D}$. The light red atom indicates its original position before the application of the external load. (Color figure online)

$$\begin{aligned}{} & {} \begin{aligned} \varvec{C}_{\textrm{H}} =&\left\langle \left[ \varvec{I} + \nabla _{\varvec{y}}^\textrm{sym} \varvec{\chi } \right] : \varvec{C}: \left[ \varvec{I} + \nabla _{\varvec{y}}^\textrm{sym} \varvec{\chi } \right] \right. \\&+\left. \frac{\partial {\varvec{\psi }}}{\partial {\varvec{y}}} \cdot \varvec{e}: \left[ \varvec{I} + \nabla _{\varvec{y}}^\textrm{sym} \varvec{\chi } \right] \right\rangle \end{aligned} \end{aligned}$$

(24)

$$\begin{aligned}{} & {} \begin{aligned} \varvec{e}_{\textrm{H}} =&\left\langle \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \cdot \varvec{e}: \left[ \varvec{I} + \nabla _{\varvec{y}}^\textrm{sym} \varvec{\chi } \right] \right. \\&-\left. \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \cdot \varvec{\epsilon }^S\cdot \frac{\partial {\varvec{\psi }}}{\partial {\varvec{y}}}\right\rangle \end{aligned} \end{aligned}$$

(25)

$$\begin{aligned}{} & {} \begin{aligned} \varvec{\epsilon }^S_{\textrm{H}} =&\left\langle - \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \cdot \varvec{e}: \frac{\partial {\varvec{\zeta }}}{\partial {\varvec{y}}} \right. \\&+\left. \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \cdot \varvec{\epsilon }^S\cdot \left[ \varvec{1} + \frac{\partial {\varvec{\eta }}}{\partial {\varvec{y}}} \right] \right\rangle \end{aligned} \end{aligned}$$

(26)

whereas by using Eqs. (20)–(23) we obtain the expressions for individual material coefficients as follows

$$\begin{aligned}{} & {} \begin{aligned} C^\mathrm{(ijkl)}_{\textrm{H}} =&\ \varvec{\varepsilon }_0^\mathrm{(ij)}: \varvec{C}_{\textrm{H}}: \varvec{\varepsilon }_0^\mathrm{(kl)}\\ =&\ \langle \varvec{\varepsilon }^\mathrm{(ij)}: \varvec{C}: \varvec{\varepsilon }^\mathrm{(kl)} - \varvec{E}^\mathrm{(ij)} \cdot \varvec{e}: \varvec{\varepsilon }^\mathrm{(kl)} \rangle \end{aligned} \end{aligned}$$

(27)

$$\begin{aligned}{} & {} \begin{aligned} e^\mathrm{(kij)}_{\textrm{H}} =&\ \varvec{E}_0^\mathrm{(k)} \cdot \varvec{e}_{\textrm{H}}: \varvec{\varepsilon }_0^\mathrm{(ij)}\\ =&\ \langle \varvec{E}^\mathrm{(k)} \cdot \varvec{e}: \varvec{\varepsilon }^\mathrm{(ij)} + \varvec{E}^\mathrm{(k)} \cdot \varvec{\epsilon }^S\cdot \varvec{E}^\mathrm{(ij)} \rangle \end{aligned} \end{aligned}$$

(28)

$$\begin{aligned}{} & {} \begin{aligned} \epsilon ^\mathrm{S \ (kl)}_{\textrm{H}} =&\ \varvec{E}_0^\mathrm{(k)} \cdot \varvec{\epsilon }^S_{\textrm{H}} \cdot \varvec{E}_0^\mathrm{(l)}\\ =&\ \langle \varvec{E}^\mathrm{(k)} \cdot \varvec{e}: \varvec{\varepsilon }^\mathrm{(l)} + \varvec{E}^\mathrm{(k)} \cdot \varvec{\epsilon }^S\cdot \varvec{E}^\mathrm{(l)} \rangle \end{aligned} \end{aligned}$$

(29)

2.3 Chosen performance measures

Two performance measures are considered in this study for the objective function: the effective Poisson’s ratio $\upsilon _{\textrm{H}}$ and the hydrostatic coupling (charge) coefficient $d^h_{\textrm{H}}$. The purpose of this work is to design structures with improved energy conversion efficiency under uniform pressure field. The works of Silva et al. (1997), Sigmund et al. (1998), Silva and Kikuchi (1999) have shown that using the hydrostatic coupling (charge) coefficient $d^h_{\textrm{H}}$ as an objective function leads to feasible, auxetic designs with good energy conversion efficiency under uniform pressure field. However, for better comparison, $\upsilon _{\textrm{H}}$ is alternatively employed in the topology optimization step with the goal of initiating the shape optimization step with an auxetic structure. The choice of $\upsilon _{\textrm{H}}$ or $d^h_{\textrm{H}}$ is then examined to determine the differences in the topologies obtained using these objective functions. However, only $d^h_{\textrm{H}}$ is used for the shape optimization step since the initial structure is already auxetic at this point.

Before proceeding with the definition of those measures, let us reformulate the effective material tensors using Voigt notation. The effective elastic tensor takes the form of a 3x3 matrix in two dimensions (for the two-dimensional case in this work we consider the direction 2 (Y) as the polarized direction) and for an orthotropic material is given by

$$\begin{aligned} \varvec{C}_{\textrm{H}} = \begin{bmatrix} C^{11}_{\textrm{H}} &{}\quad C^{12}_{\textrm{H}} &{}\quad 0\\ C^{12}_{\textrm{H}} &{}\quad C^{22}_{\textrm{H}} &{}\quad 0\\ 0 &{}\quad 0 &{}\quad C^{33}_{\textrm{H}} \end{bmatrix} \end{aligned}$$

(30)

Whereas the effective piezoelectric tensor $\varvec{e}_{\textrm{H}}$ in the stress-charge form in the Voigt notation is defined as

$$\begin{aligned} \varvec{e}_{\textrm{H}} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad e^{13}_{\textrm{H}}\\ e^{21}_{\textrm{H}} &{}\quad e^{22}_{\textrm{H}} &{}\quad 0 \end{bmatrix} \end{aligned}$$

(31)

Analogously, the effective piezoelectric tensor in the strain-charge form is given by

$$\begin{aligned} \varvec{d}_{\textrm{H}} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad d^{13}_{\textrm{H}}\\ d^{21}_{\textrm{H}} &{}\quad d^{22}_{\textrm{H}} &{}\quad 0 \end{bmatrix} \end{aligned}$$

(32)

Now, the effective Poisson’s ratio is formulated by assuming a transversely isotropic behaviour on the macroscale

$$\begin{aligned} \upsilon _{\textrm{H}} = \upsilon _{\textrm{H}}^{12} = \frac{2 C_{\textrm{H}}^{12}}{C_{\textrm{H}}^{11} + C_{\textrm{H}}^{22}} \end{aligned}$$

(33)

The second quantity of interest is the effective hydrostatic (charge) coupling coefficient $d_{\textrm{H}}^h$ which relates the electric displacement in the polarization direction $D^2_{\textrm{H}}$ with a hydrostatic pressure P under zero electric field $\varvec{E}_0 = \varvec{0}$. For clarity, we recall Eq. (9) under the assumption of zero electric field

$$\begin{aligned} \varvec{D}_{\textrm{H}} = \varvec{d}_{\textrm{H}}: \varvec{\sigma } \end{aligned}$$

(34)

Then, by considering a hydrostatic pressure (in Voigt notation) $\varvec{\sigma }= \left[ P, P, 0\right]$ , we obtain the following relation

$$\begin{aligned} \begin{bmatrix} D^1_{\textrm{H}}\\ D^2_{\textrm{H}} \end{bmatrix} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad d^{13}_{\textrm{H}}\\ d^{21}_{\textrm{H}} &{}\quad d^{22}_{\textrm{H}} &{}\quad 0 \end{bmatrix} \begin{bmatrix} P\\ P\\ 0 \end{bmatrix} \end{aligned}$$

(35)

effectively reducing to

$$\begin{aligned} D^2_{\textrm{H}} = P \left[ d^{21}_{\textrm{H}} + d^{22}_{\textrm{H}}\right] \end{aligned}$$

(36)

since the (perpendicular to the polarization direction) electric displacement $D^1_{\textrm{H}}$ is only related to $\varvec{\sigma }_{\textrm{H}}$ by $d^{31}_{\textrm{H}}$, effectively coupling only the shear stress $\sigma ^{12}_{\textrm{H}}$ with $D^1_{\textrm{H}}$, which is set to zero in the case of hydrostatic pressure. Thus, the sum $d^h_{\textrm{H}} = d^{21}_{\textrm{H}} + d^{22}_{\textrm{H}}$ is the effective hydrostatic coupling (charge) coefficient of interest, which can be formulated as follows:

$$\begin{aligned} d_{\textrm{H}}^h = d^{21}_{\textrm{H}} + d^{22}_{\textrm{H}} = \bar{\varvec{v}} \cdot \left[ \varvec{\bar{C}}_{\textrm{H}}\right] ^{-1} \cdot \varvec{\bar{e}}_{\textrm{H}} \end{aligned}$$

(37)

where $\bar{\varvec{v}} = \left[ 1 \ 1\right]$ and $\varvec{\bar{C}}_{\textrm{H}}$, $\varvec{\bar{e}}_{\textrm{H}}$ are the simplified effective tensors, given as follows

$$\begin{aligned} \varvec{\bar{C}}_{\textrm{H}} = \begin{bmatrix} C^{11}_{\textrm{H}} &{}\quad C^{12}_{\textrm{H}}\\ C^{12}_{\textrm{H}} &{}\quad C^{22}_{\textrm{H}} \end{bmatrix} \end{aligned}$$

(38)

$$\begin{aligned} \varvec{\bar{e}}_{\textrm{H}} = \begin{bmatrix} e^{21}_{\textrm{H}}\\ e^{22}_{\textrm{H}} \end{bmatrix} \end{aligned}$$

(39)

Here, the second formulation involves only the material tensors in the stress-charge form, as required for the sensitivity analysis. For a clearer understanding of the meaning behind the hydrostatic piezoelectric coupling $d^h$, in Fig. 1, the two possible ways of invoking the direct piezoelectric effect on the example of a BTO non-centro-symmetric crystal structure are shown, the sum of which make up for the $d^h$ conversion. The slight displacement of the oxygen octahedron and the titanium atom along the direction 2 (Y) creates a dipole moment, which is further affected by the external loading and the resultant deformation. The left-hand image depicts the $d^{21}$ ($d^{31}$ in 3D) direct conversion, since the load is applied in the direction 1 (X), perpendicular to the polarization direction. In the second image, the material is loaded along the polarization direction, hence the $d^{22}$ ($d^{33}$ in 3D) conversion. However, under the hydrostatic loading conditions (as in the right-hand image), the effects of $d^{21}$ and $d^{22}$ conversion annihilate each other, effectively resulting in a very low conversion rate. This is a clear consequence of the fact that any bulk material exhibits a positive Poisson’s ratio. Hence, application of the metamaterial concept to design structures with negative Poisson’s ratio appears to be reasonable in order to account for higher conversion rates under hydrostatic loading.

3 Optimization strategy

3.1 Sequential approach

In Dev et al. (2022), a sequential approach to topology and shape optimization was outlined, specifically to tackle the challenge of designing compliant mechanisms with a focus on durable flexure hinges, which require highly refined geometry and accurate numerical models. Topology optimization is the state of the art tool for structural design and has the capability to generate complex structures without initial guess of the geometry. Moreover, it is relatively simple to implement and robust. Nevertheless, shape optimization explicitly treats the boundary of the structure, providing an opportunity to consider boundary-related measures, like, for instance, the boundary stress or the curvature. In Fig. 2, the sequential approach to topology and shape optimization is sketched together with a brief overview of the advantages of each method. For detailed discussions on the applicability of this approach, refer to Dev et al. (2022), Stankiewicz et al. (2022).

Fig. 2

Optimization strategy

3.2 Topology optimization

Density-based topology optimization employs pseudo-densities $\rho _e \in \left[ 0, 1 \right]$ as design variables, which are assigned to every finite element, hence the subindex e. For pure elastic problems, SIMP (Solid Isotropic Material with Penalization) is the standard approach to determine the material properties of a structure with intermediate densities. For piezoelectric problems, the interpolation of each material tensor is analogous $\{ \bullet \} \left( {\hat{\rho }}_e \right) = \left[ {\hat{\rho }}^{p}_e \left[ 1 - \rho _{\textrm{min}}\right] + \rho _{\textrm{min}}\right] \ \{ \bullet \}$, where $\rho _{\textrm{min}} = 10^{-9}$ and for each material tensor we employ the following penalization factors: $p = 3$ for $\varvec{C}_e$ and $p = 2$ for $\varvec{e}_e$ and $\varvec{\epsilon }^S_e$. We determined that this combination of penalization factors provides, for the case studies in this work, the best convergence towards black and white designs. The hat in ${\hat{\rho }}_e$ indicates the regularized pseudo-density. The three-step regularization consists of the density filtering for member size control and checkboard pattern avoidance (Bruns and Tortorelli 2001), smoothed Heaviside projection filter for the reduction of the grey transition region (Wang et al. 2011) and projection-based filter for symmetry constraints (Vatanabe et al. 2016).

The topology optimization step is defined according to the following setup

$$\begin{aligned} \begin{aligned} {\mathcal {F}}\left( {{\hat{\rho }}, \varvec{u}} \right) \Rightarrow&\ \ \text {Min} \\ \text {st.} \ \ \ {\mathcal {G}}_{\textrm{vol}} =&\int _{\Omega }\ {\hat{\rho }}\left( {\varvec{x}} \right) \text {\,dV}- \overline{V} \le \ 0 \\ C_{\textrm{H}}^\textrm{ij} \ge&\ \bar{C}^\textrm{ij} \ \ \ \text {for \ \ i,j} \in \{ 1,2,3 \}\\ 0 \le&\ \rho \le 1 \end{aligned} \end{aligned}$$

(40)

Both of the objective functions given in 2.3 are considered in this step: the effective Poisson’s ratio $\upsilon _{\textrm{H}}$ and the effective hydrostatic piezoelectric coupling coefficient $d_{\textrm{H}}^h$. However, ${\mathcal {F}}\left( {{\hat{\rho }}, \varvec{u}} \right)$ is not a multi-objective function, but rather, it is either $\upsilon _{\textrm{H}}$ or $d_{\textrm{H}}^h$ as shown in the results section. The minimization of the first objective is employed purely to generate an auxetic behaviour, whereas maximization of the second objective directly optimizes the electromechanical conversion efficiency under hydrostatic loading conditions. Such a choice of objectives allows to assess whether those responses can be used interchangeably, since both of the formulations should ultimately generate an auxetic structure. The effective orthotropy of the structure is ensured by the symmetry constraint imposed on the design variables. The sensitivities of above mentioned responses are obtained by means of the adjoint method, similar to Silva et al. (1997). The optimal solution is sought after by means of the Method of Moving Asymptotes (MMA) (Svanberg 1987).

3.3 Shape optimization

Following the initial generation of the structure layout using topology optimization, we utilize the shape generation tool as described in the beginning of the chapter and perform the fine tuning of the design by means of the shape optimization in the EDD setting. The EDD approach to shape optimization is outlined in detail in Riehl and Steinmann (2017) for linear elastic problems and further explored for simultaneous topology and shape optimization in Stankiewicz et al. (2021) and sequential topology and shape optimization of compliant mechanisms in Dev et al. (2022) in linear elasticity with boundary stress constraints and in Stankiewicz et al. (2022) in nonlinear elasticity with the addition of local adaptive refinement strategies and curvature constraint.

In Fig. 3, the adaptive refinement procedure for the EDD state problem in the last iteration of a sequentially optimized piezoelectric metamaterial is shown. The finite elements of the embedding domain are subdivided into the inner elements (coloured in red), boundary elements (cut by the shape, coloured in green) and the outer elements (coloured in white). A close-up of the locally refined region based on the curvature criterion is additionally shown in Fig. 4. For a detailed description on the adaptive refinement strategy for EDD, we refer to Stankiewicz et al. (2022).

Fig. 3

Uniform and local adaptive refinement steps for EDD state problem. a–c depict a uniform refinement of all boundary elements of the embedding domain (coloured in green). d shows local refinement of the boundary elements based on the curvature criterion for locally improved accuracy. (Color figure online)

Fig. 4

Close up of the locally refined base and shape in the area with high curvature of the boundary

Usage of the adaptive refinement procedure for the domain and the shape guarantees accuracy of the solution and improves computational efficiency and the resolution of the shape curvature. In Fig. 5, three different discretizations are depicted that represent the design between the topology and shape optimization steps. In Table 1, the corresponding effective properties of the structure computed using each of the discretizations are given. The differences between the material coefficients reach even 45% for the case of $C^{22}_{\textrm{H}}$. A number of factors contribute to such a large divergence in the predicted material values. On the one hand, we deal with the inherent limitations of density-based topology optimization, which are the presence of grey cells that require interpolation of the material tensor, and the staircase-like boundary. On the other hand, the shape extraction procedure is purely heuristic and does not guarantee a feasible design right away. The authors are aware that the precision of the topology optimization state problem can be improved by refining the mesh (at the cost of computational efficiency). However, specifically for the sequential approach, the current methodology is sufficient, as the shape optimization step guarantees an accurate solution to the state problem (as compared with the standard FEM) and efficient computation due to the adaptive refinement strategies.

Fig. 5

Solution comparison of the state problem within the sequential optimization framework. a The density distribution at the final, converged iteration of the topology optimization step. b An EDD problem and c a standard FEM model, both based on the directly extracted shape from the topology optimization result in a

Table 1 Selected effective properties computed after the topology optimization step

The problem definition for the shape optimization step is analogous to that of the topology optimization step, with the exception of using just $d^h_{\textrm{H}}$ as the objective function and the addition of the curvature constraint

$$\begin{aligned} \begin{aligned} {\mathcal {F}}\left( {\varvec{x}, \varvec{u}} \right) =&\ d^h_{\textrm{H}} \Rightarrow \ \ \text {Max} \\ \text {st.} \ \ \ {\mathcal {G}}_{\textrm{vol}} =&\int _{{\mathcal {B}}}\ \text {\,dV}- \overline{V} \le \ 0 \\ C_{\textrm{H}}^\textrm{ij} \ge&\ \bar{C}^\textrm{ij} \ \ \ \text {for \ \ i,j} \in \{ 1,2,3 \}\\ \kappa \le&\ {\bar{\kappa }} \end{aligned} \end{aligned}$$

(41)

where the objective function is a function of the nodal positions of the boundary and $\kappa$ stands for the curvature magnitude of the shape. The curvature magnitude is defined as the second derivative of the arclength parameterized plane curve $\varvec{x} (s)$ as follows:

$$\begin{aligned} \kappa = \partial _s^2 \varvec{x} \cdot \varvec{N} \end{aligned}$$

(42)

The discretized version of $\kappa$ is obtained by employing Taylor series expansion to approximate $\partial _s^2 \varvec{x}$, which for a node i, after some modifications, is given by

$$\begin{aligned} \kappa ^i = \frac{\varvec{N}^i \cdot \left[ \varvec{x}^{i+1} - \varvec{x}^i \right] }{ \Vert \varvec{x}^{i+1} - \varvec{x}^i \Vert ^2 } + \frac{\varvec{N}^i \cdot \left[ \varvec{x}^{i-1} - \varvec{x}^i \right] }{ \Vert \varvec{x}^{i-1} - \varvec{x}^i \Vert ^2 } \end{aligned}$$

(43)

where $\varvec{N}^i$ is the nodal normal vector calculated as a weighted average of the adjacent segment normals. For more details on the derivation of the curvature constraint, see Stankiewicz et al. (2022). Note that the implementation of the curvature constraint is integrated within the regularization procedure of the shape sensitivities. The regularization of the sensitivities, the bounding box as well as the curvature constraint are all realized via the traction method, which was first introduced in Azegami and Takeuchi (2006) and adapted to the EDD setting with the inclusion of the aforementioned geometrical constraints in Stankiewicz et al. (2022).

4 Numerical studies

In the following section, we conduct numerical experiments that evaluate the influence of various optimization setups and discuss the obtained results. First of all, we consider a well-known re-entrant honeycomb (see Fig. 6) and perform pure shape optimization to investigate its influence on an already established design and its dependency on various sets of constraints for the effective elastic tensor. Next, we conduct three studies of a complete design approach, i.e. sequential topology and shape optimization. In the first study, we compare the two different objectives for the topology optimization step, as given in the previous section. In the second study, we vary the values of the elastic coefficients and in the third study, we compare the possible symmetry conditions. Moreover, in each of the studies, we consider the influence of the constraint on the effective shear elastic coefficient.

4.1 Bulk properties

The bulk material used for the case studies in this work is barium titanate (BTO, chemical formula: $\mathrm BaTiO_3$). BTO is a lead-free alternative to the commonly used lead zirconate titanate (PZT), which finds application in piezoelectric energy harvesters. As typical for a ceramic material, BTO is highly brittle and stiff, which limits its application in terms of operational frequencies and deformation under vibrational energy harvesting or sensing. Hence, the application of the metamaterial concept to BTO offers great possibilities for improved deformability, larger operational frequency range and lightweight design. Since BTO cannot undergo large strains and for computational efficiency, only a geometrically linear model for the state problem is considered. In Appendix 1 all necessary material coefficients are given.

4.2 Shape optimization of engineered structures

The re-entrant honeycomb is the standard layout when it comes to auxetic structures. An exemplary design is shown in Fig. 6, and the effective material properties of which are the following

Fig. 6

Re-entrant honeycomb design of an auxetic structure

$$\begin{aligned}{} & {} d^\textrm{h}_{\textrm{H}} = 177 \frac{\textrm{pC}}{\textrm{N}} \end{aligned}$$

(44)

$$\begin{aligned}{} & {} \varvec{C}_{\textrm{H}} = \begin{bmatrix} 3.55 &{} -4.47 &{} 0\\ -4.47 &{} 8.54 &{} 0\\ 0 &{} 0 &{} 0.11 \end{bmatrix} {\mathrm{GPa}} \end{aligned}$$

(45)

$$\begin{aligned}{} & {} \varvec{d}_{\textrm{H}} = \begin{bmatrix} 0 &{} 0 &{} 196\\ -13 &{} 190 &{} 0 \end{bmatrix} \frac{\textrm{pC}}{\textrm{N}} \end{aligned}$$

(46)

The $d^\textrm{h}_{\textrm{H}}$ value is significantly superior as compared to the bulk material ($d^\textrm{h} = 112 \mathrm{pC/N}$). The normal stiffness in vertical direction ($C^{22}_{\textrm{H}}$, direction Y) is more than twice higher than in the horizontal direction ($C^{11}_{\textrm{H}}$, direction X). Moreover, the shear stiffness $C^{33}_{\textrm{H}}$ is relatively small which signifies a potential instability under imperfect loading. In the following, we investigated two different target stiffness configurations by means of pure shape optimization. The problem setup, i.e. the objective and the constraints, is outlined in Table 2.¹ The optimization problems are each limited to 200 iterations of the Multiplier Method (Hestenes 1969; Rockafellar 1973). Therefore, we accept slight violations of the constraints as the optimization problems often do not reach the assumed tolerances.

Table 2 Shape optimization setup

We directly maximize the effective hydrostatic piezoelectric coupling coefficient $d^h_{\textrm{H}}$ under the volume restriction of 40% of the unit cell. For the stiffness constraints, in case 1, we aim to obtain a structure with similar stiffness in both X and Y directions. In case 2, however, we try to preserve similar stiffness values as in the initial design. In both cases, we introduce the curvature constraint $\kappa \le 5 \left[ \mathrm{1/mm}\right]$ which is equivalent to a minimum radius constraint of $r > 0.2 \left[ {\text {mm}}\right]$ to avoid stress concentration.

In Fig. 7, the resultant structures from the shape optimization of the re-entrant honeycomb are shown with the effective properties listed in Table 3.

Table 3 Selected effective properties of the resultant structures of the shape optimization

To begin with, an interesting observation is that the objectives $d^\textrm{h}_{\textrm{H}}$ and $\upsilon _{\textrm{H}}$ lead to strongly similar designs, see, for instance, the cases 1, 3 and 2, 4, pairwise. Moreover, the improvement of $d^\textrm{h}_{\textrm{H}}$ is negligible and the value of $\upsilon _{\textrm{H}}$ is indeed better for the cases with lower $C^{22}_{\textrm{H}}$ stiffness. For the cases 1 and 3, the value of $C^{22}_{\textrm{H}} = 4$ has not been reached. We observe a tendency; however, the decrease of $C^{22}_{\textrm{H}}$ is realized via a change in the angle of the diagonal member in the structure. Nevertheless, the design changes during shape optimization alone are in general very small, which heavily restricts possible setups of the optimization problems. Still, we draw the main conclusion of the first study that the design of the re-entrant honeycomb is mainly a question of the desired stiffness proportion between $C^{11}_{\textrm{H}}$ and $C^{22}_{\textrm{H}}$ rather than finding the maximal $d^\textrm{h}_{\textrm{H}}$ value.

Fig. 7

Close-up of the resulting structures of the shape optimization problems from Table 2. The black design is the initial configuration

4.3 Complete design approach

In the following, we explore the design possibilities of piezoelectric metamaterials by employing the sequential topology and shape optimization routine. We analyse three case studies, each time choosing different constraints to assess their influence on the design. We start with the objective function, for which the effective Poisson’s ratio $\upsilon _{\textrm{H}}$ and effective hydrostatic piezoelectric coupling coefficient $d_{\textrm{H}}^\textrm{h}$ are compared in the initial design generation phase by topology optimization. In conjunction with this, the presence of a constraint on the shear stiffness coefficient $C^{33}_{\textrm{H}}$ is studied. In Table 4, the optimization setup is shown. The responses of interest are additionally highlighted.

Table 4 Problems setup for the comparison of different objective functions and for the presence of $C^{33}_{\textrm{H}}$ constraint

Note that the objective for the shape optimization for all cases is $d^\textrm{h}_{\textrm{H}}$, since this is the ultimate design goal. Regarding the stiffness characteristics, we require the same normal stiffness in both X and Y directions. The obtained structures together with the chosen effective properties are depicted in Table 5.

Table 5 Resulting structures and their effective properties for the optimization setup in Table 4

It is immediately apparent that significantly varying geometries were obtained for each setup, although the objective values remain practically the same. In particular, the case 5 generated the well-known re-entrant honeycomb structure. Nevertheless, all of the designs exhibit auxetic behaviour with Poisson’s ratios of $\upsilon _{\textrm{H}} \approx -\,0.8$ or lower, which confirms the unanimity of both objective functions. The structures from the upper row in Table 5 exhibit very low shear stiffness coefficient $C_{\textrm{H}}^{33}$ in relation to $C_{\textrm{H}}^{11}$ and $C_{\textrm{H}}^{22}$. By introducing the constraint $C_{\textrm{H}}^{33} > 0.50$ , visibly different designs were obtained without affecting the objective values. A notable tendency observed is a shortened diagonal member, which is most visible in case 7. The effective elastic and piezoelectric tensors for the example of case 7 are given as follows:

$$\begin{aligned}{} & {} \varvec{C}_{\textrm{H}}^{(5)} = \begin{bmatrix} 4.66 &{}\quad -4.16 &{}\quad 0\\ -4.16 &{}\quad 4.31 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 0.5 \end{bmatrix} \end{aligned}$$

(47)

$$\begin{aligned}{} & {} \varvec{d}_{\textrm{H}}^{(5)} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad 200\\ -11 &{}\quad 191 &{}\quad 0 \end{bmatrix} \end{aligned}$$

(48)

which proves the feasibility of the structure with the imposed constraints on the elastic coefficients. However, we notice that the stiffness coefficients did not reach the lower limits although the structure formed very narrow hinges. This observation once again points our attention at the difference in accuracy between the state problems in the topology and shape optimization steps. Thus, taking into consideration manufacturability of the structures, an adjusted set of stiffness constraints for the shape optimization step appears to be a reasonable choice.

The second study is devoted to the influence of stiffness requirements in general, i.e. the constraints on all diagonal elastic coefficients $C_{\textrm{H}}^\textrm{ii}$ were varied, with the requirement that the normal stiffness coefficients still remain similar $C_{\textrm{H}}^{11} \approx C_{\textrm{H}}^{22}$. A summary of the problems setup for this study is given in Table 6.

Table 6 Problems setup for the comparison of different stiffness constraints on $C^{11}_{\textrm{H}}$ and $C^{22}_{\textrm{H}}$ and for the presence of $C^{33}_{\textrm{H}}$ constraint

For each case, the objective is to maximize the hydrostatic piezoelectric coupling coefficient $d^\textrm{h}_{\textrm{H}}$, as the previous study has shown that this formulation is capable of generating auxetic structures in the topology optimization step. The cases 6 and 8 are taken from the previous study, hence the unordered numbering. Two new cases (9 and 10) consider higher stiffness requirements, whereas for the case 10, the constraint on $C^{33}_{\textrm{H}}$ is increased twice as compared to the case 8. In Table 7, the resultant designs and their selected effective properties are given.

Table 7 Resulting structures and their effective properties for the optimization setup in Table 6

Besides the obvious geometrical differences, we observe a major decrease (in the absolute sense) of the effective Poisson’s ratio $\upsilon _{\textrm{H}}$ for the stiffer structures (case 9 and 10). This highlights a trade-off between increased stiffness and the ability of a structure to exhibit auxetic deformation, although this did not affect the value of $d^\textrm{h}_{\textrm{H}}$ significantly. The addition of the $C_{\textrm{H}}^{33}$ (case 10) further drops the Poisson’s ratio to $\upsilon _{\textrm{H}} = -\,0.53$. The effective elastic and piezoelectric material tensors for the example of case 10 read

$$\begin{aligned} \varvec{C}_{\textrm{H}}^{(10)} = \begin{bmatrix} 6.29 &{}\quad -\,3.28 &{}\quad 0\\ -\,3.28 &{}\quad 6.03 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 1.02 \end{bmatrix} \end{aligned}$$

(49)

$$\begin{aligned} \varvec{d}_{\textrm{H}}^{(10)} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad 202\\ -\,15 &{}\quad 190 &{}\quad 0 \end{bmatrix} \end{aligned}$$

(50)

where the value of $C_{\textrm{H}}^{11}$ is slightly greater than the lower limit of 6 GPa. Nevertheless, the constraints were not violated which confirms the feasibility of the obtained design.

The third study concerns the possible symmetry conditions, which, for two-dimensional problems, besides the double XY-symmetry used for all previous studies, it includes single X- and Y-symmetries. The problem setups are listed in Table 8, where for all the symmetry cases, a study with and without the presence of the $C_{\textrm{H}}^{33}$ constraint was conducted, similarly to the previous comparisons. The XY-symmetry cases are the same as in Table 7, hence the unordered numbering of the cases.

Table 8 Problems setup for the comparison of different symmetry conditions and for the presence of $C^{33}_{\textrm{H}}$ constraint

The presence of the symmetry condition in general is necessary to guarantee orthotropy of the material, without the need for an explicit orthotropy constraint. The resulting structures and their effective properties are given in Table 9.

Table 9 Resulting structures and their effective properties for the optimization setup in Table 8

The structures with the single symmetry condition are notably different, since a larger design space is available for the optimizer. The case 12 (Y-symmetry) generated an infeasible design with a close to zero value of $C^{33}_{\textrm{H}}$. In this particular case, the addition of the constraint on $C^{33}_{\textrm{H}}$ (see case 14) alleviated the issue and a feasible design was generated. Both of the designs with X-symmetry (case 11 and 13) reached relatively high values of $C^{33}_{\textrm{H}}$; hence, there is practically no difference between them, except the opposite orientation. Similarly to the previous comparisons, the effective elastic and piezoelectric tensors are given for case 14

$$\begin{aligned}{} & {} \varvec{C}_{\textrm{H}}^{(14)} = \begin{bmatrix} 6.76 &{}\quad -\,3.27 &{}\quad 0\\ -\,3.27 &{}\quad 6.05 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 1.17 \end{bmatrix} \end{aligned}$$

(51)

$$\begin{aligned}{} & {} \varvec{d}_{\textrm{H}}^{(14)} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad 201\\ -\,13 &{}\quad 190 &{}\quad 0 \end{bmatrix} \end{aligned}$$

(52)

where the values of $C^{11}_{\textrm{H}}$ and $C^{33}_{\textrm{H}}$ are slightly above their lower limits, not affecting, however, the feasibility of the design.

As previously observed, the transition between the topology and the shape optimization leads to a jump in the effective material coefficients due to the reduced accuracy of the state problem in the topology optimization step, as well as due to the heuristics behind the shape generation step. Therefore, if one roughly knows the discrepancy between the effective coefficients when transitioning from topology to shape optimization, a reasonable approach is to select a different set of constraint values for each of the steps. For instance, a set of constraint values can be reduced by a certain amount for the topology optimization step and set back to the desired values only for the shape optimization step. Thus, in Table 10, we introduced adjustments to the constraint values on the stiffness coefficients for the shape optimization step on the examples of cases 5 and 7.

Table 10 Problems setup for the comparison of adjusted constraints for the shape optimization step for the cases 5 and 7

The case 5 is particularly peculiar for two reasons. First, the largest jump in the value occurred for $C^{22}_{\textrm{H}}$, leading to the initial design for shape optimization yielding $C^{22}_{\textrm{H}}$ almost two times greater than $C^{11}_{\textrm{H}}$. Second, only an increase in the volume constraint from $V \le 40\%$ to $V \le 50\%$ facilitated the fulfilment of the constraint on $C^{11}_{\textrm{H}}$. In Table 11, the $C^{11}_{\textrm{H}}$ and $C^{22}_{\textrm{H}}$ stiffness coefficients for the cases 5a and 5b are shown at different stages of the sequential optimization routine.

Table 11 $C^{11}_{\textrm{H}}$ and $C^{22}_{\textrm{H}}$ stiffness values at different stages of the sequential optimization for the cases 5a and 5b

Table 12 Resulting structures and their effective properties for the optimization setup in Table 10

In the shape optimization phase, an attempt to retrieve the similarity of $C^{11}_{\textrm{H}}$ and $C^{22}_{\textrm{H}}$ stiffness coefficients was made by choosing the constraint value of $C^\textrm{ii}_{\textrm{H}} \ge 5$ for both the coefficients, which turned out not be reachable for the case 5a. Therefore, the increase in the volume constraint enabled design changes effective enough to retrieve similar values of $C^{11}_{\textrm{H}}$ and $C^{22}_{\textrm{H}}$, both of them fulfiling the optimization constraints.

The resulting structures for the adjusted sets of constraints for the cases 5 and 7 are shown in Table 12. For the case 5, one can immediately see that an increase of the $C^\textrm{ii}_{\textrm{H}}$ constraints by $25\%$ leads to a reduction of the length of the diagonal members (case 5b), mostly because of the difficulty to fulfil the constraint on $C^{11}_{\textrm{H}}$. Significantly greater stiffness coefficients could be imposed for the case 7, which shows very thin hinges in the original setup. Here, an increase of up to $75\%$ of the stiffness coefficients’ constraints lead to a generation of hinge-free structures with greatly improved manufacturability.

In Fig. 8, a visualization of three metamaterials based on the cases 7, 10 and 14 is shown (for which the effective elastic and piezoelectric tensors were given as well). This rendering helps to notice the similarities between the designs obtained in this work. In particular, the topological equivalence between all the cases is evident in the stiffness study, the results of which are depicted in Table 7.

Fig. 8

Repeated structures of the metamaterials

4.4 Role of fine tuning with shape optimization

Finally, the evolution of a sequential optimization problem is shown in two examples to illustrate the importance of the shape optimization step. In Fig. 9, the intermediate results of the sequential optimization of cases 7(b) and 10 are shown. These two examples were chosen to show two different aspects that are common when sequential optimization is applied to the design of piezoelectric microstructures.

Fig. 9

Sequential optimization of the cases 7 and 10. From the left side: topology optimization results, generated shapes, shape optimization results. The case 7 depicts two final shapes: with unchanged constraint values (case 7, green arrow) and with modified constraint values (case 7b, orange arrow). (Color figure online)

Case 7 is presented to illustrate the accuracy issue due to material interpolation within topology optimization. For this reason, the difference between the member thickness after the shape generation and after the shape optimization is particularly noticeable. This is due to the large jump in stiffness values between the topology optimization result and the generated shape. The topology optimization result contains a relatively high number of cells with intermediate density ($23\%$ cells with pseudo-density in the range $\rho \in \left[ 0.01, 0.99\right]$) after 500 iterations using the method of moving asymptotes (MMA). Therefore, the predicted material coefficients are not accurate and the resulting structure is not optimal with respect to the problem setup. As a result, an excessively stiff shape was obtained immediately after its generation. Consequently, the final structure of case 7 contains very thin features that are not suitable for manufacturing. Adjusting the stiffness constraints, as seen in case 7b, helped to obtain a manufacturable design.

Case 10 shows a large design evolution within the shape optimization alone, again caused by the jump in computed material properties after the shape generation. The shape optimization step finally approaches the stiffness constraints after a significant shortening of the diagonal members.

In Figs. 10 and 11, the evolution of the hydrostatic coupling coefficient $d^\textrm{h}_{\textrm{H}}$ and the stiffness coefficients $C^\textrm{ii}_{\textrm{H}}$, respectively, is plotted for cases 10 and 14. The iteration count encompasses the complete sequential procedure, i.e. the shape optimization iterations are counted starting from the last topology optimization iteration. In the $d^\textrm{h}_{\textrm{H}}$ plot, one can immediately notice that the value is strongly overestimated in the early iterations due to the presence of the grey cells. Hence, the evolution of $d^\textrm{h}_{\textrm{H}}$ is difficult to analyse until the design is completely black and white. However, one can clearly distinguish the transition from topology to shape optimization by the sudden drop of $d^\textrm{h}_{\textrm{H}}$ in the late stage of optimization. In case 10, it occurred at about iteration 500 and in case 14 at iteration 300. The sudden drop of $d^\textrm{h}_{\textrm{H}}$ indicates the error that even a converged topology optimization introduces in the computation of effective material parameters. This error is even more pronounced in the $C^\textrm{ii}_{\textrm{H}}$ plot, where two spikes are shown at iteration 300 and 500, respectively. The normal stiffness coefficient $C^\textrm{11}_{\textrm{H}}$ was constrained from both the bottom (6 GPa) and the top (7 MPa) to avoid excessive stiffness that occurred in other design attempts. The jump in value of up to $70\%$ is visible for $C^\textrm{11}_{\textrm{H}}$ in case 14 after shape generation, emphasizing the importance of design adaptation with shape optimization.

Fig. 10

Evolution of the hydrostatic coupling coefficient over the course of sequential optimization for cases 10 and 14

Fig. 11

Evolution of the constrained stiffness coefficients over the course of sequential optimization for cases 10 and 14

4.5 Discussion

Negligible difference between the values of the effective hydrostatic piezoelectric coupling coefficient $d^\textrm{h}_{\textrm{H}}$, which for the cases considered in this work varies merely within the range of $d^\textrm{h}_{\textrm{H}} \in \left[ 175, 181\right]$ (the value of $d^\textrm{h}_{\textrm{H}} = 185$ for the case 12 is ignored due to the infeasibility of the design) suggests that decision about a suitable design of the piezoelectric metamaterial for hydrostatic pressure applications is not a question of maximized performance, but rather the desired stiffness characteristics. Depending on the applied loading conditions, optimal performance can be achieved by designing a structure, which deformation does not exceed the stress limits and also is compliant enough to provide a sufficient energy conversion., e.g. for energy harvesting applications.

The final design also largely depends on the requirements for the $C_{\textrm{H}}^{33}$ values, since the well-known re-entrant honeycomb structures exhibit disproportionately low value of $C_{\textrm{H}}^{33}$ as compared to $C_{\textrm{H}}^{11}$ and $C_{\textrm{H}}^{22}$, thus imposition of a $C_{\textrm{H}}^{33}$ constraint significantly alters the design, providing structures more resistant to a potentially unwanted deformation when a metamaterial with large number of repeated cells is considered. Moreover, increased $C_{\textrm{H}}^{33}$ contributes to a decrease of the absolute value of the Poisson’s ratio, without affecting, however, the value of $d^\textrm{h}_{\textrm{H}}$.

Introduction of the symmetry conditions is an efficient way to guarantee orthotropy of the resultant metamaterial, simultaneously limiting the design space and hence improving the convergence to a feasible design. By imposing a double symmetry for X and Y planes, one can obtain a design equivalent to the re-entrant honeycomb (case 5).

The discrepancy between the computed material coefficients in the topology and shape optimization steps is tackled by selecting lower stiffness constraints for the topology optimization step and then by increasing them in the shape optimization step. This way, formation of narrow hinges is avoided in the final structure.

5 Conclusions

In this contribution, we applied the sequential topology and shape optimization framework (Dev et al. 2022) to design auxetic piezoelectric metamaterials made of (lead-free) barium titanate ceramic. By combining topology and shape optimization, we benefited from the advantages of each method. In the sequential framework, first topology optimization was utilized to generate the design layouts, in which the complexity of the designs was controlled by the regularization scheme: density filtering, projection and symmetry constraint. Next, by using a simple shape generation algorithm, we performed shape optimization to fine tune the designs. By exploiting the shape optimization, we profited from: (1) accurate and efficient computations of the state problem by the EDD method with adaptive base and shape refinement; and (2) presence of crisp boundary, which enabled consideration of the curvature constraint to control manufacturability and durability of the designed structures.

The case studies have shown that the design of a piezoelectric metamaterial for hydrostatic pressure applications mainly depends on the desired stiffness characteristics. By varying the constraint values for the elastic coefficients and symmetry conditions significantly varying structures were obtained. Minimization of the effective Poisson’s ratio and maximization of the hydrostatic piezoelectric coupling coefficient both lead to generation of auxetic structures. By employing the sequential optimization framework, we obtained exact designs with accurately computed effective material properties.

5.1 Outlook

Usage of brittle ceramic material motivated us to employ a geometrically linear model to compute the effective properties. However, for less brittle materials, like, for instance, PVDF, a geometrically nonlinear model would be beneficial to design structures with optimal material properties under larger deformations (despite the increased computational cost).

The sequential topology and shape optimization framework can be extended towards the design of metamaterials based on other applications or physical models, e.g. photonic metamaterials.

Furthermore, experimental investigation is planned in the future work to verify the applicability of the proposed optimization procedure and the obtained material properties. Although experimental studies are not within the scope of this contribution, in Fig. 12, we show test manufacturing of selected samples from the cases 5,7 and 10.

Fig. 12

Manufactured samples of the geometries from the cases 5,7 and 10. The samples are shown before and after the sintering

The selected samples were fabricated experimentally by ceramic injection moulding of barium titanate (BTO). The basic process for fabricating complex geometric structures has been described in detail in previous work (Biggemann et al. 2018; Köllner et al. 2023; Hoffmann et al. 2023). However, an adaptation of the process flow for complex structures was required. In order to optimize the injection moulding process, the positive moulds of structures were designed using OpenSCAD (Kintel and Wolf 2014) from the structures exported as STL. These CAD files were stereolithographically printed with an Anycubic Photon Mono 4K (Anycubic, Shenzhen, China) using Anycubic Translucent UV Resin (Anycubic, Shenzhen, China) with a z-resolution of 50 $\upmu$m. The printed templates were then moulded with a two-part RTV-2 silicone rubber (ELASTOSIL®M 4514 Wacker Chemie AG, Munich, Germany) to obtain a negative mould. Transfer moulding was performed with a ceramic feedstock containing $55\, \mathrm vol\%$ BTO powder (Sigma Aldrich, Darmstadt, Germany), $40\, \mathrm vol\%$ paraffin wax (Granopent®P, Carl Roth GmbH, Karlsruhe, Germany) and $5\, \mathrm vol\%$ carnauba wax (Naturfarben, Carl Roth GmbH, Karlsruhe, Germany). Transfer moulding was carried out at 120 °C under vacuum (< 10 Pa). One of the difficulties in comparison with previous structures was that the areas of the struts in the centre of the structure and in the head and base areas were very thin. Typically, a minimum width of 6–8 mm is required to ensure non-destructive moulding from the silicone moulds. The debinding process, performed at 650 °C for 1 h in an Al2O3 crucible (CERAMTRADE, Bad Vilbel, Germany) with an $\mathrm ZrO_2$ (Tosoh USA, Inc., Grove City, USA) powder bed, and the subsequent sintering process, performed at 1400 °C for 6 h, were performed on $\mathrm ZrO_2$ substrates placed on porous mullite substrates (Annamullit®88, Compagnie de Saint-Gobain S.A., Courbevoie, France). Optimized heating and cooling rates of 0.1 to 5 K/min were used in these procedures.

Successful injection moulding and sintering of the selected structures prove that the geometries are manufacturable directly from the results of the sequential optimization procedure despite their geometric complexity.

Appendix 1: Bulk material properties: \(\textrm{BaTiO}_3\)

For all the examples studied in this work, we utilize BTO material properties, given as follows: Material density:

$$\begin{aligned} \rho _{mat} = 5700 \frac{{\textrm{kg}}}{{\textrm{m}}^{3}} \end{aligned}$$

(53)

Compliance tensor:

$$\begin{aligned} \varvec{S} = \begin{bmatrix} 9.1 &{}\quad -\,2.9 &{}\quad 0\\ -\,2.9 &{}\quad 9.5 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 22.8 \end{bmatrix} \cdot 10^{-12} \frac{{\textrm{m}}^{2}}{\textrm{N}} \end{aligned}$$

(54)

Piezoelectric tensor (strain-charge format):

$$\begin{aligned} \varvec{d} = \begin{bmatrix} 0 &{}\quad 0 &{}\quad 260\\ -\,78 &{}\quad 190 &{}\quad 0 \end{bmatrix} \cdot 10^{-12} \frac{{\textrm{C}}}{\textrm{N}} \end{aligned}$$

(55)

Relative permittivity tensor (strain-charge format):

$$\begin{aligned} \frac{\varvec{\epsilon }^T}{\epsilon _0} = \begin{bmatrix} 1450 &{}\quad 0\\ 0 &{}\quad 1700 \end{bmatrix} \end{aligned}$$

(56)

Elasticity tensor (inverse of the compliance tensor):

$$\begin{aligned} \varvec{C} = \begin{bmatrix} 122 &{}\quad 37 &{}\quad 0\\ 37 &{}\quad 117 &{}\quad 0\\ 0 &{}\quad 0 &{}\quad 44 \end{bmatrix} \cdot 10^{9} \frac{{\textrm{N}}}{{\textrm{m}}^{2}} \end{aligned}$$

(57)