A beam element for piezoelectric planar frames under small strains but large rotations using the total Lagrangian description

* Corresponding author https://doi.org/10.1590/1679-78255963 Abstract Small strains are consistently incorporated into a large rotation piezoelectric beam theory, where the displacement is assumed to vary in accordance with the Timoshenko assumption and the electric potential has linear variation through each piezoelectric layer thickness. A finite element for planar frames, based on the total Lagrangian description, is then developed. The displacements and rotation are linearly approximated over the element and the voltage in a piezoelectric sensor layer is taken to be constant. In addition to the well-known membrane and shear lockings, it is disclosed the presence of locking in the application of Gauss law to the piezoelectric sensor layers. All the mechanical and electrical contributions to the internal load vectors are evaluated using a reduced one-point Gaussian quadrature to make this simple element efficient. No spurious modes are introduced in this process. An incremental-iterative approach, based on the Newton-Raphson algorithm, is employed in the solutions of the numerical examples to illustrate the large rotation capability


INTRODUCTION
Sensors and actuators attached or embedded into flexible host structures are two of the most basic elements of a smart structural system. In the sensor application, strains can be determined from measurement of the induced electrical potential, whereas in actuator applications strains can be controlled through the input of appropriate electric potential. Smart structures has been used in aerospace, hydrospace, nuclear, defense, automotive and civil structural applications. A vast literature is already available mainly focused on the linear analysis of beams, plates and shells, whereas less attention has been given to the more sophisticated nonlinear analysis (Cardoso and Fonseca, 2004;Chróścielewski et al., 2019). However, consideration of geometric nonlinearities can be crucial because external excitations can cause large rotations of smart structures, still in the range of small strains, due to their inherent flexibility (Zhang et al., 2015). A linear analysis may thus be insufficient to provide accurate estimations of the voltage sensed by piezoelectric sensors in light-weight and highly flexible structures (Zhang et al., 2015;Chaudhuri, 2002, 2005), which are important data for their control by actuators. It seems to be quite realistic to assume that such structures may be subjected to large rotations, but small strains. Moreover, it must be considered that the piezoceramics like the PZTs, which are the most commonly used piezoelectric material, are fragile and unable to bear large strains (Ramya et al., 2015).
In finite element analysis, the phenomenon of locking is quite well-known for linear formulations. Various types of locking, such as membrane locking, shear locking, volumetric locking, and so on, have been identified (Cook et al., 2002;Zienkiewicz et al, 2013). In particular, linear straight Timoshenko beam finite elements are likely to suffer only of shear locking. A number of such elements have appeared in the literature (Nickel and Secor, 1972;Tessler and Dong, 1981;Prathap and Bhashyam, 1982;Reddy, 2006). They differ from each other in the choice of interpolation functions used for the transverse displacement and rotation. The 2-node Timoshenko element with linear interpolation for the axial and transverse displacements, as well as for the rotation, is the simplest one. However, the shear locking makes this element very stiff in slender beams, i.e., as the length-to-thickness ratio becomes large. Detection of locking in a nonlinear formulation seems to have appeared for the first time in Kikuchi and Aizawa (1982). Lanzoni and Tarantino (2018) propose a three-dimensional displacement field to describe the plane bending of solids under large strains that takes into account anticlastic effect. The solid cross section is assumed to remain plane during deformation, with no restriction placed on the magnitude of its rotations. Some geometrical quantities of the displacement field lose their physical meaning for slender beams (Falope et al., 2019;Lanzoni and Tarantino 2020).
Lucena Neto et al. (2019) describe the behavior of Timoshenko beams under large rotations, for which small strains are consistently incorporated into a total Lagrangian formulation with the help of the polar decomposition of the deformation gradient (Malvern, 1969;Chandrasekharaiah and Debnath, 1994). In addition to shear locking, the use of the simple 2-node Timoshenko element in this highly nonlinear formulation also activates membrane locking. Both locking are alleviated by evaluating all the internal load vectors using a reduced one-point Gaussian quadrature.
Nonlinear models of piezoelectric beams subjected to large rotations with some restriction of magnitude (Mukherjee and Chaudhuri, 2002;2005;Kapuria and Alam, 2004) are more common in the literature than those with rotations free of any restrictions (Cardoso and Fonseca, 2004). In this paper, the work of Lucena Neto et al. (2019) is modified to account for piezoelectric patches attached to the beam surface, assuming electric potential with linear variation through each patch thickness (Marinković et al., 2009). In the finite element development, the voltage in a piezoelectric sensor patch is taken to be constant. This work reveals the existence of locking associated with the application of Gauss law to piezoelectric sensor layers. Numerical results demonstrate the accuracy of this simple element for nonlinear analysis of piezoelectric planar frames.

