3. The applied numerical techniques

#### 3.1. Hierarchies of models

field-potential matrix). The specific form of the latter matrix in the case of the prismatic

The nodal electric charges vector of the element e has to be defined in a different way on the

nT e

<sup>e</sup> and c are the element shape functions vector and the scalar density of the surface

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

<sup>K</sup>Mqq, hp � <sup>K</sup>Cw<sup>r</sup>, <sup>h</sup><sup>π</sup> <sup>¼</sup> <sup>F</sup><sup>V</sup> <sup>þ</sup> <sup>F</sup>S,

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 pq and πr adaptive approximations of the vectorial displacement and scalar electric fields are proposed in (9), with the common

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

> BT e Cb e

The element contributions to the other terms of (9) are defined as before, that is, in accordance

displacements and potential fields are employed here due to the different orders of

ð�ξ2þ<sup>1</sup> 0

with (2)–(4) and (6)–(7). Note that the different shape functions matrices, N

<sup>C</sup>qq, hp <sup>þ</sup> <sup>K</sup>Ew<sup>r</sup>, <sup>h</sup><sup>π</sup> <sup>¼</sup> <sup>F</sup><sup>Q</sup>

of equations. The coupling is represented by the matrix K<sup>C</sup> in the following way

KT

kC e ¼ ð1 0 ð1 0

with C representing the piezoelectric (coupling) constants matrix.

c wspð ÞJ dξ2dξ<sup>1</sup> (7)

c wspð ÞJ dξ3dη<sup>i</sup> (8)

detð ÞJ dξ1dξ2dξ<sup>3</sup> (10)

e

and n e

, for the

(9)

ð�ξ2þ<sup>1</sup> 0

element can be found in the work [26].

and

where n

electric charges.

h-approximation.

2.3. Stationary piezoelectricity

prismatic element bases and sides, that is,

f Q e ¼ ð1 0

162 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

f Q e ¼ ð1 0 ð1 0 nT e

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 M of the 3D or 3D-based elastic models M reads:

$$\mathbf{M} \in \mathbf{M}, \quad \mathbf{M} = \{\Im D, \mathbf{M}I, \mathbf{R}M, \Im D/\mathbf{R}M, \mathbf{M}l/\mathbf{R}M\} \tag{11}$$

with 3D denoting three-dimensional elasticity, MI representing hierarchical shell models of higher order, RM being the first-order shell model corresponding to Reissner theory of shells and 3D=RM and MI=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:

$$\text{MI} = \{\text{M2}, \text{M3}, \text{M4}, \dots\},$$

$$\text{MM/RM} = \{\text{M2}/\text{RM}, \text{M3}/\text{RM}, \text{M4}/\text{RM}, \dots\} \tag{12}$$

where I represents the order of the hierarchical model MI. This order is equivalent to the order of polynomial constraints defining the transverse displacement.

Subsequently, the hierarchy E of 3D-based dielectric models E includes:

$$E \in \mathcal{E}, \quad \mathcal{E} = \{\Im D, E\} \tag{13}$$

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

$$E \mathbf{J} = \{E1, E2, E3, \dots\} \tag{14}$$

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 P of piezoelectric models P consists of the following component models:

$$P \in \mathcal{P}, \quad \mathcal{P} = \{(M, E) : M \in M, E \in E\} \tag{15}$$

The second component interpolant <sup>u</sup><sup>2</sup>ð Þ <sup>ξ</sup> , corresponding to the element edges, is equal to the product of the higher-order shape function matrices, N<sup>h</sup> and Nu, of the element horizontal and

The mentioned vectors of degrees of freedom at the horizontal and vertical nodes are as

1, 2, …, Iu, where Ih and Iu are the numbers of the horizontal and upright mid-edge nodes,

The subsequent interpolant <sup>u</sup><sup>3</sup>ð Þ <sup>ξ</sup> corresponds the higher order mid-base and mid-side nodes of the element. The function is obtained through the multiplication of the shape function

, q<sup>u</sup> ¼ …; q1,j,l

ð Þ¼ ξ Nhð Þ ξ q<sup>h</sup> þ Nuð Þ ξ q<sup>u</sup> (19)

; q3,j,l ;…

ð Þ¼ ξ Nbð Þ ξ q<sup>b</sup> þ Nsð Þ ξ q<sup>s</sup> (20)

ð Þ¼ ξ Nmð Þ ξ q<sup>m</sup> (21)

, with k standing for a dof number at this

¼ qv; qh; qu; qb; qs; q<sup>m</sup>

<sup>e</sup> (22)

� �<sup>T</sup> Eq. (17)

ð Þ ξ (23)

h i<sup>T</sup>

http://dx.doi.org/10.5772/intechopen.72265

, i ¼ 1, 2, …, Ih, j ¼

, q<sup>s</sup> ¼

165

; q2,j,l

with i ¼ 1, 2, …, Ib and j ¼ 1, 2, …, Is. Here, Ib and Is denote the num-

e

ð Þþ <sup>ξ</sup> <sup>f</sup><sup>4</sup>

ð Þ ξ of the element vertices is equal to the product of the vector n<sup>v</sup> of

ð Þ¼ ξ nvð Þ ξ w<sup>v</sup> (24)

h i<sup>T</sup>

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

vertical (upright) mid-edge nodes and the corresponding dof vectors:

u2

matrices, N<sup>b</sup> and Νs, by the corresponding dofs vectors in accordance with

where the vectors of nodal dofs are equal to: q<sup>b</sup> ¼ …; q1,i, <sup>k</sup>; q2,i, <sup>k</sup>; q3,i, <sup>k</sup>;…

bers of the mid-base and mid-side nodes, while k and l are dof numbers at these nodes.

u4

¼ ½ � Nvð Þ ξ ; Nhð Þ ξ ; Nuð Þ ξ ; Nbð Þ ξ ; Nvð Þ ξ ; Nmð Þ ξ and q

tric or piezoelectric element is also defined as a sum of four components

<sup>f</sup>ð Þ¼ <sup>ξ</sup> <sup>f</sup><sup>1</sup>

h i<sup>T</sup>

uð Þ¼ ξ N e ð Þ ξ q

ð Þþ <sup>ξ</sup> <sup>f</sup><sup>2</sup>

f1

shape functions for the vertices and the corresponding nodal vector of scalar dofs, that is,

The function f ¼ fð Þ ξ interpolating electric potential at any point ξ of the normalized dielec-

ð Þþ <sup>ξ</sup> <sup>f</sup><sup>3</sup>

The last component interpolant <sup>u</sup><sup>4</sup>ð Þ <sup>ξ</sup> , assigned to the element higher order middle node, is defined as a product of the shape function matrix N<sup>m</sup> and the corresponding dof vector:

u3

h i<sup>T</sup>

while k and l represent numbers of dofs at these nodes.

where the dof vector is: q<sup>m</sup> ¼ …; q1, <sup>k</sup>; q2, <sup>k</sup>; q3, <sup>k</sup>;…

can be written in the alternative standard form

follows: q<sup>h</sup> ¼ …; q1,i, <sup>k</sup>; q2,i, <sup>k</sup>; q3,i, <sup>k</sup>; …

