Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

Abstract

In this chapter, theoretical and implementation details of the algorithms of hierarchical modeling and hierarchical hp-approximations, residual error estimation methods and four-step adaptive procedures are considered in the context of their application to modeling and simulation of the problems of elasticity, dielectricity and piezoelectricity. In the hierarchical modeling, 3D-based hierarchical elastic and dielectric models are applied. The adaptive discretization process is based on the hierarchical shape functions and the constrained approximations. In the error estimation, the equilibrated residual method is applied, which serves the total and approximation error assessment. These errors control the model and hp-adaptivity. In the case of adaptive algorithms, four-step procedure is utilized. It includes global solutions on the initial mesh, mesh modified in order to remove some undesired numerical phenomena, the intermediate h-refined mesh and the final (or target) p-enriched mesh. Examples demonstrating the effectivity of the mentioned modeling and approximation, error estimation and adaptivity control parts of the overall simulation algorithm in the three classes of problems are presented.

1. Introduction

This chapter constitutes a continuation and extension of the previous work [1] on theoretical and implementation difficulties in application of the adaptive hierarchical modeling and $\mathit{hp}$-adaptive finite element analysis to elasticity, dielectricity and piezoelectricity. In the cited work, the 3D-based elastic, dielectric and piezoelectric hierarchies of models were elucidated. These models are based on either three-dimensional theories or reduced models polynomially constrained through the thickness. In the mentioned work, also the hierarchical approximations for the three classes of hierarchical models are presented. The rules for ordering the hierarchical models and approximations are described. Then, the a posteriori error estimation, based on the equilibrated residual method (ERM) applied to the three classes of problems, is presented. The similarities and differences between the element (local) problems necessary for the element error estimation for these three cases are addressed. Finally, the three- and four-step error-controlled adaptive procedures for the three classes of problems are proposed. The procedures require the global problem solution on the initial, modified, intermediate ($h$-refined) and final ($p$-enriched) meshes.

In this chapter, attention is paid to the effectivity of the algorithms for adaptive modeling and simulation of three considered classes of physical phenomena, that is, elasticity, dielectricity and piezoelectricity. Effectivity of hierarchical approximations within elastic, dielectric and piezoelectric media is compared. For this purpose, convergence curves for the analogous model problems within three mentioned classes of problems are generated and assessed. Also, the exemplary comparative results of the approximations are presented for these three classes of model problems. In the case of the error estimation, the global and local (element) effectivity indices for the total, approximation and modeling errors, where the latter is the difference of the former two, in the exemplary model problems of elasticity, dielectricity and piezoelectricity are calculated and compared. The exemplary distribution of the element error indicators and the global values of the error estimators for the model problems of three classes are presented and compared. In the case of the adaptive procedures, the model- and $\mathit{hpq}$-adaptive algorithms, where $h$ represents the element size parameter, while $p$ and $q$ stand for the element longitudinal and transverse orders of approximation, are of our interest. These algorithms are controlled with the estimated values of the modeling, approximation and total errors. In order to check the effectivity of these algorithms, results necessary for the obtainment of the $\mathit{hp}$-adaptive convergence curves for the mentioned three model problems of elasticity, dielectricity and piezoelectricity are produced. The convergence is assessed in the context of obtainment of the target values of the errors in subsequent steps of the adaptive calculations for three classes of problems. Also, the comparative results of the adaptive solutions of the model problems of three classes are presented.

1.1. Research objectives

The main objective of this research is to demonstrate the effectivity of our generalizing algorithms [1] adapted, modified or developed for the problems of elasticity, dielectricity and piezoelectricity. Also, the issue of comparison of the corresponding effectivities for these three classes of problems is of our interest. In relation to these objectives, the presented general approach to adaptive modeling and simulation is numerically tested in the context of the approximation algorithms, error estimation algorithms and adaptivity control algorithms as well.

1.2. State-of-the-art issues

In this brief survey, the issues of hierarchical modeling, hierarchical approximations, error estimation and adaptivity control are addressed. The survey is limited to the numerical techniques used in this chapter and the papers directly utilized for this research—no general overview of the four mentioned issues is presented. The interested readers can find such overviews in some of the publications cited below.

The 3D-based hierarchical shell models utilizing three-dimensional degrees of freedom (dofs) and conforming to higher order shell theories were firstly proposed in [2] and repeated in [3]. The conventional hierarchical shells, employing mid-surface dofs, were proposed in [4]. The 3D-based approach was extended onto the first-order shell and shell-to-shell theories in [3, 5]. The latter works also extend the 3D-based hierarchical modeling onto 3D elasticity and solid-to-shell transition models. The author of this chapter is not aware of any hierarchical models of linear dielectricity. Some hierarchic piezoelectric models were presented in [6] in the context of multilayered plate structures. Suggestions on introduction of the 3D-based hierarchical dielectric and piezoelectric models were formulated in [1, 7].

The hierarchical and constrained approximations necessary for $p$- and $h$-adaptivity, respectively, are adopted in our work and were proposed in [8]. Hierarchical approximations for conventional shells were developed in [4], for 3D-based shells in [2, 9] and for complex structures in [5]. The last paper collects partial results presented in [9, 10, 11, 12]. Classical and hierarchical approximations for piezoelectric problems were elaborated in [6, 13]. Hierarchical approximations for the complex 3D or 3D-based hierarchical models of dielectrics and piezoelectrics were proposed in the works [1, 7].

The general considerations on error estimation based on the equilibrated residual method can be found in [14]. Application of this method to 3D elasticity was described in [15]. The method was also applied to the hierarchical shells of conventional character in [16]. The analogous approach for the 3D-based first-order shells was developed in [17, 18]. The method was also utilized to error estimation in the 3D-based complex structures [19]. Application of the method to dielectric and piezoelectric problems was suggested in [1, 20].

Finally, adaptivity control by means of the three-step strategy for simple structures was presented in [21]. Within this strategy, three subsequent meshes are generated—initial, intermediate (or $h$-refined) and target (or $p$-enriched) ones. The method was applied to adaptive analysis of conventional hierarchical models of shell- and plate-like structures in [16]. In that work, the third step is split into two, that is, $q$ and $p$ enrichments are performed in sequence. The original three-step strategy was then extended in [3] by addition of the fourth step in which the mesh is modified to get rid of the numerical consequences of the improper solution limit, numerical locking and edge effect. The model adaptivity is performed along with the $h$-step and $p$ and $q$ enrichments are performed simultaneously. Such a four-step adaptive strategy is applied to modeling and simulation of complex elastic structures in [19]. Adaptive simulation in electric or electromechanical problems is less advanced. Adaptivity for simple piezoelectrics was introduced in [22]. Application of the three- or four-step strategies to the analysis of simple and complex dielectrics and piezoelectrics was suggested in [1].

1.3. Novelty of the research

The main novelty of the presented research consists in application of the chosen techniques of hierarchical modeling and approximation, error estimation and adaptivity control, effective in the adaptive modeling and simulation of the elasticity problems, to the adaptive analysis of dielectric and piezoelectric phenomena.

The novelty of this particular chapter is the direct comparison of the robustness of the modeling and simulation algorithms of the coupled problem of piezoelectricity and the problems of pure elasticity and pure dielectricity.

2. Model problems

The following model problems are considered in this chapter: the linear static problem of elasticity, the linear electrostatic problem and the linear problem of stationary piezoelectricity. For each of the model problems, the appropriate finite element formulation is presented. For this purpose, the standard engineering matrix notation is applied.

2.1. Elastostatics

Here, the problems of a three-dimensional (solid) and 3D-based shell or solid-to-shell bodies are considered. Such problems were presented in [1]. In that work, the local (strong) and variational (weak) formulations of the problems are given. These formulations take advantage of the former considerations from [2, 23, 24] and are repeated in [3]. Using the variational formulation presented therein, one can derive the global finite element equations of the problem under consideration and write them in the following form:

$${\mathit{K}}_M{\mathit{q}}^{q,\mathit{hp}}={\mathit{F}}_V+{\mathit{F}}_S$$E1

where ${\mathit{K}}_M$ is the global stiffness matrix, while ${\mathit{F}}_V$ and ${\mathit{F}}_S$ represent the global vectors of the volume and surface nodal forces. The vector ${\mathit{q}}^{q,\mathit{hp}}$ stands for the global displacement degrees of freedom (dof), corresponding to $\mathit{hpq}$ approximation, and is composed of the element (local) displacement dof vectors ${\displaystyle \stackrel{e}{\mathit{q}}}$ of the elements $e=1,2,\dots ,E$, where $E$ is the total number of elements within an elastic body. These vectors are defined later in this chapter.

The global stiffness matrix is composed (aggregated) of the element stiffness matrices of the form

