Problems of Hierarchical Modelling and hp-Adaptive Finite Element Analysis in Elasticity, Dielectricity and Piezoelectricity

Abstract

In this chapter, we consider theoretical and implementation difficulties in application of the hierarchical modelling and hp-adaptive finite element approach to elasticity, dielectricity and piezoelectricity. The main feature of the applied methodology is its generalizing character which is reflected by application of the same or analogous algorithms to three mentioned physical problems, including multi-physics problem of piezoelectricity, simple and complex physical description as well as simple and complex geometries. In contrast to the most common approaches dealing with a single physical phenomenon, described by a single physical model, within a single geometrical part, this chapter presents the ideas which brake and overcome such a simplicity. This presented chapter generalizes author’s hitherto accomplishments, in hierarchical models and hp-approximations of linear elasticity, onto dielectricity and piezoelectricity. The same refers to error estimation and adaptivity control. In this context, the main similarities and differences of three physical problems are of interest in this work.

1. Introduction

In this chapter of the book we extend our hitherto propositions concerning 3D-based hierarchical models of liner elasticity onto 3D-based linear dielectric and piezoelectric media. In the case of hierarchical models of linear elasticity we apply 3D-elasticity model, hierarchical shell models, first-order shell model and the solid-to-shell or shell-to-shell transition models. In the case of dielectricity we utilize 3D-dielectricity model, and the 3D-based hierarchical dielectric models as well. The piezoelectric case needs combination of two mentioned mechanical and electric hierarchies, so as to generate the hierarchy of 3D-based piezoelectric models. Any combination of the mentioned elastic and dielectric models is possible. As far as the hp-discretization is concerned we extend the ideas of hierarchical approximations, constrained approximations and the transition approximations of the displacement field onto electric potential field of dielectricity or the coupled electro-mechanical field of piezoelectricity. The mentioned approximations allow p-adaptivity (three-dimensional or longitudinal), q-adaptivity (transverse one), h-adaptivity (three-dimensional or two-dimensional ones) and M-adaptivity (model adaptivity). The error assessment in three classes of problems is based on the equilibrated residual methods (ERM) applied to total and approximation error estimations. The modelling error is obtained as a difference of the previous two errors. The estimated error values are utilized for adaptivity control. The adaptive procedures for dielectricity and piezoelectricity are obtained through the generalization of the three- or four-step strategies applied so far to the elasticity case. The difficulties in generalization of the above-mentioned methods of hierarchical modelling, hp-approximations, error estimation and adaptivity control onto electrical and electro-mechanical problems are addressed in this chapter.

2. Considered problems

In the chapter we consider five problems. The first three correspond to stationary problems of mechanical, electric and electro-mechanical equilibrium of the elastic, dielectric and piezoelectric media, respectively. The last two problems deal with free vibration problems of linear elasticity and linear piezoelectricity.

2.1. Elastostatic problem

Let us start with the standard formulation of the linear elasticity [1]. The problem local equations include equilibrium, constitutive and geometric relations:

$\begin{array}{l}{\sigma }_{ij,j}+f_i=\mathrm{0,}\\ {\sigma }_{ij}=D_{ijkl}{\epsilon }_{kl},\\ {\epsilon }_{ij}=1/2(u_{i,j}+u_{j,i}),x\in V\end{array}$E1

where $D_{ijkl},{\sigma }_{ij},{\epsilon }_{kl},$ i,j,k,l = 1,2,3 are the elasticity constants tensor, and the stress and strain tensors, respectively. The given vector of mass load is denoted as f_i , while u_i is the unknown vector of displacements. The above equations hold in volume V of the body.

The standard boundary conditions for stresses and displacements are:

$\begin{array}{l}{\sigma }_{ij}{n}_j=p_i,x\in S_P,\\ u_i=w_i,x\in S_W\end{array}$E2

where n_j denotes components of the normal to the body surface S, composed of its loaded S_P and supported S_W parts: $S=S_P\cup S_W$. The vectors p_i and w_i represent the given stresses and displacements on S_P and S_W, respectively.

The equivalent variational formulation results from minimization of the potential energy of the elastic body:

${\displaystyle {\int }_V}{D}_{ijkl}{v}_{i,j}{u}_{k,l}dV={\displaystyle {\int }_V}{v}_i{f}_{i}dV+{\displaystyle {\int }_{S_P}}{v}_i{p}_{i}dS$E3

where v_i represent admissible displacements conforming to the displacement boundary conditions.

The above variational functional can be utilized in the derivation of the global finite element equations of the form:

$K_M{q}^{q,hp}=F_V+F_S$E4

where $K_M,F_V,F_S$, are the stiffness matrix within the mechanical equilibrium problem, and the vectors of the nodal forces due to volume and surface loadings. The term $q^{q,hp}$ stands for the nodal displacement degrees of freedom (dofs). The applied hierarchical q,hp-approximation of displacements will be addressed in the next sections.

2.2. Electrostatic problem of dielectrics

The standard local formulation of linear dielectricity [2] consists of the Gauss law, here corresponding to the lack of volume charges, the constitutive relation and the electric field $E_j$ definition:

$\begin{array}{l}{d}_{i,i}=\mathrm{0,}\\ d_i={\gamma }_{ij}{E}_j,\\ E_j=-{\phi }_j,x\in V\end{array}$E5

Above, ${\gamma }_{ij}$, i,j = 1,2,3 stands for the dielectric constants tensor, while $d_i$ denotes the electric displacement vector. The scalar term $\varphi$ represents the electric potential field, searched in the volume $V$ of the dielectric.

The boundary conditions for the electric displacements and electric potential read:

$\begin{array}{l}{d}_i{n}_i=-c,x\in S_Q,\\ \varphi =\chi ,x\in S_F\end{array}$E6

where $c$ and $\chi$ are the given scalar values of the surface charge and electric potential on the parts $S_Q$ and $S_F$, respectively, of the surface $S$ of the dielectric body $(S=S_Q\cup S_F)$.

The corresponding variational formulation which reflects minimization of the potential electric energy can be described as follows:

${\displaystyle {\int }_V}{\gamma }_{ij}{\psi }_i{\varphi }_{j}dV={\displaystyle {\int }_{S_Q}}\psi cdS$E7

with $\psi$ being the admissible electric potential conforming to the second boundary condition of Eq. (6).

The finite element formulation can be expressed as follows:

$K_E{\phi }^{\rho ,h\pi }=F_Q$E8