…; q1,j,l

node.

With N e

; q2,j,l

h i<sup>T</sup>

The linear interpolant f<sup>1</sup>

; q3,i,l ;…

The hierarchy is composed of all combinations ð Þ M; E of the elastic models M from (11) and (12) and dielectric models E defined in (13) and (14), that is,

$$\begin{aligned} \text{P} &= \{ (\text{3D}, \text{3D}), (\text{MI}, \text{3D}), (\text{RM}, \text{3D}), (\text{3D}/\text{RM}, \text{3D}), (\text{MI}/\text{RM}, \text{3D}) \\ &\quad (\text{3D}, \text{E}), (\text{MI}, \text{E})), (\text{RM}, \text{E}), (\text{3D}/\text{RM}, \text{E}), (\text{MI}/\text{RM}, \text{E})) \} \end{aligned} \tag{16}$$

#### 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 threedimensional 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 u ¼ uð Þ ξ describing displacements u ¼ ð Þ u1; u2; u<sup>3</sup> <sup>T</sup> of any point ξ of the normalized geometry of the element. This interpolant is a sum of four component functions:

$$
\mathfrak{u}(\mathfrak{E}) = \mathfrak{u}^1(\mathfrak{E}) + \mathfrak{u}^2(\mathfrak{E}) + \mathfrak{u}^3(\mathfrak{E}) + \mathfrak{u}^4(\mathfrak{E}) \tag{17}
$$

The first component function <sup>u</sup><sup>1</sup>ð Þ <sup>ξ</sup> of the element vertices is defined as a product of the linear vertex node shape function matrix N<sup>v</sup> and the corresponding vector of nodal threesomes of directional dofs, that is,

$$\mathfrak{u}^1(\xi) = \mathcal{N}\_v(\xi)\,\mathfrak{q}\_v \tag{18}$$

where the mentioned vector is: q<sup>v</sup> ¼ …; q1,i ; q2,i ; q3,i ;… h i<sup>T</sup> , and where i ¼ 1, 2, …, Iv with Iv being the number of vertex nodes within the element.

The second component interpolant <sup>u</sup><sup>2</sup>ð Þ <sup>ξ</sup> , corresponding to the element edges, is equal to the product of the higher-order shape function matrices, N<sup>h</sup> and Nu, of the element horizontal and vertical (upright) mid-edge nodes and the corresponding dof vectors:

with J standing for the order of the hierarchical dielectric theory (the polynomial order of the

The 3D-based hierarchy P of piezoelectric models P consists of the following component models:

The hierarchy is composed of all combinations ð Þ M; E of the elastic models M from (11) and

P ¼ fð Þ 3D; 3D , MI ð Þ ; 3D , RMð Þ ; 3D ,ð Þ 3D=RM; 3D , MI ð Þ =RM; 3D ð Þ 3D; EJ , MI ð Þ ; EJ , RMð Þ ; EJ ,ð Þ 3D=RM; EJ , MI ð Þg =RM; EJ

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 threedimensional 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.

lation function u ¼ uð Þ ξ describing displacements u ¼ ð Þ u1; u2; u<sup>3</sup>

<sup>u</sup>ð Þ¼ <sup>ξ</sup> <sup>u</sup><sup>1</sup>

directional dofs, that is,

where the mentioned vector is: q<sup>v</sup> ¼ …; q1,i

being the number of vertex nodes within the element.

The displacement field of the elastic and piezoelectric elements is defined through the interpo-

The first component function <sup>u</sup><sup>1</sup>ð Þ <sup>ξ</sup> of the element vertices is defined as a product of the linear vertex node shape function matrix N<sup>v</sup> and the corresponding vector of nodal threesomes of

> ; q2,i ; q3,i ;…

h i<sup>T</sup>

ð Þþ <sup>ξ</sup> <sup>u</sup><sup>3</sup>

ð Þþ <sup>ξ</sup> <sup>u</sup><sup>4</sup>

ð Þ¼ ξ Nvð Þ ξ q<sup>v</sup> (18)

malized geometry of the element. This interpolant is a sum of four component functions:

ð Þþ <sup>ξ</sup> <sup>u</sup><sup>2</sup>

u1

P∈ P, P ¼ f g ð Þ M; E : M ∈ M; E∈E (15)

(16)

<sup>T</sup> of any point ξ of the nor-

ð Þ ξ (17)

, and where i ¼ 1, 2, …, Iv with Iv

through-thickness constraints of electric potential).

164 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

3.2. Hierarchical and constrained approximations

(12) and dielectric models E defined in (13) and (14), that is,

$$\mathfrak{u}^{2}(\xi) = \mathbf{N}\_{h}(\xi)\boldsymbol{q}\_{h} + \mathbf{N}\_{u}(\xi)\boldsymbol{q}\_{u} \tag{19}$$

The mentioned vectors of degrees of freedom at the horizontal and vertical nodes are as follows: q<sup>h</sup> ¼ …; q1,i, <sup>k</sup>; q2,i, <sup>k</sup>; q3,i, <sup>k</sup>; … h i<sup>T</sup> , q<sup>u</sup> ¼ …; q1,j,l ; q2,j,l ; q3,j,l ;… h i<sup>T</sup> , i ¼ 1, 2, …, Ih, j ¼ 1, 2, …, Iu, where Ih and Iu 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 <sup>u</sup><sup>3</sup>ð Þ <sup>ξ</sup> 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, N<sup>b</sup> and Νs, by the corresponding dofs vectors in accordance with

$$\mu^3(\xi) = \mathbf{N}\_b(\xi)\,\boldsymbol{q}\_b + \mathbf{N}\_s(\xi)\,\boldsymbol{q}\_s \tag{20}$$

where the vectors of nodal dofs are equal to: q<sup>b</sup> ¼ …; q1,i, <sup>k</sup>; q2,i, <sup>k</sup>; q3,i, <sup>k</sup>;… h i<sup>T</sup> , q<sup>s</sup> ¼ …; q1,j,l ; q2,j,l ; q3,i,l ;… h i<sup>T</sup> with i ¼ 1, 2, …, Ib and j ¼ 1, 2, …, Is. Here, Ib and Is 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 <sup>u</sup><sup>4</sup>ð Þ <sup>ξ</sup> , assigned to the element higher order middle node, is defined as a product of the shape function matrix N<sup>m</sup> and the corresponding dof vector:

$$\mathfrak{u}^{4}(\mathfrak{E}) = \mathbf{N}\_{m}(\mathfrak{E})\,\mathfrak{q}\_{m} \tag{21}$$

where the dof vector is: q<sup>m</sup> ¼ …; q1, <sup>k</sup>; q2, <sup>k</sup>; q3, <sup>k</sup>;… h i<sup>T</sup> , with k standing for a dof number at this node.

With N e ¼ ½ � Nvð Þ ξ ; Nhð Þ ξ ; Nuð Þ ξ ; Nbð Þ ξ ; Nvð Þ ξ ; Nmð Þ ξ and q e ¼ qv; qh; qu; qb; qs; q<sup>m</sup> � �<sup>T</sup> Eq. (17) can be written in the alternative standard form