$${\displaystyle \stackrel{e}{\mathit{k}}}={\int }_0^1{\int }_0^1{\int }_0^{-{\xi }_2+1}{\displaystyle \stackrel{e}{{\mathit{B}}^T}}\mathit{D}{\displaystyle \stackrel{e}{\mathit{B}}}det\left(\mathit{J}\right)d{\xi }_{1}d{\xi }_{2}d{\xi }_3$$E2

where $\mathit{D}$ denotes the elastic constants matrix, ${\displaystyle \stackrel{e}{\mathit{B}}}$ represents the strain-displacement matrix, and $det\left(\mathit{J}\right)$ is the Jacobian matrix determinant. The limits and coordinates of the integration correspond to the normalized coordinates ${\xi }_i$, $i=1,2,3$ of the prismatic elements applied in [1]. The specific forms of the strain-displacement matrix can be found in the works [9, 10, 12] for the 3D-based versions of the prismatic solid (and hierarchical shell), first-order shell and solid-to-shell (and shell-to-shell) adaptive elements, respectively.

The nodal mass forces vector can be defined in the standard way

$${{\displaystyle \stackrel{e}{\mathit{f}}}}_M={\int }_0^1{\int }_0^1{\int }_0^{-{\xi }_2+1}{\displaystyle \stackrel{e}{{\mathit{N}}^T}}\mathit{f}det\left(\mathit{J}\right)d{\xi }_{1}d{\xi }_{2}d{\xi }_3,$$E3

where ${\displaystyle \stackrel{e}{\mathit{N}}}$ and $\mathit{f}$ represent the element shape functions matrix and the mass loading vector, respectively.

The element nodal forces vector due to the surface traction $\mathit{p}$ can be defined in the following two forms

$$\begin{array}{c}{\displaystyle \stackrel{e}{{\mathit{f}}_S}}={\int }_0^1{\int }_0^{-{\xi }_2+1}{\displaystyle \stackrel{e}{{\mathit{N}}^T}}\mathit{p}\mathrm{wsp}\left(\mathit{J}\right)d{\xi }_{2}d{\xi }_1\\ {\displaystyle \stackrel{e}{{\mathit{f}}_S}}={\int }_0^1{\int }_0^1{\displaystyle \stackrel{e}{{\mathit{N}}^T}}\mathit{p}\mathrm{wsp}\left(\mathit{J}\right)d{\xi }_{3}d{\eta }_i\end{array}$$E4

corresponding to the bases and sides of the prismatic element. Above the element, bases and sides are defined with the normalized longitudinal coordinates ${\xi }_j$, $j=1,2$, or the transverse normalized coordinate ${\xi }_3$ and the coordinates ${\eta }_i$, $i=1,2,3$ tangential to the sides of the element [3, 9]. The term $\mathrm{wsp}\left(\mathit{J}\right)$ is the coefficient defined with the components of the Jacobian matrix $\mathit{J}$ (see [3, 9] again).

2.2. Electrostatics

The general formulations of the problems of electrostatics can be found in [25]. Here, classical linear dielectric models are applied to such problems. The local and variational formulations for this case was presented in [1] for any 3D or 3D-based geometry (bulky, symmetric-thickness or transition ones). The corresponding finite element equations read:

$${\mathit{K}}_E{\mathit{\phi }}^{\rho ,\mathit{h\pi }}={\mathit{F}}_Q$$E5

In Eq. (5), ${\mathit{K}}_E$ represents the global characteristic matrix of dielectricity, while ${\mathit{F}}_Q$ stands for the global characteristic electric charges nodal vector. The vector ${\mathit{\phi }}^{\rho ,\mathit{h\pi }}$ is the unknown global nodal vector of electric potentials. This vector definition results from the applied $\rho ,\mathit{h\pi }$-approximation, where $\rho$ and $\pi$ represent the transverse and longitudinal orders of approximation. The global potential vector is composed of the element potential vectors ${\displaystyle \stackrel{e}{\mathit{\phi }}}$, which are described later in this chapter. The global matrix ${\mathit{K}}_E$ is the result of summation of the element contributions

$${\displaystyle \stackrel{e}{{\mathit{k}}_E}}={\int }_0^1{\int }_0^1{\int }_0^{-{\xi }_2+1}{\displaystyle \stackrel{e}{{\mathit{b}}^T}}\gamma {\displaystyle \stackrel{e}{\mathit{b}}}det\left(\mathit{J}\right)d{\xi }_{1}d{\xi }_{2}d{\xi }_3$$E6

with $\gamma$ and ${\displaystyle \stackrel{e}{\mathit{b}}}$ denoting the electric (or permittivity) constants matrix and the matrix of the relation between the electric field components and the nodal electric potentials (or shortly field-potential matrix). The specific form of the latter matrix in the case of the prismatic element can be found in the work [26].

The nodal electric charges vector of the element $e$ has to be defined in a different way on the prismatic element bases and sides, that is,

$${\displaystyle \stackrel{e}{{\mathit{f}}_Q}}={\int }_0^1{\int }_0^{-{\xi }_2+1}{\displaystyle \stackrel{e}{{\mathit{n}}^T}}c\mathrm{wsp}\left(\mathit{J}\right)d{\xi }_{2}d{\xi }_1$$E7

and

$${\displaystyle \stackrel{e}{{\mathit{f}}_Q}}={\int }_0^1{\int }_0^1{\displaystyle \stackrel{e}{{\mathit{n}}^T}}c\mathrm{wsp}\left(\mathit{J}\right)d{\xi }_{3}d{\eta }_i$$E8

where ${\displaystyle \stackrel{e}{\mathit{n}}}$ and $c$ are the element shape functions vector and the scalar density of the surface electric charges.

2.3. Stationary piezoelectricity

The local and variational formulations of linear piezoelectricity combine our former considerations concerning the linear elasticity and linear dielectricity [13, 27]. This approach was repeated in [1]. The corresponding finite element formulation can be written in the form a coupled system of equations. The coupling is represented by the matrix ${\mathit{K}}_C$ in the following way

$$\begin{array}{l}{\mathit{K}}_M{\mathit{q}}^{q,\mathit{hp}}-{\mathit{K}}_C{\mathit{\phi }}^{\rho ,\mathit{h\pi }}={\mathit{F}}_V+{\mathit{F}}_S,\\ {\mathit{K}}_C^T{\mathit{q}}^{q,\mathit{hp}}+{\mathit{K}}_E{\mathit{\phi }}^{\rho ,\mathit{h\pi }}={\mathit{F}}_Q\end{array}$$E9

The coupling term can be called the global characteristic matrix of piezoelectricity, while the rest terms retain their previous meaning. The additional remark concerns special or simplified versions of the above equation. The inverse or direct piezoelectric problems can be considered here with the right-hand side terms equal to zero in the first and second equation, respectively. It is also worth mentioning that different $\mathit{pq}$ and $\mathit{\pi \rho }$ adaptive approximations of the vectorial displacement and scalar electric fields are proposed in (9), with the common $h$-approximation.

The global matrix of piezoelectricity introduced above can be obtained through the standard finite element summation procedure, where the following element contributions are employed

$${\displaystyle \stackrel{e}{{\mathit{k}}_C}}={\int }_0^1{\int }_0^1{\int }_0^{-{\xi }_2+1}{\displaystyle \stackrel{e}{{\mathit{B}}^T}}\mathit{C}{\displaystyle \stackrel{e}{\mathit{b}}}det\left(\mathit{J}\right)d{\xi }_{1}d{\xi }_{2}d{\xi }_3$$E10

with $\mathit{C}$ representing the piezoelectric (coupling) constants matrix.

The element contributions to the other terms of (9) are defined as before, that is, in accordance with (2)–(4) and (6)–(7). Note that the different shape functions matrices, ${\displaystyle \stackrel{e}{\mathit{N}}}$ and ${\displaystyle \stackrel{e}{\mathit{n}}}$, for the displacements and potential fields are employed here due to the different orders of approximations, $\mathit{pq}$ and $\mathit{\pi \rho }$, within both fields. Thanks to this, the corresponding adaptation processes within both fields can be performed independently.

3. The applied numerical techniques

3.1. Hierarchies of models

