# Preliminary Study of Optimum Piezoelectric Cross-Ply Composites for Energy Harvesting

Authors: David N. Betts; H. Alicia Kim; Christopher R. Bowen
Journal: Smart Materials Research
 (2012)

---

## Abstract

Energy harvesting devices based on a piezoelectric material attached to asymmetric bistable laminate plates have been shown to exhibit high levels of power extraction over a wide range of frequencies. This paper optimizes for the design of bistable composites combined with piezoelectrics for energy harvesting applications. The electrical energy generated during state-change, or “snap-through,” is maximized through variation in ply thicknesses and rectangular laminate edge lengths. The design is constrained by a bistability constraint and limits on both the magnitude of deflection and the force required for the reversible actuation. Optimum solutions are obtained for differing numbers of plies and the numerical investigation results are discussed.


---

## Body

## 1. Introduction

Energy harvesting which converts ambient mechanical vibrations into electrical energy is an area of considerable research interest and has received extensive attention in the past decade. A variety of methods have been considered including inductive, capacitive, and piezoelectric materials [1–3]. In many cases harvesting devices have been designed to operate at resonance to optimize the power generation, for example, simple linear cantilever beam configurations. However, ambient vibrations generally exhibit multiple time-dependent frequencies which can include components at relatively low frequencies. This can make typical linear systems inefficient or unsuitable; particularly if the resonant frequency of the device is higher than the frequency range of the vibrations it is attempting to harvest. In order to improve the efficiency of vibrational energy harvesters, recent work has focused on exploiting nonlinearity for broadband energy harvesting. Encouraging results [2] have been obtained using nonlinear or bistable cantilevered beams. Stanton et al. [2] modeled and experimentally validated a non-linear energy harvester using a piezoelectric cantilever. An end magnet on the oscillating cantilever interacts with oppositely poled stationary magnets, which induces softening or hardening into the system and allows the resonance frequency to be tuned. This technique was shown to outperform linear systems when excited by varying frequencies. However, such a system would require an obtrusive arrangement of external magnets and could generate unwanted electromagnetic fields. An alternative method has been recently found where a piezoelectric element is attached to bistable laminate plates with 2n plies and a total (T) layup of [0n/90n]T to induce large amplitude oscillations [3]. Such harvesting structures have been shown to exhibit high levels of power extraction over a wide range of frequencies. This arrangement can be designed to occupy a smaller space and is potentially more convenient and portable for broadband energy harvesting.Bistable composites have been extensively studied for the development of morphing or adaptive structure concepts [4–7]. When a composite laminate has an asymmetric stacking sequence the resulting mismatch in thermal expansion coefficients between plies leads to a thermally induced strain. This leads to the laminate developing a curved deformation as it is cooled from its high temperature cure cycle. Under certain geometric conditions the thermal strain can lead to the development of two stable equilibrium states (“bistability”). Such structures are of interest for shape-change applications since the “snap-through” between stable-states results in a large deflection that does not require continuous energy input to be maintained. Figure 1 shows an example of this behavior for a square [0n/90n]T laminate. For a low ratio of edge length to thickness only a saddle-shaped single stable state is observed, point A, with x- and y-curvatures of equal magnitude in opposite out-of-plane directions. As the ratio increases the solution bifurcates, point B. Beyond this point, two approximately cylindrical stable states are observed, points C and D, while the saddle state becomes unstable (dashed line). When a flexible piezoelectric material is attached to the laminate surface and the structure is repeatedly actuated between the stable states, the large shape changes have the potential to generate electrical energy by the direct piezoelectric effect.Figure 1
Stable (solid) and unstable (dashed) shapes of a[0n/90n]T laminate.Betts et al. [8] presented optimization of bistable laminates by deriving an analytical solution to modeling [0n/90n]T as well as non-[0n/90n]T layups. The primary focus of the paper was to develop shape changing structures whereby a piezoelectric strain was used to induce snap-through of the bistable laminate between stable states 1 and 2. By variation of ply orientations and using a nonuniform laminate geometry, laminates were optimized for differing stiffness characteristics whilst satisfying bistability and deflection constraints. This work was extended to include piezoelectric layers into the existing analytical model for the design of laminates constrained by actuation voltage limits [9].This paper adapts this work and introduces optimization of bistable laminates combined with piezoelectric material for energy harvesting applications. The electrical energy generated during state-change (e.g., C to D in Figure1) is maximized through variation in ply thicknesses and rectangular laminate edge lengths. A bistability constraint and limits on both the magnitude of deflection and the force required for the reversible actuation are specified to find a practical and achievable solution. Optimum solutions are obtained for differing numbers of plies and the numerical investigation results are discussed.


## 2. Design Problem

### 2.1. Problem Formulation

In this section we introduce the optimization formulation for the design of a bistable composite laminate for maximum piezoelectric energy harvesting capability as follows.Maximize: Electrical energy generated by two orthogonal piezoelectric layers.Subject to: The piezoelectric laminate must be bistable.The stroke of actuation,wdef, must be less than a maximum value.The forceF, required for state change from state 1 to state 2, and the reverse change from state 2 to state 1 must be less than a maximum specified value.Variables: Piezoelectric layer and laminate ply orientations,θp and θ respectively.Nonuniform ply thicknesses.Aspect ratio,Lx/Ly, of the laminate.


### 2.2. Laminate Configuration