$$\mathfrak{u}(\xi) = \overset{\varepsilon}{\mathbf{N}}(\xi)\overset{\varepsilon}{\mathfrak{q}} \tag{22}$$

The function f ¼ fð Þ ξ interpolating electric potential at any point ξ of the normalized dielectric or piezoelectric element is also defined as a sum of four components

$$
\phi(\xi) = \phi^1(\xi) + \phi^2(\xi) + \phi^3(\xi) + \phi^4(\xi) \tag{23}
$$

The linear interpolant f<sup>1</sup> ð Þ ξ of the element vertices is equal to the product of the vector n<sup>v</sup> of shape functions for the vertices and the corresponding nodal vector of scalar dofs, that is,

$$\phi^1(\xi) = \mathfrak{n}\_v(\xi)\,\varphi\_v \tag{24}$$

where <sup>w</sup><sup>v</sup> <sup>¼</sup> …; <sup>w</sup><sup>i</sup> ½ � ;… <sup>T</sup> and <sup>i</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, …, Jv, with Jv denoting the number of vertices within the element.

The higher order interpolant f<sup>2</sup> ð Þ ξ , corresponding to horizontal and vertical (upright) midedge nodes, is equal to a product of the shape function vectors n<sup>h</sup> and n<sup>u</sup> and the corresponding vectors of nodal dofs:

$$
\phi^2(\xi) = \mathfrak{n}\_{\hbar}(\xi)\,\varphi\_{\hbar} + \mathfrak{n}\_{\hbar}(\xi)\,\varphi\_{\hbar} \tag{25}
$$

q f <sup>¼</sup> <sup>q</sup><sup>s</sup> f

w f

two-dimensional and three-dimensional cases are described.

<sup>M</sup> standing for the element stiffness matrix, f

e

calculated as the difference of the previous two errors.

<sup>¼</sup> <sup>w</sup><sup>s</sup> f

> k e <sup>M</sup>q

2 6 4

wu f

where the parts q<sup>s</sup>

3.3. Error estimation

surface forces vectors, and f

interelement stress loadings. In (31), q

with k e

element f , q<sup>w</sup>

Cq fe

f

and q<sup>u</sup> f of q f

electric potential, the analogous relation reads:

qu f

I 0

3 5

<sup>e</sup> contains displacements of the constraining nodes of the bigger neighbor <sup>e</sup>, while

3 5 2 4 ws f

3

ww e

2 6 4

qs f

3 7

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

include the unconstrained and constrained dofs of the smaller

<sup>5</sup> (29)

http://dx.doi.org/10.5772/intechopen.72265

167

5 (30)

<sup>R</sup> (31)

<sup>S</sup> denoting the nodal mass and

qw e

2 4

0 C<sup>q</sup> fe

and I represent the constraint coefficient matrix and the unity matrix. In the case of the

I 0

2 4

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

