Multilayered plate elements with node-dependent kinematics for electro-mechanical problems

ABSTRACT In the present work, a new class of finite elements (FEs) for the analysis of composite and sandwich plates embedding piezoelectric skins and patches is proposed. By making use of node-by-node variable plate theory assumptions, the new finite element allows for the simultaneous analysis of different subregions of the problem domain with different kinematics and accuracy, in a global/local sense. As a consequence, the computational costs can be reduced drastically by assuming refined theories only in those zones/nodes of the structural domain where the resulting strain and stress states, and their electro-mechanical coupling present a complex distribution. The primary advantage is that no ad-hoc techniques and mathematical artifices are required to mix the fields coming from two different and kinematically incompatible adjacent elements, because the plate structural theory varies within the finite element itself. In other words, the structural theory of the plate element is a property of the FE node in this present approach, and the continuity between two adjacent elements is ensured by adopting the same kinematics at the interface nodes. The finite element arrays of the generic plate element are formulated in terms of fundamental nuclei, which are invariants of the theory approximation order and the modeling technique (Equivalent-Single-Layer, Layer-Wise). In this work, the attention is focused on the use of Legendre polynomial expansions to describe the through-the-thickness unknowns to develop advanced plate theories. Several numerical investigations, such as composite and sandwich multilayered plates embedding piezoelectric skins and patches with various load, boundary conditions, and piezoelectric material polarizations, are carried out to validate and demonstrate the accuracy and efficiency of the present plate element, including comparison with various closed-form and FE solutions from the literature.


Introduction
Plate structures have a predominant role in a variety of engineering applications, such as piezoelectric composite structures modeling, and in the overall design procedure for smart structures and systems. The use of piezoelectric components as electro-mechanical transducers in sensor as well as in actuator applications embedded in layered CONTACT  composite structures is complicated in practice. In some cases, structures may contain regions where three-dimensional (3D) stress fields occur. To accurately capture these localized 3D stress states, solid models or higher-order theories are necessary. Analytical solution for general smart structural problems is a very tough task, and they exist, only, for few specialized and idealized cases. Meanwhile, the finite element method (FEM) has become the most widely used technique to model various physical processes, including piezoelectricity. However, the high computational costs represent the drawback of refined plate theories or three-dimensional analyses. The implementation of innovative solutions for improving the analysis efficiency for complex geometries and assemblies, possibly in a global/local scenario is the main motivations for this work.
The fundamentals of the modeling of piezoelectric materials have been given in many contributions, in particular in the pioneering works of Mindlin [1], EerNisse [2], Tiersten and Mindlin [3], and in the monograph of Tiersten [4]. The embedding of piezoelectric layers into plates and shells sharpens the requirements of an accurate modeling of the resulting adaptive structure due to the localized electro-mechanical coupling, see e.g. the review of Saravanos and Heyliger [5]. Therefore, within the framework of two-dimensional approaches, layerwise descriptions have often been proposed either for the electric field only (see e.g. the works of Kapuria [6] and Ossadzow-David and Touratier [7]) or for both the mechanical and electrical unknowns (e.g. Heyliger et al. [8]). Ballhause et al. [9] showed that a fourth order assumption for the displacements leads to the correct closed form solution. They conclude that the analysis of local responses requires at least a layer-wise descriptions of the displacements, see also [10]. Benjeddou et al. [11] emphasized that a quadratic electric potential through the plate thickness satisfies the electric charge conservation law exactly. A layerwise mixed finite element is used for piezolaminated plates in [12], and a layerwise mixed least-squares model model is used in [13]. Some of the latest contributions to the Finite Elements (FEs) analysis of piezoelectric plates and shells that are based on exact geometry solid-shell element was developed by Kulikov et al. [14,15], composite laminates consisting of passive and multi-functional materials were analyzed in [16], therefore some important aspects of modeling piezoelectric active thin-walled structures were treated in [17] Although the enormous improvements and formulations of higher-order plate structural theories, considerable work has been recently directed towards the implementation of innovative solutions for improving the analysis efficiency for complex geometries and assemblies, possibly in a global/local scenario. In this manner, the limited computational resources can be distributed in an optimal manner to study in detail only those parts of the structure that require an accurate analysis. In general, two main approaches are available to deal with a global/local analysis: (1) refining the mesh or the FE shape functions in correspondence with the critical domain; (2) formulating multi-model methods, in which different subregions of the structure are analyzed with different mathematical models. The coupling of coarse and refined mesh discretizations, or different FE shape functions, can be addressed as single-theory or single-model methods. The h-adaption method [18] is used when the structures subregions differ in mesh size, whereas the p-adaption method [19] can be applied when the subregions vary in the polynomial order of the shape functions. Moreover, the hp-adaption [20] can allow the implementation of subregions differing in both mesh size and shape functions.
In the case of multi-theory methods, in which different subregions of the structure are analyzed with different structural theories with kinematically incompatible elements, the compatibility of displacements and equilibrium of stresses at the interface between dissimilar elements have to be achieved. A wide variety of multiple model methods has been reported in the literature. In general, multi-theory methods can be divided into sequential or multistep methods and simultaneous methods. In a sequential multimodel, the global region is analyzed with an adequate model with a cheap computational cost to determine the displacement or force boundary conditions for a subsequent analysis at the local level. The local region can be modeled with a more refined theory, or it can be modeled with 3-D finite elements, see [21][22][23][24]. The simultaneous multi-model methods are characterized by the analysis of the entire structural domain, where different subregions are modeled with different mathematical models and/or distinctly different levels of domain discretization, in a unique step. One of the simplest type of simultaneous multi-model methods for composite laminates analysis is the concept of selective ply grouping or sublaminates. Recently, the authors developed multi-model elements with variable through-the-thickness approximation by using 2-D finite elements for both local and global regions [25][26][27][28]. In this approach, the continuity of the primary variables between local and global regions was straightforwardly satisfied by employing Legendre polynomials. Another well-known method to couple incompatible kinematics in multi-model methods is the use of Lagrange multipliers, which serve as additional equations to enforce compatibility between adjacent subregions. In the three-field formulation by Brezzi and Marini [29], an additional grid at the interface is introduced. The unknowns are represented independently in each sub-domain and at the interface, where the matching is provided by suitable Lagrange multipliers. This method was recently adopted by Carrera et al. [30] to couple beam elements of different orders and, thus, to develop variable kinematic beam theories. Ben Dhia et al. [31][32][33] proposed the Arlequin method to couple different numerical models using a minimization procedure. This method was adopted by Hu et al. [34,35] for the linear and nonlinear analysis of sandwich beams modelled via one-dimensional and two-dimensional finite elements, and by Biscani et al. [36] for the mechanical analysis of beams, by Biscani et al. [37] for the mechanical analysis of plates, and by Biscani et al. [38] for the electromechanical analysis of plates. In the present work, a new simultaneous multiple-model method for 2D elements with node-dependent kinematics is developed, for the analysis of electro-mechanical problems. This node-variable capability enables one to vary the kinematic assumptions within the same finite plate element. The expansion order of the plate element is, in fact, a property of the FE node in the present approach. Therefore, between finite elements, the continuity is ensured by adopting the same expansion order in the nodes at the element interface. This node-dependent finite element has been used by the authors for the mechanical analysis in [39,40] using classical and HOTtype theories with Taylor polynomials were used with an ESL approach, and the combination of HOT-type and advanced LW theories in the same finite element. In this manner, global/local models can be formulated without the use of any mathematical artifice. As a consequence, computational costs can be reduced assuming refined models only in those zones with a quasi-three-dimensional stress field, whereas computationally cheap, low-order kinematic assumptions are used in the remaining parts of the plate structure. In this paper, the governing equations of the variable-kinematics plate element for the linear static coupled electro-mechanical analysis of composite structures are derived from the Principle of Virtual Displacement (PVD). Subsequently, FEM is adopted and the Mixed Interpolation of Tensorial Components (MITC) method [41][42][43][44] is used to contrast the shear locking. The developed methodology is, therefore, assessed and used for the analysis of composite and sandwich multilayered plates embedding piezoelectric skins and patches with various load, boundary conditions, and piezoelectric material polarizations. The results are compared with various closed-form and Fem solutions and, whenever possible, with exact solutions available from the literature.

