Layer-Wise Discontinuous Galerkin Methods for Piezoelectric Laminates

Abstract

In this work, a novel high-order formulation for multilayered piezoelectric plates based on the combination of variable-order interior penalty discontinuous Galerkin methods and general layer-wise plate theories is presented, implemented and tested. The key feature of the formulation is the possibility to tune the order of the basis functions in both the in-plane approximation and the through-the-thickness expansion of the primary variables, namely displacements and electric potential. The results obtained from the application to the considered test cases show accuracy and robustness, thus confirming the developed technique as a supplementary computational tool for the analysis and design of smart laminated devices.

1. Introduction

Multilayer composite beams, plates and shells are today widely employed in several fields of engineering, for example in the aerospace, naval, automotive and biomedical sectors, among others [1]. The reason for their success is largely related to the possibility of achieving excellent stiffness-to-weight and strength-to weight ratios and especially to the possibility of tailoring such properties with respect to the application features, e.g., following the load paths, so as to increase the overall structural efficiency in terms of weight.

The availability of techniques for manufacturing composite laminates employing layers with multi-physics properties has led to the concept of smart laminates, in which different laminae may exhibit couplings between different physical fields, e.g., mechanical, thermal, electric, magnetic or even a combination of more than two of such fields, thus enabling the design of devices with multifunctional capabilities, besides the basic structural employment [2,3,4,5]. A remarkable example is provided by piezoelectric materials and laminates that, exploiting the coupling between the mechanical and electrical fields, may be employed in several applications of engineering or scientific interest [6,7,8,9].

The width of the available design space for this kind of multi-physics multilayer materials, deriving from the variety of properties of the individual laminae and from the layup flexibility, induces the need of a preliminary careful assessment of the properties of the assembled laminate: their heterogeneous nature, the mismatch of material features between contiguous layers and the variety of external loading condition may indeed result in complex multi-physics states that, if not carefully designed, may lead to unwanted behaviors or unforeseen damage/failure patterns.

For such a reason, several analytical, numerical and computational methods have been developed for the analysis of multilayer composites. The reconstruction of the fields (mechanical, thermal, electrical, etc.) within the laminate is a fully three-dimensional problem for which analytical solutions are available only for basic configurations: in general engineering applications, the employment of numerical/computational methods is then necessary. Besides fully three-dimensional approaches, which may be employed to capture complex localized field details (see, e.g., [10]), thin and relatively thick multilayered composite structures have been investigated using suitable plate theories, which are based on specific assumptions about the behavior of some specific field components throughout the laminate thickness and allow reducing the total number of unknowns and thus the computational cost associated with the numerical models: familiar examples are provided by the Kirchhoff and Mindlin kinematic assumptions in the case of purely structural plates, while also the coupling with the electric fields must be considered in piezoelectric plates (see, e.g., [1,11]).

The structural theories for laminated plates may be classified into equivalent-single-layer (ESL) and layer-wise (LW) models, which have been used for both elastic and multi-physics plates [12,13,14,15]. ESL models are based on an initial expression of some unknown fields, e.g., displacements in elastic plates or displacements and electric potential in the case of piezoelectric plates, in terms of a sum of individual contributions given by the products of known through-the-thickness functions and unknown in-plane functions, where the through-the-thickness functions span the whole thickness of the laminate. LW models differ from ESL ones for the adoption of through-the-thickness functions restrained to individual layers, which must be then complemented by suitable inter-laminar continuity conditions. While ESL models are suitable for capturing global behaviors of the laminate, LW models provide a more accurate representation of the fields on a layer-by-layer basis.

Once an ESL or LW representation of the multi-physics behavior of the structural element is assumed, the reduced equations, expressed in terms of the unknown in-plane functions, must be numerically solved. A popular strategy for the numerical treatment of the problem is based on the employment of a suitable multi-physics, or generalized, variational statement, providing the weak formulation of the problem, followed by the application of a suitable solution technique, e.g., the finite element method [13], the Rayleigh–Ritz method [16,17,18] or other similar approaches.

In FEM-based methods, the unknown fields are approximated over each finite element in terms of shape functions and nodal values and inter-element continuity is ensured by enforcing that contiguous elements share the same nodal values. Although it appears quite natural, this aspect may represent a limiting constraint when high-order, adaptive or non-conforming meshes are adopted. More advanced FE-based methods, which have been developed also for piezo-electric plates and structures, bypass such issues and allow retrieving discontinuous kinematics by resorting to node-dependent fields expansions, in which the generalized structural theory is a function of the node and may vary within the element itself (see, e.g., [19,20,21,22]).

An alternative approach to a continuous approximation of the solution, and to the mentioned advanced FE strategies with node-dependent kinematics, is provided by the discontinuous Galerkin (dG) method [23], which admits a discontinuous interpolation of the solution among the mesh elements and naturally handles arbitrary meshes, high-order basis functions and contiguous mesh elements with different order of approximation. When compared to the FEM literature, the adoption of dG-based numerical schemes for multilayered plates appears quite limited and mostly focused on the Classical Laminated Plate Theory [24,25,26,27] or the First-order Shear Deformation Theory [28,29]. Recently, Gulizzi et al. [30,31] presented dG formulations for ESL and LW theories for elastic plates, whereas, to the best of the Authors’ knowledge, no discontinuous Galerkin formulation for LW theories of piezoelectric plates is available.

In this work, we thus present an interior penalty discontinuous Galerkin formulation for multilayered piezoelectric plates based on a suitable generalization of the Principle of Virtual Displacements, which allows a high-order resolution of the electro-elastic problem. Throughout the thickness, the plate is represented adopting a LW representation, which naturally allows selecting a high-order expansion of the unknown fields within each layer, whereas, within the plate domain, the discontinuous nature of dG schemes allows employing general high-order mesh elements.

The paper is organized as follows. The basic equations of the problem in hand are first recalled in Section 2, which reviews the generalized formulation of piezoelectricity, the layer-wise formulation for piezoelectric laminates and the related laminates governing equations. The proposed method is formulated in Section 3. Numerical results on single and multiple layers piezoelectric plates are provided by Section 4, which is followed by the Conclusions.

2. Boundary-Value Problem

In this section, the basic equations describing the mechanics of general piezoelectric laminated plates are recalled. Figure 1 shows a multilayered plate obtained by stacking $N_{\ell }$ suitably oriented piezoelectric layers. In the figure, ${\theta }^{\langle \ell \rangle }$ denotes the angle between the ℓth layer’s local reference system ${\tilde{X}}_1^{\langle \ell \rangle }{\tilde{X}}_2^{\langle \ell \rangle }{\tilde{X}}_3^{\langle \ell \rangle }$ and the plate global reference system $X_1{X}_2{X}_3$.

