The Ritz – Sublaminate Generalized Unified Formulation approach for piezoelectric composite plates

ABSTRACT

This paper extends to composite plates including piezoelectric plies the variable kinematics plate modeling approach called Sublaminate Generalized Unified Formulation (SGUF). Two-dimensional plate equations are obtained upon defining a priori the through-thickness distribution of the displacement field and electric potential. According to SGUF, independent approximations can be adopted for the four components of these generalized displacements: an Equivalent Single Layer (ESL) or Layer-Wise (LW) description over an arbitrary group of plies constituting the composite plate (the sublaminate) and the polynomial order employed in each sublaminate. The solution of the two-dimensional equations is sought in weak form by means of a Ritz method. In this work, boundary functions are used in conjunction with the domain approximation expressed by an orthogonal basis spanned by Legendre polynomials. The proposed computational tool is capable to represent electroded surfaces with equipotentiality conditions. Free-vibration problems as well as static problems involving actuator and sensor configurations are addressed. Two case studies are presented, which demonstrate the high accuracy of the proposed Ritz-SGUF approach. A model assessment is proposed for showcasing to which extent the SGUF approach allows a reduction of the number of unknowns with a controlled impact on the accuracy of the result.

KEYWORDS:

• Piezoelectric composites
• variable kinematics model
• Ritz method
• unified formulation

1. Introduction

Composite structures with embedded or bonded piezoelectric sensors and/or actuators are among the most indicated solutions for implementing an adaptive strategy for, e.g., vibration suppression applications [1]. The design and implementation of such structures faces two contrasting requirements. On the one hand, it is very important for a reliable design to have the possibility of accurately predicting the global and local electromechanical response of the structure [2]. On the other hand, a computationally efficient procedure is necessary, for instance in the framework of optimization procedures [3] or in the control algorithm for allowing a real-time adaptation [4]. The three-dimensional (solid) elements with piezoelectric coupling currently available in commercial finite element software provide an excessively cumbersome computational framework, therefore much effort has been dedicated to the formulation of accurate reduced-order models, as witnessed by the early review papers [5,6] and by the more recent developments [7–12].

Unified formulations, pioneered by Reddy [13] and systematically developed by Carrera [14], allow constructing within a single computer program a variable-kinematics modeling approach with variable accuracy and computational cost, see also [15]. The importance of such an approach, in particular for coupled multi-field problems such as composite structures with piezoelectric elements, was early recognized [16–19]. Carrera Unified Formulation (CUF) was subsequently generalized by Demasi [20] (Generalized Unified Formulation, GUF). D’Ottavio [21] further extended the flexibility of GUF to include a sublaminate approach (SGUF): in SGUF, the composite stack is subdivided into an arbitrary number of laminates grouping adjacent physical plies, and a GUF variable-kinematics model can be employed for each sublaminate independently. Recently, the authors employed the displacement-based SGUF approach in conjunction with a Ritz method previously developed by Dozio and Carrera [22]: mechanical problems involving bending, buckling and free-vibration, including visco-elastic layers, have been addressed [23–25]. The proposed Ritz method is particularly useful as it provides sufficient flexibility for a rapid analysis of configurations involving arbitrary boundary and loading conditions as well as anisotropic coupling.

This paper presents the extension of the Ritz-SGUF approach to composite plates with piezoelectric plies working in extension mode (31-mode). The resulting computational tool is based on Hamilton’s principle involving the displacements and the electric potential as primary field variables. It accounts for the presence of electroded interfaces and verifies exactly the related equipotentiality conditions, which are known to play a major role for a realistic determination of the electro-mechanical coupling associated to the piezoelectric response [26].

The paper is organized as follows. The modeling approach is outlined in Section 2, where the compact index notation typical of Unified Formulations and the adopted Ritz solution are presented up to the derivation of the discrete algebraic system to be solved. Section 3 reports results obtained for representative case studies: the electro-mechanical coupling coefficients is investigated for various boundary conditions in Section 3.1, a static analysis with actuator and sensor configurations is carried out in Section 3.2, and the modeling flexibility of the SGUF approach is illustrated by referring to a double sandwich plate with piezoelectric plies in Section 3.3. Finally, Section 4 summarizes the main conclusions and suggests the direction for further work.

2. Modeling approach

2.1. Description of the geometry

A plate is considered that occupies the volume , where is the constant thickness and is the reference surface lying on the plane. The composite plate is composed of homogeneous plies stacked along the thickness direction , each of thickness , see Figure 1. According to the SGUF approach, the laminate is subdivided into sublaminates, each of constant thickness and consisting of a stack of adjacent plies. Throughout the paper, indexes and are related to the physical plies and the sublaminate, respectively, and quantities related to the ply of the sublaminate will be indicated by . Local coordinates and are introduced along with the corresponding non-dimensional coordinates

(1a)

Figure 1. SGUF: geometry description and employed coordinates across the thickness.

These coordinates are linked through the relation

(1b)

where and are the midplane coordinates of the ply and sublaminate, respectively.

2.2. The variational approach based on generalized displacements