The equilibrated residual method of error estimation [14–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:

imated 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

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:

e

<sup>V</sup> and f e

<sup>R</sup> representing the nodal forces vector due to the equilibrated

<sup>e</sup> Q,HP is the solution displacements vector in the approx-

<sup>e</sup> Q,HP <sup>¼</sup> <sup>f</sup> e <sup>V</sup> þ f e <sup>S</sup> þ f e

0 C<sup>w</sup> fe

2 6 4

with <sup>w</sup><sup>h</sup> <sup>¼</sup> …; <sup>w</sup>i, <sup>k</sup>;… <sup>T</sup> , <sup>w</sup><sup>u</sup> <sup>¼</sup> …; <sup>w</sup>i,l;… <sup>T</sup> , i ¼ 1, 2, , …, Jh, j ¼ 1, 2, …, Ju. Above, the numbers of the horizontal and upright mid-edge nodes are equal to Jh and Ju, respectively, while k and l are numbers of the consecutive dofs at these nodes.

The next higher order interpolation function f<sup>3</sup> ð Þ ξ , dealing with the mid-base and mid-side nodes, can be calculated with the multiplication of the shape function vectors n<sup>b</sup> and n<sup>s</sup> of these nodes and the respective vectors of nodal dofs:

$$\left(\boldsymbol{\phi}^{\mathsf{T}}(\boldsymbol{\xi}) = \mathfrak{n}\_{b}(\boldsymbol{\xi})\,\boldsymbol{\varphi}\_{b} + \mathfrak{n}\_{s}(\boldsymbol{\xi})\,\boldsymbol{\varphi}\_{s} \tag{26}$$

while <sup>w</sup><sup>b</sup> <sup>¼</sup> …; <sup>w</sup>i, <sup>k</sup>;… <sup>T</sup> and <sup>w</sup><sup>s</sup> <sup>¼</sup> …; <sup>w</sup>i,l; … <sup>T</sup> . Additionally, i ¼ 1, 2, …, Jb and j ¼ 1, 2, …, Js with Jb and Js 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 f<sup>4</sup> ð Þ ξ , assigned for the middle node, needs multiplication of the shape function vector n<sup>m</sup> of the node and the corresponding vector of the nodal dofs

$$\phi^4(\xi) = \mathfrak{n}\_m(\xi)\,\varphi\_m \tag{27}$$

where <sup>w</sup><sup>m</sup> <sup>¼</sup> …; <sup>w</sup><sup>k</sup> ½ � ; … <sup>T</sup>, and <sup>k</sup> represents the number of a hierarchical dof at the middle node.

Note that when n e ¼ ½ � nvð Þ ξ ; nhð Þ ξ ; nuð Þ ξ ; nbð Þ ξ ; nsð Þ ξ ; nmð Þ ξ and w <sup>e</sup> <sup>¼</sup> <sup>w</sup>v;wh;wu; <sup>w</sup>b; <sup>w</sup><sup>s</sup> ½ � ;w<sup>m</sup> T, Eq. (23) can be written in the well-known general form

$$
\phi(\xi) = \overset{\varepsilon}{\mathfrak{n}}(\xi)\overset{\varepsilon}{\varphi} \tag{28}
$$

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:

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems http://dx.doi.org/10.5772/intechopen.72265 167

$$\begin{aligned} \stackrel{f}{\dot{\boldsymbol{q}}} = \begin{bmatrix} \stackrel{f}{\dot{\boldsymbol{q}}}\_{s} \\ \stackrel{f}{\dot{\boldsymbol{q}}}\_{u} \end{bmatrix} = \begin{bmatrix} \boldsymbol{I} & \mathbf{0} \\ \mathbf{0} & \stackrel{\dot{\boldsymbol{\rho}}}{\dot{\mathbf{C}}}\_{q} \end{bmatrix} \begin{bmatrix} \stackrel{f}{\dot{\boldsymbol{q}}}\_{s} \\ \stackrel{\boldsymbol{\rho}}{\dot{\boldsymbol{q}}}\_{w} \end{bmatrix} \tag{29} \end{aligned} \tag{29}$$

where the parts q<sup>s</sup> f and q<sup>u</sup> f of q f include the unconstrained and constrained dofs of the smaller element f , q<sup>w</sup> <sup>e</sup> contains displacements of the constraining nodes of the bigger neighbor <sup>e</sup>, while Cq fe and I represent the constraint coefficient matrix and the unity matrix. In the case of the electric potential, the analogous relation reads:

$$
\stackrel{\circ}{\boldsymbol{\varphi}} = \begin{bmatrix} \stackrel{\circ}{\boldsymbol{\varphi}}\_{\boldsymbol{s}} \\ \stackrel{\circ}{\boldsymbol{\varphi}}\_{\boldsymbol{u}} \end{bmatrix} = \begin{bmatrix} \boldsymbol{I} & \mathbf{0} \\ \mathbf{0} & \stackrel{\circ}{\boldsymbol{C}}\_{\boldsymbol{\varphi}} \end{bmatrix} \begin{bmatrix} \stackrel{\circ}{\boldsymbol{\varphi}}\_{\boldsymbol{s}} \\ \stackrel{\circ}{\boldsymbol{\varphi}}\_{\boldsymbol{w}} \end{bmatrix} \tag{30}
$$

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

where <sup>w</sup><sup>v</sup> <sup>¼</sup> …; <sup>w</sup><sup>i</sup> ½ � ;… <sup>T</sup> and <sup>i</sup> <sup>¼</sup> <sup>1</sup>, <sup>2</sup>, …, Jv, with Jv denoting the number of vertices within the

edge nodes, is equal to a product of the shape function vectors n<sup>h</sup> and n<sup>u</sup> and the

of the horizontal and upright mid-edge nodes are equal to Jh and Ju, respectively, while k and l

nodes, can be calculated with the multiplication of the shape function vectors n<sup>b</sup> and n<sup>s</sup> of these

with Jb and Js standing for the numbers of the mid-base and mid-side nodes, respectively, and

where <sup>w</sup><sup>m</sup> <sup>¼</sup> …; <sup>w</sup><sup>k</sup> ½ � ; … <sup>T</sup>, and <sup>k</sup> represents the number of a hierarchical dof at the middle node.

¼ ½ � nvð Þ ξ ; nhð Þ ξ ; nuð Þ ξ ; nbð Þ ξ ; nsð Þ ξ ; nmð Þ ξ and w

fð Þ¼ ξ n e ð Þ ξ w

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

shape function vector n<sup>m</sup> of the node and the corresponding vector of the nodal dofs

f4

and vectors. The constraining relation for the case of displacements reads:

f2

166 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

, <sup>w</sup><sup>u</sup> <sup>¼</sup> …; <sup>w</sup>i,l;… <sup>T</sup>

f3

are numbers of the consecutive dofs at these nodes.

The next higher order interpolation function f<sup>3</sup>

nodes and the respective vectors of nodal dofs:

while <sup>w</sup><sup>b</sup> <sup>¼</sup> …; <sup>w</sup>i, <sup>k</sup>;… <sup>T</sup> and <sup>w</sup><sup>s</sup> <sup>¼</sup> …; <sup>w</sup>i,l; … <sup>T</sup>

k, l being dofs numbers at these nodes.

e

Eq. (23) can be written in the well-known general form

The last component function f<sup>4</sup>

Note that when n

ð Þ ξ , corresponding to horizontal and vertical (upright) mid-

ð Þ¼ ξ nhð Þ ξ w<sup>h</sup> þ nuð Þ ξ w<sup>u</sup> (25)

ð Þ¼ ξ nbð Þ ξ w<sup>b</sup> þ nsð Þ ξ w<sup>s</sup> (26)

ð Þ ξ , assigned for the middle node, needs multiplication of the

ð Þ¼ ξ nmð Þ ξ w<sup>m</sup> (27)

<sup>e</sup> <sup>¼</sup> <sup>w</sup>v;wh;wu; <sup>w</sup>b; <sup>w</sup><sup>s</sup> ½ � ;w<sup>m</sup>

<sup>e</sup> (28)

T,

, i ¼ 1, 2, , …, Jh, j ¼ 1, 2, …, Ju. Above, the numbers

ð Þ ξ , dealing with the mid-base and mid-side

. Additionally, i ¼ 1, 2, …, Jb and j ¼ 1, 2, …, Js

element.

The higher order interpolant f<sup>2</sup>

with <sup>w</sup><sup>h</sup> <sup>¼</sup> …; <sup>w</sup>i, <sup>k</sup>;… <sup>T</sup>

corresponding vectors of nodal dofs:

The equilibrated residual method of error estimation [14–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:

$$
\hat{\mathbf{k}}\_M \stackrel{\epsilon}{\mathfrak{q}}^{Q, HP} = \stackrel{\epsilon}{\tilde{\mathbf{f}}}\_V + \stackrel{\epsilon}{\tilde{\mathbf{f}}}\_S + \stackrel{\epsilon}{\tilde{\mathbf{f}}}\_R \tag{31}
$$

with k e <sup>M</sup> standing for the element stiffness matrix, f e <sup>V</sup> and f e <sup>S</sup> denoting the nodal mass and surface forces vectors, and f e <sup>R</sup> representing the nodal forces vector due to the equilibrated interelement stress loadings. In (31), q <sup>e</sup> Q,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:

$$
\stackrel{\varepsilon}{\mathbf{k}}\_E \stackrel{\varepsilon}{\mathfrak{o}}^{P,H\Pi} = \stackrel{\varepsilon}{\mathbf{f}}\_Q + \stackrel{\varepsilon}{\mathbf{f}}\_H \tag{32}
$$

where Se and S are the element and body surfaces, while r<sup>e</sup> uhpq and he f<sup>h</sup>π<sup>r</sup> 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 Sef between the elements. Some more details on how to calculate the equilibrated stresses and charge can be found in

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

169

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 hpadaptivity 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 enrich-

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 analo-

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 nIM in the intermediate mesh (denoted with index I),

<sup>0</sup><sup>M</sup> EI

is the strain energy norm of the solution from the initial mesh, while the

(37)

γ2 a,I u<sup>0</sup> 2 U 

where η<sup>0</sup><sup>M</sup> is the estimated value of the approximation error from the initial mesh, μ<sup>0</sup><sup>M</sup> represents the known or assumed h-convergence rate, d denotes dimensionality (equal to 2 or 3) of the adapted geometrical part. Additionally, EI is the total number of elements in the interme-

coefficient γa,I determines the expected relative value of the global approximation error within

ment). The problem solution on this mesh is followed by the error estimation again.

gous comments for the case of elasticity can be found in [29, 30].

n <sup>2</sup>μ0<sup>M</sup> <sup>=</sup>dþ<sup>1</sup> IM <sup>¼</sup> <sup>η</sup><sup>2</sup>

which replace an element of the initial mesh, is equal to:

[3, 19, 26].

3.4. Adaptive strategy

diate mesh, u<sup>0</sup>

the intermediate mesh.

2 U 

In the above-mentioned equation, the term k e <sup>E</sup> stands for the element characteristic matrix of dielectricity, f e <sup>Q</sup> is the element nodal vector due to surface charges, while f e <sup>H</sup> represents nodal forces due to the equilibrated interelement charges. The vector w <sup>e</sup> <sup>P</sup>,H<sup>Π</sup> is the solution vector of electric potential in the local problem. The solution is approximate and the corresponding discretization parameters H, Π, P (the element size and longitudinal and transverse approximation orders) can be different to their global counterparts h, π, r. In the cases of the approximation and total errors, the values of H ¼ h, Π ¼ π þ 1, P ¼ r and H ¼ h, Π ¼ π þ 1, P ¼ r þ 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{aligned} \stackrel{\varepsilon}{\mathbf{k}}\_{\mathcal{M}} \stackrel{\varepsilon}{\mathfrak{q}}^{\mathcal{Q},HP} - \stackrel{\varepsilon}{\mathbf{k}}\_{\mathcal{C}} \stackrel{\varepsilon}{\mathfrak{q}}^{\mathcal{P},H\Pi} &= \stackrel{\varepsilon}{\mathfrak{f}}\_{V} + \stackrel{\varepsilon}{\mathfrak{f}}\_{\mathcal{S}} + \stackrel{\varepsilon}{\mathfrak{f}}\_{R} \\ \stackrel{\varepsilon}{\mathbf{k}}\_{\mathcal{C}}^{T} \stackrel{\varepsilon}{\mathfrak{q}}^{\mathcal{Q},HP} + \stackrel{\varepsilon}{\mathbf{k}}\_{E} \stackrel{\varepsilon}{\mathfrak{p}}^{\mathcal{P},H\Pi} &= \stackrel{\varepsilon}{\mathfrak{f}}\_{Q} + \stackrel{\varepsilon}{\mathfrak{f}}\_{H} \end{aligned} \tag{33}$$

The coupled local solutions q <sup>e</sup> Q,HP and w <sup>e</sup> <sup>P</sup>,H<sup>Π</sup> 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{aligned} \stackrel{\epsilon}{\mathbf{k}}\_{M} \stackrel{\epsilon}{\mathfrak{q}}^{Q,HP} &= \stackrel{\epsilon}{\mathbf{k}}\_{C} \stackrel{\epsilon}{\mathfrak{q}}^{p,ln} + \stackrel{\epsilon}{\mathfrak{f}}\_{V} + \stackrel{\epsilon}{\mathfrak{f}}\_{S} + \stackrel{\epsilon}{\mathfrak{f}}\_{R} \\ \stackrel{\epsilon}{\mathbf{k}}\_{E} \stackrel{\epsilon}{\mathfrak{q}}^{P,H\Pi} &= -\stackrel{\epsilon}{\mathbf{k}}\_{C} \stackrel{\epsilon}{\mathfrak{q}}^{q,hp} + \stackrel{\epsilon}{\mathfrak{f}}\_{Q} + \stackrel{\epsilon}{\mathfrak{f}}\_{H} \end{aligned} \tag{34}$$

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:

$$\boldsymbol{\hat{f}}\_{R}^{\varepsilon} = \int\_{S\_{\varepsilon} \backslash \mathcal{S}} \boldsymbol{\mathsf{N}}^{\varepsilon} \left< \boldsymbol{r}\_{\varepsilon} \left( \boldsymbol{\mathsf{u}}^{\mathrm{lup}} \right) \right> d\mathcal{S}\_{\varepsilon} = \sum\_{f} \int\_{S\_{\mathcal{f}}} \boldsymbol{\mathsf{N}}^{\varepsilon} \left< \boldsymbol{r}\_{\varepsilon} \left( \boldsymbol{\mathsf{u}}^{\mathrm{lup}} \right) \right> d\mathcal{S}\_{\varepsilon'} \tag{35}$$

and

$$\oint \stackrel{\epsilon}{H} = \int\_{S\_{\epsilon} \cup \mathcal{S}} \mathfrak{n}^{T} \left< h\_{\epsilon} \left( \phi^{\hbar \pi \rho} \right) \right> d\mathcal{S}\_{\epsilon} = \sum\_{f} \int\_{S \! f} \mathfrak{n}^{T} \left< h\_{\epsilon} \left( \phi^{\hbar \pi \rho} \right) \right> d\mathcal{S}\_{\epsilon f} \tag{36}$$

where Se and S are the element and body surfaces, while r<sup>e</sup> uhpq and he f<sup>h</sup>π<sup>r</sup> 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 Sef 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

k e <sup>E</sup>w

forces due to the equilibrated interelement charges. The vector w

168 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

k e <sup>M</sup> q

> k e T Cq

k e <sup>M</sup> q

nodal forces and charge are defined in accordance with:

f R e ¼ ð Se\S N<sup>T</sup> e

f H e ¼ ð Se\S nT e

k e <sup>E</sup>w

<sup>e</sup> Q,HP � <sup>k</sup> e <sup>C</sup> w

<sup>e</sup> Q,HP <sup>þ</sup> <sup>k</sup> e <sup>E</sup>w

residual-based global estimators obtained from the local solutions of (33).

<sup>e</sup> Q,HP <sup>¼</sup> <sup>k</sup> e <sup>C</sup> w <sup>e</sup> <sup>r</sup>,h<sup>π</sup> <sup>þ</sup> <sup>f</sup> e <sup>V</sup> þ f e <sup>S</sup> þ f e R

<sup>e</sup> <sup>P</sup>,H<sup>Π</sup> ¼ �<sup>k</sup>

<sup>e</sup> Q,HP and w

In the above-mentioned equation, the term k

dielectricity, f

e

The coupled local solutions q

alternative:

and

<sup>e</sup> <sup>P</sup>,H<sup>Π</sup> <sup>¼</sup> <sup>f</sup>

e

<sup>Q</sup> is the element nodal vector due to surface charges, while f

imation errors, and the modeling error is defined as the difference of these two errors.

electric potential in the local problem. The solution is approximate and the corresponding discretization parameters H, Π, P (the element size and longitudinal and transverse approximation orders) can be different to their global counterparts h, π, r. In the cases of the approximation and total errors, the values of H ¼ h, Π ¼ π þ 1, P ¼ r and H ¼ h, Π ¼ π þ 1, P ¼ r þ 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 approx-

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:

<sup>e</sup> <sup>P</sup>,H<sup>Π</sup> <sup>¼</sup> <sup>f</sup>

<sup>e</sup> <sup>P</sup>,H<sup>Π</sup> <sup>¼</sup> <sup>f</sup>

the lack of the upper-bound property of the total, approximation and modeling errors by the

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

e T Cq <sup>e</sup> q,hp <sup>þ</sup> <sup>f</sup> e <sup>Q</sup> þ f e H

In Eqs. (31) and (32) and in the above two sets of equations, the vectors of the equilibrated

f

f

ð Sef N<sup>T</sup> e

ð Sef nT e

<sup>r</sup><sup>e</sup> <sup>u</sup>hpq � � � � dSe <sup>¼</sup> <sup>X</sup>

he <sup>f</sup><sup>h</sup>π<sup>r</sup> � � � � dSe <sup>¼</sup> <sup>X</sup>

e <sup>V</sup> þ f e <sup>S</sup> þ f e R

e <sup>Q</sup> þ f e H

<sup>e</sup> <sup>P</sup>,H<sup>Π</sup> of the above set have one disadvantage. It lies in

r<sup>e</sup> uhpq � � � � dSef (35)

he f<sup>h</sup>π<sup>r</sup> � � � � dSef (36)

e <sup>Q</sup> þ f e

<sup>H</sup> (32)

e

<sup>e</sup> <sup>P</sup>,H<sup>Π</sup> is the solution vector of

<sup>H</sup> represents nodal

(33)

(34)

<sup>E</sup> stands for the element characteristic matrix of

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 hpadaptivity 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 nIM in the intermediate mesh (denoted with index I), which replace an element of the initial mesh, is equal to:

$$m\_{I\_M}^{2\mu\_{0\_M}/d+1} = \frac{\eta\_{0\_M}^2 E\_I}{\mathcal{V}\_{a,I}^2 ||\mathfrak{u}\_0||\_{II}^2} \tag{37}$$

where η<sup>0</sup><sup>M</sup> is the estimated value of the approximation error from the initial mesh, μ<sup>0</sup><sup>M</sup> represents the known or assumed h-convergence rate, d denotes dimensionality (equal to 2 or 3) of the adapted geometrical part. Additionally, EI is the total number of elements in the intermediate mesh, u<sup>0</sup> 2 U is the strain energy norm of the solution from the initial mesh, while the coefficient γ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 pT 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:

$$p\_T^{2v\_{0M}} = \frac{p\_0^{2v\_{0M}} \eta\_{I\_M}^2}{\gamma\_{a,T}^2 ||\mu\_0||\_{II}^2} \tag{38}$$

π 2ν0<sup>M</sup> <sup>T</sup> <sup>¼</sup> <sup>π</sup>

mined. The following formula can be proposed for this purpose [26]:

<sup>r</sup><sup>T</sup> <sup>¼</sup> <sup>r</sup><sup>0</sup> � <sup>1</sup>

modeling error estimated value of a finite element of the intermediate mesh.

order within the initial mesh.

4. Numerical examples

4.1. Model structures

thin piezoelectric structure in the electromechanical case.

2ν0<sup>E</sup> <sup>0</sup><sup>E</sup> <sup>η</sup><sup>2</sup> IE EI

2 W 

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

θ2 IE EI

γ2 m,t f<sup>0</sup> 2 W 

(42)

171

(43)

γ2 a,T f<sup>0</sup>

where ηIE is the estimated approximation error value from the intermediate mesh, ν<sup>0</sup><sup>E</sup> is the assumed problem-dependant p-convergence rate and π<sup>0</sup> is the longitudinal approximation

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 deter-

<sup>2</sup> log <sup>t</sup>=2<sup>l</sup>

where r<sup>0</sup> 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 θIE is the

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.

Here, the same domain geometry is considered in the problems of elasticity, dielectricity and piezoelectricity. Its square longitudinal dimensions are equal to 2<sup>l</sup> <sup>¼</sup> <sup>3</sup>:<sup>1415</sup> � <sup>10</sup>�<sup>2</sup> m, while its thickness equals <sup>t</sup> <sup>¼</sup> <sup>0</sup>:<sup>01</sup> � <sup>10</sup>�<sup>2</sup> 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

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].

with ηIM representing the estimated value of the approximation error from the intermediate mesh, ν<sup>0</sup><sup>M</sup> being the given p-convergence rate, p<sup>0</sup> denoting the longitudinal approximation order from the initial mesh and γ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 qT 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 ||\mathfrak{u}\_0||\_{\mathcal{U}}^2} \tag{39}$$

Above, q<sup>0</sup> 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, θIM is the estimated value of the modeling error from the intermediate mesh and γ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, nIE , 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:

$$m\_{l\_{\mathbb{E}}}^{2\mu\_{0\_{\mathbb{E}}}/d+1} = \frac{\eta\_{0\_{\mathbb{E}}}^2 E\_I}{\gamma\_{a,I}^2 \|\big| \phi\_0\big|\_{\mathbb{W}}^2} \tag{40}$$

Above, the quantity η<sup>0</sup><sup>E</sup> stands for the estimated value of the approximation error in the initial mesh, the exponent μ<sup>0</sup><sup>E</sup> is the assumed h-convergence rate and the norm f<sup>0</sup> 2 W 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

$$m\_l = \max\{n\_{l\_M}, n\_{l\_\mathbb{E}}\} \tag{41}$$

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

The value of the longitudinal (or overall) approximation order π<sup>T</sup> within the electric field of an element of the target (final) mesh can be calculated from [26]:

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems http://dx.doi.org/10.5772/intechopen.72265 171

$$
\pi\_T^{2v\_{0\_M}} = \frac{\pi\_{0\_E}^{2v\_{0\_E}} \eta\_{I\_E}^2 E\_I}{\gamma\_{a,T}^2 ||\phi\_0||\_W^2} \tag{42}
$$

where ηIE is the estimated approximation error value from the intermediate mesh, ν<sup>0</sup><sup>E</sup> is the assumed problem-dependant p-convergence rate and π<sup>0</sup> 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\_{l\_E}^2 E\_I}{\mathcal{V}\_{m,t}^2 ||\phi\_0||\_W^2} \tag{43}$$