FUNDAMENTALS
In the development of the mathematical model, a beam with attached piezoelectric layers is considered. While the displacement field is assumed to vary in accordance with the Timoshenko assumption, the electric potential has linear variation through each piezoelectric layer thickness. Small strains are consistently incorporated.

Displacement, strain, electric potential and electric field
Suppose that a point ( , , ) of the beam shown in Figure 1 in the initial configuration moves to the point ( , , ) in the current configuration with no motion in the direction. Under the Timoshenko assumption (Reddy, Latin American Journal of Solids and Structures, 2020, 17(8 Thematic Section), e311 3/18 2002), the cross section is taken to be infinitely rigid in its own plane, remaining plane but not necessarily perpendicular to the beam axis during the motion. The points have their orthogonal Cartesian coordinates related by = + − sin = + cos = . (1) The displacement components ( ) and ( ) refer to the point ( , 0,0) on the beam axis, which moves to the point ( , 0,0), and ( ) is the cross-section rotation about the axis. The polar decomposition theorem (Malvern, 1969;Chandrasekharaiah and Debnath, 1994) ensures that the deformation gradient can be uniquely written in the form = where is orthogonal and fully responsible for the rigid body rotation, is symmetric and fully responsible for the strain and primes denote derivatives with respect to . If is the angle of rigid body rotation about the axis, and The symmetry of must be imposed to identify and complete the decomposition (3). However, the explicit value of to be derived from such an operation is not of interest here and, thus, will not be addressed.
One may define the displacement gradient associated with by means of to rewrite the Green strain tensor in the form Since the components of are necessarily small for small strains, the following simplified expression may be adopted in this case: Moreover, the smallness of the component sin( − ) of implies ≈ and thus leads to ϵ ≈ ϵ 0 − ′ ϵ ≈ 0 γ ≈ −(1 + ′ ) sin + ′ cos (9) where Figure 2 shows a host beam of length , having a rectangular cross section of width and thickness 2ℎ, with two piezoelectric layers bonded on its bottom and top surfaces. The lower layer has width 1 and thickness 1 , and the upper layer has width 3 and thickness 3 . In electrostatics, the electric field vector is irrotational and can therefore be derived from a scalar function ( , , ), known as electric potential, in accordance with The electric potential is here assumed to be independent of (beam with no motion in the direction) and the gradient / through the thickness of a piezoelectric layer is the only to be considered (Lammering, 1991;Almajid et al., 2001;Trindade et al., 2001;Sze et al., 2004): On the assumption that the electric potential varies linearly through the thickness of a piezoelectric layer , between the values at the bottom ( = ) and +1 at the top ( = + ) of the layer, The nonzero component of reduces to where � = +1 − is the voltage (difference of electric potential) in the piezoelectric layer . Marinković et al. (2007) and Sulbhewar and Raveendranath (2014) show that an electric potential consistent with a first-order shear deformation theory, as the one based on the Timoshenko assumption, should have an additional quadratic term in related to the curvature. Marinković et al. (2007;2009) also show that the use of (13) is sufficiently accurate for typically thin piezoelectric layers.
When a sensor (actuator) layer has the bottom surface grounded, the quantity stands for the electric potential read at (applied to) the sensor (actuator) top layer surface. On the other hand, when the sensor (actuator) layer has the top surface grounded, stands for the opposite electric potential read at (applied to) the sensor (actuator) bottom layer surface.

Principle of virtual displacements
The principle of virtual displacements for the piezoelectric beam, shown schematically in Figure 3, takes the form (Cardoso and Fonseca, 2004): The normal stress and transverse shear stress refer to the second Piola-Kirchhoff stress tensor; is the transverse electric displacement; and are the distributed axial and transverse forces (force per unit length in the initial configuration); , and are the end loads; , and are the end displacements; and is the piezoelectric beam cross-section area.