The laminate geometry has an orthogonal layup adhering to the following set of design rules as illustrated in Figure2:(1)
two piezoelectric layers, poled in the longitudinal direction (referred to as thed33 direction), on the top and bottom surfaces of the laminate positioned 90° apart;(2)
an even number of plies, with pairs of plies about the laminate midplane positioned 90° apart;(3)
equal thicknesses,tp, for both piezoelectric layers;(4)
equal ply thicknesses about the laminate midplane,t1,t2,…,tn;(5)
nonuniform edge lengths,Lx and Ly, common to both the laminate and piezoelectric layers;(6)
each ply of the laminate is made from the same material.Figure 2
An example 4-ply piezoelectric (grey) and laminate (white) geometry.The laminate geometry above is selected for three reasons. Firstly the piezoelectric layers are positioned orthogonally to align with the major curvatures of the two stable shapes, maximizing the stress in the piezoelectric elements and, hence, the mechanical to electrical energy conversion. Secondly the laminate stacking sequence allows maximizing the useful deflection between states while providing scope for tailoring the directional properties. Finally, the restricted orthogonal stacking sequence reduces the design search space by excluding general asymmetric layups which do not exhibit the high deflection bistable characteristics. It is noted that the laminate (excluding the piezoelectric layers) can be defined in terms of plies above the midplane, referred to here as the half-laminate, (θ1 and θ2 in Figure 2).We initially consider a rectangular cross-ply design,[0n/90n]T, as studied experimentally by Arrieta et al. [3] This laminate exhibits equal out-of-plane displacement at all four corners and is actuated in the configuration shown in Figure 3. The laminate is held at all four corners (zero z-displacement). A mechanical actuator is attached to the centre point of the laminate surface such that displacement applied in the z-direction (out-of-plane) induces a state-change, simulating the ambient mechanical motion. The out-of-plane displacement, referred to here as the “stroke of actuation,” is the change in maximum out-of-plane deflection from one unloaded stable state to the second, wdef, measured at the geometric centre of the laminate.Figure 3
Arrangement for actuation of a rectangular cross-ply laminate.


### 2.3. Objective Function

The objective of this study is to maximize the electrical energy generated by the mechanical deflection of a bistable piezoelectric-composite laminate. We model this based on the static states of the system. The objective function is thus the electrical energy,U, induced by stress in the piezoelectric material induced in each stable state:
(1)U=12CV2,
where C is the capacitance and V is the open circuit output voltage, as defined by the following equations:(2)C=QV=dijσAV,V=Et=-gijσtp,
where dij is the effective piezoelectric strain constant (charge per unit force), gij is the effective piezoelectric voltage constant (electric field per unit stress), σ is the stress, A is the surface area of the piezoelectrics, and tp is the thickness of the layer. Substituting (2) into (1) gives:(3)U=12(dijgij)⋅σ2⋅(Atp),
where the terms have been split into the fixed material properties of the piezoelectric (dij and gij), the stress in the material which is a function of the laminate curvature, and the material geometry (A and tp). When attached to the surface of a curved laminate the stress in the piezoelectric material varies across the surface area as a function of the strains induced by the mismatched thermal properties. The piezoelectric material, poled in the longitudinal direction in this work, is strained in both the longitudinal direction (referred to as the 33 direction) and the transverse direction (31 direction). We therefore consider the contributions of both components. The stresses in the x- and y-directions are defined by the following equations,(4)σx=Q¯11εx+Q¯12εy,σy=Q¯12εx+Q¯22εy,
where the Q̅’s are transformed stiffness terms dependent on the material orientations and ε’s are the in-plane strains. These values vary across the x-y plane and are therefore integrated over the piezoelectric surface area and substituted into (3). We then consider the total electrical energy induced by a single stable state to be the sum of the longitudinal and transverse components of both the top and bottom piezoelectric layers. Each of these terms uses the associated 33 and 31 material properties as illustrated in Figure 4.Figure 4
Longitudinal and transverse components of the total electrical energy,U. Note: subscript denotes the associated direction; superscript denotes the top (t) and bottom (b) piezoelectric layer.The total electrical energy, defining the objective function is therefore, the sum of the four terms associated with state 1, and the four terms associated with state 2, (5):(5)Maximize(Uxb+Uxt+Uyb+Uyt)state1+(Uxb+Uxt+Uyb+Uyt)state2.


### 2.4. Constraints on the Design

#### 2.4.1. Bistability

The asymmetry of laminates does not guarantee bistability. When a laminate geometry has an edge length to thickness ratio below a critical value there exists only a single saddle solution [10]. This critical ratio is a function of the stacking sequence, temperature, material properties, and aspect ratio. Furthermore, the critical value is increased by the presence of the piezoelectric layers due to their added stiffness. We introduce a constraint to ensure bistability. Noting that the x- and y-curvatures of the monostable saddle shape have equal magnitude in opposite out-of-plane directions, the sum of these curvatures is zero. Beyond the bifurcation point, see Figure 1, the sum of the curvatures is either positive or negative depending on the associated cylindrical state. We therefore impose the constraint a+b>0, where a and b are the x- and y-curvature coefficients, respectively, such that only bistable composites are considered feasible:(6)Bistability constraint:a+b>0.Close to the bifurcation point the solution is highly sensitive to manufacturing imperfections and uncertainties in material data [11]. This issue can easily be addressed by adding a small constant to the constraint, that is, a+b>X, where X is a “distance” away from the bifurcation point. Our previous studies have shown that this constraint could also be used as a minimum out-of-plane deflection constraint [8].


#### 2.4.2. Stroke of Mechanical Actuation

