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

Abstract

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.

Keywords:
Coupled field; finite elements; induced potential; piezoelectric; smart structures

1 INTRODUCTION

Beam formulations are widely used for the numerical analysis of one dimensional piezoelectric structure (Marinkovic and Marinkovic, 2012). Analytical formulations (Crawley and de Luis, 1987; Abramovich and Pletner, 1997; Crawley and Anderson, 1990) and finite elements (Bendary et al., 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, 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 and Lezgy-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 between field variables are used to define coupled quadratic polynomials for axial displacement (u_o ) and section rotation (θ), having contributions from the assumed cubic polynomial for transverse displacement (w_o ) and assumed linear polynomials for layerwise electric potential variables . The shape functions based on these polynomials efficiently handle change in stiffness due to the 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.

2 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.

Figure 1
Geometry of a general multilayered piezoelectric smart beam.

2.1 Reduced Constitutive Relations

For a general piezoelectric smart structure, the elastic Cij (i, j = 1...6), piezoelectric ekj(k = 1,2,3, j=1....6) and dielectric ∈k (k = 1,2,3) 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):

(1)

where (i=1....number of piezoelectric layers), (k=1.....n). σ, τ, ε, γ, D and E denote the axial stress (N/m2), shear stress (N/m2), normal strain, shear strain, electric displacement (C/m2) and electric field (V/m), respectively. The constants i = 1,5), ẽ and N/m ²⁾, piezoelectric (C/m ²⁾ and dielectric (F/m) properties, respectively. denote reduced elastic ( (

2.2 Mechanical Displacements and Strains

The mechanical displacement fields in the longitudinal and transverse directions for FSDT are given as (Narayanan and Balamurugan, 2003):

(2)

(3)

u₀(x) and w₀(x) are the centroidal axial and transverse displacements, respectively. θ is the section 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:

(4)

(5)

where ()' denotes derivative with respect to x.

2.3 Electric Potential and Electric Field

The layerwise two dimensional electric potential φ_i (x,z) takes the values of Φ_{i + 1} (x) and Φ_i (x) at the top and bottom faces of ^_ith piezolayer, respectively as shown in Figure 1. The through-thickness distribution of electric potential φ_i (x,z) consistent with FSDT is used (Sulbhewar and Raveendranath, 2014b):

(6)

where

The first two terms of expression (6) describe the conventional linear part in which ^_ith piezoelectric layer. The quadratic term represents the coupling between curvature strain and electric potential which constitutes induced potential. and are the mean and difference, respectively, of the top and bottom surface potentials of the

The layerwise electric field (_^Eiz ) is obtained from equation (6) as (Sulbhewar and Raveendranath, 2014b):

(7)

3 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):

(8)

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.

3.1 Variation of Electromechanical/Strain Energy

The electromechanical/strain energy variation of the piezoelectric smart beam is given as (Chee et al., 1999):

(9)

Substituting values of axial strain (ε_x ), shear strain (γ_xz ), electric field (_^Eiz ) 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:

(10)

where _^Ikq = bq + 1)./(

3.2 Variation of Kinetic Energy

Total kinetic energy of the beam is given as (Chee et al., 1999):

(11)

where ρ_k is the mass density of ^_kth layer in kgm ^-3 and (k =1...n). Substituting values of u and w from equations (2) and (3) and applying variation, to derive at:

(12)

where ^(˙) denotes ∂/∂t.

3.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 (Chee et al., 1999):

(13)

in which ^{_{fV, fS, fC}} are volume, surface and point forces, respectively. q ₀ and S _φ are the charge density and area on which charge is applied.

4 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):

(14)

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

(15)

(16)

(17)

where

Assuming that the higher order continuous derivatives of variables appearing in the governing equation (17) exist, we can write:

(18)

Using equations (15) and (18), we can write the relationship of axial displacement (u ₀) with transverse displacement (w ₀) and electric potential variable as:

(19)

where β₁ = _{^Au 2}/_{^Au 1} and βⁱ ₂ = _{^Aui 3}/_{^Au 1}.

From equations (16)-(19), we can write the relationship of section rotation (θ) with transverse displacement (w ₀) and electric potential variable as:

(20)

where . and

These equations for u"0 and θ are used in the next Section to derive coupled polynomial expressions for the field variables. It is clear that the coupling coefficients βj(j = 1,2, 3, 4) 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 that the constants Anm(m = 1;2 3; 4) (n = u, θ w) appearing in equations (15)-(17), which are used to define the coupling coefficients βj(j = 1,2, 3, 4) are different from those given in the (Sulbhewar and Raveendranath, 2015). The constants Anm 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 ^_N/m2 ) as of elastic modulus which bears the same unit (. Hence, the quantity may be termed as 'induced modulus'.

5 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, u ₀, w ₀and θ and layerwise electric potential variables i =1.....number of piezoelectric layers in the beam). where (

The equations (19) and (20), derived using the governing equilibrium equations, demand continuous third order derivative of w ₀ and first order derivative of (21a) and (21b), respectively. The transformation between the local coordinate ξ and the global coordinate x along the length of the beam is given as ξ = [2(x - x ₁) / (x ₂ - x ₁)] - 1 and (x ₂ - x ₁) - l, length of the beam element.. Hence, in terms of the natural coordinate ξ, a cubic polynomial for transverse displacement w0 and linear polynomials for layerwise electric potential variable are assumed as given in equations

(21a)

(21b)

Using equations (21a) and (21b) in equation (19) and integrating with respect to ξ, we get the coupled polynomial for axial displacement u ₀ as:

(22)

It is noted that the coupled quadratic term in equation (22) contains contributions from w ₀ and fields and does not bring in any additional generalized degree of freedom.

Substituting equations (21a) and (21b) in equation (20), the coupled polynomial expression for the section rotation θ is derived as:

(23)

Equation (23) interpolates θ by purely coupled terms with contributions from w ₀ and 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 [_^Num ](m = 1..8), [_^Nwm ](m = 1..6), [N ^θ _m ](m = 1..6) and [N ^{φ i}_m ](m = 1..2) 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.

(24)

As noted from the equation (24), while employing quadratic polynomials for axial displacement 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:

(25)

where M is 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. The matrix equations are now solved according to electrical conditions (open/closed circuit), configuration (actuator/sensor) and type of analysis (static/dynamic).

6 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 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 w ₀ given by equation (21a), coupled quadratic for u ₀ given by equation (22), coupled quadratic for θ given by equation (23) and linear for (21b)) for interpolation of field variables and layerwise consistent through-thickness potential (coupled quadratic approximation in z direction given by equation (6)). given by equation