The plate lies in the domain $V\equiv \mathsf{\Omega }\times [z_b,z_t]$, where $z_b$ identifies the height of the plate bottom surface, $z_t$ is the height of its top surface and $\mathsf{\Omega }$ denotes the modelling domain of the plate spanned by the coordinates $(X_1,X_2)$. The layers are stacked in such a way that $z_t^{\langle \ell \rangle }=z_b^{\langle \ell +1\rangle }$, for $\ell =1,\cdots ,N_{\ell }-1$, $z_b=z_b^{\langle 1\rangle }$ and $z_t=z_t^{\langle N_{\ell }\rangle }$, being $z_b^{\langle \ell \rangle }$ and $z_t^{\langle \ell \rangle }$ the coordinates of the ℓth layer’s bottom and top surfaces, respectively. Henceforth, quantities related to the layer ℓ are denoted by the superscript $\langle \ell \rangle$.

The electro-mechanical behavior of the multilayered piezoelectric plates is described in terms of mechanical displacements ${\mathit{d}}^{\langle \ell \rangle }={\{d_1,d_2,d_3\}}^{\langle \ell \rangle }$ and electric potential ${\varphi }^{\langle \ell \rangle }$, $\forall \ell =1,\cdots ,N_{\ell }$, which represent the primary variables of the problem.

The gradients of the primary variables, namely the strains ${\mathit{\gamma }}^{\langle \ell \rangle }\equiv {\{{\gamma }_{11},{\gamma }_{22},{\gamma }_{33},{\gamma }_{23},{\gamma }_{13},{\gamma }_{12}\}}^{\langle \ell \rangle }$ and the electric field ${\mathit{E}}^{\langle \ell \rangle }\equiv {\{E_1,E_2,E_3\}}^{\langle \ell \rangle }$, may be expressed as

$${\mathit{\gamma }}^{\langle \ell \rangle }={\mathit{I}}_{\alpha }^{\gamma }\frac{{\partial \mathit{d}}^{\langle \ell \rangle }}{{\partial X}_{\alpha }}+{\mathit{I}}_3^{\gamma }\frac{{\partial \mathit{d}}^{\langle \ell \rangle }}{{\partial X}_3}\mathrm{and}{\mathit{E}}^{\langle \ell \rangle }=-{\mathit{I}}_{\alpha }^E\frac{{\partial \varphi }^{\langle \ell \rangle }}{{\partial X}_{\alpha }}-{\mathit{I}}_3^E\frac{{\partial \varphi }^{\langle \ell \rangle }}{{\partial X}_3},$$

(1)

where the matrices ${\mathit{I}}_1^{\gamma }$, ${\mathit{I}}_2^{\gamma }$, ${\mathit{I}}_3^{\gamma }$, ${\mathit{I}}_1^E$, ${\mathit{I}}_2^E$, ${\mathit{I}}_3^E$ are defined as

$${\mathit{I}}_1^{\gamma }\equiv \left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right],{\mathit{I}}_2^{\gamma }\equiv \left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\\ 0& 0& 1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right],{\mathit{I}}_3^{\gamma }\equiv \left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\\ 0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right]$$

(2)

and

$${\mathit{I}}_1^E\equiv \left[\begin{array}{c}1\\ 0\\ 0\end{array}\right],{\mathit{I}}_2^E\equiv \left[\begin{array}{c}0\\ 1\\ 0\end{array}\right],{\mathit{I}}_3^E\equiv \left[\begin{array}{c}0\\ 0\\ 1\end{array}\right].$$

(3)

It is worth noting that, in Equation (1) and the continuation of the paper, Greek subscripts take value in $\{1,2\}$ and the Einstein’s implicit summation convention over repeated subscripts is assumed.

The problem description is completed by assigning the constitutive relation to its piezoelectric layers, which are assumed to be homogeneous and governed by linear piezoelectricity. Then, the constitutive behavior of the ℓth layer can be expressed as

$${\mathit{\sigma }}^{\langle \ell \rangle }={\mathit{c}}^{\langle \ell \rangle }{\mathit{\gamma }}^{\langle \ell \rangle }-{\mathit{e}}^{\langle \ell \rangle \top }{\mathit{E}}^{\langle \ell \rangle }\mathrm{and}{\mathit{D}}^{\langle \ell \rangle }={\mathit{e}}^{\langle \ell \rangle }{\mathit{\gamma }}^{\langle \ell \rangle }+{\mathit{\kappa }}^{\langle \ell \rangle }{\mathit{E}}^{\langle \ell \rangle }$$

(4)

where ${\mathit{\sigma }}^{\langle \ell \rangle }$ collects the stress components in the same order as that used for the strain components in ${\mathit{\gamma }}^{\langle \ell \rangle }$; ${\mathit{D}}^{\langle \ell \rangle }$ collects the components of the electric displacement in the same order as that used for the electric field components in ${\mathit{E}}^{\langle \ell \rangle }$; and the matrices ${\mathit{c}}^{\langle \ell \rangle }$, ${\mathit{\kappa }}^{\langle \ell \rangle }$ and ${\mathit{e}}^{\langle \ell \rangle }$ denote, respectively, the elasticity, dielectric and piezoelectric constitutive matrices of the ℓth layer. Equations (1) and (4) can be straightforwardly combined to link the stress ${\mathit{\sigma }}^{\langle \ell \rangle }$ and electric displacement ${\mathit{D}}^{\langle \ell \rangle }$ with the primary variables. It is worth noting that, when the piezoelectric constants are set to zero, decoupled elastic and dielectric behaviors are retrieved.

2.1. Generalized Formulation for Piezoelectricity

The form of Equations (1)–(4) suggests the introduction of generalized quantities (see, e.g., [32,33,34]). More specifically, for the sake of the present formulation, it is convenient to define the generalized displacement vector ${\mathit{u}}^{\langle \ell \rangle }$, the generalized strain vector ${\mathit{ϵ}}^{\langle \ell \rangle }$ and the generalized stress vector ${\mathit{\tau }}^{\langle \ell \rangle }$ as

$${\mathit{u}}^{\langle \ell \rangle }\equiv {\left\{\begin{array}{c}\mathit{d}\\ \varphi \end{array}\right\}}^{\langle \ell \rangle },{\mathit{ϵ}}^{\langle \ell \rangle }\equiv {\left\{\begin{array}{c}\mathit{\gamma }\\ -\mathit{E}\end{array}\right\}}^{\langle \ell \rangle }\mathrm{and}{\mathit{\tau }}^{\langle \ell \rangle }\equiv {\left\{\begin{array}{c}\mathit{\sigma }\\ \mathit{D}\end{array}\right\}}^{\langle \ell \rangle }.$$

(5)

Then, using Equation (5), Equation (1) can be rewritten as

$${\mathit{ϵ}}^{\langle \ell \rangle }={\mathit{I}}_{\alpha }\frac{{\partial \mathit{u}}^{\langle \ell \rangle }}{{\partial X}_{\alpha }}+{\mathit{I}}_3\frac{{\partial \mathit{u}}^{\langle \ell \rangle }}{{\partial X}_3},\mathrm{with}{\mathit{I}}_k\equiv \left[\begin{array}{cc}{\mathit{I}}_k^{\gamma }& \mathbf{0}\\ \mathbf{0}& {\mathit{I}}_k^E\end{array}\right],k=1,2,3,$$

(6)

