Abbas Homayouni-Amlashi

Thomas Schlinquer

Abdenbi Mohand-Ousaid

Micky Rakotondrabe

2D topology optimization MATLAB codes for piezoelectric actuators and energy harvesters

Keywords: Topology optimization, MATLAB code, Piezoelectric actuator, Piezoelectric energy harvester

In this paper, two separate topology optimization MATLAB codes are proposed for a piezoelectric plate in actuation and energy harvesting. The codes are written for one-layer piezoelectric plate based on 2D finite element modeling. As such, all forces and displacements are confined in the plane of the piezoelectric plate. For the material interpolation scheme, the extension of solid isotropic material with penalization approach known as PEMAP-P (piezoelectric material with penalization and polarization) which considers the density and polarization direction as optimization variables is employed. The optimality criteria and method of moving asymptotes (MMA) are utilized as optimization algorithms to update the optimization variables in each iteration. To reduce the numerical instabilities during optimization iterations, finite element equations are normalized. The efficiencies of the codes are illustrated numerically by illustrating some basic examples of actuation and energy harvesting. It is straightforward to extend the codes for various problem formulations in actuation, energy harvesting and sensing. The finite element modeling, problem formulation and MATLAB codes are explained in detail to make them appropriate for newcomers and researchers in the field of topology optimization of piezoelectric material.

Introduction

Topology optimization (TO) is a methodology to distribute the material within a design domain in an optimal way while there is no prior knowledge about the final layout of the material [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF]. This main specification of TO provides a great degree of freedom in terms of designing innovative structures to satisfy predefined engineering goals. Historically, minimization of mechanical deformation of a structure under application of different loading conditions was a classical engineering goal [START_REF] Schmit | Structural design by systematic synthesis[END_REF]. Aiming for this goal, the work of [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF] paves the way for a methodology known today as topology optimization. The general idea of this methodology is the combination of finite element method and optimization to maximize or minimize an objective function. In this regard, the design domain is discretized by a finite number of elements and design variables for the optimization problem are the variables attributed to each of these elements. Different approaches are introduced in the literature to implement the TO method [START_REF] Sigmund | Topology optimization approaches[END_REF]). Among these approaches, homogenization approach was proposed to optimize the porous elements as unit cells or microstructures within the design domain to obtain the final layout [START_REF] Suzuki | The method of moving asymptotes-a new method for structural optimization[END_REF]Bendsøe and Sigmund 1995). However, the number of optimization variables is high in homogenization method which can make the optimization cumbersome. The other famous and popular approach is the SIMP approach which stands for solid isotropic material with penalization. In this approach, the elements in the design domain can have intermediate densities [START_REF] Bendsøe | Optimal shape design as a material distribution problem[END_REF][START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF]). This will let the elements to be gray in addition to black (material) and white (void). However, due to practical constraints, it is desired that the optimization finally converges to a black and white layout. To do so, a penalization factor is defined for intermediate densities. One of the reasons for the popularity of this approach is its simplicity of implementation in comparison with other approaches. Similar to SIMP approach, there is evolutionary structural optimization (ESO) [START_REF] Xie | A simple evolutionary procedure for structural optimization[END_REF] or the more general form bi-directional evolutionary structural optimization (BESO) [START_REF] Xia | Bi-directional evolutionary structural optimization on advanced structures and materials: a comprehensive review[END_REF]) which is about removing or adding the elements inside the design domain during optimization iterations. In addition to the aforementioned approaches, there are other approaches including level set method [START_REF] Van Dijk | Levelset methods for structural topology optimization: a review[END_REF][START_REF] Andreasen | Level set topology and shape optimization by density methods using cut elements with length scale control[END_REF]) and method of moving morphable components (MMC) [START_REF] Guo | Doing topology optimization explicitly and geometrically-a new moving morphable components based framework[END_REF][START_REF] Yaghmaei | Explicit three dimensional topology optimization via moving morphable void (mmv) approach[END_REF]) which is a geometrical approach. For a detailed review and comparison between these approaches, one can refer to the following review papers: [START_REF] Sigmund | Topology optimization approaches[END_REF] and [START_REF] Cook | Topology optimization of bimorph piezoelectric energy harvesters considering variable electrode location[END_REF].

Due to the success of TO methodology, several implementation codes in different software are published in the literature. Using the SIMP approach, [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF] published the 99 lines of MATLAB code for 2D topology optimization of compliance problems. [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF] published the 88 lines of MATLAB code which was an improvement of Sigmund's 99 lines of code while having much faster speed in each iteration thanks to introducing the connectivity matrix that facilitates the assembly procedure of elemental matrices. [START_REF] Liu | An efficient 3d topology optimization code written in matlab[END_REF] published the TOP3D MATLAB code by extension of the 88 lines of code for topology optimization of 3D structures. [START_REF] Chen | A 213-line topology optimization code for geometrically nonlinear structures[END_REF] published 213 lines of MATLAB code for 2D topology optimization of geometrically nonlinear structures. There are other published codes using other TO approaches like level set method [START_REF] Challis | A discrete level-set topology optimization code written in matlab[END_REF][START_REF] Wei | An 88-line matlab code for the parameterized level set method based topology optimization using radial basis functions[END_REF][START_REF] Yaghmaei | Explicit three dimensional topology optimization via moving morphable void (mmv) approach[END_REF], BESO [START_REF] Xia | Bi-directional evolutionary structural optimization on advanced structures and materials: a comprehensive review[END_REF], and projection method [START_REF] Smith | A matlab code for topology optimization using the geometry projection method[END_REF]. These published codes facilitate the implementation of TO methodology for various applications. For this reason, the application of TO can be seen in solving different problems including the compliance problems [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF], compliant mechanism problems [START_REF] Zhu | Design of compliant mechanisms using continuum topology optimization: a review[END_REF], heat conduction [START_REF] Gersborg-Hansen | Topology optimization of heat conduction problems using the finite volume method[END_REF], and smart materials in particular the piezoelectric materials [START_REF] Sigmund | Design of smart composite materials using topology optimization[END_REF].

Due to their electromechanical coupling effect, piezoelectric materials have applications in actuation, sensing and energy harvesting. Plenty of methods can be found in the literature to analyze and improve the performance of the piezoelectric structures, whether actuators, energy harvesters or sensors such as geometrical and size opti-mization [START_REF] Schlinquer | Optimal design of a unimorph piezoelectric cantilever devoted to energy harvesting to supply animal tracking devices[END_REF][START_REF] Liseli | Optimal design of piezoelectric cantilevered actuators for charge-based self-sensing application[END_REF]Homayouni-Amlashi et al. 2020a), shape optimization [START_REF] Muthalif | Optimal piezoelectric beam shape for single and broadband vibration energy harvesting: Modeling, simulation and experimental results[END_REF], layers number optimization [START_REF] Rabenorosoa | Performances analysis of piezoelectric cantilever based energy harvester devoted to mesoscale intra-body robot Rakotondrabe M (2011) Performances inclusion for stable interval systems[END_REF], or parameters sub-optimization [START_REF] Rakotondrabe | Design of piezoelectric actuators with guaranteed performances using the performances inclusion theorem and interval tools[END_REF][START_REF] Khadraoui | Optimal design of piezoelectric cantilevered actuators with guaranteed performances by using interval techniques[END_REF] with interval techniques (Rakotondrabe 2011). After development of TO methodology, it is extended to different physics [START_REF] Alexandersen | A review of topology optimisation for fluid-based problems[END_REF][START_REF] Cook | Topology optimization of bimorph piezoelectric energy harvesters considering variable electrode location[END_REF] including the piezoelectricity. Primarily, the homogenization approach is used [START_REF] Silva | Optimal design of piezoelectric microstructures[END_REF][START_REF] Sigmund | On the design of 1-3 piezocomposites using topology optimization[END_REF]. Afterwards, other approaches including SIMP [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF], BESO (de Almeida 2019) or level set method [START_REF] Chen | A level set approach for optimal design of smart energy harvesters[END_REF]) are also explored. By defining proper objective functions, TO methodology is applied to piezoelectric actuators [START_REF] Moretti | Topology optimization of piezoelectric bimaterial actuators with velocity feedback control[END_REF][START_REF] Gonċalves | Simultaneous optimization of piezoelectric actuator topology and polarization[END_REF], sensors [START_REF] Menuzzi | Piezoelectric sensor location by the observability gramian maximization using topology optimization[END_REF]) and energy harvesters (Homayouni-Amlashi et al. 2020b;[START_REF] Homayouni-Amlashi | Analytical modelling and optimization of a piezoelectric cantilever energy harvester with in-span attachment[END_REF][START_REF] Townsend | Topology optimization of vibrational piezoelectric energy harvesters for structural health monitoring applications[END_REF]. The publications considered different types of system modeling including the static [START_REF] Zheng | Topology optimization of energy harvesting devices using piezoelectric materials[END_REF], dynamic [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF][START_REF] Wein | Topology optimization of a piezoelectric-mechanical actuator with single-and multiple-frequency excitation[END_REF], modal [START_REF] Wang | Topology optimization of piezocomposite resonator for maximizing excitation strength and synthesizing desired eigenmodes[END_REF]) and electrical circuit coupling [START_REF] Salas | A topology optimization formulation for transient design of multi-entry laminated piezocomposite energy harvesting devices coupled with electrical circuit[END_REF][START_REF] Rupp | Design of piezoelectric energy harvesting systems: a topology optimization approach based on multilayer plates and shells[END_REF]. Different types of problem formulation can be found as well such as optimization with stress constraints [START_REF] Wein | Topology optimization of a cantilevered piezoelectric energy harvester using stress norm constraints[END_REF]. Although the application of TO methodology to piezoelectric materials is well established in the literature, no implementation code is published yet.

