Design and testing for shape control of piezoelectric structures using topology optimization

Engineering Structures 97 (2015) 90–104

Contents lists available at ScienceDirect

Engineering Structures journal homepage: www.elsevier.com/locate/engstruct

Design and testing for shape control of piezoelectric structures using topology optimization Quantian Luo, Liyong Tong ⇑ School of Aerospace, Mechanical and Mechatronic Engineering, The University of Sydney, Sydney, NSW 2006, Australia

a r t i c l e

i n f o

Article history: Received 11 July 2014 Revised 1 April 2015 Accepted 1 April 2015

Keywords: Shape control Optimization Mutual strain energy Moving iso-surface threshold

a b s t r a c t This paper presents a numerical and experimental investigation into optimum topological design of morphing piezoelectric structures using a moving iso-surface threshold method. Proposed first is a novel formulation of the response function suitable for minimizing the error norm between desired and achieved shapes in the form of combined mutual strain energy densities. A design updating algorithm with an enhanced efficiency is then developed in searching for the optimum iso-surface threshold during iterations. Numerical results are presented for optimum topological material distribution for shape morphing of piezoelectric plates subjected to mechanical loading, one channel electrical voltage and a combination of both mechanical and electrical loadings. Experiments for selected optimal topological designs are conducted to validate the present response function and updating algorithm as well as the numerical results. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction Piezoelectric materials and structures have been widely studied due to the coupling effects between electrical and mechanical fields. Flexible structures with piezoelectric actuators and sensors have found a wide range of applications including structural vibration control, structural health monitoring and structural shape control. Shape control or morphing represents an important application [1–4] and there is a growing interest in morphing flexible structures using piezoelectric actuators and sensors. A desired structural shape of piezoelectric (PZT) structures may be achieved by optimal input voltages [1,3,5], optimal designs of piezoelectric actuators [6–8] and material distributions [9,10]. Finite element based topology optimization is useful for shape morphing to acquire better structural performance as it can be used to optimize voltage input, piezoelectric material distribution and structural material topology [4,11–13]. In this paper, optimal designs of structural materials for morphing piezoelectric structures will be investigated by using a moving iso-surface threshold method (MIST). Finite element based topology optimization has been substantially studied and many methods have been proposed. A homogenization method [14] is to obtain optimal structural design by homogenizing anisotropic composite with micro-scale voids. In ⇑ Corresponding author. Tel.: +61 2 93516949; fax: +61 2 93514841. E-mail address: [email protected] (L. Tong). http://dx.doi.org/10.1016/j.engstruct.2015.04.006 0141-0296/Ó 2015 Elsevier Ltd. All rights reserved.

the solid isotropic material with penalization method (SIMP) [15–17], only one variable (0 or 1 in a white-black design) is required to describe one element. In the evolutionary structural optimization method (ESO) [18,19], elements with lower physical responses are gradually removed. In the level set method (LSM) [20–22], design boundary is implicitly expressed using the zero level set of a higher-order surface function. In the moving isosurface threshold method (MIST) [23,24], a physical response for an objective function is used; each element is depicted using one variable and a moving iso-surface threshold is used to evolve design boundary. That is, features of ESO, SIMP and LSM methods are integrated into MIST. This method has been used to maximize structural stiffness and to design gripping mechanisms and pressurized actuators [23,24]. It will be further studied and used to investigate shape morphing for piezoelectric structures in this paper. A displacement component at one point can be obtained by applying a unit dummy force based on the principle of virtual work. By using this principle, optimization to design compliant mechanism has been studied by many authors [23,25,26]. In the mechanism design, an objective function can be defined as a difference between the target and realized displacements at one point where a flexible spring is often used. In the present investigation, structural displacements are firstly expressed by using mutual strain energy and the plate shapes are represented by using deflections at observation points. An objective function is defined by a sum of squared errors between the

91

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

desired and computed deflections at these points. A novel formulation of a response function is derived on the basis of mutual strain energy density and is then used to find optimal topologies of structural and piezoelectric materials as well as voltages for morphing PZT structures. A novel algorithm is also developed to efficiently find an iso-surface level in an iterative process. Numerical results computed by using MIST are presented for morphing bending and twisting shapes of plates subjected to mechanical and/or electrical loadings. A combination of the two typical shapes can comprise different shapes and thus they are used to show effectiveness of the present algorithm. Specimens for the obtained optimal designs are fabricated and experimental tests are conducted to verify the present computations. Overshoot issues caused by applied voltages in piezoelectric structures are experimentally investigated, which are important to achieve precise shape control in practice.

2. Problem statement and formulation 2.1. Basic problem statement Topology optimization for shape morphing of piezoelectric structures can be described as [8,12,27] to find the optimal material distribution that minimizes the selected shape error norm Er subject to relevant equilibrium equations for PZT structures and material volume constraint. In addition, constraints on electric energy consumption or electric field strength should also be introduced [8,12,28]. The shape morphing problem is to minimize an objective function defined as a sum of squared errors between the desired and achieved shapes [1,3,12]:

Minimize :

n X Er ¼ ðui  udi Þ2

ð1aÞ

½K uu fU k g ¼ fF k g  ½K uu fVg

ð1bÞ

½K uu fU i g ¼ fF i g

ð1cÞ

Ne X ðxe v e Þ 6 v f v 0

ð1dÞ

e¼1

maxfV e g 6 V max

ð1eÞ

0 < xmin 6 xe 6 1

ð1fÞ

where ui and udi denote the achieved (or computed) and desired (or specified) displacements at point i; ½K uu  and ½K uu ð¼ ½K uu T Þ are the structural stiffness and piezoelectric coupling matrices; fF k g and fF i g denote the applied and unit dummy force vectors; fU k g and fU i g are the displacement vectors due to the applied loads and the unit dummy force; {V} is the voltage vector; Ve (e = 1, 2, . . . , Ne) is the voltage of the eth PZT element and Ne is the total element number; Vmax denotes the maximum voltage; ve and v0 denote the volume of an element and entire design domain; vf denote the prescribed volume fraction; {xe} is the design variable vector and xe ranges from 0 to 1; xmin (= 103 normally) is the minimum value of the design variables. In a black–white design, xe = 1 or ?0. Eq. (1b) is a simplified equation for piezoelectric structures with converse piezoelectric effect only. When direct and converse piezoelectric effects are considered, equilibrium equations in finite element analysis can be expressed as [1]:

"

K uu K Tuu

# K uu  U   F  ¼ K uu V Q

1

fUg ¼ ½K uu  fFg

ð2bÞ

where

½K uu  ¼ ½K uu   ½K uu ½K uu 1 ½K Tuu 

ð2cÞ

½F ¼ fFg  ½K uu ½K uu 1 fQ g

ð2dÞ

If the electro-mechanical fields are weakly coupled, which is true for many commercially-available piezoelectric materials, the actuation equation can be simplified as [1]:

  fUg ¼ ½K uu 1 fFg  ½K uu fVg

ð2eÞ

In shape morphing for piezoelectric structures, this equation is normally adopted [1,2,29] and will also be used in this paper, When Eq. (1) is used in topology optimization, the structural stiffness and piezoelectric coupling matrices in an iterative process are calculated by [8,12]:

½K uu  ¼

Ne  X

str

pzt

xpe ½ke  þ xqe ½ke 



ð3aÞ

e¼1

½K uu  ¼

Ne X xre ½kðuuÞe 

ð3bÞ

e¼1 str

i¼1

Subject to :

where [Kuu] is the dielectric matrix and {Q} denotes the electric charges. In finite element analysis for piezoelectric structures, the actuation equation can be obtained by [1]:

ð2aÞ

pzt

in which [k(uu)e] is the piezoelectric coupling matrix; ½ke  and ½ke  are the stiffness matrices for host structural and piezoelectric materials with solid materials; subscript e represents the eth element; p, q and r are the penalty parameters. p is the penalty of structural materials; q and r are the penalty parameters for the elastic and piezoelectric constants of piezoelectric materials. In a twodimensional finite element analysis, penalties p and q can be chosen as [8,12]: p = 3 or p = 0; q = 3 or q = 0; r = 1. For PZT laminated plates, if p = 3 and q = 3, topology of the PZT material is the same as that of the host structural material and the host structural materials or the PZT materials will not be removed when p = 0 or q = 0. Morphing PZT structures can be achieved by optimal voltage channel designs [1,3,28] with an energy consumption and/or maximum voltage constraint. As indicated in [12], an arbitrary spatially distributed voltage field for the PZT actuators may be hardly realized as an optimal voltage distribution usually demands a complicated electric channel configuration so that it is too expensive to implement. Therefore, one channel actuation voltage was used in [12]. In this case, the maximum voltage constraint in Eq. (1e) is used. In this paper, one channel voltage input is also considered; topology optimization for morphing piezoelectric structures will be investigated by using MIST where an objective function is derived on the basis of mutual energy density and the numerical results with experimental validation are presented. 2.2. MIST formulation MIST [23,24] comprises of three key points: (1) expressing both objective function and constraint in the form of an integral over the design domain, (2) selecting the response (U) function, and (3) searching for optimal solutions for a sequential approximate optimization formulations based on structural responses of previous approximations.

92

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

2.2.1. Objective function It is well-known that displacement ui at point i can be expressed in terms of the mutual potential energy [23,30,31] calculated by applying a unit dummy force at point i along the direction of displacement ui. Therefore, by the unit load method, the objective function in Eq. (1a) can be expressed as:

Subject to :

  ½K uu k xpe ; k xqe fk U k g ¼ fF k g  ½K uu ðk xre ÞfVg

ð7bÞ

 K uu ðk xpe ; k xqe Þ fk U i g ¼ fF i g

ð7cÞ

Z

Hðt; k UÞdX 6 v f v 0

ð7dÞ

X

Er ¼

n X

n X

i¼1

i¼1

 i ui ¼ w

 i ðEmus Þi w

ð4aÞ

in which