The two-dimensional (2D) model of the composite piezoelectric plate is obtained upon introducing a priori assumptions for the through-the-thickness behavior of the displacements and the electrostatic potential in a variationally consistent manner. Unless otherwise stated, Latin indexes range in while Greek indexes take values in . Adopting tensor notation and the repeated index summation convention, the weak form of the mechanical and electrical equilibrium equations is expressed by the following integral statement

(2)

where , , and represent the strain and stress tensors and the electric field and electric displacement vectors, respectively, is the mass density, indicates an admissible virtual variation, is the vector of external traction imposed at the boundary and is the imposed electric charge density imposed at the boundary . In this paper it shall be assumed without loss of generality that mechanical load and electric charge densities are imposed at the top and/or bottom surfaces of the plate, viz. , while essential boundary conditions specifying displacements and electric potential are defined at . Eq.(2) can be hence re-expressed upon separating in-plane and through-the-thickness integrals and introducing the definition of sublaminates:

(3)

The variations in Eq.(2) and Eq.(3) are to be taken under the following subsidiary conditions:

• the linearized gradient relations for small perturbations:

(4)

where the notation stays for partial derivative with respect to .

• the piezoelectric constitutive equations in each ply:

(5)

where is the stiffness coefficient at constant electric field, the dielectric coefficient at constant strain and the piezoelectric stress coefficient. Elastic plies are represented by Eq.(5) with nil piezoelectric coupling, i.e., . Voigt compact notation is introduced based on the symmetry of the strain and stress tensors according to and [27]. Eq.(5) is thus expressed in terms of vectors and matrices as:

(6)

where uppercase Latin indexes . Note that Eq.(5) and Eq.(6) are considered to hold in the global Cartesian reference frame employed for describing the plate geometry, i.e., the relations include the angle describing the orientation of the principal material axis with respect to the structural axis , see also [27].

Note that the electric field is considered to be quasi-static compared to the mechanical wave propagation velocity.

2.3. Variable kinematics modeling in SGUF

Following the SGUF approach, each sublaminate is associated to a specific kinematic approximation, which is here extended for including the approximation for the electric potential . According to the compact index notation used in Unified Formulations, the generalized displacements are regrouped in one vector

(7)

In SGUF, the variation along the thickness direction of each unknown function is defined by a polynomial expansion of order and by specifying whether the plies of the sublaminate are described in an Equivalent Single Layer (ESL) or Layer-Wise (LW) manner. Therefore, the plate model is expressed in SGUF notation as

(8)

where for a displacement component and for the electric potential, see Eq.(7). The thickness functions are expressed in terms of the ply-coordinates if the plies of the sublaminate are described in a LW description; the approximation for the whole sublaminate is then obtained through the assembly of all plies upon enforcing the continuity of the function at the interface between adjacent plies. However, if an ESL description is adopted for the sublaminate , the thickness functions vary along and the ply index can be omitted in Eq.(8). If a constant approximation is used for a variable within a sublaminate, one has and . If a linear approximation is used, and the thickness functions correspond to the linear Lagrange polynomials that interpolate the values at the top and bottom of the ply (LW) or sublaminate (ESL). Higher-order approximations with are implemented through a hierarchic enrichment of the polynomial basis and employ Legendre polynomials. More details can be found elsewhere [21,23].

In a sublaminate whose plies are described in an ESL manner, a slope discontinuity can be introduced in the kinematics at the interface between adjacent plies by superposing a Zig-Zag function to the polynomial expansion given by Eq.(8). The Zig-Zag function employed in SGUF is a modification of the classical Murakami’s ZigZag Function (MZZF) [28] such that it vanishes at the top and bottom planes of the sublaminate:

(9)

This way, the continuity conditions to be imposed between adjacent sublaminates for building the matrices for the whole composite stack involve only the unknown displacements corresponding to the top and bottom of each sublaminate.

2.4. Ritz approximate solution

The weak form solution of the two-dimensional problem is sought within the space spanned by the global approximation (Ritz) functions, which can be expressed as

(10)

A map between the physical plane and the computational domain is introduced, where . Furthermore, separation of the in-plane variables is invoked thus expressing the th two-dimensional Ritz function as the product of two one-dimensional (‘beam’) functions:

(11)

where , (note that and may be, in general, different) and the relation between the th function and the beam functions and has been chosen as . Various sets of Ritz functions can be used; adopting orthogonal polynomials along with boundary functions allows to efficiently handle various boundary conditions, such as clamped, free or simply-supported edges, as well as to provide stable solutions for high values of [22]. This work employs Legendre polynomials for spanning the field solution in conjunction with boundary functions for imposing the essential boundary conditions at the edges:

(12a)

where for the employed functions are

(12b)

(12c)

The boundary conditions on the displacements can thus represent clamped, free or simply-supported edges, while those on the electric potential can describe a grounded edge () or an edge at which the potential is let free.