In this paper, 2D topology optimization codes are proposed for actuation and energy harvester by using the extension of SIMP approach known as PEMAP-P (piezoelectric material with penalization and polarization). The codes are written based on the 88 lines of MATLAB code written by [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF], except the code is extended and modified considerably to consider the electromechanical coupling effect of piezoelectric material and problem formulation. In Section 2, the finite element modeling of one-layer piezoelectric plate is presented by using the plane stress assumption. The finite element equation is derived for both actuation and energy harvesting. A normalization is applied to the finite element equation which significantly reduces the numerical instabilities in optimization iterations. Hence, the proposed codes in this paper work smoothly in both actuation and energy harvesting. In Section 3, first, the material interpolation scheme for piezoelectric material is explained. Then, the optimization problem is formulated for both actuation and energy harvesters. The problem formulations are basic for educational purposes. The sensitivity analysis is performed, and finally, the optimization algorithms are explained. In Section 4, the MATLAB codes are explained part by part in detail. The explanations in this part help the readers to extend and implement the codes for their own purposes. In Section 5, different numerical examples are illustrated and the modification to the original codes to implement those numerical examples are expressed as well.

Finite element modeling

Constitutive equation

The general linearly coupled mechanical and electrical constitutive equation of piezoelectric materials by neglecting the thermal coupling can be written as [START_REF] Lerch | Simulation of piezoelectric devices by two-and threedimensional finite elements[END_REF])

T = c E S -e Ē D = e T S + ε S Ē (1)
In ( 1), T and S are the vectors of mechanical stress and strain while c E is the stiffness tensor in constant electrical field. D and Ē are the vectors of electrical displacement and electrical field. e is the piezoelectric matrix, ε S is the matrix of permittivity in constant mechanical strain and T shows the matrix transpose.

The 4-mm tetragonal crystal class piezoelectric material [START_REF] Piefort | Finite element modelling of piezoelectric active structures[END_REF] which has orthotropic anisotropy is considered to derive the corresponding model. This class includes most of the piezoelectric material in particular the wellknown PZT (lead zirconate titanate) materials. By this consideration, the mechanical stiffness tensor, piezoelectric matrix and permittivity for full 3D modeling are

c E = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ c E 11 c E 12 c E 13 0 0 0 c E 12 c E 11 c E 13 0 0 0 c E 13 c E 13 c E 33 0 0 0 0 0 0 c E 44 0 0 0 0 0 0 c E 44 0 0 0 0 0 0 c E 66 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ e T = ⎡ ⎣
0 0 0 0 e 15 0 0 0 0 e 15 0 0 e 31 e 31 e 33 0 0 0

⎤ ⎦ ε S = ⎡ ⎣ ε S 11 0 0 0 ε S 11 0 0 0 ε S 33 ⎤ ⎦ (2)
Now, a piezoelectric plate sandwiched between two electrodes as shown in Fig. 1 is considered. Without loss of generality, several assumptions are considered for this configuration.

-the thickness to length ratio of the piezoelectric plate is less than 1/10, -the piezoelectric plate is confined to have planar movement and it is subjected to loading only in the xy plane, -the thickness of the electrodes is negligible in comparison with thickness of piezoelectric plate, -the electromechanical system is assumed to be linear, -the electrodes are perfectly conductive, -the polarization direction is perpendicular to the plate in parallel to → z axis, -the electrical field is uniform in the direction of thickness aligned with the polling direction, -the variation of the potential in the direction of the thickness is linear,

The first two assumptions let us use the plane stress assumption modeling for the piezoelectric plate [START_REF] Hutton | Fundamentals of finite element analysis[END_REF]. In this case, any stress in the direction of z will be zero. By considering transversely isotropic piezoelectric material and considering plane stress assumption, the piezoelectric plate has in-plane isotropic behavior. Furthermore, by poling the piezoelectric material in the z direction, the only non-zero electric field will be in the z direction. In this case, the reduced (2D) form of piezoelectric constitutive equation can be derived in the following form [START_REF] Junior | An electromechanical finite element model for piezoelectric energy harvester plates[END_REF]) The components of the reduced constitutive equation can be written as [START_REF] Junior | An electromechanical finite element model for piezoelectric energy harvester plates[END_REF])

⎡ ⎢ ⎢ ⎣ T 1 T 2 T 3 D 3 ⎤ ⎥ ⎥ ⎦ = ⎡ ⎢ ⎢ ⎣ c * 11 c * 12 0 -e *
⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ S 1 S 2 S 3 E 3 ⎤ ⎥ ⎥ ⎦ (3)
c * E 11 = c E 11 - c E 13 2 c E 33 , c * 12 = c E 12 - c E 13 2 c E 33 , c * 33 = c 66 e * 31 = e 31 -e 33 c E 13 c E 33 , ε * 33 = ε S 33 + e 2 33 c E 33 (4)
The obtained constitutive equation in (3) will be used in the finite element modeling of the piezoelectric plate which will be discussed in the next section.

Piezoelectric finite element model

In this section, to derive the finite element (FE) formulation, the piezoelectric plate is discretized by rectangular elements which are particular form of the more general 2D elements called "bilinear quadrilateral element" [START_REF] Hutton | Fundamentals of finite element analysis[END_REF][START_REF] Kattan | MATLAB guide to finite elements: an interactive approach[END_REF]). It should be noted that since the thickness of electrodes is negligible in comparison with the thickness of the piezoelectric plate, its structural effects are neglected in the modeling. Hence, only the piezoelectric plate is discretized by the finite number of elements. The schematic form of this discretization is illustrated in Fig. 2a. As can be seen in this figure, The piezoelectric plate is discretized as 3 by 4 elements. It is clear that a finer discretization will be used for the numerical optimization. It can be seen in the figure that each rectangular element has 4 nodes and each node has 2 in-plane mechanical degrees of freedom regarding the displacement in x and y directions. The rectangular element shown in Fig. 2b can have arbitrary length l e and width w e . In fact, with the method of finite element modeling which is used to write the optimization code, the dimensions of plate and number of elements can be defined separately. This freedom will have two advantages: first, for a predefined geometry of a piezoelectric plate, higher number of elements can be defined to have better results in terms of having small detail. Second, it is possible to define lower number of elements when reducing the computation time is necessary. As it is illustrated in Fig. 2c, the rectangular element is mapped to a parent element which is a square element with natural coordinates ξ and η. The displacement of every point within the element will be expressed by the displacement of the nodes through the interpolation functions in the following format [START_REF] Hutton | Fundamentals of finite element analysis[END_REF].

x = x 1 n 1 + x 2 n 2 + x 3 n 3 + x 4 n 4 y = y 1 n 1 + y 2 n 2 + y 3 n 3 + y 4 n 4 (5)
and the interpolation functions can be written based on the natural coordinates

n 1 = 1 4 (1 -ξ)(1 -η), n 2 = 1 4 (1 + ξ)(1 -η) n 3 = 1 4 (1 + ξ)(1 + η), n 4 = 1 4 (1 -ξ)(1 + η) (6)
where the matrix of interpolation function can be written as,

N = n 1 0 n 2 0 n 3 0 n 4 0 0 n 1 0 n 2 0 n 3 0 n 4 (7)
So far, for each element, just mechanical degrees of freedom are considered. However, in piezoelectric material, the mechanical and electrical fields are coupled. Therefore, the electrical degree of freedom should be modeled as well. The general approach for this case is to consider one electrical degree of freedom for each node in addition to mechanical degrees of freedom as explained in [START_REF] Lerch | Simulation of piezoelectric devices by two-and threedimensional finite elements[END_REF]. However, here by assuming that conductive electrodes are placed on top and bottom of the piezoelectric plates as shown in Fig. 1, the electrical potential over each electrode is constant. This condition is known as equipotential condition. Furthermore, by considering the bottom electrode as ground electrode, the whole piezoelectric plate will have one electrical degree of freedom. On the other hand, for the purpose of elemental sensitivity analysis which will be explained in Section 3, for each element, one electrical degree of freedom is considered as it is shown in yellow in Fig. 2. The global equipotential condition will be imposed after assembling the global matrices. Now, the strain and electrical field of each element can be expressed with the help of mechanical and electrical degrees of freedom

S = B u u, Ē = B φ φ (8)
In ( 8), u and φ are the vectors of mechanical displacement and scalar value of electric potential respectively. B u is the strain displacement matrix which is written as follows [START_REF] Kattan | MATLAB guide to finite elements: an interactive approach[END_REF],

B u = 1 |J | B 1 B 2 B 3 B 4 B i = ⎡ ⎢ ⎣ a ∂n i ∂ξ -b ∂n i ∂η 0 0 c ∂n i ∂η -d ∂n i ∂ξ c ∂n i ∂η -d ∂n i ∂η a ∂n i ∂ξ -b ∂n i ∂η ⎤ ⎥ ⎦ (9)
and the parameters a, b, c and d are given by [START_REF] Kattan | MATLAB guide to finite elements: an interactive approach[END_REF],

a = 1 4 [y 1 (ξ -1) + y 2 (-1 -ξ) + y 3 (1 + ξ) + y 4 (1 -ξ)] b = 1 4 [y 1 (η -1) + y 2 (1 -η) + y 3 (1 + η) + y 4 (-1 -η)] c = 1 4 [x 1 (η -1) + x 2 (1 -η) + x 3 (1 + η) + x 4 (-1 -η)] d = 1 4 [x 1 (ξ -1) + x 2 (-1 -ξ) + x 3 (1 + ξ) + x 4 (1 -ξ)] ( 10
)
where x i and y i are the coordinates of the nodes in the rectangular element before mapping. The determinant of Jacobian matrix J , which transfers the natural coordinates to the generalized coordinates, is

|J | = 1 8 x 1 x 2 x 3 x 4 × ⎡ ⎢ ⎢ ⎣ 0 1-η η -ξ ξ -1 η -1 0 ξ + 1 -ξ -η ξ -η -ξ -1 0 η + 1 1 -ξ ξ + η -η -1 0 ⎤ ⎥ ⎥ ⎦ ⎡ ⎢ ⎢ ⎣ y 1 y 2 y 3 y 4 ⎤ ⎥ ⎥ ⎦ (11)
By considering the last two assumptions, B φ is [START_REF] Junior | An electromechanical finite element model for piezoelectric energy harvester plates[END_REF])