( i ¼ w

where

2wi ðui  2udi þ u2di =ui Þ when ui – 0 0

when ui ¼ 0

ði ¼ 1; 2; . . . ; nÞ ð4bÞ

where wi (= 1, 2, . . . , n) are the weighted factors of designer’s choice; and ðEmus Þi denotes the mutual strain energy associated with ui calculated by [23,30,31]:

ui ¼ 2ðEmus Þi ¼

Z 

Z    frr gT fei g dX ¼ fri gT fer g dX ði ¼ 1;2;...; nÞ

X

where X denotes the design domain; frr g and fer g are the stress and strain fields created by the real applied loading; and fri g and fei g are the stress and strain fields produced by the applied unit dummy force at point i along ui direction. By substituting Eq. (4c) into Eq. (B1), we can re-define an objective function in the following integral format:

Z X

! ! Z n n 1X 1X  i frr gT fei gi dX ¼  i fri gT fer gi dX w w 2 i¼1 X 2 i¼1

ð5Þ

Structural displacement vector {Uk} caused by the applied forces and voltages is found by using Eq. (1b) and virtual displacement vector fU i g is calculated by using (1c). Displacements ui (i = 1, 2, . . . , n) are computed by using Eq. (4c) or mutual strain energy. The desired structural shape may be depicted by using displacements udi (i = 1, 2, . . . , n) of some points, e.g., a shape of plate deformations can be described by using deflections at some typical points. Therefore, Er defined in Eq. (5) can be used as an objective function for shape morphing and the optimization is to minimize Er. 2.2.2. Selection of the U function In MIST [23,24], an appropriate selection of the response function (U) plays a critical role in updating the design variable xe (e = 1, 2, . . . , Ne). If we consider the material at a point as solid / P t and void otherwise (or vice versa) and define the design variable xe as a ratio of the volume of the solid component to the total one of element e, all the design variables can be determined by the chosen response function (U) and the iso-surface threshold t. To minimize the squared error sum or an objective function given in Eq. (5), the U function can be constructed by using mutual strain energy density Emusd:

U¼

k

U¼

n k   1X   1X T T k  i fk ri g fk er g wi fk rr g fk ei g ¼ w 2 i¼1 2 i¼1n

ð7eÞ

ð7fÞ

The superscript k on the left top corner of the mathematical symbols represents the quantities of the kth iteration. By comparing Eqs. (7a)–(7f) with (1a)–(1f), it is seen that the key points in MIST are to find the response (U) function and an iso-surface threshold value to update design variables. 3. MIST algorithm

X

ð4cÞ

Er ¼

0 < xmin 6 xe 6 1

n n   1X   1X  i frr gT fei g ¼  i fri gT fer g w w 2 i¼1 2 i¼1

ð6Þ

In MIST, an iso-value surface of the response function U is sought to define the structural boundary by meeting the volume constraint and also to minimize an objective function [23,24]. The MIST algorithm based on standard finite element analysis without accessing the source code can be described as follows. 3.1. Step 1: Initialization The initialization mainly includes two sets of input data. One contains all necessary data for finite element analysis that are not changed in an iterative process and the other requires selection of initial values of the design variables or weighting factors {xe}0 with reference to the volume constraint that will be updated in each iteration. The weighting factors are similar to element based material density variables used in SIMP. A uniform material distribution is often assumed as an initial condition. Parameters such as a volume constraint, penalty, move limit and a filter radius are also specified in the initialization. 3.2. Step 2: Finite element analysis using current fxe gk As discussed for Eq. (2), the simplified equilibrium equation of finite element analysis (FEA) is often used in optimization for piezoelectric structures and it will also be used in this paper. In present FEA computations, the material properties for each element are calculated based on either the initial or the updated design variables or weighting factors {xe}k by using an appropriate material model in SIMP while all other FEA input data remain unchanged. An input data file for FEA is then created in the 1st iteration or modified in subsequent iterations by using the updated weighting factors and the other data required for FEA. By using the initialized or updated input data file, FEA can be conducted by solving Eqs. (7b) and (7c) and the output data file containing displacements, strains and stresses are generated. In the present computation, an in-house FEA program [3] for piezoelectric structures is used. 3.3. Step 3: Construction of the U function

2.2.3. MIST problem formulation By using the selected U function in Eq. (6), the problem given in Eq. (1) can be re-formulated in a MIST format [23] as follows: Find {xe} governed by the selected U function and the iso-surface threshold so as to:

Minimize :

Er ¼

Z X

k

UHðt; k UÞdX

ð7aÞ

3.3.1. Nodal value calculation The surface of the chosen U function for two-dimensional analysis can be constructed by the calculated U values using Eq. (6) (as z coordinate) at each node (located by its x and y coordinates in Cartesian coordinate system) by using the FEA results of the initial or previous iteration. In the FEA output data, physical responses, such as stress, strain and strain energy density, may

93

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

be given at Gauss points or element centres. When the physical quantities are output at Gauss points, the U function surface can be constructed by evaluating their nodal values through averaging the values at the surrounding Gauss points [23,24]. When the responses are output at element centres, their node values can be calculated by the Lagrange interpolation as discussed in Appendix A. 3.3.2. Filtering Topology optimization can be mesh dependent [17]. To solve this issue, filtering element sensitivity could be used for gradient-based topology optimization. In MIST [23,24], the U function can be filtered in a similar approach that in SIMP [17] and the U values are modified based on a linear weighted average within a spatial neighborhood. As values of the U function are calculated at nodal locations, the center of the spatial filter for a specific node is the nodal location and node-to-node distances of neighboring nodes are computed [23,24]. 3.3.3. Normalization To easily calculate the iso-surface level, the U function will be further normalized by:

U¼

U  Um Ua

ð8Þ

where

1 2

1 2

Um ¼ ½maxfUg þ minfUg; Ua ¼ ½maxfUg  minfUg

ð9Þ

By using Eqs. (8) and (9), the U function is normalized so that 1 6 U 6 1. 3.4. Step 4: Update of design variables or weighting factors 3.4.1. Incremental direction of the weighting factors When the U function is selected, calculation of the iso-surface level or threshold value is another key task in an iterative process. The iso-surface level of the chosen U function determined from the FEA results based on {xe}k at iteration k can be calculated using a bisection method by considering the volume constraint [23,24]. In this paper, a novel algorithm is developed to determine the iso-surface level t. In this algorithm, t may be simply found by counting nodal numbers based on the volume constraint. The details are given in Appendix B. An element with all nodal values of the U function above the iso-surface becomes solid (xe = 1) and that with all nodal values below it denotes void (xe ? 0); an element with the nodal U values above and below the iso-value plane, an area fraction with solid materials is the project area above it to this element area. Weighting factors {xe}k+1 are then obtained and incremental direction D{xe}k+1 (={xe}k+1  {xe}k) at the kth iteration is determined. 3.4.2. Incremental magnitude of the weighting factors Generally, {xe}k+1 should not be directly used for FEA at the next iterative step [23,24] owing to unstable or divergent iterations and thus move limit is introduced to reduce increments, similar to the move limit used in optimality criteria [14]. Constant or variable move limit may be used. In this study, the dynamic move limit given in [24] is used, where the move limit is halved when the objective function starts to oscillate. Oscillation may be avoided by choosing very small values of the move limit. However, this will slow down convergent rates and may even lead to a local minimum. In the iterative process, the weighing factors xe (e = 1, 2, . . . , Ne) must be in a range of (0 < xmin 6 xe 6 1) to avoid possible numerical singularity (xe < xmin) and unacceptable design (xe > 1). The

updating scheme for the weighting factors in MIST can be defined as:

Dfxe gkþ1 ¼ fxe gð1Þ kþ1  fxe gk

ð10aÞ

ð2Þ

fxe gkþ1 ¼ fxe gk þ kmv fDxe gkþ1

ð10bÞ

n o n o ð3Þ ð2Þ ð3Þ ð2Þ fxe gkþ1 ¼ max fxe gkþ1 ; xmin and fxe gkþ1 ¼ min fxe gkþ1 ; 1 ð10cÞ ð3Þ fxe gkþ1

where kmv is the move limit, 0 < kmv 6 1; weighting factors used for the next iteration.

are the updated

3.5. Step 5: Convergence test Convergence parameters may be defined in terms of the changes in the weighting factors, objective function and the response surface, etc. at element or overall domain level. The iterative process may be deemed as convergent when one or more convergence parameters are less than given tolerances. Setting the convergent tolerance is normally problem-dependent and the maximum iterative number (e.g., 200) can also be set to avoid the iteration without ending. When the convergence criteria are met, the iteration terminates otherwise go to Step 2. 4. Numerical results Numerical results are presented for plate structures in Fig. 1. Fig. 1(a) is an aluminum plate and Fig. 1(b) is an aluminum plate with bonded piezoelectric sheets on one side. The material properties are listed in Table 1 and dimensions are given in Fig. 1. The cantilever plate in Fig. 1 is often used as a part of real structures, in which bending and twisting shapes are two basic deformations, which may be defined by [1,29]:

wdb ðx; yÞ ¼

1  cosðx=aÞ Gb

wdt ðx; yÞ ¼

½coshðx=aÞ  1 sinðy=bÞ Gt

ðPlate bendingÞ

ðPlate twistingÞ

ð11aÞ

ð11bÞ

where a (= 240 mm) and b (= 200 mm) are the plate length and width; Gb and Gt are the scaling factors. Eqs. (11a) and (11b) may typically represent bending and twisting shapes [1,29]. In the present computation, shape morphing is conducted for the case of ‘desired shape specified’ [1]. The present optimization for shape morphing is to minimize the squared errors between the computed and desired shapes via optimal material distribution. In the present computations, the desired shapes are described by deflections at the chosen observation points as illustrated in Fig. 1. 4.1. A plate subjected to mechanical forces In this subsection, the MIST method is applied to find optimal topology for a plate of Fig. 1(a) subjected to mechanical forces to demonstrate the present MIST algorithm for shape morphing. In the computation, vfarc = 0.5 and p = 3. 4.1.1. Morphing bending shape When Gb = 0.4 (mm1) in Eq. (11a), desired deflections of the plate in Fig. 1(a) at x = 120 and 240 (mm) are 0.306 and 1.149 (mm). It is assumed that the cantilever plate is subjected to a transverse force F = (0, 0, 1.5) (N) at point E2 (240, 0) (mm).

94

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

y

M3

y

E3

M3

E3

M2

E2

M1

E1

73 73 PZT M2

E2

200

x

200

x

73 73 M1

E1

PZT

(a) A plate with isotropic material

(b) A PZT plate Not scaled to size : observation point

Fig. 1. (a) An aluminum plate (thickness = 2 mm); (b) an aluminum plate with PZT sheets bonded on one side (PZT thickness = 0.191; host plate thickness = 2 or 1 mm).

Table 1 Materials properties. Material

Piezoelectric ceramics PSI-5H4E

Aluminum 5000 series 5005-H14

Young’s modulus (N/m2) Poisson’s ratio Density (kg/m3) Piezoelectric constants d33 (m/v) Piezoelectric constant d31 (m/v) Coupling coefficient k33 Coupling coefficient k31 Polarizing field (v/m) Coercive field (v/m) Relative Dielectric Constant at 1000 (Hz) Operating voltage (v/mm)

6.2  1010 0.3 7800 650  1012

6.89  1010 0.33 2700

320  1012 0.75 0.44 >1.5  106 8  105 3800 +1181/472

The data are provided by the material suppliers.

The U function and iso-threshold-value surfaces for bending shape morphing at iterations 1, 10, 30 and 200 are depicted in Fig. 2(a)–(d). In MIST, an area with the U function value being larger than the threshold is treated as a solid (xe = 1) and that with the U function value being less than the threshold is void (xe = xmin ? 0). Intersection of the U function surface and the iso-surface forms the topological boundary contour at iterations 1, 10, 30 and 200 as shown in Fig. 3(a)–(d). The weighting factors of the elements in the boundary are calculated by the project area. The new weighting factors are then updated by using Eqs. 10a, 10b, 10c, which represent the material distributions as in SIMP. At iterations 1, 10, 30 and 200, material distributions are indicated in Fig. 4(a)–(d). The desired and computed shapes for morphing plate bending at these iterative steps are illustrated in Fig. 5(a)–(d). It is noted that the achieved shapes are computed by using the weighting factors at iterations 0, 9, 29 and 199. Fig. 3(d) is the optimal topology with a neat and smooth boundary. This structure is fabricated more easily than that of Fig. 4(d) as there is no grey material in the optimal topology. Fig. 5(d) shows that the desired bending shape is achieved. In the iterative computations, the error function can vary with iteration and its variation curve is plotted in Fig. 6 (Bending-F). At iterative steps 1 and 200, the error functions are 64.04 and 0.0205 (mm2), respectively. It can be seen that the squared error between the desired and computed shapes is convergent and the converged value is only 0.0321% of the initial value.

Fig. 2(c) and (d) show that it may be difficult to find iso-surface threshold by calculating an area with solid materials as the U function surface is very steep (c ? 0) in some regions. Therefore, an approach by counting node numbers is developed to find the iso-value surface level, as presented in Appendix B. 4.1.2. Morphing twisting shape When Gt = 0.3 (mm1) in Eq. (11b), desired deflections of the plate in Fig. 1(a) at points (120, ±100), (120/240, 0), (240, ±100) (mm) are ±0.204, 0, ±0.868 (mm), respectively. It is assumed that the plate is subjected to forces (0, 0, ±1.4) at points (240, ±100) (mm). The U function and iso-threshold value surfaces for the twisting shape morphing at iterations 1, 10, 30 and 200 are plotted in Fig. 7(a)–(d). The topologies at these iterations are directly obtained by intersection of the U function surface and the isosurface level t or the project of (U P t), as shown in Fig. 8(a)–(d). Fig. 9(a)–(d) illustrates the weighting factors or material distributions at iterations 1, 10, 30 and 200 and Fig. 10(a)–(d) are the desired and achieved twisting shapes at theses iterations. The curve of the error function versus iteration is also given in Fig. 6 (Twisting-F). At iterations 1 and 200, the errors are 11.03 and 0.0037 (mm2), respectively. In fact, the squared error is equal to 0.0037 (mm2) at iterations 75–200. Also, topologies and computed shapes are not changed in these iterative steps. That is, optimal design for this twisting shape morphing is obtained after 75 iterations. Fig. 7(b)–(d) also indicate that it may be difficult to determine the iso-surface level t by calculating an area with solid materials as the U function surface becomes rathe flat (c ? p/2) in the region near y = 0 and thus the new approach in Appendix B is used to find t effectively. 4.1.3. The numerical results computed by using NASTRAN for optimal topologies The optimal topologies with solid materials for morphing bending and twisting shapes shown in Figs. 3(d) and 8(d) are imported to FEA software MSC PATRAN/NASTRAN via software SolidWorks. As the topology figures are directly transferred to the FEA geometry used in NASTRAN, there are square (U P 0 at 4 nodes) and triangle elements (U P 0 at 1–3 nodes) in the FEA mesh, as shown in Fig. 11(a) and (b). Fig. 11(c) and (d) illustrate plate bending and twisting deformations and the load scenarios. It is noted that materials at some observation points have been removed in the optimal topologies but displacements at these points can be computed

95

1

Normalized phi-function

Normalized phi-function

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

Phi-function Iso-surface

0.5 0 -0.5 -1 100 0

Axis y (mm) -100 0

60

240 120 180

Axis x (mm)

0.5 0 -0.5 -1 100 0

Axis y (mm) -100 0

Phi-function Iso-surface

0.5 0 -0.5 -1 100 0

Axis y (mm)

-100 0

60

120

60

240 120 180

Axis x (mm)

(b) Niter. = 10 Normalized phi-function

Normalized phi-function

(a) Niter. = 1 1

Phi-function iso-surface

1

180

0.5 0 -0.5 -1 100 0

240

Axis x (mm)

Phi-function Iso-surface

1

Axis y (mm)

(c) Niter. = 30

-100 0

60

180

120

240

Axis x (mm)

(d) Niter. = 200

100

100

50

50

Axis y (mm)

Axis y (mm)

Fig. 2. Phi-functions and iso-surfaces for plate bending at iterations 1, 10, 30 and 200.

0

-50

-50

60

120

180

-100 0

240

120

Axis x (mm)

(a) Niter. = 1

(b) Niter. = 10

100

100

50

50

0

180200

0

-50

-50

-100

60

Axis x (mm)

Axis y (mm)

Axis y (mm)

-100 0

0

0

60

120

180

240

Axis x (mm)

(c) Niter. = 30

-100

0

60

120

180

240

Axis x (mm)

(d) Niter. = 200

Fig. 3. Topologies for the plate bending at iterations 1, 10, 30 and 200.

in the iterations as xmin > 0. Table 2 lists the specified deflections and the computed ones by using NASTRAN at the observation points. In Table 2, the data measured in the experiment are also given, which will be discussed in Section 5. In bending shape

morphing, the deflections at M1 and M3 predicted by NASTRAN are 4.9% less than the specified ones and the deflection at E2 is 7.6% larger than the specified one. In twisting shape morphing, the differences between the computed (using NASTRAN) and

96

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

(a) Niter. = 1

(b) Niter. = 10

(c) Niter. = 30

(d) Niter. = 200

Fig. 4. Material distributions for the plate bending at iterations 1, 10, 30 and 200.

(a) Niter. = 1

(c) Niter. = 30

(b) Niter. = 10

(d) Niter. = 200

Fig. 5. Specified and computed shapes of the plate bending at iterations 1, 10, 30 and 200.

the desired deflections at E1 and E3 are less than 1%. At M2, the deflection is 1.851  105 mm, which is 3 orders lower than that at (100, ±5) (mm). Fig. 11(c) and (d) indicate that the optimal

topologies obtained using MIST can be used for morphing plate bending and twisting shapes discussed in Sections 4.1.1 and 4.1.2.

97

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

and q = 3, the topology of piezoelectric materials is the same as that of the host material in the bonding region. Fig. 12(d) shows that the achieved shape correlates with the desired one. The normalized error function versus iteration is plotted in Fig. 6 (curve Bending-FV), which indicates the iteration is convergent. In the iterations 144–200, we note that Er = 0.01605 (mm2), which is less than 0.0205 (mm2). It is evident that the desired shape can be better achieved by applying both the electrical and mechanical loadings.

Normalized error (Ern /Er1 )

1

0.8

Bending-F Twisting-F Bending-FV Bending-V

0.6

0.4

0.2

4.3. A PZT plate subjected to an electrical field

0 0

20

40

60

80

100

Iteration Fig. 6. Error function versus iteration: Er1 – error at iteration 1; Ern – error at iteration n.

4.2. A PZT plate subjected to electrical and mechanical loadings

Phi-function Iso-surface

1

Normalized phi-function

Normalized phi-function

A plate with bonded piezoelectric sheets is shown in Fig. 1(b). Consider the same desired shape as that in Section 4.1. MIST will be used to find optimal design of the PZT plate subjected to electrical and mechanical loadings: V = 80 (v) and F = 1.3 (N). The host plate thickness of 2 (mm) is also considered. In all the iterations, we chose: vfarc = 0.5; p = 3; q = 3 and r = 1. Fig. 12(a) illustrates the U function and iso-surface at iteration 200. The optimal topology and the material distribution are given in Fig. 12(b) and (c), respectively, which are almost the same as those of the aluminum plate subjected to force F = 1.5 (N). This is understandable because the piezoelectric actuation is relatively small in this case. When different ratios of voltage V and force F are used, different optimal topologies can be obtained. As p = 3

0.5 0 -0.5 -1 100 0

Axis y (mm) -100 0

60

When a PZT plate in Fig. 1(b) is subjected to an electrical field only, the host plate with thickness of 1 (mm) is used due to the relatively lower piezoelectric actuation. As the maximum operation voltage is lower in anti-poling direction as listed in Table 1, a plate with bonded PZT sheets on one side is considered. The desired shape in Eq. (11a) is specified and Gb = 2/3 (mm1). At points Mi and Ei (i = 1, 2, 3), the desired deflections are 0.184 and 0.69 (mm), respectively. MIST is used to find the optimal topology for morphing bending shape of this PZT plate subjected to voltage of 150 (v). In all the iterations, we chose: vfarc = 0.6; p = 3; q = 0 and r = 1. The U function and iso-surface at iteration 200 are shown in Fig. 13(a). The optimal topology in Fig. 13(b) is directly obtained by projecting the U function surface with values of (U P 0) onto the iso-surface. The optimal material distribution is shown in Fig. 13(c). Piezoelectric materials will not be removed in the iterative processes as p = 3 and q = 0. By comparing the optimal materials distribution of Fig. 13(c) with that computed by using SIMP [32], the material distribution features obtained by the two methods (SIMP and MIST) are similar. The error function versus iteration is given in Fig. 6 (curve Bending-V) and Er = 0.01181 (mm2) in iterations 77–200, which clearly illustrates the convergence. Therefore, the desired bending shape is achieved. In the optimization for morphing a bending

120

180

240

Phi-function Iso-surface

1 0.5 0 -0.5 -1 100 0

Axis y (mm) -100 0

60

Phi-function Iso-surface

0 -0.5 -1 100 0

Axis y (mm) -100 0

60

120

180

240

Axis x (mm)

(c) Niter. = 30

240

(b) Niter. = 10 Normalized phi-function

Normalized phi-function

(a) Niter. = 1

0.5

180

Axis x (mm)

Axis x (mm)

1

120

Phi-function Iso-surface

1 0.5 0 -0.5 -1 100 0

Axis y (mm) -100 0

60

120

180

240

Axis x (mm)

(d) Niter. = 200

Fig. 7. Phi-functions and iso-surfaces for plate twisting at iterations 1, 10, 30 and 200.

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

100

100

50

50

Axis y (mm)

Axis y (mm)

98

0

-50

-50

-100 0

0

60

120

180

-100 0

240

60

Axis x (mm)

180

240

180

240

(b) Niter. = 10

100

100

50

50

Axis y (mm)

Axis y (mm)

(a) Niter. = 1

0

0

-50

-50

-100 0

120

Axis x (mm)

60

120

Axis x (mm)

(c) Niter = 30

180

240

-100 0

60

120

Axis x (mm)

(d) Niter = 200

Fig. 8. Topologies for the plate twisting at iterations 1, 10, 30 and 200.

(a) Niter. = 1

(b) Niter. = 10

(c) Niter. = 30

(d) Niter. = 200

Fig. 9. Material distributions for the plate twisting at iterations 1, 10, 30 and 200.

99

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

(a) Niter. = 1

(b) Niter. = 10

(c) Niter. = 30

(d) Niter. = 200

Fig. 10. Specified and computed shapes of the plate twisting at iterations 1, 10, 30 and 200.

(a) Mesh for bending topology

(b) Mesh for twisting topology

(c) Bending shape and load scenario

(d) Twisting shape and load scenario

Fig. 11. Nastran verification of the optimal topologies for morphing bending and twisting shapes of a plate subjected to mechanical forces.

Table 2 Deflections at observation points predicted by using NASTRAN and measured in the experiment. Shape

Observation point

Specified

Computed

Bending

M1 M3 E2

0.306 0.306 1.149

0.291 0.291 1.236

M2 E1 E3

0 0.868 0.868

1.851  105 0.863 0.863

Twisting

Plate deflection (mm)

Difference (%) Measured

Computed

Measured

0.32 1.18

4.9 4.9 7.57

4.79 2.28

– 0.576 0.576

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

100

Phi-function Iso-surface

1 0.5

Axis y (mm)

Normalized phi-function

100

0 -0.5 -1 100

0 Axis y (mm) -100 0

120 180 240

60

50 0 -50

-100

Axis x (mm)

(a) The

0

60

120

180

240

Axis x (mm)

(b) Optimal topology (Niter.= 200)

function (Niter.= 200)

(d) Specified & computed shapes (Niter.= 200)

(c) Material distribution (Niter.= 200)

Phi-function Iso-surface

1 0.5

50

Axis y (mm)

Normalized phi-function

Fig. 12. Optimal design for morphing bending plate subjected force and electrical field in piezoelectric materials.

0 -0.5 -1 100

-50 0

Axis y (mm) -100 0

60

120

180

240 0

Axis x (mm)

(a) The

0

120

180

240

Axis x (mm)

(b) Optimal topology (Niter.= 200)

function (Niter.= 200)

(c) Material distribution (Niter.= 200)

60

(d) Specified & computed shapes (Niter.= 200)

Fig. 13. Optimal design for morphing bending plate using fixed piezoelectric materials with one channel voltage.

shape of the PZT plate studied in this section, one channel voltage is used and piezoelectric materials are not removed, which would be relatively easy to realize in practice. 5. Experimental test 5.1. Test setup Optimal topologies obtained for plate bending in Section 4 are fabricated and experimental tests are conducted. Fig. 14(a)–(d)

indicates three specimens subjected to a mechanical force, a combination of mechanical and electrical loads and an electrical field for morphing plate bending shapes. Fig. 15 shows the experimental instruments and setup. In Fig. 15, the DC power supply drives the displacement sensors that are used to measure plate defections; the data logger transfers the measuring signals to computer; the signal amplifier generates electric voltage applied to piezoelectric device and the signal control is used to control electric voltages. The mechanical force is applied by hanging weight shown in Fig. 14(b) and (c). The measurement is

101

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

(a) A plate subjected to mechanical loading

(b) A plate with electrical & mechanical loads (-)

(c) A plate with electrical & mechanical loads (+)

(d) A plate subjected to electrical loading

Fig. 14. Specimens (a) optimal topology of the plate subjected to mechanical loading; (b) and (c) optimal topology of the plate subjected to mechanical and electrical loads; (d) optimal topology subjected to electrical loading.

4 6

2 3 5

1 Fig. 15. Experimental setup: (1) anti-vibration table (KS2210); (2) DC power supply (BK PECISION 1761); (3) data logger (NI cDAQ 9172/9221); (4) signal amplifier (PI HVPZT); (5) displacement sensor (Indikon 100–400 Series); (6) signal control (PCI-6723 CB-68LPR).

programmed by using LabView 8.0. In the experimental tests, the deflections at the selected observation points are measured for the specimens in Fig. 14(a)–(d).

Deflection (mm)

1.5 M3 test 1 M3 test 2 M3 test 3 E2 test 1 E2 test 2 E2 test 3

1

5.2. Testing results and comparison

0.5

0 0

2

4

6

8

10

Time (s) Fig. 16. The measured deflections at points M3 and E2 for the specimen of Fig. 14(a) in the 3 tests.

The deflection versus time for the specimen of Fig. 14(a) predicted by experiment is illustrated in Fig. 16, where the measured data in 3 repeated tests are given. The average deflections at points M3 and E2 are listed in Table 2, which are 4.79% and 2.28% higher than the specified ones. It is noted that the specified values may be deemed as the computed ones in the iteration because the error function has been minimized, as illustrated in Figs. 5(d), 10(d), 12(d) and 13(d). The measured values are slightly less than those predicted by NASTRAN as the force is not applied at the free edge. In the iteration, there are three types of materials in elements: solid, grey and void ones, whilst only the solid material

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

is considered in the experiment and the FEA by using NASTRAN. Therefore, there exists a correlation between the present computations using MIST and measured results. In the present optimization for the piezoelectric plate subjected to both electrical and mechanical loadings, the PZT material distribution is the same as that of host material in the bonding region. Fig. 14(b) and (c) show that some extra piezoelectric materials are kept in the specimen. Since these materials are less than 1% of the total piezoelectric materials and it is hard to cut piezoelectric sheets into a complex shape, they are left in the specimen. When the specimen of Fig. 14(b) is used, an electrical field should be applied in the anti-poling direction so that the deflection directions caused by the electrical and mechanical loadings are the same. When the specimen of Fig. 14(c) is tested, the applied voltage is in the poling direction. The measured deflections at points M3 and E2 for specimens of Fig. 14(b) and (c) are illustrated in Fig. 17, which show that the results are almost the same for the two loading cases. In the experiment, the weight is hanged firstly and then voltage is applied. The measured deflections at points M3 and E2 for specimens of Fig. 14(b) are 0.32 and 1.20, and those for Fig. 14(c) are 0.33 and 1.19 (mm), respectively, which correlate with the specified/computed deflections. The measured data in the 3 tests for the specimen of Fig. 14(d) are given in Fig. 18. The average deflections at points M3 and E2 are 0.18 and 0.7 (mm). The differences between the measured and

Deflection (mm)

0.2

0.1

0 0

30

LS = 6 (v/s) LS = 10 (v/s) LS = 15 (v/s) LS = 30 (v/s) LS = 150 (v/s) Switch on/off Computed

60

90

-0.1

Time (s) Fig. 19. The measured deflection at point M2 of the specimen in Fig. 14(d) for different loading speeds. LS: electrical loading/unloading speed.

specified/computed results are less than 5%. Figs. 16–18 show that the present optimal designs can be effectively used for morphing plate shapes for each considered case. 5.3. Overshoot and quasi-static state of piezoelectric actuation Deflection overshoots are observed in Figs. 17 and 18 where the voltages are suddenly applied by or equivalent to ‘switch on’. The overshoot can be harmful to shape morphing for piezoelectric structures, particularly for dynamic shape tracking. In static analysis, quasi-static actuation is generally required. Therefore, this phenomenon is further measured in the present experiments. Fig. 19 illustrates deflections at point M2 of the specimen of Fig. 14(d) subjected to Vmax = 150 (v) with different loading/unloading rates, where the desired/computed value is also marked. When voltage is abruptly loading (0 ? Vmax) or unloading (Vmax ? 0), overshoots are significant. When the loading rate is less than 15 (v/s) or unloading one is less than 10 (v/s), the overshoots are not observed. It is also observed that overshooting occurs more easily in unloading than in loading.

1.5

1 M3(+) E2(+) M3(-) E2(-)

0.5

0.3

Deflection (mm)

102

0 0

20

40

60

80

100

6. Conclusion

Time (s) Fig. 17. The measured deflections at points M3 and E2 for the specimens of Fig. 14(b) and (c). (+): the applied voltage in poling direction; (): the applied voltage in anti-poling direction.

1

M3 test 1 M3 test 2 M3 test 3 E2 test 1 E2 test 2 E2 test 3

Deflection (mm)

0.8

0.6

0.4

0.2

The novel formulations and algorithms of topology optimization using MIST for morphing structural shapes are derived and applied to piezoelectric plates. Numerical computations indicate that the present response function based on mutual strain energy density is effective to minimize the error norm and the present algorithm based on the nodal count is efficient to determine the iso-surface threshold value. By using MIST, black-white designs with clear boundary contours are obtainable, and thus their conversion to practically manufacturable products becomes possible. Numerical results and experimental testing indicate that the present MIST algorithm can yield optimal topologies of piezoelectric plates for shape morphing by applying a mechanical force, one channel electrical voltage and combination of mechanical and electrical loadings. The present study on the overshoot is helpful for precise shape morphing of piezoelectric structures.

Acknowledgment

0 0

2

4

6

8

10

Time (s) Fig. 18. The measured deflections at points M3 and E2 for the specimen of Fig. 14(c) in the 3 tests.

The authors would like to acknowledge the support of the Australian Research Council via Discovery-Project Grants (DP110104123 and DP140104408).

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

Appendix A. Calculation of the U function The U function plays a critical role in MIST method [23,24], which is a surface constructed by physical response values at element nodes. When this method is used for shape morphing, the U function is defined by mutual energy densities at nodes, which can be calculated by Eq. (5). When Eq. (5) is used, mutual strain energy density at (n, g) in the eth element can be calculated by:

1 e T fu g fBgT ½D½Bfuek g 2 i 1 T ¼ fuek g fBgT ½D½Bfuei g 2

ðeÞ

ðEmusd Þi ðn; gÞ ¼

ðA1Þ

where fuei g and fuek g are the displacement vectors of the eth element caused by a unit dummy force at point i and the applied loading; ½B and ½D are the strain and elastic matrix. By using Eq. (A1), mutual strain energy densities at the Gauss points can be calculated element by element and the U function can be constructed by evaluating its nodal value by averaging the values at the surrounding Gaussian points [23,24]. Alternatively, the mutual strain energy density can also be found by using the element stiffness matrix and node displacements. Mutual strain energy of the eth element caused by a unit dummy force at point i can be calculated by: ðeÞ

ðEmus Þi ¼

Z Xe

1 e T T fu g ½B ½D½Bfuek gdXe 2 i

ðA2Þ

Mutual strain energy density at the eth element centre is approximately given by: ðeÞ

ðEmusd Þi ð0; 0Þ ¼

1 T fue g ½ke fuek g 2v e i

ðA3Þ

where ve and [ke] are the volume and stiffness matrix of the eth element. By using Eq. (A3), mutual energy density at the element centre can be readily calculated element by element and its values at nodes can be found by using Lagrangian interpolation method and then U function can be constructed. In the present finite element analysis, 4-node elements are used to mesh the design domain and the 2nd order interpolation functions is used to find the node values over the 3 adjacent elements. Appendix B. Calculation of the iso-surface level The iso-surface level can be found by an iterative bisection method where the volume fraction is calculated in each iteration [23,24]. If the U function surface in some region is steep (see Section 4.1.1), the same volume fraction may be obtained in a wide range of threshold values. When the U function in some areas becomes very flat, a volume fraction may be hard to find (see Section 4.1.2). Angle c may be defined to describe the situation:

1

c ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 0 2 xÞ

1 þ ðU

þ ðU0y Þ

2

ðB1Þ

In the iterative process, when the iso-surface intersects the U function surface with (c ? 0) or (c ? p/2), the iso-surface threshold value would be difficult to be determined. To overcome the difficulties, the iso-surface level will be firstly found by counting nodes firstly in the present computations. The following MatLab code can be used to approximately find the iso-surface level t using node values of the U function:

103

a1 = 1.0; b1 = 1.0; % 1 6 U 6 1 after normalization using Eqs. (8) and (9) IsoL = (a1 + b1)/2; count = 0; for ib = 1:length(PhiF) if PhiF(ib) > IsoL count = count + 1; end end toler = 1.0; while toler > 0.0008; % tolerance, relying on the node number and volume constraint if count/length(PhiF) > volfrac b1 = IsoL; else a1 = IsoL; end IsoL = (a1 + b1)/2; count = 0; for ib = 1:length(PhiF) if PhiF(ib) > IsoL count = count + 1; end end toler = abs(count/length(PhiF)-volfrac)/volfrac; end

The found iso-surface level tn may satisfy the volume constraint when the mesh is fine enough. If not, further iterations to calculate the volume fraction as in [23,24] are required. In this case, the upper and low bounds (tn ± e) can be set, e.g., e = 0.01. The iso-surface level can be readily determined in (tn  e) < t < (tn ± +e) by using the bisection method in [23,24]. References [1] Koconis DB, Kollar LP, Springer GS. Shape control of composite plates and shells with embedded actuators. II. Desired shape specified. J Compos Mater 1994;28(3):262–85. [2] Irschik H. A review on static and dynamic shape control of structures by piezoelectric actuation. Eng Struct 2002;24(1):5–11. [3] Luo QT, Tong LY. High precision shape control of plates using orthotropic piezoelectric actuators. Finite Elem Anal Des 2006;42(11):1009–20. [4] Irschik H, Nader M. Actuator placement in static bending of smart beams utilizing Mohr’s analogy. Eng Struct 2009;31(8):1698–706. [5] Liu S, Tong L, Lin Z. Simultaneous optimization of control parameters and configurations of PZT actuators for morphing structural shapes. Finite Elem Anal Des 2008;44(6–7):417–24. [6] Agrawal BN, Treanor KE. Shape control of a beam using piezoelectric actuators. Smart Mater Struct 1999;8(6):729. [7] Sun D, Tong L. Design optimization of piezoelectric actuator patterns for static shape control of smart plates. Smart Mater Struct 2005;14(6):1353. [8] Mukherjee A, Joshi S. Piezoelectric sensor and actuator spatial design for shape control of piezolaminated plates. AIAA J 2002;40(6):1204–10. [9] Kögl M, Silva ECN. Topology optimization of smart structures: design of piezoelectric plate and shell actuators. Smart Mater Struct 2005;14(2):387. [10] Uysal H, Gul R, Uzman U. Optimum shape design of shell structures. Eng Struct 2007;29(1):80–7. [11] Carbonari RC, Silva ECN, Nishiwaki S. Optimum placement of piezoelectric material in piezoactuator design. Smart Mater Struct 2007;16(1):207. [12] Kang Z, Tong L. Topology optimization-based distribution design of actuation voltage in static shape control of plates. Comput Struct 2008;86(19– 20):1885–93. [13] Liu T, Veidt M, Kitipornchai S. Modelling the input–output behaviour of piezoelectric structural health monitoring systems for composite plates. Smart Mater Struct 2003;12(5):836–44. [14] Bendsoe MP, Kikuchi N. Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988;71(2):197–224. [15] Bendsøe MP. Optimal shape design as a material distribution problem. Struct Optim 1989;1(4):193–202.

104

Q. Luo, L. Tong / Engineering Structures 97 (2015) 90–104

[16] Zhou M, Rozvany GIN. The COC algorithm, Part II: Topological, geometrical and generalized shape optimization. Comput Methods Appl Mech Eng 1991;89(1– 3):309–36. [17] Bendsoe MP, Sigmund O. Topology optimization: theory, methods and applications. Berlin, New York: Springer; 2003. [18] Xie YM, Steven GP. A simple evolutionary procedure for structural optimization. Comput Struct 1993;49(5):885–96. [19] Querin OM, Steven GP, Xie YM. Evolutionary structural optimisation (ESO) using a bidirectional algorithm. Eng Comput 1998;15(8):1031–48. [20] Sethian JA. Fast marching methods. Siam Rev 1999;41(2):199–235. [21] Osher SJ, Santosa F. Level set methods for optimization problems involving geometry and constraints. I. Frequencies of a two-density inhomogeneous drum. J Comput Phys 2001;171(1):272–88. [22] Wang MY, Wang XM, Guo DM. A level set method for structural topology optimization. Comput Methods Appl Mech Eng 2003;192(1–2):227–46. [23] Tong LY, Lin JZ. Structural topology optimization with implicit design variableoptimality and algorithm. Finite Elem Anal Des 2011;47(8):922–32. [24] Vasista S, Tong LY. Design and testing of pressurized cellular planar morphing structures. AIAA J 2012;50(6):1328–38.

[25] Saxena A, Ananthasuresh GK. On an optimal property of compliant topologies. Struct Multidiscip Optim 2000;19(1):36–49. [26] Allaire G, Jouve F, Toader A-M. A level-set method for shape optimization. Compt Rend Math 2002;334(12):1125–30. [27] Tovar A, Khandelwal K. Topology optimization for minimum compliance using a control strategy. Eng Struct 2013;48:674–82. [28] Bourdin B, Chambolle A. Design-dependent loads in topology optimization. ESAIM: Contr Op Ca Va 2003;9:19–48. [29] Chee C, Tong L, Steven GP. Piezoelectric actuator orientation optimization for static shape control of composite plates. Compos Struct 2002;55(2):169–84. [30] Sigmund O. On the design of compliant mechanisms using topology optimization. Mech Struct Mach 1997;25(4):493–524. [31] Canfield S, Frecker M. Topology optimization of compliant mechanical amplifiers for piezoelectric actuators. Struct Multidiscip Optim 2000;20(4):269–79. [32] Luo Z, Luo QT, Tong LY, Gao W, Song CM. Shape morphing of laminated composite structures with photostrictive actuators via topology optimization. Compos Struct 2011;93(2):406–18.