The stroke of actuation refers to the distance travelled by the mechanical actuator in thez-direction, Figure 3 and a practical limit is imposed. This is modeled by fixing the displacement at the centre with the maximum deflection being at the corners of the laminate. Noting that the out-of-plane displacement, w of the curved laminate shape is approximated by quadratic polynomial of (7):(7)w(x,y)=12(ax2+by2+cxy)
the corner deflection between states, wdef, at x=Lx/2, y=Ly/2, can be expressed as (8) where the subscript denotes the associated stable state.(8)wdef=w1-w2=18(a1Lx2+b1Ly2+c1LxLy)-18(a2Lx2+b2Ly2+c2LxLy).Due to the orthogonality of the laminate the deflection between states is the same at all four corners. Furthermore, the curvaturesa, b, and c can be expressed in terms of the first state alone due to the two states having curvatures of equal magnitude. The constraint on the stroke, less than wmax, can therefore be reduced to (9):(9)Actuation stroke constraint:18(a1+b1)(Lx2+Ly2)<wmax⁡.


#### 2.4.3. Force of Mechanical Actuation

Similarly, the mechanical force required for actuation is constrained. The force for actuation is considered separately for the transition from state 1 to state 2,F1, and the reverse from state 2 to state 1, F2. These values are not necessarily equal for nonsquare laminates which can result in two shapes with greatly different out-of-plane displacement. The force constraint ensures that both of these values are below the available maximum force, Fmax:(10)Actuation force constraint:F1<Fmax⁡,F2<Fmax⁡.


## 2.1. Problem Formulation

In this section we introduce the optimization formulation for the design of a bistable composite laminate for maximum piezoelectric energy harvesting capability as follows.Maximize: Electrical energy generated by two orthogonal piezoelectric layers.Subject to: The piezoelectric laminate must be bistable.The stroke of actuation,wdef, must be less than a maximum value.The forceF, required for state change from state 1 to state 2, and the reverse change from state 2 to state 1 must be less than a maximum specified value.Variables: Piezoelectric layer and laminate ply orientations,θp and θ respectively.Nonuniform ply thicknesses.Aspect ratio,Lx/Ly, of the laminate.


## 2.2. Laminate Configuration

The laminate geometry has an orthogonal layup adhering to the following set of design rules as illustrated in Figure2:(1)
two piezoelectric layers, poled in the longitudinal direction (referred to as thed33 direction), on the top and bottom surfaces of the laminate positioned 90° apart;(2)
an even number of plies, with pairs of plies about the laminate midplane positioned 90° apart;(3)
equal thicknesses,tp, for both piezoelectric layers;(4)
equal ply thicknesses about the laminate midplane,t1,t2,…,tn;(5)
nonuniform edge lengths,Lx and Ly, common to both the laminate and piezoelectric layers;(6)
each ply of the laminate is made from the same material.Figure 2
An example 4-ply piezoelectric (grey) and laminate (white) geometry.The laminate geometry above is selected for three reasons. Firstly the piezoelectric layers are positioned orthogonally to align with the major curvatures of the two stable shapes, maximizing the stress in the piezoelectric elements and, hence, the mechanical to electrical energy conversion. Secondly the laminate stacking sequence allows maximizing the useful deflection between states while providing scope for tailoring the directional properties. Finally, the restricted orthogonal stacking sequence reduces the design search space by excluding general asymmetric layups which do not exhibit the high deflection bistable characteristics. It is noted that the laminate (excluding the piezoelectric layers) can be defined in terms of plies above the midplane, referred to here as the half-laminate, (θ1 and θ2 in Figure 2).We initially consider a rectangular cross-ply design,[0n/90n]T, as studied experimentally by Arrieta et al. [3] This laminate exhibits equal out-of-plane displacement at all four corners and is actuated in the configuration shown in Figure 3. The laminate is held at all four corners (zero z-displacement). A mechanical actuator is attached to the centre point of the laminate surface such that displacement applied in the z-direction (out-of-plane) induces a state-change, simulating the ambient mechanical motion. The out-of-plane displacement, referred to here as the “stroke of actuation,” is the change in maximum out-of-plane deflection from one unloaded stable state to the second, wdef, measured at the geometric centre of the laminate.Figure 3
Arrangement for actuation of a rectangular cross-ply laminate.


## 2.3. Objective Function

The objective of this study is to maximize the electrical energy generated by the mechanical deflection of a bistable piezoelectric-composite laminate. We model this based on the static states of the system. The objective function is thus the electrical energy,U, induced by stress in the piezoelectric material induced in each stable state:
(1)U=12CV2,
where C is the capacitance and V is the open circuit output voltage, as defined by the following equations:(2)C=QV=dijσAV,V=Et=-gijσtp,
where dij is the effective piezoelectric strain constant (charge per unit force), gij is the effective piezoelectric voltage constant (electric field per unit stress), σ is the stress, A is the surface area of the piezoelectrics, and tp is the thickness of the layer. Substituting (2) into (1) gives:(3)U=12(dijgij)⋅σ2⋅(Atp),
where the terms have been split into the fixed material properties of the piezoelectric (dij and gij), the stress in the material which is a function of the laminate curvature, and the material geometry (A and tp). When attached to the surface of a curved laminate the stress in the piezoelectric material varies across the surface area as a function of the strains induced by the mismatched thermal properties. The piezoelectric material, poled in the longitudinal direction in this work, is strained in both the longitudinal direction (referred to as the 33 direction) and the transverse direction (31 direction). We therefore consider the contributions of both components. The stresses in the x- and y-directions are defined by the following equations,(4)σx=Q¯11εx+Q¯12εy,σy=Q¯12εx+Q¯22εy,
where the Q̅’s are transformed stiffness terms dependent on the material orientations and ε’s are the in-plane strains. These values vary across the x-y plane and are therefore integrated over the piezoelectric surface area and substituted into (3). We then consider the total electrical energy induced by a single stable state to be the sum of the longitudinal and transverse components of both the top and bottom piezoelectric layers. Each of these terms uses the associated 33 and 31 material properties as illustrated in Figure 4.Figure 4
Longitudinal and transverse components of the total electrical energy,U. Note: subscript denotes the associated direction; superscript denotes the top (t) and bottom (b) piezoelectric layer.The total electrical energy, defining the objective function is therefore, the sum of the four terms associated with state 1, and the four terms associated with state 2, (5):(5)Maximize(Uxb+Uxt+Uyb+Uyt)state1+(Uxb+Uxt+Uyb+Uyt)state2.