The presented elastic, dielectric and piezoelectric models are all based on the 3D-based approach, which results in the application of the three-dimensional or 3D-based degrees of freedom (dofs) only. The mechanical shell and transition models are also equipped with such dofs. This means that mid-surface degrees of freedom of the conventional shell and transition models are not applied. The so-called through-thickness dofs are employed instead. Also, some constraints are imposed on the three-dimensional displacements field of the shell and transition models so as to obtain the equivalence of the conventional and 3D-based descriptions. The related issues are presented in detail in the works [3, 5]. Analogously, in [7], the 3D-based hierarchy of dielectric models was proposed. It includes the three-dimensional and symmetric-thickness hierarchical models. Three-dimensional and 3D-based through-thickness dofs are employed in these models. In the latter work, also the 3D-based mechanical and dielectric models were combined, so as to obtain the 3D-based hierarchy of the piezoelectric models. This idea was also recalled in [1]. Note that all the presented 3D-based models, either elastic, dielectric or piezoelectric ones, can be treated as the 3D models polynomially constrained through the thickness.

The mechanical hierarchy $\mathrm{M}$ of the 3D or 3D-based elastic models $M$ reads:

$$M\in \mathrm{M},\mathrm{M}=\left\{3D,\mathit{MI},\mathit{RM},3D/\mathit{RM},\mathit{MI}/\mathit{RM}\right\}$$E11

with $3D$ denoting three-dimensional elasticity, $\mathit{MI}$ representing hierarchical shell models of higher order, $\mathit{RM}$ being the first-order shell model corresponding to Reissner theory of shells and $3D/\mathit{RM}$ and $\mathit{MI}/\mathit{RM}$ standing for the transition models of solid-to-shell or shell-to-shell character. The hierarchical shell and shell-to-shell models constitute two sub-hierarchies:

$$\begin{array}{c}\mathit{MI}=\left\{M2,M3,M4,\dots \right\},\\ \mathit{MI}/\mathit{RM}=\left\{M2/\mathit{RM},M3/\mathit{RM},M4/\mathit{RM},\dots \right\}\end{array}$$E12

where $I$ represents the order of the hierarchical model $\mathit{MI}$. This order is equivalent to the order of polynomial constraints defining the transverse displacement.

Subsequently, the hierarchy $\mathrm{E}$ of 3D-based dielectric models $E$ includes:

$$E\in \mathrm{E},\mathrm{E}=\left\{3D,\mathit{EJ}\right\}$$E13

where $3D$ represents three-dimensional theory of dielectricity, while $\mathit{EJ}$ denotes the 3D-based hierarchical models. The latter models constitute the following subhierarchy:

$$\mathit{EJ}=\left\{E1,E2,E3,\dots \right\}$$E14

with $J$ standing for the order of the hierarchical dielectric theory (the polynomial order of the through-thickness constraints of electric potential).

The 3D-based hierarchy $\mathrm{P}$ of piezoelectric models $P$ consists of the following component models:

$$P\in \mathrm{P},\mathrm{P}=\left\{\left(M,E\right):M\in M,E\in E\right\}$$E15

The hierarchy is composed of all combinations $\left(M,E\right)$ of the elastic models $M$ from (11) and (12) and dielectric models $E$ defined in (13) and (14), that is,

$$\begin{array}{c}\mathrm{P}=\{\left(3D,3D\right),\left(\mathit{MI},3D\right),\left(\mathit{RM},3D\right),\left(3D/\mathit{RM},3D\right),\left(\mathit{MI}/\mathit{RM},3D\right)\\ \left(3D,\mathit{EJ}\right),\left(\mathit{MI},\mathit{EJ}\right),\left(\mathit{RM},\mathit{EJ}\right),\left(3D/\mathit{RM},\mathit{EJ}\right),\left(\mathit{MI}/\mathit{RM},\mathit{EJ}\right)\}\end{array}$$E16

3.2. Hierarchical and constrained approximations

In the proposed approach, each of the 3D or 3D-based elastic, dielectric and piezoelectric models is approximated with the three-dimensional hierarchical shape functions. The functions applied in this work are originated from [8]. Their main feature is that they allow different orders of approximation on each of the element edges and sides. Such different orders are necessary for the local (element) $q$- and $p$-adaptivity. These different orders are obtained due to the shape function definition based on tensor products of the directional (longitudinal and transverse) shape functions of different orders. The specific form of the directional and three-dimensional functions for the case of the 3D solid and 3D-based hierarchical shell elements was presented in [3, 9]. The case of the solid-to-shell and shell-to-shell elements is addressed in [3, 12], while the first-order shell element shape functions are shown in [3, 10]. The analogous functions for the three-dimensional and hierarchical symmetric-thickness dielectric elements are given in [26]. In the case of the piezoelectric elements, the idea is to combine the elastic elements of various mechanical models with the dielectric elements. This idea is implemented in [1, 26]. Some details concerning shape functions of the component (elastic and dielectric) and combined (piezoelectric) elements are presented in the following paragraphs.

The displacement field of the elastic and piezoelectric elements is defined through the interpolation function $\mathit{u}=\mathit{u}\left(\mathit{\xi }\right)$ describing displacements $\mathit{u}={\left(u_1,u_2,u_3\right)}^T$ of any point $\mathit{\xi }$ of the normalized geometry of the element. This interpolant is a sum of four component functions:

$$\mathit{u}\left(\mathit{\xi }\right)={\mathit{u}}^1\left(\mathit{\xi }\right)+{\mathit{u}}^2\left(\mathit{\xi }\right)+{\mathit{u}}^3\left(\mathit{\xi }\right)+{\mathit{u}}^4\left(\mathit{\xi }\right)$$E17

The first component function ${\mathit{u}}^1\left(\mathit{\xi }\right)$ of the element vertices is defined as a product of the linear vertex node shape function matrix ${\mathit{N}}_v$ and the corresponding vector of nodal threesomes of directional dofs, that is,

$${\mathit{u}}^1\left(\mathit{\xi }\right)={\mathit{N}}_v\left(\mathit{\xi }\right){\mathit{q}}_v$$E18

where the mentioned vector is: ${\mathit{q}}_v={\left[\dots ,q_{1,i},q_{2,i},q_{3,i},\dots \right]}^T$, and where $i=1,2,\dots ,I_v$ with $I_v$ being the number of vertex nodes within the element.

The second component interpolant ${\mathit{u}}^2\left(\mathit{\xi }\right)$, corresponding to the element edges, is equal to the product of the higher-order shape function matrices, ${\mathit{N}}_h$ and ${\mathit{N}}_u$, of the element horizontal and vertical (upright) mid-edge nodes and the corresponding dof vectors:

$${\mathit{u}}^2\left(\mathit{\xi }\right)={\mathit{N}}_h\left(\mathit{\xi }\right){\mathit{q}}_h+{\mathit{N}}_u\left(\mathit{\xi }\right){\mathit{q}}_u$$E19

The mentioned vectors of degrees of freedom at the horizontal and vertical nodes are as follows: ${\mathit{q}}_h={\left[\dots ,q_{1,i,k},q_{2,i,k},q_{3,i,k},\dots \right]}^T$, ${\mathit{q}}_u={\left[\dots ,q_{1,j,l},q_{2,j,l},q_{3,j,l},\dots \right]}^T$, $i=1,2,\dots ,I_h$, $j=1,2,\dots ,I_u$, where $I_h$ and $I_u$ are the numbers of the horizontal and upright mid-edge nodes, while $k$ and $l$ represent numbers of dofs at these nodes.

The subsequent interpolant ${\mathit{u}}^3\left(\mathit{\xi }\right)$ corresponds the higher order mid-base and mid-side nodes of the element. The function is obtained through the multiplication of the shape function matrices, ${\mathit{N}}_b$ and ${\mathit{{\rm N}}}_s$, by the corresponding dofs vectors in accordance with

$${\mathit{u}}^3\left(\mathit{\xi }\right)={\mathit{N}}_b\left(\mathit{\xi }\right){\mathit{q}}_b+{\mathit{N}}_s\left(\mathit{\xi }\right){\mathit{q}}_s$$E20

where the vectors of nodal dofs are equal to: ${\mathit{q}}_b={\left[\dots ,q_{1,i,k},q_{2,i,k},q_{3,i,k},\dots \right]}^T$, ${\mathit{q}}_s={\left[\dots ,q_{1,j,l},q_{2,j,l},q_{3,i,l},\dots \right]}^T$ with $i=1,2,\dots ,I_b$ and $j=1,2,\dots ,I_s$. Here, $I_b$ and $I_s$ denote the numbers of the mid-base and mid-side nodes, while $k$ and $l$ are dof numbers at these nodes.

The last component interpolant ${\mathit{u}}^4\left(\mathit{\xi }\right)$, assigned to the element higher order middle node, is defined as a product of the shape function matrix ${\mathit{N}}_m$ and the corresponding dof vector:

$${\mathit{u}}^4\left(\mathit{\xi }\right)={\mathit{N}}_m\left(\mathit{\xi }\right){\mathit{q}}_m$$E21

where the dof vector is: ${\mathit{q}}_m={\left[\dots ,q_{1,k},q_{2,k},q_{3,k},\dots \right]}^T$, with $k$ standing for a dof number at this node.

