## Abstract

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

### Keywords

- adaptivity
- modeling
- simulation
- finite elements
- hierarchical models
- hierarchical approximations
- error estimation
- adaptivity control
- algorithms
- effectivity
- elasticity
- dielectricity
- piezoelectricity

## 1. Introduction

This chapter constitutes a continuation and extension of the previous work [1] on theoretical and implementation difficulties in application of the adaptive hierarchical modeling and

In this chapter, attention is paid to the effectivity of the algorithms for adaptive modeling and simulation of three considered classes of physical phenomena, that is, elasticity, dielectricity and piezoelectricity. Effectivity of hierarchical approximations within elastic, dielectric and piezoelectric media is compared. For this purpose, convergence curves for the analogous model problems within three mentioned classes of problems are generated and assessed. Also, the exemplary comparative results of the approximations are presented for these three classes of model problems. In the case of the error estimation, the global and local (element) effectivity indices for the total, approximation and modeling errors, where the latter is the difference of the former two, in the exemplary model problems of elasticity, dielectricity and piezoelectricity are calculated and compared. The exemplary distribution of the element error indicators and the global values of the error estimators for the model problems of three classes are presented and compared. In the case of the adaptive procedures, the model- and

### 1.1. Research objectives

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

### 1.2. State-of-the-art issues

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

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

The hierarchical and constrained approximations necessary for

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

Finally, adaptivity control by means of the three-step strategy for simple structures was presented in [21]. Within this strategy, three subsequent meshes are generated—initial, intermediate (or

### 1.3. Novelty of the research

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

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

## 2. Model problems

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

### 2.1. Elastostatics

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

where

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

where

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

where

The element nodal forces vector due to the surface traction

corresponding to the bases and sides of the prismatic element. Above the element, bases and sides are defined with the normalized longitudinal coordinates

### 2.2. Electrostatics

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

In Eq. (5),

with

The nodal electric charges vector of the element

and

where

### 2.3. Stationary piezoelectricity

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

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

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

with

The element contributions to the other terms of (9) are defined as before, that is, in accordance with (2)–(4) and (6)–(7). Note that the different shape functions matrices,

## 3. The applied numerical techniques

### 3.1. Hierarchies of models

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

The mechanical hierarchy

with

where

Subsequently, the hierarchy

where

with

The 3D-based hierarchy

The hierarchy is composed of all combinations

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

The displacement field of the elastic and piezoelectric elements is defined through the interpolation function

The first component function

where the mentioned vector is:

The second component interpolant

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

The subsequent interpolant

where the vectors of nodal dofs are equal to:

The last component interpolant

where the dof vector is:

With

The function

The linear interpolant

where

The higher order interpolant

with

The next higher order interpolation function

while

The last component function

where

Note that when

Here, we discuss on the constrained approximation. Such an approximation is necessary for

where the parts

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

### 3.3. Error estimation

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

with

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:

In the above-mentioned equation, the term

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:

The coupled local solutions

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

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:

and

where *Se* and *S* are the element and body surfaces, while

### 3.4. Adaptive strategy

The adaptive strategy applied in this chapter takes advantage of the Texas three-step strategy [21]. The original strategy is assigned for structures of simple geometry (a single geometrical part of one type) and simple physical description (one model). Such a strategy is designed for

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

Description of the error-controlled adaptivity for elastic, dielectric and piezoelectric structures is started with the

where

In the case of the mechanical field, local (element)

with

In the case of

Above,

The

Above, the quantity

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

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

The value of the longitudinal (or overall) approximation order

where

Finally, in the case of

where

## 4. Numerical examples

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

### 4.1. Model structures

Here, the same domain geometry is considered in the problems of elasticity, dielectricity and piezoelectricity. Its square longitudinal dimensions are equal to

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

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 (

The

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

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

### 4.3. Effectivity of error estimation

In this section, results concerning application of the equilibrated residual method to the analogous model problems of elastostatics, electrostatics and stationary piezoelectricity are presented. The data of the problems are taken from Section 4.1. The thickness of the analyzed domains is now equal to *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

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 *avr*), for

The analogous estimated error distributions and the analogous plots of the effectivity indices versus

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

### 4.4. Convergence of hp -adaptivity

In Figures 19–24, the results illustrating both

In the first two figures, the purely elastic case is presented. Figure 19 presents the mesh obtained in the three-step adaptive process. Both the changes of the mesh density and the element changes of the longitudinal order of approximation can be seen in Figure 19. The initial mesh (not displayed) corresponds to the discretization parameters

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

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

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

In the case of the pure dielectricity, the relatively fine

In the case of the piezoelectricity, the changes in the common

## 5. Conclusions

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

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

The observations from the tests concerning hierarchical

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 (

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

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

## Acknowledgments

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