B φ = 1/h ( 12
)
where h is the thickness of the piezoelectric plate. Now by using the Hamilton's variational principle and neglecting the damping effect, the linear differential equation for one single element can be written in the following form [START_REF] Lerch | Simulation of piezoelectric devices by two-and threedimensional finite elements[END_REF])

m 0 0 0 ü φ + k uu k uφ k φu -k φφ u φ = f q (13)
in which m is the mass matrix, k uu is the mechanical stiffness matrix, k uφ is the piezoelectric coupling matrix, k φφ is the dielectric stiffness matrix, f is the external mechanical force and q is the charge. These components of linear differential equation are derived in the following form

k uu = h A B T u c E B u |J | dξ dη k uφ = h A B T u e T B φ |J | dξ dη k φφ = h A B T φ ε S B φ |J | dξ dη m = ρh A N T N |J | dξ dη ( 14
)
where A is the top surface area of the element and ρ is the density of the material. In fact, ( 14) illustrates the analytical calculations of elemental matrices. However, for numerical implementation in MATLAB, two-point Gauss quadrature method [START_REF] Hutton | Fundamentals of finite element analysis[END_REF]) is utilized for calculation of the elemental matrices numerically, which gives the exact values. The implementation procedure is explained later in Section 4.

To have the global FEM equation for a whole piezoelectric plate, the elemental matrices in (13) should be assembled, which is a general procedure in the FEM methodology and which will also be explained in Section 4. After assembling the elemental matrices, the global finite element equation for the whole design domain can be written as

M 0 0 0 Ü Φ + K uu K uφ K φu -K φφ U Φ = F Q (15)
Now, for two cases of actuation and energy harvesting, the global FEM (15) can be interpreted in different ways.

Here, we focus on static actuation so that the dynamics will not be considered. Therefore, the global FEM equation for the actuation can be written as

K uu U + K uφ Φ = F (16)
This equation will be used to calculate the mechanical displacement due to applied potential.

For the energy harvesting case, the external charge (Q) is considered to be zero. In addition, the external force is considered to be a harmonic excitation of frequency Ω. In this case, by considering a linear electromechanical system, the force and the response of the system can be stated as

F = f 0 e iΩt U = u 0 e iΩt , Φ = φ 0 e iΩt ( 17
)
where f 0 , u 0 and φ 0 are the amplitude of harmonic force, displacement and potential. By substituting the (17) in (15), the global FEM equation for the energy harvesting case can be written as

-Ω 2 M 0 0 0 U Φ + K uu K uφ K φu -K φφ U Φ = F 0 (18)
Equation ( 18) can also be written in following form

K uu -MΩ 2 K uφ K φu -K φφ U Φ = F 0 (19)
To solve the FEM (19), the mechanical boundary condition and equipotential condition should be applied. This will be explained in detail in Section 4.

Normalization

Here, the critical point is that the scale difference between the piezoelectric matrices including the mechanical stiffness matrices (k uu ) and (k uφ ) and the dielectric stiffness matrix (k φφ ) is huge. This huge scale difference can bring numerical instabilities in form of singularities in solving the final global FEM equation during the optimization loops. To eliminate the scale difference, a normalization is suggested (Homayouni-Amlashi 2019; Homayouni-Amlashi et al. 2020b) by factorizing the highest value of each elemental matrix which can be expressed in the following format

kuu = k uu /k 0 , kuφ = k uφ /α 0 kφφ = k φφ /β 0 , m = m/m 0 (20)
Starting by this normalization of elemental matrices, the actuation FEM ( 16) can be rewritten as

Kuu Ũ + Kuφ Φ = F (21)
in which

F = F /f 0 , Ũ = U/u 0 , Φ = Φ/φ 0 u 0 = f 0 /k 0 , φ 0 = f 0 /α 0 ( 22
)
The same normalization can be performed on the energy harvesting FEM ( 19)

Kuu -M Ω2 Kuφ Kφu -γ Kφφ Ũ Φ = F 0 ( 23
)
where

Ω2 = Ω 2 m 0 /k 0 , γ = k 0 β 0 /α 2 0 ( 24
)
Here γ is a normalization factor which keeps the solution of the system equal before and after applying the normalization. This normalization factor is having the scale of 10 1 and in this way the scale difference between the piezoelectric matrices is eliminated. The proof of normalization is provided in the Appendix. Now, by having the FEM equations ( 21) and ( 23), it is possible to enter the optimization phase. In the upcoming sections, optimization of actuator and energy harvester is separated.

As explained before, topology optimization is about distribution of the material within a design domain while there is no prior knowledge of the final optimized layout of the structures [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF]. There are several approaches for topology optimization method [START_REF] Maute | Topology optimization approaches: a comparative review[END_REF]. However, the density-based approach is chosen for this paper since its efficiency is already established in many researches specially in the area of piezoelectric actuators [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF][START_REF] Ruiz | Optimal design of a microgripper-type actuator based on aln/si heterogeneous bimorph[END_REF][START_REF] Moretti | Topology optimization of piezoelectric bimaterial actuators with velocity feedback control[END_REF] or energy harvesters (Homayouni-Amlashi 2019; [START_REF] Zheng | Topology optimization of energy harvesting devices using piezoelectric materials[END_REF][START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF].

Piezoelectric material interpolation scheme

One of the famous material interpolation schemes is the density-based approach which has been introduced to relax the optimization from the binary (void-material) problem. In this approach, for a discretized design domain by finite number of elements, the material properties of each element are related to element's density through a power law interpolation function. For passive isotropic material, this interpolation function is referred to as solid isotropic material with penalization (SIMP) which relates the element's Young modulus of elasticity to its density. For active non-isotropic piezoelectric material, the interpolation function is the extension of SIMP scheme which can be written as follows [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF])

kuu (x) = E min + x p uu (E 0 -E min ) kuu kuφ (x, P ) = (e min + x p uφ (e 0 -e min ))(2P -1) p P kuφ kφφ (x) = (ε min + x p φφ (ε 0 -ε min )) kφφ m(x) = x m (25)
where E min , e min and ε min are small numbers to define the minimum values for stiffness, coupling and dielectric matrices while E 0 , e 0 and ε 0 are equal to 1 to define the maximum values of the respected matrices. The definition of minimum values is provided to avoid the singularities during the optimization iterations. x is the density ratio of each element which has a value between 0 and 1. P is the polarization variable which also has the value between 0 and 1 and determines the direction of polarization. p uu , p uφ , p φφ and p P are penalization coefficients for the stiffness, coupling, dielectric matrices and polarization value respectively. It is obvious that in (25), the normalized form of piezoelectric matrices is used. However, the interpolation function is true for non-normalized matrices as well.

The introduced material interpolation scheme is known as PEMAP-P (piezoelectric material with penalization and polarization) [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF], which is the extension of the SIMP approach. Although some projections are defined for SIMP to have a robust topology optimization [START_REF] Wang | On projection methods, convergence and robust formulations in topology optimization[END_REF], there are alternative interpolation functions like RAMP (rational approximation of material properties) [START_REF] Stolpe | An alternative interpolation scheme for minimum compliance topology optimization[END_REF], or newer interpolation function introduced in [START_REF] Clausen | Topology optimization of coated structures and material interface problems[END_REF]. This latter one is also used in topology optimization of piezoelectric transducers [START_REF] Donoso | Topology optimization of piezo modal transducers with null-polarity phases[END_REF]. In fact, the study on preference of these interpolation functions is not in the scope of this paper. But, with the proposed MATLAB code, it will be easy for implementation of other interpolation functions.

After establishing the material interpolation scheme, the rest of this section will be divided into two parts: actuation and energy harvesting.

Actuation

Problem formulation

Following the classic approach for compliant mechanisms reported in [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF], the optimization of a planar piezoelectric actuator can be defined simply as displacement optimization or minimization of the following objective function,

minimize J act = -L T Ũ Subject to V (x) = NE i=1 x i v i ≤ V 0 < x i ≤ 1 0 ≤ P i ≤ 1 ( 26
)
where L is a vector with a value of 1 that corresponds to the output displacement node and 0 otherwise. In addition, a constraint is defined on the final volume of the optimized design. V is the target volume which is a fraction of the overall volume of the design domain while v i is the volume of each element and NE is the total number of elements while i is the number of each element.

Sensitivity analysis

For applying a gradient-based optimization, the sensitivity of objective function with respect to the optimization variables should be calculated. As such, the sensitivity of objective function with respect to the x i can be derived as

∂J ∂x i = ∂ ∂x i -L T Ũ + Λ T Kuu Ũ + Kuφ Φ -F = ∂ ∂x i -L T + Λ T Kuu Ũ + Λ T Kuφ Φ -Λ T F (27)
Through using the procedure known as adjoint method, Λ is introduced to avoid taking the derivative of displacement with respect to design variable, i.e., ∂ ũi ∂x . Therefore, the following adjoint equation should be solved

-L T + Λ T Kuu = 0 (28)
where Λ is the global adjoint vector. By solving the adjoint (28), the sensitivity values can be obtained as

∂J ∂x i = λ T i ∂ kuu ∂x i ũi + λ T i ∂ kuφ ∂x i φi (29)
where λ i is the elemental format of global adjoint vector Λ. By using the same procedure, sensitivity analysis with respect to polarization P is derived in the following form as well

∂J ∂P i = λ T i ∂ kuφ ∂P i φi ( 30
)
where λ i is the same adjoint vector which is already calculated in (28).

Based on ( 29) and ( 30), the derivative of piezoelectric stiffness and coupling matrices with respect to design variables is required which can be derived with the help of (25) as

∂ kuu ∂x i = p uu (E 0 -E min )x p uu -1 i kuu ∂ kuφ ∂x i = p uφ (e 0 -e min )x p uφ -1 i (2P i -1) p P kuφ (31) ∂ kuφ ∂P i = 2p P (e 0 -e min )(2P i -1) p P -1 x p uφ i kuφ (32)
After performing the sensitivity analysis, and defining the constraint, the optimization of variables should be done by the optimization algorithms which will be discussed later in this section.