With ${\displaystyle \stackrel{e}{\mathit{N}}}=\left[{\mathit{N}}_v\left(\mathit{\xi }\right),{\mathit{N}}_h\left(\mathit{\xi }\right),{\mathit{N}}_u\left(\mathit{\xi }\right),{\mathit{N}}_b\left(\mathit{\xi }\right),{\mathit{N}}_v\left(\mathit{\xi }\right),{\mathit{N}}_m\left(\mathit{\xi }\right)\right]$ and ${\displaystyle \stackrel{e}{\mathit{q}}}={\left[{\mathit{q}}_v,{\mathit{q}}_h,{\mathit{q}}_u,{\mathit{q}}_b,{\mathit{q}}_s,{\mathit{q}}_m\right]}^T$ Eq. (17) can be written in the alternative standard form

$$\mathit{u}\left(\mathit{\xi }\right)={\displaystyle \stackrel{e}{\mathit{N}}}\left(\mathit{\xi }\right){\displaystyle \stackrel{e}{\mathit{q}}}$$E22

The function $\varphi =\varphi \left(\mathit{\xi }\right)$ interpolating electric potential at any point $\mathit{\xi }$ of the normalized dielectric or piezoelectric element is also defined as a sum of four components

$$\varphi \left(\mathit{\xi }\right)={\varphi }^1\left(\mathit{\xi }\right)+{\varphi }^2\left(\mathit{\xi }\right)+{\varphi }^3\left(\mathit{\xi }\right)+{\varphi }^4\left(\mathit{\xi }\right)$$E23

The linear interpolant ${\varphi }^1\left(\mathit{\xi }\right)$ of the element vertices is equal to the product of the vector ${\mathit{n}}_v$ of shape functions for the vertices and the corresponding nodal vector of scalar dofs, that is,

$${\varphi }^1\left(\mathit{\xi }\right)={\mathit{n}}_v\left(\mathit{\xi }\right){\mathit{\phi }}_v$$E24

where ${\mathit{\phi }}_v={\left[\dots ,{\phi }_i,\dots \right]}^T$ and $i=1,2,\dots ,J_v$, with $J_v$ denoting the number of vertices within the element.

The higher order interpolant ${\varphi }^2\left(\mathit{\xi }\right)$, corresponding to horizontal and vertical (upright) mid-edge nodes, is equal to a product of the shape function vectors ${\mathit{n}}_h$ and ${\mathit{n}}_u$ and the corresponding vectors of nodal dofs:

$${\varphi }^2\left(\mathit{\xi }\right)={\mathit{n}}_h\left(\mathit{\xi }\right){\mathit{\phi }}_h+{\mathit{n}}_u\left(\mathit{\xi }\right){\mathit{\phi }}_u$$E25

with ${\mathit{\phi }}_h={\left[\dots ,{\phi }_{i,k},\dots \right]}^T$, ${\mathit{\phi }}_u={\left[\dots ,{\phi }_{i,l},\dots \right]}^T$, $i=1,2,,\dots ,J_h$, $j=1,2,\dots ,J_u$. Above, the numbers of the horizontal and upright mid-edge nodes are equal to $J_h$ and $J_u$, respectively, while $k$ and $l$ are numbers of the consecutive dofs at these nodes.

The next higher order interpolation function ${\varphi }^3\left(\mathit{\xi }\right)$, dealing with the mid-base and mid-side nodes, can be calculated with the multiplication of the shape function vectors ${\mathit{n}}_b$ and ${\mathit{n}}_s$ of these nodes and the respective vectors of nodal dofs:

$${\varphi }^3\left(\mathit{\xi }\right)={\mathit{n}}_b\left(\mathit{\xi }\right){\mathit{\phi }}_b+{\mathit{n}}_s\left(\mathit{\xi }\right){\mathit{\phi }}_s$$E26

while ${\mathit{\phi }}_b={\left[\dots ,{\phi }_{i,k},\dots \right]}^T$ and ${\mathit{\phi }}_s={\left[\dots ,{\phi }_{i,l},\dots \right]}^T$. Additionally, $i=1,2,\dots ,J_b$ and $j=1,2,\dots ,J_s$ with $J_b$ and $J_s$ standing for the numbers of the mid-base and mid-side nodes, respectively, and $k$, $l$ being dofs numbers at these nodes.

The last component function ${\varphi }^4\left(\mathit{\xi }\right)$, assigned for the middle node, needs multiplication of the shape function vector ${\mathit{n}}_m$ of the node and the corresponding vector of the nodal dofs

$${\varphi }^4\left(\mathit{\xi }\right)={\mathit{n}}_m\left(\mathit{\xi }\right){\mathit{\phi }}_m$$E27

where ${\mathit{\phi }}_m={\left[\dots ,{\phi }_k,\dots \right]}^T$, and $k$ represents the number of a hierarchical dof at the middle node.

Note that when ${\displaystyle \stackrel{e}{\mathit{n}}}=\left[{\mathit{n}}_v\left(\mathit{\xi }\right),{\mathit{n}}_h\left(\mathit{\xi }\right),{\mathit{n}}_u\left(\mathit{\xi }\right),{\mathit{n}}_b\left(\mathit{\xi }\right),{\mathit{n}}_s\left(\mathit{\xi }\right),{\mathit{n}}_m\left(\mathit{\xi }\right)\right]$ and ${\displaystyle \stackrel{e}{\mathit{\phi }}}={\left[{\mathit{\phi }}_v,{\mathit{\phi }}_h,{\mathit{\phi }}_u,{\mathit{\phi }}_b,{\mathit{\phi }}_s,{\mathit{\phi }}_m\right]}^T$, Eq. (23) can be written in the well-known general form

$$\varphi \left(\mathit{\xi }\right)={\displaystyle \stackrel{e}{\mathit{n}}}\left(\mathit{\xi }\right){\displaystyle \stackrel{e}{\mathit{\phi }}}$$E28

Here, we discuss on the constrained approximation. Such an approximation is necessary for $h$-adaptivity, which results in neighborhood of the element of different sizes, that is, the undivided elements $e$ of the initial mesh are adjacent to the divided elements $f$ of the $h$-adapted mesh. A further consequence of the different sizes is the constrained (or hanging) nodes of the smaller elements, which do not possess their counterparts in the neighboring bigger elements. In order to assure continuity of the field of displacements and the electric potential field between such elements, the constraining relations have to be introduced to the contributions of the smaller elements $f$ to Eqs. (1), (5) and (9), before the assemblage of the global matrices and vectors. The constraining relation for the case of displacements reads:

$${\displaystyle \stackrel{f}{\mathit{q}}}=\left[\begin{array}{c}{\displaystyle \stackrel{f}{{\mathit{q}}_s}}\\ {\displaystyle \stackrel{f}{{\mathit{q}}_u}}\end{array}\right]=\left[\begin{array}{cc}\mathit{I}& \mathbf{0}\\ \mathbf{0}& {\displaystyle \stackrel{\mathit{fe}}{{\mathit{C}}_q}}\end{array}\right]\left[\begin{array}{c}{\displaystyle \stackrel{f}{{\mathit{q}}_s}}\\ {\displaystyle \stackrel{e}{{\mathit{q}}_w}}\end{array}\right]$$E29

where the parts ${\displaystyle \stackrel{f}{{\mathit{q}}_s}}$ and ${\displaystyle \stackrel{f}{{\mathit{q}}_u}}$ of ${\displaystyle \stackrel{f}{\mathit{q}}}$ include the unconstrained and constrained dofs of the smaller element $f$, ${\displaystyle \stackrel{e}{{\mathit{q}}_w}}$ contains displacements of the constraining nodes of the bigger neighbor $e$, while ${\displaystyle \stackrel{\mathit{fe}}{{\mathit{C}}_q}}$ and $\mathit{I}$ represent the constraint coefficient matrix and the unity matrix. In the case of the electric potential, the analogous relation reads:

$${\displaystyle \stackrel{f}{\mathit{\phi }}}=\left[\begin{array}{c}{\displaystyle \stackrel{f}{{\mathit{\phi }}_s}}\\ {\displaystyle \stackrel{f}{{\mathit{\phi }}_u}}\end{array}\right]=\left[\begin{array}{cc}\mathit{I}& \mathbf{0}\\ \mathbf{0}& {\displaystyle \stackrel{\mathit{fe}}{{\mathit{C}}_{\phi }}}\end{array}\right]\left[\begin{array}{c}{\displaystyle \stackrel{f}{{\mathit{\phi }}_s}}\\ {\displaystyle \stackrel{e}{{\mathit{\phi }}_w}}\end{array}\right]$$E30