Preliminaries for electro-mechanical problems of plates
Plates are bi-dimensional structures in which one dimension (in general the thickness in the z direction) is negligible with respect to the other two dimensions. The geometry and the reference system are indicated in Figure 1.
The constitutive equations for coupled electro-mechanical problems permit to relate the mechanical stresses σ ¼ σ xx ; σ yy ; σ xy ; σ xz ; σ yz ; σ zz Â Ã , and the electric displacements to the mechanical strains ¼ xx ; yy ; xy ; xz ; yz ; zz Â Ã , and the electric , for each layer k, in the following compact form: where C is the matrix of the material stiffness coefficients, e is the matrix of the piezoelectric stiffness coefficients, and ε is the matrix of the permittivity coefficients. The mechanical strains and the electric field E are related to the mechanical displacements u ¼ u; v; w ½ and the electric potential Φ via the geometrical relations as follows: Figure 1. Geometry and reference system of the multilayered plate, including piezoelectric skins and patches.
where D g and D eg are the vectors containing the differential operators defined as follows: (3) The matrix of the material stiffness coefficients and the permittivity coefficients for orthotropic materials are defined as follows: For the sake of brevity, the expressions that relate the material coefficients C ij to the Young's moduli E 1 , E 2 , E 3 , the shear moduli G 12 , G 13 , G 23 and Poisson ratios ν 12 , ν 13 , ν 23 , ν 21 , ν 31 , ν 32 are not given here, they can be found in [45]. The piezoelectric material is characterized by the piezoelectric coefficients e ij and the permittivity coefficients ε ij , more details can be found in the book of Rogacheva [46].

Actuation modes
The piezoelectric effect is made available in polarisable crystalline materials through the application of an intense electric field which imparts a net polarization of the crystal cells. Depending on the mutual direction of the polarization and the applied loading, two important actuation modes are taken into account in this work.
2.1.1. Actuation in 3-1 mode Transverse extension mode (3-1 mode), the applied electric field is aligned with the polarization axis, but the major deformation occurs in the transverse plane due to the thinness of the piezoelectric sheet, see Figure 2. The piezoelectric coefficients for the 3-1 mode are defined as follows:

Actuation in 1-5 mode
Shear mode (1-5 mode), the applied electric field is perpendicular to the polarization direction, and the principal mechanical effect is associated with a shear deformation, see Figure 2. The piezoelectric coefficients for the 1-5 mode are defined as follows:

Refined and hierarchical plate theories for electro-mechanical problems
This work proposes a class of new finite elements which allows employing different kinematic assumptions in different subregions of the problem domain. To highlight the capabilities of the novel formulation, a four-node plate elements with nodedependent kinematics is shown in Figure 3. The element proposed in this example makes use of fourth order Layer-Wise models for the local zones that includes piezoelectric patches, a linear Layer-Wise model for the short free-edge tip, and Equivalent-Single-Layer models for the global multilayered zones without the piezoelectric patches. As it will be clear later in this paper, thanks to the hierarchical capabilities of Unified Formulation, the choice of the nodal plate theory is arbitrary, and variable-kinematics plate elements will be used to implement multi-model methods for global-local analysis.
Classical and Higher Order Theories. Classical plate models grant good results when small thickness, homogeneous structures are considered. On the other hand, the analysis of thick plates and multilayered structures may require more sophisticated theories to achieve sufficiently accurate results. As a general guideline, it is clear that the richer the kinematics of the theory, the more accurate the 2D model becomes. To overcome the limitations of classical theories, a large variety of plate higher-order theories (HOT) have been proposed in the past and recent literature. Eventually, higher-order theories can be expressed by making use of Taylor-like expansions of the generalized unknowns along the thickness. The fundamentals of the modeling of piezoelectric materials, with classical plate models, have been given in many contributions, in particular in the pioneering works of Mindlin [1], EerNisse [2], Tiersten and Mindlin [3], and in the monograph of Tiersten [4]. The localized electro-mechanical coupling, due to the use of piezoelectric layers in multilayered structures, often lead to the formulation of layerwise descriptions either for the electric field only (see e.g. the works of Kapuria [6] and of Ossadzow-David and Touratier [7]) or for both the mechanical and electrical unknowns (e.g. Heyliger et al. [8]). In this work, the attention is focused on the use of layer-wise description of both mechanical and electrical variables, with a general expansion of N terms.

Unified formulation for plates
The Unified Formulation has the capability to expand each displacement variable at any desired order. Each variable can be treated independently from the others, according to the required accuracy. This procedure becomes extremely useful when multifield problems are investigated such as thermoelastic and piezoelectric applications [27,28,47,48]. According to the UF by Carrera [49][50][51], the displacement field and the electric potential can be written as follows: In compact form: where ðx; y; zÞ is the general reference system (see Figure 1), the displacement vector u ¼ u; v; w f g and the electric potential Φ have their components expressed in this system. δ is the virtual variation associated to the virtual work, and k identifies the layer. F τ and F s are the thickness functions depending only on z: τ and s are sum indexes and N is the number of terms of the expansion in the thickness direction assumed for the displacements. For the sake of clarity, the superscript k is omitted in the definition of the Legendre polynomials.

Legendre-like polynomial expansions
The limitations, due to expressing the unknown variables in function of the midplane position of the shell, can be overcome in several ways. A possible solution can be found employing the Legendre polynomials. They permit to express the unknown variables in function of the top and bottom position of a part of the shell thickness. In the case of Legendre-like polynomial expansion models, the displacements and the electric potential are defined as follows: s ¼ 0; 1; r ; r ¼ 2; :::; N: in which P j ¼ P j ðζÞ is the Legendre polynomial of j-order defined in the ζ-domain: For the Layer-Wise (LW) models, the Legendre polynomials and the relative top and bottom position are defined for each layer.

Finite elements with node-dependent kinematics
Thanks to Unified Formulation, FEM arrays of classical to higher-order plate theories can be formulated in a straightforward and unified manner by employing a recursive index notation. By utilizing an FEM approximation, the generalized displacements of Equation (12) can be expressed as a linear combination of the shape functions to have u s ðx; yÞ ¼ N j ðx; yÞu s j j ¼ 1; :::; ðnodes per elementÞ Φ s ðx; yÞ ¼ N j ðx; yÞΦ s j j ¼ 1; :::; ðnodes per elementÞ (16) where u s j and Φ s j are the vectors of the mechanical and electrical, respectively, generalized nodal unknowns and N j can be the usual Lagrange shape functions. j denotes a summation on the element nodes. Since the principle of virtual displacements in used in this paper to obtain the elemental FE matrices, it is useful to introduce the finite element approximation of the virtual variation of the generalized displacement vector δu τ , In Equation (17), δ denotes the virtual variation, whereas indexes τ and i are used instead of s and j, respectively, for the sake of convenience. In this work, and according to Equations (12), (16) and (17), the thickness functions F s and F τ , which determine the plate theory order, are independent variables and may change for each node within the plate element. Namely, the three-dimensional displacement field and the related virtual variation can be expressed to address FE nodedependent plate kinematics as follows: uðx; y; zÞ ¼ F j s ðzÞN j ðx; yÞu s j s ¼ 0; 1; :::; N j j ¼ 1; :::; ðnodes per elementÞ Φðx; y; zÞ ¼ F j s ðzÞN j ðx; yÞΦ s j s ¼ 0; 1; :::; N j j ¼ 1; :::; ðnodes per elementÞ δuðx; y; zÞ ¼ F i τ ðzÞN i ðx; yÞδu τ i τ ¼ 0; 1; :::; N i i ¼ 1; :::; ðnodes per elementÞ δΦðx; y; zÞ ¼ F i τ ðzÞN i ðx; yÞδΦ τ i τ ¼ 0; 1; :::; N i i ¼ 1; :::; ðnodes per elementÞ (18) where the subscripts τ, s, i, and j denote summation. Superscripts i and j denote node dependency, such that for example F i τ is the thickness expanding function and N i is the number of expansion terms at node i, respectively. For the sake of clarity, the displacement and electric fields of a variable kinematic plate element, for example the fourn node element represented in Figure 4, is described in detail hereafter. The global displacement field of the element is approximated as follows: • Node 1 Plate Theory = HOT with N 1 = 2 • Node 2 Plate Theory = LW with N 2 = 2 • Node 3 Plate Theory = HOT with N 3 = 4 • Node 4 Plate Theory = HOT with N 4 = 3 According to Equation (18), it is easy to verify that the displacements at a generic point belonging to the plate element can be expressed as given in Equation (19). In this equation, only the displacement component along z-axis is given for simplicity reasons: It is intended that, due to node-variable expansion theory order, the assembling procedure of each finite element increases in complexity with respect to classical mono-theory finite elements. In the present FE approach, in fact, it is clear that both rectangular and square arrays are handled and opportunely assembled for obtaining the final elemental matrices.