where $K_E$ is the characteristic matrix of dielectricity, $F_Q$ denotes the characteristic nodal vector of electric charges and ${\phi }^{\rho ,h\pi }$ stands for the unknown nodal vector of electric potentials. The hierarchical ρ,hπ-approximation of the potential will be explained later on in this chapter.

2.3. Stationary electro-mechanical problem

Formally, the local formulation of the piezoelectric problem of electro-mechanical equilibrium [3] can be treated as a combination of the linear elasticity and linear dielectricity Eqs. (1) and (5):

$\begin{array}{l}{\sigma }_{ij,j}+f_i=\mathrm{0,}\\ d_{i,i}=\mathrm{0,}\\ {\sigma }_{ij}=D_{ijkl}{\epsilon }_{kl}-C_{kij}{E}_k,\\ d_i=C_{ikl}{\epsilon }_{kl}+{\gamma }_{ij}{E}_j,\\ {\epsilon }_{ij}=1/2(u_{i,j}+u_{j,i}),\\ E_j=-{\varphi }_j,x\in V\end{array}$E9

What couples both sets of equations are the modified constitutive relations, where the coupling piezoelectric constants tensor $C_{kij}$ appears.

The boundary conditions of the coupled problem are:

$\begin{array}{ll}{\sigma }_{ij}{n}_j=p_i,\hfill & x\in S_P,\hfill \\ d_i{n}_i=-c,\hfill & x\in S_Q,\hfill \\ u_i=w_i,\hfill & x\in S_W,\hfill \\ \phi =\chi ,\hfill & x\in S_F\hfill \end{array}$E10

The variational functional of the electro-mechanical problem consists of the terms of functionals (3) and (7) completed with the terms describing the piezoelectric coupling through the tensor $C_{ijk}$, i.e.

$\begin{array}{l}-{\displaystyle {\int }_V}{D}_{ijkl}{v}_{i,j}{u}_{k,l}dV+{\displaystyle {\int }_V}{C}_{ijk}{v}_{i,j}{\varphi }_{k}dV+{\displaystyle {\int }_V}{v}_i{f}_{i}dV+{\displaystyle {\int }_{S_P}}{v}_i{p}_{i}dS\\ +{\displaystyle {\int }_V}{\gamma }_{ik}{\psi }_i{\varphi }_{k}dV+{\displaystyle {\int }_V}{C}_{ikl}{\psi }_i{u}_{k,l}dV+{\displaystyle {\int }_{S_Q}}\psi cdS=0\end{array}$E11

The above variational formulation leads to the following finite element equations

$\begin{array}{c}{K}_M{q}^{q,hp}-K_C{\phi }^{\rho ,h\pi }=F_V+F_S,\\ K_C^T{q}^{q,hp}+K_E{\phi }^{\rho ,h\pi }=F_Q\end{array}$E12

where $K_C$ is the characteristic matrix of piezoelectric coupling.

The above finite element equations correspond to a very general case when both the direct and inverse piezoelectric phenomena are present. Substitution of the second Eq. (12) into the first one leads to the single combined equation from which the nodal displacements $q^{q,hp}$ can be calculated. The opposite substitution gives the combined equation from which the nodal electric potentials ${\phi }^{\rho ,h\pi }$ can be extracted. The first situation corresponds to the so-called actuation action of the piezoelectric, while the second one to the sensing action of the piezoelectric body. These two modes of action can be associated with the direct and inverse piezoelectric phenomena, respectively.

2.4. Mechanical problem of free vibration

The local formulation of the free vibration problem of linear elasticity is composed of the following equations

$\begin{array}{l}{\sigma }_{ij,j}=\varrho {{\displaystyle \ddot{u}}}_i,\hfill \\ {\sigma }_{ij}=D_{ijkl}{\epsilon }_{kl},\hfill \\ \begin{array}{cc}{\epsilon }_{ij}=1/2(u_{i,j}+u_{j,i}),& x\in V\end{array}\hfill \end{array}$E13

where

is a density of the elastic body, while ${{\displaystyle \ddot{u}}}_i$ represents the acceleration vector. The displacements are of harmonic character, i.e. u_i = a_i sin ωt, with $\omega$ standing for the unknown natural frequencies of the body and $a_i$ denoting the searched displacement amplitudes.

The boundary conditions are:

$\begin{array}{l}{\sigma }_{ij}{n}_j=\mathrm{0,}x\in S_P,\\ u_i=\mathrm{0,}x\in S_W\end{array}$E14

The variational formulation of the free vibration problem takes advantage of the Hamilton’s principle and reads

${\displaystyle {\int }_V}{D}_{ijkl}{v}_{i,j}{a}_{k,l}dV-{\omega }^2{\displaystyle {\int }_V}\varrho v_i{a}_{i}dV=0$E15

The finite element formulation derived from the above variational functional represents a set of uniform algebraic equations. Such a set possesses a solution if the following characteristic equation is fulfilled:

$det(K_M-{\omega }^{2}M)=0$E16

From this equation $N$ natural frequencies ${\omega }_n,n=\mathrm{1,2,}\dots \mathit{,N}$ can be calculated, where $N$ is the total number of degrees of freedom of the vibrating elastic body. Above, M represents the mass (or inertia) matrix.

For each natural frequency ${\omega }_n$ the nodal vector of displacement amplitudes $q_n^{q,hp}$ can be determined with use of the below finite element equations:

$\begin{array}{c}(K_M-{\omega }_n^{2}M)q_n^{q,hp}=0,\\ q_n^{q,h{p}^T}M{q}_n^{q,hp}=1\end{array}$E17

The second relation above is the normalization condition, completing $N-1$ geometrically independent finite element equations of the first relation.

2.5. Coupled problem of free vibration

We start here with the local (strong) formulation of the undamped vibration problem of the piezoelectric medium

$\begin{array}{l}{\sigma }_{ij,j}+f_i=\varrho {{\displaystyle \ddot{u}}}_i,\hfill \\ d_{i,i}=\mathrm{0,}\hfill \\ {\sigma }_{ij}=D_{ijkl}{\epsilon }_{kl}-C_{kij}{E}_k,\hfill \\ d_i=C_{ikl}{\epsilon }_{kl}+{\gamma }_{ij}{E}_j,\hfill \\ {\epsilon }_{ij}=1/2(u_{i,j}+u_{j,i}),\hfill \\ \begin{array}{cc}{E}_j=-{\phi }_j,& x\in V\end{array}\hfill \end{array}$E18

completed with the following boundary conditions of the coupled electromechanical field