while the constitutive relationships given in Equation (4) become

$${\mathit{\tau }}^{\langle \ell \rangle }={\mathit{p}}^{\langle \ell \rangle }{\mathit{ϵ}}^{\langle \ell \rangle },\mathrm{where}{\mathit{p}}^{\langle \ell \rangle }\equiv {\left[\begin{array}{cc}\mathit{c}& {\mathit{e}}^{\top }\\ \mathit{e}& -\mathit{\kappa }\end{array}\right]}^{\langle \ell \rangle }.$$

(7)

2.2. Layer-Wise Formulation for Multilayered Plates

Layer-wise formulations for multilayered plates are based on the introduction of suitable assumptions regarding the through-the-thickness behavior of the primary variables [1,13]. In this work, the primary variables are collected in the vector ${\mathit{u}}^{\langle \ell \rangle }$, $\ell =1,\cdots ,N_{\ell }$, which is assumed to be give by a linear combination of known thickness functions of the coordinate $X_3$ and unknown in-plane functions of the coordinates $(X_1,X_2)$. More specifically, the ith generalized displacement component $u_i^{\langle \ell \rangle }$ of ${\mathit{u}}^{\langle \ell \rangle }$ is written as

$$u_i^{\langle \ell \rangle }(X_1,X_2,X_3)=\mathit{\sum }_{k=0}^{N_i^{\langle \ell \rangle }}{f}_{ik}^{\langle \ell \rangle }\left(X_3\right)U_{ik}^{\langle \ell \rangle }(X_1,X_2),\forall \ell =1,\cdots ,N_{\ell },$$

(8)

where $N_i^{\langle \ell \rangle }$ denotes the order of expansion, $f_{ik}^{\langle \ell \rangle }\left(X_3\right)$ is kth known thickness function and $U_{ik}^{\langle \ell \rangle }(X_1,X_2)$ is the unknown in-plane function associated to ith generalized displacement component of the ℓth layer.

Consistently with the notation introduced in [30,31], Equation (8) can be expressed in matrix form as

$${\mathit{u}}^{\langle \ell \rangle }={\mathit{F}}^{\langle \ell \rangle }\left(X_3\right){\mathit{U}}^{\langle \ell \rangle }(X_1,X_2),\forall \ell =1,\cdots ,N_{\ell },$$

(9)

where ${\mathit{U}}^{\langle \ell \rangle }(X_1,X_2)$ is a $N_U^{\langle \ell \rangle }\times 1$ vector collecting the in-plane functions and ${\mathit{F}}^{\langle \ell \rangle }\left(X_3\right)$ is a $4\times N_U^{\langle \ell \rangle }$ matrix collecting the corresponding thickness functions being $N_U^{\langle \ell \rangle }\equiv (N_1^{\langle \ell \rangle }+N_2^{\langle \ell \rangle }+N_3^{\langle \ell \rangle }+N_4^{\langle \ell \rangle }+4)$. Further details about the explicit expression of the thickness functions and how the in-plane functions are ordered within the matrix ${\mathit{U}}^{\langle \ell \rangle }(X_1,X_2)$ to satisfy the continuity at the layers’ interfaces can be found in [31].

Eventually, substituting Equation (9) into Equations (6) and (7) allows expressing the generalized strain vector ${\mathit{ϵ}}^{\langle \ell \rangle }$ and the generalized stress vector ${\mathit{\tau }}^{\langle \ell \rangle }$ in terms of ${\mathit{U}}^{\langle \ell \rangle }$.

2.3. Governing Equations

The partial differential equations governing the electro-mechanical behavior of a multilayered piezoelectric plate, modelled according to the LW formulation derived in the previous section, can be obtained from the Principle of Virtual Displacements (PVD) for piezoelectric media (see, e.g., [35]). For the generic ℓth layer, the PVD reads

$${\int }_{V^{\langle \ell \rangle }}(\delta {\mathit{\gamma }}^{\langle \ell \rangle \top }{\mathit{\sigma }}^{\langle \ell \rangle }-\delta {\mathit{E}}^{\langle \ell \rangle \top }{\mathit{D}}^{\langle \ell \rangle })\mathrm{d}V={\int }_{V^{\langle \ell \rangle }}(\delta {\mathit{d}}^{\langle \ell \rangle \top }{\mathit{b}}_m^{\langle \ell \rangle }-\delta {\varphi }^{\langle \ell \rangle }{\rho }_e^{\langle \ell \rangle })\mathrm{d}V+{\int }_{{\partial V}^{\langle \ell \rangle }}(\delta {\mathit{d}}^{\langle \ell \rangle \top }{\mathit{t}}_m^{\langle \ell \rangle }+\delta {\varphi }^{\langle \ell \rangle }{\sigma }_e^{\langle \ell \rangle })\mathrm{d}\partial V,$$

(10)

where ${\mathit{b}}_m^{\langle \ell \rangle }$ is the vector of the mechanical volume forces, ${\rho }_e^{\langle \ell \rangle }$ is the body electric charge, ${\mathit{t}}_m^{\langle \ell \rangle }$ is the vector of the mechanical tractions and ${\sigma }_e^{\langle \ell \rangle }$ is the surface electric charge. Upon using the generalized notation introduced in Equation (5), the PVD can be written in a more compact form as

$${\int }_{V^{\langle \ell \rangle }}\delta {\mathit{ϵ}}^{\langle \ell \rangle \top }{\mathit{\tau }}^{\langle \ell \rangle }\mathrm{d}V={\int }_{V^{\langle \ell \rangle }}\delta {\mathit{u}}^{\langle \ell \rangle \top }{\mathit{b}}^{\langle \ell \rangle }\mathrm{d}V+{\int }_{{\partial V}^{\langle \ell \rangle }}\delta {\mathit{u}}^{\langle \ell \rangle \top }{\mathit{t}}^{\langle \ell \rangle }\mathrm{d}\partial V,$$

(11)

where ${\mathit{b}}^{\langle \ell \rangle }\equiv {\{{\mathit{b}}_m^{\langle \ell \rangle },-{\rho }_e^{\langle \ell \rangle }\}}^{\top }$ denotes the generalized volume loads and ${\mathit{t}}^{\langle \ell \rangle }\equiv {\{{\mathit{t}}_m^{\langle \ell \rangle },{\sigma }_e^{\langle \ell \rangle }\}}^{\top }$ denotes the generalized surface tractions.

