Topology optimization of piezoelectric bi-material actuators with velocity feedback control

Abstract

In recent years, the new technologies and discoveries on manufacturing materials have encouraged researchers to investigate the appearance of material properties that are not naturally available. Materials featuring a specific stiffness, or structures that combine non-structural and structural functions are applied in the aerospace, electronics and medical industry fields. Particularly, structures designed for dynamic actuation with reduced vibration response are the focus of this work. The bi-material and multifunctional concepts are considered for the design of a controlled piezoelectric actuator with vibration suppression by means of the topology optimization method (TOM). The bi-material piezoelectric actuator (BPEA) has its metallic host layer designed by the TOM, which defines the structural function, and the electric function is given by two piezo-ceramic layers that act as a sensor and an actuator coupled with a constant gain active velocity feedback control (AVFC). The AVFC, provided by the piezoelectric layers, affects the structural damping of the system through the velocity state variables readings in time domain. The dynamic equation analyzed throughout the optimization procedure is fully elaborated and implemented. The dynamic response for the rectangular four-noded finite element analysis is obtained by the Newmark’s time-integration method, which is applied to the physical and the adjoint systems, given that the adjoint formulation is needed for the sensitivity analysis. A gradient-based optimization method is applied to minimize the displacement energy output measured at a predefined degree-of-freedom of the BPEA when a transient mechanical load is applied. Results are obtained for different control gain values to evaluate their influence on the final topology.

Keywords

topology optimization method / bi-material piezoactuator / active velocity feedback control / time-domain transient analysis / host structure design / vibration suppression

1 Introduction

In this work, the topology optimization (TO) extended to a bi-material metallic design of piezoactuators in time-domain analysis, seeks the structural vibration suppression needed in modern electronic devices applications that require lightweight, fast-response and fast stabilizing actuators [1]. The bi-material piezoelectric actuator (BPEA) structure designed in this work, features non-structural functions such as the piezoelectric effect, for sensing and actuation of a controlled transducer, and optimized structural functions such as stiffness and damping of a bi-metallic host layer. Some recent applications for multifunction structures, such as the described transducer, and multifunction materials, are in morphing aircraft wings and structurally integrated electronic components [2]. For the design of the BPEA with an active velocity feedback control (AVFC) for vibration attenuation, the topology optimization method (TOM) is the most suited technique once it allows a continuous material properties variation within the microstructure according to the objective of the design.

The design of controlled piezoactuators using the TOM has been vastly explored for transient applications in frequency domain. The structural design for a spacecraft subjected to an extreme thermal environment has been studied in Ref. [3], where a membrane thermal surface load is carried out, which assumes the form of a time-harmonic external load. The thermal force is design-dependent since it is related to a thermal stress coefficient. Piezoelectric actuators have been studied in frequency domain by means of the model reduction method in Ref. [4], where piezoelectric sensors and actuator integrated with an AVFC are distributed over a host layer aiming vibration suppression, and in Ref. [5], that also distributes piezoelectric material but in order to maximize the controllability for a given vibration mode. In this last case, the structural model is written into the state-space representation and limites the control spillover effects, caused by the residual modes, by including a spillover constraint. Yet in Ref. [6], authors optimized the surface electrode distribution over piezoelectric sensors and actuators which are integrated with an AVFC for reducing the sound radiation in an unbounded acoustic domain. However, literature survey shows that fewer works involve the TOM apllied to the design of a piezoelectric structure under a time-domain transient load.

The TOM procedure associated with time-domain analysis is a sensitive problem towards convergence achievement in terms of the optimization problem parameters setup. The number of time steps used in time integration for the dynamic FE (FE) analysis and the time varying rate used for the dynamic input load are examples of parameters that have a significant influence on the final topological result [7]. For this reason, when it comes to the design of piezoelectric actuators for time domain dynamic applications, researchers mostly deal with plate structures, given a fixed host layer structure, and optimize the piezoelectric material distribution for a transient objective function [8]. However, inspired by the flextensional piezoelectric actuators, the objective of this work is to optimize the host layer topology by taking advantage of the new structural property given by two metals optimally distributed withing the domain.

The active vibration control scheme in piezoelectric systems, that is stablished by the extra voltage or electric charge supplied to the piezoceramic material [9], and that is applied to the design of the controlled BPEA therein proposed, has been studied for flexible plates [10–12], for cylindrical laminated composite shells [13], for beams [14,15], and for precision positioning of hard disk drives [16,17]. Tzou and Tseng [10] implemented the constant-gain feedback control and the constant-amplitude feedback control. Wang et al. [11] also studied a time-domain formulation with the velocity feedback control, and evaluated the system stability according to its piezoceramics placement. Besides the velocity feedback, linear quadratic regulator (LQR) and proportion integration differentiation (PID) control schemes have been studied for different excitation types applyed to flexible structures [18]. Among those works, the predefined placement of sensors and actuators associated with a controller have resulted in vibration attenuations of the system.

Therefore, this paper treats the design of an optimal topology for the bi-material host layer in time-domain analysis for fixed ceramic layers and a control gain. The proposed structure involves the clamped-free host layer to which a transient disturbance is applied. The piezoelectric sensor layer embedded to one of the surfaces monitors the structural deformation and generates an electrical charge vector ${\mathit{Q}}_{\mathrm{p}}^{\mathrm{s}}$, that differentiated in time, ${\dot{\mathit{Q}}}_{\mathrm{p}}^{\mathrm{s}}$ amplified by the control gain K and the resulting voltage ${\mathbf{\Phi }}_{\mathrm{p}}^{\mathrm{a}}$ is fed back into the actuator piezoelectric layer attached to the opposite surface of the host structure, as illustrated in Fig. 1.

Fig.1 Bi-material piezoactuator

Full size|PPT slide

To illustrate the method, the bidimensional domain is discretized in bilinear FEs with three degrees-of-freedom per node. The TOM implementation is based on a density material model, the solid isotropic material with penalization (SIMP), and on the projection filter proposed by Ref. [19]. The sensitivity analysis is based on the adjoint method [20] and the optimization problem is solved by the linear programming (LP) algorithm. The AVFC influence in the optimized topology results for vibration suppression is analyzed.

This paper is organized as follows. In Section 2, the evolving steps to obtain the FE formulation of the piezoelectric actuator equation with the AVFC is described. In Section 3, the TO problem based on the time-domain analysis is stated, by explaining the material model chosen and deriving the adjoint sensitivity formulation. In Section 4, a flowchart illustrating the step-wise iterations for the problem solved is presented. In Section 5, the BPEA design for some gain values are obtained and discussed. The results are analyzed based on different feedback control gains towards a more efficient vibration suppression of the optimized piezoelectric structure. In Section 6, some conclusions and the perspectives for future studies are inferred.

2 Controlled FE model for a BPEA

2.1 Dynamic piezoelectric equations

The direct and converse piezoelectric effects considered for the sensing and actuation of the BPEA, respectively, are given by the constitutive equations:

(1)$\{\begin{array}{l}\mathit{T}={\mathit{c}}^{\mathrm{E}}\mathit{S}-\mathit{e}\mathit{E}\\ \mathit{D}=\mathit{e}\mathit{S}+{\mathit{\epsilon }}^{\mathrm{S}}\mathit{E}\end{array},$

where $\mathit{T}$ and $\mathit{S}$ are the stress and strain tensors, respectively, and $\mathit{D}$ is the electric displacement vector. The piezoelectric elastic matrix ${\mathit{c}}^E$, evaluated under a constant electric field $\mathit{E}$, the dielectric tensor ${\mathit{\epsilon }}^{\mathrm{S}}$, evaluated under constant strain $\mathit{S}$, and the coupling piezoelectric tensor $\mathit{e}$, are the piezoelectric material constituents.

The system dynamics is approximated by a 2D-solid bilinear FE where the displacement field $\mathit{u}$ and the electric potential $\varphi$ are interpolated by means of the shape functions ${\mathit{N}}_u$ and ${\mathit{N}}_{\varphi }$, respectively. The displacement field $\mathit{U}={\left\{\begin{array}{ll}\mathit{u}& \mathit{\Phi }\end{array}\right\}}^{\mathrm{T}}$ is given by two mechanical components $\mathit{u}={\left[\begin{array}{ll}{\mathit{u}}_x& {\mathit{u}}_y\end{array}\right]}^{\mathrm{T}}$ for deformations in directions $x$ and $y$, and an electric component $\Phi =\left[\begin{array}{ll}{\Phi }_{\mathrm{f}}& {\Phi }_{\mathrm{p}}\end{array}\right]$ subdivided into free voltages, ${\Phi }_{\mathrm{f}}$, and prescribed voltages, ${\Phi }_{\mathrm{p}}$. Therefore, the dynamic FE equilibrium equations for the piezoelectric system are written [21]:

(2)$\left[\begin{array}{lll}{\mathit{M}}_{uu}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right]\left\{\begin{array}{l}\ddot{u}\\ {\ddot{\Phi }}_{\mathrm{f}}\\ {\ddot{\Phi }}_{\mathrm{p}}\end{array}\right\}+\left[\begin{array}{lll}{\mathit{C}}_{uu}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}& \mathbf{0}\end{array}\right]\left\{\begin{array}{l}\dot{u}\\ {\dot{\Phi }}_{\mathrm{f}}\\ {\dot{\Phi }}_{\mathrm{p}}\end{array}\right\}+\left[\begin{array}{lll}{\mathit{K}}_{uu}& {\mathit{K}}_{u{\varphi }_{\mathrm{f}}}& {\mathit{K}}_{u{\varphi }_{\mathrm{p}}}\\ {\mathit{K}}_{u{\varphi }_{\mathrm{f}}}^{\mathrm{T}}& -{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{f}}}& -{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{p}}}\\ {\mathit{K}}_{u{\varphi }_{\mathrm{p}}}^{\mathrm{T}}& -{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{p}}}^{\mathrm{T}}& -{\mathit{K}}_{{\varphi }_p{\varphi }_{\mathrm{p}}}\end{array}\right]\left\{\begin{array}{l}u\\ {\Phi }_{\mathrm{f}}\\ {\Phi }_{\mathrm{p}}\end{array}\right\}=\left\{\begin{array}{l}\mathit{F}\\ {\mathit{Q}}_{\mathrm{f}}\\ {\mathit{Q}}_{\mathrm{p}}\end{array}\right\},$

where the global mass matrix ${\mathit{M}}_{uu}$, stiffness matrix ${\mathit{K}}_{uu}$, piezoelectric coupling matrix, ${\mathit{K}}_{u\varphi }={\mathit{K}}_{\varphi u}$, dielectric permittivity matrix ${\mathit{K}}_{\varphi \varphi }$, external nodal force vector ${\mathit{K}}_u$ , where ${\mathit{f}}_{\mathrm{b}}$ is the body force vector and $\overline{\mathit{t}}$ is the surface traction vector, and the free charge vector ${\mathit{Q}}_{\mathrm{p}}$, where $\bar{q}$ is the free charge density per unit surface area, are defined from the elementary matrices (superscript e) as bellow:

(3)${\mathit{M}}_{uu}^e={\displaystyle {\int }_{{\mathrm{\Omega }}_e}\rho N_u^{\mathrm{T}}{N}_u^{}}\hspace{0.17em}\mathrm{d}{\mathrm{\Omega }}_e,\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{\mathit{M}}_{uu}={\displaystyle \sum _e\hspace{0.17em}{\mathit{M}}^e},$

(4)${\mathit{K}}_{uu}^e={\displaystyle {\int }_{{\mathrm{\Omega }}_e}{\mathit{B}}_u^{\mathrm{T}}{c}^{\mathrm{E}}{\mathit{B}}_u^{}\mathrm{d}{\mathrm{\Omega }}_e},\hspace{0.17em}\hspace{0.17em}{\mathit{K}}_{uu}={\displaystyle \sum _e{\mathit{K}}_{uu}^e},$

(4)

(5)${\mathit{K}}_{u\varphi }^e={\displaystyle {\int }_{{\mathrm{\Omega }}_e}{\mathit{B}}_u^{\mathrm{T}}}{e}^{\mathrm{T}}{\mathit{B}}_{\varphi }\hspace{0.17em}d{\mathrm{\Omega }}_e,\hspace{0.17em}\hspace{0.17em}{\mathit{K}}_{u\varphi }={\displaystyle \sum _e{\mathit{K}}_{u\varphi }^e},$

(5)

(6)$\hspace{0.17em}\hspace{0.17em}{\mathit{K}}_{\varphi \varphi }^e=-{\displaystyle {\int }_{{\mathrm{\Omega }}_e}{\mathit{B}}_{\varphi }^{\mathrm{T}}}{\epsilon }^{\mathrm{S}}{\mathit{B}}_{\varphi }^{}\hspace{0.17em}\mathrm{d}{\mathrm{\Omega }}_e,\hspace{0.17em}\hspace{0.17em}{\mathit{K}}_{\varphi \varphi }={\displaystyle \sum _e{\mathit{K}}_{\varphi \varphi }^e},$

(6)

(7)${{\displaystyle \mathit{F}}}_u^e={\displaystyle {\int }_{{{\displaystyle \mathrm{\Omega }}}_e}{{\displaystyle N}}_u^{\mathrm{T}}}{{\displaystyle \mathit{F}}}_b\mathrm{d}{{\displaystyle \mathrm{\Omega }}}_e+{\displaystyle {\int }_{{{\displaystyle \mathrm{\Gamma }}}_T}{{\displaystyle N}}_u^{\mathrm{T}}}\overline{\mathit{t}}\mathrm{d}{{\displaystyle \mathrm{\Gamma }}}_T,$

(7)

(8)${\mathit{Q}}_{\mathrm{p}}^e=-{\displaystyle {\int }_{{\mathrm{\Gamma }}_D}{N}_{\varphi }^{\mathrm{T}}}\overline{q}\hspace{0.17em}\mathrm{d}{\mathrm{\Gamma }}_D,$

where the boundary surface is composed by the prescribed traction region, ${{\displaystyle \mathrm{\Gamma }}}_T$ , and the electroded region, ${{\displaystyle \mathrm{\Gamma }}}_D$ . The subscript $\mathit{u}\mathit{u}$ denotes structural states, the subscript $\varphi$ denotes electrical states, and ${{\displaystyle \mathit{B}}}_u$ and ${{\displaystyle \mathit{B}}}_{\varphi }$ are, respectively, the mechanical and electrical shape functions derivatives [22]. The proportional damping matrix ${\mathit{C}}_{uu}$ is defined by

(9)${\mathit{C}}_{uu}=\alpha {\mathit{M}}_{uu}+\beta {\mathit{K}}_{uu},$

