1 INTRODUCTION

The piezoelectric smart structures have dominated the current state of the art shape and
vibration control technology due to their high reliability (^{Crawley and de Luis, 1987}). Electromechanical coupling present in the
piezoelectric materials adds the ability to respond to the subjected forced environment.
Piezoelectric beams render a wide class of smart structures (^{Benjeddou et al., 1997}). Their numerical analysis plays an important role in the
design of piezoelectric beam based control systems. The analysis of piezoelectric beam
requires efficient and accurate modeling of both mechanical and electric responses (^{Sulbhewar and Raveendranath, 2014a}; ^{Zhou et al., 2005}).

As most of the smart beams are thin (^{Marinkovic and
Marinkovic, 2012}), Euler-Bernoulli beam theory has been favored by many researchers
for their formulations. Analytical models based on the Euler-Bernoulli theory were presented
by ^{Crawley and de Luis (1987)}, ^{Crawley and Anderson (1990)} and ^{Abramovich
and Pletner (1997)} for static and dynamic analyses of piezoelectric beams. ^{Crawley and de Luis (1987)} carried out the parametric
studies to show the effectiveness of piezoelectric actuators in transmitting the strain to
the substrate. ^{Crawley and Anderson (1990)} extended
this analysis to establish the range over which the simpler analytical solutions are valid.
^{Abramovich and Pletner (1997)} proposed closed form
solutions for induced curvature and axial strain of an adaptive sandwich structure.

Also, the classical laminate theory based piezoelectric plate finite elements (^{Hwang and Park, 1993}; ^{Lam
et al., 1997}; ^{Wang et al., 1997}) and beam
elements (^{Balamurugan and Narayanan, 2002}; ^{Bendary et al., 2010}; ^{Bruant et al., 2001}; ^{Carpenter, 1997}; ^{Gaudenzi et al., 2000}; ^{Hanagud et al., 1992}; ^{Kumar and Narayanan,
2008}; ^{Manjunath and Bandyopadhyay, 2004};
^{Robbins and Reddy, 1991}; ^{Sadilek and Zemcik, 2010}; ^{Stavroulakis et
al., 2005}; ^{Zemcik and Sadilek, 2007}) are
available in the literature for the analysis of extension mode piezoelectric beams.

A four-noded classical plate bending element with an elemental electric potential degree of
freedom was used by ^{Hwang and Park (1993)} for the
vibration control analysis of composite plates with piezoelectric sensors and actuators,
while with nodal electric potential degrees of freedom for static shape control study by
^{Wang et al. (1997)}. ^{Lam et al. (1997)} developed a four-noded classical plate finite element with
negative velocity feedback control to study the active dynamic response of integrated
structures.

A very early finite element based on a classical beam theory proposed by ^{Robbins and Reddy (1991)}, is without electrical nodal
degrees of freedom and considered the strain induced by piezoelectric material as an applied
force. The active vibration control of beams using piezoelectric material was studied using
Euler-Bernoulli piezoelectric beam finite elements by ^{Balamurugan and Narayanan (2002)}, ^{Bruant et al.
(2001)}, ^{Gaudenzi et al. (2000)}, ^{Hanagud et al. (1992)}, ^{Kumar and Narayanan (2008)}, ^{Manjunath and
Bandyopadhyay (2004)} and ^{Stavroulakis et al.
(2005)}. ^{Carpenter (1997)} used energy methods
to derive Euler-Bernoulli beam finite element for analyzing elastic beams with piezoelectric
materials. Euler-Bernoulli beam element with assumed bilinear electric potential was used by
^{Zemcik and Sadilek (2007)} and ^{Sadilek and Zemcik (2010)} for modal and frequency response analysis of
piezoelectric smart beams, respectively. ^{Bendary et al.
(2010)} carried out the static and dynamic analyses of composite beams with
piezoelectric materials using a classical laminate theory based beam element with layerwise
linear through-thickness distribution of nodal electric potential.