where r<sup>0</sup> 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 θIE is the modeling error estimated value of a finite element of the intermediate mesh.

## 4. Numerical examples

In the case of the mechanical field, local (element) p-enrichments are controlled with the longitudinal (or overall) element approximation order pT 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

> 2ν0<sup>M</sup> <sup>0</sup> η<sup>2</sup> IM EI

2 U 

> θ2 IM EI

2 U 

γ2 m,T u<sup>0</sup>

<sup>0</sup><sup>E</sup> EI

(38)

(39)

(40)

2 W 

(41)

represents the

γ2 a,T u<sup>0</sup>

with ηIM representing the estimated value of the approximation error from the intermediate mesh, ν<sup>0</sup><sup>M</sup> being the given p-convergence rate, p<sup>0</sup> denoting the longitudinal approximation order from the initial mesh and γa,T standing for the expected relative value of the global

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

<sup>2</sup> log <sup>t</sup>=2<sup>l</sup>

Above, q<sup>0</sup> 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, θIM is the estimated value of the modeling error from the intermediate mesh and γm,T denotes the expected relative value of the modeling

The h-adaptivity within the electric field needs determination of the new number (denoted with index E) of elements, nIE , 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:

> γ2 a,I f<sup>0</sup> 2 W

Above, the quantity η<sup>0</sup><sup>E</sup> stands for the estimated value of the approximation error in the initial

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