Energy harvesting

Problem formulation

In the case of energy harvesting optimization, electrical boundary condition should be applied in addition to mechanical boundary condition. As mentioned previously, by considering perfectly conductive electrodes, the equipotential boundary condition can be applied by using the Boolean matrix in the following form [START_REF] Cook | Topology optimization of bimorph piezoelectric energy harvesters considering variable electrode location[END_REF])

Φ = BV p (33)
in which V p is the voltage of the top electrode while the bottom electrode considered as ground. For a general case of multi-layer piezoelectric plates, Boolean matrix B is having the dimension of N e × N P where N e is the number of nodes and N P is the number of electrodes. However, in the case of the one piezoelectric plate of this paper, B will be a vector of ones.

It is worth to note that the equipotential condition given in (33) will be applied only to the energy harvesting FEM (19). For the case of actuation, the equipotential condition will be applied automatically by defining equal applied voltage for all of the elements.

By applying the equipotential condition in (33), the energy harvesting FEM (19) can be rewritten as

K uu K uφ K φu -K φφ Ũ V p = F 0 ( 34
)
where

K uu = Kuu -M Ω2 bc K uφ = Kuφ B bc K φφ = γ B T Kφφ B ( 35
)
in which ([ ] bc ) shows the application of mechanical boundary condition. Now, the objective function for the energy harvesting application should be defined. Generally, the energy conversion ratio [START_REF] Zheng | Topology optimization of energy harvesting devices using piezoelectric materials[END_REF][START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF] or electromechanical coupling coefficient (de Almeida 2019) which is equivalent mathematically is chosen as the objective function. However, this format of objective function suffers from numerical instabilities during optimization iterations where it is suggested to penalize the mechanical energy as suggested in de Almeida (2019) and [START_REF] Salas | A topology optimization formulation for transient design of multi-entry laminated piezocomposite energy harvesting devices coupled with electrical circuit[END_REF]. Therefore, to avoid the numerical instabilities here, a classical format of objective function is defined. The optimization is defined as minimization of the weighted sum of the mechanical and electrical energy of the system, minimize

J EH = w j Π S -(1 -w j )Π E Subject to V (x) = NE i=1 x i v i ≤ V 0 < x i ≤ 1 0 ≤ P i ≤ 1 (36)
Π E and Π S are electrical and mechanical energies respectively which are defined in the following form [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF][START_REF] Zheng | Topology optimization of energy harvesting devices using piezoelectric materials[END_REF])

Π S = 1 2 Ũ T K uu Ũ, Π E = 1 2 V T p K φφ V p (37)
In optimization of (36), w j is the weighing factor which has the value between 0 and 1. Choosing the value of 1 for w j will make the optimization problem, a minimum compliance problem in which the goal is to minimize the mechanical deflection of the system under the applied force. By decreasing the value of w j , more weight will be given to maximize the output electrical energy. However, choosing very small values for w j will result in mechanically unstable layouts. Therefore, the value for w j will be found by using trial and error approach. The basis for choosing this value can be the maximum energy conversion factor of the plate under the same force. For example, in Homayouni-Amlashi et al. (2020b), the maximum energy conversion factor for an optimized piezoelectric plate under planar excitation is found to be 0.03 while in [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF] for a twolayer optimized piezoelectric plate under the bending force this ratio is 0.1. Therefore, the initial value of w j for the trial error approach can be considered between 0.01 and 0.1. The final chosen value of w j depends on the maximum stress and strain induced by the defined mechanical input to the structure which can be revealed by the post processing analysis. On the other hand, stress and strain constraints can also be considered in the optimization problem as it is investigated by [START_REF] Wein | Topology optimization of a cantilevered piezoelectric energy harvester using stress norm constraints[END_REF].

In fact, the advantage of the objective function defined here is that the optimization algorithm converges very smoothly to the final result. However, the drawback of this objective function is that the obtained result can be suboptimal depending on the chosen value for the w j . On the other hand, the problem of sub-optimal results exists in other formats of the objective function. For example, [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF] showed that by considering the energy conversion factor (Π E /Π S ) as objective function, different values of penalization factors can produce different results.

Eventually, we believe that the chosen objective function suits the educational purpose of this paper. Indeed, with the help of provided MATLAB code, the readers can easily change the code to implement other objective functions.

Sensitivity analysis

Similar to the actuation case, the next step after defining the objective function is sensitivity analysis. Since the objective function in (36) consists of mechanical and electrical energies, the sensitivity of each energy with respect to density ratio x can be found as [START_REF] Zheng | Topology optimization of energy harvesting devices using piezoelectric materials[END_REF]; Homayouni-Amlashi 2019; Homayouni-Amlashi et al. 2020b)

∂Π S ∂x i = 1 2 ũT i + λ T 1,i ∂( kuu -m Ω2 ) ∂x i ũi +λ T 1,i ∂ kuφ ∂x i φi + μ T 1,i ∂ kφu ∂x i ũi -μ T 1,i γ ∂ kφφ ∂x i φi (38) ∂Π E ∂x i = 1 2 φT i γ ∂ kφφ ∂x i φi -μ T 2,i γ ∂ kφφ ∂x i φi +λ T 2,i ∂( kuu -m Ω2 ) ∂x i u i + λ T 2,i ∂ kuφ ∂x i φi + μ T 2,i ∂ kφu ∂x i ũi (39)
in which μ and λ are the elemental adjoint vectors which are calculated by the following global coupled system

K uu K uφ K φu -K φφ Λ 1 Υ 1 = -K uu Ũ 0 K uu K uφ K φu -K φφ Λ 2 Υ 2 = 0 -K φφ V p (40)
where Λ and Υ are the global adjoint vectors which need to be disassembled to form the elemental adjoint vectors

[λ 1 ] bc = Λ 1 , [λ 2 ] bc = Λ 2 , [μ 1 ] = BΥ 1 , [μ 2 ] = BΥ 2 (41)
Now, the sensitivities with respect to polarization (P ) are calculated as well (Homayouni-Amlashi et al. 2020b; Homayouni-Amlashi 2019)

∂Π S ∂P i = λ T 1,i ∂ kuφ ∂P i φi + μ T 1,i ∂ kφu ∂P i ũi ∂Π E ∂P i = λ T 2,i ∂ kuφ ∂P i φi + μ T 2,i ∂ kφu ∂P i ũi (42)
Based on sensitivity equations in ( 39) and ( 42), the derivative of all piezoelectric matrices with respect to the design variables is required. The derivative of stiffness and coupling matrices is found in ( 31) and ( 32). Here, the derivative of dielectric matrix and mass matrix is also required which are

∂ kφφ ∂x i = p φφ (ε 0 -ε min )x p φφ -1 i kφφ ∂ m ∂x i = mi (43)
In addition to derivative of piezoelectric matrices with respect to density, derivation of the piezoelectric coupling matrix with respect to polarization variable is also required

∂ kuφ ∂P i = 2p P (2P i -1) p P -1 x p uφ i kuφ (44)
After calculation of sensitivities, the optimization variables can be updated in each iteration of optimization with the help of optimization algorithm which is the subject of the next section.

Optimization algorithms

For solving the optimization problem, there are several optimization algorithm like sequential linear programming (SLP), sequential quadratic programming (SQP), method of moving asymptotes (MMA) or the optimality criteria (OC) method. This latter one is more historical than the other methods and its application is more simple. However, MMA is more powerful in terms of solving multi-variable and multi-constraint optimization problems. In addition, the convergence in MMA is more assured due to consideration of two past successive iterations during optimization.

Therefore, in this paper to solve the actuation problem, the OC method is implemented so the proposed code is self-working. However, the OC method has convergence problem in energy harvesting code due to coupling effect of piezoelectric material. In the upcoming two sections, the OC method and MMA and their implementation codes will be explained.

Optimality criteria method

Optimality criteria is a heuristic method to update the design variables in each element of design domain during each iteration of optimization. Here, there are two design variables for each element including the density and polarization. It is common for structural optimization with optimality criteria that the mutual influence of the design variables on each other and from element to element can be ignored [START_REF] Hassani | A review of homogenization and topology optimization iii-topology optimization using optimality criteria[END_REF]. Therefore, they can be updated separately from each other in each iteration of optimization. In fact, this is similar to the case of topology optimization by homogenization method where the length, width and rotation angle of the hole in a microstructure are optimized by OC separately in each iteration of optimization [START_REF] Suzuki | The method of moving asymptotes-a new method for structural optimization[END_REF].

Here, the densities should be chosen such a way to respect the volume constraints. In this case, the following Karush-Kuhn-Tucker (KKT) condition for intermediate densities (0 < x i < 1) should be satisfied to guarantee the convergence

∂J ∂x i + λ ∂V ∂x i = 0 (45)
in which λ is the Lagrange multiplier to augment the volume constraint in the optimization. To solve the optimization problem mentioned in KKT (45), the OC algorithm given by Bendsøe and Sigmund (1995) and [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF] is used here which can be written as

x new i = ⎧ ⎨ ⎩ max(0, x i -move) min(1, x i + move) x i β η i if x i β η i ≤ max(0, x i -move) if x i β η i ≥ min(1, x i -move) otherwise (46)
where move parameter is the maximum amount of density change in each iteration of optimization, η is a numerical damping coefficient and

β i = - ∂J ∂x i λ ∂V ∂x i -1 (47)
Similar to the classical compliant problems of passive material [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF], the values of move and η are considered to be 0.2 and 0.3 respectively.

For optimization of polarization there is no volume constraint. In fact, since the polarization of each element will be optimized separately, it is not necessary to force the optimization algorithm to push the polarization value to zero for the elements with low density values. Therefore, the optimization problem can be simply defined as follows

∂J ∂P i = 0 (48)
The algorithm for this simple optimization can be obtained by modifying the OC algorithm given in (46) in the following form