• FSDT: The conventional FSDT formulation of (Narayanan and Balamurugan, 2003) which uses independent polynomials for field interpolation (linear for u _0, w _0, θ and ) and layerwise 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.

6.1 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 (h = 10mm, L = 100mm).

PVDF (Sun and Huang, 2000):

E = 2Gpa, υ = 0.29, e - 0.046Cm-2, ∈ = 0.1062 × 10-9Fm-1, ρ = 1880kgm-3

PZT 2 (EFunda.com, 2014)

{C11, C12, C13, C22, C23, C33, C44, C55, C66} = {134.87, 67.89, 68.09, 134.87, 68.09, 113.30, 22.22, 22.22, 33.44}GPa

{e31, e32, e33} = {-1.8160, -1.8160, 90506}Cm-2 {∈1, ∈2, ∈3} = {8.7655, 8.7655, 3.9843} × 10-9Fm-1

ρ = 7600kgm-3

PZT 4 (EFunda.com, 2014)

{C11, C12, C13, C22, C23, C33, C44, C55, C66} = {139, 77.84, 74.28, 139, 74.28, 115.41, 25.64, 25.64, 30.58}GPa

{e31, e32, e33} = {-5.2028, -5.2028, 15.0804}Cm-2 {∈1, ∈2, ∈3} = {1.3060, 1.3060, 1.1510} × 10-8Fm-1