nI ¼ max nIM ; nIE

The value of the longitudinal (or overall) approximation order π<sup>T</sup> within the electric field of an

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

p<sup>2</sup>ν0<sup>M</sup> <sup>T</sup> <sup>¼</sup> <sup>p</sup>

170 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

qT <sup>¼</sup> <sup>q</sup><sup>0</sup> � <sup>1</sup>

n <sup>2</sup>μ0<sup>E</sup> <sup>=</sup>dþ<sup>1</sup> IE <sup>¼</sup> <sup>η</sup><sup>2</sup>

mesh, the exponent μ<sup>0</sup><sup>E</sup> is the assumed h-convergence rate and the norm f<sup>0</sup>

electrostatic energy corresponding to the initial mesh.

element of the target (final) mesh can be calculated from [26]:

this parameter:

approximation error in the target mesh.

error in the target mesh.

with

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 2<sup>l</sup> <sup>¼</sup> <sup>3</sup>:<sup>1415</sup> � <sup>10</sup>�<sup>2</sup> m, while its thickness equals <sup>t</sup> <sup>¼</sup> <sup>0</sup>:<sup>01</sup> � <sup>10</sup>�<sup>2</sup> 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 <sup>E</sup> <sup>¼</sup> <sup>0</sup>:<sup>5</sup> � 1011 <sup>N</sup>=m2, while the applied Poisson'<sup>s</sup> ratio is ν ¼ 0:294. The isotropic dielectric properties of the dielectric and piezoelectric are characterized with the permittivity equal to <sup>δ</sup> <sup>¼</sup> <sup>0</sup>:<sup>1593</sup> � <sup>10</sup>�<sup>7</sup> F/m. The nonzero anisotropic piezoelectric constants are <sup>d</sup><sup>13</sup> <sup>¼</sup> <sup>d</sup><sup>23</sup> ¼ �0:<sup>15</sup> � <sup>10</sup>�<sup>9</sup> C/N, <sup>d</sup><sup>33</sup> <sup>¼</sup> <sup>0</sup>:<sup>3</sup> � <sup>10</sup>�<sup>9</sup> C/N and <sup>d</sup><sup>52</sup> <sup>¼</sup> <sup>d</sup><sup>61</sup> <sup>¼</sup> <sup>0</sup>:<sup>5</sup> � <sup>10</sup>�<sup>9</sup> C/N (compare [27] again). The surface load of magnitude <sup>p</sup> <sup>¼</sup> <sup>0</sup>:<sup>4</sup> � 106 <sup>N</sup>=m2 is applied to the upper surface of the elastic and piezoelectric structures. Furthermore, the surface electric charge of density <sup>c</sup> <sup>¼</sup> <sup>0</sup>:<sup>2</sup> � <sup>10</sup>�<sup>1</sup> <sup>C</sup>=m2 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 (sef) of a stress tensor (effective stress) and the electric displacement vector magnitude (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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