$\begin{array}{ll}{\sigma }_{ij}{n}_j=p_i,\hfill & x\in S_P,\hfill \\ d_i{n}_i=-c,\hfill & x\in S_Q,\hfill \\ u_i=w_i,\hfill & x\in S_W,\hfill \\ \phi =\chi ,\hfill & x\in S_F\hfill \end{array}$E19

In the case of stationary mass $f_i=f_i(x)$ and surface $p_i=p_i(x)$ loadings and charges $c=c(x)$, and the stationary displacement $w_i=w_i(x)$ and electric potential $\chi =\chi (x)$ boundary conditions as well, the problem converts into two independent ones. The first of them is exactly the stationary task of the electro-mechanical equilibrium defined with the local formulation (9) and (10). The corresponding variational and finite element formulations are exactly described with Eqs. (11) and (12), respectively. The solution in displacements $u_i=u_i(x)$ to this stationary problem determines the equilibrium state around which the free vibrations of the piezoelectric are performed. This solution allows the determination of the initial stresses ${\varsigma }_{ij}={\varsigma }_{ij}$(u) which are taken into account in the second problem of free vibration.

The local formulation of the mentioned free vibration problem of the piezoelectric can now be determined in the following way:

$\begin{array}{l}{\sigma }_{ij,j}=\varrho {{\displaystyle \ddot{u}}}_i,\hfill \\ d_{i,i}=\mathrm{0,}\hfill \\ {\sigma }_{ij}={\varsigma }_{ij}+D_{ijkl}{\epsilon }_{kl}-C_{kij}{E}_k,\hfill \\ d_i=C_{ikl}{\epsilon }_{kl}+{\gamma }_{ij}{E}_j,\hfill \\ {\epsilon }_{ij}=1/2(u_{i,j}+u_{j,i}),\hfill \\ \begin{array}{cc}{E}_j=-{\phi }_j,& x\in V\end{array}\hfill \end{array}$E20

where the displacement and the coupled potential fields are: u_i = a_i sin ωt and ϕ = α sin ωt, respectively, with $a_i$ and $\alpha$ standing for the displacement and potential amplitudes.

The above set of differential equations has to be completed with the boundary conditions of the form

$\begin{array}{ll}{\sigma }_{ij}{n}_j=\mathrm{0,}\hfill & x\in S_P,\hfill \\ d_i{n}_i=\mathrm{0,}\hfill & x\in S_Q,\hfill \\ u_i=w_i,\hfill & x\in S_W,\hfill \\ \phi =\chi ,\hfill & x\in S_F\hfill \end{array}$E21

The equivalent variational formulation of the problem reads

$\begin{array}{l}{\displaystyle {\int }_V}{\varsigma }_{ij}{v}_{k,i}{a}_{k,j}dV-{\displaystyle {\int }_V}{D}_{ijkl}{v}_{i,j}{a}_{k,l}dV+{\displaystyle {\int }_V}{C}_{ijk}{v}_{i,j}{\varphi }_{k}dV+{\omega }^2{\displaystyle {\int }_V}\varrho v_i{a}_{i}dV\\ +{\displaystyle {\int }_V}{\gamma }_{ik}{\psi }_i{\varphi }_{k}dV+{\displaystyle {\int }_V}{C}_{ikl}{\psi }_{,i}{a}_{k,l}dV=0\end{array}$E22

while the corresponding finite element equation describing free vibration of the initially stressed piezoelectric medium is

$\begin{array}{c}(K_M+K_G-{\omega }^{2}M)q^{q,hp}-K_C{\phi }^{\rho ,h\pi }=0\\ K_c^T{q}^{q,hp}+K_E{\phi }^{\rho ,h\pi }=0\end{array}$E23

where $q^{q,hp}$ and ${\phi }^{\rho ,h\pi }$ represent the nodal, displacement and electric potential, amplitude degrees of freedom, while $K_G$ stands for the so-called geometric stiffness matrix due to the initial stresses.

As the above set of linear algebraic equations is homogeneous, the solution to it can be obtained if and only if the following characteristic equation:

$det(K_M+K_G+K_C{K}_E^{-1}{K}_C^T-{\omega }^{2}M)=0$E24

is fulfilled. This equation has been obtained after substitution of the second relation (23) into the first one. This allows to remove electric potential amplitudes ${\phi }^{\rho ,h\pi }$ from the combined equation. From this equation $n=\mathrm{1,2,}\dots \mathit{,N}$ natural frequencies ${\omega }_n$ can be obtained, where $N$is the total number of degrees of freedom within the mechanical field.

The corresponding $N$ normalized mode shapes can be obtained from

$\begin{array}{c}(K_M+K_G+K_c{K}_E^{-1}{K}_C^T-{\omega }_n^{2}M)q_n^{q,hp}=0,\\ q_n^{q,h{p}^T}M{q}_n^{q,hp}=1\end{array}$E25

where the normalization condition, the same as in (17), has been applied.

3. Complexity of the modelling

There are three types of complexity considered in this chapter. The first one deals with physical complexity which consists in the presence of more than one physical phenomenon in the problem. The second complexity refers to geometry of the domain under consideration. The geometry is regarded as a complex one if more than one type of geometry is applied. One may deal with a three-dimensional geometry, thin-walled geometry and transition geometry, for example. The third type of complexity is model complexity. In this case, one employs more than one model for description of at least one physical phenomenon under consideration. Combination of these three types of complexity can be regarded as a unique feature of the presented research.

The examples of such complex modelling are electro-mechanical systems composed of geometrically complex elastic structures, joined with the geometrically complex piezoelectric actuators or sensors. In the general case of arbitrary geometry, such systems may require complex mechanical and electro-mechanical description.

3.1. Physical complexity

There are two physical sub-systems present in the considered electro-mechanical systems. The first sub-system concerns bodies subject to elastic deformation and representing structural or machine elements, while the second one concerns piezoelectric bodies acting as actuators or sensors, where the direct or inverse piezoelectric phenomena take place.

3.2. Geometrical complexity

In both, mechanical and piezoelectric, sub-domains we deal with the complex geometry of the structural and piezoelectric elements. In the case of the structural elements, they can be three-dimensional bodies, bounded with surfaces, thin- or thick-shell bodies [4] and solid-to-shell transition bodies. In the case of the piezoelectric members, they can be of three-dimensional, symmetric-thickness or transition character. The shell or symmetric-thickness elements are defined by means of the mid-surface and thickness concepts. In the case of both transition members, we deal with three-dimensional geometry, bounded with surfaces, apart from the boundary to be joined with the shell or symmetric-thickness elements. On this superficial boundary part the mid-surface and the symmetric thickness function have to be defined.