Fundamental nucleus of the stiffness matrix
Given Unified Formulation and FE approximation, the governing equations for the static response analysis of the multi-layer plate structure can be obtained by using the principle of virtual displacements, which states: where the term on the left-hand side represents the virtual variation of the strain energy; Ω and A are the integration domains in the plane and the thickness direction, respectively; ε and σ are the vector of the strain and stress components; and δL e is the virtual variation of the external loadings. By substituting the constitutive equations for composite elastic materials Equation 1, the linear geometrical relations Equation 2 as well as Equation (18) into Equation (20), the linear algebraic system in the form of governing equations is obtained in the following matrix expression: where K τsij and P τi are the element stiffness and load FE arrays written in the form of fundamental nuclei. The mechanical part K kτsij uu is a 3 Â 3 matrix, the coupling matrices K kτsij uΦ , K kτsij Φu have dimension 3 Â 1 and 1 Â 3 respectively, and the electrical part K kτsij ΦΦ is a 1 Â 1 matrix. The stiffness matrix nucleus K ktsij uu , and the mechanical external load vector nucleus P kti u are the same defined for the pure mechanical problems, reader can refers to the work [39]. The explicit expression of the stiffness electro-mechanical coupling matrices and the pure electric nucleus are given below with classical FEM method. The pure electric contribution is defined as follows: The stiffness electro-mechanical coupling matrices K ktsij uΦ and K ktsij Φu are defined as follows: where comma denote partial derivatives respect to the spatial directions. The fundamental nucleus as given above is the basic building block for the construction of the element stiffness matrix of classical, refined and variable-kinematic theories. In fact, given these nine components, element stiffness matrices of arbitrary plate models can be obtained in an automatic manner by expanding the fundamental nucleus versus the indexes τ, s, i, and j: In the development of ESL theories as in the case of this paper, the fundamental nucleus of the stiffness matrix is evaluated at the layer level and then assembled as shown in Figure 5. This figure, in particular, illustrates the expansion of the fundamental nucleus in the case of a 9-node Lagrange finite element with nodeby-node variable kinematics, as in the case of this paper. It must be added that, in this work, an MITC technique is used to overcome the shear locking phenomenon, see [48]. However, for more details about the explicit formulation of the Unified Formulation fundamental nuclei, interested readers are referred to the recent book by Carrera et al. [52].

Numerical results
Some problems have been considered to assess the capabilities of the proposed variable-kinematics plate elements and related global/local analysis. These analysis cases comprise both composite and sandwich laminated plate structures with different boundary conditions and loadings. Whenever possible, the proposed multi-theory models are compared to single-theory refined elements. According to Unified Formulation terminology, the latter models are referred to as LWN, where LW stands for Layer-Wise models (LW), and N is the theory approximation order. Eventually, 3D exact or finite element models, and Arlequin multi-model solutions are used for comparisons and, in those cases, the opportune notation is mentioned case by case. If a Navier-type closed form solution is employed instead of FEM, the subscript ðaÞ is used. On the contrary, for the sake of clarity, multi-model theories are opportunely described for each numerical case considered.