P new i = max(0, P i -move) min(1, P i + move) if ∂J ∂P i ≥ 0 if ∂J ∂P i < 0 (49)
By defining the optimization algorithm for polarization based on (49), the polarity value (P i ) for all of the elements will be steered to -1 or +1 even for the elements with very low density. However, this will not affect the optimization results since based on (25) and with maximum value of polarization (P ), low density (x) will push the coupling matrix to its minimum value.

The method of moving asymptotes

The method of moving asymptotes (MMA) is a structural optimization method proposed by Svanberg (1987). The problem formulation of this methodology is as follows [START_REF] Svanberg | Mma and gcmma-two methods for nonlinear optimization[END_REF])

minimize f 0 (χ) + a 0 z + m c i=1 c i y i + 1 2 d i y 2 i subject to f i (χ) -a i z -y i ≤ 0 i = 1, ..., m c χ ∈ X, y i ≥ 0, z ≥ 0 (50) in which X = χ ∈ n var |χ min i ≤ 0 ≤ χ max i
, while n var is the number of design variables, χ is the vector of all design variables, y i and z are the artificial optimization variables, f 0 (χ) is the cost function to be minimized, m c is the number of constraints and a 0 , a i , c i , d i are the coefficients which have to be determined to match the optimization problem mentioned in (50) to different types of optimization problems.

For the optimization problem of this paper, there is only one constraint defined on the maximum volume of the optimized design. Therefore, based on the description given in [START_REF] Svanberg | Mma and gcmma-two methods for nonlinear optimization[END_REF], by considering a 0 = 1 and a i = 0 for all i, then z = 0 in any optimal solution, and by considering d i = 0 and c i = "a large number," then the variables y i = 0 for all i and the optimization problem mentioned in (50) will be matched to the optimization problem of this paper.

Finally, the OC and MMA optimization algorithms can update the optimization variables in each iteration of optimization. To do so, the implementation code for both of OC and MMA will be given in Section 4.

Filtering

Like other problems in structural optimization, topology optimization of piezoelectric materials also suffers from numerical instabilities such as checker board problem or mesh dependency. To remedy, the solutions need to be filtered in each iteration of optimization. So far, many filtering methods are suggested in the literature [START_REF] Sigmund | Morphology-based black and white filters for topology optimization[END_REF]). Among the proposed filtering methods, sensitivity filter and density filter proved their success in overcoming the aforementioned numerical instabilities. The sensitivity filter is used in the 99 lines of topology optimization code written by [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF] and density filter is used as an alternative option in the 88 lines of topology optimization code written by [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]. As such, in this paper, the density and sensitivity filters are used as two available options. To implement these filters, the same lines of codes written in 88 lines of code by [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF] are employed. Therefore, the concepts behind these filters will not be explained here since the related explanations can be found in the mentioned references.

The important point here is that the density filter shows more promising performance in comparison with sensitivity filter specially in the case of energy harvesting problem with using the MMA as the optimization algorithm.

Finally, the filters will be applied only to the densities and filtering of the polarization variable is not necessary.

Choosing proper penalization factors

After preparing the ingredients of topology optimization algorithm, the important question would be how to tune the parameters of optimization? in addition to parameters which are the same for optimization of passive and active materials like radius of filtering and volume fraction, there are other parameters which belong to multiphysics nature of piezoelectric topology optimization like penalization coefficients of the piezoelectric matrices, i.e., p uφ , p φφ and p P . In fact by choosing the penalization coefficients, we are pursuing two goals: (1) guaranteeing the convergence to perfect void-material in the final obtained layout, (2) avoiding local optima. For choosing p uu , it is already proven that the value of 3 is the best choice [START_REF] Bendsoe | Generating optimal topologies in structural design using a homogenization method[END_REF] to reach a perfect void/material in the final layout.

There are several studies which focused on defining the criteria for choosing other penalization coefficients, i.e., p uφ , p φφ and p P [START_REF] Kim | Topology optimization of piezoelectric materials and application to the cantilever beams for vibration energy harvesting[END_REF][START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF][START_REF] Kim | Multi-physics interpolation for the topology optimization of piezoelectric systems[END_REF]. In particular, [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF] chose the penalization factors randomly and obtained different density layouts. The final conclusion was that the penalization factors have extreme effect on the final topology. But, no criteria or rule is presented in this study about the method of choosing the penalization factors. On the other hand, a detailed study is presented by [START_REF] Kim | Multi-physics interpolation for the topology optimization of piezoelectric systems[END_REF] for choosing the penalization factors. Based on this study, necessary condition for choosing the penalization factors is that the electromechanical coupling coefficient (EMCC) should be increased when the density (x) increases and vice versa. Based on this condition, the following intrinsic condition which is independent from objective function is proposed for choosing the penalization factors [START_REF] Kim | Multi-physics interpolation for the topology optimization of piezoelectric systems[END_REF] 2p uφp uu + p φφ > 0 (51)

In addition to this condition, [START_REF] Kim | Multi-physics interpolation for the topology optimization of piezoelectric systems[END_REF] proposed other objective dependent criteria for the actuation and energy harvesting application. In particular, for the actuation objective function in (26), the following condition is proposed [START_REF] Kim | Multi-physics interpolation for the topology optimization of piezoelectric systems[END_REF],

p uφ -p uu > 0 (52)
and for energy harvesting approach, when the objective function is energy conversion factor then the following conditions should be satisfied [START_REF] Kim | Multi-physics interpolation for the topology optimization of piezoelectric systems[END_REF],

p uφ -p uu > 0 p φφ -p uu > 0 (53)
It should be noted that these conditions are proposed to guarantee the convergence of the final topology and yet there is no study on the methods to define these penalization factors to avoid the local optima.

For polarization penalization, the value of 1 seems the best choice as it is suggested by [START_REF] Kögl | Topology optimization of smart structures: design of piezoelectric plate and shell actuators[END_REF]. Indeed, it is not necessary to penalize the polarization.

MATLAB implementation codes

In this section, the goal is to establish the MATLAB code for the piezoelectric optimization methodology explained in the previous section. Two MATLAB codes are mentioned in the Appendix of the paper. The first code is for actuation and the second code is for energy harvesting. Each MATLAB code is partitioned so the readers can have a perception of the goal of each part. The important lines of codes are also labeled to provide a connection between the line and the analytical calculation. Again, the rest of this section will be divided for actuation and energy harvesting codes and in each section different parts of the codes will be explained more in detail.

Actuation

General definition

The first part of the code is GENERAL DEFINITION. This part consists of geometrical dimensions of the piezoelectric plate, the resolution of the mesh, i.e., the number of elements in the → x and → y axes, the penalization factors, the type of filter, the filter radius and maximum number of iterations and the stiffness of the attached spring are defined. It should be noted that two stopping criteria are defined by the codes. The classical criteria is the density change between two last successive iterations. The other criteria is the maximum number of iterations. Satisfying either of these will stop the optimization. The reason behind adding the second criteria is that in contrast to pure mechanical problems, here we have different types of material and different types of objective functions. As such, the density change will not stop the optimization generally or it needs very high number of iterations while there are no significant oscillations in the objective function's value.

Material properties

The second part of the code is MATERIAL PROPERTIES. In this part of the code, the properties of the chosen piezoelectric material are given. These properties are the density, e * 31 coupling coefficient, * 33 permittivity coefficient and the elements of the stiffness tensor all after applying the plane stress assumption as mentioned in (4). The chosen material for the code is the PZT 4 which is popular in the literature. There are PZT materials with lower and higher coupling coefficients. If one wants to investigate the results of other piezoelectric materials, the values of this part can be changed.

It should be noted that choosing different PZT materials in this section will not affect the final layout obtained by optimization algorithm. In fact the resulted layout of the optimization is independent from Piezoelectric coefficients.

Finite element model

In the section called PREPARE FINITE ELEMENT ANALYSIS, the proposed code for the finite element model of the piezoelectric plate is the extension of the MATLAB code provided by [START_REF] Kattan | MATLAB guide to finite elements: an interactive approach[END_REF] for the bilinear quadrilateral elements of the passive materials. However, instead of using analytical calculations, elemental matrices are calculated by using the two-point Gauss quadrature method. To do so, first the geometrical sizes of each element are calculated based on the defined geometry of the design domain and desired resolution of mesh. Then, Gauss quadrature points are defined in matrix GP. Afterwards, the elemental matrices are found inside a loop with the help of Gauss points as follows,

In this method of elemental matrix calculation, there is no analytical integration, the calculation time is very fast and at the same time the code is flexible in terms of different geometrical dimensions and mesh resolution and aspect ratios.

After calculating the elemental matrices, normalization is applied based on (20) and normalization factors are obtained with the following lines of code It should be noted that for simplicity, the non-normalized matrices are replaced by normalized matrices. Since the normalization factors are saved, it is easy to find the real values after the calculation of the final results.

So far, the elemental matrices are calculated, normalized and the normalization factors are derived. For assembling the elemental matrices, the element connectivity matrix known as edofMat [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]) is built with the help of the numbering format shown in Fig. 2. In this figure, the coarse discretization is shown with the numbering of elements, nodes and mechanical and electrical degrees of freedom. This numbering format is similar to the one presented in 99 [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]) and 88 [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]) lines of MATLAB codes for structural problems. The numberings of the nodes and mechanical degrees of freedom are from top to bottom and from left to right. It should be noted that inside each element, the numbering order is different as it is shown in Fig. 2b. The numbering inside each element is counterclockwise and it starts from bottom left and the sequence of numbers in each row of edofMat matrix is following this numbering order. Furthermore, since the system here is having additional degree of freedom as potential, a vector of potential connectivity is defined as edofMatPZT. In fact, by considering one potential degree of freedom for each node, potential connectivity would be a matrix with rows containing the node IDs. However, due to equipotential condition, one potential degree of freedom is considered for each element and the potential connectivity is a vector. At last, the connectivity matrix edofMat and the potential connectivity vector edofMatPZT is defined as

The implementation lines to create the edofMat matrix and edofMatPZT vector with the help of nodeIDs are where nele is the number of elements, nodenrs is the matrix of node numbers and edofVec is the vector of nodeIDs containing the first nodeID of each element. Now, with the edofMat and edofMatPZT, it is possible to assemble the elemental matrices and build the global matrices with following lines of code