ρ = 7500kgm-3

PZT-5H (Kapuria and Hagedorn, 2007):

{C11, C12, C13, C22, C23, C33, C44, C55, C66} = {126, 79.5, 84.1, 126, 84.1, 117, 23, 23, 23.25}GPa

{e31, e32, e33} = {-6.5, -6.5, 23.3}Cm-2 {∈1, ∈2, ∈3} = {1.503, 1.503, 1.3} × 10-8Fm-1

ρ = 7500kgm-3

PZT 5A (EFunda.com, 2014)

{C11, C12, C13, C22, C23, C33, C44, C55, C66} = {120.35, 75.18, 75.09, 120.35, 75.09, 110.87, 21.05, 21.05, 22.57}GPa

{e31, e32, e33} = {-5.3512, -5.3512, 15.7835}Cm-2 {∈1, ∈2, ∈3} = {1.5317, 1.5317, 1.5052} × 10-8Fm-1

ρ = 7550kgm-3

PZT 8 (EFunda.com, 2014)

{C11, C12, C13, C22, C23, C33, C44, C55, C66} = {146.88, 81.09, 81.05, 146.88, 81.05, 131.71, 31.35, 31.35, 32.89}GPa

{e31, e32, e33} = {-3.8754, -3.8754, 13.9108}Cm-2 {∈1, ∈2, ∈3} = {1.1422, 1.1422, 0.8854} × 10-8Fm-1

ρ = 7600kgm-3

G1195N (Peng et al., 1998): E = 63GPa, υ = 0.3, d31 = 254×10-12mV-1, e31 = E×d31Cm-2, ∈3 = 15×10-9Fm-1, ρ = 7600kgm-3

Figure 2
Example 1: Geometry of a bimorph cantilever beam.

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 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. (the ratio of induced modulus to elastic modulus) of the material. The present FSDT-

Table 1
Example 1: Absolute tip deflection (μm) of 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.)

Materials	Modulus ratio	FSDT (Narayanan and Balamurugan, 2003)	ANSYS 2D	FSDT-Coupled
PVDF	0.0100	0.0690 (0.291 %)	0.0688	0.0688 (0.000 %)
PZT-2	0.0589	0.1800 (1.408 %)	0.1775	0.1774 (0.056 %)
PZT-8	0.0687	0.2910 (1.677 %)	0.2862	0.2861 (0.035 %)
PZT-4	0.0773	0.3690 (1.906 %)	0.3621	0.3620 (0.028 %)
PZT-5A	0.0849	0.5130 (2.090 %)	0.5025	0.5023 (0.040 %)
PZT-5H	0.1751	0.8244 (4.328 %)	0.7902	0.7899 (0.038 %)
G1195N	0.2710	0.7620 (6.723 %)	0.7140	0.7137 (0.042 %)

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.)

Materials	Modulus ratio	FSDT (Narayanan and Balamurugan, 2003)	ANSYS 2D	FSDT-Coupled
PVDF	0.0100	0.0460 (0.648 %)	0.0463	0.0463 (0.000 %)
PZT-2	0.0589	5.1724 (4.174 %)	5.3977	5.3977 (0.000 %)
PZT-8	0.0687	8.4348 (4.810 %)	8.8610	8.8618 (0.009 %)
PZT-4	0.0773	10.000 (5.366 %)	10.567	10.569 (0.019 %)
PZT-5A	0.0849	10.427 (5.868 %)	11.077	11.077 (0.000 %)
PZT-5H	0.1751	16.492 (11.18 %)	18.567	18.567 (0.000 %)
G1195N	0.2710	16.002 (15.99 %)	19.048	19.048 (0.000 %)

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 ANSYS 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.

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.)