where $\alpha$ and $\beta$ are Rayleigh’s coefficients, that are non-negative and small in the sense that the effective damping coefficients is $\ll \mathrm{1}$ [23], and are considered to be design-independent [24,25].

2.2 Active velocity feedback control

The problem of vibration suppression on smart structures is addressed by performing the velocity feedback control, which is derived from the electrical charge readings on the sensor layer ${\mathit{Q}}_{\mathrm{p}}^{\mathrm{s}}$, by imposing null voltage on the sensor electrodes ${\Phi }_{\mathrm{p}}^{\mathrm{s}}=\mathrm{0}$ when the system is under a mechanical transient load. Knowing that the electric charges at the piezoelectric internal nodes ${\mathit{Q}}_{\mathrm{f}}$ are null and that the actuator and sensor electrodes do not share degrees of freedom, ${\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_{\mathrm{p}}^{\mathrm{a}}}=\mathrm{0}$, the expanded Eq. (2) is as follows:

(10)${\mathit{M}}_{uu}\ddot{\mathit{u}}+{\mathit{C}}_{uu}\dot{u}+{\mathit{K}}_{uu}u+{\mathit{K}}_{u{\varphi }_{\mathrm{f}}}{\Phi }_{\mathrm{f}}+{\mathit{K}}_{u{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\Phi }_{\mathrm{p}}^{\mathrm{a}}=\mathit{F},$

(11)${\mathit{K}}_{{\varphi }_{f}u}\mathit{u}-{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{f}}}{\Phi }_{\mathrm{f}}-{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\Phi }_{\mathrm{p}}^{\mathrm{a}}=\mathrm{0},$

(12)${\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u}\mathit{u}-{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_f}{\Phi }_{\mathrm{f}}={\mathit{Q}}_{\mathrm{p}}^{\mathrm{s}}$

(13)${\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}u}\mathit{u}-{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}{\varphi }_{\mathrm{f}}}{\Phi }_{\mathrm{f}}-{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\Phi }_{\mathrm{p}}^{\mathrm{a}}={\mathit{Q}}_{\mathrm{p}}^{\mathrm{a}}.$

The derivative of the electrical charge Eq. (12) results in the nodal electric current vector ${\mathit{I}}^{\mathrm{s}}$ on the external sensor electrode:

(14)${\mathit{I}}^{\mathrm{s}}={\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u}\dot{\mathit{u}}-{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_{\mathrm{f}}}{\dot{\Phi }}_{\mathrm{f}},$

that can be added by pre-multiplying ${\mathit{I}}^{\mathrm{s}}$ by the unity vector ${\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}=\left\{\begin{array}{lll}\mathrm{1}& \cdots & \mathrm{1}\end{array}\right\}$, and the measured current in the top electrode of the sensor layer is obtained:

(15)$I^{\mathrm{s}}={\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{I}}^{\mathrm{s}}={\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}\left({\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u}\dot{u}-{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_{\mathrm{f}}}{\dot{\Phi }}_{\mathrm{f}}\right).$

To establish a velocity feedback control law, the actuator input voltage ${\Phi }_{\mathrm{p}}^{\mathrm{a}}$ is defined as a current amplifier:

(16)${\Phi }_{\mathrm{p}}^{\mathrm{a}}=K{I}^{\mathrm{s}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}},$

where ${\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}={\left\{\begin{array}{lll}\mathrm{1}& \cdots & \mathrm{1}\end{array}\right\}}^{\mathrm{T}}$ and K is the feedback gain.

The nodal voltage input at the actuator layer, Eq. (16), is then rewritten by substituting Eq. (15) into Eq. (16):

(17)${\mathrm{\Phi }}_{\mathrm{p}}^{\mathrm{a}}=K{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}u}\dot{\mathit{u}}-K{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^s}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^s{\varphi }_{\mathrm{f}}}{\dot{\Phi }}_{\mathrm{f}}.$

Now, substituting Eq. (17) into Eqs. (10), (11) and (13), we have

(18)${\mathit{M}}_{uu}\ddot{\mathit{u}}+({\mathit{C}}_{uu}+K{\mathit{K}}_{u{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u})\dot{\mathit{u}}-K{\mathit{K}}_{u{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_f}{\dot{\mathrm{\Phi }}}_{\mathrm{f}}+{\mathit{K}}_{uu}\mathit{u}+{\mathit{K}}_{u{\varphi }_{\mathrm{f}}}{\Phi }_{\mathrm{f}}=\mathit{F},$

(19)$-K{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}u}\dot{\mathit{u}}+K{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}{\varphi }_{\mathrm{f}}}{\dot{\Phi }}_{\mathrm{f}}+{\mathit{K}}_{{\varphi }_{\mathrm{f}}u}u-{\mathit{K}}_{{\varphi }_f{\varphi }_{\mathrm{f}}}{\mathbf{\Phi }}_{\mathrm{f}}=\mathrm{0},$

(20)$K{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u}\dot{\mathit{u}}-K{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_{\mathrm{f}}}{\dot{\Phi }}_{\mathrm{f}}+{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}u}\mathit{u}-{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}{\varphi }_{\mathrm{f}}}{\Phi }_{\mathrm{f}}={\mathit{Q}}_{\mathrm{p}}^{\mathrm{a}},$

and the coupled system to be solved is given by Eqs. (18) and (19):

(21)$\left[\begin{array}{ll}{\mathit{M}}_{uu}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right]\left\{\begin{array}{l}\ddot{\mathit{u}}\\ {\ddot{\Phi }}_f\end{array}\right\}+\left[\begin{array}{ll}{\mathit{C}}_{uu}+{\mathit{C}}_{\mathrm{11}}& {\mathit{C}}_{\mathrm{12}}\\ {\mathit{C}}_{\mathrm{21}}& {\mathit{C}}_{\mathrm{22}}\end{array}\right]\left\{\begin{array}{l}\dot{\mathit{u}}\\ {\dot{\Phi }}_{\mathrm{f}}\end{array}\right\}+\left[\begin{array}{ll}{\mathit{K}}_{uu}& {\mathit{K}}_{u{\varphi }_{\mathrm{f}}}\\ {\mathit{K}}_{u{\varphi }_{\mathrm{f}}}^{\mathrm{T}}& -{\mathit{K}}_{{\varphi }_{\mathrm{f}}{\varphi }_{\mathrm{f}}}\end{array}\right]\left\{\begin{array}{l}\mathit{u}\\ {\Phi }_{\mathrm{f}}\end{array}\right\}=\left\{\begin{array}{l}\mathit{F}\\ \mathbf{0}\end{array}\right\},$

where ${\mathit{C}}_{uu}$ is given by Eq. (9),

(22)${\mathit{C}}_{\mathrm{11}}=K{\mathit{K}}_{u{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u},$

(23)${\mathit{C}}_{12}=K{\mathit{K}}_{u{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}\mathrm{f}},$

(24)${\mathit{C}}_{\mathrm{21}}=-K{\mathit{K}}_{{\varphi }_{\mathrm{f}}{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}u},$

(25)${\mathit{C}}_{\mathrm{22}}=K{\mathit{K}}_{{\varphi }_{\mathrm{f}}{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{a}}}{\mathit{I}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}}{\mathit{K}}_{{\varphi }_{\mathrm{p}}^{\mathrm{s}}{\varphi }_{\mathrm{f}}}.$

Equation (13) can be used for verification purposes.

For the sake of simplifying the representation of the dynamic system above, Eq. (21) is rewritten as below:

(26)$\mathit{M}\ddot{\mathit{U}}+{\mathit{C}}_{\mathrm{fb}}\dot{\mathit{U}}+\mathit{K}\mathit{U}=\mathit{F}(t),$

for

(27)$\mathit{M}=\left[\begin{array}{ll}{\mathit{M}}_{uu}& \mathbf{0}\\ \mathbf{0}& \mathbf{0}\end{array}\right],$

(28)${\mathit{C}}_{\mathrm{fb}}=\left[\begin{array}{ll}{\mathit{C}}_{uu}+{\mathit{C}}_{\mathrm{11}}& {\mathit{C}}_{\mathrm{12}}\\ {\mathit{C}}_{\mathrm{21}}& {\mathit{C}}_{\mathrm{22}}\end{array}\right],$

(29)$\mathit{K}=\left[\begin{array}{ll}{\mathit{K}}_{uu}& {\mathit{K}}_{u{\varphi }_{\mathrm{f}}}\\ {\mathit{K}}_{u{\varphi }_{\mathrm{f}}}^{\mathrm{T}}& -{\mathit{K}}_{{\varphi }_{\mathrm{f}}{\varphi }_{\mathrm{f}}}\end{array}\right],$

(30)$\mathit{U}=\left\{\begin{array}{l}\mathit{u}\\ {\Phi }_{\mathrm{f}}\end{array}\right\},$

(31)$\mathit{F}=\left\{\begin{array}{l}{\mathit{F}}_u\\ \mathbf{0}\end{array}\right\},$

(32)${\mathit{F}}_u=\left\{\begin{array}{l}\mathbf{0}\\ ⋮\\ {\mathit{F}}_{\mathrm{dof}}\\ ⋮\\ \mathbf{0}\end{array}\right\}.$

The modified damping matrix ${\mathit{C}}_{\mathrm{fb}}$, Eq. (28), contains extra terms given the velocity readings from the piezoelectric sensor layer, as shown by Eq. (17), that provides an active damping effect to the structure.

3 Topology optimization problem

3.1 Problem formulation

The TOM employs a material model concept [26] to define the microstructure property by distributing different materials within a design domain, aiming to extremize a cost function, and uses the FE method for systematic structure analysis. Targeting the vibrabration suppression of a BPEA, the objective function defined in this work considers the minimization of a energy functional involving the displacement of a predefined degree-of-freedom as its squared integral spanned over a time interval $[\mathrm{0},t_{\mathrm{f}}]$ [8],

(33)$f(\mathit{\rho })={\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{u}_{\mathrm{dof}}^{\mathrm{2}}}\hspace{0.17em}\mathrm{d}t,$

where the $\mathit{\rho }$-vector elements are the design variables, and $u_{\mathrm{dof}}^{\mathrm{2}}={\mathit{U}}^{\mathrm{T}}\mathit{B}\mathit{U}$, where $\mathit{B}={\mathit{B}}^{\mathrm{T}}$ is

(34)$\mathit{B}=\left[\begin{array}{lllllll}\mathrm{0}& \dots & \dots & \dots & \dots & \dots & \mathrm{0}\\ ⋮& \ddots & & & & & ⋮\\ ⋮& & \mathrm{0}& & & & ⋮\\ ⋮& & & {\mathrm{1}}_{\mathrm{dof}}& & & ⋮\\ ⋮& & & & \mathrm{0}& & ⋮\\ ⋮& & & & & \ddots & ⋮\\ \mathrm{0}& \dots & \dots & \dots & \dots & \dots & \mathrm{0}\end{array}\right].$

As for the initial conditions, this dynamic problem depends on the magnitude of the input load. For a point-wise mechanical load with the profile shown in Fig. 2, the initial displacement vector ${\mathit{U}}_{\mathrm{0}}$ is given by the static analysis $\mathit{K}{\mathit{U}}_{\mathrm{0}}=\mathit{F}(\mathrm{0})$ while the initial velocity vector ${\mathit{V}}_{\mathrm{0}}$ is null. A volume constraint $V_{\max }$ limits the distribution of one of the materials. The optimization problem is therefore stated as bellow:

(35)$\begin{array}{l}\underset{\rho }{\min }\hspace{0.17em}\hspace{0.17em}f(\mathit{\rho })={\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{u}_{\mathrm{dof}}^{\mathrm{2}}\hspace{0.17em}}\mathrm{d}t,\\ \mathrm{s}.\mathrm{t}.\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\{\begin{array}{l}\mathit{M}(\mathit{\rho })\ddot{\mathit{U}}+C_{\mathrm{fb}}(\mathit{\rho })\dot{\mathit{U}}+K(\mathit{\rho })\mathit{U}=\mathit{F}(t),\hfill \\ \mathit{K}{\mathit{U}}_{\mathrm{0}}=\mathit{F}(\mathrm{0}),\hfill \\ \dot{\mathit{U}}{|}_{t=\mathrm{0}}={\mathbf{V}}_{\mathrm{0}},\hfill \\ \mathit{U}{|}_{t=\mathrm{0}}={\mathit{U}}_{\mathrm{0}},\hfill \\ {\displaystyle \sum _{n=\mathrm{1}}^N{\rho }_n}{V}_n\le V_{\max },\hfill \\ {\mathrm{10}}^{-\mathrm{15}}\le {\rho }_n\le \mathrm{1.}\hfill \end{array}\end{array}$

Fig.2 Transient load profile for the example proposed

Full size|PPT slide

As stated in optimizatioin problem Eq. (35), its constraints also includes the dynamic equilibrium Eq. (26), where $\mathit{F}(t)=\mathit{F}y(t)$, and boundary values for each element design variable ${\mathit{\rho }}_n$ of a total of $N$ elements. The metallic materials distribution is such that it minimizes the cost function Eq. (33).

3.2 Material model

In this work, for the transition from one metallic material to the other, it is adopted the SIMP, which has been proved to be robust when applied to linear elastic problems [27–29]. It allows the properties at each element to assume intermediate values, and though intermediate properties of the bi-material, since the gray scale is allowed. Smooth increments to the penalization factor $p$ are applied at each iteration of the TOM what is denominated the continuation approach. This procedure prevents a premature convergence to a local minima [30].

Based on SIMP, the elasticity matrix, ${\mathit{c}}^e(\rho )$, and the specific weight, $\zeta$, specified for the elastic domain are functions of the pseudo-density $\rho$ defined for each FE $e$:

(36)${\zeta }^e(\rho )={\rho }^{p_{\mathrm{1}}}(x_{\mathrm{c}},y_{\mathrm{c}}){\zeta }_{\mathrm{mat1}}+[\mathrm{1}-{\rho }^{p_{\mathrm{1}}}(x_{\mathrm{c}},y_{\mathrm{c}})]{\zeta }_{\mathrm{mat2}},$

(37)$c^e(\rho )={\rho }^{p_{\mathrm{2}}}(x_{\mathrm{c}},y_{\mathrm{c}})c_{\mathrm{mat1}}+[\mathrm{1}-{\rho }^{p_{\mathrm{2}}}(x_{\mathrm{c}},y_{\mathrm{c}})]c_{\mathrm{mat2}},$

where $(x_{\mathrm{c}},y_{\mathrm{c}})$ is the centroid cartesian coordinate pair of each element, and ${\mathit{c}}_{\mathrm{mat1}}$ and ${\mathit{c}}_{\mathrm{mat2}}$ are the elastic tensors for the first and second materials, respectively.

3.3 Sensitivity analysis

Considering the design problem as an optimization problem in the pseudo-densities only, the sensitivity analysis consists on finding the objective function derivatives with respect to the vector of design variables $\mathit{\rho }$. These derivatives, or the sensitivity analysis is given by the adjoint method which involves an adjoint vector $\mathit{\lambda }(t)$ that is chosen as the solution to the adjoint equations derived from a Lagrangian function $L(\mathit{U},\mathit{\lambda })$ [20]:

(38)$L(\mathit{U},\mathit{\lambda })={\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\mathit{U}}^{T}B\mathit{U}\hspace{0.17em}\mathrm{d}t}+{\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\mathit{\lambda }}^{\mathrm{T}}(t)\left[\mathit{M}(\mathit{\rho })\ddot{\mathit{U}}+{\mathit{C}}_{\mathrm{fb}}(\mathit{\rho })\dot{\mathit{U}}+\mathit{K}(\mathit{\rho })\mathit{U}-\mathit{F}(t)\right]\hspace{0.17em}\hspace{0.17em}\mathrm{d}t}.$