Substituting Equations (6), (7) and (9) into Equation (11), one obtains

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}\left[\frac{\partial \delta {\mathit{U}}^{\langle \ell \rangle \top }}{{\partial X}_{\alpha }}\left({\mathit{Q}}_{\alpha \beta }^{\langle \ell \rangle }\frac{{\partial \mathit{U}}^{\langle \ell \rangle }}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}^{\langle \ell \rangle }{\mathit{U}}^{\langle \ell \rangle }\right)+\delta {\mathit{U}}^{\langle \ell \rangle \top }\left({\mathit{R}}_{\alpha 3}^{\langle \ell \rangle \top }\frac{{\partial \mathit{U}}^{\langle \ell \rangle }}{{\partial X}_{\alpha }}+{\mathit{S}}_{33}^{\langle \ell \rangle }{\mathit{U}}^{\langle \ell \rangle }\right)\right]\mathrm{d}\mathsf{\Omega }=\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\int }_{\mathsf{\Omega }}\delta {\mathit{U}}^{\langle \ell \rangle \top }{\mathit{B}}^{\langle \ell \rangle }\mathrm{d}\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\delta {\mathit{U}}^{\langle \ell \rangle \top }{\mathit{T}}_b^{\langle \ell \rangle }\mathrm{d}\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\delta {\mathit{U}}^{\langle \ell \rangle \top }{\mathit{T}}_t^{\langle \ell \rangle }\mathrm{d}\mathsf{\Omega }+{\int }_{\partial \mathsf{\Omega }}\delta {\mathit{U}}^{\langle \ell \rangle \top }{\mathit{T}}^{\langle \ell \rangle }\mathrm{d}\partial \mathsf{\Omega },\end{array}$$

(12)

which represents the PVD for the ℓth layer of a piezoelectric plate modelled via LW theories. In Equation (12), the matrices ${\mathit{Q}}_{\alpha \beta }^{\langle \ell \rangle }$, ${\mathit{R}}_{\alpha 3}^{\langle \ell \rangle }$ and ${\mathit{S}}_{33}^{\langle \ell \rangle }$ are obtained by the following integrals defined over the thickness of the layer

$${\mathit{Q}}_{\alpha \beta }^{\langle \ell \rangle }\equiv {\int }_{z_b^{\langle \ell \rangle }}^{z_t^{\langle \ell \rangle }}{\mathit{F}}^{\langle \ell \rangle \top }{\mathit{I}}_{\alpha }^{\top }{\mathit{p}}^{\langle \ell \rangle }{\mathit{I}}_{\beta }{\mathit{F}}^{\langle \ell \rangle }\mathrm{d}{X}_3,$$

(13a)

$${\mathit{R}}_{\alpha 3}^{\langle \ell \rangle }\equiv {\int }_{z_b^{\langle \ell \rangle }}^{z_t^{\langle \ell \rangle }}{\mathit{F}}^{\langle \ell \rangle \top }{\mathit{I}}_{\alpha }^{\top }{\mathit{p}}^{\langle \ell \rangle }{\mathit{I}}_3\frac{{\mathrm{d}\mathit{F}}^{\langle \ell \rangle }}{{\mathrm{d}X}_3}\mathrm{d}{X}_3,$$

(13b)

$${\mathit{S}}_{33}^{\langle \ell \rangle }\equiv {\int }_{z_b^{\langle \ell \rangle }}^{z_t^{\langle \ell \rangle }}\frac{{\mathrm{d}\mathit{F}}^{\langle \ell \rangle \top }}{{\mathrm{d}X}_3}{\mathit{I}}_3^{\top }{\mathit{p}}^{\langle \ell \rangle }{\mathit{I}}_3\frac{{\mathrm{d}\mathit{F}}^{\langle \ell \rangle }}{{\mathrm{d}X}_3}\mathrm{d}{X}_3.$$

(13c)

On the other hand, ${\mathit{B}}^{\langle \ell \rangle }$ stems from the generalized volume loads ${\mathit{b}}^{\langle \ell \rangle }$, ${\mathit{T}}_b^{\langle \ell \rangle }$ from the generalized surface tractions defined over the bottom surface of the layer and ${\mathit{T}}_t^{\langle \ell \rangle }$ from the generalized surface tractions defined over the top surface of the layer, i.e.,

$${\mathit{B}}^{\langle \ell \rangle }\equiv {\int }_{z_b^{\langle \ell \rangle }}^{z_t^{\langle \ell \rangle }}{\mathit{F}}^{\langle \ell \rangle \top }{\mathit{b}}^{\langle \ell \rangle }\mathrm{d}{X}_3,$$

(14a)

$${\mathit{T}}_b^{\langle \ell \rangle }\equiv {\mathit{F}}^{\langle \ell \rangle \top }(X_3=z_b^{\langle \ell \rangle }){\mathit{t}}^{\langle \ell \rangle }\mathrm{and}{\mathit{T}}_t^{\langle \ell \rangle }\equiv {\mathit{F}}^{\langle \ell \rangle \top }(X_3=z_t^{\langle \ell \rangle }){\mathit{t}}^{\langle \ell \rangle }.$$

(14b)

Similarly, the generalized surface tractions defined over the lateral surface of the ℓth layer contribute to the definition of the boundary term

$${\mathit{T}}^{\langle \ell \rangle }\equiv {\int }_{z_b^{\langle \ell \rangle }}^{z_t^{\langle \ell \rangle }}{\mathit{F}}^{\langle \ell \rangle \top }{\mathit{t}}^{\langle \ell \rangle }\mathrm{d}{X}_3.$$

(15)

Equation (12) represents the layer building block allowing to state the generalized PVD for the whole piezoelectric plate. In fact, upon summing over the total number of layers and considering the equilibrium condition ${\mathit{t}}_m^{\langle \ell \rangle }(X_1,X_2,z_t^{\langle \ell \rangle })+{\mathit{t}}_m^{\langle \ell +1\rangle }(X_1,X_2,z_b^{\langle \ell +1\rangle })=\mathbf{0}$ and the zero net surface charge condition ${\sigma }_e^{\langle \ell \rangle }(X_1,X_2,z_t^{\langle \ell \rangle })+{\sigma }_e^{\langle \ell +1\rangle }(X_1,X_2,z_b^{\langle \ell +1\rangle })=0$, for $\ell =1,\cdots ,N_{\ell }-1$, the PVD for the assembled plate is written as

$$\begin{array}{c}{\int }_{\mathsf{\Omega }}\left[\frac{\partial \delta {\mathit{U}}^{\top }}{{\partial X}_{\alpha }}\left({\mathit{Q}}_{\alpha \beta }\frac{\partial \mathit{U}}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}\mathit{U}\right)+\delta {\mathit{U}}^{\top }\left({\mathit{R}}_{\alpha 3}^{\top }\frac{\partial \mathit{U}}{{\partial X}_{\alpha }}+{\mathit{S}}_{33}\mathit{U}\right)\right]\mathrm{d}\mathsf{\Omega }=\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}={\int }_{\mathsf{\Omega }}\delta {\mathit{U}}^{\top }\mathit{B}\mathrm{d}\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\delta {\mathit{U}}^{\top }{\mathit{T}}_b\mathrm{d}\mathsf{\Omega }+{\int }_{\mathsf{\Omega }}\delta {\mathit{U}}^{\top }{\mathit{T}}_t\mathrm{d}\mathsf{\Omega }+{\int }_{\partial \mathsf{\Omega }}\delta {\mathit{U}}^{\top }\mathit{T}\mathrm{d}\partial \mathsf{\Omega },\end{array}$$

(16)