Simply-supported cross-ply composite plates with piezoelectric skins
To assess this new finite element, a four-layer cross-ply square plate with a cross-ply Gr/ Ep composite core ½0 /90 and PZT-4 piezoelectric external skins is analyzed, see Figure 6. The square plate has the following geometrical data: a ¼ b ¼ 4:0, and h tot ¼ 1:0. In respect to the total thickness, a single piezoelectric skin is thick h p ¼ 0:1h tot , while the single core layer is thick h c ¼ 0:4h tot . The static analysis of the plate structure is evaluated in sensor and actuator configuration.
For the sensor case, a bi-sinusoidal transverse normal pressure is applied to the top surface of the plate: with amplitude p o z ¼ 1 and wave numbers m ¼ 1, n ¼ 1. The potential at top and bottom position is imposed Φ t ¼ Φ b ¼ 0. For the actuator case, a bi-sinusoidal electric potential is imposed at top surface: with amplitude ϕ o z ¼ 1 and wave numbers m ¼ 1, n ¼ 1. The potential at bottom position is imposed Φ b ¼ 0. No mechanical load is applied. The material properties of the plate are given in Table 1.
The plate has simply-supported boundary conditions. Due to the symmetry of both the geometry and the load, a quarter of the plate is analyzed and the following symmetry and boundary conditions (simply-supported) are applied: The boundary condition of the electric potential is taken into account to compare the results with the analytical solution, see [9], where the electric potential has the following Navier-type assumptions Φ ¼Φsinðmπx=aÞsinðnπy=bÞ.
To compare the results with other solutions present in literature [38](Arlequin method), the mid-plane domain of the plane structure was subdivided into two zones along the axes x and y, as shown in Figure 7, and multi-theory models Case A, CaseB and Case C are depicted on the FE discretization of a quarter of the plate. The FE mesh on a quarter plate is 10 Â 10 elements, the accounted mesh is the same used in the reference    [38]; the difference is that, in the reference solutions, it is used a four-node element and a Reissner Mixed Variational Theorem applied to the transverse electric displacement D z (RMVT À D z ), the transverse electric displacement is a priori modelled with the mechanical displacements. It has to be noticed that the mesh discretization of the present multi-models CaseA and CaseC is the same of the reference Arlequin solutions named ðLW1 À LWM3Þ A and ðLW2 À LWM3Þ C .
Some results of the transverse mechanical displacement w, in-plane stress σ xx , transverse shear stress σ xz , transverse normal stress σ zz , electric potential Φ, and transvserse electric displacement D z evaluated along the plate thickness are given in tabular form, see Tables 2 and 3, for the sensor case and the actuator case respectively. Monotheory models are compared with those from the present multi-model approach, furthermore, the exact 3D solutions provided by Heyliger [53], the analytical solution with layer-wise mono-models [9], and the multi-model solutions via the Arlequin method [38] are given. For the transverse shear and normal stresses, the evaluation point is located between the composite layers, their results are given for the upper layer with upscript þ and for the lower layer with upscript À . It is clear for both sensor and actuator cases that at least a second order expansion order is needed to get good results for all the considered variables.
Some results are given in terms of transverse displacement w, transverse shear stress σ xz and electric transverse displacement D z along the plate thickness. For the transverse displacement w, in sensor case configuration, the differences between single and multimodels are negligible, see Figure 8(a). On the contrary for the actuator case, see Figure 8 (b), remarkable differences are present between LW4 and LW1 single-model solutions; extending the comparison to the multi-models, it is clear that a more refined expansion is needed in the center plate where the load is higher, see CaseA and CaseC, differently Case B shows an accuracy very close to LW1 single-model.
The transverse shear stress σ xz is represented in Figure 9(a) for the actuator configuration. The stress is evaluated in its maximum shear stress value point ðx; yÞ ¼ ð0; b=2Þ. The mono-model LW4 is able to predict the correct behavior as the 3DExact reference solution [53]. The higher-order multi-models CaseA and CaseC are not sufficient to depict the good representation of the shear stress. Differently CaseB shows a good accuracy solution due to the higher-order representation in the evaluation zone of the shear stress.
Regarding the electric transverse displacement D z , represented in Figure 9(b) for the sensor configuration, it is evaluated in the center plate. The mono-model LW4 can predict the correct behavior as the 3DExact reference solution [53]. In this case, the    Results in terms of in-plane stress σ xx ðy; zÞ ¼ σ xx , and transverse electric displacement D z ¼ ð10 10 Þ Â D z along the in-plane x axis are represented in Figure 10(a,b) respectively. For both the depicted variables, the multi-model CaseC is able to reproduce the correct     . Composite plate with piezoelectric skins. In-plane stress σ xx ðy; zÞ ¼ σ xx ðb=2; þh=2Þ for the Actuator Case, and transverse electric displacement D z ðy; zÞ ¼ ð10 10 Þ Â D z ðb=2; þh=2Þ for the Sensor Case, along the in-plane direction X, the axis X is expressed in ½mm. Single and Multi-theory models.
behavior along the X axis as the reference solution LW4. The linear single model LW1 and multi-models CaseA and Case B show an incorrect solution in the zone where a linear approximation is used. In the transition elements small oscillations are present, these oscillations are due to the coarse mesh, if the mesh is more refined the oscillations tend to fade. Figure 11(a,b) show the three-dimensional distributions of the transverse normal stress σ zz , in Sensor Case configuration, of the single-model LW4 and the variable kinematic multi-model CaseA respectively. The results show the enhanced global/local capabilities of the CaseA model, which is able to predict correctly the stress state in the center zone where the loading is bigger. Moreover the three-dimensional distributions of the transverse shear stress σ xz , in Actuator Case configuration, is depicted in Figure 12(a,b) for the singlemodel LW4 and the variable kinematic multi-model CaseB respectively. The results show the enhanced global/local capabilities of the CaseB model, it represents accurately the maximum shear stress in boundary zone of the plate structure.