The matrices iK and jK corresponds to (i, j ) entry of stiffness matrix for each element. With the same strategy, iKup and jKup are written for the ith and jth entry of K uφ matrix. To avoid redundancy, there is no iKpp and jKpp. In fact, these latter are just equal to the edofMatPZT.

The material interpolation scheme mentioned in (25), is applied in sKuu and sKup matrices. xPhys(:) and pol(:) are the vectorized physical density and polarization matrices respectively which will be updated in each iteration of optimization.

Boundary condition

The mechanical boundary conditions are defined in the part DEFINITION OF BOUNDARY CONDITION. In this part, fixeddofs and freedofs are defined to contain the fixed degrees of freedom and the free degrees of freedom respectively. They are defined by using the method reported in 88 lines [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]) and 99 lines [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF]) of code. Therefore, for applying the clamped boundary condition on the left side of the design domain, the following lines of code can be used, By considering that for the coarse mesh in Fig. 2, nelx = 4 and nely = 3, fixeddofs in the aforementioned implementation line produce the numbers from 1 to 8 which are the mechanical degrees of freedom in the left side of the design domain.

Output displacement definition

In this part of the code for actuation, the goal is to define the particular point of the design domain where the maximization of displacement in a particular direction is desired. To do so, the desired mechanical degree of freedom should be defined. The variable DMDOF which is the abbreviation of desired mechanical degree of freedom is defined for this purpose. Thereafter, the vector L is created. This vector will be used later in the objective function and sensitivity analysis for optimization.

A spring is attached at DMDOF which simulates the reaction force exerted by an imaginary object. The modeled spring will modify the piezoelectric stiffness matrix with this line of code where Ks is the stiffness of the modeled spring. Actually, by changing the stiffness of the modeled spring, it is possible to determine whether more force is desired or more displacement. Indeed, since the stiffness matrix is normalized, the stiffness of the modeled spring can be determined with respect to the stiffness of the piezoelectric plate. For example, by putting the stiffness of the spring equal to 1, the stiffness will be equal to the highest value of the piezoelectric stiffness matrix. In this case, the piezoelectric layout will be optimized to produce maximum possible of force. However, by defining very low values of stiffness for the spring (i.e., 0.01), then the piezoelectric layout will be optimized for maximum possible of deflection.

Objective function

Calculation of objective function for actuator is a routine procedure to calculate the desired displacement in each iteration. To do so, first the mechanical displacement due to applied voltage should be calculated. Therefore, based on (21), the mechanical displacement vector is calculated by the following line of code then based on (36), the objective function can be calculated by

Sensitivity analysis

For sensitivity analysis, the first step is the calculation of adjoint vectors. It is calculated based on (28), then the calculation of sensitivities based on ( 29) and ( 30) starts afterwards where dc and dp are the sensitivities of the objective function with respect to x and P and dv is the sensitivities of the volume constraint with respect to x. It is obvious from ( 26) that the sensitivity of volume constraint with respect to P is zero.

As can be seen in the aforementioned lines of the code, the sensitivities containing the stiffness (kuu) and coupling matrix (kup) are calculated separately. This will help us to define different penalization coefficients for each of these matrices.

Optimization algorithm

The implementation lines of code for the OC update of densities are the same as what is mentioned in 88 lines of code [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]). In addition to densities, for updating the polarization based on algorithm (49), there is no need for the bi-sectioning loop and just one line of code can optimize the polarization in each iteration of optimization as follows

Plot densities and polarization

In part (PLOT DENSITIES & POLARIZATION), the densities and polarization will be plotted in two figures separately. As it is common, in the density figure, white area means no material while the black area means material.

Different color spectrum is chosen for polarization profile in which the red and blue color shows opposite direction of optimization while the green color shows the neutral material with no polarization. The important point here is that the OC will steer the polarization value to 0 or 1 even for the elements with lower density. Therefore, here to eliminate the confusion, the matrix of densities is multiplied to the polarization. In this way, for element with minimum density (no material), the polarization turns to green as well.

Filtering

The density filtering lines of code are similar to 88 lines of code [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]. Two types of density and sensitivity filter can be chosen. In general, the density filter is more recommended. However, since the code is written for both of these filters, the best choice will be up to readers.

Energy harvesting

General definition

In general definition part of the energy harvesting code, in addition to what is mentioned for the actuation code, the objective function's weighing factor (wj) and excitation frequency (omega) are defined.

Material properties

This part of the energy harvesting code is similar to the actuation code. However, the results of actuation are free from the PZT coefficients while in the case of energy harvesting for different materials and different piezoelectric coefficients, the results will be changed. PZT materials with higher coupling coefficients produce more electrical energy in comparison with ones with lower coupling coefficients. At the same time, PZT materials with higher coupling coefficients have more coupling effect which affect the optimization during the iteration and can change the final layout. Similar to actuation code, PZT 4 is chosen for the code.

Finite element model

In comparison with FE model presented in the actuation code, the FEM part of the energy harvesting is the same with some additional calculations. Indeed, the elemental matrices for mass and dielectric matrices are calculated as well in the following lines of code, For assembling the global mass matrix M, the same matrices of iK and jK as described in the actuation code can be used. The assembling lines of mass matrix and dielectric matrix Kφφ can be found in energy harvesting code as where the related material interpolation scheme is applied in sM and sKpp.

Furthermore, the stiffness matrix is modified by the dynamic matrix -omega*M which represents (-M Ω2 ) in ( 23). This modification is done by the following line of code After building the global piezoelectric matrices, the equipotential condition is applied with the following line of code where KupEqui and KppEqui are the piezoelectric coupling and dielectric matrices after applying the mechanical and equipotential condition which represent K uφ and K φφ .

Objective function

To find the value of objective function, the mechanical displacement and resulted potential due to applied force should be calculated based on (23). But, first a total matrix is built based on the left-hand side of ( 23). The symmetry is guaranteed afterward, and finally, the mechanical displacements and potential are calculated After calculation of mechanical displacement and potential, for calculation of objective function, Wm and We are defined for mechanical and electrical energies respectively.

It is important to note that in the written code, the mechanical and electrical energies are calculated by summing the energies in each element. That is why the mechanical displacement and electrical potentials are multiplied to elemental stiffness matrix kuu and elemental dielectric matrix kpp. Then, the objective function for the energy harvester is the weighted sum of the energies while considering the weighting factor (wj).

Sensitivity analysis

Sensitivity analysis starts by calculating adjoint vectors with the help of (40)

Here, B is the Boolean matrix as defined in (41). In fact, mu1 and mu2 are vectors of equal values for all of the elements in the design domain due to equipotential condition. For sensitivity analysis, the adjoint vectors related to mechanical and electrical states should be separated. For this reason, ADJ resolved to lambda and mu.

After calculation of adjoint vectors, the sensitivity analysis is performed by the following lines of code

The code is written in a general form to consider nf number of load cases. Therefore, for each number of load cases the mechanical displacement, potential and adjoint vectors are separated by defining the Uu i, Up i, lambda1 i, lambda2 i, mu1 i and mu2 i. In addition, mechanical displacement, potential and adjoint vectors are calculated in global format while the sensitivity analysis should be performed on the elemental scale. To do so, edofMat is utilized to convert the vector of mechanical displacement to the matrix Uu i(edofMat) in which each row is related to one element while the columns represent the mechanical displacement of that element. The same idea is also applied to the adjoint vectors, i.e., lambda1 i(edofMat) and lambda2 i(edofMat). On the other hand, due to equipotential condition, it is obvious that the potential of each element is equal to Up, and as mentioned before, values of mu1 and mu2 are equal for all elements.

Similar to the actuation code, here the sensitivities containing each piezoelectric matrix are calculated separately to facilitate the definition of different penalization factors.

Optimization algorithm

In the energy harvesting code, the MMA optimization algorithm is used. However, this code cannot be executed without having the external MMA code. The MMA implementation MATLAB code can be obtained by contacting Prof. Krister Svanberg. In this way, the interested readers will obtain two MATLAB codes to implement the MMA optimization code which includes mmasub.m and subsolv.m. The version 2007 of these codes is used in this paper. Supposing that these two MATLAB codes are already available, to implement the MMA, first initial parameters are defined in MMA Preparation part based on what is mentioned in Section 3.4.2. Then, optimization of variables will be done in MMA OPTIMIZATION OF DESIGN VARIABLES part. The important point in this part is that, due to the scale difference between the density sensitivities and polarization sensitivity, it is suggested here to normalize the polarization sensitivity with the following line of code Then, the main lines of MMA optimization implementation are

The output xmmaa is the vector of all optimization variables. Therefore, the upper half of this vector is the updated density and lower half is the polarization values. c, e, g, i andd, f, h, j respectively present the density and polarization profile of the design for the specified output stiffness after convergence. For the polarization profile, blue, red and green represent respectively negative, positive and null polarization After updating the optimization variables and plotting the results, this iteration of optimization will be finished and the optimization will be started from the beginning of the loop for the next iteration.

Numerical examples

In this section, the goal is to investigate the performances of the codes in different application cases of actuation and energy harvesting with different configurations. First, different examples of actuation will be investigated. Thereafter, different configurations of energy harvesting are explored.

Actuation

Pusher

The first example in the actuation part is a simple pusher as can be seen in Fig. 3a. The gray area shows the design domain which can be optimized by the optimization algorithm. The actuation optimization code which is mentioned in the Appendix is written for this example. As it can be seen from the code, the aspect ratio of the elements in the x and y directions is following the aspect ratio between the length and width of the plate which produces square elements for discretization of the design domain. This is not mandatory, but it is known that increasing the aspect ratio of width and length of each element increases the inaccuracy of the finite element model [START_REF] Logan | A first course in the finite element method using Algor[END_REF]. Therefore, it is recommended to follow the aspect ratio of the plate in defining the number of elements in x and y directions.

The chosen penalization factors are p uu = 3 and p uφ = 4, which satisfy the conditions mentioned in (52). These penalization factors are the same for all of the actuation examples.