3.3. Model complexity

In the case of the mechanical sub-system complex geometry, the mechanical description may include: three-dimensional elasticity model, hierarchical shell models, the first-order shell model and the solid-to-shell or shell-to-shell transition models. The latter two models allow joining the first-order shell domains with the 3D elasticity and hierarchical shell ones, respectively. In the case of the piezoelectric sub-system, the dielectric model can either represent three-dimensional dielectricity or hierarchical symmetric-thickness dielectric models. The piezoelectric model can be any combination of the listed elastic and dielelectric models.

4. The applied methodology

There are five related aspects of the presented methodology of adaptive hierarchical modelling and adaptive $hp$-finite element analysis of elastic, dielectric and piezoelectric bodies. The first of them is the 3D-based approach proposed in [5]. The second issue, i.e. hierarchical mechanical models were initiated in [6], further developed in [7] and finalized in 3D-based version in [8], while the electric models were introduced in [9] for laminated piezoelectrics and in [10] for dielectrics. General rules of the next aspect of hierarchical approximations were given in [11]. Such approximations for hierarchical shells were developed in [7] and for complex structures in [8]. The latter work generalizes the former attempts given in [12–15]. Approximations for piezoelectric problems were elaborated in [9, 16]. The next issue of error estimation by the equilibrated residual method [17] for 3D elasticity was addressed in [18], for the hierarchical shells in [19], for the 3D-based first-order shells in [20, 21] and for the 3D-based complex structures in [22]. Application of the method in electric problems was proposed in [23]. Finally, adaptivity control with the three-step strategy for simple structures was presented in [24], while adaptivity for simple piezoelectrics was introduced in [25]. The first work was extended onto the 3D-based complex structures in [22].

4.1. 3D-based approach

The applied 3D-based approach [5, 8] lies in application of only three-dimensional degrees of freedom (dofs) regardless of the applied mechanical or electric models. This means that the conventional mid-surface dofs of the shell models, i.e. mid-surface displacements, rotations and other generalized displacement dofs of the mid-surface, are replaced with the equivalent through-thickness displacement dofs similar to the three-dimensional dofs of the 3D elasticity model. Also the mid-surface dofs of the two-dimensional dielectric theory are replaced with the through-thickness electric potential dofs of the three-dimensional dielectrics.

The equivalence of the displacement mid-surface dofs and the through-thickness dofs can be expressed through:

${u^{\prime }}_j={\displaystyle {\displaystyle \sum _{m=0}^I}}{{\xi }^{\prime }}_3^m{s^{\prime }}_j^m={\displaystyle {\displaystyle \sum _{m=0}^I}}{t}_m({{\xi }^{\prime }}_3){u^{\prime }}_j^m$E26

where

represents the local, tangent $(j=\mathrm{1,2})$ and normal $(j=3)$, displacement fields. The terms

and t_m

stand for the mth power of the local, normal coordinate

(

= 0 on the mid-surface) and the mth polynomial through-thickness function of this coordinate, respectively. The mid-surface and through-thickness displacement dof functions are denoted as

and

, respectively, with m = 1,2,…,I and $I$ being the order of the shell theory.

In the case of the electric potential field, the analogous equivalence can be seen in:

$\varphi ={\displaystyle {\displaystyle \sum _{m=0}^J}}{{\xi }^{\prime }}_3^m{\delta }^m={\displaystyle {\displaystyle \sum _{m=0}^J}}{t}_m({{\xi }^{\prime }}_3){\varphi }^m$E27

where $\varphi$ denotes the scalar field of the potential, while ${\delta }^m$ and ${\varphi }^m$ represent mth (m = 1,2,…,J) mid-surface and through-thickness potential dof functions, while $J$ is the order of the two-dimensional dielectric theory.

4.2. Hierarchical models

In the case of the mechanical elastic models the hierarchy M of the 3D or 3D-based models $M$ is [5, 8]:

$M\in M,M=\{3D,MI,RM,3D/RM,MI/RM\}$E28

where 3D represents three-dimensional elasticity, $MI$ denotes hierarchical higher-order shell models, $RM$ is the first-order Reissner-Mindlin shell model, while $3D/RM$ and $MI/RM$ denote the solid-to-shell or shell-to-shell transition models. It is worth mentioning that the second and last models form the following sub-hierarchies:

$\begin{array}{c}MI=\{M\mathrm{2,}M\mathrm{3,}M\mathrm{4,}\dots \},\\ MI/RM=\{M2/RM,M3/RM,M4/RM,\dots \}\end{array}$E29

with $I$ being the order of each particular theory.

The hierarchical character of the above models results from the mentioned order $I$ of the applied 3D-based theories, i.e.

$\begin{array}{l}M=RM\Rightarrow I=\mathrm{1,}K\equiv I\\ M=M2/RM,M3/RM,M4/RM,\dots \Rightarrow I=\mathrm{2,3,4,}\dots ,K=1\\ M=M\mathrm{2,}M\mathrm{3,}M\mathrm{4,}\dots \Rightarrow I=\mathrm{2,3,4,}\dots ,K\equiv I\\ M=3D/RM\Rightarrow I\to \infty ,K=1\\ M=3D\Rightarrow I\to \infty ,K\equiv I\end{array}$E30

Note that for the pure models $(RM,MI,3D)$ $K\equiv I$, while for the transition ones $(MI/RM,3D/RM)$ $K\ne I$.

The hierarchy of 3D-based models possesses the following property:

${\displaystyle \underset{K=\mathrm{1,}I}{{\displaystyle \lim }}}\left({\displaystyle \underset{I\to \infty }{{\displaystyle \lim }}}{\Vert u^{I/K(M)}\Vert }_{U,V}\right)={\Vert u^{3D}\Vert }_{U,V}$E31

guaranteeing that the solutions $u^{I/K(M)}$ obtained with the subsequent models tend in the limit to the solution $u^{3D}$ of the three-dimensional elasticity (the highest model of the hierarchy), when $I\to \infty$. In the above relation the norm of the strain energy $U$ is applied in order to compare the solutions. The norm is defined as follows:

${\Vert u^{I/K(M)}\Vert }_{u,V}=\frac{1}{2}{\displaystyle {\int }_V}{\sigma }^T(u^{I/K(M)})\epsilon (u^{I/K(M)})dV$E32

where $\sigma$ and $\epsilon$ are the stress and strain vectors, respectively.