Generally, the Euler-Bernoulli theory based beam elements give good convergence of finite
element results (^{Reddy, 1997}). However, when the beam
cross-section consists of asymmetrical material distribution, these elements suffer from
*material-locking* due to presence of *extension-bending*
coupling (^{Raveendranath et al., 2000}). As most of the
piezoelectric beam structures contain a piezoelectric material layer asymmetrically bonded
to the host structure, it activates the extension-bending coupling and hence deteriorates
the finite element convergence. One of the ways to get rid of material-locking is the use of
a higher-order polynomial interpolation for the axial displacement along the length of the
beam (^{Raveendranath et al., 2000}; ^{Sulbhewar and Raveendranath, 2014b}). A more efficient way
to remove ill-effects of material-locking is the use of coupled polynomial field
interpolations which incorporate the extension-bending coupling in an appropriate manner,
without increasing nodal degrees of freedom (^{Raveendranath
et al., 2000}; ^{Sulbhewar and Raveendranath,
2014c}; ^{Sulbhewar and Raveendranath, 2014d}).
For Euler-Bernoulli piezoelectric beam finite element, electromechanical coupling also has
to be accommodated properly along with the extension-bending coupling in a coupled
polynomial field based formulation.

Here in this study, a novel interpolation scheme, based on coupled polynomials derived from governing equilibrium equations, is proposed to eliminate material-locking in the Euler-Bernoulli piezoelectric beam finite element. The variational formulation is used to derive governing equilibrium equations which are used to establish the relationship between field variables. The polynomial derived here for the axial displacement field for the beam contains an additional coupled term along with conventional linear terms, having contributions from an assumed cubic polynomial for transverse displacement and assumed linear polynomials for layerwise electric potential. This coupled term facilitates the use of higher order polynomial for axial displacement without introducing any additional generalized degrees of freedom in the formulation. A set of coupled shape functions is obtained from these coupled field polynomials which accommodates extension-bending and piezoelectric couplings in a variationally consistent manner. The accuracy of the present formulation is shown by comparing results with those from ANSYS 2D simulation and conventional formulations. Convergence studies are carried out to prove the efficiency of the present Euler-Bernoulli piezoelectric beam formulation with coupled shape functions over the conventional formulations with assumed independent shape functions.

2 THEORETICAL FORMULATION

The formulation is based on the equivalent single layer (ESL) Euler-Bernoulli theory for
mechanical field and a layerwise linear model for electric potential (*φ*).
The multilayered piezoelectric smart beam is considered, as shown in Figure 1. The layer(s) can be host layers of conventional/composite
material or bonded/embedded layer(s) of piezoelectric material. The beam layers are assumed
to be made up of isotropic or specially orthotropic materials, with perfect bonding among
them. The top and bottom faces of the piezoelectric layers are fully covered with
electrodes. Mechanical and electrical quantities are assumed to be small enough to apply
linear theories of elasticity and piezoelectricity.

2.1 Mechanical Displacements and Strain

The Euler-Bernoulli theory is used for which the field expressions are given as (^{Bendary et al., 2010}):

where *u*(*x,z*) and
*w*(*x,z*) are the displacements in the longitudinal and
transverse directions of the beam, respectively. *u*
_{0}(*x*) and *w*
_{0}(*x*) are the centroidal axial and transverse displacements,
respectively. ()' denotes derivative with respect to *x*. Dimensions
*L,b,h* denote the length, width and the total thickness of the beam,
respectively.

Axial strain field is derived using usual strain-displacement relation as:

2.2 Electric Potential and Electric Field

The through-thickness profile of the electric potential *φ*
(*x,z*) in a piezoelectric layer is taken as linear (^{Bendary et al., 2010}). As shown in Figure 1, the two dimensional electric potential takes the values of ϕ*
_{i+1}
*(

*x*) and ϕ

*(*

_{i+1}*x*) at the top and bottom faces of the

*i*

*th*piezoelectric layer, respectively. The electric field in the transverse direction (

*E*

*) can be derived from electric potential as (*

_{z}^{Bendary et al., 2010}):

where =
*ϕ*
*
_{i+1}
* -

*ϕ*

*, is the difference of potentials at the top and bottom faces of the*

_{i+1}*i*

*th*piezoelectric layer with thickness

*h*

*.*

_{i}3 REDUCED CONSTITUTIVE RELATIONS