Materials	Modulus ratio	FSDT (Narayanan and Balamurugan, 2003)	ANSYS 2D	FSDT-Coupled
PVDF	0.0100	2000.0 (0.090 %)	1998.2	2000.0 (0.090 %)
PZT-2	0.0589	44.885 (1.456 %)	44.241	44.268 (0.061 %)
PZT-8	0.0687	44.065 (1.373 %)	43.468	43.364 (0.239 %)
PZT-4	0.0773	46.894 (1.661 %)	46.128	46.061 (0.145 %)
PZT-5A	0.0849	62.147 (1.752 %)	61.077	60.942 (0.221 %)
PZT-5H	0.1751	59.349 (3.757 %)	57.200	57.156 (0.077 %)
G1195N	0.2710	53.180 (5.247 %)	50.529	50.368 (0.319 %)

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.)

Materials	Modulus ratio	FSDT (Narayanan and Balamurugan, 2003)	ANSYS 2D	FSDT-Coupled
PVDF	0.0100	1612.2 (0.261 %)	1608.0	1608.3 (0.019 %)
PZT-2	0.0589	81.813 (1.411 %)	80.675	80.675 (0.000 %)
PZT-8	0.0687	58.055 (1.642 %)	57.117	57.123 (0.011 %)
PZT-4	0.0773	54.778 (1.835 %)	53.791	53.795 (0.007 %)
PZT-5A	0.0849	57.410 (2.001 %)	56.284	56.287 (0.005 %)
PZT-5H	0.1751	70.381 (3.871 %)	67.758	67.758 (0.000 %)
G1195N	0.2710	105.55 (5.628 %)	99.926	99.924 (0.002 %)

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.)

Materials	Modulus ratio	FSDT (Narayanan and Balamurugan, 2003)	ANSYS 2D	FSDT-Coupled
PVDF	0.0100	5.955 (0.750 %)	6.000	5.993 (0.117 %)
PZT-2	0.0589	5.886 (2.517 %)	6.038	5.993 (0.745 %)
PZT-8	0.0687	5.873 (2.861 %)	6.046	5.994 (0.860 %)
PZT-4	0.0773	5.862 (3.011 %)	6.044	5.994 (0.827 %)
PZT-5A	0.0849	5.852 (3.305 %)	6.052	5.994 (0.958 %)
PZT-5H	0.1751	5.741 (5.076 %)	6.048	5.994 (0.893 %)
G1195N	0.2710	5.636 (6.812 %)	6.048	5.994 (0.893 %)

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.

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
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.

Modal Analysis: The accuracy and efficiency of the FSDT elements in predicting the natural frequencies of the bimorph cantilever beam shown in Figure 2 are compared here. The natural frequencies are evaluated for closed and open circuit electrical boundary conditions, with different materials. For open circuit, only the interface of the bimorph is grounded while, for closed circuit all the faces of bimorph are grounded. The results tabulated in Table 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.

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.)

Materials	Modulus ratio	FSDT (Narayanan and Balamurugan, 2003)	ANSYS 2D	FSDT-Coupled
Open circuit
PVDF	0.0100	169.74 (0.135 %)	169.97	169.94 (0.018 %)
PZT-2	0.0589	550.61 (0.736 %)	554.69	554.41 (0.051 %)
PZT-8	0.0687	556.12 (0.822 %)	560.73	560.59 (0.025 %)
PZT-4	0.0773	542.14 (0.852 %)	546.80	547.38 (0.106 %)
PZT-5A	0.0849	463.68 (1.012 %)	468.42	468.22 (0.043 %)
PZT-5H	0.1751	482.37 (1.814 %)	491.28	491.50 (0.045 %)
G1195N	0.2710	506.16 (2.628 %)	519.82	520.03 (0.040 %)
Closed circuit
PVDF	0.0100	169.11 (0.142 %)	169.35	169.32 (0.018 %)
PZT-2	0.0589	539.01 (0.739 %)	543.02	542.90 (0.022 %)
PZT-8	0.0687	542.48 (0.960 %)	547.74	547.06 (0.124 %)
PZT-4	0.0773	527.63 (0.991 %)	532.91	532.63 (0.053 %)
PZT-5A	0.0849	449.75 (1.189 %)	455.16	454.44 (0.158 %)
PZT-5H	0.1751	453.81 (2.270 %)	464.35	463.53 (0.177 %)
G1195N	0.2710	461.90 (3.206 %)	477.20	477.12 (0.017 %)

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
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
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.

