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

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.


Introduction
In this work, the topology optimization (TO) extended to a bi-material metallic design of piezoactuators in timedomain 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 nonstructural 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 finite element (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][11][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 constantgain 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 optimized 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 Q s p , that is differentiated in time, _ Q s p is amplified by the control gain K, and the resulting voltage Φ a p is fed back into the actuator piezoelectric layer attached to the opposite surface of the host structure, as illustrated in Fig. 1.
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.

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: where T and S are the stress and strain tensors, respectively, and D is the electric displacement vector. The piezoelectric elastic matrix c E , evaluated under a The system dynamics is approximated by a 2D-solid bilinear FE where the displacement field u and the electric potential f are interpolated by means of the shape functions N u and N f , respectively. The displacement subdivided into free voltages, Φ f , and prescribed voltages, Φ p . Therefore, the dynamic FE equilibrium equations for the piezoelectric system are written [21]: where the boundary surface is composed by the prescribed traction region, Γ T , and the electroded region, Γ D . The subscript uu denotes structural states, the subscript f denotes electrical states, and B u and B f are, respectively, the mechanical and electrical shape functions derivatives [22]. The proportional damping matrix C uu is defined by where α and β are Rayleigh's coefficients, that are nonnegative and small in the sense that the effective damping coefficients is < < 1 [23], and are considered to be designindependent [24,25].

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 Q s p , by imposing null voltage on the sensor electrodes Φ s p ¼ 0 when the system is under a mechanical transient load. Knowing that the electric charges at the piezoelectric internal nodes Q f are null and that the actuator and sensor electrodes do not share degrees of freedom, The derivative of the electrical charge Eq. (12) results in the nodal electric current vector I s on the external sensor electrode: that can be added by pre-multiplying I s by the unity vector I f s p ¼ 1 ⋯ 1 f g , and the measured current in the top electrode of the sensor layer is obtained: To establish a velocity feedback control law, the actuator input voltage Φ a p is defined as a current amplifier: where I f a p ¼ 1 ⋯ 1 f g 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): Now, substituting Eq. (17) into Eqs. (10), (11) and (13), we have and the coupled system to be solved is given by Eqs. (18) and (19): where C uu is given by Eq. (9), 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: for The modified damping matrix C 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.

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 an energy functional involving the displacement of a predefined degree-offreedom as its squared integral spanned over a time interval ½0,t f [8], where the ρ-vector elements are the design variables, and As for the initial conditions, this dynamic problem depends on the magnitude of the input load. For a pointwise mechanical load with the profile shown in Fig. 2, the initial displacement vector U 0 is given by the static analysis KU 0 ¼ Fð0Þ while the initial velocity vector V 0 is null. A volume constraint V max limits the distribution of one of the materials. The optimization problem is therefore stated as bellow: s:t: 10 -15 £ n £1: As stated in optimizatioin problem Eq. (35), its constraints also includes the dynamic equilibrium Eq. (26), where FðtÞ ¼ FyðtÞ, and boundary values for each element design variable n of a total of N elements. The metallic materials distribution is such that it minimizes the cost function Eq. (33).

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][28][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, c e ðÞ, and the specific weight, ζ , specified for the elastic domain are functions of the pseudo-density defined for each FE e: where ðx c ,y c Þ is the centroid cartesian coordinate pair of each element, and c mat1 and c mat2 are the elastic tensors for the first and second materials, respectively.

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 ρ. These derivatives, or the sensitivity analysis is given by the adjoint method which involves an adjoint vector lðtÞ that is chosen as the solution to the adjoint equations derived from a Lagrangian function LðU,lÞ [20]: Taking the derivative of the Lagrangian Eq. (38) with respect to the design variables ρ we get the expression that is solved by integration by parts: MðρÞl T ðtÞ t f 0 : To eliminate the state variables derivatives ∂Uðt f Þ ∂ e and ∂ _ U ðt f Þ ∂ e from Eq. (43), we define the adjoint variable l 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 C fb _ l. Given that BUðtÞ þ U T ðtÞB ¼ 2u dof ðtÞ, the adjoint system that needs to be solved is In order to redefine the system Eq. (44) as an initial value primal problem, we apply the change of variables τðtÞ ¼ t ft, and the adjoint system to be solved is as follows: where ΛðτðtÞÞ ¼ lðtÞ and _ ΛðτðtÞÞ ¼ -_ lðtÞ. Therefore, the sensitivity expression for the transient problem reduces to ∂LðU,ΛÞ ∂ e

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. 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 Δt implemented depend on the first resonance frequency ω 1 (in Hz) of the topology at each optimization iteration, as follows: 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: where l and 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 ρ is obtained and updated in the design domain.

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 Table 2.
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 bilinear FEs with three degrees-of-freedom per node. Implementation has been performed with MATLAB® software using its embedded LP optimizer.
The penalization coefficients (p 1 and p 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.
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 3 shows that as K-value increases, the   (f) post-processed topology for K = 3.6Â10 5 ; (g) optimized topology for K = 7.2Â10 5 ; (h) post-processed topology for K = 7.2Â10 5 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 4 , and K = 3.6Â10 5 .
The topologies obtained for K = 0 and K = 3.6Â10 5 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 3. 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.
When the topologies of Fig. 6 are post-processed for the feedback gains (K 1 = 0, K 2 = 9Â10 4 , K 3 = 3.6Â10 5 , K 4 = 7.2Â10 5 ), 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 4 -topology post-processed with K 4 -gain. As it should be, its value is smaller than the objective function value of K 3 -topology post-processed with K 4 -gain.  Table 3 Objective function values for each optimized topology K

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.