As far as the electric boundary conditions are concerned, the present implementation is also capable to cope with a prescribed distribution , which is imposed at an arbitrary interface . For this, the through-thickness model must be chosen so that a through-thickness DOF will be associated to the considered surface, which is straightforward within the sublaminate approach and a LW description: indeed, with reference to Eq.(8), one has to satisfy in order to have . For this variable, the Ritz approximation in Eq.(10) is then simply replaced by a one-term expression for the prescribed distribution . Electroded surfaces with a constant electric potential or a particular distribution of an applied electric potential are thus exactly represented as essential boundary conditions. It is worth noticing that only the essential boundary conditions are exactly satisfied by prescribing the values for the generalized displacements , while the conditions on the electro-mechanical fluxes (stresses and electric displacements) are let free to be satisfied in weak sense.

2.5. The discrete system

The weak form of the equilibrium equation is formally expressed in terms of the generalized displacements by substituting Eq.(4) into Eq.(6), and both into Eq.(3). The through-the-thickness approximation and the Ritz solution are then introduced for the virtual and actual field variables according to following notation:

(13)

where according to Eq.(7), and according to Eq.(8), and according to Eq.(10). All integrals along are explicitly carried out and are identified by the following notation:

(14)

In a similar manner, the in-plane integrals of the Ritz functions are evaluated according to the following notation:

(15)

The discrete expression of Eq.(3) can be finally written as

The left-hand side is the virtual internal work and includes the inertial terms, for which . is one of the material parameters of Eq.(6) depending on the nature of the generic variables and :

• if , i.e. and , then is an elastic stiffness coefficient;

• if and , i.e. and , then is a piezoelectric stress coefficient; analogously, if and , then

• if , i.e. , then is a dielectric permittivity coefficient.

The right-hand side represents the virtual work done by the external loading applied at given interfaces :

• for , a mechanical loading with a prescribed traction distribution :

• for , an electrical loading with a prescribed charge distribution

After cycling over all indexes of the thickness assumptions, assembly over all plies and sublaminates , and finally cycling over all Ritz expansion indexes , the following discrete system is obtained:

(17)

Electrical open-circuit and short-circuit conditions can be taken into account as detailed hereafter [29]. The electric potential DOF are subdivided in ‘internal’ electric DOF and those associated to an electroded surface, denoted . No electric charge is imposed in correspondance of the ‘internal’ electric potential DOF, which comprise the electric potential DOF at non-electroded (NE) surfaces also. These can hence be statically condensed out, which leads to the following discrete system:

(18a)

where

(18b)

In the short-circuit configuration (SC), one sets : the system to be solved is composed of only the first hyper-row of Eq.(18a) and its last hyper-row can be used to recover the electric charge vector within a post-processing step. In the open-circuit configuration (OC), one sets : a further static condensation of the electric potential DOF at the electrodes is carried out, which transforms the first hyper-row of Eq.(18a) into

(19)

The electric potential DOF at the open-circuited electrodes are then obtained from the second hyper-row in a post-processing step.

Static problems are obtained by dropping out the time-dependency, while free-vibration responses are obtained by setting along with and , where and is the eigenfrequency.

3. Numerical results

Before proceeding to the presentation of some numerical examples, the naming convention is introduced for identifying the variable kinematics models that can be implemented. For each sublaminate, the in-plane displacements , the transverse displacement and the electric potential can be described in an ESL (E) or LW (L) manner, and with a specific polynomial order , and . These informations are gathered for each sublaminate according to the following acronyms:

If a sublaminate does not host the electric potential variable , the corresponding index is dropped off. In fact, thanks to the sublaminate approach, the electric DOF can be retained only in those sublaminates with a piezoelectric response and omitted from those characterized by a merely elastic response.

For the whole stack of sublaminates, the acronyms are listed separated by a slash starting from the top sublaminate. For example, consider a composite layup subdivided into 3 sublaminates, all described by a FSDT kinematics; if the first (top) sublaminate is piezoelectric with a LW quadratic assumption for the electric potential, the second (central) sublaminate is purely elastic and the third (bottom) sublaminate is again piezoelectric with an ESL constant approximation for the electric potential, the resulting model is defined by the following acronym: .

3.1. Electromechanical coupling coefficients for a piezoelectric plate (T1)

This first case study, denoted as T1 and taken from [30], aims at illustrating the capabilities of the proposed approach to handle different electro-mechanical boundary conditions. In fact, these are known to play a relevant role in the possibility of the structure to convert electric and deformation energy by means of the electro-mechanical coupling. The first five bending modes are considered of a freely vibrating moderately thick plate () made of 2 equally thick plies of PZT-5 material (data in Table 1). The analysis is carried out in a two-dimensional plane strain setting. Four different sets of mechanical boundary conditions are considered at the plate edges , which can be both simply-supported (SS), both clamped (CC), clamped-free (CF) and clamped-simply-supported (CS). The simple support condition amounts to set . The interface between the two plies is electrically grounded and different electrical conditions are applied to the top and bottom surfaces of the plate, which may be covered by a thin conducting electrode or let without electrode. The electrode is modeled as an equipotential surface without mechanical stiffness at which two electric conditions can be applied: a nil electric potential is imposed, which yields a short-circuit condition (SC), or the potential value is let free, which yields an open-circuit condition (OC). If no electrode is present (NE), the electric potential is let free to vary along at the outer surfaces. In all cases, the vertical edges are electrode-free. Krommer provides analytical solutions obtained with an FSDT model (with a shear correction factor of 5/6) for the SC and OC configurations [30]. Further reference solutions are obtained by FEM computations with the commercial software Abaqus employing a mesh of 100 10 quadratic elements CPE8E.