6.2 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): E = 210GPa, υ = 0.3, ρ = 7850kgm ^-3

G1195N (Peng et al., 1998): E = 63GPa, υ = 0.3, d ₃₁ = 254×10^-12 mV ^-1, e ₃₁ = E×d ₃₁ Cm ^-2, ∈₃ = 15×10^-9 Fm ^-1, ρ = 7600kgm ^-3.

Figure 7
Example 2: Geometry of a two-layer asymmetric piezoelectric cantilever beam.

The length and total height of the beam are fixed (L = 100mm, h = 5mm), while thicknesses of the piezoelectric layer (h_p ) and the host layer (h_c ) are varied. The performance of the FSDT-based piezoelectric beam finite elements is evaluated over a wide range of the piezoelectric material proportion in the total beam thickness (thickness ratio:r = h_p / h). For a comparative evaluation of accuracy of various FSDT based formulations, converged results from a refined mesh of 40 elements have been used. The converged results from an ANSYS 2D simulation with a mesh of 200×20 elements are used as a benchmark.

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.

Table 7
Example 2: Transverse tip deflection (mm) 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.)

Thickness ratio (r)	ANSYS 2D	FSDT-Coupled	FSDT (Narayanan and Balamurugan, 2003)
0.05	-0.1671	-0.1671 (0.000 %)	-0.1671 (0.000 %)
0.1	-0.1819	-0.1819 (0.000 %)	-0.1819 (0.000 %)
0.2	-0.2109	-0.2108 (0.047 %)	-0.2110 (0.047 %)
0.3	-0.2355	-0.2354 (0.042 %)	-0.2362 (0.297 %)
0.4	-0.2523	-0.2521 (0.079 %)	-0.2543 (0.793 %)
0.5	-0.2604	-0.2601 (0.115 %)	-0.2647 (1.651 %)
0.6	-0.2622	-0.2618 (0.153 %)	-0.2699 (2.937 %)
0.7	-0.2630	-0.2625 (0.190 %)	-0.2757 (4.829 %)
0.8	-0.2704	-0.2700 (0.148 %)	-0.2914 (7.766 %)
0.9	-0.2990	-0.2987 (0.100 %)	-0.3378 (12.98 %)
1.0	-0.4008	-0.4005 (0.075 %)	-0.5088 (26.95 %)

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.)

Thickness ratio (r)	ANSYS 2D	FSDT-Coupled	FSDT (Narayanan and Balamurugan, 2003)
0.05	0.1896	0.1897 (0.053 %)	0.1897 (0.053 %)
0.1	0.4040	0.4042 (0.050 %)	0.4043 (0.074 %)
0.2	0.8915	0.8918 (0.034 %)	0.8926 (0.123 %)
0.3	1.4065	1.4063 (0.014 %)	1.4111 (0.327 %)
0.4	1.8648	1.8631 (0.091 %)	1.8792 (0.772 %)
0.5	2.1898	2.1812 (0.393 %)	2.2197 (1.365 %)
0.6	2.3170	2.3148 (0.095 %)	2.3867 (3.008 %)
0.7	2.2523	2.2522 (0.004 %)	2.3656 (5.030 %)
0.8	1.9790	1.9803 (0.065 %)	2.1376 (8.014 %)
0.9	1.4033	1.4044 (0.078 %)	1.5885 (13.20 %)
1.0	0.0000	0.0000 (0.000 %)	0.0000 (0.000 %)

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.)