## 2.4. Constraints on the Design

### 2.4.1. Bistability

The asymmetry of laminates does not guarantee bistability. When a laminate geometry has an edge length to thickness ratio below a critical value there exists only a single saddle solution [10]. This critical ratio is a function of the stacking sequence, temperature, material properties, and aspect ratio. Furthermore, the critical value is increased by the presence of the piezoelectric layers due to their added stiffness. We introduce a constraint to ensure bistability. Noting that the x- and y-curvatures of the monostable saddle shape have equal magnitude in opposite out-of-plane directions, the sum of these curvatures is zero. Beyond the bifurcation point, see Figure 1, the sum of the curvatures is either positive or negative depending on the associated cylindrical state. We therefore impose the constraint a+b>0, where a and b are the x- and y-curvature coefficients, respectively, such that only bistable composites are considered feasible:(6)Bistability constraint:a+b>0.Close to the bifurcation point the solution is highly sensitive to manufacturing imperfections and uncertainties in material data [11]. This issue can easily be addressed by adding a small constant to the constraint, that is, a+b>X, where X is a “distance” away from the bifurcation point. Our previous studies have shown that this constraint could also be used as a minimum out-of-plane deflection constraint [8].


### 2.4.2. Stroke of Mechanical Actuation

The stroke of actuation refers to the distance travelled by the mechanical actuator in thez-direction, Figure 3 and a practical limit is imposed. This is modeled by fixing the displacement at the centre with the maximum deflection being at the corners of the laminate. Noting that the out-of-plane displacement, w of the curved laminate shape is approximated by quadratic polynomial of (7):(7)w(x,y)=12(ax2+by2+cxy)
the corner deflection between states, wdef, at x=Lx/2, y=Ly/2, can be expressed as (8) where the subscript denotes the associated stable state.(8)wdef=w1-w2=18(a1Lx2+b1Ly2+c1LxLy)-18(a2Lx2+b2Ly2+c2LxLy).Due to the orthogonality of the laminate the deflection between states is the same at all four corners. Furthermore, the curvaturesa, b, and c can be expressed in terms of the first state alone due to the two states having curvatures of equal magnitude. The constraint on the stroke, less than wmax, can therefore be reduced to (9):(9)Actuation stroke constraint:18(a1+b1)(Lx2+Ly2)<wmax⁡.


### 2.4.3. Force of Mechanical Actuation

Similarly, the mechanical force required for actuation is constrained. The force for actuation is considered separately for the transition from state 1 to state 2,F1, and the reverse from state 2 to state 1, F2. These values are not necessarily equal for nonsquare laminates which can result in two shapes with greatly different out-of-plane displacement. The force constraint ensures that both of these values are below the available maximum force, Fmax:(10)Actuation force constraint:F1<Fmax⁡,F2<Fmax⁡.


## 2.4.1. Bistability

The asymmetry of laminates does not guarantee bistability. When a laminate geometry has an edge length to thickness ratio below a critical value there exists only a single saddle solution [10]. This critical ratio is a function of the stacking sequence, temperature, material properties, and aspect ratio. Furthermore, the critical value is increased by the presence of the piezoelectric layers due to their added stiffness. We introduce a constraint to ensure bistability. Noting that the x- and y-curvatures of the monostable saddle shape have equal magnitude in opposite out-of-plane directions, the sum of these curvatures is zero. Beyond the bifurcation point, see Figure 1, the sum of the curvatures is either positive or negative depending on the associated cylindrical state. We therefore impose the constraint a+b>0, where a and b are the x- and y-curvature coefficients, respectively, such that only bistable composites are considered feasible:(6)Bistability constraint:a+b>0.Close to the bifurcation point the solution is highly sensitive to manufacturing imperfections and uncertainties in material data [11]. This issue can easily be addressed by adding a small constant to the constraint, that is, a+b>X, where X is a “distance” away from the bifurcation point. Our previous studies have shown that this constraint could also be used as a minimum out-of-plane deflection constraint [8].


## 2.4.2. Stroke of Mechanical Actuation

The stroke of actuation refers to the distance travelled by the mechanical actuator in thez-direction, Figure 3 and a practical limit is imposed. This is modeled by fixing the displacement at the centre with the maximum deflection being at the corners of the laminate. Noting that the out-of-plane displacement, w of the curved laminate shape is approximated by quadratic polynomial of (7):(7)w(x,y)=12(ax2+by2+cxy)
the corner deflection between states, wdef, at x=Lx/2, y=Ly/2, can be expressed as (8) where the subscript denotes the associated stable state.(8)wdef=w1-w2=18(a1Lx2+b1Ly2+c1LxLy)-18(a2Lx2+b2Ly2+c2LxLy).Due to the orthogonality of the laminate the deflection between states is the same at all four corners. Furthermore, the curvaturesa, b, and c can be expressed in terms of the first state alone due to the two states having curvatures of equal magnitude. The constraint on the stroke, less than wmax, can therefore be reduced to (9):(9)Actuation stroke constraint:18(a1+b1)(Lx2+Ly2)<wmax⁡.


## 2.4.3. Force of Mechanical Actuation

Similarly, the mechanical force required for actuation is constrained. The force for actuation is considered separately for the transition from state 1 to state 2,F1, and the reverse from state 2 to state 1, F2. These values are not necessarily equal for nonsquare laminates which can result in two shapes with greatly different out-of-plane displacement. The force constraint ensures that both of these values are below the available maximum force, Fmax:(10)Actuation force constraint:F1<Fmax⁡,F2<Fmax⁡.