Table 1. T1: Material data for the transversely isotropic PZT-5 material polarized along the direction. Superscripts and specify that the coefficients are given at constant electric field and at constant strain, respectively. Plate length [m], thickness [m].

[GPa]	121	[Cm^−2]	−5.46
[GPa]	75.4	[Cm^−2]	15.8
[GPa]	75.2	[Cm^−2]	12.32
[GPa]	111	[nFm^−1]	15.317
[GPa]	21.8	[nFm^−1]	15.052
[GPa]	22.8	[kgm^−3]	7750

Before proceeding further, it is useful to present a convergence study of the present Ritz solution. The first and the sixth circular eigenfrequencies are considered for a CC and a CF configuration. The electrical boundary conditions are the following: the mid-surface is electroded and grounded (), the bottom surface is electroded and in open-circuit ( free), the top surface is not electroded (NE). The results of the convergence study are reported in Figure 2 for a number of Ritz expansion terms that range from 4 up to 64 (cylindrical bending is assured by and ). The plate is modeled with an FSDT kinematics and a LW quadratic expansion of the electric potential, i.e., an model. A shear correction factor of 5/6 is used. Based on the convergence analysis, an expansion order will be used in the subsequent computations.

Figure 2. T1: Convergence analysis for first and sixth eigenfrequencies of the CC and CF plate.

Table 2 – Table 5 report the first 6 circular eigenfrequencies obtained by the present Ritz-SGUF approach for the four different mechanical boundary conditions. For each case, the three electrical configurations SC, OC and NE are considered. Two different models are used: the FSDT model (5/6 as shear correction factor) and the high-order model. Note that both models retain a quadratic approximation for the electric potential in order to well capture the electric field induced by the piezoelectric coupling. The eigenfrequencies obtained by the plane strain solid Abaqus model are in perfect agreement with those of the more accurate model model and are not reported in the tables for the sake of conciseness. Furthermore, the present model identically recovers the analytical results obtained with the FSDT model by Krommer, whose numerical values can be found in [30].

Table 2. T1: Circular eigenfrequencies [] for the hinged-hinged (SS) plate.

	elec BC	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
	SC	862.97	3279.3	6862.6	11,228	16,092	21,272
	OC	912.02	3279.3	6900.1	11,228	16,118	21,272
	NE	921.78	3456.8	7124.9	11,505	16,330	21,454
	SC	863.83	3291.1	6911.1	11,350	16,325	21,654
	OC	912.59	3291.1	6946.9	11,350	16,348	21,654
	NE	923.20	3475.2	7194.8	11,665	16,612	21,878

Table 3. T1: Circular eigenfrequencies [] for the clamped-clamped (CC) plate.

	elec BC	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
	SC	1849.3	4701.2	8444.0	12,766	17,484	22,480
	OC	1849.3	4701.2	8444.0	12,766	17,484	22,480
	NE	1938.1	4844.5	8599.4	12,903	17,597	22,579
	SC	1871.6	4767.6	8582.4	13,006	17,853	23,007
	OC	1871.6	4767.6	8582.5	13,006	17,853	23,007
	NE	1966.2	4935.2	8785.1	13,207	18,033	23,157

Table 4. T1: Circular eigenfrequencies [] for the clamped-free (CF) plate.

	elec BC	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
	SC	310.44	1849.6	4832.0	8708.1	13,189	18,060
	OC	323.73	1875.0	4854.4	8728.6	13,208	18,078
	NE	332.44	1952.9	5026.2	8937.0	13,405	18,240
	SC	312.11	1862.3	4876.5	8813.6	13,388	18,385
	OC	325.23	1888.4	4898.3	8835.1	13,406	18,404
	NE	334.03	1970.1	5091.7	9088.0	13,671	18,641

Table 5. T1: Circular eigenfrequencies [] for the clamped-hinged (CS) plate.

	elec BC	Mode 1	Mode 2	Mode 3	Mode 4	Mode 5	Mode 6
	SC	1313.8	3983.5	7662.0	12,012	16,803	21,888
	OC	1326.9	3994.1	7671.0	12,020	16,809	21,893
	NE	1390.6	4150.9	7873.7	12,218	16,975	22,025
	SC	1871.6	4767.6	9551.6	13,006	18,998	23,007
	OC	1871.6	4767.6	9551.6	13,006	18,999	23,007
	NE	1966.2	4935.2	10,478	13,207	20,826	23,157