Thickness ratio (r)	ANSYS 2D	FSDT-Coupled	FSDT (Narayanan and Balamurugan, 2003)
0.05	16.356	16.355 (0.006 %)	16.355 (0.006 %)
0.1	34.840	34.848 (0.023 %)	34.851 (0.032 %)
0.2	76.880	76.880 (0.000 %)	76.950 (0.091 %)
0.3	121.25	121.24 (0.008 %)	121.66 (0.338 %)
0.4	160.62	160.62 (0.000 %)	162.01 (0.865 %)
0.5	188.05	188.05 (0.000 %)	191.36 (1.761 %)
0.6	199.57	199.57 (0.000 %)	205.76 (3.102 %)
0.7	194.17	194.16 (0.005 %)	203.94 (5.032 %)
0.8	170.73	170.72 (0.006 %)	184.29 (7.942 %)
0.9	121.08	121.07 (0.008 %)	136.95 (13.11 %)
1.0	000.00	000.00 (0.000 %)	000.00 (0.000 %)

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.

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 9
Example 2: Sensor configuration: Convergence characteristics of FSDT-based formulations to predict the potential developed at root of the two-layer asymmetric piezoelectric cantilever beam (r=0.5) subjected to a tip load of -1000 N.

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.

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.)

Thickness ratio (r)	ANSYS 2D	FSDT-Coupled	FSDT (Narayanan and Balamurugan, 2003)
0.05	9.9834	9.9843 (0.009 %)	9.9846 (0.012 %)
0.1	10.839	10.838 (0.009 %)	10.842 (0.028 %)
0.2	12.463	12.463 (0.000 %)	12.475 (0.096 %)
0.3	13.673	13.674 (0.007 %)	13.725 (0.380 %)
0.4	14.099	14.101 (0.014 %)	14.238 (0.986 %)
0.5	13.557	13.558 (0.007 %)	13.828 (1.999 %)
0.6	12.156	12.156 (0.000 %)	12.582 (3.504 %)
0.7	10.180	10.180 (0.000 %)	10.747 (5.570 %)
0.8	7.8347	7.8349 (0.002 %)	8.5019 (8.516 %)
0.9	4.9635	4.9638 (0.006 %)	5.6344 (13.52 %)
1.0	0.0000	0.0000 (0.000 %)	0.0000 (0.000 %)

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.)

Thickness ratio (r)	ANSYS 2D	FSDT-Coupled	FSDT (Narayanan and Balamurugan, 2003)
0.05	-0.1752	-0.1752 (0.000 %)	-0.1752 (0.000 %)
0.1	-0.2004	-0.2005 (0.049 %)	-0.2005 (0.049 %)
0.2	-0.2583	-0.2584 (0.039 %)	-0.2585 (0.077 %)
0.3	-0.3198	-0.3201 (0.094 %)	-0.3206 (0.250 %)
0.4	-0.3763	-0.3762 (0.027 %)	-0.3778 (0.399 %)
0.5	-0.4168	-0.4170 (0.048 %)	-0.4206 (0.912 %)
0.6	-0.4378	-0.4388 (0.228 %)	-0.4450 (1.645 %)
0.7	-0.4445	-0.4455 (0.225 %)	-0.4537 (2.069 %)
0.8	-0.4459	-0.4461 (0.045 %)	-0.4546 (1.951 %)
0.9	-0.4535	-0.4542 (0.154 %)	-0.4599 (1.411 %)
1.0	-0.5080	-0.5080 (0.000 %)	-0.5080 (0.000 %)

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.

Table 12
Example 2: Natural frequencies of the asymmetric piezoelectric cantilever beam (Steel/G1195N) in open and closed circuit electrical boundary conditions. (The absolute errors in percentage are given with respect to ANSYS 2D simulation.)