where the LW assembly technique described in [31,36] is employed to define the vectors $\mathit{U}$, $\mathit{B}$, ${\mathit{T}}_b$, ${\mathit{T}}_t$ and $\mathit{T}$ and the matrices ${\mathit{Q}}_{\alpha \beta }$, ${\mathit{R}}_{\alpha 3}$ and ${\mathit{S}}_{33}$. It is worth noting that $\mathit{U}$ collects the $N_U$ in-plane functions, which are the unknowns of the present LW formulation of piezoelectric plates.

Equation (16) represents the weak form of the governing equations for the considered problem; eventually, through standard variational calculus, it is possible to obtain the corresponding strong from of the governing equations and the relative boundary conditions as follows

$$\left\{\begin{array}{cc}-\frac{\partial }{{\partial X}_{\alpha }}\left({\mathit{Q}}_{\alpha \beta }\frac{\partial \mathit{U}}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}\mathit{U}\right)+{\mathit{R}}_{\alpha 3}^{\top }\frac{\partial \mathit{U}}{{\partial X}_{\alpha }}+{\mathit{S}}_{33}\mathit{U}=\overline{\mathit{B}},\hfill & \mathrm{in}\mathsf{\Omega }\hfill \\ \mathit{U}=\overline{\mathit{U}},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_D\hfill \\ n_{\alpha }\left({\mathit{Q}}_{\alpha \beta }\frac{\partial \mathit{U}}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}\mathit{U}\right)=\overline{\mathit{T}},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_N\hfill \end{array}\right.,$$

(17)

where $n_{\alpha }$ is the $\alpha$th component of the outer unit normal of the boundary $\partial \mathsf{\Omega }$, $\overline{\mathit{B}}\equiv \mathit{B}+{\mathit{T}}_b+{\mathit{T}}_t$ are the plate’s volume loads, $\overline{\mathit{U}}$ are the plate’s prescribed boundary displacements, $\overline{\mathit{T}}$ are the plate’s prescribed boundary tractions and $\partial {\mathsf{\Omega }}_D$ and $\partial {\mathsf{\Omega }}_N$ denote non-overlapping subsets of $\partial \mathsf{\Omega }$ over which Dirichlet and Neumann boundary conditions apply.

3. Discontinuous Galerkin Formulation

The equations given in (17) denote a system of second-order elliptic partial differential equations with suitably defined Dirichlet and Neumann boundary conditions. The derivation of a solution scheme based on the discontinuous Galerkin method generally involves [23]: (i) the introduction of an auxiliary variable that allows to rewrite the second-order system as a first-order system; and (ii) a weak statement of said first-order system for each element of the mesh discretizing the considered domain. Here, the derivation proposed by Gulizzi et al. [30,31] is employed.

The form of Equation (17) suggests defining as auxiliary variable the quantity

$${\mathit{\sum }}_{\alpha }\equiv {\mathit{Q}}_{\alpha \beta }\frac{\partial \mathit{U}}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}\mathit{U},$$

(18)

so that the system of equations given in Equation (17) may be rewritten as

$$\left\{\begin{array}{cc}-\frac{{\partial \mathit{\sum }}_{\alpha }}{{\partial X}_{\alpha }}+{\mathit{R}}_{\alpha 3}^{\top }\frac{\partial \mathit{U}}{{\partial X}_{\alpha }}+{\mathit{S}}_{33}\mathit{U}=\overline{\mathit{B}},\hfill & \mathrm{in}\mathsf{\Omega }\hfill \\ {\mathit{\sum }}_{\alpha }={\mathit{Q}}_{\alpha \beta }\frac{\partial \mathit{U}}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}\mathit{U},\hfill & \mathrm{in}\mathsf{\Omega }\hfill \\ \mathit{U}=\overline{\mathit{U}},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_D\hfill \\ {\mathit{\sum }}_{\alpha }{n}_{\alpha }=\overline{\mathit{T}},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_N\hfill \end{array}\right..$$

(19)

Upon introducing the generic mesh element ${\mathsf{\Omega }}^{\left(e\right)}$ and the test functions $\mathit{V},{\mathit{\Gamma }}_{\alpha }\in {\mathcal{V}}_{hp}^{N_U}$, where ${\mathcal{V}}_{hp}^{N_U}$ is the space of discontinuous vector fields for the considered mesh, the weak statements of the first two equations in (19) can be written as

$${\int }_{{\mathsf{\Omega }}^{\left(e\right)}}\left[\frac{{\partial \mathit{V}}^{\top }}{{\partial X}_{\alpha }}{\mathit{\sum }}_{h\alpha }+{\mathit{V}}^{\top }\left({\mathit{R}}_{\alpha 3}^{\top }\frac{{\partial \mathit{U}}_h}{{\partial X}_{\alpha }}+{\mathit{S}}_{33}{\mathit{U}}_h\right)\right]\mathrm{d}\mathsf{\Omega }={\int }_{{\partial \mathsf{\Omega }}^{\left(e\right)}}{\mathit{V}}^{\top }{\hat{\mathit{\sum }}}_{\alpha }{n}_{\alpha }\mathrm{d}\partial \mathsf{\Omega }+{\int }_{{\mathsf{\Omega }}^{\left(e\right)}}{\mathit{V}}^{\top }\overline{\mathit{B}}\mathrm{d}\mathsf{\Omega },$$

(20a)

and

$${\int }_{{\mathsf{\Omega }}^{\left(e\right)}}{\mathit{\Gamma }}_{\alpha }^{\top }{\mathit{\sum }}_{h\alpha }\mathrm{d}\mathsf{\Omega }={\int }_{{\mathsf{\Omega }}^{\left(e\right)}}{\mathit{\Gamma }}_{\alpha }^{\top }\left({\mathit{Q}}_{\alpha \beta }\frac{{\partial \mathit{U}}_h}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}{\mathit{U}}_h\right)\mathrm{d}\mathsf{\Omega }+{\int }_{{\partial \mathsf{\Omega }}^{\left(e\right)}}({\mathit{\Gamma }}_{\alpha }^{\top }{\mathit{Q}}_{\alpha \beta }+{\mathit{V}}^{\top }{\mathit{R}}_{\beta 3}^{\top })(\hat{\mathit{U}}-{\mathit{U}}_h)n_{\beta }\mathrm{d}\partial \mathsf{\Omega }.$$

(20b)

In Equation (20a,b), $\partial {\mathsf{\Omega }}^{\left(e\right)}$ denotes the boundary of ${\mathsf{\Omega }}^{\left(e\right)}$; ${\mathit{U}}_h$ and ${\mathit{\sum }}_{h\alpha }$ denote the dG solutions and are in general approximations of the exact solutions $\mathit{U}$ and ${\mathit{\sum }}_{\alpha }$, respectively; and $\hat{\mathit{U}}$ and ${{\displaystyle \hat{\mathit{\sum }}}}_{\alpha }$ are the so-called numerical fluxes, which allow recovering the continuity of the solution at the elements’ boundaries.