## 3. Analytical Model

### 3.1. Unloaded Laminate Shapes

The analytical model used to calculate the unloaded shapes of asymmetric laminates of arbitrary layup was introduced by Dano and Hyer [5] and is a nonlinear extension to classical laminated plate theory. The co-ordinate system used is that defined in Figure 3, where the geometric centre of the laminate sits at the origin and the ply orientations are measured from the x-axis. The out-of-plane displacement in the z-direction, w is assumed to be of the form of (7). The midplane strains, including geometrical nonlinearity according to the von Karman hypothesis, are defined as(11)εx0=∂u0∂x+12(∂w∂x)2,εy0=∂v0∂y+12(∂w∂y)2,εxy0=∂u0∂y+∂v0∂x+12∂w∂x∂w∂y,
where u0 and v0 are the in-plane displacements in the x- and y-directions, respectively. The midplane strains are approximated by the third order polynomials. Dano and Hyer [5] found that terms with powers of x and y that sum to an odd number are always zero. Therefore the form of the midplane strains can be reduced to the polynomials of (12):(12)εx0=e1+e2x2+e3xy+e4y2,εy0=e5+e6x2+e7xy+e8y2.Using (7), (11), and (12) and introducing the additional shape coefficients e9–11 resulting from integration of the midplane strains, expressions for the in-plane displacements u0 and v0 can be determined:(13)u0(x,y)=e1x+e9y+12(e3-12ac)x2y+(e4-c28)xy2+13(e2-12a2)x3+13e11y3,v0(x,y)=e9x+e5y+12(e7-12bc)xy2+(e6-c28)x2y+13(e8-12b2)y3+13e10x3.The total strain energy of the laminate,W, can then be expressed as the integral of strain energy density over the volume of the laminate.(14)W=∫-Lx/2Lx/2∫-Ly/2Ly/2∫-H/2H/212cijklεijεkl-α̂ijεijΔTdxdydz,
where cijkl is elastic constant, α̂ij is constant relating to the thermal expansion coefficients, Lx and Ly are the planform side lengths of the laminate, H is the total laminate thickness, ΔT is the temperature change from the cure temperature, and εij’s and εkl’s are the total strains defined as:(15)εx=εx0-za,εy=εy0-zb,εxy=εxy0-zc,
where z is the out-of-plane distance from the laminate midplane and a, b, and c are out-of-plane shape coefficients defined by (7). Expansion of (14) results in an expression for the total energy which is a function of the material and geometric properties, the temperature change from the laminate cure and the set of shape coefficients a, b, c, e1…e11. For equilibrium, the minimum energy states require:(16)fi=∂W∂ki=0,i=1,…,14,
where the ki’s are the set of shape coefficients. For the laminate stacking sequence of Figure 2, with square laminate edge lengths, Betts et al. [9] derived an analytical solution to (16) with unloaded orthogonal piezoelectric layers. This provides a fast and reliable method of calculating all stable shapes of asymmetric square laminates, necessary for optimization. In this work we use this analytical solution for a square laminate of edge lengths Lx as the initial solution for rectangular laminate shapes and the Newton-Raphson method is employed to solve the system of (16).


### 3.2. Actuation

In the arrangement shown in Figure3, as used by Arrieta et al. [3], the out-of-plane force F is applied at the centre of the laminate with the four corners held. We model this by fixing the origin and the maximum deflection is at the corners where one quarter of the force is applied at each corner, Figure 5.Figure 5
Force-deflection behavior.The out-of-plane deflection,δw is evaluated using (7) where wa denotes unloaded out-of-plane deflection and wb the deformed out-of-plane deflection:(17)δw=wa-wb=12(δaLx2+δbLy2+δcLxLy).The work done in deflecting the laminate,WF, is hence(18)WF=-4⋅F4⋅12(δaLx2+δbLy2+δcLxLy).The total energy,Wtot, (19) is minimized, where W represents the strain energy of (14):(19)Wtot=W-WF,(20)δWtot=∑i=114∂∂kiW+F2(δaLx2+δbLy2+δcLxLy)=0.Equation (20) is solved using a simple iterative scheme. The exact unloaded laminate shape is obtained first as described in Section 3.1 and the force is incrementally increased until the state transition. Figure 6 exhibits an example of determining a snap-through force where the applied force and the corresponding curvatures are plotted to show the state transition. Starting with the unloaded laminate at A there is a large positive x-curvature for state 1. Between points A and B the force is increased. From point B the curvature is seen to change dramatically with x-curvature becoming very small and the y-curvature becoming large and negative at point C. This point marks the transition from state 1 to state 2. As the load is then incrementally removed the laminate curvatures tend towards point D which represents the unloaded second state. This process can be repeated in the opposite out-of-plane direction to reset the laminate and calculate the second snap force. These two values define the force constraints for optimization.Changing curvature and snap-through behavior of a[0n/90n]T laminate with applied force.
(a)(b)


## 3.1. Unloaded Laminate Shapes