Thickness ratio (r)	ANSYS 2D	FSDT-Coupled	FSDT (Narayanan and Balamurugan, 2003)
Open circuit
0.05	399.07	398.93 (0.035 %)	398.93 (0.035 %)
0.1	382.81	382.67 (0.037 %)	382.66 (0.039 %)
0.2	356.16	356.03 (0.037 %)	355.87 (0.081 %)
0.3	337.57	337.45 (0.036 %)	336.89 (0.201 %)
0.4	326.64	326.57 (0.021 %)	325.17 (0.450 %)
0.5	322.02	322.03 (0.003 %)	319.24 (0.863 %)
0.6	321.37	321.48 (0.034 %)	316.61 (1.481 %)
0.7	321.38	321.54 (0.050 %)	313.76 (2.371 %)
0.8	317.72	317.58 (0.044 %)	305.71 (3.780 %)
0.9	302.33	302.41 (0.026 %)	284.40 (5.931 %)
1.0	261.63	261.62 (0.003 %)	232.14 (11.27 %)
Closed circuit
0.05	396.50	396.36 (0.035 %)	396.36 (0.035 %)
0.1	377.63	377.49 (0.037 %)	377.47 (0.042 %)
0.2	345.86	345.71 (0.043 %)	345.55 (0.090 %)
0.3	322.79	322.65 (0.043 %)	322.06 (0.226 %)
0.4	308.73	308.62 (0.035 %)	307.13 (0.518 %)
0.5	302.91	302.86 (0.016 %)	299.88 (1.000 %)
0.6	303.27	303.30 (0.010 %)	298.13 (1.695 %)
0.7	306.43	306.53 (0.033 %)	298.35 (2.637 %)
0.8	307.46	307.59 (0.042 %)	295.31 (3.952 %)
0.9	298.35	298.43 (0.027 %)	280.17 (6.094 %)
1.0	261.63	261.62 (0.003 %)	232.14 (11.27 %)

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.

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
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.

7 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 bending-extension 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.

REFERENCES

• Abramovich, H. (1998). Deflection control of laminated composite beams with piezoceramic layers-closed form solutions, Composite Structures 43: 217-231.
• Abramovich, H., Pletner, B. (1997). Actuation and sensing of piezolaminated sandwich type structures, Composite Structures 38: 17-27.
• Aldraihem, O.J., Khdeir, A.A. (2000). Smart beams with extension and thickness-shear piezoelectric actuators, Smart Materials and Structures, 9: 1-9.
• Beheshti-Aval, S.B., Lezgy-Nazargah, M. (2012). A coupled refined high-order global-local theory and finite element model for static electromechnaical response of smart maltilayered/sandwich beams, Archieve of Applied Mechanics 82: 1709-1752.
• Beheshti-Aval, S.B., Lezgy-Nazargah, M. (2013). Coupled refined layerwise theory for dynamic free and forced response of piezoelectric laminated composite and sandwich beams, Meccanica 48: 1479-1500.
• Bendary, I.M., Elshafei, M.A., Riad, A.M. (2010). Finite element model of smart beams with distributed piezoelectric actuators, J. Intelligent Material Systems and Structures 21: 747-758.
• Benjeddou, A., Trindade, M.A., Ohayon, R. (1997). A unified beam finite element model for extension and shear piezoelectric actuation mechanisms, J. Intelligent Material Systems and Structures 8: 1012-1025.
• Benjeddou, A., Trindade, M.A., Ohayon, R. (2000). Piezoelectric actuation mechanisms for intelligent sandwich structures, Smart Materials and Structures 9: 328-335.
• Carrera, E., Brischetto, S. (2008). Analysis of thickness locking in classical, refined and mixed multilayered plate theories, Composite Structures 82: 549-562.
• Chee, C.Y.K., Tong, L., Steven, G.P. (1999). A mixed model for composite beams with piezoelectric actuators and sensors, Smart Materials and Structures 8: 417-432.
• Crawley E.F., Anderson, E.H. (1990). Detailed models of piezoceramic actuation of beams, J. Intelligent Material Systems and Structures 1: 4-25.
• Crawley, E.F., de Luis, J. (1987). Use of piezoelectric actuators as elements of intelligent structures, AIAA J. 25: 1373-1385.
• Efunda.com, (2014). Efunda:Piezo Material Data, [online] Available at: Available at: http://www.efunda.com/materials/piezo/material_data/matdata_index.cfm [Accessed 22 Aug. 2014].
  » http://www.efunda.com/materials/piezo/material_data/matdata_index.cfm