Figure 3 reports the ratio for the first 6 bending modes, where is the eigenfrequency of the elastic plate when the piezoelectric coupling is neglected and refers to the different configurations SC, OC and NE. Figure 4 illustrates the modal electro-mechanical coupling coefficient EMCC, defined as , where for electroded (equipotential) surfaces and if a electric potential distribution can develop at the outer surfaces. A perfect match is obtained between the present results and those obtained by Krommer [30] and with the solid Abaqus model. It is interesting to note that, despite the rather large differences between the eigenfrequencies computed by the FSDT and kinematics, see Table 2–5, Figures 3 and 4 show that the discrepancies in the frequency ratios and are small for the NE configuration and practically nil for the SC and OC configurations.

Figure 3. T1: Ratio between the eigenfrequencies of the piezoelectric plate and those of the equivalent elastic plate for different mechanical boundary conditions and electric circuits.

Figure 4. T1: Modal Electro-Mechanical Coupling Coefficient for different mechanical boundary conditions with and without equipotential electrode conditions.

These results confirm the important role of the electro-mechanical boundary conditions on the coupled response of piezoelectric structures and the following comments can be made. Figure 3 shows the stiffening effect induced by the piezoelectric coupling for all considered electric boundary conditions. The largest increase is of about 10% and is obtained by the electrode-free configuration (NE configuration), while the SC configuration leads to the smallest increase. The stiffening effect of the SC configuration constantly increases from about 2.5% for mode 1 to about 4.5% for mode 6, irrespective of the mechanical boundary conditions. On the contrary, the stiffening effect associated to the NE and OC has its maximum for mode 1, with the NE configuration being alway stiffer than the OC configuration. As the mode number increases, the OC curve merge the SC curve, whereas the NE curve provides always slightly higher eigenfrequencies. The EMCC reported in Figure 4 directly quantifies the piezoelectric coupling effect. It is obvious that the NE configuration induces always a larger piezoelectric coupling than the OC configuration, an effect explained by the charge cancellation effect at the equipotential surface [31]. The EMCC is maximum for the first mode and decreases for higher modes. Furthermore, statically determined (SS and CF) plates show a higher EMCC than those with more rigid end conditions (CC and CS). For the OC configuration (equipotential electrodes), the EMCC of the SS plate vanishes for modes with an even number of half-waves, and it is identically nil for the CC plate.

3.2. Static analysis of an orthotropic laminate with a piezoelectric ply (T2)

This second case study, identified as T2, is particularly interesting because Vel and Batra provided in [32] reference 3D solutions for static actuator and sensor problems involving electro-mechanical boundary conditions different from the classical simply supported and electrically grounded edges.

The problem deals with a thick, square plate ( [m], ) consisting of a PZT-5A ply (thickness ) bonded on top of an elastic substrate composed of two Graphite-Epoxy (GrEp) plies (thickness each). The non-zero coefficients of these materials are reported in Table 6. The fibers of the bottom GrEp ply are aligned with the direction while those of the top GrEp ply are aligned with the direction. The two edges are simply-supported () and electrically grounded (); at the edges the two following sets of electro-mechanical boundary conditions are considered:

Table 6. T2: Material data for PZT-5A (transverse isotropy in plane) and for GrEp (transverse isotropy in plane). Superscripts and specify that the coefficients are given at constant electric field and at constant strain, respectively.

		PZT-5A			GrEp
[GPa]	99.201	[]	−7.209	[GPa]	183.443	[]	0
[GPa]	54.016	[]	15.118	[GPa]	4.363	[]	0
[GPa]	50.778	[]	12.322	[GPa]	3.918	[]	0
[GPa]	86.856	[]	15.30	[GPa]	11.662	[]	15.30
[GPa]	21.100	[]	15.00	[GPa]	2.870	[]	15.30
[GPa]	22.593	[]	7750	[GPa]	7.170	[]	1590

SP configuration:

CD configuration:

Therefore, the edges are simply-supported and electrically grounded in the SP configuration and are clamped and electrically insulated in the CD configuration.

The interface between the elastic substrate and the piezoelectric ply is electroded and grounded while the bottom surface of the plate is traction-free. A sensor and an actuator configuration are considered, for which the boundary conditions at the top surface and the non-dimensional output quantities are defined as follows:

sensor configuration:

(20)

actuator configuration:

(21)

The coefficients employed for the non-dimensional quantities defined above are [GPa], []

Table 7 summarizes the reference values reported in [32] with present solutions obtained by a quasi-3D model adopting an ED model for the piezoelectric ply, an LD model for the elastic substrate and Ritz expansion orders . The non-dimensional local response quantities defined in Eq.(20) (sensor configuration) and in Eq.(21) (actuator configuration) are reported for both, the SP and CD boundary conditions. For the SP boundary conditions, the present Ritz solution converges very rapidly (the reported values were already obtained with ) and exactly matches the analytical 3D solution. For the CD configuration the convergence is definitively slower, as already remarked by Vel and Batra [32], due to the difficulty in satisfying the electric boundary condition of the insulated edges of the piezoelectric ply. Nevertheless, the present results are shown to well agree with the exact 3D solution even in the CD configuration.