In the case of the dielectric theories the hierarchy E of the 3D-based models $E$:

$E\in E,E=\{3D,EJ\}$E33

is composed of the three-dimensional theory 3D and the set:

$EJ=\{E\mathrm{1,}E\mathrm{2,}E\mathrm{3,}\dots \}$E34

of the 3D-based hierarchical dielectric models $EJ$, where $J$ is the order of the corresponding dielectric theory.

The hierarchy can be ordered with respect to the order $J$ in the following way:

$\begin{array}{c}E=E\mathrm{1,}E\mathrm{2,}E\mathrm{3,}\dots \Rightarrow J=\mathrm{1,2,3,}\dots ,L\equiv J\\ E=3D\Rightarrow J\to \infty ,L\equiv J\end{array}$E35

and is characterized with the following property:

${\displaystyle \underset{L=\mathrm{1,}J}{{\displaystyle \lim }}}\left({\displaystyle \underset{J\to \infty }{{\displaystyle \lim }}}{\Vert {\varphi }^{J/L(E)}\Vert }_{W,V}\right)={\Vert {\varphi }^{3D}\Vert }_{W,V}$E36

which says that the solutions ϕ^J/L(E) based on the subsequent models of the hierarchy give in the limit $(J\to \infty )$ the solution ϕ^3D of the highest model of the hierarchy, i.e. the model of three-dimensional dielectricity. In the case of the applied pure models, $EJ$ and 3D, L ≡ J.

The norm applied for the model comparisons is based on the electrostatic energy $W$, i.e.

${\Vert {\varphi }^{J/L(E)}\Vert }_{W,V}=\frac{1}{2}{\displaystyle {\int }_V}{d}^T({\varphi }^{J/L(E)})E({\varphi }^{J/L(E)})dV$E37

with $d$ and E being the electric displacement and electric field vectors, respectively.

4.3. Hierarchical approximations

The hierarchical approximations [6, 11] applied to the hierarchy of the 3D-based elastic models can be defined as follows [5, 8]:

$\begin{array}{l}M=RM,I\equiv K=1\Rightarrow u^{q(M),hp}=u^{hp}\\ M\in MI/RM,I\ge \mathrm{2,}K=1\Rightarrow u^{q(M),hp}=u^{hpq/hp}\\ M\in MI,I\equiv K\ge 2\Rightarrow u^{q(M),hp}=u^{hpq}\\ M=3D/RM,I\to \infty ,K=1\Rightarrow u^{q(M),hp}=u^{hpp/hp}\\ M=3D,I\equiv K\to \infty \Rightarrow u^{q(M),hp}=u^{hpp}\end{array}$E38

Above we applied the general notation, $q(M)$,$hp$, valid for any hierarchical model, and the equivalent symbols of the approximations for each particular model (on the right side of the above relations), where $h$ represents the generalized size of the element, $p$ denotes the longitudinal order of approximation and $q\equiv q(M)$ is the transverse order of approximation, equivalent to the order of the theory, i.e. $q(M)\equiv I(M)$. The specific approximations are either two-dimensional $(hp)$ or three-dimensional $(hpq,hpp)$ or mixed $(hpq/hp,hpp/hp)$.

The above approximations lead in the limit $(1/h\to \infty ,p\to \infty )$ to the exact solution of the appropriate mechanical model $M$ of the order $I/K$, i.e.

${\displaystyle \underset{1/h,p\to \infty }{{\displaystyle \lim }}}{\Vert u^{q(M),hp}\Vert }_{U,\cup {\displaystyle \stackrel{e}{V}}}={\Vert u^{I/K(M)}\Vert }_{U,V}$E39

Combination of the models defined with (30) and the corresponding approximations of (38) leads to the hierarchical numerical models of the following property:

${\displaystyle \underset{K=1,I}{\lim }}\left[{\displaystyle \underset{I\to \infty }{\lim }}\left({\displaystyle \underset{1/h,p\to \infty }{\lim }}{\Vert u^{q(M),hp}\Vert }_{U,\cup {\displaystyle \stackrel{e}{V}}}\right)\right]={\displaystyle \underset{K=1,I}{\lim }}\left({\displaystyle \underset{I\to \infty }{\lim }}{\Vert u^{I/K(M)}\Vert }_{U,V}\right)={\Vert u^{3D}\Vert }_{U,V}$E40

As it can be seen above, the solutions to these models tend in the limit to the exact solution of the highest mechanical model $(3D)$.

By analogy, the numerical approximations of the 3D-based dielectric models are:

$\begin{array}{l}E\in EJ,J\equiv L\ge 2\Rightarrow {\varphi }^{\rho (E),h\pi }={\varphi }^{h\pi \rho }\\ E=3D,J\equiv L\to \infty \Rightarrow {\varphi }^{\rho (E),h\pi }={\varphi }^{h\pi \pi }\end{array}$E41

Above, for the sake of further considerations, the longitudinal approximation order $p$ is denoted as $\pi$ and the transverse one $q$ as $\rho \equiv \rho (E)$, where $\rho (E)\equiv J(E)$. The principal property of the above electric potential approximations reads

${\displaystyle \underset{1/h,\pi \to \infty }{{\displaystyle lim}}}{\Vert {\varphi }^{\rho (E),h\pi }\Vert }_{w,\cup {\displaystyle \stackrel{e}{V}}}={\Vert {\varphi }^{J/L(E)}\Vert }_{W,V}$E42

while for the hierarchical numerical dielectric models, being the combination of the electric models (35) and the above approximations (41), one has

${\displaystyle \underset{L=\mathrm{1,}J}{\lim }}\left[{\displaystyle \underset{J\to \infty }{\lim }}\left({\displaystyle \underset{1/h,\pi \to \infty }{\lim }}{\Vert {\varphi }^{\rho (E),h\pi }\Vert }_{W,\cup {\displaystyle \stackrel{e}{V}}}\right)\right]={\displaystyle \underset{L=\mathrm{1,}J}{\lim }}\left({\displaystyle \underset{J\to \infty }{\lim }}{\Vert {\varphi }^{J/L(E)}\Vert }_{W,V}\right)={\Vert {\varphi }^{3D}\Vert }_{W,V}$E43

where the limit solutions

and

^3D are present.

4.4. Error estimation

It is well known [17–19, 22] that the equilibrated residual methods of error estimation require solution of the following approximated local (element) problems of mechanical equilibrium:

$\stackrel{e}{\text{B}}(u^{Q,HP}{v}^{Q,HP})-\stackrel{e}{\text{L}}(v^{Q,HP})-{{\displaystyle {\int }_{\stackrel{e}{S}\backslash (S_P\cup S_D)}(v^{Q,HP})}}^T⟨\stackrel{e}{r}(u^{q,hp})⟩\stackrel{e}{dS}=0$E44

where ${\displaystyle \stackrel{e}{\text{B}}}$ and ${\displaystyle \stackrel{e}{\text{L}}}$ are the bilinear and linear forms representing virtual strain energy and the virtual work of the external forces. The last term stands for the virtual work of the inter-element stress fluxes

. The element size in the local problems is denoted as $H$, the longitudinal approximation order is $P$, while the transverse one is $Q\equiv Q(M)$.

The above three terms can be defined as follows

E45

with $D$ being the elasticity constants matrix. Thus, Eq. (44) can be written in the language of finite elements:

${{\displaystyle \stackrel{e}{k}}}_M{{\displaystyle \stackrel{e}{q}}}^{Q,HP}={{\displaystyle \stackrel{e}{f}}}_V+{{\displaystyle \stackrel{e}{f}}}_S+{{\displaystyle \stackrel{e}{f}}}_R$E46

where

denote element stiffness matrix, and element forces vectors due to mass, surface and inter-element stress loadings, while

and

represent solution and admissible displacement dof vectors in the local problem. The solution vector is then utilized in the error estimation. The collection of the element solutions obtained this way leads to the global error estimate which constitutes the upper bound of the true error [18, 19, 22].

The above approach can be extended onto the mechanical problem of free vibration by replacing the linear form of (44) by the virtual work of the inertia forces [26]. These forces are:

, where

is the element mass (or inertia) matrix, while the natural frequencies ${\omega }_n$ and the vector of the element displacement amplitudes

are taken from the global problem solution. These forces replace

and

in (45) and (46). Note that now the element solutions, analogous to this obtained from (46), do not guarantee the upper boundedness of the global error by its residual estimate.

In the case of electrostatic dielectricity, the local problems of the equilibrated residual method take the form [23]:

${\displaystyle \stackrel{e}{\text{b}}}({\varphi }^{P,H\text{Π}}{\psi }^{P,H\text{Π}})-{\displaystyle \stackrel{e}{1}}({\psi }^{P,H\text{Π}})-{\displaystyle {\int }_{{\displaystyle \stackrel{e}{S}}\backslash (S_Q\cup S_F)}}{\psi }^{P,H{\text{Π}}^T}⟨{\displaystyle \stackrel{e}{h}}({\varphi }^{\rho ,h\pi })⟩{\displaystyle \stackrel{e}{dS}}=0$E47

with the bilinear and linear forms,

and

, representing virtual electrostatic energy and the virtual work of external charges. The right term above is equal to the virtual work of the equilibrated charge flux

. The indices $H,\text{Π},\text{P}\equiv \text{P}(E)$ stand for the element size and the element longitudinal and transverse approximation orders.

Taking into considerations the below definitions:

E48

where γ is the dielectricity constants matrix, one can transform (47) into the following local finite element equation

${{\displaystyle \stackrel{e}{k}}}_E{{\displaystyle \stackrel{e}{\phi }}}^{P,H\text{Π}}={{\displaystyle \stackrel{e}{f}}}_Q+{{\displaystyle \stackrel{e}{f}}}_H$E49

Above

denote element dielectricity matrix, and element vectors due to surface and inter-element charges, while

and

are the solution and admissible electric potential dof vectors of the local problem. This element solution is applied to the error estimation of the dielectric problems. The global error estimate obtained with use of the above element solutions upper-bounds the true global error.

4.5. Adaptive strategy

The adaptive strategy for the complex problems of elasticity is based on the Texas Three-Step Strategy [24]. The latter strategy consists in solution of the global problem thrice on the so-called initial, intermediate and final meshes. The intermediate mesh is obtained from the initial one based on the initial mesh estimated values of element errors and the hp convergence theory relating these errors to the discretization parameter $h$. Thanks to this the h-adaptation is performed. The final mesh is obtained from the intermediate one in the analogous way through p-refinement. This process takes advantage of the intermediate mesh estimated errors and the relation between these errors and the discretization parameter p. The original strategy can be enriched with the fourth step [5, 27], called the modification one, which is performed on the initial mesh and is applied when the unpleasant numerical phenomena due to the improper solution limit, numerical locking or boundary layers appear in the mechanical problem. The purpose of this additional step is to get rid of the numerical consequences of the mentioned phenomena before the error-controlled hp-adaptivity is started. The control of this step is based on the sensitivity analysis of the local (element) problems solutions to these three phenomena.

The above four-step strategy can be easily extended onto the complex dielectric problems. It has to account only for the boundary layer phenomenon within the modification step, as two other phenomena do not appear in the problems of dielectrics. Other steps do not change.

5. Problems within the methodology

In this section of the chapter we consider the main difficulties in generalization of the hierarchical models, hierarchical approximations, equilibration residual method of error estimation and three- or four-step adaptive strategies, presented in the previous sections for the problems of elasticity and dielectricity, onto coupled problems of piezoelectricity.

5.1. Hierarchical model and approximation issues

The first task here is to compose the elastic and dielectric models, defined with (28) and (33), respectively, into one consistent hierarchy of the piezoelectric media. Our proposition on how to perform this task follows from the main feature of both component hierarchies which are characterized with the independently changing orders $I$ and $J$ of the elastic and dielectric models, as shown in (30) and (35). Because of this we can propose the definition

$P\in P,P=\{(M,E):M\in M,E\in E\}$E50

which determines the hierarchy P of piezoelectric models $P$ as composed of all combinations (M,E) of the elastic $M$ and dielectric $E$ models. Even though the following property:

${\displaystyle \underset{K=\mathrm{1,}I}{\lim }}\left\{{\displaystyle \underset{I\to \infty }{\lim }}\left[{\displaystyle \underset{L=\mathrm{1,}J}{\lim }}\left({\displaystyle \underset{J\to \infty }{\lim }}{\Vert (u^{I/K(M)},{\varphi }^{J/L(E)})\Vert }_{Z,V}\right)\right]\right\}={\Vert (u^{3D},{\varphi }^{3D})\Vert }_{Z,V}$E51

is valid in the limit, the ordering of the coupled electro-mechanical solutions measured in the energy norm is not unique due to different signs of the strain (or elastostatic) $U$, electrostatic $W$ and coupling $C$ parts of the co-energy: $Z=W-U+C$, i.e.