Sandwich cantilever plate shear-actuated in mode-15
A cantilever sandiwch plate is analysed as second example and it is shown in Figure 13. F=m; a small PZT-5H piezoelectric patches is introduced in the foam layer with dimension: a p ¼ 10 mm, b p ¼ 20 mm, h p ¼ 2 mm, the PZT-5H material has the following properties: The material PZT-5H is polarized in the x-direction, or in mode À 15. The structure is loaded at the upper and lower surfaces of the piezoelectric patch with a constant uniform electric potential equal to Φ t=b ¼ Ç10:0V: First, a convergence study on a single-theory plate model was performed. As far as an LW4 model is concerned and as shown in Table 4, a non-uniform mesh grid of 56 Â 8 elements, see Figure 14 is enough to ensure convergent results. This structural problem has become popular and it is a good benchmark problem for its selectivity. The present plate element model is compared with various solutions from the literature, including those of Zhang and Sun [54,55], Benjeddou et al. [56]. Some results are given varying the actuator position along the x-axis, the deflection at the free edge is investigated for each position of the piezoelectric patch, see Figure 15, and compared with the literature solutions [54,56], here named Sun & Zhang, and Benjeddou et al., respectively.
Various node-variable kinematic models have been used to perform the global/local analysis of the proposed plate structure, with the center of the piezoelectric patch fixed at d ¼ 85 mm. The mid-plane domain of the plate structure was subdivided into higherorder and lower-order zones along the axes x and y and they are depicted in Figure 16. The mesh discretization of the present multi-models is arbitrary. The idea behind the Table 4. Convergence study versus the number of elements of the LW4 single-theory model of the sandwich cantilever plate. Transverse displacement w ¼ 10 8 Â wða; b=2; þh=2Þ, electric potential Φ ¼ Φðd; b=2; Àh=2Þ, in-plane principal stress σ xx ¼ σ xx ðd; b=2; Àh=2Þ, transverse shear stress σ xz ¼ σ xz ðd; b=2; 0Þ, with the center of the piezoelectric patch placed at d ¼ 85 mm.  Figure 14. Mesh zone subdivisions of the sandwich plate with piezoelectric patch, for the convergence study.
discretization consists into the study of the effect of the transition zones position respect with to the variables evaluation points. Some results of the transverse mechanical displacement w, in-plane stress σ xx , transverse shear stress σ xz , transverse normal stress σ zz , electric potential Φ, and transverse electric displacement D z evaluated along the plate thickness are given in tabular form, see Table 5. Mono-theory models are compared with those from the present multi-model approach, furthermore the FEM 3D solution provided by 3D Abaqus C3D20RE element is given. Some results in terms of transverse displacement w, and electric potential Φ along the thickness are represented in Figure 17(a,b), furthermore three-dimensional view of the electric potential Φ is given in Figure 18(a,b). Some more comments can be made: • The through-the-thickness distribution of the transverse displacement w at the free tip, as shown in Figure 17(a), is correctly predicted by higher-order single-models LW3 and LW4. The same accuracy is reached by the Case A multi-model and little losses in accuracy are present in the Case B and Case C multi-models. • The behavior of the electric potential Φ along the thickness, depicted in Figure 16, is well described for every single and multi-models. Furthermore the three-dimensional view of the electric potential Φ, on deformed structure, is given by the finite element 3D Abaqus C3D20RE, see Figure 18(a), and the present mono-model LW4, see Figure 18(b). It has noticeable that the electric potential calculated by the commercial 3D Abaqus C3D20RE finite element does not tend to zero at the top and bottom positions of the inserted patch zone, differently the present LW4 model well describe the electric potential behavior without imposing any boundary Figure 15. Sandwich plate with piezoelectric patch, tip transverse displacement wðx; y; zÞ ¼ 10 8 Â wða; b=2; þh=2Þ for several position of the piezoelectric patch along the x-axis direction. Figure 16. Graphical representation of the multi-theory models, based on layer-wise models, of the sandwich plate with piezoelectric patch, for the node-variable kinematic study. Table 5. Single-theory and multi-theory models of the sandwich cantilever plate. Transverse displacement w ¼ 10 8 Â wða; b=2; þh=2Þ, electric potential Φ ¼ Φðd; b=2; Àh=2Þ, in-plane principal stress σ xx ¼ σ xx ðd; b=2; Àh=2Þ, transverse shear stress σ xz ¼ σ xz ðd; b=2; 0Þ, transverse normal stress σ zz ¼ σ zz ðd; b=2; Àh=2Þ, transverse electric displacement D z ¼ 10 4 Â D z ðd; b=2; 0Þ, with the center of the piezoelectric patch placed at d ¼ 85 mm. conditions at the top and bottom positions. Furthermore, some results in terms of mechanical stresses are given for the in-plane stress σ xx in Figure 19(a), transverse normal stress σ zz in Figure 19(b), and transverse shear stress σ xz in Figure 20(a,b) and 21(a,b). Some more comments can be made: • The through-the-thickness distribution of the in-plane stress σ xx , as shown in Figure 19(a), is correctly predicted by higher-order single-models LW3 and LW4.
The same accuracy is reached by all the considered multi-models, where the evaluation point is described by the LW4 theory. It has to be noticed that the CaseB show a little loss in accuracy due to the short distance of the evaluation point from the transition elements.  Figure 18. Sandwich cantilever plate, three-dimensional view of the electric potential Φ, on deformed structure. 3D Abaqus C3D20RE and mono-model LW4.
• The transverse normal stress σ zz , as shown in Figure 19(b), is well described by the LW4 single-model. The same accuracy is reached by all the considered multimodels, where the evaluation point is described by the LW4 theory. As mentioned for the in-plane stress σ xx , the CaseB show a little loss in accuracy due to the short distance of the evaluation point from the transition elements. • The three-dimensional view of the transverse shear stress σ xz , on undeformed structure, is given by the finite element 3D Abaqus C3D20RE in Figure 20(a), by the present mono-model LW4 in Figure 20(b), and by multi-model Case C in Figure 21(a), in which it is possible to notice the differences, at the clamped boundaries and at the transition elements close to the patch zone, respect to  Figure 19. Sandwich cantilever plate, in-plane stress σ xx ðx; yÞ ¼ σ xx ðd; b=2Þ, and transverse normal stress σ zz ðx; yÞ ¼ σ zz ðd; b=2Þ. Single and Multi-theory models. Step: Step−1 Increment 1: Step Time = 1.000 Primary Var: S, S13 the LW4 solution. Therefore, the through-the-thickness distribution of the transverse shear stress σ xz is given in Figure 21(b). All the single-model considered are not able to fulfill the interlaminar continuity condition of the shear stress.
Regarding the multi-models, they show an accuracy very close to the LW4 solution, and as mentioned for the in-plane stress σ xx , the Case B show a little loss in accuracy due to the short distance of the evaluation point from the transition elements.
By the evaluation of the various node-variable kinematic models, it is clear that an accurate representation of the stresses in localized zones is possible with DOFs reduction if an accurate distribution of the higher-order kinematic capabilities is performed. Differently, the primary variables, mechanical displacements and electric potential, are dependent on the global approximation over the whole structure. The DOFs reduction can be moderate or stronger, depending on the structure and the load case configuration.