Different choices of the numerical fluxes lead to different discontinuous Galerkin formulations [23]. However, prior to giving their explicit expression, it is useful to introduce $\partial {\mathsf{\Omega }}_D^{\left(e\right)}$ and $\partial {\mathsf{\Omega }}_N^{\left(e\right)}$ as the subsets of $\partial {\mathsf{\Omega }}^{\left(e\right)}$ where Dirichlet and Neumann boundary conditions, respectively, are to be enforced; furthermore, considering that the mesh discretization introduces a subdivision of the plate domain $\mathsf{\Omega }$, we define $\partial {\mathsf{\Omega }}_I^{\left(e\right)}$ as the part of $\partial {\mathsf{\Omega }}^{\left(e\right)}$ that the element e shares with its neighboring elements. Then, upon employing the interior penalty dG formulation introduced in [30,31], the numerical fluxes are chosen as follows

$$\tilde{\mathit{U}}=\left\{\begin{array}{cc}\left\{{\mathit{U}}_h\right\},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_I^{\left(e\right)}\hfill \\ \overline{\mathit{U}},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_D^{\left(e\right)}\hfill \\ {\mathit{U}}_h,\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_N^{\left(e\right)}\hfill \end{array}\right.$$

(21)

and

$$\left\{\begin{array}{cc}{{\displaystyle \hat{\mathit{\sum }}}}_{\alpha }=\{{\mathit{Q}}_{\alpha \beta }\frac{{\partial \mathit{U}}_h}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}{\mathit{U}}_h\}-\mu {⟦{\mathit{U}}_h⟧}_{\alpha },\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_I^{\left(e\right)}\hfill \\ {{\displaystyle \hat{\mathit{\sum }}}}_{\alpha }={\mathit{Q}}_{\alpha \beta }\frac{{\partial \mathit{U}}_h}{{\partial X}_{\beta }}+{\mathit{R}}_{\alpha 3}{\mathit{U}}_h-\mu ({\mathit{U}}_h-\overline{\mathit{U}})n_{\alpha },\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_D^{\left(e\right)}\hfill \\ {{\displaystyle \hat{\mathit{\sum }}}}_{\alpha }{n}_{\alpha }=\overline{\mathit{T}},\hfill & \mathrm{on}\partial {\mathsf{\Omega }}_N^{\left(e\right)}\hfill \end{array}\right.,$$

(22)

where $\{•\}\equiv \frac{1}{2}\left({•}^{\left(e\right)}+{•}^{\left(e^{\prime }\right)}\right)$ and ${⟦•⟧}_{\alpha }\equiv n_{\alpha }^{\left(e\right)}{•}^{\left(e\right)}+n_{\alpha }^{\left(e^{\prime }\right)}{•}^{\left(e^{\prime }\right)}$ are, respectively, the so-called average and jump operators, which are defined at the interface between the element e and its generic neighbour $e^{\prime }$. It is noted that the numerical flux ${{\displaystyle \hat{\mathit{\sum }}}}_{\alpha }$ given in (22) depends on the so-called penalty parameter $\mu$, whose influence on the numerical solution is addressed in Section 4.

Eventually, upon substituting Equations (21) and (22) into Equation (20a,b), respectively, setting ${\mathit{\Gamma }}_{\alpha }\equiv \partial \mathit{V}/\partial X_{\alpha }$ and summing over all the mesh elements, Equation (20a,b) can be combined to obtain the interior penalty dG variational statement for multilayered piezoelectric plates:

$$B(\mathit{V},{\mathit{U}}_h)=F(\mathit{V},\overline{\mathit{B}},\overline{\mathit{T}},\overline{\mathit{U}}),$$

(23)

where

$$\begin{array}{c}B(\mathit{V},{\mathit{U}}_h)={\int }_{{\mathsf{\Omega }}_h}\left[\frac{\partial {\mathit{V}}^{\top }}{\partial X_{\alpha }}\left({\mathit{Q}}_{\alpha \beta }\frac{\partial {\mathit{U}}_h}{\partial X_{\beta }}+{\mathit{R}}_{\alpha 3}{\mathit{U}}_h\right)+{\mathit{V}}^{\top }\left({\mathit{R}}_{\alpha 3}^{\top }\frac{\partial {\mathit{U}}_h}{\partial X_{\alpha }}+{\mathit{S}}_{33}{\mathit{U}}_h\right)\right]\mathrm{d}\mathsf{\Omega }+\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\int }_{{\partial \mathsf{\Omega }}_{hI}}\left[{⟦\mathit{V}⟧}_{\alpha }^{\top }\left\{{\mathit{Q}}_{\alpha \beta }\frac{\partial {\mathit{U}}_h}{\partial X_{\beta }}+{\mathit{R}}_{\alpha 3}{\mathit{U}}_h\right\}+\left\{\frac{\partial {\mathit{V}}^{\top }}{\partial X_{\alpha }}{\mathit{Q}}_{\alpha \beta }+{\mathit{V}}^{\top }{\mathit{R}}_{\beta 3}^{\top }\right\}{⟦{\mathit{U}}_h⟧}_{\beta }\right]\mathrm{d}\partial \mathsf{\Omega }+\\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-{\int }_{{\partial \mathsf{\Omega }}_{hD}}\left[n_{\alpha }{\mathit{V}}^{\top }\left({\mathit{Q}}_{\alpha \beta }\frac{\partial {\mathit{U}}_h}{\partial X_{\beta }}+{\mathit{R}}_{\alpha 3}{\mathit{U}}_h\right)+\left(\frac{{\partial \mathit{V}}^{\top }}{\partial X_{\alpha }}{\mathit{Q}}_{\alpha \beta }+{\mathit{V}}^{\top }{\mathit{R}}_{\beta 3}^{\top }\right){\mathit{U}}_h{n}_{\beta }\right]\mathrm{d}\partial \mathsf{\Omega }+\\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\int }_{{\partial \mathsf{\Omega }}_{hI}}\mu {⟦\mathit{V}⟧}_{\alpha }^{\top }{⟦{\mathit{U}}_h⟧}_{\alpha }\mathrm{d}\partial \mathsf{\Omega }+{\int }_{{\partial \mathsf{\Omega }}_{hD}}\mu {\mathit{V}}^{\top }{\mathit{U}}_h\mathrm{d}\partial \mathsf{\Omega }\end{array}$$

(24)

and

$$\begin{array}{c}F(\mathit{V},\overline{\mathit{B}},\overline{\mathit{T}},\overline{\mathit{U}})={\int }_{{\mathsf{\Omega }}_h}{\mathit{V}}^{\top }\overline{\mathit{B}}\mathrm{d}\mathsf{\Omega }+{\int }_{{\partial \mathsf{\Omega }}_{hN}}{\mathit{V}}^{\top }\overline{\mathit{T}}\mathrm{d}\partial \mathsf{\Omega }+\hfill \\ \hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} -{\int }_{{\partial \mathsf{\Omega }}_{hD}}\left(\frac{{\partial \mathit{V}}^{\top }}{{\partial X}_{\alpha }}{\mathit{Q}}_{\alpha \beta }+{\mathit{V}}^{\top }{\mathit{R}}_{\beta 3}^{\top }\right)\overline{\mathit{U}}{n}_{\beta }\mathrm{d}\partial \mathsf{\Omega }+{\int }_{{\partial \mathsf{\Omega }}_{hD}}\mu {\mathit{V}}^{\top }\overline{\mathit{U}}\mathrm{d}\partial \mathsf{\Omega }.\end{array}$$