173

dielectric example due to the influence of the coupled mechanical displacements field.

Figure 2. Magnitude of electric displacement in the purely electric case.

Figure 3. Effective stress in the coupled problem.

Figure 1. Effective stress in the purely mechanical case.

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems http://dx.doi.org/10.5772/intechopen.72265 173

Figure 2. Magnitude of electric displacement in the purely electric case.

Young's modulus is assumed to be equal to <sup>E</sup> <sup>¼</sup> <sup>0</sup>:<sup>5</sup> � 1011 <sup>N</sup>=m2, while the applied Poisson'<sup>s</sup> ratio is ν ¼ 0:294. The isotropic dielectric properties of the dielectric and piezoelectric are characterized with the permittivity equal to <sup>δ</sup> <sup>¼</sup> <sup>0</sup>:<sup>1593</sup> � <sup>10</sup>�<sup>7</sup> F/m. The nonzero anisotropic piezoelectric constants are <sup>d</sup><sup>13</sup> <sup>¼</sup> <sup>d</sup><sup>23</sup> ¼ �0:<sup>15</sup> � <sup>10</sup>�<sup>9</sup> C/N, <sup>d</sup><sup>33</sup> <sup>¼</sup> <sup>0</sup>:<sup>3</sup> � <sup>10</sup>�<sup>9</sup> C/N and <sup>d</sup><sup>52</sup> <sup>¼</sup> <sup>d</sup><sup>61</sup> <sup>¼</sup> <sup>0</sup>:<sup>5</sup> � <sup>10</sup>�<sup>9</sup> C/N (compare [27] again). The surface load of magnitude <sup>p</sup> <sup>¼</sup> <sup>0</sup>:<sup>4</sup> � 106 <sup>N</sup>=m2 is applied to the upper surface of the elastic and piezoelectric structures. Furthermore, the surface electric charge of density <sup>c</sup> <sup>¼</sup> <sup>0</sup>:<sup>2</sup> � <sup>10</sup>�<sup>1</sup> <sup>C</sup>=m2 is applied to the upper surface of the dielectric and

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

Figures 1 and 2 illustrate distributions of the effective value (sef) of a stress tensor (effective stress) and the electric displacement vector magnitude (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

of the applied geometry, load and charge distributions and boundary conditions.

4.2. Convergence of hierarchical approximations

172 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 1. Effective stress in the purely mechanical case.

piezoelectric domains.

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 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 π-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 hconvergence) or the longitudinal approximation order p or π (the case of p-convergence). The applied values of the discretization parameters of the uniform meshes are m ¼ 1, 2, …, 8, p ¼ π ¼ 1, 2, …, 6 while the transverse orders of approximation are kept constant and equal to q � I ¼ 2 or/and r � J ¼ 2, with π � pi and r � 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 ¼ π ¼ 9, h ¼ 9, q ¼ r ¼ 2.

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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

175

the former case.

Figure 6. p-Convergence in the purely elastic problem.

Figure 5. h-Convergence in the purely elastic problem.

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems http://dx.doi.org/10.5772/intechopen.72265 175

Figure 5. h-Convergence in the purely elastic problem.

The h- and p-convergence curves for the elastostatic case are presented in Figures 5 and 6. The electrostatic problem solution h- and π-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 hconvergence) or the longitudinal approximation order p or π (the case of p-convergence). The applied values of the discretization parameters of the uniform meshes are m ¼ 1, 2, …, 8, p ¼ π ¼ 1, 2, …, 6 while the transverse orders of approximation are kept constant and equal to q � I ¼ 2 or/and r � J ¼ 2, with π � pi and r � 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

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

ones obtained from the meshes of p ¼ π ¼ 9, h ¼ 9, q ¼ r ¼ 2.

Figure 4. Magnitude of electric displacement in the coupled problem.

174 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 6. p-Convergence in the purely elastic problem.

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.

Figure 7. h-Convergence in the purely dielectric problem.

Figure 8. p-Convergence in the purely dielectric problem.

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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

177

pure problems, with a tendency to be closer to the purely mechanical case.

Figure 9. h-Convergence in the coupled piezoelectric problem.

Figure 10. p-Convergence (π ¼ p) in the coupled problem.

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems http://dx.doi.org/10.5772/intechopen.72265 177

Figure 9. h-Convergence in the coupled piezoelectric problem.

Figure 10. p-Convergence (π ¼ p) in the coupled problem.

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

Figure 7. h-Convergence in the purely dielectric problem.

176 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 8. p-Convergence in the purely dielectric problem.

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 <sup>t</sup> <sup>¼</sup> <sup>0</sup>:<sup>15</sup> � <sup>10</sup>�<sup>2</sup> m. The surface traction equals now <sup>p</sup> <sup>¼</sup> <sup>4</sup>:<sup>0</sup> � 106 <sup>N</sup>=m<sup>2</sup> . 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 ð Þ M 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,

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