Constitutive relations
With polarization in the positive -axis and the remaining material principal directions aligned with the -andaxes, the following constitutive relations cover a wide range of piezoelectrics (Yang, 2018): Here 's and 's refer to the components of the second Piola-Kirchhoff stress tensor and electric displacement vector; ϵ's and ã's refer to the components of the Green strain tensor; and , and are the elastic, piezoelectric and dielectric coefficients. It is timely to say that the adopted measure for the electric displacement and electric field vectors, like the second Piola-Kirchhoff stress and Green strain tensors, refer to the initial configuration and are unaffected by rigid-body motion (Cardoso and Fonseca, 2004). This latter feature justifies the use of the relations (18), identical to those of the corresponding linear theory, as explained in Bathe (1996).
Special stress assumptions, commonly made in beam theory, require modifications of relations (18) to reduce them to the simple form where is the shear correction factor (Hutchinson, 2001;Puchegger et al., 2003). The coefficients � 11 , ̃3 1 and ̃3 depend on the type of special assumption made. Following Butz et al. (2008), Khdeir et al. (2012), Elshafei et al. (2014) and Sulbhewar and Raveendranath (2014) Substitution of (9), (14) and (19) into (17) gives where 11 = � 111 1 + � 11 + � 113 3 Here � 111 and 1 refer to the coefficient � 11 and the cross-section area of the bottom piezoelectric layer; � 11 and refer to the coefficient � 11 and the cross-section area of the host beam; and so on.

