Plane Thermoelastic Waves in Infinite Half-Space Caused FACTA UNIVERSITATIS Series: Mechanical Engineering Vol. 17, No 2, 2019, pp. 269 - 283 https://doi.org/10.22190/FUME190530030M © 2019 by University of Niš, Serbia | Creative Commons License: CC BY-NC-ND Original scientific paper ABAQUS IMPLEMENTATION OF A COROTATIONAL PIEZOELECTRIC 3-NODE SHELL ELEMENT WITH DRILLING DEGREE OF FREEDOM Dragan Marinković, Gil Rama, Manfred Zehn Department of Structural Mechanics, Berlin Institute of Technology, Germany Abstract. Integration of classical, passive structures and active elements based on multi- functional materials resulted in a novel structural concept denoted as active structures. The new structural systems are characterized by self-sensing and actuation. Coupling the two distinctive features by means of a controller enables a number of exquisite functionalities such as vibration suppression, noise attenuation, shape control, structural health monitoring, etc. Reliable, accurate and highly efficient modeling tools are an important ingredient of the active structure design. This paper addresses the Abaqus implementation of a recently developed piezoelectric 3-node shell element. The element uses co-rotational formulation to cover geometric nonlinearities. Special techniques are used to address the issues originating from low-order interpolation functions. The discrete shear gap is used to resolve the shear locking, while the assumed natural deviatoric strain technique improves the membrane behavior. Examples are computed in Abaqus upon implementation of the developed element. Key Words: Abaqus, Corotational FEM, Piezoelectric shell element, Discrete shear gap, Assumed natural deviatoric strain 1. INTRODUCTION Over two decades ago, a novel structural concept denoted as ‘active’, ‘adaptive’, or ‘smart’ structures has seen the light of day [1]. It features integration of classical, passive structures and active elements based on multi-functional materials. In this manner, artificial systems are given the ability of self-sensing and actuation coupled by control capabilities. The concept is obviously the result of mimicking the natural systems that Received May 30, 2019 / Accepted July 15, 2019 Corresponding author: Dragan Marinkovic Department of Structural Mechanics, TU Berlin, Str. d. 17. Juni 135, 10623 Berlin E-mail: dragan.marinkovic@tu-berlin.de 270 D. MARINKOVIĆ, G. RAMA, M. ZEHN actively react to environmental stimuli in order to protect their integrity and maintain optimal functionality. This provides them with a number of exquisite functionalities such as vibration suppression [2], noise attenuation [3], structural health monitoring [4], shape control [5], etc. In this manner, active systems offer significant benefits over passive ones, including improved safety, ergonomics, robustness, product lifespan, etc. Such an enticing research field has attracted many researchers, whereby, generally speaking, two major directions of work can be distinguished. One group of researchers focused their work on resolving some practical problems such as design of active elements, types of materials applied, ways of embedding active elements into passive structures, etc. The other group dedicated their work to development of modeling tools thus enabling reliable simulation and therewith less expensive and successful design. The Finite Element Method (FEM), as the most powerful method in the field of structural analysis, is addressed and many authors proposed various types of elements that enable modeling of active structures. In case of thin-walled structures, active elements are typically in the form of small, rather thin patches made of piezoelectric materials that operate based on the e31-effect. This actually means that the patches couple the electric field acting across the thickness to the in-plane strains. The same patches are used as both sensors and actuators with the essential difference only in the boundary conditions. When exposed to mechanical deformation, sensor patches deliver electric signals, i.e. voltages, proportional to the average in-plane strains in the area covered by the piezoelectric patch. Oppositely, when supplied with electric voltage, actuator patches produce distributed in-plane forces proportional to the electric voltage and acting perpendicularly to the edges of the patch. The patches are usually placed at the maximum offset distance from the mid-plane in order to produce the maximum bending moment uniformly distributed over the patch edges. They are commonly glued onto the outer shell surfaces. The effect may be amplified by using a pair of such patches collocated with respect to the structure’s mid- surface by placing one of them on the upper while the other one on the lower surface, and they are additionally oppositely polarized in order to produce only a bending moment when exposed to the same electric voltage. A large number of shell type of finite elements were proposed for piezoelectric active shell structures. Benjeddou [6] gave a thorough survey of the development in the 90’s and the development in the beginning of the new millennium retained the same intensity and innovativeness. Various approaches to modeling were investigated. One should notice that the application of piezopatches onto a passive shell structure results in a multilayered structure, even if the passive structure consists of a single layer of material. However, not rarely the passive structure is actually made of a composite laminate. The mainstream developments implemented the equivalent single-layer approach. Some of them were based on the classical laminate theory (Kirchhoff-Love kinematics) [7, 8], while others used the first-order shear deformation theory (Mindlin-Reissner kinematics) [9, 10]. The latter was more frequently used. One reason for this choice is to be sought in the greater generality as the formulation includes the transverse shear effects, but an equally important reason is the attractiveness of the C0-continuity needed from the element shape functions, whereas the Kirchhoff-Love elements demand the C1-continuity. To improve the accuracy of shell formulations at the cost of somewhat higher numerical effort, layerwise theories were also addressed. One of the elements based on this approach is the Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 271 9-node plate element proposed by Cinefra et al. [11] and which implements the technique of Mixed Interpolation of Tensor Components and variable through-the-thickness layer- wise kinematics. The approach was recently extended by Carrera et al. [12, 13] in order to locally increase the accuracy by means of node-dependent kinematics. Furthermore, both geometrically [14, 15, 16] and materially [17] nonlinear problems in the behavior of piezoelectric thin-walled structures were also considered in the work of many researchers. As it is of utmost importance to enable users to apply developed numerical tools in modeling and simulation, this paper addresses Abaqus implementation of the recently developed piezoelectric 3-node shell element. Similar works were already reported in available literature [18, 19]. In what follows, the most important aspects of the element development are briefly presented together with the results of several test cases computed by using the element implemented in Abaqus. 2. PIEZOELECTRIC 3-NODE SHELL ELEMENT WITH DRILLING DEGREE OF FREEDOM The developed element combines features of already developed and in literature available shell elements. It implements the mechanical field of the element proposed by Rama et al. [20] that was also extended to composite laminates by Rama et al. [21], and the electric field, including the both way coupling between the two, as described in the work by Marinkovic and Rama [22]. As the mentioned references provide detailed derivations, the element description given in this section will not go into all the details, but for the sake of completeness, some major aspects of the development will be presented together with the most important matrices. As is the case with most shell elements, the formulation of this one strongly relies on the use of local coordinate frame x, y, z, which is defined so as to have one of its axis, the z-axis, perpendicular to the element, while the other two axes lie in the element’s plain, Fig. 1. In the case of an isotropic material, their orientation is mainly of importance for proper evaluation of the results, as stresses and strains are determined and typically given with respect to this coordinate system. On the other hand, in the case of a composite laminate, it is of crucial importance for defining the material elastic properties. The same coordinate system comes also quite handy in the presence of piezoelectric layers, as the mathematical description of the considered e31-piezoelectric effect is straightforward with respect to it. Fig. 1 3-node shell element with the local frame and mechanical degrees of freedom 272 D. MARINKOVIĆ, G. RAMA, M. ZEHN Since the 3-node element is a flat element, its mechanical field is effectively a superposition of a plate and a membrane element, Fig. 2. Fig. 2 3-node shell element as a superposition of plate and membrane elements According to the Mindlin-Reissner kinematics, the displacement field with respect to the local coordinate frame is given as follows: 3 3 1 12 0 i i yi i i i i x i i i u u h v N v N t w w                                        , (1) where ui, vi, and wi , i=13, are the nodal displacements along the x-, y- and z-axis, respectively, xi and yi, are the nodal rotations about the x- and y-axis, respectively, Ni are the typical nodal shape functions of a linear triangular element, hi are the nodal values of shell thickness and t is the natural coordinate in the thickness direction. The plate behavior of the element is determined by plate stiffness and transverse shear stiffness. The plate stiffness is obtained by deriving the bending strains directly from the displacement field. The resulting bending strain-displacement matrix reads [20]: 23 31 12 32 13 21 23 32 31 13 12 21 0 0 0 0 0 0 1 0 0 0 0 0 0 2 0 0 0 pb y y y B x x x A y x y x y x                             , (2) where A is the element area, while xij= xi- xj and yij= yi- yj. Since shell elements are notorious for shear locking effects, particularly when low- order shape functions are used, special measures are typically needed to mitigate the effect. The developed element implements the discrete shear gap technique originally proposed by Bletzinger et al. [23] and improved by the cell strain smoothing technique suggested by Ngyen-Thoi et al. [24]. The resulting transverse shear strain-displacement matrix for a triangular element is given by [20]: 32 13 1 3 21 2 3 23 31 4 2 12 4 1 01 02 s x A x a a x a a B y A y a a y a aA                  , (3) where: Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 273 1 12 13 2 31 21 3 21 13 4 12 31 1 1 1 1 , , , 2 2 2 2 a y x a y x a x x a y y           . (4) In order to avoid the dependence of the strain-displacement matrix on the node numbering and to improve the accuracy, the strain smoothing technique [24] is applied. It implies that the element is divided into three sub-elements by an additional node at the element centroid. Equation (3) is then applied to each of the sub-elements and, at the same time, the assumption is used that the displacement at the element centroid is given as the average value of the displacements at the original three nodes. In this manner, the central node is condensed out. Finally, the resulting transverse shear strain-displacement matrix is simply given by averaging the matrices of all three sub-elements, [Bs i] : 3 1 1 3 i ps s i B B           . (5) Since the strain-displacement matrices are constant over the element area, the plate stiffness is obtained in the following manner:      T T p pb ps pb pb ps ps K K K A B D B B F B                               , (6) where [D] and [F] are the laminate bending and transverse shear stiffness matrices obtained by integrating the corresponding constants from the Hooke’s matrix across the thickness. A particular weakness of low-order 3-node shell elements resides in their membrane behavior. ANDES formulation by Felipa and Militello [25] aims to resolve this weakness and to improve the behavior of the element to the level of at least a quadratic element. It represents a combination of the free formulation (FF) proposed by Bergan and Nygard [26] and a modification of the assumed natural strain (ANS) formulation derived by Park and Stanley [27]. The formulation abandons the discretized displacement field (Eq. (1)) and instead, represents the displacements as a superposition of carefully chosen linearly independent modes consisting of rigid-body, constant-strain and linear-strain modes. The first two groups are denoted as basic modes, while the third group as higher-order modes. Accordingly, the stiffness is divided into basic and higher-order stiffness, each determining the behavior of the corresponding modes. For a unique transformation between the modal and nodal degrees of freedom, the number of modes is the same as the number of nodal degrees of freedom. The so-called drilling degree of freedom, i.e. the rotation around the local thickness axis (z-axis) plays the crucial role in the definition of the modes. On the other hand, one should notice that this degree of freedom is not a part of the discretized displacement field given by Eq. (1). Already this important difference speaks in favor of the new quality offered by ANDES formulation. The basic membrane stiffness is the same as in the FF formulation [26]:       T / ( ) mb K L A L Ah , (7) where [A] is the membrane stiffness of a laminate, i.e. a part of the ABD matrix, while [L] reads: 274 D. MARINKOVIĆ, G. RAMA, M. ZEHN   23 32 32 23 23 32 12 32 31 12 31 13 12 21 31 13 13 31 31 32 12 13 12 23 12 21 23 32 12 21 21 12 12 32 12 21 0 0 ( ) / 6 ( ) / 6 ( ) / 3 0 0 2 ( ) / 6 ( ) / 6 ( ) / 3 0 0 ( ) / 6 ( y x x y y y y x x x x y x y y x h L x y y y y x x x x y x y y x x y y y y x x                                                  23 31 23 32 31 13 ) / 6 ( ) / 3x x y x y                                   , (8) with  denoting a non-dimensional parameter. The higher-order modes cover the in-plane bending, as illustrated by the three modes on the left in Fig. 3. The linear dependency of those three modes is the reason to introduce the fourth, torsional mode (the same rotation around the z-axis at all three nodes) shown on the right in Fig. 3. Fig. 3 Higher-order modes – three bending modes (left) and a torsional mode (right) Without interpretation of single matrices that appear below, the higher-order membrane stiffness is given by [20]: T 0 9 4 mh u u K T K T              , (9) where 0 is a non-dimensional parameter, and furthermore [20]: 32 32 13 13 21 21 32 32 13 13 21 21 32 32 13 13 21 21 4 0 0 1 0 4 0 4 0 0 4 u x y A x y x y T x y x y A x y A x y x y x y A                             , (10)      T TT 4 4 5 5 6 6 1 3 nat nat nat K Q A Q Q A Q Q A Q                              , (11)   T nat nat nat A T A T          , (12) Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 275       2 2 2 23 13 21 23 13 21 23 31 32 13 21 2 2 2 31 21 32 31 21 32 31 12 13 21 322 2 2 2 12 32 13 12 32 13 12 23 21 32 13 1 4 nat y y l x x l y x x y l T y y l x x l y x x y l A y y l x x l y x x y l                                    , (13) with lij denoting the length of the element edge between nodes i and j. Furthermore:             4 1 2 5 2 3 6 1 3 1 1 1 , , 2 2 2 Q Q Q Q Q Q Q Q Q                    , (14)     3 9 7 9 5 61 2 4 2 2 2 2 2 2 2 2 2 21 21 21 21 21 21 21 21 21 5 6 3 8 94 1 2 1 2 32 2 2 2 2 2 2 32 32 32 32 32 32 32 3 7 8 9 6 54 2 2 2 2 2 2 13 13 13 13 13 13 2 2 2 , , 3 3 3 l l l l l l l l l A A A Q Q Q l l l l l l l l l l l l l l                                                                 7 2 2 2 32 32 1 2 2 2 13 13 13 l l l l                       . (15) There are different possibilities for the choice of parameters  and 0-9. This aspect is elaborated in [21] and the following choice explained: =1/8, 0= 2/4, 1=1, 2=2, 3=1, 4=0, 5=1, 6= –1, 7= –1, 8= –1, 9= –2. Upon the very compact description of the mechanical field in the element, the attention needs to be turned to the electric field in the piezoelectric layers. The function describing the electric field distribution across the thickness of piezolayers is supposed to be consistent with the Maxwell’s equations for dielectrics and with the element kinematics. In the case of Mindlin-Reissner kinematics, the consistent electric field is linear across the thickness [28]. However, it has also been shown [28] that the classical assumption of constant electric field over the thickness provides sufficient accuracy for typical, rather thin piezopatches. The electric field of the kth piezolayer is then given by: k k k E h   , (16) where k is the difference of electric potentials between the electrodes of the k th piezolayer and hk denotes the thickness of the piezolayer, leading to the typical diagonal form of the electric field – electric potential matrix:   1 k B h               , (17) The dielectric stiffness matrix is then defined by:       T V K B d B dV          , (18) 276 D. MARINKOVIĆ, G. RAMA, M. ZEHN where [d] is the matrix of dielectric constants at constant strain. Finally, the piezoelectric stiffness matrix reads:    TT u u mf V K K B e B dV               , (19) where [Bmf] collects the strain-displacement terms that define membrane and flexural strains, because the formulation considers the e31-piezoelectric effects that couples the in- plane strains with the thickness-oriented electric field. In Eq. (19) [e] is the matrix of piezoelectric constants. The integration in Eqs. (18) and (19) needs to cover all the piezoelectric layers across the thickness of the structure. In order to cover geometric nonlinearities, the element implements corotational FE formulation [29]. Beside the updated and total Lagrangian formulations [30] that both represent standard solutions in major commercially available FE codes, the corotational FE formulation has gained lately in importance in many applications ranging from real- time simulations [31] to engineering solutions [32]. All the details are available in the literature (including above mentioned references), and only the basic idea is briefly outlined here. The general idea is to introduce a new, corotational coordinate frame attached to the element that performs the same rigid-body motion as the element itself and to measure all necessary quantities, such as strains and stresses, with respect to it (Fig. 4). With the presented shell element, the local frame is used as corotational. Small strains are assumed despite finite displacements and rotations. Observing the problem in this way allows for further simplifications, such as to use linear dependency between deformational displacement (i.e. with the rigid-body motion removed from the overall displacements) and strains. With the rigid-body motion extracted from the overall motion, one may further determine all other mechanical and electric quantities in a relatively straightforward manner. For the sake of brevity, the equations are not presented here, but an interested reader may consult the above references for more details. Fig. 4 The idea behind the corotational FE formulation Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 277 3. NUMERICAL EXAMPLES The element was implemented in Abaqus by means of User Subroutine (UEL) and currently supports up to two electric degrees of freedom per element, i.e. two piezolayers that can be used both as sensors and actuators. The below given examples are of academic nature. As the focus is on numerical results, all the values are given as dimensionless. Of course, any set of consistent units corresponding to the quantities (both mechanical and electric quantities and material parameters, including dielectric and piezoelectric parameters) can be associated with them. 3.1. Bimorph beam Bimorph beam is one of the classical benchmark cases to test the developed numerical tools for piezoelectric active structures. The entire beam structure consists of two piezoelectric layers with opposite polarization. As explained in the introduction, when both layers are exposed to the same electric voltage, bending moment uniformly distributed over the edges of the structure is induced, thus causing bending. The left end of the structure is clamped (Fig. 5) and the displacement along the beam length is observed as a representative solution. The length of the beam is l=0.1, the width is w=0.005 and the overall thickness (both layers together) h=0.001. The Young’s modulus of the material is 2109, piezoelectric constant of the e31-coupling is 0.046 and the dielectric constant 0.106210 -9. The difference of electric voltages supplied to the electrodes placed on the outer surfaces of the beam is 1. Fig. 5 Bimorph beam The analytical solution for deflection [33] is obtained by assuming beam kinematics (implying Poisson’s ratio equal to 0) and it is a quadratic function in x. Fig. 6 gives a screen-shot of the deformed configuration with the contour plot of displacements obtained in Abaqus. The same case is also considered as a sensor case by exposing the structure to two transverse forces of magnitude 0.1, acting at the two free beam corners. It is assumed that the whole beam acts as a single sensor (each face completely covered by a single electrode) and the resulting electric voltage reflects the average strain in the whole beam. The obtained result for the electric voltage of 165 coincides with the analytical solution. A sensitivity analysis of the electric voltage to mesh distortion is performed. In order to summarize the results of both cases, Fig. 7a shows the diagram of deflection normalized with respect to the analytical value of the displacement at the beam tip, while Fig. 7b 278 D. MARINKOVIĆ, G. RAMA, M. ZEHN represents the voltage sensitivity to mesh distortion, by normalizing it in the same manner, i.e. with respect to the analytical value of 165. The distortion is performed by increasing the value of parameter e from 0 to 20 with the increment size of 5 (see FE mesh within the diagram in Fig. 7b). Fig. 6 Deformed bimorph beam – Abaqus contour plot for deflection Fig. 7 Bimorph beam: a) Actuator case – normalized deflection; b) Sensor case – normalized sensor voltage 3.2. Piezolaminated arch A curved structure is considered in this example. It is a semi-cylindrical arch with dimensions and boundary conditions as depicted in Fig. 7. The single layer of passive material has thickness of 5.842, the Young’s modulus is 68.95103 and the Poisson’s ratio 0.3. Each of the two outer piezolayers has the thickness of 0.254 and the following material properties: the Young’s modulus 63103 and Poisson’s ratio 0.3. Furthermore, the same piezoelectric constant e31=16.1110 -6 is assumed in all in-plane directions, and the dielectric constant is 1.6510-11. A vertical force of 200, acting upwards at the free tip, deforms the structure. The free tip displacements in the x- and y- directions and the induced sensor voltage in the inner piezolayer are observed in a geometrically nonlinear analysis performed in Abaqus. a) b) Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 279 Fig. 8 Semi-cylindrical arch – dimensions and boundary conditions The obtained results are compared with those reported by Zhang [34], who originally proposed the example. Fig. 9 depicts the deformed configuration computed using the developed 3-node shell element implemented in Abaqus. The contour plot in Fig. 9, left, corresponds to the displacement in the global x-direction (U1), while the one in Fig. 9, right, to the displacement in the global y-direction (U2). Fig. 10 shows the time history of displacements in the x- and y-directions with the increasing force, computed by the implemented 3-node shell element, Abaqus 3-node shell element and the results by Zhang [34]. Fig. 9 Abaqus contour plots for displacements in x- (left) and y-direction (right) Fig. 10 Displacements in x- and y-directions versus force 280 D. MARINKOVIĆ, G. RAMA, M. ZEHN Fig. 11 compares the results obtained for the sensor voltage by the present element with the results reported by Zhang [34]. Fig. 11 Sensor voltage versus force 3.3. Modal analysis of a piezolaminated composite plate under various electric boundary conditions The third example considers the change of natural frequency of a piezolaminated composite plate with the change of electric boundary conditions. The modal analysis is performed with short-circuited and open electrodes. Short-circuited electrodes imply zero voltage so that the behavior of the structure is purely mechanical (as if there was no piezoelectric effect present). Oppositely, with the electrodes left open, electric voltage is induced as a consequence of deformation, which further gives rise to mechanical stresses and, therewith, changes in the natural frequencies of the structure. The considered structure is a composite plate with the sequence of layers [p/0/90/0/p] with respect to the global x-axis (Fig. 11), where ‘p’ stands for piezoelectric layer, while the composite layers are made of graphite-epoxy and have the following mechanical properties with respect to the principal material orientations: Young’s moduli E11=1.32410 5 and E22=1.0810 4, Poisson’s ratio 12=0.33 and shear modulus G23=6.610 3. The piezoelectric layers have the following mechanical properties: E11=8.1310 4 and E22=8.1310 4, Poisson’s ratio 12=0.33 and shear modulus G23=2.5610 4. The piezoelectric constant is e31=14.810 -6 and the dielectric constant 1.150510-11. Fig. 12 Piezolaminated plate with: short-circuited (left) and open electrodes (right) Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 281 The span of the plate is a=200. The thickness of each composite layer is 1.068 and of each piezolayer 0.4. The modal analysis is performed for all edges simply supported. To demonstrate the influence of electric boundary conditions onto the resulting natural frequency, only the first eigenmode of the plate is observed. It is depicted in Fig. 13. Fig. 13 The first eigenmode of the piezolaminated composite plate For the purely mechanical case, the results by Abaqus linear shell element (S3) are also provided, while the result obtained using a fine mesh of Abaqus quadratic shell elements (S8) is used as a reference solution. This example was originally proposed by Saravanos et al. [35]. Hence, their results are used as a reference for the open circuit case that is affected by the piezoelectric coupling. All the results are summarized in Table 1. Table 1 Convergence analysis – normalized first natural frequency Short-circuited fref = 22915 Hz (Ref: Abaqus S8 2424 mesh) Open circuit fref = 24594 Hz (Ref: Saravanos et al. [35]) Elements Present Abaqus S3 [35] Present [35] 32 1.190 1.220 1.090 1.107 1.109 128 1.031 1.050 1.034 1.028 1.054 288 1.008 1.024 1.023 1.005 1.044 4. CONCLUSIONS Development of reliable, accurate and highly efficient numerical tools for modeling and simulation of thin-walled piezoelectric active structures is an important prerequisite for successful design of those modern structural systems. The presented shell element combines features of already available elements. Since it is a 3-node element, it offers high meshing ability and numerical efficiency. Geometric nonlinearities are accounted for based on the corotational FE formulation. The major weaknesses typical for low-order interpolation elements are addressed by means of various existing techniques in order to offer a versatile shell type finite element that also covers the electro-mechanical coupling. In order to bring the element closer to the end user, it was implemented in Abaqus. Several test cases were considered to demonstrate applicability of the element to both sensor and actuator cases in linear and geometrically nonlinear analysis. 282 D. MARINKOVIĆ, G. RAMA, M. ZEHN In further work, the element should be extended to cover directional in-plane polarization. This would allow to model composites with piezoelectric fibers. A further extension to cover nonlinearities in the piezoelectric coupling would allow to accurately model the cases that involve strong electric fields. REFERENCES 1. Gandhi, M.V., Thompson, B.S., 1992, Smart Materials and Structures, Chapman & Hall, London. 2. Bendine, K., Satla, Z., Boukhoulda, F.B., Nouari, M., 2018, Active vibration damping of smart composite beams based on system identification technique, Curved and Layered Structures, 5(1), pp. 43-48. 3. Gabbert, U., Duvigneau, F., Ringwelski, S., 2017, Noise control of vehicle systems, Facta Universitatis- Series Mechanical Engineering, 15(2), pp. 183-200. 4. Goyal, D., Pabla, B.S., 2016, The vibration monitoring methods and signal processing techniques for structural health monitoring: a review, Archives of Computational Methods in Engineering, 23(4), pp. 585-594. 5. Susheel, C.K., Kumar, R., Chauhan, V.S., 2017, Active shape and vibration control of functionally graded thin plate using functionally graded piezoelectric material, Journal of Intelligent Material Systems and Structures, 28(13), pp. 1782-1802. 6. Benjeddou, A., 2000, Advances in piezoelectric finite element modeling of adaptive structural elements: a survey, Computers and Structures, 76(1-3), pp.347-363. 7. Gabbert, U., Köppe, H., Seeger, F., Berger, H. 2002, Modeling of smart composite shell structures, Journal of Theoretical and Applied Mechanics, 3(40), pp. 575-593. 8. Jrad, H., Mallek, H., Wali, M., Dammak, F., 2018, Finite element formulation for active functionally graded thin-walled structures, Comptes Rendus Mécanique, 346(12), pp. 1159-1178. 9. Marinkovic, D, Köppe, H., Gabbert, U., 2006, Numerically efficient finite element formulation for modeling active composite laminates, Mechanics of Advanced Materials and Structures, 13(5), pp. 379–392. 10. Rama, G., Marinkovic, D., Zehn, M.W., 2017, Linear shell elements for active piezoelectric laminates, Smart Structures and Systems, 20(6), pp. 729-737. 11. Cinefra, M., Valvano, S., Carrera, E., 2015, A layer-wise MITC9 finite element for the free-vibration analysis of plates with piezo-patches, International Journal of Smart and Nano Materials, 6(2), pp. 84–104. 12. Carrera, E., Zappino, E., Li, G., 2018, Analysis of beams with piezo-patches by node-dependent kinematic finite element method models, Journal of Intelligent Material Systems and Structures, 29(7), pp. 1379-1393. 13. Carrera, E., Valvano, S., Kulikov, G.M., 2018, Multilayered plate elements with node-dependent kinematics for electro-mechanical problems, International Journal of Smart and Nano Materials, 9(4), pp. 279-317. 14. Marinkovic, D, Köppe, H, Gabbert, U., 2008, Degenerated shell element for geometrically nonlinear analysis of thin-walled piezoelectric active structures, Smart Materials and Structures, 17(1), pp. 1-10. 15. Rama, G., Marinkovic, D., Zehn, M.W., 2018, Efficient three-node finite shell element for linear and geometrically nonlinear analyses of piezoelectric laminated structures, Journal of Intelligent Material Systems and Structures, 29(3), pp. 345-357. 16. Mallek, H., Jrad, H., Wali, M., Dammak, F., 2019, Geometrically nonlinear finite element simulation of smart laminated shells using a modified first-order shear deformation theory, Journal of Intelligent Material Systems and Structures, 30(4), pp. 517-535. 17. Rao, M.N., Schmidt, R., Schröder, K.U., 2018, Static and dynamic FE analysis of piezolaminated composite shells considering electric field nonlinearity under thermo-electro-mechanical loads, Acta Mechanica, 229(12), pp. 5093-5120. 18. Nestorovic, T., Marinkovic, D., Chandrashekar, G., Marinkovic, Z., Trajkov, M., 2012, Implementation of a user defined piezoelectric shell element for analysis of active structures, Finite Elements in Analysis and Design, 52, pp. 11-22. 19. Nestorovic, T, Shabadi, S, Marinkovic, D., Trajkov, M., 2014, User defined finite element for modeling and analysis of active piezoelectric shell structures, Meccanica, 49(8), pp. 1763-1774. 20. Rama, G., Marinkovic, D., Zehn, M.W., 2018, A three-node shell element based on the discrete shear gap and assumed natural deviatoric strain approaches, Journal of the Brazilian Society of Mechanical Sciences and Engineering, 40(7), 356. Abaqus Implementation of a Corotational Piezoelectric 3-Node Shell Element with Drilling Degree... 283 21. Rama, G., Marinkovic, D., Zehn, M.W., 2018, High performance 3-node shell element for linear and geometrically nonlinear analysis of composite laminates, Composites Part B: Engineering, 151, pp. 118-126. 22. Marinkovic, D., Rama, G., 2017, Co-rotational shell element for numerical analysis of laminated piezoelectric composite structure, Composites part B: Engineering, 125, pp. 144-156. 23. Bletzinger, K.U., Bischoff, M., Ramm, E., 2000, A unified approach for shear-locking-free triangular and rectangular shell finite elements, Computers & Structures, 75(3), pp. 321–334. 24. Nguyen-Thoi, T., Phung-Van, P.,Nguyen-Xuan, H., Thai-Hoang, C., 2013, A cell-based smoothed discrete shear gap method (CSDSG3) using triangular elements for static and free vibration analyses of shell structures, International Journal of Mechanical Sciences, 74, pp. 32–45. 25. Felippa, C.A., Militello, C., 1992, Membrane triangles with corner drilling freedoms ii. The ANDES element, Finite Elements in Analysis and Design, 12, pp. 189–201. 26. Bergan, P., Nygard, M., 1984, Finite elements with increased freedom in choosing shape functions, International Journal for Numerical Methods in Engineering, 20(4), pp. 643–663. 27. Park, K., Stanley, G., 1988, Strain interpolations for a 4-node ANS shell element, In: Atluri, S.N., Yagawa, G. (Eds.), Computational Mechanics ‘88, pp. 747–750, Springer, Berlin, Heidelberg. 28. Marinkovic, D., Köppe, H., Gabbert, U., 2007, Accurate modeling of the electric field within piezoelectric layers for active composite structures, Journal of Intelligent Material Systems and Structures, 18(5), pp. 503–513. 29. Felippa, C. A., Haugen, B., 2005, Unifield formulation of small-strain corotational finite elements: I. Theory, Center for aerospace structures, College of engineering, University of Colorado. 30. Bathe, K.J., 1996, Finite Element Procedures, Prentice Hall, New York. 31. Marinkovic, D., Zehn, M.W., Rama, G., 2018, Towards real-time simulation of deformable structures by means of co-rotational finite element formulation, Meccanica, 53(11-12), pp. 3123-3136. 32. Kim, M., Im, S., 2017, A plate model for multilayer graphene sheets and its finite element implementation via corotational formulation, Computer Methods in Applied Mechanics and Engineering, 325, pp. 102-138. 33. Tzou, H.S., 1993, Piezoelectric Shells: Distributed Sensing and Control of Continua, Kluwer Academic Publishers, Amsterdam. 34. Zhang, S., 2014, Nonlinear FE Simulation and Active Vibration Control of Piezoelectric Laminated Thin- Walled Smart Structures, PhD thesis, Institute of General Mechanics, RWTH Aachen University. 35. Saravanos, D.A., Heyliger, P.R, Hopkins, D.A., 1997, Layerwise mechanics and finite element for the dynamic analysis of piezoelectric composite plates, International Journal of Solids and Structures, 34(3), pp. 359–378.