${\Vert (u^{I/K(M)},{\varphi }^{J/L(E)})\Vert }_{Z,V}=\frac{1}{2}{\displaystyle {\int }_V}|E^T\gamma E-{\epsilon }^{T}D\epsilon +2{\epsilon }^{T}CE|dV$E52

where $\epsilon =\epsilon (u$^I/K(M)) and E = E(ϕ^J/L(E)), while C is the coupling (piezoelectric) constants matrix.

The most general case of the $hp$-approximation within the coupled field of displacements and electric potential may include totally independent $h$ and $p$ approximations within the mechanical and electric fields. However, such an approach requires the vice-versa projections of the displacements and potential solutions between the independent $h$-meshes of each field. The approach results in the additional projection error which should be included in the error estimation. Because of this we propose the simplified approach which consists in application of the common $h$ mesh and the independent $p$ and $\pi$ approximations within the displacements and electric potential fields. This assumption leads to the following limit property of the coupled solution

${\displaystyle \underset{1/h,p\to \infty }{\lim }}\left({\displaystyle \underset{1/h,\pi \to \infty }{\lim }}{\Vert (u^{q(M),hp},{\varphi }^{\rho (E),h\pi })\Vert }_{Z,\cup {\displaystyle \stackrel{e}{V}}}\right)={\Vert (u^{I/K(M)},{\varphi }^{J/L(E)})\Vert }_{Z,V}$E53

It should be noticed that the monotonic character of the co-energy of the consecutive approximate solutions is not guaranteed here due to the coupled character of the electro-mechanical field.

Practical realization of the above concepts of hierarchical modelling and approximations is implemented by means of the $hpq$- and $h\pi \rho$-adaptive piezoelectric finite elements, presented in Figure 1. As said before, $h$ represents the assumed characteristic element size, common for both fields, while $p$ and $\pi$ denote independent longitudinal approximation orders within the displacements and electric potential fields. The transverse approximation orders $q$ and $\rho$ of the mechanical and electric fields, respectively, are equivalent to the independent hierarchical orders $I$ and $J$ of the elastic and dielectric models, i.e. $q\equiv I,\rho \equiv J$.

The normalized versions of the prismatic adaptive elements are presented above, where the 3D-based solid (or hierarchical shell), transition (an example of) and first-order shell mechanical elements are combined with the three-dimensional (or 3D-based symmetric-thickness) dielectric elements. In the figure the normalized coordinates are defined as ${\xi }_{\mathrm{1,}}{\xi }_{\mathrm{2,}}{\xi }_3\in (\mathrm{0,1})$, while the vertex, mid-edge, mid-base, mid-side and middle nodes of mechanical character are denoted as either $a_1,a_2,\dots ,a_{21}$ or $a_1,a_2,\dots ,a_{18}$ or $a_1,a_2,\dots ,a_{14}$, while the corresponding electric ones are marked with $b_1,b_2,\dots ,b_{21}$.

Figure 1.

Piezoelectric solid/three-dimensional (top), transition/three-dimensional (middle), shell/three-dimensional (bottom) elements.

5.2. Problems within error estimation

In order to generalize the equilibrated residual method for piezoelectricity the local elastic (44) and dielectric (47) problems have to be replaced with the coupled stationary problems describing the electro-mechanical equilibrium:

E54

where

and

, while the above coupling forms are determined as follows:

E55

Definitions (45), (48) and (55) allow to rewrite (54) in the finite element language:

E56

The main disadvantage of the local solutions

and

obtained from the above set of equations is that they do not lead to the upper bound of the global approximation error of the global problem (12) due to the coupled character of piezoelectricity. Note that the upper bound property is present in the cases of pure elasticity and pure dielectricity. In these circumstances we propose the simplified approach which consists in decoupling of the local mechanical and electric fields:

E57

The simplified solutions of the above two decoupled problems is suggested for serving the approximated error assessment of the coupled field.

In the case of the coupled free-vibration problem we propose to apply the same approach as for the purely elastic vibrations. This means that the inertia forces introduced in Subsection 4.4 have to replace the volume and surface forces in the first equations (56) and (57). The element constitutive stiffness matrix

in both equations has to be changed for the following sum:

 + 

, where

stands for the element geometric stiffness matrix due to the initial stresses which correspond to the electro-mechanical equilibrium. Because of this one may skip the electric charges

in the second equations (56) and (57).

5.3. Adaptivity control matters

The main problem with the piezoelectricity in the context of adaptivity control of the three- or four-step $hp$-adaptive procedure is that the convergence theorem for the finite element approximation of the coupled piezoelectric field is not at one’s disposal. In these circumstances we propose to use the $hp$-convergence exponents of the elasticity and dielectricity problems for the mechanical and electric components of the coupled field.

With this assumption in mind we generalize the adaptive scheme applied successfully for the complex elastic structures [28]. It can also be adopted for the dielectricity problems. The generalization is presented in Figure 2 in the form of a block diagram, where $i=\mathrm{1,2,3,4}$ represent the consecutive steps of the algorithm—the initial, modification, intermediate and final ones.

Figure 2.

Four-step hp-adaptive scheme for complex piezoelectric or piezoelectric-elastic problems.

6. Conclusions

In the presented chapter we showed how to generalize hierarchical modelling, $hp$-adaptive hierarchical approximations, the equilibrated residual method of error estimation and the four-step adaptive procedure, originally applied to elasticity and possible also in dielectricity, for the coupled pieozoelectric problems.

We suggest to combine the hierarchical models for complex elastic structures with the analogous models of dielectricity in order to obtain all coupled combinations of the component mechanical and electric models. In the case of the approximations, we suggest to apply the common $h$-mesh and independent p and $\pi$ approximations for the displacement and electric potential fields.

In the error estimation, we propose to decouple the mechanical and electric fields in the local problems of the residual approach.

Our four-step $hp$-adaptive algorithm is based on the convergence theorems for the purely elastic and purely dielectric problems. These theorems are applied to the intermediate and final meshes generation.

Acknowledgments

The support of the Polish Scientific Research Committee (now the National Science Centre) under the research grant no. N N504 5153 040 is thankfully acknowledged.

1. $a_i$ – displacement amplitude vector components $(i=\mathrm{1,2,3})$,

2. $a_j$ – a node with displacement degrees of freedom (dofs), $j=\mathrm{1,2,\dots ,21,}$