Table 7. T2: Comparison between local non-dimensional response obtained with the quasi-3D model ED/LD with and the analytical 3D solution.

	Sensor configuration				Actuator configuration
	SP configuration		CD configuration		SP configuration		CD configuration
	present	3D [32]	present	3D [32]	present	3D [32]	present	3D [32]
	−1.933	−1.933	−1.081	−1.082	10.160	10.161	4.773	4.774
	14.325	14.325	10.849	10.851	−25.861	−25.862	−14.198	−14.205
	9.329	9.329	6.683	6.652	15.517	15.517	11.114	11.115
	0.972	0.972	1.068	1.056	1.042	1.042	−0.021	0.057
	0.384	0.384	0.095	0.096	0.836	0.836	2.184	2.184
	0.419	0.419	0.462	0.462	−0.119	−0.119	−0.011	−0.012
	−3.668	−3.668	−3.019	−3.020	0.505	0.505	0.503	0.502
	21.563	21.563	13.240	13.301	−9.878	−9.878	−9.431	−9.440

The influence of electro-mechanical boundary conditions is illustrated in Figure 5 for the sensor configuration and in Figure 6 for the actuator configuration. The through-the-thickness distributions are reported of the non-dimensional in-plane and transverse displacements and (Figure 5(a) and 6(a)), as well as bending and shear stresses and (Figure 5(b) and 6(b)). The results are obtained with the quasi-3D model ED/LD with .

Figure 5. T2-Sensor configuration: influence of the boundary conditions on the mechanical response of the PZT/GrEp laminate subjected to an external pressure load.

Figure 6. T2-Actuator configuration: influence of the boundary conditions on the mechanical response of the PZT/GrEp laminate under the action of an electric potential difference in the PZT ply.

Under the action of an external pressure load (sensor case), the in-plane displacement shows a strong zig-zag distribution, which is more accentuated for the CD boundary conditions than for the SP boundary conditions, see Figure 5(a) left. Due to the clamped edges, the CD configuration reduces the deflection of the plate (Figure 5(a) right) and, hence, the bending stress (Figure 5(b) left). This stress is largest at the outer surfaces . The transverse shear stress is quadratic in each ply and shows its maximum value within the bottom GrEp ply, whose fibers are oriented along the axis (Figure 5 right). This stress is approximately 1.5 times larger in the CD configuration than in the SP configuration.

In the actuator configuration, the zig-zag shape of the in-plane displacement is barely appreciable, see Figure 6(a) left. Figure 6(a) right shows that the transverse displacement is clearly non-constant across the plate thickness, with a rather steep linear variation across the piezoelectric actuator. Obviously, the clamped edges in the CD configuration reduce the deflection induced by the actuation. The through-thickness distribution of the bending stress , illustrated in Figure 6(b) left, is sensibly different from that of the sensor configuration: on the one hand, it shows its maximum values at the interfaces between the outer plies and the central GrEp ply, and not at the outer surfaces; on the other hand, the bending stress in the piezoelectric actuator is larger in the CD configuration than in the SP configuration. The transverse shear stress is again quadratic in each ply and has its largest value within the bottom GrEp ply, but it shows an important peak also at the interface between the PZT ply and the elastic substrate.

A final comment is made concerning the fulfilment of the transverse stress boundary conditions at the top/bottom faces of the plate as well as their continuity at interfaces between dissimilar material plies, see the rightmost graphics of Figure 5(b) and 6(b). Despite the present model is a displacement-based one, the high order expansion employed in each material ply (polynomials of degree 6) yields in fact very accurate solutions even for the transverse stress field, including the fulfilment of the natural boundary conditions. In order to satisfy these conditions with low-order models, one may consider the use of (partially) mixed formulations, see, e.g. [33,34].

3.3. Free-vibration analysis of a double sandwich plate with bonded piezoelectric plies (T3)

The third test case (T3) addresses the free-vibration response of an elastic substrate consisting of a double sandwich plate with two piezoelectric plies bonded at the top and bottom surfaces. The plate is square, of total thickness and simply-supported at its four edges. Figure 7(a) illustrates the composite stack and the electric boundary conditions: the double sandwich consists of two cores sandwiched between two composite laminates and separated by a third composite laminate; each composite laminate has three plies with stacking sequence [0/90/0]; the top surface as well as the interfaces between the elastic substrate and the piezoelectric plies are electroded and grounded, whereas the bottom surface is not electroded. The properties of the piezoelectric material (PZT-5A) correspond to those listed in Table 6, the mechanical properties of the core and face ply materials are listed in Table 8. This configuration, taken from Kapuria and Nath [35], is here employed for highlighting how the proposed SGUF approach can optimize the number of DOF for reaching a desired accuracy. For this, models for the whole stack will be first proposed without resorting to the definition of sublaminates (i.e., only one sublaminate represents the whole stack) and then by subdividing the stack into 3 and 7 sublaminates, as indicated in Figure 7(b and c), respectively.

Table 8. T3: Material data for the face and core materials.

	face ply	core