The smart beam consisting layers of conventional/composite/piezoelectric material with
isotropic or specially orthotropic properties is considered here. It has axes of material
symmetry parallel to beam axes, in which electric field is applied in the transverse
direction. For extension mode beam, the transversely poled piezoelectric material is
subjected to the transverse electric field. The elastic, piezoelectric and dielectric
constants are denoted by *C*
*
_{ij}
*,

*e*

*(*

_{kj}*i*,

*j*= 1.....) and ∈

*(*

_{k}*k*= 1,2,3), respectively. The constitutive relations for such a system are given as (

^{Sulbhewar and Raveendranath, 2014c}):

where (*i *=1....number of piezoelectric layers), (*k*
=1.....total number of layers in beam). The constants , and denote elastic
(*N*/*m*
^{2)}, piezoelectric (*C*/*m*
^{2)} and dielectric (*F*/*m*) properties reduced for
the beam geometry (^{Sulbhewar and Raveendranath,
2014c}).

4 VARIATIONAL FORMULATION

Hamilton's principle is used to formulate the piezoelectric smart beam. It is expressed as
(^{Sulbhewar and Raveendranath, 2014 c}):

where, *K*=kinetic energy, *H*=electric enthalpy density
function for piezoelectric material and mechanical strain energy for the linear elastic
material and *W*=external work done.

4.1 Electromechanical and Strain Energy Variations

For the *j*
*th* conventional/composite material layer, the mechanical strain energy
variation is given as:

The electromechanical strain energy variation of the *i*
*th* piezoelectric layer is given as:

Substituting values of axial strain (*ε*
*
_{x}
*), electric field ( ) form equations (3), (4) and
using them along with constitutive relations given by equation (5) in expressions (7) and (8); the total variation
on the potential energy of the smart beam is given as:

where (*i*=1....number of piezoelectric layers),
(*k*=1.....total number of layers in beam) and

4.2 Variation of Kinetic Energy

Total kinetic energy of the beam is given as (^{Sulbhewar
and Raveendranath, 2014c}):

where *ρ*
*
_{k}
* is volumic mass density of

*k*

*th*layer in

*kgm*

^{-3}and (

*k*=1...total number of layers in the beam). Substituting values of

*u*and

*w*from equations (1) and (2) and applying variation, to derive at:

4.3 Variation of Work of External Forces

Total virtual work of the structure can be defined as the product of virtual
displacements with forces for the mechanical work and the product of the virtual electric
potential with the charges for the electrical work. The variation of total work done by
external mechanical and electrical loading is given by (^{Sulbhewar and Raveendranath, 2014c}):

in which *ƒ*
*V*
*, ƒ*
*S*
*, ƒ*
*C* are volume, surface and point forces, respectively. *q*
_{0} and *S*
*
_{φ}
* are the charge density and area on which charge is applied.

5 DERIVATION OF COUPLED FIELD RELATIONS

The relationship between field variables is established here using static governing
equations. For static conditions without any external loading, the variational principle
given in equation (6) reduces to (^{Sulbhewar and Raveendranath, 2014c}):

Applying variation to the basic variables in equation (9), the static governing equations are obtained as:

From equation (14), we can write the
relationship of axial displacement (*u*
_{0}) with transverse displacement (*w*
_{0}) and electric potential (*φ*) as:

where and . The constants which relate the field variables are the functions of geometric and material properties of the beam. This relation is used in the next section to derive coupled field polynomials.

6 FINITE ELEMENT FORMULATION

Using the variational formulation described above, a finite element model has been
developed here. The model consists of two mechanical variables (*u*
_{0} and *w*
_{0}) and layerwise electrical variables ( ) where
(*i*=1.....number of piezoelectric layers in the beam). In terms of natural
coordinate *ξ*, a cubic polynomial for transverse displacement
(*w*
_{0}) and linear polynomials for are assumed as given in equations 17 (a) and 17 (b), respectively. The transformation between
*ξ* and global coordinate *x* along the length of the beam
is given as *ξ* = 2(*x*-*x*
_{1})/(*x*
_{2}-*x*
_{1}) - 1 with (*x*
_{2}-*x*
_{1})=1, being the length of the beam element.

Using these polynomials for *w*
_{0 }and in
equation (16) and integrating with respect
to *ξ*, we get the quadratic polynomial for axial displacement
*u*
_{0} as:

It is noteworthy that equation (18) contains
a quadratic term in addition to the conventional linear terms. The formulation which use
linear interpolation for *u*
_{0} (^{Bendary et al., 2010}) is known to
exhibit material-locking when applied to beams with asymmetric cross-section (^{Sulbhewar and Raveendranath, 2014b}). A higher-order
interpolation for *u*
_{0} as given in equation (18), is expected to relieve material-locking and improve
the rate of convergence (^{Sulbhewar and Raveendranath,
2014b}). It takes care of extension-bending and electromechanical couplings in a
variationally consistent manner with the help of the coupled quadratic term which has
contributions from transverse displacement and layerwise electric potential. The
takes care of the
bending effects, while
takes care of electromechanical effects on the axial displacement. As seen from equations (17) and (18), while employing a quadratic polynomial for interpolation of axial
displacement, the generalized degrees of freedom are efficiently maintained the same as of
the conventional formulation of ^{Bendary et al. (2010)}.

Using equations (17) and (18), the coupled shape functions in equation (19) are derived by usual method.

The expressions for these shape functions in the natural coordinate system are given as:

Now the variation on basic mechanical and electrical variables can be transferred to nodal degrees of freedom. Substituting equation (19) in equations (9), (11), (12) and using them in equation (6), the following discretized form of the model is obtained:

where *M* is the mass matrix, *K*
*
_{uu}
*,

*K*

*,*

_{uφ}*K*

*,*

_{φu}*K*

*are global stiffness sub-matrices.*

_{φφ}*U*, Φ are the global nodal mechanical displacement and electric potential degrees of freedom vectors, respectively.

*F*and

*Q*are global nodal mechanical and electrical force vectors, respectively. Now the general formulation has been converted to matrix equation which can be solved according to electrical conditions (open/closed circuit), configuration (actuator/sensor) and type of analysis (static/dynamic).

7 NUMERICAL EXAMPLES AND DISCUSSIONS

The developed formulation is tested here for the accuracy and efficiency to predict the
structural behaviour of piezoelectric extension mode smart beams. The formulation is
validated for static (actuator/sensor configuration) and modal analyses. The present
Euler-Bernoulli smart beam formulation which uses coupled polynomial based interpolation
(hereafter designated as EB-CPI) is compared for performance against conventional
Euler-Bernoulli smart beam formulations with assumed independent polynomial based
interpolation as used by ^{Bendary et al. (2010)}
(designated hereafter as EB-IPI).

The particular test problem chosen here is the host layer made up of aluminum or 0/90 graphite-epoxy composite material with surface bonded piezoelectric PVDF material layer as shown in Figure 2. The material and geometric properties of the beam are:

Aluminum (^{Bendary et al., 2010}):
*E*=68.9*Gpa*; *ν*=0.25;
*ρ*=2769*kgm*
^{-3}.

Graphite-epoxy (^{Beheshti-Aval and Lezgy-Nazargah,
2013}): (*E*
*
_{L}
*,

*E*

*,*

_{T}*G*

*,*

_{LT}*G*

*) = (181,10.3,7.17,287)*

_{TT}*Gpa*(

*ν*

_{LT}*, ν*

*) = (0.28,0.33);*

_{TT}*ρ*= 1578

*kgm*

^{-3}.

PVDF (^{Sun and Huang, 2000}): *E* =
2*Gpa*; *ν*=0.29; *e*
_{31} = 0.046*Cm*
^{-2}; ∈_{3} = 0.1062×10^{-9}
*Fm*
^{-1 }ρ = 1800*kgm*
^{-3 }
*h*
*
_{c}
* = 4

*mm*;

*h*

*= 1*

_{p}*mm*; L = 100

*mm*.

The set of reference values is obtained from the numerical simulation of test problems with the ANSYS 2D simulation, for which a mesh of PLANE 183 element for host layer with size 100×9 and a mesh of PLANE 223 element for piezoelectric layer with size 100×2 is used.

7.1 Static Analysis: Sensor Configuration

For sensor configuration, the piezoelectric smart cantilever shown in Figure 2 is subjected to a tip load of 1000
*N*. The structural response of the beam is evaluated by present EB-CPI,
conventional EB-IPI of ^{Bendary et al. (2010)} and
ANSYS 2D simulation. The variations of transverse and axial deflections and potential
developed across the piezoelectric layer along the length of the 0/90/PVDF and Al/PVDF
beams are plotted in Figures 3, 4 and 5, respectively. As seen from
the plots, EB-CPI gives accurate predictions as of EB-IPI of ^{Bendary et al. (2010)} and ANSYS 2D simulation. The through-thickness
variations of axial stresses developed at the roots of the 0/90/PVDF and Al/PVDF beams are
plotted in Figure 6, which validate the use of
present coupled polynomial interpolation for calculation of derived quantities also.

Now the efficiency of the present coupled polynomial based interpolation scheme over the
conventional independent polynomial based interpolation is proved by the convergence
graphs plotted in Figures 7 and 8 for tip deflection and potential developed at the root of the beam,
respectively. As seen from the graphs, EB-IPI of ^{Bendary et
al. (2010)} exhibit material-locking and takes a large number of elements to
converge to the accurate predictions as given by 2D simulation. The present EB-CPI is free
from material-locking effects and shows single element convergence. This improvement in
the performance can be attributed to the coupled term present in the description of the
axial displacement given by equation (18).
The performance of the conventional formulation (^{Bendary et
al., 2010}) may be improved by using an assumed independent quadratic
interpolation for the axial displacement; however it will increase the number of elemental
degrees of freedom and computational efforts.

7.2 Static Analysis: Actuator Configuration

For actuator configuration, the voltage of 10 KV is applied to the top surface of
piezoelectric layer of the smart cantilever shown in Figure
2. The variations of transverse and axial deflections along the length of the
0/90/PVDF and Al/PVDF cantilever beams, obtained by EB-CPI, EB-IPI of ^{Bendary et al. (2010)} and ANSYS 2D simulation are
plotted in Figures 9 and 10, respectively. The through-thickness variations of axial stresses
developed at the roots of the 0/90/PVDF and Al/PVDF beams are plotted in Figure 11. As seen from the deflection and stress plots,
results by the present EB-CPI closely match with the results given by 2D simulation and
EB-IPI of ^{Bendary et al. (2010)}, which prove the
accuracy of the present formulation in actuator configuration.

7.3 Modal Analysis

The present EB-CPI is tested here for the accuracy and efficiency to predict the natural
frequencies of the piezoelectric smart cantilever shown in Figure 2. The converged values of first three natural frequencies of the
0/90/PVDF and Al/PVDF beams evaluated by EB-CPI, EB-IPI of ^{Bendary et al. (2010)} and ANSYS 2D simulation are tabulated in Table 1. As seen from the results, EB-CPI gives the
accurate predictions for the natural frequencies as predicted by 2D simulation and EB-IPI.
This validates the use of consistent mass matrices generated by the present coupled
polynomial field based interpolation scheme.

Now, the efficiency of the present interpolation scheme over conventional one is proved
by the convergence graphs plotted for the first three natural frequencies in Figures 12, 13 and
14. These convergence graphs once again prove the
merit of present EB-CPI over EB-IPI of ^{Bendary et al.
(2010)}. Also, the first three natural bending modes of the 0/90/PVDF
piezoelectric smart cantilever are plotted in Figure
15.

8 CONCLUSION

An efficient Euler-Bernoulli piezoelectric beam finite element formulation has been proposed. The efficiency has been achieved by using the coupled polynomial field interpolation scheme. The relationship between field variables established using governing equations, is used to generate a coupled polynomial for the axial displacement field for the beam. Consequently, the shape functions derived are also coupled which handle extension-bending and piezoelectric coupling in an efficient manner. From the numerical analysis, it was found that:

The present formulation gives accurate predictions of finite element results in both static (sensing/actuation) and modal analyses as given by conventional formulation and 2D finite element simulation, which validate the use of present coupled polynomial field based interpolation.

The present interpolation scheme is free from any locking effect unlike conventional independent polynomial based interpolation which suffers from material-locking. The convergence studies carried out prove the merit of the present coupled polynomial based interpolation over the conventional independent polynomial based interpolations.

The present formulation proves to be the most efficient way to model piezoelectric smart beam with Euler-Bernoulli beam theory.