A Timoshenko Piezoelectric Beam Finite Element with Consistent Performance Irrespective of Geometric and Material Configurations

The conventional Timoshenko piezoelectric beam finite elements based on First-order Shear Deformation Theory (FSDT) do not maintain the accuracy and convergence consistently over the applicable range of material and geometric properties. In these elements, the inaccuracy arises due to the induced potential effects in the transverse direction and inefficiency arises due to the use of independently assumed linear polynomial interpolation of the field variables in the longitudinal direction. In this work, a novel FSDT-based piezoelectric beam finite element is proposed which is devoid of these deficiencies. A variational formulation with consistent through-thickness potential is developed. The governing equilibrium equations are used to derive the coupled field relations. These relations are used to develop a polynomial interpolation scheme which properly accommodates the bending-extension, bending-shear and induced potential couplings to produce accurate results in an efficient manner. It is noteworthy that this consistently accurate and efficient beam finite element uses the same nodal variables as of conventional FSDT formulations available in the literature. Comparison of numerical results proves the consistent accuracy and efficiency of the proposed formulation irrespective of geometric and material configurations, unlike the conventional formulations.

Latin American Journal of Solids and Structures 13 (2016) 992-1015 2010; Kumar and Narayanan, 2008;Sulbhewar and Raveendranath, 2014a) based on the Euler-Bernoulli beam theory can be effectively used for the analysis of thin and slender piezoelectric smart beams.However, Euler-Bernoulli theory neglects the deformation due to shear and hence not suitable for thick and short beams.Sandwich Beam Theory (SBT) based analytical formulation (Zhang and Sun, 1996) and finite element formulations (Benjeddou et al., 1997(Benjeddou et al., , 2000;;Raja et al., 2002) considered the thicker core as a Timoshenko beam and the relatively thinner faces as Euler-Bernoulli beams.However, SBT is not suitable for short beams with thick piezoelectric layers.Timoshenko beam formulations based on the First-order Shear Deformation Theory (FSDT) consider constant shear strain across the beam cross-section.Analytical formulations (Abramovich, 1998;Aldraihem and Khdeir, 2000;Khdeir and Aldraihem, 2001) and finite elements (Robbins and Reddy, 1991;Shen, 1995;Narayanan and Balamurugan, 2003;Ray and Mallik, 2004;Neto et al., 2009) based on FSDT are widely used in the literature for the analysis of piezoelectric smart structures.
Accuracy of the conventional FSDT-based piezoelectric beam finite elements (Shen, 1995;Narayanan and Balamurugan, 2003;Ray and Malik, 2004;Neto et al., 2009) is adversely affected by the induced potential effects.These elements consider linear through-thickness distribution of electric potential which is actually nonlinear by virtue of the induced potential.The accuracy can be improved using assumed higher-order approximation of through-thickness electric potential (Jiang and Li, 2007;Kapuria and Hagedorn, 2007;Wang et al., 2007;Beheshti-Aval andLezgy-Nazargah, 2012, 2013).However, assumed higher-order potential distribution in the formulation introduces additional nodal electrical degrees of freedom in the transverse direction and hence increases the computational cost.An alternate efficient way to include the higher-order induced potential in FSDT-based formulation is to use the consistent through-thickness potential distribution derived from the electrostatic equilibrium equation (Sulbhewar and Raveendranath, 2014b).
Also, convergence of the conventional two-noded isoparametric FSDT-based piezoelectric beam element (Narayanan and Balamurugan, 2003) depends on the extent of the extension-bending and bending-shear couplings.Recently, Sulbhewar and Raveendranath (2015) proposed a novel FSDT piezoelectric beam finite element based on coupled polynomials for field variables which showed improved convergence.However, this element is not consistently accurate as the governing equations used to define coupled shape functions in this formulation are based on the assumed linear through-thickness potential.The associate errors are prominent for beams with piezoelectrically dominant cross-sections and/or with higher piezoelectric coefficients.An ideal FSDT-based formulation which is accurate and efficient over all geometric and material configurations of the piezoelectric beam should incorporate the induced potential coupling along with other mechanical couplings at the field interpolation level itself.
In the present work, an attempt is made to develop a novel FSDT piezoelectric beam formulation which is consistently accurate and efficient throughout the applicable range of geometric and material properties.The governing equations are derived using the variational formulation based on FSDT in conjunction with the consistent through-thickness potential.The relations established tween field variables are used to define coupled quadratic polynomials for axial displacement ( 0 u ) and section rotation ( ), having contributions from the assumed cubic polynomial for transverse displacement ( 0 w ) and assumed linear polynomials for layerwise electric potential variables ( i   ).
The shape functions based on these polynomials efficiently handle change in stiffness due to the Latin American Journal of Solids and Structures 13 ( 2016) 992-1015 induced potential along with bending-extension and bending-shear couplings, in an efficient manner.
Comparison of results from the present and the conventional formulations against ANSYS 2D benchmark simulation results proves the improved accuracy of the present formulation over the conventional formulations.Convergence studies are carried out to prove the improved convergence characteristics of the present FSDT element over the conventional isoparametric FSDT beam elements.It is noteworthy that owing to the fully coupled polynomial representation for section rotation and coupled quadratic term in the interpolation polynomial for axial displacement and transverse electric potential, the improved performance has been achieved with the same number of nodal degrees of freedom as of conventional two-noded isoparametric FSDT-based piezoelectric beam element.

THEORETICAL FORMULATION
An equivalent single layer (ESL) FSDT model for mechanical fields and a layerwise model for electric potential ( ) are employed for the proposed formulation.Consider a general multilayered extension mode piezoelectric smart beam with total number of layers n, as shown in Figure 1.The layers can be host layer(s) of conventional material or bonded/embedded layers of piezoelectric material.The beam layers are assumed to be made up of isotropic or specially orthotropic materials with perfect bonding among them.Top and bottom faces of 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 and assumptions of beam theory apply.

Reduced Constitutive Relations
For a general piezoelectric smart structure, the elastic constants relate the mechanical and electrical variables through the three-dimensional constitutive relations.An extension mode piezoelectric smart beam with axes of material symmetry parallel to the beam axes is considered here.For extension mode, the transversely poled piezoelectric material is subjected to the transverse electric field.For such a beam, the general constitutive relations are reduced according to beam theory, which are given as (Sulbhewar and Raveendranath 2014 b): where ( i =1….number of piezoelectric layers), ( k =1…..n)., , , ,D

Mechanical Displacements and Strains
The mechanical displacement fields in the longitudinal and transverse directions for FSDT are given as (Narayanan and Balamurugan, 2003): 0 ( ) u x and 0 ( ) w x are the centroidal axial and transverse displacements, respectively. is the sec- tion rotation of the beam.Dimensions , , L b h denote the length, width and the total thickness of the beam, respectively.
Axial and shear strain fields are derived using usual strain-displacement relations as: where ()' denotes derivative with respect to x .

Electric Potential and Electric Field
The layerwise two dimensional electric potential ( , )  (Sulbhewar and Raveendranath, 2014b): where The first two terms of expression (6) describe the conventional linear part in which i  and i   are the mean and difference, respectively, of the top and bottom surface potentials of the th i piezoe- lectric layer.The quadratic term represents the coupling between curvature strain and electric potential which constitutes induced potential.

VARIATIONAL FORMULATION
The formulation is based on Hamilton's principle which implicitly takes care of natural boundary conditions.It is expressed as (Chee et al., 1999): 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.

Variation of Electromechanical/Strain Energy
The electromechanical/strain energy variation of the piezoelectric smart beam is given as (Chee et al., 1999): Latin American Journal of Solids and Structures 13 (2016) 992-1015 Substituting values of axial strain ( x  ), shear strain ( xz  ), electric field ( i z E ) from equations (4), ( 5), (7) and using them along with constitutive relations given by equation (1) in expression (9); the total variation on the potential energy of the smart beam is given as: where

Variation of Kinetic Energy
Total kinetic energy of the beam is given as (Chee et al., 1999): where k  is the mass density of th k layer in 3 kgm  and ( k =1…n).Substituting values of u and w from equations (2) and (3) and applying variation, to derive at: where .

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 (Chee et al., 1999): 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 the field variables is established here using static governing equations.For static conditions without any external loading, the variational principle given in equation ( 8) reduces to (Sulbhewar and Raveendranath, 2015): Applying variation to the basic variables in equation ( 10), the static governing equations are obtained as: where Assuming that the higher order continuous derivatives of variables appearing in the governing equation ( 17) exist, we can write: Using equations ( 15) and ( 18), we can write the relationship of axial displacement ( 0 u ) with transverse displacement ( 0 w ) and electric potential variable ( i   ) as: where From equations ( 16)-( 19), we can write the relationship of section rotation ( ) with transverse displacement ( 0 w ) and electric potential variable ( i   ) as: where Latin American Journal of Solids and Structures 13 (2016) 992-1015 These equations for '' 0 u and  are used in the next Section to derive coupled polynomial ex- pressions for the field variables.It is clear that the coupling coefficients ( 1,2,3,4) j j   which depend on geometric and material properties of the beam, relate all the field variables by properly accommodating bending-extension, bending-shear and induced potential couplings.It is noteworthy appearing in equations ( 15)-( 17), which are used to define the coupling coefficients ( 1,2,3,4) j j   are different from those given in the Sulbhewar and Raveendranath (2015).The constants n m A in the present formulation contain additional stiffness terms (shown in curly braces) due to the induced potential effects.It may be noted that this induced stiffness is proportional to  which bears the same unit (N/m 2 ) as of elastic modulus 11 Q  .Hence, the quantity  may be termed as 'induced modulus'.

FINITE ELEMENT FORMULATION
Using the variational formulation described above, a finite element model is developed here.The two-noded beam element considered here is based on FSDT with layerwise electric potential in the transverse direction.There are three mechanical variables in the formulation namely, 0 0 , u w and  and layerwise electric potential variables i   where ( i =1.....number of piezoelectric layers in the beam).
The equations ( 19) and ( 20), derived using the governing equilibrium equations, demand continuous third order derivative of 0 w and first order derivative of i   .Hence, in terms of the natural coordinate  , a cubic polynomial for transverse displacement 0 w and linear polynomials for layerwise electric potential variable i   are assumed as given in equations (21 a) and (21 b), respectively.
The transformation between the local coordinate  and the global coordinate x along the length of the beam is given as  , length of the beam element.
Using equations (21 a) and (21 b) in equation ( 19) and integrating with respect to  , we get the coupled polynomial for axial displacement 0 u as: It is noted that the coupled quadratic term in equation ( 22) contains contributions from 0 w and i  fields and does not bring in any additional generalized degree of freedom.
Substituting equations (21 a) and (21 b) in equation ( 20), the coupled polynomial expression for the section rotation  is derived as: Latin American Journal of Solids and Structures 13 (2016) 992-1015 Equation ( 23) interpolates  by purely coupled terms with contributions from 0 w and i  fields.
It is noteworthy that equations ( 22) and ( 23) take care of extension-bending, bending-shear and induced potential couplings in a variationally consistent manner with the help of coupled terms present in the description of axial displacement and section rotation.
Using equations ( 21)-( 23), the coupled shape functions [ ] ( which relate the field variables to their nodal values as given in equation ( 24) are derived by usual method.The expressions for shape functions are given in Appendix.
As noted from the equation ( 24), while employing quadratic polynomials for axial displacement 0 u and section rotation  in the present FSDT formulation, the number and type of nodal variables are maintained the same as of the conventional isoparametric FSDT formulation.The variation on basic mechanical and electrical variables can now be transferred to nodal degrees of freedom.Substituting equation (24) in equations ( 10), ( 12), (13) and using them in equation (8), the following discretized form of the model is obtained: where M is mass matrix, , , ,   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.The matrix equations are now solved according to electrical conditions (open/closed circuit), configuration (actuator/sensor) and type of analysis (static/dynamic).

NUMERICAL EXAMPLES AND DISCUSSIONS
The software implementation of the present formulation has been carried out in MATLAB environment.The accuracy and efficiency of the proposed FSDT finite element are tested here for static (actuation/sensing) and modal (open/closed circuit) analyses and its performance is compared Latin American Journal of Solids and Structures 13 (2016) 992-1015 against the conventional two-noded isoparametric FSDT piezoelectric beam finite element available in the literature.The following designations are used:  FSDT-Coupled: The present formulation which uses the coupled polynomials (cubic for 0 w given by equation (21 a), coupled quadratic for 0 u given by equation ( 22), coupled quadratic for  given by equation ( 23) and linear for i   given by equation (21 b)) for inter- polation of field variables and layerwise consistent through-thickness potential ( coupled quadratic approximation in z direction given by equation ( 6)). FSDT: The conventional FSDT formulation of Narayanan and Balamurugan (2003) which uses independent polynomials for field interpolation (linear for 0 0 , , u w  and i   ) and layer- wise assumed linear through-thickness potential. ANSYS 2D: For a comparative evaluation of the above FSDT formulations, benchmark solutions have been obtained from a refined two-dimensional analysis using ANSYS finite element software, for which PLANE 183 elements are used to mesh conventional material layers, while PLANE 223 elements are used to mesh piezoelectric material layers.

Example 1: A Symmetric Bimorph Beam
A bimorph cantilever beam with oppositely poled piezoelectric layers as shown in Figure 2 is considered here.In order to study the effect of material properties on the performance of the FSDT elements, the following materials are used while the geometry is fixed ( 10 , 100 PVDF (Sun and Huang, 2000):  (Kapuria and Hagedorn, 2007):   For a comparative evaluation of accuracy of various FSDT-based formulations, converged results from a refined mesh of 40 elements have been used.Converged results from an ANSYS 2D simulation with a mesh of 100 10  elements are used as a benchmark.Static Analysis-Actuator Configuration: In this configuration, the interface of the bimorph is grounded and the voltages of 10volts  are applied on the free surfaces.Table 1 and 2 show the results for the tip deflection and the maximum axial stress developed in the bimorphs of different materials, respectively.Also, the associated absolute errors calculated with respect to ANSYS 2D benchmark solutions are presented in brackets.As seen from the tables, the conventional FSDT formulation (Narayanan and Balamurugan, 2003) fails to produce consistently accurate results.The percentage errors increase with the modulus ratio   (the ratio of induced modulus to elastic modulus) of the material.The present FSDT-Coupled element predicts accurate results for all the bimorphs, irrespective of the modulus ratio.This consistent performance of the present formulation can be attributed to the accommodation of induced potential effects through the coupled interpolation polynomials.Static Analysis-Sensor Configuration: Here, the beam shown in Figure 2 is subjected to a tip load of 1000 N. The results for the tip deflection, potential developed at the mid-span and the maximum axial stress developed at the root of the bimorphs of different materials are tabulated in the Tables 3, 4 and 5, respectively.The associated absolute errors (in percentage) with respect to AN-SYS 2D benchmark solutions are presented in brackets.As seen from the tables, the present FSDT-Coupled consistently reproduces the ANSYS 2D simulation results for all the bimorphs, unlike the conventional formulation.The accuracy of the present FSDT-Coupled formulation is practically insensitive to the material properties of the beam.The convergence graphs plotted in Figures 3 and 4 for the tip deflection and potential developed at the root, respectively, compare the efficiency of the FSDT-based piezoelectric beam finite elements.The G1195N bimorph which has the highest modulus ratio among the chosen materials is taken as a particular example for this study.The FSDT-Coupled shows single-element convergence, closely reproducing the ANSYS-2D solutions for both the tip deflection and the potential developed.The conventional FSDT (Narayanan and Balamurugan, 2003) overestimates the response as it neglects induced potential effects.
Latin American Journal of Solids and Structures 13 (2016) 992-1015   6 reveal the inability of conventional FSDT formulation to maintain the accuracy over the different bimorph materials.The consistent accuracy of the present FSDT-Coupled results validates the use of coupled polynomial shape functions in generating the element mass matrix consistent with the element stiffness matrix.

Example 2: A Two-Layer Asymmetric Piezoelectric Beam
A two-layer asymmetric piezoelectric cantilever beam having a steel host layer with a surface bonded piezoelectric layer of G1195N at the top, as shown in Figure 7 is considered here.The material properties used are: Steel (Carrera and Brischetto, 2008):  Static Analysis-Sensor Configuration: In this configuration, the beam shown in Figure 7 is subjected to a tip load of -1000 N. The results for transverse tip deflection, axial tip deflection and potential developed across the piezoelectric layer at the mid-span of the beam for various thickness ratios are tabulated in Tables 7, 8 and 9, respectively.The present FSDT-Coupled formulation proves its versatility, yielding consistently accurate predictions over the entire range of thickness ratio.The conventional FSDT formulation (Narayanan and Balamurugan, 2003) does not maintain the consistent accuracy, as it neglects the induced potential coupling.The associated error increases rapidly in the higher thickness ratio regimes.The convergence graphs plotted in Figures 8 and 9 for the transverse tip deflection and the potential developed at the root, respectively, prove the consistent efficiency of the present FSDT-Coupled formulation, which exhibits single-element convergence to ANSYS-2D solutions.FSDT (Narayanan and Balamurugan, 2003) model shows very slow convergence to the inaccurate results, due to induced potential effects.Static Analysis-Actuator Configuration: Here, the beam shown in Figure 7 is actuated by 100 volts.The variations of transverse and axial deflections at the tip, with thickness ratio are tabulated in Tables 10 and 11, respectively.The FSDT-Coupled formulation consistently gives accurate predictions of results over the entire range of thickness ratio as given by ANSYS 2D simulation.The conventional FSDT formulation (Narayanan and Balamurugan, 2003) does not yield consistently accurate results.Modal Analysis: The FSDT-Coupled formulation is evaluated here for its accuracy and efficiency to predict the natural frequencies of the asymmetric piezoelectric smart beam.The first natural frequency of the asymmetric Steel/G1195N beam shown in Figure 7 is computed for both open and closed circuit electrical boundary conditions.Table 12 shows the variation of first natural frequencies with thickness ratio.The results of FSDT-Coupled formulation agree very well with the ANSYS 2D simulation results.The results of the conventional FSDT formulation (Narayanan and Balamurugan, 2003), show significant deviation in the higher thickness ratio regimes where the induced potential effects are predominant.
The consistent efficiency of the present FSDT-Coupled is revealed by the convergence graphs for first natural frequency in both open and closed circuit electrical boundary conditions plotted in Figures 10 and 11, respectively.As seen from the figures, the FSDT-Coupled gives fast convergence, while FSDT (Narayanan and Balamurugan, 2003) model shows very slow convergence to the inaccurate results, due to induced potential effects.

CONCLUSIONS
Based on coupled polynomial field interpolations in conjunction with a consistent through-thickness electric potential, a novel FSDT based extension mode piezoelectric beam finite element has been presented.The derived set of coupled shape functions handles bending-extension, bending-shear and induced potential couplings in a variationally consistent manner.Numerical evaluation has proven the merits of the present formulation over the conventional formulations available in the literature, in terms of accuracy and efficiency.From the numerical analysis, it was found that:  The performance of the conventional FSDT piezoelectric beam finite elements depends on the geometric and material parameters of the beam.The accuracy depends on the proportion of piezoelectric material in the total beam thickness (thickness ratio) and the ratio of induced modulus to the elastic modulus (modulus ratio).The convergence rate of conventional FSDT elements is deteriorated by the presence of bending-shear and bendingextension couplings. The performance of the proposed FSDT-Coupled formulation proves to be insensitive to the material and geometric configuration of the beam cross-section.It consistently maintains the level of accuracy and efficiency for all modulus and thickness ratios.

Figure 1 :
Figure 1: Geometry of a general multilayered piezoelectric smart beam.

Figure 3 :
Figure 3: Example 1: Sensor configuration: Convergence characteristics of the FSDT-based piezoelectric beam finite elements to predict the tip deflection of the G1195N bimorph cantilever beam subjected to a tip load of 1000 N.

Figure 4 :
Figure 4: Example 1: Sensor configuration: Convergence characteristics of the FSDT-based piezoelectric beam finite elements to predict the potential developed at the root of the G1195N bimorph cantilever beam subjected to a tip load of 1000 N.

Figures 5
Figures 5 and 6 show the comparison of convergence characteristics of FSDT-based piezoelectric beam finite element formulations to predict the first natural frequency of the G1195N bimorph in open and closed circuit conditions, respectively.FSDT-Coupled shows quick convergence, closely reproducing the ANSYS-2D solutions for both open and closed circuit conditions.The conventional FSDT (Narayanan and Balamurugan, 2003) model underestimates the response as it neglects induced potential effects.

Figure 5 :
Figure 5: Example 1: Modal analysis: Convergence characteristics of the FSDT-based piezoelectric beam finite elements to predict the first natural frequency of the G1195N bimorph cantilever beam in open circuit electrical boundary condition.

Figure 6 :
Figure 6: Example 1: Modal analysis: Convergence characteristics of the FSDT-based piezoelectric beam finite elements to predict the first natural frequency of the G1195N bimorph cantilever beam in closed circuit electrical boundary condition.

Figure 8 :Figure 9 :
Figure 8: Example 2: Sensor configuration: Convergence characteristics of FSDT-based formulations to predict the transverse tip deflection of the two-layer asymmetric piezoelectric cantilever beam (r=0.5)subjected to a tip load of -1000 N.

Figure 10 :
Figure 10: Example 2: Modal Analysis: Convergence characteristics of FSDT-based formulations to predict the first natural frequency of the two-layer asymmetric piezoelectric cantilever beam (r=0.5) in open circuit electrical boundary condition.

Figure 11 :
Figure 11: Example 2: Modal Analysis: Convergence characteristics of FSDT-based formulations to predict the first natural frequency of the two-layer asymmetric piezoelectric cantilever beam (r=0.5) in closed circuit electrical boundary condition.

Table 2 :
Example 1: Absolute maximum axial stress developed (kPa) in the bimorph cantilever beams of different piezoelectric materials actuated by ±10 volts.(The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Table 3 :
Example 1: Absolute tip deflection ( m  ) of the bimorph cantilever beams of different piezoelectric materials subjected to a tip load of 1000 N. (The absolute errors in percentage are given with respect to ANSYS 2D simulation.)Latin American Journal of Solids and Structures 13 (2016) 992-1015

Table 4 :
Example 1: Potential developed (volts) at the mid-span of the bimorph cantilever beams of different piezoelectric materials subjected to a tip load of 1000 N. (The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Table 5 :
Example 1: Absolute maximum axial stress developed (MPa) at the root of the bimorph cantilever beams of different piezoelectric materials subjected to a tip load of 1000 N. (The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Table 6 :
Example 1: First natural frequency (Hz) of the bimorph cantilever beams of different piezoelectric materials with open and closed circuit electrical boundary conditions.(The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Table 8 :
Example 2: Axial tip deflection ( m  ) of the asymmetric piezoelectric cantilever beam (Steel/G1195N) subjected to a tip load of -1000 N. (The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Table 9 :
Example 2: Potential developed (volts) at the mid-span of the asymmetric piezoelectric cantilever beam (Steel/G1195N) subjected to a tip load of -1000 N. (The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Table 10 :
Example 2: Transverse tip deflection ( m  ) of the asymmetric piezoelectric cantilever beam (Steel/G1195N) actuated by 100 volts.(The absolute errors in percentage are given with respect to ANSYS 2D simulation.)Latin American Journal of Solids and Structures 13 (2016) 992-1015

Table 11 :
Example 2: Axial tip deflection ( m  ) of the asymmetric piezoelectric cantilever beam (Steel/G1195N) actuated by 100 volts.(The absolute errors in percentage are given with respect to ANSYS 2D simulation.)