[GPa]	133.931	1.976
[GPa]	4.423	1.773
[GPa]	4.423	1.475
[GPa]	9.226	1.791
[GPa]	4.595	1.406
[GPa]	9.226	2.761
[GPa]	2.332	0.4554
[GPa]	2.332	0.5451
[GPa]	3.558	0.01656
[]	1000	70

Figure 7. T3: (a) double sandwich plate with bonded piezoelectric plies; (b) subdivision in sublaminates; (c) subdivision in sublaminates.

The results are reported in terms of non-dimensional circular eigenfrequencies

(22)

A thick and a moderately thin plate are considered and the attention is limited to the flexural modes (1,1), (2,1) and (3,1), which correspond to vibration modes with 1, 2 and 3 halfwaves along the direction, respectively. Ritz expansion orders are used which ensure converged solutions. Table 9 and Table 10 report the percentage difference between the 3D solutions, denoted by and provided in [35], and those obtained by some models of the proposed variable kinematics approach.

Table 9. T3: assessment of models with one, three and seven sublaminates for a thick plate .

		mode (1,1)		mode (2,1)		mode (3,1)
[–]		5.6641		9.9708		14.963
1 sublaminate ()
model	DOF	[–]	[%]	[–]	[%]	[–]	[%]
	134	5.6641	0.00	9.9708	0.00	14.964	0.01
	108	5.6641	0.00	9.9709	0.00	14.964	0.01
	56	5.6437	−0.01	10.015	0.44	15.130	1.11
	41	6.8754	21.82	12.452	24.88	18.672	24.79
	38	6.9477	23.10	12.662	26.99	19.028	27.17
3 sublaminates (): /SL2/
SL2	DOF	[–]	[%]	[–]	[%]	[–]	[%]
	100	5.6393	−0.09	9.9587	−0.12	14.938	−0.16
	34	5.6430	−0.02	10.012	0.42	15.123	1.07
	23	6.7114	18.91	11.901	19.36	17.528	17.14
	20	6.8407	21.20	12.340	23.76	18.267	22.08
7 sublaminates (): /SL-F/SL-C/SL-F/SL-C/SL-F/
SL-F/SL-C	DOF	[–]	[%]	[–]	[%]	[–]	[%]
	32	5.6347	−0.17	9.9661	−0.05	14.987	0.16
	26	5.6404	−0.07	9.9850	0.14	15.039	0.51
	22	5.6495	0.09	10.032	0.61	15.176	1.42

Table 10. T3: assessment of models with with one, three and seven sublaminates for a moderately thin plate .

		mode (1,1)		mode (2,1)		mode (3,1)
[–]		11.059		26.013		44.708
1 sublaminate ()
model	DOF	[–]	[%]	[–]	[%]	[–]	[%]
	134	11.059	0.00	26.013	0.00	44.712	0.01
	108	11.059	0.00	26.013	0.00	44.712	0.01
	56	11.053	−0.06	25.989	−0.09	44.671	−0.08
	41	11.574	4.65	28.572	9.84	51.381	14.93
	38	11.599	4.88	28.808	10.75	52.098	16.53
3 sublaminates (): /SL2/
SL2	DOF	[–]	[%]	[–]	[%]	[–]	[%]
	100	11.057	−0.02	26.005	−0.03	44.691	−0.04
	34	11.052	−0.06	25.989	−0.09	44.669	−0.09
	23	11.513	4.11	28.242	8.57	50.450	12.84
	20	11.556	4.50	28.646	10.12	51.667	15.57
7 sublaminates (): /SL-F/SL-C/SL-F/SL-C/SL-F/
SL-F/SL-C	DOF	[–]	[%]	[–]	[%]	[–]	[%]
	32	11.053	−0.06	25.988	−0.01	44.657	−0.11
	26	11.055	−0.03	26.000	−0.05	44.690	−0.04
	22	11.055	−0.03	26.004	−0.04	44.709	0.00