FINITE ELEMENT MODEL
In view of the constitutive relations for and in (22), expression (16) reads − � [( 11 ϵ 0 + 12 ′ + 14 � 1 + 15 � 3 ) ϵ 0 + ( 12 ϵ 0 + 22 Introducing the constant = 0 or = 1 in order to identify the piezoelectric layer as actuator or sensor, respectively, the above expression is rewritten as − � [( 11 ϵ 0 + 12 ′ + 1 14 � 1 + 3 15 � 3 ) ϵ 0 + ( 12 ϵ 0 + 22 ′ + 1 24 � 1 + 3 25 � 3 ) ′ 0 + ã + 1 1 � 1 + 3 3 � 3 ] The displacements and rotation are linearly approximated over the element, whereas the voltage is taken to be constant in a piezoelectric sensor layer: where collects the nodal values of , and , and the element constant values assumed for � 1 and � 3 . The shape functions are given by The adoption of constant voltages � 1 and � 3 is consistent with the fact that a piezoelectric layer contains metallic electrodes on its top and bottom surfaces so that each electrode forms an equipotential region (Yasin et al., 2010). Substitution of (28) into (25) The internal load vectors , , , 1 and 3 are associated with the axial force , bending moment , shear force , and the electrical quantities 1 and 3 , respectively. The external load vector is work equivalent to the distributed forces and , and the voltages � 1 and � 3 applied to the piezoelectric actuator layers. The unknown reactions are denoted by . As (31) must hold for arbitrary , the term in parentheses must vanhish to result in the eight element equations ( ) = + + + 1 + 3 + + = .
The first six equations represent the element nodal equilibrium equations, while the last two equations are the application of Gauss law to the piezoelectric layers.

COMPUTATIONAL ASPECTS
If does not satisfy the system of nonlinear equations (33) for a given externally applied load or voltage, i.e., ( ) ≠ , corrections Δ is found using the Newton-Raphson algorithm (Zienkiewicz et al, 2013) replacing (33) by the linear system where Δ = ( ) is the unbalanced load vector. In this iterative process, the element tangent stiffness matrix is updated at each iteration, replacing by + Δ , until a state of convergence is achieved.
The linear system (34) of each element must be written first in the global system to then obtain the assembled linear system of the entire structure from all element contributions. In the notation of Figure 4, the increments Δ and Δ in the local system are related to the increments Δ and Δ in the global system by The fact that is orthogonal implies that the element tangent stiffness matrix in both systems are related by = . (38) The contribution of ϵ to the vector is clearly of follower nature because it is dependent on the structure displacements (Hibbitt, 1979). Marinković et al. (2009), Rama et al. (2016) and Marinković and Rama (2017) devote special attention in the analysis of similar contribution present in their models of piezoelectric shells. One must point out here,however, that the work done by the contribution of ϵ is not path dependent because it does not affect the symmetry of and, consequently, it is definitely conservative (Bažant and Cedolin, 2010).
Because the magnitudes of the coefficients , and are quite different, the condition number of the global tangent stiffness matrix can be very high. To overcome this problem, Qi et al. (1997) proposed to balance the magnitudes of the constitutive coefficients where a multiple of the force unity 1 N ′ ← 10 N is used with ∈ . Doing this, [N/m 2 ] and ξ [C 2 /Nm 2 ] in (18) are changed to ′ ← × 10 − and ξ ′ ← ξ × 10 . Also, the voltage � [Nm/C] is changed to � ′ ← � × 10 − . For the PZT-5H given in Table 1, the value = 9 is chosen so that the orders of 10 1 exhibited by ′ and 10 0 exhibited by ξ ′ are well balanced with the order of 10 0 exhibited by .
The internal load vectors and may be rewritten as are mechanical contributions and are electrical contributions from piezoelectric sensor layers. The magnitudes of and are very small in comparison with those of and , respectively. As reported by Laulusa and Reddy (2004), locking generates large internal load vectors that in turn may stiffen the structure. Lucena Neto et al. (2019) show the importance of evaluating and using a reduced one-point Gaussian quadrature to alleviate membrane and shear locking in purely mechanical beams. The contributions of , and are also here evaluated in this manner. While is always exactly obtained with just one Gauss point, would require more points for nonsymmetric cross sections. Reduced one-point Gaussian quadrature is also important in the evaluation of 1 and 3 because of locking, and in the evaluation of the contribution of ϵ to the vector because of computational efficiency. The use of reduced integration may introduce spurious modes into the element stiffness matrix (Cook et al., 2002;Zienkiewicz et al., 2013). Fortunately, an eigenvalue analysis of the stiffness matrix (35) in the initial configuration and obtained with just one Gauss point reveals that the element possesses only three null eigenvalues associated with the expected physical rigid body motions.

RESULTS
According to the mathematical model, the piezoelectric layers solely operate in extension mode since both polarization and electric field are in the -direction. In the numerical evaluation of the constitutive parameters � 11 , ̃3 1 and ̃3 in (20), the host structure is assumed to be made of aluminum with Young´s modulus = 70.3 GPa and Poisson´s ratio = 0.345, and the piezoelectric layers are assumed to be made of PZT-5H with properties indicated in Table 1 (Yang, 2018). As a transversely isotropic material (isotropic in the plane normal to -axis), the PZT-5H presents 23 = 13 , 22 = 11 , 55 = 44 , 66 = ( 11 − 12 )/2, 32 = 31 , 24 = 15 and 2 = 1 .
The accuracy and effectiveness of the proposed element under mechanical or electrical loading is demonstrated by means of nonlinear static problems. A code in MATLAB language has been written to carry out the numerical tests. For convergence control, displacement and unbalanced load criteria are selected (Cook et al., 2002). The Newton-Raphson iterations cease when The quantity represents the ratio between the Euclidean norms of the displacement increment vector ∆ at the end of iteration and the accumulated displacement vector +1 . The quantity represents the ratio between the Euclidean norms of the unbalanced load vector ∆ at the end of iteration and the sum of the equivalent external load vector and the load vector applied directly to the nodes. Table 2 Values of , , Θ, Φ for the cantilever beam with transverse tip force and distributed sensing according to Santos (2018).

Cantilever beam with transverse tip force and distributed sensing
The cantilever beam of Figure 5 is subjected to a transverse tip force, with two symmetrical piezoelectric layers continuously distributed as sensors. The following geometry is adopted (units in mm): Table 2 shows the solutions s , s and Θ s found in Santos (2018) For γ = 0 (no transverse shear strain) and ϵ 0 = 0 (inextensibility of the beam axis), the Santos' solutions are "exact" in the sense that they are analytically stated in a closed form. Moreover, these solutions can be widely used for a fixed value of 2 / 3 . In the sequel, the performance of the proposed element is accessed in terms of the parameters Lucena Neto et al. (2019) have shown that a reduced one-point Gaussian quadrature must be employed in the evaluation of mechanical contributions to the internal load vector to properly alleviate the membrane and shear locking effects. One should now investigate the importance of a reduced quadrature in the evaluation of 1 and 3 in order to detect the presence of locking in the application of Gauss law to the piezoelectric sensors.
Taking 2 / 3 = 4 at 10 increments and a coarse regular mesh of 4 elements, Figure 6 plots 1 and 3 against the length-to-thickness ratio / , where = 2ℎ + 1 + 3 , varying the number of Gauss points in the evaluation of 1 and 3 . Results of , , and are not shown here as they are indistinguishable from those of Santos (2018) at the scale of the figure. Whenever the quadrature rule includes two or more points, the locking effect is clearly triggered, having almost the same intensity for a given ratio / and being particularly pronounced for slender beams. A reduced one-point Gaussian quadrature alleviates the locking effects regardless of the value of / . Therefore, to make this simple element efficient, just one Gauss point will be adopted from now on in the numerical integrations of the internal load vectors.
A convergence analysis with respect to the mesh refinement is now carried out. Regular meshes of 4, 8, 16, 32 and 64 elements and tip forces of 2 / 3 = 1, 4, 10, each applied at 10 increments, are considered. Table 3 shows that the obtained results converge quickly to the solutions given by Santos (2018). It should be noted that: (a) the tip displacements and rotation converge faster than the voltages at the clamped end of the beam. In fact, the voltages � 1 (0) and � 3 (0) given by Santos (2018) are compared with the constant values measured over the sensor patches attached to the element placed next to the clamped end; (b) the convergences of the displacements and sensed   voltage are monotonic from below, while the convergence of the rotation is monotonic from above; (c) the induced voltage measured by the bottom sensor patch is, in general, slightly smaller than that measured by the top one. This can be physically understood if one imagines that the contraction of the bottom patch is also, in general, slightly smaller than the extension of the top patch.
Although the table results are obtained using / = 50, nearly no changes are observed for / = 400 which reveals the effective locking alleviation given by the one-point reduced integration scheme. However, when the code is running for / = 4 (see tabled results in parentheses) larger differences are observed. This is a clear consequence of the relevant contribution given by the transverse shear strain. It is worth mentioning that for all ratios / the obtained results are in close agreement with those found by Meireles (2019), which uses finite element corotational models based on a first-order shear deformation theory.
Deformed shapes of the cantilever beam at successive load levels 2 / 3 = 1, 4, 10 are depicted in Figure 7. As shown, excellent results are quickly reached with the mesh refinement. This same comment also applies to the results ploted in Figure 8, where the distributions of � 1 and � 3 along the beam are obtained for 2 / 3 = 10 using a mesh of 8 elements.

Cantilever beam with distributed actuation
With respect to the previous example, the following changes are carried out: (a) the applied tip load is removed; (b) the electric potentials � 1 = 10 and � 3 = 0 (units in volt) are applied to the bottom and top piezoelectric layers,respectively; (c) the thickness of the bottom piezoelectric layer is set at a higher value of 1 = 2.5 mm. Santos (2018) has reported the "exact" results s = −5.1638 × 10 −8 mm s = −5.9929 × 10 −6 mm Θ s = −5.9929 × 10 −5 rad (48) for this problem.
Changes in the number of Gauss points in the evaluation of the contribution of ϵ to the vector have led to nearly the same results. Adopting for simplicity a reduced one-point Gaussian quadrature in that evaluation, the convergence in terms of , and is attained, with at least 3 significant figures, even for a coarse mesh of 2 elements and 2 load increments. This remarkably result with a reduced number of increments is mainly due to the properly account of ⁄ in the tangent stiffness matrix, whose approach is apparently not followed by Marinković et al. (2009), Rama et al. (2016) and Marinković and Rama (2017) in their nonlinear analyses of piezoelectric shells. Figure 8: Distribution of the sensed voltage Φ 1 = 44 � 1 / 24 1 and Φ 3 = 55 � 3 / 25 3 for 2 / 3 = 10 using 8 elements.

Modified Williams toggle frame
The toggle frame of Williams (1964) has been used by several authors (Wood and Zienkiewicz, 1977;Nanakorn and Vu, 2006;Lee et al., 2014) to verify the performance and efficiency of nonlinear finite elements. The nonlinear response of the frame is dictated by the geometry and inclination of the members. The original version of the toggle frame involves clamped supports. Herein, in order to intensify the nonlinear effects, one has modified the original geometry and changed the support conditions to simply-supported ones (see Figure 9). Moreover, two piezoelectric sensor patches, with 25 mm of width and 1 mm of thichness, are bonded to the top surface of each frame leg which has 25 mm of width and 2 mm of thichness. Due to the presence of snap through, the vertical displacement ∆= 6 mm applied at 10 increments, rather than the load , is prescribed at the apex (displacement control). In a more sophisticated manner, the arc-length technique could also be used to trace the whole equilibrium path in a slightly expensive way (Crisfield, 1991).
Assuming that the structure deformation takes place symmetrically, only half of the frame needs to be meshed. Figure 10 shows the responses (∆), obtained with a regular mesh of eight elements per frame leg, for the structure with and without sensor patches. For the purely mechanical case, converged results obtained from a Nastran model of CBEAM elements are also included in the figure (marked with small circle), which are in very good agreement with the developed finite element. Finally, the voltages � 3 predicted by the sensor patches 1 and 2 along the loading history are plotted in Figure 11. Both patches exhibit highly nonlinear responses which can be useful to access the limit point of the structure. It is worth mentioning that the presence of turning points in the space � 3 , as shown in the figure, does not bring any numerical difficulties because the solution control is guided by the vertical displacement ∆ at the apex (see space ∆ in Figure 10).