Sandwich cantilever plate under mechanical loading
A cantilever sandiwch plate is analysed as shown in Figure 22. The geometrical dimensions are: a ¼ b ¼ 20 mm, h total ¼ 6 mm. The upper and lower layers are made of Aluminum. The central layer is made of Foam, and a PZT-5H piezoelectric patches is introduced in the foam layer with dimension: a p ¼ 10 mm, b p ¼ 20 mm, h p ¼ 2 mm, it is centered at x ¼ 10 mm. The material properties are the same of the previous numerical example. The three layers have the same thickness h ¼ 2 mm. The structure is loaded at the free tip ðx; y; zÞ ¼ ða; b=2; þh=2Þ with a concentrated transverse mechanical load The structure analyzed in this numerical section is taken from the work of Sun and Zhang [54]. The present single-and multi-model solutions are compared with a calculated three-dimensional FEM ABAQUS solution. A non-uniform mesh grid of 60 Â 16 elements ensures the convergence of the solution with a LW4 single-model, see Figure 23. For the sake of brevity, the study of the convergence is here omitted. The adopted refined mesh is necessary to study the behavior of the mechanical and electrical variables along the whole plate domain, and not in one single point. The difficult task is to obtain a good behavior of the mechanical stresses, electric potential, and electric displacements, along with the inplane directions close to the interfaces of the piezoelectric patch, avoiding strange oscillations due to the changing of the element size.
Various node-variable kinematic models have been used to perform the global/local analysis of the proposed plate structure. The mid-plane domain of the plate structure was subdivided into different higher-and lower-order zones along the axes x and y and they are depicted in Figure 23. The mesh discretization of the present multi-models is arbitrary. The idea behind the discretization consists into the study of the effect of the transition zones position respect with to the variables evaluation points, and respect with to the material change in the particular case of the foam core with the inserted piezoelectric patch. Some results of the transverse mechanical displacement w, in-plane stress σ xx , transverse shear stress σ xz , transverse normal stress σ zz , electric potential Φ, and transverse electric displacement D z evaluated along the plate thickness are given in tabular form, see Table 6. Mono-theory models are compared with those from the present multi-model approach, furthermore the FEM 3D solution provided by 3D Abaqus C3D20RE element is given.
Some results are given in terms of transverse displacement w, transverse shear stress σ xz , electric potential Φ and electric transverse displacement D z along the plate thickness. For the transverse displacement w the differences between single and multimodels are negligible in the lower part of the multilayer, see Figure 24(a). On the contrary, in the upper part of the multilayer close to the applied concentrated load, remarkable differences are present between LW4 and CaseA respect to the other single and multi-model solutions.
The transverse shear stress σ xz is represented in Figure 24(b). The stress is evaluated in the center patch ðx; yÞ ¼ ða=2; b=2Þ. The mono-model LW4 can predict the correct Figure 23. Non-uniform adopted mesh and graphical representation of the multi-model cases, for the sandwich plate.
behavior satisfying the interlaminar continuity condition, and its accuracy is almost the same of the Abaqus C3D20RE finite element solution. The lower single-models LW2 and LW1 are not able to represent correctly the stress behavior. The higher-order multimodels show a good accuracy solution due to the higher-order representation in the evaluation zone of the shear stress.
Higher-order single-models are needed to well describe to non-linear behavior of the electric potential and to capture its maximum value located at the interfaces corner ð3a=4; b=2; þh=6Þ. The top and bottom position values tend naturally to zero without imposing any boundary conditions. The multi-model solutions have almost the same accuracy of the LW4 solution, except for the Case B multi-model which shows an increase of the maximum value at the interfaces corner, this is due to the influence of the transition zone with the LW1 zone elements, as shown in Figure 23. It has to be noticeable that the present solutions are compared with the Abaqus C3D20RE finite element solution which shows a comparable electric potential description in the center part of the thickness multilayer, on the contrary the top and bottom values do not naturally tend to zero.
Regarding the electric transverse displacement D z , represented in Figure 25(b), it is evaluated in the center patch. The same considerations of the shear stress can be made here. The mono-model LW4 is able to predict the correct behavior, its accuracy is almost the same of the Abaqus C3D20RE finite element solution. The lower single-models LW2 and LW1 are not able to represent correctly the electric displacement behavior. The higher-order multi-models show a good accuracy solution due to the higher-order representation in the evaluation zone.
Results in terms of transverse shear and normal stresses σ xz ðy; zÞ ¼ ð10 À5 Þ Â σ xz , σ zz ðy; zÞ ¼ ð10 À5 Þ Â σ zz , electric potential Φðy; zÞ ¼ Φ, and transverse electric displacement D z ¼ ð10 5 Þ Â D z , along the in-plane x axis at the interface between the upper skin and the sandwich core, are represented in Figure 26(a,b) and 27(a,b) respectively. For both the transverse stress variables, see Figure 26(a,b), the LW4 single-model and higher-order multi-models show the same behavior and accuracy. Higher peak values are noticeable at the side-edges of the piezoelectric patch x ¼ 5; 15 mm. The multi- model Case B show an increase of the maximum peak value at x ¼ 15 mm, this is due to the transition zone between LW4 and LW1 models, as shown in Figure 23. The linear single-model LW1 completely underestimate the stresses description.  Figure 26. Sandwich cantilever plate under concentrated mechanical load. Transverse shear stress σ xz ðy; zÞ ¼ 10 À5 Â σ xz ðb=2; þh=6Þ, and transverse normal stress σ zz ðy; zÞ ¼ 10 À5 Â σ zz ðb=2; þh=6Þ along the in-plane x-axis. Single and Multi-theory models. The electric potential is well depicted by all the single and multi-models, as shown in Figure 27(a). As mentioned before, the multi-model Case B show an increase of the maximum peak value at x ¼ 15 mm, this is due to the transition zone between LW4 and LW1 models.
Regarding the transverse electric displacement D z , the single LW4 and all the multi-models, as shown in Figure 27(b), show a good description along with the inplane direction with some small oscillations in the zones close to the side-edge of the patch at x ¼ 5; 15mm. It has to be noticed that the linear single-model LW1 is completely not able to correct describe the transverse electric displacement, at x ¼ 5 mm the peak values show an inverse, positive, sign respect to the other single and multi-models, and at x ¼ 15 mm the maximum peak value is almost double respect to the other models.
Finally in Figure 28(a,b) the three-dimensional distributions of the electric potential Φ, obtained with the Abaqus 3D finite element C3D20RE and the present LW4 singlemodel, respectively are depicted on the entire plate structure, represented with an initial section at y ¼ 0 and the middle section at y ¼ b=2. It has to be noticed that the present LW4 single model well describes the phenomena without imposing any boundary condition, the electric potential tend naturally to zero.
The electric in-plane and transverse displacements D x and D z are depicted in Figure 29(a,b) and 30(a,b), respectively. The present LW4 single model and the Abaqus 3D finite element C3D20RE are in good agreement on the whole plate structure.