• Jiang, J.P., Li, D.X. (2007). A new finite element model for piezothermoelastic composite beam, J. Sound and Vibration 306: 849-864.
• Kapuria, S., Hagedorn, P. (2007). Unified efficient layerwise theory for smart beams with segmented extension/shear mode, piezoelectric actuators and sensors, J. Mechanics of Materials and Structures 2: 1267-1298.
• Khdeir, A.A., Aldraihem, O.J. (2001). Deflection analysis of beams with extension and shear piezoelectric patches using discontinuity functions, Smart Materials and Structures 10: 212-220.
• Kumar, K.R., Narayanan, S. (2008). Active vibration control of beams with optimal placement of piezoelectric sensor/actuator pairs, Smart Materials and Structures 17: 55008-55022.
• Marinkovic, D., Marinkovic, Z. (2012). On FEM modeling of piezoelectric actuators and sensors for thin-walled structures, Smart Structures and Systems 9: 411-426.
• Narayanan, S., Balamurugan, V. (2003). Finite element modelling of piezolaminated smart structures for active vibration control with distributed sensors and actuators, J. Sound and Vibration 262: 529-562.
• Neto, M. A., Yu, W., Roy, S. (2009). Two finite elements for general composite beams with piezoelectric actuators and sensors, Finite Elements in Analysis and Design 45: 295-304.
• Peng, X.Q., Lam, K.Y., Liu, G.R. (1998). Active vibration control of composite beams with piezoelectrics: a finite element model with third order theory, J. Sound and Vibration 209: 635-650.
• Raja, S., Prathap, G., Sinha, P.K. (2002). Active vibration control of composite sandwich beams with piezoelectric extension-bending and shear actuators, Smart Materials and Structures 11: 63-71.
• Ray, M.C., Mallik, N. (2004). Active control of laminated composite beams using a piezoelectric fiber reinforced composite layer, Smart Materials and Structures 13: 146-152.
• Robbins, D.H., Reddy, J.N. (1991). Analysis of piezoelectrically actuated beams using a layer-wise displacement theory, Computers and Structures 41: 265-279.
• Shen, M.H.H. (1995). A new modeling technique for piezoelectrically actuated beams, Computers and Structures 57: 361-366.
• Sulbhewar, L.N., Ravenndranath, P. (2014a). A numerically accurate and efficient coupled polynomial field interpolation for Euler-Bernoulli piezoelectric beam finite element with induced potential effect, J. Intelligent Material Systems and Structures, DOI: 1045389X14544149, 26(12):1539-1550.
• Sulbhewar, L.N., Raveendranath, P. (2014b). An accurate novel coupled field Timoshenko piezoelectric beam finite element with induced potential effects, Latin American J. Solids and Structures 11: 1628-1650.
• Sulbhewar, L.N., Raveendranath, P. (2015). A locking-free coupled polynomial Timoshenko piezoelectric beam finite element, Engineering Computations, 32(5):1251-1274.
• Sun, B., Huang, D. (2000). Analytical vibration suppression analysis of composite beams with piezoelectric laminae, Smart Materials and Structures 9: 751-760.
• Wang, F., Tang, G.J., Li, D.K. (2007). Accuarate modeling of a piezoelectric composite beam, Smart Materials and Structures 16: 1595-1602.
• Zhang, X.D., Sun, C.T. (1996). Formulation of an adaptive sandwich beam, Smart Materials and Structures 5: 814-823.

Appendix: Coupled Shape Functions

Publication Dates

• Publication in this collection
  May 2016

History

• Received
  03 Dec 2014
• Reviewed
  04 Jan 2016
• Accepted
  16 Feb 2016