(25)

The integrals appearing in Equations (24) and (25) are referred to as the broken integrals and are defined as

$${\int }_{{\mathsf{\Omega }}_h}•\mathrm{d}\mathsf{\Omega }\equiv {\mathit{\sum }}_{e=1}^{N_e}{\int }_{{\mathsf{\Omega }}^{\left(e\right)}}{•}^{\left(e\right)}\mathrm{d}\mathsf{\Omega },{\int }_{{\partial \mathsf{\Omega }}_{hD}}•\mathrm{d}\partial \mathsf{\Omega }\equiv {\mathit{\sum }}_{e=1}^{N_e}{\int }_{{\partial \mathsf{\Omega }}_D^{\left(e\right)}}{•}^{\left(e\right)}\mathrm{d}\partial \mathsf{\Omega },{\int }_{{\partial \mathsf{\Omega }}_{hN}}•\mathrm{d}\partial \mathsf{\Omega }\equiv {\mathit{\sum }}_{e=1}^{N_e}{\int }_{{\partial \mathsf{\Omega }}_N^{\left(e\right)}}{•}^{\left(e\right)}\mathrm{d}\partial \mathsf{\Omega }$$

(26)

and

$${\int }_{{\partial \mathsf{\Omega }}_{hI}}•\mathrm{d}\partial \mathsf{\Omega }\equiv \mathit{\sum }_{i=1}^{N_i}{\int }_{I^{\left(i\right)}}{•}^{\left(i\right)}\mathrm{d}\partial \mathsf{\Omega },$$

(27)

where $N_e$ is the total number of mesh elements, $N_i$ is the total number of inter-element interfaces that have been introduced along with the mesh discretization and $I^i$ is the generic ith interface between the neighboring elements e and $e^{\prime }$, i.e., $I^{\left(i\right)}={\partial \mathsf{\Omega }}^{\left(e\right)}\cap {\partial \mathsf{\Omega }}^{\left(e^{\prime }\right)}$, ${\mathsf{\Omega }}_h\equiv {\cup }_{e=1}^{N_e}{\mathsf{\Omega }}^{\left(e\right)}$, ${\partial \mathsf{\Omega }}_{hD}\equiv {\cup }_{e=1}^{N_e}\partial {\mathsf{\Omega }}_D^{\left(e\right)}$, ${\partial \mathsf{\Omega }}_{hN}\equiv {\cup }_{e=1}^{N_e}\partial {\mathsf{\Omega }}_N^{\left(e\right)}$ and $\partial {\mathsf{\Omega }}_{hI}\equiv {\cup }_{i=1}^{N_i}{I}^{\left(i\right)}$.

It is underlined that the bilinear form given in Equation (24), stemming from the choice of the numerical fluxes given in Equations (21) and (22), is symmetric and verifies the Galerkin orthogonality. As shown in Section 4, this allows avoiding super-penalization [23] and obtaining optimal $hp$-convergence with values of the penalty parameter $\mu$ that have the same order of magnitude of the constants of the differential operators involved in the governing equations.

4. Numerical Results

In this section, the present formulation is employed to analyze the following test cases: (i) a simply-supported square piezoelectric layer subjected to a bi-sinusoidal load and zero surface charge at the top and bottom surfaces; and (ii) a simply-supported square smart plate subjected to a bi-sinusoidal load and with grounded top and bottom surfaces. The materials and the stacking sequences of the considered numerical tests are reported in

Table 1. Properties of the considered piezoelectric materials, as taken from Ballhause et al. [36].

Material ID	Property	Component	Value	Unit
PZT-4	Young’s moduli	$E_1$, $E_2$	81.3	GPa
		$E_3$	64.5	GPa
	Poisson’s ratios	${\nu }_{23}$, ${\nu }_{13}$	0.432	-
		${\nu }_{12}$	0.329	-
	Shear moduli	$G_{23}$, $G_{13}$	25.6	GPa
		$G_{12}$	30.6	GPa
	Dielectric constants	${\kappa }_1$, ${\kappa }_2$	$6.46$×${10}^{-9}$	C/(V·m)
		${\kappa }_3$	$5.62$×${10}^{-9}$	C/(V·m)
	Piezoelectric constants	$e_{15}$, $e_{24}$	$12.72$	C/m${}^2$
		$e_{31}$, $e_{32}$	$-5.2$	C/m${}^2$
		$e_{33}$	$15.08$	C/m${}^2$
Gr/Ep	Young’s moduli	$E_1$	132.38	GPa
		$E_2$, $E_3$	10.756	GPa
	Poisson’s ratios	${\nu }_{23}$	0.49	-
		${\nu }_{12}$, ${\nu }_{13}$	0.24	-
	Shear moduli	$G_{23}$	3.606	GPa
		$G_{12}$, $G_{13}$	5.6537	GPa
	Dielectric constants	${\kappa }_1$	$1.53$×${10}^{-11}$	C/(V·m)
		${\kappa }_2$, ${\kappa }_3$	$1.31$×${10}^{-11}$	C/(V·m)

and

Table 2. Properties of the considered plate stacking sequence.

Plate ID	Material	Layup	Layers Thickness [m]
${\mathrm{P}}_1$	[PZT-4]	$\left[0\right]$	[0.025]
${\mathrm{P}}_2$	[PZT-4,Gr/Ep,Gr/Ep,PZT-4]	$[0,0,90,0]$	[0.025,0.1,0.1,0.025]

, respectively. The same order of expansion is used for each layer and for each component of the generalized displacement field ${\mathit{u}}^{\langle \ell \rangle }$, i.e., $N_1^{\langle \ell \rangle }=N_2^{\langle \ell \rangle }=N_3^{\langle \ell \rangle }=N_4^{\langle \ell \rangle }=N$, $\forall \ell =1,\cdots ,N_{\ell }$, and the corresponding theory is denoted by ${\mathrm{LW}}_N$. Legendre polynomials are used as basis functions.

The numerical results are reported in terms of non-dimensional quantities, denoted by the hat $\widehat{·}$ and defined as

$$\widehat{\mathit{d}}\equiv \frac{\mathit{d}}{l},\widehat{\varphi }\equiv \frac{\varphi }{l\sqrt{c/\kappa }},\widehat{\sigma }\equiv \frac{\sigma }{c},\mathrm{and}\widehat{\mathit{D}}\equiv \frac{\mathit{D}}{\sqrt{\kappa c}},$$

(28)

where l, c and $\kappa$ denote, respectively, a characteristic length and the stiffness and dielectric permittivity constants of the considered problem.

4.1. Piezoelectric Layer