Taking the derivative of the Lagrangian Eq. (38) with respect to the design variables $\mathit{\rho }$ we get the expression

(39)$\begin{array}{l}{\displaystyle \frac{\partial L(\mathit{U},\lambda )}{\partial {\rho }_e}}={\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\displaystyle \frac{\partial }{\partial {\rho }_e}}({\mathit{U}}^{\mathrm{T}}\mathit{B}\mathit{U})\hspace{0.17em}\hspace{0.17em}\mathrm{d}t}\\ +{\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\mathit{\lambda }}^{\mathrm{T}}(t)\left({\displaystyle \frac{\partial \mathit{M}(\mathit{\rho })}{\partial {\rho }_e}}\ddot{\mathit{U}}+{\displaystyle \frac{\partial {\mathit{C}}_{\mathrm{fb}}(\mathit{\rho })}{\partial {\rho }_e}}\dot{\mathit{U}}+{\displaystyle \frac{\partial \mathit{K}(\mathit{\rho })}{\partial {\rho }_e}}\mathit{U}-{\displaystyle \frac{\partial \mathit{F}(t)}{\partial {\rho }_e}}\right)\hspace{0.17em}\hspace{0.17em}\mathrm{d}t}\\ +{\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\mathit{\lambda }}^{\mathrm{T}}(t)\left(\mathit{M}(\mathit{\rho }){\displaystyle \frac{\partial \ddot{\mathit{U}}}{\partial {\rho }_e}}+{\mathit{C}}_{\mathrm{fb}}(\mathit{\rho }){\displaystyle \frac{\partial \dot{\mathit{U}}}{\partial {\rho }_e}}+\mathit{K}(\mathit{\rho }){\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}\right)\hspace{0.17em}\hspace{0.17em}\mathrm{d}t},\end{array}$

that is solved by integration by parts:

(40)

(41)$$\begin{gathered} {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}\mathit{M}(\mathit{\rho })\lambda (t){\displaystyle \frac{\partial \dot{\mathit{U}}}{\partial {\rho }_e}}\mathrm{d}t=\mathit{M}({\mathit{\rho }}_e)\frac{\mathrm{d}}{\mathrm{d}t}\left({\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}\right){\lambda (t)|}_{\mathrm{0}}^{t_{\mathrm{f}}}-{\int }_{\mathrm{0}}^{t_{\mathrm{f}}}\mathit{M}(\mathit{\rho })\dot{\mathit{\lambda }}(t)\frac{\mathrm{d}}{\mathrm{d}t}\left({\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}\right)\mathrm{d}t \\ =\mathit{M}(\mathit{\rho })\lambda (t){\frac{\mathrm{d}}{\mathrm{d}t}\left({\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}\right)|}_{\mathrm{0}}^{t_{\mathrm{f}}}-\mathit{M}(\mathit{\rho })\dot{\mathit{\lambda }}(t){{\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}|}_{\mathrm{0}}^{t_{\mathrm{f}}}+{\int }_{\mathrm{0}}^{t_{\mathrm{f}}}\mathit{M}(\mathit{\rho })\ddot{\mathit{\lambda }}(t){\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}\mathrm{d}t, \end{gathered}$$

(42)

Rearranging Eqs. (40)–(42), Eq. (39) is then rewritten:

(43)$\begin{array}{l}{\displaystyle \frac{\partial L(\mathit{U},\mathit{\lambda })}{\partial {\rho }_e}}={\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\lambda }^{\mathrm{T}}(t)\left({\displaystyle \frac{\partial \mathit{M}(\mathit{\rho })}{\partial {\rho }_e}}\ddot{\mathit{U}}+{\displaystyle \frac{\partial {\mathit{C}}_{\mathrm{fb}}(\mathit{\rho })}{\partial {\rho }_e}}\dot{\mathit{U}}+{\displaystyle \frac{\partial \mathit{K}(\mathit{\rho })}{\partial {\rho }_e}}\mathit{U}-{\displaystyle \frac{\partial \mathit{F}(t)}{\partial {\rho }_e}}\right)\hspace{0.17em}\hspace{0.17em}\mathrm{d}t}\\ +{\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}\left(\mathit{M}(\mathit{\rho }){\ddot{\mathit{\lambda }}}^{\mathrm{T}}(t)-{\mathit{C}}_{\mathrm{fb}}(\mathit{\rho }){\dot{\mathit{\lambda }}}^{\mathrm{T}}(t)+K(\mathit{\rho }){\mathit{\lambda }}^{\mathrm{T}}(t)+\mathit{B}\mathit{U}(t)+{\mathit{U}}^{\mathrm{T}}(t)\mathit{B}\right)}{\displaystyle \frac{\partial \mathit{U}}{\partial {\rho }_e}}\hspace{0.17em}\mathrm{d}t\\ +{\displaystyle \frac{\partial U}{\partial {\rho }_e}}{\left|-\mathit{M}(\mathit{\rho }){\dot{\mathit{\lambda }}}^{\mathrm{T}}(t)+{\mathit{C}}_{\mathrm{fb}}(\mathit{\rho }){\mathit{\lambda }}^{\mathrm{T}}(t)\right|}_{\mathrm{0}}^{t_{\mathrm{f}}}+{\displaystyle \frac{\partial \dot{\mathit{U}}}{\partial {\rho }_e}}{\left|M(\mathit{\rho }){\mathit{\lambda }}^{\mathrm{T}}(t)\right|}_{\mathrm{0}}^{t_{\mathrm{f}}},\end{array}$

To eliminate the state variables derivatives ${\displaystyle \frac{\partial \mathit{U}(t_{\mathrm{f}})}{\partial {\rho }_e}}$ and ${\displaystyle \frac{\partial \dot{\mathit{U}}(t_{\mathrm{f}})}{\partial {\rho }_e}}$ from Eq. (43), we define the adjoint variable ${\lambda }^T$ such that it annihilates the expression in square brackets of the second term given the initial conditions in the last two terms, resulting in the system Eq. (43). The derived adjoint ordinary differential equation is similar to the primal problem, except for a sign change on the adjoint velocity term ${\mathit{C}}_{\mathrm{fb}}\dot{\mathit{\lambda }}$. Given that $\mathit{B}\mathit{U}(t)+{\mathit{U}}^{\mathrm{T}}(t)\mathit{B}=\mathrm{2}{u}_{\mathrm{dof}}(t)$, the adjoint system that needs to be solved is:

(44)$\{\begin{array}{l}\mathit{M}{\ddot{\mathit{\lambda }}}^{\mathrm{T}}(t)-C_{\mathrm{fb}}{\dot{\mathit{\lambda }}}^{\mathrm{T}}(t)+\mathit{K}{\mathit{\lambda }}^{\mathrm{T}}(t)=-\mathrm{2}{u}_{\mathrm{dof}}(t),\hfill \\ {\mathit{\lambda }}^{\mathrm{T}}(t_{\mathrm{f}})=\mathbf{0},\hfill \\ {\dot{\mathit{\lambda }}}^{\mathrm{T}}(t_{\mathrm{f}})=\mathbf{0}.\hfill \end{array}$

In order to redefine the system Eq. (44) as an initial value primal problem, we apply the change of variables $\tau (t)=t_f-t$, and the adjoint system to be solved is as follows:

(45)$\{\begin{array}{l}\mathit{M}{\ddot{\Lambda }}^{\mathrm{T}}(\tau (t))+{\mathit{C}}_{\mathrm{fb}}{\dot{\Lambda }}^{\mathrm{T}}(\tau (t))+\mathit{K}{\Lambda }^{\mathrm{T}}(\tau (t))=-\mathrm{2}{u}_{\mathrm{dof}}(\tau (t)),\hfill \\ {\mathbf{\Lambda }}^{\mathrm{T}}(\mathrm{0})=\mathbf{0},\hfill \\ {\dot{\Lambda }}^{\mathrm{T}}(\mathrm{0})=\mathbf{0},\hfill \end{array}$

where $\Lambda (\tau (t))=\mathbf{\lambda }(t)$ and $\dot{\Lambda }(\tau (t))=-\dot{\mathit{\lambda }}(t)$. Therefore, the sensitivity expression for the transient problem reduces to

(46)${\displaystyle \frac{\partial L(\mathit{U},\Lambda )}{\partial {\rho }_e}}={\displaystyle {\int }_{\mathrm{0}}^{t_{\mathrm{f}}}{\Lambda }^{\mathrm{T}}}(t_{\mathrm{f}}-t)\left({\displaystyle \frac{\partial \mathit{M}(\mathit{\rho })}{\partial {\rho }_e}}\ddot{\mathit{U}}+{\displaystyle \frac{\partial C_{\mathrm{fb}}(\mathit{\rho })}{\partial {\rho }_e}}\dot{\mathit{U}}+{\displaystyle \frac{\partial \mathit{K}(\mathit{\rho })}{\partial {\rho }_e}}\mathit{U}-{\displaystyle \frac{\partial \mathit{F}(t)}{\partial {\rho }_e}}\right)\mathrm{d}t.$

4 Numerical implementation

The TO problem implementation given the input data that includes the nodal Cartesian coordinates, their connectivity, boundary conditions and initial material distribution, begins with the dynamic FE method analysis for which the objective function is calculated. If the convergence is not achieved, the sensitivities are calculated and the optimizatioin problem is solved through the LP algorithm. The physical and adjoint dynamic systems, Eqs. (26) and (45), respectively, are solved by the α-Newmark time-integration method. The design variables are updated and the next iteration loop initiates, as indicated in the flowchart of Fig. 3.

Fig.3 Flowchart for the implemented optimization procedure

Full size|PPT slide

In this work, the approach proposed by Ref. [19] is used to calculate the pseudo-densities of each FE by applying a linear projection filter to the design variables, which are defined for the centroid coordinate of the element. This linear projection technique solves the undesired mesh dependency problem that arises when intermediate material densities are penalized.

Regarding the Newmark’s time-integration, the sampling rate fs and the time-step $\mathrm{\Delta }t$ implemented depend on the first resonance frequency ${\omega }_{\mathrm{1}}$ (in Hz) of the topology at each optimization iteration, as follows:

$fs=\mathrm{20}{\omega }_{n_{\mathrm{1}}},$

(47)$\mathrm{\Delta }t={\displaystyle \frac{\mathrm{1}}{fs}}.$

The optimization solver based in LP is implemented with iterative adjustments of the design variables within a percentage of their original values, what is called the moving limits and which are dependent on the objective function value history:

(48)${\rho }_l=\{\begin{array}{l}\mathrm{0.85}{\rho }_{\mathrm{it}},f({\rho }_{\mathrm{it}})-f({\rho }_{\mathrm{it}-1})\ge \mathrm{0.6}(f({\rho }_{\mathrm{it}-1})-f({\rho }_{\mathrm{it}-\mathrm{2}})),\hfill \\ \mathrm{0.90}{\rho }_{\mathrm{it}},\mathrm{0.3}(f({\rho }_{\mathrm{it}-1})-f({\rho }_{\mathrm{it}-\mathrm{2}}))\le f({\rho }_{\mathrm{it}})-f({\rho }_{\mathrm{it}-\mathrm{1}})\le \mathrm{0.6}(f({\rho }_{\mathrm{it}-\mathrm{1}})-f({\rho }_{\mathrm{it}-\mathrm{2}})),\hfill \\ \mathrm{0.95}{\rho }_{\mathrm{it}},f({\rho }_{\mathrm{it}})-f({\rho }_{\mathrm{it}-1})\le \mathrm{0.3}(f({\rho }_{\mathrm{it}-\mathrm{1}})-f({\rho }_{\mathrm{it}-\mathrm{2}})),\hfill \end{array}$

(49)${\rho }_u=\{\begin{array}{l}\mathrm{1.15}{\rho }_{\mathrm{it}},f({\rho }_{\mathrm{it}})-f({\rho }_{\mathrm{it}-\mathrm{1}})\ge \mathrm{0.6}(f({\rho }_{\mathrm{it}-\mathrm{1}})-f({\rho }_{\mathrm{it}-\mathrm{2}})),\hfill \\ \mathrm{1.10}{\rho }_{\mathrm{it}},\mathrm{0.3}(f({\rho }_{\mathrm{it}-\mathrm{1}})-f({\rho }_{\mathrm{it}-\mathrm{2}}))\le f({\rho }_{\mathrm{it}})-f({\rho }_{\mathrm{it}-\mathrm{1}})\le \mathrm{0.6}(f({\rho }_{\mathrm{it}-\mathrm{1}})-f({\rho }_{\mathrm{it}-\mathrm{2}})),\hfill \\ \mathrm{1.05}{\rho }_{\mathrm{it}},f({\rho }_{\mathrm{it}})-f({\rho }_{\mathrm{it}-\mathrm{1}})\le \mathrm{0.3}(f({\rho }_{\mathrm{it}-\mathrm{1}})-f({\rho }_{\mathrm{it}-\mathrm{2}})),\hfill \end{array}$

where ${\rho }_l$ and ${\rho }_u$ are, respectively, the lower and upper bounds of the design variables, and “it” is the current iteration step. After each optimization step, the design variables vector $\mathit{\rho }$is obtained and updated in the design domain.

5 Numerical results

In this section, a numerical example is presented to demonstrate the methodology described in this text. The optimized topology of the host structure is obtained for two distinct metallic materials for the predefined locations of the piezoelectric sensor and actuator layers. The domain characteristics are given in Table 1. The material properties of PZT-5A [31] are: $c_{11}$=139 GPa, $c_{12}$=77.84 GPa, $c_{13}$=74.28 GPa, $c_{33}$=115.41 GPa, $c_{44}$=25.64 GPa, $c_{66}$=30.58 GPa, $e_{31}$=‒14.8 C/m², $e_{33}$=15.7835 C/m², $e_{15}$=12.67 C/m², ${\epsilon }_{11}$=1474.36 pF/m, ${\epsilon }_{33}$=1299.38 pF/m, ${\epsilon }_0$=8.854188 pF/m, $\zeta$=7500 kg/m³. The material properties of mat1 and mat2 in Eqs. (36) and (37) are shown in Table 2.

Tab.1 Design domain definitions for the vibration attenuation problem

Parameter	Symbol	Value	Unit
Domain size (width, heigh, depth)	$(w_{\mathrm{d}},h_{\mathrm{d}},d_{\mathrm{d}})$	(20,4,0.2)	cm
Substrate size (width, heigh, depth)	$(w_{\mathrm{s}},h_{\mathrm{s}},d_{\mathrm{s}})$	(20,2.67,0.2)	cm
Piezoceramic size (width, heigh, depth)	$(w_{\mathrm{p}},h_{\mathrm{p}},d_{\mathrm{p}})$	(20,0.67,0.2)	cm
Substrate mesh size	($x_{\mathrm{s}}$,$y_{\mathrm{s}}$)	(120,16)	‒
Piezo ceramic mesh size	($x_{\mathrm{p}}$,$y_{\mathrm{p}}$)	(120,4)	‒
Mechanical load amplitude	$\Vert \mathit{F}\Vert$	1	N
Initial material volume	$V_{\mathrm{ini}}$	30%	‒
Material volume constraint	$V_{\max }$	70%	‒
Proportional damping	(α,β)	(0,1×10‒5)	‒

Tab.2 Material properties (mat1 and mat2 in Eqs. (36) and (37))

Material	$E$/GPa	$\zeta$/(kg∙m−3)	$\nu$
Fictitious (mat1)	300	1	0.24
Alumina (mat2)	392	4000	0.24

The proposed two-dimensional BPEA has its metallic nodes clamped in x and y directions at one end and has a transient load applied at its free opposite end (Fig. 4). The AVFC aims to minimize the vibration response measured at the loaded tip degree-of-freedom. The fixed design domain is discretized into four-noded bi-linear FEs with three degrees-of-freedom per node. Implementation has been performed with MATLAB® software using its embedded LP optimizer.

The penalization coefficients ($p_{\mathrm{1}}$ and $p_{\mathrm{2}}$) update is given by the curves shown in Fig. 5. The radius value of the projection filter applied is such that the circumference encloses five elements diagonally.

Fig.4 Boundaries and transient load

Full size|PPT slide

Fig.5 Penalization coefficients update

Full size|PPT slide

After 35 iteration steps the convergence is achieved, resulting in topologies with a shadowy gray scale region for the feedback gain values K (K₁=0, K₂=9×10⁴, K₃=3.6×10⁵, K₄=7.2×10⁵). The resulting designs are shown in Fig. 6. The pseudo-densities $\rho$ are plotted for the obtained topologies in Figs. 6(a), 6(c), 6(e) and 6(g). In Figs. 6(b), 6(d), 6(f) and 6(h), those resulted material distribution are interpreted as material Type 1 (black) for design variables greater than or equal to 0.7, as material Type 2 (white) for design variables smaller than 0.7, and the ceramics are the gray regions. The influence of the different constant gains to the final topology is investigated.

The respective convergence curves during the optimization steps are shown in Fig. 7, and the volume constraint curves are shown in Fig. 8.

Fig.6 Optimized topologies for different feedback gain values. (a) Optimized topology for $K=\mathrm{0}$; (b) post-processed topology for $K=\mathrm{0}$; (c) optimized topology for K=9×104; (d) Post-processed Topology for K=9×104; (e) optimized topology for K=3.6×105; (f) post-processed topology for K=3.6×105; (g) optimized topology for K=7.2×105; (h) post-processed topology for K=7.2×105

Full size|PPT slide

For a homogeneous initial material distribution of 30%-material Type 1, the optimization managed to minimize the vibration energy function by increasing the concentration of the less stiff material (mat1). The volume constraint of 70%-material Type 1 is inactive for all simulated feedback gains, as shown in Fig. 8. Regarding the influence of the control gain, Table 4 shows that as K-value increases, the concentration of material Type 2 increases, which has greater specific weight and consequently the first natural frequency decreases. As a natural frequency constraint is not specified in the optimization problem formulation, the increase in the feedback gain compensates for the lower energy dissipation within the specified time interval caused by a lower natural frequency. This described tendency can be seen for all of the feedback gains except for K=0, where the absence of additional damping effect forces the proportional structural damping to be larger by increasing the concentration of mat2, comparing to K=9×10⁴, and K=3.6×10⁵.

Fig.7 Convergence comparison evaluated by Eq. (33) for $K=K_i\times {\mathrm{10}}^{\mathrm{5}}$ (i=1,2,3,4)

Full size|PPT slide

Fig.8 Volume constraint comparison for $K=K_i\times {\mathrm{10}}^{\mathrm{5}}$ (i=1,2,3,4)

Full size|PPT slide

The topologies obtained for K=0 and K=3.6×10⁵ have both approximately the same material volume but the first one shows a symmetric material distribution, what contrasts with the antisymmetric material distribution shown in the second case, Figs. 6(b) and 6(f), respectively.

Therefore, as seen in Fig. 7, as the constant gain value increases, the objective function converges at a lower value, what is quantitative shown in Table 4. The transient responses in Fig. 9 shows that the faster the oscillation decays, the higher the transient voltage peak (Fig. 10) inputted into the actuator, what demostrates the AVFC influence.

Tab.3 Objective function values for each optimized topology

K		Optimized			Post-processed
	$f(\mathit{\rho })$		$V_{\mathrm{mat1}}$	${\omega }_{n_{\mathrm{1}}}$/Hz	$f(\mathit{\rho })$	$V_{\mathrm{mat1}}$	${\omega }_{n_{\mathrm{1}}}$/Hz
0.00	3.6508×10−15		60.02%	1217.7312	3.6638×10−15	53.96%	1220.31
9.00×104	3.4135×10−15		60.95%	1217.7312	3.4270×10−15	55.05%	1222.25
3.60×105	3.0654×10−15		60.17%	1214.3580	3.0773×10−15	54.01%	1218.80
7.20×105	2.8837×10−15		58.15%	1214.3580	2.8947×10−15	51.67%	1210.43

Fig.9 Transient response for the optimized topologies $K=K_i\times {\mathrm{10}}^{\mathrm{5}}$ (i=1,2,3,4)

Full size|PPT slide

Fig.10 Transient feedback voltage for the optimized topologies $K=K_i\times {\mathrm{10}}^{\mathrm{5}}$ (i=1,2,3,4)

Full size|PPT slide

When the topologies of Fig. 6 are post-processed for the feedback gains (K₁=0, K₂=9×10⁴, K₃=3.6×10⁵, K₄=7.2×10⁵), the objective function values shown in Fig. 11 are obtained. The effect of designing the BPEA structure by applying the TOM with a feedback gain can be seen by comparing the objective function values of a topology obtained with a higher constant gain with the topology obtained with a lower constant gain but post-processed with the higher gain. The last statement can be exemplified by taking the objective function value of the K₄-topology post-processed with K₄-gain. As it should be, its value is smaller than the objective function value of K₃-topology post-processed with K₄-gain.

Fig.11 Analysis of the objective function values for the optimized topologies subjected to the proposed gains $K=K_i\times {\mathrm{10}}^{\mathrm{5}}$ (i=1,2,3,4)

Full size|PPT slide

6 Conclusions

The TO applied to the design of the host structure of BPEAs with active vibration suppression is carried out in this paper. The extra terms added to the damping matrix given by the piezoelectric sensor layer readings in time domain, has been shown to interfere to the metallic material distribution within the design domain such that a slower response decay, given by a lower natural frequency of the structure, is compensated by the electrical damping terms, what attenuates the vibration response of the piezoelectric system. It has been shown that as the constant gain increases, the objective function achieves a lower minimum and the concentration between the metallic materials are rebalanced. The final topology is though dependent on the feedback control gain, what shows the importance and viability of solving the TO of a controlled piezoelectric system in time domain, with no loss to the transient response information. Besides the fact that the material distribution within the domain does not vary significantly as the feedback gain increases, the loss of symmetry is evident when the constant gain goes from zero to the maximum value considered. Other combination of material types, with lower elastic modulus for instance, might be considered so the stiffness of the structure depends more on the material distribution when the feedback gain varies than on the presence of the piezoceramics.

References

[1] Ou J S, Kikuchi N. Integrated optimal structural and vibration control design. Structural Optimization, 1996, 12(4): 209–216

[2] Gibson R F. A review of recent research on mechanics of multifunctional composite materials and structures. Composite Structures, 2010, 92(12): 2793–2810

[3] Yang X, Li Y. Topology optimization to minimize the dynamic compliance of a bi-material plate in a thermal environment. Structural and Multidisciplinary Optimization, 2013, 47(3): 399–408

[4] Zhang X, Kang Z. Topology optimization of piezoelectric layers in plates with active vibration control. Journal of Intelligent Material Systems and Structures, 2014, 25(6): 697–712

[5] Gonçalves J F, De Leon D M, Perondi E A. Topology optimization of embedded piezoelectric actuators considering control spillover effects. Journal of Sound and Vibration, 2017, 388: 20–41

[6] Zhang X, Kang Z, Li M. Topology optimization of electrode coverage of piezoelectric thin-walled structures with CGVF control for minimizing sound radiation. Structural and Multidisciplinary Optimization, 2014, 50(5): 799–814

[7] Zhao J, Wang C. Dynamic response topology optimization in the time domain using model reduction method. Structural and Multidisciplinary Optimization, 2016, 53(1): 101–114

[8] Zhang X, Kang Z. Dynamic topology optimization of piezoelectric structures with active control for reducing transient response. Computer Methods in Applied Mechanics and Engineering, 2014, 281: 200–219

[9] Takezawa A, Makihara K, Kogiso N, et al. Layout optimization methodology of piezoelectric transducers in energy-recycling semi-active vibration control systems. Journal of Sound and Vibration, 2014, 333(2): 327–344

[10] Tzou H, Tseng C. Distributed piezoelectric sensor/actuator design for dynamic measurement/control of distributed parameter systems: A piezoelectric finite element approach. Journal of Sound and Vibration, 1990, 138(1): 17–34

[11] Wang S, Quek S, Ang K. Vibration control of smart piezoelectric composite plates. Smart Materials and Structures, 2001, 10(4): 637–644

[12] Caruso G, Galeani S, Menini L. Active vibration control of an elastic plate using multiple piezoelectric sensors and actuators. Simulation Modelling Practice and Theory, 2003, 11(5‒6): 403–419

[13] Ray M, Reddy J. Optimal control of thin circular cylindrical laminated composite shells using active constrained layer damping treatment. Smart Materials and Structures, 2004, 13(1): 64–72

[14] Vasques C M A, Dias Rodrigues J. Active vibration control of smart piezoelectric beams: Comparison of classical and optimal feedback control strategies. Computers & Structures, 2006, 84(22‒23): 1402–1414

[15] Kim S M, Wang S, Brennan M J. Comparison of negative and positive position feedback control of a exible structure. Smart Materials and Structures, 2011, 20(1): 015011

[16] Choi S B, Kim H S, Park J S. Multi-mode vibration reduction of a CD-ROM drive base using a piezoelectric shunt circuit. Journal of Sound and Vibration, 2007, 300(1–2): 160–175

[17] Chan K W, Liao W H. Precision positioning of hard disk drives using piezoelectric actuators with passive damping. In: Proceedings of 2006 IEEE International Conference on Mechatronics and Automation. Luoyang: IEEE, 2006, 1269–1274

[18] Zhang S, Schmidt R, Qin X. Active vibration control of piezoelectric bonded smart structures using PID algorithm. Chinese Journal of Aeronautics, 2015, 28(1): 305–313

[19] Bruns T E, Tortorelli D A. Topology optimization of non-linear elastic structures and compliant mechanisms. Computer Methods in Applied Mechanics and Engineering, 2001, 190(26–27): 3443–3459

[20] Jensen J S, Nakshatrala P B, Tortorelli D A. On the consistency of adjoint sensitivity analysis for structural optimization of linear dynamic problems. Structural and Multidisciplinary Optimization, 2014, 49(5): 831–837

[21] Moetakef M A, Lawrence K L, Joshi S P, Closed-form expressions for higher order electroelastic tetrahedral elements. AIAA Journal, 1995, 33(1): 136–142

[22] Reddy J N. An Introduction to the Finite Element Method. Vol. 2. New York: McGraw-Hill, 1993

[23] Dong X, Ye L, Peng Z, A study on controller structure interaction of piezoelectric smart structures based on finite element method. Journal of Intelligent Material Systems and Structures, 2014, 25(12): 1401–1413

[24] Tcherniak D. Topology optimization of resonating structures using simp method. International Journal for Numerical Methods in Engineering, 2002, 54(11): 1605–1622

[25] Zhao J, Wang C. Topology optimization for minimizing the maximum dynamic response in the time domain using aggregation functional method. Computers & Structures, 2017, 190: 41–60

[26] Rozvany G. Aims, scope, methods, history and unified terminology of computer-aided topology optimization in structural mechanics. Structural and Multidisciplinary Optimization, 2001, 21(2): 90–108

[27] Bendsøe M P, Sigmund O. Topology optimization: Theory, Methods and Applications. Berlin: Springer, 2003

[28] Zargham S, Ward T A, Ramli R, Topology optimization: A review for structural designs under vibration problems. Structural and Multidisciplinary Optimization, 2016, 53(6): 1157–1177

[29] Rozvany G I. A critical review of established methods of structural topology optimization. Structural and Multidisciplinary Optimization, 2009, 37(3): 217–237

[30] Yin L, Ananthasuresh G. Topology optimization of compliant mechanisms with multiple materials using a peak function material interpolation scheme. Structural and Multidisciplinary Optimization, 2001, 23(1): 49–62

[31] Balamurugan V, Narayanan S. A piezolaminated composite degenerated shell finite element for active control of structures with distributed piezosensors and actuators. Smart Materials and Structures, 2008, 17(3): 035031

Acknowledgements

The first author thanks CAPES (Coordination for Improvement of Higher Education Personnel) for the doctoral fellowship. The second author is pleased to acknowledge the support by CNPq (National Council for Scientific and Technological Development) under grant 302658/2018-1 and CAPES project number A023-2013.

Open Access

This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

RIGHTS & PERMISSIONS

2019 The Author(s) 2019. This article is published with open access at link.springer.com and journal.hep.com.cn