3. ${\displaystyle \stackrel{e}{\text{b}}}$ – a bilinear form representing the element virtual electrostatic energy,

4. $b_j$ – a node with electric potential dofs, $j=\mathrm{1,2,\dots ,21,}$

5. ${\displaystyle \stackrel{e}{\text{B}}}$ – a bilinear form representing the element virtual strain energy,

6. c – the electric surface charge,

7. 
   – a bilinear form representing the element coupling electro-mechanical energy,

8. $C_{ikl}$ – components of the piezoelectricity tensor (i,k,l = 1,2,3),

9. $C$ – the piezoelectric constants matrix,

10. $d_i$ – components of the electric displacement vector $(i=\mathrm{1,2,3})$,

11. $d$ – the electric displacement vector,

12. $D_{ijkl}$ – elasticity tensor components, i,j,k,l = 1,2,3,

13. $D$ – the elasticity constants matrix,

14. $e$ – the number of an element,

15. $E-$a dielectric model,

16. E – the hierarchy of dielectric models,

17. $E_j$ – electric field vector components $(j=\mathrm{1,2,3})$,

18. $E$ – the electric field vector,

19. $F_i$ – mass load vector components $(i=\mathrm{1,2,3})$,

20. $f$ – the mass load vector,

21. 
    – the inter-element charge vector,

22. 
    – the element mass forces vector,

23. 
    – the element charge vector,

24. 
    – the inter-element stress forces vector,

25. 
    – the element surface forces vector,

26. $F_M$ – the global mass forces vector,

27. $F_S$ – the global surface forces vector,

28. $F_Q$ – the global charge vector,

29. $h$ – measure of the element size,

30. 
    – the equilibrated inter-element charge,

31. $H$ – a measure of the element size in the local problem,

32. $I$ – the order of the hierarchical mechanical model,

33. $J$ – the order of the hierarchical dielectric model,

34. 
    – the element piezoelectric (coupling) matrix,

35. 
    – the element dielectric matrix,

36. 
    – the element geometric stiffness matrix,

37. 
    – the element (constitutive) stiffness matrix,

38. $K_C$ – the global piezoelectric (coupling) matrix,

39. $K_E$ – the global dielectric matrix,

40. $K_G$ – the global geometric stiffness matrix,

41. $K_M$ – the global (constitutive) stiffness matrix,

42. 
    – a linear form equal to the element virtual work of the external charges,

43. 
    – a linear form representing the element virtual work of the external forces,

44. 
    – the element mass (inertia) matrix,

45. $m$ – the number of a dof function,

46. $M$ – mechanical model,

47. $\text{M}$ – the hierarchy of mechanical models,

48. $M$ – the global mass (inertia) matrix,

49. $n$ – the number of a natural frequency,

50. $n_j$ – vector components $(j=\mathrm{1,2,3})$ of the normal to the body surface,

51. $N$ – the global number of dofs (degrees of freedom),

52. $p$ – the element approximation order or the longitudinal order of approximation for displacements,

53. $p_i$ – surface load vector components,

54. $P$ – the longitudinal order of approximation for element displacements,

55. $P$ – piezoelectric model,

56. $\text{P}$ – the hierarchy of piezoelectric models,

57. $p$ – the surface load vector,

58. $q$ – the transverse order of approximation for displacements,

59. $q^{q,hp}$ – the global displacement (or displacement amplitude) dofs vector,

60. 
    – the element displacement (or displacement amplitude) dofs vector,

61. $Q$ – the local problem transverse order of approximation for displacements,

62. 
    – a vector of the equilibrated inter-element stresses,

63. 
    – a mid-surface displacement dof function,

64. $S$ – the surface of a body (or medium),

65. $S_F$ – the body surface part with given electric potential,

66. $S_P$ – the loaded part of the body surface,

67. $S_Q$ – the electrically charged part of the body (or medium) surface,

68. $S_W$ – the body surface part with given displacements,

69. t – time,

70. $t_m$ – the mth polynomial through-thickness function,

71. $u$ – the vector of global displacements,

72. $u_i$ – the global displacement components, $i=\mathrm{1,2,3}$,

73. $u^{Q,HP}$ – the displacement vector in the local (element) problem,

74. 
    – local components $(j=\mathrm{1,2,3})$ of the displacement vector,

75. 
    – a through-thickness displacement dof function,

76. ${{\displaystyle \ddot{u}}}_i$ – acceleration vector components, $i=\mathrm{1,2,3}$,

77. $v_i$ – admissible displacement (or its amplitude) vector components $(i=\mathrm{1,2,3})$,

78. $v^{Q,HP}$ – the admissible displacement vector in the local (element) problem,

79. 
    – the admissible displacement dof vector of an element,

80. $V$ – volume of a body (or medium),

81. $w_i$ – given displacement vector components $(i=\mathrm{1,2,3})$ on the body surface,

82. $x$ – global Cartesian coordinate vector,

83. ${\gamma }_{ij}$ – components of the dielectricity tensor $(i\mathit{,j}=\mathrm{1,2,3})$,

84. γ – the dielectricity constants matrix,

85. ${\delta }^m$ – the $m\text{th}$ mid-surface electric potential dof function,

86. ${\epsilon }_{ij}$ – strain tensor components $(i\mathit{,j}=\mathrm{1,2,3})$,

87. $\epsilon$ – the strain vector,

88. ${\xi }_i$ – normalized coordinates of an element $(i=\mathrm{1,2,3})$,

89. 
    – the local normal coordinate,

90. $\pi$ – the longitudinal approximation order of electric potential,

91. $\text{Π}$ – the element longitudinal approximation order of electric potential,

92. $\rho$ – the transverse approximation order of electric potential,

93. 
    – mass density of the body,

94. $\text{P}$ – the element transverse approximation order of electric potential,

95. ${\varsigma }_{ij}$ – initial stress tensor components (i,j = 1,2,3),

96. ${\sigma }_{ij}$ – stress tensor components (i,j = 1,2,3),

97. $\sigma$ – the stress vector,

98. ${\phi }^{\rho ,h\pi }$ – the global electric potential (or potential amplitude) dofs vector,

99. 
    – the local (element) electric potential dofs vector,

100. $\varphi$ – the electric potential,

101. ${\varphi }^m$ – the mth through-thickness electric potential dof function,

102. $\chi$ – the given electric potential on the body surface,

103. $\psi$ – the admissible electric potential,

104. 
     – the admissible electric potential dofs vector of an element,

105. $\omega$ – a natural frequency,

106. ${\omega }_n$ – the nth natural frequency