In the first test case, the plate labeled as ${\mathrm{P}}_1$ in Table 2 is considered. The plate is shown in Figure 2a and consists of a single piezoelectric layer that occupies the volume $V=[0,L]\times [0,L]\times [-z/2,z/2]$, where $L=1$ m and $z=0.025$ m. For this problem, the non-dimensionalization constants are chosen as $l=1$ m, $c=81.3$ GPa and $\kappa =6.46$ × ${10}^{-9}$ C/(V·m). The plate is subjected to simply-supported boundary conditions over the lateral surfaces, zero tractions and zero surface charges over the bottom surface and a bisinusoidal vertical traction and zero surface charge over the top surface; more specifically, the mechanical tractions over the top surface have the following expression

$${\overline{\mathit{t}}}_m=qsin\left(\frac{\pi X_1}{L}\right)sin\left(\frac{\pi X_2}{L}\right){\hat{\mathbf{e}}}_3,\mathrm{at}{X}_3=z/2,$$

(29)

where q denotes the amplitude and ${\hat{\mathbf{e}}}_3$ is the unit vector along the $X_3$ direction.

Given the chosen material, geometry and set of boundary conditions, this problem admits Navier-type solutions (see, e.g., [36]), which are used to assess the accuracy of the proposed formulation as a function of the basis functions order p, the penalty parameter $\mu$ and the mesh size $h=L/n$, being n the number of mesh elements per side of the plate.

Eventually, the following error measure is introduced:

$$e\left({\mathit{U}}_h\right)\equiv \frac{\left|\right|{\mathit{U}}_h-{\mathit{U}}_{\mathrm{ref}}{\left|\right|}_2}{\left|\right|{\mathit{U}}_{\mathrm{ref}}{\left|\right|}_2},$$

(30)

where ${\left|\right|\mathit{U}\left|\right|}_2$ is the $L^2\left(\mathsf{\Omega }\right)$ norm, ${\mathit{U}}_{\mathrm{ref}}$ is the Navier solution and ${\mathit{U}}_h$ is the solution computed using the present formulation. It is worth stressing that the solution fields ${\mathit{U}}_h={\mathit{U}}_h(X_1,X_2)$ provided by the developed dG scheme and assessed through the error definition in Equation (30) are those appearing in the generalized kinematical model in Equation (9), which identifies them as a generalization of the mid-plane displacements and rotations employed in classical plate theories.

Figure 3a shows the effect of the penalty parameter $\mu$ introduced in Equation (24) on the error measure $e\left({\mathit{U}}_h\right)$ as a function of the basis functions order p and for a mesh grid $2\times 2$ and a layer-wise theory LW${}_3$. In the figure, $p_{11}\equiv {\left(\mathit{p}\right)}_{11}$, where $\mathit{p}$ is the generalized constitutive matrix of the considered piezoelectric material as introduced in Equation (7). Consistently with the results obtained in [30,31], the magnitude of the penalty parameter has little influence on the error values as long as it is chosen sufficiently large. It is also worth noting that, unlike super-penalized techniques, $\mu$ can be chosen one or two orders of magnitude larger than $p_{11}$, thus not affecting the conditioning of the system of equations; this is a typical feature of the interior penalty dG formulation [23].

Then, upon choosing $\mu /p_{11}=10/h$ as common in interior penalty dG formulations, the error $e\left({\mathit{U}}_h\right)$ is computed as a function of the mesh size h and the polynomial order p; the obtained results are presented in Figure 3b and show the optimal convergence of the present formulation.

For the sake of completeness, it is interesting to show a comparison between the LW theories and the 3D solution, which can be obtained in closed form [6]. Figure 4 shows the through-the-thickness behavior of some selected mechanical and electrical quantities, which have been evaluated using the exact 3D solution and the present formulation with two different layer-wise theories, namely LW${}_2$, LW${}_3$, LW${}_4$ and LW${}_5$, a polynomial order $p=7$ and a mesh grid $2\times 2$. As shown in the figure, a low-order LW theory is able to capture accurately the values of the primary variables, i.e., the displacement components and the electric potential; however, high-order LW theories are required to reproduce the through-the-thickness stress and electric displacement distribution, especially when the out-of-plane components are considered.

4.2. Smart Plate

In the second test case, the dG formulation is employed to evaluate the electro-mechanical response of the smart plate analyzed in [36], which corresponds to the plate labeled as ${\mathrm{P}}_2$ in Table 2 and shown in Figure 2b. The plate occupies the volume $V=[0,L]\times [0,L]\times [-z/2,z/2]$, where $L=1$ m and $z=0.25$ m, and consists of four layers: two outer piezoelectric layers with thickness $z^{\langle 1\rangle }=z^{\langle 4\rangle }=0.025$ m made of PZT-4 and two structural layers with thickness $z^{\langle 2\rangle }=z^{\langle 3\rangle }=0.1$ m made of Gr/Ep in the [0/90] configuration. For this problem, the non-dimensionalization constants are chosen as $l=1$ m, $c=132.38$ GPa and $\kappa =6.46$×${10}^{-9}$ C/(V·m).

The plate is subjected to the same mechanical boundary conditions of the piezoelectric layer considered in the previous test case, i.e., a bisinusoidal normal traction over the top surface. However, from an electrical standpoint, two sets of boundary conditions are considered, namely closed-circuit and open-circuit, which correspond to zero potential and zero surface charge, respectively, at the bottom and top surfaces. These sets of boundary conditions aim at reproducing realistic situations when smart plates are employed as actuator or sensor devices.

In addition, in this test case, the 3D solution can be computed in closed form [6]. Figure 5 shows the through-the-thickness behavior of some selected mechanical and electrical quantities in the short-circuit configuration, whereas Figure 6 shows the same quantities in the open-circuit configuration. The curves have been evaluated using the exact 3D solution and the present formulation using the layer-wise theory LW${}_7$, a polynomial order $p=7$ and a mesh grid $2\times 2$. As shown in the figures, a perfect agreement is observed between the exact 3D results and the results obtained using the present formulation.

5. Conclusions

In this work, we present a novel high-order formulation for multilayered piezoelectric plates. The method is based on the combination of a variable-order interior penalty discontinuous Galerkin method and general layer-wise plate theories. One of the main advantages of the developed formulation is that it allows tuning the order of the numerical approximation throughout both the thickness and the modelling plane of the considered multilayered piezoelectric plates, thus enabling high-order layer-wise resolution of the electromechanical fields, which is typically obtainable via fully three-dimensional models. Additionally, the use of high-order basis functions allows to attain a certain level of accuracy with a smaller number of mesh elements, and therefore of total degrees of freedom, with respect to those needed by other popular low-order methods. The readers interested in these aspects are referred to [30], which presents a quantitative comparison between the present dG approach and standard FEM in terms of accuracy versus total degrees of freedom for the case of purely elastic plates. Numerical tests were performed on a single piezoelectric layer and on a smart multilayered plate consisting of a structural two-layer core and two outer piezoelectric layers, and the effect of the penalty parameter $\mu$ was investigated. The obtained results show that the present formulation offers optimal $hp$-convergence rate and that it is able to resolve with high-order accuracy both the in-plane and the through-the-thickness behavior of the primary variables as well as the behavior of the derived variables, i.e., mechanical stresses and electric displacements.