Conclusions
In this paper a new methodology for global/local analysis of composite and sandwich plate structure embedding piezoelectric skins and patches has been introduced. This approach makes use of advanced finite plate elements with node-dependent kinematics, which are formulated in the domain of the Unified Formulation. In fact, the Step: Step−1 Increment 1: Step finite element arrays of the generic plate element are formulated in terms of fundamental nuclei, which are invariants of the theory approximation order and the modeling technique (ESL, LW). In this manner, the plate theory can vary within the same finite elements with no difficulties. Thus, given a finite element model, the theory approximation accuracy can be enriched locally in a very straightforward manner by enforcing the same kinematics at the interface nodes between kinematically incompatible plate elements. The resulting global/local approach is very efficient because it does not employ any mathematical artifice to enforce the displacement/ stress continuity, such as those methods based on Lagrange multipliers or overlapping regions. Step: Step−1 Increment 1: Step  Step: Step−1 Increment 1: Step The present node-dependent variable kinematic model allows to locally improve the solution. Two main aspects can be highlighted: a reduction of computational costs with respect to Layer-Wise single-model solutions, and a simultaneous multi-models globallocal analysis can be performed in one-single analysis step. An accurate representation of secondary variables (mechanical stresses and electric displacements) in localized zones is possible with DOFs reduction if an accurate distribution of the higher-order kinematic capabilities is performed. On the contrary, the accuracy of the solution in terms of primary variables (mechanical displacements and electric potential) values depends on the global approximation over the whole structure. The efficacy of the node-dependent variable kinematic and global/local models, thus, depends on the characteristics of the problem under consideration as well as on the required analysis type. The proposed methodology has been widely assessed in this paper by analyzing composite and sandwich plates embedding piezoelectric skins and patches, in sensor and actuator configurations, and shear actuated configuration, and by comparison with solutions from the literature and those from finite element commercial tools. Future developments will deal with the extension of this global/local methodology to hierarchical shell theories and thermo-mechanical problems.