The analytical model used to calculate the unloaded shapes of asymmetric laminates of arbitrary layup was introduced by Dano and Hyer [5] and is a nonlinear extension to classical laminated plate theory. The co-ordinate system used is that defined in Figure 3, where the geometric centre of the laminate sits at the origin and the ply orientations are measured from the x-axis. The out-of-plane displacement in the z-direction, w is assumed to be of the form of (7). The midplane strains, including geometrical nonlinearity according to the von Karman hypothesis, are defined as(11)εx0=∂u0∂x+12(∂w∂x)2,εy0=∂v0∂y+12(∂w∂y)2,εxy0=∂u0∂y+∂v0∂x+12∂w∂x∂w∂y,
where u0 and v0 are the in-plane displacements in the x- and y-directions, respectively. The midplane strains are approximated by the third order polynomials. Dano and Hyer [5] found that terms with powers of x and y that sum to an odd number are always zero. Therefore the form of the midplane strains can be reduced to the polynomials of (12):(12)εx0=e1+e2x2+e3xy+e4y2,εy0=e5+e6x2+e7xy+e8y2.Using (7), (11), and (12) and introducing the additional shape coefficients e9–11 resulting from integration of the midplane strains, expressions for the in-plane displacements u0 and v0 can be determined:(13)u0(x,y)=e1x+e9y+12(e3-12ac)x2y+(e4-c28)xy2+13(e2-12a2)x3+13e11y3,v0(x,y)=e9x+e5y+12(e7-12bc)xy2+(e6-c28)x2y+13(e8-12b2)y3+13e10x3.The total strain energy of the laminate,W, can then be expressed as the integral of strain energy density over the volume of the laminate.(14)W=∫-Lx/2Lx/2∫-Ly/2Ly/2∫-H/2H/212cijklεijεkl-α̂ijεijΔTdxdydz,
where cijkl is elastic constant, α̂ij is constant relating to the thermal expansion coefficients, Lx and Ly are the planform side lengths of the laminate, H is the total laminate thickness, ΔT is the temperature change from the cure temperature, and εij’s and εkl’s are the total strains defined as:(15)εx=εx0-za,εy=εy0-zb,εxy=εxy0-zc,
where z is the out-of-plane distance from the laminate midplane and a, b, and c are out-of-plane shape coefficients defined by (7). Expansion of (14) results in an expression for the total energy which is a function of the material and geometric properties, the temperature change from the laminate cure and the set of shape coefficients a, b, c, e1…e11. For equilibrium, the minimum energy states require:(16)fi=∂W∂ki=0,i=1,…,14,
where the ki’s are the set of shape coefficients. For the laminate stacking sequence of Figure 2, with square laminate edge lengths, Betts et al. [9] derived an analytical solution to (16) with unloaded orthogonal piezoelectric layers. This provides a fast and reliable method of calculating all stable shapes of asymmetric square laminates, necessary for optimization. In this work we use this analytical solution for a square laminate of edge lengths Lx as the initial solution for rectangular laminate shapes and the Newton-Raphson method is employed to solve the system of (16).


## 3.2. Actuation

In the arrangement shown in Figure3, as used by Arrieta et al. [3], the out-of-plane force F is applied at the centre of the laminate with the four corners held. We model this by fixing the origin and the maximum deflection is at the corners where one quarter of the force is applied at each corner, Figure 5.Figure 5
Force-deflection behavior.The out-of-plane deflection,δw is evaluated using (7) where wa denotes unloaded out-of-plane deflection and wb the deformed out-of-plane deflection:(17)δw=wa-wb=12(δaLx2+δbLy2+δcLxLy).The work done in deflecting the laminate,WF, is hence(18)WF=-4⋅F4⋅12(δaLx2+δbLy2+δcLxLy).The total energy,Wtot, (19) is minimized, where W represents the strain energy of (14):(19)Wtot=W-WF,(20)δWtot=∑i=114∂∂kiW+F2(δaLx2+δbLy2+δcLxLy)=0.Equation (20) is solved using a simple iterative scheme. The exact unloaded laminate shape is obtained first as described in Section 3.1 and the force is incrementally increased until the state transition. Figure 6 exhibits an example of determining a snap-through force where the applied force and the corresponding curvatures are plotted to show the state transition. Starting with the unloaded laminate at A there is a large positive x-curvature for state 1. Between points A and B the force is increased. From point B the curvature is seen to change dramatically with x-curvature becoming very small and the y-curvature becoming large and negative at point C. This point marks the transition from state 1 to state 2. As the load is then incrementally removed the laminate curvatures tend towards point D which represents the unloaded second state. This process can be repeated in the opposite out-of-plane direction to reset the laminate and calculate the second snap force. These two values define the force constraints for optimization.Changing curvature and snap-through behavior of a[0n/90n]T laminate with applied force.
(a)(b)


## 4. Results

This section presents the results of a design study. The optimization problem is solved using Matlab’s sequential quadratic programming, fmincon. All examples use M21/T800 material properties [12] for the laminate and M8557-P1 MFC properties for the piezoelectric layers [13].


### 4.1. Influence of Aspect Ratio for a Single Laminate

Prior to optimization, we first investigate the influence of aspect ratio on the strain energy components of (5) using [0/0/90/90]T laminates of aspect ratio 1 and 5 (Figure 7). A breakdown of the objective function values is shown in Table 1. For the laminate with an aspect ratio of 1 (Figure 7(a)), two equal and opposite shapes are observed and each state contributes 50% of the total objective function. Of each half, the vast majority is contributed by the longitudinal direction of the piezoelectric aligned with the major curvature, Uxt for state 1 and Uyb for state 2. Interestingly the same pattern is observed for the aspect ratio of 5 where the two stable states have very different displacements, Figure 7(b). Despite the appearance, states 1 and 2 exhibit the same magnitude of curvature in opposite out-of-plane directions, thus the vast majority of the electrical energy is provided by the longitudinal direction of the piezoelectric layer aligned with the major curvature. Our detailed investigation revealed that this is consistently observed for all cross-ply laminates.Table 1
Breakdown of electrical energy for two[0/0/90/90]T laminates.
Component% Total energy, AR = 1% Total energy, AR = 5State 1Uxt46.2647.31Uyt0.220.03Uxb2.632.53Uyb0.890.13Total50.0050.00State 2Uxt0.890.13Uyt2.632.53Uxb0.220.03Uyb46.2647.31Total50.0050.00Two stable states of a[0/0/90/90]T laminate with aspect ratio of (a) 1 and (b) 5.
(a)(b)