The general rules for the constrained approximation are presented in [8]. These rules are applied in [3, 14, 26] where the methods of obtainment of the constraint coefficients for the two-dimensional and three-dimensional cases are described.

3.3. Error estimation

The equilibrated residual method of error estimation [14, 15, 16, 19] applied to solid mechanics is based on the solution of the approximated local (element) problems of mechanical equilibrium. The corresponding equilibrium condition, written in the language of finite elements, takes the form:

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

with ${{\displaystyle \stackrel{e}{\mathit{k}}}}_M$ standing for the element stiffness matrix, ${{\displaystyle \stackrel{e}{\mathit{f}}}}_V$ and ${{\displaystyle \stackrel{e}{\mathit{f}}}}_S$ denoting the nodal mass and surface forces vectors, and ${{\displaystyle \stackrel{e}{\mathit{f}}}}_R$ representing the nodal forces vector due to the equilibrated interelement stress loadings. In (31), ${{\displaystyle \stackrel{e}{\mathit{q}}}}^{Q,\mathit{HP}}$ is the solution displacements vector in the approximated local problem. Note that the discretization parameters $H$, $P$ and $Q$ (the element size, and the longitudinal and transverse approximation orders in the local problems) can be different to their global counterparts $h$, $p$ and $q$. The element solutions from the above relation give the global error estimate, which upper-bounds the true error [15, 16, 19]. In the presented approach, the above relation is applied to both the approximation and total errors estimation. However, different values of the discretization parameters are applied in both cases, that is, $H=h$, $P=p+1$, $Q=q$ and $H=h$, $P=p+1$, $Q=q+1$, respectively. The modeling error is calculated as the difference of the previous two errors.

As demonstrated in the work [20], application of the equilibrated residual methods to dielectric problems needs solution of the local (element) electric equilibrium problems. Such equilibrium, expressed in the language of finite elements, can be written in the following way:

$${{\displaystyle \stackrel{e}{\mathit{k}}}}_E{{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\mathrm{P},H\mathrm{\Pi }}={{\displaystyle \stackrel{e}{\mathit{f}}}}_Q+{{\displaystyle \stackrel{e}{\mathit{f}}}}_H$$E32

In the above-mentioned equation, the term ${{\displaystyle \stackrel{e}{\mathit{k}}}}_E$ stands for the element characteristic matrix of dielectricity, ${{\displaystyle \stackrel{e}{\mathit{f}}}}_Q$ is the element nodal vector due to surface charges, while ${{\displaystyle \stackrel{e}{\mathit{f}}}}_H$ represents nodal forces due to the equilibrated interelement charges. The vector ${{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\mathrm{P},H\mathrm{\Pi }}$ is the solution vector of electric potential in the local problem. The solution is approximate and the corresponding discretization parameters $H$, $\mathrm{\Pi }$, $\mathrm{P}$ (the element size and longitudinal and transverse approximation orders) can be different to their global counterparts $h$, $\pi$, $\rho$. In the cases of the approximation and total errors, the values of $H=h$, $\mathrm{\Pi }=\pi +1$, $\mathrm{P}=\rho$ and $H=h$, $\mathrm{\Pi }=\pi +1$, $\mathrm{P}=\rho +1$ are applied, respectively. Again, the collection of the local solutions obtained by means of (32) leads to the error estimate, which upper-bounds the true total and true approximation errors, and the modeling error is defined as the difference of these two errors.

Generalization of the equilibrated residual methods onto piezoelectric problems was proposed and developed in [1, 20, 26]. In accordance with this proposition, the above mechanical and electrical local equilibria Eqs. (31) and (32) have to be replaced by the coupled equations describing electromechanical equilibrium. Such equations take the following finite element form:

$$\begin{array}{c}{{\displaystyle \stackrel{e}{\mathit{k}}}}_M{{\displaystyle \stackrel{e}{\mathit{q}}}}^{Q,\mathit{HP}}-{{\displaystyle \stackrel{e}{\mathit{k}}}}_C{{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\mathrm{P},H\mathrm{\Pi }}={{\displaystyle \stackrel{e}{\mathit{f}}}}_V+{{\displaystyle \stackrel{e}{\mathit{f}}}}_S+{{\displaystyle \stackrel{e}{\mathit{f}}}}_R\\ {{\displaystyle \stackrel{e}{\mathit{k}}}}_C^T{{\displaystyle \stackrel{e}{\mathit{q}}}}^{Q,\mathit{HP}}+{{\displaystyle \stackrel{e}{\mathit{k}}}}_E{{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\mathrm{P},H\mathrm{\Pi }}={{\displaystyle \stackrel{e}{\mathit{f}}}}_Q+{{\displaystyle \stackrel{e}{\mathit{f}}}}_H\end{array}$$E33

The coupled local solutions ${{\displaystyle \stackrel{e}{\mathit{q}}}}^{Q,\mathit{HP}}$ and ${{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\mathrm{P},H\mathrm{\Pi }}$ of the above set have one disadvantage. It lies in the lack of the upper-bound property of the total, approximation and modeling errors by the residual-based global estimators obtained from the local solutions of (33).

In the works [1, 26], the decoupled version of the above set was also proposed as a simpler alternative:

$$\begin{array}{c}{{\displaystyle \stackrel{e}{\mathit{k}}}}_M{{\displaystyle \stackrel{e}{\mathit{q}}}}^{Q,\mathit{HP}}={{\displaystyle \stackrel{e}{\mathit{k}}}}_C{{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\rho ,h\mathrm{\pi }}+{{\displaystyle \stackrel{e}{\mathit{f}}}}_V+{{\displaystyle \stackrel{e}{\mathit{f}}}}_S+{{\displaystyle \stackrel{e}{\mathit{f}}}}_R\\ {{\displaystyle \stackrel{e}{\mathit{k}}}}_E{{\displaystyle \stackrel{e}{\mathit{\phi }}}}^{\mathrm{P},H\mathrm{\Pi }}=-{{\displaystyle \stackrel{e}{\mathit{k}}}}_C^T{{\displaystyle \stackrel{e}{\mathit{q}}}}^{q,\mathit{hp}}+{{\displaystyle \stackrel{e}{\mathit{f}}}}_Q+{{\displaystyle \stackrel{e}{\mathit{f}}}}_H\end{array}$$E34

In Eqs. (31) and (32) and in the above two sets of equations, the vectors of the equilibrated nodal forces and charge are defined in accordance with:

$${\displaystyle \stackrel{e}{{\mathit{f}}_R}}={\int }_{S_e}{\displaystyle \stackrel{e}{{\mathit{N}}^T}}〈{\mathit{r}}_e\left({\mathit{u}}^{\mathit{hpq}}\right)〉{\mathit{dS}}_e=\sum _f{\int }_{S_{\mathit{ef}}}{\displaystyle \stackrel{e}{{\mathit{N}}^T}}〈{\mathit{r}}_e\left({\mathit{u}}^{\mathit{hpq}}\right)〉{\mathit{dS}}_{\mathit{ef}}$$E35

and

$${\displaystyle \stackrel{e}{{\mathit{f}}_H}}={\int }_{S_e\backslash S}{\displaystyle \stackrel{e}{{\mathit{n}}^T}}〈h_e\left({\varphi }^{\mathit{h\pi \rho }}\right)〉{\mathit{dS}}_e=\sum _f{\int }_{\mathit{Sef}}{\displaystyle \stackrel{e}{{\mathit{n}}^T}}〈h_e\left({\varphi }^{\mathit{h\pi \rho }}\right)〉{\mathit{dS}}_{\mathit{ef}}$$E36

where S_e and S are the element and body surfaces, while $〈{\mathit{r}}_e\left({\mathit{u}}^{\mathit{hpq}}\right)〉$ and $〈h_e\left({\varphi }^{\mathit{h\pi \rho }}\right)〉$ denote the equilibrated interelement (between element $e$ and any of its neighbors $f$) stress vector and the equilibrated interelement charge density acting on the common sides $S_{\mathit{ef}}$ between the elements. Some more details on how to calculate the equilibrated stresses and charge can be found in [3, 19, 26].

3.4. Adaptive strategy

The adaptive strategy applied in this chapter takes advantage of the Texas three-step strategy [21]. The original strategy is assigned for structures of simple geometry (a single geometrical part of one type) and simple physical description (one model). Such a strategy is designed for $\mathit{hp}$-adaptivity and consists of three steps: initial, intermediate and final ones, where three subsequent global problem solutions are obtained. In this strategy, the solution on the initial mesh is followed by the error estimation and error-controlled $h$-adaptation (element refinement). In this way, the intermediate mesh is formed. The solution on this mesh is then followed by the error estimation and the error-controlled $p$-adaptivity (element approximation order enrichment). The problem solution on this mesh is followed by the error estimation again.

The original strategy was extended in three ways [3, 26]. First, structures of complex geometry and complex physical description can be analyzed in the presented approach. Second, the $h$-step of the adaptation is enriched with the model adaptivity, while the $p$-step of the strategy is completed with $q$-adaptivity (enrichment of the element transverse approximation order) [19] in the way different to the earlier proposition of [4]. Third, the strategy is enhanced through the addition of one more adaptation step called the modification one. In this step, the initial mesh is modified so as to remove the undesired numerical phenomena, such as the improper solution limit, numerical locking or boundary layers [28]. The extended strategy can be applied to elastic [3, 19], dielectric [1] and piezoelectric [1, 26] problems. Some difficulties in the application of the strategy to the cases of piezoelectricity are described in [1]. The analogous comments for the case of elasticity can be found in [29, 30].

Description of the error-controlled adaptivity for elastic, dielectric and piezoelectric structures is started with the $h$-adaptivity within the mechanical field (hence index $M$). In accordance with [19, 21], the new number of elements $n_{I_M}$ in the intermediate mesh (denoted with index $I$), which replace an element of the initial mesh, is equal to:

$$n_{I_M}^{2{\mu }_{0_M}/d+1}=\frac{{\eta }_{0_M}^2{E}_I}{{\gamma }_{a,I}^2‖{}_{}^{\mathrm{}}{\mathit{u}}_0‖{}_U^2}$$E37

where ${\eta }_{0_M}$ is the estimated value of the approximation error from the initial mesh, ${\mu }_{0_M}$ represents the known or assumed $h$-convergence rate, $d$ denotes dimensionality (equal to 2 or 3) of the adapted geometrical part. Additionally, $E_I$ is the total number of elements in the intermediate mesh, $‖{\mathit{u}}_0‖{}_U^2$ is the strain energy norm of the solution from the initial mesh, while the coefficient ${\gamma }_{a,I}$ determines the expected relative value of the global approximation error within the intermediate mesh.

In the case of the mechanical field, local (element) $p$-enrichments are controlled with the longitudinal (or overall) element approximation order $p_T$ of the final (or target) mesh (marked with the index $T$). The longitudinal order corresponds to thin-walled structures, while the overall order to three-dimensional bodies. The following common formula [19, 21] determines this parameter:

$$\begin{array}{l}{p}_T^{2{\nu }_{0M}}=\frac{p_0^{2{\nu }_{0_M}}{\eta }_{I_M}^2{E}_I}{{\gamma }_{a,T}^2‖{}_{}^{\mathrm{}}{\mathit{u}}_0‖{}_U^2}\\ \end{array}$$E38

with ${\eta }_{I_M}$ representing the estimated value of the approximation error from the intermediate mesh, ${\nu }_{0_M}$ being the given $p$-convergence rate, $p_0$ denoting the longitudinal approximation order from the initial mesh and ${\gamma }_{a,T}$ standing for the expected relative value of the global approximation error in the target mesh.

In the case of $q$-adaptivity within the mechanical field of the thin-walled structures, the target values of the element transverse approximation orders $q_T$ can be defined in accordance with [3, 19]:

$$q_T=q_0-\frac{1}{2}{log}_{t/2l}\frac{{\theta }_{I_M}^2{E}_I}{{\gamma }_{m,T}^2‖{}_{}^{\mathrm{}}{\mathit{u}}_0‖{}_U^2}$$E39

Above, $q_0$ is the element transverse order from the initial mesh, $t$ and $2l$ represent the thickness and length of the thin-walled part of the structure, ${\theta }_{I_M}$ is the estimated value of the modeling error from the intermediate mesh and ${\gamma }_{m,T}$ denotes the expected relative value of the modeling error in the target mesh.

The $h$-adaptivity within the electric field needs determination of the new number (denoted with index $E$) of elements, $n_{I_E}$, in the intermediate mesh (marked with index $I$ again). Such a number has to be determined for each element of the initial mesh. This number is equal to:

$$n_{I_E}^{2{\mu }_{0_E}/d+1}=\frac{{\eta }_{0_E}^2{E}_I}{{\gamma }_{a,I}^2‖{}_{}^{\mathrm{}}{\varphi }_0‖{}_W^2}$$E40

Above, the quantity ${\eta }_{0_E}$ stands for the estimated value of the approximation error in the initial mesh, the exponent ${\mu }_{0_E}$ is the assumed $h$-convergence rate and the norm $‖{\varphi }_0‖{}_W^2$ represents the electrostatic energy corresponding to the initial mesh.

Note that in the case of piezoelectricity, for each element, the final subdivision is determined with

$$n_I=max\left(n_{I_M},n_{I_E}\right)$$E41

as the common mesh division is applied for the mechanical and electric fields (see [1]).

The value of the longitudinal (or overall) approximation order ${\pi }_T$ within the electric field of an element of the target (final) mesh can be calculated from [26]:

$${\pi }_T^{2{\nu }_{0_M}}=\frac{{\pi }_{0_E}^{2{\nu }_{0_E}}{\eta }_{I_E}^2{E}_I}{{\gamma }_{a,T}^2‖{}_{}^{\mathrm{}}{\varphi }_0‖{}_W^2}$$E42

where ${\eta }_{I_E}$ is the estimated approximation error value from the intermediate mesh, ${\nu }_{0_E}$ is the assumed problem-dependant $p$-convergence rate and ${\pi }_0$ is the longitudinal approximation order within the initial mesh.

Finally, in the case of $q$-adaptivity of thin dielectric or piezoelectric body, the element transverse order of approximation of the electric potential field in the final mesh has to be determined. The following formula can be proposed for this purpose [26]:

$${\rho }_T={\rho }_0-\frac{1}{2}{log}_{t/2l}\frac{{\theta }_{I_E}^2{E}_I}{{\gamma }_{m,t}^2‖{\varphi }_0‖{}_W^2}$$E43

where ${\rho }_0$ is the element transverse approximation order applied in the initial mesh, $t$ and $2l$ are the transverse and longitudinal dimensions of a thin member (body or part) and ${\theta }_{I_E}$ is the modeling error estimated value of a finite element of the intermediate mesh.

4. Numerical examples

In this section, some comparative examples for problems of elasticity, dielectricity and piezoelectricity are presented. Attention is focused on the comparison of effectiveness of three main algorithms applied in our generalizing approach to adaptive modeling and analysis of the problems of three mentioned classes. The tested numerical procedures include hierarchical approximations, error estimation and error-controlled adaptivity. In the first case of hierarchical approximations, convergence curves for the three classes of problems are compared. In the second case of error estimation with the equilibrated residual method, effectivities of the global estimators in three problems are presented. In the third case of adaptive procedure, effectivity of the adaptation is checked through the comparison of the adaptive convergence curves. Also, the ability to reach the assumed admissible error in problems of three types is assessed.

4.1. Model structures

Here, the same domain geometry is considered in the problems of elasticity, dielectricity and piezoelectricity. Its square longitudinal dimensions are equal to $2l=3.1415\cdot {10}^{-2}$ m, while its thickness equals $t=0.01\cdot {10}^{-2}$ m. The domain thickness may change if necessary. This domain can represent a plate structure in the mechanical case, a thin dielectric in the electric case and a thin piezoelectric structure in the electromechanical case.

So as to be able to compare three different physical problems, physical properties of the materials and the external load and charge are assumed such that the induced mechanical and electric potential energies are of the same order. The isotropic mechanical properties of the plate structure and thin piezoelectric correspond to a typical piezoelectric material and are taken from [27]. Young’s modulus is assumed to be equal to $E=0.5\cdot {10}^{11}$ $\mathrm{N}/{\mathrm{m}}^2$, while the applied Poisson’s ratio is $\nu =0.294$. The isotropic dielectric properties of the dielectric and piezoelectric are characterized with the permittivity equal to $\delta =0.1593\cdot {10}^{-7}$ F/m. The nonzero anisotropic piezoelectric constants are $d_{13}=d_{23}=-0.15\cdot {10}^{-9}$ C/N, $d_{33}=0.3\cdot {10}^{-9}$ C/N and $d_{52}=d_{61}=0.5\cdot {10}^{-9}$ C/N (compare [27] again). The surface load of magnitude $p=0.4\cdot {10}^6$ $\mathrm{N}/{\mathrm{m}}^2$ is applied to the upper surface of the elastic and piezoelectric structures. Furthermore, the surface electric charge of density $c=0.2\cdot {10}^{-1}$ $\mathrm{C}/{\mathrm{m}}^2$ is applied to the upper surface of the dielectric and piezoelectric domains.

The kinematic boundary conditions within the mechanical field of displacements of the elastic and piezoelectric structures assume all edges (lateral sides) clamped—no displacements on these edges are present. In the case of the electric field of potential within the dielectric and piezoelectric domains, grounding is assumed around the domain (zero potential on the lateral sides).

In the next sections, only symmetric quarters of the structures are shown due to the symmetry of the applied geometry, load and charge distributions and boundary conditions.

4.2. Convergence of hierarchical approximations

Figures 1 and 2 illustrate distributions of the effective value ($\mathit{se}f$) of a stress tensor (effective stress) and the electric displacement vector magnitude ($\mathit{dm}$) for the purely mechanical and electric problems, while Figures 3 and 4 present the same quantities corresponding to the electromechanical problem. The displayed values are obtained due to solution of the corresponding global problems, either (1) or (5) or (9). Comparing Figure 1 with Figure 3 as well as Figures 2 and 4, one can notice that the stresses in the piezoelectric case are changed with respect to the purely elastic case due to the presence of the electromechanical coupling. Similarly, the electric displacements of the piezoelectric case look different to those of the purely dielectric example due to the influence of the coupled mechanical displacements field.

Figure 1.

Effective stress in the purely mechanical case.

Figure 2.

Magnitude of electric displacement in the purely electric case.

Figure 3.

Effective stress in the coupled problem.

Figure 4.

Magnitude of electric displacement in the coupled problem.

The $h$- and $p$-convergence curves for the elastostatic case are presented in Figures 5 and 6. The electrostatic problem solution $h$- and $\pi$-convergence curves are displayed in Figures 7 and 8. The solution convergence curves for the stationary piezoelectricity are shown in Figures 9 and 10. In these figures, the absolute values of the approximation errors are plotted versus the number $N$ of the applied degrees of freedom. The number of degrees of freedom changes due to the increase in either the number of subdivisions $m=l/h$ (the case of $h$-convergence) or the longitudinal approximation order $p$ or $\pi$ (the case of $p$-convergence). The applied values of the discretization parameters of the uniform meshes are $m=1,2,\dots ,8$, $p=\pi =1,2,\dots ,6$ while the transverse orders of approximation are kept constant and equal to $q\equiv I=2$ or/and $\rho \equiv J=2$, with $\pi \equiv \mathit{pi}$ and $\rho \equiv \mathit{rho}$. The error is calculated as a square root of the difference of the potential energy of the numerical solution and the exact value of this energy in three cases: mechanical, electric or electromechanical. As the exact values of the solutions to three problems are not known, these values are replaced by the best numerical ones obtained from the meshes of $p=\pi =9$, $h=9$, $q=\rho =2$.

Figure 5.

$h$-Convergence in the purely elastic problem.

Figure 6.

$p$-Convergence in the purely elastic problem.

Figure 7.

$h$-Convergence in the purely dielectric problem.

Figure 8.

$p$-Convergence in the purely dielectric problem.

Figure 9.

$h$-Convergence in the coupled piezoelectric problem.

Figure 10.

$p$-Convergence ($\pi =p$) in the coupled problem.

The following findings can be formulated based on the analysis and comparison of the drawings. The convergence curves for the purely mechanical problem consist of three parts. The first part, almost horizontal and flat, corresponds to the presence of the numerical locking. The second part of the highest slope corresponds to the so-called asymptotic convergence. The third part of worse convergence is affected by the influence of the boundary layer. Such a picture of convergence is typical for elastic problems and displacement finite element formulation (compare [5]). In the case of the pure dielectricity, the curves consists of two parts only—the second and third ones. No locking is observed for this problem. In the case of the coupled problem of piezoelectricity, three parts of the curves appear again. As far as the asymptotic convergence range is concerned, it should be noticed that in the purely electric case, the convergence is much higher (slopes are more steep) than in the purely mechanical problem. Additionally, the boundary layer effect in the dielectric problem is less severe than in the mechanical one—the slopes of the third parts of the curves are higher in the former case.

As it comes to the piezoelectricity, the following observations can be seen. The first parts of the convergence curves in the case of the coupled problem are not flat but are bent, that is, the monotonic character of these parts is not retained. This observation reflects the fact of change of the sign of the potential energy, which is composed of mechanical, electric and coupling contributions of different signs. The solution convergence curves of the piezoelectric problem, in the asymptotic and boundary-layer ranges, lie just between the corresponding curves of the pure problems, with a tendency to be closer to the purely mechanical case.

4.3. Effectivity of error estimation

In this section, results concerning application of the equilibrated residual method to the analogous model problems of elastostatics, electrostatics and stationary piezoelectricity are presented. The data of the problems are taken from Section 4.1. The thickness of the analyzed domains is now equal to $t=0.15\cdot {10}^{-2}$ m. The surface traction equals now $p=4.0\cdot {10}^6$ $\mathrm{N}/{\mathrm{m}}^2$. This value has changed due to the thickness change so as to assure the same order of the electric and mechanical potential energies, as well as the same order of the electric and mechanical parts of this energy, in the tests concerning three mentioned problems of elasticity, dielectricity and piezoelectricity. Presentation of the results starts with the purely elastic case. In Figure 11, the chosen example of the estimated total error distribution is presented for the uniform mesh discretization parameters $m=3$, $p=4$ and $q=2$. The level of the relative estimated total errors in elements (denoted as $\left(M\right)\mathit{nt}$) for the mechanical problem can be seen in the figure, as well as the average global value of the error estimator, marked as avr. The estimated local errors represent the square root of the difference between the mechanical potential energies of the numerical global solution under consideration, determined by (1), and the solution obtained from the local problems (31). In the case of a local error indicator, these energies are limited to a single element, while in the case of the global error estimator, the energies of all elements are taken into account. Figure 12 presents effectivity indices of the global error estimators as a function of the longitudinal approximation order $p$. The effectivity indices are calculated as ratios of the global estimators by the global values of the true total, approximation and modeling errors. The global values of the errors are equal to the square root of the difference of the appropriate energies. As the exact values of the solutions, necessary for the exact energy values, are not known, they are replaced with the numerical values obtained by using (1), from the finest and richest meshes possible, that is, with $m=9$, $p=9$ and $q=2$ or $q=6$ in the cases of the approximation and total errors calculations, respectively.

Figure 11.

Estimated total errors in the purely mechanical case.

Figure 12.

Effectivities of the estimators in the purely mechanical case.

Figures 13 and 14 display the analogous results for the electrostatic case. The first figure shows the exemplary distribution of the locally estimated relative total errors (denoted as $\left(E\right)\mathit{nt}$), and also the average error estimate (marked with avr), for $m=3$, $\pi \equiv \mathit{pi}=4$, $\rho \equiv \mathit{rho}=2$. The global solution to the problem (5) and the solutions to the local problems (32) are employed for the determination of two electric potential energies. These energies correspond to the numerical solution and the estimate of the exact solution. As far as the second figure is concerned, it illustrates the change of the effectivity indices of estimation of the total, approximation and modeling errors as a function of the longitudinal order of approximation $\pi$. In order to obtain the necessary exact values of the energies, the global problem (5) was applied again, with $m=9$, $\pi =9$ and $\rho =2$ or $\rho =6$, for the calculation cases of the approximation and total errors, respectively.

Figure 13.

Estimated total errors in the purely electric case.

Figure 14.

Effectivities of the estimators in the purely electric case.

The analogous estimated error distributions and the analogous plots of the effectivity indices versus $p=\pi$ for the piezoelectricity case are shown in Figures 15–18. Here, the numerical and estimated values of the mechanical and electric parts of potential energies, necessary for the exemplary local and average estimated error values determination, are obtained from (9), with $m=3$, $p=\pi =4$, $q=\rho =2$ and (33), respectively. The exact energy values necessary for effectivity calculations are obtained through the global numerical approximation (9), with $m=9$, $p=\pi =9$ and $q=\rho =2$ or $q=\rho =6$ and $\pi \equiv \mathit{pi}$, $\rho \equiv \mathit{rho}$.

Figure 15.

Mechanical parts of the estimated total errors (piezoelectricity).

Comparing Figures 11, 13, 15 and 17, one can notice that the estimated total error level for the cases of elasticity and dielectricity are not the same for the corresponding model problems, as the latter problem produces the lower local and global error estimates. The average relative values of the total error estimates for these two cases are equal to 0.128 and 0.013, respectively. In the case of the piezoelectricity, the estimated local and global error values are higher than in the previous two cases. The mechanical and electric parts of the average total error estimate are equal to 0.221 and 0.018, respectively. As far as the effectivities for three model problems from Figures 12, 14, 16, 18 are concerned, one can see that the pure problems deliver very similar effectivities—in both problems close to the desired values of 1. The values are almost everywhere larger than 1, and this result is consistent with the upper bound property of the estimation present in the pure problems (compare [19]). In the case of the coupled piezoelectric problem, global effectivities are much greater than 1, that is, the estimated global error values are much overestimated. This suggests that some additional techniques should be applied in the error estimation by the equilibrated residual method if one wants the effectivities to become closer to 1. A remedy can be the application of the higher order equilibration [14] in the piezoelectricity problems, instead of the linear equilibration used in this work and usually applied to pure problems.

Figure 16.

Effectivities of estimators’ mechanical parts (piezoelectricity).

Figure 17.

Electrical parts of the estimated total errors (piezoelectricity).

Figure 18.

Effectivities of estimators’ electrical parts (piezoelectricity).

4.4. Convergence of $\mathit{\text{hp}}$-adaptivity

In Figures 19–24, the results illustrating both $\mathit{hp}$-adapted meshes and convergence of the corresponding adaptation processes are presented for three model problems of elasticity, dielectricity and piezoelectricity. In these problems, the original data are recalled from Section 4.1. The thickness of the domains is equal to $t=0.15\cdot {10}^{-2}$ m, as in the previous test. In this way, the influence of the error estimation, performed in the previous subsection, on effectiveness of the mesh adaptation presented here can be assessed. Also, it is worth noticing that for the applied thickness value, the locking and boundary layer phenomena do not influence convergence very much. Note also that the presented adaptation is controlled with the estimated values of the element approximation errors within the displacements and electric potential fields by means of the formulas (37), (38) or/and (40), (42).

Figure 19.

Final $\mathit{hp}$-adapted mesh in the purely mechanical case.

Figure 20.

$\mathit{hp}$-Adaptive convergence in the purely mechanical case.

Figure 21.

Final $\mathit{h\pi }$-adapted mesh in the purely electric case.

Figure 22.

$\mathit{h\pi }$-Adaptive convergence in the purely electric case.

Figure 23.

Final $\mathit{hp}$-adapted mesh in the coupled problem.

Figure 24.

$\mathit{hp}$-Adaptive convergence in the coupled problem.

In the first two figures, the purely elastic case is presented. Figure 19 presents the mesh obtained in the three-step adaptive process. Both the changes of the mesh density and the element changes of the longitudinal order of approximation can be seen in Figure 19. The initial mesh (not displayed) corresponds to the discretization parameters $m=3$, $p=4$ and $q=2$. In this mesh, the longitudinal approximation order equal to $p=4$ is applied in order to remove the influence of the locking phenomenon present in purely mechanical problems. The target approximation error for the final mesh is assumed to be 0.1 with the ratio of the errors from the intermediate and final meshes equal to 2. The next Figure 20 displays the $\mathit{hp}$-convergence of the solution during the adaptation process. In the error calculations, the exact value of the energy is replaced with the value obtained for the best numerical discretization of the second-order ($q=2$) hierarchical shell model, where $m=9$ and $p=9$. The final true error value corresponds to the lowest (third) point of the convergence curve. This value can be compared with the horizontal dotted line corresponding to the admissible error level.

Figures 21 and 22 present the similar results for the purely electric case. The first of them displays the final mesh resulting from the three-step adaptation. The only difference within the applied discretization parameters is $\pi =2$ within the electric potential field of the initial mesh. This value replaces $p=4$ in the displacements field of the previous example. The assumption of $\pi =2$ results from lower error level within the former field and the lack of the numerical locking within dielectric problems. Subsequently, the second figure illustrates the $\mathit{h\pi }$-convergence curve of the adapted solution. The curve can be compared with the admissible error level again. The exact solution, necessary for the error calculations, is approximated by the numerical solution corresponding to $m=9$, $\pi =9$ and $\rho =2$.

The next two couples, Figures 23 and 24, as well as Figures 25 and 26, present exactly the same results, that is, final meshes and adaptive convergence curves for the displacements and electric potential fields, respectively, in the case of the coupled problem of piezoelectricity. For both fields, exactly the same initial mesh and error control parameters are applied as for the pure problems of elasticity and dielectricity. The presence of locking in piezoelectric problems is taken into account, hence $p=4$ and $\pi =2$ are set within the initial meshes of the corresponding fields. In the case of the adaptive convergence curves, the admissible error levels are marked with the dotted lines again. The numerical approximation of the exact solution entering the error calculations is obtained from the discretization based on $m=9$, $p=\pi =9$, $q=\rho =2$.

Figure 25.

Final $\mathit{h\pi }$-adapted mesh in the coupled problem.

Figure 26.

$\mathit{h\pi }$-Adaptive convergence in the coupled problem.

Discussion of the obtained numerical results can be concluded in the following way. In the case of the pure elasticity, the admissible error value is reached in three steps. The estimated average approximation error relative value for the initial mesh is 0.084. A relatively rare $h$-mesh is generated as the error is relatively low in the initial mesh and overestimation for $p=4$ does not occur (compare Figure 12). Then the mesh is well $p$-enriched.

In the case of the pure dielectricity, the relatively fine $h$-mesh is produced because of the error overestimation for $\pi =2$ (see Figure 14) present in the initial mesh. The estimated average approximation error relative value for the initial mesh is moderate and equal to 0.124. As a result, the $p$-enrichment is reduced (barely visible). The admissible error is reached, however.

In the case of the piezoelectricity, the changes in the common $h$-mesh come from the relatively large errors of the displacements field in the initial mesh. The corresponding average error value for the displacements field is equal to 0.227, while for the electric potential field, it equals 0.127. The displacements field errors in elements are larger than those of the electric potential field. As a result, the common $h$-mesh is too fine for the electric potential. This effect is enforced by the overestimation for $\pi =2$ (Figure 18). The following $\pi$-enrichment is weak. However, the admissible error level is reached for the potential field. This demonstrates that the idea of the common $h$-mesh for both fields works well also in the case of the field of lower errors. As far as the displacement field is concerned, one can notice that the error overestimation for $p=4$ (see Figure 16) in the common $h$-mesh produces too rich $p$-mesh for displacements and thus the corresponding admissible error value is far exceeded. The obtained solution for displacements is better than expected.

5. Conclusions

In this chapter, the algorithms of hierarchical modeling, hierarchical approximations, error estimation and adaptivity control, so far utilized successfully for elastic problems, are applied to the problems of dielectricity and piezoelectricity. To the best knowledge, no examples of such application had been reported in the existing literature.

This chapter shows how to assess effectivity of the algorithms of hierarchical approximations, equilibrated residual method of error estimation and three-step adaptive procedure, originally applied to elasticity and possible also in dielectricity and the coupled piezoelectric problems.

The observations from the tests concerning hierarchical $\mathit{hp}$-approximations allow for the following generalizations. The applied hierarchical approximations can be effective in modeling and simulation of all three classes of problems. The $h$- and $p$-convergence rate is the highest for the purely dielectric problems and the lowest for the analogous purely elastic ones. The convergence of the piezoelectric problems is located between the convergences of the corresponding pure problems. In the case of the pure problems, the convergence curves are monotonic, while in the case of the coupled problems of piezoelectricity, the loss of monotonicity can happen because of the sign change of the electromechanical potential energy.

The general conclusions concerning the applied error estimation method are as follows. In the pure problem of elasticity, the effectivities of the modeling, approximation and total errors are all close to 1.0, except for the cases of poor discretization ($p=1,2$). In the pure dielectric problems, the values of the effectivities are between 0.9 and 1.0. For the cases of poor discretization, these values are close to 2.0. In the case of piezoelectricity, the indices are close to their values from the pure problems for low values of the approximation order, that is, for $p\le 3$ or/and $\pi \le 3$. For higher values of the approximation order, the corresponding effectivities are much higher than 1.0 for the approximation error. A bit smaller overestimation can be seen in the case of the total error. It can be concluded that for rich discretizations, some additional techniques are necessary in the case of piezoelectricity, so as to enforce values of the effectivities closer to 1.0.

Generalizations related to the applied adaptivity control algorithms can be formulated in the following way. For the analogous mechanical, electric and electromechanical problems, the three-step adaptive strategy leads to similar convergence results as it comes to the final mesh true error level, even though the $\mathit{hp}$-path in each of three cases can be different. The final error values are usually smaller or close to the admissible error level. In the case of piezoelectricity, too fine or/and too rich meshes may be generated, if overestimation, coming from the error estimation procedure, occurs.

The results concerning the effectivity of the application of the mentioned three algorithms to the dielectric and piezoelectric problems as well as the comparative effectivity analysis of the analogous mechanical, electric and electromechanical problems are unique—no other examples of such results can be found in the related literature.

Acknowledgments

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