For having a completely symmetrical response with respect to the horizontal dotted line in Fig. 3a and to decrease the number of elements in the design domain, the defined design domain in the code is the upper half of the piezoelectric plate. Therefore, in the symmetry line of the design domain, the roller mechanical boundary condition is applied with the following line of code, With this mechanical boundary condition, the nodes connected to the symmetry line can have displacement in the x direction but not in the y direction.

As can be seen in Fig. 3, the results of topology optimization for different spring stiffness are plotted. The upper row of the figure shows the density layout while the lower row shows the polarization profile. By mirroring the obtained result with respect to the symmetry line, the results are illustrated for whole piezoelectric plate.

The numerical results for different spring stiffness are also reported in Table 1. The numerical results show the amplification ratio of the optimized design with respect to the full plate under application of same value of voltage. The objective value which is reported by the code is not showing this amplification ratio. To calculate the amplification ratio, the final value of the objective function after finishing the optimization should be divided by the objective function value of the full plate. To find the objective function of the full plate, it is possible to define the initial values of density equal to 1. To do so, the following line of code put 1 instead of volfrac, then stop the code after the calculation of objective function. In this case, the value of objective function for the full plate is obtained.

From the plots in Fig. 3c andd, it is clear that when the spring stiffness is one, the optimal layout is a very simple lumped design with uniform polarization profile. Based on Table 1, for this case, the full plate is having more displacement. On the other hand, by considering very low stiffness (Ks = 0.005) and then based on Fig. 3i andj, the density layout is more complicated and polarization profile is not uniform anymore. In fact, it is obvious that the blue region in the polarization profile will have extension while the red part will have compression (shrinkage). The combination of this extension and compression will produce an amplification ratio with respect to full plate equal to 2.75 as reported in Table 1.

Gripper

The second example of piezoelectric actuation is a gripper, which is similar to the case discussed in [START_REF] Ruiz | Optimal design of a microgripper-type actuator based on aln/si heterogeneous bimorph[END_REF]. The goal is to design a gripper to grab an object as it is shown in Fig. 4a andb. To do so, some modifications should be done to the actuation code in the Appendix. First of all, the OUTPUT DISPLACEMENT DEFINITION part should be completely changed by replacing the following lines where L grip and W grip are the length and width of the empty box in the piezoelectric plate as shown in Fig. 4a andb.

Next, to enforce zero material in the desired box of the design domain, the passive elements should be defined. The strategy is the same as in 99 lines [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF] and in 88 lines [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]) of code. The following part should be added after the part INITIALIZE ITERATION and before the part START ITERATION Then to apply the passive material in each iteration, the following line should be added after the OC update line Now by executing the code, the results of Fig. 4 for different spring stiffness will be obtained and the numerical results are reported in Table 1. It is interesting to note that for the gripper, the amplification ratios in optimized designs are much higher than the amplification ratios for optimized pushers. Indeed, the polarization optimization plays a major role in designing the gripper. That is why for any chosen values of Ks, the polarization profile is not uniform and the Fig. 4 Topology optimization of a piezoelectric gripper for different stiffness of the modeled spring. The panels follow the same presentation as Fig. 3 gripper needs the combination of expansion and retraction for increasing the amplification ratio.

Energy harvesting

Lateral force

The first example of the energy harvesting code is a plate under a lateral force excitation as it is shown in Fig. 5. The code in the Appendix is written for this case. Here, the goal is to maximize the output electrical energy while minimizing the mechanical energy of the system. For this case, the problem is static and the excitation frequency is considered to be zero. The chosen penalization factors for the energy harvesting code in contrast to actuation part are not the same for all cases. As such in Table 2, the penalization factors are reported for each case. But, for all cases, the chosen penalization factors are satisfying the conditions mentioned in ( 51) and ( 52). The reason for different penalization factors for each case is that the energy harvesting optimization is more complicated than actuation due to existence of the coupling effect. This coupling effect is highly affected by the chosen penalization factors in particular p uφ and p φφ . By choosing proper penalization factors, it is possible to avoid the nonsymmetric results or to improve the convergence.

The results of optimization for different values of the weighting factor (w j ) are illustrated in Fig. 5. For the first case, the weighting (w j ) is equal to 1. As such, the problem is now a compliance problem in which minimization of deflection is the target. In this case, the optimization is done without polarization optimization. To do optimization without polarization, it is possible to simply put the penalization factor for the polarization equal to 0, i.e., penalPol = 0 in GENERAL DEFINITIONS part of the code. As can be seen in Fig. 5b, the obtained density layout for this case is similar to the results of the topology optimization of passive materials as reported by 99 lines [START_REF] Sigmund | A 99 line topology optimization code written in matlab[END_REF] or 88 lines [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF]) of MATLAB code. This was expected since the PZT materials have the plane isotropic behavior. The numerical results of the optimizations are given in Table 2. It is reported for the aforementioned case that the output electrical energy is zero which is due to the charge cancellation. In fact, lateral force induces tension and compression in different parts of the piezoelectric plate which produces voltages with opposite sign on the surface of the electrode. The opposite signs of voltages nullify each other. For the next case, polarization is also optimized by putting penalPol = 1. In Fig. 5c andd, it is obvious that the density layout did not change and the polarization profile is not uniform anymore. By this polarization optimization, based on Table 2, not only the mechanical energy of the piezoelectric plate is reduced in comparison with the first case, but also the problem of charge cancellation is suppressed and we have a non-zero electrical energy.

In the next case, the goal is to increase the electrical energy due to the same amount of force. To do so, the weighting factor is decreased to 0.01. As can be seen in Fig. 5e andf, the density layout is changed and the polarization profile is changed accordingly as well. By observing the obtained numerical results in Table 2, the electrical energy is increased in comparison with the previous cases with the cost of increasing the mechanical energy as well. By further increasing the weighting factor, the results of Fig. 5g and h will be obtained in which a jump in both mechanical energy and electrical energy of the design can be seen in Table 2.

It can be noticed that in the optimized polarization profiles, there are areas with null polarity (green color) while there are materials (non-void). This null polarity is mostly at places where there is a transition between the polarization direction. It is possible to take these null polarity areas into consideration in the optimization problem formulation as discussed in Donoso and Sigmund (2016).

In the next example, the piezoelectric plate with the same configuration of boundary and load condition is considered while the force is considered to be harmonic. In Fig. 6, the results of topology optimization under harmonic force can be seen. The related numerical results are also reported in Table 2. In this figure, the excitation frequency is primarily considered to be 1 kHz. At this frequency, the optimization converges to the final layout for different values of w j . However, by increasing the frequency of excitation, convergence problems begin. The problems are parasitic effects of the material layout and the force will be disjointed from the material. These problems are due to the fact that the excitation frequency is close to the resonance and antiresonance frequency of the piezoelectric plate. The resonances are the natural frequencies when the electrodes are short-circuited and antiresonances are the ones when the electrodes are in open-circuit condition [START_REF] Lerch | Simulation of piezoelectric devices by two-and threedimensional finite elements[END_REF]). During the optimization, these frequencies are changing in each iteration. Therefore, it is possible that the excitation frequency comes close to the resonance frequency during the optimization which will introduce singularity in the FEM ( 23) and the amplitude of the displacement vector reaches to infinity. Even by defining the damping, still the jump in the displacement vector will result in numerical instabilities and disjoint problem in which the force will be disconnected from material as it is reported by [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF]. One solution to reduce the numerical instabilities is to restrict the move limit of the MMA code. Indeed, the move limit of the MMA optimization code is by default set to 1 which means the density of material can jump from 0 (void) to 1 (material) in a single iteration. To modify this move limit, in the mmasub.m, one needs to change the value of move from 1 to 0.1. With this restriction on the move limit, convergence to a black and white final layout is achieved for 3.5 kHz excitation frequencies as it is shown in Fig. 6f andg. However, by reducing the MMA move limit, the final result can be trapped in the local optima. This is the maximum excitation frequency that the convergence can be achieved. After this frequency, again the force will be disjointed from the material and the convergence problem can be seen. To overcome the challenges of dynamic topology optimization of mechanical structures, several methods can be found in the literature. For example, [START_REF] Olhoff | Topological design of continuum structures subjected to forced vibration[END_REF] suggest that the excitation frequency can be increased gradually during the optimization iterations. [START_REF] Liu | A comparative study of dynamic analysis methods for structural topology optimization under harmonic force excitations[END_REF] and [START_REF] Jensen | Topology optimization of dynamics problems with padé approximants[END_REF] modeled the damping in the dynamic system and optimized the structure for an interval of frequencies including the resonance ones. For mechanical structures, the goal is to reduce the displacement or stored mechanical energy in the system. For piezoelectric energy harvesters in which the maximization of electrical output regarding the mechanical input is desired, [START_REF] Noh | Topology optimization of piezoelectric energy harvesting devices considering static and harmonic dynamic loads[END_REF] solved the disjoint problem of force and material in dynamic topology optimization by defining a constraint on the mechanical energy of the piezoelectric structure and they defined the objective function as maximization of the electrical output of the system.

Reducing the weighting factor w j brings the resonance frequency closer to the excitation frequency to some extent. However, a very low weighing factor can introduce convergence problems as well. On the other hand, for piezoelectric energy harvesters, matching the resonance frequency and excitation frequency is favorable. To do so, the problem formulation should be changed. For example, [START_REF] Kim | Topology optimization of piezoelectric materials and application to the cantilever beams for vibration energy harvesting[END_REF] did eigenfrequency optimization to increase the electromechanical coupling coefficient of the design and to match the resonance frequency and the excitation frequency. [START_REF] Wang | Topology optimization of piezocomposite resonator for maximizing excitation strength and synthesizing desired eigenmodes[END_REF] and [START_REF] Nakasone | Dynamic design of piezoelectric laminated sensors and actuators using topology optimization[END_REF]

Two loads case