### 4.2. Optimization of Cross-Ply Laminates

This section presents optimization of a[0/0/90/90]T laminate to determine ply thickness, t, and the laminate aspect ratio, Lx/Ly. The piezoelectric layers are oriented at 0° and 90° and have a fixed thickness of 0.2 mm. The surface area of the laminate and the piezoelectric layers is fixed to 0.04 m2. The full design space for this problem constrained by bistability, a maximum deflection of 25 mm and a maximum reversible snap force of 5 N is shown in Figure 8. The colored contours represent the objective function value bounded by each of the constraints.Figure 8
Full design space for optimization of a[0/0/90/90]T laminate. Colored contours represent the objective function, bounded by bistability, deflection and force constraints. The white circle marks a local optimum, the black circle the global optimum, and the white squares the local solutions constrained by discrete ply thicknesses.The outer limits of the design space (the white “monostable” area) are due to the bistability constraint and are observed to be at a maximum half-laminate thickness of 0.866 mm (with an aspect ratio 1), and a maximum aspect ratio of 6.85 (at a half-laminate thickness of 0.196 mm), Figure8. The left hand boundary of the monostable region defines the minimum laminate thickness to achieve bistability. As the thickness of the piezoelectric is fixed, the thickness ratio of the laminates to piezoelectrics becomes too small beyond this minimum limit and the stiffness contribution of the piezoelectrics to the piezo-composite laminate becomes significant, leading to the loss of bistability. The right hand boundary represents the bistability bifurcation limit due to low edge to thickness ratio. As shown in Figure 1, there is only one saddle equilibrium state for low edge to thickness laminates. For the optimization study presented here, the laminate area is fixed which limits the maximum achievable edge length. As the laminate thickness increases on the x-axis in Figure 8, the edge to thickness ratio decreases and reaches the right hand boundary beyond which the laminates can only have one equilibrium state.There are three grey areas marking the infeasible regions. The lightest region indicates laminates which violates the actuation stroke constraint, that is, they require a deflection greater than the allowable maximum. This constraint violation, however, does not influence the optimum designs as the objective function values of these solutions are relatively low.The other two constraints, snap-through forcesF1 and F2, required to induce state transition, however, are critical and determine the optimum solutions. These are two darker grey regions in Figure 8 where the objective function values are maximum. The laminates in these regions require a greater force to induce the state change than the available force. For the case of a [0/0/90/90]T laminate, the snap force from state 1 to 2 is completely enclosed by the snap force from state 2 to 1. This is a result of the aspect ratio where the displacements between two states can vary significantly for high aspect ratio (e.g., Figure 7). The difference between the snap force requirements is more prominent for low thickness laminates in Figure 8 and indeed, these regions meet at aspect ratio 1 indicating that the force requirement from states 1 to 2 and 2 to 1 are identical for square laminates. The lower limit of the constraint boundary indicates that the snap force requirement is reduced as the thickness of the laminates decreases. As the thickness increases, the increased stiffness of the laminates reduces the curvature thus requiring less snap force and displacement. This gives rise to the upper constraint boundary.Two optimal solutions are obtained. The first, indicated as a white circle in Figure8 is a local optimum with an objective function value of 0.012 J. This laminate has an aspect ratio of 1 and half-laminate thickness of 0.245 mm. The second, the black circle, has an aspect ratio of 1 and half-laminate thickness of 0.731 mm, giving an objective function value of 0.0295 J. The second solution is the global optimum.


### 4.3. Discrete Ply Thicknesses

The problem so far assumed that the ply thickness variation is continuous. However in practice, the number of plies and, thus, the laminate thickness are discrete. For a typical single ply thickness of 0.125 mm, we can reconsider the design space of Figure8 as a [0n/90n]T laminate with n ranging from 1 to 6, that is, 2- to 12-ply laminates. With this added constraint the discrete optimum solution for each number of plies are marked by white squares in Figure 8 and summarized in Table 2. The solutions tend to be for laminates with aspect ratio 1 unless the square laminates violate a constraint. The active constraint in this case was the mechanical force limit. The global solution was found to be a square laminate with n=6, whose objective function is an order of magnitude greater than of the n=1 solution.Table 2
Optimum solution for discrete laminate thicknesses.
No. of plies,nAspect ratioObjective function value (J)110.0027321.870.012132.500.020442.390.025051.830.0281610.0285


## 4.1. Influence of Aspect Ratio for a Single Laminate

Prior to optimization, we first investigate the influence of aspect ratio on the strain energy components of (5) using [0/0/90/90]T laminates of aspect ratio 1 and 5 (Figure 7). A breakdown of the objective function values is shown in Table 1. For the laminate with an aspect ratio of 1 (Figure 7(a)), two equal and opposite shapes are observed and each state contributes 50% of the total objective function. Of each half, the vast majority is contributed by the longitudinal direction of the piezoelectric aligned with the major curvature, Uxt for state 1 and Uyb for state 2. Interestingly the same pattern is observed for the aspect ratio of 5 where the two stable states have very different displacements, Figure 7(b). Despite the appearance, states 1 and 2 exhibit the same magnitude of curvature in opposite out-of-plane directions, thus the vast majority of the electrical energy is provided by the longitudinal direction of the piezoelectric layer aligned with the major curvature. Our detailed investigation revealed that this is consistently observed for all cross-ply laminates.Table 1
Breakdown of electrical energy for two[0/0/90/90]T laminates.
Component% Total energy, AR = 1% Total energy, AR = 5State 1Uxt46.2647.31Uyt0.220.03Uxb2.632.53Uyb0.890.13Total50.0050.00State 2Uxt0.890.13Uyt2.632.53Uxb0.220.03Uyb46.2647.31Total50.0050.00Two stable states of a[0/0/90/90]T laminate with aspect ratio of (a) 1 and (b) 5.
(a)(b)


