An Efficient Coupled Polynomial Interpolation Scheme to Eliminate Material-locking in the Euler-Bernoulli Piezoelectric Beam Finite Element

The convergence characteristic of the conventional two-noded Euler-Bernoulli piezoelectric beam finite element depends on the configuration of the beam cross-section. The element shows slower convergence for the asymmetric material distribution in the beam cross-section due to ‘material-locking’ caused by extension-bending coupling. Hence, the use of conventional Euler-Bernoulli beam finite element to analyze piezoelectric beams which are generally made of the host layer with asymmetrically surface bonded piezoelectric layers/patches, leads to increased computational effort to yield converged results. Here, an efficient coupled polynomial interpolation scheme is proposed to improve the convergence of the Euler-Bernoulli piezoelectric beam finite elements, by eliminating ill-effects of material-locking. The equilibrium equations, derived using a variational formulation, are used to establish relationships between field variables. These relations are used to find a coupled quadratic polynomial for axial displacement, having contributions from an assumed cubic polynomial for transverse displacement and assumed linear polynomials for layerwise electric potentials. A set of coupled shape functions derived using these polynomials efficiently handles extension-bending and electromechanical couplings at the field interpolation level itself in a variationally consistent manner, without increasing the number of nodal degrees of freedom. The comparison of results obtained from numerical simulation of test problems shows that the convergence characteristic of the proposed element is insensitive to the material configuration of the beam cross-section.


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, 2014 a;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.
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 fournoded 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 andNarayanan (2002), Bruent 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.
Latin American Journal of Solids and Structures 12 (2015) 153-172 (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 extensionbending 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, 2014 c;2014 d).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 extensionbending 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.

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 Latin American Journal of Solids and Structures 12 (2015) 153-172 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.

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.0 () ux and 0 () wxare 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:

Electric Potential and Electric Field
The through-thickness profile of the electric potential ( , ) xz  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 1 ()  at the top and bottom faces of the th i piezoelectric layer, respectively.The electric field in the transverse direction ( z E ) can be derived from electric potential as (Bendary et al., 2010): Latin American Journal of Solids and Structures 12 (2015) 153-172   ( where , is the difference of potentials at the top and bottom faces of the th i piezoelectric layer with thickness i h .

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 , ( , 1.....6) tively.The constitutive relations for such a system are given as (Sulbhewar and Raveendranath, 2014 c): where ( i =1….number of piezoelectric layers), ( k =1…..total number of layers in beam).The con- stants , Qeand denote elastic ( 2Nm ), piezoelectric ( 2Cm ) and dielectric ( Fm) properties reduced for the beam geometry (Sulbhewar and Raveendranath, 2014 c).

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.

Electromechanical and Strain Energy Variations
For the th j conventional/composite material layer, the mechanical strain energy variation is given as: The electromechanical strain energy variation of the th i piezoelectric layer is given as: Latin American Journal of Solids and Structures 12 (2015) 153-172 Substituting values of axial strain ( x  ), electric field ( i z E ) 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

Variation of Kinetic Energy
Total kinetic energy of the beam is given as (Sulbhewar and Raveendranath, 2014 c): where k  is volumic mass density of th k layer in 3 kgm  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 de- rive at: () denotes t .

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, 2014 c): .
Latin American Journal of Solids and Structures 12 (2015) 153-172 in which ,, f f f are volume, surface and point forces, respectively.0 q and S  are the charge density and area on which charge is applied.

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, 2014 c): 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 ( 0 u ) with transverse displacement ( 0 w ) and electric potential ( ) as: where 1 11 1 11 0 ( ) / ( ) 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.

FINITE ELEMENT FORMULATION
Using the variational formulation described above, a finite element model has been developed here.
The model consists of two mechanical variables ( 0 u and 0 w ) 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 ( 0 w ) and linear polynomials for i  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 Using these polynomials for 0 w and i  in equation ( 16) and integrating with respect to  , we get the quadratic polynomial for axial displacement 0 u as: Latin American Journal of Solids and Structures 12 (2015) 153-172 It is noteworthy that equation ( 18) contains a quadratic term in addition to the conventional linear terms.The formulation which use linear interpolation for 0 u (Bendary et al., 2010) is known to exhibit material-locking when applied to beams with asymmetric cross-section (Sulbhewar and Raveendranath, 2014 b).A higher-order interpolation for 0 u as given in equation ( 18), is expected to relieve material-locking and improve the rate of convergence (Sulbhewar and Raveendranath, 2014 b).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 1  u takes care of the bending effects, while 2 ui  takes care of electromechanical effects on the axial displacement.As seen from equa- tions ( 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: Latin American Journal of Solids and Structures 12 (2015) 153-172 where M is the mass matrix, ,,,  are global stiffness sub-matrices., U  are the glob- al 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 for- mulation 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).

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 PVDF (Sun and Huang, 2000):  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 8  and a mesh of PLANE 223 element for piezoelectric layer with size 100 2  is used.

Static Analysis: Sensor Configuration
For sensor configuration, the piezoelectric smart cantilever shown in Figure 2    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.

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.

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: Latin American Journal of Solids and Structures 12 (2015) 153-172  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.

Figure 1 :
Figure 1: Geometry of a general multilayered piezoelectric smart beam.
being the length of the beam element.
. The material and geometric properties of the beam are: Aluminum(Bendary et al., 2010):

Figure 2 :
Figure 2: Geometry of a extension mode piezoelectric smart cantilever beam.
is subjected to a tip load of 1000 N .The structural response of the beam is evaluated by present EB-CPI, conventional EB-IPI ofBendary 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 Figures3, 4 and 5, respectively.As seen from the plots, EB-CPI gives accurate predictions as of EB-IPI ofBendary 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 Figure6, which validate the use of present coupled polynomial interpolation for calculation of derived quantities also.

Figure 3 :
Figure 3: Sensor configuration: Variation of the transverse deflection along the length of the piezoelectric smart cantilever subjected to a tip load of 1000 N.

Figure 4 :
Figure 4: Sensor configuration: Variation of the axial deflection along the length of the piezoelectric smart cantilever subjected to a tip load of 1000 N.

Figure 7 :Figure 8 :Figure 9 :
Figure 7: Sensor configuration: Convergence characteristics of the Euler-Bernoulli based piezoelectric beam finite elements to predict the tip deflection of the piezoelectric smart cantilever subjected to a tip load of 1000 N.

Figure 10 :
Figure 10: Actuator configuration: Variation of the axial deflection along the length of the piezoelectric smart cantilever actuated by 10 KV.

Figure 12 :Figure 13 :
Figure 12: Modal analysis: Convergence characteristics of the Euler-Bernoulli based piezoelectric beam finite elements to predict the first natural frequency of the piezoelectric smart cantilever.

Table 1 :
Natural frequencies in Hz for piezoelectric smart cantilever beam (Refer Figure2).