In this example, the goal is to optimize piezoelectric energy harvesters for two loads case as shown in Fig. 7a. In fact, the piezoelectric plate will be optimized for in-plane forces that can come from different directions as it is discussed in Homayouni-Amlashi et al. (2020b) since any in-plane force can be decomposed to the load cases shown in Fig. 7a. The mechanical boundary condition of this figure is proposed to make the harvested energy symmetric with respect to the forces in each direction as it is possible.

To implement the two loads case, the following changes should be made to the energy harvesting code mentioned in the Appendix. The parameters in the GENERAL DEFINITIONS part should be changed by putting Lp = 3e-2, Wp = 3e-2, nelx = 100, and nely = 100. To define the boundary condition, the DEFINITION OF BOUNDARY CONDITION part will be changed completely as follows BCTratio is the ratio of the clamped part to the total length of the edge. To define the two load case the part FORCE DEFINITION will be changed with the following lines:

In contrast to previous case, the initial values for the densities (x {0} ) are not equal to volume fraction (volfrac). It is observed that by changing the initial values to one, the obtained results are more symmetric. Therefore, the following line of energy harvesting code is updated,

The results for the defined load and boundary condition are illustrated in Fig. 7. Figure 7b shows the compliance result due to static force without polarization optimization. In Fig. 7c andd the excitation frequency is zero, w j is decreased and polarization is also optimized. Other cases are optimized for different excitation frequencies and weighting factors. The numerical results are reported in Table 2. For defined geometry, load and boundary condition of this example, the highest frequency that the convergence to black and white is achieved is 1.5 kHz. This is due to the fact that with the defined geometry and boundary condition, the first resonance frequency is lower than the previous example of Fig. 6.

Discussion

In the actuation code, the optimization algorithm is OC. However, it is very simple to apply the MMA to the actuation code. One has just to copy the MMA Preparation part from the energy harvesting code and paste it after the INITIALIZE ITERATION part of the actuation code. Then, one has to remove the OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES part and substitute the MMA OPTIMIZATION OF DESIGN VARIABLES part from the energy harvesting code. It will be observed that the obtained result remains the same.

It is possible to implement the OC in the energy harvesting code as well. This can be simply done by substituting the MMA OPTIMIZATION OF DESIGN VARIABLES part with the OPTIMALITY CRITERIA UPDATE OF DESIGN VARIABLES part of the actuation code. However, due to coupling effect, sometimes the convergence problem appears.

It should be noted that the results presented in this paper are not the best results that can be obtained from the codes. In fact by changing the penalization factors, volume fraction, filter radius and initial values, different results can be obtained. However, finding the best results or performing the parameters analysis on the final results is not the subject of this educational paper.

There are similar aspects between the proposed codes here and the 88 lines of MATLAB code [START_REF] Andreassen | Efficient topology optimization in matlab using 88 lines of code[END_REF], such as the use of different filtering methods, different boundary conditions, and different load cases. These aspects are not discussed here since the procedures of implementation are similar.

Extension of the codes for different goals such as different objective functions, multi-material optimization, and 3D finite element modeling is straightforward. For example, the 169 lines of code for 3D topology optimization [START_REF] Liu | An efficient 3d topology optimization code written in matlab[END_REF] of passive material are the extension of the 88 lines of code. The same strategy can be used for the proposed code here to extend the code for 3D finite elements to consider multi-layer piezoelectric plates and out of plane forces.

Conclusion

Two MATLAB codes are proposed for topology optimization of piezoelectric actuators and energy harvesters. The codes are developed based on the finite element modeling of piezoelectric materials. The PEMAP-P as an extension of SIMP approaches is used for material interpolation scheme. Optimality criteria and method of moving asymptotes are used for optimization of element's density and polarization direction. Different parts of the codes are explained in detail to make the implementation and extension of the codes straightforward. Some basic and general examples are chosen to show the effectiveness of the codes. The aim of the codes is to help the students and newcomers in the field of topology optimization of smart materials in particular the piezoelectric material. While the codes have been and can be used for energy harvesting and actuation applications, perspective works include their application to optimization of piezoelectric sensors as well as piezoelectric sensors-actuators, also named as self-sensing (Bafumba [START_REF] Liseli | Optimal design of piezoelectric cantilevered actuators for charge-based self-sensing application[END_REF][START_REF] Rakotondrabe | Combining self-sensing with an unkowninput-observer to estimate the displacement, the force and the state in piezoelectric cantilevered actuator[END_REF][START_REF] Aljanaideh | Observer and robust h-inf control of a 2dof piezoelectric actuator equiped with self-measurement[END_REF]).

Fig. 1

 1 Fig. 1 Piezoelectric plate sandwiched between two electrodes. a Isometric view. b Side view

Fig. 2

 2 Fig. 2 Finite element discretization of design domain. Panel (a) is the numbering format inside the design domain. Panel (b) is the numbering format inside each element. Panel (c) is the parent element in the natural coordinates

Fig. 3

 3 Fig. 3 Topology optimization of a piezoelectric actuator (pusher) for different stiffness of the modeled spring. Panel a presents the specification for a practical use and panel b is the mechanical model for implementation in the finite element software. Panels c, e, g, i

Fig. 5

 5 Fig. 5 Topology optimization of a piezoelectric energy harvester under application of a lateral static force for different values of w j . Panel a presents the specification for a practical use and panel b is the result that can be obtained through classical compliance optimization. Panels

Fig. 6

 6 Fig. 6 Topology optimization of a piezoelectric energy harvester under application of a lateral harmonic force for different values of w j . Panel a presents the specification for a practical use and panels b, d, f,

  formulated the optimization problem to optimize the eigenmodes in addition to optimization of the eigenfrequency. It is worth mentioning that for these problem formulations, alternative material interpolation functions should be used to avoid the artificial local modes in the low-density regions (Pedersen 2000). The alternative material interpolation functions can be the one introduced by Huang et al. (2010) for the stiffness matrix or the RAMP interpolation function which is used by Nakasone and Silva (2010) in combination with the PEMAP-P interpolation function.

Fig. 7

 7 Fig. 7 Topology optimization of a piezoelectric energy harvester under application of 2 load cases for different values of w j and Ω. Panel a presents the specification for a practical use and panel b is the result that can be obtained through classical compliance optimization. Panels

Table 1

 1 Displacement amplification ratio of optimized actuators with respect to full plate

		Pusher		Gripper
		k s	Amplification	k s	Amplification
			ratio		ratio
	CASE (1)	1	0.76	1	3.95
	CASE (2)	0.1	0.99	0.1	11.42
	CASE (3)	0.01	1.54	0.01	34.04
	CASE (4)	0.005	2.75	0.003	60.11

Table 2

 2 Output energies of the optimized designs for different weighting factors

		w j	Π S	Π E	p uu	p uφ	p φφ	p P	Ω	x {0}	MMA move
	Lateral force
	Compliance	1	116.66	0.00	3	6	4	0	0	volfrac	1
	CASE (1)	1	99.90	1.25	3	6	4	1	0	volfrac	1
	CASE (2)	0.01	152.87	1.70	3	4	4	1	0	volfrac	1
	CASE (3)	0.005	741.09	7.77	3	4	4	1	0	volfrac	1
	CASE (4)	1	106.04	1.38	3	6	6	1	1 kHz	volfrac	1
	CASE (5)	0.01	123.15	1.43	3	6	6	1	1 kHz	volfrac	0.1
	CASE (6)	0.02	240.50	2.63	3	6	4	1	3.5 kHz	volfrac	0.1
	2 load case
	Compliance	1	106	0.03	3	6	6	0	0	1	1
	CASE (1)	0.005	168.54	1.20	3	6	6	1	0	1	1
	CASE (2)	1	107.88	0.86	3	6	6	1	1 kHz	1	0.1
	CASE (3)	0.005	163.43	1.04	3	6	4	1	1 kHz	1	0.1
	CASE (4)	0.05	120.06	0.95	3	6	6	1	1.5 kHz	1	0.1

Funding This work has been supported by the national CODE-TRACK project (ANR-17-CE05-0014-01, Control theory tools for optimal design of piezoelectric energy harvesters devoted to birds tracking devices). This work has also been partially supported by the Bourgogne Franche-Comté region project COMPACT.

Compliance with ethical standards

Conflict of interest

The authors declare that they have no conflict of interest.

Replication of results

All the results presented in this work can be reproduced with the MATLAB code available within the paper and by adding reference of this latter.

Appendix: Proof of normalization

Here the goal is to prove the normalization which is proposed by authors in Homayouni-Amlashi (2019) and Homayouni-Amlashi et al. (2020b). To normalize the global FEM equation of the piezoelectric plate

the first step of normalization is to factorize the values defined in (20) and ( 22)

It should be noted that since the responses of the system, i.e., U and Φ are unknown, their factorization values will be defined later.

With the help of factorization defined in (A.2), (A.1) can be rewritten as follows,

This is two linearly coupled equations. The second equation can be written as

Then, the Φ will be

The first linear equation from (A.3) can be written as

By substituting (A.5) to (A.6) and dividing the resulted equation by k 0 u 0 , one will have

Now, it is possible to define,

By substituting (A.8) to (A.7)

Now, the normalization factor can be defined as

Again, the (A.9) can be written as,

Now, it can be considered that

If one write (A.12) as Kφu Ũγ Kφφ Φ = 0 (A.13) then the following linearly coupled equation can be written with the help of (A.11) and (A.13)

which is the same as ( 23) of the paper. However, right now we have two equations for Φ in (A.12) and (A.5). By equating these two equations results in

We can then simplify both sides and replace the γ coefficient

We will find the following equation for φ 0 ,

But from (A.8) u 0 can be replaced and the final value for φ 0 will be

Now, by finding the values of γ , u 0 and φ 0 , the normalization of global FEM equations of piezoelectric material is finished. Actuation is the sub problem for the aforementioned normalization. Only ( 16) should be normalized. This equation with the help of the normalization factors defined in ( 20) and ( 24) can be written as

By dividing each side of equation by f 0 ,

Now by considering u 0 and φ 0 equal to values in (A.8) and (A.18), the normalized (21) will be obtained.

Kuu Ũ + Kuφ Φ = F

(A.21)