## 4.2. Optimization of Cross-Ply Laminates

This section presents optimization of a[0/0/90/90]T laminate to determine ply thickness, t, and the laminate aspect ratio, Lx/Ly. The piezoelectric layers are oriented at 0° and 90° and have a fixed thickness of 0.2 mm. The surface area of the laminate and the piezoelectric layers is fixed to 0.04 m2. The full design space for this problem constrained by bistability, a maximum deflection of 25 mm and a maximum reversible snap force of 5 N is shown in Figure 8. The colored contours represent the objective function value bounded by each of the constraints.Figure 8
Full design space for optimization of a[0/0/90/90]T laminate. Colored contours represent the objective function, bounded by bistability, deflection and force constraints. The white circle marks a local optimum, the black circle the global optimum, and the white squares the local solutions constrained by discrete ply thicknesses.The outer limits of the design space (the white “monostable” area) are due to the bistability constraint and are observed to be at a maximum half-laminate thickness of 0.866 mm (with an aspect ratio 1), and a maximum aspect ratio of 6.85 (at a half-laminate thickness of 0.196 mm), Figure8. The left hand boundary of the monostable region defines the minimum laminate thickness to achieve bistability. As the thickness of the piezoelectric is fixed, the thickness ratio of the laminates to piezoelectrics becomes too small beyond this minimum limit and the stiffness contribution of the piezoelectrics to the piezo-composite laminate becomes significant, leading to the loss of bistability. The right hand boundary represents the bistability bifurcation limit due to low edge to thickness ratio. As shown in Figure 1, there is only one saddle equilibrium state for low edge to thickness laminates. For the optimization study presented here, the laminate area is fixed which limits the maximum achievable edge length. As the laminate thickness increases on the x-axis in Figure 8, the edge to thickness ratio decreases and reaches the right hand boundary beyond which the laminates can only have one equilibrium state.There are three grey areas marking the infeasible regions. The lightest region indicates laminates which violates the actuation stroke constraint, that is, they require a deflection greater than the allowable maximum. This constraint violation, however, does not influence the optimum designs as the objective function values of these solutions are relatively low.The other two constraints, snap-through forcesF1 and F2, required to induce state transition, however, are critical and determine the optimum solutions. These are two darker grey regions in Figure 8 where the objective function values are maximum. The laminates in these regions require a greater force to induce the state change than the available force. For the case of a [0/0/90/90]T laminate, the snap force from state 1 to 2 is completely enclosed by the snap force from state 2 to 1. This is a result of the aspect ratio where the displacements between two states can vary significantly for high aspect ratio (e.g., Figure 7). The difference between the snap force requirements is more prominent for low thickness laminates in Figure 8 and indeed, these regions meet at aspect ratio 1 indicating that the force requirement from states 1 to 2 and 2 to 1 are identical for square laminates. The lower limit of the constraint boundary indicates that the snap force requirement is reduced as the thickness of the laminates decreases. As the thickness increases, the increased stiffness of the laminates reduces the curvature thus requiring less snap force and displacement. This gives rise to the upper constraint boundary.Two optimal solutions are obtained. The first, indicated as a white circle in Figure8 is a local optimum with an objective function value of 0.012 J. This laminate has an aspect ratio of 1 and half-laminate thickness of 0.245 mm. The second, the black circle, has an aspect ratio of 1 and half-laminate thickness of 0.731 mm, giving an objective function value of 0.0295 J. The second solution is the global optimum.


## 4.3. Discrete Ply Thicknesses

The problem so far assumed that the ply thickness variation is continuous. However in practice, the number of plies and, thus, the laminate thickness are discrete. For a typical single ply thickness of 0.125 mm, we can reconsider the design space of Figure8 as a [0n/90n]T laminate with n ranging from 1 to 6, that is, 2- to 12-ply laminates. With this added constraint the discrete optimum solution for each number of plies are marked by white squares in Figure 8 and summarized in Table 2. The solutions tend to be for laminates with aspect ratio 1 unless the square laminates violate a constraint. The active constraint in this case was the mechanical force limit. The global solution was found to be a square laminate with n=6, whose objective function is an order of magnitude greater than of the n=1 solution.Table 2
Optimum solution for discrete laminate thicknesses.
No. of plies,nAspect ratioObjective function value (J)110.0027321.870.012132.500.020442.390.025051.830.0281610.0285


## 5. Conclusions

This paper presented a modeling and optimization formulation of bistable piezo-composite laminates for energy harvesting applications. The method allows design parameters such as harvester aspect ratio and device thickness to be optimized based on realistic design constraints such as snap-through forces, displacements, and the influence of the stiffness of the piezoelectric elements on bistability of the laminate. This preliminary investigation of optimum cross-ply laminates has found that despite the high deflection of high aspect ratio laminates, the square laminates are the optimum with highest electrical energy. The mechanical force required to induce state change was found to be the critical factor in determining the optimum solution whilst the deflection constraint was observed to be inactive. The discrete ply thickness and its influence on the solutions were also examined. Based on this preliminary study, continuing research is extending the optimization methodology to include additional design constraints such as the maximum stress in the piezoelectric material to prevent depolarization. Further design parameters are also considered such as choice of piezoelectric material, piezoelectric area, and configurations and stacking sequence of non-cross-ply laminates.



---

*Source: 621364-2012-04-09.xml*