179

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 ð Þ E nt), and also the average error estimate (marked with avr), for m ¼ 3, π � pi ¼ 4, r � 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 π. In order to obtain the necessary exact values of the energies, the global problem (5) was applied again, with m ¼ 9, π ¼ 9 and r ¼ 2 or r ¼ 6, for the calculation cases of the approximation and

The analogous estimated error distributions and the analogous plots of the effectivity indices versus p ¼ π 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 ¼ π ¼ 4, q ¼ r ¼ 2 and (33), respectively. The exact energy values necessary for effectivity calculations are obtained through the global numerical approximation (9), with

m ¼ 9, p ¼ π ¼ 9 and q ¼ r ¼ 2 or q ¼ r ¼ 6 and π � pi, r � rho.

Figure 12. Effectivities of the estimators in the purely mechanical case.

respectively.

total errors, respectively.

Figure 11. Estimated total errors in the purely mechanical case.

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.

4.3. Effectivity of error estimation

178 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 11. Estimated total errors in the purely mechanical case.

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 <sup>t</sup> <sup>¼</sup> <sup>0</sup>:<sup>15</sup> � <sup>10</sup>�<sup>2</sup> m. The surface traction equals now <sup>p</sup> <sup>¼</sup> <sup>4</sup>:<sup>0</sup> � 106 <sup>N</sup>=m<sup>2</sup>

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 ð Þ M 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,

.

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 ð Þ E nt), and also the average error estimate (marked with avr), for m ¼ 3, π � pi ¼ 4, r � 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 π. In order to obtain the necessary exact values of the energies, the global problem (5) was applied again, with m ¼ 9, π ¼ 9 and r ¼ 2 or r ¼ 6, for the calculation cases of the approximation and total errors, respectively.

The analogous estimated error distributions and the analogous plots of the effectivity indices versus p ¼ π 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 ¼ π ¼ 4, q ¼ r ¼ 2 and (33), respectively. The exact energy values necessary for effectivity calculations are obtained through the global numerical approximation (9), with m ¼ 9, p ¼ π ¼ 9 and q ¼ r ¼ 2 or q ¼ r ¼ 6 and π � pi, r � rho.

Figure 12. Effectivities of the estimators in the purely mechanical case.

Figure 13. Estimated total errors in the purely electric case.

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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

181

Figure 15. Mechanical parts of the estimated total errors (piezoelectricity).

Figure 16. Effectivities of estimators' mechanical parts (piezoelectricity).

Figure 14. Effectivities of the estimators in the purely electric case.

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

Figure 13. Estimated total errors in the purely electric case.

180 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 14. Effectivities of the estimators in the purely electric case.

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

Figure 16. Effectivities of estimators' mechanical parts (piezoelectricity).

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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

183

In Figures 19–24, the results illustrating both 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 <sup>t</sup> <sup>¼</sup> <sup>0</sup>:<sup>15</sup> � <sup>10</sup>�<sup>2</sup> 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

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

problems.

4.4. Convergence of hp-adaptivity

fields by means of the formulas (37), (38) or/and (40), (42).

Figure 19. Final hp-adapted mesh in the purely mechanical case.

Figure 17. Electrical parts of the estimated total errors (piezoelectricity).

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

Figure 18. Effectivities of estimators' electrical parts (piezoelectricity).

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.

### 4.4. Convergence of hp-adaptivity

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

Figure 17. Electrical parts of the estimated total errors (piezoelectricity).

182 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 18. Effectivities of estimators' electrical parts (piezoelectricity).

In Figures 19–24, the results illustrating both 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 <sup>t</sup> <sup>¼</sup> <sup>0</sup>:<sup>15</sup> � <sup>10</sup>�<sup>2</sup> 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).

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

Figure 19. Final hp-adapted mesh in the purely mechanical case.

Figure 20. hp-Adaptive convergence in the purely mechanical case.

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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

185

Figure 22. hπ-Adaptive convergence in the purely electric case.

Figure 23. Final hp-adapted mesh in the coupled problem.

Figure 21. Final hπ-adapted mesh in the purely electric case.

Figure 22. hπ-Adaptive convergence in the purely electric case.

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

Figure 20. hp-Adaptive convergence in the purely mechanical case.

184 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

Figure 21. Final hπ-adapted mesh in the purely electric case.

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

Figure 23. Final hp-adapted mesh in the coupled problem.

Figure 24. hp-Adaptive convergence in the coupled problem.

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.

again. The numerical approximation of the exact solution entering the error calculations is

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

187

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

In the case of the pure dielectricity, the relatively fine h-mesh is produced because of the error overestimation for π ¼ 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 π ¼ 2 (Figure 18). The following π-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

obtained from the discretization based on m ¼ 9, p ¼ π ¼ 9, q ¼ r ¼ 2.

Figure 25. Final hπ-adapted mesh in the coupled problem.

p ¼ 4 does not occur (compare Figure 12). Then the mesh is well p-enriched.

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 π ¼ 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 π ¼ 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 hπ-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, π ¼ 9 and r ¼ 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 π ¼ 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

Figure 25. Final hπ-adapted mesh in the coupled problem.

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

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 π ¼ 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 π ¼ 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 hπ-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

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 π ¼ 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

compared with the horizontal dotted line corresponding to the admissible error level.

numerical solution corresponding to m ¼ 9, π ¼ 9 and r ¼ 2.

Figure 24. hp-Adaptive convergence in the coupled problem.

186 Finite Element Method - Simulation, Numerical Analysis and Solution Techniques

again. The numerical approximation of the exact solution entering the error calculations is obtained from the discretization based on m ¼ 9, p ¼ π ¼ 9, q ¼ r ¼ 2.

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 π ¼ 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 π ¼ 2 (Figure 18). The following π-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

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 ≤ 3 or/and π ≤ 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

Adaptive Modeling and Simulation of Elastic, Dielectric and Piezoelectric Problems

http://dx.doi.org/10.5772/intechopen.72265

189

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 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

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 exam-

Partial financial support of the Polish Scientific Research Committee (now the National Science

[1] Zboiński G. Problems of hierarchical modelling and hp-adaptive finite element analysis in elasticity, dielectricity and piezoelectricity. In: Petrova R, editor. Perusal of the Finite Element

Centre) under the research grant no. N N504 5153 040 is thankfully acknowledged.

1 Institute of Fluid Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland 2 Faculty of Technical Sciences, University of Warmia and Mazury, Olsztyn, Poland

case of piezoelectricity, so as to enforce values of the effectivities closer to 1.0.

estimation procedure, occurs.

Acknowledgements

Author details

References

Grzegorz Zboiński1,2\*

ples of such results can be found in the related literature.

\*Address all correspondence to: zboi@imp.gda.pl

Method. Rijeka: InTech; 2016. Ch. 1. pp. 1-29

Figure 26. hπ-Adaptive convergence in the coupled problem.

corresponding admissible error value is far exceeded. The obtained solution for displacements is better than expected.