By using only one sublaminate, , the electric potential is always LW in order to meet the electric boundary conditions, while the kinematics may be chosen as LW or ESL. The high-order LW model and the classical CUF model LD2 both provide the exact 3D solution, but these models need at least 100 DOF. Adopting the FSDT kinematics in each ply reduces the number of DOF to 56; in this case, an error larger than 1% is found only for the shortest wavelength (3,1) mode of the thick plate (. Models that rely upon an ESL kinematics for the whole composite stack are seen to be inadequate, even if the Zig-Zag function is introduced. In fact, the Zig-Zag function contributes only marginally because it is incapable of coping with the different interfaces that characterize the stack, e.g., the interface between the PZT-5A ply and the face ply and that between the face ply and the core. It is worth noticing that models adopting the plane stress condition for the piezoelectric constitutive law, which is the case for the FSDT kinematics, may yield lower eigenfrequencies than the 3D reference. In fact, low-order kinematics in conjunction with the plane stress assumption comes along with an underestimation of the electrical energy induced by the electro-mechanical coupling, which in turn results in an overly compliant response of the piezoelectric structure. Similar results can be found, e.g., in [35,36].

Models with sublaminates are defined upon regrouping all plies constituting the double sandwich elastic plate into the second sublaminate (SL2), see Figure 7(b). The results in Tables 9 and 10 are thus obtained with several representative models for the double sandwich structure SL2. The piezoelectric plies are modeled with an FSDT kinematics in conjunction with an electric potential distribution that is quadratic for the short-circuited ply at the top (SL1) and linear for the open-circuited ply at the bottom (SL3). With respect to the case with only one sublaminate, the number of DOF is already reduced by because no electric potential DOF are associated to the elastic basis structure. Moreover, since the piezoelectric plies and the SL2 are assembled in an LW sense, the accuracy is in general higher than that of models with only one sublaminate.

An even more efficient modeling strategy consists in regrouping into different sublaminates the two piezoelectric plies, the two core plies and the three laminated faces of the sandwich structure, see Figure 7(c). In Tables 9 and 10, this representation involving 7 sublaminates is identified by specifying the sublaminate models for the cores (SL3 = SL5 = SL-C) and those for the laminated faces (SL2 = SL4 = SL6 = SL-F). The models and are again employed for the sublaminates SL1 and SL7 corresponding to the top and bottom piezoelectric plies, respectively. In general, it can be seen that the LW assembly of sublaminates allows to reduce the order of the approximations inside each sublaminate without affecting the accuracy. So, a simple FSDT kinematics with or without Zig-Zag function for the in-plane displacements can be used for the laminated faces. The core model may be based on FSDT also, or it may be enhanced upon retaining the through-thickness stretch ( model for SL-C). This latter effect is seen to play a certain role only in the (3,1) mode of the thick plate (. If an error of 1.4% is acceptable for this mode, the less expensive model has only 22 DOF, otherwise a model with 26 DOF can be used, which ensures a maximum error of 0.5% for all considered modes. This level of accuracy can be achieved at this relatively low number of DOF thanks to the subdivision of the whole stack into several sublaminates.

In closure, the number of DOF required by the present Ritz-SGUF approach is compared with respect to a solid FEM approach available in commercial packages for the piezoelectric analysis. Only the number of DOF along the thickness coordinate shall be here considered, the discretization employed to resolve the gradients in the plane won’t be taken into account. Nevertheless, it is worth recalling that the FE mesh should avoid extreme aspect ratios of the 3D elements in order to ensure a proper conditioning, which means that the in-plane mesh density should depend on the number of elements employed across the thickness. In order to capture the piezoelectric coupling inside the top ply, at least 2 linear finite elements are required across the thickness for approximating the electric potential, i.e., 3 electric potential DOF. Since in solid elements the same approximation is used for all unknowns, a total of 12 DOF should be used across the thickness of the piezoelectric ply. On the contrary, since the present variable kinematics approach can introduce different approximations for the field variables, the employed model ED uses only 8 DOF for the piezoelectric ply. Similar considerations hold for all models in which the order of expansion is different for the different field variables. It is finally mentioned that some commercial FEM packages allow to model sublaminates by means of laminated solid-shell elements. In this case, only FSDT kinematics can be used and no piezoelectric coupling is present.

4. Conclusions and outlook

A variable-kinematics modeling approach based on the Sublaminate Generalized Unified Formulation (SGUF) has been extended to composite plates with bonded piezoelectric plies working in extension mode. Ritz method has been used to solve in weak form the resulting two-dimensional differential equations, which may include various types of essential boundary conditions, in particular equipotential surfaces. The proposed results have demonstrated that the resulting SGUF-Ritz approach is a valuable tool for investigating piezoelectric plate problems, including static actuator or sensor applications and free-vibration analysis with various electro-mechanical boundary conditions; assessing a large number of models in view of a minimization of the number of DOF with a controlled accuracy, i.e., for identifying the best compromise between computational cost and results’ accuracy. As a general guideline, the accuracy is greatly enhanced by adopting separate sublaminates for groups of plies characterised by a strong mismatch with respect to adjacent plies, i.e., by using separate sublaminates for the piezoelectric plies, the stiff composite laminates and the soft core plies.

The effect of model assumptions on the accuracy can be quantified in a systematic manner with respect to the associated computational cost by referring to the so-called mixed axiomatic/asymptotic approach [37–39]: Best Theory Diagrams have been derived in the framework of analytical Navier-type and FEM solutions for sample problems involving homogeneous metallic plates as well as laminated and sandwich plates. This technique has been applied to analytical Navier-type solutions for piezoelectric plates in [40]. The implementation of this methodology on the basis of the Ritz-SGUF approach could be an interesting subject for a further study.

Further work will be dedicated to the inclusion of viscoelastic elements for the design and analysis of integrated active-passive vibration damping systems. A finite element implementation of the SGUF approach with a consistent global-local submodeling approach related to the 2D mesh shall allow an efficient investigation of more complex configurations, e.g. realistic structures with multiple active and/or passive damping patches.

Disclosure statement

No potential conflict of interest was reported by the authors.